Don't like ads? PRO users don't see any ads ;-)
Public Pastes
Guest

red/black sor

By: a guest on Jan 24th, 2012  |  syntax: C  |  size: 3.29 KB  |  hits: 138  |  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.  
  7. #define TOLERANCE 0.00002        /* termination criterion */
  8. #define MAGIC 0.8                /* magic factor */
  9.  
  10. int even (int i)
  11. {
  12.         return !(i & 1);
  13. }
  14.  
  15.  
  16. double stencil (double** G, int row, int col)
  17. {
  18.         return (G[row-1][col] + G[row+1][col] +
  19.                 G[row][col-1] + G[row][col+1] ) / 4.0;
  20. }
  21.  
  22.  
  23. void alloc_grid(double ***Gptr, int N)
  24. {
  25.         int i;
  26.         double** G = (double**)malloc(N*sizeof(double*));
  27.         if ( G == 0 ) {
  28.                 fprintf(stderr, "malloc failed\n");
  29.                 exit(42);
  30.         }
  31.  
  32.         for (i = 0; i<N; i++) { /* malloc the own range plus one more line */
  33.                 /* of overlap on each side */
  34.                 G[i] = (double*)malloc(N*sizeof(double));
  35.                 if ( G[i] == 0 ) {
  36.                         fprintf(stderr, "malloc failed\n");
  37.                         exit(42);
  38.                 }
  39.         }
  40.  
  41.         *Gptr = G;
  42. }
  43.  
  44.  
  45. void init_grid(double **G, int N)
  46. {
  47.         int i, j;
  48.  
  49.         /* initialize the grid */
  50.         for (i = 0; i < N; i++){
  51.                 for (j = 0; j < N; j++){
  52.                         if (i == 0) G[i][j] = 4.56;
  53.                         else if (i == N-1) G[i][j] = 9.85;
  54.                         else if (j == 0) G[i][j] = 7.32;
  55.                         else if (j == N-1) G[i][j] = 6.88;
  56.                         else G[i][j] = 0.0;
  57.                 }
  58.         }
  59. }
  60.  
  61.  
  62. void print_grid(double** G, int N)
  63. {
  64.         int i, j;
  65.  
  66.         for ( i = 1 ; i < N-1 ; i++ ) {
  67.                 for ( j = 1 ; j < N-1 ; j++ ) {
  68.                         printf("%10.3f ", G[i][j]);
  69.                 }
  70.                 printf("\n");
  71.         }
  72. }
  73.  
  74.  
  75. int main (int argc, char *argv[])
  76. {
  77.         int N;                     /* problem size */
  78.         int ncol,nrow;             /* number of rows and columns */
  79.         double Gnew,r,omega,       /* some float values */
  80.            stopdiff,maxdiff,diff;  /* differences btw grid points in iters */
  81.         double **G;                /* the grid */
  82.         int i,j,phase,iteration;   /* counters */
  83.         struct timeval start;
  84.         struct timeval end;
  85.         double time;
  86.         int print = 0;
  87.  
  88.         /* set up problem size */
  89.         N = 1000;
  90.  
  91.         for(i=1; i<argc; i++) {
  92.                 if(!strcmp(argv[i], "-print")) {
  93.                         print = 1;
  94.                 } else {
  95.                         N = atoi(argv[i]);
  96.                 }
  97.         }
  98.  
  99.         fprintf(stderr, "Running %d x %d SOR\n", N, N);
  100.  
  101.         N += 2; /* add the two border lines */
  102.         /* finally N*N (from argv) array points will be computed */
  103.  
  104.         /* set up a quadratic grid */
  105.         ncol = nrow = N;
  106.         r        = 0.5 * ( cos( M_PI / ncol ) + cos( M_PI / nrow ) );
  107.         omega    = 2.0 / ( 1 + sqrt( 1 - r * r ) );
  108.         stopdiff = TOLERANCE / ( 2.0 - omega );
  109.         omega   *= MAGIC;
  110.  
  111.         alloc_grid(&G, N);
  112.         init_grid(G, N);
  113.  
  114.         if(gettimeofday(&start, 0) != 0) {
  115.                 fprintf(stderr, "could not do timing\n");
  116.                 exit(1);
  117.         }
  118.  
  119.         /* now  the "real" computation */
  120.         iteration = 0;
  121.         do {
  122.                 maxdiff = 0.0;
  123.                 for ( phase = 0; phase < 2 ; phase++){
  124.                         for ( i = 1 ; i < N-1 ; i++ ){
  125.                                 for ( j = 1 + (even(i) ^ phase); j < N-1 ; j += 2 ){
  126.                                         Gnew = stencil(G,i,j);
  127.                                         diff = fabs(Gnew - G[i][j]);
  128.                                         if ( diff > maxdiff )
  129.                                                 maxdiff = diff;
  130.                                         G[i][j] = G[i][j] + omega * (Gnew-G[i][j]);
  131.                                 }
  132.                         }
  133.                 }
  134.                 iteration++;
  135.         } while (maxdiff > stopdiff);
  136.  
  137.         if(gettimeofday(&end, 0) != 0) {
  138.                 fprintf(stderr, "could not do timing\n");
  139.                 exit(1);
  140.         }
  141.  
  142.         time = (end.tv_sec + (end.tv_usec / 1000000.0)) -
  143.                 (start.tv_sec + (start.tv_usec / 1000000.0));
  144.  
  145.         fprintf(stderr, "SOR took %10.3f seconds\n", time);
  146.  
  147.         printf("Used %5d iterations, diff is %10.6f, allowed diff is %10.6f\n",
  148.                iteration,maxdiff,stopdiff);
  149.  
  150.         if(print == 1) {
  151.                 print_grid(G, N);
  152.         }
  153.  
  154.         return 0;
  155. }
create a new version of this paste RAW Paste Data