冬Blog

醉心技术、醉心生活
  博客园  :: 首页  :: 新随笔  :: 订阅 订阅  :: 管理

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

Posted on 2007-10-24 09:12  冬冬  阅读(5536)  评论(13编辑  收藏  举报
  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