lu.c

这是用并行程序语言MPI做的关于矩阵LU分解的程序

#include "stdio.h"
#include "stdlib.h"
#include "mpi.h"
#define a(x,y) a[(x)*M+y]
/*A为M*M矩阵*/
#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,z;
    int i1,i2;
    int v,w;
	float h;
    float *a,*f,*l,*u;
    FILE *fdA;
	double starttime;
	double time1,time2;

    MPI_Init(&argc,&argv);
    MPI_Comm_size(MPI_COMM_WORLD,&group_size);
    MPI_Comm_rank(MPI_COMM_WORLD,&my_rank);
	starttime=MPI_Wtime();
    p=group_size;

    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("the original matrix is:\n");
		for(i=0;i<M;i++){
			for(j=0;j<M;j++)
				printf("%f ",A(i,j) );
			printf("\n");
		}
		printf("\n");
*/
    }

    /*0号进程将M广播给所有进程*/
    MPI_Bcast(&M,1,MPI_INT,0,MPI_COMM_WORLD);
    m=M/p;
    if (M%p!=0) m++;

    /*分配至各进程的子矩阵大小为m*M*/
    a=(float*)malloc(floatsize*m*M);

    /*各进程为主行元素建立发送和接收缓冲区*/
    f=(float*)malloc(floatsize*M);

    /*0号进程为l和u矩阵分配内存，以分离出经过变换后的A矩阵中的l和u矩阵*/
    if (my_rank==0)
    {
        l=(float*)malloc(floatsize*M*M);
        u=(float*)malloc(floatsize*M*M);
    }

    /*0号进程采用行交叉划分将矩阵A划分为大小m*M的p块子矩阵，依次发送给1至p-1号进程*/
    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);
		if (M%p==0) 
		{
			for(i=0;i<M;i++)
				if ((i/m)!=0)
				{
					i1=i/m;
					i2=i%m+1;
					MPI_Send(&A(i,0),M,MPI_FLOAT,i1,i2,MPI_COMM_WORLD);
				}
		}
		else
        {
			for(i=0;i<(M%p)*m;i++)
			{
				if ((i/m)!=0)
				{
					i1=i/m;
					i2=i%m+1;
					MPI_Send(&A(i,0),M,MPI_FLOAT,i1,i2,MPI_COMM_WORLD);
				}
			}
			for(i=(M%p)*m;i<M;i++)
			{
				i1 = (i-(M%p)*m)/(m-1)+M%p;
				i2 = (i-(M%p)*m)%(m-1)+1;
				MPI_Send(&A(i,0),M,MPI_FLOAT,i1,i2,MPI_COMM_WORLD);
			}
        }
/*		printf("the rank %d :\n",my_rank);
		for(i=0;i<m;i++){
			for(j=0;j<M;j++)
				printf("%f ",a(i,j));
			printf("\n");
		}
*/	
    }
    else
    {
		if (my_rank<(M%p)||M%p==0) {
			for(i=0;i<m;i++)
				MPI_Recv(&a(i,0),M,MPI_FLOAT,0,i+1,MPI_COMM_WORLD,&status);
		}
		else
		{
			for(i=0;i<m-1;i++)
				MPI_Recv(&a(i,0),M,MPI_FLOAT,0,i+1,MPI_COMM_WORLD,&status);
		}
/*		printf("the rank %d :\n",my_rank);
		for(i=0;i<m;i++){
			for(j=0;j<M;j++)
				printf("%f ",a(i,j));
			printf("\n");
		}	
*/
    }

	time1=MPI_Wtime();

	for(i=0;i<p;i++)
	{
		for(j=0;j<m;j++)
		{
			if(my_rank == i)
			{
				for(k=0;k<M;k++)
				{
					f[k]=a(j,k);
				}
				MPI_Bcast(f,M,MPI_FLOAT,my_rank,MPI_COMM_WORLD );
				for(k=j;k<m-1;k++)
				{
					h = a(k+1,i*m+j)/f[i*m+j];
					for(z=i*m+j;z<M;z++){
						a(k+1,z)+=f[z]*(-1)*h;
					}
					a(k+1,i*m+j) = h;
				}
			}
			else
			{
				MPI_Bcast(f,M,MPI_FLOAT,i,MPI_COMM_WORLD);
			}
			if(my_rank>i)
			{
				for(k=0;k<m;k++)
				{
					h=a(k,i*m+j)/f[i*m+j];					
					for(z=i*m+j;z<M;z++)
					{
						a(k,z)+=-1*f[z]*h;
					}
					a(k,i*m+j) = h;
				}
			}
		}
	}


    /*0号进程从其余各进程中接收子矩阵a，得到经过变换的矩阵A*/
    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+1,MPI_COMM_WORLD);
    }
    else
    {
        for(i=1;i<p;i++)
            for(j=0;j<m;j++)
        {
            MPI_Recv(f,M,MPI_FLOAT,i,j+1,MPI_COMM_WORLD,&status);
            for(k=0;k<M;k++)
			{
				z=i*m+j;
                A(z,k)=f[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("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\t",A(i,j));
            printf("\n");
        }
        printf("\nOutput of LU operation\n");
        printf("Matrix L:\n");
        for(i=0;i<M;i++)
        {
            for(j=0;j<M;j++)
                printf("%f\t",l(i,j));
            printf("\n");
        }
        printf("Matrix U:\n");
        for(i=0;i<M;i++)
        {
            for(j=0;j<M;j++)
                printf("%f\t",u(i,j));
            printf("\n");
        }
    }
	time2=MPI_Wtime();
    if (my_rank==0)
    {
        printf("\n");
        printf("Whole running time    = %f seconds\n",time2-starttime);
        printf("Distribute data time  = %f seconds\n",time1-starttime);
        printf("Parallel compute time = %f seconds\n",time2-time1);

    }
    MPI_Finalize();
    Environment_Finalize(a,f);
    return(0);
}