Ответ 1
Я согласен с транспозицией для лучшего кэширования (но см. мои комментарии к этому в конце), и там больше делать, поэтому давайте посмотрим, что мы можем сделать с полной функцией...
Оригинальная функция, для справки (с некоторыми правилами для моего удобства):
void MultiDiagonalSymmetricMatrix::CholeskyBackSolve(float *x, float *b){
//We want to solve L D Lt x = b where D is a diagonal matrix described by Diagonals[0] and L is a unit lower triagular matrix described by the rest of the diagonals.
//Let D Lt x = y. Then, first solve L y = b.
float *y = new float[n];
float **d = IncompleteCholeskyFactorization->Diagonals;
unsigned int *s = IncompleteCholeskyFactorization->StartRows;
unsigned int M = IncompleteCholeskyFactorization->m;
unsigned int N = IncompleteCholeskyFactorization->n;
unsigned int i, j;
for(j = 0; j != N; j++){
float sub = 0;
for(i = 1; i != M; i++){
int c = (int)j - (int)s[i];
if(c < 0) break;
if(c==j) {
sub += d[i][c]*b[c];
} else {
sub += d[i][c]*y[c];
}
}
y[j] = b[j] - sub;
}
//Now, solve x from D Lt x = y -> Lt x = D^-1 y
// Took this one out of the while, so it can be parallelized now, which speeds up, because division is expensive
#pragma omp parallel for
for(j = 0; j < N; j++){
x[j] = y[j]/d[0][j];
}
while(j-- != 0){
float sub = 0;
for(i = 1; i != M; i++){
if(j + s[i] >= N) break;
sub += d[i][j]*x[j + s[i]];
}
x[j] -= sub;
}
delete[] y;
}
Из-за комментария о параллельном делении, дающем ускорение скорости (несмотря на то, что это только O (N)), я предполагаю, что сама функция получает много имен. Так зачем выделять память? Просто отметьте x как __restrict__ и измените y на x всюду (__restrict__ является расширением GCC, взятым из C99. Возможно, вы захотите использовать для него define. Возможно, в библиотеке уже есть один).
Точно так же, хотя, я думаю, вы не можете изменить подпись, вы можете заставить функцию принять только один параметр и изменить его. b никогда не используется, если x или y установлены. Это также означает, что вы можете избавиться от ветки в первом цикле, который запускает ~ N * M раз. Используйте memcpy в начале, если у вас должно быть 2 параметра.
И почему d массив указателей? Должно быть? Это кажется слишком глубоким в исходном коде, поэтому я не буду его трогать, но если есть возможность сгладить сохраненный массив, это будет ускорение скорости, даже если вы не можете его транспонировать (умножить, добавить, разыгрывать быстрее чем разыменование, добавление, разыменование).
Итак, новый код:
void MultiDiagonalSymmetricMatrix::CholeskyBackSolve(float *__restrict__ x){
// comments removed so that suggestions are more visible. Don't remove them in the real code!
// these definitions got long. Feel free to remove const; it does nothing for the optimiser
const float *const __restrict__ *const __restrict__ d = IncompleteCholeskyFactorization->Diagonals;
const unsigned int *const __restrict__ s = IncompleteCholeskyFactorization->StartRows;
const unsigned int M = IncompleteCholeskyFactorization->m;
const unsigned int N = IncompleteCholeskyFactorization->n;
unsigned int i;
unsigned int j;
for(j = 0; j < N; j++){ // don't use != as an optimisation; compilers can do more with <
float sub = 0;
for(i = 1; i < M && j >= s[i]; i++){
const unsigned int c = j - s[i];
sub += d[i][c]*x[c];
}
x[j] -= sub;
}
// Consider using processor-specific optimisations for this
#pragma omp parallel for
for(j = 0; j < N; j++){
x[j] /= d[0][j];
}
for( j = N; (j --) > 0; ){ // changed for clarity
float sub = 0;
for(i = 1; i < M && j + s[i] < N; i++){
sub += d[i][j]*x[j + s[i]];
}
x[j] -= sub;
}
}
Хорошо, что он выглядит более аккуратным, а нехватка памяти и сокращение разветвления, если не что иное, является стимулом. Если вы можете изменить s, чтобы добавить дополнительное значение UINT_MAX в конце, вы можете удалить больше ветвей (как проверки i<M, которые снова запускают ~ N * M раз).
Теперь мы не можем сделать больше циклов параллельно, и мы не можем комбинировать циклы. Усиление теперь будет, как было предложено в другом ответе, перестроить d. За исключением того, что работа, требуемая для переупорядочения d, имеет точно такие же проблемы с кешем, как и работа над циклом. И для этого потребуется выделенная память. Нехорошо. Единственными опциями для оптимизации являются следующие: изменить структуру самого IncompleteCholeskyFactorization->Diagonals, что, вероятно, будет означать много изменений или найти другой алгоритм, который лучше работает с данными в этом порядке.
Если вы хотите идти дальше, ваши оптимизации должны будут влиять на довольно много кода (не так уж плохо, если нет веской причины для Diagonals, являющегося массивом указателей, похоже, что это может быть связано с рефактору).