Actual source code: mhypre.c

  1: /*
  2:     Creates hypre ijmatrix from PETSc matrix
  3: */

  5: #include <petscpkg_version.h>
  6: #include <petsc/private/petschypre.h>
  7: #include <petscmathypre.h>
  8: #include <petsc/private/matimpl.h>
  9: #include <petsc/private/deviceimpl.h>
 10: #include <../src/mat/impls/hypre/mhypre.h>
 11: #include <../src/mat/impls/aij/mpi/mpiaij.h>
 12: #include <../src/vec/vec/impls/hypre/vhyp.h>
 13: #include <HYPRE.h>
 14: #include <HYPRE_utilities.h>
 15: #include <_hypre_parcsr_ls.h>
 16: #include <_hypre_sstruct_ls.h>

 18: #if PETSC_PKG_HYPRE_VERSION_LT(2, 18, 0)
 19:   #define hypre_ParCSRMatrixClone(A, B) hypre_ParCSRMatrixCompleteClone(A)
 20: #endif

 22: #if PETSC_PKG_HYPRE_VERSION_GE(2, 15, 0)
 23:   #define HYPRE_AssumedPartitionCheck() 1
 24: #endif

 26: static PetscErrorCode MatHYPRE_CreateFromMat(Mat, Mat_HYPRE *);
 27: static PetscErrorCode MatHYPRE_IJMatrixPreallocate(Mat, Mat, HYPRE_IJMatrix);
 28: static PetscErrorCode MatHYPRE_IJMatrixCopyIJ_MPIAIJ(Mat, HYPRE_IJMatrix);
 29: static PetscErrorCode MatHYPRE_IJMatrixCopyIJ_SeqAIJ(Mat, HYPRE_IJMatrix);
 30: static PetscErrorCode MatHYPRE_MultKernel_Private(Mat, HYPRE_Complex, Vec, HYPRE_Complex, Vec, PetscBool);
 31: static PetscErrorCode MatSetValues_HYPRE(Mat, PetscInt, const PetscInt[], PetscInt, const PetscInt[], const PetscScalar[], InsertMode ins);

 33: static PetscErrorCode MatHYPRE_IJMatrixPreallocate(Mat A_d, Mat A_o, HYPRE_IJMatrix ij)
 34: {
 35:   PetscInt        i, n_d, n_o;
 36:   const PetscInt *ia_d, *ia_o;
 37:   PetscBool       done_d = PETSC_FALSE, done_o = PETSC_FALSE;
 38:   HYPRE_Int      *nnz_d = NULL, *nnz_o = NULL;

 40:   PetscFunctionBegin;
 41:   if (A_d) { /* determine number of nonzero entries in local diagonal part */
 42:     PetscCall(MatGetRowIJ(A_d, 0, PETSC_FALSE, PETSC_FALSE, &n_d, &ia_d, NULL, &done_d));
 43:     if (done_d) {
 44:       PetscCall(PetscMalloc1(n_d, &nnz_d));
 45:       for (i = 0; i < n_d; i++) nnz_d[i] = ia_d[i + 1] - ia_d[i];
 46:     }
 47:     PetscCall(MatRestoreRowIJ(A_d, 0, PETSC_FALSE, PETSC_FALSE, NULL, &ia_d, NULL, &done_d));
 48:   }
 49:   if (A_o) { /* determine number of nonzero entries in local off-diagonal part */
 50:     PetscCall(MatGetRowIJ(A_o, 0, PETSC_FALSE, PETSC_FALSE, &n_o, &ia_o, NULL, &done_o));
 51:     if (done_o) {
 52:       PetscCall(PetscMalloc1(n_o, &nnz_o));
 53:       for (i = 0; i < n_o; i++) nnz_o[i] = ia_o[i + 1] - ia_o[i];
 54:     }
 55:     PetscCall(MatRestoreRowIJ(A_o, 0, PETSC_FALSE, PETSC_FALSE, &n_o, &ia_o, NULL, &done_o));
 56:   }
 57:   if (done_d) {    /* set number of nonzeros in HYPRE IJ matrix */
 58:     if (!done_o) { /* only diagonal part */
 59:       PetscCall(PetscCalloc1(n_d, &nnz_o));
 60:     }
 61: #if PETSC_PKG_HYPRE_VERSION_GE(2, 16, 0)
 62:     { /* If we don't do this, the columns of the matrix will be all zeros! */
 63:       hypre_AuxParCSRMatrix *aux_matrix;
 64:       aux_matrix = (hypre_AuxParCSRMatrix *)hypre_IJMatrixTranslator(ij);
 65:       hypre_AuxParCSRMatrixDestroy(aux_matrix);
 66:       hypre_IJMatrixTranslator(ij) = NULL;
 67:       PetscCallExternal(HYPRE_IJMatrixSetDiagOffdSizes, ij, nnz_d, nnz_o);
 68:       /* it seems they partially fixed it in 2.19.0 */
 69:   #if PETSC_PKG_HYPRE_VERSION_LT(2, 19, 0)
 70:       aux_matrix                               = (hypre_AuxParCSRMatrix *)hypre_IJMatrixTranslator(ij);
 71:       hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 1;
 72:   #endif
 73:     }
 74: #else
 75:     PetscCallExternal(HYPRE_IJMatrixSetDiagOffdSizes, ij, nnz_d, nnz_o);
 76: #endif
 77:     PetscCall(PetscFree(nnz_d));
 78:     PetscCall(PetscFree(nnz_o));
 79:   }
 80:   PetscFunctionReturn(PETSC_SUCCESS);
 81: }

 83: static PetscErrorCode MatHYPRE_CreateFromMat(Mat A, Mat_HYPRE *hA)
 84: {
 85:   PetscInt rstart, rend, cstart, cend;

 87:   PetscFunctionBegin;
 88:   PetscCall(PetscLayoutSetUp(A->rmap));
 89:   PetscCall(PetscLayoutSetUp(A->cmap));
 90:   rstart = A->rmap->rstart;
 91:   rend   = A->rmap->rend;
 92:   cstart = A->cmap->rstart;
 93:   cend   = A->cmap->rend;
 94:   PetscHYPREInitialize();
 95:   if (hA->ij) {
 96:     if (!hA->inner_free) hypre_IJMatrixObject(hA->ij) = NULL;
 97:     PetscCallExternal(HYPRE_IJMatrixDestroy, hA->ij);
 98:   }
 99:   PetscCallExternal(HYPRE_IJMatrixCreate, hA->comm, rstart, rend - 1, cstart, cend - 1, &hA->ij);
100:   PetscCallExternal(HYPRE_IJMatrixSetObjectType, hA->ij, HYPRE_PARCSR);
101:   {
102:     PetscBool       same;
103:     Mat             A_d, A_o;
104:     const PetscInt *colmap;
105:     PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIAIJ, &same));
106:     if (same) {
107:       PetscCall(MatMPIAIJGetSeqAIJ(A, &A_d, &A_o, &colmap));
108:       PetscCall(MatHYPRE_IJMatrixPreallocate(A_d, A_o, hA->ij));
109:       PetscFunctionReturn(PETSC_SUCCESS);
110:     }
111:     PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIBAIJ, &same));
112:     if (same) {
113:       PetscCall(MatMPIBAIJGetSeqBAIJ(A, &A_d, &A_o, &colmap));
114:       PetscCall(MatHYPRE_IJMatrixPreallocate(A_d, A_o, hA->ij));
115:       PetscFunctionReturn(PETSC_SUCCESS);
116:     }
117:     PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATSEQAIJ, &same));
118:     if (same) {
119:       PetscCall(MatHYPRE_IJMatrixPreallocate(A, NULL, hA->ij));
120:       PetscFunctionReturn(PETSC_SUCCESS);
121:     }
122:     PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATSEQBAIJ, &same));
123:     if (same) {
124:       PetscCall(MatHYPRE_IJMatrixPreallocate(A, NULL, hA->ij));
125:       PetscFunctionReturn(PETSC_SUCCESS);
126:     }
127:   }
128:   PetscFunctionReturn(PETSC_SUCCESS);
129: }

131: static PetscErrorCode MatHYPRE_IJMatrixCopyIJ(Mat A, HYPRE_IJMatrix ij)
132: {
133:   PetscBool flg;

135:   PetscFunctionBegin;
136: #if PETSC_PKG_HYPRE_VERSION_LT(2, 19, 0)
137:   PetscCallExternal(HYPRE_IJMatrixInitialize, ij);
138: #else
139:   PetscCallExternal(HYPRE_IJMatrixInitialize_v2, ij, HYPRE_MEMORY_HOST);
140: #endif
141:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIAIJ, &flg));
142:   if (flg) {
143:     PetscCall(MatHYPRE_IJMatrixCopyIJ_MPIAIJ(A, ij));
144:     PetscFunctionReturn(PETSC_SUCCESS);
145:   }
146:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATSEQAIJ, &flg));
147:   if (flg) {
148:     PetscCall(MatHYPRE_IJMatrixCopyIJ_SeqAIJ(A, ij));
149:     PetscFunctionReturn(PETSC_SUCCESS);
150:   }
151:   PetscCheck(PETSC_FALSE, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "No support for matrix type %s", ((PetscObject)A)->type_name);
152:   PetscFunctionReturn(PETSC_SUCCESS);
153: }

155: static PetscErrorCode MatHYPRE_IJMatrixCopyIJ_SeqAIJ(Mat A, HYPRE_IJMatrix ij)
156: {
157:   Mat_SeqAIJ            *pdiag = (Mat_SeqAIJ *)A->data;
158:   HYPRE_Int              type;
159:   hypre_ParCSRMatrix    *par_matrix;
160:   hypre_AuxParCSRMatrix *aux_matrix;
161:   hypre_CSRMatrix       *hdiag;
162:   PetscBool              sameint = (PetscBool)(sizeof(PetscInt) == sizeof(HYPRE_Int));

164:   PetscFunctionBegin;
165:   PetscCallExternal(HYPRE_IJMatrixGetObjectType, ij, &type);
166:   PetscCheck(type == HYPRE_PARCSR, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Only HYPRE_PARCSR is supported");
167:   PetscCallExternal(HYPRE_IJMatrixGetObject, ij, (void **)&par_matrix);
168:   hdiag = hypre_ParCSRMatrixDiag(par_matrix);
169:   /*
170:        this is the Hack part where we monkey directly with the hypre datastructures
171:   */
172:   if (sameint) {
173:     PetscCall(PetscArraycpy(hdiag->i, pdiag->i, A->rmap->n + 1));
174:     PetscCall(PetscArraycpy(hdiag->j, pdiag->j, pdiag->nz));
175:   } else {
176:     PetscInt i;

178:     for (i = 0; i < A->rmap->n + 1; i++) hdiag->i[i] = (HYPRE_Int)pdiag->i[i];
179:     for (i = 0; i < pdiag->nz; i++) hdiag->j[i] = (HYPRE_Int)pdiag->j[i];
180:   }

182:   aux_matrix                               = (hypre_AuxParCSRMatrix *)hypre_IJMatrixTranslator(ij);
183:   hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 0;
184:   PetscFunctionReturn(PETSC_SUCCESS);
185: }

187: static PetscErrorCode MatHYPRE_IJMatrixCopyIJ_MPIAIJ(Mat A, HYPRE_IJMatrix ij)
188: {
189:   Mat_MPIAIJ            *pA = (Mat_MPIAIJ *)A->data;
190:   Mat_SeqAIJ            *pdiag, *poffd;
191:   PetscInt               i, *garray = pA->garray, *jj, cstart, *pjj;
192:   HYPRE_Int             *hjj, type;
193:   hypre_ParCSRMatrix    *par_matrix;
194:   hypre_AuxParCSRMatrix *aux_matrix;
195:   hypre_CSRMatrix       *hdiag, *hoffd;
196:   PetscBool              sameint = (PetscBool)(sizeof(PetscInt) == sizeof(HYPRE_Int));

198:   PetscFunctionBegin;
199:   pdiag = (Mat_SeqAIJ *)pA->A->data;
200:   poffd = (Mat_SeqAIJ *)pA->B->data;
201:   /* cstart is only valid for square MPIAIJ laid out in the usual way */
202:   PetscCall(MatGetOwnershipRange(A, &cstart, NULL));

204:   PetscCallExternal(HYPRE_IJMatrixGetObjectType, ij, &type);
205:   PetscCheck(type == HYPRE_PARCSR, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Only HYPRE_PARCSR is supported");
206:   PetscCallExternal(HYPRE_IJMatrixGetObject, ij, (void **)&par_matrix);
207:   hdiag = hypre_ParCSRMatrixDiag(par_matrix);
208:   hoffd = hypre_ParCSRMatrixOffd(par_matrix);

210:   if (sameint) {
211:     PetscCall(PetscArraycpy(hdiag->i, pdiag->i, pA->A->rmap->n + 1));
212:   } else {
213:     for (i = 0; i < pA->A->rmap->n + 1; i++) hdiag->i[i] = (HYPRE_Int)pdiag->i[i];
214:   }

216:   hjj = hdiag->j;
217:   pjj = pdiag->j;
218: #if PETSC_PKG_HYPRE_VERSION_GE(2, 16, 0)
219:   for (i = 0; i < pdiag->nz; i++) hjj[i] = pjj[i];
220: #else
221:   for (i = 0; i < pdiag->nz; i++) hjj[i] = cstart + pjj[i];
222: #endif
223:   if (sameint) {
224:     PetscCall(PetscArraycpy(hoffd->i, poffd->i, pA->A->rmap->n + 1));
225:   } else {
226:     for (i = 0; i < pA->A->rmap->n + 1; i++) hoffd->i[i] = (HYPRE_Int)poffd->i[i];
227:   }

229:   jj = (PetscInt *)hoffd->j;
230: #if PETSC_PKG_HYPRE_VERSION_GE(2, 16, 0)
231:   PetscCallExternal(hypre_CSRMatrixBigInitialize, hoffd);
232:   jj = (PetscInt *)hoffd->big_j;
233: #endif
234:   pjj = poffd->j;
235:   for (i = 0; i < poffd->nz; i++) jj[i] = garray[pjj[i]];

237:   aux_matrix                               = (hypre_AuxParCSRMatrix *)hypre_IJMatrixTranslator(ij);
238:   hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 0;
239:   PetscFunctionReturn(PETSC_SUCCESS);
240: }

242: static PetscErrorCode MatConvert_HYPRE_IS(Mat A, MatType mtype, MatReuse reuse, Mat *B)
243: {
244:   Mat_HYPRE             *mhA = (Mat_HYPRE *)A->data;
245:   Mat                    lA;
246:   ISLocalToGlobalMapping rl2g, cl2g;
247:   IS                     is;
248:   hypre_ParCSRMatrix    *hA;
249:   hypre_CSRMatrix       *hdiag, *hoffd;
250:   MPI_Comm               comm;
251:   HYPRE_Complex         *hdd, *hod, *aa;
252:   PetscScalar           *data;
253:   HYPRE_BigInt          *col_map_offd;
254:   HYPRE_Int             *hdi, *hdj, *hoi, *hoj;
255:   PetscInt              *ii, *jj, *iptr, *jptr;
256:   PetscInt               cum, dr, dc, oc, str, stc, nnz, i, jd, jo, M, N;
257:   HYPRE_Int              type;
258:   MatType                lmattype   = NULL;
259:   PetscBool              freeparcsr = PETSC_FALSE;

261:   PetscFunctionBegin;
262:   comm = PetscObjectComm((PetscObject)A);
263:   PetscCallExternal(HYPRE_IJMatrixGetObjectType, mhA->ij, &type);
264:   PetscCheck(type == HYPRE_PARCSR, comm, PETSC_ERR_SUP, "Only HYPRE_PARCSR is supported");
265:   PetscCallExternal(HYPRE_IJMatrixGetObject, mhA->ij, (void **)&hA);
266: #if defined(PETSC_HAVE_HYPRE_DEVICE)
267:   if (HYPRE_MEMORY_DEVICE == hypre_IJMatrixMemoryLocation(mhA->ij)) {
268:     /* Support by copying back on the host and copy to GPU
269:        Kind of inefficient, but this is the best we can do now */
270:   #if defined(HYPRE_USING_HIP)
271:     lmattype = MATSEQAIJHIPSPARSE;
272:   #elif defined(HYPRE_USING_CUDA)
273:     lmattype = MATSEQAIJCUSPARSE;
274:   #endif
275:     hA         = hypre_ParCSRMatrixClone_v2(hA, 1, HYPRE_MEMORY_HOST);
276:     freeparcsr = PETSC_TRUE;
277:   }
278: #endif
279:   M     = hypre_ParCSRMatrixGlobalNumRows(hA);
280:   N     = hypre_ParCSRMatrixGlobalNumCols(hA);
281:   str   = hypre_ParCSRMatrixFirstRowIndex(hA);
282:   stc   = hypre_ParCSRMatrixFirstColDiag(hA);
283:   hdiag = hypre_ParCSRMatrixDiag(hA);
284:   hoffd = hypre_ParCSRMatrixOffd(hA);
285:   dr    = hypre_CSRMatrixNumRows(hdiag);
286:   dc    = hypre_CSRMatrixNumCols(hdiag);
287:   nnz   = hypre_CSRMatrixNumNonzeros(hdiag);
288:   hdi   = hypre_CSRMatrixI(hdiag);
289:   hdj   = hypre_CSRMatrixJ(hdiag);
290:   hdd   = hypre_CSRMatrixData(hdiag);
291:   oc    = hypre_CSRMatrixNumCols(hoffd);
292:   nnz += hypre_CSRMatrixNumNonzeros(hoffd);
293:   hoi = hypre_CSRMatrixI(hoffd);
294:   hoj = hypre_CSRMatrixJ(hoffd);
295:   hod = hypre_CSRMatrixData(hoffd);
296:   if (reuse != MAT_REUSE_MATRIX) {
297:     PetscInt *aux;

299:     /* generate l2g maps for rows and cols */
300:     PetscCall(ISCreateStride(comm, dr, str, 1, &is));
301:     PetscCall(ISLocalToGlobalMappingCreateIS(is, &rl2g));
302:     PetscCall(ISDestroy(&is));
303:     col_map_offd = hypre_ParCSRMatrixColMapOffd(hA);
304:     PetscCall(PetscMalloc1(dc + oc, &aux));
305:     for (i = 0; i < dc; i++) aux[i] = i + stc;
306:     for (i = 0; i < oc; i++) aux[i + dc] = col_map_offd[i];
307:     PetscCall(ISCreateGeneral(comm, dc + oc, aux, PETSC_OWN_POINTER, &is));
308:     PetscCall(ISLocalToGlobalMappingCreateIS(is, &cl2g));
309:     PetscCall(ISDestroy(&is));
310:     /* create MATIS object */
311:     PetscCall(MatCreate(comm, B));
312:     PetscCall(MatSetSizes(*B, dr, dc, M, N));
313:     PetscCall(MatSetType(*B, MATIS));
314:     PetscCall(MatSetLocalToGlobalMapping(*B, rl2g, cl2g));
315:     PetscCall(ISLocalToGlobalMappingDestroy(&rl2g));
316:     PetscCall(ISLocalToGlobalMappingDestroy(&cl2g));

318:     /* allocate CSR for local matrix */
319:     PetscCall(PetscMalloc1(dr + 1, &iptr));
320:     PetscCall(PetscMalloc1(nnz, &jptr));
321:     PetscCall(PetscMalloc1(nnz, &data));
322:   } else {
323:     PetscInt  nr;
324:     PetscBool done;
325:     PetscCall(MatISGetLocalMat(*B, &lA));
326:     PetscCall(MatGetRowIJ(lA, 0, PETSC_FALSE, PETSC_FALSE, &nr, (const PetscInt **)&iptr, (const PetscInt **)&jptr, &done));
327:     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);
328:     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);
329:     PetscCall(MatSeqAIJGetArrayWrite(lA, &data));
330:   }
331:   /* merge local matrices */
332:   ii  = iptr;
333:   jj  = jptr;
334:   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++ */
335:   *ii = *(hdi++) + *(hoi++);
336:   for (jd = 0, jo = 0, cum = 0; *ii < nnz; cum++) {
337:     PetscScalar *aold = (PetscScalar *)aa;
338:     PetscInt    *jold = jj, nc = jd + jo;
339:     for (; jd < *hdi; jd++) {
340:       *jj++ = *hdj++;
341:       *aa++ = *hdd++;
342:     }
343:     for (; jo < *hoi; jo++) {
344:       *jj++ = *hoj++ + dc;
345:       *aa++ = *hod++;
346:     }
347:     *(++ii) = *(hdi++) + *(hoi++);
348:     PetscCall(PetscSortIntWithScalarArray(jd + jo - nc, jold, aold));
349:   }
350:   for (; cum < dr; cum++) *(++ii) = nnz;
351:   if (reuse != MAT_REUSE_MATRIX) {
352:     Mat_SeqAIJ *a;

354:     PetscCall(MatCreateSeqAIJWithArrays(PETSC_COMM_SELF, dr, dc + oc, iptr, jptr, data, &lA));
355:     /* hack SeqAIJ */
356:     a          = (Mat_SeqAIJ *)lA->data;
357:     a->free_a  = PETSC_TRUE;
358:     a->free_ij = PETSC_TRUE;
359:     if (lmattype) PetscCall(MatConvert(lA, lmattype, MAT_INPLACE_MATRIX, &lA));
360:     PetscCall(MatISSetLocalMat(*B, lA));
361:     PetscCall(MatDestroy(&lA));
362:   } else {
363:     PetscCall(MatSeqAIJRestoreArrayWrite(lA, &data));
364:   }
365:   PetscCall(MatAssemblyBegin(*B, MAT_FINAL_ASSEMBLY));
366:   PetscCall(MatAssemblyEnd(*B, MAT_FINAL_ASSEMBLY));
367:   if (reuse == MAT_INPLACE_MATRIX) PetscCall(MatHeaderReplace(A, B));
368:   if (freeparcsr) PetscCallExternal(hypre_ParCSRMatrixDestroy, hA);
369:   PetscFunctionReturn(PETSC_SUCCESS);
370: }

372: static PetscErrorCode MatHYPRE_DestroyCOOMat(Mat mat)
373: {
374:   Mat_HYPRE *hA = (Mat_HYPRE *)mat->data;

376:   PetscFunctionBegin;
377:   if (hA->cooMat) { /* If cooMat is present we need to destroy the column indices */
378:     PetscCall(MatDestroy(&hA->cooMat));
379:     if (hA->cooMatAttached) {
380:       hypre_CSRMatrix     *csr;
381:       hypre_ParCSRMatrix  *parcsr;
382:       HYPRE_MemoryLocation mem;

384:       PetscCallExternal(HYPRE_IJMatrixGetObject, hA->ij, (void **)&parcsr);
385:       csr = hypre_ParCSRMatrixDiag(parcsr);
386:       if (csr) {
387:         mem = hypre_CSRMatrixMemoryLocation(csr);
388:         PetscStackCallExternalVoid("hypre_TFree", hypre_TFree(hypre_CSRMatrixJ(csr), mem));
389:         PetscStackCallExternalVoid("hypre_TFree", hypre_TFree(hypre_CSRMatrixBigJ(csr), mem));
390:       }
391:       csr = hypre_ParCSRMatrixOffd(parcsr);
392:       if (csr) {
393:         mem = hypre_CSRMatrixMemoryLocation(csr);
394:         PetscStackCallExternalVoid("hypre_TFree", hypre_TFree(hypre_CSRMatrixJ(csr), mem));
395:         PetscStackCallExternalVoid("hypre_TFree", hypre_TFree(hypre_CSRMatrixBigJ(csr), mem));
396:       }
397:     }
398:   }
399:   hA->cooMatAttached = PETSC_FALSE;
400:   PetscFunctionReturn(PETSC_SUCCESS);
401: }

403: static PetscErrorCode MatHYPRE_CreateCOOMat(Mat mat)
404: {
405:   MPI_Comm    comm;
406:   PetscMPIInt size;
407:   PetscLayout rmap, cmap;
408:   Mat_HYPRE  *hmat    = (Mat_HYPRE *)mat->data;
409:   MatType     matType = MATAIJ; /* default type of cooMat */

411:   PetscFunctionBegin;
412:   /* Build an agent matrix cooMat with AIJ format
413:      It has the same sparsity pattern as mat, and also shares the data array with mat. We use cooMat to do the COO work.
414:    */
415:   PetscCall(PetscObjectGetComm((PetscObject)mat, &comm));
416:   PetscCallMPI(MPI_Comm_size(comm, &size));
417:   PetscCall(PetscLayoutSetUp(mat->rmap));
418:   PetscCall(PetscLayoutSetUp(mat->cmap));
419:   PetscCall(MatGetLayouts(mat, &rmap, &cmap));

421: #if defined(PETSC_HAVE_HYPRE_DEVICE)
422:   if (!mat->boundtocpu) { /* mat will be on device, so will cooMat */
423:   #if defined(HYPRE_USING_HIP)
424:     matType = MATAIJHIPSPARSE;
425:   #elif defined(HYPRE_USING_CUDA)
426:     matType = MATAIJCUSPARSE;
427:   #elif defined(HYPRE_USING_SYCL) && defined(PETSC_HAVE_KOKKOS_KERNELS)
428:     matType = MATAIJKOKKOS;
429:   #else
430:     SETERRQ(comm, PETSC_ERR_SUP, "No HYPRE device available. Suggest re-installing with Kokkos Kernels");
431:   #endif
432:   }
433: #endif

435:   /* Do COO preallocation through cooMat */
436:   PetscCall(MatHYPRE_DestroyCOOMat(mat));
437:   PetscCall(MatCreate(comm, &hmat->cooMat));
438:   PetscCall(MatSetType(hmat->cooMat, matType));
439:   PetscCall(MatSetLayouts(hmat->cooMat, rmap, cmap));

441:   /* allocate local matrices if needed */
442:   PetscCall(MatMPIAIJSetPreallocation(hmat->cooMat, 0, NULL, 0, NULL));
443:   PetscFunctionReturn(PETSC_SUCCESS);
444: }

446: /* Attach cooMat data array to hypre matrix.
447:    When AIJCUPMSPARSE will support raw device pointers and not THRUSTARRAY
448:    we should swap the arrays: i.e., attach hypre matrix array to cooMat
449:    This is because hypre should be in charge of handling the memory,
450:    cooMat is only a way to reuse PETSc COO code.
451:    attaching the memory will then be done at MatSetValuesCOO time and it will dynamically
452:    support hypre matrix migrating to host.
453: */
454: static PetscErrorCode MatHYPRE_AttachCOOMat(Mat mat)
455: {
456:   Mat_HYPRE           *hmat = (Mat_HYPRE *)mat->data;
457:   hypre_CSRMatrix     *diag, *offd;
458:   hypre_ParCSRMatrix  *parCSR;
459:   HYPRE_MemoryLocation hmem = HYPRE_MEMORY_HOST;
460:   PetscMemType         pmem;
461:   Mat                  A, B;
462:   PetscScalar         *a;
463:   PetscMPIInt          size;
464:   MPI_Comm             comm;

466:   PetscFunctionBegin;
467:   PetscCheck(hmat->cooMat, PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "HYPRE COO delegate matrix has not been created yet");
468:   if (hmat->cooMatAttached) PetscFunctionReturn(PETSC_SUCCESS);
469:   PetscCheck(hmat->cooMat->preallocated, PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "HYPRE COO delegate matrix is not preallocated");
470:   PetscCall(PetscObjectSetName((PetscObject)hmat->cooMat, "_internal_COO_mat_for_hypre"));
471:   PetscCall(PetscObjectGetComm((PetscObject)mat, &comm));
472:   PetscCallMPI(MPI_Comm_size(comm, &size));

474:   /* Alias cooMat's data array to IJMatrix's */
475:   PetscCallExternal(HYPRE_IJMatrixGetObject, hmat->ij, (void **)&parCSR);
476:   diag = hypre_ParCSRMatrixDiag(parCSR);
477:   offd = hypre_ParCSRMatrixOffd(parCSR);

479:   A = (size == 1) ? hmat->cooMat : ((Mat_MPIAIJ *)hmat->cooMat->data)->A;
480:   B = (size == 1) ? NULL : ((Mat_MPIAIJ *)hmat->cooMat->data)->B;

482:   PetscCall(PetscObjectSetName((PetscObject)A, "_internal_COO_mat_for_hypre"));
483:   hmem = hypre_CSRMatrixMemoryLocation(diag);
484:   PetscCall(MatSeqAIJGetCSRAndMemType(A, NULL, NULL, &a, &pmem));
485:   PetscAssert((PetscMemTypeHost(pmem) && hmem == HYPRE_MEMORY_HOST) || (PetscMemTypeDevice(pmem) && hmem == HYPRE_MEMORY_DEVICE), comm, PETSC_ERR_PLIB, "PETSc and hypre's memory types mismatch");
486:   PetscStackCallExternalVoid("hypre_TFree", hypre_TFree(hypre_CSRMatrixData(diag), hmem));
487:   hypre_CSRMatrixData(diag)     = (HYPRE_Complex *)a;
488:   hypre_CSRMatrixOwnsData(diag) = 0; /* Take ownership of (j,a) away from hypre. As a result, we need to free them on our own */

490:   if (B) {
491:     hmem = hypre_CSRMatrixMemoryLocation(offd);
492:     PetscCall(MatSeqAIJGetCSRAndMemType(B, NULL, NULL, &a, &pmem));
493:     PetscAssert((PetscMemTypeHost(pmem) && hmem == HYPRE_MEMORY_HOST) || (PetscMemTypeDevice(pmem) && hmem == HYPRE_MEMORY_DEVICE), comm, PETSC_ERR_PLIB, "PETSc and hypre's memory types mismatch");
494:     PetscStackCallExternalVoid("hypre_TFree", hypre_TFree(hypre_CSRMatrixData(offd), hmem));
495:     hypre_CSRMatrixData(offd)     = (HYPRE_Complex *)a;
496:     hypre_CSRMatrixOwnsData(offd) = 0;
497:   }
498:   hmat->cooMatAttached = PETSC_TRUE;
499:   PetscFunctionReturn(PETSC_SUCCESS);
500: }

502: // Build COO's coordinate list i[], j[] based on CSR's i[], j[] arrays and the number of local rows 'n'
503: static PetscErrorCode CSRtoCOO_Private(PetscInt n, const PetscInt ii[], const PetscInt jj[], PetscCount *ncoo, PetscInt **coo_i, PetscInt **coo_j)
504: {
505:   PetscInt *cooi, *cooj;

507:   PetscFunctionBegin;
508:   *ncoo = ii[n];
509:   PetscCall(PetscMalloc2(*ncoo, &cooi, *ncoo, &cooj));
510:   for (PetscInt i = 0; i < n; i++) {
511:     for (PetscInt j = ii[i]; j < ii[i + 1]; j++) cooi[j] = i;
512:   }
513:   PetscCall(PetscArraycpy(cooj, jj, *ncoo));
514:   *coo_i = cooi;
515:   *coo_j = cooj;
516:   PetscFunctionReturn(PETSC_SUCCESS);
517: }

519: // Similar to CSRtoCOO_Private, but the CSR's i[], j[] are of type HYPRE_Int
520: static PetscErrorCode CSRtoCOO_HYPRE_Int_Private(PetscInt n, const HYPRE_Int ii[], const HYPRE_Int jj[], PetscCount *ncoo, PetscInt **coo_i, PetscInt **coo_j)
521: {
522:   PetscInt *cooi, *cooj;

524:   PetscFunctionBegin;
525:   *ncoo = ii[n];
526:   PetscCall(PetscMalloc2(*ncoo, &cooi, *ncoo, &cooj));
527:   for (PetscInt i = 0; i < n; i++) {
528:     for (HYPRE_Int j = ii[i]; j < ii[i + 1]; j++) cooi[j] = i;
529:   }
530:   for (PetscCount i = 0; i < *ncoo; i++) cooj[i] = jj[i];
531:   *coo_i = cooi;
532:   *coo_j = cooj;
533:   PetscFunctionReturn(PETSC_SUCCESS);
534: }

536: // Build a COO data structure for the seqaij matrix, as if the nonzeros are laid out in the same order as in the CSR
537: static PetscErrorCode MatSeqAIJGetCOO_Private(Mat A, PetscCount *ncoo, PetscInt **coo_i, PetscInt **coo_j)
538: {
539:   PetscInt        n;
540:   const PetscInt *ii, *jj;
541:   PetscBool       done;

543:   PetscFunctionBegin;
544:   PetscCall(MatGetRowIJ(A, 0, PETSC_FALSE, PETSC_FALSE, &n, &ii, &jj, &done));
545:   PetscCheck(done, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Failure for MatGetRowIJ");
546:   PetscCall(CSRtoCOO_Private(n, ii, jj, ncoo, coo_i, coo_j));
547:   PetscCall(MatRestoreRowIJ(A, 0, PETSC_FALSE, PETSC_FALSE, &n, &ii, &jj, &done));
548:   PetscCheck(done, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Failure for MatRestoreRowIJ");
549:   PetscFunctionReturn(PETSC_SUCCESS);
550: }

552: // Build a COO data structure for the hypreCSRMatrix, as if the nonzeros are laid out in the same order as in the hypreCSRMatrix
553: static PetscErrorCode hypreCSRMatrixGetCOO_Private(hypre_CSRMatrix *A, PetscCount *ncoo, PetscInt **coo_i, PetscInt **coo_j)
554: {
555:   PetscInt             n = hypre_CSRMatrixNumRows(A);
556:   HYPRE_Int           *ii, *jj;
557:   HYPRE_MemoryLocation mem = HYPRE_MEMORY_HOST;

559:   PetscFunctionBegin;
560: #if defined(PETSC_HAVE_HYPRE_DEVICE)
561:   mem = hypre_CSRMatrixMemoryLocation(A);
562:   if (mem != HYPRE_MEMORY_HOST) {
563:     PetscCount nnz = hypre_CSRMatrixNumNonzeros(A);
564:     PetscCall(PetscMalloc2(n + 1, &ii, nnz, &jj));
565:     hypre_TMemcpy(ii, hypre_CSRMatrixI(A), HYPRE_Int, n + 1, HYPRE_MEMORY_HOST, mem);
566:     hypre_TMemcpy(jj, hypre_CSRMatrixJ(A), HYPRE_Int, nnz, HYPRE_MEMORY_HOST, mem);
567:   } else {
568: #else
569:   {
570: #endif
571:     ii = hypre_CSRMatrixI(A);
572:     jj = hypre_CSRMatrixJ(A);
573:   }
574:   PetscCall(CSRtoCOO_HYPRE_Int_Private(n, ii, jj, ncoo, coo_i, coo_j));
575:   if (mem != HYPRE_MEMORY_HOST) PetscCall(PetscFree2(ii, jj));
576:   PetscFunctionReturn(PETSC_SUCCESS);
577: }

579: static PetscErrorCode MatSetValuesCOOFromCSRMatrix_Private(Mat A, hypre_CSRMatrix *H)
580: {
581:   PetscBool            iscpu = PETSC_TRUE;
582:   PetscScalar         *a;
583:   HYPRE_MemoryLocation mem = HYPRE_MEMORY_HOST;

585:   PetscFunctionBegin;
586: #if defined(PETSC_HAVE_HYPRE_DEVICE)
587:   mem = hypre_CSRMatrixMemoryLocation(H);
588:   PetscCall(PetscObjectTypeCompare((PetscObject)A, MATSEQAIJ, &iscpu));
589: #endif
590:   if (iscpu && mem != HYPRE_MEMORY_HOST) {
591:     PetscCount nnz = hypre_CSRMatrixNumNonzeros(H);
592:     PetscCall(PetscMalloc1(nnz, &a));
593:     hypre_TMemcpy(a, hypre_CSRMatrixData(H), PetscScalar, nnz, HYPRE_MEMORY_HOST, mem);
594:   } else {
595:     a = (PetscScalar *)hypre_CSRMatrixData(H);
596:   }
597:   PetscCall(MatSetValuesCOO(A, a, INSERT_VALUES));
598:   if (iscpu && mem != HYPRE_MEMORY_HOST) PetscCall(PetscFree(a));
599:   PetscFunctionReturn(PETSC_SUCCESS);
600: }

602: PETSC_INTERN PetscErrorCode MatConvert_AIJ_HYPRE(Mat A, MatType type, MatReuse reuse, Mat *B)
603: {
604:   MPI_Comm     comm = PetscObjectComm((PetscObject)A);
605:   Mat          M = NULL, dH = NULL, oH = NULL, dA = NULL, oA = NULL;
606:   PetscBool    ismpiaij, issbaij, isbaij, boundtocpu = PETSC_TRUE;
607:   Mat_HYPRE   *hA;
608:   PetscMemType memtype = PETSC_MEMTYPE_HOST;

610:   PetscFunctionBegin;
611:   if (PetscDefined(HAVE_HYPRE_DEVICE)) {
612:     PetscCall(MatGetCurrentMemType(A, &memtype));
613:     PetscHYPREInitialize();
614:     boundtocpu = PetscMemTypeHost(memtype) ? PETSC_TRUE : PETSC_FALSE;
615:     PetscCallExternal(HYPRE_SetMemoryLocation, boundtocpu ? HYPRE_MEMORY_HOST : HYPRE_MEMORY_DEVICE);
616:   }

618:   PetscCall(PetscObjectTypeCompareAny((PetscObject)A, &issbaij, MATSEQSBAIJ, MATMPIBAIJ, ""));
619:   PetscCall(PetscObjectTypeCompareAny((PetscObject)A, &isbaij, MATSEQBAIJ, MATMPIBAIJ, ""));
620:   if (isbaij || issbaij) { /* handle BAIJ and SBAIJ */
621:     PetscBool ismpi;
622:     MatType   newtype;

624:     PetscCall(PetscObjectTypeCompareAny((PetscObject)A, &ismpi, MATMPISBAIJ, MATMPIBAIJ, ""));
625:     newtype = ismpi ? MATMPIAIJ : MATSEQAIJ;
626:     if (reuse == MAT_REUSE_MATRIX) {
627:       PetscCall(MatConvert(*B, newtype, MAT_INPLACE_MATRIX, B));
628:       PetscCall(MatConvert(A, newtype, MAT_REUSE_MATRIX, B));
629:       PetscCall(MatConvert(*B, MATHYPRE, MAT_INPLACE_MATRIX, B));
630:     } else if (reuse == MAT_INITIAL_MATRIX) {
631:       PetscCall(MatConvert(A, newtype, MAT_INITIAL_MATRIX, B));
632:       PetscCall(MatConvert(*B, MATHYPRE, MAT_INPLACE_MATRIX, B));
633:     } else {
634:       PetscCall(MatConvert(A, newtype, MAT_INPLACE_MATRIX, &A));
635:       PetscCall(MatConvert(A, MATHYPRE, MAT_INPLACE_MATRIX, &A));
636:     }
637: #if defined(PETSC_HAVE_DEVICE)
638:     (*B)->boundtocpu = boundtocpu;
639: #endif
640:     PetscFunctionReturn(PETSC_SUCCESS);
641:   }

643:   dA = A;
644:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIAIJ, &ismpiaij));
645:   if (ismpiaij) PetscCall(MatMPIAIJGetSeqAIJ(A, &dA, &oA, NULL));

647:   if (reuse != MAT_REUSE_MATRIX) {
648:     PetscCount coo_n;
649:     PetscInt  *coo_i, *coo_j;

651:     PetscCall(MatCreate(comm, &M));
652:     PetscCall(MatSetType(M, MATHYPRE));
653:     PetscCall(MatSetSizes(M, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N));
654:     PetscCall(MatSetOption(M, MAT_SORTED_FULL, PETSC_TRUE));
655:     PetscCall(MatSetOption(M, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE));

657:     hA = (Mat_HYPRE *)M->data;
658:     PetscCall(MatHYPRE_CreateFromMat(A, hA));
659:     PetscCall(MatHYPRE_IJMatrixCopyIJ(A, hA->ij));

661:     PetscCall(MatHYPRE_CreateCOOMat(M));

663:     dH = hA->cooMat;
664:     PetscCall(PetscObjectBaseTypeCompare((PetscObject)hA->cooMat, MATMPIAIJ, &ismpiaij));
665:     if (ismpiaij) PetscCall(MatMPIAIJGetSeqAIJ(hA->cooMat, &dH, &oH, NULL));

667:     PetscCall(PetscObjectSetName((PetscObject)dH, "_internal_COO_mat_for_hypre"));
668:     PetscCall(MatSeqAIJGetCOO_Private(dA, &coo_n, &coo_i, &coo_j));
669:     PetscCall(MatSetPreallocationCOO(dH, coo_n, coo_i, coo_j));
670:     PetscCall(PetscFree2(coo_i, coo_j));
671:     if (oH) {
672:       PetscCall(PetscLayoutDestroy(&oH->cmap));
673:       PetscCall(PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)oH), oA->cmap->n, oA->cmap->n, 1, &oH->cmap));
674:       PetscCall(MatSeqAIJGetCOO_Private(oA, &coo_n, &coo_i, &coo_j));
675:       PetscCall(MatSetPreallocationCOO(oH, coo_n, coo_i, coo_j));
676:       PetscCall(PetscFree2(coo_i, coo_j));
677:     }
678:     hA->cooMat->assembled = PETSC_TRUE;

680:     M->preallocated = PETSC_TRUE;
681:     PetscCall(MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY));
682:     PetscCall(MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY));

684:     PetscCall(MatHYPRE_AttachCOOMat(M));
685:     if (reuse == MAT_INITIAL_MATRIX) *B = M;
686:   } else M = *B;

688:   hA = (Mat_HYPRE *)M->data;
689:   PetscCheck(hA->cooMat, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "HYPRE COO delegate matrix has not been created yet");

691:   dH = hA->cooMat;
692:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)hA->cooMat, MATMPIAIJ, &ismpiaij));
693:   if (ismpiaij) PetscCall(MatMPIAIJGetSeqAIJ(hA->cooMat, &dH, &oH, NULL));

695:   PetscScalar *a;
696:   PetscCall(MatSeqAIJGetCSRAndMemType(dA, NULL, NULL, &a, NULL));
697:   PetscCall(MatSetValuesCOO(dH, a, INSERT_VALUES));
698:   if (oH) {
699:     PetscCall(MatSeqAIJGetCSRAndMemType(oA, NULL, NULL, &a, NULL));
700:     PetscCall(MatSetValuesCOO(oH, a, INSERT_VALUES));
701:   }

703:   if (reuse == MAT_INPLACE_MATRIX) PetscCall(MatHeaderReplace(A, &M));
704: #if defined(PETSC_HAVE_DEVICE)
705:   (*B)->boundtocpu = boundtocpu;
706: #endif
707:   PetscFunctionReturn(PETSC_SUCCESS);
708: }

710: static PetscErrorCode MatConvert_HYPRE_AIJ(Mat A, MatType mtype, MatReuse reuse, Mat *B)
711: {
712:   Mat                 M, dA = NULL, oA = NULL;
713:   hypre_ParCSRMatrix *parcsr;
714:   hypre_CSRMatrix    *dH, *oH;
715:   MPI_Comm            comm;
716:   PetscBool           ismpiaij, isseqaij;

718:   PetscFunctionBegin;
719:   comm = PetscObjectComm((PetscObject)A);
720:   if (reuse == MAT_REUSE_MATRIX) {
721:     PetscCall(PetscObjectBaseTypeCompare((PetscObject)*B, MATMPIAIJ, &ismpiaij));
722:     PetscCall(PetscObjectBaseTypeCompare((PetscObject)*B, MATSEQAIJ, &isseqaij));
723:     PetscCheck(ismpiaij || isseqaij, comm, PETSC_ERR_SUP, "Only MATMPIAIJ or MATSEQAIJ base types are supported");
724:   }
725:   PetscCall(MatHYPREGetParCSR(A, &parcsr));
726: #if defined(PETSC_HAVE_HYPRE_DEVICE)
727:   if (HYPRE_MEMORY_DEVICE == hypre_ParCSRMatrixMemoryLocation(parcsr)) {
728:     PetscBool isaij;

730:     PetscCall(PetscStrcmp(mtype, MATAIJ, &isaij));
731:     if (isaij) {
732:       PetscMPIInt size;

734:       PetscCallMPI(MPI_Comm_size(comm, &size));
735:   #if defined(HYPRE_USING_HIP)
736:       mtype = size > 1 ? MATMPIAIJHIPSPARSE : MATSEQAIJHIPSPARSE;
737:   #elif defined(HYPRE_USING_CUDA)
738:       mtype = size > 1 ? MATMPIAIJCUSPARSE : MATSEQAIJCUSPARSE;
739:   #else
740:       mtype = size > 1 ? MATMPIAIJ : MATSEQAIJ;
741:   #endif
742:     }
743:   }
744: #endif
745:   dH = hypre_ParCSRMatrixDiag(parcsr);
746:   oH = hypre_ParCSRMatrixOffd(parcsr);
747:   if (reuse != MAT_REUSE_MATRIX) {
748:     PetscCount coo_n;
749:     PetscInt  *coo_i, *coo_j;

751:     PetscCall(MatCreate(comm, &M));
752:     PetscCall(MatSetType(M, mtype));
753:     PetscCall(MatSetSizes(M, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N));
754:     PetscCall(MatMPIAIJSetPreallocation(M, 0, NULL, 0, NULL));

756:     dA = M;
757:     PetscCall(PetscObjectBaseTypeCompare((PetscObject)M, MATMPIAIJ, &ismpiaij));
758:     if (ismpiaij) PetscCall(MatMPIAIJGetSeqAIJ(M, &dA, &oA, NULL));

760:     PetscCall(hypreCSRMatrixGetCOO_Private(dH, &coo_n, &coo_i, &coo_j));
761:     PetscCall(MatSetPreallocationCOO(dA, coo_n, coo_i, coo_j));
762:     PetscCall(PetscFree2(coo_i, coo_j));
763:     if (ismpiaij) {
764:       HYPRE_Int nc = hypre_CSRMatrixNumCols(oH);

766:       PetscCall(PetscLayoutDestroy(&oA->cmap));
767:       PetscCall(PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)oA), nc, nc, 1, &oA->cmap));
768:       PetscCall(hypreCSRMatrixGetCOO_Private(oH, &coo_n, &coo_i, &coo_j));
769:       PetscCall(MatSetPreallocationCOO(oA, coo_n, coo_i, coo_j));
770:       PetscCall(PetscFree2(coo_i, coo_j));

772:       /* garray */
773:       Mat_MPIAIJ   *aij    = (Mat_MPIAIJ *)M->data;
774:       HYPRE_BigInt *harray = hypre_ParCSRMatrixColMapOffd(parcsr);
775:       PetscInt     *garray;

777:       PetscCall(PetscFree(aij->garray));
778:       PetscCall(PetscMalloc1(nc, &garray));
779:       for (HYPRE_Int i = 0; i < nc; i++) garray[i] = (PetscInt)harray[i];
780:       aij->garray = garray;
781:       PetscCall(MatSetUpMultiply_MPIAIJ(M));
782:     }
783:     if (reuse == MAT_INITIAL_MATRIX) *B = M;
784:   } else M = *B;

786:   dA = M;
787:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)M, MATMPIAIJ, &ismpiaij));
788:   if (ismpiaij) PetscCall(MatMPIAIJGetSeqAIJ(M, &dA, &oA, NULL));
789:   PetscCall(MatSetValuesCOOFromCSRMatrix_Private(dA, dH));
790:   if (oA) PetscCall(MatSetValuesCOOFromCSRMatrix_Private(oA, oH));
791:   M->assembled = PETSC_TRUE;
792:   if (reuse == MAT_INPLACE_MATRIX) PetscCall(MatHeaderReplace(A, &M));
793:   PetscFunctionReturn(PETSC_SUCCESS);
794: }

796: static PetscErrorCode MatAIJGetParCSR_Private(Mat A, hypre_ParCSRMatrix **hA)
797: {
798:   hypre_ParCSRMatrix *tA;
799:   hypre_CSRMatrix    *hdiag, *hoffd;
800:   Mat_SeqAIJ         *diag, *offd;
801:   PetscInt           *garray, i, noffd, dnnz, onnz, *row_starts, *col_starts;
802:   MPI_Comm            comm = PetscObjectComm((PetscObject)A);
803:   PetscBool           ismpiaij, isseqaij;
804:   PetscBool           sameint = (PetscBool)(sizeof(PetscInt) == sizeof(HYPRE_Int));
805:   HYPRE_Int          *hdi = NULL, *hdj = NULL, *hoi = NULL, *hoj = NULL;
806:   PetscInt           *pdi = NULL, *pdj = NULL, *poi = NULL, *poj = NULL;
807:   PetscBool           iscuda, iship;
808: #if defined(PETSC_HAVE_DEVICE) && defined(PETSC_HAVE_HYPRE_DEVICE)
809:   PetscBool boundtocpu = A->boundtocpu;
810: #else
811:   PetscBool boundtocpu = PETSC_TRUE;
812: #endif

814:   PetscFunctionBegin;
815:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIAIJ, &ismpiaij));
816:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATSEQAIJ, &isseqaij));
817:   PetscCheck(ismpiaij || isseqaij, comm, PETSC_ERR_SUP, "Unsupported type %s", ((PetscObject)A)->type_name);
818:   PetscCall(PetscObjectTypeCompareAny((PetscObject)A, &iscuda, MATSEQAIJCUSPARSE, MATMPIAIJCUSPARSE, ""));
819:   PetscCall(PetscObjectTypeCompareAny((PetscObject)A, &iship, MATSEQAIJHIPSPARSE, MATMPIAIJHIPSPARSE, ""));
820:   PetscHYPREInitialize();
821:   if (ismpiaij) {
822:     Mat_MPIAIJ *a = (Mat_MPIAIJ *)A->data;

824:     diag = (Mat_SeqAIJ *)a->A->data;
825:     offd = (Mat_SeqAIJ *)a->B->data;
826:     if (!boundtocpu && (iscuda || iship)) {
827: #if defined(HYPRE_USING_CUDA) && defined(PETSC_HAVE_CUDA)
828:       if (iscuda) {
829:         sameint = PETSC_TRUE;
830:         PetscCall(MatSeqAIJCUSPARSEGetIJ(a->A, PETSC_FALSE, (const HYPRE_Int **)&hdi, (const HYPRE_Int **)&hdj));
831:         PetscCall(MatSeqAIJCUSPARSEGetIJ(a->B, PETSC_FALSE, (const HYPRE_Int **)&hoi, (const HYPRE_Int **)&hoj));
832:       }
833: #endif
834: #if defined(HYPRE_USING_HIP) && defined(PETSC_HAVE_HIP)
835:       if (iship) {
836:         sameint = PETSC_TRUE;
837:         PetscCall(MatSeqAIJHIPSPARSEGetIJ(a->A, PETSC_FALSE, (const HYPRE_Int **)&hdi, (const HYPRE_Int **)&hdj));
838:         PetscCall(MatSeqAIJHIPSPARSEGetIJ(a->B, PETSC_FALSE, (const HYPRE_Int **)&hoi, (const HYPRE_Int **)&hoj));
839:       }
840: #endif
841:     } else {
842:       boundtocpu = PETSC_TRUE;
843:       pdi        = diag->i;
844:       pdj        = diag->j;
845:       poi        = offd->i;
846:       poj        = offd->j;
847:       if (sameint) {
848:         hdi = (HYPRE_Int *)pdi;
849:         hdj = (HYPRE_Int *)pdj;
850:         hoi = (HYPRE_Int *)poi;
851:         hoj = (HYPRE_Int *)poj;
852:       }
853:     }
854:     garray = a->garray;
855:     noffd  = a->B->cmap->N;
856:     dnnz   = diag->nz;
857:     onnz   = offd->nz;
858:   } else {
859:     diag = (Mat_SeqAIJ *)A->data;
860:     offd = NULL;
861:     if (!boundtocpu && (iscuda || iship)) {
862: #if defined(HYPRE_USING_CUDA) && defined(PETSC_HAVE_CUDA)
863:       if (iscuda) {
864:         sameint = PETSC_TRUE;
865:         PetscCall(MatSeqAIJCUSPARSEGetIJ(A, PETSC_FALSE, (const HYPRE_Int **)&hdi, (const HYPRE_Int **)&hdj));
866:       }
867: #endif
868: #if defined(HYPRE_USING_HIP) && defined(PETSC_HAVE_HIP)
869:       if (iship) {
870:         sameint = PETSC_TRUE;
871:         PetscCall(MatSeqAIJHIPSPARSEGetIJ(A, PETSC_FALSE, (const HYPRE_Int **)&hdi, (const HYPRE_Int **)&hdj));
872:       }
873: #endif
874:     } else {
875:       boundtocpu = PETSC_TRUE;
876:       pdi        = diag->i;
877:       pdj        = diag->j;
878:       if (sameint) {
879:         hdi = (HYPRE_Int *)pdi;
880:         hdj = (HYPRE_Int *)pdj;
881:       }
882:     }
883:     garray = NULL;
884:     noffd  = 0;
885:     dnnz   = diag->nz;
886:     onnz   = 0;
887:   }

889:   /* create a temporary ParCSR */
890:   if (HYPRE_AssumedPartitionCheck()) {
891:     PetscMPIInt myid;

893:     PetscCallMPI(MPI_Comm_rank(comm, &myid));
894:     row_starts = A->rmap->range + myid;
895:     col_starts = A->cmap->range + myid;
896:   } else {
897:     row_starts = A->rmap->range;
898:     col_starts = A->cmap->range;
899:   }
900:   tA = hypre_ParCSRMatrixCreate(comm, A->rmap->N, A->cmap->N, (HYPRE_BigInt *)row_starts, (HYPRE_BigInt *)col_starts, noffd, dnnz, onnz);
901: #if defined(hypre_ParCSRMatrixOwnsRowStarts)
902:   hypre_ParCSRMatrixSetRowStartsOwner(tA, 0);
903:   hypre_ParCSRMatrixSetColStartsOwner(tA, 0);
904: #endif

906:   /* set diagonal part */
907:   hdiag = hypre_ParCSRMatrixDiag(tA);
908:   if (!sameint) { /* malloc CSR pointers */
909:     PetscCall(PetscMalloc2(A->rmap->n + 1, &hdi, dnnz, &hdj));
910:     for (i = 0; i < A->rmap->n + 1; i++) hdi[i] = (HYPRE_Int)pdi[i];
911:     for (i = 0; i < dnnz; i++) hdj[i] = (HYPRE_Int)pdj[i];
912:   }
913:   hypre_CSRMatrixI(hdiag)           = hdi;
914:   hypre_CSRMatrixJ(hdiag)           = hdj;
915:   hypre_CSRMatrixData(hdiag)        = (HYPRE_Complex *)diag->a;
916:   hypre_CSRMatrixNumNonzeros(hdiag) = diag->nz;
917:   hypre_CSRMatrixSetDataOwner(hdiag, 0);

919:   /* set off-diagonal part */
920:   hoffd = hypre_ParCSRMatrixOffd(tA);
921:   if (offd) {
922:     if (!sameint) { /* malloc CSR pointers */
923:       PetscCall(PetscMalloc2(A->rmap->n + 1, &hoi, onnz, &hoj));
924:       for (i = 0; i < A->rmap->n + 1; i++) hoi[i] = (HYPRE_Int)poi[i];
925:       for (i = 0; i < onnz; i++) hoj[i] = (HYPRE_Int)poj[i];
926:     }
927:     hypre_CSRMatrixI(hoffd)           = hoi;
928:     hypre_CSRMatrixJ(hoffd)           = hoj;
929:     hypre_CSRMatrixData(hoffd)        = (HYPRE_Complex *)offd->a;
930:     hypre_CSRMatrixNumNonzeros(hoffd) = offd->nz;
931:     hypre_CSRMatrixSetDataOwner(hoffd, 0);
932:   }
933: #if defined(PETSC_HAVE_HYPRE_DEVICE)
934:   PetscCallExternal(hypre_ParCSRMatrixInitialize_v2, tA, !boundtocpu ? HYPRE_MEMORY_DEVICE : HYPRE_MEMORY_HOST);
935: #else
936:   #if PETSC_PKG_HYPRE_VERSION_LT(2, 18, 0)
937:   PetscCallExternal(hypre_ParCSRMatrixInitialize, tA);
938:   #else
939:   PetscCallExternal(hypre_ParCSRMatrixInitialize_v2, tA, HYPRE_MEMORY_HOST);
940:   #endif
941: #endif

943:   /* MatrixSetRownnz comes after MatrixInitialize, so the first uses the right memory location */
944:   hypre_CSRMatrixSetRownnz(hdiag);
945:   if (offd) hypre_CSRMatrixSetRownnz(hoffd);

947:   hypre_TFree(hypre_ParCSRMatrixColMapOffd(tA), HYPRE_MEMORY_HOST);
948:   hypre_ParCSRMatrixSetNumNonzeros(tA);
949:   hypre_ParCSRMatrixColMapOffd(tA) = (HYPRE_BigInt *)garray;
950:   if (!hypre_ParCSRMatrixCommPkg(tA)) PetscCallExternal(hypre_MatvecCommPkgCreate, tA);
951:   *hA = tA;
952:   PetscFunctionReturn(PETSC_SUCCESS);
953: }

955: static PetscErrorCode MatAIJRestoreParCSR_Private(Mat A, hypre_ParCSRMatrix **hA)
956: {
957:   hypre_CSRMatrix *hdiag, *hoffd;
958:   PetscBool        ismpiaij, sameint = (PetscBool)(sizeof(PetscInt) == sizeof(HYPRE_Int));
959:   PetscBool        iscuda, iship;

961:   PetscFunctionBegin;
962:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIAIJ, &ismpiaij));
963:   PetscCall(PetscObjectTypeCompareAny((PetscObject)A, &iscuda, MATSEQAIJCUSPARSE, MATMPIAIJCUSPARSE, ""));
964:   PetscCall(PetscObjectTypeCompareAny((PetscObject)A, &iship, MATSEQAIJHIPSPARSE, MATMPIAIJHIPSPARSE, ""));
965: #if defined(HYPRE_USING_CUDA) && defined(PETSC_HAVE_CUDA)
966:   if (iscuda) sameint = PETSC_TRUE;
967: #elif defined(HYPRE_USING_HIP) && defined(PETSC_HAVE_HIP)
968:   if (iship) sameint = PETSC_TRUE;
969: #endif
970:   hdiag = hypre_ParCSRMatrixDiag(*hA);
971:   hoffd = hypre_ParCSRMatrixOffd(*hA);
972:   /* free temporary memory allocated by PETSc
973:      set pointers to NULL before destroying tA */
974:   if (!sameint) {
975:     HYPRE_Int *hi, *hj;

977:     hi = hypre_CSRMatrixI(hdiag);
978:     hj = hypre_CSRMatrixJ(hdiag);
979:     PetscCall(PetscFree2(hi, hj));
980:     if (ismpiaij) {
981:       hi = hypre_CSRMatrixI(hoffd);
982:       hj = hypre_CSRMatrixJ(hoffd);
983:       PetscCall(PetscFree2(hi, hj));
984:     }
985:   }
986:   hypre_CSRMatrixI(hdiag)    = NULL;
987:   hypre_CSRMatrixJ(hdiag)    = NULL;
988:   hypre_CSRMatrixData(hdiag) = NULL;
989:   if (ismpiaij) {
990:     hypre_CSRMatrixI(hoffd)    = NULL;
991:     hypre_CSRMatrixJ(hoffd)    = NULL;
992:     hypre_CSRMatrixData(hoffd) = NULL;
993:   }
994:   hypre_ParCSRMatrixColMapOffd(*hA) = NULL;
995:   hypre_ParCSRMatrixDestroy(*hA);
996:   *hA = NULL;
997:   PetscFunctionReturn(PETSC_SUCCESS);
998: }

1000: /* calls RAP from BoomerAMG:
1001:    the resulting ParCSR will not own the column and row starts
1002:    It looks like we don't need to have the diagonal entries ordered first */
1003: static PetscErrorCode MatHYPRE_ParCSR_RAP(hypre_ParCSRMatrix *hR, hypre_ParCSRMatrix *hA, hypre_ParCSRMatrix *hP, hypre_ParCSRMatrix **hRAP)
1004: {
1005: #if defined(hypre_ParCSRMatrixOwnsRowStarts)
1006:   HYPRE_Int P_owns_col_starts, R_owns_row_starts;
1007: #endif

1009:   PetscFunctionBegin;
1010: #if defined(hypre_ParCSRMatrixOwnsRowStarts)
1011:   P_owns_col_starts = hypre_ParCSRMatrixOwnsColStarts(hP);
1012:   R_owns_row_starts = hypre_ParCSRMatrixOwnsRowStarts(hR);
1013: #endif
1014:   /* can be replaced by version test later */
1015: #if defined(PETSC_HAVE_HYPRE_DEVICE)
1016:   PetscStackPushExternal("hypre_ParCSRMatrixRAP");
1017:   *hRAP = hypre_ParCSRMatrixRAP(hR, hA, hP);
1018:   PetscStackPop;
1019: #else
1020:   PetscCallExternal(hypre_BoomerAMGBuildCoarseOperator, hR, hA, hP, hRAP);
1021:   PetscCallExternal(hypre_ParCSRMatrixSetNumNonzeros, *hRAP);
1022: #endif
1023:   /* hypre_BoomerAMGBuildCoarseOperator steals the col_starts from P and the row_starts from R */
1024: #if defined(hypre_ParCSRMatrixOwnsRowStarts)
1025:   hypre_ParCSRMatrixSetRowStartsOwner(*hRAP, 0);
1026:   hypre_ParCSRMatrixSetColStartsOwner(*hRAP, 0);
1027:   if (P_owns_col_starts) hypre_ParCSRMatrixSetColStartsOwner(hP, 1);
1028:   if (R_owns_row_starts) hypre_ParCSRMatrixSetRowStartsOwner(hR, 1);
1029: #endif
1030:   PetscFunctionReturn(PETSC_SUCCESS);
1031: }

1033: static PetscErrorCode MatPtAPNumeric_AIJ_AIJ_wHYPRE(Mat A, Mat P, Mat C)
1034: {
1035:   Mat                 B;
1036:   hypre_ParCSRMatrix *hA, *hP, *hPtAP = NULL;
1037:   Mat_Product        *product = C->product;

1039:   PetscFunctionBegin;
1040:   PetscCall(MatAIJGetParCSR_Private(A, &hA));
1041:   PetscCall(MatAIJGetParCSR_Private(P, &hP));
1042:   PetscCall(MatHYPRE_ParCSR_RAP(hP, hA, hP, &hPtAP));
1043:   PetscCall(MatCreateFromParCSR(hPtAP, MATAIJ, PETSC_OWN_POINTER, &B));

1045:   PetscCall(MatHeaderMerge(C, &B));
1046:   C->product = product;

1048:   PetscCall(MatAIJRestoreParCSR_Private(A, &hA));
1049:   PetscCall(MatAIJRestoreParCSR_Private(P, &hP));
1050:   PetscFunctionReturn(PETSC_SUCCESS);
1051: }

1053: PETSC_INTERN PetscErrorCode MatPtAPSymbolic_AIJ_AIJ_wHYPRE(Mat A, Mat P, PetscReal fill, Mat C)
1054: {
1055:   PetscFunctionBegin;
1056:   PetscCall(MatSetType(C, MATAIJ));
1057:   C->ops->ptapnumeric    = MatPtAPNumeric_AIJ_AIJ_wHYPRE;
1058:   C->ops->productnumeric = MatProductNumeric_PtAP;
1059:   PetscFunctionReturn(PETSC_SUCCESS);
1060: }

1062: static PetscErrorCode MatPtAPNumeric_AIJ_HYPRE(Mat A, Mat P, Mat C)
1063: {
1064:   Mat                 B;
1065:   Mat_HYPRE          *hP;
1066:   hypre_ParCSRMatrix *hA = NULL, *Pparcsr, *ptapparcsr = NULL;
1067:   HYPRE_Int           type;
1068:   MPI_Comm            comm = PetscObjectComm((PetscObject)A);
1069:   PetscBool           ishypre;

1071:   PetscFunctionBegin;
1072:   PetscCall(PetscObjectTypeCompare((PetscObject)P, MATHYPRE, &ishypre));
1073:   PetscCheck(ishypre, comm, PETSC_ERR_USER, "P should be of type %s", MATHYPRE);
1074:   hP = (Mat_HYPRE *)P->data;
1075:   PetscCallExternal(HYPRE_IJMatrixGetObjectType, hP->ij, &type);
1076:   PetscCheck(type == HYPRE_PARCSR, comm, PETSC_ERR_SUP, "Only HYPRE_PARCSR is supported");
1077:   PetscCallExternal(HYPRE_IJMatrixGetObject, hP->ij, (void **)&Pparcsr);

1079:   PetscCall(MatAIJGetParCSR_Private(A, &hA));
1080:   PetscCall(MatHYPRE_ParCSR_RAP(Pparcsr, hA, Pparcsr, &ptapparcsr));
1081:   PetscCall(MatAIJRestoreParCSR_Private(A, &hA));

1083:   /* create temporary matrix and merge to C */
1084:   PetscCall(MatCreateFromParCSR(ptapparcsr, ((PetscObject)C)->type_name, PETSC_OWN_POINTER, &B));
1085:   PetscCall(MatHeaderMerge(C, &B));
1086:   PetscFunctionReturn(PETSC_SUCCESS);
1087: }

1089: static PetscErrorCode MatPtAPNumeric_HYPRE_HYPRE(Mat A, Mat P, Mat C)
1090: {
1091:   Mat                 B;
1092:   hypre_ParCSRMatrix *Aparcsr, *Pparcsr, *ptapparcsr = NULL;
1093:   Mat_HYPRE          *hA, *hP;
1094:   PetscBool           ishypre;
1095:   HYPRE_Int           type;

1097:   PetscFunctionBegin;
1098:   PetscCall(PetscObjectTypeCompare((PetscObject)P, MATHYPRE, &ishypre));
1099:   PetscCheck(ishypre, PetscObjectComm((PetscObject)P), PETSC_ERR_USER, "P should be of type %s", MATHYPRE);
1100:   PetscCall(PetscObjectTypeCompare((PetscObject)A, MATHYPRE, &ishypre));
1101:   PetscCheck(ishypre, PetscObjectComm((PetscObject)A), PETSC_ERR_USER, "A should be of type %s", MATHYPRE);
1102:   hA = (Mat_HYPRE *)A->data;
1103:   hP = (Mat_HYPRE *)P->data;
1104:   PetscCallExternal(HYPRE_IJMatrixGetObjectType, hA->ij, &type);
1105:   PetscCheck(type == HYPRE_PARCSR, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Only HYPRE_PARCSR is supported");
1106:   PetscCallExternal(HYPRE_IJMatrixGetObjectType, hP->ij, &type);
1107:   PetscCheck(type == HYPRE_PARCSR, PetscObjectComm((PetscObject)P), PETSC_ERR_SUP, "Only HYPRE_PARCSR is supported");
1108:   PetscCallExternal(HYPRE_IJMatrixGetObject, hA->ij, (void **)&Aparcsr);
1109:   PetscCallExternal(HYPRE_IJMatrixGetObject, hP->ij, (void **)&Pparcsr);
1110:   PetscCall(MatHYPRE_ParCSR_RAP(Pparcsr, Aparcsr, Pparcsr, &ptapparcsr));
1111:   PetscCall(MatCreateFromParCSR(ptapparcsr, MATHYPRE, PETSC_OWN_POINTER, &B));
1112:   PetscCall(MatHeaderMerge(C, &B));
1113:   PetscFunctionReturn(PETSC_SUCCESS);
1114: }

1116: /* calls hypre_ParMatmul
1117:    hypre_ParMatMul uses hypre_ParMatrixCreate with the communicator of hA
1118:    hypre_ParMatrixCreate does not duplicate the communicator
1119:    It looks like we don't need to have the diagonal entries ordered first */
1120: static PetscErrorCode MatHYPRE_ParCSR_MatMatMult(hypre_ParCSRMatrix *hA, hypre_ParCSRMatrix *hB, hypre_ParCSRMatrix **hAB)
1121: {
1122:   PetscFunctionBegin;
1123:   /* can be replaced by version test later */
1124: #if defined(PETSC_HAVE_HYPRE_DEVICE)
1125:   PetscStackPushExternal("hypre_ParCSRMatMat");
1126:   *hAB = hypre_ParCSRMatMat(hA, hB);
1127: #else
1128:   PetscStackPushExternal("hypre_ParMatmul");
1129:   *hAB = hypre_ParMatmul(hA, hB);
1130: #endif
1131:   PetscStackPop;
1132:   PetscFunctionReturn(PETSC_SUCCESS);
1133: }

1135: static PetscErrorCode MatMatMultNumeric_AIJ_AIJ_wHYPRE(Mat A, Mat B, Mat C)
1136: {
1137:   Mat                 D;
1138:   hypre_ParCSRMatrix *hA, *hB, *hAB = NULL;
1139:   Mat_Product        *product = C->product;

1141:   PetscFunctionBegin;
1142:   PetscCall(MatAIJGetParCSR_Private(A, &hA));
1143:   PetscCall(MatAIJGetParCSR_Private(B, &hB));
1144:   PetscCall(MatHYPRE_ParCSR_MatMatMult(hA, hB, &hAB));
1145:   PetscCall(MatCreateFromParCSR(hAB, MATAIJ, PETSC_OWN_POINTER, &D));

1147:   PetscCall(MatHeaderMerge(C, &D));
1148:   C->product = product;

1150:   PetscCall(MatAIJRestoreParCSR_Private(A, &hA));
1151:   PetscCall(MatAIJRestoreParCSR_Private(B, &hB));
1152:   PetscFunctionReturn(PETSC_SUCCESS);
1153: }

1155: PETSC_INTERN PetscErrorCode MatMatMultSymbolic_AIJ_AIJ_wHYPRE(Mat A, Mat B, PetscReal fill, Mat C)
1156: {
1157:   PetscFunctionBegin;
1158:   PetscCall(MatSetType(C, MATAIJ));
1159:   C->ops->matmultnumeric = MatMatMultNumeric_AIJ_AIJ_wHYPRE;
1160:   C->ops->productnumeric = MatProductNumeric_AB;
1161:   PetscFunctionReturn(PETSC_SUCCESS);
1162: }

1164: static PetscErrorCode MatMatMultNumeric_HYPRE_HYPRE(Mat A, Mat B, Mat C)
1165: {
1166:   Mat                 D;
1167:   hypre_ParCSRMatrix *Aparcsr, *Bparcsr, *ABparcsr = NULL;
1168:   Mat_HYPRE          *hA, *hB;
1169:   PetscBool           ishypre;
1170:   HYPRE_Int           type;
1171:   Mat_Product        *product;

1173:   PetscFunctionBegin;
1174:   PetscCall(PetscObjectTypeCompare((PetscObject)B, MATHYPRE, &ishypre));
1175:   PetscCheck(ishypre, PetscObjectComm((PetscObject)B), PETSC_ERR_USER, "B should be of type %s", MATHYPRE);
1176:   PetscCall(PetscObjectTypeCompare((PetscObject)A, MATHYPRE, &ishypre));
1177:   PetscCheck(ishypre, PetscObjectComm((PetscObject)A), PETSC_ERR_USER, "A should be of type %s", MATHYPRE);
1178:   hA = (Mat_HYPRE *)A->data;
1179:   hB = (Mat_HYPRE *)B->data;
1180:   PetscCallExternal(HYPRE_IJMatrixGetObjectType, hA->ij, &type);
1181:   PetscCheck(type == HYPRE_PARCSR, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Only HYPRE_PARCSR is supported");
1182:   PetscCallExternal(HYPRE_IJMatrixGetObjectType, hB->ij, &type);
1183:   PetscCheck(type == HYPRE_PARCSR, PetscObjectComm((PetscObject)B), PETSC_ERR_SUP, "Only HYPRE_PARCSR is supported");
1184:   PetscCallExternal(HYPRE_IJMatrixGetObject, hA->ij, (void **)&Aparcsr);
1185:   PetscCallExternal(HYPRE_IJMatrixGetObject, hB->ij, (void **)&Bparcsr);
1186:   PetscCall(MatHYPRE_ParCSR_MatMatMult(Aparcsr, Bparcsr, &ABparcsr));
1187:   PetscCall(MatCreateFromParCSR(ABparcsr, MATHYPRE, PETSC_OWN_POINTER, &D));

1189:   /* need to use HeaderReplace because HeaderMerge messes up with the communicator */
1190:   product    = C->product; /* save it from MatHeaderReplace() */
1191:   C->product = NULL;
1192:   PetscCall(MatHeaderReplace(C, &D));
1193:   C->product             = product;
1194:   C->ops->matmultnumeric = MatMatMultNumeric_HYPRE_HYPRE;
1195:   C->ops->productnumeric = MatProductNumeric_AB;
1196:   PetscFunctionReturn(PETSC_SUCCESS);
1197: }

1199: PETSC_INTERN PetscErrorCode MatTransposeMatMatMultNumeric_AIJ_AIJ_AIJ_wHYPRE(Mat A, Mat B, Mat C, Mat D)
1200: {
1201:   Mat                 E;
1202:   hypre_ParCSRMatrix *hA, *hB, *hC, *hABC = NULL;

1204:   PetscFunctionBegin;
1205:   PetscCall(MatAIJGetParCSR_Private(A, &hA));
1206:   PetscCall(MatAIJGetParCSR_Private(B, &hB));
1207:   PetscCall(MatAIJGetParCSR_Private(C, &hC));
1208:   PetscCall(MatHYPRE_ParCSR_RAP(hA, hB, hC, &hABC));
1209:   PetscCall(MatCreateFromParCSR(hABC, MATAIJ, PETSC_OWN_POINTER, &E));
1210:   PetscCall(MatHeaderMerge(D, &E));
1211:   PetscCall(MatAIJRestoreParCSR_Private(A, &hA));
1212:   PetscCall(MatAIJRestoreParCSR_Private(B, &hB));
1213:   PetscCall(MatAIJRestoreParCSR_Private(C, &hC));
1214:   PetscFunctionReturn(PETSC_SUCCESS);
1215: }

1217: PETSC_INTERN PetscErrorCode MatTransposeMatMatMultSymbolic_AIJ_AIJ_AIJ_wHYPRE(Mat A, Mat B, Mat C, PetscReal fill, Mat D)
1218: {
1219:   PetscFunctionBegin;
1220:   PetscCall(MatSetType(D, MATAIJ));
1221:   PetscFunctionReturn(PETSC_SUCCESS);
1222: }

1224: static PetscErrorCode MatProductSymbolic_AB_HYPRE(Mat C)
1225: {
1226:   PetscFunctionBegin;
1227:   C->ops->productnumeric = MatProductNumeric_AB;
1228:   PetscFunctionReturn(PETSC_SUCCESS);
1229: }

1231: static PetscErrorCode MatProductSetFromOptions_HYPRE_AB(Mat C)
1232: {
1233:   Mat_Product *product = C->product;
1234:   PetscBool    Ahypre;

1236:   PetscFunctionBegin;
1237:   PetscCall(PetscObjectTypeCompare((PetscObject)product->A, MATHYPRE, &Ahypre));
1238:   if (Ahypre) { /* A is a Hypre matrix */
1239:     PetscCall(MatSetType(C, MATHYPRE));
1240:     C->ops->productsymbolic = MatProductSymbolic_AB_HYPRE;
1241:     C->ops->matmultnumeric  = MatMatMultNumeric_HYPRE_HYPRE;
1242:     PetscFunctionReturn(PETSC_SUCCESS);
1243:   }
1244:   PetscFunctionReturn(PETSC_SUCCESS);
1245: }

1247: static PetscErrorCode MatProductSymbolic_PtAP_HYPRE(Mat C)
1248: {
1249:   PetscFunctionBegin;
1250:   C->ops->productnumeric = MatProductNumeric_PtAP;
1251:   PetscFunctionReturn(PETSC_SUCCESS);
1252: }

1254: static PetscErrorCode MatProductSetFromOptions_HYPRE_PtAP(Mat C)
1255: {
1256:   Mat_Product *product = C->product;
1257:   PetscBool    flg;
1258:   PetscInt     type        = 0;
1259:   const char  *outTypes[4] = {"aij", "seqaij", "mpiaij", "hypre"};
1260:   PetscInt     ntype       = 4;
1261:   Mat          A           = product->A;
1262:   PetscBool    Ahypre;

1264:   PetscFunctionBegin;
1265:   PetscCall(PetscObjectTypeCompare((PetscObject)A, MATHYPRE, &Ahypre));
1266:   if (Ahypre) { /* A is a Hypre matrix */
1267:     PetscCall(MatSetType(C, MATHYPRE));
1268:     C->ops->productsymbolic = MatProductSymbolic_PtAP_HYPRE;
1269:     C->ops->ptapnumeric     = MatPtAPNumeric_HYPRE_HYPRE;
1270:     PetscFunctionReturn(PETSC_SUCCESS);
1271:   }

1273:   /* A is AIJ, P is Hypre, C = PtAP can be either AIJ or Hypre format */
1274:   /* Get runtime option */
1275:   if (product->api_user) {
1276:     PetscOptionsBegin(PetscObjectComm((PetscObject)C), ((PetscObject)C)->prefix, "MatPtAP_HYPRE", "Mat");
1277:     PetscCall(PetscOptionsEList("-matptap_hypre_outtype", "MatPtAP outtype", "MatPtAP outtype", outTypes, ntype, outTypes[type], &type, &flg));
1278:     PetscOptionsEnd();
1279:   } else {
1280:     PetscOptionsBegin(PetscObjectComm((PetscObject)C), ((PetscObject)C)->prefix, "MatProduct_PtAP_HYPRE", "Mat");
1281:     PetscCall(PetscOptionsEList("-mat_product_algorithm_hypre_outtype", "MatProduct_PtAP outtype", "MatProduct_PtAP", outTypes, ntype, outTypes[type], &type, &flg));
1282:     PetscOptionsEnd();
1283:   }

1285:   if (type == 0 || type == 1 || type == 2) {
1286:     PetscCall(MatSetType(C, MATAIJ));
1287:   } else if (type == 3) {
1288:     PetscCall(MatSetType(C, MATHYPRE));
1289:   } else SETERRQ(PetscObjectComm((PetscObject)C), PETSC_ERR_SUP, "MatPtAP outtype is not supported");
1290:   C->ops->productsymbolic = MatProductSymbolic_PtAP_HYPRE;
1291:   C->ops->ptapnumeric     = MatPtAPNumeric_AIJ_HYPRE;
1292:   PetscFunctionReturn(PETSC_SUCCESS);
1293: }

1295: static PetscErrorCode MatProductSetFromOptions_HYPRE(Mat C)
1296: {
1297:   Mat_Product *product = C->product;

1299:   PetscFunctionBegin;
1300:   switch (product->type) {
1301:   case MATPRODUCT_AB:
1302:     PetscCall(MatProductSetFromOptions_HYPRE_AB(C));
1303:     break;
1304:   case MATPRODUCT_PtAP:
1305:     PetscCall(MatProductSetFromOptions_HYPRE_PtAP(C));
1306:     break;
1307:   default:
1308:     break;
1309:   }
1310:   PetscFunctionReturn(PETSC_SUCCESS);
1311: }

1313: static PetscErrorCode MatMultTranspose_HYPRE(Mat A, Vec x, Vec y)
1314: {
1315:   PetscFunctionBegin;
1316:   PetscCall(MatHYPRE_MultKernel_Private(A, 1.0, x, 0.0, y, PETSC_TRUE));
1317:   PetscFunctionReturn(PETSC_SUCCESS);
1318: }

1320: static PetscErrorCode MatMult_HYPRE(Mat A, Vec x, Vec y)
1321: {
1322:   PetscFunctionBegin;
1323:   PetscCall(MatHYPRE_MultKernel_Private(A, 1.0, x, 0.0, y, PETSC_FALSE));
1324:   PetscFunctionReturn(PETSC_SUCCESS);
1325: }

1327: static PetscErrorCode MatMultAdd_HYPRE(Mat A, Vec x, Vec y, Vec z)
1328: {
1329:   PetscFunctionBegin;
1330:   if (y != z) PetscCall(VecCopy(y, z));
1331:   PetscCall(MatHYPRE_MultKernel_Private(A, 1.0, x, 1.0, z, PETSC_FALSE));
1332:   PetscFunctionReturn(PETSC_SUCCESS);
1333: }

1335: static PetscErrorCode MatMultTransposeAdd_HYPRE(Mat A, Vec x, Vec y, Vec z)
1336: {
1337:   PetscFunctionBegin;
1338:   if (y != z) PetscCall(VecCopy(y, z));
1339:   PetscCall(MatHYPRE_MultKernel_Private(A, 1.0, x, 1.0, z, PETSC_TRUE));
1340:   PetscFunctionReturn(PETSC_SUCCESS);
1341: }

1343: /* y = a * A * x + b * y or y = a * A^t * x + b * y depending on trans */
1344: static PetscErrorCode MatHYPRE_MultKernel_Private(Mat A, HYPRE_Complex a, Vec x, HYPRE_Complex b, Vec y, PetscBool trans)
1345: {
1346:   Mat_HYPRE          *hA = (Mat_HYPRE *)A->data;
1347:   hypre_ParCSRMatrix *parcsr;
1348:   hypre_ParVector    *hx, *hy;

1350:   PetscFunctionBegin;
1351:   if (trans) {
1352:     PetscCall(VecHYPRE_IJVectorPushVecRead(hA->b, x));
1353:     if (b != 0.0) PetscCall(VecHYPRE_IJVectorPushVec(hA->x, y));
1354:     else PetscCall(VecHYPRE_IJVectorPushVecWrite(hA->x, y));
1355:     PetscCallExternal(HYPRE_IJVectorGetObject, hA->b->ij, (void **)&hx);
1356:     PetscCallExternal(HYPRE_IJVectorGetObject, hA->x->ij, (void **)&hy);
1357:   } else {
1358:     PetscCall(VecHYPRE_IJVectorPushVecRead(hA->x, x));
1359:     if (b != 0.0) PetscCall(VecHYPRE_IJVectorPushVec(hA->b, y));
1360:     else PetscCall(VecHYPRE_IJVectorPushVecWrite(hA->b, y));
1361:     PetscCallExternal(HYPRE_IJVectorGetObject, hA->x->ij, (void **)&hx);
1362:     PetscCallExternal(HYPRE_IJVectorGetObject, hA->b->ij, (void **)&hy);
1363:   }
1364:   PetscCallExternal(HYPRE_IJMatrixGetObject, hA->ij, (void **)&parcsr);
1365:   if (trans) {
1366:     PetscCallExternal(hypre_ParCSRMatrixMatvecT, a, parcsr, hx, b, hy);
1367:   } else {
1368:     PetscCallExternal(hypre_ParCSRMatrixMatvec, a, parcsr, hx, b, hy);
1369:   }
1370:   PetscCall(VecHYPRE_IJVectorPopVec(hA->x));
1371:   PetscCall(VecHYPRE_IJVectorPopVec(hA->b));
1372:   PetscFunctionReturn(PETSC_SUCCESS);
1373: }

1375: static PetscErrorCode MatDestroy_HYPRE(Mat A)
1376: {
1377:   Mat_HYPRE *hA = (Mat_HYPRE *)A->data;

1379:   PetscFunctionBegin;
1380:   PetscCall(VecHYPRE_IJVectorDestroy(&hA->x));
1381:   PetscCall(VecHYPRE_IJVectorDestroy(&hA->b));
1382:   PetscCall(MatHYPRE_DestroyCOOMat(A)); /* must be called before destroying the individual CSR */
1383:   if (hA->ij) {
1384:     if (!hA->inner_free) hypre_IJMatrixObject(hA->ij) = NULL;
1385:     PetscCallExternal(HYPRE_IJMatrixDestroy, hA->ij);
1386:   }
1387:   if (hA->comm) PetscCall(PetscCommRestoreComm(PetscObjectComm((PetscObject)A), &hA->comm));

1389:   PetscCall(MatStashDestroy_Private(&A->stash));
1390:   PetscCall(PetscFree(hA->array));
1391:   if (hA->rows_d) PetscStackCallExternalVoid("hypre_Free", hypre_Free(hA->rows_d, HYPRE_MEMORY_DEVICE));

1393:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_hypre_aij_C", NULL));
1394:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_hypre_is_C", NULL));
1395:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_seqaij_hypre_C", NULL));
1396:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_mpiaij_hypre_C", NULL));
1397:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_seqaijhipsparse_hypre_C", NULL));
1398:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_mpiaijhipsparse_hypre_C", NULL));
1399:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_seqaijcusparse_hypre_C", NULL));
1400:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_mpiaijcusparse_hypre_C", NULL));
1401:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHYPRESetPreallocation_C", NULL));
1402:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHYPREGetParCSR_C", NULL));
1403:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSetPreallocationCOO_C", NULL));
1404:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSetValuesCOO_C", NULL));
1405:   PetscCall(PetscFree(A->data));
1406:   PetscFunctionReturn(PETSC_SUCCESS);
1407: }

1409: static PetscErrorCode MatSetUp_HYPRE(Mat A)
1410: {
1411:   PetscFunctionBegin;
1412:   if (!A->preallocated) PetscCall(MatHYPRESetPreallocation(A, PETSC_DEFAULT, NULL, PETSC_DEFAULT, NULL));
1413:   PetscFunctionReturn(PETSC_SUCCESS);
1414: }

1416: //TODO FIX hypre_CSRMatrixMatvecOutOfPlace
1417: #if defined(PETSC_HAVE_HYPRE_DEVICE)
1418: static PetscErrorCode MatBindToCPU_HYPRE(Mat A, PetscBool bind)
1419: {
1420:   Mat_HYPRE           *hA   = (Mat_HYPRE *)A->data;
1421:   HYPRE_MemoryLocation hmem = bind ? HYPRE_MEMORY_HOST : HYPRE_MEMORY_DEVICE;

1423:   PetscFunctionBegin;
1424:   A->boundtocpu = bind;
1425:   if (hA->ij && hypre_IJMatrixAssembleFlag(hA->ij) && hmem != hypre_IJMatrixMemoryLocation(hA->ij)) {
1426:     hypre_ParCSRMatrix *parcsr;
1427:     PetscCallExternal(HYPRE_IJMatrixGetObject, hA->ij, (void **)&parcsr);
1428:     PetscCallExternal(hypre_ParCSRMatrixMigrate, parcsr, hmem);
1429:   }
1430:   if (hA->x) PetscCall(VecHYPRE_IJBindToCPU(hA->x, bind));
1431:   if (hA->b) PetscCall(VecHYPRE_IJBindToCPU(hA->b, bind));
1432:   PetscFunctionReturn(PETSC_SUCCESS);
1433: }
1434: #endif

1436: static PetscErrorCode MatAssemblyEnd_HYPRE(Mat A, MatAssemblyType mode)
1437: {
1438:   Mat_HYPRE   *hA = (Mat_HYPRE *)A->data;
1439:   PetscMPIInt  n;
1440:   PetscInt     i, j, rstart, ncols, flg;
1441:   PetscInt    *row, *col;
1442:   PetscScalar *val;

1444:   PetscFunctionBegin;
1445:   PetscCheck(mode != MAT_FLUSH_ASSEMBLY, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "MAT_FLUSH_ASSEMBLY currently not supported with MATHYPRE");

1447:   if (!A->nooffprocentries) {
1448:     while (1) {
1449:       PetscCall(MatStashScatterGetMesg_Private(&A->stash, &n, &row, &col, &val, &flg));
1450:       if (!flg) break;

1452:       for (i = 0; i < n;) {
1453:         /* Now identify the consecutive vals belonging to the same row */
1454:         for (j = i, rstart = row[j]; j < n; j++) {
1455:           if (row[j] != rstart) break;
1456:         }
1457:         if (j < n) ncols = j - i;
1458:         else ncols = n - i;
1459:         /* Now assemble all these values with a single function call */
1460:         PetscCall(MatSetValues_HYPRE(A, 1, row + i, ncols, col + i, val + i, A->insertmode));

1462:         i = j;
1463:       }
1464:     }
1465:     PetscCall(MatStashScatterEnd_Private(&A->stash));
1466:   }

1468:   PetscCallExternal(HYPRE_IJMatrixAssemble, hA->ij);
1469:   /* The assembly routine destroys the aux_matrix, we recreate it here by calling HYPRE_IJMatrixInitialize */
1470:   /* 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 */
1471:   if (!A->sortedfull) {
1472:     hypre_AuxParCSRMatrix *aux_matrix;

1474:     /* call destroy just to make sure we do not leak anything */
1475:     aux_matrix = (hypre_AuxParCSRMatrix *)hypre_IJMatrixTranslator(hA->ij);
1476:     PetscCallExternal(hypre_AuxParCSRMatrixDestroy, aux_matrix);
1477:     hypre_IJMatrixTranslator(hA->ij) = NULL;

1479:     /* Initialize with assembled flag -> it only recreates the aux_par_matrix */
1480:     PetscCallExternal(HYPRE_IJMatrixInitialize, hA->ij);
1481:     aux_matrix = (hypre_AuxParCSRMatrix *)hypre_IJMatrixTranslator(hA->ij);
1482:     if (aux_matrix) {
1483:       hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 1; /* see comment in MatHYPRESetPreallocation_HYPRE */
1484: #if PETSC_PKG_HYPRE_VERSION_LT(2, 19, 0)
1485:       PetscCallExternal(hypre_AuxParCSRMatrixInitialize, aux_matrix);
1486: #else
1487:       PetscCallExternal(hypre_AuxParCSRMatrixInitialize_v2, aux_matrix, HYPRE_MEMORY_HOST);
1488: #endif
1489:     }
1490:   }
1491:   {
1492:     hypre_ParCSRMatrix *parcsr;

1494:     PetscCallExternal(HYPRE_IJMatrixGetObject, hA->ij, (void **)&parcsr);
1495:     if (!hypre_ParCSRMatrixCommPkg(parcsr)) PetscCallExternal(hypre_MatvecCommPkgCreate, parcsr);
1496:   }
1497:   if (!hA->x) PetscCall(VecHYPRE_IJVectorCreate(A->cmap, &hA->x));
1498:   if (!hA->b) PetscCall(VecHYPRE_IJVectorCreate(A->rmap, &hA->b));
1499: #if defined(PETSC_HAVE_HYPRE_DEVICE)
1500:   PetscCall(MatBindToCPU_HYPRE(A, A->boundtocpu));
1501: #endif
1502:   PetscFunctionReturn(PETSC_SUCCESS);
1503: }

1505: static PetscErrorCode MatGetArray_HYPRE(Mat A, PetscInt size, void **array)
1506: {
1507:   Mat_HYPRE *hA = (Mat_HYPRE *)A->data;

1509:   PetscFunctionBegin;
1510:   PetscCheck(hA->array_available, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "Temporary space is in use");

1512:   if (hA->array_size >= size) {
1513:     *array = hA->array;
1514:   } else {
1515:     PetscCall(PetscFree(hA->array));
1516:     hA->array_size = size;
1517:     PetscCall(PetscMalloc(hA->array_size, &hA->array));
1518:     *array = hA->array;
1519:   }

1521:   hA->array_available = PETSC_FALSE;
1522:   PetscFunctionReturn(PETSC_SUCCESS);
1523: }

1525: static PetscErrorCode MatRestoreArray_HYPRE(Mat A, void **array)
1526: {
1527:   Mat_HYPRE *hA = (Mat_HYPRE *)A->data;

1529:   PetscFunctionBegin;
1530:   *array              = NULL;
1531:   hA->array_available = PETSC_TRUE;
1532:   PetscFunctionReturn(PETSC_SUCCESS);
1533: }

1535: static PetscErrorCode MatSetValues_HYPRE(Mat A, PetscInt nr, const PetscInt rows[], PetscInt nc, const PetscInt cols[], const PetscScalar v[], InsertMode ins)
1536: {
1537:   Mat_HYPRE     *hA   = (Mat_HYPRE *)A->data;
1538:   PetscScalar   *vals = (PetscScalar *)v;
1539:   HYPRE_Complex *sscr;
1540:   PetscInt      *cscr[2];
1541:   PetscInt       i, nzc;
1542:   PetscInt       rst = A->rmap->rstart, ren = A->rmap->rend;
1543:   void          *array = NULL;

1545:   PetscFunctionBegin;
1546:   PetscCall(MatGetArray_HYPRE(A, sizeof(PetscInt) * (2 * nc) + sizeof(HYPRE_Complex) * nc * nr, &array));
1547:   cscr[0] = (PetscInt *)array;
1548:   cscr[1] = ((PetscInt *)array) + nc;
1549:   sscr    = (HYPRE_Complex *)(((PetscInt *)array) + nc * 2);
1550:   for (i = 0, nzc = 0; i < nc; i++) {
1551:     if (cols[i] >= 0) {
1552:       cscr[0][nzc]   = cols[i];
1553:       cscr[1][nzc++] = i;
1554:     }
1555:   }
1556:   if (!nzc) {
1557:     PetscCall(MatRestoreArray_HYPRE(A, &array));
1558:     PetscFunctionReturn(PETSC_SUCCESS);
1559:   }

1561: #if 0 //defined(PETSC_HAVE_HYPRE_DEVICE)
1562:   if (HYPRE_MEMORY_HOST != hypre_IJMatrixMemoryLocation(hA->ij)) {
1563:     hypre_ParCSRMatrix *parcsr;

1565:     PetscCallExternal(HYPRE_IJMatrixGetObject,hA->ij,(void**)&parcsr);
1566:     PetscCallExternal(hypre_ParCSRMatrixMigrate,parcsr, HYPRE_MEMORY_HOST);
1567:   }
1568: #endif

1570:   if (ins == ADD_VALUES) {
1571:     for (i = 0; i < nr; i++) {
1572:       if (rows[i] >= 0) {
1573:         PetscInt  j;
1574:         HYPRE_Int hnc = (HYPRE_Int)nzc;

1576:         if (!nzc) continue;
1577:         /* nonlocal values */
1578:         if (rows[i] < rst || rows[i] >= ren) {
1579:           PetscCheck(!A->nooffprocentries, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Setting off process row %" PetscInt_FMT " even though MatSetOption(,MAT_NO_OFF_PROC_ENTRIES,PETSC_TRUE) was set", rows[i]);
1580:           if (hA->donotstash) continue;
1581:         }
1582:         PetscCheck((PetscInt)hnc == nzc, PETSC_COMM_SELF, PETSC_ERR_SUP, "Hypre overflow! number of columns %" PetscInt_FMT " for row %" PetscInt_FMT, nzc, rows[i]);
1583:         for (j = 0; j < nzc; j++) PetscCall(PetscHYPREScalarCast(vals[cscr[1][j]], &sscr[j]));
1584:         PetscCallExternal(HYPRE_IJMatrixAddToValues, hA->ij, 1, &hnc, (HYPRE_BigInt *)(rows + i), (HYPRE_BigInt *)cscr[0], sscr);
1585:       }
1586:       vals += nc;
1587:     }
1588:   } else { /* INSERT_VALUES */
1589:     for (i = 0; i < nr; i++) {
1590:       if (rows[i] >= 0) {
1591:         PetscInt  j;
1592:         HYPRE_Int hnc = (HYPRE_Int)nzc;

1594:         if (!nzc) continue;
1595:         PetscCheck((PetscInt)hnc == nzc, PETSC_COMM_SELF, PETSC_ERR_SUP, "Hypre overflow! number of columns %" PetscInt_FMT " for row %" PetscInt_FMT, nzc, rows[i]);
1596:         for (j = 0; j < nzc; j++) PetscCall(PetscHYPREScalarCast(vals[cscr[1][j]], &sscr[j]));
1597:         /* nonlocal values */
1598:         if (rows[i] < rst || rows[i] >= ren) {
1599:           PetscCheck(!A->nooffprocentries, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Setting off process row %" PetscInt_FMT " even though MatSetOption(,MAT_NO_OFF_PROC_ENTRIES,PETSC_TRUE) was set", rows[i]);
1600:           if (!hA->donotstash) PetscCall(MatStashValuesRow_Private(&A->stash, rows[i], nzc, cscr[0], (PetscScalar *)sscr, PETSC_FALSE));
1601:         }
1602:         /* local values */
1603:         else
1604:           PetscCallExternal(HYPRE_IJMatrixSetValues, hA->ij, 1, &hnc, (HYPRE_BigInt *)(rows + i), (HYPRE_BigInt *)cscr[0], sscr);
1605:       }
1606:       vals += nc;
1607:     }
1608:   }

1610:   PetscCall(MatRestoreArray_HYPRE(A, &array));
1611:   PetscFunctionReturn(PETSC_SUCCESS);
1612: }

1614: static PetscErrorCode MatHYPRESetPreallocation_HYPRE(Mat A, PetscInt dnz, const PetscInt dnnz[], PetscInt onz, const PetscInt onnz[])
1615: {
1616:   Mat_HYPRE  *hA = (Mat_HYPRE *)A->data;
1617:   HYPRE_Int  *hdnnz, *honnz;
1618:   PetscInt    i, rs, re, cs, ce, bs;
1619:   PetscMPIInt size;

1621:   PetscFunctionBegin;
1622:   PetscCall(PetscLayoutSetUp(A->rmap));
1623:   PetscCall(PetscLayoutSetUp(A->cmap));
1624:   rs = A->rmap->rstart;
1625:   re = A->rmap->rend;
1626:   cs = A->cmap->rstart;
1627:   ce = A->cmap->rend;
1628:   if (!hA->ij) {
1629:     PetscCallExternal(HYPRE_IJMatrixCreate, hA->comm, rs, re - 1, cs, ce - 1, &hA->ij);
1630:     PetscCallExternal(HYPRE_IJMatrixSetObjectType, hA->ij, HYPRE_PARCSR);
1631:   } else {
1632:     HYPRE_BigInt hrs, hre, hcs, hce;
1633:     PetscCallExternal(HYPRE_IJMatrixGetLocalRange, hA->ij, &hrs, &hre, &hcs, &hce);
1634:     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);
1635:     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);
1636:   }
1637:   PetscCall(MatHYPRE_DestroyCOOMat(A));
1638:   PetscCall(MatGetBlockSize(A, &bs));
1639:   if (dnz == PETSC_DEFAULT || dnz == PETSC_DECIDE) dnz = 10 * bs;
1640:   if (onz == PETSC_DEFAULT || onz == PETSC_DECIDE) onz = 10 * bs;

1642:   if (!dnnz) {
1643:     PetscCall(PetscMalloc1(A->rmap->n, &hdnnz));
1644:     for (i = 0; i < A->rmap->n; i++) hdnnz[i] = dnz;
1645:   } else {
1646:     hdnnz = (HYPRE_Int *)dnnz;
1647:   }
1648:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
1649:   if (size > 1) {
1650:     hypre_AuxParCSRMatrix *aux_matrix;
1651:     if (!onnz) {
1652:       PetscCall(PetscMalloc1(A->rmap->n, &honnz));
1653:       for (i = 0; i < A->rmap->n; i++) honnz[i] = onz;
1654:     } else honnz = (HYPRE_Int *)onnz;
1655:     /* SetDiagOffdSizes sets hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 0, since it seems
1656:        they assume the user will input the entire row values, properly sorted
1657:        In PETSc, we don't make such an assumption and set this flag to 1,
1658:        unless the option MAT_SORTED_FULL is set to true.
1659:        Also, to avoid possible memory leaks, we destroy and recreate the translator
1660:        This has to be done here, as HYPRE_IJMatrixInitialize will properly initialize
1661:        the IJ matrix for us */
1662:     aux_matrix = (hypre_AuxParCSRMatrix *)hypre_IJMatrixTranslator(hA->ij);
1663:     hypre_AuxParCSRMatrixDestroy(aux_matrix);
1664:     hypre_IJMatrixTranslator(hA->ij) = NULL;
1665:     PetscCallExternal(HYPRE_IJMatrixSetDiagOffdSizes, hA->ij, hdnnz, honnz);
1666:     aux_matrix                               = (hypre_AuxParCSRMatrix *)hypre_IJMatrixTranslator(hA->ij);
1667:     hypre_AuxParCSRMatrixNeedAux(aux_matrix) = !A->sortedfull;
1668:   } else {
1669:     honnz = NULL;
1670:     PetscCallExternal(HYPRE_IJMatrixSetRowSizes, hA->ij, hdnnz);
1671:   }

1673:   /* reset assembled flag and call the initialize method */
1674:   hypre_IJMatrixAssembleFlag(hA->ij) = 0;
1675: #if PETSC_PKG_HYPRE_VERSION_LT(2, 19, 0)
1676:   PetscCallExternal(HYPRE_IJMatrixInitialize, hA->ij);
1677: #else
1678:   PetscCallExternal(HYPRE_IJMatrixInitialize_v2, hA->ij, HYPRE_MEMORY_HOST);
1679: #endif
1680:   if (!dnnz) PetscCall(PetscFree(hdnnz));
1681:   if (!onnz && honnz) PetscCall(PetscFree(honnz));
1682:   /* Match AIJ logic */
1683:   A->preallocated = PETSC_TRUE;
1684:   A->assembled    = PETSC_FALSE;
1685:   PetscFunctionReturn(PETSC_SUCCESS);
1686: }

1688: /*@C
1689:   MatHYPRESetPreallocation - Preallocates memory for a sparse parallel matrix in HYPRE IJ format

1691:   Collective

1693:   Input Parameters:
1694: + A    - the matrix
1695: . dnz  - number of nonzeros per row in DIAGONAL portion of local submatrix
1696:           (same value is used for all local rows)
1697: . dnnz - array containing the number of nonzeros in the various rows of the
1698:           DIAGONAL portion of the local submatrix (possibly different for each row)
1699:           or `NULL` (`PETSC_NULL_INTEGER` in Fortran), if `d_nz` is used to specify the nonzero structure.
1700:           The size of this array is equal to the number of local rows, i.e `m`.
1701:           For matrices that will be factored, you must leave room for (and set)
1702:           the diagonal entry even if it is zero.
1703: . onz  - number of nonzeros per row in the OFF-DIAGONAL portion of local
1704:           submatrix (same value is used for all local rows).
1705: - onnz - array containing the number of nonzeros in the various rows of the
1706:           OFF-DIAGONAL portion of the local submatrix (possibly different for
1707:           each row) or `NULL` (`PETSC_NULL_INTEGER` in Fortran), if `o_nz` is used to specify the nonzero
1708:           structure. The size of this array is equal to the number
1709:           of local rows, i.e `m`.

1711:   Level: intermediate

1713:   Note:
1714:   If the *nnz parameter is given then the *nz parameter is ignored; for sequential matrices, `onz` and `onnz` are ignored.

1716: .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MatMPIAIJSetPreallocation()`, `MATHYPRE`, `MATAIJ`
1717: @*/
1718: PetscErrorCode MatHYPRESetPreallocation(Mat A, PetscInt dnz, const PetscInt dnnz[], PetscInt onz, const PetscInt onnz[])
1719: {
1720:   PetscFunctionBegin;
1723:   PetscTryMethod(A, "MatHYPRESetPreallocation_C", (Mat, PetscInt, const PetscInt[], PetscInt, const PetscInt[]), (A, dnz, dnnz, onz, onnz));
1724:   PetscFunctionReturn(PETSC_SUCCESS);
1725: }

1727: /*@C
1728:   MatCreateFromParCSR - Creates a `Mat` from a `hypre_ParCSRMatrix`

1730:   Collective

1732:   Input Parameters:
1733: + parcsr   - the pointer to the `hypre_ParCSRMatrix`
1734: . mtype    - matrix type to be created. Currently `MATAIJ`, `MATIS` and `MATHYPRE` are supported.
1735: - copymode - PETSc copying options, see  `PetscCopyMode`

1737:   Output Parameter:
1738: . A - the matrix

1740:   Level: intermediate

1742: .seealso: [](ch_matrices), `Mat`, `MATHYPRE`, `PetscCopyMode`
1743: @*/
1744: PETSC_EXTERN PetscErrorCode MatCreateFromParCSR(hypre_ParCSRMatrix *parcsr, MatType mtype, PetscCopyMode copymode, Mat *A)
1745: {
1746:   Mat        T;
1747:   Mat_HYPRE *hA;
1748:   MPI_Comm   comm;
1749:   PetscInt   rstart, rend, cstart, cend, M, N;
1750:   PetscBool  isseqaij, isseqaijmkl, ismpiaij, isaij, ishyp, isis;

1752:   PetscFunctionBegin;
1753:   comm = hypre_ParCSRMatrixComm(parcsr);
1754:   PetscCall(PetscStrcmp(mtype, MATSEQAIJ, &isseqaij));
1755:   PetscCall(PetscStrcmp(mtype, MATSEQAIJMKL, &isseqaijmkl));
1756:   PetscCall(PetscStrcmp(mtype, MATMPIAIJ, &ismpiaij));
1757:   PetscCall(PetscStrcmp(mtype, MATAIJ, &isaij));
1758:   PetscCall(PetscStrcmp(mtype, MATHYPRE, &ishyp));
1759:   PetscCall(PetscStrcmp(mtype, MATIS, &isis));
1760:   isaij = (PetscBool)(isseqaij || isseqaijmkl || ismpiaij || isaij);
1761:   /* TODO */
1762:   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);
1763:   /* access ParCSRMatrix */
1764:   rstart = hypre_ParCSRMatrixFirstRowIndex(parcsr);
1765:   rend   = hypre_ParCSRMatrixLastRowIndex(parcsr);
1766:   cstart = hypre_ParCSRMatrixFirstColDiag(parcsr);
1767:   cend   = hypre_ParCSRMatrixLastColDiag(parcsr);
1768:   M      = hypre_ParCSRMatrixGlobalNumRows(parcsr);
1769:   N      = hypre_ParCSRMatrixGlobalNumCols(parcsr);

1771:   /* create PETSc matrix with MatHYPRE */
1772:   PetscCall(MatCreate(comm, &T));
1773:   PetscCall(MatSetSizes(T, PetscMax(rend - rstart + 1, 0), PetscMax(cend - cstart + 1, 0), M, N));
1774:   PetscCall(MatSetType(T, MATHYPRE));
1775:   hA = (Mat_HYPRE *)T->data;

1777:   /* create HYPRE_IJMatrix */
1778:   PetscCallExternal(HYPRE_IJMatrixCreate, hA->comm, rstart, rend, cstart, cend, &hA->ij);
1779:   PetscCallExternal(HYPRE_IJMatrixSetObjectType, hA->ij, HYPRE_PARCSR);

1781:   /* create new ParCSR object if needed */
1782:   if (ishyp && copymode == PETSC_COPY_VALUES) {
1783:     hypre_ParCSRMatrix *new_parcsr;
1784: #if PETSC_PKG_HYPRE_VERSION_LT(2, 18, 0)
1785:     hypre_CSRMatrix *hdiag, *hoffd, *ndiag, *noffd;

1787:     new_parcsr = hypre_ParCSRMatrixClone(parcsr, 0);
1788:     hdiag      = hypre_ParCSRMatrixDiag(parcsr);
1789:     hoffd      = hypre_ParCSRMatrixOffd(parcsr);
1790:     ndiag      = hypre_ParCSRMatrixDiag(new_parcsr);
1791:     noffd      = hypre_ParCSRMatrixOffd(new_parcsr);
1792:     PetscCall(PetscArraycpy(hypre_CSRMatrixData(ndiag), hypre_CSRMatrixData(hdiag), hypre_CSRMatrixNumNonzeros(hdiag)));
1793:     PetscCall(PetscArraycpy(hypre_CSRMatrixData(noffd), hypre_CSRMatrixData(hoffd), hypre_CSRMatrixNumNonzeros(hoffd)));
1794: #else
1795:     new_parcsr = hypre_ParCSRMatrixClone(parcsr, 1);
1796: #endif
1797:     parcsr   = new_parcsr;
1798:     copymode = PETSC_OWN_POINTER;
1799:   }

1801:   /* set ParCSR object */
1802:   hypre_IJMatrixObject(hA->ij) = parcsr;
1803:   T->preallocated              = PETSC_TRUE;

1805:   /* set assembled flag */
1806:   hypre_IJMatrixAssembleFlag(hA->ij) = 1;
1807: #if 0
1808:   PetscCallExternal(HYPRE_IJMatrixInitialize,hA->ij);
1809: #endif
1810:   if (ishyp) {
1811:     PetscMPIInt myid = 0;

1813:     /* make sure we always have row_starts and col_starts available */
1814:     if (HYPRE_AssumedPartitionCheck()) PetscCallMPI(MPI_Comm_rank(comm, &myid));
1815: #if defined(hypre_ParCSRMatrixOwnsRowStarts)
1816:     if (!hypre_ParCSRMatrixOwnsColStarts(parcsr)) {
1817:       PetscLayout map;

1819:       PetscCall(MatGetLayouts(T, NULL, &map));
1820:       PetscCall(PetscLayoutSetUp(map));
1821:       hypre_ParCSRMatrixColStarts(parcsr) = (HYPRE_BigInt *)(map->range + myid);
1822:     }
1823:     if (!hypre_ParCSRMatrixOwnsRowStarts(parcsr)) {
1824:       PetscLayout map;

1826:       PetscCall(MatGetLayouts(T, &map, NULL));
1827:       PetscCall(PetscLayoutSetUp(map));
1828:       hypre_ParCSRMatrixRowStarts(parcsr) = (HYPRE_BigInt *)(map->range + myid);
1829:     }
1830: #endif
1831:     /* prevent from freeing the pointer */
1832:     if (copymode == PETSC_USE_POINTER) hA->inner_free = PETSC_FALSE;
1833:     *A = T;
1834:     PetscCall(MatSetOption(*A, MAT_SORTED_FULL, PETSC_TRUE));
1835:     PetscCall(MatAssemblyBegin(*A, MAT_FINAL_ASSEMBLY));
1836:     PetscCall(MatAssemblyEnd(*A, MAT_FINAL_ASSEMBLY));
1837:   } else if (isaij) {
1838:     if (copymode != PETSC_OWN_POINTER) {
1839:       /* prevent from freeing the pointer */
1840:       hA->inner_free = PETSC_FALSE;
1841:       PetscCall(MatConvert_HYPRE_AIJ(T, MATAIJ, MAT_INITIAL_MATRIX, A));
1842:       PetscCall(MatDestroy(&T));
1843:     } else { /* AIJ return type with PETSC_OWN_POINTER */
1844:       PetscCall(MatConvert_HYPRE_AIJ(T, MATAIJ, MAT_INPLACE_MATRIX, &T));
1845:       *A = T;
1846:     }
1847:   } else if (isis) {
1848:     PetscCall(MatConvert_HYPRE_IS(T, MATIS, MAT_INITIAL_MATRIX, A));
1849:     if (copymode != PETSC_OWN_POINTER) hA->inner_free = PETSC_FALSE;
1850:     PetscCall(MatDestroy(&T));
1851:   }
1852:   PetscFunctionReturn(PETSC_SUCCESS);
1853: }

1855: static PetscErrorCode MatHYPREGetParCSR_HYPRE(Mat A, hypre_ParCSRMatrix **parcsr)
1856: {
1857:   Mat_HYPRE *hA = (Mat_HYPRE *)A->data;
1858:   HYPRE_Int  type;

1860:   PetscFunctionBegin;
1861:   PetscCheck(hA->ij, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "HYPRE_IJMatrix not present");
1862:   PetscCallExternal(HYPRE_IJMatrixGetObjectType, hA->ij, &type);
1863:   PetscCheck(type == HYPRE_PARCSR, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "HYPRE_IJMatrix is not of type HYPRE_PARCSR");
1864:   PetscCallExternal(HYPRE_IJMatrixGetObject, hA->ij, (void **)parcsr);
1865:   PetscFunctionReturn(PETSC_SUCCESS);
1866: }

1868: /*@C
1869:   MatHYPREGetParCSR - Gets the pointer to the ParCSR matrix

1871:   Not Collective, No Fortran Support

1873:   Input Parameter:
1874: . A - the `MATHYPRE` object

1876:   Output Parameter:
1877: . parcsr - the pointer to the `hypre_ParCSRMatrix`

1879:   Level: intermediate

1881: .seealso: [](ch_matrices), `Mat`, `MATHYPRE`, `PetscCopyMode`
1882: @*/
1883: PetscErrorCode MatHYPREGetParCSR(Mat A, hypre_ParCSRMatrix **parcsr)
1884: {
1885:   PetscFunctionBegin;
1888:   PetscUseMethod(A, "MatHYPREGetParCSR_C", (Mat, hypre_ParCSRMatrix **), (A, parcsr));
1889:   PetscFunctionReturn(PETSC_SUCCESS);
1890: }

1892: static PetscErrorCode MatMissingDiagonal_HYPRE(Mat A, PetscBool *missing, PetscInt *dd)
1893: {
1894:   hypre_ParCSRMatrix *parcsr;
1895:   hypre_CSRMatrix    *ha;
1896:   PetscInt            rst;

1898:   PetscFunctionBegin;
1899:   PetscCheck(A->rmap->n == A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not implemented with non-square diagonal blocks");
1900:   PetscCall(MatGetOwnershipRange(A, &rst, NULL));
1901:   PetscCall(MatHYPREGetParCSR_HYPRE(A, &parcsr));
1902:   if (missing) *missing = PETSC_FALSE;
1903:   if (dd) *dd = -1;
1904:   ha = hypre_ParCSRMatrixDiag(parcsr);
1905:   if (ha) {
1906:     PetscInt   size, i;
1907:     HYPRE_Int *ii, *jj;

1909:     size = hypre_CSRMatrixNumRows(ha);
1910:     ii   = hypre_CSRMatrixI(ha);
1911:     jj   = hypre_CSRMatrixJ(ha);
1912:     for (i = 0; i < size; i++) {
1913:       PetscInt  j;
1914:       PetscBool found = PETSC_FALSE;

1916:       for (j = ii[i]; j < ii[i + 1] && !found; j++) found = (jj[j] == i) ? PETSC_TRUE : PETSC_FALSE;

1918:       if (!found) {
1919:         PetscCall(PetscInfo(A, "Matrix is missing local diagonal entry %" PetscInt_FMT "\n", i));
1920:         if (missing) *missing = PETSC_TRUE;
1921:         if (dd) *dd = i + rst;
1922:         PetscFunctionReturn(PETSC_SUCCESS);
1923:       }
1924:     }
1925:     if (!size) {
1926:       PetscCall(PetscInfo(A, "Matrix has no diagonal entries therefore is missing diagonal\n"));
1927:       if (missing) *missing = PETSC_TRUE;
1928:       if (dd) *dd = rst;
1929:     }
1930:   } else {
1931:     PetscCall(PetscInfo(A, "Matrix has no diagonal entries therefore is missing diagonal\n"));
1932:     if (missing) *missing = PETSC_TRUE;
1933:     if (dd) *dd = rst;
1934:   }
1935:   PetscFunctionReturn(PETSC_SUCCESS);
1936: }

1938: static PetscErrorCode MatScale_HYPRE(Mat A, PetscScalar s)
1939: {
1940:   hypre_ParCSRMatrix *parcsr;
1941: #if PETSC_PKG_HYPRE_VERSION_LT(2, 19, 0)
1942:   hypre_CSRMatrix *ha;
1943: #endif
1944:   HYPRE_Complex hs;

1946:   PetscFunctionBegin;
1947:   PetscCall(PetscHYPREScalarCast(s, &hs));
1948:   PetscCall(MatHYPREGetParCSR_HYPRE(A, &parcsr));
1949: #if PETSC_PKG_HYPRE_VERSION_GE(2, 19, 0)
1950:   PetscCallExternal(hypre_ParCSRMatrixScale, parcsr, hs);
1951: #else /* diagonal part */
1952:   ha = hypre_ParCSRMatrixDiag(parcsr);
1953:   if (ha) {
1954:     PetscInt       size, i;
1955:     HYPRE_Int     *ii;
1956:     HYPRE_Complex *a;

1958:     size = hypre_CSRMatrixNumRows(ha);
1959:     a    = hypre_CSRMatrixData(ha);
1960:     ii   = hypre_CSRMatrixI(ha);
1961:     for (i = 0; i < ii[size]; i++) a[i] *= hs;
1962:   }
1963:   /* off-diagonal part */
1964:   ha = hypre_ParCSRMatrixOffd(parcsr);
1965:   if (ha) {
1966:     PetscInt       size, i;
1967:     HYPRE_Int     *ii;
1968:     HYPRE_Complex *a;

1970:     size = hypre_CSRMatrixNumRows(ha);
1971:     a    = hypre_CSRMatrixData(ha);
1972:     ii   = hypre_CSRMatrixI(ha);
1973:     for (i = 0; i < ii[size]; i++) a[i] *= hs;
1974:   }
1975: #endif
1976:   PetscFunctionReturn(PETSC_SUCCESS);
1977: }

1979: static PetscErrorCode MatZeroRowsColumns_HYPRE(Mat A, PetscInt numRows, const PetscInt rows[], PetscScalar diag, Vec x, Vec b)
1980: {
1981:   hypre_ParCSRMatrix *parcsr;
1982:   HYPRE_Int          *lrows;
1983:   PetscInt            rst, ren, i;

1985:   PetscFunctionBegin;
1986:   PetscCheck(!x && !b, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "To be implemented");
1987:   PetscCall(MatHYPREGetParCSR_HYPRE(A, &parcsr));
1988:   PetscCall(PetscMalloc1(numRows, &lrows));
1989:   PetscCall(MatGetOwnershipRange(A, &rst, &ren));
1990:   for (i = 0; i < numRows; i++) {
1991:     PetscCheck(rows[i] >= rst && rows[i] < ren, PETSC_COMM_SELF, PETSC_ERR_SUP, "Non-local rows not yet supported");
1992:     lrows[i] = rows[i] - rst;
1993:   }
1994:   PetscCallExternal(hypre_ParCSRMatrixEliminateRowsCols, parcsr, numRows, lrows);
1995:   PetscCall(PetscFree(lrows));
1996:   PetscFunctionReturn(PETSC_SUCCESS);
1997: }

1999: static PetscErrorCode MatZeroEntries_HYPRE_CSRMatrix(hypre_CSRMatrix *ha)
2000: {
2001:   PetscFunctionBegin;
2002:   if (ha) {
2003:     HYPRE_Int     *ii, size;
2004:     HYPRE_Complex *a;

2006:     size = hypre_CSRMatrixNumRows(ha);
2007:     a    = hypre_CSRMatrixData(ha);
2008:     ii   = hypre_CSRMatrixI(ha);

2010:     if (a) PetscCall(PetscArrayzero(a, ii[size]));
2011:   }
2012:   PetscFunctionReturn(PETSC_SUCCESS);
2013: }

2015: static PetscErrorCode MatZeroEntries_HYPRE(Mat A)
2016: {
2017:   Mat_HYPRE *hA = (Mat_HYPRE *)A->data;

2019:   PetscFunctionBegin;
2020:   if (HYPRE_MEMORY_DEVICE == hypre_IJMatrixMemoryLocation(hA->ij)) {
2021:     PetscCallExternal(HYPRE_IJMatrixSetConstantValues, hA->ij, 0.0);
2022:   } else {
2023:     hypre_ParCSRMatrix *parcsr;

2025:     PetscCall(MatHYPREGetParCSR_HYPRE(A, &parcsr));
2026:     PetscCall(MatZeroEntries_HYPRE_CSRMatrix(hypre_ParCSRMatrixDiag(parcsr)));
2027:     PetscCall(MatZeroEntries_HYPRE_CSRMatrix(hypre_ParCSRMatrixOffd(parcsr)));
2028:   }
2029:   PetscFunctionReturn(PETSC_SUCCESS);
2030: }

2032: static PetscErrorCode MatZeroRows_HYPRE_CSRMatrix(hypre_CSRMatrix *hA, PetscInt N, const PetscInt rows[], HYPRE_Complex diag)
2033: {
2034:   PetscInt       ii;
2035:   HYPRE_Int     *i, *j;
2036:   HYPRE_Complex *a;

2038:   PetscFunctionBegin;
2039:   if (!hA) PetscFunctionReturn(PETSC_SUCCESS);

2041:   i = hypre_CSRMatrixI(hA);
2042:   j = hypre_CSRMatrixJ(hA);
2043:   a = hypre_CSRMatrixData(hA);
2044: #if defined(PETSC_HAVE_HYPRE_DEVICE)
2045:   if (HYPRE_MEMORY_DEVICE == hypre_CSRMatrixMemoryLocation(hA)) {
2046:   #if defined(HYPRE_USING_CUDA)
2047:     MatZeroRows_CUDA(N, rows, i, j, a, diag);
2048:   #elif defined(HYPRE_USING_HIP)
2049:     MatZeroRows_HIP(N, rows, i, j, a, diag);
2050:   #elif defined(PETSC_HAVE_KOKKOS)
2051:     MatZeroRows_Kokkos(N, rows, i, j, a, diag);
2052:   #else
2053:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for MatZeroRows on a hypre matrix in this memory location");
2054:   #endif
2055:   } else
2056: #endif
2057:   {
2058:     for (ii = 0; ii < N; ii++) {
2059:       HYPRE_Int jj, ibeg, iend, irow;

2061:       irow = rows[ii];
2062:       ibeg = i[irow];
2063:       iend = i[irow + 1];
2064:       for (jj = ibeg; jj < iend; jj++)
2065:         if (j[jj] == irow) a[jj] = diag;
2066:         else a[jj] = 0.0;
2067:     }
2068:   }
2069:   PetscFunctionReturn(PETSC_SUCCESS);
2070: }

2072: static PetscErrorCode MatZeroRows_HYPRE(Mat A, PetscInt N, const PetscInt rows[], PetscScalar diag, Vec x, Vec b)
2073: {
2074:   hypre_ParCSRMatrix *parcsr;
2075:   PetscInt           *lrows, len, *lrows2;
2076:   HYPRE_Complex       hdiag;

2078:   PetscFunctionBegin;
2079:   PetscCheck(!x && !b, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Does not support to modify the solution and the right hand size");
2080:   PetscCall(PetscHYPREScalarCast(diag, &hdiag));
2081:   /* retrieve the internal matrix */
2082:   PetscCall(MatHYPREGetParCSR_HYPRE(A, &parcsr));
2083:   /* get locally owned rows */
2084:   PetscCall(MatZeroRowsMapLocal_Private(A, N, rows, &len, &lrows));

2086: #if defined(PETSC_HAVE_HYPRE_DEVICE)
2087:   if (HYPRE_MEMORY_DEVICE == hypre_CSRMatrixMemoryLocation(hypre_ParCSRMatrixDiag(parcsr))) {
2088:     Mat_HYPRE *hA = (Mat_HYPRE *)A->data;
2089:     PetscInt   m;
2090:     PetscCall(MatGetLocalSize(A, &m, NULL));
2091:     if (!hA->rows_d) {
2092:       hA->rows_d = hypre_TAlloc(PetscInt, m, HYPRE_MEMORY_DEVICE);
2093:       if (m) PetscCheck(hA->rows_d, PETSC_COMM_SELF, PETSC_ERR_MEM, "HYPRE_TAlloc failed");
2094:     }
2095:     PetscCheck(len <= m, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Too many rows in rows[]");
2096:     PetscStackCallExternalVoid("hypre_Memcpy", hypre_Memcpy(hA->rows_d, lrows, sizeof(PetscInt) * len, HYPRE_MEMORY_DEVICE, HYPRE_MEMORY_HOST));
2097:     lrows2 = hA->rows_d;
2098:   } else
2099: #endif
2100:   {
2101:     lrows2 = lrows;
2102:   }

2104:   /* zero diagonal part */
2105:   PetscCall(MatZeroRows_HYPRE_CSRMatrix(hypre_ParCSRMatrixDiag(parcsr), len, lrows2, hdiag));
2106:   /* zero off-diagonal part */
2107:   PetscCall(MatZeroRows_HYPRE_CSRMatrix(hypre_ParCSRMatrixOffd(parcsr), len, lrows2, 0.0));

2109:   PetscCall(PetscFree(lrows));
2110:   PetscFunctionReturn(PETSC_SUCCESS);
2111: }

2113: static PetscErrorCode MatAssemblyBegin_HYPRE(Mat mat, MatAssemblyType mode)
2114: {
2115:   PetscFunctionBegin;
2116:   if (mat->nooffprocentries) PetscFunctionReturn(PETSC_SUCCESS);

2118:   PetscCall(MatStashScatterBegin_Private(mat, &mat->stash, mat->rmap->range));
2119:   PetscFunctionReturn(PETSC_SUCCESS);
2120: }

2122: static PetscErrorCode MatGetRow_HYPRE(Mat A, PetscInt row, PetscInt *nz, PetscInt **idx, PetscScalar **v)
2123: {
2124:   hypre_ParCSRMatrix *parcsr;
2125:   HYPRE_Int           hnz;

2127:   PetscFunctionBegin;
2128:   /* retrieve the internal matrix */
2129:   PetscCall(MatHYPREGetParCSR_HYPRE(A, &parcsr));
2130:   /* call HYPRE API */
2131:   PetscCallExternal(HYPRE_ParCSRMatrixGetRow, parcsr, row, &hnz, (HYPRE_BigInt **)idx, (HYPRE_Complex **)v);
2132:   if (nz) *nz = (PetscInt)hnz;
2133:   PetscFunctionReturn(PETSC_SUCCESS);
2134: }

2136: static PetscErrorCode MatRestoreRow_HYPRE(Mat A, PetscInt row, PetscInt *nz, PetscInt **idx, PetscScalar **v)
2137: {
2138:   hypre_ParCSRMatrix *parcsr;
2139:   HYPRE_Int           hnz;

2141:   PetscFunctionBegin;
2142:   /* retrieve the internal matrix */
2143:   PetscCall(MatHYPREGetParCSR_HYPRE(A, &parcsr));
2144:   /* call HYPRE API */
2145:   hnz = nz ? (HYPRE_Int)(*nz) : 0;
2146:   PetscCallExternal(HYPRE_ParCSRMatrixRestoreRow, parcsr, row, &hnz, (HYPRE_BigInt **)idx, (HYPRE_Complex **)v);
2147:   PetscFunctionReturn(PETSC_SUCCESS);
2148: }

2150: static PetscErrorCode MatGetValues_HYPRE(Mat A, PetscInt m, const PetscInt idxm[], PetscInt n, const PetscInt idxn[], PetscScalar v[])
2151: {
2152:   Mat_HYPRE *hA = (Mat_HYPRE *)A->data;
2153:   PetscInt   i;

2155:   PetscFunctionBegin;
2156:   if (!m || !n) PetscFunctionReturn(PETSC_SUCCESS);
2157:   /* Ignore negative row indices
2158:    * And negative column indices should be automatically ignored in hypre
2159:    * */
2160:   for (i = 0; i < m; i++) {
2161:     if (idxm[i] >= 0) {
2162:       HYPRE_Int hn = (HYPRE_Int)n;
2163:       PetscCallExternal(HYPRE_IJMatrixGetValues, hA->ij, 1, &hn, (HYPRE_BigInt *)&idxm[i], (HYPRE_BigInt *)idxn, (HYPRE_Complex *)(v + i * n));
2164:     }
2165:   }
2166:   PetscFunctionReturn(PETSC_SUCCESS);
2167: }

2169: static PetscErrorCode MatSetOption_HYPRE(Mat A, MatOption op, PetscBool flg)
2170: {
2171:   Mat_HYPRE *hA = (Mat_HYPRE *)A->data;

2173:   PetscFunctionBegin;
2174:   switch (op) {
2175:   case MAT_NO_OFF_PROC_ENTRIES:
2176:     if (flg) PetscCallExternal(HYPRE_IJMatrixSetMaxOffProcElmts, hA->ij, 0);
2177:     break;
2178:   case MAT_IGNORE_OFF_PROC_ENTRIES:
2179:     hA->donotstash = flg;
2180:     break;
2181:   default:
2182:     break;
2183:   }
2184:   PetscFunctionReturn(PETSC_SUCCESS);
2185: }

2187: static PetscErrorCode MatView_HYPRE(Mat A, PetscViewer view)
2188: {
2189:   PetscViewerFormat format;

2191:   PetscFunctionBegin;
2192:   PetscCall(PetscViewerGetFormat(view, &format));
2193:   if (format == PETSC_VIEWER_ASCII_FACTOR_INFO || format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) PetscFunctionReturn(PETSC_SUCCESS);
2194:   if (format != PETSC_VIEWER_NATIVE) {
2195:     Mat                 B;
2196:     hypre_ParCSRMatrix *parcsr;
2197:     PetscErrorCode (*mview)(Mat, PetscViewer) = NULL;

2199:     PetscCall(MatHYPREGetParCSR_HYPRE(A, &parcsr));
2200:     PetscCall(MatCreateFromParCSR(parcsr, MATAIJ, PETSC_USE_POINTER, &B));
2201:     PetscCall(MatGetOperation(B, MATOP_VIEW, (PetscErrorCodeFn **)&mview));
2202:     PetscCheck(mview, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing view operation");
2203:     PetscCall((*mview)(B, view));
2204:     PetscCall(MatDestroy(&B));
2205:   } else {
2206:     Mat_HYPRE  *hA = (Mat_HYPRE *)A->data;
2207:     PetscMPIInt size;
2208:     PetscBool   isascii;
2209:     const char *filename;

2211:     /* HYPRE uses only text files */
2212:     PetscCall(PetscObjectTypeCompare((PetscObject)view, PETSCVIEWERASCII, &isascii));
2213:     PetscCheck(isascii, PetscObjectComm((PetscObject)view), PETSC_ERR_SUP, "PetscViewerType %s: native HYPRE format needs PETSCVIEWERASCII", ((PetscObject)view)->type_name);
2214:     PetscCall(PetscViewerFileGetName(view, &filename));
2215:     PetscCallExternal(HYPRE_IJMatrixPrint, hA->ij, filename);
2216:     PetscCallMPI(MPI_Comm_size(hA->comm, &size));
2217:     if (size > 1) {
2218:       PetscCall(PetscViewerASCIIPrintf(view, "Matrix files: %s.%05d ... %s.%05d\n", filename, 0, filename, size - 1));
2219:     } else {
2220:       PetscCall(PetscViewerASCIIPrintf(view, "Matrix file: %s.%05d\n", filename, 0));
2221:     }
2222:   }
2223:   PetscFunctionReturn(PETSC_SUCCESS);
2224: }

2226: static PetscErrorCode MatCopy_HYPRE(Mat A, Mat B, MatStructure str)
2227: {
2228:   hypre_ParCSRMatrix *acsr, *bcsr;

2230:   PetscFunctionBegin;
2231:   if (str == SAME_NONZERO_PATTERN && A->ops->copy == B->ops->copy) {
2232:     PetscCall(MatHYPREGetParCSR_HYPRE(A, &acsr));
2233:     PetscCall(MatHYPREGetParCSR_HYPRE(B, &bcsr));
2234:     PetscCallExternal(hypre_ParCSRMatrixCopy, acsr, bcsr, 1);
2235:     PetscCall(MatSetOption(B, MAT_SORTED_FULL, PETSC_TRUE)); /* "perfect" preallocation, so no need for hypre_AuxParCSRMatrixNeedAux */
2236:     PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
2237:     PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
2238:   } else {
2239:     PetscCall(MatCopy_Basic(A, B, str));
2240:   }
2241:   PetscFunctionReturn(PETSC_SUCCESS);
2242: }

2244: static PetscErrorCode MatGetDiagonal_HYPRE(Mat A, Vec d)
2245: {
2246:   hypre_ParCSRMatrix *parcsr;
2247:   hypre_CSRMatrix    *dmat;
2248:   HYPRE_Complex      *a;
2249:   PetscBool           cong;

2251:   PetscFunctionBegin;
2252:   PetscCall(MatHasCongruentLayouts(A, &cong));
2253:   PetscCheck(cong, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Only for square matrices with same local distributions of rows and columns");
2254:   PetscCall(MatHYPREGetParCSR_HYPRE(A, &parcsr));
2255:   dmat = hypre_ParCSRMatrixDiag(parcsr);
2256:   if (dmat) {
2257: #if defined(PETSC_HAVE_HYPRE_DEVICE)
2258:     HYPRE_MemoryLocation mem = hypre_CSRMatrixMemoryLocation(dmat);
2259: #else
2260:     HYPRE_MemoryLocation mem = HYPRE_MEMORY_HOST;
2261: #endif

2263:     if (mem != HYPRE_MEMORY_HOST) PetscCall(VecGetArrayWriteAndMemType(d, (PetscScalar **)&a, NULL));
2264:     else PetscCall(VecGetArrayWrite(d, (PetscScalar **)&a));
2265:     hypre_CSRMatrixExtractDiagonal(dmat, a, 0);
2266:     if (mem != HYPRE_MEMORY_HOST) PetscCall(VecRestoreArrayWriteAndMemType(d, (PetscScalar **)&a));
2267:     else PetscCall(VecRestoreArrayWrite(d, (PetscScalar **)&a));
2268:   }
2269:   PetscFunctionReturn(PETSC_SUCCESS);
2270: }

2272: #include <petscblaslapack.h>

2274: static PetscErrorCode MatAXPY_HYPRE(Mat Y, PetscScalar a, Mat X, MatStructure str)
2275: {
2276:   PetscFunctionBegin;
2277: #if defined(PETSC_HAVE_HYPRE_DEVICE)
2278:   {
2279:     Mat                 B;
2280:     hypre_ParCSRMatrix *x, *y, *z;

2282:     PetscCall(MatHYPREGetParCSR(Y, &y));
2283:     PetscCall(MatHYPREGetParCSR(X, &x));
2284:     PetscCallExternal(hypre_ParCSRMatrixAdd, 1.0, y, 1.0, x, &z);
2285:     PetscCall(MatCreateFromParCSR(z, MATHYPRE, PETSC_OWN_POINTER, &B));
2286:     PetscCall(MatHeaderMerge(Y, &B));
2287:   }
2288: #else
2289:   if (str == SAME_NONZERO_PATTERN) {
2290:     hypre_ParCSRMatrix *x, *y;
2291:     hypre_CSRMatrix    *xloc, *yloc;
2292:     PetscInt            xnnz, ynnz;
2293:     HYPRE_Complex      *xarr, *yarr;
2294:     PetscBLASInt        one = 1, bnz;

2296:     PetscCall(MatHYPREGetParCSR(Y, &y));
2297:     PetscCall(MatHYPREGetParCSR(X, &x));

2299:     /* diagonal block */
2300:     xloc = hypre_ParCSRMatrixDiag(x);
2301:     yloc = hypre_ParCSRMatrixDiag(y);
2302:     xnnz = 0;
2303:     ynnz = 0;
2304:     xarr = NULL;
2305:     yarr = NULL;
2306:     if (xloc) {
2307:       xarr = hypre_CSRMatrixData(xloc);
2308:       xnnz = hypre_CSRMatrixNumNonzeros(xloc);
2309:     }
2310:     if (yloc) {
2311:       yarr = hypre_CSRMatrixData(yloc);
2312:       ynnz = hypre_CSRMatrixNumNonzeros(yloc);
2313:     }
2314:     PetscCheck(xnnz == ynnz, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Different number of nonzeros in diagonal block %" PetscInt_FMT " != %" PetscInt_FMT, xnnz, ynnz);
2315:     PetscCall(PetscBLASIntCast(xnnz, &bnz));
2316:     PetscCallBLAS("BLASaxpy", BLASaxpy_(&bnz, &a, (PetscScalar *)xarr, &one, (PetscScalar *)yarr, &one));

2318:     /* off-diagonal block */
2319:     xloc = hypre_ParCSRMatrixOffd(x);
2320:     yloc = hypre_ParCSRMatrixOffd(y);
2321:     xnnz = 0;
2322:     ynnz = 0;
2323:     xarr = NULL;
2324:     yarr = NULL;
2325:     if (xloc) {
2326:       xarr = hypre_CSRMatrixData(xloc);
2327:       xnnz = hypre_CSRMatrixNumNonzeros(xloc);
2328:     }
2329:     if (yloc) {
2330:       yarr = hypre_CSRMatrixData(yloc);
2331:       ynnz = hypre_CSRMatrixNumNonzeros(yloc);
2332:     }
2333:     PetscCheck(xnnz == ynnz, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Different number of nonzeros in off-diagonal block %" PetscInt_FMT " != %" PetscInt_FMT, xnnz, ynnz);
2334:     PetscCall(PetscBLASIntCast(xnnz, &bnz));
2335:     PetscCallBLAS("BLASaxpy", BLASaxpy_(&bnz, &a, (PetscScalar *)xarr, &one, (PetscScalar *)yarr, &one));
2336:   } else if (str == SUBSET_NONZERO_PATTERN) {
2337:     PetscCall(MatAXPY_Basic(Y, a, X, str));
2338:   } else {
2339:     Mat B;

2341:     PetscCall(MatAXPY_Basic_Preallocate(Y, X, &B));
2342:     PetscCall(MatAXPY_BasicWithPreallocation(B, Y, a, X, str));
2343:     PetscCall(MatHeaderReplace(Y, &B));
2344:   }
2345: #endif
2346:   PetscFunctionReturn(PETSC_SUCCESS);
2347: }

2349: static PetscErrorCode MatDuplicate_HYPRE(Mat A, MatDuplicateOption op, Mat *B)
2350: {
2351:   hypre_ParCSRMatrix *parcsr = NULL;
2352:   PetscCopyMode       cpmode;
2353:   Mat_HYPRE          *hA;

2355:   PetscFunctionBegin;
2356:   PetscCall(MatHYPREGetParCSR_HYPRE(A, &parcsr));
2357:   if (op == MAT_DO_NOT_COPY_VALUES || op == MAT_SHARE_NONZERO_PATTERN) {
2358:     parcsr = hypre_ParCSRMatrixClone(parcsr, 0);
2359:     cpmode = PETSC_OWN_POINTER;
2360:   } else {
2361:     cpmode = PETSC_COPY_VALUES;
2362:   }
2363:   PetscCall(MatCreateFromParCSR(parcsr, MATHYPRE, cpmode, B));
2364:   hA = (Mat_HYPRE *)A->data;
2365:   if (hA->cooMat) {
2366:     Mat_HYPRE *hB = (Mat_HYPRE *)((*B)->data);
2367:     op            = (op == MAT_DO_NOT_COPY_VALUES) ? op : MAT_COPY_VALUES;
2368:     /* Cannot simply increase the reference count of hA->cooMat, since B needs to share cooMat's data array */
2369:     PetscCall(MatDuplicate(hA->cooMat, op, &hB->cooMat));
2370:     PetscCall(MatHYPRE_AttachCOOMat(*B));
2371:   }
2372:   PetscFunctionReturn(PETSC_SUCCESS);
2373: }

2375: static PetscErrorCode MatSetPreallocationCOO_HYPRE(Mat mat, PetscCount coo_n, PetscInt coo_i[], PetscInt coo_j[])
2376: {
2377:   Mat_HYPRE *hmat = (Mat_HYPRE *)mat->data;

2379:   PetscFunctionBegin;
2380:   /* Build an agent matrix cooMat with AIJ format
2381:      It has the same sparsity pattern as mat, and also shares the data array with mat. We use cooMat to do the COO work.
2382:    */
2383:   PetscCall(MatHYPRE_CreateCOOMat(mat));
2384:   PetscCall(MatSetOption(hmat->cooMat, MAT_IGNORE_OFF_PROC_ENTRIES, hmat->donotstash));
2385:   PetscCall(MatSetOption(hmat->cooMat, MAT_NO_OFF_PROC_ENTRIES, mat->nooffprocentries));

2387:   /* MatSetPreallocationCOO_SeqAIJ and MatSetPreallocationCOO_MPIAIJ uses this specific
2388:      name to automatically put the diagonal entries first */
2389:   PetscCall(PetscObjectSetName((PetscObject)hmat->cooMat, "_internal_COO_mat_for_hypre"));
2390:   PetscCall(MatSetPreallocationCOO(hmat->cooMat, coo_n, coo_i, coo_j));
2391:   hmat->cooMat->assembled = PETSC_TRUE;

2393:   /* Copy the sparsity pattern from cooMat to hypre IJMatrix hmat->ij */
2394:   PetscCall(MatSetOption(mat, MAT_SORTED_FULL, PETSC_TRUE));
2395:   PetscCall(MatHYPRE_CreateFromMat(hmat->cooMat, hmat));      /* Create hmat->ij and preallocate it */
2396:   PetscCall(MatHYPRE_IJMatrixCopyIJ(hmat->cooMat, hmat->ij)); /* Copy A's (i,j) to hmat->ij */

2398:   mat->preallocated = PETSC_TRUE;
2399:   PetscCall(MatAssemblyBegin(mat, MAT_FINAL_ASSEMBLY));
2400:   PetscCall(MatAssemblyEnd(mat, MAT_FINAL_ASSEMBLY)); /* Migrate mat to device if it is bound to. Hypre builds its own SpMV context here */

2402:   /* Attach cooMat to mat */
2403:   PetscCall(MatHYPRE_AttachCOOMat(mat));
2404:   PetscFunctionReturn(PETSC_SUCCESS);
2405: }

2407: static PetscErrorCode MatSetValuesCOO_HYPRE(Mat mat, const PetscScalar v[], InsertMode imode)
2408: {
2409:   Mat_HYPRE *hmat = (Mat_HYPRE *)mat->data;

2411:   PetscFunctionBegin;
2412:   PetscCheck(hmat->cooMat, PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "HYPRE COO delegate matrix has not been created yet");
2413:   PetscCall(MatSetValuesCOO(hmat->cooMat, v, imode));
2414:   PetscCall(MatViewFromOptions(hmat->cooMat, (PetscObject)mat, "-cooMat_view"));
2415:   PetscFunctionReturn(PETSC_SUCCESS);
2416: }

2418: static PetscErrorCode MatGetCurrentMemType_HYPRE(Mat A, PetscMemType *m)
2419: {
2420:   PetscBool petsconcpu;

2422:   PetscFunctionBegin;
2423:   PetscCall(MatBoundToCPU(A, &petsconcpu));
2424:   *m = petsconcpu ? PETSC_MEMTYPE_HOST : PETSC_MEMTYPE_DEVICE;
2425:   PetscFunctionReturn(PETSC_SUCCESS);
2426: }

2428: /*MC
2429:    MATHYPRE - "hypre" - A matrix type to be used for sequential and parallel sparse matrices
2430:           based on the hypre IJ interface.

2432:    Level: intermediate

2434: .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MatHYPRESetPreallocation`
2435: M*/
2436: PETSC_EXTERN PetscErrorCode MatCreate_HYPRE(Mat B)
2437: {
2438:   Mat_HYPRE *hB;
2439: #if defined(PETSC_HAVE_HYPRE_DEVICE)
2440:   HYPRE_MemoryLocation memory_location;
2441: #endif

2443:   PetscFunctionBegin;
2444:   PetscHYPREInitialize();
2445:   PetscCall(PetscNew(&hB));

2447:   hB->inner_free      = PETSC_TRUE;
2448:   hB->array_available = PETSC_TRUE;

2450:   B->data = (void *)hB;

2452:   PetscCall(PetscMemzero(B->ops, sizeof(struct _MatOps)));
2453:   B->ops->mult                  = MatMult_HYPRE;
2454:   B->ops->multtranspose         = MatMultTranspose_HYPRE;
2455:   B->ops->multadd               = MatMultAdd_HYPRE;
2456:   B->ops->multtransposeadd      = MatMultTransposeAdd_HYPRE;
2457:   B->ops->setup                 = MatSetUp_HYPRE;
2458:   B->ops->destroy               = MatDestroy_HYPRE;
2459:   B->ops->assemblyend           = MatAssemblyEnd_HYPRE;
2460:   B->ops->assemblybegin         = MatAssemblyBegin_HYPRE;
2461:   B->ops->setvalues             = MatSetValues_HYPRE;
2462:   B->ops->missingdiagonal       = MatMissingDiagonal_HYPRE;
2463:   B->ops->scale                 = MatScale_HYPRE;
2464:   B->ops->zerorowscolumns       = MatZeroRowsColumns_HYPRE;
2465:   B->ops->zeroentries           = MatZeroEntries_HYPRE;
2466:   B->ops->zerorows              = MatZeroRows_HYPRE;
2467:   B->ops->getrow                = MatGetRow_HYPRE;
2468:   B->ops->restorerow            = MatRestoreRow_HYPRE;
2469:   B->ops->getvalues             = MatGetValues_HYPRE;
2470:   B->ops->setoption             = MatSetOption_HYPRE;
2471:   B->ops->duplicate             = MatDuplicate_HYPRE;
2472:   B->ops->copy                  = MatCopy_HYPRE;
2473:   B->ops->view                  = MatView_HYPRE;
2474:   B->ops->getdiagonal           = MatGetDiagonal_HYPRE;
2475:   B->ops->axpy                  = MatAXPY_HYPRE;
2476:   B->ops->productsetfromoptions = MatProductSetFromOptions_HYPRE;
2477:   B->ops->getcurrentmemtype     = MatGetCurrentMemType_HYPRE;
2478: #if defined(PETSC_HAVE_HYPRE_DEVICE)
2479:   B->ops->bindtocpu = MatBindToCPU_HYPRE;
2480:   /* Get hypre's default memory location. Users can control this using the corresponding HYPRE_SetMemoryLocation API */
2481:   PetscCallExternal(HYPRE_GetMemoryLocation, &memory_location);
2482:   B->boundtocpu = (memory_location == HYPRE_MEMORY_HOST) ? PETSC_TRUE : PETSC_FALSE;
2483: #endif

2485:   /* build cache for off array entries formed */
2486:   PetscCall(MatStashCreate_Private(PetscObjectComm((PetscObject)B), 1, &B->stash));

2488:   PetscCall(PetscCommGetComm(PetscObjectComm((PetscObject)B), &hB->comm));
2489:   PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATHYPRE));
2490:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_hypre_aij_C", MatConvert_HYPRE_AIJ));
2491:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_hypre_is_C", MatConvert_HYPRE_IS));
2492:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqaij_hypre_C", MatProductSetFromOptions_HYPRE));
2493:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_mpiaij_hypre_C", MatProductSetFromOptions_HYPRE));
2494:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatHYPRESetPreallocation_C", MatHYPRESetPreallocation_HYPRE));
2495:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatHYPREGetParCSR_C", MatHYPREGetParCSR_HYPRE));
2496:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSetPreallocationCOO_C", MatSetPreallocationCOO_HYPRE));
2497:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSetValuesCOO_C", MatSetValuesCOO_HYPRE));
2498: #if defined(PETSC_HAVE_HYPRE_DEVICE)
2499:   #if defined(HYPRE_USING_HIP)
2500:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqaijhipsparse_hypre_C", MatProductSetFromOptions_HYPRE));
2501:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_mpiaijhipsparse_hypre_C", MatProductSetFromOptions_HYPRE));
2502:   PetscCall(PetscDeviceInitialize(PETSC_DEVICE_HIP));
2503:   PetscCall(MatSetVecType(B, VECHIP));
2504:   #endif
2505:   #if defined(HYPRE_USING_CUDA)
2506:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqaijcusparse_hypre_C", MatProductSetFromOptions_HYPRE));
2507:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_mpiaijcusparse_hypre_C", MatProductSetFromOptions_HYPRE));
2508:   PetscCall(PetscDeviceInitialize(PETSC_DEVICE_CUDA));
2509:   PetscCall(MatSetVecType(B, VECCUDA));
2510:   #endif
2511: #endif
2512:   PetscFunctionReturn(PETSC_SUCCESS);
2513: }