001 package org.ujmp.core.doublematrix.calculation.general.decomposition; 002 003 import org.ujmp.core.Matrix; 004 import org.ujmp.core.MatrixFactory; 005 import org.ujmp.core.doublematrix.DenseDoubleMatrix2D; 006 import org.ujmp.core.util.DecompositionOps; 007 import org.ujmp.core.util.UJMPSettings; 008 009 /** 010 * LU Decomposition. 011 * <P> 012 * For an m-by-n matrix A with m >= n, the LU decomposition is an m-by-n unit 013 * lower triangular matrix L, an n-by-n upper triangular matrix U, and a 014 * permutation vector piv of length m so that A(piv,:) = L*U. If m < n, then L 015 * is m-by-m and U is m-by-n. 016 * <P> 017 * The LU decompostion with pivoting always exists, even if the matrix is 018 * singular, so the constructor will never fail. The primary use of the LU 019 * decomposition is in the solution of square systems of simultaneous linear 020 * equations. This will fail if isNonsingular() returns false. 021 */ 022 023 public interface LU<T> { 024 025 public T[] calc(T source); 026 027 public T solve(T source, T b); 028 029 public static int THRESHOLD = 100; 030 031 public static final LU<Matrix> MATRIX = new LU<Matrix>() { 032 033 public final Matrix[] calc(Matrix source) { 034 if (UJMPSettings.getNumberOfThreads() == 1) { 035 if (source.getRowCount() >= THRESHOLD && source.getColumnCount() >= THRESHOLD) { 036 return MATRIXLARGESINGLETHREADED.calc(source); 037 } else { 038 return MATRIXSMALLSINGLETHREADED.calc(source); 039 } 040 } else { 041 if (source.getRowCount() >= THRESHOLD && source.getColumnCount() >= THRESHOLD) { 042 return MATRIXLARGEMULTITHREADED.calc(source); 043 } else { 044 return MATRIXSMALLMULTITHREADED.calc(source); 045 } 046 } 047 } 048 049 public final Matrix solve(Matrix source, Matrix b) { 050 if (UJMPSettings.getNumberOfThreads() == 1) { 051 if (source.getRowCount() >= THRESHOLD && source.getColumnCount() >= THRESHOLD) { 052 return MATRIXLARGESINGLETHREADED.solve(source, b); 053 } else { 054 return MATRIXSMALLSINGLETHREADED.solve(source, b); 055 } 056 } else { 057 if (source.getRowCount() >= THRESHOLD && source.getColumnCount() >= THRESHOLD) { 058 return MATRIXLARGEMULTITHREADED.solve(source, b); 059 } else { 060 return MATRIXSMALLMULTITHREADED.solve(source, b); 061 } 062 } 063 } 064 }; 065 066 public static final LU<Matrix> INSTANCE = MATRIX; 067 068 public static final LU<Matrix> UJMP = new LU<Matrix>() { 069 070 public final Matrix[] calc(Matrix source) { 071 LUMatrix lu = new LUMatrix(source); 072 return new Matrix[] { lu.getL(), lu.getU(), lu.getP() }; 073 } 074 075 public final Matrix solve(Matrix source, Matrix b) { 076 LUMatrix lu = new LUMatrix(source); 077 return lu.solve(b); 078 } 079 }; 080 081 public static final LU<Matrix> MATRIXSMALLMULTITHREADED = UJMP; 082 083 public static final LU<Matrix> MATRIXSMALLSINGLETHREADED = UJMP; 084 085 public static final LU<Matrix> MATRIXLARGESINGLETHREADED = new LU<Matrix>() { 086 public final Matrix[] calc(Matrix source) { 087 LU<Matrix> lu = null; 088 if (UJMPSettings.isUseJBlas()) { 089 lu = DecompositionOps.LU_JBLAS; 090 } 091 if (lu == null) { 092 lu = UJMP; 093 } 094 return lu.calc(source); 095 } 096 097 public final Matrix solve(Matrix source, Matrix b) { 098 LU<Matrix> lu = null; 099 if (UJMPSettings.isUseJBlas()) { 100 lu = DecompositionOps.LU_JBLAS; 101 } 102 if (lu == null && UJMPSettings.isUseOjalgo()) { 103 lu = DecompositionOps.LU_OJALGO; 104 } 105 if (lu == null) { 106 lu = UJMP; 107 } 108 return lu.solve(source, b); 109 } 110 }; 111 112 public static final LU<Matrix> MATRIXLARGEMULTITHREADED = new LU<Matrix>() { 113 public final Matrix[] calc(Matrix source) { 114 LU<Matrix> lu = null; 115 if (UJMPSettings.isUseJBlas()) { 116 lu = DecompositionOps.LU_JBLAS; 117 } 118 if (lu == null && UJMPSettings.isUseOjalgo()) { 119 lu = DecompositionOps.LU_OJALGO; 120 } 121 if (lu == null) { 122 lu = UJMP; 123 } 124 return lu.calc(source); 125 } 126 127 public final Matrix solve(Matrix source, Matrix b) { 128 LU<Matrix> lu = null; 129 if (UJMPSettings.isUseJBlas()) { 130 lu = DecompositionOps.LU_JBLAS; 131 } 132 if (lu == null && UJMPSettings.isUseOjalgo()) { 133 lu = DecompositionOps.LU_OJALGO; 134 } 135 if (lu == null) { 136 lu = UJMP; 137 } 138 return lu.solve(source, b); 139 } 140 }; 141 142 final class LUMatrix { 143 /** 144 * Array for internal storage of decomposition. 145 * 146 * @serial internal array storage. 147 */ 148 private final double[][] LU; 149 150 /** 151 * Row and column dimensions, and pivot sign. 152 * 153 * @serial column dimension. 154 * @serial row dimension. 155 * @serial pivot sign. 156 */ 157 private final int m, n; 158 private int pivsign; 159 160 /** 161 * Internal storage of pivot vector. 162 * 163 * @serial pivot vector. 164 */ 165 private final int[] piv; 166 167 /* 168 * ------------------------ Constructor ------------------------ 169 */ 170 171 /** 172 * LU Decomposition 173 * 174 * @param A 175 * Rectangular matrix 176 * @return Structure to access L, U and piv. 177 */ 178 public LUMatrix(Matrix A) { 179 180 // Use a "left-looking", dot-product, Crout/Doolittle algorithm. 181 182 LU = A.toDoubleArray(); 183 m = (int) A.getRowCount(); 184 n = (int) A.getColumnCount(); 185 piv = new int[m]; 186 for (int i = 0; i < m; i++) { 187 piv[i] = i; 188 } 189 pivsign = 1; 190 191 final double[] LUcolj = new double[m]; 192 double[] LUrowi = null; 193 194 // Outer loop. 195 196 for (int j = 0; j < n; j++) { 197 198 // Make a copy of the j-th column to localize references. 199 for (int i = 0; i < m; i++) { 200 LUcolj[i] = LU[i][j]; 201 } 202 203 // Apply previous transformations. 204 205 for (int i = 0; i < m; i++) { 206 LUrowi = LU[i]; 207 208 // Most of the time is spent in the following dot product. 209 210 final int kmax = Math.min(i, j); 211 double s = 0.0; 212 for (int k = 0; k < kmax; k++) { 213 s += LUrowi[k] * LUcolj[k]; 214 } 215 216 LUrowi[j] = LUcolj[i] -= s; 217 } 218 219 int p = j; 220 for (int i = j + 1; i < m; i++) { 221 if (Math.abs(LUcolj[i]) > Math.abs(LUcolj[p])) { 222 p = i; 223 } 224 } 225 if (p != j) { 226 for (int k = 0; k < n; k++) { 227 double t = LU[p][k]; 228 LU[p][k] = LU[j][k]; 229 LU[j][k] = t; 230 } 231 int k = piv[p]; 232 piv[p] = piv[j]; 233 piv[j] = k; 234 pivsign = -pivsign; 235 } 236 237 // Compute multipliers. 238 239 // http://cio.nist.gov/esd/emaildir/lists/jama/msg01498.html 240 if (j < m && LU[j][j] != 0.0) { 241 for (int i = j + 1; i < m; i++) { 242 LU[i][j] /= LU[j][j]; 243 } 244 } 245 } 246 } 247 248 /* 249 * ------------------------ Temporary, experimental code. 250 * ------------------------ *\ 251 * 252 * \** LU Decomposition, computed by Gaussian elimination. <P> This 253 * constructor computes L and U with the "daxpy"-based elimination 254 * algorithm used in LINPACK and MATLAB. In Java, we suspect the 255 * dot-product, Crout algorithm will be faster. We have temporarily 256 * included this constructor until timing experiments confirm this 257 * suspicion. <P> 258 * 259 * @param A Rectangular matrix 260 * 261 * @param linpackflag Use Gaussian elimination. Actual value ignored. 262 * 263 * @return Structure to access L, U and piv.\ 264 * 265 * public LUDecomposition (Matrix A, int linpackflag) { // Initialize. 266 * LU = A.getArrayCopy(); m = A.getRowDimension(); n = 267 * A.getColumnDimension(); piv = new int[m]; for (int i = 0; i < m; i++) 268 * { piv[i] = i; } pivsign = 1; // Main loop. for (int k = 0; k < n; 269 * k++) { // Find pivot. int p = k; for (int i = k+1; i < m; i++) { if 270 * (Math.abs(LU[i][k]) > Math.abs(LU[p][k])) { p = i; } } // Exchange if 271 * necessary. if (p != k) { for (int j = 0; j < n; j++) { double t = 272 * LU[p][j]; LU[p][j] = LU[k][j]; LU[k][j] = t; } int t = piv[p]; piv[p] 273 * = piv[k]; piv[k] = t; pivsign = -pivsign; } // Compute multipliers 274 * and eliminate k-th column. if (LU[k][k] != 0.0) { for (int i = k+1; i 275 * < m; i++) { LU[i][k] /= LU[k][k]; for (int j = k+1; j < n; j++) { 276 * LU[i][j] -= LU[i][k]*LU[k][j]; } } } } } 277 * 278 * \* ------------------------ End of temporary code. 279 * ------------------------ 280 */ 281 282 /* 283 * ------------------------ Public Methods ------------------------ 284 */ 285 286 /** 287 * Is the matrix nonsingular? 288 * 289 * @return true if U, and hence A, is nonsingular. 290 */ 291 292 public final boolean isNonsingular() { 293 for (int j = 0; j < n; j++) { 294 if (LU[j][j] == 0) 295 return false; 296 } 297 return true; 298 } 299 300 /** 301 * Return lower triangular factor 302 * 303 * @return L 304 */ 305 306 public final Matrix getL() { 307 final int min = Math.min(m, n); 308 final double[][] L = new double[m][min]; 309 for (int i = 0; i < m; i++) { 310 for (int j = 0; j < min; j++) { 311 if (i > j) { 312 L[i][j] = LU[i][j]; 313 } else if (i == j) { 314 L[i][j] = 1.0; 315 } 316 } 317 } 318 return MatrixFactory.linkToArray(L); 319 } 320 321 /** 322 * Return upper triangular factor 323 * 324 * @return U 325 */ 326 327 public final Matrix getU() { 328 final int min = Math.min(m, n); 329 final double[][] U = new double[min][n]; 330 for (int i = 0; i < min; i++) { 331 for (int j = 0; j < n; j++) { 332 if (i <= j) { 333 U[i][j] = LU[i][j]; 334 } 335 } 336 } 337 return MatrixFactory.linkToArray(U); 338 } 339 340 /** 341 * Return pivot permutation vector 342 * 343 * @return piv 344 */ 345 346 public final int[] getPivot() { 347 final int[] p = new int[m]; 348 for (int i = 0; i < m; i++) { 349 p[i] = piv[i]; 350 } 351 return p; 352 } 353 354 public final Matrix getP() { 355 final DenseDoubleMatrix2D p = DenseDoubleMatrix2D.factory.zeros(m, m); 356 for (int i = 0; i < m; i++) { 357 p.setDouble(1, i, piv[i]); 358 } 359 return p; 360 } 361 362 /** 363 * Return pivot permutation vector as a one-dimensional double array 364 * 365 * @return (double) piv 366 */ 367 368 public final double[] getDoublePivot() { 369 final double[] vals = new double[m]; 370 for (int i = 0; i < m; i++) { 371 vals[i] = (double) piv[i]; 372 } 373 return vals; 374 } 375 376 /** 377 * Determinant 378 * 379 * @return det(A) 380 * @exception IllegalArgumentException 381 * Matrix must be square 382 */ 383 384 public final double det() { 385 if (m != n) { 386 throw new IllegalArgumentException("Matrix must be square."); 387 } 388 double d = (double) pivsign; 389 for (int j = 0; j < n; j++) { 390 d *= LU[j][j]; 391 } 392 return d; 393 } 394 395 /** 396 * Solve A*X = B 397 * 398 * @param B 399 * A Matrix with as many rows as A and any number of columns. 400 * @return X so that L*U*X = B(piv,:) 401 * @exception IllegalArgumentException 402 * Matrix row dimensions must agree. 403 * @exception RuntimeException 404 * Matrix is singular. 405 */ 406 407 public final Matrix solve(Matrix B) { 408 if (B.getRowCount() != m) { 409 throw new IllegalArgumentException("Matrix row dimensions must agree."); 410 } 411 if (!this.isNonsingular()) { 412 throw new RuntimeException("Matrix is singular."); 413 } 414 415 // Copy right hand side with pivoting 416 final int nx = (int) B.getColumnCount(); 417 // final Matrix Xmat = B.selectRows(Ret.NEW, 418 // MathUtil.toLongArray(piv)); 419 420 final double[][] X = new double[piv.length][(int) B.getColumnCount()]; 421 if (B instanceof DenseDoubleMatrix2D) { 422 DenseDoubleMatrix2D m = (DenseDoubleMatrix2D) B; 423 for (int c = (int) B.getColumnCount(); --c >= 0;) { 424 for (int r = piv.length; --r >= 0;) { 425 X[r][c] = m.getDouble(piv[r], c); 426 } 427 } 428 } else { 429 for (int c = (int) B.getColumnCount(); --c >= 0;) { 430 for (int r = piv.length; --r >= 0;) { 431 X[r][c] = B.getAsDouble(piv[r], c); 432 } 433 } 434 } 435 436 // Solve L*Y = B(piv,:) 437 for (int k = 0; k < n; k++) { 438 for (int i = k + 1; i < n; i++) { 439 for (int j = 0; j < nx; j++) { 440 X[i][j] -= X[k][j] * LU[i][k]; 441 } 442 } 443 } 444 445 for (int k = n - 1; k >= 0; k--) { 446 for (int j = 0; j < nx; j++) { 447 X[k][j] /= LU[k][k]; 448 } 449 for (int i = 0; i < k; i++) { 450 for (int j = 0; j < nx; j++) { 451 X[i][j] -= X[k][j] * LU[i][k]; 452 } 453 } 454 } 455 return MatrixFactory.linkToArray(X); 456 } 457 458 }; 459 }