1 //=================================================================
  2 // Name : 使用雅可比迭代法求解線性方程組
  3 // Author : 袁冬（2107020046）
  4 // Copyright : 中國海洋大學(xué)
  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 //調(diào)試符號
 13 
 14 #define TRUE 1
 15 #define FALSE 0
 16 #define bool int
 17 
 18 #define MAX_N 100 //允許的最大未知數(shù)個(gè)數(shù)
 19 #define MAX_A (MAX_N * MAX_N) //允許最大的系數(shù)的個(gè)數(shù)
 20 
 21 #define MAX_ITERATION 10000 //最大迭代次數(shù)
 22 #define TOLERANCE 0.001 //誤差
 23 
 24 #include "mpi.h"
 25 #include <stdio.h>
 26 #include <stdlib.h>
 27 
 28 int pID, pSize; //pID：當(dāng)前進(jìn)程ID，pSize：總的進(jìn)程數(shù)
 29 int n, iteration = 0; //n：未知數(shù)的個(gè)數(shù)，iternation：迭代的次數(shù)
 30 float x[MAX_N], new_x[MAX_N], result_x[MAX_N]; //x：表示上一次迭代的結(jié)果，new_x：表示這一次迭代的結(jié)果，result_x：表示歸約之后得到的總的結(jié)果
 31 float a[MAX_N][MAX_N]; //系數(shù)
 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     //輸入系數(shù)
 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 //輸出結(jié)果
 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     //有一個(gè)不滿足誤差的，則返回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 //入口，主函數(shù)
 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     //每個(gè)進(jìn)程對應(yīng)一個(gè)未知數(shù)
100     n = pSize; 
101 
102     //根進(jìn)程負(fù)責(zé)輸入
103     if(!pID) input();
104         
105     //廣播系數(shù)
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     //當(dāng)前進(jìn)程處理的元素為該進(jìn)程的ID
125     i = pID;
126     
127     //迭代求x[i]的值
128     do{
129         //迭代次數(shù)加1
130         iteration++;
131         
132         //將上一輪迭代的結(jié)果作為本輪迭代的起始值，并設(shè)置新的迭代結(jié)果為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         //打印本輪迭代的結(jié)果
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         //使用歸約的方法，同步所有計(jì)算結(jié)果
167         MPI_Allreduce(new_x, result_x, n, MPI_FLOAT, MPI_SUM, MPI_COMM_WORLD);
168         
169 #ifdef DEBUG
170         //打印歸約后的結(jié)果
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         //如果迭代次數(shù)超過了最大迭代次數(shù)則退出
178         if(iteration > MAX_ITERATION) {
179             break;
180         }
181     } while(!tolerance()); //精度不滿足要求繼續(xù)迭代
182     
183     //根進(jìn)程負(fù)責(zé)輸出結(jié)果
184     if(!pID) output();
185     
186     //結(jié)束
187     MPI_Finalize();
188     return 0;
189 }
190