/****************************************************************************
程序名：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;
}