Оптимизированное матричное умножение в 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

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

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

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

Answer 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];

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

Answer 4

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

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

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

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

Answer 5

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

Блокировка

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

Настройка

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

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

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

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

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

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

Answer 6

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

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

Answer 7

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

Answer 8

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

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. Поскольку умножение обходится дороже, чем сложение, вы хотите, чтобы машина распараллеливала его как можно больше, поэтому сохранение ваших остановок для сложения означает, что вы тратите меньше времени на ожидание в цикле сложения, чем на цикл умножения.

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

Answer 9

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

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

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

Answer 10

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

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

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