在 CUDA 中将矢量化二维方阵和压缩三对角矩阵相乘答案

【问题标题】：Multiply vectorized 2D square matrix and compressed tridiagonal matrix in CUDA在 CUDA 中将矢量化二维方阵和压缩三对角矩阵相乘
【发布时间】：2015-06-26 13:03:16
【问题描述】：

我有两个矩阵

#define MATRIX_SIZE 20
#define BLOCK_SIZE 2
#define TILE_SIZE  2

double** A
double** B

矩阵 A 是稠密的，矩阵 B 是三对角的。我创建了 A 的矢量化表示

/* sz = A.rowlen = B.rowlen = A.collen = B.collen */
double* A1d = matrix_to_vector(sz, A);

我还使用以下函数创建了 B 的压缩表示

double* l_array = new double(sz - 1);
double* m_array = new double(sz);
double* r_array = new double(sz-1);
int current_l_idx = 0;
int current_m_idx = 0;
int current_r_idx = 0;

for (int i = 0; i < sz; i++) {
  for (int j = 0; j < sz; j++) {
    if ((i == j+1) || (i-1 == j)) {
      l_array[current_l_idx] = B[i][j];
      current_l_idx++;
    }
    else if ((i == j-1) || (i+1 == j)) {
      r_array[current_r_idx] = B[i][j];
      current_r_idx++;
    }
    else if (i == j) {
      m_array[current_m_idx] = B[i][j];
      current_m_idx++;
    }
  }
}

然后，我为CUDA 创建一个空的二维向量化矩阵 E 以及我的所有对象

double* E1d = matrix_to_vector(sz, E);

double* d_A
double* d_B_l;
double* d_B_m;
double* d_B_r;
double* d_E;

size_t sizeA = sz * sz * sizeof(double);
size_t sizeB_lr = (sz - 1) * sizeof(double);
size_t sizeB_m = sz * sizeof(double);

cudaMalloc(&d_A, sizeA);
cudaMalloc(&d_B_l. sizeB_lr);
cudaMalloc(&d_B_m, sizeB_m);
cudaMalloc(&d_B_r, sizeB_lr);
cudaMalloc(&d_E, sizeA);

cudaMemcpy(d_A, A1d, sizeA, cudaMemcpyHostToDevice);
cudaMemcpy(d_B_l, l_array, sizeB_lr, cudaMemcpyHostToDevice);
cudaMemcpy(d_B_m, m_array, sizeB_m, cudaMemcpyHostToDevice);
cudaMemcpy(d_B_r, r_array, sizeB_lr, cudaMemcpyHostToDevice);
cudaMemcpy(d_E, E1d, sizeA, cudaMemcpyHostToDevice);

dim3 threads(BLOCK_SIZE, BLOCK_SIZE);
dim3 grid(MATRIX_SIZE / threads.x, MATRIX_SIZE / threads.y);

cudakernel<<<grid, threads>>>(sz, d_A, d_B_l, d_B_m, d_B_r, d_E);

我可以连续执行这个乘法，但不幸的是，我不知道如何在 CUDA 设备上实现它

假设

A 和 B 总是正方形的
sz 总是能被 BLOCK_SIZE 和 TILE_SIZE 整除
BLOCK_SIZE 将始终等于 TILE_SIZE

【问题讨论】：

如果你愿意将你的三对角矩阵转换为 CSR 存储，你可以使用 CUSPARSE csrmm 因为你的矩阵是正方形的，你想要的操作（密集 x 稀疏）可以通过转置来实现操作并反转 csrmm 的顺序（稀疏 x 密集）。
@Robert CSR 是可以接受的，但不幸的是我不能使用库。 :(
矩阵乘法计算输出矩阵 (C) 中每个点的行 (A) 乘以 (B) 列的总和。由于 B 是您的三对角矩阵，这意味着您在 B 的列中对于每个输出元素只有 3 个元素（或者对于边缘情况为 2 个）。一个简单的方法可能是适当地限制 for 循环（只做 3 次乘法）并提出适当的索引来索引到您的 B 向量中。
@Robert 我同意。这就是向我建议的方法。我只是不太擅长设备上的矩阵迭代。 :( 你能提供如何使用这种方法编写这个内核吗？
这是做作业的吗？如果您了解如何进行密集矩阵-矩阵乘法，您应该能够弄清楚在三对角情况下每个元素需要哪些 3 个乘积项。弄清楚这一点（C 中的每个结果都需要 A 和 B 的 3 个乘积项）实际上与 GPU 无关。

标签： c++ matrix cuda sparse-matrix matrix-multiplication

【解决方案1】：

根据您的设置代码，我怀疑您正在寻找这种矩阵乘法的平铺共享内存方法，而我并不想为您做功课，所以我将演示一个示例不使用共享内存。

如果您了解矩阵乘法的工作原理，并且您还了解如何创建普通的shared memory GPU matrix multiply kernel，那么将以下代码转换为使用共享内存应该相对简单：

#include <stdio.h>

#define DSIZE 256
#define BSIZE 32
#define TOL 0.0001
typedef  double mytype;

#define cudaCheckErrors(msg) \
    do { \
        cudaError_t __err = cudaGetLastError(); \
        if (__err != cudaSuccess) { \
            fprintf(stderr, "Fatal error: %s (%s at %s:%d)\n", \
                msg, cudaGetErrorString(__err), \
                __FILE__, __LINE__); \
            fprintf(stderr, "*** FAILED - ABORTING\n"); \
            exit(1); \
        } \
    } while (0)




// C = A x B
// A,B,C are all dense
template <typename T>
__global__ void mm(const T * __restrict__ A, const T * __restrict__ B, T * __restrict__ C, const int sz){

  int idx = threadIdx.x+blockDim.x*blockIdx.x;
  int idy = threadIdx.y+blockDim.y*blockIdx.y;

  if ((idx < sz) && (idy < sz)){
    T temp = 0;
    for (int i = 0; i < sz; i++)
      temp += A[idy*sz+i]*B[i*sz+idx];
    C[idy*sz+idx] = temp;}
}


// C = A x B
// A,C are dense, B is tridiagonal
template <typename T>
__global__ void mmt(const T * __restrict__ A, const T * __restrict__ B_l, const T * __restrict__ B_m, const T * __restrict__ B_r, T * __restrict__ C, const int sz){

  int idx = threadIdx.x+blockDim.x*blockIdx.x;
  int idy = threadIdx.y+blockDim.y*blockIdx.y;

  if ((idx < sz) && (idy < sz)){
    T temp = 0;
    if (idx > 0)      temp += A[idy*sz+(idx-1)]*B_r[idx-1];
                      temp += A[idy*sz+(idx)  ]*B_m[idx];
    if (idx < (sz-1)) temp += A[idy*sz+(idx+1)]*B_l[idx];
    C[idy*sz+idx] = temp;}
}



int main(){

  mytype *d_A, *h_A, *d_B, *h_B, *d_C, *h_Cd, *h_Cs, *d_B_l, *h_B_l, *d_B_m, *h_B_m, *d_B_r, *h_B_r;
  size_t msz  = DSIZE*DSIZE;
  size_t mszb = msz*sizeof(mytype);
// host side allocations
  h_A = (mytype *)malloc(mszb);
  h_B = (mytype *)malloc(mszb);
  h_Cd =(mytype *)malloc(mszb);
  h_Cs =(mytype *)malloc(mszb);
  h_B_l = (mytype *)malloc((DSIZE-1)*sizeof(mytype));
  h_B_r = (mytype *)malloc((DSIZE-1)*sizeof(mytype));
  h_B_m = (mytype *)malloc(    DSIZE*sizeof(mytype));
  if (!h_A || !h_B || !h_Cd || !h_Cs || !h_B_l || !h_B_r || !h_B_m) {printf("malloc fail\n"); return -1;}
// device side allocations
  cudaMalloc(&d_A, mszb);
  cudaMalloc(&d_B, mszb);
  cudaMalloc(&d_C, mszb);
  cudaMalloc(&d_B_l, (DSIZE-1)*sizeof(mytype));
  cudaMalloc(&d_B_r, (DSIZE-1)*sizeof(mytype));
  cudaMalloc(&d_B_m,     DSIZE*sizeof(mytype));
  cudaCheckErrors("cudaMalloc fail");
// prepare A, B matrices

/*
       |1 1 1 ...|
   A = |2 2 2 ...|
       |3 3 3 ...|
       |4 4 4 ...|
       |...      |


       |2 1 0 ...|   B_l = left/lower subdiagonal    (i.e. all 3's)
   B = |3 2 1 ...|   B_m = middle/main diagonal      (i.e. all 2's)
       |0 3 2 ...|   B_r = right/upper superdiagonal (i.e. all 1's)
       |0 0 3 ...|
       |...      |
*/



  for (int i = 0; i < DSIZE; i++){
    if (i < DSIZE-1){
      h_B_r[i] = 1;
      h_B_l[i] = 3;}
    h_B_m[i] = 2;
    for (int j = 0; j < DSIZE; j++){
      h_A[i*DSIZE+j] = i+1;
      if (j==i+1)      h_B[i*DSIZE+j] = 1;
      else if (j==i)   h_B[i*DSIZE+j] = 2;
      else if (j==i-1) h_B[i*DSIZE+j] = 3;
      else             h_B[i*DSIZE+j] = 0;}}

// copy data to device

  cudaMemcpy(d_A, h_A, mszb, cudaMemcpyHostToDevice);
  cudaMemcpy(d_B, h_B, mszb, cudaMemcpyHostToDevice);
  cudaMemcpy(d_B_l, h_B_l, (DSIZE-1)*sizeof(mytype), cudaMemcpyHostToDevice);
  cudaMemcpy(d_B_r, h_B_r, (DSIZE-1)*sizeof(mytype), cudaMemcpyHostToDevice);
  cudaMemcpy(d_B_m, h_B_m,     DSIZE*sizeof(mytype), cudaMemcpyHostToDevice);
  cudaCheckErrors("cudaMemcpy1 fail");

// perform dense-dense multiply

  dim3 block(BSIZE,BSIZE);
  dim3  grid((DSIZE+block.x-1)/block.x, (DSIZE+block.y-1)/block.y);

  cudaMemset(d_C, 0, mszb);
  mm<<<grid, block>>>(d_A, d_B, d_C, DSIZE);
  cudaMemcpy(h_Cd, d_C, mszb, cudaMemcpyDeviceToHost);
  cudaCheckErrors("cudaMemcpy 2/kernel fail");

// perform dense-sparse multiply

  cudaMemset(d_C, 0, mszb);
  mmt<<<grid, block>>>(d_A, d_B_l, d_B_m, d_B_r, d_C, DSIZE);
  cudaMemcpy(h_Cs, d_C, mszb, cudaMemcpyDeviceToHost);
  cudaCheckErrors("cudaMemcpy 3/kernel fail");

// compare results

  for (int i = 0; i < DSIZE; i++)
    for (int j = 0; j < DSIZE; j++)
      if (abs(h_Cs[i*DSIZE+j] - h_Cd[i*DSIZE+j]) > TOL) {printf("results mismatch at (%d, %d) dense: %f sparse: %f\n", i, j, h_Cd[i*DSIZE+j], h_Cs[i*DSIZE+j]); return -1;}
  printf("Success!\n");

  return 0;
}

注意事项：

mmt 内核中的所有全局内存访问（即A、B 向量和C）都应该正确地跨线程合并。因此，使用共享内存的转换也应该很容易产生对共享内存的非银行冲突访问。
虽然研究此代码可能对学习有用，但我建议使用 CUSPARSE 的例程（例如 csrmm）完成任何严重的稀疏密集矩阵乘法。它几乎肯定会比上述代码更有效（更快），并且可能也比上述代码的任何共享内存转换更快。

【讨论】：