Actual source code: mhyp.c
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.
324: Level: intermediate
326: Notes:
327: Unlike hypre's general semi-struct object consisting of a collection of structured-grid objects and unstructured
328: grid objects, we restrict the semi-struct objects to consist of only structured-grid components.
330: 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
331: be defined by a DMDA.
333: The matrix needs a DMDA associated with it by either a call to MatSetDM() or if the matrix is obtained from DMCreateMatrix()
335: M*/
337: PetscErrorCode MatSetValuesLocal_HYPRESStruct_3d(Mat mat,PetscInt nrow,const PetscInt irow[],PetscInt ncol,const PetscInt icol[],const PetscScalar y[],InsertMode addv)
338: {
339: PetscErrorCode ierr;
340: HYPRE_Int index[3],*entries;
341: PetscInt i,j,stencil;
342: HYPRE_Complex *values = (HYPRE_Complex*)y;
343: Mat_HYPRESStruct *ex = (Mat_HYPRESStruct*) mat->data;
345: PetscInt part = 0; /* Petsc sstruct interface only allows 1 part */
346: PetscInt ordering;
347: PetscInt grid_rank, to_grid_rank;
348: PetscInt var_type, to_var_type;
349: PetscInt to_var_entry = 0;
350: PetscInt nvars= ex->nvars;
351: PetscInt row;
354: PetscMalloc1(7*nvars,&entries);
355: ordering = ex->dofs_order; /* ordering= 0 nodal ordering
356: 1 variable ordering */
357: /* stencil entries are ordered 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] = (HYPRE_Int)(ex->xs + (row % ex->nx));
390: index[1] = (HYPRE_Int)(ex->ys + ((row/ex->nx) % ex->ny));
391: index[2] = (HYPRE_Int)(ex->zs + (row/(ex->nxny)));
393: if (addv == ADD_VALUES) {
394: PetscStackCallStandard(HYPRE_SStructMatrixAddToValues,(ex->ss_mat,part,index,var_type,ncol,entries,values));
395: } else {
396: PetscStackCallStandard(HYPRE_SStructMatrixSetValues,(ex->ss_mat,part,index,var_type,ncol,entries,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] = (HYPRE_Int)(ex->xs + (row % ex->nx));
432: index[1] = (HYPRE_Int)(ex->ys + ((row/ex->nx) % ex->ny));
433: index[2] = (HYPRE_Int)(ex->zs + (row/(ex->nxny)));
435: if (addv == ADD_VALUES) {
436: PetscStackCallStandard(HYPRE_SStructMatrixAddToValues,(ex->ss_mat,part,index,var_type,ncol,entries,values));
437: } else {
438: PetscStackCallStandard(HYPRE_SStructMatrixSetValues,(ex->ss_mat,part,index,var_type,ncol,entries,values));
439: }
440: values += ncol;
441: }
442: }
443: PetscFree(entries);
444: return(0);
445: }
447: PetscErrorCode MatZeroRowsLocal_HYPRESStruct_3d(Mat mat,PetscInt nrow,const PetscInt irow[],PetscScalar d,Vec x,Vec b)
448: {
449: PetscErrorCode ierr;
450: HYPRE_Int index[3],*entries;
451: PetscInt i;
452: HYPRE_Complex **values;
453: Mat_HYPRESStruct *ex = (Mat_HYPRESStruct*) mat->data;
455: PetscInt part = 0; /* Petsc sstruct interface only allows 1 part */
456: PetscInt ordering = ex->dofs_order;
457: PetscInt grid_rank;
458: PetscInt var_type;
459: PetscInt nvars= ex->nvars;
460: PetscInt row;
463: if (x && b) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"No support");
464: PetscMalloc1(7*nvars,&entries);
466: PetscMalloc1(nvars,&values);
467: PetscMalloc1(7*nvars*nvars,&values[0]);
468: for (i=1; i<nvars; i++) {
469: values[i] = values[i-1] + nvars*7;
470: }
472: for (i=0; i< nvars; i++) {
473: PetscArrayzero(values[i],nvars*7*sizeof(HYPRE_Complex));
474: PetscHYPREScalarCast(d,values[i]+3);
475: }
477: for (i=0; i< nvars*7; i++) entries[i] = i;
479: if (!ordering) {
480: for (i=0; i<nrow; i++) {
481: grid_rank = irow[i] / nvars;
482: var_type = (irow[i] % nvars);
484: row = ex->gindices[grid_rank] - ex->rstart;
485: index[0] = (HYPRE_Int)(ex->xs + (row % ex->nx));
486: index[1] = (HYPRE_Int)(ex->ys + ((row/ex->nx) % ex->ny));
487: index[2] = (HYPRE_Int)(ex->zs + (row/(ex->nxny)));
488: PetscStackCallStandard(HYPRE_SStructMatrixSetValues,(ex->ss_mat,part,index,var_type,7*nvars,entries,values[var_type]));
489: }
490: } else {
491: for (i=0; i<nrow; i++) {
492: var_type = irow[i]/(ex->gnxgnygnz);
493: grid_rank= irow[i] - var_type*(ex->gnxgnygnz);
495: row = ex->gindices[grid_rank] - ex->rstart;
496: index[0] = (HYPRE_Int)(ex->xs + (row % ex->nx));
497: index[1] = (HYPRE_Int)(ex->ys + ((row/ex->nx) % ex->ny));
498: index[2] = (HYPRE_Int)(ex->zs + (row/(ex->nxny)));
499: PetscStackCallStandard(HYPRE_SStructMatrixSetValues,(ex->ss_mat,part,index,var_type,7*nvars,entries,values[var_type]));
500: }
501: }
502: PetscStackCallStandard(HYPRE_SStructMatrixAssemble,(ex->ss_mat));
503: PetscFree(values[0]);
504: PetscFree(values);
505: PetscFree(entries);
506: return(0);
507: }
509: PetscErrorCode MatZeroEntries_HYPRESStruct_3d(Mat mat)
510: {
511: PetscErrorCode ierr;
512: Mat_HYPRESStruct *ex = (Mat_HYPRESStruct*) mat->data;
513: PetscInt nvars= ex->nvars;
514: PetscInt size;
515: PetscInt part= 0; /* only one part */
518: 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);
519: {
520: HYPRE_Int i,*entries,iupper[3],ilower[3];
521: HYPRE_Complex *values;
523: for (i= 0; i< 3; i++) {
524: ilower[i] = ex->hbox.imin[i];
525: iupper[i] = ex->hbox.imax[i];
526: }
528: PetscMalloc2(nvars*7,&entries,nvars*7*size,&values);
529: for (i= 0; i< nvars*7; i++) entries[i] = i;
530: PetscArrayzero(values,nvars*7*size);
532: for (i= 0; i< nvars; i++) {
533: PetscStackCallStandard(HYPRE_SStructMatrixSetBoxValues,(ex->ss_mat,part,ilower,iupper,i,nvars*7,entries,values));
534: }
535: PetscFree2(entries,values);
536: }
537: PetscStackCallStandard(HYPRE_SStructMatrixAssemble,(ex->ss_mat));
538: return(0);
539: }
541: static PetscErrorCode MatSetUp_HYPRESStruct(Mat mat)
542: {
543: PetscErrorCode ierr;
544: Mat_HYPRESStruct *ex = (Mat_HYPRESStruct*) mat->data;
545: PetscInt dim,dof,sw[3],nx,ny,nz;
546: PetscInt ilower[3],iupper[3],ssize,i;
547: DMBoundaryType px,py,pz;
548: DMDAStencilType st;
549: PetscInt nparts= 1; /* assuming only one part */
550: PetscInt part = 0;
551: ISLocalToGlobalMapping ltog;
552: DM da;
555: MatGetDM(mat,(DM*)&da);
556: ex->da = da;
557: PetscObjectReference((PetscObject)da);
559: DMDAGetInfo(ex->da,&dim,0,0,0,0,0,0,&dof,&sw[0],&px,&py,&pz,&st);
560: DMDAGetCorners(ex->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);
561: iupper[0] += ilower[0] - 1;
562: iupper[1] += ilower[1] - 1;
563: iupper[2] += ilower[2] - 1;
564: /* the hypre_Box is used to zero out the matrix entries in MatZeroValues() */
565: ex->hbox.imin[0] = ilower[0];
566: ex->hbox.imin[1] = ilower[1];
567: ex->hbox.imin[2] = ilower[2];
568: ex->hbox.imax[0] = iupper[0];
569: ex->hbox.imax[1] = iupper[1];
570: ex->hbox.imax[2] = iupper[2];
572: ex->dofs_order = 0;
574: /* assuming that the same number of dofs on each gridpoint. Also assume all cell-centred based */
575: ex->nvars= dof;
577: /* create the hypre grid object and set its information */
578: if (px || py || pz) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Ask us to add periodic support by calling HYPRE_SStructGridSetPeriodic()");
579: PetscStackCallStandard(HYPRE_SStructGridCreate,(ex->hcomm,dim,nparts,&ex->ss_grid));
580: PetscStackCallStandard(HYPRE_SStructGridSetExtents,(ex->ss_grid,part,ex->hbox.imin,ex->hbox.imax));
581: {
582: HYPRE_SStructVariable *vartypes;
583: PetscMalloc1(ex->nvars,&vartypes);
584: for (i= 0; i< ex->nvars; i++) vartypes[i]= HYPRE_SSTRUCT_VARIABLE_CELL;
585: PetscStackCallStandard(HYPRE_SStructGridSetVariables,(ex->ss_grid, part, ex->nvars,vartypes));
586: PetscFree(vartypes);
587: }
588: PetscStackCallStandard(HYPRE_SStructGridAssemble,(ex->ss_grid));
590: sw[1] = sw[0];
591: sw[2] = sw[1];
592: /* PetscStackCallStandard(HYPRE_SStructGridSetNumGhost,(ex->ss_grid,sw)); */
594: /* create the hypre stencil object and set its information */
595: if (sw[0] > 1) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Ask us to add support for wider stencils");
596: if (st == DMDA_STENCIL_BOX) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Ask us to add support for box stencils");
598: if (dim == 1) {
599: HYPRE_Int offsets[3][1] = {{-1},{0},{1}};
600: PetscInt j, cnt;
602: ssize = 3*(ex->nvars);
603: PetscStackCallStandard(HYPRE_SStructStencilCreate,(dim,ssize,&ex->ss_stencil));
604: cnt= 0;
605: for (i = 0; i < (ex->nvars); i++) {
606: for (j = 0; j < 3; j++) {
607: PetscStackCallStandard(HYPRE_SStructStencilSetEntry,(ex->ss_stencil, cnt, offsets[j], i));
608: cnt++;
609: }
610: }
611: } else if (dim == 2) {
612: HYPRE_Int offsets[5][2] = {{0,-1},{-1,0},{0,0},{1,0},{0,1}};
613: PetscInt j, cnt;
615: ssize = 5*(ex->nvars);
616: PetscStackCallStandard(HYPRE_SStructStencilCreate,(dim,ssize,&ex->ss_stencil));
617: cnt= 0;
618: for (i= 0; i< (ex->nvars); i++) {
619: for (j= 0; j< 5; j++) {
620: PetscStackCallStandard(HYPRE_SStructStencilSetEntry,(ex->ss_stencil, cnt, offsets[j], i));
621: cnt++;
622: }
623: }
624: } else if (dim == 3) {
625: 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}};
626: PetscInt j, cnt;
628: ssize = 7*(ex->nvars);
629: PetscStackCallStandard(HYPRE_SStructStencilCreate,(dim,ssize,&ex->ss_stencil));
630: cnt= 0;
631: for (i= 0; i< (ex->nvars); i++) {
632: for (j= 0; j< 7; j++) {
633: PetscStackCallStandard(HYPRE_SStructStencilSetEntry,(ex->ss_stencil, cnt, offsets[j], i));
634: cnt++;
635: }
636: }
637: }
639: /* create the HYPRE graph */
640: PetscStackCallStandard(HYPRE_SStructGraphCreate,(ex->hcomm, ex->ss_grid, &(ex->ss_graph)));
642: /* set the stencil graph. Note that each variable has the same graph. This means that each
643: variable couples to all the other variable and with the same stencil pattern. */
644: for (i= 0; i< (ex->nvars); i++) {
645: PetscStackCallStandard(HYPRE_SStructGraphSetStencil,(ex->ss_graph,part,i,ex->ss_stencil));
646: }
647: PetscStackCallStandard(HYPRE_SStructGraphAssemble,(ex->ss_graph));
649: /* create the HYPRE sstruct vectors for rhs and solution */
650: PetscStackCallStandard(HYPRE_SStructVectorCreate,(ex->hcomm,ex->ss_grid,&ex->ss_b));
651: PetscStackCallStandard(HYPRE_SStructVectorCreate,(ex->hcomm,ex->ss_grid,&ex->ss_x));
652: PetscStackCallStandard(HYPRE_SStructVectorInitialize,(ex->ss_b));
653: PetscStackCallStandard(HYPRE_SStructVectorInitialize,(ex->ss_x));
654: PetscStackCallStandard(HYPRE_SStructVectorAssemble,(ex->ss_b));
655: PetscStackCallStandard(HYPRE_SStructVectorAssemble,(ex->ss_x));
657: /* create the hypre matrix object and set its information */
658: PetscStackCallStandard(HYPRE_SStructMatrixCreate,(ex->hcomm,ex->ss_graph,&ex->ss_mat));
659: PetscStackCallStandard(HYPRE_SStructGridDestroy,(ex->ss_grid));
660: PetscStackCallStandard(HYPRE_SStructStencilDestroy,(ex->ss_stencil));
661: if (ex->needsinitialization) {
662: PetscStackCallStandard(HYPRE_SStructMatrixInitialize,(ex->ss_mat));
663: ex->needsinitialization = PETSC_FALSE;
664: }
666: /* set the global and local sizes of the matrix */
667: DMDAGetCorners(da,0,0,0,&nx,&ny,&nz);
668: MatSetSizes(mat,dof*nx*ny*nz,dof*nx*ny*nz,PETSC_DECIDE,PETSC_DECIDE);
669: PetscLayoutSetBlockSize(mat->rmap,dof);
670: PetscLayoutSetBlockSize(mat->cmap,dof);
671: PetscLayoutSetUp(mat->rmap);
672: PetscLayoutSetUp(mat->cmap);
673: mat->preallocated = PETSC_TRUE;
675: if (dim == 3) {
676: mat->ops->setvalueslocal = MatSetValuesLocal_HYPRESStruct_3d;
677: mat->ops->zerorowslocal = MatZeroRowsLocal_HYPRESStruct_3d;
678: mat->ops->zeroentries = MatZeroEntries_HYPRESStruct_3d;
680: /* MatZeroEntries_HYPRESStruct_3d(mat); */
681: } else SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Only support for 3d DMDA currently");
683: /* get values that will be used repeatedly in MatSetValuesLocal() and MatZeroRowsLocal() repeatedly */
684: MatGetOwnershipRange(mat,&ex->rstart,NULL);
685: DMGetLocalToGlobalMapping(ex->da,<og);
686: ISLocalToGlobalMappingGetIndices(ltog, (const PetscInt **) &ex->gindices);
687: DMDAGetGhostCorners(ex->da,0,0,0,&ex->gnx,&ex->gnxgny,&ex->gnxgnygnz);
689: ex->gnxgny *= ex->gnx;
690: ex->gnxgnygnz *= ex->gnxgny;
692: DMDAGetCorners(ex->da,&ex->xs,&ex->ys,&ex->zs,&ex->nx,&ex->ny,&ex->nz);
694: ex->nxny = ex->nx*ex->ny;
695: ex->nxnynz = ex->nz*ex->nxny;
696: return(0);
697: }
699: PetscErrorCode MatMult_HYPRESStruct(Mat A,Vec x,Vec y)
700: {
701: PetscErrorCode ierr;
702: const PetscScalar *xx;
703: PetscScalar *yy;
704: HYPRE_Int hlower[3],hupper[3];
705: PetscInt ilower[3],iupper[3];
706: Mat_HYPRESStruct *mx = (Mat_HYPRESStruct*)(A->data);
707: PetscInt ordering= mx->dofs_order;
708: PetscInt nvars = mx->nvars;
709: PetscInt part = 0;
710: PetscInt size;
711: PetscInt i;
714: DMDAGetCorners(mx->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);
716: /* when HYPRE_MIXEDINT is defined, sizeof(HYPRE_Int) == 32 */
717: iupper[0] += ilower[0] - 1;
718: iupper[1] += ilower[1] - 1;
719: iupper[2] += ilower[2] - 1;
720: hlower[0] = (HYPRE_Int)ilower[0];
721: hlower[1] = (HYPRE_Int)ilower[1];
722: hlower[2] = (HYPRE_Int)ilower[2];
723: hupper[0] = (HYPRE_Int)iupper[0];
724: hupper[1] = (HYPRE_Int)iupper[1];
725: hupper[2] = (HYPRE_Int)iupper[2];
727: size= 1;
728: for (i = 0; i < 3; i++) size *= (iupper[i]-ilower[i]+1);
730: /* copy x values over to hypre for variable ordering */
731: if (ordering) {
732: PetscStackCallStandard(HYPRE_SStructVectorSetConstantValues,(mx->ss_b,0.0));
733: VecGetArrayRead(x,&xx);
734: for (i= 0; i< nvars; i++) {
735: PetscStackCallStandard(HYPRE_SStructVectorSetBoxValues,(mx->ss_b,part,hlower,hupper,i,(HYPRE_Complex*)(xx+(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,hlower,hupper,i,(HYPRE_Complex*)(yy+(size*i))));
745: }
746: VecRestoreArray(y,&yy);
747: } else { /* nodal ordering must be mapped to variable ordering for sys_pfmg */
748: PetscScalar *z;
749: PetscInt j, k;
751: PetscMalloc1(nvars*size,&z);
752: PetscStackCallStandard(HYPRE_SStructVectorSetConstantValues,(mx->ss_b,0.0));
753: VecGetArrayRead(x,&xx);
755: /* transform nodal to hypre's variable ordering for sys_pfmg */
756: for (i= 0; i< size; i++) {
757: k= i*nvars;
758: for (j= 0; j< nvars; j++) z[j*size+i]= xx[k+j];
759: }
760: for (i= 0; i< nvars; i++) {
761: PetscStackCallStandard(HYPRE_SStructVectorSetBoxValues,(mx->ss_b,part,hlower,hupper,i,(HYPRE_Complex*)(z+(size*i))));
762: }
763: VecRestoreArrayRead(x,&xx);
764: PetscStackCallStandard(HYPRE_SStructVectorAssemble,(mx->ss_b));
765: PetscStackCallStandard(HYPRE_SStructMatrixMatvec,(1.0,mx->ss_mat,mx->ss_b,0.0,mx->ss_x));
767: /* copy solution values back to PETSc */
768: VecGetArray(y,&yy);
769: for (i= 0; i< nvars; i++) {
770: PetscStackCallStandard(HYPRE_SStructVectorGetBoxValues,(mx->ss_x,part,hlower,hupper,i,(HYPRE_Complex*)(z+(size*i))));
771: }
772: /* transform hypre's variable ordering for sys_pfmg to nodal ordering */
773: for (i= 0; i< size; i++) {
774: k= i*nvars;
775: for (j= 0; j< nvars; j++) yy[k+j]= z[j*size+i];
776: }
777: VecRestoreArray(y,&yy);
778: PetscFree(z);
779: }
780: return(0);
781: }
783: PetscErrorCode MatAssemblyEnd_HYPRESStruct(Mat mat,MatAssemblyType mode)
784: {
785: Mat_HYPRESStruct *ex = (Mat_HYPRESStruct*) mat->data;
788: PetscStackCallStandard(HYPRE_SStructMatrixAssemble,(ex->ss_mat));
789: return(0);
790: }
792: PetscErrorCode MatZeroEntries_HYPRESStruct(Mat mat)
793: {
795: /* before the DMDA is set to the matrix the zero doesn't need to do anything */
796: return(0);
797: }
799: PetscErrorCode MatDestroy_HYPRESStruct(Mat mat)
800: {
801: Mat_HYPRESStruct *ex = (Mat_HYPRESStruct*) mat->data;
802: PetscErrorCode ierr;
803: ISLocalToGlobalMapping ltog;
806: DMGetLocalToGlobalMapping(ex->da,<og);
807: ISLocalToGlobalMappingRestoreIndices(ltog, (const PetscInt **) &ex->gindices);
808: PetscStackCallStandard(HYPRE_SStructGraphDestroy,(ex->ss_graph));
809: PetscStackCallStandard(HYPRE_SStructMatrixDestroy,(ex->ss_mat));
810: PetscStackCallStandard(HYPRE_SStructVectorDestroy,(ex->ss_x));
811: PetscStackCallStandard(HYPRE_SStructVectorDestroy,(ex->ss_b));
812: PetscObjectDereference((PetscObject)ex->da);
813: MPI_Comm_free(&(ex->hcomm));
814: PetscFree(ex);
815: return(0);
816: }
818: PETSC_EXTERN PetscErrorCode MatCreate_HYPRESStruct(Mat B)
819: {
820: Mat_HYPRESStruct *ex;
821: PetscErrorCode ierr;
824: PetscNewLog(B,&ex);
825: B->data = (void*)ex;
826: B->rmap->bs = 1;
827: B->assembled = PETSC_FALSE;
829: B->insertmode = NOT_SET_VALUES;
831: B->ops->assemblyend = MatAssemblyEnd_HYPRESStruct;
832: B->ops->mult = MatMult_HYPRESStruct;
833: B->ops->zeroentries = MatZeroEntries_HYPRESStruct;
834: B->ops->destroy = MatDestroy_HYPRESStruct;
835: B->ops->setup = MatSetUp_HYPRESStruct;
837: ex->needsinitialization = PETSC_TRUE;
839: MPI_Comm_dup(PetscObjectComm((PetscObject)B),&(ex->hcomm));
840: PetscObjectChangeTypeName((PetscObject)B,MATHYPRESSTRUCT);
841: return(0);
842: }