LU分解连续分块的实现

墨蓝 2022-05-14 06:10 201阅读 0赞

#include "string.h"
    #include "stdio.h"
    #include "stdlib.h"
    #include "mpi.h"
    
    #define a(x,y) a[x*M+y]
    
    #define A(x,y) A[x*M+y]
    #define l(x,y) l[x*M+y]
    #define u(x,y) u[x*M+y]
    #define floatsize sizeof(float)
    #define intsize sizeof(int)
    
    int M,N;
    int m;
    float *A;
    int my_rank;
    int p;
    MPI_Status status;
    
    void fatal(char *message)
    {
        printf("%s\n",message);
        exit(1);
    }
    
    void Environment_Finalize(float *a,float *f)
    {
        free(a);
        free(f);
    }
    
    int main(int argc, char **argv)
    {
        int i,j,k,my_rank,group_size;
        int i1,i2;
        int v,w;
        float *a,*f,*l,*u;
        FILE *fdA;
        double time1,time2;
    
        MPI_Init(&argc,&argv);
        MPI_Comm_size(MPI_COMM_WORLD,&group_size);
        MPI_Comm_rank(MPI_COMM_WORLD,&my_rank);
    
        p=group_size;
        time1=MPI_Wtime();
        if (my_rank==0)
        {
            fdA=fopen("dataIn.txt","r");
            fscanf(fdA,"%d %d", &M, &N);
            if(M != N)
            {
                puts("The input is error!");
                exit(0);
            }
            A=(float *)malloc(floatsize*M*M);
            for(i = 0; i < M; i ++)
                for(j = 0; j < M; j ++)
                    fscanf(fdA, "%f", A+i*M+j);
            fclose(fdA);
            printf("Input of file \"dataIn.txt\"\n");
            printf("%d\t %d\n",M, N);
            for(i=0;i<M;i++)
            {
                for(j=0;j<N;j++)
                    printf("%f ",A(i,j));
                printf("\n");
            }
        }
    
       
        MPI_Bcast(&M,1,MPI_INT,0,MPI_COMM_WORLD);
        m=M/p;
        if (M%p!=0) m++;
    
        
        a=(float*)malloc(floatsize*m*M);
    
     
        f=(float*)malloc(floatsize*M);
    
       
        if (my_rank==0)
        {
            l=(float*)malloc(floatsize*M*M);
            u=(float*)malloc(floatsize*M*M);
        }
    
       
        if (a==NULL) fatal("allocate error\n");
    
        if (my_rank==0)
        {
            for(i=0;i<m;i++)
                for(j=0;j<M;j++)
                    a(i,j)=A(i,j);
            for(i=0;i<M;i++)
                if (i>=M/p)
            {
                  i1 = i/m;
                  i2 = i%m+1;
            //    i1=i%p;
            //    i2=i/p+1;
                MPI_Send(&A(i,0),M,MPI_FLOAT,i1,i2,MPI_COMM_WORLD);
            }
        }
        else
        {
            for(i=0;i<m;i++)
                MPI_Recv(&a(i,0),M,MPI_FLOAT,0,i+1,MPI_COMM_WORLD,&status);
        }
    
        for(i=0;i<p;i++)
            for(j=0;j<m;j++)        
        {
            /*i*/
            if (my_rank==i)
    //        if (my_rank==j)
            {
    //            v=i*p+j;
                  v = m*i +j;
                for (k=v;k<M;k++)
                    f[k]=a(j,k);
    
                MPI_Bcast(f,M,MPI_FLOAT,my_rank,MPI_COMM_WORLD);
            }
            else
            {
    //          v=i*p+j;
                v = m*i +j;
                MPI_Bcast(f,M,MPI_FLOAT,i,MPI_COMM_WORLD);
            }
    
            /**/
    //        if (my_rank<=i)
            if (my_rank==i)
                for(k=j+1;k<m;k++)
            {
                a(k,v)=a(k,v)/f[v];
                for(w=v+1;w<M;w++)
                    a(k,w)=a(k,w)-f[w]*a(k,v);
            }
    
            /**/
            if (my_rank>i)
    //        if (my_rank>j)
                for(k=0;k<m;k++)
            {
                a(k,v)=a(k,v)/f[v];
                for(w=v+1;w<M;w++)
                    a(k,w)=a(k,w)-f[w]*a(k,v);
            }
        }
    
        /**/
        if (my_rank==0)
        {
            for(i=0;i<m;i++)
                for(j=0;j<M;j++)
                    A(i,j)=a(i,j);
        }
        if (my_rank!=0)
        {
            for(i=0;i<m;i++)
                MPI_Send(&a(i,0),M,MPI_FLOAT,0,i,MPI_COMM_WORLD);
        }
        else
        {
            for(i=1;i<p;i++)
                for(j=0;j<m;j++)
            {
                MPI_Recv(&a(j,0),M,MPI_FLOAT,i,j,MPI_COMM_WORLD,&status);
                for(k=0;k<M;k++)
                    A((i*m+j),k)=a(j,k);
            }
        }
    
        if (my_rank==0)
        {
            for(i=0;i<M;i++)
                for(j=0;j<M;j++)
                    u(i,j)=0.0;
            for(i=0;i<M;i++)
                for(j=0;j<M;j++)
                    if (i==j)
                        l(i,j)=1.0;
            else
                l(i,j)=0.0;
            for(i=0;i<M;i++)
                for(j=0;j<M;j++)
                    if (i>j)
                        l(i,j)=A(i,j);
            else
                u(i,j)=A(i,j);
            printf("\nOutput of LU operation\n");
            printf("Matrix L:\n");
            for(i=0;i<M;i++)
            {
                for(j=0;j<M;j++)
                    printf("%f ",l(i,j));
                printf("\n");
            }
            printf("Matrix U:\n");
            for(i=0;i<M;i++)
            {
                for(j=0;j<M;j++)
                    printf("%f ",u(i,j));
                printf("\n");
            }
        }
        time2=MPI_Wtime();
        printf("The whole time is %f\n",time2-time1);
        MPI_Finalize();
        Environment_Finalize(a,f);
        return(0);
    }