• MPI学习笔记（二）：矩阵相乘的两种实现方法

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

            最近领导让把之前安装的软件lapack、blas里的dgemm运算提取出来独立作为一套程序，然后把这段程序改为并行的，并测试一下进程规模扩展到128时的并行效率。
            我发现这个是dgemm.f文件，里面主要是对C=αAB+βC的实现，因此在此总结一下MPI的矩阵乘法使用。
            其主要思想：是把相乘的矩阵按行分解（任务分解），分别分给不同的进程，然后在汇总到一个进程上，在程序上实现则用到了主从模式，人为的把进程分为主进程和从进程，主进程负责对原始矩阵初始化赋值，并把数据均匀分发（为了负载均衡）到从进程上进行相乘运算，主要用到的知识是MPI点对点通信和组通信的机制。

    一、使用简单的MPI_Send和MPI_Recv实现

    ?
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    47
    48
    49
    50
    51
    52
    53
    54
    55
    56
    57
    58
    59
    60
    61
    62
    63
    64
    65
    66
    67
    68
    69
    70
    71
    72
    73
    74
    75
    76
    77
    78
    79
    80
    81
    82
    83
    84
    85
    86
    87
    88
    89
    90
    91
    92
    93
    94
    95
    96
    97
    98
    99
    100
    101
    102
    103
    104
    105
    106
    107
    108
    109
    110
    111
    #include
    #include "mpi.h"
    #include
    #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
             for(int j=0;j
                Matrix_A[i*N+j]=i+1;
          }
          for(int i=0;i
             for(int j=0;j
                Matrix_B[i*K+j]=j+1;
          }
          for(int i=0;i
             for(int j=0;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
             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
             MPI_Recv(ans,line*K,MPI_DOUBLE,p,33,MPI_COMM_WORLD,&status);
             for(int i=0;i
                for(int j=0;j
                   Matrix_C[((p-1)*line+i)*K+j]=ans[i*K+j];
          }
     
          // 计算A剩下的行数据
          for(int i=(comm_sz-1)*line;i
             for(int j=0;j
                double temp=0;
                for(int p=0;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
             for(int j=0;j
                double temp=0;
                for(int p=0;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实现

    ?
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    47
    48
    49
    50
    51
    52
    53
    54
    55
    56
    57
    58
    59
    60
    61
    62
    63
    64
    65
    66
    67
    68
    69
    70
    71
    72
    73
    74
    75
    76
    77
    78
    79
    80
    81
    82
    83
    84
    85
    86
    87
    88
    89
    90
    91
    92
    93
    94
    95
    96
    97
    98
    99
    #include
    #include "mpi.h"
    #include
    #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
             for(int j=0;j
                Matrix_A[i*N+j]=i+1;
             for(int p=0;p
                Matrix_C[i*K+p]=1;
          }
          for(int i=0;i
             for(int j=0;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
          for(int j=0;j
             double temp=0;
             for(int p=0;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
                for(int j=0;j
                   double temp=0;
                   for(int p=0;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内容

    ?
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    /********** 输出函数 **********/
    void Matrix_print(double *A,int M,int N)
    {
       for(int i=0;i
          for(int j=0;j
             printf("%.1f ",A[i*N+j]);
          printf("\n");
       }
          printf("\n");
    }

     五、对以上两种程序的不同写法

    1、MPI_Send和MPI_Recv

    ?
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    47
    48
    49
    50
    51
    52
    53
    54
    55
    56
    57
    58
    59
    60
    61
    62
    63
    64
    65
    66
    67
    68
    69
    70
    71
    72
    73
    74
    75
    76
    77
    78
    79
    80
    81
    82
    83
    84
    85
    86
    87
    88
    89
    90
    91
    92
    93
    94
    95
    96
    97
    98
    99
    100
    101
    102
    103
    104
    105
    106
    107
    108
    109
    110
    111
    112
    113
    114
    #include
    #include "mpi.h"
    #include
    #include "dgemm_1.h"
     
    #define M 1200 // 矩阵维度
    #define N 1000
    #define K 1300
     
    int main(int argc, char **argv)
    {
       int my_rank,comm_sz,line;
       double start, stop; //计时时间
       double alpha=2,beta=2; // 系数C=aA*B+bC
       MPI_Status status;
     
       MPI_Init(&argc,&argv);
       MPI_Comm_size(MPI_COMM_WORLD, &comm_sz);
       MPI_Comm_rank(MPI_COMM_WORLD,&my_rank);
     
       line=M/comm_sz; // 每个进程分多少行数据
     
       // 给矩阵A B,C赋值
       if(my_rank==0){
          double *Matrix_A,*Matrix_B,*Matrix_C,*result_Matrix;
          Matrix_A=(double*)malloc(M*N*sizeof(double));
          Matrix_B=(double*)malloc(N*K*sizeof(double));
          Matrix_C=(double*)malloc(M*K*sizeof(double));
          result_Matrix=(double*)malloc(M*K*sizeof(double)); // 保存数据计算结果
     
          start=MPI_Wtime();
          for(int i=0;i
             for(int j=0;j
                Matrix_A[i*N+j]=i+1;
          }
          for(int i=0;i
             for(int j=0;j
                Matrix_B[i*K+j]=j+1;
          }
          for(int i=0;i
             for(int j=0;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
             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 i=0;i
             MPI_Recv(result_Matrix+i*line*K,line*K,MPI_DOUBLE,i+1,33,MPI_COMM_WORLD,&status);
          }
     
          // 计算A剩下的行数据
          //matMulti(Matrix_A+(comm_sz-1)*line,Matrix_B,result_Matrix,M,N,K,alpha,beta);
          for(int i=(comm_sz-1)*line;i
             for(int j=0;j
                double temp=0;
                for(int p=0;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(result_Matrix);
       }
       else{
          double *Matrix_B,*buffer_A,*buffer_C,*ans;
          Matrix_B=(double*)malloc(N*K*sizeof(double));
          buffer_A=(double*)malloc(line*N*sizeof(double)); // A的均分行的数据
          buffer_C=(double*)malloc(line*K*sizeof(double)); // C的均分行的数据
     
          //接收广播的数据
          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);
     
          //计算乘积结果,并将结果发送给主进程
          //matMulti(buffer_A,Matrix_B,buffer_C,line,N,K,alpha,beta);
          for(int i=0;i
             for(int j=0;j
                double temp=0;
                for(int p=0;p
                   temp+=buffer_A[i*N+p]*Matrix_B[p*K+j];
                }
                buffer_C[i*line+j]=alpha*temp+beta*buffer_C[i*K+j];
             }
          }
          MPI_Send(buffer_C,line*K,MPI_DOUBLE,0,33,MPI_COMM_WORLD);
     
          free(Matrix_B);
          free(buffer_A);
          free(buffer_C);
       }
     
       MPI_Finalize();
       return 0;
    }

    2、MPI_Scatter和MPI_Gather

    在上面的MPI_Scatter和MPI_Gather的程序中，矩阵的行数据并不能平均的分配给每个子进程，当不能矩阵的行数据不能与总进程数整除时，就需要单独计算一下矩阵剩余的行数据。

    ?
    1
    2
    3
    4
    if(M%comm_sz!=0){
        M-=M%comm_sz;
        M+=comm_sz;
    }

    　　为了使矩阵的行数据能平均的分配给每个子进程，因此记comm_sz为进程数，矩阵的维度需要是comm_sz的倍数。我们将矩阵的维度扩展到comm_sz的倍数，多余的部分用0填充，保证正确性。

    ?
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    47
    48
    49
    50
    51
    52
    53
    54
    55
    56
    57
    58
    59
    60
    61
    62
    63
    64
    65
    66
    67
    68
    69
    70
    71
    72
    73
    74
    75
    76
    77
    78
    79
    80
    81
    82
    83
    84
    85
    86
    87
    88
    89
    90
    91
    92
    93
    94
    95
    96
    97
    98
    99
    100
    101
    102
    103
    104
    105
    106
    107
    108
    109
    110
    111
    112
    113
    114
    115
    116
    117
    118
    119
    120
    121
    122
    123
    124
    125
    126
    127
    128
    129
    130
    #include
    #include "mpi.h"
    #include
    #include "dgemm_1.h"
     
    //#define M 3 // 矩阵维度
    #define N 1000
    #define K 1300
     
    int main(int argc, char **argv)
    {
       int M=1400;
       int my_rank,comm_sz,line;
       double start, stop; //计时时间
       MPI_Status status;
       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);
     
       if(comm_sz==1){
          double *Matrix_A,*Matrix_B,*Matrix_C;
     
          Matrix_A=(double*)malloc(M*N*sizeof(double));
          Matrix_B=(double*)malloc(N*K*sizeof(double));
          Matrix_C=(double*)malloc(M*K*sizeof(double));
          start=MPI_Wtime();
          // 给矩阵A B,C赋值
          for(int i=0;i
             for(int j=0;j
                Matrix_A[i*N+j]=i+1;
             for(int p=0;p
                Matrix_C[i*K+p]=1;
          }
          for(int i=0;i
             for(int j=0;j
                Matrix_B[i*K+j]=j+1;
          }
     
          // 计算结果
          for(int i=0;i
             for(int j=0;j
                double temp=0;
                for(int p=0;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);
       }
       else{
          double *Matrix_A,*Matrix_B,*Matrix_C,*ans,*buffer_A,*buffer_C,*result_Matrix;
          int saveM=M;
          if(M%comm_sz!=0){
             M-=M%comm_sz;
             M+=comm_sz;
          }
          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
                for(int j=0;j
                   if(i
                      Matrix_A[i*N+j]=i+1;
                   else
                      Matrix_A[i*N+j]=0;
                for(int p=0;p
                      Matrix_C[i*K+p]=1;
             }
             for(int i=0;i
                for(int j=0;j
                   Matrix_B[i*K+j]=j+1;
             }
     
          // 输出A,B,C
          //Matrix_print(Matrix_A,saveM,N);
          //Matrix_print(Matrix_B,N,K);
          //Matrix_print(Matrix_C,saveM,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
             for(int j=0;j
                double temp=0;
                for(int p=0;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);
     
          if(my_rank==0){
             //Matrix_print(result_Matrix,saveM,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;
    }

    　　

     

     

     

    结束。

     

  • 相关阅读:
    【网络原理】TCP/IP协议
    JS点击按钮获取列表信息，纯JS代码创建表格，如何引入CSS框架，＜svg＞标签
    第三章微服务配置中心
    Vue2动态显示时间
    智能小车开发
    Linux调试器-gdb使用
    读书笔记：彼得·德鲁克《认识管理》第8章 战略规划：企业家技能
    MySQL数据库索引的数据结构
    [附源码]计算机毕业设计JAVA购买车票系统
    逐字节讲解 Redis 持久化（RDB 和 AOF）的文件格式
  • 原文地址：https://www.cnblogs.com/babyclass/p/16518726.html