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.math3.linear; 019 020 import java.lang.reflect.Array; 021 022 import org.apache.commons.math3.Field; 023 import org.apache.commons.math3.FieldElement; 024 import org.apache.commons.math3.exception.DimensionMismatchException; 025 026 /** 027 * Calculates the LUP-decomposition of a square matrix. 028 * <p>The LUP-decomposition of a matrix A consists of three matrices 029 * L, U and P that satisfy: PA = LU, L is lower triangular, and U is 030 * upper triangular and P is a permutation matrix. All matrices are 031 * m×m.</p> 032 * <p>Since {@link FieldElement field elements} do not provide an ordering 033 * operator, the permutation matrix is computed here only in order to avoid 034 * a zero pivot element, no attempt is done to get the largest pivot 035 * element.</p> 036 * <p>This class is based on the class with similar name from the 037 * <a href="http://math.nist.gov/javanumerics/jama/">JAMA</a> library.</p> 038 * <ul> 039 * <li>a {@link #getP() getP} method has been added,</li> 040 * <li>the {@code det} method has been renamed as {@link #getDeterminant() 041 * getDeterminant},</li> 042 * <li>the {@code getDoublePivot} method has been removed (but the int based 043 * {@link #getPivot() getPivot} method has been kept),</li> 044 * <li>the {@code solve} and {@code isNonSingular} methods have been replaced 045 * by a {@link #getSolver() getSolver} method and the equivalent methods 046 * provided by the returned {@link DecompositionSolver}.</li> 047 * </ul> 048 * 049 * @param <T> the type of the field elements 050 * @see <a href="http://mathworld.wolfram.com/LUDecomposition.html">MathWorld</a> 051 * @see <a href="http://en.wikipedia.org/wiki/LU_decomposition">Wikipedia</a> 052 * @version $Id: FieldLUDecomposition.java 1416643 2012-12-03 19:37:14Z tn $ 053 * @since 2.0 (changed to concrete class in 3.0) 054 */ 055 public class FieldLUDecomposition<T extends FieldElement<T>> { 056 057 /** Field to which the elements belong. */ 058 private final Field<T> field; 059 060 /** Entries of LU decomposition. */ 061 private T[][] lu; 062 063 /** Pivot permutation associated with LU decomposition. */ 064 private int[] pivot; 065 066 /** Parity of the permutation associated with the LU decomposition. */ 067 private boolean even; 068 069 /** Singularity indicator. */ 070 private boolean singular; 071 072 /** Cached value of L. */ 073 private FieldMatrix<T> cachedL; 074 075 /** Cached value of U. */ 076 private FieldMatrix<T> cachedU; 077 078 /** Cached value of P. */ 079 private FieldMatrix<T> cachedP; 080 081 /** 082 * Calculates the LU-decomposition of the given matrix. 083 * @param matrix The matrix to decompose. 084 * @throws NonSquareMatrixException if matrix is not square 085 */ 086 public FieldLUDecomposition(FieldMatrix<T> matrix) { 087 if (!matrix.isSquare()) { 088 throw new NonSquareMatrixException(matrix.getRowDimension(), 089 matrix.getColumnDimension()); 090 } 091 092 final int m = matrix.getColumnDimension(); 093 field = matrix.getField(); 094 lu = matrix.getData(); 095 pivot = new int[m]; 096 cachedL = null; 097 cachedU = null; 098 cachedP = null; 099 100 // Initialize permutation array and parity 101 for (int row = 0; row < m; row++) { 102 pivot[row] = row; 103 } 104 even = true; 105 singular = false; 106 107 // Loop over columns 108 for (int col = 0; col < m; col++) { 109 110 T sum = field.getZero(); 111 112 // upper 113 for (int row = 0; row < col; row++) { 114 final T[] luRow = lu[row]; 115 sum = luRow[col]; 116 for (int i = 0; i < row; i++) { 117 sum = sum.subtract(luRow[i].multiply(lu[i][col])); 118 } 119 luRow[col] = sum; 120 } 121 122 // lower 123 int nonZero = col; // permutation row 124 for (int row = col; row < m; row++) { 125 final T[] luRow = lu[row]; 126 sum = luRow[col]; 127 for (int i = 0; i < col; i++) { 128 sum = sum.subtract(luRow[i].multiply(lu[i][col])); 129 } 130 luRow[col] = sum; 131 132 if (lu[nonZero][col].equals(field.getZero())) { 133 // try to select a better permutation choice 134 ++nonZero; 135 } 136 } 137 138 // Singularity check 139 if (nonZero >= m) { 140 singular = true; 141 return; 142 } 143 144 // Pivot if necessary 145 if (nonZero != col) { 146 T tmp = field.getZero(); 147 for (int i = 0; i < m; i++) { 148 tmp = lu[nonZero][i]; 149 lu[nonZero][i] = lu[col][i]; 150 lu[col][i] = tmp; 151 } 152 int temp = pivot[nonZero]; 153 pivot[nonZero] = pivot[col]; 154 pivot[col] = temp; 155 even = !even; 156 } 157 158 // Divide the lower elements by the "winning" diagonal elt. 159 final T luDiag = lu[col][col]; 160 for (int row = col + 1; row < m; row++) { 161 final T[] luRow = lu[row]; 162 luRow[col] = luRow[col].divide(luDiag); 163 } 164 } 165 166 } 167 168 /** 169 * Returns the matrix L of the decomposition. 170 * <p>L is a lower-triangular matrix</p> 171 * @return the L matrix (or null if decomposed matrix is singular) 172 */ 173 public FieldMatrix<T> getL() { 174 if ((cachedL == null) && !singular) { 175 final int m = pivot.length; 176 cachedL = new Array2DRowFieldMatrix<T>(field, m, m); 177 for (int i = 0; i < m; ++i) { 178 final T[] luI = lu[i]; 179 for (int j = 0; j < i; ++j) { 180 cachedL.setEntry(i, j, luI[j]); 181 } 182 cachedL.setEntry(i, i, field.getOne()); 183 } 184 } 185 return cachedL; 186 } 187 188 /** 189 * Returns the matrix U of the decomposition. 190 * <p>U is an upper-triangular matrix</p> 191 * @return the U matrix (or null if decomposed matrix is singular) 192 */ 193 public FieldMatrix<T> getU() { 194 if ((cachedU == null) && !singular) { 195 final int m = pivot.length; 196 cachedU = new Array2DRowFieldMatrix<T>(field, m, m); 197 for (int i = 0; i < m; ++i) { 198 final T[] luI = lu[i]; 199 for (int j = i; j < m; ++j) { 200 cachedU.setEntry(i, j, luI[j]); 201 } 202 } 203 } 204 return cachedU; 205 } 206 207 /** 208 * Returns the P rows permutation matrix. 209 * <p>P is a sparse matrix with exactly one element set to 1.0 in 210 * each row and each column, all other elements being set to 0.0.</p> 211 * <p>The positions of the 1 elements are given by the {@link #getPivot() 212 * pivot permutation vector}.</p> 213 * @return the P rows permutation matrix (or null if decomposed matrix is singular) 214 * @see #getPivot() 215 */ 216 public FieldMatrix<T> getP() { 217 if ((cachedP == null) && !singular) { 218 final int m = pivot.length; 219 cachedP = new Array2DRowFieldMatrix<T>(field, m, m); 220 for (int i = 0; i < m; ++i) { 221 cachedP.setEntry(i, pivot[i], field.getOne()); 222 } 223 } 224 return cachedP; 225 } 226 227 /** 228 * Returns the pivot permutation vector. 229 * @return the pivot permutation vector 230 * @see #getP() 231 */ 232 public int[] getPivot() { 233 return pivot.clone(); 234 } 235 236 /** 237 * Return the determinant of the matrix. 238 * @return determinant of the matrix 239 */ 240 public T getDeterminant() { 241 if (singular) { 242 return field.getZero(); 243 } else { 244 final int m = pivot.length; 245 T determinant = even ? field.getOne() : field.getZero().subtract(field.getOne()); 246 for (int i = 0; i < m; i++) { 247 determinant = determinant.multiply(lu[i][i]); 248 } 249 return determinant; 250 } 251 } 252 253 /** 254 * Get a solver for finding the A × X = B solution in exact linear sense. 255 * @return a solver 256 */ 257 public FieldDecompositionSolver<T> getSolver() { 258 return new Solver<T>(field, lu, pivot, singular); 259 } 260 261 /** Specialized solver. */ 262 private static class Solver<T extends FieldElement<T>> implements FieldDecompositionSolver<T> { 263 264 /** Field to which the elements belong. */ 265 private final Field<T> field; 266 267 /** Entries of LU decomposition. */ 268 private final T[][] lu; 269 270 /** Pivot permutation associated with LU decomposition. */ 271 private final int[] pivot; 272 273 /** Singularity indicator. */ 274 private final boolean singular; 275 276 /** 277 * Build a solver from decomposed matrix. 278 * @param field field to which the matrix elements belong 279 * @param lu entries of LU decomposition 280 * @param pivot pivot permutation associated with LU decomposition 281 * @param singular singularity indicator 282 */ 283 private Solver(final Field<T> field, final T[][] lu, 284 final int[] pivot, final boolean singular) { 285 this.field = field; 286 this.lu = lu; 287 this.pivot = pivot; 288 this.singular = singular; 289 } 290 291 /** {@inheritDoc} */ 292 public boolean isNonSingular() { 293 return !singular; 294 } 295 296 /** {@inheritDoc} */ 297 public FieldVector<T> solve(FieldVector<T> b) { 298 try { 299 return solve((ArrayFieldVector<T>) b); 300 } catch (ClassCastException cce) { 301 302 final int m = pivot.length; 303 if (b.getDimension() != m) { 304 throw new DimensionMismatchException(b.getDimension(), m); 305 } 306 if (singular) { 307 throw new SingularMatrixException(); 308 } 309 310 @SuppressWarnings("unchecked") // field is of type T 311 final T[] bp = (T[]) Array.newInstance(field.getRuntimeClass(), m); 312 313 // Apply permutations to b 314 for (int row = 0; row < m; row++) { 315 bp[row] = b.getEntry(pivot[row]); 316 } 317 318 // Solve LY = b 319 for (int col = 0; col < m; col++) { 320 final T bpCol = bp[col]; 321 for (int i = col + 1; i < m; i++) { 322 bp[i] = bp[i].subtract(bpCol.multiply(lu[i][col])); 323 } 324 } 325 326 // Solve UX = Y 327 for (int col = m - 1; col >= 0; col--) { 328 bp[col] = bp[col].divide(lu[col][col]); 329 final T bpCol = bp[col]; 330 for (int i = 0; i < col; i++) { 331 bp[i] = bp[i].subtract(bpCol.multiply(lu[i][col])); 332 } 333 } 334 335 return new ArrayFieldVector<T>(field, bp, false); 336 337 } 338 } 339 340 /** Solve the linear equation A × X = B. 341 * <p>The A matrix is implicit here. It is </p> 342 * @param b right-hand side of the equation A × X = B 343 * @return a vector X such that A × X = B 344 * @throws DimensionMismatchException if the matrices dimensions do not match. 345 * @throws SingularMatrixException if the decomposed matrix is singular. 346 */ 347 public ArrayFieldVector<T> solve(ArrayFieldVector<T> b) { 348 final int m = pivot.length; 349 final int length = b.getDimension(); 350 if (length != m) { 351 throw new DimensionMismatchException(length, m); 352 } 353 if (singular) { 354 throw new SingularMatrixException(); 355 } 356 357 @SuppressWarnings("unchecked") 358 // field is of type T 359 final T[] bp = (T[]) Array.newInstance(field.getRuntimeClass(), 360 m); 361 362 // Apply permutations to b 363 for (int row = 0; row < m; row++) { 364 bp[row] = b.getEntry(pivot[row]); 365 } 366 367 // Solve LY = b 368 for (int col = 0; col < m; col++) { 369 final T bpCol = bp[col]; 370 for (int i = col + 1; i < m; i++) { 371 bp[i] = bp[i].subtract(bpCol.multiply(lu[i][col])); 372 } 373 } 374 375 // Solve UX = Y 376 for (int col = m - 1; col >= 0; col--) { 377 bp[col] = bp[col].divide(lu[col][col]); 378 final T bpCol = bp[col]; 379 for (int i = 0; i < col; i++) { 380 bp[i] = bp[i].subtract(bpCol.multiply(lu[i][col])); 381 } 382 } 383 384 return new ArrayFieldVector<T>(bp, false); 385 } 386 387 /** {@inheritDoc} */ 388 public FieldMatrix<T> solve(FieldMatrix<T> b) { 389 final int m = pivot.length; 390 if (b.getRowDimension() != m) { 391 throw new DimensionMismatchException(b.getRowDimension(), m); 392 } 393 if (singular) { 394 throw new SingularMatrixException(); 395 } 396 397 final int nColB = b.getColumnDimension(); 398 399 // Apply permutations to b 400 @SuppressWarnings("unchecked") // field is of type T 401 final T[][] bp = (T[][]) Array.newInstance(field.getRuntimeClass(), new int[] { m, nColB }); 402 for (int row = 0; row < m; row++) { 403 final T[] bpRow = bp[row]; 404 final int pRow = pivot[row]; 405 for (int col = 0; col < nColB; col++) { 406 bpRow[col] = b.getEntry(pRow, col); 407 } 408 } 409 410 // Solve LY = b 411 for (int col = 0; col < m; col++) { 412 final T[] bpCol = bp[col]; 413 for (int i = col + 1; i < m; i++) { 414 final T[] bpI = bp[i]; 415 final T luICol = lu[i][col]; 416 for (int j = 0; j < nColB; j++) { 417 bpI[j] = bpI[j].subtract(bpCol[j].multiply(luICol)); 418 } 419 } 420 } 421 422 // Solve UX = Y 423 for (int col = m - 1; col >= 0; col--) { 424 final T[] bpCol = bp[col]; 425 final T luDiag = lu[col][col]; 426 for (int j = 0; j < nColB; j++) { 427 bpCol[j] = bpCol[j].divide(luDiag); 428 } 429 for (int i = 0; i < col; i++) { 430 final T[] bpI = bp[i]; 431 final T luICol = lu[i][col]; 432 for (int j = 0; j < nColB; j++) { 433 bpI[j] = bpI[j].subtract(bpCol[j].multiply(luICol)); 434 } 435 } 436 } 437 438 return new Array2DRowFieldMatrix<T>(field, bp, false); 439 440 } 441 442 /** {@inheritDoc} */ 443 public FieldMatrix<T> getInverse() { 444 final int m = pivot.length; 445 final T one = field.getOne(); 446 FieldMatrix<T> identity = new Array2DRowFieldMatrix<T>(field, m, m); 447 for (int i = 0; i < m; ++i) { 448 identity.setEntry(i, i, one); 449 } 450 return solve(identity); 451 } 452 } 453 }