矩阵计算
矩阵计算问题有很多种类型，例如：
求解线性代数方程组 Ax = b
线性最小二乘问题 given b in R^m， for x in R^n，minimize ||Ax - b||^2
矩阵特征值问题 Ax = λx
矩阵奇异值分解 A = U∑V^T
很多矩阵计算问题都有并行的计算方法，例如矩阵乘法，我们现在来学习并且用代码实现他们，从而更深地理解并行计算的思想。
并行矩阵乘法
并行计算，就是多个进程并行协作，完成特定的任务。现在我们假定一个并行系统，包含了p个处理机，每个处理机一个进程，我们可以分别用字符“0”，“1”，...，“p-1”来引用它们，或者为了清晰，我们用 Pi 来引用它们，i 表示一个进程的进程号，进程之间可以相互传递消息，所谓消息，指的是一个数据结构。设用x引用一个消息，那么函数 send(x, i)表示当前进程(处理机)把消息x发送给序号等于i的进程，recv(x, i)表示当前进程从序号为i的进程接收消息x。在并行编程中，我们用程序代码定义好一个过程，每个进程都将运行这段程序代码定义的过程，也就是说，代码必须是通用的。并行编程的是一个定义各个局部之间相互关系的过程，全局的发展方向不仅仅依赖于局部与局部之间的这些关系，而且对初始状态也是极其依赖的，类似于数学和物理学中的混沌现象，我们需要选择合适的初始状态和局部之间的相互作用关系，使得整体得到希望的发展方向。
接下来我们来说明用这p个处理机来进行并行算法。给出两个矩阵A和B，要C = A×B，我们可以用矩阵乘法的标准定义，还可以用分块矩阵来先进行变换。
算法1：行列划分算法
我们对A进行行分块，对B进行列分块：
这两个矩阵的子矩阵分别两两相乘得到 Ai×Bj = Ci,j ，我们可以得到p×p个矩阵，这些矩阵拼接起来，就得到了结果C，亦即，C = [Ci,j]p×p。我们可以用举世闻名的数学归纳法来证明如此分块的正确性，不过这并不是本文重点，不再赘述。
由此我们得到如下伪代码：
PROCEDURE 1 INPUT: A , B . A and B are both in processor 0. OUTPUT: C = A × B in processor 0. BEGIN:
#数据散发阶段，进程0把矩阵A和B分发给各个进程 if myid == 0: for i in range(1,p): send(Ai, i), send(Bi, i) barrier() else: recv(Amyid, 0), recv(Bmyid, 0) barrier()
#并行计算各个分块矩阵的阶段 left = (myid - 1)%p , right = (myid + 1)%p for i in range(0, p): k = (i + myid)%p C[myid][k] = Amyid * Bk if i != p-1 : send(Bk,left), recv(Bk+1,right)
#聚集阶段，每个节点把自己的计算结果发送给进程0 if myid == 0: for i in range(1, p): j = (i - 1)%p start_new_thread(recv, (C[j],i)) else: j = (myid - 1)%p send(C[j],0) return C END.
在初始状态时，矩阵A和B都在节点0，我们通过数据散发操作使得通信器达到这样的状态：我们用数字 i 来引用通信器在每次循环通信的状态，用数字 k 引用在此状态下进程在自己的内存中保留着的分块矩阵 Ak 和 Bk。 Ai 和 Bi相乘可以得到C[i][i]，计算和循环通信完成后，B的分块矩阵在节点中完成一次轮换，通信器进入下一个状态( 用数字 i+1 来引用这个状态)，继续计算得到C[i][i+1]，直到计算出每个C的分块，最后每个进程把自己计算出来的分块发送给节点0，完成拼接得到C。
考虑这个算法的计算时间，如果是串行执行乘法，需要的时间是 O(mnk)，p个处理机同时计算的计算时间当然就是O(mnk/p)，注意计算总量仍旧是不变的。
再考虑通信时间。串行算法没有通信时间，对于并行算法，忽略每次通信的启动时间，数据散发阶段的通信的时间复杂度O( (m-m0)k+(n-n0)k )，再计算阶段，共要完成p-1次循环传输，每次总共传输nk个数据，时间复杂度O( (p-1)nk )，聚集阶段的时间复杂度O( (m-m0)n )。再没有采用多线程同时进行计算和传输的情况下，算法的运行时间复杂度就是计算时间与通信时间的和 O( mnk/p + (m-m0)k+(n-n0)k + (p-1)nk + (m-m0)n )。
这个算法是否真的更快，就有待实践检验了。当然上述算法还有很多可以优化的地方，这个算法只是一个思路。比如由于矩阵B已经存在于0引用的那个节点，所以不需要把节点0纳入循环；每个节点计算矩阵的过程中，可以每次计算完一列，就可以采用非阻塞通信或者多线的方式，把这一列发送出去，从而同时进行计算和传输。
MPI来实现行列划分算法
MPI是一组接口标准，调用这些接口可以实现并行编程中的各种动作，如初始化或者结束并行编程环境，进程间传递消息等。
通信器提供进程间通信的基本环境,MPI 程序中所有通信都必须在特定的通信器中完成。 MPI 环境在初始化时(也就是会设置一些全局变量)会自动创建两个通信器,一个称为 MPI_COMM_WORLD,它包含程序中的所有进程,另一个称为 MPI_COMM_SELF,它是每个进程独自构成的、仅包含自己的通信器。MPI 系统提供了一个特殊进程号 MPI_PROC_NULL,它代表空进程 (不存在的进程),与 MPI_PROC_NULL 进行通信相当于一个空操 作,对程序的运行没有任何影响。具体的API可以查阅书籍或者看 http://www.cnblogs.com/xinchrome/p/4859119.html
下面上代码：
#include <mpi.h> #include <stdio.h> #include <stdlib.h> #include <pthread.h> #define M 1000 #define N 1000 #define K 1000 /* 职责：读取数据到矩阵中 */ void read_to_mat(int mat[][],int m,int n){ int i = 0,j = 0; char filename[80]; FILE *fp = NULL; fscanf(stdin,"%s",filename); fp = fopen(filename,"r"); for(i = 0;i < m;i++){ for(j = 0;j < n;j++){ if(EOF == fscanf(fp,"%d",&mat[i][j])){ printf("[ERROR] - filename: %s",filename); exit(1); } } } fclose(fp); return; } void write_to_file(int mat[][],int m,int n){ char filename[] = "result.txt" FILE *fp = fopen(filename,"w"); for(int i = 0;i < m;i++){ if(i > 0) fprintf(fp," "); for(int j = 0;j < n;j++){ fprintf(fp,"%d",mat[i][j]); if(j != n-1) fprintf(fp," "); } } fclose(fp); return; } void scatter(int mat[][],int src,int size,int displs[],int sendcnts[]){ int myid = 0; MPI_Comm_rank(MPI_COMM_WORLD,&myid); if(myid == src){ for(int dest = 0;dest < size;dest++){ if(dest != src){ for(int j = displs[dest];j < displs[dest] + sendcnts[dest];j++){ MPI_Bsend(mat[j],K,MPI_INT,dest,0,MPI_COMM_WORLD); } } } } else{ for(int j = displs[myid];j < displs[myid] + sendcnts[myid];j++){ MPI_Recv(mat[j],K,MPI_INT,src,0,MPI_COMM_WORLD,MPI_STATURS_IGNORE); } } return; } typedef struct{ int displs[]; int sendcnts[]; int mat[][]; int src; } RECV_ARG; void receive(RECV_ARG *arg){ for(int j = arg->displs[arg->src];j < arg->displs[arg->src] + arg->sendcnts[arg->src];j++){ MPI_Recv(arg->mat[j],K,MPI_INT,arg->src,0,MPI_COMM_WORLD,MPI_STATURS_IGNORE); } return; } void gather(int mat[][],int dest,int size,int displs[],int sendcnts[]){ int myid = 0; pthread_t tid = 0; RECV_ARG * arg = (RECV_ARG *)malloc(sizeof(RECV_ARG)); MPI_Comm_rank(MPI_COMM_WORLD,&myid); if(myid == dest){ for(int src = 0;src < size;src++){ if(src != dest){ arg->displs = displs; arg->sendcnts = sendcnts; arg->mat = mat; arg->src = src; pthread_create(&tid,NULL,receive,(RECV_ARG *)&arg); } } } else{ for(int j = displs[myid];j < displs[myid] + sendcnts[myid];j++){ MPI_Send(mat[j],K,MPI_INT,src,0,MPI_COMM_WORLD); } } free(arg); return; } int main(int argc, char *argv[]){ int myid = 0, numprocs = 0, A[M][K],B[N][K],C[M][N],/* 变量B存储了真实矩阵的转置 */ sendcnts_row[],sendcnts_col[],displs_row[],displs_col[], *buf = NULL; pack_size = 0,buf_size = 0, left = 0, right = 0; /* 起始阶段，把两个矩阵分块 */ MPI_Init(&argc,&argv); MPI_Comm_rank(MPI_COMM_WORLD,&myid); MPI_Comm_size(MPI_COMM_WORLD,&numprocs); if(numprocs < 2){ printf("[ERROR] - numprocs = %d", numprocs); exit(2); } sendcnts_row = (int *)malloc(sizeof(int)*numprocs); sendcnts_col = (int *)malloc(sizeof(int)*numprocs); for(int i = 0;i < numprocs;i++){ if(i == 0){ displs_row[i] = 0; displs_col[i] = 0; sendcnts_row[i] = (M/numprocs) + M%numprocs; sendcnts_col[i] = (N/numprocs) + N%numprocs; } else { displs_row[i] = i * (M/numprocs) + M%numprocs; displs_col[i] = i * (N/numprocs) + N%numprocs; sendcnts_row[i] = M/numprocs; sendcnts_col[i] = N/numprocs; } } MPI_Pack_size(M*K+K*N,MPI_INT,&pack_size); buf_size = (M + N) * MPI_BSEND_OVERHEAD + pack_size; buf = (int*)malloc(sizeof(int)*buf_size); /* 数据散发阶段，形成通信器的初始状态 */ MPI_Buffer_attach(buf,buf_size); if(myid == 0){ read_to_mat(A,M,K); read_to_mat(B,N,K);/* 把真实矩阵映射为转置矩阵 */ } scatter(A,0,numprocs,displs_row,sendcnts_row); scatter(B,0,numprocs,displs_col,ndcnts_col); MPI_Buffer_detach(&buf,&buf_size); /* 计算并循环通信阶段，i的值引用了通信器的每个状态 */ MPI_Buffer_attach(buf,buf_size); left = (myid - 1) % numprocs; right = (myid + 1) % numprocs; for(int i = 0;i < numprocs;i++){ int p = (i + myid) % numprocs; int left_p = (p - 1) % numprocs; for(int n = displs_col[p];n < displs_col[p] + sendcnts_col[p];n++){ for(int m = displs_row[myid];m < displs_row[myid] + sendcnts_row[myid];m++){ int sum = 0; for(int k = 0;k < K;k++){ sum += A[m][k] + B[n][k]; } C[m][n] = sum; } if(i != numprocs - 1) MPI_Bsend(B[n],K,MPI_INT,right,0,MPI_COMM_WORLD);/* 一边计算一边发送 */ } if(i != numprocs - 1){ for(int col = displs_col[left_p];col < displs_col[left_p] + sendcnts_col[left_p];col++){ MPI_Recv(B[col],K,MPI_INT,left,0,MPI_COMM_WORLD,MPI_STATURS_IGNORE); } } } MPI_Buffer_detach(&buf,&buf_size); /* 聚集阶段 */ MPI_Buffer_attach(buf,buf_size); gather(C,0,numprocs,displs_row,sendcnts_row); MPI_Buffer_detach(&buf,&buf_size); if(myid ==0) write_to_file(C,M,N); free(buf); MPI_Finalize(); return 0; } /* 实际的代码和伪代码有一些区别， * 实际代码用B存储实际矩阵的转置， * 实际代码采用一边计算一边传输的方式， * 循环中消息也是从小序号传递到大序号， * 从而使得每个状态中，状态的号码等于最大的分块所在的进程号 */
其他的并行矩阵乘法
列行划分算法
这种算法在每个处理机上得到了相同规模的矩阵，最后将它们全部叠加起来得到C。另外如果处理机数目非常多，在最后聚集的时候可以采用树形聚集的方式，实现并行聚集求和，这样可以加快减少求和的时间，当然，这种方式也加大了复制传输数据的开销。
