1 /* 2 * Licensed to the Apache Software Foundation (ASF) under one or more 3 * contributor license agreements. See the NOTICE file distributed with 4 * this work for additional information regarding copyright ownership. 5 * The ASF licenses this file to You under the Apache License, Version 2.0 6 * (the "License"); you may not use this file except in compliance with 7 * the License. You may obtain a copy of the License at 8 * 9 * http://www.apache.org/licenses/LICENSE-2.0 10 * 11 * Unless required by applicable law or agreed to in writing, software 12 * distributed under the License is distributed on an "AS IS" BASIS, 13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 14 * See the License for the specific language governing permissions and 15 * limitations under the License. 16 */ 17 18 package org.apache.commons.math.linear; 19 20 import org.apache.commons.math.MathRuntimeException; 21 import org.apache.commons.math.exception.util.LocalizedFormats; 22 import org.apache.commons.math.util.FastMath; 23 24 /** 25 * Calculates the LUP-decomposition of a square matrix. 26 * <p>The LUP-decomposition of a matrix A consists of three matrices 27 * L, U and P that satisfy: PA = LU, L is lower triangular, and U is 28 * upper triangular and P is a permutation matrix. All matrices are 29 * m×m.</p> 30 * <p>As shown by the presence of the P matrix, this decomposition is 31 * implemented using partial pivoting.</p> 32 * 33 * @version $Revision: 990655 $ $Date: 2010-08-29 23:49:40 +0200 (dim. 29 aot 2010) $ 34 * @since 2.0 35 */ 36 public class LUDecompositionImpl implements LUDecomposition { 37 38 /** Default bound to determine effective singularity in LU decomposition */ 39 private static final double DEFAULT_TOO_SMALL = 10E-12; 40 41 /** Entries of LU decomposition. */ 42 private double lu[][]; 43 44 /** Pivot permutation associated with LU decomposition */ 45 private int[] pivot; 46 47 /** Parity of the permutation associated with the LU decomposition */ 48 private boolean even; 49 50 /** Singularity indicator. */ 51 private boolean singular; 52 53 /** Cached value of L. */ 54 private RealMatrix cachedL; 55 56 /** Cached value of U. */ 57 private RealMatrix cachedU; 58 59 /** Cached value of P. */ 60 private RealMatrix cachedP; 61 62 /** 63 * Calculates the LU-decomposition of the given matrix. 64 * @param matrix The matrix to decompose. 65 * @exception InvalidMatrixException if matrix is not square 66 */ 67 public LUDecompositionImpl(RealMatrix matrix) 68 throws InvalidMatrixException { 69 this(matrix, DEFAULT_TOO_SMALL); 70 } 71 72 /** 73 * Calculates the LU-decomposition of the given matrix. 74 * @param matrix The matrix to decompose. 75 * @param singularityThreshold threshold (based on partial row norm) 76 * under which a matrix is considered singular 77 * @exception NonSquareMatrixException if matrix is not square 78 */ 79 public LUDecompositionImpl(RealMatrix matrix, double singularityThreshold) 80 throws NonSquareMatrixException { 81 82 if (!matrix.isSquare()) { 83 throw new NonSquareMatrixException(matrix.getRowDimension(), matrix.getColumnDimension()); 84 } 85 86 final int m = matrix.getColumnDimension(); 87 lu = matrix.getData(); 88 pivot = new int[m]; 89 cachedL = null; 90 cachedU = null; 91 cachedP = null; 92 93 // Initialize permutation array and parity 94 for (int row = 0; row < m; row++) { 95 pivot[row] = row; 96 } 97 even = true; 98 singular = false; 99 100 // Loop over columns 101 for (int col = 0; col < m; col++) { 102 103 double sum = 0; 104 105 // upper 106 for (int row = 0; row < col; row++) { 107 final double[] luRow = lu[row]; 108 sum = luRow[col]; 109 for (int i = 0; i < row; i++) { 110 sum -= luRow[i] * lu[i][col]; 111 } 112 luRow[col] = sum; 113 } 114 115 // lower 116 int max = col; // permutation row 117 double largest = Double.NEGATIVE_INFINITY; 118 for (int row = col; row < m; row++) { 119 final double[] luRow = lu[row]; 120 sum = luRow[col]; 121 for (int i = 0; i < col; i++) { 122 sum -= luRow[i] * lu[i][col]; 123 } 124 luRow[col] = sum; 125 126 // maintain best permutation choice 127 if (FastMath.abs(sum) > largest) { 128 largest = FastMath.abs(sum); 129 max = row; 130 } 131 } 132 133 // Singularity check 134 if (FastMath.abs(lu[max][col]) < singularityThreshold) { 135 singular = true; 136 return; 137 } 138 139 // Pivot if necessary 140 if (max != col) { 141 double tmp = 0; 142 final double[] luMax = lu[max]; 143 final double[] luCol = lu[col]; 144 for (int i = 0; i < m; i++) { 145 tmp = luMax[i]; 146 luMax[i] = luCol[i]; 147 luCol[i] = tmp; 148 } 149 int temp = pivot[max]; 150 pivot[max] = pivot[col]; 151 pivot[col] = temp; 152 even = !even; 153 } 154 155 // Divide the lower elements by the "winning" diagonal elt. 156 final double luDiag = lu[col][col]; 157 for (int row = col + 1; row < m; row++) { 158 lu[row][col] /= luDiag; 159 } 160 } 161 162 } 163 164 /** {@inheritDoc} */ 165 public RealMatrix getL() { 166 if ((cachedL == null) && !singular) { 167 final int m = pivot.length; 168 cachedL = MatrixUtils.createRealMatrix(m, m); 169 for (int i = 0; i < m; ++i) { 170 final double[] luI = lu[i]; 171 for (int j = 0; j < i; ++j) { 172 cachedL.setEntry(i, j, luI[j]); 173 } 174 cachedL.setEntry(i, i, 1.0); 175 } 176 } 177 return cachedL; 178 } 179 180 /** {@inheritDoc} */ 181 public RealMatrix getU() { 182 if ((cachedU == null) && !singular) { 183 final int m = pivot.length; 184 cachedU = MatrixUtils.createRealMatrix(m, m); 185 for (int i = 0; i < m; ++i) { 186 final double[] luI = lu[i]; 187 for (int j = i; j < m; ++j) { 188 cachedU.setEntry(i, j, luI[j]); 189 } 190 } 191 } 192 return cachedU; 193 } 194 195 /** {@inheritDoc} */ 196 public RealMatrix getP() { 197 if ((cachedP == null) && !singular) { 198 final int m = pivot.length; 199 cachedP = MatrixUtils.createRealMatrix(m, m); 200 for (int i = 0; i < m; ++i) { 201 cachedP.setEntry(i, pivot[i], 1.0); 202 } 203 } 204 return cachedP; 205 } 206 207 /** {@inheritDoc} */ 208 public int[] getPivot() { 209 return pivot.clone(); 210 } 211 212 /** {@inheritDoc} */ 213 public double getDeterminant() { 214 if (singular) { 215 return 0; 216 } else { 217 final int m = pivot.length; 218 double determinant = even ? 1 : -1; 219 for (int i = 0; i < m; i++) { 220 determinant *= lu[i][i]; 221 } 222 return determinant; 223 } 224 } 225 226 /** {@inheritDoc} */ 227 public DecompositionSolver getSolver() { 228 return new Solver(lu, pivot, singular); 229 } 230 231 /** Specialized solver. */ 232 private static class Solver implements DecompositionSolver { 233 234 /** Entries of LU decomposition. */ 235 private final double lu[][]; 236 237 /** Pivot permutation associated with LU decomposition. */ 238 private final int[] pivot; 239 240 /** Singularity indicator. */ 241 private final boolean singular; 242 243 /** 244 * Build a solver from decomposed matrix. 245 * @param lu entries of LU decomposition 246 * @param pivot pivot permutation associated with LU decomposition 247 * @param singular singularity indicator 248 */ 249 private Solver(final double[][] lu, final int[] pivot, final boolean singular) { 250 this.lu = lu; 251 this.pivot = pivot; 252 this.singular = singular; 253 } 254 255 /** {@inheritDoc} */ 256 public boolean isNonSingular() { 257 return !singular; 258 } 259 260 /** {@inheritDoc} */ 261 public double[] solve(double[] b) 262 throws IllegalArgumentException, InvalidMatrixException { 263 264 final int m = pivot.length; 265 if (b.length != m) { 266 throw MathRuntimeException.createIllegalArgumentException( 267 LocalizedFormats.VECTOR_LENGTH_MISMATCH, b.length, m); 268 } 269 if (singular) { 270 throw new SingularMatrixException(); 271 } 272 273 final double[] bp = new double[m]; 274 275 // Apply permutations to b 276 for (int row = 0; row < m; row++) { 277 bp[row] = b[pivot[row]]; 278 } 279 280 // Solve LY = b 281 for (int col = 0; col < m; col++) { 282 final double bpCol = bp[col]; 283 for (int i = col + 1; i < m; i++) { 284 bp[i] -= bpCol * lu[i][col]; 285 } 286 } 287 288 // Solve UX = Y 289 for (int col = m - 1; col >= 0; col--) { 290 bp[col] /= lu[col][col]; 291 final double bpCol = bp[col]; 292 for (int i = 0; i < col; i++) { 293 bp[i] -= bpCol * lu[i][col]; 294 } 295 } 296 297 return bp; 298 299 } 300 301 /** {@inheritDoc} */ 302 public RealVector solve(RealVector b) 303 throws IllegalArgumentException, InvalidMatrixException { 304 try { 305 return solve((ArrayRealVector) b); 306 } catch (ClassCastException cce) { 307 308 final int m = pivot.length; 309 if (b.getDimension() != m) { 310 throw MathRuntimeException.createIllegalArgumentException( 311 LocalizedFormats.VECTOR_LENGTH_MISMATCH, b.getDimension(), m); 312 } 313 if (singular) { 314 throw new SingularMatrixException(); 315 } 316 317 final double[] bp = new double[m]; 318 319 // Apply permutations to b 320 for (int row = 0; row < m; row++) { 321 bp[row] = b.getEntry(pivot[row]); 322 } 323 324 // Solve LY = b 325 for (int col = 0; col < m; col++) { 326 final double bpCol = bp[col]; 327 for (int i = col + 1; i < m; i++) { 328 bp[i] -= bpCol * lu[i][col]; 329 } 330 } 331 332 // Solve UX = Y 333 for (int col = m - 1; col >= 0; col--) { 334 bp[col] /= lu[col][col]; 335 final double bpCol = bp[col]; 336 for (int i = 0; i < col; i++) { 337 bp[i] -= bpCol * lu[i][col]; 338 } 339 } 340 341 return new ArrayRealVector(bp, false); 342 343 } 344 } 345 346 /** Solve the linear equation A × X = B. 347 * <p>The A matrix is implicit here. It is </p> 348 * @param b right-hand side of the equation A × X = B 349 * @return a vector X such that A × X = B 350 * @exception IllegalArgumentException if matrices dimensions don't match 351 * @exception InvalidMatrixException if decomposed matrix is singular 352 */ 353 public ArrayRealVector solve(ArrayRealVector b) 354 throws IllegalArgumentException, InvalidMatrixException { 355 return new ArrayRealVector(solve(b.getDataRef()), false); 356 } 357 358 /** {@inheritDoc} */ 359 public RealMatrix solve(RealMatrix b) 360 throws IllegalArgumentException, InvalidMatrixException { 361 362 final int m = pivot.length; 363 if (b.getRowDimension() != m) { 364 throw MathRuntimeException.createIllegalArgumentException( 365 LocalizedFormats.DIMENSIONS_MISMATCH_2x2, 366 b.getRowDimension(), b.getColumnDimension(), m, "n"); 367 } 368 if (singular) { 369 throw new SingularMatrixException(); 370 } 371 372 final int nColB = b.getColumnDimension(); 373 374 // Apply permutations to b 375 final double[][] bp = new double[m][nColB]; 376 for (int row = 0; row < m; row++) { 377 final double[] bpRow = bp[row]; 378 final int pRow = pivot[row]; 379 for (int col = 0; col < nColB; col++) { 380 bpRow[col] = b.getEntry(pRow, col); 381 } 382 } 383 384 // Solve LY = b 385 for (int col = 0; col < m; col++) { 386 final double[] bpCol = bp[col]; 387 for (int i = col + 1; i < m; i++) { 388 final double[] bpI = bp[i]; 389 final double luICol = lu[i][col]; 390 for (int j = 0; j < nColB; j++) { 391 bpI[j] -= bpCol[j] * luICol; 392 } 393 } 394 } 395 396 // Solve UX = Y 397 for (int col = m - 1; col >= 0; col--) { 398 final double[] bpCol = bp[col]; 399 final double luDiag = lu[col][col]; 400 for (int j = 0; j < nColB; j++) { 401 bpCol[j] /= luDiag; 402 } 403 for (int i = 0; i < col; i++) { 404 final double[] bpI = bp[i]; 405 final double luICol = lu[i][col]; 406 for (int j = 0; j < nColB; j++) { 407 bpI[j] -= bpCol[j] * luICol; 408 } 409 } 410 } 411 412 return new Array2DRowRealMatrix(bp, false); 413 414 } 415 416 /** {@inheritDoc} */ 417 public RealMatrix getInverse() throws InvalidMatrixException { 418 return solve(MatrixUtils.createRealIdentityMatrix(pivot.length)); 419 } 420 421 } 422 423 } 424
