Actual source code: mhypre.c

petsc-3.10.5 2019-03-28
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,Vec,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_HYPRE      *hB;
328:   MPI_Comm       comm = PetscObjectComm((PetscObject)A);

332:   if (reuse == MAT_INPLACE_MATRIX) SETERRQ(comm,PETSC_ERR_SUP,"Unsupported MAT_INPLACE_MATRIX");
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,B);
343:     MatSetType(*B,MATHYPRE);
344:     MatSetSizes(*B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
345:     hB   = (Mat_HYPRE*)((*B)->data);
346:   }
347:   MatHYPRE_CreateFromMat(A,hB);
348:   MatHYPRE_IJMatrixCopy(A,hB->ij);
349:   (*B)->preallocated = PETSC_TRUE;
350:   MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
351:   MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
352:   return(0);
353: }

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

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

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

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

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

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

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

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

514:       PetscContainerCreate(comm,&c);
515:       PetscContainerSetPointer(c,ptrs[i]);
516:       PetscContainerSetUserDestroy(c,hypre_array_destroy);
517:       PetscObjectCompose((PetscObject)(*B),names[i],(PetscObject)c);
518:       PetscContainerDestroy(&c);
519:     }
520:   }
521:   return(0);
522: }

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

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

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

556:   /* create a temporary ParCSR */
557:   if (HYPRE_AssumedPartitionCheck()) {
558:     PetscMPIInt myid;

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

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

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

596: static PetscErrorCode MatAIJRestoreParCSR_Private(Mat A, hypre_ParCSRMatrix **hA)
597: {
598:   hypre_CSRMatrix    *hdiag,*hoffd;

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

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

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

639: static PetscErrorCode MatPtAPNumeric_AIJ_AIJ_wHYPRE(Mat A,Mat P,Mat C)
640: {
641:   Mat                B;
642:   hypre_ParCSRMatrix *hA,*hP,*hPtAP;
643:   PetscErrorCode     ierr;

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

656: PETSC_INTERN PetscErrorCode MatPtAPSymbolic_AIJ_AIJ_wHYPRE(Mat A,Mat P,PetscReal fill,Mat *C)
657: {

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

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

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

685:   MatAIJGetParCSR_Private(A,&hA);
686:   MatHYPRE_ParCSR_RAP(Pparcsr,hA,Pparcsr,&ptapparcsr);
687:   MatAIJRestoreParCSR_Private(A,&hA);

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

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

700:   if (scall == MAT_INITIAL_MATRIX) {
701:     const char *deft = MATAIJ;
702:     char       type[256];
703:     PetscBool  flg;

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

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

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

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

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

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

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

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

802: PETSC_INTERN PetscErrorCode MatMatMultSymbolic_AIJ_AIJ_wHYPRE(Mat A,Mat B,PetscReal fill,Mat *C)
803: {

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

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

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

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

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

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

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

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

885:   MatCreate(PetscObjectComm((PetscObject)A),D);
886:   MatSetType(*D,MATAIJ);
887:   return(0);
888: }

890: static PetscErrorCode MatMultTranspose_HYPRE(Mat A, Vec x, Vec y)
891: {

895:   MatHYPRE_MultKernel_Private(A,x,y,PETSC_TRUE);
896:   return(0);
897: }

899: static PetscErrorCode MatMult_HYPRE(Mat A, Vec x, Vec y)
900: {

904:   MatHYPRE_MultKernel_Private(A,x,y,PETSC_FALSE);
905:   return(0);
906: }

908: static PetscErrorCode MatHYPRE_MultKernel_Private(Mat A, Vec x, Vec y, PetscBool trans)
909: {
910:   Mat_HYPRE          *hA = (Mat_HYPRE*)A->data;
911:   hypre_ParCSRMatrix *parcsr;
912:   hypre_ParVector    *hx,*hy;
913:   PetscScalar        *ax,*ay,*sax,*say;
914:   PetscErrorCode     ierr;

917:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hA->ij,(void**)&parcsr));
918:   PetscStackCallStandard(HYPRE_IJVectorGetObject,(hA->x,(void**)&hx));
919:   PetscStackCallStandard(HYPRE_IJVectorGetObject,(hA->b,(void**)&hy));
920:   VecGetArrayRead(x,(const PetscScalar**)&ax);
921:   VecGetArray(y,&ay);
922:   if (trans) {
923:     VecHYPRE_ParVectorReplacePointer(hA->x,ay,say);
924:     VecHYPRE_ParVectorReplacePointer(hA->b,ax,sax);
925:     hypre_ParCSRMatrixMatvecT(1.,parcsr,hy,0.,hx);
926:     VecHYPRE_ParVectorReplacePointer(hA->x,say,ay);
927:     VecHYPRE_ParVectorReplacePointer(hA->b,sax,ax);
928:   } else {
929:     VecHYPRE_ParVectorReplacePointer(hA->x,ax,sax);
930:     VecHYPRE_ParVectorReplacePointer(hA->b,ay,say);
931:     hypre_ParCSRMatrixMatvec(1.,parcsr,hx,0.,hy);
932:     VecHYPRE_ParVectorReplacePointer(hA->x,sax,ax);
933:     VecHYPRE_ParVectorReplacePointer(hA->b,say,ay);
934:   }
935:   VecRestoreArrayRead(x,(const PetscScalar**)&ax);
936:   VecRestoreArray(y,&ay);
937:   return(0);
938: }

940: static PetscErrorCode MatDestroy_HYPRE(Mat A)
941: {
942:   Mat_HYPRE      *hA = (Mat_HYPRE*)A->data;

946:   if (hA->x) PetscStackCallStandard(HYPRE_IJVectorDestroy,(hA->x));
947:   if (hA->b) PetscStackCallStandard(HYPRE_IJVectorDestroy,(hA->b));
948:   if (hA->ij) {
949:     if (!hA->inner_free) hypre_IJMatrixObject(hA->ij) = NULL;
950:     PetscStackCallStandard(HYPRE_IJMatrixDestroy,(hA->ij));
951:   }
952:   if (hA->comm) { MPI_Comm_free(&hA->comm); }

954:   MatStashDestroy_Private(&A->stash);

956:   PetscFree(hA->array);

958:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_hypre_aij_C",NULL);
959:   PetscObjectComposeFunction((PetscObject)A,"MatConvert_hypre_is_C",NULL);
960:   PetscObjectComposeFunction((PetscObject)A,"MatPtAP_seqaij_hypre_C",NULL);
961:   PetscObjectComposeFunction((PetscObject)A,"MatPtAP_mpiaij_hypre_C",NULL);
962:   PetscObjectComposeFunction((PetscObject)A,"MatHYPRESetPreallocation_C",NULL);
963:   PetscObjectComposeFunction((PetscObject)A,"MatHYPREGetParCSR_C",NULL);
964:   PetscObjectComposeFunction((PetscObject)A,"MatPtAP_is_hypre_C",NULL);
965:   PetscFree(A->data);
966:   return(0);
967: }

969: static PetscErrorCode MatSetUp_HYPRE(Mat A)
970: {

974:   MatHYPRESetPreallocation(A,PETSC_DEFAULT,NULL,PETSC_DEFAULT,NULL);
975:   return(0);
976: }


979: static PetscErrorCode MatAssemblyEnd_HYPRE(Mat A, MatAssemblyType mode)
980: {
981:   Mat_HYPRE          *hA = (Mat_HYPRE*)A->data;
982:   Vec                x,b;
983:   PetscMPIInt        n;
984:   PetscInt           i,j,rstart,ncols,flg;
985:   PetscInt           *row,*col;
986:   PetscScalar        *val;
987:   PetscErrorCode     ierr;

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

992:   if (!A->nooffprocentries) {
993:     while (1) {
994:       MatStashScatterGetMesg_Private(&A->stash,&n,&row,&col,&val,&flg);
995:       if (!flg) break;

997:       for (i=0; i<n; ) {
998:         /* Now identify the consecutive vals belonging to the same row */
999:         for (j=i,rstart=row[j]; j<n; j++) {
1000:           if (row[j] != rstart) break;
1001:         }
1002:         if (j < n) ncols = j-i;
1003:         else       ncols = n-i;
1004:         /* Now assemble all these values with a single function call */
1005:         MatSetValues_HYPRE(A,1,row+i,ncols,col+i,val+i,A->insertmode);

1007:         i = j;
1008:       }
1009:     }
1010:     MatStashScatterEnd_Private(&A->stash);
1011:   }

1013:   PetscStackCallStandard(HYPRE_IJMatrixAssemble,(hA->ij));
1014:   if (hA->x) return(0);
1015:   PetscLayoutSetUp(A->rmap);
1016:   PetscLayoutSetUp(A->cmap);
1017:   VecCreateMPIWithArray(PetscObjectComm((PetscObject)A),1,A->cmap->n,A->cmap->N,NULL,&x);
1018:   VecCreateMPIWithArray(PetscObjectComm((PetscObject)A),1,A->rmap->n,A->rmap->N,NULL,&b);
1019:   VecHYPRE_IJVectorCreate(x,&hA->x);
1020:   VecHYPRE_IJVectorCreate(b,&hA->b);
1021:   VecDestroy(&x);
1022:   VecDestroy(&b);
1023:   return(0);
1024: }

1026: static PetscErrorCode MatGetArray_HYPRE(Mat A, PetscInt size, void **array)
1027: {
1028:   Mat_HYPRE          *hA = (Mat_HYPRE*)A->data;
1029:   PetscErrorCode     ierr;

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

1034:   if (hA->size >= size) *array = hA->array;
1035:   else {
1036:     PetscFree(hA->array);
1037:     hA->size = size;
1038:     PetscMalloc(hA->size,&hA->array);
1039:     *array = hA->array;
1040:   }

1042:   hA->available = PETSC_FALSE;
1043:   return(0);
1044: }

1046: static PetscErrorCode MatRestoreArray_HYPRE(Mat A, void **array)
1047: {
1048:   Mat_HYPRE          *hA = (Mat_HYPRE*)A->data;

1051:   *array = NULL;
1052:   hA->available = PETSC_TRUE;
1053:   return(0);
1054: }


1057: PetscErrorCode MatSetValues_HYPRE(Mat A, PetscInt nr, const PetscInt rows[], PetscInt nc, const PetscInt cols[], const PetscScalar v[], InsertMode ins)
1058: {
1059:   Mat_HYPRE          *hA = (Mat_HYPRE*)A->data;
1060:   PetscScalar        *vals = (PetscScalar *)v;
1061:   PetscScalar        *sscr;
1062:   PetscInt           *cscr[2];
1063:   PetscInt           i,nzc;
1064:   void               *array = NULL;
1065:   PetscErrorCode     ierr;

1068:   MatGetArray_HYPRE(A,sizeof(PetscInt)*(2*nc)+sizeof(PetscScalar)*nc*nr,&array);
1069:   cscr[0] = (PetscInt*)array;
1070:   cscr[1] = ((PetscInt*)array)+nc;
1071:   sscr = (PetscScalar*)(((PetscInt*)array)+nc*2);
1072:   for (i=0,nzc=0;i<nc;i++) {
1073:     if (cols[i] >= 0) {
1074:       cscr[0][nzc  ] = cols[i];
1075:       cscr[1][nzc++] = i;
1076:     }
1077:   }
1078:   if (!nzc) {
1079:     MatRestoreArray_HYPRE(A,&array);
1080:     return(0);
1081:   }

1083:   if (ins == ADD_VALUES) {
1084:     for (i=0;i<nr;i++) {
1085:       if (rows[i] >= 0 && nzc) {
1086:         PetscInt j;
1087:         for (j=0;j<nzc;j++) sscr[j] = vals[cscr[1][j]];
1088:         PetscStackCallStandard(HYPRE_IJMatrixAddToValues,(hA->ij,1,(HYPRE_Int*)&nzc,(HYPRE_Int*)(rows+i),(HYPRE_Int*)cscr[0],sscr));
1089:       }
1090:       vals += nc;
1091:     }
1092:   } else { /* INSERT_VALUES */

1094:     PetscInt rst,ren;
1095:     MatGetOwnershipRange(A,&rst,&ren);

1097:     for (i=0;i<nr;i++) {
1098:       if (rows[i] >= 0 && nzc) {
1099:         PetscInt j;
1100:         for (j=0;j<nzc;j++) sscr[j] = vals[cscr[1][j]];
1101:         /* nonlocal values */
1102:         if (rows[i] < rst || rows[i] >= ren) { MatStashValuesRow_Private(&A->stash,rows[i],nzc,cscr[0],sscr, PETSC_FALSE); }
1103:         /* local values */
1104:         else PetscStackCallStandard(HYPRE_IJMatrixSetValues,(hA->ij,1,(HYPRE_Int*)&nzc,(HYPRE_Int*)(rows+i),(HYPRE_Int*)cscr[0],sscr));
1105:       }
1106:       vals += nc;
1107:     }
1108:   }

1110:   MatRestoreArray_HYPRE(A,&array);
1111:   return(0);
1112: }

1114: static PetscErrorCode MatHYPRESetPreallocation_HYPRE(Mat A, PetscInt dnz, const PetscInt dnnz[], PetscInt onz, const PetscInt onnz[])
1115: {
1116:   Mat_HYPRE          *hA = (Mat_HYPRE*)A->data;
1117:   HYPRE_Int          *hdnnz,*honnz;
1118:   PetscInt           i,rs,re,cs,ce,bs;
1119:   PetscMPIInt        size;
1120:   PetscErrorCode     ierr;

1123:   MatGetBlockSize(A,&bs);
1124:   PetscLayoutSetUp(A->rmap);
1125:   PetscLayoutSetUp(A->cmap);
1126:   rs   = A->rmap->rstart;
1127:   re   = A->rmap->rend;
1128:   cs   = A->cmap->rstart;
1129:   ce   = A->cmap->rend;
1130:   if (!hA->ij) {
1131:     PetscStackCallStandard(HYPRE_IJMatrixCreate,(hA->comm,rs,re-1,cs,ce-1,&hA->ij));
1132:     PetscStackCallStandard(HYPRE_IJMatrixSetObjectType,(hA->ij,HYPRE_PARCSR));
1133:   } else {
1134:     HYPRE_Int hrs,hre,hcs,hce;
1135:     PetscStackCallStandard(HYPRE_IJMatrixGetLocalRange,(hA->ij,&hrs,&hre,&hcs,&hce));
1136:     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);
1137:     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);
1138:   }
1139:   PetscStackCallStandard(HYPRE_IJMatrixInitialize,(hA->ij));

1141:   if (dnz == PETSC_DEFAULT || dnz == PETSC_DECIDE) dnz = 10*bs;
1142:   if (onz == PETSC_DEFAULT || onz == PETSC_DECIDE) onz = 10*bs;

1144:   if (!dnnz) {
1145:     PetscMalloc1(A->rmap->n,&hdnnz);
1146:     for (i=0;i<A->rmap->n;i++) hdnnz[i] = dnz;
1147:   } else {
1148:     hdnnz = (HYPRE_Int*)dnnz;
1149:   }
1150:   MPI_Comm_size(PetscObjectComm((PetscObject)A),&size);
1151:   if (size > 1) {
1152:     if (!onnz) {
1153:       PetscMalloc1(A->rmap->n,&honnz);
1154:       for (i=0;i<A->rmap->n;i++) honnz[i] = onz;
1155:     } else {
1156:       honnz = (HYPRE_Int*)onnz;
1157:     }
1158:     PetscStackCallStandard(HYPRE_IJMatrixSetDiagOffdSizes,(hA->ij,hdnnz,honnz));
1159:   } else {
1160:     honnz = NULL;
1161:     PetscStackCallStandard(HYPRE_IJMatrixSetRowSizes,(hA->ij,hdnnz));
1162:   }
1163:   if (!dnnz) {
1164:     PetscFree(hdnnz);
1165:   }
1166:   if (!onnz && honnz) {
1167:     PetscFree(honnz);
1168:   }
1169:   A->preallocated = PETSC_TRUE;

1171:   /* SetDiagOffdSizes sets hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 0 */
1172:   {
1173:     hypre_AuxParCSRMatrix *aux_matrix;
1174:     aux_matrix = (hypre_AuxParCSRMatrix*)hypre_IJMatrixTranslator(hA->ij);
1175:     hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 1;
1176:   }
1177:   return(0);
1178: }

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

1183:    Collective on Mat

1185:    Input Parameters:
1186: +  A - the matrix
1187: .  dnz  - number of nonzeros per row in DIAGONAL portion of local submatrix
1188:           (same value is used for all local rows)
1189: .  dnnz - array containing the number of nonzeros in the various rows of the
1190:           DIAGONAL portion of the local submatrix (possibly different for each row)
1191:           or NULL (PETSC_NULL_INTEGER in Fortran), if d_nz is used to specify the nonzero structure.
1192:           The size of this array is equal to the number of local rows, i.e 'm'.
1193:           For matrices that will be factored, you must leave room for (and set)
1194:           the diagonal entry even if it is zero.
1195: .  onz  - number of nonzeros per row in the OFF-DIAGONAL portion of local
1196:           submatrix (same value is used for all local rows).
1197: -  onnz - array containing the number of nonzeros in the various rows of the
1198:           OFF-DIAGONAL portion of the local submatrix (possibly different for
1199:           each row) or NULL (PETSC_NULL_INTEGER in Fortran), if o_nz is used to specify the nonzero
1200:           structure. The size of this array is equal to the number
1201:           of local rows, i.e 'm'.

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

1206:    Level: intermediate

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

1210: .seealso: MatCreate(), MatMPIAIJSetPreallocation, MATHYPRE
1211: @*/
1212: PetscErrorCode MatHYPRESetPreallocation(Mat A, PetscInt dnz, const PetscInt dnnz[], PetscInt onz, const PetscInt onnz[])
1213: {

1219:   PetscTryMethod(A,"MatHYPRESetPreallocation_C",(Mat,PetscInt,const PetscInt[],PetscInt,const PetscInt[]),(A,dnz,dnnz,onz,onnz));
1220:   return(0);
1221: }

1223: /*
1224:    MatCreateFromParCSR - Creates a matrix from a hypre_ParCSRMatrix

1226:    Collective

1228:    Input Parameters:
1229: +  vparcsr  - the pointer to the hypre_ParCSRMatrix
1230: .  mtype    - matrix type to be created. Currently MATAIJ, MATIS and MATHYPRE are supported.
1231: -  copymode - PETSc copying options

1233:    Output Parameter:
1234: .  A  - the matrix

1236:    Level: intermediate

1238: .seealso: MatHYPRE, PetscCopyMode
1239: */
1240: PETSC_EXTERN PetscErrorCode MatCreateFromParCSR(hypre_ParCSRMatrix *vparcsr, MatType mtype, PetscCopyMode copymode, Mat* A)
1241: {
1242:   Mat                   T;
1243:   Mat_HYPRE             *hA;
1244:   hypre_ParCSRMatrix    *parcsr;
1245:   MPI_Comm              comm;
1246:   PetscInt              rstart,rend,cstart,cend,M,N;
1247:   PetscBool             isseqaij,ismpiaij,isaij,ishyp,isis;
1248:   PetscErrorCode        ierr;

1251:   parcsr = (hypre_ParCSRMatrix *)vparcsr;
1252:   comm   = hypre_ParCSRMatrixComm(parcsr);
1253:   PetscStrcmp(mtype,MATSEQAIJ,&isseqaij);
1254:   PetscStrcmp(mtype,MATMPIAIJ,&ismpiaij);
1255:   PetscStrcmp(mtype,MATAIJ,&isaij);
1256:   PetscStrcmp(mtype,MATHYPRE,&ishyp);
1257:   PetscStrcmp(mtype,MATIS,&isis);
1258:   isaij  = (PetscBool)(isseqaij || ismpiaij || isaij);
1259:   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);
1260:   if (ishyp && copymode == PETSC_COPY_VALUES) SETERRQ(comm,PETSC_ERR_SUP,"Unsupported copymode PETSC_COPY_VALUES");

1262:   /* access ParCSRMatrix */
1263:   rstart = hypre_ParCSRMatrixFirstRowIndex(parcsr);
1264:   rend   = hypre_ParCSRMatrixLastRowIndex(parcsr);
1265:   cstart = hypre_ParCSRMatrixFirstColDiag(parcsr);
1266:   cend   = hypre_ParCSRMatrixLastColDiag(parcsr);
1267:   M      = hypre_ParCSRMatrixGlobalNumRows(parcsr);
1268:   N      = hypre_ParCSRMatrixGlobalNumCols(parcsr);

1270:   /* fix for empty local rows/columns */
1271:   if (rend < rstart) rend = rstart;
1272:   if (cend < cstart) cend = cstart;

1274:   /* create PETSc matrix with MatHYPRE */
1275:   MatCreate(comm,&T);
1276:   MatSetSizes(T,rend-rstart+1,cend-cstart+1,M,N);
1277:   MatSetType(T,MATHYPRE);
1278:   hA   = (Mat_HYPRE*)(T->data);

1280:   /* create HYPRE_IJMatrix */
1281:   PetscStackCallStandard(HYPRE_IJMatrixCreate,(hA->comm,rstart,rend-1,cstart,cend-1,&hA->ij));

1283:   /* set ParCSR object */
1284:   PetscStackCallStandard(HYPRE_IJMatrixSetObjectType,(hA->ij,HYPRE_PARCSR));
1285:   hypre_IJMatrixObject(hA->ij) = parcsr;
1286:   T->preallocated = PETSC_TRUE;

1288:   /* set assembled flag */
1289:   hypre_IJMatrixAssembleFlag(hA->ij) = 1;
1290:   PetscStackCallStandard(HYPRE_IJMatrixInitialize,(hA->ij));
1291:   if (ishyp) {
1292:     PetscMPIInt myid = 0;

1294:     /* make sure we always have row_starts and col_starts available */
1295:     if (HYPRE_AssumedPartitionCheck()) {
1296:       MPI_Comm_rank(comm,&myid);
1297:     }
1298:     if (!hypre_ParCSRMatrixOwnsColStarts(parcsr)) {
1299:       PetscLayout map;

1301:       MatGetLayouts(T,NULL,&map);
1302:       PetscLayoutSetUp(map);
1303:       hypre_ParCSRMatrixColStarts(parcsr) = (HYPRE_Int*)(map->range + myid);
1304:     }
1305:     if (!hypre_ParCSRMatrixOwnsRowStarts(parcsr)) {
1306:       PetscLayout map;

1308:       MatGetLayouts(T,&map,NULL);
1309:       PetscLayoutSetUp(map);
1310:       hypre_ParCSRMatrixRowStarts(parcsr) = (HYPRE_Int*)(map->range + myid);
1311:     }
1312:     /* prevent from freeing the pointer */
1313:     if (copymode == PETSC_USE_POINTER) hA->inner_free = PETSC_FALSE;
1314:     *A   = T;
1315:     MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY);
1316:     MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY);
1317:   } else if (isaij) {
1318:     if (copymode != PETSC_OWN_POINTER) {
1319:       /* prevent from freeing the pointer */
1320:       hA->inner_free = PETSC_FALSE;
1321:       MatConvert_HYPRE_AIJ(T,MATAIJ,MAT_INITIAL_MATRIX,A);
1322:       MatDestroy(&T);
1323:     } else { /* AIJ return type with PETSC_OWN_POINTER */
1324:       MatConvert_HYPRE_AIJ(T,MATAIJ,MAT_INPLACE_MATRIX,&T);
1325:       *A   = T;
1326:     }
1327:   } else if (isis) {
1328:     MatConvert_HYPRE_IS(T,MATIS,MAT_INITIAL_MATRIX,A);
1329:     if (copymode != PETSC_OWN_POINTER) hA->inner_free = PETSC_FALSE;
1330:     MatDestroy(&T);
1331:   }
1332:   return(0);
1333: }

1335: static PetscErrorCode MatHYPREGetParCSR_HYPRE(Mat A, hypre_ParCSRMatrix **parcsr)
1336: {
1337:   Mat_HYPRE*            hA = (Mat_HYPRE*)A->data;
1338:   HYPRE_Int             type;
1339:   PetscErrorCode        ierr;

1342:   if (!hA->ij) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"HYPRE_IJMatrix not present");
1343:   PetscStackCallStandard(HYPRE_IJMatrixGetObjectType,(hA->ij,&type));
1344:   if (type != HYPRE_PARCSR) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"HYPRE_IJMatrix is not of type HYPRE_PARCSR");
1345:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hA->ij,(void**)parcsr));
1346:   return(0);
1347: }

1349: /*
1350:    MatHYPREGetParCSR - Gets the pointer to the ParCSR matrix

1352:    Not collective

1354:    Input Parameters:
1355: +  A  - the MATHYPRE object

1357:    Output Parameter:
1358: .  parcsr  - the pointer to the hypre_ParCSRMatrix

1360:    Level: intermediate

1362: .seealso: MatHYPRE, PetscCopyMode
1363: */
1364: PetscErrorCode MatHYPREGetParCSR(Mat A, hypre_ParCSRMatrix **parcsr)
1365: {

1371:   PetscUseMethod(A,"MatHYPREGetParCSR_C",(Mat,hypre_ParCSRMatrix**),(A,parcsr));
1372:   return(0);
1373: }

1375: static PetscErrorCode MatMissingDiagonal_HYPRE(Mat A, PetscBool *missing, PetscInt *dd)
1376: {
1377:   hypre_ParCSRMatrix *parcsr;
1378:   hypre_CSRMatrix    *ha;
1379:   PetscInt           rst;
1380:   PetscErrorCode     ierr;

1383:   if (A->rmap->n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not implemented with non-square diagonal blocks");
1384:   MatGetOwnershipRange(A,&rst,NULL);
1385:   MatHYPREGetParCSR_HYPRE(A,&parcsr);
1386:   if (missing) *missing = PETSC_FALSE;
1387:   if (dd) *dd = -1;
1388:   ha = hypre_ParCSRMatrixDiag(parcsr);
1389:   if (ha) {
1390:     PetscInt  size,i;
1391:     HYPRE_Int *ii,*jj;

1393:     size = hypre_CSRMatrixNumRows(ha);
1394:     ii   = hypre_CSRMatrixI(ha);
1395:     jj   = hypre_CSRMatrixJ(ha);
1396:     for (i = 0; i < size; i++) {
1397:       PetscInt  j;
1398:       PetscBool found = PETSC_FALSE;

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

1403:       if (!found) {
1404:         PetscInfo1(A,"Matrix is missing local diagonal entry %D\n",i);
1405:         if (missing) *missing = PETSC_TRUE;
1406:         if (dd) *dd = i+rst;
1407:         return(0);
1408:       }
1409:     }
1410:     if (!size) {
1411:       PetscInfo(A,"Matrix has no diagonal entries therefore is missing diagonal\n");
1412:       if (missing) *missing = PETSC_TRUE;
1413:       if (dd) *dd = rst;
1414:     }
1415:   } else {
1416:     PetscInfo(A,"Matrix has no diagonal entries therefore is missing diagonal\n");
1417:     if (missing) *missing = PETSC_TRUE;
1418:     if (dd) *dd = rst;
1419:   }
1420:   return(0);
1421: }

1423: static PetscErrorCode MatScale_HYPRE(Mat A, PetscScalar s)
1424: {
1425:   hypre_ParCSRMatrix *parcsr;
1426:   hypre_CSRMatrix    *ha;
1427:   PetscErrorCode     ierr;

1430:   MatHYPREGetParCSR_HYPRE(A,&parcsr);
1431:   /* diagonal part */
1432:   ha = hypre_ParCSRMatrixDiag(parcsr);
1433:   if (ha) {
1434:     PetscInt    size,i;
1435:     HYPRE_Int   *ii;
1436:     PetscScalar *a;

1438:     size = hypre_CSRMatrixNumRows(ha);
1439:     a    = hypre_CSRMatrixData(ha);
1440:     ii   = hypre_CSRMatrixI(ha);
1441:     for (i = 0; i < ii[size]; i++) a[i] *= s;
1442:   }
1443:   /* offdiagonal part */
1444:   ha = hypre_ParCSRMatrixOffd(parcsr);
1445:   if (ha) {
1446:     PetscInt    size,i;
1447:     HYPRE_Int   *ii;
1448:     PetscScalar *a;

1450:     size = hypre_CSRMatrixNumRows(ha);
1451:     a    = hypre_CSRMatrixData(ha);
1452:     ii   = hypre_CSRMatrixI(ha);
1453:     for (i = 0; i < ii[size]; i++) a[i] *= s;
1454:   }
1455:   return(0);
1456: }

1458: static PetscErrorCode MatZeroRowsColumns_HYPRE(Mat A, PetscInt numRows, const PetscInt rows[], PetscScalar diag, Vec x, Vec b)
1459: {
1460:   hypre_ParCSRMatrix *parcsr;
1461:   HYPRE_Int          *lrows;
1462:   PetscInt           rst,ren,i;
1463:   PetscErrorCode     ierr;

1466:   if (x || b) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"To be implemented");
1467:   MatHYPREGetParCSR_HYPRE(A,&parcsr);
1468:   PetscMalloc1(numRows,&lrows);
1469:   MatGetOwnershipRange(A,&rst,&ren);
1470:   for (i=0;i<numRows;i++) {
1471:     if (rows[i] < rst || rows[i] >= ren)
1472:       SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Non-local rows not yet supported");
1473:     lrows[i] = rows[i] - rst;
1474:   }
1475:   PetscStackCallStandard(hypre_ParCSRMatrixEliminateRowsCols,(parcsr,numRows,lrows));
1476:   PetscFree(lrows);
1477:   return(0);
1478: }

1480: static PetscErrorCode MatZeroEntries_HYPRE_CSRMatrix(hypre_CSRMatrix *ha)
1481: {
1482:   PetscErrorCode      ierr;

1485:   if (ha) {
1486:     HYPRE_Int     *ii, size;
1487:     HYPRE_Complex *a;

1489:     size = hypre_CSRMatrixNumRows(ha);
1490:     a    = hypre_CSRMatrixData(ha);
1491:     ii   = hypre_CSRMatrixI(ha);

1493:     if (a) { PetscMemzero(a,(ii[size])*sizeof(HYPRE_Complex)); }
1494:   }
1495:   return(0);
1496: }

1498: PetscErrorCode MatZeroEntries_HYPRE(Mat A)
1499: {
1500:   hypre_ParCSRMatrix  *parcsr;
1501:   PetscErrorCode      ierr;

1504:   MatHYPREGetParCSR_HYPRE(A,&parcsr);
1505:   /* diagonal part */
1506:   MatZeroEntries_HYPRE_CSRMatrix(hypre_ParCSRMatrixDiag(parcsr));
1507:   /* off-diagonal part */
1508:   MatZeroEntries_HYPRE_CSRMatrix(hypre_ParCSRMatrixOffd(parcsr));
1509:   return(0);
1510: }


1513: static PetscErrorCode MatZeroRows_HYPRE_CSRMatrix(hypre_CSRMatrix *hA,PetscInt N,const PetscInt rows[],PetscScalar diag)
1514: {
1515:   PetscInt        ii, jj, ibeg, iend, irow;
1516:   PetscInt        *i, *j;
1517:   PetscScalar     *a;


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

1523:   i = (PetscInt*) hypre_CSRMatrixI(hA);
1524:   j = (PetscInt*) hypre_CSRMatrixJ(hA);
1525:   a = hypre_CSRMatrixData(hA);

1527:   for (ii = 0; ii < N; ii++) {
1528:     irow = rows[ii];
1529:     ibeg = i[irow];
1530:     iend = i[irow+1];
1531:     for (jj = ibeg; jj < iend; jj++)
1532:       if (j[jj] == irow) a[jj] = diag;
1533:       else a[jj] = 0.0;
1534:    }

1536:    return(0);
1537: }

1539: PetscErrorCode MatZeroRows_HYPRE(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
1540: {
1541:   hypre_ParCSRMatrix  *parcsr;
1542:   PetscInt            *lrows,len;
1543:   PetscErrorCode      ierr;

1546:   if (x || b) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Does not support to modify the solution and the right hand size \n");
1547:   /* retrieve the internal matrix */
1548:   MatHYPREGetParCSR_HYPRE(A,&parcsr);
1549:   /* get locally owned rows */
1550:   MatZeroRowsMapLocal_Private(A,N,rows,&len,&lrows);
1551:   /* zero diagonal part */
1552:   MatZeroRows_HYPRE_CSRMatrix(hypre_ParCSRMatrixDiag(parcsr),len,lrows,diag);
1553:   /* zero off-diagonal part */
1554:   MatZeroRows_HYPRE_CSRMatrix(hypre_ParCSRMatrixOffd(parcsr),len,lrows,0.0);

1556:   PetscFree(lrows);
1557:   return(0);
1558: }


1561: PetscErrorCode MatAssemblyBegin_HYPRE(Mat mat,MatAssemblyType mode)
1562: {

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

1568:   MatStashScatterBegin_Private(mat,&mat->stash,mat->rmap->range);
1569:   return(0);
1570: }

1572: PetscErrorCode MatGetRow_HYPRE(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v)
1573: {
1574:   hypre_ParCSRMatrix  *parcsr;
1575:   PetscErrorCode      ierr;

1578:   /* retrieve the internal matrix */
1579:   MatHYPREGetParCSR_HYPRE(A,&parcsr);
1580:   /* call HYPRE API */
1581:   PetscStackCallStandard(HYPRE_ParCSRMatrixGetRow,(parcsr,row,(HYPRE_Int*)nz,(HYPRE_Int**)idx,v));
1582:   return(0);
1583: }

1585: PetscErrorCode MatRestoreRow_HYPRE(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v)
1586: {
1587:   hypre_ParCSRMatrix  *parcsr;
1588:   PetscErrorCode      ierr;

1591:   /* retrieve the internal matrix */
1592:   MatHYPREGetParCSR_HYPRE(A,&parcsr);
1593:   /* call HYPRE API */
1594:   PetscStackCallStandard(HYPRE_ParCSRMatrixRestoreRow,(parcsr,row,(HYPRE_Int*)nz,(HYPRE_Int**)idx,v));
1595:   return(0);
1596: }



1600: PetscErrorCode MatGetValues_HYPRE(Mat A,PetscInt m,const PetscInt idxm[],PetscInt n,const PetscInt idxn[],PetscScalar v[])
1601: {
1602:   HYPRE_IJMatrix     *hIJ = (HYPRE_IJMatrix*)A->data;
1603:   PetscErrorCode      ierr;
1604:   PetscInt            i;
1606:   if (!m || !n) return(0);

1608:   /* Ignore negative row indices
1609:    * And negative column indices should be automatically ignored in hypre
1610:    * */
1611:   for (i=0; i<m; i++)
1612:     if (idxm[i] >= 0) PetscStackCallStandard(HYPRE_IJMatrixGetValues,(*hIJ,1,(HYPRE_Int*)&n,(HYPRE_Int*)&idxm[i],(HYPRE_Int*)idxn,&v[i*n]));

1614:   return(0);
1615: }


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

1622:    Level: intermediate

1624: .seealso: MatCreate()
1625: M*/

1627: PETSC_EXTERN PetscErrorCode MatCreate_HYPRE(Mat B)
1628: {
1629:   Mat_HYPRE      *hB;

1633:   PetscNewLog(B,&hB);
1634:   hB->inner_free = PETSC_TRUE;
1635:   hB->available  = PETSC_TRUE;
1636:   hB->size       = 0;
1637:   hB->array      = NULL;

1639:   B->data       = (void*)hB;
1640:   B->rmap->bs   = 1;
1641:   B->assembled  = PETSC_FALSE;

1643:   PetscMemzero(B->ops,sizeof(struct _MatOps));
1644:   B->ops->mult            = MatMult_HYPRE;
1645:   B->ops->multtranspose   = MatMultTranspose_HYPRE;
1646:   B->ops->setup           = MatSetUp_HYPRE;
1647:   B->ops->destroy         = MatDestroy_HYPRE;

1649:   /* build cache for off array entries formed */
1650:   MatStashCreate_Private(PetscObjectComm((PetscObject)B),1,&B->stash);
1651:   B->ops->assemblyend     = MatAssemblyEnd_HYPRE;
1652:   B->ops->assemblybegin   = MatAssemblyBegin_HYPRE;

1654:   B->ops->ptap            = MatPtAP_HYPRE_HYPRE;
1655:   B->ops->matmult         = MatMatMult_HYPRE_HYPRE;
1656:   B->ops->setvalues       = MatSetValues_HYPRE;
1657:   B->ops->missingdiagonal = MatMissingDiagonal_HYPRE;
1658:   B->ops->scale           = MatScale_HYPRE;
1659:   B->ops->zerorowscolumns = MatZeroRowsColumns_HYPRE;
1660:   B->ops->zeroentries     = MatZeroEntries_HYPRE;
1661:   B->ops->zerorows        = MatZeroRows_HYPRE;
1662:   B->ops->getrow          = MatGetRow_HYPRE;
1663:   B->ops->restorerow      = MatRestoreRow_HYPRE;
1664:   B->ops->getvalues       = MatGetValues_HYPRE;

1666:   MPI_Comm_dup(PetscObjectComm((PetscObject)B),&hB->comm);
1667:   PetscObjectChangeTypeName((PetscObject)B,MATHYPRE);
1668:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_hypre_aij_C",MatConvert_HYPRE_AIJ);
1669:   PetscObjectComposeFunction((PetscObject)B,"MatConvert_hypre_is_C",MatConvert_HYPRE_IS);
1670:   PetscObjectComposeFunction((PetscObject)B,"MatPtAP_seqaij_hypre_C",MatPtAP_AIJ_HYPRE);
1671:   PetscObjectComposeFunction((PetscObject)B,"MatPtAP_mpiaij_hypre_C",MatPtAP_AIJ_HYPRE);
1672:   PetscObjectComposeFunction((PetscObject)B,"MatHYPRESetPreallocation_C",MatHYPRESetPreallocation_HYPRE);
1673:   PetscObjectComposeFunction((PetscObject)B,"MatHYPREGetParCSR_C",MatHYPREGetParCSR_HYPRE);
1674:   PetscObjectComposeFunction((PetscObject)B,"MatPtAP_is_hypre_C",MatPtAP_IS_XAIJ);
1675:   return(0);
1676: }

1678: static PetscErrorCode hypre_array_destroy(void *ptr)
1679: {
1681:    hypre_TFree(ptr,HYPRE_MEMORY_HOST);
1682:    return(0);
1683: }