Re: OT: negative zeros with ATLAS 3.8.4.x

[Marco Caliari]

On Mon, 11 Apr 2011, John W. Eaton wrote:

> On 11-Apr-2011, Marco Caliari wrote:
>
> | Dear maintainers,
> |
> | the following code
> |
> | clear all
> | m = 45;
> | y = linspace(0,1,m+1)';
> | y = y(1:m);
> | F = [0:m-1].'*2*pi*y.';
> | F(:,1)
> |
> | produces some strange negative zeros (in random positions) with Octave
> | 3.4.0 and ATLAS 3.8.4.x (x=1 or x=2, Core232SSE3), whereas everything is
> | fine with ATLAS 3.8.3. So, I suspect a bug in ATLAS (uninitialized
> | values?). Maybe someone is able to write a neat C code to submit to ATLAS
> | developers and/or to check that there is surely no fault in Octave.
>
> Can you try to make the expression a little simpler?
>
> Do you get the same result for
>
>  m = 45;
>  y = linspace(0,1,m+1);
>  y = y(1:m);
>  F = [0:m-1].'*2*pi*y;
>  F(:,1)
>
> ?  What about
>
>  m = 45;
>  y = linspace(0,1,m+1);
>  y = y(1:m);
>  t = (0:m-1)'*2*pi;
>  F = t*y;
>  F(:,1)
>
> ?  Or
>
>  m = 45;
>  y = linspace(0,1,m+1);
>  y = y(1:m);
>  t = (0:m-1)';
>  F = t*y;
>  F(:,1)
>
> ?  If you can narrow it down to one matrix multiplication operation
> then it should be simple to write the equivalent Fortran program
> calling LAPACK.

I discovered that the problem is given by the fix described here:

http://math-atlas.sourceforge.net/errata.html#scal0

In few words, if you call the BLAS *scal functions with the scalar set 
to zero, ATLAS set the vector to zeros, whereas standard BLAS multiply the 
vector by zero (and -0*0=-0, NaN*0=NaN, etc.). If I initialize retval with

for (int j = 0; j < a_nr; j++)
  for (int i = 0; i < b_nc; i++)
    retval.xelem (j,i) = 0.;

in dMatrix.cc before the call to dgemm (around line 3180), everything is 
fine with ATLAS 3.8.4.x.

Best regards,

Marco