"Fossies" - the Fresh Open Source Software Archive

Member "hpcc-1.5.0/FFT/tstfft.c" (18 Mar 2016, 3871 Bytes) of package /linux/privat/hpcc-1.5.0.tar.gz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "tstfft.c" see the Fossies "Dox" file reference documentation and the latest Fossies "Diffs" side-by-side code changes report: 1.5.0b_vs_1.5.0.

    1 /* -*- mode: C; tab-width: 2; indent-tabs-mode: nil; fill-column: 79; coding: iso-latin-1-unix -*- */
    2 /* tstfft.c
    3  */
    4 
    5 #include <hpcc.h>
    6 
    7 #include "hpccfft.h"
    8 
    9 static int
   10 TestFFT1(HPCC_Params *params, int doIO, FILE *outFile, double *UGflops, int *Un, int *Ufailure) {
   11   fftw_complex *in, *out;
   12   fftw_plan p;
   13   hpcc_fftw_plan ip;
   14   double Gflops = -1.0;
   15   double maxErr, tmp1, tmp2, tmp3, t0, t1, t2, t3;
   16   int i, n, flags, failure = 1;
   17   double deps = HPL_dlamch( HPL_MACH_EPS );
   18 
   19 #ifdef HPCC_FFT_235
   20   int f[3];
   21 
   22   /* Need 2 vectors for input and output and 1 vector of scratch spaces */
   23   n = HPCC_LocalVectorSize( params, 3, sizeof(fftw_complex), 0 );
   24 
   25   /* Adjust local size for factors */
   26   for ( ; HPCC_factor235( n, f ); n--)
   27     ; /* EMPTY */
   28 #else
   29   /* Need 2 vectors and vectors' sizes as power of 2 */
   30   n = HPCC_LocalVectorSize( params, 2, sizeof(fftw_complex), 1 );
   31 #endif
   32 
   33   /* need to use fftw_malloc() so that the returned pointers will be aligned properly for SSE
   34      instructions on Intel/AMD systems */
   35   in  = (fftw_complex *)HPCC_fftw_malloc( (sizeof *in)  * n );
   36   out = (fftw_complex *)HPCC_fftw_malloc( (sizeof *out) * n );
   37 
   38   if (! in || ! out) goto comp_end;
   39 
   40   /* Make sure that `inout' and `work' are initialized in parallel if using
   41      Open MP: this will ensure better placement of pages if first-touch policy
   42      is used by a distrubuted shared memory machine. */
   43 #ifdef _OPENMP
   44 #pragma omp parallel for
   45   for (i = 0; i < n; ++i) {
   46     c_re( in[i] ) = c_re( out[i] ) = 0.0;
   47     c_re( in[i] ) = c_im( out[i] ) = 0.0;
   48   }
   49 #endif
   50 
   51   t0 = -MPI_Wtime();
   52   HPCC_bcnrand( 2*n, 0, in );
   53   t0 += MPI_Wtime();
   54 
   55 #ifdef HPCC_FFTW_ESTIMATE
   56   flags = FFTW_ESTIMATE;
   57 #else
   58   flags = FFTW_MEASURE;
   59 #endif
   60 
   61   t1 = -MPI_Wtime();
   62   p = fftw_create_plan( n, FFTW_FORWARD, flags );
   63   t1 += MPI_Wtime();
   64 
   65   if (! p) goto comp_end;
   66 
   67   t2 = -MPI_Wtime();
   68   fftw_one( p, in, out );
   69   t2 += MPI_Wtime();
   70 
   71   fftw_destroy_plan(p);
   72 
   73   ip = HPCC_fftw_create_plan( n, FFTW_BACKWARD, FFTW_ESTIMATE );
   74 
   75   if (ip) {
   76     t3 = -MPI_Wtime();
   77     HPCC_fftw_one( ip, out, in );
   78     t3 += MPI_Wtime();
   79 
   80     HPCC_fftw_destroy_plan( ip );
   81   }
   82 
   83   HPCC_bcnrand( 2*(s64Int)n, 0, out ); /* regenerate data */
   84   maxErr = 0.0;
   85   for (i = 0; i < n; i++) {
   86     tmp1 = c_re( in[i] ) - c_re( out[i] );
   87     tmp2 = c_im( in[i] ) - c_im( out[i] );
   88     tmp3 = sqrt( tmp1*tmp1 + tmp2*tmp2 );
   89     maxErr = maxErr >= tmp3 ? maxErr : tmp3;
   90   }
   91 
   92   if (maxErr / log(n) / deps < params->test.thrsh) failure = 0;
   93 
   94   if (doIO) {
   95     fprintf( outFile, "Vector size: %d\n", n );
   96     fprintf( outFile, "Generation time: %9.3f\n", t0 );
   97     fprintf( outFile, "Tuning: %9.3f\n", t1 );
   98     fprintf( outFile, "Computing: %9.3f\n", t2 );
   99     fprintf( outFile, "Inverse FFT: %9.3f\n", t3 );
  100     fprintf( outFile, "max(|x-x0|): %9.3e\n", maxErr );
  101   }
  102 
  103   if (t2 > 0.0) Gflops = 1e-9 * (5.0 * n * log(n) / log(2.0)) / t2;
  104 
  105   comp_end:
  106 
  107   if (out) HPCC_fftw_free( out );
  108   if (in)  HPCC_fftw_free( in );
  109 
  110   *UGflops = Gflops;
  111   *Un = n;
  112   *Ufailure = failure;
  113 
  114   return 0;
  115 }
  116 
  117 int
  118 HPCC_TestFFT(HPCC_Params *params, int doIO, double *UGflops, int *Un, int *Ufailure) {
  119   int rv, n, failure = 1;
  120   double Gflops;
  121   FILE *outFile = NULL;
  122 
  123   if (doIO) {
  124     outFile = fopen( params->outFname, "a" );
  125     if (! outFile) {
  126       outFile = stderr;
  127       fprintf( outFile, "Cannot open output file.\n" );
  128       return 1;
  129     }
  130   } else {
  131     outFile = fopen( "/dev/null", "w" ); /* special filename Unix file systems */
  132     if (! outFile) {
  133       outFile = fopen( "nul", "w"); /* special filename on Windows, produces no output */
  134     }
  135   }
  136 
  137   n = 0;
  138   Gflops = -1.0;
  139   rv = TestFFT1( params, doIO, outFile, &Gflops, &n, &failure );
  140 
  141   if (doIO) {
  142     fflush( outFile );
  143   }
  144 
  145   if (outFile)
  146     fclose( outFile );
  147 
  148   if (UGflops) *UGflops = Gflops;
  149   if (Un) *Un = n;
  150   if (Ufailure) *Ufailure = failure;
  151 
  152   return rv;
  153 }