Panel użytkownika
Nazwa użytkownika:
Hasło:
Nie masz jeszcze konta?

Algorytm Strassena - rekurencja w wielu wywołaniach

Ostatnio zmodyfikowano 2016-12-07 14:42
Autor Wiadomość
Korybut
Temat założony przez niniejszego użytkownika
Algorytm Strassena - rekurencja w wielu wywołaniach
» 2016-12-05 10:59:57
Cześć!

Ostatnio realizuje mnożenie "dużych" macierzy 4k na 4k.. postanowiłem więc wykorzystać algorytm Strassena aby uciąć jedno mnożenie. Program działał prawidłowo, zwracał poprawny wynik, ale problem zaczął się gdy spróbowałem rekurencji. Nie wiem czy może gdzieś pamięć ginie, lub wyniki przechowywane z wyższego poziomu rekurencji są nadpisywane, ale macierz końcowa to zbiór abstrakcyjnych liczb. Może zamiast opisywać dokładnie problem to zaprezentuje kod - na pewno więcej wyjaśni.

header.h
C/C++
#ifndef HEADER_H_
#define HEADER_H_
#include <iostream>
#include <time.h>
#include <stdlib.h>

using namespace std;

int ** allocMatrixEX( int, int );
void deleteMatrixEX( int **&, int size );

void generateMX( int **& matrix, int size );
void copyMX( int **& matrix, int **& MX0, int size );

void subdivide( int ** MX0, int **& matrix, int i, int j, int size );

int ** add( int **& MX0, int **& MX1, int **& MX2, int size );
int ** subtract( int **& MX0, int **& MX1, int **& MX2, int size );
int ** multiply( int **& MX0, int **& MX1, int **& MX2, int size );
int ** strassen( int **& MXC, int **& MXA, int **& MXB, int value );
void summary( int **& MOC, int **& C11, int **& C12, int **& C21, int **& C22, int size );

void showM( int **& MX0, int size );

#endif

main.cpp
C/C++
#include <iostream>
#include <time.h>
#include <stdlib.h>
#include "header.h"

int main( void ) {
    srand( time( NULL ) );
    int value = 4096;
   
    // tu powinna zacząć się funkcja strassen
    int size = value;
   
    int ** MA = allocMatrixEX( size, size );
    int ** MB = allocMatrixEX( size, size );
    int ** MC = allocMatrixEX( size, size );
   
    generateMX( MA, size );
    generateMX( MB, size );
    //showM(MA, size);
    // showM(MB, size);
    // system("pause");
   
    MC = strassen( MC, MA, MB, size );
    // showM(MC, size);
   
    deleteMatrixEX( MA, size );
    deleteMatrixEX( MB, size );
    deleteMatrixEX( MC, size );
   
    return 0;
}

definition.cpp
C/C++
#include <iostream>
#include <time.h>
#include <stdlib.h>
#include "header.h"

// ALOKACJA PAMIĘCI DLA MACIERZY
int ** allocMatrixEX( int height, int width ) {
    int ** matrix = new int *[ height ];
    for( int i = 0; i < height; ++i )
         matrix[ i ] = new int[ width ];
   
    return matrix;
}

// USUWANIE PAMIĘCI PO MACIERZY
void deleteMatrixEX( int **& matrix, int height ) {
    for( int i = 0; i < height; ++i )
         delete matrix[ i ];
   
    delete matrix;
}

// GENEROWANIE MACIERZY
void generateMX( int **& matrix, int size ) {
    for( int i = 0; i < size; i++ )
    for( int j = 0; j < size; j++ )
         matrix[ i ][ j ] = rand() % 3;
   
}

// KOPIOWANIE MACIERZY
void copyMX( int **& matrix, int **& MX0, int size ) {
    for( int i = 0; i < size; i++ )
    for( int j = 0; j < size; j++ )
         matrix[ i ][ j ] = MX0[ i ][ j ];
   
}

// MNOŻENIE MACIERZY
int ** multiply( int **& MX0, int **& MX1, int **& MX2, int size ) {
    for( int y = 0; y < size / 2; y++ ) {
        for( int j = 0; j < size / 2; j++ ) {
            for( int i = 0; i < size / 2; i++ ) {
                MX0[ y ][ j ] += MX1[ y ][ i ] * MX2[ i ][ j ];
            }
        }
        //cout << "mnozenie kolumny " << y << " macierzy " << MX0 << endl;
    }
    return MX0;
}

// DODAWANIE MACIERZY
int ** add( int **& MX0, int **& MX1, int **& MX2, int size ) {
    for( int i = 0; i < size / 2; i++ ) {
        for( int j = 0; j < size / 2; j++ )
             MX0[ i ][ j ] = MX1[ i ][ j ] + MX2[ i ][ j ];
        //cout << "dodawanie kolumny " << i << " macierzy " << MX0 << endl;
    }
    return MX0;
}

// ODEJMOWANIE MACIERZY
int ** subtract( int **& MX0, int **& MX1, int **& MX2, int size ) {
    for( int i = 0; i < size / 2; i++ ) {
        for( int j = 0; j < size / 2; j++ )
             MX0[ i ][ j ] = MX1[ i ][ j ] - MX2[ i ][ j ];
        //cout << "odejmowanie kolumny " << i << " macierzy " << MX0 << endl;
    }
    return MX0;
}

// PODZIAŁ GŁOWNYCH MACIERZY NA ĆWIARTKI
void subdivide( int ** MX0, int **& matrix, int i, int j, int size ) {
    int x = i;
    for( int a = 0; a < size; a++ )
    {
        int y = j;
        for( int b = 0; b < size; b++ )
        {
            MX0[ a ][ b ] = matrix[ x ][ y ];
            y++;
        }
        x++;
    }
}

// KOŃCOWE SUMOWANIE MACIERZY WYNIKOWEJ
void summary( int **& MOC, int **& c11, int **& c12, int **& c21, int **& c22, int size ) {
    int x = 0;
    for( int a = 0; a < size / 2; a++ )
    {
        int y = 0;
        for( int b = 0; b < size / 2; b++ )
        {
            MOC[ a ][ b ] = c11[ x ][ y ];
            y++;
        }
        x++;
    }
    x = 0;
    for( int a = 0; a < size / 2; a++ )
    {
        int y = 0;
        for( int b = size / 2; b < size; b++ )
        {
            MOC[ a ][ b ] = c12[ x ][ y ];
            y++;
        }
        x++;
    }
    x = 0;
    for( int a = size / 2; a < size; a++ )
    {
        int y = 0;
        for( int b = 0; b < size / 2; b++ )
        {
            MOC[ a ][ b ] = c21[ x ][ y ];
            y++;
        }
        x++;
    }
    x = 0;
    for( int a = size / 2; a < size; a++ )
    {
        int y = 0;
        for( int b = size / 2; b < size; b++ )
        {
            MOC[ a ][ b ] = c22[ x ][ y ];
            y++;
        }
        x++;
    }
}

// WYŚWIETLANIE MACIERZY
void showM( int **& MX0, int size ) {
    for( int i = 0; i < size; i++ ) {
        for( int j = 0; j < size; j++ ) {
            cout << MX0[ i ][ j ] << " | ";
        }
        cout << endl;
    }
    cout << "**********************************************\n";
}

// FUNKCJA STRASSEN

int ** strassen( int **& MXC, int **& MXA, int **& MXB, int value ) {
    int size = value;
   
    /* ALOKACJA PAMIĘCI NA POTRZEBY MACIERZY OPERACYJNYCH */
    int ** matrixA = allocMatrixEX( size, size );
    int ** matrixB = allocMatrixEX( size, size );
    int ** matrixC = allocMatrixEX( size, size );
   
    int ** A11 = allocMatrixEX( size / 2, size / 2 );
    int ** A12 = allocMatrixEX( size / 2, size / 2 );
    int ** A21 = allocMatrixEX( size / 2, size / 2 );
    int ** A22 = allocMatrixEX( size / 2, size / 2 );
    int ** B11 = allocMatrixEX( size / 2, size / 2 );
    int ** B12 = allocMatrixEX( size / 2, size / 2 );
    int ** B21 = allocMatrixEX( size / 2, size / 2 );
    int ** B22 = allocMatrixEX( size / 2, size / 2 );
   
    int ** M1 = allocMatrixEX( size / 2, size / 2 );
    int ** M2 = allocMatrixEX( size / 2, size / 2 );
    int ** M3 = allocMatrixEX( size / 2, size / 2 );
    int ** M4 = allocMatrixEX( size / 2, size / 2 );
    int ** M5 = allocMatrixEX( size / 2, size / 2 );
    int ** M6 = allocMatrixEX( size / 2, size / 2 );
    int ** M7 = allocMatrixEX( size / 2, size / 2 );
    /* KONIEC ALOKACJI */
   
    /* PRZYPISANIE WARTOŚCI DO MACIERZY I ICH PODZIAŁ NA ĆWIARTKI */
   
    copyMX( matrixA, MXA, size );
    subdivide( A11, matrixA, 0, 0, size / 2 );
    subdivide( A12, matrixA, 0, size / 2, size / 2 );
    subdivide( A21, matrixA, size / 2, 0, size / 2 );
    subdivide( A22, matrixA, size / 2, size / 2, size / 2 );
   
    //showM(A11, size/2);
    //showM(A12, size/2);
    //showM(A21, size/2);
    //showM(A22, size/2);
    // system("pause");
    deleteMatrixEX( matrixA, size );
   
    copyMX( matrixB, MXB, size );
    subdivide( B11, matrixB, 0, 0, size / 2 );
    subdivide( B12, matrixB, 0, size / 2, size / 2 );
    subdivide( B21, matrixB, size / 2, 0, size / 2 );
    subdivide( B22, matrixB, size / 2, size / 2, size / 2 );
   
    //showM(B11, size/2);
    //showM(B12, size/2);
    //showM(B21, size/2);
    //showM(B22, size/2);
    // system("pause");
    deleteMatrixEX( matrixB, size );
    /* KONIEC PRZYPISYWANIA I PODZIAŁU MACIERZY WEJŚCIOWYCH */
   
    /* OBLICZANIE MACIERZY OPERACYJNYCH */
    // jeśli ma dzielić dalej to mnożenie odbywa się za pomocą funkcji strassen(), w innym wypadku za pomocą funkcji multiply() //
   
    int ** temp0 = allocMatrixEX( size / 2, size / 2 );
    int ** temp1 = allocMatrixEX( size / 2, size / 2 );
   
    if( size > 256 ) {
        //krok 1
        temp0 = add( temp0, A11, A22, size );
        temp1 = add( temp1, B11, B22, size );
        M1 = strassen( M1, temp0, temp1, size / 2 );
        //krok 2
        temp0 = add( temp0, A21, A22, size );
        M2 = strassen( M2, temp0, B11, size / 2 );
        //krok 3
        temp1 = subtract( temp1, B12, B22, size );
        M3 = strassen( M3, A11, temp1, size / 2 );
        //krok 4
        temp1 = subtract( temp1, B21, B11, size );
        M4 = strassen( M4, A22, temp1, size / 2 );
        //krok 5
        temp1 = add( temp1, A11, A12, size );
        M5 = strassen( M5, temp1, B22, size / 2 );
        //krok 6
        temp0 = subtract( temp0, A21, A11, size );
        temp1 = add( temp1, B11, B12, size );
        M6 = strassen( M6, temp0, temp1, size / 2 );
        //krok 7
        temp0 = subtract( temp0, A12, A22, size );
        temp1 = add( temp1, B21, B22, size );
        M7 = strassen( M7, temp0, temp1, size / 2 );
    } else {
        //krok 1
        temp0 = add( temp0, A11, A22, size );
        temp1 = add( temp1, B11, B22, size );
        M1 = multiply( M1, temp0, temp1, size / 2 );
        //krok 2
        temp0 = add( temp0, A21, A22, size );
        M2 = multiply( M2, temp0, B11, size / 2 );
        //krok 3
        temp1 = subtract( temp1, B12, B22, size );
        M3 = multiply( M3, A11, temp1, size / 2 );
        //krok 4
        temp1 = subtract( temp1, B21, B11, size );
        M4 = multiply( M4, A22, temp1, size / 2 );
        //krok 5
        temp1 = add( temp1, A11, A12, size );
        M5 = multiply( M5, temp1, B22, size / 2 );
        //krok 6
        temp0 = subtract( temp0, A21, A11, size );
        temp1 = add( temp1, B11, B12, size );
        M6 = multiply( M6, temp0, temp1, size / 2 );
        //krok 7
        temp0 = subtract( temp0, A12, A22, size );
        temp1 = add( temp1, B21, B22, size );
        M7 = multiply( M7, temp0, temp1, size / 2 );
    }
   
    /* KONIEC OBLICZEŃ OPERACYJNYCH */
   
    deleteMatrixEX( A11, size / 2 );
    deleteMatrixEX( A12, size / 2 );
    deleteMatrixEX( A21, size / 2 );
    deleteMatrixEX( A22, size / 2 );
    deleteMatrixEX( B11, size / 2 );
    deleteMatrixEX( B12, size / 2 );
    deleteMatrixEX( B21, size / 2 );
    deleteMatrixEX( B22, size / 2 );
    deleteMatrixEX( temp0, size / 2 );
    deleteMatrixEX( temp1, size / 2 );
   
    int ** C11 = allocMatrixEX( size / 2, size / 2 );
    int ** C12 = allocMatrixEX( size / 2, size / 2 );
    int ** C21 = allocMatrixEX( size / 2, size / 2 );
    int ** C22 = allocMatrixEX( size / 2, size / 2 );
    int ** temp2 = allocMatrixEX( size / 2, size / 2 );
    int ** temp3 = allocMatrixEX( size / 2, size / 2 );
   
    /* OSTATECZNE SUMOWANIE WARTOŚCI WYJŚCIOWYCH */
    temp2 = add( temp2, M1, M4, size );
    temp3 = subtract( temp3, temp2, M5, size );
    C11 = add( C11, temp3, M7, size );
    C12 = add( C12, M3, M5, size );
    C21 = add( C21, M2, M4, size );
    temp2 = subtract( temp2, M1, M2, size );
    temp3 = add( temp3, temp2, M3, size );
    C22 = add( C22, temp3, M6, size );
    /* KONIEC SUMOWANIA */
   
    deleteMatrixEX( M1, size / 2 );
    deleteMatrixEX( M2, size / 2 );
    deleteMatrixEX( M3, size / 2 );
    deleteMatrixEX( M4, size / 2 );
    deleteMatrixEX( M5, size / 2 );
    deleteMatrixEX( M6, size / 2 );
    deleteMatrixEX( M7, size / 2 );
    deleteMatrixEX( temp2, size / 2 );
    deleteMatrixEX( temp3, size / 2 );
   
    /* SKŁADANIE MACIERZY WYJŚCIOWEJ */
    summary( MXC, C11, C12, C21, C22, size );
   
    deleteMatrixEX( C11, size / 2 );
    deleteMatrixEX( C12, size / 2 );
    deleteMatrixEX( C21, size / 2 );
    deleteMatrixEX( C22, size / 2 );
   
    cout << "zaraz zwroce macierz wyjsciowa M o wyniku:\n\n";
    //showM(MXC, size);
    copyMX( matrixC, MXC, size );
    //deleteMatrixEX(MXC, size);
    return matrixC;
}
P-154527
DejaVu
» 2016-12-06 16:03:55
https://pl.wikipedia.org/wiki​/Algorytm_Strassena

Tu jest implementacja w C.

Możesz więc:
1. porównać wyniki swoje z wynikami jakie zwraca ten program
2. poszukać możliwie małego scenariusza testowego na którym będziesz w stanie stwierdzić, że wyniki się różnią
3. uruchomić debugger i poszukać na tym małym przykładzie jakie obliczenia są wykonywane przez jeden i przez drugi program.
P-154586
darko202
» 2016-12-07 14:42:23
1.
Moim zdaniem (po analizie kodu)

popełniasz klasyczny błąd :
* rezerwujesz pamięć na zmienną
* tworzysz wskaźnik do tej zmiennej
* zwalniasz pamięć
* odwołujesz się poprzez wskaźnik do zwolnionego już obszaru pamięci


np.
w
int ** strassen(...){
...
deleteMatrixEX( matrixA, size );  // tu zwalniasz pamięć
...
// tu gdzieś odwołujesz się do zwolnionej wyżej pamięci
...
}

przepraszam, ale nie mam w tej chwili więcej czasu na dokładną analizę :(
nie jestem pewny (100%), ale tak to wygląda.


2.
>> "spróbowałem rekurencji"
to co przedstawiłeś nie jest rekurencją

np. funkcja rekurencyjna
int F(int n)
{
...
 if (n>1)
   return F(n-1) * 5*n;
 else
   return 3;
}
P-154616
« 1 »
  Strona 1 z 1