red/black sor

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <sys/time.h>

#define TOLERANCE 0.00002        /* termination criterion */
#define MAGIC 0.8                /* magic factor */

int even (int i)
{
    return !(i & 1);
}


double stencil (double** G, int row, int col)
{
    return (G[row-1][col] + G[row+1][col] +
        G[row][col-1] + G[row][col+1] ) / 4.0;
}


void alloc_grid(double ***Gptr, int N)
{
    int i;
    double** G = (double**)malloc(N*sizeof(double*));
    if ( G == 0 ) {
        fprintf(stderr, "malloc failed\n");
        exit(42);
    }

    for (i = 0; i<N; i++) { /* malloc the own range plus one more line */
        /* of overlap on each side */
        G[i] = (double*)malloc(N*sizeof(double));
        if ( G[i] == 0 ) {
            fprintf(stderr, "malloc failed\n");
            exit(42);
        }
    }

    *Gptr = G;
}


void init_grid(double **G, int N)
{
    int i, j;

    /* initialize the grid */
    for (i = 0; i < N; i++){
        for (j = 0; j < N; j++){
            if (i == 0) G[i][j] = 4.56;
            else if (i == N-1) G[i][j] = 9.85;
            else if (j == 0) G[i][j] = 7.32;
            else if (j == N-1) G[i][j] = 6.88;
            else G[i][j] = 0.0;
        }
    }
}


void print_grid(double** G, int N)
{
    int i, j;

    for ( i = 1 ; i < N-1 ; i++ ) {
        for ( j = 1 ; j < N-1 ; j++ ) {
            printf("%10.3f ", G[i][j]);
        }
        printf("\n");
    }
}


int main (int argc, char *argv[])
{
    int N;                     /* problem size */
    int ncol,nrow;             /* number of rows and columns */
    double Gnew,r,omega,       /* some float values */
       stopdiff,maxdiff,diff;  /* differences btw grid points in iters */
    double **G;                /* the grid */
    int i,j,phase,iteration;   /* counters */
    struct timeval start;
    struct timeval end;
    double time;
    int print = 0;

    /* set up problem size */
    N = 1000;

    for(i=1; i<argc; i++) {
        if(!strcmp(argv[i], "-print")) {
            print = 1;
        } else {
            N = atoi(argv[i]);
        }
    }

    fprintf(stderr, "Running %d x %d SOR\n", N, N);

    N += 2; /* add the two border lines */
    /* finally N*N (from argv) array points will be computed */

    /* set up a quadratic grid */
    ncol = nrow = N;
    r        = 0.5 * ( cos( M_PI / ncol ) + cos( M_PI / nrow ) );
    omega    = 2.0 / ( 1 + sqrt( 1 - r * r ) );
    stopdiff = TOLERANCE / ( 2.0 - omega );
    omega   *= MAGIC;

    alloc_grid(&G, N);
    init_grid(G, N);

    if(gettimeofday(&start, 0) != 0) {
        fprintf(stderr, "could not do timing\n");
        exit(1);
    }

    /* now  the "real" computation */
    iteration = 0;
    do {
        maxdiff = 0.0;
        for ( phase = 0; phase < 2 ; phase++){
            for ( i = 1 ; i < N-1 ; i++ ){
                for ( j = 1 + (even(i) ^ phase); j < N-1 ; j += 2 ){
                    Gnew = stencil(G,i,j);
                    diff = fabs(Gnew - G[i][j]);
                    if ( diff > maxdiff )
                        maxdiff = diff;
                    G[i][j] = G[i][j] + omega * (Gnew-G[i][j]);
                }
            }
        }
        iteration++;
    } while (maxdiff > stopdiff);

    if(gettimeofday(&end, 0) != 0) {
        fprintf(stderr, "could not do timing\n");
        exit(1);
    }

    time = (end.tv_sec + (end.tv_usec / 1000000.0)) -
        (start.tv_sec + (start.tv_usec / 1000000.0));

    fprintf(stderr, "SOR took %10.3f seconds\n", time);

    printf("Used %5d iterations, diff is %10.6f, allowed diff is %10.6f\n",
           iteration,maxdiff,stopdiff);

    if(print == 1) {
        print_grid(G, N);
    }

    return 0;
}