edu.mines.jtk.la.DMatrixLud Maven / Gradle / Ivy

Go to download
Show more of this group Show more artifacts with this name
Show all versions of edu-mines-jtk Show documentation
Java packages for science and engineering
The newest version!
/****************************************************************************
Copyright 2006, Colorado School of Mines and others.
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 edu.mines.jtk.la;

import static java.lang.Math.abs;
import static java.lang.Math.min;

import edu.mines.jtk.util.Check;

/**
 * LU decomposition (with pivoting) of a matrix A. 
 * For an m-by-n matrix A, with m>=n, the LU decomposition is 
 * A(piv,:) = L*U, where L is an m-by-n unit lower triangular matrix, 
 * U is an n-by-n upper-triangular matrix, and piv is a permutation
 * vector of length m.
 * 
 * The LU decomposition with pivoting always exists, even for singular
 * matrices A. The primary use of LU decomposition is in the solution of
 * square systems of simultaneous linear equations. These solutions will
 * fila if the matrix A is singular.
 * 

 * This class was adapted from the package Jama, which was developed by 
 * Joe Hicklin, Cleve Moler, and Peter Webb of The MathWorks, Inc., and by
 * Ronald Boisvert, Bruce Miller, Roldan Pozo, and Karin Remington of the
 * National Institue of Standards and Technology.
 * @author Dave Hale, Colorado School of Mines
 * @version 2006.09.15
 */
public class DMatrixLud {

  /** 
   * Constructs an LU decomposition for the specified matrix A.
   * @param a the matrix A.
   */
  public DMatrixLud(DMatrix a) {
    int m = _m = a.getM();
    int n = _n = a.getN();
    double[][] lu = _lu = a.get();
    _piv = new int[m];
    for (int i=0; iabs(lucolj[p]))
          p = i;
      }
      if (p!=j) {
        for (int k=0; kj) {
          l[i][j] = _lu[i][j];
        } else if (i==j) {
          l[i][j] = 1.0;
        } else {
          l[i][j] = 0.0;
        }
      }
    }
    return new DMatrix(_m,_n,l);
  }

  /** 
   * Gets the n-by-n upper triangular matrix factor U.
   * @return the n-by-n matrix factor U.
   */
  public DMatrix getU() {
    double[][] u = new double[_n][_n];
    for (int i=0; i<_n; ++i) {
      for (int j=0; j<_n; ++j) {
        if (i<=j) {
          u[i][j] = _lu[i][j];
        } else {
          u[i][j] = 0.0;
        }
      }
    }
    return new DMatrix(_n,_n,u);
  }

  /**
   * Gets the pivot vector, an array of length m.
   * @return the pivot vector.
   */
  public int[] getPivot() {
    int[] p = new int[_m];
    for (int i=0; i<_m; ++i)
      p[i] = _piv[i];
    return p;
  }

  /**
   * Returns the solution X of the system A*X = B.
   * This solution exists only if the matrix A is non-singular.
   * @param b a matrix of right-hand-side vectors. This matrix must
   *  have the same number (m) of rows as the matrix A, but may have
   *  any number of columns.
   * @return the matrix solution X.
   * @exception IllegalStateException if A is singular.
   */
  public DMatrix solve(DMatrix b) {
    Check.argument(b.getM()==_m,"A and B have the same number of rows");
    Check.state(this.isNonSingular(),"A is non-singular");

    // Copy of right-hand side with pivoting.
    int nx = b.getN();
    DMatrix xx = b.get(_piv,0,nx-1);
    double[][] x = xx.getArray();

    // Solve L*Y = B(piv,:).
    for (int k=0; k<_n; ++k) {
      for (int i=k+1; i<_n; ++i) {
        for (int j=0; j=0; --k) {
      for (int j=0; j