Actual source code: mhypre.c

petsc-master 2020-10-23
Report Typos and Errors

  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: 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:   PetscObjectBaseTypeCompare((PetscObject)A,MATMPIAIJ,&flg);
139:   MatGetSize(A,&nr,&nc);
140:   if (flg && nr == nc) {
141:     MatHYPRE_IJMatrixFastCopy_MPIAIJ(A,ij);
142:     return(0);
143:   }
144:   PetscObjectBaseTypeCompare((PetscObject)A,MATSEQAIJ,&flg);
145:   if (flg) {
146:     MatHYPRE_IJMatrixFastCopy_SeqAIJ(A,ij);
147:     return(0);
148:   }

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

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

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

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

190:     for (i=0;i<A->rmap->n + 1;i++) hdiag->i[i] = (HYPRE_Int)pdiag->i[i];
191:     for (i=0;i<pdiag->nz;i++)      hdiag->j[i] = (HYPRE_Int)pdiag->j[i];
192:   }
193:   PetscArraycpy(hdiag->data,pdiag->a,pdiag->nz);

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

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

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

218:   PetscStackCallStandard(HYPRE_IJMatrixInitialize,(ij));
219:   PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(ij,&type));
220:   if (type != HYPRE_PARCSR) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported");
221:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(ij,(void**)&par_matrix));
222:   hdiag = hypre_ParCSRMatrixDiag(par_matrix);
223:   hoffd = hypre_ParCSRMatrixOffd(par_matrix);

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

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

259:   PetscArraycpy(hoffd->data,poffd->a,poffd->nz);

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

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

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

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

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

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

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

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

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

425:   comm = PetscObjectComm((PetscObject)A);
426:   PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(hA->ij,&type));
427:   if (type != HYPRE_PARCSR) SETERRQ(comm,PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported");
428:   if (reuse == MAT_REUSE_MATRIX) {
429:     PetscBool ismpiaij,isseqaij;
430:     PetscObjectBaseTypeCompare((PetscObject)*B,MATMPIAIJ,&ismpiaij);
431:     PetscObjectBaseTypeCompare((PetscObject)*B,MATSEQAIJ,&isseqaij);
432:     if (!ismpiaij && !isseqaij) SETERRQ(comm,PETSC_ERR_SUP,"Only MATMPIAIJ or MATSEQAIJ are supported");
433:   }
434:   MPI_Comm_size(comm,&size);

436:   PetscStackCallStandard(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);
454:       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);
455:       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);
456:       MatSeqAIJGetArray(b->A,&da);
457:     } else {
458:       MatGetRowIJ(*B,0,PETSC_FALSE,PETSC_FALSE,&nr,(const PetscInt**)&dii,(const PetscInt**)&djj,&done);
459:       if (nr != m) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot reuse mat: invalid number of local rows! %D != %D",nr,m);
460:       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);
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:     for (i=0;i<m+1;i++)  dii[i] = (PetscInt)(hypre_CSRMatrixI(hdiag)[i]);
476:     for (i=0;i<dnnz;i++) djj[i] = (PetscInt)(hypre_CSRMatrixJ(hdiag)[i]);
477:   } else {
478:     PetscArraycpy(dii,hypre_CSRMatrixI(hdiag),m+1);
479:     PetscArraycpy(djj,hypre_CSRMatrixJ(hdiag),dnnz);
480:   }
481:   PetscArraycpy(da,hypre_CSRMatrixData(hdiag),dnnz);
482:   iptr = djj;
483:   aptr = da;
484:   for (i=0; i<m; i++) {
485:     PetscInt nc = dii[i+1]-dii[i];
486:     PetscSortIntWithScalarArray(nc,iptr,aptr);
487:     iptr += nc;
488:     aptr += nc;
489:   }
490:   if (size > 1) {
491:     HYPRE_BigInt *coffd;
492:     HYPRE_Int    *offdj;

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

503:       MatGetRowIJ(b->B,0,PETSC_FALSE,PETSC_FALSE,&nr,(const PetscInt**)&oii,(const PetscInt**)&ojj,&done);
504:       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);
505:       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);
506:       MatSeqAIJGetArray(b->B,&oa);
507:     } else { /* MAT_INPLACE_MATRIX */
508:       if (!sameint) {
509:         PetscMalloc1(m+1,&oii);
510:         PetscMalloc1(onnz,&ojj);
511:       } else {
512:         oii = (PetscInt*)hypre_CSRMatrixI(hoffd);
513:         ojj = (PetscInt*)hypre_CSRMatrixJ(hoffd);
514:       }
515:       oa = (PetscScalar*)hypre_CSRMatrixData(hoffd);
516:     }
517:     if (!sameint) {
518:       for (i=0;i<m+1;i++) oii[i] = (PetscInt)(hypre_CSRMatrixI(hoffd)[i]);
519:     } else {
520:       PetscArraycpy(oii,hypre_CSRMatrixI(hoffd),m+1);
521:     }
522:     offdj = hypre_CSRMatrixJ(hoffd);
523:     coffd = hypre_ParCSRMatrixColMapOffd(parcsr);
524:     for (i=0; i<onnz; i++) ojj[i] = coffd[offdj[i]];
525:     PetscArraycpy(oa,hypre_CSRMatrixData(hoffd),onnz);
526:     iptr = ojj;
527:     aptr = oa;
528:     for (i=0; i<m; i++) {
529:        PetscInt nc = oii[i+1]-oii[i];
530:        PetscSortIntWithScalarArray(nc,iptr,aptr);
531:        iptr += nc;
532:        aptr += nc;
533:     }
534:     if (reuse == MAT_INITIAL_MATRIX) {
535:       Mat_MPIAIJ *b;
536:       Mat_SeqAIJ *d,*o;

538:       MatCreateMPIAIJWithSplitArrays(comm,m,n,PETSC_DECIDE,PETSC_DECIDE,dii,djj,da,oii,ojj,oa,B);
539:       /* hack MPIAIJ */
540:       b          = (Mat_MPIAIJ*)((*B)->data);
541:       d          = (Mat_SeqAIJ*)b->A->data;
542:       o          = (Mat_SeqAIJ*)b->B->data;
543:       d->free_a  = PETSC_TRUE;
544:       d->free_ij = PETSC_TRUE;
545:       o->free_a  = PETSC_TRUE;
546:       o->free_ij = PETSC_TRUE;
547:     } else if (reuse == MAT_INPLACE_MATRIX) {
548:       Mat T;

550:       MatCreateMPIAIJWithSplitArrays(comm,m,n,PETSC_DECIDE,PETSC_DECIDE,dii,djj,da,oii,ojj,oa,&T);
551:       if (sameint) { /* ownership of CSR pointers is transferred to PETSc */
552:         hypre_CSRMatrixI(hdiag) = NULL;
553:         hypre_CSRMatrixJ(hdiag) = NULL;
554:         hypre_CSRMatrixI(hoffd) = NULL;
555:         hypre_CSRMatrixJ(hoffd) = NULL;
556:       } else { /* Hack MPIAIJ -> free ij but not a */
557:         Mat_MPIAIJ *b = (Mat_MPIAIJ*)(T->data);
558:         Mat_SeqAIJ *d = (Mat_SeqAIJ*)(b->A->data);
559:         Mat_SeqAIJ *o = (Mat_SeqAIJ*)(b->B->data);

561:         d->free_ij = PETSC_TRUE;
562:         o->free_ij = PETSC_TRUE;
563:       }
564:       hypre_CSRMatrixData(hdiag) = NULL;
565:       hypre_CSRMatrixData(hoffd) = NULL;
566:       MatHeaderReplace(A,&T);
567:     }
568:   } else {
569:     oii  = NULL;
570:     ojj  = NULL;
571:     oa   = NULL;
572:     if (reuse == MAT_INITIAL_MATRIX) {
573:       Mat_SeqAIJ* b;

575:       MatCreateSeqAIJWithArrays(comm,m,n,dii,djj,da,B);
576:       /* hack SeqAIJ */
577:       b          = (Mat_SeqAIJ*)((*B)->data);
578:       b->free_a  = PETSC_TRUE;
579:       b->free_ij = PETSC_TRUE;
580:     } else if (reuse == MAT_INPLACE_MATRIX) {
581:       Mat T;

583:       MatCreateSeqAIJWithArrays(comm,m,n,dii,djj,da,&T);
584:       if (sameint) { /* ownership of CSR pointers is transferred to PETSc */
585:         hypre_CSRMatrixI(hdiag) = NULL;
586:         hypre_CSRMatrixJ(hdiag) = NULL;
587:       } else { /* free ij but not a */
588:         Mat_SeqAIJ* b = (Mat_SeqAIJ*)(T->data);

590:         b->free_ij = PETSC_TRUE;
591:       }
592:       hypre_CSRMatrixData(hdiag) = NULL;
593:       MatHeaderReplace(A,&T);
594:     }
595:   }

597:   /* we have to use hypre_Tfree to free the HYPRE arrays
598:      that PETSc now onws */
599:   if (reuse == MAT_INPLACE_MATRIX) {
600:     PetscInt nh;
601:     void *ptrs[6] = {da,oa,dii,djj,oii,ojj};
602:     const char *names[6] = {"_hypre_csr_da",
603:                             "_hypre_csr_oa",
604:                             "_hypre_csr_dii",
605:                             "_hypre_csr_djj",
606:                             "_hypre_csr_oii",
607:                             "_hypre_csr_ojj"};
608:     nh = sameint ? 6 : 2;
609:     for (i=0; i<nh; i++) {
610:       PetscContainer c;

612:       PetscContainerCreate(comm,&c);
613:       PetscContainerSetPointer(c,ptrs[i]);
614:       PetscContainerSetUserDestroy(c,hypre_array_destroy);
615:       PetscObjectCompose((PetscObject)(*B),names[i],(PetscObject)c);
616:       PetscContainerDestroy(&c);
617:     }
618:   }
619:   return(0);
620: }

622: static PetscErrorCode MatAIJGetParCSR_Private(Mat A, hypre_ParCSRMatrix **hA)
623: {
624:   hypre_ParCSRMatrix *tA;
625:   hypre_CSRMatrix    *hdiag,*hoffd;
626:   Mat_SeqAIJ         *diag,*offd;
627:   PetscInt           *garray,i,noffd,dnnz,onnz,*row_starts,*col_starts;
628:   MPI_Comm           comm = PetscObjectComm((PetscObject)A);
629:   PetscBool          ismpiaij,isseqaij;
630:   PetscBool          sameint = (PetscBool)(sizeof(PetscInt) == sizeof(HYPRE_Int));
631:   PetscErrorCode     ierr;

634:   PetscObjectBaseTypeCompare((PetscObject)A,MATMPIAIJ,&ismpiaij);
635:   PetscObjectBaseTypeCompare((PetscObject)A,MATSEQAIJ,&isseqaij);
636:   if (!ismpiaij && !isseqaij) SETERRQ1(comm,PETSC_ERR_SUP,"Unsupported type %s",((PetscObject)A)->type_name);
637:   if (ismpiaij) {
638:     Mat_MPIAIJ *a = (Mat_MPIAIJ*)(A->data);

640:     diag   = (Mat_SeqAIJ*)a->A->data;
641:     offd   = (Mat_SeqAIJ*)a->B->data;
642:     garray = a->garray;
643:     noffd  = a->B->cmap->N;
644:     dnnz   = diag->nz;
645:     onnz   = offd->nz;
646:   } else {
647:     diag   = (Mat_SeqAIJ*)A->data;
648:     offd   = NULL;
649:     garray = NULL;
650:     noffd  = 0;
651:     dnnz   = diag->nz;
652:     onnz   = 0;
653:   }

655:   /* create a temporary ParCSR */
656:   if (HYPRE_AssumedPartitionCheck()) {
657:     PetscMPIInt myid;

659:     MPI_Comm_rank(comm,&myid);
660:     row_starts = A->rmap->range + myid;
661:     col_starts = A->cmap->range + myid;
662:   } else {
663:     row_starts = A->rmap->range;
664:     col_starts = A->cmap->range;
665:   }
666:   tA = hypre_ParCSRMatrixCreate(comm,A->rmap->N,A->cmap->N,(HYPRE_BigInt*)row_starts,(HYPRE_BigInt*)col_starts,noffd,dnnz,onnz);
667:   hypre_ParCSRMatrixSetRowStartsOwner(tA,0);
668:   hypre_ParCSRMatrixSetColStartsOwner(tA,0);

670:   /* set diagonal part */
671:   hdiag = hypre_ParCSRMatrixDiag(tA);
672:   if (sameint) { /* reuse CSR pointers */
673:     hypre_CSRMatrixI(hdiag) = (HYPRE_Int*)diag->i;
674:     hypre_CSRMatrixJ(hdiag) = (HYPRE_Int*)diag->j;
675:   } else { /* malloc CSR pointers */
676:     HYPRE_Int *hi,*hj;

678:     PetscMalloc2(A->rmap->n+1,&hi,dnnz,&hj);
679:     for (i = 0; i < A->rmap->n+1; i++) hi[i] = (HYPRE_Int)(diag->i[i]);
680:     for (i = 0; i < dnnz; i++)         hj[i] = (HYPRE_Int)(diag->j[i]);
681:     hypre_CSRMatrixI(hdiag) = hi;
682:     hypre_CSRMatrixJ(hdiag) = hj;
683:   }
684:   hypre_CSRMatrixData(hdiag)        = (HYPRE_Complex*)diag->a;
685:   hypre_CSRMatrixNumNonzeros(hdiag) = diag->nz;
686:   hypre_CSRMatrixSetRownnz(hdiag);
687:   hypre_CSRMatrixSetDataOwner(hdiag,0);

689:   /* set offdiagonal part */
690:   hoffd = hypre_ParCSRMatrixOffd(tA);
691:   if (offd) {
692:     if (sameint) { /* reuse CSR pointers */
693:       hypre_CSRMatrixI(hoffd) = (HYPRE_Int*)offd->i;
694:       hypre_CSRMatrixJ(hoffd) = (HYPRE_Int*)offd->j;
695:     } else { /* malloc CSR pointers */
696:       HYPRE_Int *hi,*hj;

698:       PetscMalloc2(A->rmap->n+1,&hi,onnz,&hj);
699:       for (i = 0; i < A->rmap->n+1; i++) hi[i] = (HYPRE_Int)(offd->i[i]);
700:       for (i = 0; i < onnz; i++)         hj[i] = (HYPRE_Int)(offd->j[i]);
701:       hypre_CSRMatrixI(hoffd) = hi;
702:       hypre_CSRMatrixJ(hoffd) = hj;
703:     }
704:     hypre_CSRMatrixData(hoffd)        = (HYPRE_Complex*)offd->a;
705:     hypre_CSRMatrixNumNonzeros(hoffd) = offd->nz;
706:     hypre_CSRMatrixSetRownnz(hoffd);
707:     hypre_CSRMatrixSetDataOwner(hoffd,0);
708:     hypre_ParCSRMatrixSetNumNonzeros(tA);
709:     hypre_ParCSRMatrixColMapOffd(tA) = (HYPRE_BigInt*)garray;
710:   }
711:   *hA = tA;
712:   return(0);
713: }

715: static PetscErrorCode MatAIJRestoreParCSR_Private(Mat A, hypre_ParCSRMatrix **hA)
716: {
717:   hypre_CSRMatrix *hdiag,*hoffd;
718:   PetscBool       sameint = (PetscBool)(sizeof(PetscInt) == sizeof(HYPRE_Int));
719:   PetscErrorCode  ierr;

722:   hdiag = hypre_ParCSRMatrixDiag(*hA);
723:   hoffd = hypre_ParCSRMatrixOffd(*hA);
724:   /* free temporary memory allocated by PETSc */
725:   if (!sameint) {
726:     HYPRE_Int *hi,*hj;

728:     hi = hypre_CSRMatrixI(hdiag);
729:     hj = hypre_CSRMatrixJ(hdiag);
730:     PetscFree2(hi,hj);
731:     if (hoffd) {
732:       hi = hypre_CSRMatrixI(hoffd);
733:       hj = hypre_CSRMatrixJ(hoffd);
734:       PetscFree2(hi,hj);
735:     }
736:   }
737:   /* set pointers to NULL before destroying tA */
738:   hypre_CSRMatrixI(hdiag)           = NULL;
739:   hypre_CSRMatrixJ(hdiag)           = NULL;
740:   hypre_CSRMatrixData(hdiag)        = NULL;
741:   hypre_CSRMatrixI(hoffd)           = NULL;
742:   hypre_CSRMatrixJ(hoffd)           = NULL;
743:   hypre_CSRMatrixData(hoffd)        = NULL;
744:   hypre_ParCSRMatrixColMapOffd(*hA) = NULL;
745:   hypre_ParCSRMatrixDestroy(*hA);
746:   *hA = NULL;
747:   return(0);
748: }

750: /* calls RAP from BoomerAMG:
751:    the resulting ParCSR will not own the column and row starts
752:    It looks like we don't need to have the diagonal entries
753:    ordered first in the rows of the diagonal part
754:    for boomerAMGBuildCoarseOperator to work */
755: static PetscErrorCode MatHYPRE_ParCSR_RAP(hypre_ParCSRMatrix *hR, hypre_ParCSRMatrix *hA,hypre_ParCSRMatrix *hP, hypre_ParCSRMatrix **hRAP)
756: {
757:   HYPRE_Int P_owns_col_starts,R_owns_row_starts;

760:   P_owns_col_starts = hypre_ParCSRMatrixOwnsColStarts(hP);
761:   R_owns_row_starts = hypre_ParCSRMatrixOwnsRowStarts(hR);
762:   PetscStackCallStandard(hypre_BoomerAMGBuildCoarseOperator,(hR,hA,hP,hRAP));
763:   PetscStackCallStandard(hypre_ParCSRMatrixSetNumNonzeros,(*hRAP));
764:   /* hypre_BoomerAMGBuildCoarseOperator steals the col_starts from P and the row_starts from R */
765:   hypre_ParCSRMatrixSetRowStartsOwner(*hRAP,0);
766:   hypre_ParCSRMatrixSetColStartsOwner(*hRAP,0);
767:   if (P_owns_col_starts) hypre_ParCSRMatrixSetColStartsOwner(hP,1);
768:   if (R_owns_row_starts) hypre_ParCSRMatrixSetRowStartsOwner(hR,1);
769:   return(0);
770: }

772: static PetscErrorCode MatPtAPNumeric_AIJ_AIJ_wHYPRE(Mat A,Mat P,Mat C)
773: {
774:   Mat                B;
775:   hypre_ParCSRMatrix *hA,*hP,*hPtAP;
776:   PetscErrorCode     ierr;
777:   Mat_Product        *product=C->product;

780:   MatAIJGetParCSR_Private(A,&hA);
781:   MatAIJGetParCSR_Private(P,&hP);
782:   MatHYPRE_ParCSR_RAP(hP,hA,hP,&hPtAP);
783:   MatCreateFromParCSR(hPtAP,MATAIJ,PETSC_OWN_POINTER,&B);

785:   MatHeaderMerge(C,&B);
786:   C->product = product;

788:   MatAIJRestoreParCSR_Private(A,&hA);
789:   MatAIJRestoreParCSR_Private(P,&hP);
790:   return(0);
791: }

793: PETSC_INTERN PetscErrorCode MatPtAPSymbolic_AIJ_AIJ_wHYPRE(Mat A,Mat P,PetscReal fill,Mat C)
794: {

798:   MatSetType(C,MATAIJ);
799:   C->ops->ptapnumeric    = MatPtAPNumeric_AIJ_AIJ_wHYPRE;
800:   C->ops->productnumeric = MatProductNumeric_PtAP;
801:   return(0);
802: }

804: static PetscErrorCode MatPtAPNumeric_AIJ_HYPRE(Mat A,Mat P,Mat C)
805: {
806:   Mat                B;
807:   Mat_HYPRE          *hP;
808:   hypre_ParCSRMatrix *hA = NULL,*Pparcsr,*ptapparcsr;
809:   HYPRE_Int          type;
810:   MPI_Comm           comm = PetscObjectComm((PetscObject)A);
811:   PetscBool          ishypre;
812:   PetscErrorCode     ierr;

815:   PetscObjectTypeCompare((PetscObject)P,MATHYPRE,&ishypre);
816:   if (!ishypre) SETERRQ1(comm,PETSC_ERR_USER,"P should be of type %s",MATHYPRE);
817:   hP = (Mat_HYPRE*)P->data;
818:   PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(hP->ij,&type));
819:   if (type != HYPRE_PARCSR) SETERRQ(comm,PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported");
820:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hP->ij,(void**)&Pparcsr));

822:   MatAIJGetParCSR_Private(A,&hA);
823:   MatHYPRE_ParCSR_RAP(Pparcsr,hA,Pparcsr,&ptapparcsr);
824:   MatAIJRestoreParCSR_Private(A,&hA);

826:   /* create temporary matrix and merge to C */
827:   MatCreateFromParCSR(ptapparcsr,((PetscObject)C)->type_name,PETSC_OWN_POINTER,&B);
828:   MatHeaderMerge(C,&B);
829:   return(0);
830: }

832: static PetscErrorCode MatPtAPNumeric_HYPRE_HYPRE(Mat A,Mat P,Mat C)
833: {
834:   Mat                B;
835:   hypre_ParCSRMatrix *Aparcsr,*Pparcsr,*ptapparcsr;
836:   Mat_HYPRE          *hA,*hP;
837:   PetscBool          ishypre;
838:   HYPRE_Int          type;
839:   PetscErrorCode     ierr;

842:   PetscObjectTypeCompare((PetscObject)P,MATHYPRE,&ishypre);
843:   if (!ishypre) SETERRQ1(PetscObjectComm((PetscObject)P),PETSC_ERR_USER,"P should be of type %s",MATHYPRE);
844:   PetscObjectTypeCompare((PetscObject)A,MATHYPRE,&ishypre);
845:   if (!ishypre) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"A should be of type %s",MATHYPRE);
846:   hA = (Mat_HYPRE*)A->data;
847:   hP = (Mat_HYPRE*)P->data;
848:   PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(hA->ij,&type));
849:   if (type != HYPRE_PARCSR) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported");
850:   PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(hP->ij,&type));
851:   if (type != HYPRE_PARCSR) SETERRQ(PetscObjectComm((PetscObject)P),PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported");
852:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hA->ij,(void**)&Aparcsr));
853:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hP->ij,(void**)&Pparcsr));
854:   MatHYPRE_ParCSR_RAP(Pparcsr,Aparcsr,Pparcsr,&ptapparcsr);
855:   MatCreateFromParCSR(ptapparcsr,MATHYPRE,PETSC_OWN_POINTER,&B);
856:   MatHeaderMerge(C,&B);
857:   return(0);
858: }

860: /* calls hypre_ParMatmul
861:    hypre_ParMatMul uses hypre_ParMatrixCreate with the communicator of hA
862:    hypre_ParMatrixCreate does not duplicate the communicator
863:    It looks like we don't need to have the diagonal entries
864:    ordered first in the rows of the diagonal part
865:    for boomerAMGBuildCoarseOperator to work */
866: static PetscErrorCode MatHYPRE_ParCSR_MatMatMult(hypre_ParCSRMatrix *hA, hypre_ParCSRMatrix *hB, hypre_ParCSRMatrix **hAB)
867: {
869:   PetscStackPush("hypre_ParMatmul");
870:   *hAB = hypre_ParMatmul(hA,hB);
871:   PetscStackPop;
872:   return(0);
873: }

875: static PetscErrorCode MatMatMultNumeric_AIJ_AIJ_wHYPRE(Mat A,Mat B,Mat C)
876: {
877:   Mat                D;
878:   hypre_ParCSRMatrix *hA,*hB,*hAB = NULL;
879:   PetscErrorCode     ierr;
880:   Mat_Product        *product=C->product;

883:   MatAIJGetParCSR_Private(A,&hA);
884:   MatAIJGetParCSR_Private(B,&hB);
885:   MatHYPRE_ParCSR_MatMatMult(hA,hB,&hAB);
886:   MatCreateFromParCSR(hAB,MATAIJ,PETSC_OWN_POINTER,&D);

888:   MatHeaderMerge(C,&D);
889:   C->product = product;

891:   MatAIJRestoreParCSR_Private(A,&hA);
892:   MatAIJRestoreParCSR_Private(B,&hB);
893:   return(0);
894: }

896: PETSC_INTERN PetscErrorCode MatMatMultSymbolic_AIJ_AIJ_wHYPRE(Mat A,Mat B,PetscReal fill,Mat C)
897: {

901:   MatSetType(C,MATAIJ);
902:   C->ops->matmultnumeric = MatMatMultNumeric_AIJ_AIJ_wHYPRE;
903:   C->ops->productnumeric = MatProductNumeric_AB;
904:   return(0);
905: }

907: static PetscErrorCode MatMatMultNumeric_HYPRE_HYPRE(Mat A,Mat B,Mat C)
908: {
909:   Mat                D;
910:   hypre_ParCSRMatrix *Aparcsr,*Bparcsr,*ABparcsr = NULL;
911:   Mat_HYPRE          *hA,*hB;
912:   PetscBool          ishypre;
913:   HYPRE_Int          type;
914:   PetscErrorCode     ierr;
915:   Mat_Product        *product;

918:   PetscObjectTypeCompare((PetscObject)B,MATHYPRE,&ishypre);
919:   if (!ishypre) SETERRQ1(PetscObjectComm((PetscObject)B),PETSC_ERR_USER,"B should be of type %s",MATHYPRE);
920:   PetscObjectTypeCompare((PetscObject)A,MATHYPRE,&ishypre);
921:   if (!ishypre) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"A should be of type %s",MATHYPRE);
922:   hA = (Mat_HYPRE*)A->data;
923:   hB = (Mat_HYPRE*)B->data;
924:   PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(hA->ij,&type));
925:   if (type != HYPRE_PARCSR) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported");
926:   PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(hB->ij,&type));
927:   if (type != HYPRE_PARCSR) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported");
928:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hA->ij,(void**)&Aparcsr));
929:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hB->ij,(void**)&Bparcsr));
930:   MatHYPRE_ParCSR_MatMatMult(Aparcsr,Bparcsr,&ABparcsr);
931:   MatCreateFromParCSR(ABparcsr,MATHYPRE,PETSC_OWN_POINTER,&D);

933:   /* need to use HeaderReplace because HeaderMerge messes up with the communicator */
934:   product    = C->product;  /* save it from MatHeaderReplace() */
935:   C->product = NULL;
936:   MatHeaderReplace(C,&D);
937:   C->product = product;
938:   C->ops->matmultnumeric = MatMatMultNumeric_HYPRE_HYPRE;
939:   C->ops->productnumeric = MatProductNumeric_AB;
940:   return(0);
941: }

943: PETSC_INTERN PetscErrorCode MatTransposeMatMatMultNumeric_AIJ_AIJ_AIJ_wHYPRE(Mat A,Mat B,Mat C,Mat D)
944: {
945:   Mat                E;
946:   hypre_ParCSRMatrix *hA,*hB,*hC,*hABC;
947:   PetscErrorCode     ierr;

950:   MatAIJGetParCSR_Private(A,&hA);
951:   MatAIJGetParCSR_Private(B,&hB);
952:   MatAIJGetParCSR_Private(C,&hC);
953:   MatHYPRE_ParCSR_RAP(hA,hB,hC,&hABC);
954:   MatCreateFromParCSR(hABC,MATAIJ,PETSC_OWN_POINTER,&E);
955:   MatHeaderMerge(D,&E);
956:   MatAIJRestoreParCSR_Private(A,&hA);
957:   MatAIJRestoreParCSR_Private(B,&hB);
958:   MatAIJRestoreParCSR_Private(C,&hC);
959:   return(0);
960: }

962: PETSC_INTERN PetscErrorCode MatTransposeMatMatMultSymbolic_AIJ_AIJ_AIJ_wHYPRE(Mat A,Mat B,Mat C,PetscReal fill,Mat D)
963: {

967:   MatSetType(D,MATAIJ);
968:   return(0);
969: }

971: /* ---------------------------------------------------- */
972: static PetscErrorCode MatProductSymbolic_AB_HYPRE(Mat C)
973: {
975:   C->ops->productnumeric = MatProductNumeric_AB;
976:   return(0);
977: }

979: static PetscErrorCode MatProductSetFromOptions_HYPRE_AB(Mat C)
980: {
982:   Mat_Product    *product = C->product;
983:   PetscBool      Ahypre;

986:   PetscObjectTypeCompare((PetscObject)product->A,MATHYPRE,&Ahypre);
987:   if (Ahypre) { /* A is a Hypre matrix */
988:     MatSetType(C,MATHYPRE);
989:     C->ops->productsymbolic = MatProductSymbolic_AB_HYPRE;
990:     C->ops->matmultnumeric  = MatMatMultNumeric_HYPRE_HYPRE;
991:     return(0);
992:   }
993:   return(0);
994: }

996: static PetscErrorCode MatProductSymbolic_PtAP_HYPRE(Mat C)
997: {
999:   C->ops->productnumeric = MatProductNumeric_PtAP;
1000:   return(0);
1001: }

1003: static PetscErrorCode MatProductSetFromOptions_HYPRE_PtAP(Mat C)
1004: {
1006:   Mat_Product    *product = C->product;
1007:   PetscBool      flg;
1008:   PetscInt       type = 0;
1009:   const char     *outTypes[4] = {"aij","seqaij","mpiaij","hypre"};
1010:   PetscInt       ntype = 4;
1011:   Mat            A = product->A;
1012:   PetscBool      Ahypre;

1015:   PetscObjectTypeCompare((PetscObject)A,MATHYPRE,&Ahypre);
1016:   if (Ahypre) { /* A is a Hypre matrix */
1017:     MatSetType(C,MATHYPRE);
1018:     C->ops->productsymbolic = MatProductSymbolic_PtAP_HYPRE;
1019:     C->ops->ptapnumeric     = MatPtAPNumeric_HYPRE_HYPRE;
1020:     return(0);
1021:   }

1023:   /* A is AIJ, P is Hypre, C = PtAP can be either AIJ or Hypre format */
1024:   /* Get runtime option */
1025:   if (product->api_user) {
1026:     PetscOptionsBegin(PetscObjectComm((PetscObject)C),((PetscObject)C)->prefix,"MatPtAP_HYPRE","Mat");
1027:     PetscOptionsEList("-matptap_hypre_outtype","MatPtAP outtype","MatPtAP outtype",outTypes,ntype,outTypes[type],&type,&flg);
1028:     PetscOptionsEnd();
1029:   } else {
1030:     PetscOptionsBegin(PetscObjectComm((PetscObject)C),((PetscObject)C)->prefix,"MatProduct_PtAP_HYPRE","Mat");
1031:     PetscOptionsEList("-matproduct_ptap_hypre_outtype","MatProduct_PtAP outtype","MatProduct_PtAP",outTypes,ntype,outTypes[type],&type,&flg);
1032:     PetscOptionsEnd();
1033:   }

1035:   if (type == 0 || type == 1 || type == 2) {
1036:     MatSetType(C,MATAIJ);
1037:   } else if (type == 3) {
1038:     MatSetType(C,MATHYPRE);
1039:   } else SETERRQ(PetscObjectComm((PetscObject)C),PETSC_ERR_SUP,"MatPtAP outtype is not supported");
1040:   C->ops->productsymbolic = MatProductSymbolic_PtAP_HYPRE;
1041:   C->ops->ptapnumeric     = MatPtAPNumeric_AIJ_HYPRE;
1042:   return(0);
1043: }

1045: static PetscErrorCode MatProductSetFromOptions_HYPRE(Mat C)
1046: {
1048:   Mat_Product    *product = C->product;

1051:   switch (product->type) {
1052:   case MATPRODUCT_AB:
1053:     MatProductSetFromOptions_HYPRE_AB(C);
1054:     break;
1055:   case MATPRODUCT_PtAP:
1056:     MatProductSetFromOptions_HYPRE_PtAP(C);
1057:     break;
1058:   default:
1059:     break;
1060:   }
1061:   return(0);
1062: }

1064: /* -------------------------------------------------------- */

1066: static PetscErrorCode MatMultTranspose_HYPRE(Mat A, Vec x, Vec y)
1067: {

1071:   MatHYPRE_MultKernel_Private(A,1.0,x,0.0,y,PETSC_TRUE);
1072:   return(0);
1073: }

1075: static PetscErrorCode MatMult_HYPRE(Mat A, Vec x, Vec y)
1076: {

1080:   MatHYPRE_MultKernel_Private(A,1.0,x,0.0,y,PETSC_FALSE);
1081:   return(0);
1082: }

1084: static PetscErrorCode MatMultAdd_HYPRE(Mat A, Vec x, Vec y, Vec z)
1085: {

1089:   if (y != z) {
1090:     VecCopy(y,z);
1091:   }
1092:   MatHYPRE_MultKernel_Private(A,1.0,x,1.0,z,PETSC_FALSE);
1093:   return(0);
1094: }

1096: static PetscErrorCode MatMultTransposeAdd_HYPRE(Mat A, Vec x, Vec y, Vec z)
1097: {

1101:   if (y != z) {
1102:     VecCopy(y,z);
1103:   }
1104:   MatHYPRE_MultKernel_Private(A,1.0,x,1.0,z,PETSC_TRUE);
1105:   return(0);
1106: }

1108: /* y = a * A * x + b * y or y = a * A^t * x + b * y depending on trans */
1109: static PetscErrorCode MatHYPRE_MultKernel_Private(Mat A, HYPRE_Complex a, Vec x, HYPRE_Complex b, Vec y, PetscBool trans)
1110: {
1111:   Mat_HYPRE          *hA = (Mat_HYPRE*)A->data;
1112:   hypre_ParCSRMatrix *parcsr;
1113:   hypre_ParVector    *hx,*hy;
1114:   HYPRE_Complex      *ax,*ay,*sax,*say;
1115:   PetscErrorCode     ierr;

1118:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hA->ij,(void**)&parcsr));
1119:   PetscStackCallStandard(HYPRE_IJVectorGetObject,(hA->x,(void**)&hx));
1120:   PetscStackCallStandard(HYPRE_IJVectorGetObject,(hA->b,(void**)&hy));
1121:   VecGetArrayRead(x,(const PetscScalar**)&ax);
1122:   VecGetArray(y,(PetscScalar**)&ay);
1123:   if (trans) {
1124:     VecHYPRE_ParVectorReplacePointer(hA->x,ay,say);
1125:     VecHYPRE_ParVectorReplacePointer(hA->b,ax,sax);
1126:     hypre_ParCSRMatrixMatvecT(a,parcsr,hy,b,hx);
1127:     VecHYPRE_ParVectorReplacePointer(hA->x,say,ay);
1128:     VecHYPRE_ParVectorReplacePointer(hA->b,sax,ax);
1129:   } else {
1130:     VecHYPRE_ParVectorReplacePointer(hA->x,ax,sax);
1131:     VecHYPRE_ParVectorReplacePointer(hA->b,ay,say);
1132:     hypre_ParCSRMatrixMatvec(a,parcsr,hx,b,hy);
1133:     VecHYPRE_ParVectorReplacePointer(hA->x,sax,ax);
1134:     VecHYPRE_ParVectorReplacePointer(hA->b,say,ay);
1135:   }
1136:   VecRestoreArrayRead(x,(const PetscScalar**)&ax);
1137:   VecRestoreArray(y,(PetscScalar**)&ay);
1138:   return(0);
1139: }

1141: static PetscErrorCode MatDestroy_HYPRE(Mat A)
1142: {
1143:   Mat_HYPRE      *hA = (Mat_HYPRE*)A->data;

1147:   if (hA->x) PetscStackCallStandard(HYPRE_IJVectorDestroy,(hA->x));
1148:   if (hA->b) PetscStackCallStandard(HYPRE_IJVectorDestroy,(hA->b));
1149:   if (hA->ij) {
1150:     if (!hA->inner_free) hypre_IJMatrixObject(hA->ij) = NULL;
1151:     PetscStackCallStandard(HYPRE_IJMatrixDestroy,(hA->ij));
1152:   }
1153:   if (hA->comm) { MPI_Comm_free(&hA->comm); }

1155:   MatStashDestroy_Private(&A->stash);

1157:   PetscFree(hA->array);

1159:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_hypre_aij_C",NULL);
1160:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_hypre_is_C",NULL);
1161:   PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_seqaij_hypre_C",NULL);
1162:   PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_mpiaij_hypre_C",NULL);
1163:   PetscObjectComposeFunction((PetscObject)A,"MatHYPRESetPreallocation_C",NULL);
1164:   PetscObjectComposeFunction((PetscObject)A,"MatHYPREGetParCSR_C",NULL);
1165:   PetscFree(A->data);
1166:   return(0);
1167: }

1169: static PetscErrorCode MatSetUp_HYPRE(Mat A)
1170: {

1174:   MatHYPRESetPreallocation(A,PETSC_DEFAULT,NULL,PETSC_DEFAULT,NULL);
1175:   return(0);
1176: }

1178: static PetscErrorCode MatAssemblyEnd_HYPRE(Mat A, MatAssemblyType mode)
1179: {
1180:   Mat_HYPRE          *hA = (Mat_HYPRE*)A->data;
1181:   Vec                x,b;
1182:   PetscMPIInt        n;
1183:   PetscInt           i,j,rstart,ncols,flg;
1184:   PetscInt           *row,*col;
1185:   PetscScalar        *val;
1186:   PetscErrorCode     ierr;

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

1191:   if (!A->nooffprocentries) {
1192:     while (1) {
1193:       MatStashScatterGetMesg_Private(&A->stash,&n,&row,&col,&val,&flg);
1194:       if (!flg) break;

1196:       for (i=0; i<n;) {
1197:         /* Now identify the consecutive vals belonging to the same row */
1198:         for (j=i,rstart=row[j]; j<n; j++) {
1199:           if (row[j] != rstart) break;
1200:         }
1201:         if (j < n) ncols = j-i;
1202:         else       ncols = n-i;
1203:         /* Now assemble all these values with a single function call */
1204:         MatSetValues_HYPRE(A,1,row+i,ncols,col+i,val+i,A->insertmode);

1206:         i = j;
1207:       }
1208:     }
1209:     MatStashScatterEnd_Private(&A->stash);
1210:   }

1212:   PetscStackCallStandard(HYPRE_IJMatrixAssemble,(hA->ij));
1213:   /* The assembly routine destroys the aux_matrix, we recreate it here by calling HYPRE_IJMatrixInitialize */
1214:   /* 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 */
1215:   if (!hA->sorted_full) {
1216:     hypre_AuxParCSRMatrix *aux_matrix;

1218:     /* call destroy just to make sure we do not leak anything */
1219:     aux_matrix = (hypre_AuxParCSRMatrix*)hypre_IJMatrixTranslator(hA->ij);
1220:     PetscStackCallStandard(hypre_AuxParCSRMatrixDestroy,(aux_matrix));
1221:     hypre_IJMatrixTranslator(hA->ij) = NULL;

1223:     /* Initialize with assembled flag -> it only recreates the aux_par_matrix */
1224:     PetscStackCallStandard(HYPRE_IJMatrixInitialize,(hA->ij));
1225:     aux_matrix = (hypre_AuxParCSRMatrix*)hypre_IJMatrixTranslator(hA->ij);
1226:     hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 1; /* see comment in MatHYPRESetPreallocation_HYPRE */
1227: #if PETSC_PKG_HYPRE_VERSION_LT(2,19,0)
1228:     PetscStackCallStandard(hypre_AuxParCSRMatrixInitialize,(aux_matrix));
1229: #else
1230:     PetscStackCallStandard(hypre_AuxParCSRMatrixInitialize_v2,(aux_matrix,HYPRE_MEMORY_HOST));
1231: #endif
1232:   }
1233:   if (hA->x) return(0);
1234:   PetscLayoutSetUp(A->rmap);
1235:   PetscLayoutSetUp(A->cmap);
1236:   VecCreateMPIWithArray(PetscObjectComm((PetscObject)A),1,A->cmap->n,A->cmap->N,NULL,&x);
1237:   VecCreateMPIWithArray(PetscObjectComm((PetscObject)A),1,A->rmap->n,A->rmap->N,NULL,&b);
1238:   VecHYPRE_IJVectorCreate(x,&hA->x);
1239:   VecHYPRE_IJVectorCreate(b,&hA->b);
1240:   VecDestroy(&x);
1241:   VecDestroy(&b);
1242:   return(0);
1243: }

1245: static PetscErrorCode MatGetArray_HYPRE(Mat A, PetscInt size, void **array)
1246: {
1247:   Mat_HYPRE          *hA = (Mat_HYPRE*)A->data;
1248:   PetscErrorCode     ierr;

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

1253:   if (hA->size >= size) {
1254:     *array = hA->array;
1255:   } else {
1256:     PetscFree(hA->array);
1257:     hA->size = size;
1258:     PetscMalloc(hA->size,&hA->array);
1259:     *array = hA->array;
1260:   }

1262:   hA->available = PETSC_FALSE;
1263:   return(0);
1264: }

1266: static PetscErrorCode MatRestoreArray_HYPRE(Mat A, void **array)
1267: {
1268:   Mat_HYPRE          *hA = (Mat_HYPRE*)A->data;

1271:   *array = NULL;
1272:   hA->available = PETSC_TRUE;
1273:   return(0);
1274: }

1276: PetscErrorCode MatSetValues_HYPRE(Mat A, PetscInt nr, const PetscInt rows[], PetscInt nc, const PetscInt cols[], const PetscScalar v[], InsertMode ins)
1277: {
1278:   Mat_HYPRE          *hA = (Mat_HYPRE*)A->data;
1279:   PetscScalar        *vals = (PetscScalar *)v;
1280:   HYPRE_Complex      *sscr;
1281:   PetscInt           *cscr[2];
1282:   PetscInt           i,nzc;
1283:   void               *array = NULL;
1284:   PetscErrorCode     ierr;

1287:   MatGetArray_HYPRE(A,sizeof(PetscInt)*(2*nc)+sizeof(HYPRE_Complex)*nc*nr,&array);
1288:   cscr[0] = (PetscInt*)array;
1289:   cscr[1] = ((PetscInt*)array)+nc;
1290:   sscr = (HYPRE_Complex*)(((PetscInt*)array)+nc*2);
1291:   for (i=0,nzc=0;i<nc;i++) {
1292:     if (cols[i] >= 0) {
1293:       cscr[0][nzc  ] = cols[i];
1294:       cscr[1][nzc++] = i;
1295:     }
1296:   }
1297:   if (!nzc) {
1298:     MatRestoreArray_HYPRE(A,&array);
1299:     return(0);
1300:   }

1302:   if (ins == ADD_VALUES) {
1303:     for (i=0;i<nr;i++) {
1304:       if (rows[i] >= 0 && nzc) {
1305:         PetscInt  j;
1306:         HYPRE_Int hnc = (HYPRE_Int)nzc;

1308:         if ((PetscInt)hnc != nzc) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Hypre overflow! number of columns %D for row %D",nzc,rows[i]);
1309:         for (j=0;j<nzc;j++) { PetscHYPREScalarCast(vals[cscr[1][j]],&sscr[j]); }
1310:         PetscStackCallStandard(HYPRE_IJMatrixAddToValues,(hA->ij,1,&hnc,(HYPRE_BigInt*)(rows+i),(HYPRE_BigInt*)cscr[0],sscr));
1311:       }
1312:       vals += nc;
1313:     }
1314:   } else { /* INSERT_VALUES */
1315:     PetscInt rst,ren;

1317:     MatGetOwnershipRange(A,&rst,&ren);
1318:     for (i=0;i<nr;i++) {
1319:       if (rows[i] >= 0 && nzc) {
1320:         PetscInt  j;
1321:         HYPRE_Int hnc = (HYPRE_Int)nzc;

1323:         if ((PetscInt)hnc != nzc) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Hypre overflow! number of columns %D for row %D",nzc,rows[i]);
1324:         for (j=0;j<nzc;j++) { PetscHYPREScalarCast(vals[cscr[1][j]],&sscr[j]); }
1325:         /* nonlocal values */
1326:         if (rows[i] < rst || rows[i] >= ren) { MatStashValuesRow_Private(&A->stash,rows[i],nzc,cscr[0],(PetscScalar*)sscr,PETSC_FALSE); }
1327:         /* local values */
1328:         else PetscStackCallStandard(HYPRE_IJMatrixSetValues,(hA->ij,1,&hnc,(HYPRE_BigInt*)(rows+i),(HYPRE_BigInt*)cscr[0],sscr));
1329:       }
1330:       vals += nc;
1331:     }
1332:   }

1334:   MatRestoreArray_HYPRE(A,&array);
1335:   return(0);
1336: }

1338: static PetscErrorCode MatHYPRESetPreallocation_HYPRE(Mat A, PetscInt dnz, const PetscInt dnnz[], PetscInt onz, const PetscInt onnz[])
1339: {
1340:   Mat_HYPRE      *hA = (Mat_HYPRE*)A->data;
1341:   HYPRE_Int      *hdnnz,*honnz;
1342:   PetscInt       i,rs,re,cs,ce,bs;
1343:   PetscMPIInt    size;

1347:   MatGetBlockSize(A,&bs);
1348:   PetscLayoutSetUp(A->rmap);
1349:   PetscLayoutSetUp(A->cmap);
1350:   rs   = A->rmap->rstart;
1351:   re   = A->rmap->rend;
1352:   cs   = A->cmap->rstart;
1353:   ce   = A->cmap->rend;
1354:   if (!hA->ij) {
1355:     PetscStackCallStandard(HYPRE_IJMatrixCreate,(hA->comm,rs,re-1,cs,ce-1,&hA->ij));
1356:     PetscStackCallStandard(HYPRE_IJMatrixSetObjectType,(hA->ij,HYPRE_PARCSR));
1357:   } else {
1358:     HYPRE_BigInt hrs,hre,hcs,hce;
1359:     PetscStackCallStandard(HYPRE_IJMatrixGetLocalRange,(hA->ij,&hrs,&hre,&hcs,&hce));
1360:     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);
1361:     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);
1362:   }
1363:   if (dnz == PETSC_DEFAULT || dnz == PETSC_DECIDE) dnz = 10*bs;
1364:   if (onz == PETSC_DEFAULT || onz == PETSC_DECIDE) onz = 10*bs;

1366:   if (!dnnz) {
1367:     PetscMalloc1(A->rmap->n,&hdnnz);
1368:     for (i=0;i<A->rmap->n;i++) hdnnz[i] = dnz;
1369:   } else {
1370:     hdnnz = (HYPRE_Int*)dnnz;
1371:   }
1372:   MPI_Comm_size(PetscObjectComm((PetscObject)A),&size);
1373:   if (size > 1) {
1374:     hypre_AuxParCSRMatrix *aux_matrix;
1375:     if (!onnz) {
1376:       PetscMalloc1(A->rmap->n,&honnz);
1377:       for (i=0;i<A->rmap->n;i++) honnz[i] = onz;
1378:     } else honnz = (HYPRE_Int*)onnz;
1379:     /* SetDiagOffdSizes sets hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 0, since it seems
1380:        they assume the user will input the entire row values, properly sorted
1381:        In PETSc, we don't make such an assumption and set this flag to 1,
1382:        unless the option MAT_SORTED_FULL is set to true.
1383:        Also, to avoid possible memory leaks, we destroy and recreate the translator
1384:        This has to be done here, as HYPRE_IJMatrixInitialize will properly initialize
1385:        the IJ matrix for us */
1386:     aux_matrix = (hypre_AuxParCSRMatrix*)hypre_IJMatrixTranslator(hA->ij);
1387:     hypre_AuxParCSRMatrixDestroy(aux_matrix);
1388:     hypre_IJMatrixTranslator(hA->ij) = NULL;
1389:     PetscStackCallStandard(HYPRE_IJMatrixSetDiagOffdSizes,(hA->ij,hdnnz,honnz));
1390:     aux_matrix = (hypre_AuxParCSRMatrix*)hypre_IJMatrixTranslator(hA->ij);
1391:     hypre_AuxParCSRMatrixNeedAux(aux_matrix) = !hA->sorted_full;
1392:   } else {
1393:     honnz = NULL;
1394:     PetscStackCallStandard(HYPRE_IJMatrixSetRowSizes,(hA->ij,hdnnz));
1395:   }

1397:   /* reset assembled flag and call the initialize method */
1398:   hypre_IJMatrixAssembleFlag(hA->ij) = 0;
1399:   PetscStackCallStandard(HYPRE_IJMatrixInitialize,(hA->ij));
1400:   if (!dnnz) {
1401:     PetscFree(hdnnz);
1402:   }
1403:   if (!onnz && honnz) {
1404:     PetscFree(honnz);
1405:   }

1407:   /* Match AIJ logic */
1408:   A->preallocated = PETSC_TRUE;
1409:   A->assembled    = PETSC_FALSE;
1410:   return(0);
1411: }

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

1416:    Collective on Mat

1418:    Input Parameters:
1419: +  A - the matrix
1420: .  dnz  - number of nonzeros per row in DIAGONAL portion of local submatrix
1421:           (same value is used for all local rows)
1422: .  dnnz - array containing the number of nonzeros in the various rows of the
1423:           DIAGONAL portion of the local submatrix (possibly different for each row)
1424:           or NULL (PETSC_NULL_INTEGER in Fortran), if d_nz is used to specify the nonzero structure.
1425:           The size of this array is equal to the number of local rows, i.e 'm'.
1426:           For matrices that will be factored, you must leave room for (and set)
1427:           the diagonal entry even if it is zero.
1428: .  onz  - number of nonzeros per row in the OFF-DIAGONAL portion of local
1429:           submatrix (same value is used for all local rows).
1430: -  onnz - array containing the number of nonzeros in the various rows of the
1431:           OFF-DIAGONAL portion of the local submatrix (possibly different for
1432:           each row) or NULL (PETSC_NULL_INTEGER in Fortran), if o_nz is used to specify the nonzero
1433:           structure. The size of this array is equal to the number
1434:           of local rows, i.e 'm'.

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

1439:    Level: intermediate

1441: .seealso: MatCreate(), MatMPIAIJSetPreallocation(), MATHYPRE
1442: @*/
1443: PetscErrorCode MatHYPRESetPreallocation(Mat A, PetscInt dnz, const PetscInt dnnz[], PetscInt onz, const PetscInt onnz[])
1444: {

1450:   PetscTryMethod(A,"MatHYPRESetPreallocation_C",(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt[]),(A,dnz,dnnz,onz,onnz));
1451:   return(0);
1452: }

1454: /*
1455:    MatCreateFromParCSR - Creates a matrix from a hypre_ParCSRMatrix

1457:    Collective

1459:    Input Parameters:
1460: +  parcsr   - the pointer to the hypre_ParCSRMatrix
1461: .  mtype    - matrix type to be created. Currently MATAIJ, MATIS and MATHYPRE are supported.
1462: -  copymode - PETSc copying options

1464:    Output Parameter:
1465: .  A  - the matrix

1467:    Level: intermediate

1469: .seealso: MatHYPRE, PetscCopyMode
1470: */
1471: PETSC_EXTERN PetscErrorCode MatCreateFromParCSR(hypre_ParCSRMatrix *parcsr, MatType mtype, PetscCopyMode copymode, Mat* A)
1472: {
1473:   Mat                   T;
1474:   Mat_HYPRE             *hA;
1475:   MPI_Comm              comm;
1476:   PetscInt              rstart,rend,cstart,cend,M,N;
1477:   PetscBool             isseqaij,isseqaijmkl,ismpiaij,isaij,ishyp,isis;
1478:   PetscErrorCode        ierr;

1481:   comm   = hypre_ParCSRMatrixComm(parcsr);
1482:   PetscStrcmp(mtype,MATSEQAIJ,&isseqaij);
1483:   PetscStrcmp(mtype,MATSEQAIJMKL,&isseqaijmkl);
1484:   PetscStrcmp(mtype,MATMPIAIJ,&ismpiaij);
1485:   PetscStrcmp(mtype,MATAIJ,&isaij);
1486:   PetscStrcmp(mtype,MATHYPRE,&ishyp);
1487:   PetscStrcmp(mtype,MATIS,&isis);
1488:   isaij  = (PetscBool)(isseqaij || isseqaijmkl || ismpiaij || isaij);
1489:   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);
1490:   /* access ParCSRMatrix */
1491:   rstart = hypre_ParCSRMatrixFirstRowIndex(parcsr);
1492:   rend   = hypre_ParCSRMatrixLastRowIndex(parcsr);
1493:   cstart = hypre_ParCSRMatrixFirstColDiag(parcsr);
1494:   cend   = hypre_ParCSRMatrixLastColDiag(parcsr);
1495:   M      = hypre_ParCSRMatrixGlobalNumRows(parcsr);
1496:   N      = hypre_ParCSRMatrixGlobalNumCols(parcsr);

1498:   /* fix for empty local rows/columns */
1499:   if (rend < rstart) rend = rstart;
1500:   if (cend < cstart) cend = cstart;

1502:   /* PETSc convention */
1503:   rend++;
1504:   cend++;
1505:   rend = PetscMin(rend,M);
1506:   cend = PetscMin(cend,N);

1508:   /* create PETSc matrix with MatHYPRE */
1509:   MatCreate(comm,&T);
1510:   MatSetSizes(T,rend-rstart,cend-cstart,M,N);
1511:   MatSetType(T,MATHYPRE);
1512:   hA   = (Mat_HYPRE*)(T->data);

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

1518:   /* create new ParCSR object if needed */
1519:   if (ishyp && copymode == PETSC_COPY_VALUES) {
1520:     hypre_ParCSRMatrix *new_parcsr;
1521:     hypre_CSRMatrix    *hdiag,*hoffd,*ndiag,*noffd;

1523:     new_parcsr = hypre_ParCSRMatrixClone(parcsr,0);
1524:     hdiag      = hypre_ParCSRMatrixDiag(parcsr);
1525:     hoffd      = hypre_ParCSRMatrixOffd(parcsr);
1526:     ndiag      = hypre_ParCSRMatrixDiag(new_parcsr);
1527:     noffd      = hypre_ParCSRMatrixOffd(new_parcsr);
1528:     PetscArraycpy(hypre_CSRMatrixData(ndiag),hypre_CSRMatrixData(hdiag),hypre_CSRMatrixNumNonzeros(hdiag));
1529:     PetscArraycpy(hypre_CSRMatrixData(noffd),hypre_CSRMatrixData(hoffd),hypre_CSRMatrixNumNonzeros(hoffd));
1530:     parcsr     = new_parcsr;
1531:     copymode   = PETSC_OWN_POINTER;
1532:   }

1534:   /* set ParCSR object */
1535:   hypre_IJMatrixObject(hA->ij) = parcsr;
1536:   T->preallocated = PETSC_TRUE;

1538:   /* set assembled flag */
1539:   hypre_IJMatrixAssembleFlag(hA->ij) = 1;
1540:   PetscStackCallStandard(HYPRE_IJMatrixInitialize,(hA->ij));
1541:   if (ishyp) {
1542:     PetscMPIInt myid = 0;

1544:     /* make sure we always have row_starts and col_starts available */
1545:     if (HYPRE_AssumedPartitionCheck()) {
1546:       MPI_Comm_rank(comm,&myid);
1547:     }
1548:     if (!hypre_ParCSRMatrixOwnsColStarts(parcsr)) {
1549:       PetscLayout map;

1551:       MatGetLayouts(T,NULL,&map);
1552:       PetscLayoutSetUp(map);
1553:       hypre_ParCSRMatrixColStarts(parcsr) = (HYPRE_BigInt*)(map->range + myid);
1554:     }
1555:     if (!hypre_ParCSRMatrixOwnsRowStarts(parcsr)) {
1556:       PetscLayout map;

1558:       MatGetLayouts(T,&map,NULL);
1559:       PetscLayoutSetUp(map);
1560:       hypre_ParCSRMatrixRowStarts(parcsr) = (HYPRE_BigInt*)(map->range + myid);
1561:     }
1562:     /* prevent from freeing the pointer */
1563:     if (copymode == PETSC_USE_POINTER) hA->inner_free = PETSC_FALSE;
1564:     *A   = T;
1565:     MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);
1566:     MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);
1567:   } else if (isaij) {
1568:     if (copymode != PETSC_OWN_POINTER) {
1569:       /* prevent from freeing the pointer */
1570:       hA->inner_free = PETSC_FALSE;
1571:       MatConvert_HYPRE_AIJ(T,MATAIJ,MAT_INITIAL_MATRIX,A);
1572:       MatDestroy(&T);
1573:     } else { /* AIJ return type with PETSC_OWN_POINTER */
1574:       MatConvert_HYPRE_AIJ(T,MATAIJ,MAT_INPLACE_MATRIX,&T);
1575:       *A   = T;
1576:     }
1577:   } else if (isis) {
1578:     MatConvert_HYPRE_IS(T,MATIS,MAT_INITIAL_MATRIX,A);
1579:     if (copymode != PETSC_OWN_POINTER) hA->inner_free = PETSC_FALSE;
1580:     MatDestroy(&T);
1581:   }
1582:   return(0);
1583: }

1585: static PetscErrorCode MatHYPREGetParCSR_HYPRE(Mat A, hypre_ParCSRMatrix **parcsr)
1586: {
1587:   Mat_HYPRE *hA = (Mat_HYPRE*)A->data;
1588:   HYPRE_Int type;

1591:   if (!hA->ij) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"HYPRE_IJMatrix not present");
1592:   PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(hA->ij,&type));
1593:   if (type != HYPRE_PARCSR) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"HYPRE_IJMatrix is not of type HYPRE_PARCSR");
1594:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hA->ij,(void**)parcsr));
1595:   return(0);
1596: }

1598: /*
1599:    MatHYPREGetParCSR - Gets the pointer to the ParCSR matrix

1601:    Not collective

1603:    Input Parameters:
1604: +  A  - the MATHYPRE object

1606:    Output Parameter:
1607: .  parcsr  - the pointer to the hypre_ParCSRMatrix

1609:    Level: intermediate

1611: .seealso: MatHYPRE, PetscCopyMode
1612: */
1613: PetscErrorCode MatHYPREGetParCSR(Mat A, hypre_ParCSRMatrix **parcsr)
1614: {

1620:   PetscUseMethod(A,"MatHYPREGetParCSR_C",(Mat,hypre_ParCSRMatrix**),(A,parcsr));
1621:   return(0);
1622: }

1624: static PetscErrorCode MatMissingDiagonal_HYPRE(Mat A, PetscBool *missing, PetscInt *dd)
1625: {
1626:   hypre_ParCSRMatrix *parcsr;
1627:   hypre_CSRMatrix    *ha;
1628:   PetscInt           rst;
1629:   PetscErrorCode     ierr;

1632:   if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not implemented with non-square diagonal blocks");
1633:   MatGetOwnershipRange(A,&rst,NULL);
1634:   MatHYPREGetParCSR_HYPRE(A,&parcsr);
1635:   if (missing) *missing = PETSC_FALSE;
1636:   if (dd) *dd = -1;
1637:   ha = hypre_ParCSRMatrixDiag(parcsr);
1638:   if (ha) {
1639:     PetscInt  size,i;
1640:     HYPRE_Int *ii,*jj;

1642:     size = hypre_CSRMatrixNumRows(ha);
1643:     ii   = hypre_CSRMatrixI(ha);
1644:     jj   = hypre_CSRMatrixJ(ha);
1645:     for (i = 0; i < size; i++) {
1646:       PetscInt  j;
1647:       PetscBool found = PETSC_FALSE;

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

1652:       if (!found) {
1653:         PetscInfo1(A,"Matrix is missing local diagonal entry %D\n",i);
1654:         if (missing) *missing = PETSC_TRUE;
1655:         if (dd) *dd = i+rst;
1656:         return(0);
1657:       }
1658:     }
1659:     if (!size) {
1660:       PetscInfo(A,"Matrix has no diagonal entries therefore is missing diagonal\n");
1661:       if (missing) *missing = PETSC_TRUE;
1662:       if (dd) *dd = rst;
1663:     }
1664:   } else {
1665:     PetscInfo(A,"Matrix has no diagonal entries therefore is missing diagonal\n");
1666:     if (missing) *missing = PETSC_TRUE;
1667:     if (dd) *dd = rst;
1668:   }
1669:   return(0);
1670: }

1672: static PetscErrorCode MatScale_HYPRE(Mat A, PetscScalar s)
1673: {
1674:   hypre_ParCSRMatrix *parcsr;
1675:   hypre_CSRMatrix    *ha;
1676:   PetscErrorCode     ierr;
1677:   HYPRE_Complex      hs;

1680:   PetscHYPREScalarCast(s,&hs);
1681:   MatHYPREGetParCSR_HYPRE(A,&parcsr);
1682:   /* diagonal part */
1683:   ha = hypre_ParCSRMatrixDiag(parcsr);
1684:   if (ha) {
1685:     PetscInt      size,i;
1686:     HYPRE_Int     *ii;
1687:     HYPRE_Complex *a;

1689:     size = hypre_CSRMatrixNumRows(ha);
1690:     a    = hypre_CSRMatrixData(ha);
1691:     ii   = hypre_CSRMatrixI(ha);
1692:     for (i = 0; i < ii[size]; i++) a[i] *= hs;
1693:   }
1694:   /* offdiagonal part */
1695:   ha = hypre_ParCSRMatrixOffd(parcsr);
1696:   if (ha) {
1697:     PetscInt      size,i;
1698:     HYPRE_Int     *ii;
1699:     HYPRE_Complex *a;

1701:     size = hypre_CSRMatrixNumRows(ha);
1702:     a    = hypre_CSRMatrixData(ha);
1703:     ii   = hypre_CSRMatrixI(ha);
1704:     for (i = 0; i < ii[size]; i++) a[i] *= hs;
1705:   }
1706:   return(0);
1707: }

1709: static PetscErrorCode MatZeroRowsColumns_HYPRE(Mat A, PetscInt numRows, const PetscInt rows[], PetscScalar diag, Vec x, Vec b)
1710: {
1711:   hypre_ParCSRMatrix *parcsr;
1712:   HYPRE_Int          *lrows;
1713:   PetscInt           rst,ren,i;
1714:   PetscErrorCode     ierr;

1717:   if (x || b) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"To be implemented");
1718:   MatHYPREGetParCSR_HYPRE(A,&parcsr);
1719:   PetscMalloc1(numRows,&lrows);
1720:   MatGetOwnershipRange(A,&rst,&ren);
1721:   for (i=0;i<numRows;i++) {
1722:     if (rows[i] < rst || rows[i] >= ren)
1723:       SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Non-local rows not yet supported");
1724:     lrows[i] = rows[i] - rst;
1725:   }
1726:   PetscStackCallStandard(hypre_ParCSRMatrixEliminateRowsCols,(parcsr,numRows,lrows));
1727:   PetscFree(lrows);
1728:   return(0);
1729: }

1731: static PetscErrorCode MatZeroEntries_HYPRE_CSRMatrix(hypre_CSRMatrix *ha)
1732: {
1733:   PetscErrorCode      ierr;

1736:   if (ha) {
1737:     HYPRE_Int     *ii, size;
1738:     HYPRE_Complex *a;

1740:     size = hypre_CSRMatrixNumRows(ha);
1741:     a    = hypre_CSRMatrixData(ha);
1742:     ii   = hypre_CSRMatrixI(ha);

1744:     if (a) {PetscArrayzero(a,ii[size]);}
1745:   }
1746:   return(0);
1747: }

1749: PetscErrorCode MatZeroEntries_HYPRE(Mat A)
1750: {
1751:   hypre_ParCSRMatrix  *parcsr;
1752:   PetscErrorCode      ierr;

1755:   MatHYPREGetParCSR_HYPRE(A,&parcsr);
1756:   /* diagonal part */
1757:   MatZeroEntries_HYPRE_CSRMatrix(hypre_ParCSRMatrixDiag(parcsr));
1758:   /* off-diagonal part */
1759:   MatZeroEntries_HYPRE_CSRMatrix(hypre_ParCSRMatrixOffd(parcsr));
1760:   return(0);
1761: }

1763: static PetscErrorCode MatZeroRows_HYPRE_CSRMatrix(hypre_CSRMatrix *hA,PetscInt N,const PetscInt rows[],HYPRE_Complex diag)
1764: {
1765:   PetscInt        ii;
1766:   HYPRE_Int       *i, *j;
1767:   HYPRE_Complex   *a;

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

1772:   i = hypre_CSRMatrixI(hA);
1773:   j = hypre_CSRMatrixJ(hA);
1774:   a = hypre_CSRMatrixData(hA);

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

1779:     irow = rows[ii];
1780:     ibeg = i[irow];
1781:     iend = i[irow+1];
1782:     for (jj = ibeg; jj < iend; jj++)
1783:       if (j[jj] == irow) a[jj] = diag;
1784:       else a[jj] = 0.0;
1785:    }
1786:    return(0);
1787: }

1789: static PetscErrorCode MatZeroRows_HYPRE(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
1790: {
1791:   hypre_ParCSRMatrix  *parcsr;
1792:   PetscInt            *lrows,len;
1793:   HYPRE_Complex       hdiag;
1794:   PetscErrorCode      ierr;

1797:   if (x || b) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Does not support to modify the solution and the right hand size");
1798:   PetscHYPREScalarCast(diag,&hdiag);
1799:   /* retrieve the internal matrix */
1800:   MatHYPREGetParCSR_HYPRE(A,&parcsr);
1801:   /* get locally owned rows */
1802:   MatZeroRowsMapLocal_Private(A,N,rows,&len,&lrows);
1803:   /* zero diagonal part */
1804:   MatZeroRows_HYPRE_CSRMatrix(hypre_ParCSRMatrixDiag(parcsr),len,lrows,hdiag);
1805:   /* zero off-diagonal part */
1806:   MatZeroRows_HYPRE_CSRMatrix(hypre_ParCSRMatrixOffd(parcsr),len,lrows,0.0);

1808:   PetscFree(lrows);
1809:   return(0);
1810: }

1812: static PetscErrorCode MatAssemblyBegin_HYPRE(Mat mat,MatAssemblyType mode)
1813: {

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

1819:   MatStashScatterBegin_Private(mat,&mat->stash,mat->rmap->range);
1820:   return(0);
1821: }

1823: static PetscErrorCode MatGetRow_HYPRE(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v)
1824: {
1825:   hypre_ParCSRMatrix  *parcsr;
1826:   HYPRE_Int           hnz;
1827:   PetscErrorCode      ierr;

1830:   /* retrieve the internal matrix */
1831:   MatHYPREGetParCSR_HYPRE(A,&parcsr);
1832:   /* call HYPRE API */
1833:   PetscStackCallStandard(HYPRE_ParCSRMatrixGetRow,(parcsr,row,&hnz,(HYPRE_BigInt**)idx,(HYPRE_Complex**)v));
1834:   if (nz) *nz = (PetscInt)hnz;
1835:   return(0);
1836: }

1838: static PetscErrorCode MatRestoreRow_HYPRE(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v)
1839: {
1840:   hypre_ParCSRMatrix  *parcsr;
1841:   HYPRE_Int           hnz;
1842:   PetscErrorCode      ierr;

1845:   /* retrieve the internal matrix */
1846:   MatHYPREGetParCSR_HYPRE(A,&parcsr);
1847:   /* call HYPRE API */
1848:   hnz = nz ? (HYPRE_Int)(*nz) : 0;
1849:   PetscStackCallStandard(HYPRE_ParCSRMatrixRestoreRow,(parcsr,row,&hnz,(HYPRE_BigInt**)idx,(HYPRE_Complex**)v));
1850:   return(0);
1851: }

1853: static PetscErrorCode MatGetValues_HYPRE(Mat A,PetscInt m,const PetscInt idxm[],PetscInt n,const PetscInt idxn[],PetscScalar v[])
1854: {
1855:   Mat_HYPRE *hA = (Mat_HYPRE*)A->data;
1856:   PetscInt  i;

1859:   if (!m || !n) return(0);
1860:   /* Ignore negative row indices
1861:    * And negative column indices should be automatically ignored in hypre
1862:    * */
1863:   for (i=0; i<m; i++) {
1864:     if (idxm[i] >= 0) {
1865:       HYPRE_Int hn = (HYPRE_Int)n;
1866:       PetscStackCallStandard(HYPRE_IJMatrixGetValues,(hA->ij,1,&hn,(HYPRE_BigInt*)&idxm[i],(HYPRE_BigInt*)idxn,(HYPRE_Complex*)(v + i*n)));
1867:     }
1868:   }
1869:   return(0);
1870: }

1872: static PetscErrorCode MatSetOption_HYPRE(Mat A,MatOption op,PetscBool flg)
1873: {
1874:   Mat_HYPRE *hA = (Mat_HYPRE*)A->data;

1877:   switch (op) {
1878:   case MAT_NO_OFF_PROC_ENTRIES:
1879:     if (flg) {
1880:       PetscStackCallStandard(HYPRE_IJMatrixSetMaxOffProcElmts,(hA->ij,0));
1881:     }
1882:     break;
1883:   case MAT_SORTED_FULL:
1884:     hA->sorted_full = flg;
1885:     break;
1886:   default:
1887:     break;
1888:   }
1889:   return(0);
1890: }

1892: static PetscErrorCode MatView_HYPRE(Mat A, PetscViewer view)
1893: {
1894:   hypre_ParCSRMatrix *parcsr;
1895:   PetscErrorCode     ierr;
1896:   Mat                B;
1897:   PetscViewerFormat  format;
1898:   PetscErrorCode     (*mview)(Mat,PetscViewer) = NULL;

1901:   PetscViewerGetFormat(view,&format);
1902:   if (format != PETSC_VIEWER_NATIVE) {
1903:     MatHYPREGetParCSR_HYPRE(A,&parcsr);
1904:     MatCreateFromParCSR(parcsr,MATAIJ,PETSC_USE_POINTER,&B);
1905:     MatGetOperation(B,MATOP_VIEW,(void(**)(void))&mview);
1906:     if (!mview) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Missing view operation");
1907:     (*mview)(B,view);
1908:     MatDestroy(&B);
1909:   } else {
1910:     Mat_HYPRE  *hA = (Mat_HYPRE*)A->data;
1911:     PetscMPIInt size;
1912:     PetscBool   isascii;
1913:     const char *filename;

1915:     /* HYPRE uses only text files */
1916:     PetscObjectTypeCompare((PetscObject)view,PETSCVIEWERASCII,&isascii);
1917:     if (!isascii) SETERRQ1(PetscObjectComm((PetscObject)view),PETSC_ERR_SUP,"PetscViewerType %s: native HYPRE format needs PETSCVIEWERASCII",((PetscObject)view)->type_name);
1918:     PetscViewerFileGetName(view,&filename);
1919:     PetscStackCallStandard(HYPRE_IJMatrixPrint,(hA->ij,filename));
1920:     MPI_Comm_size(hA->comm,&size);
1921:     if (size > 1) {
1922:       PetscViewerASCIIPrintf(view,"Matrix files: %s.%05d ... %s.%05d\n",filename,0,filename,size-1);
1923:     } else {
1924:       PetscViewerASCIIPrintf(view,"Matrix file: %s.%05d\n",filename,0);
1925:     }
1926:   }
1927:   return(0);
1928: }

1930: static PetscErrorCode MatDuplicate_HYPRE(Mat A,MatDuplicateOption op, Mat *B)
1931: {
1932:   hypre_ParCSRMatrix *parcsr;
1933:   PetscErrorCode     ierr;
1934:   PetscCopyMode      cpmode;

1937:   MatHYPREGetParCSR_HYPRE(A,&parcsr);
1938:   if (op == MAT_DO_NOT_COPY_VALUES || op == MAT_SHARE_NONZERO_PATTERN) {
1939:     parcsr = hypre_ParCSRMatrixClone(parcsr,0);
1940:     cpmode = PETSC_OWN_POINTER;
1941:   } else {
1942:     cpmode = PETSC_COPY_VALUES;
1943:   }
1944:   MatCreateFromParCSR(parcsr,MATHYPRE,cpmode,B);
1945:   return(0);
1946: }

1948: static PetscErrorCode MatCopy_HYPRE(Mat A, Mat B, MatStructure str)
1949: {
1950:   hypre_ParCSRMatrix *acsr,*bcsr;
1951:   PetscErrorCode     ierr;

1954:   if (str == SAME_NONZERO_PATTERN && A->ops->copy == B->ops->copy) {
1955:     MatHYPREGetParCSR_HYPRE(A,&acsr);
1956:     MatHYPREGetParCSR_HYPRE(B,&bcsr);
1957:     PetscStackCallStandard(hypre_ParCSRMatrixCopy,(acsr,bcsr,1));
1958:     MatSetOption(B,MAT_SORTED_FULL,PETSC_TRUE); /* "perfect" preallocation, so no need for hypre_AuxParCSRMatrixNeedAux */
1959:     MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
1960:     MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
1961:   } else {
1962:     MatCopy_Basic(A,B,str);
1963:   }
1964:   return(0);
1965: }

1967: static PetscErrorCode MatGetDiagonal_HYPRE(Mat A, Vec d)
1968: {
1969:   hypre_ParCSRMatrix *parcsr;
1970:   hypre_CSRMatrix    *dmat;
1971:   HYPRE_Complex      *a;
1972:   HYPRE_Complex      *data = NULL;
1973:   HYPRE_Int          *diag = NULL;
1974:   PetscInt           i;
1975:   PetscBool          cong;
1976:   PetscErrorCode     ierr;

1979:   MatHasCongruentLayouts(A,&cong);
1980:   if (!cong) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Only for square matrices with same local distributions of rows and columns");
1981:   if (PetscDefined(USE_DEBUG)) {
1982:     PetscBool miss;
1983:     MatMissingDiagonal(A,&miss,NULL);
1984:     if (miss && A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not implemented when diagonal entries are missing");
1985:   }
1986:   MatHYPREGetParCSR_HYPRE(A,&parcsr);
1987:   dmat = hypre_ParCSRMatrixDiag(parcsr);
1988:   if (dmat) {
1989:     /* this cast fixes the clang error: implicit conversion from 'HYPRE_Complex' (aka '_Complex double') to 'double' is not permitted in C++ */
1990:     VecGetArray(d,(PetscScalar**)&a);
1991:     diag = hypre_CSRMatrixI(dmat);
1992:     data = hypre_CSRMatrixData(dmat);
1993:     for (i=0;i<A->rmap->n;i++) a[i] = data[diag[i]];
1994:     VecRestoreArray(d,(PetscScalar**)&a);
1995:   }
1996:   return(0);
1997: }

1999: #include <petscblaslapack.h>

2001: static PetscErrorCode MatAXPY_HYPRE(Mat Y,PetscScalar a,Mat X,MatStructure str)
2002: {

2006:   if (str == SAME_NONZERO_PATTERN) {
2007:     hypre_ParCSRMatrix *x,*y;
2008:     hypre_CSRMatrix    *xloc,*yloc;
2009:     PetscInt           xnnz,ynnz;
2010:     HYPRE_Complex      *xarr,*yarr;
2011:     PetscBLASInt       one=1,bnz;

2013:     MatHYPREGetParCSR(Y,&y);
2014:     MatHYPREGetParCSR(X,&x);

2016:     /* diagonal block */
2017:     xloc = hypre_ParCSRMatrixDiag(x);
2018:     yloc = hypre_ParCSRMatrixDiag(y);
2019:     xnnz = 0;
2020:     ynnz = 0;
2021:     xarr = NULL;
2022:     yarr = NULL;
2023:     if (xloc) {
2024:       xarr = hypre_CSRMatrixData(xloc);
2025:       xnnz = hypre_CSRMatrixNumNonzeros(xloc);
2026:     }
2027:     if (yloc) {
2028:       yarr = hypre_CSRMatrixData(yloc);
2029:       ynnz = hypre_CSRMatrixNumNonzeros(yloc);
2030:     }
2031:     if (xnnz != ynnz) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Different number of nonzeros in diagonal block %D != %D",xnnz,ynnz);
2032:     PetscBLASIntCast(xnnz,&bnz);
2033:     PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&bnz,&a,(PetscScalar*)xarr,&one,(PetscScalar*)yarr,&one));

2035:     /* off-diagonal block */
2036:     xloc = hypre_ParCSRMatrixOffd(x);
2037:     yloc = hypre_ParCSRMatrixOffd(y);
2038:     xnnz = 0;
2039:     ynnz = 0;
2040:     xarr = NULL;
2041:     yarr = NULL;
2042:     if (xloc) {
2043:       xarr = hypre_CSRMatrixData(xloc);
2044:       xnnz = hypre_CSRMatrixNumNonzeros(xloc);
2045:     }
2046:     if (yloc) {
2047:       yarr = hypre_CSRMatrixData(yloc);
2048:       ynnz = hypre_CSRMatrixNumNonzeros(yloc);
2049:     }
2050:     if (xnnz != ynnz) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Different number of nonzeros in off-diagonal block %D != %D",xnnz,ynnz);
2051:     PetscBLASIntCast(xnnz,&bnz);
2052:     PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&bnz,&a,(PetscScalar*)xarr,&one,(PetscScalar*)yarr,&one));
2053:   } else if (str == SUBSET_NONZERO_PATTERN) {
2054:     MatAXPY_Basic(Y,a,X,str);
2055:   } else {
2056:     Mat B;

2058:     MatAXPY_Basic_Preallocate(Y,X,&B);
2059:     MatAXPY_BasicWithPreallocation(B,Y,a,X,str);
2060:     MatHeaderReplace(Y,&B);
2061:   }
2062:   return(0);
2063: }

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

2069:    Level: intermediate

2071: .seealso: MatCreate()
2072: M*/

2074: PETSC_EXTERN PetscErrorCode MatCreate_HYPRE(Mat B)
2075: {
2076:   Mat_HYPRE      *hB;

2080:   PetscNewLog(B,&hB);
2081:   hB->inner_free = PETSC_TRUE;
2082:   hB->available  = PETSC_TRUE;
2083:   hB->sorted_full= PETSC_FALSE; /* no assumption whether column indices are sorted or not */
2084:   hB->size       = 0;
2085:   hB->array      = NULL;

2087:   B->data       = (void*)hB;
2088:   B->rmap->bs   = 1;
2089:   B->assembled  = PETSC_FALSE;

2091:   PetscMemzero(B->ops,sizeof(struct _MatOps));
2092:   B->ops->mult             = MatMult_HYPRE;
2093:   B->ops->multtranspose    = MatMultTranspose_HYPRE;
2094:   B->ops->multadd          = MatMultAdd_HYPRE;
2095:   B->ops->multtransposeadd = MatMultTransposeAdd_HYPRE;
2096:   B->ops->setup            = MatSetUp_HYPRE;
2097:   B->ops->destroy          = MatDestroy_HYPRE;
2098:   B->ops->assemblyend      = MatAssemblyEnd_HYPRE;
2099:   B->ops->assemblybegin    = MatAssemblyBegin_HYPRE;
2100:   B->ops->setvalues        = MatSetValues_HYPRE;
2101:   B->ops->missingdiagonal  = MatMissingDiagonal_HYPRE;
2102:   B->ops->scale            = MatScale_HYPRE;
2103:   B->ops->zerorowscolumns  = MatZeroRowsColumns_HYPRE;
2104:   B->ops->zeroentries      = MatZeroEntries_HYPRE;
2105:   B->ops->zerorows         = MatZeroRows_HYPRE;
2106:   B->ops->getrow           = MatGetRow_HYPRE;
2107:   B->ops->restorerow       = MatRestoreRow_HYPRE;
2108:   B->ops->getvalues        = MatGetValues_HYPRE;
2109:   B->ops->setoption        = MatSetOption_HYPRE;
2110:   B->ops->duplicate        = MatDuplicate_HYPRE;
2111:   B->ops->copy             = MatCopy_HYPRE;
2112:   B->ops->view             = MatView_HYPRE;
2113:   B->ops->getdiagonal      = MatGetDiagonal_HYPRE;
2114:   B->ops->axpy             = MatAXPY_HYPRE;
2115:   B->ops->productsetfromoptions = MatProductSetFromOptions_HYPRE;

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

2120:   MPI_Comm_dup(PetscObjectComm((PetscObject)B),&hB->comm);
2121:   PetscObjectChangeTypeName((PetscObject)B,MATHYPRE);
2122:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_hypre_aij_C",MatConvert_HYPRE_AIJ);
2123:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_hypre_is_C",MatConvert_HYPRE_IS);
2124:   PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqaij_hypre_C",MatProductSetFromOptions_HYPRE);
2125:   PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_mpiaij_hypre_C",MatProductSetFromOptions_HYPRE);
2126:   PetscObjectComposeFunction((PetscObject)B,"MatHYPRESetPreallocation_C",MatHYPRESetPreallocation_HYPRE);
2127:   PetscObjectComposeFunction((PetscObject)B,"MatHYPREGetParCSR_C",MatHYPREGetParCSR_HYPRE);
2128:   return(0);
2129: }

2131: static PetscErrorCode hypre_array_destroy(void *ptr)
2132: {
2134:    hypre_TFree(ptr,HYPRE_MEMORY_HOST);
2135:    return(0);
2136: }