高斯消去法解线性方程组（MPI）

用一上午的时间，用MPI编写了高斯消去法解线性方程组。这次只是针对单线程负责一个线程方程的求解，对于超大规模的方程组，需要按行分块，后面会在这个基础上进行修改。总结一下这次遇到的问题：

（1）MPI_Allreduce（）函数的使用；

（2）MPI_Allgather（）函数的使用；

（3）线程之间不使用通信函数进行值传递（地址传递）是没有办法使用其他线程的数据，这是设计并行程序中最容易忽视的一点。

#include "stdio.h"
#include "mpi.h"
#include "malloc.h"

#define n 4
#define BLOCK_LOW(id,p,n) ((id) * (n)/(p))
#define BLOCK_HIGH(id,p,n) (BLOCK_LOW((id)+1,p,n)-1)
#define BLOCK_SIZE(id,p,n) (BLOCK_LOW((id)+1,p,n)-BLOCK_LOW((id),p,n))
#define BLOCK_OWNER(index,p,n) (((p) * ((index)+1)-1)/(n))
#define MIN(a,b) ((a)<(b)?(a):(b))

int NotIn(int id,int *picked);

struct {
double value;
int index;
}local,global;

int main(int argc,char *argv[])
{
  int id,p;
  MPI_Init(&argc,&argv);
  MPI_Comm_rank(MPI_COMM_WORLD,&id);
  MPI_Comm_size(MPI_COMM_WORLD,&p);
  
  //double a[n][n+1] = {{4,6,2,-2,8},{2,0,5,-2,4},{-4,-3,-5,4,1},{8,18,-2,3,40}};
  double a[n][n+1] = {{2,1,-5,1,8},{1,-3,0,-6,9},{0,2,-1,2,-5},{1,4,-7,6,0}};
  
  int i,j;
  int index;
  
  int *picked;
  picked = (int *)malloc(sizeof(int) * n);   //记录已被选中的行
  for(i=0;i<n;i++)
    picked[i] = -1;
  
  for(i=0;i<n;i++)
  {
      
      if(NotIn(id,picked))    //判断该行是否已被选中，没有选择则进行下一步
      {
      local.value = a[id][i];
      local.index = id;
      }
      else
      {
      local.value = -100;   //假设各个参数最小值不小于-100
          local.index = id;
      }
      
      MPI_Allreduce(&local,&global,1,MPI_DOUBLE_INT,MPI_MAXLOC,MPI_COMM_WORLD);   // 归约最大的值，并存入到global中
     // printf(" i = %d,id =%d,value = %f,index = %d
",i,id,global.value,global.index);
      picked[i] = global.index;
      
    if(id == global.index)
     {
       MPI_Bcast(&global.index,1,MPI_INT,id,MPI_COMM_WORLD);
     }
     
     MPI_Allgather(&a[id][0],n+1,MPI_DOUBLE,a,n+1,MPI_DOUBLE,MPI_COMM_WORLD);  
         //每个线程解决的是对应行的求解，例如：线程号为0的线程仅得到0行的解，但是第1行的改动，0线程没有办法得到，只有1线程自己才知道，所以需要使用MPI_Allg          ather（）函数进行去收集，并将结果存入到各个线程中，最后各个线程得到a为最新解
    
     if(NotIn(id,picked))
     {
       double temp = 0 - a[id][i] / a[picked[i]][i];
       for(j=i;j<n+1;j++)
       {
           a[id][j] += a[picked[i]][j] * temp;
       }
     }

  }
 
  MPI_Gather(&a[id][0],n+1,MPI_DOUBLE,a,n+1,MPI_DOUBLE,0,MPI_COMM_WORLD);  //
      //解上三角形，因为都需要使用到上一线程的数值，这样造成通信开销增大，不如直接让单一线程去求解上三角形矩阵
  if(id == 0)
  {
     for(i=0;i<n;i++)
     {
    for(j=0;j<n+1;j++)
    {
       printf("%f	",a[i][j]);
     }
     printf("
");
     }
     
    /* for(i=0;i<n;i++)
     {
       printf("%d	",picked[i]);
     }
     */
     double *x;
     x = (double *)malloc(sizeof(double) * n);
     for(i=(n-1);i>=0;i--)   //这里还有一个小插曲，在这一行的后面加了”；“，导致i变成-1，使程序报错 Segmentation fault (11)，在Linux下经常遇到这个问题，大体2点：1.占用的内存超出了系统内存 2.循环中，数组越界

     {
        //printf("
 n =%d,i = %d",n,i);
        x[i] = a[picked[i]][n] / a[picked[i]][i];
        printf("x[%d] = %f
",i,x[i]);
        for(j=0;j<n;j++)
        {
            a[picked[j]][n] = a[picked[j]][n] - x[i] * a[picked[j]][i] ;
            a[picked[j]][i] = 0;
        }
     
     }
   }


  MPI_Finalize();
  return 0;
}


int NotIn(int id,int *picked)
{
  int i;
  for(i=0;i<n;i++)
  {
    if(id == picked[i])
    {
      return 0;
    }
  }
  return 1;
}