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]); 
}