Statistics
| Revision:

root / Pi / C / Hybrid / Pi_Hybrid.c @ 183

History | View | Annotate | Download (8.6 kB)

1 64 equemene
//
2 64 equemene
// Estimation of Pi using Monte Carlo exploration process
3 64 equemene
// gcc -std=c99 -O3 -o Pi Pi.c -lm
4 64 equemene
// Emmanuel Quemener <emmanuel.quemener@ens-lyon.fr>
5 64 equemene
// Cecill v2
6 64 equemene
7 64 equemene
// Needed for gethostname
8 64 equemene
#define _BSD_SOURCE
9 64 equemene
#include <sys/unistd.h>
10 64 equemene
11 64 equemene
#include <math.h>
12 64 equemene
#include <stdio.h>
13 64 equemene
#include <stdlib.h>
14 64 equemene
#include <limits.h>
15 64 equemene
#include <mpi.h>
16 64 equemene
#include <stddef.h>
17 64 equemene
18 64 equemene
#ifdef TIME
19 64 equemene
#include <sys/time.h>
20 64 equemene
#endif
21 64 equemene
22 64 equemene
// Marsaglia RNG very simple implementation
23 64 equemene
#define znew  ((z=36969*(z&65535)+(z>>16))<<16)
24 64 equemene
#define wnew  ((w=18000*(w&65535)+(w>>16))&65535)
25 64 equemene
#define MWC   (znew+wnew)
26 64 equemene
#define SHR3  (jsr=(jsr=(jsr=jsr^(jsr<<17))^(jsr>>13))^(jsr<<5))
27 64 equemene
#define CONG  (jcong=69069*jcong+1234567)
28 64 equemene
#define KISS  ((MWC^CONG)+SHR3)
29 64 equemene
30 79 equemene
#define ITERATIONS 1000000000
31 79 equemene
32 64 equemene
#define MWCfp MWC * 2.328306435454494e-10f
33 64 equemene
#define KISSfp KISS * 2.328306435454494e-10f
34 79 equemene
#define SHR3fp SHR3 * 2.328306435454494e-10f
35 79 equemene
#define CONGfp CONG * 2.328306435454494e-10f
36 64 equemene
37 64 equemene
#define PROCESS 1
38 64 equemene
39 64 equemene
#ifdef LONG
40 64 equemene
#define LENGTH long long
41 64 equemene
#else
42 64 equemene
#define LENGTH int
43 64 equemene
#endif
44 64 equemene
45 64 equemene
typedef struct compute_node {
46 183 equemene
  LENGTH iterations;
47 183 equemene
  int process;
48 64 equemene
} node;
49 64 equemene
50 183 equemene
typedef struct compute_result {
51 183 equemene
  LENGTH inside;
52 183 equemene
  long int useconds;
53 125 equemene
} result;
54 125 equemene
55 64 equemene
unsigned int rotl(unsigned int value, int shift) {
56 183 equemene
  return (value << shift) | (value >> (sizeof(value) * CHAR_BIT - shift));
57 64 equemene
}
58 183 equemene
59 64 equemene
unsigned int rotr(unsigned int value, int shift) {
60 183 equemene
  return (value >> shift) | (value << (sizeof(value) * CHAR_BIT - shift));
61 64 equemene
}
62 64 equemene
63 64 equemene
LENGTH MainLoopGlobal(LENGTH iterations,unsigned int seed_w,unsigned int seed_z)
64 64 equemene
{
65 183 equemene
66 79 equemene
#if defined TCONG
67 183 equemene
  unsigned int jcong=seed_z;
68 79 equemene
#elif defined TSHR3
69 183 equemene
  unsigned int jsr=seed_w;
70 79 equemene
#elif defined TMWC
71 183 equemene
  unsigned int z=seed_z;
72 183 equemene
  unsigned int w=seed_w;
73 79 equemene
#elif defined TKISS
74 183 equemene
  unsigned int jcong=seed_z;
75 183 equemene
  unsigned int jsr=seed_w;
76 183 equemene
  unsigned int z=seed_z;
77 183 equemene
  unsigned int w=seed_w;
78 79 equemene
#endif
79 79 equemene
80 183 equemene
  LENGTH total=0;
81 183 equemene
82 183 equemene
  for (LENGTH i=0;i<iterations;i++) {
83 183 equemene
84 79 equemene
#if defined TINT32
85 183 equemene
#define THEONE 1073741824
86 183 equemene
#if defined TCONG
87 183 equemene
    unsigned int x=CONG>>17 ;
88 183 equemene
    unsigned int y=CONG>>17 ;
89 183 equemene
#elif defined TSHR3
90 183 equemene
    unsigned int x=SHR3>>17 ;
91 183 equemene
    unsigned int y=SHR3>>17 ;
92 183 equemene
#elif defined TMWC
93 183 equemene
    unsigned int x=MWC>>17 ;
94 183 equemene
    unsigned int y=MWC>>17 ;
95 183 equemene
#elif defined TKISS
96 183 equemene
    unsigned int x=KISS>>17 ;
97 183 equemene
    unsigned int y=KISS>>17 ;
98 183 equemene
#endif
99 79 equemene
#elif defined TINT64
100 183 equemene
#define THEONE 4611686018427387904
101 183 equemene
#if defined TCONG
102 183 equemene
    unsigned long x=(unsigned long)(CONG>>1) ;
103 183 equemene
    unsigned long y=(unsigned long)(CONG>>1) ;
104 183 equemene
#elif defined TSHR3
105 183 equemene
    unsigned long x=(unsigned long)(SHR3>>1) ;
106 183 equemene
    unsigned long y=(unsigned long)(SHR3>>1) ;
107 183 equemene
#elif defined TMWC
108 183 equemene
    unsigned long x=(unsigned long)(MWC>>1) ;
109 183 equemene
    unsigned long y=(unsigned long)(MWC>>1) ;
110 183 equemene
#elif defined TKISS
111 183 equemene
    unsigned long x=(unsigned long)(KISS>>1) ;
112 183 equemene
    unsigned long y=(unsigned long)(KISS>>1) ;
113 183 equemene
#endif
114 79 equemene
#elif defined TFP32
115 183 equemene
#define THEONE 1.0f
116 183 equemene
#if defined TCONG
117 183 equemene
    float x=CONGfp ;
118 183 equemene
    float y=CONGfp ;
119 183 equemene
#elif defined TSHR3
120 183 equemene
    float x=SHR3fp ;
121 183 equemene
    float y=SHR3fp ;
122 183 equemene
#elif defined TMWC
123 183 equemene
    float x=MWCfp ;
124 183 equemene
    float y=MWCfp ;
125 183 equemene
#elif defined TKISS
126 183 equemene
    float x=KISSfp ;
127 183 equemene
    float y=KISSfp ;
128 183 equemene
#endif
129 79 equemene
#elif defined TFP64
130 183 equemene
#define THEONE 1.0f
131 183 equemene
#if defined TCONG
132 183 equemene
    double x=(double)CONGfp ;
133 183 equemene
    double y=(double)CONGfp ;
134 183 equemene
#elif defined TSHR3
135 183 equemene
    double x=(double)SHR3fp ;
136 183 equemene
    double y=(double)SHR3fp ;
137 183 equemene
#elif defined TMWC
138 183 equemene
    double x=(double)MWCfp ;
139 183 equemene
    double y=(double)MWCfp ;
140 183 equemene
#elif defined TKISS
141 183 equemene
    double x=(double)KISSfp ;
142 183 equemene
    double y=(double)KISSfp ;
143 79 equemene
#endif
144 183 equemene
#endif
145 64 equemene
146 183 equemene
    // Matching test
147 183 equemene
    unsigned long inside=((x*x+y*y) < THEONE) ? 1:0;
148 183 equemene
    total+=inside;
149 183 equemene
  }
150 183 equemene
151 183 equemene
  return(total);
152 183 equemene
153 64 equemene
}
154 64 equemene
155 64 equemene
int main(int argc, char *argv[]) {
156 183 equemene
157 64 equemene
  unsigned int seed_z=362436069,seed_w=52128862,process=PROCESS;
158 64 equemene
  // Number of NP or OpenMP processes <1024
159 125 equemene
  LENGTH iterations=ITERATIONS,insideMPI[8192],insideOpenMP[1024],
160 183 equemene
    part_inside=0,part_iterations,insides=0;
161 64 equemene
  int numtasks,rank,rc,tag=1,i;
162 64 equemene
  float pi;
163 64 equemene
164 64 equemene
  // Hostname supposed to be <128 characters
165 64 equemene
  char hostname[128];
166 183 equemene
167 64 equemene
  gethostname(hostname, sizeof hostname);
168 183 equemene
169 64 equemene
  struct timeval start,end;
170 64 equemene
  long int useconds;
171 183 equemene
172 64 equemene
  MPI_Status Stat;
173 76 equemene
174 64 equemene
  rc = MPI_Init(&argc,&argv);
175 64 equemene
  if (rc != MPI_SUCCESS) {
176 64 equemene
    printf ("Error starting MPI program. Terminating.\n");
177 64 equemene
    MPI_Abort(MPI_COMM_WORLD, rc);
178 64 equemene
  }
179 64 equemene
180 64 equemene
  MPI_Comm_size(MPI_COMM_WORLD,&numtasks);
181 125 equemene
  MPI_Comm_rank(MPI_COMM_WORLD,&rank);
182 183 equemene
183 125 equemene
  const int nitems=2;
184 125 equemene
  int blocklengths[2] = {1,1};
185 183 equemene
186 64 equemene
#ifdef LONG
187 125 equemene
  MPI_Datatype types_node[2] = {MPI_LONG, MPI_INT};
188 125 equemene
  MPI_Datatype types_result[2] = {MPI_LONG, MPI_LONG};
189 64 equemene
#else
190 125 equemene
  MPI_Datatype types_node[2] = {MPI_INT, MPI_INT};
191 125 equemene
  MPI_Datatype types_result[2] = {MPI_INT, MPI_LONG};
192 64 equemene
#endif
193 183 equemene
194 125 equemene
  MPI_Datatype mpi_node_type,mpi_result_type;
195 125 equemene
  MPI_Aint     offsets[2],offsetsr[2];
196 183 equemene
197 125 equemene
  offsets[0] = offsetof(node, iterations);
198 125 equemene
  offsets[1] = offsetof(node, process);
199 183 equemene
200 183 equemene
  MPI_Type_create_struct(nitems, blocklengths, offsets, types_node, &mpi_node_type);
201 125 equemene
  MPI_Type_commit(&mpi_node_type);
202 183 equemene
203 125 equemene
  offsetsr[0] = offsetof(result, inside);
204 125 equemene
  offsetsr[1] = offsetof(result, useconds);
205 183 equemene
206 183 equemene
  MPI_Type_create_struct(nitems, blocklengths, offsetsr, types_result, &mpi_result_type);
207 125 equemene
  MPI_Type_commit(&mpi_result_type);
208 183 equemene
209 64 equemene
  if (rank==0) {
210 64 equemene
211 64 equemene
    if (argc > 1) {
212 64 equemene
      iterations=(LENGTH)atoll(argv[1]);
213 64 equemene
      process=atoi(argv[2]);
214 64 equemene
    }
215 64 equemene
    else {
216 64 equemene
      printf("\n\tPi : Estimate Pi with Monte Carlo exploration\n\n");
217 183 equemene
      printf("\t\t#1 : number of iterations (default 1 billion)\n");
218 183 equemene
      printf("\t\t#2 : number of OpenMP processes (default 1)\n\n");
219 64 equemene
    }
220 64 equemene
221 64 equemene
    printf ("\n\tInformation about architecture:\n\n");
222 64 equemene
223 183 equemene
    printf ("\tSizeof int = %lld bytes.\n", (long long)sizeof(int));
224 183 equemene
    printf ("\tSizeof long = %lld bytes.\n", (long long)sizeof(long));
225 183 equemene
    printf ("\tSizeof long long = %lld bytes.\n", (long long)sizeof(long long));
226 64 equemene
227 183 equemene
    printf ("\tMax int = %u\n", INT_MAX);
228 183 equemene
    printf ("\tMax long = %ld\n", LONG_MAX);
229 183 equemene
    printf ("\tMax long long = %lld\n\n", LLONG_MAX);
230 64 equemene
231 183 equemene
    part_iterations=(((iterations%numtasks)%process) == 0) ? iterations/numtasks/process:iterations/numtasks/process+1 ;
232 183 equemene
233 64 equemene
    node send;
234 64 equemene
    send.iterations=part_iterations;
235 64 equemene
    send.process=process;
236 64 equemene
237 64 equemene
    // Split part of code
238 76 equemene
    for (i=1;i<numtasks;i++) {
239 124 equemene
      rc = MPI_Send(&send, 1, mpi_node_type, i, tag, MPI_COMM_WORLD);
240 64 equemene
    }
241 64 equemene
242 64 equemene
    gettimeofday(&start,(struct timezone *)0);
243 64 equemene
244 64 equemene
#pragma omp parallel for
245 183 equemene
    for (int i=0 ; i<process; i++) {
246 183 equemene
      insideOpenMP[i]=MainLoopGlobal(part_iterations,
247 183 equemene
                                     rotr(seed_w,i),
248 183 equemene
                                     rotl(seed_z,i));
249 183 equemene
      /*
250 183 equemene
        printf("\t(%s,%i) found %lld for process %i\n",hostname,0,
251 183 equemene
        (long long)insideOpenMP[i],i); */
252 183 equemene
    }
253 183 equemene
    /*
254 183 equemene
      printf("\n");
255 183 equemene
    */
256 125 equemene
257 183 equemene
    insides=0;
258 183 equemene
    for (int i=0 ; i<process; i++) {
259 183 equemene
      insides+=insideOpenMP[i];
260 183 equemene
    }
261 125 equemene
262 183 equemene
    gettimeofday(&end,(struct timezone *)0);
263 183 equemene
    useconds=(end.tv_sec-start.tv_sec)*1000000+end.tv_usec-start.tv_usec;
264 183 equemene
265 183 equemene
    printf("\tOn %s with rank #%i find %lld inside in %lu useconds.\n",
266 183 equemene
           hostname,rank,(long long)insides,useconds);
267 183 equemene
268 183 equemene
    // Join part of code
269 183 equemene
    for (i=1;i<numtasks;i++) {
270 64 equemene
271 183 equemene
      result recv;
272 64 equemene
273 183 equemene
      rc = MPI_Recv(&recv, 1, mpi_result_type, i, tag, MPI_COMM_WORLD,&Stat);
274 183 equemene
275 183 equemene
      insideMPI[i]=recv.inside;
276 183 equemene
      useconds=recv.useconds;
277 183 equemene
278 183 equemene
      printf("\tReceive from rank #%i, find %lld inside in %lu useconds\n",i,(long long)insideMPI[i],useconds);
279 183 equemene
280 183 equemene
      insides+=insideMPI[i];
281 183 equemene
    }
282 183 equemene
283 183 equemene
    pi=4.*(float)insides/(float)(part_iterations*numtasks*process);
284 183 equemene
285 183 equemene
    printf("\n\tPi=%.40f\n\twith error %.40f\n\twith %lld iterations\n\n",pi,
286 183 equemene
           fabs(pi-4*atan(1.))/pi,(long long)(part_iterations*numtasks*process));
287 183 equemene
288 64 equemene
  }
289 64 equemene
  else
290 64 equemene
    {
291 64 equemene
      // Receive information from master
292 183 equemene
293 64 equemene
      node recv;
294 183 equemene
295 124 equemene
      rc = MPI_Recv(&recv, 1, mpi_node_type, 0, tag, MPI_COMM_WORLD,&Stat);
296 183 equemene
      /*
297 64 equemene
      printf("\t(%s,%i) receive from master %lld with %i process\n",
298 183 equemene
      hostname,rank,(long long)recv.iterations,recv.process);
299 125 equemene
      */
300 183 equemene
301 64 equemene
      gettimeofday(&start,(struct timezone *)0);
302 183 equemene
303 64 equemene
#pragma omp parallel for
304 125 equemene
      for (int i=0 ; i<recv.process; i++) {
305 183 equemene
        insideOpenMP[i]=MainLoopGlobal(recv.iterations,rotr(seed_w,rank+i),rotl(seed_z,rank-i));
306 125 equemene
        /*
307 183 equemene
          printf("\t(%s,%i) found %lld for process %i\n",hostname,rank,
308 183 equemene
          (long long)insideOpenMP[i],i);
309 125 equemene
        */
310 125 equemene
      }
311 64 equemene
312 125 equemene
      /* printf("\n"); */
313 183 equemene
314 125 equemene
      for (int i=0 ; i<recv.process; i++) {
315 125 equemene
        part_inside+=insideOpenMP[i];
316 125 equemene
      }
317 125 equemene
318 64 equemene
      gettimeofday(&end,(struct timezone *)0);
319 64 equemene
      useconds=(end.tv_sec-start.tv_sec)*1000000+end.tv_usec-start.tv_usec;
320 125 equemene
      /*
321 183 equemene
        printf("\tOn %s rank %i find %lld inside in %lu useconds.\n",
322 183 equemene
        hostname,rank,(long long)part_inside,useconds);
323 125 equemene
      */
324 125 equemene
      result send;
325 125 equemene
      send.inside=part_inside;
326 125 equemene
      send.useconds=useconds;
327 64 equemene
328 125 equemene
      rc = MPI_Send(&send, 1, mpi_result_type, 0, tag, MPI_COMM_WORLD);
329 183 equemene
330 64 equemene
    }
331 64 equemene
332 64 equemene
  MPI_Type_free(&mpi_node_type);
333 125 equemene
  MPI_Type_free(&mpi_result_type);
334 183 equemene
335 64 equemene
  MPI_Finalize();
336 64 equemene
337 64 equemene
}