矩阵相乘并行算法

并行处理技术

课程设计分析报告

课程设计题目矩阵相乘并行算法设计姓名廖杰

学号M201372880

专业计算机技术

任课教师金海石宣化

所在学院计算机科学与技术学院报告提交日期2014-01-13

一、实验目的

1、学习使用集群；

2、掌握并行处理或分布计算的编程方法；

3、学会以并行处理的思想分析问题。

二、实验要求

1、自行生成矩阵作为算法的输入；

2、使用并行处理技术编程，例如：MPI、OpenMP、MR；

3、矩阵大小至少为1000*1000；

4、加速比越大成绩越高。

三、实验内容

3.1、矩阵的划分：

对于矩阵相乘的并行算法，可以有三种：对矩阵按行划分、按列划分和棋盘式分块划分。和按行或列划分相比，棋盘式划分可以开发出更高的并行度。对于一个n×n的方阵，棋盘划分最多可以使用n^2个处理器进行并行计算，但使用按行或列分解最多可以使用n个。对矩阵相乘采用棋盘式划分的算法通常称作Cannon算法。

A）行列划分

又叫带状划分（Striped Partitioning），就是将矩阵整行或者整列分成若干个组，每个组指派给一个处理器。下图所例为4个CPU，8×8矩阵的带状划分。

在带状划分情况下，每个CPU将会均匀分配到2行(列)数据。8×8矩阵变成了一个1×4或4×1的分块矩阵，每个CPU所属的分块矩阵大小为8×2或2×8。

B）棋盘划分

就是将矩阵分成若干个子矩阵，每个子矩阵指派给一个处理器，此时任一处理器均不包含整行或者整列。下图所示即为4个处理器情况下8×8矩阵的棋盘划分，其中处理器阵列为2×2，每个处理器分配到的子矩阵大小为4×4。

矩阵划分成棋盘状可以和处理器连成二维网孔相对应。对于一个n×n维矩阵和p×p的二维处理器阵列，每个处理器均匀分配有（n/p）×(n/p)=n^2/p^2个元素。使用棋盘式划分的矩阵相乘算法一般有两种，Cannon算法和Summa算法。SUMMA算法能够计算

m*l的A矩阵和l*n的B矩阵相乘（m、l、n可不相等），而cannon算法只能实现n*n 的A矩阵和n*n的B矩阵相乘，具有很大的局限性。

3.2、算法原理

A) 行划分法

假设是M*N，计算前，将矩阵N发送给所有从进程，然后将矩阵M分块，将M中数据按行分给各从进程，在从进程中计算M中部分行数据和N的乘积，最后将结果发送给主进程。这里为了方便，有多少进程，就将M分了多少块，除最后一块外的其他数据块大小都相等，最后一块是剩下的数据，大小大于等于其他数据块大小，因为矩阵行数不一定整除进程数。最后一块数据在主进程中计算，其他的在从进程中计算。

定义两个矩阵M和N，N所有进程都需要，M可以只在主进程中定义。其他的变量视主进程和从进程需要按要求定义在合适的位置。

代码参见附录部分。

B) Cannon算法

Cannon算法的基本思想可以如下表示：假设两个矩阵A和B相乘，把A和B矩阵划分成p个方块，进程的编号从到，并在最初把子矩阵和分配给。虽然第i行的每个进程需要全部的个子矩阵，但我们还是能调度第i行个进程的计算，使得

每个进程在任何时刻都是用不同的。每完成一次矩阵乘法，这些块在各进程之间被轮流使用，似的每次轮流之后每个进程都可以得到新的。对列使用同样的调度，则在任何时

刻，任何进程至多拥有每个矩阵的一个块，在所有进程中，改算法需要的总内存量为。下图为此算法中不同进程上子矩阵乘法的调度过程。

假如矩阵C=A*B,则C的的计算公式如下：

进程P 存储分块矩阵这一部分。块矩阵乘法要计算所有匹配的和，然而只有在主对

角线的才是匹配的。因此需要采用循环移动分块矩阵的方法来使每个进程都有一对可以直接相乘的匹配的块，具体方法如下：

(1)将排第i行的块循环左移i个位置，将第列．块循环上移j个位置；

(2) 进程执行乘一加运算，然后将移动得到的块循环左移1个位置，将移动得到的

块循环上移1个位置；

(3)重复第2步(一1)次，每次移动后进程执行乘一加运算。

经过以上操作后就可以得到矩阵C的解。

代码请参见附录部分

C) Summa算法

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

C= 0

for i= 0 t o k-1 step nb do

cur col = i×c/ n

cur row = i×r / m

if my col = cur rol 向本行广播A 从i mod (k/c) 列开始的nb 列, 存于A′

if my row = cur row 向本列广播B 从i mod (k/r) 行开始的nb 行, 存于B ′

C= A′×B ′

end for

SUMMA算法的核心思想是：各处理器收集来自同一行处理器中A矩阵子块的所有列和同一列处理器中B矩阵子块的所有行，然后将行列相乘后累加，形成一个C矩阵的分块矩阵。最后由rank=0的处理器将其他处理器的数据收集起来，形成最终的矩阵C。

SUMMA算法相较于cannon算法的优势只要体现在SUMMA算法能够计算m*l的A 矩阵和l*n的B矩阵相乘（m、l、n可不相等），而cannon算法只能实现n*n的A矩阵和n*n的B矩阵相乘，具有很大的局限性。

代码参见附录部分。

3.3、程序运行结果对比分析

A) 统一的实验条件

矩阵大小：1000*1000；

矩阵数字范围：0~10；

矩阵数字分布是否随机：是；

分配的进程数：9；

B) 实验进程数解释

由于Cannon算法本身局限性，要使用Cannon算法，必须满足进程数为整数的平方，比如1、4、9、16等。在本次的实验环境之下，经过多次对比分析，发现对于分行还是分块算法，进程数安排在8~15可以得到最理想的运行速度：进程数目过小则每个进程单独运算的时间过多，进程数过大则选路时间（进程与进程之间的通信时间）过长。而对比要求每个算法的进程相同，故此处选择进程数目为9.

C) 算法运行时间对比

Cannon算法运行时间如下：

分行法运行时间如下：

串行算法运行时间如下：

由于Summa算法与Cannon算法思路几乎相同，而且在算法预处理阶段要比Cannon 算法更耗时，故没有做多余的实验。

显而易见，单纯的运用分行算法所花费的时间是最短的。

D) 结果分析

Cannon算法相对于简单的行划分并行处理算法，其优势仅仅在于并行度可以更高（可达到N*N个进程，N为矩阵宽），但在并行度相同的情况下，其多出的预处理过程、矩阵发送与结果回收机制会占用更多的时间。

3.4、程序调优

A) 行划分算法优化

1. 循环优化

在预估计矩阵大小为10的倍数的基础上，对每一个步长为1的循环做处理，改为步长为10的循环，将十次循环体全部压缩在一次循环中，从而大量减少了循环的判别时间，提升循环运算速度。例如在单个线程在计算部分结果时，采用的循环为：

for(i=0;i

for(j=0;j

DATA temp=0;

for(k=0;k

temp += buffer[i*width+k]*n[j*width+k];

temp += buffer[i*width+k+1]*n[j*width+k+1];

temp += buffer[i*width+k+2]*n[j*width+k+2];

temp += buffer[i*width+k+3]*n[j*width+k+3];

temp += buffer[i*width+k+4]*n[j*width+k+4];

temp += buffer[i*width+k+5]*n[j*width+k+5];

temp += buffer[i*width+k+6]*n[j*width+k+6];

temp += buffer[i*width+k+7]*n[j*width+k+7];

temp += buffer[i*width+k+8]*n[j*width+k+8];

temp += buffer[i*width+k+9]*n[j*width+k+9];

}

ans[i*width+j] = temp;

}

}

在将循环次数压缩的同时，为了进一步减少循环的运算量，在每一个步长为10的循环之前做预处理，避免循环体中的重复运算。例如在主进程在接受其他进程时，将结果矩阵整合的过程：

for(k=1;k

{

MPI_Recv(ans,line*width,MPI_INT,k,2,MPI_COMM_WORLD,&status);

for(i=0;i

{

count=i*k*width; //将i*k*width提前算好，减少了下一步循环的重复运算

count1=i*width;

for(j=0;j

p[count+j] = ans[count1+j];

p[count+j+1] = ans[count1+j+1];

p[count+j+2] = ans[count1+j+2];

p[count+j+3] = ans[count1+j+3];

p[count+j+4] = ans[count1+j+4];

p[count+j+5] = ans[count1+j+5];

p[count+j+6] = ans[count1+j+6];

p[count+j+7] = ans[count1+j+7];

p[count+j+8] = ans[count1+j+8];

p[count+j+9] = ans[count1+j+9];

}

}

}

2. 节省空间

在进行矩阵工作量划分并传送的时候，为每一个进程开辟仅仅是自己所需要大小的空间，例如在9进程的环境下，每个进程所需要接受的缓存空间为B矩阵大小以及大约1/9大小A 矩阵。

内存开辟：

buffer = (DATA *)malloc(sizeof(DATA)*width*line);

矩阵A分块传输：

for(k=1;k

{

for(i=k;i

{

count=i/numprocs*width;

count1=i*width;

for(j=0;j

{

buffer[count+j]=m[count1+j];

buffer[count+j+1]=m[count1+j+1];

buffer[count+j+2]=m[count1+j+2];

buffer[count+j+3]=m[count1+j+3];

buffer[count+j+4]=m[count1+j+4];

buffer[count+j+5]=m[count1+j+5];

buffer[count+j+6]=m[count1+j+6];

buffer[count+j+7]=m[count1+j+7];

buffer[count+j+8]=m[count1+j+8];

buffer[count+j+9]=m[count1+j+9];

}

}

MPI_Send(buffer,line*width,MPI_INT,k,1,MPI_COMM_WORLD);

同样的方式也运用在运行空间的开辟上。

这样做不仅仅是内存空间的节约，同时也减少了进程之间的数据传输量，大大节省了进程之间的协作时间！

B) 稀疏矩阵处理

虽然程序并未对稀疏矩阵进行优化，但是还是试着对程序的输入数据模式进行更改，体验一下稀疏矩阵运算的速度提升有多快。

void magten(DATA *a,int width)

{

int i,j;

srand((unsigned)time(NULL));

for(i=0;i

for(j=0;j

a[i*width+j] = (rand()%MUL)*(rand()%2)*(rand()%3);

a[i*width+j+1] = (rand()%MUL)*(rand()%2)*(rand()%3);

a[i*width+j+2] = (rand()%MUL)*(rand()%2)*(rand()%3);

a[i*width+j+3] = (rand()%MUL)*(rand()%2)*(rand()%3);

a[i*width+j+4] = (rand()%MUL)*(rand()%2)*(rand()%3);

a[i*width+j+5] = (rand()%MUL)*(rand()%2)*(rand()%3);

a[i*width+j+6] = (rand()%MUL)*(rand()%2)*(rand()%3);

a[i*width+j+7] = (rand()%MUL)*(rand()%2)*(rand()%3);

a[i*width+j+8] = (rand()%MUL)*(rand()%2)*(rand()%3);

a[i*width+j+9] = (rand()%MUL)*(rand()%2)*(rand()%3);

}

}

}

如上面所示，对于每一个生成的数据都再一次进行乘法运算，其中乘数是0的概率为2/3.

运行结果如下：

可以看出，稀疏矩阵的乘法由0.76s变为0.75s，仅仅是短暂的提升。

提升计时显示的精度，可以看到，对于稀疏矩阵的处理要比普通矩阵快0.015s，提速约为2%

3.5、算法最优速度

对算法运用不同的进程数目运算进行了大量重复试验，最终得出在进程数大概为12的时候，本算法的运行速度最快。最终结果如下：

发出工作量时间为0.138184s

运算时间为0.495569s

接收答案时间为0.025461s

总运算时间0.659240s

四、总结分析

进程数目12 1

加速比15.80 1

效率15.80/12=1.317 1

从效率大于1上可以看出，本次课程设计做出的算法为超线性加速，这主要得益于对循环体的优化。

附录

Cannon算法

#include

#include "mpi.h"

#include

#include

#include

#include

#define MUL 10

MPI_Status status;

double **A, **B, **C; //C=A*B

double *a,*b,*c; //各个进程的缓冲区

int n; //矩阵的行列数

int np; //每个进程控制的小矩阵的行列数

int p,rank; //进程个个数、当前进程的编号,笛卡尔进程编号

double *tempa, *tempb;

void ProduceABC(); //在根处理器中生成矩阵AB,初始化矩阵C

void PrintABC();//输出结果

void ScatterAB();// 分发矩阵AB中的元素到各个进程中

void MainProcess(); //cannon算法的主过程

void collectC(); //收集结果矩阵C

void Mutiply(); //矩阵相乘

void Printab();

void Printc();

int main(int argc, char *argv[])

{

int i;

double starttime,endtime;

MPI_Init(&argc, &argv);

MPI_Comm_size(MPI_COMM_WORLD, &p);

MPI_Comm_rank(MPI_COMM_WORLD, &rank);

if(rank == 0)

{

printf("please input the raw number:n= ");

fflush(stdout);

scanf("%d", &n);

printf("\n");

}

MPI_Bcast(&n, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);

// n = atoi(argv[1]);

np = n/(int)sqrt(p);

a = (double*)malloc(np*np*sizeof(double));

b = (double*)malloc(np*np*sizeof(double));

c = (double*)malloc(np*np*sizeof(double));

memset(c, 0, np*np*sizeof(double));

tempa = (double*)malloc(np*np*sizeof(double));

tempb = (double*)malloc(np*np*sizeof(double));

if(rank == 0)

{

//在根处理器中为矩阵ABC分配空间

A = (double**)malloc(n*sizeof(double*));

B = (double**)malloc(n*sizeof(double*));

C = (double**)malloc(n*sizeof(double*));

for(i = 0; i < n; i++)

{

A[i] = (double*)malloc(n*sizeof(double));

B[i] = (double*)malloc(n*sizeof(double));

C[i] = (double*)malloc(n*sizeof(double));

}

ProduceABC(); //在根处理器中随机生成矩阵AB，初始化矩阵C

ScatterAB();// 分发矩阵AB中的元素到各个进程中

}

else

{

MPI_Recv(a, np*np, MPI_DOUBLE, 0, 1, MPI_COMM_WORLD, &status);

MPI_Recv(b, np*np, MPI_DOUBLE, 0, 2, MPI_COMM_WORLD, &status); }

starttime=MPI_Wtime(); //开始时间

MainProcess(); //cannon算法的主过程

if(rank == 0)

{

collectC(); //收集结果矩阵C

PrintABC(); //输出结果

endtime=MPI_Wtime();

printf("time used: %lf\n",endtime - starttime);

for(i = 0; i < n; i++)

{

free(A[i]);

free(B[i]);

free(C[i]);

}

free(A);

free(B);

free(C);

}

else

{

MPI_Send(c, np*np, MPI_DOUBLE, 0, 1, MPI_COMM_WORLD);

}

free(a);

free(b);

free(c);

free(tempa);

free(tempb);

MPI_Finalize();

return 0;

}

void ProduceABC()//在根处理器中生成矩阵AB

{

int i,j;

srand((unsigned)time(NULL));

for(i=0; i

for(j=0; j

{

A[i][j]=rand()%MUL;

B[i][j]=rand()%MUL;

C[i][j]=0.0;

}

}

void PrintABC()//输出结果

{

printf("A[0][0]=%f\nB[0][0]=%f\nC[0][0]=%f\n",A[0][0],B[0][0],C[0][0]);

}

void ScatterAB()// 分发矩阵AB中的元素到各个进程中

{

int imin,imax,jmin,jmax;

int sp;

int i,j,k,m;

for(k=0; k

{

/*计算相应处理器所分得的矩阵块在总矩阵中的坐标范围*/

sp = (int)sqrt(p);

imin = (k / sp) * np;

imax = imin + np - 1;

jmin = (k % sp) * np;

jmax = jmin + np -1;

/*rank=0的处理器将A,B中的相应块拷至tempa,tempb，准备向其他处理器发送*/

m = 0;

for(i=imin; i<=imax; i++)

{

for(j=jmin; j<=jmax; j++)

{

tempa[m] = A[i][j];

tempb[m] = B[j][i]; //矩阵B按列优先存储

m++;

}

}

/*根处理器将自己对应的矩阵块从tempa,tempb拷至a,b*/

if(k==0)

{

memcpy(a, tempa, np*np*sizeof(double));

memcpy(b, tempb, np*np*sizeof(double));

}

else /*根处理器向其他处理器发送tempa,tempb中相关的矩阵块*/

{

MPI_Send(tempa, np*np, MPI_DOUBLE, k, 1, MPI_COMM_WORLD);

MPI_Send(tempb, np*np, MPI_DOUBLE, k, 2, MPI_COMM_WORLD);

}

}

}

void MainProcess() //cannon算法的主过程

{

MPI_Comm comm; //笛卡尔结构通讯器

int crank;

int dims[2],periods[2], coords[2];

int source, dest, up, down, right, left;

int i;

dims[0] = dims[1] = (int)sqrt(p);

periods[0] = periods[1] = 1;

MPI_Cart_create(MPI_COMM_WORLD, 2, dims, periods, 1, &comm);

MPI_Comm_rank(comm, &crank);

MPI_Cart_coords(comm, crank, 2, coords);

MPI_Cart_shift(comm, 1, -1, &right, &left);

MPI_Cart_shift(comm, 0, -1, &down, &up);

MPI_Cart_shift(comm, 1, -coords[0], &source, &dest);

MPI_Sendrecv_replace(a, np*np, MPI_DOUBLE, dest, 1, source, 1, comm, &status);

MPI_Cart_shift(comm, 0, -coords[1], &source, &dest);

MPI_Sendrecv_replace(b, np*np, MPI_DOUBLE, dest, 1, source, 1, comm, &status);

Mutiply(); //矩阵相乘

for(i = 1; i < dims[0]; i++)

{

MPI_Sendrecv_replace(a, np*np, MPI_DOUBLE, left, 1, right, 1, comm, &status);

MPI_Sendrecv_replace(b, np*np, MPI_DOUBLE, up, 1, down, 1, comm, &status);

Mutiply(); //矩阵相乘

}

MPI_Comm_free(&comm);

}

void collectC() //收集结果矩阵C

{

int i,j,k,s,m;

int imin,imax,jmin,jmax;

int sp= (int)sqrt(p);

/* 根处理器中的c赋给总矩阵C */

for (i=0;i

{

for(j=0;j

C[i][j]=c[i*np+j];

}

for (k=1;k

{

/*根处理器从其他处理器接收相应的分块c*/

MPI_Recv(c, np*np, MPI_DOUBLE, k, 1, MPI_COMM_WORLD, &status);

//printf("rank = %d\n", k);Printc();

imin = (k / sp) * np;

imax = imin + np - 1;