程序代写 jama_cholesky.h Source File

jama_cholesky.h Source File

Main Page Namespace List Compound List File List Compound Members

jama_cholesky.h
Go to the documentation of this file.00001 #ifndef JAMA_CHOLESKY_H
00002 #define JAMA_CHOLESKY_H
00004 #include “math.h”
00005 /* needed for sqrt() below. */
00008 namespace JAMA
00011 using namespace TNT;
00046 template
00047 class Cholesky
00049 Array2D L_; // lower triangular factor
00050 int isspd; // 1 if matrix to be factored was SPD
00052 public:
00054 Cholesky();
00055 Cholesky(const Array2D &A);
00056 Array2D getL() const;
00057 Array1D solve(const Array1D &B);
00058 Array2D solve(const Array2D &B);
00059 int is_spd() const;
00063 template
00064 Cholesky::Cholesky() : L_(0,0), isspd(0) {}
00070 template
00071 int Cholesky::is_spd() const
00073 return isspd;
00079 template
00080 Array2D Cholesky::getL() const
00082 return L_;
00091 template
00092 Cholesky::Cholesky(const Array2D &A)
00096 int m = A.dim1();
00097 int n = A.dim2();
00099 isspd = (m == n);
00101 if (m != n)
00102 {
00103 L_ = Array2D(0,0);
00104 return;
00105 }
00107 L_ = Array2D(n,n);
00110 // Main loop.
00111 for (int j = 0; j < n; j++) 00112 { 00113 double d = 0.0; 00114 for (int k = 0; k < j; k++) 00115 { 00116 Real s = 0.0; 00117 for (int i = 0; i < k; i++) 00118 { 00119 s += L_[k][i]*L_[j][i]; 00120 } 00121 L_[j][k] = s = (A[j][k] - s)/L_[k][k]; 00122 d = d + s*s; 00123 isspd = isspd && (A[k][j] == A[j][k]); 00124 } 00125 d = A[j][j] - d; 00126 isspd = isspd && (d > 0.0);
00127 L_[j][j] = sqrt(d > 0.0 ? d : 0.0);
00128 for (int k = j+1; k < n; k++) 00129 { 00130 L_[j][k] = 0.0; 00131 } 00132 } 00145 template
00146 Array1D Cholesky::solve(const Array1D &b)
00148 int n = L_.dim1();
00149 if (b.dim1() != n)
00150 return Array1D();
00153 Array1D x = b.copy();
00156 // Solve L*y = b;
00157 for (int k = 0; k < n; k++) 00158 { 00159 for (int i = 0; i < k; i++) 00160 x[k] -= x[i]*L_[k][i]; 00161 x[k] /= L_[k][k]; 00163 } 00165 // Solve L'*X = Y; 00166 for (int k = n-1; k >= 0; k–)
00167 {
00168 for (int i = k+1; i < n; i++) 00169 x[k] -= x[i]*L_[i][k]; 00170 x[k] /= L_[k][k]; 00171 } 00173 return x; 00187 template
00188 Array2D Cholesky::solve(const Array2D &B)
00190 int n = L_.dim1();
00191 if (B.dim1() != n)
00192 return Array2D();
00195 Array2D X = B.copy();
00196 int nx = B.dim2();
00198 // Cleve’s original code
00199 #if 0
00200 // Solve L*Y = B;
00201 for (int k = 0; k < n; k++) { 00202 for (int i = k+1; i < n; i++) { 00203 for (int j = 0; j < nx; j++) { 00204 X[i][j] -= X[k][j]*L_[k][i]; 00205 } 00206 } 00207 for (int j = 0; j < nx; j++) { 00208 X[k][j] /= L_[k][k]; 00209 } 00210 } 00212 // Solve L'*X = Y; 00213 for (int k = n-1; k >= 0; k–) {
00214 for (int j = 0; j < nx; j++) { 00215 X[k][j] /= L_[k][k]; 00216 } 00217 for (int i = 0; i < k; i++) { 00218 for (int j = 0; j < nx; j++) { 00219 X[i][j] -= X[k][j]*L_[k][i]; 00220 } 00221 } 00222 } 00223 #endif 00226 // Solve L*y = b; 00227 for (int j=0; j< nx; j++) 00228 { 00229 for (int k = 0; k < n; k++) 00230 { 00231 for (int i = 0; i < k; i++) 00232 X[k][j] -= X[i][j]*L_[k][i]; 00233 X[k][j] /= L_[k][k]; 00234 } 00235 } 00237 // Solve L'*X = Y; 00238 for (int j=0; j= 0; k–)
00241 {
00242 for (int i = k+1; i < n; i++) 00243 X[k][j] -= X[i][j]*L_[i][k]; 00244 X[k][j] /= L_[k][k]; 00245 } 00246 } 00250 return X; 00255 // namespace JAMA 00257 #endif 00258 // JAMA_CHOLESKY_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