Actual source code: mhypre.c

petsc-3.11.3 2019-06-26
Report Typos and Errors

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

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

 16: PETSC_INTERN PetscErrorCode MatPtAP_IS_XAIJ(Mat,Mat,MatReuse,PetscReal,Mat*);

 18: static PetscErrorCode MatHYPRE_CreateFromMat(Mat,Mat_HYPRE*);
 19: static PetscErrorCode MatHYPRE_IJMatrixPreallocate(Mat,Mat,HYPRE_IJMatrix);
 20: static PetscErrorCode MatHYPRE_IJMatrixFastCopy_MPIAIJ(Mat,HYPRE_IJMatrix);
 21: static PetscErrorCode MatHYPRE_IJMatrixFastCopy_SeqAIJ(Mat,HYPRE_IJMatrix);
 22: static PetscErrorCode MatHYPRE_MultKernel_Private(Mat,PetscScalar,Vec,PetscScalar,Vec,PetscBool);
 23: static PetscErrorCode hypre_array_destroy(void*);
 24: PetscErrorCode MatSetValues_HYPRE(Mat, PetscInt,const PetscInt[],PetscInt,const PetscInt[],const PetscScalar[],InsertMode ins);

 26: static PetscErrorCode MatHYPRE_IJMatrixPreallocate(Mat A_d, Mat A_o, HYPRE_IJMatrix ij)
 27: {
 29:   PetscInt       i,n_d,n_o;
 30:   const PetscInt *ia_d,*ia_o;
 31:   PetscBool      done_d=PETSC_FALSE,done_o=PETSC_FALSE;
 32:   PetscInt       *nnz_d=NULL,*nnz_o=NULL;

 35:   if (A_d) { /* determine number of nonzero entries in local diagonal part */
 36:     MatGetRowIJ(A_d,0,PETSC_FALSE,PETSC_FALSE,&n_d,&ia_d,NULL,&done_d);
 37:     if (done_d) {
 38:       PetscMalloc1(n_d,&nnz_d);
 39:       for (i=0; i<n_d; i++) {
 40:         nnz_d[i] = ia_d[i+1] - ia_d[i];
 41:       }
 42:     }
 43:     MatRestoreRowIJ(A_d,0,PETSC_FALSE,PETSC_FALSE,NULL,&ia_d,NULL,&done_d);
 44:   }
 45:   if (A_o) { /* determine number of nonzero entries in local off-diagonal part */
 46:     MatGetRowIJ(A_o,0,PETSC_FALSE,PETSC_FALSE,&n_o,&ia_o,NULL,&done_o);
 47:     if (done_o) {
 48:       PetscMalloc1(n_o,&nnz_o);
 49:       for (i=0; i<n_o; i++) {
 50:         nnz_o[i] = ia_o[i+1] - ia_o[i];
 51:       }
 52:     }
 53:     MatRestoreRowIJ(A_o,0,PETSC_FALSE,PETSC_FALSE,&n_o,&ia_o,NULL,&done_o);
 54:   }
 55:   if (done_d) {    /* set number of nonzeros in HYPRE IJ matrix */
 56:     if (!done_o) { /* only diagonal part */
 57:       PetscMalloc1(n_d,&nnz_o);
 58:       for (i=0; i<n_d; i++) {
 59:         nnz_o[i] = 0;
 60:       }
 61:     }
 62:     PetscStackCallStandard(HYPRE_IJMatrixSetDiagOffdSizes,(ij,(HYPRE_Int *)nnz_d,(HYPRE_Int *)nnz_o));
 63:     PetscFree(nnz_d);
 64:     PetscFree(nnz_o);
 65:   }
 66:   return(0);
 67: }

 69: static PetscErrorCode MatHYPRE_CreateFromMat(Mat A, Mat_HYPRE *hA)
 70: {
 72:   PetscInt       rstart,rend,cstart,cend;

 75:   PetscLayoutSetUp(A->rmap);
 76:   PetscLayoutSetUp(A->cmap);
 77:   rstart = A->rmap->rstart;
 78:   rend   = A->rmap->rend;
 79:   cstart = A->cmap->rstart;
 80:   cend   = A->cmap->rend;
 81:   PetscStackCallStandard(HYPRE_IJMatrixCreate,(hA->comm,rstart,rend-1,cstart,cend-1,&hA->ij));
 82:   PetscStackCallStandard(HYPRE_IJMatrixSetObjectType,(hA->ij,HYPRE_PARCSR));
 83:   {
 84:     PetscBool      same;
 85:     Mat            A_d,A_o;
 86:     const PetscInt *colmap;
 87:     PetscObjectTypeCompare((PetscObject)A,MATMPIAIJ,&same);
 88:     if (same) {
 89:       MatMPIAIJGetSeqAIJ(A,&A_d,&A_o,&colmap);
 90:       MatHYPRE_IJMatrixPreallocate(A_d,A_o,hA->ij);
 91:       return(0);
 92:     }
 93:     PetscObjectTypeCompare((PetscObject)A,MATMPIBAIJ,&same);
 94:     if (same) {
 95:       MatMPIBAIJGetSeqBAIJ(A,&A_d,&A_o,&colmap);
 96:       MatHYPRE_IJMatrixPreallocate(A_d,A_o,hA->ij);
 97:       return(0);
 98:     }
 99:     PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&same);
100:     if (same) {
101:       MatHYPRE_IJMatrixPreallocate(A,NULL,hA->ij);
102:       return(0);
103:     }
104:     PetscObjectTypeCompare((PetscObject)A,MATSEQBAIJ,&same);
105:     if (same) {
106:       MatHYPRE_IJMatrixPreallocate(A,NULL,hA->ij);
107:       return(0);
108:     }
109:   }
110:   return(0);
111: }

113: static PetscErrorCode MatHYPRE_IJMatrixCopy(Mat A, HYPRE_IJMatrix ij)
114: {
115:   PetscErrorCode    ierr;
116:   PetscInt          i,rstart,rend,ncols,nr,nc;
117:   const PetscScalar *values;
118:   const PetscInt    *cols;
119:   PetscBool         flg;

122:   PetscObjectTypeCompare((PetscObject)A,MATMPIAIJ,&flg);
123:   MatGetSize(A,&nr,&nc);
124:   if (flg && nr == nc) {
125:     MatHYPRE_IJMatrixFastCopy_MPIAIJ(A,ij);
126:     return(0);
127:   }
128:   PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&flg);
129:   if (flg) {
130:     MatHYPRE_IJMatrixFastCopy_SeqAIJ(A,ij);
131:     return(0);
132:   }

134:   PetscStackCallStandard(HYPRE_IJMatrixInitialize,(ij));
135:   MatGetOwnershipRange(A,&rstart,&rend);
136:   for (i=rstart; i<rend; i++) {
137:     MatGetRow(A,i,&ncols,&cols,&values);
138:     if (ncols) {
139:       PetscStackCallStandard(HYPRE_IJMatrixSetValues,(ij,1,(HYPRE_Int *)&ncols,(HYPRE_Int *)&i,(HYPRE_Int *)cols,values));
140:     }
141:     MatRestoreRow(A,i,&ncols,&cols,&values);
142:   }
143:   return(0);
144: }

146: static PetscErrorCode MatHYPRE_IJMatrixFastCopy_SeqAIJ(Mat A, HYPRE_IJMatrix ij)
147: {
148:   PetscErrorCode        ierr;
149:   Mat_SeqAIJ            *pdiag = (Mat_SeqAIJ*)A->data;
150:   HYPRE_Int             type;
151:   hypre_ParCSRMatrix    *par_matrix;
152:   hypre_AuxParCSRMatrix *aux_matrix;
153:   hypre_CSRMatrix       *hdiag;

156:   PetscStackCallStandard(HYPRE_IJMatrixInitialize,(ij));
157:   PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(ij,&type));
158:   if (type != HYPRE_PARCSR) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported");
159:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(ij,(void**)&par_matrix));
160:   hdiag = hypre_ParCSRMatrixDiag(par_matrix);
161:   /*
162:        this is the Hack part where we monkey directly with the hypre datastructures
163:   */
164:   PetscMemcpy(hdiag->i,pdiag->i,(A->rmap->n + 1)*sizeof(PetscInt));
165:   PetscMemcpy(hdiag->j,pdiag->j,pdiag->nz*sizeof(PetscInt));
166:   PetscMemcpy(hdiag->data,pdiag->a,pdiag->nz*sizeof(PetscScalar));

168:   aux_matrix = (hypre_AuxParCSRMatrix*)hypre_IJMatrixTranslator(ij);
169:   hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 0;
170:   return(0);
171: }

173: static PetscErrorCode MatHYPRE_IJMatrixFastCopy_MPIAIJ(Mat A, HYPRE_IJMatrix ij)
174: {
175:   PetscErrorCode        ierr;
176:   Mat_MPIAIJ            *pA = (Mat_MPIAIJ*)A->data;
177:   Mat_SeqAIJ            *pdiag,*poffd;
178:   PetscInt              i,*garray = pA->garray,*jj,cstart,*pjj;
179:   HYPRE_Int             type;
180:   hypre_ParCSRMatrix    *par_matrix;
181:   hypre_AuxParCSRMatrix *aux_matrix;
182:   hypre_CSRMatrix       *hdiag,*hoffd;

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

190:   PetscStackCallStandard(HYPRE_IJMatrixInitialize,(ij));
191:   PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(ij,&type));
192:   if (type != HYPRE_PARCSR) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported");
193:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(ij,(void**)&par_matrix));
194:   hdiag = hypre_ParCSRMatrixDiag(par_matrix);
195:   hoffd = hypre_ParCSRMatrixOffd(par_matrix);

197:   /*
198:        this is the Hack part where we monkey directly with the hypre datastructures
199:   */
200:   PetscMemcpy(hdiag->i,pdiag->i,(pA->A->rmap->n + 1)*sizeof(PetscInt));
201:   /* need to shift the diag column indices (hdiag->j) back to global numbering since hypre is expecting this */
202:   jj  = (PetscInt*)hdiag->j;
203:   pjj = (PetscInt*)pdiag->j;
204:   for (i=0; i<pdiag->nz; i++) jj[i] = cstart + pjj[i];
205:   PetscMemcpy(hdiag->data,pdiag->a,pdiag->nz*sizeof(PetscScalar));
206:   PetscMemcpy(hoffd->i,poffd->i,(pA->A->rmap->n + 1)*sizeof(PetscInt));
207:   /* need to move the offd column indices (hoffd->j) back to global numbering since hypre is expecting this
208:      If we hacked a hypre a bit more we might be able to avoid this step */
209:   jj  = (PetscInt*) hoffd->j;
210:   pjj = (PetscInt*) poffd->j;
211:   for (i=0; i<poffd->nz; i++) jj[i] = garray[pjj[i]];
212:   PetscMemcpy(hoffd->data,poffd->a,poffd->nz*sizeof(PetscScalar));

214:   aux_matrix = (hypre_AuxParCSRMatrix*)hypre_IJMatrixTranslator(ij);
215:   hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 0;
216:   return(0);
217: }

219: static PetscErrorCode MatConvert_HYPRE_IS(Mat A, MatType mtype, MatReuse reuse, Mat* B)
220: {
221:   Mat_HYPRE*             mhA = (Mat_HYPRE*)(A->data);
222:   Mat                    lA;
223:   ISLocalToGlobalMapping rl2g,cl2g;
224:   IS                     is;
225:   hypre_ParCSRMatrix     *hA;
226:   hypre_CSRMatrix        *hdiag,*hoffd;
227:   MPI_Comm               comm;
228:   PetscScalar            *hdd,*hod,*aa,*data;
229:   HYPRE_Int              *col_map_offd,*hdi,*hdj,*hoi,*hoj;
230:   PetscInt               *ii,*jj,*iptr,*jptr;
231:   PetscInt               cum,dr,dc,oc,str,stc,nnz,i,jd,jo,M,N;
232:   HYPRE_Int              type;
233:   PetscErrorCode         ierr;

236:   comm = PetscObjectComm((PetscObject)A);
237:   PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(mhA->ij,&type));
238:   if (type != HYPRE_PARCSR) SETERRQ(comm,PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported");
239:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(mhA->ij,(void**)&hA));
240:   M     = hypre_ParCSRMatrixGlobalNumRows(hA);
241:   N     = hypre_ParCSRMatrixGlobalNumCols(hA);
242:   str   = hypre_ParCSRMatrixFirstRowIndex(hA);
243:   stc   = hypre_ParCSRMatrixFirstColDiag(hA);
244:   hdiag = hypre_ParCSRMatrixDiag(hA);
245:   hoffd = hypre_ParCSRMatrixOffd(hA);
246:   dr    = hypre_CSRMatrixNumRows(hdiag);
247:   dc    = hypre_CSRMatrixNumCols(hdiag);
248:   nnz   = hypre_CSRMatrixNumNonzeros(hdiag);
249:   hdi   = hypre_CSRMatrixI(hdiag);
250:   hdj   = hypre_CSRMatrixJ(hdiag);
251:   hdd   = hypre_CSRMatrixData(hdiag);
252:   oc    = hypre_CSRMatrixNumCols(hoffd);
253:   nnz  += hypre_CSRMatrixNumNonzeros(hoffd);
254:   hoi   = hypre_CSRMatrixI(hoffd);
255:   hoj   = hypre_CSRMatrixJ(hoffd);
256:   hod   = hypre_CSRMatrixData(hoffd);
257:   if (reuse != MAT_REUSE_MATRIX) {
258:     PetscInt *aux;

260:     /* generate l2g maps for rows and cols */
261:     ISCreateStride(comm,dr,str,1,&is);
262:     ISLocalToGlobalMappingCreateIS(is,&rl2g);
263:     ISDestroy(&is);
264:     col_map_offd = hypre_ParCSRMatrixColMapOffd(hA);
265:     PetscMalloc1(dc+oc,&aux);
266:     for (i=0; i<dc; i++) aux[i] = i+stc;
267:     for (i=0; i<oc; i++) aux[i+dc] = col_map_offd[i];
268:     ISCreateGeneral(comm,dc+oc,aux,PETSC_OWN_POINTER,&is);
269:     ISLocalToGlobalMappingCreateIS(is,&cl2g);
270:     ISDestroy(&is);
271:     /* create MATIS object */
272:     MatCreate(comm,B);
273:     MatSetSizes(*B,dr,dc,M,N);
274:     MatSetType(*B,MATIS);
275:     MatSetLocalToGlobalMapping(*B,rl2g,cl2g);
276:     ISLocalToGlobalMappingDestroy(&rl2g);
277:     ISLocalToGlobalMappingDestroy(&cl2g);

279:     /* allocate CSR for local matrix */
280:     PetscMalloc1(dr+1,&iptr);
281:     PetscMalloc1(nnz,&jptr);
282:     PetscMalloc1(nnz,&data);
283:   } else {
284:     PetscInt  nr;
285:     PetscBool done;
286:     MatISGetLocalMat(*B,&lA);
287:     MatGetRowIJ(lA,0,PETSC_FALSE,PETSC_FALSE,&nr,(const PetscInt**)&iptr,(const PetscInt**)&jptr,&done);
288:     if (nr != dr) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot reuse mat: invalid number of rows in local mat! %D != %D",nr,dr);
289:     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);
290:     MatSeqAIJGetArray(lA,&data);
291:   }
292:   /* merge local matrices */
293:   ii   = iptr;
294:   jj   = jptr;
295:   aa   = data;
296:   *ii  = *(hdi++) + *(hoi++);
297:   for (jd=0,jo=0,cum=0; *ii<nnz; cum++) {
298:     PetscScalar *aold = aa;
299:     PetscInt    *jold = jj,nc = jd+jo;
300:     for (; jd<*hdi; jd++) { *jj++ = *hdj++;      *aa++ = *hdd++; }
301:     for (; jo<*hoi; jo++) { *jj++ = *hoj++ + dc; *aa++ = *hod++; }
302:     *(++ii) = *(hdi++) + *(hoi++);
303:     PetscSortIntWithScalarArray(jd+jo-nc,jold,aold);
304:   }
305:   for (; cum<dr; cum++) *(++ii) = nnz;
306:   if (reuse != MAT_REUSE_MATRIX) {
307:     Mat_SeqAIJ* a;

309:     MatCreateSeqAIJWithArrays(PETSC_COMM_SELF,dr,dc+oc,iptr,jptr,data,&lA);
310:     MatISSetLocalMat(*B,lA);
311:     /* hack SeqAIJ */
312:     a          = (Mat_SeqAIJ*)(lA->data);
313:     a->free_a  = PETSC_TRUE;
314:     a->free_ij = PETSC_TRUE;
315:     MatDestroy(&lA);
316:   }
317:   MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
318:   MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
319:   if (reuse == MAT_INPLACE_MATRIX) {
320:     MatHeaderReplace(A,B);
321:   }
322:   return(0);
323: }

325: PETSC_INTERN PetscErrorCode MatConvert_AIJ_HYPRE(Mat A, MatType type, MatReuse reuse, Mat *B)
326: {
327:   Mat            M = NULL;
328:   Mat_HYPRE      *hB;
329:   MPI_Comm       comm = PetscObjectComm((PetscObject)A);

333:   if (reuse == MAT_REUSE_MATRIX) {
334:     /* always destroy the old matrix and create a new memory;
335:        hope this does not churn the memory too much. The problem
336:        is I do not know if it is possible to put the matrix back to
337:        its initial state so that we can directly copy the values
338:        the second time through. */
339:     hB = (Mat_HYPRE*)((*B)->data);
340:     PetscStackCallStandard(HYPRE_IJMatrixDestroy,(hB->ij));
341:   } else {
342:     MatCreate(comm,&M);
343:     MatSetType(M,MATHYPRE);
344:     MatSetSizes(M,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
345:     hB   = (Mat_HYPRE*)(M->data);
346:     if (reuse == MAT_INITIAL_MATRIX) *B = M;
347:   }
348:   MatHYPRE_CreateFromMat(A,hB);
349:   MatHYPRE_IJMatrixCopy(A,hB->ij);
350:   if (reuse == MAT_INPLACE_MATRIX) {
351:     MatHeaderReplace(A,&M);
352:   }
353:   (*B)->preallocated = PETSC_TRUE;
354:   MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
355:   MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
356:   return(0);
357: }

359: static PetscErrorCode MatConvert_HYPRE_AIJ(Mat A, MatType mtype, MatReuse reuse, Mat *B)
360: {
361:   Mat_HYPRE          *hA = (Mat_HYPRE*)A->data;
362:   hypre_ParCSRMatrix *parcsr;
363:   hypre_CSRMatrix    *hdiag,*hoffd;
364:   MPI_Comm           comm;
365:   PetscScalar        *da,*oa,*aptr;
366:   PetscInt           *dii,*djj,*oii,*ojj,*iptr;
367:   PetscInt           i,dnnz,onnz,m,n;
368:   HYPRE_Int          type;
369:   PetscMPIInt        size;
370:   PetscErrorCode     ierr;

373:   comm = PetscObjectComm((PetscObject)A);
374:   PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(hA->ij,&type));
375:   if (type != HYPRE_PARCSR) SETERRQ(comm,PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported");
376:   if (reuse == MAT_REUSE_MATRIX) {
377:     PetscBool ismpiaij,isseqaij;
378:     PetscObjectTypeCompare((PetscObject)*B,MATMPIAIJ,&ismpiaij);
379:     PetscObjectBaseTypeCompare((PetscObject)*B,MATSEQAIJ,&isseqaij);
380:     if (!ismpiaij && !isseqaij) SETERRQ(comm,PETSC_ERR_SUP,"Only MATMPIAIJ or MATSEQAIJ are supported");
381:   }
382:   MPI_Comm_size(comm,&size);

384:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hA->ij,(void**)&parcsr));
385:   hdiag = hypre_ParCSRMatrixDiag(parcsr);
386:   hoffd = hypre_ParCSRMatrixOffd(parcsr);
387:   m     = hypre_CSRMatrixNumRows(hdiag);
388:   n     = hypre_CSRMatrixNumCols(hdiag);
389:   dnnz  = hypre_CSRMatrixNumNonzeros(hdiag);
390:   onnz  = hypre_CSRMatrixNumNonzeros(hoffd);
391:   if (reuse == MAT_INITIAL_MATRIX) {
392:     PetscMalloc1(m+1,&dii);
393:     PetscMalloc1(dnnz,&djj);
394:     PetscMalloc1(dnnz,&da);
395:   } else if (reuse == MAT_REUSE_MATRIX) {
396:     PetscInt  nr;
397:     PetscBool done;
398:     if (size > 1) {
399:       Mat_MPIAIJ *b = (Mat_MPIAIJ*)((*B)->data);

401:       MatGetRowIJ(b->A,0,PETSC_FALSE,PETSC_FALSE,&nr,(const PetscInt**)&dii,(const PetscInt**)&djj,&done);
402:       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);
403:       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);
404:       MatSeqAIJGetArray(b->A,&da);
405:     } else {
406:       MatGetRowIJ(*B,0,PETSC_FALSE,PETSC_FALSE,&nr,(const PetscInt**)&dii,(const PetscInt**)&djj,&done);
407:       if (nr != m) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_USER,"Cannot reuse mat: invalid number of local rows! %D != %D",nr,m);
408:       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);
409:       MatSeqAIJGetArray(*B,&da);
410:     }
411:   } else { /* MAT_INPLACE_MATRIX */
412:     dii = (PetscInt*)hypre_CSRMatrixI(hdiag);
413:     djj = (PetscInt*)hypre_CSRMatrixJ(hdiag);
414:     da  = hypre_CSRMatrixData(hdiag);
415:   }
416:   PetscMemcpy(dii,hypre_CSRMatrixI(hdiag),(m+1)*sizeof(PetscInt));
417:   PetscMemcpy(djj,hypre_CSRMatrixJ(hdiag),dnnz*sizeof(PetscInt));
418:   PetscMemcpy(da,hypre_CSRMatrixData(hdiag),dnnz*sizeof(PetscScalar));
419:   iptr = djj;
420:   aptr = da;
421:   for (i=0; i<m; i++) {
422:     PetscInt nc = dii[i+1]-dii[i];
423:     PetscSortIntWithScalarArray(nc,iptr,aptr);
424:     iptr += nc;
425:     aptr += nc;
426:   }
427:   if (size > 1) {
428:     HYPRE_Int *offdj,*coffd;

430:     if (reuse == MAT_INITIAL_MATRIX) {
431:       PetscMalloc1(m+1,&oii);
432:       PetscMalloc1(onnz,&ojj);
433:       PetscMalloc1(onnz,&oa);
434:     } else if (reuse == MAT_REUSE_MATRIX) {
435:       Mat_MPIAIJ *b = (Mat_MPIAIJ*)((*B)->data);
436:       PetscInt   nr,hr = hypre_CSRMatrixNumRows(hoffd);
437:       PetscBool  done;

439:       MatGetRowIJ(b->B,0,PETSC_FALSE,PETSC_FALSE,&nr,(const PetscInt**)&oii,(const PetscInt**)&ojj,&done);
440:       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);
441:       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);
442:       MatSeqAIJGetArray(b->B,&oa);
443:     } else { /* MAT_INPLACE_MATRIX */
444:       oii = (PetscInt*)hypre_CSRMatrixI(hoffd);
445:       ojj = (PetscInt*)hypre_CSRMatrixJ(hoffd);
446:       oa  = hypre_CSRMatrixData(hoffd);
447:     }
448:     PetscMemcpy(oii,hypre_CSRMatrixI(hoffd),(m+1)*sizeof(PetscInt));
449:     offdj = hypre_CSRMatrixJ(hoffd);
450:     coffd = hypre_ParCSRMatrixColMapOffd(parcsr);
451:     for (i=0; i<onnz; i++) ojj[i] = coffd[offdj[i]];
452:     PetscMemcpy(oa,hypre_CSRMatrixData(hoffd),onnz*sizeof(PetscScalar));
453:     iptr = ojj;
454:     aptr = oa;
455:     for (i=0; i<m; i++) {
456:        PetscInt nc = oii[i+1]-oii[i];
457:        PetscSortIntWithScalarArray(nc,iptr,aptr);
458:        iptr += nc;
459:        aptr += nc;
460:     }
461:     if (reuse == MAT_INITIAL_MATRIX) {
462:       Mat_MPIAIJ *b;
463:       Mat_SeqAIJ *d,*o;

465:       MatCreateMPIAIJWithSplitArrays(comm,m,n,PETSC_DECIDE,PETSC_DECIDE,dii,djj,da,oii,ojj,oa,B);
466:       /* hack MPIAIJ */
467:       b          = (Mat_MPIAIJ*)((*B)->data);
468:       d          = (Mat_SeqAIJ*)b->A->data;
469:       o          = (Mat_SeqAIJ*)b->B->data;
470:       d->free_a  = PETSC_TRUE;
471:       d->free_ij = PETSC_TRUE;
472:       o->free_a  = PETSC_TRUE;
473:       o->free_ij = PETSC_TRUE;
474:     } else if (reuse == MAT_INPLACE_MATRIX) {
475:       Mat T;
476:       MatCreateMPIAIJWithSplitArrays(comm,m,n,PETSC_DECIDE,PETSC_DECIDE,dii,djj,da,oii,ojj,oa,&T);
477:       hypre_CSRMatrixI(hdiag)    = NULL;
478:       hypre_CSRMatrixJ(hdiag)    = NULL;
479:       hypre_CSRMatrixData(hdiag) = NULL;
480:       hypre_CSRMatrixI(hoffd)    = NULL;
481:       hypre_CSRMatrixJ(hoffd)    = NULL;
482:       hypre_CSRMatrixData(hoffd) = NULL;
483:       MatHeaderReplace(A,&T);
484:     }
485:   } else {
486:     oii  = NULL;
487:     ojj  = NULL;
488:     oa   = NULL;
489:     if (reuse == MAT_INITIAL_MATRIX) {
490:       Mat_SeqAIJ* b;
491:       MatCreateSeqAIJWithArrays(comm,m,n,dii,djj,da,B);
492:       /* hack SeqAIJ */
493:       b          = (Mat_SeqAIJ*)((*B)->data);
494:       b->free_a  = PETSC_TRUE;
495:       b->free_ij = PETSC_TRUE;
496:     } else if (reuse == MAT_INPLACE_MATRIX) {
497:       Mat T;
498:       MatCreateSeqAIJWithArrays(comm,m,n,dii,djj,da,&T);
499:       hypre_CSRMatrixI(hdiag)    = NULL;
500:       hypre_CSRMatrixJ(hdiag)    = NULL;
501:       hypre_CSRMatrixData(hdiag) = NULL;
502:       MatHeaderReplace(A,&T);
503:     }
504:   }

506:   /* we have to use hypre_Tfree to free the arrays */
507:   if (reuse == MAT_INPLACE_MATRIX) {
508:     void *ptrs[6] = {dii,djj,da,oii,ojj,oa};
509:     const char *names[6] = {"_hypre_csr_dii",
510:                             "_hypre_csr_djj",
511:                             "_hypre_csr_da",
512:                             "_hypre_csr_oii",
513:                             "_hypre_csr_ojj",
514:                             "_hypre_csr_oa"};
515:     for (i=0; i<6; i++) {
516:       PetscContainer c;

518:       PetscContainerCreate(comm,&c);
519:       PetscContainerSetPointer(c,ptrs[i]);
520:       PetscContainerSetUserDestroy(c,hypre_array_destroy);
521:       PetscObjectCompose((PetscObject)(*B),names[i],(PetscObject)c);
522:       PetscContainerDestroy(&c);
523:     }
524:   }
525:   return(0);
526: }

528: static PetscErrorCode MatAIJGetParCSR_Private(Mat A, hypre_ParCSRMatrix **hA)
529: {
530:   hypre_ParCSRMatrix *tA;
531:   hypre_CSRMatrix    *hdiag,*hoffd;
532:   Mat_SeqAIJ         *diag,*offd;
533:   PetscInt           *garray,noffd,dnnz,onnz,*row_starts,*col_starts;
534:   MPI_Comm           comm = PetscObjectComm((PetscObject)A);
535:   PetscBool          ismpiaij,isseqaij;
536:   PetscErrorCode     ierr;

539:   PetscObjectTypeCompare((PetscObject)A,MATMPIAIJ,&ismpiaij);
540:   PetscObjectBaseTypeCompare((PetscObject)A,MATSEQAIJ,&isseqaij);
541:   if (!ismpiaij && !isseqaij) SETERRQ1(comm,PETSC_ERR_SUP,"Unsupported type %s",((PetscObject)A)->type);
542:   if (ismpiaij) {
543:     Mat_MPIAIJ *a = (Mat_MPIAIJ*)(A->data);

545:     diag   = (Mat_SeqAIJ*)a->A->data;
546:     offd   = (Mat_SeqAIJ*)a->B->data;
547:     garray = a->garray;
548:     noffd  = a->B->cmap->N;
549:     dnnz   = diag->nz;
550:     onnz   = offd->nz;
551:   } else {
552:     diag   = (Mat_SeqAIJ*)A->data;
553:     offd   = NULL;
554:     garray = NULL;
555:     noffd  = 0;
556:     dnnz   = diag->nz;
557:     onnz   = 0;
558:   }

560:   /* create a temporary ParCSR */
561:   if (HYPRE_AssumedPartitionCheck()) {
562:     PetscMPIInt myid;

564:     MPI_Comm_rank(comm,&myid);
565:     row_starts = A->rmap->range + myid;
566:     col_starts = A->cmap->range + myid;
567:   } else {
568:     row_starts = A->rmap->range;
569:     col_starts = A->cmap->range;
570:   }
571:   tA = hypre_ParCSRMatrixCreate(comm,A->rmap->N,A->cmap->N,(HYPRE_Int*)row_starts,(HYPRE_Int*)col_starts,noffd,dnnz,onnz);
572:   hypre_ParCSRMatrixSetRowStartsOwner(tA,0);
573:   hypre_ParCSRMatrixSetColStartsOwner(tA,0);

575:   /* set diagonal part */
576:   hdiag = hypre_ParCSRMatrixDiag(tA);
577:   hypre_CSRMatrixI(hdiag)           = (HYPRE_Int*)diag->i;
578:   hypre_CSRMatrixJ(hdiag)           = (HYPRE_Int*)diag->j;
579:   hypre_CSRMatrixData(hdiag)        = diag->a;
580:   hypre_CSRMatrixNumNonzeros(hdiag) = diag->nz;
581:   hypre_CSRMatrixSetRownnz(hdiag);
582:   hypre_CSRMatrixSetDataOwner(hdiag,0);

584:   /* set offdiagonal part */
585:   hoffd = hypre_ParCSRMatrixOffd(tA);
586:   if (offd) {
587:     hypre_CSRMatrixI(hoffd)           = (HYPRE_Int*)offd->i;
588:     hypre_CSRMatrixJ(hoffd)           = (HYPRE_Int*)offd->j;
589:     hypre_CSRMatrixData(hoffd)        = offd->a;
590:     hypre_CSRMatrixNumNonzeros(hoffd) = offd->nz;
591:     hypre_CSRMatrixSetRownnz(hoffd);
592:     hypre_CSRMatrixSetDataOwner(hoffd,0);
593:     hypre_ParCSRMatrixSetNumNonzeros(tA);
594:     hypre_ParCSRMatrixColMapOffd(tA) = (HYPRE_Int*)garray;
595:   }
596:   *hA = tA;
597:   return(0);
598: }

600: static PetscErrorCode MatAIJRestoreParCSR_Private(Mat A, hypre_ParCSRMatrix **hA)
601: {
602:   hypre_CSRMatrix    *hdiag,*hoffd;

605:   hdiag = hypre_ParCSRMatrixDiag(*hA);
606:   hoffd = hypre_ParCSRMatrixOffd(*hA);
607:   /* set pointers to NULL before destroying tA */
608:   hypre_CSRMatrixI(hdiag)           = NULL;
609:   hypre_CSRMatrixJ(hdiag)           = NULL;
610:   hypre_CSRMatrixData(hdiag)        = NULL;
611:   hypre_CSRMatrixI(hoffd)           = NULL;
612:   hypre_CSRMatrixJ(hoffd)           = NULL;
613:   hypre_CSRMatrixData(hoffd)        = NULL;
614:   hypre_ParCSRMatrixColMapOffd(*hA) = NULL;
615:   hypre_ParCSRMatrixDestroy(*hA);
616:   *hA = NULL;
617:   return(0);
618: }

620: /* calls RAP from BoomerAMG:
621:    the resulting ParCSR will not own the column and row starts
622:    It looks like we don't need to have the diagonal entries
623:    ordered first in the rows of the diagonal part
624:    for boomerAMGBuildCoarseOperator to work */
625: static PetscErrorCode MatHYPRE_ParCSR_RAP(hypre_ParCSRMatrix *hR, hypre_ParCSRMatrix *hA,hypre_ParCSRMatrix *hP, hypre_ParCSRMatrix **hRAP)
626: {
627:   HYPRE_Int P_owns_col_starts,R_owns_row_starts;

630:   P_owns_col_starts = hypre_ParCSRMatrixOwnsColStarts(hP);
631:   R_owns_row_starts = hypre_ParCSRMatrixOwnsRowStarts(hR);
632:   PetscStackCallStandard(hypre_BoomerAMGBuildCoarseOperator,(hR,hA,hP,hRAP));
633:   PetscStackCallStandard(hypre_ParCSRMatrixSetNumNonzeros,(*hRAP));
634:   /* hypre_BoomerAMGBuildCoarseOperator steals the col_starts from P and the row_starts from R */
635:   hypre_ParCSRMatrixSetRowStartsOwner(*hRAP,0);
636:   hypre_ParCSRMatrixSetColStartsOwner(*hRAP,0);
637:   if (P_owns_col_starts) hypre_ParCSRMatrixSetColStartsOwner(hP,1);
638:   if (R_owns_row_starts) hypre_ParCSRMatrixSetRowStartsOwner(hR,1);
639:   return(0);
640: }

642: static PetscErrorCode MatPtAPNumeric_AIJ_AIJ_wHYPRE(Mat A,Mat P,Mat C)
643: {
644:   Mat                B;
645:   hypre_ParCSRMatrix *hA,*hP,*hPtAP;
646:   PetscErrorCode     ierr;

649:   MatAIJGetParCSR_Private(A,&hA);
650:   MatAIJGetParCSR_Private(P,&hP);
651:   MatHYPRE_ParCSR_RAP(hP,hA,hP,&hPtAP);
652:   MatCreateFromParCSR(hPtAP,MATAIJ,PETSC_OWN_POINTER,&B);
653:   MatHeaderMerge(C,&B);
654:   MatAIJRestoreParCSR_Private(A,&hA);
655:   MatAIJRestoreParCSR_Private(P,&hP);
656:   return(0);
657: }

659: PETSC_INTERN PetscErrorCode MatPtAPSymbolic_AIJ_AIJ_wHYPRE(Mat A,Mat P,PetscReal fill,Mat *C)
660: {

664:   MatCreate(PetscObjectComm((PetscObject)A),C);
665:   MatSetType(*C,MATAIJ);
666:   (*C)->ops->ptapnumeric = MatPtAPNumeric_AIJ_AIJ_wHYPRE;
667:   return(0);
668: }

670: static PetscErrorCode MatPtAPNumeric_AIJ_HYPRE(Mat A,Mat P,Mat C)
671: {
672:   Mat                B;
673:   Mat_HYPRE          *hP;
674:   hypre_ParCSRMatrix *hA = NULL,*Pparcsr,*ptapparcsr;
675:   HYPRE_Int          type;
676:   MPI_Comm           comm = PetscObjectComm((PetscObject)A);
677:   PetscBool          ishypre;
678:   PetscErrorCode     ierr;

681:   PetscObjectTypeCompare((PetscObject)P,MATHYPRE,&ishypre);
682:   if (!ishypre) SETERRQ1(comm,PETSC_ERR_USER,"P should be of type %s",MATHYPRE);
683:   hP = (Mat_HYPRE*)P->data;
684:   PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(hP->ij,&type));
685:   if (type != HYPRE_PARCSR) SETERRQ(comm,PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported");
686:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hP->ij,(void**)&Pparcsr));

688:   MatAIJGetParCSR_Private(A,&hA);
689:   MatHYPRE_ParCSR_RAP(Pparcsr,hA,Pparcsr,&ptapparcsr);
690:   MatAIJRestoreParCSR_Private(A,&hA);

692:   /* create temporary matrix and merge to C */
693:   MatCreateFromParCSR(ptapparcsr,((PetscObject)C)->type_name,PETSC_OWN_POINTER,&B);
694:   MatHeaderMerge(C,&B);
695:   return(0);
696: }

698: static PetscErrorCode MatPtAP_AIJ_HYPRE(Mat A,Mat P,MatReuse scall,PetscReal fill,Mat *C)
699: {

703:   if (scall == MAT_INITIAL_MATRIX) {
704:     const char *deft = MATAIJ;
705:     char       type[256];
706:     PetscBool  flg;

708:     PetscObjectOptionsBegin((PetscObject)A);
709:     PetscOptionsFList("-matptap_hypre_outtype","Matrix type",NULL,MatList,deft,type,256,&flg);
710:     PetscOptionsEnd();
711:     PetscLogEventBegin(MAT_PtAPSymbolic,A,P,0,0);
712:     MatCreate(PetscObjectComm((PetscObject)A),C);
713:     if (flg) {
714:       MatSetType(*C,type);
715:     } else {
716:       MatSetType(*C,deft);
717:     }
718:     (*C)->ops->ptapnumeric = MatPtAPNumeric_AIJ_HYPRE;
719:     PetscLogEventEnd(MAT_PtAPSymbolic,A,P,0,0);
720:   }
721:   PetscLogEventBegin(MAT_PtAPNumeric,A,P,0,0);
722:   (*(*C)->ops->ptapnumeric)(A,P,*C);
723:   PetscLogEventEnd(MAT_PtAPNumeric,A,P,0,0);
724:   return(0);
725: }

727: static PetscErrorCode MatPtAPNumeric_HYPRE_HYPRE(Mat A,Mat P,Mat C)
728: {
729:   Mat                B;
730:   hypre_ParCSRMatrix *Aparcsr,*Pparcsr,*ptapparcsr;
731:   Mat_HYPRE          *hA,*hP;
732:   PetscBool          ishypre;
733:   HYPRE_Int          type;
734:   PetscErrorCode     ierr;

737:   PetscObjectTypeCompare((PetscObject)P,MATHYPRE,&ishypre);
738:   if (!ishypre) SETERRQ1(PetscObjectComm((PetscObject)P),PETSC_ERR_USER,"P should be of type %s",MATHYPRE);
739:   PetscObjectTypeCompare((PetscObject)A,MATHYPRE,&ishypre);
740:   if (!ishypre) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"A should be of type %s",MATHYPRE);
741:   hA = (Mat_HYPRE*)A->data;
742:   hP = (Mat_HYPRE*)P->data;
743:   PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(hA->ij,&type));
744:   if (type != HYPRE_PARCSR) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported");
745:   PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(hP->ij,&type));
746:   if (type != HYPRE_PARCSR) SETERRQ(PetscObjectComm((PetscObject)P),PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported");
747:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hA->ij,(void**)&Aparcsr));
748:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hP->ij,(void**)&Pparcsr));
749:   MatHYPRE_ParCSR_RAP(Pparcsr,Aparcsr,Pparcsr,&ptapparcsr);
750:   MatCreateFromParCSR(ptapparcsr,MATHYPRE,PETSC_OWN_POINTER,&B);
751:   MatHeaderMerge(C,&B);
752:   return(0);
753: }

755: static PetscErrorCode MatPtAP_HYPRE_HYPRE(Mat A,Mat P,MatReuse scall,PetscReal fill,Mat *C)
756: {

760:   if (scall == MAT_INITIAL_MATRIX) {
761:     PetscLogEventBegin(MAT_PtAPSymbolic,A,P,0,0);
762:     MatCreate(PetscObjectComm((PetscObject)A),C);
763:     MatSetType(*C,MATHYPRE);
764:     (*C)->ops->ptapnumeric = MatPtAPNumeric_HYPRE_HYPRE;
765:     PetscLogEventEnd(MAT_PtAPSymbolic,A,P,0,0);
766:   }
767:   PetscLogEventBegin(MAT_PtAPNumeric,A,P,0,0);
768:   (*(*C)->ops->ptapnumeric)(A,P,*C);
769:   PetscLogEventEnd(MAT_PtAPNumeric,A,P,0,0);
770:   return(0);
771: }

773: /* calls hypre_ParMatmul
774:    hypre_ParMatMul uses hypre_ParMatrixCreate with the communicator of hA
775:    hypre_ParMatrixCreate does not duplicate the communicator
776:    It looks like we don't need to have the diagonal entries
777:    ordered first in the rows of the diagonal part
778:    for boomerAMGBuildCoarseOperator to work */
779: static PetscErrorCode MatHYPRE_ParCSR_MatMatMult(hypre_ParCSRMatrix *hA, hypre_ParCSRMatrix *hB, hypre_ParCSRMatrix **hAB)
780: {
782:   PetscStackPush("hypre_ParMatmul");
783:   *hAB = hypre_ParMatmul(hA,hB);
784:   PetscStackPop;
785:   return(0);
786: }

788: static PetscErrorCode MatMatMultNumeric_AIJ_AIJ_wHYPRE(Mat A,Mat B,Mat C)
789: {
790:   Mat                D;
791:   hypre_ParCSRMatrix *hA,*hB,*hAB = NULL;
792:   PetscErrorCode     ierr;

795:   MatAIJGetParCSR_Private(A,&hA);
796:   MatAIJGetParCSR_Private(B,&hB);
797:   MatHYPRE_ParCSR_MatMatMult(hA,hB,&hAB);
798:   MatCreateFromParCSR(hAB,MATAIJ,PETSC_OWN_POINTER,&D);
799:   MatHeaderMerge(C,&D);
800:   MatAIJRestoreParCSR_Private(A,&hA);
801:   MatAIJRestoreParCSR_Private(B,&hB);
802:   return(0);
803: }

805: PETSC_INTERN PetscErrorCode MatMatMultSymbolic_AIJ_AIJ_wHYPRE(Mat A,Mat B,PetscReal fill,Mat *C)
806: {

810:   MatCreate(PetscObjectComm((PetscObject)A),C);
811:   MatSetType(*C,MATAIJ);
812:   (*C)->ops->matmultnumeric = MatMatMultNumeric_AIJ_AIJ_wHYPRE;
813:   return(0);
814: }

816: static PetscErrorCode MatMatMultNumeric_HYPRE_HYPRE(Mat A,Mat B,Mat C)
817: {
818:   Mat                D;
819:   hypre_ParCSRMatrix *Aparcsr,*Bparcsr,*ABparcsr = NULL;
820:   Mat_HYPRE          *hA,*hB;
821:   PetscBool          ishypre;
822:   HYPRE_Int          type;
823:   PetscErrorCode     ierr;

826:   PetscObjectTypeCompare((PetscObject)B,MATHYPRE,&ishypre);
827:   if (!ishypre) SETERRQ1(PetscObjectComm((PetscObject)B),PETSC_ERR_USER,"B should be of type %s",MATHYPRE);
828:   PetscObjectTypeCompare((PetscObject)A,MATHYPRE,&ishypre);
829:   if (!ishypre) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"A should be of type %s",MATHYPRE);
830:   hA = (Mat_HYPRE*)A->data;
831:   hB = (Mat_HYPRE*)B->data;
832:   PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(hA->ij,&type));
833:   if (type != HYPRE_PARCSR) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported");
834:   PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(hB->ij,&type));
835:   if (type != HYPRE_PARCSR) SETERRQ(PetscObjectComm((PetscObject)B),PETSC_ERR_SUP,"Only HYPRE_PARCSR is supported");
836:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hA->ij,(void**)&Aparcsr));
837:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hB->ij,(void**)&Bparcsr));
838:   MatHYPRE_ParCSR_MatMatMult(Aparcsr,Bparcsr,&ABparcsr);
839:   MatCreateFromParCSR(ABparcsr,MATHYPRE,PETSC_OWN_POINTER,&D);
840:   /* need to use HeaderReplace because HeaderMerge messes up with the communicator */
841:   MatHeaderReplace(C,&D);
842:   C->ops->matmultnumeric = MatMatMultNumeric_HYPRE_HYPRE;
843:   return(0);
844: }

846: static PetscErrorCode MatMatMult_HYPRE_HYPRE(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C)
847: {

851:   if (scall == MAT_INITIAL_MATRIX) {
852:     PetscLogEventBegin(MAT_MatMultSymbolic,A,B,0,0);
853:     MatCreate(PetscObjectComm((PetscObject)A),C);
854:     MatSetType(*C,MATHYPRE);
855:     (*C)->ops->matmultnumeric = MatMatMultNumeric_HYPRE_HYPRE;
856:     PetscLogEventEnd(MAT_MatMultSymbolic,A,B,0,0);
857:   }
858:   PetscLogEventBegin(MAT_MatMultNumeric,A,B,0,0);
859:   (*(*C)->ops->matmultnumeric)(A,B,*C);
860:   PetscLogEventEnd(MAT_MatMultNumeric,A,B,0,0);
861:   return(0);
862: }

864: PETSC_INTERN PetscErrorCode MatTransposeMatMatMultNumeric_AIJ_AIJ_AIJ_wHYPRE(Mat A,Mat B,Mat C,Mat D)
865: {
866:   Mat                E;
867:   hypre_ParCSRMatrix *hA,*hB,*hC,*hABC;
868:   PetscErrorCode     ierr;

871:   MatAIJGetParCSR_Private(A,&hA);
872:   MatAIJGetParCSR_Private(B,&hB);
873:   MatAIJGetParCSR_Private(C,&hC);
874:   MatHYPRE_ParCSR_RAP(hA,hB,hC,&hABC);
875:   MatCreateFromParCSR(hABC,MATAIJ,PETSC_OWN_POINTER,&E);
876:   MatHeaderMerge(D,&E);
877:   MatAIJRestoreParCSR_Private(A,&hA);
878:   MatAIJRestoreParCSR_Private(B,&hB);
879:   MatAIJRestoreParCSR_Private(C,&hC);
880:   return(0);
881: }

883: PETSC_INTERN PetscErrorCode MatTransposeMatMatMultSymbolic_AIJ_AIJ_AIJ_wHYPRE(Mat A,Mat B,Mat C,PetscReal fill,Mat *D)
884: {

888:   MatCreate(PetscObjectComm((PetscObject)A),D);
889:   MatSetType(*D,MATAIJ);
890:   return(0);
891: }

893: static PetscErrorCode MatMultTranspose_HYPRE(Mat A, Vec x, Vec y)
894: {

898:   MatHYPRE_MultKernel_Private(A,1.0,x,0.0,y,PETSC_TRUE);
899:   return(0);
900: }

902: static PetscErrorCode MatMult_HYPRE(Mat A, Vec x, Vec y)
903: {

907:   MatHYPRE_MultKernel_Private(A,1.0,x,0.0,y,PETSC_FALSE);
908:   return(0);
909: }

911: static PetscErrorCode MatMultAdd_HYPRE(Mat A, Vec x, Vec y, Vec z)
912: {

916:   if (y != z) {
917:     VecCopy(y,z);
918:   }
919:   MatHYPRE_MultKernel_Private(A,1.0,x,1.0,z,PETSC_FALSE);
920:   return(0);
921: }

923: static PetscErrorCode MatMultTransposeAdd_HYPRE(Mat A, Vec x, Vec y, Vec z)
924: {

928:   if (y != z) {
929:     VecCopy(y,z);
930:   }
931:   MatHYPRE_MultKernel_Private(A,1.0,x,1.0,z,PETSC_TRUE);
932:   return(0);
933: }

935: /* y = a * A * x + b * y or y = a * A^t * x + b * y depending on trans */
936: static PetscErrorCode MatHYPRE_MultKernel_Private(Mat A, PetscScalar a, Vec x, PetscScalar b, Vec y, PetscBool trans)
937: {
938:   Mat_HYPRE          *hA = (Mat_HYPRE*)A->data;
939:   hypre_ParCSRMatrix *parcsr;
940:   hypre_ParVector    *hx,*hy;
941:   PetscScalar        *ax,*ay,*sax,*say;
942:   PetscErrorCode     ierr;

945:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hA->ij,(void**)&parcsr));
946:   PetscStackCallStandard(HYPRE_IJVectorGetObject,(hA->x,(void**)&hx));
947:   PetscStackCallStandard(HYPRE_IJVectorGetObject,(hA->b,(void**)&hy));
948:   VecGetArrayRead(x,(const PetscScalar**)&ax);
949:   VecGetArray(y,&ay);
950:   if (trans) {
951:     VecHYPRE_ParVectorReplacePointer(hA->x,ay,say);
952:     VecHYPRE_ParVectorReplacePointer(hA->b,ax,sax);
953:     hypre_ParCSRMatrixMatvecT(a,parcsr,hy,b,hx);
954:     VecHYPRE_ParVectorReplacePointer(hA->x,say,ay);
955:     VecHYPRE_ParVectorReplacePointer(hA->b,sax,ax);
956:   } else {
957:     VecHYPRE_ParVectorReplacePointer(hA->x,ax,sax);
958:     VecHYPRE_ParVectorReplacePointer(hA->b,ay,say);
959:     hypre_ParCSRMatrixMatvec(a,parcsr,hx,b,hy);
960:     VecHYPRE_ParVectorReplacePointer(hA->x,sax,ax);
961:     VecHYPRE_ParVectorReplacePointer(hA->b,say,ay);
962:   }
963:   VecRestoreArrayRead(x,(const PetscScalar**)&ax);
964:   VecRestoreArray(y,&ay);
965:   return(0);
966: }

968: static PetscErrorCode MatDestroy_HYPRE(Mat A)
969: {
970:   Mat_HYPRE      *hA = (Mat_HYPRE*)A->data;

974:   if (hA->x) PetscStackCallStandard(HYPRE_IJVectorDestroy,(hA->x));
975:   if (hA->b) PetscStackCallStandard(HYPRE_IJVectorDestroy,(hA->b));
976:   if (hA->ij) {
977:     if (!hA->inner_free) hypre_IJMatrixObject(hA->ij) = NULL;
978:     PetscStackCallStandard(HYPRE_IJMatrixDestroy,(hA->ij));
979:   }
980:   if (hA->comm) { MPI_Comm_free(&hA->comm); }

982:   MatStashDestroy_Private(&A->stash);

984:   PetscFree(hA->array);

986:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_hypre_aij_C",NULL);
987:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_hypre_is_C",NULL);
988:   PetscObjectComposeFunction((PetscObject)A,"MatPtAP_seqaij_hypre_C",NULL);
989:   PetscObjectComposeFunction((PetscObject)A,"MatPtAP_mpiaij_hypre_C",NULL);
990:   PetscObjectComposeFunction((PetscObject)A,"MatHYPRESetPreallocation_C",NULL);
991:   PetscObjectComposeFunction((PetscObject)A,"MatHYPREGetParCSR_C",NULL);
992:   PetscObjectComposeFunction((PetscObject)A,"MatPtAP_is_hypre_C",NULL);
993:   PetscFree(A->data);
994:   return(0);
995: }

997: static PetscErrorCode MatSetUp_HYPRE(Mat A)
998: {

1002:   MatHYPRESetPreallocation(A,PETSC_DEFAULT,NULL,PETSC_DEFAULT,NULL);
1003:   return(0);
1004: }

1006: static PetscErrorCode MatAssemblyEnd_HYPRE(Mat A, MatAssemblyType mode)
1007: {
1008:   Mat_HYPRE          *hA = (Mat_HYPRE*)A->data;
1009:   Vec                x,b;
1010:   PetscMPIInt        n;
1011:   PetscInt           i,j,rstart,ncols,flg;
1012:   PetscInt           *row,*col;
1013:   PetscScalar        *val;
1014:   PetscErrorCode     ierr;

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

1019:   if (!A->nooffprocentries) {
1020:     while (1) {
1021:       MatStashScatterGetMesg_Private(&A->stash,&n,&row,&col,&val,&flg);
1022:       if (!flg) break;

1024:       for (i=0; i<n; ) {
1025:         /* Now identify the consecutive vals belonging to the same row */
1026:         for (j=i,rstart=row[j]; j<n; j++) {
1027:           if (row[j] != rstart) break;
1028:         }
1029:         if (j < n) ncols = j-i;
1030:         else       ncols = n-i;
1031:         /* Now assemble all these values with a single function call */
1032:         MatSetValues_HYPRE(A,1,row+i,ncols,col+i,val+i,A->insertmode);

1034:         i = j;
1035:       }
1036:     }
1037:     MatStashScatterEnd_Private(&A->stash);
1038:   }

1040:   PetscStackCallStandard(HYPRE_IJMatrixAssemble,(hA->ij));
1041:   /* The assembly routine destroys the aux_matrix, we recreate it here by calling HYPRE_IJMatrixInitalize */
1042:   {
1043:     hypre_AuxParCSRMatrix *aux_matrix;

1045:     /* call destroy just to make sure we do not leak anything */
1046:     aux_matrix = (hypre_AuxParCSRMatrix*)hypre_IJMatrixTranslator(hA->ij);
1047:     PetscStackCallStandard(hypre_AuxParCSRMatrixDestroy,(aux_matrix));
1048:     hypre_IJMatrixTranslator(hA->ij) = NULL;

1050:     /* Initialize with assembled flag -> it only recreates the aux_par_matrix */
1051:     PetscStackCallStandard(HYPRE_IJMatrixInitialize,(hA->ij));
1052:     aux_matrix = (hypre_AuxParCSRMatrix*)hypre_IJMatrixTranslator(hA->ij);
1053:     hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 1; /* see comment in MatHYPRESetPreallocation_HYPRE */
1054:     PetscStackCallStandard(hypre_AuxParCSRMatrixInitialize,(aux_matrix));
1055:   }
1056:   if (hA->x) return(0);
1057:   PetscLayoutSetUp(A->rmap);
1058:   PetscLayoutSetUp(A->cmap);
1059:   VecCreateMPIWithArray(PetscObjectComm((PetscObject)A),1,A->cmap->n,A->cmap->N,NULL,&x);
1060:   VecCreateMPIWithArray(PetscObjectComm((PetscObject)A),1,A->rmap->n,A->rmap->N,NULL,&b);
1061:   VecHYPRE_IJVectorCreate(x,&hA->x);
1062:   VecHYPRE_IJVectorCreate(b,&hA->b);
1063:   VecDestroy(&x);
1064:   VecDestroy(&b);
1065:   return(0);
1066: }

1068: static PetscErrorCode MatGetArray_HYPRE(Mat A, PetscInt size, void **array)
1069: {
1070:   Mat_HYPRE          *hA = (Mat_HYPRE*)A->data;
1071:   PetscErrorCode     ierr;

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

1076:   if (hA->size >= size) *array = hA->array;
1077:   else {
1078:     PetscFree(hA->array);
1079:     hA->size = size;
1080:     PetscMalloc(hA->size,&hA->array);
1081:     *array = hA->array;
1082:   }

1084:   hA->available = PETSC_FALSE;
1085:   return(0);
1086: }

1088: static PetscErrorCode MatRestoreArray_HYPRE(Mat A, void **array)
1089: {
1090:   Mat_HYPRE          *hA = (Mat_HYPRE*)A->data;

1093:   *array = NULL;
1094:   hA->available = PETSC_TRUE;
1095:   return(0);
1096: }

1098: PetscErrorCode MatSetValues_HYPRE(Mat A, PetscInt nr, const PetscInt rows[], PetscInt nc, const PetscInt cols[], const PetscScalar v[], InsertMode ins)
1099: {
1100:   Mat_HYPRE          *hA = (Mat_HYPRE*)A->data;
1101:   PetscScalar        *vals = (PetscScalar *)v;
1102:   PetscScalar        *sscr;
1103:   PetscInt           *cscr[2];
1104:   PetscInt           i,nzc;
1105:   void               *array = NULL;
1106:   PetscErrorCode     ierr;

1109:   MatGetArray_HYPRE(A,sizeof(PetscInt)*(2*nc)+sizeof(PetscScalar)*nc*nr,&array);
1110:   cscr[0] = (PetscInt*)array;
1111:   cscr[1] = ((PetscInt*)array)+nc;
1112:   sscr = (PetscScalar*)(((PetscInt*)array)+nc*2);
1113:   for (i=0,nzc=0;i<nc;i++) {
1114:     if (cols[i] >= 0) {
1115:       cscr[0][nzc  ] = cols[i];
1116:       cscr[1][nzc++] = i;
1117:     }
1118:   }
1119:   if (!nzc) {
1120:     MatRestoreArray_HYPRE(A,&array);
1121:     return(0);
1122:   }

1124:   if (ins == ADD_VALUES) {
1125:     for (i=0;i<nr;i++) {
1126:       if (rows[i] >= 0 && nzc) {
1127:         PetscInt j;
1128:         for (j=0;j<nzc;j++) sscr[j] = vals[cscr[1][j]];
1129:         PetscStackCallStandard(HYPRE_IJMatrixAddToValues,(hA->ij,1,(HYPRE_Int*)&nzc,(HYPRE_Int*)(rows+i),(HYPRE_Int*)cscr[0],sscr));
1130:       }
1131:       vals += nc;
1132:     }
1133:   } else { /* INSERT_VALUES */

1135:     PetscInt rst,ren;
1136:     MatGetOwnershipRange(A,&rst,&ren);

1138:     for (i=0;i<nr;i++) {
1139:       if (rows[i] >= 0 && nzc) {
1140:         PetscInt j;
1141:         for (j=0;j<nzc;j++) sscr[j] = vals[cscr[1][j]];
1142:         /* nonlocal values */
1143:         if (rows[i] < rst || rows[i] >= ren) { MatStashValuesRow_Private(&A->stash,rows[i],nzc,cscr[0],sscr,PETSC_FALSE); }
1144:         /* local values */
1145:         else PetscStackCallStandard(HYPRE_IJMatrixSetValues,(hA->ij,1,(HYPRE_Int*)&nzc,(HYPRE_Int*)(rows+i),(HYPRE_Int*)cscr[0],sscr));
1146:       }
1147:       vals += nc;
1148:     }
1149:   }

1151:   MatRestoreArray_HYPRE(A,&array);
1152:   return(0);
1153: }

1155: static PetscErrorCode MatHYPRESetPreallocation_HYPRE(Mat A, PetscInt dnz, const PetscInt dnnz[], PetscInt onz, const PetscInt onnz[])
1156: {
1157:   Mat_HYPRE      *hA = (Mat_HYPRE*)A->data;
1158:   HYPRE_Int      *hdnnz,*honnz;
1159:   PetscInt       i,rs,re,cs,ce,bs;
1160:   PetscMPIInt    size;

1164:   MatGetBlockSize(A,&bs);
1165:   PetscLayoutSetUp(A->rmap);
1166:   PetscLayoutSetUp(A->cmap);
1167:   rs   = A->rmap->rstart;
1168:   re   = A->rmap->rend;
1169:   cs   = A->cmap->rstart;
1170:   ce   = A->cmap->rend;
1171:   if (!hA->ij) {
1172:     PetscStackCallStandard(HYPRE_IJMatrixCreate,(hA->comm,rs,re-1,cs,ce-1,&hA->ij));
1173:     PetscStackCallStandard(HYPRE_IJMatrixSetObjectType,(hA->ij,HYPRE_PARCSR));
1174:   } else {
1175:     HYPRE_Int hrs,hre,hcs,hce;
1176:     PetscStackCallStandard(HYPRE_IJMatrixGetLocalRange,(hA->ij,&hrs,&hre,&hcs,&hce));
1177:     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);
1178:     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);
1179:   }
1180:   if (dnz == PETSC_DEFAULT || dnz == PETSC_DECIDE) dnz = 10*bs;
1181:   if (onz == PETSC_DEFAULT || onz == PETSC_DECIDE) onz = 10*bs;

1183:   if (!dnnz) {
1184:     PetscMalloc1(A->rmap->n,&hdnnz);
1185:     for (i=0;i<A->rmap->n;i++) hdnnz[i] = dnz;
1186:   } else {
1187:     hdnnz = (HYPRE_Int*)dnnz;
1188:   }
1189:   MPI_Comm_size(PetscObjectComm((PetscObject)A),&size);
1190:   if (size > 1) {
1191:     hypre_AuxParCSRMatrix *aux_matrix;
1192:     if (!onnz) {
1193:       PetscMalloc1(A->rmap->n,&honnz);
1194:       for (i=0;i<A->rmap->n;i++) honnz[i] = onz;
1195:     } else {
1196:       honnz = (HYPRE_Int*)onnz;
1197:     }
1198:     /* SetDiagOffdSizes sets hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 0, since it seems
1199:        they assume the user will input the entire row values, properly sorted
1200:        In PETSc, we don't make such an assumption, and we instead set this flag to 1
1201:        Also, to avoid possible memory leaks, we destroy and recreate the translator
1202:        This has to be done here, as HYPRE_IJMatrixInitialize will properly initialize
1203:        the IJ matrix for us */
1204:     aux_matrix = (hypre_AuxParCSRMatrix*)hypre_IJMatrixTranslator(hA->ij);
1205:     hypre_AuxParCSRMatrixDestroy(aux_matrix);
1206:     hypre_IJMatrixTranslator(hA->ij) = NULL;
1207:     PetscStackCallStandard(HYPRE_IJMatrixSetDiagOffdSizes,(hA->ij,hdnnz,honnz));
1208:     aux_matrix = (hypre_AuxParCSRMatrix*)hypre_IJMatrixTranslator(hA->ij);
1209:     hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 1;
1210:   } else {
1211:     honnz = NULL;
1212:     PetscStackCallStandard(HYPRE_IJMatrixSetRowSizes,(hA->ij,hdnnz));
1213:   }

1215:   /* reset assembled flag and call the initialize method */
1216:   hypre_IJMatrixAssembleFlag(hA->ij) = 0;
1217:   PetscStackCallStandard(HYPRE_IJMatrixInitialize,(hA->ij));
1218:   if (!dnnz) {
1219:     PetscFree(hdnnz);
1220:   }
1221:   if (!onnz && honnz) {
1222:     PetscFree(honnz);
1223:   }

1225:   /* Match AIJ logic */
1226:   A->preallocated = PETSC_TRUE;
1227:   A->assembled    = PETSC_FALSE;
1228:   return(0);
1229: }

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

1234:    Collective on Mat

1236:    Input Parameters:
1237: +  A - the matrix
1238: .  dnz  - number of nonzeros per row in DIAGONAL portion of local submatrix
1239:           (same value is used for all local rows)
1240: .  dnnz - array containing the number of nonzeros in the various rows of the
1241:           DIAGONAL portion of the local submatrix (possibly different for each row)
1242:           or NULL (PETSC_NULL_INTEGER in Fortran), if d_nz is used to specify the nonzero structure.
1243:           The size of this array is equal to the number of local rows, i.e 'm'.
1244:           For matrices that will be factored, you must leave room for (and set)
1245:           the diagonal entry even if it is zero.
1246: .  onz  - number of nonzeros per row in the OFF-DIAGONAL portion of local
1247:           submatrix (same value is used for all local rows).
1248: -  onnz - array containing the number of nonzeros in the various rows of the
1249:           OFF-DIAGONAL portion of the local submatrix (possibly different for
1250:           each row) or NULL (PETSC_NULL_INTEGER in Fortran), if o_nz is used to specify the nonzero
1251:           structure. The size of this array is equal to the number
1252:           of local rows, i.e 'm'.

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

1257:    Level: intermediate

1259: .keywords: matrix, aij, compressed row, sparse, parallel

1261: .seealso: MatCreate(), MatMPIAIJSetPreallocation(), MATHYPRE
1262: @*/
1263: PetscErrorCode MatHYPRESetPreallocation(Mat A, PetscInt dnz, const PetscInt dnnz[], PetscInt onz, const PetscInt onnz[])
1264: {

1270:   PetscTryMethod(A,"MatHYPRESetPreallocation_C",(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt[]),(A,dnz,dnnz,onz,onnz));
1271:   return(0);
1272: }

1274: /*
1275:    MatCreateFromParCSR - Creates a matrix from a hypre_ParCSRMatrix

1277:    Collective

1279:    Input Parameters:
1280: +  parcsr   - the pointer to the hypre_ParCSRMatrix
1281: .  mtype    - matrix type to be created. Currently MATAIJ, MATIS and MATHYPRE are supported.
1282: -  copymode - PETSc copying options

1284:    Output Parameter:
1285: .  A  - the matrix

1287:    Level: intermediate

1289: .seealso: MatHYPRE, PetscCopyMode
1290: */
1291: PETSC_EXTERN PetscErrorCode MatCreateFromParCSR(hypre_ParCSRMatrix *parcsr, MatType mtype, PetscCopyMode copymode, Mat* A)
1292: {
1293:   Mat                   T;
1294:   Mat_HYPRE             *hA;
1295:   MPI_Comm              comm;
1296:   PetscInt              rstart,rend,cstart,cend,M,N;
1297:   PetscBool             isseqaij,ismpiaij,isaij,ishyp,isis;
1298:   PetscErrorCode        ierr;

1301:   comm   = hypre_ParCSRMatrixComm(parcsr);
1302:   PetscStrcmp(mtype,MATSEQAIJ,&isseqaij);
1303:   PetscStrcmp(mtype,MATMPIAIJ,&ismpiaij);
1304:   PetscStrcmp(mtype,MATAIJ,&isaij);
1305:   PetscStrcmp(mtype,MATHYPRE,&ishyp);
1306:   PetscStrcmp(mtype,MATIS,&isis);
1307:   isaij  = (PetscBool)(isseqaij || ismpiaij || isaij);
1308:   if (!isaij && !ishyp && !isis) SETERRQ6(comm,PETSC_ERR_SUP,"Unsupported MatType %s! Supported types are %s, %s, %s, %s, and %s",mtype,MATAIJ,MATSEQAIJ,MATMPIAIJ,MATIS,MATHYPRE);
1309:   /* access ParCSRMatrix */
1310:   rstart = hypre_ParCSRMatrixFirstRowIndex(parcsr);
1311:   rend   = hypre_ParCSRMatrixLastRowIndex(parcsr);
1312:   cstart = hypre_ParCSRMatrixFirstColDiag(parcsr);
1313:   cend   = hypre_ParCSRMatrixLastColDiag(parcsr);
1314:   M      = hypre_ParCSRMatrixGlobalNumRows(parcsr);
1315:   N      = hypre_ParCSRMatrixGlobalNumCols(parcsr);

1317:   /* fix for empty local rows/columns */
1318:   if (rend < rstart) rend = rstart;
1319:   if (cend < cstart) cend = cstart;

1321:   /* PETSc convention */
1322:   rend++;
1323:   cend++;
1324:   rend = PetscMin(rend,M);
1325:   cend = PetscMin(cend,N);

1327:   /* create PETSc matrix with MatHYPRE */
1328:   MatCreate(comm,&T);
1329:   MatSetSizes(T,rend-rstart,cend-cstart,M,N);
1330:   MatSetType(T,MATHYPRE);
1331:   hA   = (Mat_HYPRE*)(T->data);

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

1337:   /* create new ParCSR object if needed */
1338:   if (ishyp && copymode == PETSC_COPY_VALUES) {
1339:     hypre_ParCSRMatrix *new_parcsr;
1340:     hypre_CSRMatrix    *hdiag,*hoffd,*ndiag,*noffd;

1342:     new_parcsr = hypre_ParCSRMatrixCompleteClone(parcsr);
1343:     hdiag      = hypre_ParCSRMatrixDiag(parcsr);
1344:     hoffd      = hypre_ParCSRMatrixOffd(parcsr);
1345:     ndiag      = hypre_ParCSRMatrixDiag(new_parcsr);
1346:     noffd      = hypre_ParCSRMatrixOffd(new_parcsr);
1347:     PetscMemcpy(hypre_CSRMatrixData(ndiag),hypre_CSRMatrixData(hdiag),hypre_CSRMatrixNumNonzeros(hdiag)*sizeof(PetscScalar));
1348:     PetscMemcpy(hypre_CSRMatrixData(noffd),hypre_CSRMatrixData(hoffd),hypre_CSRMatrixNumNonzeros(hoffd)*sizeof(PetscScalar));
1349:     parcsr     = new_parcsr;
1350:     copymode   = PETSC_OWN_POINTER;
1351:   }

1353:   /* set ParCSR object */
1354:   hypre_IJMatrixObject(hA->ij) = parcsr;
1355:   T->preallocated = PETSC_TRUE;

1357:   /* set assembled flag */
1358:   hypre_IJMatrixAssembleFlag(hA->ij) = 1;
1359:   PetscStackCallStandard(HYPRE_IJMatrixInitialize,(hA->ij));
1360:   if (ishyp) {
1361:     PetscMPIInt myid = 0;

1363:     /* make sure we always have row_starts and col_starts available */
1364:     if (HYPRE_AssumedPartitionCheck()) {
1365:       MPI_Comm_rank(comm,&myid);
1366:     }
1367:     if (!hypre_ParCSRMatrixOwnsColStarts(parcsr)) {
1368:       PetscLayout map;

1370:       MatGetLayouts(T,NULL,&map);
1371:       PetscLayoutSetUp(map);
1372:       hypre_ParCSRMatrixColStarts(parcsr) = (HYPRE_Int*)(map->range + myid);
1373:     }
1374:     if (!hypre_ParCSRMatrixOwnsRowStarts(parcsr)) {
1375:       PetscLayout map;

1377:       MatGetLayouts(T,&map,NULL);
1378:       PetscLayoutSetUp(map);
1379:       hypre_ParCSRMatrixRowStarts(parcsr) = (HYPRE_Int*)(map->range + myid);
1380:     }
1381:     /* prevent from freeing the pointer */
1382:     if (copymode == PETSC_USE_POINTER) hA->inner_free = PETSC_FALSE;
1383:     *A   = T;
1384:     MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);
1385:     MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);
1386:   } else if (isaij) {
1387:     if (copymode != PETSC_OWN_POINTER) {
1388:       /* prevent from freeing the pointer */
1389:       hA->inner_free = PETSC_FALSE;
1390:       MatConvert_HYPRE_AIJ(T,MATAIJ,MAT_INITIAL_MATRIX,A);
1391:       MatDestroy(&T);
1392:     } else { /* AIJ return type with PETSC_OWN_POINTER */
1393:       MatConvert_HYPRE_AIJ(T,MATAIJ,MAT_INPLACE_MATRIX,&T);
1394:       *A   = T;
1395:     }
1396:   } else if (isis) {
1397:     MatConvert_HYPRE_IS(T,MATIS,MAT_INITIAL_MATRIX,A);
1398:     if (copymode != PETSC_OWN_POINTER) hA->inner_free = PETSC_FALSE;
1399:     MatDestroy(&T);
1400:   }
1401:   return(0);
1402: }

1404: static PetscErrorCode MatHYPREGetParCSR_HYPRE(Mat A, hypre_ParCSRMatrix **parcsr)
1405: {
1406:   Mat_HYPRE *hA = (Mat_HYPRE*)A->data;
1407:   HYPRE_Int type;

1410:   if (!hA->ij) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"HYPRE_IJMatrix not present");
1411:   PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(hA->ij,&type));
1412:   if (type != HYPRE_PARCSR) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"HYPRE_IJMatrix is not of type HYPRE_PARCSR");
1413:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hA->ij,(void**)parcsr));
1414:   return(0);
1415: }

1417: /*
1418:    MatHYPREGetParCSR - Gets the pointer to the ParCSR matrix

1420:    Not collective

1422:    Input Parameters:
1423: +  A  - the MATHYPRE object

1425:    Output Parameter:
1426: .  parcsr  - the pointer to the hypre_ParCSRMatrix

1428:    Level: intermediate

1430: .seealso: MatHYPRE, PetscCopyMode
1431: */
1432: PetscErrorCode MatHYPREGetParCSR(Mat A, hypre_ParCSRMatrix **parcsr)
1433: {

1439:   PetscUseMethod(A,"MatHYPREGetParCSR_C",(Mat,hypre_ParCSRMatrix**),(A,parcsr));
1440:   return(0);
1441: }

1443: static PetscErrorCode MatMissingDiagonal_HYPRE(Mat A, PetscBool *missing, PetscInt *dd)
1444: {
1445:   hypre_ParCSRMatrix *parcsr;
1446:   hypre_CSRMatrix    *ha;
1447:   PetscInt           rst;
1448:   PetscErrorCode     ierr;

1451:   if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not implemented with non-square diagonal blocks");
1452:   MatGetOwnershipRange(A,&rst,NULL);
1453:   MatHYPREGetParCSR_HYPRE(A,&parcsr);
1454:   if (missing) *missing = PETSC_FALSE;
1455:   if (dd) *dd = -1;
1456:   ha = hypre_ParCSRMatrixDiag(parcsr);
1457:   if (ha) {
1458:     PetscInt  size,i;
1459:     HYPRE_Int *ii,*jj;

1461:     size = hypre_CSRMatrixNumRows(ha);
1462:     ii   = hypre_CSRMatrixI(ha);
1463:     jj   = hypre_CSRMatrixJ(ha);
1464:     for (i = 0; i < size; i++) {
1465:       PetscInt  j;
1466:       PetscBool found = PETSC_FALSE;

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

1471:       if (!found) {
1472:         PetscInfo1(A,"Matrix is missing local diagonal entry %D\n",i);
1473:         if (missing) *missing = PETSC_TRUE;
1474:         if (dd) *dd = i+rst;
1475:         return(0);
1476:       }
1477:     }
1478:     if (!size) {
1479:       PetscInfo(A,"Matrix has no diagonal entries therefore is missing diagonal\n");
1480:       if (missing) *missing = PETSC_TRUE;
1481:       if (dd) *dd = rst;
1482:     }
1483:   } else {
1484:     PetscInfo(A,"Matrix has no diagonal entries therefore is missing diagonal\n");
1485:     if (missing) *missing = PETSC_TRUE;
1486:     if (dd) *dd = rst;
1487:   }
1488:   return(0);
1489: }

1491: static PetscErrorCode MatScale_HYPRE(Mat A, PetscScalar s)
1492: {
1493:   hypre_ParCSRMatrix *parcsr;
1494:   hypre_CSRMatrix    *ha;
1495:   PetscErrorCode     ierr;

1498:   MatHYPREGetParCSR_HYPRE(A,&parcsr);
1499:   /* diagonal part */
1500:   ha = hypre_ParCSRMatrixDiag(parcsr);
1501:   if (ha) {
1502:     PetscInt    size,i;
1503:     HYPRE_Int   *ii;
1504:     PetscScalar *a;

1506:     size = hypre_CSRMatrixNumRows(ha);
1507:     a    = hypre_CSRMatrixData(ha);
1508:     ii   = hypre_CSRMatrixI(ha);
1509:     for (i = 0; i < ii[size]; i++) a[i] *= s;
1510:   }
1511:   /* offdiagonal part */
1512:   ha = hypre_ParCSRMatrixOffd(parcsr);
1513:   if (ha) {
1514:     PetscInt    size,i;
1515:     HYPRE_Int   *ii;
1516:     PetscScalar *a;

1518:     size = hypre_CSRMatrixNumRows(ha);
1519:     a    = hypre_CSRMatrixData(ha);
1520:     ii   = hypre_CSRMatrixI(ha);
1521:     for (i = 0; i < ii[size]; i++) a[i] *= s;
1522:   }
1523:   return(0);
1524: }

1526: static PetscErrorCode MatZeroRowsColumns_HYPRE(Mat A, PetscInt numRows, const PetscInt rows[], PetscScalar diag, Vec x, Vec b)
1527: {
1528:   hypre_ParCSRMatrix *parcsr;
1529:   HYPRE_Int          *lrows;
1530:   PetscInt           rst,ren,i;
1531:   PetscErrorCode     ierr;

1534:   if (x || b) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"To be implemented");
1535:   MatHYPREGetParCSR_HYPRE(A,&parcsr);
1536:   PetscMalloc1(numRows,&lrows);
1537:   MatGetOwnershipRange(A,&rst,&ren);
1538:   for (i=0;i<numRows;i++) {
1539:     if (rows[i] < rst || rows[i] >= ren)
1540:       SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Non-local rows not yet supported");
1541:     lrows[i] = rows[i] - rst;
1542:   }
1543:   PetscStackCallStandard(hypre_ParCSRMatrixEliminateRowsCols,(parcsr,numRows,lrows));
1544:   PetscFree(lrows);
1545:   return(0);
1546: }

1548: static PetscErrorCode MatZeroEntries_HYPRE_CSRMatrix(hypre_CSRMatrix *ha)
1549: {
1550:   PetscErrorCode      ierr;

1553:   if (ha) {
1554:     HYPRE_Int     *ii, size;
1555:     HYPRE_Complex *a;

1557:     size = hypre_CSRMatrixNumRows(ha);
1558:     a    = hypre_CSRMatrixData(ha);
1559:     ii   = hypre_CSRMatrixI(ha);

1561:     if (a) { PetscMemzero(a,(ii[size])*sizeof(HYPRE_Complex)); }
1562:   }
1563:   return(0);
1564: }

1566: PetscErrorCode MatZeroEntries_HYPRE(Mat A)
1567: {
1568:   hypre_ParCSRMatrix  *parcsr;
1569:   PetscErrorCode      ierr;

1572:   MatHYPREGetParCSR_HYPRE(A,&parcsr);
1573:   /* diagonal part */
1574:   MatZeroEntries_HYPRE_CSRMatrix(hypre_ParCSRMatrixDiag(parcsr));
1575:   /* off-diagonal part */
1576:   MatZeroEntries_HYPRE_CSRMatrix(hypre_ParCSRMatrixOffd(parcsr));
1577:   return(0);
1578: }

1580: static PetscErrorCode MatZeroRows_HYPRE_CSRMatrix(hypre_CSRMatrix *hA,PetscInt N,const PetscInt rows[],PetscScalar diag)
1581: {
1582:   PetscInt        ii, jj, ibeg, iend, irow;
1583:   PetscInt        *i, *j;
1584:   PetscScalar     *a;

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

1589:   i = (PetscInt*) hypre_CSRMatrixI(hA);
1590:   j = (PetscInt*) hypre_CSRMatrixJ(hA);
1591:   a = hypre_CSRMatrixData(hA);

1593:   for (ii = 0; ii < N; ii++) {
1594:     irow = rows[ii];
1595:     ibeg = i[irow];
1596:     iend = i[irow+1];
1597:     for (jj = ibeg; jj < iend; jj++)
1598:       if (j[jj] == irow) a[jj] = diag;
1599:       else a[jj] = 0.0;
1600:    }
1601:    return(0);
1602: }

1604: static PetscErrorCode MatZeroRows_HYPRE(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
1605: {
1606:   hypre_ParCSRMatrix  *parcsr;
1607:   PetscInt            *lrows,len;
1608:   PetscErrorCode      ierr;

1611:   if (x || b) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Does not support to modify the solution and the right hand size \n");
1612:   /* retrieve the internal matrix */
1613:   MatHYPREGetParCSR_HYPRE(A,&parcsr);
1614:   /* get locally owned rows */
1615:   MatZeroRowsMapLocal_Private(A,N,rows,&len,&lrows);
1616:   /* zero diagonal part */
1617:   MatZeroRows_HYPRE_CSRMatrix(hypre_ParCSRMatrixDiag(parcsr),len,lrows,diag);
1618:   /* zero off-diagonal part */
1619:   MatZeroRows_HYPRE_CSRMatrix(hypre_ParCSRMatrixOffd(parcsr),len,lrows,0.0);

1621:   PetscFree(lrows);
1622:   return(0);
1623: }

1625: static PetscErrorCode MatAssemblyBegin_HYPRE(Mat mat,MatAssemblyType mode)
1626: {

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

1632:   MatStashScatterBegin_Private(mat,&mat->stash,mat->rmap->range);
1633:   return(0);
1634: }

1636: static PetscErrorCode MatGetRow_HYPRE(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v)
1637: {
1638:   hypre_ParCSRMatrix  *parcsr;
1639:   PetscErrorCode      ierr;

1642:   /* retrieve the internal matrix */
1643:   MatHYPREGetParCSR_HYPRE(A,&parcsr);
1644:   /* call HYPRE API */
1645:   PetscStackCallStandard(HYPRE_ParCSRMatrixGetRow,(parcsr,row,(HYPRE_Int*)nz,(HYPRE_Int**)idx,v));
1646:   return(0);
1647: }

1649: static PetscErrorCode MatRestoreRow_HYPRE(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v)
1650: {
1651:   hypre_ParCSRMatrix  *parcsr;
1652:   PetscErrorCode      ierr;

1655:   /* retrieve the internal matrix */
1656:   MatHYPREGetParCSR_HYPRE(A,&parcsr);
1657:   /* call HYPRE API */
1658:   PetscStackCallStandard(HYPRE_ParCSRMatrixRestoreRow,(parcsr,row,(HYPRE_Int*)nz,(HYPRE_Int**)idx,v));
1659:   return(0);
1660: }

1662: static PetscErrorCode MatGetValues_HYPRE(Mat A,PetscInt m,const PetscInt idxm[],PetscInt n,const PetscInt idxn[],PetscScalar v[])
1663: {
1664:   Mat_HYPRE *hA = (Mat_HYPRE*)A->data;
1665:   PetscInt  i;

1668:   if (!m || !n) return(0);
1669:   /* Ignore negative row indices
1670:    * And negative column indices should be automatically ignored in hypre
1671:    * */
1672:   for (i=0; i<m; i++)
1673:     if (idxm[i] >= 0) PetscStackCallStandard(HYPRE_IJMatrixGetValues,(hA->ij,1,(HYPRE_Int*)&n,(HYPRE_Int*)&idxm[i],(HYPRE_Int*)idxn,&v[i*n]));
1674:   return(0);
1675: }

1677: static PetscErrorCode MatSetOption_HYPRE(Mat A,MatOption op,PetscBool flg)
1678: {
1679:   Mat_HYPRE *hA = (Mat_HYPRE*)A->data;

1682:   switch (op)
1683:   {
1684:   case MAT_NO_OFF_PROC_ENTRIES:
1685:     if (flg) {
1686:       PetscStackCallStandard(HYPRE_IJMatrixSetMaxOffProcElmts,(hA->ij,0));
1687:     }
1688:     break;
1689:   default:
1690:     break;
1691:   }
1692:   return(0);
1693: }

1695: static PetscErrorCode MatView_HYPRE(Mat A, PetscViewer view)
1696: {
1697:   hypre_ParCSRMatrix *parcsr;
1698:   PetscErrorCode     ierr;
1699:   Mat                B;
1700:   PetscViewerFormat  format;
1701:   PetscErrorCode     (*mview)(Mat,PetscViewer) = NULL;

1704:   PetscViewerGetFormat(view,&format);
1705:   if (format != PETSC_VIEWER_NATIVE) {
1706:     MatHYPREGetParCSR_HYPRE(A,&parcsr);
1707:     MatCreateFromParCSR(parcsr,MATAIJ,PETSC_USE_POINTER,&B);
1708:     MatGetOperation(B,MATOP_VIEW,(void(**)(void))&mview);
1709:     if (!mview) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Missing view operation");
1710:     (*mview)(B,view);
1711:     MatDestroy(&B);
1712:   } else {
1713:     Mat_HYPRE  *hA = (Mat_HYPRE*)A->data;
1714:     PetscMPIInt size;
1715:     PetscBool   isascii;
1716:     const char *filename;

1718:     /* HYPRE uses only text files */
1719:     PetscObjectTypeCompare((PetscObject)view,PETSCVIEWERASCII,&isascii);
1720:     if (!isascii) SETERRQ1(PetscObjectComm((PetscObject)view),PETSC_ERR_SUP,"PetscViewerType %s: native HYPRE format needs PETSCVIEWERASCII",((PetscObject)view)->type_name);
1721:     PetscViewerFileGetName(view,&filename);
1722:     PetscStackCallStandard(HYPRE_IJMatrixPrint,(hA->ij,filename));
1723:     MPI_Comm_size(hA->comm,&size);
1724:     if (size > 1) {
1725:       PetscViewerASCIIPrintf(view,"Matrix files: %s.%05d ... %s.%05d\n",filename,0,filename,size-1);
1726:     } else {
1727:       PetscViewerASCIIPrintf(view,"Matrix file: %s.%05d\n",filename,0);
1728:     }
1729:   }
1730:   return(0);
1731: }

1733: static PetscErrorCode MatDuplicate_HYPRE(Mat A,MatDuplicateOption op, Mat *B)
1734: {
1735:   hypre_ParCSRMatrix *parcsr;
1736:   PetscErrorCode     ierr;
1737:   PetscCopyMode      cpmode;

1740:   MatHYPREGetParCSR_HYPRE(A,&parcsr);
1741:   if (op == MAT_DO_NOT_COPY_VALUES || op == MAT_SHARE_NONZERO_PATTERN) {
1742:     parcsr = hypre_ParCSRMatrixCompleteClone(parcsr);
1743:     cpmode = PETSC_OWN_POINTER;
1744:   } else {
1745:     cpmode = PETSC_COPY_VALUES;
1746:   }
1747:   MatCreateFromParCSR(parcsr,MATHYPRE,cpmode,B);
1748:   return(0);
1749: }

1751: static PetscErrorCode MatCopy_HYPRE(Mat A, Mat B, MatStructure str)
1752: {
1753:   hypre_ParCSRMatrix *acsr,*bcsr;
1754:   PetscErrorCode     ierr;

1757:   if (str == SAME_NONZERO_PATTERN && A->ops->copy == B->ops->copy) {
1758:     MatHYPREGetParCSR_HYPRE(A,&acsr);
1759:     MatHYPREGetParCSR_HYPRE(B,&bcsr);
1760:     PetscStackCallStandard(hypre_ParCSRMatrixCopy,(acsr,bcsr,1));
1761:     MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
1762:     MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
1763:   } else {
1764:     MatCopy_Basic(A,B,str);
1765:   }
1766:   return(0);
1767: }

1769: static PetscErrorCode MatGetDiagonal_HYPRE(Mat A, Vec d)
1770: {
1771:   hypre_ParCSRMatrix *parcsr;
1772:   hypre_CSRMatrix    *dmat;
1773:   PetscScalar        *a,*data = NULL;
1774:   PetscInt           i,*diag = NULL;
1775:   PetscBool          cong;
1776:   PetscErrorCode     ierr;

1779:   MatHasCongruentLayouts(A,&cong);
1780:   if (!cong) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Only for square matrices with same local distributions of rows and columns");
1781: #if defined(PETSC_USE_DEBUG)
1782:   {
1783:     PetscBool miss;
1784:     MatMissingDiagonal(A,&miss,NULL);
1785:     if (miss && A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not implemented when diagonal entries are missing");
1786:   }
1787: #endif
1788:   MatHYPREGetParCSR_HYPRE(A,&parcsr);
1789:   dmat = hypre_ParCSRMatrixDiag(parcsr);
1790:   if (dmat) {
1791:     VecGetArray(d,&a);
1792:     diag = (PetscInt*)(hypre_CSRMatrixI(dmat));
1793:     data = (PetscScalar*)(hypre_CSRMatrixData(dmat));
1794:     for (i=0;i<A->rmap->n;i++) a[i] = data[diag[i]];
1795:     VecRestoreArray(d,&a);
1796:   }
1797:   return(0);
1798: }

1800:  #include <petscblaslapack.h>

1802: static PetscErrorCode MatAXPY_HYPRE(Mat Y,PetscScalar a,Mat X,MatStructure str)
1803: {

1807:   if (str == SAME_NONZERO_PATTERN) {
1808:     hypre_ParCSRMatrix *x,*y;
1809:     hypre_CSRMatrix    *xloc,*yloc;
1810:     PetscInt           xnnz,ynnz;
1811:     PetscScalar        *xarr,*yarr;
1812:     PetscBLASInt       one=1,bnz;

1814:     MatHYPREGetParCSR(Y,&y);
1815:     MatHYPREGetParCSR(X,&x);

1817:     /* diagonal block */
1818:     xloc = hypre_ParCSRMatrixDiag(x);
1819:     yloc = hypre_ParCSRMatrixDiag(y);
1820:     xnnz = 0;
1821:     ynnz = 0;
1822:     xarr = NULL;
1823:     yarr = NULL;
1824:     if (xloc) {
1825:       xarr = (PetscScalar*)(hypre_CSRMatrixData(xloc));
1826:       xnnz = hypre_CSRMatrixNumNonzeros(xloc);
1827:     }
1828:     if (yloc) {
1829:       yarr = (PetscScalar*)(hypre_CSRMatrixData(yloc));
1830:       ynnz = hypre_CSRMatrixNumNonzeros(yloc);
1831:     }
1832:     if (xnnz != ynnz) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Different number of nonzeros in diagonal block %D != %D",xnnz,ynnz);
1833:     PetscBLASIntCast(xnnz,&bnz);
1834:     PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&bnz,&a,xarr,&one,yarr,&one));

1836:     /* off-diagonal block */
1837:     xloc = hypre_ParCSRMatrixOffd(x);
1838:     yloc = hypre_ParCSRMatrixOffd(y);
1839:     xnnz = 0;
1840:     ynnz = 0;
1841:     xarr = NULL;
1842:     yarr = NULL;
1843:     if (xloc) {
1844:       xarr = (PetscScalar*)(hypre_CSRMatrixData(xloc));
1845:       xnnz = hypre_CSRMatrixNumNonzeros(xloc);
1846:     }
1847:     if (yloc) {
1848:       yarr = (PetscScalar*)(hypre_CSRMatrixData(yloc));
1849:       ynnz = hypre_CSRMatrixNumNonzeros(yloc);
1850:     }
1851:     if (xnnz != ynnz) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Different number of nonzeros in off-diagonal block %D != %D",xnnz,ynnz);
1852:     PetscBLASIntCast(xnnz,&bnz);
1853:     PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&bnz,&a,xarr,&one,yarr,&one));
1854:   } else if (str == SUBSET_NONZERO_PATTERN) {
1855:     MatAXPY_Basic(Y,a,X,str);
1856:   } else {
1857:     Mat B;

1859:     MatAXPY_Basic_Preallocate(Y,X,&B);
1860:     MatAXPY_BasicWithPreallocation(B,Y,a,X,str);
1861:     MatHeaderReplace(Y,&B);
1862:   }
1863:   return(0);
1864: }

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

1870:    Level: intermediate

1872: .seealso: MatCreate()
1873: M*/

1875: PETSC_EXTERN PetscErrorCode MatCreate_HYPRE(Mat B)
1876: {
1877:   Mat_HYPRE      *hB;

1881:   PetscNewLog(B,&hB);
1882:   hB->inner_free = PETSC_TRUE;
1883:   hB->available  = PETSC_TRUE;
1884:   hB->size       = 0;
1885:   hB->array      = NULL;

1887:   B->data       = (void*)hB;
1888:   B->rmap->bs   = 1;
1889:   B->assembled  = PETSC_FALSE;

1891:   PetscMemzero(B->ops,sizeof(struct _MatOps));
1892:   B->ops->mult             = MatMult_HYPRE;
1893:   B->ops->multtranspose    = MatMultTranspose_HYPRE;
1894:   B->ops->multadd          = MatMultAdd_HYPRE;
1895:   B->ops->multtransposeadd = MatMultTransposeAdd_HYPRE;
1896:   B->ops->setup            = MatSetUp_HYPRE;
1897:   B->ops->destroy          = MatDestroy_HYPRE;
1898:   B->ops->assemblyend      = MatAssemblyEnd_HYPRE;
1899:   B->ops->assemblybegin    = MatAssemblyBegin_HYPRE;
1900:   B->ops->ptap             = MatPtAP_HYPRE_HYPRE;
1901:   B->ops->matmult          = MatMatMult_HYPRE_HYPRE;
1902:   B->ops->setvalues        = MatSetValues_HYPRE;
1903:   B->ops->missingdiagonal  = MatMissingDiagonal_HYPRE;
1904:   B->ops->scale            = MatScale_HYPRE;
1905:   B->ops->zerorowscolumns  = MatZeroRowsColumns_HYPRE;
1906:   B->ops->zeroentries      = MatZeroEntries_HYPRE;
1907:   B->ops->zerorows         = MatZeroRows_HYPRE;
1908:   B->ops->getrow           = MatGetRow_HYPRE;
1909:   B->ops->restorerow       = MatRestoreRow_HYPRE;
1910:   B->ops->getvalues        = MatGetValues_HYPRE;
1911:   B->ops->setoption        = MatSetOption_HYPRE;
1912:   B->ops->duplicate        = MatDuplicate_HYPRE;
1913:   B->ops->copy             = MatCopy_HYPRE;
1914:   B->ops->view             = MatView_HYPRE;
1915:   B->ops->getdiagonal      = MatGetDiagonal_HYPRE;
1916:   B->ops->axpy             = MatAXPY_HYPRE;

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

1921:   MPI_Comm_dup(PetscObjectComm((PetscObject)B),&hB->comm);
1922:   PetscObjectChangeTypeName((PetscObject)B,MATHYPRE);
1923:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_hypre_aij_C",MatConvert_HYPRE_AIJ);
1924:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_hypre_is_C",MatConvert_HYPRE_IS);
1925:   PetscObjectComposeFunction((PetscObject)B,"MatPtAP_seqaij_hypre_C",MatPtAP_AIJ_HYPRE);
1926:   PetscObjectComposeFunction((PetscObject)B,"MatPtAP_mpiaij_hypre_C",MatPtAP_AIJ_HYPRE);
1927:   PetscObjectComposeFunction((PetscObject)B,"MatHYPRESetPreallocation_C",MatHYPRESetPreallocation_HYPRE);
1928:   PetscObjectComposeFunction((PetscObject)B,"MatHYPREGetParCSR_C",MatHYPREGetParCSR_HYPRE);
1929:   PetscObjectComposeFunction((PetscObject)B,"MatPtAP_is_hypre_C",MatPtAP_IS_XAIJ);
1930:   return(0);
1931: }

1933: static PetscErrorCode hypre_array_destroy(void *ptr)
1934: {
1936:    hypre_TFree(ptr,HYPRE_MEMORY_HOST);
1937:    return(0);
1938: }