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 <../src/mat/impls/hypre/mhypre.h>
 11: #include <../src/mat/impls/aij/mpi/mpiaij.h>
 12: #include <../src/vec/vec/impls/hypre/vhyp.h>
 13: #include <HYPRE.h>
 14: #include <HYPRE_utilities.h>
 15: #include <_hypre_parcsr_ls.h>
 16: #include <_hypre_sstruct_ls.h>

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

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

 30: static PetscErrorCode MatHYPRE_IJMatrixPreallocate(Mat A_d, Mat A_o, HYPRE_IJMatrix ij)
 31: {
 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;

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

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

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

129: static PetscErrorCode MatHYPRE_IJMatrixCopy(Mat A, HYPRE_IJMatrix ij)
130: {
131:   PetscErrorCode    ierr;
132:   PetscInt          i,rstart,rend,ncols,nr,nc;
133:   const PetscScalar *values;
134:   const PetscInt    *cols;
135:   PetscBool         flg;

138: #if PETSC_PKG_HYPRE_VERSION_LT(2,19,0)
139:   PetscStackCallStandard(HYPRE_IJMatrixInitialize,(ij));
140: #else
141:   PetscStackCallStandard(HYPRE_IJMatrixInitialize_v2,(ij,HYPRE_MEMORY_HOST));
142: #endif
143:   PetscObjectBaseTypeCompare((PetscObject)A,MATMPIAIJ,&flg);
144:   MatGetSize(A,&nr,&nc);
145:   if (flg && nr == nc) {
146:     MatHYPRE_IJMatrixFastCopy_MPIAIJ(A,ij);
147:     return(0);
148:   }
149:   PetscObjectBaseTypeCompare((PetscObject)A,MATSEQAIJ,&flg);
150:   if (flg) {
151:     MatHYPRE_IJMatrixFastCopy_SeqAIJ(A,ij);
152:     return(0);
153:   }

155:   MatGetOwnershipRange(A,&rstart,&rend);
156:   for (i=rstart; i<rend; i++) {
157:     MatGetRow(A,i,&ncols,&cols,&values);
158:     if (ncols) {
159:       HYPRE_Int nc = (HYPRE_Int)ncols;

161:       if ((PetscInt)nc != ncols) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Hypre overflow! number of columns %D for row %D",ncols,i);
162:       PetscStackCallStandard(HYPRE_IJMatrixSetValues,(ij,1,&nc,(HYPRE_BigInt *)&i,(HYPRE_BigInt *)cols,(HYPRE_Complex *)values));
163:     }
164:     MatRestoreRow(A,i,&ncols,&cols,&values);
165:   }
166:   return(0);
167: }

169: static PetscErrorCode MatHYPRE_IJMatrixFastCopy_SeqAIJ(Mat A, HYPRE_IJMatrix ij)
170: {
171:   PetscErrorCode        ierr;
172:   Mat_SeqAIJ            *pdiag = (Mat_SeqAIJ*)A->data;
173:   HYPRE_Int             type;
174:   hypre_ParCSRMatrix    *par_matrix;
175:   hypre_AuxParCSRMatrix *aux_matrix;
176:   hypre_CSRMatrix       *hdiag;
177:   PetscBool             sameint = (PetscBool)(sizeof(PetscInt) == sizeof(HYPRE_Int));
178:   const PetscScalar     *pa;

181:   PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(ij,&type));
182:   if (type != HYPRE_PARCSR) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported");
183:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(ij,(void**)&par_matrix));
184:   hdiag = hypre_ParCSRMatrixDiag(par_matrix);
185:   /*
186:        this is the Hack part where we monkey directly with the hypre datastructures
187:   */
188:   if (sameint) {
189:     PetscArraycpy(hdiag->i,pdiag->i,A->rmap->n + 1);
190:     PetscArraycpy(hdiag->j,pdiag->j,pdiag->nz);
191:   } else {
192:     PetscInt i;

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

198:   MatSeqAIJGetArrayRead(A,&pa);
199:   PetscArraycpy(hdiag->data,pa,pdiag->nz);
200:   MatSeqAIJRestoreArrayRead(A,&pa);

202:   aux_matrix = (hypre_AuxParCSRMatrix*)hypre_IJMatrixTranslator(ij);
203:   hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 0;
204:   return(0);
205: }

207: static PetscErrorCode MatHYPRE_IJMatrixFastCopy_MPIAIJ(Mat A, HYPRE_IJMatrix ij)
208: {
209:   PetscErrorCode        ierr;
210:   Mat_MPIAIJ            *pA = (Mat_MPIAIJ*)A->data;
211:   Mat_SeqAIJ            *pdiag,*poffd;
212:   PetscInt              i,*garray = pA->garray,*jj,cstart,*pjj;
213:   HYPRE_Int             *hjj,type;
214:   hypre_ParCSRMatrix    *par_matrix;
215:   hypre_AuxParCSRMatrix *aux_matrix;
216:   hypre_CSRMatrix       *hdiag,*hoffd;
217:   PetscBool             sameint = (PetscBool)(sizeof(PetscInt) == sizeof(HYPRE_Int));
218:   const PetscScalar     *pa;

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

226:   PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(ij,&type));
227:   if (type != HYPRE_PARCSR) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported");
228:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(ij,(void**)&par_matrix));
229:   hdiag = hypre_ParCSRMatrixDiag(par_matrix);
230:   hoffd = hypre_ParCSRMatrixOffd(par_matrix);

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

257:   /* need to move the offd column indices (hoffd->j) back to global numbering since hypre is expecting this
258:      If we hacked a hypre a bit more we might be able to avoid this step */
259: #if PETSC_PKG_HYPRE_VERSION_GE(2,16,0)
260:   PetscStackCallStandard(hypre_CSRMatrixBigInitialize,(hoffd));
261:   jj  = (PetscInt*) hoffd->big_j;
262: #else
263:   jj  = (PetscInt*) hoffd->j;
264: #endif
265:   pjj = poffd->j;
266:   for (i=0; i<poffd->nz; i++) jj[i] = garray[pjj[i]];

268:   MatSeqAIJGetArrayRead(pA->B,&pa);
269:   PetscArraycpy(hoffd->data,pa,poffd->nz);
270:   MatSeqAIJRestoreArrayRead(pA->B,&pa);

272:   aux_matrix = (hypre_AuxParCSRMatrix*)hypre_IJMatrixTranslator(ij);
273:   hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 0;
274:   return(0);
275: }

277: static PetscErrorCode MatConvert_HYPRE_IS(Mat A, MatType mtype, MatReuse reuse, Mat* B)
278: {
279:   Mat_HYPRE*             mhA = (Mat_HYPRE*)(A->data);
280:   Mat                    lA;
281:   ISLocalToGlobalMapping rl2g,cl2g;
282:   IS                     is;
283:   hypre_ParCSRMatrix     *hA;
284:   hypre_CSRMatrix        *hdiag,*hoffd;
285:   MPI_Comm               comm;
286:   HYPRE_Complex          *hdd,*hod,*aa;
287:   PetscScalar            *data;
288:   HYPRE_BigInt           *col_map_offd;
289:   HYPRE_Int              *hdi,*hdj,*hoi,*hoj;
290:   PetscInt               *ii,*jj,*iptr,*jptr;
291:   PetscInt               cum,dr,dc,oc,str,stc,nnz,i,jd,jo,M,N;
292:   HYPRE_Int              type;
293:   PetscErrorCode         ierr;

296:   comm = PetscObjectComm((PetscObject)A);
297:   PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(mhA->ij,&type));
298:   if (type != HYPRE_PARCSR) SETERRQ(comm,PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported");
299:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(mhA->ij,(void**)&hA));
300:   M     = hypre_ParCSRMatrixGlobalNumRows(hA);
301:   N     = hypre_ParCSRMatrixGlobalNumCols(hA);
302:   str   = hypre_ParCSRMatrixFirstRowIndex(hA);
303:   stc   = hypre_ParCSRMatrixFirstColDiag(hA);
304:   hdiag = hypre_ParCSRMatrixDiag(hA);
305:   hoffd = hypre_ParCSRMatrixOffd(hA);
306:   dr    = hypre_CSRMatrixNumRows(hdiag);
307:   dc    = hypre_CSRMatrixNumCols(hdiag);
308:   nnz   = hypre_CSRMatrixNumNonzeros(hdiag);
309:   hdi   = hypre_CSRMatrixI(hdiag);
310:   hdj   = hypre_CSRMatrixJ(hdiag);
311:   hdd   = hypre_CSRMatrixData(hdiag);
312:   oc    = hypre_CSRMatrixNumCols(hoffd);
313:   nnz  += hypre_CSRMatrixNumNonzeros(hoffd);
314:   hoi   = hypre_CSRMatrixI(hoffd);
315:   hoj   = hypre_CSRMatrixJ(hoffd);
316:   hod   = hypre_CSRMatrixData(hoffd);
317:   if (reuse != MAT_REUSE_MATRIX) {
318:     PetscInt *aux;

320:     /* generate l2g maps for rows and cols */
321:     ISCreateStride(comm,dr,str,1,&is);
322:     ISLocalToGlobalMappingCreateIS(is,&rl2g);
323:     ISDestroy(&is);
324:     col_map_offd = hypre_ParCSRMatrixColMapOffd(hA);
325:     PetscMalloc1(dc+oc,&aux);
326:     for (i=0; i<dc; i++) aux[i] = i+stc;
327:     for (i=0; i<oc; i++) aux[i+dc] = col_map_offd[i];
328:     ISCreateGeneral(comm,dc+oc,aux,PETSC_OWN_POINTER,&is);
329:     ISLocalToGlobalMappingCreateIS(is,&cl2g);
330:     ISDestroy(&is);
331:     /* create MATIS object */
332:     MatCreate(comm,B);
333:     MatSetSizes(*B,dr,dc,M,N);
334:     MatSetType(*B,MATIS);
335:     MatSetLocalToGlobalMapping(*B,rl2g,cl2g);
336:     ISLocalToGlobalMappingDestroy(&rl2g);
337:     ISLocalToGlobalMappingDestroy(&cl2g);

339:     /* allocate CSR for local matrix */
340:     PetscMalloc1(dr+1,&iptr);
341:     PetscMalloc1(nnz,&jptr);
342:     PetscMalloc1(nnz,&data);
343:   } else {
344:     PetscInt  nr;
345:     PetscBool done;
346:     MatISGetLocalMat(*B,&lA);
347:     MatGetRowIJ(lA,0,PETSC_FALSE,PETSC_FALSE,&nr,(const PetscInt**)&iptr,(const PetscInt**)&jptr,&done);
348:     if (nr != dr) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot reuse mat: invalid number of rows in local mat! %D != %D",nr,dr);
349:     if (iptr[nr] < nnz) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot reuse mat: invalid number of nonzeros in local mat! reuse %D requested %D",iptr[nr],nnz);
350:     MatSeqAIJGetArray(lA,&data);
351:   }
352:   /* merge local matrices */
353:   ii  = iptr;
354:   jj  = jptr;
355:   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++ */
356:   *ii = *(hdi++) + *(hoi++);
357:   for (jd=0,jo=0,cum=0; *ii<nnz; cum++) {
358:     PetscScalar *aold = (PetscScalar*)aa;
359:     PetscInt    *jold = jj,nc = jd+jo;
360:     for (; jd<*hdi; jd++) { *jj++ = *hdj++;      *aa++ = *hdd++; }
361:     for (; jo<*hoi; jo++) { *jj++ = *hoj++ + dc; *aa++ = *hod++; }
362:     *(++ii) = *(hdi++) + *(hoi++);
363:     PetscSortIntWithScalarArray(jd+jo-nc,jold,aold);
364:   }
365:   for (; cum<dr; cum++) *(++ii) = nnz;
366:   if (reuse != MAT_REUSE_MATRIX) {
367:     Mat_SeqAIJ* a;

369:     MatCreateSeqAIJWithArrays(PETSC_COMM_SELF,dr,dc+oc,iptr,jptr,data,&lA);
370:     MatISSetLocalMat(*B,lA);
371:     /* hack SeqAIJ */
372:     a          = (Mat_SeqAIJ*)(lA->data);
373:     a->free_a  = PETSC_TRUE;
374:     a->free_ij = PETSC_TRUE;
375:     MatDestroy(&lA);
376:   }
377:   MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
378:   MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
379:   if (reuse == MAT_INPLACE_MATRIX) {
380:     MatHeaderReplace(A,B);
381:   }
382:   return(0);
383: }

385: PETSC_INTERN PetscErrorCode MatConvert_AIJ_HYPRE(Mat A, MatType type, MatReuse reuse, Mat *B)
386: {
387:   Mat            M = NULL;
388:   Mat_HYPRE      *hB;
389:   MPI_Comm       comm = PetscObjectComm((PetscObject)A);

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

421: static PetscErrorCode MatConvert_HYPRE_AIJ(Mat A, MatType mtype, MatReuse reuse, Mat *B)
422: {
423:   Mat_HYPRE          *hA = (Mat_HYPRE*)A->data;
424:   hypre_ParCSRMatrix *parcsr;
425:   hypre_CSRMatrix    *hdiag,*hoffd;
426:   MPI_Comm           comm;
427:   PetscScalar        *da,*oa,*aptr;
428:   PetscInt           *dii,*djj,*oii,*ojj,*iptr;
429:   PetscInt           i,dnnz,onnz,m,n;
430:   HYPRE_Int          type;
431:   PetscMPIInt        size;
432:   PetscBool          sameint = (PetscBool)(sizeof(PetscInt) == sizeof(HYPRE_Int));
433:   PetscErrorCode     ierr;

436:   comm = PetscObjectComm((PetscObject)A);
437:   PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(hA->ij,&type));
438:   if (type != HYPRE_PARCSR) SETERRQ(comm,PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported");
439:   if (reuse == MAT_REUSE_MATRIX) {
440:     PetscBool ismpiaij,isseqaij;
441:     PetscObjectBaseTypeCompare((PetscObject)*B,MATMPIAIJ,&ismpiaij);
442:     PetscObjectBaseTypeCompare((PetscObject)*B,MATSEQAIJ,&isseqaij);
443:     if (!ismpiaij && !isseqaij) SETERRQ(comm,PETSC_ERR_SUP,"Only MATMPIAIJ or MATSEQAIJ are supported");
444:   }
445: #if defined(PETSC_HAVE_HYPRE_DEVICE)
446:   if (HYPRE_MEMORY_DEVICE == hypre_IJMatrixMemoryLocation(hA->ij)) SETERRQ(comm,PETSC_ERR_SUP,"Not yet implemented");
447: #endif
448:   MPI_Comm_size(comm,&size);

450:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hA->ij,(void**)&parcsr));
451:   hdiag = hypre_ParCSRMatrixDiag(parcsr);
452:   hoffd = hypre_ParCSRMatrixOffd(parcsr);
453:   m     = hypre_CSRMatrixNumRows(hdiag);
454:   n     = hypre_CSRMatrixNumCols(hdiag);
455:   dnnz  = hypre_CSRMatrixNumNonzeros(hdiag);
456:   onnz  = hypre_CSRMatrixNumNonzeros(hoffd);
457:   if (reuse == MAT_INITIAL_MATRIX) {
458:     PetscMalloc1(m+1,&dii);
459:     PetscMalloc1(dnnz,&djj);
460:     PetscMalloc1(dnnz,&da);
461:   } else if (reuse == MAT_REUSE_MATRIX) {
462:     PetscInt  nr;
463:     PetscBool done;
464:     if (size > 1) {
465:       Mat_MPIAIJ *b = (Mat_MPIAIJ*)((*B)->data);

467:       MatGetRowIJ(b->A,0,PETSC_FALSE,PETSC_FALSE,&nr,(const PetscInt**)&dii,(const PetscInt**)&djj,&done);
468:       if (nr != m) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot reuse mat: invalid number of local rows in diag part! %D != %D",nr,m);
469:       if (dii[nr] < dnnz) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot reuse mat: invalid number of nonzeros in diag part! reuse %D hypre %D",dii[nr],dnnz);
470:       MatSeqAIJGetArray(b->A,&da);
471:     } else {
472:       MatGetRowIJ(*B,0,PETSC_FALSE,PETSC_FALSE,&nr,(const PetscInt**)&dii,(const PetscInt**)&djj,&done);
473:       if (nr != m) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot reuse mat: invalid number of local rows! %D != %D",nr,m);
474:       if (dii[nr] < dnnz) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot reuse mat: invalid number of nonzeros! reuse %D hypre %D",dii[nr],dnnz);
475:       MatSeqAIJGetArray(*B,&da);
476:     }
477:   } else { /* MAT_INPLACE_MATRIX */
478:     if (!sameint) {
479:       PetscMalloc1(m+1,&dii);
480:       PetscMalloc1(dnnz,&djj);
481:     } else {
482:       dii = (PetscInt*)hypre_CSRMatrixI(hdiag);
483:       djj = (PetscInt*)hypre_CSRMatrixJ(hdiag);
484:     }
485:     da = (PetscScalar*)hypre_CSRMatrixData(hdiag);
486:   }

488:   if (!sameint) {
489:     for (i=0;i<m+1;i++)  dii[i] = (PetscInt)(hypre_CSRMatrixI(hdiag)[i]);
490:     for (i=0;i<dnnz;i++) djj[i] = (PetscInt)(hypre_CSRMatrixJ(hdiag)[i]);
491:   } else {
492:     PetscArraycpy(dii,hypre_CSRMatrixI(hdiag),m+1);
493:     PetscArraycpy(djj,hypre_CSRMatrixJ(hdiag),dnnz);
494:   }
495:   PetscArraycpy(da,hypre_CSRMatrixData(hdiag),dnnz);
496:   iptr = djj;
497:   aptr = da;
498:   for (i=0; i<m; i++) {
499:     PetscInt nc = dii[i+1]-dii[i];
500:     PetscSortIntWithScalarArray(nc,iptr,aptr);
501:     iptr += nc;
502:     aptr += nc;
503:   }
504:   if (size > 1) {
505:     HYPRE_BigInt *coffd;
506:     HYPRE_Int    *offdj;

508:     if (reuse == MAT_INITIAL_MATRIX) {
509:       PetscMalloc1(m+1,&oii);
510:       PetscMalloc1(onnz,&ojj);
511:       PetscMalloc1(onnz,&oa);
512:     } else if (reuse == MAT_REUSE_MATRIX) {
513:       Mat_MPIAIJ *b = (Mat_MPIAIJ*)((*B)->data);
514:       PetscInt   nr,hr = hypre_CSRMatrixNumRows(hoffd);
515:       PetscBool  done;

517:       MatGetRowIJ(b->B,0,PETSC_FALSE,PETSC_FALSE,&nr,(const PetscInt**)&oii,(const PetscInt**)&ojj,&done);
518:       if (nr != hr) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot reuse mat: invalid number of local rows in offdiag part! %D != %D",nr,hr);
519:       if (oii[nr] < onnz) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot reuse mat: invalid number of nonzeros in offdiag part! reuse %D hypre %D",oii[nr],onnz);
520:       MatSeqAIJGetArray(b->B,&oa);
521:     } else { /* MAT_INPLACE_MATRIX */
522:       if (!sameint) {
523:         PetscMalloc1(m+1,&oii);
524:         PetscMalloc1(onnz,&ojj);
525:       } else {
526:         oii = (PetscInt*)hypre_CSRMatrixI(hoffd);
527:         ojj = (PetscInt*)hypre_CSRMatrixJ(hoffd);
528:       }
529:       oa = (PetscScalar*)hypre_CSRMatrixData(hoffd);
530:     }
531:     if (!sameint) {
532:       for (i=0;i<m+1;i++) oii[i] = (PetscInt)(hypre_CSRMatrixI(hoffd)[i]);
533:     } else {
534:       PetscArraycpy(oii,hypre_CSRMatrixI(hoffd),m+1);
535:     }
536:     offdj = hypre_CSRMatrixJ(hoffd);
537:     coffd = hypre_ParCSRMatrixColMapOffd(parcsr);
538:     for (i=0; i<onnz; i++) ojj[i] = coffd[offdj[i]];
539:     PetscArraycpy(oa,hypre_CSRMatrixData(hoffd),onnz);
540:     iptr = ojj;
541:     aptr = oa;
542:     for (i=0; i<m; i++) {
543:        PetscInt nc = oii[i+1]-oii[i];
544:        PetscSortIntWithScalarArray(nc,iptr,aptr);
545:        iptr += nc;
546:        aptr += nc;
547:     }
548:     if (reuse == MAT_INITIAL_MATRIX) {
549:       Mat_MPIAIJ *b;
550:       Mat_SeqAIJ *d,*o;

552:       MatCreateMPIAIJWithSplitArrays(comm,m,n,PETSC_DECIDE,PETSC_DECIDE,dii,djj,da,oii,ojj,oa,B);
553:       /* hack MPIAIJ */
554:       b          = (Mat_MPIAIJ*)((*B)->data);
555:       d          = (Mat_SeqAIJ*)b->A->data;
556:       o          = (Mat_SeqAIJ*)b->B->data;
557:       d->free_a  = PETSC_TRUE;
558:       d->free_ij = PETSC_TRUE;
559:       o->free_a  = PETSC_TRUE;
560:       o->free_ij = PETSC_TRUE;
561:     } else if (reuse == MAT_INPLACE_MATRIX) {
562:       Mat T;

564:       MatCreateMPIAIJWithSplitArrays(comm,m,n,PETSC_DECIDE,PETSC_DECIDE,dii,djj,da,oii,ojj,oa,&T);
565:       if (sameint) { /* ownership of CSR pointers is transferred to PETSc */
566:         hypre_CSRMatrixI(hdiag) = NULL;
567:         hypre_CSRMatrixJ(hdiag) = NULL;
568:         hypre_CSRMatrixI(hoffd) = NULL;
569:         hypre_CSRMatrixJ(hoffd) = NULL;
570:       } else { /* Hack MPIAIJ -> free ij but not a */
571:         Mat_MPIAIJ *b = (Mat_MPIAIJ*)(T->data);
572:         Mat_SeqAIJ *d = (Mat_SeqAIJ*)(b->A->data);
573:         Mat_SeqAIJ *o = (Mat_SeqAIJ*)(b->B->data);

575:         d->free_ij = PETSC_TRUE;
576:         o->free_ij = PETSC_TRUE;
577:       }
578:       hypre_CSRMatrixData(hdiag) = NULL;
579:       hypre_CSRMatrixData(hoffd) = NULL;
580:       MatHeaderReplace(A,&T);
581:     }
582:   } else {
583:     oii  = NULL;
584:     ojj  = NULL;
585:     oa   = NULL;
586:     if (reuse == MAT_INITIAL_MATRIX) {
587:       Mat_SeqAIJ* b;

589:       MatCreateSeqAIJWithArrays(comm,m,n,dii,djj,da,B);
590:       /* hack SeqAIJ */
591:       b          = (Mat_SeqAIJ*)((*B)->data);
592:       b->free_a  = PETSC_TRUE;
593:       b->free_ij = PETSC_TRUE;
594:     } else if (reuse == MAT_INPLACE_MATRIX) {
595:       Mat T;

597:       MatCreateSeqAIJWithArrays(comm,m,n,dii,djj,da,&T);
598:       if (sameint) { /* ownership of CSR pointers is transferred to PETSc */
599:         hypre_CSRMatrixI(hdiag) = NULL;
600:         hypre_CSRMatrixJ(hdiag) = NULL;
601:       } else { /* free ij but not a */
602:         Mat_SeqAIJ* b = (Mat_SeqAIJ*)(T->data);

604:         b->free_ij = PETSC_TRUE;
605:       }
606:       hypre_CSRMatrixData(hdiag) = NULL;
607:       MatHeaderReplace(A,&T);
608:     }
609:   }

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

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

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

653:   PetscObjectBaseTypeCompare((PetscObject)A,MATMPIAIJ,&ismpiaij);
654:   PetscObjectBaseTypeCompare((PetscObject)A,MATSEQAIJ,&isseqaij);
655:   if (!ismpiaij && !isseqaij) SETERRQ1(comm,PETSC_ERR_SUP,"Unsupported type %s",((PetscObject)A)->type_name);
656:   if (ismpiaij) {
657:     Mat_MPIAIJ *a = (Mat_MPIAIJ*)(A->data);

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

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

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

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

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

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

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

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

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

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

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

860:   MatAIJGetParCSR_Private(A,&hA);
861:   MatAIJGetParCSR_Private(P,&hP);
862:   MatHYPRE_ParCSR_RAP(hP,hA,hP,&hPtAP);
863:   MatCreateFromParCSR(hPtAP,MATAIJ,PETSC_OWN_POINTER,&B);

865:   MatHeaderMerge(C,&B);
866:   C->product = product;

868:   MatAIJRestoreParCSR_Private(A,&hA);
869:   MatAIJRestoreParCSR_Private(P,&hP);
870:   return(0);
871: }

873: PETSC_INTERN PetscErrorCode MatPtAPSymbolic_AIJ_AIJ_wHYPRE(Mat A,Mat P,PetscReal fill,Mat C)
874: {

878:   MatSetType(C,MATAIJ);
879:   C->ops->ptapnumeric    = MatPtAPNumeric_AIJ_AIJ_wHYPRE;
880:   C->ops->productnumeric = MatProductNumeric_PtAP;
881:   return(0);
882: }

884: static PetscErrorCode MatPtAPNumeric_AIJ_HYPRE(Mat A,Mat P,Mat C)
885: {
886:   Mat                B;
887:   Mat_HYPRE          *hP;
888:   hypre_ParCSRMatrix *hA = NULL,*Pparcsr,*ptapparcsr = NULL;
889:   HYPRE_Int          type;
890:   MPI_Comm           comm = PetscObjectComm((PetscObject)A);
891:   PetscBool          ishypre;
892:   PetscErrorCode     ierr;

895:   PetscObjectTypeCompare((PetscObject)P,MATHYPRE,&ishypre);
896:   if (!ishypre) SETERRQ1(comm,PETSC_ERR_USER,"P should be of type %s",MATHYPRE);
897:   hP = (Mat_HYPRE*)P->data;
898:   PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(hP->ij,&type));
899:   if (type != HYPRE_PARCSR) SETERRQ(comm,PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported");
900:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hP->ij,(void**)&Pparcsr));

902:   MatAIJGetParCSR_Private(A,&hA);
903:   MatHYPRE_ParCSR_RAP(Pparcsr,hA,Pparcsr,&ptapparcsr);
904:   MatAIJRestoreParCSR_Private(A,&hA);

906:   /* create temporary matrix and merge to C */
907:   MatCreateFromParCSR(ptapparcsr,((PetscObject)C)->type_name,PETSC_OWN_POINTER,&B);
908:   MatHeaderMerge(C,&B);
909:   return(0);
910: }

912: static PetscErrorCode MatPtAPNumeric_HYPRE_HYPRE(Mat A,Mat P,Mat C)
913: {
914:   Mat                B;
915:   hypre_ParCSRMatrix *Aparcsr,*Pparcsr,*ptapparcsr = NULL;
916:   Mat_HYPRE          *hA,*hP;
917:   PetscBool          ishypre;
918:   HYPRE_Int          type;
919:   PetscErrorCode     ierr;

922:   PetscObjectTypeCompare((PetscObject)P,MATHYPRE,&ishypre);
923:   if (!ishypre) SETERRQ1(PetscObjectComm((PetscObject)P),PETSC_ERR_USER,"P should be of type %s",MATHYPRE);
924:   PetscObjectTypeCompare((PetscObject)A,MATHYPRE,&ishypre);
925:   if (!ishypre) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"A should be of type %s",MATHYPRE);
926:   hA = (Mat_HYPRE*)A->data;
927:   hP = (Mat_HYPRE*)P->data;
928:   PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(hA->ij,&type));
929:   if (type != HYPRE_PARCSR) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported");
930:   PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(hP->ij,&type));
931:   if (type != HYPRE_PARCSR) SETERRQ(PetscObjectComm((PetscObject)P),PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported");
932:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hA->ij,(void**)&Aparcsr));
933:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hP->ij,(void**)&Pparcsr));
934:   MatHYPRE_ParCSR_RAP(Pparcsr,Aparcsr,Pparcsr,&ptapparcsr);
935:   MatCreateFromParCSR(ptapparcsr,MATHYPRE,PETSC_OWN_POINTER,&B);
936:   MatHeaderMerge(C,&B);
937:   return(0);
938: }

940: /* calls hypre_ParMatmul
941:    hypre_ParMatMul uses hypre_ParMatrixCreate with the communicator of hA
942:    hypre_ParMatrixCreate does not duplicate the communicator
943:    It looks like we don't need to have the diagonal entries ordered first */
944: static PetscErrorCode MatHYPRE_ParCSR_MatMatMult(hypre_ParCSRMatrix *hA, hypre_ParCSRMatrix *hB, hypre_ParCSRMatrix **hAB)
945: {
947:   /* can be replaced by version test later */
948: #if defined(PETSC_HAVE_HYPRE_DEVICE)
949:   PetscStackPush("hypre_ParCSRMatMat");
950:   *hAB = hypre_ParCSRMatMat(hA,hB);
951: #else
952:   PetscStackPush("hypre_ParMatmul");
953:   *hAB = hypre_ParMatmul(hA,hB);
954: #endif
955:   PetscStackPop;
956:   return(0);
957: }

959: static PetscErrorCode MatMatMultNumeric_AIJ_AIJ_wHYPRE(Mat A,Mat B,Mat C)
960: {
961:   Mat                D;
962:   hypre_ParCSRMatrix *hA,*hB,*hAB = NULL;
963:   PetscErrorCode     ierr;
964:   Mat_Product        *product=C->product;

967:   MatAIJGetParCSR_Private(A,&hA);
968:   MatAIJGetParCSR_Private(B,&hB);
969:   MatHYPRE_ParCSR_MatMatMult(hA,hB,&hAB);
970:   MatCreateFromParCSR(hAB,MATAIJ,PETSC_OWN_POINTER,&D);

972:   MatHeaderMerge(C,&D);
973:   C->product = product;

975:   MatAIJRestoreParCSR_Private(A,&hA);
976:   MatAIJRestoreParCSR_Private(B,&hB);
977:   return(0);
978: }

980: PETSC_INTERN PetscErrorCode MatMatMultSymbolic_AIJ_AIJ_wHYPRE(Mat A,Mat B,PetscReal fill,Mat C)
981: {

985:   MatSetType(C,MATAIJ);
986:   C->ops->matmultnumeric = MatMatMultNumeric_AIJ_AIJ_wHYPRE;
987:   C->ops->productnumeric = MatProductNumeric_AB;
988:   return(0);
989: }

991: static PetscErrorCode MatMatMultNumeric_HYPRE_HYPRE(Mat A,Mat B,Mat C)
992: {
993:   Mat                D;
994:   hypre_ParCSRMatrix *Aparcsr,*Bparcsr,*ABparcsr = NULL;
995:   Mat_HYPRE          *hA,*hB;
996:   PetscBool          ishypre;
997:   HYPRE_Int          type;
998:   PetscErrorCode     ierr;
999:   Mat_Product        *product;

1002:   PetscObjectTypeCompare((PetscObject)B,MATHYPRE,&ishypre);
1003:   if (!ishypre) SETERRQ1(PetscObjectComm((PetscObject)B),PETSC_ERR_USER,"B should be of type %s",MATHYPRE);
1004:   PetscObjectTypeCompare((PetscObject)A,MATHYPRE,&ishypre);
1005:   if (!ishypre) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"A should be of type %s",MATHYPRE);
1006:   hA = (Mat_HYPRE*)A->data;
1007:   hB = (Mat_HYPRE*)B->data;
1008:   PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(hA->ij,&type));
1009:   if (type != HYPRE_PARCSR) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported");
1010:   PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(hB->ij,&type));
1011:   if (type != HYPRE_PARCSR) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported");
1012:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hA->ij,(void**)&Aparcsr));
1013:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hB->ij,(void**)&Bparcsr));
1014:   MatHYPRE_ParCSR_MatMatMult(Aparcsr,Bparcsr,&ABparcsr);
1015:   MatCreateFromParCSR(ABparcsr,MATHYPRE,PETSC_OWN_POINTER,&D);

1017:   /* need to use HeaderReplace because HeaderMerge messes up with the communicator */
1018:   product    = C->product;  /* save it from MatHeaderReplace() */
1019:   C->product = NULL;
1020:   MatHeaderReplace(C,&D);
1021:   C->product = product;
1022:   C->ops->matmultnumeric = MatMatMultNumeric_HYPRE_HYPRE;
1023:   C->ops->productnumeric = MatProductNumeric_AB;
1024:   return(0);
1025: }

1027: PETSC_INTERN PetscErrorCode MatTransposeMatMatMultNumeric_AIJ_AIJ_AIJ_wHYPRE(Mat A,Mat B,Mat C,Mat D)
1028: {
1029:   Mat                E;
1030:   hypre_ParCSRMatrix *hA,*hB,*hC,*hABC = NULL;
1031:   PetscErrorCode     ierr;

1034:   MatAIJGetParCSR_Private(A,&hA);
1035:   MatAIJGetParCSR_Private(B,&hB);
1036:   MatAIJGetParCSR_Private(C,&hC);
1037:   MatHYPRE_ParCSR_RAP(hA,hB,hC,&hABC);
1038:   MatCreateFromParCSR(hABC,MATAIJ,PETSC_OWN_POINTER,&E);
1039:   MatHeaderMerge(D,&E);
1040:   MatAIJRestoreParCSR_Private(A,&hA);
1041:   MatAIJRestoreParCSR_Private(B,&hB);
1042:   MatAIJRestoreParCSR_Private(C,&hC);
1043:   return(0);
1044: }

1046: PETSC_INTERN PetscErrorCode MatTransposeMatMatMultSymbolic_AIJ_AIJ_AIJ_wHYPRE(Mat A,Mat B,Mat C,PetscReal fill,Mat D)
1047: {

1051:   MatSetType(D,MATAIJ);
1052:   return(0);
1053: }

1055: /* ---------------------------------------------------- */
1056: static PetscErrorCode MatProductSymbolic_AB_HYPRE(Mat C)
1057: {
1059:   C->ops->productnumeric = MatProductNumeric_AB;
1060:   return(0);
1061: }

1063: static PetscErrorCode MatProductSetFromOptions_HYPRE_AB(Mat C)
1064: {
1066:   Mat_Product    *product = C->product;
1067:   PetscBool      Ahypre;

1070:   PetscObjectTypeCompare((PetscObject)product->A,MATHYPRE,&Ahypre);
1071:   if (Ahypre) { /* A is a Hypre matrix */
1072:     MatSetType(C,MATHYPRE);
1073:     C->ops->productsymbolic = MatProductSymbolic_AB_HYPRE;
1074:     C->ops->matmultnumeric  = MatMatMultNumeric_HYPRE_HYPRE;
1075:     return(0);
1076:   }
1077:   return(0);
1078: }

1080: static PetscErrorCode MatProductSymbolic_PtAP_HYPRE(Mat C)
1081: {
1083:   C->ops->productnumeric = MatProductNumeric_PtAP;
1084:   return(0);
1085: }

1087: static PetscErrorCode MatProductSetFromOptions_HYPRE_PtAP(Mat C)
1088: {
1090:   Mat_Product    *product = C->product;
1091:   PetscBool      flg;
1092:   PetscInt       type = 0;
1093:   const char     *outTypes[4] = {"aij","seqaij","mpiaij","hypre"};
1094:   PetscInt       ntype = 4;
1095:   Mat            A = product->A;
1096:   PetscBool      Ahypre;

1099:   PetscObjectTypeCompare((PetscObject)A,MATHYPRE,&Ahypre);
1100:   if (Ahypre) { /* A is a Hypre matrix */
1101:     MatSetType(C,MATHYPRE);
1102:     C->ops->productsymbolic = MatProductSymbolic_PtAP_HYPRE;
1103:     C->ops->ptapnumeric     = MatPtAPNumeric_HYPRE_HYPRE;
1104:     return(0);
1105:   }

1107:   /* A is AIJ, P is Hypre, C = PtAP can be either AIJ or Hypre format */
1108:   /* Get runtime option */
1109:   if (product->api_user) {
1110:     PetscOptionsBegin(PetscObjectComm((PetscObject)C),((PetscObject)C)->prefix,"MatPtAP_HYPRE","Mat");
1111:     PetscOptionsEList("-matptap_hypre_outtype","MatPtAP outtype","MatPtAP outtype",outTypes,ntype,outTypes[type],&type,&flg);
1112:     PetscOptionsEnd();
1113:   } else {
1114:     PetscOptionsBegin(PetscObjectComm((PetscObject)C),((PetscObject)C)->prefix,"MatProduct_PtAP_HYPRE","Mat");
1115:     PetscOptionsEList("-matproduct_ptap_hypre_outtype","MatProduct_PtAP outtype","MatProduct_PtAP",outTypes,ntype,outTypes[type],&type,&flg);
1116:     PetscOptionsEnd();
1117:   }

1119:   if (type == 0 || type == 1 || type == 2) {
1120:     MatSetType(C,MATAIJ);
1121:   } else if (type == 3) {
1122:     MatSetType(C,MATHYPRE);
1123:   } else SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_SUP,"MatPtAP outtype is not supported");
1124:   C->ops->productsymbolic = MatProductSymbolic_PtAP_HYPRE;
1125:   C->ops->ptapnumeric     = MatPtAPNumeric_AIJ_HYPRE;
1126:   return(0);
1127: }

1129: static PetscErrorCode MatProductSetFromOptions_HYPRE(Mat C)
1130: {
1132:   Mat_Product    *product = C->product;

1135:   switch (product->type) {
1136:   case MATPRODUCT_AB:
1137:     MatProductSetFromOptions_HYPRE_AB(C);
1138:     break;
1139:   case MATPRODUCT_PtAP:
1140:     MatProductSetFromOptions_HYPRE_PtAP(C);
1141:     break;
1142:   default:
1143:     break;
1144:   }
1145:   return(0);
1146: }

1148: /* -------------------------------------------------------- */

1150: static PetscErrorCode MatMultTranspose_HYPRE(Mat A, Vec x, Vec y)
1151: {

1155:   MatHYPRE_MultKernel_Private(A,1.0,x,0.0,y,PETSC_TRUE);
1156:   return(0);
1157: }

1159: static PetscErrorCode MatMult_HYPRE(Mat A, Vec x, Vec y)
1160: {

1164:   MatHYPRE_MultKernel_Private(A,1.0,x,0.0,y,PETSC_FALSE);
1165:   return(0);
1166: }

1168: static PetscErrorCode MatMultAdd_HYPRE(Mat A, Vec x, Vec y, Vec z)
1169: {

1173:   if (y != z) {
1174:     VecCopy(y,z);
1175:   }
1176:   MatHYPRE_MultKernel_Private(A,1.0,x,1.0,z,PETSC_FALSE);
1177:   return(0);
1178: }

1180: static PetscErrorCode MatMultTransposeAdd_HYPRE(Mat A, Vec x, Vec y, Vec z)
1181: {

1185:   if (y != z) {
1186:     VecCopy(y,z);
1187:   }
1188:   MatHYPRE_MultKernel_Private(A,1.0,x,1.0,z,PETSC_TRUE);
1189:   return(0);
1190: }

1192: /* y = a * A * x + b * y or y = a * A^t * x + b * y depending on trans */
1193: static PetscErrorCode MatHYPRE_MultKernel_Private(Mat A, HYPRE_Complex a, Vec x, HYPRE_Complex b, Vec y, PetscBool trans)
1194: {
1195:   Mat_HYPRE          *hA = (Mat_HYPRE*)A->data;
1196:   hypre_ParCSRMatrix *parcsr;
1197:   hypre_ParVector    *hx,*hy;
1198:   PetscErrorCode     ierr;

1201:   if (trans) {
1202:     VecHYPRE_IJVectorPushVecRead(hA->b,x);
1203:     if (b != 0.0) { VecHYPRE_IJVectorPushVec(hA->x,y); }
1204:     else { VecHYPRE_IJVectorPushVecWrite(hA->x,y); }
1205:     PetscStackCallStandard(HYPRE_IJVectorGetObject,(hA->b->ij,(void**)&hx));
1206:     PetscStackCallStandard(HYPRE_IJVectorGetObject,(hA->x->ij,(void**)&hy));
1207:   } else {
1208:     VecHYPRE_IJVectorPushVecRead(hA->x,x);
1209:     if (b != 0.0) { VecHYPRE_IJVectorPushVec(hA->b,y); }
1210:     else { VecHYPRE_IJVectorPushVecWrite(hA->b,y); }
1211:     PetscStackCallStandard(HYPRE_IJVectorGetObject,(hA->x->ij,(void**)&hx));
1212:     PetscStackCallStandard(HYPRE_IJVectorGetObject,(hA->b->ij,(void**)&hy));
1213:   }
1214:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hA->ij,(void**)&parcsr));
1215:   if (trans) {
1216:     PetscStackCallStandard(hypre_ParCSRMatrixMatvecT,(a,parcsr,hx,b,hy));
1217:   } else {
1218:     PetscStackCallStandard(hypre_ParCSRMatrixMatvec,(a,parcsr,hx,b,hy));
1219:   }
1220:   VecHYPRE_IJVectorPopVec(hA->x);
1221:   VecHYPRE_IJVectorPopVec(hA->b);
1222:   return(0);
1223: }

1225: static PetscErrorCode MatDestroy_HYPRE(Mat A)
1226: {
1227:   Mat_HYPRE      *hA = (Mat_HYPRE*)A->data;

1231:   VecHYPRE_IJVectorDestroy(&hA->x);
1232:   VecHYPRE_IJVectorDestroy(&hA->b);
1233:   if (hA->ij) {
1234:     if (!hA->inner_free) hypre_IJMatrixObject(hA->ij) = NULL;
1235:     PetscStackCallStandard(HYPRE_IJMatrixDestroy,(hA->ij));
1236:   }
1237:   if (hA->comm) {MPI_Comm_free(&hA->comm);}

1239:   MatStashDestroy_Private(&A->stash);

1241:   PetscFree(hA->array);

1243:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_hypre_aij_C",NULL);
1244:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_hypre_is_C",NULL);
1245:   PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_seqaij_hypre_C",NULL);
1246:   PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_mpiaij_hypre_C",NULL);
1247:   PetscObjectComposeFunction((PetscObject)A,"MatHYPRESetPreallocation_C",NULL);
1248:   PetscObjectComposeFunction((PetscObject)A,"MatHYPREGetParCSR_C",NULL);
1249:   PetscFree(A->data);
1250:   return(0);
1251: }

1253: static PetscErrorCode MatSetUp_HYPRE(Mat A)
1254: {

1258:   MatHYPRESetPreallocation(A,PETSC_DEFAULT,NULL,PETSC_DEFAULT,NULL);
1259:   return(0);
1260: }

1262: //TODO FIX hypre_CSRMatrixMatvecOutOfPlace
1263: #if defined(PETSC_HAVE_HYPRE_DEVICE)
1264: static PetscErrorCode MatBindToCPU_HYPRE(Mat A, PetscBool bind)
1265: {
1266:   Mat_HYPRE            *hA = (Mat_HYPRE*)A->data;
1267:   HYPRE_MemoryLocation hmem = bind ? HYPRE_MEMORY_HOST : HYPRE_MEMORY_DEVICE;
1268:   PetscErrorCode       ierr;

1271:   A->boundtocpu = bind;
1272:   if (hypre_IJMatrixAssembleFlag(hA->ij) && hmem != hypre_IJMatrixMemoryLocation(hA->ij)) {
1273:     hypre_ParCSRMatrix *parcsr;
1274:     PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hA->ij,(void**)&parcsr));
1275:     PetscStackCallStandard(hypre_ParCSRMatrixMigrate,(parcsr, hmem));
1276:   }
1277:   if (hA->x) { VecHYPRE_IJBindToCPU(hA->x,bind); }
1278:   if (hA->b) { VecHYPRE_IJBindToCPU(hA->b,bind); }
1279:   return(0);
1280: }
1281: #endif

1283: static PetscErrorCode MatAssemblyEnd_HYPRE(Mat A, MatAssemblyType mode)
1284: {
1285:   Mat_HYPRE          *hA = (Mat_HYPRE*)A->data;
1286:   PetscMPIInt        n;
1287:   PetscInt           i,j,rstart,ncols,flg;
1288:   PetscInt           *row,*col;
1289:   PetscScalar        *val;
1290:   PetscErrorCode     ierr;

1293:   if (mode == MAT_FLUSH_ASSEMBLY) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"MAT_FLUSH_ASSEMBLY currently not supported with MATHYPRE");

1295:   if (!A->nooffprocentries) {
1296:     while (1) {
1297:       MatStashScatterGetMesg_Private(&A->stash,&n,&row,&col,&val,&flg);
1298:       if (!flg) break;

1300:       for (i=0; i<n;) {
1301:         /* Now identify the consecutive vals belonging to the same row */
1302:         for (j=i,rstart=row[j]; j<n; j++) {
1303:           if (row[j] != rstart) break;
1304:         }
1305:         if (j < n) ncols = j-i;
1306:         else       ncols = n-i;
1307:         /* Now assemble all these values with a single function call */
1308:         MatSetValues_HYPRE(A,1,row+i,ncols,col+i,val+i,A->insertmode);

1310:         i = j;
1311:       }
1312:     }
1313:     MatStashScatterEnd_Private(&A->stash);
1314:   }

1316:   PetscStackCallStandard(HYPRE_IJMatrixAssemble,(hA->ij));
1317:   /* The assembly routine destroys the aux_matrix, we recreate it here by calling HYPRE_IJMatrixInitialize */
1318:   /* 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 */
1319:   if (!hA->sorted_full) {
1320:     hypre_AuxParCSRMatrix *aux_matrix;

1322:     /* call destroy just to make sure we do not leak anything */
1323:     aux_matrix = (hypre_AuxParCSRMatrix*)hypre_IJMatrixTranslator(hA->ij);
1324:     PetscStackCallStandard(hypre_AuxParCSRMatrixDestroy,(aux_matrix));
1325:     hypre_IJMatrixTranslator(hA->ij) = NULL;

1327:     /* Initialize with assembled flag -> it only recreates the aux_par_matrix */
1328:     PetscStackCallStandard(HYPRE_IJMatrixInitialize,(hA->ij));
1329:     aux_matrix = (hypre_AuxParCSRMatrix*)hypre_IJMatrixTranslator(hA->ij);
1330:     if (aux_matrix) {
1331:       hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 1; /* see comment in MatHYPRESetPreallocation_HYPRE */
1332: #if PETSC_PKG_HYPRE_VERSION_LT(2,19,0)
1333:       PetscStackCallStandard(hypre_AuxParCSRMatrixInitialize,(aux_matrix));
1334: #else
1335:       PetscStackCallStandard(hypre_AuxParCSRMatrixInitialize_v2,(aux_matrix,HYPRE_MEMORY_HOST));
1336: #endif
1337:     }
1338:   }
1339:   {
1340:     hypre_ParCSRMatrix *parcsr;

1342:     PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hA->ij,(void**)&parcsr));
1343:     if (!hypre_ParCSRMatrixCommPkg(parcsr)) PetscStackCallStandard(hypre_MatvecCommPkgCreate,(parcsr));
1344:   }
1345:   if (!hA->x) { VecHYPRE_IJVectorCreate(A->cmap,&hA->x); }
1346:   if (!hA->b) { VecHYPRE_IJVectorCreate(A->rmap,&hA->b); }
1347: #if defined(PETSC_HAVE_HYPRE_DEVICE)
1348:   MatBindToCPU_HYPRE(A,A->boundtocpu);
1349: #endif
1350:   return(0);
1351: }

1353: static PetscErrorCode MatGetArray_HYPRE(Mat A, PetscInt size, void **array)
1354: {
1355:   Mat_HYPRE          *hA = (Mat_HYPRE*)A->data;
1356:   PetscErrorCode     ierr;

1359:   if (!hA->available) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_NULL,"Temporary space is in use");

1361:   if (hA->size >= size) {
1362:     *array = hA->array;
1363:   } else {
1364:     PetscFree(hA->array);
1365:     hA->size = size;
1366:     PetscMalloc(hA->size,&hA->array);
1367:     *array = hA->array;
1368:   }

1370:   hA->available = PETSC_FALSE;
1371:   return(0);
1372: }

1374: static PetscErrorCode MatRestoreArray_HYPRE(Mat A, void **array)
1375: {
1376:   Mat_HYPRE          *hA = (Mat_HYPRE*)A->data;

1379:   *array = NULL;
1380:   hA->available = PETSC_TRUE;
1381:   return(0);
1382: }

1384: static PetscErrorCode MatSetValues_HYPRE(Mat A, PetscInt nr, const PetscInt rows[], PetscInt nc, const PetscInt cols[], const PetscScalar v[], InsertMode ins)
1385: {
1386:   Mat_HYPRE      *hA = (Mat_HYPRE*)A->data;
1387:   PetscScalar    *vals = (PetscScalar *)v;
1388:   HYPRE_Complex  *sscr;
1389:   PetscInt       *cscr[2];
1390:   PetscInt       i,nzc;
1391:   void           *array = NULL;

1395:   MatGetArray_HYPRE(A,sizeof(PetscInt)*(2*nc)+sizeof(HYPRE_Complex)*nc*nr,&array);
1396:   cscr[0] = (PetscInt*)array;
1397:   cscr[1] = ((PetscInt*)array)+nc;
1398:   sscr = (HYPRE_Complex*)(((PetscInt*)array)+nc*2);
1399:   for (i=0,nzc=0;i<nc;i++) {
1400:     if (cols[i] >= 0) {
1401:       cscr[0][nzc  ] = cols[i];
1402:       cscr[1][nzc++] = i;
1403:     }
1404:   }
1405:   if (!nzc) {
1406:     MatRestoreArray_HYPRE(A,&array);
1407:     return(0);
1408:   }

1410: #if 0 //defined(PETSC_HAVE_HYPRE_DEVICE)
1411:   if (HYPRE_MEMORY_HOST != hypre_IJMatrixMemoryLocation(hA->ij)) {
1412:     hypre_ParCSRMatrix *parcsr;

1414:     PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hA->ij,(void**)&parcsr));
1415:     PetscStackCallStandard(hypre_ParCSRMatrixMigrate,(parcsr, HYPRE_MEMORY_HOST));
1416:   }
1417: #endif

1419:   if (ins == ADD_VALUES) {
1420:     for (i=0;i<nr;i++) {
1421:       if (rows[i] >= 0) {
1422:         PetscInt  j;
1423:         HYPRE_Int hnc = (HYPRE_Int)nzc;

1425:         if ((PetscInt)hnc != nzc) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Hypre overflow! number of columns %D for row %D",nzc,rows[i]);
1426:         for (j=0;j<nzc;j++) { PetscHYPREScalarCast(vals[cscr[1][j]],&sscr[j]); }
1427:         PetscStackCallStandard(HYPRE_IJMatrixAddToValues,(hA->ij,1,&hnc,(HYPRE_BigInt*)(rows+i),(HYPRE_BigInt*)cscr[0],sscr));
1428:       }
1429:       vals += nc;
1430:     }
1431:   } else { /* INSERT_VALUES */
1432:     PetscInt rst,ren;

1434:     MatGetOwnershipRange(A,&rst,&ren);
1435:     for (i=0;i<nr;i++) {
1436:       if (rows[i] >= 0) {
1437:         PetscInt  j;
1438:         HYPRE_Int hnc = (HYPRE_Int)nzc;

1440:         if ((PetscInt)hnc != nzc) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Hypre overflow! number of columns %D for row %D",nzc,rows[i]);
1441:         for (j=0;j<nzc;j++) { PetscHYPREScalarCast(vals[cscr[1][j]],&sscr[j]); }
1442:         /* nonlocal values */
1443:         if (rows[i] < rst || rows[i] >= ren) { MatStashValuesRow_Private(&A->stash,rows[i],nzc,cscr[0],(PetscScalar*)sscr,PETSC_FALSE); }
1444:         /* local values */
1445:         else PetscStackCallStandard(HYPRE_IJMatrixSetValues,(hA->ij,1,&hnc,(HYPRE_BigInt*)(rows+i),(HYPRE_BigInt*)cscr[0],sscr));
1446:       }
1447:       vals += nc;
1448:     }
1449:   }

1451:   MatRestoreArray_HYPRE(A,&array);
1452:   return(0);
1453: }

1455: static PetscErrorCode MatHYPRESetPreallocation_HYPRE(Mat A, PetscInt dnz, const PetscInt dnnz[], PetscInt onz, const PetscInt onnz[])
1456: {
1457:   Mat_HYPRE      *hA = (Mat_HYPRE*)A->data;
1458:   HYPRE_Int      *hdnnz,*honnz;
1459:   PetscInt       i,rs,re,cs,ce,bs;
1460:   PetscMPIInt    size;

1464:   MatGetBlockSize(A,&bs);
1465:   PetscLayoutSetUp(A->rmap);
1466:   PetscLayoutSetUp(A->cmap);
1467:   rs   = A->rmap->rstart;
1468:   re   = A->rmap->rend;
1469:   cs   = A->cmap->rstart;
1470:   ce   = A->cmap->rend;
1471:   if (!hA->ij) {
1472:     PetscStackCallStandard(HYPRE_IJMatrixCreate,(hA->comm,rs,re-1,cs,ce-1,&hA->ij));
1473:     PetscStackCallStandard(HYPRE_IJMatrixSetObjectType,(hA->ij,HYPRE_PARCSR));
1474:   } else {
1475:     HYPRE_BigInt hrs,hre,hcs,hce;
1476:     PetscStackCallStandard(HYPRE_IJMatrixGetLocalRange,(hA->ij,&hrs,&hre,&hcs,&hce));
1477:     if (hre-hrs+1 != re -rs) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Inconsistent local rows: IJMatrix [%D,%D), PETSc [%D,%D)",hrs,hre+1,rs,re);
1478:     if (hce-hcs+1 != ce -cs) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Inconsistent local cols: IJMatrix [%D,%D), PETSc [%D,%D)",hcs,hce+1,cs,ce);
1479:   }
1480:   if (dnz == PETSC_DEFAULT || dnz == PETSC_DECIDE) dnz = 10*bs;
1481:   if (onz == PETSC_DEFAULT || onz == PETSC_DECIDE) onz = 10*bs;

1483:   if (!dnnz) {
1484:     PetscMalloc1(A->rmap->n,&hdnnz);
1485:     for (i=0;i<A->rmap->n;i++) hdnnz[i] = dnz;
1486:   } else {
1487:     hdnnz = (HYPRE_Int*)dnnz;
1488:   }
1489:   MPI_Comm_size(PetscObjectComm((PetscObject)A),&size);
1490:   if (size > 1) {
1491:     hypre_AuxParCSRMatrix *aux_matrix;
1492:     if (!onnz) {
1493:       PetscMalloc1(A->rmap->n,&honnz);
1494:       for (i=0;i<A->rmap->n;i++) honnz[i] = onz;
1495:     } else honnz = (HYPRE_Int*)onnz;
1496:     /* SetDiagOffdSizes sets hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 0, since it seems
1497:        they assume the user will input the entire row values, properly sorted
1498:        In PETSc, we don't make such an assumption and set this flag to 1,
1499:        unless the option MAT_SORTED_FULL is set to true.
1500:        Also, to avoid possible memory leaks, we destroy and recreate the translator
1501:        This has to be done here, as HYPRE_IJMatrixInitialize will properly initialize
1502:        the IJ matrix for us */
1503:     aux_matrix = (hypre_AuxParCSRMatrix*)hypre_IJMatrixTranslator(hA->ij);
1504:     hypre_AuxParCSRMatrixDestroy(aux_matrix);
1505:     hypre_IJMatrixTranslator(hA->ij) = NULL;
1506:     PetscStackCallStandard(HYPRE_IJMatrixSetDiagOffdSizes,(hA->ij,hdnnz,honnz));
1507:     aux_matrix = (hypre_AuxParCSRMatrix*)hypre_IJMatrixTranslator(hA->ij);
1508:     hypre_AuxParCSRMatrixNeedAux(aux_matrix) = !hA->sorted_full;
1509:   } else {
1510:     honnz = NULL;
1511:     PetscStackCallStandard(HYPRE_IJMatrixSetRowSizes,(hA->ij,hdnnz));
1512:   }

1514:   /* reset assembled flag and call the initialize method */
1515:   hypre_IJMatrixAssembleFlag(hA->ij) = 0;
1516: #if PETSC_PKG_HYPRE_VERSION_LT(2,19,0)
1517:   PetscStackCallStandard(HYPRE_IJMatrixInitialize,(hA->ij));
1518: #else
1519:   PetscStackCallStandard(HYPRE_IJMatrixInitialize_v2,(hA->ij,HYPRE_MEMORY_HOST));
1520: #endif
1521:   if (!dnnz) {
1522:     PetscFree(hdnnz);
1523:   }
1524:   if (!onnz && honnz) {
1525:     PetscFree(honnz);
1526:   }
1527:   /* Match AIJ logic */
1528:   A->preallocated = PETSC_TRUE;
1529:   A->assembled    = PETSC_FALSE;
1530:   return(0);
1531: }

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

1536:    Collective on Mat

1538:    Input Parameters:
1539: +  A - the matrix
1540: .  dnz  - number of nonzeros per row in DIAGONAL portion of local submatrix
1541:           (same value is used for all local rows)
1542: .  dnnz - array containing the number of nonzeros in the various rows of the
1543:           DIAGONAL portion of the local submatrix (possibly different for each row)
1544:           or NULL (PETSC_NULL_INTEGER in Fortran), if d_nz is used to specify the nonzero structure.
1545:           The size of this array is equal to the number of local rows, i.e 'm'.
1546:           For matrices that will be factored, you must leave room for (and set)
1547:           the diagonal entry even if it is zero.
1548: .  onz  - number of nonzeros per row in the OFF-DIAGONAL portion of local
1549:           submatrix (same value is used for all local rows).
1550: -  onnz - array containing the number of nonzeros in the various rows of the
1551:           OFF-DIAGONAL portion of the local submatrix (possibly different for
1552:           each row) or NULL (PETSC_NULL_INTEGER in Fortran), if o_nz is used to specify the nonzero
1553:           structure. The size of this array is equal to the number
1554:           of local rows, i.e 'm'.

1556:    Notes:
1557:     If the *nnz parameter is given then the *nz parameter is ignored; for sequential matrices, onz and onnz are ignored.

1559:    Level: intermediate

1561: .seealso: MatCreate(), MatMPIAIJSetPreallocation(), MATHYPRE
1562: @*/
1563: PetscErrorCode MatHYPRESetPreallocation(Mat A, PetscInt dnz, const PetscInt dnnz[], PetscInt onz, const PetscInt onnz[])
1564: {

1570:   PetscTryMethod(A,"MatHYPRESetPreallocation_C",(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt[]),(A,dnz,dnnz,onz,onnz));
1571:   return(0);
1572: }

1574: /*
1575:    MatCreateFromParCSR - Creates a matrix from a hypre_ParCSRMatrix

1577:    Collective

1579:    Input Parameters:
1580: +  parcsr   - the pointer to the hypre_ParCSRMatrix
1581: .  mtype    - matrix type to be created. Currently MATAIJ, MATIS and MATHYPRE are supported.
1582: -  copymode - PETSc copying options

1584:    Output Parameter:
1585: .  A  - the matrix

1587:    Level: intermediate

1589: .seealso: MatHYPRE, PetscCopyMode
1590: */
1591: PETSC_EXTERN PetscErrorCode MatCreateFromParCSR(hypre_ParCSRMatrix *parcsr, MatType mtype, PetscCopyMode copymode, Mat* A)
1592: {
1593:   Mat            T;
1594:   Mat_HYPRE      *hA;
1595:   MPI_Comm       comm;
1596:   PetscInt       rstart,rend,cstart,cend,M,N;
1597:   PetscBool      isseqaij,isseqaijmkl,ismpiaij,isaij,ishyp,isis;

1601:   comm  = hypre_ParCSRMatrixComm(parcsr);
1602:   PetscStrcmp(mtype,MATSEQAIJ,&isseqaij);
1603:   PetscStrcmp(mtype,MATSEQAIJMKL,&isseqaijmkl);
1604:   PetscStrcmp(mtype,MATMPIAIJ,&ismpiaij);
1605:   PetscStrcmp(mtype,MATAIJ,&isaij);
1606:   PetscStrcmp(mtype,MATHYPRE,&ishyp);
1607:   PetscStrcmp(mtype,MATIS,&isis);
1608:   isaij = (PetscBool)(isseqaij || isseqaijmkl || ismpiaij || isaij);
1609:   /* TODO */
1610:   if (!isaij && !ishyp && !isis) SETERRQ7(comm,PETSC_ERR_SUP,"Unsupported MatType %s! Supported types are %s, %s, %s, %s, %s, and %s",mtype,MATAIJ,MATSEQAIJ,MATSEQAIJMKL,MATMPIAIJ,MATIS,MATHYPRE);
1611:   /* access ParCSRMatrix */
1612:   rstart = hypre_ParCSRMatrixFirstRowIndex(parcsr);
1613:   rend   = hypre_ParCSRMatrixLastRowIndex(parcsr);
1614:   cstart = hypre_ParCSRMatrixFirstColDiag(parcsr);
1615:   cend   = hypre_ParCSRMatrixLastColDiag(parcsr);
1616:   M      = hypre_ParCSRMatrixGlobalNumRows(parcsr);
1617:   N      = hypre_ParCSRMatrixGlobalNumCols(parcsr);

1619:   /* fix for empty local rows/columns */
1620:   if (rend < rstart) rend = rstart;
1621:   if (cend < cstart) cend = cstart;

1623:   /* PETSc convention */
1624:   rend++;
1625:   cend++;
1626:   rend = PetscMin(rend,M);
1627:   cend = PetscMin(cend,N);

1629:   /* create PETSc matrix with MatHYPRE */
1630:   MatCreate(comm,&T);
1631:   MatSetSizes(T,rend-rstart,cend-cstart,M,N);
1632:   MatSetType(T,MATHYPRE);
1633:   hA   = (Mat_HYPRE*)(T->data);

1635:   /* create HYPRE_IJMatrix */
1636:   PetscStackCallStandard(HYPRE_IJMatrixCreate,(hA->comm,rstart,rend-1,cstart,cend-1,&hA->ij));
1637:   PetscStackCallStandard(HYPRE_IJMatrixSetObjectType,(hA->ij,HYPRE_PARCSR));

1639: // TODO DEV
1640:   /* create new ParCSR object if needed */
1641:   if (ishyp && copymode == PETSC_COPY_VALUES) {
1642:     hypre_ParCSRMatrix *new_parcsr;
1643: #if PETSC_PKG_HYPRE_VERSION_LT(2,18,0)
1644:     hypre_CSRMatrix    *hdiag,*hoffd,*ndiag,*noffd;

1646:     new_parcsr = hypre_ParCSRMatrixClone(parcsr,0);
1647:     hdiag      = hypre_ParCSRMatrixDiag(parcsr);
1648:     hoffd      = hypre_ParCSRMatrixOffd(parcsr);
1649:     ndiag      = hypre_ParCSRMatrixDiag(new_parcsr);
1650:     noffd      = hypre_ParCSRMatrixOffd(new_parcsr);
1651:     PetscArraycpy(hypre_CSRMatrixData(ndiag),hypre_CSRMatrixData(hdiag),hypre_CSRMatrixNumNonzeros(hdiag));
1652:     PetscArraycpy(hypre_CSRMatrixData(noffd),hypre_CSRMatrixData(hoffd),hypre_CSRMatrixNumNonzeros(hoffd));
1653: #else
1654:     new_parcsr = hypre_ParCSRMatrixClone(parcsr,1);
1655: #endif
1656:     parcsr     = new_parcsr;
1657:     copymode   = PETSC_OWN_POINTER;
1658:   }

1660:   /* set ParCSR object */
1661:   hypre_IJMatrixObject(hA->ij) = parcsr;
1662:   T->preallocated = PETSC_TRUE;

1664:   /* set assembled flag */
1665:   hypre_IJMatrixAssembleFlag(hA->ij) = 1;
1666: #if 0
1667:   PetscStackCallStandard(HYPRE_IJMatrixInitialize,(hA->ij));
1668: #endif
1669:   if (ishyp) {
1670:     PetscMPIInt myid = 0;

1672:     /* make sure we always have row_starts and col_starts available */
1673:     if (HYPRE_AssumedPartitionCheck()) {
1674:       MPI_Comm_rank(comm,&myid);
1675:     }
1676: #if defined(hypre_ParCSRMatrixOwnsRowStarts)
1677:     if (!hypre_ParCSRMatrixOwnsColStarts(parcsr)) {
1678:       PetscLayout map;

1680:       MatGetLayouts(T,NULL,&map);
1681:       PetscLayoutSetUp(map);
1682:       hypre_ParCSRMatrixColStarts(parcsr) = (HYPRE_BigInt*)(map->range + myid);
1683:     }
1684:     if (!hypre_ParCSRMatrixOwnsRowStarts(parcsr)) {
1685:       PetscLayout map;

1687:       MatGetLayouts(T,&map,NULL);
1688:       PetscLayoutSetUp(map);
1689:       hypre_ParCSRMatrixRowStarts(parcsr) = (HYPRE_BigInt*)(map->range + myid);
1690:     }
1691: #endif
1692:     /* prevent from freeing the pointer */
1693:     if (copymode == PETSC_USE_POINTER) hA->inner_free = PETSC_FALSE;
1694:     *A   = T;
1695:     MatSetOption(*A,MAT_SORTED_FULL,PETSC_TRUE);
1696:     MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);
1697:     MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);
1698:   } else if (isaij) {
1699:     if (copymode != PETSC_OWN_POINTER) {
1700:       /* prevent from freeing the pointer */
1701:       hA->inner_free = PETSC_FALSE;
1702:       MatConvert_HYPRE_AIJ(T,MATAIJ,MAT_INITIAL_MATRIX,A);
1703:       MatDestroy(&T);
1704:     } else { /* AIJ return type with PETSC_OWN_POINTER */
1705:       MatConvert_HYPRE_AIJ(T,MATAIJ,MAT_INPLACE_MATRIX,&T);
1706:       *A   = T;
1707:     }
1708:   } else if (isis) {
1709:     MatConvert_HYPRE_IS(T,MATIS,MAT_INITIAL_MATRIX,A);
1710:     if (copymode != PETSC_OWN_POINTER) hA->inner_free = PETSC_FALSE;
1711:     MatDestroy(&T);
1712:   }
1713:   return(0);
1714: }

1716: static PetscErrorCode MatHYPREGetParCSR_HYPRE(Mat A, hypre_ParCSRMatrix **parcsr)
1717: {
1718:   Mat_HYPRE *hA = (Mat_HYPRE*)A->data;
1719:   HYPRE_Int type;

1722:   if (!hA->ij) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"HYPRE_IJMatrix not present");
1723:   PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(hA->ij,&type));
1724:   if (type != HYPRE_PARCSR) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"HYPRE_IJMatrix is not of type HYPRE_PARCSR");
1725:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hA->ij,(void**)parcsr));
1726:   return(0);
1727: }

1729: /*
1730:    MatHYPREGetParCSR - Gets the pointer to the ParCSR matrix

1732:    Not collective

1734:    Input Parameters:
1735: +  A  - the MATHYPRE object

1737:    Output Parameter:
1738: .  parcsr  - the pointer to the hypre_ParCSRMatrix

1740:    Level: intermediate

1742: .seealso: MatHYPRE, PetscCopyMode
1743: */
1744: PetscErrorCode MatHYPREGetParCSR(Mat A, hypre_ParCSRMatrix **parcsr)
1745: {

1751:   PetscUseMethod(A,"MatHYPREGetParCSR_C",(Mat,hypre_ParCSRMatrix**),(A,parcsr));
1752:   return(0);
1753: }

1755: static PetscErrorCode MatMissingDiagonal_HYPRE(Mat A, PetscBool *missing, PetscInt *dd)
1756: {
1757:   hypre_ParCSRMatrix *parcsr;
1758:   hypre_CSRMatrix    *ha;
1759:   PetscInt           rst;
1760:   PetscErrorCode     ierr;

1763:   if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not implemented with non-square diagonal blocks");
1764:   MatGetOwnershipRange(A,&rst,NULL);
1765:   MatHYPREGetParCSR_HYPRE(A,&parcsr);
1766:   if (missing) *missing = PETSC_FALSE;
1767:   if (dd) *dd = -1;
1768:   ha = hypre_ParCSRMatrixDiag(parcsr);
1769:   if (ha) {
1770:     PetscInt  size,i;
1771:     HYPRE_Int *ii,*jj;

1773:     size = hypre_CSRMatrixNumRows(ha);
1774:     ii   = hypre_CSRMatrixI(ha);
1775:     jj   = hypre_CSRMatrixJ(ha);
1776:     for (i = 0; i < size; i++) {
1777:       PetscInt  j;
1778:       PetscBool found = PETSC_FALSE;

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

1783:       if (!found) {
1784:         PetscInfo1(A,"Matrix is missing local diagonal entry %D\n",i);
1785:         if (missing) *missing = PETSC_TRUE;
1786:         if (dd) *dd = i+rst;
1787:         return(0);
1788:       }
1789:     }
1790:     if (!size) {
1791:       PetscInfo(A,"Matrix has no diagonal entries therefore is missing diagonal\n");
1792:       if (missing) *missing = PETSC_TRUE;
1793:       if (dd) *dd = rst;
1794:     }
1795:   } else {
1796:     PetscInfo(A,"Matrix has no diagonal entries therefore is missing diagonal\n");
1797:     if (missing) *missing = PETSC_TRUE;
1798:     if (dd) *dd = rst;
1799:   }
1800:   return(0);
1801: }

1803: static PetscErrorCode MatScale_HYPRE(Mat A, PetscScalar s)
1804: {
1805:   hypre_ParCSRMatrix *parcsr;
1806: #if PETSC_PKG_HYPRE_VERSION_LT(2,19,0)
1807:   hypre_CSRMatrix    *ha;
1808: #endif
1809:   PetscErrorCode     ierr;
1810:   HYPRE_Complex      hs;

1813:   PetscHYPREScalarCast(s,&hs);
1814:   MatHYPREGetParCSR_HYPRE(A,&parcsr);
1815: #if PETSC_PKG_HYPRE_VERSION_GE(2,19,0)
1816:   PetscStackCallStandard(hypre_ParCSRMatrixScale,(parcsr,hs));
1817: #else  /* diagonal part */
1818:   ha = hypre_ParCSRMatrixDiag(parcsr);
1819:   if (ha) {
1820:     PetscInt      size,i;
1821:     HYPRE_Int     *ii;
1822:     HYPRE_Complex *a;

1824:     size = hypre_CSRMatrixNumRows(ha);
1825:     a    = hypre_CSRMatrixData(ha);
1826:     ii   = hypre_CSRMatrixI(ha);
1827:     for (i = 0; i < ii[size]; i++) a[i] *= hs;
1828:   }
1829:   /* offdiagonal part */
1830:   ha = hypre_ParCSRMatrixOffd(parcsr);
1831:   if (ha) {
1832:     PetscInt      size,i;
1833:     HYPRE_Int     *ii;
1834:     HYPRE_Complex *a;

1836:     size = hypre_CSRMatrixNumRows(ha);
1837:     a    = hypre_CSRMatrixData(ha);
1838:     ii   = hypre_CSRMatrixI(ha);
1839:     for (i = 0; i < ii[size]; i++) a[i] *= hs;
1840:   }
1841: #endif
1842:   return(0);
1843: }

1845: static PetscErrorCode MatZeroRowsColumns_HYPRE(Mat A, PetscInt numRows, const PetscInt rows[], PetscScalar diag, Vec x, Vec b)
1846: {
1847:   hypre_ParCSRMatrix *parcsr;
1848:   HYPRE_Int          *lrows;
1849:   PetscInt           rst,ren,i;
1850:   PetscErrorCode     ierr;

1853:   if (x || b) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"To be implemented");
1854:   MatHYPREGetParCSR_HYPRE(A,&parcsr);
1855:   PetscMalloc1(numRows,&lrows);
1856:   MatGetOwnershipRange(A,&rst,&ren);
1857:   for (i=0;i<numRows;i++) {
1858:     if (rows[i] < rst || rows[i] >= ren)
1859:       SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Non-local rows not yet supported");
1860:     lrows[i] = rows[i] - rst;
1861:   }
1862:   PetscStackCallStandard(hypre_ParCSRMatrixEliminateRowsCols,(parcsr,numRows,lrows));
1863:   PetscFree(lrows);
1864:   return(0);
1865: }

1867: static PetscErrorCode MatZeroEntries_HYPRE_CSRMatrix(hypre_CSRMatrix *ha)
1868: {
1869:   PetscErrorCode      ierr;

1872:   if (ha) {
1873:     HYPRE_Int     *ii, size;
1874:     HYPRE_Complex *a;

1876:     size = hypre_CSRMatrixNumRows(ha);
1877:     a    = hypre_CSRMatrixData(ha);
1878:     ii   = hypre_CSRMatrixI(ha);

1880:     if (a) {PetscArrayzero(a,ii[size]);}
1881:   }
1882:   return(0);
1883: }

1885: PetscErrorCode MatZeroEntries_HYPRE(Mat A)
1886: {
1887:   Mat_HYPRE *hA = (Mat_HYPRE*)A->data;

1890:   if (HYPRE_MEMORY_DEVICE == hypre_IJMatrixMemoryLocation(hA->ij)) {
1891:     PetscStackCallStandard(HYPRE_IJMatrixSetConstantValues,(hA->ij,0.0));
1892:   } else {
1893:     hypre_ParCSRMatrix *parcsr;
1894:     PetscErrorCode     ierr;

1896:     MatHYPREGetParCSR_HYPRE(A,&parcsr);
1897:     MatZeroEntries_HYPRE_CSRMatrix(hypre_ParCSRMatrixDiag(parcsr));
1898:     MatZeroEntries_HYPRE_CSRMatrix(hypre_ParCSRMatrixOffd(parcsr));
1899:   }
1900:   return(0);
1901: }

1903: static PetscErrorCode MatZeroRows_HYPRE_CSRMatrix(hypre_CSRMatrix *hA,PetscInt N,const PetscInt rows[],HYPRE_Complex diag)
1904: {
1905:   PetscInt        ii;
1906:   HYPRE_Int       *i, *j;
1907:   HYPRE_Complex   *a;

1910:   if (!hA) return(0);

1912:   i = hypre_CSRMatrixI(hA);
1913:   j = hypre_CSRMatrixJ(hA);
1914:   a = hypre_CSRMatrixData(hA);

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

1919:     irow = rows[ii];
1920:     ibeg = i[irow];
1921:     iend = i[irow+1];
1922:     for (jj = ibeg; jj < iend; jj++)
1923:       if (j[jj] == irow) a[jj] = diag;
1924:       else a[jj] = 0.0;
1925:    }
1926:    return(0);
1927: }

1929: static PetscErrorCode MatZeroRows_HYPRE(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
1930: {
1931:   hypre_ParCSRMatrix  *parcsr;
1932:   PetscInt            *lrows,len;
1933:   HYPRE_Complex       hdiag;
1934:   PetscErrorCode      ierr;

1937:   if (x || b) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Does not support to modify the solution and the right hand size");
1938:   PetscHYPREScalarCast(diag,&hdiag);
1939:   /* retrieve the internal matrix */
1940:   MatHYPREGetParCSR_HYPRE(A,&parcsr);
1941:   /* get locally owned rows */
1942:   MatZeroRowsMapLocal_Private(A,N,rows,&len,&lrows);
1943:   /* zero diagonal part */
1944:   MatZeroRows_HYPRE_CSRMatrix(hypre_ParCSRMatrixDiag(parcsr),len,lrows,hdiag);
1945:   /* zero off-diagonal part */
1946:   MatZeroRows_HYPRE_CSRMatrix(hypre_ParCSRMatrixOffd(parcsr),len,lrows,0.0);

1948:   PetscFree(lrows);
1949:   return(0);
1950: }

1952: static PetscErrorCode MatAssemblyBegin_HYPRE(Mat mat,MatAssemblyType mode)
1953: {

1957:   if (mat->nooffprocentries) return(0);

1959:   MatStashScatterBegin_Private(mat,&mat->stash,mat->rmap->range);
1960:   return(0);
1961: }

1963: static PetscErrorCode MatGetRow_HYPRE(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v)
1964: {
1965:   hypre_ParCSRMatrix  *parcsr;
1966:   HYPRE_Int           hnz;
1967:   PetscErrorCode      ierr;

1970:   /* retrieve the internal matrix */
1971:   MatHYPREGetParCSR_HYPRE(A,&parcsr);
1972:   /* call HYPRE API */
1973:   PetscStackCallStandard(HYPRE_ParCSRMatrixGetRow,(parcsr,row,&hnz,(HYPRE_BigInt**)idx,(HYPRE_Complex**)v));
1974:   if (nz) *nz = (PetscInt)hnz;
1975:   return(0);
1976: }

1978: static PetscErrorCode MatRestoreRow_HYPRE(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v)
1979: {
1980:   hypre_ParCSRMatrix  *parcsr;
1981:   HYPRE_Int           hnz;
1982:   PetscErrorCode      ierr;

1985:   /* retrieve the internal matrix */
1986:   MatHYPREGetParCSR_HYPRE(A,&parcsr);
1987:   /* call HYPRE API */
1988:   hnz = nz ? (HYPRE_Int)(*nz) : 0;
1989:   PetscStackCallStandard(HYPRE_ParCSRMatrixRestoreRow,(parcsr,row,&hnz,(HYPRE_BigInt**)idx,(HYPRE_Complex**)v));
1990:   return(0);
1991: }

1993: static PetscErrorCode MatGetValues_HYPRE(Mat A,PetscInt m,const PetscInt idxm[],PetscInt n,const PetscInt idxn[],PetscScalar v[])
1994: {
1995:   Mat_HYPRE *hA = (Mat_HYPRE*)A->data;
1996:   PetscInt  i;

1999:   if (!m || !n) return(0);
2000:   /* Ignore negative row indices
2001:    * And negative column indices should be automatically ignored in hypre
2002:    * */
2003:   for (i=0; i<m; i++) {
2004:     if (idxm[i] >= 0) {
2005:       HYPRE_Int hn = (HYPRE_Int)n;
2006:       PetscStackCallStandard(HYPRE_IJMatrixGetValues,(hA->ij,1,&hn,(HYPRE_BigInt*)&idxm[i],(HYPRE_BigInt*)idxn,(HYPRE_Complex*)(v + i*n)));
2007:     }
2008:   }
2009:   return(0);
2010: }

2012: static PetscErrorCode MatSetOption_HYPRE(Mat A,MatOption op,PetscBool flg)
2013: {
2014:   Mat_HYPRE *hA = (Mat_HYPRE*)A->data;

2017:   switch (op) {
2018:   case MAT_NO_OFF_PROC_ENTRIES:
2019:     if (flg) {
2020:       PetscStackCallStandard(HYPRE_IJMatrixSetMaxOffProcElmts,(hA->ij,0));
2021:     }
2022:     break;
2023:   case MAT_SORTED_FULL:
2024:     hA->sorted_full = flg;
2025:     break;
2026:   default:
2027:     break;
2028:   }
2029:   return(0);
2030: }

2032: static PetscErrorCode MatView_HYPRE(Mat A, PetscViewer view)
2033: {
2034:   PetscErrorCode     ierr;
2035:   PetscViewerFormat  format;

2038:   PetscViewerGetFormat(view,&format);
2039:   if (format == PETSC_VIEWER_ASCII_FACTOR_INFO || format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) return(0);
2040:   if (format != PETSC_VIEWER_NATIVE) {
2041:     Mat                B;
2042:     hypre_ParCSRMatrix *parcsr;
2043:     PetscErrorCode     (*mview)(Mat,PetscViewer) = NULL;

2045:     MatHYPREGetParCSR_HYPRE(A,&parcsr);
2046:     MatCreateFromParCSR(parcsr,MATAIJ,PETSC_USE_POINTER,&B);
2047:     MatGetOperation(B,MATOP_VIEW,(void(**)(void))&mview);
2048:     if (!mview) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Missing view operation");
2049:     (*mview)(B,view);
2050:     MatDestroy(&B);
2051:   } else {
2052:     Mat_HYPRE  *hA = (Mat_HYPRE*)A->data;
2053:     PetscMPIInt size;
2054:     PetscBool   isascii;
2055:     const char *filename;

2057:     /* HYPRE uses only text files */
2058:     PetscObjectTypeCompare((PetscObject)view,PETSCVIEWERASCII,&isascii);
2059:     if (!isascii) SETERRQ1(PetscObjectComm((PetscObject)view),PETSC_ERR_SUP,"PetscViewerType %s: native HYPRE format needs PETSCVIEWERASCII",((PetscObject)view)->type_name);
2060:     PetscViewerFileGetName(view,&filename);
2061:     PetscStackCallStandard(HYPRE_IJMatrixPrint,(hA->ij,filename));
2062:     MPI_Comm_size(hA->comm,&size);
2063:     if (size > 1) {
2064:       PetscViewerASCIIPrintf(view,"Matrix files: %s.%05d ... %s.%05d\n",filename,0,filename,size-1);
2065:     } else {
2066:       PetscViewerASCIIPrintf(view,"Matrix file: %s.%05d\n",filename,0);
2067:     }
2068:   }
2069:   return(0);
2070: }

2072: static PetscErrorCode MatDuplicate_HYPRE(Mat A,MatDuplicateOption op, Mat *B)
2073: {
2074:   hypre_ParCSRMatrix *parcsr = NULL;
2075:   PetscErrorCode     ierr;
2076:   PetscCopyMode      cpmode;

2079:   MatHYPREGetParCSR_HYPRE(A,&parcsr);
2080:   if (op == MAT_DO_NOT_COPY_VALUES || op == MAT_SHARE_NONZERO_PATTERN) {
2081:     parcsr = hypre_ParCSRMatrixClone(parcsr,0);
2082:     cpmode = PETSC_OWN_POINTER;
2083:   } else {
2084:     cpmode = PETSC_COPY_VALUES;
2085:   }
2086:   MatCreateFromParCSR(parcsr,MATHYPRE,cpmode,B);
2087:   return(0);
2088: }

2090: static PetscErrorCode MatCopy_HYPRE(Mat A, Mat B, MatStructure str)
2091: {
2092:   hypre_ParCSRMatrix *acsr,*bcsr;
2093:   PetscErrorCode     ierr;

2096:   if (str == SAME_NONZERO_PATTERN && A->ops->copy == B->ops->copy) {
2097:     MatHYPREGetParCSR_HYPRE(A,&acsr);
2098:     MatHYPREGetParCSR_HYPRE(B,&bcsr);
2099:     PetscStackCallStandard(hypre_ParCSRMatrixCopy,(acsr,bcsr,1));
2100:     MatSetOption(B,MAT_SORTED_FULL,PETSC_TRUE); /* "perfect" preallocation, so no need for hypre_AuxParCSRMatrixNeedAux */
2101:     MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
2102:     MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
2103:   } else {
2104:     MatCopy_Basic(A,B,str);
2105:   }
2106:   return(0);
2107: }

2109: static PetscErrorCode MatGetDiagonal_HYPRE(Mat A, Vec d)
2110: {
2111:   hypre_ParCSRMatrix *parcsr;
2112:   hypre_CSRMatrix    *dmat;
2113:   HYPRE_Complex      *a;
2114:   HYPRE_Complex      *data = NULL;
2115:   HYPRE_Int          *diag = NULL;
2116:   PetscInt           i;
2117:   PetscBool          cong;
2118:   PetscErrorCode     ierr;

2121:   MatHasCongruentLayouts(A,&cong);
2122:   if (!cong) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Only for square matrices with same local distributions of rows and columns");
2123:   if (PetscDefined(USE_DEBUG)) {
2124:     PetscBool miss;
2125:     MatMissingDiagonal(A,&miss,NULL);
2126:     if (miss && A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not implemented when diagonal entries are missing");
2127:   }
2128:   MatHYPREGetParCSR_HYPRE(A,&parcsr);
2129:   dmat = hypre_ParCSRMatrixDiag(parcsr);
2130:   if (dmat) {
2131:     /* this cast fixes the clang error: implicit conversion from 'HYPRE_Complex' (aka '_Complex double') to 'double' is not permitted in C++ */
2132:     VecGetArray(d,(PetscScalar**)&a);
2133:     diag = hypre_CSRMatrixI(dmat);
2134:     data = hypre_CSRMatrixData(dmat);
2135:     for (i=0;i<A->rmap->n;i++) a[i] = data[diag[i]];
2136:     VecRestoreArray(d,(PetscScalar**)&a);
2137:   }
2138:   return(0);
2139: }

2141: #include <petscblaslapack.h>

2143: static PetscErrorCode MatAXPY_HYPRE(Mat Y,PetscScalar a,Mat X,MatStructure str)
2144: {

2148: #if defined(PETSC_HAVE_HYPRE_DEVICE)
2149:   {
2150:     Mat                B;
2151:     hypre_ParCSRMatrix *x,*y,*z;

2153:     MatHYPREGetParCSR(Y,&y);
2154:     MatHYPREGetParCSR(X,&x);
2155:     PetscStackCallStandard(hypre_ParCSRMatrixAdd,(1.0,y,1.0,x,&z));
2156:     MatCreateFromParCSR(z,MATHYPRE,PETSC_OWN_POINTER,&B);
2157:     MatHeaderMerge(Y,&B);
2158:   }
2159: #else
2160:   if (str == SAME_NONZERO_PATTERN) {
2161:     hypre_ParCSRMatrix *x,*y;
2162:     hypre_CSRMatrix    *xloc,*yloc;
2163:     PetscInt           xnnz,ynnz;
2164:     HYPRE_Complex      *xarr,*yarr;
2165:     PetscBLASInt       one=1,bnz;

2167:     MatHYPREGetParCSR(Y,&y);
2168:     MatHYPREGetParCSR(X,&x);

2170:     /* diagonal block */
2171:     xloc = hypre_ParCSRMatrixDiag(x);
2172:     yloc = hypre_ParCSRMatrixDiag(y);
2173:     xnnz = 0;
2174:     ynnz = 0;
2175:     xarr = NULL;
2176:     yarr = NULL;
2177:     if (xloc) {
2178:       xarr = hypre_CSRMatrixData(xloc);
2179:       xnnz = hypre_CSRMatrixNumNonzeros(xloc);
2180:     }
2181:     if (yloc) {
2182:       yarr = hypre_CSRMatrixData(yloc);
2183:       ynnz = hypre_CSRMatrixNumNonzeros(yloc);
2184:     }
2185:     if (xnnz != ynnz) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Different number of nonzeros in diagonal block %D != %D",xnnz,ynnz);
2186:     PetscBLASIntCast(xnnz,&bnz);
2187:     PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&bnz,&a,(PetscScalar*)xarr,&one,(PetscScalar*)yarr,&one));

2189:     /* off-diagonal block */
2190:     xloc = hypre_ParCSRMatrixOffd(x);
2191:     yloc = hypre_ParCSRMatrixOffd(y);
2192:     xnnz = 0;
2193:     ynnz = 0;
2194:     xarr = NULL;
2195:     yarr = NULL;
2196:     if (xloc) {
2197:       xarr = hypre_CSRMatrixData(xloc);
2198:       xnnz = hypre_CSRMatrixNumNonzeros(xloc);
2199:     }
2200:     if (yloc) {
2201:       yarr = hypre_CSRMatrixData(yloc);
2202:       ynnz = hypre_CSRMatrixNumNonzeros(yloc);
2203:     }
2204:     if (xnnz != ynnz) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Different number of nonzeros in off-diagonal block %D != %D",xnnz,ynnz);
2205:     PetscBLASIntCast(xnnz,&bnz);
2206:     PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&bnz,&a,(PetscScalar*)xarr,&one,(PetscScalar*)yarr,&one));
2207:   } else if (str == SUBSET_NONZERO_PATTERN) {
2208:     MatAXPY_Basic(Y,a,X,str);
2209:   } else {
2210:     Mat B;

2212:     MatAXPY_Basic_Preallocate(Y,X,&B);
2213:     MatAXPY_BasicWithPreallocation(B,Y,a,X,str);
2214:     MatHeaderReplace(Y,&B);
2215:   }
2216: #endif
2217:   return(0);
2218: }

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

2224:    Level: intermediate

2226: .seealso: MatCreate()
2227: M*/

2229: PETSC_EXTERN PetscErrorCode MatCreate_HYPRE(Mat B)
2230: {
2231:   Mat_HYPRE      *hB;

2235:   PetscNewLog(B,&hB);

2237:   hB->inner_free  = PETSC_TRUE;
2238:   hB->available   = PETSC_TRUE;
2239:   hB->sorted_full = PETSC_FALSE; /* no assumption whether column indices are sorted or not */
2240:   hB->size        = 0;
2241:   hB->array       = NULL;

2243:   B->data       = (void*)hB;
2244:   B->rmap->bs   = 1;
2245:   B->assembled  = PETSC_FALSE;

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

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

2280:   MPI_Comm_dup(PetscObjectComm((PetscObject)B),&hB->comm);
2281:   PetscObjectChangeTypeName((PetscObject)B,MATHYPRE);
2282:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_hypre_aij_C",MatConvert_HYPRE_AIJ);
2283:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_hypre_is_C",MatConvert_HYPRE_IS);
2284:   PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqaij_hypre_C",MatProductSetFromOptions_HYPRE);
2285:   PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_mpiaij_hypre_C",MatProductSetFromOptions_HYPRE);
2286:   PetscObjectComposeFunction((PetscObject)B,"MatHYPRESetPreallocation_C",MatHYPRESetPreallocation_HYPRE);
2287:   PetscObjectComposeFunction((PetscObject)B,"MatHYPREGetParCSR_C",MatHYPREGetParCSR_HYPRE);
2288: #if defined(PETSC_HAVE_HYPRE_DEVICE)
2289: #if defined(HYPRE_USING_HIP)
2290:   PetscHIPInitializeCheck();
2291:   MatSetVecType(B,VECHIP);
2292: #endif
2293: #if defined(HYPRE_USING_CUDA)
2294:   PetscCUDAInitializeCheck();
2295:   MatSetVecType(B,VECCUDA);
2296: #endif
2297: #endif
2298:   return(0);
2299: }

2301: static PetscErrorCode hypre_array_destroy(void *ptr)
2302: {
2304:    hypre_TFree(ptr,HYPRE_MEMORY_HOST);
2305:    return(0);
2306: }