/*
    * Copyright (C) 2012 Alberto Irurueta Carro (alberto@irurueta.com)
    *
    * Licensed 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 com.irurueta.algebra;
  
  /**
   * This class allows decomposition of matrices using LU decomposition, which
   * consists on retrieving two triangular matrices (lower triangular and upper
   * triangular) as a decomposition of provided input matrix.
   * In other words, if input matrix is A, then: A = L * U, where L is lower
   * triangular matrix and U is upper triangular matrix.
   * LU decomposition is a useful and fast way of solving systems of linear
   * equations, computing determinants or finding whether a matrix is singular.
   */
  @SuppressWarnings("DuplicatedCode")
  public class LUDecomposer extends Decomposer {
  
      /**
       * Constant defining default round error when determining singularity of
       * matrices. This value is zero by default.
       */
      public static final double DEFAULT_ROUND_ERROR = 0.0;
  
      /**
       * Constant defining minimum allowed round error value when determining
       * singularity of matrices.
       */
      public static final double MIN_ROUND_ERROR = 0.0;
  
      /**
       * Internal matrix containing results of decomposition.
       */
      private Matrix lu;
  
      /**
       * Internal array containing pivotings after decomposition.
       */
      int[] piv;
  
      /**
       * Member containing pivot sign after decomposition.
       */
      int pivSign;
  
      /**
       * Constructor of this class.
       */
      public LUDecomposer() {
          super();
          lu = null;
          piv = null;
      }
  
      /**
       * Constructor of this class.
       *
       * @param inputMatrix Reference to input matrix to be decomposed
       */
      public LUDecomposer(final Matrix inputMatrix) {
          super(inputMatrix);
          lu = null;
          piv = null;
      }
  
      /**
       * Returns decomposer type corresponding to LU decomposition
       *
       * @return Decomposer type
       */
      @Override
      public DecomposerType getDecomposerType() {
          return DecomposerType.LU_DECOMPOSITION;
      }
  
      /**
       * Sets reference to input matrix to be decomposed.
       *
       * @param inputMatrix Reference to input matrix to be decomposed.
       * @throws LockedException Exception thrown if attempting to call this
       *                         method while this instance remains locked.
       */
      @Override
      public void setInputMatrix(final Matrix inputMatrix) throws LockedException {
          super.setInputMatrix(inputMatrix);
          lu = null;
          piv = null;
      }
 
     /**
      * Returns boolean indicating whether decomposition has been computed and
      * results can be retrieved.
      * Attempting to retrieve decomposition results when not available, will
      * probably raise a NotAvailableException
      *
      * @return Boolean indicating whether decomposition has been computed and
      * results can be retrieved.
      */
     @Override
     public boolean isDecompositionAvailable() {
         return lu != null;
     }
 
     /**
      * This method computes LU matrix decomposition, which consists on
      * retrieving two triangular matrices (Lower triangular and Upper
      * triangular) as a decomposition of provided input matrix.
      * In other words, if input matrix is A, then: A = L * U
      * Note: During execution of this method, this instance will be locked,
      * and hence attempting to set some parameters might raise a
      * LockedException.
      * Note: After execution of this method, LU decomposition will be
      * available and operations such as retrieving L and U matrices or
      * computing determinants among others will be able to be done.
      * Attempting to call any of such operations before calling this method
      * will raise a NotAvailableException because they require computation of
      * LU decomposition first.
      *
      * @throws NotReadyException   Exception thrown if attempting to call this
      *                             method when this instance is not ready (i.e. no input matrix has been
      *                             provided).
      * @throws LockedException     Exception thrown if attempting to call this
      *                             method when this instance is not ready (i.e. no input matrix has been
      *                             provided).
      * @throws DecomposerException Exception thrown if for any reason
      *                             decomposition fails while executing, like when convergence of results
      *                             can not be obtained, etc.
      */
     @Override
     public void decompose() throws NotReadyException, LockedException, DecomposerException {
 
         if (!isReady()) {
             throw new NotReadyException();
         }
         if (isLocked()) {
             throw new LockedException();
         }
 
         final var rows = inputMatrix.getRows();
         final var columns = inputMatrix.getColumns();
 
         if (columns > rows) {
             throw new DecomposerException();
         }
 
         locked = true;
 
         // copy matrix contents
         lu = new Matrix(inputMatrix);
 
         piv = new int[rows];
         for (var i = 0; i < rows; i++) {
             piv[i] = i;
         }
 
         pivSign = 1;
 
         // Main loop
         for (var k = 0; k < columns; k++) {
             // Find pivot
             var p = k;
             for (var i = k + 1; i < rows; i++) {
                 p = Math.abs(lu.getElementAt(i, k)) > Math.abs(lu.getElementAt(p, k)) ? i : p;
             }
             // Exchange if necessary
             if (p != k) {
                 for (var j = 0; j < columns; j++) {
                     final var t = lu.getElementAt(p, j);
                     lu.setElementAt(p, j, lu.getElementAt(k, j));
                     lu.setElementAt(k, j, t);
                 }
                 final var t = piv[p];
                 piv[p] = piv[k];
                 piv[k] = t;
                 pivSign = -pivSign;
             }
             // Compute multipliers and eliminate k-th column
             if (lu.getElementAt(k, k) != 0.0) {
                 for (var i = k + 1; i < rows; i++) {
                     lu.setElementAt(i, k, lu.getElementAt(i, k) / lu.getElementAt(k, k));
                     for (var j = k + 1; j < columns; j++) {
                         lu.setElementAt(i, j, lu.getElementAt(i, j)
                                 - lu.getElementAt(i, k) * lu.getElementAt(k, j));
                     }
                 }
             }
         }
 
         locked = false;
     }
 
     /**
      * Return boolean indicating whether provided input matrix is singular
      * or not after computing LU decomposition. Returns true if singular and
      * false otherwise.
      * A matrix is defined as singular if its determinant is zero, hence
      * provided input matrix must be square (even though LU decomposition can be
      * computed for non-square matrices), otherwise a WrongSizeException will
      * be raised when calling this method. LU decomposition can be used to avoid
      * determinant computation by means of pivoting, because LU decomposition
      * obtains two triangular matrices, and the determinant of a triangular
      * matrix is just the product of the diagonal elements. Hence, if any
      * element on the diagonal of LU decomposition is zero, determinant will be
      * zero and input matrix will be singular.
      *
      * @return Boolean indicating whether provided input matrix is singular
      * or not.
      * @throws NotAvailableException    Exception thrown if attempting to call this
      *                                  method before computing LU decomposition. To avoid this exception call
      *                                  decompose() method first.
      * @throws WrongSizeException       Exception thrown if attempting to call this
      *                                  method using a non-square input matrix.
      * @throws IllegalArgumentException Exception thrown if provided rounding
      *                                  error is lower than minimum allowed value (MIN_ROUND_ERROR)
      * @see #decompose()
      */
     public boolean isSingular() throws NotAvailableException, WrongSizeException {
         return isSingular(DEFAULT_ROUND_ERROR);
     }
 
     /**
      * Return boolean indicating whether provided input matrix is singular
      * or not after computing LU decomposition. Returns true if singular and
      * false otherwise.
      * A matrix is defined as singular if its determinant is zero, hence
      * provided input matrix must be square (even though LU decomposition can be
      * computed for non-square matrices), otherwise a WrongSizeException will
      * be raised when calling this method. LU decomposition can be used to avoid
      * determinant computation by means of pivoting, because LU decomposition
      * obtains two triangular matrices, and the determinant of a triangular
      * matrix is just the product of the diagonal elements. Hence, if any
      * element on the diagonal of LU decomposition is zero, determinant will be
      * zero and input matrix will be singular.
      *
      * @param roundingError Determines the amount of margin given to determine
      *                      whether a matrix is singular or not due to rounding errors. If not
      *                      provided, by default rounding error is set to zero, but this value can be
      *                      relaxed if needed.
      * @return Boolean indicating whether provided input matrix is singular
      * or not.
      * @throws NotAvailableException    Exception thrown if attempting to call this
      *                                  method before computing LU decomposition. To avoid this exception call
      *                                  decompose() method first.
      * @throws WrongSizeException       Exception thrown if attempting to call this
      *                                  method using a non-square input matrix.
      * @throws IllegalArgumentException Exception thrown if provided rounding
      *                                  error is lower than minimum allowed value (MIN_ROUND_ERROR)
      * @see #decompose()
      */
     public boolean isSingular(final double roundingError) throws NotAvailableException, WrongSizeException {
 
         if (!isDecompositionAvailable()) {
             throw new NotAvailableException();
         }
         if (roundingError < MIN_ROUND_ERROR) {
             throw new IllegalArgumentException();
         }
 
         // A matrix is singular when its determinant is zero. Hence, in order to
         // compute singularity matrix must be square
         final var rows = inputMatrix.getRows();
         final var columns = inputMatrix.getColumns();
         if (rows != columns) {
             throw new WrongSizeException();
         }
 
         // Since we have computed LU decomposition into triangular matrices. The
         // determinant of a triangular matrix is the product of its diagonal
         // elements. Hence, if any element in the diagonal is zero, input matrix
         // will be singular.
         for (var j = 0; j < columns; j++) {
             if (lu.getElementAt(j, j) == 0.0) {
                 return true;
             }
         }
         return false;
     }
 
     /**
      * Fills provided matrix instance with the Lower triangular matrix
      * resulting from LU decomposition before correcting any possible pivots.
      * Hence, this matrix is only ensured to be Lower triangular.
      * In other words, this matrix does not ensure the product A = L * U, to
      * achieve this, we need to apply pivot correction.
      * A pivot corrected version of this matrix can be obtained by calling
      * method getL().
      *
      * @param pivottedL Lower triangular matrix.
      * @throws NotAvailableException Exception thrown if attempting to call
      *                               this method before computing LU decomposition. To avoid this exception
      *                               call decompose() method first.
      * @see #getL()
      * @see #decompose()
      */
     public void getPivottedL(final Matrix pivottedL) throws NotAvailableException {
 
         if (!isDecompositionAvailable()) {
             throw new NotAvailableException();
         }
 
         final var rows = lu.getRows();
         final var columns = lu.getColumns();
 
         if (pivottedL.getRows() != rows || pivottedL.getColumns() != columns) {
             try {
                 pivottedL.resize(rows, columns);
             } catch (final WrongSizeException ignore) {
                 // never happens
             }
         }
 
         for (var i = 0; i < rows; i++) {
             for (var j = 0; j < columns; j++) {
                 if (i > j) {
                     pivottedL.setElementAt(i, j, lu.getElementAt(i, j));
                 } else if (i == j) {
                     pivottedL.setElementAt(i, j, 1.0);
                 } else {
                     pivottedL.setElementAt(i, j, 0.0);
                 }
             }
         }
     }
 
     /**
      * Returns a new matrix instance containing the Lower triangular matrix
      * resulting from LU decomposition before correcting any possible pivots.
      * Hence, this matrix is only ensured to be Lower triangular.
      * In other words, this matrix does not ensure the product A = L * U, to
      * achieve this, we need to apply pivot correction.
      * A pivot corrected version of this matrix can be obtained by calling
      * method getL().
      *
      * @return Lower triangular matrix.
      * @throws NotAvailableException Exception thrown if attempting to call
      *                               this method before computing LU decomposition. To avoid this exception
      *                               call decompose() method first.
      * @see #getL()
      * @see #decompose()
      */
     public Matrix getPivottedL() throws NotAvailableException {
 
         if (!isDecompositionAvailable()) {
             throw new NotAvailableException();
         }
 
         final var rows = lu.getRows();
         final var columns = lu.getColumns();
 
         Matrix out = null;
         try {
             out = new Matrix(rows, columns);
         } catch (final WrongSizeException ignore) {
             // never happens
         }
         getPivottedL(out);
         return out;
     }
 
     /**
      * Fills provided matrix instance with the pivot corrected Lower
      * triangular matrix resulting from LU decomposition for provided input
      * matrix.
      * Since this matrix is pivot corrected, it might not be completely
      * triangular, except for some row pivotting.
      * Notice that LU decomposition obtains matrices in the form of
      * A = L * U, where A is provided input matrix, L is lower triangular
      * matrix and U is upper triangular matrix.
      *
      * @param l the Lower triangular matrix resulting from LU
      *          decomposition for provided input matrix.
      * @throws NotAvailableException Exception thrown if attempting to call
      *                               this method before computing LU decomposition. To avoid this exception
      *                               call decompose() method first.
      * @see #decompose()
      */
     public void getL(final Matrix l) throws NotAvailableException {
         if (!isDecompositionAvailable()) {
             throw new NotAvailableException();
         }
 
         final var rows = lu.getRows();
         final var columns = lu.getColumns();
 
         if (l.getRows() != rows || l.getColumns() != columns) {
             try {
                 l.resize(rows, columns);
             } catch (final WrongSizeException ignore) {
                 // never happens
             }
         }
 
         for (var i = 0; i < rows; i++) {
             for (var j = 0; j < columns; j++) {
                 if (i > j) {
                     l.setElementAt(piv[i], j, lu.getElementAt(i, j));
                 } else if (i == j) {
                     l.setElementAt(piv[i], j, 1.0);
                 } else {
                     l.setElementAt(piv[i], j, 0.0);
                 }
             }
         }
     }
 
     /**
      * Returns a new matrix instance containing the pivot corrected Lower
      * triangular matrix resulting from LU decomposition for provided input
      * matrix.
      * Since this matrix is pivot corrected, it might not be completely
      * triangular, except for some row pivotting.
      * Notice that LU decomposition obtains matrices in the form of
      * A = L * U, where A is provided input matrix, L is lower triangular
      * matrix and U is upper triangular matrix.
      *
      * @return Returns the Lower triangular matrix resulting from LU
      * decomposition for provided input matrix.
      * @throws NotAvailableException Exception thrown if attempting to call
      *                               this method before computing LU decomposition. To avoid this exception
      *                               call decompose() method first.
      * @see #decompose()
      */
     public Matrix getL() throws NotAvailableException {
         if (!isDecompositionAvailable()) {
             throw new NotAvailableException();
         }
 
         final var rows = lu.getRows();
         final var columns = lu.getColumns();
         Matrix out = null;
         try {
             out = new Matrix(rows, columns);
         } catch (final WrongSizeException ignore) {
             //never happens
         }
         getL(out);
         return out;
     }
 
     /**
      * Fills provided matrix instance with the Upper triangular matrix
      * resulting from LU decomposition for provided input matrix.
      * Notice that LU decomposition obtains matrices in the form A = L * U,
      * where A is provided input matrix, L is lower triangular matrix and U
      * is upper triangular matrix.
      *
      * @param u Returns the Upper triangular matrix resulting from LU
      *          decomposition for provided input matrix.
      * @throws NotAvailableException Exception thrown if attempting to call
      *                               this method before computing LU decomposition. To avoid this exception
      *                               call decompose() method first.
      * @see #decompose()
      */
     public void getU(final Matrix u) throws NotAvailableException {
         if (!isDecompositionAvailable()) {
             throw new NotAvailableException();
         }
 
         final var columns = lu.getColumns();
 
         if (u.getRows() != columns || u.getColumns() != columns) {
             try {
                 u.resize(columns, columns);
             } catch (final WrongSizeException ignore) {
                 // never happens
             }
         }
         for (var i = 0; i < columns; i++) {
             for (var j = 0; j < columns; j++) {
                 if (i <= j) {
                     u.setElementAt(i, j, lu.getElementAt(i, j));
                 } else {
                     u.setElementAt(i, j, 0.0);
                 }
             }
         }
     }
 
     /**
      * Returns a new matrix instance containing the Upper triangular matrix
      * resulting from LU decomposition for provided input matrix.
      * Notice that LU decomposition obtains matrices in the form A = L * U,
      * where A is provided input matrix, L is lower triangular matrix and U
      * is upper triangular matrix.
      *
      * @return Returns the Upper triangular matrix resulting from LU
      * decomposition for provided input matrix.
      * @throws NotAvailableException Exception thrown if attempting to call
      *                               this method before computing LU decomposition. To avoid this exception
      *                               call decompose() method first.
      * @see #decompose()
      */
     public Matrix getU() throws NotAvailableException {
         if (!isDecompositionAvailable()) {
             throw new NotAvailableException();
         }
 
         final var columns = lu.getColumns();
         Matrix out = null;
         try {
             out = new Matrix(columns, columns);
         } catch (final WrongSizeException ignore) {
             // never happens
         }
         getU(out);
         return out;
     }
 
     /**
      * Returns pivot permutation vector.
      *
      * @return Pivot permutation vector.
      * @throws NotAvailableException Exception thrown if attempting to call
      *                               this method before computing LU decomposition. To avoid this exception
      *                               call decompose() method first.
      */
     public int[] getPivot() throws NotAvailableException {
 
         if (!this.isDecompositionAvailable()) {
             throw new NotAvailableException();
         }
         return piv;
     }
 
     /**
      * Returns determinant of provided input matrix using LU decomposition as
      * means to obtain it.
      * Provided input matrix must be square (even though LU decomposition can
      * be computed for non-square matrices), otherwise a WrongSizeException will
      * be raised when calling this method.
      * LU decomposition can be used to avoid determinant computation using other
      * slow methods, because LU decomposition obtains two triangular matrices,
      * and the determinant of a triangular matrix is just the product of the
      * diagonal elements.
      * Since the determinant of a matrix product is the product of determinants,
      * then determinant of input matrix can be computed as the product of
      * determinants of L and U.
      * Finally, since L has ones on its diagonal, its determinant will be +-1,
      * depending on the amount of pivots done on L, and determinant of U will be
      * just the product of its diagonal elements.
      *
      * @return Determinant of provided input matrix.
      * @throws NotAvailableException Exception thrown if attempting to call
      *                               this method before computing LU decomposition. To avoid this exception
      *                               call decompose() method first.
      * @throws WrongSizeException    Exception thrown if attempting to call this
      *                               method using a non-square input matrix.
      * @see #decompose()
      */
     public double determinant() throws NotAvailableException, WrongSizeException {
 
         if (!isDecompositionAvailable()) {
             throw new NotAvailableException();
         }
 
         // Determinants can only be computed on squared matrices
         final var rows = inputMatrix.getRows();
         final var columns = inputMatrix.getColumns();
 
         if (rows != columns) {
             throw new WrongSizeException();
         }
 
         double d = pivSign;
         for (int j = 0; j < columns; j++) {
             d *= lu.getElementAt(j, j);
         }
 
         return d;
     }
 
 
     /**
      * Solves a linear system of equations of the following form:
      * A * X = B.
      * Where A is the input matrix provided for LU decomposition, X is the
      * solution to the system of equations, and B is the parameters
      * vector/matrix.
      * Note: This method can be reused for different b vectors/matrices without
      * having to recompute LU decomposition on the same input matrix.
      * Note: Provided b matrix must have the same number of rows as provided
      * input matrix A, otherwise a WrongSizeException will be raised.
      * Note: Provided input matrix A must be square, otherwise a
      * WrongSizeException will be raised as well.
      * Note: If provided input matrix A is singular, a SingularMatrixException
      * will be thrown.
      * Note: In order to execute this method, an LU decomposition must be
      * available, otherwise a NotAvailableException will be raised. In order
      * to avoid this exception call decompose() method first.
      * Note: DEFAULT_ROUND_ERROR is used as rounding error
      * Note: Solution of linear system of equations is stored in provided result
      * matrix
      *
      * @param b      Parameters matrix that determine a linear system of equations.
      *               Provided matrix must have the same number of rows as provided input
      *               matrix for LU decomposition. Besides, each column on parameters matrix
      *               will represent a new system of equations, whose solution will be returned
      *               on appropriate column as an output of this method.
      * @param result Matrix containing solution of linear system of equations on
      *               each column for each column of provided parameters matrix b.
      * @throws NotAvailableException    Exception thrown if attempting to call this
      *                                  method before computing LU decomposition. To avoid this exception call
      *                                  decompose() method first.
      * @throws WrongSizeException       Exception thrown if attempting to call this
      *                                  method using a non-square input matrix; or if provided parameters matrix
      *                                  (b) does not have the same number of rows as input matrix being LU
      *                                  decomposed.
      * @throws SingularMatrixException  Exception thrown if provided input matrix
      *                                  to be LU decomposed is singular. In this case linear system of equations
      *                                  cannot be solved.
      * @throws IllegalArgumentException Exception thrown if provided rounding
      *                                  error is lower than minimum allowed value (MIN_ROUND_ERROR).
      * @see #decompose()
      */
     public void solve(final Matrix b, final Matrix result) throws NotAvailableException, WrongSizeException,
             SingularMatrixException {
         solve(b, DEFAULT_ROUND_ERROR, result);
     }
 
     /**
      * Solves a linear system of equations of the following form:
      * A * X = B.
      * Where A is the input matrix provided for LU decomposition, X is the
      * solution to the system of equations, and B is the parameters
      * vector/matrix.
      * Note: This method can be reused for different b vectors/matrices without
      * having to recompute LU decomposition on the same input matrix.
      * Note: Provided b matrix must have the same number of rows as provided
      * input matrix A, otherwise a WrongSizeException will be raised.
      * Note: Provided input matrix A must be square, otherwise a
      * WrongSizeException will be raised as well.
      * Note: If provided input matrix A is singular, a SingularMatrixException
      * will be thrown.
      * Note: In order to execute this method, an LU decomposition must be
      * available, otherwise a NotAvailableException will be raised. In order
      * to avoid this exception call decompose() method first.
      * Note: Solution of linear system of equations is stored in provided result
      * matrix
      *
      * @param b             Parameters matrix that determine a linear system of equations.
      *                      Provided matrix must have the same number of rows as provided input
      *                      matrix for LU decomposition. Besides, each column on parameters matrix
      *                      will represent a new system of equations, whose solution will be returned
      *                      on appropriate column as an output of this method.
      * @param roundingError Determines the amount of margin given to determine
      *                      whether a matrix is singular or not due to rounding errors. If not
      *                      provided, by default rounding error is set to zero, but this value can
      *                      be relaxed if needed.
      * @param result        Matrix containing solution of linear system of equations on
      *                      each column for each column of provided parameters matrix b.
      * @throws NotAvailableException    Exception thrown if attempting to call this
      *                                  method before computing LU decomposition. To avoid this exception call
      *                                  decompose() method first.
      * @throws WrongSizeException       Exception thrown if attempting to call this
      *                                  method using a non-square input matrix; or if provided parameters matrix
      *                                  (b) does not have the same number of rows as input matrix being LU
      *                                  decomposed.
      * @throws SingularMatrixException  Exception thrown if provided input matrix
      *                                  to be LU decomposed is singular. In this case linear system of equations
      *                                  cannot be solved.
      * @throws IllegalArgumentException Exception thrown if provided rounding
      *                                  error is lower than minimum allowed value (MIN_ROUND_ERROR).
      * @see #decompose()
      */
     public void solve(final Matrix b, final double roundingError, final Matrix result)
             throws NotAvailableException, WrongSizeException, SingularMatrixException {
 
         if (!isDecompositionAvailable()) {
             throw new NotAvailableException();
         }
 
         if (b.getRows() != inputMatrix.getRows()) {
             throw new WrongSizeException();
         }
 
         if (roundingError < MIN_ROUND_ERROR) {
             throw new IllegalArgumentException();
         }
 
         // Copy right hand side with pivoting
         final var rows = lu.getRows();
         final var columns = lu.getColumns();
         final var colsB = b.getColumns();
 
         if (rows != columns) {
             throw new WrongSizeException();
         }
         if (isSingular(roundingError)) {
             throw new SingularMatrixException();
         }
 
         // resize result matrix if needed
         if (result.getRows() != columns || result.getColumns() != colsB) {
             result.resize(columns, colsB);
         }
         result.initialize(0.0);
 
         for (var i = 0; i < columns; i++) {
             for (var j = 0; j < colsB; j++) {
                 result.setElementAt(i, j, b.getElementAt(piv[i], j));
             }
         }
 
         // Solve L * Y = b(piv, :)
         for (var k = 0; k < columns; k++) {
             for (var i = k + 1; i < columns; i++) {
                 for (var j = 0; j < colsB; j++) {
                     result.setElementAt(i, j, result.getElementAt(i, j)
                             - result.getElementAt(k, j) * lu.getElementAt(i, k));
                 }
             }
         }
 
         // Solve U * X = Y
         for (var k = columns - 1; k >= 0; k--) {
             for (var j = 0; j < colsB; j++) {
                 result.setElementAt(k, j, result.getElementAt(k, j) / lu.getElementAt(k, k));
             }
             for (var i = 0; i < k; i++) {
                 for (var j = 0; j < colsB; j++) {
                     result.setElementAt(i, j, result.getElementAt(i, j)
                             - result.getElementAt(k, j) * lu.getElementAt(i, k));
                 }
             }
         }
     }
 
     /**
      * Solves a linear system of equations of the following form:
      * A * X = B.
      * Where A is the input matrix provided for LU decomposition, X is the
      * solution to the system of equations, and B is the parameters
      * vector/matrix.
      * Note: This method can be reused for different b vectors/matrices without
      * having to recompute LU decomposition on the same input matrix.
      * Note: Provided b matrix must have the same number of rows as provided
      * input matrix A, otherwise a WrongSizeException will be raised.
      * Note: Provided input matrix A must be square, otherwise a
      * WrongSizeException will be raised as well.
      * Note: If provided input matrix A is singular, a SingularMatrixException
      * will be thrown.
      * Note: In order to execute this method, an LU decomposition must be
      * available, otherwise a NotAvailableException will be raised. In order
      * to avoid this exception call decompose() method first.
      * Note: DEFAULT_ROUND_ERROR is used as rounding error
      *
      * @param b Parameters matrix that determine a linear system of equations.
      *          Provided matrix must have the same number of rows as provided input
      *          matrix for LU decomposition. Besides, each column on parameters matrix
      *          will represent a new system of equations, whose solution will be returned
      *          on appropriate column as an output of this method.
      * @return Matrix containing solution of linear system of equations on each
      * column for each column of provided parameters matrix b.
      * @throws NotAvailableException    Exception thrown if attempting to call this
      *                                  method before computing LU decomposition. To avoid this exception call
      *                                  decompose() method first.
      * @throws WrongSizeException       Exception thrown if attempting to call this
      *                                  method using a non-square input matrix; or if provided parameters matrix
      *                                  (b) does not have the same number of rows as input matrix being LU
      *                                  decomposed.
      * @throws SingularMatrixException  Exception thrown if provided input matrix
      *                                  to be LU decomposed is singular. In this case linear system of equations
      *                                  cannot be solved.
      * @throws IllegalArgumentException Exception thrown if provided rounding
      *                                  error is lower than minimum allowed value (MIN_ROUND_ERROR).
      * @see #decompose()
      */
     public Matrix solve(final Matrix b) throws NotAvailableException, WrongSizeException, SingularMatrixException {
         return solve(b, DEFAULT_ROUND_ERROR);
     }
 
     /**
      * Solves a linear system of equations of the following form:
      * A * X = B.
      * Where A is the input matrix provided for LU decomposition, X is the
      * solution to the system of equations, and B is the parameters
      * vector/matrix.
      * Note: This method can be reused for different b vectors/matrices without
      * having to recompute LU decomposition on the same input matrix.
      * Note: Provided b matrix must have the same number of rows as provided
      * input matrix A, otherwise a WrongSizeException will be raised.
      * Note: Provided input matrix A must be square, otherwise a
      * WrongSizeException will be raised as well.
      * Note: If provided input matrix A is singular, a SingularMatrixException
      * will be thrown.
      * Note: In order to execute this method, an LU decomposition must be
      * available, otherwise a NotAvailableException will be raised. In order
      * to avoid this exception call decompose() method first.
      *
      * @param b             Parameters matrix that determine a linear system of equations.
      *                      Provided matrix must have the same number of rows as provided input
      *                      matrix for LU decomposition. Besides, each column on parameters matrix
      *                      will represent a new system of equations, whose solution will be returned
      *                      on appropriate column as an output of this method.
      * @param roundingError Determines the amount of margin given to determine
      *                      whether a matrix is singular or not due to rounding errors. If not
      *                      provided, by default rounding error is set to zero, but this value can
      *                      be relaxed if needed.
      * @return Matrix containing solution of linear system of equations on each
      * column for each column of provided parameters matrix b.
      * @throws NotAvailableException    Exception thrown if attempting to call this
      *                                  method before computing LU decomposition. To avoid this exception call
      *                                  decompose() method first.
      * @throws WrongSizeException       Exception thrown if attempting to call this
      *                                  method using a non-square input matrix; or if provided parameters matrix
      *                                  (b) does not have the same number of rows as input matrix being LU
      *                                  decomposed.
      * @throws SingularMatrixException  Exception thrown if provided input matrix
      *                                  to be LU decomposed is singular. In this case linear system of equations
      *                                  cannot be solved.
      * @throws IllegalArgumentException Exception thrown if provided rounding
      *                                  error is lower than minimum allowed value (MIN_ROUND_ERROR).
      * @see #decompose()
      */
     public Matrix solve(final Matrix b, final double roundingError) throws NotAvailableException, WrongSizeException,
             SingularMatrixException {
 
         if (!isDecompositionAvailable()) {
             throw new NotAvailableException();
         }
 
         final var columns = lu.getColumns();
         final var colsB = b.getColumns();
         final var out = new Matrix(columns, colsB);
         solve(b, roundingError, out);
         return out;
     }
 }