Actual source code: mhyp.c

petsc-master 2016-12-07
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: Unlike the more general support for blocks in hypre this allows only one block per process and requires the block
 19:           be defined by a DMDA.

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

 23: .seealso: MatCreate(), PCPFMG, MatSetupDM(), DMCreateMatrix()
 24: M*/


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

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

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

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

 97: PetscErrorCode MatZeroEntries_HYPREStruct_3d(Mat mat)
 98: {
 99:   PetscErrorCode  ierr;
100:   PetscInt        indices[7] = {0,1,2,3,4,5,6};
101:   Mat_HYPREStruct *ex        = (Mat_HYPREStruct*) mat->data;

104:   /* hypre has no public interface to do this */
105:   PetscStackCallStandard(hypre_StructMatrixClearBoxValues,(ex->hmat,&ex->hbox,7,(HYPRE_Int *)indices,0,1));
106:   PetscStackCallStandard(HYPRE_StructMatrixAssemble,(ex->hmat));
107:   return(0);
108: }

112: static PetscErrorCode  MatSetupDM_HYPREStruct(Mat mat,DM da)
113: {
114:   PetscErrorCode         ierr;
115:   Mat_HYPREStruct        *ex = (Mat_HYPREStruct*) mat->data;
116:   PetscInt               dim,dof,sw[3],nx,ny,nz,ilower[3],iupper[3],ssize,i;
117:   DMBoundaryType         px,py,pz;
118:   DMDAStencilType        st;
119:   ISLocalToGlobalMapping ltog;

122:   ex->da = da;
123:   PetscObjectReference((PetscObject)da);

125:   DMDAGetInfo(ex->da,&dim,0,0,0,0,0,0,&dof,&sw[0],&px,&py,&pz,&st);
126:   DMDAGetCorners(ex->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);
127:   iupper[0] += ilower[0] - 1;
128:   iupper[1] += ilower[1] - 1;
129:   iupper[2] += ilower[2] - 1;

131:   /* the hypre_Box is used to zero out the matrix entries in MatZeroValues() */
132:   ex->hbox.imin[0] = ilower[0];
133:   ex->hbox.imin[1] = ilower[1];
134:   ex->hbox.imin[2] = ilower[2];
135:   ex->hbox.imax[0] = iupper[0];
136:   ex->hbox.imax[1] = iupper[1];
137:   ex->hbox.imax[2] = iupper[2];

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

146:   sw[1] = sw[0];
147:   sw[2] = sw[1];
148:   PetscStackCallStandard(HYPRE_StructGridSetNumGhost,(ex->hgrid,(HYPRE_Int *)sw));

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

176:   /* create the HYPRE vector for rhs and solution */
177:   PetscStackCallStandard(HYPRE_StructVectorCreate,(ex->hcomm,ex->hgrid,&ex->hb));
178:   PetscStackCallStandard(HYPRE_StructVectorCreate,(ex->hcomm,ex->hgrid,&ex->hx));
179:   PetscStackCallStandard(HYPRE_StructVectorInitialize,(ex->hb));
180:   PetscStackCallStandard(HYPRE_StructVectorInitialize,(ex->hx));
181:   PetscStackCallStandard(HYPRE_StructVectorAssemble,(ex->hb));
182:   PetscStackCallStandard(HYPRE_StructVectorAssemble,(ex->hx));

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

193:   /* set the global and local sizes of the matrix */
194:   DMDAGetCorners(da,0,0,0,&nx,&ny,&nz);
195:   MatSetSizes(mat,dof*nx*ny*nz,dof*nx*ny*nz,PETSC_DECIDE,PETSC_DECIDE);
196:   PetscLayoutSetBlockSize(mat->rmap,1);
197:   PetscLayoutSetBlockSize(mat->cmap,1);
198:   PetscLayoutSetUp(mat->rmap);
199:   PetscLayoutSetUp(mat->cmap);

201:   if (dim == 3) {
202:     mat->ops->setvalueslocal = MatSetValuesLocal_HYPREStruct_3d;
203:     mat->ops->zerorowslocal  = MatZeroRowsLocal_HYPREStruct_3d;
204:     mat->ops->zeroentries    = MatZeroEntries_HYPREStruct_3d;

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

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

222: PetscErrorCode MatMult_HYPREStruct(Mat A,Vec x,Vec y)
223: {
224:   PetscErrorCode  ierr;
225:   PetscScalar     *xx,*yy;
226:   PetscInt        ilower[3],iupper[3];
227:   Mat_HYPREStruct *mx = (Mat_HYPREStruct*)(A->data);

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

232:   iupper[0] += ilower[0] - 1;
233:   iupper[1] += ilower[1] - 1;
234:   iupper[2] += ilower[2] - 1;

236:   /* copy x values over to hypre */
237:   PetscStackCallStandard(HYPRE_StructVectorSetConstantValues,(mx->hb,0.0));
238:   VecGetArray(x,&xx);
239:   PetscStackCallStandard(HYPRE_StructVectorSetBoxValues,(mx->hb,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,xx));
240:   VecRestoreArray(x,&xx);
241:   PetscStackCallStandard(HYPRE_StructVectorAssemble,(mx->hb));
242:   PetscStackCallStandard(HYPRE_StructMatrixMatvec,(1.0,mx->hmat,mx->hb,0.0,mx->hx));

244:   /* copy solution values back to PETSc */
245:   VecGetArray(y,&yy);
246:   PetscStackCallStandard(HYPRE_StructVectorGetBoxValues,(mx->hx,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,yy));
247:   VecRestoreArray(y,&yy);
248:   return(0);
249: }

253: PetscErrorCode MatAssemblyEnd_HYPREStruct(Mat mat,MatAssemblyType mode)
254: {
255:   Mat_HYPREStruct *ex = (Mat_HYPREStruct*) mat->data;
256:   PetscErrorCode  ierr;

259:   PetscStackCallStandard(HYPRE_StructMatrixAssemble,(ex->hmat));
260:   /* PetscStackCallStandard(HYPRE_StructMatrixPrint,("dummy",ex->hmat,0)); */
261:   return(0);
262: }

266: PetscErrorCode MatZeroEntries_HYPREStruct(Mat mat)
267: {
269:   /* before the DMDA is set to the matrix the zero doesn't need to do anything */
270:   return(0);
271: }


276: PetscErrorCode MatDestroy_HYPREStruct(Mat mat)
277: {
278:   Mat_HYPREStruct *ex = (Mat_HYPREStruct*) mat->data;
279:   PetscErrorCode  ierr;

282:   PetscStackCallStandard(HYPRE_StructMatrixDestroy,(ex->hmat));
283:   PetscStackCallStandard(HYPRE_StructVectorDestroy,(ex->hx));
284:   PetscStackCallStandard(HYPRE_StructVectorDestroy,(ex->hb));
285:   return(0);
286: }


291: PETSC_EXTERN PetscErrorCode MatCreate_HYPREStruct(Mat B)
292: {
293:   Mat_HYPREStruct *ex;
294:   PetscErrorCode  ierr;

297:   PetscNewLog(B,&ex);
298:   B->data      = (void*)ex;
299:   B->rmap->bs  = 1;
300:   B->assembled = PETSC_FALSE;

302:   B->insertmode = NOT_SET_VALUES;

304:   B->ops->assemblyend = MatAssemblyEnd_HYPREStruct;
305:   B->ops->mult        = MatMult_HYPREStruct;
306:   B->ops->zeroentries = MatZeroEntries_HYPREStruct;
307:   B->ops->destroy     = MatDestroy_HYPREStruct;

309:   ex->needsinitialization = PETSC_TRUE;

311:   MPI_Comm_dup(PetscObjectComm((PetscObject)B),&(ex->hcomm));
312:   PetscObjectComposeFunction((PetscObject)B,"MatSetupDM_C",MatSetupDM_HYPREStruct);
313:   PetscObjectChangeTypeName((PetscObject)B,MATHYPRESTRUCT);
314:   return(0);
315: }

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


322:    Level: intermediate

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

327:           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
328:           be defined by a DMDA.

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

332: M*/

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

343:   PetscInt part = 0;          /* Petsc sstruct interface only allows 1 part */
344:   PetscInt ordering;
345:   PetscInt grid_rank, to_grid_rank;
346:   PetscInt var_type, to_var_type;
347:   PetscInt to_var_entry = 0;

349:   PetscInt nvars= ex->nvars;
350:   PetscInt row,*entries;

353:   PetscMalloc1(7*nvars,&entries);
354:   ordering= ex->dofs_order;  /* ordering= 0   nodal ordering
355:                                           1   variable ordering */
356:   /* stencil entries are orderer 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] = ex->xs + (row % ex->nx);
390:       index[1] = ex->ys + ((row/ex->nx) % ex->ny);
391:       index[2] = ex->zs + (row/(ex->nxny));

393:       if (addv == ADD_VALUES) {
394:         PetscStackCallStandard(HYPRE_SStructMatrixAddToValues,(ex->ss_mat,part,(HYPRE_Int *)index,var_type,ncol,(HYPRE_Int *)entries,(PetscScalar*)values));
395:       } else {
396:         PetscStackCallStandard(HYPRE_SStructMatrixSetValues,(ex->ss_mat,part,(HYPRE_Int *)index,var_type,ncol,(HYPRE_Int *)entries,(PetscScalar*)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] = ex->xs + (row % ex->nx);
432:       index[1] = ex->ys + ((row/ex->nx) % ex->ny);
433:       index[2] = ex->zs + (row/(ex->nxny));

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

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

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

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

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

473:   for (i=0; i< nvars; i++) {
474:     PetscMemzero(values[i],nvars*7*sizeof(PetscScalar));
475:     *(values[i]+3) = d;
476:   }

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

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

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

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

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

521:   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);
522:   {
523:     PetscInt    i,*entries,iupper[3],ilower[3];
524:     PetscScalar *values;


527:     for (i= 0; i< 3; i++) {
528:       ilower[i]= ex->hbox.imin[i];
529:       iupper[i]= ex->hbox.imax[i];
530:     }

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

536:     for (i= 0; i< nvars; i++) {
537:       PetscStackCallStandard(HYPRE_SStructMatrixSetBoxValues,(ex->ss_mat,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,nvars*7,(HYPRE_Int *)entries,values));
538:     }
539:     PetscFree2(entries,values);
540:   }
541:   PetscStackCallStandard(HYPRE_SStructMatrixAssemble,(ex->ss_mat));
542:   return(0);
543: }


548: static PetscErrorCode  MatSetupDM_HYPRESStruct(Mat mat,DM da)
549: {
550:   PetscErrorCode         ierr;
551:   Mat_HYPRESStruct       *ex = (Mat_HYPRESStruct*) mat->data;
552:   PetscInt               dim,dof,sw[3],nx,ny,nz;
553:   PetscInt               ilower[3],iupper[3],ssize,i;
554:   DMBoundaryType         px,py,pz;
555:   DMDAStencilType        st;
556:   PetscInt               nparts= 1;  /* assuming only one part */
557:   PetscInt               part  = 0;
558:   ISLocalToGlobalMapping ltog;
560:   ex->da = da;
561:   PetscObjectReference((PetscObject)da);

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

576:   ex->dofs_order = 0;

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

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

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

587:   {
588:     HYPRE_SStructVariable *vartypes;
589:     PetscMalloc1(ex->nvars,&vartypes);
590:     for (i= 0; i< ex->nvars; i++) vartypes[i]= HYPRE_SSTRUCT_VARIABLE_CELL;
591:     PetscStackCallStandard(HYPRE_SStructGridSetVariables,(ex->ss_grid, part, ex->nvars,vartypes));
592:     PetscFree(vartypes);
593:   }
594:   PetscStackCallStandard(HYPRE_SStructGridAssemble,(ex->ss_grid));

596:   sw[1] = sw[0];
597:   sw[2] = sw[1];
598:   /* PetscStackCallStandard(HYPRE_SStructGridSetNumGhost,(ex->ss_grid,sw)); */

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

604:   if (dim == 1) {
605:     PetscInt offsets[3][1] = {{-1},{0},{1}};
606:     PetscInt j, cnt;

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

621:     ssize = 5*(ex->nvars);
622:     PetscStackCallStandard(HYPRE_SStructStencilCreate,(dim,ssize,&ex->ss_stencil));
623:     cnt= 0;
624:     for (i= 0; i< (ex->nvars); i++) {
625:       for (j= 0; j< 5; j++) {
626:         PetscStackCallStandard(HYPRE_SStructStencilSetEntry,(ex->ss_stencil, cnt, (HYPRE_Int *)offsets[j], i));
627:         cnt++;
628:       }
629:     }
630:   } else if (dim == 3) {
631:     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}};
632:     PetscInt j, cnt;

634:     ssize = 7*(ex->nvars);
635:     PetscStackCallStandard(HYPRE_SStructStencilCreate,(dim,ssize,&ex->ss_stencil));
636:     cnt= 0;
637:     for (i= 0; i< (ex->nvars); i++) {
638:       for (j= 0; j< 7; j++) {
639:         PetscStackCallStandard(HYPRE_SStructStencilSetEntry,(ex->ss_stencil, cnt, (HYPRE_Int *)offsets[j], i));
640:         cnt++;
641:       }
642:     }
643:   }

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

648:   /* set the stencil graph. Note that each variable has the same graph. This means that each
649:      variable couples to all the other variable and with the same stencil pattern. */
650:   for (i= 0; i< (ex->nvars); i++) {
651:     PetscStackCallStandard(HYPRE_SStructGraphSetStencil,(ex->ss_graph,part,i,ex->ss_stencil));
652:   }
653:   PetscStackCallStandard(HYPRE_SStructGraphAssemble,(ex->ss_graph));

655:   /* create the HYPRE sstruct vectors for rhs and solution */
656:   PetscStackCallStandard(HYPRE_SStructVectorCreate,(ex->hcomm,ex->ss_grid,&ex->ss_b));
657:   PetscStackCallStandard(HYPRE_SStructVectorCreate,(ex->hcomm,ex->ss_grid,&ex->ss_x));
658:   PetscStackCallStandard(HYPRE_SStructVectorInitialize,(ex->ss_b));
659:   PetscStackCallStandard(HYPRE_SStructVectorInitialize,(ex->ss_x));
660:   PetscStackCallStandard(HYPRE_SStructVectorAssemble,(ex->ss_b));
661:   PetscStackCallStandard(HYPRE_SStructVectorAssemble,(ex->ss_x));

663:   /* create the hypre matrix object and set its information */
664:   PetscStackCallStandard(HYPRE_SStructMatrixCreate,(ex->hcomm,ex->ss_graph,&ex->ss_mat));
665:   PetscStackCallStandard(HYPRE_SStructGridDestroy,(ex->ss_grid));
666:   PetscStackCallStandard(HYPRE_SStructStencilDestroy,(ex->ss_stencil));
667:   if (ex->needsinitialization) {
668:     PetscStackCallStandard(HYPRE_SStructMatrixInitialize,(ex->ss_mat));
669:     ex->needsinitialization = PETSC_FALSE;
670:   }

672:   /* set the global and local sizes of the matrix */
673:   DMDAGetCorners(da,0,0,0,&nx,&ny,&nz);
674:   MatSetSizes(mat,dof*nx*ny*nz,dof*nx*ny*nz,PETSC_DECIDE,PETSC_DECIDE);
675:   PetscLayoutSetBlockSize(mat->rmap,1);
676:   PetscLayoutSetBlockSize(mat->cmap,1);
677:   PetscLayoutSetUp(mat->rmap);
678:   PetscLayoutSetUp(mat->cmap);

680:   if (dim == 3) {
681:     mat->ops->setvalueslocal = MatSetValuesLocal_HYPRESStruct_3d;
682:     mat->ops->zerorowslocal  = MatZeroRowsLocal_HYPRESStruct_3d;
683:     mat->ops->zeroentries    = MatZeroEntries_HYPRESStruct_3d;

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

688:   /* get values that will be used repeatedly in MatSetValuesLocal() and MatZeroRowsLocal() repeatedly */
689:   MatGetOwnershipRange(mat,&ex->rstart,NULL);
690:   DMGetLocalToGlobalMapping(ex->da,&ltog);
691:   ISLocalToGlobalMappingGetIndices(ltog, (const PetscInt **) &ex->gindices);
692:   DMDAGetGhostCorners(ex->da,0,0,0,&ex->gnx,&ex->gnxgny,&ex->gnxgnygnz);

694:   ex->gnxgny    *= ex->gnx;
695:   ex->gnxgnygnz *= ex->gnxgny;

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

699:   ex->nxny   = ex->nx*ex->ny;
700:   ex->nxnynz = ex->nz*ex->nxny;
701:   return(0);
702: }

706: PetscErrorCode MatMult_HYPRESStruct(Mat A,Vec x,Vec y)
707: {
708:   PetscErrorCode   ierr;
709:   PetscScalar      *xx,*yy;
710:   PetscInt         ilower[3],iupper[3];
711:   Mat_HYPRESStruct *mx     = (Mat_HYPRESStruct*)(A->data);
712:   PetscInt         ordering= mx->dofs_order;
713:   PetscInt         nvars   = mx->nvars;
714:   PetscInt         part    = 0;
715:   PetscInt         size;
716:   PetscInt         i;

719:   DMDAGetCorners(mx->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);
720:   iupper[0] += ilower[0] - 1;
721:   iupper[1] += ilower[1] - 1;
722:   iupper[2] += ilower[2] - 1;

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

727:   /* copy x values over to hypre for variable ordering */
728:   if (ordering) {
729:     PetscStackCallStandard(HYPRE_SStructVectorSetConstantValues,(mx->ss_b,0.0));
730:     VecGetArray(x,&xx);
731:     for (i= 0; i< nvars; i++) {
732:       PetscStackCallStandard(HYPRE_SStructVectorSetBoxValues,(mx->ss_b,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,xx+(size*i)));
733:     }
734:     VecRestoreArray(x,&xx);
735:     PetscStackCallStandard(HYPRE_SStructVectorAssemble,(mx->ss_b));
736:     PetscStackCallStandard(HYPRE_SStructMatrixMatvec,(1.0,mx->ss_mat,mx->ss_b,0.0,mx->ss_x));

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

748:     PetscMalloc1(nvars*size,&z);
749:     PetscStackCallStandard(HYPRE_SStructVectorSetConstantValues,(mx->ss_b,0.0));
750:     VecGetArray(x,&xx);

752:     /* transform nodal to hypre's variable ordering for sys_pfmg */
753:     for (i= 0; i< size; i++) {
754:       k= i*nvars;
755:       for (j= 0; j< nvars; j++) z[j*size+i]= xx[k+j];
756:     }
757:     for (i= 0; i< nvars; i++) {
758:       PetscStackCallStandard(HYPRE_SStructVectorSetBoxValues,(mx->ss_b,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,z+(size*i)));
759:     }
760:     VecRestoreArray(x,&xx);
761:     PetscStackCallStandard(HYPRE_SStructVectorAssemble,(mx->ss_b));
762:     PetscStackCallStandard(HYPRE_SStructMatrixMatvec,(1.0,mx->ss_mat,mx->ss_b,0.0,mx->ss_x));

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

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

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

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


804: PetscErrorCode MatDestroy_HYPRESStruct(Mat mat)
805: {
806:   Mat_HYPRESStruct *ex = (Mat_HYPRESStruct*) mat->data;
807:   PetscErrorCode   ierr;

810:   PetscStackCallStandard(HYPRE_SStructGraphDestroy,(ex->ss_graph));
811:   PetscStackCallStandard(HYPRE_SStructMatrixDestroy,(ex->ss_mat));
812:   PetscStackCallStandard(HYPRE_SStructVectorDestroy,(ex->ss_x));
813:   PetscStackCallStandard(HYPRE_SStructVectorDestroy,(ex->ss_b));
814:   return(0);
815: }

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

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

830:   B->insertmode = NOT_SET_VALUES;

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

837:   ex->needsinitialization = PETSC_TRUE;

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