Used astyle --style=kr formatted source codes.
[charm.git] / example / fft-trans / fft_ref.cpp
1 #include <mpi.h> 
2 #include <cstdlib>
3 #include "fftw3-mpi.h"
4 #include <cmath>
5 #include <limits>
6 #include "fileio.h"
7
8 using namespace std;
9
10 int main(int argc, char *argv[]) {
11   int rank, size; 
12   MPI_Init(&argc, &argv);
13   MPI_Comm_size(MPI_COMM_WORLD, &size);
14   MPI_Comm_rank(MPI_COMM_WORLD, &rank); 
15
16   fftw_plan plan;
17   fftw_complex *data;
18
19   fftw_mpi_init();
20
21   if(rank == 0) {
22     if(argc != 2) {
23       printf("Usage: ./binary <N>\n");
24       MPI_Abort(MPI_COMM_WORLD,-1);
25     }
26   }
27
28   int N = atoi(argv[1]);
29
30   ptrdiff_t local_ni=N*N/size, local_i_start = N*N/size*rank;
31   ptrdiff_t local_no=local_ni, local_o_start = local_i_start;
32
33   int b_or_f = FFTW_BACKWARD;
34
35   ptrdiff_t alloc_local =  fftw_mpi_local_size_1d(N*N, MPI_COMM_WORLD,
36       b_or_f, FFTW_ESTIMATE, &local_ni, &local_i_start,
37       &local_no, &local_o_start);
38
39   data = fftw_alloc_complex(alloc_local);
40
41   plan = fftw_mpi_plan_dft_1d(N*N, data, data, MPI_COMM_WORLD, b_or_f, FFTW_ESTIMATE);
42
43   char filename[80];
44   sprintf(filename, "%d-%d.dump%d", size, N, rank);
45   readCommFile(data, filename);
46
47   fftw_execute(plan);
48
49   double infNorm = 0.0;
50   srand48(rank);
51   for(int i = 0; i < N*N/size; i++) {
52     data[i][0] = data[i][0]/(N*N) - drand48();
53     data[i][1] = data[i][1]/(N*N) - drand48();
54
55     double mag = sqrt(pow(data[i][0], 2) + pow(data[i][1], 2));
56     if(mag > infNorm) infNorm = mag;
57   }
58
59   double my_r = infNorm / (std::numeric_limits<double>::epsilon() * log((double)N * N));
60   double r;
61
62   MPI_Reduce(&my_r, &r, 1, MPI_DOUBLE, MPI_MAX, 0, MPI_COMM_WORLD);
63
64   if(rank == 0) {
65     if(r < 16)
66       printf("r = %g, PASS!\n",r);
67     else
68       printf("r = %g, FAIL\n",r);
69   }
70
71   fftw_destroy_plan(plan);
72
73   fftw_mpi_cleanup();
74   MPI_Finalize();      
75   return 0; 
76