1519f805aSKarl Rupp #if !defined(__CUSPARSEMATIMPL) 29ae82921SPaul Mullowney #define __CUSPARSEMATIMPL 39ae82921SPaul Mullowney 4afb2bd1cSJunchao Zhang #include <petscpkg_version.h> 5303a667bSSatish Balay #include <petsc/private/cudavecimpl.h> 69ae82921SPaul Mullowney 79ae82921SPaul Mullowney #include <cusparse_v2.h> 89ae82921SPaul Mullowney 99ae82921SPaul Mullowney #include <algorithm> 109ae82921SPaul Mullowney #include <vector> 119ae82921SPaul Mullowney 12c41cb2e2SAlejandro Lamas Daviña #include <thrust/device_vector.h> 13c41cb2e2SAlejandro Lamas Daviña #include <thrust/device_ptr.h> 1450ab121bSSatish Balay #include <thrust/device_malloc_allocator.h> 15c41cb2e2SAlejandro Lamas Daviña #include <thrust/transform.h> 16c41cb2e2SAlejandro Lamas Daviña #include <thrust/functional.h> 17554b8892SKarl Rupp #include <thrust/sequence.h> 187eaca502SStefano Zampini #include <thrust/system/system_error.h> 19c41cb2e2SAlejandro Lamas Daviña 20185e5899SJunchao Zhang #if (CUSPARSE_VER_MAJOR > 10 || CUSPARSE_VER_MAJOR == 10 && CUSPARSE_VER_MINOR >= 2) /* According to cuda/10.1.168 on OLCF Summit */ 2157d48284SJunchao Zhang #define CHKERRCUSPARSE(stat)\ 2257d48284SJunchao Zhang do {\ 2357d48284SJunchao Zhang if (PetscUnlikely(stat)) {\ 2457d48284SJunchao Zhang const char *name = cusparseGetErrorName(stat);\ 2557d48284SJunchao Zhang const char *descr = cusparseGetErrorString(stat);\ 26*a952e0d5SBarry Smith if (((stat == CUSPARSE_STATUS_NOT_INITIALIZED) || (stat == CUSPARSE_STATUS_ALLOC_FAILED)) && PetscCUDAInitialized) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_GPU_RESOURCE,"cuSPARSE error %d (%s) : %s. Reports not initialized or alloc failed; this indicates the GPU has run out resources",(int)stat,name,descr); \ 27*a952e0d5SBarry Smith else SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_GPU,"cuSPARSE error %d (%s) : %s",(int)stat,name,descr);\ 2857d48284SJunchao Zhang }\ 2957d48284SJunchao Zhang } while (0) 306d22fe3dSJunchao Zhang #else 31589f383fSBarry Smith #define CHKERRCUSPARSE(stat) do {if (PetscUnlikely(stat)) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_GPU,"cusparse error %d",(int)stat);} while (0) 326d22fe3dSJunchao Zhang #endif 3357d48284SJunchao Zhang 347eaca502SStefano Zampini #define PetscStackCallThrust(body) do { \ 357eaca502SStefano Zampini try { \ 367eaca502SStefano Zampini body; \ 377eaca502SStefano Zampini } catch(thrust::system_error& e) { \ 387eaca502SStefano Zampini SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in Thrust %s",e.what());\ 397eaca502SStefano Zampini } \ 407eaca502SStefano Zampini } while (0) 417eaca502SStefano Zampini 42aa372e3fSPaul Mullowney #if defined(PETSC_USE_COMPLEX) 43aa372e3fSPaul Mullowney #if defined(PETSC_USE_REAL_SINGLE) 44ccdfe979SStefano Zampini const cuComplex PETSC_CUSPARSE_ONE = {1.0f, 0.0f}; 45ccdfe979SStefano Zampini const cuComplex PETSC_CUSPARSE_ZERO = {0.0f, 0.0f}; 46afb2bd1cSJunchao Zhang #elif defined(PETSC_USE_REAL_DOUBLE) 47afb2bd1cSJunchao Zhang const cuDoubleComplex PETSC_CUSPARSE_ONE = {1.0, 0.0}; 48afb2bd1cSJunchao Zhang const cuDoubleComplex PETSC_CUSPARSE_ZERO = {0.0, 0.0}; 49afb2bd1cSJunchao Zhang #endif 50afb2bd1cSJunchao Zhang #else 51afb2bd1cSJunchao Zhang const PetscScalar PETSC_CUSPARSE_ONE = 1.0; 52afb2bd1cSJunchao Zhang const PetscScalar PETSC_CUSPARSE_ZERO = 0.0; 53afb2bd1cSJunchao Zhang #endif 54afb2bd1cSJunchao Zhang 551b0a6780SStefano Zampini #if PETSC_PKG_CUDA_VERSION_GE(9,0,0) 56afb2bd1cSJunchao Zhang #define cusparse_create_analysis_info cusparseCreateCsrsv2Info 57afb2bd1cSJunchao Zhang #define cusparse_destroy_analysis_info cusparseDestroyCsrsv2Info 58afb2bd1cSJunchao Zhang #if defined(PETSC_USE_COMPLEX) 59afb2bd1cSJunchao Zhang #if defined(PETSC_USE_REAL_SINGLE) 60afb2bd1cSJunchao Zhang #define cusparse_get_svbuffsize(a,b,c,d,e,f,g,h,i,j) cusparseCcsrsv2_bufferSize(a,b,c,d,e,(cuComplex*)(f),g,h,i,j) 61afb2bd1cSJunchao Zhang #define cusparse_analysis(a,b,c,d,e,f,g,h,i,j,k) cusparseCcsrsv2_analysis(a,b,c,d,e,(const cuComplex*)(f),g,h,i,j,k) 62afb2bd1cSJunchao Zhang #define cusparse_solve(a,b,c,d,e,f,g,h,i,j,k,l,m,n) cusparseCcsrsv2_solve(a,b,c,d,(const cuComplex*)(e),f,(const cuComplex*)(g),h,i,j,(const cuComplex*)(k),(cuComplex*)(l),m,n) 63afb2bd1cSJunchao Zhang #elif defined(PETSC_USE_REAL_DOUBLE) 64afb2bd1cSJunchao Zhang #define cusparse_get_svbuffsize(a,b,c,d,e,f,g,h,i,j) cusparseZcsrsv2_bufferSize(a,b,c,d,e,(cuDoubleComplex*)(f),g,h,i,j) 65afb2bd1cSJunchao Zhang #define cusparse_analysis(a,b,c,d,e,f,g,h,i,j,k) cusparseZcsrsv2_analysis(a,b,c,d,e,(const cuDoubleComplex*)(f),g,h,i,j,k) 66afb2bd1cSJunchao Zhang #define cusparse_solve(a,b,c,d,e,f,g,h,i,j,k,l,m,n) cusparseZcsrsv2_solve(a,b,c,d,(const cuDoubleComplex*)(e),f,(const cuDoubleComplex*)(g),h,i,j,(const cuDoubleComplex*)(k),(cuDoubleComplex*)(l),m,n) 67afb2bd1cSJunchao Zhang #endif 68afb2bd1cSJunchao Zhang #else /* not complex */ 69afb2bd1cSJunchao Zhang #if defined(PETSC_USE_REAL_SINGLE) 70afb2bd1cSJunchao Zhang #define cusparse_get_svbuffsize cusparseScsrsv2_bufferSize 71afb2bd1cSJunchao Zhang #define cusparse_analysis cusparseScsrsv2_analysis 72afb2bd1cSJunchao Zhang #define cusparse_solve cusparseScsrsv2_solve 73afb2bd1cSJunchao Zhang #elif defined(PETSC_USE_REAL_DOUBLE) 74afb2bd1cSJunchao Zhang #define cusparse_get_svbuffsize cusparseDcsrsv2_bufferSize 75afb2bd1cSJunchao Zhang #define cusparse_analysis cusparseDcsrsv2_analysis 76afb2bd1cSJunchao Zhang #define cusparse_solve cusparseDcsrsv2_solve 77afb2bd1cSJunchao Zhang #endif 78afb2bd1cSJunchao Zhang #endif 79afb2bd1cSJunchao Zhang #else 80afb2bd1cSJunchao Zhang #define cusparse_create_analysis_info cusparseCreateSolveAnalysisInfo 81afb2bd1cSJunchao Zhang #define cusparse_destroy_analysis_info cusparseDestroySolveAnalysisInfo 82afb2bd1cSJunchao Zhang #if defined(PETSC_USE_COMPLEX) 83afb2bd1cSJunchao Zhang #if defined(PETSC_USE_REAL_SINGLE) 84ec42abe4SAlejandro Lamas Daviña #define cusparse_solve(a,b,c,d,e,f,g,h,i,j,k) cusparseCcsrsv_solve((a),(b),(c),(cuComplex*)(d),(e),(cuComplex*)(f),(g),(h),(i),(cuComplex*)(j),(cuComplex*)(k)) 85ec42abe4SAlejandro Lamas Daviña #define cusparse_analysis(a,b,c,d,e,f,g,h,i) cusparseCcsrsv_analysis((a),(b),(c),(d),(e),(cuComplex*)(f),(g),(h),(i)) 861b0a6780SStefano Zampini #elif defined(PETSC_USE_REAL_DOUBLE) 871b0a6780SStefano Zampini #define cusparse_solve(a,b,c,d,e,f,g,h,i,j,k) cusparseZcsrsv_solve((a),(b),(c),(cuDoubleComplex*)(d),(e),(cuDoubleComplex*)(f),(g),(h),(i),(cuDoubleComplex*)(j),(cuDoubleComplex*)(k)) 881b0a6780SStefano Zampini #define cusparse_analysis(a,b,c,d,e,f,g,h,i) cusparseZcsrsv_analysis((a),(b),(c),(d),(e),(cuDoubleComplex*)(f),(g),(h),(i)) 891b0a6780SStefano Zampini #endif 901b0a6780SStefano Zampini #else /* not complex */ 911b0a6780SStefano Zampini #if defined(PETSC_USE_REAL_SINGLE) 921b0a6780SStefano Zampini #define cusparse_solve cusparseScsrsv_solve 931b0a6780SStefano Zampini #define cusparse_analysis cusparseScsrsv_analysis 941b0a6780SStefano Zampini #elif defined(PETSC_USE_REAL_DOUBLE) 951b0a6780SStefano Zampini #define cusparse_solve cusparseDcsrsv_solve 961b0a6780SStefano Zampini #define cusparse_analysis cusparseDcsrsv_analysis 971b0a6780SStefano Zampini #endif 981b0a6780SStefano Zampini #endif 991b0a6780SStefano Zampini #endif 1001b0a6780SStefano Zampini 1011b0a6780SStefano Zampini #if PETSC_PKG_CUDA_VERSION_GE(11,0,0) 1021b0a6780SStefano Zampini #define cusparse_csr2csc cusparseCsr2cscEx2 1031b0a6780SStefano Zampini #if defined(PETSC_USE_COMPLEX) 1041b0a6780SStefano Zampini #if defined(PETSC_USE_REAL_SINGLE) 1051b0a6780SStefano Zampini #define cusparse_scalartype CUDA_C_32F 106039c6fbaSStefano Zampini #define cusparse_csr_spgeam(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t) cusparseCcsrgeam2(a,b,c,(cuComplex*)d,e,f,(cuComplex*)g,h,i,(cuComplex*)j,k,l,(cuComplex*)m,n,o,p,(cuComplex*)q,r,s,t) 107039c6fbaSStefano Zampini #define cusparse_csr_spgeam_bufferSize(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t) cusparseCcsrgeam2_bufferSizeExt(a,b,c,(cuComplex*)d,e,f,(cuComplex*)g,h,i,(cuComplex*)j,k,l,(cuComplex*)m,n,o,p,(cuComplex*)q,r,s,t) 1081b0a6780SStefano Zampini #elif defined(PETSC_USE_REAL_DOUBLE) 1091b0a6780SStefano Zampini #define cusparse_scalartype CUDA_C_64F 110039c6fbaSStefano Zampini #define cusparse_csr_spgeam(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t) cusparseZcsrgeam2(a,b,c,(cuDoubleComplex*)d,e,f,(cuDoubleComplex*)g,h,i,(cuDoubleComplex*)j,k,l,(cuDoubleComplex*)m,n,o,p,(cuDoubleComplex*)q,r,s,t) 111039c6fbaSStefano Zampini #define cusparse_csr_spgeam_bufferSize(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t) cusparseZcsrgeam2_bufferSizeExt(a,b,c,(cuDoubleComplex*)d,e,f,(cuDoubleComplex*)g,h,i,(cuDoubleComplex*)j,k,l,(cuDoubleComplex*)m,n,o,p,(cuDoubleComplex*)q,r,s,t) 1121b0a6780SStefano Zampini #endif 1131b0a6780SStefano Zampini #else /* not complex */ 1141b0a6780SStefano Zampini #if defined(PETSC_USE_REAL_SINGLE) 1151b0a6780SStefano Zampini #define cusparse_scalartype CUDA_R_32F 116039c6fbaSStefano Zampini #define cusparse_csr_spgeam cusparseScsrgeam2 117039c6fbaSStefano Zampini #define cusparse_csr_spgeam_bufferSize cusparseScsrgeam2_bufferSizeExt 1181b0a6780SStefano Zampini #elif defined(PETSC_USE_REAL_DOUBLE) 1191b0a6780SStefano Zampini #define cusparse_scalartype CUDA_R_64F 120039c6fbaSStefano Zampini #define cusparse_csr_spgeam cusparseDcsrgeam2 121039c6fbaSStefano Zampini #define cusparse_csr_spgeam_bufferSize cusparseDcsrgeam2_bufferSizeExt 1221b0a6780SStefano Zampini #endif 1231b0a6780SStefano Zampini #endif 1241b0a6780SStefano Zampini #else 1251b0a6780SStefano Zampini #if defined(PETSC_USE_COMPLEX) 1261b0a6780SStefano Zampini #if defined(PETSC_USE_REAL_SINGLE) 127ec42abe4SAlejandro Lamas Daviña #define cusparse_csr_spmv(a,b,c,d,e,f,g,h,i,j,k,l,m) cusparseCcsrmv((a),(b),(c),(d),(e),(cuComplex*)(f),(g),(cuComplex*)(h),(i),(j),(cuComplex*)(k),(cuComplex*)(l),(cuComplex*)(m)) 128ccdfe979SStefano Zampini #define cusparse_csr_spmm(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p) cusparseCcsrmm((a),(b),(c),(d),(e),(f),(cuComplex*)(g),(h),(cuComplex*)(i),(j),(k),(cuComplex*)(l),(m),(cuComplex*)(n),(cuComplex*)(o),(p)) 129ec42abe4SAlejandro Lamas Daviña #define cusparse_csr2csc(a,b,c,d,e,f,g,h,i,j,k,l) cusparseCcsr2csc((a),(b),(c),(d),(cuComplex*)(e),(f),(g),(cuComplex*)(h),(i),(j),(k),(l)) 130ec42abe4SAlejandro Lamas Daviña #define cusparse_hyb_spmv(a,b,c,d,e,f,g,h) cusparseChybmv((a),(b),(cuComplex*)(c),(d),(e),(cuComplex*)(f),(cuComplex*)(g),(cuComplex*)(h)) 131ec42abe4SAlejandro Lamas Daviña #define cusparse_csr2hyb(a,b,c,d,e,f,g,h,i,j) cusparseCcsr2hyb((a),(b),(c),(d),(cuComplex*)(e),(f),(g),(h),(i),(j)) 132ec42abe4SAlejandro Lamas Daviña #define cusparse_hyb2csr(a,b,c,d,e,f) cusparseChyb2csr((a),(b),(c),(cuComplex*)(d),(e),(f)) 133fcdce8c4SStefano Zampini #define cusparse_csr_spgemm(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t) cusparseCcsrgemm(a,b,c,d,e,f,g,h,(cuComplex*)i,j,k,l,m,(cuComplex*)n,o,p,q,(cuComplex*)r,s,t) 134039c6fbaSStefano Zampini #define cusparse_csr_spgeam(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s) cusparseCcsrgeam(a,b,c,(cuComplex*)d,e,f,(cuComplex*)g,h,i,(cuComplex*)j,k,l,(cuComplex*)m,n,o,p,(cuComplex*)q,r,s) 135aa372e3fSPaul Mullowney #elif defined(PETSC_USE_REAL_DOUBLE) 136ec42abe4SAlejandro Lamas Daviña #define cusparse_csr_spmv(a,b,c,d,e,f,g,h,i,j,k,l,m) cusparseZcsrmv((a),(b),(c),(d),(e),(cuDoubleComplex*)(f),(g),(cuDoubleComplex*)(h),(i),(j),(cuDoubleComplex*)(k),(cuDoubleComplex*)(l),(cuDoubleComplex*)(m)) 137ccdfe979SStefano Zampini #define cusparse_csr_spmm(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p) cusparseZcsrmm((a),(b),(c),(d),(e),(f),(cuDoubleComplex*)(g),(h),(cuDoubleComplex*)(i),(j),(k),(cuDoubleComplex*)(l),(m),(cuDoubleComplex*)(n),(cuDoubleComplex*)(o),(p)) 138ec42abe4SAlejandro Lamas Daviña #define cusparse_csr2csc(a,b,c,d,e,f,g,h,i,j,k,l) cusparseZcsr2csc((a),(b),(c),(d),(cuDoubleComplex*)(e),(f),(g),(cuDoubleComplex*)(h),(i),(j),(k),(l)) 139ec42abe4SAlejandro Lamas Daviña #define cusparse_hyb_spmv(a,b,c,d,e,f,g,h) cusparseZhybmv((a),(b),(cuDoubleComplex*)(c),(d),(e),(cuDoubleComplex*)(f),(cuDoubleComplex*)(g),(cuDoubleComplex*)(h)) 140ec42abe4SAlejandro Lamas Daviña #define cusparse_csr2hyb(a,b,c,d,e,f,g,h,i,j) cusparseZcsr2hyb((a),(b),(c),(d),(cuDoubleComplex*)(e),(f),(g),(h),(i),(j)) 141ec42abe4SAlejandro Lamas Daviña #define cusparse_hyb2csr(a,b,c,d,e,f) cusparseZhyb2csr((a),(b),(c),(cuDoubleComplex*)(d),(e),(f)) 142fcdce8c4SStefano Zampini #define cusparse_csr_spgemm(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t) cusparseZcsrgemm(a,b,c,d,e,f,g,h,(cuDoubleComplex*)i,j,k,l,m,(cuDoubleComplex*)n,o,p,q,(cuDoubleComplex*)r,s,t) 143039c6fbaSStefano Zampini #define cusparse_csr_spgeam(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s) cusparseZcsrgeam(a,b,c,(cuDoubleComplex*)d,e,f,(cuDoubleComplex*)g,h,i,(cuDoubleComplex*)j,k,l,(cuDoubleComplex*)m,n,o,p,(cuDoubleComplex*)q,r,s) 144aa372e3fSPaul Mullowney #endif 145aa372e3fSPaul Mullowney #else 146aa372e3fSPaul Mullowney #if defined(PETSC_USE_REAL_SINGLE) 147aa372e3fSPaul Mullowney #define cusparse_csr_spmv cusparseScsrmv 148ccdfe979SStefano Zampini #define cusparse_csr_spmm cusparseScsrmm 149aa372e3fSPaul Mullowney #define cusparse_csr2csc cusparseScsr2csc 150aa372e3fSPaul Mullowney #define cusparse_hyb_spmv cusparseShybmv 151aa372e3fSPaul Mullowney #define cusparse_csr2hyb cusparseScsr2hyb 152aa372e3fSPaul Mullowney #define cusparse_hyb2csr cusparseShyb2csr 153fcdce8c4SStefano Zampini #define cusparse_csr_spgemm cusparseScsrgemm 154039c6fbaSStefano Zampini #define cusparse_csr_spgeam cusparseScsrgeam 155aa372e3fSPaul Mullowney #elif defined(PETSC_USE_REAL_DOUBLE) 156aa372e3fSPaul Mullowney #define cusparse_csr_spmv cusparseDcsrmv 157ccdfe979SStefano Zampini #define cusparse_csr_spmm cusparseDcsrmm 158aa372e3fSPaul Mullowney #define cusparse_csr2csc cusparseDcsr2csc 159aa372e3fSPaul Mullowney #define cusparse_hyb_spmv cusparseDhybmv 160aa372e3fSPaul Mullowney #define cusparse_csr2hyb cusparseDcsr2hyb 161aa372e3fSPaul Mullowney #define cusparse_hyb2csr cusparseDhyb2csr 162fcdce8c4SStefano Zampini #define cusparse_csr_spgemm cusparseDcsrgemm 163039c6fbaSStefano Zampini #define cusparse_csr_spgeam cusparseDcsrgeam 164aa372e3fSPaul Mullowney #endif 165aa372e3fSPaul Mullowney #endif 166afb2bd1cSJunchao Zhang #endif 167aa372e3fSPaul Mullowney 168aa372e3fSPaul Mullowney #define THRUSTINTARRAY32 thrust::device_vector<int> 169aa372e3fSPaul Mullowney #define THRUSTINTARRAY thrust::device_vector<PetscInt> 170aa372e3fSPaul Mullowney #define THRUSTARRAY thrust::device_vector<PetscScalar> 171aa372e3fSPaul Mullowney 172aa372e3fSPaul Mullowney /* A CSR matrix structure */ 173aa372e3fSPaul Mullowney struct CsrMatrix { 174aa372e3fSPaul Mullowney PetscInt num_rows; 175aa372e3fSPaul Mullowney PetscInt num_cols; 176aa372e3fSPaul Mullowney PetscInt num_entries; 177aa372e3fSPaul Mullowney THRUSTINTARRAY32 *row_offsets; 178aa372e3fSPaul Mullowney THRUSTINTARRAY32 *column_indices; 179aa372e3fSPaul Mullowney THRUSTARRAY *values; 1809ae82921SPaul Mullowney }; 1819ae82921SPaul Mullowney 182aa372e3fSPaul Mullowney /* This is struct holding the relevant data needed to a MatSolve */ 183aa372e3fSPaul Mullowney struct Mat_SeqAIJCUSPARSETriFactorStruct { 184aa372e3fSPaul Mullowney /* Data needed for triangular solve */ 185aa372e3fSPaul Mullowney cusparseMatDescr_t descr; 186aa372e3fSPaul Mullowney cusparseOperation_t solveOp; 187aa372e3fSPaul Mullowney CsrMatrix *csrMat; 1881b0a6780SStefano Zampini #if PETSC_PKG_CUDA_VERSION_GE(9,0,0) 189afb2bd1cSJunchao Zhang csrsv2Info_t solveInfo; 190afb2bd1cSJunchao Zhang #else 191afb2bd1cSJunchao Zhang cusparseSolveAnalysisInfo_t solveInfo; 192afb2bd1cSJunchao Zhang #endif 1931b0a6780SStefano Zampini cusparseSolvePolicy_t solvePolicy; /* whether level information is generated and used */ 1941b0a6780SStefano Zampini int solveBufferSize; 1951b0a6780SStefano Zampini void *solveBuffer; 1961b0a6780SStefano Zampini size_t csr2cscBufferSize; /* to transpose the triangular factor (only used for CUDA >= 11.0) */ 1971b0a6780SStefano Zampini void *csr2cscBuffer; 1982cbc15d9SMark PetscScalar *AA_h; /* managed host buffer for moving values to the GPU */ 199aa372e3fSPaul Mullowney }; 200aa372e3fSPaul Mullowney 201afb2bd1cSJunchao Zhang /* This is a larger struct holding all the triangular factors for a solve, transpose solve, and any indices used in a reordering */ 202aa372e3fSPaul Mullowney struct Mat_SeqAIJCUSPARSETriFactors { 203aa372e3fSPaul Mullowney Mat_SeqAIJCUSPARSETriFactorStruct *loTriFactorPtr; /* pointer for lower triangular (factored matrix) on GPU */ 204aa372e3fSPaul Mullowney Mat_SeqAIJCUSPARSETriFactorStruct *upTriFactorPtr; /* pointer for upper triangular (factored matrix) on GPU */ 205aa372e3fSPaul Mullowney Mat_SeqAIJCUSPARSETriFactorStruct *loTriFactorPtrTranspose; /* pointer for lower triangular (factored matrix) on GPU for the transpose (useful for BiCG) */ 206aa372e3fSPaul Mullowney Mat_SeqAIJCUSPARSETriFactorStruct *upTriFactorPtrTranspose; /* pointer for upper triangular (factored matrix) on GPU for the transpose (useful for BiCG)*/ 207aa372e3fSPaul Mullowney THRUSTINTARRAY *rpermIndices; /* indices used for any reordering */ 208aa372e3fSPaul Mullowney THRUSTINTARRAY *cpermIndices; /* indices used for any reordering */ 209aa372e3fSPaul Mullowney THRUSTARRAY *workVector; 210aa372e3fSPaul Mullowney cusparseHandle_t handle; /* a handle to the cusparse library */ 211aa372e3fSPaul Mullowney PetscInt nnz; /* number of nonzeros ... need this for accurate logging between ICC and ILU */ 212aa372e3fSPaul Mullowney }; 213aa372e3fSPaul Mullowney 214afb2bd1cSJunchao Zhang struct Mat_CusparseSpMV { 215afb2bd1cSJunchao Zhang PetscBool initialized; /* Don't rely on spmvBuffer != NULL to test if the struct is initialized, */ 216afb2bd1cSJunchao Zhang size_t spmvBufferSize; /* since I'm not sure if smvBuffer can be NULL even after cusparseSpMV_bufferSize() */ 217afb2bd1cSJunchao Zhang void *spmvBuffer; 2189db3cbf9SStefano Zampini #if PETSC_PKG_CUDA_VERSION_GE(11,0,0) /* these are present from CUDA 10.1, but PETSc code makes use of them from CUDA 11 on */ 219afb2bd1cSJunchao Zhang cusparseDnVecDescr_t vecXDescr,vecYDescr; /* descriptor for the dense vectors in y=op(A)x */ 2209db3cbf9SStefano Zampini #endif 221afb2bd1cSJunchao Zhang }; 222afb2bd1cSJunchao Zhang 223afb2bd1cSJunchao Zhang /* This is struct holding the relevant data needed to a MatMult */ 224afb2bd1cSJunchao Zhang struct Mat_SeqAIJCUSPARSEMultStruct { 225afb2bd1cSJunchao Zhang void *mat; /* opaque pointer to a matrix. This could be either a cusparseHybMat_t or a CsrMatrix */ 226afb2bd1cSJunchao Zhang cusparseMatDescr_t descr; /* Data needed to describe the matrix for a multiply */ 227afb2bd1cSJunchao Zhang THRUSTINTARRAY *cprowIndices; /* compressed row indices used in the parallel SpMV */ 228afb2bd1cSJunchao Zhang PetscScalar *alpha_one; /* pointer to a device "scalar" storing the alpha parameter in the SpMV */ 229afb2bd1cSJunchao Zhang PetscScalar *beta_zero; /* pointer to a device "scalar" storing the beta parameter in the SpMV as zero*/ 230afb2bd1cSJunchao Zhang PetscScalar *beta_one; /* pointer to a device "scalar" storing the beta parameter in the SpMV as one */ 231afb2bd1cSJunchao Zhang #if PETSC_PKG_CUDA_VERSION_GE(11,0,0) 232afb2bd1cSJunchao Zhang cusparseSpMatDescr_t matDescr; /* descriptor for the matrix, used by SpMV and SpMM */ 233afb2bd1cSJunchao Zhang Mat_CusparseSpMV cuSpMV[3]; /* different Mat_CusparseSpMV structs for non-transpose, transpose, conj-transpose */ 234afb2bd1cSJunchao Zhang Mat_SeqAIJCUSPARSEMultStruct() : matDescr(NULL) { 235afb2bd1cSJunchao Zhang for (int i=0; i<3; i++) cuSpMV[i].initialized = PETSC_FALSE; 236afb2bd1cSJunchao Zhang } 237afb2bd1cSJunchao Zhang #endif 238afb2bd1cSJunchao Zhang }; 239afb2bd1cSJunchao Zhang 240aa372e3fSPaul Mullowney /* This is a larger struct holding all the matrices for a SpMV, and SpMV Tranpose */ 2419ae82921SPaul Mullowney struct Mat_SeqAIJCUSPARSE { 242aa372e3fSPaul Mullowney Mat_SeqAIJCUSPARSEMultStruct *mat; /* pointer to the matrix on the GPU */ 243aa372e3fSPaul Mullowney Mat_SeqAIJCUSPARSEMultStruct *matTranspose; /* pointer to the matrix on the GPU (for the transpose ... useful for BiCG) */ 244aa372e3fSPaul Mullowney THRUSTARRAY *workVector; /* pointer to a workvector to which we can copy the relevant indices of a vector we want to multiply */ 245029808eeSJunchao Zhang THRUSTINTARRAY32 *rowoffsets_gpu; /* rowoffsets on GPU in non-compressed-row format. It is used to convert CSR to CSC */ 246213423ffSJunchao Zhang PetscInt nrows; /* number of rows of the matrix seen by GPU */ 247e057df02SPaul Mullowney MatCUSPARSEStorageFormat format; /* the storage format for the matrix on the device */ 248aa372e3fSPaul Mullowney cudaStream_t stream; /* a stream for the parallel SpMV ... this is not owned and should not be deleted */ 249aa372e3fSPaul Mullowney cusparseHandle_t handle; /* a handle to the cusparse library ... this may not be owned (if we're working in parallel i.e. multiGPUs) */ 25054da937aSStefano Zampini PetscObjectState nonzerostate; /* track nonzero state to possibly recreate the GPU matrix */ 25154da937aSStefano Zampini PetscBool transgen; /* whether or not to generate explicit transpose for MatMultTranspose operations */ 252afb2bd1cSJunchao Zhang #if PETSC_PKG_CUDA_VERSION_GE(11,0,0) 253afb2bd1cSJunchao Zhang size_t csr2cscBufferSize; /* stuff used to compute the matTranspose above */ 254afb2bd1cSJunchao Zhang void *csr2cscBuffer; /* This is used as a C struct and is calloc'ed by PetscNewLog() */ 255afb2bd1cSJunchao Zhang cusparseCsr2CscAlg_t csr2cscAlg; /* algorithms can be selected from command line options */ 256afb2bd1cSJunchao Zhang cusparseSpMVAlg_t spmvAlg; 257afb2bd1cSJunchao Zhang cusparseSpMMAlg_t spmmAlg; 258afb2bd1cSJunchao Zhang #endif 2593fa6b06aSMark Adams PetscSplitCSRDataStructure *deviceMat; /* Matrix on device for, eg, assembly */ 2607e8381f9SStefano Zampini THRUSTINTARRAY *cooPerm; 2617e8381f9SStefano Zampini THRUSTINTARRAY *cooPerm_a; 2629ae82921SPaul Mullowney }; 2639ae82921SPaul Mullowney 2645a576424SJed Brown PETSC_INTERN PetscErrorCode MatCUSPARSECopyToGPU(Mat); 265b06137fdSPaul Mullowney PETSC_INTERN PetscErrorCode MatCUSPARSESetStream(Mat, const cudaStream_t stream); 266b06137fdSPaul Mullowney PETSC_INTERN PetscErrorCode MatCUSPARSESetHandle(Mat, const cusparseHandle_t handle); 267b06137fdSPaul Mullowney PETSC_INTERN PetscErrorCode MatCUSPARSEClearHandle(Mat); 268e61fc153SStefano Zampini PETSC_INTERN PetscErrorCode MatSetPreallocationCOO_SeqAIJCUSPARSE(Mat,PetscInt,const PetscInt[],const PetscInt[]); 269e61fc153SStefano Zampini PETSC_INTERN PetscErrorCode MatSetValuesCOO_SeqAIJCUSPARSE(Mat,const PetscScalar[],InsertMode); 270ed502f03SStefano Zampini PETSC_INTERN PetscErrorCode MatSeqAIJCUSPARSEGetArrayRead(Mat,const PetscScalar**); 271ed502f03SStefano Zampini PETSC_INTERN PetscErrorCode MatSeqAIJCUSPARSERestoreArrayRead(Mat,const PetscScalar**); 272ed502f03SStefano Zampini PETSC_INTERN PetscErrorCode MatSeqAIJCUSPARSEGetArrayWrite(Mat,PetscScalar**); 273ed502f03SStefano Zampini PETSC_INTERN PetscErrorCode MatSeqAIJCUSPARSERestoreArrayWrite(Mat,PetscScalar**); 274039c6fbaSStefano Zampini PETSC_INTERN PetscErrorCode MatSeqAIJCUSPARSEGetArray(Mat,PetscScalar**); 275039c6fbaSStefano Zampini PETSC_INTERN PetscErrorCode MatSeqAIJCUSPARSERestoreArray(Mat,PetscScalar**); 276ed502f03SStefano Zampini PETSC_INTERN PetscErrorCode MatSeqAIJCUSPARSEMergeMats(Mat,Mat,MatReuse,Mat*); 277ed502f03SStefano Zampini 27808391a17SStefano Zampini PETSC_STATIC_INLINE bool isCudaMem(const void *data) 27908391a17SStefano Zampini { 28008391a17SStefano Zampini cudaError_t cerr; 28108391a17SStefano Zampini struct cudaPointerAttributes attr; 28208391a17SStefano Zampini enum cudaMemoryType mtype; 28308391a17SStefano Zampini cerr = cudaPointerGetAttributes(&attr,data); /* Do not check error since before CUDA 11.0, passing a host pointer returns cudaErrorInvalidValue */ 28408391a17SStefano Zampini cudaGetLastError(); /* Reset the last error */ 28508391a17SStefano Zampini #if (CUDART_VERSION < 10000) 28608391a17SStefano Zampini mtype = attr.memoryType; 28708391a17SStefano Zampini #else 28808391a17SStefano Zampini mtype = attr.type; 28908391a17SStefano Zampini #endif 29008391a17SStefano Zampini if (cerr == cudaSuccess && mtype == cudaMemoryTypeDevice) return true; 29108391a17SStefano Zampini else return false; 29208391a17SStefano Zampini } 29308391a17SStefano Zampini 2949ae82921SPaul Mullowney #endif 295