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

  1 #include "stdio.h"
  2 #include "mpi.h"
  3 #include "malloc.h"
  4 
  5 #define n 4
  6 #define BLOCK_LOW(id,p,n) ((id) * (n)/(p))
  7 #define BLOCK_HIGH(id,p,n) (BLOCK_LOW((id)+1,p,n)-1)
  8 #define BLOCK_SIZE(id,p,n) (BLOCK_LOW((id)+1,p,n)-BLOCK_LOW((id),p,n))
  9 #define BLOCK_OWNER(index,p,n) (((p) * ((index)+1)-1)/(n))
 10 #define MIN(a,b) ((a)<(b)?(a):(b))

 11 int NotIn(int id,int *picked);

 12 struct {
 13 double value;
 14 int index;
 15 }local,global;
 16 
 17 int main(int argc,char *argv[])
 18 {
 19   int id,p;
 20   MPI_Init(&argc,&argv);
 21   MPI_Comm_rank(MPI_COMM_WORLD,&id);
 22   MPI_Comm_size(MPI_COMM_WORLD,&p);
 23   
 24   //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}};
 25   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}};
 26   
 27   int i,j;
 28   int index;
 29   
 30   int *picked;
 31   picked = (int *)malloc(sizeof(int) * n);   //记录已被选中的行
 32   for(i=0;i<n;i++)
 33     picked[i] = -1;
 34   
 35   for(i=0;i<n;i++)
 36   {
 37       
 38       if(NotIn(id,picked))    //判断该行是否已被选中，没有选择则进行下一步
 39       {
 40       local.value = a[id][i];
 41       local.index = id;
 42       }
 43       else
 44       {
 45       local.value = -100;   //假设各个参数最小值不小于-100
 46           local.index = id;
 47       }
 48       
 49       MPI_Allreduce(&local,&global,1,MPI_DOUBLE_INT,MPI_MAXLOC,MPI_COMM_WORLD);   // 归约最大的值，并存入到global中
 50      // printf(" i = %d,id =%d,value = %f,index = %d\n",i,id,global.value,global.index);
 51       picked[i] = global.index;
 52       
 53     if(id == global.index)
 54      {
 55        MPI_Bcast(&global.index,1,MPI_INT,id,MPI_COMM_WORLD);
 56      }
 57      
 58      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为最新解
 59     
 60      if(NotIn(id,picked))
 61      {
 62        double temp = 0 - a[id][i] / a[picked[i]][i];
 63        for(j=i;j<n+1;j++)
 64        {
 65            a[id][j] += a[picked[i]][j] * temp;
 66        }
 67      }
 68 
 69   }
 70  
 71   MPI_Gather(&a[id][0],n+1,MPI_DOUBLE,a,n+1,MPI_DOUBLE,0,MPI_COMM_WORLD);  //
      //解上三角形，因为都需要使用到上一线程的数值，这样造成通信开销增大，不如直接让单一线程去求解上三角形矩阵
 72   if(id == 0)
 73   {
 74      for(i=0;i<n;i++)
 75      {
 76     for(j=0;j<n+1;j++)
 77     {
 78        printf("%f\t",a[i][j]);
 79      }
 80      printf("\n");
 81      }
 82      
 83     /* for(i=0;i<n;i++)
 84      {
 85        printf("%d\t",picked[i]);
 86      }
 87      */
 88      double *x;
 89      x = (double *)malloc(sizeof(double) * n);
 90      for(i=(n-1);i>=0;i--)   //这里还有一个小插曲，在这一行的后面加了”；“，导致i变成-1，使程序报错 Segmentation fault (11)，在Linux下经常遇到这个问题，大体2点：1.占用的内存超出了系统内存 2.循环中，数组越界

 91      {
 92         //printf("\n n =%d,i = %d",n,i);
 93         x[i] = a[picked[i]][n] / a[picked[i]][i];
 94         printf("x[%d] = %f\n",i,x[i]);
 95         for(j=0;j<n;j++)
 96         {
 97             a[picked[j]][n] = a[picked[j]][n] - x[i] * a[picked[j]][i] ;
 98             a[picked[j]][i] = 0;
 99         }
100      
101      }
102    }
103 
104 
105   MPI_Finalize();
106   return 0;
107 }
108 
109 
110 int NotIn(int id,int *picked)
111 {
112   int i;
113   for(i=0;i<n;i++)
114   {
115     if(id == picked[i])
116     {
117       return 0;
118     }
119   }
120   return 1;
121 }