Actual source code: fdda.c
1: #include <petsc/private/dmdaimpl.h>
2: #include <petscmat.h>
4: extern PetscErrorCode DMCreateColoring_DA_1d_MPIAIJ(DM, ISColoringType, ISColoring *);
5: extern PetscErrorCode DMCreateColoring_DA_2d_MPIAIJ(DM, ISColoringType, ISColoring *);
6: extern PetscErrorCode DMCreateColoring_DA_2d_5pt_MPIAIJ(DM, ISColoringType, ISColoring *);
7: extern PetscErrorCode DMCreateColoring_DA_3d_MPIAIJ(DM, ISColoringType, ISColoring *);
9: /*
10: For ghost i that may be negative or greater than the upper bound this
11: maps it into the 0:m-1 range using periodicity
12: */
13: #define SetInRange(i, m) (((i) < 0) ? (m) + (i) : (((i) >= (m)) ? (i) - (m) : (i)))
15: static PetscErrorCode DMDASetBlockFills_Private(const PetscInt *dfill, PetscInt w, PetscInt **rfill)
16: {
17: PetscInt i, j, nz, *fill;
19: PetscFunctionBegin;
20: if (!dfill) PetscFunctionReturn(PETSC_SUCCESS);
22: /* count number nonzeros */
23: nz = 0;
24: for (i = 0; i < w; i++) {
25: for (j = 0; j < w; j++) {
26: if (dfill[w * i + j]) nz++;
27: }
28: }
29: PetscCall(PetscMalloc1(nz + w + 1, &fill));
30: /* construct modified CSR storage of nonzero structure */
31: /* fill[0 -- w] marks starts of each row of column indices (and end of last row)
32: so fill[1] - fill[0] gives number of nonzeros in first row etc */
33: nz = w + 1;
34: for (i = 0; i < w; i++) {
35: fill[i] = nz;
36: for (j = 0; j < w; j++) {
37: if (dfill[w * i + j]) {
38: fill[nz] = j;
39: nz++;
40: }
41: }
42: }
43: fill[w] = nz;
45: *rfill = fill;
46: PetscFunctionReturn(PETSC_SUCCESS);
47: }
49: static PetscErrorCode DMDASetBlockFillsSparse_Private(const PetscInt *dfillsparse, PetscInt w, PetscInt **rfill)
50: {
51: PetscInt nz;
53: PetscFunctionBegin;
54: if (!dfillsparse) PetscFunctionReturn(PETSC_SUCCESS);
56: /* Determine number of non-zeros */
57: nz = (dfillsparse[w] - w - 1);
59: /* Allocate space for our copy of the given sparse matrix representation. */
60: PetscCall(PetscMalloc1(nz + w + 1, rfill));
61: PetscCall(PetscArraycpy(*rfill, dfillsparse, nz + w + 1));
62: PetscFunctionReturn(PETSC_SUCCESS);
63: }
65: static PetscErrorCode DMDASetBlockFills_Private2(DM_DA *dd)
66: {
67: PetscInt i, k, cnt = 1;
69: PetscFunctionBegin;
70: /* ofillcount tracks the columns of ofill that have any nonzero in them; the value in each location is the number of
71: columns to the left with any nonzeros in them plus 1 */
72: PetscCall(PetscCalloc1(dd->w, &dd->ofillcols));
73: for (i = 0; i < dd->w; i++) {
74: for (k = dd->ofill[i]; k < dd->ofill[i + 1]; k++) dd->ofillcols[dd->ofill[k]] = 1;
75: }
76: for (i = 0; i < dd->w; i++) {
77: if (dd->ofillcols[i]) dd->ofillcols[i] = cnt++;
78: }
79: PetscFunctionReturn(PETSC_SUCCESS);
80: }
82: /*@
83: DMDASetBlockFills - Sets the fill pattern in each block for a multi-component problem
84: of the matrix returned by `DMCreateMatrix()`.
86: Logically Collective
88: Input Parameters:
89: + da - the `DMDA`
90: . dfill - the fill pattern in the diagonal block (may be `NULL`, means use dense block)
91: - ofill - the fill pattern in the off-diagonal blocks
93: Level: developer
95: Notes:
96: This only makes sense when you are doing multicomponent problems but using the
97: `MATMPIAIJ` matrix format
99: The format for `dfill` and `ofill` is a 2 dimensional dof by dof matrix with 1 entries
100: representing coupling and 0 entries for missing coupling. For example
101: .vb
102: dfill[9] = {1, 0, 0,
103: 1, 1, 0,
104: 0, 1, 1}
105: .ve
106: means that row 0 is coupled with only itself in the diagonal block, row 1 is coupled with
107: itself and row 0 (in the diagonal block) and row 2 is coupled with itself and row 1 (in the
108: diagonal block).
110: `DMDASetGetMatrix()` allows you to provide general code for those more complicated nonzero patterns then
111: can be represented in the `dfill`, `ofill` format
113: Contributed by\:
114: Glenn Hammond
116: .seealso: [](sec_struct), `DM`, `DMDA`, `DMCreateMatrix()`, `DMDASetGetMatrix()`, `DMSetMatrixPreallocateOnly()`, `DMDASetBlockFillsSparse()`
117: @*/
118: PetscErrorCode DMDASetBlockFills(DM da, const PetscInt *dfill, const PetscInt *ofill)
119: {
120: DM_DA *dd = (DM_DA *)da->data;
122: PetscFunctionBegin;
123: /* save the given dfill and ofill information */
124: PetscCall(DMDASetBlockFills_Private(dfill, dd->w, &dd->dfill));
125: PetscCall(DMDASetBlockFills_Private(ofill, dd->w, &dd->ofill));
127: /* count nonzeros in ofill columns */
128: PetscCall(DMDASetBlockFills_Private2(dd));
129: PetscFunctionReturn(PETSC_SUCCESS);
130: }
132: /*@
133: DMDASetBlockFillsSparse - Sets the fill pattern in each block for a multi-component problem
134: of the matrix returned by `DMCreateMatrix()`, using sparse representations
135: of fill patterns.
137: Logically Collective
139: Input Parameters:
140: + da - the `DMDA`
141: . dfillsparse - the sparse fill pattern in the diagonal block (may be `NULL`, means use dense block)
142: - ofillsparse - the sparse fill pattern in the off-diagonal blocks
144: Level: developer
146: Notes:
147: This only makes sense when you are doing multicomponent problems but using the
148: `MATMPIAIJ` matrix format
150: The format for `dfill` and `ofill` is a sparse representation of a
151: dof-by-dof matrix with 1 entries representing coupling and 0 entries
152: for missing coupling. The sparse representation is a 1 dimensional
153: array of length nz + dof + 1, where nz is the number of non-zeros in
154: the matrix. The first dof entries in the array give the
155: starting array indices of each row's items in the rest of the array,
156: the dof+1st item contains the value nz + dof + 1 (i.e. the entire length of the array)
157: and the remaining nz items give the column indices of each of
158: the 1s within the logical 2D matrix. Each row's items within
159: the array are the column indices of the 1s within that row
160: of the 2D matrix. PETSc developers may recognize that this is the
161: same format as that computed by the `DMDASetBlockFills_Private()`
162: function from a dense 2D matrix representation.
164: `DMDASetGetMatrix()` allows you to provide general code for those more complicated nonzero patterns then
165: can be represented in the `dfill`, `ofill` format
167: Contributed by\:
168: Philip C. Roth
170: .seealso: [](sec_struct), `DM`, `DMDA`, `DMDASetBlockFills()`, `DMCreateMatrix()`, `DMDASetGetMatrix()`, `DMSetMatrixPreallocateOnly()`
171: @*/
172: PetscErrorCode DMDASetBlockFillsSparse(DM da, const PetscInt *dfillsparse, const PetscInt *ofillsparse)
173: {
174: DM_DA *dd = (DM_DA *)da->data;
176: PetscFunctionBegin;
177: /* save the given dfill and ofill information */
178: PetscCall(DMDASetBlockFillsSparse_Private(dfillsparse, dd->w, &dd->dfill));
179: PetscCall(DMDASetBlockFillsSparse_Private(ofillsparse, dd->w, &dd->ofill));
181: /* count nonzeros in ofill columns */
182: PetscCall(DMDASetBlockFills_Private2(dd));
183: PetscFunctionReturn(PETSC_SUCCESS);
184: }
186: PetscErrorCode DMCreateColoring_DA(DM da, ISColoringType ctype, ISColoring *coloring)
187: {
188: PetscInt dim, m, n, p, nc;
189: DMBoundaryType bx, by, bz;
190: MPI_Comm comm;
191: PetscMPIInt size;
192: PetscBool isBAIJ;
193: DM_DA *dd = (DM_DA *)da->data;
195: PetscFunctionBegin;
196: /*
197: m
198: ------------------------------------------------------
199: | |
200: | |
201: | ---------------------- |
202: | | | |
203: n | yn | | |
204: | | | |
205: | .--------------------- |
206: | (xs,ys) xn |
207: | . |
208: | (gxs,gys) |
209: | |
210: -----------------------------------------------------
211: */
213: /*
214: nc - number of components per grid point
215: col - number of colors needed in one direction for single component problem
217: */
218: PetscCall(DMDAGetInfo(da, &dim, NULL, NULL, NULL, &m, &n, &p, &nc, NULL, &bx, &by, &bz, NULL));
220: PetscCall(PetscObjectGetComm((PetscObject)da, &comm));
221: PetscCallMPI(MPI_Comm_size(comm, &size));
222: if (ctype == IS_COLORING_LOCAL) {
223: PetscCheck(!((m == 1 && bx == DM_BOUNDARY_PERIODIC) || (n == 1 && by == DM_BOUNDARY_PERIODIC) || (p == 1 && bz == DM_BOUNDARY_PERIODIC)), PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "IS_COLORING_LOCAL cannot be used for periodic boundary condition having both sides of the domain on the same process");
224: }
226: /* Tell the DMDA it has 1 degree of freedom per grid point so that the coloring for BAIJ
227: matrices is for the blocks, not the individual matrix elements */
228: PetscCall(PetscStrbeginswith(da->mattype, MATBAIJ, &isBAIJ));
229: if (!isBAIJ) PetscCall(PetscStrbeginswith(da->mattype, MATMPIBAIJ, &isBAIJ));
230: if (!isBAIJ) PetscCall(PetscStrbeginswith(da->mattype, MATSEQBAIJ, &isBAIJ));
231: if (isBAIJ) {
232: dd->w = 1;
233: dd->xs = dd->xs / nc;
234: dd->xe = dd->xe / nc;
235: dd->Xs = dd->Xs / nc;
236: dd->Xe = dd->Xe / nc;
237: }
239: /*
240: We do not provide a getcoloring function in the DMDA operations because
241: the basic DMDA does not know about matrices. We think of DMDA as being
242: more low-level then matrices.
243: */
244: if (dim == 1) PetscCall(DMCreateColoring_DA_1d_MPIAIJ(da, ctype, coloring));
245: else if (dim == 2) PetscCall(DMCreateColoring_DA_2d_MPIAIJ(da, ctype, coloring));
246: else if (dim == 3) PetscCall(DMCreateColoring_DA_3d_MPIAIJ(da, ctype, coloring));
247: else SETERRQ(PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "Not done for %" PetscInt_FMT " dimension, send us mail petsc-maint@mcs.anl.gov for code", dim);
248: if (isBAIJ) {
249: dd->w = nc;
250: dd->xs = dd->xs * nc;
251: dd->xe = dd->xe * nc;
252: dd->Xs = dd->Xs * nc;
253: dd->Xe = dd->Xe * nc;
254: }
255: PetscFunctionReturn(PETSC_SUCCESS);
256: }
258: PetscErrorCode DMCreateColoring_DA_2d_MPIAIJ(DM da, ISColoringType ctype, ISColoring *coloring)
259: {
260: PetscInt xs, ys, nx, ny, i, j, ii, gxs, gys, gnx, gny, m, n, M, N, dim, s, k, nc, col;
261: PetscInt ncolors = 0;
262: MPI_Comm comm;
263: DMBoundaryType bx, by;
264: DMDAStencilType st;
265: ISColoringValue *colors;
266: DM_DA *dd = (DM_DA *)da->data;
268: PetscFunctionBegin;
269: /*
270: nc - number of components per grid point
271: col - number of colors needed in one direction for single component problem
273: */
274: PetscCall(DMDAGetInfo(da, &dim, &m, &n, NULL, &M, &N, NULL, &nc, &s, &bx, &by, NULL, &st));
275: col = 2 * s + 1;
276: PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &nx, &ny, NULL));
277: PetscCall(DMDAGetGhostCorners(da, &gxs, &gys, NULL, &gnx, &gny, NULL));
278: PetscCall(PetscObjectGetComm((PetscObject)da, &comm));
280: /* special case as taught to us by Paul Hovland */
281: if (st == DMDA_STENCIL_STAR && s == 1) {
282: PetscCall(DMCreateColoring_DA_2d_5pt_MPIAIJ(da, ctype, coloring));
283: } else {
284: if (ctype == IS_COLORING_GLOBAL) {
285: if (!dd->localcoloring) {
286: PetscCall(PetscMalloc1(nc * nx * ny, &colors));
287: ii = 0;
288: for (j = ys; j < ys + ny; j++) {
289: for (i = xs; i < xs + nx; i++) {
290: for (k = 0; k < nc; k++) PetscCall(ISColoringValueCast(k + nc * ((i % col) + col * (j % col)), colors + ii++));
291: }
292: }
293: ncolors = nc + nc * (col - 1 + col * (col - 1));
294: PetscCall(ISColoringCreate(comm, ncolors, nc * nx * ny, colors, PETSC_OWN_POINTER, &dd->localcoloring));
295: }
296: *coloring = dd->localcoloring;
297: } else if (ctype == IS_COLORING_LOCAL) {
298: if (!dd->ghostedcoloring) {
299: PetscCall(PetscMalloc1(nc * gnx * gny, &colors));
300: ii = 0;
301: for (j = gys; j < gys + gny; j++) {
302: for (i = gxs; i < gxs + gnx; i++) {
303: for (k = 0; k < nc; k++) {
304: /* the complicated stuff is to handle periodic boundaries */
305: PetscCall(ISColoringValueCast(k + nc * ((SetInRange(i, m) % col) + col * (SetInRange(j, n) % col)), colors + ii++));
306: }
307: }
308: }
309: ncolors = nc + nc * (col - 1 + col * (col - 1));
310: PetscCall(ISColoringCreate(comm, ncolors, nc * gnx * gny, colors, PETSC_OWN_POINTER, &dd->ghostedcoloring));
311: /* PetscIntView(ncolors,(PetscInt*)colors,0); */
313: PetscCall(ISColoringSetType(dd->ghostedcoloring, IS_COLORING_LOCAL));
314: }
315: *coloring = dd->ghostedcoloring;
316: } else SETERRQ(PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_WRONG, "Unknown ISColoringType %d", (int)ctype);
317: }
318: PetscCall(ISColoringReference(*coloring));
319: PetscFunctionReturn(PETSC_SUCCESS);
320: }
322: PetscErrorCode DMCreateColoring_DA_3d_MPIAIJ(DM da, ISColoringType ctype, ISColoring *coloring)
323: {
324: PetscInt xs, ys, nx, ny, i, j, gxs, gys, gnx, gny, m, n, p, dim, s, k, nc, col, zs, gzs, ii, l, nz, gnz, M, N, P;
325: PetscInt ncolors;
326: MPI_Comm comm;
327: DMBoundaryType bx, by, bz;
328: DMDAStencilType st;
329: ISColoringValue *colors;
330: DM_DA *dd = (DM_DA *)da->data;
332: PetscFunctionBegin;
333: /*
334: nc - number of components per grid point
335: col - number of colors needed in one direction for single component problem
336: */
337: PetscCall(DMDAGetInfo(da, &dim, &m, &n, &p, &M, &N, &P, &nc, &s, &bx, &by, &bz, &st));
338: col = 2 * s + 1;
339: PetscCheck(bx != DM_BOUNDARY_PERIODIC || (m % col) == 0, PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "For coloring efficiency ensure number of grid points in X is divisible by 2*stencil_width + 1");
340: PetscCheck(by != DM_BOUNDARY_PERIODIC || (n % col) == 0, PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "For coloring efficiency ensure number of grid points in Y is divisible by 2*stencil_width + 1");
341: PetscCheck(bz != DM_BOUNDARY_PERIODIC || (p % col) == 0, PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "For coloring efficiency ensure number of grid points in Z is divisible by 2*stencil_width + 1");
343: PetscCall(DMDAGetCorners(da, &xs, &ys, &zs, &nx, &ny, &nz));
344: PetscCall(DMDAGetGhostCorners(da, &gxs, &gys, &gzs, &gnx, &gny, &gnz));
345: PetscCall(PetscObjectGetComm((PetscObject)da, &comm));
347: /* create the coloring */
348: if (ctype == IS_COLORING_GLOBAL) {
349: if (!dd->localcoloring) {
350: PetscCall(PetscMalloc1(nc * nx * ny * nz, &colors));
351: ii = 0;
352: for (k = zs; k < zs + nz; k++) {
353: for (j = ys; j < ys + ny; j++) {
354: for (i = xs; i < xs + nx; i++) {
355: for (l = 0; l < nc; l++) PetscCall(ISColoringValueCast(l + nc * ((i % col) + col * (j % col) + col * col * (k % col)), colors + ii++));
356: }
357: }
358: }
359: ncolors = nc + nc * (col - 1 + col * (col - 1) + col * col * (col - 1));
360: PetscCall(ISColoringCreate(comm, ncolors, nc * nx * ny * nz, colors, PETSC_OWN_POINTER, &dd->localcoloring));
361: }
362: *coloring = dd->localcoloring;
363: } else if (ctype == IS_COLORING_LOCAL) {
364: if (!dd->ghostedcoloring) {
365: PetscCall(PetscMalloc1(nc * gnx * gny * gnz, &colors));
366: ii = 0;
367: for (k = gzs; k < gzs + gnz; k++) {
368: for (j = gys; j < gys + gny; j++) {
369: for (i = gxs; i < gxs + gnx; i++) {
370: for (l = 0; l < nc; l++) {
371: /* the complicated stuff is to handle periodic boundaries */
372: PetscCall(ISColoringValueCast(l + nc * ((SetInRange(i, m) % col) + col * (SetInRange(j, n) % col) + col * col * (SetInRange(k, p) % col)), colors + ii++));
373: }
374: }
375: }
376: }
377: ncolors = nc + nc * (col - 1 + col * (col - 1) + col * col * (col - 1));
378: PetscCall(ISColoringCreate(comm, ncolors, nc * gnx * gny * gnz, colors, PETSC_OWN_POINTER, &dd->ghostedcoloring));
379: PetscCall(ISColoringSetType(dd->ghostedcoloring, IS_COLORING_LOCAL));
380: }
381: *coloring = dd->ghostedcoloring;
382: } else SETERRQ(PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_WRONG, "Unknown ISColoringType %d", (int)ctype);
383: PetscCall(ISColoringReference(*coloring));
384: PetscFunctionReturn(PETSC_SUCCESS);
385: }
387: PetscErrorCode DMCreateColoring_DA_1d_MPIAIJ(DM da, ISColoringType ctype, ISColoring *coloring)
388: {
389: PetscInt xs, nx, i, i1, gxs, gnx, l, m, M, dim, s, nc, col;
390: PetscInt ncolors;
391: MPI_Comm comm;
392: DMBoundaryType bx;
393: ISColoringValue *colors;
394: DM_DA *dd = (DM_DA *)da->data;
396: PetscFunctionBegin;
397: /*
398: nc - number of components per grid point
399: col - number of colors needed in one direction for single component problem
400: */
401: PetscCall(DMDAGetInfo(da, &dim, &m, NULL, NULL, &M, NULL, NULL, &nc, &s, &bx, NULL, NULL, NULL));
402: col = 2 * s + 1;
403: PetscCheck(bx != DM_BOUNDARY_PERIODIC || !(m % col), PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "IS_COLORING_GLOBAL can only be used for periodic boundary conditions if the number of grid points %" PetscInt_FMT " is divisible by the number of colors %" PetscInt_FMT, m, col);
404: PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &nx, NULL, NULL));
405: PetscCall(DMDAGetGhostCorners(da, &gxs, NULL, NULL, &gnx, NULL, NULL));
406: PetscCall(PetscObjectGetComm((PetscObject)da, &comm));
408: /* create the coloring */
409: if (ctype == IS_COLORING_GLOBAL) {
410: if (!dd->localcoloring) {
411: PetscCall(PetscMalloc1(nc * nx, &colors));
412: if (dd->ofillcols) {
413: PetscInt tc = 0;
414: for (i = 0; i < nc; i++) tc += (PetscInt)(dd->ofillcols[i] > 0);
415: i1 = 0;
416: for (i = xs; i < xs + nx; i++) {
417: for (l = 0; l < nc; l++) {
418: if (dd->ofillcols[l] && (i % col)) {
419: PetscCall(ISColoringValueCast(nc - 1 + tc * ((i % col) - 1) + dd->ofillcols[l], colors + i1++));
420: } else {
421: PetscCall(ISColoringValueCast(l, colors + i1++));
422: }
423: }
424: }
425: ncolors = nc + 2 * s * tc;
426: } else {
427: i1 = 0;
428: for (i = xs; i < xs + nx; i++) {
429: for (l = 0; l < nc; l++) PetscCall(ISColoringValueCast(l + nc * (i % col), colors + i1++));
430: }
431: ncolors = nc + nc * (col - 1);
432: }
433: PetscCall(ISColoringCreate(comm, ncolors, nc * nx, colors, PETSC_OWN_POINTER, &dd->localcoloring));
434: }
435: *coloring = dd->localcoloring;
436: } else if (ctype == IS_COLORING_LOCAL) {
437: if (!dd->ghostedcoloring) {
438: PetscCall(PetscMalloc1(nc * gnx, &colors));
439: i1 = 0;
440: for (i = gxs; i < gxs + gnx; i++) {
441: for (l = 0; l < nc; l++) {
442: /* the complicated stuff is to handle periodic boundaries */
443: PetscCall(ISColoringValueCast(l + nc * (SetInRange(i, m) % col), colors + i1++));
444: }
445: }
446: ncolors = nc + nc * (col - 1);
447: PetscCall(ISColoringCreate(comm, ncolors, nc * gnx, colors, PETSC_OWN_POINTER, &dd->ghostedcoloring));
448: PetscCall(ISColoringSetType(dd->ghostedcoloring, IS_COLORING_LOCAL));
449: }
450: *coloring = dd->ghostedcoloring;
451: } else SETERRQ(PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_WRONG, "Unknown ISColoringType %d", (int)ctype);
452: PetscCall(ISColoringReference(*coloring));
453: PetscFunctionReturn(PETSC_SUCCESS);
454: }
456: PetscErrorCode DMCreateColoring_DA_2d_5pt_MPIAIJ(DM da, ISColoringType ctype, ISColoring *coloring)
457: {
458: PetscInt xs, ys, nx, ny, i, j, ii, gxs, gys, gnx, gny, m, n, dim, s, k, nc;
459: PetscInt ncolors;
460: MPI_Comm comm;
461: DMBoundaryType bx, by;
462: ISColoringValue *colors;
463: DM_DA *dd = (DM_DA *)da->data;
465: PetscFunctionBegin;
466: /*
467: nc - number of components per grid point
468: col - number of colors needed in one direction for single component problem
469: */
470: PetscCall(DMDAGetInfo(da, &dim, &m, &n, NULL, NULL, NULL, NULL, &nc, &s, &bx, &by, NULL, NULL));
471: PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &nx, &ny, NULL));
472: PetscCall(DMDAGetGhostCorners(da, &gxs, &gys, NULL, &gnx, &gny, NULL));
473: PetscCall(PetscObjectGetComm((PetscObject)da, &comm));
474: /* create the coloring */
475: if (ctype == IS_COLORING_GLOBAL) {
476: if (!dd->localcoloring) {
477: PetscCall(PetscMalloc1(nc * nx * ny, &colors));
478: ii = 0;
479: for (j = ys; j < ys + ny; j++) {
480: for (i = xs; i < xs + nx; i++) {
481: for (k = 0; k < nc; k++) PetscCall(ISColoringValueCast(k + nc * ((3 * j + i) % 5), colors + ii++));
482: }
483: }
484: ncolors = 5 * nc;
485: PetscCall(ISColoringCreate(comm, ncolors, nc * nx * ny, colors, PETSC_OWN_POINTER, &dd->localcoloring));
486: }
487: *coloring = dd->localcoloring;
488: } else if (ctype == IS_COLORING_LOCAL) {
489: if (!dd->ghostedcoloring) {
490: PetscCall(PetscMalloc1(nc * gnx * gny, &colors));
491: ii = 0;
492: for (j = gys; j < gys + gny; j++) {
493: for (i = gxs; i < gxs + gnx; i++) {
494: for (k = 0; k < nc; k++) PetscCall(ISColoringValueCast(k + nc * ((3 * SetInRange(j, n) + SetInRange(i, m)) % 5), colors + ii++));
495: }
496: }
497: ncolors = 5 * nc;
498: PetscCall(ISColoringCreate(comm, ncolors, nc * gnx * gny, colors, PETSC_OWN_POINTER, &dd->ghostedcoloring));
499: PetscCall(ISColoringSetType(dd->ghostedcoloring, IS_COLORING_LOCAL));
500: }
501: *coloring = dd->ghostedcoloring;
502: } else SETERRQ(PetscObjectComm((PetscObject)da), PETSC_ERR_ARG_WRONG, "Unknown ISColoringType %d", (int)ctype);
503: PetscFunctionReturn(PETSC_SUCCESS);
504: }
506: /* =========================================================================== */
507: extern PetscErrorCode DMCreateMatrix_DA_1d_MPIAIJ(DM, Mat);
508: extern PetscErrorCode DMCreateMatrix_DA_1d_MPIAIJ_Fill(DM, Mat);
509: extern PetscErrorCode DMCreateMatrix_DA_1d_SeqAIJ_NoPreallocation(DM, Mat);
510: extern PetscErrorCode DMCreateMatrix_DA_2d_MPIAIJ(DM, Mat);
511: extern PetscErrorCode DMCreateMatrix_DA_2d_MPIAIJ_Fill(DM, Mat);
512: extern PetscErrorCode DMCreateMatrix_DA_3d_MPIAIJ(DM, Mat);
513: extern PetscErrorCode DMCreateMatrix_DA_3d_MPIAIJ_Fill(DM, Mat);
514: extern PetscErrorCode DMCreateMatrix_DA_2d_MPIBAIJ(DM, Mat);
515: extern PetscErrorCode DMCreateMatrix_DA_3d_MPIBAIJ(DM, Mat);
516: extern PetscErrorCode DMCreateMatrix_DA_2d_MPISBAIJ(DM, Mat);
517: extern PetscErrorCode DMCreateMatrix_DA_3d_MPISBAIJ(DM, Mat);
518: extern PetscErrorCode DMCreateMatrix_DA_2d_MPISELL(DM, Mat);
519: extern PetscErrorCode DMCreateMatrix_DA_3d_MPISELL(DM, Mat);
520: extern PetscErrorCode DMCreateMatrix_DA_IS(DM, Mat);
522: static PetscErrorCode MatView_MPI_DA(Mat A, PetscViewer viewer)
523: {
524: DM da;
525: const char *prefix;
526: Mat AA, Anatural;
527: AO ao;
528: PetscInt rstart, rend, *petsc, i;
529: IS is;
530: MPI_Comm comm;
531: PetscViewerFormat format;
532: PetscBool flag;
534: PetscFunctionBegin;
535: /* Check whether we are just printing info, in which case MatView() already viewed everything we wanted to view */
536: PetscCall(PetscViewerGetFormat(viewer, &format));
537: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) PetscFunctionReturn(PETSC_SUCCESS);
539: PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
540: PetscCall(MatGetDM(A, &da));
541: PetscCheck(da, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Matrix not generated from a DMDA");
543: PetscCall(DMDAGetAO(da, &ao));
544: PetscCall(MatGetOwnershipRange(A, &rstart, &rend));
545: PetscCall(PetscMalloc1(rend - rstart, &petsc));
546: for (i = rstart; i < rend; i++) petsc[i - rstart] = i;
547: PetscCall(AOApplicationToPetsc(ao, rend - rstart, petsc));
548: PetscCall(ISCreateGeneral(comm, rend - rstart, petsc, PETSC_OWN_POINTER, &is));
550: /* call viewer on natural ordering */
551: PetscCall(PetscObjectBaseTypeCompare((PetscObject)A, MATMPISELL, &flag));
552: if (flag) {
553: PetscCall(MatConvert(A, MATAIJ, MAT_INITIAL_MATRIX, &AA));
554: A = AA;
555: }
556: PetscCall(MatCreateSubMatrix(A, is, is, MAT_INITIAL_MATRIX, &Anatural));
557: PetscCall(ISDestroy(&is));
558: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)A, &prefix));
559: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)Anatural, prefix));
560: PetscCall(PetscObjectSetName((PetscObject)Anatural, ((PetscObject)A)->name));
561: ((PetscObject)Anatural)->donotPetscObjectPrintClassNamePrefixType = PETSC_TRUE;
562: PetscCall(MatView(Anatural, viewer));
563: ((PetscObject)Anatural)->donotPetscObjectPrintClassNamePrefixType = PETSC_FALSE;
564: PetscCall(MatDestroy(&Anatural));
565: if (flag) PetscCall(MatDestroy(&AA));
566: PetscFunctionReturn(PETSC_SUCCESS);
567: }
569: static PetscErrorCode MatLoad_MPI_DA(Mat A, PetscViewer viewer)
570: {
571: DM da;
572: Mat Anatural, Aapp;
573: AO ao;
574: PetscInt rstart, rend, *app, i, m, n, M, N;
575: IS is;
576: MPI_Comm comm;
578: PetscFunctionBegin;
579: PetscCall(PetscObjectGetComm((PetscObject)A, &comm));
580: PetscCall(MatGetDM(A, &da));
581: PetscCheck(da, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_WRONG, "Matrix not generated from a DMDA");
583: /* Load the matrix in natural ordering */
584: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), &Anatural));
585: PetscCall(MatSetType(Anatural, ((PetscObject)A)->type_name));
586: PetscCall(MatGetSize(A, &M, &N));
587: PetscCall(MatGetLocalSize(A, &m, &n));
588: PetscCall(MatSetSizes(Anatural, m, n, M, N));
589: PetscCall(MatLoad(Anatural, viewer));
591: /* Map natural ordering to application ordering and create IS */
592: PetscCall(DMDAGetAO(da, &ao));
593: PetscCall(MatGetOwnershipRange(Anatural, &rstart, &rend));
594: PetscCall(PetscMalloc1(rend - rstart, &app));
595: for (i = rstart; i < rend; i++) app[i - rstart] = i;
596: PetscCall(AOPetscToApplication(ao, rend - rstart, app));
597: PetscCall(ISCreateGeneral(comm, rend - rstart, app, PETSC_OWN_POINTER, &is));
599: /* Do permutation and replace header */
600: PetscCall(MatCreateSubMatrix(Anatural, is, is, MAT_INITIAL_MATRIX, &Aapp));
601: PetscCall(MatHeaderReplace(A, &Aapp));
602: PetscCall(ISDestroy(&is));
603: PetscCall(MatDestroy(&Anatural));
604: PetscFunctionReturn(PETSC_SUCCESS);
605: }
607: PetscErrorCode DMCreateMatrix_DA(DM da, Mat *J)
608: {
609: PetscInt dim, dof, nx, ny, nz, dims[3], starts[3], M, N, P;
610: Mat A;
611: MPI_Comm comm;
612: MatType Atype;
613: MatType mtype;
614: PetscMPIInt size;
615: DM_DA *dd = (DM_DA *)da->data;
616: PetscBool aij = PETSC_FALSE, baij = PETSC_FALSE, sbaij = PETSC_FALSE, sell = PETSC_FALSE, is = PETSC_FALSE;
618: PetscFunctionBegin;
619: PetscCall(MatInitializePackage());
620: mtype = da->mattype;
622: /*
623: m
624: ------------------------------------------------------
625: | |
626: | |
627: | ---------------------- |
628: | | | |
629: n | ny | | |
630: | | | |
631: | .--------------------- |
632: | (xs,ys) nx |
633: | . |
634: | (gxs,gys) |
635: | |
636: -----------------------------------------------------
637: */
639: /*
640: nc - number of components per grid point
641: col - number of colors needed in one direction for single component problem
642: */
643: M = dd->M;
644: N = dd->N;
645: P = dd->P;
646: dim = da->dim;
647: dof = dd->w;
648: /* PetscCall(DMDAGetInfo(da,&dim,&M,&N,&P,NULL,NULL,NULL,&dof,NULL,NULL,NULL,NULL,NULL)); */
649: PetscCall(DMDAGetCorners(da, NULL, NULL, NULL, &nx, &ny, &nz));
650: PetscCall(PetscObjectGetComm((PetscObject)da, &comm));
651: PetscCall(MatCreate(comm, &A));
652: PetscCall(MatSetSizes(A, dof * nx * ny * nz, dof * nx * ny * nz, dof * M * N * P, dof * M * N * P));
653: PetscCall(MatSetType(A, mtype));
654: PetscCall(MatSetFromOptions(A));
655: if (dof * nx * ny * nz < da->bind_below) {
656: PetscCall(MatSetBindingPropagates(A, PETSC_TRUE));
657: PetscCall(MatBindToCPU(A, PETSC_TRUE));
658: }
659: PetscCall(MatSetDM(A, da));
660: if (da->structure_only) PetscCall(MatSetOption(A, MAT_STRUCTURE_ONLY, PETSC_TRUE));
661: PetscCall(MatGetType(A, &Atype));
662: /*
663: We do not provide a getmatrix function in the DMDA operations because
664: the basic DMDA does not know about matrices. We think of DMDA as being more
665: more low-level than matrices. This is kind of cheating but, cause sometimes
666: we think of DMDA has higher level than matrices.
668: We could switch based on Atype (or mtype), but we do not since the
669: specialized setting routines depend only on the particular preallocation
670: details of the matrix, not the type itself.
671: */
672: if (!da->prealloc_skip) { // Flag is likely set when user intends to use MatSetPreallocationCOO()
673: PetscCall(PetscObjectHasFunction((PetscObject)A, "MatMPIAIJSetPreallocation_C", &aij));
674: if (!aij) PetscCall(PetscObjectHasFunction((PetscObject)A, "MatSeqAIJSetPreallocation_C", &aij));
675: if (!aij) {
676: PetscCall(PetscObjectHasFunction((PetscObject)A, "MatMPIBAIJSetPreallocation_C", &baij));
677: if (!baij) PetscCall(PetscObjectHasFunction((PetscObject)A, "MatSeqBAIJSetPreallocation_C", &baij));
678: if (!baij) {
679: PetscCall(PetscObjectHasFunction((PetscObject)A, "MatMPISBAIJSetPreallocation_C", &sbaij));
680: if (!sbaij) PetscCall(PetscObjectHasFunction((PetscObject)A, "MatSeqSBAIJSetPreallocation_C", &sbaij));
681: if (!sbaij) {
682: PetscCall(PetscObjectHasFunction((PetscObject)A, "MatMPISELLSetPreallocation_C", &sell));
683: if (!sell) PetscCall(PetscObjectHasFunction((PetscObject)A, "MatSeqSELLSetPreallocation_C", &sell));
684: }
685: if (!sell) PetscCall(PetscObjectHasFunction((PetscObject)A, "MatISSetPreallocation_C", &is));
686: }
687: }
688: }
689: if (aij) {
690: if (dim == 1) {
691: if (dd->ofill) {
692: PetscCall(DMCreateMatrix_DA_1d_MPIAIJ_Fill(da, A));
693: } else {
694: DMBoundaryType bx;
695: PetscMPIInt size;
696: PetscCall(DMDAGetInfo(da, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, &bx, NULL, NULL, NULL));
697: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)da), &size));
698: if (size == 1 && bx == DM_BOUNDARY_NONE) {
699: PetscCall(DMCreateMatrix_DA_1d_SeqAIJ_NoPreallocation(da, A));
700: } else {
701: PetscCall(DMCreateMatrix_DA_1d_MPIAIJ(da, A));
702: }
703: }
704: } else if (dim == 2) {
705: if (dd->ofill) {
706: PetscCall(DMCreateMatrix_DA_2d_MPIAIJ_Fill(da, A));
707: } else {
708: PetscCall(DMCreateMatrix_DA_2d_MPIAIJ(da, A));
709: }
710: } else if (dim == 3) {
711: if (dd->ofill) {
712: PetscCall(DMCreateMatrix_DA_3d_MPIAIJ_Fill(da, A));
713: } else {
714: PetscCall(DMCreateMatrix_DA_3d_MPIAIJ(da, A));
715: }
716: }
717: } else if (baij) {
718: if (dim == 2) {
719: PetscCall(DMCreateMatrix_DA_2d_MPIBAIJ(da, A));
720: } else if (dim == 3) {
721: PetscCall(DMCreateMatrix_DA_3d_MPIBAIJ(da, A));
722: } else SETERRQ(PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "Not implemented for %" PetscInt_FMT " dimension and Matrix Type: %s in %" PetscInt_FMT " dimension! Send mail to petsc-maint@mcs.anl.gov for code", dim, Atype, dim);
723: } else if (sbaij) {
724: if (dim == 2) {
725: PetscCall(DMCreateMatrix_DA_2d_MPISBAIJ(da, A));
726: } else if (dim == 3) {
727: PetscCall(DMCreateMatrix_DA_3d_MPISBAIJ(da, A));
728: } else SETERRQ(PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "Not implemented for %" PetscInt_FMT " dimension and Matrix Type: %s in %" PetscInt_FMT " dimension! Send mail to petsc-maint@mcs.anl.gov for code", dim, Atype, dim);
729: } else if (sell) {
730: if (dim == 2) {
731: PetscCall(DMCreateMatrix_DA_2d_MPISELL(da, A));
732: } else if (dim == 3) {
733: PetscCall(DMCreateMatrix_DA_3d_MPISELL(da, A));
734: } else SETERRQ(PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "Not implemented for %" PetscInt_FMT " dimension and Matrix Type: %s in %" PetscInt_FMT " dimension! Send mail to petsc-maint@mcs.anl.gov for code", dim, Atype, dim);
735: } else if (is) {
736: PetscCall(DMCreateMatrix_DA_IS(da, A));
737: } else { // unknown type or da->prealloc_skip so structural information may be needed, but no prealloc
738: ISLocalToGlobalMapping ltog;
740: PetscCall(MatSetBlockSize(A, dof));
741: PetscCall(MatSetUp(A));
742: PetscCall(DMGetLocalToGlobalMapping(da, <og));
743: PetscCall(MatSetLocalToGlobalMapping(A, ltog, ltog));
744: }
745: PetscCall(DMDAGetGhostCorners(da, &starts[0], &starts[1], &starts[2], &dims[0], &dims[1], &dims[2]));
746: PetscCall(MatSetStencil(A, dim, dims, starts, dof));
747: PetscCall(MatSetDM(A, da));
748: PetscCallMPI(MPI_Comm_size(comm, &size));
749: if (size > 1) {
750: /* change viewer to display matrix in natural ordering */
751: PetscCall(MatSetOperation(A, MATOP_VIEW, (PetscErrorCodeFn *)MatView_MPI_DA));
752: PetscCall(MatSetOperation(A, MATOP_LOAD, (PetscErrorCodeFn *)MatLoad_MPI_DA));
753: }
754: *J = A;
755: PetscFunctionReturn(PETSC_SUCCESS);
756: }
758: PetscErrorCode DMCreateMatrix_DA_IS(DM dm, Mat J)
759: {
760: DM_DA *da = (DM_DA *)dm->data;
761: Mat lJ, P;
762: ISLocalToGlobalMapping ltog;
763: IS is;
764: PetscBT bt;
765: const PetscInt *e_loc, *idx;
766: PetscInt i, nel, nen, nv, dof, *gidx, n, N;
768: /* The l2g map of DMDA has all ghosted nodes, and e_loc is a subset of all the local nodes (including the ghosted)
769: We need to filter out the local indices that are not represented through the DMDAGetElements decomposition */
770: PetscFunctionBegin;
771: dof = da->w;
772: PetscCall(MatSetBlockSize(J, dof));
773: PetscCall(DMGetLocalToGlobalMapping(dm, <og));
775: /* flag local elements indices in local DMDA numbering */
776: PetscCall(ISLocalToGlobalMappingGetSize(ltog, &nv));
777: PetscCall(PetscBTCreate(nv / dof, &bt));
778: PetscCall(DMDAGetElements(dm, &nel, &nen, &e_loc)); /* this will throw an error if the stencil type is not DMDA_STENCIL_BOX */
779: for (i = 0; i < nel * nen; i++) PetscCall(PetscBTSet(bt, e_loc[i]));
781: /* filter out (set to -1) the global indices not used by the local elements */
782: PetscCall(PetscMalloc1(nv / dof, &gidx));
783: PetscCall(ISLocalToGlobalMappingGetBlockIndices(ltog, &idx));
784: PetscCall(PetscArraycpy(gidx, idx, nv / dof));
785: PetscCall(ISLocalToGlobalMappingRestoreBlockIndices(ltog, &idx));
786: for (i = 0; i < nv / dof; i++)
787: if (!PetscBTLookup(bt, i)) gidx[i] = -1;
788: PetscCall(PetscBTDestroy(&bt));
789: PetscCall(ISCreateBlock(PetscObjectComm((PetscObject)dm), dof, nv / dof, gidx, PETSC_OWN_POINTER, &is));
790: PetscCall(ISLocalToGlobalMappingCreateIS(is, <og));
791: PetscCall(MatSetLocalToGlobalMapping(J, ltog, ltog));
792: PetscCall(ISLocalToGlobalMappingDestroy(<og));
793: PetscCall(ISDestroy(&is));
795: /* Preallocation */
796: if (dm->prealloc_skip) {
797: PetscCall(MatSetUp(J));
798: } else {
799: PetscCall(MatISGetLocalMat(J, &lJ));
800: PetscCall(MatGetLocalToGlobalMapping(lJ, <og, NULL));
801: PetscCall(MatCreate(PetscObjectComm((PetscObject)lJ), &P));
802: PetscCall(MatSetType(P, MATPREALLOCATOR));
803: PetscCall(MatSetLocalToGlobalMapping(P, ltog, ltog));
804: PetscCall(MatGetSize(lJ, &N, NULL));
805: PetscCall(MatGetLocalSize(lJ, &n, NULL));
806: PetscCall(MatSetSizes(P, n, n, N, N));
807: PetscCall(MatSetBlockSize(P, dof));
808: PetscCall(MatSetUp(P));
809: for (i = 0; i < nel; i++) PetscCall(MatSetValuesBlockedLocal(P, nen, e_loc + i * nen, nen, e_loc + i * nen, NULL, INSERT_VALUES));
810: PetscCall(MatPreallocatorPreallocate(P, (PetscBool)!da->prealloc_only, lJ));
811: PetscCall(MatISRestoreLocalMat(J, &lJ));
812: PetscCall(DMDARestoreElements(dm, &nel, &nen, &e_loc));
813: PetscCall(MatDestroy(&P));
815: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
816: PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
817: }
818: PetscFunctionReturn(PETSC_SUCCESS);
819: }
821: PetscErrorCode DMCreateMatrix_DA_2d_MPISELL(DM da, Mat J)
822: {
823: PetscInt xs, ys, nx, ny, i, j, slot, gxs, gys, gnx, gny, m, n, dim, s, *cols = NULL, k, nc, *rows = NULL, col, cnt, l, p;
824: PetscInt lstart, lend, pstart, pend, *dnz, *onz;
825: MPI_Comm comm;
826: PetscScalar *values;
827: DMBoundaryType bx, by;
828: ISLocalToGlobalMapping ltog;
829: DMDAStencilType st;
831: PetscFunctionBegin;
832: /*
833: nc - number of components per grid point
834: col - number of colors needed in one direction for single component problem
835: */
836: PetscCall(DMDAGetInfo(da, &dim, &m, &n, NULL, NULL, NULL, NULL, &nc, &s, &bx, &by, NULL, &st));
837: col = 2 * s + 1;
838: PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &nx, &ny, NULL));
839: PetscCall(DMDAGetGhostCorners(da, &gxs, &gys, NULL, &gnx, &gny, NULL));
840: PetscCall(PetscObjectGetComm((PetscObject)da, &comm));
842: PetscCall(PetscMalloc2(nc, &rows, col * col * nc * nc, &cols));
843: PetscCall(DMGetLocalToGlobalMapping(da, <og));
845: PetscCall(MatSetBlockSize(J, nc));
846: /* determine the matrix preallocation information */
847: MatPreallocateBegin(comm, nc * nx * ny, nc * nx * ny, dnz, onz);
848: for (i = xs; i < xs + nx; i++) {
849: pstart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s, -i));
850: pend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s, m - i - 1));
852: for (j = ys; j < ys + ny; j++) {
853: slot = i - gxs + gnx * (j - gys);
855: lstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s, -j));
856: lend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s, n - j - 1));
858: cnt = 0;
859: for (k = 0; k < nc; k++) {
860: for (l = lstart; l < lend + 1; l++) {
861: for (p = pstart; p < pend + 1; p++) {
862: if (st == DMDA_STENCIL_BOX || !l || !p) { /* entries on star have either l = 0 or p = 0 */
863: cols[cnt++] = k + nc * (slot + gnx * l + p);
864: }
865: }
866: }
867: rows[k] = k + nc * slot;
868: }
869: PetscCall(MatPreallocateSetLocal(ltog, nc, rows, ltog, cnt, cols, dnz, onz));
870: }
871: }
872: PetscCall(MatSetBlockSize(J, nc));
873: PetscCall(MatSeqSELLSetPreallocation(J, 0, dnz));
874: PetscCall(MatMPISELLSetPreallocation(J, 0, dnz, 0, onz));
875: MatPreallocateEnd(dnz, onz);
877: PetscCall(MatSetLocalToGlobalMapping(J, ltog, ltog));
879: /*
880: For each node in the grid: we get the neighbors in the local (on processor ordering
881: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
882: PETSc ordering.
883: */
884: if (!da->prealloc_only) {
885: PetscCall(PetscCalloc1(col * col * nc * nc, &values));
886: for (i = xs; i < xs + nx; i++) {
887: pstart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s, -i));
888: pend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s, m - i - 1));
890: for (j = ys; j < ys + ny; j++) {
891: slot = i - gxs + gnx * (j - gys);
893: lstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s, -j));
894: lend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s, n - j - 1));
896: cnt = 0;
897: for (k = 0; k < nc; k++) {
898: for (l = lstart; l < lend + 1; l++) {
899: for (p = pstart; p < pend + 1; p++) {
900: if (st == DMDA_STENCIL_BOX || !l || !p) { /* entries on star have either l = 0 or p = 0 */
901: cols[cnt++] = k + nc * (slot + gnx * l + p);
902: }
903: }
904: }
905: rows[k] = k + nc * slot;
906: }
907: PetscCall(MatSetValuesLocal(J, nc, rows, cnt, cols, values, INSERT_VALUES));
908: }
909: }
910: PetscCall(PetscFree(values));
911: /* do not copy values to GPU since they are all zero and not yet needed there */
912: PetscCall(MatBindToCPU(J, PETSC_TRUE));
913: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
914: PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
915: PetscCall(MatBindToCPU(J, PETSC_FALSE));
916: PetscCall(MatSetOption(J, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE));
917: }
918: PetscCall(PetscFree2(rows, cols));
919: PetscFunctionReturn(PETSC_SUCCESS);
920: }
922: PetscErrorCode DMCreateMatrix_DA_3d_MPISELL(DM da, Mat J)
923: {
924: PetscInt xs, ys, nx, ny, i, j, slot, gxs, gys, gnx, gny;
925: PetscInt m, n, dim, s, *cols = NULL, k, nc, *rows = NULL, col, cnt, l, p, *dnz = NULL, *onz = NULL;
926: PetscInt istart, iend, jstart, jend, kstart, kend, zs, nz, gzs, gnz, ii, jj, kk, M, N, P;
927: MPI_Comm comm;
928: PetscScalar *values;
929: DMBoundaryType bx, by, bz;
930: ISLocalToGlobalMapping ltog;
931: DMDAStencilType st;
933: PetscFunctionBegin;
934: /*
935: nc - number of components per grid point
936: col - number of colors needed in one direction for single component problem
937: */
938: PetscCall(DMDAGetInfo(da, &dim, &m, &n, &p, &M, &N, &P, &nc, &s, &bx, &by, &bz, &st));
939: col = 2 * s + 1;
940: PetscCall(DMDAGetCorners(da, &xs, &ys, &zs, &nx, &ny, &nz));
941: PetscCall(DMDAGetGhostCorners(da, &gxs, &gys, &gzs, &gnx, &gny, &gnz));
942: PetscCall(PetscObjectGetComm((PetscObject)da, &comm));
944: PetscCall(PetscMalloc2(nc, &rows, col * col * col * nc * nc, &cols));
945: PetscCall(DMGetLocalToGlobalMapping(da, <og));
947: PetscCall(MatSetBlockSize(J, nc));
948: /* determine the matrix preallocation information */
949: MatPreallocateBegin(comm, nc * nx * ny * nz, nc * nx * ny * nz, dnz, onz);
950: for (i = xs; i < xs + nx; i++) {
951: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s, -i));
952: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s, m - i - 1));
953: for (j = ys; j < ys + ny; j++) {
954: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s, -j));
955: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s, n - j - 1));
956: for (k = zs; k < zs + nz; k++) {
957: kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s, -k));
958: kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s, p - k - 1));
960: slot = i - gxs + gnx * (j - gys) + gnx * gny * (k - gzs);
962: cnt = 0;
963: for (l = 0; l < nc; l++) {
964: for (ii = istart; ii < iend + 1; ii++) {
965: for (jj = jstart; jj < jend + 1; jj++) {
966: for (kk = kstart; kk < kend + 1; kk++) {
967: if (st == DMDA_STENCIL_BOX || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) { /* entries on star*/
968: cols[cnt++] = l + nc * (slot + ii + gnx * jj + gnx * gny * kk);
969: }
970: }
971: }
972: }
973: rows[l] = l + nc * slot;
974: }
975: PetscCall(MatPreallocateSetLocal(ltog, nc, rows, ltog, cnt, cols, dnz, onz));
976: }
977: }
978: }
979: PetscCall(MatSetBlockSize(J, nc));
980: PetscCall(MatSeqSELLSetPreallocation(J, 0, dnz));
981: PetscCall(MatMPISELLSetPreallocation(J, 0, dnz, 0, onz));
982: MatPreallocateEnd(dnz, onz);
983: PetscCall(MatSetLocalToGlobalMapping(J, ltog, ltog));
985: /*
986: For each node in the grid: we get the neighbors in the local (on processor ordering
987: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
988: PETSc ordering.
989: */
990: if (!da->prealloc_only) {
991: PetscCall(PetscCalloc1(col * col * col * nc * nc * nc, &values));
992: for (i = xs; i < xs + nx; i++) {
993: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s, -i));
994: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s, m - i - 1));
995: for (j = ys; j < ys + ny; j++) {
996: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s, -j));
997: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s, n - j - 1));
998: for (k = zs; k < zs + nz; k++) {
999: kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s, -k));
1000: kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s, p - k - 1));
1002: slot = i - gxs + gnx * (j - gys) + gnx * gny * (k - gzs);
1004: cnt = 0;
1005: for (l = 0; l < nc; l++) {
1006: for (ii = istart; ii < iend + 1; ii++) {
1007: for (jj = jstart; jj < jend + 1; jj++) {
1008: for (kk = kstart; kk < kend + 1; kk++) {
1009: if (st == DMDA_STENCIL_BOX || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) { /* entries on star*/
1010: cols[cnt++] = l + nc * (slot + ii + gnx * jj + gnx * gny * kk);
1011: }
1012: }
1013: }
1014: }
1015: rows[l] = l + nc * slot;
1016: }
1017: PetscCall(MatSetValuesLocal(J, nc, rows, cnt, cols, values, INSERT_VALUES));
1018: }
1019: }
1020: }
1021: PetscCall(PetscFree(values));
1022: /* do not copy values to GPU since they are all zero and not yet needed there */
1023: PetscCall(MatBindToCPU(J, PETSC_TRUE));
1024: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
1025: PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
1026: PetscCall(MatBindToCPU(J, PETSC_FALSE));
1027: PetscCall(MatSetOption(J, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE));
1028: }
1029: PetscCall(PetscFree2(rows, cols));
1030: PetscFunctionReturn(PETSC_SUCCESS);
1031: }
1033: PetscErrorCode DMCreateMatrix_DA_2d_MPIAIJ(DM da, Mat J)
1034: {
1035: PetscInt xs, ys, nx, ny, i, j, slot, gxs, gys, gnx, gny, m, n, dim, s, *cols = NULL, k, nc, *rows = NULL, col, cnt, l, p, M, N;
1036: PetscInt lstart, lend, pstart, pend, *dnz, *onz;
1037: MPI_Comm comm;
1038: DMBoundaryType bx, by;
1039: ISLocalToGlobalMapping ltog, mltog;
1040: DMDAStencilType st;
1041: PetscBool removedups = PETSC_FALSE, alreadyboundtocpu = PETSC_TRUE;
1043: PetscFunctionBegin;
1044: /*
1045: nc - number of components per grid point
1046: col - number of colors needed in one direction for single component problem
1047: */
1048: PetscCall(DMDAGetInfo(da, &dim, &m, &n, NULL, &M, &N, NULL, &nc, &s, &bx, &by, NULL, &st));
1049: if (bx == DM_BOUNDARY_NONE && by == DM_BOUNDARY_NONE) PetscCall(MatSetOption(J, MAT_SORTED_FULL, PETSC_TRUE));
1050: col = 2 * s + 1;
1051: /*
1052: With one processor in periodic domains in a skinny dimension the code will label nonzero columns multiple times
1053: because of "wrapping" around the end of the domain hitting an entry already counted in the other direction.
1054: */
1055: if (M == 1 && 2 * s >= m) removedups = PETSC_TRUE;
1056: if (N == 1 && 2 * s >= n) removedups = PETSC_TRUE;
1057: PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &nx, &ny, NULL));
1058: PetscCall(DMDAGetGhostCorners(da, &gxs, &gys, NULL, &gnx, &gny, NULL));
1059: PetscCall(PetscObjectGetComm((PetscObject)da, &comm));
1061: PetscCall(PetscMalloc2(nc, &rows, col * col * nc * nc, &cols));
1062: PetscCall(DMGetLocalToGlobalMapping(da, <og));
1064: PetscCall(MatSetBlockSize(J, nc));
1065: /* determine the matrix preallocation information */
1066: MatPreallocateBegin(comm, nc * nx * ny, nc * nx * ny, dnz, onz);
1067: for (i = xs; i < xs + nx; i++) {
1068: pstart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s, -i));
1069: pend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s, m - i - 1));
1071: for (j = ys; j < ys + ny; j++) {
1072: slot = i - gxs + gnx * (j - gys);
1074: lstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s, -j));
1075: lend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s, n - j - 1));
1077: cnt = 0;
1078: for (k = 0; k < nc; k++) {
1079: for (l = lstart; l < lend + 1; l++) {
1080: for (p = pstart; p < pend + 1; p++) {
1081: if (st == DMDA_STENCIL_BOX || !l || !p) { /* entries on star have either l = 0 or p = 0 */
1082: cols[cnt++] = k + nc * (slot + gnx * l + p);
1083: }
1084: }
1085: }
1086: rows[k] = k + nc * slot;
1087: }
1088: if (removedups) PetscCall(MatPreallocateSetLocalRemoveDups(ltog, nc, rows, ltog, cnt, cols, dnz, onz));
1089: else PetscCall(MatPreallocateSetLocal(ltog, nc, rows, ltog, cnt, cols, dnz, onz));
1090: }
1091: }
1092: PetscCall(MatSetBlockSize(J, nc));
1093: PetscCall(MatSeqAIJSetPreallocation(J, 0, dnz));
1094: PetscCall(MatMPIAIJSetPreallocation(J, 0, dnz, 0, onz));
1095: MatPreallocateEnd(dnz, onz);
1096: PetscCall(MatGetLocalToGlobalMapping(J, &mltog, NULL));
1097: if (!mltog) PetscCall(MatSetLocalToGlobalMapping(J, ltog, ltog));
1099: /*
1100: For each node in the grid: we get the neighbors in the local (on processor ordering
1101: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1102: PETSc ordering.
1103: */
1104: if (!da->prealloc_only) {
1105: for (i = xs; i < xs + nx; i++) {
1106: pstart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s, -i));
1107: pend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s, m - i - 1));
1109: for (j = ys; j < ys + ny; j++) {
1110: slot = i - gxs + gnx * (j - gys);
1112: lstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s, -j));
1113: lend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s, n - j - 1));
1115: cnt = 0;
1116: for (l = lstart; l < lend + 1; l++) {
1117: for (p = pstart; p < pend + 1; p++) {
1118: if (st == DMDA_STENCIL_BOX || !l || !p) { /* entries on star have either l = 0 or p = 0 */
1119: cols[cnt++] = nc * (slot + gnx * l + p);
1120: for (k = 1; k < nc; k++) {
1121: cols[cnt] = 1 + cols[cnt - 1];
1122: cnt++;
1123: }
1124: }
1125: }
1126: }
1127: for (k = 0; k < nc; k++) rows[k] = k + nc * slot;
1128: PetscCall(MatSetValuesLocal(J, nc, rows, cnt, cols, NULL, INSERT_VALUES));
1129: }
1130: }
1131: /* do not copy values to GPU since they are all zero and not yet needed there */
1132: PetscCall(MatBoundToCPU(J, &alreadyboundtocpu));
1133: PetscCall(MatBindToCPU(J, PETSC_TRUE));
1134: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
1135: PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
1136: if (!alreadyboundtocpu) PetscCall(MatBindToCPU(J, PETSC_FALSE));
1137: PetscCall(MatSetOption(J, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE));
1138: if (bx == DM_BOUNDARY_NONE && by == DM_BOUNDARY_NONE) PetscCall(MatSetOption(J, MAT_SORTED_FULL, PETSC_FALSE));
1139: }
1140: PetscCall(PetscFree2(rows, cols));
1141: PetscFunctionReturn(PETSC_SUCCESS);
1142: }
1144: PetscErrorCode DMCreateMatrix_DA_2d_MPIAIJ_Fill(DM da, Mat J)
1145: {
1146: PetscInt xs, ys, nx, ny, i, j, slot, gxs, gys, gnx, gny;
1147: PetscInt m, n, dim, s, *cols, k, nc, row, col, cnt, maxcnt = 0, l, p, M, N;
1148: PetscInt lstart, lend, pstart, pend, *dnz, *onz;
1149: DM_DA *dd = (DM_DA *)da->data;
1150: PetscInt ifill_col, *ofill = dd->ofill, *dfill = dd->dfill;
1151: MPI_Comm comm;
1152: DMBoundaryType bx, by;
1153: ISLocalToGlobalMapping ltog;
1154: DMDAStencilType st;
1155: PetscBool removedups = PETSC_FALSE;
1157: PetscFunctionBegin;
1158: /*
1159: nc - number of components per grid point
1160: col - number of colors needed in one direction for single component problem
1161: */
1162: PetscCall(DMDAGetInfo(da, &dim, &m, &n, NULL, &M, &N, NULL, &nc, &s, &bx, &by, NULL, &st));
1163: col = 2 * s + 1;
1164: /*
1165: With one processor in periodic domains in a skinny dimension the code will label nonzero columns multiple times
1166: because of "wrapping" around the end of the domain hitting an entry already counted in the other direction.
1167: */
1168: if (M == 1 && 2 * s >= m) removedups = PETSC_TRUE;
1169: if (N == 1 && 2 * s >= n) removedups = PETSC_TRUE;
1170: PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &nx, &ny, NULL));
1171: PetscCall(DMDAGetGhostCorners(da, &gxs, &gys, NULL, &gnx, &gny, NULL));
1172: PetscCall(PetscObjectGetComm((PetscObject)da, &comm));
1174: PetscCall(PetscMalloc1(col * col * nc, &cols));
1175: PetscCall(DMGetLocalToGlobalMapping(da, <og));
1177: PetscCall(MatSetBlockSize(J, nc));
1178: /* determine the matrix preallocation information */
1179: MatPreallocateBegin(comm, nc * nx * ny, nc * nx * ny, dnz, onz);
1180: for (i = xs; i < xs + nx; i++) {
1181: pstart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s, -i));
1182: pend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s, m - i - 1));
1184: for (j = ys; j < ys + ny; j++) {
1185: slot = i - gxs + gnx * (j - gys);
1187: lstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s, -j));
1188: lend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s, n - j - 1));
1190: for (k = 0; k < nc; k++) {
1191: cnt = 0;
1192: for (l = lstart; l < lend + 1; l++) {
1193: for (p = pstart; p < pend + 1; p++) {
1194: if (l || p) {
1195: if (st == DMDA_STENCIL_BOX || !l || !p) { /* entries on star */
1196: for (ifill_col = ofill[k]; ifill_col < ofill[k + 1]; ifill_col++) cols[cnt++] = ofill[ifill_col] + nc * (slot + gnx * l + p);
1197: }
1198: } else {
1199: if (dfill) {
1200: for (ifill_col = dfill[k]; ifill_col < dfill[k + 1]; ifill_col++) cols[cnt++] = dfill[ifill_col] + nc * (slot + gnx * l + p);
1201: } else {
1202: for (ifill_col = 0; ifill_col < nc; ifill_col++) cols[cnt++] = ifill_col + nc * (slot + gnx * l + p);
1203: }
1204: }
1205: }
1206: }
1207: row = k + nc * slot;
1208: maxcnt = PetscMax(maxcnt, cnt);
1209: if (removedups) PetscCall(MatPreallocateSetLocalRemoveDups(ltog, 1, &row, ltog, cnt, cols, dnz, onz));
1210: else PetscCall(MatPreallocateSetLocal(ltog, 1, &row, ltog, cnt, cols, dnz, onz));
1211: }
1212: }
1213: }
1214: PetscCall(MatSeqAIJSetPreallocation(J, 0, dnz));
1215: PetscCall(MatMPIAIJSetPreallocation(J, 0, dnz, 0, onz));
1216: MatPreallocateEnd(dnz, onz);
1217: PetscCall(MatSetLocalToGlobalMapping(J, ltog, ltog));
1219: /*
1220: For each node in the grid: we get the neighbors in the local (on processor ordering
1221: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1222: PETSc ordering.
1223: */
1224: if (!da->prealloc_only) {
1225: for (i = xs; i < xs + nx; i++) {
1226: pstart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s, -i));
1227: pend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s, m - i - 1));
1229: for (j = ys; j < ys + ny; j++) {
1230: slot = i - gxs + gnx * (j - gys);
1232: lstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s, -j));
1233: lend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s, n - j - 1));
1235: for (k = 0; k < nc; k++) {
1236: cnt = 0;
1237: for (l = lstart; l < lend + 1; l++) {
1238: for (p = pstart; p < pend + 1; p++) {
1239: if (l || p) {
1240: if (st == DMDA_STENCIL_BOX || !l || !p) { /* entries on star */
1241: for (ifill_col = ofill[k]; ifill_col < ofill[k + 1]; ifill_col++) cols[cnt++] = ofill[ifill_col] + nc * (slot + gnx * l + p);
1242: }
1243: } else {
1244: if (dfill) {
1245: for (ifill_col = dfill[k]; ifill_col < dfill[k + 1]; ifill_col++) cols[cnt++] = dfill[ifill_col] + nc * (slot + gnx * l + p);
1246: } else {
1247: for (ifill_col = 0; ifill_col < nc; ifill_col++) cols[cnt++] = ifill_col + nc * (slot + gnx * l + p);
1248: }
1249: }
1250: }
1251: }
1252: row = k + nc * slot;
1253: PetscCall(MatSetValuesLocal(J, 1, &row, cnt, cols, NULL, INSERT_VALUES));
1254: }
1255: }
1256: }
1257: /* do not copy values to GPU since they are all zero and not yet needed there */
1258: PetscCall(MatBindToCPU(J, PETSC_TRUE));
1259: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
1260: PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
1261: PetscCall(MatBindToCPU(J, PETSC_FALSE));
1262: PetscCall(MatSetOption(J, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE));
1263: }
1264: PetscCall(PetscFree(cols));
1265: PetscFunctionReturn(PETSC_SUCCESS);
1266: }
1268: PetscErrorCode DMCreateMatrix_DA_3d_MPIAIJ(DM da, Mat J)
1269: {
1270: PetscInt xs, ys, nx, ny, i, j, slot, gxs, gys, gnx, gny;
1271: PetscInt m, n, dim, s, *cols = NULL, k, nc, *rows = NULL, col, cnt, l, p, *dnz = NULL, *onz = NULL;
1272: PetscInt istart, iend, jstart, jend, kstart, kend, zs, nz, gzs, gnz, ii, jj, kk, M, N, P;
1273: MPI_Comm comm;
1274: DMBoundaryType bx, by, bz;
1275: ISLocalToGlobalMapping ltog, mltog;
1276: DMDAStencilType st;
1277: PetscBool removedups = PETSC_FALSE;
1279: PetscFunctionBegin;
1280: /*
1281: nc - number of components per grid point
1282: col - number of colors needed in one direction for single component problem
1283: */
1284: PetscCall(DMDAGetInfo(da, &dim, &m, &n, &p, &M, &N, &P, &nc, &s, &bx, &by, &bz, &st));
1285: if (bx == DM_BOUNDARY_NONE && by == DM_BOUNDARY_NONE && bz == DM_BOUNDARY_NONE) PetscCall(MatSetOption(J, MAT_SORTED_FULL, PETSC_TRUE));
1286: col = 2 * s + 1;
1288: /*
1289: With one processor in periodic domains in a skinny dimension the code will label nonzero columns multiple times
1290: because of "wrapping" around the end of the domain hitting an entry already counted in the other direction.
1291: */
1292: if (M == 1 && 2 * s >= m) removedups = PETSC_TRUE;
1293: if (N == 1 && 2 * s >= n) removedups = PETSC_TRUE;
1294: if (P == 1 && 2 * s >= p) removedups = PETSC_TRUE;
1296: PetscCall(DMDAGetCorners(da, &xs, &ys, &zs, &nx, &ny, &nz));
1297: PetscCall(DMDAGetGhostCorners(da, &gxs, &gys, &gzs, &gnx, &gny, &gnz));
1298: PetscCall(PetscObjectGetComm((PetscObject)da, &comm));
1300: PetscCall(PetscMalloc2(nc, &rows, col * col * col * nc * nc, &cols));
1301: PetscCall(DMGetLocalToGlobalMapping(da, <og));
1303: PetscCall(MatSetBlockSize(J, nc));
1304: /* determine the matrix preallocation information */
1305: MatPreallocateBegin(comm, nc * nx * ny * nz, nc * nx * ny * nz, dnz, onz);
1306: for (i = xs; i < xs + nx; i++) {
1307: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s, -i));
1308: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s, m - i - 1));
1309: for (j = ys; j < ys + ny; j++) {
1310: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s, -j));
1311: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s, n - j - 1));
1312: for (k = zs; k < zs + nz; k++) {
1313: kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s, -k));
1314: kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s, p - k - 1));
1316: slot = i - gxs + gnx * (j - gys) + gnx * gny * (k - gzs);
1318: cnt = 0;
1319: for (l = 0; l < nc; l++) {
1320: for (ii = istart; ii < iend + 1; ii++) {
1321: for (jj = jstart; jj < jend + 1; jj++) {
1322: for (kk = kstart; kk < kend + 1; kk++) {
1323: if (st == DMDA_STENCIL_BOX || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) { /* entries on star*/
1324: cols[cnt++] = l + nc * (slot + ii + gnx * jj + gnx * gny * kk);
1325: }
1326: }
1327: }
1328: }
1329: rows[l] = l + nc * slot;
1330: }
1331: if (removedups) PetscCall(MatPreallocateSetLocalRemoveDups(ltog, nc, rows, ltog, cnt, cols, dnz, onz));
1332: else PetscCall(MatPreallocateSetLocal(ltog, nc, rows, ltog, cnt, cols, dnz, onz));
1333: }
1334: }
1335: }
1336: PetscCall(MatSetBlockSize(J, nc));
1337: PetscCall(MatSeqAIJSetPreallocation(J, 0, dnz));
1338: PetscCall(MatMPIAIJSetPreallocation(J, 0, dnz, 0, onz));
1339: MatPreallocateEnd(dnz, onz);
1340: PetscCall(MatGetLocalToGlobalMapping(J, &mltog, NULL));
1341: if (!mltog) PetscCall(MatSetLocalToGlobalMapping(J, ltog, ltog));
1343: /*
1344: For each node in the grid: we get the neighbors in the local (on processor ordering
1345: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1346: PETSc ordering.
1347: */
1348: if (!da->prealloc_only) {
1349: for (i = xs; i < xs + nx; i++) {
1350: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s, -i));
1351: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s, m - i - 1));
1352: for (j = ys; j < ys + ny; j++) {
1353: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s, -j));
1354: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s, n - j - 1));
1355: for (k = zs; k < zs + nz; k++) {
1356: kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s, -k));
1357: kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s, p - k - 1));
1359: slot = i - gxs + gnx * (j - gys) + gnx * gny * (k - gzs);
1361: cnt = 0;
1362: for (kk = kstart; kk < kend + 1; kk++) {
1363: for (jj = jstart; jj < jend + 1; jj++) {
1364: for (ii = istart; ii < iend + 1; ii++) {
1365: if (st == DMDA_STENCIL_BOX || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) { /* entries on star*/
1366: cols[cnt++] = nc * (slot + ii + gnx * jj + gnx * gny * kk);
1367: for (l = 1; l < nc; l++) {
1368: cols[cnt] = 1 + cols[cnt - 1];
1369: cnt++;
1370: }
1371: }
1372: }
1373: }
1374: }
1375: rows[0] = nc * slot;
1376: for (l = 1; l < nc; l++) rows[l] = 1 + rows[l - 1];
1377: PetscCall(MatSetValuesLocal(J, nc, rows, cnt, cols, NULL, INSERT_VALUES));
1378: }
1379: }
1380: }
1381: /* do not copy values to GPU since they are all zero and not yet needed there */
1382: PetscCall(MatBindToCPU(J, PETSC_TRUE));
1383: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
1384: PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
1385: if (bx == DM_BOUNDARY_NONE && by == DM_BOUNDARY_NONE && bz == DM_BOUNDARY_NONE) PetscCall(MatSetOption(J, MAT_SORTED_FULL, PETSC_FALSE));
1386: PetscCall(MatBindToCPU(J, PETSC_FALSE));
1387: PetscCall(MatSetOption(J, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE));
1388: }
1389: PetscCall(PetscFree2(rows, cols));
1390: PetscFunctionReturn(PETSC_SUCCESS);
1391: }
1393: PetscErrorCode DMCreateMatrix_DA_1d_MPIAIJ_Fill(DM da, Mat J)
1394: {
1395: DM_DA *dd = (DM_DA *)da->data;
1396: PetscInt xs, nx, i, j, gxs, gnx, row, k, l;
1397: PetscInt m, dim, s, *cols = NULL, nc, cnt, maxcnt = 0, *ocols;
1398: PetscInt *ofill = dd->ofill, *dfill = dd->dfill;
1399: DMBoundaryType bx;
1400: ISLocalToGlobalMapping ltog;
1401: PetscMPIInt rank, size;
1403: PetscFunctionBegin;
1404: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)da), &rank));
1405: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)da), &size));
1407: /*
1408: nc - number of components per grid point
1409: */
1410: PetscCall(DMDAGetInfo(da, &dim, &m, NULL, NULL, NULL, NULL, NULL, &nc, &s, &bx, NULL, NULL, NULL));
1411: PetscCheck(s <= 1, PetscObjectComm((PetscObject)da), PETSC_ERR_SUP, "Matrix creation for 1d not implemented correctly for stencil width larger than 1");
1412: PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &nx, NULL, NULL));
1413: PetscCall(DMDAGetGhostCorners(da, &gxs, NULL, NULL, &gnx, NULL, NULL));
1415: PetscCall(MatSetBlockSize(J, nc));
1416: PetscCall(PetscCalloc2(nx * nc, &cols, nx * nc, &ocols));
1418: /*
1419: note should be smaller for first and last process with no periodic
1420: does not handle dfill
1421: */
1422: cnt = 0;
1423: /* coupling with process to the left */
1424: for (i = 0; i < s; i++) {
1425: for (j = 0; j < nc; j++) {
1426: ocols[cnt] = ((rank == 0) ? 0 : (s - i) * (ofill[j + 1] - ofill[j]));
1427: cols[cnt] = dfill[j + 1] - dfill[j] + (s + i) * (ofill[j + 1] - ofill[j]);
1428: if (rank == 0 && (dd->bx == DM_BOUNDARY_PERIODIC)) {
1429: if (size > 1) ocols[cnt] += (s - i) * (ofill[j + 1] - ofill[j]);
1430: else cols[cnt] += (s - i) * (ofill[j + 1] - ofill[j]);
1431: }
1432: maxcnt = PetscMax(maxcnt, ocols[cnt] + cols[cnt]);
1433: cnt++;
1434: }
1435: }
1436: for (i = s; i < nx - s; i++) {
1437: for (j = 0; j < nc; j++) {
1438: cols[cnt] = dfill[j + 1] - dfill[j] + 2 * s * (ofill[j + 1] - ofill[j]);
1439: maxcnt = PetscMax(maxcnt, ocols[cnt] + cols[cnt]);
1440: cnt++;
1441: }
1442: }
1443: /* coupling with process to the right */
1444: for (i = nx - s; i < nx; i++) {
1445: for (j = 0; j < nc; j++) {
1446: ocols[cnt] = ((rank == (size - 1)) ? 0 : (i - nx + s + 1) * (ofill[j + 1] - ofill[j]));
1447: cols[cnt] = dfill[j + 1] - dfill[j] + (s + nx - i - 1) * (ofill[j + 1] - ofill[j]);
1448: if (rank == size - 1 && dd->bx == DM_BOUNDARY_PERIODIC) {
1449: if (size > 1) ocols[cnt] += (i - nx + s + 1) * (ofill[j + 1] - ofill[j]);
1450: else cols[cnt] += (i - nx + s + 1) * (ofill[j + 1] - ofill[j]);
1451: }
1452: maxcnt = PetscMax(maxcnt, ocols[cnt] + cols[cnt]);
1453: cnt++;
1454: }
1455: }
1457: PetscCall(MatSeqAIJSetPreallocation(J, 0, cols));
1458: PetscCall(MatMPIAIJSetPreallocation(J, 0, cols, 0, ocols));
1459: PetscCall(PetscFree2(cols, ocols));
1461: PetscCall(DMGetLocalToGlobalMapping(da, <og));
1462: PetscCall(MatSetLocalToGlobalMapping(J, ltog, ltog));
1464: /*
1465: For each node in the grid: we get the neighbors in the local (on processor ordering
1466: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1467: PETSc ordering.
1468: */
1469: if (!da->prealloc_only) {
1470: PetscCall(PetscMalloc1(maxcnt, &cols));
1471: row = xs * nc;
1472: /* coupling with process to the left */
1473: for (i = xs; i < xs + s; i++) {
1474: for (j = 0; j < nc; j++) {
1475: cnt = 0;
1476: if (rank) {
1477: for (l = 0; l < s; l++) {
1478: for (k = ofill[j]; k < ofill[j + 1]; k++) cols[cnt++] = (i - s + l) * nc + ofill[k];
1479: }
1480: }
1481: if (rank == 0 && (dd->bx == DM_BOUNDARY_PERIODIC)) {
1482: for (l = 0; l < s; l++) {
1483: for (k = ofill[j]; k < ofill[j + 1]; k++) cols[cnt++] = (m + i - s - l) * nc + ofill[k];
1484: }
1485: }
1486: if (dfill) {
1487: for (k = dfill[j]; k < dfill[j + 1]; k++) cols[cnt++] = i * nc + dfill[k];
1488: } else {
1489: for (k = 0; k < nc; k++) cols[cnt++] = i * nc + k;
1490: }
1491: for (l = 0; l < s; l++) {
1492: for (k = ofill[j]; k < ofill[j + 1]; k++) cols[cnt++] = (i + s - l) * nc + ofill[k];
1493: }
1494: PetscCall(MatSetValues(J, 1, &row, cnt, cols, NULL, INSERT_VALUES));
1495: row++;
1496: }
1497: }
1498: for (i = xs + s; i < xs + nx - s; i++) {
1499: for (j = 0; j < nc; j++) {
1500: cnt = 0;
1501: for (l = 0; l < s; l++) {
1502: for (k = ofill[j]; k < ofill[j + 1]; k++) cols[cnt++] = (i - s + l) * nc + ofill[k];
1503: }
1504: if (dfill) {
1505: for (k = dfill[j]; k < dfill[j + 1]; k++) cols[cnt++] = i * nc + dfill[k];
1506: } else {
1507: for (k = 0; k < nc; k++) cols[cnt++] = i * nc + k;
1508: }
1509: for (l = 0; l < s; l++) {
1510: for (k = ofill[j]; k < ofill[j + 1]; k++) cols[cnt++] = (i + s - l) * nc + ofill[k];
1511: }
1512: PetscCall(MatSetValues(J, 1, &row, cnt, cols, NULL, INSERT_VALUES));
1513: row++;
1514: }
1515: }
1516: /* coupling with process to the right */
1517: for (i = xs + nx - s; i < xs + nx; i++) {
1518: for (j = 0; j < nc; j++) {
1519: cnt = 0;
1520: for (l = 0; l < s; l++) {
1521: for (k = ofill[j]; k < ofill[j + 1]; k++) cols[cnt++] = (i - s + l) * nc + ofill[k];
1522: }
1523: if (dfill) {
1524: for (k = dfill[j]; k < dfill[j + 1]; k++) cols[cnt++] = i * nc + dfill[k];
1525: } else {
1526: for (k = 0; k < nc; k++) cols[cnt++] = i * nc + k;
1527: }
1528: if (rank < size - 1) {
1529: for (l = 0; l < s; l++) {
1530: for (k = ofill[j]; k < ofill[j + 1]; k++) cols[cnt++] = (i + s - l) * nc + ofill[k];
1531: }
1532: }
1533: if (rank == size - 1 && dd->bx == DM_BOUNDARY_PERIODIC) {
1534: for (l = 0; l < s; l++) {
1535: for (k = ofill[j]; k < ofill[j + 1]; k++) cols[cnt++] = (i - s - l - m + 2) * nc + ofill[k];
1536: }
1537: }
1538: PetscCall(MatSetValues(J, 1, &row, cnt, cols, NULL, INSERT_VALUES));
1539: row++;
1540: }
1541: }
1542: PetscCall(PetscFree(cols));
1543: /* do not copy values to GPU since they are all zero and not yet needed there */
1544: PetscCall(MatBindToCPU(J, PETSC_TRUE));
1545: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
1546: PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
1547: PetscCall(MatBindToCPU(J, PETSC_FALSE));
1548: PetscCall(MatSetOption(J, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE));
1549: }
1550: PetscFunctionReturn(PETSC_SUCCESS);
1551: }
1553: PetscErrorCode DMCreateMatrix_DA_1d_MPIAIJ(DM da, Mat J)
1554: {
1555: PetscInt xs, nx, i, i1, slot, gxs, gnx;
1556: PetscInt m, dim, s, *cols = NULL, nc, *rows = NULL, col, cnt, l;
1557: PetscInt istart, iend;
1558: DMBoundaryType bx;
1559: ISLocalToGlobalMapping ltog, mltog;
1561: PetscFunctionBegin;
1562: /*
1563: nc - number of components per grid point
1564: col - number of colors needed in one direction for single component problem
1565: */
1566: PetscCall(DMDAGetInfo(da, &dim, &m, NULL, NULL, NULL, NULL, NULL, &nc, &s, &bx, NULL, NULL, NULL));
1567: if (bx == DM_BOUNDARY_NONE) PetscCall(MatSetOption(J, MAT_SORTED_FULL, PETSC_TRUE));
1568: col = 2 * s + 1;
1570: PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &nx, NULL, NULL));
1571: PetscCall(DMDAGetGhostCorners(da, &gxs, NULL, NULL, &gnx, NULL, NULL));
1573: PetscCall(MatSetBlockSize(J, nc));
1574: PetscCall(MatSeqAIJSetPreallocation(J, col * nc, NULL));
1575: PetscCall(MatMPIAIJSetPreallocation(J, col * nc, NULL, col * nc, NULL));
1577: PetscCall(DMGetLocalToGlobalMapping(da, <og));
1578: PetscCall(MatGetLocalToGlobalMapping(J, &mltog, NULL));
1579: if (!mltog) PetscCall(MatSetLocalToGlobalMapping(J, ltog, ltog));
1581: /*
1582: For each node in the grid: we get the neighbors in the local (on processor ordering
1583: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1584: PETSc ordering.
1585: */
1586: if (!da->prealloc_only) {
1587: PetscCall(PetscMalloc2(nc, &rows, col * nc * nc, &cols));
1588: for (i = xs; i < xs + nx; i++) {
1589: istart = PetscMax(-s, gxs - i);
1590: iend = PetscMin(s, gxs + gnx - i - 1);
1591: slot = i - gxs;
1593: cnt = 0;
1594: for (i1 = istart; i1 < iend + 1; i1++) {
1595: cols[cnt++] = nc * (slot + i1);
1596: for (l = 1; l < nc; l++) {
1597: cols[cnt] = 1 + cols[cnt - 1];
1598: cnt++;
1599: }
1600: }
1601: rows[0] = nc * slot;
1602: for (l = 1; l < nc; l++) rows[l] = 1 + rows[l - 1];
1603: PetscCall(MatSetValuesLocal(J, nc, rows, cnt, cols, NULL, INSERT_VALUES));
1604: }
1605: /* do not copy values to GPU since they are all zero and not yet needed there */
1606: PetscCall(MatBindToCPU(J, PETSC_TRUE));
1607: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
1608: PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
1609: if (bx == DM_BOUNDARY_NONE) PetscCall(MatSetOption(J, MAT_SORTED_FULL, PETSC_FALSE));
1610: PetscCall(MatBindToCPU(J, PETSC_FALSE));
1611: PetscCall(MatSetOption(J, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE));
1612: PetscCall(PetscFree2(rows, cols));
1613: }
1614: PetscFunctionReturn(PETSC_SUCCESS);
1615: }
1617: PetscErrorCode DMCreateMatrix_DA_1d_SeqAIJ_NoPreallocation(DM da, Mat J)
1618: {
1619: PetscInt xs, nx, i, i1, slot, gxs, gnx;
1620: PetscInt m, dim, s, *cols = NULL, nc, *rows = NULL, col, cnt, l;
1621: PetscInt istart, iend;
1622: DMBoundaryType bx;
1623: ISLocalToGlobalMapping ltog, mltog;
1625: PetscFunctionBegin;
1626: /*
1627: nc - number of components per grid point
1628: col - number of colors needed in one direction for single component problem
1629: */
1630: PetscCall(DMDAGetInfo(da, &dim, &m, NULL, NULL, NULL, NULL, NULL, &nc, &s, &bx, NULL, NULL, NULL));
1631: col = 2 * s + 1;
1633: PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &nx, NULL, NULL));
1634: PetscCall(DMDAGetGhostCorners(da, &gxs, NULL, NULL, &gnx, NULL, NULL));
1636: PetscCall(MatSetBlockSize(J, nc));
1637: PetscCall(MatSeqAIJSetTotalPreallocation(J, nx * nc * col * nc));
1639: PetscCall(DMGetLocalToGlobalMapping(da, <og));
1640: PetscCall(MatGetLocalToGlobalMapping(J, &mltog, NULL));
1641: if (!mltog) PetscCall(MatSetLocalToGlobalMapping(J, ltog, ltog));
1643: /*
1644: For each node in the grid: we get the neighbors in the local (on processor ordering
1645: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1646: PETSc ordering.
1647: */
1648: if (!da->prealloc_only) {
1649: PetscCall(PetscMalloc2(nc, &rows, col * nc * nc, &cols));
1650: for (i = xs; i < xs + nx; i++) {
1651: istart = PetscMax(-s, gxs - i);
1652: iend = PetscMin(s, gxs + gnx - i - 1);
1653: slot = i - gxs;
1655: cnt = 0;
1656: for (i1 = istart; i1 < iend + 1; i1++) {
1657: cols[cnt++] = nc * (slot + i1);
1658: for (l = 1; l < nc; l++) {
1659: cols[cnt] = 1 + cols[cnt - 1];
1660: cnt++;
1661: }
1662: }
1663: rows[0] = nc * slot;
1664: for (l = 1; l < nc; l++) rows[l] = 1 + rows[l - 1];
1665: PetscCall(MatSetValuesLocal(J, nc, rows, cnt, cols, NULL, INSERT_VALUES));
1666: }
1667: /* do not copy values to GPU since they are all zero and not yet needed there */
1668: PetscCall(MatBindToCPU(J, PETSC_TRUE));
1669: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
1670: PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
1671: if (bx == DM_BOUNDARY_NONE) PetscCall(MatSetOption(J, MAT_SORTED_FULL, PETSC_FALSE));
1672: PetscCall(MatBindToCPU(J, PETSC_FALSE));
1673: PetscCall(MatSetOption(J, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE));
1674: PetscCall(PetscFree2(rows, cols));
1675: }
1676: PetscCall(MatSetOption(J, MAT_SORTED_FULL, PETSC_FALSE));
1677: PetscFunctionReturn(PETSC_SUCCESS);
1678: }
1680: PetscErrorCode DMCreateMatrix_DA_2d_MPIBAIJ(DM da, Mat J)
1681: {
1682: PetscInt xs, ys, nx, ny, i, j, slot, gxs, gys, gnx, gny;
1683: PetscInt m, n, dim, s, *cols, nc, col, cnt, *dnz, *onz;
1684: PetscInt istart, iend, jstart, jend, ii, jj;
1685: MPI_Comm comm;
1686: PetscScalar *values;
1687: DMBoundaryType bx, by;
1688: DMDAStencilType st;
1689: ISLocalToGlobalMapping ltog;
1691: PetscFunctionBegin;
1692: /*
1693: nc - number of components per grid point
1694: col - number of colors needed in one direction for single component problem
1695: */
1696: PetscCall(DMDAGetInfo(da, &dim, &m, &n, NULL, NULL, NULL, NULL, &nc, &s, &bx, &by, NULL, &st));
1697: col = 2 * s + 1;
1699: PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &nx, &ny, NULL));
1700: PetscCall(DMDAGetGhostCorners(da, &gxs, &gys, NULL, &gnx, &gny, NULL));
1701: PetscCall(PetscObjectGetComm((PetscObject)da, &comm));
1703: PetscCall(PetscMalloc1(col * col * nc * nc, &cols));
1705: PetscCall(DMGetLocalToGlobalMapping(da, <og));
1707: /* determine the matrix preallocation information */
1708: MatPreallocateBegin(comm, nx * ny, nx * ny, dnz, onz);
1709: for (i = xs; i < xs + nx; i++) {
1710: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s, -i));
1711: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s, m - i - 1));
1712: for (j = ys; j < ys + ny; j++) {
1713: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s, -j));
1714: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s, n - j - 1));
1715: slot = i - gxs + gnx * (j - gys);
1717: /* Find block columns in block row */
1718: cnt = 0;
1719: for (ii = istart; ii < iend + 1; ii++) {
1720: for (jj = jstart; jj < jend + 1; jj++) {
1721: if (st == DMDA_STENCIL_BOX || !ii || !jj) { /* BOX or on the STAR */
1722: cols[cnt++] = slot + ii + gnx * jj;
1723: }
1724: }
1725: }
1726: PetscCall(MatPreallocateSetLocalBlock(ltog, 1, &slot, ltog, cnt, cols, dnz, onz));
1727: }
1728: }
1729: PetscCall(MatSeqBAIJSetPreallocation(J, nc, 0, dnz));
1730: PetscCall(MatMPIBAIJSetPreallocation(J, nc, 0, dnz, 0, onz));
1731: MatPreallocateEnd(dnz, onz);
1733: PetscCall(MatSetLocalToGlobalMapping(J, ltog, ltog));
1735: /*
1736: For each node in the grid: we get the neighbors in the local (on processor ordering
1737: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1738: PETSc ordering.
1739: */
1740: if (!da->prealloc_only) {
1741: PetscCall(PetscCalloc1(col * col * nc * nc, &values));
1742: for (i = xs; i < xs + nx; i++) {
1743: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s, -i));
1744: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s, m - i - 1));
1745: for (j = ys; j < ys + ny; j++) {
1746: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s, -j));
1747: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s, n - j - 1));
1748: slot = i - gxs + gnx * (j - gys);
1749: cnt = 0;
1750: for (ii = istart; ii < iend + 1; ii++) {
1751: for (jj = jstart; jj < jend + 1; jj++) {
1752: if (st == DMDA_STENCIL_BOX || !ii || !jj) { /* BOX or on the STAR */
1753: cols[cnt++] = slot + ii + gnx * jj;
1754: }
1755: }
1756: }
1757: PetscCall(MatSetValuesBlockedLocal(J, 1, &slot, cnt, cols, values, INSERT_VALUES));
1758: }
1759: }
1760: PetscCall(PetscFree(values));
1761: /* do not copy values to GPU since they are all zero and not yet needed there */
1762: PetscCall(MatBindToCPU(J, PETSC_TRUE));
1763: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
1764: PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
1765: PetscCall(MatBindToCPU(J, PETSC_FALSE));
1766: PetscCall(MatSetOption(J, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE));
1767: }
1768: PetscCall(PetscFree(cols));
1769: PetscFunctionReturn(PETSC_SUCCESS);
1770: }
1772: PetscErrorCode DMCreateMatrix_DA_3d_MPIBAIJ(DM da, Mat J)
1773: {
1774: PetscInt xs, ys, nx, ny, i, j, slot, gxs, gys, gnx, gny;
1775: PetscInt m, n, dim, s, *cols, k, nc, col, cnt, p, *dnz, *onz;
1776: PetscInt istart, iend, jstart, jend, kstart, kend, zs, nz, gzs, gnz, ii, jj, kk;
1777: MPI_Comm comm;
1778: PetscScalar *values;
1779: DMBoundaryType bx, by, bz;
1780: DMDAStencilType st;
1781: ISLocalToGlobalMapping ltog;
1783: PetscFunctionBegin;
1784: /*
1785: nc - number of components per grid point
1786: col - number of colors needed in one direction for single component problem
1787: */
1788: PetscCall(DMDAGetInfo(da, &dim, &m, &n, &p, NULL, NULL, NULL, &nc, &s, &bx, &by, &bz, &st));
1789: col = 2 * s + 1;
1791: PetscCall(DMDAGetCorners(da, &xs, &ys, &zs, &nx, &ny, &nz));
1792: PetscCall(DMDAGetGhostCorners(da, &gxs, &gys, &gzs, &gnx, &gny, &gnz));
1793: PetscCall(PetscObjectGetComm((PetscObject)da, &comm));
1795: PetscCall(PetscMalloc1(col * col * col, &cols));
1797: PetscCall(DMGetLocalToGlobalMapping(da, <og));
1799: /* determine the matrix preallocation information */
1800: MatPreallocateBegin(comm, nx * ny * nz, nx * ny * nz, dnz, onz);
1801: for (i = xs; i < xs + nx; i++) {
1802: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s, -i));
1803: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s, m - i - 1));
1804: for (j = ys; j < ys + ny; j++) {
1805: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s, -j));
1806: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s, n - j - 1));
1807: for (k = zs; k < zs + nz; k++) {
1808: kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s, -k));
1809: kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s, p - k - 1));
1811: slot = i - gxs + gnx * (j - gys) + gnx * gny * (k - gzs);
1813: /* Find block columns in block row */
1814: cnt = 0;
1815: for (ii = istart; ii < iend + 1; ii++) {
1816: for (jj = jstart; jj < jend + 1; jj++) {
1817: for (kk = kstart; kk < kend + 1; kk++) {
1818: if (st == DMDA_STENCIL_BOX || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) { /* entries on star*/
1819: cols[cnt++] = slot + ii + gnx * jj + gnx * gny * kk;
1820: }
1821: }
1822: }
1823: }
1824: PetscCall(MatPreallocateSetLocalBlock(ltog, 1, &slot, ltog, cnt, cols, dnz, onz));
1825: }
1826: }
1827: }
1828: PetscCall(MatSeqBAIJSetPreallocation(J, nc, 0, dnz));
1829: PetscCall(MatMPIBAIJSetPreallocation(J, nc, 0, dnz, 0, onz));
1830: MatPreallocateEnd(dnz, onz);
1832: PetscCall(MatSetLocalToGlobalMapping(J, ltog, ltog));
1834: /*
1835: For each node in the grid: we get the neighbors in the local (on processor ordering
1836: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1837: PETSc ordering.
1838: */
1839: if (!da->prealloc_only) {
1840: PetscCall(PetscCalloc1(col * col * col * nc * nc, &values));
1841: for (i = xs; i < xs + nx; i++) {
1842: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s, -i));
1843: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s, m - i - 1));
1844: for (j = ys; j < ys + ny; j++) {
1845: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s, -j));
1846: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s, n - j - 1));
1847: for (k = zs; k < zs + nz; k++) {
1848: kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s, -k));
1849: kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s, p - k - 1));
1851: slot = i - gxs + gnx * (j - gys) + gnx * gny * (k - gzs);
1853: cnt = 0;
1854: for (ii = istart; ii < iend + 1; ii++) {
1855: for (jj = jstart; jj < jend + 1; jj++) {
1856: for (kk = kstart; kk < kend + 1; kk++) {
1857: if (st == DMDA_STENCIL_BOX || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) { /* entries on star*/
1858: cols[cnt++] = slot + ii + gnx * jj + gnx * gny * kk;
1859: }
1860: }
1861: }
1862: }
1863: PetscCall(MatSetValuesBlockedLocal(J, 1, &slot, cnt, cols, values, INSERT_VALUES));
1864: }
1865: }
1866: }
1867: PetscCall(PetscFree(values));
1868: /* do not copy values to GPU since they are all zero and not yet needed there */
1869: PetscCall(MatBindToCPU(J, PETSC_TRUE));
1870: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
1871: PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
1872: PetscCall(MatBindToCPU(J, PETSC_FALSE));
1873: PetscCall(MatSetOption(J, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE));
1874: }
1875: PetscCall(PetscFree(cols));
1876: PetscFunctionReturn(PETSC_SUCCESS);
1877: }
1879: /*
1880: This helper is for of SBAIJ preallocation, to discard the lower-triangular values which are difficult to
1881: identify in the local ordering with periodic domain.
1882: */
1883: static PetscErrorCode L2GFilterUpperTriangular(ISLocalToGlobalMapping ltog, PetscInt *row, PetscInt *cnt, PetscInt col[])
1884: {
1885: PetscInt i, n;
1887: PetscFunctionBegin;
1888: PetscCall(ISLocalToGlobalMappingApplyBlock(ltog, 1, row, row));
1889: PetscCall(ISLocalToGlobalMappingApplyBlock(ltog, *cnt, col, col));
1890: for (i = 0, n = 0; i < *cnt; i++) {
1891: if (col[i] >= *row) col[n++] = col[i];
1892: }
1893: *cnt = n;
1894: PetscFunctionReturn(PETSC_SUCCESS);
1895: }
1897: PetscErrorCode DMCreateMatrix_DA_2d_MPISBAIJ(DM da, Mat J)
1898: {
1899: PetscInt xs, ys, nx, ny, i, j, slot, gxs, gys, gnx, gny;
1900: PetscInt m, n, dim, s, *cols, nc, col, cnt, *dnz, *onz;
1901: PetscInt istart, iend, jstart, jend, ii, jj;
1902: MPI_Comm comm;
1903: PetscScalar *values;
1904: DMBoundaryType bx, by;
1905: DMDAStencilType st;
1906: ISLocalToGlobalMapping ltog;
1908: PetscFunctionBegin;
1909: /*
1910: nc - number of components per grid point
1911: col - number of colors needed in one direction for single component problem
1912: */
1913: PetscCall(DMDAGetInfo(da, &dim, &m, &n, NULL, NULL, NULL, NULL, &nc, &s, &bx, &by, NULL, &st));
1914: col = 2 * s + 1;
1916: PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &nx, &ny, NULL));
1917: PetscCall(DMDAGetGhostCorners(da, &gxs, &gys, NULL, &gnx, &gny, NULL));
1918: PetscCall(PetscObjectGetComm((PetscObject)da, &comm));
1920: PetscCall(PetscMalloc1(col * col * nc * nc, &cols));
1922: PetscCall(DMGetLocalToGlobalMapping(da, <og));
1924: /* determine the matrix preallocation information */
1925: MatPreallocateBegin(comm, nx * ny, nx * ny, dnz, onz);
1926: for (i = xs; i < xs + nx; i++) {
1927: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s, -i));
1928: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s, m - i - 1));
1929: for (j = ys; j < ys + ny; j++) {
1930: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s, -j));
1931: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s, n - j - 1));
1932: slot = i - gxs + gnx * (j - gys);
1934: /* Find block columns in block row */
1935: cnt = 0;
1936: for (ii = istart; ii < iend + 1; ii++) {
1937: for (jj = jstart; jj < jend + 1; jj++) {
1938: if (st == DMDA_STENCIL_BOX || !ii || !jj) cols[cnt++] = slot + ii + gnx * jj;
1939: }
1940: }
1941: PetscCall(L2GFilterUpperTriangular(ltog, &slot, &cnt, cols));
1942: PetscCall(MatPreallocateSymmetricSetBlock(slot, cnt, cols, dnz, onz));
1943: }
1944: }
1945: PetscCall(MatSeqSBAIJSetPreallocation(J, nc, 0, dnz));
1946: PetscCall(MatMPISBAIJSetPreallocation(J, nc, 0, dnz, 0, onz));
1947: MatPreallocateEnd(dnz, onz);
1949: PetscCall(MatSetLocalToGlobalMapping(J, ltog, ltog));
1951: /*
1952: For each node in the grid: we get the neighbors in the local (on processor ordering
1953: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
1954: PETSc ordering.
1955: */
1956: if (!da->prealloc_only) {
1957: PetscCall(PetscCalloc1(col * col * nc * nc, &values));
1958: for (i = xs; i < xs + nx; i++) {
1959: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s, -i));
1960: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s, m - i - 1));
1961: for (j = ys; j < ys + ny; j++) {
1962: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s, -j));
1963: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s, n - j - 1));
1964: slot = i - gxs + gnx * (j - gys);
1966: /* Find block columns in block row */
1967: cnt = 0;
1968: for (ii = istart; ii < iend + 1; ii++) {
1969: for (jj = jstart; jj < jend + 1; jj++) {
1970: if (st == DMDA_STENCIL_BOX || !ii || !jj) cols[cnt++] = slot + ii + gnx * jj;
1971: }
1972: }
1973: PetscCall(L2GFilterUpperTriangular(ltog, &slot, &cnt, cols));
1974: PetscCall(MatSetValuesBlocked(J, 1, &slot, cnt, cols, values, INSERT_VALUES));
1975: }
1976: }
1977: PetscCall(PetscFree(values));
1978: /* do not copy values to GPU since they are all zero and not yet needed there */
1979: PetscCall(MatBindToCPU(J, PETSC_TRUE));
1980: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
1981: PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
1982: PetscCall(MatBindToCPU(J, PETSC_FALSE));
1983: PetscCall(MatSetOption(J, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE));
1984: }
1985: PetscCall(PetscFree(cols));
1986: PetscFunctionReturn(PETSC_SUCCESS);
1987: }
1989: PetscErrorCode DMCreateMatrix_DA_3d_MPISBAIJ(DM da, Mat J)
1990: {
1991: PetscInt xs, ys, nx, ny, i, j, slot, gxs, gys, gnx, gny;
1992: PetscInt m, n, dim, s, *cols, k, nc, col, cnt, p, *dnz, *onz;
1993: PetscInt istart, iend, jstart, jend, kstart, kend, zs, nz, gzs, gnz, ii, jj, kk;
1994: MPI_Comm comm;
1995: PetscScalar *values;
1996: DMBoundaryType bx, by, bz;
1997: DMDAStencilType st;
1998: ISLocalToGlobalMapping ltog;
2000: PetscFunctionBegin;
2001: /*
2002: nc - number of components per grid point
2003: col - number of colors needed in one direction for single component problem
2004: */
2005: PetscCall(DMDAGetInfo(da, &dim, &m, &n, &p, NULL, NULL, NULL, &nc, &s, &bx, &by, &bz, &st));
2006: col = 2 * s + 1;
2008: PetscCall(DMDAGetCorners(da, &xs, &ys, &zs, &nx, &ny, &nz));
2009: PetscCall(DMDAGetGhostCorners(da, &gxs, &gys, &gzs, &gnx, &gny, &gnz));
2010: PetscCall(PetscObjectGetComm((PetscObject)da, &comm));
2012: /* create the matrix */
2013: PetscCall(PetscMalloc1(col * col * col, &cols));
2015: PetscCall(DMGetLocalToGlobalMapping(da, <og));
2017: /* determine the matrix preallocation information */
2018: MatPreallocateBegin(comm, nx * ny * nz, nx * ny * nz, dnz, onz);
2019: for (i = xs; i < xs + nx; i++) {
2020: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s, -i));
2021: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s, m - i - 1));
2022: for (j = ys; j < ys + ny; j++) {
2023: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s, -j));
2024: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s, n - j - 1));
2025: for (k = zs; k < zs + nz; k++) {
2026: kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s, -k));
2027: kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s, p - k - 1));
2029: slot = i - gxs + gnx * (j - gys) + gnx * gny * (k - gzs);
2031: /* Find block columns in block row */
2032: cnt = 0;
2033: for (ii = istart; ii < iend + 1; ii++) {
2034: for (jj = jstart; jj < jend + 1; jj++) {
2035: for (kk = kstart; kk < kend + 1; kk++) {
2036: if (st == DMDA_STENCIL_BOX || (!ii && !jj) || (!jj && !kk) || (!ii && !kk)) cols[cnt++] = slot + ii + gnx * jj + gnx * gny * kk;
2037: }
2038: }
2039: }
2040: PetscCall(L2GFilterUpperTriangular(ltog, &slot, &cnt, cols));
2041: PetscCall(MatPreallocateSymmetricSetBlock(slot, cnt, cols, dnz, onz));
2042: }
2043: }
2044: }
2045: PetscCall(MatSeqSBAIJSetPreallocation(J, nc, 0, dnz));
2046: PetscCall(MatMPISBAIJSetPreallocation(J, nc, 0, dnz, 0, onz));
2047: MatPreallocateEnd(dnz, onz);
2049: PetscCall(MatSetLocalToGlobalMapping(J, ltog, ltog));
2051: /*
2052: For each node in the grid: we get the neighbors in the local (on processor ordering
2053: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
2054: PETSc ordering.
2055: */
2056: if (!da->prealloc_only) {
2057: PetscCall(PetscCalloc1(col * col * col * nc * nc, &values));
2058: for (i = xs; i < xs + nx; i++) {
2059: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s, -i));
2060: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s, m - i - 1));
2061: for (j = ys; j < ys + ny; j++) {
2062: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s, -j));
2063: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s, n - j - 1));
2064: for (k = zs; k < zs + nz; k++) {
2065: kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s, -k));
2066: kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s, p - k - 1));
2068: slot = i - gxs + gnx * (j - gys) + gnx * gny * (k - gzs);
2070: cnt = 0;
2071: for (ii = istart; ii < iend + 1; ii++) {
2072: for (jj = jstart; jj < jend + 1; jj++) {
2073: for (kk = kstart; kk < kend + 1; kk++) {
2074: if (st == DMDA_STENCIL_BOX || (!ii && !jj) || (!jj && !kk) || (!ii && !kk)) cols[cnt++] = slot + ii + gnx * jj + gnx * gny * kk;
2075: }
2076: }
2077: }
2078: PetscCall(L2GFilterUpperTriangular(ltog, &slot, &cnt, cols));
2079: PetscCall(MatSetValuesBlocked(J, 1, &slot, cnt, cols, values, INSERT_VALUES));
2080: }
2081: }
2082: }
2083: PetscCall(PetscFree(values));
2084: /* do not copy values to GPU since they are all zero and not yet needed there */
2085: PetscCall(MatBindToCPU(J, PETSC_TRUE));
2086: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
2087: PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
2088: PetscCall(MatBindToCPU(J, PETSC_FALSE));
2089: PetscCall(MatSetOption(J, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE));
2090: }
2091: PetscCall(PetscFree(cols));
2092: PetscFunctionReturn(PETSC_SUCCESS);
2093: }
2095: PetscErrorCode DMCreateMatrix_DA_3d_MPIAIJ_Fill(DM da, Mat J)
2096: {
2097: PetscInt xs, ys, nx, ny, i, j, slot, gxs, gys, gnx, gny;
2098: PetscInt m, n, dim, s, *cols, k, nc, row, col, cnt, maxcnt = 0, l, p, *dnz, *onz;
2099: PetscInt istart, iend, jstart, jend, kstart, kend, zs, nz, gzs, gnz, ii, jj, kk, M, N, P;
2100: DM_DA *dd = (DM_DA *)da->data;
2101: PetscInt ifill_col, *dfill = dd->dfill, *ofill = dd->ofill;
2102: MPI_Comm comm;
2103: PetscScalar *values;
2104: DMBoundaryType bx, by, bz;
2105: ISLocalToGlobalMapping ltog;
2106: DMDAStencilType st;
2107: PetscBool removedups = PETSC_FALSE;
2109: PetscFunctionBegin;
2110: /*
2111: nc - number of components per grid point
2112: col - number of colors needed in one direction for single component problem
2113: */
2114: PetscCall(DMDAGetInfo(da, &dim, &m, &n, &p, &M, &N, &P, &nc, &s, &bx, &by, &bz, &st));
2115: col = 2 * s + 1;
2117: /*
2118: With one processor in periodic domains in a skinny dimension the code will label nonzero columns multiple times
2119: because of "wrapping" around the end of the domain hitting an entry already counted in the other direction.
2120: */
2121: if (M == 1 && 2 * s >= m) removedups = PETSC_TRUE;
2122: if (N == 1 && 2 * s >= n) removedups = PETSC_TRUE;
2123: if (P == 1 && 2 * s >= p) removedups = PETSC_TRUE;
2125: PetscCall(DMDAGetCorners(da, &xs, &ys, &zs, &nx, &ny, &nz));
2126: PetscCall(DMDAGetGhostCorners(da, &gxs, &gys, &gzs, &gnx, &gny, &gnz));
2127: PetscCall(PetscObjectGetComm((PetscObject)da, &comm));
2129: PetscCall(PetscMalloc1(col * col * col * nc, &cols));
2130: PetscCall(DMGetLocalToGlobalMapping(da, <og));
2132: /* determine the matrix preallocation information */
2133: MatPreallocateBegin(comm, nc * nx * ny * nz, nc * nx * ny * nz, dnz, onz);
2135: PetscCall(MatSetBlockSize(J, nc));
2136: for (i = xs; i < xs + nx; i++) {
2137: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s, -i));
2138: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s, m - i - 1));
2139: for (j = ys; j < ys + ny; j++) {
2140: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s, -j));
2141: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s, n - j - 1));
2142: for (k = zs; k < zs + nz; k++) {
2143: kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s, -k));
2144: kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s, p - k - 1));
2146: slot = i - gxs + gnx * (j - gys) + gnx * gny * (k - gzs);
2148: for (l = 0; l < nc; l++) {
2149: cnt = 0;
2150: for (ii = istart; ii < iend + 1; ii++) {
2151: for (jj = jstart; jj < jend + 1; jj++) {
2152: for (kk = kstart; kk < kend + 1; kk++) {
2153: if (ii || jj || kk) {
2154: if (st == DMDA_STENCIL_BOX || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) { /* entries on star*/
2155: for (ifill_col = ofill[l]; ifill_col < ofill[l + 1]; ifill_col++) cols[cnt++] = ofill[ifill_col] + nc * (slot + ii + gnx * jj + gnx * gny * kk);
2156: }
2157: } else {
2158: if (dfill) {
2159: for (ifill_col = dfill[l]; ifill_col < dfill[l + 1]; ifill_col++) cols[cnt++] = dfill[ifill_col] + nc * (slot + ii + gnx * jj + gnx * gny * kk);
2160: } else {
2161: for (ifill_col = 0; ifill_col < nc; ifill_col++) cols[cnt++] = ifill_col + nc * (slot + ii + gnx * jj + gnx * gny * kk);
2162: }
2163: }
2164: }
2165: }
2166: }
2167: row = l + nc * slot;
2168: maxcnt = PetscMax(maxcnt, cnt);
2169: if (removedups) PetscCall(MatPreallocateSetLocalRemoveDups(ltog, 1, &row, ltog, cnt, cols, dnz, onz));
2170: else PetscCall(MatPreallocateSetLocal(ltog, 1, &row, ltog, cnt, cols, dnz, onz));
2171: }
2172: }
2173: }
2174: }
2175: PetscCall(MatSeqAIJSetPreallocation(J, 0, dnz));
2176: PetscCall(MatMPIAIJSetPreallocation(J, 0, dnz, 0, onz));
2177: MatPreallocateEnd(dnz, onz);
2178: PetscCall(MatSetLocalToGlobalMapping(J, ltog, ltog));
2180: /*
2181: For each node in the grid: we get the neighbors in the local (on processor ordering
2182: that includes the ghost points) then MatSetValuesLocal() maps those indices to the global
2183: PETSc ordering.
2184: */
2185: if (!da->prealloc_only) {
2186: PetscCall(PetscCalloc1(maxcnt, &values));
2187: for (i = xs; i < xs + nx; i++) {
2188: istart = (bx == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s, -i));
2189: iend = (bx == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s, m - i - 1));
2190: for (j = ys; j < ys + ny; j++) {
2191: jstart = (by == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s, -j));
2192: jend = (by == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s, n - j - 1));
2193: for (k = zs; k < zs + nz; k++) {
2194: kstart = (bz == DM_BOUNDARY_PERIODIC) ? -s : (PetscMax(-s, -k));
2195: kend = (bz == DM_BOUNDARY_PERIODIC) ? s : (PetscMin(s, p - k - 1));
2197: slot = i - gxs + gnx * (j - gys) + gnx * gny * (k - gzs);
2199: for (l = 0; l < nc; l++) {
2200: cnt = 0;
2201: for (ii = istart; ii < iend + 1; ii++) {
2202: for (jj = jstart; jj < jend + 1; jj++) {
2203: for (kk = kstart; kk < kend + 1; kk++) {
2204: if (ii || jj || kk) {
2205: if (st == DMDA_STENCIL_BOX || ((!ii && !jj) || (!jj && !kk) || (!ii && !kk))) { /* entries on star*/
2206: for (ifill_col = ofill[l]; ifill_col < ofill[l + 1]; ifill_col++) cols[cnt++] = ofill[ifill_col] + nc * (slot + ii + gnx * jj + gnx * gny * kk);
2207: }
2208: } else {
2209: if (dfill) {
2210: for (ifill_col = dfill[l]; ifill_col < dfill[l + 1]; ifill_col++) cols[cnt++] = dfill[ifill_col] + nc * (slot + ii + gnx * jj + gnx * gny * kk);
2211: } else {
2212: for (ifill_col = 0; ifill_col < nc; ifill_col++) cols[cnt++] = ifill_col + nc * (slot + ii + gnx * jj + gnx * gny * kk);
2213: }
2214: }
2215: }
2216: }
2217: }
2218: row = l + nc * slot;
2219: PetscCall(MatSetValuesLocal(J, 1, &row, cnt, cols, values, INSERT_VALUES));
2220: }
2221: }
2222: }
2223: }
2224: PetscCall(PetscFree(values));
2225: /* do not copy values to GPU since they are all zero and not yet needed there */
2226: PetscCall(MatBindToCPU(J, PETSC_TRUE));
2227: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
2228: PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
2229: PetscCall(MatBindToCPU(J, PETSC_FALSE));
2230: PetscCall(MatSetOption(J, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE));
2231: }
2232: PetscCall(PetscFree(cols));
2233: PetscFunctionReturn(PETSC_SUCCESS);
2234: }