A 和 B 之间的矩阵-矩阵乘法，其中 AN 是一个随机矩阵答案

【问题标题】：Matrix-matrix multiplication between A and B, with A a random matrixA 和 B 之间的矩阵-矩阵乘法，其中 AN 是一个随机矩阵
【发布时间】：2013-09-11 19:59:19
【问题描述】：

我正在 CUDA 中研究一个特殊的矩阵-矩阵乘法 (AxB)，其中 A 是随机的 MxM 矩阵，B 是 MxN 矩阵。在下面的代码中，M 只是2000，但在实际情况下，它将被替换为大数字，以便为A 提供超过2GB 的矩阵。实际上，A 的所有元素都是随机的，并且限制在某个范围内，因此将由 randomize 函数生成。

我编写了下面的代码，其中A 的每个元素都是从数组中随机抽取的，因此原始AxB 将被修改为长度为M 乘以B 的向量。这是我编写代码的方式，但它似乎不起作用

#include <iostream>
#include <cusp/complex.h>

using namespace std;

#define M 2000
#define N 300

typedef cusp::complex<double> Complex;

__global__ void MVult(Complex* ad, Complex* bd, Complex* cd, int m1, int n1, int n2) 
{
  int x = (blockIdx.x * blockDim.x) + threadIdx.x;
  int y = (blockIdx.y * blockDim.y) + threadIdx.y;

  if(x < n2 && y < m1) 
  {   
    Complex sum = Complex(0.0, 0.0);
    int ridx = (rand()%(M-1)); // here I randomize the starting ridx 
    for(int i=0; i<n1; i++) sum += ad[ridx + i] * bd[i * n2 + x];
    cd[y * n2 + x] = v;
  }
}

int main(int argc, char *argv[])
{
  std::vector< Complex > _A(2*M+1);
  std::vector< Complex > _B(M*N);
  Complex *A, *B, *C;

  cudaMalloc((void**)&A, (2*M+1)*sizeof(Complex));
  cudaMalloc((void**)&B, M*N*sizeof(Complex));
  cudaMalloc((void**)&C, M*N*sizeof(Complex));

  for (int i=0; i<2*M+1; i++) _A[i] = Complex((double)i, (double)i);
  for (int i=0; i<M*N; i++) _B[i] = Complex(1.0, 0.0);

  cudaMemcpy( A, &_A[0], (2*M+1)*sizeof(Complex), cudaMemcpyHostToDevice );
  cudaMemcpy( B, &_B[0], (M*N)*sizeof(Complex), cudaMemcpyHostToDevice );

  dim3 block(32, 32);           
  dim3 grid((N+31)/32, (M+31)/32);

  MVult<<<grid, block>>>(A, B, C, M, M, N);
  cudaMemcpy(&_B[0], &C[0], (M*N)*sizeof(Complex), cudaMemcpyDeviceToHost);

  cudaFree(A);
  cudaFree(B);
  cudaFree(C);

  return 0;
}

我尝试使用 CPU 循环将其循环 M 次，并且每次运行向量和矩阵乘法（在 CUDA 中完成），但它太慢了。我正在寻找一种更快的方法来解决问题。

【问题讨论】：

标签： cuda matrix-multiplication

【解决方案1】：

您的代码会很慢，主要有两个原因：

在构造矩阵A 的方式中，您是在随机访问全局内存，因此防止合并访问；
正如@talonmies 在下面的评论中指出的那样，您正在实现自己的矩阵向量乘法例程，这肯定会比高度优化的例程如cuBLAS 慢。

为了加速你的代码，而不是使用你的__global__函数MVult，你可以/应该使用

cuRAND 用随机数填充矩阵A；
cuBLAS 执行A 和B 之间的矩阵乘法，特别是cublasCgemm() 用于单精度复杂计算。

如果矩阵A 太大，那么可以尝试将A*B 的计算分成更小的块，然后使用cuBLAS 的批处理功能（使用cublasSetStream()）尝试使用CUDA 实现并发执行流。

您可能还希望查看以下使用也推力的示例：

Matrix multiplication on GPU using CUDA with CUBLAS, CURAND and Thrust

正如@talonmies 所建议的那样，您可能还希望重新考虑您的方法。例如，如果A 是一个随机矩阵，那么A*B 也将是随机的。是否有可能利用A 的统计数据以及可能对B 的先验知识，通过随机方法直接构造矩阵A*B，而无需使用矩阵乘法？

【讨论】：

我已经尝试了 cuBLAS 和 cublasZgemm，但还是太慢了。我对 cublasSetStream 了解不多，有什么可以效仿的例子吗？我只在网上看到文档，没有看到完整的例子
@user12854197：您将无法编写比 CUBLAS 中已有的更快的例程。如果 CUBLAS 的 Zgemm 太慢，那么您需要重新考虑解决问题的方法，而不是尝试编写自己的复杂 gemm 内核
@user1285419 “太慢”是什么意思？它比你自己的方法慢吗？当A 是2Gb 矩阵时，您希望将A 和B 相乘的执行时间是多少，即浮点精度大约为20000x20000 元素的矩阵？
@user1285419 talonmies 还建议重新考虑您的方法。如果A 是一个随机矩阵，那么A*B 也将是随机的。是否有可能利用A 的统计数据和可能的先验知识B 直接通过随机方法构造矩阵A*B，而无需使用矩阵乘法？跨度>