Essentia  2.1-beta5-dev
jama_cholesky.h
Go to the documentation of this file.
1 #ifndef JAMA_CHOLESKY_H
2 #define JAMA_CHOLESKY_H
3 
4 #include "math.h"
5  /* needed for sqrt() below. */
6 
7 
8 namespace JAMA
9 {
10 
11 using namespace TNT;
12 
46 template <class Real>
47 class Cholesky
48 {
49  Array2D<Real> L_; // lower triangular factor
50  int isspd; // 1 if matrix to be factored was SPD
51 
52 public:
53 
54  Cholesky();
55  Cholesky(const Array2D<Real> &A);
56  Array2D<Real> getL() const;
57  Array1D<Real> solve(const Array1D<Real> &B);
58  Array2D<Real> solve(const Array2D<Real> &B);
59  int is_spd() const;
60 
61 };
62 
63 template <class Real>
64 Cholesky<Real>::Cholesky() : L_(0,0), isspd(0) {}
65 
70 template <class Real>
72 {
73  return isspd;
74 }
75 
79 template <class Real>
81 {
82  return L_;
83 }
84 
91 template <class Real>
93 {
94 
95 
96  int m = A.dim1();
97  int n = A.dim2();
98 
99  isspd = (m == n);
100 
101  if (m != n)
102  {
103  L_ = Array2D<Real>(0,0);
104  return;
105  }
106 
107  L_ = Array2D<Real>(n,n);
108 
109 
110  // Main loop.
111  for (int j = 0; j < n; j++)
112  {
113  double d = 0.0;
114  for (int k = 0; k < j; k++)
115  {
116  Real s = 0.0;
117  for (int i = 0; i < k; i++)
118  {
119  s += L_[k][i]*L_[j][i];
120  }
121  L_[j][k] = s = (A[j][k] - s)/L_[k][k];
122  d = d + s*s;
123  isspd = isspd && (A[k][j] == A[j][k]);
124  }
125  d = A[j][j] - d;
126  isspd = isspd && (d > 0.0);
127  L_[j][j] = sqrt(d > 0.0 ? d : 0.0);
128  for (int k = j+1; k < n; k++)
129  {
130  L_[j][k] = 0.0;
131  }
132  }
133 }
134 
145 template <class Real>
147 {
148  int n = L_.dim1();
149  if (b.dim1() != n)
150  return Array1D<Real>();
151 
152 
153  Array1D<Real> x = b.copy();
154 
155 
156  // Solve L*y = b;
157  for (int k = 0; k < n; k++)
158  {
159  for (int i = 0; i < k; i++)
160  x[k] -= x[i]*L_[k][i];
161  x[k] /= L_[k][k];
162 
163  }
164 
165  // Solve L'*X = Y;
166  for (int k = n-1; k >= 0; k--)
167  {
168  for (int i = k+1; i < n; i++)
169  x[k] -= x[i]*L_[i][k];
170  x[k] /= L_[k][k];
171  }
172 
173  return x;
174 }
175 
176 
187 template <class Real>
189 {
190  int n = L_.dim1();
191  if (B.dim1() != n)
192  return Array2D<Real>();
193 
194 
195  Array2D<Real> X = B.copy();
196  int nx = B.dim2();
197 
198 // Cleve's original code
199 #if 0
200  // Solve L*Y = B;
201  for (int k = 0; k < n; k++) {
202  for (int i = k+1; i < n; i++) {
203  for (int j = 0; j < nx; j++) {
204  X[i][j] -= X[k][j]*L_[k][i];
205  }
206  }
207  for (int j = 0; j < nx; j++) {
208  X[k][j] /= L_[k][k];
209  }
210  }
211 
212  // Solve L'*X = Y;
213  for (int k = n-1; k >= 0; k--) {
214  for (int j = 0; j < nx; j++) {
215  X[k][j] /= L_[k][k];
216  }
217  for (int i = 0; i < k; i++) {
218  for (int j = 0; j < nx; j++) {
219  X[i][j] -= X[k][j]*L_[k][i];
220  }
221  }
222  }
223 #endif
224 
225 
226  // Solve L*y = b;
227  for (int j=0; j< nx; j++)
228  {
229  for (int k = 0; k < n; k++)
230  {
231  for (int i = 0; i < k; i++)
232  X[k][j] -= X[i][j]*L_[k][i];
233  X[k][j] /= L_[k][k];
234  }
235  }
236 
237  // Solve L'*X = Y;
238  for (int j=0; j<nx; j++)
239  {
240  for (int k = n-1; k >= 0; k--)
241  {
242  for (int i = k+1; i < n; i++)
243  X[k][j] -= X[i][j]*L_[i][k];
244  X[k][j] /= L_[k][k];
245  }
246  }
247 
248 
249 
250  return X;
251 }
252 
253 
254 }
255 // namespace JAMA
256 
257 #endif
258 // JAMA_CHOLESKY_H
int dim1() const
Definition: tnt_array2d.h:231
int is_spd() const
Definition: jama_cholesky.h:71
Array2D copy() const
Definition: tnt_array2d.h:180
Definition: tnt_array1d.h:35
Definition: jama_cholesky.h:47
Array1D< Real > solve(const Array1D< Real > &B)
Definition: jama_cholesky.h:146
Definition: jama_cholesky.h:8
Cholesky()
Definition: jama_cholesky.h:64
int dim2() const
Definition: tnt_array2d.h:234
int dim1() const
Definition: tnt_array1d.h:218
float Real
Definition: types.h:68
Array2D< Real > getL() const
Definition: jama_cholesky.h:80
Array2D< Real > L_
Definition: jama_cholesky.h:49
Array1D copy() const
Definition: tnt_array1d.h:176
int isspd
Definition: jama_cholesky.h:50