PRACE NPT

#include <stdio.h>

#include <stdlib.h>

#include <math.h>

#include <mpi.h>

#include "fundefs.h"

#define  MAXNOD 1000000

int me, nodes;

int offset[MAXNOD][2], sizes[MAXNOD][2];

int main( int argc, char* argv[] )

{ 

   MPI_Comm comm = MPI_COMM_WORLD;

   int      m, n, nrep;

   int      mflint, mfltrn, ok, gok;

   int      i, irep, m1, m2, mc, mr, nc, nr, nx;

   double   **arr1r, **arr1i, **arr2r, **arr2i, **carr, **cari;

   double   *ur, *ui, *wr, *wi;

   double   corr, err, frac, mflops, time1, gtime1, time2, gtime2;

   FILE     *inl;

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

   MPI_Init( &argc, &argv );

   MPI_Comm_rank( comm, &me );

   MPI_Comm_size( comm, &nodes );

   if ( me == 0 ) {

      state( "mod2f" );

      prthead( nodes );

   }

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

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

      m   = ilog2( n );

      m2  = m/2;

      m1 = m%2 == 0 ? m2 : m2 + 1;

      nc = pow( 2, m2 );

      nr = pow( 2, m1 );

      sizoff( nr, nc );

      nx  = nc > nr ? nc : nr;

      mr  = nr/nodes;

      mc  = nc/nodes;

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

// --- Check that No. of processes matches the problem size.

      tstinp( mc, mr, nc, nr, me, nodes );

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

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

// --- Allocate partitioned arrays and generate data.

      carr  = makmat( mr, nc ); cari  = makmat( mr, nc );

      arr1r = makmat( mr, nc ); arr1i = makmat( mr, nc );

      arr2r = makmat( mc, nr ); arr2i = makmat( mc, nr );

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

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

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

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

      gendat( mr, nc, carr, cari );

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

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

      time1 = MPI_Wtime();

      time2 = 0.0;

      for( irep = 1; irep <= nrep; irep++ ) {

         cp_arr2d( mr, nc, carr, arr1r );

         cp_arr2d( mr, nc, cari, arr1i );

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

//--- Do 1st transposition.

         gtrans( nr, nc, sizes[me][0], sizes[me][1], 

                 arr1r, arr1i, arr2r, arr2i, 0, &time2 );

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

//--- Do 1st pass of MC NR-length FFTs per processor.

         cfft4( 0, m1, ur, ui, arr2r[0], arr2i[0], wr, wi );

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

            cfft4( 1, m1, ur, ui, arr2r[i], arr2i[i], wr, wi );

         }

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

//--- Multiply with twiddle factors.

         twiddle( mc, nr, arr2r, arr2i );

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

//--- Do 2nd transposition.

         gtrans( nc, nr, sizes[me][1], sizes[me][0],

                 arr2r, arr2i, arr1r, arr1i, 1, &time2 );

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

//--- Do 2nd pass of MR NC-length FFTs per processor.

         cfft4( 0, m2, ur, ui, arr1r[0], arr1i[0], wr, wi );

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

            cfft4( 1, m2, ur, ui, arr1r[i], arr1i[i], wr, wi );

         } 

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

//--- Do 3rd transposition.

         gtrans( nr, nc, sizes[me][0], sizes[me][1], 

                 arr1r, arr1i, arr2r, arr2i, 0, &time2 );

      }

      time1 = MPI_Wtime() - time1;

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

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

      ok   = check( mc, nr, arr2r, arr2i, err );

      corr = MPI_Wtime();

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

         cp_arr2d( mr, nc, carr, arr1r );

         cp_arr2d( mr, nc, cari, arr1i );

      }

      corr  = MPI_Wtime() - corr;

      time1 = time1 - corr;

      MPI_Reduce( &time1, &gtime1, 1, MPI_DOUBLE, MPI_MAX, 0, comm );

      time1 = gtime1/(double)nrep;

      MPI_Reduce( &time2, &gtime2, 1, MPI_DOUBLE, MPI_MAX, 0, comm );

      time2 = gtime2/(double)nrep;

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

// --- Calculate Mflop rates.

      if ( me == 0 ) {

         nflops( m, &mflint, &mfltrn );

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

         frac   = 100.0*(time2/time1);

         prtspeed( n, time1, mflops, time2, frac, ok );

      }

      free( wi ); free( wr );

      free( ui ); free( ur );

      delmat( mc, arr2i ); delmat( mc, arr2r );

      delmat( mr, arr1i ); delmat( mr, arr1r );

      delmat( mr, cari ) ; delmat( mr, carr );

   }

   if ( me == 0 ) {

      printf( "--------------------------------------------------" );

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

   }

   MPI_Finalize();

}