jacobi.c

利用MPI实现并行化的jacobi算法

/****************************************************************************
程序名：jacobi.c
程序功能：
1、利用MPI实现利用jacobi算法实现n*n线性方程组求解的并行计算
2、在主节点进行串行jacobi算法，对同一方程组求解
3、对比执行时间
****************************************************************************/
//需引用头文件：
#include "mpi.h"
#include <stdio.h>
#include <stdlib.h>

//常量定义区
#define MAX_Node 100 //允许的最大未知数个数
#define MAX_BufSize (MAX_Node * MAX_Node) //允许最大的系数的个数

#define MAX_loops 10000 //最大迭代次数
#define TOLERANCE 0.0001 //误差

//变量定义区
int pID, pSize; //pID：当前进程ID，pSize：总的进程数
int size,n, iteration = 0; //n：未知数的个数，iternation：迭代的次数
float x[MAX_Node], new_x[MAX_Node], result_x[MAX_Node]; //x：表示上一次迭代的结果，new_x：表示这一次迭代的结果，result_x：表示归约之后得到的总的结果
float a[MAX_Node][MAX_Node]; //系数
float b[MAX_Node]; 
FILE *finput;//输入数据

//函数定义区
float norm_inf(float x[],int n)
{
 float norm;
 int i;
 norm=fabs(x[0]);
 for(i=1;i<n;i++)
 {
  if(fabs(x[i])>norm)
   norm=fabs(x[i]);
 }
 return norm;
} 
//串行算法实现函数：
void serial_jacobi(float a[MAX_Node][MAX_Node],float g[MAX_Node],int n)
{
 float b[MAX_Node][MAX_Node]={0},x0[MAX_Node],x1[MAX_Node],x1_x0[MAX_Node],norm,temp;
 int i,j,k;
 //初始化结果集：
 for(i=0;i<n;i++)
 {
  g[i]=g[i]/a[i][i];
  for(j=0;j<n;j++)
  {
   if(i==j)
    continue;
   b[i][j]=-a[i][j]/a[i][i];
  }
 }
 
 for(i=0;i<n;i++)
 {
  x0[i]=0;
  x1[i]=1;
  x1_x0[i]=x1[i]-x0[i];
 }
 k=0;
 norm=norm_inf(x1_x0,n);
 
 //叠代计算新的结果集
 while((norm>=TOLERANCE)&&(k<MAX_loops))
 {
  for(i=0;i<n;i++)
   x0[i]=x1[i];
  for(i=0;i<n;i++)
  {
   temp=0;
   for(j=0;j<n;j++)
    temp=temp+b[i][j]*x0[j];
   x1[i]=temp+g[i];
   x1_x0[i]=x1[i]-x0[i];
  }
  norm=norm_inf(x1_x0,n);
  k++;
 }
 //结果输出
 for(i=0;i<n;i++)
  printf("x[%d]=%lf\n",i,x1[i]);
 printf("Serial Total Iteration : %d\n", k);

}
//读取输入文件
void input(){
	
    
    int i, j;
    printf("The n is %d\n", n);
    fflush(stdout);
        
		fscanf(finput,"%d", &size);
    // 判断是否是合法输入,如果不是,程序退出
    if(size != n)
    {
        puts("The input is error!");
        MPI_Finalize();
        exit(0);
    }
    //读入A的系数和b,并保存到aij和bj中
    for(i = 0; i < size; i ++)
    {
        for(j = 0; j < size; j ++) fscanf(finput, "%f", &a[i][j]);
        fscanf(finput, "%f", &b[i]);
    }
    
    //输出待测内容：
    for(i = 0; i < size; i ++)
    {
    		for(j = 0; j < size; j ++) 
    		{
    			if (j==size-1) printf("%f*X[%d]\t",a[i][j],j);
    				else printf("%f*X[%d] + ",a[i][j],j);
    		}
     		printf(" = %f\n",b[i]); 			
  	}
    fclose(finput);
}
//输出结果
void p_output(){
    
    int i;
    printf("parallel jacobi：\n");
        if(iteration > MAX_loops) printf(", that is out of the limit of iteration!\n");
    for(i = 0; i < n; i++)
        printf("x[%d] is %f\n", i, result_x[i]);
    printf("Parallel Total Iteration : %d\n", iteration);
}
//判断精度是否满足要求，满足要求则返回1，否则，返回0
char tolerance(){
    
    int i;
    
    //有一个不满足误差的，则返回0
    for(i = 0; i < n; i++)
        if(result_x[i] - x[i] > TOLERANCE || x[i] - result_x[i] > TOLERANCE)
            return 0;
    
    //全部满足误差，返回1
    return 1;
}

//入口，主函数
int main(int argc, char* argv[]) {
    
    MPI_Status status; 
    int i, j;
    float sum;
    double starttime, endtime; 
    char	CfgFileName[100];
    
    //输入文件：默认为input，如加参数可指定读入文件名，但必须和可执行文件处于同一目录下
    if (argc == 1)
			finput=fopen("input","r");
		else if (argc == 2)
			sprintf(CfgFileName, "./%s", argv[1]);
    finput=fopen(CfgFileName,"r");
     
    //初始化
    MPI_Init(&argc, &argv); 
    MPI_Comm_rank(MPI_COMM_WORLD, &pID); 
    MPI_Comm_size(MPI_COMM_WORLD, &pSize); 
   
   	    
    //每个进程对应一个未知数
    n = pSize; 

    //根进程负责输入
    if(!pID) input();
    
    //开始计时
    starttime = MPI_Wtime(); 
        
    //广播系数
    MPI_Bcast(a, MAX_BufSize, MPI_FLOAT, 0, MPI_COMM_WORLD);
    //广播b
    MPI_Bcast(b, MAX_Node, MPI_FLOAT, 0, MPI_COMM_WORLD);
    
    //初始化x的值
    for(i = 0; i < n; i++) {
        x[i] = b[i];
        new_x[i] = b[i];
    }
    
       //当前进程处理的元素为该进程的ID
    i = pID;
    
    //迭代求x[i]的值
    do{
        //迭代次数加1
        iteration++;
        
        //将上一轮迭代的结果作为本轮迭代的起始值，并设置新的迭代结果为0
        for(j = 0; j < n; j++)
        {
            x[j] = result_x[j];
            new_x[j] = 0;
            result_x[j] = 0;
        }            
        
        //同步等待
        MPI_Barrier(MPI_COMM_WORLD);
          
        //求和
        sum = - a[i][i] * x[i];
        for(j = 0; j < n; j++) sum += a[i][j] * x[j];
        
        //求新的x[i]
        new_x[i] = (b[i] - sum) / a[i][i];
          
        //使用归约的方法，同步所有计算结果
        MPI_Allreduce(new_x, result_x, n, MPI_FLOAT, MPI_SUM, MPI_COMM_WORLD);
        
        //如果迭代次数超过了最大迭代次数则退出
        if(iteration > MAX_loops) {
            break;
        }
    } while(!tolerance()); //精度不满足要求继续迭代
    
    //根进程负责输出结果
    if(!pID) {
    p_output();
        endtime   = MPI_Wtime(); 
   	printf("parallel took %f seconds\n",endtime-starttime); 
  	}
    //并行结束
   
   	//主进程运算串行算法
    if(!pID) {
    printf("serial jacobi：\n");

    starttime = MPI_Wtime();
    serial_jacobi(a,b,n);
    endtime   = MPI_Wtime(); 

   	printf("serial took %f seconds\n",endtime-starttime); 
  	}
    
    MPI_Finalize(); 
    return 0;
}