Re: BiCGstab

Posted by Peter Caspers-4 on Nov 20, 2015; 6:15am
URL: http://quantlib.414.s1.nabble.com/BiCGstab-tp17097p17101.html

Hi Klaus,

thank you. Actually the breakdown test case was born because I forgot
to provide the guess in the solve(...) method (and then 0 is taken).
With the guess provided (and the fix discussed earlier) I did not run
into any problems up to now, also in bigger and more realistic / less
trivial cases for A and b.

In the Numerical Recipes (2.7) I found this statement

"While there is no guarantee that this whole procedure will not break
down or become
unstable for general A, in practice this is rare."

So everything's fine, sorry for my slight disbelief in the
implementation, it was based on ignorance :-)

Thanks again for your help
Peter


On 19 November 2015 at 21:31, Klaus Spanderen <[hidden email]> wrote:

> 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,  <[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
>

------------------------------------------------------------------------------
_______________________________________________
QuantLib-dev mailing list
[hidden email]
https://lists.sourceforge.net/lists/listinfo/quantlib-dev