使用 CUBLAS 进行批量 QR 分解

Question

我一直在尝试与 CUDA 并行执行许多小矩阵的 QR 分解。 因此，我在 Cublas 中使用了 cublasDgeqrfBatched 函数。我找不到上述功能的工作示例，并且在调用它的文档中发现了一些歧义。 事实上，我尝试在Wikipedia 的 Householder 反射部分的示例中测试cublasDgeqrfBatched，因为cublasDgeqrfBatched 正在使用相同的方法。 2个输入小矩阵是相同的，分别是：

A= 12 -51   4
   6   167 -68
  -4   24  -41

根据文档，Aarray 是一个指向维度为mxn 的矩阵的指针数组，而 TauArray 是一个指向维度至少为max (1, min(m, n). 的向量的指针数组

cublasDgeqrfBatched 对每个Aarray[i] 进行 QR 分解

i =0, ...,batchSize-1

每个矩阵Q[i]存储在每个Aarray[i]的下部

我使用下面的代码来调用这个函数：

#include "cuda_runtime.h"
#include "device_launch_paraMeters.h"

#include<stdlib.h>
#include<stdio.h>
#include<assert.h>
#include <cublas.h>
#include "cublas_v2.h"
#include "Utilities.cuh"
#include <helper_cuda.h>


/********/
/* MAIN */
/********/
int main(){

//mxn: size of Array[i]

const int m = 3;
const int n = 3;
double h_A[3*3*2]={12, -51, 4, 6, 167, -68, -4, 24, -41, 12, -51, 4, 6, 167, -68, -4, 24, -41};// two 3x3 identical matrices for test


const int batchSize=2;//2 small matrices 
const int ltau=3; //ltau = max(1,min(m,n))

// --- CUBLAS initialization
    cublasHandle_t cublas_handle;
    cublasStatus_t stat;
    cublasSafeCall(cublasCreate(&cublas_handle));


// --- CUDA batched QR initialization

    double *d_A, *d_TAU;


checkCudaErrors(cudaMalloc((void**)&d_A, m*n*batchSize*sizeof(double)));  
checkCudaErrors(cudaMalloc((void**)&d_TAU, ltau*batchSize*sizeof(double))); 

checkCudaErrors(cudaMemcpy(d_A,h_A,m*n*batchSize*sizeof(double),cudaMemcpyHostToDevice));

double *d_Aarray[batchSize],*d_TauArray[batchSize];

for (int i = 0; i < batchSize; i++)
{ 
    d_Aarray[i] = d_A+ i*m*n;
    d_TauArray[i] = d_TAU + i*ltau;
}

int lda=3;
int info;

stat=cublasDgeqrfBatched(cublas_handle, m, n, d_Aarray, lda, d_TauArray, &info, batchSize);
if (stat != CUBLAS_STATUS_SUCCESS) 
    printf("\n cublasDgeqrfBatched failed");


double *A0,*A1;
A0=(double*)malloc(m*n*batchSize*sizeof(double));
A1=(double*)malloc(m*n*sizeof(double));

checkCudaErrors(cudaMemcpy(A0,d_Aarray[0],m*n*sizeof(double),cudaMemcpyDeviceToHost));
checkCudaErrors(cudaMemcpy(A1,d_Aarray[1],m*n*sizeof(double),cudaMemcpyDeviceToHost));

}

但是，出现错误"CUDA error batched_QR/kernel.cu:64 code=4(cudaErrorLaunchFailure) "cudaMemcpy(A0,d_Aarray[0],m*n*sizeof(double),cudaMemcpyDeviceToHost)"

我认为指针的使用存在错误，但我无法更正。请问问题出在哪里？

编辑：

为了使 d_Aarray 和 d_TauArray 设备阵列符合建议的要求，我添加了以下内容：

     double *d_A, *d_TAU;


checkCudaErrors(cudaMalloc((void**)&d_A, m*n*batchSize*sizeof(*d_A)));  
checkCudaErrors(cudaMalloc((void**)&d_TAU, ltau*batchSize*sizeof(*d_TAU))); 

checkCudaErrors(cudaMemcpy(d_A,h_A,m*n*batchSize*sizeof(double),cudaMemcpyHostToDevice));
checkCudaErrors(cudaMemset(d_TAU, 0, ltau*batchSize* sizeof(*d_TAU)));

但将结果复制回主机时总是出现相同的错误。

Answer 1

我认为使用指针有错误

你是对的。您传递给cublasDgeqrfBatched 的设备指针数组是主机数组而不是设备数组：

double *d_Aarray[batchSize],*d_TauArray[batchSize];

for (int i = 0; i < batchSize; i++)
{ 
    d_Aarray[i] = d_A+ i*m*n;
    d_TauArray[i] = d_TAU + i*ltau;
}

您必须将d_Aarray 和d_TauArray 复制到设备并将设备副本的地址传递给cublasDgeqrfBatched 才能正常工作。像这样的：

double *d_Aarray[batchSize],*d_TauArray[batchSize];

for (int i = 0; i < batchSize; i++)
{ 
    d_Aarray[i] = d_A+ i*m*n;
    d_TauArray[i] = d_TAU + i*ltau;
}

double ** d_Aarray_, ** d_TauArray_;
cudaMalloc((void **)&d_Aarray_, sizeof(d_Aarray));
cudaMalloc((void **)&d_TauArray_, sizeof(d_TauArray));

cudaMemcpy(d_Aarray_, d_Aarray, sizeof(d_Aarray), cudaMemcpyHostToDevice);
cudaMemcpy(d_TauArray_, d_TauArray, sizeof(d_TauArray), cudaMemcpyHostToDevice);

stat = cublasDgeqrfBatched(cublas_handle, m, n, d_Aarray_, lda, d_TauArray_, &info, batchSize)

[免责声明：用浏览器编写]

这里的d_Aarray_ 和d_TauArray_ 是d_Aarray 和d_TauArray 的设备 内存副本。