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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

639:     PetscCall(MatHYPRE_CreateCOOMat(M));

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1662:   Collective

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

1682:   Level: intermediate

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

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

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

1701:   Collective

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

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

1711:   Level: intermediate

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1852:   Not Collective

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

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

1860:   Level: intermediate

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2253: #include <petscblaslapack.h>

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2403:    Level: intermediate

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

2408: PETSC_EXTERN PetscErrorCode MatCreate_HYPRE(Mat B)
2409: {
2410:   Mat_HYPRE *hB;

2412:   PetscFunctionBegin;
2413:   PetscCall(PetscNew(&hB));

2415:   hB->inner_free      = PETSC_TRUE;
2416:   hB->array_available = PETSC_TRUE;

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

2420:   PetscCall(PetscMemzero(B->ops, sizeof(struct _MatOps)));
2421:   B->ops->mult                  = MatMult_HYPRE;
2422:   B->ops->multtranspose         = MatMultTranspose_HYPRE;
2423:   B->ops->multadd               = MatMultAdd_HYPRE;
2424:   B->ops->multtransposeadd      = MatMultTransposeAdd_HYPRE;
2425:   B->ops->setup                 = MatSetUp_HYPRE;
2426:   B->ops->destroy               = MatDestroy_HYPRE;
2427:   B->ops->assemblyend           = MatAssemblyEnd_HYPRE;
2428:   B->ops->assemblybegin         = MatAssemblyBegin_HYPRE;
2429:   B->ops->setvalues             = MatSetValues_HYPRE;
2430:   B->ops->missingdiagonal       = MatMissingDiagonal_HYPRE;
2431:   B->ops->scale                 = MatScale_HYPRE;
2432:   B->ops->zerorowscolumns       = MatZeroRowsColumns_HYPRE;
2433:   B->ops->zeroentries           = MatZeroEntries_HYPRE;
2434:   B->ops->zerorows              = MatZeroRows_HYPRE;
2435:   B->ops->getrow                = MatGetRow_HYPRE;
2436:   B->ops->restorerow            = MatRestoreRow_HYPRE;
2437:   B->ops->getvalues             = MatGetValues_HYPRE;
2438:   B->ops->setoption             = MatSetOption_HYPRE;
2439:   B->ops->duplicate             = MatDuplicate_HYPRE;
2440:   B->ops->copy                  = MatCopy_HYPRE;
2441:   B->ops->view                  = MatView_HYPRE;
2442:   B->ops->getdiagonal           = MatGetDiagonal_HYPRE;
2443:   B->ops->axpy                  = MatAXPY_HYPRE;
2444:   B->ops->productsetfromoptions = MatProductSetFromOptions_HYPRE;
2445: #if defined(PETSC_HAVE_HYPRE_DEVICE)
2446:   B->ops->bindtocpu = MatBindToCPU_HYPRE;
2447:   B->boundtocpu     = PETSC_FALSE;
2448: #endif

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

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