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
00003
00004 #include "math.h"
00005         /* needed for sqrt() below. */
00006
00007
00008 namespace JAMA
00009 {
00010
00011 using namespace TNT;
00012
00045
00046 template <class Real>
00047 class Cholesky
00048 {
00049         Array2D<Real> L_;               // lower triangular factor
00050         int isspd;                              // 1 if matrix to be factored was SPD
00051
00052 public:
00053
00054         Cholesky();
00055         Cholesky(const Array2D<Real> &A);
00056         Array2D<Real> getL() const;
00057         Array1D<Real> solve(const Array1D<Real> &B);
00058         Array2D<Real> solve(const Array2D<Real> &B);
00059         int is_spd() const;
00060
00061 };
00062
00063 template <class Real>
00064 Cholesky<Real>::Cholesky() : L_(0,0), isspd(0) {}
00065
00070 template <class Real>
00071 int Cholesky<Real>::is_spd() const
00072 {
00073         return isspd;
00074 }
00075
00079 template <class Real>
00080 Array2D<Real> Cholesky<Real>::getL() const
00081 {
00082         return L_;
00083 }
00084
00091 template <class Real>
00092 Cholesky<Real>::Cholesky(const Array2D<Real> &A)
00093 {
00094
00095
00096         int m = A.dim1();
00097         int n = A.dim2();
00098
00099         isspd = (m == n);
00100
00101         if (m != n)
00102         {
00103                 L_ = Array2D<Real>(0,0);
00104                 return;
00105         }
00106
00107         L_ = Array2D<Real>(n,n);
00108
00109
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         }
00133 }
00134
00145 template <class Real>
00146 Array1D<Real> Cholesky<Real>::solve(const Array1D<Real> &b)
00147 {
00148         int n = L_.dim1();
00149         if (b.dim1() != n)
00150                 return Array1D<Real>();
00151
00152
00153         Array1D<Real> x = b.copy();
00154
00155
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];
00162
00163       }
00164
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       }
00172
00173         return x;
00174 }
00175
00176
00187 template <class Real>
00188 Array2D<Real> Cholesky<Real>::solve(const Array2D<Real> &B)
00189 {
00190         int n = L_.dim1();
00191         if (B.dim1() != n)
00192                 return Array2D<Real>();
00193
00194
00195         Array2D<Real> X = B.copy();
00196         int nx = B.dim2();
00197
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       }
00211
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
00224
00225
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       }
00236
00237       // Solve L'*X = Y;
00238      for (int j=0; j<nx; j++)
00239          {
00240         for (int k = n-1; k >= 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       }
00247
00248
00249
00250         return X;
00251 }
00252
00253
00254 }
00255 // namespace JAMA
00256
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