94 MIMD共享存储模型的并行算法

9.4MIMD共享存储模型的并行算法(一)异步枚举排序算法枚举排序(EnumerationSort)是一种最简单的排序算法，通常也称为秩排序(RankSort)。对于基于枚举比较原理的异步排序算法，为了排序n个数的序列S=(x1,x2，...,xn)，算法要生成n个进程。进程i(IWiWn)将xi与S中其余元素进行比较，并且使用局部变量k及下所有小于xi的元素的数目。当所有的比较都完成时，就将xi置入排序序列中的k+1的位置上，，因此每个进程都可能彼此独立地执行，而无通信的要求。令X是存在共享存储器中长度为n的数组，开始时放入被排序的序列，当算法结束时，结果置于共享存储器中的T数组内，变量i,j,k是算法生成的每个进程的局部变量［1。。7］。该并行算法的描述如下：算法9.4.1MIMD-CREW模型上的异步枚举排序算法输入：S=(x1,x2,.,xn)置于共享存储器的X数组中输出：已排序的序列置于共享存储器的T数组中ENUERSORTTOC\o"1-5"\h\z{for(i=1;i<=n;i++){CreateProcessi;}Processi:k=0;for(j=1;j<=n;j++){if(X(i)>X(j)){k^k+1;}elseif((X(i)==X(j)&&i>j)){k^k+1;}}T(k+1)=X(i);} 在该并行算法中有，由于每个处理器至多确定数组X中的「n/p］个元素的位置，而每确定一个元素的位置需O(n)时间，因此对X进行排序需「n/p】XO(n)时间。例9_8设S=(8,6,6,9,7)，p(n)=2。创建进程的过程中，假定由P1生成5个进程，此后所有的进程均等待开始。假定使用“先进先出”的调度策略。进程1和2分别由P1和P2执行，它们同时计算元素8和6的次第，然后将其分别置于T数组的相应位置上，如图9_19(a)所示。欲知下一进程有哪一个处理器启动执行，需研究进程1和2的执行时间。假定比较操作和赋值操作大致用相同的时间，当执行

X(i)>X(j),X(i)=X(j)和i>j以及k=0,k=k+1和T(k+1)=X(i)时，进程1和2分别需要14和17各时间步，如图9_19(b)所示。所以进程3应由P1启动执行，之后3个时间步P2启动执行进程4。这两个进程负责将元素6和9分别置于数组T的相应位置上，此时进程3需要18个时间步，进程4需要13个时间步，所以进程4比进程3早结束。下一进程5由P2启动。之后当进程3执行完后，因无等待的进程，所以P1变为空闲。当进程5执行15个时间步后，元素7便放到数组T的相应位置上。算法总共需要45个时间步。Pi进程1 进程Pi进程1 进程3□ 14 32P2进程2讲程4进程50 17 30 45(b)进度调度(a)数组T(b)进度调度图9_19异步枚举排序数组变化及进程调度(二)单源点最短路径并行算法假设一有向图各边赋于非负整数，某条路径的长度就是沿该路径所有边权之和。最短路径问题，就是对每一点对i和j，丘照期间任何最短长度的路径。单源点最短路径，即在一个图中寻找从一个指定顶点到所有其他定点的最短路径。(a)单处理机上的Moore算法。在Moore算法中，设源点为s^V。从s到其他各顶点的最短路径长度用一维数组dist存储。首先置dist(s)=0,dist(v)=8，其中v*，且vEV。算法使用了一个队列queue。开始执行时队列中仅含有源点s；以后每次只要队列非空，就将排头顶点u从队列中移去，令其作为本次搜索的顶点，然后检查u的所有射出边(u,v)EE：若dist(u)+w(u,v)<dist(v)，则此时就找到了一条从s到v的更短路径，置dist(v)=dist(u)+w(u,v)；若v不再搜索队列中，则把v加到队尾。如此重复进行，直到整个待搜索队列空时，算法就终止。算法9・4・2_1SISD机器上的单源点最短路径算法输入：有向加权图G(V,E)的权矩阵输出：从源s到其余顶点i(i/s)的最短路径dist(i),iEVSISDSHORTESTPATH{(1)for(i=1;i<=n;i++)/*初始化*/3{INITIALIZE(i);}(2)将s插入队列queue中；(3)while(队列非空){SEARCH;}8}StatusSEARCH{(3.1)dequeuevertexu;/*从队列中删去顶点u*/(3.2)for(eachedge(u,v)EE)TOC\o"1-5"\h\z{new_dist6dist(u)+w(u,v);if(new_dist<dist(v)){dist(v)6new_dist;}if(visnotinthequeue){enqueuevertexv;}}}StatusINITIALIZE{(1.1)dist(s)60;(1.2)dist(v)。8;(v#s,vEV)(1.3)queue60;}(b)Moore算法的并行化Deo等人基于MIMD紧耦合共享存储模型实现了Moore算法的并行化。直观上讲，算法9.4.2_1有两处可并行化的地方：一是SEARCH的第(3.2)步；二是主算法的第(3)步。前者，任何一个顶点均可能有多条射出边，它们都可并行地被检查；后者，在任何时候都可能有多个顶点在队列中，因此有可能每次检查多个顶点的射出边。由于后者可产生较大的加速，而且当G是个稀疏图时，并行度受定点射出变得影响，所以选用后者。首先，队列用源点初始化。然后创建了许多异步进程，每个进程都从队列中删除一个顶点，检查其射出边，将已发现有更短路径的顶点加入到队列。算法的第(1)步采用预调度方法很容易并行化。而第(3)步的while循环需作适当的修改，以能反映并执行SEARCH过程时一些异步进程的存在。显然，当一个进程发现队列为空时，就停止执行是不合适的。因而必须采用两个变量联合使用的办法，以决定什么时候无工作可做：第一个是数组变量waiting，它记住哪一个进程正处于等待状态；第二个是布尔变量halt，仅当队列为空和所有进程处于等待状态时为真。INITIALIZE过程置数组waiting中的第一个元素为假。SEARCH过程也必须作适当的修改。因为队列的插入、删除操作不是原子操作，所以执行上述操作时必须给队列上锁：其次，在一个进程将刚找到的v路径new_dist与当前的最短路径dist(v)比较之前，变量dist(v)也必须上锁，否则两个进程有可能同时修改它；最后，若一个进程发现队列为空时，则Swaiting中的相应元素为真。若进程1处于等待状态，则它要检查每个进程是否都处在等待状态，如果是，则halt置为真，而进程1检查每个进程是否处于等待状态时，队列也必须上锁。该并行算法的描述如下：算法9.4.2MIMD-SM模型上单源点最短路径算法

1234567891011121314151617181920212223242526272829303132333435363738394041输入：有向加权图G的权矩阵W输出：从源s到其余顶点i(iNs)的最短路径dist(i),iEVMIMDSHORTESTPATH{parfor(i=1;i<=p;i++)/*p为进程数*/{for(j=i;j<=n;j=j+p)/*初始化*/{INITIALIZE(j);}}enqueues;/*s入队*/halt=FALSE;parfor(i=1;i<=p;i++){while(nothalt){SEARCH(i);}}}StatusSEARCH(i){lockthequeue;/*队列上锁*/if(queueisempty)/*队列空时，等待进程为真*/{waiting(i)6TRUE;if(i==1)/*进程1等待时，其它进程均须等待*/{halt=waiting(2)八waiting(3)八•••八waiting(p);}unlockthequeue;/*队列开锁*/}else{dequeueu;/*从队列中删除u*/waiting(i)6FALSE;unlockthequeue;/*队列开锁*/for(everyedge(u,v)inGraph)/*检查每条射出边*/{new_dist。dist(u)+w(u,v);lock(dist(v))/*dist(v)上锁*/if(new_dist<dist(v)){

42dist(v)6new_dist;/*更新new_dist*/43unlock(dist(v));/*dist(v)开锁*/44if(v史queue)45{46lockthequeue;/*队列上锁*/47enqueuev;/*v入队*/48unlockthequeue;/*队列开锁*/49}50}51else{unlock(dist(v));}/*dist(v)开锁*/52}53 }54}(二)最小生成树并行算法一个无向连通图G的生成树是指满足如下条件的G的子图T：G和T具有相同的顶点数；在T中有足够的边能连接G所有的顶点，但不出现回路。最小生成树，即最小代价生成树(MinimumCostSpanningTree)，它是加权连通无向图的所有生成树中权值最小的生成树。(a)Sollin顺序求MST算法算法开始时，图的n个孤立顶点视为一片森林，而每个顶点均视为一棵树；算法共迭代logn次，每次迭代时，森林中的每棵树，同时决定其最小的邻边，并将它们加入到森林中(即合并各棵树)；此过程重复到森林中只剩下一棵树。因为森林中树的数目，每次都以2的倍数减少，所以迭代至多需要hogn】次就可找到MST；而每次迭代时，找顶点的最小邻边至多执行O(n2)次比较。因此Sollin的顺序算法的复杂度为O(n2logn)。函数FIND(v)找顶点v所在的树的名字，UNION(v,u)合并包括顶点v和u的两棵树。算法9・4・3_1SISD机器上的SollinMST算法输入：无向图G的加权矩阵W输出：G的最小生成树T(树以边的形式存储)SOLLIN-MSTTOC\o"1-5"\h\z{(1)for(i=1;i<=n;i++)/*初始化*/{vertexiisinitiallyinseti;}(2)T60;/*T就是MST*/

7(3)while(|T|<n-1)8{9(3.1)for(eachtreei)/*每棵树i找不在同一树中的最小边权*/10{11closest(i)<8;12}13(3.2)for(each(v,u)EE)14{15if(FIND(v)^FIND(u))/*v和u属于不同的连通片*/16{17if(w(v,u)<closest(FIND(v)))18{19closest(FIND(v))6w(v,u);20edge(FIND(v))<(v,u);21}22}23}24(3.3)for(eachtreei)25{26(3.3.1)(v,u)^edge(i);27(3.3.2)if(FIND(v)丰FIND(u))28{29(3.3.2.1)T6TU{(v,u)};/*加入新的树边*/30(3.3.2.2)UNION(v,u);/*合并成较大的树*/31}32}33}34}(b)Quinn的并行化算法在共享存储的MIMD机系统上，Sollin算法可按如下思路并行化：首先，应设法对第(3)步while循环进行并行化，但遗憾的是因循环之间存在着先后制约关系，即在第k+1次循环之前，第k次循环的子树必须同一个有最小权值的边关联的另一棵子树归并，所以while语句的并行化是有限的。其次，考虑循环体内的并行化。算法的第(3.1)步通过适当的预调度可以使其并行化，设第t次循环时已有nt棵子树:若nt>p，则把nt棵子树较均匀的分配到p个处理器中，每个处理器约有［nt/p〕棵子树；否则，这nt棵子树就分配给nt个处理器。算法的第(3.2)步并行化的最有效做法是：首先，每个处理器检查它内部的顶点的边，然后再检查不在同一处理器上树之间的边。算法的第(3.3)步并行化稍微复杂些。图9_20示出了此情况，假定一个处理器企图将树B与其最近的树A进行归并，变量edge(A)有一条权值为k的边(vA,uA)，变量edge(B)也有一条权值为k的边(vB,uB)。然而假定两处理器都在实行UNION之前执行了第(ii)步的测试，则(vA,uA)和(vB,uB)都将加入到树T中，结果形成一个环，显然这是错的。因此，欲使第3.3)步并行化，必须在执行第(3.3.2)步之前上锁，执行完第(3.3.2.2)步后再开锁。每次仅允许一个处理器进入临界区。为了避免死锁的产生，当有多个请求上锁的进程申请临界区时，仅仅让标号最小的子树上锁。在MIMD-SM模型上，所法所需的时间为~匚 l.n算法9・4・4MIMD-CREW模型上的算法9・4・4MIMD-CREW模型上的Gauss-Seidel算法输入：AnXn，向量b，初值x0(i=1,...,n)，精度c(0<c<1)输出：满足精度要求的x的近似解GUASS-SEIDEL { for(i=1;i<=n;i++)3{O(logn(一+—+n+p))，处理器数为O(p)[ioo7]。PP树A 树B图9_20并行化Sollin算法所引起的复杂情况i.Gauss-Seidel算法Gauss-Seidel算法是求解Ax=b的一种算法。先将系数矩阵A分解为A=E+D+F，其中D、E和F均为nXn的矩阵，它们的元素分别定义如下：d」％,eW卜jf<jij[0,否则，j[0,否则j[0,否则这样，Ax=(D+E+F)x=b，从而Dx=b-Ex-Fx。给定初始值x0，则第k次迭代为Dxk=b-Exk-Fxk-1对于某一k和给定的误差允许值c，如果下式成立，则认为迭代时收敛的。/|xk+1一Xk|<Ci=1Gauss-Seidel异步并行算法，设有N个处理器(NWn)，算法生成n个进程，每一个进程计算x的一个分量。这些进程由N个处理器异步的执行。令x0代表初值，Iold；代表旧值，new；代表当前值，c为精度[1007]。该并行算法的描述如下：old£xo；7CreateProcessi;/*生成i个进程*/8}9Processi:10do11{old.<new.;ii12睥可£尝归外"d)f(⑶k=1 k=i+113}while(!(81new-old.\<c))14=115xi^newi;}6new£气。;xoldk))«例9_9使用算法9.4.4求解下述方程：|2x1+x2=3"1+2x2=4设有2个处理器(N=2)，取初值x0=1/2,x0=3/2,c=0.02。设置old]=1/2,neW]=1/2，生成进程processl；设置old2=3/2,new2=3/2，生成进程process2。process1：设置old1=1/2；计算new1=(1/2)(3-3/2)=3/4；process2：设置old2=3/2；计算new2=(1/2)(4-1/2)=7/4o重复计算如下：new1=5/8, new2=13/8；new1=11/16,new2=27/16；new1=21/32,new2=53/32；new1=43/64,new2=107/64；new1=85/128,new2=213/128o因为|43/64-85/128|+|107/64-213/128|<0.02,所以迭代结束。牛顿求根并行算法牛顿求根法，设f(x)在(a,b)上连续，且f’(x)N0,f”(x)N0,f(a)•f(b)<0,则方程f(x)=0的根z的近似值可用如下迭代公式计算，直到误差|xn+1-xn|<c，其中c为给定的精度： nxn+1=xn-f(xn)/f’(xn)并取初始值如下：枷,当广3)<0X=< 0lb,当广⑴>0牛顿方法的几何解释如图9_21所示(假定f”(x)>0)。下一次迭代值xn+1就是曲线f(x)在xn处的切线与x轴的交点。MIMD机器上的牛顿求根算法，将包含f(x)的零点z的区间(a,b)(设a<b)划分成N+1等分(NN2)。各分点取为z的近似值。计算由N个进程组成，每一进程开始用各分点值进行牛顿法求根。各进程并发异步的执行。一旦某一进程收敛，它就向共享存储单元ROOT写它所求得的值。开始时根ROOT被置成8,一旦ROOT值被某一进程修改，所有进程将随之结束。如果两个进程同时收敛，就会造成写冲突，于是按最小号码的进程优先的策略来裁决，其它进程都被拒绝。如果在预定的迭代数目r之内某一进程不收敛，则它就被挂起［1007］。该并行算法的描述如下: 算法9.4.5MIMD-CREW模型上的牛顿求根算法输入：f(x)，f’(x)，(a,b)，精度c，最大容许迭代次数r输出：返回答案在ROOT中ROOTs((b-a)/(N+1);for(k=1;k<=N;k++){createprocessk;/*生成进程*/}ROOT。8;Processk:xold(a+ks;iteration。0;while(iteration<r&&ROOT==8){iteration++;xnew。xoldfgl/政1xnew-xold3{17ROOT。x;new18}19x,,6x ;old new20}21}在牛顿求根并行算法中，令N为可用处理器数，如果N足够大，则有一个起始点可充分接近z。如果f(x),f’(x),f”(x)在区间(a,b)内连续有界，那么其中有一个处理器将会在O(logm)时间内收敛，其中m为所希望的精度的位数。例9_10令f(x)=x3-4x-5。因此f’(x)=3x2-4。在区间(-3,3)存在着f(x)的零点，令N=5,区间被6等分，各分点的值依次为-2,-1,0,1,2。算法生成5个进程，令e=9-10，假定5个处理器同时执行5个进程。在此情况下，进程5最先收敛到根，其值ROOT=2.456678。b)MIMD异步通信模型的并行算法i.快速排序并行算法快速排序(QuickSort)是一种最基本的排序算法，它的基本思想是，在当前无序序列R[1,n]中取一个记录作为比较的基准，用此基准将当前的无序序列R[1,n]划分成两个无序序列R[1,i-1]和R[i,n](1WiWn)，且R[1,i-1]中记录的所有关键字均小于等于基准的关键字，R[i,n]中记录的所有关键字均大于等于基准的关键字；当R[1,i-1]和R[i,n]非空时，分别对他们重复上述的划分过程。对每次划分过后所得到的两个序列分别使用两个处理器完成递归排序。例如对一个长为n的序列，首先划分得到两个长为n/2的序列，将其交给两个处理器分别处理；而后进一步划分得到4个长为n/4的序列，在分别交给4个处理器处理；如此递归下去最终得到排序号的序列[1001]。该并行算法的描述如下: 算法9・5・1快速排序并行算法输入：无序数组data[1,n]，使用的处理器个数2m输出：有序数组data[1,n]SORT{para_quicksort(data,1,n,m,0);}para_quicksort(data,i,j,m,id)4{5 if((j-i)Wk||m==0)6{Pid:quicksort(data,i,j);8}else{Pid:r=partition(data,i,j);Pidsenddata[r+q,m-1]toPid+2m-1para_quicksort(data,I,r-1,m-1,id);paraquicksort(data,r+1,j,m-1,id+2m-1);P，d+2m-Senddata[r+1,m-1]toP，d；} 算法9・5・算法9・5・2矩阵转置并行算法输入：矩阵A*nxn输出：矩阵Anxn的转置矩阵ATnxnTRANSPOSEDMATRIX {/*对所有处理器my_rank(my_rank=0,...,p-1)*/v(my_rank/sqrt(p);/*计算子块的行号*/u(my_rankmodsqrt(p);/*计算子块的列号*/if(u<v)/*对存放下三角块的处理器*/6{send所存的子块to其对角块所在的处理器；receive其对角块所在的处理器中发来的子块；}}在最优的情况下该并行算法形成一个高度为logn的排序树，其计算复杂度为O(n),通信复杂度为O(n)；在最坏情况下其计算复杂度降为O(n2)；正常情况下该算法的计算复杂度为O(n)。ii. 二维网孔上的矩阵转置并行算法网孔上的矩阵转置并行算法思路是，实现矩阵转置时，若处理器个数为p,且它们的编号依次是0,1,...,p-1,则将n阶矩阵A分成p个大小为mXm的子块，m=「n/p〕。p个字块组成一个.£FXv-7的子块阵列。记其中第i行第j列的子块为Aj它含有A的第(i-1)m+1至第im行中的第(j-1)m+1至第jm列的所有元素。对每一处理器按行主方式赋以二维下标，记编号为i的处理器的而为下标为(v,u)，其中v=£jp」u=imodjp，将A的子块存入下标为(v,u)表示的对应处理器中。这样，转置过程分两步进行：第一步，子块转置，具体过程如图9_22所示；第二步，处理器内部局部转置。图9_22子块转置为了避免对应子块交换数据时处理器发生死锁，可令下三角子块先向与之对应的上三角子块发送数据，然后从上三角子块接收数据；上三角子块先将数据存放在缓冲区buffer中，然后从与之对应的下三角子块接收数据；最后再将缓冲区中的数据发送给下三角子块[1001]。该并行算法的描述如下：else/*对存放上三角块的处理器*/TOC\o"1-5"\h\z{将所存的子块在缓冲区buffer中做备份；receive其对角块所在的处理器中发来的子块；sendbuffer中所存子块to其对角块所在的处理器；}for(i=1;i<=m;i++) /*处理器内部局部转置*/{for(j=1;j<=i;j++){交换a[i,j]和a[j,i];}}}若记ts为发送启动时间，tw为单位数据传输时间，th为处理器间的延迟时间，则第一步由于每个子块有n2/p"个元素，又由于通信过程中为了避免死锁，错开下三角子块与上三角子块的发送顺序，因此子块的交换时间为n2 n22-+2tA；P+n2 n22-+2tA；P+2八丁+《\：p。2p p部转置时间为n2/2p。因此所需的并行计算时间7p=矩阵并行分块乘法算法矩阵并行分块乘法的基本思想：将矩阵A按行划分为p块(p为处理器个数)，设u=「m/p】，每块含有连续的u行向量，这些行块依次记为A。,％...▲-]，分别存放在标号为0,1,...,p-1的处理器中。对矩阵B按列划分为p块，记v=k/p，每块含有连续的v列向量，这些列块依次记为B0,B],...,Bp1，分别存放在标号为0,1,...,p-1的处理器中。将结果矩阵C也相应的同时进行行、列划分，得到pXp个大小为uXv的子矩阵，记第i行第j列的子矩阵为C”则有C^AiXBj，其中Ai大小为uXn，Bj大小为nXv。开始，各处理器的存储内容如图9_23(a)所示。此时各处理器并行计算C厂A|XB|，其中i,j=0,1,...,p-1，此后第i号处理器将所存储的B的列块送至第i-1号处理器(第0号处理器将B的列块送至第p-1号处理器中，形成循环传送)，各处理器中存储内容如图9_23(b)所示。它们再次并行计算Cij.=AiXBj，这里j=(i+1)modp。B的列块在处理器中以这样的方式循环传送p-1次，并做p次子矩阵相乘运算，就生成了矩阵C的所有子矩阵。编号为i的处理器的内部存储器存有子矩阵Ci0,Ci1,.,Ci(p1)O为了避免在通信过程中发生死锁，奇数号及偶数号处理器的收发顺序被错开，使偶数号处理器先发送后接收；而奇数号处理器先将B的列块存于缓冲区buffer中，然后接受编号在其后面的处理器所发送的B的列块，最后再将缓冲区的原矩阵B的列块发送给编号在其前面的处理器"01】。该并行算法的描述如下：——处理器编号 存储器内容处理器编号 存储器内容121F12345678910111213141516171819202122232425262728(a) (b)图9_23矩阵相乘并行算法中的数据交换算法9・5・3矩阵并行分块乘法算法输入：矩阵AZBnXk输出：矩阵C灯mXkBLOCK-MATRIXMUTIPLICATION{/*对所有处理器my_rank(my_rank=0,…,p-1)*/K(i+my_rank)%p;/*计算C的子块号*/for(z=0;z<=u-1;z++){c[l,z,j](0；for(j=0;j<=v-1;j++){c[l,z,j](c[l,z,j]+a[z,s]*b[s,j];}}mm1((p+my_rank-1)%p;/*计算左邻处理器的标号*/mm1((my_rank+1)%p;/*计算右邻处理器的标号*/if(iJp-1){if(my_rank%2==0)/*编号为偶数的处理器*/{将所存的B的子块发送到其左邻处理器中接收其右邻处理器中发来的B的子块}if(my_rank%2/0)/*编号为奇数的处理器*/{将所存的B子块在缓冲区buffer中做备份；receive其右邻处理器中发来的B的子块；sendbuffer中所存的B的子块to其左邻处理器;}}}矩阵并行分块乘法中，设一次乘法和加法运算时间为一个单位时间，由于每个处理器计算p个uXn与nXv阶的子矩阵相乘，因此计算时间为uXnXvXp；所有处理器交换数据p-1次，每次的通信量为vXn，通信过程中为了避免死锁，错开奇数号及偶数号处理器的收发顺序，通信时间为2(p-1)(ts+nvtw)=O(nk)，所以并

行计算时间T=uvnp+2(p-1)(t+nvt)=mnk/p+2(p-1)(t+nvt)。PJ. '•A /'e wt^ A '•A /'e wt^sw sw分布式矩阵求逆的并行算法矩阵求逆的并行算法的基本思想：在矩阵求逆的过程中，一次利用主行i(i=0,1,...,n-1)对其余各行j(j^i)作初等变换，由于各行计算之间没有数据相关关系，因此对矩阵A按行划分来实现并行计算。考虑到在计算过程中处理器之间的负载均衡，对A采用行交叉划分:设处理器个数为p,矩阵A的阶数为n,m=「n/p】，对矩阵A行交叉划分后，编号为i(i=0,1,...,p-1)的处理器存有A的第i,i+p,...,i+(m-1)p行。在计算中，依次将第0,1,...,n-1行作为主行，将其广播给所有处理器，这实际上是各处理器轮流选出主行并广播。发送主行数据的处理器利用主行对其主行之外的m-1行行向量做行变换，其余处理器则利用主行对其m行行向量做行变换［1001］。该并行算法的描述如下：算法9・5・4矩阵求逆并行算法输入：矩阵A*nxn输出：矩阵A-1nxnINVERSEMATRIX1{2/*对所有处理器my_rank(my_rank=0,...,p-1)*/3for(i=0;i<=m-1;i++)4{5for(j=0;j<=p-1;j++)6{7if(my_rank=j)/*主兀素在本处理器中*/89{v6i*p+j;/*v为王兀素的行号*/10a[i,v]61/a[i,v];11for(k=0;k<=n-1;k++)12{13if(k^v)14{15a[i,k]6a[i,k]*a[i,v];16}17}18for(k=0;k<=n-1;k++)19{20f[k]6a[i,k];21}2223broadcast变换后的主彳丁兀素(存于数组f中)to所有处理器:■ 3、・-1-7-^ J..LI-rmRDr-l-r“/else/*主兀素不在本处理器中*/2425{v6i*p+j;/*v为主兀素行号*/26reveive主彳丁兀素存于数组中；27}

282930313233343536373839404142434445464748495051525354555657585960616263646566676869 }70 }if(my_rank勺)/*主行元素不在本处理器中*/{for(k=0;k<=m-1;k++)/*处理非主行、非主列元素*/{for(w=0;w<=n-1;w++){if(w^v){a[k,w]<a[k,w]-f[w]*a[k,v];}}}for(k=0;k<=m-1;k++)/*处理主行元素*/{a[k,v]<-f[v]*a[k,v];}}else/*处理主行所在处理器中的其它元素*/{for(k=0;k<=m-1;k++){if(k^i){for(w=0;w<=n-1;w++){if(w^v){a[k,w]<a[k,w]-f[w]*a[k,v];}}}}for(k=0;k<=m-1;k++){if(k^i){a[k,v]6-f[v]*a[k,v];}}}}矩阵求逆的并行算法中，若一次乘法和加法运算或一次除法运算时间为一个单

位时间，则所需的计算时间为mn2；又由于共有n行数据依次作为主行被广播，其通信时间为n(t+nt)logp，所以该算法并行计算时间为T=mn2+n(t+nt)logp。SW p sw分布式高斯消去并行算法高斯消去法是利用主行i对其余各行j(j>i)，作初等变换，各行计算之间没有数据相关关系，因此可以对矩阵A按行划分。考虑到在计算过程中处理器之间的负载均衡，对A采用行交叉划分。设处理器个数为p，矩阵A的阶数为n,m=「n/p】，对矩阵A行交叉划分后，编号为i(i=0,1,...,p-1)的处理器含有A的第i,i+p,...,i+(m-1)p行和向量B的第i,i+p,...,i+(m-1)p一共m个元素。消去过程的并行是依次以第0,1,...n-1行为主行进行消去计算，由于对行的交叉划分与分布，这实际上是由各处理器轮流选出主行。在每次消去计算前，各处理器并行求其局部存储器中右下角子阵的最大元。若以编号为my_rank的处理器的第i行作为主行，则编号在my_rank后面的处理器(包括my_rank本身)求其局部存储器中第i行至第m-1行元素的最大元，并记录其行号、列号及所在处理器编号；编号在my_rank前面的处理器求其局部存储器中第i+1行至第m-1行元素的最大元，并记录其行号、列号及所在处理器编号。然后通过扩展收集操作，将局部存储器中最大元按处理器编号连接起来，并广播给所有处理器，各处理器以此求得整个矩阵右下角阶子阵的最大元maxvalue及其所在行号、列号和处理器编号。若maxvalue的列号不是原主元素akk的列号，则交换第k列与maxvalue所在列的两列数据；若maxvalue的处理器编号不是原主元素a^的处理器编号，则在处理期间的进行行交换;若maxvalue的处理器编号是原主元素a^的处理器编号，但行号不是原主元素a^的行号，则在处理其内部进行行交换。在消去计算中，首先对主行元素作归一化 操作akj=akj/akk，然后将主行广播给所有处理器，各处理器利用接收到的主行元素对其部分行向量做行变换。若以编号为my_rank的处理器的第i行作为主行，并将它播送给所有的处理器。则编号在my_rank前面的处理器(包括my_rank)本身利用主行对其第i+1,...,m-1行数据和子向量B做行变换。编号在my_rank后面的处理器利用主行对其第i,...,m-1行数据和子向量B做行变换。回代过程的并行是按xn,xn-1,.x1的顺序由各处理器依次计算x(i*p+myrank)，一旦x(i*p+myrank)被计算出来就立即广播给所有处理器，用于与对应项相乘并作求和计算I*：。该并行算法的描述如下：算法9・5・5全主元高斯消去法过程的并行算法输入：系数矩阵Anxn，常数向量bnx1输出：解向量xnx1GAUSSELIMINATIONTOC\o"1-5"\h\z{/*对所有处理器my_rank(my_rank=0,...,p-1)*//*消去过程*/for(i=0;i<=m-1;i++){for(j=0;j<=p-1;j++){if(my_rank<j)/*对于主行后面的块*/{vi*p+j;/*v为主元素的行号*/

1112131415161718192021222324252627282930313233343536373839404142434445464748495051525354/*确定本处理器所存未消元部分的最大元及其位置，存于数组lmax中*/lmax[0]6a[i+1,v];for(k=i+1;k<=m-1;k++){for(t=v;t<=n-1;t++){if(|a[k,t]|>lmax[0]){lmax[0]6a[k,t];lmax[1]6k;lmax⑵6t;lmax[3]6my_rank;}}}}if(my_rankNj)/*对于主行后面的块*/{v6i*p+j;/*v为主元素的行号*//*确定本处理器所存未消元部分的最大元及其位置，存于数组lmax中*/lmax[0]6a[i,v];for(k=i;k<=m-1;k++){for(t=v;t<=n-1;t++){if(|a[k,t]|>lmax[0]){lmax[0]6a[k,t];lmax[1]6k;lmax⑵6t;lmax[3]6my_ran;}}}}broadcast处理器中的lmax元素to其他所有处理器；/*确定最大元及其位置*/maxvalue6getpivort(max);maxrow6getrow(max);maxcolumn6getcolumn(max);maxrank6getrank(max);/*列交换*/

if(maxcolumn#v){for(t=0;t<=m;t++){交换a[t,v]与a[t,maxcolumn];}交换shift[v]与shift[maxcolumn];}/*行交换*/if(my_rank=j){if(maxcolumn#a[i,v]){if((maxrank==j)&&(i#maxrow)){innerexchangerow();/*处理器内部换行*/}if(maxrank#j){outerexchangerow();/*处理器之间换行*/}}}if(maxrank#j)/*主行所在的处理器*//*对主行作归一化运算*/{for(k=v+1;k<=n-1;k++){a[i,k]6a[i,k]/a[i,v];}b[i]6b[i]/a[i,v];a[i,v]J;/*将主行元素赋于数组f*/for(k=v+1;k<=n-1;k++){f[k]6a[i,k];}f[n]6b[i];broadcast主行to所有处理器；}else/*非主行所在的处理器*/{receive主行元素存于数组f中；5556575859606162636465666768697071727374757677787980818283848586878889909192939495969798}99100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142if(my_rankWj){for(k=i+1;k<=m-1;k++){for(w=v+1;w<=n-1;w++){a[k,w]<a[k,w]-f[w]*a[k,v];}b[k]<b[k]-f[n]*a[k,v];}}if(my_rank>j){f