Оптимизированное матричное умножение в C

Question

Я пытаюсь сравнить разные методы умножения матриц. Первый нормальный метод:

do
{
    for (j = 0; j < i; j++)
    {
        for (k = 0; k < i; k++)
        {
            suma = 0;
            for (l = 0; l < i; l++)
                suma += MatrixA[j][l]*MatrixB[l][k];
                MatrixR[j][k] = suma;
            }
        }
    }
    c++;
} while (c<iteraciones);

Второй состоит из транспонирования матрицы B, а затем умножения на строки:

int f, co;
for (f = 0; f < i; f++) {
    for ( co = 0; co < i; co++) {
        MatrixB[f][co] = MatrixB[co][f];
    }
}

c = 0;
do
{
    for (j = 0; j < i; j++)
    {
        for (k = 0; k < i; k++)
        {
            suma = 0;
            for (l = 0; l < i; l++)
                suma += MatrixA[j][l]*MatrixB[k][l];
                MatrixR[j][k] = suma;
            }
        }
     }
     c++;
} while (c<iteraciones);

Второй метод должен быть намного быстрее, потому что мы получаем доступ к смежным слотам памяти, но я не получаю существенного улучшения производительности. Я что-то не так делаю?

Я могу опубликовать полный код, но я думаю, что это не нужно.

Answer 1

Вообще говоря, транспонирование B должно оказаться гораздо быстрее, чем наивная реализация, но за счет потери другой памяти NxN. Я только что провел неделю, разбираясь с оптимизацией умножения матриц, и пока абсолютный победитель в рукопашном бою таков:

for (int i = 0; i < N; i++)
    for (int k = 0; k < N; k++)
        for (int j = 0; j < N; j++)
            if (likely(k)) /* #define likely(x) __builtin_expect(!!(x), 1) */
                C[i][j] += A[i][k] * B[k][j];
            else
                C[i][j] = A[i][k] * B[k][j];

Это даже лучше, чем метод Дреппера, упомянутый в предыдущем комментарии, поскольку он работает оптимально, независимо от свойств кэша базового ЦП. Хитрость заключается в переупорядочении циклов так, чтобы все три матрицы были доступны в главном порядке строк.

Answer 2

Очень старый вопрос, но здесь моя текущая реализация для моих проектов opengl:

typedef float matN[N][N];

inline void matN_mul(matN dest, matN src1, matN src2)
{
    unsigned int i;
    for(i = 0; i < N^2; i++)
    {
        unsigned int row = (int) i / 4, col = i % 4;
        dest[row][col] = src1[row][0] * src2[0][col] +
                         src1[row][1] * src2[1][col] +
                         ....
                         src[row][N-1] * src3[N-1][col];
    }
}

Где N заменяется размером матрицы. Так что если вы умножаете матрицы 4х4, то вы используете:

typedef float mat4[4][4];    

inline void mat4_mul(mat4 dest, mat4 src1, mat4 src2)
{
    unsigned int i;
    for(i = 0; i < 16; i++)
    {
        unsigned int row = (int) i / 4, col = i % 4;
        dest[row][col] = src1[row][0] * src2[0][col] +
                         src1[row][1] * src2[1][col] +
                         src1[row][2] * src2[2][col] +
                         src1[row][3] * src2[3][col];
    }
}

Эта функция в основном сводит к минимуму циклы, но модуль может облагаться налогом ... На моем компьютере эта функция выполнялась примерно на 50% быстрее, чем функция умножения цикла на тройную.

Минусы:

• Требуется много кода (например, различные функции для mat3 x mat3, mat5 x mat5 ...)

• Твики, необходимые для нерегулярного умножения (например, mat3 x mat4) .....

Answer 3

Если вы работаете с небольшими числами, то упомянутое вами улучшение незначительно. Кроме того, производительность будет зависеть от оборудования, на котором вы работаете. Но если вы работаете с числами в миллионах, то это повлияет. Приходя в программу, можете ли вы вставить написанную вами программу.