1 #if !defined(CUSPARSEMATIMPL) 2 #define CUSPARSEMATIMPL 3 4 #include <petscpkg_version.h> 5 #include <petsc/private/cudavecimpl.h> 6 #include <petscaijdevice.h> 7 8 #include <cusparse_v2.h> 9 10 #include <algorithm> 11 #include <vector> 12 13 #include <thrust/device_vector.h> 14 #include <thrust/device_ptr.h> 15 #include <thrust/device_malloc_allocator.h> 16 #include <thrust/transform.h> 17 #include <thrust/functional.h> 18 #include <thrust/sequence.h> 19 #include <thrust/system/system_error.h> 20 21 #define PetscStackCallThrust(body) do { \ 22 try { \ 23 body; \ 24 } catch(thrust::system_error& e) { \ 25 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in Thrust %s",e.what());\ 26 } \ 27 } while (0) 28 29 #if defined(PETSC_USE_COMPLEX) 30 #if defined(PETSC_USE_REAL_SINGLE) 31 const cuComplex PETSC_CUSPARSE_ONE = {1.0f, 0.0f}; 32 const cuComplex PETSC_CUSPARSE_ZERO = {0.0f, 0.0f}; 33 #define cusparseXcsrilu02_bufferSize(a,b,c,d,e,f,g,h,i) cusparseCcsrilu02_bufferSize(a,b,c,d,(cuComplex*)e,f,g,h,i) 34 #define cusparseXcsrilu02_analysis(a,b,c,d,e,f,g,h,i,j) cusparseCcsrilu02_analysis(a,b,c,d,(cuComplex*)e,f,g,h,i,j) 35 #define cusparseXcsrilu02(a,b,c,d,e,f,g,h,i,j) cusparseCcsrilu02(a,b,c,d,(cuComplex*)e,f,g,h,i,j) 36 #define cusparseXcsric02_bufferSize(a,b,c,d,e,f,g,h,i) cusparseCcsric02_bufferSize(a,b,c,d,(cuComplex*)e,f,g,h,i) 37 #define cusparseXcsric02_analysis(a,b,c,d,e,f,g,h,i,j) cusparseCcsric02_analysis(a,b,c,d,(cuComplex*)e,f,g,h,i,j) 38 #define cusparseXcsric02(a,b,c,d,e,f,g,h,i,j) cusparseCcsric02(a,b,c,d,(cuComplex*)e,f,g,h,i,j) 39 #elif defined(PETSC_USE_REAL_DOUBLE) 40 const cuDoubleComplex PETSC_CUSPARSE_ONE = {1.0, 0.0}; 41 const cuDoubleComplex PETSC_CUSPARSE_ZERO = {0.0, 0.0}; 42 #define cusparseXcsrilu02_bufferSize(a,b,c,d,e,f,g,h,i) cusparseZcsrilu02_bufferSize(a,b,c,d,(cuDoubleComplex*)e,f,g,h,i) 43 #define cusparseXcsrilu02_analysis(a,b,c,d,e,f,g,h,i,j) cusparseZcsrilu02_analysis(a,b,c,d,(cuDoubleComplex*)e,f,g,h,i,j) 44 #define cusparseXcsrilu02(a,b,c,d,e,f,g,h,i,j) cusparseZcsrilu02(a,b,c,d,(cuDoubleComplex*)e,f,g,h,i,j) 45 #define cusparseXcsric02_bufferSize(a,b,c,d,e,f,g,h,i) cusparseZcsric02_bufferSize(a,b,c,d,(cuDoubleComplex*)e,f,g,h,i) 46 #define cusparseXcsric02_analysis(a,b,c,d,e,f,g,h,i,j) cusparseZcsric02_analysis(a,b,c,d,(cuDoubleComplex*)e,f,g,h,i,j) 47 #define cusparseXcsric02(a,b,c,d,e,f,g,h,i,j) cusparseZcsric02(a,b,c,d,(cuDoubleComplex*)e,f,g,h,i,j) 48 #endif 49 #else 50 const PetscScalar PETSC_CUSPARSE_ONE = 1.0; 51 const PetscScalar PETSC_CUSPARSE_ZERO = 0.0; 52 #if defined(PETSC_USE_REAL_SINGLE) 53 #define cusparseXcsrilu02_bufferSize cusparseScsrilu02_bufferSize 54 #define cusparseXcsrilu02_analysis cusparseScsrilu02_analysis 55 #define cusparseXcsrilu02 cusparseScsrilu02 56 #define cusparseXcsric02_bufferSize cusparseScsric02_bufferSize 57 #define cusparseXcsric02_analysis cusparseScsric02_analysis 58 #define cusparseXcsric02 cusparseScsric02 59 #elif defined(PETSC_USE_REAL_DOUBLE) 60 #define cusparseXcsrilu02_bufferSize cusparseDcsrilu02_bufferSize 61 #define cusparseXcsrilu02_analysis cusparseDcsrilu02_analysis 62 #define cusparseXcsrilu02 cusparseDcsrilu02 63 #define cusparseXcsric02_bufferSize cusparseDcsric02_bufferSize 64 #define cusparseXcsric02_analysis cusparseDcsric02_analysis 65 #define cusparseXcsric02 cusparseDcsric02 66 #endif 67 #endif 68 69 #if PETSC_PKG_CUDA_VERSION_GE(9,0,0) 70 #define csrsvInfo_t csrsv2Info_t 71 #define cusparseCreateCsrsvInfo cusparseCreateCsrsv2Info 72 #define cusparseDestroyCsrsvInfo cusparseDestroyCsrsv2Info 73 #if defined(PETSC_USE_COMPLEX) 74 #if defined(PETSC_USE_REAL_SINGLE) 75 #define cusparseXcsrsv_buffsize(a,b,c,d,e,f,g,h,i,j) cusparseCcsrsv2_bufferSize(a,b,c,d,e,(cuComplex*)(f),g,h,i,j) 76 #define cusparseXcsrsv_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) 77 #define cusparseXcsrsv_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) 78 #elif defined(PETSC_USE_REAL_DOUBLE) 79 #define cusparseXcsrsv_buffsize(a,b,c,d,e,f,g,h,i,j) cusparseZcsrsv2_bufferSize(a,b,c,d,e,(cuDoubleComplex*)(f),g,h,i,j) 80 #define cusparseXcsrsv_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) 81 #define cusparseXcsrsv_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) 82 #endif 83 #else /* not complex */ 84 #if defined(PETSC_USE_REAL_SINGLE) 85 #define cusparseXcsrsv_buffsize cusparseScsrsv2_bufferSize 86 #define cusparseXcsrsv_analysis cusparseScsrsv2_analysis 87 #define cusparseXcsrsv_solve cusparseScsrsv2_solve 88 #elif defined(PETSC_USE_REAL_DOUBLE) 89 #define cusparseXcsrsv_buffsize cusparseDcsrsv2_bufferSize 90 #define cusparseXcsrsv_analysis cusparseDcsrsv2_analysis 91 #define cusparseXcsrsv_solve cusparseDcsrsv2_solve 92 #endif 93 #endif 94 #else 95 #define csrsvInfo_t cusparseSolveAnalysisInfo_t 96 #define cusparseCreateCsrsvInfo cusparseCreateSolveAnalysisInfo 97 #define cusparseDestroyCsrsvInfo cusparseDestroySolveAnalysisInfo 98 #if defined(PETSC_USE_COMPLEX) 99 #if defined(PETSC_USE_REAL_SINGLE) 100 #define cusparseXcsrsv_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)) 101 #define cusparseXcsrsv_analysis(a,b,c,d,e,f,g,h,i) cusparseCcsrsv_analysis((a),(b),(c),(d),(e),(cuComplex*)(f),(g),(h),(i)) 102 #elif defined(PETSC_USE_REAL_DOUBLE) 103 #define cusparseXcsrsv_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)) 104 #define cusparseXcsrsv_analysis(a,b,c,d,e,f,g,h,i) cusparseZcsrsv_analysis((a),(b),(c),(d),(e),(cuDoubleComplex*)(f),(g),(h),(i)) 105 #endif 106 #else /* not complex */ 107 #if defined(PETSC_USE_REAL_SINGLE) 108 #define cusparseXcsrsv_solve cusparseScsrsv_solve 109 #define cusparseXcsrsv_analysis cusparseScsrsv_analysis 110 #elif defined(PETSC_USE_REAL_DOUBLE) 111 #define cusparseXcsrsv_solve cusparseDcsrsv_solve 112 #define cusparseXcsrsv_analysis cusparseDcsrsv_analysis 113 #endif 114 #endif 115 #endif 116 117 #if PETSC_PKG_CUDA_VERSION_GE(11,0,0) 118 #define cusparse_csr2csc cusparseCsr2cscEx2 119 #if defined(PETSC_USE_COMPLEX) 120 #if defined(PETSC_USE_REAL_SINGLE) 121 #define cusparse_scalartype CUDA_C_32F 122 #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) 123 #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) 124 #elif defined(PETSC_USE_REAL_DOUBLE) 125 #define cusparse_scalartype CUDA_C_64F 126 #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) 127 #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) 128 #endif 129 #else /* not complex */ 130 #if defined(PETSC_USE_REAL_SINGLE) 131 #define cusparse_scalartype CUDA_R_32F 132 #define cusparse_csr_spgeam cusparseScsrgeam2 133 #define cusparse_csr_spgeam_bufferSize cusparseScsrgeam2_bufferSizeExt 134 #elif defined(PETSC_USE_REAL_DOUBLE) 135 #define cusparse_scalartype CUDA_R_64F 136 #define cusparse_csr_spgeam cusparseDcsrgeam2 137 #define cusparse_csr_spgeam_bufferSize cusparseDcsrgeam2_bufferSizeExt 138 #endif 139 #endif 140 #else 141 #if defined(PETSC_USE_COMPLEX) 142 #if defined(PETSC_USE_REAL_SINGLE) 143 #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)) 144 #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)) 145 #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)) 146 #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)) 147 #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)) 148 #define cusparse_hyb2csr(a,b,c,d,e,f) cusparseChyb2csr((a),(b),(c),(cuComplex*)(d),(e),(f)) 149 #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) 150 #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) 151 #elif defined(PETSC_USE_REAL_DOUBLE) 152 #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)) 153 #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)) 154 #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)) 155 #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)) 156 #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)) 157 #define cusparse_hyb2csr(a,b,c,d,e,f) cusparseZhyb2csr((a),(b),(c),(cuDoubleComplex*)(d),(e),(f)) 158 #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) 159 #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) 160 #endif 161 #else 162 #if defined(PETSC_USE_REAL_SINGLE) 163 #define cusparse_csr_spmv cusparseScsrmv 164 #define cusparse_csr_spmm cusparseScsrmm 165 #define cusparse_csr2csc cusparseScsr2csc 166 #define cusparse_hyb_spmv cusparseShybmv 167 #define cusparse_csr2hyb cusparseScsr2hyb 168 #define cusparse_hyb2csr cusparseShyb2csr 169 #define cusparse_csr_spgemm cusparseScsrgemm 170 #define cusparse_csr_spgeam cusparseScsrgeam 171 #elif defined(PETSC_USE_REAL_DOUBLE) 172 #define cusparse_csr_spmv cusparseDcsrmv 173 #define cusparse_csr_spmm cusparseDcsrmm 174 #define cusparse_csr2csc cusparseDcsr2csc 175 #define cusparse_hyb_spmv cusparseDhybmv 176 #define cusparse_csr2hyb cusparseDcsr2hyb 177 #define cusparse_hyb2csr cusparseDhyb2csr 178 #define cusparse_csr_spgemm cusparseDcsrgemm 179 #define cusparse_csr_spgeam cusparseDcsrgeam 180 #endif 181 #endif 182 #endif 183 184 #define THRUSTINTARRAY32 thrust::device_vector<int> 185 #define THRUSTINTARRAY thrust::device_vector<PetscInt> 186 #define THRUSTARRAY thrust::device_vector<PetscScalar> 187 188 /* A CSR matrix structure */ 189 struct CsrMatrix { 190 PetscInt num_rows; 191 PetscInt num_cols; 192 PetscInt num_entries; 193 THRUSTINTARRAY32 *row_offsets; 194 THRUSTINTARRAY32 *column_indices; 195 THRUSTARRAY *values; 196 }; 197 198 /* This is struct holding the relevant data needed to a MatSolve */ 199 struct Mat_SeqAIJCUSPARSETriFactorStruct { 200 /* Data needed for triangular solve */ 201 cusparseMatDescr_t descr; 202 cusparseOperation_t solveOp; 203 CsrMatrix *csrMat; 204 csrsvInfo_t solveInfo; 205 cusparseSolvePolicy_t solvePolicy; /* whether level information is generated and used */ 206 int solveBufferSize; 207 void *solveBuffer; 208 size_t csr2cscBufferSize; /* to transpose the triangular factor (only used for CUDA >= 11.0) */ 209 void *csr2cscBuffer; 210 PetscScalar *AA_h; /* managed host buffer for moving values to the GPU */ 211 }; 212 213 /* This is a larger struct holding all the triangular factors for a solve, transpose solve, and any indices used in a reordering */ 214 struct Mat_SeqAIJCUSPARSETriFactors { 215 Mat_SeqAIJCUSPARSETriFactorStruct *loTriFactorPtr; /* pointer for lower triangular (factored matrix) on GPU */ 216 Mat_SeqAIJCUSPARSETriFactorStruct *upTriFactorPtr; /* pointer for upper triangular (factored matrix) on GPU */ 217 Mat_SeqAIJCUSPARSETriFactorStruct *loTriFactorPtrTranspose; /* pointer for lower triangular (factored matrix) on GPU for the transpose (useful for BiCG) */ 218 Mat_SeqAIJCUSPARSETriFactorStruct *upTriFactorPtrTranspose; /* pointer for upper triangular (factored matrix) on GPU for the transpose (useful for BiCG)*/ 219 THRUSTINTARRAY *rpermIndices; /* indices used for any reordering */ 220 THRUSTINTARRAY *cpermIndices; /* indices used for any reordering */ 221 THRUSTARRAY *workVector; 222 cusparseHandle_t handle; /* a handle to the cusparse library */ 223 PetscInt nnz; /* number of nonzeros ... need this for accurate logging between ICC and ILU */ 224 PetscScalar *a_band_d; /* GPU data for banded CSR LU factorization matrix diag(L)=1 */ 225 int *i_band_d; /* this could be optimized away */ 226 cudaDeviceProp dev_prop; 227 PetscBool init_dev_prop; 228 229 /* csrilu0/csric0 appeared in cusparse-8.0, but we use it along with cusparseSpSV, 230 which first appeared in cusparse-11.5 with cuda-11.3. 231 */ 232 PetscBool factorizeOnDevice; /* Do factorization on device or not */ 233 #if CUSPARSE_VERSION >= 11500 234 PetscScalar *csrVal; 235 int *csrRowPtr,*csrColIdx; /* a,i,j of M. Using int since some cusparse APIs only support 32-bit indices */ 236 237 /* Mixed mat descriptor types? yes, different cusparse APIs use different types */ 238 cusparseMatDescr_t matDescr_M; 239 cusparseSpMatDescr_t spMatDescr_L,spMatDescr_U; 240 cusparseSpSVDescr_t spsvDescr_L,spsvDescr_Lt,spsvDescr_U,spsvDescr_Ut; 241 242 cusparseDnVecDescr_t dnVecDescr_X,dnVecDescr_Y; 243 PetscScalar *X,*Y; /* data array of dnVec X and Y */ 244 245 /* Mixed size types? yes, CUDA-11.7.0 declared cusparseDcsrilu02_bufferSizeExt() that returns size_t but did not implement it! */ 246 int factBufferSize_M; /* M ~= LU or LLt */ 247 size_t spsvBufferSize_L,spsvBufferSize_Lt,spsvBufferSize_U,spsvBufferSize_Ut; 248 /* cusparse needs various buffers for factorization and solve of L, U, Lt, or Ut. 249 So save memory, we share the factorization buffer with one of spsvBuffer_L/U. 250 */ 251 void *factBuffer_M,*spsvBuffer_L,*spsvBuffer_U,*spsvBuffer_Lt,*spsvBuffer_Ut; 252 253 csrilu02Info_t ilu0Info_M; 254 csric02Info_t ic0Info_M; 255 int structural_zero,numerical_zero; 256 cusparseSolvePolicy_t policy_M; 257 258 /* In MatSolveTranspose() for ILU0, we use the two flags to do on-demand solve */ 259 PetscBool createdTransposeSpSVDescr; /* Have we created SpSV descriptors for Lt, Ut? */ 260 PetscBool updatedTransposeSpSVAnalysis; /* Have we updated SpSV analysis with the latest L, U values? */ 261 262 PetscLogDouble numericFactFlops; /* Estimated FLOPs in ILU0/ICC0 numeric factorization */ 263 #endif 264 }; 265 266 struct Mat_CusparseSpMV { 267 PetscBool initialized; /* Don't rely on spmvBuffer != NULL to test if the struct is initialized, */ 268 size_t spmvBufferSize; /* since I'm not sure if smvBuffer can be NULL even after cusparseSpMV_bufferSize() */ 269 void *spmvBuffer; 270 #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 */ 271 cusparseDnVecDescr_t vecXDescr,vecYDescr; /* descriptor for the dense vectors in y=op(A)x */ 272 #endif 273 }; 274 275 /* This is struct holding the relevant data needed to a MatMult */ 276 struct Mat_SeqAIJCUSPARSEMultStruct { 277 void *mat; /* opaque pointer to a matrix. This could be either a cusparseHybMat_t or a CsrMatrix */ 278 cusparseMatDescr_t descr; /* Data needed to describe the matrix for a multiply */ 279 THRUSTINTARRAY *cprowIndices; /* compressed row indices used in the parallel SpMV */ 280 PetscScalar *alpha_one; /* pointer to a device "scalar" storing the alpha parameter in the SpMV */ 281 PetscScalar *beta_zero; /* pointer to a device "scalar" storing the beta parameter in the SpMV as zero*/ 282 PetscScalar *beta_one; /* pointer to a device "scalar" storing the beta parameter in the SpMV as one */ 283 #if PETSC_PKG_CUDA_VERSION_GE(11,0,0) 284 cusparseSpMatDescr_t matDescr; /* descriptor for the matrix, used by SpMV and SpMM */ 285 Mat_CusparseSpMV cuSpMV[3]; /* different Mat_CusparseSpMV structs for non-transpose, transpose, conj-transpose */ 286 Mat_SeqAIJCUSPARSEMultStruct() : matDescr(NULL) { 287 for (int i=0; i<3; i++) cuSpMV[i].initialized = PETSC_FALSE; 288 } 289 #endif 290 }; 291 292 /* This is a larger struct holding all the matrices for a SpMV, and SpMV Transpose */ 293 struct Mat_SeqAIJCUSPARSE { 294 Mat_SeqAIJCUSPARSEMultStruct *mat; /* pointer to the matrix on the GPU */ 295 Mat_SeqAIJCUSPARSEMultStruct *matTranspose; /* pointer to the matrix on the GPU (for the transpose ... useful for BiCG) */ 296 THRUSTARRAY *workVector; /* pointer to a workvector to which we can copy the relevant indices of a vector we want to multiply */ 297 THRUSTINTARRAY32 *rowoffsets_gpu; /* rowoffsets on GPU in non-compressed-row format. It is used to convert CSR to CSC */ 298 PetscInt nrows; /* number of rows of the matrix seen by GPU */ 299 MatCUSPARSEStorageFormat format; /* the storage format for the matrix on the device */ 300 PetscBool use_cpu_solve; /* Use AIJ_Seq (I)LU solve */ 301 cudaStream_t stream; /* a stream for the parallel SpMV ... this is not owned and should not be deleted */ 302 cusparseHandle_t handle; /* a handle to the cusparse library ... this may not be owned (if we're working in parallel i.e. multiGPUs) */ 303 PetscObjectState nonzerostate; /* track nonzero state to possibly recreate the GPU matrix */ 304 #if PETSC_PKG_CUDA_VERSION_GE(11,0,0) 305 size_t csr2cscBufferSize; /* stuff used to compute the matTranspose above */ 306 void *csr2cscBuffer; /* This is used as a C struct and is calloc'ed by PetscNewLog() */ 307 cusparseCsr2CscAlg_t csr2cscAlg; /* algorithms can be selected from command line options */ 308 cusparseSpMVAlg_t spmvAlg; 309 cusparseSpMMAlg_t spmmAlg; 310 #endif 311 THRUSTINTARRAY *csr2csc_i; 312 PetscSplitCSRDataStructure deviceMat; /* Matrix on device for, eg, assembly */ 313 314 /* Stuff for basic COO support */ 315 THRUSTINTARRAY *cooPerm; /* permutation array that sorts the input coo entris by row and col */ 316 THRUSTINTARRAY *cooPerm_a; /* ordered array that indicate i-th nonzero (after sorting) is the j-th unique nonzero */ 317 318 /* Stuff for extended COO support */ 319 PetscBool use_extended_coo; /* Use extended COO format */ 320 PetscCount *jmap_d; /* perm[disp+jmap[i]..disp+jmap[i+1]) gives indices of entries in v[] associated with i-th nonzero of the matrix */ 321 PetscCount *perm_d; 322 323 Mat_SeqAIJCUSPARSE() : use_extended_coo(PETSC_FALSE), perm_d(NULL), jmap_d(NULL) {} 324 }; 325 326 PETSC_INTERN PetscErrorCode MatSeqAIJCUSPARSECopyToGPU(Mat); 327 PETSC_INTERN PetscErrorCode MatSetPreallocationCOO_SeqAIJCUSPARSE_Basic(Mat,PetscCount,const PetscInt[],const PetscInt[]); 328 PETSC_INTERN PetscErrorCode MatSetValuesCOO_SeqAIJCUSPARSE_Basic(Mat,const PetscScalar[],InsertMode); 329 PETSC_INTERN PetscErrorCode MatSeqAIJCUSPARSEMergeMats(Mat,Mat,MatReuse,Mat*); 330 PETSC_INTERN PetscErrorCode MatSeqAIJCUSPARSETriFactors_Reset(Mat_SeqAIJCUSPARSETriFactors_p*); 331 332 static inline bool isCudaMem(const void *data) 333 { 334 cudaError_t cerr; 335 struct cudaPointerAttributes attr; 336 enum cudaMemoryType mtype; 337 cerr = cudaPointerGetAttributes(&attr,data); /* Do not check error since before CUDA 11.0, passing a host pointer returns cudaErrorInvalidValue */ 338 cudaGetLastError(); /* Reset the last error */ 339 #if (CUDART_VERSION < 10000) 340 mtype = attr.memoryType; 341 #else 342 mtype = attr.type; 343 #endif 344 if (cerr == cudaSuccess && mtype == cudaMemoryTypeDevice) return true; 345 else return false; 346 } 347 348 #endif 349