important. If b_i and (A*guess)_i are zero then BICGSTAB is likely to fail if
(A^(-1)*b)_i is not zero.
implementation.
> 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, <
[hidden email]> 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" <
[hidden email]>
> > An: "QuantLib Mailing Lists" <
[hidden email]>
> > 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
> >
[hidden email]
> >
https://lists.sourceforge.net/lists/listinfo/quantlib-dev>
> ----------------------------------------------------------------------------
> -- _______________________________________________
> QuantLib-dev mailing list
>
[hidden email]
>
https://lists.sourceforge.net/lists/listinfo/quantlib-dev