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) {
120: PetscCall(MatGetDiagonalBlock(pc->pmat, &pmat));
121: } else pmat = mat;
122: }
124: /*
125: Setup code depends on the number of blocks
126: */
127: if (jac->n_local == 1) {
128: PetscCall(PCSetUp_BJacobi_Singleblock(pc, mat, pmat));
129: } else {
130: PetscCall(PCSetUp_BJacobi_Multiblock(pc, mat, pmat));
131: }
132: PetscFunctionReturn(PETSC_SUCCESS);
133: }
135: /* Default destroy, if it has never been setup */
136: static PetscErrorCode PCDestroy_BJacobi(PC pc)
137: {
138: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
140: PetscFunctionBegin;
141: PetscCall(PetscFree(jac->g_lens));
142: PetscCall(PetscFree(jac->l_lens));
143: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBJacobiGetSubKSP_C", NULL));
144: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBJacobiSetTotalBlocks_C", NULL));
145: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBJacobiGetTotalBlocks_C", NULL));
146: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBJacobiSetLocalBlocks_C", NULL));
147: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBJacobiGetLocalBlocks_C", NULL));
148: PetscCall(PetscFree(pc->data));
149: PetscFunctionReturn(PETSC_SUCCESS);
150: }
152: static PetscErrorCode PCSetFromOptions_BJacobi(PC pc, PetscOptionItems PetscOptionsObject)
153: {
154: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
155: PetscInt blocks, i;
156: PetscBool flg;
158: PetscFunctionBegin;
159: PetscOptionsHeadBegin(PetscOptionsObject, "Block Jacobi options");
160: PetscCall(PetscOptionsInt("-pc_bjacobi_blocks", "Total number of blocks", "PCBJacobiSetTotalBlocks", jac->n, &blocks, &flg));
161: if (flg) PetscCall(PCBJacobiSetTotalBlocks(pc, blocks, NULL));
162: PetscCall(PetscOptionsInt("-pc_bjacobi_local_blocks", "Local number of blocks", "PCBJacobiSetLocalBlocks", jac->n_local, &blocks, &flg));
163: if (flg) PetscCall(PCBJacobiSetLocalBlocks(pc, blocks, NULL));
164: if (jac->ksp) {
165: /* The sub-KSP has already been set up (e.g., PCSetUp_BJacobi_Singleblock), but KSPSetFromOptions was not called
166: * unless we had already been called. */
167: for (i = 0; i < jac->n_local; i++) PetscCall(KSPSetFromOptions(jac->ksp[i]));
168: }
169: PetscOptionsHeadEnd();
170: PetscFunctionReturn(PETSC_SUCCESS);
171: }
173: #include <petscdraw.h>
174: static PetscErrorCode PCView_BJacobi(PC pc, PetscViewer viewer)
175: {
176: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
177: PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc *)jac->data;
178: PetscMPIInt rank;
179: PetscInt i;
180: PetscBool iascii, isstring, isdraw;
181: PetscViewer sviewer;
182: PetscViewerFormat format;
183: const char *prefix;
185: PetscFunctionBegin;
186: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
187: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSTRING, &isstring));
188: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
189: if (iascii) {
190: if (pc->useAmat) PetscCall(PetscViewerASCIIPrintf(viewer, " using Amat local matrix, number of blocks = %" PetscInt_FMT "\n", jac->n));
191: PetscCall(PetscViewerASCIIPrintf(viewer, " number of blocks = %" PetscInt_FMT "\n", jac->n));
192: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)pc), &rank));
193: PetscCall(PetscViewerGetFormat(viewer, &format));
194: if (format != PETSC_VIEWER_ASCII_INFO_DETAIL) {
195: PetscCall(PetscViewerASCIIPrintf(viewer, " Local solver information for first block is in the following KSP and PC objects on rank 0:\n"));
196: PetscCall(PCGetOptionsPrefix(pc, &prefix));
197: PetscCall(PetscViewerASCIIPrintf(viewer, " Use -%sksp_view ::ascii_info_detail to display information for all blocks\n", prefix ? prefix : ""));
198: if (jac->ksp && !jac->psubcomm) {
199: PetscCall(PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &sviewer));
200: if (rank == 0) {
201: PetscCall(PetscViewerASCIIPushTab(sviewer));
202: PetscCall(KSPView(jac->ksp[0], sviewer));
203: PetscCall(PetscViewerASCIIPopTab(sviewer));
204: }
205: PetscCall(PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &sviewer));
206: /* extra call needed because of the two calls to PetscViewerASCIIPushSynchronized() in PetscViewerGetSubViewer() */
207: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
208: } else if (mpjac && jac->ksp && mpjac->psubcomm) {
209: PetscCall(PetscViewerGetSubViewer(viewer, mpjac->psubcomm->child, &sviewer));
210: if (!mpjac->psubcomm->color) {
211: PetscCall(PetscViewerASCIIPushTab(sviewer));
212: PetscCall(KSPView(*jac->ksp, sviewer));
213: PetscCall(PetscViewerASCIIPopTab(sviewer));
214: }
215: PetscCall(PetscViewerRestoreSubViewer(viewer, mpjac->psubcomm->child, &sviewer));
216: /* extra call needed because of the two calls to PetscViewerASCIIPushSynchronized() in PetscViewerGetSubViewer() */
217: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
218: }
219: } else {
220: PetscInt n_global;
221: PetscCallMPI(MPIU_Allreduce(&jac->n_local, &n_global, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)pc)));
222: PetscCall(PetscViewerASCIIPushSynchronized(viewer));
223: PetscCall(PetscViewerASCIIPrintf(viewer, " Local solver information for each block is in the following KSP and PC objects:\n"));
224: PetscCall(PetscViewerASCIIPushTab(viewer));
225: PetscCall(PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &sviewer));
226: PetscCall(PetscViewerASCIIPrintf(sviewer, "[%d] number of local blocks = %" PetscInt_FMT ", first local block number = %" PetscInt_FMT "\n", rank, jac->n_local, jac->first_local));
227: for (i = 0; i < jac->n_local; i++) {
228: PetscCall(PetscViewerASCIIPrintf(sviewer, "[%d] local block number %" PetscInt_FMT "\n", rank, i));
229: PetscCall(KSPView(jac->ksp[i], sviewer));
230: PetscCall(PetscViewerASCIIPrintf(sviewer, "- - - - - - - - - - - - - - - - - -\n"));
231: }
232: PetscCall(PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &sviewer));
233: PetscCall(PetscViewerASCIIPopTab(viewer));
234: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
235: }
236: } else if (isstring) {
237: PetscCall(PetscViewerStringSPrintf(viewer, " blks=%" PetscInt_FMT, jac->n));
238: PetscCall(PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &sviewer));
239: if (jac->ksp) PetscCall(KSPView(jac->ksp[0], sviewer));
240: PetscCall(PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &sviewer));
241: } else if (isdraw) {
242: PetscDraw draw;
243: char str[25];
244: PetscReal x, y, bottom, h;
246: PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
247: PetscCall(PetscDrawGetCurrentPoint(draw, &x, &y));
248: PetscCall(PetscSNPrintf(str, 25, "Number blocks %" PetscInt_FMT, jac->n));
249: PetscCall(PetscDrawStringBoxed(draw, x, y, PETSC_DRAW_RED, PETSC_DRAW_BLACK, str, NULL, &h));
250: bottom = y - h;
251: PetscCall(PetscDrawPushCurrentPoint(draw, x, bottom));
252: /* warning the communicator on viewer is different then on ksp in parallel */
253: if (jac->ksp) PetscCall(KSPView(jac->ksp[0], viewer));
254: PetscCall(PetscDrawPopCurrentPoint(draw));
255: }
256: PetscFunctionReturn(PETSC_SUCCESS);
257: }
259: static PetscErrorCode PCBJacobiGetSubKSP_BJacobi(PC pc, PetscInt *n_local, PetscInt *first_local, KSP **ksp)
260: {
261: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
263: PetscFunctionBegin;
264: PetscCheck(pc->setupcalled, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Must call KSPSetUp() or PCSetUp() first");
266: if (n_local) *n_local = jac->n_local;
267: if (first_local) *first_local = jac->first_local;
268: if (ksp) *ksp = jac->ksp;
269: PetscFunctionReturn(PETSC_SUCCESS);
270: }
272: static PetscErrorCode PCBJacobiSetTotalBlocks_BJacobi(PC pc, PetscInt blocks, const PetscInt *lens)
273: {
274: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
276: PetscFunctionBegin;
277: PetscCheck(pc->setupcalled <= 0 || jac->n == blocks, PetscObjectComm((PetscObject)pc), PETSC_ERR_ORDER, "Cannot alter number of blocks after PCSetUp()/KSPSetUp() has been called");
278: jac->n = blocks;
279: if (!lens) jac->g_lens = NULL;
280: else {
281: PetscCall(PetscMalloc1(blocks, &jac->g_lens));
282: PetscCall(PetscArraycpy(jac->g_lens, lens, blocks));
283: }
284: PetscFunctionReturn(PETSC_SUCCESS);
285: }
287: static PetscErrorCode PCBJacobiGetTotalBlocks_BJacobi(PC pc, PetscInt *blocks, const PetscInt *lens[])
288: {
289: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
291: PetscFunctionBegin;
292: *blocks = jac->n;
293: if (lens) *lens = jac->g_lens;
294: PetscFunctionReturn(PETSC_SUCCESS);
295: }
297: static PetscErrorCode PCBJacobiSetLocalBlocks_BJacobi(PC pc, PetscInt blocks, const PetscInt lens[])
298: {
299: PC_BJacobi *jac;
301: PetscFunctionBegin;
302: jac = (PC_BJacobi *)pc->data;
304: jac->n_local = blocks;
305: if (!lens) jac->l_lens = NULL;
306: else {
307: PetscCall(PetscMalloc1(blocks, &jac->l_lens));
308: PetscCall(PetscArraycpy(jac->l_lens, lens, blocks));
309: }
310: PetscFunctionReturn(PETSC_SUCCESS);
311: }
313: static PetscErrorCode PCBJacobiGetLocalBlocks_BJacobi(PC pc, PetscInt *blocks, const PetscInt *lens[])
314: {
315: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
317: PetscFunctionBegin;
318: *blocks = jac->n_local;
319: if (lens) *lens = jac->l_lens;
320: PetscFunctionReturn(PETSC_SUCCESS);
321: }
323: /*@C
324: PCBJacobiGetSubKSP - Gets the local `KSP` contexts for all blocks on
325: this processor.
327: Not Collective
329: Input Parameter:
330: . pc - the preconditioner context
332: Output Parameters:
333: + n_local - the number of blocks on this processor, or NULL
334: . first_local - the global number of the first block on this processor, or NULL
335: - ksp - the array of KSP contexts
337: Level: advanced
339: Notes:
340: After `PCBJacobiGetSubKSP()` the array of `KSP` contexts is not to be freed.
342: Currently for some matrix implementations only 1 block per processor
343: is supported.
345: You must call `KSPSetUp()` or `PCSetUp()` before calling `PCBJacobiGetSubKSP()`.
347: .seealso: [](ch_ksp), `PCBJACOBI`, `PCASM`, `PCASMGetSubKSP()`
348: @*/
349: PetscErrorCode PCBJacobiGetSubKSP(PC pc, PetscInt *n_local, PetscInt *first_local, KSP *ksp[])
350: {
351: PetscFunctionBegin;
353: PetscUseMethod(pc, "PCBJacobiGetSubKSP_C", (PC, PetscInt *, PetscInt *, KSP **), (pc, n_local, first_local, ksp));
354: PetscFunctionReturn(PETSC_SUCCESS);
355: }
357: /*@
358: PCBJacobiSetTotalBlocks - Sets the global number of blocks for the block
359: Jacobi preconditioner.
361: Collective
363: Input Parameters:
364: + pc - the preconditioner context
365: . blocks - the number of blocks
366: - lens - [optional] integer array containing the size of each block
368: Options Database Key:
369: . -pc_bjacobi_blocks <blocks> - Sets the number of global blocks
371: Level: intermediate
373: Note:
374: Currently only a limited number of blocking configurations are supported.
375: All processors sharing the `PC` must call this routine with the same data.
377: .seealso: [](ch_ksp), `PCBJACOBI`, `PCSetUseAmat()`, `PCBJacobiSetLocalBlocks()`
378: @*/
379: PetscErrorCode PCBJacobiSetTotalBlocks(PC pc, PetscInt blocks, const PetscInt lens[])
380: {
381: PetscFunctionBegin;
383: PetscCheck(blocks > 0, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Must have positive blocks");
384: PetscTryMethod(pc, "PCBJacobiSetTotalBlocks_C", (PC, PetscInt, const PetscInt[]), (pc, blocks, lens));
385: PetscFunctionReturn(PETSC_SUCCESS);
386: }
388: /*@C
389: PCBJacobiGetTotalBlocks - Gets the global number of blocks for the block
390: Jacobi, `PCBJACOBI`, preconditioner.
392: Not Collective
394: Input Parameter:
395: . pc - the preconditioner context
397: Output Parameters:
398: + blocks - the number of blocks
399: - lens - integer array containing the size of each block
401: Level: intermediate
403: .seealso: [](ch_ksp), `PCBJACOBI`, `PCSetUseAmat()`, `PCBJacobiGetLocalBlocks()`
404: @*/
405: PetscErrorCode PCBJacobiGetTotalBlocks(PC pc, PetscInt *blocks, const PetscInt *lens[])
406: {
407: PetscFunctionBegin;
409: PetscAssertPointer(blocks, 2);
410: PetscUseMethod(pc, "PCBJacobiGetTotalBlocks_C", (PC, PetscInt *, const PetscInt *[]), (pc, blocks, lens));
411: PetscFunctionReturn(PETSC_SUCCESS);
412: }
414: /*@
415: PCBJacobiSetLocalBlocks - Sets the local number of blocks for the block
416: Jacobi, `PCBJACOBI`, preconditioner.
418: Not Collective
420: Input Parameters:
421: + pc - the preconditioner context
422: . blocks - the number of blocks
423: - lens - [optional] integer array containing size of each block
425: Options Database Key:
426: . -pc_bjacobi_local_blocks <blocks> - Sets the number of local blocks
428: Level: intermediate
430: Note:
431: Currently only a limited number of blocking configurations are supported.
433: .seealso: [](ch_ksp), `PCBJACOBI`, `PCSetUseAmat()`, `PCBJacobiSetTotalBlocks()`
434: @*/
435: PetscErrorCode PCBJacobiSetLocalBlocks(PC pc, PetscInt blocks, const PetscInt lens[])
436: {
437: PetscFunctionBegin;
439: PetscCheck(blocks >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must have nonegative blocks");
440: PetscTryMethod(pc, "PCBJacobiSetLocalBlocks_C", (PC, PetscInt, const PetscInt[]), (pc, blocks, lens));
441: PetscFunctionReturn(PETSC_SUCCESS);
442: }
444: /*@C
445: PCBJacobiGetLocalBlocks - Gets the local number of blocks for the block
446: Jacobi, `PCBJACOBI`, preconditioner.
448: Not Collective
450: Input Parameters:
451: + pc - the preconditioner context
452: . blocks - the number of blocks
453: - lens - [optional] integer array containing size of each block
455: Level: intermediate
457: Note:
458: Currently only a limited number of blocking configurations are supported.
460: .seealso: [](ch_ksp), `PCBJACOBI`, `PCSetUseAmat()`, `PCBJacobiGetTotalBlocks()`
461: @*/
462: PetscErrorCode PCBJacobiGetLocalBlocks(PC pc, PetscInt *blocks, const PetscInt *lens[])
463: {
464: PetscFunctionBegin;
466: PetscAssertPointer(blocks, 2);
467: PetscUseMethod(pc, "PCBJacobiGetLocalBlocks_C", (PC, PetscInt *, const PetscInt *[]), (pc, blocks, lens));
468: PetscFunctionReturn(PETSC_SUCCESS);
469: }
471: /*MC
472: PCBJACOBI - Use block Jacobi preconditioning, each block is (approximately) solved with
473: its own `KSP` object.
475: Options Database Keys:
476: + -pc_use_amat - use Amat to apply block of operator in inner Krylov method
477: - -pc_bjacobi_blocks <n> - use n total blocks
479: Level: beginner
481: Notes:
482: See `PCJACOBI` for diagonal Jacobi, `PCVPBJACOBI` for variable point block, and `PCPBJACOBI` for fixed size point block
484: Each processor can have one or more blocks, or a single block can be shared by several processes. Defaults to one block per processor.
486: To set options on the solvers for each block append -sub_ to all the `KSP` and `PC`
487: options database keys. For example, -sub_pc_type ilu -sub_pc_factor_levels 1 -sub_ksp_type preonly
489: To set the options on the solvers separate for each block call `PCBJacobiGetSubKSP()`
490: and set the options directly on the resulting `KSP` object (you can access its `PC`
491: `KSPGetPC()`)
493: For GPU-based vectors (`VECCUDA`, `VECViennaCL`) it is recommended to use exactly one block per MPI process for best
494: performance. Different block partitioning may lead to additional data transfers
495: between host and GPU that lead to degraded performance.
497: When multiple processes share a single block, each block encompasses exactly all the unknowns owned its set of processes.
499: .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCType`,
500: `PCASM`, `PCSetUseAmat()`, `PCGetUseAmat()`, `PCBJacobiGetSubKSP()`, `PCBJacobiSetTotalBlocks()`,
501: `PCBJacobiSetLocalBlocks()`, `PCSetModifySubMatrices()`, `PCJACOBI`, `PCVPBJACOBI`, `PCPBJACOBI`
502: M*/
504: PETSC_EXTERN PetscErrorCode PCCreate_BJacobi(PC pc)
505: {
506: PetscMPIInt rank;
507: PC_BJacobi *jac;
509: PetscFunctionBegin;
510: PetscCall(PetscNew(&jac));
511: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)pc), &rank));
513: pc->ops->apply = NULL;
514: pc->ops->matapply = NULL;
515: pc->ops->applytranspose = NULL;
516: pc->ops->setup = PCSetUp_BJacobi;
517: pc->ops->destroy = PCDestroy_BJacobi;
518: pc->ops->setfromoptions = PCSetFromOptions_BJacobi;
519: pc->ops->view = PCView_BJacobi;
520: pc->ops->applyrichardson = NULL;
522: pc->data = (void *)jac;
523: jac->n = -1;
524: jac->n_local = -1;
525: jac->first_local = rank;
526: jac->ksp = NULL;
527: jac->g_lens = NULL;
528: jac->l_lens = NULL;
529: jac->psubcomm = NULL;
531: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBJacobiGetSubKSP_C", PCBJacobiGetSubKSP_BJacobi));
532: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBJacobiSetTotalBlocks_C", PCBJacobiSetTotalBlocks_BJacobi));
533: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBJacobiGetTotalBlocks_C", PCBJacobiGetTotalBlocks_BJacobi));
534: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBJacobiSetLocalBlocks_C", PCBJacobiSetLocalBlocks_BJacobi));
535: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBJacobiGetLocalBlocks_C", PCBJacobiGetLocalBlocks_BJacobi));
536: PetscFunctionReturn(PETSC_SUCCESS);
537: }
539: /*
540: These are for a single block per processor; works for AIJ, BAIJ; Seq and MPI
541: */
542: static PetscErrorCode PCReset_BJacobi_Singleblock(PC pc)
543: {
544: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
545: PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock *)jac->data;
547: PetscFunctionBegin;
548: PetscCall(KSPReset(jac->ksp[0]));
549: PetscCall(VecDestroy(&bjac->x));
550: PetscCall(VecDestroy(&bjac->y));
551: PetscFunctionReturn(PETSC_SUCCESS);
552: }
554: static PetscErrorCode PCDestroy_BJacobi_Singleblock(PC pc)
555: {
556: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
557: PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock *)jac->data;
559: PetscFunctionBegin;
560: PetscCall(PCReset_BJacobi_Singleblock(pc));
561: PetscCall(KSPDestroy(&jac->ksp[0]));
562: PetscCall(PetscFree(jac->ksp));
563: PetscCall(PetscFree(bjac));
564: PetscCall(PCDestroy_BJacobi(pc));
565: PetscFunctionReturn(PETSC_SUCCESS);
566: }
568: static PetscErrorCode PCSetUpOnBlocks_BJacobi_Singleblock(PC pc)
569: {
570: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
571: KSP subksp = jac->ksp[0];
572: KSPConvergedReason reason;
574: PetscFunctionBegin;
575: PetscCall(KSPSetUp(subksp));
576: PetscCall(KSPGetConvergedReason(subksp, &reason));
577: if (reason == KSP_DIVERGED_PC_FAILED) pc->failedreason = PC_SUBPC_ERROR;
578: PetscFunctionReturn(PETSC_SUCCESS);
579: }
581: static PetscErrorCode PCApply_BJacobi_Singleblock(PC pc, Vec x, Vec y)
582: {
583: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
584: PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock *)jac->data;
586: PetscFunctionBegin;
587: PetscCall(VecGetLocalVectorRead(x, bjac->x));
588: PetscCall(VecGetLocalVector(y, bjac->y));
589: /* Since the inner KSP matrix may point directly to the diagonal block of an MPI matrix the inner
590: matrix may change even if the outer KSP/PC has not updated the preconditioner, this will trigger a rebuild
591: of the inner preconditioner automatically unless we pass down the outer preconditioners reuse flag.*/
592: PetscCall(KSPSetReusePreconditioner(jac->ksp[0], pc->reusepreconditioner));
593: PetscCall(PetscLogEventBegin(PC_ApplyOnBlocks, jac->ksp[0], bjac->x, bjac->y, 0));
594: PetscCall(KSPSolve(jac->ksp[0], bjac->x, bjac->y));
595: PetscCall(KSPCheckSolve(jac->ksp[0], pc, bjac->y));
596: PetscCall(PetscLogEventEnd(PC_ApplyOnBlocks, jac->ksp[0], bjac->x, bjac->y, 0));
597: PetscCall(VecRestoreLocalVectorRead(x, bjac->x));
598: PetscCall(VecRestoreLocalVector(y, bjac->y));
599: PetscFunctionReturn(PETSC_SUCCESS);
600: }
602: static PetscErrorCode PCMatApply_BJacobi_Singleblock(PC pc, Mat X, Mat Y)
603: {
604: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
605: Mat sX, sY;
607: PetscFunctionBegin;
608: /* Since the inner KSP matrix may point directly to the diagonal block of an MPI matrix the inner
609: matrix may change even if the outer KSP/PC has not updated the preconditioner, this will trigger a rebuild
610: of the inner preconditioner automatically unless we pass down the outer preconditioners reuse flag.*/
611: PetscCall(KSPSetReusePreconditioner(jac->ksp[0], pc->reusepreconditioner));
612: PetscCall(MatDenseGetLocalMatrix(X, &sX));
613: PetscCall(MatDenseGetLocalMatrix(Y, &sY));
614: PetscCall(KSPMatSolve(jac->ksp[0], sX, sY));
615: PetscFunctionReturn(PETSC_SUCCESS);
616: }
618: static PetscErrorCode PCApplySymmetricLeft_BJacobi_Singleblock(PC pc, Vec x, Vec y)
619: {
620: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
621: PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock *)jac->data;
622: PetscScalar *y_array;
623: const PetscScalar *x_array;
624: PC subpc;
626: PetscFunctionBegin;
627: /*
628: The VecPlaceArray() is to avoid having to copy the
629: y vector into the bjac->x vector. The reason for
630: the bjac->x vector is that we need a sequential vector
631: for the sequential solve.
632: */
633: PetscCall(VecGetArrayRead(x, &x_array));
634: PetscCall(VecGetArray(y, &y_array));
635: PetscCall(VecPlaceArray(bjac->x, x_array));
636: PetscCall(VecPlaceArray(bjac->y, y_array));
637: /* apply the symmetric left portion of the inner PC operator */
638: /* note this bypasses the inner KSP and its options completely */
639: PetscCall(KSPGetPC(jac->ksp[0], &subpc));
640: PetscCall(PCApplySymmetricLeft(subpc, bjac->x, bjac->y));
641: PetscCall(VecResetArray(bjac->x));
642: PetscCall(VecResetArray(bjac->y));
643: PetscCall(VecRestoreArrayRead(x, &x_array));
644: PetscCall(VecRestoreArray(y, &y_array));
645: PetscFunctionReturn(PETSC_SUCCESS);
646: }
648: static PetscErrorCode PCApplySymmetricRight_BJacobi_Singleblock(PC pc, Vec x, Vec y)
649: {
650: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
651: PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock *)jac->data;
652: PetscScalar *y_array;
653: const PetscScalar *x_array;
654: PC subpc;
656: PetscFunctionBegin;
657: /*
658: The VecPlaceArray() is to avoid having to copy the
659: y vector into the bjac->x vector. The reason for
660: the bjac->x vector is that we need a sequential vector
661: for the sequential solve.
662: */
663: PetscCall(VecGetArrayRead(x, &x_array));
664: PetscCall(VecGetArray(y, &y_array));
665: PetscCall(VecPlaceArray(bjac->x, x_array));
666: PetscCall(VecPlaceArray(bjac->y, y_array));
668: /* apply the symmetric right portion of the inner PC operator */
669: /* note this bypasses the inner KSP and its options completely */
671: PetscCall(KSPGetPC(jac->ksp[0], &subpc));
672: PetscCall(PCApplySymmetricRight(subpc, bjac->x, bjac->y));
674: PetscCall(VecResetArray(bjac->x));
675: PetscCall(VecResetArray(bjac->y));
676: PetscCall(VecRestoreArrayRead(x, &x_array));
677: PetscCall(VecRestoreArray(y, &y_array));
678: PetscFunctionReturn(PETSC_SUCCESS);
679: }
681: static PetscErrorCode PCApplyTranspose_BJacobi_Singleblock(PC pc, Vec x, Vec y)
682: {
683: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
684: PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock *)jac->data;
685: PetscScalar *y_array;
686: const PetscScalar *x_array;
688: PetscFunctionBegin;
689: /*
690: The VecPlaceArray() is to avoid having to copy the
691: y vector into the bjac->x vector. The reason for
692: the bjac->x vector is that we need a sequential vector
693: for the sequential solve.
694: */
695: PetscCall(VecGetArrayRead(x, &x_array));
696: PetscCall(VecGetArray(y, &y_array));
697: PetscCall(VecPlaceArray(bjac->x, x_array));
698: PetscCall(VecPlaceArray(bjac->y, y_array));
699: PetscCall(PetscLogEventBegin(PC_ApplyTransposeOnBlocks, jac->ksp[0], bjac->x, bjac->y, 0));
700: PetscCall(KSPSolveTranspose(jac->ksp[0], bjac->x, bjac->y));
701: PetscCall(KSPCheckSolve(jac->ksp[0], pc, bjac->y));
702: PetscCall(PetscLogEventEnd(PC_ApplyTransposeOnBlocks, jac->ksp[0], bjac->x, bjac->y, 0));
703: PetscCall(VecResetArray(bjac->x));
704: PetscCall(VecResetArray(bjac->y));
705: PetscCall(VecRestoreArrayRead(x, &x_array));
706: PetscCall(VecRestoreArray(y, &y_array));
707: PetscFunctionReturn(PETSC_SUCCESS);
708: }
710: static PetscErrorCode PCSetUp_BJacobi_Singleblock(PC pc, Mat mat, Mat pmat)
711: {
712: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
713: PetscInt m;
714: KSP ksp;
715: PC_BJacobi_Singleblock *bjac;
716: PetscBool wasSetup = PETSC_TRUE;
717: VecType vectype;
718: const char *prefix;
720: PetscFunctionBegin;
721: if (!pc->setupcalled) {
722: if (!jac->ksp) {
723: PetscInt nestlevel;
725: wasSetup = PETSC_FALSE;
727: PetscCall(KSPCreate(PETSC_COMM_SELF, &ksp));
728: PetscCall(PCGetKSPNestLevel(pc, &nestlevel));
729: PetscCall(KSPSetNestLevel(ksp, nestlevel + 1));
730: PetscCall(KSPSetErrorIfNotConverged(ksp, pc->erroriffailure));
731: PetscCall(PetscObjectIncrementTabLevel((PetscObject)ksp, (PetscObject)pc, 1));
732: PetscCall(KSPSetType(ksp, KSPPREONLY));
733: PetscCall(PCGetOptionsPrefix(pc, &prefix));
734: PetscCall(KSPSetOptionsPrefix(ksp, prefix));
735: PetscCall(KSPAppendOptionsPrefix(ksp, "sub_"));
737: pc->ops->reset = PCReset_BJacobi_Singleblock;
738: pc->ops->destroy = PCDestroy_BJacobi_Singleblock;
739: pc->ops->apply = PCApply_BJacobi_Singleblock;
740: pc->ops->matapply = PCMatApply_BJacobi_Singleblock;
741: pc->ops->applysymmetricleft = PCApplySymmetricLeft_BJacobi_Singleblock;
742: pc->ops->applysymmetricright = PCApplySymmetricRight_BJacobi_Singleblock;
743: pc->ops->applytranspose = PCApplyTranspose_BJacobi_Singleblock;
744: pc->ops->setuponblocks = PCSetUpOnBlocks_BJacobi_Singleblock;
746: PetscCall(PetscMalloc1(1, &jac->ksp));
747: jac->ksp[0] = ksp;
749: PetscCall(PetscNew(&bjac));
750: jac->data = (void *)bjac;
751: } else {
752: ksp = jac->ksp[0];
753: bjac = (PC_BJacobi_Singleblock *)jac->data;
754: }
756: /*
757: The reason we need to generate these vectors is to serve
758: as the right-hand side and solution vector for the solve on the
759: block. We do not need to allocate space for the vectors since
760: that is provided via VecPlaceArray() just before the call to
761: KSPSolve() on the block.
762: */
763: PetscCall(MatGetSize(pmat, &m, &m));
764: PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, m, NULL, &bjac->x));
765: PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, m, NULL, &bjac->y));
766: PetscCall(MatGetVecType(pmat, &vectype));
767: PetscCall(VecSetType(bjac->x, vectype));
768: PetscCall(VecSetType(bjac->y, vectype));
769: } else {
770: ksp = jac->ksp[0];
771: bjac = (PC_BJacobi_Singleblock *)jac->data;
772: }
773: PetscCall(KSPGetOptionsPrefix(ksp, &prefix));
774: if (pc->useAmat) {
775: PetscCall(KSPSetOperators(ksp, mat, pmat));
776: PetscCall(MatSetOptionsPrefix(mat, prefix));
777: } else {
778: PetscCall(KSPSetOperators(ksp, pmat, pmat));
779: }
780: PetscCall(MatSetOptionsPrefix(pmat, prefix));
781: if (!wasSetup && pc->setfromoptionscalled) {
782: /* If PCSetFromOptions_BJacobi is called later, KSPSetFromOptions will be called at that time. */
783: PetscCall(KSPSetFromOptions(ksp));
784: }
785: PetscFunctionReturn(PETSC_SUCCESS);
786: }
788: static PetscErrorCode PCReset_BJacobi_Multiblock(PC pc)
789: {
790: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
791: PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock *)jac->data;
792: PetscInt i;
794: PetscFunctionBegin;
795: if (bjac && bjac->pmat) {
796: PetscCall(MatDestroyMatrices(jac->n_local, &bjac->pmat));
797: if (pc->useAmat) PetscCall(MatDestroyMatrices(jac->n_local, &bjac->mat));
798: }
800: for (i = 0; i < jac->n_local; i++) {
801: PetscCall(KSPReset(jac->ksp[i]));
802: if (bjac && bjac->x) {
803: PetscCall(VecDestroy(&bjac->x[i]));
804: PetscCall(VecDestroy(&bjac->y[i]));
805: PetscCall(ISDestroy(&bjac->is[i]));
806: }
807: }
808: PetscCall(PetscFree(jac->l_lens));
809: PetscCall(PetscFree(jac->g_lens));
810: PetscFunctionReturn(PETSC_SUCCESS);
811: }
813: static PetscErrorCode PCDestroy_BJacobi_Multiblock(PC pc)
814: {
815: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
816: PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock *)jac->data;
817: PetscInt i;
819: PetscFunctionBegin;
820: PetscCall(PCReset_BJacobi_Multiblock(pc));
821: if (bjac) {
822: PetscCall(PetscFree2(bjac->x, bjac->y));
823: PetscCall(PetscFree(bjac->starts));
824: PetscCall(PetscFree(bjac->is));
825: }
826: PetscCall(PetscFree(jac->data));
827: for (i = 0; i < jac->n_local; i++) PetscCall(KSPDestroy(&jac->ksp[i]));
828: PetscCall(PetscFree(jac->ksp));
829: PetscCall(PCDestroy_BJacobi(pc));
830: PetscFunctionReturn(PETSC_SUCCESS);
831: }
833: static PetscErrorCode PCSetUpOnBlocks_BJacobi_Multiblock(PC pc)
834: {
835: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
836: PetscInt i, n_local = jac->n_local;
837: KSPConvergedReason reason;
839: PetscFunctionBegin;
840: for (i = 0; i < n_local; i++) {
841: PetscCall(KSPSetUp(jac->ksp[i]));
842: PetscCall(KSPGetConvergedReason(jac->ksp[i], &reason));
843: if (reason == KSP_DIVERGED_PC_FAILED) pc->failedreason = PC_SUBPC_ERROR;
844: }
845: PetscFunctionReturn(PETSC_SUCCESS);
846: }
848: static PetscErrorCode PCApply_BJacobi_Multiblock(PC pc, Vec x, Vec y)
849: {
850: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
851: PetscInt i, n_local = jac->n_local;
852: PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock *)jac->data;
853: PetscScalar *yin;
854: const PetscScalar *xin;
856: PetscFunctionBegin;
857: PetscCall(VecGetArrayRead(x, &xin));
858: PetscCall(VecGetArray(y, &yin));
859: for (i = 0; i < n_local; i++) {
860: /*
861: To avoid copying the subvector from x into a workspace we instead
862: make the workspace vector array point to the subpart of the array of
863: the global vector.
864: */
865: PetscCall(VecPlaceArray(bjac->x[i], xin + bjac->starts[i]));
866: PetscCall(VecPlaceArray(bjac->y[i], yin + bjac->starts[i]));
868: PetscCall(PetscLogEventBegin(PC_ApplyOnBlocks, jac->ksp[i], bjac->x[i], bjac->y[i], 0));
869: PetscCall(KSPSolve(jac->ksp[i], bjac->x[i], bjac->y[i]));
870: PetscCall(KSPCheckSolve(jac->ksp[i], pc, bjac->y[i]));
871: PetscCall(PetscLogEventEnd(PC_ApplyOnBlocks, jac->ksp[i], bjac->x[i], bjac->y[i], 0));
873: PetscCall(VecResetArray(bjac->x[i]));
874: PetscCall(VecResetArray(bjac->y[i]));
875: }
876: PetscCall(VecRestoreArrayRead(x, &xin));
877: PetscCall(VecRestoreArray(y, &yin));
878: PetscFunctionReturn(PETSC_SUCCESS);
879: }
881: static PetscErrorCode PCApplySymmetricLeft_BJacobi_Multiblock(PC pc, Vec x, Vec y)
882: {
883: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
884: PetscInt i, n_local = jac->n_local;
885: PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock *)jac->data;
886: PetscScalar *yin;
887: const PetscScalar *xin;
888: PC subpc;
890: PetscFunctionBegin;
891: PetscCall(VecGetArrayRead(x, &xin));
892: PetscCall(VecGetArray(y, &yin));
893: for (i = 0; i < n_local; i++) {
894: /*
895: To avoid copying the subvector from x into a workspace we instead
896: make the workspace vector array point to the subpart of the array of
897: the global vector.
898: */
899: PetscCall(VecPlaceArray(bjac->x[i], xin + bjac->starts[i]));
900: PetscCall(VecPlaceArray(bjac->y[i], yin + bjac->starts[i]));
902: PetscCall(PetscLogEventBegin(PC_ApplyOnBlocks, jac->ksp[i], bjac->x[i], bjac->y[i], 0));
903: /* apply the symmetric left portion of the inner PC operator */
904: /* note this bypasses the inner KSP and its options completely */
905: PetscCall(KSPGetPC(jac->ksp[i], &subpc));
906: PetscCall(PCApplySymmetricLeft(subpc, bjac->x[i], bjac->y[i]));
907: PetscCall(PetscLogEventEnd(PC_ApplyOnBlocks, jac->ksp[i], bjac->x[i], bjac->y[i], 0));
909: PetscCall(VecResetArray(bjac->x[i]));
910: PetscCall(VecResetArray(bjac->y[i]));
911: }
912: PetscCall(VecRestoreArrayRead(x, &xin));
913: PetscCall(VecRestoreArray(y, &yin));
914: PetscFunctionReturn(PETSC_SUCCESS);
915: }
917: static PetscErrorCode PCApplySymmetricRight_BJacobi_Multiblock(PC pc, Vec x, Vec y)
918: {
919: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
920: PetscInt i, n_local = jac->n_local;
921: PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock *)jac->data;
922: PetscScalar *yin;
923: const PetscScalar *xin;
924: PC subpc;
926: PetscFunctionBegin;
927: PetscCall(VecGetArrayRead(x, &xin));
928: PetscCall(VecGetArray(y, &yin));
929: for (i = 0; i < n_local; i++) {
930: /*
931: To avoid copying the subvector from x into a workspace we instead
932: make the workspace vector array point to the subpart of the array of
933: the global vector.
934: */
935: PetscCall(VecPlaceArray(bjac->x[i], xin + bjac->starts[i]));
936: PetscCall(VecPlaceArray(bjac->y[i], yin + bjac->starts[i]));
938: PetscCall(PetscLogEventBegin(PC_ApplyOnBlocks, jac->ksp[i], bjac->x[i], bjac->y[i], 0));
939: /* apply the symmetric left portion of the inner PC operator */
940: /* note this bypasses the inner KSP and its options completely */
941: PetscCall(KSPGetPC(jac->ksp[i], &subpc));
942: PetscCall(PCApplySymmetricRight(subpc, bjac->x[i], bjac->y[i]));
943: PetscCall(PetscLogEventEnd(PC_ApplyOnBlocks, jac->ksp[i], bjac->x[i], bjac->y[i], 0));
945: PetscCall(VecResetArray(bjac->x[i]));
946: PetscCall(VecResetArray(bjac->y[i]));
947: }
948: PetscCall(VecRestoreArrayRead(x, &xin));
949: PetscCall(VecRestoreArray(y, &yin));
950: PetscFunctionReturn(PETSC_SUCCESS);
951: }
953: static PetscErrorCode PCApplyTranspose_BJacobi_Multiblock(PC pc, Vec x, Vec y)
954: {
955: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
956: PetscInt i, n_local = jac->n_local;
957: PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock *)jac->data;
958: PetscScalar *yin;
959: const PetscScalar *xin;
961: PetscFunctionBegin;
962: PetscCall(VecGetArrayRead(x, &xin));
963: PetscCall(VecGetArray(y, &yin));
964: for (i = 0; i < n_local; i++) {
965: /*
966: To avoid copying the subvector from x into a workspace we instead
967: make the workspace vector array point to the subpart of the array of
968: the global vector.
969: */
970: PetscCall(VecPlaceArray(bjac->x[i], xin + bjac->starts[i]));
971: PetscCall(VecPlaceArray(bjac->y[i], yin + bjac->starts[i]));
973: PetscCall(PetscLogEventBegin(PC_ApplyTransposeOnBlocks, jac->ksp[i], bjac->x[i], bjac->y[i], 0));
974: PetscCall(KSPSolveTranspose(jac->ksp[i], bjac->x[i], bjac->y[i]));
975: PetscCall(KSPCheckSolve(jac->ksp[i], pc, bjac->y[i]));
976: PetscCall(PetscLogEventEnd(PC_ApplyTransposeOnBlocks, jac->ksp[i], bjac->x[i], bjac->y[i], 0));
978: PetscCall(VecResetArray(bjac->x[i]));
979: PetscCall(VecResetArray(bjac->y[i]));
980: }
981: PetscCall(VecRestoreArrayRead(x, &xin));
982: PetscCall(VecRestoreArray(y, &yin));
983: PetscFunctionReturn(PETSC_SUCCESS);
984: }
986: static PetscErrorCode PCSetUp_BJacobi_Multiblock(PC pc, Mat mat, Mat pmat)
987: {
988: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
989: PetscInt m, n_local, N, M, start, i;
990: const char *prefix;
991: KSP ksp;
992: Vec x, y;
993: PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock *)jac->data;
994: PC subpc;
995: IS is;
996: MatReuse scall;
997: VecType vectype;
998: MatNullSpace *nullsp_mat = NULL, *nullsp_pmat = NULL;
1000: PetscFunctionBegin;
1001: PetscCall(MatGetLocalSize(pc->pmat, &M, &N));
1003: n_local = jac->n_local;
1005: if (pc->useAmat) {
1006: PetscBool same;
1007: PetscCall(PetscObjectTypeCompare((PetscObject)mat, ((PetscObject)pmat)->type_name, &same));
1008: PetscCheck(same, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "Matrices not of same type");
1009: }
1011: if (!pc->setupcalled) {
1012: PetscInt nestlevel;
1014: scall = MAT_INITIAL_MATRIX;
1016: if (!jac->ksp) {
1017: pc->ops->reset = PCReset_BJacobi_Multiblock;
1018: pc->ops->destroy = PCDestroy_BJacobi_Multiblock;
1019: pc->ops->apply = PCApply_BJacobi_Multiblock;
1020: pc->ops->matapply = NULL;
1021: pc->ops->applysymmetricleft = PCApplySymmetricLeft_BJacobi_Multiblock;
1022: pc->ops->applysymmetricright = PCApplySymmetricRight_BJacobi_Multiblock;
1023: pc->ops->applytranspose = PCApplyTranspose_BJacobi_Multiblock;
1024: pc->ops->setuponblocks = PCSetUpOnBlocks_BJacobi_Multiblock;
1026: PetscCall(PetscNew(&bjac));
1027: PetscCall(PetscMalloc1(n_local, &jac->ksp));
1028: PetscCall(PetscMalloc2(n_local, &bjac->x, n_local, &bjac->y));
1029: PetscCall(PetscMalloc1(n_local, &bjac->starts));
1031: jac->data = (void *)bjac;
1032: PetscCall(PetscMalloc1(n_local, &bjac->is));
1034: for (i = 0; i < n_local; i++) {
1035: PetscCall(KSPCreate(PETSC_COMM_SELF, &ksp));
1036: PetscCall(PCGetKSPNestLevel(pc, &nestlevel));
1037: PetscCall(KSPSetNestLevel(ksp, nestlevel + 1));
1038: PetscCall(KSPSetErrorIfNotConverged(ksp, pc->erroriffailure));
1039: PetscCall(PetscObjectIncrementTabLevel((PetscObject)ksp, (PetscObject)pc, 1));
1040: PetscCall(KSPSetType(ksp, KSPPREONLY));
1041: PetscCall(KSPGetPC(ksp, &subpc));
1042: PetscCall(PCGetOptionsPrefix(pc, &prefix));
1043: PetscCall(KSPSetOptionsPrefix(ksp, prefix));
1044: PetscCall(KSPAppendOptionsPrefix(ksp, "sub_"));
1046: jac->ksp[i] = ksp;
1047: }
1048: } else {
1049: bjac = (PC_BJacobi_Multiblock *)jac->data;
1050: }
1052: start = 0;
1053: PetscCall(MatGetVecType(pmat, &vectype));
1054: for (i = 0; i < n_local; i++) {
1055: m = jac->l_lens[i];
1056: /*
1057: The reason we need to generate these vectors is to serve
1058: as the right-hand side and solution vector for the solve on the
1059: block. We do not need to allocate space for the vectors since
1060: that is provided via VecPlaceArray() just before the call to
1061: KSPSolve() on the block.
1063: */
1064: PetscCall(VecCreateSeq(PETSC_COMM_SELF, m, &x));
1065: PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, m, NULL, &y));
1066: PetscCall(VecSetType(x, vectype));
1067: PetscCall(VecSetType(y, vectype));
1069: bjac->x[i] = x;
1070: bjac->y[i] = y;
1071: bjac->starts[i] = start;
1073: PetscCall(ISCreateStride(PETSC_COMM_SELF, m, start, 1, &is));
1074: bjac->is[i] = is;
1076: start += m;
1077: }
1078: } else {
1079: bjac = (PC_BJacobi_Multiblock *)jac->data;
1080: /*
1081: Destroy the blocks from the previous iteration
1082: */
1083: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
1084: PetscCall(MatGetNullSpaces(n_local, bjac->pmat, &nullsp_pmat));
1085: PetscCall(MatDestroyMatrices(n_local, &bjac->pmat));
1086: if (pc->useAmat) {
1087: PetscCall(MatGetNullSpaces(n_local, bjac->mat, &nullsp_mat));
1088: PetscCall(MatDestroyMatrices(n_local, &bjac->mat));
1089: }
1090: scall = MAT_INITIAL_MATRIX;
1091: } else scall = MAT_REUSE_MATRIX;
1092: }
1094: PetscCall(MatCreateSubMatrices(pmat, n_local, bjac->is, bjac->is, scall, &bjac->pmat));
1095: if (nullsp_pmat) PetscCall(MatRestoreNullSpaces(n_local, bjac->pmat, &nullsp_pmat));
1096: if (pc->useAmat) {
1097: PetscCall(MatCreateSubMatrices(mat, n_local, bjac->is, bjac->is, scall, &bjac->mat));
1098: if (nullsp_mat) PetscCall(MatRestoreNullSpaces(n_local, bjac->mat, &nullsp_mat));
1099: }
1100: /* Return control to the user so that the submatrices can be modified (e.g., to apply
1101: different boundary conditions for the submatrices than for the global problem) */
1102: PetscCall(PCModifySubMatrices(pc, n_local, bjac->is, bjac->is, bjac->pmat, pc->modifysubmatricesP));
1104: for (i = 0; i < n_local; i++) {
1105: PetscCall(KSPGetOptionsPrefix(jac->ksp[i], &prefix));
1106: if (pc->useAmat) {
1107: PetscCall(KSPSetOperators(jac->ksp[i], bjac->mat[i], bjac->pmat[i]));
1108: PetscCall(MatSetOptionsPrefix(bjac->mat[i], prefix));
1109: } else {
1110: PetscCall(KSPSetOperators(jac->ksp[i], bjac->pmat[i], bjac->pmat[i]));
1111: }
1112: PetscCall(MatSetOptionsPrefix(bjac->pmat[i], prefix));
1113: if (pc->setfromoptionscalled) PetscCall(KSPSetFromOptions(jac->ksp[i]));
1114: }
1115: PetscFunctionReturn(PETSC_SUCCESS);
1116: }
1118: /*
1119: These are for a single block with multiple processes
1120: */
1121: static PetscErrorCode PCSetUpOnBlocks_BJacobi_Multiproc(PC pc)
1122: {
1123: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
1124: KSP subksp = jac->ksp[0];
1125: KSPConvergedReason reason;
1127: PetscFunctionBegin;
1128: PetscCall(KSPSetUp(subksp));
1129: PetscCall(KSPGetConvergedReason(subksp, &reason));
1130: if (reason == KSP_DIVERGED_PC_FAILED) pc->failedreason = PC_SUBPC_ERROR;
1131: PetscFunctionReturn(PETSC_SUCCESS);
1132: }
1134: static PetscErrorCode PCReset_BJacobi_Multiproc(PC pc)
1135: {
1136: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
1137: PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc *)jac->data;
1139: PetscFunctionBegin;
1140: PetscCall(VecDestroy(&mpjac->ysub));
1141: PetscCall(VecDestroy(&mpjac->xsub));
1142: PetscCall(MatDestroy(&mpjac->submats));
1143: if (jac->ksp) PetscCall(KSPReset(jac->ksp[0]));
1144: PetscFunctionReturn(PETSC_SUCCESS);
1145: }
1147: static PetscErrorCode PCDestroy_BJacobi_Multiproc(PC pc)
1148: {
1149: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
1150: PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc *)jac->data;
1152: PetscFunctionBegin;
1153: PetscCall(PCReset_BJacobi_Multiproc(pc));
1154: PetscCall(KSPDestroy(&jac->ksp[0]));
1155: PetscCall(PetscFree(jac->ksp));
1156: PetscCall(PetscSubcommDestroy(&mpjac->psubcomm));
1158: PetscCall(PetscFree(mpjac));
1159: PetscCall(PCDestroy_BJacobi(pc));
1160: PetscFunctionReturn(PETSC_SUCCESS);
1161: }
1163: static PetscErrorCode PCApply_BJacobi_Multiproc(PC pc, Vec x, Vec y)
1164: {
1165: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
1166: PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc *)jac->data;
1167: PetscScalar *yarray;
1168: const PetscScalar *xarray;
1169: KSPConvergedReason reason;
1171: PetscFunctionBegin;
1172: /* place x's and y's local arrays into xsub and ysub */
1173: PetscCall(VecGetArrayRead(x, &xarray));
1174: PetscCall(VecGetArray(y, &yarray));
1175: PetscCall(VecPlaceArray(mpjac->xsub, xarray));
1176: PetscCall(VecPlaceArray(mpjac->ysub, yarray));
1178: /* apply preconditioner on each matrix block */
1179: PetscCall(PetscLogEventBegin(PC_ApplyOnBlocks, jac->ksp[0], mpjac->xsub, mpjac->ysub, 0));
1180: PetscCall(KSPSolve(jac->ksp[0], mpjac->xsub, mpjac->ysub));
1181: PetscCall(KSPCheckSolve(jac->ksp[0], pc, mpjac->ysub));
1182: PetscCall(PetscLogEventEnd(PC_ApplyOnBlocks, jac->ksp[0], mpjac->xsub, mpjac->ysub, 0));
1183: PetscCall(KSPGetConvergedReason(jac->ksp[0], &reason));
1184: if (reason == KSP_DIVERGED_PC_FAILED) pc->failedreason = PC_SUBPC_ERROR;
1186: PetscCall(VecResetArray(mpjac->xsub));
1187: PetscCall(VecResetArray(mpjac->ysub));
1188: PetscCall(VecRestoreArrayRead(x, &xarray));
1189: PetscCall(VecRestoreArray(y, &yarray));
1190: PetscFunctionReturn(PETSC_SUCCESS);
1191: }
1193: static PetscErrorCode PCMatApply_BJacobi_Multiproc(PC pc, Mat X, Mat Y)
1194: {
1195: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
1196: KSPConvergedReason reason;
1197: Mat sX, sY;
1198: const PetscScalar *x;
1199: PetscScalar *y;
1200: PetscInt m, N, lda, ldb;
1202: PetscFunctionBegin;
1203: /* apply preconditioner on each matrix block */
1204: PetscCall(MatGetLocalSize(X, &m, NULL));
1205: PetscCall(MatGetSize(X, NULL, &N));
1206: PetscCall(MatDenseGetLDA(X, &lda));
1207: PetscCall(MatDenseGetLDA(Y, &ldb));
1208: PetscCall(MatDenseGetArrayRead(X, &x));
1209: PetscCall(MatDenseGetArrayWrite(Y, &y));
1210: PetscCall(MatCreateDense(PetscObjectComm((PetscObject)jac->ksp[0]), m, PETSC_DECIDE, PETSC_DECIDE, N, (PetscScalar *)x, &sX));
1211: PetscCall(MatCreateDense(PetscObjectComm((PetscObject)jac->ksp[0]), m, PETSC_DECIDE, PETSC_DECIDE, N, y, &sY));
1212: PetscCall(MatDenseSetLDA(sX, lda));
1213: PetscCall(MatDenseSetLDA(sY, ldb));
1214: PetscCall(PetscLogEventBegin(PC_ApplyOnBlocks, jac->ksp[0], X, Y, 0));
1215: PetscCall(KSPMatSolve(jac->ksp[0], sX, sY));
1216: PetscCall(KSPCheckSolve(jac->ksp[0], pc, NULL));
1217: PetscCall(PetscLogEventEnd(PC_ApplyOnBlocks, jac->ksp[0], X, Y, 0));
1218: PetscCall(MatDestroy(&sY));
1219: PetscCall(MatDestroy(&sX));
1220: PetscCall(MatDenseRestoreArrayWrite(Y, &y));
1221: PetscCall(MatDenseRestoreArrayRead(X, &x));
1222: PetscCall(KSPGetConvergedReason(jac->ksp[0], &reason));
1223: if (reason == KSP_DIVERGED_PC_FAILED) pc->failedreason = PC_SUBPC_ERROR;
1224: PetscFunctionReturn(PETSC_SUCCESS);
1225: }
1227: static PetscErrorCode PCSetUp_BJacobi_Multiproc(PC pc)
1228: {
1229: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
1230: PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc *)jac->data;
1231: PetscInt m, n;
1232: MPI_Comm comm, subcomm = 0;
1233: const char *prefix;
1234: PetscBool wasSetup = PETSC_TRUE;
1235: VecType vectype;
1237: PetscFunctionBegin;
1238: PetscCall(PetscObjectGetComm((PetscObject)pc, &comm));
1239: PetscCheck(jac->n_local <= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Only a single block in a subcommunicator is supported");
1240: jac->n_local = 1; /* currently only a single block is supported for a subcommunicator */
1241: if (!pc->setupcalled) {
1242: PetscInt nestlevel;
1244: wasSetup = PETSC_FALSE;
1245: PetscCall(PetscNew(&mpjac));
1246: jac->data = (void *)mpjac;
1248: /* initialize datastructure mpjac */
1249: if (!jac->psubcomm) {
1250: /* Create default contiguous subcommunicatiors if user does not provide them */
1251: PetscCall(PetscSubcommCreate(comm, &jac->psubcomm));
1252: PetscCall(PetscSubcommSetNumber(jac->psubcomm, jac->n));
1253: PetscCall(PetscSubcommSetType(jac->psubcomm, PETSC_SUBCOMM_CONTIGUOUS));
1254: }
1255: mpjac->psubcomm = jac->psubcomm;
1256: subcomm = PetscSubcommChild(mpjac->psubcomm);
1258: /* Get matrix blocks of pmat */
1259: PetscCall(MatGetMultiProcBlock(pc->pmat, subcomm, MAT_INITIAL_MATRIX, &mpjac->submats));
1261: /* create a new PC that processors in each subcomm have copy of */
1262: PetscCall(PetscMalloc1(1, &jac->ksp));
1263: PetscCall(KSPCreate(subcomm, &jac->ksp[0]));
1264: PetscCall(PCGetKSPNestLevel(pc, &nestlevel));
1265: PetscCall(KSPSetNestLevel(jac->ksp[0], nestlevel + 1));
1266: PetscCall(KSPSetErrorIfNotConverged(jac->ksp[0], pc->erroriffailure));
1267: PetscCall(PetscObjectIncrementTabLevel((PetscObject)jac->ksp[0], (PetscObject)pc, 1));
1268: PetscCall(KSPSetOperators(jac->ksp[0], mpjac->submats, mpjac->submats));
1269: PetscCall(KSPGetPC(jac->ksp[0], &mpjac->pc));
1271: PetscCall(PCGetOptionsPrefix(pc, &prefix));
1272: PetscCall(KSPSetOptionsPrefix(jac->ksp[0], prefix));
1273: PetscCall(KSPAppendOptionsPrefix(jac->ksp[0], "sub_"));
1274: PetscCall(KSPGetOptionsPrefix(jac->ksp[0], &prefix));
1275: PetscCall(MatSetOptionsPrefix(mpjac->submats, prefix));
1277: /* create dummy vectors xsub and ysub */
1278: PetscCall(MatGetLocalSize(mpjac->submats, &m, &n));
1279: PetscCall(VecCreateMPIWithArray(subcomm, 1, n, PETSC_DECIDE, NULL, &mpjac->xsub));
1280: PetscCall(VecCreateMPIWithArray(subcomm, 1, m, PETSC_DECIDE, NULL, &mpjac->ysub));
1281: PetscCall(MatGetVecType(mpjac->submats, &vectype));
1282: PetscCall(VecSetType(mpjac->xsub, vectype));
1283: PetscCall(VecSetType(mpjac->ysub, vectype));
1285: pc->ops->setuponblocks = PCSetUpOnBlocks_BJacobi_Multiproc;
1286: pc->ops->reset = PCReset_BJacobi_Multiproc;
1287: pc->ops->destroy = PCDestroy_BJacobi_Multiproc;
1288: pc->ops->apply = PCApply_BJacobi_Multiproc;
1289: pc->ops->matapply = PCMatApply_BJacobi_Multiproc;
1290: } else { /* pc->setupcalled */
1291: subcomm = PetscSubcommChild(mpjac->psubcomm);
1292: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
1293: /* destroy old matrix blocks, then get new matrix blocks */
1294: if (mpjac->submats) PetscCall(MatDestroy(&mpjac->submats));
1295: PetscCall(MatGetMultiProcBlock(pc->pmat, subcomm, MAT_INITIAL_MATRIX, &mpjac->submats));
1296: } else {
1297: PetscCall(MatGetMultiProcBlock(pc->pmat, subcomm, MAT_REUSE_MATRIX, &mpjac->submats));
1298: }
1299: PetscCall(KSPSetOperators(jac->ksp[0], mpjac->submats, mpjac->submats));
1300: }
1302: if (!wasSetup && pc->setfromoptionscalled) PetscCall(KSPSetFromOptions(jac->ksp[0]));
1303: PetscFunctionReturn(PETSC_SUCCESS);
1304: }