MPI 实现 SUMMA 矩阵乘法

SUMMA 算法

SUMMA 算法和Fox 算法一样,将A , B 和C 划分为相同大小的矩阵,对应放在r×c 二维 mesh 上. 但
SUMMA 算法将矩阵乘法分解为一系列的秩nb 修正,即各处理器中的A 和B 分别被分解为nb 大小的列块
和行块进行相乘, nb≤min( k/ r , k / c) , 前面所说的分块尺寸就是指nb 的大小. 算法中, 广播实现为逻辑处理
器行环或列环上的流水线传送, 达到了计算与通信的重叠.具体描述如算法1所示.

算法1. 二维mesh 上的SUMMA 算法
C= 0
for i= 0 t o k- 1 step nb　do
　cur col = i×c/ n
　cur r ow = i×r / m
　if my col = curr ol 向本行广播 A 从 i mo d ( k/ c) 列开始的 nb 列,存于 A′
　if my r ow = curr ow 向本列广播 B 从 i mod ( k / r )行开始的 nb 行,存于 B′
　C= C+ A′ ×B′
end fo r
SUMMA 算法的计算量为 ( n^3/ p ) . 算法中通信部分采用了流水线广播的技术, 其通信开销应为( k/ nb+ 2c- 3) ( ts + m×nb tw / r) + ( k / nb+ 2r - 3) ( t s+ n×nbtw / c) ,当m= k= n 且r = c= p 时通信开销为2( n/ nb+ 2 p - 3) ( ts+ n×nb tw / p ) .由于每个处理器上需要存放A , B 和C 三个子矩阵,再加上接收消息所需空间,总共需要的空间为 ( 3^n2/ p + 2nb×n/ p ) .
SUMMA 算法的计算复杂度和Fo x 算法相当,而所需的辅助空间小于Fox 算法, 在处理器数目相同时,
SUMMA 算法可以求解更大规模的问题.另外又由于SUMMA 算法采用了流水线广播的技术,其通信开销
中不含lo gP 因子,故SUMMA 算法具有更好的可扩放性,对于大规模并行机, SUMMA 算法要优于Fox 算法.

实现代码：

[cpp]
/****************summa***************************************/
/********** 2012.11.17   YingfengChen    *******************/
/************************************************************/
#include <mpi.h>
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <time.h>
#include <sys/time.h>

void PrintMatrixForVector(int * matrix,int high,int len)
{
    int i;
    for(i=0;i<high*len;i++)
    {

        printf("%6d ",matrix[i]);

        if(i%len==len-1&&i!=0)
            printf("\n");


    }


}

/*******可以利用open MP并行优化**********/
void MatrixMultiply(int * A,int *B,int *C,unsigned m,unsigned n,unsigned p)
{   int i,j,k;

    /*printf("A: \n");
    PrintMatrixForVector(A,m,n);

    printf("B: \n");
    PrintMatrixForVector(B,n,p);*/

    for(i=0;i<m;i++)
       for(j=0;j<p;j++)

       {
           int result=0;
           for(k=0;k<n;k++)
            {
               result=A[i*n+k]*B[k*p+j]+result;

            }
           C[i*p+j]=result;


       }


}

/*******可以利用open MP并行优化**********/
void MatrixAdd(int * A,int *B,unsigned m,unsigned n) //the result remain in A
{ int i,j;
   for(i=0;i<m;i++)
      for(j=0;j<n;j++)
        {
            A[i*n+j]=A[i*n+j]+B[i*n+j];
        }

}

void PrintMatrix(int ** matrix,int high,int len)
{
    int i,j;
    for(i=0;i<high;i++)
    {
        for(j=0;j<len;j++)
        {
            printf("%6d ",matrix[i][j]);

        }
        printf("\n");
    }

}

/****随机数据以待计算****/
void RandomMatrix(int *matrix,int len)
{
   struct timeval tpstart; //使数据更均匀
   gettimeofday(&tpstart,NULL);
   srand(tpstart.tv_usec);
   int i=0;
   for(i=0;i<len;i++)
   matrix[i]=rand()%8;

}
int main(int argc,char **argv)
{

    int rank;
    MPI_Status status;
    MPI_Init(&argc,&argv);
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);

    int nodeNum; //节点数，必须为可开平方

    int matrixHighA;// 矩阵A行数
    int matrixLenA;//矩阵A列数


    int matrixHighB;
    int matrixLenB;

    /****   参数检查参数错误用默认参数*****/
    if(argc!=5&&rank==0)
    {

         printf("The para is wrong!using default para\n");

         nodeNum=4;

         matrixHighA=6;
         matrixLenA=8;


         matrixHighB=8;
&nb

补充：软件开发 , C++ ,