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