プログラム

このプログラムではノードが１つだと何もできない。 rankが0のものは制御用、残りが計算用となっている。 m = (N + myid) / (np -1)とすることで、うまく各ノードにデータを配布している。

#include <stdio.h>
#include <malloc.h>
#include <math.h>
#include "mpi.h"

#define N 1000

int main(int argc, char **argv){
  int myid, np;
  double *x, *y, *z;
  double *parx, *pary, *parz;          // each PEs data
  double wtime, timer;
  MPI_Status st;

  int i, j, k, m, offset;

  MPI_Init(&argc, &argv);
  MPI_Comm_rank(MPI_COMM_WORLD, &myid);
  MPI_Comm_size(MPI_COMM_WORLD, &np);

  /* Data Initialize */
  if(myid == 0){
    x = (double *)malloc(sizeof(double) * N);
    y = (double *)malloc(sizeof(double) * N);
    z = (double *)malloc(sizeof(double) * N);

    for(i = 0; i < N; i++)
      x[i] = y[i] = 0.00001 * i;
  }

  /* timer start */
  MPI_Barrier(MPI_COMM_WORLD);
  if(myid == 0){
  wtime = MPI_Wtime();
  }
  /* Calc datasize of each PEs by itself */
  m = (N + myid) / (np -1);
  parx = (double *)malloc(sizeof(double) * m);
  pary = (double *)malloc(sizeof(double) * m);
  parz = (double *)malloc(sizeof(double) * m);

  /* send data from PE0 to each PEs */
  if(myid == 0){
    /* selfcopy PE0 to PE0 */

    memcpy(parx, x, sizeof(double) * m);
    memcpy(pary, y, sizeof(double) * m);


    offset = 0;
    for(j = 1; j < np; j++){
      /* PE0 calcs other PEs datasizes */
      k = (N + j) / (np - 1);
      MPI_Send(&x[offset], k, MPI_DOUBLE, j, 1000, MPI_COMM_WORLD);
      MPI_Send(&y[offset], k, MPI_DOUBLE, j, 2000, MPI_COMM_WORLD);
      offset = offset + k;
    }
    /*
    for(i = 0; i < m; i++){
      z[i] = sqrt(parx[i] * parx[i] + pary[i] * pary[i]);
    }
    */
    offset = 0;
    for(j = 1; j < np; j++){
      k = (N + j) / (np - 1);
      MPI_Recv(&z[offset], k, MPI_DOUBLE, j, 3000, MPI_COMM_WORLD, &st);
      offset = offset + k;
    }

  }
  else{
    /* Each PE knows itself datasize at m */
    MPI_Recv(parx, m, MPI_DOUBLE, 0, 1000, MPI_COMM_WORLD, &st);
    MPI_Recv(pary, m, MPI_DOUBLE, 0, 2000, MPI_COMM_WORLD, &st);

    /* calculate norm */
    for(i = 0;i < m;i++){
      parz[i] = sqrt(parx[i] * parx[i] + pary[i] * pary[i]);
    }

    MPI_Send(parz, m, MPI_DOUBLE, 0, 3000, MPI_COMM_WORLD);
  }



  /* for debug */
  /*
  for(i = 0; i < np; i++){
    if(myid == i){
      for(k = 0; k < m; k++)
        printf("%d: (x, y) = (%10.9f, %10.9f)\n", myid, parx[k], pary[k]);
      fflush(stdout);
    }
  }
  */
  if(myid == 0){
    timer = MPI_Wtime() - wtime;
    /*
    for(i=0;i<N;i++){
      printf("norm:%d = %f\n",i,z[i]);
      fflush(stdout);
      }
    */
    printf("time:%f\n",timer);
  }


  MPI_Finalize();

}