src/smpi/sample/matrix.c

   1 /*
   2  * Mark Stillwell
   3  * ICS691: High Performance Computing
   4  * Fall 2006
   5  * Homework 3, Exercise 2, Step 1
   6  */
   7 #include <stdio.h>
   8 #include <stdlib.h>
   9 #include <libgen.h>
  10 #include <mpi.h>
  11
  12 #define ITERATIONS         10
  13 #define STEPS              1
  14 #define STEP_SIZE          0
  15
  16 #define USAGE_ERROR        1
  17 #define MALLOC_ERROR       2
  18 #define GETTIMEOFDAY_ERROR 3
  19
  20 void * checked_malloc(int rank, const char * varname, size_t size) {
  21   void * ptr;
  22   ptr = malloc(size);
  23   if (NULL == ptr) {
  24     printf("node %d could not malloc memory for %s.\n", rank, varname);
  25     MPI_Abort(MPI_COMM_WORLD, MALLOC_ERROR);
  26     exit(MALLOC_ERROR);
  27   }
  28   return ptr;
  29 }
  30
  31 int main(int argc, char* argv[]) {
  32
  33   // timing/system variables
  34   int iteration, iterations = ITERATIONS;
  35   int step, steps = STEPS, step_size = STEP_SIZE;
  36   long usecs, total_usecs;
  37   struct timeval *start_time, *stop_time;
  38   char *program;
  39
  40   // mpi/communications variables
  41   int rank;
  42   int row, col;
  43   MPI_Comm row_comm, col_comm;
  44
  45   // algorithm variables
  46   int N_start, N, P;
  47   int *A, *A_t, *B, *C, *D, *a, *b, *abuf, *bbuf;
  48   int n, i, j, k, I, J;
  49
  50   MPI_Init(&argc, &argv);
  51
  52   MPI_Comm_rank(MPI_COMM_WORLD, &rank);
  53
  54   if (0 == rank) {
  55     int size;
  56     MPI_Comm_size(MPI_COMM_WORLD, &size);
  57
  58     program = basename(argv[0]);
  59
  60     // root node parses cmdline args
  61     /*
  62     if (3 > argc || !isdigit(*argv[1]) || !isdigit(*argv[2])) {
  63       printf("usage:\n%s <N> <P> [<iterations>]\n", program);
  64       MPI_Abort(MPI_COMM_WORLD, USAGE_ERROR);
  65       exit(USAGE_ERROR);
  66     }
  67     */
  68
  69     //N_start = atoi(argv[1]);
  70     //P = atoi(argv[2]);
  71     N_start = 100;
  72     P = 2;
  73
  74     /*
  75     if (4 <= argc && isdigit(*argv[3])) {
  76       iterations = atoi(argv[3]);
  77     }
  78
  79     if (5 <= argc && isdigit(*argv[4])) {
  80       steps = atoi(argv[4]);
  81     }
  82
  83     if (6 <= argc && isdigit(*argv[5])) {
  84       step_size = atoi(argv[5]);
  85     }
  86     */
  87
  88     if (P*P != size) {
  89       printf("P^2 must equal size.\n");
  90       MPI_Abort(MPI_COMM_WORLD, USAGE_ERROR);
  91       exit(USAGE_ERROR);
  92     }
  93
  94     start_time = (struct timeval *)checked_malloc(rank, "start_time", sizeof(struct timeval));
  95     stop_time  = (struct timeval *)checked_malloc(rank, "stop_time",  sizeof(struct timeval));
  96
  97   }
  98
  99   // send command line parameters except N, since it can vary
 100   MPI_Bcast(&P, 1, MPI_INT, 0, MPI_COMM_WORLD);
 101   MPI_Bcast(&iterations, 1, MPI_INT, 0, MPI_COMM_WORLD);
 102   MPI_Bcast(&steps, 1, MPI_INT, 0, MPI_COMM_WORLD);
 103   MPI_Bcast(&step_size, 1, MPI_INT, 0, MPI_COMM_WORLD);
 104
 105   row = rank / P;
 106   col = rank % P;
 107
 108   // create row/column communicators
 109   MPI_Comm_split(MPI_COMM_WORLD, row, col, &row_comm);
 110   MPI_Comm_split(MPI_COMM_WORLD, col, row, &col_comm);
 111
 112   for (step = 0; step < steps; step++) {
 113
 114     total_usecs = 0;
 115
 116     if (0 == rank) {
 117       N = N_start + step * step_size;
 118       if ((N/P)*P != N) {
 119         printf("P must divide N and %d does not divide %d.\n", P, N);
 120         N = -1;
 121       }
 122     }
 123
 124     MPI_Bcast(&N, 1, MPI_INT, 0, MPI_COMM_WORLD);
 125
 126     // if root passes N = -1, skip this round
 127     if (-1 == N) continue;
 128
 129     n = N / P;
 130
 131     // initialize matrix components
 132     A   = (int *)checked_malloc(rank, "A",   n*n*sizeof(int));
 133     A_t = (int *)checked_malloc(rank, "A_t", n*n*sizeof(int));
 134     B   = (int *)checked_malloc(rank, "B",   n*n*sizeof(int));
 135     C   = (int *)checked_malloc(rank, "C",   n*n*sizeof(int));
 136     D   = (int *)checked_malloc(rank, "D",   n*n*sizeof(int));
 137
 138     for (i = 0; i < n; i++) {
 139       for (j = 0; j < n; j++) {
 140
 141         I = n*row+i;
 142         J = n*col+j;
 143
 144         A[n*i+j] = I+J;
 145         B[n*i+j] = I;
 146
 147         // d is the check matrix
 148         D[n*i+j] = 0;
 149         for (k = 0; k < N; k++) {
 150           // A[I,k] = I+k
 151           // B[k,J] = k
 152           D[n*i+j] += (I+k) * k;
 153         }
 154
 155       }
 156     }
 157
 158     // buffers
 159     abuf = (int *)checked_malloc(rank, "abuf", n*sizeof(int));
 160     bbuf = (int *)checked_malloc(rank, "bbuf", n*sizeof(int));
 161
 162     for (iteration = 0; iteration < iterations; iteration++) {
 163
 164       for (i = 0; i < n*n; i++) {
 165         C[i] = 0;
 166       }
 167
 168       // node zero sets start time
 169       if (0 == rank && -1 == gettimeofday(start_time, NULL)) {
 170         printf("couldn't set start_time on node 0!\n");
 171         MPI_Abort(MPI_COMM_WORLD, GETTIMEOFDAY_ERROR);
 172         exit(GETTIMEOFDAY_ERROR);
 173       }
 174
 175       // populate transpose of A
 176       for (i = 0; i < n; i++) {
 177         for (j = 0; j < n; j++) {
 178           A_t[n*i+j] = A[n*j+i];
 179         }
 180       }
 181
 182       // perform calculations
 183       for (k = 0; k < N; k++) {
 184
 185         if (k/n == col) {
 186           a = A_t + n*(k%n);
 187         } else {
 188           a = abuf;
 189         }
 190
 191         if (k/n == row) {
 192           b = B + n*(k%n);
 193         } else {
 194           b = bbuf;
 195         }
 196
 197         MPI_Bcast(a, n, MPI_INT, k/n, row_comm);
 198         MPI_Bcast(b, n, MPI_INT, k/n, col_comm);
 199
 200         for (i = 0; i < n; i++) {
 201           for (j = 0; j < n; j++) {
 202             C[n*i+j] += a[i] * b[j];
 203           }
 204         }
 205
 206       } // for k
 207
 208       // node zero sets stop time
 209       if (0 == rank && -1 == gettimeofday(stop_time, NULL)) {
 210         printf("couldn't set stop_time on node 0!\n");
 211         MPI_Abort(MPI_COMM_WORLD, GETTIMEOFDAY_ERROR);
 212         exit(GETTIMEOFDAY_ERROR);
 213       }
 214
 215       // check calculation
 216       for (i = 0; i < n*n && C[i] == D[i]; i++);
 217       j = (n*n == i);
 218       MPI_Reduce(&j, &k, 1, MPI_INT, MPI_LAND, 0, MPI_COMM_WORLD);
 219
 220       // node zero prints stats
 221       if (0 == rank) {
 222         usecs = (stop_time->tv_sec*1000000+stop_time->tv_usec) - (start_time->tv_sec*1000000+start_time->tv_usec);
 223         printf("prog: %s, N: %d, P: %d, procs: %d, time: %d us, check: %d\n", program, N, P, P*P, usecs, k);
 224         total_usecs += usecs;
 225       }
 226
 227     }
 228
 229     // node 0 prints final stats
 230     if (0 == rank) {
 231       printf("prog: %s, N: %d, P: %d, procs: %d, iterations: %d, avg. time: %d us\n",
 232           program, N, P, P*P, iterations, total_usecs / iterations);
 233     }
 234
 235     // free data structures
 236     free(A);
 237     free(A_t);
 238     free(B);
 239     free(C);
 240     free(D);
 241     free(abuf);
 242     free(bbuf);
 243
 244   }
 245
 246   if (0 == rank) {
 247     free(start_time);
 248     free(stop_time);
 249   }
 250
 251   MPI_Finalize();
 252
 253   return 0;
 254 }