/*
* /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(¤t_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;
}