1 /*
2 Creates hypre ijmatrix from PETSc matrix
3 */
4
5 #include <petscpkg_version.h>
6 #include <petsc/private/petschypre.h>
7 #include <petscmathypre.h>
8 #include <petsc/private/matimpl.h>
9 #include <petsc/private/deviceimpl.h>
10 #include <../src/mat/impls/hypre/mhypre.h>
11 #include <../src/mat/impls/aij/mpi/mpiaij.h>
12 #include <../src/vec/vec/impls/hypre/vhyp.h>
13 #include <HYPRE.h>
14 #include <HYPRE_utilities.h>
15 #include <_hypre_parcsr_ls.h>
16 #include <_hypre_sstruct_ls.h>
17
18 #if PETSC_PKG_HYPRE_VERSION_LT(2, 18, 0)
19 #define hypre_ParCSRMatrixClone(A, B) hypre_ParCSRMatrixCompleteClone(A)
20 #endif
21
22 #if PETSC_PKG_HYPRE_VERSION_GE(2, 15, 0)
23 #define HYPRE_AssumedPartitionCheck() 1
24 #endif
25
26 static PetscErrorCode MatHYPRE_CreateFromMat(Mat, Mat_HYPRE *);
27 static PetscErrorCode MatHYPRE_IJMatrixPreallocate(Mat, Mat, HYPRE_IJMatrix);
28 static PetscErrorCode MatHYPRE_IJMatrixCopyIJ_MPIAIJ(Mat, HYPRE_IJMatrix);
29 static PetscErrorCode MatHYPRE_IJMatrixCopyIJ_SeqAIJ(Mat, HYPRE_IJMatrix);
30 static PetscErrorCode MatHYPRE_MultKernel_Private(Mat, HYPRE_Complex, Vec, HYPRE_Complex, Vec, PetscBool);
31 static PetscErrorCode MatSetValues_HYPRE(Mat, PetscInt, const PetscInt[], PetscInt, const PetscInt[], const PetscScalar[], InsertMode ins);
32
MatHYPRE_IJMatrixPreallocate(Mat A_d,Mat A_o,HYPRE_IJMatrix ij)33 static PetscErrorCode MatHYPRE_IJMatrixPreallocate(Mat A_d, Mat A_o, HYPRE_IJMatrix ij)
34 {
35 PetscInt i, n_d, n_o;
36 const PetscInt *ia_d, *ia_o;
37 PetscBool done_d = PETSC_FALSE, done_o = PETSC_FALSE;
38 HYPRE_Int *nnz_d = NULL, *nnz_o = NULL;
39
40 PetscFunctionBegin;
41 if (A_d) { /* determine number of nonzero entries in local diagonal part */
42 PetscCall(MatGetRowIJ(A_d, 0, PETSC_FALSE, PETSC_FALSE, &n_d, &ia_d, NULL, &done_d));
43 if (done_d) {
44 PetscCall(PetscMalloc1(n_d, &nnz_d));
45 for (i = 0; i < n_d; i++) nnz_d[i] = (HYPRE_Int)(ia_d[i + 1] - ia_d[i]);
46 }
47 PetscCall(MatRestoreRowIJ(A_d, 0, PETSC_FALSE, PETSC_FALSE, NULL, &ia_d, NULL, &done_d));
48 }
49 if (A_o) { /* determine number of nonzero entries in local off-diagonal part */
50 PetscCall(MatGetRowIJ(A_o, 0, PETSC_FALSE, PETSC_FALSE, &n_o, &ia_o, NULL, &done_o));
51 if (done_o) {
52 PetscCall(PetscMalloc1(n_o, &nnz_o));
53 for (i = 0; i < n_o; i++) nnz_o[i] = (HYPRE_Int)(ia_o[i + 1] - ia_o[i]);
54 }
55 PetscCall(MatRestoreRowIJ(A_o, 0, PETSC_FALSE, PETSC_FALSE, &n_o, &ia_o, NULL, &done_o));
56 }
57 if (done_d) { /* set number of nonzeros in HYPRE IJ matrix */
58 if (!done_o) { /* only diagonal part */
59 PetscCall(PetscCalloc1(n_d, &nnz_o));
60 }
61 #if PETSC_PKG_HYPRE_VERSION_GE(2, 16, 0)
62 { /* If we don't do this, the columns of the matrix will be all zeros! */
63 hypre_AuxParCSRMatrix *aux_matrix;
64 aux_matrix = (hypre_AuxParCSRMatrix *)hypre_IJMatrixTranslator(ij);
65 hypre_AuxParCSRMatrixDestroy(aux_matrix);
66 hypre_IJMatrixTranslator(ij) = NULL;
67 PetscCallHYPRE(HYPRE_IJMatrixSetDiagOffdSizes(ij, nnz_d, nnz_o));
68 /* it seems they partially fixed it in 2.19.0 */
69 #if PETSC_PKG_HYPRE_VERSION_LT(2, 19, 0)
70 aux_matrix = (hypre_AuxParCSRMatrix *)hypre_IJMatrixTranslator(ij);
71 hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 1;
72 #endif
73 }
74 #else
75 PetscCallHYPRE(HYPRE_IJMatrixSetDiagOffdSizes(ij, nnz_d, nnz_o));
76 #endif
77 PetscCall(PetscFree(nnz_d));
78 PetscCall(PetscFree(nnz_o));
79 }
80 PetscFunctionReturn(PETSC_SUCCESS);
81 }
82
MatHYPRE_CreateFromMat(Mat A,Mat_HYPRE * hA)83 static PetscErrorCode MatHYPRE_CreateFromMat(Mat A, Mat_HYPRE *hA)
84 {
85 HYPRE_Int rstart, rend, cstart, cend;
86
87 PetscFunctionBegin;
88 PetscCall(PetscLayoutSetUp(A->rmap));
89 PetscCall(PetscLayoutSetUp(A->cmap));
90 rstart = (HYPRE_Int)A->rmap->rstart;
91 rend = (HYPRE_Int)A->rmap->rend;
92 cstart = (HYPRE_Int)A->cmap->rstart;
93 cend = (HYPRE_Int)A->cmap->rend;
94 PetscCall(PetscHYPREInitialize());
95 if (hA->ij) {
96 if (!hA->inner_free) hypre_IJMatrixObject(hA->ij) = NULL;
97 PetscCallHYPRE(HYPRE_IJMatrixDestroy(hA->ij));
98 }
99 PetscCallHYPRE(HYPRE_IJMatrixCreate(hA->comm, rstart, rend - 1, cstart, cend - 1, &hA->ij));
100 PetscCallHYPRE(HYPRE_IJMatrixSetObjectType(hA->ij, HYPRE_PARCSR));
101 {
102 PetscBool same;
103 Mat A_d, A_o;
104 const PetscInt *colmap;
105 PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIAIJ, &same));
106 if (same) {
107 PetscCall(MatMPIAIJGetSeqAIJ(A, &A_d, &A_o, &colmap));
108 PetscCall(MatHYPRE_IJMatrixPreallocate(A_d, A_o, hA->ij));
109 PetscFunctionReturn(PETSC_SUCCESS);
110 }
111 PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIBAIJ, &same));
112 if (same) {
113 PetscCall(MatMPIBAIJGetSeqBAIJ(A, &A_d, &A_o, &colmap));
114 PetscCall(MatHYPRE_IJMatrixPreallocate(A_d, A_o, hA->ij));
115 PetscFunctionReturn(PETSC_SUCCESS);
116 }
117 PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATSEQAIJ, &same));
118 if (same) {
119 PetscCall(MatHYPRE_IJMatrixPreallocate(A, NULL, hA->ij));
120 PetscFunctionReturn(PETSC_SUCCESS);
121 }
122 PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATSEQBAIJ, &same));
123 if (same) {
124 PetscCall(MatHYPRE_IJMatrixPreallocate(A, NULL, hA->ij));
125 PetscFunctionReturn(PETSC_SUCCESS);
126 }
127 }
128 PetscFunctionReturn(PETSC_SUCCESS);
129 }
130
MatHYPRE_IJMatrixCopyIJ(Mat A,HYPRE_IJMatrix ij)131 static PetscErrorCode MatHYPRE_IJMatrixCopyIJ(Mat A, HYPRE_IJMatrix ij)
132 {
133 PetscBool flg;
134
135 PetscFunctionBegin;
136 #if PETSC_PKG_HYPRE_VERSION_LT(2, 19, 0)
137 PetscCallHYPRE(HYPRE_IJMatrixInitialize(ij));
138 #else
139 PetscCallHYPRE(HYPRE_IJMatrixInitialize_v2(ij, HYPRE_MEMORY_HOST));
140 #endif
141 PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIAIJ, &flg));
142 if (flg) {
143 PetscCall(MatHYPRE_IJMatrixCopyIJ_MPIAIJ(A, ij));
144 PetscFunctionReturn(PETSC_SUCCESS);
145 }
146 PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATSEQAIJ, &flg));
147 if (flg) {
148 PetscCall(MatHYPRE_IJMatrixCopyIJ_SeqAIJ(A, ij));
149 PetscFunctionReturn(PETSC_SUCCESS);
150 }
151 PetscCheck(PETSC_FALSE, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "No support for matrix type %s", ((PetscObject)A)->type_name);
152 PetscFunctionReturn(PETSC_SUCCESS);
153 }
154
MatHYPRE_IJMatrixCopyIJ_SeqAIJ(Mat A,HYPRE_IJMatrix ij)155 static PetscErrorCode MatHYPRE_IJMatrixCopyIJ_SeqAIJ(Mat A, HYPRE_IJMatrix ij)
156 {
157 Mat_SeqAIJ *pdiag = (Mat_SeqAIJ *)A->data;
158 HYPRE_Int type;
159 hypre_ParCSRMatrix *par_matrix;
160 hypre_AuxParCSRMatrix *aux_matrix;
161 hypre_CSRMatrix *hdiag;
162 PetscBool sameint = (PetscBool)(sizeof(PetscInt) == sizeof(HYPRE_Int));
163
164 PetscFunctionBegin;
165 PetscCallHYPRE(HYPRE_IJMatrixGetObjectType(ij, &type));
166 PetscCheck(type == HYPRE_PARCSR, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Only HYPRE_PARCSR is supported");
167 PetscCallHYPRE(HYPRE_IJMatrixGetObject(ij, (void **)&par_matrix));
168 hdiag = hypre_ParCSRMatrixDiag(par_matrix);
169 /*
170 this is the Hack part where we monkey directly with the hypre datastructures
171 */
172 if (sameint) {
173 PetscCall(PetscArraycpy(hdiag->i, pdiag->i, A->rmap->n + 1));
174 PetscCall(PetscArraycpy(hdiag->j, pdiag->j, pdiag->nz));
175 } else {
176 PetscInt i;
177
178 for (i = 0; i < A->rmap->n + 1; i++) hdiag->i[i] = (HYPRE_Int)pdiag->i[i];
179 for (i = 0; i < pdiag->nz; i++) hdiag->j[i] = (HYPRE_Int)pdiag->j[i];
180 }
181
182 aux_matrix = (hypre_AuxParCSRMatrix *)hypre_IJMatrixTranslator(ij);
183 hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 0;
184 PetscFunctionReturn(PETSC_SUCCESS);
185 }
186
MatHYPRE_IJMatrixCopyIJ_MPIAIJ(Mat A,HYPRE_IJMatrix ij)187 static PetscErrorCode MatHYPRE_IJMatrixCopyIJ_MPIAIJ(Mat A, HYPRE_IJMatrix ij)
188 {
189 Mat_MPIAIJ *pA = (Mat_MPIAIJ *)A->data;
190 Mat_SeqAIJ *pdiag, *poffd;
191 PetscInt i, *garray = pA->garray, *jj, cstart, *pjj;
192 HYPRE_Int *hjj, type;
193 hypre_ParCSRMatrix *par_matrix;
194 hypre_AuxParCSRMatrix *aux_matrix;
195 hypre_CSRMatrix *hdiag, *hoffd;
196 PetscBool sameint = (PetscBool)(sizeof(PetscInt) == sizeof(HYPRE_Int));
197
198 PetscFunctionBegin;
199 pdiag = (Mat_SeqAIJ *)pA->A->data;
200 poffd = (Mat_SeqAIJ *)pA->B->data;
201 /* cstart is only valid for square MPIAIJ laid out in the usual way */
202 PetscCall(MatGetOwnershipRange(A, &cstart, NULL));
203
204 PetscCallHYPRE(HYPRE_IJMatrixGetObjectType(ij, &type));
205 PetscCheck(type == HYPRE_PARCSR, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Only HYPRE_PARCSR is supported");
206 PetscCallHYPRE(HYPRE_IJMatrixGetObject(ij, (void **)&par_matrix));
207 hdiag = hypre_ParCSRMatrixDiag(par_matrix);
208 hoffd = hypre_ParCSRMatrixOffd(par_matrix);
209
210 if (sameint) {
211 PetscCall(PetscArraycpy(hdiag->i, pdiag->i, pA->A->rmap->n + 1));
212 } else {
213 for (i = 0; i < pA->A->rmap->n + 1; i++) hdiag->i[i] = (HYPRE_Int)pdiag->i[i];
214 }
215
216 hjj = hdiag->j;
217 pjj = pdiag->j;
218 #if PETSC_PKG_HYPRE_VERSION_GE(2, 16, 0)
219 for (i = 0; i < pdiag->nz; i++) hjj[i] = (HYPRE_Int)pjj[i];
220 #else
221 for (i = 0; i < pdiag->nz; i++) hjj[i] = cstart + pjj[i];
222 #endif
223 if (sameint) {
224 PetscCall(PetscArraycpy(hoffd->i, poffd->i, pA->A->rmap->n + 1));
225 } else {
226 for (i = 0; i < pA->A->rmap->n + 1; i++) hoffd->i[i] = (HYPRE_Int)poffd->i[i];
227 }
228
229 jj = (PetscInt *)hoffd->j;
230 #if PETSC_PKG_HYPRE_VERSION_GE(2, 16, 0)
231 PetscCallHYPRE(hypre_CSRMatrixBigInitialize(hoffd));
232 jj = (PetscInt *)hoffd->big_j;
233 #endif
234 pjj = poffd->j;
235 for (i = 0; i < poffd->nz; i++) jj[i] = garray[pjj[i]];
236
237 aux_matrix = (hypre_AuxParCSRMatrix *)hypre_IJMatrixTranslator(ij);
238 hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 0;
239 PetscFunctionReturn(PETSC_SUCCESS);
240 }
241
MatConvert_HYPRE_IS(Mat A,MatType mtype,MatReuse reuse,Mat * B)242 static PetscErrorCode MatConvert_HYPRE_IS(Mat A, MatType mtype, MatReuse reuse, Mat *B)
243 {
244 Mat_HYPRE *mhA = (Mat_HYPRE *)A->data;
245 Mat lA;
246 ISLocalToGlobalMapping rl2g, cl2g;
247 IS is;
248 hypre_ParCSRMatrix *hA;
249 hypre_CSRMatrix *hdiag, *hoffd;
250 MPI_Comm comm;
251 HYPRE_Complex *hdd, *hod, *aa;
252 PetscScalar *data;
253 HYPRE_BigInt *col_map_offd;
254 HYPRE_Int *hdi, *hdj, *hoi, *hoj;
255 PetscInt *ii, *jj, *iptr, *jptr;
256 PetscInt cum, dr, dc, oc, str, stc, nnz, i, jd, jo, M, N;
257 HYPRE_Int type;
258 MatType lmattype = NULL;
259 PetscBool freeparcsr = PETSC_FALSE;
260
261 PetscFunctionBegin;
262 comm = PetscObjectComm((PetscObject)A);
263 PetscCallHYPRE(HYPRE_IJMatrixGetObjectType(mhA->ij, &type));
264 PetscCheck(type == HYPRE_PARCSR, comm, PETSC_ERR_SUP, "Only HYPRE_PARCSR is supported");
265 PetscCallHYPRE(HYPRE_IJMatrixGetObject(mhA->ij, (void **)&hA));
266 #if defined(PETSC_HAVE_HYPRE_DEVICE)
267 if (HYPRE_MEMORY_DEVICE == hypre_IJMatrixMemoryLocation(mhA->ij)) {
268 /* Support by copying back on the host and copy to GPU
269 Kind of inefficient, but this is the best we can do now */
270 #if defined(HYPRE_USING_HIP)
271 lmattype = MATSEQAIJHIPSPARSE;
272 #elif defined(HYPRE_USING_CUDA)
273 lmattype = MATSEQAIJCUSPARSE;
274 #endif
275 hA = hypre_ParCSRMatrixClone_v2(hA, 1, HYPRE_MEMORY_HOST);
276 freeparcsr = PETSC_TRUE;
277 }
278 #endif
279 M = hypre_ParCSRMatrixGlobalNumRows(hA);
280 N = hypre_ParCSRMatrixGlobalNumCols(hA);
281 str = hypre_ParCSRMatrixFirstRowIndex(hA);
282 stc = hypre_ParCSRMatrixFirstColDiag(hA);
283 hdiag = hypre_ParCSRMatrixDiag(hA);
284 hoffd = hypre_ParCSRMatrixOffd(hA);
285 dr = hypre_CSRMatrixNumRows(hdiag);
286 dc = hypre_CSRMatrixNumCols(hdiag);
287 nnz = hypre_CSRMatrixNumNonzeros(hdiag);
288 hdi = hypre_CSRMatrixI(hdiag);
289 hdj = hypre_CSRMatrixJ(hdiag);
290 hdd = hypre_CSRMatrixData(hdiag);
291 oc = hypre_CSRMatrixNumCols(hoffd);
292 nnz += hypre_CSRMatrixNumNonzeros(hoffd);
293 hoi = hypre_CSRMatrixI(hoffd);
294 hoj = hypre_CSRMatrixJ(hoffd);
295 hod = hypre_CSRMatrixData(hoffd);
296 if (reuse != MAT_REUSE_MATRIX) {
297 PetscInt *aux;
298
299 /* generate l2g maps for rows and cols */
300 PetscCall(ISCreateStride(comm, dr, str, 1, &is));
301 PetscCall(ISLocalToGlobalMappingCreateIS(is, &rl2g));
302 PetscCall(ISDestroy(&is));
303 col_map_offd = hypre_ParCSRMatrixColMapOffd(hA);
304 PetscCall(PetscMalloc1(dc + oc, &aux));
305 for (i = 0; i < dc; i++) aux[i] = i + stc;
306 for (i = 0; i < oc; i++) aux[i + dc] = col_map_offd[i];
307 PetscCall(ISCreateGeneral(comm, dc + oc, aux, PETSC_OWN_POINTER, &is));
308 PetscCall(ISLocalToGlobalMappingCreateIS(is, &cl2g));
309 PetscCall(ISDestroy(&is));
310 /* create MATIS object */
311 PetscCall(MatCreate(comm, B));
312 PetscCall(MatSetSizes(*B, dr, dc, M, N));
313 PetscCall(MatSetType(*B, MATIS));
314 PetscCall(MatSetLocalToGlobalMapping(*B, rl2g, cl2g));
315 PetscCall(ISLocalToGlobalMappingDestroy(&rl2g));
316 PetscCall(ISLocalToGlobalMappingDestroy(&cl2g));
317
318 /* allocate CSR for local matrix */
319 PetscCall(PetscMalloc1(dr + 1, &iptr));
320 PetscCall(PetscMalloc1(nnz, &jptr));
321 PetscCall(PetscMalloc1(nnz, &data));
322 } else {
323 PetscInt nr;
324 PetscBool done;
325 PetscCall(MatISGetLocalMat(*B, &lA));
326 PetscCall(MatGetRowIJ(lA, 0, PETSC_FALSE, PETSC_FALSE, &nr, (const PetscInt **)&iptr, (const PetscInt **)&jptr, &done));
327 PetscCheck(nr == dr, PETSC_COMM_SELF, PETSC_ERR_USER, "Cannot reuse mat: invalid number of rows in local mat! %" PetscInt_FMT " != %" PetscInt_FMT, nr, dr);
328 PetscCheck(iptr[nr] >= nnz, PETSC_COMM_SELF, PETSC_ERR_USER, "Cannot reuse mat: invalid number of nonzeros in local mat! reuse %" PetscInt_FMT " requested %" PetscInt_FMT, iptr[nr], nnz);
329 PetscCall(MatSeqAIJGetArrayWrite(lA, &data));
330 }
331 /* merge local matrices */
332 ii = iptr;
333 jj = jptr;
334 aa = (HYPRE_Complex *)data; /* this cast fixes the clang error when doing the assignments below: implicit conversion from 'HYPRE_Complex' (aka '_Complex double') to 'double' is not permitted in C++ */
335 *ii = *(hdi++) + *(hoi++);
336 for (jd = 0, jo = 0, cum = 0; *ii < nnz; cum++) {
337 PetscScalar *aold = (PetscScalar *)aa;
338 PetscInt *jold = jj, nc = jd + jo;
339 for (; jd < *hdi; jd++) {
340 *jj++ = *hdj++;
341 *aa++ = *hdd++;
342 }
343 for (; jo < *hoi; jo++) {
344 *jj++ = *hoj++ + dc;
345 *aa++ = *hod++;
346 }
347 *(++ii) = *(hdi++) + *(hoi++);
348 PetscCall(PetscSortIntWithScalarArray(jd + jo - nc, jold, aold));
349 }
350 for (; cum < dr; cum++) *(++ii) = nnz;
351 if (reuse != MAT_REUSE_MATRIX) {
352 Mat_SeqAIJ *a;
353
354 PetscCall(MatCreateSeqAIJWithArrays(PETSC_COMM_SELF, dr, dc + oc, iptr, jptr, data, &lA));
355 /* hack SeqAIJ */
356 a = (Mat_SeqAIJ *)lA->data;
357 a->free_a = PETSC_TRUE;
358 a->free_ij = PETSC_TRUE;
359 if (lmattype) PetscCall(MatConvert(lA, lmattype, MAT_INPLACE_MATRIX, &lA));
360 PetscCall(MatISSetLocalMat(*B, lA));
361 PetscCall(MatDestroy(&lA));
362 } else {
363 PetscCall(MatSeqAIJRestoreArrayWrite(lA, &data));
364 }
365 PetscCall(MatAssemblyBegin(*B, MAT_FINAL_ASSEMBLY));
366 PetscCall(MatAssemblyEnd(*B, MAT_FINAL_ASSEMBLY));
367 if (reuse == MAT_INPLACE_MATRIX) PetscCall(MatHeaderReplace(A, B));
368 if (freeparcsr) PetscCallHYPRE(hypre_ParCSRMatrixDestroy(hA));
369 PetscFunctionReturn(PETSC_SUCCESS);
370 }
371
MatHYPRE_DestroyCOOMat(Mat mat)372 static PetscErrorCode MatHYPRE_DestroyCOOMat(Mat mat)
373 {
374 Mat_HYPRE *hA = (Mat_HYPRE *)mat->data;
375
376 PetscFunctionBegin;
377 if (hA->cooMat) { /* If cooMat is present we need to destroy the column indices */
378 PetscCall(MatDestroy(&hA->cooMat));
379 if (hA->cooMatAttached) {
380 hypre_CSRMatrix *csr;
381 hypre_ParCSRMatrix *parcsr;
382 HYPRE_MemoryLocation mem;
383
384 PetscCallHYPRE(HYPRE_IJMatrixGetObject(hA->ij, (void **)&parcsr));
385 csr = hypre_ParCSRMatrixDiag(parcsr);
386 if (csr) {
387 mem = hypre_CSRMatrixMemoryLocation(csr);
388 PetscStackCallExternalVoid("hypre_TFree", hypre_TFree(hypre_CSRMatrixJ(csr), mem));
389 PetscStackCallExternalVoid("hypre_TFree", hypre_TFree(hypre_CSRMatrixBigJ(csr), mem));
390 }
391 csr = hypre_ParCSRMatrixOffd(parcsr);
392 if (csr) {
393 mem = hypre_CSRMatrixMemoryLocation(csr);
394 PetscStackCallExternalVoid("hypre_TFree", hypre_TFree(hypre_CSRMatrixJ(csr), mem));
395 PetscStackCallExternalVoid("hypre_TFree", hypre_TFree(hypre_CSRMatrixBigJ(csr), mem));
396 }
397 }
398 }
399 hA->cooMatAttached = PETSC_FALSE;
400 PetscFunctionReturn(PETSC_SUCCESS);
401 }
402
MatHYPRE_CreateCOOMat(Mat mat)403 static PetscErrorCode MatHYPRE_CreateCOOMat(Mat mat)
404 {
405 MPI_Comm comm;
406 PetscMPIInt size;
407 PetscLayout rmap, cmap;
408 Mat_HYPRE *hmat = (Mat_HYPRE *)mat->data;
409 MatType matType = MATAIJ; /* default type of cooMat */
410
411 PetscFunctionBegin;
412 /* Build an agent matrix cooMat with AIJ format
413 It has the same sparsity pattern as mat, and also shares the data array with mat. We use cooMat to do the COO work.
414 */
415 PetscCall(PetscObjectGetComm((PetscObject)mat, &comm));
416 PetscCallMPI(MPI_Comm_size(comm, &size));
417 PetscCall(PetscLayoutSetUp(mat->rmap));
418 PetscCall(PetscLayoutSetUp(mat->cmap));
419 PetscCall(MatGetLayouts(mat, &rmap, &cmap));
420
421 #if defined(PETSC_HAVE_HYPRE_DEVICE)
422 if (!mat->boundtocpu) { /* mat will be on device, so will cooMat */
423 #if defined(HYPRE_USING_HIP)
424 matType = MATAIJHIPSPARSE;
425 #elif defined(HYPRE_USING_CUDA)
426 matType = MATAIJCUSPARSE;
427 #elif defined(HYPRE_USING_SYCL) && defined(PETSC_HAVE_KOKKOS_KERNELS)
428 matType = MATAIJKOKKOS;
429 #else
430 SETERRQ(comm, PETSC_ERR_SUP, "No HYPRE device available. Suggest re-installing with Kokkos Kernels");
431 #endif
432 }
433 #endif
434
435 /* Do COO preallocation through cooMat */
436 PetscCall(MatHYPRE_DestroyCOOMat(mat));
437 PetscCall(MatCreate(comm, &hmat->cooMat));
438 PetscCall(MatSetType(hmat->cooMat, matType));
439 PetscCall(MatSetLayouts(hmat->cooMat, rmap, cmap));
440
441 /* allocate local matrices if needed */
442 PetscCall(MatMPIAIJSetPreallocation(hmat->cooMat, 0, NULL, 0, NULL));
443 PetscFunctionReturn(PETSC_SUCCESS);
444 }
445
446 /* Attach cooMat data array to hypre matrix.
447 When AIJCUPMSPARSE will support raw device pointers and not THRUSTARRAY
448 we should swap the arrays: i.e., attach hypre matrix array to cooMat
449 This is because hypre should be in charge of handling the memory,
450 cooMat is only a way to reuse PETSc COO code.
451 attaching the memory will then be done at MatSetValuesCOO time and it will dynamically
452 support hypre matrix migrating to host.
453 */
MatHYPRE_AttachCOOMat(Mat mat)454 static PetscErrorCode MatHYPRE_AttachCOOMat(Mat mat)
455 {
456 Mat_HYPRE *hmat = (Mat_HYPRE *)mat->data;
457 hypre_CSRMatrix *diag, *offd;
458 hypre_ParCSRMatrix *parCSR;
459 HYPRE_MemoryLocation hmem = HYPRE_MEMORY_HOST;
460 PetscMemType pmem;
461 Mat A, B;
462 PetscScalar *a;
463 PetscMPIInt size;
464 MPI_Comm comm;
465
466 PetscFunctionBegin;
467 PetscCheck(hmat->cooMat, PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "HYPRE COO delegate matrix has not been created yet");
468 if (hmat->cooMatAttached) PetscFunctionReturn(PETSC_SUCCESS);
469 PetscCheck(hmat->cooMat->preallocated, PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "HYPRE COO delegate matrix is not preallocated");
470 PetscCall(PetscObjectSetName((PetscObject)hmat->cooMat, "_internal_COO_mat_for_hypre"));
471 PetscCall(PetscObjectGetComm((PetscObject)mat, &comm));
472 PetscCallMPI(MPI_Comm_size(comm, &size));
473
474 /* Alias cooMat's data array to IJMatrix's */
475 PetscCallHYPRE(HYPRE_IJMatrixGetObject(hmat->ij, (void **)&parCSR));
476 diag = hypre_ParCSRMatrixDiag(parCSR);
477 offd = hypre_ParCSRMatrixOffd(parCSR);
478
479 A = (size == 1) ? hmat->cooMat : ((Mat_MPIAIJ *)hmat->cooMat->data)->A;
480 B = (size == 1) ? NULL : ((Mat_MPIAIJ *)hmat->cooMat->data)->B;
481
482 PetscCall(PetscObjectSetName((PetscObject)A, "_internal_COO_mat_for_hypre"));
483 hmem = hypre_CSRMatrixMemoryLocation(diag);
484 PetscCall(MatSeqAIJGetCSRAndMemType(A, NULL, NULL, &a, &pmem));
485 PetscAssert((PetscMemTypeHost(pmem) && hmem == HYPRE_MEMORY_HOST) || (PetscMemTypeDevice(pmem) && hmem == HYPRE_MEMORY_DEVICE), comm, PETSC_ERR_PLIB, "PETSc and hypre's memory types mismatch");
486 PetscStackCallExternalVoid("hypre_TFree", hypre_TFree(hypre_CSRMatrixData(diag), hmem));
487 hypre_CSRMatrixData(diag) = (HYPRE_Complex *)a;
488 hypre_CSRMatrixOwnsData(diag) = 0; /* Take ownership of (j,a) away from hypre. As a result, we need to free them on our own */
489
490 if (B) {
491 hmem = hypre_CSRMatrixMemoryLocation(offd);
492 PetscCall(MatSeqAIJGetCSRAndMemType(B, NULL, NULL, &a, &pmem));
493 PetscAssert((PetscMemTypeHost(pmem) && hmem == HYPRE_MEMORY_HOST) || (PetscMemTypeDevice(pmem) && hmem == HYPRE_MEMORY_DEVICE), comm, PETSC_ERR_PLIB, "PETSc and hypre's memory types mismatch");
494 PetscStackCallExternalVoid("hypre_TFree", hypre_TFree(hypre_CSRMatrixData(offd), hmem));
495 hypre_CSRMatrixData(offd) = (HYPRE_Complex *)a;
496 hypre_CSRMatrixOwnsData(offd) = 0;
497 }
498 hmat->cooMatAttached = PETSC_TRUE;
499 PetscFunctionReturn(PETSC_SUCCESS);
500 }
501
502 // Build COO's coordinate list i[], j[] based on CSR's i[], j[] arrays and the number of local rows 'n'
CSRtoCOO_Private(PetscInt n,const PetscInt ii[],const PetscInt jj[],PetscCount * ncoo,PetscInt ** coo_i,PetscInt ** coo_j)503 static PetscErrorCode CSRtoCOO_Private(PetscInt n, const PetscInt ii[], const PetscInt jj[], PetscCount *ncoo, PetscInt **coo_i, PetscInt **coo_j)
504 {
505 PetscInt *cooi, *cooj;
506
507 PetscFunctionBegin;
508 *ncoo = ii[n];
509 PetscCall(PetscMalloc2(*ncoo, &cooi, *ncoo, &cooj));
510 for (PetscInt i = 0; i < n; i++) {
511 for (PetscInt j = ii[i]; j < ii[i + 1]; j++) cooi[j] = i;
512 }
513 PetscCall(PetscArraycpy(cooj, jj, *ncoo));
514 *coo_i = cooi;
515 *coo_j = cooj;
516 PetscFunctionReturn(PETSC_SUCCESS);
517 }
518
519 // Similar to CSRtoCOO_Private, but the CSR's i[], j[] are of type HYPRE_Int
CSRtoCOO_HYPRE_Int_Private(PetscInt n,const HYPRE_Int ii[],const HYPRE_Int jj[],PetscCount * ncoo,PetscInt ** coo_i,PetscInt ** coo_j)520 static PetscErrorCode CSRtoCOO_HYPRE_Int_Private(PetscInt n, const HYPRE_Int ii[], const HYPRE_Int jj[], PetscCount *ncoo, PetscInt **coo_i, PetscInt **coo_j)
521 {
522 PetscInt *cooi, *cooj;
523
524 PetscFunctionBegin;
525 *ncoo = ii[n];
526 PetscCall(PetscMalloc2(*ncoo, &cooi, *ncoo, &cooj));
527 for (PetscInt i = 0; i < n; i++) {
528 for (HYPRE_Int j = ii[i]; j < ii[i + 1]; j++) cooi[j] = i;
529 }
530 for (PetscCount i = 0; i < *ncoo; i++) cooj[i] = jj[i];
531 *coo_i = cooi;
532 *coo_j = cooj;
533 PetscFunctionReturn(PETSC_SUCCESS);
534 }
535
536 // Build a COO data structure for the seqaij matrix, as if the nonzeros are laid out in the same order as in the CSR
MatSeqAIJGetCOO_Private(Mat A,PetscCount * ncoo,PetscInt ** coo_i,PetscInt ** coo_j)537 static PetscErrorCode MatSeqAIJGetCOO_Private(Mat A, PetscCount *ncoo, PetscInt **coo_i, PetscInt **coo_j)
538 {
539 PetscInt n;
540 const PetscInt *ii, *jj;
541 PetscBool done;
542
543 PetscFunctionBegin;
544 PetscCall(MatGetRowIJ(A, 0, PETSC_FALSE, PETSC_FALSE, &n, &ii, &jj, &done));
545 PetscCheck(done, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Failure for MatGetRowIJ");
546 PetscCall(CSRtoCOO_Private(n, ii, jj, ncoo, coo_i, coo_j));
547 PetscCall(MatRestoreRowIJ(A, 0, PETSC_FALSE, PETSC_FALSE, &n, &ii, &jj, &done));
548 PetscCheck(done, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Failure for MatRestoreRowIJ");
549 PetscFunctionReturn(PETSC_SUCCESS);
550 }
551
552 // Build a COO data structure for the hypreCSRMatrix, as if the nonzeros are laid out in the same order as in the hypreCSRMatrix
hypreCSRMatrixGetCOO_Private(hypre_CSRMatrix * A,PetscCount * ncoo,PetscInt ** coo_i,PetscInt ** coo_j)553 static PetscErrorCode hypreCSRMatrixGetCOO_Private(hypre_CSRMatrix *A, PetscCount *ncoo, PetscInt **coo_i, PetscInt **coo_j)
554 {
555 PetscInt n = hypre_CSRMatrixNumRows(A);
556 HYPRE_Int *ii, *jj;
557 HYPRE_MemoryLocation mem = HYPRE_MEMORY_HOST;
558
559 PetscFunctionBegin;
560 #if defined(PETSC_HAVE_HYPRE_DEVICE)
561 mem = hypre_CSRMatrixMemoryLocation(A);
562 if (mem != HYPRE_MEMORY_HOST) {
563 PetscCount nnz = hypre_CSRMatrixNumNonzeros(A);
564 PetscCall(PetscMalloc2(n + 1, &ii, nnz, &jj));
565 hypre_TMemcpy(ii, hypre_CSRMatrixI(A), HYPRE_Int, n + 1, HYPRE_MEMORY_HOST, mem);
566 hypre_TMemcpy(jj, hypre_CSRMatrixJ(A), HYPRE_Int, nnz, HYPRE_MEMORY_HOST, mem);
567 } else {
568 #else
569 {
570 #endif
571 ii = hypre_CSRMatrixI(A);
572 jj = hypre_CSRMatrixJ(A);
573 }
574 PetscCall(CSRtoCOO_HYPRE_Int_Private(n, ii, jj, ncoo, coo_i, coo_j));
575 if (mem != HYPRE_MEMORY_HOST) PetscCall(PetscFree2(ii, jj));
576 PetscFunctionReturn(PETSC_SUCCESS);
577 }
578
579 static PetscErrorCode MatSetValuesCOOFromCSRMatrix_Private(Mat A, hypre_CSRMatrix *H)
580 {
581 PetscBool iscpu = PETSC_TRUE;
582 PetscScalar *a;
583 HYPRE_MemoryLocation mem = HYPRE_MEMORY_HOST;
584
585 PetscFunctionBegin;
586 #if defined(PETSC_HAVE_HYPRE_DEVICE)
587 mem = hypre_CSRMatrixMemoryLocation(H);
588 PetscCall(PetscObjectTypeCompare((PetscObject)A, MATSEQAIJ, &iscpu));
589 #endif
590 if (iscpu && mem != HYPRE_MEMORY_HOST) {
591 PetscCount nnz = hypre_CSRMatrixNumNonzeros(H);
592 PetscCall(PetscMalloc1(nnz, &a));
593 hypre_TMemcpy(a, hypre_CSRMatrixData(H), PetscScalar, nnz, HYPRE_MEMORY_HOST, mem);
594 } else {
595 a = (PetscScalar *)hypre_CSRMatrixData(H);
596 }
597 PetscCall(MatSetValuesCOO(A, a, INSERT_VALUES));
598 if (iscpu && mem != HYPRE_MEMORY_HOST) PetscCall(PetscFree(a));
599 PetscFunctionReturn(PETSC_SUCCESS);
600 }
601
602 PETSC_INTERN PetscErrorCode MatConvert_AIJ_HYPRE(Mat A, MatType type, MatReuse reuse, Mat *B)
603 {
604 MPI_Comm comm = PetscObjectComm((PetscObject)A);
605 Mat M = NULL, dH = NULL, oH = NULL, dA = NULL, oA = NULL;
606 PetscBool ismpiaij, issbaij, isbaij, boundtocpu = PETSC_TRUE;
607 Mat_HYPRE *hA;
608 PetscMemType memtype = PETSC_MEMTYPE_HOST;
609
610 PetscFunctionBegin;
611 if (PetscDefined(HAVE_HYPRE_DEVICE)) {
612 PetscCall(MatGetCurrentMemType(A, &memtype));
613 PetscCall(PetscHYPREInitialize());
614 boundtocpu = PetscMemTypeHost(memtype) ? PETSC_TRUE : PETSC_FALSE;
615 PetscCallHYPRE(HYPRE_SetMemoryLocation(boundtocpu ? HYPRE_MEMORY_HOST : HYPRE_MEMORY_DEVICE));
616 }
617
618 PetscCall(PetscObjectTypeCompareAny((PetscObject)A, &issbaij, MATSEQSBAIJ, MATMPIBAIJ, ""));
619 PetscCall(PetscObjectTypeCompareAny((PetscObject)A, &isbaij, MATSEQBAIJ, MATMPIBAIJ, ""));
620 if (isbaij || issbaij) { /* handle BAIJ and SBAIJ */
621 PetscBool ismpi;
622 MatType newtype;
623
624 PetscCall(PetscObjectTypeCompareAny((PetscObject)A, &ismpi, MATMPISBAIJ, MATMPIBAIJ, ""));
625 newtype = ismpi ? MATMPIAIJ : MATSEQAIJ;
626 if (reuse == MAT_REUSE_MATRIX) {
627 PetscCall(MatConvert(*B, newtype, MAT_INPLACE_MATRIX, B));
628 PetscCall(MatConvert(A, newtype, MAT_REUSE_MATRIX, B));
629 PetscCall(MatConvert(*B, MATHYPRE, MAT_INPLACE_MATRIX, B));
630 } else if (reuse == MAT_INITIAL_MATRIX) {
631 PetscCall(MatConvert(A, newtype, MAT_INITIAL_MATRIX, B));
632 PetscCall(MatConvert(*B, MATHYPRE, MAT_INPLACE_MATRIX, B));
633 } else {
634 PetscCall(MatConvert(A, newtype, MAT_INPLACE_MATRIX, &A));
635 PetscCall(MatConvert(A, MATHYPRE, MAT_INPLACE_MATRIX, &A));
636 }
637 #if defined(PETSC_HAVE_DEVICE)
638 (*B)->boundtocpu = boundtocpu;
639 #endif
640 PetscFunctionReturn(PETSC_SUCCESS);
641 }
642
643 dA = A;
644 PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIAIJ, &ismpiaij));
645 if (ismpiaij) PetscCall(MatMPIAIJGetSeqAIJ(A, &dA, &oA, NULL));
646
647 if (reuse != MAT_REUSE_MATRIX) {
648 PetscCount coo_n;
649 PetscInt *coo_i, *coo_j;
650
651 PetscCall(MatCreate(comm, &M));
652 PetscCall(MatSetType(M, MATHYPRE));
653 PetscCall(MatSetSizes(M, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N));
654 PetscCall(MatSetOption(M, MAT_SORTED_FULL, PETSC_TRUE));
655 PetscCall(MatSetOption(M, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE));
656
657 hA = (Mat_HYPRE *)M->data;
658 PetscCall(MatHYPRE_CreateFromMat(A, hA));
659 PetscCall(MatHYPRE_IJMatrixCopyIJ(A, hA->ij));
660
661 PetscCall(MatHYPRE_CreateCOOMat(M));
662
663 dH = hA->cooMat;
664 PetscCall(PetscObjectBaseTypeCompare((PetscObject)hA->cooMat, MATMPIAIJ, &ismpiaij));
665 if (ismpiaij) PetscCall(MatMPIAIJGetSeqAIJ(hA->cooMat, &dH, &oH, NULL));
666
667 PetscCall(PetscObjectSetName((PetscObject)dH, "_internal_COO_mat_for_hypre"));
668 PetscCall(MatSeqAIJGetCOO_Private(dA, &coo_n, &coo_i, &coo_j));
669 PetscCall(MatSetPreallocationCOO(dH, coo_n, coo_i, coo_j));
670 PetscCall(PetscFree2(coo_i, coo_j));
671 if (oH) {
672 PetscCall(PetscLayoutDestroy(&oH->cmap));
673 PetscCall(PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)oH), oA->cmap->n, oA->cmap->n, 1, &oH->cmap));
674 PetscCall(MatSeqAIJGetCOO_Private(oA, &coo_n, &coo_i, &coo_j));
675 PetscCall(MatSetPreallocationCOO(oH, coo_n, coo_i, coo_j));
676 PetscCall(PetscFree2(coo_i, coo_j));
677 }
678 hA->cooMat->assembled = PETSC_TRUE;
679
680 M->preallocated = PETSC_TRUE;
681 PetscCall(MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY));
682 PetscCall(MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY));
683
684 PetscCall(MatHYPRE_AttachCOOMat(M));
685 if (reuse == MAT_INITIAL_MATRIX) *B = M;
686 } else M = *B;
687
688 hA = (Mat_HYPRE *)M->data;
689 PetscCheck(hA->cooMat, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "HYPRE COO delegate matrix has not been created yet");
690
691 dH = hA->cooMat;
692 PetscCall(PetscObjectBaseTypeCompare((PetscObject)hA->cooMat, MATMPIAIJ, &ismpiaij));
693 if (ismpiaij) PetscCall(MatMPIAIJGetSeqAIJ(hA->cooMat, &dH, &oH, NULL));
694
695 PetscScalar *a;
696 PetscCall(MatSeqAIJGetCSRAndMemType(dA, NULL, NULL, &a, NULL));
697 PetscCall(MatSetValuesCOO(dH, a, INSERT_VALUES));
698 if (oH) {
699 PetscCall(MatSeqAIJGetCSRAndMemType(oA, NULL, NULL, &a, NULL));
700 PetscCall(MatSetValuesCOO(oH, a, INSERT_VALUES));
701 }
702
703 if (reuse == MAT_INPLACE_MATRIX) PetscCall(MatHeaderReplace(A, &M));
704 #if defined(PETSC_HAVE_DEVICE)
705 (*B)->boundtocpu = boundtocpu;
706 #endif
707 PetscFunctionReturn(PETSC_SUCCESS);
708 }
709
710 static PetscErrorCode MatConvert_HYPRE_AIJ(Mat A, MatType mtype, MatReuse reuse, Mat *B)
711 {
712 Mat M, dA = NULL, oA = NULL;
713 hypre_ParCSRMatrix *parcsr;
714 hypre_CSRMatrix *dH, *oH;
715 MPI_Comm comm;
716 PetscBool ismpiaij, isseqaij;
717
718 PetscFunctionBegin;
719 comm = PetscObjectComm((PetscObject)A);
720 if (reuse == MAT_REUSE_MATRIX) {
721 PetscCall(PetscObjectBaseTypeCompare((PetscObject)*B, MATMPIAIJ, &ismpiaij));
722 PetscCall(PetscObjectBaseTypeCompare((PetscObject)*B, MATSEQAIJ, &isseqaij));
723 PetscCheck(ismpiaij || isseqaij, comm, PETSC_ERR_SUP, "Only MATMPIAIJ or MATSEQAIJ base types are supported");
724 }
725 PetscCall(MatHYPREGetParCSR(A, &parcsr));
726 #if defined(PETSC_HAVE_HYPRE_DEVICE)
727 if (HYPRE_MEMORY_DEVICE == hypre_ParCSRMatrixMemoryLocation(parcsr)) {
728 PetscBool isaij;
729
730 PetscCall(PetscStrcmp(mtype, MATAIJ, &isaij));
731 if (isaij) {
732 PetscMPIInt size;
733
734 PetscCallMPI(MPI_Comm_size(comm, &size));
735 #if defined(HYPRE_USING_HIP)
736 mtype = size > 1 ? MATMPIAIJHIPSPARSE : MATSEQAIJHIPSPARSE;
737 #elif defined(HYPRE_USING_CUDA)
738 mtype = size > 1 ? MATMPIAIJCUSPARSE : MATSEQAIJCUSPARSE;
739 #else
740 mtype = size > 1 ? MATMPIAIJ : MATSEQAIJ;
741 #endif
742 }
743 }
744 #endif
745 dH = hypre_ParCSRMatrixDiag(parcsr);
746 oH = hypre_ParCSRMatrixOffd(parcsr);
747 if (reuse != MAT_REUSE_MATRIX) {
748 PetscCount coo_n;
749 PetscInt *coo_i, *coo_j;
750
751 PetscCall(MatCreate(comm, &M));
752 PetscCall(MatSetType(M, mtype));
753 PetscCall(MatSetSizes(M, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N));
754 PetscCall(MatMPIAIJSetPreallocation(M, 0, NULL, 0, NULL));
755
756 dA = M;
757 PetscCall(PetscObjectBaseTypeCompare((PetscObject)M, MATMPIAIJ, &ismpiaij));
758 if (ismpiaij) PetscCall(MatMPIAIJGetSeqAIJ(M, &dA, &oA, NULL));
759
760 PetscCall(hypreCSRMatrixGetCOO_Private(dH, &coo_n, &coo_i, &coo_j));
761 PetscCall(MatSetPreallocationCOO(dA, coo_n, coo_i, coo_j));
762 PetscCall(PetscFree2(coo_i, coo_j));
763 if (ismpiaij) {
764 HYPRE_Int nc = hypre_CSRMatrixNumCols(oH);
765
766 PetscCall(PetscLayoutDestroy(&oA->cmap));
767 PetscCall(PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)oA), nc, nc, 1, &oA->cmap));
768 PetscCall(hypreCSRMatrixGetCOO_Private(oH, &coo_n, &coo_i, &coo_j));
769 PetscCall(MatSetPreallocationCOO(oA, coo_n, coo_i, coo_j));
770 PetscCall(PetscFree2(coo_i, coo_j));
771
772 /* garray */
773 Mat_MPIAIJ *aij = (Mat_MPIAIJ *)M->data;
774 HYPRE_BigInt *harray = hypre_ParCSRMatrixColMapOffd(parcsr);
775 PetscInt *garray;
776
777 PetscCall(PetscFree(aij->garray));
778 PetscCall(PetscMalloc1(nc, &garray));
779 for (HYPRE_Int i = 0; i < nc; i++) garray[i] = (PetscInt)harray[i];
780 aij->garray = garray;
781 PetscCall(MatSetUpMultiply_MPIAIJ(M));
782 }
783 if (reuse == MAT_INITIAL_MATRIX) *B = M;
784 } else M = *B;
785
786 dA = M;
787 PetscCall(PetscObjectBaseTypeCompare((PetscObject)M, MATMPIAIJ, &ismpiaij));
788 if (ismpiaij) PetscCall(MatMPIAIJGetSeqAIJ(M, &dA, &oA, NULL));
789 PetscCall(MatSetValuesCOOFromCSRMatrix_Private(dA, dH));
790 if (oA) PetscCall(MatSetValuesCOOFromCSRMatrix_Private(oA, oH));
791 M->assembled = PETSC_TRUE;
792 if (reuse == MAT_INPLACE_MATRIX) PetscCall(MatHeaderReplace(A, &M));
793 PetscFunctionReturn(PETSC_SUCCESS);
794 }
795
796 static PetscErrorCode MatAIJGetParCSR_Private(Mat A, hypre_ParCSRMatrix **hA)
797 {
798 hypre_ParCSRMatrix *tA;
799 hypre_CSRMatrix *hdiag, *hoffd;
800 Mat_SeqAIJ *diag, *offd;
801 PetscInt *garray, i, noffd, dnnz, onnz, *row_starts, *col_starts;
802 MPI_Comm comm = PetscObjectComm((PetscObject)A);
803 PetscBool ismpiaij, isseqaij;
804 PetscBool sameint = (PetscBool)(sizeof(PetscInt) == sizeof(HYPRE_Int));
805 HYPRE_Int *hdi = NULL, *hdj = NULL, *hoi = NULL, *hoj = NULL;
806 PetscInt *pdi = NULL, *pdj = NULL, *poi = NULL, *poj = NULL;
807 PetscBool iscuda, iship;
808 #if defined(PETSC_HAVE_DEVICE) && defined(PETSC_HAVE_HYPRE_DEVICE)
809 PetscBool boundtocpu = A->boundtocpu;
810 #else
811 PetscBool boundtocpu = PETSC_TRUE;
812 #endif
813
814 PetscFunctionBegin;
815 PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIAIJ, &ismpiaij));
816 PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATSEQAIJ, &isseqaij));
817 PetscCheck(ismpiaij || isseqaij, comm, PETSC_ERR_SUP, "Unsupported type %s", ((PetscObject)A)->type_name);
818 PetscCall(PetscObjectTypeCompareAny((PetscObject)A, &iscuda, MATSEQAIJCUSPARSE, MATMPIAIJCUSPARSE, ""));
819 PetscCall(PetscObjectTypeCompareAny((PetscObject)A, &iship, MATSEQAIJHIPSPARSE, MATMPIAIJHIPSPARSE, ""));
820 PetscCall(PetscHYPREInitialize());
821 if (ismpiaij) {
822 Mat_MPIAIJ *a = (Mat_MPIAIJ *)A->data;
823
824 diag = (Mat_SeqAIJ *)a->A->data;
825 offd = (Mat_SeqAIJ *)a->B->data;
826 if (!boundtocpu && (iscuda || iship)) {
827 #if defined(HYPRE_USING_CUDA) && defined(PETSC_HAVE_CUDA)
828 if (iscuda) {
829 sameint = PETSC_TRUE;
830 PetscCall(MatSeqAIJCUSPARSEGetIJ(a->A, PETSC_FALSE, (const HYPRE_Int **)&hdi, (const HYPRE_Int **)&hdj));
831 PetscCall(MatSeqAIJCUSPARSEGetIJ(a->B, PETSC_FALSE, (const HYPRE_Int **)&hoi, (const HYPRE_Int **)&hoj));
832 }
833 #endif
834 #if defined(HYPRE_USING_HIP) && defined(PETSC_HAVE_HIP)
835 if (iship) {
836 sameint = PETSC_TRUE;
837 PetscCall(MatSeqAIJHIPSPARSEGetIJ(a->A, PETSC_FALSE, (const HYPRE_Int **)&hdi, (const HYPRE_Int **)&hdj));
838 PetscCall(MatSeqAIJHIPSPARSEGetIJ(a->B, PETSC_FALSE, (const HYPRE_Int **)&hoi, (const HYPRE_Int **)&hoj));
839 }
840 #endif
841 } else {
842 boundtocpu = PETSC_TRUE;
843 pdi = diag->i;
844 pdj = diag->j;
845 poi = offd->i;
846 poj = offd->j;
847 if (sameint) {
848 hdi = (HYPRE_Int *)pdi;
849 hdj = (HYPRE_Int *)pdj;
850 hoi = (HYPRE_Int *)poi;
851 hoj = (HYPRE_Int *)poj;
852 }
853 }
854 garray = a->garray;
855 noffd = a->B->cmap->N;
856 dnnz = diag->nz;
857 onnz = offd->nz;
858 } else {
859 diag = (Mat_SeqAIJ *)A->data;
860 offd = NULL;
861 if (!boundtocpu && (iscuda || iship)) {
862 #if defined(HYPRE_USING_CUDA) && defined(PETSC_HAVE_CUDA)
863 if (iscuda) {
864 sameint = PETSC_TRUE;
865 PetscCall(MatSeqAIJCUSPARSEGetIJ(A, PETSC_FALSE, (const HYPRE_Int **)&hdi, (const HYPRE_Int **)&hdj));
866 }
867 #endif
868 #if defined(HYPRE_USING_HIP) && defined(PETSC_HAVE_HIP)
869 if (iship) {
870 sameint = PETSC_TRUE;
871 PetscCall(MatSeqAIJHIPSPARSEGetIJ(A, PETSC_FALSE, (const HYPRE_Int **)&hdi, (const HYPRE_Int **)&hdj));
872 }
873 #endif
874 } else {
875 boundtocpu = PETSC_TRUE;
876 pdi = diag->i;
877 pdj = diag->j;
878 if (sameint) {
879 hdi = (HYPRE_Int *)pdi;
880 hdj = (HYPRE_Int *)pdj;
881 }
882 }
883 garray = NULL;
884 noffd = 0;
885 dnnz = diag->nz;
886 onnz = 0;
887 }
888
889 /* create a temporary ParCSR */
890 if (HYPRE_AssumedPartitionCheck()) {
891 PetscMPIInt myid;
892
893 PetscCallMPI(MPI_Comm_rank(comm, &myid));
894 row_starts = A->rmap->range + myid;
895 col_starts = A->cmap->range + myid;
896 } else {
897 row_starts = A->rmap->range;
898 col_starts = A->cmap->range;
899 }
900 tA = hypre_ParCSRMatrixCreate(comm, A->rmap->N, A->cmap->N, (HYPRE_BigInt *)row_starts, (HYPRE_BigInt *)col_starts, (HYPRE_Int)noffd, (HYPRE_Int)dnnz, (HYPRE_Int)onnz);
901 #if defined(hypre_ParCSRMatrixOwnsRowStarts)
902 hypre_ParCSRMatrixSetRowStartsOwner(tA, 0);
903 hypre_ParCSRMatrixSetColStartsOwner(tA, 0);
904 #endif
905
906 /* set diagonal part */
907 hdiag = hypre_ParCSRMatrixDiag(tA);
908 if (!sameint) { /* malloc CSR pointers */
909 PetscCall(PetscMalloc2(A->rmap->n + 1, &hdi, dnnz, &hdj));
910 for (i = 0; i < A->rmap->n + 1; i++) hdi[i] = (HYPRE_Int)pdi[i];
911 for (i = 0; i < dnnz; i++) hdj[i] = (HYPRE_Int)pdj[i];
912 }
913 hypre_CSRMatrixI(hdiag) = hdi;
914 hypre_CSRMatrixJ(hdiag) = hdj;
915 hypre_CSRMatrixData(hdiag) = (HYPRE_Complex *)diag->a;
916 hypre_CSRMatrixNumNonzeros(hdiag) = (HYPRE_Int)diag->nz;
917 hypre_CSRMatrixSetDataOwner(hdiag, 0);
918
919 /* set off-diagonal part */
920 hoffd = hypre_ParCSRMatrixOffd(tA);
921 if (offd) {
922 if (!sameint) { /* malloc CSR pointers */
923 PetscCall(PetscMalloc2(A->rmap->n + 1, &hoi, onnz, &hoj));
924 for (i = 0; i < A->rmap->n + 1; i++) hoi[i] = (HYPRE_Int)poi[i];
925 for (i = 0; i < onnz; i++) hoj[i] = (HYPRE_Int)poj[i];
926 }
927 hypre_CSRMatrixI(hoffd) = hoi;
928 hypre_CSRMatrixJ(hoffd) = hoj;
929 hypre_CSRMatrixData(hoffd) = (HYPRE_Complex *)offd->a;
930 hypre_CSRMatrixNumNonzeros(hoffd) = (HYPRE_Int)offd->nz;
931 hypre_CSRMatrixSetDataOwner(hoffd, 0);
932 }
933 #if defined(PETSC_HAVE_HYPRE_DEVICE)
934 PetscCallHYPRE(hypre_ParCSRMatrixInitialize_v2(tA, !boundtocpu ? HYPRE_MEMORY_DEVICE : HYPRE_MEMORY_HOST));
935 #else
936 #if PETSC_PKG_HYPRE_VERSION_LT(2, 18, 0)
937 PetscCallHYPRE(hypre_ParCSRMatrixInitialize(tA));
938 #else
939 PetscCallHYPRE(hypre_ParCSRMatrixInitialize_v2(tA, HYPRE_MEMORY_HOST));
940 #endif
941 #endif
942
943 /* MatrixSetRownnz comes after MatrixInitialize, so the first uses the right memory location */
944 hypre_CSRMatrixSetRownnz(hdiag);
945 if (offd) hypre_CSRMatrixSetRownnz(hoffd);
946
947 hypre_TFree(hypre_ParCSRMatrixColMapOffd(tA), HYPRE_MEMORY_HOST);
948 hypre_ParCSRMatrixSetNumNonzeros(tA);
949 hypre_ParCSRMatrixColMapOffd(tA) = (HYPRE_BigInt *)garray;
950 if (!hypre_ParCSRMatrixCommPkg(tA)) PetscCallHYPRE(hypre_MatvecCommPkgCreate(tA));
951 *hA = tA;
952 PetscFunctionReturn(PETSC_SUCCESS);
953 }
954
955 static PetscErrorCode MatAIJRestoreParCSR_Private(Mat A, hypre_ParCSRMatrix **hA)
956 {
957 hypre_CSRMatrix *hdiag, *hoffd;
958 PetscBool ismpiaij, sameint = (PetscBool)(sizeof(PetscInt) == sizeof(HYPRE_Int));
959 PetscBool iscuda, iship;
960
961 PetscFunctionBegin;
962 PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIAIJ, &ismpiaij));
963 PetscCall(PetscObjectTypeCompareAny((PetscObject)A, &iscuda, MATSEQAIJCUSPARSE, MATMPIAIJCUSPARSE, ""));
964 PetscCall(PetscObjectTypeCompareAny((PetscObject)A, &iship, MATSEQAIJHIPSPARSE, MATMPIAIJHIPSPARSE, ""));
965 #if defined(HYPRE_USING_CUDA) && defined(PETSC_HAVE_CUDA)
966 if (iscuda) sameint = PETSC_TRUE;
967 #elif defined(HYPRE_USING_HIP) && defined(PETSC_HAVE_HIP)
968 if (iship) sameint = PETSC_TRUE;
969 #endif
970 hdiag = hypre_ParCSRMatrixDiag(*hA);
971 hoffd = hypre_ParCSRMatrixOffd(*hA);
972 /* free temporary memory allocated by PETSc
973 set pointers to NULL before destroying tA */
974 if (!sameint) {
975 HYPRE_Int *hi, *hj;
976
977 hi = hypre_CSRMatrixI(hdiag);
978 hj = hypre_CSRMatrixJ(hdiag);
979 PetscCall(PetscFree2(hi, hj));
980 if (ismpiaij) {
981 hi = hypre_CSRMatrixI(hoffd);
982 hj = hypre_CSRMatrixJ(hoffd);
983 PetscCall(PetscFree2(hi, hj));
984 }
985 }
986 hypre_CSRMatrixI(hdiag) = NULL;
987 hypre_CSRMatrixJ(hdiag) = NULL;
988 hypre_CSRMatrixData(hdiag) = NULL;
989 if (ismpiaij) {
990 hypre_CSRMatrixI(hoffd) = NULL;
991 hypre_CSRMatrixJ(hoffd) = NULL;
992 hypre_CSRMatrixData(hoffd) = NULL;
993 }
994 hypre_ParCSRMatrixColMapOffd(*hA) = NULL;
995 hypre_ParCSRMatrixDestroy(*hA);
996 *hA = NULL;
997 PetscFunctionReturn(PETSC_SUCCESS);
998 }
999
1000 /* calls RAP from BoomerAMG:
1001 the resulting ParCSR will not own the column and row starts
1002 It looks like we don't need to have the diagonal entries ordered first */
1003 static PetscErrorCode MatHYPRE_ParCSR_RAP(hypre_ParCSRMatrix *hR, hypre_ParCSRMatrix *hA, hypre_ParCSRMatrix *hP, hypre_ParCSRMatrix **hRAP)
1004 {
1005 #if defined(hypre_ParCSRMatrixOwnsRowStarts)
1006 HYPRE_Int P_owns_col_starts, R_owns_row_starts;
1007 #endif
1008
1009 PetscFunctionBegin;
1010 #if defined(hypre_ParCSRMatrixOwnsRowStarts)
1011 P_owns_col_starts = hypre_ParCSRMatrixOwnsColStarts(hP);
1012 R_owns_row_starts = hypre_ParCSRMatrixOwnsRowStarts(hR);
1013 #endif
1014 /* can be replaced by version test later */
1015 #if defined(PETSC_HAVE_HYPRE_DEVICE)
1016 PetscStackPushExternal("hypre_ParCSRMatrixRAP");
1017 *hRAP = hypre_ParCSRMatrixRAP(hR, hA, hP);
1018 PetscStackPop;
1019 #else
1020 PetscCallHYPRE(hypre_BoomerAMGBuildCoarseOperator(hR, hA, hP, hRAP));
1021 PetscCallHYPRE(hypre_ParCSRMatrixSetNumNonzeros(*hRAP));
1022 #endif
1023 /* hypre_BoomerAMGBuildCoarseOperator steals the col_starts from P and the row_starts from R */
1024 #if defined(hypre_ParCSRMatrixOwnsRowStarts)
1025 hypre_ParCSRMatrixSetRowStartsOwner(*hRAP, 0);
1026 hypre_ParCSRMatrixSetColStartsOwner(*hRAP, 0);
1027 if (P_owns_col_starts) hypre_ParCSRMatrixSetColStartsOwner(hP, 1);
1028 if (R_owns_row_starts) hypre_ParCSRMatrixSetRowStartsOwner(hR, 1);
1029 #endif
1030 PetscFunctionReturn(PETSC_SUCCESS);
1031 }
1032
1033 static PetscErrorCode MatPtAPNumeric_AIJ_AIJ_wHYPRE(Mat A, Mat P, Mat C)
1034 {
1035 Mat B;
1036 hypre_ParCSRMatrix *hA, *hP, *hPtAP = NULL;
1037 Mat_Product *product = C->product;
1038
1039 PetscFunctionBegin;
1040 PetscCall(MatAIJGetParCSR_Private(A, &hA));
1041 PetscCall(MatAIJGetParCSR_Private(P, &hP));
1042 PetscCall(MatHYPRE_ParCSR_RAP(hP, hA, hP, &hPtAP));
1043 PetscCall(MatCreateFromParCSR(hPtAP, MATAIJ, PETSC_OWN_POINTER, &B));
1044
1045 PetscCall(MatHeaderMerge(C, &B));
1046 C->product = product;
1047
1048 PetscCall(MatAIJRestoreParCSR_Private(A, &hA));
1049 PetscCall(MatAIJRestoreParCSR_Private(P, &hP));
1050 PetscFunctionReturn(PETSC_SUCCESS);
1051 }
1052
1053 PETSC_INTERN PetscErrorCode MatPtAPSymbolic_AIJ_AIJ_wHYPRE(Mat A, Mat P, PetscReal fill, Mat C)
1054 {
1055 PetscFunctionBegin;
1056 PetscCall(MatSetType(C, MATAIJ));
1057 C->ops->ptapnumeric = MatPtAPNumeric_AIJ_AIJ_wHYPRE;
1058 C->ops->productnumeric = MatProductNumeric_PtAP;
1059 PetscFunctionReturn(PETSC_SUCCESS);
1060 }
1061
1062 static PetscErrorCode MatPtAPNumeric_AIJ_HYPRE(Mat A, Mat P, Mat C)
1063 {
1064 Mat B;
1065 Mat_HYPRE *hP;
1066 hypre_ParCSRMatrix *hA = NULL, *Pparcsr, *ptapparcsr = NULL;
1067 HYPRE_Int type;
1068 MPI_Comm comm = PetscObjectComm((PetscObject)A);
1069 PetscBool ishypre;
1070
1071 PetscFunctionBegin;
1072 PetscCall(PetscObjectTypeCompare((PetscObject)P, MATHYPRE, &ishypre));
1073 PetscCheck(ishypre, comm, PETSC_ERR_USER, "P should be of type %s", MATHYPRE);
1074 hP = (Mat_HYPRE *)P->data;
1075 PetscCallHYPRE(HYPRE_IJMatrixGetObjectType(hP->ij, &type));
1076 PetscCheck(type == HYPRE_PARCSR, comm, PETSC_ERR_SUP, "Only HYPRE_PARCSR is supported");
1077 PetscCallHYPRE(HYPRE_IJMatrixGetObject(hP->ij, (void **)&Pparcsr));
1078
1079 PetscCall(MatAIJGetParCSR_Private(A, &hA));
1080 PetscCall(MatHYPRE_ParCSR_RAP(Pparcsr, hA, Pparcsr, &ptapparcsr));
1081 PetscCall(MatAIJRestoreParCSR_Private(A, &hA));
1082
1083 /* create temporary matrix and merge to C */
1084 PetscCall(MatCreateFromParCSR(ptapparcsr, ((PetscObject)C)->type_name, PETSC_OWN_POINTER, &B));
1085 PetscCall(MatHeaderMerge(C, &B));
1086 PetscFunctionReturn(PETSC_SUCCESS);
1087 }
1088
1089 static PetscErrorCode MatPtAPNumeric_HYPRE_HYPRE(Mat A, Mat P, Mat C)
1090 {
1091 Mat B;
1092 hypre_ParCSRMatrix *Aparcsr, *Pparcsr, *ptapparcsr = NULL;
1093 Mat_HYPRE *hA, *hP;
1094 PetscBool ishypre;
1095 HYPRE_Int type;
1096
1097 PetscFunctionBegin;
1098 PetscCall(PetscObjectTypeCompare((PetscObject)P, MATHYPRE, &ishypre));
1099 PetscCheck(ishypre, PetscObjectComm((PetscObject)P), PETSC_ERR_USER, "P should be of type %s", MATHYPRE);
1100 PetscCall(PetscObjectTypeCompare((PetscObject)A, MATHYPRE, &ishypre));
1101 PetscCheck(ishypre, PetscObjectComm((PetscObject)A), PETSC_ERR_USER, "A should be of type %s", MATHYPRE);
1102 hA = (Mat_HYPRE *)A->data;
1103 hP = (Mat_HYPRE *)P->data;
1104 PetscCallHYPRE(HYPRE_IJMatrixGetObjectType(hA->ij, &type));
1105 PetscCheck(type == HYPRE_PARCSR, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Only HYPRE_PARCSR is supported");
1106 PetscCallHYPRE(HYPRE_IJMatrixGetObjectType(hP->ij, &type));
1107 PetscCheck(type == HYPRE_PARCSR, PetscObjectComm((PetscObject)P), PETSC_ERR_SUP, "Only HYPRE_PARCSR is supported");
1108 PetscCallHYPRE(HYPRE_IJMatrixGetObject(hA->ij, (void **)&Aparcsr));
1109 PetscCallHYPRE(HYPRE_IJMatrixGetObject(hP->ij, (void **)&Pparcsr));
1110 PetscCall(MatHYPRE_ParCSR_RAP(Pparcsr, Aparcsr, Pparcsr, &ptapparcsr));
1111 PetscCall(MatCreateFromParCSR(ptapparcsr, MATHYPRE, PETSC_OWN_POINTER, &B));
1112 PetscCall(MatHeaderMerge(C, &B));
1113 PetscFunctionReturn(PETSC_SUCCESS);
1114 }
1115
1116 /* calls hypre_ParMatmul
1117 hypre_ParMatMul uses hypre_ParMatrixCreate with the communicator of hA
1118 hypre_ParMatrixCreate does not duplicate the communicator
1119 It looks like we don't need to have the diagonal entries ordered first */
1120 static PetscErrorCode MatHYPRE_ParCSR_MatMatMult(hypre_ParCSRMatrix *hA, hypre_ParCSRMatrix *hB, hypre_ParCSRMatrix **hAB)
1121 {
1122 PetscFunctionBegin;
1123 /* can be replaced by version test later */
1124 #if defined(PETSC_HAVE_HYPRE_DEVICE)
1125 PetscStackPushExternal("hypre_ParCSRMatMat");
1126 *hAB = hypre_ParCSRMatMat(hA, hB);
1127 #else
1128 PetscStackPushExternal("hypre_ParMatmul");
1129 *hAB = hypre_ParMatmul(hA, hB);
1130 #endif
1131 PetscStackPop;
1132 PetscFunctionReturn(PETSC_SUCCESS);
1133 }
1134
1135 static PetscErrorCode MatMatMultNumeric_AIJ_AIJ_wHYPRE(Mat A, Mat B, Mat C)
1136 {
1137 Mat D;
1138 hypre_ParCSRMatrix *hA, *hB, *hAB = NULL;
1139 Mat_Product *product = C->product;
1140
1141 PetscFunctionBegin;
1142 PetscCall(MatAIJGetParCSR_Private(A, &hA));
1143 PetscCall(MatAIJGetParCSR_Private(B, &hB));
1144 PetscCall(MatHYPRE_ParCSR_MatMatMult(hA, hB, &hAB));
1145 PetscCall(MatCreateFromParCSR(hAB, MATAIJ, PETSC_OWN_POINTER, &D));
1146
1147 PetscCall(MatHeaderMerge(C, &D));
1148 C->product = product;
1149
1150 PetscCall(MatAIJRestoreParCSR_Private(A, &hA));
1151 PetscCall(MatAIJRestoreParCSR_Private(B, &hB));
1152 PetscFunctionReturn(PETSC_SUCCESS);
1153 }
1154
1155 PETSC_INTERN PetscErrorCode MatMatMultSymbolic_AIJ_AIJ_wHYPRE(Mat A, Mat B, PetscReal fill, Mat C)
1156 {
1157 PetscFunctionBegin;
1158 PetscCall(MatSetType(C, MATAIJ));
1159 C->ops->matmultnumeric = MatMatMultNumeric_AIJ_AIJ_wHYPRE;
1160 C->ops->productnumeric = MatProductNumeric_AB;
1161 PetscFunctionReturn(PETSC_SUCCESS);
1162 }
1163
1164 static PetscErrorCode MatMatMultNumeric_HYPRE_HYPRE(Mat A, Mat B, Mat C)
1165 {
1166 Mat D;
1167 hypre_ParCSRMatrix *Aparcsr, *Bparcsr, *ABparcsr = NULL;
1168 Mat_HYPRE *hA, *hB;
1169 PetscBool ishypre;
1170 HYPRE_Int type;
1171 Mat_Product *product;
1172
1173 PetscFunctionBegin;
1174 PetscCall(PetscObjectTypeCompare((PetscObject)B, MATHYPRE, &ishypre));
1175 PetscCheck(ishypre, PetscObjectComm((PetscObject)B), PETSC_ERR_USER, "B should be of type %s", MATHYPRE);
1176 PetscCall(PetscObjectTypeCompare((PetscObject)A, MATHYPRE, &ishypre));
1177 PetscCheck(ishypre, PetscObjectComm((PetscObject)A), PETSC_ERR_USER, "A should be of type %s", MATHYPRE);
1178 hA = (Mat_HYPRE *)A->data;
1179 hB = (Mat_HYPRE *)B->data;
1180 PetscCallHYPRE(HYPRE_IJMatrixGetObjectType(hA->ij, &type));
1181 PetscCheck(type == HYPRE_PARCSR, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Only HYPRE_PARCSR is supported");
1182 PetscCallHYPRE(HYPRE_IJMatrixGetObjectType(hB->ij, &type));
1183 PetscCheck(type == HYPRE_PARCSR, PetscObjectComm((PetscObject)B), PETSC_ERR_SUP, "Only HYPRE_PARCSR is supported");
1184 PetscCallHYPRE(HYPRE_IJMatrixGetObject(hA->ij, (void **)&Aparcsr));
1185 PetscCallHYPRE(HYPRE_IJMatrixGetObject(hB->ij, (void **)&Bparcsr));
1186 PetscCall(MatHYPRE_ParCSR_MatMatMult(Aparcsr, Bparcsr, &ABparcsr));
1187 PetscCall(MatCreateFromParCSR(ABparcsr, MATHYPRE, PETSC_OWN_POINTER, &D));
1188
1189 /* need to use HeaderReplace because HeaderMerge messes up with the communicator */
1190 product = C->product; /* save it from MatHeaderReplace() */
1191 C->product = NULL;
1192 PetscCall(MatHeaderReplace(C, &D));
1193 C->product = product;
1194 C->ops->matmultnumeric = MatMatMultNumeric_HYPRE_HYPRE;
1195 C->ops->productnumeric = MatProductNumeric_AB;
1196 PetscFunctionReturn(PETSC_SUCCESS);
1197 }
1198
1199 PETSC_INTERN PetscErrorCode MatTransposeMatMatMultNumeric_AIJ_AIJ_AIJ_wHYPRE(Mat A, Mat B, Mat C, Mat D)
1200 {
1201 Mat E;
1202 hypre_ParCSRMatrix *hA, *hB, *hC, *hABC = NULL;
1203
1204 PetscFunctionBegin;
1205 PetscCall(MatAIJGetParCSR_Private(A, &hA));
1206 PetscCall(MatAIJGetParCSR_Private(B, &hB));
1207 PetscCall(MatAIJGetParCSR_Private(C, &hC));
1208 PetscCall(MatHYPRE_ParCSR_RAP(hA, hB, hC, &hABC));
1209 PetscCall(MatCreateFromParCSR(hABC, MATAIJ, PETSC_OWN_POINTER, &E));
1210 PetscCall(MatHeaderMerge(D, &E));
1211 PetscCall(MatAIJRestoreParCSR_Private(A, &hA));
1212 PetscCall(MatAIJRestoreParCSR_Private(B, &hB));
1213 PetscCall(MatAIJRestoreParCSR_Private(C, &hC));
1214 PetscFunctionReturn(PETSC_SUCCESS);
1215 }
1216
1217 PETSC_INTERN PetscErrorCode MatTransposeMatMatMultSymbolic_AIJ_AIJ_AIJ_wHYPRE(Mat A, Mat B, Mat C, PetscReal fill, Mat D)
1218 {
1219 PetscFunctionBegin;
1220 PetscCall(MatSetType(D, MATAIJ));
1221 PetscFunctionReturn(PETSC_SUCCESS);
1222 }
1223
1224 static PetscErrorCode MatProductSymbolic_AB_HYPRE(Mat C)
1225 {
1226 PetscFunctionBegin;
1227 C->ops->productnumeric = MatProductNumeric_AB;
1228 PetscFunctionReturn(PETSC_SUCCESS);
1229 }
1230
1231 static PetscErrorCode MatProductSetFromOptions_HYPRE_AB(Mat C)
1232 {
1233 Mat_Product *product = C->product;
1234 PetscBool Ahypre;
1235
1236 PetscFunctionBegin;
1237 PetscCall(PetscObjectTypeCompare((PetscObject)product->A, MATHYPRE, &Ahypre));
1238 if (Ahypre) { /* A is a Hypre matrix */
1239 PetscCall(MatSetType(C, MATHYPRE));
1240 C->ops->productsymbolic = MatProductSymbolic_AB_HYPRE;
1241 C->ops->matmultnumeric = MatMatMultNumeric_HYPRE_HYPRE;
1242 PetscFunctionReturn(PETSC_SUCCESS);
1243 }
1244 PetscFunctionReturn(PETSC_SUCCESS);
1245 }
1246
1247 static PetscErrorCode MatProductSymbolic_PtAP_HYPRE(Mat C)
1248 {
1249 PetscFunctionBegin;
1250 C->ops->productnumeric = MatProductNumeric_PtAP;
1251 PetscFunctionReturn(PETSC_SUCCESS);
1252 }
1253
1254 static PetscErrorCode MatProductSetFromOptions_HYPRE_PtAP(Mat C)
1255 {
1256 Mat_Product *product = C->product;
1257 PetscBool flg;
1258 PetscInt type = 0;
1259 const char *outTypes[4] = {"aij", "seqaij", "mpiaij", "hypre"};
1260 PetscInt ntype = 4;
1261 Mat A = product->A;
1262 PetscBool Ahypre;
1263
1264 PetscFunctionBegin;
1265 PetscCall(PetscObjectTypeCompare((PetscObject)A, MATHYPRE, &Ahypre));
1266 if (Ahypre) { /* A is a Hypre matrix */
1267 PetscCall(MatSetType(C, MATHYPRE));
1268 C->ops->productsymbolic = MatProductSymbolic_PtAP_HYPRE;
1269 C->ops->ptapnumeric = MatPtAPNumeric_HYPRE_HYPRE;
1270 PetscFunctionReturn(PETSC_SUCCESS);
1271 }
1272
1273 /* A is AIJ, P is Hypre, C = PtAP can be either AIJ or Hypre format */
1274 /* Get runtime option */
1275 if (product->api_user) {
1276 PetscOptionsBegin(PetscObjectComm((PetscObject)C), ((PetscObject)C)->prefix, "MatPtAP_HYPRE", "Mat");
1277 PetscCall(PetscOptionsEList("-matptap_hypre_outtype", "MatPtAP outtype", "MatPtAP outtype", outTypes, ntype, outTypes[type], &type, &flg));
1278 PetscOptionsEnd();
1279 } else {
1280 PetscOptionsBegin(PetscObjectComm((PetscObject)C), ((PetscObject)C)->prefix, "MatProduct_PtAP_HYPRE", "Mat");
1281 PetscCall(PetscOptionsEList("-mat_product_algorithm_hypre_outtype", "MatProduct_PtAP outtype", "MatProduct_PtAP", outTypes, ntype, outTypes[type], &type, &flg));
1282 PetscOptionsEnd();
1283 }
1284
1285 if (type == 0 || type == 1 || type == 2) {
1286 PetscCall(MatSetType(C, MATAIJ));
1287 } else if (type == 3) {
1288 PetscCall(MatSetType(C, MATHYPRE));
1289 } else SETERRQ(PetscObjectComm((PetscObject)C), PETSC_ERR_SUP, "MatPtAP outtype is not supported");
1290 C->ops->productsymbolic = MatProductSymbolic_PtAP_HYPRE;
1291 C->ops->ptapnumeric = MatPtAPNumeric_AIJ_HYPRE;
1292 PetscFunctionReturn(PETSC_SUCCESS);
1293 }
1294
1295 static PetscErrorCode MatProductSetFromOptions_HYPRE(Mat C)
1296 {
1297 Mat_Product *product = C->product;
1298
1299 PetscFunctionBegin;
1300 switch (product->type) {
1301 case MATPRODUCT_AB:
1302 PetscCall(MatProductSetFromOptions_HYPRE_AB(C));
1303 break;
1304 case MATPRODUCT_PtAP:
1305 PetscCall(MatProductSetFromOptions_HYPRE_PtAP(C));
1306 break;
1307 default:
1308 break;
1309 }
1310 PetscFunctionReturn(PETSC_SUCCESS);
1311 }
1312
1313 static PetscErrorCode MatMultTranspose_HYPRE(Mat A, Vec x, Vec y)
1314 {
1315 PetscFunctionBegin;
1316 PetscCall(MatHYPRE_MultKernel_Private(A, 1.0, x, 0.0, y, PETSC_TRUE));
1317 PetscFunctionReturn(PETSC_SUCCESS);
1318 }
1319
1320 static PetscErrorCode MatMult_HYPRE(Mat A, Vec x, Vec y)
1321 {
1322 PetscFunctionBegin;
1323 PetscCall(MatHYPRE_MultKernel_Private(A, 1.0, x, 0.0, y, PETSC_FALSE));
1324 PetscFunctionReturn(PETSC_SUCCESS);
1325 }
1326
1327 static PetscErrorCode MatMultAdd_HYPRE(Mat A, Vec x, Vec y, Vec z)
1328 {
1329 PetscFunctionBegin;
1330 if (y != z) PetscCall(VecCopy(y, z));
1331 PetscCall(MatHYPRE_MultKernel_Private(A, 1.0, x, 1.0, z, PETSC_FALSE));
1332 PetscFunctionReturn(PETSC_SUCCESS);
1333 }
1334
1335 static PetscErrorCode MatMultTransposeAdd_HYPRE(Mat A, Vec x, Vec y, Vec z)
1336 {
1337 PetscFunctionBegin;
1338 if (y != z) PetscCall(VecCopy(y, z));
1339 PetscCall(MatHYPRE_MultKernel_Private(A, 1.0, x, 1.0, z, PETSC_TRUE));
1340 PetscFunctionReturn(PETSC_SUCCESS);
1341 }
1342
1343 /* y = a * A * x + b * y or y = a * A^t * x + b * y depending on trans */
1344 static PetscErrorCode MatHYPRE_MultKernel_Private(Mat A, HYPRE_Complex a, Vec x, HYPRE_Complex b, Vec y, PetscBool trans)
1345 {
1346 Mat_HYPRE *hA = (Mat_HYPRE *)A->data;
1347 hypre_ParCSRMatrix *parcsr;
1348 hypre_ParVector *hx, *hy;
1349
1350 PetscFunctionBegin;
1351 if (trans) {
1352 PetscCall(VecHYPRE_IJVectorPushVecRead(hA->b, x));
1353 if (b != 0.0) PetscCall(VecHYPRE_IJVectorPushVec(hA->x, y));
1354 else PetscCall(VecHYPRE_IJVectorPushVecWrite(hA->x, y));
1355 PetscCallHYPRE(HYPRE_IJVectorGetObject(hA->b->ij, (void **)&hx));
1356 PetscCallHYPRE(HYPRE_IJVectorGetObject(hA->x->ij, (void **)&hy));
1357 } else {
1358 PetscCall(VecHYPRE_IJVectorPushVecRead(hA->x, x));
1359 if (b != 0.0) PetscCall(VecHYPRE_IJVectorPushVec(hA->b, y));
1360 else PetscCall(VecHYPRE_IJVectorPushVecWrite(hA->b, y));
1361 PetscCallHYPRE(HYPRE_IJVectorGetObject(hA->x->ij, (void **)&hx));
1362 PetscCallHYPRE(HYPRE_IJVectorGetObject(hA->b->ij, (void **)&hy));
1363 }
1364 PetscCallHYPRE(HYPRE_IJMatrixGetObject(hA->ij, (void **)&parcsr));
1365 if (trans) {
1366 PetscCallHYPRE(hypre_ParCSRMatrixMatvecT(a, parcsr, hx, b, hy));
1367 } else {
1368 PetscCallHYPRE(hypre_ParCSRMatrixMatvec(a, parcsr, hx, b, hy));
1369 }
1370 PetscCall(VecHYPRE_IJVectorPopVec(hA->x));
1371 PetscCall(VecHYPRE_IJVectorPopVec(hA->b));
1372 PetscFunctionReturn(PETSC_SUCCESS);
1373 }
1374
1375 static PetscErrorCode MatDestroy_HYPRE(Mat A)
1376 {
1377 Mat_HYPRE *hA = (Mat_HYPRE *)A->data;
1378
1379 PetscFunctionBegin;
1380 PetscCall(VecHYPRE_IJVectorDestroy(&hA->x));
1381 PetscCall(VecHYPRE_IJVectorDestroy(&hA->b));
1382 PetscCall(MatHYPRE_DestroyCOOMat(A)); /* must be called before destroying the individual CSR */
1383 if (hA->ij) {
1384 if (!hA->inner_free) hypre_IJMatrixObject(hA->ij) = NULL;
1385 PetscCallHYPRE(HYPRE_IJMatrixDestroy(hA->ij));
1386 }
1387 if (hA->comm) PetscCall(PetscCommRestoreComm(PetscObjectComm((PetscObject)A), &hA->comm));
1388
1389 PetscCall(MatStashDestroy_Private(&A->stash));
1390 PetscCall(PetscFree(hA->array));
1391 if (hA->rows_d) PetscStackCallExternalVoid("hypre_Free", hypre_Free(hA->rows_d, HYPRE_MEMORY_DEVICE));
1392
1393 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_hypre_aij_C", NULL));
1394 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_hypre_is_C", NULL));
1395 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_seqaij_hypre_C", NULL));
1396 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_mpiaij_hypre_C", NULL));
1397 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_seqaijhipsparse_hypre_C", NULL));
1398 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_mpiaijhipsparse_hypre_C", NULL));
1399 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_seqaijcusparse_hypre_C", NULL));
1400 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_mpiaijcusparse_hypre_C", NULL));
1401 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHYPRESetPreallocation_C", NULL));
1402 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHYPREGetParCSR_C", NULL));
1403 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSetPreallocationCOO_C", NULL));
1404 PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSetValuesCOO_C", NULL));
1405 PetscCall(PetscFree(A->data));
1406 PetscFunctionReturn(PETSC_SUCCESS);
1407 }
1408
1409 static PetscErrorCode MatSetUp_HYPRE(Mat A)
1410 {
1411 PetscFunctionBegin;
1412 if (!A->preallocated) PetscCall(MatHYPRESetPreallocation(A, PETSC_DEFAULT, NULL, PETSC_DEFAULT, NULL));
1413 PetscFunctionReturn(PETSC_SUCCESS);
1414 }
1415
1416 //TODO FIX hypre_CSRMatrixMatvecOutOfPlace
1417 #if defined(PETSC_HAVE_HYPRE_DEVICE)
1418 static PetscErrorCode MatBindToCPU_HYPRE(Mat A, PetscBool bind)
1419 {
1420 Mat_HYPRE *hA = (Mat_HYPRE *)A->data;
1421 HYPRE_MemoryLocation hmem = bind ? HYPRE_MEMORY_HOST : HYPRE_MEMORY_DEVICE;
1422
1423 PetscFunctionBegin;
1424 A->boundtocpu = bind;
1425 if (hA->ij && hypre_IJMatrixAssembleFlag(hA->ij) && hmem != hypre_IJMatrixMemoryLocation(hA->ij)) {
1426 hypre_ParCSRMatrix *parcsr;
1427 PetscCallHYPRE(HYPRE_IJMatrixGetObject(hA->ij, (void **)&parcsr));
1428 PetscCallHYPRE(hypre_ParCSRMatrixMigrate(parcsr, hmem));
1429 }
1430 if (hA->x) PetscCall(VecHYPRE_IJBindToCPU(hA->x, bind));
1431 if (hA->b) PetscCall(VecHYPRE_IJBindToCPU(hA->b, bind));
1432 PetscFunctionReturn(PETSC_SUCCESS);
1433 }
1434 #endif
1435
1436 static PetscErrorCode MatAssemblyEnd_HYPRE(Mat A, MatAssemblyType mode)
1437 {
1438 Mat_HYPRE *hA = (Mat_HYPRE *)A->data;
1439 PetscMPIInt n;
1440 PetscInt i, j, rstart, ncols, flg;
1441 PetscInt *row, *col;
1442 PetscScalar *val;
1443
1444 PetscFunctionBegin;
1445 PetscCheck(mode != MAT_FLUSH_ASSEMBLY, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "MAT_FLUSH_ASSEMBLY currently not supported with MATHYPRE");
1446
1447 if (!A->nooffprocentries) {
1448 while (1) {
1449 PetscCall(MatStashScatterGetMesg_Private(&A->stash, &n, &row, &col, &val, &flg));
1450 if (!flg) break;
1451
1452 for (i = 0; i < n;) {
1453 /* Now identify the consecutive vals belonging to the same row */
1454 for (j = i, rstart = row[j]; j < n; j++) {
1455 if (row[j] != rstart) break;
1456 }
1457 if (j < n) ncols = j - i;
1458 else ncols = n - i;
1459 /* Now assemble all these values with a single function call */
1460 PetscCall(MatSetValues_HYPRE(A, 1, row + i, ncols, col + i, val + i, A->insertmode));
1461
1462 i = j;
1463 }
1464 }
1465 PetscCall(MatStashScatterEnd_Private(&A->stash));
1466 }
1467
1468 PetscCallHYPRE(HYPRE_IJMatrixAssemble(hA->ij));
1469 /* The assembly routine destroys the aux_matrix, we recreate it here by calling HYPRE_IJMatrixInitialize */
1470 /* If the option MAT_SORTED_FULL is set to true, the indices and values can be passed to hypre directly, so we don't need the aux_matrix */
1471 if (!A->sortedfull) {
1472 hypre_AuxParCSRMatrix *aux_matrix;
1473
1474 /* call destroy just to make sure we do not leak anything */
1475 aux_matrix = (hypre_AuxParCSRMatrix *)hypre_IJMatrixTranslator(hA->ij);
1476 PetscCallHYPRE(hypre_AuxParCSRMatrixDestroy(aux_matrix));
1477 hypre_IJMatrixTranslator(hA->ij) = NULL;
1478
1479 /* Initialize with assembled flag -> it only recreates the aux_par_matrix */
1480 PetscCallHYPRE(HYPRE_IJMatrixInitialize(hA->ij));
1481 aux_matrix = (hypre_AuxParCSRMatrix *)hypre_IJMatrixTranslator(hA->ij);
1482 if (aux_matrix) {
1483 hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 1; /* see comment in MatHYPRESetPreallocation_HYPRE */
1484 #if PETSC_PKG_HYPRE_VERSION_LT(2, 19, 0)
1485 PetscCallHYPRE(hypre_AuxParCSRMatrixInitialize(aux_matrix));
1486 #else
1487 PetscCallHYPRE(hypre_AuxParCSRMatrixInitialize_v2(aux_matrix, HYPRE_MEMORY_HOST));
1488 #endif
1489 }
1490 }
1491 {
1492 hypre_ParCSRMatrix *parcsr;
1493
1494 PetscCallHYPRE(HYPRE_IJMatrixGetObject(hA->ij, (void **)&parcsr));
1495 if (!hypre_ParCSRMatrixCommPkg(parcsr)) PetscCallHYPRE(hypre_MatvecCommPkgCreate(parcsr));
1496 }
1497 if (!hA->x) PetscCall(VecHYPRE_IJVectorCreate(A->cmap, &hA->x));
1498 if (!hA->b) PetscCall(VecHYPRE_IJVectorCreate(A->rmap, &hA->b));
1499 #if defined(PETSC_HAVE_HYPRE_DEVICE)
1500 PetscCall(MatBindToCPU_HYPRE(A, A->boundtocpu));
1501 #endif
1502 PetscFunctionReturn(PETSC_SUCCESS);
1503 }
1504
1505 static PetscErrorCode MatGetArray_HYPRE(Mat A, PetscInt size, void **array)
1506 {
1507 Mat_HYPRE *hA = (Mat_HYPRE *)A->data;
1508
1509 PetscFunctionBegin;
1510 PetscCheck(hA->array_available, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "Temporary space is in use");
1511
1512 if (hA->array_size >= size) {
1513 *array = hA->array;
1514 } else {
1515 PetscCall(PetscFree(hA->array));
1516 hA->array_size = size;
1517 PetscCall(PetscMalloc(hA->array_size, &hA->array));
1518 *array = hA->array;
1519 }
1520
1521 hA->array_available = PETSC_FALSE;
1522 PetscFunctionReturn(PETSC_SUCCESS);
1523 }
1524
1525 static PetscErrorCode MatRestoreArray_HYPRE(Mat A, void **array)
1526 {
1527 Mat_HYPRE *hA = (Mat_HYPRE *)A->data;
1528
1529 PetscFunctionBegin;
1530 *array = NULL;
1531 hA->array_available = PETSC_TRUE;
1532 PetscFunctionReturn(PETSC_SUCCESS);
1533 }
1534
1535 static PetscErrorCode MatSetValues_HYPRE(Mat A, PetscInt nr, const PetscInt rows[], PetscInt nc, const PetscInt cols[], const PetscScalar v[], InsertMode ins)
1536 {
1537 Mat_HYPRE *hA = (Mat_HYPRE *)A->data;
1538 PetscScalar *vals = (PetscScalar *)v;
1539 HYPRE_Complex *sscr;
1540 PetscInt *cscr[2];
1541 PetscInt i, nzc;
1542 PetscInt rst = A->rmap->rstart, ren = A->rmap->rend;
1543 void *array = NULL;
1544
1545 PetscFunctionBegin;
1546 PetscCall(MatGetArray_HYPRE(A, sizeof(PetscInt) * (2 * nc) + sizeof(HYPRE_Complex) * nc * nr, &array));
1547 cscr[0] = (PetscInt *)array;
1548 cscr[1] = ((PetscInt *)array) + nc;
1549 sscr = (HYPRE_Complex *)(((PetscInt *)array) + nc * 2);
1550 for (i = 0, nzc = 0; i < nc; i++) {
1551 if (cols[i] >= 0) {
1552 cscr[0][nzc] = cols[i];
1553 cscr[1][nzc++] = i;
1554 }
1555 }
1556 if (!nzc) {
1557 PetscCall(MatRestoreArray_HYPRE(A, &array));
1558 PetscFunctionReturn(PETSC_SUCCESS);
1559 }
1560
1561 #if 0 //defined(PETSC_HAVE_HYPRE_DEVICE)
1562 if (HYPRE_MEMORY_HOST != hypre_IJMatrixMemoryLocation(hA->ij)) {
1563 hypre_ParCSRMatrix *parcsr;
1564
1565 PetscCallHYPRE(HYPRE_IJMatrixGetObject(hA->ij,(void**)&parcsr));
1566 PetscCallHYPRE(hypre_ParCSRMatrixMigrate(parcsr, HYPRE_MEMORY_HOST));
1567 }
1568 #endif
1569
1570 if (ins == ADD_VALUES) {
1571 for (i = 0; i < nr; i++) {
1572 if (rows[i] >= 0) {
1573 PetscInt j;
1574 HYPRE_Int hnc = (HYPRE_Int)nzc;
1575
1576 if (!nzc) continue;
1577 /* nonlocal values */
1578 if (rows[i] < rst || rows[i] >= ren) {
1579 PetscCheck(!A->nooffprocentries, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Setting off process row %" PetscInt_FMT " even though MatSetOption(,MAT_NO_OFF_PROC_ENTRIES,PETSC_TRUE) was set", rows[i]);
1580 if (hA->donotstash) continue;
1581 }
1582 PetscCheck((PetscInt)hnc == nzc, PETSC_COMM_SELF, PETSC_ERR_SUP, "Hypre overflow! number of columns %" PetscInt_FMT " for row %" PetscInt_FMT, nzc, rows[i]);
1583 for (j = 0; j < nzc; j++) PetscCall(PetscHYPREScalarCast(vals[cscr[1][j]], &sscr[j]));
1584 PetscCallHYPRE(HYPRE_IJMatrixAddToValues(hA->ij, 1, &hnc, (HYPRE_BigInt *)(rows + i), (HYPRE_BigInt *)cscr[0], sscr));
1585 }
1586 vals += nc;
1587 }
1588 } else { /* INSERT_VALUES */
1589 for (i = 0; i < nr; i++) {
1590 if (rows[i] >= 0) {
1591 PetscInt j;
1592 HYPRE_Int hnc = (HYPRE_Int)nzc;
1593
1594 if (!nzc) continue;
1595 PetscCheck((PetscInt)hnc == nzc, PETSC_COMM_SELF, PETSC_ERR_SUP, "Hypre overflow! number of columns %" PetscInt_FMT " for row %" PetscInt_FMT, nzc, rows[i]);
1596 for (j = 0; j < nzc; j++) PetscCall(PetscHYPREScalarCast(vals[cscr[1][j]], &sscr[j]));
1597 /* nonlocal values */
1598 if (rows[i] < rst || rows[i] >= ren) {
1599 PetscCheck(!A->nooffprocentries, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Setting off process row %" PetscInt_FMT " even though MatSetOption(,MAT_NO_OFF_PROC_ENTRIES,PETSC_TRUE) was set", rows[i]);
1600 if (!hA->donotstash) PetscCall(MatStashValuesRow_Private(&A->stash, rows[i], nzc, cscr[0], (PetscScalar *)sscr, PETSC_FALSE));
1601 }
1602 /* local values */
1603 else
1604 PetscCallHYPRE(HYPRE_IJMatrixSetValues(hA->ij, 1, &hnc, (HYPRE_BigInt *)(rows + i), (HYPRE_BigInt *)cscr[0], sscr));
1605 }
1606 vals += nc;
1607 }
1608 }
1609
1610 PetscCall(MatRestoreArray_HYPRE(A, &array));
1611 PetscFunctionReturn(PETSC_SUCCESS);
1612 }
1613
1614 static PetscErrorCode MatHYPRESetPreallocation_HYPRE(Mat A, PetscInt dnz, const PetscInt dnnz[], PetscInt onz, const PetscInt onnz[])
1615 {
1616 Mat_HYPRE *hA = (Mat_HYPRE *)A->data;
1617 HYPRE_Int *hdnnz, *honnz;
1618 PetscInt i, rs, re, cs, ce, bs;
1619 PetscMPIInt size;
1620
1621 PetscFunctionBegin;
1622 PetscCall(PetscLayoutSetUp(A->rmap));
1623 PetscCall(PetscLayoutSetUp(A->cmap));
1624 rs = A->rmap->rstart;
1625 re = A->rmap->rend;
1626 cs = A->cmap->rstart;
1627 ce = A->cmap->rend;
1628 if (!hA->ij) {
1629 PetscCallHYPRE(HYPRE_IJMatrixCreate(hA->comm, rs, re - 1, cs, ce - 1, &hA->ij));
1630 PetscCallHYPRE(HYPRE_IJMatrixSetObjectType(hA->ij, HYPRE_PARCSR));
1631 } else {
1632 HYPRE_BigInt hrs, hre, hcs, hce;
1633 PetscCallHYPRE(HYPRE_IJMatrixGetLocalRange(hA->ij, &hrs, &hre, &hcs, &hce));
1634 PetscCheck(hre - hrs + 1 == re - rs, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent local rows: IJMatrix [%" PetscHYPRE_BigInt_FMT ",%" PetscHYPRE_BigInt_FMT "), PETSc [%" PetscInt_FMT ",%" PetscInt_FMT ")", hrs, hre + 1, rs, re);
1635 PetscCheck(hce - hcs + 1 == ce - cs, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent local cols: IJMatrix [%" PetscHYPRE_BigInt_FMT ",%" PetscHYPRE_BigInt_FMT "), PETSc [%" PetscInt_FMT ",%" PetscInt_FMT ")", hcs, hce + 1, cs, ce);
1636 }
1637 PetscCall(MatHYPRE_DestroyCOOMat(A));
1638 PetscCall(MatGetBlockSize(A, &bs));
1639 if (dnz == PETSC_DEFAULT || dnz == PETSC_DECIDE) dnz = 10 * bs;
1640 if (onz == PETSC_DEFAULT || onz == PETSC_DECIDE) onz = 10 * bs;
1641
1642 if (!dnnz) {
1643 PetscCall(PetscMalloc1(A->rmap->n, &hdnnz));
1644 for (i = 0; i < A->rmap->n; i++) hdnnz[i] = (HYPRE_Int)dnz;
1645 } else {
1646 hdnnz = (HYPRE_Int *)dnnz;
1647 }
1648 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
1649 if (size > 1) {
1650 hypre_AuxParCSRMatrix *aux_matrix;
1651 if (!onnz) {
1652 PetscCall(PetscMalloc1(A->rmap->n, &honnz));
1653 for (i = 0; i < A->rmap->n; i++) honnz[i] = (HYPRE_Int)onz;
1654 } else honnz = (HYPRE_Int *)onnz;
1655 /* SetDiagOffdSizes sets hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 0, since it seems
1656 they assume the user will input the entire row values, properly sorted
1657 In PETSc, we don't make such an assumption and set this flag to 1,
1658 unless the option MAT_SORTED_FULL is set to true.
1659 Also, to avoid possible memory leaks, we destroy and recreate the translator
1660 This has to be done here, as HYPRE_IJMatrixInitialize will properly initialize
1661 the IJ matrix for us */
1662 aux_matrix = (hypre_AuxParCSRMatrix *)hypre_IJMatrixTranslator(hA->ij);
1663 hypre_AuxParCSRMatrixDestroy(aux_matrix);
1664 hypre_IJMatrixTranslator(hA->ij) = NULL;
1665 PetscCallHYPRE(HYPRE_IJMatrixSetDiagOffdSizes(hA->ij, hdnnz, honnz));
1666 aux_matrix = (hypre_AuxParCSRMatrix *)hypre_IJMatrixTranslator(hA->ij);
1667 hypre_AuxParCSRMatrixNeedAux(aux_matrix) = !A->sortedfull;
1668 } else {
1669 honnz = NULL;
1670 PetscCallHYPRE(HYPRE_IJMatrixSetRowSizes(hA->ij, hdnnz));
1671 }
1672
1673 /* reset assembled flag and call the initialize method */
1674 hypre_IJMatrixAssembleFlag(hA->ij) = 0;
1675 #if PETSC_PKG_HYPRE_VERSION_LT(2, 19, 0)
1676 PetscCallHYPRE(HYPRE_IJMatrixInitialize(hA->ij));
1677 #else
1678 PetscCallHYPRE(HYPRE_IJMatrixInitialize_v2(hA->ij, HYPRE_MEMORY_HOST));
1679 #endif
1680 if (!dnnz) PetscCall(PetscFree(hdnnz));
1681 if (!onnz && honnz) PetscCall(PetscFree(honnz));
1682 /* Match AIJ logic */
1683 A->preallocated = PETSC_TRUE;
1684 A->assembled = PETSC_FALSE;
1685 PetscFunctionReturn(PETSC_SUCCESS);
1686 }
1687
1688 /*@C
1689 MatHYPRESetPreallocation - Preallocates memory for a sparse parallel matrix in HYPRE IJ format
1690
1691 Collective
1692
1693 Input Parameters:
1694 + A - the matrix
1695 . dnz - number of nonzeros per row in DIAGONAL portion of local submatrix
1696 (same value is used for all local rows)
1697 . dnnz - array containing the number of nonzeros in the various rows of the
1698 DIAGONAL portion of the local submatrix (possibly different for each row)
1699 or `NULL` (`PETSC_NULL_INTEGER` in Fortran), if `d_nz` is used to specify the nonzero structure.
1700 The size of this array is equal to the number of local rows, i.e `m`.
1701 For matrices that will be factored, you must leave room for (and set)
1702 the diagonal entry even if it is zero.
1703 . onz - number of nonzeros per row in the OFF-DIAGONAL portion of local
1704 submatrix (same value is used for all local rows).
1705 - onnz - array containing the number of nonzeros in the various rows of the
1706 OFF-DIAGONAL portion of the local submatrix (possibly different for
1707 each row) or `NULL` (`PETSC_NULL_INTEGER` in Fortran), if `o_nz` is used to specify the nonzero
1708 structure. The size of this array is equal to the number
1709 of local rows, i.e `m`.
1710
1711 Level: intermediate
1712
1713 Note:
1714 If the *nnz parameter is given then the *nz parameter is ignored; for sequential matrices, `onz` and `onnz` are ignored.
1715
1716 .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MatMPIAIJSetPreallocation()`, `MATHYPRE`, `MATAIJ`
1717 @*/
1718 PetscErrorCode MatHYPRESetPreallocation(Mat A, PetscInt dnz, const PetscInt dnnz[], PetscInt onz, const PetscInt onnz[])
1719 {
1720 PetscFunctionBegin;
1721 PetscValidHeaderSpecific(A, MAT_CLASSID, 1);
1722 PetscValidType(A, 1);
1723 PetscTryMethod(A, "MatHYPRESetPreallocation_C", (Mat, PetscInt, const PetscInt[], PetscInt, const PetscInt[]), (A, dnz, dnnz, onz, onnz));
1724 PetscFunctionReturn(PETSC_SUCCESS);
1725 }
1726
1727 /*@C
1728 MatCreateFromParCSR - Creates a `Mat` from a `hypre_ParCSRMatrix`
1729
1730 Collective
1731
1732 Input Parameters:
1733 + parcsr - the pointer to the `hypre_ParCSRMatrix`
1734 . mtype - matrix type to be created. Currently `MATAIJ`, `MATIS` and `MATHYPRE` are supported.
1735 - copymode - PETSc copying options, see `PetscCopyMode`
1736
1737 Output Parameter:
1738 . A - the matrix
1739
1740 Level: intermediate
1741
1742 .seealso: [](ch_matrices), `Mat`, `MATHYPRE`, `PetscCopyMode`
1743 @*/
1744 PETSC_EXTERN PetscErrorCode MatCreateFromParCSR(hypre_ParCSRMatrix *parcsr, MatType mtype, PetscCopyMode copymode, Mat *A)
1745 {
1746 Mat T;
1747 Mat_HYPRE *hA;
1748 MPI_Comm comm;
1749 PetscInt rstart, rend, cstart, cend, M, N;
1750 PetscBool isseqaij, isseqaijmkl, ismpiaij, isaij, ishyp, isis;
1751
1752 PetscFunctionBegin;
1753 comm = hypre_ParCSRMatrixComm(parcsr);
1754 PetscCall(PetscStrcmp(mtype, MATSEQAIJ, &isseqaij));
1755 PetscCall(PetscStrcmp(mtype, MATSEQAIJMKL, &isseqaijmkl));
1756 PetscCall(PetscStrcmp(mtype, MATMPIAIJ, &ismpiaij));
1757 PetscCall(PetscStrcmp(mtype, MATAIJ, &isaij));
1758 PetscCall(PetscStrcmp(mtype, MATHYPRE, &ishyp));
1759 PetscCall(PetscStrcmp(mtype, MATIS, &isis));
1760 isaij = (PetscBool)(isseqaij || isseqaijmkl || ismpiaij || isaij);
1761 /* TODO */
1762 PetscCheck(isaij || ishyp || isis, comm, PETSC_ERR_SUP, "Unsupported MatType %s! Supported types are %s, %s, %s, %s, %s, and %s", mtype, MATAIJ, MATSEQAIJ, MATSEQAIJMKL, MATMPIAIJ, MATIS, MATHYPRE);
1763 /* access ParCSRMatrix */
1764 rstart = hypre_ParCSRMatrixFirstRowIndex(parcsr);
1765 rend = hypre_ParCSRMatrixLastRowIndex(parcsr);
1766 cstart = hypre_ParCSRMatrixFirstColDiag(parcsr);
1767 cend = hypre_ParCSRMatrixLastColDiag(parcsr);
1768 M = hypre_ParCSRMatrixGlobalNumRows(parcsr);
1769 N = hypre_ParCSRMatrixGlobalNumCols(parcsr);
1770
1771 /* create PETSc matrix with MatHYPRE */
1772 PetscCall(MatCreate(comm, &T));
1773 PetscCall(MatSetSizes(T, PetscMax(rend - rstart + 1, 0), PetscMax(cend - cstart + 1, 0), M, N));
1774 PetscCall(MatSetType(T, MATHYPRE));
1775 hA = (Mat_HYPRE *)T->data;
1776
1777 /* create HYPRE_IJMatrix */
1778 PetscCallHYPRE(HYPRE_IJMatrixCreate(hA->comm, rstart, rend, cstart, cend, &hA->ij));
1779 PetscCallHYPRE(HYPRE_IJMatrixSetObjectType(hA->ij, HYPRE_PARCSR));
1780
1781 /* create new ParCSR object if needed */
1782 if (ishyp && copymode == PETSC_COPY_VALUES) {
1783 hypre_ParCSRMatrix *new_parcsr;
1784 #if PETSC_PKG_HYPRE_VERSION_LT(2, 18, 0)
1785 hypre_CSRMatrix *hdiag, *hoffd, *ndiag, *noffd;
1786
1787 new_parcsr = hypre_ParCSRMatrixClone(parcsr, 0);
1788 hdiag = hypre_ParCSRMatrixDiag(parcsr);
1789 hoffd = hypre_ParCSRMatrixOffd(parcsr);
1790 ndiag = hypre_ParCSRMatrixDiag(new_parcsr);
1791 noffd = hypre_ParCSRMatrixOffd(new_parcsr);
1792 PetscCall(PetscArraycpy(hypre_CSRMatrixData(ndiag), hypre_CSRMatrixData(hdiag), hypre_CSRMatrixNumNonzeros(hdiag)));
1793 PetscCall(PetscArraycpy(hypre_CSRMatrixData(noffd), hypre_CSRMatrixData(hoffd), hypre_CSRMatrixNumNonzeros(hoffd)));
1794 #else
1795 new_parcsr = hypre_ParCSRMatrixClone(parcsr, 1);
1796 #endif
1797 parcsr = new_parcsr;
1798 copymode = PETSC_OWN_POINTER;
1799 }
1800
1801 /* set ParCSR object */
1802 hypre_IJMatrixObject(hA->ij) = parcsr;
1803 T->preallocated = PETSC_TRUE;
1804
1805 /* set assembled flag */
1806 hypre_IJMatrixAssembleFlag(hA->ij) = 1;
1807 #if 0
1808 PetscCallHYPRE(HYPRE_IJMatrixInitialize(hA->ij));
1809 #endif
1810 if (ishyp) {
1811 PetscMPIInt myid = 0;
1812
1813 /* make sure we always have row_starts and col_starts available */
1814 if (HYPRE_AssumedPartitionCheck()) PetscCallMPI(MPI_Comm_rank(comm, &myid));
1815 #if defined(hypre_ParCSRMatrixOwnsRowStarts)
1816 if (!hypre_ParCSRMatrixOwnsColStarts(parcsr)) {
1817 PetscLayout map;
1818
1819 PetscCall(MatGetLayouts(T, NULL, &map));
1820 PetscCall(PetscLayoutSetUp(map));
1821 hypre_ParCSRMatrixColStarts(parcsr) = (HYPRE_BigInt *)(map->range + myid);
1822 }
1823 if (!hypre_ParCSRMatrixOwnsRowStarts(parcsr)) {
1824 PetscLayout map;
1825
1826 PetscCall(MatGetLayouts(T, &map, NULL));
1827 PetscCall(PetscLayoutSetUp(map));
1828 hypre_ParCSRMatrixRowStarts(parcsr) = (HYPRE_BigInt *)(map->range + myid);
1829 }
1830 #endif
1831 /* prevent from freeing the pointer */
1832 if (copymode == PETSC_USE_POINTER) hA->inner_free = PETSC_FALSE;
1833 *A = T;
1834 PetscCall(MatSetOption(*A, MAT_SORTED_FULL, PETSC_TRUE));
1835 PetscCall(MatAssemblyBegin(*A, MAT_FINAL_ASSEMBLY));
1836 PetscCall(MatAssemblyEnd(*A, MAT_FINAL_ASSEMBLY));
1837 } else if (isaij) {
1838 if (copymode != PETSC_OWN_POINTER) {
1839 /* prevent from freeing the pointer */
1840 hA->inner_free = PETSC_FALSE;
1841 PetscCall(MatConvert_HYPRE_AIJ(T, MATAIJ, MAT_INITIAL_MATRIX, A));
1842 PetscCall(MatDestroy(&T));
1843 } else { /* AIJ return type with PETSC_OWN_POINTER */
1844 PetscCall(MatConvert_HYPRE_AIJ(T, MATAIJ, MAT_INPLACE_MATRIX, &T));
1845 *A = T;
1846 }
1847 } else if (isis) {
1848 PetscCall(MatConvert_HYPRE_IS(T, MATIS, MAT_INITIAL_MATRIX, A));
1849 if (copymode != PETSC_OWN_POINTER) hA->inner_free = PETSC_FALSE;
1850 PetscCall(MatDestroy(&T));
1851 }
1852 PetscFunctionReturn(PETSC_SUCCESS);
1853 }
1854
1855 static PetscErrorCode MatHYPREGetParCSR_HYPRE(Mat A, hypre_ParCSRMatrix **parcsr)
1856 {
1857 Mat_HYPRE *hA = (Mat_HYPRE *)A->data;
1858 HYPRE_Int type;
1859
1860 PetscFunctionBegin;
1861 PetscCheck(hA->ij, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "HYPRE_IJMatrix not present");
1862 PetscCallHYPRE(HYPRE_IJMatrixGetObjectType(hA->ij, &type));
1863 PetscCheck(type == HYPRE_PARCSR, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "HYPRE_IJMatrix is not of type HYPRE_PARCSR");
1864 PetscCallHYPRE(HYPRE_IJMatrixGetObject(hA->ij, (void **)parcsr));
1865 PetscFunctionReturn(PETSC_SUCCESS);
1866 }
1867
1868 /*@C
1869 MatHYPREGetParCSR - Gets the pointer to the ParCSR matrix
1870
1871 Not Collective, No Fortran Support
1872
1873 Input Parameter:
1874 . A - the `MATHYPRE` object
1875
1876 Output Parameter:
1877 . parcsr - the pointer to the `hypre_ParCSRMatrix`
1878
1879 Level: intermediate
1880
1881 .seealso: [](ch_matrices), `Mat`, `MATHYPRE`, `PetscCopyMode`
1882 @*/
1883 PetscErrorCode MatHYPREGetParCSR(Mat A, hypre_ParCSRMatrix **parcsr)
1884 {
1885 PetscFunctionBegin;
1886 PetscValidHeaderSpecific(A, MAT_CLASSID, 1);
1887 PetscValidType(A, 1);
1888 PetscUseMethod(A, "MatHYPREGetParCSR_C", (Mat, hypre_ParCSRMatrix **), (A, parcsr));
1889 PetscFunctionReturn(PETSC_SUCCESS);
1890 }
1891
1892 static PetscErrorCode MatScale_HYPRE(Mat A, PetscScalar s)
1893 {
1894 hypre_ParCSRMatrix *parcsr;
1895 #if PETSC_PKG_HYPRE_VERSION_LT(2, 19, 0)
1896 hypre_CSRMatrix *ha;
1897 #endif
1898 HYPRE_Complex hs;
1899
1900 PetscFunctionBegin;
1901 PetscCall(PetscHYPREScalarCast(s, &hs));
1902 PetscCall(MatHYPREGetParCSR_HYPRE(A, &parcsr));
1903 #if PETSC_PKG_HYPRE_VERSION_GE(2, 19, 0)
1904 PetscCallHYPRE(hypre_ParCSRMatrixScale(parcsr, hs));
1905 #else /* diagonal part */
1906 ha = hypre_ParCSRMatrixDiag(parcsr);
1907 if (ha) {
1908 PetscInt size, i;
1909 HYPRE_Int *ii;
1910 HYPRE_Complex *a;
1911
1912 size = hypre_CSRMatrixNumRows(ha);
1913 a = hypre_CSRMatrixData(ha);
1914 ii = hypre_CSRMatrixI(ha);
1915 for (i = 0; i < ii[size]; i++) a[i] *= hs;
1916 }
1917 /* off-diagonal part */
1918 ha = hypre_ParCSRMatrixOffd(parcsr);
1919 if (ha) {
1920 PetscInt size, i;
1921 HYPRE_Int *ii;
1922 HYPRE_Complex *a;
1923
1924 size = hypre_CSRMatrixNumRows(ha);
1925 a = hypre_CSRMatrixData(ha);
1926 ii = hypre_CSRMatrixI(ha);
1927 for (i = 0; i < ii[size]; i++) a[i] *= hs;
1928 }
1929 #endif
1930 PetscFunctionReturn(PETSC_SUCCESS);
1931 }
1932
1933 static PetscErrorCode MatZeroRowsColumns_HYPRE(Mat A, PetscInt numRows, const PetscInt rows[], PetscScalar diag, Vec x, Vec b)
1934 {
1935 hypre_ParCSRMatrix *parcsr;
1936 HYPRE_Int *lrows;
1937 PetscInt rst, ren, i;
1938
1939 PetscFunctionBegin;
1940 PetscCheck(!x && !b, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "To be implemented");
1941 PetscCall(MatHYPREGetParCSR_HYPRE(A, &parcsr));
1942 PetscCall(PetscMalloc1(numRows, &lrows));
1943 PetscCall(MatGetOwnershipRange(A, &rst, &ren));
1944 for (i = 0; i < numRows; i++) {
1945 PetscCheck(rows[i] >= rst && rows[i] < ren, PETSC_COMM_SELF, PETSC_ERR_SUP, "Non-local rows not yet supported");
1946 lrows[i] = (HYPRE_Int)(rows[i] - rst);
1947 }
1948 PetscCallHYPRE(hypre_ParCSRMatrixEliminateRowsCols(parcsr, (HYPRE_Int)numRows, lrows));
1949 PetscCall(PetscFree(lrows));
1950 PetscFunctionReturn(PETSC_SUCCESS);
1951 }
1952
1953 static PetscErrorCode MatZeroEntries_HYPRE_CSRMatrix(hypre_CSRMatrix *ha)
1954 {
1955 PetscFunctionBegin;
1956 if (ha) {
1957 HYPRE_Int *ii, size;
1958 HYPRE_Complex *a;
1959
1960 size = hypre_CSRMatrixNumRows(ha);
1961 a = hypre_CSRMatrixData(ha);
1962 ii = hypre_CSRMatrixI(ha);
1963
1964 if (a) PetscCall(PetscArrayzero(a, ii[size]));
1965 }
1966 PetscFunctionReturn(PETSC_SUCCESS);
1967 }
1968
1969 static PetscErrorCode MatZeroEntries_HYPRE(Mat A)
1970 {
1971 Mat_HYPRE *hA = (Mat_HYPRE *)A->data;
1972
1973 PetscFunctionBegin;
1974 if (HYPRE_MEMORY_DEVICE == hypre_IJMatrixMemoryLocation(hA->ij)) {
1975 PetscCallHYPRE(HYPRE_IJMatrixSetConstantValues(hA->ij, 0.0));
1976 } else {
1977 hypre_ParCSRMatrix *parcsr;
1978
1979 PetscCall(MatHYPREGetParCSR_HYPRE(A, &parcsr));
1980 PetscCall(MatZeroEntries_HYPRE_CSRMatrix(hypre_ParCSRMatrixDiag(parcsr)));
1981 PetscCall(MatZeroEntries_HYPRE_CSRMatrix(hypre_ParCSRMatrixOffd(parcsr)));
1982 }
1983 PetscFunctionReturn(PETSC_SUCCESS);
1984 }
1985
1986 static PetscErrorCode MatZeroRows_HYPRE_CSRMatrix(hypre_CSRMatrix *hA, PetscInt N, const PetscInt rows[], HYPRE_Complex diag)
1987 {
1988 PetscInt ii;
1989 HYPRE_Int *i, *j;
1990 HYPRE_Complex *a;
1991
1992 PetscFunctionBegin;
1993 if (!hA) PetscFunctionReturn(PETSC_SUCCESS);
1994
1995 i = hypre_CSRMatrixI(hA);
1996 j = hypre_CSRMatrixJ(hA);
1997 a = hypre_CSRMatrixData(hA);
1998 #if defined(PETSC_HAVE_HYPRE_DEVICE)
1999 if (HYPRE_MEMORY_DEVICE == hypre_CSRMatrixMemoryLocation(hA)) {
2000 #if defined(HYPRE_USING_CUDA)
2001 PetscCall(MatZeroRows_CUDA(N, rows, i, j, a, diag));
2002 #elif defined(HYPRE_USING_HIP)
2003 PetscCall(MatZeroRows_HIP(N, rows, i, j, a, diag));
2004 #elif defined(PETSC_HAVE_KOKKOS)
2005 PetscCall(MatZeroRows_Kokkos(N, rows, i, j, a, diag));
2006 #else
2007 SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for MatZeroRows on a hypre matrix in this memory location");
2008 #endif
2009 } else
2010 #endif
2011 {
2012 for (ii = 0; ii < N; ii++) {
2013 HYPRE_Int jj, ibeg, iend, irow;
2014
2015 irow = (HYPRE_Int)rows[ii];
2016 ibeg = i[irow];
2017 iend = i[irow + 1];
2018 for (jj = ibeg; jj < iend; jj++)
2019 if (j[jj] == irow) a[jj] = diag;
2020 else a[jj] = 0.0;
2021 }
2022 }
2023 PetscFunctionReturn(PETSC_SUCCESS);
2024 }
2025
2026 static PetscErrorCode MatZeroRows_HYPRE(Mat A, PetscInt N, const PetscInt rows[], PetscScalar diag, Vec x, Vec b)
2027 {
2028 hypre_ParCSRMatrix *parcsr;
2029 PetscInt *lrows, len, *lrows2;
2030 HYPRE_Complex hdiag;
2031
2032 PetscFunctionBegin;
2033 PetscCheck(!x && !b, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Does not support to modify the solution and the right hand size");
2034 PetscCall(PetscHYPREScalarCast(diag, &hdiag));
2035 /* retrieve the internal matrix */
2036 PetscCall(MatHYPREGetParCSR_HYPRE(A, &parcsr));
2037 /* get locally owned rows */
2038 PetscCall(MatZeroRowsMapLocal_Private(A, N, rows, &len, &lrows));
2039
2040 #if defined(PETSC_HAVE_HYPRE_DEVICE)
2041 if (HYPRE_MEMORY_DEVICE == hypre_CSRMatrixMemoryLocation(hypre_ParCSRMatrixDiag(parcsr))) {
2042 Mat_HYPRE *hA = (Mat_HYPRE *)A->data;
2043 PetscInt m;
2044 PetscCall(MatGetLocalSize(A, &m, NULL));
2045 if (!hA->rows_d) {
2046 hA->rows_d = hypre_TAlloc(PetscInt, m, HYPRE_MEMORY_DEVICE);
2047 if (m) PetscCheck(hA->rows_d, PETSC_COMM_SELF, PETSC_ERR_MEM, "HYPRE_TAlloc failed");
2048 }
2049 PetscCheck(len <= m, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Too many rows in rows[]");
2050 PetscStackCallExternalVoid("hypre_Memcpy", hypre_Memcpy(hA->rows_d, lrows, sizeof(PetscInt) * len, HYPRE_MEMORY_DEVICE, HYPRE_MEMORY_HOST));
2051 lrows2 = hA->rows_d;
2052 } else
2053 #endif
2054 {
2055 lrows2 = lrows;
2056 }
2057
2058 /* zero diagonal part */
2059 PetscCall(MatZeroRows_HYPRE_CSRMatrix(hypre_ParCSRMatrixDiag(parcsr), len, lrows2, hdiag));
2060 /* zero off-diagonal part */
2061 PetscCall(MatZeroRows_HYPRE_CSRMatrix(hypre_ParCSRMatrixOffd(parcsr), len, lrows2, 0.0));
2062
2063 PetscCall(PetscFree(lrows));
2064 PetscFunctionReturn(PETSC_SUCCESS);
2065 }
2066
2067 static PetscErrorCode MatAssemblyBegin_HYPRE(Mat mat, MatAssemblyType mode)
2068 {
2069 PetscFunctionBegin;
2070 if (mat->nooffprocentries) PetscFunctionReturn(PETSC_SUCCESS);
2071
2072 PetscCall(MatStashScatterBegin_Private(mat, &mat->stash, mat->rmap->range));
2073 PetscFunctionReturn(PETSC_SUCCESS);
2074 }
2075
2076 static PetscErrorCode MatGetRow_HYPRE(Mat A, PetscInt row, PetscInt *nz, PetscInt **idx, PetscScalar **v)
2077 {
2078 hypre_ParCSRMatrix *parcsr;
2079 HYPRE_Int hnz;
2080
2081 PetscFunctionBegin;
2082 /* retrieve the internal matrix */
2083 PetscCall(MatHYPREGetParCSR_HYPRE(A, &parcsr));
2084 /* call HYPRE API */
2085 PetscCallHYPRE(HYPRE_ParCSRMatrixGetRow(parcsr, row, &hnz, (HYPRE_BigInt **)idx, (HYPRE_Complex **)v));
2086 if (nz) *nz = (PetscInt)hnz;
2087 PetscFunctionReturn(PETSC_SUCCESS);
2088 }
2089
2090 static PetscErrorCode MatRestoreRow_HYPRE(Mat A, PetscInt row, PetscInt *nz, PetscInt **idx, PetscScalar **v)
2091 {
2092 hypre_ParCSRMatrix *parcsr;
2093 HYPRE_Int hnz;
2094
2095 PetscFunctionBegin;
2096 /* retrieve the internal matrix */
2097 PetscCall(MatHYPREGetParCSR_HYPRE(A, &parcsr));
2098 /* call HYPRE API */
2099 hnz = nz ? (HYPRE_Int)(*nz) : 0;
2100 PetscCallHYPRE(HYPRE_ParCSRMatrixRestoreRow(parcsr, row, &hnz, (HYPRE_BigInt **)idx, (HYPRE_Complex **)v));
2101 PetscFunctionReturn(PETSC_SUCCESS);
2102 }
2103
2104 static PetscErrorCode MatGetValues_HYPRE(Mat A, PetscInt m, const PetscInt idxm[], PetscInt n, const PetscInt idxn[], PetscScalar v[])
2105 {
2106 Mat_HYPRE *hA = (Mat_HYPRE *)A->data;
2107 PetscInt i;
2108
2109 PetscFunctionBegin;
2110 if (!m || !n) PetscFunctionReturn(PETSC_SUCCESS);
2111 /* Ignore negative row indices
2112 * And negative column indices should be automatically ignored in hypre
2113 * */
2114 for (i = 0; i < m; i++) {
2115 if (idxm[i] >= 0) {
2116 HYPRE_Int hn = (HYPRE_Int)n;
2117 PetscCallHYPRE(HYPRE_IJMatrixGetValues(hA->ij, 1, &hn, (HYPRE_BigInt *)&idxm[i], (HYPRE_BigInt *)idxn, (HYPRE_Complex *)(v + i * n)));
2118 }
2119 }
2120 PetscFunctionReturn(PETSC_SUCCESS);
2121 }
2122
2123 static PetscErrorCode MatSetOption_HYPRE(Mat A, MatOption op, PetscBool flg)
2124 {
2125 Mat_HYPRE *hA = (Mat_HYPRE *)A->data;
2126
2127 PetscFunctionBegin;
2128 switch (op) {
2129 case MAT_NO_OFF_PROC_ENTRIES:
2130 if (flg) PetscCallHYPRE(HYPRE_IJMatrixSetMaxOffProcElmts(hA->ij, 0));
2131 break;
2132 case MAT_IGNORE_OFF_PROC_ENTRIES:
2133 hA->donotstash = flg;
2134 break;
2135 default:
2136 break;
2137 }
2138 PetscFunctionReturn(PETSC_SUCCESS);
2139 }
2140
2141 static PetscErrorCode MatView_HYPRE(Mat A, PetscViewer view)
2142 {
2143 PetscViewerFormat format;
2144
2145 PetscFunctionBegin;
2146 PetscCall(PetscViewerGetFormat(view, &format));
2147 if (format == PETSC_VIEWER_ASCII_FACTOR_INFO || format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) PetscFunctionReturn(PETSC_SUCCESS);
2148 if (format != PETSC_VIEWER_NATIVE) {
2149 Mat B;
2150 hypre_ParCSRMatrix *parcsr;
2151 PetscErrorCode (*mview)(Mat, PetscViewer) = NULL;
2152
2153 PetscCall(MatHYPREGetParCSR_HYPRE(A, &parcsr));
2154 PetscCall(MatCreateFromParCSR(parcsr, MATAIJ, PETSC_USE_POINTER, &B));
2155 PetscCall(MatGetOperation(B, MATOP_VIEW, (PetscErrorCodeFn **)&mview));
2156 PetscCheck(mview, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing view operation");
2157 PetscCall((*mview)(B, view));
2158 PetscCall(MatDestroy(&B));
2159 } else {
2160 Mat_HYPRE *hA = (Mat_HYPRE *)A->data;
2161 PetscMPIInt size;
2162 PetscBool isascii;
2163 const char *filename;
2164
2165 /* HYPRE uses only text files */
2166 PetscCall(PetscObjectTypeCompare((PetscObject)view, PETSCVIEWERASCII, &isascii));
2167 PetscCheck(isascii, PetscObjectComm((PetscObject)view), PETSC_ERR_SUP, "PetscViewerType %s: native HYPRE format needs PETSCVIEWERASCII", ((PetscObject)view)->type_name);
2168 PetscCall(PetscViewerFileGetName(view, &filename));
2169 PetscCallHYPRE(HYPRE_IJMatrixPrint(hA->ij, filename));
2170 PetscCallMPI(MPI_Comm_size(hA->comm, &size));
2171 if (size > 1) {
2172 PetscCall(PetscViewerASCIIPrintf(view, "Matrix files: %s.%05d ... %s.%05d\n", filename, 0, filename, size - 1));
2173 } else {
2174 PetscCall(PetscViewerASCIIPrintf(view, "Matrix file: %s.%05d\n", filename, 0));
2175 }
2176 }
2177 PetscFunctionReturn(PETSC_SUCCESS);
2178 }
2179
2180 static PetscErrorCode MatCopy_HYPRE(Mat A, Mat B, MatStructure str)
2181 {
2182 hypre_ParCSRMatrix *acsr, *bcsr;
2183
2184 PetscFunctionBegin;
2185 if (str == SAME_NONZERO_PATTERN && A->ops->copy == B->ops->copy) {
2186 PetscCall(MatHYPREGetParCSR_HYPRE(A, &acsr));
2187 PetscCall(MatHYPREGetParCSR_HYPRE(B, &bcsr));
2188 PetscCallHYPRE(hypre_ParCSRMatrixCopy(acsr, bcsr, 1));
2189 PetscCall(MatSetOption(B, MAT_SORTED_FULL, PETSC_TRUE)); /* "perfect" preallocation, so no need for hypre_AuxParCSRMatrixNeedAux */
2190 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
2191 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
2192 } else {
2193 PetscCall(MatCopy_Basic(A, B, str));
2194 }
2195 PetscFunctionReturn(PETSC_SUCCESS);
2196 }
2197
2198 static PetscErrorCode MatGetDiagonal_HYPRE(Mat A, Vec d)
2199 {
2200 hypre_ParCSRMatrix *parcsr;
2201 hypre_CSRMatrix *dmat;
2202 HYPRE_Complex *a;
2203 PetscBool cong;
2204
2205 PetscFunctionBegin;
2206 PetscCall(MatHasCongruentLayouts(A, &cong));
2207 PetscCheck(cong, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Only for square matrices with same local distributions of rows and columns");
2208 PetscCall(MatHYPREGetParCSR_HYPRE(A, &parcsr));
2209 dmat = hypre_ParCSRMatrixDiag(parcsr);
2210 if (dmat) {
2211 #if defined(PETSC_HAVE_HYPRE_DEVICE)
2212 HYPRE_MemoryLocation mem = hypre_CSRMatrixMemoryLocation(dmat);
2213 #else
2214 HYPRE_MemoryLocation mem = HYPRE_MEMORY_HOST;
2215 #endif
2216
2217 if (mem != HYPRE_MEMORY_HOST) PetscCall(VecGetArrayWriteAndMemType(d, (PetscScalar **)&a, NULL));
2218 else PetscCall(VecGetArrayWrite(d, (PetscScalar **)&a));
2219 hypre_CSRMatrixExtractDiagonal(dmat, a, 0);
2220 if (mem != HYPRE_MEMORY_HOST) PetscCall(VecRestoreArrayWriteAndMemType(d, (PetscScalar **)&a));
2221 else PetscCall(VecRestoreArrayWrite(d, (PetscScalar **)&a));
2222 }
2223 PetscFunctionReturn(PETSC_SUCCESS);
2224 }
2225
2226 #include <petscblaslapack.h>
2227
2228 static PetscErrorCode MatAXPY_HYPRE(Mat Y, PetscScalar a, Mat X, MatStructure str)
2229 {
2230 PetscFunctionBegin;
2231 #if defined(PETSC_HAVE_HYPRE_DEVICE)
2232 {
2233 Mat B;
2234 hypre_ParCSRMatrix *x, *y, *z;
2235
2236 PetscCall(MatHYPREGetParCSR(Y, &y));
2237 PetscCall(MatHYPREGetParCSR(X, &x));
2238 PetscCallHYPRE(hypre_ParCSRMatrixAdd(1.0, y, 1.0, x, &z));
2239 PetscCall(MatCreateFromParCSR(z, MATHYPRE, PETSC_OWN_POINTER, &B));
2240 PetscCall(MatHeaderMerge(Y, &B));
2241 }
2242 #else
2243 if (str == SAME_NONZERO_PATTERN) {
2244 hypre_ParCSRMatrix *x, *y;
2245 hypre_CSRMatrix *xloc, *yloc;
2246 PetscInt xnnz, ynnz;
2247 HYPRE_Complex *xarr, *yarr;
2248 PetscBLASInt one = 1, bnz;
2249
2250 PetscCall(MatHYPREGetParCSR(Y, &y));
2251 PetscCall(MatHYPREGetParCSR(X, &x));
2252
2253 /* diagonal block */
2254 xloc = hypre_ParCSRMatrixDiag(x);
2255 yloc = hypre_ParCSRMatrixDiag(y);
2256 xnnz = 0;
2257 ynnz = 0;
2258 xarr = NULL;
2259 yarr = NULL;
2260 if (xloc) {
2261 xarr = hypre_CSRMatrixData(xloc);
2262 xnnz = hypre_CSRMatrixNumNonzeros(xloc);
2263 }
2264 if (yloc) {
2265 yarr = hypre_CSRMatrixData(yloc);
2266 ynnz = hypre_CSRMatrixNumNonzeros(yloc);
2267 }
2268 PetscCheck(xnnz == ynnz, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Different number of nonzeros in diagonal block %" PetscInt_FMT " != %" PetscInt_FMT, xnnz, ynnz);
2269 PetscCall(PetscBLASIntCast(xnnz, &bnz));
2270 PetscCallBLAS("BLASaxpy", BLASaxpy_(&bnz, &a, (PetscScalar *)xarr, &one, (PetscScalar *)yarr, &one));
2271
2272 /* off-diagonal block */
2273 xloc = hypre_ParCSRMatrixOffd(x);
2274 yloc = hypre_ParCSRMatrixOffd(y);
2275 xnnz = 0;
2276 ynnz = 0;
2277 xarr = NULL;
2278 yarr = NULL;
2279 if (xloc) {
2280 xarr = hypre_CSRMatrixData(xloc);
2281 xnnz = hypre_CSRMatrixNumNonzeros(xloc);
2282 }
2283 if (yloc) {
2284 yarr = hypre_CSRMatrixData(yloc);
2285 ynnz = hypre_CSRMatrixNumNonzeros(yloc);
2286 }
2287 PetscCheck(xnnz == ynnz, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Different number of nonzeros in off-diagonal block %" PetscInt_FMT " != %" PetscInt_FMT, xnnz, ynnz);
2288 PetscCall(PetscBLASIntCast(xnnz, &bnz));
2289 PetscCallBLAS("BLASaxpy", BLASaxpy_(&bnz, &a, (PetscScalar *)xarr, &one, (PetscScalar *)yarr, &one));
2290 } else if (str == SUBSET_NONZERO_PATTERN) {
2291 PetscCall(MatAXPY_Basic(Y, a, X, str));
2292 } else {
2293 Mat B;
2294
2295 PetscCall(MatAXPY_Basic_Preallocate(Y, X, &B));
2296 PetscCall(MatAXPY_BasicWithPreallocation(B, Y, a, X, str));
2297 PetscCall(MatHeaderReplace(Y, &B));
2298 }
2299 #endif
2300 PetscFunctionReturn(PETSC_SUCCESS);
2301 }
2302
2303 static PetscErrorCode MatDuplicate_HYPRE(Mat A, MatDuplicateOption op, Mat *B)
2304 {
2305 hypre_ParCSRMatrix *parcsr = NULL;
2306 PetscCopyMode cpmode;
2307 Mat_HYPRE *hA;
2308
2309 PetscFunctionBegin;
2310 PetscCall(MatHYPREGetParCSR_HYPRE(A, &parcsr));
2311 if (op == MAT_DO_NOT_COPY_VALUES || op == MAT_SHARE_NONZERO_PATTERN) {
2312 parcsr = hypre_ParCSRMatrixClone(parcsr, 0);
2313 cpmode = PETSC_OWN_POINTER;
2314 } else {
2315 cpmode = PETSC_COPY_VALUES;
2316 }
2317 PetscCall(MatCreateFromParCSR(parcsr, MATHYPRE, cpmode, B));
2318 hA = (Mat_HYPRE *)A->data;
2319 if (hA->cooMat) {
2320 Mat_HYPRE *hB = (Mat_HYPRE *)((*B)->data);
2321 op = (op == MAT_DO_NOT_COPY_VALUES) ? op : MAT_COPY_VALUES;
2322 /* Cannot simply increase the reference count of hA->cooMat, since B needs to share cooMat's data array */
2323 PetscCall(MatDuplicate(hA->cooMat, op, &hB->cooMat));
2324 PetscCall(MatHYPRE_AttachCOOMat(*B));
2325 }
2326 PetscFunctionReturn(PETSC_SUCCESS);
2327 }
2328
2329 static PetscErrorCode MatSetPreallocationCOO_HYPRE(Mat mat, PetscCount coo_n, PetscInt coo_i[], PetscInt coo_j[])
2330 {
2331 Mat_HYPRE *hmat = (Mat_HYPRE *)mat->data;
2332
2333 PetscFunctionBegin;
2334 /* Build an agent matrix cooMat with AIJ format
2335 It has the same sparsity pattern as mat, and also shares the data array with mat. We use cooMat to do the COO work.
2336 */
2337 PetscCall(MatHYPRE_CreateCOOMat(mat));
2338 PetscCall(MatSetOption(hmat->cooMat, MAT_IGNORE_OFF_PROC_ENTRIES, hmat->donotstash));
2339 PetscCall(MatSetOption(hmat->cooMat, MAT_NO_OFF_PROC_ENTRIES, mat->nooffprocentries));
2340
2341 /* MatSetPreallocationCOO_SeqAIJ and MatSetPreallocationCOO_MPIAIJ uses this specific
2342 name to automatically put the diagonal entries first */
2343 PetscCall(PetscObjectSetName((PetscObject)hmat->cooMat, "_internal_COO_mat_for_hypre"));
2344 PetscCall(MatSetPreallocationCOO(hmat->cooMat, coo_n, coo_i, coo_j));
2345 hmat->cooMat->assembled = PETSC_TRUE;
2346
2347 /* Copy the sparsity pattern from cooMat to hypre IJMatrix hmat->ij */
2348 PetscCall(MatSetOption(mat, MAT_SORTED_FULL, PETSC_TRUE));
2349 PetscCall(MatHYPRE_CreateFromMat(hmat->cooMat, hmat)); /* Create hmat->ij and preallocate it */
2350 PetscCall(MatHYPRE_IJMatrixCopyIJ(hmat->cooMat, hmat->ij)); /* Copy A's (i,j) to hmat->ij */
2351
2352 mat->preallocated = PETSC_TRUE;
2353 PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY));
2354 PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY)); /* Migrate mat to device if it is bound to. Hypre builds its own SpMV context here */
2355
2356 /* Attach cooMat to mat */
2357 PetscCall(MatHYPRE_AttachCOOMat(mat));
2358 PetscFunctionReturn(PETSC_SUCCESS);
2359 }
2360
2361 static PetscErrorCode MatSetValuesCOO_HYPRE(Mat mat, const PetscScalar v[], InsertMode imode)
2362 {
2363 Mat_HYPRE *hmat = (Mat_HYPRE *)mat->data;
2364
2365 PetscFunctionBegin;
2366 PetscCheck(hmat->cooMat, PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "HYPRE COO delegate matrix has not been created yet");
2367 PetscCall(MatSetValuesCOO(hmat->cooMat, v, imode));
2368 PetscCall(MatViewFromOptions(hmat->cooMat, (PetscObject)mat, "-cooMat_view"));
2369 PetscFunctionReturn(PETSC_SUCCESS);
2370 }
2371
2372 static PetscErrorCode MatGetCurrentMemType_HYPRE(Mat A, PetscMemType *m)
2373 {
2374 PetscBool petsconcpu;
2375
2376 PetscFunctionBegin;
2377 PetscCall(MatBoundToCPU(A, &petsconcpu));
2378 *m = petsconcpu ? PETSC_MEMTYPE_HOST : PETSC_MEMTYPE_DEVICE;
2379 PetscFunctionReturn(PETSC_SUCCESS);
2380 }
2381
2382 /*MC
2383 MATHYPRE - "hypre" - A matrix type to be used for sequential and parallel sparse matrices
2384 based on the hypre IJ interface.
2385
2386 Level: intermediate
2387
2388 .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MatHYPRESetPreallocation`
2389 M*/
2390 PETSC_EXTERN PetscErrorCode MatCreate_HYPRE(Mat B)
2391 {
2392 Mat_HYPRE *hB;
2393 #if defined(PETSC_HAVE_HYPRE_DEVICE)
2394 HYPRE_MemoryLocation memory_location;
2395 #endif
2396
2397 PetscFunctionBegin;
2398 PetscCall(PetscHYPREInitialize());
2399 PetscCall(PetscNew(&hB));
2400
2401 hB->inner_free = PETSC_TRUE;
2402 hB->array_available = PETSC_TRUE;
2403
2404 B->data = (void *)hB;
2405
2406 PetscCall(PetscMemzero(B->ops, sizeof(struct _MatOps)));
2407 B->ops->mult = MatMult_HYPRE;
2408 B->ops->multtranspose = MatMultTranspose_HYPRE;
2409 B->ops->multadd = MatMultAdd_HYPRE;
2410 B->ops->multtransposeadd = MatMultTransposeAdd_HYPRE;
2411 B->ops->setup = MatSetUp_HYPRE;
2412 B->ops->destroy = MatDestroy_HYPRE;
2413 B->ops->assemblyend = MatAssemblyEnd_HYPRE;
2414 B->ops->assemblybegin = MatAssemblyBegin_HYPRE;
2415 B->ops->setvalues = MatSetValues_HYPRE;
2416 B->ops->scale = MatScale_HYPRE;
2417 B->ops->zerorowscolumns = MatZeroRowsColumns_HYPRE;
2418 B->ops->zeroentries = MatZeroEntries_HYPRE;
2419 B->ops->zerorows = MatZeroRows_HYPRE;
2420 B->ops->getrow = MatGetRow_HYPRE;
2421 B->ops->restorerow = MatRestoreRow_HYPRE;
2422 B->ops->getvalues = MatGetValues_HYPRE;
2423 B->ops->setoption = MatSetOption_HYPRE;
2424 B->ops->duplicate = MatDuplicate_HYPRE;
2425 B->ops->copy = MatCopy_HYPRE;
2426 B->ops->view = MatView_HYPRE;
2427 B->ops->getdiagonal = MatGetDiagonal_HYPRE;
2428 B->ops->axpy = MatAXPY_HYPRE;
2429 B->ops->productsetfromoptions = MatProductSetFromOptions_HYPRE;
2430 B->ops->getcurrentmemtype = MatGetCurrentMemType_HYPRE;
2431 #if defined(PETSC_HAVE_HYPRE_DEVICE)
2432 B->ops->bindtocpu = MatBindToCPU_HYPRE;
2433 /* Get hypre's default memory location. Users can control this using the corresponding HYPRE_SetMemoryLocation API */
2434 PetscCallHYPRE(HYPRE_GetMemoryLocation(&memory_location));
2435 B->boundtocpu = (memory_location == HYPRE_MEMORY_HOST) ? PETSC_TRUE : PETSC_FALSE;
2436 #endif
2437
2438 /* build cache for off array entries formed */
2439 PetscCall(MatStashCreate_Private(PetscObjectComm((PetscObject)B), 1, &B->stash));
2440
2441 PetscCall(PetscCommGetComm(PetscObjectComm((PetscObject)B), &hB->comm));
2442 PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATHYPRE));
2443 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_hypre_aij_C", MatConvert_HYPRE_AIJ));
2444 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_hypre_is_C", MatConvert_HYPRE_IS));
2445 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqaij_hypre_C", MatProductSetFromOptions_HYPRE));
2446 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_mpiaij_hypre_C", MatProductSetFromOptions_HYPRE));
2447 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatHYPRESetPreallocation_C", MatHYPRESetPreallocation_HYPRE));
2448 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatHYPREGetParCSR_C", MatHYPREGetParCSR_HYPRE));
2449 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSetPreallocationCOO_C", MatSetPreallocationCOO_HYPRE));
2450 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSetValuesCOO_C", MatSetValuesCOO_HYPRE));
2451 #if defined(PETSC_HAVE_HYPRE_DEVICE)
2452 #if defined(HYPRE_USING_HIP)
2453 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqaijhipsparse_hypre_C", MatProductSetFromOptions_HYPRE));
2454 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_mpiaijhipsparse_hypre_C", MatProductSetFromOptions_HYPRE));
2455 PetscCall(PetscDeviceInitialize(PETSC_DEVICE_HIP));
2456 PetscCall(MatSetVecType(B, VECHIP));
2457 #endif
2458 #if defined(HYPRE_USING_CUDA)
2459 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqaijcusparse_hypre_C", MatProductSetFromOptions_HYPRE));
2460 PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_mpiaijcusparse_hypre_C", MatProductSetFromOptions_HYPRE));
2461 PetscCall(PetscDeviceInitialize(PETSC_DEVICE_CUDA));
2462 PetscCall(MatSetVecType(B, VECCUDA));
2463 #endif
2464 #endif
2465 PetscFunctionReturn(PETSC_SUCCESS);
2466 }
2467