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: 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 SETERRQ(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] = (HYPRE_Int)(ex->xs + (row % ex->nx));
57: index[1] = (HYPRE_Int)(ex->ys + ((row/ex->nx) % ex->ny));
58: index[2] = (HYPRE_Int)(ex->zs + (row/(ex->nxny)));
59: if (addv == ADD_VALUES) {
60: PetscStackCallStandard(HYPRE_StructMatrixAddToValues,ex->hmat,index,ncol,entries,values);
61: } else {
62: PetscStackCallStandard(HYPRE_StructMatrixSetValues,ex->hmat,index,ncol,entries,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: HYPRE_Int index[3],entries[7] = {0,1,2,3,4,5,6};
72: PetscInt row,i;
73: HYPRE_Complex values[7];
74: Mat_HYPREStruct *ex = (Mat_HYPREStruct*) mat->data;
77: PetscArrayzero(values,7);
78: PetscHYPREScalarCast(d,&values[3]);
79: for (i=0; i<nrow; i++) {
80: row = ex->gindices[irow[i]] - ex->rstart;
81: index[0] = (HYPRE_Int)(ex->xs + (row % ex->nx));
82: index[1] = (HYPRE_Int)(ex->ys + ((row/ex->nx) % ex->ny));
83: index[2] = (HYPRE_Int)(ex->zs + (row/(ex->nxny)));
84: PetscStackCallStandard(HYPRE_StructMatrixSetValues,ex->hmat,index,7,entries,values);
85: }
86: PetscStackCallStandard(HYPRE_StructMatrixAssemble,ex->hmat);
87: return 0;
88: }
90: PetscErrorCode MatZeroEntries_HYPREStruct_3d(Mat mat)
91: {
92: HYPRE_Int indices[7] = {0,1,2,3,4,5,6};
93: Mat_HYPREStruct *ex = (Mat_HYPREStruct*) mat->data;
95: /* hypre has no public interface to do this */
96: PetscStackCallStandard(hypre_StructMatrixClearBoxValues,ex->hmat,&ex->hbox,7,indices,0,1);
97: PetscStackCallStandard(HYPRE_StructMatrixAssemble,ex->hmat);
98: return 0;
99: }
101: static PetscErrorCode MatSetUp_HYPREStruct(Mat mat)
102: {
103: Mat_HYPREStruct *ex = (Mat_HYPREStruct*) mat->data;
104: HYPRE_Int sw[6];
105: HYPRE_Int hlower[3],hupper[3];
106: PetscInt dim,dof,psw,nx,ny,nz,ilower[3],iupper[3],ssize,i;
107: DMBoundaryType px,py,pz;
108: DMDAStencilType st;
109: ISLocalToGlobalMapping ltog;
110: DM da;
112: MatGetDM(mat,(DM*)&da);
113: ex->da = da;
114: PetscObjectReference((PetscObject)da);
116: DMDAGetInfo(ex->da,&dim,0,0,0,0,0,0,&dof,&psw,&px,&py,&pz,&st);
117: DMDAGetCorners(ex->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);
119: /* when HYPRE_MIXEDINT is defined, sizeof(HYPRE_Int) == 32 */
120: iupper[0] += ilower[0] - 1;
121: iupper[1] += ilower[1] - 1;
122: iupper[2] += ilower[2] - 1;
123: hlower[0] = (HYPRE_Int)ilower[0];
124: hlower[1] = (HYPRE_Int)ilower[1];
125: hlower[2] = (HYPRE_Int)ilower[2];
126: hupper[0] = (HYPRE_Int)iupper[0];
127: hupper[1] = (HYPRE_Int)iupper[1];
128: hupper[2] = (HYPRE_Int)iupper[2];
130: /* the hypre_Box is used to zero out the matrix entries in MatZeroValues() */
131: ex->hbox.imin[0] = hlower[0];
132: ex->hbox.imin[1] = hlower[1];
133: ex->hbox.imin[2] = hlower[2];
134: ex->hbox.imax[0] = hupper[0];
135: ex->hbox.imax[1] = hupper[1];
136: ex->hbox.imax[2] = hupper[2];
138: /* create the hypre grid object and set its information */
141: PetscStackCallStandard(HYPRE_StructGridCreate,ex->hcomm,dim,&ex->hgrid);
142: PetscStackCallStandard(HYPRE_StructGridSetExtents,ex->hgrid,hlower,hupper);
143: PetscStackCallStandard(HYPRE_StructGridAssemble,ex->hgrid);
145: sw[5] = sw[4] = sw[3] = sw[2] = sw[1] = sw[0] = psw;
146: PetscStackCallStandard(HYPRE_StructGridSetNumGhost,ex->hgrid,sw);
148: /* create the hypre stencil object and set its information */
151: if (dim == 1) {
152: HYPRE_Int offsets[3][1] = {{-1},{0},{1}};
153: ssize = 3;
154: PetscStackCallStandard(HYPRE_StructStencilCreate,dim,ssize,&ex->hstencil);
155: for (i=0; i<ssize; i++) {
156: PetscStackCallStandard(HYPRE_StructStencilSetElement,ex->hstencil,i,offsets[i]);
157: }
158: } else if (dim == 2) {
159: HYPRE_Int offsets[5][2] = {{0,-1},{-1,0},{0,0},{1,0},{0,1}};
160: ssize = 5;
161: PetscStackCallStandard(HYPRE_StructStencilCreate,dim,ssize,&ex->hstencil);
162: for (i=0; i<ssize; i++) {
163: PetscStackCallStandard(HYPRE_StructStencilSetElement,ex->hstencil,i,offsets[i]);
164: }
165: } else if (dim == 3) {
166: 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}};
167: ssize = 7;
168: PetscStackCallStandard(HYPRE_StructStencilCreate,dim,ssize,&ex->hstencil);
169: for (i=0; i<ssize; i++) {
170: PetscStackCallStandard(HYPRE_StructStencilSetElement,ex->hstencil,i,offsets[i]);
171: }
172: }
174: /* create the HYPRE vector for rhs and solution */
175: PetscStackCallStandard(HYPRE_StructVectorCreate,ex->hcomm,ex->hgrid,&ex->hb);
176: PetscStackCallStandard(HYPRE_StructVectorCreate,ex->hcomm,ex->hgrid,&ex->hx);
177: PetscStackCallStandard(HYPRE_StructVectorInitialize,ex->hb);
178: PetscStackCallStandard(HYPRE_StructVectorInitialize,ex->hx);
179: PetscStackCallStandard(HYPRE_StructVectorAssemble,ex->hb);
180: PetscStackCallStandard(HYPRE_StructVectorAssemble,ex->hx);
182: /* create the hypre matrix object and set its information */
183: PetscStackCallStandard(HYPRE_StructMatrixCreate,ex->hcomm,ex->hgrid,ex->hstencil,&ex->hmat);
184: PetscStackCallStandard(HYPRE_StructGridDestroy,ex->hgrid);
185: PetscStackCallStandard(HYPRE_StructStencilDestroy,ex->hstencil);
186: if (ex->needsinitialization) {
187: PetscStackCallStandard(HYPRE_StructMatrixInitialize,ex->hmat);
188: ex->needsinitialization = PETSC_FALSE;
189: }
191: /* set the global and local sizes of the matrix */
192: DMDAGetCorners(da,0,0,0,&nx,&ny,&nz);
193: MatSetSizes(mat,dof*nx*ny*nz,dof*nx*ny*nz,PETSC_DECIDE,PETSC_DECIDE);
194: PetscLayoutSetBlockSize(mat->rmap,1);
195: PetscLayoutSetBlockSize(mat->cmap,1);
196: PetscLayoutSetUp(mat->rmap);
197: PetscLayoutSetUp(mat->cmap);
198: mat->preallocated = PETSC_TRUE;
200: if (dim == 3) {
201: mat->ops->setvalueslocal = MatSetValuesLocal_HYPREStruct_3d;
202: mat->ops->zerorowslocal = MatZeroRowsLocal_HYPREStruct_3d;
203: mat->ops->zeroentries = MatZeroEntries_HYPREStruct_3d;
205: /* MatZeroEntries_HYPREStruct_3d(mat); */
206: } else SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Only support for 3d DMDA currently");
208: /* get values that will be used repeatedly in MatSetValuesLocal() and MatZeroRowsLocal() repeatedly */
209: MatGetOwnershipRange(mat,&ex->rstart,NULL);
210: DMGetLocalToGlobalMapping(ex->da,<og);
211: ISLocalToGlobalMappingGetIndices(ltog, (const PetscInt **) &ex->gindices);
212: DMDAGetGhostCorners(ex->da,0,0,0,&ex->gnx,&ex->gnxgny,0);
213: ex->gnxgny *= ex->gnx;
214: DMDAGetCorners(ex->da,&ex->xs,&ex->ys,&ex->zs,&ex->nx,&ex->ny,0);
215: ex->nxny = ex->nx*ex->ny;
216: return 0;
217: }
219: PetscErrorCode MatMult_HYPREStruct(Mat A,Vec x,Vec y)
220: {
221: const PetscScalar *xx;
222: PetscScalar *yy;
223: PetscInt ilower[3],iupper[3];
224: HYPRE_Int hlower[3],hupper[3];
225: Mat_HYPREStruct *mx = (Mat_HYPREStruct*)(A->data);
227: DMDAGetCorners(mx->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);
228: /* when HYPRE_MIXEDINT is defined, sizeof(HYPRE_Int) == 32 */
229: iupper[0] += ilower[0] - 1;
230: iupper[1] += ilower[1] - 1;
231: iupper[2] += ilower[2] - 1;
232: hlower[0] = (HYPRE_Int)ilower[0];
233: hlower[1] = (HYPRE_Int)ilower[1];
234: hlower[2] = (HYPRE_Int)ilower[2];
235: hupper[0] = (HYPRE_Int)iupper[0];
236: hupper[1] = (HYPRE_Int)iupper[1];
237: hupper[2] = (HYPRE_Int)iupper[2];
239: /* copy x values over to hypre */
240: PetscStackCallStandard(HYPRE_StructVectorSetConstantValues,mx->hb,0.0);
241: VecGetArrayRead(x,&xx);
242: PetscStackCallStandard(HYPRE_StructVectorSetBoxValues,mx->hb,hlower,hupper,(HYPRE_Complex*)xx);
243: VecRestoreArrayRead(x,&xx);
244: PetscStackCallStandard(HYPRE_StructVectorAssemble,mx->hb);
245: PetscStackCallStandard(HYPRE_StructMatrixMatvec,1.0,mx->hmat,mx->hb,0.0,mx->hx);
247: /* copy solution values back to PETSc */
248: VecGetArray(y,&yy);
249: PetscStackCallStandard(HYPRE_StructVectorGetBoxValues,mx->hx,hlower,hupper,(HYPRE_Complex*)yy);
250: VecRestoreArray(y,&yy);
251: return 0;
252: }
254: PetscErrorCode MatAssemblyEnd_HYPREStruct(Mat mat,MatAssemblyType mode)
255: {
256: Mat_HYPREStruct *ex = (Mat_HYPREStruct*) mat->data;
258: PetscStackCallStandard(HYPRE_StructMatrixAssemble,ex->hmat);
259: /* PetscStackCallStandard(HYPRE_StructMatrixPrint,"dummy",ex->hmat,0); */
260: return 0;
261: }
263: PetscErrorCode MatZeroEntries_HYPREStruct(Mat mat)
264: {
265: /* before the DMDA is set to the matrix the zero doesn't need to do anything */
266: return 0;
267: }
269: PetscErrorCode MatDestroy_HYPREStruct(Mat mat)
270: {
271: Mat_HYPREStruct *ex = (Mat_HYPREStruct*) mat->data;
273: PetscStackCallStandard(HYPRE_StructMatrixDestroy,ex->hmat);
274: PetscStackCallStandard(HYPRE_StructVectorDestroy,ex->hx);
275: PetscStackCallStandard(HYPRE_StructVectorDestroy,ex->hb);
276: PetscObjectDereference((PetscObject)ex->da);
277: MPI_Comm_free(&(ex->hcomm));
278: PetscFree(ex);
279: return 0;
280: }
282: PETSC_EXTERN PetscErrorCode MatCreate_HYPREStruct(Mat B)
283: {
284: Mat_HYPREStruct *ex;
286: PetscNewLog(B,&ex);
287: B->data = (void*)ex;
288: B->rmap->bs = 1;
289: B->assembled = PETSC_FALSE;
291: B->insertmode = NOT_SET_VALUES;
293: B->ops->assemblyend = MatAssemblyEnd_HYPREStruct;
294: B->ops->mult = MatMult_HYPREStruct;
295: B->ops->zeroentries = MatZeroEntries_HYPREStruct;
296: B->ops->destroy = MatDestroy_HYPREStruct;
297: B->ops->setup = MatSetUp_HYPREStruct;
299: ex->needsinitialization = PETSC_TRUE;
301: MPI_Comm_dup(PetscObjectComm((PetscObject)B),&(ex->hcomm));
302: PetscObjectChangeTypeName((PetscObject)B,MATHYPRESTRUCT);
303: return 0;
304: }
306: /*MC
307: MATHYPRESSTRUCT - MATHYPRESSTRUCT = "hypresstruct" - A matrix type to be used for parallel sparse matrices
308: based on the hypre HYPRE_SStructMatrix.
310: Level: intermediate
312: Notes:
313: Unlike hypre's general semi-struct object consisting of a collection of structured-grid objects and unstructured
314: grid objects, we restrict the semi-struct objects to consist of only structured-grid components.
316: 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
317: be defined by a DMDA.
319: The matrix needs a DMDA associated with it by either a call to MatSetDM() or if the matrix is obtained from DMCreateMatrix()
321: M*/
323: PetscErrorCode MatSetValuesLocal_HYPRESStruct_3d(Mat mat,PetscInt nrow,const PetscInt irow[],PetscInt ncol,const PetscInt icol[],const PetscScalar y[],InsertMode addv)
324: {
325: HYPRE_Int index[3],*entries;
326: PetscInt i,j,stencil;
327: HYPRE_Complex *values = (HYPRE_Complex*)y;
328: Mat_HYPRESStruct *ex = (Mat_HYPRESStruct*) mat->data;
330: PetscInt part = 0; /* Petsc sstruct interface only allows 1 part */
331: PetscInt ordering;
332: PetscInt grid_rank, to_grid_rank;
333: PetscInt var_type, to_var_type;
334: PetscInt to_var_entry = 0;
335: PetscInt nvars= ex->nvars;
336: PetscInt row;
338: PetscMalloc1(7*nvars,&entries);
339: ordering = ex->dofs_order; /* ordering= 0 nodal ordering
340: 1 variable ordering */
341: /* stencil entries are ordered by variables: var0_stencil0, var0_stencil1, ..., var0_stencil6, var1_stencil0, var1_stencil1, ... */
342: if (!ordering) {
343: for (i=0; i<nrow; i++) {
344: grid_rank= irow[i]/nvars;
345: var_type = (irow[i] % nvars);
347: for (j=0; j<ncol; j++) {
348: to_grid_rank= icol[j]/nvars;
349: to_var_type = (icol[j] % nvars);
351: to_var_entry= to_var_entry*7;
352: entries[j] = to_var_entry;
354: stencil = to_grid_rank-grid_rank;
355: if (!stencil) {
356: entries[j] += 3;
357: } else if (stencil == -1) {
358: entries[j] += 2;
359: } else if (stencil == 1) {
360: entries[j] += 4;
361: } else if (stencil == -ex->gnx) {
362: entries[j] += 1;
363: } else if (stencil == ex->gnx) {
364: entries[j] += 5;
365: } else if (stencil == -ex->gnxgny) {
366: entries[j] += 0;
367: } else if (stencil == ex->gnxgny) {
368: entries[j] += 6;
369: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Local row %D local column %D have bad stencil %D",irow[i],icol[j],stencil);
370: }
372: row = ex->gindices[grid_rank] - ex->rstart;
373: index[0] = (HYPRE_Int)(ex->xs + (row % ex->nx));
374: index[1] = (HYPRE_Int)(ex->ys + ((row/ex->nx) % ex->ny));
375: index[2] = (HYPRE_Int)(ex->zs + (row/(ex->nxny)));
377: if (addv == ADD_VALUES) {
378: PetscStackCallStandard(HYPRE_SStructMatrixAddToValues,ex->ss_mat,part,index,var_type,ncol,entries,values);
379: } else {
380: PetscStackCallStandard(HYPRE_SStructMatrixSetValues,ex->ss_mat,part,index,var_type,ncol,entries,values);
381: }
382: values += ncol;
383: }
384: } else {
385: for (i=0; i<nrow; i++) {
386: var_type = irow[i]/(ex->gnxgnygnz);
387: grid_rank= irow[i] - var_type*(ex->gnxgnygnz);
389: for (j=0; j<ncol; j++) {
390: to_var_type = icol[j]/(ex->gnxgnygnz);
391: to_grid_rank= icol[j] - to_var_type*(ex->gnxgnygnz);
393: to_var_entry= to_var_entry*7;
394: entries[j] = to_var_entry;
396: stencil = to_grid_rank-grid_rank;
397: if (!stencil) {
398: entries[j] += 3;
399: } else if (stencil == -1) {
400: entries[j] += 2;
401: } else if (stencil == 1) {
402: entries[j] += 4;
403: } else if (stencil == -ex->gnx) {
404: entries[j] += 1;
405: } else if (stencil == ex->gnx) {
406: entries[j] += 5;
407: } else if (stencil == -ex->gnxgny) {
408: entries[j] += 0;
409: } else if (stencil == ex->gnxgny) {
410: entries[j] += 6;
411: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Local row %D local column %D have bad stencil %D",irow[i],icol[j],stencil);
412: }
414: row = ex->gindices[grid_rank] - ex->rstart;
415: index[0] = (HYPRE_Int)(ex->xs + (row % ex->nx));
416: index[1] = (HYPRE_Int)(ex->ys + ((row/ex->nx) % ex->ny));
417: index[2] = (HYPRE_Int)(ex->zs + (row/(ex->nxny)));
419: if (addv == ADD_VALUES) {
420: PetscStackCallStandard(HYPRE_SStructMatrixAddToValues,ex->ss_mat,part,index,var_type,ncol,entries,values);
421: } else {
422: PetscStackCallStandard(HYPRE_SStructMatrixSetValues,ex->ss_mat,part,index,var_type,ncol,entries,values);
423: }
424: values += ncol;
425: }
426: }
427: PetscFree(entries);
428: return 0;
429: }
431: PetscErrorCode MatZeroRowsLocal_HYPRESStruct_3d(Mat mat,PetscInt nrow,const PetscInt irow[],PetscScalar d,Vec x,Vec b)
432: {
433: HYPRE_Int index[3],*entries;
434: PetscInt i;
435: HYPRE_Complex **values;
436: Mat_HYPRESStruct *ex = (Mat_HYPRESStruct*) mat->data;
438: PetscInt part = 0; /* Petsc sstruct interface only allows 1 part */
439: PetscInt ordering = ex->dofs_order;
440: PetscInt grid_rank;
441: PetscInt var_type;
442: PetscInt nvars= ex->nvars;
443: PetscInt row;
446: PetscMalloc1(7*nvars,&entries);
448: PetscMalloc1(nvars,&values);
449: PetscMalloc1(7*nvars*nvars,&values[0]);
450: for (i=1; i<nvars; i++) {
451: values[i] = values[i-1] + nvars*7;
452: }
454: for (i=0; i< nvars; i++) {
455: PetscArrayzero(values[i],nvars*7*sizeof(HYPRE_Complex));
456: PetscHYPREScalarCast(d,values[i]+3);
457: }
459: for (i=0; i< nvars*7; i++) entries[i] = i;
461: if (!ordering) {
462: for (i=0; i<nrow; i++) {
463: grid_rank = irow[i] / nvars;
464: var_type = (irow[i] % nvars);
466: row = ex->gindices[grid_rank] - ex->rstart;
467: index[0] = (HYPRE_Int)(ex->xs + (row % ex->nx));
468: index[1] = (HYPRE_Int)(ex->ys + ((row/ex->nx) % ex->ny));
469: index[2] = (HYPRE_Int)(ex->zs + (row/(ex->nxny)));
470: PetscStackCallStandard(HYPRE_SStructMatrixSetValues,ex->ss_mat,part,index,var_type,7*nvars,entries,values[var_type]);
471: }
472: } else {
473: for (i=0; i<nrow; i++) {
474: var_type = irow[i]/(ex->gnxgnygnz);
475: grid_rank= irow[i] - var_type*(ex->gnxgnygnz);
477: row = ex->gindices[grid_rank] - ex->rstart;
478: index[0] = (HYPRE_Int)(ex->xs + (row % ex->nx));
479: index[1] = (HYPRE_Int)(ex->ys + ((row/ex->nx) % ex->ny));
480: index[2] = (HYPRE_Int)(ex->zs + (row/(ex->nxny)));
481: PetscStackCallStandard(HYPRE_SStructMatrixSetValues,ex->ss_mat,part,index,var_type,7*nvars,entries,values[var_type]);
482: }
483: }
484: PetscStackCallStandard(HYPRE_SStructMatrixAssemble,ex->ss_mat);
485: PetscFree(values[0]);
486: PetscFree(values);
487: PetscFree(entries);
488: return 0;
489: }
491: PetscErrorCode MatZeroEntries_HYPRESStruct_3d(Mat mat)
492: {
493: Mat_HYPRESStruct *ex = (Mat_HYPRESStruct*) mat->data;
494: PetscInt nvars= ex->nvars;
495: PetscInt size;
496: PetscInt part= 0; /* only one part */
498: 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);
499: {
500: HYPRE_Int i,*entries,iupper[3],ilower[3];
501: HYPRE_Complex *values;
503: for (i= 0; i< 3; i++) {
504: ilower[i] = ex->hbox.imin[i];
505: iupper[i] = ex->hbox.imax[i];
506: }
508: PetscMalloc2(nvars*7,&entries,nvars*7*size,&values);
509: for (i= 0; i< nvars*7; i++) entries[i] = i;
510: PetscArrayzero(values,nvars*7*size);
512: for (i= 0; i< nvars; i++) {
513: PetscStackCallStandard(HYPRE_SStructMatrixSetBoxValues,ex->ss_mat,part,ilower,iupper,i,nvars*7,entries,values);
514: }
515: PetscFree2(entries,values);
516: }
517: PetscStackCallStandard(HYPRE_SStructMatrixAssemble,ex->ss_mat);
518: return 0;
519: }
521: static PetscErrorCode MatSetUp_HYPRESStruct(Mat mat)
522: {
523: Mat_HYPRESStruct *ex = (Mat_HYPRESStruct*) mat->data;
524: PetscInt dim,dof,sw[3],nx,ny,nz;
525: PetscInt ilower[3],iupper[3],ssize,i;
526: DMBoundaryType px,py,pz;
527: DMDAStencilType st;
528: PetscInt nparts= 1; /* assuming only one part */
529: PetscInt part = 0;
530: ISLocalToGlobalMapping ltog;
531: DM da;
533: MatGetDM(mat,(DM*)&da);
534: ex->da = da;
535: PetscObjectReference((PetscObject)da);
537: DMDAGetInfo(ex->da,&dim,0,0,0,0,0,0,&dof,&sw[0],&px,&py,&pz,&st);
538: DMDAGetCorners(ex->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);
539: iupper[0] += ilower[0] - 1;
540: iupper[1] += ilower[1] - 1;
541: iupper[2] += ilower[2] - 1;
542: /* the hypre_Box is used to zero out the matrix entries in MatZeroValues() */
543: ex->hbox.imin[0] = ilower[0];
544: ex->hbox.imin[1] = ilower[1];
545: ex->hbox.imin[2] = ilower[2];
546: ex->hbox.imax[0] = iupper[0];
547: ex->hbox.imax[1] = iupper[1];
548: ex->hbox.imax[2] = iupper[2];
550: ex->dofs_order = 0;
552: /* assuming that the same number of dofs on each gridpoint. Also assume all cell-centred based */
553: ex->nvars= dof;
555: /* create the hypre grid object and set its information */
557: PetscStackCallStandard(HYPRE_SStructGridCreate,ex->hcomm,dim,nparts,&ex->ss_grid);
558: PetscStackCallStandard(HYPRE_SStructGridSetExtents,ex->ss_grid,part,ex->hbox.imin,ex->hbox.imax);
559: {
560: HYPRE_SStructVariable *vartypes;
561: PetscMalloc1(ex->nvars,&vartypes);
562: for (i= 0; i< ex->nvars; i++) vartypes[i]= HYPRE_SSTRUCT_VARIABLE_CELL;
563: PetscStackCallStandard(HYPRE_SStructGridSetVariables,ex->ss_grid, part, ex->nvars,vartypes);
564: PetscFree(vartypes);
565: }
566: PetscStackCallStandard(HYPRE_SStructGridAssemble,ex->ss_grid);
568: sw[1] = sw[0];
569: sw[2] = sw[1];
570: /* PetscStackCallStandard(HYPRE_SStructGridSetNumGhost,ex->ss_grid,sw); */
572: /* create the hypre stencil object and set its information */
576: if (dim == 1) {
577: HYPRE_Int offsets[3][1] = {{-1},{0},{1}};
578: PetscInt j, cnt;
580: ssize = 3*(ex->nvars);
581: PetscStackCallStandard(HYPRE_SStructStencilCreate,dim,ssize,&ex->ss_stencil);
582: cnt= 0;
583: for (i = 0; i < (ex->nvars); i++) {
584: for (j = 0; j < 3; j++) {
585: PetscStackCallStandard(HYPRE_SStructStencilSetEntry,ex->ss_stencil, cnt, offsets[j], i);
586: cnt++;
587: }
588: }
589: } else if (dim == 2) {
590: HYPRE_Int offsets[5][2] = {{0,-1},{-1,0},{0,0},{1,0},{0,1}};
591: PetscInt j, cnt;
593: ssize = 5*(ex->nvars);
594: PetscStackCallStandard(HYPRE_SStructStencilCreate,dim,ssize,&ex->ss_stencil);
595: cnt= 0;
596: for (i= 0; i< (ex->nvars); i++) {
597: for (j= 0; j< 5; j++) {
598: PetscStackCallStandard(HYPRE_SStructStencilSetEntry,ex->ss_stencil, cnt, offsets[j], i);
599: cnt++;
600: }
601: }
602: } else if (dim == 3) {
603: 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}};
604: PetscInt j, cnt;
606: ssize = 7*(ex->nvars);
607: PetscStackCallStandard(HYPRE_SStructStencilCreate,dim,ssize,&ex->ss_stencil);
608: cnt= 0;
609: for (i= 0; i< (ex->nvars); i++) {
610: for (j= 0; j< 7; j++) {
611: PetscStackCallStandard(HYPRE_SStructStencilSetEntry,ex->ss_stencil, cnt, offsets[j], i);
612: cnt++;
613: }
614: }
615: }
617: /* create the HYPRE graph */
618: PetscStackCallStandard(HYPRE_SStructGraphCreate,ex->hcomm, ex->ss_grid, &(ex->ss_graph));
620: /* set the stencil graph. Note that each variable has the same graph. This means that each
621: variable couples to all the other variable and with the same stencil pattern. */
622: for (i= 0; i< (ex->nvars); i++) {
623: PetscStackCallStandard(HYPRE_SStructGraphSetStencil,ex->ss_graph,part,i,ex->ss_stencil);
624: }
625: PetscStackCallStandard(HYPRE_SStructGraphAssemble,ex->ss_graph);
627: /* create the HYPRE sstruct vectors for rhs and solution */
628: PetscStackCallStandard(HYPRE_SStructVectorCreate,ex->hcomm,ex->ss_grid,&ex->ss_b);
629: PetscStackCallStandard(HYPRE_SStructVectorCreate,ex->hcomm,ex->ss_grid,&ex->ss_x);
630: PetscStackCallStandard(HYPRE_SStructVectorInitialize,ex->ss_b);
631: PetscStackCallStandard(HYPRE_SStructVectorInitialize,ex->ss_x);
632: PetscStackCallStandard(HYPRE_SStructVectorAssemble,ex->ss_b);
633: PetscStackCallStandard(HYPRE_SStructVectorAssemble,ex->ss_x);
635: /* create the hypre matrix object and set its information */
636: PetscStackCallStandard(HYPRE_SStructMatrixCreate,ex->hcomm,ex->ss_graph,&ex->ss_mat);
637: PetscStackCallStandard(HYPRE_SStructGridDestroy,ex->ss_grid);
638: PetscStackCallStandard(HYPRE_SStructStencilDestroy,ex->ss_stencil);
639: if (ex->needsinitialization) {
640: PetscStackCallStandard(HYPRE_SStructMatrixInitialize,ex->ss_mat);
641: ex->needsinitialization = PETSC_FALSE;
642: }
644: /* set the global and local sizes of the matrix */
645: DMDAGetCorners(da,0,0,0,&nx,&ny,&nz);
646: MatSetSizes(mat,dof*nx*ny*nz,dof*nx*ny*nz,PETSC_DECIDE,PETSC_DECIDE);
647: PetscLayoutSetBlockSize(mat->rmap,dof);
648: PetscLayoutSetBlockSize(mat->cmap,dof);
649: PetscLayoutSetUp(mat->rmap);
650: PetscLayoutSetUp(mat->cmap);
651: mat->preallocated = PETSC_TRUE;
653: if (dim == 3) {
654: mat->ops->setvalueslocal = MatSetValuesLocal_HYPRESStruct_3d;
655: mat->ops->zerorowslocal = MatZeroRowsLocal_HYPRESStruct_3d;
656: mat->ops->zeroentries = MatZeroEntries_HYPRESStruct_3d;
658: /* MatZeroEntries_HYPRESStruct_3d(mat); */
659: } else SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Only support for 3d DMDA currently");
661: /* get values that will be used repeatedly in MatSetValuesLocal() and MatZeroRowsLocal() repeatedly */
662: MatGetOwnershipRange(mat,&ex->rstart,NULL);
663: DMGetLocalToGlobalMapping(ex->da,<og);
664: ISLocalToGlobalMappingGetIndices(ltog, (const PetscInt **) &ex->gindices);
665: DMDAGetGhostCorners(ex->da,0,0,0,&ex->gnx,&ex->gnxgny,&ex->gnxgnygnz);
667: ex->gnxgny *= ex->gnx;
668: ex->gnxgnygnz *= ex->gnxgny;
670: DMDAGetCorners(ex->da,&ex->xs,&ex->ys,&ex->zs,&ex->nx,&ex->ny,&ex->nz);
672: ex->nxny = ex->nx*ex->ny;
673: ex->nxnynz = ex->nz*ex->nxny;
674: return 0;
675: }
677: PetscErrorCode MatMult_HYPRESStruct(Mat A,Vec x,Vec y)
678: {
679: const PetscScalar *xx;
680: PetscScalar *yy;
681: HYPRE_Int hlower[3],hupper[3];
682: PetscInt ilower[3],iupper[3];
683: Mat_HYPRESStruct *mx = (Mat_HYPRESStruct*)(A->data);
684: PetscInt ordering= mx->dofs_order;
685: PetscInt nvars = mx->nvars;
686: PetscInt part = 0;
687: PetscInt size;
688: PetscInt i;
690: DMDAGetCorners(mx->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);
692: /* when HYPRE_MIXEDINT is defined, sizeof(HYPRE_Int) == 32 */
693: iupper[0] += ilower[0] - 1;
694: iupper[1] += ilower[1] - 1;
695: iupper[2] += ilower[2] - 1;
696: hlower[0] = (HYPRE_Int)ilower[0];
697: hlower[1] = (HYPRE_Int)ilower[1];
698: hlower[2] = (HYPRE_Int)ilower[2];
699: hupper[0] = (HYPRE_Int)iupper[0];
700: hupper[1] = (HYPRE_Int)iupper[1];
701: hupper[2] = (HYPRE_Int)iupper[2];
703: size= 1;
704: for (i = 0; i < 3; i++) size *= (iupper[i]-ilower[i]+1);
706: /* copy x values over to hypre for variable ordering */
707: if (ordering) {
708: PetscStackCallStandard(HYPRE_SStructVectorSetConstantValues,mx->ss_b,0.0);
709: VecGetArrayRead(x,&xx);
710: for (i= 0; i< nvars; i++) {
711: PetscStackCallStandard(HYPRE_SStructVectorSetBoxValues,mx->ss_b,part,hlower,hupper,i,(HYPRE_Complex*)(xx+(size*i)));
712: }
713: VecRestoreArrayRead(x,&xx);
714: PetscStackCallStandard(HYPRE_SStructVectorAssemble,mx->ss_b);
715: PetscStackCallStandard(HYPRE_SStructMatrixMatvec,1.0,mx->ss_mat,mx->ss_b,0.0,mx->ss_x);
717: /* copy solution values back to PETSc */
718: VecGetArray(y,&yy);
719: for (i= 0; i< nvars; i++) {
720: PetscStackCallStandard(HYPRE_SStructVectorGetBoxValues,mx->ss_x,part,hlower,hupper,i,(HYPRE_Complex*)(yy+(size*i)));
721: }
722: VecRestoreArray(y,&yy);
723: } else { /* nodal ordering must be mapped to variable ordering for sys_pfmg */
724: PetscScalar *z;
725: PetscInt j, k;
727: PetscMalloc1(nvars*size,&z);
728: PetscStackCallStandard(HYPRE_SStructVectorSetConstantValues,mx->ss_b,0.0);
729: VecGetArrayRead(x,&xx);
731: /* transform nodal to hypre's variable ordering for sys_pfmg */
732: for (i= 0; i< size; i++) {
733: k= i*nvars;
734: for (j= 0; j< nvars; j++) z[j*size+i]= xx[k+j];
735: }
736: for (i= 0; i< nvars; i++) {
737: PetscStackCallStandard(HYPRE_SStructVectorSetBoxValues,mx->ss_b,part,hlower,hupper,i,(HYPRE_Complex*)(z+(size*i)));
738: }
739: VecRestoreArrayRead(x,&xx);
740: PetscStackCallStandard(HYPRE_SStructVectorAssemble,mx->ss_b);
741: PetscStackCallStandard(HYPRE_SStructMatrixMatvec,1.0,mx->ss_mat,mx->ss_b,0.0,mx->ss_x);
743: /* copy solution values back to PETSc */
744: VecGetArray(y,&yy);
745: for (i= 0; i< nvars; i++) {
746: PetscStackCallStandard(HYPRE_SStructVectorGetBoxValues,mx->ss_x,part,hlower,hupper,i,(HYPRE_Complex*)(z+(size*i)));
747: }
748: /* transform hypre's variable ordering for sys_pfmg to nodal ordering */
749: for (i= 0; i< size; i++) {
750: k= i*nvars;
751: for (j= 0; j< nvars; j++) yy[k+j]= z[j*size+i];
752: }
753: VecRestoreArray(y,&yy);
754: PetscFree(z);
755: }
756: return 0;
757: }
759: PetscErrorCode MatAssemblyEnd_HYPRESStruct(Mat mat,MatAssemblyType mode)
760: {
761: Mat_HYPRESStruct *ex = (Mat_HYPRESStruct*) mat->data;
763: PetscStackCallStandard(HYPRE_SStructMatrixAssemble,ex->ss_mat);
764: return 0;
765: }
767: PetscErrorCode MatZeroEntries_HYPRESStruct(Mat mat)
768: {
769: /* before the DMDA is set to the matrix the zero doesn't need to do anything */
770: return 0;
771: }
773: PetscErrorCode MatDestroy_HYPRESStruct(Mat mat)
774: {
775: Mat_HYPRESStruct *ex = (Mat_HYPRESStruct*) mat->data;
776: ISLocalToGlobalMapping ltog;
778: DMGetLocalToGlobalMapping(ex->da,<og);
779: ISLocalToGlobalMappingRestoreIndices(ltog, (const PetscInt **) &ex->gindices);
780: PetscStackCallStandard(HYPRE_SStructGraphDestroy,ex->ss_graph);
781: PetscStackCallStandard(HYPRE_SStructMatrixDestroy,ex->ss_mat);
782: PetscStackCallStandard(HYPRE_SStructVectorDestroy,ex->ss_x);
783: PetscStackCallStandard(HYPRE_SStructVectorDestroy,ex->ss_b);
784: PetscObjectDereference((PetscObject)ex->da);
785: MPI_Comm_free(&(ex->hcomm));
786: PetscFree(ex);
787: return 0;
788: }
790: PETSC_EXTERN PetscErrorCode MatCreate_HYPRESStruct(Mat B)
791: {
792: Mat_HYPRESStruct *ex;
794: PetscNewLog(B,&ex);
795: B->data = (void*)ex;
796: B->rmap->bs = 1;
797: B->assembled = PETSC_FALSE;
799: B->insertmode = NOT_SET_VALUES;
801: B->ops->assemblyend = MatAssemblyEnd_HYPRESStruct;
802: B->ops->mult = MatMult_HYPRESStruct;
803: B->ops->zeroentries = MatZeroEntries_HYPRESStruct;
804: B->ops->destroy = MatDestroy_HYPRESStruct;
805: B->ops->setup = MatSetUp_HYPRESStruct;
807: ex->needsinitialization = PETSC_TRUE;
809: MPI_Comm_dup(PetscObjectComm((PetscObject)B),&(ex->hcomm));
810: PetscObjectChangeTypeName((PetscObject)B,MATHYPRESSTRUCT);
811: return 0;
812: }