使用 cublasSgetriBatched 在 gpu 上求逆两个矩阵

Posted

技术标签:

【中文标题】使用 cublasSgetriBatched 在 gpu 上求逆两个矩阵【英文标题】:Inversion of two matrices on a gpu using cublasSgetriBatched 【发布时间】:2021-07-26 04:50:46 【问题描述】:

我是 cublas 的新手。我想在 GPU 上并行计算两个矩阵的逆。矩阵是 [4 8;3 9] 和 [5 2;1 7]。是否可以使用 cublasSgetriBatched 来做到这一点?这是我的代码,我得到的结果不正确。这里我取了 2x2 矩阵,但我想找到一种方法来解决多个 mxm 矩阵的这个问题。

#include <stdio.h>

#include "cuda_runtime.h" 
#include "device_launch_parameters.h"

#include "cublas_v2.h"



int main() 

const unsigned int N = 2; 

const unsigned int Nmatrices = 2;

cublasHandle_t handle;
cublasCreate(&handle);

// --- Matrices to be inverted 
float *h_A = new float[N*N*Nmatrices];
float *r_A = new float[N*N*Nmatrices];//result 

h_A[0] = 4.f;  
h_A[1] = 3.f;
h_A[2] = 8.f;
h_A[3] = 9.f;
h_A[4] = 5.f; 
h_A[5] = 1.f; 
h_A[6] = 2.f; 
h_A[7] = 7.f;


// --- Allocate device matrices 
float *d_A; cudaMalloc((void**)&d_A, N*N*Nmatrices*sizeof(float));
float *c_A; cudaMalloc((void**)&c_A, N*N*Nmatrices*sizeof(float));

// --- Move the matrix to be inverted from host to device
cudaMemcpy(d_A,h_A,N*N*Nmatrices*sizeof(float),cudaMemcpyHostToDevice);

// --- Creating the array of pointers needed as input to the batched getrf
float **h_inout_pointers = (float **)malloc(Nmatrices*sizeof(float *));
//for (int i=0; i<Nmatrices; i++) h_inout_pointers[i]=(float *)((char*)d_A+i*((size_t)N*N)*sizeof(float));
*h_inout_pointers=d_A;

float **d_inout_pointers;
cudaMalloc((void**)&d_inout_pointers, Nmatrices*sizeof(float *));
cudaMemcpy(d_inout_pointers,h_inout_pointers,Nmatrices*sizeof(float *),cudaMemcpyHostToDevice);
//free(h_inout_pointers);

float **r_inout_pointers = (float **)malloc(Nmatrices*sizeof(float *));
//for (int i=0; i<Nmatrices; i++) h_inout_pointers[i]=(float *)((char*)d_A+i*((size_t)N*N)*sizeof(float));
*r_inout_pointers=c_A;

float **rd_inout_pointers;
cudaMalloc((void**)&rd_inout_pointers, Nmatrices*sizeof(float *));
cudaMemcpy(rd_inout_pointers,r_inout_pointers,Nmatrices*sizeof(float *),cudaMemcpyHostToDevice);

int *d_PivotArray; cudaMalloc((void**)&d_PivotArray, N*Nmatrices*sizeof(int));
int *d_InfoArray;  cudaMalloc((void**)&d_InfoArray,  Nmatrices*sizeof(int));

int *h_PivotArray = (int *)malloc(N*Nmatrices*sizeof(int));
int *h_InfoArray  = (int *)malloc(  Nmatrices*sizeof(int));

cublasSgetrfBatched(handle, N, d_inout_pointers, N, d_PivotArray, d_InfoArray, Nmatrices);
//cublasSafeCall(cublasSgetrfBatched(handle, N, d_inout_pointers, N, NULL, d_InfoArray, Nmatrices));

//gpuErrchk(cudaMemcpy(h_InfoArray,d_InfoArray,Nmatrices*sizeof(int),cudaMemcpyDeviceToHost));
cublasSgetriBatched(handle, N, d_inout_pointers, N, d_PivotArray, rd_inout_pointers, N, d_InfoArray, 
Nmatrices);


cudaMemcpy(h_A,d_A,N*N*sizeof(float),cudaMemcpyDeviceToHost);
cudaMemcpy(r_A,c_A,N*N*sizeof(float),cudaMemcpyDeviceToHost);
//gpuErrchk(cudaMemcpy(h_PivotArray,d_PivotArray,N*Nmatrices*sizeof(int),cudaMemcpyDeviceToHost));

for (int i=0; i<N*N*Nmatrices; i++) printf("A[%i]=%f\n", i, r_A[i]);


return 0;

【问题讨论】:

我的 cublas 知识已经生疏(12 年前),但要小心一些原始的东西,一些 API(如果不是全部的话)在内存中使用 column-major 数据表示。 【参考方案1】:

您的代码至少有 3 个使用问题:

    您没有正确设置 d_A 的指针数组。 您没有正确设置c_A 的指针数组。 您对r_A 的最终cudaMemcpy 语句仅传输一个矩阵值的数据 (N*N)。

以下代码已修复项目,并且运行时没有运行时错误,产生非零结果。如果您认为这些结果是错误的,您应该说明原因以及您认为正确结果应该是什么。

#include <stdio.h>

#include <cublas_v2.h>



int main() 

const unsigned int N = 2;

const unsigned int Nmatrices = 2;

cublasHandle_t handle;
cublasCreate(&handle);

// --- Matrices to be inverted
float *h_A = new float[N*N*Nmatrices];
float *r_A = new float[N*N*Nmatrices];//result

h_A[0] = 4.f;
h_A[1] = 3.f;
h_A[2] = 8.f;
h_A[3] = 9.f;
h_A[4] = 5.f;
h_A[5] = 1.f;
h_A[6] = 2.f;
h_A[7] = 7.f;


// --- Allocate device matrices
float *d_A; cudaMalloc((void**)&d_A, N*N*Nmatrices*sizeof(float));
float *c_A; cudaMalloc((void**)&c_A, N*N*Nmatrices*sizeof(float));

// --- Move the matrix to be inverted from host to device
cudaMemcpy(d_A,h_A,N*N*Nmatrices*sizeof(float),cudaMemcpyHostToDevice);

// --- Creating the array of pointers needed as input to the batched getrf
float **h_inout_pointers = (float **)malloc(Nmatrices*sizeof(float *));
//for (int i=0; i<Nmatrices; i++) h_inout_pointers[i]=(float *)((char*)d_A+i*((size_t)N*N)*sizeof(float));
h_inout_pointers[0]=d_A;
h_inout_pointers[1]=d_A+N*N;

float **d_inout_pointers;
cudaMalloc((void**)&d_inout_pointers, Nmatrices*sizeof(float *));
cudaMemcpy(d_inout_pointers,h_inout_pointers,Nmatrices*sizeof(float *),cudaMemcpyHostToDevice);
//free(h_inout_pointers);

float **r_inout_pointers = (float **)malloc(Nmatrices*sizeof(float *));
//for (int i=0; i<Nmatrices; i++) h_inout_pointers[i]=(float *)((char*)d_A+i*((size_t)N*N)*sizeof(float));
r_inout_pointers[0]=c_A;
r_inout_pointers[1]=c_A+N*N;

float **rd_inout_pointers;
cudaMalloc((void**)&rd_inout_pointers, Nmatrices*sizeof(float *));
cudaMemcpy(rd_inout_pointers,r_inout_pointers,Nmatrices*sizeof(float *),cudaMemcpyHostToDevice);

int *d_PivotArray; cudaMalloc((void**)&d_PivotArray, N*Nmatrices*sizeof(int));
int *d_InfoArray;  cudaMalloc((void**)&d_InfoArray,  Nmatrices*sizeof(int));

int *h_PivotArray = (int *)malloc(N*Nmatrices*sizeof(int));
int *h_InfoArray  = (int *)malloc(  Nmatrices*sizeof(int));

cublasSgetrfBatched(handle, N, d_inout_pointers, N, d_PivotArray, d_InfoArray, Nmatrices);
//cublasSafeCall(cublasSgetrfBatched(handle, N, d_inout_pointers, N, NULL, d_InfoArray, Nmatrices));

//gpuErrchk(cudaMemcpy(h_InfoArray,d_InfoArray,Nmatrices*sizeof(int),cudaMemcpyDeviceToHost));
cublasSgetriBatched(handle, N, d_inout_pointers, N, d_PivotArray, rd_inout_pointers, N, d_InfoArray,
Nmatrices);


cudaMemcpy(h_A,d_A,N*N*sizeof(float),cudaMemcpyDeviceToHost);
cudaMemcpy(r_A,c_A,Nmatrices*N*N*sizeof(float),cudaMemcpyDeviceToHost);
//gpuErrchk(cudaMemcpy(h_PivotArray,d_PivotArray,N*Nmatrices*sizeof(int),cudaMemcpyDeviceToHost));

for (int i=0; i<N*N*Nmatrices; i++) printf("A[%i]=%f\n", i, r_A[i]);


return 0;

$ nvcc -o t133 t133.cu -lcublas
$ cuda-memcheck ./t133
========= CUDA-MEMCHECK
A[0]=0.750000
A[1]=-0.250000
A[2]=-0.666667
A[3]=0.333333
A[4]=0.212121
A[5]=-0.030303
A[6]=-0.060606
A[7]=0.151515
========= ERROR SUMMARY: 0 errors
$

(FWIW我相信结果是正确的,用计算器测试here)

【讨论】:

以上是关于使用 cublasSgetriBatched 在 gpu 上求逆两个矩阵的主要内容,如果未能解决你的问题,请参考以下文章

在哪里使用 CORBA 以及在哪里使用 SNMP 进行监控?

为啥在使用 unicode 时我不能在 :before :after 内容之后使用空格

在哪里使用 callable 以及在哪里使用 Runnable Interface?

在 Observable RxSwift 中使用 'asPromise()' 可以在 PromiseKit Promise 中使用吗?

可以在 SELECT 查询中使用 IF() 但不能在 UPDATE 中使用

使用 React,在使用 react-transition-group 时,在 StrictMode 中不推荐使用 findDOMNode 作为警告