Это должно во многом зависеть от точного разреженности матрицы и используемой платформы. Я протестировал несколько вещей с gcc 8.3.0 и флагами компилятора -O3 -march=native (что составляет -march=skylake на моем ЦП) в нижнем треугольнике этой матрицы измерения 3006 с 19554 ненулевыми записями. Надеюсь, это несколько близко к вашей настройке, но в любом случае я надеюсь, что они могут дать вам представление о том, с чего начать.
Для определения времени я использовал google / benchmark с this исходный файл . Он определяет benchBacksolveBaseline, который измеряет реализацию, представленную в вопросе, и benchBacksolveOptimized, который сравнивает предложенные "оптимизированные" реализации. Существует также benchFillRhs, который отдельно сравнивает функцию, которая используется в обоих для генерации некоторых не совсем тривиальных значений для правой части. Чтобы получить время «чистых» обратных процессов, нужно вычесть время, которое занимает benchFillRhs.
1. Итерация строго назад
Внешний l oop в вашей реализации перебирает столбцы в обратном направлении, в то время как внутренний l oop перебирает текущий столбец вперед. Похоже, было бы более последовательным итерировать каждый столбец в обратном направлении:
for (int i=n-1; i>=0; --i) {
for (int j=Lp[i+1]-1; j>=Lp[i]; --j) {
x[i] -= Lx[j] * x[Li[j]];
}
}
Это едва меняет сборку (https://godbolt.org/z/CBZAT5), но временные параметры тестирования показывают заметное улучшение :
------------------------------------------------------------------
Benchmark Time CPU Iterations
------------------------------------------------------------------
benchFillRhs 2737 ns 2734 ns 5120000
benchBacksolveBaseline 17412 ns 17421 ns 829630
benchBacksolveOptimized 16046 ns 16040 ns 853333
Я предполагаю, что это вызвано каким-то более предсказуемым доступом к кешу, но я не стал вдаваться в подробности.
2. Меньше нагрузок / накоплений во внутренней l oop
Поскольку A меньше tri angular, у нас есть i < Li[j]. Поэтому мы знаем, что x[Li[j]] не изменится из-за изменений в x[i] во внутреннем l oop. Мы можем применить эти знания в нашей реализации, используя временную переменную:
for (int i=n-1; i>=0; --i) {
double xi_temp = x[i];
for (int j=Lp[i+1]-1; j>=Lp[i]; --j) {
xi_temp -= Lx[j] * x[Li[j]];
}
x[i] = xi_temp;
}
Это заставит gcc 8.3.0 переместить хранилище в память из внутренней l oop прямо после его окончания ( https://godbolt.org/z/vM4gPD). Тест для тестовой матрицы в моей системе показывает небольшое улучшение:
------------------------------------------------------------------
Benchmark Time CPU Iterations
------------------------------------------------------------------
benchFillRhs 2737 ns 2740 ns 5120000
benchBacksolveBaseline 17410 ns 17418 ns 814545
benchBacksolveOptimized 15155 ns 15147 ns 887129
3. Разверните l oop
Пока clang уже начинает раскручивать l oop после первого предлагаемого изменения кода, gcc 8.3.0 по-прежнему нет. Итак, давайте попробуем, дополнительно передав -funroll-loops.
------------------------------------------------------------------
Benchmark Time CPU Iterations
------------------------------------------------------------------
benchFillRhs 2733 ns 2734 ns 5120000
benchBacksolveBaseline 15079 ns 15081 ns 953191
benchBacksolveOptimized 14392 ns 14385 ns 963441
Обратите внимание, что базовая линия также улучшается, так как l oop в этой реализации также развернут. Наша оптимизированная версия также выигрывает от развертывания l oop, но, возможно, не так, как нам хотелось бы. Глядя на сгенерированную сборку (https://godbolt.org/z/_LJC5f), кажется, что gcc, возможно, ушло немного далеко с 8 развертываниями. Для моей настройки, я могу сделать немного лучше, просто развернув вручную. Поэтому снова сбросьте флаг -funroll-loops и выполните развертывание примерно так:
for (int i=n-1; i>=0; --i) {
const int col_begin = Lp[i];
const int col_end = Lp[i+1];
const bool is_col_nnz_odd = (col_end - col_begin) & 1;
double xi_temp = x[i];
int j = col_end - 1;
if (is_col_nnz_odd) {
xi_temp -= Lx[j] * x[Li[j]];
--j;
}
for (; j >= col_begin; j -= 2) {
xi_temp -= Lx[j - 0] * x[Li[j - 0]] +
Lx[j - 1] * x[Li[j - 1]];
}
x[i] = xi_temp;
}
С этим я измеряю:
------------------------------------------------------------------
Benchmark Time CPU Iterations
------------------------------------------------------------------
benchFillRhs 2728 ns 2729 ns 5090909
benchBacksolveBaseline 17451 ns 17449 ns 822018
benchBacksolveOptimized 13440 ns 13443 ns 1018182
Другие алгоритмы
Все эти версии по-прежнему используют одну и ту же простую реализацию обратного решения для структуры разреженной матрицы. По сути, работа с такими разреженными матричными структурами может иметь значительные проблемы с трафиком памяти c. По крайней мере, для матричных факторизаций существуют более сложные методы, которые работают на плотных подматрицах, собранных из разреженной структуры. Примерами являются супернодальные и мультифронтальные методы. Я немного размышляю над этим, но я думаю, что такие методы также применят эту идею для компоновки и используют плотные матричные операции для более низких три angular обратных решений (например, для факторизаций типа Холецкого). Поэтому, возможно, стоит взглянуть на такие методы, если вы не обязаны придерживаться простого метода, который напрямую работает с разреженной структурой. См. Например этот опрос Дэвиса.