nurbz
02-14-2008, 04:19 AM
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
davekw7x
02-14-2008, 11:42 AM
...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