Index: test-suite/trunk/MultiSource/Benchmarks/MiBench/telecomm-FFT/CMakeLists.txt =================================================================== --- test-suite/trunk/MultiSource/Benchmarks/MiBench/telecomm-FFT/CMakeLists.txt +++ test-suite/trunk/MultiSource/Benchmarks/MiBench/telecomm-FFT/CMakeLists.txt @@ -3,4 +3,5 @@ set(RUN_OPTIONS 8 32768 -i) set(FP_ABSTOLERANCE 0.00001) set(HASH_PROGRAM_OUTPUT 1) +add_definitions(-DFP_ABSTOLERANCE=1e-5) llvm_multisource() Index: test-suite/trunk/MultiSource/Benchmarks/MiBench/telecomm-FFT/Makefile =================================================================== --- test-suite/trunk/MultiSource/Benchmarks/MiBench/telecomm-FFT/Makefile +++ test-suite/trunk/MultiSource/Benchmarks/MiBench/telecomm-FFT/Makefile @@ -2,6 +2,7 @@ PROG = telecomm-fft LDFLAGS = -lm +CFLAGS += -DFP_ABSTOLERANCE=1e-5 RUN_OPTIONS = 8 32768 -i FP_ABSTOLERANCE := 0.00001 HASH_PROGRAM_OUTPUT = 1 Index: test-suite/trunk/MultiSource/Benchmarks/MiBench/telecomm-FFT/fourier.h =================================================================== --- test-suite/trunk/MultiSource/Benchmarks/MiBench/telecomm-FFT/fourier.h +++ test-suite/trunk/MultiSource/Benchmarks/MiBench/telecomm-FFT/fourier.h @@ -45,6 +45,13 @@ 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 ); Index: test-suite/trunk/MultiSource/Benchmarks/MiBench/telecomm-FFT/fourierf.c =================================================================== --- test-suite/trunk/MultiSource/Benchmarks/MiBench/telecomm-FFT/fourierf.c +++ test-suite/trunk/MultiSource/Benchmarks/MiBench/telecomm-FFT/fourierf.c @@ -146,5 +146,116 @@ } } +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 ---*/ Index: test-suite/trunk/MultiSource/Benchmarks/MiBench/telecomm-FFT/main.c =================================================================== --- test-suite/trunk/MultiSource/Benchmarks/MiBench/telecomm-FFT/main.c +++ test-suite/trunk/MultiSource/Benchmarks/MiBench/telecomm-FFT/main.c @@ -21,6 +21,20 @@ } // 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; @@ -29,6 +43,8 @@ float *ImagIn; float *RealOut; float *ImagOut; + float *RealOut_StrictFP; + float *ImagOut_StrictFP; float *coeff; float *amp; int invfft=0; @@ -51,6 +67,8 @@ 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); @@ -81,21 +99,30 @@ /* 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