c++ - Cuda，计算3d对象之间的距离矩阵

Question

我有一个 3D 连接的 N 个对象(原子)的“字符串”(分子)(每个原子都有一个坐标)。我需要计算分子中每对原子之间的距离(参见下面的伪代码)。用 CUDA 怎么办？我应该将 2 个 3D 数组传递给内核函数吗？或者 3 个坐标数组:X[N]、Y[N]、Z[N]？谢谢。

结构原子 { 双 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用于存储坐标。

您想生成一个内存访问模式，其中包含尽可能多的 union 读写。为此，安排由每个 warp 中的 32 个线程生成的地址或索引尽可能彼此靠近(稍微简化)。

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

我在下面为您的算法提供了一个完整的示例。在示例中，i 索引是从 threadIdx.x 派生的，因此，为了检查扭曲是否会生成合并的读取和写入，我将检查代码，同时插入一些i 的连续值，例如 0、1 和 2，并检查生成的索引是否也是连续的。

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