[SciPy-Dev] ANN: SciPy 0.11.0 release candidate 1

[Rüdiger Kessel]

Andreas Hilboll <lists <at> hilboll.de> writes:

> 
> > 20.07.2012 11:34, Andreas Hilboll kirjoitti:
> > [clip]
> >> ======================================================================
> >> FAIL: test_basic.TestNorm.test_stable
> >> ----------------------------------------------------------------------
> > [clip]
> >
> > Broken BLAS or something? Try this fortran program:
> >
> >
> > !!! test.f90
> > program main
> >     implicit none
> >     real a(10001)
> >     real b, snrm2
> >     external snrm2
> >
> >     a(1) = 1e4
> >     a(2:10001) = 1.0
> >
> >     b = snrm2(10001, a, 1)
> >     write(*,*) b   ! should give 10000.5
> > end program main
> >
> >
> 
> % gfortran -lblas test_blas.f90
> % ./a.out
>    10000.0000
> 
> So apparently, something's wrong. But what?
> 
> Cheers, A.
> 

Hi,

I think BLAS is right.

I implemented the core of the original BLAS routine
(http://www.netlib.org/blas/snrm2.f) in python:

def snrm2(X,INCX=1):
    import math
    from numpy import float32
    N=len(X)
    SCALE = float32(0.0)
    SSQ = float32(0.0)
    for ix in range(0,N,INCX):
        if X[ix]!=0.0:
            ABSXI = float32(abs(X[ix]))
            if SCALE < ABSXI:
                SSQ = float32(1.0) + SSQ * float32(float32(SCALE/ABSXI)**2)
                SCALE = ABSXI
            else:
                SSQ = SSQ + float32(float32(ABSXI/SCALE)**2)
    NORM = SCALE*float32(math.sqrt(SSQ))
    return NORM

It shows the same behavior.

The reason is that the ratio between 10000**2 / 1**2 is 1e8. So the loop is
adding a small value to a large value in single precision.

float32(1e8+1) = 1e8

This can only be improved if snrm2 would internally use double precision. But
this was obviously not the intention.

Greetings
Rüdiger