Actual source code: mhyp.c

petsc-master 2018-02-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: 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 MatSetDM() or if the matrix is obtained from DMCreateMatrix()

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


 27: PetscErrorCode  MatSetValuesLocal_HYPREStruct_3d(Mat mat,PetscInt nrow,const PetscInt irow[],PetscInt ncol,const PetscInt icol[],const PetscScalar y[],InsertMode addv)
 28: {
 29:   PetscErrorCode    ierr;
 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:   PetscErrorCode  ierr;
 94:   PetscInt        indices[7] = {0,1,2,3,4,5,6};
 95:   Mat_HYPREStruct *ex        = (Mat_HYPREStruct*) mat->data;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

288:   B->insertmode = NOT_SET_VALUES;

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

296:   ex->needsinitialization = PETSC_TRUE;

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

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


308:    Level: intermediate

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

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

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

318: M*/

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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


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

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

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


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

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

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

555:   ex->dofs_order = 0;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

763:   PetscStackCallStandard(HYPRE_SStructMatrixAssemble,(ex->ss_mat));
764:   return(0);
765: }

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

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

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

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

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

804:   B->insertmode = NOT_SET_VALUES;

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

812:   ex->needsinitialization = PETSC_TRUE;

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