I've discovered an error in the algorithm for the (bi)conjugate gradient method. This is equation 2.7.32 on p. 84 of the C version, 2nd ed.

Specifically, the expression for p(k+1) should be
p(k+1) = r(k+1) + alpha(k)*p(k)
and similarly in the equation for "p-bar"

That is, p(k+1) should use the ***most recently computed*** value of r, namely r(k+1) and NOT r(k) as written in the book.

I've tested the algorithm on some simple 2x2 matrices. With this correction, my program converges to the solution in 2 steps as expected. Otherwise, the solution is not reached even after 100 iterations through the loop.

Regards,

Mark Widmer

Dear Mark,

The linbcg algorithm is correct in the book. There may be some misunderstanding about what is "k" and what is "k+1". If you look at the code, you will see that the initial r is computed and used to compute p. Then a new r is computed, and so on. In this sense, the current value of r is always being used to compute p.

Note that the algorithm is useless in practice without preconditioning. The equations that are actually implemented are 2.7.43.

Saul Teukolsky

Saul,

Thanks for replying.

I wasn't speaking directly about the linbcg routine, rather about the equations that describe the underlying algorithm: 2.7.32 and for that matter 2.7.42&43 also. I am programming in VB, but own the "C" book, so I have used the description of the algorithm to write my own version of code.

Whatever k and k+1 are, they should be self-consistent within each set of equations, whether it's 2.7.32, or 2.7.32 plus the substitutions given by 2.7.42&43.

I maintain my position, that the expressions for p[k+1] (in 2.7.32 and 2.7.43) should use k+1 as the subscript for r or z on the right-hand-side. r[k+1] and z[k+1] are the "most recently computed" versions of r & z, according to the notation used in these equations.

I guess we do agree that the most recently computed r's (and z's, for that matter) are to be used when computing p.

Thanks again,

Mark

p.s. Just to add to the confusion, I see that the linbcg routine begins the loop by computing "beta", "p", and "p-bar", whereas the equations 2.7.32 & 2.7.43 begin the loop computing "alpha", "r", and "r-bar".

Dear Mark,
I have NR in C++ 2nd ed. Equation (2.7.43) is printed as p_{k+1}=z_{k+1}+\beta_k * p_k. Same as what you've pointed out. I've read the linbcg algorithm and does not find anything wrong with it.
However, I do find the linbcg algorithm converges VERY slowly. Take this code for example, it takes 100+ iterations to converge at the correct solution, whereas MATLAB bicg function only takes two iterations, which is expected.
DP a[4] = {4, 1, 1, 3};
Mat_DP A(a,2,2);
Vec_DP b(2), x(2);
sa_p = new Vec_DP(5);
ija_p = new Vec_INT(5);
b[0] = 1; b[1] = 2;
x[0] = 2; x[1] = 1;
NR::sprsin(A, 1e-16, *sa_p, *ija_p);
NR::linbcg(b, x, 1, tol, 2, iter, err);

Nonetheless, if I change linbcg to conjugate gradient method by changing the lines
for ( bknum=0.0,j=0; j<n; j++ ) bknum+=z[j]*rr[j];
for (akden=0.0,j=0; j<n; j++) akden+=z[j]*pp[j];
to
for ( bknum=0.0,j=0; j<n; j++ ) bknum+=z[j]*r[j];
for (akden=0.0,j=0; j<n; j++) akden+=z[j]*p[j];
the changed linbcg converges with two iterations.

I know you found this convergence problem in linbcg? How did you solve it?

Lianqing