Dobra, po kolei.
- Myślałem, że jak się bierzesz za mnożenie macierzy to przynajmniej znasz podstawy. BLAS to "Basic Linear Algebra System" i w nim masz procedurę mnożenia macierzy. Na typach zmiennoprzecinkowych, bo ja osobiście nie znam żadnej aplikacji, która robiłaby mnożenie przez siebie dwóch macierzy na typie innym niż zmiennoprzecinkowy. Mnie na studiach ostrzegali, żeby nie robić samemu mnożenia macierzy, bo raz, że będzie wolno, dwa, że zmarnuję całe tygodnie na klepanie rzeczy już dawno poprawnie zaimplementowanych. Ciekawe, że w przypadku web developmentu nikt nie zaczyna od klepania ręcznej obsługi gniazdek BSD, tylko wszyscy używają porządnych bibliotek / frameworków. Za to każdy, kto widzi macierz od razu zaczyna od klepania mnożenia macierzy własnoręcznie. Fuck logic.
- "Co do kodu jest napisany tak dobrze , że ciężko mi jest się połapać jakie kroki są wykonywane po kolei." No to witamy w dziedzinie, gdzie kod rozumieją tylko specjaliści. Nie rozumiesz kodu - zajmij się czymś innym.
Masz tutaj poniżej wyczyszczony, dobrze skomentowany kod.
/*
* 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.