Как эффективно вычислять функции в трехмерной сетке с помощью распараллеливания графического процессора? ⇐ C++

[Anonymous]

Я работаю над проектом, в котором мне нужно вычислить функцию в каждой точке 3D-сетки. Изначально я использовал вложенный цикл для перебора всех точек и выполнения вычислений. Этот подход был многопоточным, в котором работа распределялась по оси Z, где каждый поток обрабатывал часть сетки:

Код: Выделить всё

// Iterate over the Z-axis of the 3D grid
for (int z = 0; z < dim[2]; ++z) {
// Iterate over the Y-axis of the 3D grid
for (int y = 0; y < dim[1]; ++y) {
// Iterate over the X-axis of the 3D grid
for (int x = 0; x < dim[0]; ++x) {
// Create a 3D point (pt) to represent the current position in the grid
vec3 pt;
// Calculate the X coordinate of the point as a normalized value between 0 and 1
pt[0] = (static_cast(x) / static_cast(dim[0]));
// Calculate the Y coordinate of the point as a normalized value between 0 and 1
pt[1] = (static_cast(y) / static_cast(dim[1]));
// Calculate the Z coordinate of the point as a normalized value between 0 and 1
pt[2] = (static_cast(z) / static_cast(dim[2]));
// Compute the value of the function at the current point
double val = computeFunction(pt);
// Store the computed value in the results vector
results.push_back(val);
}
}
}

Теперь, конечно, функцию ComputeFunction() необходимо перенести на CUDA, но я хочу сосредоточиться на логике эффективного разделения трехмерной сетки точек на блоки точек в параллельно.


Я думал об этой реализации, которая возвращает одномерный массив:

Код: Выделить всё

__device__ double computeFunction(double x, double y, double z) {
return x * y * z; // Example function
}

__global__ void compute3dGrid(double* results, int dimX, int dimY, int dimZ) {
int x = blockIdx.x * blockDim.x + threadIdx.x;
int y = blockIdx.y * blockDim.y + threadIdx.y;
int z = blockIdx.z * blockDim.z + threadIdx.z;

if (x < dimX && y < dimY && z < dimZ) {
double pt_x = static_cast(x) / static_cast(dimX);
double pt_y = static_cast(y) / static_cast(dimY);
double pt_z = static_cast(z) / static_cast(dimZ);

double val = computeFunction(pt_x, pt_y, pt_z);

int idx = z * dimX * dimY + y * dimX + x;
results[idx] = val;
}
}

void launchGridComputation(double* h_results, int dimX, int dimY, int dimZ) {
double* d_results;
size_t size = dimX * dimY * dimZ * sizeof(double);

cudaMalloc(&d_results, size);

dim3 blockSize(8, 8, 8);
dim3 gridSize(
(dimX + blockSize.x - 1) / blockSize.x,
(dimY + blockSize.y - 1) / blockSize.y,
(dimZ + blockSize.z - 1) / blockSize.z
);

compute3dGrid(d_results, dimX, dimY, dimZ);

cudaMemcpy(h_results, d_results, size, cudaMemcpyDeviceToHost);

cudaFree(d_results);
}

В конечном итоге эта проблема может потребовать метода проб и ошибок для оптимизации размеров блоков?
Есть ли лучший способ сделать это?< /strong>
ОБНОВЛЕНИЕ:
Я хотел бы как можно больше ответить на различные вопросы, которые были опубликованы в комментариях к исходному сообщению. :
Вот код, который выполняет вычисления в каждой точке:

Код: Выделить всё

double Molecule::getFastNeutralSphericalDensity(QList < Atom* >* list, vec3 point)
{
if (list == nullptr) return(0.0);
if (list->isEmpty()) return(0.0);

Function_Result ret;                                // Multi Data declaration
ret.Data_1d = 0.0;                                  // Multi data Initialisation
ret.Data_2d = 0.0;
ret.Data_3d = 0.0;
ret = getFastNeutralSphericalDensity_MultiData(list, point);    // Multi data function call
return(ret.Data_1d);
}

Function_Result Molecule::getFastNeutralSphericalDensity_MultiData(QList < Atom* >* list, vec3 point)
{
Function_Result ret;                                // Multi Data declaration
ret.Data_1d = 0.0;                                  // Multi data Initialisation
ret.Data_2d = 0.0;
ret.Data_3d = 0.0;

if (list == nullptr) return(ret);
if (list->isEmpty()) return(ret);

double rho, rhototal, imin, imax, valmin, valmax, disttonuclei, dmin, incr;
Atom *atm = nullptr;
QList radialdensity;
int chem;
int rsize;

dmin = 5.0;
incr = 0.05;
rhototal = 0.0;
double occmult, mult;

QList::iterator it;
for (it = list->begin(); it != list->end(); ++it)
{
atm = *it;

if (atm->isDummy()) continue;

disttonuclei = arma::norm(point-atm->getCoords());
//if (disttonuclei > dmin) continue;

chem = atm->getZDiff();
mult = atm->getMult();
if (mult getOccupancy()/mult;
radialdensity = mol_radial_densities.value(chem);

if (radialdensity.isEmpty())
{
qWarning() = rsize) continue;
valmin = radialdensity.at(static_cast(imin));
valmax = radialdensity.at(static_cast(imax));
rho = interpol_lineaire(valmin, valmax, disttonuclei, imin, imax);
}
rhototal += occmult*rho;
}

ret.Data_1d = rhototal;                                       // Multi data result
ret.Data_2d = 0.0;
ret.Data_3d = 0.0;
return(ret);
}

Функция getFastNeutralSphericalDensity и ее помощник getFastNeutralSphericalDensity_MultiData вычисляют плотность в заданной точке, суммируя вклады каждого атома в заданном диапазоне. Этот тип вычислений часто используется в компьютерной химии и моделировании молекулярной динамики для моделирования распределения электронной плотности или других скалярных полей вокруг атомов в молекуле.
Например, одна из молекул, которые мы работают над созданием сетки из 91 миллиона точек, выполнение которой на многопоточном процессоре занимает несколько часов.
Я успешно создал реализацию графического процессора, которая выполняется примерно за 60–80 секунд. 91 млн баллов. НО[/b] У меня возникла проблема с использованием памяти:



Платформа
Обработанные баллы
Время выполненияИспользование памяти (ОЗУ)




< td>ЦП
91M
6 часов
4G (постоянно)


ГП
91M
70,77 с
До 16G в системе 16G



Что самое странное, так это то, что я почти не вижу использования графического процессора в диспетчере задач, а также использования памяти графического процессора. Использование оперативной памяти зашкаливает, однако при использовании метода графического процессора оно составляет около 16 ГБ и остается постоянным на уровне около 7,4 ГБ после завершения вычислений.
Я все еще исследую, почему использование памяти настолько велико .



Подробнее здесь: https://stackoverflow.com/questions/787 ... lelization