Dobra, po kolei.
Masz tutaj poniżej wyczyszczony, dobrze skomentowany kod.
Kopiuj
/*
* Poniższy program oblicza iloczyn:
*
* A x B = C
*
* Korzystając z biblioteki cblas oraz MPI. Idea zrównoleglania jest
* przedstawiona poniżej. Każdy z 2 workerów otrzymuje połowę macierzy A, całą
* macierz B i oblicza połowę macierzy C. Staramy się zachować zasadę prawie
* równego podziału pracy między workery.
*
* worker
* -- -- -- -- -- --
* w1 | a11 a12 a13 a14 | | b11 b12 b13 b14 | | c11 c12 c13 c14 |
* _____| a21 a22 a23 a24 | \/ | b21 b22 b23 b24 | -- | c21 c22 c23 c24 |____
* w2 | a31 a32 a33 a34 | /\ | b31 b32 b33 b34 | -- | c31 c32 c33 c34 |
* | a41 a42 a43 a44 | | b41 b42 b43 b44 | | c41 c42 c43 c44 |
* -- -- -- -- -- --
*/
#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <cblas.h>
#include <mpi.h>
#define IDX(i, j, n) ((i) * (n) + (j))
#define CACHE_LINE_SIZE (64)
static double *
matrix_alloc(int m, int n)
{
double *tmp;
int ret;
/*
* Kilka słów wyjaśnienia, czemu używamy posix_memalign zamiast malloc
* albo calloc. Funkcja posix_memalign gwarantuje utrzymanie wyrównania
* do zadanej potęgi dwójki. My chcemy, aby wyrównanie było do linijki
* cache. Wynika to z dwóch powodów:
*
* 1. Instrukcje ładowania danych do rejestrów dostępne w SSE/AVX/AVX2
* są bardzo czułe na wyrównanie danych. W szczególności wyrównanie
* do rozmiaru linijki cache gwarantuje od razu optymalne wyrównanie
* dla instrukcji ładowania.
* 2. Drugą kwestią jest gwarancja braku tzw. false sharing (to ma sens
* przede wszystkim dla programowania wielowątkowego). Ale dla
* porządku możemy też utrzymać w naszym wieloprocesowym kodzie ten
* sam mechanizm, bo raz, że niewiele kosztuje, dwa - będzie bardziej
* przenośny między środowiskami wieloprocesowymi i wielowątkowymi.
*/
ret = posix_memalign((void **)&tmp, CACHE_LINE_SIZE, m * n *
sizeof(*tmp));
assert(ret == 0);
memset(tmp, 0, n * m * sizeof(*tmp));
return (tmp);
}
static void
matrix_fill_rand(double *M, int m, int n)
{
int i, j;
/* Wypełnij macierz wartościami pseudolosowymi z przedziału [0, 1[. */
for (i = 0; i < m; i++) {
for (j = 0; j < n; j++) {
M[IDX(i, j, m)] = 1.0f * rand() / RAND_MAX;
}
}
}
#ifdef DEBUG
/*
* Procedura testująca, czy wyliczony wynik się zgadza. Wykonywana przez
* workera o id = 0, przeliczająca iloczyn A x B = D i porównująca macierze
* C i D.
*/
static void
check_correctness(double *A, double *B, double *C, int n)
{
double *D;
int i, j, ret;
ret = posix_memalign((void **)&D, 64, n * n * sizeof(*D));
assert(ret == 0);
cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
n, n, n,
1.0,
A, n,
B, n,
0.0,
D, n);
for (i = 0; i < n; i++) {
for (j = 0; j < n; j++) {
assert(C[IDX(i, j, n)] == D[IDX(i, j, n)]);
}
}
free(D);
}
#endif
int
main(int argc, char **argv)
{
double *A, *B, *C;
int i, n, ret;
int comm_size, comm_id;
int *offsets, *sizes;
double t1, t2;
/*
* Dla ćwiczenia zrobimy sobie nowy typ danych w MPI - jeden wiersz
* macierzy. W ramach ćwiczeń - z kodu do mnożenia macierzy kwadratowych
* zrób mnożenie macierzy o dowolnych wymiarach. De facto przeróbek nie
* powinno być wiele.
*/
MPI_Datatype row;
MPI_Status status;
MPI_Init(&argc, &argv);
t1 = t2 = 0;
n = 8;
if (argc > 1)
n = atoi(argv[1]);
MPI_Comm_size(MPI_COMM_WORLD, &comm_size);
MPI_Comm_rank(MPI_COMM_WORLD, &comm_id);
printf("Run as %d among %d processes\n", comm_id, comm_size);
/*
* Utwórzmy specjalny typ danych dla wiersza macierzy. Równie dobrze
* można by tu uzyć typu continuous. Poeksperymentuj.
*/
ret = MPI_Type_vector(1, n, n, MPI_DOUBLE, &row);
assert(ret == MPI_SUCCESS);
ret = MPI_Type_commit(&row);
assert(ret == MPI_SUCCESS);
/*
* Here magic happens. Każdy worker musi wiedzieć, jak duży fragment
* macierzy ma mnożyć. Dzielimy macierz A na comm_size workerów tak,
* aby różnica liczby wierszy przydzielonych do dowolnych 2 workerów
* była nie większa niż 1. Da się to osiągnąć wzorkiem poniżej.
*
* Obie tablice przechowują liczby wierszy:
* sizes[i] - rozmiar slice'a przydzielony do i-tego procesu
* offsets[i] - początek slice'a dla i-tego procesu
*
* Ta metoda ma też taką miłą właściwość, że jeden kod obsługuje zarówno
* przypadek sekwencyjny jak i wieloprocesowy.
*/
sizes = calloc(comm_size, sizeof(*sizes));
assert(sizes != NULL);
offsets = calloc(comm_size, sizeof(*offsets));
assert(offsets != NULL);
for (i = 0; i < comm_size; i++) {
sizes[i] = n / comm_size + (i < (n % comm_size) ? 1 : 0);
if (i > 0)
offsets[i] = offsets[i - 1] + sizes[i - 1];
}
printf("Comm %d has size=%d offset=%d\n", comm_id, sizes[comm_id],
offsets[comm_id]);
B = matrix_alloc(n, n);
if (comm_id == 0) {
A = matrix_alloc(n, n);
C = matrix_alloc(n, n);
matrix_fill_rand(A, n, n);
matrix_fill_rand(B, n, n);
} else {
A = matrix_alloc(sizes[comm_id], n);
C = matrix_alloc(sizes[comm_id], n);
}
t1 = MPI_Wtime();
/* Macierz B jest współdzielona przez wszystkich. */
ret = MPI_Bcast(B, n, row, 0, MPI_COMM_WORLD);
assert(ret == MPI_SUCCESS);
/*
* Rozsyłamy chunki macierzy A z procesu głównego do wszystkich
* pozostałych. Korzystamy z wcześniej przygotowanych tablic sizes
* i offsets.
*/
if (comm_id == 0) {
for (i = 1; i < comm_size; i++) {
ret = MPI_Send(A + IDX(offsets[i], 0, n), sizes[i], row,
i, 0, MPI_COMM_WORLD);
assert(ret == MPI_SUCCESS);
}
} else {
ret = MPI_Recv(A, sizes[comm_id], row, 0, 0, MPI_COMM_WORLD,
&status);
assert(ret == MPI_SUCCESS);
}
/*
* Dany proces oblicza fragment macierzy:
* C[offset:offset+size][0:n] = A[offset:offset+size][0:n] x B[0:n][0:n]
*/
printf("Starting computation in %d\n [%d x %d] * [%d x %d]\n", comm_id,
sizes[comm_id], n, n, n);
/*
* Elegancka metoda mnożenia. Użycie tablic offsetów i rozmiarów
* powoduje, że możemy elegancko zapisać mnożenie. Wszystkie procesy
* robią dokładnie tą samą operację DGEMM z tymi samymi lokalnymi
* parametrami.
*/
cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
sizes[comm_id], /* A ma rozmiar sizes[comm_id] x n */
n, n, /* B ma na sztywno rozmiar n x n */
1.0,
A, n,
B, n,
0.0,
C, n);
/*
* Odbieramy fragmenty macierzy C od workerów. Znowu przydają nam się
* tablice rozmiarów i offsetów.
*/
if (comm_id == 0) {
for (i = 1; i < comm_size; i++) {
ret = MPI_Recv(C + IDX(offsets[i], 0, n), sizes[i], row,
i, 0, MPI_COMM_WORLD, &status);
assert(ret == MPI_SUCCESS);
}
} else {
ret = MPI_Send(C, sizes[comm_id], row, 0, 0, MPI_COMM_WORLD);
assert(ret == MPI_SUCCESS);
}
t2 = MPI_Wtime();
/*
* Mamy w pełni złożoną macierz C. Możemy sprawdzić wyniki, zapisać je
* gdzieś, cokolwiek.
*/
if (comm_id == 0) {
#ifdef DEBUG
/* TDD */
check_correctness(A, B, C, n);
#endif
printf("Total time: %lf[s]\n", t2 - t1);
}
free(offsets);
free(sizes);
free(A);
free(B);
free(C);
MPI_Finalize();
return (0);
}
Zadanie: przerobić to na dowolne rozmiary macierzy (nie tylko kwadraty), zastanowić się nad modelem scatter/gather, zastanowić się nad podobną metodą, gdzie dystrybuuje się kolumny macierzy B, całą macierz A oraz otrzymujemy kolumny macierzy C, zastanowić się, która metoda (wierszowa, czy kolumnowa) jest lepsza i dlaczego.
Jak sobie odpowiesz na wszystkie pytania powyżej będziesz dobrze rozumiał na czym polega mnożenie macierzy, nauczysz się MPI i może nawet błyśniesz przed prowadzącym.
