Actual source code: mhyp.c
petsc-3.9.4 2018-09-11
2: /*
3: Creates hypre ijmatrix from PETSc matrix
4: */
5: #include <petscsys.h>
6: #include <petsc/private/matimpl.h>
7: #include <petscdmda.h>
8: #include <../src/dm/impls/da/hypre/mhyp.h>
10: /* -----------------------------------------------------------------------------------------------------------------*/
12: /*MC
13: MATHYPRESTRUCT - MATHYPRESTRUCT = "hyprestruct" - A matrix type to be used for parallel sparse matrices
14: based on the hypre HYPRE_StructMatrix.
16: Level: intermediate
18: Notes: Unlike the more general support for blocks in hypre this allows only one block per process and requires the block
19: be defined by a DMDA.
21: The matrix needs a DMDA associated with it by either a call to MatSetDM() or if the matrix is obtained from DMCreateMatrix()
23: .seealso: MatCreate(), PCPFMG, MatSetDM(), DMCreateMatrix()
24: M*/
27: PetscErrorCode MatSetValuesLocal_HYPREStruct_3d(Mat mat,PetscInt nrow,const PetscInt irow[],PetscInt ncol,const PetscInt icol[],const PetscScalar y[],InsertMode addv)
28: {
29: PetscErrorCode ierr;
30: PetscInt i,j,stencil,index[3],row,entries[7];
31: const PetscScalar *values = y;
32: Mat_HYPREStruct *ex = (Mat_HYPREStruct*) mat->data;
35: if (PetscUnlikely(ncol > 7)) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"ncol %D > 7 too large",ncol);
36: for (i=0; i<nrow; i++) {
37: for (j=0; j<ncol; j++) {
38: stencil = icol[j] - irow[i];
39: if (!stencil) {
40: entries[j] = 3;
41: } else if (stencil == -1) {
42: entries[j] = 2;
43: } else if (stencil == 1) {
44: entries[j] = 4;
45: } else if (stencil == -ex->gnx) {
46: entries[j] = 1;
47: } else if (stencil == ex->gnx) {
48: entries[j] = 5;
49: } else if (stencil == -ex->gnxgny) {
50: entries[j] = 0;
51: } else if (stencil == ex->gnxgny) {
52: entries[j] = 6;
53: } else SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Local row %D local column %D have bad stencil %D",irow[i],icol[j],stencil);
54: }
55: row = ex->gindices[irow[i]] - ex->rstart;
56: index[0] = ex->xs + (row % ex->nx);
57: index[1] = ex->ys + ((row/ex->nx) % ex->ny);
58: index[2] = ex->zs + (row/(ex->nxny));
59: if (addv == ADD_VALUES) {
60: PetscStackCallStandard(HYPRE_StructMatrixAddToValues,(ex->hmat,(HYPRE_Int *)index,ncol,(HYPRE_Int *)entries,(PetscScalar*)values));
61: } else {
62: PetscStackCallStandard(HYPRE_StructMatrixSetValues,(ex->hmat,(HYPRE_Int *)index,ncol,(HYPRE_Int *)entries,(PetscScalar*)values));
63: }
64: values += ncol;
65: }
66: return(0);
67: }
69: PetscErrorCode MatZeroRowsLocal_HYPREStruct_3d(Mat mat,PetscInt nrow,const PetscInt irow[],PetscScalar d,Vec x,Vec b)
70: {
71: PetscErrorCode ierr;
72: PetscInt i,index[3],row,entries[7] = {0,1,2,3,4,5,6};
73: PetscScalar values[7];
74: Mat_HYPREStruct *ex = (Mat_HYPREStruct*) mat->data;
77: if (x && b) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"No support");
78: PetscMemzero(values,7*sizeof(PetscScalar));
79: values[3] = d;
80: for (i=0; i<nrow; i++) {
81: row = ex->gindices[irow[i]] - ex->rstart;
82: index[0] = ex->xs + (row % ex->nx);
83: index[1] = ex->ys + ((row/ex->nx) % ex->ny);
84: index[2] = ex->zs + (row/(ex->nxny));
85: PetscStackCallStandard(HYPRE_StructMatrixSetValues,(ex->hmat,(HYPRE_Int *)index,7,(HYPRE_Int *)entries,values));
86: }
87: PetscStackCallStandard(HYPRE_StructMatrixAssemble,(ex->hmat));
88: return(0);
89: }
91: PetscErrorCode MatZeroEntries_HYPREStruct_3d(Mat mat)
92: {
93: PetscErrorCode ierr;
94: PetscInt indices[7] = {0,1,2,3,4,5,6};
95: Mat_HYPREStruct *ex = (Mat_HYPREStruct*) mat->data;
98: /* hypre has no public interface to do this */
99: PetscStackCallStandard(hypre_StructMatrixClearBoxValues,(ex->hmat,&ex->hbox,7,(HYPRE_Int *)indices,0,1));
100: PetscStackCallStandard(HYPRE_StructMatrixAssemble,(ex->hmat));
101: return(0);
102: }
104: static PetscErrorCode MatSetUp_HYPREStruct(Mat mat)
105: {
106: PetscErrorCode ierr;
107: Mat_HYPREStruct *ex = (Mat_HYPREStruct*) mat->data;
108: PetscInt dim,dof,sw[6],nx,ny,nz,ilower[3],iupper[3],ssize,i;
109: DMBoundaryType px,py,pz;
110: DMDAStencilType st;
111: ISLocalToGlobalMapping ltog;
112: DM da;
115: MatGetDM(mat,(DM*)&da);
116: ex->da = da;
117: PetscObjectReference((PetscObject)da);
119: DMDAGetInfo(ex->da,&dim,0,0,0,0,0,0,&dof,&sw[0],&px,&py,&pz,&st);
120: DMDAGetCorners(ex->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);
121: iupper[0] += ilower[0] - 1;
122: iupper[1] += ilower[1] - 1;
123: iupper[2] += ilower[2] - 1;
125: /* the hypre_Box is used to zero out the matrix entries in MatZeroValues() */
126: ex->hbox.imin[0] = ilower[0];
127: ex->hbox.imin[1] = ilower[1];
128: ex->hbox.imin[2] = ilower[2];
129: ex->hbox.imax[0] = iupper[0];
130: ex->hbox.imax[1] = iupper[1];
131: ex->hbox.imax[2] = iupper[2];
133: /* create the hypre grid object and set its information */
134: if (dof > 1) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Currently only support for scalar problems");
135: if (px || py || pz) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Ask us to add periodic support by calling HYPRE_StructGridSetPeriodic()");
136: PetscStackCallStandard(HYPRE_StructGridCreate,(ex->hcomm,dim,&ex->hgrid));
137: PetscStackCallStandard(HYPRE_StructGridSetExtents,(ex->hgrid,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper));
138: PetscStackCallStandard(HYPRE_StructGridAssemble,(ex->hgrid));
140: sw[5] = sw[4] = sw[3] = sw[2] = sw[1] = sw[0];
141: PetscStackCallStandard(HYPRE_StructGridSetNumGhost,(ex->hgrid,(HYPRE_Int *)sw));
143: /* create the hypre stencil object and set its information */
144: if (sw[0] > 1) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Ask us to add support for wider stencils");
145: if (st == DMDA_STENCIL_BOX) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Ask us to add support for box stencils");
146: if (dim == 1) {
147: PetscInt offsets[3][1] = {{-1},{0},{1}};
148: ssize = 3;
149: PetscStackCallStandard(HYPRE_StructStencilCreate,(dim,ssize,&ex->hstencil));
150: for (i=0; i<ssize; i++) {
151: PetscStackCallStandard(HYPRE_StructStencilSetElement,(ex->hstencil,i,(HYPRE_Int *)offsets[i]));
152: }
153: } else if (dim == 2) {
154: PetscInt offsets[5][2] = {{0,-1},{-1,0},{0,0},{1,0},{0,1}};
155: ssize = 5;
156: PetscStackCallStandard(HYPRE_StructStencilCreate,(dim,ssize,&ex->hstencil));
157: for (i=0; i<ssize; i++) {
158: PetscStackCallStandard(HYPRE_StructStencilSetElement,(ex->hstencil,i,(HYPRE_Int *)offsets[i]));
159: }
160: } else if (dim == 3) {
161: PetscInt offsets[7][3] = {{0,0,-1},{0,-1,0},{-1,0,0},{0,0,0},{1,0,0},{0,1,0},{0,0,1}};
162: ssize = 7;
163: PetscStackCallStandard(HYPRE_StructStencilCreate,(dim,ssize,&ex->hstencil));
164: for (i=0; i<ssize; i++) {
165: PetscStackCallStandard(HYPRE_StructStencilSetElement,(ex->hstencil,i,(HYPRE_Int *)offsets[i]));
166: }
167: }
169: /* create the HYPRE vector for rhs and solution */
170: PetscStackCallStandard(HYPRE_StructVectorCreate,(ex->hcomm,ex->hgrid,&ex->hb));
171: PetscStackCallStandard(HYPRE_StructVectorCreate,(ex->hcomm,ex->hgrid,&ex->hx));
172: PetscStackCallStandard(HYPRE_StructVectorInitialize,(ex->hb));
173: PetscStackCallStandard(HYPRE_StructVectorInitialize,(ex->hx));
174: PetscStackCallStandard(HYPRE_StructVectorAssemble,(ex->hb));
175: PetscStackCallStandard(HYPRE_StructVectorAssemble,(ex->hx));
177: /* create the hypre matrix object and set its information */
178: PetscStackCallStandard(HYPRE_StructMatrixCreate,(ex->hcomm,ex->hgrid,ex->hstencil,&ex->hmat));
179: PetscStackCallStandard(HYPRE_StructGridDestroy,(ex->hgrid));
180: PetscStackCallStandard(HYPRE_StructStencilDestroy,(ex->hstencil));
181: if (ex->needsinitialization) {
182: PetscStackCallStandard(HYPRE_StructMatrixInitialize,(ex->hmat));
183: ex->needsinitialization = PETSC_FALSE;
184: }
186: /* set the global and local sizes of the matrix */
187: DMDAGetCorners(da,0,0,0,&nx,&ny,&nz);
188: MatSetSizes(mat,dof*nx*ny*nz,dof*nx*ny*nz,PETSC_DECIDE,PETSC_DECIDE);
189: PetscLayoutSetBlockSize(mat->rmap,1);
190: PetscLayoutSetBlockSize(mat->cmap,1);
191: PetscLayoutSetUp(mat->rmap);
192: PetscLayoutSetUp(mat->cmap);
193: mat->preallocated = PETSC_TRUE;
195: if (dim == 3) {
196: mat->ops->setvalueslocal = MatSetValuesLocal_HYPREStruct_3d;
197: mat->ops->zerorowslocal = MatZeroRowsLocal_HYPREStruct_3d;
198: mat->ops->zeroentries = MatZeroEntries_HYPREStruct_3d;
200: /* MatZeroEntries_HYPREStruct_3d(mat); */
201: } else SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Only support for 3d DMDA currently");
203: /* get values that will be used repeatedly in MatSetValuesLocal() and MatZeroRowsLocal() repeatedly */
204: MatGetOwnershipRange(mat,&ex->rstart,NULL);
205: DMGetLocalToGlobalMapping(ex->da,<og);
206: ISLocalToGlobalMappingGetIndices(ltog, (const PetscInt **) &ex->gindices);
207: DMDAGetGhostCorners(ex->da,0,0,0,&ex->gnx,&ex->gnxgny,0);
208: ex->gnxgny *= ex->gnx;
209: DMDAGetCorners(ex->da,&ex->xs,&ex->ys,&ex->zs,&ex->nx,&ex->ny,0);
210: ex->nxny = ex->nx*ex->ny;
211: return(0);
212: }
214: PetscErrorCode MatMult_HYPREStruct(Mat A,Vec x,Vec y)
215: {
216: PetscErrorCode ierr;
217: const PetscScalar *xx;
218: PetscScalar *yy;
219: PetscInt ilower[3],iupper[3];
220: Mat_HYPREStruct *mx = (Mat_HYPREStruct*)(A->data);
223: DMDAGetCorners(mx->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);
225: iupper[0] += ilower[0] - 1;
226: iupper[1] += ilower[1] - 1;
227: iupper[2] += ilower[2] - 1;
229: /* copy x values over to hypre */
230: PetscStackCallStandard(HYPRE_StructVectorSetConstantValues,(mx->hb,0.0));
231: VecGetArrayRead(x,&xx);
232: PetscStackCallStandard(HYPRE_StructVectorSetBoxValues,(mx->hb,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,(PetscScalar*)xx));
233: VecRestoreArrayRead(x,&xx);
234: PetscStackCallStandard(HYPRE_StructVectorAssemble,(mx->hb));
235: PetscStackCallStandard(HYPRE_StructMatrixMatvec,(1.0,mx->hmat,mx->hb,0.0,mx->hx));
237: /* copy solution values back to PETSc */
238: VecGetArray(y,&yy);
239: PetscStackCallStandard(HYPRE_StructVectorGetBoxValues,(mx->hx,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,yy));
240: VecRestoreArray(y,&yy);
241: return(0);
242: }
244: PetscErrorCode MatAssemblyEnd_HYPREStruct(Mat mat,MatAssemblyType mode)
245: {
246: Mat_HYPREStruct *ex = (Mat_HYPREStruct*) mat->data;
247: PetscErrorCode ierr;
250: PetscStackCallStandard(HYPRE_StructMatrixAssemble,(ex->hmat));
251: /* PetscStackCallStandard(HYPRE_StructMatrixPrint,("dummy",ex->hmat,0)); */
252: return(0);
253: }
255: PetscErrorCode MatZeroEntries_HYPREStruct(Mat mat)
256: {
258: /* before the DMDA is set to the matrix the zero doesn't need to do anything */
259: return(0);
260: }
262: PetscErrorCode MatDestroy_HYPREStruct(Mat mat)
263: {
264: Mat_HYPREStruct *ex = (Mat_HYPREStruct*) mat->data;
265: PetscErrorCode ierr;
268: PetscStackCallStandard(HYPRE_StructMatrixDestroy,(ex->hmat));
269: PetscStackCallStandard(HYPRE_StructVectorDestroy,(ex->hx));
270: PetscStackCallStandard(HYPRE_StructVectorDestroy,(ex->hb));
271: PetscObjectDereference((PetscObject)ex->da);
272: MPI_Comm_free(&(ex->hcomm));
273: PetscFree(ex);
274: return(0);
275: }
277: PETSC_EXTERN PetscErrorCode MatCreate_HYPREStruct(Mat B)
278: {
279: Mat_HYPREStruct *ex;
280: PetscErrorCode ierr;
283: PetscNewLog(B,&ex);
284: B->data = (void*)ex;
285: B->rmap->bs = 1;
286: B->assembled = PETSC_FALSE;
288: B->insertmode = NOT_SET_VALUES;
290: B->ops->assemblyend = MatAssemblyEnd_HYPREStruct;
291: B->ops->mult = MatMult_HYPREStruct;
292: B->ops->zeroentries = MatZeroEntries_HYPREStruct;
293: B->ops->destroy = MatDestroy_HYPREStruct;
294: B->ops->setup = MatSetUp_HYPREStruct;
296: ex->needsinitialization = PETSC_TRUE;
298: MPI_Comm_dup(PetscObjectComm((PetscObject)B),&(ex->hcomm));
299: PetscObjectChangeTypeName((PetscObject)B,MATHYPRESTRUCT);
300: return(0);
301: }
303: /*MC
304: MATHYPRESSTRUCT - MATHYPRESSTRUCT = "hypresstruct" - A matrix type to be used for parallel sparse matrices
305: based on the hypre HYPRE_SStructMatrix.
308: Level: intermediate
310: Notes: Unlike hypre's general semi-struct object consisting of a collection of structured-grid objects and unstructured
311: grid objects, we restrict the semi-struct objects to consist of only structured-grid components.
313: Unlike the more general support for parts and blocks in hypre this allows only one part, and one block per process and requires the block
314: be defined by a DMDA.
316: The matrix needs a DMDA associated with it by either a call to MatSetDM() or if the matrix is obtained from DMCreateMatrix()
318: M*/
320: PetscErrorCode MatSetValuesLocal_HYPRESStruct_3d(Mat mat,PetscInt nrow,const PetscInt irow[],PetscInt ncol,const PetscInt icol[],const PetscScalar y[],InsertMode addv)
321: {
322: PetscErrorCode ierr;
323: PetscInt i,j,stencil,index[3];
324: const PetscScalar *values = y;
325: Mat_HYPRESStruct *ex = (Mat_HYPRESStruct*) mat->data;
327: PetscInt part = 0; /* Petsc sstruct interface only allows 1 part */
328: PetscInt ordering;
329: PetscInt grid_rank, to_grid_rank;
330: PetscInt var_type, to_var_type;
331: PetscInt to_var_entry = 0;
332: PetscInt nvars= ex->nvars;
333: PetscInt row,*entries;
336: PetscMalloc1(7*nvars,&entries);
337: ordering = ex->dofs_order; /* ordering= 0 nodal ordering
338: 1 variable ordering */
339: /* stencil entries are ordered by variables: var0_stencil0, var0_stencil1, ..., var0_stencil6, var1_stencil0, var1_stencil1, ... */
340: if (!ordering) {
341: for (i=0; i<nrow; i++) {
342: grid_rank= irow[i]/nvars;
343: var_type = (irow[i] % nvars);
345: for (j=0; j<ncol; j++) {
346: to_grid_rank= icol[j]/nvars;
347: to_var_type = (icol[j] % nvars);
349: to_var_entry= to_var_entry*7;
350: entries[j] = to_var_entry;
352: stencil = to_grid_rank-grid_rank;
353: if (!stencil) {
354: entries[j] += 3;
355: } else if (stencil == -1) {
356: entries[j] += 2;
357: } else if (stencil == 1) {
358: entries[j] += 4;
359: } else if (stencil == -ex->gnx) {
360: entries[j] += 1;
361: } else if (stencil == ex->gnx) {
362: entries[j] += 5;
363: } else if (stencil == -ex->gnxgny) {
364: entries[j] += 0;
365: } else if (stencil == ex->gnxgny) {
366: entries[j] += 6;
367: } else SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Local row %D local column %D have bad stencil %D",irow[i],icol[j],stencil);
368: }
370: row = ex->gindices[grid_rank] - ex->rstart;
371: index[0] = ex->xs + (row % ex->nx);
372: index[1] = ex->ys + ((row/ex->nx) % ex->ny);
373: index[2] = ex->zs + (row/(ex->nxny));
375: if (addv == ADD_VALUES) {
376: PetscStackCallStandard(HYPRE_SStructMatrixAddToValues,(ex->ss_mat,part,(HYPRE_Int *)index,var_type,ncol,(HYPRE_Int *)entries,(PetscScalar*)values));
377: } else {
378: PetscStackCallStandard(HYPRE_SStructMatrixSetValues,(ex->ss_mat,part,(HYPRE_Int *)index,var_type,ncol,(HYPRE_Int *)entries,(PetscScalar*)values));
379: }
380: values += ncol;
381: }
382: } else {
383: for (i=0; i<nrow; i++) {
384: var_type = irow[i]/(ex->gnxgnygnz);
385: grid_rank= irow[i] - var_type*(ex->gnxgnygnz);
387: for (j=0; j<ncol; j++) {
388: to_var_type = icol[j]/(ex->gnxgnygnz);
389: to_grid_rank= icol[j] - to_var_type*(ex->gnxgnygnz);
391: to_var_entry= to_var_entry*7;
392: entries[j] = to_var_entry;
394: stencil = to_grid_rank-grid_rank;
395: if (!stencil) {
396: entries[j] += 3;
397: } else if (stencil == -1) {
398: entries[j] += 2;
399: } else if (stencil == 1) {
400: entries[j] += 4;
401: } else if (stencil == -ex->gnx) {
402: entries[j] += 1;
403: } else if (stencil == ex->gnx) {
404: entries[j] += 5;
405: } else if (stencil == -ex->gnxgny) {
406: entries[j] += 0;
407: } else if (stencil == ex->gnxgny) {
408: entries[j] += 6;
409: } else SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Local row %D local column %D have bad stencil %D",irow[i],icol[j],stencil);
410: }
412: row = ex->gindices[grid_rank] - ex->rstart;
413: index[0] = ex->xs + (row % ex->nx);
414: index[1] = ex->ys + ((row/ex->nx) % ex->ny);
415: index[2] = ex->zs + (row/(ex->nxny));
417: if (addv == ADD_VALUES) {
418: PetscStackCallStandard(HYPRE_SStructMatrixAddToValues,(ex->ss_mat,part,(HYPRE_Int *)index,var_type,ncol,(HYPRE_Int *)entries,(PetscScalar*)values));
419: } else {
420: PetscStackCallStandard(HYPRE_SStructMatrixSetValues,(ex->ss_mat,part,(HYPRE_Int *)index,var_type,ncol,(HYPRE_Int *)entries,(PetscScalar*)values));
421: }
422: values += ncol;
423: }
424: }
425: PetscFree(entries);
426: return(0);
427: }
429: PetscErrorCode MatZeroRowsLocal_HYPRESStruct_3d(Mat mat,PetscInt nrow,const PetscInt irow[],PetscScalar d,Vec x,Vec b)
430: {
431: PetscErrorCode ierr;
432: PetscInt i,index[3];
433: PetscScalar **values;
434: Mat_HYPRESStruct *ex = (Mat_HYPRESStruct*) mat->data;
436: PetscInt part = 0; /* Petsc sstruct interface only allows 1 part */
437: PetscInt ordering = ex->dofs_order;
438: PetscInt grid_rank;
439: PetscInt var_type;
440: PetscInt nvars= ex->nvars;
441: PetscInt row,*entries;
444: if (x && b) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"No support");
445: PetscMalloc1(7*nvars,&entries);
447: PetscMalloc1(nvars,&values);
448: PetscMalloc1(7*nvars*nvars,&values[0]);
449: for (i=1; i<nvars; i++) {
450: values[i] = values[i-1] + nvars*7;
451: }
453: for (i=0; i< nvars; i++) {
454: PetscMemzero(values[i],nvars*7*sizeof(PetscScalar));
455: *(values[i]+3) = d;
456: }
458: for (i= 0; i< nvars*7; i++) entries[i] = i;
460: if (!ordering) {
461: for (i=0; i<nrow; i++) {
462: grid_rank = irow[i] / nvars;
463: var_type =(irow[i] % nvars);
465: row = ex->gindices[grid_rank] - ex->rstart;
466: index[0] = ex->xs + (row % ex->nx);
467: index[1] = ex->ys + ((row/ex->nx) % ex->ny);
468: index[2] = ex->zs + (row/(ex->nxny));
469: PetscStackCallStandard(HYPRE_SStructMatrixSetValues,(ex->ss_mat,part,(HYPRE_Int *)index,var_type,7*nvars,(HYPRE_Int *)entries,values[var_type]));
470: }
471: } else {
472: for (i=0; i<nrow; i++) {
473: var_type = irow[i]/(ex->gnxgnygnz);
474: grid_rank= irow[i] - var_type*(ex->gnxgnygnz);
476: row = ex->gindices[grid_rank] - ex->rstart;
477: index[0] = ex->xs + (row % ex->nx);
478: index[1] = ex->ys + ((row/ex->nx) % ex->ny);
479: index[2] = ex->zs + (row/(ex->nxny));
480: PetscStackCallStandard(HYPRE_SStructMatrixSetValues,(ex->ss_mat,part,(HYPRE_Int *)index,var_type,7*nvars,(HYPRE_Int *)entries,values[var_type]));
481: }
482: }
483: PetscStackCallStandard(HYPRE_SStructMatrixAssemble,(ex->ss_mat));
484: PetscFree(values[0]);
485: PetscFree(values);
486: PetscFree(entries);
487: return(0);
488: }
490: PetscErrorCode MatZeroEntries_HYPRESStruct_3d(Mat mat)
491: {
492: PetscErrorCode ierr;
493: Mat_HYPRESStruct *ex = (Mat_HYPRESStruct*) mat->data;
494: PetscInt nvars= ex->nvars;
495: PetscInt size;
496: PetscInt part= 0; /* only one part */
499: size = ((ex->hbox.imax[0])-(ex->hbox.imin[0])+1)*((ex->hbox.imax[1])-(ex->hbox.imin[1])+1)*((ex->hbox.imax[2])-(ex->hbox.imin[2])+1);
500: {
501: PetscInt i,*entries,iupper[3],ilower[3];
502: PetscScalar *values;
505: for (i= 0; i< 3; i++) {
506: ilower[i]= ex->hbox.imin[i];
507: iupper[i]= ex->hbox.imax[i];
508: }
510: PetscMalloc2(nvars*7,&entries,nvars*7*size,&values);
511: for (i= 0; i< nvars*7; i++) entries[i]= i;
512: PetscMemzero(values,nvars*7*size*sizeof(PetscScalar));
514: for (i= 0; i< nvars; i++) {
515: PetscStackCallStandard(HYPRE_SStructMatrixSetBoxValues,(ex->ss_mat,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,nvars*7,(HYPRE_Int *)entries,values));
516: }
517: PetscFree2(entries,values);
518: }
519: PetscStackCallStandard(HYPRE_SStructMatrixAssemble,(ex->ss_mat));
520: return(0);
521: }
524: static PetscErrorCode MatSetUp_HYPRESStruct(Mat mat)
525: {
526: PetscErrorCode ierr;
527: Mat_HYPRESStruct *ex = (Mat_HYPRESStruct*) mat->data;
528: PetscInt dim,dof,sw[3],nx,ny,nz;
529: PetscInt ilower[3],iupper[3],ssize,i;
530: DMBoundaryType px,py,pz;
531: DMDAStencilType st;
532: PetscInt nparts= 1; /* assuming only one part */
533: PetscInt part = 0;
534: ISLocalToGlobalMapping ltog;
535: DM da;
538: MatGetDM(mat,(DM*)&da);
539: ex->da = da;
540: PetscObjectReference((PetscObject)da);
542: DMDAGetInfo(ex->da,&dim,0,0,0,0,0,0,&dof,&sw[0],&px,&py,&pz,&st);
543: DMDAGetCorners(ex->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);
544: iupper[0] += ilower[0] - 1;
545: iupper[1] += ilower[1] - 1;
546: iupper[2] += ilower[2] - 1;
547: /* the hypre_Box is used to zero out the matrix entries in MatZeroValues() */
548: ex->hbox.imin[0] = ilower[0];
549: ex->hbox.imin[1] = ilower[1];
550: ex->hbox.imin[2] = ilower[2];
551: ex->hbox.imax[0] = iupper[0];
552: ex->hbox.imax[1] = iupper[1];
553: ex->hbox.imax[2] = iupper[2];
555: ex->dofs_order = 0;
557: /* assuming that the same number of dofs on each gridpoint. Also assume all cell-centred based */
558: ex->nvars= dof;
560: /* create the hypre grid object and set its information */
561: if (px || py || pz) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Ask us to add periodic support by calling HYPRE_SStructGridSetPeriodic()");
562: PetscStackCallStandard(HYPRE_SStructGridCreate,(ex->hcomm,dim,nparts,&ex->ss_grid));
563: PetscStackCallStandard(HYPRE_SStructGridSetExtents,(ex->ss_grid,part,ex->hbox.imin,ex->hbox.imax));
564: {
565: HYPRE_SStructVariable *vartypes;
566: PetscMalloc1(ex->nvars,&vartypes);
567: for (i= 0; i< ex->nvars; i++) vartypes[i]= HYPRE_SSTRUCT_VARIABLE_CELL;
568: PetscStackCallStandard(HYPRE_SStructGridSetVariables,(ex->ss_grid, part, ex->nvars,vartypes));
569: PetscFree(vartypes);
570: }
571: PetscStackCallStandard(HYPRE_SStructGridAssemble,(ex->ss_grid));
573: sw[1] = sw[0];
574: sw[2] = sw[1];
575: /* PetscStackCallStandard(HYPRE_SStructGridSetNumGhost,(ex->ss_grid,sw)); */
577: /* create the hypre stencil object and set its information */
578: if (sw[0] > 1) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Ask us to add support for wider stencils");
579: if (st == DMDA_STENCIL_BOX) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Ask us to add support for box stencils");
581: if (dim == 1) {
582: PetscInt offsets[3][1] = {{-1},{0},{1}};
583: PetscInt j, cnt;
585: ssize = 3*(ex->nvars);
586: PetscStackCallStandard(HYPRE_SStructStencilCreate,(dim,ssize,&ex->ss_stencil));
587: cnt= 0;
588: for (i = 0; i < (ex->nvars); i++) {
589: for (j = 0; j < 3; j++) {
590: PetscStackCallStandard(HYPRE_SStructStencilSetEntry,(ex->ss_stencil, cnt, (HYPRE_Int *)offsets[j], i));
591: cnt++;
592: }
593: }
594: } else if (dim == 2) {
595: PetscInt offsets[5][2] = {{0,-1},{-1,0},{0,0},{1,0},{0,1}};
596: PetscInt j, cnt;
598: ssize = 5*(ex->nvars);
599: PetscStackCallStandard(HYPRE_SStructStencilCreate,(dim,ssize,&ex->ss_stencil));
600: cnt= 0;
601: for (i= 0; i< (ex->nvars); i++) {
602: for (j= 0; j< 5; j++) {
603: PetscStackCallStandard(HYPRE_SStructStencilSetEntry,(ex->ss_stencil, cnt, (HYPRE_Int *)offsets[j], i));
604: cnt++;
605: }
606: }
607: } else if (dim == 3) {
608: PetscInt offsets[7][3] = {{0,0,-1},{0,-1,0},{-1,0,0},{0,0,0},{1,0,0},{0,1,0},{0,0,1}};
609: PetscInt j, cnt;
611: ssize = 7*(ex->nvars);
612: PetscStackCallStandard(HYPRE_SStructStencilCreate,(dim,ssize,&ex->ss_stencil));
613: cnt= 0;
614: for (i= 0; i< (ex->nvars); i++) {
615: for (j= 0; j< 7; j++) {
616: PetscStackCallStandard(HYPRE_SStructStencilSetEntry,(ex->ss_stencil, cnt, (HYPRE_Int *)offsets[j], i));
617: cnt++;
618: }
619: }
620: }
622: /* create the HYPRE graph */
623: PetscStackCallStandard(HYPRE_SStructGraphCreate,(ex->hcomm, ex->ss_grid, &(ex->ss_graph)));
625: /* set the stencil graph. Note that each variable has the same graph. This means that each
626: variable couples to all the other variable and with the same stencil pattern. */
627: for (i= 0; i< (ex->nvars); i++) {
628: PetscStackCallStandard(HYPRE_SStructGraphSetStencil,(ex->ss_graph,part,i,ex->ss_stencil));
629: }
630: PetscStackCallStandard(HYPRE_SStructGraphAssemble,(ex->ss_graph));
632: /* create the HYPRE sstruct vectors for rhs and solution */
633: PetscStackCallStandard(HYPRE_SStructVectorCreate,(ex->hcomm,ex->ss_grid,&ex->ss_b));
634: PetscStackCallStandard(HYPRE_SStructVectorCreate,(ex->hcomm,ex->ss_grid,&ex->ss_x));
635: PetscStackCallStandard(HYPRE_SStructVectorInitialize,(ex->ss_b));
636: PetscStackCallStandard(HYPRE_SStructVectorInitialize,(ex->ss_x));
637: PetscStackCallStandard(HYPRE_SStructVectorAssemble,(ex->ss_b));
638: PetscStackCallStandard(HYPRE_SStructVectorAssemble,(ex->ss_x));
640: /* create the hypre matrix object and set its information */
641: PetscStackCallStandard(HYPRE_SStructMatrixCreate,(ex->hcomm,ex->ss_graph,&ex->ss_mat));
642: PetscStackCallStandard(HYPRE_SStructGridDestroy,(ex->ss_grid));
643: PetscStackCallStandard(HYPRE_SStructStencilDestroy,(ex->ss_stencil));
644: if (ex->needsinitialization) {
645: PetscStackCallStandard(HYPRE_SStructMatrixInitialize,(ex->ss_mat));
646: ex->needsinitialization = PETSC_FALSE;
647: }
649: /* set the global and local sizes of the matrix */
650: DMDAGetCorners(da,0,0,0,&nx,&ny,&nz);
651: MatSetSizes(mat,dof*nx*ny*nz,dof*nx*ny*nz,PETSC_DECIDE,PETSC_DECIDE);
652: PetscLayoutSetBlockSize(mat->rmap,dof);
653: PetscLayoutSetBlockSize(mat->cmap,dof);
654: PetscLayoutSetUp(mat->rmap);
655: PetscLayoutSetUp(mat->cmap);
656: mat->preallocated = PETSC_TRUE;
658: if (dim == 3) {
659: mat->ops->setvalueslocal = MatSetValuesLocal_HYPRESStruct_3d;
660: mat->ops->zerorowslocal = MatZeroRowsLocal_HYPRESStruct_3d;
661: mat->ops->zeroentries = MatZeroEntries_HYPRESStruct_3d;
663: /* MatZeroEntries_HYPRESStruct_3d(mat); */
664: } else SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Only support for 3d DMDA currently");
666: /* get values that will be used repeatedly in MatSetValuesLocal() and MatZeroRowsLocal() repeatedly */
667: MatGetOwnershipRange(mat,&ex->rstart,NULL);
668: DMGetLocalToGlobalMapping(ex->da,<og);
669: ISLocalToGlobalMappingGetIndices(ltog, (const PetscInt **) &ex->gindices);
670: DMDAGetGhostCorners(ex->da,0,0,0,&ex->gnx,&ex->gnxgny,&ex->gnxgnygnz);
672: ex->gnxgny *= ex->gnx;
673: ex->gnxgnygnz *= ex->gnxgny;
675: DMDAGetCorners(ex->da,&ex->xs,&ex->ys,&ex->zs,&ex->nx,&ex->ny,&ex->nz);
677: ex->nxny = ex->nx*ex->ny;
678: ex->nxnynz = ex->nz*ex->nxny;
679: return(0);
680: }
682: PetscErrorCode MatMult_HYPRESStruct(Mat A,Vec x,Vec y)
683: {
684: PetscErrorCode ierr;
685: const PetscScalar *xx;
686: PetscScalar *yy;
687: PetscInt ilower[3],iupper[3];
688: Mat_HYPRESStruct *mx = (Mat_HYPRESStruct*)(A->data);
689: PetscInt ordering= mx->dofs_order;
690: PetscInt nvars = mx->nvars;
691: PetscInt part = 0;
692: PetscInt size;
693: PetscInt i;
696: DMDAGetCorners(mx->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);
697: iupper[0] += ilower[0] - 1;
698: iupper[1] += ilower[1] - 1;
699: iupper[2] += ilower[2] - 1;
701: size= 1;
702: for (i = 0; i < 3; i++) size *= (iupper[i]-ilower[i]+1);
704: /* copy x values over to hypre for variable ordering */
705: if (ordering) {
706: PetscStackCallStandard(HYPRE_SStructVectorSetConstantValues,(mx->ss_b,0.0));
707: VecGetArrayRead(x,&xx);
708: for (i= 0; i< nvars; i++) {
709: PetscStackCallStandard(HYPRE_SStructVectorSetBoxValues,(mx->ss_b,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,(PetscScalar*)xx+(size*i)));
710: }
711: VecRestoreArrayRead(x,&xx);
712: PetscStackCallStandard(HYPRE_SStructVectorAssemble,(mx->ss_b));
713: PetscStackCallStandard(HYPRE_SStructMatrixMatvec,(1.0,mx->ss_mat,mx->ss_b,0.0,mx->ss_x));
715: /* copy solution values back to PETSc */
716: VecGetArray(y,&yy);
717: for (i= 0; i< nvars; i++) {
718: PetscStackCallStandard(HYPRE_SStructVectorGetBoxValues,(mx->ss_x,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,yy+(size*i)));
719: }
720: VecRestoreArray(y,&yy);
721: } else { /* nodal ordering must be mapped to variable ordering for sys_pfmg */
722: PetscScalar *z;
723: PetscInt j, k;
725: PetscMalloc1(nvars*size,&z);
726: PetscStackCallStandard(HYPRE_SStructVectorSetConstantValues,(mx->ss_b,0.0));
727: VecGetArrayRead(x,&xx);
729: /* transform nodal to hypre's variable ordering for sys_pfmg */
730: for (i= 0; i< size; i++) {
731: k= i*nvars;
732: for (j= 0; j< nvars; j++) z[j*size+i]= xx[k+j];
733: }
734: for (i= 0; i< nvars; i++) {
735: PetscStackCallStandard(HYPRE_SStructVectorSetBoxValues,(mx->ss_b,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,z+(size*i)));
736: }
737: VecRestoreArrayRead(x,&xx);
738: PetscStackCallStandard(HYPRE_SStructVectorAssemble,(mx->ss_b));
739: PetscStackCallStandard(HYPRE_SStructMatrixMatvec,(1.0,mx->ss_mat,mx->ss_b,0.0,mx->ss_x));
741: /* copy solution values back to PETSc */
742: VecGetArray(y,&yy);
743: for (i= 0; i< nvars; i++) {
744: PetscStackCallStandard(HYPRE_SStructVectorGetBoxValues,(mx->ss_x,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,z+(size*i)));
745: }
746: /* transform hypre's variable ordering for sys_pfmg to nodal ordering */
747: for (i= 0; i< size; i++) {
748: k= i*nvars;
749: for (j= 0; j< nvars; j++) yy[k+j]= z[j*size+i];
750: }
751: VecRestoreArray(y,&yy);
752: PetscFree(z);
753: }
754: return(0);
755: }
757: PetscErrorCode MatAssemblyEnd_HYPRESStruct(Mat mat,MatAssemblyType mode)
758: {
759: Mat_HYPRESStruct *ex = (Mat_HYPRESStruct*) mat->data;
760: PetscErrorCode ierr;
763: PetscStackCallStandard(HYPRE_SStructMatrixAssemble,(ex->ss_mat));
764: return(0);
765: }
767: PetscErrorCode MatZeroEntries_HYPRESStruct(Mat mat)
768: {
770: /* before the DMDA is set to the matrix the zero doesn't need to do anything */
771: return(0);
772: }
774: PetscErrorCode MatDestroy_HYPRESStruct(Mat mat)
775: {
776: Mat_HYPRESStruct *ex = (Mat_HYPRESStruct*) mat->data;
777: PetscErrorCode ierr;
778: ISLocalToGlobalMapping ltog;
781: DMGetLocalToGlobalMapping(ex->da,<og);
782: ISLocalToGlobalMappingRestoreIndices(ltog, (const PetscInt **) &ex->gindices);
783: PetscStackCallStandard(HYPRE_SStructGraphDestroy,(ex->ss_graph));
784: PetscStackCallStandard(HYPRE_SStructMatrixDestroy,(ex->ss_mat));
785: PetscStackCallStandard(HYPRE_SStructVectorDestroy,(ex->ss_x));
786: PetscStackCallStandard(HYPRE_SStructVectorDestroy,(ex->ss_b));
787: PetscObjectDereference((PetscObject)ex->da);
788: MPI_Comm_free(&(ex->hcomm));
789: PetscFree(ex);
790: return(0);
791: }
793: PETSC_EXTERN PetscErrorCode MatCreate_HYPRESStruct(Mat B)
794: {
795: Mat_HYPRESStruct *ex;
796: PetscErrorCode ierr;
799: PetscNewLog(B,&ex);
800: B->data = (void*)ex;
801: B->rmap->bs = 1;
802: B->assembled = PETSC_FALSE;
804: B->insertmode = NOT_SET_VALUES;
806: B->ops->assemblyend = MatAssemblyEnd_HYPRESStruct;
807: B->ops->mult = MatMult_HYPRESStruct;
808: B->ops->zeroentries = MatZeroEntries_HYPRESStruct;
809: B->ops->destroy = MatDestroy_HYPRESStruct;
810: B->ops->setup = MatSetUp_HYPRESStruct;
812: ex->needsinitialization = PETSC_TRUE;
814: MPI_Comm_dup(PetscObjectComm((PetscObject)B),&(ex->hcomm));
815: PetscObjectChangeTypeName((PetscObject)B,MATHYPRESSTRUCT);
816: return(0);
817: }