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

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

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

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

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

Ответ 1

Что каждый программист должен знать о памяти (pdf-ссылка) Ульриха Дреппера имеет много хороших идей об эффективности памяти, но, в частности, он использует матричное умножение в качестве примера того, как знание о памяти и использование этого знания могут ускорить этот процесс. Посмотрите приложение A.1 в своей статье и прочитайте раздел 6.2.1. В таблице 6.2 в документе показано, что он может получить свое время работы на 10% от наивного времени реализации для матрицы 1000x1000.

Конечно, его окончательный код довольно волосатый и использует множество системных вещей и настройку времени компиляции, но все же, если вам действительно нужна скорость, чтение этой бумаги и чтение его реализации определенно стоит того.

Ответ 2

Получение этого права может быть нетривиальным. Одна из оптимизаций, которая имеет особое значение для больших матриц, - это мозаика умножения для сохранения содержимого в кэше. Однажды я измерил разницу в производительности в 12 раз, но я специально выбрал размер матрицы, который потреблял кратные значения моего кеша (около 97 года, поэтому кеш был маленьким).

Там много литературы на эту тему. Отправной точкой является:

http://en.wikipedia.org/wiki/Loop_tiling

Для более глубокого изучения могут быть полезны следующие ссылки, особенно книги Банерджи:

[Ban93] Banerjee, Utpal, Loop Transformations для реструктурируемых компиляторов: фонды, Kluwer Academic Publishers, Norwell, MA, 1993.

[Ban94] Banerjee, Utpal, Loop Parallelization, Kluwer Academic Publishers, Norwell, MA, 1994.

[BGS93] Бэкон, Дэвид Ф., Сьюзен Л. Грэм и Оливер Шарп, Преобразования компиляторов для высокопроизводительных вычислений, Отдел компьютерных наук, Университет Калифорнии, Беркли, Калифорния, Технический отчет № UCB/CSD-93-781.

[LRW91] Лэм, Моника С., Эдвард Э. Ротберг и Майкл Э. Вольф. Производительность кэша и оптимизация заблокированных алгоритмов, в 4-й Международной конференции по архитектурной поддержке языков программирования, состоявшейся в Санта-Кларе, Калифорния, апрель 1991 г., 63-74.

[LW91] Лэм, Моника С. и Майкл Э. Вольф. Теория циклического преобразования и алгоритм максимизации параллелизма, в транзакциях IEEE на параллельных и распределенных системах, 1991, 2 (4): 452-471.

[PW86] Падуя, Дэвид А. и Майкл Дж. Вулф, Усовершенствованная оптимизация компиляторов для суперкомпьютеров, In Communications of ACM, 29 (12): 1184-1201, 1986.

[Wolfe89] Wolfe, Michael J. Оптимизация суперкомпиляторов для суперкомпьютеров, The MIT Press, Cambridge, MA, 1989.

[Wolfe96] Wolfe, Michael J., Высокопроизводительные компиляторы для параллельных вычислений, Addison-Wesley, CA, 1996.

Ответ 3

ВНИМАНИЕ: У вас есть ошибка в вашей второй реализации

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

Когда вы делаете f = 0, c = 1

        MatrixB[0][1] = MatrixB[1][0];

вы перезаписываете MatrixB[0][1] и теряете это значение! Когда петля попадает в f = 1, c = 0

        MatrixB[1][0] = MatrixB[0][1];

значение, скопированное, совпадает с тем, которое уже было там.

Ответ 4

Если матрица недостаточно велика или вы не повторяете операции, большое количество раз вы не увидите заметных различий.

Если матрица, скажем, 1,000x1000, вы увидите улучшения, но я бы сказал, что если она ниже 100x100, вы не должны беспокоиться об этом.

Кроме того, любое "улучшение" может быть порядка миллисекунд, если только yoy не работают с чрезвычайно большими матрицами или не повторяют операцию тысячи раз.

Наконец, если вы измените компьютер, который вы используете для более быстрого, различия будут еще более узкими!

Ответ 5

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

Не забывайте при измерении времени выполнения включить время для переноса матрицы B.

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

Ответ 6

Как большие улучшения, которые вы получите, будут зависеть от:

Размер кеша
Размер строки кэша
Степень ассоциативности кеша

Для небольших матричных размеров и современных процессоров весьма вероятно, что данные, которые будут отображаться как в MatrixA, так и в MatrixB, будут почти полностью сохранены в кеше после первого касания.

Ответ 7

Просто для вас попробовать (но это будет иметь значение только для больших матриц): отделите логику добавления от логики умножения во внутреннем цикле следующим образом:

for (k = 0; k < i; k++)
{
    int sums[i];//I know this size declaration is illegal in C. consider 
            //this pseudo-code.
    for (l = 0; l < i; l++)
        sums[l] = MatrixA[j][l]*MatrixB[k][l];

    int suma = 0;
    for(int s = 0; s < i; s++)
       suma += sums[s];
}

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

Это только моя логика. Другие, имеющие больше знаний в этой области, могут не согласиться.

Ответ 8

Сложность вычисления умножения двух матриц N * N равна O (N ^ 3). Эффективность будет значительно улучшена, если вы используете алгоритм O (N ^ 2.73), который, вероятно, был принят MATLAB. Если вы установили MATLAB, попробуйте умножить две матрицы 1024 * 1024. На моем компьютере MATLAB завершает его в 0,7, но реализация наивного алгоритма C\С++, подобного вашей, занимает 20 секунд. Если вы действительно заботитесь о производительности, обратитесь к более сложным алгоритмам. Я слышал, что существует O (N ^ 2.4) алгоритм, однако ему нужна очень большая матрица, так что другими манипуляциями можно пренебречь.

Ответ 9

Вы не должны писать умножение матрицы. Вы должны зависеть от внешних библиотек. В частности, вы должны использовать подпрограмму GEMM из библиотеки BLAS. GEMM часто обеспечивает следующие оптимизации

Блокировка

Эффективное умножение матрицы зависит от блокировки вашей матрицы и выполнения нескольких меньших блокированных умножений. В идеале размер каждого блока выбирается так, чтобы он хорошо вписывался в кеш, что значительно улучшало производительность.

Настройка

Идеальный размер блока зависит от базовой иерархии памяти (насколько большой кеш?). В результате библиотеки должны быть настроены и скомпилированы для каждой конкретной машины. Это делается, в частности, реализацией ATLAS BLAS.

Оптимизация уровня сборки

Матричная мультипликация настолько распространена, что разработчики будут ее оптимизировать вручную. В частности это делается в GotoBLAS.

Гетерогенные вычисления (GPU)

Matrix Multiply очень интенсивно работает на FLOP/compute, что делает его идеальным кандидатом для работы на графических процессорах. cuBLAS и MAGMA являются хорошими кандидатами для этого.

Короче говоря, плотная линейная алгебра является хорошо изученной темой. Люди посвящают свои жизни усовершенствованию этих алгоритмов. Вы должны использовать их работу; это сделает их счастливыми.

Ответ 10

не столь особенный, но лучше:

    c = 0;
do
{
    for (j = 0; j < i; j++)
    {
        for (k = 0; k < i; k++)
        {
            sum = 0; sum_ = 0;
            for (l = 0; l < i; l++) {
                MatrixB[j][k] = MatrixB[k][j];
                sum += MatrixA[j][l]*MatrixB[k][l];
                l++;
                MatrixB[j][k] = MatrixB[k][j];
                sum_ += MatrixA[j][l]*MatrixB[k][l];

                sum += sum_;
            }
            MatrixR[j][k] = sum;
        }
     }
     c++;
} while (c<iteraciones);

Ответ 11

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

Ответ 12

Очень старый вопрос, но это моя текущая реализация для моих проектов 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 заменяется размером матрицы. Поэтому, если вы умножаете матрицы 4x4, вы используете:

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

Ответ 13

Вообще говоря, перенос 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];

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

Ответ 14

Вас может заинтересовать этот ответ: fooobar.com/info/2305837/...