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
#include
#include
#include
#include
#include
void PrintMatrixForVector(int * matrix,int high,int len)
{
int i;
for(i=0;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
{
int result=0;
for(k=0;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
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
for(j=0;j
printf("%6d ",matrix[i][j]);
}
printf("\n");
}
}
/****随机数据以待计算****/
void RandomMatrix(int *matrix,int len)
{
struct timeva l tpstart; //使数据更均匀
gettimeofday(&tpstart,NULL);
srand(tpstart.tv_usec);
int i=0;
for(i=0;i
}
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;
matrixLenB=10;
}
else
{
nodeNum=atoi(argv[1]);
matrixHighA=atoi(argv[2]);
matrixLenA=atoi(argv[3]);
matrixHighB=atoi(argv[3]);
matrixLenB=atoi(argv[4