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, i;
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 (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: PetscInt i;
176: PetscBool isascii, isstring, isdraw;
177: PetscViewer sviewer;
178: PetscViewerFormat format;
179: const char *prefix;
181: PetscFunctionBegin;
182: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
183: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSTRING, &isstring));
184: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERDRAW, &isdraw));
185: if (isascii) {
186: if (pc->useAmat) PetscCall(PetscViewerASCIIPrintf(viewer, " using Amat local matrix, number of blocks = %" PetscInt_FMT "\n", jac->n));
187: PetscCall(PetscViewerASCIIPrintf(viewer, " number of blocks = %" PetscInt_FMT "\n", jac->n));
188: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)pc), &rank));
189: PetscCall(PetscViewerGetFormat(viewer, &format));
190: if (format != PETSC_VIEWER_ASCII_INFO_DETAIL) {
191: PetscCall(PetscViewerASCIIPrintf(viewer, " Local solver information for first block is in the following KSP and PC objects on rank 0:\n"));
192: PetscCall(PCGetOptionsPrefix(pc, &prefix));
193: PetscCall(PetscViewerASCIIPrintf(viewer, " Use -%sksp_view ::ascii_info_detail to display information for all blocks\n", prefix ? prefix : ""));
194: if (jac->ksp && !jac->psubcomm) {
195: PetscCall(PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &sviewer));
196: if (rank == 0) {
197: PetscCall(PetscViewerASCIIPushTab(sviewer));
198: PetscCall(KSPView(jac->ksp[0], sviewer));
199: PetscCall(PetscViewerASCIIPopTab(sviewer));
200: }
201: PetscCall(PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &sviewer));
202: /* extra call needed because of the two calls to PetscViewerASCIIPushSynchronized() in PetscViewerGetSubViewer() */
203: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
204: } else if (mpjac && jac->ksp && mpjac->psubcomm) {
205: PetscCall(PetscViewerGetSubViewer(viewer, mpjac->psubcomm->child, &sviewer));
206: if (!mpjac->psubcomm->color) {
207: PetscCall(PetscViewerASCIIPushTab(sviewer));
208: PetscCall(KSPView(*jac->ksp, sviewer));
209: PetscCall(PetscViewerASCIIPopTab(sviewer));
210: }
211: PetscCall(PetscViewerRestoreSubViewer(viewer, mpjac->psubcomm->child, &sviewer));
212: /* extra call needed because of the two calls to PetscViewerASCIIPushSynchronized() in PetscViewerGetSubViewer() */
213: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
214: }
215: } else {
216: PetscInt n_global;
217: PetscCallMPI(MPIU_Allreduce(&jac->n_local, &n_global, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)pc)));
218: PetscCall(PetscViewerASCIIPushSynchronized(viewer));
219: PetscCall(PetscViewerASCIIPrintf(viewer, " Local solver information for each block is in the following KSP and PC objects:\n"));
220: PetscCall(PetscViewerASCIIPushTab(viewer));
221: PetscCall(PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &sviewer));
222: PetscCall(PetscViewerASCIIPrintf(sviewer, "[%d] number of local blocks = %" PetscInt_FMT ", first local block number = %" PetscInt_FMT "\n", rank, jac->n_local, jac->first_local));
223: for (i = 0; i < jac->n_local; i++) {
224: PetscCall(PetscViewerASCIIPrintf(sviewer, "[%d] local block number %" PetscInt_FMT "\n", rank, i));
225: PetscCall(KSPView(jac->ksp[i], sviewer));
226: PetscCall(PetscViewerASCIIPrintf(sviewer, "- - - - - - - - - - - - - - - - - -\n"));
227: }
228: PetscCall(PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &sviewer));
229: PetscCall(PetscViewerASCIIPopTab(viewer));
230: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
231: }
232: } else if (isstring) {
233: PetscCall(PetscViewerStringSPrintf(viewer, " blks=%" PetscInt_FMT, jac->n));
234: PetscCall(PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &sviewer));
235: if (jac->ksp) PetscCall(KSPView(jac->ksp[0], sviewer));
236: PetscCall(PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &sviewer));
237: } else if (isdraw) {
238: PetscDraw draw;
239: char str[25];
240: PetscReal x, y, bottom, h;
242: PetscCall(PetscViewerDrawGetDraw(viewer, 0, &draw));
243: PetscCall(PetscDrawGetCurrentPoint(draw, &x, &y));
244: PetscCall(PetscSNPrintf(str, 25, "Number blocks %" PetscInt_FMT, jac->n));
245: PetscCall(PetscDrawStringBoxed(draw, x, y, PETSC_DRAW_RED, PETSC_DRAW_BLACK, str, NULL, &h));
246: bottom = y - h;
247: PetscCall(PetscDrawPushCurrentPoint(draw, x, bottom));
248: /* warning the communicator on viewer is different then on ksp in parallel */
249: if (jac->ksp) PetscCall(KSPView(jac->ksp[0], viewer));
250: PetscCall(PetscDrawPopCurrentPoint(draw));
251: }
252: PetscFunctionReturn(PETSC_SUCCESS);
253: }
255: static PetscErrorCode PCBJacobiGetSubKSP_BJacobi(PC pc, PetscInt *n_local, PetscInt *first_local, KSP **ksp)
256: {
257: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
259: PetscFunctionBegin;
260: PetscCheck(pc->setupcalled, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Must call KSPSetUp() or PCSetUp() first");
262: if (n_local) *n_local = jac->n_local;
263: if (first_local) *first_local = jac->first_local;
264: if (ksp) *ksp = jac->ksp;
265: PetscFunctionReturn(PETSC_SUCCESS);
266: }
268: static PetscErrorCode PCBJacobiSetTotalBlocks_BJacobi(PC pc, PetscInt blocks, const PetscInt *lens)
269: {
270: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
272: PetscFunctionBegin;
273: PetscCheck(!pc->setupcalled || jac->n == blocks, PetscObjectComm((PetscObject)pc), PETSC_ERR_ORDER, "Cannot alter number of blocks after PCSetUp()/KSPSetUp() has been called");
274: jac->n = blocks;
275: if (!lens) jac->g_lens = NULL;
276: else {
277: PetscCall(PetscMalloc1(blocks, &jac->g_lens));
278: PetscCall(PetscArraycpy(jac->g_lens, lens, blocks));
279: }
280: PetscFunctionReturn(PETSC_SUCCESS);
281: }
283: static PetscErrorCode PCBJacobiGetTotalBlocks_BJacobi(PC pc, PetscInt *blocks, const PetscInt *lens[])
284: {
285: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
287: PetscFunctionBegin;
288: *blocks = jac->n;
289: if (lens) *lens = jac->g_lens;
290: PetscFunctionReturn(PETSC_SUCCESS);
291: }
293: static PetscErrorCode PCBJacobiSetLocalBlocks_BJacobi(PC pc, PetscInt blocks, const PetscInt lens[])
294: {
295: PC_BJacobi *jac;
297: PetscFunctionBegin;
298: jac = (PC_BJacobi *)pc->data;
300: jac->n_local = blocks;
301: if (!lens) jac->l_lens = NULL;
302: else {
303: PetscCall(PetscMalloc1(blocks, &jac->l_lens));
304: PetscCall(PetscArraycpy(jac->l_lens, lens, blocks));
305: }
306: PetscFunctionReturn(PETSC_SUCCESS);
307: }
309: static PetscErrorCode PCBJacobiGetLocalBlocks_BJacobi(PC pc, PetscInt *blocks, const PetscInt *lens[])
310: {
311: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
313: PetscFunctionBegin;
314: *blocks = jac->n_local;
315: if (lens) *lens = jac->l_lens;
316: PetscFunctionReturn(PETSC_SUCCESS);
317: }
319: /*@C
320: PCBJacobiGetSubKSP - Gets the local `KSP` contexts for all blocks on
321: this processor.
323: Not Collective
325: Input Parameter:
326: . pc - the preconditioner context
328: Output Parameters:
329: + n_local - the number of blocks on this processor, or NULL
330: . first_local - the global number of the first block on this processor, or NULL
331: - ksp - the array of KSP contexts
333: Level: advanced
335: Notes:
336: After `PCBJacobiGetSubKSP()` the array of `KSP` contexts is not to be freed.
338: Currently for some matrix implementations only 1 block per processor
339: is supported.
341: You must call `KSPSetUp()` or `PCSetUp()` before calling `PCBJacobiGetSubKSP()`.
343: Fortran Note:
344: Call `PCBJacobiRestoreSubKSP()` when you no longer need access to the array of `KSP`
346: .seealso: [](ch_ksp), `PCBJACOBI`, `PCASM`, `PCASMGetSubKSP()`
347: @*/
348: PetscErrorCode PCBJacobiGetSubKSP(PC pc, PetscInt *n_local, PetscInt *first_local, KSP *ksp[])
349: {
350: PetscFunctionBegin;
352: PetscUseMethod(pc, "PCBJacobiGetSubKSP_C", (PC, PetscInt *, PetscInt *, KSP **), (pc, n_local, first_local, ksp));
353: PetscFunctionReturn(PETSC_SUCCESS);
354: }
356: /*@
357: PCBJacobiSetTotalBlocks - Sets the global number of blocks for the block
358: Jacobi preconditioner.
360: Collective
362: Input Parameters:
363: + pc - the preconditioner context
364: . blocks - the number of blocks
365: - lens - [optional] integer array containing the size of each block
367: Options Database Key:
368: . -pc_bjacobi_blocks blocks - Sets the number of global blocks
370: Level: intermediate
372: Note:
373: Currently only a limited number of blocking configurations are supported.
374: All processors sharing the `PC` must call this routine with the same data.
376: .seealso: [](ch_ksp), `PCBJACOBI`, `PCSetUseAmat()`, `PCBJacobiSetLocalBlocks()`
377: @*/
378: PetscErrorCode PCBJacobiSetTotalBlocks(PC pc, PetscInt blocks, const PetscInt lens[])
379: {
380: PetscFunctionBegin;
382: PetscCheck(blocks > 0, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Must have positive blocks");
383: PetscTryMethod(pc, "PCBJacobiSetTotalBlocks_C", (PC, PetscInt, const PetscInt[]), (pc, blocks, lens));
384: PetscFunctionReturn(PETSC_SUCCESS);
385: }
387: /*@C
388: PCBJacobiGetTotalBlocks - Gets the global number of blocks for the block
389: Jacobi, `PCBJACOBI`, preconditioner.
391: Not Collective
393: Input Parameter:
394: . pc - the preconditioner context
396: Output Parameters:
397: + blocks - the number of blocks
398: - lens - integer array containing the size of each block
400: Level: intermediate
402: .seealso: [](ch_ksp), `PCBJACOBI`, `PCSetUseAmat()`, `PCBJacobiGetLocalBlocks()`
403: @*/
404: PetscErrorCode PCBJacobiGetTotalBlocks(PC pc, PetscInt *blocks, const PetscInt *lens[])
405: {
406: PetscFunctionBegin;
408: PetscAssertPointer(blocks, 2);
409: PetscUseMethod(pc, "PCBJacobiGetTotalBlocks_C", (PC, PetscInt *, const PetscInt *[]), (pc, blocks, lens));
410: PetscFunctionReturn(PETSC_SUCCESS);
411: }
413: /*@
414: PCBJacobiSetLocalBlocks - Sets the local number of blocks for the block
415: Jacobi, `PCBJACOBI`, preconditioner.
417: Not Collective
419: Input Parameters:
420: + pc - the preconditioner context
421: . blocks - the number of blocks
422: - lens - [optional] integer array containing size of each block
424: Options Database Key:
425: . -pc_bjacobi_local_blocks blocks - Sets the number of local blocks
427: Level: intermediate
429: Note:
430: Currently only a limited number of blocking configurations are supported.
432: .seealso: [](ch_ksp), `PCBJACOBI`, `PCSetUseAmat()`, `PCBJacobiSetTotalBlocks()`
433: @*/
434: PetscErrorCode PCBJacobiSetLocalBlocks(PC pc, PetscInt blocks, const PetscInt lens[])
435: {
436: PetscFunctionBegin;
438: PetscCheck(blocks >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must have nonegative blocks");
439: PetscTryMethod(pc, "PCBJacobiSetLocalBlocks_C", (PC, PetscInt, const PetscInt[]), (pc, blocks, lens));
440: PetscFunctionReturn(PETSC_SUCCESS);
441: }
443: /*@C
444: PCBJacobiGetLocalBlocks - Gets the local number of blocks for the block
445: Jacobi, `PCBJACOBI`, preconditioner.
447: Not Collective
449: Input Parameters:
450: + pc - the preconditioner context
451: . blocks - the number of blocks
452: - lens - [optional] integer array containing size of each block
454: Level: intermediate
456: Note:
457: Currently only a limited number of blocking configurations are supported.
459: .seealso: [](ch_ksp), `PCBJACOBI`, `PCSetUseAmat()`, `PCBJacobiGetTotalBlocks()`
460: @*/
461: PetscErrorCode PCBJacobiGetLocalBlocks(PC pc, PetscInt *blocks, const PetscInt *lens[])
462: {
463: PetscFunctionBegin;
465: PetscAssertPointer(blocks, 2);
466: PetscUseMethod(pc, "PCBJacobiGetLocalBlocks_C", (PC, PetscInt *, const PetscInt *[]), (pc, blocks, lens));
467: PetscFunctionReturn(PETSC_SUCCESS);
468: }
470: /*MC
471: PCBJACOBI - Use block Jacobi preconditioning, each block is (approximately) solved with its own `KSP` object.
473: Options Database Keys:
474: + -pc_use_amat (true|false) - use `Amat` to apply block of operator in inner Krylov method
475: - -pc_bjacobi_blocks n - use n total blocks
477: Level: beginner
479: Notes:
480: See `PCJACOBI` for diagonal Jacobi, `PCVPBJACOBI` for variable point block, and `PCPBJACOBI` for fixed size point block
482: Each processor can have one or more blocks, or a single block can be shared by several processes. Defaults to one block per processor.
484: To set options on the solvers for each block append -sub_ to all the `KSP` and `PC`
485: options database keys. For example, -sub_pc_type ilu -sub_pc_factor_levels 1 -sub_ksp_type preonly
487: To set the options on the solvers separate for each block call `PCBJacobiGetSubKSP()`
488: and set the options directly on the resulting `KSP` object (you can access its `PC`
489: `KSPGetPC()`)
491: For GPU-based vectors (`VECCUDA`, `VECViennaCL`) it is recommended to use exactly one block per MPI process for best
492: performance. Different block partitioning may lead to additional data transfers
493: between host and GPU that lead to degraded performance.
495: When multiple processes share a single block, each block encompasses exactly all the unknowns owned its set of processes.
497: .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`,
498: `PCASM`, `PCSetUseAmat()`, `PCGetUseAmat()`, `PCBJacobiGetSubKSP()`, `PCBJacobiSetTotalBlocks()`,
499: `PCBJacobiSetLocalBlocks()`, `PCSetModifySubMatrices()`, `PCJACOBI`, `PCVPBJACOBI`, `PCPBJACOBI`
500: M*/
502: PETSC_EXTERN PetscErrorCode PCCreate_BJacobi(PC pc)
503: {
504: PetscMPIInt rank;
505: PC_BJacobi *jac;
507: PetscFunctionBegin;
508: PetscCall(PetscNew(&jac));
509: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)pc), &rank));
511: pc->ops->apply = NULL;
512: pc->ops->matapply = NULL;
513: pc->ops->applytranspose = NULL;
514: pc->ops->matapplytranspose = NULL;
515: pc->ops->setup = PCSetUp_BJacobi;
516: pc->ops->destroy = PCDestroy_BJacobi;
517: pc->ops->setfromoptions = PCSetFromOptions_BJacobi;
518: pc->ops->view = PCView_BJacobi;
519: pc->ops->applyrichardson = NULL;
521: pc->data = (void *)jac;
522: jac->n = -1;
523: jac->n_local = -1;
524: jac->first_local = rank;
525: jac->ksp = NULL;
526: jac->g_lens = NULL;
527: jac->l_lens = NULL;
528: jac->psubcomm = NULL;
530: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBJacobiGetSubKSP_C", PCBJacobiGetSubKSP_BJacobi));
531: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBJacobiSetTotalBlocks_C", PCBJacobiSetTotalBlocks_BJacobi));
532: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBJacobiGetTotalBlocks_C", PCBJacobiGetTotalBlocks_BJacobi));
533: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBJacobiSetLocalBlocks_C", PCBJacobiSetLocalBlocks_BJacobi));
534: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBJacobiGetLocalBlocks_C", PCBJacobiGetLocalBlocks_BJacobi));
535: PetscFunctionReturn(PETSC_SUCCESS);
536: }
538: /*
539: These are for a single block per processor; works for AIJ, BAIJ; Seq and MPI
540: */
541: static PetscErrorCode PCReset_BJacobi_Singleblock(PC pc)
542: {
543: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
544: PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock *)jac->data;
546: PetscFunctionBegin;
547: PetscCall(KSPReset(jac->ksp[0]));
548: PetscCall(VecDestroy(&bjac->x));
549: PetscCall(VecDestroy(&bjac->y));
550: PetscFunctionReturn(PETSC_SUCCESS);
551: }
553: static PetscErrorCode PCDestroy_BJacobi_Singleblock(PC pc)
554: {
555: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
556: PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock *)jac->data;
558: PetscFunctionBegin;
559: PetscCall(PCReset_BJacobi_Singleblock(pc));
560: PetscCall(KSPDestroy(&jac->ksp[0]));
561: PetscCall(PetscFree(jac->ksp));
562: PetscCall(PetscFree(bjac));
563: PetscCall(PCDestroy_BJacobi(pc));
564: PetscFunctionReturn(PETSC_SUCCESS);
565: }
567: static PetscErrorCode PCSetUpOnBlocks_BJacobi_Singleblock(PC pc)
568: {
569: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
570: KSP subksp = jac->ksp[0];
571: KSPConvergedReason reason;
573: PetscFunctionBegin;
574: PetscCall(KSPSetUp(subksp));
575: PetscCall(KSPGetConvergedReason(subksp, &reason));
576: if (reason == KSP_DIVERGED_PC_FAILED) pc->failedreason = PC_SUBPC_ERROR;
577: PetscFunctionReturn(PETSC_SUCCESS);
578: }
580: static PetscErrorCode PCApply_BJacobi_Singleblock(PC pc, Vec x, Vec y)
581: {
582: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
583: PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock *)jac->data;
585: PetscFunctionBegin;
586: PetscCall(VecGetLocalVectorRead(x, bjac->x));
587: PetscCall(VecGetLocalVector(y, bjac->y));
588: /* Since the inner KSP matrix may point directly to the diagonal block of an MPI matrix the inner
589: matrix may change even if the outer KSP/PC has not updated the preconditioner, this will trigger a rebuild
590: of the inner preconditioner automatically unless we pass down the outer preconditioners reuse flag.*/
591: PetscCall(KSPSetReusePreconditioner(jac->ksp[0], pc->reusepreconditioner));
592: PetscCall(PetscLogEventBegin(PC_ApplyOnBlocks, jac->ksp[0], bjac->x, bjac->y, 0));
593: PetscCall(KSPSolve(jac->ksp[0], bjac->x, bjac->y));
594: PetscCall(KSPCheckSolve(jac->ksp[0], pc, bjac->y));
595: PetscCall(PetscLogEventEnd(PC_ApplyOnBlocks, jac->ksp[0], bjac->x, bjac->y, 0));
596: PetscCall(VecRestoreLocalVectorRead(x, bjac->x));
597: PetscCall(VecRestoreLocalVector(y, bjac->y));
598: PetscFunctionReturn(PETSC_SUCCESS);
599: }
601: static PetscErrorCode PCMatApply_BJacobi_Singleblock_Private(PC pc, Mat X, Mat Y, PetscBool transpose)
602: {
603: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
604: Mat sX, sY;
606: PetscFunctionBegin;
607: /* Since the inner KSP matrix may point directly to the diagonal block of an MPI matrix the inner
608: matrix may change even if the outer KSP/PC has not updated the preconditioner, this will trigger a rebuild
609: of the inner preconditioner automatically unless we pass down the outer preconditioners reuse flag.*/
610: PetscCall(KSPSetReusePreconditioner(jac->ksp[0], pc->reusepreconditioner));
611: PetscCall(MatDenseGetLocalMatrix(X, &sX));
612: PetscCall(MatDenseGetLocalMatrix(Y, &sY));
613: if (!transpose) {
614: PetscCall(PetscLogEventBegin(PC_ApplyOnBlocks, jac->ksp[0], sX, sY, 0));
615: PetscCall(KSPMatSolve(jac->ksp[0], sX, sY));
616: PetscCall(PetscLogEventEnd(PC_ApplyOnBlocks, jac->ksp[0], sX, sY, 0));
617: } else {
618: PetscCall(PetscLogEventBegin(PC_ApplyTransposeOnBlocks, jac->ksp[0], sX, sY, 0));
619: PetscCall(KSPMatSolveTranspose(jac->ksp[0], sX, sY));
620: PetscCall(PetscLogEventEnd(PC_ApplyTransposeOnBlocks, jac->ksp[0], sX, sY, 0));
621: }
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;
1192: KSPConvergedReason reason;
1194: PetscFunctionBegin;
1195: /* place x's and y's local arrays into xsub and ysub */
1196: PetscCall(VecGetArrayRead(x, &xarray));
1197: PetscCall(VecGetArray(y, &yarray));
1198: PetscCall(VecPlaceArray(mpjac->xsub, xarray));
1199: PetscCall(VecPlaceArray(mpjac->ysub, yarray));
1201: /* apply preconditioner on each matrix block */
1202: PetscCall(PetscLogEventBegin(PC_ApplyOnBlocks, jac->ksp[0], mpjac->xsub, mpjac->ysub, 0));
1203: PetscCall(KSPSolve(jac->ksp[0], mpjac->xsub, mpjac->ysub));
1204: PetscCall(KSPCheckSolve(jac->ksp[0], pc, mpjac->ysub));
1205: PetscCall(PetscLogEventEnd(PC_ApplyOnBlocks, jac->ksp[0], mpjac->xsub, mpjac->ysub, 0));
1206: PetscCall(KSPGetConvergedReason(jac->ksp[0], &reason));
1207: if (reason == KSP_DIVERGED_PC_FAILED) pc->failedreason = PC_SUBPC_ERROR;
1209: PetscCall(VecResetArray(mpjac->xsub));
1210: PetscCall(VecResetArray(mpjac->ysub));
1211: PetscCall(VecRestoreArrayRead(x, &xarray));
1212: PetscCall(VecRestoreArray(y, &yarray));
1213: PetscFunctionReturn(PETSC_SUCCESS);
1214: }
1216: static PetscErrorCode PCMatApply_BJacobi_Multiproc(PC pc, Mat X, Mat Y)
1217: {
1218: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
1219: KSPConvergedReason reason;
1220: Mat sX, sY;
1221: const PetscScalar *x;
1222: PetscScalar *y;
1223: PetscInt m, N, lda, ldb;
1225: PetscFunctionBegin;
1226: /* apply preconditioner on each matrix block */
1227: PetscCall(MatGetLocalSize(X, &m, NULL));
1228: PetscCall(MatGetSize(X, NULL, &N));
1229: PetscCall(MatDenseGetLDA(X, &lda));
1230: PetscCall(MatDenseGetLDA(Y, &ldb));
1231: PetscCall(MatDenseGetArrayRead(X, &x));
1232: PetscCall(MatDenseGetArrayWrite(Y, &y));
1233: PetscCall(MatCreateDense(PetscObjectComm((PetscObject)jac->ksp[0]), m, PETSC_DECIDE, PETSC_DECIDE, N, (PetscScalar *)x, &sX));
1234: PetscCall(MatCreateDense(PetscObjectComm((PetscObject)jac->ksp[0]), m, PETSC_DECIDE, PETSC_DECIDE, N, y, &sY));
1235: PetscCall(MatDenseSetLDA(sX, lda));
1236: PetscCall(MatDenseSetLDA(sY, ldb));
1237: PetscCall(PetscLogEventBegin(PC_ApplyOnBlocks, jac->ksp[0], X, Y, 0));
1238: PetscCall(KSPMatSolve(jac->ksp[0], sX, sY));
1239: PetscCall(KSPCheckSolve(jac->ksp[0], pc, NULL));
1240: PetscCall(PetscLogEventEnd(PC_ApplyOnBlocks, jac->ksp[0], X, Y, 0));
1241: PetscCall(MatDestroy(&sY));
1242: PetscCall(MatDestroy(&sX));
1243: PetscCall(MatDenseRestoreArrayWrite(Y, &y));
1244: PetscCall(MatDenseRestoreArrayRead(X, &x));
1245: PetscCall(KSPGetConvergedReason(jac->ksp[0], &reason));
1246: if (reason == KSP_DIVERGED_PC_FAILED) pc->failedreason = PC_SUBPC_ERROR;
1247: PetscFunctionReturn(PETSC_SUCCESS);
1248: }
1250: static PetscErrorCode PCSetUp_BJacobi_Multiproc(PC pc)
1251: {
1252: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
1253: PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc *)jac->data;
1254: PetscInt m, n;
1255: MPI_Comm comm, subcomm = 0;
1256: const char *prefix;
1257: PetscBool wasSetup = PETSC_TRUE;
1258: VecType vectype;
1260: PetscFunctionBegin;
1261: PetscCall(PetscObjectGetComm((PetscObject)pc, &comm));
1262: PetscCheck(jac->n_local <= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Only a single block in a subcommunicator is supported");
1263: jac->n_local = 1; /* currently only a single block is supported for a subcommunicator */
1264: if (!pc->setupcalled) {
1265: PetscInt nestlevel;
1267: wasSetup = PETSC_FALSE;
1268: PetscCall(PetscNew(&mpjac));
1269: jac->data = (void *)mpjac;
1271: /* initialize datastructure mpjac */
1272: if (!jac->psubcomm) {
1273: /* Create default contiguous subcommunicatiors if user does not provide them */
1274: PetscCall(PetscSubcommCreate(comm, &jac->psubcomm));
1275: PetscCall(PetscSubcommSetNumber(jac->psubcomm, jac->n));
1276: PetscCall(PetscSubcommSetType(jac->psubcomm, PETSC_SUBCOMM_CONTIGUOUS));
1277: }
1278: mpjac->psubcomm = jac->psubcomm;
1279: subcomm = PetscSubcommChild(mpjac->psubcomm);
1281: /* Get matrix blocks of pmat */
1282: PetscCall(MatGetMultiProcBlock(pc->pmat, subcomm, MAT_INITIAL_MATRIX, &mpjac->submats));
1284: /* create a new PC that processors in each subcomm have copy of */
1285: PetscCall(PetscMalloc1(1, &jac->ksp));
1286: PetscCall(KSPCreate(subcomm, &jac->ksp[0]));
1287: PetscCall(PCGetKSPNestLevel(pc, &nestlevel));
1288: PetscCall(KSPSetNestLevel(jac->ksp[0], nestlevel + 1));
1289: PetscCall(KSPSetErrorIfNotConverged(jac->ksp[0], pc->erroriffailure));
1290: PetscCall(PetscObjectIncrementTabLevel((PetscObject)jac->ksp[0], (PetscObject)pc, 1));
1291: PetscCall(KSPSetOperators(jac->ksp[0], mpjac->submats, mpjac->submats));
1292: PetscCall(KSPGetPC(jac->ksp[0], &mpjac->pc));
1294: PetscCall(PCGetOptionsPrefix(pc, &prefix));
1295: PetscCall(KSPSetOptionsPrefix(jac->ksp[0], prefix));
1296: PetscCall(KSPAppendOptionsPrefix(jac->ksp[0], "sub_"));
1297: PetscCall(KSPGetOptionsPrefix(jac->ksp[0], &prefix));
1298: PetscCall(MatSetOptionsPrefix(mpjac->submats, prefix));
1300: PetscCall(MatGetLocalSize(mpjac->submats, &m, &n));
1301: PetscCall(VecCreateMPIWithArray(subcomm, 1, n, PETSC_DECIDE, NULL, &mpjac->xsub));
1302: PetscCall(VecCreateMPIWithArray(subcomm, 1, m, PETSC_DECIDE, NULL, &mpjac->ysub));
1303: PetscCall(MatGetVecType(mpjac->submats, &vectype));
1304: PetscCall(VecSetType(mpjac->xsub, vectype));
1305: PetscCall(VecSetType(mpjac->ysub, vectype));
1307: pc->ops->setuponblocks = PCSetUpOnBlocks_BJacobi_Multiproc;
1308: pc->ops->reset = PCReset_BJacobi_Multiproc;
1309: pc->ops->destroy = PCDestroy_BJacobi_Multiproc;
1310: pc->ops->apply = PCApply_BJacobi_Multiproc;
1311: pc->ops->matapply = PCMatApply_BJacobi_Multiproc;
1312: } else { /* pc->setupcalled */
1313: subcomm = PetscSubcommChild(mpjac->psubcomm);
1314: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
1315: /* destroy old matrix blocks, then get new matrix blocks */
1316: if (mpjac->submats) PetscCall(MatDestroy(&mpjac->submats));
1317: PetscCall(MatGetMultiProcBlock(pc->pmat, subcomm, MAT_INITIAL_MATRIX, &mpjac->submats));
1318: } else {
1319: PetscCall(MatGetMultiProcBlock(pc->pmat, subcomm, MAT_REUSE_MATRIX, &mpjac->submats));
1320: }
1321: PetscCall(KSPSetOperators(jac->ksp[0], mpjac->submats, mpjac->submats));
1322: }
1324: if (!wasSetup && pc->setfromoptionscalled) PetscCall(KSPSetFromOptions(jac->ksp[0]));
1325: PetscFunctionReturn(PETSC_SUCCESS);
1326: }