https://blog.csdn.net/kewang93/article/details/118702824?spm=1001.2014.3001.5502
这一篇为复数矩阵求逆
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>
#include <math.h>
#include <cuda_runtime.h>
#include <cublas_v2.h>
#include <helper_cuda.h>
//#include "batchCUBLAS.h"
#define T_ELEM cuComplex
#define CUBLASTEST_FAILED 1
#define MATDEM 2
const char *sSDKname = "batchCUBLAS";
static __inline__ int imax(int x, int y)
{
return (x > y) ? x : y;
}
static __inline__ unsigned cuRand(void)
{
/* George Marsaglia's fast inline random number generator */
#define CUDA_ZNEW (cuda_z=36969*(cuda_z&65535)+(cuda_z>>16))
#define CUDA_WNEW (cuda_w=18000*(cuda_w&65535)+(cuda_w>>16))
#define CUDA_MWC ((CUDA_ZNEW<<16)+CUDA_WNEW)
#define CUDA_SHR3 (cuda_jsr=cuda_jsr^(cuda_jsr<<17), \
cuda_jsr = cuda_jsr ^ (cuda_jsr >> 13), \
cuda_jsr = cuda_jsr ^ (cuda_jsr << 5))
#define CUDA_CONG (cuda_jcong=69069*cuda_jcong+1234567)
#define KISS ((CUDA_MWC^CUDA_CONG)+CUDA_SHR3)
static unsigned int cuda_z = 362436069, cuda_w = 521288629;
static unsigned int cuda_jsr = 123456789, cuda_jcong = 380116160;
return KISS;
}
#include <windows.h>
static __inline__ double second(void)
{
LARGE_INTEGER t;
static double oofreq;
static int checkedForHighResTimer;
static BOOL hasHighResTimer;
if (!checkedForHighResTimer)
{
hasHighResTimer = QueryPerformanceFrequency(&t);
oofreq = 1.0 / (double)t.QuadPart;
checkedForHighResTimer = 1;
}
if (hasHighResTimer)
{
QueryPerformanceCounter(&t);
return (double)t.QuadPart * oofreq;
}
else
{
return (double)GetTickCount() / 1000.0;
}
}
struct gemmOpts
{
int m;
int n;
int k;
//testMethod test_method;
char *elem_type;
int N; // number of multiplications
};
struct gemmTestParams
{
cublasOperation_t transa;
cublasOperation_t transb;
int m;
int n;
int k;
T_ELEM alpha;
T_ELEM beta;
};
#define CLEANUP() \
do { \
if (A) free (A); \
if (B) free (B); \
if (C) free (C); \
for(int i = 0; i < opts.N; ++i) { \
if(devPtrA[i]) cudaFree(devPtrA[i]);\
if(devPtrB[i]) cudaFree(devPtrB[i]);\
if(devPtrC[i]) cudaFree(devPtrC[i]);\
} \
if (devPtrA) free(devPtrA); \
if (devPtrB) free(devPtrB); \
if (devPtrC) free(devPtrC); \
if (devPtrA_dev) cudaFree(devPtrA_dev); \
if (devPtrB_dev) cudaFree(devPtrB_dev); \
if (devPtrC_dev) cudaFree(devPtrC_dev); \
fflush (stdout); \
} while (0)
void fillupMatrixDebug(T_ELEM *A, int lda, int rows, int cols)
{
for (int j = 0; j < cols; j++)
{
for (int i = 0; i < rows; i++)
{
A[i + lda*j].x = (i + rand()%13);
A[i + lda*j].y = (j + rand() % 17);
}
}
}
void PrintMatrixDebug(T_ELEM *A, int lda, int rows, int cols)
{
printf("========The matrix is \n");
for (int j = 0; j < cols; j++)
{
for (int i = 0; i < rows; i++)
{
printf("%+.04f%+.04fi ", A[i + lda*j].x, A[i + lda*j].y);
//A[i + lda*j] = (i + j);
}
printf(" \n");
}
}
int main(int argc, char *argv[])
{
printf("%s Starting...\n\n", sSDKname);
int dev = findCudaDevice(argc, (const char **)argv);
if (dev == -1)
{
return CUBLASTEST_FAILED;
}
cublasHandle_t handle;
if (cublasCreate(&handle) != CUBLAS_STATUS_SUCCESS)
{
fprintf(stdout, "CUBLAS initialization failed!\n");
exit(EXIT_FAILURE);
}
cublasStatus_t status1, status2, status3;
T_ELEM *A = NULL;
T_ELEM *B = NULL;
T_ELEM *C = NULL;
T_ELEM **devPtrA = 0;
T_ELEM **devPtrB = 0;
T_ELEM **devPtrC = 0;
T_ELEM **devPtrA_dev = NULL;
T_ELEM **devPtrB_dev = NULL;
T_ELEM **devPtrC_dev = NULL;
int matrixM, matrixN, matrixK;
int rowsA, rowsB, rowsC;
int colsA, colsB, colsC;
int m, n, k;
int matrixSizeA, matrixSizeB, matrixSizeC;
int errors;
double start, stop;
gemmOpts opts;
opts.N = 20;
gemmTestParams params;
printf("Testing Batch INV Cublas\n");
matrixM = MATDEM;
matrixN = MATDEM;
matrixK = MATDEM;
int * info;//用于记录LU分解是否成功
int * pivo;//用于记录LU分解的信息
cudaMalloc((void **)& info, sizeof(int)* opts.N);
cudaMalloc((void **)& pivo, sizeof(int)* matrixM * opts.N);
rowsA = imax(1, matrixM);
colsA = imax(1, matrixK);
rowsB = imax(1, matrixK);
colsB = imax(1, matrixN);
rowsC = imax(1, matrixM);
colsC = imax(1, matrixN);
matrixSizeA = rowsA * colsA;
matrixSizeB = rowsB * colsB;
matrixSizeC = rowsC * colsC;
params.transa = CUBLAS_OP_N;
params.transb = CUBLAS_OP_N;
m = matrixM;
n = matrixN;
k = matrixK;
params.m = m;
params.n = n;
params.k = k;
params.alpha.x = 1.0; params.alpha.y = 0.0;
params.beta.x = 0.0; params.beta.y = 0.0;
devPtrA = (T_ELEM **)malloc(opts.N * sizeof(T_ELEM *));
devPtrB = (T_ELEM **)malloc(opts.N * sizeof(*devPtrB));
devPtrC = (T_ELEM **)malloc(opts.N * sizeof(*devPtrC));
for (int i = 0; i < opts.N; i++)
{
cudaError_t err1 = cudaMalloc((void **)&devPtrA[i], matrixSizeA * sizeof(T_ELEM ));
cudaError_t err2 = cudaMalloc((void **)&devPtrB[i], matrixSizeB * sizeof(devPtrB[0][0]));
cudaError_t err3 = cudaMalloc((void **)&devPtrC[i], matrixSizeC * sizeof(devPtrC[0][0]));
}
// For batched processing we need those arrays on the device
cudaError_t err1 = cudaMalloc((void **)&devPtrA_dev, opts.N * sizeof(*devPtrA));
cudaError_t err2 = cudaMalloc((void **)&devPtrB_dev, opts.N * sizeof(*devPtrB));
cudaError_t err3 = cudaMalloc((void **)&devPtrC_dev, opts.N * sizeof(*devPtrC));
err1 = cudaMemcpy(devPtrA_dev, devPtrA, opts.N * sizeof(*devPtrA), cudaMemcpyHostToDevice);
err2 = cudaMemcpy(devPtrB_dev, devPtrB, opts.N * sizeof(*devPtrB), cudaMemcpyHostToDevice);
err3 = cudaMemcpy(devPtrC_dev, devPtrC, opts.N * sizeof(*devPtrC), cudaMemcpyHostToDevice);
A = (T_ELEM *)malloc(matrixSizeA * sizeof(A[0]));
B = (T_ELEM *)malloc(matrixSizeB * sizeof(B[0]));
C = (T_ELEM *)malloc(matrixSizeC * sizeof(C[0]));
printf("#### args: lda=%d ldb=%d ldc=%d\n", rowsA, rowsB, rowsC);
m = cuRand() % matrixM;
n = cuRand() % matrixN;
k = cuRand() % matrixK;
memset(A, 0xFF, matrixSizeA* sizeof(A[0]));
fillupMatrixDebug(A, rowsA, rowsA, rowsA);
memset(B, 0xFF, matrixSizeB* sizeof(B[0]));
fillupMatrixDebug(B, rowsB, rowsA, rowsA);
for (int i = 0; i < opts.N; i++)
{
status1 = cublasSetMatrix(rowsA, colsA, sizeof(A[0]), A, rowsA, devPtrA[i], rowsA);
status2 = cublasSetMatrix(rowsB, colsB, sizeof(B[0]), B, rowsB, devPtrB[i], rowsB);
status3 = cublasSetMatrix(rowsC, colsC, sizeof(C[0]), C, rowsC, devPtrC[i], rowsC);
}
start = second();
status1 = cublasCgemmBatched(handle, params.transa, params.transb, params.m, params.n,
params.k, ¶ms.alpha, (const T_ELEM **)devPtrA_dev, rowsA,
(const T_ELEM **)devPtrB_dev, rowsB, ¶ms.beta, devPtrC_dev, rowsC, opts.N);
cublasCgetrfBatched(handle, params.m, devPtrA_dev, params.m, pivo, info, opts.N);//第四个参数是矩阵的主导维,由于这里假设数据在内存中的存放是连续的,所以是size
cublasCgetriBatched(handle, params.m, devPtrA_dev, params.m, pivo, devPtrC_dev, params.m, info, opts.N);
//cublasCgetrfBatched(cublasHandle_t handle,
// int n,
// cuComplex *Aarray[],
// int lda,
// int *PivotArray,
// int *infoArray,
// int batchSize);
//cudaMemcpy(devPtrC_dev, devPtrC, opts.N * sizeof(*devPtrC), cudaMemcpyHostToDevice);
for (int i = 0; i < opts.N; i++)
{
//status1 = cublasSetMatrix(rowsA, colsA, sizeof(A[0]), A, rowsA, devPtrA[i], rowsA);
//status2 = cublasSetMatrix(rowsB, colsB, sizeof(B[0]), B, rowsB, devPtrB[i], rowsB);
status3 = cublasGetMatrix(rowsC, colsC, sizeof(C[0]), devPtrC[i], rowsC, C, rowsC);
}
PrintMatrixDebug(A, params.m, params.n, params.k);
PrintMatrixDebug(B, params.m, params.n, params.k);
PrintMatrixDebug(C, params.m, params.n, params.k);
stop = second();
fprintf(stdout, "^^^^ elapsed = %10.8f sec \n", (stop - start));
CLEANUP();
cublasDestroy(handle);
printf("\nTest Summary\n");
system("pause");
return 0;
}