[Quantlib-dev] BiCGstab

[Peter C. <pca...@gm...>]

Hi Ralph, Klaus,

I am using your biconjugate gradient stabilized solver and experience
some issues in certain cases. A small test case would be the system
Ax=b with

A=
| 1 0 0 0 0 0 0 0 0 |
| 0 1 0 0 0 0 0 0 0 |
| 0 0 1 0 0 0 0 0 0 |
| 0 0 0 1 0 0 0 0 0 |
| 0 0 0 0 1 0 0 0 0 |
| 0 0 -0.5 0 0 1 0 0 -0.5 |
| 0 0 0 -0.5 0 0 1 -0.5 0 |
| 0 0 0 0 0 0 0 1 0 |
| 0 0 0 0 0 -0.5 0 -0.5 1 |

b = [ 5; 5; 5; 5; 5; 0; 0; 5; 0 ]

With the guess [5,5,...,5], which is the exact solution, the solver
throws an exception "could not converge". This is because "error" is
initialized to Null<Real>() in L52, and the iteration loop is exited
at L58 already (because rho=0.0), so the Null<Real>() value is checked
in L90. I think this can be fixed by simply changing L52 to

        Real error=norm2(r)/bnorm2;

With the guess [0,0,....,0] on the other hand, the same exception is
thrown. We have

x = [ 5; 5; 5; 5; 5; 2.77778; 5.55556; 5; 2.77778 ]
r = [ 0; 0; 0; 0; 0; 1.11111; -0.555556; 0; 1.11111 ]
rTld = [ 5; 5; 5; 5; 5; 0; 0; 5; 0 ]

at L86 after the first iteration, so rho = 0.0 again, but this time
exiting (at L57) with a huge error of course. Is that expected ? While
I have not looked at the code very closely, a quick observation is
that rTld is never updated during the iterations, which doesn't feel
right intuitively.

Can you help?

Thanks
Peter