Re: [Quantlib-dev] BiCGstab

[Klaus S. <kl...@sp...>]

Hi Peter,

rTld  does not get updated, which makes the choice of the guess even more 
important. If b_i and (A*guess)_i are zero then BICGSTAB is likely to fail if 
(A^(-1)*b)_i is not zero. 

MatLab also fails to solve your second example using MatLab's own BICGSTAB 
implementation. 

regards
Klaus

On Thursday, November 19, 2015 03:55:54 PM Peter Caspers wrote:
> 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
> 
> ----------------------------------------------------------------------------
> -- _______________________________________________
> QuantLib-dev mailing list
> Qua...@li...
> https://lists.sourceforge.net/lists/listinfo/quantlib-dev