实现 CUDA 矢量对角化的最佳方法

Question

我想要做的是输入我的 m x n 矩阵，并并行地为矩阵的每一列构造 n 个方对角矩阵，对每个方对角矩阵执行操作，然后重新组合结果。我该怎么做呢？

到目前为止，我从一个 m x n 矩阵开始；先前矩阵计算的结果，其中每个元素都是使用函数 y = f(g(x)) 计算的。

这给了我一个包含 n 列元素 [f1, f2...fn] 的矩阵，其中每个 fn 代表一个高度为 m 的列向量。

从这里开始，我想区分矩阵的每一列关于 g(x)。区分 fn(x) w.r.t. g(x) 产生一个包含元素 f'(x) 的方阵。在约束条件下，该方阵简化为雅可比行列式，每行的元素沿方阵的对角线，且等于 fn'，所有其他元素均为零。

这就是为什么需要为每个向量行 fn 构造对角线的原因。

为此，我采用定义为 A(hA x 1) 的目标向量，该向量是从较大的 A(m x n) 矩阵中提取的。然后我准备了一个定义为 C(hA x hA) 的归零矩阵，用于保存对角线。

目的是将向量 A 对角化为一个方阵，其中 A 的每个元素都位于 C 的对角线上，其他所有元素都为零。

可能有更有效的方法可以使用一些预先构建的例程来完成此任务，而无需构建全新的内核，但请注意，出于这些目的，此方法是必要的。

完成此操作的内核代码（有效）如下所示：

_cudaDiagonalizeTest << <5, 1 >> >(d_A, matrix_size.uiWA, matrix_size.uiHA, d_C, matrix_size.uiWC, matrix_size.uiHC);

__global__ void _cudaDiagonalizeTest(float *A, int wA, int hA, float *C, int wC, int hC)
{
    int ix, iy, idx;

    ix = blockIdx.x * blockDim.x + threadIdx.x;
    iy = blockIdx.y * blockDim.y + threadIdx.y;

    idx = iy * wA + ix;

    C[idx * (wC + 1)] = A[idx];

}

我有点怀疑这是一种非常幼稚的解决方案，并且想知道是否有人可以举例说明我如何使用

a) 减少

b) 推力

对于大行大小的向量，我希望能够使用 GPU 的多线程功能将任务分成小作业，并在最后将每个结果与 __syncthreads() 组合。

下图显示了想要的结果。

我已阅读NVIDIA's article on reduction，但未能达到预期的效果。

非常欢迎任何帮助或解释。

谢谢。

矩阵 A 是具有 4 列的目标。我想获取每一列，并将其元素作为对角线复制到矩阵 B 中，遍历每一列。

Answer 1

我创建了一个基于推力的简单示例。它使用列优先顺序将矩阵存储在thrust::device_vector 中。它应该可以很好地适应更大的行/列数。

另一种方法可能基于thrust strided_range example。

这个例子做你想做的（根据输入向量填充对角线）。但是，根据您如何将生成的矩阵进行到“微分”步骤，可能仍然值得研究稀疏存储（没有所有零条目）是否可能，因为这将减少内存消耗并简化迭代。

#include <thrust/device_vector.h>
#include <thrust/scatter.h>
#include <thrust/sequence.h>
#include <thrust/iterator/transform_iterator.h>
#include <thrust/iterator/counting_iterator.h>
#include <thrust/functional.h>
#include <iostream>


template<typename V>
void print_matrix(const V& mat, int rows, int cols)
{
   for(int i = 0; i < rows; ++i)
   {
     for(int j = 0; j < cols; ++j)
     {
      std::cout << mat[i + j*rows] << "\t";
     }
     std::cout << std::endl;
   }
}

struct diag_index : public thrust::unary_function<int,int>
{
  diag_index(int rows) : rows(rows){}

  __host__ __device__
  int operator()(const int index) const
  {
      return (index*rows + (index%rows));
  }

  const int rows;
};

int main()
{
  const int rows = 5; 
  const int cols = 4;

  // allocate memory and fill with demo data
  // we use column-major order
  thrust::device_vector<int> A(rows*cols);
  thrust::sequence(A.begin(), A.end());

  thrust::device_vector<int> B(rows*rows*cols, 0);

  // fill diagonal matrix
  thrust::scatter(A.begin(), A.end(), thrust::make_transform_iterator(thrust::make_counting_iterator(0),diag_index(rows)), B.begin());

  print_matrix(A, rows, cols);
  std::cout << std::endl;
  print_matrix(B, rows, rows*cols);
  return 0;
}

这个例子会输出：

0    5    10    15    
1    6    11    16    
2    7    12    17    
3    8    13    18    
4    9    14    19    

0    0    0    0    0    5    0    0    0    0    10    0    0    0    0    15    0    0    0    0    
0    1    0    0    0    0    6    0    0    0    0    11    0    0    0    0    16    0    0    0    
0    0    2    0    0    0    0    7    0    0    0    0    12    0    0    0    0    17    0    0    
0    0    0    3    0    0    0    0    8    0    0    0    0    13    0    0    0    0    18    0    
0    0    0    0    4    0    0    0    0    9    0    0    0    0    14    0    0    0    0    19    

Answer 2

不使用推力的替代答案如下：

_cudaMatrixTest << <5, 5 >> >(d_A, matrix_size.uiWA, matrix_size.uiHA, d_C, matrix_size.uiWC, matrix_size.uiHC);

__global__ void _cudaMatrixTest(float *A, int wA, int hA, float *C, int wC, int hC)
{
    int ix, iy, idx;

    ix = blockIdx.x * blockDim.x + threadIdx.x;
    iy = blockIdx.y * blockDim.y + threadIdx.y;

    idx = iy * wA + ix;

    C[idx * wC + (idx % wC)] = A[threadIdx.x * wA + (ix / wC)];
}

d_A 在哪里

0    5    10    15    
1    6    11    16    
2    7    12    17    
3    8    13    18    
4    9    14    19    

这两个答案都是可行的解决方案。问题是，哪个更好/更快？