00001 #ifndef JAMA_CHOLESKY_H
00002 #define JAMA_CHOLESKY_H
00003
00004 #include "math.h"
00005
00006
00007
00008 namespace JAMA
00009 {
00010
00043
00044 template <class Real>
00045 class Cholesky
00046 {
00047 TNT::Array2D<Real> L;
00048 int isspd;
00049
00050
00051 public:
00052
00053 Cholesky();
00054 Cholesky(const TNT::Array2D<Real> &A);
00055 TNT::Array2D<Real> getL() const;
00056 int is_spd() const;
00057
00058 };
00059
00060 template <class Real>
00061 Cholesky<Real>::Cholesky() : L(0,0), isspd(0) {}
00062
00067 template <class Real>
00068 int Cholesky<Real>::is_spd() const
00069 {
00070 return isspd;
00071 }
00072
00076 template <class Real>
00077 TNT::Array2D<Real> Cholesky<Real>::getL() const
00078 {
00079 return L;
00080 }
00081
00088 template <class Real>
00089 Cholesky<Real>::Cholesky(const TNT::Array2D<Real> &A)
00090 {
00091
00092
00093 int m = A.dim1();
00094 int n = A.dim2();
00095
00096 isspd = (m == n);
00097
00098 if (m != n)
00099 {
00100 L = TNT::Array2D<Real>(0,0);
00101 return;
00102 }
00103
00104 L = TNT::Array2D<Real>(n,n);
00105
00106
00107
00108 for (int j = 0; j < n; j++)
00109 {
00110 double d = 0.0;
00111 for (int k = 0; k < j; k++)
00112 {
00113 Real s = 0.0;
00114 for (int i = 0; i < k; i++)
00115 {
00116 s += L[k][i]*L[j][i];
00117 }
00118 L[j][k] = s = (A[j][k] - s)/L[k][k];
00119 d = d + s*s;
00120 isspd = isspd && (A[k][j] == A[j][k]);
00121 }
00122 d = A[j][j] - d;
00123 isspd = isspd && (d > 0.0);
00124 L[j][j] = sqrt(d > 0.0 ? d : 0.0);
00125 for (int k = j+1; k < n; k++)
00126 {
00127 L[j][k] = 0.0;
00128 }
00129 }
00130 }
00131
00132 }
00133
00134
00135 #endif
00136