Actual source code: mhyp.c
petsc-3.11.4 2019-09-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: 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: PetscInt indices[7] = {0,1,2,3,4,5,6};
94: Mat_HYPREStruct *ex = (Mat_HYPREStruct*) mat->data;
97: /* hypre has no public interface to do this */
98: PetscStackCallStandard(hypre_StructMatrixClearBoxValues,(ex->hmat,&ex->hbox,7,(HYPRE_Int *)indices,0,1));
99: PetscStackCallStandard(HYPRE_StructMatrixAssemble,(ex->hmat));
100: return(0);
101: }
103: static PetscErrorCode MatSetUp_HYPREStruct(Mat mat)
104: {
105: PetscErrorCode ierr;
106: Mat_HYPREStruct *ex = (Mat_HYPREStruct*) mat->data;
107: PetscInt dim,dof,sw[6],nx,ny,nz,ilower[3],iupper[3],ssize,i;
108: DMBoundaryType px,py,pz;
109: DMDAStencilType st;
110: ISLocalToGlobalMapping ltog;
111: DM da;
114: MatGetDM(mat,(DM*)&da);
115: ex->da = da;
116: PetscObjectReference((PetscObject)da);
118: DMDAGetInfo(ex->da,&dim,0,0,0,0,0,0,&dof,&sw[0],&px,&py,&pz,&st);
119: DMDAGetCorners(ex->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);
120: iupper[0] += ilower[0] - 1;
121: iupper[1] += ilower[1] - 1;
122: iupper[2] += ilower[2] - 1;
124: /* the hypre_Box is used to zero out the matrix entries in MatZeroValues() */
125: ex->hbox.imin[0] = ilower[0];
126: ex->hbox.imin[1] = ilower[1];
127: ex->hbox.imin[2] = ilower[2];
128: ex->hbox.imax[0] = iupper[0];
129: ex->hbox.imax[1] = iupper[1];
130: ex->hbox.imax[2] = iupper[2];
132: /* create the hypre grid object and set its information */
133: if (dof > 1) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Currently only support for scalar problems");
134: if (px || py || pz) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Ask us to add periodic support by calling HYPRE_StructGridSetPeriodic()");
135: PetscStackCallStandard(HYPRE_StructGridCreate,(ex->hcomm,dim,&ex->hgrid));
136: PetscStackCallStandard(HYPRE_StructGridSetExtents,(ex->hgrid,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper));
137: PetscStackCallStandard(HYPRE_StructGridAssemble,(ex->hgrid));
139: sw[5] = sw[4] = sw[3] = sw[2] = sw[1] = sw[0];
140: PetscStackCallStandard(HYPRE_StructGridSetNumGhost,(ex->hgrid,(HYPRE_Int *)sw));
142: /* create the hypre stencil object and set its information */
143: if (sw[0] > 1) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Ask us to add support for wider stencils");
144: if (st == DMDA_STENCIL_BOX) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Ask us to add support for box stencils");
145: if (dim == 1) {
146: PetscInt offsets[3][1] = {{-1},{0},{1}};
147: ssize = 3;
148: PetscStackCallStandard(HYPRE_StructStencilCreate,(dim,ssize,&ex->hstencil));
149: for (i=0; i<ssize; i++) {
150: PetscStackCallStandard(HYPRE_StructStencilSetElement,(ex->hstencil,i,(HYPRE_Int *)offsets[i]));
151: }
152: } else if (dim == 2) {
153: PetscInt offsets[5][2] = {{0,-1},{-1,0},{0,0},{1,0},{0,1}};
154: ssize = 5;
155: PetscStackCallStandard(HYPRE_StructStencilCreate,(dim,ssize,&ex->hstencil));
156: for (i=0; i<ssize; i++) {
157: PetscStackCallStandard(HYPRE_StructStencilSetElement,(ex->hstencil,i,(HYPRE_Int *)offsets[i]));
158: }
159: } else if (dim == 3) {
160: 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}};
161: ssize = 7;
162: PetscStackCallStandard(HYPRE_StructStencilCreate,(dim,ssize,&ex->hstencil));
163: for (i=0; i<ssize; i++) {
164: PetscStackCallStandard(HYPRE_StructStencilSetElement,(ex->hstencil,i,(HYPRE_Int *)offsets[i]));
165: }
166: }
168: /* create the HYPRE vector for rhs and solution */
169: PetscStackCallStandard(HYPRE_StructVectorCreate,(ex->hcomm,ex->hgrid,&ex->hb));
170: PetscStackCallStandard(HYPRE_StructVectorCreate,(ex->hcomm,ex->hgrid,&ex->hx));
171: PetscStackCallStandard(HYPRE_StructVectorInitialize,(ex->hb));
172: PetscStackCallStandard(HYPRE_StructVectorInitialize,(ex->hx));
173: PetscStackCallStandard(HYPRE_StructVectorAssemble,(ex->hb));
174: PetscStackCallStandard(HYPRE_StructVectorAssemble,(ex->hx));
176: /* create the hypre matrix object and set its information */
177: PetscStackCallStandard(HYPRE_StructMatrixCreate,(ex->hcomm,ex->hgrid,ex->hstencil,&ex->hmat));
178: PetscStackCallStandard(HYPRE_StructGridDestroy,(ex->hgrid));
179: PetscStackCallStandard(HYPRE_StructStencilDestroy,(ex->hstencil));
180: if (ex->needsinitialization) {
181: PetscStackCallStandard(HYPRE_StructMatrixInitialize,(ex->hmat));
182: ex->needsinitialization = PETSC_FALSE;
183: }
185: /* set the global and local sizes of the matrix */
186: DMDAGetCorners(da,0,0,0,&nx,&ny,&nz);
187: MatSetSizes(mat,dof*nx*ny*nz,dof*nx*ny*nz,PETSC_DECIDE,PETSC_DECIDE);
188: PetscLayoutSetBlockSize(mat->rmap,1);
189: PetscLayoutSetBlockSize(mat->cmap,1);
190: PetscLayoutSetUp(mat->rmap);
191: PetscLayoutSetUp(mat->cmap);
192: mat->preallocated = PETSC_TRUE;
194: if (dim == 3) {
195: mat->ops->setvalueslocal = MatSetValuesLocal_HYPREStruct_3d;
196: mat->ops->zerorowslocal = MatZeroRowsLocal_HYPREStruct_3d;
197: mat->ops->zeroentries = MatZeroEntries_HYPREStruct_3d;
199: /* MatZeroEntries_HYPREStruct_3d(mat); */
200: } else SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Only support for 3d DMDA currently");
202: /* get values that will be used repeatedly in MatSetValuesLocal() and MatZeroRowsLocal() repeatedly */
203: MatGetOwnershipRange(mat,&ex->rstart,NULL);
204: DMGetLocalToGlobalMapping(ex->da,<og);
205: ISLocalToGlobalMappingGetIndices(ltog, (const PetscInt **) &ex->gindices);
206: DMDAGetGhostCorners(ex->da,0,0,0,&ex->gnx,&ex->gnxgny,0);
207: ex->gnxgny *= ex->gnx;
208: DMDAGetCorners(ex->da,&ex->xs,&ex->ys,&ex->zs,&ex->nx,&ex->ny,0);
209: ex->nxny = ex->nx*ex->ny;
210: return(0);
211: }
213: PetscErrorCode MatMult_HYPREStruct(Mat A,Vec x,Vec y)
214: {
215: PetscErrorCode ierr;
216: const PetscScalar *xx;
217: PetscScalar *yy;
218: PetscInt ilower[3],iupper[3];
219: Mat_HYPREStruct *mx = (Mat_HYPREStruct*)(A->data);
222: DMDAGetCorners(mx->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);
224: iupper[0] += ilower[0] - 1;
225: iupper[1] += ilower[1] - 1;
226: iupper[2] += ilower[2] - 1;
228: /* copy x values over to hypre */
229: PetscStackCallStandard(HYPRE_StructVectorSetConstantValues,(mx->hb,0.0));
230: VecGetArrayRead(x,&xx);
231: PetscStackCallStandard(HYPRE_StructVectorSetBoxValues,(mx->hb,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,(PetscScalar*)xx));
232: VecRestoreArrayRead(x,&xx);
233: PetscStackCallStandard(HYPRE_StructVectorAssemble,(mx->hb));
234: PetscStackCallStandard(HYPRE_StructMatrixMatvec,(1.0,mx->hmat,mx->hb,0.0,mx->hx));
236: /* copy solution values back to PETSc */
237: VecGetArray(y,&yy);
238: PetscStackCallStandard(HYPRE_StructVectorGetBoxValues,(mx->hx,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,yy));
239: VecRestoreArray(y,&yy);
240: return(0);
241: }
243: PetscErrorCode MatAssemblyEnd_HYPREStruct(Mat mat,MatAssemblyType mode)
244: {
245: Mat_HYPREStruct *ex = (Mat_HYPREStruct*) mat->data;
248: PetscStackCallStandard(HYPRE_StructMatrixAssemble,(ex->hmat));
249: /* PetscStackCallStandard(HYPRE_StructMatrixPrint,("dummy",ex->hmat,0)); */
250: return(0);
251: }
253: PetscErrorCode MatZeroEntries_HYPREStruct(Mat mat)
254: {
256: /* before the DMDA is set to the matrix the zero doesn't need to do anything */
257: return(0);
258: }
260: PetscErrorCode MatDestroy_HYPREStruct(Mat mat)
261: {
262: Mat_HYPREStruct *ex = (Mat_HYPREStruct*) mat->data;
263: PetscErrorCode ierr;
266: PetscStackCallStandard(HYPRE_StructMatrixDestroy,(ex->hmat));
267: PetscStackCallStandard(HYPRE_StructVectorDestroy,(ex->hx));
268: PetscStackCallStandard(HYPRE_StructVectorDestroy,(ex->hb));
269: PetscObjectDereference((PetscObject)ex->da);
270: MPI_Comm_free(&(ex->hcomm));
271: PetscFree(ex);
272: return(0);
273: }
275: PETSC_EXTERN PetscErrorCode MatCreate_HYPREStruct(Mat B)
276: {
277: Mat_HYPREStruct *ex;
278: PetscErrorCode ierr;
281: PetscNewLog(B,&ex);
282: B->data = (void*)ex;
283: B->rmap->bs = 1;
284: B->assembled = PETSC_FALSE;
286: B->insertmode = NOT_SET_VALUES;
288: B->ops->assemblyend = MatAssemblyEnd_HYPREStruct;
289: B->ops->mult = MatMult_HYPREStruct;
290: B->ops->zeroentries = MatZeroEntries_HYPREStruct;
291: B->ops->destroy = MatDestroy_HYPREStruct;
292: B->ops->setup = MatSetUp_HYPREStruct;
294: ex->needsinitialization = PETSC_TRUE;
296: MPI_Comm_dup(PetscObjectComm((PetscObject)B),&(ex->hcomm));
297: PetscObjectChangeTypeName((PetscObject)B,MATHYPRESTRUCT);
298: return(0);
299: }
301: /*MC
302: MATHYPRESSTRUCT - MATHYPRESSTRUCT = "hypresstruct" - A matrix type to be used for parallel sparse matrices
303: based on the hypre HYPRE_SStructMatrix.
306: Level: intermediate
308: Notes:
309: Unlike hypre's general semi-struct object consisting of a collection of structured-grid objects and unstructured
310: grid objects, we restrict the semi-struct objects to consist of only structured-grid components.
312: 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
313: be defined by a DMDA.
315: The matrix needs a DMDA associated with it by either a call to MatSetDM() or if the matrix is obtained from DMCreateMatrix()
317: M*/
319: PetscErrorCode MatSetValuesLocal_HYPRESStruct_3d(Mat mat,PetscInt nrow,const PetscInt irow[],PetscInt ncol,const PetscInt icol[],const PetscScalar y[],InsertMode addv)
320: {
321: PetscErrorCode ierr;
322: PetscInt i,j,stencil,index[3];
323: const PetscScalar *values = y;
324: Mat_HYPRESStruct *ex = (Mat_HYPRESStruct*) mat->data;
326: PetscInt part = 0; /* Petsc sstruct interface only allows 1 part */
327: PetscInt ordering;
328: PetscInt grid_rank, to_grid_rank;
329: PetscInt var_type, to_var_type;
330: PetscInt to_var_entry = 0;
331: PetscInt nvars= ex->nvars;
332: PetscInt row,*entries;
335: PetscMalloc1(7*nvars,&entries);
336: ordering = ex->dofs_order; /* ordering= 0 nodal ordering
337: 1 variable ordering */
338: /* stencil entries are ordered by variables: var0_stencil0, var0_stencil1, ..., var0_stencil6, var1_stencil0, var1_stencil1, ... */
339: if (!ordering) {
340: for (i=0; i<nrow; i++) {
341: grid_rank= irow[i]/nvars;
342: var_type = (irow[i] % nvars);
344: for (j=0; j<ncol; j++) {
345: to_grid_rank= icol[j]/nvars;
346: to_var_type = (icol[j] % nvars);
348: to_var_entry= to_var_entry*7;
349: entries[j] = to_var_entry;
351: stencil = to_grid_rank-grid_rank;
352: if (!stencil) {
353: entries[j] += 3;
354: } else if (stencil == -1) {
355: entries[j] += 2;
356: } else if (stencil == 1) {
357: entries[j] += 4;
358: } else if (stencil == -ex->gnx) {
359: entries[j] += 1;
360: } else if (stencil == ex->gnx) {
361: entries[j] += 5;
362: } else if (stencil == -ex->gnxgny) {
363: entries[j] += 0;
364: } else if (stencil == ex->gnxgny) {
365: entries[j] += 6;
366: } else SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Local row %D local column %D have bad stencil %D",irow[i],icol[j],stencil);
367: }
369: row = ex->gindices[grid_rank] - ex->rstart;
370: index[0] = ex->xs + (row % ex->nx);
371: index[1] = ex->ys + ((row/ex->nx) % ex->ny);
372: index[2] = ex->zs + (row/(ex->nxny));
374: if (addv == ADD_VALUES) {
375: PetscStackCallStandard(HYPRE_SStructMatrixAddToValues,(ex->ss_mat,part,(HYPRE_Int *)index,var_type,ncol,(HYPRE_Int *)entries,(PetscScalar*)values));
376: } else {
377: PetscStackCallStandard(HYPRE_SStructMatrixSetValues,(ex->ss_mat,part,(HYPRE_Int *)index,var_type,ncol,(HYPRE_Int *)entries,(PetscScalar*)values));
378: }
379: values += ncol;
380: }
381: } else {
382: for (i=0; i<nrow; i++) {
383: var_type = irow[i]/(ex->gnxgnygnz);
384: grid_rank= irow[i] - var_type*(ex->gnxgnygnz);
386: for (j=0; j<ncol; j++) {
387: to_var_type = icol[j]/(ex->gnxgnygnz);
388: to_grid_rank= icol[j] - to_var_type*(ex->gnxgnygnz);
390: to_var_entry= to_var_entry*7;
391: entries[j] = to_var_entry;
393: stencil = to_grid_rank-grid_rank;
394: if (!stencil) {
395: entries[j] += 3;
396: } else if (stencil == -1) {
397: entries[j] += 2;
398: } else if (stencil == 1) {
399: entries[j] += 4;
400: } else if (stencil == -ex->gnx) {
401: entries[j] += 1;
402: } else if (stencil == ex->gnx) {
403: entries[j] += 5;
404: } else if (stencil == -ex->gnxgny) {
405: entries[j] += 0;
406: } else if (stencil == ex->gnxgny) {
407: entries[j] += 6;
408: } else SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Local row %D local column %D have bad stencil %D",irow[i],icol[j],stencil);
409: }
411: row = ex->gindices[grid_rank] - ex->rstart;
412: index[0] = ex->xs + (row % ex->nx);
413: index[1] = ex->ys + ((row/ex->nx) % ex->ny);
414: index[2] = ex->zs + (row/(ex->nxny));
416: if (addv == ADD_VALUES) {
417: PetscStackCallStandard(HYPRE_SStructMatrixAddToValues,(ex->ss_mat,part,(HYPRE_Int *)index,var_type,ncol,(HYPRE_Int *)entries,(PetscScalar*)values));
418: } else {
419: PetscStackCallStandard(HYPRE_SStructMatrixSetValues,(ex->ss_mat,part,(HYPRE_Int *)index,var_type,ncol,(HYPRE_Int *)entries,(PetscScalar*)values));
420: }
421: values += ncol;
422: }
423: }
424: PetscFree(entries);
425: return(0);
426: }
428: PetscErrorCode MatZeroRowsLocal_HYPRESStruct_3d(Mat mat,PetscInt nrow,const PetscInt irow[],PetscScalar d,Vec x,Vec b)
429: {
430: PetscErrorCode ierr;
431: PetscInt i,index[3];
432: PetscScalar **values;
433: Mat_HYPRESStruct *ex = (Mat_HYPRESStruct*) mat->data;
435: PetscInt part = 0; /* Petsc sstruct interface only allows 1 part */
436: PetscInt ordering = ex->dofs_order;
437: PetscInt grid_rank;
438: PetscInt var_type;
439: PetscInt nvars= ex->nvars;
440: PetscInt row,*entries;
443: if (x && b) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"No support");
444: PetscMalloc1(7*nvars,&entries);
446: PetscMalloc1(nvars,&values);
447: PetscMalloc1(7*nvars*nvars,&values[0]);
448: for (i=1; i<nvars; i++) {
449: values[i] = values[i-1] + nvars*7;
450: }
452: for (i=0; i< nvars; i++) {
453: PetscMemzero(values[i],nvars*7*sizeof(PetscScalar));
454: *(values[i]+3) = d;
455: }
457: for (i= 0; i< nvars*7; i++) entries[i] = i;
459: if (!ordering) {
460: for (i=0; i<nrow; i++) {
461: grid_rank = irow[i] / nvars;
462: var_type =(irow[i] % nvars);
464: row = ex->gindices[grid_rank] - ex->rstart;
465: index[0] = ex->xs + (row % ex->nx);
466: index[1] = ex->ys + ((row/ex->nx) % ex->ny);
467: index[2] = ex->zs + (row/(ex->nxny));
468: PetscStackCallStandard(HYPRE_SStructMatrixSetValues,(ex->ss_mat,part,(HYPRE_Int *)index,var_type,7*nvars,(HYPRE_Int *)entries,values[var_type]));
469: }
470: } else {
471: for (i=0; i<nrow; i++) {
472: var_type = irow[i]/(ex->gnxgnygnz);
473: grid_rank= irow[i] - var_type*(ex->gnxgnygnz);
475: row = ex->gindices[grid_rank] - ex->rstart;
476: index[0] = ex->xs + (row % ex->nx);
477: index[1] = ex->ys + ((row/ex->nx) % ex->ny);
478: index[2] = ex->zs + (row/(ex->nxny));
479: PetscStackCallStandard(HYPRE_SStructMatrixSetValues,(ex->ss_mat,part,(HYPRE_Int *)index,var_type,7*nvars,(HYPRE_Int *)entries,values[var_type]));
480: }
481: }
482: PetscStackCallStandard(HYPRE_SStructMatrixAssemble,(ex->ss_mat));
483: PetscFree(values[0]);
484: PetscFree(values);
485: PetscFree(entries);
486: return(0);
487: }
489: PetscErrorCode MatZeroEntries_HYPRESStruct_3d(Mat mat)
490: {
491: PetscErrorCode ierr;
492: Mat_HYPRESStruct *ex = (Mat_HYPRESStruct*) mat->data;
493: PetscInt nvars= ex->nvars;
494: PetscInt size;
495: 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: PetscInt i,*entries,iupper[3],ilower[3];
501: PetscScalar *values;
504: for (i= 0; i< 3; i++) {
505: ilower[i]= ex->hbox.imin[i];
506: iupper[i]= ex->hbox.imax[i];
507: }
509: PetscMalloc2(nvars*7,&entries,nvars*7*size,&values);
510: for (i= 0; i< nvars*7; i++) entries[i]= i;
511: PetscMemzero(values,nvars*7*size*sizeof(PetscScalar));
513: for (i= 0; i< nvars; i++) {
514: PetscStackCallStandard(HYPRE_SStructMatrixSetBoxValues,(ex->ss_mat,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,nvars*7,(HYPRE_Int *)entries,values));
515: }
516: PetscFree2(entries,values);
517: }
518: PetscStackCallStandard(HYPRE_SStructMatrixAssemble,(ex->ss_mat));
519: return(0);
520: }
523: static PetscErrorCode MatSetUp_HYPRESStruct(Mat mat)
524: {
525: PetscErrorCode ierr;
526: Mat_HYPRESStruct *ex = (Mat_HYPRESStruct*) mat->data;
527: PetscInt dim,dof,sw[3],nx,ny,nz;
528: PetscInt ilower[3],iupper[3],ssize,i;
529: DMBoundaryType px,py,pz;
530: DMDAStencilType st;
531: PetscInt nparts= 1; /* assuming only one part */
532: PetscInt part = 0;
533: ISLocalToGlobalMapping ltog;
534: DM da;
537: MatGetDM(mat,(DM*)&da);
538: ex->da = da;
539: PetscObjectReference((PetscObject)da);
541: DMDAGetInfo(ex->da,&dim,0,0,0,0,0,0,&dof,&sw[0],&px,&py,&pz,&st);
542: DMDAGetCorners(ex->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);
543: iupper[0] += ilower[0] - 1;
544: iupper[1] += ilower[1] - 1;
545: iupper[2] += ilower[2] - 1;
546: /* the hypre_Box is used to zero out the matrix entries in MatZeroValues() */
547: ex->hbox.imin[0] = ilower[0];
548: ex->hbox.imin[1] = ilower[1];
549: ex->hbox.imin[2] = ilower[2];
550: ex->hbox.imax[0] = iupper[0];
551: ex->hbox.imax[1] = iupper[1];
552: ex->hbox.imax[2] = iupper[2];
554: ex->dofs_order = 0;
556: /* assuming that the same number of dofs on each gridpoint. Also assume all cell-centred based */
557: ex->nvars= dof;
559: /* create the hypre grid object and set its information */
560: if (px || py || pz) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Ask us to add periodic support by calling HYPRE_SStructGridSetPeriodic()");
561: PetscStackCallStandard(HYPRE_SStructGridCreate,(ex->hcomm,dim,nparts,&ex->ss_grid));
562: PetscStackCallStandard(HYPRE_SStructGridSetExtents,(ex->ss_grid,part,ex->hbox.imin,ex->hbox.imax));
563: {
564: HYPRE_SStructVariable *vartypes;
565: PetscMalloc1(ex->nvars,&vartypes);
566: for (i= 0; i< ex->nvars; i++) vartypes[i]= HYPRE_SSTRUCT_VARIABLE_CELL;
567: PetscStackCallStandard(HYPRE_SStructGridSetVariables,(ex->ss_grid, part, ex->nvars,vartypes));
568: PetscFree(vartypes);
569: }
570: PetscStackCallStandard(HYPRE_SStructGridAssemble,(ex->ss_grid));
572: sw[1] = sw[0];
573: sw[2] = sw[1];
574: /* PetscStackCallStandard(HYPRE_SStructGridSetNumGhost,(ex->ss_grid,sw)); */
576: /* create the hypre stencil object and set its information */
577: if (sw[0] > 1) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Ask us to add support for wider stencils");
578: if (st == DMDA_STENCIL_BOX) SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Ask us to add support for box stencils");
580: if (dim == 1) {
581: PetscInt offsets[3][1] = {{-1},{0},{1}};
582: PetscInt j, cnt;
584: ssize = 3*(ex->nvars);
585: PetscStackCallStandard(HYPRE_SStructStencilCreate,(dim,ssize,&ex->ss_stencil));
586: cnt= 0;
587: for (i = 0; i < (ex->nvars); i++) {
588: for (j = 0; j < 3; j++) {
589: PetscStackCallStandard(HYPRE_SStructStencilSetEntry,(ex->ss_stencil, cnt, (HYPRE_Int *)offsets[j], i));
590: cnt++;
591: }
592: }
593: } else if (dim == 2) {
594: PetscInt offsets[5][2] = {{0,-1},{-1,0},{0,0},{1,0},{0,1}};
595: PetscInt j, cnt;
597: ssize = 5*(ex->nvars);
598: PetscStackCallStandard(HYPRE_SStructStencilCreate,(dim,ssize,&ex->ss_stencil));
599: cnt= 0;
600: for (i= 0; i< (ex->nvars); i++) {
601: for (j= 0; j< 5; j++) {
602: PetscStackCallStandard(HYPRE_SStructStencilSetEntry,(ex->ss_stencil, cnt, (HYPRE_Int *)offsets[j], i));
603: cnt++;
604: }
605: }
606: } else if (dim == 3) {
607: 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}};
608: PetscInt j, cnt;
610: ssize = 7*(ex->nvars);
611: PetscStackCallStandard(HYPRE_SStructStencilCreate,(dim,ssize,&ex->ss_stencil));
612: cnt= 0;
613: for (i= 0; i< (ex->nvars); i++) {
614: for (j= 0; j< 7; j++) {
615: PetscStackCallStandard(HYPRE_SStructStencilSetEntry,(ex->ss_stencil, cnt, (HYPRE_Int *)offsets[j], i));
616: cnt++;
617: }
618: }
619: }
621: /* create the HYPRE graph */
622: PetscStackCallStandard(HYPRE_SStructGraphCreate,(ex->hcomm, ex->ss_grid, &(ex->ss_graph)));
624: /* set the stencil graph. Note that each variable has the same graph. This means that each
625: variable couples to all the other variable and with the same stencil pattern. */
626: for (i= 0; i< (ex->nvars); i++) {
627: PetscStackCallStandard(HYPRE_SStructGraphSetStencil,(ex->ss_graph,part,i,ex->ss_stencil));
628: }
629: PetscStackCallStandard(HYPRE_SStructGraphAssemble,(ex->ss_graph));
631: /* create the HYPRE sstruct vectors for rhs and solution */
632: PetscStackCallStandard(HYPRE_SStructVectorCreate,(ex->hcomm,ex->ss_grid,&ex->ss_b));
633: PetscStackCallStandard(HYPRE_SStructVectorCreate,(ex->hcomm,ex->ss_grid,&ex->ss_x));
634: PetscStackCallStandard(HYPRE_SStructVectorInitialize,(ex->ss_b));
635: PetscStackCallStandard(HYPRE_SStructVectorInitialize,(ex->ss_x));
636: PetscStackCallStandard(HYPRE_SStructVectorAssemble,(ex->ss_b));
637: PetscStackCallStandard(HYPRE_SStructVectorAssemble,(ex->ss_x));
639: /* create the hypre matrix object and set its information */
640: PetscStackCallStandard(HYPRE_SStructMatrixCreate,(ex->hcomm,ex->ss_graph,&ex->ss_mat));
641: PetscStackCallStandard(HYPRE_SStructGridDestroy,(ex->ss_grid));
642: PetscStackCallStandard(HYPRE_SStructStencilDestroy,(ex->ss_stencil));
643: if (ex->needsinitialization) {
644: PetscStackCallStandard(HYPRE_SStructMatrixInitialize,(ex->ss_mat));
645: ex->needsinitialization = PETSC_FALSE;
646: }
648: /* set the global and local sizes of the matrix */
649: DMDAGetCorners(da,0,0,0,&nx,&ny,&nz);
650: MatSetSizes(mat,dof*nx*ny*nz,dof*nx*ny*nz,PETSC_DECIDE,PETSC_DECIDE);
651: PetscLayoutSetBlockSize(mat->rmap,dof);
652: PetscLayoutSetBlockSize(mat->cmap,dof);
653: PetscLayoutSetUp(mat->rmap);
654: PetscLayoutSetUp(mat->cmap);
655: mat->preallocated = PETSC_TRUE;
657: if (dim == 3) {
658: mat->ops->setvalueslocal = MatSetValuesLocal_HYPRESStruct_3d;
659: mat->ops->zerorowslocal = MatZeroRowsLocal_HYPRESStruct_3d;
660: mat->ops->zeroentries = MatZeroEntries_HYPRESStruct_3d;
662: /* MatZeroEntries_HYPRESStruct_3d(mat); */
663: } else SETERRQ(PetscObjectComm((PetscObject)da),PETSC_ERR_SUP,"Only support for 3d DMDA currently");
665: /* get values that will be used repeatedly in MatSetValuesLocal() and MatZeroRowsLocal() repeatedly */
666: MatGetOwnershipRange(mat,&ex->rstart,NULL);
667: DMGetLocalToGlobalMapping(ex->da,<og);
668: ISLocalToGlobalMappingGetIndices(ltog, (const PetscInt **) &ex->gindices);
669: DMDAGetGhostCorners(ex->da,0,0,0,&ex->gnx,&ex->gnxgny,&ex->gnxgnygnz);
671: ex->gnxgny *= ex->gnx;
672: ex->gnxgnygnz *= ex->gnxgny;
674: DMDAGetCorners(ex->da,&ex->xs,&ex->ys,&ex->zs,&ex->nx,&ex->ny,&ex->nz);
676: ex->nxny = ex->nx*ex->ny;
677: ex->nxnynz = ex->nz*ex->nxny;
678: return(0);
679: }
681: PetscErrorCode MatMult_HYPRESStruct(Mat A,Vec x,Vec y)
682: {
683: PetscErrorCode ierr;
684: const PetscScalar *xx;
685: PetscScalar *yy;
686: PetscInt ilower[3],iupper[3];
687: Mat_HYPRESStruct *mx = (Mat_HYPRESStruct*)(A->data);
688: PetscInt ordering= mx->dofs_order;
689: PetscInt nvars = mx->nvars;
690: PetscInt part = 0;
691: PetscInt size;
692: PetscInt i;
695: DMDAGetCorners(mx->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);
696: iupper[0] += ilower[0] - 1;
697: iupper[1] += ilower[1] - 1;
698: iupper[2] += ilower[2] - 1;
700: size= 1;
701: for (i = 0; i < 3; i++) size *= (iupper[i]-ilower[i]+1);
703: /* copy x values over to hypre for variable ordering */
704: if (ordering) {
705: PetscStackCallStandard(HYPRE_SStructVectorSetConstantValues,(mx->ss_b,0.0));
706: VecGetArrayRead(x,&xx);
707: for (i= 0; i< nvars; i++) {
708: PetscStackCallStandard(HYPRE_SStructVectorSetBoxValues,(mx->ss_b,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,(PetscScalar*)xx+(size*i)));
709: }
710: VecRestoreArrayRead(x,&xx);
711: PetscStackCallStandard(HYPRE_SStructVectorAssemble,(mx->ss_b));
712: PetscStackCallStandard(HYPRE_SStructMatrixMatvec,(1.0,mx->ss_mat,mx->ss_b,0.0,mx->ss_x));
714: /* copy solution values back to PETSc */
715: VecGetArray(y,&yy);
716: for (i= 0; i< nvars; i++) {
717: PetscStackCallStandard(HYPRE_SStructVectorGetBoxValues,(mx->ss_x,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,yy+(size*i)));
718: }
719: VecRestoreArray(y,&yy);
720: } else { /* nodal ordering must be mapped to variable ordering for sys_pfmg */
721: PetscScalar *z;
722: PetscInt j, k;
724: PetscMalloc1(nvars*size,&z);
725: PetscStackCallStandard(HYPRE_SStructVectorSetConstantValues,(mx->ss_b,0.0));
726: VecGetArrayRead(x,&xx);
728: /* transform nodal to hypre's variable ordering for sys_pfmg */
729: for (i= 0; i< size; i++) {
730: k= i*nvars;
731: for (j= 0; j< nvars; j++) z[j*size+i]= xx[k+j];
732: }
733: for (i= 0; i< nvars; i++) {
734: PetscStackCallStandard(HYPRE_SStructVectorSetBoxValues,(mx->ss_b,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,z+(size*i)));
735: }
736: VecRestoreArrayRead(x,&xx);
737: PetscStackCallStandard(HYPRE_SStructVectorAssemble,(mx->ss_b));
738: PetscStackCallStandard(HYPRE_SStructMatrixMatvec,(1.0,mx->ss_mat,mx->ss_b,0.0,mx->ss_x));
740: /* copy solution values back to PETSc */
741: VecGetArray(y,&yy);
742: for (i= 0; i< nvars; i++) {
743: PetscStackCallStandard(HYPRE_SStructVectorGetBoxValues,(mx->ss_x,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,z+(size*i)));
744: }
745: /* transform hypre's variable ordering for sys_pfmg to nodal ordering */
746: for (i= 0; i< size; i++) {
747: k= i*nvars;
748: for (j= 0; j< nvars; j++) yy[k+j]= z[j*size+i];
749: }
750: VecRestoreArray(y,&yy);
751: PetscFree(z);
752: }
753: return(0);
754: }
756: PetscErrorCode MatAssemblyEnd_HYPRESStruct(Mat mat,MatAssemblyType mode)
757: {
758: Mat_HYPRESStruct *ex = (Mat_HYPRESStruct*) mat->data;
761: PetscStackCallStandard(HYPRE_SStructMatrixAssemble,(ex->ss_mat));
762: return(0);
763: }
765: PetscErrorCode MatZeroEntries_HYPRESStruct(Mat mat)
766: {
768: /* before the DMDA is set to the matrix the zero doesn't need to do anything */
769: return(0);
770: }
772: PetscErrorCode MatDestroy_HYPRESStruct(Mat mat)
773: {
774: Mat_HYPRESStruct *ex = (Mat_HYPRESStruct*) mat->data;
775: PetscErrorCode ierr;
776: ISLocalToGlobalMapping ltog;
779: DMGetLocalToGlobalMapping(ex->da,<og);
780: ISLocalToGlobalMappingRestoreIndices(ltog, (const PetscInt **) &ex->gindices);
781: PetscStackCallStandard(HYPRE_SStructGraphDestroy,(ex->ss_graph));
782: PetscStackCallStandard(HYPRE_SStructMatrixDestroy,(ex->ss_mat));
783: PetscStackCallStandard(HYPRE_SStructVectorDestroy,(ex->ss_x));
784: PetscStackCallStandard(HYPRE_SStructVectorDestroy,(ex->ss_b));
785: PetscObjectDereference((PetscObject)ex->da);
786: MPI_Comm_free(&(ex->hcomm));
787: PetscFree(ex);
788: return(0);
789: }
791: PETSC_EXTERN PetscErrorCode MatCreate_HYPRESStruct(Mat B)
792: {
793: Mat_HYPRESStruct *ex;
794: PetscErrorCode ierr;
797: PetscNewLog(B,&ex);
798: B->data = (void*)ex;
799: B->rmap->bs = 1;
800: B->assembled = PETSC_FALSE;
802: B->insertmode = NOT_SET_VALUES;
804: B->ops->assemblyend = MatAssemblyEnd_HYPRESStruct;
805: B->ops->mult = MatMult_HYPRESStruct;
806: B->ops->zeroentries = MatZeroEntries_HYPRESStruct;
807: B->ops->destroy = MatDestroy_HYPRESStruct;
808: B->ops->setup = MatSetUp_HYPRESStruct;
810: ex->needsinitialization = PETSC_TRUE;
812: MPI_Comm_dup(PetscObjectComm((PetscObject)B),&(ex->hcomm));
813: PetscObjectChangeTypeName((PetscObject)B,MATHYPRESSTRUCT);
814: return(0);
815: }