Laboratoire ICAR » Plateforme TXM

package JamaPlus;

   /** Cholesky Decomposition.

   <P>

   For a symmetric, positive definite matrix A, the Cholesky decomposition

   is an lower triangular matrix L so that A = L*L'.

   <P>

   If the matrix is not symmetric or positive definite, the constructor

   returns a partial decomposition and sets an internal flag that may

   be queried by the isSPD() method.

   */

public class CholeskyDecomposition implements java.io.Serializable {

/* ------------------------

   Class variables

 * ------------------------ */

   /** Array for internal storage of decomposition.

   @serial internal array storage.

   */

   private double[][] L;

   /** Row and column dimension (square matrix).

   @serial matrix dimension.

   */

   private int n;

   /** Symmetric and positive definite flag.

   @serial is symmetric and positive definite flag.

   */

   private boolean isspd;

/* ------------------------

   Constructor

 * ------------------------ */

   /** Cholesky algorithm for symmetric and positive definite matrix.

   @param  A   Square, symmetric matrix.

   @return     Structure to access L and isspd flag.

   */

   public CholeskyDecomposition (Matrix Arg) {

     // Initialize.

      double[][] A = Arg.getArray();

      n = Arg.getRowDimension();

      L = new double[n][n];

      isspd = (Arg.getColumnDimension() == n);

      // Main loop.

      for (int j = 0; j < n; j++) {

         double[] Lrowj = L[j];

         double d = 0.0;

         for (int k = 0; k < j; k++) {

            double[] Lrowk = L[k];

            double s = 0.0;

            for (int i = 0; i < k; i++) {

               s += Lrowk[i]*Lrowj[i];

            }

            Lrowj[k] = s = (A[j][k] - s)/L[k][k];

            d = d + s*s;

            isspd = isspd & (A[k][j] == A[j][k]);

         }

         d = A[j][j] - d;

         isspd = isspd & (d > 0.0);

         L[j][j] = Math.sqrt(Math.max(d,0.0));

         for (int k = j+1; k < n; k++) {

            L[j][k] = 0.0;

         }

      }

   }

/* ------------------------

   Temporary, experimental code.

 * ------------------------ *\

   \** Right Triangular Cholesky Decomposition.

   <P>

   For a symmetric, positive definite matrix A, the Right Cholesky

   decomposition is an upper triangular matrix R so that A = R'*R.

   This constructor computes R with the Fortran inspired column oriented

   algorithm used in LINPACK and MATLAB.  In Java, we suspect a row oriented,

   lower triangular decomposition is faster.  We have temporarily included

   this constructor here until timing experiments confirm this suspicion.

   *\

   \** Array for internal storage of right triangular decomposition. **\

   private transient double[][] R;

   \** Cholesky algorithm for symmetric and positive definite matrix.

   @param  A           Square, symmetric matrix.

   @param  rightflag   Actual value ignored.

   @return             Structure to access R and isspd flag.

   *\

   public CholeskyDecomposition (Matrix Arg, int rightflag) {

      // Initialize.

      double[][] A = Arg.getArray();

      n = Arg.getColumnDimension();

      R = new double[n][n];

      isspd = (Arg.getColumnDimension() == n);

      // Main loop.

      for (int j = 0; j < n; j++) {

         double d = 0.0;

         for (int k = 0; k < j; k++) {

            double s = A[k][j];

            for (int i = 0; i < k; i++) {

               s = s - R[i][k]*R[i][j];

            }

            R[k][j] = s = s/R[k][k];

            d = d + s*s;

            isspd = isspd & (A[k][j] == A[j][k]);

         }

         d = A[j][j] - d;

         isspd = isspd & (d > 0.0);

         R[j][j] = Math.sqrt(Math.max(d,0.0));

         for (int k = j+1; k < n; k++) {

            R[k][j] = 0.0;

         }

      }

   }

   \** Return upper triangular factor.

   @return     R

   *\

   public Matrix getR () {

      return new Matrix(R,n,n);

   }

\* ------------------------

   End of temporary code.

 * ------------------------ */

/* ------------------------

   Public Methods

 * ------------------------ */

   /** Is the matrix symmetric and positive definite?

   @return     true if A is symmetric and positive definite.

   */

   public boolean isSPD () {

      return isspd;

   }

   /** Return triangular factor.

   @return     L

   */

   public Matrix getL () {

      return new Matrix(L,n,n);

   }

   /** Solve A*X = B

   @param  B   A Matrix with as many rows as A and any number of columns.

   @return     X so that L*L'*X = B

   @exception  IllegalArgumentException  Matrix row dimensions must agree.

   @exception  RuntimeException  Matrix is not symmetric positive definite.

   */

   public Matrix solve (Matrix B) {

      if (B.getRowDimension() != n) {

         throw new IllegalArgumentException("Matrix row dimensions must agree.");

      }

      if (!isspd) {

         throw new RuntimeException("Matrix is not symmetric positive definite.");

      }

      // Copy right hand side.

      double[][] X = B.getArrayCopy();

      int nx = B.getColumnDimension();

              // Solve L*Y = B;

              for (int k = 0; k < n; k++) {

                for (int j = 0; j < nx; j++) {

                   for (int i = 0; i < k ; i++) {

                       X[k][j] -= X[i][j]*L[k][i];

                   }

                   X[k][j] /= L[k][k];

                }

              }

              // Solve L'*X = Y;

              for (int k = n-1; k >= 0; k--) {

                for (int j = 0; j < nx; j++) {

                   for (int i = k+1; i < n ; i++) {

                       X[k][j] -= X[i][j]*L[i][k];

                   }

                   X[k][j] /= L[k][k];

                }

              }

      return new Matrix(X,n,nx);

   }

}