Re: [Quantlib-dev] BiCGstab

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

Hi Ralph,

thanks a lot. I will send Luigi a PR with the fix for the error initialization.

Best regards
Peter


On 19 November 2015 at 15:03,  <ori...@we...> wrote:
> Hi Peter,
>
> I've just had a quick look ... The proper initialization of "error" seems to
> be missing and changing L52 to
>
> Real error=norm2(r)/bnorm2;
>
> as you suggested would be the correct bugfix.
>
> On the other issue, the intermediate results rho == 0.0 or omega == 0.0 are
> so-called breakdown cases for the algorithm. To cope with that, one would
> have to combine the BiCG with other methods ... In the BiCGstab algorithm,
> rTld is not updated indeed.
>
>
> Best regards
>
> Ralph
>
>
>
>
> Gesendet: Mittwoch, 18. November 2015 um 19:53 Uhr
> Von: "Peter Caspers" <pca...@gm...>
> An: "QuantLib Mailing Lists" <qua...@li...>
> Betreff: [Quantlib-dev] BiCGstab
> 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
>
> ------------------------------------------------------------------------------
> _______________________________________________
> QuantLib-dev mailing list
> Qua...@li...
> https://lists.sourceforge.net/lists/listinfo/quantlib-dev