Cuda, вычислить матрицу расстояний между 3d объектами

Question

У меня есть «строка» (молекула) связанных N объектов (атомов) в 3D (каждый атом имеет координаты). И мне нужно рассчитать расстояние между каждой парой атомов в молекуле (см. Ниже псевдокод). Как это можно сделать с CUDA? Должен ли я перейти к функции ядра 2 3D-массива? Или 3 массива с координатами: X [N], Y [N], Z [N]? Благодарю.

struct atom { double x, y, z; }

int main() 
{ 

//N number of atoms in a molecule 
double DistanceMatrix[N][N]; 
double d; 
atom Atoms[N]; 

for (int i = 0; i < N; i ++) 
    for (int j = 0; j < N; j++) 
     DistanceMatrix[i][j] = (atoms[i].x -atoms[j].x)*(atoms[i].x -atoms[j].x) + 
        (atoms[i].y -atoms[j].y)* (atoms[i].y -atoms[j].y) + (atoms[i].z -atoms[j].z)* (atoms[i].z -atoms[j].z; 

} 

Answer 1

Если вы работаете с очень большими молекулами, то, вероятно, не будет достаточно работы, чтобы держать GPU занят, поэтому расчеты будут быстрее с процессором.

Если Вы хотите вычислить Euclidean distance, то ваш расчет некорректный. Вам нужна 3D-версия теоремы Пифагора.

Я бы использовал SoA для хранения координат.

Вы хотите создать шаблон доступа к памяти с таким количеством коалесцированных чтений и записи, насколько это возможно. Для этого упорядочьте адреса или индексы, сгенерированные 32 потоками в каждом warp, как можно ближе друг к другу (немного упрощенное).

threadIdx обозначает индексы потоков внутри блока и blockIdx обозначает индексы блоков в сетке. blockIdx всегда одинаково для всех потоков в деформации. Только threadIdx изменяется в пределах потоков в блоке. Чтобы визуализировать, как 3 измерения threadIdx назначены потокам, считайте их вложенными циклами, где x является внутренним циклом, а z - это внешний цикл. Таким образом, потоки с соседними значениями x наиболее вероятно будут находиться в пределах одного и того же варпа, и если x делится на 32, только теки, имеющие одинаковое значение x/32, находятся в одном и том же варпе.

Я включил полный пример для вашего алгоритма ниже. В этом примере индекс i получен из , поэтому, чтобы проверить, что перекосы будут генерировать коалесцированные чтения и записи, я бы рассмотрел код, вставляя несколько последовательных значений, таких как 0, 1 и 2, для i и проверки того, что сгенерированные индексы также будут последовательными.

Адрес, сгенерированный из индекса j менее важен, как j происходит от threadIdx.y и поэтому менее вероятно, изменяться в пределах основы (и никогда не будет изменяться, если threadIdx.x делится на 32).

#include "cuda_runtime.h" 
#include <iostream> 

using namespace std; 

const int N(20); 

#define check(ans) { _check((ans), __FILE__, __LINE__); } 
inline void _check(cudaError_t code, char *file, int line) 
{ 
    if (code != cudaSuccess) { 
    fprintf(stderr,"CUDA Error: %s %s %d\n", cudaGetErrorString(code), file, line); 
    exit(code); 
    } 
} 

int div_up(int a, int b) { 
    return ((a % b) != 0) ? (a/b + 1) : (a/b); 
} 

__global__ void calc_distances(double* distances, 
    double* atoms_x, double* atoms_y, double* atoms_z); 

int main(int argc, char **argv) 
{ 
    double* atoms_x_h; 
    check(cudaMallocHost(&atoms_x_h, N * sizeof(double))); 

    double* atoms_y_h; 
    check(cudaMallocHost(&atoms_y_h, N * sizeof(double))); 

    double* atoms_z_h; 
    check(cudaMallocHost(&atoms_z_h, N * sizeof(double))); 

    for (int i(0); i < N; ++i) { 
    atoms_x_h[i] = i; 
    atoms_y_h[i] = i; 
    atoms_z_h[i] = i; 
    } 

    double* atoms_x_d; 
    check(cudaMalloc(&atoms_x_d, N * sizeof(double))); 

    double* atoms_y_d; 
    check(cudaMalloc(&atoms_y_d, N * sizeof(double))); 

    double* atoms_z_d; 
    check(cudaMalloc(&atoms_z_d, N * sizeof(double))); 

    check(cudaMemcpy(atoms_x_d, atoms_x_h, N * sizeof(double), cudaMemcpyHostToDevice)); 
    check(cudaMemcpy(atoms_y_d, atoms_y_h, N * sizeof(double), cudaMemcpyHostToDevice)); 
    check(cudaMemcpy(atoms_z_d, atoms_z_h, N * sizeof(double), cudaMemcpyHostToDevice)); 

    double* distances_d; 
    check(cudaMalloc(&distances_d, N * N * sizeof(double))); 

    const int threads_per_block(256); 
    dim3 n_blocks(div_up(N, threads_per_block)); 

    calc_distances<<<n_blocks, threads_per_block>>>(distances_d, atoms_x_d, atoms_y_d, atoms_z_d); 

    check(cudaPeekAtLastError()); 
    check(cudaDeviceSynchronize()); 

    double* distances_h; 
    check(cudaMallocHost(&distances_h, N * N * sizeof(double))); 

    check(cudaMemcpy(distances_h, distances_d, N * N * sizeof(double), cudaMemcpyDeviceToHost)); 

    for (int i(0); i < N; ++i) { 
    for (int j(0); j < N; ++j) { 
     cout << "(" << i << "," << j << "): " << distances_h[i + N * j] << endl; 
    } 
    } 

    check(cudaFree(distances_d)); 
    check(cudaFreeHost(distances_h)); 
    check(cudaFree(atoms_x_d)); 
    check(cudaFreeHost(atoms_x_h)); 
    check(cudaFree(atoms_y_d)); 
    check(cudaFreeHost(atoms_y_h)); 
    check(cudaFree(atoms_z_d)); 
    check(cudaFreeHost(atoms_z_h)); 

    return 0; 
} 

__global__ void calc_distances(double* distances, 
    double* atoms_x, double* atoms_y, double* atoms_z) 
{ 
    int i(threadIdx.x + blockIdx.x * blockDim.x); 
    int j(threadIdx.y + blockIdx.y * blockDim.y); 

    if (i >= N || j >= N) { 
    return; 
    } 

    distances[i + N * j] = 
    (atoms_x[i] - atoms_x[j]) * (atoms_x[i] - atoms_x[j]) + 
    (atoms_y[i] - atoms_y[j]) * (atoms_y[i] - atoms_y[j]) + 
    (atoms_z[i] - atoms_z[j]) * (atoms_z[i] - atoms_z[j]); 
}