Pastebin launched a little side project called VERYVIRAL.com, check it out ;-) Want more features on Pastebin? Sign Up, it's FREE!
Guest

mpi red black sor

By: a guest on Jan 24th, 2012  |  syntax: C  |  size: 5.96 KB  |  views: 431  |  expires: Never
download  |  raw  |  embed  |  report abuse  |  print
Text below is selected. Please press Ctrl+C to copy to your clipboard. (⌘+C on Mac)
  1. #include <stdio.h>
  2. #include <stdlib.h>
  3. #include <math.h>
  4. #include <string.h>
  5. #include <sys/time.h>
  6. #include <mpi.h>
  7.  
  8. #define MIN(X,Y) ((X) < (Y) ? (X) : (Y))
  9.  
  10. #define TOLERANCE 0.00002        /* termination criterion */
  11. #define MAGIC 0.8                /* magic factor */
  12.  
  13. int N; /* problem size */
  14. double **G; /* the grid */
  15. double **buff;
  16. MPI_Status istatus;
  17. MPI_Request ireqs1, ireqs2, ireqr1, ireqr2;
  18. int mynode, totalnodes;
  19. int jsta, jend, inext, iprev;
  20. double *bufs1, *bufr1;
  21. double *bufs2, *bufr2;
  22. double stopdiff, maxdiff, diff; /* differences btw grid points in iters */
  23.  
  24. int even (int i)
  25. {
  26.         return !(i & 1);
  27. }
  28.  
  29. double stencil (double** G, int row, int col)
  30. {
  31.         return (G[row-1][col] + G[row+1][col] + G[row][col-1] + G[row][col+1] ) / 4.0;
  32. }
  33.  
  34. void alloc_grid(double ***Gptr, int N)
  35. {
  36.         int i;
  37.         double** G = (double**)malloc(N*sizeof(double*));
  38.         if ( G == 0 ) {
  39.                 fprintf(stderr, "malloc failed\n");
  40.                 exit(42);
  41.         }
  42.  
  43.         for (i = 0; i<N; i++) { /* malloc the own range plus one more line */
  44.                 /* of overlap on each side */
  45.                 G[i] = (double*)malloc(N*sizeof(double));
  46.                 if ( G[i] == 0 ) {
  47.                         fprintf(stderr, "malloc failed\n");
  48.                         exit(42);
  49.                 }
  50.         }
  51.  
  52.         *Gptr = G;
  53. }
  54.  
  55. void init_grid(double **G, int N)
  56. {
  57.         int i, j;
  58.  
  59.         /* initialize the grid */
  60.         for (i = 0; i < N; i++){
  61.                 for (j = 0; j < N; j++){
  62.                         if (i == 0) G[i][j] = 4.56;
  63.                         else if (i == N-1) G[i][j] = 9.85;
  64.                         else if (j == 0) G[i][j] = 7.32;
  65.                         else if (j == N-1) G[i][j] = 6.88;
  66.                         else G[i][j] = 0.0;
  67.                 }
  68.         }
  69. }
  70.  
  71. void print_grid(double** G, int N)
  72. {
  73.         int i, j;
  74.  
  75.         for ( i = 1 ; i < N-1 ; i++ ) {
  76.                 for ( j = 1 ; j < N-1 ; j++ ) {
  77.                         printf("%10.3f ", G[i][j]);
  78.                 }
  79.                 printf("\n");
  80.         }
  81. }
  82.  
  83. void range(int n1, int n2, int nprocs, int irank, int *ista, int *iend) {
  84.     float iwork1 = (n2 - n1 + 1) / nprocs;
  85.     float iwork2 = (n2 - n1 + 1) % nprocs;
  86.     int start, end;
  87.  
  88.     start = (int) (irank * iwork1 + n1 + MIN(irank, iwork2));
  89.     end = (int) (start + iwork1 - 1);
  90.  
  91.     if (iwork2 > irank) {
  92.         end = end + 1;
  93.     }
  94.  
  95.     *ista = start;
  96.     *iend = end;
  97. }
  98.  
  99. void exchange(int phase) {
  100.     int is1 = ((jsta + phase) % 2) + 1;
  101.     int is2 = ((jend + phase) % 2) + 1;
  102.     int ir1 = fabs(3 - is1);
  103.     int ir2 = fabs(3 - is2);
  104.     int icnt = 0;
  105.     int icnt1 = 0, icnt2 = 0;
  106.     int m = N - 1;
  107.     int i;
  108.  
  109.     if (mynode != 0) {
  110.         icnt1 = 0;
  111.         for (i = is1; i <= m; i += 2) {
  112.             bufs1[icnt1] = G[i][jsta];
  113.             icnt1 = icnt1 + 1;
  114.         }
  115.     }
  116.  
  117.     if (mynode != (totalnodes - 1)) {
  118.         icnt2 = 0;
  119.         for (i = is2; i <= m; i += 2) {
  120.             bufs2[icnt2] = G[i][jend];
  121.             icnt2 = icnt2 + 1;
  122.         }
  123.     }
  124.  
  125.     MPI_Isend(bufs1, icnt1, MPI_DOUBLE, iprev, 1, MPI_COMM_WORLD, &ireqs1);
  126.     MPI_Isend(bufs2, icnt2, MPI_DOUBLE, inext, 1, MPI_COMM_WORLD, &ireqs2);
  127.  
  128.     MPI_Irecv(bufr1, N, MPI_DOUBLE, iprev, 1, MPI_COMM_WORLD, &ireqr1);
  129.     MPI_Irecv(bufr2, N, MPI_DOUBLE, inext, 1, MPI_COMM_WORLD, &ireqr2);
  130.  
  131.     MPI_Wait(&ireqs1, &istatus);
  132.     MPI_Wait(&ireqs2, &istatus);
  133.     MPI_Wait(&ireqr1, &istatus);
  134.     MPI_Wait(&ireqr2, &istatus);
  135.  
  136.     if (mynode != 0) {
  137.         icnt = 0;
  138.         for (i = ir1; i <= m; i += 2) {
  139.             G[i][jsta - 1] = bufr1[icnt];
  140.             icnt = icnt + 1;
  141.         }
  142.     }
  143.  
  144.     if (mynode != (totalnodes - 1)) {
  145.         icnt = 0;
  146.         for (i = ir2; i <= m; i += 2) {
  147.             G[i][jend + 1] = bufr2[icnt];
  148.             icnt = icnt + 1;
  149.         }
  150.     }
  151.  
  152. }
  153.  
  154. int main (int argc, char *argv[])
  155. {
  156.         int ncol,nrow;             /* number of rows and columns */
  157.         double Gnew,r,omega,       /* some float values */
  158.            stopdiff,maxdiff,diff;  /* differences btw grid points in iters */
  159.         int i,j,phase,iteration;   /* counters */
  160.         int print = 0;
  161.         double inittime, totaltime;
  162.  
  163.         /* Initializing MPI */
  164.     MPI_Init(&argc, &argv);
  165.     MPI_Comm_size(MPI_COMM_WORLD, &totalnodes);
  166.     MPI_Comm_rank(MPI_COMM_WORLD, &mynode);
  167.  
  168.         /* set up problem size */
  169.         N = 1000;
  170.  
  171.         for(i=1; i<argc; i++) {
  172.                 if(!strcmp(argv[i], "-print")) {
  173.                         print = 1;
  174.                 } else {
  175.                         N = atoi(argv[i]);
  176.                 }
  177.         }
  178.  
  179.         if(mynode == 0) {
  180.                 fprintf(stderr, "Running %d x %d SOR\n", N, N);
  181.         }
  182.  
  183.         N += 2; /* add the two border lines */
  184.         /* finally N*N (from argv) array points will be computed */
  185.  
  186.         /* set up a quadratic grid */
  187.         ncol = nrow = N;
  188.         r        = 0.5 * ( cos( M_PI / ncol ) + cos( M_PI / nrow ) );
  189.         omega    = 2.0 / ( 1 + sqrt( 1 - r * r ) );
  190.         stopdiff = TOLERANCE / ( 2.0 - omega );
  191.         omega   *= MAGIC;
  192.  
  193.         alloc_grid(&G, N);
  194.         alloc_grid(&buff, N);
  195.         init_grid(G, N);
  196.  
  197.         // send receive buffers for left neighbor processes (iprev)
  198.     bufs1 = (double *) malloc(N * sizeof (double));
  199.     bufr1 = (double *) malloc(N * sizeof (double));
  200.  
  201.     // send receive buffers for right neighbor processes (inext)
  202.     bufs2 = (double *) malloc(N * sizeof (double));
  203.     bufr2 = (double *) malloc(N * sizeof (double));
  204.  
  205.     // start - end matrix array cols for each processes
  206.     range(1, N - 1, totalnodes, mynode, &jsta, &jend);
  207.  
  208.     inext = mynode + 1;
  209.     iprev = mynode - 1;
  210.  
  211.     if (inext == totalnodes)
  212.         inext = MPI_PROC_NULL;
  213.  
  214.     if (iprev == -1)
  215.         iprev = MPI_PROC_NULL;
  216.  
  217.         inittime = MPI_Wtime();
  218.  
  219.         /* now do the "real" computation */
  220.         iteration = 0;
  221.         do {
  222.                 maxdiff = 0.0;
  223.                 for ( phase = 0; phase < 2 ; phase++){
  224.  
  225.                         exchange(phase);
  226.  
  227.                         for ( i = 1 ; i < N-1 ; i++ ){
  228.                                 for ( j = jsta + (even(i) ^ phase); j <= jend ; j += 2 ){
  229.                                         Gnew = stencil(G,i,j);
  230.                                         diff = fabs(Gnew - G[i][j]);
  231.                                         if ( diff > maxdiff )
  232.                                                 maxdiff = diff;
  233.                                         G[i][j] = G[i][j] + omega * (Gnew-G[i][j]);
  234.                                 }
  235.                         }
  236.                 }
  237.  
  238.                 MPI_Allreduce(&maxdiff, &diff, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
  239.                 maxdiff = diff;
  240.  
  241.                 iteration++;
  242.         } while (maxdiff > stopdiff);
  243.  
  244.         totaltime = MPI_Wtime() - inittime;
  245.  
  246.         MPI_Barrier(MPI_COMM_WORLD);
  247.  
  248.         if(print == 1) {
  249.                 printf("Node: %d\n", mynode);
  250.                 print_grid(G, N);
  251.                  printf("\n");
  252.         }
  253.  
  254.         MPI_Finalize();
  255.  
  256.         return 0;
  257. }