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