1 #ifndef JAMA_LU_H 2 #define JAMA_LU_H 3 4 #include "tnt.h" 5 #include <algorithm> 6 //for min(), max() below 7 8 using namespace TNT; 9 using namespace std; 10 11 namespace JAMA 12 { 13 14 /** LU Decomposition. 15 <P> 16 For an m-by-n matrix A with m >= n, the LU decomposition is an m-by-n 17 unit lower triangular matrix L, an n-by-n upper triangular matrix U, 18 and a permutation vector piv of length m so that A(piv,:) = L*U. 19 If m < n, then L is m-by-m and U is m-by-n. 20 <P> 21 The LU decompostion with pivoting always exists, even if the matrix is 22 singular, so the constructor will never fail. The primary use of the 23 LU decomposition is in the solution of square systems of simultaneous 24 linear equations. This will fail if isNonsingular() returns false. 25 */ 26 template <class Real> 27 class LU 28 { 29 30 31 32 /* Array for internal storage of decomposition. */ 33 Array2D<Real> LU_; 34 int m, n, pivsign; 35 Array1D<int> piv; 36 37 permute_copy(const Array2D<Real> & A,const Array1D<int> & piv,int j0,int j1)38 Array2D<Real> permute_copy(const Array2D<Real> &A, 39 const Array1D<int> &piv, int j0, int j1) 40 { 41 int piv_length = piv.dim(); 42 43 Array2D<Real> X(piv_length, j1-j0+1); 44 45 46 for (int i = 0; i < piv_length; i++) 47 for (int j = j0; j <= j1; j++) 48 X[i][j-j0] = A[piv[i]][j]; 49 50 return X; 51 } 52 permute_copy(const Array1D<Real> & A,const Array1D<int> & piv)53 Array1D<Real> permute_copy(const Array1D<Real> &A, 54 const Array1D<int> &piv) 55 { 56 int piv_length = piv.dim(); 57 if (piv_length != A.dim()) 58 return Array1D<Real>(); 59 60 Array1D<Real> x(piv_length); 61 62 63 for (int i = 0; i < piv_length; i++) 64 x[i] = A[piv[i]]; 65 66 return x; 67 } 68 69 70 public : 71 72 /** LU Decomposition 73 @param A Rectangular matrix 74 @return LU Decomposition object to access L, U and piv. 75 */ 76 LU(const Array2D<Real> & A)77 LU (const Array2D<Real> &A) : LU_(A.copy()), m(A.dim1()), n(A.dim2()), 78 piv(A.dim1()) 79 80 { 81 82 // Use a "left-looking", dot-product, Crout/Doolittle algorithm. 83 84 85 for (int i = 0; i < m; i++) { 86 piv[i] = i; 87 } 88 pivsign = 1; 89 Real *LUrowi = 0;; 90 Array1D<Real> LUcolj(m); 91 92 // Outer loop. 93 94 for (int j = 0; j < n; j++) { 95 96 // Make a copy of the j-th column to localize references. 97 98 for (int i = 0; i < m; i++) { 99 LUcolj[i] = LU_[i][j]; 100 } 101 102 // Apply previous transformations. 103 104 for (int i = 0; i < m; i++) { 105 LUrowi = LU_[i]; 106 107 // Most of the time is spent in the following dot product. 108 109 int kmax = min(i,j); 110 double s = 0.0; 111 for (int k = 0; k < kmax; k++) { 112 s += LUrowi[k]*LUcolj[k]; 113 } 114 115 LUrowi[j] = LUcolj[i] -= s; 116 } 117 118 // Find pivot and exchange if necessary. 119 120 int p = j; 121 for (int i = j+1; i < m; i++) { 122 if (abs(LUcolj[i]) > abs(LUcolj[p])) { 123 p = i; 124 } 125 } 126 if (p != j) { 127 int k=0; 128 for (k = 0; k < n; k++) { 129 double t = LU_[p][k]; 130 LU_[p][k] = LU_[j][k]; 131 LU_[j][k] = t; 132 } 133 k = piv[p]; 134 piv[p] = piv[j]; 135 piv[j] = k; 136 pivsign = -pivsign; 137 } 138 139 // Compute multipliers. 140 141 if ((j < m) && (LU_[j][j] != 0.0)) { 142 for (int i = j+1; i < m; i++) { 143 LU_[i][j] /= LU_[j][j]; 144 } 145 } 146 } 147 } 148 149 150 /** Is the matrix nonsingular? 151 @return 1 (true) if upper triangular factor U (and hence A) 152 is nonsingular, 0 otherwise. 153 */ 154 isNonsingular()155 int isNonsingular () { 156 for (int j = 0; j < n; j++) { 157 if (LU_[j][j] == 0) 158 return 0; 159 } 160 return 1; 161 } 162 163 /** Return lower triangular factor 164 @return L 165 */ 166 getL()167 Array2D<Real> getL () { 168 Array2D<Real> L_(m,n); 169 for (int i = 0; i < m; i++) { 170 for (int j = 0; j < n; j++) { 171 if (i > j) { 172 L_[i][j] = LU_[i][j]; 173 } else if (i == j) { 174 L_[i][j] = 1.0; 175 } else { 176 L_[i][j] = 0.0; 177 } 178 } 179 } 180 return L_; 181 } 182 183 /** Return upper triangular factor 184 @return U portion of LU factorization. 185 */ 186 getU()187 Array2D<Real> getU () { 188 Array2D<Real> U_(n,n); 189 for (int i = 0; i < n; i++) { 190 for (int j = 0; j < n; j++) { 191 if (i <= j) { 192 U_[i][j] = LU_[i][j]; 193 } else { 194 U_[i][j] = 0.0; 195 } 196 } 197 } 198 return U_; 199 } 200 201 /** Return pivot permutation vector 202 @return piv 203 */ 204 getPivot()205 Array1D<int> getPivot () { 206 return piv; 207 } 208 209 210 /** Compute determinant using LU factors. 211 @return determinant of A, or 0 if A is not square. 212 */ 213 det()214 Real det () { 215 if (m != n) { 216 return Real(0); 217 } 218 Real d = Real(pivsign); 219 for (int j = 0; j < n; j++) { 220 d *= LU_[j][j]; 221 } 222 return d; 223 } 224 225 /** Solve A*X = B 226 @param B A Matrix with as many rows as A and any number of columns. 227 @return X so that L*U*X = B(piv,:), if B is nonconformant, returns 228 0x0 (null) array. 229 */ 230 solve(const Array2D<Real> & B)231 Array2D<Real> solve (const Array2D<Real> &B) 232 { 233 234 /* Dimensions: A is mxn, X is nxk, B is mxk */ 235 236 if (B.dim1() != m) { 237 return Array2D<Real>(0,0); 238 } 239 if (!isNonsingular()) { 240 return Array2D<Real>(0,0); 241 } 242 243 // Copy right hand side with pivoting 244 int nx = B.dim2(); 245 246 247 Array2D<Real> X = permute_copy(B, piv, 0, nx-1); 248 249 // Solve L*Y = B(piv,:) 250 for (int k = 0; k < n; k++) { 251 for (int i = k+1; i < n; i++) { 252 for (int j = 0; j < nx; j++) { 253 X[i][j] -= X[k][j]*LU_[i][k]; 254 } 255 } 256 } 257 // Solve U*X = Y; 258 for (int k = n-1; k >= 0; k--) { 259 for (int j = 0; j < nx; j++) { 260 X[k][j] /= LU_[k][k]; 261 } 262 for (int i = 0; i < k; i++) { 263 for (int j = 0; j < nx; j++) { 264 X[i][j] -= X[k][j]*LU_[i][k]; 265 } 266 } 267 } 268 return X; 269 } 270 271 272 /** Solve A*x = b, where x and b are vectors of length equal 273 to the number of rows in A. 274 275 @param b a vector (Array1D> of length equal to the first dimension 276 of A. 277 @return x a vector (Array1D> so that L*U*x = b(piv), if B is nonconformant, 278 returns 0x0 (null) array. 279 */ 280 solve(const Array1D<Real> & b)281 Array1D<Real> solve (const Array1D<Real> &b) 282 { 283 284 /* Dimensions: A is mxn, X is nxk, B is mxk */ 285 286 if (b.dim1() != m) { 287 return Array1D<Real>(); 288 } 289 if (!isNonsingular()) { 290 return Array1D<Real>(); 291 } 292 293 294 Array1D<Real> x = permute_copy(b, piv); 295 296 // Solve L*Y = B(piv) 297 for (int k = 0; k < n; k++) { 298 for (int i = k+1; i < n; i++) { 299 x[i] -= x[k]*LU_[i][k]; 300 } 301 } 302 303 // Solve U*X = Y; 304 for (int k = n-1; k >= 0; k--) { 305 x[k] /= LU_[k][k]; 306 for (int i = 0; i < k; i++) 307 x[i] -= x[k]*LU_[i][k]; 308 } 309 310 311 return x; 312 } 313 314 }; /* class LU */ 315 316 } /* namespace JAMA */ 317 318 #endif 319 /* JAMA_LU_H */ 320