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;
matrixLenB=10;
}
else
{
nodeNum=atoi(argv[1]);
matrixHighA=atoi(argv[2]);
matrixLenA=atoi(argv[3]);
matrixHighB=atoi(argv[3]);
matrixLenB=atoi(argv[4]);
}
int p=sqrt(nodeNum);
/*************計算local A B行列數*************************/
int localHighA=matrixHighA/p;
int localLenA=matrixLenA/p;
int localHighB=matrixHighB/p;
int localLenB=matrixLenB/p;
int i;
int j;
int k;
int l;
int * A=(int *)malloc(localLenA*localHighA*sizeof(int));
RandomMatrix(A,localHighA*localLenA);
int * B=(int *)malloc(localLenB*localHighB*sizeof(int));
RandomMatrix(B,localHighB*localLenB);
int * C=(int *)malloc(localHighA*localLenB*sizeof(int));
for(i=0;i<localHighA*localLenB;i++)C[i]=0;
/* printf("%d local A:\n",rank);
PrintMatrixForVector(A,localHighA,localLenA);
printf("%d local B:\n",rank);
PrintMatrixForVector(B,localHighB,localLenB);*/
int myRow=rank/p;//計算當前節點在mesh中的行列值
int myCol=rank%p;
//將數據發送到0號進程收集顯示
MPI_Send(A,localHighA*localLenA,MPI_INT,0,rank+100,MPI_COMM_WORLD);
MPI_Send(B,localHighB*localLenB,MPI_INT,0,rank+200,MPI_COMM_WORLD);
if(rank==0)
{
int **matrixA=(int **)malloc(matrixHighA*sizeof(int *));
for (i=0;i<matrixHighA;i++)
matrixA[i]=(int *)malloc(matrixLenA*sizeof(int));
int **matrixB=(int **)malloc(matrixHighB*sizeof(int *));
for (i=0;i<matrixHighB;i++)
matrixB[i]=(int *)malloc(matrixLenB*sizeof(int));
for(i=0;i<nodeNum;i++)
{
int *receiveATemp=(int *)malloc(localLenA*localHighA*sizeof(int));
int *receiveBTemp=(int *)malloc(localLenB*localHighB*sizeof(int));
MPI_Recv(receiveATemp,localHighA*localLenA,MPI_INT,i,i+100,MPI_COMM_WORLD,&status);
MPI_Recv(receiveBTemp,localHighB*localLenB,MPI_INT,i,i+200,MPI_COMM_WORLD,&status);
l=0;
for(j=0;j<localHighA;j++)
for(k=0;k<localLenA;k++)
{
matrixA[j+(int)(i/p)*localHighA][k+(int)(i%p)*localLenA]=receiveATemp[l++];
}
l=0;
for(j=0;j<localHighB;j++)
for(k=0;k<localLenB;k++)
{
matrixB[j+(int)(i/p)*localHighB][k+(int)(i%p)*localLenB]=receiveBTemp[l++];
}
free(receiveATemp);
free(receiveBTemp);
}
printf("A:\n");
PrintMatrix(matrixA,matrixHighA,matrixLenA);
printf("B:\n");
PrintMatrix(matrixB,matrixHighB,matrixLenB);
for (i=0;i<matrixHighA;i++)
free(matrixA[i]);
for (i=0;i<matrixHighB;i++)
free(matrixB[i]);
free(matrixA);
free(matrixB);
}
for(i=0;i<p;i++)//每個節點向同行同列發送局部數據 A B
{
// if(myCol!=i)
{
MPI_Send(A,localHighA*localLenA,MPI_INT,myRow*p+i,1,MPI_COMM_WORLD);
MPI_Send(B,localHighB*localLenB,MPI_INT,myRow*p+i,2,MPI_COMM_WORLD);
}
// if(myRow!=i)
{
MPI_Send(A,localHighA*localLenA,MPI_INT,i*p+myCol,1,MPI_COMM_WORLD);
MPI_Send(B,localHighB*localLenB,MPI_INT,i*p+myCol,2,MPI_COMM_WORLD);
}
}
int *receiveA=(int *)malloc(localLenA*localHighA*sizeof(int));
int *receiveB=(int *)malloc(localLenB*localHighB*sizeof(int));
int *resultC= (int *)malloc(localHighA*localLenB*sizeof(int));
for(i=0;i<localHighA*localLenB;i++)resultC[i]=0;
/*********************計算矩陣的值*****************************/
for(i=0;i<p;i++)
{
MPI_Recv(receiveA,localHighA*localLenA,MPI_INT,myRow*p+i,1,MPI_COMM_WORLD,&status);
MPI_Recv(receiveB,localHighB*localLenB,MPI_INT,i*p+myCol,2,MPI_COMM_WORLD,&status);
MatrixMultiply(receiveA,receiveB,resultC,localHighA,localLenA,localLenB);
/* if(rank==0)
{
printf("%d A:\n",i);
PrintMatrixForVector(receiveA,localHighA,localLenA);
printf("%d B:\n",i);
PrintMatrixForVector(receiveB,localHighB,localLenB);
printf("%d C:\n",i);
PrintMatrixForVector(resultC,localHighA,localLenB);
}*/
MatrixAdd(C,resultC,localHighA,localLenB);
}
MPI_Send(C,localHighA*localLenB,MPI_INT,0,rank+400,MPI_COMM_WORLD);//將局部結果C發送至0號收集
if(rank==0)//收集數據並且在後面顯示
{
int **matrixC=(int **)malloc(matrixHighA*sizeof(int *));
for (i=0;i<matrixHighA;i++)
matrixC[i]=(int *)malloc(matrixLenB*sizeof(int));
for(i=0;i<nodeNum;i++)
{
int *receiveCTemp=(int *)malloc(localLenB*localHighA*sizeof(int));
MPI_Recv(receiveCTemp,localHighA*localLenB,MPI_INT,i,i+400,MPI_COMM_WORLD,&status);
l=0;
for(j=0;j<localHighA;j++)
for(k=0;k<localLenB;k++)
{
matrixC[j+(int)(i/p)*localHighA][k+(int)(i%p)*localLenB]=receiveCTemp[l++];
}
free(receiveCTemp);
}
printf("C:\n");
PrintMatrix(matrixC,matrixHighA,matrixLenB);
}
MPI_Finalize();
return 0;
}
編譯 mpicc -o summa summa.c -lm
運行(需要指定處理節點的個數一般為可開方整數)
mpirun -np 處理器個數 ./summa 處理器個數 A行數 A列數 B列數
例如運行
結果如下
