Actual source code: mhypre.c


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

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

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

 23: static PetscErrorCode MatHYPRE_CreateFromMat(Mat, Mat_HYPRE *);
 24: static PetscErrorCode MatHYPRE_IJMatrixPreallocate(Mat, Mat, HYPRE_IJMatrix);
 25: static PetscErrorCode MatHYPRE_IJMatrixFastCopy_MPIAIJ(Mat, HYPRE_IJMatrix);
 26: static PetscErrorCode MatHYPRE_IJMatrixFastCopy_SeqAIJ(Mat, HYPRE_IJMatrix);
 27: static PetscErrorCode MatHYPRE_MultKernel_Private(Mat, HYPRE_Complex, Vec, HYPRE_Complex, Vec, PetscBool);
 28: static PetscErrorCode hypre_array_destroy(void *);
 29: static PetscErrorCode MatSetValues_HYPRE(Mat, PetscInt, const PetscInt[], PetscInt, const PetscInt[], const PetscScalar[], InsertMode ins);

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

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

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

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

122: static PetscErrorCode MatHYPRE_IJMatrixCopy(Mat A, HYPRE_IJMatrix ij)
123: {
124:   PetscInt           i, rstart, rend, ncols, nr, nc;
125:   const PetscScalar *values;
126:   const PetscInt    *cols;
127:   PetscBool          flg;

129: #if PETSC_PKG_HYPRE_VERSION_LT(2, 19, 0)
130:   PetscCallExternal(HYPRE_IJMatrixInitialize, ij);
131: #else
132:   PetscCallExternal(HYPRE_IJMatrixInitialize_v2, ij, HYPRE_MEMORY_HOST);
133: #endif
134:   PetscObjectBaseTypeCompare((PetscObject)A, MATMPIAIJ, &flg);
135:   MatGetSize(A, &nr, &nc);
136:   if (flg && nr == nc) {
137:     MatHYPRE_IJMatrixFastCopy_MPIAIJ(A, ij);
138:     return 0;
139:   }
140:   PetscObjectBaseTypeCompare((PetscObject)A, MATSEQAIJ, &flg);
141:   if (flg) {
142:     MatHYPRE_IJMatrixFastCopy_SeqAIJ(A, ij);
143:     return 0;
144:   }

146:   /* Do not need Aux since we have done precise i[],j[] allocation in MatHYPRE_CreateFromMat() */
147:   hypre_AuxParCSRMatrixNeedAux((hypre_AuxParCSRMatrix *)hypre_IJMatrixTranslator(ij)) = 0;

149:   MatGetOwnershipRange(A, &rstart, &rend);
150:   for (i = rstart; i < rend; i++) {
151:     MatGetRow(A, i, &ncols, &cols, &values);
152:     if (ncols) {
153:       HYPRE_Int nc = (HYPRE_Int)ncols;

156:       PetscCallExternal(HYPRE_IJMatrixSetValues, ij, 1, &nc, (HYPRE_BigInt *)&i, (HYPRE_BigInt *)cols, (HYPRE_Complex *)values);
157:     }
158:     MatRestoreRow(A, i, &ncols, &cols, &values);
159:   }
160:   return 0;
161: }

163: static PetscErrorCode MatHYPRE_IJMatrixFastCopy_SeqAIJ(Mat A, HYPRE_IJMatrix ij)
164: {
165:   Mat_SeqAIJ            *pdiag = (Mat_SeqAIJ *)A->data;
166:   HYPRE_Int              type;
167:   hypre_ParCSRMatrix    *par_matrix;
168:   hypre_AuxParCSRMatrix *aux_matrix;
169:   hypre_CSRMatrix       *hdiag;
170:   PetscBool              sameint = (PetscBool)(sizeof(PetscInt) == sizeof(HYPRE_Int));
171:   const PetscScalar     *pa;

173:   PetscCallExternal(HYPRE_IJMatrixGetObjectType, ij, &type);
175:   PetscCallExternal(HYPRE_IJMatrixGetObject, ij, (void **)&par_matrix);
176:   hdiag = hypre_ParCSRMatrixDiag(par_matrix);
177:   /*
178:        this is the Hack part where we monkey directly with the hypre datastructures
179:   */
180:   if (sameint) {
181:     PetscArraycpy(hdiag->i, pdiag->i, A->rmap->n + 1);
182:     PetscArraycpy(hdiag->j, pdiag->j, pdiag->nz);
183:   } else {
184:     PetscInt i;

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

190:   MatSeqAIJGetArrayRead(A, &pa);
191:   PetscArraycpy(hdiag->data, pa, pdiag->nz);
192:   MatSeqAIJRestoreArrayRead(A, &pa);

194:   aux_matrix                               = (hypre_AuxParCSRMatrix *)hypre_IJMatrixTranslator(ij);
195:   hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 0;
196:   return 0;
197: }

199: static PetscErrorCode MatHYPRE_IJMatrixFastCopy_MPIAIJ(Mat A, HYPRE_IJMatrix ij)
200: {
201:   Mat_MPIAIJ            *pA = (Mat_MPIAIJ *)A->data;
202:   Mat_SeqAIJ            *pdiag, *poffd;
203:   PetscInt               i, *garray = pA->garray, *jj, cstart, *pjj;
204:   HYPRE_Int             *hjj, type;
205:   hypre_ParCSRMatrix    *par_matrix;
206:   hypre_AuxParCSRMatrix *aux_matrix;
207:   hypre_CSRMatrix       *hdiag, *hoffd;
208:   PetscBool              sameint = (PetscBool)(sizeof(PetscInt) == sizeof(HYPRE_Int));
209:   const PetscScalar     *pa;

211:   pdiag = (Mat_SeqAIJ *)pA->A->data;
212:   poffd = (Mat_SeqAIJ *)pA->B->data;
213:   /* cstart is only valid for square MPIAIJ layed out in the usual way */
214:   MatGetOwnershipRange(A, &cstart, NULL);

216:   PetscCallExternal(HYPRE_IJMatrixGetObjectType, ij, &type);
218:   PetscCallExternal(HYPRE_IJMatrixGetObject, ij, (void **)&par_matrix);
219:   hdiag = hypre_ParCSRMatrixDiag(par_matrix);
220:   hoffd = hypre_ParCSRMatrixOffd(par_matrix);

222:   /*
223:        this is the Hack part where we monkey directly with the hypre datastructures
224:   */
225:   if (sameint) {
226:     PetscArraycpy(hdiag->i, pdiag->i, pA->A->rmap->n + 1);
227:   } else {
228:     for (i = 0; i < pA->A->rmap->n + 1; i++) hdiag->i[i] = (HYPRE_Int)(pdiag->i[i]);
229:   }
230:   /* need to shift the diag column indices (hdiag->j) back to global numbering since hypre is expecting this */
231:   hjj = hdiag->j;
232:   pjj = pdiag->j;
233: #if PETSC_PKG_HYPRE_VERSION_GE(2, 16, 0)
234:   for (i = 0; i < pdiag->nz; i++) hjj[i] = pjj[i];
235: #else
236:   for (i = 0; i < pdiag->nz; i++) hjj[i] = cstart + pjj[i];
237: #endif
238:   MatSeqAIJGetArrayRead(pA->A, &pa);
239:   PetscArraycpy(hdiag->data, pa, pdiag->nz);
240:   MatSeqAIJRestoreArrayRead(pA->A, &pa);
241:   if (sameint) {
242:     PetscArraycpy(hoffd->i, poffd->i, pA->A->rmap->n + 1);
243:   } else {
244:     for (i = 0; i < pA->A->rmap->n + 1; i++) hoffd->i[i] = (HYPRE_Int)(poffd->i[i]);
245:   }

247:   /* need to move the offd column indices (hoffd->j) back to global numbering since hypre is expecting this
248:      If we hacked a hypre a bit more we might be able to avoid this step */
249: #if PETSC_PKG_HYPRE_VERSION_GE(2, 16, 0)
250:   PetscCallExternal(hypre_CSRMatrixBigInitialize, hoffd);
251:   jj = (PetscInt *)hoffd->big_j;
252: #else
253:   jj = (PetscInt *)hoffd->j;
254: #endif
255:   pjj = poffd->j;
256:   for (i = 0; i < poffd->nz; i++) jj[i] = garray[pjj[i]];

258:   MatSeqAIJGetArrayRead(pA->B, &pa);
259:   PetscArraycpy(hoffd->data, pa, poffd->nz);
260:   MatSeqAIJRestoreArrayRead(pA->B, &pa);

262:   aux_matrix                               = (hypre_AuxParCSRMatrix *)hypre_IJMatrixTranslator(ij);
263:   hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 0;
264:   return 0;
265: }

267: static PetscErrorCode MatConvert_HYPRE_IS(Mat A, MatType mtype, MatReuse reuse, Mat *B)
268: {
269:   Mat_HYPRE             *mhA = (Mat_HYPRE *)(A->data);
270:   Mat                    lA;
271:   ISLocalToGlobalMapping rl2g, cl2g;
272:   IS                     is;
273:   hypre_ParCSRMatrix    *hA;
274:   hypre_CSRMatrix       *hdiag, *hoffd;
275:   MPI_Comm               comm;
276:   HYPRE_Complex         *hdd, *hod, *aa;
277:   PetscScalar           *data;
278:   HYPRE_BigInt          *col_map_offd;
279:   HYPRE_Int             *hdi, *hdj, *hoi, *hoj;
280:   PetscInt              *ii, *jj, *iptr, *jptr;
281:   PetscInt               cum, dr, dc, oc, str, stc, nnz, i, jd, jo, M, N;
282:   HYPRE_Int              type;

284:   comm = PetscObjectComm((PetscObject)A);
285:   PetscCallExternal(HYPRE_IJMatrixGetObjectType, mhA->ij, &type);
287:   PetscCallExternal(HYPRE_IJMatrixGetObject, mhA->ij, (void **)&hA);
288:   M     = hypre_ParCSRMatrixGlobalNumRows(hA);
289:   N     = hypre_ParCSRMatrixGlobalNumCols(hA);
290:   str   = hypre_ParCSRMatrixFirstRowIndex(hA);
291:   stc   = hypre_ParCSRMatrixFirstColDiag(hA);
292:   hdiag = hypre_ParCSRMatrixDiag(hA);
293:   hoffd = hypre_ParCSRMatrixOffd(hA);
294:   dr    = hypre_CSRMatrixNumRows(hdiag);
295:   dc    = hypre_CSRMatrixNumCols(hdiag);
296:   nnz   = hypre_CSRMatrixNumNonzeros(hdiag);
297:   hdi   = hypre_CSRMatrixI(hdiag);
298:   hdj   = hypre_CSRMatrixJ(hdiag);
299:   hdd   = hypre_CSRMatrixData(hdiag);
300:   oc    = hypre_CSRMatrixNumCols(hoffd);
301:   nnz += hypre_CSRMatrixNumNonzeros(hoffd);
302:   hoi = hypre_CSRMatrixI(hoffd);
303:   hoj = hypre_CSRMatrixJ(hoffd);
304:   hod = hypre_CSRMatrixData(hoffd);
305:   if (reuse != MAT_REUSE_MATRIX) {
306:     PetscInt *aux;

308:     /* generate l2g maps for rows and cols */
309:     ISCreateStride(comm, dr, str, 1, &is);
310:     ISLocalToGlobalMappingCreateIS(is, &rl2g);
311:     ISDestroy(&is);
312:     col_map_offd = hypre_ParCSRMatrixColMapOffd(hA);
313:     PetscMalloc1(dc + oc, &aux);
314:     for (i = 0; i < dc; i++) aux[i] = i + stc;
315:     for (i = 0; i < oc; i++) aux[i + dc] = col_map_offd[i];
316:     ISCreateGeneral(comm, dc + oc, aux, PETSC_OWN_POINTER, &is);
317:     ISLocalToGlobalMappingCreateIS(is, &cl2g);
318:     ISDestroy(&is);
319:     /* create MATIS object */
320:     MatCreate(comm, B);
321:     MatSetSizes(*B, dr, dc, M, N);
322:     MatSetType(*B, MATIS);
323:     MatSetLocalToGlobalMapping(*B, rl2g, cl2g);
324:     ISLocalToGlobalMappingDestroy(&rl2g);
325:     ISLocalToGlobalMappingDestroy(&cl2g);

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

363:     MatCreateSeqAIJWithArrays(PETSC_COMM_SELF, dr, dc + oc, iptr, jptr, data, &lA);
364:     MatISSetLocalMat(*B, lA);
365:     /* hack SeqAIJ */
366:     a          = (Mat_SeqAIJ *)(lA->data);
367:     a->free_a  = PETSC_TRUE;
368:     a->free_ij = PETSC_TRUE;
369:     MatDestroy(&lA);
370:   }
371:   MatAssemblyBegin(*B, MAT_FINAL_ASSEMBLY);
372:   MatAssemblyEnd(*B, MAT_FINAL_ASSEMBLY);
373:   if (reuse == MAT_INPLACE_MATRIX) MatHeaderReplace(A, B);
374:   return 0;
375: }

377: PETSC_INTERN PetscErrorCode MatConvert_AIJ_HYPRE(Mat A, MatType type, MatReuse reuse, Mat *B)
378: {
379:   Mat        M = NULL;
380:   Mat_HYPRE *hB;
381:   MPI_Comm   comm = PetscObjectComm((PetscObject)A);

383:   if (reuse == MAT_REUSE_MATRIX) {
384:     /* always destroy the old matrix and create a new memory;
385:        hope this does not churn the memory too much. The problem
386:        is I do not know if it is possible to put the matrix back to
387:        its initial state so that we can directly copy the values
388:        the second time through. */
389:     hB = (Mat_HYPRE *)((*B)->data);
390:     PetscCallExternal(HYPRE_IJMatrixDestroy, hB->ij);
391:   } else {
392:     MatCreate(comm, &M);
393:     MatSetType(M, MATHYPRE);
394:     MatSetSizes(M, A->rmap->n, A->cmap->n, A->rmap->N, A->cmap->N);
395:     hB = (Mat_HYPRE *)(M->data);
396:     if (reuse == MAT_INITIAL_MATRIX) *B = M;
397:   }
398:   MatSetOption(*B, MAT_SORTED_FULL, PETSC_TRUE); /* "perfect" preallocation, so no need for hypre_AuxParCSRMatrixNeedAux */
399:   MatSetOption(*B, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE);
400:   MatHYPRE_CreateFromMat(A, hB);
401:   MatHYPRE_IJMatrixCopy(A, hB->ij);
402:   if (reuse == MAT_INPLACE_MATRIX) MatHeaderReplace(A, &M);
403:   (*B)->preallocated = PETSC_TRUE;
404:   MatAssemblyBegin(*B, MAT_FINAL_ASSEMBLY);
405:   MatAssemblyEnd(*B, MAT_FINAL_ASSEMBLY);
406:   return 0;
407: }

409: static PetscErrorCode MatConvert_HYPRE_AIJ(Mat A, MatType mtype, MatReuse reuse, Mat *B)
410: {
411:   Mat_HYPRE          *hA = (Mat_HYPRE *)A->data;
412:   hypre_ParCSRMatrix *parcsr;
413:   hypre_CSRMatrix    *hdiag, *hoffd;
414:   MPI_Comm            comm;
415:   PetscScalar        *da, *oa, *aptr;
416:   PetscInt           *dii, *djj, *oii, *ojj, *iptr;
417:   PetscInt            i, dnnz, onnz, m, n;
418:   HYPRE_Int           type;
419:   PetscMPIInt         size;
420:   PetscBool           sameint = (PetscBool)(sizeof(PetscInt) == sizeof(HYPRE_Int));

422:   comm = PetscObjectComm((PetscObject)A);
423:   PetscCallExternal(HYPRE_IJMatrixGetObjectType, hA->ij, &type);
425:   if (reuse == MAT_REUSE_MATRIX) {
426:     PetscBool ismpiaij, isseqaij;
427:     PetscObjectBaseTypeCompare((PetscObject)*B, MATMPIAIJ, &ismpiaij);
428:     PetscObjectBaseTypeCompare((PetscObject)*B, MATSEQAIJ, &isseqaij);
430:   }
431: #if defined(PETSC_HAVE_HYPRE_DEVICE)
433: #endif
434:   MPI_Comm_size(comm, &size);

436:   PetscCallExternal(HYPRE_IJMatrixGetObject, hA->ij, (void **)&parcsr);
437:   hdiag = hypre_ParCSRMatrixDiag(parcsr);
438:   hoffd = hypre_ParCSRMatrixOffd(parcsr);
439:   m     = hypre_CSRMatrixNumRows(hdiag);
440:   n     = hypre_CSRMatrixNumCols(hdiag);
441:   dnnz  = hypre_CSRMatrixNumNonzeros(hdiag);
442:   onnz  = hypre_CSRMatrixNumNonzeros(hoffd);
443:   if (reuse == MAT_INITIAL_MATRIX) {
444:     PetscMalloc1(m + 1, &dii);
445:     PetscMalloc1(dnnz, &djj);
446:     PetscMalloc1(dnnz, &da);
447:   } else if (reuse == MAT_REUSE_MATRIX) {
448:     PetscInt  nr;
449:     PetscBool done;
450:     if (size > 1) {
451:       Mat_MPIAIJ *b = (Mat_MPIAIJ *)((*B)->data);

453:       MatGetRowIJ(b->A, 0, PETSC_FALSE, PETSC_FALSE, &nr, (const PetscInt **)&dii, (const PetscInt **)&djj, &done);
456:       MatSeqAIJGetArray(b->A, &da);
457:     } else {
458:       MatGetRowIJ(*B, 0, PETSC_FALSE, PETSC_FALSE, &nr, (const PetscInt **)&dii, (const PetscInt **)&djj, &done);
461:       MatSeqAIJGetArray(*B, &da);
462:     }
463:   } else { /* MAT_INPLACE_MATRIX */
464:     if (!sameint) {
465:       PetscMalloc1(m + 1, &dii);
466:       PetscMalloc1(dnnz, &djj);
467:     } else {
468:       dii = (PetscInt *)hypre_CSRMatrixI(hdiag);
469:       djj = (PetscInt *)hypre_CSRMatrixJ(hdiag);
470:     }
471:     da = (PetscScalar *)hypre_CSRMatrixData(hdiag);
472:   }

474:   if (!sameint) {
475:     if (reuse != MAT_REUSE_MATRIX) {
476:       for (i = 0; i < m + 1; i++) dii[i] = (PetscInt)(hypre_CSRMatrixI(hdiag)[i]);
477:     }
478:     for (i = 0; i < dnnz; i++) djj[i] = (PetscInt)(hypre_CSRMatrixJ(hdiag)[i]);
479:   } else {
480:     if (reuse != MAT_REUSE_MATRIX) PetscArraycpy(dii, hypre_CSRMatrixI(hdiag), m + 1);
481:     PetscArraycpy(djj, hypre_CSRMatrixJ(hdiag), dnnz);
482:   }
483:   PetscArraycpy(da, hypre_CSRMatrixData(hdiag), dnnz);
484:   iptr = djj;
485:   aptr = da;
486:   for (i = 0; i < m; i++) {
487:     PetscInt nc = dii[i + 1] - dii[i];
488:     PetscSortIntWithScalarArray(nc, iptr, aptr);
489:     iptr += nc;
490:     aptr += nc;
491:   }
492:   if (size > 1) {
493:     HYPRE_BigInt *coffd;
494:     HYPRE_Int    *offdj;

496:     if (reuse == MAT_INITIAL_MATRIX) {
497:       PetscMalloc1(m + 1, &oii);
498:       PetscMalloc1(onnz, &ojj);
499:       PetscMalloc1(onnz, &oa);
500:     } else if (reuse == MAT_REUSE_MATRIX) {
501:       Mat_MPIAIJ *b = (Mat_MPIAIJ *)((*B)->data);
502:       PetscInt    nr, hr = hypre_CSRMatrixNumRows(hoffd);
503:       PetscBool   done;

505:       MatGetRowIJ(b->B, 0, PETSC_FALSE, PETSC_FALSE, &nr, (const PetscInt **)&oii, (const PetscInt **)&ojj, &done);
508:       MatSeqAIJGetArray(b->B, &oa);
509:     } else { /* MAT_INPLACE_MATRIX */
510:       if (!sameint) {
511:         PetscMalloc1(m + 1, &oii);
512:         PetscMalloc1(onnz, &ojj);
513:       } else {
514:         oii = (PetscInt *)hypre_CSRMatrixI(hoffd);
515:         ojj = (PetscInt *)hypre_CSRMatrixJ(hoffd);
516:       }
517:       oa = (PetscScalar *)hypre_CSRMatrixData(hoffd);
518:     }
519:     if (reuse != MAT_REUSE_MATRIX) {
520:       if (!sameint) {
521:         for (i = 0; i < m + 1; i++) oii[i] = (PetscInt)(hypre_CSRMatrixI(hoffd)[i]);
522:       } else {
523:         PetscArraycpy(oii, hypre_CSRMatrixI(hoffd), m + 1);
524:       }
525:     }
526:     PetscArraycpy(oa, hypre_CSRMatrixData(hoffd), onnz);

528:     offdj = hypre_CSRMatrixJ(hoffd);
529:     coffd = hypre_ParCSRMatrixColMapOffd(parcsr);
530:     /* we only need the permutation to be computed properly, I don't know if HYPRE
531:        messes up with the ordering. Just in case, allocate some memory and free it
532:        later */
533:     if (reuse == MAT_REUSE_MATRIX) {
534:       Mat_MPIAIJ *b = (Mat_MPIAIJ *)((*B)->data);
535:       PetscInt    mnz;

537:       MatSeqAIJGetMaxRowNonzeros(b->B, &mnz);
538:       PetscMalloc1(mnz, &ojj);
539:     } else
540:       for (i = 0; i < onnz; i++) ojj[i] = coffd[offdj[i]];
541:     iptr = ojj;
542:     aptr = oa;
543:     for (i = 0; i < m; i++) {
544:       PetscInt nc = oii[i + 1] - oii[i];
545:       if (reuse == MAT_REUSE_MATRIX) {
546:         PetscInt j;

548:         iptr = ojj;
549:         for (j = 0; j < nc; j++) iptr[j] = coffd[offdj[oii[i] + j]];
550:       }
551:       PetscSortIntWithScalarArray(nc, iptr, aptr);
552:       iptr += nc;
553:       aptr += nc;
554:     }
555:     if (reuse == MAT_REUSE_MATRIX) PetscFree(ojj);
556:     if (reuse == MAT_INITIAL_MATRIX) {
557:       Mat_MPIAIJ *b;
558:       Mat_SeqAIJ *d, *o;

560:       MatCreateMPIAIJWithSplitArrays(comm, m, n, PETSC_DECIDE, PETSC_DECIDE, dii, djj, da, oii, ojj, oa, B);
561:       /* hack MPIAIJ */
562:       b          = (Mat_MPIAIJ *)((*B)->data);
563:       d          = (Mat_SeqAIJ *)b->A->data;
564:       o          = (Mat_SeqAIJ *)b->B->data;
565:       d->free_a  = PETSC_TRUE;
566:       d->free_ij = PETSC_TRUE;
567:       o->free_a  = PETSC_TRUE;
568:       o->free_ij = PETSC_TRUE;
569:     } else if (reuse == MAT_INPLACE_MATRIX) {
570:       Mat T;

572:       MatCreateMPIAIJWithSplitArrays(comm, m, n, PETSC_DECIDE, PETSC_DECIDE, dii, djj, da, oii, ojj, oa, &T);
573:       if (sameint) { /* ownership of CSR pointers is transferred to PETSc */
574:         hypre_CSRMatrixI(hdiag) = NULL;
575:         hypre_CSRMatrixJ(hdiag) = NULL;
576:         hypre_CSRMatrixI(hoffd) = NULL;
577:         hypre_CSRMatrixJ(hoffd) = NULL;
578:       } else { /* Hack MPIAIJ -> free ij but not a */
579:         Mat_MPIAIJ *b = (Mat_MPIAIJ *)(T->data);
580:         Mat_SeqAIJ *d = (Mat_SeqAIJ *)(b->A->data);
581:         Mat_SeqAIJ *o = (Mat_SeqAIJ *)(b->B->data);

583:         d->free_ij = PETSC_TRUE;
584:         o->free_ij = PETSC_TRUE;
585:       }
586:       hypre_CSRMatrixData(hdiag) = NULL;
587:       hypre_CSRMatrixData(hoffd) = NULL;
588:       MatHeaderReplace(A, &T);
589:     }
590:   } else {
591:     oii = NULL;
592:     ojj = NULL;
593:     oa  = NULL;
594:     if (reuse == MAT_INITIAL_MATRIX) {
595:       Mat_SeqAIJ *b;

597:       MatCreateSeqAIJWithArrays(comm, m, n, dii, djj, da, B);
598:       /* hack SeqAIJ */
599:       b          = (Mat_SeqAIJ *)((*B)->data);
600:       b->free_a  = PETSC_TRUE;
601:       b->free_ij = PETSC_TRUE;
602:     } else if (reuse == MAT_INPLACE_MATRIX) {
603:       Mat T;

605:       MatCreateSeqAIJWithArrays(comm, m, n, dii, djj, da, &T);
606:       if (sameint) { /* ownership of CSR pointers is transferred to PETSc */
607:         hypre_CSRMatrixI(hdiag) = NULL;
608:         hypre_CSRMatrixJ(hdiag) = NULL;
609:       } else { /* free ij but not a */
610:         Mat_SeqAIJ *b = (Mat_SeqAIJ *)(T->data);

612:         b->free_ij = PETSC_TRUE;
613:       }
614:       hypre_CSRMatrixData(hdiag) = NULL;
615:       MatHeaderReplace(A, &T);
616:     }
617:   }

619:   /* we have to use hypre_Tfree to free the HYPRE arrays
620:      that PETSc now onws */
621:   if (reuse == MAT_INPLACE_MATRIX) {
622:     PetscInt    nh;
623:     void       *ptrs[6]  = {da, oa, dii, djj, oii, ojj};
624:     const char *names[6] = {"_hypre_csr_da", "_hypre_csr_oa", "_hypre_csr_dii", "_hypre_csr_djj", "_hypre_csr_oii", "_hypre_csr_ojj"};
625:     nh                   = sameint ? 6 : 2;
626:     for (i = 0; i < nh; i++) {
627:       PetscContainer c;

629:       PetscContainerCreate(comm, &c);
630:       PetscContainerSetPointer(c, ptrs[i]);
631:       PetscContainerSetUserDestroy(c, hypre_array_destroy);
632:       PetscObjectCompose((PetscObject)(*B), names[i], (PetscObject)c);
633:       PetscContainerDestroy(&c);
634:     }
635:   }
636:   return 0;
637: }

639: static PetscErrorCode MatAIJGetParCSR_Private(Mat A, hypre_ParCSRMatrix **hA)
640: {
641:   hypre_ParCSRMatrix *tA;
642:   hypre_CSRMatrix    *hdiag, *hoffd;
643:   Mat_SeqAIJ         *diag, *offd;
644:   PetscInt           *garray, i, noffd, dnnz, onnz, *row_starts, *col_starts;
645:   MPI_Comm            comm = PetscObjectComm((PetscObject)A);
646:   PetscBool           ismpiaij, isseqaij;
647:   PetscBool           sameint = (PetscBool)(sizeof(PetscInt) == sizeof(HYPRE_Int));
648:   HYPRE_Int          *hdi = NULL, *hdj = NULL, *hoi = NULL, *hoj = NULL;
649:   PetscInt           *pdi = NULL, *pdj = NULL, *poi = NULL, *poj = NULL;
650: #if defined(PETSC_HAVE_HYPRE_DEVICE)
651:   PetscBool iscuda = PETSC_FALSE;
652: #endif

654:   PetscObjectBaseTypeCompare((PetscObject)A, MATMPIAIJ, &ismpiaij);
655:   PetscObjectBaseTypeCompare((PetscObject)A, MATSEQAIJ, &isseqaij);
657:   if (ismpiaij) {
658:     Mat_MPIAIJ *a = (Mat_MPIAIJ *)(A->data);

660:     diag = (Mat_SeqAIJ *)a->A->data;
661:     offd = (Mat_SeqAIJ *)a->B->data;
662: #if defined(PETSC_HAVE_CUDA) && defined(PETSC_HAVE_HYPRE_DEVICE) && defined(HYPRE_USING_CUDA)
663:     PetscObjectTypeCompare((PetscObject)A, MATMPIAIJCUSPARSE, &iscuda);
664:     if (iscuda && !A->boundtocpu) {
665:       sameint = PETSC_TRUE;
666:       MatSeqAIJCUSPARSEGetIJ(a->A, PETSC_FALSE, (const HYPRE_Int **)&hdi, (const HYPRE_Int **)&hdj);
667:       MatSeqAIJCUSPARSEGetIJ(a->B, PETSC_FALSE, (const HYPRE_Int **)&hoi, (const HYPRE_Int **)&hoj);
668:     } else {
669: #else
670:     {
671: #endif
672:       pdi = diag->i;
673:       pdj = diag->j;
674:       poi = offd->i;
675:       poj = offd->j;
676:       if (sameint) {
677:         hdi = (HYPRE_Int *)pdi;
678:         hdj = (HYPRE_Int *)pdj;
679:         hoi = (HYPRE_Int *)poi;
680:         hoj = (HYPRE_Int *)poj;
681:       }
682:     }
683:     garray = a->garray;
684:     noffd  = a->B->cmap->N;
685:     dnnz   = diag->nz;
686:     onnz   = offd->nz;
687:   } else {
688:     diag = (Mat_SeqAIJ *)A->data;
689:     offd = NULL;
690: #if defined(PETSC_HAVE_CUDA) && defined(PETSC_HAVE_HYPRE_DEVICE)
691:     PetscObjectTypeCompare((PetscObject)A, MATSEQAIJCUSPARSE, &iscuda);
692:     if (iscuda && !A->boundtocpu) {
693:       sameint = PETSC_TRUE;
694:       MatSeqAIJCUSPARSEGetIJ(A, PETSC_FALSE, (const HYPRE_Int **)&hdi, (const HYPRE_Int **)&hdj);
695:     } else {
696: #else
697:     {
698: #endif
699:       pdi = diag->i;
700:       pdj = diag->j;
701:       if (sameint) {
702:         hdi = (HYPRE_Int *)pdi;
703:         hdj = (HYPRE_Int *)pdj;
704:       }
705:     }
706:     garray = NULL;
707:     noffd  = 0;
708:     dnnz   = diag->nz;
709:     onnz   = 0;
710:   }

712:   /* create a temporary ParCSR */
713:   if (HYPRE_AssumedPartitionCheck()) {
714:     PetscMPIInt myid;

716:     MPI_Comm_rank(comm, &myid);
717:     row_starts = A->rmap->range + myid;
718:     col_starts = A->cmap->range + myid;
719:   } else {
720:     row_starts = A->rmap->range;
721:     col_starts = A->cmap->range;
722:   }
723:   tA = hypre_ParCSRMatrixCreate(comm, A->rmap->N, A->cmap->N, (HYPRE_BigInt *)row_starts, (HYPRE_BigInt *)col_starts, noffd, dnnz, onnz);
724: #if defined(hypre_ParCSRMatrixOwnsRowStarts)
725:   hypre_ParCSRMatrixSetRowStartsOwner(tA, 0);
726:   hypre_ParCSRMatrixSetColStartsOwner(tA, 0);
727: #endif

729:   /* set diagonal part */
730:   hdiag = hypre_ParCSRMatrixDiag(tA);
731:   if (!sameint) { /* malloc CSR pointers */
732:     PetscMalloc2(A->rmap->n + 1, &hdi, dnnz, &hdj);
733:     for (i = 0; i < A->rmap->n + 1; i++) hdi[i] = (HYPRE_Int)(pdi[i]);
734:     for (i = 0; i < dnnz; i++) hdj[i] = (HYPRE_Int)(pdj[i]);
735:   }
736:   hypre_CSRMatrixI(hdiag)           = hdi;
737:   hypre_CSRMatrixJ(hdiag)           = hdj;
738:   hypre_CSRMatrixData(hdiag)        = (HYPRE_Complex *)diag->a;
739:   hypre_CSRMatrixNumNonzeros(hdiag) = diag->nz;
740:   hypre_CSRMatrixSetRownnz(hdiag);
741:   hypre_CSRMatrixSetDataOwner(hdiag, 0);

743:   /* set offdiagonal part */
744:   hoffd = hypre_ParCSRMatrixOffd(tA);
745:   if (offd) {
746:     if (!sameint) { /* malloc CSR pointers */
747:       PetscMalloc2(A->rmap->n + 1, &hoi, onnz, &hoj);
748:       for (i = 0; i < A->rmap->n + 1; i++) hoi[i] = (HYPRE_Int)(poi[i]);
749:       for (i = 0; i < onnz; i++) hoj[i] = (HYPRE_Int)(poj[i]);
750:     }
751:     hypre_CSRMatrixI(hoffd)           = hoi;
752:     hypre_CSRMatrixJ(hoffd)           = hoj;
753:     hypre_CSRMatrixData(hoffd)        = (HYPRE_Complex *)offd->a;
754:     hypre_CSRMatrixNumNonzeros(hoffd) = offd->nz;
755:     hypre_CSRMatrixSetRownnz(hoffd);
756:     hypre_CSRMatrixSetDataOwner(hoffd, 0);
757:   }
758: #if defined(PETSC_HAVE_HYPRE_DEVICE)
759:   PetscCallExternal(hypre_ParCSRMatrixInitialize_v2, tA, iscuda ? HYPRE_MEMORY_DEVICE : HYPRE_MEMORY_HOST);
760: #else
761:   #if PETSC_PKG_HYPRE_VERSION_LT(2, 18, 0)
762:   PetscCallExternal(hypre_ParCSRMatrixInitialize, tA);
763:   #else
764:   PetscCallExternal(hypre_ParCSRMatrixInitialize_v2, tA, HYPRE_MEMORY_HOST);
765:   #endif
766: #endif
767:   hypre_TFree(hypre_ParCSRMatrixColMapOffd(tA), HYPRE_MEMORY_HOST);
768:   hypre_ParCSRMatrixSetNumNonzeros(tA);
769:   hypre_ParCSRMatrixColMapOffd(tA) = (HYPRE_BigInt *)garray;
770:   if (!hypre_ParCSRMatrixCommPkg(tA)) PetscCallExternal(hypre_MatvecCommPkgCreate, tA);
771:   *hA = tA;
772:   return 0;
773: }

775: static PetscErrorCode MatAIJRestoreParCSR_Private(Mat A, hypre_ParCSRMatrix **hA)
776: {
777:   hypre_CSRMatrix *hdiag, *hoffd;
778:   PetscBool        ismpiaij, sameint = (PetscBool)(sizeof(PetscInt) == sizeof(HYPRE_Int));
779: #if defined(PETSC_HAVE_HYPRE_DEVICE)
780:   PetscBool iscuda = PETSC_FALSE;
781: #endif

783:   PetscObjectBaseTypeCompare((PetscObject)A, MATMPIAIJ, &ismpiaij);
784: #if defined(PETSC_HAVE_HYPRE_DEVICE)
785:   PetscObjectTypeCompareAny((PetscObject)A, &iscuda, MATSEQAIJCUSPARSE, MATMPIAIJCUSPARSE, "");
786:   if (iscuda) sameint = PETSC_TRUE;
787: #endif
788:   hdiag = hypre_ParCSRMatrixDiag(*hA);
789:   hoffd = hypre_ParCSRMatrixOffd(*hA);
790:   /* free temporary memory allocated by PETSc
791:      set pointers to NULL before destroying tA */
792:   if (!sameint) {
793:     HYPRE_Int *hi, *hj;

795:     hi = hypre_CSRMatrixI(hdiag);
796:     hj = hypre_CSRMatrixJ(hdiag);
797:     PetscFree2(hi, hj);
798:     if (ismpiaij) {
799:       hi = hypre_CSRMatrixI(hoffd);
800:       hj = hypre_CSRMatrixJ(hoffd);
801:       PetscFree2(hi, hj);
802:     }
803:   }
804:   hypre_CSRMatrixI(hdiag)    = NULL;
805:   hypre_CSRMatrixJ(hdiag)    = NULL;
806:   hypre_CSRMatrixData(hdiag) = NULL;
807:   if (ismpiaij) {
808:     hypre_CSRMatrixI(hoffd)    = NULL;
809:     hypre_CSRMatrixJ(hoffd)    = NULL;
810:     hypre_CSRMatrixData(hoffd) = NULL;
811:   }
812:   hypre_ParCSRMatrixColMapOffd(*hA) = NULL;
813:   hypre_ParCSRMatrixDestroy(*hA);
814:   *hA = NULL;
815:   return 0;
816: }

818: /* calls RAP from BoomerAMG:
819:    the resulting ParCSR will not own the column and row starts
820:    It looks like we don't need to have the diagonal entries ordered first */
821: static PetscErrorCode MatHYPRE_ParCSR_RAP(hypre_ParCSRMatrix *hR, hypre_ParCSRMatrix *hA, hypre_ParCSRMatrix *hP, hypre_ParCSRMatrix **hRAP)
822: {
823: #if defined(hypre_ParCSRMatrixOwnsRowStarts)
824:   HYPRE_Int P_owns_col_starts, R_owns_row_starts;
825: #endif

827: #if defined(hypre_ParCSRMatrixOwnsRowStarts)
828:   P_owns_col_starts = hypre_ParCSRMatrixOwnsColStarts(hP);
829:   R_owns_row_starts = hypre_ParCSRMatrixOwnsRowStarts(hR);
830: #endif
831:   /* can be replaced by version test later */
832: #if defined(PETSC_HAVE_HYPRE_DEVICE)
833:   PetscStackPushExternal("hypre_ParCSRMatrixRAP");
834:   *hRAP = hypre_ParCSRMatrixRAP(hR, hA, hP);
835:   PetscStackPop;
836: #else
837:   PetscCallExternal(hypre_BoomerAMGBuildCoarseOperator, hR, hA, hP, hRAP);
838:   PetscCallExternal(hypre_ParCSRMatrixSetNumNonzeros, *hRAP);
839: #endif
840:   /* hypre_BoomerAMGBuildCoarseOperator steals the col_starts from P and the row_starts from R */
841: #if defined(hypre_ParCSRMatrixOwnsRowStarts)
842:   hypre_ParCSRMatrixSetRowStartsOwner(*hRAP, 0);
843:   hypre_ParCSRMatrixSetColStartsOwner(*hRAP, 0);
844:   if (P_owns_col_starts) hypre_ParCSRMatrixSetColStartsOwner(hP, 1);
845:   if (R_owns_row_starts) hypre_ParCSRMatrixSetRowStartsOwner(hR, 1);
846: #endif
847:   return 0;
848: }

850: static PetscErrorCode MatPtAPNumeric_AIJ_AIJ_wHYPRE(Mat A, Mat P, Mat C)
851: {
852:   Mat                 B;
853:   hypre_ParCSRMatrix *hA, *hP, *hPtAP = NULL;
854:   Mat_Product        *product = C->product;

856:   MatAIJGetParCSR_Private(A, &hA);
857:   MatAIJGetParCSR_Private(P, &hP);
858:   MatHYPRE_ParCSR_RAP(hP, hA, hP, &hPtAP);
859:   MatCreateFromParCSR(hPtAP, MATAIJ, PETSC_OWN_POINTER, &B);

861:   MatHeaderMerge(C, &B);
862:   C->product = product;

864:   MatAIJRestoreParCSR_Private(A, &hA);
865:   MatAIJRestoreParCSR_Private(P, &hP);
866:   return 0;
867: }

869: PETSC_INTERN PetscErrorCode MatPtAPSymbolic_AIJ_AIJ_wHYPRE(Mat A, Mat P, PetscReal fill, Mat C)
870: {
871:   MatSetType(C, MATAIJ);
872:   C->ops->ptapnumeric    = MatPtAPNumeric_AIJ_AIJ_wHYPRE;
873:   C->ops->productnumeric = MatProductNumeric_PtAP;
874:   return 0;
875: }

877: static PetscErrorCode MatPtAPNumeric_AIJ_HYPRE(Mat A, Mat P, Mat C)
878: {
879:   Mat                 B;
880:   Mat_HYPRE          *hP;
881:   hypre_ParCSRMatrix *hA = NULL, *Pparcsr, *ptapparcsr = NULL;
882:   HYPRE_Int           type;
883:   MPI_Comm            comm = PetscObjectComm((PetscObject)A);
884:   PetscBool           ishypre;

886:   PetscObjectTypeCompare((PetscObject)P, MATHYPRE, &ishypre);
888:   hP = (Mat_HYPRE *)P->data;
889:   PetscCallExternal(HYPRE_IJMatrixGetObjectType, hP->ij, &type);
891:   PetscCallExternal(HYPRE_IJMatrixGetObject, hP->ij, (void **)&Pparcsr);

893:   MatAIJGetParCSR_Private(A, &hA);
894:   MatHYPRE_ParCSR_RAP(Pparcsr, hA, Pparcsr, &ptapparcsr);
895:   MatAIJRestoreParCSR_Private(A, &hA);

897:   /* create temporary matrix and merge to C */
898:   MatCreateFromParCSR(ptapparcsr, ((PetscObject)C)->type_name, PETSC_OWN_POINTER, &B);
899:   MatHeaderMerge(C, &B);
900:   return 0;
901: }

903: static PetscErrorCode MatPtAPNumeric_HYPRE_HYPRE(Mat A, Mat P, Mat C)
904: {
905:   Mat                 B;
906:   hypre_ParCSRMatrix *Aparcsr, *Pparcsr, *ptapparcsr = NULL;
907:   Mat_HYPRE          *hA, *hP;
908:   PetscBool           ishypre;
909:   HYPRE_Int           type;

911:   PetscObjectTypeCompare((PetscObject)P, MATHYPRE, &ishypre);
913:   PetscObjectTypeCompare((PetscObject)A, MATHYPRE, &ishypre);
915:   hA = (Mat_HYPRE *)A->data;
916:   hP = (Mat_HYPRE *)P->data;
917:   PetscCallExternal(HYPRE_IJMatrixGetObjectType, hA->ij, &type);
919:   PetscCallExternal(HYPRE_IJMatrixGetObjectType, hP->ij, &type);
921:   PetscCallExternal(HYPRE_IJMatrixGetObject, hA->ij, (void **)&Aparcsr);
922:   PetscCallExternal(HYPRE_IJMatrixGetObject, hP->ij, (void **)&Pparcsr);
923:   MatHYPRE_ParCSR_RAP(Pparcsr, Aparcsr, Pparcsr, &ptapparcsr);
924:   MatCreateFromParCSR(ptapparcsr, MATHYPRE, PETSC_OWN_POINTER, &B);
925:   MatHeaderMerge(C, &B);
926:   return 0;
927: }

929: /* calls hypre_ParMatmul
930:    hypre_ParMatMul uses hypre_ParMatrixCreate with the communicator of hA
931:    hypre_ParMatrixCreate does not duplicate the communicator
932:    It looks like we don't need to have the diagonal entries ordered first */
933: static PetscErrorCode MatHYPRE_ParCSR_MatMatMult(hypre_ParCSRMatrix *hA, hypre_ParCSRMatrix *hB, hypre_ParCSRMatrix **hAB)
934: {
935:   /* can be replaced by version test later */
936: #if defined(PETSC_HAVE_HYPRE_DEVICE)
937:   PetscStackPushExternal("hypre_ParCSRMatMat");
938:   *hAB = hypre_ParCSRMatMat(hA, hB);
939: #else
940:   PetscStackPushExternal("hypre_ParMatmul");
941:   *hAB = hypre_ParMatmul(hA, hB);
942: #endif
943:   PetscStackPop;
944:   return 0;
945: }

947: static PetscErrorCode MatMatMultNumeric_AIJ_AIJ_wHYPRE(Mat A, Mat B, Mat C)
948: {
949:   Mat                 D;
950:   hypre_ParCSRMatrix *hA, *hB, *hAB = NULL;
951:   Mat_Product        *product = C->product;

953:   MatAIJGetParCSR_Private(A, &hA);
954:   MatAIJGetParCSR_Private(B, &hB);
955:   MatHYPRE_ParCSR_MatMatMult(hA, hB, &hAB);
956:   MatCreateFromParCSR(hAB, MATAIJ, PETSC_OWN_POINTER, &D);

958:   MatHeaderMerge(C, &D);
959:   C->product = product;

961:   MatAIJRestoreParCSR_Private(A, &hA);
962:   MatAIJRestoreParCSR_Private(B, &hB);
963:   return 0;
964: }

966: PETSC_INTERN PetscErrorCode MatMatMultSymbolic_AIJ_AIJ_wHYPRE(Mat A, Mat B, PetscReal fill, Mat C)
967: {
968:   MatSetType(C, MATAIJ);
969:   C->ops->matmultnumeric = MatMatMultNumeric_AIJ_AIJ_wHYPRE;
970:   C->ops->productnumeric = MatProductNumeric_AB;
971:   return 0;
972: }

974: static PetscErrorCode MatMatMultNumeric_HYPRE_HYPRE(Mat A, Mat B, Mat C)
975: {
976:   Mat                 D;
977:   hypre_ParCSRMatrix *Aparcsr, *Bparcsr, *ABparcsr = NULL;
978:   Mat_HYPRE          *hA, *hB;
979:   PetscBool           ishypre;
980:   HYPRE_Int           type;
981:   Mat_Product        *product;

983:   PetscObjectTypeCompare((PetscObject)B, MATHYPRE, &ishypre);
985:   PetscObjectTypeCompare((PetscObject)A, MATHYPRE, &ishypre);
987:   hA = (Mat_HYPRE *)A->data;
988:   hB = (Mat_HYPRE *)B->data;
989:   PetscCallExternal(HYPRE_IJMatrixGetObjectType, hA->ij, &type);
991:   PetscCallExternal(HYPRE_IJMatrixGetObjectType, hB->ij, &type);
993:   PetscCallExternal(HYPRE_IJMatrixGetObject, hA->ij, (void **)&Aparcsr);
994:   PetscCallExternal(HYPRE_IJMatrixGetObject, hB->ij, (void **)&Bparcsr);
995:   MatHYPRE_ParCSR_MatMatMult(Aparcsr, Bparcsr, &ABparcsr);
996:   MatCreateFromParCSR(ABparcsr, MATHYPRE, PETSC_OWN_POINTER, &D);

998:   /* need to use HeaderReplace because HeaderMerge messes up with the communicator */
999:   product    = C->product; /* save it from MatHeaderReplace() */
1000:   C->product = NULL;
1001:   MatHeaderReplace(C, &D);
1002:   C->product             = product;
1003:   C->ops->matmultnumeric = MatMatMultNumeric_HYPRE_HYPRE;
1004:   C->ops->productnumeric = MatProductNumeric_AB;
1005:   return 0;
1006: }

1008: PETSC_INTERN PetscErrorCode MatTransposeMatMatMultNumeric_AIJ_AIJ_AIJ_wHYPRE(Mat A, Mat B, Mat C, Mat D)
1009: {
1010:   Mat                 E;
1011:   hypre_ParCSRMatrix *hA, *hB, *hC, *hABC = NULL;

1013:   MatAIJGetParCSR_Private(A, &hA);
1014:   MatAIJGetParCSR_Private(B, &hB);
1015:   MatAIJGetParCSR_Private(C, &hC);
1016:   MatHYPRE_ParCSR_RAP(hA, hB, hC, &hABC);
1017:   MatCreateFromParCSR(hABC, MATAIJ, PETSC_OWN_POINTER, &E);
1018:   MatHeaderMerge(D, &E);
1019:   MatAIJRestoreParCSR_Private(A, &hA);
1020:   MatAIJRestoreParCSR_Private(B, &hB);
1021:   MatAIJRestoreParCSR_Private(C, &hC);
1022:   return 0;
1023: }

1025: PETSC_INTERN PetscErrorCode MatTransposeMatMatMultSymbolic_AIJ_AIJ_AIJ_wHYPRE(Mat A, Mat B, Mat C, PetscReal fill, Mat D)
1026: {
1027:   MatSetType(D, MATAIJ);
1028:   return 0;
1029: }

1031: /* ---------------------------------------------------- */
1032: static PetscErrorCode MatProductSymbolic_AB_HYPRE(Mat C)
1033: {
1034:   C->ops->productnumeric = MatProductNumeric_AB;
1035:   return 0;
1036: }

1038: static PetscErrorCode MatProductSetFromOptions_HYPRE_AB(Mat C)
1039: {
1040:   Mat_Product *product = C->product;
1041:   PetscBool    Ahypre;

1043:   PetscObjectTypeCompare((PetscObject)product->A, MATHYPRE, &Ahypre);
1044:   if (Ahypre) { /* A is a Hypre matrix */
1045:     MatSetType(C, MATHYPRE);
1046:     C->ops->productsymbolic = MatProductSymbolic_AB_HYPRE;
1047:     C->ops->matmultnumeric  = MatMatMultNumeric_HYPRE_HYPRE;
1048:     return 0;
1049:   }
1050:   return 0;
1051: }

1053: static PetscErrorCode MatProductSymbolic_PtAP_HYPRE(Mat C)
1054: {
1055:   C->ops->productnumeric = MatProductNumeric_PtAP;
1056:   return 0;
1057: }

1059: static PetscErrorCode MatProductSetFromOptions_HYPRE_PtAP(Mat C)
1060: {
1061:   Mat_Product *product = C->product;
1062:   PetscBool    flg;
1063:   PetscInt     type        = 0;
1064:   const char  *outTypes[4] = {"aij", "seqaij", "mpiaij", "hypre"};
1065:   PetscInt     ntype       = 4;
1066:   Mat          A           = product->A;
1067:   PetscBool    Ahypre;

1069:   PetscObjectTypeCompare((PetscObject)A, MATHYPRE, &Ahypre);
1070:   if (Ahypre) { /* A is a Hypre matrix */
1071:     MatSetType(C, MATHYPRE);
1072:     C->ops->productsymbolic = MatProductSymbolic_PtAP_HYPRE;
1073:     C->ops->ptapnumeric     = MatPtAPNumeric_HYPRE_HYPRE;
1074:     return 0;
1075:   }

1077:   /* A is AIJ, P is Hypre, C = PtAP can be either AIJ or Hypre format */
1078:   /* Get runtime option */
1079:   if (product->api_user) {
1080:     PetscOptionsBegin(PetscObjectComm((PetscObject)C), ((PetscObject)C)->prefix, "MatPtAP_HYPRE", "Mat");
1081:     PetscOptionsEList("-matptap_hypre_outtype", "MatPtAP outtype", "MatPtAP outtype", outTypes, ntype, outTypes[type], &type, &flg);
1082:     PetscOptionsEnd();
1083:   } else {
1084:     PetscOptionsBegin(PetscObjectComm((PetscObject)C), ((PetscObject)C)->prefix, "MatProduct_PtAP_HYPRE", "Mat");
1085:     PetscOptionsEList("-mat_product_algorithm_hypre_outtype", "MatProduct_PtAP outtype", "MatProduct_PtAP", outTypes, ntype, outTypes[type], &type, &flg);
1086:     PetscOptionsEnd();
1087:   }

1089:   if (type == 0 || type == 1 || type == 2) {
1090:     MatSetType(C, MATAIJ);
1091:   } else if (type == 3) {
1092:     MatSetType(C, MATHYPRE);
1093:   } else SETERRQ(PetscObjectComm((PetscObject)C), PETSC_ERR_SUP, "MatPtAP outtype is not supported");
1094:   C->ops->productsymbolic = MatProductSymbolic_PtAP_HYPRE;
1095:   C->ops->ptapnumeric     = MatPtAPNumeric_AIJ_HYPRE;
1096:   return 0;
1097: }

1099: static PetscErrorCode MatProductSetFromOptions_HYPRE(Mat C)
1100: {
1101:   Mat_Product *product = C->product;

1103:   switch (product->type) {
1104:   case MATPRODUCT_AB:
1105:     MatProductSetFromOptions_HYPRE_AB(C);
1106:     break;
1107:   case MATPRODUCT_PtAP:
1108:     MatProductSetFromOptions_HYPRE_PtAP(C);
1109:     break;
1110:   default:
1111:     break;
1112:   }
1113:   return 0;
1114: }

1116: /* -------------------------------------------------------- */

1118: static PetscErrorCode MatMultTranspose_HYPRE(Mat A, Vec x, Vec y)
1119: {
1120:   MatHYPRE_MultKernel_Private(A, 1.0, x, 0.0, y, PETSC_TRUE);
1121:   return 0;
1122: }

1124: static PetscErrorCode MatMult_HYPRE(Mat A, Vec x, Vec y)
1125: {
1126:   MatHYPRE_MultKernel_Private(A, 1.0, x, 0.0, y, PETSC_FALSE);
1127:   return 0;
1128: }

1130: static PetscErrorCode MatMultAdd_HYPRE(Mat A, Vec x, Vec y, Vec z)
1131: {
1132:   if (y != z) VecCopy(y, z);
1133:   MatHYPRE_MultKernel_Private(A, 1.0, x, 1.0, z, PETSC_FALSE);
1134:   return 0;
1135: }

1137: static PetscErrorCode MatMultTransposeAdd_HYPRE(Mat A, Vec x, Vec y, Vec z)
1138: {
1139:   if (y != z) VecCopy(y, z);
1140:   MatHYPRE_MultKernel_Private(A, 1.0, x, 1.0, z, PETSC_TRUE);
1141:   return 0;
1142: }

1144: /* y = a * A * x + b * y or y = a * A^t * x + b * y depending on trans */
1145: static PetscErrorCode MatHYPRE_MultKernel_Private(Mat A, HYPRE_Complex a, Vec x, HYPRE_Complex b, Vec y, PetscBool trans)
1146: {
1147:   Mat_HYPRE          *hA = (Mat_HYPRE *)A->data;
1148:   hypre_ParCSRMatrix *parcsr;
1149:   hypre_ParVector    *hx, *hy;

1151:   if (trans) {
1152:     VecHYPRE_IJVectorPushVecRead(hA->b, x);
1153:     if (b != 0.0) VecHYPRE_IJVectorPushVec(hA->x, y);
1154:     else VecHYPRE_IJVectorPushVecWrite(hA->x, y);
1155:     PetscCallExternal(HYPRE_IJVectorGetObject, hA->b->ij, (void **)&hx);
1156:     PetscCallExternal(HYPRE_IJVectorGetObject, hA->x->ij, (void **)&hy);
1157:   } else {
1158:     VecHYPRE_IJVectorPushVecRead(hA->x, x);
1159:     if (b != 0.0) VecHYPRE_IJVectorPushVec(hA->b, y);
1160:     else VecHYPRE_IJVectorPushVecWrite(hA->b, y);
1161:     PetscCallExternal(HYPRE_IJVectorGetObject, hA->x->ij, (void **)&hx);
1162:     PetscCallExternal(HYPRE_IJVectorGetObject, hA->b->ij, (void **)&hy);
1163:   }
1164:   PetscCallExternal(HYPRE_IJMatrixGetObject, hA->ij, (void **)&parcsr);
1165:   if (trans) {
1166:     PetscCallExternal(hypre_ParCSRMatrixMatvecT, a, parcsr, hx, b, hy);
1167:   } else {
1168:     PetscCallExternal(hypre_ParCSRMatrixMatvec, a, parcsr, hx, b, hy);
1169:   }
1170:   VecHYPRE_IJVectorPopVec(hA->x);
1171:   VecHYPRE_IJVectorPopVec(hA->b);
1172:   return 0;
1173: }

1175: static PetscErrorCode MatDestroy_HYPRE(Mat A)
1176: {
1177:   Mat_HYPRE *hA = (Mat_HYPRE *)A->data;

1179:   VecHYPRE_IJVectorDestroy(&hA->x);
1180:   VecHYPRE_IJVectorDestroy(&hA->b);
1181:   if (hA->ij) {
1182:     if (!hA->inner_free) hypre_IJMatrixObject(hA->ij) = NULL;
1183:     PetscCallExternal(HYPRE_IJMatrixDestroy, hA->ij);
1184:   }
1185:   if (hA->comm) PetscCommRestoreComm(PetscObjectComm((PetscObject)A), &hA->comm);

1187:   MatStashDestroy_Private(&A->stash);
1188:   PetscFree(hA->array);

1190:   if (hA->cooMat) {
1191:     MatDestroy(&hA->cooMat);
1192:     PetscStackCallExternalVoid("hypre_TFree", hypre_TFree(hA->diagJ, hA->memType));
1193:     PetscStackCallExternalVoid("hypre_TFree", hypre_TFree(hA->offdJ, hA->memType));
1194:     PetscStackCallExternalVoid("hypre_TFree", hypre_TFree(hA->diag, hA->memType));
1195:   }

1197:   PetscObjectComposeFunction((PetscObject)A, "MatConvert_hypre_aij_C", NULL);
1198:   PetscObjectComposeFunction((PetscObject)A, "MatConvert_hypre_is_C", NULL);
1199:   PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_seqaij_hypre_C", NULL);
1200:   PetscObjectComposeFunction((PetscObject)A, "MatProductSetFromOptions_mpiaij_hypre_C", NULL);
1201:   PetscObjectComposeFunction((PetscObject)A, "MatHYPRESetPreallocation_C", NULL);
1202:   PetscObjectComposeFunction((PetscObject)A, "MatHYPREGetParCSR_C", NULL);
1203:   PetscObjectComposeFunction((PetscObject)A, "MatSetPreallocationCOO_C", NULL);
1204:   PetscObjectComposeFunction((PetscObject)A, "MatSetValuesCOO_C", NULL);
1205:   PetscFree(A->data);
1206:   return 0;
1207: }

1209: static PetscErrorCode MatSetUp_HYPRE(Mat A)
1210: {
1211:   MatHYPRESetPreallocation(A, PETSC_DEFAULT, NULL, PETSC_DEFAULT, NULL);
1212:   return 0;
1213: }

1215: //TODO FIX hypre_CSRMatrixMatvecOutOfPlace
1216: #if defined(PETSC_HAVE_HYPRE_DEVICE)
1217: static PetscErrorCode MatBindToCPU_HYPRE(Mat A, PetscBool bind)
1218: {
1219:   Mat_HYPRE           *hA   = (Mat_HYPRE *)A->data;
1220:   HYPRE_MemoryLocation hmem = bind ? HYPRE_MEMORY_HOST : HYPRE_MEMORY_DEVICE;

1222:   A->boundtocpu = bind;
1223:   if (hA->ij && hypre_IJMatrixAssembleFlag(hA->ij) && hmem != hypre_IJMatrixMemoryLocation(hA->ij)) {
1224:     hypre_ParCSRMatrix *parcsr;
1225:     PetscCallExternal(HYPRE_IJMatrixGetObject, hA->ij, (void **)&parcsr);
1226:     PetscCallExternal(hypre_ParCSRMatrixMigrate, parcsr, hmem);
1227:   }
1228:   if (hA->x) VecHYPRE_IJBindToCPU(hA->x, bind);
1229:   if (hA->b) VecHYPRE_IJBindToCPU(hA->b, bind);
1230:   return 0;
1231: }
1232: #endif

1234: static PetscErrorCode MatAssemblyEnd_HYPRE(Mat A, MatAssemblyType mode)
1235: {
1236:   Mat_HYPRE   *hA = (Mat_HYPRE *)A->data;
1237:   PetscMPIInt  n;
1238:   PetscInt     i, j, rstart, ncols, flg;
1239:   PetscInt    *row, *col;
1240:   PetscScalar *val;


1244:   if (!A->nooffprocentries) {
1245:     while (1) {
1246:       MatStashScatterGetMesg_Private(&A->stash, &n, &row, &col, &val, &flg);
1247:       if (!flg) break;

1249:       for (i = 0; i < n;) {
1250:         /* Now identify the consecutive vals belonging to the same row */
1251:         for (j = i, rstart = row[j]; j < n; j++) {
1252:           if (row[j] != rstart) break;
1253:         }
1254:         if (j < n) ncols = j - i;
1255:         else ncols = n - i;
1256:         /* Now assemble all these values with a single function call */
1257:         MatSetValues_HYPRE(A, 1, row + i, ncols, col + i, val + i, A->insertmode);

1259:         i = j;
1260:       }
1261:     }
1262:     MatStashScatterEnd_Private(&A->stash);
1263:   }

1265:   PetscCallExternal(HYPRE_IJMatrixAssemble, hA->ij);
1266:   /* The assembly routine destroys the aux_matrix, we recreate it here by calling HYPRE_IJMatrixInitialize */
1267:   /* 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 */
1268:   if (!hA->sorted_full) {
1269:     hypre_AuxParCSRMatrix *aux_matrix;

1271:     /* call destroy just to make sure we do not leak anything */
1272:     aux_matrix = (hypre_AuxParCSRMatrix *)hypre_IJMatrixTranslator(hA->ij);
1273:     PetscCallExternal(hypre_AuxParCSRMatrixDestroy, aux_matrix);
1274:     hypre_IJMatrixTranslator(hA->ij) = NULL;

1276:     /* Initialize with assembled flag -> it only recreates the aux_par_matrix */
1277:     PetscCallExternal(HYPRE_IJMatrixInitialize, hA->ij);
1278:     aux_matrix = (hypre_AuxParCSRMatrix *)hypre_IJMatrixTranslator(hA->ij);
1279:     if (aux_matrix) {
1280:       hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 1; /* see comment in MatHYPRESetPreallocation_HYPRE */
1281: #if PETSC_PKG_HYPRE_VERSION_LT(2, 19, 0)
1282:       PetscCallExternal(hypre_AuxParCSRMatrixInitialize, aux_matrix);
1283: #else
1284:       PetscCallExternal(hypre_AuxParCSRMatrixInitialize_v2, aux_matrix, HYPRE_MEMORY_HOST);
1285: #endif
1286:     }
1287:   }
1288:   {
1289:     hypre_ParCSRMatrix *parcsr;

1291:     PetscCallExternal(HYPRE_IJMatrixGetObject, hA->ij, (void **)&parcsr);
1292:     if (!hypre_ParCSRMatrixCommPkg(parcsr)) PetscCallExternal(hypre_MatvecCommPkgCreate, parcsr);
1293:   }
1294:   if (!hA->x) VecHYPRE_IJVectorCreate(A->cmap, &hA->x);
1295:   if (!hA->b) VecHYPRE_IJVectorCreate(A->rmap, &hA->b);
1296: #if defined(PETSC_HAVE_HYPRE_DEVICE)
1297:   MatBindToCPU_HYPRE(A, A->boundtocpu);
1298: #endif
1299:   return 0;
1300: }

1302: static PetscErrorCode MatGetArray_HYPRE(Mat A, PetscInt size, void **array)
1303: {
1304:   Mat_HYPRE *hA = (Mat_HYPRE *)A->data;


1308:   if (hA->size >= size) {
1309:     *array = hA->array;
1310:   } else {
1311:     PetscFree(hA->array);
1312:     hA->size = size;
1313:     PetscMalloc(hA->size, &hA->array);
1314:     *array = hA->array;
1315:   }

1317:   hA->available = PETSC_FALSE;
1318:   return 0;
1319: }

1321: static PetscErrorCode MatRestoreArray_HYPRE(Mat A, void **array)
1322: {
1323:   Mat_HYPRE *hA = (Mat_HYPRE *)A->data;

1325:   *array        = NULL;
1326:   hA->available = PETSC_TRUE;
1327:   return 0;
1328: }

1330: static PetscErrorCode MatSetValues_HYPRE(Mat A, PetscInt nr, const PetscInt rows[], PetscInt nc, const PetscInt cols[], const PetscScalar v[], InsertMode ins)
1331: {
1332:   Mat_HYPRE     *hA   = (Mat_HYPRE *)A->data;
1333:   PetscScalar   *vals = (PetscScalar *)v;
1334:   HYPRE_Complex *sscr;
1335:   PetscInt      *cscr[2];
1336:   PetscInt       i, nzc;
1337:   void          *array = NULL;

1339:   MatGetArray_HYPRE(A, sizeof(PetscInt) * (2 * nc) + sizeof(HYPRE_Complex) * nc * nr, &array);
1340:   cscr[0] = (PetscInt *)array;
1341:   cscr[1] = ((PetscInt *)array) + nc;
1342:   sscr    = (HYPRE_Complex *)(((PetscInt *)array) + nc * 2);
1343:   for (i = 0, nzc = 0; i < nc; i++) {
1344:     if (cols[i] >= 0) {
1345:       cscr[0][nzc]   = cols[i];
1346:       cscr[1][nzc++] = i;
1347:     }
1348:   }
1349:   if (!nzc) {
1350:     MatRestoreArray_HYPRE(A, &array);
1351:     return 0;
1352:   }

1354: #if 0 //defined(PETSC_HAVE_HYPRE_DEVICE)
1355:   if (HYPRE_MEMORY_HOST != hypre_IJMatrixMemoryLocation(hA->ij)) {
1356:     hypre_ParCSRMatrix *parcsr;

1358:     PetscCallExternal(HYPRE_IJMatrixGetObject,hA->ij,(void**)&parcsr);
1359:     PetscCallExternal(hypre_ParCSRMatrixMigrate,parcsr, HYPRE_MEMORY_HOST);
1360:   }
1361: #endif

1363:   if (ins == ADD_VALUES) {
1364:     for (i = 0; i < nr; i++) {
1365:       if (rows[i] >= 0) {
1366:         PetscInt  j;
1367:         HYPRE_Int hnc = (HYPRE_Int)nzc;

1370:         for (j = 0; j < nzc; j++) PetscHYPREScalarCast(vals[cscr[1][j]], &sscr[j]);
1371:         PetscCallExternal(HYPRE_IJMatrixAddToValues, hA->ij, 1, &hnc, (HYPRE_BigInt *)(rows + i), (HYPRE_BigInt *)cscr[0], sscr);
1372:       }
1373:       vals += nc;
1374:     }
1375:   } else { /* INSERT_VALUES */
1376:     PetscInt rst, ren;

1378:     MatGetOwnershipRange(A, &rst, &ren);
1379:     for (i = 0; i < nr; i++) {
1380:       if (rows[i] >= 0) {
1381:         PetscInt  j;
1382:         HYPRE_Int hnc = (HYPRE_Int)nzc;

1385:         for (j = 0; j < nzc; j++) PetscHYPREScalarCast(vals[cscr[1][j]], &sscr[j]);
1386:         /* nonlocal values */
1387:         if (rows[i] < rst || rows[i] >= ren) MatStashValuesRow_Private(&A->stash, rows[i], nzc, cscr[0], (PetscScalar *)sscr, PETSC_FALSE);
1388:         /* local values */
1389:         else PetscCallExternal(HYPRE_IJMatrixSetValues, hA->ij, 1, &hnc, (HYPRE_BigInt *)(rows + i), (HYPRE_BigInt *)cscr[0], sscr);
1390:       }
1391:       vals += nc;
1392:     }
1393:   }

1395:   MatRestoreArray_HYPRE(A, &array);
1396:   return 0;
1397: }

1399: static PetscErrorCode MatHYPRESetPreallocation_HYPRE(Mat A, PetscInt dnz, const PetscInt dnnz[], PetscInt onz, const PetscInt onnz[])
1400: {
1401:   Mat_HYPRE  *hA = (Mat_HYPRE *)A->data;
1402:   HYPRE_Int  *hdnnz, *honnz;
1403:   PetscInt    i, rs, re, cs, ce, bs;
1404:   PetscMPIInt size;

1406:   PetscLayoutSetUp(A->rmap);
1407:   PetscLayoutSetUp(A->cmap);
1408:   rs = A->rmap->rstart;
1409:   re = A->rmap->rend;
1410:   cs = A->cmap->rstart;
1411:   ce = A->cmap->rend;
1412:   if (!hA->ij) {
1413:     PetscCallExternal(HYPRE_IJMatrixCreate, hA->comm, rs, re - 1, cs, ce - 1, &hA->ij);
1414:     PetscCallExternal(HYPRE_IJMatrixSetObjectType, hA->ij, HYPRE_PARCSR);
1415:   } else {
1416:     HYPRE_BigInt hrs, hre, hcs, hce;
1417:     PetscCallExternal(HYPRE_IJMatrixGetLocalRange, hA->ij, &hrs, &hre, &hcs, &hce);
1420:   }
1421:   MatGetBlockSize(A, &bs);
1422:   if (dnz == PETSC_DEFAULT || dnz == PETSC_DECIDE) dnz = 10 * bs;
1423:   if (onz == PETSC_DEFAULT || onz == PETSC_DECIDE) onz = 10 * bs;

1425:   if (!dnnz) {
1426:     PetscMalloc1(A->rmap->n, &hdnnz);
1427:     for (i = 0; i < A->rmap->n; i++) hdnnz[i] = dnz;
1428:   } else {
1429:     hdnnz = (HYPRE_Int *)dnnz;
1430:   }
1431:   MPI_Comm_size(PetscObjectComm((PetscObject)A), &size);
1432:   if (size > 1) {
1433:     hypre_AuxParCSRMatrix *aux_matrix;
1434:     if (!onnz) {
1435:       PetscMalloc1(A->rmap->n, &honnz);
1436:       for (i = 0; i < A->rmap->n; i++) honnz[i] = onz;
1437:     } else honnz = (HYPRE_Int *)onnz;
1438:     /* SetDiagOffdSizes sets hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 0, since it seems
1439:        they assume the user will input the entire row values, properly sorted
1440:        In PETSc, we don't make such an assumption and set this flag to 1,
1441:        unless the option MAT_SORTED_FULL is set to true.
1442:        Also, to avoid possible memory leaks, we destroy and recreate the translator
1443:        This has to be done here, as HYPRE_IJMatrixInitialize will properly initialize
1444:        the IJ matrix for us */
1445:     aux_matrix = (hypre_AuxParCSRMatrix *)hypre_IJMatrixTranslator(hA->ij);
1446:     hypre_AuxParCSRMatrixDestroy(aux_matrix);
1447:     hypre_IJMatrixTranslator(hA->ij) = NULL;
1448:     PetscCallExternal(HYPRE_IJMatrixSetDiagOffdSizes, hA->ij, hdnnz, honnz);
1449:     aux_matrix                               = (hypre_AuxParCSRMatrix *)hypre_IJMatrixTranslator(hA->ij);
1450:     hypre_AuxParCSRMatrixNeedAux(aux_matrix) = !hA->sorted_full;
1451:   } else {
1452:     honnz = NULL;
1453:     PetscCallExternal(HYPRE_IJMatrixSetRowSizes, hA->ij, hdnnz);
1454:   }

1456:   /* reset assembled flag and call the initialize method */
1457:   hypre_IJMatrixAssembleFlag(hA->ij) = 0;
1458: #if PETSC_PKG_HYPRE_VERSION_LT(2, 19, 0)
1459:   PetscCallExternal(HYPRE_IJMatrixInitialize, hA->ij);
1460: #else
1461:   PetscCallExternal(HYPRE_IJMatrixInitialize_v2, hA->ij, HYPRE_MEMORY_HOST);
1462: #endif
1463:   if (!dnnz) PetscFree(hdnnz);
1464:   if (!onnz && honnz) PetscFree(honnz);
1465:   /* Match AIJ logic */
1466:   A->preallocated = PETSC_TRUE;
1467:   A->assembled    = PETSC_FALSE;
1468:   return 0;
1469: }

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

1474:    Collective

1476:    Input Parameters:
1477: +  A - the matrix
1478: .  dnz  - number of nonzeros per row in DIAGONAL portion of local submatrix
1479:           (same value is used for all local rows)
1480: .  dnnz - array containing the number of nonzeros in the various rows of the
1481:           DIAGONAL portion of the local submatrix (possibly different for each row)
1482:           or NULL (`PETSC_NULL_INTEGER` in Fortran), if d_nz is used to specify the nonzero structure.
1483:           The size of this array is equal to the number of local rows, i.e 'm'.
1484:           For matrices that will be factored, you must leave room for (and set)
1485:           the diagonal entry even if it is zero.
1486: .  onz  - number of nonzeros per row in the OFF-DIAGONAL portion of local
1487:           submatrix (same value is used for all local rows).
1488: -  onnz - array containing the number of nonzeros in the various rows of the
1489:           OFF-DIAGONAL portion of the local submatrix (possibly different for
1490:           each row) or NULL (`PETSC_NULL_INTEGER` in Fortran), if o_nz is used to specify the nonzero
1491:           structure. The size of this array is equal to the number
1492:           of local rows, i.e 'm'.

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

1497:    Level: intermediate

1499: .seealso: `MatCreate()`, `MatMPIAIJSetPreallocation()`, `MATHYPRE`
1500: @*/
1501: PetscErrorCode MatHYPRESetPreallocation(Mat A, PetscInt dnz, const PetscInt dnnz[], PetscInt onz, const PetscInt onnz[])
1502: {
1505:   PetscTryMethod(A, "MatHYPRESetPreallocation_C", (Mat, PetscInt, const PetscInt[], PetscInt, const PetscInt[]), (A, dnz, dnnz, onz, onnz));
1506:   return 0;
1507: }

1509: /*
1510:    MatCreateFromParCSR - Creates a matrix from a hypre_ParCSRMatrix

1512:    Collective

1514:    Input Parameters:
1515: +  parcsr   - the pointer to the hypre_ParCSRMatrix
1516: .  mtype    - matrix type to be created. Currently MATAIJ, MATIS and MATHYPRE are supported.
1517: -  copymode - PETSc copying options

1519:    Output Parameter:
1520: .  A  - the matrix

1522:    Level: intermediate

1524: .seealso: `MatHYPRE`, `PetscCopyMode`
1525: */
1526: PETSC_EXTERN PetscErrorCode MatCreateFromParCSR(hypre_ParCSRMatrix *parcsr, MatType mtype, PetscCopyMode copymode, Mat *A)
1527: {
1528:   Mat        T;
1529:   Mat_HYPRE *hA;
1530:   MPI_Comm   comm;
1531:   PetscInt   rstart, rend, cstart, cend, M, N;
1532:   PetscBool  isseqaij, isseqaijmkl, ismpiaij, isaij, ishyp, isis;

1534:   comm = hypre_ParCSRMatrixComm(parcsr);
1535:   PetscStrcmp(mtype, MATSEQAIJ, &isseqaij);
1536:   PetscStrcmp(mtype, MATSEQAIJMKL, &isseqaijmkl);
1537:   PetscStrcmp(mtype, MATMPIAIJ, &ismpiaij);
1538:   PetscStrcmp(mtype, MATAIJ, &isaij);
1539:   PetscStrcmp(mtype, MATHYPRE, &ishyp);
1540:   PetscStrcmp(mtype, MATIS, &isis);
1541:   isaij = (PetscBool)(isseqaij || isseqaijmkl || ismpiaij || isaij);
1542:   /* TODO */
1544:   /* access ParCSRMatrix */
1545:   rstart = hypre_ParCSRMatrixFirstRowIndex(parcsr);
1546:   rend   = hypre_ParCSRMatrixLastRowIndex(parcsr);
1547:   cstart = hypre_ParCSRMatrixFirstColDiag(parcsr);
1548:   cend   = hypre_ParCSRMatrixLastColDiag(parcsr);
1549:   M      = hypre_ParCSRMatrixGlobalNumRows(parcsr);
1550:   N      = hypre_ParCSRMatrixGlobalNumCols(parcsr);

1552:   /* fix for empty local rows/columns */
1553:   if (rend < rstart) rend = rstart;
1554:   if (cend < cstart) cend = cstart;

1556:   /* PETSc convention */
1557:   rend++;
1558:   cend++;
1559:   rend = PetscMin(rend, M);
1560:   cend = PetscMin(cend, N);

1562:   /* create PETSc matrix with MatHYPRE */
1563:   MatCreate(comm, &T);
1564:   MatSetSizes(T, rend - rstart, cend - cstart, M, N);
1565:   MatSetType(T, MATHYPRE);
1566:   hA = (Mat_HYPRE *)(T->data);

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

1572:   // TODO DEV
1573:   /* create new ParCSR object if needed */
1574:   if (ishyp && copymode == PETSC_COPY_VALUES) {
1575:     hypre_ParCSRMatrix *new_parcsr;
1576: #if PETSC_PKG_HYPRE_VERSION_LT(2, 18, 0)
1577:     hypre_CSRMatrix *hdiag, *hoffd, *ndiag, *noffd;

1579:     new_parcsr = hypre_ParCSRMatrixClone(parcsr, 0);
1580:     hdiag      = hypre_ParCSRMatrixDiag(parcsr);
1581:     hoffd      = hypre_ParCSRMatrixOffd(parcsr);
1582:     ndiag      = hypre_ParCSRMatrixDiag(new_parcsr);
1583:     noffd      = hypre_ParCSRMatrixOffd(new_parcsr);
1584:     PetscArraycpy(hypre_CSRMatrixData(ndiag), hypre_CSRMatrixData(hdiag), hypre_CSRMatrixNumNonzeros(hdiag));
1585:     PetscArraycpy(hypre_CSRMatrixData(noffd), hypre_CSRMatrixData(hoffd), hypre_CSRMatrixNumNonzeros(hoffd));
1586: #else
1587:     new_parcsr = hypre_ParCSRMatrixClone(parcsr, 1);
1588: #endif
1589:     parcsr   = new_parcsr;
1590:     copymode = PETSC_OWN_POINTER;
1591:   }

1593:   /* set ParCSR object */
1594:   hypre_IJMatrixObject(hA->ij) = parcsr;
1595:   T->preallocated              = PETSC_TRUE;

1597:   /* set assembled flag */
1598:   hypre_IJMatrixAssembleFlag(hA->ij) = 1;
1599: #if 0
1600:   PetscCallExternal(HYPRE_IJMatrixInitialize,hA->ij);
1601: #endif
1602:   if (ishyp) {
1603:     PetscMPIInt myid = 0;

1605:     /* make sure we always have row_starts and col_starts available */
1606:     if (HYPRE_AssumedPartitionCheck()) MPI_Comm_rank(comm, &myid);
1607: #if defined(hypre_ParCSRMatrixOwnsRowStarts)
1608:     if (!hypre_ParCSRMatrixOwnsColStarts(parcsr)) {
1609:       PetscLayout map;

1611:       MatGetLayouts(T, NULL, &map);
1612:       PetscLayoutSetUp(map);
1613:       hypre_ParCSRMatrixColStarts(parcsr) = (HYPRE_BigInt *)(map->range + myid);
1614:     }
1615:     if (!hypre_ParCSRMatrixOwnsRowStarts(parcsr)) {
1616:       PetscLayout map;

1618:       MatGetLayouts(T, &map, NULL);
1619:       PetscLayoutSetUp(map);
1620:       hypre_ParCSRMatrixRowStarts(parcsr) = (HYPRE_BigInt *)(map->range + myid);
1621:     }
1622: #endif
1623:     /* prevent from freeing the pointer */
1624:     if (copymode == PETSC_USE_POINTER) hA->inner_free = PETSC_FALSE;
1625:     *A = T;
1626:     MatSetOption(*A, MAT_SORTED_FULL, PETSC_TRUE);
1627:     MatAssemblyBegin(*A, MAT_FINAL_ASSEMBLY);
1628:     MatAssemblyEnd(*A, MAT_FINAL_ASSEMBLY);
1629:   } else if (isaij) {
1630:     if (copymode != PETSC_OWN_POINTER) {
1631:       /* prevent from freeing the pointer */
1632:       hA->inner_free = PETSC_FALSE;
1633:       MatConvert_HYPRE_AIJ(T, MATAIJ, MAT_INITIAL_MATRIX, A);
1634:       MatDestroy(&T);
1635:     } else { /* AIJ return type with PETSC_OWN_POINTER */
1636:       MatConvert_HYPRE_AIJ(T, MATAIJ, MAT_INPLACE_MATRIX, &T);
1637:       *A = T;
1638:     }
1639:   } else if (isis) {
1640:     MatConvert_HYPRE_IS(T, MATIS, MAT_INITIAL_MATRIX, A);
1641:     if (copymode != PETSC_OWN_POINTER) hA->inner_free = PETSC_FALSE;
1642:     MatDestroy(&T);
1643:   }
1644:   return 0;
1645: }

1647: static PetscErrorCode MatHYPREGetParCSR_HYPRE(Mat A, hypre_ParCSRMatrix **parcsr)
1648: {
1649:   Mat_HYPRE *hA = (Mat_HYPRE *)A->data;
1650:   HYPRE_Int  type;

1653:   PetscCallExternal(HYPRE_IJMatrixGetObjectType, hA->ij, &type);
1655:   PetscCallExternal(HYPRE_IJMatrixGetObject, hA->ij, (void **)parcsr);
1656:   return 0;
1657: }

1659: /*
1660:    MatHYPREGetParCSR - Gets the pointer to the ParCSR matrix

1662:    Not collective

1664:    Input Parameters:
1665: +  A  - the MATHYPRE object

1667:    Output Parameter:
1668: .  parcsr  - the pointer to the hypre_ParCSRMatrix

1670:    Level: intermediate

1672: .seealso: `MatHYPRE`, `PetscCopyMode`
1673: */
1674: PetscErrorCode MatHYPREGetParCSR(Mat A, hypre_ParCSRMatrix **parcsr)
1675: {
1678:   PetscUseMethod(A, "MatHYPREGetParCSR_C", (Mat, hypre_ParCSRMatrix **), (A, parcsr));
1679:   return 0;
1680: }

1682: static PetscErrorCode MatMissingDiagonal_HYPRE(Mat A, PetscBool *missing, PetscInt *dd)
1683: {
1684:   hypre_ParCSRMatrix *parcsr;
1685:   hypre_CSRMatrix    *ha;
1686:   PetscInt            rst;

1689:   MatGetOwnershipRange(A, &rst, NULL);
1690:   MatHYPREGetParCSR_HYPRE(A, &parcsr);
1691:   if (missing) *missing = PETSC_FALSE;
1692:   if (dd) *dd = -1;
1693:   ha = hypre_ParCSRMatrixDiag(parcsr);
1694:   if (ha) {
1695:     PetscInt   size, i;
1696:     HYPRE_Int *ii, *jj;

1698:     size = hypre_CSRMatrixNumRows(ha);
1699:     ii   = hypre_CSRMatrixI(ha);
1700:     jj   = hypre_CSRMatrixJ(ha);
1701:     for (i = 0; i < size; i++) {
1702:       PetscInt  j;
1703:       PetscBool found = PETSC_FALSE;

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

1707:       if (!found) {
1708:         PetscInfo(A, "Matrix is missing local diagonal entry %" PetscInt_FMT "\n", i);
1709:         if (missing) *missing = PETSC_TRUE;
1710:         if (dd) *dd = i + rst;
1711:         return 0;
1712:       }
1713:     }
1714:     if (!size) {
1715:       PetscInfo(A, "Matrix has no diagonal entries therefore is missing diagonal\n");
1716:       if (missing) *missing = PETSC_TRUE;
1717:       if (dd) *dd = rst;
1718:     }
1719:   } else {
1720:     PetscInfo(A, "Matrix has no diagonal entries therefore is missing diagonal\n");
1721:     if (missing) *missing = PETSC_TRUE;
1722:     if (dd) *dd = rst;
1723:   }
1724:   return 0;
1725: }

1727: static PetscErrorCode MatScale_HYPRE(Mat A, PetscScalar s)
1728: {
1729:   hypre_ParCSRMatrix *parcsr;
1730: #if PETSC_PKG_HYPRE_VERSION_LT(2, 19, 0)
1731:   hypre_CSRMatrix *ha;
1732: #endif
1733:   HYPRE_Complex hs;

1735:   PetscHYPREScalarCast(s, &hs);
1736:   MatHYPREGetParCSR_HYPRE(A, &parcsr);
1737: #if PETSC_PKG_HYPRE_VERSION_GE(2, 19, 0)
1738:   PetscCallExternal(hypre_ParCSRMatrixScale, parcsr, hs);
1739: #else /* diagonal part */
1740:   ha = hypre_ParCSRMatrixDiag(parcsr);
1741:   if (ha) {
1742:     PetscInt size, i;
1743:     HYPRE_Int *ii;
1744:     HYPRE_Complex *a;

1746:     size = hypre_CSRMatrixNumRows(ha);
1747:     a = hypre_CSRMatrixData(ha);
1748:     ii = hypre_CSRMatrixI(ha);
1749:     for (i = 0; i < ii[size]; i++) a[i] *= hs;
1750:   }
1751:   /* offdiagonal part */
1752:   ha = hypre_ParCSRMatrixOffd(parcsr);
1753:   if (ha) {
1754:     PetscInt size, i;
1755:     HYPRE_Int *ii;
1756:     HYPRE_Complex *a;

1758:     size = hypre_CSRMatrixNumRows(ha);
1759:     a = hypre_CSRMatrixData(ha);
1760:     ii = hypre_CSRMatrixI(ha);
1761:     for (i = 0; i < ii[size]; i++) a[i] *= hs;
1762:   }
1763: #endif
1764:   return 0;
1765: }

1767: static PetscErrorCode MatZeroRowsColumns_HYPRE(Mat A, PetscInt numRows, const PetscInt rows[], PetscScalar diag, Vec x, Vec b)
1768: {
1769:   hypre_ParCSRMatrix *parcsr;
1770:   HYPRE_Int          *lrows;
1771:   PetscInt            rst, ren, i;

1774:   MatHYPREGetParCSR_HYPRE(A, &parcsr);
1775:   PetscMalloc1(numRows, &lrows);
1776:   MatGetOwnershipRange(A, &rst, &ren);
1777:   for (i = 0; i < numRows; i++) {
1779:     lrows[i] = rows[i] - rst;
1780:   }
1781:   PetscCallExternal(hypre_ParCSRMatrixEliminateRowsCols, parcsr, numRows, lrows);
1782:   PetscFree(lrows);
1783:   return 0;
1784: }

1786: static PetscErrorCode MatZeroEntries_HYPRE_CSRMatrix(hypre_CSRMatrix *ha)
1787: {
1788:   if (ha) {
1789:     HYPRE_Int     *ii, size;
1790:     HYPRE_Complex *a;

1792:     size = hypre_CSRMatrixNumRows(ha);
1793:     a    = hypre_CSRMatrixData(ha);
1794:     ii   = hypre_CSRMatrixI(ha);

1796:     if (a) PetscArrayzero(a, ii[size]);
1797:   }
1798:   return 0;
1799: }

1801: PetscErrorCode MatZeroEntries_HYPRE(Mat A)
1802: {
1803:   Mat_HYPRE *hA = (Mat_HYPRE *)A->data;

1805:   if (HYPRE_MEMORY_DEVICE == hypre_IJMatrixMemoryLocation(hA->ij)) {
1806:     PetscCallExternal(HYPRE_IJMatrixSetConstantValues, hA->ij, 0.0);
1807:   } else {
1808:     hypre_ParCSRMatrix *parcsr;

1810:     MatHYPREGetParCSR_HYPRE(A, &parcsr);
1811:     MatZeroEntries_HYPRE_CSRMatrix(hypre_ParCSRMatrixDiag(parcsr));
1812:     MatZeroEntries_HYPRE_CSRMatrix(hypre_ParCSRMatrixOffd(parcsr));
1813:   }
1814:   return 0;
1815: }

1817: static PetscErrorCode MatZeroRows_HYPRE_CSRMatrix(hypre_CSRMatrix *hA, PetscInt N, const PetscInt rows[], HYPRE_Complex diag)
1818: {
1819:   PetscInt       ii;
1820:   HYPRE_Int     *i, *j;
1821:   HYPRE_Complex *a;

1823:   if (!hA) return 0;

1825:   i = hypre_CSRMatrixI(hA);
1826:   j = hypre_CSRMatrixJ(hA);
1827:   a = hypre_CSRMatrixData(hA);

1829:   for (ii = 0; ii < N; ii++) {
1830:     HYPRE_Int jj, ibeg, iend, irow;

1832:     irow = rows[ii];
1833:     ibeg = i[irow];
1834:     iend = i[irow + 1];
1835:     for (jj = ibeg; jj < iend; jj++)
1836:       if (j[jj] == irow) a[jj] = diag;
1837:       else a[jj] = 0.0;
1838:   }
1839:   return 0;
1840: }

1842: static PetscErrorCode MatZeroRows_HYPRE(Mat A, PetscInt N, const PetscInt rows[], PetscScalar diag, Vec x, Vec b)
1843: {
1844:   hypre_ParCSRMatrix *parcsr;
1845:   PetscInt           *lrows, len;
1846:   HYPRE_Complex       hdiag;

1849:   PetscHYPREScalarCast(diag, &hdiag);
1850:   /* retrieve the internal matrix */
1851:   MatHYPREGetParCSR_HYPRE(A, &parcsr);
1852:   /* get locally owned rows */
1853:   MatZeroRowsMapLocal_Private(A, N, rows, &len, &lrows);
1854:   /* zero diagonal part */
1855:   MatZeroRows_HYPRE_CSRMatrix(hypre_ParCSRMatrixDiag(parcsr), len, lrows, hdiag);
1856:   /* zero off-diagonal part */
1857:   MatZeroRows_HYPRE_CSRMatrix(hypre_ParCSRMatrixOffd(parcsr), len, lrows, 0.0);

1859:   PetscFree(lrows);
1860:   return 0;
1861: }

1863: static PetscErrorCode MatAssemblyBegin_HYPRE(Mat mat, MatAssemblyType mode)
1864: {
1865:   if (mat->nooffprocentries) return 0;

1867:   MatStashScatterBegin_Private(mat, &mat->stash, mat->rmap->range);
1868:   return 0;
1869: }

1871: static PetscErrorCode MatGetRow_HYPRE(Mat A, PetscInt row, PetscInt *nz, PetscInt **idx, PetscScalar **v)
1872: {
1873:   hypre_ParCSRMatrix *parcsr;
1874:   HYPRE_Int           hnz;

1876:   /* retrieve the internal matrix */
1877:   MatHYPREGetParCSR_HYPRE(A, &parcsr);
1878:   /* call HYPRE API */
1879:   PetscCallExternal(HYPRE_ParCSRMatrixGetRow, parcsr, row, &hnz, (HYPRE_BigInt **)idx, (HYPRE_Complex **)v);
1880:   if (nz) *nz = (PetscInt)hnz;
1881:   return 0;
1882: }

1884: static PetscErrorCode MatRestoreRow_HYPRE(Mat A, PetscInt row, PetscInt *nz, PetscInt **idx, PetscScalar **v)
1885: {
1886:   hypre_ParCSRMatrix *parcsr;
1887:   HYPRE_Int           hnz;

1889:   /* retrieve the internal matrix */
1890:   MatHYPREGetParCSR_HYPRE(A, &parcsr);
1891:   /* call HYPRE API */
1892:   hnz = nz ? (HYPRE_Int)(*nz) : 0;
1893:   PetscCallExternal(HYPRE_ParCSRMatrixRestoreRow, parcsr, row, &hnz, (HYPRE_BigInt **)idx, (HYPRE_Complex **)v);
1894:   return 0;
1895: }

1897: static PetscErrorCode MatGetValues_HYPRE(Mat A, PetscInt m, const PetscInt idxm[], PetscInt n, const PetscInt idxn[], PetscScalar v[])
1898: {
1899:   Mat_HYPRE *hA = (Mat_HYPRE *)A->data;
1900:   PetscInt   i;

1902:   if (!m || !n) return 0;
1903:   /* Ignore negative row indices
1904:    * And negative column indices should be automatically ignored in hypre
1905:    * */
1906:   for (i = 0; i < m; i++) {
1907:     if (idxm[i] >= 0) {
1908:       HYPRE_Int hn = (HYPRE_Int)n;
1909:       PetscCallExternal(HYPRE_IJMatrixGetValues, hA->ij, 1, &hn, (HYPRE_BigInt *)&idxm[i], (HYPRE_BigInt *)idxn, (HYPRE_Complex *)(v + i * n));
1910:     }
1911:   }
1912:   return 0;
1913: }

1915: static PetscErrorCode MatSetOption_HYPRE(Mat A, MatOption op, PetscBool flg)
1916: {
1917:   Mat_HYPRE *hA = (Mat_HYPRE *)A->data;

1919:   switch (op) {
1920:   case MAT_NO_OFF_PROC_ENTRIES:
1921:     if (flg) PetscCallExternal(HYPRE_IJMatrixSetMaxOffProcElmts, hA->ij, 0);
1922:     break;
1923:   case MAT_SORTED_FULL:
1924:     hA->sorted_full = flg;
1925:     break;
1926:   default:
1927:     break;
1928:   }
1929:   return 0;
1930: }

1932: static PetscErrorCode MatView_HYPRE(Mat A, PetscViewer view)
1933: {
1934:   PetscViewerFormat format;

1936:   PetscViewerGetFormat(view, &format);
1937:   if (format == PETSC_VIEWER_ASCII_FACTOR_INFO || format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) return 0;
1938:   if (format != PETSC_VIEWER_NATIVE) {
1939:     Mat                 B;
1940:     hypre_ParCSRMatrix *parcsr;
1941:     PetscErrorCode (*mview)(Mat, PetscViewer) = NULL;

1943:     MatHYPREGetParCSR_HYPRE(A, &parcsr);
1944:     MatCreateFromParCSR(parcsr, MATAIJ, PETSC_USE_POINTER, &B);
1945:     MatGetOperation(B, MATOP_VIEW, (void (**)(void)) & mview);
1947:     (*mview)(B, view);
1948:     MatDestroy(&B);
1949:   } else {
1950:     Mat_HYPRE  *hA = (Mat_HYPRE *)A->data;
1951:     PetscMPIInt size;
1952:     PetscBool   isascii;
1953:     const char *filename;

1955:     /* HYPRE uses only text files */
1956:     PetscObjectTypeCompare((PetscObject)view, PETSCVIEWERASCII, &isascii);
1958:     PetscViewerFileGetName(view, &filename);
1959:     PetscCallExternal(HYPRE_IJMatrixPrint, hA->ij, filename);
1960:     MPI_Comm_size(hA->comm, &size);
1961:     if (size > 1) {
1962:       PetscViewerASCIIPrintf(view, "Matrix files: %s.%05d ... %s.%05d\n", filename, 0, filename, size - 1);
1963:     } else {
1964:       PetscViewerASCIIPrintf(view, "Matrix file: %s.%05d\n", filename, 0);
1965:     }
1966:   }
1967:   return 0;
1968: }

1970: static PetscErrorCode MatDuplicate_HYPRE(Mat A, MatDuplicateOption op, Mat *B)
1971: {
1972:   hypre_ParCSRMatrix *parcsr = NULL;
1973:   PetscCopyMode       cpmode;

1975:   MatHYPREGetParCSR_HYPRE(A, &parcsr);
1976:   if (op == MAT_DO_NOT_COPY_VALUES || op == MAT_SHARE_NONZERO_PATTERN) {
1977:     parcsr = hypre_ParCSRMatrixClone(parcsr, 0);
1978:     cpmode = PETSC_OWN_POINTER;
1979:   } else {
1980:     cpmode = PETSC_COPY_VALUES;
1981:   }
1982:   MatCreateFromParCSR(parcsr, MATHYPRE, cpmode, B);
1983:   return 0;
1984: }

1986: static PetscErrorCode MatCopy_HYPRE(Mat A, Mat B, MatStructure str)
1987: {
1988:   hypre_ParCSRMatrix *acsr, *bcsr;

1990:   if (str == SAME_NONZERO_PATTERN && A->ops->copy == B->ops->copy) {
1991:     MatHYPREGetParCSR_HYPRE(A, &acsr);
1992:     MatHYPREGetParCSR_HYPRE(B, &bcsr);
1993:     PetscCallExternal(hypre_ParCSRMatrixCopy, acsr, bcsr, 1);
1994:     MatSetOption(B, MAT_SORTED_FULL, PETSC_TRUE); /* "perfect" preallocation, so no need for hypre_AuxParCSRMatrixNeedAux */
1995:     MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
1996:     MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
1997:   } else {
1998:     MatCopy_Basic(A, B, str);
1999:   }
2000:   return 0;
2001: }

2003: static PetscErrorCode MatGetDiagonal_HYPRE(Mat A, Vec d)
2004: {
2005:   hypre_ParCSRMatrix *parcsr;
2006:   hypre_CSRMatrix    *dmat;
2007:   HYPRE_Complex      *a;
2008:   HYPRE_Complex      *data = NULL;
2009:   HYPRE_Int          *diag = NULL;
2010:   PetscInt            i;
2011:   PetscBool           cong;

2013:   MatHasCongruentLayouts(A, &cong);
2015:   if (PetscDefined(USE_DEBUG)) {
2016:     PetscBool miss;
2017:     MatMissingDiagonal(A, &miss, NULL);
2019:   }
2020:   MatHYPREGetParCSR_HYPRE(A, &parcsr);
2021:   dmat = hypre_ParCSRMatrixDiag(parcsr);
2022:   if (dmat) {
2023:     /* this cast fixes the clang error: implicit conversion from 'HYPRE_Complex' (aka '_Complex double') to 'double' is not permitted in C++ */
2024:     VecGetArray(d, (PetscScalar **)&a);
2025:     diag = hypre_CSRMatrixI(dmat);
2026:     data = hypre_CSRMatrixData(dmat);
2027:     for (i = 0; i < A->rmap->n; i++) a[i] = data[diag[i]];
2028:     VecRestoreArray(d, (PetscScalar **)&a);
2029:   }
2030:   return 0;
2031: }

2033: #include <petscblaslapack.h>

2035: static PetscErrorCode MatAXPY_HYPRE(Mat Y, PetscScalar a, Mat X, MatStructure str)
2036: {
2037: #if defined(PETSC_HAVE_HYPRE_DEVICE)
2038:   {
2039:     Mat                 B;
2040:     hypre_ParCSRMatrix *x, *y, *z;

2042:     MatHYPREGetParCSR(Y, &y);
2043:     MatHYPREGetParCSR(X, &x);
2044:     PetscCallExternal(hypre_ParCSRMatrixAdd, 1.0, y, 1.0, x, &z);
2045:     MatCreateFromParCSR(z, MATHYPRE, PETSC_OWN_POINTER, &B);
2046:     MatHeaderMerge(Y, &B);
2047:   }
2048: #else
2049:   if (str == SAME_NONZERO_PATTERN) {
2050:     hypre_ParCSRMatrix *x, *y;
2051:     hypre_CSRMatrix *xloc, *yloc;
2052:     PetscInt xnnz, ynnz;
2053:     HYPRE_Complex *xarr, *yarr;
2054:     PetscBLASInt one = 1, bnz;

2056:     MatHYPREGetParCSR(Y, &y);
2057:     MatHYPREGetParCSR(X, &x);

2059:     /* diagonal block */
2060:     xloc = hypre_ParCSRMatrixDiag(x);
2061:     yloc = hypre_ParCSRMatrixDiag(y);
2062:     xnnz = 0;
2063:     ynnz = 0;
2064:     xarr = NULL;
2065:     yarr = NULL;
2066:     if (xloc) {
2067:       xarr = hypre_CSRMatrixData(xloc);
2068:       xnnz = hypre_CSRMatrixNumNonzeros(xloc);
2069:     }
2070:     if (yloc) {
2071:       yarr = hypre_CSRMatrixData(yloc);
2072:       ynnz = hypre_CSRMatrixNumNonzeros(yloc);
2073:     }
2075:     PetscBLASIntCast(xnnz, &bnz);
2076:     PetscCallBLAS("BLASaxpy", BLASaxpy_(&bnz, &a, (PetscScalar *)xarr, &one, (PetscScalar *)yarr, &one));

2078:     /* off-diagonal block */
2079:     xloc = hypre_ParCSRMatrixOffd(x);
2080:     yloc = hypre_ParCSRMatrixOffd(y);
2081:     xnnz = 0;
2082:     ynnz = 0;
2083:     xarr = NULL;
2084:     yarr = NULL;
2085:     if (xloc) {
2086:       xarr = hypre_CSRMatrixData(xloc);
2087:       xnnz = hypre_CSRMatrixNumNonzeros(xloc);
2088:     }
2089:     if (yloc) {
2090:       yarr = hypre_CSRMatrixData(yloc);
2091:       ynnz = hypre_CSRMatrixNumNonzeros(yloc);
2092:     }
2094:     PetscBLASIntCast(xnnz, &bnz);
2095:     PetscCallBLAS("BLASaxpy", BLASaxpy_(&bnz, &a, (PetscScalar *)xarr, &one, (PetscScalar *)yarr, &one));
2096:   } else if (str == SUBSET_NONZERO_PATTERN) {
2097:     MatAXPY_Basic(Y, a, X, str);
2098:   } else {
2099:     Mat B;

2101:     MatAXPY_Basic_Preallocate(Y, X, &B);
2102:     MatAXPY_BasicWithPreallocation(B, Y, a, X, str);
2103:     MatHeaderReplace(Y, &B);
2104:   }
2105: #endif
2106:   return 0;
2107: }

2109: static PetscErrorCode MatSetPreallocationCOO_HYPRE(Mat mat, PetscCount coo_n, PetscInt coo_i[], PetscInt coo_j[])
2110: {
2111:   MPI_Comm             comm;
2112:   PetscMPIInt          size;
2113:   PetscLayout          rmap, cmap;
2114:   Mat_HYPRE           *hmat;
2115:   hypre_ParCSRMatrix  *parCSR;
2116:   hypre_CSRMatrix     *diag, *offd;
2117:   Mat                  A, B, cooMat;
2118:   PetscScalar         *Aa, *Ba;
2119:   HYPRE_MemoryLocation hypreMemtype = HYPRE_MEMORY_HOST;
2120:   PetscMemType         petscMemtype;
2121:   MatType              matType = MATAIJ; /* default type of cooMat */

2123:   /* Build an agent matrix cooMat whose type is either MATAIJ or MATAIJKOKKOS.
2124:      It has the same sparsity pattern as mat, and also shares the data array with mat. We use cooMat to do the COO work.
2125:    */
2126:   PetscObjectGetComm((PetscObject)mat, &comm);
2127:   MPI_Comm_size(comm, &size);
2128:   PetscLayoutSetUp(mat->rmap);
2129:   PetscLayoutSetUp(mat->cmap);
2130:   MatGetLayouts(mat, &rmap, &cmap);

2132:   /* I do not know how hypre_ParCSRMatrix stores diagonal elements for non-square matrices, so I just give up now */

2135: #if defined(PETSC_HAVE_DEVICE)
2136:   if (!mat->boundtocpu) { /* mat will be on device, so will cooMat */
2137:   #if defined(PETSC_HAVE_KOKKOS)
2138:     matType = MATAIJKOKKOS;
2139:   #else
2140:     SETERRQ(comm, PETSC_ERR_SUP, "To support MATHYPRE COO assembly on device, we need Kokkos, e.g., --download-kokkos --download-kokkos-kernels");
2141:   #endif
2142:   }
2143: #endif

2145:   /* Do COO preallocation through cooMat */
2146:   hmat = (Mat_HYPRE *)mat->data;
2147:   MatDestroy(&hmat->cooMat);
2148:   MatCreate(comm, &cooMat);
2149:   MatSetType(cooMat, matType);
2150:   MatSetLayouts(cooMat, rmap, cmap);
2151:   MatSetPreallocationCOO(cooMat, coo_n, coo_i, coo_j);

2153:   /* Copy the sparsity pattern from cooMat to hypre IJMatrix hmat->ij */
2154:   MatSetOption(mat, MAT_SORTED_FULL, PETSC_TRUE);
2155:   MatSetOption(mat, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE);
2156:   MatHYPRE_CreateFromMat(cooMat, hmat);    /* Create hmat->ij and preallocate it */
2157:   MatHYPRE_IJMatrixCopy(cooMat, hmat->ij); /* Copy A's (a,i,j) to hmat->ij. To reuse code. Copying 'a' is not really needed */

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

2163:   /* Alias cooMat's data array to IJMatrix's */
2164:   PetscCallExternal(HYPRE_IJMatrixGetObject, hmat->ij, (void **)&parCSR);
2165:   diag = hypre_ParCSRMatrixDiag(parCSR);
2166:   offd = hypre_ParCSRMatrixOffd(parCSR);

2168:   hypreMemtype = hypre_CSRMatrixMemoryLocation(diag);
2169:   A            = (size == 1) ? cooMat : ((Mat_MPIAIJ *)cooMat->data)->A;
2170:   MatSeqAIJGetCSRAndMemType(A, NULL, NULL, &Aa, &petscMemtype);
2171:   PetscAssert((PetscMemTypeHost(petscMemtype) && hypreMemtype == HYPRE_MEMORY_HOST) || (PetscMemTypeDevice(petscMemtype) && hypreMemtype == HYPRE_MEMORY_DEVICE), comm, PETSC_ERR_PLIB, "PETSc and hypre's memory types mismatch");

2173:   hmat->diagJ = hypre_CSRMatrixJ(diag);
2174:   PetscStackCallExternalVoid("hypre_TFree", hypre_TFree(hypre_CSRMatrixData(diag), hypreMemtype));
2175:   hypre_CSRMatrixData(diag)     = (HYPRE_Complex *)Aa;
2176:   hypre_CSRMatrixOwnsData(diag) = 0; /* Take ownership of (j,a) away from hypre. As a result, we need to free them on our own */

2178:   /* Copy diagonal pointers of A to device to facilitate MatSeqAIJMoveDiagonalValuesFront_SeqAIJKokkos */
2179:   if (hypreMemtype == HYPRE_MEMORY_DEVICE) {
2180:     PetscStackCallExternalVoid("hypre_TAlloc", hmat->diag = hypre_TAlloc(PetscInt, rmap->n, hypreMemtype));
2181:     MatMarkDiagonal_SeqAIJ(A); /* We need updated diagonal positions */
2182:     PetscStackCallExternalVoid("hypre_TMemcpy", hypre_TMemcpy(hmat->diag, ((Mat_SeqAIJ *)A->data)->diag, PetscInt, rmap->n, hypreMemtype, HYPRE_MEMORY_HOST));
2183:   }

2185:   if (size > 1) {
2186:     B = ((Mat_MPIAIJ *)cooMat->data)->B;
2187:     MatSeqAIJGetCSRAndMemType(B, NULL, NULL, &Ba, &petscMemtype);
2188:     hmat->offdJ = hypre_CSRMatrixJ(offd);
2189:     PetscStackCallExternalVoid("hypre_TFree", hypre_TFree(hypre_CSRMatrixData(offd), hypreMemtype));
2190:     hypre_CSRMatrixData(offd)     = (HYPRE_Complex *)Ba;
2191:     hypre_CSRMatrixOwnsData(offd) = 0;
2192:   }

2194:   /* Record cooMat for use in MatSetValuesCOO_HYPRE */
2195:   hmat->cooMat  = cooMat;
2196:   hmat->memType = hypreMemtype;
2197:   return 0;
2198: }

2200: static PetscErrorCode MatSetValuesCOO_HYPRE(Mat mat, const PetscScalar v[], InsertMode imode)
2201: {
2202:   Mat_HYPRE  *hmat = (Mat_HYPRE *)mat->data;
2203:   PetscMPIInt size;
2204:   Mat         A;

2207:   MPI_Comm_size(hmat->comm, &size);
2208:   MatSetValuesCOO(hmat->cooMat, v, imode);

2210:   /* Move diagonal elements of the diagonal block to the front of their row, as needed by ParCSRMatrix. So damn hacky */
2211:   A = (size == 1) ? hmat->cooMat : ((Mat_MPIAIJ *)hmat->cooMat->data)->A;
2212:   if (hmat->memType == HYPRE_MEMORY_HOST) {
2213:     Mat_SeqAIJ  *aij = (Mat_SeqAIJ *)A->data;
2214:     PetscInt     i, m, *Ai = aij->i, *Adiag = aij->diag;
2215:     PetscScalar *Aa = aij->a, tmp;

2217:     MatGetSize(A, &m, NULL);
2218:     for (i = 0; i < m; i++) {
2219:       if (Adiag[i] >= Ai[i] && Adiag[i] < Ai[i + 1]) { /* Digonal element of this row exists in a[] and j[] */
2220:         tmp          = Aa[Ai[i]];
2221:         Aa[Ai[i]]    = Aa[Adiag[i]];
2222:         Aa[Adiag[i]] = tmp;
2223:       }
2224:     }
2225:   } else {
2226: #if defined(PETSC_HAVE_KOKKOS_KERNELS)
2227:     MatSeqAIJMoveDiagonalValuesFront_SeqAIJKokkos(A, hmat->diag);
2228: #endif
2229:   }
2230:   return 0;
2231: }

2233: /*MC
2234:    MATHYPRE - MATHYPRE = "hypre" - A matrix type to be used for sequential and parallel sparse matrices
2235:           based on the hypre IJ interface.

2237:    Level: intermediate

2239: .seealso: `MatCreate()`, `MatHYPRESetPreallocation`
2240: M*/

2242: PETSC_EXTERN PetscErrorCode MatCreate_HYPRE(Mat B)
2243: {
2244:   Mat_HYPRE *hB;

2246:   PetscNew(&hB);

2248:   hB->inner_free  = PETSC_TRUE;
2249:   hB->available   = PETSC_TRUE;
2250:   hB->sorted_full = PETSC_FALSE; /* no assumption whether column indices are sorted or not */
2251:   hB->size        = 0;
2252:   hB->array       = NULL;

2254:   B->data      = (void *)hB;
2255:   B->assembled = PETSC_FALSE;

2257:   PetscMemzero(B->ops, sizeof(struct _MatOps));
2258:   B->ops->mult                  = MatMult_HYPRE;
2259:   B->ops->multtranspose         = MatMultTranspose_HYPRE;
2260:   B->ops->multadd               = MatMultAdd_HYPRE;
2261:   B->ops->multtransposeadd      = MatMultTransposeAdd_HYPRE;
2262:   B->ops->setup                 = MatSetUp_HYPRE;
2263:   B->ops->destroy               = MatDestroy_HYPRE;
2264:   B->ops->assemblyend           = MatAssemblyEnd_HYPRE;
2265:   B->ops->assemblybegin         = MatAssemblyBegin_HYPRE;
2266:   B->ops->setvalues             = MatSetValues_HYPRE;
2267:   B->ops->missingdiagonal       = MatMissingDiagonal_HYPRE;
2268:   B->ops->scale                 = MatScale_HYPRE;
2269:   B->ops->zerorowscolumns       = MatZeroRowsColumns_HYPRE;
2270:   B->ops->zeroentries           = MatZeroEntries_HYPRE;
2271:   B->ops->zerorows              = MatZeroRows_HYPRE;
2272:   B->ops->getrow                = MatGetRow_HYPRE;
2273:   B->ops->restorerow            = MatRestoreRow_HYPRE;
2274:   B->ops->getvalues             = MatGetValues_HYPRE;
2275:   B->ops->setoption             = MatSetOption_HYPRE;
2276:   B->ops->duplicate             = MatDuplicate_HYPRE;
2277:   B->ops->copy                  = MatCopy_HYPRE;
2278:   B->ops->view                  = MatView_HYPRE;
2279:   B->ops->getdiagonal           = MatGetDiagonal_HYPRE;
2280:   B->ops->axpy                  = MatAXPY_HYPRE;
2281:   B->ops->productsetfromoptions = MatProductSetFromOptions_HYPRE;
2282: #if defined(PETSC_HAVE_HYPRE_DEVICE)
2283:   B->ops->bindtocpu = MatBindToCPU_HYPRE;
2284:   B->boundtocpu     = PETSC_FALSE;
2285: #endif

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

2290:   PetscCommGetComm(PetscObjectComm((PetscObject)B), &hB->comm);
2291:   PetscObjectChangeTypeName((PetscObject)B, MATHYPRE);
2292:   PetscObjectComposeFunction((PetscObject)B, "MatConvert_hypre_aij_C", MatConvert_HYPRE_AIJ);
2293:   PetscObjectComposeFunction((PetscObject)B, "MatConvert_hypre_is_C", MatConvert_HYPRE_IS);
2294:   PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_seqaij_hypre_C", MatProductSetFromOptions_HYPRE);
2295:   PetscObjectComposeFunction((PetscObject)B, "MatProductSetFromOptions_mpiaij_hypre_C", MatProductSetFromOptions_HYPRE);
2296:   PetscObjectComposeFunction((PetscObject)B, "MatHYPRESetPreallocation_C", MatHYPRESetPreallocation_HYPRE);
2297:   PetscObjectComposeFunction((PetscObject)B, "MatHYPREGetParCSR_C", MatHYPREGetParCSR_HYPRE);
2298:   PetscObjectComposeFunction((PetscObject)B, "MatSetPreallocationCOO_C", MatSetPreallocationCOO_HYPRE);
2299:   PetscObjectComposeFunction((PetscObject)B, "MatSetValuesCOO_C", MatSetValuesCOO_HYPRE);
2300: #if defined(PETSC_HAVE_HYPRE_DEVICE)
2301:   #if defined(HYPRE_USING_HIP)
2302:   PetscDeviceInitialize(PETSC_DEVICE_HIP);
2303:   MatSetVecType(B, VECHIP);
2304:   #endif
2305:   #if defined(HYPRE_USING_CUDA)
2306:   PetscDeviceInitialize(PETSC_DEVICE_CUDA);
2307:   MatSetVecType(B, VECCUDA);
2308:   #endif
2309: #endif
2310:   return 0;
2311: }

2313: static PetscErrorCode hypre_array_destroy(void *ptr)
2314: {
2315:   hypre_TFree(ptr, HYPRE_MEMORY_HOST);
2316:   return 0;
2317: }