home nsu nsulab NGU-2009 7201 ikuznecov nsulab smp4x64 NGU-2009 7201 i

/*
 * /home/nsu/nsulab/NGU-2009/7201/ikuznecov/1.c
 *
 * nsulab@smp4x64:~/NGU-2009/7201/ikuznecov$ mpicc 1.c -DDEBUG
 * nsulab@smp4x64:~/NGU-2009/7201/ikuznecov$ mpirun -np 4 ./a.out
 * A0: 0.00 1.00 2.00 3.00 
 * A1: 4.00 5.00 6.00 7.00 
 * A3: 12.00 13.00 14.00 15.00 
 * A2: 8.00 9.00 10.00 11.00 
 * 6.00 6.00 6.00 6.00 
 * 22.00 22.00 22.00 22.00 
 * 38.00 38.00 38.00 38.00 
 * 54.00 54.00 54.00 54.00 
*/

#include <stdio.h>
#include <stdlib.h>
#include <string.h> /* memset */
#include <mpi.h>
#include <sys/time.h>

/* матрица NxN */
#define N 4
#define DIMS 1


int make_limits (int size, int *limits)
{
        int i, n, max = 0;

        memset(limits, 0, sizeof(int) * size);

        for (i = 0; i < size; i++)
        {
                limits[i] = N / size; /* всем нодам одинаковое количество строк */
        }
        
        /* раздаём остаток */
        n = N % size;

        while (n > 0)
        {
                for (i = 0; i < size && n > 0; i++)
                {
                        limits[i]++;
                        n--;
                }
        }
        
        /* ищем размер самого большого куска, по нему будем выделять память для B */
        for (i = 1; i < size; i++)
        {
                if (limits[i] > limits[max])
                        max = i;
        }

        return limits[max];
}

int main (int argc, char *argv[])
{
        int i, j, k, l;
        
        int offset = 0;
        int current_lim = 0, max_lim = 0;
        
        int *limits;
        double *A, *B, *C;
        
        int size, rank;
        int dims[DIMS], periods[DIMS];
        int src = 0, dst = 0;
        MPI_Status st;
        MPI_Comm ring;

        MPI_Init(&argc, &argv);

        MPI_Comm_size(MPI_COMM_WORLD, &size);
        MPI_Comm_rank(MPI_COMM_WORLD, &rank);

        for (i = 0; i < DIMS; i++)
        {
                dims[i] = 0;
                periods[i] = 1;
        }

        MPI_Dims_create(size, DIMS, dims);
        MPI_Cart_create(MPI_COMM_WORLD, DIMS, dims, periods, 0, &ring); /* создаём коммуникатор с топологией "кольцо" */
        MPI_Comm_rank(ring, &rank);

        limits = malloc(sizeof(int)*size);

        max_lim = make_limits(size, limits);

        A = malloc(sizeof(double) * N * limits[rank]);
        B = malloc(sizeof(double) * N * max_lim);
        C = malloc(sizeof(double) * N * limits[rank]);

        /* заполняем свои части (limits[rank] x N)
         * матриц A (0.0, 1.0, 2.0, ...), B (1.0, 1.0, 1.0, ...) и C (0.0, 0.0, 0.0, ...)
        */
        for(i = 0; i < limits[rank]; i++)
        {
                for(j = 0; j < N; j++)
                {
                        A[i*N+j] = N*limits[rank]*rank + i*N + j;
                        B[i*N+j] = 1.;
                        C[i*N+j] = 0.;
                }
        }

#ifdef DEBUG
        for (i = 0; i < limits[rank]; i++)
        {
                printf ("A%d: ", rank);
                for (j = 0; j < N; j++)
                        printf("%.2f ", A[i*N+j]);
                printf("\n");
        }
#endif

        for (i = 0, offset = 0; i < rank; i++)
                offset += limits[i]; /* сдвиг в матрице до текущей ноды */

        current_lim = limits[rank]; /* нужен в отдельной памяти, т.к. будет передаваться */
        for (l = 0; l < size; l++) /* по всем нодам, считаем произведение текущих кусков A и B; передаем матрицу B и её размер по кругу */
        {
                for (i = 0; i < limits[rank]; i++)
                        for (j = 0; j < current_lim; j++)
                                for (k = 0; k < N; k++)
                                        C[i*N+j+offset] += A[i*N+k] * B[j*N+k];

                offset = (offset + current_lim) % N;

                MPI_Cart_shift(ring, 0, -1, &src, &dst); /* узнаём ранк предыдущей и следующей по кольцу нод */
                MPI_Sendrecv_replace(B, max_lim*N, MPI_DOUBLE, dst, 1, src, 1, ring, &st);
                MPI_Sendrecv_replace(&current_lim, 1, MPI_INT, dst, 2, src, 2, ring, &st);
        }

        for(k = 0; k < size; k++)
        {
                MPI_Barrier(MPI_COMM_WORLD); /* ждём пока все ноды дойдут до барьера */
                if (k == rank) /* каждая нода по порядку выводит свою часть матрицы C */
                {
                        for (i = 0; i < limits[k]; i++)
                        {
                                for (j = 0; j < N; j++)
                                        printf ("%.2f ", C[i*N+j]);
                                printf ("\n");
                        }
                }
                fflush(stdout);
                MPI_Barrier(MPI_COMM_WORLD);
        }

        free(A);
        free(B);
        free(C);
        free(limits);

        MPI_Comm_free(&ring);
        MPI_Finalize();
        return 0;
}