Actual source code: mhypre.c

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    318:     /* allocate CSR for local matrix */
    319:     PetscCall(PetscMalloc1(dr + 1, &iptr));
    320:     PetscCall(PetscMalloc1(nnz, &jptr));
    321:     PetscCall(PetscMalloc1(nnz, &data));
    322:   } else {
    323:     PetscInt  nr;
    324:     PetscBool done;
    325:     PetscCall(MatISGetLocalMat(*B, &lA));
    326:     PetscCall(MatGetRowIJ(lA, 0, PETSC_FALSE, PETSC_FALSE, &nr, (const PetscInt **)&iptr, (const PetscInt **)&jptr, &done));
    327:     PetscCheck(nr == dr, PETSC_COMM_SELF, PETSC_ERR_USER, "Cannot reuse mat: invalid number of rows in local mat! %" PetscInt_FMT " != %" PetscInt_FMT, nr, dr);
    328:     PetscCheck(iptr[nr] >= nnz, PETSC_COMM_SELF, PETSC_ERR_USER, "Cannot reuse mat: invalid number of nonzeros in local mat! reuse %" PetscInt_FMT " requested %" PetscInt_FMT, iptr[nr], nnz);
    329:     PetscCall(MatSeqAIJGetArrayWrite(lA, &data));
    330:   }
    331:   /* merge local matrices */
    332:   ii  = iptr;
    333:   jj  = jptr;
    334:   aa  = (HYPRE_Complex *)data; /* this cast fixes the clang error when doing the assignments below: implicit conversion from 'HYPRE_Complex' (aka '_Complex double') to 'double' is not permitted in C++ */
    335:   *ii = *(hdi++) + *(hoi++);
    336:   for (jd = 0, jo = 0, cum = 0; *ii < nnz; cum++) {
    337:     PetscScalar *aold = (PetscScalar *)aa;
    338:     PetscInt    *jold = jj, nc = jd + jo;
    339:     for (; jd < *hdi; jd++) {
    340:       *jj++ = *hdj++;
    341:       *aa++ = *hdd++;
    342:     }
    343:     for (; jo < *hoi; jo++) {
    344:       *jj++ = *hoj++ + dc;
    345:       *aa++ = *hod++;
    346:     }
    347:     *(++ii) = *(hdi++) + *(hoi++);
    348:     PetscCall(PetscSortIntWithScalarArray(jd + jo - nc, jold, aold));
    349:   }
    350:   for (; cum < dr; cum++) *(++ii) = nnz;
    351:   if (reuse != MAT_REUSE_MATRIX) {
    352:     Mat_SeqAIJ *a;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    661:     PetscCall(MatHYPRE_CreateCOOMat(M));

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    1468:   PetscCallExternal(HYPRE_IJMatrixAssemble, hA->ij);
    1469:   /* The assembly routine destroys the aux_matrix, we recreate it here by calling HYPRE_IJMatrixInitialize */
    1470:   /* If the option MAT_SORTED_FULL is set to true, the indices and values can be passed to hypre directly, so we don't need the aux_matrix */
    1471:   if (!A->sortedfull) {
    1472:     hypre_AuxParCSRMatrix *aux_matrix;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    1621:   PetscFunctionBegin;
    1622:   PetscCall(PetscLayoutSetUp(A->rmap));
    1623:   PetscCall(PetscLayoutSetUp(A->cmap));
    1624:   rs = A->rmap->rstart;
    1625:   re = A->rmap->rend;
    1626:   cs = A->cmap->rstart;
    1627:   ce = A->cmap->rend;
    1628:   if (!hA->ij) {
    1629:     PetscCallExternal(HYPRE_IJMatrixCreate, hA->comm, rs, re - 1, cs, ce - 1, &hA->ij);
    1630:     PetscCallExternal(HYPRE_IJMatrixSetObjectType, hA->ij, HYPRE_PARCSR);
    1631:   } else {
    1632:     HYPRE_BigInt hrs, hre, hcs, hce;
    1633:     PetscCallExternal(HYPRE_IJMatrixGetLocalRange, hA->ij, &hrs, &hre, &hcs, &hce);
    1634:     PetscCheck(hre - hrs + 1 == re - rs, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent local rows: IJMatrix [%" PetscHYPRE_BigInt_FMT ",%" PetscHYPRE_BigInt_FMT "), PETSc [%" PetscInt_FMT ",%" PetscInt_FMT ")", hrs, hre + 1, rs, re);
    1635:     PetscCheck(hce - hcs + 1 == ce - cs, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Inconsistent local cols: IJMatrix [%" PetscHYPRE_BigInt_FMT ",%" PetscHYPRE_BigInt_FMT "), PETSc [%" PetscInt_FMT ",%" PetscInt_FMT ")", hcs, hce + 1, cs, ce);
    1636:   }
    1637:   PetscCall(MatHYPRE_DestroyCOOMat(A));
    1638:   PetscCall(MatGetBlockSize(A, &bs));
    1639:   if (dnz == PETSC_DEFAULT || dnz == PETSC_DECIDE) dnz = 10 * bs;
    1640:   if (onz == PETSC_DEFAULT || onz == PETSC_DECIDE) onz = 10 * bs;

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

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

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

    1691:   Collective

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

    1711:   Level: intermediate

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

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

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

    1730:   Collective

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

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

    1740:   Level: intermediate

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

    1752:   PetscFunctionBegin;
    1753:   comm = hypre_ParCSRMatrixComm(parcsr);
    1754:   PetscCall(PetscStrcmp(mtype, MATSEQAIJ, &isseqaij));
    1755:   PetscCall(PetscStrcmp(mtype, MATSEQAIJMKL, &isseqaijmkl));
    1756:   PetscCall(PetscStrcmp(mtype, MATMPIAIJ, &ismpiaij));
    1757:   PetscCall(PetscStrcmp(mtype, MATAIJ, &isaij));
    1758:   PetscCall(PetscStrcmp(mtype, MATHYPRE, &ishyp));
    1759:   PetscCall(PetscStrcmp(mtype, MATIS, &isis));
    1760:   isaij = (PetscBool)(isseqaij || isseqaijmkl || ismpiaij || isaij);
    1761:   /* TODO */
    1762:   PetscCheck(isaij || ishyp || isis, comm, PETSC_ERR_SUP, "Unsupported MatType %s! Supported types are %s, %s, %s, %s, %s, and %s", mtype, MATAIJ, MATSEQAIJ, MATSEQAIJMKL, MATMPIAIJ, MATIS, MATHYPRE);
    1763:   /* access ParCSRMatrix */
    1764:   rstart = hypre_ParCSRMatrixFirstRowIndex(parcsr);
    1765:   rend   = hypre_ParCSRMatrixLastRowIndex(parcsr);
    1766:   cstart = hypre_ParCSRMatrixFirstColDiag(parcsr);
    1767:   cend   = hypre_ParCSRMatrixLastColDiag(parcsr);
    1768:   M      = hypre_ParCSRMatrixGlobalNumRows(parcsr);
    1769:   N      = hypre_ParCSRMatrixGlobalNumCols(parcsr);

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

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

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

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

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

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

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

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

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

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

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

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

    1871:   Not Collective, No Fortran Support

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

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

    1879:   Level: intermediate

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    2272: #include <petscblaslapack.h>

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    2432:    Level: intermediate

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

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

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

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

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

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

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