/*
    * Licensed to the Apache Software Foundation (ASF) under one or more
    * contributor license agreements.  See the NOTICE file distributed with
    * this work for additional information regarding copyright ownership.
    * The ASF licenses this file to You under the Apache License, Version 2.0
    * (the "License"); you may not use this file except in compliance with
    * the License.  You may obtain a copy of the License at
    *
    *      http://www.apache.org/licenses/LICENSE-2.0
   *
   * Unless required by applicable law or agreed to in writing, software
   * distributed under the License is distributed on an "AS IS" BASIS,
   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
   * See the License for the specific language governing permissions and
   * limitations under the License.
   */
  
  package org.apache.commons.math4.legacy.linear;
  
  import org.apache.commons.math4.legacy.exception.DimensionMismatchException;
  import org.apache.commons.math4.core.jdkmath.JdkMath;
  
  /**
   * Calculates the LUP-decomposition of a square matrix.
   * <p>The LUP-decomposition of a matrix A consists of three matrices L, U and
   * P that satisfy: P&times;A = L&times;U. L is lower triangular (with unit
   * diagonal terms), U is upper triangular and P is a permutation matrix. All
   * matrices are m&times;m.</p>
   * <p>As shown by the presence of the P matrix, this decomposition is
   * implemented using partial pivoting.</p>
   * <p>This class is based on the class with similar name from the
   * <a href="http://math.nist.gov/javanumerics/jama/">JAMA</a> library.</p>
   * <ul>
   *   <li>a {@link #getP() getP} method has been added,</li>
   *   <li>the {@code det} method has been renamed as {@link #getDeterminant()
   *   getDeterminant},</li>
   *   <li>the {@code getDoublePivot} method has been removed (but the int based
   *   {@link #getPivot() getPivot} method has been kept),</li>
   *   <li>the {@code solve} and {@code isNonSingular} methods have been replaced
   *   by a {@link #getSolver() getSolver} method and the equivalent methods
   *   provided by the returned {@link DecompositionSolver}.</li>
   * </ul>
   *
   * @see <a href="http://mathworld.wolfram.com/LUDecomposition.html">MathWorld</a>
   * @see <a href="http://en.wikipedia.org/wiki/LU_decomposition">Wikipedia</a>
   * @since 2.0 (changed to concrete class in 3.0)
   */
  public class LUDecomposition {
      /** Default bound to determine effective singularity in LU decomposition. */
      private static final double DEFAULT_TOO_SMALL = 1e-11;
      /** Entries of LU decomposition. */
      private final double[][] lu;
      /** Pivot permutation associated with LU decomposition. */
      private final int[] pivot;
      /** Parity of the permutation associated with the LU decomposition. */
      private boolean even;
      /** Singularity indicator. */
      private boolean singular;
      /** Cached value of L. */
      private RealMatrix cachedL;
      /** Cached value of U. */
      private RealMatrix cachedU;
      /** Cached value of P. */
      private RealMatrix cachedP;
  
      /**
       * Calculates the LU-decomposition of the given matrix.
       * This constructor uses 1e-11 as default value for the singularity
       * threshold.
       *
       * @param matrix Matrix to decompose.
       * @throws NonSquareMatrixException if matrix is not square.
       */
      public LUDecomposition(RealMatrix matrix) {
          this(matrix, DEFAULT_TOO_SMALL);
      }
  
      /**
       * Calculates the LU-decomposition of the given matrix.
       * @param matrix The matrix to decompose.
       * @param singularityThreshold threshold (based on partial row norm)
       * under which a matrix is considered singular
       * @throws NonSquareMatrixException if matrix is not square
       */
      public LUDecomposition(RealMatrix matrix, double singularityThreshold) {
          if (!matrix.isSquare()) {
              throw new NonSquareMatrixException(matrix.getRowDimension(),
                                                 matrix.getColumnDimension());
          }
  
          final int m = matrix.getColumnDimension();
          lu = matrix.getData();
          pivot = new int[m];
          cachedL = null;
          cachedU = null;
          cachedP = null;
  
          // Initialize permutation array and parity
          for (int row = 0; row < m; row++) {
             pivot[row] = row;
         }
         even     = true;
         singular = false;
 
         // Loop over columns
         for (int col = 0; col < m; col++) {
 
             // upper
             for (int row = 0; row < col; row++) {
                 final double[] luRow = lu[row];
                 double sum = luRow[col];
                 for (int i = 0; i < row; i++) {
                     sum -= luRow[i] * lu[i][col];
                 }
                 luRow[col] = sum;
             }
 
             // lower
             int max = col; // permutation row
             double largest = Double.NEGATIVE_INFINITY;
             for (int row = col; row < m; row++) {
                 final double[] luRow = lu[row];
                 double sum = luRow[col];
                 for (int i = 0; i < col; i++) {
                     sum -= luRow[i] * lu[i][col];
                 }
                 luRow[col] = sum;
 
                 // maintain best permutation choice
                 if (JdkMath.abs(sum) > largest) {
                     largest = JdkMath.abs(sum);
                     max = row;
                 }
             }
 
             // Singularity check
             if (JdkMath.abs(lu[max][col]) < singularityThreshold) {
                 singular = true;
                 return;
             }
 
             // Pivot if necessary
             if (max != col) {
                 double tmp = 0;
                 final double[] luMax = lu[max];
                 final double[] luCol = lu[col];
                 for (int i = 0; i < m; i++) {
                     tmp = luMax[i];
                     luMax[i] = luCol[i];
                     luCol[i] = tmp;
                 }
                 int temp = pivot[max];
                 pivot[max] = pivot[col];
                 pivot[col] = temp;
                 even = !even;
             }
 
             // Divide the lower elements by the "winning" diagonal elt.
             final double luDiag = lu[col][col];
             for (int row = col + 1; row < m; row++) {
                 lu[row][col] /= luDiag;
             }
         }
     }
 
     /**
      * Returns the matrix L of the decomposition.
      * <p>L is a lower-triangular matrix</p>
      * @return the L matrix (or null if decomposed matrix is singular)
      */
     public RealMatrix getL() {
         if (cachedL == null && !singular) {
             final int m = pivot.length;
             cachedL = MatrixUtils.createRealMatrix(m, m);
             for (int i = 0; i < m; ++i) {
                 final double[] luI = lu[i];
                 for (int j = 0; j < i; ++j) {
                     cachedL.setEntry(i, j, luI[j]);
                 }
                 cachedL.setEntry(i, i, 1.0);
             }
         }
         return cachedL;
     }
 
     /**
      * Returns the matrix U of the decomposition.
      * <p>U is an upper-triangular matrix</p>
      * @return the U matrix (or null if decomposed matrix is singular)
      */
     public RealMatrix getU() {
         if (cachedU == null && !singular) {
             final int m = pivot.length;
             cachedU = MatrixUtils.createRealMatrix(m, m);
             for (int i = 0; i < m; ++i) {
                 final double[] luI = lu[i];
                 for (int j = i; j < m; ++j) {
                     cachedU.setEntry(i, j, luI[j]);
                 }
             }
         }
         return cachedU;
     }
 
     /**
      * Returns the P rows permutation matrix.
      * <p>P is a sparse matrix with exactly one element set to 1.0 in
      * each row and each column, all other elements being set to 0.0.</p>
      * <p>The positions of the 1 elements are given by the {@link #getPivot()
      * pivot permutation vector}.</p>
      * @return the P rows permutation matrix (or null if decomposed matrix is singular)
      * @see #getPivot()
      */
     public RealMatrix getP() {
         if (cachedP == null && !singular) {
             final int m = pivot.length;
             cachedP = MatrixUtils.createRealMatrix(m, m);
             for (int i = 0; i < m; ++i) {
                 cachedP.setEntry(i, pivot[i], 1.0);
             }
         }
         return cachedP;
     }
 
     /**
      * Returns the pivot permutation vector.
      * @return the pivot permutation vector
      * @see #getP()
      */
     public int[] getPivot() {
         return pivot.clone();
     }
 
     /**
      * Return the determinant of the matrix.
      * @return determinant of the matrix
      */
     public double getDeterminant() {
         if (singular) {
             return 0;
         } else {
             final int m = pivot.length;
             double determinant = even ? 1 : -1;
             for (int i = 0; i < m; i++) {
                 determinant *= lu[i][i];
             }
             return determinant;
         }
     }
 
     /**
      * Get a solver for finding the A &times; X = B solution in exact linear
      * sense.
      * @return a solver
      */
     public DecompositionSolver getSolver() {
         return new Solver(lu, pivot, singular);
     }
 
     /** Specialized solver. */
     private static final class Solver implements DecompositionSolver {
 
         /** Entries of LU decomposition. */
         private final double[][] lu;
 
         /** Pivot permutation associated with LU decomposition. */
         private final int[] pivot;
 
         /** Singularity indicator. */
         private final boolean singular;
 
         /**
          * Build a solver from decomposed matrix.
          * @param lu entries of LU decomposition
          * @param pivot pivot permutation associated with LU decomposition
          * @param singular singularity indicator
          */
         private Solver(final double[][] lu, final int[] pivot, final boolean singular) {
             this.lu       = lu;
             this.pivot    = pivot;
             this.singular = singular;
         }
 
         /** {@inheritDoc} */
         @Override
         public boolean isNonSingular() {
             return !singular;
         }
 
         /** {@inheritDoc} */
         @Override
         public RealVector solve(RealVector b) {
             final int m = pivot.length;
             if (b.getDimension() != m) {
                 throw new DimensionMismatchException(b.getDimension(), m);
             }
             if (singular) {
                 throw new SingularMatrixException();
             }
 
             final double[] bp = new double[m];
 
             // Apply permutations to b
             for (int row = 0; row < m; row++) {
                 bp[row] = b.getEntry(pivot[row]);
             }
 
             // Solve LY = b
             for (int col = 0; col < m; col++) {
                 final double bpCol = bp[col];
                 for (int i = col + 1; i < m; i++) {
                     bp[i] -= bpCol * lu[i][col];
                 }
             }
 
             // Solve UX = Y
             for (int col = m - 1; col >= 0; col--) {
                 bp[col] /= lu[col][col];
                 final double bpCol = bp[col];
                 for (int i = 0; i < col; i++) {
                     bp[i] -= bpCol * lu[i][col];
                 }
             }
 
             return new ArrayRealVector(bp, false);
         }
 
         /** {@inheritDoc} */
         @Override
         public RealMatrix solve(RealMatrix b) {
 
             final int m = pivot.length;
             if (b.getRowDimension() != m) {
                 throw new DimensionMismatchException(b.getRowDimension(), m);
             }
             if (singular) {
                 throw new SingularMatrixException();
             }
 
             final int nColB = b.getColumnDimension();
 
             // Apply permutations to b
             final double[][] bp = new double[m][nColB];
             for (int row = 0; row < m; row++) {
                 final double[] bpRow = bp[row];
                 final int pRow = pivot[row];
                 for (int col = 0; col < nColB; col++) {
                     bpRow[col] = b.getEntry(pRow, col);
                 }
             }
 
             // Solve LY = b
             for (int col = 0; col < m; col++) {
                 final double[] bpCol = bp[col];
                 for (int i = col + 1; i < m; i++) {
                     final double[] bpI = bp[i];
                     final double luICol = lu[i][col];
                     for (int j = 0; j < nColB; j++) {
                         bpI[j] -= bpCol[j] * luICol;
                     }
                 }
             }
 
             // Solve UX = Y
             for (int col = m - 1; col >= 0; col--) {
                 final double[] bpCol = bp[col];
                 final double luDiag = lu[col][col];
                 for (int j = 0; j < nColB; j++) {
                     bpCol[j] /= luDiag;
                 }
                 for (int i = 0; i < col; i++) {
                     final double[] bpI = bp[i];
                     final double luICol = lu[i][col];
                     for (int j = 0; j < nColB; j++) {
                         bpI[j] -= bpCol[j] * luICol;
                     }
                 }
             }
 
             return new Array2DRowRealMatrix(bp, false);
         }
 
         /**
          * Get the inverse of the decomposed matrix.
          *
          * @return the inverse matrix.
          * @throws SingularMatrixException if the decomposed matrix is singular.
          */
         @Override
         public RealMatrix getInverse() {
             return solve(MatrixUtils.createRealIdentityMatrix(pivot.length));
         }
     }
 }