Actual source code: mhypre.c

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

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

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

 22: static PetscErrorCode MatHYPRE_CreateFromMat(Mat, Mat_HYPRE *);
 23: static PetscErrorCode MatHYPRE_IJMatrixPreallocate(Mat, Mat, HYPRE_IJMatrix);
 24: static PetscErrorCode MatHYPRE_IJMatrixCopyIJ_MPIAIJ(Mat, HYPRE_IJMatrix);
 25: static PetscErrorCode MatHYPRE_IJMatrixCopyIJ_SeqAIJ(Mat, HYPRE_IJMatrix);
 26: static PetscErrorCode MatHYPRE_MultKernel_Private(Mat, HYPRE_Complex, Vec, HYPRE_Complex, Vec, PetscBool);
 27: static PetscErrorCode MatSetValues_HYPRE(Mat, PetscInt, const PetscInt[], PetscInt, const PetscInt[], const PetscScalar[], InsertMode ins);

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

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

 79: static PetscErrorCode MatHYPRE_CreateFromMat(Mat A, Mat_HYPRE *hA)
 80: {
 81:   PetscInt rstart, rend, cstart, cend;

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

127: static PetscErrorCode MatHYPRE_IJMatrixCopyIJ(Mat A, HYPRE_IJMatrix ij)
128: {
129:   PetscBool flg;

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

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

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

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

178:   aux_matrix                               = (hypre_AuxParCSRMatrix *)hypre_IJMatrixTranslator(ij);
179:   hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 0;
180:   PetscFunctionReturn(PETSC_SUCCESS);
181: }

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

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

200:   PetscCallExternal(HYPRE_IJMatrixGetObjectType, ij, &type);
201:   PetscCheck(type == HYPRE_PARCSR, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Only HYPRE_PARCSR is supported");
202:   PetscCallExternal(HYPRE_IJMatrixGetObject, ij, (void **)&par_matrix);
203:   hdiag = hypre_ParCSRMatrixDiag(par_matrix);
204:   hoffd = hypre_ParCSRMatrixOffd(par_matrix);

206:   if (sameint) {
207:     PetscCall(PetscArraycpy(hdiag->i, pdiag->i, pA->A->rmap->n + 1));
208:   } else {
209:     for (i = 0; i < pA->A->rmap->n + 1; i++) hdiag->i[i] = (HYPRE_Int)pdiag->i[i];
210:   }

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

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

233:   aux_matrix                               = (hypre_AuxParCSRMatrix *)hypre_IJMatrixTranslator(ij);
234:   hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 0;
235:   PetscFunctionReturn(PETSC_SUCCESS);
236: }

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

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

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

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

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

368: static PetscErrorCode MatHYPRE_DestroyCOOMat(Mat mat)
369: {
370:   Mat_HYPRE *hA = (Mat_HYPRE *)mat->data;

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

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

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

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

417: #if defined(PETSC_HAVE_HYPRE_DEVICE)
418:   if (!mat->boundtocpu) { /* mat will be on device, so will cooMat */
419:   #if defined(HYPRE_USING_HIP)
420:     matType = MATAIJHIPSPARSE;
421:   #elif defined(HYPRE_USING_CUDA)
422:     matType = MATAIJCUSPARSE;
423:   #else
424:     SETERRQ(comm, PETSC_ERR_SUP, "Do not know the HYPRE device");
425:   #endif
426:   }
427: #endif

429:   /* Do COO preallocation through cooMat */
430:   PetscCall(MatHYPRE_DestroyCOOMat(mat));
431:   PetscCall(MatCreate(comm, &hmat->cooMat));
432:   PetscCall(MatSetType(hmat->cooMat, matType));
433:   PetscCall(MatSetLayouts(hmat->cooMat, rmap, cmap));

435:   /* allocate local matrices if needed */
436:   PetscCall(MatMPIAIJSetPreallocation(hmat->cooMat, 0, NULL, 0, NULL));
437:   PetscFunctionReturn(PETSC_SUCCESS);
438: }

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

460:   PetscFunctionBegin;
461:   PetscCheck(hmat->cooMat, PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "HYPRE COO delegate matrix has not been created yet");
462:   if (hmat->cooMatAttached) PetscFunctionReturn(PETSC_SUCCESS);
463:   PetscCheck(hmat->cooMat->preallocated, PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "HYPRE COO delegate matrix is not preallocated");
464:   PetscCall(PetscObjectSetName((PetscObject)hmat->cooMat, "_internal_COO_mat_for_hypre"));
465:   PetscCall(PetscObjectGetComm((PetscObject)mat, &comm));
466:   PetscCallMPI(MPI_Comm_size(comm, &size));

468:   /* Alias cooMat's data array to IJMatrix's */
469:   PetscCallExternal(HYPRE_IJMatrixGetObject, hmat->ij, (void **)&parCSR);
470:   diag = hypre_ParCSRMatrixDiag(parCSR);
471:   offd = hypre_ParCSRMatrixOffd(parCSR);

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

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

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

496: static PetscErrorCode CSRtoCOO_Private(PetscInt n, const PetscInt ii[], const PetscInt jj[], PetscCount *ncoo, PetscInt **coo_i, PetscInt **coo_j)
497: {
498:   PetscInt *cooi, *cooj;

500:   PetscFunctionBegin;
501:   *ncoo = ii[n];
502:   PetscCall(PetscMalloc2(*ncoo, &cooi, *ncoo, &cooj));
503:   for (PetscInt i = 0; i < n; i++) {
504:     for (PetscInt j = ii[i]; j < ii[i + 1]; j++) cooi[j] = i;
505:   }
506:   PetscCall(PetscArraycpy(cooj, jj, *ncoo));
507:   *coo_i = cooi;
508:   *coo_j = cooj;
509:   PetscFunctionReturn(PETSC_SUCCESS);
510: }

512: static PetscErrorCode CSRtoCOO_HYPRE_Int_Private(PetscInt n, const HYPRE_Int ii[], const HYPRE_Int jj[], PetscCount *ncoo, PetscInt **coo_i, PetscInt **coo_j)
513: {
514:   PetscInt *cooi, *cooj;

516:   PetscFunctionBegin;
517:   *ncoo = ii[n];
518:   PetscCall(PetscMalloc2(*ncoo, &cooi, *ncoo, &cooj));
519:   for (PetscInt i = 0; i < n; i++) {
520:     for (HYPRE_Int j = ii[i]; j < ii[i + 1]; j++) cooi[j] = i;
521:   }
522:   for (PetscCount i = 0; i < *ncoo; i++) cooj[i] = jj[i];
523:   *coo_i = cooi;
524:   *coo_j = cooj;
525:   PetscFunctionReturn(PETSC_SUCCESS);
526: }

528: static PetscErrorCode MatSeqAIJGetCOO_Private(Mat A, PetscCount *ncoo, PetscInt **coo_i, PetscInt **coo_j)
529: {
530:   PetscInt        n;
531:   const PetscInt *ii, *jj;
532:   PetscBool       done;

534:   PetscFunctionBegin;
535:   PetscCall(MatGetRowIJ(A, 0, PETSC_FALSE, PETSC_FALSE, &n, &ii, &jj, &done));
536:   PetscCheck(done, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Failure for MatGetRowIJ");
537:   PetscCall(CSRtoCOO_Private(n, ii, jj, ncoo, coo_i, coo_j));
538:   PetscCall(MatRestoreRowIJ(A, 0, PETSC_FALSE, PETSC_FALSE, &n, &ii, &jj, &done));
539:   PetscCheck(done, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Failure for MatRestoreRowIJ");
540:   PetscFunctionReturn(PETSC_SUCCESS);
541: }

543: static PetscErrorCode hypreCSRMatrixGetCOO_Private(hypre_CSRMatrix *A, PetscCount *ncoo, PetscInt **coo_i, PetscInt **coo_j)
544: {
545:   PetscInt             n = hypre_CSRMatrixNumRows(A);
546:   HYPRE_Int           *ii, *jj;
547:   HYPRE_MemoryLocation mem = HYPRE_MEMORY_HOST;

549:   PetscFunctionBegin;
550: #if defined(PETSC_HAVE_HYPRE_DEVICE)
551:   mem = hypre_CSRMatrixMemoryLocation(A);
552:   if (mem != HYPRE_MEMORY_HOST) {
553:     PetscCount nnz = hypre_CSRMatrixNumNonzeros(A);
554:     PetscCall(PetscMalloc2(n + 1, &ii, nnz, &jj));
555:     hypre_TMemcpy(ii, hypre_CSRMatrixI(A), HYPRE_Int, n + 1, HYPRE_MEMORY_HOST, mem);
556:     hypre_TMemcpy(jj, hypre_CSRMatrixJ(A), HYPRE_Int, nnz, HYPRE_MEMORY_HOST, mem);
557:   } else {
558: #else
559:   {
560: #endif
561:     ii = hypre_CSRMatrixI(A);
562:     jj = hypre_CSRMatrixJ(A);
563:   }
564:   PetscCall(CSRtoCOO_HYPRE_Int_Private(n, ii, jj, ncoo, coo_i, coo_j));
565:   if (mem != HYPRE_MEMORY_HOST) PetscCall(PetscFree2(ii, jj));
566:   PetscFunctionReturn(PETSC_SUCCESS);
567: }

569: static PetscErrorCode MatSetValuesCOOFromCSRMatrix_Private(Mat A, hypre_CSRMatrix *H)
570: {
571:   PetscBool            iscpu = PETSC_TRUE;
572:   PetscScalar         *a;
573:   HYPRE_MemoryLocation mem = HYPRE_MEMORY_HOST;

575:   PetscFunctionBegin;
576: #if defined(PETSC_HAVE_HYPRE_DEVICE)
577:   mem = hypre_CSRMatrixMemoryLocation(H);
578:   PetscCall(PetscObjectTypeCompare((PetscObject)A, MATSEQAIJ, &iscpu));
579: #endif
580:   if (iscpu && mem != HYPRE_MEMORY_HOST) {
581:     PetscCount nnz = hypre_CSRMatrixNumNonzeros(H);
582:     PetscCall(PetscMalloc1(nnz, &a));
583:     hypre_TMemcpy(a, hypre_CSRMatrixData(H), PetscScalar, nnz, HYPRE_MEMORY_HOST, mem);
584:   } else {
585:     a = (PetscScalar *)hypre_CSRMatrixData(H);
586:   }
587:   PetscCall(MatSetValuesCOO(A, a, INSERT_VALUES));
588:   if (iscpu && mem != HYPRE_MEMORY_HOST) PetscCall(PetscFree(a));
589:   PetscFunctionReturn(PETSC_SUCCESS);
590: }

592: PETSC_INTERN PetscErrorCode MatConvert_AIJ_HYPRE(Mat A, MatType type, MatReuse reuse, Mat *B)
593: {
594:   MPI_Comm   comm = PetscObjectComm((PetscObject)A);
595:   Mat        M = NULL, dH = NULL, oH = NULL, dA = NULL, oA = NULL;
596:   PetscBool  ismpiaij, issbaij, isbaij;
597:   Mat_HYPRE *hA;

599:   PetscFunctionBegin;
600:   PetscCall(PetscObjectTypeCompareAny((PetscObject)A, &issbaij, MATSEQSBAIJ, MATMPIBAIJ, ""));
601:   PetscCall(PetscObjectTypeCompareAny((PetscObject)A, &isbaij, MATSEQBAIJ, MATMPIBAIJ, ""));
602:   if (isbaij || issbaij) { /* handle BAIJ and SBAIJ */
603:     PetscBool ismpi;
604:     MatType   newtype;

606:     PetscCall(PetscObjectTypeCompareAny((PetscObject)A, &ismpi, MATMPISBAIJ, MATMPIBAIJ, ""));
607:     newtype = ismpi ? MATMPIAIJ : MATSEQAIJ;
608:     if (reuse == MAT_REUSE_MATRIX) {
609:       PetscCall(MatConvert(*B, newtype, MAT_INPLACE_MATRIX, B));
610:       PetscCall(MatConvert(A, newtype, MAT_REUSE_MATRIX, B));
611:       PetscCall(MatConvert(*B, MATHYPRE, MAT_INPLACE_MATRIX, B));
612:     } else if (reuse == MAT_INITIAL_MATRIX) {
613:       PetscCall(MatConvert(A, newtype, MAT_INITIAL_MATRIX, B));
614:       PetscCall(MatConvert(*B, MATHYPRE, MAT_INPLACE_MATRIX, B));
615:     } else {
616:       PetscCall(MatConvert(A, newtype, MAT_INPLACE_MATRIX, &A));
617:       PetscCall(MatConvert(A, MATHYPRE, MAT_INPLACE_MATRIX, &A));
618:     }
619:     PetscFunctionReturn(PETSC_SUCCESS);
620:   }

622:   dA = A;
623:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIAIJ, &ismpiaij));
624:   if (ismpiaij) PetscCall(MatMPIAIJGetSeqAIJ(A, &dA, &oA, NULL));

626:   if (reuse != MAT_REUSE_MATRIX) {
627:     PetscCount coo_n;
628:     PetscInt  *coo_i, *coo_j;

630:     PetscCall(MatCreate(comm, &M));
631:     PetscCall(MatSetType(M, MATHYPRE));
632:     PetscCall(MatSetSizes(M, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N));
633:     PetscCall(MatSetOption(M, MAT_SORTED_FULL, PETSC_TRUE));
634:     PetscCall(MatSetOption(M, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE));

636:     hA = (Mat_HYPRE *)M->data;
637:     PetscCall(MatHYPRE_CreateFromMat(A, hA));
638:     PetscCall(MatHYPRE_IJMatrixCopyIJ(A, hA->ij));

640:     PetscCall(MatHYPRE_CreateCOOMat(M));

642:     dH = hA->cooMat;
643:     PetscCall(PetscObjectBaseTypeCompare((PetscObject)hA->cooMat, MATMPIAIJ, &ismpiaij));
644:     if (ismpiaij) PetscCall(MatMPIAIJGetSeqAIJ(hA->cooMat, &dH, &oH, NULL));

646:     PetscCall(PetscObjectSetName((PetscObject)dH, "_internal_COO_mat_for_hypre"));
647:     PetscCall(MatSeqAIJGetCOO_Private(dA, &coo_n, &coo_i, &coo_j));
648:     PetscCall(MatSetPreallocationCOO(dH, coo_n, coo_i, coo_j));
649:     PetscCall(PetscFree2(coo_i, coo_j));
650:     if (oH) {
651:       PetscCall(PetscLayoutDestroy(&oH->cmap));
652:       PetscCall(PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)oH), oA->cmap->n, oA->cmap->n, 1, &oH->cmap));
653:       PetscCall(MatSeqAIJGetCOO_Private(oA, &coo_n, &coo_i, &coo_j));
654:       PetscCall(MatSetPreallocationCOO(oH, coo_n, coo_i, coo_j));
655:       PetscCall(PetscFree2(coo_i, coo_j));
656:     }
657:     hA->cooMat->assembled = PETSC_TRUE;

659:     M->preallocated = PETSC_TRUE;
660:     PetscCall(MatAssemblyBegin(M, MAT_FINAL_ASSEMBLY));
661:     PetscCall(MatAssemblyEnd(M, MAT_FINAL_ASSEMBLY));

663:     PetscCall(MatHYPRE_AttachCOOMat(M));
664:     if (reuse == MAT_INITIAL_MATRIX) *B = M;
665:   } else M = *B;

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

670:   dH = hA->cooMat;
671:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)hA->cooMat, MATMPIAIJ, &ismpiaij));
672:   if (ismpiaij) PetscCall(MatMPIAIJGetSeqAIJ(hA->cooMat, &dH, &oH, NULL));

674:   PetscScalar *a;
675:   PetscCall(MatSeqAIJGetCSRAndMemType(dA, NULL, NULL, &a, NULL));
676:   PetscCall(MatSetValuesCOO(dH, a, INSERT_VALUES));
677:   if (oH) {
678:     PetscCall(MatSeqAIJGetCSRAndMemType(oA, NULL, NULL, &a, NULL));
679:     PetscCall(MatSetValuesCOO(oH, a, INSERT_VALUES));
680:   }

682:   if (reuse == MAT_INPLACE_MATRIX) PetscCall(MatHeaderReplace(A, &M));
683:   PetscFunctionReturn(PETSC_SUCCESS);
684: }

686: static PetscErrorCode MatConvert_HYPRE_AIJ(Mat A, MatType mtype, MatReuse reuse, Mat *B)
687: {
688:   Mat                 M, dA = NULL, oA = NULL;
689:   hypre_ParCSRMatrix *parcsr;
690:   hypre_CSRMatrix    *dH, *oH;
691:   MPI_Comm            comm;
692:   PetscBool           ismpiaij, isseqaij;

694:   PetscFunctionBegin;
695:   comm = PetscObjectComm((PetscObject)A);
696:   if (reuse == MAT_REUSE_MATRIX) {
697:     PetscCall(PetscObjectBaseTypeCompare((PetscObject)*B, MATMPIAIJ, &ismpiaij));
698:     PetscCall(PetscObjectBaseTypeCompare((PetscObject)*B, MATSEQAIJ, &isseqaij));
699:     PetscCheck(ismpiaij || isseqaij, comm, PETSC_ERR_SUP, "Only MATMPIAIJ or MATSEQAIJ base types are supported");
700:   }
701:   PetscCall(MatHYPREGetParCSR(A, &parcsr));
702: #if defined(PETSC_HAVE_HYPRE_DEVICE)
703:   if (HYPRE_MEMORY_DEVICE == hypre_ParCSRMatrixMemoryLocation(parcsr)) {
704:     PetscBool isaij;

706:     PetscCall(PetscStrcmp(mtype, MATAIJ, &isaij));
707:     if (isaij) {
708:       PetscMPIInt size;

710:       PetscCallMPI(MPI_Comm_size(comm, &size));
711:   #if defined(HYPRE_USING_HIP)
712:       mtype = size > 1 ? MATMPIAIJHIPSPARSE : MATSEQAIJHIPSPARSE;
713:   #elif defined(HYPRE_USING_CUDA)
714:       mtype = size > 1 ? MATMPIAIJCUSPARSE : MATSEQAIJCUSPARSE;
715:   #else
716:       mtype = size > 1 ? MATMPIAIJ : MATSEQAIJ;
717:   #endif
718:     }
719:   }
720: #endif
721:   dH = hypre_ParCSRMatrixDiag(parcsr);
722:   oH = hypre_ParCSRMatrixOffd(parcsr);
723:   if (reuse != MAT_REUSE_MATRIX) {
724:     PetscCount coo_n;
725:     PetscInt  *coo_i, *coo_j;

727:     PetscCall(MatCreate(comm, &M));
728:     PetscCall(MatSetType(M, mtype));
729:     PetscCall(MatSetSizes(M, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N));
730:     PetscCall(MatMPIAIJSetPreallocation(M, 0, NULL, 0, NULL));

732:     dA = M;
733:     PetscCall(PetscObjectBaseTypeCompare((PetscObject)M, MATMPIAIJ, &ismpiaij));
734:     if (ismpiaij) PetscCall(MatMPIAIJGetSeqAIJ(M, &dA, &oA, NULL));

736:     PetscCall(hypreCSRMatrixGetCOO_Private(dH, &coo_n, &coo_i, &coo_j));
737:     PetscCall(MatSetPreallocationCOO(dA, coo_n, coo_i, coo_j));
738:     PetscCall(PetscFree2(coo_i, coo_j));
739:     if (ismpiaij) {
740:       HYPRE_Int nc = hypre_CSRMatrixNumCols(oH);

742:       PetscCall(PetscLayoutDestroy(&oA->cmap));
743:       PetscCall(PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)oA), nc, nc, 1, &oA->cmap));
744:       PetscCall(hypreCSRMatrixGetCOO_Private(oH, &coo_n, &coo_i, &coo_j));
745:       PetscCall(MatSetPreallocationCOO(oA, coo_n, coo_i, coo_j));
746:       PetscCall(PetscFree2(coo_i, coo_j));

748:       /* garray */
749:       Mat_MPIAIJ   *aij    = (Mat_MPIAIJ *)M->data;
750:       HYPRE_BigInt *harray = hypre_ParCSRMatrixColMapOffd(parcsr);
751:       PetscInt     *garray;

753:       PetscCall(PetscFree(aij->garray));
754:       PetscCall(PetscMalloc1(nc, &garray));
755:       for (HYPRE_Int i = 0; i < nc; i++) garray[i] = (PetscInt)harray[i];
756:       aij->garray = garray;
757:       PetscCall(MatSetUpMultiply_MPIAIJ(M));
758:     }
759:     if (reuse == MAT_INITIAL_MATRIX) *B = M;
760:   } else M = *B;

762:   dA = M;
763:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)M, MATMPIAIJ, &ismpiaij));
764:   if (ismpiaij) PetscCall(MatMPIAIJGetSeqAIJ(M, &dA, &oA, NULL));
765:   PetscCall(MatSetValuesCOOFromCSRMatrix_Private(dA, dH));
766:   if (oA) PetscCall(MatSetValuesCOOFromCSRMatrix_Private(oA, oH));
767:   M->assembled = PETSC_TRUE;
768:   if (reuse == MAT_INPLACE_MATRIX) PetscCall(MatHeaderReplace(A, &M));
769:   PetscFunctionReturn(PETSC_SUCCESS);
770: }

772: static PetscErrorCode MatAIJGetParCSR_Private(Mat A, hypre_ParCSRMatrix **hA)
773: {
774:   hypre_ParCSRMatrix *tA;
775:   hypre_CSRMatrix    *hdiag, *hoffd;
776:   Mat_SeqAIJ         *diag, *offd;
777:   PetscInt           *garray, i, noffd, dnnz, onnz, *row_starts, *col_starts;
778:   MPI_Comm            comm = PetscObjectComm((PetscObject)A);
779:   PetscBool           ismpiaij, isseqaij;
780:   PetscBool           sameint = (PetscBool)(sizeof(PetscInt) == sizeof(HYPRE_Int));
781:   HYPRE_Int          *hdi = NULL, *hdj = NULL, *hoi = NULL, *hoj = NULL;
782:   PetscInt           *pdi = NULL, *pdj = NULL, *poi = NULL, *poj = NULL;
783:   PetscBool           iscuda, iship;
784: #if defined(PETSC_HAVE_DEVICE) && defined(PETSC_HAVE_HYPRE_DEVICE)
785:   PetscBool boundtocpu = A->boundtocpu;
786: #else
787:   PetscBool boundtocpu = PETSC_TRUE;
788: #endif

790:   PetscFunctionBegin;
791:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIAIJ, &ismpiaij));
792:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATSEQAIJ, &isseqaij));
793:   PetscCheck(ismpiaij || isseqaij, comm, PETSC_ERR_SUP, "Unsupported type %s", ((PetscObject)A)->type_name);
794:   PetscCall(PetscObjectTypeCompareAny((PetscObject)A, &iscuda, MATSEQAIJHIPSPARSE, MATMPIAIJCUSPARSE, ""));
795:   PetscCall(PetscObjectTypeCompareAny((PetscObject)A, &iship, MATSEQAIJCUSPARSE, MATMPIAIJHIPSPARSE, ""));
796:   PetscHYPREInitialize();
797:   if (ismpiaij) {
798:     Mat_MPIAIJ *a = (Mat_MPIAIJ *)A->data;

800:     diag = (Mat_SeqAIJ *)a->A->data;
801:     offd = (Mat_SeqAIJ *)a->B->data;
802:     if (!boundtocpu && (iscuda || iship)) {
803: #if defined(HYPRE_USING_CUDA) && defined(PETSC_HAVE_CUDA)
804:       if (iscuda) {
805:         sameint = PETSC_TRUE;
806:         PetscCall(MatSeqAIJCUSPARSEGetIJ(a->A, PETSC_FALSE, (const HYPRE_Int **)&hdi, (const HYPRE_Int **)&hdj));
807:         PetscCall(MatSeqAIJCUSPARSEGetIJ(a->B, PETSC_FALSE, (const HYPRE_Int **)&hoi, (const HYPRE_Int **)&hoj));
808:       }
809: #endif
810: #if defined(HYPRE_USING_HIP) && defined(PETSC_HAVE_HIP)
811:       if (iship) {
812:         sameint = PETSC_TRUE;
813:         PetscCall(MatSeqAIJHIPSPARSEGetIJ(a->A, PETSC_FALSE, (const HYPRE_Int **)&hdi, (const HYPRE_Int **)&hdj));
814:         PetscCall(MatSeqAIJHIPSPARSEGetIJ(a->B, PETSC_FALSE, (const HYPRE_Int **)&hoi, (const HYPRE_Int **)&hoj));
815:       }
816: #endif
817:     } else {
818:       boundtocpu = PETSC_TRUE;
819:       pdi        = diag->i;
820:       pdj        = diag->j;
821:       poi        = offd->i;
822:       poj        = offd->j;
823:       if (sameint) {
824:         hdi = (HYPRE_Int *)pdi;
825:         hdj = (HYPRE_Int *)pdj;
826:         hoi = (HYPRE_Int *)poi;
827:         hoj = (HYPRE_Int *)poj;
828:       }
829:     }
830:     garray = a->garray;
831:     noffd  = a->B->cmap->N;
832:     dnnz   = diag->nz;
833:     onnz   = offd->nz;
834:   } else {
835:     diag = (Mat_SeqAIJ *)A->data;
836:     offd = NULL;
837:     if (!boundtocpu && (iscuda || iship)) {
838: #if defined(HYPRE_USING_CUDA) && defined(PETSC_HAVE_CUDA)
839:       if (iscuda) {
840:         sameint = PETSC_TRUE;
841:         PetscCall(MatSeqAIJCUSPARSEGetIJ(A, PETSC_FALSE, (const HYPRE_Int **)&hdi, (const HYPRE_Int **)&hdj));
842:       }
843: #endif
844: #if defined(HYPRE_USING_HIP) && defined(PETSC_HAVE_HIP)
845:       if (iship) {
846:         sameint = PETSC_TRUE;
847:         PetscCall(MatSeqAIJHIPSPARSEGetIJ(A, PETSC_FALSE, (const HYPRE_Int **)&hdi, (const HYPRE_Int **)&hdj));
848:       }
849: #endif
850:     } else {
851:       boundtocpu = PETSC_TRUE;
852:       pdi        = diag->i;
853:       pdj        = diag->j;
854:       if (sameint) {
855:         hdi = (HYPRE_Int *)pdi;
856:         hdj = (HYPRE_Int *)pdj;
857:       }
858:     }
859:     garray = NULL;
860:     noffd  = 0;
861:     dnnz   = diag->nz;
862:     onnz   = 0;
863:   }

865:   /* create a temporary ParCSR */
866:   if (HYPRE_AssumedPartitionCheck()) {
867:     PetscMPIInt myid;

869:     PetscCallMPI(MPI_Comm_rank(comm, &myid));
870:     row_starts = A->rmap->range + myid;
871:     col_starts = A->cmap->range + myid;
872:   } else {
873:     row_starts = A->rmap->range;
874:     col_starts = A->cmap->range;
875:   }
876:   tA = hypre_ParCSRMatrixCreate(comm, A->rmap->N, A->cmap->N, (HYPRE_BigInt *)row_starts, (HYPRE_BigInt *)col_starts, noffd, dnnz, onnz);
877: #if defined(hypre_ParCSRMatrixOwnsRowStarts)
878:   hypre_ParCSRMatrixSetRowStartsOwner(tA, 0);
879:   hypre_ParCSRMatrixSetColStartsOwner(tA, 0);
880: #endif

882:   /* set diagonal part */
883:   hdiag = hypre_ParCSRMatrixDiag(tA);
884:   if (!sameint) { /* malloc CSR pointers */
885:     PetscCall(PetscMalloc2(A->rmap->n + 1, &hdi, dnnz, &hdj));
886:     for (i = 0; i < A->rmap->n + 1; i++) hdi[i] = (HYPRE_Int)pdi[i];
887:     for (i = 0; i < dnnz; i++) hdj[i] = (HYPRE_Int)pdj[i];
888:   }
889:   hypre_CSRMatrixI(hdiag)           = hdi;
890:   hypre_CSRMatrixJ(hdiag)           = hdj;
891:   hypre_CSRMatrixData(hdiag)        = (HYPRE_Complex *)diag->a;
892:   hypre_CSRMatrixNumNonzeros(hdiag) = diag->nz;
893:   hypre_CSRMatrixSetRownnz(hdiag);
894:   hypre_CSRMatrixSetDataOwner(hdiag, 0);

896:   /* set off-diagonal part */
897:   hoffd = hypre_ParCSRMatrixOffd(tA);
898:   if (offd) {
899:     if (!sameint) { /* malloc CSR pointers */
900:       PetscCall(PetscMalloc2(A->rmap->n + 1, &hoi, onnz, &hoj));
901:       for (i = 0; i < A->rmap->n + 1; i++) hoi[i] = (HYPRE_Int)poi[i];
902:       for (i = 0; i < onnz; i++) hoj[i] = (HYPRE_Int)poj[i];
903:     }
904:     hypre_CSRMatrixI(hoffd)           = hoi;
905:     hypre_CSRMatrixJ(hoffd)           = hoj;
906:     hypre_CSRMatrixData(hoffd)        = (HYPRE_Complex *)offd->a;
907:     hypre_CSRMatrixNumNonzeros(hoffd) = offd->nz;
908:     hypre_CSRMatrixSetRownnz(hoffd);
909:     hypre_CSRMatrixSetDataOwner(hoffd, 0);
910:   }
911: #if defined(PETSC_HAVE_HYPRE_DEVICE)
912:   PetscCallExternal(hypre_ParCSRMatrixInitialize_v2, tA, !boundtocpu ? HYPRE_MEMORY_DEVICE : HYPRE_MEMORY_HOST);
913: #else
914:   #if PETSC_PKG_HYPRE_VERSION_LT(2, 18, 0)
915:   PetscCallExternal(hypre_ParCSRMatrixInitialize, tA);
916:   #else
917:   PetscCallExternal(hypre_ParCSRMatrixInitialize_v2, tA, HYPRE_MEMORY_HOST);
918:   #endif
919: #endif
920:   hypre_TFree(hypre_ParCSRMatrixColMapOffd(tA), HYPRE_MEMORY_HOST);
921:   hypre_ParCSRMatrixSetNumNonzeros(tA);
922:   hypre_ParCSRMatrixColMapOffd(tA) = (HYPRE_BigInt *)garray;
923:   if (!hypre_ParCSRMatrixCommPkg(tA)) PetscCallExternal(hypre_MatvecCommPkgCreate, tA);
924:   *hA = tA;
925:   PetscFunctionReturn(PETSC_SUCCESS);
926: }

928: static PetscErrorCode MatAIJRestoreParCSR_Private(Mat A, hypre_ParCSRMatrix **hA)
929: {
930:   hypre_CSRMatrix *hdiag, *hoffd;
931:   PetscBool        ismpiaij, sameint = (PetscBool)(sizeof(PetscInt) == sizeof(HYPRE_Int));
932: #if defined(PETSC_HAVE_HYPRE_DEVICE)
933:   PetscBool iscuda = PETSC_FALSE;
934: #endif

936:   PetscFunctionBegin;
937:   PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPIAIJ, &ismpiaij));
938: #if defined(PETSC_HAVE_HYPRE_DEVICE)
939:   PetscCall(PetscObjectTypeCompareAny((PetscObject)A, &iscuda, MATSEQAIJCUSPARSE, MATMPIAIJCUSPARSE, ""));
940:   if (iscuda) sameint = PETSC_TRUE;
941: #endif
942:   hdiag = hypre_ParCSRMatrixDiag(*hA);
943:   hoffd = hypre_ParCSRMatrixOffd(*hA);
944:   /* free temporary memory allocated by PETSc
945:      set pointers to NULL before destroying tA */
946:   if (!sameint) {
947:     HYPRE_Int *hi, *hj;

949:     hi = hypre_CSRMatrixI(hdiag);
950:     hj = hypre_CSRMatrixJ(hdiag);
951:     PetscCall(PetscFree2(hi, hj));
952:     if (ismpiaij) {
953:       hi = hypre_CSRMatrixI(hoffd);
954:       hj = hypre_CSRMatrixJ(hoffd);
955:       PetscCall(PetscFree2(hi, hj));
956:     }
957:   }
958:   hypre_CSRMatrixI(hdiag)    = NULL;
959:   hypre_CSRMatrixJ(hdiag)    = NULL;
960:   hypre_CSRMatrixData(hdiag) = NULL;
961:   if (ismpiaij) {
962:     hypre_CSRMatrixI(hoffd)    = NULL;
963:     hypre_CSRMatrixJ(hoffd)    = NULL;
964:     hypre_CSRMatrixData(hoffd) = NULL;
965:   }
966:   hypre_ParCSRMatrixColMapOffd(*hA) = NULL;
967:   hypre_ParCSRMatrixDestroy(*hA);
968:   *hA = NULL;
969:   PetscFunctionReturn(PETSC_SUCCESS);
970: }

972: /* calls RAP from BoomerAMG:
973:    the resulting ParCSR will not own the column and row starts
974:    It looks like we don't need to have the diagonal entries ordered first */
975: static PetscErrorCode MatHYPRE_ParCSR_RAP(hypre_ParCSRMatrix *hR, hypre_ParCSRMatrix *hA, hypre_ParCSRMatrix *hP, hypre_ParCSRMatrix **hRAP)
976: {
977: #if defined(hypre_ParCSRMatrixOwnsRowStarts)
978:   HYPRE_Int P_owns_col_starts, R_owns_row_starts;
979: #endif

981:   PetscFunctionBegin;
982: #if defined(hypre_ParCSRMatrixOwnsRowStarts)
983:   P_owns_col_starts = hypre_ParCSRMatrixOwnsColStarts(hP);
984:   R_owns_row_starts = hypre_ParCSRMatrixOwnsRowStarts(hR);
985: #endif
986:   /* can be replaced by version test later */
987: #if defined(PETSC_HAVE_HYPRE_DEVICE)
988:   PetscStackPushExternal("hypre_ParCSRMatrixRAP");
989:   *hRAP = hypre_ParCSRMatrixRAP(hR, hA, hP);
990:   PetscStackPop;
991: #else
992:   PetscCallExternal(hypre_BoomerAMGBuildCoarseOperator, hR, hA, hP, hRAP);
993:   PetscCallExternal(hypre_ParCSRMatrixSetNumNonzeros, *hRAP);
994: #endif
995:   /* hypre_BoomerAMGBuildCoarseOperator steals the col_starts from P and the row_starts from R */
996: #if defined(hypre_ParCSRMatrixOwnsRowStarts)
997:   hypre_ParCSRMatrixSetRowStartsOwner(*hRAP, 0);
998:   hypre_ParCSRMatrixSetColStartsOwner(*hRAP, 0);
999:   if (P_owns_col_starts) hypre_ParCSRMatrixSetColStartsOwner(hP, 1);
1000:   if (R_owns_row_starts) hypre_ParCSRMatrixSetRowStartsOwner(hR, 1);
1001: #endif
1002:   PetscFunctionReturn(PETSC_SUCCESS);
1003: }

1005: static PetscErrorCode MatPtAPNumeric_AIJ_AIJ_wHYPRE(Mat A, Mat P, Mat C)
1006: {
1007:   Mat                 B;
1008:   hypre_ParCSRMatrix *hA, *hP, *hPtAP = NULL;
1009:   Mat_Product        *product = C->product;

1011:   PetscFunctionBegin;
1012:   PetscCall(MatAIJGetParCSR_Private(A, &hA));
1013:   PetscCall(MatAIJGetParCSR_Private(P, &hP));
1014:   PetscCall(MatHYPRE_ParCSR_RAP(hP, hA, hP, &hPtAP));
1015:   PetscCall(MatCreateFromParCSR(hPtAP, MATAIJ, PETSC_OWN_POINTER, &B));

1017:   PetscCall(MatHeaderMerge(C, &B));
1018:   C->product = product;

1020:   PetscCall(MatAIJRestoreParCSR_Private(A, &hA));
1021:   PetscCall(MatAIJRestoreParCSR_Private(P, &hP));
1022:   PetscFunctionReturn(PETSC_SUCCESS);
1023: }

1025: PETSC_INTERN PetscErrorCode MatPtAPSymbolic_AIJ_AIJ_wHYPRE(Mat A, Mat P, PetscReal fill, Mat C)
1026: {
1027:   PetscFunctionBegin;
1028:   PetscCall(MatSetType(C, MATAIJ));
1029:   C->ops->ptapnumeric    = MatPtAPNumeric_AIJ_AIJ_wHYPRE;
1030:   C->ops->productnumeric = MatProductNumeric_PtAP;
1031:   PetscFunctionReturn(PETSC_SUCCESS);
1032: }

1034: static PetscErrorCode MatPtAPNumeric_AIJ_HYPRE(Mat A, Mat P, Mat C)
1035: {
1036:   Mat                 B;
1037:   Mat_HYPRE          *hP;
1038:   hypre_ParCSRMatrix *hA = NULL, *Pparcsr, *ptapparcsr = NULL;
1039:   HYPRE_Int           type;
1040:   MPI_Comm            comm = PetscObjectComm((PetscObject)A);
1041:   PetscBool           ishypre;

1043:   PetscFunctionBegin;
1044:   PetscCall(PetscObjectTypeCompare((PetscObject)P, MATHYPRE, &ishypre));
1045:   PetscCheck(ishypre, comm, PETSC_ERR_USER, "P should be of type %s", MATHYPRE);
1046:   hP = (Mat_HYPRE *)P->data;
1047:   PetscCallExternal(HYPRE_IJMatrixGetObjectType, hP->ij, &type);
1048:   PetscCheck(type == HYPRE_PARCSR, comm, PETSC_ERR_SUP, "Only HYPRE_PARCSR is supported");
1049:   PetscCallExternal(HYPRE_IJMatrixGetObject, hP->ij, (void **)&Pparcsr);

1051:   PetscCall(MatAIJGetParCSR_Private(A, &hA));
1052:   PetscCall(MatHYPRE_ParCSR_RAP(Pparcsr, hA, Pparcsr, &ptapparcsr));
1053:   PetscCall(MatAIJRestoreParCSR_Private(A, &hA));

1055:   /* create temporary matrix and merge to C */
1056:   PetscCall(MatCreateFromParCSR(ptapparcsr, ((PetscObject)C)->type_name, PETSC_OWN_POINTER, &B));
1057:   PetscCall(MatHeaderMerge(C, &B));
1058:   PetscFunctionReturn(PETSC_SUCCESS);
1059: }

1061: static PetscErrorCode MatPtAPNumeric_HYPRE_HYPRE(Mat A, Mat P, Mat C)
1062: {
1063:   Mat                 B;
1064:   hypre_ParCSRMatrix *Aparcsr, *Pparcsr, *ptapparcsr = NULL;
1065:   Mat_HYPRE          *hA, *hP;
1066:   PetscBool           ishypre;
1067:   HYPRE_Int           type;

1069:   PetscFunctionBegin;
1070:   PetscCall(PetscObjectTypeCompare((PetscObject)P, MATHYPRE, &ishypre));
1071:   PetscCheck(ishypre, PetscObjectComm((PetscObject)P), PETSC_ERR_USER, "P should be of type %s", MATHYPRE);
1072:   PetscCall(PetscObjectTypeCompare((PetscObject)A, MATHYPRE, &ishypre));
1073:   PetscCheck(ishypre, PetscObjectComm((PetscObject)A), PETSC_ERR_USER, "A should be of type %s", MATHYPRE);
1074:   hA = (Mat_HYPRE *)A->data;
1075:   hP = (Mat_HYPRE *)P->data;
1076:   PetscCallExternal(HYPRE_IJMatrixGetObjectType, hA->ij, &type);
1077:   PetscCheck(type == HYPRE_PARCSR, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Only HYPRE_PARCSR is supported");
1078:   PetscCallExternal(HYPRE_IJMatrixGetObjectType, hP->ij, &type);
1079:   PetscCheck(type == HYPRE_PARCSR, PetscObjectComm((PetscObject)P), PETSC_ERR_SUP, "Only HYPRE_PARCSR is supported");
1080:   PetscCallExternal(HYPRE_IJMatrixGetObject, hA->ij, (void **)&Aparcsr);
1081:   PetscCallExternal(HYPRE_IJMatrixGetObject, hP->ij, (void **)&Pparcsr);
1082:   PetscCall(MatHYPRE_ParCSR_RAP(Pparcsr, Aparcsr, Pparcsr, &ptapparcsr));
1083:   PetscCall(MatCreateFromParCSR(ptapparcsr, MATHYPRE, PETSC_OWN_POINTER, &B));
1084:   PetscCall(MatHeaderMerge(C, &B));
1085:   PetscFunctionReturn(PETSC_SUCCESS);
1086: }

1088: /* calls hypre_ParMatmul
1089:    hypre_ParMatMul uses hypre_ParMatrixCreate with the communicator of hA
1090:    hypre_ParMatrixCreate does not duplicate the communicator
1091:    It looks like we don't need to have the diagonal entries ordered first */
1092: static PetscErrorCode MatHYPRE_ParCSR_MatMatMult(hypre_ParCSRMatrix *hA, hypre_ParCSRMatrix *hB, hypre_ParCSRMatrix **hAB)
1093: {
1094:   PetscFunctionBegin;
1095:   /* can be replaced by version test later */
1096: #if defined(PETSC_HAVE_HYPRE_DEVICE)
1097:   PetscStackPushExternal("hypre_ParCSRMatMat");
1098:   *hAB = hypre_ParCSRMatMat(hA, hB);
1099: #else
1100:   PetscStackPushExternal("hypre_ParMatmul");
1101:   *hAB = hypre_ParMatmul(hA, hB);
1102: #endif
1103:   PetscStackPop;
1104:   PetscFunctionReturn(PETSC_SUCCESS);
1105: }

1107: static PetscErrorCode MatMatMultNumeric_AIJ_AIJ_wHYPRE(Mat A, Mat B, Mat C)
1108: {
1109:   Mat                 D;
1110:   hypre_ParCSRMatrix *hA, *hB, *hAB = NULL;
1111:   Mat_Product        *product = C->product;

1113:   PetscFunctionBegin;
1114:   PetscCall(MatAIJGetParCSR_Private(A, &hA));
1115:   PetscCall(MatAIJGetParCSR_Private(B, &hB));
1116:   PetscCall(MatHYPRE_ParCSR_MatMatMult(hA, hB, &hAB));
1117:   PetscCall(MatCreateFromParCSR(hAB, MATAIJ, PETSC_OWN_POINTER, &D));

1119:   PetscCall(MatHeaderMerge(C, &D));
1120:   C->product = product;

1122:   PetscCall(MatAIJRestoreParCSR_Private(A, &hA));
1123:   PetscCall(MatAIJRestoreParCSR_Private(B, &hB));
1124:   PetscFunctionReturn(PETSC_SUCCESS);
1125: }

1127: PETSC_INTERN PetscErrorCode MatMatMultSymbolic_AIJ_AIJ_wHYPRE(Mat A, Mat B, PetscReal fill, Mat C)
1128: {
1129:   PetscFunctionBegin;
1130:   PetscCall(MatSetType(C, MATAIJ));
1131:   C->ops->matmultnumeric = MatMatMultNumeric_AIJ_AIJ_wHYPRE;
1132:   C->ops->productnumeric = MatProductNumeric_AB;
1133:   PetscFunctionReturn(PETSC_SUCCESS);
1134: }

1136: static PetscErrorCode MatMatMultNumeric_HYPRE_HYPRE(Mat A, Mat B, Mat C)
1137: {
1138:   Mat                 D;
1139:   hypre_ParCSRMatrix *Aparcsr, *Bparcsr, *ABparcsr = NULL;
1140:   Mat_HYPRE          *hA, *hB;
1141:   PetscBool           ishypre;
1142:   HYPRE_Int           type;
1143:   Mat_Product        *product;

1145:   PetscFunctionBegin;
1146:   PetscCall(PetscObjectTypeCompare((PetscObject)B, MATHYPRE, &ishypre));
1147:   PetscCheck(ishypre, PetscObjectComm((PetscObject)B), PETSC_ERR_USER, "B should be of type %s", MATHYPRE);
1148:   PetscCall(PetscObjectTypeCompare((PetscObject)A, MATHYPRE, &ishypre));
1149:   PetscCheck(ishypre, PetscObjectComm((PetscObject)A), PETSC_ERR_USER, "A should be of type %s", MATHYPRE);
1150:   hA = (Mat_HYPRE *)A->data;
1151:   hB = (Mat_HYPRE *)B->data;
1152:   PetscCallExternal(HYPRE_IJMatrixGetObjectType, hA->ij, &type);
1153:   PetscCheck(type == HYPRE_PARCSR, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Only HYPRE_PARCSR is supported");
1154:   PetscCallExternal(HYPRE_IJMatrixGetObjectType, hB->ij, &type);
1155:   PetscCheck(type == HYPRE_PARCSR, PetscObjectComm((PetscObject)B), PETSC_ERR_SUP, "Only HYPRE_PARCSR is supported");
1156:   PetscCallExternal(HYPRE_IJMatrixGetObject, hA->ij, (void **)&Aparcsr);
1157:   PetscCallExternal(HYPRE_IJMatrixGetObject, hB->ij, (void **)&Bparcsr);
1158:   PetscCall(MatHYPRE_ParCSR_MatMatMult(Aparcsr, Bparcsr, &ABparcsr));
1159:   PetscCall(MatCreateFromParCSR(ABparcsr, MATHYPRE, PETSC_OWN_POINTER, &D));

1161:   /* need to use HeaderReplace because HeaderMerge messes up with the communicator */
1162:   product    = C->product; /* save it from MatHeaderReplace() */
1163:   C->product = NULL;
1164:   PetscCall(MatHeaderReplace(C, &D));
1165:   C->product             = product;
1166:   C->ops->matmultnumeric = MatMatMultNumeric_HYPRE_HYPRE;
1167:   C->ops->productnumeric = MatProductNumeric_AB;
1168:   PetscFunctionReturn(PETSC_SUCCESS);
1169: }

1171: PETSC_INTERN PetscErrorCode MatTransposeMatMatMultNumeric_AIJ_AIJ_AIJ_wHYPRE(Mat A, Mat B, Mat C, Mat D)
1172: {
1173:   Mat                 E;
1174:   hypre_ParCSRMatrix *hA, *hB, *hC, *hABC = NULL;

1176:   PetscFunctionBegin;
1177:   PetscCall(MatAIJGetParCSR_Private(A, &hA));
1178:   PetscCall(MatAIJGetParCSR_Private(B, &hB));
1179:   PetscCall(MatAIJGetParCSR_Private(C, &hC));
1180:   PetscCall(MatHYPRE_ParCSR_RAP(hA, hB, hC, &hABC));
1181:   PetscCall(MatCreateFromParCSR(hABC, MATAIJ, PETSC_OWN_POINTER, &E));
1182:   PetscCall(MatHeaderMerge(D, &E));
1183:   PetscCall(MatAIJRestoreParCSR_Private(A, &hA));
1184:   PetscCall(MatAIJRestoreParCSR_Private(B, &hB));
1185:   PetscCall(MatAIJRestoreParCSR_Private(C, &hC));
1186:   PetscFunctionReturn(PETSC_SUCCESS);
1187: }

1189: PETSC_INTERN PetscErrorCode MatTransposeMatMatMultSymbolic_AIJ_AIJ_AIJ_wHYPRE(Mat A, Mat B, Mat C, PetscReal fill, Mat D)
1190: {
1191:   PetscFunctionBegin;
1192:   PetscCall(MatSetType(D, MATAIJ));
1193:   PetscFunctionReturn(PETSC_SUCCESS);
1194: }

1196: static PetscErrorCode MatProductSymbolic_AB_HYPRE(Mat C)
1197: {
1198:   PetscFunctionBegin;
1199:   C->ops->productnumeric = MatProductNumeric_AB;
1200:   PetscFunctionReturn(PETSC_SUCCESS);
1201: }

1203: static PetscErrorCode MatProductSetFromOptions_HYPRE_AB(Mat C)
1204: {
1205:   Mat_Product *product = C->product;
1206:   PetscBool    Ahypre;

1208:   PetscFunctionBegin;
1209:   PetscCall(PetscObjectTypeCompare((PetscObject)product->A, MATHYPRE, &Ahypre));
1210:   if (Ahypre) { /* A is a Hypre matrix */
1211:     PetscCall(MatSetType(C, MATHYPRE));
1212:     C->ops->productsymbolic = MatProductSymbolic_AB_HYPRE;
1213:     C->ops->matmultnumeric  = MatMatMultNumeric_HYPRE_HYPRE;
1214:     PetscFunctionReturn(PETSC_SUCCESS);
1215:   }
1216:   PetscFunctionReturn(PETSC_SUCCESS);
1217: }

1219: static PetscErrorCode MatProductSymbolic_PtAP_HYPRE(Mat C)
1220: {
1221:   PetscFunctionBegin;
1222:   C->ops->productnumeric = MatProductNumeric_PtAP;
1223:   PetscFunctionReturn(PETSC_SUCCESS);
1224: }

1226: static PetscErrorCode MatProductSetFromOptions_HYPRE_PtAP(Mat C)
1227: {
1228:   Mat_Product *product = C->product;
1229:   PetscBool    flg;
1230:   PetscInt     type        = 0;
1231:   const char  *outTypes[4] = {"aij", "seqaij", "mpiaij", "hypre"};
1232:   PetscInt     ntype       = 4;
1233:   Mat          A           = product->A;
1234:   PetscBool    Ahypre;

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

1245:   /* A is AIJ, P is Hypre, C = PtAP can be either AIJ or Hypre format */
1246:   /* Get runtime option */
1247:   if (product->api_user) {
1248:     PetscOptionsBegin(PetscObjectComm((PetscObject)C), ((PetscObject)C)->prefix, "MatPtAP_HYPRE", "Mat");
1249:     PetscCall(PetscOptionsEList("-matptap_hypre_outtype", "MatPtAP outtype", "MatPtAP outtype", outTypes, ntype, outTypes[type], &type, &flg));
1250:     PetscOptionsEnd();
1251:   } else {
1252:     PetscOptionsBegin(PetscObjectComm((PetscObject)C), ((PetscObject)C)->prefix, "MatProduct_PtAP_HYPRE", "Mat");
1253:     PetscCall(PetscOptionsEList("-mat_product_algorithm_hypre_outtype", "MatProduct_PtAP outtype", "MatProduct_PtAP", outTypes, ntype, outTypes[type], &type, &flg));
1254:     PetscOptionsEnd();
1255:   }

1257:   if (type == 0 || type == 1 || type == 2) {
1258:     PetscCall(MatSetType(C, MATAIJ));
1259:   } else if (type == 3) {
1260:     PetscCall(MatSetType(C, MATHYPRE));
1261:   } else SETERRQ(PetscObjectComm((PetscObject)C), PETSC_ERR_SUP, "MatPtAP outtype is not supported");
1262:   C->ops->productsymbolic = MatProductSymbolic_PtAP_HYPRE;
1263:   C->ops->ptapnumeric     = MatPtAPNumeric_AIJ_HYPRE;
1264:   PetscFunctionReturn(PETSC_SUCCESS);
1265: }

1267: static PetscErrorCode MatProductSetFromOptions_HYPRE(Mat C)
1268: {
1269:   Mat_Product *product = C->product;

1271:   PetscFunctionBegin;
1272:   switch (product->type) {
1273:   case MATPRODUCT_AB:
1274:     PetscCall(MatProductSetFromOptions_HYPRE_AB(C));
1275:     break;
1276:   case MATPRODUCT_PtAP:
1277:     PetscCall(MatProductSetFromOptions_HYPRE_PtAP(C));
1278:     break;
1279:   default:
1280:     break;
1281:   }
1282:   PetscFunctionReturn(PETSC_SUCCESS);
1283: }

1285: static PetscErrorCode MatMultTranspose_HYPRE(Mat A, Vec x, Vec y)
1286: {
1287:   PetscFunctionBegin;
1288:   PetscCall(MatHYPRE_MultKernel_Private(A, 1.0, x, 0.0, y, PETSC_TRUE));
1289:   PetscFunctionReturn(PETSC_SUCCESS);
1290: }

1292: static PetscErrorCode MatMult_HYPRE(Mat A, Vec x, Vec y)
1293: {
1294:   PetscFunctionBegin;
1295:   PetscCall(MatHYPRE_MultKernel_Private(A, 1.0, x, 0.0, y, PETSC_FALSE));
1296:   PetscFunctionReturn(PETSC_SUCCESS);
1297: }

1299: static PetscErrorCode MatMultAdd_HYPRE(Mat A, Vec x, Vec y, Vec z)
1300: {
1301:   PetscFunctionBegin;
1302:   if (y != z) PetscCall(VecCopy(y, z));
1303:   PetscCall(MatHYPRE_MultKernel_Private(A, 1.0, x, 1.0, z, PETSC_FALSE));
1304:   PetscFunctionReturn(PETSC_SUCCESS);
1305: }

1307: static PetscErrorCode MatMultTransposeAdd_HYPRE(Mat A, Vec x, Vec y, Vec z)
1308: {
1309:   PetscFunctionBegin;
1310:   if (y != z) PetscCall(VecCopy(y, z));
1311:   PetscCall(MatHYPRE_MultKernel_Private(A, 1.0, x, 1.0, z, PETSC_TRUE));
1312:   PetscFunctionReturn(PETSC_SUCCESS);
1313: }

1315: /* y = a * A * x + b * y or y = a * A^t * x + b * y depending on trans */
1316: static PetscErrorCode MatHYPRE_MultKernel_Private(Mat A, HYPRE_Complex a, Vec x, HYPRE_Complex b, Vec y, PetscBool trans)
1317: {
1318:   Mat_HYPRE          *hA = (Mat_HYPRE *)A->data;
1319:   hypre_ParCSRMatrix *parcsr;
1320:   hypre_ParVector    *hx, *hy;

1322:   PetscFunctionBegin;
1323:   if (trans) {
1324:     PetscCall(VecHYPRE_IJVectorPushVecRead(hA->b, x));
1325:     if (b != 0.0) PetscCall(VecHYPRE_IJVectorPushVec(hA->x, y));
1326:     else PetscCall(VecHYPRE_IJVectorPushVecWrite(hA->x, y));
1327:     PetscCallExternal(HYPRE_IJVectorGetObject, hA->b->ij, (void **)&hx);
1328:     PetscCallExternal(HYPRE_IJVectorGetObject, hA->x->ij, (void **)&hy);
1329:   } else {
1330:     PetscCall(VecHYPRE_IJVectorPushVecRead(hA->x, x));
1331:     if (b != 0.0) PetscCall(VecHYPRE_IJVectorPushVec(hA->b, y));
1332:     else PetscCall(VecHYPRE_IJVectorPushVecWrite(hA->b, y));
1333:     PetscCallExternal(HYPRE_IJVectorGetObject, hA->x->ij, (void **)&hx);
1334:     PetscCallExternal(HYPRE_IJVectorGetObject, hA->b->ij, (void **)&hy);
1335:   }
1336:   PetscCallExternal(HYPRE_IJMatrixGetObject, hA->ij, (void **)&parcsr);
1337:   if (trans) {
1338:     PetscCallExternal(hypre_ParCSRMatrixMatvecT, a, parcsr, hx, b, hy);
1339:   } else {
1340:     PetscCallExternal(hypre_ParCSRMatrixMatvec, a, parcsr, hx, b, hy);
1341:   }
1342:   PetscCall(VecHYPRE_IJVectorPopVec(hA->x));
1343:   PetscCall(VecHYPRE_IJVectorPopVec(hA->b));
1344:   PetscFunctionReturn(PETSC_SUCCESS);
1345: }

1347: static PetscErrorCode MatDestroy_HYPRE(Mat A)
1348: {
1349:   Mat_HYPRE *hA = (Mat_HYPRE *)A->data;

1351:   PetscFunctionBegin;
1352:   PetscCall(VecHYPRE_IJVectorDestroy(&hA->x));
1353:   PetscCall(VecHYPRE_IJVectorDestroy(&hA->b));
1354:   PetscCall(MatHYPRE_DestroyCOOMat(A)); /* must be called before destroying the individual CSR */
1355:   if (hA->ij) {
1356:     if (!hA->inner_free) hypre_IJMatrixObject(hA->ij) = NULL;
1357:     PetscCallExternal(HYPRE_IJMatrixDestroy, hA->ij);
1358:   }
1359:   if (hA->comm) PetscCall(PetscCommRestoreComm(PetscObjectComm((PetscObject)A), &hA->comm));

1361:   PetscCall(MatStashDestroy_Private(&A->stash));
1362:   PetscCall(PetscFree(hA->array));
1363:   if (hA->rows_d) PetscStackCallExternalVoid("hypre_Free", hypre_Free(hA->rows_d, HYPRE_MEMORY_DEVICE));

1365:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_hypre_aij_C", NULL));
1366:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatConvert_hypre_is_C", NULL));
1367:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_seqaij_hypre_C", NULL));
1368:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_mpiaij_hypre_C", NULL));
1369:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_seqaijhipsparse_hypre_C", NULL));
1370:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_mpiaijhipsparse_hypre_C", NULL));
1371:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_seqaijcusparse_hypre_C", NULL));
1372:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_mpiaijcusparse_hypre_C", NULL));
1373:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHYPRESetPreallocation_C", NULL));
1374:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatHYPREGetParCSR_C", NULL));
1375:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSetPreallocationCOO_C", NULL));
1376:   PetscCall(PetscObjectComposeFunction((PetscObject)A, "MatSetValuesCOO_C", NULL));
1377:   PetscCall(PetscFree(A->data));
1378:   PetscFunctionReturn(PETSC_SUCCESS);
1379: }

1381: static PetscErrorCode MatSetUp_HYPRE(Mat A)
1382: {
1383:   PetscFunctionBegin;
1384:   if (!A->preallocated) PetscCall(MatHYPRESetPreallocation(A, PETSC_DEFAULT, NULL, PETSC_DEFAULT, NULL));
1385:   PetscFunctionReturn(PETSC_SUCCESS);
1386: }

1388: //TODO FIX hypre_CSRMatrixMatvecOutOfPlace
1389: #if defined(PETSC_HAVE_HYPRE_DEVICE)
1390: static PetscErrorCode MatBindToCPU_HYPRE(Mat A, PetscBool bind)
1391: {
1392:   Mat_HYPRE           *hA   = (Mat_HYPRE *)A->data;
1393:   HYPRE_MemoryLocation hmem = bind ? HYPRE_MEMORY_HOST : HYPRE_MEMORY_DEVICE;

1395:   PetscFunctionBegin;
1396:   A->boundtocpu = bind;
1397:   if (hA->ij && hypre_IJMatrixAssembleFlag(hA->ij) && hmem != hypre_IJMatrixMemoryLocation(hA->ij)) {
1398:     hypre_ParCSRMatrix *parcsr;
1399:     PetscCallExternal(HYPRE_IJMatrixGetObject, hA->ij, (void **)&parcsr);
1400:     PetscCallExternal(hypre_ParCSRMatrixMigrate, parcsr, hmem);
1401:   }
1402:   if (hA->x) PetscCall(VecHYPRE_IJBindToCPU(hA->x, bind));
1403:   if (hA->b) PetscCall(VecHYPRE_IJBindToCPU(hA->b, bind));
1404:   PetscFunctionReturn(PETSC_SUCCESS);
1405: }
1406: #endif

1408: static PetscErrorCode MatAssemblyEnd_HYPRE(Mat A, MatAssemblyType mode)
1409: {
1410:   Mat_HYPRE   *hA = (Mat_HYPRE *)A->data;
1411:   PetscMPIInt  n;
1412:   PetscInt     i, j, rstart, ncols, flg;
1413:   PetscInt    *row, *col;
1414:   PetscScalar *val;

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

1419:   if (!A->nooffprocentries) {
1420:     while (1) {
1421:       PetscCall(MatStashScatterGetMesg_Private(&A->stash, &n, &row, &col, &val, &flg));
1422:       if (!flg) break;

1424:       for (i = 0; i < n;) {
1425:         /* Now identify the consecutive vals belonging to the same row */
1426:         for (j = i, rstart = row[j]; j < n; j++) {
1427:           if (row[j] != rstart) break;
1428:         }
1429:         if (j < n) ncols = j - i;
1430:         else ncols = n - i;
1431:         /* Now assemble all these values with a single function call */
1432:         PetscCall(MatSetValues_HYPRE(A, 1, row + i, ncols, col + i, val + i, A->insertmode));

1434:         i = j;
1435:       }
1436:     }
1437:     PetscCall(MatStashScatterEnd_Private(&A->stash));
1438:   }

1440:   PetscCallExternal(HYPRE_IJMatrixAssemble, hA->ij);
1441:   /* The assembly routine destroys the aux_matrix, we recreate it here by calling HYPRE_IJMatrixInitialize */
1442:   /* 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 */
1443:   if (!A->sortedfull) {
1444:     hypre_AuxParCSRMatrix *aux_matrix;

1446:     /* call destroy just to make sure we do not leak anything */
1447:     aux_matrix = (hypre_AuxParCSRMatrix *)hypre_IJMatrixTranslator(hA->ij);
1448:     PetscCallExternal(hypre_AuxParCSRMatrixDestroy, aux_matrix);
1449:     hypre_IJMatrixTranslator(hA->ij) = NULL;

1451:     /* Initialize with assembled flag -> it only recreates the aux_par_matrix */
1452:     PetscCallExternal(HYPRE_IJMatrixInitialize, hA->ij);
1453:     aux_matrix = (hypre_AuxParCSRMatrix *)hypre_IJMatrixTranslator(hA->ij);
1454:     if (aux_matrix) {
1455:       hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 1; /* see comment in MatHYPRESetPreallocation_HYPRE */
1456: #if PETSC_PKG_HYPRE_VERSION_LT(2, 19, 0)
1457:       PetscCallExternal(hypre_AuxParCSRMatrixInitialize, aux_matrix);
1458: #else
1459:       PetscCallExternal(hypre_AuxParCSRMatrixInitialize_v2, aux_matrix, HYPRE_MEMORY_HOST);
1460: #endif
1461:     }
1462:   }
1463:   {
1464:     hypre_ParCSRMatrix *parcsr;

1466:     PetscCallExternal(HYPRE_IJMatrixGetObject, hA->ij, (void **)&parcsr);
1467:     if (!hypre_ParCSRMatrixCommPkg(parcsr)) PetscCallExternal(hypre_MatvecCommPkgCreate, parcsr);
1468:   }
1469:   if (!hA->x) PetscCall(VecHYPRE_IJVectorCreate(A->cmap, &hA->x));
1470:   if (!hA->b) PetscCall(VecHYPRE_IJVectorCreate(A->rmap, &hA->b));
1471: #if defined(PETSC_HAVE_HYPRE_DEVICE)
1472:   PetscCall(MatBindToCPU_HYPRE(A, A->boundtocpu));
1473: #endif
1474:   PetscFunctionReturn(PETSC_SUCCESS);
1475: }

1477: static PetscErrorCode MatGetArray_HYPRE(Mat A, PetscInt size, void **array)
1478: {
1479:   Mat_HYPRE *hA = (Mat_HYPRE *)A->data;

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

1484:   if (hA->array_size >= size) {
1485:     *array = hA->array;
1486:   } else {
1487:     PetscCall(PetscFree(hA->array));
1488:     hA->array_size = size;
1489:     PetscCall(PetscMalloc(hA->array_size, &hA->array));
1490:     *array = hA->array;
1491:   }

1493:   hA->array_available = PETSC_FALSE;
1494:   PetscFunctionReturn(PETSC_SUCCESS);
1495: }

1497: static PetscErrorCode MatRestoreArray_HYPRE(Mat A, void **array)
1498: {
1499:   Mat_HYPRE *hA = (Mat_HYPRE *)A->data;

1501:   PetscFunctionBegin;
1502:   *array              = NULL;
1503:   hA->array_available = PETSC_TRUE;
1504:   PetscFunctionReturn(PETSC_SUCCESS);
1505: }

1507: static PetscErrorCode MatSetValues_HYPRE(Mat A, PetscInt nr, const PetscInt rows[], PetscInt nc, const PetscInt cols[], const PetscScalar v[], InsertMode ins)
1508: {
1509:   Mat_HYPRE     *hA   = (Mat_HYPRE *)A->data;
1510:   PetscScalar   *vals = (PetscScalar *)v;
1511:   HYPRE_Complex *sscr;
1512:   PetscInt      *cscr[2];
1513:   PetscInt       i, nzc;
1514:   PetscInt       rst = A->rmap->rstart, ren = A->rmap->rend;
1515:   void          *array = NULL;

1517:   PetscFunctionBegin;
1518:   PetscCall(MatGetArray_HYPRE(A, sizeof(PetscInt) * (2 * nc) + sizeof(HYPRE_Complex) * nc * nr, &array));
1519:   cscr[0] = (PetscInt *)array;
1520:   cscr[1] = ((PetscInt *)array) + nc;
1521:   sscr    = (HYPRE_Complex *)(((PetscInt *)array) + nc * 2);
1522:   for (i = 0, nzc = 0; i < nc; i++) {
1523:     if (cols[i] >= 0) {
1524:       cscr[0][nzc]   = cols[i];
1525:       cscr[1][nzc++] = i;
1526:     }
1527:   }
1528:   if (!nzc) {
1529:     PetscCall(MatRestoreArray_HYPRE(A, &array));
1530:     PetscFunctionReturn(PETSC_SUCCESS);
1531:   }

1533: #if 0 //defined(PETSC_HAVE_HYPRE_DEVICE)
1534:   if (HYPRE_MEMORY_HOST != hypre_IJMatrixMemoryLocation(hA->ij)) {
1535:     hypre_ParCSRMatrix *parcsr;

1537:     PetscCallExternal(HYPRE_IJMatrixGetObject,hA->ij,(void**)&parcsr);
1538:     PetscCallExternal(hypre_ParCSRMatrixMigrate,parcsr, HYPRE_MEMORY_HOST);
1539:   }
1540: #endif

1542:   if (ins == ADD_VALUES) {
1543:     for (i = 0; i < nr; i++) {
1544:       if (rows[i] >= 0) {
1545:         PetscInt  j;
1546:         HYPRE_Int hnc = (HYPRE_Int)nzc;

1548:         if (!nzc) continue;
1549:         /* nonlocal values */
1550:         if (rows[i] < rst || rows[i] >= ren) {
1551:           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]);
1552:           if (hA->donotstash) continue;
1553:         }
1554:         PetscCheck((PetscInt)hnc == nzc, PETSC_COMM_SELF, PETSC_ERR_SUP, "Hypre overflow! number of columns %" PetscInt_FMT " for row %" PetscInt_FMT, nzc, rows[i]);
1555:         for (j = 0; j < nzc; j++) PetscCall(PetscHYPREScalarCast(vals[cscr[1][j]], &sscr[j]));
1556:         PetscCallExternal(HYPRE_IJMatrixAddToValues, hA->ij, 1, &hnc, (HYPRE_BigInt *)(rows + i), (HYPRE_BigInt *)cscr[0], sscr);
1557:       }
1558:       vals += nc;
1559:     }
1560:   } else { /* INSERT_VALUES */
1561:     for (i = 0; i < nr; i++) {
1562:       if (rows[i] >= 0) {
1563:         PetscInt  j;
1564:         HYPRE_Int hnc = (HYPRE_Int)nzc;

1566:         if (!nzc) continue;
1567:         PetscCheck((PetscInt)hnc == nzc, PETSC_COMM_SELF, PETSC_ERR_SUP, "Hypre overflow! number of columns %" PetscInt_FMT " for row %" PetscInt_FMT, nzc, rows[i]);
1568:         for (j = 0; j < nzc; j++) PetscCall(PetscHYPREScalarCast(vals[cscr[1][j]], &sscr[j]));
1569:         /* nonlocal values */
1570:         if (rows[i] < rst || rows[i] >= ren) {
1571:           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]);
1572:           if (!hA->donotstash) PetscCall(MatStashValuesRow_Private(&A->stash, rows[i], nzc, cscr[0], (PetscScalar *)sscr, PETSC_FALSE));
1573:         }
1574:         /* local values */
1575:         else
1576:           PetscCallExternal(HYPRE_IJMatrixSetValues, hA->ij, 1, &hnc, (HYPRE_BigInt *)(rows + i), (HYPRE_BigInt *)cscr[0], sscr);
1577:       }
1578:       vals += nc;
1579:     }
1580:   }

1582:   PetscCall(MatRestoreArray_HYPRE(A, &array));
1583:   PetscFunctionReturn(PETSC_SUCCESS);
1584: }

1586: static PetscErrorCode MatHYPRESetPreallocation_HYPRE(Mat A, PetscInt dnz, const PetscInt dnnz[], PetscInt onz, const PetscInt onnz[])
1587: {
1588:   Mat_HYPRE  *hA = (Mat_HYPRE *)A->data;
1589:   HYPRE_Int  *hdnnz, *honnz;
1590:   PetscInt    i, rs, re, cs, ce, bs;
1591:   PetscMPIInt size;

1593:   PetscFunctionBegin;
1594:   PetscCall(PetscLayoutSetUp(A->rmap));
1595:   PetscCall(PetscLayoutSetUp(A->cmap));
1596:   rs = A->rmap->rstart;
1597:   re = A->rmap->rend;
1598:   cs = A->cmap->rstart;
1599:   ce = A->cmap->rend;
1600:   if (!hA->ij) {
1601:     PetscCallExternal(HYPRE_IJMatrixCreate, hA->comm, rs, re - 1, cs, ce - 1, &hA->ij);
1602:     PetscCallExternal(HYPRE_IJMatrixSetObjectType, hA->ij, HYPRE_PARCSR);
1603:   } else {
1604:     HYPRE_BigInt hrs, hre, hcs, hce;
1605:     PetscCallExternal(HYPRE_IJMatrixGetLocalRange, hA->ij, &hrs, &hre, &hcs, &hce);
1606:     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);
1607:     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);
1608:   }
1609:   PetscCall(MatHYPRE_DestroyCOOMat(A));
1610:   PetscCall(MatGetBlockSize(A, &bs));
1611:   if (dnz == PETSC_DEFAULT || dnz == PETSC_DECIDE) dnz = 10 * bs;
1612:   if (onz == PETSC_DEFAULT || onz == PETSC_DECIDE) onz = 10 * bs;

1614:   if (!dnnz) {
1615:     PetscCall(PetscMalloc1(A->rmap->n, &hdnnz));
1616:     for (i = 0; i < A->rmap->n; i++) hdnnz[i] = dnz;
1617:   } else {
1618:     hdnnz = (HYPRE_Int *)dnnz;
1619:   }
1620:   PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
1621:   if (size > 1) {
1622:     hypre_AuxParCSRMatrix *aux_matrix;
1623:     if (!onnz) {
1624:       PetscCall(PetscMalloc1(A->rmap->n, &honnz));
1625:       for (i = 0; i < A->rmap->n; i++) honnz[i] = onz;
1626:     } else honnz = (HYPRE_Int *)onnz;
1627:     /* SetDiagOffdSizes sets hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 0, since it seems
1628:        they assume the user will input the entire row values, properly sorted
1629:        In PETSc, we don't make such an assumption and set this flag to 1,
1630:        unless the option MAT_SORTED_FULL is set to true.
1631:        Also, to avoid possible memory leaks, we destroy and recreate the translator
1632:        This has to be done here, as HYPRE_IJMatrixInitialize will properly initialize
1633:        the IJ matrix for us */
1634:     aux_matrix = (hypre_AuxParCSRMatrix *)hypre_IJMatrixTranslator(hA->ij);
1635:     hypre_AuxParCSRMatrixDestroy(aux_matrix);
1636:     hypre_IJMatrixTranslator(hA->ij) = NULL;
1637:     PetscCallExternal(HYPRE_IJMatrixSetDiagOffdSizes, hA->ij, hdnnz, honnz);
1638:     aux_matrix                               = (hypre_AuxParCSRMatrix *)hypre_IJMatrixTranslator(hA->ij);
1639:     hypre_AuxParCSRMatrixNeedAux(aux_matrix) = !A->sortedfull;
1640:   } else {
1641:     honnz = NULL;
1642:     PetscCallExternal(HYPRE_IJMatrixSetRowSizes, hA->ij, hdnnz);
1643:   }

1645:   /* reset assembled flag and call the initialize method */
1646:   hypre_IJMatrixAssembleFlag(hA->ij) = 0;
1647: #if PETSC_PKG_HYPRE_VERSION_LT(2, 19, 0)
1648:   PetscCallExternal(HYPRE_IJMatrixInitialize, hA->ij);
1649: #else
1650:   PetscCallExternal(HYPRE_IJMatrixInitialize_v2, hA->ij, HYPRE_MEMORY_HOST);
1651: #endif
1652:   if (!dnnz) PetscCall(PetscFree(hdnnz));
1653:   if (!onnz && honnz) PetscCall(PetscFree(honnz));
1654:   /* Match AIJ logic */
1655:   A->preallocated = PETSC_TRUE;
1656:   A->assembled    = PETSC_FALSE;
1657:   PetscFunctionReturn(PETSC_SUCCESS);
1658: }

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

1663:   Collective

1665:   Input Parameters:
1666: + A    - the matrix
1667: . dnz  - number of nonzeros per row in DIAGONAL portion of local submatrix
1668:           (same value is used for all local rows)
1669: . dnnz - array containing the number of nonzeros in the various rows of the
1670:           DIAGONAL portion of the local submatrix (possibly different for each row)
1671:           or `NULL` (`PETSC_NULL_INTEGER` in Fortran), if `d_nz` is used to specify the nonzero structure.
1672:           The size of this array is equal to the number of local rows, i.e `m`.
1673:           For matrices that will be factored, you must leave room for (and set)
1674:           the diagonal entry even if it is zero.
1675: . onz  - number of nonzeros per row in the OFF-DIAGONAL portion of local
1676:           submatrix (same value is used for all local rows).
1677: - onnz - array containing the number of nonzeros in the various rows of the
1678:           OFF-DIAGONAL portion of the local submatrix (possibly different for
1679:           each row) or `NULL` (`PETSC_NULL_INTEGER` in Fortran), if `o_nz` is used to specify the nonzero
1680:           structure. The size of this array is equal to the number
1681:           of local rows, i.e `m`.

1683:   Level: intermediate

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

1688: .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MatMPIAIJSetPreallocation()`, `MATHYPRE`, `MATAIJ`
1689: @*/
1690: PetscErrorCode MatHYPRESetPreallocation(Mat A, PetscInt dnz, const PetscInt dnnz[], PetscInt onz, const PetscInt onnz[])
1691: {
1692:   PetscFunctionBegin;
1695:   PetscTryMethod(A, "MatHYPRESetPreallocation_C", (Mat, PetscInt, const PetscInt[], PetscInt, const PetscInt[]), (A, dnz, dnnz, onz, onnz));
1696:   PetscFunctionReturn(PETSC_SUCCESS);
1697: }

1699: /*@C
1700:   MatCreateFromParCSR - Creates a `Mat` from a `hypre_ParCSRMatrix`

1702:   Collective

1704:   Input Parameters:
1705: + parcsr   - the pointer to the `hypre_ParCSRMatrix`
1706: . mtype    - matrix type to be created. Currently `MATAIJ`, `MATIS` and `MATHYPRE` are supported.
1707: - copymode - PETSc copying options, see  `PetscCopyMode`

1709:   Output Parameter:
1710: . A - the matrix

1712:   Level: intermediate

1714: .seealso: [](ch_matrices), `Mat`, `MatHYPRE`, `PetscCopyMode`
1715: @*/
1716: PETSC_EXTERN PetscErrorCode MatCreateFromParCSR(hypre_ParCSRMatrix *parcsr, MatType mtype, PetscCopyMode copymode, Mat *A)
1717: {
1718:   Mat        T;
1719:   Mat_HYPRE *hA;
1720:   MPI_Comm   comm;
1721:   PetscInt   rstart, rend, cstart, cend, M, N;
1722:   PetscBool  isseqaij, isseqaijmkl, ismpiaij, isaij, ishyp, isis;

1724:   PetscFunctionBegin;
1725:   comm = hypre_ParCSRMatrixComm(parcsr);
1726:   PetscCall(PetscStrcmp(mtype, MATSEQAIJ, &isseqaij));
1727:   PetscCall(PetscStrcmp(mtype, MATSEQAIJMKL, &isseqaijmkl));
1728:   PetscCall(PetscStrcmp(mtype, MATMPIAIJ, &ismpiaij));
1729:   PetscCall(PetscStrcmp(mtype, MATAIJ, &isaij));
1730:   PetscCall(PetscStrcmp(mtype, MATHYPRE, &ishyp));
1731:   PetscCall(PetscStrcmp(mtype, MATIS, &isis));
1732:   isaij = (PetscBool)(isseqaij || isseqaijmkl || ismpiaij || isaij);
1733:   /* TODO */
1734:   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);
1735:   /* access ParCSRMatrix */
1736:   rstart = hypre_ParCSRMatrixFirstRowIndex(parcsr);
1737:   rend   = hypre_ParCSRMatrixLastRowIndex(parcsr);
1738:   cstart = hypre_ParCSRMatrixFirstColDiag(parcsr);
1739:   cend   = hypre_ParCSRMatrixLastColDiag(parcsr);
1740:   M      = hypre_ParCSRMatrixGlobalNumRows(parcsr);
1741:   N      = hypre_ParCSRMatrixGlobalNumCols(parcsr);

1743:   /* fix for empty local rows/columns */
1744:   if (rend < rstart) rend = rstart;
1745:   if (cend < cstart) cend = cstart;

1747:   /* PETSc convention */
1748:   rend++;
1749:   cend++;
1750:   rend = PetscMin(rend, M);
1751:   cend = PetscMin(cend, N);

1753:   /* create PETSc matrix with MatHYPRE */
1754:   PetscCall(MatCreate(comm, &T));
1755:   PetscCall(MatSetSizes(T, rend - rstart, cend - cstart, M, N));
1756:   PetscCall(MatSetType(T, MATHYPRE));
1757:   hA = (Mat_HYPRE *)T->data;

1759:   /* create HYPRE_IJMatrix */
1760:   PetscCallExternal(HYPRE_IJMatrixCreate, hA->comm, rstart, rend - 1, cstart, cend - 1, &hA->ij);
1761:   PetscCallExternal(HYPRE_IJMatrixSetObjectType, hA->ij, HYPRE_PARCSR);

1763:   /* create new ParCSR object if needed */
1764:   if (ishyp && copymode == PETSC_COPY_VALUES) {
1765:     hypre_ParCSRMatrix *new_parcsr;
1766: #if PETSC_PKG_HYPRE_VERSION_LT(2, 18, 0)
1767:     hypre_CSRMatrix *hdiag, *hoffd, *ndiag, *noffd;

1769:     new_parcsr = hypre_ParCSRMatrixClone(parcsr, 0);
1770:     hdiag      = hypre_ParCSRMatrixDiag(parcsr);
1771:     hoffd      = hypre_ParCSRMatrixOffd(parcsr);
1772:     ndiag      = hypre_ParCSRMatrixDiag(new_parcsr);
1773:     noffd      = hypre_ParCSRMatrixOffd(new_parcsr);
1774:     PetscCall(PetscArraycpy(hypre_CSRMatrixData(ndiag), hypre_CSRMatrixData(hdiag), hypre_CSRMatrixNumNonzeros(hdiag)));
1775:     PetscCall(PetscArraycpy(hypre_CSRMatrixData(noffd), hypre_CSRMatrixData(hoffd), hypre_CSRMatrixNumNonzeros(hoffd)));
1776: #else
1777:     new_parcsr = hypre_ParCSRMatrixClone(parcsr, 1);
1778: #endif
1779:     parcsr   = new_parcsr;
1780:     copymode = PETSC_OWN_POINTER;
1781:   }

1783:   /* set ParCSR object */
1784:   hypre_IJMatrixObject(hA->ij) = parcsr;
1785:   T->preallocated              = PETSC_TRUE;

1787:   /* set assembled flag */
1788:   hypre_IJMatrixAssembleFlag(hA->ij) = 1;
1789: #if 0
1790:   PetscCallExternal(HYPRE_IJMatrixInitialize,hA->ij);
1791: #endif
1792:   if (ishyp) {
1793:     PetscMPIInt myid = 0;

1795:     /* make sure we always have row_starts and col_starts available */
1796:     if (HYPRE_AssumedPartitionCheck()) PetscCallMPI(MPI_Comm_rank(comm, &myid));
1797: #if defined(hypre_ParCSRMatrixOwnsRowStarts)
1798:     if (!hypre_ParCSRMatrixOwnsColStarts(parcsr)) {
1799:       PetscLayout map;

1801:       PetscCall(MatGetLayouts(T, NULL, &map));
1802:       PetscCall(PetscLayoutSetUp(map));
1803:       hypre_ParCSRMatrixColStarts(parcsr) = (HYPRE_BigInt *)(map->range + myid);
1804:     }
1805:     if (!hypre_ParCSRMatrixOwnsRowStarts(parcsr)) {
1806:       PetscLayout map;

1808:       PetscCall(MatGetLayouts(T, &map, NULL));
1809:       PetscCall(PetscLayoutSetUp(map));
1810:       hypre_ParCSRMatrixRowStarts(parcsr) = (HYPRE_BigInt *)(map->range + myid);
1811:     }
1812: #endif
1813:     /* prevent from freeing the pointer */
1814:     if (copymode == PETSC_USE_POINTER) hA->inner_free = PETSC_FALSE;
1815:     *A = T;
1816:     PetscCall(MatSetOption(*A, MAT_SORTED_FULL, PETSC_TRUE));
1817:     PetscCall(MatAssemblyBegin(*A, MAT_FINAL_ASSEMBLY));
1818:     PetscCall(MatAssemblyEnd(*A, MAT_FINAL_ASSEMBLY));
1819:   } else if (isaij) {
1820:     if (copymode != PETSC_OWN_POINTER) {
1821:       /* prevent from freeing the pointer */
1822:       hA->inner_free = PETSC_FALSE;
1823:       PetscCall(MatConvert_HYPRE_AIJ(T, MATAIJ, MAT_INITIAL_MATRIX, A));
1824:       PetscCall(MatDestroy(&T));
1825:     } else { /* AIJ return type with PETSC_OWN_POINTER */
1826:       PetscCall(MatConvert_HYPRE_AIJ(T, MATAIJ, MAT_INPLACE_MATRIX, &T));
1827:       *A = T;
1828:     }
1829:   } else if (isis) {
1830:     PetscCall(MatConvert_HYPRE_IS(T, MATIS, MAT_INITIAL_MATRIX, A));
1831:     if (copymode != PETSC_OWN_POINTER) hA->inner_free = PETSC_FALSE;
1832:     PetscCall(MatDestroy(&T));
1833:   }
1834:   PetscFunctionReturn(PETSC_SUCCESS);
1835: }

1837: static PetscErrorCode MatHYPREGetParCSR_HYPRE(Mat A, hypre_ParCSRMatrix **parcsr)
1838: {
1839:   Mat_HYPRE *hA = (Mat_HYPRE *)A->data;
1840:   HYPRE_Int  type;

1842:   PetscFunctionBegin;
1843:   PetscCheck(hA->ij, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "HYPRE_IJMatrix not present");
1844:   PetscCallExternal(HYPRE_IJMatrixGetObjectType, hA->ij, &type);
1845:   PetscCheck(type == HYPRE_PARCSR, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "HYPRE_IJMatrix is not of type HYPRE_PARCSR");
1846:   PetscCallExternal(HYPRE_IJMatrixGetObject, hA->ij, (void **)parcsr);
1847:   PetscFunctionReturn(PETSC_SUCCESS);
1848: }

1850: /*@C
1851:   MatHYPREGetParCSR - Gets the pointer to the ParCSR matrix

1853:   Not Collective, No Fortran Support

1855:   Input Parameter:
1856: . A - the `MATHYPRE` object

1858:   Output Parameter:
1859: . parcsr - the pointer to the `hypre_ParCSRMatrix`

1861:   Level: intermediate

1863: .seealso: [](ch_matrices), `Mat`, `MatHYPRE`, `PetscCopyMode`
1864: @*/
1865: PetscErrorCode MatHYPREGetParCSR(Mat A, hypre_ParCSRMatrix **parcsr)
1866: {
1867:   PetscFunctionBegin;
1870:   PetscUseMethod(A, "MatHYPREGetParCSR_C", (Mat, hypre_ParCSRMatrix **), (A, parcsr));
1871:   PetscFunctionReturn(PETSC_SUCCESS);
1872: }

1874: static PetscErrorCode MatMissingDiagonal_HYPRE(Mat A, PetscBool *missing, PetscInt *dd)
1875: {
1876:   hypre_ParCSRMatrix *parcsr;
1877:   hypre_CSRMatrix    *ha;
1878:   PetscInt            rst;

1880:   PetscFunctionBegin;
1881:   PetscCheck(A->rmap->n == A->cmap->n, PETSC_COMM_SELF, PETSC_ERR_SUP, "Not implemented with non-square diagonal blocks");
1882:   PetscCall(MatGetOwnershipRange(A, &rst, NULL));
1883:   PetscCall(MatHYPREGetParCSR_HYPRE(A, &parcsr));
1884:   if (missing) *missing = PETSC_FALSE;
1885:   if (dd) *dd = -1;
1886:   ha = hypre_ParCSRMatrixDiag(parcsr);
1887:   if (ha) {
1888:     PetscInt   size, i;
1889:     HYPRE_Int *ii, *jj;

1891:     size = hypre_CSRMatrixNumRows(ha);
1892:     ii   = hypre_CSRMatrixI(ha);
1893:     jj   = hypre_CSRMatrixJ(ha);
1894:     for (i = 0; i < size; i++) {
1895:       PetscInt  j;
1896:       PetscBool found = PETSC_FALSE;

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

1900:       if (!found) {
1901:         PetscCall(PetscInfo(A, "Matrix is missing local diagonal entry %" PetscInt_FMT "\n", i));
1902:         if (missing) *missing = PETSC_TRUE;
1903:         if (dd) *dd = i + rst;
1904:         PetscFunctionReturn(PETSC_SUCCESS);
1905:       }
1906:     }
1907:     if (!size) {
1908:       PetscCall(PetscInfo(A, "Matrix has no diagonal entries therefore is missing diagonal\n"));
1909:       if (missing) *missing = PETSC_TRUE;
1910:       if (dd) *dd = rst;
1911:     }
1912:   } else {
1913:     PetscCall(PetscInfo(A, "Matrix has no diagonal entries therefore is missing diagonal\n"));
1914:     if (missing) *missing = PETSC_TRUE;
1915:     if (dd) *dd = rst;
1916:   }
1917:   PetscFunctionReturn(PETSC_SUCCESS);
1918: }

1920: static PetscErrorCode MatScale_HYPRE(Mat A, PetscScalar s)
1921: {
1922:   hypre_ParCSRMatrix *parcsr;
1923: #if PETSC_PKG_HYPRE_VERSION_LT(2, 19, 0)
1924:   hypre_CSRMatrix *ha;
1925: #endif
1926:   HYPRE_Complex hs;

1928:   PetscFunctionBegin;
1929:   PetscCall(PetscHYPREScalarCast(s, &hs));
1930:   PetscCall(MatHYPREGetParCSR_HYPRE(A, &parcsr));
1931: #if PETSC_PKG_HYPRE_VERSION_GE(2, 19, 0)
1932:   PetscCallExternal(hypre_ParCSRMatrixScale, parcsr, hs);
1933: #else /* diagonal part */
1934:   ha = hypre_ParCSRMatrixDiag(parcsr);
1935:   if (ha) {
1936:     PetscInt       size, i;
1937:     HYPRE_Int     *ii;
1938:     HYPRE_Complex *a;

1940:     size = hypre_CSRMatrixNumRows(ha);
1941:     a    = hypre_CSRMatrixData(ha);
1942:     ii   = hypre_CSRMatrixI(ha);
1943:     for (i = 0; i < ii[size]; i++) a[i] *= hs;
1944:   }
1945:   /* off-diagonal part */
1946:   ha = hypre_ParCSRMatrixOffd(parcsr);
1947:   if (ha) {
1948:     PetscInt       size, i;
1949:     HYPRE_Int     *ii;
1950:     HYPRE_Complex *a;

1952:     size = hypre_CSRMatrixNumRows(ha);
1953:     a    = hypre_CSRMatrixData(ha);
1954:     ii   = hypre_CSRMatrixI(ha);
1955:     for (i = 0; i < ii[size]; i++) a[i] *= hs;
1956:   }
1957: #endif
1958:   PetscFunctionReturn(PETSC_SUCCESS);
1959: }

1961: static PetscErrorCode MatZeroRowsColumns_HYPRE(Mat A, PetscInt numRows, const PetscInt rows[], PetscScalar diag, Vec x, Vec b)
1962: {
1963:   hypre_ParCSRMatrix *parcsr;
1964:   HYPRE_Int          *lrows;
1965:   PetscInt            rst, ren, i;

1967:   PetscFunctionBegin;
1968:   PetscCheck(!x && !b, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "To be implemented");
1969:   PetscCall(MatHYPREGetParCSR_HYPRE(A, &parcsr));
1970:   PetscCall(PetscMalloc1(numRows, &lrows));
1971:   PetscCall(MatGetOwnershipRange(A, &rst, &ren));
1972:   for (i = 0; i < numRows; i++) {
1973:     PetscCheck(rows[i] >= rst && rows[i] < ren, PETSC_COMM_SELF, PETSC_ERR_SUP, "Non-local rows not yet supported");
1974:     lrows[i] = rows[i] - rst;
1975:   }
1976:   PetscCallExternal(hypre_ParCSRMatrixEliminateRowsCols, parcsr, numRows, lrows);
1977:   PetscCall(PetscFree(lrows));
1978:   PetscFunctionReturn(PETSC_SUCCESS);
1979: }

1981: static PetscErrorCode MatZeroEntries_HYPRE_CSRMatrix(hypre_CSRMatrix *ha)
1982: {
1983:   PetscFunctionBegin;
1984:   if (ha) {
1985:     HYPRE_Int     *ii, size;
1986:     HYPRE_Complex *a;

1988:     size = hypre_CSRMatrixNumRows(ha);
1989:     a    = hypre_CSRMatrixData(ha);
1990:     ii   = hypre_CSRMatrixI(ha);

1992:     if (a) PetscCall(PetscArrayzero(a, ii[size]));
1993:   }
1994:   PetscFunctionReturn(PETSC_SUCCESS);
1995: }

1997: static PetscErrorCode MatZeroEntries_HYPRE(Mat A)
1998: {
1999:   Mat_HYPRE *hA = (Mat_HYPRE *)A->data;

2001:   PetscFunctionBegin;
2002:   if (HYPRE_MEMORY_DEVICE == hypre_IJMatrixMemoryLocation(hA->ij)) {
2003:     PetscCallExternal(HYPRE_IJMatrixSetConstantValues, hA->ij, 0.0);
2004:   } else {
2005:     hypre_ParCSRMatrix *parcsr;

2007:     PetscCall(MatHYPREGetParCSR_HYPRE(A, &parcsr));
2008:     PetscCall(MatZeroEntries_HYPRE_CSRMatrix(hypre_ParCSRMatrixDiag(parcsr)));
2009:     PetscCall(MatZeroEntries_HYPRE_CSRMatrix(hypre_ParCSRMatrixOffd(parcsr)));
2010:   }
2011:   PetscFunctionReturn(PETSC_SUCCESS);
2012: }

2014: static PetscErrorCode MatZeroRows_HYPRE_CSRMatrix(hypre_CSRMatrix *hA, PetscInt N, const PetscInt rows[], HYPRE_Complex diag)
2015: {
2016:   PetscInt       ii;
2017:   HYPRE_Int     *i, *j;
2018:   HYPRE_Complex *a;

2020:   PetscFunctionBegin;
2021:   if (!hA) PetscFunctionReturn(PETSC_SUCCESS);

2023:   i = hypre_CSRMatrixI(hA);
2024:   j = hypre_CSRMatrixJ(hA);
2025:   a = hypre_CSRMatrixData(hA);
2026: #if defined(PETSC_HAVE_HYPRE_DEVICE)
2027:   if (HYPRE_MEMORY_DEVICE == hypre_CSRMatrixMemoryLocation(hA)) {
2028:   #if defined(HYPRE_USING_CUDA)
2029:     MatZeroRows_CUDA(N, rows, i, j, a, diag);
2030:   #elif defined(HYPRE_USING_HIP)
2031:     MatZeroRows_HIP(N, rows, i, j, a, diag);
2032:   #elif defined(PETSC_HAVE_KOKKOS)
2033:     MatZeroRows_Kokkos(N, rows, i, j, a, diag);
2034:   #else
2035:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for MatZeroRows on a hypre matrix in this memory location");
2036:   #endif
2037:   } else
2038: #endif
2039:   {
2040:     for (ii = 0; ii < N; ii++) {
2041:       HYPRE_Int jj, ibeg, iend, irow;

2043:       irow = rows[ii];
2044:       ibeg = i[irow];
2045:       iend = i[irow + 1];
2046:       for (jj = ibeg; jj < iend; jj++)
2047:         if (j[jj] == irow) a[jj] = diag;
2048:         else a[jj] = 0.0;
2049:     }
2050:   }
2051:   PetscFunctionReturn(PETSC_SUCCESS);
2052: }

2054: static PetscErrorCode MatZeroRows_HYPRE(Mat A, PetscInt N, const PetscInt rows[], PetscScalar diag, Vec x, Vec b)
2055: {
2056:   hypre_ParCSRMatrix *parcsr;
2057:   PetscInt           *lrows, len, *lrows2;
2058:   HYPRE_Complex       hdiag;

2060:   PetscFunctionBegin;
2061:   PetscCheck(!x && !b, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Does not support to modify the solution and the right hand size");
2062:   PetscCall(PetscHYPREScalarCast(diag, &hdiag));
2063:   /* retrieve the internal matrix */
2064:   PetscCall(MatHYPREGetParCSR_HYPRE(A, &parcsr));
2065:   /* get locally owned rows */
2066:   PetscCall(MatZeroRowsMapLocal_Private(A, N, rows, &len, &lrows));

2068: #if defined(PETSC_HAVE_HYPRE_DEVICE)
2069:   if (HYPRE_MEMORY_DEVICE == hypre_CSRMatrixMemoryLocation(hypre_ParCSRMatrixDiag(parcsr))) {
2070:     Mat_HYPRE *hA = (Mat_HYPRE *)A->data;
2071:     PetscInt   m;
2072:     PetscCall(MatGetLocalSize(A, &m, NULL));
2073:     if (!hA->rows_d) {
2074:       hA->rows_d = hypre_TAlloc(PetscInt, m, HYPRE_MEMORY_DEVICE);
2075:       if (m) PetscCheck(hA->rows_d, PETSC_COMM_SELF, PETSC_ERR_MEM, "HYPRE_TAlloc failed");
2076:     }
2077:     PetscCheck(len <= m, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Too many rows in rows[]");
2078:     PetscStackCallExternalVoid("hypre_Memcpy", hypre_Memcpy(hA->rows_d, lrows, sizeof(PetscInt) * len, HYPRE_MEMORY_DEVICE, HYPRE_MEMORY_HOST));
2079:     lrows2 = hA->rows_d;
2080:   } else
2081: #endif
2082:   {
2083:     lrows2 = lrows;
2084:   }

2086:   /* zero diagonal part */
2087:   PetscCall(MatZeroRows_HYPRE_CSRMatrix(hypre_ParCSRMatrixDiag(parcsr), len, lrows2, hdiag));
2088:   /* zero off-diagonal part */
2089:   PetscCall(MatZeroRows_HYPRE_CSRMatrix(hypre_ParCSRMatrixOffd(parcsr), len, lrows2, 0.0));

2091:   PetscCall(PetscFree(lrows));
2092:   PetscFunctionReturn(PETSC_SUCCESS);
2093: }

2095: static PetscErrorCode MatAssemblyBegin_HYPRE(Mat mat, MatAssemblyType mode)
2096: {
2097:   PetscFunctionBegin;
2098:   if (mat->nooffprocentries) PetscFunctionReturn(PETSC_SUCCESS);

2100:   PetscCall(MatStashScatterBegin_Private(mat, &mat->stash, mat->rmap->range));
2101:   PetscFunctionReturn(PETSC_SUCCESS);
2102: }

2104: static PetscErrorCode MatGetRow_HYPRE(Mat A, PetscInt row, PetscInt *nz, PetscInt **idx, PetscScalar **v)
2105: {
2106:   hypre_ParCSRMatrix *parcsr;
2107:   HYPRE_Int           hnz;

2109:   PetscFunctionBegin;
2110:   /* retrieve the internal matrix */
2111:   PetscCall(MatHYPREGetParCSR_HYPRE(A, &parcsr));
2112:   /* call HYPRE API */
2113:   PetscCallExternal(HYPRE_ParCSRMatrixGetRow, parcsr, row, &hnz, (HYPRE_BigInt **)idx, (HYPRE_Complex **)v);
2114:   if (nz) *nz = (PetscInt)hnz;
2115:   PetscFunctionReturn(PETSC_SUCCESS);
2116: }

2118: static PetscErrorCode MatRestoreRow_HYPRE(Mat A, PetscInt row, PetscInt *nz, PetscInt **idx, PetscScalar **v)
2119: {
2120:   hypre_ParCSRMatrix *parcsr;
2121:   HYPRE_Int           hnz;

2123:   PetscFunctionBegin;
2124:   /* retrieve the internal matrix */
2125:   PetscCall(MatHYPREGetParCSR_HYPRE(A, &parcsr));
2126:   /* call HYPRE API */
2127:   hnz = nz ? (HYPRE_Int)(*nz) : 0;
2128:   PetscCallExternal(HYPRE_ParCSRMatrixRestoreRow, parcsr, row, &hnz, (HYPRE_BigInt **)idx, (HYPRE_Complex **)v);
2129:   PetscFunctionReturn(PETSC_SUCCESS);
2130: }

2132: static PetscErrorCode MatGetValues_HYPRE(Mat A, PetscInt m, const PetscInt idxm[], PetscInt n, const PetscInt idxn[], PetscScalar v[])
2133: {
2134:   Mat_HYPRE *hA = (Mat_HYPRE *)A->data;
2135:   PetscInt   i;

2137:   PetscFunctionBegin;
2138:   if (!m || !n) PetscFunctionReturn(PETSC_SUCCESS);
2139:   /* Ignore negative row indices
2140:    * And negative column indices should be automatically ignored in hypre
2141:    * */
2142:   for (i = 0; i < m; i++) {
2143:     if (idxm[i] >= 0) {
2144:       HYPRE_Int hn = (HYPRE_Int)n;
2145:       PetscCallExternal(HYPRE_IJMatrixGetValues, hA->ij, 1, &hn, (HYPRE_BigInt *)&idxm[i], (HYPRE_BigInt *)idxn, (HYPRE_Complex *)(v + i * n));
2146:     }
2147:   }
2148:   PetscFunctionReturn(PETSC_SUCCESS);
2149: }

2151: static PetscErrorCode MatSetOption_HYPRE(Mat A, MatOption op, PetscBool flg)
2152: {
2153:   Mat_HYPRE *hA = (Mat_HYPRE *)A->data;

2155:   PetscFunctionBegin;
2156:   switch (op) {
2157:   case MAT_NO_OFF_PROC_ENTRIES:
2158:     if (flg) PetscCallExternal(HYPRE_IJMatrixSetMaxOffProcElmts, hA->ij, 0);
2159:     break;
2160:   case MAT_IGNORE_OFF_PROC_ENTRIES:
2161:     hA->donotstash = flg;
2162:     break;
2163:   default:
2164:     break;
2165:   }
2166:   PetscFunctionReturn(PETSC_SUCCESS);
2167: }

2169: static PetscErrorCode MatView_HYPRE(Mat A, PetscViewer view)
2170: {
2171:   PetscViewerFormat format;

2173:   PetscFunctionBegin;
2174:   PetscCall(PetscViewerGetFormat(view, &format));
2175:   if (format == PETSC_VIEWER_ASCII_FACTOR_INFO || format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) PetscFunctionReturn(PETSC_SUCCESS);
2176:   if (format != PETSC_VIEWER_NATIVE) {
2177:     Mat                 B;
2178:     hypre_ParCSRMatrix *parcsr;
2179:     PetscErrorCode (*mview)(Mat, PetscViewer) = NULL;

2181:     PetscCall(MatHYPREGetParCSR_HYPRE(A, &parcsr));
2182:     PetscCall(MatCreateFromParCSR(parcsr, MATAIJ, PETSC_USE_POINTER, &B));
2183:     PetscCall(MatGetOperation(B, MATOP_VIEW, (void (**)(void))&mview));
2184:     PetscCheck(mview, PetscObjectComm((PetscObject)A), PETSC_ERR_PLIB, "Missing view operation");
2185:     PetscCall((*mview)(B, view));
2186:     PetscCall(MatDestroy(&B));
2187:   } else {
2188:     Mat_HYPRE  *hA = (Mat_HYPRE *)A->data;
2189:     PetscMPIInt size;
2190:     PetscBool   isascii;
2191:     const char *filename;

2193:     /* HYPRE uses only text files */
2194:     PetscCall(PetscObjectTypeCompare((PetscObject)view, PETSCVIEWERASCII, &isascii));
2195:     PetscCheck(isascii, PetscObjectComm((PetscObject)view), PETSC_ERR_SUP, "PetscViewerType %s: native HYPRE format needs PETSCVIEWERASCII", ((PetscObject)view)->type_name);
2196:     PetscCall(PetscViewerFileGetName(view, &filename));
2197:     PetscCallExternal(HYPRE_IJMatrixPrint, hA->ij, filename);
2198:     PetscCallMPI(MPI_Comm_size(hA->comm, &size));
2199:     if (size > 1) {
2200:       PetscCall(PetscViewerASCIIPrintf(view, "Matrix files: %s.%05d ... %s.%05d\n", filename, 0, filename, size - 1));
2201:     } else {
2202:       PetscCall(PetscViewerASCIIPrintf(view, "Matrix file: %s.%05d\n", filename, 0));
2203:     }
2204:   }
2205:   PetscFunctionReturn(PETSC_SUCCESS);
2206: }

2208: static PetscErrorCode MatCopy_HYPRE(Mat A, Mat B, MatStructure str)
2209: {
2210:   hypre_ParCSRMatrix *acsr, *bcsr;

2212:   PetscFunctionBegin;
2213:   if (str == SAME_NONZERO_PATTERN && A->ops->copy == B->ops->copy) {
2214:     PetscCall(MatHYPREGetParCSR_HYPRE(A, &acsr));
2215:     PetscCall(MatHYPREGetParCSR_HYPRE(B, &bcsr));
2216:     PetscCallExternal(hypre_ParCSRMatrixCopy, acsr, bcsr, 1);
2217:     PetscCall(MatSetOption(B, MAT_SORTED_FULL, PETSC_TRUE)); /* "perfect" preallocation, so no need for hypre_AuxParCSRMatrixNeedAux */
2218:     PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
2219:     PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
2220:   } else {
2221:     PetscCall(MatCopy_Basic(A, B, str));
2222:   }
2223:   PetscFunctionReturn(PETSC_SUCCESS);
2224: }

2226: static PetscErrorCode MatGetDiagonal_HYPRE(Mat A, Vec d)
2227: {
2228:   hypre_ParCSRMatrix *parcsr;
2229:   hypre_CSRMatrix    *dmat;
2230:   HYPRE_Complex      *a;
2231:   PetscBool           cong;

2233:   PetscFunctionBegin;
2234:   PetscCall(MatHasCongruentLayouts(A, &cong));
2235:   PetscCheck(cong, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Only for square matrices with same local distributions of rows and columns");
2236:   PetscCall(MatHYPREGetParCSR_HYPRE(A, &parcsr));
2237:   dmat = hypre_ParCSRMatrixDiag(parcsr);
2238:   if (dmat) {
2239: #if defined(PETSC_HAVE_HYPRE_DEVICE)
2240:     HYPRE_MemoryLocation mem = hypre_CSRMatrixMemoryLocation(dmat);
2241: #else
2242:     HYPRE_MemoryLocation mem = HYPRE_MEMORY_HOST;
2243: #endif

2245:     if (mem != HYPRE_MEMORY_HOST) PetscCall(VecGetArrayWriteAndMemType(d, (PetscScalar **)&a, NULL));
2246:     else PetscCall(VecGetArrayWrite(d, (PetscScalar **)&a));
2247:     hypre_CSRMatrixExtractDiagonal(dmat, a, 0);
2248:     if (mem != HYPRE_MEMORY_HOST) PetscCall(VecRestoreArrayWriteAndMemType(d, (PetscScalar **)&a));
2249:     else PetscCall(VecRestoreArrayWrite(d, (PetscScalar **)&a));
2250:   }
2251:   PetscFunctionReturn(PETSC_SUCCESS);
2252: }

2254: #include <petscblaslapack.h>

2256: static PetscErrorCode MatAXPY_HYPRE(Mat Y, PetscScalar a, Mat X, MatStructure str)
2257: {
2258:   PetscFunctionBegin;
2259: #if defined(PETSC_HAVE_HYPRE_DEVICE)
2260:   {
2261:     Mat                 B;
2262:     hypre_ParCSRMatrix *x, *y, *z;

2264:     PetscCall(MatHYPREGetParCSR(Y, &y));
2265:     PetscCall(MatHYPREGetParCSR(X, &x));
2266:     PetscCallExternal(hypre_ParCSRMatrixAdd, 1.0, y, 1.0, x, &z);
2267:     PetscCall(MatCreateFromParCSR(z, MATHYPRE, PETSC_OWN_POINTER, &B));
2268:     PetscCall(MatHeaderMerge(Y, &B));
2269:   }
2270: #else
2271:   if (str == SAME_NONZERO_PATTERN) {
2272:     hypre_ParCSRMatrix *x, *y;
2273:     hypre_CSRMatrix    *xloc, *yloc;
2274:     PetscInt            xnnz, ynnz;
2275:     HYPRE_Complex      *xarr, *yarr;
2276:     PetscBLASInt        one = 1, bnz;

2278:     PetscCall(MatHYPREGetParCSR(Y, &y));
2279:     PetscCall(MatHYPREGetParCSR(X, &x));

2281:     /* diagonal block */
2282:     xloc = hypre_ParCSRMatrixDiag(x);
2283:     yloc = hypre_ParCSRMatrixDiag(y);
2284:     xnnz = 0;
2285:     ynnz = 0;
2286:     xarr = NULL;
2287:     yarr = NULL;
2288:     if (xloc) {
2289:       xarr = hypre_CSRMatrixData(xloc);
2290:       xnnz = hypre_CSRMatrixNumNonzeros(xloc);
2291:     }
2292:     if (yloc) {
2293:       yarr = hypre_CSRMatrixData(yloc);
2294:       ynnz = hypre_CSRMatrixNumNonzeros(yloc);
2295:     }
2296:     PetscCheck(xnnz == ynnz, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Different number of nonzeros in diagonal block %" PetscInt_FMT " != %" PetscInt_FMT, xnnz, ynnz);
2297:     PetscCall(PetscBLASIntCast(xnnz, &bnz));
2298:     PetscCallBLAS("BLASaxpy", BLASaxpy_(&bnz, &a, (PetscScalar *)xarr, &one, (PetscScalar *)yarr, &one));

2300:     /* off-diagonal block */
2301:     xloc = hypre_ParCSRMatrixOffd(x);
2302:     yloc = hypre_ParCSRMatrixOffd(y);
2303:     xnnz = 0;
2304:     ynnz = 0;
2305:     xarr = NULL;
2306:     yarr = NULL;
2307:     if (xloc) {
2308:       xarr = hypre_CSRMatrixData(xloc);
2309:       xnnz = hypre_CSRMatrixNumNonzeros(xloc);
2310:     }
2311:     if (yloc) {
2312:       yarr = hypre_CSRMatrixData(yloc);
2313:       ynnz = hypre_CSRMatrixNumNonzeros(yloc);
2314:     }
2315:     PetscCheck(xnnz == ynnz, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Different number of nonzeros in off-diagonal block %" PetscInt_FMT " != %" PetscInt_FMT, xnnz, ynnz);
2316:     PetscCall(PetscBLASIntCast(xnnz, &bnz));
2317:     PetscCallBLAS("BLASaxpy", BLASaxpy_(&bnz, &a, (PetscScalar *)xarr, &one, (PetscScalar *)yarr, &one));
2318:   } else if (str == SUBSET_NONZERO_PATTERN) {
2319:     PetscCall(MatAXPY_Basic(Y, a, X, str));
2320:   } else {
2321:     Mat B;

2323:     PetscCall(MatAXPY_Basic_Preallocate(Y, X, &B));
2324:     PetscCall(MatAXPY_BasicWithPreallocation(B, Y, a, X, str));
2325:     PetscCall(MatHeaderReplace(Y, &B));
2326:   }
2327: #endif
2328:   PetscFunctionReturn(PETSC_SUCCESS);
2329: }

2331: static PetscErrorCode MatDuplicate_HYPRE(Mat A, MatDuplicateOption op, Mat *B)
2332: {
2333:   hypre_ParCSRMatrix *parcsr = NULL;
2334:   PetscCopyMode       cpmode;
2335:   Mat_HYPRE          *hA;

2337:   PetscFunctionBegin;
2338:   PetscCall(MatHYPREGetParCSR_HYPRE(A, &parcsr));
2339:   if (op == MAT_DO_NOT_COPY_VALUES || op == MAT_SHARE_NONZERO_PATTERN) {
2340:     parcsr = hypre_ParCSRMatrixClone(parcsr, 0);
2341:     cpmode = PETSC_OWN_POINTER;
2342:   } else {
2343:     cpmode = PETSC_COPY_VALUES;
2344:   }
2345:   PetscCall(MatCreateFromParCSR(parcsr, MATHYPRE, cpmode, B));
2346:   hA = (Mat_HYPRE *)A->data;
2347:   if (hA->cooMat) {
2348:     Mat_HYPRE *hB = (Mat_HYPRE *)((*B)->data);
2349:     op            = (op == MAT_DO_NOT_COPY_VALUES) ? op : MAT_COPY_VALUES;
2350:     /* Cannot simply increase the reference count of hA->cooMat, since B needs to share cooMat's data array */
2351:     PetscCall(MatDuplicate(hA->cooMat, op, &hB->cooMat));
2352:     PetscCall(MatHYPRE_AttachCOOMat(*B));
2353:   }
2354:   PetscFunctionReturn(PETSC_SUCCESS);
2355: }

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

2361:   PetscFunctionBegin;
2362:   /* Build an agent matrix cooMat with AIJ format
2363:      It has the same sparsity pattern as mat, and also shares the data array with mat. We use cooMat to do the COO work.
2364:    */
2365:   PetscCall(MatHYPRE_CreateCOOMat(mat));
2366:   PetscCall(MatSetOption(hmat->cooMat, MAT_IGNORE_OFF_PROC_ENTRIES, hmat->donotstash));
2367:   PetscCall(MatSetOption(hmat->cooMat, MAT_NO_OFF_PROC_ENTRIES, mat->nooffprocentries));

2369:   /* MatSetPreallocationCOO_SeqAIJ and MatSetPreallocationCOO_MPIAIJ uses this specific
2370:      name to automatically put the diagonal entries first */
2371:   PetscCall(PetscObjectSetName((PetscObject)hmat->cooMat, "_internal_COO_mat_for_hypre"));
2372:   PetscCall(MatSetPreallocationCOO(hmat->cooMat, coo_n, coo_i, coo_j));
2373:   hmat->cooMat->assembled = PETSC_TRUE;

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

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

2384:   /* Attach cooMat to mat */
2385:   PetscCall(MatHYPRE_AttachCOOMat(mat));
2386:   PetscFunctionReturn(PETSC_SUCCESS);
2387: }

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

2393:   PetscFunctionBegin;
2394:   PetscCheck(hmat->cooMat, PetscObjectComm((PetscObject)mat), PETSC_ERR_PLIB, "HYPRE COO delegate matrix has not been created yet");
2395:   PetscCall(MatSetValuesCOO(hmat->cooMat, v, imode));
2396:   PetscCall(MatViewFromOptions(hmat->cooMat, (PetscObject)mat, "-cooMat_view"));
2397:   PetscFunctionReturn(PETSC_SUCCESS);
2398: }

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

2404:    Level: intermediate

2406: .seealso: [](ch_matrices), `Mat`, `MatCreate()`, `MatHYPRESetPreallocation`
2407: M*/

2409: PETSC_EXTERN PetscErrorCode MatCreate_HYPRE(Mat B)
2410: {
2411:   Mat_HYPRE *hB;
2412: #if defined(PETSC_HAVE_HYPRE_DEVICE)
2413:   HYPRE_MemoryLocation memory_location;
2414: #endif

2416:   PetscFunctionBegin;
2417:   PetscHYPREInitialize();
2418:   PetscCall(PetscNew(&hB));

2420:   hB->inner_free      = PETSC_TRUE;
2421:   hB->array_available = PETSC_TRUE;

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

2425:   PetscCall(PetscMemzero(B->ops, sizeof(struct _MatOps)));
2426:   B->ops->mult                  = MatMult_HYPRE;
2427:   B->ops->multtranspose         = MatMultTranspose_HYPRE;
2428:   B->ops->multadd               = MatMultAdd_HYPRE;
2429:   B->ops->multtransposeadd      = MatMultTransposeAdd_HYPRE;
2430:   B->ops->setup                 = MatSetUp_HYPRE;
2431:   B->ops->destroy               = MatDestroy_HYPRE;
2432:   B->ops->assemblyend           = MatAssemblyEnd_HYPRE;
2433:   B->ops->assemblybegin         = MatAssemblyBegin_HYPRE;
2434:   B->ops->setvalues             = MatSetValues_HYPRE;
2435:   B->ops->missingdiagonal       = MatMissingDiagonal_HYPRE;
2436:   B->ops->scale                 = MatScale_HYPRE;
2437:   B->ops->zerorowscolumns       = MatZeroRowsColumns_HYPRE;
2438:   B->ops->zeroentries           = MatZeroEntries_HYPRE;
2439:   B->ops->zerorows              = MatZeroRows_HYPRE;
2440:   B->ops->getrow                = MatGetRow_HYPRE;
2441:   B->ops->restorerow            = MatRestoreRow_HYPRE;
2442:   B->ops->getvalues             = MatGetValues_HYPRE;
2443:   B->ops->setoption             = MatSetOption_HYPRE;
2444:   B->ops->duplicate             = MatDuplicate_HYPRE;
2445:   B->ops->copy                  = MatCopy_HYPRE;
2446:   B->ops->view                  = MatView_HYPRE;
2447:   B->ops->getdiagonal           = MatGetDiagonal_HYPRE;
2448:   B->ops->axpy                  = MatAXPY_HYPRE;
2449:   B->ops->productsetfromoptions = MatProductSetFromOptions_HYPRE;
2450: #if defined(PETSC_HAVE_HYPRE_DEVICE)
2451:   B->ops->bindtocpu = MatBindToCPU_HYPRE;
2452:   /* Get hypre's default memory location. Users can control this using the corresponding HYPRE_SetMemoryLocation API */
2453:   PetscCallExternal(HYPRE_GetMemoryLocation, &memory_location);
2454:   B->boundtocpu = (memory_location == HYPRE_MEMORY_HOST) ? PETSC_TRUE : PETSC_FALSE;
2455: #endif

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

2460:   PetscCall(PetscCommGetComm(PetscObjectComm((PetscObject)B), &hB->comm));
2461:   PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATHYPRE));
2462:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_hypre_aij_C", MatConvert_HYPRE_AIJ));
2463:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatConvert_hypre_is_C", MatConvert_HYPRE_IS));
2464:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqaij_hypre_C", MatProductSetFromOptions_HYPRE));
2465:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_mpiaij_hypre_C", MatProductSetFromOptions_HYPRE));
2466:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatHYPRESetPreallocation_C", MatHYPRESetPreallocation_HYPRE));
2467:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatHYPREGetParCSR_C", MatHYPREGetParCSR_HYPRE));
2468:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSetPreallocationCOO_C", MatSetPreallocationCOO_HYPRE));
2469:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatSetValuesCOO_C", MatSetValuesCOO_HYPRE));
2470: #if defined(PETSC_HAVE_HYPRE_DEVICE)
2471:   #if defined(HYPRE_USING_HIP)
2472:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqaijhipsparse_hypre_C", MatProductSetFromOptions_HYPRE));
2473:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_mpiaijhipsparse_hypre_C", MatProductSetFromOptions_HYPRE));
2474:   PetscCall(PetscDeviceInitialize(PETSC_DEVICE_HIP));
2475:   PetscCall(MatSetVecType(B, VECHIP));
2476:   #endif
2477:   #if defined(HYPRE_USING_CUDA)
2478:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqaijcusparse_hypre_C", MatProductSetFromOptions_HYPRE));
2479:   PetscCall(PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_mpiaijcusparse_hypre_C", MatProductSetFromOptions_HYPRE));
2480:   PetscCall(PetscDeviceInitialize(PETSC_DEVICE_CUDA));
2481:   PetscCall(MatSetVecType(B, VECCUDA));
2482:   #endif
2483: #endif
2484:   PetscFunctionReturn(PETSC_SUCCESS);
2485: }