Алгоритм Штрассена. Оптимизация

Привет, Хабр.
Имеется алгоритм Штрассена умножения матриц, но скоростью он не блещет. Буду рад увидеть ваши замечания по оптимизации кода.

Для матриц размером 512*512 дал такие результаты:
Normal mode 1538911
Strassen mode 784078

Надеюсь код читаемый :)
int Strassen(int N, int **MatrixA, int **MatrixB, int **MatrixC)
{
    
    int HalfSize = N/2;
    int newSize = N/2;
    
    if ( N <= 32 )//choosing the threshold is extremely important, try N<=2 to see the result
    {
        MUL(MatrixA,MatrixB,MatrixC,N);
    }
    else
    {
        int** A11;
        int** A12;
        int** A21;
        int** A22;
        
        int** B11;
        int** B12;
        int** B21;
        int** B22;
        
        int** C11;
        int** C12;
        int** C21;
        int** C22;
        
        int** M1;
        int** M2;
        int** M3;
        int** M4;
        int** M5;
        int** M6;
        int** M7;
        int** AResult;
        int** BResult;
        
        //making a 1 diminsional pointer based array.
        A11 = new int *[newSize];
        A12 = new int *[newSize];
        A21 = new int *[newSize];
        A22 = new int *[newSize];
        
        B11 = new int *[newSize];
        B12 = new int *[newSize];
        B21 = new int *[newSize];
        B22 = new int *[newSize];
        
        C11 = new int *[newSize];
        C12 = new int *[newSize];
        C21 = new int *[newSize];
        C22 = new int *[newSize];
        
        M1 = new int *[newSize];
        M2 = new int *[newSize];
        M3 = new int *[newSize];
        M4 = new int *[newSize];
        M5 = new int *[newSize];
        M6 = new int *[newSize];
        M7 = new int *[newSize];
        
        AResult = new int *[newSize];
        BResult = new int *[newSize];
        
        int newLength = newSize;
        
        //making that 1 diminsional pointer based array , a 2D pointer based array
        for ( int i = 0; i < newSize; i++)
        {
            A11[i] = new int[newLength];
            A12[i] = new int[newLength];
            A21[i] = new int[newLength];
            A22[i] = new int[newLength];
            
            B11[i] = new int[newLength];
            B12[i] = new int[newLength];
            B21[i] = new int[newLength];
            B22[i] = new int[newLength];
            
            C11[i] = new int[newLength];
            C12[i] = new int[newLength];
            C21[i] = new int[newLength];
            C22[i] = new int[newLength];
            
            M1[i] = new int[newLength];
            M2[i] = new int[newLength];
            M3[i] = new int[newLength];
            M4[i] = new int[newLength];
            M5[i] = new int[newLength];
            M6[i] = new int[newLength];
            M7[i] = new int[newLength];
            
            AResult[i] = new int[newLength];
            BResult[i] = new int[newLength];
            
            
        }
        //splitting input Matrixes, into 4 submatrices each.
        for (int i = 0; i < N / 2; i++)
        {
            for (int j = 0; j < N / 2; j++)
            {
                A11[i][j] = MatrixA[i][j];
                A12[i][j] = MatrixA[i][j + N / 2];
                A21[i][j] = MatrixA[i + N / 2][j];
                A22[i][j] = MatrixA[i + N / 2][j + N / 2];
                
                B11[i][j] = MatrixB[i][j];
                B12[i][j] = MatrixB[i][j + N / 2];
                B21[i][j] = MatrixB[i + N / 2][j];
                B22[i][j] = MatrixB[i + N / 2][j + N / 2];
                
            }
        }
        
        //here we calculate M1..M7 matrices .
        //M1[][]
        ADD( A11,A22,AResult, HalfSize);
        ADD( B11,B22,BResult, HalfSize);
        Strassen( HalfSize, AResult, BResult, M1 ); //now that we need to multiply this , we use the strassen itself .
        
        
        //M2[][]
        ADD( A21,A22,AResult, HalfSize);              //M2=(A21+A22)B11
        Strassen(HalfSize, AResult, B11, M2);       //Mul(AResult,B11,M2);
        
        //M3[][]
        SUB( B12,B22,BResult, HalfSize);              //M3=A11(B12-B22)
        Strassen(HalfSize, A11, BResult, M3);       //Mul(A11,BResult,M3);
        
        //M4[][]
        SUB( B21, B11, BResult, HalfSize);           //M4=A22(B21-B11)
        Strassen(HalfSize, A22, BResult, M4);       //Mul(A22,BResult,M4);
        
        //M5[][]
        ADD( A11, A12, AResult, HalfSize);           //M5=(A11+A12)B22
        Strassen(HalfSize, AResult, B22, M5);       //Mul(AResult,B22,M5);
        
        
        //M6[][]
        SUB( A21, A11, AResult, HalfSize);
        ADD( B11, B12, BResult, HalfSize);             //M6=(A21-A11)(B11+B12)
        Strassen( HalfSize, AResult, BResult, M6);    //Mul(AResult,BResult,M6);
        
        //M7[][]
        SUB(A12, A22, AResult, HalfSize);
        ADD(B21, B22, BResult, HalfSize);             //M7=(A12-A22)(B21+B22)
        Strassen(HalfSize, AResult, BResult, M7);     //Mul(AResult,BResult,M7);
        
        //C11 = M1 + M4 - M5 + M7;
        ADD( M1, M4, AResult, HalfSize);
        SUB( M7, M5, BResult, HalfSize);
        ADD( AResult, BResult, C11, HalfSize);
        
        //C12 = M3 + M5;
        ADD( M3, M5, C12, HalfSize);
        
        //C21 = M2 + M4;
        ADD( M2, M4, C21, HalfSize);
        
        //C22 = M1 + M3 - M2 + M6;
        ADD( M1, M3, AResult, HalfSize);
        SUB( M6, M2, BResult, HalfSize);
        ADD( AResult, BResult, C22, HalfSize);
        
        
        //at this point , we have calculated the c11..c22 matrices, and now we are going to
        //put them together and make a unit matrix which would describe our resulting Matrix.
        for (int i = 0; i < N/2 ; i++)
        {
            for (int j = 0 ; j < N/2 ; j++)
            {
                MatrixC[i][j] = C11[i][j];
                MatrixC[i][j + N / 2] = C12[i][j];
                MatrixC[i + N / 2][j] = C21[i][j];
                MatrixC[i + N / 2][j + N / 2] = C22[i][j];
            }
        }
        
        // dont forget to free the space we alocated for matrices,
        for (int i = 0; i < newLength; i++)
        {
            delete[] A11[i];delete[] A12[i];delete[] A21[i];
            delete[] A22[i];
            
            delete[] B11[i];delete[] B12[i];delete[] B21[i];
            delete[] B22[i];
            delete[] C11[i];delete[] C12[i];delete[] C21[i];
            delete[] C22[i];
            delete[] M1[i];delete[] M2[i];delete[] M3[i];delete[] M4[i];
            delete[] M5[i];delete[] M6[i];delete[] M7[i];
            delete[] AResult[i];delete[] BResult[i] ;
        }
        delete[] A11;delete[] A12;delete[] A21;delete[] A22;
        delete[] B11;delete[] B12;delete[] B21;delete[] B22;
        delete[] C11;delete[] C12;delete[] C21;delete[] C22;
        delete[] M1;delete[] M2;delete[] M3;delete[] M4;delete[] M5;
        delete[] M6;delete[] M7;
        delete[] AResult;
        delete[] BResult ;
        
        
    }//end of else
    
    
	return 0;
}

int ADD(int** MatrixA, int** MatrixB, int** MatrixResult, int MatrixSize )
{
    for ( int i = 0; i < MatrixSize; i++)
    {
        for ( int j = 0; j < MatrixSize; j++)
        {
            MatrixResult[i][j] =  MatrixA[i][j] + MatrixB[i][j];
        }
    }
	return 0;
}

int SUB(int** MatrixA, int** MatrixB, int** MatrixResult, int MatrixSize )
{
    for ( int i = 0; i < MatrixSize; i++)
    {
        for ( int j = 0; j < MatrixSize; j++)
        {
            MatrixResult[i][j] =  MatrixA[i][j] - MatrixB[i][j];
        }
    }
	return 0;
}

int MUL( int** MatrixA, int** MatrixB, int** MatrixResult, int MatrixSize )
{
    for (int i=0;i<MatrixSize ;i++)
    {
        for (int j=0;j<MatrixSize ;j++)
        {
            MatrixResult[i][j]=0;
            for (int k=0;k<MatrixSize ;k++)
            {
                MatrixResult[i][j]=MatrixResult[i][j]+MatrixA[i][k]*MatrixB[k][j];
            }
        }
    }
	return 0;
}



Спасибо!
  • Вопрос задан
  • 10814 просмотров
Решения вопроса 1
@gribozavr
Для начала преобразуйте все матрицы в плоские массивы int[N * N], а обращения как M[j * N + i].
Ответ написан
Пригласить эксперта
Ответы на вопрос 7
@BrainHacker
Блин, у вас столько динамического выделения памяти, что аж старшно становиться! Попробуйте хотя бы выделять память большой пачкой. А лучше вообще на старте работы алгоритма выделить всю необходимую память. Один раз.
Ответ написан
kefiijrw
@kefiijrw
У вас, по-моему, ошибка в функции MUL. Вспоминать алгоритм лень, поэтому предположение.

В строке
     MatrixResult[i][j]=MatrixResult[i][j]+MatrixA[i][k]*MatrixB[k][j];


по смыслу вместо = должно быть += или что-то вроде того.

Ну и может распараллелить?
Ответ написан
@bayandin
Можно почитать zealint.ru/fast-matrix-multiplication-results.html.
Позволю себе небольшую цитату из этой записи:

«Во-первых, сначала я написал самую тупую версию, которая умножает строку на столбец в тройном цикле (C [ i ] [ j ] += A [ i ] [ k ] * B [ k ] [ j ], циклы идут в порядке i-j-k). Это и были те самые 800 с. (надо же с чего-то начинать). Здесь же можно поиграть перестановкой циклов местами: самое быстрое получается, если выбрать i-k-j, в этом случае все три матрицы сканируются по строкам. Такая версия уже укладывалась в 3 минуты (что почти в 5 раз быстрее). Можно просто транспонировать вторую матрицу прямо на вводе и поменять местами индексы k и j при обращении к матрице B. Это ещё быстрее (150 с.). Что ещё можно сделать без вмешательства ассемблера? — поменять int (32 бита) на short (16 бит) для входящих матриц и скомпилировать родным для процессора компилятором с опциями максимального ускорения. Это дает 94 секунды. До 55 секунд можно дойти, если развернуть матрицы в массив. Тогда будут получаться записи вида A [ i * n + k ], но все такие умножения нужно либо вынести за внутренний цикл, либо вообще заменить сложением указателей с числом n в нужных местах. В этом случае компилятор оптимизирует гораздо лучше. Все, без ассемблера можно сделать только ещё одну вещь: разбить матрицу на блоки, которые помещаются в кэш, но я решил, что сделаю это сразу на ассемблере.»
Ответ написан
@Eddy_Em
Вам нужно именно реализовать свой метод (т.е. это — какая-то лабораторная) или же вы изобретаете велосипед?
Если второе, то советую воспользоваться, к примеру, библиотеками BLAS или GSL.
Ответ написан
@IlVin
Прочитайте оптимизации Алексея Тутубалина — он ещё в лохматом 2007 году перемножением матриц процессоры «щупал»: blog.lexa.ru/2007/01/04/o_peremnozhenii_matric_i_prochix_arxitekturnix_zamorochkax.html
Ответ написан
Nomad1
@Nomad1
Дам Вам необычный совет: напишите сначала код на Managed языке программирования (C#, Java, или вообще что-либо высокого уровня, без выделения памяти, указателей и пр.). Затем примените простейшие оптимизации (одномерные массивы, вынос всяких N/2 в отдельные переменные). Полученную программу назовите прототипом, пересоберите в релизном режиме, замерьте скорость и заморозьте код.
Затем, любая реализация на другом языке должна сверяться с прототипом по результату и скорости работы. К примеру, приведенный в вопросе C++ код будет выполняться медленнее C# кода за счет более правильной работы с памятью и хорошей оптимизируемости Jit-ом у второго. Любое ответвление, где новый код стал медленее прототипа надо выкидывать, например с не-оптимальным использованием SIMD инструкций, динамическим выделением памяти в цикле, неудачным распараллеливанием.
Такой подход буквально к любой time-critical задаче (прототипирование на языке с высоким уровнем абстракции, а затем реализация) позволит не уходить в дебри и не оптимизировать вручную раз то, что отлично могут сделать компилятор на пару с оптимизатором.
Ответ написан
Ваш ответ на вопрос

Войдите, чтобы написать ответ

Войти через центр авторизации
Похожие вопросы