Logo AND Algorithmique Numérique Distribuée

Public GIT Repository
only ask one version of the standard for pybind11
[simgrid.git] / examples / smpi / NAS / ep.c
1 #include <stdlib.h>
2 #include <stdio.h>
3 #include <string.h>
4 #include <math.h>
5 #include "smpi/mpi.h"
6 #include "nas_common.h"
7 #include "simgrid/instr.h" //TRACE_
8
9 char class;
10 int nprocs;
11
12 #define TRUE 1
13 #define FALSE 0
14
15 int main(int argc, char **argv) {
16   double dum[3] = {1.,1.,1.};
17   double x1, x2, sx, sy, tm, an, gc;
18   XBT_ATTRIB_UNUSED double tt;
19   double Mops;
20   double epsilon=1.0E-8, a = 1220703125., s=271828183.;
21   double t1, t2, t3, t4;
22   double sx_verify_value, sy_verify_value, sx_err, sy_err;
23
24   int    m;
25   int    mk=16;
26   int    nk = (int)(pow(2,mk)),
27          nq=10,
28          np, node, no_nodes, i, ik, kk, l, k, nit, no_large_nodes, np_add, k_offset;
29   int    verified;
30   char   size[500]; // mind the size of the string to represent a big number
31
32   double *x = (double *) SMPI_SHARED_MALLOC (2*nk*sizeof(double));
33   double *q = (double *) SMPI_SHARED_MALLOC (nq*sizeof(double));
34
35   MPI_Init( &argc, &argv );
36   MPI_Comm_size( MPI_COMM_WORLD, &no_nodes);
37   MPI_Comm_rank( MPI_COMM_WORLD, &node);
38
39   TRACE_smpi_set_category ("start");
40
41   get_info(argc, argv, &nprocs, &class);
42   check_info(EP, nprocs, class);
43
44   if      (class == 'S')
45     { m = 24; }
46   else if (class == 'W')
47     { m = 25; }
48   else if (class == 'A')
49     { m = 28; }
50   else if (class == 'B')
51     { m = 30; }
52   else if (class == 'C')
53     { m = 32; }
54   else if (class == 'D')
55     { m = 36; }
56   else if (class == 'E')
57     { m = 40; }
58   else {
59     printf("EP: Internal error: invalid class type %c\n", class);
60     exit(1);
61   }
62
63   int mm = m -mk;
64   int nn = (int)(pow(2,mm));
65
66   int root = 0;
67   if (node == root ) {
68     /* Because the size of the problem is too large to store in a 32-bit integer for some classes, we put it into a
69      * string (for printing). Have to strip off the decimal point put in there by the floating point print statement
70      * (internal file) */
71     fprintf(stdout," NAS Parallel Benchmarks 3.2 -- EP Benchmark");
72     snprintf(size,500,"%lu",(unsigned long)pow(2,m+1));
73     fprintf(stdout," Number of random numbers generated: %s\n",size);
74     fprintf(stdout," Number of active processes: %d\n",no_nodes);
75   }
76
77   /* Compute the number of "batches" of random number pairs generated per processor. Adjust if the number of processors
78    * does not evenly divide the total number  */
79   np = nn / no_nodes;
80   no_large_nodes = nn % no_nodes;
81   if (node < no_large_nodes)
82     np_add = 1;
83   else
84     np_add = 0;
85   np = np + np_add;
86
87   if (np == 0) {
88     fprintf(stdout,"Too many nodes: %d  %d",no_nodes,nn);
89     MPI_Abort(MPI_COMM_WORLD,1);
90     exit(0);
91   }
92
93   /* Call the random number generator functions and initialize the x-array to reduce the effects of paging the timings.
94    Also, call all mathematical functions that are used. Make sure initializations cannot be eliminated as dead code. */
95   vranlc(0, dum[0], dum[1], &(dum[2]));
96
97   dum[0] = randlc(&(dum[1]),&(dum[2]));
98   for (i=0;i<2*nk;i++) {
99     x[i] = -1e99;
100   }
101
102   /* Synchronize before placing time stamp */
103   MPI_Barrier( MPI_COMM_WORLD );
104
105   TRACE_smpi_set_category ("ep");
106
107   timer_clear(1);
108   timer_clear(2);
109   timer_clear(3);
110   timer_start(1);
111
112   t1 = a;
113   t1 = vranlc(0, t1, a, x);
114
115   /* Compute AN = A ^ (2 * NK) (mod 2^46). */
116   t1 = a;
117   for (i=1; i <= mk+1; i++) {
118     t2 = randlc(&t1, &t1);
119   }
120   an = t1;
121   tt = s;
122   gc = 0;
123   tt = 0.;
124   sx = 0.;
125   sy = 0.;
126   for (i=0; i < nq ; i++) {
127     q[i] = 0.;
128   }
129
130 /* Each instance of this loop may be performed independently. We compute the k offsets separately to take into account
131  * the fact that some nodes have more numbers to generate than others */
132   if (np_add == 1)
133     k_offset = node * np -1;
134   else
135     k_offset = no_large_nodes*(np+1) + (node-no_large_nodes)*np -1;
136
137   int stop = FALSE;
138   for(k = 1; k <= np; k++) { SMPI_SAMPLE_GLOBAL(0.25 * np, 0.03) {
139     stop = FALSE;
140     kk = k_offset + k ;
141     t1 = s;
142     t2 = an;
143
144     //       Find starting seed t1 for this kk.
145     for (i=1;i<=100 && !stop;i++) {
146       ik = kk / 2;
147       if (2 * ik != kk)  {
148         randlc(&t1, &t2);
149       }
150       if (ik==0)
151         stop = TRUE;
152       else {
153         randlc(&t2, &t2);
154         kk = ik;
155       }
156     }
157     //       Compute uniform pseudorandom numbers.
158
159     timer_start(3);
160
161     t1 = vranlc(2 * nk, t1, a, x);
162
163     timer_stop(3);
164
165     /* Compute Gaussian deviates by acceptance-rejection method and tally counts in concentric square annuli.
166      * This loop is not vectorizable. */
167     timer_start(2);
168
169     for(i=1; i<=nk;i++) {
170       x1 = 2. * x[2*i-2] -1.0;
171       x2 = 2. * x[2*i-1] - 1.0;
172       t1 = x1*x1 + x2*x2;
173       if (t1 <= 1.) {
174         t2   = sqrt(-2. * log(t1) / t1);
175         t3   = (x1 * t2);
176         t4   = (x2 * t2);
177         l    = (int)(fabs(t3) > fabs(t4) ? fabs(t3) : fabs(t4));
178         q[l] = q[l] + 1.;
179         sx   = sx + t3;
180         sy   = sy + t4;
181       }
182     }
183     timer_stop(2);
184     } }
185
186   TRACE_smpi_set_category ("finalize");
187
188   MPI_Allreduce(&sx, x, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
189   sx = x[0];
190   MPI_Allreduce(&sy, x, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
191   sy = x[0];
192   MPI_Allreduce(q, x, nq, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
193
194   for(i = 0; i < nq; i++) {
195     q[i] = x[i];
196   }
197   for(i = 0; i < nq; i++) {
198     gc += q[i];
199   }
200
201   timer_stop(1);
202   tm = timer_read(1);
203   MPI_Allreduce(&tm, x, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
204   tm = x[0];
205
206   SMPI_SHARED_FREE(x);
207
208   if(node == root) {
209     nit = 0;
210     verified = TRUE;
211
212     if(m == 24) {
213       sx_verify_value = -3.247834652034740E3;
214       sy_verify_value = -6.958407078382297E3;
215     } else if(m == 25) {
216       sx_verify_value = -2.863319731645753E3;
217       sy_verify_value = -6.320053679109499E3;
218     } else if(m == 28) {
219       sx_verify_value = -4.295875165629892E3;
220       sy_verify_value = -1.580732573678431E4;
221     } else if(m == 30) {
222       sx_verify_value =  4.033815542441498E4;
223       sy_verify_value = -2.660669192809235E4;
224     } else if(m == 32) {
225       sx_verify_value =  4.764367927995374E4;
226       sy_verify_value = -8.084072988043731E4;
227     } else if(m == 36) {
228       sx_verify_value =  1.982481200946593E5;
229       sy_verify_value = -1.020596636361769E5;
230     } else {
231       verified = FALSE;
232     }
233
234     if(verified) {
235       sx_err = fabs((sx - sx_verify_value)/sx_verify_value);
236       sy_err = fabs((sy - sy_verify_value)/sy_verify_value);
237       verified = ((sx_err < epsilon) && (sy_err < epsilon));
238     }
239
240     Mops = (pow(2.0, m+1))/tm/1000;
241
242     fprintf(stdout,"EP Benchmark Results:\n");
243     fprintf(stdout,"CPU Time=%d\n",(int) tm);
244     fprintf(stdout,"N = 2^%d\n",m);
245     fprintf(stdout,"No. Gaussain Pairs =%d\n",(int) gc);
246     fprintf(stdout,"Sum = %f %ld\n",sx,(long) sy);
247     fprintf(stdout,"Count:");
248     for(i = 0; i < nq; i++) {
249       fprintf(stdout,"%d\t %ld\n",i,(long) q[i]);
250     }
251     c_print_results("EP", class, m+1, 0, 0, nit, nprocs, no_nodes, tm, Mops, "Random number generated",verified);
252
253     fprintf(stdout,"Total time:     %f\n",(timer_read(1)/1000));
254     fprintf(stdout,"Gaussian pairs: %f\n",(timer_read(2)/1000));
255     fprintf(stdout,"Random numbers: %f\n",(timer_read(3)/1000));
256   }
257
258   SMPI_SHARED_FREE(q);
259
260   MPI_Finalize();
261   return 0;
262 }