Хотя известно, что использование вложенных std::vector для представления матриц - плохая идея , давайте пока воспользуемся им, поскольку он гибкий и многие существующие функции могут обрабатывать std::vector.
Я думал, что в небольших случаях разницу в скорости можно игнорировать. Но оказалось, что vector<vector<double>> в 10 + раз медленнее , чем numpy.dot().
Пусть A и B - матрицы с размером size x size , Предполагая, квадратные матрицы просто для простоты. (Мы не намерены ограничивать обсуждение случаем квадратных матриц.) Мы инициализируем каждую матрицу детерминистическим c способом и, наконец, вычисляем C = A * B.
Мы определяем «время вычисления» как время прошло только для расчета C = A * B. Другими словами, различные накладные расходы не включены.
Python3 код
import numpy as np
import time
import sys
if (len(sys.argv) != 2):
print("Pass `size` as an argument.", file = sys.stderr);
sys.exit(1);
size = int(sys.argv[1]);
A = np.ndarray((size, size));
B = np.ndarray((size, size));
for i in range(size):
for j in range(size):
A[i][j] = i * 3.14 + j
B[i][j] = i * 3.14 - j
start = time.time()
C = np.dot(A, B);
print("{:.3e}".format(time.time() - start), file = sys.stderr);
C ++ код
using namespace std;
#include <iostream>
#include <vector>
#include <chrono>
int main(int argc, char **argv) {
if (argc != 2) {
cerr << "Pass `size` as an argument.\n";
return 1;
}
const unsigned size = atoi(argv[1]);
vector<vector<double>> A(size, vector<double>(size));
vector<vector<double>> B(size, vector<double>(size));
for (int i = 0; i < size; ++i) {
for (int j = 0; j < size; ++j) {
A[i][j] = i * 3.14 + j;
B[i][j] = i * 3.14 - j;
}
}
auto start = chrono::system_clock::now();
vector<vector<double>> C(size, vector<double>(size, /* initial_value = */ 0));
for (int i = 0; i < size; ++i) {
for (int j = 0; j < size; ++j) {
for (int k = 0; k < size; ++k) {
C[i][j] += A[i][k] * B[k][j];
}
}
}
cerr << scientific;
cerr.precision(3);
cerr << chrono::duration<double>(chrono::system_clock::now() - start).count() << "\n";
}
код C ++ (многопоточный)
Мы также написали многопоточную версию кода C ++, поскольку numpy.dot() автоматически рассчитывается параллельно .
Вы можете получить все коды из GitHub .
Результат
C++ версия в 10+ раз медленнее, чем Python 3 (с numpy) версия.
matrix_size: 200x200
--------------- Time in seconds ---------------
C++ (not multithreaded): 8.45e-03
C++ (1 thread): 8.66e-03
C++ (2 threads): 4.68e-03
C++ (3 threads): 3.14e-03
C++ (4 threads): 2.43e-03
Python 3: 4.07e-04
-----------------------------------------------
matrix_size: 400x400
--------------- Time in seconds ---------------
C++ (not multithreaded): 7.011e-02
C++ (1 thread): 6.985e-02
C++ (2 threads): 3.647e-02
C++ (3 threads): 2.462e-02
C++ (4 threads): 1.915e-02
Python 3: 1.466e-03
-----------------------------------------------
Вопрос
Есть ли способ ускорить реализацию C ++?
Оптимизации, которые я пробовал
порядок вычисления свопа -> максимум в 3,5 раза быстрее (не код numpy, но код C ++)
оптимизация 1 плюс частичная развертка -> максимум в 4,5 раза быстрее, , но это может быть сделано только тогда, когда заранее известно size Нет. Как указано в этот комментарий , size не является обязательным закончили быть известным. Мы можем просто ограничить максимальное значение l oop переменных развернутых циклов и обработать оставшиеся элементы обычными циклами. См., Например, моя реализация .
оптимизация 2, плюс минимизация вызова C[i][j] путем введения простой переменной sum -> максимально в 5,2 раза быстрее , Реализация здесь . Этот результат подразумевает, что std::vector::operator[] является неузнаваемо медленным.
оптимизация 3, плюс флаг g ++ -march=native -> максимум в 6,2 раза быстрее (кстати, мы используем -O3 конечно.)
Оптимизация 3, плюс уменьшение вызова оператора [] путем введения указателя на элемент A, поскольку элементы A последовательно доступны в развернутый л oop. -> Максимум в 6,2 раза быстрее и немного быстрее, чем оптимизация 4. Код показан ниже.
g ++ -funroll-loops флаг для развертывания for циклов -> нет изменить
g ++ #pragma GCC unroll n -> без изменений
g ++ -flto флаг для включения оптимизации времени соединения -> без изменений
блочный алгоритм -> без изменений
транспонировать B, чтобы избежать кэширования пропустить -> без изменений
длинная линейная std::vector вместо вложенной std::vector<std::vector>, порядок вычисления свопов, алгоритм блокировки и частичная развертка -> максимум в 2,2 раза быстрее
Оптимизация 1, плюс P GO (оптимизация по профилю) -> в 4,7 раза быстрее
Оптимизация 3 , плюс P GO -> то же самое, что Оптимизация 3
Оптимизация 3, плюс спецификация g ++ c __builtin_prefetch() -> то же самое, что Оптимизация 3
Текущий статус
(о * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *) * * * * * * * * 11 * * 11 * * 11 * * 11 * *. Но давайте приведем некоторые коды, все из которых являются функциями, вызываемыми из многопоточной версии кода C ++.
Оригинальный код ( GitHub )
void f(const vector<vector<double>> &A, const vector<vector<double>> &B, vector<vector<double>> &C, unsigned row_start, unsigned row_end) {
const unsigned j_max = B[0].size();
const unsigned k_max = B.size();
for (int i = row_start; i < row_end; ++i) {
for (int j = 0; j < j_max; ++j) {
for (int k = 0; k < k_max; ++k) {
C[i][j] += A[i][k] * B[k][j];
}
}
}
}
Текущий лучший код ( GitHub )
Это реализация описанной выше оптимизации 5.
void f(const vector<vector<double>> &A, const vector<vector<double>> &B, vector<vector<double>> &C, unsigned row_start, unsigned row_end) {
static const unsigned num_unroll = 5;
const unsigned j_max = B[0].size();
const unsigned k_max_for_unrolled_loop = B.size() / num_unroll * num_unroll;
const unsigned k_max = B.size();
for (int i = row_start; i < row_end; ++i) {
for (int k = 0; k < k_max_for_unrolled_loop; k += num_unroll) {
for (int j = 0; j < j_max; ++j) {
const double *p = A[i].data() + k;
double sum;
sum = *p++ * B[k][j];
sum += *p++ * B[k+1][j];
sum += *p++ * B[k+2][j];
sum += *p++ * B[k+3][j];
sum += *p++ * B[k+4][j];
C[i][j] += sum;
}
}
for (int k = k_max_for_unrolled_loop; k < k_max; ++k) {
const double a = A[i][k];
for (int j = 0; j < j_max; ++j) {
C[i][j] += a * B[k][j];
}
}
}
}
Мы перепробовали много оптимизаций с тех пор, как впервые опубликовали этот вопрос. Мы потратили целых два дня на борьбу с этой проблемой и, наконец, пришли к тому, что у нас больше нет идей, как оптимизировать текущий лучший код. Мы сомневаемся, что более сложные алгоритмы, такие как Strassen's , сделают это лучше, поскольку случаи, которые мы обрабатываем, не велики, и каждая операция на std::vector настолько дорогая, что, как мы видели, просто сокращает вызов [] хорошо улучшили производительность.
Мы (хотим) верить, что мы можем улучшить ее, хотя.