使用 CUDA 并行的特征值求解器答案

【问题标题】：Eigenvalue solver in parallel using CUDA使用 CUDA 并行的特征值求解器
【发布时间】：2016-06-30 00:45:35
【问题描述】：

我一直在网上搜索和搜索，但似乎找不到我正在寻找的答案。我有一个特殊的问题。

我正在编辑这个以简化问题并希望它更具可读性和可理解性。

假设我有 5000 个 20x20 对称密集矩阵。我想在 CUDA 中创建一个内核，让每个线程负责计算每个对称矩阵的特征值。

如果可能的话，CUDA 内核的示例代码会很棒。

我们将不胜感激任何和所有的帮助/建议！

谢谢，

乔纳森

【问题讨论】：

Jacobi 也是迭代的。也许问题是您最初的猜测。如果您的猜测正确，您将需要更少的迭代来收敛到一个解决方案。你是只对特征值感兴趣，还是需要特征向量？
我只需要特征值。有关于并行执行雅可比的论文。这就是我上面所说的。对困惑感到抱歉。为了更好地解释我的情况，我事先进行了其他计算，以获取每个点的某些特征。然后我根据这些特征计算特征值。迭代部分只是计算那些特征值。我希望这有帮助。让我知道这是否有帮助！
我不确定我是否理解“点”的含义。这是方阵中行/列数的同义词吗？更大的矩阵意味着更多的计算——没有什么能改变这一点。
就好像我在处理 500 x 3 矩阵一样。 500 是点数（行），3 是维度（列）。本质上，我想将其扩展到 >100,000 x 3 矩阵或更多。
@Johnathan：请将所有这些内容编辑到您的问题中。正如所写的那样，它冗长而杂乱无章，很难理解。如果您对 N 的值有所了解，它可能也会有所帮助

标签： c++ c algorithm parallel-processing cuda

【解决方案1】：

我想在 CUDA 中创建一个内核，让每个线程负责计算每个对称矩阵的特征值。

我怀疑这是否是最快的方法，但它可能适用于非常小的矩阵。即使在这种情况下，也可能会进行一些数据存储优化（跨线程交错全局数据），但这会使事情复杂化。

如上所述，该请求可以映射为“令人尴尬的并行”算法，其中每个线程处理一个完全独立的问题。我们只需要找到合适的单线程“捐助者代码”。快速谷歌搜索后，我遇到了this。修改代码以这种独立于线程的方式运行是相当简单的。我们只需要借用 3 个例程（jacobi_eigenvalue、r8mat_diag_get_vector 和 r8mat_identity），并用__host__ __device__ 适当地装饰这些例程以在 GPU 上使用，同时没有其他更改。

有问题的代码似乎是佛罗里达州立大学的 J Burkardt 许可的 GNU LGPL。因此，考虑到这一点，并遵循conventional wisdom，我没有在这个答案中包含任何大量的代码。但是您应该能够使用我给出的说明通过实验重建我的结果。

注意：我不确定使用此代码的法律后果是什么，该代码声称已获得 GNU LGPL 许可。如果您选择使用此代码或其中的一部分，您应该确保遵守any necessary requirements。我在这里使用它的主要目的是演示单线程问题解决程序的相对微不足道的“令人尴尬的并行”扩展的概念。

通过转到 here 并将 3 个指示的函数复制粘贴到剩余代码骨架中指示的位置来重构我的完整代码应该是微不足道的。但这不会改变前面提到的任何通知/免责声明。使用它需要您自担风险。

同样，从性能的角度来看，不做其他更改可能不是最好的主意，但它会产生微不足道的工作量，并且可以作为一个可能有用的起点。一些可能的优化可能是：

寻找数据交错策略，使相邻线程更有可能读取相邻数据
从线程代码中消除new和delete函数，并用固定分配替换它（这很容易做到）
删除不必要的代码 - 例如计算和排序特征向量的代码（如果不需要该数据）

无论如何，使用上面修饰过的捐助者代码，我们只需要在它周围包裹一个简单的内核（je），以启动每个线程对单独的数据集（即矩阵）进行操作，并且每个线程产生自己的集合特征值（和特征向量 - 对于这个特定的代码库）。

出于测试目的，我将其设计为仅使用 3 个线程和 3 个 4x4 矩阵，但将它扩展到您希望的任意数量的矩阵/线程应该是微不足道的。

为了简洁起见，我已经放弃了the usual error checking，但我建议您使用它，或者如果您进行任何修改，至少使用cuda-memcheck 运行您的代码。

我还构建了代码来向上调整设备堆大小以适应内核内new 操作，具体取决于矩阵（即线程）的数量和矩阵尺寸。如果您进行了上面提到的第二个优化，您可能会删除它。

t1177.cu:

#include <stdio.h>
#include <iostream>
const int num_mat = 3; // total number of matrices = total number of threads
const int N = 4;   // square symmetric matrix dimension
const int nTPB = 256;  // threads per block

// test symmetric matrices

  double a1[N*N] = {
      4.0,  -30.0,    60.0,   -35.0, 
    -30.0,  300.0,  -675.0,   420.0, 
     60.0, -675.0,  1620.0, -1050.0, 
    -35.0,  420.0, -1050.0,   700.0 };

  double a2[N*N] = {
    4.0, 0.0, 0.0, 0.0, 
    0.0, 1.0, 0.0, 0.0, 
    0.0, 0.0, 3.0, 0.0, 
    0.0, 0.0, 0.0, 2.0 };

  double a3[N*N] = {
    -2.0,   1.0,   0.0,   0.0,
     1.0,  -2.0,   1.0,   0.0,
     0.0,   1.0,  -2.0,   1.0,
     0.0,   0.0,   1.0,  -2.0 }; 


/* ---------------------------------------------------------------- */
//
// the following functions come from here:
//
// https://people.sc.fsu.edu/~jburkardt/cpp_src/jacobi_eigenvalue/jacobi_eigenvalue.cpp
//
// attributed to j. burkardt, FSU
// they are unmodified except to add __host__ __device__ decorations
//
//****************************************************************************80
__host__ __device__
void r8mat_diag_get_vector ( int n, double a[], double v[] )
/* PASTE IN THE CODE HERE, FROM THE ABOVE LINK, FOR THIS FUNCTION */
//****************************************************************************80
__host__ __device__
void r8mat_identity ( int n, double a[] )
/* PASTE IN THE CODE HERE, FROM THE ABOVE LINK, FOR THIS FUNCTION */
//****************************************************************************80
__host__ __device__
void jacobi_eigenvalue ( int n, double a[], int it_max, double v[], 
  double d[], int &it_num, int &rot_num )
/* PASTE IN THE CODE HERE, FROM THE ABOVE LINK, FOR THIS FUNCTION */

// end of FSU code
/* ---------------------------------------------------------------- */

__global__ void je(int num_matr, int n, double *a, int it_max, double *v, double *d){

  int idx = threadIdx.x+blockDim.x*blockIdx.x;
  int it_num;
  int rot_num;
  if (idx < num_matr){
    jacobi_eigenvalue(n, a+(idx*n*n), it_max, v+(idx*n*n), d+(idx*n), it_num, rot_num);
  }
}

void initialize_matrix(int mat_id, int n, double *mat, double *v){

  for (int i = 0; i < n*n; i++) *(v+(mat_id*n*n)+i) = mat[i];
}

void print_vec(int vec_id, int n, double *d){

  std::cout << "matrix " << vec_id << " eigenvalues: " << std::endl;
  for (int i = 0; i < n; i++) std::cout << i << ": " << *(d+(n*vec_id)+i) << std::endl;
  std::cout << std::endl;
}
int main(){
// make sure device heap has enough space for in-kernel new allocations
  const int heapsize = num_mat*N*sizeof(double)*2;
  const int chunks = heapsize/(8192*1024) + 1;
  cudaError_t cudaStatus = cudaDeviceSetLimit(cudaLimitMallocHeapSize, (8192*1024) * chunks);
  if (cudaStatus != cudaSuccess) {
        fprintf(stderr, "set device heap limit failed!");
    }
  const int max_iter = 1000;
  double *h_a, *d_a, *h_v, *d_v, *h_d, *d_d;
  h_a = (double *)malloc(num_mat*N*N*sizeof(double));
  h_v = (double *)malloc(num_mat*N*N*sizeof(double));
  h_d = (double *)malloc(num_mat*  N*sizeof(double));
  cudaMalloc(&d_a, num_mat*N*N*sizeof(double));
  cudaMalloc(&d_v, num_mat*N*N*sizeof(double));
  cudaMalloc(&d_d, num_mat*  N*sizeof(double));
  memset(h_a, 0, num_mat*N*N*sizeof(double));
  memset(h_v, 0, num_mat*N*N*sizeof(double));
  memset(h_d, 0, num_mat*  N*sizeof(double));
  initialize_matrix(0, N, a1, h_a);
  initialize_matrix(1, N, a2, h_a);
  initialize_matrix(2, N, a3, h_a);
  cudaMemcpy(d_a, h_a, num_mat*N*N*sizeof(double), cudaMemcpyHostToDevice);
  cudaMemcpy(d_v, h_v, num_mat*N*N*sizeof(double), cudaMemcpyHostToDevice);
  cudaMemcpy(d_d, h_d, num_mat*  N*sizeof(double), cudaMemcpyHostToDevice);
  je<<<(num_mat+nTPB-1)/nTPB, nTPB>>>(num_mat, N, d_a, max_iter, d_v, d_d);
  cudaMemcpy(h_d, d_d, num_mat*N*sizeof(double), cudaMemcpyDeviceToHost);
  print_vec(0, N, h_d);
  print_vec(1, N, h_d);
  print_vec(2, N, h_d);
  return 0;
}

编译和示例运行：

$ nvcc -o t1177 t1177.cu
$ cuda-memcheck ./t1177
========= CUDA-MEMCHECK
matrix 0 eigenvalues:
0: 0.166643
1: 1.47805
2: 37.1015
3: 2585.25

matrix 1 eigenvalues:
0: 1
1: 2
2: 3
3: 4

matrix 2 eigenvalues:
0: -3.61803
1: -2.61803
2: -1.38197
3: -0.381966

========= ERROR SUMMARY: 0 errors
$

输出对我来说似乎是合理的，大部分与输出 here 匹配。

【讨论】：