Actual source code: mhyp.c

petsc-master 2018-06-16
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:   PetscErrorCode    ierr;
 31:   PetscInt          i,j,stencil,index[3],row,entries[7];
 32:   const PetscScalar *values = 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] = ex->xs + (row % ex->nx);
 58:     index[1] = ex->ys + ((row/ex->nx) % ex->ny);
 59:     index[2] = ex->zs + (row/(ex->nxny));
 60:     if (addv == ADD_VALUES) {
 61:       PetscStackCallStandard(HYPRE_StructMatrixAddToValues,(ex->hmat,(HYPRE_Int *)index,ncol,(HYPRE_Int *)entries,(PetscScalar*)values));
 62:     } else {
 63:       PetscStackCallStandard(HYPRE_StructMatrixSetValues,(ex->hmat,(HYPRE_Int *)index,ncol,(HYPRE_Int *)entries,(PetscScalar*)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:   PetscInt        i,index[3],row,entries[7] = {0,1,2,3,4,5,6};
 74:   PetscScalar     values[7];
 75:   Mat_HYPREStruct *ex = (Mat_HYPREStruct*) mat->data;

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

 92: PetscErrorCode MatZeroEntries_HYPREStruct_3d(Mat mat)
 93: {
 94:   PetscErrorCode  ierr;
 95:   PetscInt        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,(HYPRE_Int *)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:   PetscInt               dim,dof,sw[6],nx,ny,nz,ilower[3],iupper[3],ssize,i;
110:   DMBoundaryType         px,py,pz;
111:   DMDAStencilType        st;
112:   ISLocalToGlobalMapping ltog;
113:   DM                     da;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

226:   iupper[0] += ilower[0] - 1;
227:   iupper[1] += ilower[1] - 1;
228:   iupper[2] += ilower[2] - 1;

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

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

245: PetscErrorCode MatAssemblyEnd_HYPREStruct(Mat mat,MatAssemblyType mode)
246: {
247:   Mat_HYPREStruct *ex = (Mat_HYPREStruct*) mat->data;
248:   PetscErrorCode  ierr;

251:   PetscStackCallStandard(HYPRE_StructMatrixAssemble,(ex->hmat));
252:   /* PetscStackCallStandard(HYPRE_StructMatrixPrint,("dummy",ex->hmat,0)); */
253:   return(0);
254: }

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

263: PetscErrorCode MatDestroy_HYPREStruct(Mat mat)
264: {
265:   Mat_HYPREStruct *ex = (Mat_HYPREStruct*) mat->data;
266:   PetscErrorCode  ierr;

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

278: PETSC_EXTERN PetscErrorCode MatCreate_HYPREStruct(Mat B)
279: {
280:   Mat_HYPREStruct *ex;
281:   PetscErrorCode  ierr;

284:   PetscNewLog(B,&ex);
285:   B->data      = (void*)ex;
286:   B->rmap->bs  = 1;
287:   B->assembled = PETSC_FALSE;

289:   B->insertmode = NOT_SET_VALUES;

291:   B->ops->assemblyend = MatAssemblyEnd_HYPREStruct;
292:   B->ops->mult        = MatMult_HYPREStruct;
293:   B->ops->zeroentries = MatZeroEntries_HYPREStruct;
294:   B->ops->destroy     = MatDestroy_HYPREStruct;
295:   B->ops->setup       = MatSetUp_HYPREStruct;

297:   ex->needsinitialization = PETSC_TRUE;

299:   MPI_Comm_dup(PetscObjectComm((PetscObject)B),&(ex->hcomm));
300:   PetscObjectChangeTypeName((PetscObject)B,MATHYPRESTRUCT);
301:   return(0);
302: }

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


309:    Level: intermediate

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

315:           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
316:           be defined by a DMDA.

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

320: M*/

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

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

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

347:       for (j=0; j<ncol; j++) {
348:         to_grid_rank= icol[j]/nvars;
349:         to_var_type = (icol[j] % nvars);

351:         to_var_entry= to_var_entry*7;
352:         entries[j]  = to_var_entry;

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

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

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

389:       for (j=0; j<ncol; j++) {
390:         to_var_type = icol[j]/(ex->gnxgnygnz);
391:         to_grid_rank= icol[j] - to_var_type*(ex->gnxgnygnz);

393:         to_var_entry= to_var_entry*7;
394:         entries[j]  = to_var_entry;

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

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

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

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

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

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

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

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

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

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

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

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

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

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


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

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

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


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

540:   MatGetDM(mat,(DM*)&da);
541:   ex->da = da;
542:   PetscObjectReference((PetscObject)da);

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

557:   ex->dofs_order = 0;

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

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

575:   sw[1] = sw[0];
576:   sw[2] = sw[1];
577:   /* PetscStackCallStandard(HYPRE_SStructGridSetNumGhost,(ex->ss_grid,sw)); */

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

583:   if (dim == 1) {
584:     PetscInt offsets[3][1] = {{-1},{0},{1}};
585:     PetscInt j, cnt;

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

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

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

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

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

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

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

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

660:   if (dim == 3) {
661:     mat->ops->setvalueslocal = MatSetValuesLocal_HYPRESStruct_3d;
662:     mat->ops->zerorowslocal  = MatZeroRowsLocal_HYPRESStruct_3d;
663:     mat->ops->zeroentries    = MatZeroEntries_HYPRESStruct_3d;

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

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

674:   ex->gnxgny    *= ex->gnx;
675:   ex->gnxgnygnz *= ex->gnxgny;

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

679:   ex->nxny   = ex->nx*ex->ny;
680:   ex->nxnynz = ex->nz*ex->nxny;
681:   return(0);
682: }

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

698:   DMDAGetCorners(mx->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);
699:   iupper[0] += ilower[0] - 1;
700:   iupper[1] += ilower[1] - 1;
701:   iupper[2] += ilower[2] - 1;

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

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

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

727:     PetscMalloc1(nvars*size,&z);
728:     PetscStackCallStandard(HYPRE_SStructVectorSetConstantValues,(mx->ss_b,0.0));
729:     VecGetArrayRead(x,&xx);

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

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

759: PetscErrorCode MatAssemblyEnd_HYPRESStruct(Mat mat,MatAssemblyType mode)
760: {
761:   Mat_HYPRESStruct *ex = (Mat_HYPRESStruct*) mat->data;
762:   PetscErrorCode   ierr;

765:   PetscStackCallStandard(HYPRE_SStructMatrixAssemble,(ex->ss_mat));
766:   return(0);
767: }

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

776: PetscErrorCode MatDestroy_HYPRESStruct(Mat mat)
777: {
778:   Mat_HYPRESStruct       *ex = (Mat_HYPRESStruct*) mat->data;
779:   PetscErrorCode         ierr;
780:   ISLocalToGlobalMapping ltog;

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

795: PETSC_EXTERN PetscErrorCode MatCreate_HYPRESStruct(Mat B)
796: {
797:   Mat_HYPRESStruct *ex;
798:   PetscErrorCode   ierr;

801:   PetscNewLog(B,&ex);
802:   B->data      = (void*)ex;
803:   B->rmap->bs  = 1;
804:   B->assembled = PETSC_FALSE;

806:   B->insertmode = NOT_SET_VALUES;

808:   B->ops->assemblyend = MatAssemblyEnd_HYPRESStruct;
809:   B->ops->mult        = MatMult_HYPRESStruct;
810:   B->ops->zeroentries = MatZeroEntries_HYPRESStruct;
811:   B->ops->destroy     = MatDestroy_HYPRESStruct;
812:   B->ops->setup       = MatSetUp_HYPRESStruct;

814:   ex->needsinitialization = PETSC_TRUE;

816:   MPI_Comm_dup(PetscObjectComm((PetscObject)B),&(ex->hcomm));
817:   PetscObjectChangeTypeName((PetscObject)B,MATHYPRESSTRUCT);
818:   return(0);
819: }