Перемножение матриц в С++

[Andrey_Kireew]

В моё приложении на С++ основное время затрачивается на перемножение квадратных матриц,
которое реализуется следующим образом

1. ....................................... 
2. for (ii=0; ii<NStr; ii++) 
3.       for (j=0; j<NStr; j++) 
4.       {   ind=j*NStr; 
5.           *(D+ii+ind)=-*(A+ii)**(B+ind); 
6.           for (k=1; k<NStr; k++) 
7.               *(D+ii+ind)-=*(A+ii+k*NStr)**(B+k+ind); 
8.       } 
9. ...................................... 

Матрицы размерами 1000х1000 перемножаются примерно за 10 секунд (проверял с помощью функции clock()),
однако проведя аналогичную операцию в системе Matlab , я был мягко говоря озадачен

1. >> A=randn(1000,1000); 
2. >> B=randn(1000,1000); 
3. >> tic;D=A*B;toc; 
4. Elapsed time is 0.122922 seconds. 
5. >>  

Matlab выполняет перемножение примерно в 100 раз быстрее.
Не могу понять в чём дело? как можно оптимизировать алгоритм?
вообще программы на С++ обычно работают быстрее чем в Matlab,
кроме того, насколько я понял никаких быстрых алгоритмов перемножения матриц такого размера не существует.
Подскажите пожалуйста как лучше организовать вычисления?

[iifat]

А почему там -= вместо +=?
А ** — это умножить на число по адресу? Или ++ понимает как возведение в степень?
Можно попробовать ввести ещё несколько указателей для исключения умножения при взятии элементов матриц по индексам.
Возможно, функция Matlab написана на ассемблере, а это ещё быстрее.

[aa_dav]

Посмотрите в википедии - в статье "умножение матриц" целый раздел выделен про "быстрое умножение матриц".
P.S.
Скорее всего вы неправильно произвели оценку. Между 1000 в степени 3 и 1000 в степени 2,3727 значительная разница.

[mihaild]

Поменяйте порядок циклов так, чтобы самый внутренний шел по памяти подряд (если я правильно понял, это будет цикл по ii). Это потребует предварительно инициализировать D нулями, но не страшно.

Потом посмотрите, справляется ли компилятор с векторизацией и попробуйте скомпилировать с флагами -Ofast -march=native -mtune=native

Алгоритм Копперсмита-Винограда на практике интереса не представляет.

[Xaositect]

Andrey_Kireew в сообщении #1225268 писал(а):

В моё приложении на С++ основное время затрачивается на перемножение квадратных матриц,

Возьмите BLAS.

aa_dav в сообщении #1225284 писал(а):

Посмотрите в википедии - в статье "умножение матриц" целый раздел выделен про "быстрое умножение матриц".
P.S.
Скорее всего вы неправильно произвели оценку. Между 1000 в степени 3 и 1000 в степени 2,3727 значительная разница.

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

[iakovk]

Andrey_Kireew в сообщении #1225268 писал(а):

Матрицы размерами 1000х1000 перемножаются примерно за 10 секунд (проверял с помощью функции clock()),

А у меня получилось 8 сек. Без особых ухищрений.
Visual Studio 2015, Debug x86 (все по умолчанию); Тип данных double.
Если немножко пошаманить, то меньше одной сек - но до матлаба не дотягивает.
"Пошаманить" значит: включить Release конфигурацию, включить оптимизацию Full и favor for speed. А также использовать тип данных float.
Как ни странно, если все это сделать, то

mihaild в сообщении #1225288 писал(а):

Поменяйте порядок циклов так, чтобы самый внутренний шел по памяти подряд

приводит к замедлению.
Код:

Используется синтаксис C++

#include <stdlib.h>
#include <chrono>
#include <iostream>

void main()
{
        typedef float elemtype;
        static elemtype a[1000][1000];
        static elemtype b[1000][1000];
        static elemtype c[1000][1000];
        for (int j1 = 0; j1 < 1000; j1++)
        {
                for (int j2 = 0; j2 < 1000; j2++)
                {
                        a[j1][j2] = elemtype(rand());
                        b[j1][j2] = elemtype(rand());
                        c[j1][j2] = 0;
                }
        }
        std::cout << "Multiplication started" << std::endl;
        auto startTime = std::chrono::system_clock::now();
        for (int j1 = 0; j1 < 1000; j1++)
        {
                for (int j2 = 0; j2 < 1000; j2++)
                {
                        for (int j3 = 0; j3 < 1000; j3++)
                        {
                                c[j1][j2] += a[j1][j3] * b[j3][j2];
                        }
                }
        }
        std::chrono::duration<double> delay = std::chrono::system_clock::now() - startTime;
        std::cout << "Multiplication ended: delay = " << delay.count() << std::endl;
}
 

[Pphantom]

Xaositect в сообщении #1225303 писал(а):

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

Так $O(n^{\approx 2.37})$ - это не Штрассен, это Копперсмит-Виноград (который и при любых реальных размерах проиграет из-за константы).

[mihaild]

iakovk в сообщении #1225338 писал(а):

Visual Studio 2015, Debug x86

Смотреть на производительность в отладочной сборке - бесполезно, там получается случайное число.

iakovk в сообщении #1225338 писал(а):

приводит к замедлению.

Ну у вас же не подряд идет. Внутренний цикл у вас по j3, который идет старшим индексом у b. В вашем коде правильный порядок - j1, j3, j2. На моей машине ваш код отрабатывает за 6.5с, после перестановки циклов - за 0.22с.

[Andrey_Kireew]

Спасибо за многочисленные комментарии,
отвечаю на вопросы по порядку:

iifat в сообщении #1225282 писал(а):

А почему там -= вместо +=?

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

iifat в сообщении #1225282 писал(а):

А ** — это умножить на число по адресу?

совершенно верно

iifat в сообщении #1225282 писал(а):

Можно попробовать ввести ещё несколько указателей для исключения умножения при взятии элементов матриц по индексам.

я пробовал это делать - но производительность не меняется, переменная ind введена по тем же соображениям - однако если её убрать и в скобках написать

1.  *(D+ii+j*NStr)-=*(A+ii+k*NStr)**(B+k+j*NStr);  

производительность остаётся такой же, так что видимо это не самое узкое место

iifat в сообщении #1225282 писал(а):

Возможно, функция Matlab написана на ассемблере, а это ещё быстрее

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

aa_dav в сообщении #1225284 писал(а):

Посмотрите в википедии

это я всё смотрел, где то в сети написано, что в Matlab реализован эффективный алгоритм перемножения, но что за алгоритм конкретно - не сказано.

Цитата:

Скорее всего вы неправильно произвели оценку. Между 1000 в степени 3 и 1000 в степени 2,3727 значительная разница.

что я 0.1 секунду от 10 неправильно отличил? показатель степени как таковой меня интересует мало, хотя согласен с последующими комментариями 2,3727 - это вряд ли, так как это практически бесполезный алгоритм. Скорее всего там алгоритм Штрассена с логарифмом 2.8 и при этом маленькая константа, но это всё мои догадки ...

Xaositect в сообщении #1225303 писал(а):

Возьмите BLAS

спасибо, попробую

iakovk в сообщении #1225338 писал(а):

А у меня получилось 8 сек. Без особых ухищрений

Это только подтверждает, что никаких особых ошибок у меня нет.
Замедление на пару сек. может быть из за изменения загрузки компьютера, у меня каждый раз получаются разные цифры, так что на отклонения менее 10-15% думаю ориентироваться не стоит. Прирост должен быть в разы, а в идеале - на порядки.

mihaild в сообщении #1225352 писал(а):

На моей машине ваш код отрабатывает за 6.5с, после перестановки циклов - за 0.22с.

спасибо, попробую, это уже кое что

[dsge]

Вы можете скомпилировать м-файл и потом посмотреть С++ коды, какой получится алгоритм.

[mihaild]

Еще можно спросить компилятор, векторизует ли он код (это тоже влияет на производительность):

$ clang++ -std=c++14 -Ofast -march=native -mtune=native me.cpp -Rpass=loop-vectorize -Rpass-missed=loop-vectorize  -Rpass-analysis=loop-vectorize
me.cpp:24:4: remark: vectorized loop (vectorization width: 8, interleaved count: 4) [-Rpass=loop-vectorize]
                        for (int j2 = 0; j2 < 1000; j2++) {
                        ^

Если я правильно помню, на матрицах такого размера MKL обгонял наивное перемножение (естественно с линеаризованными матрицами) с оптимизациями компилятора примерно на 30% (и думаю это можно с хорошей точностью считать лучшим возможным результатом).

[mustitz]

А MATLAB в один поток считает?

[Dmitriy40]

По моему у Вас сильно недооптимизированы вычисления (по вине компилятора), за то же время можно выполнить в тысячу раз больше операций, если использовать современный процессор и FMA расширения системы команд:

SY1234 в сообщении #1043753 писал(а):

Haswell - очень удачная архитектура: позволяет, например, на 6-ти ядерном процессоре с частотой 3.3 GHz перемножить две квадратные матрицы double с n=10000 всего за 7.0 с для оси x64.

Матрицы вдесятеро больше, а время примерно то же! Я тогда прикидывал по тактам, выходило практически пиковая теоретическая производительность процессора по умножениям достигалась. На немного более старом процессоре (лишь с AVX расширением) скорость умножений упадёт ровно вдвое. Без даже AVX совсем грустно (ещё ровно вдвое медленнее).

(Прикидка теоретической скорости)

Можно и снова прикинуть, учитываем только умножения (сложения выполняются бесплатно, в параллельном потоке на том же ядре процессора), скорость работы с памятью считаем уже оптимизированной и не тормозящей процесс (именно для этого циклы переставляются). FMA команда выполняет 4 умножения за полтакта в потоке, т.е. 8 умножений за такт. AVX на более старом процессоре позволяет умножать 4 числа за такт.
Для перемножения матриц 1000х1000 надо выполнить суммарно $10^9$ умножений. Время перемножения будет не меньше $T=\dfrac{10^9}{8|4 \cdot N_\text{core} \cdot \text{Freq}}= 83|42 \text{мс}$ на каждый поток для 3ГГц процессора (83 для более старого процессора без FMA). На четырёх потоках время может быть 10 (FMA) - 20 (AVX) мс.
Повторю, с использованием оптимизированных библиотек (от той же Intel) или с оптимизацией на асме вполне достигается до 90% этой теоретической скорости.
Разумеется это лишь чистое перемножение матриц, без дополнительных операций над ними.

[Andrey_Kireew]

mustitz в сообщении #1225403 писал(а):

А MATLAB в один поток считает?

да в один поток,
кстати экспонента алгоритма перемножения matlab примерно 2.69, это лучше, чем алгоритм Штрассена

[Andrey_Kireew]

mihaild в сообщении #1225288 писал(а):

Поменяйте порядок циклов так, чтобы самый внутренний шел по памяти подряд (если я правильно понял, это будет цикл по ii). Это потребует предварительно инициализировать D нулями, но не страшно.

попробовал как Вы советовали, с небольшими усовершенствованиями

1. indD=D; 
2.    for (j=0; j<NStr; j++) 
3.    { indB=B; 
4.       for (ii=0; ii<NStr; ii++) 
5.       { 
6.          *indD=*indB**(M+j*NStr); 
7.            indD++; indB++; 
8.       } 
9.    } 
10.     
11.   indB=B; 
12.   for (k=1; k<NStr; k++) 
13.   {   indD=D; 
14.       for (j=0; j<NStr; j++) 
15.       { val=*(M+k+j*NStr);   
16.           for (ii=0; ii<NStr; ii++) 
17.           { 
18.               *indD+=*indB*val; 
19.                indD++; indB++; 
20.           } 
21.                indB-=NStr; 
22.       } 
23.       indB++; 
24.   } 

теперь цикл выполняется за 1.5 секунды (вместо 10), ускорение в 6 с лишним раз