MPI 并行解上三角矩阵表示的方程组

MPI 并行解上三角矩阵表示的方程组

MPI最基本的回代法实现

1. #include "mpi.h"   
2. #include <stdio.h>   
3. #include <stdlib.h>   
4.   
5. #define ROOT 0   
6. #define TAG 0   
7. int main(int argc,char *argv[]) {  
8.     float matrix[][4]={{2,3,4,5},//the upper trangular matrix representing a equation set   
9.                 {0,2,3,4},  
10.                 {0,0,2,3}};  
11.     int xcnt=3;//cout of x   
12.     int self,size,tag=0;  
13.     MPI_Init(&argc,&argv);  
14.     MPI_Comm_rank(MPI_COMM_WORLD,&self);  
15.     MPI_Comm_size(MPI_COMM_WORLD,&size);  
16.   
17.     MPI_Request r;  
18.     MPI_Status s;  
19.   
20.     MPI_Datatype MPI_VEC;  
21.     MPI_Type_vector(xcnt+1,1,1,MPI_FLOAT,&MPI_VEC);  
22.     MPI_Type_commit(&MPI_VEC);  
23.     float* equation=(float*)malloc((xcnt+1)*sizeof(float));  
24.     float xs,xr;  
25.     if(0==self) {//send each process the correspond equation   
26.         //parellel send and copy   
27.         MPI_Issend(matrix[1],1,MPI_VEC,1,TAG,MPI_COMM_WORLD,&r);  
28.         for(int i=0;i<=xcnt;++i) {  
29.             equation[i]=matrix[0][i];  
30.         }  
31.         MPI_Wait(&r,&s);  
32.         for(int i=2;i<size;++i) {  
33.             MPI_Ssend(matrix[i],1,MPI_VEC,i,TAG,MPI_COMM_WORLD);  
34.         }  
35.     } else {  
36.         MPI_Recv(equation,1,MPI_VEC,ROOT,TAG,MPI_COMM_WORLD,&s);  
37.     }  
38.       
39.     for(int i=xcnt-1;i>=0;--i) {  
40.         if(i==self) {  
41.             xs=equation[xcnt]/equation[i];  
42.             printf("x%d = %f \n",self,xs);  
43.             MPI_Bcast(&xs,1,MPI_FLOAT,i,MPI_COMM_WORLD);  
44.         } else {  
45.             MPI_Bcast(&xs,1,MPI_FLOAT,i,MPI_COMM_WORLD);  
46.             if(i>self) {  
47.                 equation[xcnt]-=equation[i]*xs;  
48.             }  
49.         }  
50.     }  
51.     free(equation);  
52.     MPI_Finalize();  
53.     return 0;  
54. }