Essentia  2.1-beta6-dev
jama_lu.h
Go to the documentation of this file.
1 #ifndef JAMA_LU_H
2 #define JAMA_LU_H
3
4 #include "tnt.h"
5 #include <algorithm>
6 //for min(), max() below
7
8 using namespace TNT;
9 using namespace std;
10
11 namespace JAMA
12 {
13
26 template <class Real>
27 class LU
28 {
29
30
31
32  /* Array for internal storage of decomposition. */
34  int m, n, pivsign;
36
37
39  const Array1D<int> &piv, int j0, int j1)
40  {
41  int piv_length = piv.dim();
42
43  Array2D<Real> X(piv_length, j1-j0+1);
44
45
46  for (int i = 0; i < piv_length; i++)
47  for (int j = j0; j <= j1; j++)
48  X[i][j-j0] = A[piv[i]][j];
49
50  return X;
51  }
52
54  const Array1D<int> &piv)
55  {
56  int piv_length = piv.dim();
57  if (piv_length != A.dim())
58  return Array1D<Real>();
59
60  Array1D<Real> x(piv_length);
61
62
63  for (int i = 0; i < piv_length; i++)
64  x[i] = A[piv[i]];
65
66  return x;
67  }
68
69
70  public :
71
77  LU (const Array2D<Real> &A) : LU_(A.copy()), m(A.dim1()), n(A.dim2()),
78  piv(A.dim1())
79
80  {
81
82  // Use a "left-looking", dot-product, Crout/Doolittle algorithm.
83
84
85  for (int i = 0; i < m; i++) {
86  piv[i] = i;
87  }
88  pivsign = 1;
89  Real *LUrowi = 0;;
90  Array1D<Real> LUcolj(m);
91
92  // Outer loop.
93
94  for (int j = 0; j < n; j++) {
95
96  // Make a copy of the j-th column to localize references.
97
98  for (int i = 0; i < m; i++) {
99  LUcolj[i] = LU_[i][j];
100  }
101
102  // Apply previous transformations.
103
104  for (int i = 0; i < m; i++) {
105  LUrowi = LU_[i];
106
107  // Most of the time is spent in the following dot product.
108
109  int kmax = min(i,j);
110  double s = 0.0;
111  for (int k = 0; k < kmax; k++) {
112  s += LUrowi[k]*LUcolj[k];
113  }
114
115  LUrowi[j] = LUcolj[i] -= s;
116  }
117
118  // Find pivot and exchange if necessary.
119
120  int p = j;
121  for (int i = j+1; i < m; i++) {
122  if (abs(LUcolj[i]) > abs(LUcolj[p])) {
123  p = i;
124  }
125  }
126  if (p != j) {
127  int k=0;
128  for (k = 0; k < n; k++) {
129  double t = LU_[p][k];
130  LU_[p][k] = LU_[j][k];
131  LU_[j][k] = t;
132  }
133  k = piv[p];
134  piv[p] = piv[j];
135  piv[j] = k;
136  pivsign = -pivsign;
137  }
138
139  // Compute multipliers.
140
141  if ((j < m) && (LU_[j][j] != 0.0)) {
142  for (int i = j+1; i < m; i++) {
143  LU_[i][j] /= LU_[j][j];
144  }
145  }
146  }
147  }
148
149
155  int isNonsingular () {
156  for (int j = 0; j < n; j++) {
157  if (LU_[j][j] == 0)
158  return 0;
159  }
160  return 1;
161  }
162
168  Array2D<Real> L_(m,n);
169  for (int i = 0; i < m; i++) {
170  for (int j = 0; j < n; j++) {
171  if (i > j) {
172  L_[i][j] = LU_[i][j];
173  } else if (i == j) {
174  L_[i][j] = 1.0;
175  } else {
176  L_[i][j] = 0.0;
177  }
178  }
179  }
180  return L_;
181  }
182
188  Array2D<Real> U_(n,n);
189  for (int i = 0; i < n; i++) {
190  for (int j = 0; j < n; j++) {
191  if (i <= j) {
192  U_[i][j] = LU_[i][j];
193  } else {
194  U_[i][j] = 0.0;
195  }
196  }
197  }
198  return U_;
199  }
200
206  return piv;
207  }
208
209
214  Real det () {
215  if (m != n) {
216  return Real(0);
217  }
218  Real d = Real(pivsign);
219  for (int j = 0; j < n; j++) {
220  d *= LU_[j][j];
221  }
222  return d;
223  }
224
232  {
233
234  /* Dimensions: A is mxn, X is nxk, B is mxk */
235
236  if (B.dim1() != m) {
237  return Array2D<Real>(0,0);
238  }
239  if (!isNonsingular()) {
240  return Array2D<Real>(0,0);
241  }
242
243  // Copy right hand side with pivoting
244  int nx = B.dim2();
245
246
247  Array2D<Real> X = permute_copy(B, piv, 0, nx-1);
248
249  // Solve L*Y = B(piv,:)
250  for (int k = 0; k < n; k++) {
251  for (int i = k+1; i < n; i++) {
252  for (int j = 0; j < nx; j++) {
253  X[i][j] -= X[k][j]*LU_[i][k];
254  }
255  }
256  }
257  // Solve U*X = Y;
258  for (int k = n-1; k >= 0; k--) {
259  for (int j = 0; j < nx; j++) {
260  X[k][j] /= LU_[k][k];
261  }
262  for (int i = 0; i < k; i++) {
263  for (int j = 0; j < nx; j++) {
264  X[i][j] -= X[k][j]*LU_[i][k];
265  }
266  }
267  }
268  return X;
269  }
270
271
282  {
283
284  /* Dimensions: A is mxn, X is nxk, B is mxk */
285
286  if (b.dim1() != m) {
287  return Array1D<Real>();
288  }
289  if (!isNonsingular()) {
290  return Array1D<Real>();
291  }
292
293
294  Array1D<Real> x = permute_copy(b, piv);
295
296  // Solve L*Y = B(piv)
297  for (int k = 0; k < n; k++) {
298  for (int i = k+1; i < n; i++) {
299  x[i] -= x[k]*LU_[i][k];
300  }
301  }
302
303  // Solve U*X = Y;
304  for (int k = n-1; k >= 0; k--) {
305  x[k] /= LU_[k][k];
306  for (int i = 0; i < k; i++)
307  x[i] -= x[k]*LU_[i][k];
308  }
309
310
311  return x;
312  }
313
314 }; /* class LU */
315
316 } /* namespace JAMA */
317
318 #endif
319 /* JAMA_LU_H */
Array1D< int > getPivot()
Definition: jama_lu.h:205
int dim1() const
Definition: tnt_array2d.h:231
Real det()
Definition: jama_lu.h:214
LU(const Array2D< Real > &A)
Definition: jama_lu.h:77
Array1D< Real > solve(const Array1D< Real > &b)
Definition: jama_lu.h:281
Definition: tnt_array1d.h:35
Definition: jama_cholesky.h:8
int pivsign
Definition: jama_lu.h:34
Array2D< Real > solve(const Array2D< Real > &B)
Definition: jama_lu.h:231
Definition: jama_lu.h:27
Array2D< Real > getL()
Definition: jama_lu.h:167
int dim2() const
Definition: tnt_array2d.h:234
int dim1() const
Definition: tnt_array1d.h:218
float Real
Definition: types.h:69
Array2D< Real > LU_
Definition: jama_lu.h:33
int dim() const
Definition: tnt_array1d.h:221
Array2D< Real > getU()
Definition: jama_lu.h:187
Array2D< Real > permute_copy(const Array2D< Real > &A, const Array1D< int > &piv, int j0, int j1)
Definition: jama_lu.h:38
Array1D< Real > permute_copy(const Array1D< Real > &A, const Array1D< int > &piv)
Definition: jama_lu.h:53
int isNonsingular()
Definition: jama_lu.h:155
Array1D< int > piv
Definition: jama_lu.h:35