Actual source code: mhyp.c
petsc-3.13.6 2020-09-29
2: /*
3: Creates hypre ijmatrix from PETSc matrix
4: */
5: #include <petscsys.h>
6: #include <petsc/private/petschypre.h>
7: #include <petsc/private/matimpl.h>
8: #include <petscdmda.h>
9: #include <../src/dm/impls/da/hypre/mhyp.h>
11: /* -----------------------------------------------------------------------------------------------------------------*/
13: /*MC
14: MATHYPRESTRUCT - MATHYPRESTRUCT = "hyprestruct" - A matrix type to be used for parallel sparse matrices
15: based on the hypre HYPRE_StructMatrix.
17: Level: intermediate
19: Notes:
20: Unlike the more general support for blocks in hypre this allows only one block per process and requires the block
21: be defined by a DMDA.
23: The matrix needs a DMDA associated with it by either a call to MatSetDM() or if the matrix is obtained from DMCreateMatrix()
25: .seealso: MatCreate(), PCPFMG, MatSetDM(), DMCreateMatrix()
26: M*/
28: PetscErrorCode MatSetValuesLocal_HYPREStruct_3d(Mat mat,PetscInt nrow,const PetscInt irow[],PetscInt ncol,const PetscInt icol[],const PetscScalar y[],InsertMode addv)
29: {
30: HYPRE_Int index[3],entries[7];
31: PetscInt i,j,stencil,row;
32: HYPRE_Complex *values = (HYPRE_Complex*)y;
33: Mat_HYPREStruct *ex = (Mat_HYPREStruct*) mat->data;
36: if (PetscUnlikely(ncol > 7)) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"ncol %D > 7 too large",ncol);
37: for (i=0; i<nrow; i++) {
38: for (j=0; j<ncol; j++) {
39: stencil = icol[j] - irow[i];
40: if (!stencil) {
41: entries[j] = 3;
42: } else if (stencil == -1) {
43: entries[j] = 2;
44: } else if (stencil == 1) {
45: entries[j] = 4;
46: } else if (stencil == -ex->gnx) {
47: entries[j] = 1;
48: } else if (stencil == ex->gnx) {
49: entries[j] = 5;
50: } else if (stencil == -ex->gnxgny) {
51: entries[j] = 0;
52: } else if (stencil == ex->gnxgny) {
53: entries[j] = 6;
54: } else SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Local row %D local column %D have bad stencil %D",irow[i],icol[j],stencil);
55: }
56: row = ex->gindices[irow[i]] - ex->rstart;
57: index[0] = (HYPRE_Int)(ex->xs + (row % ex->nx));
58: index[1] = (HYPRE_Int)(ex->ys + ((row/ex->nx) % ex->ny));
59: index[2] = (HYPRE_Int)(ex->zs + (row/(ex->nxny)));
60: if (addv == ADD_VALUES) {
61: PetscStackCallStandard(HYPRE_StructMatrixAddToValues,(ex->hmat,index,ncol,entries,values));
62: } else {
63: PetscStackCallStandard(HYPRE_StructMatrixSetValues,(ex->hmat,index,ncol,entries,values));
64: }
65: values += ncol;
66: }
67: return(0);
68: }
70: PetscErrorCode MatZeroRowsLocal_HYPREStruct_3d(Mat mat,PetscInt nrow,const PetscInt irow[],PetscScalar d,Vec x,Vec b)
71: {
72: PetscErrorCode ierr;
73: HYPRE_Int index[3],entries[7] = {0,1,2,3,4,5,6};
74: PetscInt row,i;
75: HYPRE_Complex values[7];
76: Mat_HYPREStruct *ex = (Mat_HYPREStruct*) mat->data;
79: if (x && b) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"No support");
80: PetscArrayzero(values,7);
81: PetscHYPREScalarCast(d,&values[3]);
82: for (i=0; i<nrow; i++) {
83: row = ex->gindices[irow[i]] - ex->rstart;
84: index[0] = (HYPRE_Int)(ex->xs + (row % ex->nx));
85: index[1] = (HYPRE_Int)(ex->ys + ((row/ex->nx) % ex->ny));
86: index[2] = (HYPRE_Int)(ex->zs + (row/(ex->nxny)));
87: PetscStackCallStandard(HYPRE_StructMatrixSetValues,(ex->hmat,index,7,entries,values));
88: }
89: PetscStackCallStandard(HYPRE_StructMatrixAssemble,(ex->hmat));
90: return(0);
91: }
93: PetscErrorCode MatZeroEntries_HYPREStruct_3d(Mat mat)
94: {
95: HYPRE_Int indices[7] = {0,1,2,3,4,5,6};
96: Mat_HYPREStruct *ex = (Mat_HYPREStruct*) mat->data;
99: /* hypre has no public interface to do this */
100: PetscStackCallStandard(hypre_StructMatrixClearBoxValues,(ex->hmat,&ex->hbox,7,indices,0,1));
101: PetscStackCallStandard(HYPRE_StructMatrixAssemble,(ex->hmat));
102: return(0);
103: }
105: static PetscErrorCode MatSetUp_HYPREStruct(Mat mat)
106: {
107: PetscErrorCode ierr;
108: Mat_HYPREStruct *ex = (Mat_HYPREStruct*) mat->data;
109: HYPRE_Int sw[6];
110: HYPRE_Int hlower[3],hupper[3];
111: PetscInt dim,dof,psw,nx,ny,nz,ilower[3],iupper[3],ssize,i;
112: DMBoundaryType px,py,pz;
113: DMDAStencilType st;
114: ISLocalToGlobalMapping ltog;
115: DM da;
118: MatGetDM(mat,(DM*)&da);
119: ex->da = da;
120: PetscObjectReference((PetscObject)da);
122: DMDAGetInfo(ex->da,&dim,0,0,0,0,0,0,&dof,&psw,&px,&py,&pz,&st);
123: DMDAGetCorners(ex->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);
125: /* when HYPRE_MIXEDINT is defined, sizeof(HYPRE_Int) == 32 */
126: iupper[0] += ilower[0] - 1;
127: iupper[1] += ilower[1] - 1;
128: iupper[2] += ilower[2] - 1;
129: hlower[0] = (HYPRE_Int)ilower[0];
130: hlower[1] = (HYPRE_Int)ilower[1];
131: hlower[2] = (HYPRE_Int)ilower[2];
132: hupper[0] = (HYPRE_Int)iupper[0];
133: hupper[1] = (HYPRE_Int)iupper[1];
134: hupper[2] = (HYPRE_Int)iupper[2];
136: /* the hypre_Box is used to zero out the matrix entries in MatZeroValues() */
137: ex->hbox.imin[0] = hlower[0];
138: ex->hbox.imin[1] = hlower[1];
139: ex->hbox.imin[2] = hlower[2];
140: ex->hbox.imax[0] = hupper[0];
141: ex->hbox.imax[1] = hupper[1];
142: ex->hbox.imax[2] = hupper[2];
144: /* create the hypre grid object and set its information */
145: if (dof > 1) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Currently only support for scalar problems");
146: if (px || py || pz) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Ask us to add periodic support by calling HYPRE_StructGridSetPeriodic()");
147: PetscStackCallStandard(HYPRE_StructGridCreate,(ex->hcomm,dim,&ex->hgrid));
148: PetscStackCallStandard(HYPRE_StructGridSetExtents,(ex->hgrid,hlower,hupper));
149: PetscStackCallStandard(HYPRE_StructGridAssemble,(ex->hgrid));
151: sw[5] = sw[4] = sw[3] = sw[2] = sw[1] = sw[0] = psw;
152: PetscStackCallStandard(HYPRE_StructGridSetNumGhost,(ex->hgrid,sw));
154: /* create the hypre stencil object and set its information */
155: if (sw[0] > 1) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Ask us to add support for wider stencils");
156: if (st == DMDA_STENCIL_BOX) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Ask us to add support for box stencils");
157: if (dim == 1) {
158: HYPRE_Int offsets[3][1] = {{-1},{0},{1}};
159: ssize = 3;
160: PetscStackCallStandard(HYPRE_StructStencilCreate,(dim,ssize,&ex->hstencil));
161: for (i=0; i<ssize; i++) {
162: PetscStackCallStandard(HYPRE_StructStencilSetElement,(ex->hstencil,i,offsets[i]));
163: }
164: } else if (dim == 2) {
165: HYPRE_Int offsets[5][2] = {{0,-1},{-1,0},{0,0},{1,0},{0,1}};
166: ssize = 5;
167: PetscStackCallStandard(HYPRE_StructStencilCreate,(dim,ssize,&ex->hstencil));
168: for (i=0; i<ssize; i++) {
169: PetscStackCallStandard(HYPRE_StructStencilSetElement,(ex->hstencil,i,offsets[i]));
170: }
171: } else if (dim == 3) {
172: HYPRE_Int offsets[7][3] = {{0,0,-1},{0,-1,0},{-1,0,0},{0,0,0},{1,0,0},{0,1,0},{0,0,1}};
173: ssize = 7;
174: PetscStackCallStandard(HYPRE_StructStencilCreate,(dim,ssize,&ex->hstencil));
175: for (i=0; i<ssize; i++) {
176: PetscStackCallStandard(HYPRE_StructStencilSetElement,(ex->hstencil,i,offsets[i]));
177: }
178: }
180: /* create the HYPRE vector for rhs and solution */
181: PetscStackCallStandard(HYPRE_StructVectorCreate,(ex->hcomm,ex->hgrid,&ex->hb));
182: PetscStackCallStandard(HYPRE_StructVectorCreate,(ex->hcomm,ex->hgrid,&ex->hx));
183: PetscStackCallStandard(HYPRE_StructVectorInitialize,(ex->hb));
184: PetscStackCallStandard(HYPRE_StructVectorInitialize,(ex->hx));
185: PetscStackCallStandard(HYPRE_StructVectorAssemble,(ex->hb));
186: PetscStackCallStandard(HYPRE_StructVectorAssemble,(ex->hx));
188: /* create the hypre matrix object and set its information */
189: PetscStackCallStandard(HYPRE_StructMatrixCreate,(ex->hcomm,ex->hgrid,ex->hstencil,&ex->hmat));
190: PetscStackCallStandard(HYPRE_StructGridDestroy,(ex->hgrid));
191: PetscStackCallStandard(HYPRE_StructStencilDestroy,(ex->hstencil));
192: if (ex->needsinitialization) {
193: PetscStackCallStandard(HYPRE_StructMatrixInitialize,(ex->hmat));
194: ex->needsinitialization = PETSC_FALSE;
195: }
197: /* set the global and local sizes of the matrix */
198: DMDAGetCorners(da,0,0,0,&nx,&ny,&nz);
199: MatSetSizes(mat,dof*nx*ny*nz,dof*nx*ny*nz,PETSC_DECIDE,PETSC_DECIDE);
200: PetscLayoutSetBlockSize(mat->rmap,1);
201: PetscLayoutSetBlockSize(mat->cmap,1);
202: PetscLayoutSetUp(mat->rmap);
203: PetscLayoutSetUp(mat->cmap);
204: mat->preallocated = PETSC_TRUE;
206: if (dim == 3) {
207: mat->ops->setvalueslocal = MatSetValuesLocal_HYPREStruct_3d;
208: mat->ops->zerorowslocal = MatZeroRowsLocal_HYPREStruct_3d;
209: mat->ops->zeroentries = MatZeroEntries_HYPREStruct_3d;
211: /* MatZeroEntries_HYPREStruct_3d(mat); */
212: } else SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Only support for 3d DMDA currently");
214: /* get values that will be used repeatedly in MatSetValuesLocal() and MatZeroRowsLocal() repeatedly */
215: MatGetOwnershipRange(mat,&ex->rstart,NULL);
216: DMGetLocalToGlobalMapping(ex->da,<og);
217: ISLocalToGlobalMappingGetIndices(ltog, (const PetscInt **) &ex->gindices);
218: DMDAGetGhostCorners(ex->da,0,0,0,&ex->gnx,&ex->gnxgny,0);
219: ex->gnxgny *= ex->gnx;
220: DMDAGetCorners(ex->da,&ex->xs,&ex->ys,&ex->zs,&ex->nx,&ex->ny,0);
221: ex->nxny = ex->nx*ex->ny;
222: return(0);
223: }
225: PetscErrorCode MatMult_HYPREStruct(Mat A,Vec x,Vec y)
226: {
227: PetscErrorCode ierr;
228: const PetscScalar *xx;
229: PetscScalar *yy;
230: PetscInt ilower[3],iupper[3];
231: HYPRE_Int hlower[3],hupper[3];
232: Mat_HYPREStruct *mx = (Mat_HYPREStruct*)(A->data);
235: DMDAGetCorners(mx->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);
236: /* when HYPRE_MIXEDINT is defined, sizeof(HYPRE_Int) == 32 */
237: iupper[0] += ilower[0] - 1;
238: iupper[1] += ilower[1] - 1;
239: iupper[2] += ilower[2] - 1;
240: hlower[0] = (HYPRE_Int)ilower[0];
241: hlower[1] = (HYPRE_Int)ilower[1];
242: hlower[2] = (HYPRE_Int)ilower[2];
243: hupper[0] = (HYPRE_Int)iupper[0];
244: hupper[1] = (HYPRE_Int)iupper[1];
245: hupper[2] = (HYPRE_Int)iupper[2];
247: /* copy x values over to hypre */
248: PetscStackCallStandard(HYPRE_StructVectorSetConstantValues,(mx->hb,0.0));
249: VecGetArrayRead(x,&xx);
250: PetscStackCallStandard(HYPRE_StructVectorSetBoxValues,(mx->hb,hlower,hupper,(HYPRE_Complex*)xx));
251: VecRestoreArrayRead(x,&xx);
252: PetscStackCallStandard(HYPRE_StructVectorAssemble,(mx->hb));
253: PetscStackCallStandard(HYPRE_StructMatrixMatvec,(1.0,mx->hmat,mx->hb,0.0,mx->hx));
255: /* copy solution values back to PETSc */
256: VecGetArray(y,&yy);
257: PetscStackCallStandard(HYPRE_StructVectorGetBoxValues,(mx->hx,hlower,hupper,(HYPRE_Complex*)yy));
258: VecRestoreArray(y,&yy);
259: return(0);
260: }
262: PetscErrorCode MatAssemblyEnd_HYPREStruct(Mat mat,MatAssemblyType mode)
263: {
264: Mat_HYPREStruct *ex = (Mat_HYPREStruct*) mat->data;
267: PetscStackCallStandard(HYPRE_StructMatrixAssemble,(ex->hmat));
268: /* PetscStackCallStandard(HYPRE_StructMatrixPrint,("dummy",ex->hmat,0)); */
269: return(0);
270: }
272: PetscErrorCode MatZeroEntries_HYPREStruct(Mat mat)
273: {
275: /* before the DMDA is set to the matrix the zero doesn't need to do anything */
276: return(0);
277: }
279: PetscErrorCode MatDestroy_HYPREStruct(Mat mat)
280: {
281: Mat_HYPREStruct *ex = (Mat_HYPREStruct*) mat->data;
282: PetscErrorCode ierr;
285: PetscStackCallStandard(HYPRE_StructMatrixDestroy,(ex->hmat));
286: PetscStackCallStandard(HYPRE_StructVectorDestroy,(ex->hx));
287: PetscStackCallStandard(HYPRE_StructVectorDestroy,(ex->hb));
288: PetscObjectDereference((PetscObject)ex->da);
289: MPI_Comm_free(&(ex->hcomm));
290: PetscFree(ex);
291: return(0);
292: }
294: PETSC_EXTERN PetscErrorCode MatCreate_HYPREStruct(Mat B)
295: {
296: Mat_HYPREStruct *ex;
297: PetscErrorCode ierr;
300: PetscNewLog(B,&ex);
301: B->data = (void*)ex;
302: B->rmap->bs = 1;
303: B->assembled = PETSC_FALSE;
305: B->insertmode = NOT_SET_VALUES;
307: B->ops->assemblyend = MatAssemblyEnd_HYPREStruct;
308: B->ops->mult = MatMult_HYPREStruct;
309: B->ops->zeroentries = MatZeroEntries_HYPREStruct;
310: B->ops->destroy = MatDestroy_HYPREStruct;
311: B->ops->setup = MatSetUp_HYPREStruct;
313: ex->needsinitialization = PETSC_TRUE;
315: MPI_Comm_dup(PetscObjectComm((PetscObject)B),&(ex->hcomm));
316: PetscObjectChangeTypeName((PetscObject)B,MATHYPRESTRUCT);
317: return(0);
318: }
320: /*MC
321: MATHYPRESSTRUCT - MATHYPRESSTRUCT = "hypresstruct" - A matrix type to be used for parallel sparse matrices
322: based on the hypre HYPRE_SStructMatrix.
325: Level: intermediate
327: Notes:
328: Unlike hypre's general semi-struct object consisting of a collection of structured-grid objects and unstructured
329: grid objects, we restrict the semi-struct objects to consist of only structured-grid components.
331: 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
332: be defined by a DMDA.
334: The matrix needs a DMDA associated with it by either a call to MatSetDM() or if the matrix is obtained from DMCreateMatrix()
336: M*/
338: PetscErrorCode MatSetValuesLocal_HYPRESStruct_3d(Mat mat,PetscInt nrow,const PetscInt irow[],PetscInt ncol,const PetscInt icol[],const PetscScalar y[],InsertMode addv)
339: {
340: PetscErrorCode ierr;
341: HYPRE_Int index[3],*entries;
342: PetscInt i,j,stencil;
343: HYPRE_Complex *values = (HYPRE_Complex*)y;
344: Mat_HYPRESStruct *ex = (Mat_HYPRESStruct*) mat->data;
346: PetscInt part = 0; /* Petsc sstruct interface only allows 1 part */
347: PetscInt ordering;
348: PetscInt grid_rank, to_grid_rank;
349: PetscInt var_type, to_var_type;
350: PetscInt to_var_entry = 0;
351: PetscInt nvars= ex->nvars;
352: PetscInt row;
355: PetscMalloc1(7*nvars,&entries);
356: ordering = ex->dofs_order; /* ordering= 0 nodal ordering
357: 1 variable ordering */
358: /* stencil entries are ordered by variables: var0_stencil0, var0_stencil1, ..., var0_stencil6, var1_stencil0, var1_stencil1, ... */
359: if (!ordering) {
360: for (i=0; i<nrow; i++) {
361: grid_rank= irow[i]/nvars;
362: var_type = (irow[i] % nvars);
364: for (j=0; j<ncol; j++) {
365: to_grid_rank= icol[j]/nvars;
366: to_var_type = (icol[j] % nvars);
368: to_var_entry= to_var_entry*7;
369: entries[j] = to_var_entry;
371: stencil = to_grid_rank-grid_rank;
372: if (!stencil) {
373: entries[j] += 3;
374: } else if (stencil == -1) {
375: entries[j] += 2;
376: } else if (stencil == 1) {
377: entries[j] += 4;
378: } else if (stencil == -ex->gnx) {
379: entries[j] += 1;
380: } else if (stencil == ex->gnx) {
381: entries[j] += 5;
382: } else if (stencil == -ex->gnxgny) {
383: entries[j] += 0;
384: } else if (stencil == ex->gnxgny) {
385: entries[j] += 6;
386: } else SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Local row %D local column %D have bad stencil %D",irow[i],icol[j],stencil);
387: }
389: row = ex->gindices[grid_rank] - ex->rstart;
390: index[0] = (HYPRE_Int)(ex->xs + (row % ex->nx));
391: index[1] = (HYPRE_Int)(ex->ys + ((row/ex->nx) % ex->ny));
392: index[2] = (HYPRE_Int)(ex->zs + (row/(ex->nxny)));
394: if (addv == ADD_VALUES) {
395: PetscStackCallStandard(HYPRE_SStructMatrixAddToValues,(ex->ss_mat,part,index,var_type,ncol,entries,values));
396: } else {
397: PetscStackCallStandard(HYPRE_SStructMatrixSetValues,(ex->ss_mat,part,index,var_type,ncol,entries,values));
398: }
399: values += ncol;
400: }
401: } else {
402: for (i=0; i<nrow; i++) {
403: var_type = irow[i]/(ex->gnxgnygnz);
404: grid_rank= irow[i] - var_type*(ex->gnxgnygnz);
406: for (j=0; j<ncol; j++) {
407: to_var_type = icol[j]/(ex->gnxgnygnz);
408: to_grid_rank= icol[j] - to_var_type*(ex->gnxgnygnz);
410: to_var_entry= to_var_entry*7;
411: entries[j] = to_var_entry;
413: stencil = to_grid_rank-grid_rank;
414: if (!stencil) {
415: entries[j] += 3;
416: } else if (stencil == -1) {
417: entries[j] += 2;
418: } else if (stencil == 1) {
419: entries[j] += 4;
420: } else if (stencil == -ex->gnx) {
421: entries[j] += 1;
422: } else if (stencil == ex->gnx) {
423: entries[j] += 5;
424: } else if (stencil == -ex->gnxgny) {
425: entries[j] += 0;
426: } else if (stencil == ex->gnxgny) {
427: entries[j] += 6;
428: } else SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Local row %D local column %D have bad stencil %D",irow[i],icol[j],stencil);
429: }
431: row = ex->gindices[grid_rank] - ex->rstart;
432: index[0] = (HYPRE_Int)(ex->xs + (row % ex->nx));
433: index[1] = (HYPRE_Int)(ex->ys + ((row/ex->nx) % ex->ny));
434: index[2] = (HYPRE_Int)(ex->zs + (row/(ex->nxny)));
436: if (addv == ADD_VALUES) {
437: PetscStackCallStandard(HYPRE_SStructMatrixAddToValues,(ex->ss_mat,part,index,var_type,ncol,entries,values));
438: } else {
439: PetscStackCallStandard(HYPRE_SStructMatrixSetValues,(ex->ss_mat,part,index,var_type,ncol,entries,values));
440: }
441: values += ncol;
442: }
443: }
444: PetscFree(entries);
445: return(0);
446: }
448: PetscErrorCode MatZeroRowsLocal_HYPRESStruct_3d(Mat mat,PetscInt nrow,const PetscInt irow[],PetscScalar d,Vec x,Vec b)
449: {
450: PetscErrorCode ierr;
451: HYPRE_Int index[3],*entries;
452: PetscInt i;
453: HYPRE_Complex **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;
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: PetscArrayzero(values[i],nvars*7*sizeof(HYPRE_Complex));
475: PetscHYPREScalarCast(d,values[i]+3);
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] = (HYPRE_Int)(ex->xs + (row % ex->nx));
487: index[1] = (HYPRE_Int)(ex->ys + ((row/ex->nx) % ex->ny));
488: index[2] = (HYPRE_Int)(ex->zs + (row/(ex->nxny)));
489: PetscStackCallStandard(HYPRE_SStructMatrixSetValues,(ex->ss_mat,part,index,var_type,7*nvars,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] = (HYPRE_Int)(ex->xs + (row % ex->nx));
498: index[1] = (HYPRE_Int)(ex->ys + ((row/ex->nx) % ex->ny));
499: index[2] = (HYPRE_Int)(ex->zs + (row/(ex->nxny)));
500: PetscStackCallStandard(HYPRE_SStructMatrixSetValues,(ex->ss_mat,part,index,var_type,7*nvars,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: }
510: PetscErrorCode MatZeroEntries_HYPRESStruct_3d(Mat mat)
511: {
512: PetscErrorCode ierr;
513: Mat_HYPRESStruct *ex = (Mat_HYPRESStruct*) mat->data;
514: PetscInt nvars= ex->nvars;
515: PetscInt size;
516: PetscInt part= 0; /* only one part */
519: 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);
520: {
521: HYPRE_Int i,*entries,iupper[3],ilower[3];
522: HYPRE_Complex *values;
525: for (i= 0; i< 3; i++) {
526: ilower[i] = ex->hbox.imin[i];
527: iupper[i] = ex->hbox.imax[i];
528: }
530: PetscMalloc2(nvars*7,&entries,nvars*7*size,&values);
531: for (i= 0; i< nvars*7; i++) entries[i] = i;
532: PetscArrayzero(values,nvars*7*size);
534: for (i= 0; i< nvars; i++) {
535: PetscStackCallStandard(HYPRE_SStructMatrixSetBoxValues,(ex->ss_mat,part,ilower,iupper,i,nvars*7,entries,values));
536: }
537: PetscFree2(entries,values);
538: }
539: PetscStackCallStandard(HYPRE_SStructMatrixAssemble,(ex->ss_mat));
540: return(0);
541: }
544: static PetscErrorCode MatSetUp_HYPRESStruct(Mat mat)
545: {
546: PetscErrorCode ierr;
547: Mat_HYPRESStruct *ex = (Mat_HYPRESStruct*) mat->data;
548: PetscInt dim,dof,sw[3],nx,ny,nz;
549: PetscInt ilower[3],iupper[3],ssize,i;
550: DMBoundaryType px,py,pz;
551: DMDAStencilType st;
552: PetscInt nparts= 1; /* assuming only one part */
553: PetscInt part = 0;
554: ISLocalToGlobalMapping ltog;
555: DM da;
558: MatGetDM(mat,(DM*)&da);
559: ex->da = da;
560: PetscObjectReference((PetscObject)da);
562: DMDAGetInfo(ex->da,&dim,0,0,0,0,0,0,&dof,&sw[0],&px,&py,&pz,&st);
563: DMDAGetCorners(ex->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);
564: iupper[0] += ilower[0] - 1;
565: iupper[1] += ilower[1] - 1;
566: iupper[2] += ilower[2] - 1;
567: /* the hypre_Box is used to zero out the matrix entries in MatZeroValues() */
568: ex->hbox.imin[0] = ilower[0];
569: ex->hbox.imin[1] = ilower[1];
570: ex->hbox.imin[2] = ilower[2];
571: ex->hbox.imax[0] = iupper[0];
572: ex->hbox.imax[1] = iupper[1];
573: ex->hbox.imax[2] = iupper[2];
575: ex->dofs_order = 0;
577: /* assuming that the same number of dofs on each gridpoint. Also assume all cell-centred based */
578: ex->nvars= dof;
580: /* create the hypre grid object and set its information */
581: if (px || py || pz) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Ask us to add periodic support by calling HYPRE_SStructGridSetPeriodic()");
582: PetscStackCallStandard(HYPRE_SStructGridCreate,(ex->hcomm,dim,nparts,&ex->ss_grid));
583: PetscStackCallStandard(HYPRE_SStructGridSetExtents,(ex->ss_grid,part,ex->hbox.imin,ex->hbox.imax));
584: {
585: HYPRE_SStructVariable *vartypes;
586: PetscMalloc1(ex->nvars,&vartypes);
587: for (i= 0; i< ex->nvars; i++) vartypes[i]= HYPRE_SSTRUCT_VARIABLE_CELL;
588: PetscStackCallStandard(HYPRE_SStructGridSetVariables,(ex->ss_grid, part, ex->nvars,vartypes));
589: PetscFree(vartypes);
590: }
591: PetscStackCallStandard(HYPRE_SStructGridAssemble,(ex->ss_grid));
593: sw[1] = sw[0];
594: sw[2] = sw[1];
595: /* PetscStackCallStandard(HYPRE_SStructGridSetNumGhost,(ex->ss_grid,sw)); */
597: /* create the hypre stencil object and set its information */
598: if (sw[0] > 1) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Ask us to add support for wider stencils");
599: if (st == DMDA_STENCIL_BOX) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Ask us to add support for box stencils");
601: if (dim == 1) {
602: HYPRE_Int offsets[3][1] = {{-1},{0},{1}};
603: PetscInt j, cnt;
605: ssize = 3*(ex->nvars);
606: PetscStackCallStandard(HYPRE_SStructStencilCreate,(dim,ssize,&ex->ss_stencil));
607: cnt= 0;
608: for (i = 0; i < (ex->nvars); i++) {
609: for (j = 0; j < 3; j++) {
610: PetscStackCallStandard(HYPRE_SStructStencilSetEntry,(ex->ss_stencil, cnt, offsets[j], i));
611: cnt++;
612: }
613: }
614: } else if (dim == 2) {
615: HYPRE_Int offsets[5][2] = {{0,-1},{-1,0},{0,0},{1,0},{0,1}};
616: PetscInt j, cnt;
618: ssize = 5*(ex->nvars);
619: PetscStackCallStandard(HYPRE_SStructStencilCreate,(dim,ssize,&ex->ss_stencil));
620: cnt= 0;
621: for (i= 0; i< (ex->nvars); i++) {
622: for (j= 0; j< 5; j++) {
623: PetscStackCallStandard(HYPRE_SStructStencilSetEntry,(ex->ss_stencil, cnt, offsets[j], i));
624: cnt++;
625: }
626: }
627: } else if (dim == 3) {
628: HYPRE_Int offsets[7][3] = {{0,0,-1},{0,-1,0},{-1,0,0},{0,0,0},{1,0,0},{0,1,0},{0,0,1}};
629: PetscInt j, cnt;
631: ssize = 7*(ex->nvars);
632: PetscStackCallStandard(HYPRE_SStructStencilCreate,(dim,ssize,&ex->ss_stencil));
633: cnt= 0;
634: for (i= 0; i< (ex->nvars); i++) {
635: for (j= 0; j< 7; j++) {
636: PetscStackCallStandard(HYPRE_SStructStencilSetEntry,(ex->ss_stencil, cnt, offsets[j], i));
637: cnt++;
638: }
639: }
640: }
642: /* create the HYPRE graph */
643: PetscStackCallStandard(HYPRE_SStructGraphCreate,(ex->hcomm, ex->ss_grid, &(ex->ss_graph)));
645: /* set the stencil graph. Note that each variable has the same graph. This means that each
646: variable couples to all the other variable and with the same stencil pattern. */
647: for (i= 0; i< (ex->nvars); i++) {
648: PetscStackCallStandard(HYPRE_SStructGraphSetStencil,(ex->ss_graph,part,i,ex->ss_stencil));
649: }
650: PetscStackCallStandard(HYPRE_SStructGraphAssemble,(ex->ss_graph));
652: /* create the HYPRE sstruct vectors for rhs and solution */
653: PetscStackCallStandard(HYPRE_SStructVectorCreate,(ex->hcomm,ex->ss_grid,&ex->ss_b));
654: PetscStackCallStandard(HYPRE_SStructVectorCreate,(ex->hcomm,ex->ss_grid,&ex->ss_x));
655: PetscStackCallStandard(HYPRE_SStructVectorInitialize,(ex->ss_b));
656: PetscStackCallStandard(HYPRE_SStructVectorInitialize,(ex->ss_x));
657: PetscStackCallStandard(HYPRE_SStructVectorAssemble,(ex->ss_b));
658: PetscStackCallStandard(HYPRE_SStructVectorAssemble,(ex->ss_x));
660: /* create the hypre matrix object and set its information */
661: PetscStackCallStandard(HYPRE_SStructMatrixCreate,(ex->hcomm,ex->ss_graph,&ex->ss_mat));
662: PetscStackCallStandard(HYPRE_SStructGridDestroy,(ex->ss_grid));
663: PetscStackCallStandard(HYPRE_SStructStencilDestroy,(ex->ss_stencil));
664: if (ex->needsinitialization) {
665: PetscStackCallStandard(HYPRE_SStructMatrixInitialize,(ex->ss_mat));
666: ex->needsinitialization = PETSC_FALSE;
667: }
669: /* set the global and local sizes of the matrix */
670: DMDAGetCorners(da,0,0,0,&nx,&ny,&nz);
671: MatSetSizes(mat,dof*nx*ny*nz,dof*nx*ny*nz,PETSC_DECIDE,PETSC_DECIDE);
672: PetscLayoutSetBlockSize(mat->rmap,dof);
673: PetscLayoutSetBlockSize(mat->cmap,dof);
674: PetscLayoutSetUp(mat->rmap);
675: PetscLayoutSetUp(mat->cmap);
676: mat->preallocated = PETSC_TRUE;
678: if (dim == 3) {
679: mat->ops->setvalueslocal = MatSetValuesLocal_HYPRESStruct_3d;
680: mat->ops->zerorowslocal = MatZeroRowsLocal_HYPRESStruct_3d;
681: mat->ops->zeroentries = MatZeroEntries_HYPRESStruct_3d;
683: /* MatZeroEntries_HYPRESStruct_3d(mat); */
684: } else SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Only support for 3d DMDA currently");
686: /* get values that will be used repeatedly in MatSetValuesLocal() and MatZeroRowsLocal() repeatedly */
687: MatGetOwnershipRange(mat,&ex->rstart,NULL);
688: DMGetLocalToGlobalMapping(ex->da,<og);
689: ISLocalToGlobalMappingGetIndices(ltog, (const PetscInt **) &ex->gindices);
690: DMDAGetGhostCorners(ex->da,0,0,0,&ex->gnx,&ex->gnxgny,&ex->gnxgnygnz);
692: ex->gnxgny *= ex->gnx;
693: ex->gnxgnygnz *= ex->gnxgny;
695: DMDAGetCorners(ex->da,&ex->xs,&ex->ys,&ex->zs,&ex->nx,&ex->ny,&ex->nz);
697: ex->nxny = ex->nx*ex->ny;
698: ex->nxnynz = ex->nz*ex->nxny;
699: return(0);
700: }
702: PetscErrorCode MatMult_HYPRESStruct(Mat A,Vec x,Vec y)
703: {
704: PetscErrorCode ierr;
705: const PetscScalar *xx;
706: PetscScalar *yy;
707: HYPRE_Int hlower[3],hupper[3];
708: PetscInt ilower[3],iupper[3];
709: Mat_HYPRESStruct *mx = (Mat_HYPRESStruct*)(A->data);
710: PetscInt ordering= mx->dofs_order;
711: PetscInt nvars = mx->nvars;
712: PetscInt part = 0;
713: PetscInt size;
714: PetscInt i;
717: DMDAGetCorners(mx->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);
719: /* when HYPRE_MIXEDINT is defined, sizeof(HYPRE_Int) == 32 */
720: iupper[0] += ilower[0] - 1;
721: iupper[1] += ilower[1] - 1;
722: iupper[2] += ilower[2] - 1;
723: hlower[0] = (HYPRE_Int)ilower[0];
724: hlower[1] = (HYPRE_Int)ilower[1];
725: hlower[2] = (HYPRE_Int)ilower[2];
726: hupper[0] = (HYPRE_Int)iupper[0];
727: hupper[1] = (HYPRE_Int)iupper[1];
728: hupper[2] = (HYPRE_Int)iupper[2];
730: size= 1;
731: for (i = 0; i < 3; i++) size *= (iupper[i]-ilower[i]+1);
733: /* copy x values over to hypre for variable ordering */
734: if (ordering) {
735: PetscStackCallStandard(HYPRE_SStructVectorSetConstantValues,(mx->ss_b,0.0));
736: VecGetArrayRead(x,&xx);
737: for (i= 0; i< nvars; i++) {
738: PetscStackCallStandard(HYPRE_SStructVectorSetBoxValues,(mx->ss_b,part,hlower,hupper,i,(HYPRE_Complex*)(xx+(size*i))));
739: }
740: VecRestoreArrayRead(x,&xx);
741: PetscStackCallStandard(HYPRE_SStructVectorAssemble,(mx->ss_b));
742: PetscStackCallStandard(HYPRE_SStructMatrixMatvec,(1.0,mx->ss_mat,mx->ss_b,0.0,mx->ss_x));
744: /* copy solution values back to PETSc */
745: VecGetArray(y,&yy);
746: for (i= 0; i< nvars; i++) {
747: PetscStackCallStandard(HYPRE_SStructVectorGetBoxValues,(mx->ss_x,part,hlower,hupper,i,(HYPRE_Complex*)(yy+(size*i))));
748: }
749: VecRestoreArray(y,&yy);
750: } else { /* nodal ordering must be mapped to variable ordering for sys_pfmg */
751: PetscScalar *z;
752: PetscInt j, k;
754: PetscMalloc1(nvars*size,&z);
755: PetscStackCallStandard(HYPRE_SStructVectorSetConstantValues,(mx->ss_b,0.0));
756: VecGetArrayRead(x,&xx);
758: /* transform nodal to hypre's variable ordering for sys_pfmg */
759: for (i= 0; i< size; i++) {
760: k= i*nvars;
761: for (j= 0; j< nvars; j++) z[j*size+i]= xx[k+j];
762: }
763: for (i= 0; i< nvars; i++) {
764: PetscStackCallStandard(HYPRE_SStructVectorSetBoxValues,(mx->ss_b,part,hlower,hupper,i,(HYPRE_Complex*)(z+(size*i))));
765: }
766: VecRestoreArrayRead(x,&xx);
767: PetscStackCallStandard(HYPRE_SStructVectorAssemble,(mx->ss_b));
768: PetscStackCallStandard(HYPRE_SStructMatrixMatvec,(1.0,mx->ss_mat,mx->ss_b,0.0,mx->ss_x));
770: /* copy solution values back to PETSc */
771: VecGetArray(y,&yy);
772: for (i= 0; i< nvars; i++) {
773: PetscStackCallStandard(HYPRE_SStructVectorGetBoxValues,(mx->ss_x,part,hlower,hupper,i,(HYPRE_Complex*)(z+(size*i))));
774: }
775: /* transform hypre's variable ordering for sys_pfmg to nodal ordering */
776: for (i= 0; i< size; i++) {
777: k= i*nvars;
778: for (j= 0; j< nvars; j++) yy[k+j]= z[j*size+i];
779: }
780: VecRestoreArray(y,&yy);
781: PetscFree(z);
782: }
783: return(0);
784: }
786: PetscErrorCode MatAssemblyEnd_HYPRESStruct(Mat mat,MatAssemblyType mode)
787: {
788: Mat_HYPRESStruct *ex = (Mat_HYPRESStruct*) mat->data;
791: PetscStackCallStandard(HYPRE_SStructMatrixAssemble,(ex->ss_mat));
792: return(0);
793: }
795: PetscErrorCode MatZeroEntries_HYPRESStruct(Mat mat)
796: {
798: /* before the DMDA is set to the matrix the zero doesn't need to do anything */
799: return(0);
800: }
802: PetscErrorCode MatDestroy_HYPRESStruct(Mat mat)
803: {
804: Mat_HYPRESStruct *ex = (Mat_HYPRESStruct*) mat->data;
805: PetscErrorCode ierr;
806: ISLocalToGlobalMapping ltog;
809: DMGetLocalToGlobalMapping(ex->da,<og);
810: ISLocalToGlobalMappingRestoreIndices(ltog, (const PetscInt **) &ex->gindices);
811: PetscStackCallStandard(HYPRE_SStructGraphDestroy,(ex->ss_graph));
812: PetscStackCallStandard(HYPRE_SStructMatrixDestroy,(ex->ss_mat));
813: PetscStackCallStandard(HYPRE_SStructVectorDestroy,(ex->ss_x));
814: PetscStackCallStandard(HYPRE_SStructVectorDestroy,(ex->ss_b));
815: PetscObjectDereference((PetscObject)ex->da);
816: MPI_Comm_free(&(ex->hcomm));
817: PetscFree(ex);
818: return(0);
819: }
821: PETSC_EXTERN PetscErrorCode MatCreate_HYPRESStruct(Mat B)
822: {
823: Mat_HYPRESStruct *ex;
824: PetscErrorCode ierr;
827: PetscNewLog(B,&ex);
828: B->data = (void*)ex;
829: B->rmap->bs = 1;
830: B->assembled = PETSC_FALSE;
832: B->insertmode = NOT_SET_VALUES;
834: B->ops->assemblyend = MatAssemblyEnd_HYPRESStruct;
835: B->ops->mult = MatMult_HYPRESStruct;
836: B->ops->zeroentries = MatZeroEntries_HYPRESStruct;
837: B->ops->destroy = MatDestroy_HYPRESStruct;
838: B->ops->setup = MatSetUp_HYPRESStruct;
840: ex->needsinitialization = PETSC_TRUE;
842: MPI_Comm_dup(PetscObjectComm((PetscObject)B),&(ex->hcomm));
843: PetscObjectChangeTypeName((PetscObject)B,MATHYPRESSTRUCT);
844: return(0);
845: }