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