001 /* 002 * Copyright (C) 2008-2010 by Holger Arndt 003 * 004 * This file is part of the Universal Java Matrix Package (UJMP). 005 * See the NOTICE file distributed with this work for additional 006 * information regarding copyright ownership and licensing. 007 * 008 * UJMP is free software; you can redistribute it and/or modify 009 * it under the terms of the GNU Lesser General Public License as 010 * published by the Free Software Foundation; either version 2 011 * of the License, or (at your option) any later version. 012 * 013 * UJMP is distributed in the hope that it will be useful, 014 * but WITHOUT ANY WARRANTY; without even the implied warranty of 015 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 016 * GNU Lesser General Public License for more details. 017 * 018 * You should have received a copy of the GNU Lesser General Public 019 * License along with UJMP; if not, write to the 020 * Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, 021 * Boston, MA 02110-1301 USA 022 */ 023 024 package org.ujmp.core.doublematrix.calculation.general.decomposition; 025 026 import org.ujmp.core.Matrix; 027 import org.ujmp.core.MatrixFactory; 028 import org.ujmp.core.util.DecompositionOps; 029 import org.ujmp.core.util.UJMPSettings; 030 031 /** 032 * Cholesky Decomposition. 033 * <P> 034 * For a symmetric, positive definite matrix A, the Cholesky decomposition is an 035 * lower triangular matrix L so that A = L*L'. 036 * <P> 037 * If the matrix is not symmetric or positive definite, the constructor returns 038 * a partial decomposition and sets an internal flag that may be queried by the 039 * isSPD() method. 040 */ 041 042 public interface Chol<T> { 043 044 public T calc(T source); 045 046 public T solve(T source, T b); 047 048 public static int THRESHOLD = 100; 049 050 public static final Chol<Matrix> MATRIX = new Chol<Matrix>() { 051 052 public final Matrix calc(Matrix source) { 053 if (UJMPSettings.getNumberOfThreads() == 1) { 054 if (source.getRowCount() >= THRESHOLD && source.getColumnCount() >= THRESHOLD) { 055 return MATRIXLARGESINGLETHREADED.calc(source); 056 } else { 057 return MATRIXSMALLSINGLETHREADED.calc(source); 058 } 059 } else { 060 if (source.getRowCount() >= THRESHOLD && source.getColumnCount() >= THRESHOLD) { 061 return MATRIXLARGEMULTITHREADED.calc(source); 062 } else { 063 return MATRIXSMALLMULTITHREADED.calc(source); 064 } 065 } 066 } 067 068 public final Matrix solve(Matrix source, Matrix b) { 069 if (UJMPSettings.getNumberOfThreads() == 1) { 070 if (source.getRowCount() >= THRESHOLD && source.getColumnCount() >= THRESHOLD) { 071 return MATRIXLARGESINGLETHREADED.solve(source, b); 072 } else { 073 return MATRIXSMALLSINGLETHREADED.solve(source, b); 074 } 075 } else { 076 if (source.getRowCount() >= THRESHOLD && source.getColumnCount() >= THRESHOLD) { 077 return MATRIXLARGEMULTITHREADED.solve(source, b); 078 } else { 079 return MATRIXSMALLMULTITHREADED.solve(source, b); 080 } 081 } 082 } 083 }; 084 085 public static final Chol<Matrix> INSTANCE = MATRIX; 086 087 public static final Chol<Matrix> UJMP = new Chol<Matrix>() { 088 089 public final Matrix calc(Matrix source) { 090 CholMatrix chol = new CholMatrix(source); 091 return chol.getL(); 092 } 093 094 public final Matrix solve(Matrix source, Matrix b) { 095 CholMatrix chol = new CholMatrix(source); 096 return chol.solve(b); 097 } 098 }; 099 100 public static final Chol<Matrix> MATRIXSMALLMULTITHREADED = UJMP; 101 102 public static final Chol<Matrix> MATRIXSMALLSINGLETHREADED = UJMP; 103 104 public static final Chol<Matrix> MATRIXLARGESINGLETHREADED = new Chol<Matrix>() { 105 public Matrix calc(Matrix source) { 106 Chol<Matrix> chol = null; 107 if (UJMPSettings.isUseJBlas()) { 108 chol = DecompositionOps.CHOL_JBLAS; 109 } 110 if (chol == null) { 111 chol = UJMP; 112 } 113 return chol.calc(source); 114 } 115 116 public final Matrix solve(Matrix source, Matrix b) { 117 Chol<Matrix> chol = null; 118 if (UJMPSettings.isUseJBlas()) { 119 chol = DecompositionOps.CHOL_JBLAS; 120 } 121 if (chol == null) { 122 chol = UJMP; 123 } 124 return chol.solve(source, b); 125 } 126 }; 127 128 public static final Chol<Matrix> MATRIXLARGEMULTITHREADED = new Chol<Matrix>() { 129 public Matrix calc(Matrix source) { 130 Chol<Matrix> chol = null; 131 if (UJMPSettings.isUseJBlas()) { 132 chol = DecompositionOps.CHOL_JBLAS; 133 } 134 if (chol == null && UJMPSettings.isUseOjalgo()) { 135 chol = DecompositionOps.CHOL_OJALGO; 136 } 137 if (chol == null) { 138 chol = UJMP; 139 } 140 return chol.calc(source); 141 } 142 143 public final Matrix solve(Matrix source, Matrix b) { 144 Chol<Matrix> chol = null; 145 if (UJMPSettings.isUseJBlas()) { 146 chol = DecompositionOps.CHOL_JBLAS; 147 } 148 if (chol == null && UJMPSettings.isUseOjalgo()) { 149 chol = DecompositionOps.CHOL_OJALGO; 150 } 151 if (chol == null) { 152 chol = UJMP; 153 } 154 return chol.solve(source, b); 155 } 156 }; 157 158 final class CholMatrix { 159 private static final long serialVersionUID = 400514872358216115L; 160 161 /** 162 * Array for internal storage of decomposition. 163 * 164 * @serial internal array storage. 165 */ 166 private final double[][] L; 167 168 /** 169 * Row and column dimension (square matrix). 170 * 171 * @serial matrix dimension. 172 */ 173 private final int n; 174 175 /** 176 * Symmetric and positive definite flag. 177 * 178 * @serial is symmetric and positive definite flag. 179 */ 180 private boolean isspd; 181 182 /* 183 * ------------------------ Constructor ------------------------ 184 */ 185 186 /** 187 * Cholesky algorithm for symmetric and positive definite matrix. 188 * 189 * @param A 190 * Square, symmetric matrix. 191 * @return Structure to access L and isspd flag. 192 */ 193 194 public CholMatrix(Matrix Arg) { 195 final double[][] A = Arg.toDoubleArray(); 196 n = (int) Arg.getRowCount(); 197 L = new double[n][n]; 198 isspd = (Arg.getColumnCount() == n); 199 // Main loop. 200 double[] Lrowj = null; 201 double[] Lrowk = null; 202 double[] Aj = null; 203 for (int j = 0; j < n; j++) { 204 Lrowj = L[j]; 205 Aj = A[j]; 206 double d = 0.0; 207 for (int k = 0; k < j; k++) { 208 Lrowk = L[k]; 209 double s = 0.0; 210 for (int i = 0; i < k; i++) { 211 s += Lrowk[i] * Lrowj[i]; 212 } 213 Lrowj[k] = s = (Aj[k] - s) / Lrowk[k]; 214 d = d + s * s; 215 isspd = isspd & (A[k][j] == Aj[k]); 216 } 217 d = Aj[j] - d; 218 isspd = isspd & (d > 0.0); 219 Lrowj[j] = Math.sqrt(Math.max(d, 0.0)); 220 for (int k = j + 1; k < n; k++) { 221 Lrowj[k] = 0.0; 222 } 223 } 224 } 225 226 /* 227 * ------------------------ Temporary, experimental code. 228 * ------------------------ *\ 229 * 230 * \** Right Triangular Cholesky Decomposition. <P> For a symmetric, 231 * positive definite matrix A, the Right Cholesky decomposition is an 232 * upper triangular matrix R so that A = R'*R. This constructor computes 233 * R with the Fortran inspired column oriented algorithm used in LINPACK 234 * and MATLAB. In Java, we suspect a row oriented, lower triangular 235 * decomposition is faster. We have temporarily included this 236 * constructor here until timing experiments confirm this suspicion.\ 237 * 238 * \** Array for internal storage of right triangular decomposition. **\ 239 * private transient double[][] R; 240 * 241 * \** Cholesky algorithm for symmetric and positive definite matrix. 242 * 243 * @param A Square, symmetric matrix. 244 * 245 * @param rightflag Actual value ignored. 246 * 247 * @return Structure to access R and isspd flag.\ 248 * 249 * public CholeskyDecomposition (Matrix Arg, int rightflag) { // 250 * Initialize. double[][] A = Arg.getArray(); n = 251 * Arg.getColumnDimension(); R = new double[n][n]; isspd = 252 * (Arg.getColumnDimension() == n); // Main loop. for (int j = 0; j < n; 253 * j++) { double d = 0.0; for (int k = 0; k < j; k++) { double s = 254 * A[k][j]; for (int i = 0; i < k; i++) { s = s - R[i][k]*R[i][j]; } 255 * R[k][j] = s = s/R[k][k]; d = d + s*s; isspd = isspd & (A[k][j] == 256 * A[j][k]); } d = A[j][j] - d; isspd = isspd & (d > 0.0); R[j][j] = 257 * Math.sqrt(Math.max(d,0.0)); for (int k = j+1; k < n; k++) { R[k][j] = 258 * 0.0; } } } 259 * 260 * \** Return upper triangular factor. 261 * 262 * @return R\ 263 * 264 * public Matrix getR () { return new Matrix(R,n,n); } 265 * 266 * \* ------------------------ End of temporary code. 267 * ------------------------ 268 */ 269 270 /* 271 * ------------------------ Public Methods ------------------------ 272 */ 273 274 /** 275 * Is the matrix symmetric and positive definite? 276 * 277 * @return true if A is symmetric and positive definite. 278 */ 279 280 public final boolean isSPD() { 281 return isspd; 282 } 283 284 /** 285 * Return triangular factor. 286 * 287 * @return L 288 */ 289 290 public final Matrix getL() { 291 return MatrixFactory.linkToArray(L); 292 } 293 294 /** 295 * Solve A*X = B 296 * 297 * @param B 298 * A Matrix with as many rows as A and any number of columns. 299 * @return X so that L*L'*X = B 300 * @exception IllegalArgumentException 301 * Matrix row dimensions must agree. 302 * @exception RuntimeException 303 * Matrix is not symmetric positive definite. 304 */ 305 306 public final Matrix solve(Matrix B) { 307 if (B.getRowCount() != n) { 308 throw new IllegalArgumentException("Matrix row dimensions must agree."); 309 } 310 if (!isspd) { 311 throw new RuntimeException("Matrix is not symmetric positive definite."); 312 } 313 314 // Copy right hand side. 315 final double[][] X = B.toDoubleArray(); 316 final int nx = (int) B.getColumnCount(); 317 318 // Solve L*Y = B; 319 for (int k = 0; k < n; k++) { 320 for (int j = 0; j < nx; j++) { 321 for (int i = 0; i < k; i++) { 322 X[k][j] -= X[i][j] * L[k][i]; 323 } 324 X[k][j] /= L[k][k]; 325 } 326 } 327 328 // Solve L'*X = Y; 329 for (int k = n - 1; k >= 0; k--) { 330 for (int j = 0; j < nx; j++) { 331 for (int i = k + 1; i < n; i++) { 332 X[k][j] -= X[i][j] * L[i][k]; 333 } 334 X[k][j] /= L[k][k]; 335 } 336 } 337 338 return MatrixFactory.linkToArray(X); 339 } 340 341 } 342 }