Actual source code: mhyp.c
petsc-3.10.5 2019-03-28
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:
19: Unlike the more general support for blocks in hypre this allows only one block per process and requires the block
20: be defined by a DMDA.
22: The matrix needs a DMDA associated with it by either a call to MatSetDM() or if the matrix is obtained from DMCreateMatrix()
24: .seealso: MatCreate(), PCPFMG, MatSetDM(), DMCreateMatrix()
25: 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: PetscErrorCode ierr;
31: PetscInt i,j,stencil,index[3],row,entries[7];
32: const PetscScalar *values = 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] = ex->xs + (row % ex->nx);
58: index[1] = ex->ys + ((row/ex->nx) % ex->ny);
59: index[2] = ex->zs + (row/(ex->nxny));
60: if (addv == ADD_VALUES) {
61: PetscStackCallStandard(HYPRE_StructMatrixAddToValues,(ex->hmat,(HYPRE_Int *)index,ncol,(HYPRE_Int *)entries,(PetscScalar*)values));
62: } else {
63: PetscStackCallStandard(HYPRE_StructMatrixSetValues,(ex->hmat,(HYPRE_Int *)index,ncol,(HYPRE_Int *)entries,(PetscScalar*)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: PetscInt i,index[3],row,entries[7] = {0,1,2,3,4,5,6};
74: PetscScalar values[7];
75: Mat_HYPREStruct *ex = (Mat_HYPREStruct*) mat->data;
78: if (x && b) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"No support");
79: PetscMemzero(values,7*sizeof(PetscScalar));
80: values[3] = d;
81: for (i=0; i<nrow; i++) {
82: row = ex->gindices[irow[i]] - ex->rstart;
83: index[0] = ex->xs + (row % ex->nx);
84: index[1] = ex->ys + ((row/ex->nx) % ex->ny);
85: index[2] = ex->zs + (row/(ex->nxny));
86: PetscStackCallStandard(HYPRE_StructMatrixSetValues,(ex->hmat,(HYPRE_Int *)index,7,(HYPRE_Int *)entries,values));
87: }
88: PetscStackCallStandard(HYPRE_StructMatrixAssemble,(ex->hmat));
89: return(0);
90: }
92: PetscErrorCode MatZeroEntries_HYPREStruct_3d(Mat mat)
93: {
94: PetscErrorCode ierr;
95: PetscInt 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,(HYPRE_Int *)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: PetscInt dim,dof,sw[6],nx,ny,nz,ilower[3],iupper[3],ssize,i;
110: DMBoundaryType px,py,pz;
111: DMDAStencilType st;
112: ISLocalToGlobalMapping ltog;
113: DM da;
116: MatGetDM(mat,(DM*)&da);
117: ex->da = da;
118: PetscObjectReference((PetscObject)da);
120: DMDAGetInfo(ex->da,&dim,0,0,0,0,0,0,&dof,&sw[0],&px,&py,&pz,&st);
121: DMDAGetCorners(ex->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);
122: iupper[0] += ilower[0] - 1;
123: iupper[1] += ilower[1] - 1;
124: iupper[2] += ilower[2] - 1;
126: /* the hypre_Box is used to zero out the matrix entries in MatZeroValues() */
127: ex->hbox.imin[0] = ilower[0];
128: ex->hbox.imin[1] = ilower[1];
129: ex->hbox.imin[2] = ilower[2];
130: ex->hbox.imax[0] = iupper[0];
131: ex->hbox.imax[1] = iupper[1];
132: ex->hbox.imax[2] = iupper[2];
134: /* create the hypre grid object and set its information */
135: if (dof > 1) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Currently only support for scalar problems");
136: if (px || py || pz) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Ask us to add periodic support by calling HYPRE_StructGridSetPeriodic()");
137: PetscStackCallStandard(HYPRE_StructGridCreate,(ex->hcomm,dim,&ex->hgrid));
138: PetscStackCallStandard(HYPRE_StructGridSetExtents,(ex->hgrid,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper));
139: PetscStackCallStandard(HYPRE_StructGridAssemble,(ex->hgrid));
141: sw[5] = sw[4] = sw[3] = sw[2] = sw[1] = sw[0];
142: PetscStackCallStandard(HYPRE_StructGridSetNumGhost,(ex->hgrid,(HYPRE_Int *)sw));
144: /* create the hypre stencil object and set its information */
145: if (sw[0] > 1) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Ask us to add support for wider stencils");
146: if (st == DMDA_STENCIL_BOX) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Ask us to add support for box stencils");
147: if (dim == 1) {
148: PetscInt offsets[3][1] = {{-1},{0},{1}};
149: ssize = 3;
150: PetscStackCallStandard(HYPRE_StructStencilCreate,(dim,ssize,&ex->hstencil));
151: for (i=0; i<ssize; i++) {
152: PetscStackCallStandard(HYPRE_StructStencilSetElement,(ex->hstencil,i,(HYPRE_Int *)offsets[i]));
153: }
154: } else if (dim == 2) {
155: PetscInt offsets[5][2] = {{0,-1},{-1,0},{0,0},{1,0},{0,1}};
156: ssize = 5;
157: PetscStackCallStandard(HYPRE_StructStencilCreate,(dim,ssize,&ex->hstencil));
158: for (i=0; i<ssize; i++) {
159: PetscStackCallStandard(HYPRE_StructStencilSetElement,(ex->hstencil,i,(HYPRE_Int *)offsets[i]));
160: }
161: } else if (dim == 3) {
162: 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}};
163: ssize = 7;
164: PetscStackCallStandard(HYPRE_StructStencilCreate,(dim,ssize,&ex->hstencil));
165: for (i=0; i<ssize; i++) {
166: PetscStackCallStandard(HYPRE_StructStencilSetElement,(ex->hstencil,i,(HYPRE_Int *)offsets[i]));
167: }
168: }
170: /* create the HYPRE vector for rhs and solution */
171: PetscStackCallStandard(HYPRE_StructVectorCreate,(ex->hcomm,ex->hgrid,&ex->hb));
172: PetscStackCallStandard(HYPRE_StructVectorCreate,(ex->hcomm,ex->hgrid,&ex->hx));
173: PetscStackCallStandard(HYPRE_StructVectorInitialize,(ex->hb));
174: PetscStackCallStandard(HYPRE_StructVectorInitialize,(ex->hx));
175: PetscStackCallStandard(HYPRE_StructVectorAssemble,(ex->hb));
176: PetscStackCallStandard(HYPRE_StructVectorAssemble,(ex->hx));
178: /* create the hypre matrix object and set its information */
179: PetscStackCallStandard(HYPRE_StructMatrixCreate,(ex->hcomm,ex->hgrid,ex->hstencil,&ex->hmat));
180: PetscStackCallStandard(HYPRE_StructGridDestroy,(ex->hgrid));
181: PetscStackCallStandard(HYPRE_StructStencilDestroy,(ex->hstencil));
182: if (ex->needsinitialization) {
183: PetscStackCallStandard(HYPRE_StructMatrixInitialize,(ex->hmat));
184: ex->needsinitialization = PETSC_FALSE;
185: }
187: /* set the global and local sizes of the matrix */
188: DMDAGetCorners(da,0,0,0,&nx,&ny,&nz);
189: MatSetSizes(mat,dof*nx*ny*nz,dof*nx*ny*nz,PETSC_DECIDE,PETSC_DECIDE);
190: PetscLayoutSetBlockSize(mat->rmap,1);
191: PetscLayoutSetBlockSize(mat->cmap,1);
192: PetscLayoutSetUp(mat->rmap);
193: PetscLayoutSetUp(mat->cmap);
194: mat->preallocated = PETSC_TRUE;
196: if (dim == 3) {
197: mat->ops->setvalueslocal = MatSetValuesLocal_HYPREStruct_3d;
198: mat->ops->zerorowslocal = MatZeroRowsLocal_HYPREStruct_3d;
199: mat->ops->zeroentries = MatZeroEntries_HYPREStruct_3d;
201: /* MatZeroEntries_HYPREStruct_3d(mat); */
202: } else SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Only support for 3d DMDA currently");
204: /* get values that will be used repeatedly in MatSetValuesLocal() and MatZeroRowsLocal() repeatedly */
205: MatGetOwnershipRange(mat,&ex->rstart,NULL);
206: DMGetLocalToGlobalMapping(ex->da,<og);
207: ISLocalToGlobalMappingGetIndices(ltog, (const PetscInt **) &ex->gindices);
208: DMDAGetGhostCorners(ex->da,0,0,0,&ex->gnx,&ex->gnxgny,0);
209: ex->gnxgny *= ex->gnx;
210: DMDAGetCorners(ex->da,&ex->xs,&ex->ys,&ex->zs,&ex->nx,&ex->ny,0);
211: ex->nxny = ex->nx*ex->ny;
212: return(0);
213: }
215: PetscErrorCode MatMult_HYPREStruct(Mat A,Vec x,Vec y)
216: {
217: PetscErrorCode ierr;
218: const PetscScalar *xx;
219: PetscScalar *yy;
220: PetscInt ilower[3],iupper[3];
221: Mat_HYPREStruct *mx = (Mat_HYPREStruct*)(A->data);
224: DMDAGetCorners(mx->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);
226: iupper[0] += ilower[0] - 1;
227: iupper[1] += ilower[1] - 1;
228: iupper[2] += ilower[2] - 1;
230: /* copy x values over to hypre */
231: PetscStackCallStandard(HYPRE_StructVectorSetConstantValues,(mx->hb,0.0));
232: VecGetArrayRead(x,&xx);
233: PetscStackCallStandard(HYPRE_StructVectorSetBoxValues,(mx->hb,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,(PetscScalar*)xx));
234: VecRestoreArrayRead(x,&xx);
235: PetscStackCallStandard(HYPRE_StructVectorAssemble,(mx->hb));
236: PetscStackCallStandard(HYPRE_StructMatrixMatvec,(1.0,mx->hmat,mx->hb,0.0,mx->hx));
238: /* copy solution values back to PETSc */
239: VecGetArray(y,&yy);
240: PetscStackCallStandard(HYPRE_StructVectorGetBoxValues,(mx->hx,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,yy));
241: VecRestoreArray(y,&yy);
242: return(0);
243: }
245: PetscErrorCode MatAssemblyEnd_HYPREStruct(Mat mat,MatAssemblyType mode)
246: {
247: Mat_HYPREStruct *ex = (Mat_HYPREStruct*) mat->data;
248: PetscErrorCode ierr;
251: PetscStackCallStandard(HYPRE_StructMatrixAssemble,(ex->hmat));
252: /* PetscStackCallStandard(HYPRE_StructMatrixPrint,("dummy",ex->hmat,0)); */
253: return(0);
254: }
256: PetscErrorCode MatZeroEntries_HYPREStruct(Mat mat)
257: {
259: /* before the DMDA is set to the matrix the zero doesn't need to do anything */
260: return(0);
261: }
263: PetscErrorCode MatDestroy_HYPREStruct(Mat mat)
264: {
265: Mat_HYPREStruct *ex = (Mat_HYPREStruct*) mat->data;
266: PetscErrorCode ierr;
269: PetscStackCallStandard(HYPRE_StructMatrixDestroy,(ex->hmat));
270: PetscStackCallStandard(HYPRE_StructVectorDestroy,(ex->hx));
271: PetscStackCallStandard(HYPRE_StructVectorDestroy,(ex->hb));
272: PetscObjectDereference((PetscObject)ex->da);
273: MPI_Comm_free(&(ex->hcomm));
274: PetscFree(ex);
275: return(0);
276: }
278: PETSC_EXTERN PetscErrorCode MatCreate_HYPREStruct(Mat B)
279: {
280: Mat_HYPREStruct *ex;
281: PetscErrorCode ierr;
284: PetscNewLog(B,&ex);
285: B->data = (void*)ex;
286: B->rmap->bs = 1;
287: B->assembled = PETSC_FALSE;
289: B->insertmode = NOT_SET_VALUES;
291: B->ops->assemblyend = MatAssemblyEnd_HYPREStruct;
292: B->ops->mult = MatMult_HYPREStruct;
293: B->ops->zeroentries = MatZeroEntries_HYPREStruct;
294: B->ops->destroy = MatDestroy_HYPREStruct;
295: B->ops->setup = MatSetUp_HYPREStruct;
297: ex->needsinitialization = PETSC_TRUE;
299: MPI_Comm_dup(PetscObjectComm((PetscObject)B),&(ex->hcomm));
300: PetscObjectChangeTypeName((PetscObject)B,MATHYPRESTRUCT);
301: return(0);
302: }
304: /*MC
305: MATHYPRESSTRUCT - MATHYPRESSTRUCT = "hypresstruct" - A matrix type to be used for parallel sparse matrices
306: based on the hypre HYPRE_SStructMatrix.
309: Level: intermediate
311: Notes:
312: Unlike hypre's general semi-struct object consisting of a collection of structured-grid objects and unstructured
313: grid objects, we restrict the semi-struct objects to consist of only structured-grid components.
315: 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
316: be defined by a DMDA.
318: The matrix needs a DMDA associated with it by either a call to MatSetDM() or if the matrix is obtained from DMCreateMatrix()
320: M*/
322: PetscErrorCode MatSetValuesLocal_HYPRESStruct_3d(Mat mat,PetscInt nrow,const PetscInt irow[],PetscInt ncol,const PetscInt icol[],const PetscScalar y[],InsertMode addv)
323: {
324: PetscErrorCode ierr;
325: PetscInt i,j,stencil,index[3];
326: const PetscScalar *values = y;
327: Mat_HYPRESStruct *ex = (Mat_HYPRESStruct*) mat->data;
329: PetscInt part = 0; /* Petsc sstruct interface only allows 1 part */
330: PetscInt ordering;
331: PetscInt grid_rank, to_grid_rank;
332: PetscInt var_type, to_var_type;
333: PetscInt to_var_entry = 0;
334: PetscInt nvars= ex->nvars;
335: PetscInt row,*entries;
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 SETERRQ3(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] = ex->xs + (row % ex->nx);
374: index[1] = ex->ys + ((row/ex->nx) % ex->ny);
375: index[2] = ex->zs + (row/(ex->nxny));
377: if (addv == ADD_VALUES) {
378: PetscStackCallStandard(HYPRE_SStructMatrixAddToValues,(ex->ss_mat,part,(HYPRE_Int *)index,var_type,ncol,(HYPRE_Int *)entries,(PetscScalar*)values));
379: } else {
380: PetscStackCallStandard(HYPRE_SStructMatrixSetValues,(ex->ss_mat,part,(HYPRE_Int *)index,var_type,ncol,(HYPRE_Int *)entries,(PetscScalar*)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 SETERRQ3(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] = ex->xs + (row % ex->nx);
416: index[1] = ex->ys + ((row/ex->nx) % ex->ny);
417: index[2] = ex->zs + (row/(ex->nxny));
419: if (addv == ADD_VALUES) {
420: PetscStackCallStandard(HYPRE_SStructMatrixAddToValues,(ex->ss_mat,part,(HYPRE_Int *)index,var_type,ncol,(HYPRE_Int *)entries,(PetscScalar*)values));
421: } else {
422: PetscStackCallStandard(HYPRE_SStructMatrixSetValues,(ex->ss_mat,part,(HYPRE_Int *)index,var_type,ncol,(HYPRE_Int *)entries,(PetscScalar*)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: PetscErrorCode ierr;
434: PetscInt i,index[3];
435: PetscScalar **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,*entries;
446: if (x && b) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"No support");
447: PetscMalloc1(7*nvars,&entries);
449: PetscMalloc1(nvars,&values);
450: PetscMalloc1(7*nvars*nvars,&values[0]);
451: for (i=1; i<nvars; i++) {
452: values[i] = values[i-1] + nvars*7;
453: }
455: for (i=0; i< nvars; i++) {
456: PetscMemzero(values[i],nvars*7*sizeof(PetscScalar));
457: *(values[i]+3) = d;
458: }
460: for (i= 0; i< nvars*7; i++) entries[i] = i;
462: if (!ordering) {
463: for (i=0; i<nrow; i++) {
464: grid_rank = irow[i] / nvars;
465: var_type =(irow[i] % nvars);
467: row = ex->gindices[grid_rank] - ex->rstart;
468: index[0] = ex->xs + (row % ex->nx);
469: index[1] = ex->ys + ((row/ex->nx) % ex->ny);
470: index[2] = ex->zs + (row/(ex->nxny));
471: PetscStackCallStandard(HYPRE_SStructMatrixSetValues,(ex->ss_mat,part,(HYPRE_Int *)index,var_type,7*nvars,(HYPRE_Int *)entries,values[var_type]));
472: }
473: } else {
474: for (i=0; i<nrow; i++) {
475: var_type = irow[i]/(ex->gnxgnygnz);
476: grid_rank= irow[i] - var_type*(ex->gnxgnygnz);
478: row = ex->gindices[grid_rank] - ex->rstart;
479: index[0] = ex->xs + (row % ex->nx);
480: index[1] = ex->ys + ((row/ex->nx) % ex->ny);
481: index[2] = ex->zs + (row/(ex->nxny));
482: PetscStackCallStandard(HYPRE_SStructMatrixSetValues,(ex->ss_mat,part,(HYPRE_Int *)index,var_type,7*nvars,(HYPRE_Int *)entries,values[var_type]));
483: }
484: }
485: PetscStackCallStandard(HYPRE_SStructMatrixAssemble,(ex->ss_mat));
486: PetscFree(values[0]);
487: PetscFree(values);
488: PetscFree(entries);
489: return(0);
490: }
492: PetscErrorCode MatZeroEntries_HYPRESStruct_3d(Mat mat)
493: {
494: PetscErrorCode ierr;
495: Mat_HYPRESStruct *ex = (Mat_HYPRESStruct*) mat->data;
496: PetscInt nvars= ex->nvars;
497: PetscInt size;
498: PetscInt part= 0; /* only one part */
501: 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);
502: {
503: PetscInt i,*entries,iupper[3],ilower[3];
504: PetscScalar *values;
507: for (i= 0; i< 3; i++) {
508: ilower[i]= ex->hbox.imin[i];
509: iupper[i]= ex->hbox.imax[i];
510: }
512: PetscMalloc2(nvars*7,&entries,nvars*7*size,&values);
513: for (i= 0; i< nvars*7; i++) entries[i]= i;
514: PetscMemzero(values,nvars*7*size*sizeof(PetscScalar));
516: for (i= 0; i< nvars; i++) {
517: PetscStackCallStandard(HYPRE_SStructMatrixSetBoxValues,(ex->ss_mat,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,nvars*7,(HYPRE_Int *)entries,values));
518: }
519: PetscFree2(entries,values);
520: }
521: PetscStackCallStandard(HYPRE_SStructMatrixAssemble,(ex->ss_mat));
522: return(0);
523: }
526: static PetscErrorCode MatSetUp_HYPRESStruct(Mat mat)
527: {
528: PetscErrorCode ierr;
529: Mat_HYPRESStruct *ex = (Mat_HYPRESStruct*) mat->data;
530: PetscInt dim,dof,sw[3],nx,ny,nz;
531: PetscInt ilower[3],iupper[3],ssize,i;
532: DMBoundaryType px,py,pz;
533: DMDAStencilType st;
534: PetscInt nparts= 1; /* assuming only one part */
535: PetscInt part = 0;
536: ISLocalToGlobalMapping ltog;
537: DM da;
540: MatGetDM(mat,(DM*)&da);
541: ex->da = da;
542: PetscObjectReference((PetscObject)da);
544: DMDAGetInfo(ex->da,&dim,0,0,0,0,0,0,&dof,&sw[0],&px,&py,&pz,&st);
545: DMDAGetCorners(ex->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);
546: iupper[0] += ilower[0] - 1;
547: iupper[1] += ilower[1] - 1;
548: iupper[2] += ilower[2] - 1;
549: /* the hypre_Box is used to zero out the matrix entries in MatZeroValues() */
550: ex->hbox.imin[0] = ilower[0];
551: ex->hbox.imin[1] = ilower[1];
552: ex->hbox.imin[2] = ilower[2];
553: ex->hbox.imax[0] = iupper[0];
554: ex->hbox.imax[1] = iupper[1];
555: ex->hbox.imax[2] = iupper[2];
557: ex->dofs_order = 0;
559: /* assuming that the same number of dofs on each gridpoint. Also assume all cell-centred based */
560: ex->nvars= dof;
562: /* create the hypre grid object and set its information */
563: if (px || py || pz) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Ask us to add periodic support by calling HYPRE_SStructGridSetPeriodic()");
564: PetscStackCallStandard(HYPRE_SStructGridCreate,(ex->hcomm,dim,nparts,&ex->ss_grid));
565: PetscStackCallStandard(HYPRE_SStructGridSetExtents,(ex->ss_grid,part,ex->hbox.imin,ex->hbox.imax));
566: {
567: HYPRE_SStructVariable *vartypes;
568: PetscMalloc1(ex->nvars,&vartypes);
569: for (i= 0; i< ex->nvars; i++) vartypes[i]= HYPRE_SSTRUCT_VARIABLE_CELL;
570: PetscStackCallStandard(HYPRE_SStructGridSetVariables,(ex->ss_grid, part, ex->nvars,vartypes));
571: PetscFree(vartypes);
572: }
573: PetscStackCallStandard(HYPRE_SStructGridAssemble,(ex->ss_grid));
575: sw[1] = sw[0];
576: sw[2] = sw[1];
577: /* PetscStackCallStandard(HYPRE_SStructGridSetNumGhost,(ex->ss_grid,sw)); */
579: /* create the hypre stencil object and set its information */
580: if (sw[0] > 1) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Ask us to add support for wider stencils");
581: if (st == DMDA_STENCIL_BOX) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Ask us to add support for box stencils");
583: if (dim == 1) {
584: PetscInt offsets[3][1] = {{-1},{0},{1}};
585: PetscInt j, cnt;
587: ssize = 3*(ex->nvars);
588: PetscStackCallStandard(HYPRE_SStructStencilCreate,(dim,ssize,&ex->ss_stencil));
589: cnt= 0;
590: for (i = 0; i < (ex->nvars); i++) {
591: for (j = 0; j < 3; j++) {
592: PetscStackCallStandard(HYPRE_SStructStencilSetEntry,(ex->ss_stencil, cnt, (HYPRE_Int *)offsets[j], i));
593: cnt++;
594: }
595: }
596: } else if (dim == 2) {
597: PetscInt offsets[5][2] = {{0,-1},{-1,0},{0,0},{1,0},{0,1}};
598: PetscInt j, cnt;
600: ssize = 5*(ex->nvars);
601: PetscStackCallStandard(HYPRE_SStructStencilCreate,(dim,ssize,&ex->ss_stencil));
602: cnt= 0;
603: for (i= 0; i< (ex->nvars); i++) {
604: for (j= 0; j< 5; j++) {
605: PetscStackCallStandard(HYPRE_SStructStencilSetEntry,(ex->ss_stencil, cnt, (HYPRE_Int *)offsets[j], i));
606: cnt++;
607: }
608: }
609: } else if (dim == 3) {
610: 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}};
611: PetscInt j, cnt;
613: ssize = 7*(ex->nvars);
614: PetscStackCallStandard(HYPRE_SStructStencilCreate,(dim,ssize,&ex->ss_stencil));
615: cnt= 0;
616: for (i= 0; i< (ex->nvars); i++) {
617: for (j= 0; j< 7; j++) {
618: PetscStackCallStandard(HYPRE_SStructStencilSetEntry,(ex->ss_stencil, cnt, (HYPRE_Int *)offsets[j], i));
619: cnt++;
620: }
621: }
622: }
624: /* create the HYPRE graph */
625: PetscStackCallStandard(HYPRE_SStructGraphCreate,(ex->hcomm, ex->ss_grid, &(ex->ss_graph)));
627: /* set the stencil graph. Note that each variable has the same graph. This means that each
628: variable couples to all the other variable and with the same stencil pattern. */
629: for (i= 0; i< (ex->nvars); i++) {
630: PetscStackCallStandard(HYPRE_SStructGraphSetStencil,(ex->ss_graph,part,i,ex->ss_stencil));
631: }
632: PetscStackCallStandard(HYPRE_SStructGraphAssemble,(ex->ss_graph));
634: /* create the HYPRE sstruct vectors for rhs and solution */
635: PetscStackCallStandard(HYPRE_SStructVectorCreate,(ex->hcomm,ex->ss_grid,&ex->ss_b));
636: PetscStackCallStandard(HYPRE_SStructVectorCreate,(ex->hcomm,ex->ss_grid,&ex->ss_x));
637: PetscStackCallStandard(HYPRE_SStructVectorInitialize,(ex->ss_b));
638: PetscStackCallStandard(HYPRE_SStructVectorInitialize,(ex->ss_x));
639: PetscStackCallStandard(HYPRE_SStructVectorAssemble,(ex->ss_b));
640: PetscStackCallStandard(HYPRE_SStructVectorAssemble,(ex->ss_x));
642: /* create the hypre matrix object and set its information */
643: PetscStackCallStandard(HYPRE_SStructMatrixCreate,(ex->hcomm,ex->ss_graph,&ex->ss_mat));
644: PetscStackCallStandard(HYPRE_SStructGridDestroy,(ex->ss_grid));
645: PetscStackCallStandard(HYPRE_SStructStencilDestroy,(ex->ss_stencil));
646: if (ex->needsinitialization) {
647: PetscStackCallStandard(HYPRE_SStructMatrixInitialize,(ex->ss_mat));
648: ex->needsinitialization = PETSC_FALSE;
649: }
651: /* set the global and local sizes of the matrix */
652: DMDAGetCorners(da,0,0,0,&nx,&ny,&nz);
653: MatSetSizes(mat,dof*nx*ny*nz,dof*nx*ny*nz,PETSC_DECIDE,PETSC_DECIDE);
654: PetscLayoutSetBlockSize(mat->rmap,dof);
655: PetscLayoutSetBlockSize(mat->cmap,dof);
656: PetscLayoutSetUp(mat->rmap);
657: PetscLayoutSetUp(mat->cmap);
658: mat->preallocated = PETSC_TRUE;
660: if (dim == 3) {
661: mat->ops->setvalueslocal = MatSetValuesLocal_HYPRESStruct_3d;
662: mat->ops->zerorowslocal = MatZeroRowsLocal_HYPRESStruct_3d;
663: mat->ops->zeroentries = MatZeroEntries_HYPRESStruct_3d;
665: /* MatZeroEntries_HYPRESStruct_3d(mat); */
666: } else SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Only support for 3d DMDA currently");
668: /* get values that will be used repeatedly in MatSetValuesLocal() and MatZeroRowsLocal() repeatedly */
669: MatGetOwnershipRange(mat,&ex->rstart,NULL);
670: DMGetLocalToGlobalMapping(ex->da,<og);
671: ISLocalToGlobalMappingGetIndices(ltog, (const PetscInt **) &ex->gindices);
672: DMDAGetGhostCorners(ex->da,0,0,0,&ex->gnx,&ex->gnxgny,&ex->gnxgnygnz);
674: ex->gnxgny *= ex->gnx;
675: ex->gnxgnygnz *= ex->gnxgny;
677: DMDAGetCorners(ex->da,&ex->xs,&ex->ys,&ex->zs,&ex->nx,&ex->ny,&ex->nz);
679: ex->nxny = ex->nx*ex->ny;
680: ex->nxnynz = ex->nz*ex->nxny;
681: return(0);
682: }
684: PetscErrorCode MatMult_HYPRESStruct(Mat A,Vec x,Vec y)
685: {
686: PetscErrorCode ierr;
687: const PetscScalar *xx;
688: PetscScalar *yy;
689: PetscInt ilower[3],iupper[3];
690: Mat_HYPRESStruct *mx = (Mat_HYPRESStruct*)(A->data);
691: PetscInt ordering= mx->dofs_order;
692: PetscInt nvars = mx->nvars;
693: PetscInt part = 0;
694: PetscInt size;
695: PetscInt i;
698: DMDAGetCorners(mx->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);
699: iupper[0] += ilower[0] - 1;
700: iupper[1] += ilower[1] - 1;
701: iupper[2] += ilower[2] - 1;
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,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,(PetscScalar*)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,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,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,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,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,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,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;
762: PetscErrorCode ierr;
765: PetscStackCallStandard(HYPRE_SStructMatrixAssemble,(ex->ss_mat));
766: return(0);
767: }
769: PetscErrorCode MatZeroEntries_HYPRESStruct(Mat mat)
770: {
772: /* before the DMDA is set to the matrix the zero doesn't need to do anything */
773: return(0);
774: }
776: PetscErrorCode MatDestroy_HYPRESStruct(Mat mat)
777: {
778: Mat_HYPRESStruct *ex = (Mat_HYPRESStruct*) mat->data;
779: PetscErrorCode ierr;
780: ISLocalToGlobalMapping ltog;
783: DMGetLocalToGlobalMapping(ex->da,<og);
784: ISLocalToGlobalMappingRestoreIndices(ltog, (const PetscInt **) &ex->gindices);
785: PetscStackCallStandard(HYPRE_SStructGraphDestroy,(ex->ss_graph));
786: PetscStackCallStandard(HYPRE_SStructMatrixDestroy,(ex->ss_mat));
787: PetscStackCallStandard(HYPRE_SStructVectorDestroy,(ex->ss_x));
788: PetscStackCallStandard(HYPRE_SStructVectorDestroy,(ex->ss_b));
789: PetscObjectDereference((PetscObject)ex->da);
790: MPI_Comm_free(&(ex->hcomm));
791: PetscFree(ex);
792: return(0);
793: }
795: PETSC_EXTERN PetscErrorCode MatCreate_HYPRESStruct(Mat B)
796: {
797: Mat_HYPRESStruct *ex;
798: PetscErrorCode ierr;
801: PetscNewLog(B,&ex);
802: B->data = (void*)ex;
803: B->rmap->bs = 1;
804: B->assembled = PETSC_FALSE;
806: B->insertmode = NOT_SET_VALUES;
808: B->ops->assemblyend = MatAssemblyEnd_HYPRESStruct;
809: B->ops->mult = MatMult_HYPRESStruct;
810: B->ops->zeroentries = MatZeroEntries_HYPRESStruct;
811: B->ops->destroy = MatDestroy_HYPRESStruct;
812: B->ops->setup = MatSetUp_HYPRESStruct;
814: ex->needsinitialization = PETSC_TRUE;
816: MPI_Comm_dup(PetscObjectComm((PetscObject)B),&(ex->hcomm));
817: PetscObjectChangeTypeName((PetscObject)B,MATHYPRESSTRUCT);
818: return(0);
819: }