В настоящее время я использую треугольный решатель для разреженных матриц и пытаюсь ускорить работу, используя директивы OpenACC.Учитывая мои матричные коэффициенты LU в разреженном формате CSR, OpenACC удалось правильно решить коэффициент L, но коэффициент U дает совершенно неверные результаты по сравнению с истинным решением приложения.Вот код ускоренного ядра для задачи обратной замены:
#pragma acc kernels deviceptr( ia, ja, factorValsU, y, x )
{
for ( int idx = size; idx > 0; idx-- )
{
double temp = 0.0;
int rowInit = ia[ idx - 1];
int rowEnd = ia[ idx ];
#pragma acc loop vector reduction( + : temp)
for ( int k = rowInit + 1; k < rowEnd; k++ )
{
temp += factorValsU[ k ] * x[ ja[ k ] ];
}
x[ idx ] = (y[ idx ] - temp) / factorValsU[ rowInit ];
}
}
Я не знаю, почему это ядро дает неверный результат.Я уже пробовал другую версию для ядра, где матрица сохраняется в обратном направлении, то есть снизу вверх, что в принципе можно решить с помощью следующего ядра:
#pragma acc kernels deviceptr( ia, ja, factorValsU, y, x )
{
for ( int idx = 0; idx < size; idx++ )
{
double temp = 0.0;
int rowInit = ia[ idx ];
int rowEnd = ia[ idx + 1 ];
#pragma acc loop vector reduction( + : temp)
for ( int k = rowInit + 1; k < rowEnd; k++ )
{
temp += factorValsU[ k ] * x[ ja[ k ] ];
}
x[ size - idx ] = (y[ size - idx ] - temp) / factorValsU[ rowInit ];
}
}
Но результат всегда неверен,Я пропустил что-то фундаментальное в декорировании обычного кода с помощью директив OpenACC для достижения правильных результатов?
Как уже упоминалось ранее, прямая замена L-фактора работает правильно, поэтому для полноты я опубликую код здесь.
#pragma acc kernels deviceptr( ia, ja, factorValsL, y, x )
{
for ( int idx = 0; idx < size; idx++ )
{
double temp = 0.0;
int rowInit = ia[ idx ];
int rowEnd = ia[ idx + 1 ];
#pragma acc loop vector reduction( + : temp)
for ( int k = rowInit; k < rowEnd; k++ )
{
temp += factorValsL[ k ] * x[ ja[ k ] ];
}
x[ idx ] = y[ idx ] - temp;
}
}
Обратите внимание на тонкую разницу между ядром для прямой замены (работает) и обратной замены (обе не работают), это область памяти, в которой сохранен результат:
x[ idx ] = y[ idx ] - temp for the L factor
x[ size - idx ] = (y[ size - idx ] - temp) / factorValsU[ rowInit ] for the U factor;
Есть ликакая-то причина для U-фактора, чтобы вычислить ошибочные результаты, причина порядка, в котором выполняется назначение (и лекция) в памяти?
Для полноты информации, предоставленной компилятором pgi18.4 относительно ядра:
triangularSolverU_acc(int, const int *, const int *, const double *, const double *, double *, bool):
614, Complex loop carried dependence of y->,x->,factorVals-> prevents parallelization
Loop carried dependence of x-> prevents parallelization
Loop carried backward dependence of x-> prevents vectorization
Accelerator kernel generated
Generating Tesla code
614, #pragma acc loop seq
621, #pragma acc loop vector(128) /* threadIdx.x */
Generating reduction(+:temp)
621, Loop is parallelizable
, который показывает, что внешний цикл был сериализован, а внутренний цикл является сокращением.