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