xref: /petsc/src/mat/impls/aij/seq/seqcusparse/cusparsematimpl.h (revision a952e0d5b3c7341b1f0ba71fc417c6e83b2e888e)
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