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)

Even simpler:

clear all
m = 45;
y = (0:m-1);
t = (0:m-1)';
F = t*y;
F(:,1)

I observe *almost* always two negative zeros (last and third-last
position). But with the following Fortran code

      program test
      implicit none
      integer i
      double precision a(45,1),b(1,45),c(45,45)
      do i = 1,45
         a(i,1) = (i-1) * 1.0d0
         b(1,i) = (i-1) * 1.0d0
      enddo
      call dgemm('N','N',45,45,1,1.0d0,a,45,b,1,0.0d0,c,45)
      do i = 1,45
         write(6,*) c(i,1)
      enddo
      stop
      end

(linked to the threaded ATLAS libraries and with the same compiler options 
I used to build Octave) I don't see negative zeros.

Marco