[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