Actual source code: mhyp.c

petsc-3.5.2 2014-09-08
Report Typos and Errors
  2: /*
  3:     Creates hypre ijmatrix from PETSc matrix
  4: */
  5: #include <petscsys.h>
  6: #include <petsc-private/matimpl.h>
  7: #include <petscdmda.h>                /*I "petscdmda.h" I*/
  8: #include <../src/dm/impls/da/hypre/mhyp.h>

 12: PetscErrorCode MatHYPRE_IJMatrixPreallocate(Mat A_d, Mat A_o,HYPRE_IJMatrix ij)
 13: {
 15:   PetscInt       i,n_d,n_o;
 16:   const PetscInt *ia_d,*ia_o;
 17:   PetscBool      done_d=PETSC_FALSE,done_o=PETSC_FALSE;
 18:   PetscInt       *nnz_d=NULL,*nnz_o=NULL;

 21:   if (A_d) { /* determine number of nonzero entries in local diagonal part */
 22:     MatGetRowIJ(A_d,0,PETSC_FALSE,PETSC_FALSE,&n_d,&ia_d,NULL,&done_d);
 23:     if (done_d) {
 24:       PetscMalloc1(n_d,&nnz_d);
 25:       for (i=0; i<n_d; i++) {
 26:         nnz_d[i] = ia_d[i+1] - ia_d[i];
 27:       }
 28:     }
 29:     MatRestoreRowIJ(A_d,0,PETSC_FALSE,PETSC_FALSE,NULL,&ia_d,NULL,&done_d);
 30:   }
 31:   if (A_o) { /* determine number of nonzero entries in local off-diagonal part */
 32:     MatGetRowIJ(A_o,0,PETSC_FALSE,PETSC_FALSE,&n_o,&ia_o,NULL,&done_o);
 33:     if (done_o) {
 34:       PetscMalloc1(n_o,&nnz_o);
 35:       for (i=0; i<n_o; i++) {
 36:         nnz_o[i] = ia_o[i+1] - ia_o[i];
 37:       }
 38:     }
 39:     MatRestoreRowIJ(A_o,0,PETSC_FALSE,PETSC_FALSE,&n_o,&ia_o,NULL,&done_o);
 40:   }
 41:   if (done_d) {    /* set number of nonzeros in HYPRE IJ matrix */
 42:     if (!done_o) { /* only diagonal part */
 43:       PetscMalloc1(n_d,&nnz_o);
 44:       for (i=0; i<n_d; i++) {
 45:         nnz_o[i] = 0;
 46:       }
 47:     }
 48:     PetscStackCallStandard(HYPRE_IJMatrixSetDiagOffdSizes,(ij,(HYPRE_Int *)nnz_d,(HYPRE_Int *)nnz_o));
 49:     PetscFree(nnz_d);
 50:     PetscFree(nnz_o);
 51:   }
 52:   return(0);
 53: }

 57: PetscErrorCode MatHYPRE_IJMatrixCreate(Mat A,HYPRE_IJMatrix *ij)
 58: {
 60:   PetscInt       rstart,rend,cstart,cend;

 66:   MatSetUp(A);
 67:   rstart = A->rmap->rstart;
 68:   rend   = A->rmap->rend;
 69:   cstart = A->cmap->rstart;
 70:   cend   = A->cmap->rend;
 71:   PetscStackCallStandard(HYPRE_IJMatrixCreate,(PetscObjectComm((PetscObject)A),rstart,rend-1,cstart,cend-1,ij));
 72:   PetscStackCallStandard(HYPRE_IJMatrixSetObjectType,(*ij,HYPRE_PARCSR));
 73:   {
 74:     PetscBool      same;
 75:     Mat            A_d,A_o;
 76:     const PetscInt *colmap;
 77:     PetscObjectTypeCompare((PetscObject)A,MATMPIAIJ,&same);
 78:     if (same) {
 79:       MatMPIAIJGetSeqAIJ(A,&A_d,&A_o,&colmap);
 80:       MatHYPRE_IJMatrixPreallocate(A_d,A_o,*ij);
 81:       return(0);
 82:     }
 83:     PetscObjectTypeCompare((PetscObject)A,MATMPIBAIJ,&same);
 84:     if (same) {
 85:       MatMPIBAIJGetSeqBAIJ(A,&A_d,&A_o,&colmap);
 86:       MatHYPRE_IJMatrixPreallocate(A_d,A_o,*ij);
 87:       return(0);
 88:     }
 89:     PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&same);
 90:     if (same) {
 91:       MatHYPRE_IJMatrixPreallocate(A,NULL,*ij);
 92:       return(0);
 93:     }
 94:     PetscObjectTypeCompare((PetscObject)A,MATSEQBAIJ,&same);
 95:     if (same) {
 96:       MatHYPRE_IJMatrixPreallocate(A,NULL,*ij);
 97:       return(0);
 98:     }
 99:   }
100:   return(0);
101: }

103: extern PetscErrorCode MatHYPRE_IJMatrixFastCopy_MPIAIJ(Mat,HYPRE_IJMatrix);
104: extern PetscErrorCode MatHYPRE_IJMatrixFastCopy_SeqAIJ(Mat,HYPRE_IJMatrix);
105: /*
106:     Copies the data over (column indices, numerical values) to hypre matrix
107: */

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

120:   PetscObjectTypeCompare((PetscObject)A,MATMPIAIJ,&flg);
121:   if (flg) {
122:     MatHYPRE_IJMatrixFastCopy_MPIAIJ(A,ij);
123:     return(0);
124:   }
125:   PetscObjectTypeCompare((PetscObject)A,MATSEQAIJ,&flg);
126:   if (flg) {
127:     MatHYPRE_IJMatrixFastCopy_SeqAIJ(A,ij);
128:     return(0);
129:   }

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

144: /*
145:     This copies the CSR format directly from the PETSc data structure to the hypre
146:     data structure without calls to MatGetRow() or hypre's set values.

148: */
149: #include <_hypre_IJ_mv.h>
150: #include <HYPRE_IJ_mv.h>
151: #include <../src/mat/impls/aij/mpi/mpiaij.h>

155: PetscErrorCode MatHYPRE_IJMatrixFastCopy_SeqAIJ(Mat A,HYPRE_IJMatrix ij)
156: {
158:   Mat_SeqAIJ     *pdiag = (Mat_SeqAIJ*)A->data;

160:   hypre_ParCSRMatrix    *par_matrix;
161:   hypre_AuxParCSRMatrix *aux_matrix;
162:   hypre_CSRMatrix       *hdiag;


169:   PetscLogEventBegin(MAT_Convert,A,0,0,0);
170:   PetscStackCallStandard(HYPRE_IJMatrixInitialize,(ij));
171:   par_matrix = (hypre_ParCSRMatrix*)hypre_IJMatrixObject(ij);
172:   aux_matrix = (hypre_AuxParCSRMatrix*)hypre_IJMatrixTranslator(ij);
173:   hdiag      = hypre_ParCSRMatrixDiag(par_matrix);

175:   /*
176:        this is the Hack part where we monkey directly with the hypre datastructures
177:   */
178:   PetscMemcpy(hdiag->i,pdiag->i,(A->rmap->n + 1)*sizeof(PetscInt));
179:   PetscMemcpy(hdiag->j,pdiag->j,pdiag->nz*sizeof(PetscInt));
180:   PetscMemcpy(hdiag->data,pdiag->a,pdiag->nz*sizeof(PetscScalar));

182:   hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 0;
183:   PetscStackCallStandard(HYPRE_IJMatrixAssemble,(ij));
184:   PetscLogEventEnd(MAT_Convert,A,0,0,0);
185:   return(0);
186: }

190: PetscErrorCode MatHYPRE_IJMatrixFastCopy_MPIAIJ(Mat A,HYPRE_IJMatrix ij)
191: {
193:   Mat_MPIAIJ     *pA = (Mat_MPIAIJ*)A->data;
194:   Mat_SeqAIJ     *pdiag,*poffd;
195:   PetscInt       i,*garray = pA->garray,*jj,cstart,*pjj;

197:   hypre_ParCSRMatrix    *par_matrix;
198:   hypre_AuxParCSRMatrix *aux_matrix;
199:   hypre_CSRMatrix       *hdiag,*hoffd;

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

210:   PetscLogEventBegin(MAT_Convert,A,0,0,0);
211:   PetscStackCallStandard(HYPRE_IJMatrixInitialize,(ij));
212:   par_matrix = (hypre_ParCSRMatrix*)hypre_IJMatrixObject(ij);
213:   aux_matrix = (hypre_AuxParCSRMatrix*)hypre_IJMatrixTranslator(ij);
214:   hdiag      = hypre_ParCSRMatrixDiag(par_matrix);
215:   hoffd      = hypre_ParCSRMatrixOffd(par_matrix);

217:   /*
218:        this is the Hack part where we monkey directly with the hypre datastructures
219:   */
220:   PetscMemcpy(hdiag->i,pdiag->i,(pA->A->rmap->n + 1)*sizeof(PetscInt));
221:   /* need to shift the diag column indices (hdiag->j) back to global numbering since hypre is expecting this */
222:   jj  = (PetscInt*)hdiag->j;
223:   pjj = (PetscInt*)pdiag->j;
224:   for (i=0; i<pdiag->nz; i++) jj[i] = cstart + pjj[i];
225:   PetscMemcpy(hdiag->data,pdiag->a,pdiag->nz*sizeof(PetscScalar));
226:   PetscMemcpy(hoffd->i,poffd->i,(pA->A->rmap->n + 1)*sizeof(PetscInt));
227:   /* need to move the offd column indices (hoffd->j) back to global numbering since hypre is expecting this
228:      If we hacked a hypre a bit more we might be able to avoid this step */
229:   jj  = (PetscInt*) hoffd->j;
230:   pjj = (PetscInt*) poffd->j;
231:   for (i=0; i<poffd->nz; i++) jj[i] = garray[pjj[i]];
232:   PetscMemcpy(hoffd->data,poffd->a,poffd->nz*sizeof(PetscScalar));

234:   hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 0;
235:   PetscStackCallStandard(HYPRE_IJMatrixAssemble,(ij));
236:   PetscLogEventEnd(MAT_Convert,A,0,0,0);
237:   return(0);
238: }

240: /*
241:     Does NOT copy the data over, instead uses DIRECTLY the pointers from the PETSc MPIAIJ format

243:     This is UNFINISHED and does NOT work! The problem is that hypre puts the diagonal entry first
244:     which will corrupt the PETSc data structure if we did this. Need a work around to this problem.
245: */
246: #include <_hypre_IJ_mv.h>
247: #include <HYPRE_IJ_mv.h>

251: PetscErrorCode MatHYPRE_IJMatrixLink(Mat A,HYPRE_IJMatrix *ij)
252: {
253:   PetscErrorCode        ierr;
254:   PetscInt              rstart,rend,cstart,cend;
255:   PetscBool             flg;
256:   hypre_AuxParCSRMatrix *aux_matrix;

262:   PetscObjectTypeCompare((PetscObject)A,MATMPIAIJ,&flg);
263:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Can only use with PETSc MPIAIJ matrices");
264:   MatSetUp(A);

266:   PetscLogEventBegin(MAT_Convert,A,0,0,0);
267:   rstart = A->rmap->rstart;
268:   rend   = A->rmap->rend;
269:   cstart = A->cmap->rstart;
270:   cend   = A->cmap->rend;
271:   PetscStackCallStandard(HYPRE_IJMatrixCreate,(PetscObjectComm((PetscObject)A),rstart,rend-1,cstart,cend-1,ij));
272:   PetscStackCallStandard(HYPRE_IJMatrixSetObjectType,(*ij,HYPRE_PARCSR));

274:   PetscStackCallStandard(HYPRE_IJMatrixInitialize,(*ij));
275:   PetscStackCall("hypre_IJMatrixTranslator",aux_matrix = (hypre_AuxParCSRMatrix*)hypre_IJMatrixTranslator(*ij));

277:   hypre_AuxParCSRMatrixNeedAux(aux_matrix) = 0;

279:   /* this is the Hack part where we monkey directly with the hypre datastructures */

281:   PetscStackCallStandard(HYPRE_IJMatrixAssemble,(*ij));
282:   PetscLogEventEnd(MAT_Convert,A,0,0,0);
283:   return(0);
284: }

286: /* -----------------------------------------------------------------------------------------------------------------*/

288: /*MC
289:    MATHYPRESTRUCT - MATHYPRESTRUCT = "hyprestruct" - A matrix type to be used for parallel sparse matrices
290:           based on the hypre HYPRE_StructMatrix.

292:    Level: intermediate

294:    Notes: Unlike the more general support for blocks in hypre this allows only one block per process and requires the block
295:           be defined by a DMDA.

297:           The matrix needs a DMDA associated with it by either a call to MatSetupDM() or if the matrix is obtained from DMCreateMatrix()

299: .seealso: MatCreate(), PCPFMG, MatSetupDM(), DMCreateMatrix()
300: M*/


305: PetscErrorCode  MatSetValuesLocal_HYPREStruct_3d(Mat mat,PetscInt nrow,const PetscInt irow[],PetscInt ncol,const PetscInt icol[],const PetscScalar y[],InsertMode addv)
306: {
307:   PetscErrorCode    ierr;
308:   PetscInt          i,j,stencil,index[3],row,entries[7];
309:   const PetscScalar *values = y;
310:   Mat_HYPREStruct   *ex     = (Mat_HYPREStruct*) mat->data;

313:   for (i=0; i<nrow; i++) {
314:     for (j=0; j<ncol; j++) {
315:       stencil = icol[j] - irow[i];
316:       if (!stencil) {
317:         entries[j] = 3;
318:       } else if (stencil == -1) {
319:         entries[j] = 2;
320:       } else if (stencil == 1) {
321:         entries[j] = 4;
322:       } else if (stencil == -ex->gnx) {
323:         entries[j] = 1;
324:       } else if (stencil == ex->gnx) {
325:         entries[j] = 5;
326:       } else if (stencil == -ex->gnxgny) {
327:         entries[j] = 0;
328:       } else if (stencil == ex->gnxgny) {
329:         entries[j] = 6;
330:       } else SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Local row %D local column %D have bad stencil %D",irow[i],icol[j],stencil);
331:     }
332:     row      = ex->gindices[irow[i]] - ex->rstart;
333:     index[0] = ex->xs + (row % ex->nx);
334:     index[1] = ex->ys + ((row/ex->nx) % ex->ny);
335:     index[2] = ex->zs + (row/(ex->nxny));
336:     if (addv == ADD_VALUES) {
337:       PetscStackCallStandard(HYPRE_StructMatrixAddToValues,(ex->hmat,(HYPRE_Int *)index,ncol,(HYPRE_Int *)entries,(PetscScalar*)values));
338:     } else {
339:       PetscStackCallStandard(HYPRE_StructMatrixSetValues,(ex->hmat,(HYPRE_Int *)index,ncol,(HYPRE_Int *)entries,(PetscScalar*)values));
340:     }
341:     values += ncol;
342:   }
343:   return(0);
344: }

348: PetscErrorCode  MatZeroRowsLocal_HYPREStruct_3d(Mat mat,PetscInt nrow,const PetscInt irow[],PetscScalar d,Vec x,Vec b)
349: {
350:   PetscErrorCode  ierr;
351:   PetscInt        i,index[3],row,entries[7] = {0,1,2,3,4,5,6};
352:   PetscScalar     values[7];
353:   Mat_HYPREStruct *ex = (Mat_HYPREStruct*) mat->data;

356:   if (x && b) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"No support");
357:   PetscMemzero(values,7*sizeof(PetscScalar));
358:   values[3] = d;
359:   for (i=0; i<nrow; i++) {
360:     row      = ex->gindices[irow[i]] - ex->rstart;
361:     index[0] = ex->xs + (row % ex->nx);
362:     index[1] = ex->ys + ((row/ex->nx) % ex->ny);
363:     index[2] = ex->zs + (row/(ex->nxny));
364:     PetscStackCallStandard(HYPRE_StructMatrixSetValues,(ex->hmat,(HYPRE_Int *)index,7,(HYPRE_Int *)entries,values));
365:   }
366:   PetscStackCallStandard(HYPRE_StructMatrixAssemble,(ex->hmat));
367:   return(0);
368: }

372: PetscErrorCode MatZeroEntries_HYPREStruct_3d(Mat mat)
373: {
374:   PetscErrorCode  ierr;
375:   PetscInt        indices[7] = {0,1,2,3,4,5,6};
376:   Mat_HYPREStruct *ex        = (Mat_HYPREStruct*) mat->data;

379:   /* hypre has no public interface to do this */
380:   PetscStackCallStandard(hypre_StructMatrixClearBoxValues,(ex->hmat,&ex->hbox,7,(HYPRE_Int *)indices,0,1));
381:   PetscStackCallStandard(HYPRE_StructMatrixAssemble,(ex->hmat));
382:   return(0);
383: }

387: static PetscErrorCode  MatSetupDM_HYPREStruct(Mat mat,DM da)
388: {
389:   PetscErrorCode         ierr;
390:   Mat_HYPREStruct        *ex = (Mat_HYPREStruct*) mat->data;
391:   PetscInt               dim,dof,sw[3],nx,ny,nz,ilower[3],iupper[3],ssize,i;
392:   DMBoundaryType         px,py,pz;
393:   DMDAStencilType        st;
394:   ISLocalToGlobalMapping ltog;

397:   ex->da = da;
398:   PetscObjectReference((PetscObject)da);

400:   DMDAGetInfo(ex->da,&dim,0,0,0,0,0,0,&dof,&sw[0],&px,&py,&pz,&st);
401:   DMDAGetCorners(ex->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);
402:   iupper[0] += ilower[0] - 1;
403:   iupper[1] += ilower[1] - 1;
404:   iupper[2] += ilower[2] - 1;

406:   /* the hypre_Box is used to zero out the matrix entries in MatZeroValues() */
407:   ex->hbox.imin[0] = ilower[0];
408:   ex->hbox.imin[1] = ilower[1];
409:   ex->hbox.imin[2] = ilower[2];
410:   ex->hbox.imax[0] = iupper[0];
411:   ex->hbox.imax[1] = iupper[1];
412:   ex->hbox.imax[2] = iupper[2];

414:   /* create the hypre grid object and set its information */
415:   if (dof > 1) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Currently only support for scalar problems");
416:   if (px || py || pz) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Ask us to add periodic support by calling HYPRE_StructGridSetPeriodic()");
417:   PetscStackCallStandard(HYPRE_StructGridCreate,(ex->hcomm,dim,&ex->hgrid));
418:   PetscStackCallStandard(HYPRE_StructGridSetExtents,(ex->hgrid,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper));
419:   PetscStackCallStandard(HYPRE_StructGridAssemble,(ex->hgrid));

421:   sw[1] = sw[0];
422:   sw[2] = sw[1];
423:   PetscStackCallStandard(HYPRE_StructGridSetNumGhost,(ex->hgrid,(HYPRE_Int *)sw));

425:   /* create the hypre stencil object and set its information */
426:   if (sw[0] > 1) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Ask us to add support for wider stencils");
427:   if (st == DMDA_STENCIL_BOX) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Ask us to add support for box stencils");
428:   if (dim == 1) {
429:     PetscInt offsets[3][1] = {{-1},{0},{1}};
430:     ssize = 3;
431:     PetscStackCallStandard(HYPRE_StructStencilCreate,(dim,ssize,&ex->hstencil));
432:     for (i=0; i<ssize; i++) {
433:       PetscStackCallStandard(HYPRE_StructStencilSetElement,(ex->hstencil,i,(HYPRE_Int *)offsets[i]));
434:     }
435:   } else if (dim == 2) {
436:     PetscInt offsets[5][2] = {{0,-1},{-1,0},{0,0},{1,0},{0,1}};
437:     ssize = 5;
438:     PetscStackCallStandard(HYPRE_StructStencilCreate,(dim,ssize,&ex->hstencil));
439:     for (i=0; i<ssize; i++) {
440:       PetscStackCallStandard(HYPRE_StructStencilSetElement,(ex->hstencil,i,(HYPRE_Int *)offsets[i]));
441:     }
442:   } else if (dim == 3) {
443:     PetscInt offsets[7][3] = {{0,0,-1},{0,-1,0},{-1,0,0},{0,0,0},{1,0,0},{0,1,0},{0,0,1}};
444:     ssize = 7;
445:     PetscStackCallStandard(HYPRE_StructStencilCreate,(dim,ssize,&ex->hstencil));
446:     for (i=0; i<ssize; i++) {
447:       PetscStackCallStandard(HYPRE_StructStencilSetElement,(ex->hstencil,i,(HYPRE_Int *)offsets[i]));
448:     }
449:   }

451:   /* create the HYPRE vector for rhs and solution */
452:   PetscStackCallStandard(HYPRE_StructVectorCreate,(ex->hcomm,ex->hgrid,&ex->hb));
453:   PetscStackCallStandard(HYPRE_StructVectorCreate,(ex->hcomm,ex->hgrid,&ex->hx));
454:   PetscStackCallStandard(HYPRE_StructVectorInitialize,(ex->hb));
455:   PetscStackCallStandard(HYPRE_StructVectorInitialize,(ex->hx));
456:   PetscStackCallStandard(HYPRE_StructVectorAssemble,(ex->hb));
457:   PetscStackCallStandard(HYPRE_StructVectorAssemble,(ex->hx));

459:   /* create the hypre matrix object and set its information */
460:   PetscStackCallStandard(HYPRE_StructMatrixCreate,(ex->hcomm,ex->hgrid,ex->hstencil,&ex->hmat));
461:   PetscStackCallStandard(HYPRE_StructGridDestroy,(ex->hgrid));
462:   PetscStackCallStandard(HYPRE_StructStencilDestroy,(ex->hstencil));
463:   if (ex->needsinitialization) {
464:     PetscStackCallStandard(HYPRE_StructMatrixInitialize,(ex->hmat));
465:     ex->needsinitialization = PETSC_FALSE;
466:   }

468:   /* set the global and local sizes of the matrix */
469:   DMDAGetCorners(da,0,0,0,&nx,&ny,&nz);
470:   MatSetSizes(mat,dof*nx*ny*nz,dof*nx*ny*nz,PETSC_DECIDE,PETSC_DECIDE);
471:   PetscLayoutSetBlockSize(mat->rmap,1);
472:   PetscLayoutSetBlockSize(mat->cmap,1);
473:   PetscLayoutSetUp(mat->rmap);
474:   PetscLayoutSetUp(mat->cmap);

476:   if (dim == 3) {
477:     mat->ops->setvalueslocal = MatSetValuesLocal_HYPREStruct_3d;
478:     mat->ops->zerorowslocal  = MatZeroRowsLocal_HYPREStruct_3d;
479:     mat->ops->zeroentries    = MatZeroEntries_HYPREStruct_3d;

481:     MatZeroEntries_HYPREStruct_3d(mat);
482:   } else SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Only support for 3d DMDA currently");

484:   /* get values that will be used repeatedly in MatSetValuesLocal() and MatZeroRowsLocal() repeatedly */
485:   MatGetOwnershipRange(mat,&ex->rstart,NULL);
486:   DMGetLocalToGlobalMapping(ex->da,&ltog);
487:   ISLocalToGlobalMappingGetIndices(ltog, (const PetscInt **) &ex->gindices);
488:   DMDAGetGhostCorners(ex->da,0,0,0,&ex->gnx,&ex->gnxgny,0);
489:   ex->gnxgny *= ex->gnx;
490:   DMDAGetCorners(ex->da,&ex->xs,&ex->ys,&ex->zs,&ex->nx,&ex->ny,0);
491:   ex->nxny    = ex->nx*ex->ny;
492:   return(0);
493: }

497: PetscErrorCode MatMult_HYPREStruct(Mat A,Vec x,Vec y)
498: {
499:   PetscErrorCode  ierr;
500:   PetscScalar     *xx,*yy;
501:   PetscInt        ilower[3],iupper[3];
502:   Mat_HYPREStruct *mx = (Mat_HYPREStruct*)(A->data);

505:   DMDAGetCorners(mx->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);

507:   iupper[0] += ilower[0] - 1;
508:   iupper[1] += ilower[1] - 1;
509:   iupper[2] += ilower[2] - 1;

511:   /* copy x values over to hypre */
512:   PetscStackCallStandard(HYPRE_StructVectorSetConstantValues,(mx->hb,0.0));
513:   VecGetArray(x,&xx);
514:   PetscStackCallStandard(HYPRE_StructVectorSetBoxValues,(mx->hb,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,xx));
515:   VecRestoreArray(x,&xx);
516:   PetscStackCallStandard(HYPRE_StructVectorAssemble,(mx->hb));
517:   PetscStackCallStandard(HYPRE_StructMatrixMatvec,(1.0,mx->hmat,mx->hb,0.0,mx->hx));

519:   /* copy solution values back to PETSc */
520:   VecGetArray(y,&yy);
521:   PetscStackCallStandard(HYPRE_StructVectorGetBoxValues,(mx->hx,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,yy));
522:   VecRestoreArray(y,&yy);
523:   return(0);
524: }

528: PetscErrorCode MatAssemblyEnd_HYPREStruct(Mat mat,MatAssemblyType mode)
529: {
530:   Mat_HYPREStruct *ex = (Mat_HYPREStruct*) mat->data;
531:   PetscErrorCode  ierr;

534:   PetscStackCallStandard(HYPRE_StructMatrixAssemble,(ex->hmat));
535:   /* PetscStackCallStandard(HYPRE_StructMatrixPrint,("dummy",ex->hmat,0)); */
536:   return(0);
537: }

541: PetscErrorCode MatZeroEntries_HYPREStruct(Mat mat)
542: {
544:   /* before the DMDA is set to the matrix the zero doesn't need to do anything */
545:   return(0);
546: }


551: PetscErrorCode MatDestroy_HYPREStruct(Mat mat)
552: {
553:   Mat_HYPREStruct *ex = (Mat_HYPREStruct*) mat->data;
554:   PetscErrorCode  ierr;

557:   PetscStackCallStandard(HYPRE_StructMatrixDestroy,(ex->hmat));
558:   PetscStackCallStandard(HYPRE_StructVectorDestroy,(ex->hx));
559:   PetscStackCallStandard(HYPRE_StructVectorDestroy,(ex->hb));
560:   return(0);
561: }


566: PETSC_EXTERN PetscErrorCode MatCreate_HYPREStruct(Mat B)
567: {
568:   Mat_HYPREStruct *ex;
569:   PetscErrorCode  ierr;

572:   PetscNewLog(B,&ex);
573:   B->data      = (void*)ex;
574:   B->rmap->bs  = 1;
575:   B->assembled = PETSC_FALSE;

577:   B->insertmode = NOT_SET_VALUES;

579:   B->ops->assemblyend = MatAssemblyEnd_HYPREStruct;
580:   B->ops->mult        = MatMult_HYPREStruct;
581:   B->ops->zeroentries = MatZeroEntries_HYPREStruct;
582:   B->ops->destroy     = MatDestroy_HYPREStruct;

584:   ex->needsinitialization = PETSC_TRUE;

586:   MPI_Comm_dup(PetscObjectComm((PetscObject)B),&(ex->hcomm));
587:   PetscObjectComposeFunction((PetscObject)B,"MatSetupDM_C",MatSetupDM_HYPREStruct);
588:   PetscObjectChangeTypeName((PetscObject)B,MATHYPRESTRUCT);
589:   return(0);
590: }

592: /*MC
593:    MATHYPRESSTRUCT - MATHYPRESSTRUCT = "hypresstruct" - A matrix type to be used for parallel sparse matrices
594:           based on the hypre HYPRE_SStructMatrix.


597:    Level: intermediate

599:    Notes: Unlike hypre's general semi-struct object consisting of a collection of structured-grid objects and unstructured
600:           grid objects, we will restrict the semi-struct objects to consist of only structured-grid components.

602:           Unlike the more general support for parts and blocks in hypre this allows only one part, and one block per process and requires the block
603:           be defined by a DMDA.

605:           The matrix needs a DMDA associated with it by either a call to MatSetupDM() or if the matrix is obtained from DMCreateMatrix()

607: M*/

611: PetscErrorCode  MatSetValuesLocal_HYPRESStruct_3d(Mat mat,PetscInt nrow,const PetscInt irow[],PetscInt ncol,const PetscInt icol[],const PetscScalar y[],InsertMode addv)
612: {
613:   PetscErrorCode    ierr;
614:   PetscInt          i,j,stencil,index[3];
615:   const PetscScalar *values = y;
616:   Mat_HYPRESStruct  *ex     = (Mat_HYPRESStruct*) mat->data;

618:   PetscInt part = 0;          /* Petsc sstruct interface only allows 1 part */
619:   PetscInt ordering;
620:   PetscInt grid_rank, to_grid_rank;
621:   PetscInt var_type, to_var_type;
622:   PetscInt to_var_entry = 0;

624:   PetscInt nvars= ex->nvars;
625:   PetscInt row,*entries;

628:   PetscMalloc1(7*nvars,&entries);
629:   ordering= ex->dofs_order;  /* ordering= 0   nodal ordering
630:                                           1   variable ordering */
631:   /* stencil entries are orderer by variables: var0_stencil0, var0_stencil1, ..., var0_stencil6, var1_stencil0, var1_stencil1, ...  */

633:   if (!ordering) {
634:     for (i=0; i<nrow; i++) {
635:       grid_rank= irow[i]/nvars;
636:       var_type = (irow[i] % nvars);

638:       for (j=0; j<ncol; j++) {
639:         to_grid_rank= icol[j]/nvars;
640:         to_var_type = (icol[j] % nvars);

642:         to_var_entry= to_var_entry*7;
643:         entries[j]  = to_var_entry;

645:         stencil = to_grid_rank-grid_rank;
646:         if (!stencil) {
647:           entries[j] += 3;
648:         } else if (stencil == -1) {
649:           entries[j] += 2;
650:         } else if (stencil == 1) {
651:           entries[j] += 4;
652:         } else if (stencil == -ex->gnx) {
653:           entries[j] += 1;
654:         } else if (stencil == ex->gnx) {
655:           entries[j] += 5;
656:         } else if (stencil == -ex->gnxgny) {
657:           entries[j] += 0;
658:         } else if (stencil == ex->gnxgny) {
659:           entries[j] += 6;
660:         } else SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Local row %D local column %D have bad stencil %D",irow[i],icol[j],stencil);
661:       }

663:       row      = ex->gindices[grid_rank] - ex->rstart;
664:       index[0] = ex->xs + (row % ex->nx);
665:       index[1] = ex->ys + ((row/ex->nx) % ex->ny);
666:       index[2] = ex->zs + (row/(ex->nxny));

668:       if (addv == ADD_VALUES) {
669:         PetscStackCallStandard(HYPRE_SStructMatrixAddToValues,(ex->ss_mat,part,(HYPRE_Int *)index,var_type,ncol,(HYPRE_Int *)entries,(PetscScalar*)values));
670:       } else {
671:         PetscStackCallStandard(HYPRE_SStructMatrixSetValues,(ex->ss_mat,part,(HYPRE_Int *)index,var_type,ncol,(HYPRE_Int *)entries,(PetscScalar*)values));
672:       }
673:       values += ncol;
674:     }
675:   } else {
676:     for (i=0; i<nrow; i++) {
677:       var_type = irow[i]/(ex->gnxgnygnz);
678:       grid_rank= irow[i] - var_type*(ex->gnxgnygnz);

680:       for (j=0; j<ncol; j++) {
681:         to_var_type = icol[j]/(ex->gnxgnygnz);
682:         to_grid_rank= icol[j] - to_var_type*(ex->gnxgnygnz);

684:         to_var_entry= to_var_entry*7;
685:         entries[j]  = to_var_entry;

687:         stencil = to_grid_rank-grid_rank;
688:         if (!stencil) {
689:           entries[j] += 3;
690:         } else if (stencil == -1) {
691:           entries[j] += 2;
692:         } else if (stencil == 1) {
693:           entries[j] += 4;
694:         } else if (stencil == -ex->gnx) {
695:           entries[j] += 1;
696:         } else if (stencil == ex->gnx) {
697:           entries[j] += 5;
698:         } else if (stencil == -ex->gnxgny) {
699:           entries[j] += 0;
700:         } else if (stencil == ex->gnxgny) {
701:           entries[j] += 6;
702:         } else SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Local row %D local column %D have bad stencil %D",irow[i],icol[j],stencil);
703:       }

705:       row      = ex->gindices[grid_rank] - ex->rstart;
706:       index[0] = ex->xs + (row % ex->nx);
707:       index[1] = ex->ys + ((row/ex->nx) % ex->ny);
708:       index[2] = ex->zs + (row/(ex->nxny));

710:       if (addv == ADD_VALUES) {
711:         PetscStackCallStandard(HYPRE_SStructMatrixAddToValues,(ex->ss_mat,part,(HYPRE_Int *)index,var_type,ncol,(HYPRE_Int *)entries,(PetscScalar*)values));
712:       } else {
713:         PetscStackCallStandard(HYPRE_SStructMatrixSetValues,(ex->ss_mat,part,(HYPRE_Int *)index,var_type,ncol,(HYPRE_Int *)entries,(PetscScalar*)values));
714:       }
715:       values += ncol;
716:     }
717:   }
718:   PetscFree(entries);
719:   return(0);
720: }

724: PetscErrorCode  MatZeroRowsLocal_HYPRESStruct_3d(Mat mat,PetscInt nrow,const PetscInt irow[],PetscScalar d,Vec x,Vec b)
725: {
726:   PetscErrorCode   ierr;
727:   PetscInt         i,index[3];
728:   PetscScalar      **values;
729:   Mat_HYPRESStruct *ex = (Mat_HYPRESStruct*) mat->data;

731:   PetscInt part     = 0;      /* Petsc sstruct interface only allows 1 part */
732:   PetscInt ordering = ex->dofs_order;
733:   PetscInt grid_rank;
734:   PetscInt var_type;
735:   PetscInt nvars= ex->nvars;
736:   PetscInt row,*entries;

739:   if (x && b) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"No support");
740:   PetscMalloc1(7*nvars,&entries);

742:   PetscMalloc1(nvars,&values);
743:   PetscMalloc1(7*nvars*nvars,&values[0]);
744:   for (i=1; i<nvars; i++) {
745:     values[i] = values[i-1] + nvars*7;
746:   }

748:   for (i=0; i< nvars; i++) {
749:     PetscMemzero(values[i],nvars*7*sizeof(PetscScalar));
750:     *(values[i]+3) = d;
751:   }

753:   for (i= 0; i< nvars*7; i++) entries[i] = i;

755:   if (!ordering) {
756:     for (i=0; i<nrow; i++) {
757:       grid_rank = irow[i] / nvars;
758:       var_type  =(irow[i] % nvars);

760:       row      = ex->gindices[grid_rank] - ex->rstart;
761:       index[0] = ex->xs + (row % ex->nx);
762:       index[1] = ex->ys + ((row/ex->nx) % ex->ny);
763:       index[2] = ex->zs + (row/(ex->nxny));
764:       PetscStackCallStandard(HYPRE_SStructMatrixSetValues,(ex->ss_mat,part,(HYPRE_Int *)index,var_type,7*nvars,(HYPRE_Int *)entries,values[var_type]));
765:     }
766:   } else {
767:     for (i=0; i<nrow; i++) {
768:       var_type = irow[i]/(ex->gnxgnygnz);
769:       grid_rank= irow[i] - var_type*(ex->gnxgnygnz);

771:       row      = ex->gindices[grid_rank] - ex->rstart;
772:       index[0] = ex->xs + (row % ex->nx);
773:       index[1] = ex->ys + ((row/ex->nx) % ex->ny);
774:       index[2] = ex->zs + (row/(ex->nxny));
775:       PetscStackCallStandard(HYPRE_SStructMatrixSetValues,(ex->ss_mat,part,(HYPRE_Int *)index,var_type,7*nvars,(HYPRE_Int *)entries,values[var_type]));
776:     }
777:   }
778:   PetscStackCallStandard(HYPRE_SStructMatrixAssemble,(ex->ss_mat));
779:   PetscFree(values[0]);
780:   PetscFree(values);
781:   PetscFree(entries);
782:   return(0);
783: }

787: PetscErrorCode MatZeroEntries_HYPRESStruct_3d(Mat mat)
788: {
789:   PetscErrorCode   ierr;
790:   Mat_HYPRESStruct *ex  = (Mat_HYPRESStruct*) mat->data;
791:   PetscInt         nvars= ex->nvars;
792:   PetscInt         size;
793:   PetscInt         part= 0;   /* only one part */

796:   size = ((ex->hbox.imax[0])-(ex->hbox.imin[0])+1)*((ex->hbox.imax[1])-(ex->hbox.imin[1])+1)*((ex->hbox.imax[2])-(ex->hbox.imin[2])+1);
797:   {
798:     PetscInt    i,*entries,iupper[3],ilower[3];
799:     PetscScalar *values;


802:     for (i= 0; i< 3; i++) {
803:       ilower[i]= ex->hbox.imin[i];
804:       iupper[i]= ex->hbox.imax[i];
805:     }

807:     PetscMalloc2(nvars*7,&entries,nvars*7*size,&values);
808:     for (i= 0; i< nvars*7; i++) entries[i]= i;
809:     PetscMemzero(values,nvars*7*size*sizeof(PetscScalar));

811:     for (i= 0; i< nvars; i++) {
812:       PetscStackCallStandard(HYPRE_SStructMatrixSetBoxValues,(ex->ss_mat,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,nvars*7,(HYPRE_Int *)entries,values));
813:     }
814:     PetscFree2(entries,values);
815:   }
816:   PetscStackCallStandard(HYPRE_SStructMatrixAssemble,(ex->ss_mat));
817:   return(0);
818: }


823: static PetscErrorCode  MatSetupDM_HYPRESStruct(Mat mat,DM da)
824: {
825:   PetscErrorCode         ierr;
826:   Mat_HYPRESStruct       *ex = (Mat_HYPRESStruct*) mat->data;
827:   PetscInt               dim,dof,sw[3],nx,ny,nz;
828:   PetscInt               ilower[3],iupper[3],ssize,i;
829:   DMBoundaryType         px,py,pz;
830:   DMDAStencilType        st;
831:   PetscInt               nparts= 1;  /* assuming only one part */
832:   PetscInt               part  = 0;
833:   ISLocalToGlobalMapping ltog;
835:   ex->da = da;
836:   PetscObjectReference((PetscObject)da);

838:   DMDAGetInfo(ex->da,&dim,0,0,0,0,0,0,&dof,&sw[0],&px,&py,&pz,&st);
839:   DMDAGetCorners(ex->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);
840:   iupper[0] += ilower[0] - 1;
841:   iupper[1] += ilower[1] - 1;
842:   iupper[2] += ilower[2] - 1;
843:   /* the hypre_Box is used to zero out the matrix entries in MatZeroValues() */
844:   ex->hbox.imin[0] = ilower[0];
845:   ex->hbox.imin[1] = ilower[1];
846:   ex->hbox.imin[2] = ilower[2];
847:   ex->hbox.imax[0] = iupper[0];
848:   ex->hbox.imax[1] = iupper[1];
849:   ex->hbox.imax[2] = iupper[2];

851:   ex->dofs_order = 0;

853:   /* assuming that the same number of dofs on each gridpoint. Also assume all cell-centred based */
854:   ex->nvars= dof;

856:   /* create the hypre grid object and set its information */
857:   if (px || py || pz) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Ask us to add periodic support by calling HYPRE_SStructGridSetPeriodic()");
858:   PetscStackCallStandard(HYPRE_SStructGridCreate,(ex->hcomm,dim,nparts,&ex->ss_grid));

860:   PetscStackCallStandard(HYPRE_SStructGridSetExtents,(ex->ss_grid,part,ex->hbox.imin,ex->hbox.imax));

862:   {
863:     HYPRE_SStructVariable *vartypes;
864:     PetscMalloc1(ex->nvars,&vartypes);
865:     for (i= 0; i< ex->nvars; i++) vartypes[i]= HYPRE_SSTRUCT_VARIABLE_CELL;
866:     PetscStackCallStandard(HYPRE_SStructGridSetVariables,(ex->ss_grid, part, ex->nvars,vartypes));
867:     PetscFree(vartypes);
868:   }
869:   PetscStackCallStandard(HYPRE_SStructGridAssemble,(ex->ss_grid));

871:   sw[1] = sw[0];
872:   sw[2] = sw[1];
873:   /* PetscStackCallStandard(HYPRE_SStructGridSetNumGhost,(ex->ss_grid,sw)); */

875:   /* create the hypre stencil object and set its information */
876:   if (sw[0] > 1) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Ask us to add support for wider stencils");
877:   if (st == DMDA_STENCIL_BOX) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Ask us to add support for box stencils");

879:   if (dim == 1) {
880:     PetscInt offsets[3][1] = {{-1},{0},{1}};
881:     PetscInt j, cnt;

883:     ssize = 3*(ex->nvars);
884:     PetscStackCallStandard(HYPRE_SStructStencilCreate,(dim,ssize,&ex->ss_stencil));
885:     cnt= 0;
886:     for (i = 0; i < (ex->nvars); i++) {
887:       for (j = 0; j < 3; j++) {
888:         PetscStackCallStandard(HYPRE_SStructStencilSetEntry,(ex->ss_stencil, cnt, (HYPRE_Int *)offsets[j], i));
889:         cnt++;
890:       }
891:     }
892:   } else if (dim == 2) {
893:     PetscInt offsets[5][2] = {{0,-1},{-1,0},{0,0},{1,0},{0,1}};
894:     PetscInt j, cnt;

896:     ssize = 5*(ex->nvars);
897:     PetscStackCallStandard(HYPRE_SStructStencilCreate,(dim,ssize,&ex->ss_stencil));
898:     cnt= 0;
899:     for (i= 0; i< (ex->nvars); i++) {
900:       for (j= 0; j< 5; j++) {
901:         PetscStackCallStandard(HYPRE_SStructStencilSetEntry,(ex->ss_stencil, cnt, (HYPRE_Int *)offsets[j], i));
902:         cnt++;
903:       }
904:     }
905:   } else if (dim == 3) {
906:     PetscInt offsets[7][3] = {{0,0,-1},{0,-1,0},{-1,0,0},{0,0,0},{1,0,0},{0,1,0},{0,0,1}};
907:     PetscInt j, cnt;

909:     ssize = 7*(ex->nvars);
910:     PetscStackCallStandard(HYPRE_SStructStencilCreate,(dim,ssize,&ex->ss_stencil));
911:     cnt= 0;
912:     for (i= 0; i< (ex->nvars); i++) {
913:       for (j= 0; j< 7; j++) {
914:         PetscStackCallStandard(HYPRE_SStructStencilSetEntry,(ex->ss_stencil, cnt, (HYPRE_Int *)offsets[j], i));
915:         cnt++;
916:       }
917:     }
918:   }

920:   /* create the HYPRE graph */
921:   PetscStackCallStandard(HYPRE_SStructGraphCreate,(ex->hcomm, ex->ss_grid, &(ex->ss_graph)));

923:   /* set the stencil graph. Note that each variable has the same graph. This means that each
924:      variable couples to all the other variable and with the same stencil pattern. */
925:   for (i= 0; i< (ex->nvars); i++) {
926:     PetscStackCallStandard(HYPRE_SStructGraphSetStencil,(ex->ss_graph,part,i,ex->ss_stencil));
927:   }
928:   PetscStackCallStandard(HYPRE_SStructGraphAssemble,(ex->ss_graph));

930:   /* create the HYPRE sstruct vectors for rhs and solution */
931:   PetscStackCallStandard(HYPRE_SStructVectorCreate,(ex->hcomm,ex->ss_grid,&ex->ss_b));
932:   PetscStackCallStandard(HYPRE_SStructVectorCreate,(ex->hcomm,ex->ss_grid,&ex->ss_x));
933:   PetscStackCallStandard(HYPRE_SStructVectorInitialize,(ex->ss_b));
934:   PetscStackCallStandard(HYPRE_SStructVectorInitialize,(ex->ss_x));
935:   PetscStackCallStandard(HYPRE_SStructVectorAssemble,(ex->ss_b));
936:   PetscStackCallStandard(HYPRE_SStructVectorAssemble,(ex->ss_x));

938:   /* create the hypre matrix object and set its information */
939:   PetscStackCallStandard(HYPRE_SStructMatrixCreate,(ex->hcomm,ex->ss_graph,&ex->ss_mat));
940:   PetscStackCallStandard(HYPRE_SStructGridDestroy,(ex->ss_grid));
941:   PetscStackCallStandard(HYPRE_SStructStencilDestroy,(ex->ss_stencil));
942:   if (ex->needsinitialization) {
943:     PetscStackCallStandard(HYPRE_SStructMatrixInitialize,(ex->ss_mat));
944:     ex->needsinitialization = PETSC_FALSE;
945:   }

947:   /* set the global and local sizes of the matrix */
948:   DMDAGetCorners(da,0,0,0,&nx,&ny,&nz);
949:   MatSetSizes(mat,dof*nx*ny*nz,dof*nx*ny*nz,PETSC_DECIDE,PETSC_DECIDE);
950:   PetscLayoutSetBlockSize(mat->rmap,1);
951:   PetscLayoutSetBlockSize(mat->cmap,1);
952:   PetscLayoutSetUp(mat->rmap);
953:   PetscLayoutSetUp(mat->cmap);

955:   if (dim == 3) {
956:     mat->ops->setvalueslocal = MatSetValuesLocal_HYPRESStruct_3d;
957:     mat->ops->zerorowslocal  = MatZeroRowsLocal_HYPRESStruct_3d;
958:     mat->ops->zeroentries    = MatZeroEntries_HYPRESStruct_3d;

960:     MatZeroEntries_HYPRESStruct_3d(mat);
961:   } else SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Only support for 3d DMDA currently");

963:   /* get values that will be used repeatedly in MatSetValuesLocal() and MatZeroRowsLocal() repeatedly */
964:   MatGetOwnershipRange(mat,&ex->rstart,NULL);
965:   DMGetLocalToGlobalMapping(ex->da,&ltog);
966:   ISLocalToGlobalMappingGetIndices(ltog, (const PetscInt **) &ex->gindices);
967:   DMDAGetGhostCorners(ex->da,0,0,0,&ex->gnx,&ex->gnxgny,&ex->gnxgnygnz);

969:   ex->gnxgny    *= ex->gnx;
970:   ex->gnxgnygnz *= ex->gnxgny;

972:   DMDAGetCorners(ex->da,&ex->xs,&ex->ys,&ex->zs,&ex->nx,&ex->ny,&ex->nz);

974:   ex->nxny   = ex->nx*ex->ny;
975:   ex->nxnynz = ex->nz*ex->nxny;
976:   return(0);
977: }

981: PetscErrorCode MatMult_HYPRESStruct(Mat A,Vec x,Vec y)
982: {
983:   PetscErrorCode   ierr;
984:   PetscScalar      *xx,*yy;
985:   PetscInt         ilower[3],iupper[3];
986:   Mat_HYPRESStruct *mx     = (Mat_HYPRESStruct*)(A->data);
987:   PetscInt         ordering= mx->dofs_order;
988:   PetscInt         nvars   = mx->nvars;
989:   PetscInt         part    = 0;
990:   PetscInt         size;
991:   PetscInt         i;

994:   DMDAGetCorners(mx->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);
995:   iupper[0] += ilower[0] - 1;
996:   iupper[1] += ilower[1] - 1;
997:   iupper[2] += ilower[2] - 1;

999:   size= 1;
1000:   for (i = 0; i < 3; i++) size *= (iupper[i]-ilower[i]+1);

1002:   /* copy x values over to hypre for variable ordering */
1003:   if (ordering) {
1004:     PetscStackCallStandard(HYPRE_SStructVectorSetConstantValues,(mx->ss_b,0.0));
1005:     VecGetArray(x,&xx);
1006:     for (i= 0; i< nvars; i++) {
1007:       PetscStackCallStandard(HYPRE_SStructVectorSetBoxValues,(mx->ss_b,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,xx+(size*i)));
1008:     }
1009:     VecRestoreArray(x,&xx);
1010:     PetscStackCallStandard(HYPRE_SStructVectorAssemble,(mx->ss_b));
1011:     PetscStackCallStandard(HYPRE_SStructMatrixMatvec,(1.0,mx->ss_mat,mx->ss_b,0.0,mx->ss_x));

1013:     /* copy solution values back to PETSc */
1014:     VecGetArray(y,&yy);
1015:     for (i= 0; i< nvars; i++) {
1016:       PetscStackCallStandard(HYPRE_SStructVectorGetBoxValues,(mx->ss_x,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,yy+(size*i)));
1017:     }
1018:     VecRestoreArray(y,&yy);
1019:   } else {      /* nodal ordering must be mapped to variable ordering for sys_pfmg */
1020:     PetscScalar *z;
1021:     PetscInt    j, k;

1023:     PetscMalloc1(nvars*size,&z);
1024:     PetscStackCallStandard(HYPRE_SStructVectorSetConstantValues,(mx->ss_b,0.0));
1025:     VecGetArray(x,&xx);

1027:     /* transform nodal to hypre's variable ordering for sys_pfmg */
1028:     for (i= 0; i< size; i++) {
1029:       k= i*nvars;
1030:       for (j= 0; j< nvars; j++) z[j*size+i]= xx[k+j];
1031:     }
1032:     for (i= 0; i< nvars; i++) {
1033:       PetscStackCallStandard(HYPRE_SStructVectorSetBoxValues,(mx->ss_b,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,z+(size*i)));
1034:     }
1035:     VecRestoreArray(x,&xx);
1036:     PetscStackCallStandard(HYPRE_SStructVectorAssemble,(mx->ss_b));
1037:     PetscStackCallStandard(HYPRE_SStructMatrixMatvec,(1.0,mx->ss_mat,mx->ss_b,0.0,mx->ss_x));

1039:     /* copy solution values back to PETSc */
1040:     VecGetArray(y,&yy);
1041:     for (i= 0; i< nvars; i++) {
1042:       PetscStackCallStandard(HYPRE_SStructVectorGetBoxValues,(mx->ss_x,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,z+(size*i)));
1043:     }
1044:     /* transform hypre's variable ordering for sys_pfmg to nodal ordering */
1045:     for (i= 0; i< size; i++) {
1046:       k= i*nvars;
1047:       for (j= 0; j< nvars; j++) yy[k+j]= z[j*size+i];
1048:     }
1049:     VecRestoreArray(y,&yy);
1050:     PetscFree(z);
1051:   }
1052:   return(0);
1053: }

1057: PetscErrorCode MatAssemblyEnd_HYPRESStruct(Mat mat,MatAssemblyType mode)
1058: {
1059:   Mat_HYPRESStruct *ex = (Mat_HYPRESStruct*) mat->data;
1060:   PetscErrorCode   ierr;

1063:   PetscStackCallStandard(HYPRE_SStructMatrixAssemble,(ex->ss_mat));
1064:   return(0);
1065: }

1069: PetscErrorCode MatZeroEntries_HYPRESStruct(Mat mat)
1070: {
1072:   /* before the DMDA is set to the matrix the zero doesn't need to do anything */
1073:   return(0);
1074: }


1079: PetscErrorCode MatDestroy_HYPRESStruct(Mat mat)
1080: {
1081:   Mat_HYPRESStruct *ex = (Mat_HYPRESStruct*) mat->data;
1082:   PetscErrorCode   ierr;

1085:   PetscStackCallStandard(HYPRE_SStructGraphDestroy,(ex->ss_graph));
1086:   PetscStackCallStandard(HYPRE_SStructMatrixDestroy,(ex->ss_mat));
1087:   PetscStackCallStandard(HYPRE_SStructVectorDestroy,(ex->ss_x));
1088:   PetscStackCallStandard(HYPRE_SStructVectorDestroy,(ex->ss_b));
1089:   return(0);
1090: }

1094: PETSC_EXTERN PetscErrorCode MatCreate_HYPRESStruct(Mat B)
1095: {
1096:   Mat_HYPRESStruct *ex;
1097:   PetscErrorCode   ierr;

1100:   PetscNewLog(B,&ex);
1101:   B->data      = (void*)ex;
1102:   B->rmap->bs  = 1;
1103:   B->assembled = PETSC_FALSE;

1105:   B->insertmode = NOT_SET_VALUES;

1107:   B->ops->assemblyend = MatAssemblyEnd_HYPRESStruct;
1108:   B->ops->mult        = MatMult_HYPRESStruct;
1109:   B->ops->zeroentries = MatZeroEntries_HYPRESStruct;
1110:   B->ops->destroy     = MatDestroy_HYPRESStruct;

1112:   ex->needsinitialization = PETSC_TRUE;

1114:   MPI_Comm_dup(PetscObjectComm((PetscObject)B),&(ex->hcomm));
1115:   PetscObjectComposeFunction((PetscObject)B,"MatSetupDM_C",MatSetupDM_HYPRESStruct);
1116:   PetscObjectChangeTypeName((PetscObject)B,MATHYPRESSTRUCT);
1117:   return(0);
1118: }