Actual source code: mhyp.c


  2: /*
  3:     Creates hypre ijmatrix from PETSc matrix
  4: */
  5: #include <petscsys.h>
  6: #include <petsc/private/petschypre.h>
  7: #include <petsc/private/matimpl.h>
  8: #include <petscdmda.h>
  9: #include <../src/dm/impls/da/hypre/mhyp.h>

 11: /* -----------------------------------------------------------------------------------------------------------------*/

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

 17:    Level: intermediate

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

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

 25: .seealso: MatCreate(), PCPFMG, MatSetDM(), DMCreateMatrix()
 26: M*/

 28: PetscErrorCode  MatSetValuesLocal_HYPREStruct_3d(Mat mat,PetscInt nrow,const PetscInt irow[],PetscInt ncol,const PetscInt icol[],const PetscScalar y[],InsertMode addv)
 29: {
 30:   HYPRE_Int       index[3],entries[7];
 31:   PetscInt        i,j,stencil,row;
 32:   HYPRE_Complex   *values = (HYPRE_Complex*)y;
 33:   Mat_HYPREStruct *ex     = (Mat_HYPREStruct*) mat->data;

 36:   if (PetscUnlikely(ncol > 7)) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"ncol %D > 7 too large",ncol);
 37:   for (i=0; i<nrow; i++) {
 38:     for (j=0; j<ncol; j++) {
 39:       stencil = icol[j] - irow[i];
 40:       if (!stencil) {
 41:         entries[j] = 3;
 42:       } else if (stencil == -1) {
 43:         entries[j] = 2;
 44:       } else if (stencil == 1) {
 45:         entries[j] = 4;
 46:       } else if (stencil == -ex->gnx) {
 47:         entries[j] = 1;
 48:       } else if (stencil == ex->gnx) {
 49:         entries[j] = 5;
 50:       } else if (stencil == -ex->gnxgny) {
 51:         entries[j] = 0;
 52:       } else if (stencil == ex->gnxgny) {
 53:         entries[j] = 6;
 54:       } else SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Local row %D local column %D have bad stencil %D",irow[i],icol[j],stencil);
 55:     }
 56:     row      = ex->gindices[irow[i]] - ex->rstart;
 57:     index[0] = (HYPRE_Int)(ex->xs + (row % ex->nx));
 58:     index[1] = (HYPRE_Int)(ex->ys + ((row/ex->nx) % ex->ny));
 59:     index[2] = (HYPRE_Int)(ex->zs + (row/(ex->nxny)));
 60:     if (addv == ADD_VALUES) {
 61:       PetscStackCallStandard(HYPRE_StructMatrixAddToValues,(ex->hmat,index,ncol,entries,values));
 62:     } else {
 63:       PetscStackCallStandard(HYPRE_StructMatrixSetValues,(ex->hmat,index,ncol,entries,values));
 64:     }
 65:     values += ncol;
 66:   }
 67:   return(0);
 68: }

 70: PetscErrorCode  MatZeroRowsLocal_HYPREStruct_3d(Mat mat,PetscInt nrow,const PetscInt irow[],PetscScalar d,Vec x,Vec b)
 71: {
 72:   PetscErrorCode  ierr;
 73:   HYPRE_Int       index[3],entries[7] = {0,1,2,3,4,5,6};
 74:   PetscInt        row,i;
 75:   HYPRE_Complex   values[7];
 76:   Mat_HYPREStruct *ex = (Mat_HYPREStruct*) mat->data;

 79:   if (x && b) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"No support");
 80:   PetscArrayzero(values,7);
 81:   PetscHYPREScalarCast(d,&values[3]);
 82:   for (i=0; i<nrow; i++) {
 83:     row      = ex->gindices[irow[i]] - ex->rstart;
 84:     index[0] = (HYPRE_Int)(ex->xs + (row % ex->nx));
 85:     index[1] = (HYPRE_Int)(ex->ys + ((row/ex->nx) % ex->ny));
 86:     index[2] = (HYPRE_Int)(ex->zs + (row/(ex->nxny)));
 87:     PetscStackCallStandard(HYPRE_StructMatrixSetValues,(ex->hmat,index,7,entries,values));
 88:   }
 89:   PetscStackCallStandard(HYPRE_StructMatrixAssemble,(ex->hmat));
 90:   return(0);
 91: }

 93: PetscErrorCode MatZeroEntries_HYPREStruct_3d(Mat mat)
 94: {
 95:   HYPRE_Int       indices[7] = {0,1,2,3,4,5,6};
 96:   Mat_HYPREStruct *ex        = (Mat_HYPREStruct*) mat->data;

 99:   /* hypre has no public interface to do this */
100:   PetscStackCallStandard(hypre_StructMatrixClearBoxValues,(ex->hmat,&ex->hbox,7,indices,0,1));
101:   PetscStackCallStandard(HYPRE_StructMatrixAssemble,(ex->hmat));
102:   return(0);
103: }

105: static PetscErrorCode  MatSetUp_HYPREStruct(Mat mat)
106: {
107:   PetscErrorCode         ierr;
108:   Mat_HYPREStruct        *ex = (Mat_HYPREStruct*) mat->data;
109:   HYPRE_Int              sw[6];
110:   HYPRE_Int              hlower[3],hupper[3];
111:   PetscInt               dim,dof,psw,nx,ny,nz,ilower[3],iupper[3],ssize,i;
112:   DMBoundaryType         px,py,pz;
113:   DMDAStencilType        st;
114:   ISLocalToGlobalMapping ltog;
115:   DM                     da;

118:   MatGetDM(mat,(DM*)&da);
119:   ex->da = da;
120:   PetscObjectReference((PetscObject)da);

122:   DMDAGetInfo(ex->da,&dim,0,0,0,0,0,0,&dof,&psw,&px,&py,&pz,&st);
123:   DMDAGetCorners(ex->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);

125:   /* when HYPRE_MIXEDINT is defined, sizeof(HYPRE_Int) == 32 */
126:   iupper[0] += ilower[0] - 1;
127:   iupper[1] += ilower[1] - 1;
128:   iupper[2] += ilower[2] - 1;
129:   hlower[0]  = (HYPRE_Int)ilower[0];
130:   hlower[1]  = (HYPRE_Int)ilower[1];
131:   hlower[2]  = (HYPRE_Int)ilower[2];
132:   hupper[0]  = (HYPRE_Int)iupper[0];
133:   hupper[1]  = (HYPRE_Int)iupper[1];
134:   hupper[2]  = (HYPRE_Int)iupper[2];

136:   /* the hypre_Box is used to zero out the matrix entries in MatZeroValues() */
137:   ex->hbox.imin[0] = hlower[0];
138:   ex->hbox.imin[1] = hlower[1];
139:   ex->hbox.imin[2] = hlower[2];
140:   ex->hbox.imax[0] = hupper[0];
141:   ex->hbox.imax[1] = hupper[1];
142:   ex->hbox.imax[2] = hupper[2];

144:   /* create the hypre grid object and set its information */
145:   if (dof > 1) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Currently only support for scalar problems");
146:   if (px || py || pz) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Ask us to add periodic support by calling HYPRE_StructGridSetPeriodic()");
147:   PetscStackCallStandard(HYPRE_StructGridCreate,(ex->hcomm,dim,&ex->hgrid));
148:   PetscStackCallStandard(HYPRE_StructGridSetExtents,(ex->hgrid,hlower,hupper));
149:   PetscStackCallStandard(HYPRE_StructGridAssemble,(ex->hgrid));

151:   sw[5] = sw[4] = sw[3] = sw[2] = sw[1] = sw[0] = psw;
152:   PetscStackCallStandard(HYPRE_StructGridSetNumGhost,(ex->hgrid,sw));

154:   /* create the hypre stencil object and set its information */
155:   if (sw[0] > 1) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Ask us to add support for wider stencils");
156:   if (st == DMDA_STENCIL_BOX) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Ask us to add support for box stencils");
157:   if (dim == 1) {
158:     HYPRE_Int offsets[3][1] = {{-1},{0},{1}};
159:     ssize = 3;
160:     PetscStackCallStandard(HYPRE_StructStencilCreate,(dim,ssize,&ex->hstencil));
161:     for (i=0; i<ssize; i++) {
162:       PetscStackCallStandard(HYPRE_StructStencilSetElement,(ex->hstencil,i,offsets[i]));
163:     }
164:   } else if (dim == 2) {
165:     HYPRE_Int offsets[5][2] = {{0,-1},{-1,0},{0,0},{1,0},{0,1}};
166:     ssize = 5;
167:     PetscStackCallStandard(HYPRE_StructStencilCreate,(dim,ssize,&ex->hstencil));
168:     for (i=0; i<ssize; i++) {
169:       PetscStackCallStandard(HYPRE_StructStencilSetElement,(ex->hstencil,i,offsets[i]));
170:     }
171:   } else if (dim == 3) {
172:     HYPRE_Int offsets[7][3] = {{0,0,-1},{0,-1,0},{-1,0,0},{0,0,0},{1,0,0},{0,1,0},{0,0,1}};
173:     ssize = 7;
174:     PetscStackCallStandard(HYPRE_StructStencilCreate,(dim,ssize,&ex->hstencil));
175:     for (i=0; i<ssize; i++) {
176:       PetscStackCallStandard(HYPRE_StructStencilSetElement,(ex->hstencil,i,offsets[i]));
177:     }
178:   }

180:   /* create the HYPRE vector for rhs and solution */
181:   PetscStackCallStandard(HYPRE_StructVectorCreate,(ex->hcomm,ex->hgrid,&ex->hb));
182:   PetscStackCallStandard(HYPRE_StructVectorCreate,(ex->hcomm,ex->hgrid,&ex->hx));
183:   PetscStackCallStandard(HYPRE_StructVectorInitialize,(ex->hb));
184:   PetscStackCallStandard(HYPRE_StructVectorInitialize,(ex->hx));
185:   PetscStackCallStandard(HYPRE_StructVectorAssemble,(ex->hb));
186:   PetscStackCallStandard(HYPRE_StructVectorAssemble,(ex->hx));

188:   /* create the hypre matrix object and set its information */
189:   PetscStackCallStandard(HYPRE_StructMatrixCreate,(ex->hcomm,ex->hgrid,ex->hstencil,&ex->hmat));
190:   PetscStackCallStandard(HYPRE_StructGridDestroy,(ex->hgrid));
191:   PetscStackCallStandard(HYPRE_StructStencilDestroy,(ex->hstencil));
192:   if (ex->needsinitialization) {
193:     PetscStackCallStandard(HYPRE_StructMatrixInitialize,(ex->hmat));
194:     ex->needsinitialization = PETSC_FALSE;
195:   }

197:   /* set the global and local sizes of the matrix */
198:   DMDAGetCorners(da,0,0,0,&nx,&ny,&nz);
199:   MatSetSizes(mat,dof*nx*ny*nz,dof*nx*ny*nz,PETSC_DECIDE,PETSC_DECIDE);
200:   PetscLayoutSetBlockSize(mat->rmap,1);
201:   PetscLayoutSetBlockSize(mat->cmap,1);
202:   PetscLayoutSetUp(mat->rmap);
203:   PetscLayoutSetUp(mat->cmap);
204:   mat->preallocated = PETSC_TRUE;

206:   if (dim == 3) {
207:     mat->ops->setvalueslocal = MatSetValuesLocal_HYPREStruct_3d;
208:     mat->ops->zerorowslocal  = MatZeroRowsLocal_HYPREStruct_3d;
209:     mat->ops->zeroentries    = MatZeroEntries_HYPREStruct_3d;

211:     /*        MatZeroEntries_HYPREStruct_3d(mat); */
212:   } else SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Only support for 3d DMDA currently");

214:   /* get values that will be used repeatedly in MatSetValuesLocal() and MatZeroRowsLocal() repeatedly */
215:   MatGetOwnershipRange(mat,&ex->rstart,NULL);
216:   DMGetLocalToGlobalMapping(ex->da,&ltog);
217:   ISLocalToGlobalMappingGetIndices(ltog, (const PetscInt **) &ex->gindices);
218:   DMDAGetGhostCorners(ex->da,0,0,0,&ex->gnx,&ex->gnxgny,0);
219:   ex->gnxgny *= ex->gnx;
220:   DMDAGetCorners(ex->da,&ex->xs,&ex->ys,&ex->zs,&ex->nx,&ex->ny,0);
221:   ex->nxny    = ex->nx*ex->ny;
222:   return(0);
223: }

225: PetscErrorCode MatMult_HYPREStruct(Mat A,Vec x,Vec y)
226: {
227:   PetscErrorCode    ierr;
228:   const PetscScalar *xx;
229:   PetscScalar       *yy;
230:   PetscInt          ilower[3],iupper[3];
231:   HYPRE_Int         hlower[3],hupper[3];
232:   Mat_HYPREStruct   *mx = (Mat_HYPREStruct*)(A->data);

235:   DMDAGetCorners(mx->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);
236:   /* when HYPRE_MIXEDINT is defined, sizeof(HYPRE_Int) == 32 */
237:   iupper[0] += ilower[0] - 1;
238:   iupper[1] += ilower[1] - 1;
239:   iupper[2] += ilower[2] - 1;
240:   hlower[0]  = (HYPRE_Int)ilower[0];
241:   hlower[1]  = (HYPRE_Int)ilower[1];
242:   hlower[2]  = (HYPRE_Int)ilower[2];
243:   hupper[0]  = (HYPRE_Int)iupper[0];
244:   hupper[1]  = (HYPRE_Int)iupper[1];
245:   hupper[2]  = (HYPRE_Int)iupper[2];

247:   /* copy x values over to hypre */
248:   PetscStackCallStandard(HYPRE_StructVectorSetConstantValues,(mx->hb,0.0));
249:   VecGetArrayRead(x,&xx);
250:   PetscStackCallStandard(HYPRE_StructVectorSetBoxValues,(mx->hb,hlower,hupper,(HYPRE_Complex*)xx));
251:   VecRestoreArrayRead(x,&xx);
252:   PetscStackCallStandard(HYPRE_StructVectorAssemble,(mx->hb));
253:   PetscStackCallStandard(HYPRE_StructMatrixMatvec,(1.0,mx->hmat,mx->hb,0.0,mx->hx));

255:   /* copy solution values back to PETSc */
256:   VecGetArray(y,&yy);
257:   PetscStackCallStandard(HYPRE_StructVectorGetBoxValues,(mx->hx,hlower,hupper,(HYPRE_Complex*)yy));
258:   VecRestoreArray(y,&yy);
259:   return(0);
260: }

262: PetscErrorCode MatAssemblyEnd_HYPREStruct(Mat mat,MatAssemblyType mode)
263: {
264:   Mat_HYPREStruct *ex = (Mat_HYPREStruct*) mat->data;

267:   PetscStackCallStandard(HYPRE_StructMatrixAssemble,(ex->hmat));
268:   /* PetscStackCallStandard(HYPRE_StructMatrixPrint,("dummy",ex->hmat,0)); */
269:   return(0);
270: }

272: PetscErrorCode MatZeroEntries_HYPREStruct(Mat mat)
273: {
275:   /* before the DMDA is set to the matrix the zero doesn't need to do anything */
276:   return(0);
277: }

279: PetscErrorCode MatDestroy_HYPREStruct(Mat mat)
280: {
281:   Mat_HYPREStruct *ex = (Mat_HYPREStruct*) mat->data;
282:   PetscErrorCode  ierr;

285:   PetscStackCallStandard(HYPRE_StructMatrixDestroy,(ex->hmat));
286:   PetscStackCallStandard(HYPRE_StructVectorDestroy,(ex->hx));
287:   PetscStackCallStandard(HYPRE_StructVectorDestroy,(ex->hb));
288:   PetscObjectDereference((PetscObject)ex->da);
289:   MPI_Comm_free(&(ex->hcomm));
290:   PetscFree(ex);
291:   return(0);
292: }

294: PETSC_EXTERN PetscErrorCode MatCreate_HYPREStruct(Mat B)
295: {
296:   Mat_HYPREStruct *ex;
297:   PetscErrorCode  ierr;

300:   PetscNewLog(B,&ex);
301:   B->data      = (void*)ex;
302:   B->rmap->bs  = 1;
303:   B->assembled = PETSC_FALSE;

305:   B->insertmode = NOT_SET_VALUES;

307:   B->ops->assemblyend = MatAssemblyEnd_HYPREStruct;
308:   B->ops->mult        = MatMult_HYPREStruct;
309:   B->ops->zeroentries = MatZeroEntries_HYPREStruct;
310:   B->ops->destroy     = MatDestroy_HYPREStruct;
311:   B->ops->setup       = MatSetUp_HYPREStruct;

313:   ex->needsinitialization = PETSC_TRUE;

315:   MPI_Comm_dup(PetscObjectComm((PetscObject)B),&(ex->hcomm));
316:   PetscObjectChangeTypeName((PetscObject)B,MATHYPRESTRUCT);
317:   return(0);
318: }

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

324:    Level: intermediate

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

330:           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
331:           be defined by a DMDA.

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

335: M*/

337: PetscErrorCode  MatSetValuesLocal_HYPRESStruct_3d(Mat mat,PetscInt nrow,const PetscInt irow[],PetscInt ncol,const PetscInt icol[],const PetscScalar y[],InsertMode addv)
338: {
339:   PetscErrorCode    ierr;
340:   HYPRE_Int         index[3],*entries;
341:   PetscInt          i,j,stencil;
342:   HYPRE_Complex     *values = (HYPRE_Complex*)y;
343:   Mat_HYPRESStruct  *ex     = (Mat_HYPRESStruct*) mat->data;

345:   PetscInt          part = 0;          /* Petsc sstruct interface only allows 1 part */
346:   PetscInt          ordering;
347:   PetscInt          grid_rank, to_grid_rank;
348:   PetscInt          var_type, to_var_type;
349:   PetscInt          to_var_entry = 0;
350:   PetscInt          nvars= ex->nvars;
351:   PetscInt          row;

354:   PetscMalloc1(7*nvars,&entries);
355:   ordering = ex->dofs_order;  /* ordering= 0   nodal ordering
356:                                            1   variable ordering */
357:   /* stencil entries are ordered by variables: var0_stencil0, var0_stencil1, ..., var0_stencil6, var1_stencil0, var1_stencil1, ...  */
358:   if (!ordering) {
359:     for (i=0; i<nrow; i++) {
360:       grid_rank= irow[i]/nvars;
361:       var_type = (irow[i] % nvars);

363:       for (j=0; j<ncol; j++) {
364:         to_grid_rank= icol[j]/nvars;
365:         to_var_type = (icol[j] % nvars);

367:         to_var_entry= to_var_entry*7;
368:         entries[j]  = to_var_entry;

370:         stencil = to_grid_rank-grid_rank;
371:         if (!stencil) {
372:           entries[j] += 3;
373:         } else if (stencil == -1) {
374:           entries[j] += 2;
375:         } else if (stencil == 1) {
376:           entries[j] += 4;
377:         } else if (stencil == -ex->gnx) {
378:           entries[j] += 1;
379:         } else if (stencil == ex->gnx) {
380:           entries[j] += 5;
381:         } else if (stencil == -ex->gnxgny) {
382:           entries[j] += 0;
383:         } else if (stencil == ex->gnxgny) {
384:           entries[j] += 6;
385:         } else SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Local row %D local column %D have bad stencil %D",irow[i],icol[j],stencil);
386:       }

388:       row      = ex->gindices[grid_rank] - ex->rstart;
389:       index[0] = (HYPRE_Int)(ex->xs + (row % ex->nx));
390:       index[1] = (HYPRE_Int)(ex->ys + ((row/ex->nx) % ex->ny));
391:       index[2] = (HYPRE_Int)(ex->zs + (row/(ex->nxny)));

393:       if (addv == ADD_VALUES) {
394:         PetscStackCallStandard(HYPRE_SStructMatrixAddToValues,(ex->ss_mat,part,index,var_type,ncol,entries,values));
395:       } else {
396:         PetscStackCallStandard(HYPRE_SStructMatrixSetValues,(ex->ss_mat,part,index,var_type,ncol,entries,values));
397:       }
398:       values += ncol;
399:     }
400:   } else {
401:     for (i=0; i<nrow; i++) {
402:       var_type = irow[i]/(ex->gnxgnygnz);
403:       grid_rank= irow[i] - var_type*(ex->gnxgnygnz);

405:       for (j=0; j<ncol; j++) {
406:         to_var_type = icol[j]/(ex->gnxgnygnz);
407:         to_grid_rank= icol[j] - to_var_type*(ex->gnxgnygnz);

409:         to_var_entry= to_var_entry*7;
410:         entries[j]  = to_var_entry;

412:         stencil = to_grid_rank-grid_rank;
413:         if (!stencil) {
414:           entries[j] += 3;
415:         } else if (stencil == -1) {
416:           entries[j] += 2;
417:         } else if (stencil == 1) {
418:           entries[j] += 4;
419:         } else if (stencil == -ex->gnx) {
420:           entries[j] += 1;
421:         } else if (stencil == ex->gnx) {
422:           entries[j] += 5;
423:         } else if (stencil == -ex->gnxgny) {
424:           entries[j] += 0;
425:         } else if (stencil == ex->gnxgny) {
426:           entries[j] += 6;
427:         } else SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Local row %D local column %D have bad stencil %D",irow[i],icol[j],stencil);
428:       }

430:       row      = ex->gindices[grid_rank] - ex->rstart;
431:       index[0] = (HYPRE_Int)(ex->xs + (row % ex->nx));
432:       index[1] = (HYPRE_Int)(ex->ys + ((row/ex->nx) % ex->ny));
433:       index[2] = (HYPRE_Int)(ex->zs + (row/(ex->nxny)));

435:       if (addv == ADD_VALUES) {
436:         PetscStackCallStandard(HYPRE_SStructMatrixAddToValues,(ex->ss_mat,part,index,var_type,ncol,entries,values));
437:       } else {
438:         PetscStackCallStandard(HYPRE_SStructMatrixSetValues,(ex->ss_mat,part,index,var_type,ncol,entries,values));
439:       }
440:       values += ncol;
441:     }
442:   }
443:   PetscFree(entries);
444:   return(0);
445: }

447: PetscErrorCode  MatZeroRowsLocal_HYPRESStruct_3d(Mat mat,PetscInt nrow,const PetscInt irow[],PetscScalar d,Vec x,Vec b)
448: {
449:   PetscErrorCode   ierr;
450:   HYPRE_Int        index[3],*entries;
451:   PetscInt         i;
452:   HYPRE_Complex    **values;
453:   Mat_HYPRESStruct *ex = (Mat_HYPRESStruct*) mat->data;

455:   PetscInt         part     = 0;      /* Petsc sstruct interface only allows 1 part */
456:   PetscInt         ordering = ex->dofs_order;
457:   PetscInt         grid_rank;
458:   PetscInt         var_type;
459:   PetscInt         nvars= ex->nvars;
460:   PetscInt         row;

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

466:   PetscMalloc1(nvars,&values);
467:   PetscMalloc1(7*nvars*nvars,&values[0]);
468:   for (i=1; i<nvars; i++) {
469:     values[i] = values[i-1] + nvars*7;
470:   }

472:   for (i=0; i< nvars; i++) {
473:     PetscArrayzero(values[i],nvars*7*sizeof(HYPRE_Complex));
474:     PetscHYPREScalarCast(d,values[i]+3);
475:   }

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

479:   if (!ordering) {
480:     for (i=0; i<nrow; i++) {
481:       grid_rank = irow[i] / nvars;
482:       var_type  = (irow[i] % nvars);

484:       row      = ex->gindices[grid_rank] - ex->rstart;
485:       index[0] = (HYPRE_Int)(ex->xs + (row % ex->nx));
486:       index[1] = (HYPRE_Int)(ex->ys + ((row/ex->nx) % ex->ny));
487:       index[2] = (HYPRE_Int)(ex->zs + (row/(ex->nxny)));
488:       PetscStackCallStandard(HYPRE_SStructMatrixSetValues,(ex->ss_mat,part,index,var_type,7*nvars,entries,values[var_type]));
489:     }
490:   } else {
491:     for (i=0; i<nrow; i++) {
492:       var_type = irow[i]/(ex->gnxgnygnz);
493:       grid_rank= irow[i] - var_type*(ex->gnxgnygnz);

495:       row      = ex->gindices[grid_rank] - ex->rstart;
496:       index[0] = (HYPRE_Int)(ex->xs + (row % ex->nx));
497:       index[1] = (HYPRE_Int)(ex->ys + ((row/ex->nx) % ex->ny));
498:       index[2] = (HYPRE_Int)(ex->zs + (row/(ex->nxny)));
499:       PetscStackCallStandard(HYPRE_SStructMatrixSetValues,(ex->ss_mat,part,index,var_type,7*nvars,entries,values[var_type]));
500:     }
501:   }
502:   PetscStackCallStandard(HYPRE_SStructMatrixAssemble,(ex->ss_mat));
503:   PetscFree(values[0]);
504:   PetscFree(values);
505:   PetscFree(entries);
506:   return(0);
507: }

509: PetscErrorCode MatZeroEntries_HYPRESStruct_3d(Mat mat)
510: {
511:   PetscErrorCode   ierr;
512:   Mat_HYPRESStruct *ex  = (Mat_HYPRESStruct*) mat->data;
513:   PetscInt         nvars= ex->nvars;
514:   PetscInt         size;
515:   PetscInt         part= 0;   /* only one part */

518:   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);
519:   {
520:     HYPRE_Int     i,*entries,iupper[3],ilower[3];
521:     HYPRE_Complex *values;

523:     for (i= 0; i< 3; i++) {
524:       ilower[i] = ex->hbox.imin[i];
525:       iupper[i] = ex->hbox.imax[i];
526:     }

528:     PetscMalloc2(nvars*7,&entries,nvars*7*size,&values);
529:     for (i= 0; i< nvars*7; i++) entries[i] = i;
530:     PetscArrayzero(values,nvars*7*size);

532:     for (i= 0; i< nvars; i++) {
533:       PetscStackCallStandard(HYPRE_SStructMatrixSetBoxValues,(ex->ss_mat,part,ilower,iupper,i,nvars*7,entries,values));
534:     }
535:     PetscFree2(entries,values);
536:   }
537:   PetscStackCallStandard(HYPRE_SStructMatrixAssemble,(ex->ss_mat));
538:   return(0);
539: }

541: static PetscErrorCode  MatSetUp_HYPRESStruct(Mat mat)
542: {
543:   PetscErrorCode         ierr;
544:   Mat_HYPRESStruct       *ex = (Mat_HYPRESStruct*) mat->data;
545:   PetscInt               dim,dof,sw[3],nx,ny,nz;
546:   PetscInt               ilower[3],iupper[3],ssize,i;
547:   DMBoundaryType         px,py,pz;
548:   DMDAStencilType        st;
549:   PetscInt               nparts= 1;  /* assuming only one part */
550:   PetscInt               part  = 0;
551:   ISLocalToGlobalMapping ltog;
552:   DM                     da;

555:   MatGetDM(mat,(DM*)&da);
556:   ex->da = da;
557:   PetscObjectReference((PetscObject)da);

559:   DMDAGetInfo(ex->da,&dim,0,0,0,0,0,0,&dof,&sw[0],&px,&py,&pz,&st);
560:   DMDAGetCorners(ex->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);
561:   iupper[0] += ilower[0] - 1;
562:   iupper[1] += ilower[1] - 1;
563:   iupper[2] += ilower[2] - 1;
564:   /* the hypre_Box is used to zero out the matrix entries in MatZeroValues() */
565:   ex->hbox.imin[0] = ilower[0];
566:   ex->hbox.imin[1] = ilower[1];
567:   ex->hbox.imin[2] = ilower[2];
568:   ex->hbox.imax[0] = iupper[0];
569:   ex->hbox.imax[1] = iupper[1];
570:   ex->hbox.imax[2] = iupper[2];

572:   ex->dofs_order = 0;

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

577:   /* create the hypre grid object and set its information */
578:   if (px || py || pz) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Ask us to add periodic support by calling HYPRE_SStructGridSetPeriodic()");
579:   PetscStackCallStandard(HYPRE_SStructGridCreate,(ex->hcomm,dim,nparts,&ex->ss_grid));
580:   PetscStackCallStandard(HYPRE_SStructGridSetExtents,(ex->ss_grid,part,ex->hbox.imin,ex->hbox.imax));
581:   {
582:     HYPRE_SStructVariable *vartypes;
583:     PetscMalloc1(ex->nvars,&vartypes);
584:     for (i= 0; i< ex->nvars; i++) vartypes[i]= HYPRE_SSTRUCT_VARIABLE_CELL;
585:     PetscStackCallStandard(HYPRE_SStructGridSetVariables,(ex->ss_grid, part, ex->nvars,vartypes));
586:     PetscFree(vartypes);
587:   }
588:   PetscStackCallStandard(HYPRE_SStructGridAssemble,(ex->ss_grid));

590:   sw[1] = sw[0];
591:   sw[2] = sw[1];
592:   /* PetscStackCallStandard(HYPRE_SStructGridSetNumGhost,(ex->ss_grid,sw)); */

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

598:   if (dim == 1) {
599:     HYPRE_Int offsets[3][1] = {{-1},{0},{1}};
600:     PetscInt  j, cnt;

602:     ssize = 3*(ex->nvars);
603:     PetscStackCallStandard(HYPRE_SStructStencilCreate,(dim,ssize,&ex->ss_stencil));
604:     cnt= 0;
605:     for (i = 0; i < (ex->nvars); i++) {
606:       for (j = 0; j < 3; j++) {
607:         PetscStackCallStandard(HYPRE_SStructStencilSetEntry,(ex->ss_stencil, cnt, offsets[j], i));
608:         cnt++;
609:       }
610:     }
611:   } else if (dim == 2) {
612:     HYPRE_Int offsets[5][2] = {{0,-1},{-1,0},{0,0},{1,0},{0,1}};
613:     PetscInt  j, cnt;

615:     ssize = 5*(ex->nvars);
616:     PetscStackCallStandard(HYPRE_SStructStencilCreate,(dim,ssize,&ex->ss_stencil));
617:     cnt= 0;
618:     for (i= 0; i< (ex->nvars); i++) {
619:       for (j= 0; j< 5; j++) {
620:         PetscStackCallStandard(HYPRE_SStructStencilSetEntry,(ex->ss_stencil, cnt, offsets[j], i));
621:         cnt++;
622:       }
623:     }
624:   } else if (dim == 3) {
625:     HYPRE_Int offsets[7][3] = {{0,0,-1},{0,-1,0},{-1,0,0},{0,0,0},{1,0,0},{0,1,0},{0,0,1}};
626:     PetscInt  j, cnt;

628:     ssize = 7*(ex->nvars);
629:     PetscStackCallStandard(HYPRE_SStructStencilCreate,(dim,ssize,&ex->ss_stencil));
630:     cnt= 0;
631:     for (i= 0; i< (ex->nvars); i++) {
632:       for (j= 0; j< 7; j++) {
633:         PetscStackCallStandard(HYPRE_SStructStencilSetEntry,(ex->ss_stencil, cnt, offsets[j], i));
634:         cnt++;
635:       }
636:     }
637:   }

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

642:   /* set the stencil graph. Note that each variable has the same graph. This means that each
643:      variable couples to all the other variable and with the same stencil pattern. */
644:   for (i= 0; i< (ex->nvars); i++) {
645:     PetscStackCallStandard(HYPRE_SStructGraphSetStencil,(ex->ss_graph,part,i,ex->ss_stencil));
646:   }
647:   PetscStackCallStandard(HYPRE_SStructGraphAssemble,(ex->ss_graph));

649:   /* create the HYPRE sstruct vectors for rhs and solution */
650:   PetscStackCallStandard(HYPRE_SStructVectorCreate,(ex->hcomm,ex->ss_grid,&ex->ss_b));
651:   PetscStackCallStandard(HYPRE_SStructVectorCreate,(ex->hcomm,ex->ss_grid,&ex->ss_x));
652:   PetscStackCallStandard(HYPRE_SStructVectorInitialize,(ex->ss_b));
653:   PetscStackCallStandard(HYPRE_SStructVectorInitialize,(ex->ss_x));
654:   PetscStackCallStandard(HYPRE_SStructVectorAssemble,(ex->ss_b));
655:   PetscStackCallStandard(HYPRE_SStructVectorAssemble,(ex->ss_x));

657:   /* create the hypre matrix object and set its information */
658:   PetscStackCallStandard(HYPRE_SStructMatrixCreate,(ex->hcomm,ex->ss_graph,&ex->ss_mat));
659:   PetscStackCallStandard(HYPRE_SStructGridDestroy,(ex->ss_grid));
660:   PetscStackCallStandard(HYPRE_SStructStencilDestroy,(ex->ss_stencil));
661:   if (ex->needsinitialization) {
662:     PetscStackCallStandard(HYPRE_SStructMatrixInitialize,(ex->ss_mat));
663:     ex->needsinitialization = PETSC_FALSE;
664:   }

666:   /* set the global and local sizes of the matrix */
667:   DMDAGetCorners(da,0,0,0,&nx,&ny,&nz);
668:   MatSetSizes(mat,dof*nx*ny*nz,dof*nx*ny*nz,PETSC_DECIDE,PETSC_DECIDE);
669:   PetscLayoutSetBlockSize(mat->rmap,dof);
670:   PetscLayoutSetBlockSize(mat->cmap,dof);
671:   PetscLayoutSetUp(mat->rmap);
672:   PetscLayoutSetUp(mat->cmap);
673:   mat->preallocated = PETSC_TRUE;

675:   if (dim == 3) {
676:     mat->ops->setvalueslocal = MatSetValuesLocal_HYPRESStruct_3d;
677:     mat->ops->zerorowslocal  = MatZeroRowsLocal_HYPRESStruct_3d;
678:     mat->ops->zeroentries    = MatZeroEntries_HYPRESStruct_3d;

680:     /* MatZeroEntries_HYPRESStruct_3d(mat); */
681:   } else SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Only support for 3d DMDA currently");

683:   /* get values that will be used repeatedly in MatSetValuesLocal() and MatZeroRowsLocal() repeatedly */
684:   MatGetOwnershipRange(mat,&ex->rstart,NULL);
685:   DMGetLocalToGlobalMapping(ex->da,&ltog);
686:   ISLocalToGlobalMappingGetIndices(ltog, (const PetscInt **) &ex->gindices);
687:   DMDAGetGhostCorners(ex->da,0,0,0,&ex->gnx,&ex->gnxgny,&ex->gnxgnygnz);

689:   ex->gnxgny    *= ex->gnx;
690:   ex->gnxgnygnz *= ex->gnxgny;

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

694:   ex->nxny   = ex->nx*ex->ny;
695:   ex->nxnynz = ex->nz*ex->nxny;
696:   return(0);
697: }

699: PetscErrorCode MatMult_HYPRESStruct(Mat A,Vec x,Vec y)
700: {
701:   PetscErrorCode    ierr;
702:   const PetscScalar *xx;
703:   PetscScalar       *yy;
704:   HYPRE_Int         hlower[3],hupper[3];
705:   PetscInt          ilower[3],iupper[3];
706:   Mat_HYPRESStruct  *mx     = (Mat_HYPRESStruct*)(A->data);
707:   PetscInt          ordering= mx->dofs_order;
708:   PetscInt          nvars   = mx->nvars;
709:   PetscInt          part    = 0;
710:   PetscInt          size;
711:   PetscInt          i;

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

716:   /* when HYPRE_MIXEDINT is defined, sizeof(HYPRE_Int) == 32 */
717:   iupper[0] += ilower[0] - 1;
718:   iupper[1] += ilower[1] - 1;
719:   iupper[2] += ilower[2] - 1;
720:   hlower[0]  = (HYPRE_Int)ilower[0];
721:   hlower[1]  = (HYPRE_Int)ilower[1];
722:   hlower[2]  = (HYPRE_Int)ilower[2];
723:   hupper[0]  = (HYPRE_Int)iupper[0];
724:   hupper[1]  = (HYPRE_Int)iupper[1];
725:   hupper[2]  = (HYPRE_Int)iupper[2];

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

730:   /* copy x values over to hypre for variable ordering */
731:   if (ordering) {
732:     PetscStackCallStandard(HYPRE_SStructVectorSetConstantValues,(mx->ss_b,0.0));
733:     VecGetArrayRead(x,&xx);
734:     for (i= 0; i< nvars; i++) {
735:       PetscStackCallStandard(HYPRE_SStructVectorSetBoxValues,(mx->ss_b,part,hlower,hupper,i,(HYPRE_Complex*)(xx+(size*i))));
736:     }
737:     VecRestoreArrayRead(x,&xx);
738:     PetscStackCallStandard(HYPRE_SStructVectorAssemble,(mx->ss_b));
739:     PetscStackCallStandard(HYPRE_SStructMatrixMatvec,(1.0,mx->ss_mat,mx->ss_b,0.0,mx->ss_x));

741:     /* copy solution values back to PETSc */
742:     VecGetArray(y,&yy);
743:     for (i= 0; i< nvars; i++) {
744:       PetscStackCallStandard(HYPRE_SStructVectorGetBoxValues,(mx->ss_x,part,hlower,hupper,i,(HYPRE_Complex*)(yy+(size*i))));
745:     }
746:     VecRestoreArray(y,&yy);
747:   } else {      /* nodal ordering must be mapped to variable ordering for sys_pfmg */
748:     PetscScalar *z;
749:     PetscInt    j, k;

751:     PetscMalloc1(nvars*size,&z);
752:     PetscStackCallStandard(HYPRE_SStructVectorSetConstantValues,(mx->ss_b,0.0));
753:     VecGetArrayRead(x,&xx);

755:     /* transform nodal to hypre's variable ordering for sys_pfmg */
756:     for (i= 0; i< size; i++) {
757:       k= i*nvars;
758:       for (j= 0; j< nvars; j++) z[j*size+i]= xx[k+j];
759:     }
760:     for (i= 0; i< nvars; i++) {
761:       PetscStackCallStandard(HYPRE_SStructVectorSetBoxValues,(mx->ss_b,part,hlower,hupper,i,(HYPRE_Complex*)(z+(size*i))));
762:     }
763:     VecRestoreArrayRead(x,&xx);
764:     PetscStackCallStandard(HYPRE_SStructVectorAssemble,(mx->ss_b));
765:     PetscStackCallStandard(HYPRE_SStructMatrixMatvec,(1.0,mx->ss_mat,mx->ss_b,0.0,mx->ss_x));

767:     /* copy solution values back to PETSc */
768:     VecGetArray(y,&yy);
769:     for (i= 0; i< nvars; i++) {
770:       PetscStackCallStandard(HYPRE_SStructVectorGetBoxValues,(mx->ss_x,part,hlower,hupper,i,(HYPRE_Complex*)(z+(size*i))));
771:     }
772:     /* transform hypre's variable ordering for sys_pfmg to nodal ordering */
773:     for (i= 0; i< size; i++) {
774:       k= i*nvars;
775:       for (j= 0; j< nvars; j++) yy[k+j]= z[j*size+i];
776:     }
777:     VecRestoreArray(y,&yy);
778:     PetscFree(z);
779:   }
780:   return(0);
781: }

783: PetscErrorCode MatAssemblyEnd_HYPRESStruct(Mat mat,MatAssemblyType mode)
784: {
785:   Mat_HYPRESStruct *ex = (Mat_HYPRESStruct*) mat->data;

788:   PetscStackCallStandard(HYPRE_SStructMatrixAssemble,(ex->ss_mat));
789:   return(0);
790: }

792: PetscErrorCode MatZeroEntries_HYPRESStruct(Mat mat)
793: {
795:   /* before the DMDA is set to the matrix the zero doesn't need to do anything */
796:   return(0);
797: }

799: PetscErrorCode MatDestroy_HYPRESStruct(Mat mat)
800: {
801:   Mat_HYPRESStruct       *ex = (Mat_HYPRESStruct*) mat->data;
802:   PetscErrorCode         ierr;
803:   ISLocalToGlobalMapping ltog;

806:   DMGetLocalToGlobalMapping(ex->da,&ltog);
807:   ISLocalToGlobalMappingRestoreIndices(ltog, (const PetscInt **) &ex->gindices);
808:   PetscStackCallStandard(HYPRE_SStructGraphDestroy,(ex->ss_graph));
809:   PetscStackCallStandard(HYPRE_SStructMatrixDestroy,(ex->ss_mat));
810:   PetscStackCallStandard(HYPRE_SStructVectorDestroy,(ex->ss_x));
811:   PetscStackCallStandard(HYPRE_SStructVectorDestroy,(ex->ss_b));
812:   PetscObjectDereference((PetscObject)ex->da);
813:   MPI_Comm_free(&(ex->hcomm));
814:   PetscFree(ex);
815:   return(0);
816: }

818: PETSC_EXTERN PetscErrorCode MatCreate_HYPRESStruct(Mat B)
819: {
820:   Mat_HYPRESStruct *ex;
821:   PetscErrorCode   ierr;

824:   PetscNewLog(B,&ex);
825:   B->data      = (void*)ex;
826:   B->rmap->bs  = 1;
827:   B->assembled = PETSC_FALSE;

829:   B->insertmode = NOT_SET_VALUES;

831:   B->ops->assemblyend = MatAssemblyEnd_HYPRESStruct;
832:   B->ops->mult        = MatMult_HYPRESStruct;
833:   B->ops->zeroentries = MatZeroEntries_HYPRESStruct;
834:   B->ops->destroy     = MatDestroy_HYPRESStruct;
835:   B->ops->setup       = MatSetUp_HYPRESStruct;

837:   ex->needsinitialization = PETSC_TRUE;

839:   MPI_Comm_dup(PetscObjectComm((PetscObject)B),&(ex->hcomm));
840:   PetscObjectChangeTypeName((PetscObject)B,MATHYPRESSTRUCT);
841:   return(0);
842: }