Форум программистов [Powered by Invision Power Board]

Код

void main(int argc, char **argv)
{
    int  id, N, i, j, s, k, l, index, m, n, p, locN, t;
    double starttime, endtime, max, c;
    double **a, *b, *buf;

    MPI_Init (&argc, &argv);
    starttime = MPI_Wtime ();

    //n = atoi (argv[1]);
    n = 16;     // matrix size

    MPI_Comm_size (MPI_COMM_WORLD,  &N);
    MPI_Comm_rank (MPI_COMM_WORLD,  &id);

    if (n % N)
    {
        MPI_Finalize ();
        printf ("Matrix size (%i) can't be divided to the number of processors (%i)\n",  n,  N);
        exit (1);
    }

    locN = n/N;
    if (id == 0)
        printf ("locN = %i\n\n", locN);

    a  =  (double**) malloc (locN * sizeof (double));       // create matrix a (for each process)
    if (id == 0)
        b = (double*) malloc (n * sizeof (double));         // create vector b (for bm only)

    buf = (double*) malloc (n * sizeof (double));

    // fill matrix a[locN x n]
    for (i=0; i<locN; i++)
    { 
        a[i] = (double*) malloc (n * sizeof (double));
        for (j=0; j<n; j++)
        {
            if (j == i + id * locN)
                a[i][j] = j+1;
            else if (j == i + id * locN + 1)
                a[i][j] = 2;
            else
                a[i][j] = 0;
        }
    }

    // fill vector b
    if (id == 0)
        for (i = 0;i<n;i++)
            b[i] = (i+1) * 2;

    if (id == 0)
    {
        printf("Vector:\n");

        for (i=0; i<n; i++)
            printf("% 7.1f", b[i]);

        printf ("\n\n");
    }

    if (id == 0)
    {
        printf ("Matrix:\n");
        for (i=0; i<locN; i++)
        {
            for (j=0; j<n; j++)
            {
                printf ("% 7.1f", a[i][j]);
            }
            printf ("\n");
        }
        printf ("\n");
    }
   
    for (p=0; p<N; p++)
    {
        for (i=0; i<locN; i++)
        {
            // find max element in matrix a for the current process
            if (id == p)
            {
                max = a[i][i + p*locN];
                index = i + p*locN;

                for (j=i + p*locN + 1; j<n; j++)
                {
                    if (a[i][j] > max)
                    {
                        max = a[i][j];
                        index = j;
                    }
                }
            }

            // broadcast to other processes
            MPI_Bcast (&index, 1, MPI_INT, p, MPI_COMM_WORLD);
            MPI_Bcast (&max, 1, MPI_DOUBLE, p, MPI_COMM_WORLD);
           
            /*
            printf("id: %d   max: %f  index: %d\n", id, max, index);   
            printf("Matrix before subst ... \n");

            for (k = 0;k<locN;k++)
            {
            for (s = 0;s<n;s++)
            printf("id: %d  a: %.2f ", id, a[k][s]);
            printf("\n");
            }
            */
           
            if (index != i + p*locN)
            {   
                for (l=0; l<locN; l++)
                {
                    c = a[l][i + p*locN];
                    a[l][i + p*locN] = a[l][index];
                    a[l][index] = c;
                }

                c = b[i + p*locN];
                b[i + p*locN] = b[index];
                b[index] = c;
            }


            /*    printf("Matrix after subst ... \n");
            for (k = 0;k<locN;k++)
            {
            for (s = 0;s<n;s++)
            printf("id: %d  a: %.2f ", id, a[k][s]);
            printf("\n");
            }
            */       
            if (id == p)
            {
                for (k = 0;k<n;k++)
                {
                    buf[k] = a[i][k];
                    /* printf("id: %d a[%d][%d]: %.2f\n", id, i, k, a[i][k]);*/
                }
            }

            MPI_Bcast(buf, n, MPI_DOUBLE, p, MPI_COMM_WORLD);
            /*    for (s = 0;s<n;s++)
            printf("id: %d  buf: %f\n", id, buf[s]);

            printf("New matrix before ... \n");
            for (s = 0;s<locN;s++)
            {
            for (k = 0;k<n;k++)
            printf("id: %d  a: %.2f ", id, a[s][k]);
            printf("\n");
            }
            */       
            for (s = 0;s<locN;s++)
                for (t = 0;t<n;t++)
                    if (t!= i+p*locN)
                        a[s][t] -= a[s][i+p*locN]*(buf[t]/buf[i+p*locN]);     

            if (id == 0)
                for (k = 0;k<n;k++)
                    if (k!= i+p*locN)
                        b[k] -= b[i+p*locN]*(buf[k]/buf[i+p*locN]);

            MPI_Barrier(MPI_COMM_WORLD);
            /*    printf("New matrix after next step... \n", i+1);
            for (k = 0;k<locN;k++)
            {
            for (j = 0;j<n;j++)
            printf("id: %d  a: %.2f ", id, a[k][j]);
            printf("\n");
            }
            */
        }         
    }           

    MPI_Barrier(MPI_COMM_WORLD);
    /* if (id == 0) printf("END MATRIX ... \n");
    for (i = 0;i<locN;i++)
    {
    for (j = 0;j<n;j++)
    printf("id: %d  a: %.2f ", id, a[i][j]);
    printf("\n");
    }
    if (id == 0)
    {     
    printf("END VECTOR ... \n");
    for (k = 0;k<n;k++)
    printf("%f\n", b[k]);
    }
    */
    MPI_Barrier(MPI_COMM_WORLD);
    MPI_Bcast(b, n, MPI_DOUBLE, 0, MPI_COMM_WORLD);
    /*for (i = 0;i<n;i++)
    printf("id: %d   b: %f\n", id, b[i]);
    if (id == 0)
    printf("So,  we find the solution of SLAE ... \n");
    for (p = 0;p<locN;p++)
    printf("id: %d res: %f\n", id, b[p+id*locN]/a[p][p+id*locN]);
    */   
    MPI_Barrier(MPI_COMM_WORLD);
    endtime = MPI_Wtime()-starttime;     
    if (id == 0) 
        printf("Time %f\n", endtime);

    for (i = 0;i<locN;i++)
        free(a[i]);
    free(a);
    free(b);
    free(buf);
    MPI_Finalize();
    exit(1);
}