00001 #ifndef JAMA_CHOLESKY_H
00002 #define JAMA_CHOLESKY_H
00003
00004 #include "math.h"
00005
00006
00007
00008 namespace JAMA
00009 {
00010
00011 using namespace TNT;
00012
00046 template <class Real>
00047 class Cholesky
00048 {
00049 Array2D<Real> L_;
00050 int isspd;
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
00111 for (int j = 0; j < n; j++)
00112 {
00113 Real 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
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
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
00199 #if 0
00200
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
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
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
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
00256
00257 #endif
00258