[es] - Množenje matrica - pomoć oko Strassena

Iljo
Marko Ilić
Varazdin-Zagreb

Član broj: 65402
Poruke: 35
*.adsl.net.t-com.hr.

Profil

Množenje matrica - pomoć oko Strassena

^{09.06.2007. u 00:35 - pre 205 meseci}

http://en.wikipedia.org/wiki/Strassen_algorithm
Evo mog koda napisanog u C-u, na Linuxu a imam ako treba i za Windowse.

Code:

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <sys/time.h>

#define DIM     512        //dimenzija matrice
#define NISTA    2
#define RANDOM    1
#define NULAMA    0
#define DODAJ    1
#define NOVA    0

void napravi_matricu (int ***M, int dim_matrice, int nacin);
void ispisi_matricu (int **M);
void dealociraj_matricu (int **M, int dim_matrice);
void suma_matrica (int **A, int **B, int **C, int dim_matrice, int nacin);    //C = A + B
void suma1_matrica (int **A, int i1, int j1, int **B, int i2, int j2, int **C, int dim_matrice, int nacin);
void razlika_matrica (int **A, int **B, int **C, int dim_matrice, int nacin);    //C = A - B
void razlika1_matrica (int **A, int i1, int j1, int **B, int i2, int j2, int **C, int dim_matrice, int nacin);
void strassen_matrice (int **A, int **B, int **C, int dim_matrice);

int main (void)
{
    int **A, **B, **C;
    double start, end;
    struct timeval vrijeme;

    srand( 0 );
    napravi_matricu( &A, DIM, RANDOM );
    napravi_matricu( &B, DIM, RANDOM );
    napravi_matricu( &C, DIM, NULAMA );

    gettimeofday( &vrijeme, NULL );
    start = vrijeme.tv_sec + (vrijeme.tv_usec/1000000.0);

    strassen_matrice( A, B, C, DIM );

    gettimeofday( &vrijeme, NULL );
    end = vrijeme.tv_sec + (vrijeme.tv_usec/1000000.0);

/*    printf( "--- A ---\n" ); ispisi_matricu( A );
    printf( "--- B ---\n" ); ispisi_matricu( B );
    printf( "--- C ---\n" ); ispisi_matricu( C );
*/    dealociraj_matricu( A, DIM );
    dealociraj_matricu( B, DIM );
    dealociraj_matricu( C, DIM );
    printf( "%d - Vrijeme: %fs\n", DIM, end-start );

    return 0;
}

void napravi_matricu (int ***M, int dim_matrice, int nacin)
{
    int i, j;

    (*M) = (int **)malloc( dim_matrice * sizeof(int *) );
    if ((*M) == NULL) {
        printf( "Greska: Alokacija memorije nije uspjela!\n" );
        exit( -1 );
    }
    for (i = 0; i < dim_matrice; i++) {
        (*M)[i] = (int *)malloc( dim_matrice * sizeof(int) );
        if ((*M)[i] == NULL) {
            printf( "Greska 2: Alokacija memorije nije uspjela!\n" );
            exit( -1 );
        }
    }

    if (nacin == RANDOM)
        for (i = 0; i < dim_matrice; i++)
            for (j = 0; j < dim_matrice; j++)
                (*M)[i][j] = rand() % 101;
    else if (nacin == NULAMA)
        for (i = 0; i < dim_matrice; i++)
            for (j = 0; j < dim_matrice; j++)
                (*M)[i][j] = 0;

    return;
}

void ispisi_matricu (int **M)
{
    int i, j;

    for (i = 0; i < DIM; i++) {
        for (j = 0; j < DIM; j++)
            printf( "%7d", M[i][j] );
        printf( "\n" );
    }

    return;
}

void dealociraj_matricu (int **M, int dim_matrice)
{
    int i;

    for (i = 0; i < dim_matrice; i++)
        free( M[i] );
    free( M );

    return;
}

void suma_matrica (int **A, int **B, int **C, int dim_matrice, int nacin)
{
    int i, j;

    if (nacin == DODAJ)
        for (i = 0; i < dim_matrice; i++)
            for (j = 0; j < dim_matrice; j++)
                C[i][j] += A[i][j] + B[i][j];
    else
        for (i = 0; i < dim_matrice; i++)
            for (j = 0; j < dim_matrice; j++)
                C[i][j] = A[i][j] + B[i][j];

    return;
}

void razlika_matrica (int **A, int **B, int **C, int dim_matrice, int nacin)
{
    int i, j;

    if (nacin == DODAJ)
        for (i = 0; i < dim_matrice; i++)
            for (j = 0; j < dim_matrice; j++)
                C[i][j] += A[i][j] - B[i][j];
    else
        for (i = 0; i < dim_matrice; i++)
            for (j = 0; j < dim_matrice; j++)
                C[i][j] = A[i][j] - B[i][j];

    return;
}

void suma1_matrica (int **A, int i1, int j1, int **B, int i2, int j2, int **C, int dim_matrice, int nacin)
{
    int i, j;

    if (nacin == DODAJ)
        for (i = 0; i < dim_matrice; i++)
            for (j = 0; j < dim_matrice; j++)
                C[i][j] += A[i1+i][j1+j] + B[i2+i][j2+j];
    else
        for (i = 0; i < dim_matrice; i++)
            for (j = 0; j < dim_matrice; j++)
                C[i][j] = A[i1+i][j1+j] + B[i2+i][j2+j];

    return;
}

void razlika1_matrica (int **A, int i1, int j1, int **B, int i2, int j2, int **C, int dim_matrice, int nacin)
{
    int i, j;

    if (nacin == DODAJ)
        for (i = 0; i < dim_matrice; i++)
            for (j = 0; j < dim_matrice; j++)
                C[i][j] += A[i1+i][j1+j] - B[i2+i][j2+j];
    else
        for (i = 0; i < dim_matrice; i++)
            for (j = 0; j < dim_matrice; j++)
                C[i][j] = A[i1+i][j1+j] - B[i2+i][j2+j];

    return;
}

void strassen_matrice (int **A, int **B, int **C, int dim_matrice)
{
    int     i, j, k,
        dim_submatrice = dim_matrice/2,
        **Temp1, **Temp2, **Temp3,
        **M1, **M2, **M3, **M4, **M5, **M6, **M7;

    if (dim_matrice <= 1) {
        C[0][0] = A[0][0] * B[0][0];
        return;
    }

    napravi_matricu( &Temp1, dim_submatrice, NISTA ); napravi_matricu( &Temp2, dim_submatrice, NISTA );
    napravi_matricu( &Temp3, dim_submatrice, NISTA );
    napravi_matricu( &M1, dim_submatrice, NISTA ); napravi_matricu( &M2, dim_submatrice, NISTA );
    napravi_matricu( &M3, dim_submatrice, NISTA ); napravi_matricu( &M4, dim_submatrice, NISTA );
    napravi_matricu( &M5, dim_submatrice, NISTA ); napravi_matricu( &M6, dim_submatrice, NISTA );
    napravi_matricu( &M7, dim_submatrice, NISTA );

    suma1_matrica( A, 0, 0, A, dim_submatrice, dim_submatrice, Temp1, dim_submatrice, NOVA );
    suma1_matrica( B, 0, 0, B, dim_submatrice, dim_submatrice, Temp2, dim_submatrice, NOVA );
    strassen_matrice( Temp1, Temp2, M1, dim_submatrice );

    suma1_matrica( A, dim_submatrice, 0, A, dim_submatrice, dim_submatrice, Temp1, dim_submatrice, NOVA );
    for (i = 0; i < dim_submatrice; i++)
        for (j = 0; j < dim_submatrice; j++)
            Temp3[i][j] = B[i][j];
    strassen_matrice( Temp1, Temp3, M2, dim_submatrice );

    razlika1_matrica( B, 0, dim_submatrice, B, dim_submatrice, dim_submatrice, Temp2, dim_submatrice, NOVA );
    for (i = 0; i < dim_submatrice; i++)
        for (j = 0; j < dim_submatrice; j++)
            Temp3[i][j] = A[i][j];
    strassen_matrice( Temp3, Temp2, M3, dim_submatrice );

    razlika1_matrica( B, dim_submatrice, 0, B, 0, 0, Temp2, dim_submatrice, NOVA );
    for (i = dim_submatrice; i < dim_matrice; i++)
        for (j = dim_submatrice; j < dim_matrice; j++)
            Temp3[(i-dim_submatrice)][(j-dim_submatrice)] = A[i][j];
    strassen_matrice( Temp3, Temp2, M4, dim_submatrice );

    suma1_matrica( A, 0, 0, A, 0, dim_submatrice, Temp1, dim_submatrice, NOVA );
    for (i = dim_submatrice; i < dim_matrice; i++)
        for (j = dim_submatrice; j < dim_matrice; j++)
            Temp3[(i-dim_submatrice)][(j-dim_submatrice)] = B[i][j];
    strassen_matrice( Temp1, Temp3, M5, dim_submatrice );

    razlika1_matrica( A, dim_submatrice, 0, A, 0, 0, Temp1, dim_submatrice, NOVA );
    suma1_matrica( B, 0, 0, B, 0, dim_submatrice, Temp2, dim_submatrice, NOVA );
    strassen_matrice( Temp1, Temp2, M6, dim_submatrice );

    razlika1_matrica( A, 0, dim_submatrice, A, dim_submatrice, dim_submatrice, Temp1, dim_submatrice, NOVA );
    suma1_matrica( B, dim_submatrice, 0, B, dim_submatrice, dim_submatrice, Temp2, dim_submatrice, NOVA );
    strassen_matrice( Temp1, Temp2, M7, dim_submatrice );


    razlika_matrica( M4, M5, Temp1, dim_submatrice, NOVA );
    suma_matrica( M1, M7, Temp1, dim_submatrice, DODAJ );
    for (i = 0; i < dim_submatrice; i++)
        for (j = 0; j < dim_submatrice; j++){
            C[i][j] = Temp1[i][j];
        }
    suma_matrica( M3, M5, Temp1, dim_submatrice, NOVA );
    for (i = 0; i < dim_submatrice; i++)
        for (j = 0; j < dim_submatrice; j++){
            C[i][(j+dim_submatrice)] = Temp1[i][j];
        }
    suma_matrica( M2, M4, Temp1, dim_submatrice, NOVA );
    for (i = 0; i < dim_submatrice; i++)
        for (j = 0; j < dim_submatrice; j++){
            C[(i+dim_submatrice)][j] = Temp1[i][j];
        }
    razlika_matrica( M1, M2, Temp1, dim_submatrice, NOVA );
    suma_matrica( M3, M6, Temp1, dim_submatrice, DODAJ );
    for (i = 0; i < dim_submatrice; i++)
        for (j = 0; j < dim_submatrice; j++){
            C[(i+dim_submatrice)][(j+dim_submatrice)] = Temp1[i][j];
        }

    dealociraj_matricu( Temp1, dim_submatrice ); dealociraj_matricu( Temp2, dim_submatrice );
    dealociraj_matricu( Temp3, dim_submatrice );
    dealociraj_matricu( M1, dim_submatrice ); dealociraj_matricu( M2, dim_submatrice );
    dealociraj_matricu( M3, dim_submatrice ); dealociraj_matricu( M4, dim_submatrice );
    dealociraj_matricu( M5, dim_submatrice ); dealociraj_matricu( M6, dim_submatrice );
    dealociraj_matricu( M7, dim_submatrice );
    return;
}

E sad on je puno sporiji nego kad napišem običan algoritam za množenje, a trebao bi biti brži.
Znam da bi trebao biti brži za velike matrice, ali čak i kad stavim da mi je red matrice 2048 on je još uvijek sporiji.

Radio sam neka poboljšanja u programu da bude malo brži od ovog al je i taj puno sporiji od običnog. (nisam ga stavio na forum jer je teško čitljiv)
Probao sam maknut pozivanja funkcije razlika1_matrica i suma1_matrica u funkciji strassen_matrice, i stavit taj kod direkno tamo gdje se one pozivaju, al poboljšanje je gotovo nikakvo.
Još sam probao srezat doljnje rekurzije tako da sam ovaj dio koda:

Code:

    if (dim_matrice <= 1) {
        C[0][0] = A[0][0] * B[0][0];
        return;
    }

zamjenio sa:

Code:

    if (dim_matrice <= 64) {
        for (i = 0; i < dim_matrice; i++)
            for (k = 0; k < dim_matrice; k++)
                for (j = 0; j < dim_matrice; j++)
                    C[i][j] = 0;
        for (i = 0; i < dim_matrice; i++)
            for (k = 0; k < dim_matrice; k++)
                for (j = 0; j < dim_matrice; j++)
                    C[i][j] += A[i][k] * B[k][j];
        return;
    }

i sada je strassen brži (zapravo nije ni pravi strassen sada), al tek kada je red matrice >= 2048. (brži je za 15 sec).

Bio bih zahvalan kad bi mi netko pogledao program i rekao kako da ga ubrzam, ili dao neki svoj program za Strassena napisanog u C/C++, može i nekom drugom sličnom jeziku.

Odgovor na temu