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: }