[Insight-developers] Numerical issues with new code

[Bill Lorensen]

Actually, you don't need the borland compiler to see the issues. Look
at the output of one of the tests and you will see many nan's.

I'm curious, who converted the fortran to c++, someone at Stanford?

Bill

On Wed, Apr 14, 2010 at 12:16 PM, Luis Ibanez <luis.ibanez at kitware.com> wrote:
>
> Bill,
>
> Good call !
>
> The Valgrind errors were indeed very enlightening.
>
>
> Here is what seems to happen:
>
> The array "w[n]" is a "working space", declared in
> line 338, of lsqrBase.cxx:
>
>                 double * w = new double[n];
>
> and only gets to be initialized if the variable "alpha"
> is larger than zero, in lines 364-368:
>
>   if( alpha > zero )
>     {
>     this->Scale( n, ( one / alpha ), v );
>     CopyVector( n, v, w );
>     }
>
> where "alpha" is the norm of he vector "v".
>
> In the original Fortran code, this is:
>
>     if (alpha > zero) then
>        call dscal (n, (one/alpha), v, 1)
>        w = v
>     end if
>
>
> To put this in context, from the paper:
> http://www.stanford.edu/group/SOL/software/lsqr/lsqr-toms82a.pdf
>
> The solver is dealing with the problem
>
>                           A*x = b
>
> where the matrix A is known,
> the vector b is known,
> and we are looking for the vector x.
>
> One of the first things that the solver does
> is to scale the problem by defining a vector "u"
> such that
>
>                      beta * u = b,
>
> and a vector "v" such that
>
>                  alpha * v = A' * u.
>
> where A' is the transpose of A.
>
> So,
>
> beta is the norm of the input vector "b",
>
> while
>
> alpha is the norm of the vector resulting
> from the product : A' * u.
>
>
> ---
>
> This also shows that I mis-translated a condition
> for early termination, that corresponds to the case
> of "alpha == 0" or "beta == 0".
>
>
> I'm correction that section at this point....
>
>
>
>        Luis
>
>
> -------------------------------------------------------------------------------------------------
> On Wed, Apr 14, 2010 at 8:27 AM, Bill Lorensen <bill.lorensen at gmail.com>
> wrote:
>>
>> Luis,
>>
>> Maybe these valgrind errors will help:
>>
>> http://www.cdash.org/CDash/viewDynamicAnalysis.php?buildid=586308
>>
>> Bill
>>
>> On Tue, Apr 13, 2010 at 4:18 PM, Luis Ibanez <luis.ibanez at kitware.com>
>> wrote:
>> >
>> > Hi Bill,
>> >
>> > Yeap, that seems to be the case.
>> >
>> > I'm updating my Borland build now...
>> >
>> >
>> >     Thanks for pointing this out.
>> >
>> >
>> >           Luis
>> >
>> > ---------------------------------------------------
>> > On Mon, Apr 12, 2010 at 3:15 PM, Bill Lorensen <bill.lorensen at gmail.com>
>> > wrote:
>> >>
>> >> Luis,
>> >>
>> >> Looks like there are some numerical problems with the new
>> >> vnl_lsqr code.
>> >>
>> >> Borland is reporting
>> >> http://www.cdash.org/CDash/viewTest.php?onlyfailed&buildid=584274
>> >>
>> >> I think Borland is the only compiler that currently catches floating
>> >> point divide by zero, etc.
>> >>
>> >> This is probably due to a real error in the new code.
>> >>
>> >> Bill
>> >
>> >
>
>