CUDA，计算3D对象之间的距离矩阵

Question

我有一个3D连接的N个对象（原子）的“串”（分子）（每个原子都有一个坐标）。我需要计算分子中每对原子之间的距离（参见下面的伪代码）。 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忙，所以计算会随着CPU速度更快。

如果您打算计算欧几里德距离，您的计算是不正确的。你需要毕达哥拉斯定理的3D版本。

我会使用SoA来存储坐标。

您想要生成一个内存访问模式，尽可能多的合并读写操作。为此，将每个warp中由32个线程生成的地址或索引尽可能地彼此靠近（稍微简化一点）。

threadIdx指定块内的线索索引并且blockIdx指定网格内的块索引。 blockIdx对于变形中的所有线程总是相同的。只有threadIdx在块内的线程内变化。要想象如何将threadIdx的3个维度分配给线程，请将它们视为嵌套循环，其中x是内部循环，而z是外部循环。因此，具有相邻x值的线程最有可能位于同一翘曲内，并且如果x可被32整除，则只有共享相同x/32值的线程处于相同翘曲内。

我在下面为您的算法包含了一个完整的示例。在这个例子中，i索引是从threadIdx.x派生的，因此，为了检查经线是否会产生合并读写，我会在插入一些连续值（例如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]); 
}