examples/smpi/mvmul.c

   1 /* Copyright (c) 2009-2010, 2013-2014. The SimGrid Team.
   2  * All rights reserved.                                                     */
   3
   4 /* This program is free software; you can redistribute it and/or modify it
   5  * under the terms of the license (GNU LGPL) which comes with this package. */
   6
   7 #include <stdio.h>
   8 #include <stdlib.h>
   9 #include <ctype.h>
  10 #include <mpi.h>
  11 #include <xbt/str.h>
  12
  13 #define ITERATIONS         10
  14 #define USAGE_ERROR        1
  15 #define SANITY_ERROR       2
  16 #define GETTIMEOFDAY_ERROR 3
  17
  18 int main(int argc, char *argv[])
  19 {
  20   int size, rank;
  21   int N, n, i, j, k, current_iteration, successful_iterations = 0;
  22   double *matrix = NULL, *vector = NULL, *vcalc, *vcheck;
  23   MPI_Status status;
  24   struct timeval *start_time = NULL, *stop_time = NULL;
  25   long parallel_usecs, parallel_usecs_total = 0, sequential_usecs, sequential_usecs_total = 0;
  26
  27   MPI_Init(&argc, &argv);
  28
  29   MPI_Comm_rank(MPI_COMM_WORLD, &rank);
  30   MPI_Comm_size(MPI_COMM_WORLD, &size);
  31
  32   if (0 == rank) {
  33     // root node parses cmdline args
  34     if (2 > argc || !isdigit(*argv[1])) {
  35       printf("usage:\n%s <size>\n", argv[0]);
  36       MPI_Abort(MPI_COMM_WORLD, USAGE_ERROR);
  37       exit(USAGE_ERROR);
  38     }
  39
  40     N = xbt_str_parse_int(argv[1], "Invalid size: %s");
  41
  42     start_time = (struct timeval *) malloc(sizeof(struct timeval));
  43     stop_time = (struct timeval *) malloc(sizeof(struct timeval));
  44   }
  45
  46   for (current_iteration = 0; current_iteration < ITERATIONS; current_iteration++) {
  47     if (0 == rank) {
  48       matrix = (double *) malloc(N * N * sizeof(double));
  49       vector = (double *) malloc(N * sizeof(double));
  50
  51       for (i = 0; i < N * N; i++) {
  52         matrix[i] = (double) rand() / ((double) RAND_MAX + 1);
  53       }
  54
  55       for (i = 0; i < N; i++) {
  56         vector[i] = (double) rand() / ((double) RAND_MAX + 1);
  57       }
  58
  59       // for the sake of argument, the parallel algorithm begins when the root node begins to transmit the matrix to the
  60       // workers.
  61       if (-1 == gettimeofday(start_time, NULL)) {
  62         printf("couldn't set start_time on node 0!\n");
  63         MPI_Abort(MPI_COMM_WORLD, GETTIMEOFDAY_ERROR);
  64         exit(GETTIMEOFDAY_ERROR);
  65       }
  66
  67       for (i = 1; i < size; i++) {
  68         MPI_Send(&N, 1, MPI_INT, i, 0, MPI_COMM_WORLD);
  69       }
  70     } else {
  71       MPI_Recv(&N, 1, MPI_INT, 0, 0, MPI_COMM_WORLD, &status);
  72     }
  73
  74     // this algorithm uses at most N processors...
  75     if (rank < N) {
  76       if (size > N)
  77         size = N;
  78       n = N / size + ((rank < (N % size)) ? 1 : 0);
  79
  80       if (0 == rank) {
  81         for (i = 1, j = n; i < size && j < N; i++, j += k) {
  82           k = N / size + ((i < (N % size)) ? 1 : 0);
  83           MPI_Send(matrix + N * j, N * k, MPI_DOUBLE, i, 0, MPI_COMM_WORLD);
  84           MPI_Send(vector, N, MPI_DOUBLE, i, 0, MPI_COMM_WORLD);
  85         }
  86
  87         // sanity check
  88 #ifdef DEBUG
  89         if (i != size || j != N) {
  90           printf("index calc error: i = %d, size = %d, j = %d, N = %d\n", i, size, j, N);
  91           MPI_Abort(MPI_COMM_WORLD, SANITY_ERROR);
  92           exit(SANITY_ERROR);
  93         }
  94 #endif
  95
  96         vcalc = (double *) malloc(N * sizeof(double));
  97       } else {
  98         matrix = (double *) malloc(N * n * sizeof(double));
  99         vector = (double *) malloc(N * sizeof(double));
 100         vcalc = (double *) malloc(n * sizeof(double));
 101
 102         MPI_Recv(matrix, N * n, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD, &status);
 103         MPI_Recv(vector, N, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD, &status);
 104       }
 105
 106       for (i = 0; i < n; i++) {
 107         vcalc[i] = 0.0;
 108         for (j = 0; j < N; j++) {
 109           vcalc[i] += matrix[N * i + j] * vector[j];
 110         }
 111       }
 112
 113       if (0 != rank) {
 114         MPI_Send(vcalc, n, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD);
 115       } else {
 116         for (i = 1, j = n; i < size && j < N; i++, j += k) {
 117           k = N / size + ((i < (N % size)) ? 1 : 0);
 118           MPI_Recv(vcalc + j, k, MPI_DOUBLE, i, 0, MPI_COMM_WORLD, &status);
 119         }
 120
 121         // sanity check
 122 #ifdef DEBUG
 123         if (i != size || j != N) {
 124           printf("index calc error 2: i = %d, size = %d, j = %d, N = %d\n", i, size, j, N);
 125           MPI_Abort(MPI_COMM_WORLD, SANITY_ERROR);
 126           exit(SANITY_ERROR);
 127         }
 128 #endif
 129
 130         if (-1 == gettimeofday(stop_time, NULL)) {
 131           printf("couldn't set stop_time on node 0!\n");
 132           MPI_Abort(MPI_COMM_WORLD, GETTIMEOFDAY_ERROR);
 133           exit(GETTIMEOFDAY_ERROR);
 134         }
 135
 136         parallel_usecs = (stop_time->tv_sec * 1000000 + stop_time->tv_usec) -
 137                          (start_time->tv_sec * 1000000 + start_time->tv_usec);
 138
 139         if (-1 == gettimeofday(start_time, NULL)) {
 140           printf("couldn't set start_time on node 0!\n");
 141           MPI_Abort(MPI_COMM_WORLD, GETTIMEOFDAY_ERROR);
 142           exit(GETTIMEOFDAY_ERROR);
 143         }
 144         // calculate serially
 145         vcheck = (double *) malloc(N * sizeof(double));
 146         for (i = 0; i < N; i++) {
 147           vcheck[i] = 0.0;
 148           for (j = 0; j < N; j++) {
 149             vcheck[i] += matrix[N * i + j] * vector[j];
 150           }
 151         }
 152
 153         if (-1 == gettimeofday(stop_time, NULL)) {
 154           printf("couldn't set stop_time on node 0!\n");
 155           MPI_Abort(MPI_COMM_WORLD, GETTIMEOFDAY_ERROR);
 156           exit(GETTIMEOFDAY_ERROR);
 157         }
 158
 159         sequential_usecs = (stop_time->tv_sec * 1000000 + stop_time->tv_usec) -
 160                            (start_time->tv_sec * 1000000 + start_time->tv_usec);
 161
 162         // verify correctness
 163         for (i = 0; i < N && vcalc[i] == vcheck[i]; i++);
 164
 165         printf("prog: blocking, i: %d ", current_iteration);
 166
 167         if (i == N) {
 168           printf("ptime: %ld us, stime: %ld us, speedup: %.3f, nodes: %d, efficiency: %.3f\n",
 169                parallel_usecs, sequential_usecs, (double) sequential_usecs / (double) parallel_usecs, size,
 170                (double) sequential_usecs / ((double) parallel_usecs * (double) size));
 171
 172           parallel_usecs_total += parallel_usecs;
 173           sequential_usecs_total += sequential_usecs;
 174           successful_iterations++;
 175         } else {
 176           printf("parallel calc != serial calc, ");
 177         }
 178
 179         free(vcheck);
 180       }
 181
 182       free(matrix);
 183       free(vector);
 184       free(vcalc);
 185     }
 186   }
 187
 188   if (0 == rank) {
 189     printf("prog: blocking, ");
 190     if (0 < successful_iterations) {
 191       printf("iterations: %d, avg. ptime: %.3f us, avg. stime: %.3f us, avg. speedup: %.3f, nodes: %d, avg. efficiency: %.3f\n",
 192            successful_iterations, (double) parallel_usecs_total / (double) successful_iterations,
 193            (double) sequential_usecs_total / (double) successful_iterations,
 194            (double) sequential_usecs_total / (double) parallel_usecs_total, size,
 195            (double) sequential_usecs_total / ((double) parallel_usecs_total * (double) size));
 196     } else {
 197       printf("no successful iterations!\n");
 198     }
 199
 200     free(start_time);
 201     free(stop_time);
 202   }
 203
 204   MPI_Finalize();
 205
 206   return 0;
 207 }