MIDAPACK - MIcrowave Data Analysis PACKage  1.1b
Parallel software tools for high performance CMB DA analysis
 All Data Structures Files Functions Variables Typedefs Groups Pages
butterfly_extra.c
Go to the documentation of this file.
1 
9 #ifdef W_MPI
10 #include <mpi.h>
11 #include <stdlib.h>
12 #include <string.h>
13 
14 
37 int butterfly_init(int *indices, int count, int **R, int *nR, int **S, int *nS, int **com_indices, int *com_count, int steps, MPI_Comm comm){
38 
39  int i, k, p2k;
40  int rank, size, rk, sk;
41  int tag;
42  MPI_Request s_request, r_request;
43  int nbuf, *buf;
44  int **I, *nI;
45  int **J, *nJ;
46 
47  MPI_Comm_size(comm, &size);
48  MPI_Comm_rank(comm, &rank);
49 
50  I = (int **) malloc(steps * sizeof(int*));
51  nI = (int *) malloc(steps * sizeof(int));
52  tag=0;
53  p2k=size/2;
54 
55  for(k=0; k<steps; k++){ //butterfly first pass : bottom up (fill tabs nI and I)
56  sk=(rank+size-p2k)%size;
57  rk=(rank+p2k)%size;
58 
59  if(k==0){ //S^0 := A
60  nS[k] = count;
61  S[k] = (int *) malloc(nS[k] * sizeof(int));
62  memcpy( S[k], indices, nS[k]*sizeof(int));
63  }
64  else{ //S^k := S^{k-1} \cup R^{k-1}
65  nS[k] = card_or(S[k-1], nS[k-1], I[steps-k], nI[steps-k]);
66  S[k] = (int *) malloc(nS[k] * sizeof(int));
67  set_or(S[k-1], nS[k-1], I[steps-k], nI[steps-k], S[k]);
68  }
69 
70  MPI_Irecv(&nI[steps-k-1], 1, MPI_INT, rk, tag, comm, &r_request); //receive number of indices
71  MPI_Isend(&nS[k], 1, MPI_INT, sk, tag, comm, &s_request); //send number of indices
72  MPI_Wait(&r_request, MPI_STATUS_IGNORE);
73  MPI_Wait(&s_request, MPI_STATUS_IGNORE);
74 
75  I[steps-k-1]= (int *) malloc(nI[steps-k-1] * sizeof(int));
76 
77  tag++;
78  MPI_Irecv(I[steps-k-1], nI[steps-k-1], MPI_INT, rk, tag, comm, &r_request); //receive indices
79  MPI_Isend(S[k], nS[k], MPI_INT, sk, tag, comm, &s_request); //send indices
80  MPI_Wait(&r_request, MPI_STATUS_IGNORE);
81  MPI_Wait(&s_request, MPI_STATUS_IGNORE);
82 
83  p2k/=2;
84  tag++;
85  }
86 
87  J = (int **) malloc(steps * sizeof(int*));
88  nJ = (int *) malloc(steps * sizeof(int));
89 
90  tag=0;
91  p2k=1;
92  for(k=0; k<steps; k++){ //buuterfly second pass : top down (fill tabs nJ and J)
93  free(S[k]);
94  sk=(rank+p2k)%size;
95  rk=(rank+size-p2k)%size;
96  if(k==0){
97  nJ[k] = count;
98  J[k] = (int *) malloc(nJ[k] * sizeof(int));
99  memcpy( J[k], indices, nJ[k]*sizeof(int));
100  }
101  else{
102  nJ[k] = card_or(J[k-1], nJ[k-1], R[k-1], nR[k-1]);
103  J[k] = (int *) malloc(nJ[k] * sizeof(int));
104  set_or(J[k-1], nJ[k-1], R[k-1], nR[k-1], J[k]); //J^k=R^k-1 \cup J^k-1
105  free(R[k-1]);
106  }
107  if(k!=steps-1){
108  MPI_Irecv(&nR[k], 1, MPI_INT, rk, tag, comm, &r_request);
109  MPI_Isend(&nJ[k], 1, MPI_INT, sk, tag, comm, &s_request);
110  MPI_Wait(&r_request, MPI_STATUS_IGNORE);
111  MPI_Wait(&s_request, MPI_STATUS_IGNORE);
112 
113  R[k]= (int *) malloc( nR[k] * sizeof(int));
114  tag++;
115 
116  MPI_Irecv(R[k], nR[k], MPI_INT, rk, tag, comm, &r_request);
117  MPI_Isend(J[k], nJ[k], MPI_INT, sk, tag, comm, &s_request);
118  MPI_Wait(&r_request, MPI_STATUS_IGNORE);
119  MPI_Wait(&s_request, MPI_STATUS_IGNORE);
120  }
121  p2k*=2;
122  tag++;
123  }
124 
125 
126  tag=0;
127  p2k=1;
128  for(k=0; k<steps; k++){ //butterfly last pass : know that Sending tab is S = I \cap J, so send S and we'll get R
129  sk=(rank+p2k)%size;
130  rk=(rank+size-p2k)%size;
131 
132  nS[k] = card_and(I[k], nI[k], J[k], nJ[k]);
133  S[k] = (int *) malloc(nJ[k] *sizeof(int));
134  set_and( I[k], nI[k], J[k], nJ[k], S[k]); //S^k=I^k \cap J^k
135 
136  free(I[k]);
137  free(J[k]);
138 
139  MPI_Irecv(&nR[k],1, MPI_INT, rk, tag, comm, &r_request); //receive size
140  MPI_Isend(&nS[k], 1, MPI_INT, sk, tag, comm, &s_request); //send size
141  MPI_Wait(&r_request, MPI_STATUS_IGNORE);
142  MPI_Wait(&s_request, MPI_STATUS_IGNORE);
143 
144  R[k]= (int *) malloc( nR[k] * sizeof(int));
145  tag++;
146 
147  MPI_Irecv(R[k], nR[k], MPI_INT, rk, tag, comm, &r_request); //receive indices
148  MPI_Isend(S[k], nS[k], MPI_INT, sk, tag, comm, &s_request); //send indices
149  MPI_Wait(&r_request, MPI_STATUS_IGNORE);
150  MPI_Wait(&s_request, MPI_STATUS_IGNORE);
151 
152  p2k*=2;
153  tag++;
154  }
155 
156  //Now we work locally
157  int **USR, *nUSR, **U, *nU;
158 
159  USR = (int **) malloc(steps*sizeof(int *));
160  nUSR = (int *) malloc(steps*sizeof(int));
161  U = (int **) malloc(steps*sizeof(int *));
162  nU = (int *) malloc(steps*sizeof(int));
163 
164  for(k=0; k<steps; k++){
165  nUSR[k] = card_or(S[k], nS[k], R[k], nR[k]);
166  USR[k] = (int *) malloc(nUSR[k]*sizeof(int));
167  set_or(S[k], nS[k], R[k], nR[k], USR[k]);
168  }
169  for(k=0; k<steps; k++){
170  if(k==0){
171  nU[k]=nUSR[k];
172  U[k] = (int *) malloc(nU[k] * sizeof(int));
173  memcpy( U[k], USR[k], nU[k]*sizeof(int));
174  }
175  else{
176  nU[k] = card_or(U[k-1], nU[k-1], USR[k], nUSR[k]);
177  U[k] = (int *) malloc(nU[k]*sizeof(int *));
178  set_or(U[k-1], nU[k-1], USR[k], nUSR[k], U[k]);
179  }
180  }
181  *com_count=nU[steps-1];
182  *com_indices = (int *) malloc(*com_count * sizeof(int));
183  memcpy(*com_indices, U[steps-1], *com_count * sizeof(int));
184  //====================================================================
185 
186  for(k=0; k<steps; k++){
187  subset2map(*com_indices, *com_count, S[k], nS[k]);
188  subset2map(*com_indices, *com_count, R[k], nR[k]);
189  }
190  free(USR);
191  free(U);
192 
193  return 0;
194 }
195 
196 
209 double butterfly_reduce(int **R, int *nR, int nRmax, int **S, int *nS, int nSmax, double *val, int steps, MPI_Comm comm){
210  double st, t;
211  t=0.0;
212  int k, p2k, tag;
213  int rank, size, rk, sk;
214  MPI_Request s_request, r_request;
215  double *sbuf, *rbuf;
216 
217  MPI_Comm_size(comm, &size);
218  MPI_Comm_rank(comm, &rank);
219 
220  sbuf = (double *) malloc(nSmax * sizeof(double));
221  rbuf = (double *) malloc(nRmax * sizeof(double));
222  tag=0;
223  p2k=1;
224 
225  for(k=0; k<steps; k++){
226  sk=(rank+p2k)%size;
227  rk=(rank+size-p2k)%size;
228 
229  m2s(val, sbuf, S[k], nS[k]); //fill the sending buffer
230 
231  st=MPI_Wtime();
232  MPI_Irecv(rbuf, nR[k], MPI_DOUBLE, rk, tag, comm, &r_request);
233  MPI_Isend(sbuf, nS[k], MPI_DOUBLE, sk, tag, comm, &s_request);
234 
235  MPI_Wait(&r_request, MPI_STATUS_IGNORE);
236  MPI_Wait(&s_request, MPI_STATUS_IGNORE);
237 
238  t=t+MPI_Wtime()-st;
239 
240  s2m_sum(val, rbuf, R[k], nR[k]); //sum receive buffer into values
241 
242  p2k*=2;
243  tag++;
244 
245  }
246  free(sbuf);
247  free(rbuf);
248  return t;
249 }
250 
251 int truebutterfly_init(int *indices, int count, int **R, int *nR, int **S, int *nS, int **com_indices, int *com_count, int steps, MPI_Comm comm){
252 
253  int i, k, p2k, p2k1;
254  int rank, size, rk, sk;
255  int tag;
256  MPI_Request s_request, r_request;
257  int nbuf, *buf;
258  int **I, *nI;
259  int **J, *nJ;
260 
261  MPI_Comm_size(comm, &size);
262  MPI_Comm_rank(comm, &rank);
263 
264  I = (int **) malloc(steps * sizeof(int*));
265  nI = (int *) malloc(steps * sizeof(int));
266  tag=0;
267  p2k=size/2;
268  p2k1=2*p2k;
269 
270  for(k=0; k<steps; k++){ //butterfly first pass : bottom up (fill tabs nI and I)
271 
272  if( rank%p2k1 < p2k) sk=rk=rank+p2k; else sk=rk=rank-p2k;
273 
274  if(k==0){ //S^0 := A
275  nS[k] = count;
276  S[k] = (int *) malloc(nS[k] * sizeof(int));
277  memcpy( S[k], indices, nS[k]*sizeof(int));
278  }
279  else{ //S^k := S^{k-1} \cup R^{k-1}
280  nS[k] = card_or(S[k-1], nS[k-1], I[steps-k], nI[steps-k]);
281  S[k] = (int *) malloc(nS[k] * sizeof(int));
282  set_or(S[k-1], nS[k-1], I[steps-k], nI[steps-k], S[k]);
283  }
284 
285  MPI_Irecv(&nI[steps-k-1], 1, MPI_INT, rk, tag, comm, &r_request); //receive number of indices
286  MPI_Isend(&nS[k], 1, MPI_INT, sk, tag, comm, &s_request); //send number of indices
287  MPI_Wait(&r_request, MPI_STATUS_IGNORE);
288  MPI_Wait(&s_request, MPI_STATUS_IGNORE);
289 
290  I[steps-k-1]= (int *) malloc(nI[steps-k-1] * sizeof(int));
291 
292  tag++;
293  MPI_Irecv(I[steps-k-1], nI[steps-k-1], MPI_INT, rk, tag, comm, &r_request); //receive indices
294  MPI_Isend(S[k], nS[k], MPI_INT, sk, tag, comm, &s_request); //send indices
295  MPI_Wait(&r_request, MPI_STATUS_IGNORE);
296  MPI_Wait(&s_request, MPI_STATUS_IGNORE);
297 
298  p2k/=2;
299  p2k1/=2;
300  tag++;
301  }
302 
303  J = (int **) malloc(steps * sizeof(int*));
304  nJ = (int *) malloc(steps * sizeof(int));
305 
306  tag=0;
307  p2k=1;
308  p2k1=p2k*2;
309  for(k=0; k<steps; k++){ //buuterfly second pass : top down (fill tabs nJ and J)
310  free(S[k]);
311 
312  if( rank%p2k1 < p2k) sk=rk=rank+p2k; else sk=rk=rank-p2k;
313 
314  if(k==0){
315  nJ[k] = count;
316  J[k] = (int *) malloc(nJ[k] * sizeof(int));
317  memcpy( J[k], indices, nJ[k]*sizeof(int));
318  }
319  else{
320  nJ[k] = card_or(J[k-1], nJ[k-1], R[k-1], nR[k-1]);
321  J[k] = (int *) malloc(nJ[k] * sizeof(int));
322  set_or(J[k-1], nJ[k-1], R[k-1], nR[k-1], J[k]); //J^k=R^k-1 \cup J^k-1
323  free(R[k-1]);
324  }
325  if(k!=steps-1){
326  MPI_Irecv(&nR[k], 1, MPI_INT, rk, tag, comm, &r_request);
327  MPI_Isend(&nJ[k], 1, MPI_INT, sk, tag, comm, &s_request);
328  MPI_Wait(&r_request, MPI_STATUS_IGNORE);
329  MPI_Wait(&s_request, MPI_STATUS_IGNORE);
330 
331  R[k]= (int *) malloc( nR[k] * sizeof(int));
332  tag++;
333 
334  MPI_Irecv(R[k], nR[k], MPI_INT, rk, tag, comm, &r_request);
335  MPI_Isend(J[k], nJ[k], MPI_INT, sk, tag, comm, &s_request);
336  MPI_Wait(&r_request, MPI_STATUS_IGNORE);
337  MPI_Wait(&s_request, MPI_STATUS_IGNORE);
338  }
339  p2k*=2;
340  p2k1*=2;
341  tag++;
342  }
343 
344 
345  tag=0;
346  p2k=1;
347  p2k1=p2k*2;
348  for(k=0; k<steps; k++){ //butterfly last pass : know that Sending tab is S = I \cap J, so send S and we'll get R
349 
350  if( rank%p2k1 < p2k) sk=rk=rank+p2k; else sk=rk=rank-p2k;
351 
352  nS[k] = card_and(I[k], nI[k], J[k], nJ[k]);
353  S[k] = (int *) malloc(nJ[k] *sizeof(int));
354  set_and( I[k], nI[k], J[k], nJ[k], S[k]); //S^k=I^k \cap J^k
355 
356  free(I[k]);
357  free(J[k]);
358 
359  MPI_Irecv(&nR[k],1, MPI_INT, rk, tag, comm, &r_request); //receive size
360  MPI_Isend(&nS[k], 1, MPI_INT, sk, tag, comm, &s_request); //send size
361  MPI_Wait(&r_request, MPI_STATUS_IGNORE);
362  MPI_Wait(&s_request, MPI_STATUS_IGNORE);
363 
364  R[k]= (int *) malloc( nR[k] * sizeof(int));
365  tag++;
366 
367  MPI_Irecv(R[k], nR[k], MPI_INT, rk, tag, comm, &r_request); //receive indices
368  MPI_Isend(S[k], nS[k], MPI_INT, sk, tag, comm, &s_request); //send indices
369  MPI_Wait(&r_request, MPI_STATUS_IGNORE);
370  MPI_Wait(&s_request, MPI_STATUS_IGNORE);
371 
372  p2k*=2;
373  p2k1*=2;
374  tag++;
375  }
376 
377  //Now we work locally
378  int **USR, *nUSR, **U, *nU;
379 
380  USR = (int **) malloc(steps*sizeof(int *));
381  nUSR = (int *) malloc(steps*sizeof(int));
382  U = (int **) malloc(steps*sizeof(int *));
383  nU = (int *) malloc(steps*sizeof(int));
384 
385  for(k=0; k<steps; k++){
386  nUSR[k] = card_or(S[k], nS[k], R[k], nR[k]);
387  USR[k] = (int *) malloc(nUSR[k]*sizeof(int));
388  set_or(S[k], nS[k], R[k], nR[k], USR[k]);
389  }
390  for(k=0; k<steps; k++){
391  if(k==0){
392  nU[k]=nUSR[k];
393  U[k] = (int *) malloc(nU[k] * sizeof(int));
394  memcpy( U[k], USR[k], nU[k]*sizeof(int));
395  }
396  else{
397  nU[k] = card_or(U[k-1], nU[k-1], USR[k], nUSR[k]);
398  U[k] = (int *) malloc(nU[k]*sizeof(int *));
399  set_or(U[k-1], nU[k-1], USR[k], nUSR[k], U[k]);
400  }
401  }
402  *com_count=nU[steps-1];
403  *com_indices = (int *) malloc(*com_count * sizeof(int));
404  memcpy(*com_indices, U[steps-1], *com_count * sizeof(int));
405  //====================================================================
406 
407  for(k=0; k<steps; k++){
408  subset2map(*com_indices, *com_count, S[k], nS[k]);
409  subset2map(*com_indices, *com_count, R[k], nR[k]);
410  }
411  free(USR);
412  free(U);
413 
414  return 0;
415 }
416 
417 
430 double truebutterfly_reduce(int **R, int *nR, int nRmax, int **S, int *nS, int nSmax, double *val, int steps, MPI_Comm comm){
431  double st, t;
432  t=0.0;
433  int k, p2k, p2k1, tag;
434  int rank, size, rk, sk;
435  MPI_Status status;
436  MPI_Request s_request, r_request;
437  double *sbuf, *rbuf;
438 
439  MPI_Comm_size(comm, &size);
440  MPI_Comm_rank(comm, &rank);
441 
442  sbuf = (double *) malloc(nSmax * sizeof(double));
443  rbuf = (double *) malloc(nRmax * sizeof(double));
444  tag=0;
445  p2k=1;
446  p2k1=p2k*2;
447 
448  for(k=0; k<steps; k++){
449 
450  if( rank%p2k1 < p2k){
451 
452  sk=rk=rank+p2k;
453 
454  st=MPI_Wtime();
455 
456  // MPI_Sendrecv(sbuf, nS[k], MPI_DOUBLE, sk, tag, rbuf, nR[k], MPI_DOUBLE, rk, tag, comm, &status);
457 
458  m2s(val, sbuf, S[k], nS[k]); //fill the sending buffer
459  MPI_Isend(sbuf, nS[k], MPI_DOUBLE, sk, tag, comm, &s_request);
460  MPI_Irecv(rbuf, nR[k], MPI_DOUBLE, rk, tag, comm, &r_request);
461 
462  MPI_Wait(&s_request, MPI_STATUS_IGNORE);
463  MPI_Wait(&r_request, MPI_STATUS_IGNORE);
464  s2m_sum(val, rbuf, R[k], nR[k]); //sum receive buffer into values
465 
466 
467  t=t+MPI_Wtime()-st;
468 
469  } else {
470 
471  sk=rk=rank-p2k;
472 
473  st=MPI_Wtime();
474 
475  MPI_Irecv(rbuf, nR[k], MPI_DOUBLE, rk, tag, comm, &r_request);
476  m2s(val, sbuf, S[k], nS[k]); //fill the sending buffer
477  MPI_Isend(sbuf, nS[k], MPI_DOUBLE, sk, tag, comm, &s_request);
478 
479  MPI_Wait(&r_request, MPI_STATUS_IGNORE);
480  s2m_sum(val, rbuf, R[k], nR[k]); //sum receive buffer into values
481 
482  MPI_Wait(&s_request, MPI_STATUS_IGNORE);
483 
484  // MPI_Sendrecv(sbuf, nS[k], MPI_DOUBLE, sk, tag, rbuf, nR[k], MPI_DOUBLE, rk, tag, comm, &status);
485 
486  t=t+MPI_Wtime()-st;
487 
488  }
489 
490  p2k*=2;
491  p2k1*=2;
492  tag++;
493 
494  }
495  free(sbuf);
496  free(rbuf);
497  return t;
498 }
499 
500 #endif
501 
502