使用 CUBLAS 进行批量 QR 分解

Batched QR decomposition using CUBLAS

我一直在尝试与 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)));

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

I think there is an error in the use of pointers

你是对的。您传递给 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_Aarrayd_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_Aarrayd_TauArray