Mnożenie macierzy z wykorzystaniem MPI

0

Witam!
Napisałem sobie o to taki prosty programik do mnożenia macierzy.


#include <stdio.h>
#include <stdlib.h> //funkcja rand()
#include <time.h> // zarzadzanie czasem
//#include "mpi.h"
int main(int argc, char *argv[])
{
        int n, m, j, k; //deklaracja wymiarow macierzy

        printf("Podaj wymiary dwoch macierzy:"); // tylko 3 wymiary bo liczba kolumn pierwszej musi byc rowna liczbie  wierszy drugiej
        scanf("%d%d%d", &n, &m, &k);

        j = m ;

        int macierzA[n][m];
        int macierzB[j][k];
        int wynik[n][k];

        printf("Pierwsza macierz: n=%d, m=%d\n", n, m);
        printf("Druga macierz: j= %d, k=%d\n", j, k);

        printf("\n");
        printf("Pierwsza macierz: \n");

        int i, p, u; //deklaracja licznikow
        for (i = 0; i < n; i++) //wpisywanie losowych wartosci wierszami
                for (p = 0;p < m; p++)
                       macierzA[i][p] = rand() % 10; // generowanie wartosci od 0 do 9

        for (i = 0; i < n; i++) // wypisywanie tablicy wierszami
        {
                printf("\n");
                for (p = 0; p < m; p++)
                        printf("%d\t", macierzA[i][p]);
        }

        printf("\n\n");
        printf("Druga macierz:\n");

        for (i = 0; i < j; i++) //wpisywanie losowych wartosci wierszami
                for (p = 0; p < k; p++)
                        macierzB[i][p] = rand() % 10; // generuje wartosci od 0 do 9

        for (i = 0;i < j; i++) //wypisywanie tablicy wierszami
        {
                printf("\n");
        for (p = 0; p < k; p++)
                printf("%d\t", macierzB[i][p]);
        }

        printf("\n\n");

        for (i = 0; i < n; i++) // mnozenie macierzy
        for (p = 0; p < k; p++)
        {
                wynik[i][p] = 0; // przygotowanie macierzy wynikowej "wynik"
                for (u= 0; u < j; u++)
                        wynik[i][p]= wynik[i][p] + macierzA[i][u] * macierzB[u][p]; // mnozenie
        }

        // wypisanie macierzy
        printf("Macierz wynikowa: \n");
        for (i = 0; i < n; i++)
        {
                printf("\n");
                for (p = 0; p < k; p++)
                        printf("%d\t", wynik[i][p]);
        }
        return 0;
}

Problem w tym , że chciałbym teraz przerobić go tak aby wykonywał się przy użyciu załóżmy 4 PC czyli przy użyciu MPI.
Rozumiem ,że musiał bym zacząć od podzielenia macierzy wynikowej na 4 części tak aby każdy z PC otrzymał do wyliczenia 1/4 macierzy. Następnie każdy PC wylicza daną część macierzy
zwracają one wynik i następuje połączenie. Ze strony teoretycznej wiem jak to powinno wyglądać. Problem w tym ,że nie bardzo wiem jak skonstruować pętlę, która dzieliła i liczyła by macierz o rozmiarze podawanym dynamicznie ( weźmy na początek tylko macierze kwadratowe typu 2x2, 4x4, 8x8, 16x16 itp). Prosił bym o pomoc w tej sprawie :-)

0

Zrób inaczej, ponieważ MPI ma mądre metody do wysyłki i odbierania danych, zwłaszcza w postaci macierzowej, to za ich pomocą odpowiednio roześlij dane, a w poszczególnych procesach odbierz je i przetwarzaj jak by to były zupełnie niezależne macierze. Dopiero przy odbiorze bloków poskładaj je tak, aby otrzymać prawidłowy wynik.

P.s. macierz typu int? Robisz równoległe przetwarzanie plików graficznych czy jak :>?

Edit:
Jeśli chodzi o samo dzielenie, to nie musisz robić żadnej pętli, po prostu dzielisz liczbę kolumn (n) i wierszy macierzy (m) na c (c = liczba cpu) i do pierwszych (n % c) kolumn i (m % c) wierszy dodajesz 1, aby zrobić taki mały load-balancing.

0

Ależ skąd.. :) jestem studentem i miałem ostatnio wykład z MPI. A korzystając z tego, że mam dostęp do sieci uczelnianej i skonfigurowanego MPI to chciałem sprawdzić jak to działa i napisać sprawko na zaliczenie z danego tematu :) Czyli rozumiem , że powinienem rozesłać całą macierz A i B do każdego z procesu, następnie zrobić tak aby każdy z procesów obliczył
daną część macierzy wynikowej . Na końcu składam macierze w jedną całość ?

Example:


 macierzA                             macierz B
 a11  a12  a13  a14            b11  b12  b13  b14
 a21  a22  a23  a24            b21  b22  b23  b24
 a31  a32  a33  a34            b31  b32  b33  b34
 a41  a42  a43  a44            b41  b42  b43  b44

obydwie całe wysyłam do procesów o nr: 0 1 2 3

proces 0 oblicza lewa górna cześć macierzy wynikowej
ab11  ab12    
ab21  ab22

proces 1
ab13  ab14
ab23  ab24

proces 2
ab31  ab32
ab41  ab42

proces 3
ab33  ab34
ab43  ab44

Po obliczeniu zwracają wyliczone macierze
i składam je w jedną macierz wynikowa

wynikowa:
ab11  ab12  ab13  ab14
ab21  ab22  ab23  ab24
ab31  ab32  ab33  ab34
ab41  ab42  ab43  ab44

Dobrze rozumiem ? Czy zupełnie chodzi o coś innego ?

0

Możesz wysłać całą macierz, a możesz użyć MPI_Type_vector do wysłania do poszczególnych procesów tylko wierszy i kolumn, które są potrzebne do wyliczenia danego fragmentu macierzy C.

Edit: tak, dobrze rozumiesz. Polecam Ci jeszcze dwa usprawnienia. Zmień typ macierzy na double i użyj BLAS level 3 (dokładniej: dgemm).

Edit2: Używasz błędnej definicji mnożenia macierzy o ile zauważyłem. Żeby wyliczyć block C[offset:offset + size][offset:offset + size] musisz mieć wiersze A[offset:offset+size][0:size] oraz kolumny B[0:size][offset:offset+size].

2

Na szybkiego napisałem Ci PoC-a jak takie mnożenie mogłoby wyglądać.
Ja po prostu rozbiłem macierz A wierszami tak, żeby każdy z procesów miał mniej więcej po tyle samo pracy.
Teoretycznie dałoby się to jescze bardzo mocno usprawnić (broadcast macierzy B, nieblokujący send, itp.), ale na początek wystarczy.

#include <assert.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include <cblas.h>
#include <mpi.h>

/* Row-major order */
#define IDX(i, j, n)    (i * n + j)

int
main(int argc, char **argv)
{
    double *A, *B, *C;

    int i, j, n, ret;
    int comm_size, comm_id;

    int *offsets, *sizes;

    double t1, t2;

    MPI_Datatype row;
    MPI_Status status;

    MPI_Init(&argc, &argv);

    t1 = t2 = 0;
    n = 8;
    if (argc > 1)
        n = atoi(argv[1]);

#ifdef DEBUG
    printf("n=%d\n", n);
#endif
    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);

    ret = MPI_Type_vector(1, n, n, MPI_DOUBLE, &row);
    assert(ret == MPI_SUCCESS);
    ret = MPI_Type_commit(&row);
    assert(ret == MPI_SUCCESS);

    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]);

    if (comm_id == 0) {
        ret = posix_memalign((void **)&A, 64, n * n * sizeof(*A));
    } else {
        ret = posix_memalign((void **)&A, 64, sizes[comm_id] * n *
            sizeof(*A));
    }
    assert(ret == 0);

    ret = posix_memalign((void **)&B, 64, n * n * sizeof(*B));
    assert(ret == 0);

    if (comm_id == 0) {
        ret = posix_memalign((void **)&C, 64, n * n * sizeof(*C));
    } else {
        ret = posix_memalign((void **)&C, 64, sizes[comm_id] * n *
            sizeof(*C));
    }
    assert(ret == 0);

    if (comm_id == 0)
        memset(C, 0, n * n * sizeof(*C));

    if (comm_id == 0) {

        for (i = 0; i < n; i++) {
            for (j = 0; j < n; j++) {
                A[IDX(i, j, n)] = 1.0f * rand() / RAND_MAX;
                B[IDX(i, j, n)] = 1.0f * rand() / RAND_MAX;
            }
        }
        if (comm_id == 0)
            t1 = MPI_Wtime();

        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);

            ret = MPI_Send(B, n, 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);

        ret = MPI_Recv(B, n, row, 0, 0, MPI_COMM_WORLD, &status);
        assert(ret == MPI_SUCCESS);
    }

    /*
     * Given process computes block of C:
     * 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);

    cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans,
        sizes[comm_id], /* A is sizes[comm_id] x n */
        n, n, /* B is n x n */
        1.0,
        A, n,
        B, n,
        0.0,
        C, n);

    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);
    }

    if (comm_id == 0)
        t2 = MPI_Wtime();

    /*
     * In C we have right now fully assembled matrix. Check for
     * correctness.
     */

#ifdef DEBUG
    if (comm_id == 0) {
        double *D;

        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)]);

                printf("%lf ", C[IDX(i, j, n)]);
            }
            printf("\n");
        }
        free(D);
    }
#endif

    if (comm_id == 0)
        printf("Total time: %lf[s]\n", t2 - t1);

    free(A);
    free(B);
    free(C);
    MPI_Finalize();
    return (0);
}
0

Dziękuję ślicznie :) muszę sobie ten kod na spokojnie przeanalizować :)

0

Kurcze podczas kompilacji tego kodu przez:


mpicc -o matrix matrix.c

wyrzuca błąd :


/tmp/ccZuhMF0.o: In function `main':
matrix.c:(.text+0x733): undefined reference to `cblas_dgemm'
collect2: error: ld returned 1 exit status

Pewnie z powodu braku biblioteki na komputerach uczelni.

 #include <cblas.h>

Co do kodu jest napisany tak dobrze , że ciężko mi jest się połapać jakie kroki są wykonywane po kolei.

0

Dobra, po kolei.

  1. 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.
  2. "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.

0
cuadam1992 napisał(a):

Witam!
Napisałem sobie o to taki prosty programik do mnożenia macierzy.

[...]
int macierzA[n][m];

    for (i = 0; i < n; i++) //wpisywanie losowych wartosci wierszami
            for (p = 0;p < m; p++)
                   macierzA[i][p] = rand() % 10; // generowanie wartosci od 0 do 9

[...]



@cuadam1992
Na moje inżynierskie oko, tam jeszcze gdzieśkolwiek powinien się malloc po drodze znaleźć.

1 użytkowników online, w tym zalogowanych: 0, gości: 1, botów: 0