AMP C++ вычисляет максимальное значение по массиву

Question

Я хочу сравнить максимальную абсолютную разницу между двумя трехмерными массивами с C++ AMP.

С OpenMP это просто. Учитывая 2 массивов

float*** u, ***uOld; 

Код является:

double residual = 0; 
#pragma omp parallel for schedule(static) collapse(3) reduction(max : residual) 
for (int i = 0; i < nbX; i++) 
    for (int j = 0; j < nbY; j++) 
     for (int k = 0; k < nbTheta; k++) 
      residual = std::max(residual, fabs(u[i][j][k] - uOld[i][j][k])); 

Было бы легко использовать max_element от AMP алгоритмов, но это не реализовано. Я думаю, что-то вроде этого, но сокращение требуется на уровне внешнего контура:

extent<1> extTheta(nbTheta); 
parallel_for_each(extTheta, [=, &u_, &uOld_](index<1> iTheta) restrict(amp) 
{ 
    type residual = 0; 
    for (int iX = 0; iX < nbX; iX++) 
    for (int iY = 0; iY < nbY; iY++) 
    residual = fast_math::fmax(residual, fast_math::fabs(u_[iX][iY][iTheta] - uOld_[iX][iY][iTheta])); 
}) 

Данные на GPU, и я не хочу его транзит на GPU по соображениям эффективности. Как это сделать эффективно?

Answer 1

Вот решение, вдохновленные из MSDN блога: https://blogs.msdn.microsoft.com/nativeconcurrency/2012/03/08/parallel-reduction-using-c-amp

parallel_for_each(extent<3>(nbTheta, nbX, nbY), [=, &u_, &uOld_](index<3> idx) restrict(amp) 
{ 
    uOld_[idx[0]][idx[1]][idx[2]] = abs1(u_[idx[0]][idx[1]][idx[2]] - uOld_[idx[0]][idx[1]][idx[2]]); 
}); 


array_view<float, 1> residualReduce_ = uOld_.view_as<1>(extent<1>(nbTheta*nbX*nbY)); 
array_view<float, 1> residual_ = residualReduce_.section(index<1>(0), extent<1>(1)); 
for (unsigned shift = nbTheta*nbX*nbY/2; shift > 0; shift /= 2) 
{ 
    parallel_for_each(extent<1>(shift), [=](index<1> idx) restrict(amp) 
    { 
     residualReduce_[idx[0]] = fast_math::fmax(residualReduce_[idx[0]], residualReduce_[idx[0] + shift]); 
     if (shift % 2){ //If odd, each thread includes a shifted entry. One will match the end of the queue 
      residualReduce_[idx[0]] = fast_math::fmax(residualReduce_[idx[0]], residualReduce_[idx[0] + shift + 1]); 
     } 
    }); 
} 
concurrency::copy(residual_, &residual); 
parallel_for_each(extent<3>(nbTheta, nbX, nbY), [=, &u_, &uOld_](index<3> idx) restrict(amp) 
{ 
    uOld_[idx[0]][idx[1]][idx[2]] = u_[idx[0]][idx[1]][idx[2]]; 
}) 

В отличие от сниппета в вопросе, решение включает в себя обновление Уолд на U.

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