MPI學習筆記（二）：矩陣相乘的兩種實現方法

2022 年 7 月 25 日
筆記
MPI

mpi矩陣乘法（C=αAB+βC）

最近領導讓把之前安裝的軟體lapack、blas里的dgemm運算提取出來獨立作為一套程式，然後把這段程式改為並行的，並測試一下進程規模擴展到128時的並行效率。
我發現這個是dgemm.f文件，裡面主要是對C=αAB+βC的實現，因此在此總結一下MPI的矩陣乘法使用。
其主要思想：是把相乘的矩陣按行分解（任務分解），分別分給不同的進程，然後在匯總到一個進程上，在程式上實現則用到了主從模式，人為的把進程分為主進程和從進程，主進程負責對原始矩陣初始化賦值，並把數據均勻分發（為了負載均衡）到從進程上進行相乘運算，主要用到的知識是MPI點對點通訊和組通訊的機制。

一、使用簡單的MPI_Send和MPI_Recv實現

#include <stdio.h>
#include "mpi.h"
#include <stdlib.h>
#include "functions.h"

#define M 1000 // 矩陣維度
#define N 1100
#define K 900

int main(int argc, char **argv)
{
   int my_rank,comm_sz,line;
   double start, stop; //計時時間
   MPI_Status status;
   char Processorname[20];

   double *Matrix_A,*Matrix_B,*Matrix_C,*ans,*buffer_A,*buffer_C;
   double alpha=2,beta=2; // 係數C=aA*B+bC

   MPI_Init(&argc,&argv);
   MPI_Comm_size(MPI_COMM_WORLD, &comm_sz);
   MPI_Comm_rank(MPI_COMM_WORLD,&my_rank);

   line=M/comm_sz; // 每個進程分多少行數據
   Matrix_A=(double*)malloc(M*N*sizeof(double));
   Matrix_B=(double*)malloc(N*K*sizeof(double));
   Matrix_C=(double*)malloc(M*K*sizeof(double));
   buffer_A=(double*)malloc(line*N*sizeof(double)); // A的均分行的數據
   buffer_C=(double*)malloc(line*K*sizeof(double)); // C的均分行的數據
   ans=(double*)malloc(line*K*sizeof(double)); // 臨時保存部分數據計算結果

   // 給矩陣A B,C賦值
   if(my_rank==0){
      start=MPI_Wtime();
      for(int i=0;i<M;i++){
         for(int j=0;j<N;j++)
            Matrix_A[i*N+j]=i+1;
      }
      for(int i=0;i<N;i++){
         for(int j=0;j<K;j++)
            Matrix_B[i*K+j]=j+1;
      }
      for(int i=0;i<M;i++){
         for(int j=0;j<K;j++)
            Matrix_C[i*K+j]=1;
      }

      // 輸出A,B,C
      /*Matrix_print(Matrix_A,M,N);
      Matrix_print(Matrix_B,N,K);
      Matrix_print(Matrix_C,M,K);
      */
      /*將矩陣廣播出去*/
      for(int i=1;i<comm_sz;i++){
         MPI_Send(Matrix_A+(i-1)*line*N,line*N,MPI_DOUBLE,i,66,MPI_COMM_WORLD);
         MPI_Send(Matrix_C+(i-1)*line*K,line*K,MPI_DOUBLE,i,99,MPI_COMM_WORLD);
      }
      MPI_Bcast(Matrix_B,N*K,MPI_DOUBLE,0,MPI_COMM_WORLD);

      // 接收從進程的計算結果
      for(int p=1;p<comm_sz;p++){
         MPI_Recv(ans,line*K,MPI_DOUBLE,p,33,MPI_COMM_WORLD,&status);
         for(int i=0;i<line;i+=comm_sz)
            for(int j=0;j<K;j++)
               Matrix_C[((p-1)*line+i)*K+j]=ans[i*K+j];
      }

      // 計算A剩下的行數據
      for(int i=(comm_sz-1)*line;i<M;i++){
         for(int j=0;j<K;j++){
            double temp=0;
            for(int p=0;p<N;p++)
               temp+=Matrix_A[i*N+p]*Matrix_B[p*K+j];
            Matrix_C[i*K+j]=alpha*temp+beta*Matrix_C[i*K+j];
         }
      }

      //Matrix_print(Matrix_C,M,K);
      stop=MPI_Wtime();

      printf("rank:%d time:%lfs\n",my_rank,stop-start);

      free(Matrix_A);
      free(Matrix_B);
      free(Matrix_C);
      free(buffer_A);
      free(buffer_C);
      free(ans);
   }
   else{
      //接收廣播的數據
      MPI_Recv(buffer_A,line*N,MPI_DOUBLE,0,66,MPI_COMM_WORLD,&status);
      MPI_Recv(buffer_C,line*K,MPI_DOUBLE,0,99,MPI_COMM_WORLD,&status);
      MPI_Bcast(Matrix_B,N*K,MPI_DOUBLE,0,MPI_COMM_WORLD);

      //計算乘積結果,並將結果發送給主進程
      for(int i=0;i<line;i++){
         for(int j=0;j<K;j++){
            double temp=0;
            for(int p=0;p<N;p++){
               temp+=buffer_A[i*N+p]*Matrix_B[p*K+j];
            }
            ans[i*line+j]=alpha*temp+beta*buffer_C[i*K+j];
         }
      }
      MPI_Send(ans,line*K,MPI_DOUBLE,0,33,MPI_COMM_WORLD);
   }

   MPI_Finalize();
   return 0;
}

二、使用較高級的MPI_Scatter和MPI_Gather實現

#include <stdio.h>
#include "mpi.h"
#include <stdlib.h>
#include "functions.h"

#define M 1200 // 矩陣維度
#define N 1000
#define K 1100

int main(int argc, char **argv)
{
   int my_rank,comm_sz,line;
   double start, stop; //計時時間
   MPI_Status status;

   double *Matrix_A,*Matrix_B,*Matrix_C,*ans,*buffer_A,*buffer_C,*result_Matrix;
   double alpha=2,beta=2; // 係數C=aA*B+bC

   MPI_Init(&argc,&argv);
   MPI_Comm_size(MPI_COMM_WORLD, &comm_sz);
   MPI_Comm_rank(MPI_COMM_WORLD,&my_rank);

   line=M/comm_sz; // 每個進程分多少行數據
   Matrix_A=(double*)malloc(M*N*sizeof(double));
   Matrix_B=(double*)malloc(N*K*sizeof(double));
   Matrix_C=(double*)malloc(M*K*sizeof(double));
   buffer_A=(double*)malloc(line*N*sizeof(double)); // A的均分行的數據
   buffer_C=(double*)malloc(line*K*sizeof(double)); // C的均分行的數據
   ans=(double*)malloc(line*K*sizeof(double)); // 保存部分數據計算結果
   result_Matrix=(double*)malloc(M*K*sizeof(double)); // 保存數據計算結果

   // 給矩陣A B,C賦值
   if(my_rank==0){
      start=MPI_Wtime();
      for(int i=0;i<M;i++){
         for(int j=0;j<N;j++)
            Matrix_A[i*N+j]=i+1;
         for(int p=0;p<K;p++)
            Matrix_C[i*K+p]=1;
      }
      for(int i=0;i<N;i++){
         for(int j=0;j<K;j++)
            Matrix_B[i*K+j]=j+1;
      }

      // 輸出A,B,C
      //Matrix_print(Matrix_A,M,N);
      //Matrix_print(Matrix_B,N,K);
      //Matrix_print(Matrix_C,M,K);
   }

   // 數據分發
   MPI_Scatter(Matrix_A,line*N,MPI_DOUBLE,buffer_A,line*N,MPI_DOUBLE,0,MPI_COMM_WORLD);
   MPI_Scatter(Matrix_C,line*K,MPI_DOUBLE,buffer_C,line*K,MPI_DOUBLE,0,MPI_COMM_WORLD);
   // 數據廣播
   MPI_Bcast(Matrix_B,N*K,MPI_DOUBLE,0,MPI_COMM_WORLD);

   // 計算 結果
   for(int i=0;i<line;i++){
      for(int j=0;j<K;j++){
         double temp=0;
         for(int p=0;p<N;p++)
            temp+=buffer_A[i*N+p]*Matrix_B[p*K+j];
         ans[i*K+j]=alpha*temp+beta*buffer_C[i*K+j];
      }
   }
   // 結果聚集
   MPI_Gather(ans,line*K,MPI_DOUBLE,result_Matrix,line*K,MPI_DOUBLE,0,MPI_COMM_WORLD);

   // 計算A剩下的行數據
   if(my_rank==0){
      int rest=M%comm_sz;
      if(rest!=0){
         for(int i=M-rest-1;i<M;i++)
            for(int j=0;j<K;j++){
               double temp=0;
               for(int p=0;p<N;p++)
                  temp+=Matrix_A[i*N+p]*Matrix_B[p*K+j];
               result_Matrix[i*K+j]=alpha*temp+beta*Matrix_C[i*K+j];
            }
      }

      //Matrix_print(result_Matrix,M,K);
      stop=MPI_Wtime();

      printf("rank:%d time:%lfs\n",my_rank,stop-start);
   }

   free(Matrix_A);
   free(Matrix_B);
   free(Matrix_C);
   free(ans);
   free(buffer_A);
   free(buffer_C);
   free(result_Matrix);

   MPI_Finalize();
   return 0;
}

三、結果分析

下圖為上面兩種方法的耗時：

1、執行時間分析：
並行時，隨著進程數目的增多，並行計算的時間越來越短；當達到一定的進程數時，執行時間小到最小值；然後再隨著進程數的增多，執行時間反而越來越長。
2、加速比分析：
隨著進程數的增大，加速比也是逐漸增大到最大值；再隨著進程數的增大，加速比逐漸減小。
3、執行效率分析：
隨著進程數的增大，程式執行效率不斷降低

由於消息傳遞需要成本，而且不是每個進程都同時開始和結束，所以隨著進程數的上升，平均每進程的效率下降

四、頭文件functions.h內容

/********** 輸出函數 **********/
void Matrix_print(double *A,int M,int N)
{
   for(int i=0;i<M;i++){
      for(int j=0;j<N;j++)
         printf("%.1f ",A[i*N+j]);
      printf("\n");
   }
      printf("\n");
}

　　

結束。

Tags: MPI

Previous post

華為補貼活動上線：買手機、智慧屏便宜了

Next post

成立20多年曾經中國第一電商網站倒閉：原因一聲嘆息