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

• 2019-01-07
• 2,300
• 2
• 1

代码是在 https://blog.csdn.net/a337875954/article/details/79042062 的基础上修改的

原文存在错误，在计算上三角矩阵时，只对系数矩阵进行了运算，而不是对增广矩阵进行运算。

我只是对这个错误进行了修改，并优化、精简了输出。

代码：

#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;
}