[test-suite] fix telecomm-FFT for CFLAGS="-ffp-contract=on"

[sebpop]

Summary

This patch fixes MultiSource/Benchmarks/MiBench/telecomm-FFT when compiled with -ffp-contract=on.

[rengolin]

Same as previous, LGTM. Thanks!

test-suite/trunk/MultiSource/Benchmarks/MiBench/telecomm-FFT/CMakeLists.txt

 set(PROG telecomm-fft)
 list(APPEND LDFLAGS -lm)
 set(RUN_OPTIONS 8 32768 -i)
 set(FP_ABSTOLERANCE 0.00001)
 set(HASH_PROGRAM_OUTPUT 1)
+add_definitions(-DFP_ABSTOLERANCE=1e-5)
 llvm_multisource()

test-suite/trunk/MultiSource/Benchmarks/MiBench/telecomm-FFT/Makefile

 LEVEL = ../../../..
 
 PROG     = telecomm-fft
 LDFLAGS  = -lm
+CFLAGS += -DFP_ABSTOLERANCE=1e-5
 RUN_OPTIONS = 8 32768 -i
 FP_ABSTOLERANCE :=  0.00001
 HASH_PROGRAM_OUTPUT = 1
 
 include $(LEVEL)/MultiSource/Makefile.multisrc

test-suite/trunk/MultiSource/Benchmarks/MiBench/telecomm-FFT/fourier.h

 void fft_float (
     unsigned  NumSamples,          /* must be a power of 2 */
     int       InverseTransform,    /* 0=forward FFT, 1=inverse FFT */
     float    *RealIn,              /* array of input's real samples */
     float    *ImaginaryIn,         /* array of input's imag samples */
     float    *RealOut,             /* array of output's reals */
     float    *ImaginaryOut );      /* array of output's imaginaries */
 
+void fft_float_StrictFP (
+    unsigned  NumSamples,          /* must be a power of 2 */
+    int       InverseTransform,    /* 0=forward FFT, 1=inverse FFT */
+    float    *RealIn,              /* array of input's real samples */
+    float    *ImaginaryIn,         /* array of input's imag samples */
+    float    *RealOut,             /* array of output's reals */
+    float    *ImaginaryOut );      /* array of output's imaginaries */
 
 int IsPowerOfTwo ( unsigned x );
 unsigned NumberOfBitsNeeded ( unsigned PowerOfTwo );
 unsigned ReverseBits ( unsigned index, unsigned NumBits );
 
 /*
 **   The following function returns an "abstract frequency" of a
 **   given index into a buffer with a given number of frequency samples.

test-suite/trunk/MultiSource/Benchmarks/MiBench/telecomm-FFT/fourierf.c

if ( InverseTransform )
         for ( i=0; i < NumSamples; i++ )
         {
             RealOut[i] /= denom;
             ImagOut[i] /= denom;
         }
     }
 }
 
+void fft_float_StrictFP (
+    unsigned  NumSamples,
+    int       InverseTransform,
+    float    *RealIn,
+    float    *ImagIn,
+    float    *RealOut,
+    float    *ImagOut )
+{
+#pragma STDC FP_CONTRACT OFF
+    unsigned NumBits;    /* Number of bits needed to store indices */
+    unsigned i, j, k, n;
+    unsigned BlockSize, BlockEnd;
+
+    double angle_numerator = 2.0 * DDC_PI;
+    double tr, ti;     /* temp real, temp imaginary */
+
+    if ( !IsPowerOfTwo(NumSamples) )
+    {
+        fprintf (
+            stderr,
+            "Error in fft():  NumSamples=%u is not power of two\n",
+            NumSamples );
+
+        exit(1);
+    }
+
+    if ( InverseTransform )
+        angle_numerator = -angle_numerator;
+
+    CHECKPOINTER ( RealIn );
+    CHECKPOINTER ( RealOut );
+    CHECKPOINTER ( ImagOut );
+
+    NumBits = NumberOfBitsNeeded ( NumSamples );
+
+    /*
+    **   Do simultaneous data copy and bit-reversal ordering into outputs...
+    */
+
+    for ( i=0; i < NumSamples; i++ )
+    {
+        j = ReverseBits ( i, NumBits );
+        RealOut[j] = RealIn[i];
+        ImagOut[j] = (ImagIn == NULL) ? 0.0 : ImagIn[i];
+    }
+
+    /*
+    **   Do the FFT itself...
+    */
+
+    BlockEnd = 1;
+    for ( BlockSize = 2; BlockSize <= NumSamples; BlockSize <<= 1 )
+    {
+        double delta_angle = angle_numerator / (double)BlockSize;
+        double sm2 = sin ( -2 * delta_angle );
+        double sm1 = sin ( -delta_angle );
+        double cm2 = cos ( -2 * delta_angle );
+        double cm1 = cos ( -delta_angle );
+        double w = 2 * cm1;
+        double ar[3], ai[3];
+        double temp;
+
+        for ( i=0; i < NumSamples; i += BlockSize )
+        {
+            ar[2] = cm2;
+            ar[1] = cm1;
+
+            ai[2] = sm2;
+            ai[1] = sm1;
+
+            for ( j=i, n=0; n < BlockEnd; j++, n++ )
+            {
+                ar[0] = w*ar[1] - ar[2];
+                ar[2] = ar[1];
+                ar[1] = ar[0];
+
+                ai[0] = w*ai[1] - ai[2];
+                ai[2] = ai[1];
+                ai[1] = ai[0];
+
+                k = j + BlockEnd;
+                tr = ar[0]*RealOut[k] - ai[0]*ImagOut[k];
+                ti = ar[0]*ImagOut[k] + ai[0]*RealOut[k];
+
+                RealOut[k] = RealOut[j] - tr;
+                ImagOut[k] = ImagOut[j] - ti;
+
+                RealOut[j] += tr;
+                ImagOut[j] += ti;
+            }
+        }
+
+        BlockEnd = BlockSize;
+    }
+
+    /*
+    **   Need to normalize if inverse transform...
+    */
+
+    if ( InverseTransform )
+    {
+        double denom = (double)NumSamples;
+
+        for ( i=0; i < NumSamples; i++ )
+        {
+            RealOut[i] /= denom;
+            ImagOut[i] /= denom;
+        }
+    }
+}
+
 
 /*--- end of file fourierf.c ---*/

test-suite/trunk/MultiSource/Benchmarks/MiBench/telecomm-FFT/main.c

int rand( void ) {
     return (unsigned int)(next / 65536) % RAND_MAX+1;
 }
 
 void srand( unsigned int seed ) {
     next = seed;
 }
 // End of RNG implementation
 
+/* Return 0 when V1 and V2 do not match within the allowed FP_ABSTOLERANCE.  */
+static inline int
+check_FP(float V1, float V2) {
+  double AbsTolerance = FP_ABSTOLERANCE;
+  double Diff = fabs(V1 - V2);
+  if (Diff > AbsTolerance) {
+    fprintf(stderr, "A = %lf and B = %lf differ more than"
+                  " FP_ABSTOLERANCE = %lf\n", V1, V2, AbsTolerance);
+    return 0;
+  }
+
+  return 1;
+}
+
 int main(int argc, char *argv[]) {
   unsigned MAXSIZE;
   unsigned MAXWAVES;
   unsigned i,j;
   float *RealIn;
   float *ImagIn;
   float *RealOut;
   float *ImagOut;
+  float *RealOut_StrictFP;
+  float *ImagOut_StrictFP;
   float *coeff;
   float *amp;
   int invfft=0;
 
   if (argc<3)
   {
     printf("Usage: fft <waves> <length> -i\n");
     printf("-i performs an inverse fft\n");
     printf("make <waves> random sinusoids");
     printf("<length> is the number of samples\n");
     exit(-1);
   }
   else if (argc==4)
     invfft = !strncmp(argv[3],"-i",2);
   MAXSIZE=atoi(argv[2]);
   MAXWAVES=atoi(argv[1]);
   srand(1);
     
  RealIn=(float*)malloc(sizeof(float)*MAXSIZE);
  ImagIn=(float*)malloc(sizeof(float)*MAXSIZE);
  RealOut=(float*)malloc(sizeof(float)*MAXSIZE);
  ImagOut=(float*)malloc(sizeof(float)*MAXSIZE);
+ RealOut_StrictFP=(float*)malloc(sizeof(float)*MAXSIZE);
+ ImagOut_StrictFP=(float*)malloc(sizeof(float)*MAXSIZE);
  coeff=(float*)malloc(sizeof(float)*MAXWAVES);
  amp=(float*)malloc(sizeof(float)*MAXWAVES);
 
  /* Makes MAXWAVES waves of random amplitude and period */
   for(i=0;i<MAXWAVES;i++) 
   {
     coeff[i] = i%1000;
     amp[i] = i%1000;
else
        RealIn[i]+=coeff[j]*sin(amp[j]*i);
      }
      ImagIn[i]=0;
    }
  }
 
  /* regular*/
  fft_float (MAXSIZE,invfft,RealIn,ImagIn,RealOut,ImagOut);
+ fft_float_StrictFP (MAXSIZE,invfft,RealIn,ImagIn,RealOut_StrictFP,ImagOut_StrictFP);
  
  printf("RealOut:\n");
- for (i=0;i<MAXSIZE;i++)
-   printf("%f \t", RealOut[i]);
+ for (i=0;i<MAXSIZE;i++) {
+   if (!check_FP(RealOut[i], RealOut_StrictFP[i]))
+     return 1;
+   printf("%f \t", RealOut_StrictFP[i]);
+ }
  printf("\n");
 
 printf("ImagOut:\n");
- for (i=0;i<MAXSIZE;i++)
-   printf("%f \t", ImagOut[i]);
-   printf("\n");
+ for (i=0;i<MAXSIZE;i++) {
+   if (!check_FP(ImagOut[i], ImagOut_StrictFP[i]))
+     return 1;
+   printf("%f \t", ImagOut_StrictFP[i]);
+ }
 
  free(RealIn);
  free(ImagIn);
  free(RealOut);
  free(ImagOut);
+ free(RealOut_StrictFP);
+ free(ImagOut_StrictFP);
  free(coeff);
  free(amp);
  exit(0);
 
 
 }