Оптимизированное объединение для 3D конечных разностей с использованием CUDA

Я написал ядро CUDA для выполнения многочисленных операций (включая конечные разности) на трехмерной прямоугольной карте с размером сетки [N_x, N_y, N_z] и равномерным шагом сетки. Я разделяю область на трехмерные блоки размером [B_x, B_y, B_z] и группирую рабочие элементы так, что каждая рабочая группа вычисляет конечные разности для 3D блока. Я прочитал ряд источников о том, как лучше подготовить свое ядро, чтобы использовать слияние и общий доступ к памяти для повышения пропускной способности. На этих данных должно быть выполнено несколько ядер, поэтому я придумал идею, как упорядочить данные, чтобы максимизировать слияние между потоками. Мне хотелось бы услышать ваше мнение о том, насколько эта идея абсурдна.

Суть идеи заключается в том, чтобы упорядочить данные в формате “блока”, так что, когда я вызываю (например, FD) ядро на этих данных, импорт в общую память рабочей группы будет слиянием. Ядро для преобразования из стандартного формата строк в блочное упорядочение данных выглядит следующим образом:

__constant__ int NX, NY, NZ, NT;      // Глобальные размеры области
__constant__ int BX, BY, BZ, BT;      // Размеры блока
__constant__ int NBX, NBY, NBZ;       // Количество блоков в области

// gid_rm(...) принимает глобальные индексы x, y, z и вычисляет "монолитный" ид узла
__device__\n
int gid_rm(const int i,const int j, const int k, const int DX,const int DY, const int DZ){ 
return i*DY*DZ + j*DZ + k;
} 

// gid_bl(...) принимает глобальные индексы x, y, z и вычисляет "блоковый" ид узла
__device__
int gid_bl(const int im,const int jm, const int km){ 
const int ib = im/BX, jb = jm/BY, kb = km/BZ;
const int gl = gid(im-ib*BX, jm-jb*BY, km-kb*BZ, BX, BY, BZ);
const int bls = gid(ib,jb,kb,NBX,NBY,NBZ);
return gl + BX*BY*BZ*bls; 
}

template<typename T>
__global__
void Data2BlockFormat(const T* src, T* dst) {
    const int gx = threadIdx.x + blockIdx.x*BX; 
    const int gy = threadIdx.y + blockIdx.y*BY; 
    const int gz = threadIdx.z + blockIdx.z*BZ; 
    const int mid = gid_rm(gx,gy,gz,NX,NY,NZ);          
    const int bid = gid_bl(gx,gy,gz); 
    dst[bid]       = src[mid];
    dst[bid+NT]    = src[mid+NT];
    dst[bid+2*NT]  = src[mid+2*NT];
}";

Я упрощенно представил все входные массивы для повышения эффективности копирования данных. Операция над данными с блочным упорядочением будет выглядеть примерно так:

template<int NFDGrid, typename T>
__global__
void Div(const T* src, const int* halo, T* grad) {

// Объявление массивов общей памяти (паддинг массивов для учета гало)
__shared__ T u[NFDGrid];
__shared__ T v[NFDGrid];
__shared__ T w[NFDGrid];

// Подготовка индексов
const int gx = threadIdx.x + gx0;       // Глобальные индексы 
const int gy = threadIdx.y + gy0; 
const int gz = threadIdx.z + gz0; 
const int txh = threadIdx.x + Halo;     // Локальные индексы в паддинговом массиве
const int tyh = threadIdx.y + Halo;
const int tzh = threadIdx.z + Halo;

const int bl_id = gidb(gx,gy,gz);

// Загрузка "ядерных" значений блока в общую память (должно срежиссировать!)
u[pid] = f[bl_id];
v[pid] = f[bl_id+NT];
w[pid] = f[bl_id+2*NT];
__syncthreads();

// Загрузка значений гало в общую память (не может быть срежиссировано)
(...)

// Выполнение операций над общей памятью
T dudx=0., dvdy=0., dwdz=0.;
(...)

// Установка выходных данных в блочном упорядочении данных
grad[bl_id]      = dudx;
grad[bl_id+NT]   = dvdy;
grad[bl_id+2*NT] = dwdz;

Мой вопрос следующий: Если я выберу блок размером, например, [B_x, B_y, B_z] = [8,8,8], произойдет ли идеальное слияние для 512 элементов блока, или данные будут загружаться всего лишь в последовательности по 32 бита за варп? Имеет ли более смысл структурировать импорт данных как 32-битные карандаши и импортировать их с помощью for цикла? Если есть общий лучший подход к этой проблеме, я был бы рад узнать о нем!

Я пытался выполнить ту же операцию с карандашами вдоль наиболее быстро вращающегося индекса, и мой подход по-прежнему кажется быстрее, однако это может быть связано с наивной реализацией моих ядер FD.