实现 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 中。它应该可以很好地扩展 row/column 计数。

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

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

#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

两个答案都是可行的解决方案。问题是，better/faster?

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

Best way to achieve CUDA Vector Diagonalization

cuda

matrix