2011年8月3日水曜日

メモ:C=A*B MPI

N*Nの行列計算C = A・Bを行うためのプログラムにMPIを組み込んだもの。

これは最終的に出来上がったプログラムだが、一つ前のステップとして、各MPIプロセスはA,Bの行列データの全てを持ち(N*Nのデータ)、それぞれのプロセスでそのデータの一部を使い、Cの別々の要素について計算を行うというものを作った。このプログラムでは、各プロセスの持つデータも分割するようにしている。勿論この方が無駄なデータを保持する必要が無く最適化されたプログラムである。

また、並列処理以外にも、行列積のプログラムではループ交換法やタイリング(ブロック化)法により、アルゴリズムを最適化する必要がある。ループ交換法では行列の積を内積形式にするか、外積形式にするか、中間積形式にするかで、言語により性能が違ってくる。このプログラムでは内積形式を使っているが、行方向のアクセスはC言語では連続メモリアクセスとなり、Fortranでは不連続アクセスとなる(列方向は反対)。タイリングの仕方も、一次元なのか、二次元なのか、サイクリックなのか否かで性能が変わる(データサイズや計算負荷のバランス、プロセス数によっても変わる)。

このプログラムでは各プロセスはAの[N/プロセス数]*[N]、Bの[N]*[N/プロセス数]のデータを持つので、各プロセスでCの[N/プロセス数]*[N]を計算するためには、通信を使ってBのデータを更新する必要がある。このプログラムではBのデータの更新のため、循環左シフト転送を実装している。ただこのアルゴリズムでMPI実装することを考えると、全てのプロセスが一斉にメッセージの送信と受信(MPI_Send,MPI_Recv)をすると送信バッファが永遠にクリアできず、デットロックが起きてしまう。そのため、プロセス番号を2で割る条件で、あるプロセスは最初に送信し受信、一方のプロセスは受信し送信するという風に通信をずらしている(タグ番号もずらすことが大事)。

ローカルな行列積では、Cの列方向の要素の場所をプロセスごとにずれるようにしなければ最終的な計算結果(N*Nの行列C)を得られないため、jstart変数を定義している。

--------------C=A*B-------------

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include <mpi.h>

#define  N        128
#define  NPROCS   64

#define  DEBUG  1
#define  EPS    1.0e-18

double  A[N/NPROCS][N];
double  B[N][N/NPROCS];
double  C[N/NPROCS][N];

double  B_T[N][N/NPROCS];

int     myid, numprocs;

void MyMatMat(double [][], double [][], double [][], int); 

void main(int argc, char* argv[]) {

     double  t0, t1, t2, t_w;
     double  dc_inv, d_mflops;

     int     ierr;
     int     i, j;      
     int     iflag, iflag_t;

     ierr = MPI_Init(&argc, &argv);
     ierr = MPI_Comm_rank(MPI_COMM_WORLD, &myid);
     ierr = MPI_Comm_size(MPI_COMM_WORLD, &numprocs);

     /* matrix generation --------------------------*/
     if (DEBUG == 1) {
       for(j=0; j<N/NPROCS; j++) {
         for(i=0; i<N; i++) {
           A[j][i] = 1.0;
           C[j][i] = 0.0;
         }
       }
       for(j=0; j<N; j++) {
         for(i=0; i<N/NPROCS; i++) {
           B[j][i] = 1.0;
         }
       }

     } else {
       srand(myid);
       dc_inv = 1.0/(double)RAND_MAX;
 
      for(j=0; j<N/NPROCS; j++) {
         for(i=0; i<N; i++) {
           A[j][i] = rand()*dc_inv;
           C[j][i] = 0.0;
         }
       }
      for(j=0; j<N; j++) {
         for(i=0; i<N/NPROCS; i++) {
           B[j][i] = rand()*dc_inv;
         }
       }

     } /* end of matrix generation --------------------------*/

     /* Start of mat-vec routine ----------------------------*/
     ierr = MPI_Barrier(MPI_COMM_WORLD);
     t1 = MPI_Wtime();

     MyMatMat(C, A, B, N);

     ierr = MPI_Barrier(MPI_COMM_WORLD);
     t2 = MPI_Wtime();
     t0 =  t2 - t1; 
     ierr = MPI_Reduce(&t0, &t_w, 1, MPI_DOUBLE, MPI_MAX, 0, MPI_COMM_WORLD);
     /* End of mat-vec routine --------------------------- */

     if (myid == 0) {

       printf("N  = %d \n",N);
       printf("Mat-Mat time  = %lf [sec.] \n",t_w);

       d_mflops = 2.0*(double)N*(double)N*(double)N/t_w;
       d_mflops = d_mflops * 1.0e-6;
       printf(" %lf [MFLOPS] \n", d_mflops);
     }

     if (DEBUG == 1) {
       /* Verification routine ----------------- */
       iflag = 0;
       for(j=0; j<N/NPROCS; j++) { 
         for(i=0; i<N; i++) { 
           if (fabs(C[j][i] - (double)N) > EPS) {
             printf(" Error! in ( %d , %d )-th argument in PE %d \n",j, i, myid);
             iflag = 1;
             exit(1);
           } 
         }
       }
       /* ------------------------------------- */

       MPI_Reduce(&iflag, &iflag_t, 1, MPI_INT, MPI_SUM, 0, MPI_COMM_WORLD);
       if (myid == 0) {
         if (iflag_t == 0) printf(" OK! \n");
       }

     }

     ierr = MPI_Finalize();

     exit(0);
}

void MyMatMat(double C[N/NPROCS][N], 
              double A[N/NPROCS][N], double B[N][N/NPROCS], int n) 
{
     int  i, j, k, p, q;
     int  ib;
     int  ierr;
     int  iloop; 
     int  jstart; 
     int  isendPE, irecvPE;
   
     MPI_Status istatus;

     /* Information of Send and recv PEs */
     isendPE = myid - 1;
     irecvPE = myid + 1;
     if (myid == 0) isendPE = numprocs - 1;
     if (myid == numprocs - 1) irecvPE = 0;   
    
     /* Block Length */
     ib = n/numprocs;

     for(iloop=0;iloop<=NPROCS-1;iloop++){
        jstart = ib*((myid+iloop)%NPROCS);
        for(i=0; i<ib; i++) {
          for(j=0; j<ib; j++) {
            for(k=0; k<n; k++) {
              C[i][jstart+j] += A[i][k] * B[k][j];
            }
          }
        }

        if (iloop != (numprocs-1)) {
           if(myid % 2 == 0){
              ierr = MPI_Send(B,ib*n,MPI_DOUBLE,isendPE,iloop,MPI_COMM_WORLD);
              ierr = MPI_Recv(B_T,ib*n,MPI_DOUBLE,irecvPE,iloop+numprocs,MPI_COMM_WORLD,&istatus);
           }
           else {
              ierr = MPI_Recv(B_T,ib*n,MPI_DOUBLE,irecvPE,iloop,MPI_COMM_WORLD,&istatus);
              ierr = MPI_Send(B,ib*n,MPI_DOUBLE,isendPE,iloop+numprocs,MPI_COMM_WORLD);
           }
           for(p=0;p<=ib-1;p++){
             for(q=0;q<=n-1;q++){
                B[p][q] = B_T[p][q];
             }
           }
        }

     }
}

0 件のコメント: