[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 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 */

#includestdio.h
#includestdlib.h
#includesys/types.h
#includesys/time.h
#includeunistd.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;ji;j++)
 {
  *(ai+j)=0.0;
 }
 an=0.0;
 for(j=i;jZA;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;jZA;j++)
 {
  an+= *(ai+j) * *(ai+j);
 }

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

 for(l=0;lSA;l++)
 {
  for(k=0;kZA;k++)
  {
   *(b1+k)=0.0;
for(j=0;jk;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;jZA;j++)
{
 e= -2.0* *(ai+j) * *(ai+k);
 *(b1+k)+=e * *(R+j*SA+l);
}
  }

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

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

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

 for(k=0;kZA;k++)
 {
   for(j=0;jk;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;jZA;j++)
   {
e= -2.0* *(ai+j) * *(ai+k);
*(b1+k)+=e * *(Qinvb+j);
   }
 }
 for(k=0;kZA;k++)
 {
  *(Qinvb+k)=*(b1+k);
 }



for(i=1;iSA;i++)
{
 for(j=0;ji;j++)
 {
  *(ai+j)=0.0;
 }
 an=0.0;
 for(j=i;jZA;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;jZA;j++)
 {
  an+= *(ai+j) * *(ai+j);
 }

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

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

 for(l=i;lSA;l++)
 {
  for(k=0;kZA;k++)
  {
   *(b1+k)=0.0;
   for(j=0;jk;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;jZA;j++)
   {
e= -2.0* *(ai+j) * *(ai+k);
*(b1+k)+=e * *(R+j*SA+l);
   }
  }
  for(k=0;kZA;k++)
  {
   *(R+k*SA+l)=*(b1+k);
  }
 }


  for(k=0;kZA;k++)
  {
   *(b1+k)=0.0;
   for(j=0;jk;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;jZA;j++)
   {
e= -2.0* *(ai+j) * *(ai+k);
*(b1+k)+=e * *(Qinvb+j);
   }
  }
  for(k=0;kZA;k++)
  {
   *(Qinvb+k)=*(b1+k);
  }

 for(k=i+1;kZA;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

[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