[Bug c/38848] New: Optimizer -O2 doesn't work on linear algebra code on double data type

[Bug c/38848] New: Optimizer -O2 doesn't work on linear algebra code on double data type

[rbenedik at fsmat dot htu dot tuwien dot ac dot at]

/* The program doesn't terminate with -O2 but works with no optimization on
double Datatype */

#include<stdio.h>
#include<stdlib.h>
#include<sys/types.h>
#include<sys/time.h>
#include<unistd.h>

#define STRMAX 1000

#ifdef FLOAT
#define FTYPE(function) function##_f
#define DTYPE float
#define SIZE 4
#endif

#ifdef DOUBLE
#define FTYPE(function) function##_d
#define DTYPE double
#define SIZE 8
#endif

#ifdef EXTENDED
#define FTYPE(function) function##_e
#define DTYPE long double
#define SIZE 16
#endif

void FTYPE(LEAST_SQUARE_QR)(DTYPE *A, DTYPE *x, DTYPE *b, int ZA, int SA)
{
 DTYPE *R, *Qinvb;
 int i,j,k,l;
 DTYPE *ai,an,e,*b1;

 R=(DTYPE *)malloc(ZA*SA*SIZE);
 Qinvb=(DTYPE *)malloc(ZA*SIZE);
 ai=(DTYPE *)malloc(ZA*SIZE);
 b1=(DTYPE *)malloc(ZA*SIZE);

 FTYPE(COPY)(A,R,ZA*SA);
 FTYPE(COPY)(b,Qinvb,ZA);

 i=0;

 for(j=0;j<i;j++)
 {
  *(ai+j)=0.0;
 }
 an=0.0;
 for(j=i;j<ZA;j++)
 {
  *(ai+j)=*(R+j*SA+i);
  an+=*(ai+j) * *(ai+j);
 }
 an=FTYPE(SQUARE_ROOT)(an);

 if(*(ai+i)>0)
 {
  *(ai+i)+=an;
 }
 else
 {
  *(ai+i)-=an;
 }

 an=0.0;
 for(j=i;j<ZA;j++)
 {
  an+= *(ai+j) * *(ai+j);
 }

 an=FTYPE(SQUARE_ROOT)(an);
 an=1.0/an;
for(j=i;j<ZA;j++)
 {
  *(ai+j)=*(ai+j)*an;
 }

 for(l=0;l<SA;l++)
 {
  for(k=0;k<ZA;k++)
  {
   *(b1+k)=0.0;
    for(j=0;j<k;j++)
    {
     e= -2.0* *(ai+j) * *(ai+k);
     *(b1+k)+=e * *(R+j*SA+l);
    }
    e=1.0 - 2.0* *(ai+j) * *(ai+k);
    *(b1+k)+=e * *(R+j*SA+l);
    for(j=k+1;j<ZA;j++)
    {
     e= -2.0* *(ai+j) * *(ai+k);
     *(b1+k)+=e * *(R+j*SA+l);
    }
  }

  for(k=0;k<ZA;k++)
  {
   *(R+k*SA+l)=*(b1+k);
  }
 }

 for(k=1;k<ZA;k++)
 {
  *(R+k*SA)=0.0;
 }

 for(k=0;k<ZA;k++)
 {
  *(b1+k)=0.0;
 }

 for(k=0;k<ZA;k++)
 {
   for(j=0;j<k;j++)
   {
    e = -2.0* *(ai+j) * *(ai+k);
    *(b1+k) +=e * *(Qinvb+j);
   }
   e=1.0 - 2.0* *(ai+j) * *(ai+k);
   *(b1+k)+=e * *(Qinvb+j);
   for(j=k+1;j<ZA;j++)
   {
    e= -2.0* *(ai+j) * *(ai+k);
    *(b1+k)+=e * *(Qinvb+j);
   }
 }
 for(k=0;k<ZA;k++)
 {
  *(Qinvb+k)=*(b1+k);
 }



for(i=1;i<SA;i++)
{
 for(j=0;j<i;j++)
 {
  *(ai+j)=0.0;
 }
 an=0.0;
 for(j=i;j<ZA;j++)
 {
  *(ai+j)=*(R+j*SA+i);
  an+=*(ai+j) * *(ai+j);
 }
 an=FTYPE(SQUARE_ROOT)(an);

 if(*(ai+i)>0)
 {
  *(ai+i)+=an;
 }
 else
 {
  *(ai+i)-=an;
 }




 an=0.0;
 for(j=i;j<ZA;j++)
 {
  an+= *(ai+j) * *(ai+j);
 }

 an=FTYPE(SQUARE_ROOT)(an);
 an=1.0/an;

 for(j=i;j<ZA;j++)
 {
  *(ai+j)=*(ai+j)*an;
 }

 for(l=i;l<SA;l++)
 {
  for(k=0;k<ZA;k++)
  {
   *(b1+k)=0.0;
   for(j=0;j<k;j++)
   {
    e= -2.0* *(ai+j) * *(ai+k);
    *(b1+k)+=e * *(R+j*SA+l);
   }
   e=1.0 - 2.0* *(ai+j) * *(ai+k);
   *(b1+k)+=e * *(R+j*SA+l);
   for(j=k+1;j<ZA;j++)
   {
    e= -2.0* *(ai+j) * *(ai+k);
    *(b1+k)+=e * *(R+j*SA+l);
   }
  }
  for(k=0;k<ZA;k++)
  {
   *(R+k*SA+l)=*(b1+k);
  }
 }


  for(k=0;k<ZA;k++)
  {
   *(b1+k)=0.0;
   for(j=0;j<k;j++)
   {
    e= -2.0* *(ai+j) * *(ai+k);
    *(b1+k)+=e * *(Qinvb+j);
   }
   e=1.0 - 2.0* *(ai+j) * *(ai+k);
   *(b1+k)+=e * *(Qinvb+j);
   for(j=k+1;j<ZA;j++)
   {
    e= -2.0* *(ai+j) * *(ai+k);
    *(b1+k)+=e * *(Qinvb+j);
   }
  }
  for(k=0;k<ZA;k++)
  {
   *(Qinvb+k)=*(b1+k);
  }

 for(k=i+1;k<ZA;k++)
 {
  *(R+k*SA+i)=0.0;
 }
}

for(i=SA;i>=1;i--)
{
 *(x+i-1)=*(Qinvb+i-1);
 for(j=i+1;j<=SA;j++)
 {
  *(x+i-1)-= *(x+j-1) * *(R+SA*(i-1)+j-1);
 }
 *(x+i-1)=*(x+i-1)/ *(R+SA*(i-1)+i-1);
}

free(b1);
free(R);
free(Qinvb);
free(ai);
return;
}


-- 
           Summary: Optimizer -O2 doesn't work on linear algebra code on
                    double data type
           Product: gcc
           Version: 4.3.2
            Status: UNCONFIRMED
          Severity: major
          Priority: P3
         Component: c
        AssignedTo: unassigned at gcc dot gnu dot org
        ReportedBy: rbenedik at fsmat dot htu dot tuwien dot ac dot at
 GCC build triplet: gcc-Version 4.3.2 20081105 (Red Hat 4.3.2-7) (GCC)
  GCC host triplet: x86_64-redhat-linux - Fedora 10.0
GCC target triplet: x86_64-redhat-linux


http://gcc.gnu.org/bugzilla/show_bug.cgi?id=38848

[Bug middle-end/38848] Optimizer -O2 doesn't work on linear algebra code on double data type

[pinskia at gcc dot gnu dot org]

-- 

pinskia at gcc dot gnu dot org changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
           Severity|major                       |normal
          Component|c                           |middle-end


http://gcc.gnu.org/bugzilla/show_bug.cgi?id=38848

[Bug middle-end/38848] Optimizer -O2 doesn't work on linear algebra code on double data type

[rguenth at gcc dot gnu dot org]

------- Comment #1 from rguenth at gcc dot gnu dot org  2009-01-14 21:05 -------
I cannot even compile or link this in any obvious way.

gcc -o t t.c -DDOUBLE
/usr/lib/gcc/i486-linux-gnu/4.3.2/../../../../lib/crt1.o: In function `_start':
(.text+0x18): undefined reference to `main'
/tmp/ccNVE5Zw.o: In function `LEAST_SQUARE_QR_d':
t.c:(.text+0x67): undefined reference to `COPY_d'
t.c:(.text+0x80): undefined reference to `COPY_d'
t.c:(.text+0x111): undefined reference to `SQUARE_ROOT_d'
t.c:(.text+0x1ae): undefined reference to `SQUARE_ROOT_d'
t.c:(.text+0x5db): undefined reference to `SQUARE_ROOT_d'
t.c:(.text+0x678): undefined reference to `SQUARE_ROOT_d'
collect2: ld returned 1 exit status


-- 


http://gcc.gnu.org/bugzilla/show_bug.cgi?id=38848

[Bug middle-end/38848] Optimizer -O2 doesn't work on linear algebra code on double data type

[rbenedik at fsmat dot htu dot tuwien dot ac dot at]

------- Comment #2 from rbenedik at fsmat dot htu dot tuwien dot ac dot at  2009-01-14 21:50 -------
I tried to make a single file from the code but the bug does not show up there.
You can download the software at:

http://fsmat.htu.tuwien.ac.at/~rbenedik/gauss-2.06.tar

Run 

make clean; make
and run

lsq_perf_d 10

in the benchmark directory; float and extended do work on the System: e5200
@2.5 Ghz Fedora 10.0 64 bit double does not.


-- 


http://gcc.gnu.org/bugzilla/show_bug.cgi?id=38848

[Bug middle-end/38848] Optimizer -O2 doesn't work on linear algebra code on double data type

[rguenth at gcc dot gnu dot org]

------- Comment #3 from rguenth at gcc dot gnu dot org  2009-01-14 22:25 -------
You need to do some work to reduce this to a testcase.


-- 

rguenth at gcc dot gnu dot org changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
             Status|UNCONFIRMED                 |WAITING


http://gcc.gnu.org/bugzilla/show_bug.cgi?id=38848

[Bug middle-end/38848] Optimizer -O2 doesn't work on linear algebra code on double data type

[rbenedik at fsmat dot htu dot tuwien dot ac dot at]

------- Comment #4 from rbenedik at fsmat dot htu dot tuwien dot ac dot at  2009-01-15 12:12 -------
I traced the error back to a iteration scheme which is not stable when optmized
with -O2. Sorry this was my fault, so theres no Bug in -O2:

DTYPE FTYPE(SQUARE_ROOT)(DTYPE A)
{
 DTYPE X_N1;
 DTYPE X_N2;

 X_N1=A;

 while(1)
 {
  X_N2 = (DTYPE)(1/2.0) * (X_N1+A/X_N1);
  if(X_N2 == X_N1) break;
  X_N1 = X_N2;
 }

 return X_N2;
}


-- 


http://gcc.gnu.org/bugzilla/show_bug.cgi?id=38848

[Bug middle-end/38848] Optimizer -O2 doesn't work on linear algebra code on double data type

[rguenth at gcc dot gnu dot org]

------- Comment #5 from rguenth at gcc dot gnu dot org  2009-01-15 14:35 -------
Thus invalid.


-- 

rguenth at gcc dot gnu dot org changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
             Status|WAITING                     |RESOLVED
         Resolution|                            |INVALID


http://gcc.gnu.org/bugzilla/show_bug.cgi?id=38848