PRACE NPT

#include <stdio.h>

#include <stdlib.h>

#include <math.h>

double cclock( void );

int    check( double a[], double b[], int m, double err );

void   cfft4( int, int, double ur[], double ui[], double xr[], double xi[],

              double wr[], double wi[] );

void   gendat( double *a, int );

void   nflops( int, int *, int * );

void   prthead( void );

void   prtspeed( int, double, double, int );

void   state( char * );

int main()

{ 

   int      m, n, nrep;

   int      mflint, mfltrn, ok;

   int      i, j;

   double   *a, *b, *c, *ur, *ui, *wr, *wi;

   double   corr, err, mflops, time;

   FILE     *inl;

// ------------------------------------------------------------------------

   state( "mod2f" );

   prthead();

   inl = fopen( "mod2f.in", "r" );

   while( ( fscanf( inl, "%d%d\n", &m, &nrep ) != EOF ) ){

      n   = pow( 2, m );

      err = ( 10.0*n*m )* 1.0e-10;         // --- Allowed error tolerance.

      a   = calloc( n, sizeof( double ) );

      b   = calloc( n, sizeof( double ) );

      c   = calloc( n, sizeof( double ) );

      ur  = calloc( n, sizeof( double ) );

      ui  = calloc( n, sizeof( double ) );

      wr  = calloc( n, sizeof( double ) );

      wi  = calloc( n, sizeof( double ) );

// -------------------------------------------------------------------------

// --- Generate data.

      gendat( c, n );

// -------------------------------------------------------------------------

// --- Repeat FFT 'nrep' times for this problem size.

      time = cclock();

      for( i = 0; i < nrep; i++ ) {

         for( j = 0; j < n; j++ ) {

            a[j] = c[j];

            b[j] = 0.0;

         }

         cfft4( 0, m, ur, ui, a, b, wr, wi );

         cfft4( 1, m, ur, ui, a, b, wr, wi );

      }

      time = cclock() - time;

// -------------------------------------------------------------------------

// --- Check for errors and correct timing for filling of arrays.

      ok   = check( a, b, n, err );

      corr = cclock();

      for( i = 0; i < nrep; i++ ) {

         for( j = 0; j < n; j++ ) {

            a[j] = c[j];

            b[j] = 0.0;

         }

      }

      corr = cclock() - corr;

      time = ( time - corr )/(double)nrep;

// -------------------------------------------------------------------------

// --- Calculate Mflop rates.

      nflops( m, &mflint, &mfltrn );

      mflops = 1.0e-6*( mflint + mfltrn )/time;

      prtspeed( n, time, mflops, ok );

      free( wi ); free( wr ); free( ui ); free( ur );

      free( c );  free( b ); free( a );

   }

   printf( "-----------------------------------------------\n" );

}