001 /* 002 * Licensed to the Apache Software Foundation (ASF) under one or more 003 * contributor license agreements. See the NOTICE file distributed with 004 * this work for additional information regarding copyright ownership. 005 * The ASF licenses this file to You under the Apache License, Version 2.0 006 * (the "License"); you may not use this file except in compliance with 007 * the License. You may obtain a copy of the License at 008 * 009 * http://www.apache.org/licenses/LICENSE-2.0 010 * 011 * Unless required by applicable law or agreed to in writing, software 012 * distributed under the License is distributed on an "AS IS" BASIS, 013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 014 * See the License for the specific language governing permissions and 015 * limitations under the License. 016 */ 017 018 package org.apache.commons.math.linear; 019 020 import org.apache.commons.math.MathRuntimeException; 021 022 /** 023 * Calculates the LUP-decomposition of a square matrix. 024 * <p>The LUP-decomposition of a matrix A consists of three matrices 025 * L, U and P that satisfy: PA = LU, L is lower triangular, and U is 026 * upper triangular and P is a permutation matrix. All matrices are 027 * m×m.</p> 028 * <p>As shown by the presence of the P matrix, this decomposition is 029 * implemented using partial pivoting.</p> 030 * 031 * @version $Revision: 799857 $ $Date: 2009-08-01 09:07:12 -0400 (Sat, 01 Aug 2009) $ 032 * @since 2.0 033 */ 034 public class LUDecompositionImpl implements LUDecomposition { 035 036 /** Entries of LU decomposition. */ 037 private double lu[][]; 038 039 /** Pivot permutation associated with LU decomposition */ 040 private int[] pivot; 041 042 /** Parity of the permutation associated with the LU decomposition */ 043 private boolean even; 044 045 /** Singularity indicator. */ 046 private boolean singular; 047 048 /** Cached value of L. */ 049 private RealMatrix cachedL; 050 051 /** Cached value of U. */ 052 private RealMatrix cachedU; 053 054 /** Cached value of P. */ 055 private RealMatrix cachedP; 056 057 /** Default bound to determine effective singularity in LU decomposition */ 058 private static final double DEFAULT_TOO_SMALL = 10E-12; 059 060 /** 061 * Calculates the LU-decomposition of the given matrix. 062 * @param matrix The matrix to decompose. 063 * @exception InvalidMatrixException if matrix is not square 064 */ 065 public LUDecompositionImpl(RealMatrix matrix) 066 throws InvalidMatrixException { 067 this(matrix, DEFAULT_TOO_SMALL); 068 } 069 070 /** 071 * Calculates the LU-decomposition of the given matrix. 072 * @param matrix The matrix to decompose. 073 * @param singularityThreshold threshold (based on partial row norm) 074 * under which a matrix is considered singular 075 * @exception NonSquareMatrixException if matrix is not square 076 */ 077 public LUDecompositionImpl(RealMatrix matrix, double singularityThreshold) 078 throws NonSquareMatrixException { 079 080 if (!matrix.isSquare()) { 081 throw new NonSquareMatrixException(matrix.getRowDimension(), matrix.getColumnDimension()); 082 } 083 084 final int m = matrix.getColumnDimension(); 085 lu = matrix.getData(); 086 pivot = new int[m]; 087 cachedL = null; 088 cachedU = null; 089 cachedP = null; 090 091 // Initialize permutation array and parity 092 for (int row = 0; row < m; row++) { 093 pivot[row] = row; 094 } 095 even = true; 096 singular = false; 097 098 // Loop over columns 099 for (int col = 0; col < m; col++) { 100 101 double sum = 0; 102 103 // upper 104 for (int row = 0; row < col; row++) { 105 final double[] luRow = lu[row]; 106 sum = luRow[col]; 107 for (int i = 0; i < row; i++) { 108 sum -= luRow[i] * lu[i][col]; 109 } 110 luRow[col] = sum; 111 } 112 113 // lower 114 int max = col; // permutation row 115 double largest = Double.NEGATIVE_INFINITY; 116 for (int row = col; row < m; row++) { 117 final double[] luRow = lu[row]; 118 sum = luRow[col]; 119 for (int i = 0; i < col; i++) { 120 sum -= luRow[i] * lu[i][col]; 121 } 122 luRow[col] = sum; 123 124 // maintain best permutation choice 125 if (Math.abs(sum) > largest) { 126 largest = Math.abs(sum); 127 max = row; 128 } 129 } 130 131 // Singularity check 132 if (Math.abs(lu[max][col]) < singularityThreshold) { 133 singular = true; 134 return; 135 } 136 137 // Pivot if necessary 138 if (max != col) { 139 double tmp = 0; 140 final double[] luMax = lu[max]; 141 final double[] luCol = lu[col]; 142 for (int i = 0; i < m; i++) { 143 tmp = luMax[i]; 144 luMax[i] = luCol[i]; 145 luCol[i] = tmp; 146 } 147 int temp = pivot[max]; 148 pivot[max] = pivot[col]; 149 pivot[col] = temp; 150 even = !even; 151 } 152 153 // Divide the lower elements by the "winning" diagonal elt. 154 final double luDiag = lu[col][col]; 155 for (int row = col + 1; row < m; row++) { 156 lu[row][col] /= luDiag; 157 } 158 } 159 160 } 161 162 /** {@inheritDoc} */ 163 public RealMatrix getL() { 164 if ((cachedL == null) && !singular) { 165 final int m = pivot.length; 166 cachedL = MatrixUtils.createRealMatrix(m, m); 167 for (int i = 0; i < m; ++i) { 168 final double[] luI = lu[i]; 169 for (int j = 0; j < i; ++j) { 170 cachedL.setEntry(i, j, luI[j]); 171 } 172 cachedL.setEntry(i, i, 1.0); 173 } 174 } 175 return cachedL; 176 } 177 178 /** {@inheritDoc} */ 179 public RealMatrix getU() { 180 if ((cachedU == null) && !singular) { 181 final int m = pivot.length; 182 cachedU = MatrixUtils.createRealMatrix(m, m); 183 for (int i = 0; i < m; ++i) { 184 final double[] luI = lu[i]; 185 for (int j = i; j < m; ++j) { 186 cachedU.setEntry(i, j, luI[j]); 187 } 188 } 189 } 190 return cachedU; 191 } 192 193 /** {@inheritDoc} */ 194 public RealMatrix getP() { 195 if ((cachedP == null) && !singular) { 196 final int m = pivot.length; 197 cachedP = MatrixUtils.createRealMatrix(m, m); 198 for (int i = 0; i < m; ++i) { 199 cachedP.setEntry(i, pivot[i], 1.0); 200 } 201 } 202 return cachedP; 203 } 204 205 /** {@inheritDoc} */ 206 public int[] getPivot() { 207 return pivot.clone(); 208 } 209 210 /** {@inheritDoc} */ 211 public double getDeterminant() { 212 if (singular) { 213 return 0; 214 } else { 215 final int m = pivot.length; 216 double determinant = even ? 1 : -1; 217 for (int i = 0; i < m; i++) { 218 determinant *= lu[i][i]; 219 } 220 return determinant; 221 } 222 } 223 224 /** {@inheritDoc} */ 225 public DecompositionSolver getSolver() { 226 return new Solver(lu, pivot, singular); 227 } 228 229 /** Specialized solver. */ 230 private static class Solver implements DecompositionSolver { 231 232 /** Entries of LU decomposition. */ 233 private final double lu[][]; 234 235 /** Pivot permutation associated with LU decomposition. */ 236 private final int[] pivot; 237 238 /** Singularity indicator. */ 239 private final boolean singular; 240 241 /** 242 * Build a solver from decomposed matrix. 243 * @param lu entries of LU decomposition 244 * @param pivot pivot permutation associated with LU decomposition 245 * @param singular singularity indicator 246 */ 247 private Solver(final double[][] lu, final int[] pivot, final boolean singular) { 248 this.lu = lu; 249 this.pivot = pivot; 250 this.singular = singular; 251 } 252 253 /** {@inheritDoc} */ 254 public boolean isNonSingular() { 255 return !singular; 256 } 257 258 /** {@inheritDoc} */ 259 public double[] solve(double[] b) 260 throws IllegalArgumentException, InvalidMatrixException { 261 262 final int m = pivot.length; 263 if (b.length != m) { 264 throw MathRuntimeException.createIllegalArgumentException( 265 "vector length mismatch: got {0} but expected {1}", 266 b.length, m); 267 } 268 if (singular) { 269 throw new SingularMatrixException(); 270 } 271 272 final double[] bp = new double[m]; 273 274 // Apply permutations to b 275 for (int row = 0; row < m; row++) { 276 bp[row] = b[pivot[row]]; 277 } 278 279 // Solve LY = b 280 for (int col = 0; col < m; col++) { 281 final double bpCol = bp[col]; 282 for (int i = col + 1; i < m; i++) { 283 bp[i] -= bpCol * lu[i][col]; 284 } 285 } 286 287 // Solve UX = Y 288 for (int col = m - 1; col >= 0; col--) { 289 bp[col] /= lu[col][col]; 290 final double bpCol = bp[col]; 291 for (int i = 0; i < col; i++) { 292 bp[i] -= bpCol * lu[i][col]; 293 } 294 } 295 296 return bp; 297 298 } 299 300 /** {@inheritDoc} */ 301 public RealVector solve(RealVector b) 302 throws IllegalArgumentException, InvalidMatrixException { 303 try { 304 return solve((ArrayRealVector) b); 305 } catch (ClassCastException cce) { 306 307 final int m = pivot.length; 308 if (b.getDimension() != m) { 309 throw MathRuntimeException.createIllegalArgumentException( 310 "vector length mismatch: got {0} but expected {1}", 311 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 "dimensions mismatch: got {0}x{1} but expected {2}x{3}", 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 }