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