Hi,

great book, but unclear sometimes. For instance, Cholesky decomposition, the solve function, it is not mentioned anywhere that x is first set equal to b. Rather important!

A

...it is not mentioned anywhere that x is first set equal to b...

Well, it would be important if it were true, but I do not find that to be necessary:

// Whatever forms of includes you need here

#include "../code/nr3.h"
#include "../code/cholesky.h"

// Use Cholesky class to solve equations where
// the matrix is positive semidefinite

int main(void)
{
// Generic loop counters
int i, j;

// Size of the system
const int N = 3;

// An array to feed the MatDoub constructor
const Doub a_d[N * N] = {
1.0e+2, 1.5e+1, 1.0e-2,
1.5e+1, 2.3e+0, 1.0e-2,
1.0e-2, 1.0e-2, 1.0e+0
};

// An array to feed the VecDoub constructor
const Doub b_d[N] = {
4.0e-1,
2.0e-2,
9.9e+1
};

// Declare and initialize the array to feed the Cholesky constructor
MatDoub a(N, N, a_d);

// The constructor does the decomposition
Cholesky chol(a);

// Declare and initialize the array for the right-hand side
// to feed to the solve function
VecDoub b(N, b_d);

// Array to hold the solution
VecDoub x(N);

chol.solve(b, x);

cout << "Original matrix and right-hand side:" << endl;

cout << scientific << setprecision(6) << showpos;
for (i = 0; i < N; i++) {
for (j = 0; j < N; j++) {
cout << setw(15) << a_d[i * N + j] ;
}
cout << " |" << setw(15) << b[i] << endl;
}
cout << endl;

cout << "Product of original matrix times solution from solve:" << endl;

for (i = 0; i < N; i++) {
Doub sum = 0;
for (j = 0; j < N; j++) {
sum += a_d[i * N + j] * x[j];
}
cout << setw(16) << sum << endl;
}
cout << endl;

return 0;
}


Output:


Original matrix and right-hand side:
+1.000000e+02 +1.500000e+01 +1.000000e-02 | +4.000000e-01
+1.500000e+01 +2.300000e+00 +1.000000e-02 | +2.000000e-02
+1.000000e-02 +1.000000e-02 +1.000000e+00 | +9.900000e+01

Product of original matrix times solution from solve:
+4.000000e-01
+2.000000e-02
+9.900000e+01


Or maybe I am missing your point???

Regards,

Dave