001/* 002 * BioJava development code 003 * 004 * This code may be freely distributed and modified under the 005 * terms of the GNU Lesser General Public Licence. This should 006 * be distributed with the code. If you do not have a copy, 007 * see: 008 * 009 * http://www.gnu.org/copyleft/lesser.html 010 * 011 * Copyright for this code is held jointly by the individual 012 * authors. These should be listed in @author doc comments. 013 * 014 * For more information on the BioJava project and its aims, 015 * or to join the biojava-l mailing list, visit the home page 016 * at: 017 * 018 * http://www.biojava.org/ 019 * 020 */ 021package org.biojava.nbio.structure.jama; 022 023 /** LU Decomposition. 024 <P> 025 For an m-by-n matrix A with m >= n, the LU decomposition is an m-by-n 026 unit lower triangular matrix L, an n-by-n upper triangular matrix U, 027 and a permutation vector piv of length m so that A(piv,:) = L*U. 028 If m < n, then L is m-by-m and U is m-by-n. 029 <P> 030 The LU decompostion with pivoting always exists, even if the matrix is 031 singular, so the constructor will never fail. The primary use of the 032 LU decomposition is in the solution of square systems of simultaneous 033 linear equations. This will fail if isNonsingular() returns false. 034 */ 035 036public class LUDecomposition implements java.io.Serializable { 037 038 static final long serialVersionUID = 9271028462937843l; 039 040/* ------------------------ 041 Class variables 042 * ------------------------ */ 043 044 /** Array for internal storage of decomposition. 045 @serial internal array storage. 046 */ 047 private double[][] LU; 048 049 /** Row and column dimensions, and pivot sign. 050 @serial column dimension. 051 @serial row dimension. 052 @serial pivot sign. 053 */ 054 private int m, n, pivsign; 055 056 /** Internal storage of pivot vector. 057 @serial pivot vector. 058 */ 059 private int[] piv; 060 061/* ------------------------ 062 Constructor 063 * ------------------------ */ 064 065 /** LU Decomposition provides a data structure to access L, U and piv. 066 @param A Rectangular matrix 067 */ 068 069 public LUDecomposition (Matrix A) { 070 071 // Use a "left-looking", dot-product, Crout/Doolittle algorithm. 072 073 LU = A.getArrayCopy(); 074 m = A.getRowDimension(); 075 n = A.getColumnDimension(); 076 piv = new int[m]; 077 for (int i = 0; i < m; i++) { 078 piv[i] = i; 079 } 080 pivsign = 1; 081 double[] LUrowi; 082 double[] LUcolj = new double[m]; 083 084 // Outer loop. 085 086 for (int j = 0; j < n; j++) { 087 088 // Make a copy of the j-th column to localize references. 089 090 for (int i = 0; i < m; i++) { 091 LUcolj[i] = LU[i][j]; 092 } 093 094 // Apply previous transformations. 095 096 for (int i = 0; i < m; i++) { 097 LUrowi = LU[i]; 098 099 // Most of the time is spent in the following dot product. 100 101 int kmax = Math.min(i,j); 102 double s = 0.0; 103 for (int k = 0; k < kmax; k++) { 104 s += LUrowi[k]*LUcolj[k]; 105 } 106 107 LUrowi[j] = LUcolj[i] -= s; 108 } 109 110 // Find pivot and exchange if necessary. 111 112 int p = j; 113 for (int i = j+1; i < m; i++) { 114 if (Math.abs(LUcolj[i]) > Math.abs(LUcolj[p])) { 115 p = i; 116 } 117 } 118 if (p != j) { 119 for (int k = 0; k < n; k++) { 120 double t = LU[p][k]; LU[p][k] = LU[j][k]; LU[j][k] = t; 121 } 122 int k = piv[p]; piv[p] = piv[j]; piv[j] = k; 123 pivsign = -pivsign; 124 } 125 126 // Compute multipliers. 127 128 if (j < m && LU[j][j] != 0.0) { 129 for (int i = j+1; i < m; i++) { 130 LU[i][j] /= LU[j][j]; 131 } 132 } 133 } 134 } 135 136/* ------------------------ 137 Temporary, experimental code. 138 ------------------------ *\ 139 140 \** LU Decomposition, computed by Gaussian elimination. 141 <P> 142 This constructor computes L and U with the "daxpy"-based elimination 143 algorithm used in LINPACK and MATLAB. In Java, we suspect the dot-product, 144 Crout algorithm will be faster. We have temporarily included this 145 constructor until timing experiments confirm this suspicion. 146 <P> 147 @param A Rectangular matrix 148 @param linpackflag Use Gaussian elimination. Actual value ignored. 149 @return Structure to access L, U and piv. 150 *\ 151 152 public LUDecomposition (Matrix A, int linpackflag) { 153 // Initialize. 154 LU = A.getArrayCopy(); 155 m = A.getRowDimension(); 156 n = A.getColumnDimension(); 157 piv = new int[m]; 158 for (int i = 0; i < m; i++) { 159 piv[i] = i; 160 } 161 pivsign = 1; 162 // Main loop. 163 for (int k = 0; k < n; k++) { 164 // Find pivot. 165 int p = k; 166 for (int i = k+1; i < m; i++) { 167 if (Math.abs(LU[i][k]) > Math.abs(LU[p][k])) { 168 p = i; 169 } 170 } 171 // Exchange if necessary. 172 if (p != k) { 173 for (int j = 0; j < n; j++) { 174 double t = LU[p][j]; LU[p][j] = LU[k][j]; LU[k][j] = t; 175 } 176 int t = piv[p]; piv[p] = piv[k]; piv[k] = t; 177 pivsign = -pivsign; 178 } 179 // Compute multipliers and eliminate k-th column. 180 if (LU[k][k] != 0.0) { 181 for (int i = k+1; i < m; i++) { 182 LU[i][k] /= LU[k][k]; 183 for (int j = k+1; j < n; j++) { 184 LU[i][j] -= LU[i][k]*LU[k][j]; 185 } 186 } 187 } 188 } 189 } 190 191\* ------------------------ 192 End of temporary code. 193 * ------------------------ */ 194 195/* ------------------------ 196 Public Methods 197 * ------------------------ */ 198 199 /** Is the matrix nonsingular? 200 @return true if U, and hence A, is nonsingular. 201 */ 202 203 public boolean isNonsingular () { 204 for (int j = 0; j < n; j++) { 205 if (LU[j][j] == 0) 206 return false; 207 } 208 return true; 209 } 210 211 /** Return lower triangular factor 212 @return L 213 */ 214 215 public Matrix getL () { 216 Matrix X = new Matrix(m,n); 217 double[][] L = X.getArray(); 218 for (int i = 0; i < m; i++) { 219 for (int j = 0; j < n; j++) { 220 if (i > j) { 221 L[i][j] = LU[i][j]; 222 } else if (i == j) { 223 L[i][j] = 1.0; 224 } else { 225 L[i][j] = 0.0; 226 } 227 } 228 } 229 return X; 230 } 231 232 /** Return upper triangular factor 233 @return U 234 */ 235 236 public Matrix getU () { 237 Matrix X = new Matrix(n,n); 238 double[][] U = X.getArray(); 239 for (int i = 0; i < n; i++) { 240 for (int j = 0; j < n; j++) { 241 if (i <= j) { 242 U[i][j] = LU[i][j]; 243 } else { 244 U[i][j] = 0.0; 245 } 246 } 247 } 248 return X; 249 } 250 251 /** Return pivot permutation vector 252 @return piv 253 */ 254 255 public int[] getPivot () { 256 int[] p = new int[m]; 257 for (int i = 0; i < m; i++) { 258 p[i] = piv[i]; 259 } 260 return p; 261 } 262 263 /** Return pivot permutation vector as a one-dimensional double array 264 @return (double) piv 265 */ 266 267 public double[] getDoublePivot () { 268 double[] vals = new double[m]; 269 for (int i = 0; i < m; i++) { 270 vals[i] = piv[i]; 271 } 272 return vals; 273 } 274 275 /** Determinant 276 @return det(A) 277 @exception IllegalArgumentException Matrix must be square 278 */ 279 280 public double det () { 281 if (m != n) { 282 throw new IllegalArgumentException("Matrix must be square."); 283 } 284 double d = pivsign; 285 for (int j = 0; j < n; j++) { 286 d *= LU[j][j]; 287 } 288 return d; 289 } 290 291 /** Solve A*X = B 292 @param B A Matrix with as many rows as A and any number of columns. 293 @return X so that L*U*X = B(piv,:) 294 @exception IllegalArgumentException Matrix row dimensions must agree. 295 @exception RuntimeException Matrix is singular. 296 */ 297 298 public Matrix solve (Matrix B) { 299 if (B.getRowDimension() != m) { 300 throw new IllegalArgumentException("Matrix row dimensions must agree."); 301 } 302 if (!this.isNonsingular()) { 303 throw new RuntimeException("Matrix is singular."); 304 } 305 306 // Copy right hand side with pivoting 307 int nx = B.getColumnDimension(); 308 Matrix Xmat = B.getMatrix(piv,0,nx-1); 309 double[][] X = Xmat.getArray(); 310 311 // Solve L*Y = B(piv,:) 312 for (int k = 0; k < n; k++) { 313 for (int i = k+1; i < n; i++) { 314 for (int j = 0; j < nx; j++) { 315 X[i][j] -= X[k][j]*LU[i][k]; 316 } 317 } 318 } 319 // Solve U*X = Y; 320 for (int k = n-1; k >= 0; k--) { 321 for (int j = 0; j < nx; j++) { 322 X[k][j] /= LU[k][k]; 323 } 324 for (int i = 0; i < k; i++) { 325 for (int j = 0; j < nx; j++) { 326 X[i][j] -= X[k][j]*LU[i][k]; 327 } 328 } 329 } 330 return Xmat; 331 } 332}