cjacobi.c

高教出版社的《并行算法实践》一书的源代码

#include "stdio.h"
#include "stdlib.h"
#include "mpi.h"
#include "math.h"
#include "string.h"
#define E 0.00001
#define min -1
#define intsize sizeof(int)
#define charsize sizeof(char)
#define floatsize sizeof(float)
#define A(x,y) A[x*N+y]
#define I(x,y) I[x*N+y]
#define a(x,y) a[x*N+y]
#define e(x,y) e[x*N+y]
#define b(x) b[x]
#define buffer(x,y) buffer[x*N+y]
#define buffee(x,y) buffee[x*N+y]
int M,N;
int *b;
int m,p;
int myid,group_size;
float *A,*I;
float *a,*e;
float max;
float sum;
double starttime,time1,time2,time3,endtime;
MPI_Status status;

float sgn(float v)
{
    float f;
    if (v>=0) f=1;
    else f=-1;
    return f;
}


void read_fileA()
{
    int i,j;

    FILE *fdA;
    time1=MPI_Wtime();

    fdA=fopen("dataIn.txt","r");
    fscanf(fdA,"%d %d", &M, &N);

    if(M != N)
    {
        puts("The input is error!");
        exit(0);
    }

    m=N/p; if (N%p!=0) m++;
    A=(float*)malloc(floatsize*N*m*p);
    I=(float*)malloc(floatsize*N*m*p);

    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\t",A(i,j));
        printf("\n");
    }

    for(i=0;i<N;i++)
    {
        for(j=0;j<N;j++)
        {
            if (i==j) I(i,j)=1;
            else I(i,j)=0;
        }
    }
}


void send_a()
{
    int i;
    for(i=1;i<p;i++)
    {
        MPI_Send(&(A(m*i,0)),m*N,MPI_FLOAT,i,i,MPI_COMM_WORLD);
        MPI_Send(&(I(m*i,0)),m*N,MPI_FLOAT,i,i,MPI_COMM_WORLD);
    }
    free(A);
    free(I);
}


void get_a()
{
    int i,j;
    if (myid==0)
    {
        for(i=0;i<m;i++)
            for(j=0;j<N;j++)
        {
            a(i,j)=A(i,j);
            e(i,j)=I(i,j);
        }
    }
    else
    {
        MPI_Recv(a,m*N,MPI_FLOAT,0,myid,MPI_COMM_WORLD,&status);
        MPI_Recv(e,m*N,MPI_FLOAT,0,myid,MPI_COMM_WORLD,&status);
    }
}


int main(int argc,char **argv)
{
    float *c;
    int k;
    int loop;
    int i,j,v,z,r,t,y;
    int i1,j1;
    float f,g,h;
    float sin1,sin2,cos1;
    float s1,c1;
    float *br,*bt,*bi,*bj,*zi,*zj;
    float *temp1,*temp2,*buffer,*buffee;
    int counter,current;
    int *buf;
    int mml,mpl;
    float bpp,bqq,bpq,bqp;
    float lmax;

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

    max=10.0;
    loop=0;

    if (myid==0) read_fileA();

    MPI_Bcast(&N,1,MPI_INT,0,MPI_COMM_WORLD);

    m=N/p;
    if (N%p!=0) m++;

    a=(float*)malloc(floatsize*m*N);
    e=(float*)malloc(floatsize*m*N);
    b=(int*)malloc(intsize*m);
    br=(float*)malloc(floatsize*N);
    bt=(float*)malloc(floatsize*N);
    bi=(float*)malloc(floatsize*m);
    bj=(float*)malloc(floatsize*m);
    zi=(float*)malloc(floatsize*m);
    zj=(float*)malloc(floatsize*m);
    temp1=(float*)malloc(floatsize*4);
    temp2=(float*)malloc(floatsize*4*p);

    if ((myid==p-1)&&(myid%2!=0))
    {
        buffer=(float*)malloc(floatsize*m/2*N);
        buffee=(float*)malloc(floatsize*m/2*N);
        buf=(int*)malloc(intsize*m/2);            /* store i,j */
    }

    if ((myid%2!=0)&&(myid!=p-1))
    {
        buffer=(float*)malloc(floatsize*m*N);
        buffee=(float*)malloc(floatsize*m*N);
        buf=(int*)malloc(intsize*m);
    }

    if (myid==0)
    {
        get_a();
        send_a();
    }
    else
    {
        get_a();
    }

    if (myid==0)                                  /* start computing now */
        time2=MPI_Wtime();

    for(i=0;i<m;i++)
        b(i)=myid*m+i;

    while (fabs(max)>E)
    {
        loop++;
        for(i=myid*m;i<(myid+1)*m-1;i++)
            for(j=i+1;j<=(myid+1)*m-1;j++)
        {
            r=i%m; t=j%m;                         /* to make a(r,j) zero use a(r,i) a(t,j) */
            if(a(r,j)!=0)
            {
                f=(-a(r,j));
                g=(a(t,j)-a(r,i))/2;
                h=sgn(g)*f/sqrt(f*f+g*g);
                sin2=h;
                sin1=h/sqrt(2*(1+sqrt(1-h*h)));
                cos1=sqrt(1-sin1*sin1);

                bpp=a(r,i)*cos1*cos1+a(t,j)*sin1*sin1+a(r,j)*sin2;
                bqq=a(r,i)*sin1*sin1+a(t,j)*cos1*cos1-a(r,j)*sin2;
                bpq=0; bqp=0;

                for(v=0;v<N;v++)                  /* compute row */
                    if ((v!=i)&&(v!=j))
                {
                    br[v]=a(r,v)*cos1+a(t,v)*sin1;
                    a(t,v)=(-a(r,v))*sin1+a(t,v)*cos1;
                }

                for(v=0;v<N;v++)                  /* row */
                    if ((v!=i)&&(v!=j))
                {
                    a(r,v)=br[v];
                }

                for(v=0;v<m;v++)
                    br[v]=e(v,i)*cos1+e(v,j)*sin1;

                for(v=0;v<m;v++)
                    e(v,j)=e(v,i)*(-sin1)+e(v,j)*cos1;

                for(v=0;v<m;v++)
                    e(v,i)=br[v];

                for(v=0;v<m;v++)                  /* compute col */
                    if ((v!=r)&&(v!=t))
                {
                    bi[v]=a(v,i)*cos1+a(v,j)*sin1;
                    a(v,j)=(-a(v,i))*sin1+a(v,j)*cos1;
                }

                for(v=0;v<m;v++)                  /* col */
                    if ((v!=r)&&(v!=t))
                        a(v,i)=bi[v];

                a(r,i)=bpp;
                a(t,j)=bqq;
                a(r,j)=bpq;
                a(t,i)=bqp;

                temp1[0]=sin1;
                temp1[1]=cos1;
                temp1[2]=(float)i;
                temp1[3]=(float)j;
            }
            else
            {
                temp1[0]=0.0;
                temp1[1]=0.0;
                temp1[2]=0.0;
                temp1[3]=0.0;
            }                                     /* end of if */

            MPI_Allgather(temp1,4,MPI_FLOAT,temp2,4,MPI_FLOAT,MPI_COMM_WORLD);
            current=0;

            for(v=1;v<=p;v++)
            {
                s1=temp2[(v-1)*4+0];
                c1=temp2[(v-1)*4+1];
                i1=(int)temp2[(v-1)*4+2];
                j1=(int)temp2[(v-1)*4+3];

                if ((s1!=0.0)||(c1!=0.0)||(i1!=0)||(j1!=0))
                {
                    if (myid!=current)
                    {

                        for(z=0;z<m;z++)
                        {
                            zi[z]=a(z,i1)*c1 + a(z,j1)*s1;
                            a(z,j1)=-a(z,i1)*s1 + a(z,j1)*c1;
                        }

                        for(z=0;z<m;z++)
                            a(z,i1)=zi[z];

                        for(z=0;z<m;z++)
                            zi[z]=e(z,i1)*c1+e(z,j1)*s1;

                        for(z=0;z<m;z++)
                            e(z,j1)=-e(z,i1)*s1+e(z,j1)*c1;

                        for(z=0;z<m;z++)
                            e(z,i1)=zi[z];

                    }                             /* if myid!=current */

                }                                 /* if */

                current=current+1;
            }                                     /* for v */
        }                                         /* for i,j */

        for(counter=1;counter<=2*p-2;counter++)
        {

            if (myid==0)
            {
                MPI_Send(&(a(m/2,0)),m/2*N,MPI_FLOAT,myid+1,myid+1,MPI_COMM_WORLD);
                MPI_Send(&(e(m/2,0)),m/2*N,MPI_FLOAT,myid+1,myid+1,MPI_COMM_WORLD);
                MPI_Send(&b(m/2),m/2,MPI_INT,myid+1,myid+1,MPI_COMM_WORLD);

                MPI_Recv(&(a(m/2,0)),m/2*N,MPI_FLOAT,myid+1,myid,MPI_COMM_WORLD,&status);
                MPI_Recv(&(e(m/2,0)),m/2*N,MPI_FLOAT,myid+1,myid,MPI_COMM_WORLD,&status);
                MPI_Recv(&b(m/2),m/2,MPI_INT,myid+1,myid,MPI_COMM_WORLD,&status);
            }

            if ((myid==p-1)&&(myid%2!=0))
            {
                for(i=m/2;i<m;i++)
                    for(j=0;j<N;j++)
                        buffer((i-m/2),j)=a(i,j);

                for(i=m/2;i<m;i++)
                    for(j=0;j<N;j++)
                        buffee((i-m/2),j)=e(i,j);

                for(i=m/2;i<m;i++)
                    buf[i-m/2]=b(i);

                for(i=0;i<m/2;i++)
                    for(j=0;j<N;j++)
                        a((i+m/2),j)=a(i,j);

                for(i=0;i<m/2;i++)
                    for(j=0;j<N;j++)
                        e((i+m/2),j)=e(i,j);

                for(i=0;i<m/2;i++)
                    b(m/2+i)=b(i);

                MPI_Recv(&(a(0,0)),m/2*N,MPI_FLOAT,myid-1,myid,MPI_COMM_WORLD,&status);
                MPI_Recv(&(e(0,0)),m/2*N,MPI_FLOAT,myid-1,myid,MPI_COMM_WORLD,&status);
                MPI_Recv(&b(0),m/2,MPI_INT,myid-1,myid,MPI_COMM_WORLD,&status);
                MPI_Send(buffer,m/2*N,MPI_FLOAT,myid-1,myid-1,MPI_COMM_WORLD);
                MPI_Send(buffee,m/2*N,MPI_FLOAT,myid-1,myid-1,MPI_COMM_WORLD);
                MPI_Send(buf,m/2,MPI_INT,myid-1,myid-1,MPI_COMM_WORLD);
            }

            if ((myid==p-1)&&(myid%2==0))
            {
                MPI_Send(&(a(m/2,0)),m/2*N,MPI_FLOAT,myid-1,myid-1,MPI_COMM_WORLD);
                MPI_Send(&(e(m/2,0)),m/2*N,MPI_FLOAT,myid-1,myid-1,MPI_COMM_WORLD);
                MPI_Send(&b(m/2),m/2,MPI_INT,myid-1,myid-1,MPI_COMM_WORLD);

                for(i=0;i<m/2;i++)
                    for(j=0;j<N;j++)
                        a((i+m/2),j)=a(i,j);

                for(i=0;i<m/2;i++)
                    for(j=0;j<N;j++)