Actual source code: bjacobi.c
1: /*
2: Defines a block Jacobi preconditioner.
3: */
5: #include <../src/ksp/pc/impls/bjacobi/bjacobi.h>
7: static PetscErrorCode PCSetUp_BJacobi_Singleblock(PC, Mat, Mat);
8: static PetscErrorCode PCSetUp_BJacobi_Multiblock(PC, Mat, Mat);
9: static PetscErrorCode PCSetUp_BJacobi_Multiproc(PC);
11: static PetscErrorCode PCSetUp_BJacobi(PC pc)
12: {
13: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
14: Mat mat = pc->mat, pmat = pc->pmat;
15: PetscBool hasop;
16: PetscInt N, M, start, i, sum, end;
17: PetscInt bs, i_start = -1, i_end = -1;
18: PetscMPIInt rank, size;
20: PetscFunctionBegin;
21: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)pc), &rank));
22: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size));
23: PetscCall(MatGetLocalSize(pc->pmat, &M, &N));
24: PetscCall(MatGetBlockSize(pc->pmat, &bs));
26: if (jac->n > 0 && jac->n < size) {
27: PetscCall(PCSetUp_BJacobi_Multiproc(pc));
28: PetscFunctionReturn(PETSC_SUCCESS);
29: }
31: /* Determines the number of blocks assigned to each processor */
32: /* local block count given */
33: if (jac->n_local > 0 && jac->n < 0) {
34: PetscCallMPI(MPIU_Allreduce(&jac->n_local, &jac->n, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)pc)));
35: if (jac->l_lens) { /* check that user set these correctly */
36: sum = 0;
37: for (i = 0; i < jac->n_local; i++) {
38: PetscCheck(jac->l_lens[i] / bs * bs == jac->l_lens[i], PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Mat blocksize doesn't match block Jacobi layout");
39: sum += jac->l_lens[i];
40: }
41: PetscCheck(sum == M, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Local lens set incorrectly");
42: } else {
43: PetscCall(PetscMalloc1(jac->n_local, &jac->l_lens));
44: for (i = 0; i < jac->n_local; i++) jac->l_lens[i] = bs * ((M / bs) / jac->n_local + (((M / bs) % jac->n_local) > i));
45: }
46: } else if (jac->n > 0 && jac->n_local < 0) { /* global block count given */
47: /* global blocks given: determine which ones are local */
48: if (jac->g_lens) {
49: /* check if the g_lens is has valid entries */
50: for (i = 0; i < jac->n; i++) {
51: PetscCheck(jac->g_lens[i], PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Zero block not allowed");
52: PetscCheck(jac->g_lens[i] / bs * bs == jac->g_lens[i], PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Mat blocksize doesn't match block Jacobi layout");
53: }
54: if (size == 1) {
55: jac->n_local = jac->n;
56: PetscCall(PetscMalloc1(jac->n_local, &jac->l_lens));
57: PetscCall(PetscArraycpy(jac->l_lens, jac->g_lens, jac->n_local));
58: /* check that user set these correctly */
59: sum = 0;
60: for (i = 0; i < jac->n_local; i++) sum += jac->l_lens[i];
61: PetscCheck(sum == M, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Global lens set incorrectly");
62: } else {
63: PetscCall(MatGetOwnershipRange(pc->pmat, &start, &end));
64: /* loop over blocks determining first one owned by me */
65: sum = 0;
66: for (i = 0; i < jac->n + 1; i++) {
67: if (sum == start) {
68: i_start = i;
69: goto start_1;
70: }
71: if (i < jac->n) sum += jac->g_lens[i];
72: }
73: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Block sizes used in PCBJacobiSetTotalBlocks()\nare not compatible with parallel matrix layout");
74: start_1:
75: for (i = i_start; i < jac->n + 1; i++) {
76: if (sum == end) {
77: i_end = i;
78: goto end_1;
79: }
80: if (i < jac->n) sum += jac->g_lens[i];
81: }
82: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Block sizes used in PCBJacobiSetTotalBlocks()\nare not compatible with parallel matrix layout");
83: end_1:
84: jac->n_local = i_end - i_start;
85: PetscCall(PetscMalloc1(jac->n_local, &jac->l_lens));
86: PetscCall(PetscArraycpy(jac->l_lens, jac->g_lens + i_start, jac->n_local));
87: }
88: } else { /* no global blocks given, determine then using default layout */
89: jac->n_local = jac->n / size + ((jac->n % size) > rank);
90: PetscCall(PetscMalloc1(jac->n_local, &jac->l_lens));
91: for (i = 0; i < jac->n_local; i++) {
92: jac->l_lens[i] = ((M / bs) / jac->n_local + (((M / bs) % jac->n_local) > i)) * bs;
93: PetscCheck(jac->l_lens[i], PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Too many blocks given");
94: }
95: }
96: } else if (jac->n < 0 && jac->n_local < 0) { /* no blocks given */
97: jac->n = size;
98: jac->n_local = 1;
99: PetscCall(PetscMalloc1(1, &jac->l_lens));
100: jac->l_lens[0] = M;
101: } else { /* jac->n > 0 && jac->n_local > 0 */
102: if (!jac->l_lens) {
103: PetscCall(PetscMalloc1(jac->n_local, &jac->l_lens));
104: for (i = 0; i < jac->n_local; i++) jac->l_lens[i] = bs * ((M / bs) / jac->n_local + (((M / bs) % jac->n_local) > i));
105: }
106: }
107: PetscCheck(jac->n_local >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of blocks is less than number of processors");
109: /* Determines mat and pmat */
110: PetscCall(MatHasOperation(pc->mat, MATOP_GET_DIAGONAL_BLOCK, &hasop));
111: if (!hasop && size == 1) {
112: mat = pc->mat;
113: pmat = pc->pmat;
114: } else {
115: if (pc->useAmat) {
116: /* use block from Amat matrix, not Pmat for local MatMult() */
117: PetscCall(MatGetDiagonalBlock(pc->mat, &mat));
118: }
119: if (pc->pmat != pc->mat || !pc->useAmat) PetscCall(MatGetDiagonalBlock(pc->pmat, &pmat));
120: else pmat = mat;
121: }
123: /*
124: Setup code depends on the number of blocks
125: */
126: if (jac->n_local == 1) PetscCall(PCSetUp_BJacobi_Singleblock(pc, mat, pmat));
127: else PetscCall(PCSetUp_BJacobi_Multiblock(pc, mat, pmat));
128: PetscFunctionReturn(PETSC_SUCCESS);
129: }
131: /* Default destroy, if it has never been setup */
132: static PetscErrorCode PCDestroy_BJacobi(PC pc)
133: {
134: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
136: PetscFunctionBegin;
137: PetscCall(PetscFree(jac->g_lens));
138: PetscCall(PetscFree(jac->l_lens));
139: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBJacobiGetSubKSP_C", NULL));
140: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBJacobiSetTotalBlocks_C", NULL));
141: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBJacobiGetTotalBlocks_C", NULL));
142: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBJacobiSetLocalBlocks_C", NULL));
143: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBJacobiGetLocalBlocks_C", NULL));
144: PetscCall(PetscFree(pc->data));
145: PetscFunctionReturn(PETSC_SUCCESS);
146: }
148: static PetscErrorCode PCSetFromOptions_BJacobi(PC pc, PetscOptionItems PetscOptionsObject)
149: {
150: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
151: PetscInt blocks;
152: PetscBool flg;
154: PetscFunctionBegin;
155: PetscOptionsHeadBegin(PetscOptionsObject, "Block Jacobi options");
156: PetscCall(PetscOptionsInt("-pc_bjacobi_blocks", "Total number of blocks", "PCBJacobiSetTotalBlocks", jac->n, &blocks, &flg));
157: if (flg) PetscCall(PCBJacobiSetTotalBlocks(pc, blocks, NULL));
158: PetscCall(PetscOptionsInt("-pc_bjacobi_local_blocks", "Local number of blocks", "PCBJacobiSetLocalBlocks", jac->n_local, &blocks, &flg));
159: if (flg) PetscCall(PCBJacobiSetLocalBlocks(pc, blocks, NULL));
160: if (jac->ksp) {
161: /* The sub-KSP has already been set up (e.g., PCSetUp_BJacobi_Singleblock), but KSPSetFromOptions was not called
162: * unless we had already been called. */
163: for (PetscInt i = 0; i < jac->n_local; i++) PetscCall(KSPSetFromOptions(jac->ksp[i]));
164: }
165: PetscOptionsHeadEnd();
166: PetscFunctionReturn(PETSC_SUCCESS);
167: }
169: #include <petscdraw.h>
170: static PetscErrorCode PCView_BJacobi(PC pc, PetscViewer viewer)
171: {
172: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
173: PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc *)jac->data;
174: PetscMPIInt rank;
175: PetscBool isascii, isstring, isdraw;
176: PetscViewer sviewer;
177: PetscViewerFormat format;
178: const char *prefix;
180: PetscFunctionBegin;
181: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
182: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSTRING, &isstring));
183: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
184: if (isascii) {
185: if (pc->useAmat) PetscCall(PetscViewerASCIIPrintf(viewer, " using Amat local matrix, number of blocks = %" PetscInt_FMT "\n", jac->n));
186: PetscCall(PetscViewerASCIIPrintf(viewer, " number of blocks = %" PetscInt_FMT "\n", jac->n));
187: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)pc), &rank));
188: PetscCall(PetscViewerGetFormat(viewer, &format));
189: if (format != PETSC_VIEWER_ASCII_INFO_DETAIL) {
190: PetscCall(PetscViewerASCIIPrintf(viewer, " Local solver information for first block is in the following KSP and PC objects on rank 0:\n"));
191: PetscCall(PCGetOptionsPrefix(pc, &prefix));
192: PetscCall(PetscViewerASCIIPrintf(viewer, " Use -%sksp_view ::ascii_info_detail to display information for all blocks\n", prefix ? prefix : ""));
193: if (jac->ksp && !jac->psubcomm) {
194: PetscCall(PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &sviewer));
195: if (rank == 0) {
196: PetscCall(PetscViewerASCIIPushTab(sviewer));
197: PetscCall(KSPView(jac->ksp[0], sviewer));
198: PetscCall(PetscViewerASCIIPopTab(sviewer));
199: }
200: PetscCall(PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &sviewer));
201: /* extra call needed because of the two calls to PetscViewerASCIIPushSynchronized() in PetscViewerGetSubViewer() */
202: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
203: } else if (mpjac && jac->ksp && mpjac->psubcomm) {
204: PetscCall(PetscViewerGetSubViewer(viewer, mpjac->psubcomm->child, &sviewer));
205: if (!mpjac->psubcomm->color) {
206: PetscCall(PetscViewerASCIIPushTab(sviewer));
207: PetscCall(KSPView(*jac->ksp, sviewer));
208: PetscCall(PetscViewerASCIIPopTab(sviewer));
209: }
210: PetscCall(PetscViewerRestoreSubViewer(viewer, mpjac->psubcomm->child, &sviewer));
211: /* extra call needed because of the two calls to PetscViewerASCIIPushSynchronized() in PetscViewerGetSubViewer() */
212: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
213: }
214: } else {
215: PetscInt n_global;
216: PetscCallMPI(MPIU_Allreduce(&jac->n_local, &n_global, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)pc)));
217: PetscCall(PetscViewerASCIIPushSynchronized(viewer));
218: PetscCall(PetscViewerASCIIPrintf(viewer, " Local solver information for each block is in the following KSP and PC objects:\n"));
219: PetscCall(PetscViewerASCIIPushTab(viewer));
220: PetscCall(PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &sviewer));
221: PetscCall(PetscViewerASCIIPrintf(sviewer, "[%d] number of local blocks = %" PetscInt_FMT ", first local block number = %" PetscInt_FMT "\n", rank, jac->n_local, jac->first_local));
222: for (PetscInt i = 0; i < jac->n_local; i++) {
223: PetscCall(PetscViewerASCIIPrintf(sviewer, "[%d] local block number %" PetscInt_FMT "\n", rank, i));
224: PetscCall(KSPView(jac->ksp[i], sviewer));
225: PetscCall(PetscViewerASCIIPrintf(sviewer, "- - - - - - - - - - - - - - - - - -\n"));
226: }
227: PetscCall(PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &sviewer));
228: PetscCall(PetscViewerASCIIPopTab(viewer));
229: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
230: }
231: } else if (isstring) {
232: PetscCall(PetscViewerStringSPrintf(viewer, " blks=%" PetscInt_FMT, jac->n));
233: PetscCall(PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &sviewer));
234: if (jac->ksp) PetscCall(KSPView(jac->ksp[0], sviewer));
235: PetscCall(PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &sviewer));
236: } else if (isdraw) {
237: PetscDraw draw;
238: char str[25];
239: PetscReal x, y, bottom, h;
241: PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
242: PetscCall(PetscDrawGetCurrentPoint(draw, &x, &y));
243: PetscCall(PetscSNPrintf(str, 25, "Number blocks %" PetscInt_FMT, jac->n));
244: PetscCall(PetscDrawStringBoxed(draw, x, y, PETSC_DRAW_RED, PETSC_DRAW_BLACK, str, NULL, &h));
245: bottom = y - h;
246: PetscCall(PetscDrawPushCurrentPoint(draw, x, bottom));
247: /* warning the communicator on viewer is different then on ksp in parallel */
248: if (jac->ksp) PetscCall(KSPView(jac->ksp[0], viewer));
249: PetscCall(PetscDrawPopCurrentPoint(draw));
250: }
251: PetscFunctionReturn(PETSC_SUCCESS);
252: }
254: static PetscErrorCode PCBJacobiGetSubKSP_BJacobi(PC pc, PetscInt *n_local, PetscInt *first_local, KSP **ksp)
255: {
256: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
258: PetscFunctionBegin;
259: PetscCheck(pc->setupcalled, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Must call KSPSetUp() or PCSetUp() first");
261: if (n_local) *n_local = jac->n_local;
262: if (first_local) *first_local = jac->first_local;
263: if (ksp) *ksp = jac->ksp;
264: PetscFunctionReturn(PETSC_SUCCESS);
265: }
267: static PetscErrorCode PCBJacobiSetTotalBlocks_BJacobi(PC pc, PetscInt blocks, const PetscInt *lens)
268: {
269: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
271: PetscFunctionBegin;
272: PetscCheck(!pc->setupcalled || jac->n == blocks, PetscObjectComm((PetscObject)pc), PETSC_ERR_ORDER, "Cannot alter number of blocks after PCSetUp()/KSPSetUp() has been called");
273: jac->n = blocks;
274: if (!lens) jac->g_lens = NULL;
275: else {
276: PetscCall(PetscMalloc1(blocks, &jac->g_lens));
277: PetscCall(PetscArraycpy(jac->g_lens, lens, blocks));
278: }
279: PetscFunctionReturn(PETSC_SUCCESS);
280: }
282: static PetscErrorCode PCBJacobiGetTotalBlocks_BJacobi(PC pc, PetscInt *blocks, const PetscInt *lens[])
283: {
284: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
286: PetscFunctionBegin;
287: *blocks = jac->n;
288: if (lens) *lens = jac->g_lens;
289: PetscFunctionReturn(PETSC_SUCCESS);
290: }
292: static PetscErrorCode PCBJacobiSetLocalBlocks_BJacobi(PC pc, PetscInt blocks, const PetscInt lens[])
293: {
294: PC_BJacobi *jac;
296: PetscFunctionBegin;
297: jac = (PC_BJacobi *)pc->data;
299: jac->n_local = blocks;
300: if (!lens) jac->l_lens = NULL;
301: else {
302: PetscCall(PetscMalloc1(blocks, &jac->l_lens));
303: PetscCall(PetscArraycpy(jac->l_lens, lens, blocks));
304: }
305: PetscFunctionReturn(PETSC_SUCCESS);
306: }
308: static PetscErrorCode PCBJacobiGetLocalBlocks_BJacobi(PC pc, PetscInt *blocks, const PetscInt *lens[])
309: {
310: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
312: PetscFunctionBegin;
313: *blocks = jac->n_local;
314: if (lens) *lens = jac->l_lens;
315: PetscFunctionReturn(PETSC_SUCCESS);
316: }
318: /*@
319: PCBJacobiGetSubKSP - Gets the local `KSP` contexts for all blocks on
320: this processor.
322: Not Collective
324: Input Parameter:
325: . pc - the preconditioner context
327: Output Parameters:
328: + n_local - the number of blocks on this processor, or NULL
329: . first_local - the global number of the first block on this processor, or NULL
330: - ksp - the array of KSP contexts
332: Level: advanced
334: Notes:
335: After `PCBJacobiGetSubKSP()` the array of `KSP` contexts is not to be freed.
337: Currently for some matrix implementations only 1 block per processor
338: is supported.
340: You must call `KSPSetUp()` or `PCSetUp()` before calling `PCBJacobiGetSubKSP()`.
342: Fortran Note:
343: Call `PCBJacobiRestoreSubKSP()` when you no longer need access to the array of `KSP`
345: .seealso: [](ch_ksp), `PCBJACOBI`, `PCASM`, `PCASMGetSubKSP()`
346: @*/
347: PetscErrorCode PCBJacobiGetSubKSP(PC pc, PetscInt *n_local, PetscInt *first_local, KSP *ksp[])
348: {
349: PetscFunctionBegin;
351: PetscUseMethod(pc, "PCBJacobiGetSubKSP_C", (PC, PetscInt *, PetscInt *, KSP **), (pc, n_local, first_local, ksp));
352: PetscFunctionReturn(PETSC_SUCCESS);
353: }
355: /*@
356: PCBJacobiSetTotalBlocks - Sets the global number of blocks for the block
357: Jacobi preconditioner.
359: Collective
361: Input Parameters:
362: + pc - the preconditioner context
363: . blocks - the number of blocks
364: - lens - [optional] integer array containing the size of each block
366: Options Database Key:
367: . -pc_bjacobi_blocks blocks - Sets the number of global blocks
369: Level: intermediate
371: Note:
372: Currently only a limited number of blocking configurations are supported.
373: All processors sharing the `PC` must call this routine with the same data.
375: .seealso: [](ch_ksp), `PCBJACOBI`, `PCSetUseAmat()`, `PCBJacobiSetLocalBlocks()`
376: @*/
377: PetscErrorCode PCBJacobiSetTotalBlocks(PC pc, PetscInt blocks, const PetscInt lens[])
378: {
379: PetscFunctionBegin;
381: PetscCheck(blocks > 0, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Must have positive blocks");
382: PetscTryMethod(pc, "PCBJacobiSetTotalBlocks_C", (PC, PetscInt, const PetscInt[]), (pc, blocks, lens));
383: PetscFunctionReturn(PETSC_SUCCESS);
384: }
386: /*@
387: PCBJacobiGetTotalBlocks - Gets the global number of blocks for the block
388: Jacobi, `PCBJACOBI`, preconditioner.
390: Not Collective
392: Input Parameter:
393: . pc - the preconditioner context
395: Output Parameters:
396: + blocks - the number of blocks
397: - lens - integer array containing the size of each block
399: Level: intermediate
401: .seealso: [](ch_ksp), `PCBJACOBI`, `PCSetUseAmat()`, `PCBJacobiGetLocalBlocks()`
402: @*/
403: PetscErrorCode PCBJacobiGetTotalBlocks(PC pc, PetscInt *blocks, const PetscInt *lens[])
404: {
405: PetscFunctionBegin;
407: PetscAssertPointer(blocks, 2);
408: PetscUseMethod(pc, "PCBJacobiGetTotalBlocks_C", (PC, PetscInt *, const PetscInt *[]), (pc, blocks, lens));
409: PetscFunctionReturn(PETSC_SUCCESS);
410: }
412: /*@
413: PCBJacobiSetLocalBlocks - Sets the local number of blocks for the block
414: Jacobi, `PCBJACOBI`, preconditioner.
416: Not Collective
418: Input Parameters:
419: + pc - the preconditioner context
420: . blocks - the number of blocks
421: - lens - [optional] integer array containing size of each block
423: Options Database Key:
424: . -pc_bjacobi_local_blocks blocks - Sets the number of local blocks
426: Level: intermediate
428: Note:
429: Currently only a limited number of blocking configurations are supported.
431: .seealso: [](ch_ksp), `PCBJACOBI`, `PCSetUseAmat()`, `PCBJacobiSetTotalBlocks()`
432: @*/
433: PetscErrorCode PCBJacobiSetLocalBlocks(PC pc, PetscInt blocks, const PetscInt lens[])
434: {
435: PetscFunctionBegin;
437: PetscCheck(blocks >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must have nonegative blocks");
438: PetscTryMethod(pc, "PCBJacobiSetLocalBlocks_C", (PC, PetscInt, const PetscInt[]), (pc, blocks, lens));
439: PetscFunctionReturn(PETSC_SUCCESS);
440: }
442: /*@
443: PCBJacobiGetLocalBlocks - Gets the local number of blocks for the block
444: Jacobi, `PCBJACOBI`, preconditioner.
446: Not Collective
448: Input Parameters:
449: + pc - the preconditioner context
450: . blocks - the number of blocks
451: - lens - [optional] integer array containing size of each block
453: Level: intermediate
455: Note:
456: Currently only a limited number of blocking configurations are supported.
458: .seealso: [](ch_ksp), `PCBJACOBI`, `PCSetUseAmat()`, `PCBJacobiGetTotalBlocks()`
459: @*/
460: PetscErrorCode PCBJacobiGetLocalBlocks(PC pc, PetscInt *blocks, const PetscInt *lens[])
461: {
462: PetscFunctionBegin;
464: PetscAssertPointer(blocks, 2);
465: PetscUseMethod(pc, "PCBJacobiGetLocalBlocks_C", (PC, PetscInt *, const PetscInt *[]), (pc, blocks, lens));
466: PetscFunctionReturn(PETSC_SUCCESS);
467: }
469: /*MC
470: PCBJACOBI - Use block Jacobi preconditioning, each block is (approximately) solved with its own `KSP` object.
472: Options Database Keys:
473: + -pc_use_amat (true|false) - use `Amat` to apply block of operator in inner Krylov method
474: - -pc_bjacobi_blocks n - use n total blocks
476: Level: beginner
478: Notes:
479: See `PCJACOBI` for diagonal Jacobi, `PCVPBJACOBI` for variable point block, and `PCPBJACOBI` for fixed size point block
481: Each processor can have one or more blocks, or a single block can be shared by several processes. Defaults to one block per processor.
483: To set options on the solvers for each block append -sub_ to all the `KSP` and `PC`
484: options database keys. For example, -sub_pc_type ilu -sub_pc_factor_levels 1 -sub_ksp_type preonly
486: To set the options on the solvers separate for each block call `PCBJacobiGetSubKSP()`
487: and set the options directly on the resulting `KSP` object (you can access its `PC`
488: `KSPGetPC()`)
490: For GPU-based vectors (`VECCUDA`, `VECViennaCL`) it is recommended to use exactly one block per MPI process for best
491: performance. Different block partitioning may lead to additional data transfers
492: between host and GPU that lead to degraded performance.
494: When multiple processes share a single block, each block encompasses exactly all the unknowns owned its set of processes.
496: .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`,
497: `PCASM`, `PCSetUseAmat()`, `PCGetUseAmat()`, `PCBJacobiGetSubKSP()`, `PCBJacobiSetTotalBlocks()`,
498: `PCBJacobiSetLocalBlocks()`, `PCSetModifySubMatrices()`, `PCJACOBI`, `PCVPBJACOBI`, `PCPBJACOBI`
499: M*/
501: PETSC_EXTERN PetscErrorCode PCCreate_BJacobi(PC pc)
502: {
503: PetscMPIInt rank;
504: PC_BJacobi *jac;
506: PetscFunctionBegin;
507: PetscCall(PetscNew(&jac));
508: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)pc), &rank));
510: pc->ops->apply = NULL;
511: pc->ops->matapply = NULL;
512: pc->ops->applytranspose = NULL;
513: pc->ops->matapplytranspose = NULL;
514: pc->ops->setup = PCSetUp_BJacobi;
515: pc->ops->destroy = PCDestroy_BJacobi;
516: pc->ops->setfromoptions = PCSetFromOptions_BJacobi;
517: pc->ops->view = PCView_BJacobi;
518: pc->ops->applyrichardson = NULL;
520: pc->data = (void *)jac;
521: jac->n = -1;
522: jac->n_local = -1;
523: jac->first_local = rank;
524: jac->ksp = NULL;
525: jac->g_lens = NULL;
526: jac->l_lens = NULL;
527: jac->psubcomm = NULL;
529: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBJacobiGetSubKSP_C", PCBJacobiGetSubKSP_BJacobi));
530: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBJacobiSetTotalBlocks_C", PCBJacobiSetTotalBlocks_BJacobi));
531: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBJacobiGetTotalBlocks_C", PCBJacobiGetTotalBlocks_BJacobi));
532: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBJacobiSetLocalBlocks_C", PCBJacobiSetLocalBlocks_BJacobi));
533: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBJacobiGetLocalBlocks_C", PCBJacobiGetLocalBlocks_BJacobi));
534: PetscFunctionReturn(PETSC_SUCCESS);
535: }
537: /*
538: These are for a single block per processor; works for AIJ, BAIJ; Seq and MPI
539: */
540: static PetscErrorCode PCReset_BJacobi_Singleblock(PC pc)
541: {
542: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
543: PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock *)jac->data;
545: PetscFunctionBegin;
546: PetscCall(KSPReset(jac->ksp[0]));
547: PetscCall(VecDestroy(&bjac->x));
548: PetscCall(VecDestroy(&bjac->y));
549: PetscFunctionReturn(PETSC_SUCCESS);
550: }
552: static PetscErrorCode PCDestroy_BJacobi_Singleblock(PC pc)
553: {
554: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
555: PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock *)jac->data;
557: PetscFunctionBegin;
558: PetscCall(PCReset_BJacobi_Singleblock(pc));
559: PetscCall(KSPDestroy(&jac->ksp[0]));
560: PetscCall(PetscFree(jac->ksp));
561: PetscCall(PetscFree(bjac));
562: PetscCall(PCDestroy_BJacobi(pc));
563: PetscFunctionReturn(PETSC_SUCCESS);
564: }
566: static PetscErrorCode PCSetUpOnBlocks_BJacobi_Singleblock(PC pc)
567: {
568: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
569: KSP subksp = jac->ksp[0];
570: KSPConvergedReason reason;
572: PetscFunctionBegin;
573: PetscCall(KSPSetUp(subksp));
574: PetscCall(KSPGetConvergedReason(subksp, &reason));
575: if (reason == KSP_DIVERGED_PC_FAILED) pc->failedreason = PC_SUBPC_ERROR;
576: PetscFunctionReturn(PETSC_SUCCESS);
577: }
579: static PetscErrorCode PCApply_BJacobi_Singleblock(PC pc, Vec x, Vec y)
580: {
581: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
582: PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock *)jac->data;
584: PetscFunctionBegin;
585: PetscCall(VecGetLocalVectorRead(x, bjac->x));
586: PetscCall(VecGetLocalVector(y, bjac->y));
587: /* Since the inner KSP matrix may point directly to the diagonal block of an MPI matrix the inner
588: matrix may change even if the outer KSP/PC has not updated the preconditioner, this will trigger a rebuild
589: of the inner preconditioner automatically unless we pass down the outer preconditioners reuse flag.*/
590: PetscCall(KSPSetReusePreconditioner(jac->ksp[0], pc->reusepreconditioner));
591: PetscCall(PetscLogEventBegin(PC_ApplyOnBlocks, jac->ksp[0], bjac->x, bjac->y, 0));
592: PetscCall(KSPSolve(jac->ksp[0], bjac->x, bjac->y));
593: PetscCall(KSPCheckSolve(jac->ksp[0], pc, bjac->y));
594: PetscCall(PetscLogEventEnd(PC_ApplyOnBlocks, jac->ksp[0], bjac->x, bjac->y, 0));
595: PetscCall(VecRestoreLocalVectorRead(x, bjac->x));
596: PetscCall(VecRestoreLocalVector(y, bjac->y));
597: PetscFunctionReturn(PETSC_SUCCESS);
598: }
600: static PetscErrorCode PCMatApply_BJacobi_Singleblock_Private(PC pc, Mat X, Mat Y, PetscBool transpose)
601: {
602: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
603: Mat sX, sY;
605: PetscFunctionBegin;
606: /* Since the inner KSP matrix may point directly to the diagonal block of an MPI matrix the inner
607: matrix may change even if the outer KSP/PC has not updated the preconditioner, this will trigger a rebuild
608: of the inner preconditioner automatically unless we pass down the outer preconditioners reuse flag.*/
609: PetscCall(KSPSetReusePreconditioner(jac->ksp[0], pc->reusepreconditioner));
610: PetscCall(MatDenseGetLocalMatrix(X, &sX));
611: PetscCall(MatDenseGetLocalMatrix(Y, &sY));
612: if (!transpose) {
613: PetscCall(PetscLogEventBegin(PC_ApplyOnBlocks, jac->ksp[0], sX, sY, 0));
614: PetscCall(KSPMatSolve(jac->ksp[0], sX, sY));
615: PetscCall(PetscLogEventEnd(PC_ApplyOnBlocks, jac->ksp[0], sX, sY, 0));
616: } else {
617: PetscCall(PetscLogEventBegin(PC_ApplyTransposeOnBlocks, jac->ksp[0], sX, sY, 0));
618: PetscCall(KSPMatSolveTranspose(jac->ksp[0], sX, sY));
619: PetscCall(PetscLogEventEnd(PC_ApplyTransposeOnBlocks, jac->ksp[0], sX, sY, 0));
620: }
621: PetscCall(KSPCheckMatSolve(jac->ksp[0], pc, Y));
622: PetscFunctionReturn(PETSC_SUCCESS);
623: }
625: static PetscErrorCode PCMatApply_BJacobi_Singleblock(PC pc, Mat X, Mat Y)
626: {
627: PetscFunctionBegin;
628: PetscCall(PCMatApply_BJacobi_Singleblock_Private(pc, X, Y, PETSC_FALSE));
629: PetscFunctionReturn(PETSC_SUCCESS);
630: }
632: static PetscErrorCode PCMatApplyTranspose_BJacobi_Singleblock(PC pc, Mat X, Mat Y)
633: {
634: PetscFunctionBegin;
635: PetscCall(PCMatApply_BJacobi_Singleblock_Private(pc, X, Y, PETSC_TRUE));
636: PetscFunctionReturn(PETSC_SUCCESS);
637: }
639: static PetscErrorCode PCApplySymmetricLeft_BJacobi_Singleblock(PC pc, Vec x, Vec y)
640: {
641: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
642: PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock *)jac->data;
643: PetscScalar *y_array;
644: const PetscScalar *x_array;
645: PC subpc;
647: PetscFunctionBegin;
648: /*
649: The VecPlaceArray() is to avoid having to copy the
650: y vector into the bjac->x vector. The reason for
651: the bjac->x vector is that we need a sequential vector
652: for the sequential solve.
653: */
654: PetscCall(VecGetArrayRead(x, &x_array));
655: PetscCall(VecGetArray(y, &y_array));
656: PetscCall(VecPlaceArray(bjac->x, x_array));
657: PetscCall(VecPlaceArray(bjac->y, y_array));
658: /* apply the symmetric left portion of the inner PC operator */
659: /* note this bypasses the inner KSP and its options completely */
660: PetscCall(KSPGetPC(jac->ksp[0], &subpc));
661: PetscCall(PCApplySymmetricLeft(subpc, bjac->x, bjac->y));
662: PetscCall(VecResetArray(bjac->x));
663: PetscCall(VecResetArray(bjac->y));
664: PetscCall(VecRestoreArrayRead(x, &x_array));
665: PetscCall(VecRestoreArray(y, &y_array));
666: PetscFunctionReturn(PETSC_SUCCESS);
667: }
669: static PetscErrorCode PCApplySymmetricRight_BJacobi_Singleblock(PC pc, Vec x, Vec y)
670: {
671: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
672: PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock *)jac->data;
673: PetscScalar *y_array;
674: const PetscScalar *x_array;
675: PC subpc;
677: PetscFunctionBegin;
678: /*
679: The VecPlaceArray() is to avoid having to copy the
680: y vector into the bjac->x vector. The reason for
681: the bjac->x vector is that we need a sequential vector
682: for the sequential solve.
683: */
684: PetscCall(VecGetArrayRead(x, &x_array));
685: PetscCall(VecGetArray(y, &y_array));
686: PetscCall(VecPlaceArray(bjac->x, x_array));
687: PetscCall(VecPlaceArray(bjac->y, y_array));
689: /* apply the symmetric right portion of the inner PC operator */
690: /* note this bypasses the inner KSP and its options completely */
692: PetscCall(KSPGetPC(jac->ksp[0], &subpc));
693: PetscCall(PCApplySymmetricRight(subpc, bjac->x, bjac->y));
695: PetscCall(VecResetArray(bjac->x));
696: PetscCall(VecResetArray(bjac->y));
697: PetscCall(VecRestoreArrayRead(x, &x_array));
698: PetscCall(VecRestoreArray(y, &y_array));
699: PetscFunctionReturn(PETSC_SUCCESS);
700: }
702: static PetscErrorCode PCApplyTranspose_BJacobi_Singleblock(PC pc, Vec x, Vec y)
703: {
704: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
705: PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock *)jac->data;
706: PetscScalar *y_array;
707: const PetscScalar *x_array;
709: PetscFunctionBegin;
710: /*
711: The VecPlaceArray() is to avoid having to copy the
712: y vector into the bjac->x vector. The reason for
713: the bjac->x vector is that we need a sequential vector
714: for the sequential solve.
715: */
716: PetscCall(VecGetArrayRead(x, &x_array));
717: PetscCall(VecGetArray(y, &y_array));
718: PetscCall(VecPlaceArray(bjac->x, x_array));
719: PetscCall(VecPlaceArray(bjac->y, y_array));
720: PetscCall(PetscLogEventBegin(PC_ApplyTransposeOnBlocks, jac->ksp[0], bjac->x, bjac->y, 0));
721: PetscCall(KSPSolveTranspose(jac->ksp[0], bjac->x, bjac->y));
722: PetscCall(KSPCheckSolve(jac->ksp[0], pc, bjac->y));
723: PetscCall(PetscLogEventEnd(PC_ApplyTransposeOnBlocks, jac->ksp[0], bjac->x, bjac->y, 0));
724: PetscCall(VecResetArray(bjac->x));
725: PetscCall(VecResetArray(bjac->y));
726: PetscCall(VecRestoreArrayRead(x, &x_array));
727: PetscCall(VecRestoreArray(y, &y_array));
728: PetscFunctionReturn(PETSC_SUCCESS);
729: }
731: static PetscErrorCode PCSetUp_BJacobi_Singleblock(PC pc, Mat mat, Mat pmat)
732: {
733: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
734: PetscInt m;
735: KSP ksp;
736: PC_BJacobi_Singleblock *bjac;
737: PetscBool wasSetup = PETSC_TRUE;
738: VecType vectype;
739: const char *prefix;
741: PetscFunctionBegin;
742: if (!pc->setupcalled) {
743: if (!jac->ksp) {
744: PetscInt nestlevel;
746: wasSetup = PETSC_FALSE;
748: PetscCall(KSPCreate(PETSC_COMM_SELF, &ksp));
749: PetscCall(PCGetKSPNestLevel(pc, &nestlevel));
750: PetscCall(KSPSetNestLevel(ksp, nestlevel + 1));
751: PetscCall(KSPSetErrorIfNotConverged(ksp, pc->erroriffailure));
752: PetscCall(PetscObjectIncrementTabLevel((PetscObject)ksp, (PetscObject)pc, 1));
753: PetscCall(KSPSetType(ksp, KSPPREONLY));
754: PetscCall(PCGetOptionsPrefix(pc, &prefix));
755: PetscCall(KSPSetOptionsPrefix(ksp, prefix));
756: PetscCall(KSPAppendOptionsPrefix(ksp, "sub_"));
758: pc->ops->reset = PCReset_BJacobi_Singleblock;
759: pc->ops->destroy = PCDestroy_BJacobi_Singleblock;
760: pc->ops->apply = PCApply_BJacobi_Singleblock;
761: pc->ops->matapply = PCMatApply_BJacobi_Singleblock;
762: pc->ops->matapplytranspose = PCMatApplyTranspose_BJacobi_Singleblock;
763: pc->ops->applysymmetricleft = PCApplySymmetricLeft_BJacobi_Singleblock;
764: pc->ops->applysymmetricright = PCApplySymmetricRight_BJacobi_Singleblock;
765: pc->ops->applytranspose = PCApplyTranspose_BJacobi_Singleblock;
766: pc->ops->setuponblocks = PCSetUpOnBlocks_BJacobi_Singleblock;
768: PetscCall(PetscMalloc1(1, &jac->ksp));
769: jac->ksp[0] = ksp;
771: PetscCall(PetscNew(&bjac));
772: jac->data = (void *)bjac;
773: } else {
774: ksp = jac->ksp[0];
775: bjac = (PC_BJacobi_Singleblock *)jac->data;
776: }
778: /*
779: The reason we need to generate these vectors is to serve
780: as the right-hand side and solution vector for the solve on the
781: block. We do not need to allocate space for the vectors since
782: that is provided via VecPlaceArray() just before the call to
783: KSPSolve() on the block.
784: */
785: PetscCall(MatGetSize(pmat, &m, &m));
786: PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, m, NULL, &bjac->x));
787: PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, m, NULL, &bjac->y));
788: PetscCall(MatGetVecType(pmat, &vectype));
789: PetscCall(VecSetType(bjac->x, vectype));
790: PetscCall(VecSetType(bjac->y, vectype));
791: } else {
792: ksp = jac->ksp[0];
793: bjac = (PC_BJacobi_Singleblock *)jac->data;
794: }
795: PetscCall(KSPGetOptionsPrefix(ksp, &prefix));
796: if (pc->useAmat) {
797: PetscCall(KSPSetOperators(ksp, mat, pmat));
798: PetscCall(MatSetOptionsPrefix(mat, prefix));
799: } else {
800: PetscCall(KSPSetOperators(ksp, pmat, pmat));
801: }
802: PetscCall(MatSetOptionsPrefix(pmat, prefix));
803: if (!wasSetup && pc->setfromoptionscalled) {
804: /* If PCSetFromOptions_BJacobi is called later, KSPSetFromOptions will be called at that time. */
805: PetscCall(KSPSetFromOptions(ksp));
806: }
807: PetscFunctionReturn(PETSC_SUCCESS);
808: }
810: static PetscErrorCode PCReset_BJacobi_Multiblock(PC pc)
811: {
812: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
813: PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock *)jac->data;
814: PetscInt i;
816: PetscFunctionBegin;
817: if (bjac && bjac->pmat) {
818: PetscCall(MatDestroyMatrices(jac->n_local, &bjac->pmat));
819: if (pc->useAmat) PetscCall(MatDestroyMatrices(jac->n_local, &bjac->mat));
820: }
822: for (i = 0; i < jac->n_local; i++) {
823: PetscCall(KSPReset(jac->ksp[i]));
824: if (bjac && bjac->x) {
825: PetscCall(VecDestroy(&bjac->x[i]));
826: PetscCall(VecDestroy(&bjac->y[i]));
827: PetscCall(ISDestroy(&bjac->is[i]));
828: }
829: }
830: PetscCall(PetscFree(jac->l_lens));
831: PetscCall(PetscFree(jac->g_lens));
832: PetscFunctionReturn(PETSC_SUCCESS);
833: }
835: static PetscErrorCode PCDestroy_BJacobi_Multiblock(PC pc)
836: {
837: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
838: PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock *)jac->data;
839: PetscInt i;
841: PetscFunctionBegin;
842: PetscCall(PCReset_BJacobi_Multiblock(pc));
843: if (bjac) {
844: PetscCall(PetscFree2(bjac->x, bjac->y));
845: PetscCall(PetscFree(bjac->starts));
846: PetscCall(PetscFree(bjac->is));
847: }
848: PetscCall(PetscFree(jac->data));
849: for (i = 0; i < jac->n_local; i++) PetscCall(KSPDestroy(&jac->ksp[i]));
850: PetscCall(PetscFree(jac->ksp));
851: PetscCall(PCDestroy_BJacobi(pc));
852: PetscFunctionReturn(PETSC_SUCCESS);
853: }
855: static PetscErrorCode PCSetUpOnBlocks_BJacobi_Multiblock(PC pc)
856: {
857: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
858: PetscInt i, n_local = jac->n_local;
859: KSPConvergedReason reason;
861: PetscFunctionBegin;
862: for (i = 0; i < n_local; i++) {
863: PetscCall(KSPSetUp(jac->ksp[i]));
864: PetscCall(KSPGetConvergedReason(jac->ksp[i], &reason));
865: if (reason == KSP_DIVERGED_PC_FAILED) pc->failedreason = PC_SUBPC_ERROR;
866: }
867: PetscFunctionReturn(PETSC_SUCCESS);
868: }
870: static PetscErrorCode PCApply_BJacobi_Multiblock(PC pc, Vec x, Vec y)
871: {
872: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
873: PetscInt i, n_local = jac->n_local;
874: PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock *)jac->data;
875: PetscScalar *yin;
876: const PetscScalar *xin;
878: PetscFunctionBegin;
879: PetscCall(VecGetArrayRead(x, &xin));
880: PetscCall(VecGetArray(y, &yin));
881: for (i = 0; i < n_local; i++) {
882: /*
883: To avoid copying the subvector from x into a workspace we instead
884: make the workspace vector array point to the subpart of the array of
885: the global vector.
886: */
887: PetscCall(VecPlaceArray(bjac->x[i], xin + bjac->starts[i]));
888: PetscCall(VecPlaceArray(bjac->y[i], yin + bjac->starts[i]));
890: PetscCall(PetscLogEventBegin(PC_ApplyOnBlocks, jac->ksp[i], bjac->x[i], bjac->y[i], 0));
891: PetscCall(KSPSolve(jac->ksp[i], bjac->x[i], bjac->y[i]));
892: PetscCall(KSPCheckSolve(jac->ksp[i], pc, bjac->y[i]));
893: PetscCall(PetscLogEventEnd(PC_ApplyOnBlocks, jac->ksp[i], bjac->x[i], bjac->y[i], 0));
895: PetscCall(VecResetArray(bjac->x[i]));
896: PetscCall(VecResetArray(bjac->y[i]));
897: }
898: PetscCall(VecRestoreArrayRead(x, &xin));
899: PetscCall(VecRestoreArray(y, &yin));
900: PetscFunctionReturn(PETSC_SUCCESS);
901: }
903: static PetscErrorCode PCApplySymmetricLeft_BJacobi_Multiblock(PC pc, Vec x, Vec y)
904: {
905: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
906: PetscInt i, n_local = jac->n_local;
907: PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock *)jac->data;
908: PetscScalar *yin;
909: const PetscScalar *xin;
910: PC subpc;
912: PetscFunctionBegin;
913: PetscCall(VecGetArrayRead(x, &xin));
914: PetscCall(VecGetArray(y, &yin));
915: for (i = 0; i < n_local; i++) {
916: /*
917: To avoid copying the subvector from x into a workspace we instead
918: make the workspace vector array point to the subpart of the array of
919: the global vector.
920: */
921: PetscCall(VecPlaceArray(bjac->x[i], xin + bjac->starts[i]));
922: PetscCall(VecPlaceArray(bjac->y[i], yin + bjac->starts[i]));
924: PetscCall(PetscLogEventBegin(PC_ApplyOnBlocks, jac->ksp[i], bjac->x[i], bjac->y[i], 0));
925: /* apply the symmetric left portion of the inner PC operator */
926: /* note this bypasses the inner KSP and its options completely */
927: PetscCall(KSPGetPC(jac->ksp[i], &subpc));
928: PetscCall(PCApplySymmetricLeft(subpc, bjac->x[i], bjac->y[i]));
929: PetscCall(PetscLogEventEnd(PC_ApplyOnBlocks, jac->ksp[i], bjac->x[i], bjac->y[i], 0));
931: PetscCall(VecResetArray(bjac->x[i]));
932: PetscCall(VecResetArray(bjac->y[i]));
933: }
934: PetscCall(VecRestoreArrayRead(x, &xin));
935: PetscCall(VecRestoreArray(y, &yin));
936: PetscFunctionReturn(PETSC_SUCCESS);
937: }
939: static PetscErrorCode PCApplySymmetricRight_BJacobi_Multiblock(PC pc, Vec x, Vec y)
940: {
941: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
942: PetscInt i, n_local = jac->n_local;
943: PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock *)jac->data;
944: PetscScalar *yin;
945: const PetscScalar *xin;
946: PC subpc;
948: PetscFunctionBegin;
949: PetscCall(VecGetArrayRead(x, &xin));
950: PetscCall(VecGetArray(y, &yin));
951: for (i = 0; i < n_local; i++) {
952: /*
953: To avoid copying the subvector from x into a workspace we instead
954: make the workspace vector array point to the subpart of the array of
955: the global vector.
956: */
957: PetscCall(VecPlaceArray(bjac->x[i], xin + bjac->starts[i]));
958: PetscCall(VecPlaceArray(bjac->y[i], yin + bjac->starts[i]));
960: PetscCall(PetscLogEventBegin(PC_ApplyOnBlocks, jac->ksp[i], bjac->x[i], bjac->y[i], 0));
961: /* apply the symmetric left portion of the inner PC operator */
962: /* note this bypasses the inner KSP and its options completely */
963: PetscCall(KSPGetPC(jac->ksp[i], &subpc));
964: PetscCall(PCApplySymmetricRight(subpc, bjac->x[i], bjac->y[i]));
965: PetscCall(PetscLogEventEnd(PC_ApplyOnBlocks, jac->ksp[i], bjac->x[i], bjac->y[i], 0));
967: PetscCall(VecResetArray(bjac->x[i]));
968: PetscCall(VecResetArray(bjac->y[i]));
969: }
970: PetscCall(VecRestoreArrayRead(x, &xin));
971: PetscCall(VecRestoreArray(y, &yin));
972: PetscFunctionReturn(PETSC_SUCCESS);
973: }
975: static PetscErrorCode PCApplyTranspose_BJacobi_Multiblock(PC pc, Vec x, Vec y)
976: {
977: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
978: PetscInt i, n_local = jac->n_local;
979: PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock *)jac->data;
980: PetscScalar *yin;
981: const PetscScalar *xin;
983: PetscFunctionBegin;
984: PetscCall(VecGetArrayRead(x, &xin));
985: PetscCall(VecGetArray(y, &yin));
986: for (i = 0; i < n_local; i++) {
987: /*
988: To avoid copying the subvector from x into a workspace we instead
989: make the workspace vector array point to the subpart of the array of
990: the global vector.
991: */
992: PetscCall(VecPlaceArray(bjac->x[i], xin + bjac->starts[i]));
993: PetscCall(VecPlaceArray(bjac->y[i], yin + bjac->starts[i]));
995: PetscCall(PetscLogEventBegin(PC_ApplyTransposeOnBlocks, jac->ksp[i], bjac->x[i], bjac->y[i], 0));
996: PetscCall(KSPSolveTranspose(jac->ksp[i], bjac->x[i], bjac->y[i]));
997: PetscCall(KSPCheckSolve(jac->ksp[i], pc, bjac->y[i]));
998: PetscCall(PetscLogEventEnd(PC_ApplyTransposeOnBlocks, jac->ksp[i], bjac->x[i], bjac->y[i], 0));
1000: PetscCall(VecResetArray(bjac->x[i]));
1001: PetscCall(VecResetArray(bjac->y[i]));
1002: }
1003: PetscCall(VecRestoreArrayRead(x, &xin));
1004: PetscCall(VecRestoreArray(y, &yin));
1005: PetscFunctionReturn(PETSC_SUCCESS);
1006: }
1008: static PetscErrorCode PCSetUp_BJacobi_Multiblock(PC pc, Mat mat, Mat pmat)
1009: {
1010: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
1011: PetscInt m, n_local, N, M, start, i;
1012: const char *prefix;
1013: KSP ksp;
1014: Vec x, y;
1015: PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock *)jac->data;
1016: PC subpc;
1017: IS is;
1018: MatReuse scall;
1019: VecType vectype;
1020: MatNullSpace *nullsp_mat = NULL, *nullsp_pmat = NULL;
1022: PetscFunctionBegin;
1023: PetscCall(MatGetLocalSize(pc->pmat, &M, &N));
1025: n_local = jac->n_local;
1027: if (pc->useAmat) {
1028: PetscBool same;
1029: PetscCall(PetscObjectTypeCompare((PetscObject)mat, ((PetscObject)pmat)->type_name, &same));
1030: PetscCheck(same, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "Matrices not of same type");
1031: }
1033: if (!pc->setupcalled) {
1034: PetscInt nestlevel;
1036: scall = MAT_INITIAL_MATRIX;
1038: if (!jac->ksp) {
1039: pc->ops->reset = PCReset_BJacobi_Multiblock;
1040: pc->ops->destroy = PCDestroy_BJacobi_Multiblock;
1041: pc->ops->apply = PCApply_BJacobi_Multiblock;
1042: pc->ops->matapply = NULL;
1043: pc->ops->matapplytranspose = NULL;
1044: pc->ops->applysymmetricleft = PCApplySymmetricLeft_BJacobi_Multiblock;
1045: pc->ops->applysymmetricright = PCApplySymmetricRight_BJacobi_Multiblock;
1046: pc->ops->applytranspose = PCApplyTranspose_BJacobi_Multiblock;
1047: pc->ops->setuponblocks = PCSetUpOnBlocks_BJacobi_Multiblock;
1049: PetscCall(PetscNew(&bjac));
1050: PetscCall(PetscMalloc1(n_local, &jac->ksp));
1051: PetscCall(PetscMalloc2(n_local, &bjac->x, n_local, &bjac->y));
1052: PetscCall(PetscMalloc1(n_local, &bjac->starts));
1054: jac->data = (void *)bjac;
1055: PetscCall(PetscMalloc1(n_local, &bjac->is));
1057: for (i = 0; i < n_local; i++) {
1058: PetscCall(KSPCreate(PETSC_COMM_SELF, &ksp));
1059: PetscCall(PCGetKSPNestLevel(pc, &nestlevel));
1060: PetscCall(KSPSetNestLevel(ksp, nestlevel + 1));
1061: PetscCall(KSPSetErrorIfNotConverged(ksp, pc->erroriffailure));
1062: PetscCall(PetscObjectIncrementTabLevel((PetscObject)ksp, (PetscObject)pc, 1));
1063: PetscCall(KSPSetType(ksp, KSPPREONLY));
1064: PetscCall(KSPGetPC(ksp, &subpc));
1065: PetscCall(PCGetOptionsPrefix(pc, &prefix));
1066: PetscCall(KSPSetOptionsPrefix(ksp, prefix));
1067: PetscCall(KSPAppendOptionsPrefix(ksp, "sub_"));
1069: jac->ksp[i] = ksp;
1070: }
1071: } else {
1072: bjac = (PC_BJacobi_Multiblock *)jac->data;
1073: }
1075: start = 0;
1076: PetscCall(MatGetVecType(pmat, &vectype));
1077: for (i = 0; i < n_local; i++) {
1078: m = jac->l_lens[i];
1079: /*
1080: The reason we need to generate these vectors is to serve
1081: as the right-hand side and solution vector for the solve on the
1082: block. We do not need to allocate space for the vectors since
1083: that is provided via VecPlaceArray() just before the call to
1084: KSPSolve() on the block.
1086: */
1087: PetscCall(VecCreateSeq(PETSC_COMM_SELF, m, &x));
1088: PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, m, NULL, &y));
1089: PetscCall(VecSetType(x, vectype));
1090: PetscCall(VecSetType(y, vectype));
1092: bjac->x[i] = x;
1093: bjac->y[i] = y;
1094: bjac->starts[i] = start;
1096: PetscCall(ISCreateStride(PETSC_COMM_SELF, m, start, 1, &is));
1097: bjac->is[i] = is;
1099: start += m;
1100: }
1101: } else {
1102: bjac = (PC_BJacobi_Multiblock *)jac->data;
1103: /*
1104: Destroy the blocks from the previous iteration
1105: */
1106: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
1107: PetscCall(MatGetNullSpaces(n_local, bjac->pmat, &nullsp_pmat));
1108: PetscCall(MatDestroyMatrices(n_local, &bjac->pmat));
1109: if (pc->useAmat) {
1110: PetscCall(MatGetNullSpaces(n_local, bjac->mat, &nullsp_mat));
1111: PetscCall(MatDestroyMatrices(n_local, &bjac->mat));
1112: }
1113: scall = MAT_INITIAL_MATRIX;
1114: } else scall = MAT_REUSE_MATRIX;
1115: }
1117: PetscCall(MatCreateSubMatrices(pmat, n_local, bjac->is, bjac->is, scall, &bjac->pmat));
1118: if (nullsp_pmat) PetscCall(MatRestoreNullSpaces(n_local, bjac->pmat, &nullsp_pmat));
1119: if (pc->useAmat) {
1120: PetscCall(MatCreateSubMatrices(mat, n_local, bjac->is, bjac->is, scall, &bjac->mat));
1121: if (nullsp_mat) PetscCall(MatRestoreNullSpaces(n_local, bjac->mat, &nullsp_mat));
1122: }
1123: /* Return control to the user so that the submatrices can be modified (e.g., to apply
1124: different boundary conditions for the submatrices than for the global problem) */
1125: PetscCall(PCModifySubMatrices(pc, n_local, bjac->is, bjac->is, bjac->pmat, pc->modifysubmatricesP));
1127: for (i = 0; i < n_local; i++) {
1128: PetscCall(KSPGetOptionsPrefix(jac->ksp[i], &prefix));
1129: if (pc->useAmat) {
1130: PetscCall(KSPSetOperators(jac->ksp[i], bjac->mat[i], bjac->pmat[i]));
1131: PetscCall(MatSetOptionsPrefix(bjac->mat[i], prefix));
1132: } else {
1133: PetscCall(KSPSetOperators(jac->ksp[i], bjac->pmat[i], bjac->pmat[i]));
1134: }
1135: PetscCall(MatSetOptionsPrefix(bjac->pmat[i], prefix));
1136: if (pc->setfromoptionscalled) PetscCall(KSPSetFromOptions(jac->ksp[i]));
1137: }
1138: PetscFunctionReturn(PETSC_SUCCESS);
1139: }
1141: /*
1142: These are for a single block with multiple processes
1143: */
1144: static PetscErrorCode PCSetUpOnBlocks_BJacobi_Multiproc(PC pc)
1145: {
1146: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
1147: KSP subksp = jac->ksp[0];
1148: KSPConvergedReason reason;
1150: PetscFunctionBegin;
1151: PetscCall(KSPSetUp(subksp));
1152: PetscCall(KSPGetConvergedReason(subksp, &reason));
1153: if (reason == KSP_DIVERGED_PC_FAILED) pc->failedreason = PC_SUBPC_ERROR;
1154: PetscFunctionReturn(PETSC_SUCCESS);
1155: }
1157: static PetscErrorCode PCReset_BJacobi_Multiproc(PC pc)
1158: {
1159: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
1160: PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc *)jac->data;
1162: PetscFunctionBegin;
1163: PetscCall(VecDestroy(&mpjac->ysub));
1164: PetscCall(VecDestroy(&mpjac->xsub));
1165: PetscCall(MatDestroy(&mpjac->submats));
1166: if (jac->ksp) PetscCall(KSPReset(jac->ksp[0]));
1167: PetscFunctionReturn(PETSC_SUCCESS);
1168: }
1170: static PetscErrorCode PCDestroy_BJacobi_Multiproc(PC pc)
1171: {
1172: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
1173: PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc *)jac->data;
1175: PetscFunctionBegin;
1176: PetscCall(PCReset_BJacobi_Multiproc(pc));
1177: PetscCall(KSPDestroy(&jac->ksp[0]));
1178: PetscCall(PetscFree(jac->ksp));
1179: PetscCall(PetscSubcommDestroy(&mpjac->psubcomm));
1181: PetscCall(PetscFree(mpjac));
1182: PetscCall(PCDestroy_BJacobi(pc));
1183: PetscFunctionReturn(PETSC_SUCCESS);
1184: }
1186: static PetscErrorCode PCApply_BJacobi_Multiproc(PC pc, Vec x, Vec y)
1187: {
1188: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
1189: PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc *)jac->data;
1190: PetscScalar *yarray;
1191: const PetscScalar *xarray;
1193: PetscFunctionBegin;
1194: /* place x's and y's local arrays into xsub and ysub */
1195: PetscCall(VecGetArrayRead(x, &xarray));
1196: PetscCall(VecGetArray(y, &yarray));
1197: PetscCall(VecPlaceArray(mpjac->xsub, xarray));
1198: PetscCall(VecPlaceArray(mpjac->ysub, yarray));
1200: /* apply preconditioner on each matrix block */
1201: PetscCall(PetscLogEventBegin(PC_ApplyOnBlocks, jac->ksp[0], mpjac->xsub, mpjac->ysub, 0));
1202: PetscCall(KSPSolve(jac->ksp[0], mpjac->xsub, mpjac->ysub));
1203: PetscCall(KSPCheckSolve(jac->ksp[0], pc, mpjac->ysub));
1204: PetscCall(PetscLogEventEnd(PC_ApplyOnBlocks, jac->ksp[0], mpjac->xsub, mpjac->ysub, 0));
1206: PetscCall(VecResetArray(mpjac->xsub));
1207: PetscCall(VecResetArray(mpjac->ysub));
1208: PetscCall(VecRestoreArrayRead(x, &xarray));
1209: PetscCall(VecRestoreArray(y, &yarray));
1210: PetscFunctionReturn(PETSC_SUCCESS);
1211: }
1213: static PetscErrorCode PCMatApply_BJacobi_Multiproc(PC pc, Mat X, Mat Y)
1214: {
1215: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
1216: Mat sX, sY;
1217: const PetscScalar *x;
1218: PetscScalar *y;
1219: PetscInt m, N, lda, ldb;
1221: PetscFunctionBegin;
1222: /* apply preconditioner on each matrix block */
1223: PetscCall(MatGetLocalSize(X, &m, NULL));
1224: PetscCall(MatGetSize(X, NULL, &N));
1225: PetscCall(MatDenseGetLDA(X, &lda));
1226: PetscCall(MatDenseGetLDA(Y, &ldb));
1227: PetscCall(MatDenseGetArrayRead(X, &x));
1228: PetscCall(MatDenseGetArrayWrite(Y, &y));
1229: PetscCall(MatCreateDense(PetscObjectComm((PetscObject)jac->ksp[0]), m, PETSC_DECIDE, PETSC_DECIDE, N, (PetscScalar *)x, &sX));
1230: PetscCall(MatCreateDense(PetscObjectComm((PetscObject)jac->ksp[0]), m, PETSC_DECIDE, PETSC_DECIDE, N, y, &sY));
1231: PetscCall(MatDenseSetLDA(sX, lda));
1232: PetscCall(MatDenseSetLDA(sY, ldb));
1233: PetscCall(PetscLogEventBegin(PC_ApplyOnBlocks, jac->ksp[0], X, Y, 0));
1234: PetscCall(KSPMatSolve(jac->ksp[0], sX, sY));
1235: PetscCall(KSPCheckMatSolve(jac->ksp[0], pc, sY));
1236: PetscCall(PetscLogEventEnd(PC_ApplyOnBlocks, jac->ksp[0], X, Y, 0));
1237: PetscCall(MatDestroy(&sY));
1238: PetscCall(MatDestroy(&sX));
1239: PetscCall(MatDenseRestoreArrayWrite(Y, &y));
1240: PetscCall(MatDenseRestoreArrayRead(X, &x));
1241: PetscFunctionReturn(PETSC_SUCCESS);
1242: }
1244: static PetscErrorCode PCSetUp_BJacobi_Multiproc(PC pc)
1245: {
1246: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
1247: PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc *)jac->data;
1248: PetscInt m, n;
1249: MPI_Comm comm, subcomm = 0;
1250: const char *prefix;
1251: PetscBool wasSetup = PETSC_TRUE;
1252: VecType vectype;
1254: PetscFunctionBegin;
1255: PetscCall(PetscObjectGetComm((PetscObject)pc, &comm));
1256: PetscCheck(jac->n_local <= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Only a single block in a subcommunicator is supported");
1257: jac->n_local = 1; /* currently only a single block is supported for a subcommunicator */
1258: if (!pc->setupcalled) {
1259: PetscInt nestlevel;
1261: wasSetup = PETSC_FALSE;
1262: PetscCall(PetscNew(&mpjac));
1263: jac->data = (void *)mpjac;
1265: /* initialize datastructure mpjac */
1266: if (!jac->psubcomm) {
1267: /* Create default contiguous subcommunicatiors if user does not provide them */
1268: PetscCall(PetscSubcommCreate(comm, &jac->psubcomm));
1269: PetscCall(PetscSubcommSetNumber(jac->psubcomm, jac->n));
1270: PetscCall(PetscSubcommSetType(jac->psubcomm, PETSC_SUBCOMM_CONTIGUOUS));
1271: }
1272: mpjac->psubcomm = jac->psubcomm;
1273: subcomm = PetscSubcommChild(mpjac->psubcomm);
1275: /* Get matrix blocks of pmat */
1276: PetscCall(MatGetMultiProcBlock(pc->pmat, subcomm, MAT_INITIAL_MATRIX, &mpjac->submats));
1278: /* create a new PC that processors in each subcomm have copy of */
1279: PetscCall(PetscMalloc1(1, &jac->ksp));
1280: PetscCall(KSPCreate(subcomm, &jac->ksp[0]));
1281: PetscCall(PCGetKSPNestLevel(pc, &nestlevel));
1282: PetscCall(KSPSetNestLevel(jac->ksp[0], nestlevel + 1));
1283: PetscCall(KSPSetErrorIfNotConverged(jac->ksp[0], pc->erroriffailure));
1284: PetscCall(PetscObjectIncrementTabLevel((PetscObject)jac->ksp[0], (PetscObject)pc, 1));
1285: PetscCall(KSPSetOperators(jac->ksp[0], mpjac->submats, mpjac->submats));
1286: PetscCall(KSPGetPC(jac->ksp[0], &mpjac->pc));
1288: PetscCall(PCGetOptionsPrefix(pc, &prefix));
1289: PetscCall(KSPSetOptionsPrefix(jac->ksp[0], prefix));
1290: PetscCall(KSPAppendOptionsPrefix(jac->ksp[0], "sub_"));
1291: PetscCall(KSPGetOptionsPrefix(jac->ksp[0], &prefix));
1292: PetscCall(MatSetOptionsPrefix(mpjac->submats, prefix));
1294: PetscCall(MatGetLocalSize(mpjac->submats, &m, &n));
1295: PetscCall(VecCreateMPIWithArray(subcomm, 1, n, PETSC_DECIDE, NULL, &mpjac->xsub));
1296: PetscCall(VecCreateMPIWithArray(subcomm, 1, m, PETSC_DECIDE, NULL, &mpjac->ysub));
1297: PetscCall(MatGetVecType(mpjac->submats, &vectype));
1298: PetscCall(VecSetType(mpjac->xsub, vectype));
1299: PetscCall(VecSetType(mpjac->ysub, vectype));
1301: pc->ops->setuponblocks = PCSetUpOnBlocks_BJacobi_Multiproc;
1302: pc->ops->reset = PCReset_BJacobi_Multiproc;
1303: pc->ops->destroy = PCDestroy_BJacobi_Multiproc;
1304: pc->ops->apply = PCApply_BJacobi_Multiproc;
1305: pc->ops->matapply = PCMatApply_BJacobi_Multiproc;
1306: } else { /* pc->setupcalled */
1307: subcomm = PetscSubcommChild(mpjac->psubcomm);
1308: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
1309: /* destroy old matrix blocks, then get new matrix blocks */
1310: PetscCall(MatDestroy(&mpjac->submats));
1311: PetscCall(MatGetMultiProcBlock(pc->pmat, subcomm, MAT_INITIAL_MATRIX, &mpjac->submats));
1312: } else {
1313: PetscCall(MatGetMultiProcBlock(pc->pmat, subcomm, MAT_REUSE_MATRIX, &mpjac->submats));
1314: }
1315: PetscCall(KSPSetOperators(jac->ksp[0], mpjac->submats, mpjac->submats));
1316: }
1318: if (!wasSetup && pc->setfromoptionscalled) PetscCall(KSPSetFromOptions(jac->ksp[0]));
1319: PetscFunctionReturn(PETSC_SUCCESS);
1320: }