RTXI  2.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
cholesky.cpp
Go to the documentation of this file.
1 //
2 // File = cholesky.cpp
3 //
4 
5 #include "cholesky.h"
6 #include "complex.h"
7 #include <fstream>
8 
9 #ifdef _DEBUG
10 extern std::ofstream DebugFile;
11 #endif
12 
13 int
14 CholeskyDecomp(int ord, matrix<complex>* ax, complex* bx, double epsilon)
15 {
16  double* dx;
17  int i, j, k;
18  complex* yx;
19  matrix<complex> lx(1, ord, 1, ord);
20 
21  dx = new double[ord + 1];
22  yx = new complex[ord + 1];
23 
24  dx[1] = real((*ax)[1][1]);
25  for (i = 2; i <= ord; i++) {
26  for (j = 1; j < i; j++) {
27  lx[i][j] = conj((*ax)[j][i]) / dx[j];
28  if (j == 1)
29  continue;
30 
31  for (k = 1; k < j; k++) {
32  lx[i][j] -= lx[i][k] * dx[k] * conj(lx[j][k]) / dx[j];
33  } // end of loop over k
34  } // end of loop over j
35  dx[i] = real((*ax)[i][i]);
36 
37  for (k = 1; k < i; k++) {
38  dx[i] -= dx[k] * mag_sqrd(lx[i][k]);
39  }
40  //
41  // Non-positive dx[i] is an error condition most likely
42  // caused by ill-conditioned ax matrix
43 
44  if (dx[i] < 0.0)
45  // if( dx[i] < epsilon )
46  {
47  delete[] dx;
48  delete[] yx;
49  return (-1);
50  }
51  } // end of loop over i
52 
53  //- - - - - - - - - - - - - - - - - -
54  // Solve for y vector in Eq. (31.x)
55  //
56  yx[1] = bx[1];
57  for (k = 2; k <= ord; k++) {
58  yx[k] = bx[k];
59 
60  for (j = 1; j < k; j++) {
61  yx[k] -= lx[k][j] * yx[j];
62  }
63  } // end of loop over k
64 
65  //- - - - - - - - - - - - - - - - - - -
66  // Solve for x vector in Eq. (31.x)
67  //
68  bx[ord] = yx[ord] / dx[ord];
69 
70  for (k = ord - 1; k >= 1; k--) {
71  bx[k] = yx[k] / dx[k];
72 
73  for (j = k + 1; j <= ord; j++) {
74  bx[k] -= conj(lx[j][k]) * bx[j];
75  }
76  }
77  delete[] dx;
78  delete[] yx;
79  return (0);
80 }
complex conj(const complex _z)
Definition: complex.h:108
int CholeskyDecomp(int ord, matrix< complex > *ax, complex *bx, double epsilon)
Definition: cholesky.cpp:14
std::ofstream DebugFile
double mag_sqrd(const complex _z)
Definition: complex.h:128
double real(const complex &_z)
Definition: complex.h:98