# MPI高斯消元求解线性方程组

• 2019-01-07
• 3,416
• 2
• 1

### 代码：

``````#include "mpi.h"
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
int main(int argc,char *argv[])
{
int self,size;
MPI_Status s;
//equivalents
float equs[][11]={{1,2,6,3,5,23,14,6,21,14,2},
{2,8,6,10,21,3,14,21,15,23,4},
{6,6,9,21,14,23,6,54,1,14,6},
{21,1,2,32,2,36,21,4,2,15,2},
{2,23,54,3,32,11,12,1,44,1,41},
{21,12,14,1,23,24,14,21,12,11,21},
{1,24,3,36,6,21,4,5,12,21,21},
{21,41,45,12,32,21,41,12,4,6,4},
{6,12,32,21,14,21,4,5,6,21,6},
{21,12,41,8,18,32,16,14,7,17,27}
};  // 修改后的equs数组储存的是增广矩阵
MPI_Init(&argc,&argv);
MPI_Comm_rank(MPI_COMM_WORLD,&self);
MPI_Comm_size(MPI_COMM_WORLD,&size);
float *equivalent=(float*)malloc((size+1)*sizeof(float));
float *recv=(float*)malloc((size+1)*sizeof(float));
float *x=(float*)malloc((size+1)*sizeof(float));
float B[10]={2,4,6,2,41,21,21,4,6,27};  // B数组是记录Ax=B中的B
float temp=0.0;
int marked=0;
int i, j, k;
int flag=0;
//distribute
if(0==self)
{
for(i=0;i<size+1;i++)
{
equivalent[i]=equs[0][i];
}
for(j=1;j<size;j++)
{
MPI_Send(equs[j],size+1,MPI_FLOAT,j,0,MPI_COMM_WORLD);
}
}
else
{
MPI_Recv(equivalent,size+1,MPI_FLOAT,0,0,MPI_COMM_WORLD,&s);
}

//pivot
for(i=0;i<size-1+1;i++)
{
if(self==i)
{
for(k=0;k<size+1;k++)
{
recv[k]=equivalent[k];
}
marked=1;
}
MPI_Bcast(recv,size+1,MPI_FLOAT,i,MPI_COMM_WORLD);
if(0==marked)
{
float ratio=equivalent[i]/recv[i];
for(j=i;j<size+1;j++)
{
equivalent[j]-=recv[j]*ratio;
}
}
}
//backtracking
for (i=size-1; i>=0; i--)
{
if(self==i)
{
j=size-i-1;
k=size-1;
while(j--)
{
temp+=x[k]*equivalent[k];
k--;
}
x[i]=(equivalent[size]-temp)/equivalent[i];
B[i]=equivalent[size];
}
MPI_Bcast(x,size,MPI_FLOAT,i,MPI_COMM_WORLD);
MPI_Bcast(B,size,MPI_FLOAT,i,MPI_COMM_WORLD);
MPI_Barrier(MPI_COMM_WORLD);
}
MPI_Barrier(MPI_COMM_WORLD);
if(size-1==self)
{
if(equivalent[size-1]==0)
{
//            printf("This matix is unsolvable!\n");
flag = 1;
}
else
{
//            printf("This matix is solvable!\n");
}
}
if (flag)
{
free(recv);
free(equivalent);
free(x);
MPI_Finalize();
return 0;
}
if(0==self)
{
//        printf(" The anwer is :\n");
for (i=0; i<size; i++)
{
printf ("%.3f ", x[i]);
}
putchar('\n');
}
free(recv);
free(equivalent);
free(x);
MPI_Finalize();
return 0;
}
``````

#### 评论

• ###### 梁明晗 2020-06-09 10:13回复

MPI_Send(equs[j],size+1,MPI_FLOAT,j,0,MPI_COMM_WORLD);
MPI_Recv(equivalent,size+1,MPI_FLOAT,0,0,MPI_COMM_WORLD,&s);
这里传输的大小为什么和size有关呢

• ###### zjw1111 2020-06-09 10:44回复

MPI_Send/MPI_Recv的第二个参数是需要进行MPI通信的元素数量，equs[j]这一行中有size+1个元素