[SciPy-Dev] TestNorm.test_stable failure [was ANN: SciPy 0.11.0 release candidate 2]

[Rüdiger Kessel]

Don't we talk about snrm2() for float32? So it should always return
float32 and not float64, doesn't it?

Ruediger

2012/9/22 Ralf Gommers <ralf.gommers at gmail.com>:
> 57711ac8b
>
>
>
> On Sat, Sep 22, 2012 at 6:14 PM, Rüdiger Kessel <ruediger.kessel at gmail.com>
> wrote:
>>
>> Hi,
>>
>> test_basic.TestNorm.test_stable is very sensitive against the
>> implementation of the basic BLAS routine snrm2() as I have explained
>> before.
>> Its good to tell implementations apart, but it is not very robust.
>>
>> A pure float32 implementation of snrm2() will give 10000.0 for norm(a).
>> An internal double implementation with rounding to float32 at the end
>> gives 10000.5 for norm(a).
>> A mixed implementation with rounding to float32 before rescaling gives
>> 10000.4990234375 for norm(a).
>>
>> I have not looked into the implementation in ATLAS. I just played with
>> my python implementation to see if I can find an explanation.
>>
>> I think the problem here is that for matrix a some data was chosen
>> where a pure float32 implementation of snrm2() will not give a precise
>> result. So any implementation with gives something in the range from
>> 10000.0 to 10000.5 need to be considered as correct.
>>
>> Maybe it would be better to chose data where rounded double and
>> float32 implementation will give the same result like:
>>
>> a = array(range(90000,100000), dtype=float32)/10.0
>>
>> norm(a) should give 950434.0 for any implementation.
>
>
> I think you mean 950433.5, which is what float64 gives (approximately). On
> win32 both scipy and numpy give 950433.5, on OS X scipy gives 950433.5 and
> numpy 950433.56. The original purpose of the test was to check that the
> scipy version was more stable than the numpy one, so the input you suggest
> may work.
>
> Thanks,
> Ralf
>
>
>
> _______________________________________________
> SciPy-Dev mailing list
> SciPy-Dev at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-dev
>