Re: odd behaviour with 2x2 matrix floating-point multiplication

[Marco Atzeri]

On 12/16/2014 5:05 PM, Jordi Gutiérrez Hermoso wrote:
> On Tue, 2014-12-16 at 06:37 +0100, Marco Atzeri wrote:
> > additional, for both blas's using the 64 bit versions:
> >
> > ii =  10000
> > count_err =
> >
> >      0
> >      0
> >
> > ans =
> >
> >      0
> >      0
>
> This seems like a dead-giveaway that some of the BLASes is using
> shorter x86 registers, or perhaps is splitting 64-bit doubles to fit
> into 32-bit registers. The x86 FPU uses 80 bits for its internal
> registers, but values in memory are stored as 64-bit doubles. When
> moving values back and forth, precision could be getting lost.

I suspect that on 32 bit there is a mismatch between octave
and blas claculation, while on 64 bit there is no mismatch.
From numerical point of view we can argue forevar which rounding
is the right one..

> What does the following C++ program output when you compile it with
> 32-bit and 64-bit MinGW?
>
>      #include <iostream>
>
>      int main(){
>        double eps = 1.0;
>        while (1 + eps > 1) {
>          eps /= 2;
>        }
>        std::cout << "Machine epsilon is " << 2*eps << std::endl;
>      }

I am not a mingw user. ;-) Just the cygwin package manager.

Neither the original complainer...

> - Jordi G. H.