Actual source code: mhyp.c

petsc-master 2019-05-22
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>
  8:  #include <../src/dm/impls/da/hypre/mhyp.h>

 10: /* -----------------------------------------------------------------------------------------------------------------*/

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

 16:    Level: intermediate

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

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

 24: .seealso: MatCreate(), PCPFMG, MatSetDM(), DMCreateMatrix()
 25: 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:   PetscInt          i,j,stencil,index[3],row,entries[7];
 31:   const PetscScalar *values = y;
 32:   Mat_HYPREStruct   *ex     = (Mat_HYPREStruct*) mat->data;

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

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

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

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

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

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

114:   MatGetDM(mat,(DM*)&da);
115:   ex->da = da;
116:   PetscObjectReference((PetscObject)da);

118:   DMDAGetInfo(ex->da,&dim,0,0,0,0,0,0,&dof,&sw[0],&px,&py,&pz,&st);
119:   DMDAGetCorners(ex->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);
120:   iupper[0] += ilower[0] - 1;
121:   iupper[1] += ilower[1] - 1;
122:   iupper[2] += ilower[2] - 1;

124:   /* the hypre_Box is used to zero out the matrix entries in MatZeroValues() */
125:   ex->hbox.imin[0] = ilower[0];
126:   ex->hbox.imin[1] = ilower[1];
127:   ex->hbox.imin[2] = ilower[2];
128:   ex->hbox.imax[0] = iupper[0];
129:   ex->hbox.imax[1] = iupper[1];
130:   ex->hbox.imax[2] = iupper[2];

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

139:   sw[5] = sw[4] = sw[3] = sw[2] = sw[1] = sw[0];
140:   PetscStackCallStandard(HYPRE_StructGridSetNumGhost,(ex->hgrid,(HYPRE_Int *)sw));

142:   /* create the hypre stencil object and set its information */
143:   if (sw[0] > 1) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Ask us to add support for wider stencils");
144:   if (st == DMDA_STENCIL_BOX) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Ask us to add support for box stencils");
145:   if (dim == 1) {
146:     PetscInt offsets[3][1] = {{-1},{0},{1}};
147:     ssize = 3;
148:     PetscStackCallStandard(HYPRE_StructStencilCreate,(dim,ssize,&ex->hstencil));
149:     for (i=0; i<ssize; i++) {
150:       PetscStackCallStandard(HYPRE_StructStencilSetElement,(ex->hstencil,i,(HYPRE_Int *)offsets[i]));
151:     }
152:   } else if (dim == 2) {
153:     PetscInt offsets[5][2] = {{0,-1},{-1,0},{0,0},{1,0},{0,1}};
154:     ssize = 5;
155:     PetscStackCallStandard(HYPRE_StructStencilCreate,(dim,ssize,&ex->hstencil));
156:     for (i=0; i<ssize; i++) {
157:       PetscStackCallStandard(HYPRE_StructStencilSetElement,(ex->hstencil,i,(HYPRE_Int *)offsets[i]));
158:     }
159:   } else if (dim == 3) {
160:     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}};
161:     ssize = 7;
162:     PetscStackCallStandard(HYPRE_StructStencilCreate,(dim,ssize,&ex->hstencil));
163:     for (i=0; i<ssize; i++) {
164:       PetscStackCallStandard(HYPRE_StructStencilSetElement,(ex->hstencil,i,(HYPRE_Int *)offsets[i]));
165:     }
166:   }

168:   /* create the HYPRE vector for rhs and solution */
169:   PetscStackCallStandard(HYPRE_StructVectorCreate,(ex->hcomm,ex->hgrid,&ex->hb));
170:   PetscStackCallStandard(HYPRE_StructVectorCreate,(ex->hcomm,ex->hgrid,&ex->hx));
171:   PetscStackCallStandard(HYPRE_StructVectorInitialize,(ex->hb));
172:   PetscStackCallStandard(HYPRE_StructVectorInitialize,(ex->hx));
173:   PetscStackCallStandard(HYPRE_StructVectorAssemble,(ex->hb));
174:   PetscStackCallStandard(HYPRE_StructVectorAssemble,(ex->hx));

176:   /* create the hypre matrix object and set its information */
177:   PetscStackCallStandard(HYPRE_StructMatrixCreate,(ex->hcomm,ex->hgrid,ex->hstencil,&ex->hmat));
178:   PetscStackCallStandard(HYPRE_StructGridDestroy,(ex->hgrid));
179:   PetscStackCallStandard(HYPRE_StructStencilDestroy,(ex->hstencil));
180:   if (ex->needsinitialization) {
181:     PetscStackCallStandard(HYPRE_StructMatrixInitialize,(ex->hmat));
182:     ex->needsinitialization = PETSC_FALSE;
183:   }

185:   /* set the global and local sizes of the matrix */
186:   DMDAGetCorners(da,0,0,0,&nx,&ny,&nz);
187:   MatSetSizes(mat,dof*nx*ny*nz,dof*nx*ny*nz,PETSC_DECIDE,PETSC_DECIDE);
188:   PetscLayoutSetBlockSize(mat->rmap,1);
189:   PetscLayoutSetBlockSize(mat->cmap,1);
190:   PetscLayoutSetUp(mat->rmap);
191:   PetscLayoutSetUp(mat->cmap);
192:   mat->preallocated = PETSC_TRUE;

194:   if (dim == 3) {
195:     mat->ops->setvalueslocal = MatSetValuesLocal_HYPREStruct_3d;
196:     mat->ops->zerorowslocal  = MatZeroRowsLocal_HYPREStruct_3d;
197:     mat->ops->zeroentries    = MatZeroEntries_HYPREStruct_3d;

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

202:   /* get values that will be used repeatedly in MatSetValuesLocal() and MatZeroRowsLocal() repeatedly */
203:   MatGetOwnershipRange(mat,&ex->rstart,NULL);
204:   DMGetLocalToGlobalMapping(ex->da,&ltog);
205:   ISLocalToGlobalMappingGetIndices(ltog, (const PetscInt **) &ex->gindices);
206:   DMDAGetGhostCorners(ex->da,0,0,0,&ex->gnx,&ex->gnxgny,0);
207:   ex->gnxgny *= ex->gnx;
208:   DMDAGetCorners(ex->da,&ex->xs,&ex->ys,&ex->zs,&ex->nx,&ex->ny,0);
209:   ex->nxny    = ex->nx*ex->ny;
210:   return(0);
211: }

213: PetscErrorCode MatMult_HYPREStruct(Mat A,Vec x,Vec y)
214: {
215:   PetscErrorCode    ierr;
216:   const PetscScalar *xx;
217:   PetscScalar       *yy;
218:   PetscInt          ilower[3],iupper[3];
219:   Mat_HYPREStruct   *mx = (Mat_HYPREStruct*)(A->data);

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

224:   iupper[0] += ilower[0] - 1;
225:   iupper[1] += ilower[1] - 1;
226:   iupper[2] += ilower[2] - 1;

228:   /* copy x values over to hypre */
229:   PetscStackCallStandard(HYPRE_StructVectorSetConstantValues,(mx->hb,0.0));
230:   VecGetArrayRead(x,&xx);
231:   PetscStackCallStandard(HYPRE_StructVectorSetBoxValues,(mx->hb,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,(PetscScalar*)xx));
232:   VecRestoreArrayRead(x,&xx);
233:   PetscStackCallStandard(HYPRE_StructVectorAssemble,(mx->hb));
234:   PetscStackCallStandard(HYPRE_StructMatrixMatvec,(1.0,mx->hmat,mx->hb,0.0,mx->hx));

236:   /* copy solution values back to PETSc */
237:   VecGetArray(y,&yy);
238:   PetscStackCallStandard(HYPRE_StructVectorGetBoxValues,(mx->hx,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,yy));
239:   VecRestoreArray(y,&yy);
240:   return(0);
241: }

243: PetscErrorCode MatAssemblyEnd_HYPREStruct(Mat mat,MatAssemblyType mode)
244: {
245:   Mat_HYPREStruct *ex = (Mat_HYPREStruct*) mat->data;

248:   PetscStackCallStandard(HYPRE_StructMatrixAssemble,(ex->hmat));
249:   /* PetscStackCallStandard(HYPRE_StructMatrixPrint,("dummy",ex->hmat,0)); */
250:   return(0);
251: }

253: PetscErrorCode MatZeroEntries_HYPREStruct(Mat mat)
254: {
256:   /* before the DMDA is set to the matrix the zero doesn't need to do anything */
257:   return(0);
258: }

260: PetscErrorCode MatDestroy_HYPREStruct(Mat mat)
261: {
262:   Mat_HYPREStruct *ex = (Mat_HYPREStruct*) mat->data;
263:   PetscErrorCode  ierr;

266:   PetscStackCallStandard(HYPRE_StructMatrixDestroy,(ex->hmat));
267:   PetscStackCallStandard(HYPRE_StructVectorDestroy,(ex->hx));
268:   PetscStackCallStandard(HYPRE_StructVectorDestroy,(ex->hb));
269:   PetscObjectDereference((PetscObject)ex->da);
270:   MPI_Comm_free(&(ex->hcomm));
271:   PetscFree(ex);
272:   return(0);
273: }

275: PETSC_EXTERN PetscErrorCode MatCreate_HYPREStruct(Mat B)
276: {
277:   Mat_HYPREStruct *ex;
278:   PetscErrorCode  ierr;

281:   PetscNewLog(B,&ex);
282:   B->data      = (void*)ex;
283:   B->rmap->bs  = 1;
284:   B->assembled = PETSC_FALSE;

286:   B->insertmode = NOT_SET_VALUES;

288:   B->ops->assemblyend = MatAssemblyEnd_HYPREStruct;
289:   B->ops->mult        = MatMult_HYPREStruct;
290:   B->ops->zeroentries = MatZeroEntries_HYPREStruct;
291:   B->ops->destroy     = MatDestroy_HYPREStruct;
292:   B->ops->setup       = MatSetUp_HYPREStruct;

294:   ex->needsinitialization = PETSC_TRUE;

296:   MPI_Comm_dup(PetscObjectComm((PetscObject)B),&(ex->hcomm));
297:   PetscObjectChangeTypeName((PetscObject)B,MATHYPRESTRUCT);
298:   return(0);
299: }

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


306:    Level: intermediate

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

312:           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
313:           be defined by a DMDA.

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

317: M*/

319: PetscErrorCode  MatSetValuesLocal_HYPRESStruct_3d(Mat mat,PetscInt nrow,const PetscInt irow[],PetscInt ncol,const PetscInt icol[],const PetscScalar y[],InsertMode addv)
320: {
321:   PetscErrorCode    ierr;
322:   PetscInt          i,j,stencil,index[3];
323:   const PetscScalar *values = y;
324:   Mat_HYPRESStruct  *ex     = (Mat_HYPRESStruct*) mat->data;

326:   PetscInt          part = 0;          /* Petsc sstruct interface only allows 1 part */
327:   PetscInt          ordering;
328:   PetscInt          grid_rank, to_grid_rank;
329:   PetscInt          var_type, to_var_type;
330:   PetscInt          to_var_entry = 0;
331:   PetscInt          nvars= ex->nvars;
332:   PetscInt          row,*entries;

335:   PetscMalloc1(7*nvars,&entries);
336:   ordering = ex->dofs_order;  /* ordering= 0   nodal ordering
337:                                            1   variable ordering */
338:   /* stencil entries are ordered by variables: var0_stencil0, var0_stencil1, ..., var0_stencil6, var1_stencil0, var1_stencil1, ...  */
339:   if (!ordering) {
340:     for (i=0; i<nrow; i++) {
341:       grid_rank= irow[i]/nvars;
342:       var_type = (irow[i] % nvars);

344:       for (j=0; j<ncol; j++) {
345:         to_grid_rank= icol[j]/nvars;
346:         to_var_type = (icol[j] % nvars);

348:         to_var_entry= to_var_entry*7;
349:         entries[j]  = to_var_entry;

351:         stencil = to_grid_rank-grid_rank;
352:         if (!stencil) {
353:           entries[j] += 3;
354:         } else if (stencil == -1) {
355:           entries[j] += 2;
356:         } else if (stencil == 1) {
357:           entries[j] += 4;
358:         } else if (stencil == -ex->gnx) {
359:           entries[j] += 1;
360:         } else if (stencil == ex->gnx) {
361:           entries[j] += 5;
362:         } else if (stencil == -ex->gnxgny) {
363:           entries[j] += 0;
364:         } else if (stencil == ex->gnxgny) {
365:           entries[j] += 6;
366:         } else SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Local row %D local column %D have bad stencil %D",irow[i],icol[j],stencil);
367:       }

369:       row      = ex->gindices[grid_rank] - ex->rstart;
370:       index[0] = ex->xs + (row % ex->nx);
371:       index[1] = ex->ys + ((row/ex->nx) % ex->ny);
372:       index[2] = ex->zs + (row/(ex->nxny));

374:       if (addv == ADD_VALUES) {
375:         PetscStackCallStandard(HYPRE_SStructMatrixAddToValues,(ex->ss_mat,part,(HYPRE_Int *)index,var_type,ncol,(HYPRE_Int *)entries,(PetscScalar*)values));
376:       } else {
377:         PetscStackCallStandard(HYPRE_SStructMatrixSetValues,(ex->ss_mat,part,(HYPRE_Int *)index,var_type,ncol,(HYPRE_Int *)entries,(PetscScalar*)values));
378:       }
379:       values += ncol;
380:     }
381:   } else {
382:     for (i=0; i<nrow; i++) {
383:       var_type = irow[i]/(ex->gnxgnygnz);
384:       grid_rank= irow[i] - var_type*(ex->gnxgnygnz);

386:       for (j=0; j<ncol; j++) {
387:         to_var_type = icol[j]/(ex->gnxgnygnz);
388:         to_grid_rank= icol[j] - to_var_type*(ex->gnxgnygnz);

390:         to_var_entry= to_var_entry*7;
391:         entries[j]  = to_var_entry;

393:         stencil = to_grid_rank-grid_rank;
394:         if (!stencil) {
395:           entries[j] += 3;
396:         } else if (stencil == -1) {
397:           entries[j] += 2;
398:         } else if (stencil == 1) {
399:           entries[j] += 4;
400:         } else if (stencil == -ex->gnx) {
401:           entries[j] += 1;
402:         } else if (stencil == ex->gnx) {
403:           entries[j] += 5;
404:         } else if (stencil == -ex->gnxgny) {
405:           entries[j] += 0;
406:         } else if (stencil == ex->gnxgny) {
407:           entries[j] += 6;
408:         } else SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Local row %D local column %D have bad stencil %D",irow[i],icol[j],stencil);
409:       }

411:       row      = ex->gindices[grid_rank] - ex->rstart;
412:       index[0] = ex->xs + (row % ex->nx);
413:       index[1] = ex->ys + ((row/ex->nx) % ex->ny);
414:       index[2] = ex->zs + (row/(ex->nxny));

416:       if (addv == ADD_VALUES) {
417:         PetscStackCallStandard(HYPRE_SStructMatrixAddToValues,(ex->ss_mat,part,(HYPRE_Int *)index,var_type,ncol,(HYPRE_Int *)entries,(PetscScalar*)values));
418:       } else {
419:         PetscStackCallStandard(HYPRE_SStructMatrixSetValues,(ex->ss_mat,part,(HYPRE_Int *)index,var_type,ncol,(HYPRE_Int *)entries,(PetscScalar*)values));
420:       }
421:       values += ncol;
422:     }
423:   }
424:   PetscFree(entries);
425:   return(0);
426: }

428: PetscErrorCode  MatZeroRowsLocal_HYPRESStruct_3d(Mat mat,PetscInt nrow,const PetscInt irow[],PetscScalar d,Vec x,Vec b)
429: {
430:   PetscErrorCode   ierr;
431:   PetscInt         i,index[3];
432:   PetscScalar      **values;
433:   Mat_HYPRESStruct *ex = (Mat_HYPRESStruct*) mat->data;

435:   PetscInt         part     = 0;      /* Petsc sstruct interface only allows 1 part */
436:   PetscInt         ordering = ex->dofs_order;
437:   PetscInt         grid_rank;
438:   PetscInt         var_type;
439:   PetscInt         nvars= ex->nvars;
440:   PetscInt         row,*entries;

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

446:   PetscMalloc1(nvars,&values);
447:   PetscMalloc1(7*nvars*nvars,&values[0]);
448:   for (i=1; i<nvars; i++) {
449:     values[i] = values[i-1] + nvars*7;
450:   }

452:   for (i=0; i< nvars; i++) {
453:     PetscMemzero(values[i],nvars*7*sizeof(PetscScalar));
454:     *(values[i]+3) = d;
455:   }

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

459:   if (!ordering) {
460:     for (i=0; i<nrow; i++) {
461:       grid_rank = irow[i] / nvars;
462:       var_type  =(irow[i] % nvars);

464:       row      = ex->gindices[grid_rank] - ex->rstart;
465:       index[0] = ex->xs + (row % ex->nx);
466:       index[1] = ex->ys + ((row/ex->nx) % ex->ny);
467:       index[2] = ex->zs + (row/(ex->nxny));
468:       PetscStackCallStandard(HYPRE_SStructMatrixSetValues,(ex->ss_mat,part,(HYPRE_Int *)index,var_type,7*nvars,(HYPRE_Int *)entries,values[var_type]));
469:     }
470:   } else {
471:     for (i=0; i<nrow; i++) {
472:       var_type = irow[i]/(ex->gnxgnygnz);
473:       grid_rank= irow[i] - var_type*(ex->gnxgnygnz);

475:       row      = ex->gindices[grid_rank] - ex->rstart;
476:       index[0] = ex->xs + (row % ex->nx);
477:       index[1] = ex->ys + ((row/ex->nx) % ex->ny);
478:       index[2] = ex->zs + (row/(ex->nxny));
479:       PetscStackCallStandard(HYPRE_SStructMatrixSetValues,(ex->ss_mat,part,(HYPRE_Int *)index,var_type,7*nvars,(HYPRE_Int *)entries,values[var_type]));
480:     }
481:   }
482:   PetscStackCallStandard(HYPRE_SStructMatrixAssemble,(ex->ss_mat));
483:   PetscFree(values[0]);
484:   PetscFree(values);
485:   PetscFree(entries);
486:   return(0);
487: }

489: PetscErrorCode MatZeroEntries_HYPRESStruct_3d(Mat mat)
490: {
491:   PetscErrorCode   ierr;
492:   Mat_HYPRESStruct *ex  = (Mat_HYPRESStruct*) mat->data;
493:   PetscInt         nvars= ex->nvars;
494:   PetscInt         size;
495:   PetscInt         part= 0;   /* only one part */

498:   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);
499:   {
500:     PetscInt    i,*entries,iupper[3],ilower[3];
501:     PetscScalar *values;


504:     for (i= 0; i< 3; i++) {
505:       ilower[i]= ex->hbox.imin[i];
506:       iupper[i]= ex->hbox.imax[i];
507:     }

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

513:     for (i= 0; i< nvars; i++) {
514:       PetscStackCallStandard(HYPRE_SStructMatrixSetBoxValues,(ex->ss_mat,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,nvars*7,(HYPRE_Int *)entries,values));
515:     }
516:     PetscFree2(entries,values);
517:   }
518:   PetscStackCallStandard(HYPRE_SStructMatrixAssemble,(ex->ss_mat));
519:   return(0);
520: }


523: static PetscErrorCode  MatSetUp_HYPRESStruct(Mat mat)
524: {
525:   PetscErrorCode         ierr;
526:   Mat_HYPRESStruct       *ex = (Mat_HYPRESStruct*) mat->data;
527:   PetscInt               dim,dof,sw[3],nx,ny,nz;
528:   PetscInt               ilower[3],iupper[3],ssize,i;
529:   DMBoundaryType         px,py,pz;
530:   DMDAStencilType        st;
531:   PetscInt               nparts= 1;  /* assuming only one part */
532:   PetscInt               part  = 0;
533:   ISLocalToGlobalMapping ltog;
534:   DM                     da;

537:   MatGetDM(mat,(DM*)&da);
538:   ex->da = da;
539:   PetscObjectReference((PetscObject)da);

541:   DMDAGetInfo(ex->da,&dim,0,0,0,0,0,0,&dof,&sw[0],&px,&py,&pz,&st);
542:   DMDAGetCorners(ex->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);
543:   iupper[0] += ilower[0] - 1;
544:   iupper[1] += ilower[1] - 1;
545:   iupper[2] += ilower[2] - 1;
546:   /* the hypre_Box is used to zero out the matrix entries in MatZeroValues() */
547:   ex->hbox.imin[0] = ilower[0];
548:   ex->hbox.imin[1] = ilower[1];
549:   ex->hbox.imin[2] = ilower[2];
550:   ex->hbox.imax[0] = iupper[0];
551:   ex->hbox.imax[1] = iupper[1];
552:   ex->hbox.imax[2] = iupper[2];

554:   ex->dofs_order = 0;

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

559:   /* create the hypre grid object and set its information */
560:   if (px || py || pz) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Ask us to add periodic support by calling HYPRE_SStructGridSetPeriodic()");
561:   PetscStackCallStandard(HYPRE_SStructGridCreate,(ex->hcomm,dim,nparts,&ex->ss_grid));
562:   PetscStackCallStandard(HYPRE_SStructGridSetExtents,(ex->ss_grid,part,ex->hbox.imin,ex->hbox.imax));
563:   {
564:     HYPRE_SStructVariable *vartypes;
565:     PetscMalloc1(ex->nvars,&vartypes);
566:     for (i= 0; i< ex->nvars; i++) vartypes[i]= HYPRE_SSTRUCT_VARIABLE_CELL;
567:     PetscStackCallStandard(HYPRE_SStructGridSetVariables,(ex->ss_grid, part, ex->nvars,vartypes));
568:     PetscFree(vartypes);
569:   }
570:   PetscStackCallStandard(HYPRE_SStructGridAssemble,(ex->ss_grid));

572:   sw[1] = sw[0];
573:   sw[2] = sw[1];
574:   /* PetscStackCallStandard(HYPRE_SStructGridSetNumGhost,(ex->ss_grid,sw)); */

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

580:   if (dim == 1) {
581:     PetscInt offsets[3][1] = {{-1},{0},{1}};
582:     PetscInt j, cnt;

584:     ssize = 3*(ex->nvars);
585:     PetscStackCallStandard(HYPRE_SStructStencilCreate,(dim,ssize,&ex->ss_stencil));
586:     cnt= 0;
587:     for (i = 0; i < (ex->nvars); i++) {
588:       for (j = 0; j < 3; j++) {
589:         PetscStackCallStandard(HYPRE_SStructStencilSetEntry,(ex->ss_stencil, cnt, (HYPRE_Int *)offsets[j], i));
590:         cnt++;
591:       }
592:     }
593:   } else if (dim == 2) {
594:     PetscInt offsets[5][2] = {{0,-1},{-1,0},{0,0},{1,0},{0,1}};
595:     PetscInt j, cnt;

597:     ssize = 5*(ex->nvars);
598:     PetscStackCallStandard(HYPRE_SStructStencilCreate,(dim,ssize,&ex->ss_stencil));
599:     cnt= 0;
600:     for (i= 0; i< (ex->nvars); i++) {
601:       for (j= 0; j< 5; j++) {
602:         PetscStackCallStandard(HYPRE_SStructStencilSetEntry,(ex->ss_stencil, cnt, (HYPRE_Int *)offsets[j], i));
603:         cnt++;
604:       }
605:     }
606:   } else if (dim == 3) {
607:     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}};
608:     PetscInt j, cnt;

610:     ssize = 7*(ex->nvars);
611:     PetscStackCallStandard(HYPRE_SStructStencilCreate,(dim,ssize,&ex->ss_stencil));
612:     cnt= 0;
613:     for (i= 0; i< (ex->nvars); i++) {
614:       for (j= 0; j< 7; j++) {
615:         PetscStackCallStandard(HYPRE_SStructStencilSetEntry,(ex->ss_stencil, cnt, (HYPRE_Int *)offsets[j], i));
616:         cnt++;
617:       }
618:     }
619:   }

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

624:   /* set the stencil graph. Note that each variable has the same graph. This means that each
625:      variable couples to all the other variable and with the same stencil pattern. */
626:   for (i= 0; i< (ex->nvars); i++) {
627:     PetscStackCallStandard(HYPRE_SStructGraphSetStencil,(ex->ss_graph,part,i,ex->ss_stencil));
628:   }
629:   PetscStackCallStandard(HYPRE_SStructGraphAssemble,(ex->ss_graph));

631:   /* create the HYPRE sstruct vectors for rhs and solution */
632:   PetscStackCallStandard(HYPRE_SStructVectorCreate,(ex->hcomm,ex->ss_grid,&ex->ss_b));
633:   PetscStackCallStandard(HYPRE_SStructVectorCreate,(ex->hcomm,ex->ss_grid,&ex->ss_x));
634:   PetscStackCallStandard(HYPRE_SStructVectorInitialize,(ex->ss_b));
635:   PetscStackCallStandard(HYPRE_SStructVectorInitialize,(ex->ss_x));
636:   PetscStackCallStandard(HYPRE_SStructVectorAssemble,(ex->ss_b));
637:   PetscStackCallStandard(HYPRE_SStructVectorAssemble,(ex->ss_x));

639:   /* create the hypre matrix object and set its information */
640:   PetscStackCallStandard(HYPRE_SStructMatrixCreate,(ex->hcomm,ex->ss_graph,&ex->ss_mat));
641:   PetscStackCallStandard(HYPRE_SStructGridDestroy,(ex->ss_grid));
642:   PetscStackCallStandard(HYPRE_SStructStencilDestroy,(ex->ss_stencil));
643:   if (ex->needsinitialization) {
644:     PetscStackCallStandard(HYPRE_SStructMatrixInitialize,(ex->ss_mat));
645:     ex->needsinitialization = PETSC_FALSE;
646:   }

648:   /* set the global and local sizes of the matrix */
649:   DMDAGetCorners(da,0,0,0,&nx,&ny,&nz);
650:   MatSetSizes(mat,dof*nx*ny*nz,dof*nx*ny*nz,PETSC_DECIDE,PETSC_DECIDE);
651:   PetscLayoutSetBlockSize(mat->rmap,dof);
652:   PetscLayoutSetBlockSize(mat->cmap,dof);
653:   PetscLayoutSetUp(mat->rmap);
654:   PetscLayoutSetUp(mat->cmap);
655:   mat->preallocated = PETSC_TRUE;

657:   if (dim == 3) {
658:     mat->ops->setvalueslocal = MatSetValuesLocal_HYPRESStruct_3d;
659:     mat->ops->zerorowslocal  = MatZeroRowsLocal_HYPRESStruct_3d;
660:     mat->ops->zeroentries    = MatZeroEntries_HYPRESStruct_3d;

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

665:   /* get values that will be used repeatedly in MatSetValuesLocal() and MatZeroRowsLocal() repeatedly */
666:   MatGetOwnershipRange(mat,&ex->rstart,NULL);
667:   DMGetLocalToGlobalMapping(ex->da,&ltog);
668:   ISLocalToGlobalMappingGetIndices(ltog, (const PetscInt **) &ex->gindices);
669:   DMDAGetGhostCorners(ex->da,0,0,0,&ex->gnx,&ex->gnxgny,&ex->gnxgnygnz);

671:   ex->gnxgny    *= ex->gnx;
672:   ex->gnxgnygnz *= ex->gnxgny;

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

676:   ex->nxny   = ex->nx*ex->ny;
677:   ex->nxnynz = ex->nz*ex->nxny;
678:   return(0);
679: }

681: PetscErrorCode MatMult_HYPRESStruct(Mat A,Vec x,Vec y)
682: {
683:   PetscErrorCode    ierr;
684:   const PetscScalar *xx;
685:   PetscScalar       *yy;
686:   PetscInt          ilower[3],iupper[3];
687:   Mat_HYPRESStruct  *mx     = (Mat_HYPRESStruct*)(A->data);
688:   PetscInt          ordering= mx->dofs_order;
689:   PetscInt          nvars   = mx->nvars;
690:   PetscInt          part    = 0;
691:   PetscInt          size;
692:   PetscInt          i;

695:   DMDAGetCorners(mx->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);
696:   iupper[0] += ilower[0] - 1;
697:   iupper[1] += ilower[1] - 1;
698:   iupper[2] += ilower[2] - 1;

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

703:   /* copy x values over to hypre for variable ordering */
704:   if (ordering) {
705:     PetscStackCallStandard(HYPRE_SStructVectorSetConstantValues,(mx->ss_b,0.0));
706:     VecGetArrayRead(x,&xx);
707:     for (i= 0; i< nvars; i++) {
708:       PetscStackCallStandard(HYPRE_SStructVectorSetBoxValues,(mx->ss_b,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,(PetscScalar*)xx+(size*i)));
709:     }
710:     VecRestoreArrayRead(x,&xx);
711:     PetscStackCallStandard(HYPRE_SStructVectorAssemble,(mx->ss_b));
712:     PetscStackCallStandard(HYPRE_SStructMatrixMatvec,(1.0,mx->ss_mat,mx->ss_b,0.0,mx->ss_x));

714:     /* copy solution values back to PETSc */
715:     VecGetArray(y,&yy);
716:     for (i= 0; i< nvars; i++) {
717:       PetscStackCallStandard(HYPRE_SStructVectorGetBoxValues,(mx->ss_x,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,yy+(size*i)));
718:     }
719:     VecRestoreArray(y,&yy);
720:   } else {      /* nodal ordering must be mapped to variable ordering for sys_pfmg */
721:     PetscScalar *z;
722:     PetscInt    j, k;

724:     PetscMalloc1(nvars*size,&z);
725:     PetscStackCallStandard(HYPRE_SStructVectorSetConstantValues,(mx->ss_b,0.0));
726:     VecGetArrayRead(x,&xx);

728:     /* transform nodal to hypre's variable ordering for sys_pfmg */
729:     for (i= 0; i< size; i++) {
730:       k= i*nvars;
731:       for (j= 0; j< nvars; j++) z[j*size+i]= xx[k+j];
732:     }
733:     for (i= 0; i< nvars; i++) {
734:       PetscStackCallStandard(HYPRE_SStructVectorSetBoxValues,(mx->ss_b,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,z+(size*i)));
735:     }
736:     VecRestoreArrayRead(x,&xx);
737:     PetscStackCallStandard(HYPRE_SStructVectorAssemble,(mx->ss_b));
738:     PetscStackCallStandard(HYPRE_SStructMatrixMatvec,(1.0,mx->ss_mat,mx->ss_b,0.0,mx->ss_x));

740:     /* copy solution values back to PETSc */
741:     VecGetArray(y,&yy);
742:     for (i= 0; i< nvars; i++) {
743:       PetscStackCallStandard(HYPRE_SStructVectorGetBoxValues,(mx->ss_x,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,z+(size*i)));
744:     }
745:     /* transform hypre's variable ordering for sys_pfmg to nodal ordering */
746:     for (i= 0; i< size; i++) {
747:       k= i*nvars;
748:       for (j= 0; j< nvars; j++) yy[k+j]= z[j*size+i];
749:     }
750:     VecRestoreArray(y,&yy);
751:     PetscFree(z);
752:   }
753:   return(0);
754: }

756: PetscErrorCode MatAssemblyEnd_HYPRESStruct(Mat mat,MatAssemblyType mode)
757: {
758:   Mat_HYPRESStruct *ex = (Mat_HYPRESStruct*) mat->data;

761:   PetscStackCallStandard(HYPRE_SStructMatrixAssemble,(ex->ss_mat));
762:   return(0);
763: }

765: PetscErrorCode MatZeroEntries_HYPRESStruct(Mat mat)
766: {
768:   /* before the DMDA is set to the matrix the zero doesn't need to do anything */
769:   return(0);
770: }

772: PetscErrorCode MatDestroy_HYPRESStruct(Mat mat)
773: {
774:   Mat_HYPRESStruct       *ex = (Mat_HYPRESStruct*) mat->data;
775:   PetscErrorCode         ierr;
776:   ISLocalToGlobalMapping ltog;

779:   DMGetLocalToGlobalMapping(ex->da,&ltog);
780:   ISLocalToGlobalMappingRestoreIndices(ltog, (const PetscInt **) &ex->gindices);
781:   PetscStackCallStandard(HYPRE_SStructGraphDestroy,(ex->ss_graph));
782:   PetscStackCallStandard(HYPRE_SStructMatrixDestroy,(ex->ss_mat));
783:   PetscStackCallStandard(HYPRE_SStructVectorDestroy,(ex->ss_x));
784:   PetscStackCallStandard(HYPRE_SStructVectorDestroy,(ex->ss_b));
785:   PetscObjectDereference((PetscObject)ex->da);
786:   MPI_Comm_free(&(ex->hcomm));
787:   PetscFree(ex);
788:   return(0);
789: }

791: PETSC_EXTERN PetscErrorCode MatCreate_HYPRESStruct(Mat B)
792: {
793:   Mat_HYPRESStruct *ex;
794:   PetscErrorCode   ierr;

797:   PetscNewLog(B,&ex);
798:   B->data      = (void*)ex;
799:   B->rmap->bs  = 1;
800:   B->assembled = PETSC_FALSE;

802:   B->insertmode = NOT_SET_VALUES;

804:   B->ops->assemblyend = MatAssemblyEnd_HYPRESStruct;
805:   B->ops->mult        = MatMult_HYPRESStruct;
806:   B->ops->zeroentries = MatZeroEntries_HYPRESStruct;
807:   B->ops->destroy     = MatDestroy_HYPRESStruct;
808:   B->ops->setup       = MatSetUp_HYPRESStruct;

810:   ex->needsinitialization = PETSC_TRUE;

812:   MPI_Comm_dup(PetscObjectComm((PetscObject)B),&(ex->hcomm));
813:   PetscObjectChangeTypeName((PetscObject)B,MATHYPRESSTRUCT);
814:   return(0);
815: }