CS代写 jama_lu.h Source File

jama_lu.h Source File

Main Page Namespace List Compound List File List Compound Members

Go to the documentation of this file.00001 #ifndef JAMA_LU_H
00002 #define JAMA_LU_H
00004 #include “tnt.h”
00006 using namespace TNT;
00009 namespace JAMA
00024 template
00025 class LU
00030 /* Array for internal storage of decomposition. */
00031 Array2D LU_;
00032 int m, n, pivsign;
00033 Array1D piv;
00036 Array2D permute_copy(const Array2D &A,
00037 const Array1D &piv, int j0, int j1)
00038 {
00039 int piv_length = piv.dim();
00041 Array2D X(piv_length, j1-j0+1);
00044 for (int i = 0; i < piv_length; i++) 00045 for (int j = j0; j <= j1; j++) 00046 X[i][j-j0] = A[piv[i]][j]; 00048 return X; 00049 } 00051 Array1D permute_copy(const Array1D &A,
00052 const Array1D &piv)
00053 {
00054 int piv_length = piv.dim();
00055 if (piv_length != A.dim())
00056 return Array1D();
00058 Array1D x(piv_length);
00061 for (int i = 0; i < piv_length; i++) 00062 x[i] = A[piv[i]]; 00064 return x; 00065 } 00068 public : 00075 LU (const Array2D &A) : LU_(A.copy()), m(A.dim1()), n(A.dim2()),
00076 piv(A.dim1())
00078 {
00080 // Use a “left-looking”, dot-product, Crout/Doolittle algorithm.
00082 int i=0;
00083 int j=0;
00084 int k=0;
00086 for (i = 0; i < m; i++) { 00087 piv[i] = i; 00088 } 00089 pivsign = 1; 00090 Real *LUrowi = 0;; 00091 Array1D LUcolj(m);
00093 // Outer loop.
00095 for (j = 0; j < n; j++) { 00097 // Make a copy of the j-th column to localize references. 00099 for (i = 0; i < m; i++) { 00100 LUcolj[i] = LU_[i][j]; 00101 } 00103 // Apply previous transformations. 00105 for (int i = 0; i < m; i++) { 00106 LUrowi = LU_[i]; 00108 // Most of the time is spent in the following dot product. 00110 int kmax = min(i,j); 00111 double s = 0.0; 00112 for (k = 0; k < kmax; k++) { 00113 s += LUrowi[k]*LUcolj[k]; 00114 } 00116 LUrowi[j] = LUcolj[i] -= s; 00117 } 00119 // Find pivot and exchange if necessary. 00121 int p = j; 00122 for (int i = j+1; i < m; i++) { 00123 if (abs(LUcolj[i]) > abs(LUcolj[p])) {
00124 p = i;
00125 }
00126 }
00127 if (p != j) {
00128 for (k = 0; k < n; k++) { 00129 double t = LU_[p][k]; 00130 LU_[p][k] = LU_[j][k]; 00131 LU_[j][k] = t; 00132 } 00133 k = piv[p]; 00134 piv[p] = piv[j]; 00135 piv[j] = k; 00136 pivsign = -pivsign; 00137 } 00139 // Compute multipliers. 00141 if ((j < m) && (LU_[j][j] != 0.0)) { 00142 for (i = j+1; i < m; i++) { 00143 LU_[i][j] /= LU_[j][j]; 00144 } 00145 } 00146 } 00147 } 00155 int isNonsingular () { 00156 for (int j = 0; j < n; j++) { 00157 if (LU_[j][j] == 0) 00158 return 0; 00159 } 00160 return 1; 00161 } 00167 Array2D getL () {
00168 Array2D L_(m,n);
00169 for (int i = 0; i < m; i++) { 00170 for (int j = 0; j < n; j++) { 00171 if (i > j) {
00172 L_[i][j] = LU_[i][j];
00173 } else if (i == j) {
00174 L_[i][j] = 1.0;
00175 } else {
00176 L_[i][j] = 0.0;
00177 }
00178 }
00179 }
00180 return L_;
00181 }
00187 Array2D getU () {
00188 Array2D U_(n,n);
00189 for (int i = 0; i < n; i++) { 00190 for (int j = 0; j < n; j++) { 00191 if (i <= j) { 00192 U_[i][j] = LU_[i][j]; 00193 } else { 00194 U_[i][j] = 0.0; 00195 } 00196 } 00197 } 00198 return U_; 00199 } 00205 Array1D getPivot () {
00206 return p;
00207 }
00214 Real det () {
00215 if (m != n) {
00216 return Real(0);
00217 }
00218 Real d = Real(pivsign);
00219 for (int j = 0; j < n; j++) { 00220 d *= LU_[j][j]; 00221 } 00222 return d; 00223 } 00231 Array2D solve (const Array2D &B)
00232 {
00234 /* Dimensions: A is mxn, X is nxk, B is mxk */
00236 if (B.dim1() != m) {
00237 return Array2D(0,0);
00238 }
00239 if (!isNonsingular()) {
00240 return Array2D(0,0);
00241 }
00243 // Copy right hand side with pivoting
00244 int nx = B.dim2();
00247 Array2D X = permute_copy(B, piv, 0, nx-1);
00249 // Solve L*Y = B(piv,:)
00250 for (int k = 0; k < n; k++) { 00251 for (int i = k+1; i < n; i++) { 00252 for (int j = 0; j < nx; j++) { 00253 X[i][j] -= X[k][j]*LU_[i][k]; 00254 } 00255 } 00256 } 00257 // Solve U*X = Y; 00258 for (int k = n-1; k >= 0; k–) {
00259 for (int j = 0; j < nx; j++) { 00260 X[k][j] /= LU_[k][k]; 00261 } 00262 for (int i = 0; i < k; i++) { 00263 for (int j = 0; j < nx; j++) { 00264 X[i][j] -= X[k][j]*LU_[i][k]; 00265 } 00266 } 00267 } 00268 return X; 00269 } 00281 Array1D solve (const Array1D &b)
00282 {
00284 /* Dimensions: A is mxn, X is nxk, B is mxk */
00286 if (b.dim1() != m) {
00287 return Array1D();
00288 }
00289 if (!isNonsingular()) {
00290 return Array1D();
00291 }
00294 Array1D x = permute_copy(b, piv);
00296 // Solve L*Y = B(piv)
00297 for (int k = 0; k < n; k++) { 00298 for (int i = k+1; i < n; i++) { 00299 x[i] -= x[k]*LU_[i][k]; 00300 } 00301 } 00303 // Solve U*X = Y; 00304 for (int k = n-1; k >= 0; k–) {
00305 x[k] /= LU_[k][k];
00306 for (int i = 0; i < k; i++) 00307 x[i] -= x[k]*LU_[i][k]; 00308 } 00311 return x; 00312 } 00314 }; /* class LU */ 00316 } /* namespace JAMA */ 00318 #endif 00319 /* JAMA_LU_H */ Generated at Mon Jan 20 07:47:17 2003 for JAMA/C++ by 1.2.5 written by Dimitri van Heesch, © 1997-2001 程序代写 CS代考加微信: powcoder QQ: 1823890830 Email: powcoder@163.com

Related Posts