用雅可比迭代法求线性方程组的解的并行算法（MPI）

 1  // =================================================================
  2  //  Name : 使用雅可比迭代法求解线性方程组
  3  //  Author : 袁冬（2107020046）
  4  //  Copyright : 中国海洋大学
  5  //  LastUpdate : 2007.10.18    
  6  //  Develop Evrionment : Windows Server 2003 SP2
  7  //                         + Eclipse 3.3.1 + CDT 4.0.1 
  8  //                         + MinGW 3.81 + GDB 6.6 
  9  //                         + mpich2-1.0.6-win32-ia32
 10  // =================================================================
 11 
 12  // #define DEBUG 1  // 调试符号
 13 
 14  #define TRUE 1
 15  #define FALSE 0
 16  #define bool int
 17 
 18  #define MAX_N 100  // 允许的最大未知数个数
 19  #define MAX_A (MAX_N * MAX_N)  // 允许最大的系数的个数
 20 
 21  #define MAX_ITERATION 10000  // 最大迭代次数
 22  #define TOLERANCE 0.001  // 误差
 23 
 24 #include "mpi.h"
 25 #include <stdio.h>
 26 #include <stdlib.h>
 27 
 28  int pID, pSize;  // pID：当前进程ID，pSize：总的进程数
 29  int n, iteration = 0;  // n：未知数的个数，iternation：迭代的次数
 30  float x[MAX_N], new_x[MAX_N], result_x[MAX_N];  // x：表示上一次迭代的结果，new_x：表示这一次迭代的结果，result_x：表示归约之后得到的总的结果
 31  float a[MAX_N][MAX_N];  // 系数
 32  float b[MAX_N]; 
 33 
 34  // 输入
 35  void input(){
 36     
 37      int i, j;
 38     
 39     printf("The n is %d\n", n);
 40     fflush(stdout);
 41     
 42      // 输入系数
 43      for(i = 0; i < n; i++) {
 44         printf("Input a[%d][0] to a[%d][n-1]:\n", i, i);
 45         fflush(stdout);
 46          for(j = 0; j < n; j++)
 47             scanf("%f", &a[i][j]);
 48     }
 49     
 50      // 输入b
 51     printf("Input b[0] to b[n-1]:\n");
 52     fflush(stdout);
 53      for(j = 0; j < n; j++)
 54         scanf("%f", &b[j]);
 55 }
 56  // 输出结果
 57  void output(){
 58     
 59      int i;
 60     
 61     printf("Total iteration : %d", iteration);
 62      if(iteration > MAX_ITERATION) printf(", that is out of the limit of iteration!");
 63     printf("\n");
 64     
 65      for(i = 0; i < n; i++)
 66         printf("x[%d] is %f\n", i, result_x[i]);
 67 }
 68  // 判断精度是否满足要求，满足要求则返回TRUE，否则，返回FALSE
 69  bool tolerance(){
 70     
 71      int i;
 72     
 73      // 有一个不满足误差的，则返回FALSE
 74      for(i = 0; i < n; i++)
 75          if(result_x[i] - x[i] > TOLERANCE || x[i] - result_x[i] > TOLERANCE)
 76              return FALSE;
 77     
 78 #ifdef DEBUG
 79     printf("TRUE From %d\n", pID);
 80     fflush(stdout);
 81  #endif
 82     
 83      // 全部满足误差，返回TRUE
 84      return TRUE;
 85 }
 86 
 87  // 入口，主函数
 88  int main( int argc,  char* argv[]) {
 89     
 90     MPI_Status status; 
 91      int i, j;
 92      float sum;
 93     
 94      // 初始化
 95     MPI_Init(&argc, &argv); 
 96     MPI_Comm_rank(MPI_COMM_WORLD, &pID); 
 97     MPI_Comm_size(MPI_COMM_WORLD, &pSize); 
 98     
 99      // 每个进程对应一个未知数
100     n = pSize; 
101 
102      // 根进程负责输入
103      if(!pID) input();
104         
105      // 广播系数
106     MPI_Bcast(a, MAX_A, MPI_FLOAT, 0, MPI_COMM_WORLD);
107      // 广播b
108     MPI_Bcast(b, MAX_N, MPI_FLOAT, 0, MPI_COMM_WORLD);
109     
110 #ifdef DEBUG
111      // 打印a, b
112      for(j = 0; j < n; j++)
113         printf("%f ", b[j]);
114     printf(" From %d\n", pID);
115     fflush(stdout);
116  #endif
117     
118      // 初始化x的值
119      for(i = 0; i < n; i++) {
120         x[i] = b[i];
121         new_x[i] = b[i];
122     }
123     
124      // 当前进程处理的元素为该进程的ID
125     i = pID;
126     
127      // 迭代求x[i]的值
128      do{
129          // 迭代次数加1
130         iteration++;
131         
132          // 将上一轮迭代的结果作为本轮迭代的起始值，并设置新的迭代结果为0
133          for(j = 0; j < n; j++)
134         {
135             x[j] = result_x[j];
136             new_x[j] = 0;
137             result_x[j] = 0;
138         }            
139         
140          // 同步等待
141         MPI_Barrier(MPI_COMM_WORLD);
142         
143 #ifdef DEBUG
144          // 打印本轮迭代的初始值
145          for(j = 0; j < n; j++)
146             printf("%f ", x[j]);
147         printf(" From %d\n", pID);
148         fflush(stdout);
149  #endif
150         
151          // 求和
152         sum = - a[i][i] * x[i];
153          for(j = 0; j < n; j++) sum += a[i][j] * x[j];
154         
155          // 求新的x[i]
156         new_x[i] = (b[i] - sum) / a[i][i];
157         
158 #ifdef DEBUG
159          // 打印本轮迭代的结果
160          for(j = 0; j < n; j++)
161             printf("%f ", new_x[j]);
162         printf(" From %d\n", pID);
163         fflush(stdout);
164  #endif        
165         
166          // 使用归约的方法，同步所有计算结果
167         MPI_Allreduce(new_x, result_x, n, MPI_FLOAT, MPI_SUM, MPI_COMM_WORLD);
168         
169 #ifdef DEBUG
170          // 打印归约后的结果
171          for(j = 0; j < n; j++)
172             printf("%f ", result_x[j]);
173         printf(" From %d\n", pID);
174         fflush(stdout);
175  #endif        
176         
177          // 如果迭代次数超过了最大迭代次数则退出
178          if(iteration > MAX_ITERATION) {
179              break;
180         }
181     }  while(!tolerance());  // 精度不满足要求继续迭代
182      
183       // 根进程负责输出结果
184      if(!pID) output();
185     
186      // 结束
187     MPI_Finalize();
188      return 0;
189 }
190