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
475: its own `KSP` object.
477: Options Database Keys:
478: + -pc_use_amat - use Amat to apply block of operator in inner Krylov method
479: - -pc_bjacobi_blocks <n> - use n total blocks
481: Level: beginner
483: Notes:
484: See `PCJACOBI` for diagonal Jacobi, `PCVPBJACOBI` for variable point block, and `PCPBJACOBI` for fixed size point block
486: Each processor can have one or more blocks, or a single block can be shared by several processes. Defaults to one block per processor.
488: To set options on the solvers for each block append -sub_ to all the `KSP` and `PC`
489: options database keys. For example, -sub_pc_type ilu -sub_pc_factor_levels 1 -sub_ksp_type preonly
491: To set the options on the solvers separate for each block call `PCBJacobiGetSubKSP()`
492: and set the options directly on the resulting `KSP` object (you can access its `PC`
493: `KSPGetPC()`)
495: For GPU-based vectors (`VECCUDA`, `VECViennaCL`) it is recommended to use exactly one block per MPI process for best
496: performance. Different block partitioning may lead to additional data transfers
497: between host and GPU that lead to degraded performance.
499: When multiple processes share a single block, each block encompasses exactly all the unknowns owned its set of processes.
501: .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCType`,
502: `PCASM`, `PCSetUseAmat()`, `PCGetUseAmat()`, `PCBJacobiGetSubKSP()`, `PCBJacobiSetTotalBlocks()`,
503: `PCBJacobiSetLocalBlocks()`, `PCSetModifySubMatrices()`, `PCJACOBI`, `PCVPBJACOBI`, `PCPBJACOBI`
504: M*/
506: PETSC_EXTERN PetscErrorCode PCCreate_BJacobi(PC pc)
507: {
508: PetscMPIInt rank;
509: PC_BJacobi *jac;
511: PetscFunctionBegin;
512: PetscCall(PetscNew(&jac));
513: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)pc), &rank));
515: pc->ops->apply = NULL;
516: pc->ops->matapply = NULL;
517: pc->ops->applytranspose = NULL;
518: pc->ops->matapplytranspose = NULL;
519: pc->ops->setup = PCSetUp_BJacobi;
520: pc->ops->destroy = PCDestroy_BJacobi;
521: pc->ops->setfromoptions = PCSetFromOptions_BJacobi;
522: pc->ops->view = PCView_BJacobi;
523: pc->ops->applyrichardson = NULL;
525: pc->data = (void *)jac;
526: jac->n = -1;
527: jac->n_local = -1;
528: jac->first_local = rank;
529: jac->ksp = NULL;
530: jac->g_lens = NULL;
531: jac->l_lens = NULL;
532: jac->psubcomm = NULL;
534: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBJacobiGetSubKSP_C", PCBJacobiGetSubKSP_BJacobi));
535: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBJacobiSetTotalBlocks_C", PCBJacobiSetTotalBlocks_BJacobi));
536: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBJacobiGetTotalBlocks_C", PCBJacobiGetTotalBlocks_BJacobi));
537: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBJacobiSetLocalBlocks_C", PCBJacobiSetLocalBlocks_BJacobi));
538: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCBJacobiGetLocalBlocks_C", PCBJacobiGetLocalBlocks_BJacobi));
539: PetscFunctionReturn(PETSC_SUCCESS);
540: }
542: /*
543: These are for a single block per processor; works for AIJ, BAIJ; Seq and MPI
544: */
545: static PetscErrorCode PCReset_BJacobi_Singleblock(PC pc)
546: {
547: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
548: PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock *)jac->data;
550: PetscFunctionBegin;
551: PetscCall(KSPReset(jac->ksp[0]));
552: PetscCall(VecDestroy(&bjac->x));
553: PetscCall(VecDestroy(&bjac->y));
554: PetscFunctionReturn(PETSC_SUCCESS);
555: }
557: static PetscErrorCode PCDestroy_BJacobi_Singleblock(PC pc)
558: {
559: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
560: PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock *)jac->data;
562: PetscFunctionBegin;
563: PetscCall(PCReset_BJacobi_Singleblock(pc));
564: PetscCall(KSPDestroy(&jac->ksp[0]));
565: PetscCall(PetscFree(jac->ksp));
566: PetscCall(PetscFree(bjac));
567: PetscCall(PCDestroy_BJacobi(pc));
568: PetscFunctionReturn(PETSC_SUCCESS);
569: }
571: static PetscErrorCode PCSetUpOnBlocks_BJacobi_Singleblock(PC pc)
572: {
573: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
574: KSP subksp = jac->ksp[0];
575: KSPConvergedReason reason;
577: PetscFunctionBegin;
578: PetscCall(KSPSetUp(subksp));
579: PetscCall(KSPGetConvergedReason(subksp, &reason));
580: if (reason == KSP_DIVERGED_PC_FAILED) pc->failedreason = PC_SUBPC_ERROR;
581: PetscFunctionReturn(PETSC_SUCCESS);
582: }
584: static PetscErrorCode PCApply_BJacobi_Singleblock(PC pc, Vec x, Vec y)
585: {
586: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
587: PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock *)jac->data;
589: PetscFunctionBegin;
590: PetscCall(VecGetLocalVectorRead(x, bjac->x));
591: PetscCall(VecGetLocalVector(y, bjac->y));
592: /* Since the inner KSP matrix may point directly to the diagonal block of an MPI matrix the inner
593: matrix may change even if the outer KSP/PC has not updated the preconditioner, this will trigger a rebuild
594: of the inner preconditioner automatically unless we pass down the outer preconditioners reuse flag.*/
595: PetscCall(KSPSetReusePreconditioner(jac->ksp[0], pc->reusepreconditioner));
596: PetscCall(PetscLogEventBegin(PC_ApplyOnBlocks, jac->ksp[0], bjac->x, bjac->y, 0));
597: PetscCall(KSPSolve(jac->ksp[0], bjac->x, bjac->y));
598: PetscCall(KSPCheckSolve(jac->ksp[0], pc, bjac->y));
599: PetscCall(PetscLogEventEnd(PC_ApplyOnBlocks, jac->ksp[0], bjac->x, bjac->y, 0));
600: PetscCall(VecRestoreLocalVectorRead(x, bjac->x));
601: PetscCall(VecRestoreLocalVector(y, bjac->y));
602: PetscFunctionReturn(PETSC_SUCCESS);
603: }
605: static PetscErrorCode PCMatApply_BJacobi_Singleblock_Private(PC pc, Mat X, Mat Y, PetscBool transpose)
606: {
607: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
608: Mat sX, sY;
610: PetscFunctionBegin;
611: /* Since the inner KSP matrix may point directly to the diagonal block of an MPI matrix the inner
612: matrix may change even if the outer KSP/PC has not updated the preconditioner, this will trigger a rebuild
613: of the inner preconditioner automatically unless we pass down the outer preconditioners reuse flag.*/
614: PetscCall(KSPSetReusePreconditioner(jac->ksp[0], pc->reusepreconditioner));
615: PetscCall(MatDenseGetLocalMatrix(X, &sX));
616: PetscCall(MatDenseGetLocalMatrix(Y, &sY));
617: if (!transpose) {
618: PetscCall(PetscLogEventBegin(PC_ApplyOnBlocks, jac->ksp[0], sX, sY, 0));
619: PetscCall(KSPMatSolve(jac->ksp[0], sX, sY));
620: PetscCall(PetscLogEventEnd(PC_ApplyOnBlocks, jac->ksp[0], sX, sY, 0));
621: } else {
622: PetscCall(PetscLogEventBegin(PC_ApplyTransposeOnBlocks, jac->ksp[0], sX, sY, 0));
623: PetscCall(KSPMatSolveTranspose(jac->ksp[0], sX, sY));
624: PetscCall(PetscLogEventEnd(PC_ApplyTransposeOnBlocks, jac->ksp[0], sX, sY, 0));
625: }
626: PetscFunctionReturn(PETSC_SUCCESS);
627: }
629: static PetscErrorCode PCMatApply_BJacobi_Singleblock(PC pc, Mat X, Mat Y)
630: {
631: PetscFunctionBegin;
632: PetscCall(PCMatApply_BJacobi_Singleblock_Private(pc, X, Y, PETSC_FALSE));
633: PetscFunctionReturn(PETSC_SUCCESS);
634: }
636: static PetscErrorCode PCMatApplyTranspose_BJacobi_Singleblock(PC pc, Mat X, Mat Y)
637: {
638: PetscFunctionBegin;
639: PetscCall(PCMatApply_BJacobi_Singleblock_Private(pc, X, Y, PETSC_TRUE));
640: PetscFunctionReturn(PETSC_SUCCESS);
641: }
643: static PetscErrorCode PCApplySymmetricLeft_BJacobi_Singleblock(PC pc, Vec x, Vec y)
644: {
645: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
646: PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock *)jac->data;
647: PetscScalar *y_array;
648: const PetscScalar *x_array;
649: PC subpc;
651: PetscFunctionBegin;
652: /*
653: The VecPlaceArray() is to avoid having to copy the
654: y vector into the bjac->x vector. The reason for
655: the bjac->x vector is that we need a sequential vector
656: for the sequential solve.
657: */
658: PetscCall(VecGetArrayRead(x, &x_array));
659: PetscCall(VecGetArray(y, &y_array));
660: PetscCall(VecPlaceArray(bjac->x, x_array));
661: PetscCall(VecPlaceArray(bjac->y, y_array));
662: /* apply the symmetric left portion of the inner PC operator */
663: /* note this bypasses the inner KSP and its options completely */
664: PetscCall(KSPGetPC(jac->ksp[0], &subpc));
665: PetscCall(PCApplySymmetricLeft(subpc, bjac->x, bjac->y));
666: PetscCall(VecResetArray(bjac->x));
667: PetscCall(VecResetArray(bjac->y));
668: PetscCall(VecRestoreArrayRead(x, &x_array));
669: PetscCall(VecRestoreArray(y, &y_array));
670: PetscFunctionReturn(PETSC_SUCCESS);
671: }
673: static PetscErrorCode PCApplySymmetricRight_BJacobi_Singleblock(PC pc, Vec x, Vec y)
674: {
675: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
676: PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock *)jac->data;
677: PetscScalar *y_array;
678: const PetscScalar *x_array;
679: PC subpc;
681: PetscFunctionBegin;
682: /*
683: The VecPlaceArray() is to avoid having to copy the
684: y vector into the bjac->x vector. The reason for
685: the bjac->x vector is that we need a sequential vector
686: for the sequential solve.
687: */
688: PetscCall(VecGetArrayRead(x, &x_array));
689: PetscCall(VecGetArray(y, &y_array));
690: PetscCall(VecPlaceArray(bjac->x, x_array));
691: PetscCall(VecPlaceArray(bjac->y, y_array));
693: /* apply the symmetric right portion of the inner PC operator */
694: /* note this bypasses the inner KSP and its options completely */
696: PetscCall(KSPGetPC(jac->ksp[0], &subpc));
697: PetscCall(PCApplySymmetricRight(subpc, bjac->x, bjac->y));
699: PetscCall(VecResetArray(bjac->x));
700: PetscCall(VecResetArray(bjac->y));
701: PetscCall(VecRestoreArrayRead(x, &x_array));
702: PetscCall(VecRestoreArray(y, &y_array));
703: PetscFunctionReturn(PETSC_SUCCESS);
704: }
706: static PetscErrorCode PCApplyTranspose_BJacobi_Singleblock(PC pc, Vec x, Vec y)
707: {
708: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
709: PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock *)jac->data;
710: PetscScalar *y_array;
711: const PetscScalar *x_array;
713: PetscFunctionBegin;
714: /*
715: The VecPlaceArray() is to avoid having to copy the
716: y vector into the bjac->x vector. The reason for
717: the bjac->x vector is that we need a sequential vector
718: for the sequential solve.
719: */
720: PetscCall(VecGetArrayRead(x, &x_array));
721: PetscCall(VecGetArray(y, &y_array));
722: PetscCall(VecPlaceArray(bjac->x, x_array));
723: PetscCall(VecPlaceArray(bjac->y, y_array));
724: PetscCall(PetscLogEventBegin(PC_ApplyTransposeOnBlocks, jac->ksp[0], bjac->x, bjac->y, 0));
725: PetscCall(KSPSolveTranspose(jac->ksp[0], bjac->x, bjac->y));
726: PetscCall(KSPCheckSolve(jac->ksp[0], pc, bjac->y));
727: PetscCall(PetscLogEventEnd(PC_ApplyTransposeOnBlocks, jac->ksp[0], bjac->x, bjac->y, 0));
728: PetscCall(VecResetArray(bjac->x));
729: PetscCall(VecResetArray(bjac->y));
730: PetscCall(VecRestoreArrayRead(x, &x_array));
731: PetscCall(VecRestoreArray(y, &y_array));
732: PetscFunctionReturn(PETSC_SUCCESS);
733: }
735: static PetscErrorCode PCSetUp_BJacobi_Singleblock(PC pc, Mat mat, Mat pmat)
736: {
737: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
738: PetscInt m;
739: KSP ksp;
740: PC_BJacobi_Singleblock *bjac;
741: PetscBool wasSetup = PETSC_TRUE;
742: VecType vectype;
743: const char *prefix;
745: PetscFunctionBegin;
746: if (!pc->setupcalled) {
747: if (!jac->ksp) {
748: PetscInt nestlevel;
750: wasSetup = PETSC_FALSE;
752: PetscCall(KSPCreate(PETSC_COMM_SELF, &ksp));
753: PetscCall(PCGetKSPNestLevel(pc, &nestlevel));
754: PetscCall(KSPSetNestLevel(ksp, nestlevel + 1));
755: PetscCall(KSPSetErrorIfNotConverged(ksp, pc->erroriffailure));
756: PetscCall(PetscObjectIncrementTabLevel((PetscObject)ksp, (PetscObject)pc, 1));
757: PetscCall(KSPSetType(ksp, KSPPREONLY));
758: PetscCall(PCGetOptionsPrefix(pc, &prefix));
759: PetscCall(KSPSetOptionsPrefix(ksp, prefix));
760: PetscCall(KSPAppendOptionsPrefix(ksp, "sub_"));
762: pc->ops->reset = PCReset_BJacobi_Singleblock;
763: pc->ops->destroy = PCDestroy_BJacobi_Singleblock;
764: pc->ops->apply = PCApply_BJacobi_Singleblock;
765: pc->ops->matapply = PCMatApply_BJacobi_Singleblock;
766: pc->ops->matapplytranspose = PCMatApplyTranspose_BJacobi_Singleblock;
767: pc->ops->applysymmetricleft = PCApplySymmetricLeft_BJacobi_Singleblock;
768: pc->ops->applysymmetricright = PCApplySymmetricRight_BJacobi_Singleblock;
769: pc->ops->applytranspose = PCApplyTranspose_BJacobi_Singleblock;
770: pc->ops->setuponblocks = PCSetUpOnBlocks_BJacobi_Singleblock;
772: PetscCall(PetscMalloc1(1, &jac->ksp));
773: jac->ksp[0] = ksp;
775: PetscCall(PetscNew(&bjac));
776: jac->data = (void *)bjac;
777: } else {
778: ksp = jac->ksp[0];
779: bjac = (PC_BJacobi_Singleblock *)jac->data;
780: }
782: /*
783: The reason we need to generate these vectors is to serve
784: as the right-hand side and solution vector for the solve on the
785: block. We do not need to allocate space for the vectors since
786: that is provided via VecPlaceArray() just before the call to
787: KSPSolve() on the block.
788: */
789: PetscCall(MatGetSize(pmat, &m, &m));
790: PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, m, NULL, &bjac->x));
791: PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, m, NULL, &bjac->y));
792: PetscCall(MatGetVecType(pmat, &vectype));
793: PetscCall(VecSetType(bjac->x, vectype));
794: PetscCall(VecSetType(bjac->y, vectype));
795: } else {
796: ksp = jac->ksp[0];
797: bjac = (PC_BJacobi_Singleblock *)jac->data;
798: }
799: PetscCall(KSPGetOptionsPrefix(ksp, &prefix));
800: if (pc->useAmat) {
801: PetscCall(KSPSetOperators(ksp, mat, pmat));
802: PetscCall(MatSetOptionsPrefix(mat, prefix));
803: } else {
804: PetscCall(KSPSetOperators(ksp, pmat, pmat));
805: }
806: PetscCall(MatSetOptionsPrefix(pmat, prefix));
807: if (!wasSetup && pc->setfromoptionscalled) {
808: /* If PCSetFromOptions_BJacobi is called later, KSPSetFromOptions will be called at that time. */
809: PetscCall(KSPSetFromOptions(ksp));
810: }
811: PetscFunctionReturn(PETSC_SUCCESS);
812: }
814: static PetscErrorCode PCReset_BJacobi_Multiblock(PC pc)
815: {
816: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
817: PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock *)jac->data;
818: PetscInt i;
820: PetscFunctionBegin;
821: if (bjac && bjac->pmat) {
822: PetscCall(MatDestroyMatrices(jac->n_local, &bjac->pmat));
823: if (pc->useAmat) PetscCall(MatDestroyMatrices(jac->n_local, &bjac->mat));
824: }
826: for (i = 0; i < jac->n_local; i++) {
827: PetscCall(KSPReset(jac->ksp[i]));
828: if (bjac && bjac->x) {
829: PetscCall(VecDestroy(&bjac->x[i]));
830: PetscCall(VecDestroy(&bjac->y[i]));
831: PetscCall(ISDestroy(&bjac->is[i]));
832: }
833: }
834: PetscCall(PetscFree(jac->l_lens));
835: PetscCall(PetscFree(jac->g_lens));
836: PetscFunctionReturn(PETSC_SUCCESS);
837: }
839: static PetscErrorCode PCDestroy_BJacobi_Multiblock(PC pc)
840: {
841: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
842: PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock *)jac->data;
843: PetscInt i;
845: PetscFunctionBegin;
846: PetscCall(PCReset_BJacobi_Multiblock(pc));
847: if (bjac) {
848: PetscCall(PetscFree2(bjac->x, bjac->y));
849: PetscCall(PetscFree(bjac->starts));
850: PetscCall(PetscFree(bjac->is));
851: }
852: PetscCall(PetscFree(jac->data));
853: for (i = 0; i < jac->n_local; i++) PetscCall(KSPDestroy(&jac->ksp[i]));
854: PetscCall(PetscFree(jac->ksp));
855: PetscCall(PCDestroy_BJacobi(pc));
856: PetscFunctionReturn(PETSC_SUCCESS);
857: }
859: static PetscErrorCode PCSetUpOnBlocks_BJacobi_Multiblock(PC pc)
860: {
861: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
862: PetscInt i, n_local = jac->n_local;
863: KSPConvergedReason reason;
865: PetscFunctionBegin;
866: for (i = 0; i < n_local; i++) {
867: PetscCall(KSPSetUp(jac->ksp[i]));
868: PetscCall(KSPGetConvergedReason(jac->ksp[i], &reason));
869: if (reason == KSP_DIVERGED_PC_FAILED) pc->failedreason = PC_SUBPC_ERROR;
870: }
871: PetscFunctionReturn(PETSC_SUCCESS);
872: }
874: static PetscErrorCode PCApply_BJacobi_Multiblock(PC pc, Vec x, Vec y)
875: {
876: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
877: PetscInt i, n_local = jac->n_local;
878: PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock *)jac->data;
879: PetscScalar *yin;
880: const PetscScalar *xin;
882: PetscFunctionBegin;
883: PetscCall(VecGetArrayRead(x, &xin));
884: PetscCall(VecGetArray(y, &yin));
885: for (i = 0; i < n_local; i++) {
886: /*
887: To avoid copying the subvector from x into a workspace we instead
888: make the workspace vector array point to the subpart of the array of
889: the global vector.
890: */
891: PetscCall(VecPlaceArray(bjac->x[i], xin + bjac->starts[i]));
892: PetscCall(VecPlaceArray(bjac->y[i], yin + bjac->starts[i]));
894: PetscCall(PetscLogEventBegin(PC_ApplyOnBlocks, jac->ksp[i], bjac->x[i], bjac->y[i], 0));
895: PetscCall(KSPSolve(jac->ksp[i], bjac->x[i], bjac->y[i]));
896: PetscCall(KSPCheckSolve(jac->ksp[i], pc, bjac->y[i]));
897: PetscCall(PetscLogEventEnd(PC_ApplyOnBlocks, jac->ksp[i], bjac->x[i], bjac->y[i], 0));
899: PetscCall(VecResetArray(bjac->x[i]));
900: PetscCall(VecResetArray(bjac->y[i]));
901: }
902: PetscCall(VecRestoreArrayRead(x, &xin));
903: PetscCall(VecRestoreArray(y, &yin));
904: PetscFunctionReturn(PETSC_SUCCESS);
905: }
907: static PetscErrorCode PCApplySymmetricLeft_BJacobi_Multiblock(PC pc, Vec x, Vec y)
908: {
909: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
910: PetscInt i, n_local = jac->n_local;
911: PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock *)jac->data;
912: PetscScalar *yin;
913: const PetscScalar *xin;
914: PC subpc;
916: PetscFunctionBegin;
917: PetscCall(VecGetArrayRead(x, &xin));
918: PetscCall(VecGetArray(y, &yin));
919: for (i = 0; i < n_local; i++) {
920: /*
921: To avoid copying the subvector from x into a workspace we instead
922: make the workspace vector array point to the subpart of the array of
923: the global vector.
924: */
925: PetscCall(VecPlaceArray(bjac->x[i], xin + bjac->starts[i]));
926: PetscCall(VecPlaceArray(bjac->y[i], yin + bjac->starts[i]));
928: PetscCall(PetscLogEventBegin(PC_ApplyOnBlocks, jac->ksp[i], bjac->x[i], bjac->y[i], 0));
929: /* apply the symmetric left portion of the inner PC operator */
930: /* note this bypasses the inner KSP and its options completely */
931: PetscCall(KSPGetPC(jac->ksp[i], &subpc));
932: PetscCall(PCApplySymmetricLeft(subpc, bjac->x[i], bjac->y[i]));
933: PetscCall(PetscLogEventEnd(PC_ApplyOnBlocks, jac->ksp[i], bjac->x[i], bjac->y[i], 0));
935: PetscCall(VecResetArray(bjac->x[i]));
936: PetscCall(VecResetArray(bjac->y[i]));
937: }
938: PetscCall(VecRestoreArrayRead(x, &xin));
939: PetscCall(VecRestoreArray(y, &yin));
940: PetscFunctionReturn(PETSC_SUCCESS);
941: }
943: static PetscErrorCode PCApplySymmetricRight_BJacobi_Multiblock(PC pc, Vec x, Vec y)
944: {
945: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
946: PetscInt i, n_local = jac->n_local;
947: PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock *)jac->data;
948: PetscScalar *yin;
949: const PetscScalar *xin;
950: PC subpc;
952: PetscFunctionBegin;
953: PetscCall(VecGetArrayRead(x, &xin));
954: PetscCall(VecGetArray(y, &yin));
955: for (i = 0; i < n_local; i++) {
956: /*
957: To avoid copying the subvector from x into a workspace we instead
958: make the workspace vector array point to the subpart of the array of
959: the global vector.
960: */
961: PetscCall(VecPlaceArray(bjac->x[i], xin + bjac->starts[i]));
962: PetscCall(VecPlaceArray(bjac->y[i], yin + bjac->starts[i]));
964: PetscCall(PetscLogEventBegin(PC_ApplyOnBlocks, jac->ksp[i], bjac->x[i], bjac->y[i], 0));
965: /* apply the symmetric left portion of the inner PC operator */
966: /* note this bypasses the inner KSP and its options completely */
967: PetscCall(KSPGetPC(jac->ksp[i], &subpc));
968: PetscCall(PCApplySymmetricRight(subpc, bjac->x[i], bjac->y[i]));
969: PetscCall(PetscLogEventEnd(PC_ApplyOnBlocks, jac->ksp[i], bjac->x[i], bjac->y[i], 0));
971: PetscCall(VecResetArray(bjac->x[i]));
972: PetscCall(VecResetArray(bjac->y[i]));
973: }
974: PetscCall(VecRestoreArrayRead(x, &xin));
975: PetscCall(VecRestoreArray(y, &yin));
976: PetscFunctionReturn(PETSC_SUCCESS);
977: }
979: static PetscErrorCode PCApplyTranspose_BJacobi_Multiblock(PC pc, Vec x, Vec y)
980: {
981: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
982: PetscInt i, n_local = jac->n_local;
983: PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock *)jac->data;
984: PetscScalar *yin;
985: const PetscScalar *xin;
987: PetscFunctionBegin;
988: PetscCall(VecGetArrayRead(x, &xin));
989: PetscCall(VecGetArray(y, &yin));
990: for (i = 0; i < n_local; i++) {
991: /*
992: To avoid copying the subvector from x into a workspace we instead
993: make the workspace vector array point to the subpart of the array of
994: the global vector.
995: */
996: PetscCall(VecPlaceArray(bjac->x[i], xin + bjac->starts[i]));
997: PetscCall(VecPlaceArray(bjac->y[i], yin + bjac->starts[i]));
999: PetscCall(PetscLogEventBegin(PC_ApplyTransposeOnBlocks, jac->ksp[i], bjac->x[i], bjac->y[i], 0));
1000: PetscCall(KSPSolveTranspose(jac->ksp[i], bjac->x[i], bjac->y[i]));
1001: PetscCall(KSPCheckSolve(jac->ksp[i], pc, bjac->y[i]));
1002: PetscCall(PetscLogEventEnd(PC_ApplyTransposeOnBlocks, jac->ksp[i], bjac->x[i], bjac->y[i], 0));
1004: PetscCall(VecResetArray(bjac->x[i]));
1005: PetscCall(VecResetArray(bjac->y[i]));
1006: }
1007: PetscCall(VecRestoreArrayRead(x, &xin));
1008: PetscCall(VecRestoreArray(y, &yin));
1009: PetscFunctionReturn(PETSC_SUCCESS);
1010: }
1012: static PetscErrorCode PCSetUp_BJacobi_Multiblock(PC pc, Mat mat, Mat pmat)
1013: {
1014: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
1015: PetscInt m, n_local, N, M, start, i;
1016: const char *prefix;
1017: KSP ksp;
1018: Vec x, y;
1019: PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock *)jac->data;
1020: PC subpc;
1021: IS is;
1022: MatReuse scall;
1023: VecType vectype;
1024: MatNullSpace *nullsp_mat = NULL, *nullsp_pmat = NULL;
1026: PetscFunctionBegin;
1027: PetscCall(MatGetLocalSize(pc->pmat, &M, &N));
1029: n_local = jac->n_local;
1031: if (pc->useAmat) {
1032: PetscBool same;
1033: PetscCall(PetscObjectTypeCompare((PetscObject)mat, ((PetscObject)pmat)->type_name, &same));
1034: PetscCheck(same, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "Matrices not of same type");
1035: }
1037: if (!pc->setupcalled) {
1038: PetscInt nestlevel;
1040: scall = MAT_INITIAL_MATRIX;
1042: if (!jac->ksp) {
1043: pc->ops->reset = PCReset_BJacobi_Multiblock;
1044: pc->ops->destroy = PCDestroy_BJacobi_Multiblock;
1045: pc->ops->apply = PCApply_BJacobi_Multiblock;
1046: pc->ops->matapply = NULL;
1047: pc->ops->matapplytranspose = NULL;
1048: pc->ops->applysymmetricleft = PCApplySymmetricLeft_BJacobi_Multiblock;
1049: pc->ops->applysymmetricright = PCApplySymmetricRight_BJacobi_Multiblock;
1050: pc->ops->applytranspose = PCApplyTranspose_BJacobi_Multiblock;
1051: pc->ops->setuponblocks = PCSetUpOnBlocks_BJacobi_Multiblock;
1053: PetscCall(PetscNew(&bjac));
1054: PetscCall(PetscMalloc1(n_local, &jac->ksp));
1055: PetscCall(PetscMalloc2(n_local, &bjac->x, n_local, &bjac->y));
1056: PetscCall(PetscMalloc1(n_local, &bjac->starts));
1058: jac->data = (void *)bjac;
1059: PetscCall(PetscMalloc1(n_local, &bjac->is));
1061: for (i = 0; i < n_local; i++) {
1062: PetscCall(KSPCreate(PETSC_COMM_SELF, &ksp));
1063: PetscCall(PCGetKSPNestLevel(pc, &nestlevel));
1064: PetscCall(KSPSetNestLevel(ksp, nestlevel + 1));
1065: PetscCall(KSPSetErrorIfNotConverged(ksp, pc->erroriffailure));
1066: PetscCall(PetscObjectIncrementTabLevel((PetscObject)ksp, (PetscObject)pc, 1));
1067: PetscCall(KSPSetType(ksp, KSPPREONLY));
1068: PetscCall(KSPGetPC(ksp, &subpc));
1069: PetscCall(PCGetOptionsPrefix(pc, &prefix));
1070: PetscCall(KSPSetOptionsPrefix(ksp, prefix));
1071: PetscCall(KSPAppendOptionsPrefix(ksp, "sub_"));
1073: jac->ksp[i] = ksp;
1074: }
1075: } else {
1076: bjac = (PC_BJacobi_Multiblock *)jac->data;
1077: }
1079: start = 0;
1080: PetscCall(MatGetVecType(pmat, &vectype));
1081: for (i = 0; i < n_local; i++) {
1082: m = jac->l_lens[i];
1083: /*
1084: The reason we need to generate these vectors is to serve
1085: as the right-hand side and solution vector for the solve on the
1086: block. We do not need to allocate space for the vectors since
1087: that is provided via VecPlaceArray() just before the call to
1088: KSPSolve() on the block.
1090: */
1091: PetscCall(VecCreateSeq(PETSC_COMM_SELF, m, &x));
1092: PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, m, NULL, &y));
1093: PetscCall(VecSetType(x, vectype));
1094: PetscCall(VecSetType(y, vectype));
1096: bjac->x[i] = x;
1097: bjac->y[i] = y;
1098: bjac->starts[i] = start;
1100: PetscCall(ISCreateStride(PETSC_COMM_SELF, m, start, 1, &is));
1101: bjac->is[i] = is;
1103: start += m;
1104: }
1105: } else {
1106: bjac = (PC_BJacobi_Multiblock *)jac->data;
1107: /*
1108: Destroy the blocks from the previous iteration
1109: */
1110: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
1111: PetscCall(MatGetNullSpaces(n_local, bjac->pmat, &nullsp_pmat));
1112: PetscCall(MatDestroyMatrices(n_local, &bjac->pmat));
1113: if (pc->useAmat) {
1114: PetscCall(MatGetNullSpaces(n_local, bjac->mat, &nullsp_mat));
1115: PetscCall(MatDestroyMatrices(n_local, &bjac->mat));
1116: }
1117: scall = MAT_INITIAL_MATRIX;
1118: } else scall = MAT_REUSE_MATRIX;
1119: }
1121: PetscCall(MatCreateSubMatrices(pmat, n_local, bjac->is, bjac->is, scall, &bjac->pmat));
1122: if (nullsp_pmat) PetscCall(MatRestoreNullSpaces(n_local, bjac->pmat, &nullsp_pmat));
1123: if (pc->useAmat) {
1124: PetscCall(MatCreateSubMatrices(mat, n_local, bjac->is, bjac->is, scall, &bjac->mat));
1125: if (nullsp_mat) PetscCall(MatRestoreNullSpaces(n_local, bjac->mat, &nullsp_mat));
1126: }
1127: /* Return control to the user so that the submatrices can be modified (e.g., to apply
1128: different boundary conditions for the submatrices than for the global problem) */
1129: PetscCall(PCModifySubMatrices(pc, n_local, bjac->is, bjac->is, bjac->pmat, pc->modifysubmatricesP));
1131: for (i = 0; i < n_local; i++) {
1132: PetscCall(KSPGetOptionsPrefix(jac->ksp[i], &prefix));
1133: if (pc->useAmat) {
1134: PetscCall(KSPSetOperators(jac->ksp[i], bjac->mat[i], bjac->pmat[i]));
1135: PetscCall(MatSetOptionsPrefix(bjac->mat[i], prefix));
1136: } else {
1137: PetscCall(KSPSetOperators(jac->ksp[i], bjac->pmat[i], bjac->pmat[i]));
1138: }
1139: PetscCall(MatSetOptionsPrefix(bjac->pmat[i], prefix));
1140: if (pc->setfromoptionscalled) PetscCall(KSPSetFromOptions(jac->ksp[i]));
1141: }
1142: PetscFunctionReturn(PETSC_SUCCESS);
1143: }
1145: /*
1146: These are for a single block with multiple processes
1147: */
1148: static PetscErrorCode PCSetUpOnBlocks_BJacobi_Multiproc(PC pc)
1149: {
1150: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
1151: KSP subksp = jac->ksp[0];
1152: KSPConvergedReason reason;
1154: PetscFunctionBegin;
1155: PetscCall(KSPSetUp(subksp));
1156: PetscCall(KSPGetConvergedReason(subksp, &reason));
1157: if (reason == KSP_DIVERGED_PC_FAILED) pc->failedreason = PC_SUBPC_ERROR;
1158: PetscFunctionReturn(PETSC_SUCCESS);
1159: }
1161: static PetscErrorCode PCReset_BJacobi_Multiproc(PC pc)
1162: {
1163: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
1164: PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc *)jac->data;
1166: PetscFunctionBegin;
1167: PetscCall(VecDestroy(&mpjac->ysub));
1168: PetscCall(VecDestroy(&mpjac->xsub));
1169: PetscCall(MatDestroy(&mpjac->submats));
1170: if (jac->ksp) PetscCall(KSPReset(jac->ksp[0]));
1171: PetscFunctionReturn(PETSC_SUCCESS);
1172: }
1174: static PetscErrorCode PCDestroy_BJacobi_Multiproc(PC pc)
1175: {
1176: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
1177: PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc *)jac->data;
1179: PetscFunctionBegin;
1180: PetscCall(PCReset_BJacobi_Multiproc(pc));
1181: PetscCall(KSPDestroy(&jac->ksp[0]));
1182: PetscCall(PetscFree(jac->ksp));
1183: PetscCall(PetscSubcommDestroy(&mpjac->psubcomm));
1185: PetscCall(PetscFree(mpjac));
1186: PetscCall(PCDestroy_BJacobi(pc));
1187: PetscFunctionReturn(PETSC_SUCCESS);
1188: }
1190: static PetscErrorCode PCApply_BJacobi_Multiproc(PC pc, Vec x, Vec y)
1191: {
1192: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
1193: PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc *)jac->data;
1194: PetscScalar *yarray;
1195: const PetscScalar *xarray;
1196: KSPConvergedReason reason;
1198: PetscFunctionBegin;
1199: /* place x's and y's local arrays into xsub and ysub */
1200: PetscCall(VecGetArrayRead(x, &xarray));
1201: PetscCall(VecGetArray(y, &yarray));
1202: PetscCall(VecPlaceArray(mpjac->xsub, xarray));
1203: PetscCall(VecPlaceArray(mpjac->ysub, yarray));
1205: /* apply preconditioner on each matrix block */
1206: PetscCall(PetscLogEventBegin(PC_ApplyOnBlocks, jac->ksp[0], mpjac->xsub, mpjac->ysub, 0));
1207: PetscCall(KSPSolve(jac->ksp[0], mpjac->xsub, mpjac->ysub));
1208: PetscCall(KSPCheckSolve(jac->ksp[0], pc, mpjac->ysub));
1209: PetscCall(PetscLogEventEnd(PC_ApplyOnBlocks, jac->ksp[0], mpjac->xsub, mpjac->ysub, 0));
1210: PetscCall(KSPGetConvergedReason(jac->ksp[0], &reason));
1211: if (reason == KSP_DIVERGED_PC_FAILED) pc->failedreason = PC_SUBPC_ERROR;
1213: PetscCall(VecResetArray(mpjac->xsub));
1214: PetscCall(VecResetArray(mpjac->ysub));
1215: PetscCall(VecRestoreArrayRead(x, &xarray));
1216: PetscCall(VecRestoreArray(y, &yarray));
1217: PetscFunctionReturn(PETSC_SUCCESS);
1218: }
1220: static PetscErrorCode PCMatApply_BJacobi_Multiproc(PC pc, Mat X, Mat Y)
1221: {
1222: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
1223: KSPConvergedReason reason;
1224: Mat sX, sY;
1225: const PetscScalar *x;
1226: PetscScalar *y;
1227: PetscInt m, N, lda, ldb;
1229: PetscFunctionBegin;
1230: /* apply preconditioner on each matrix block */
1231: PetscCall(MatGetLocalSize(X, &m, NULL));
1232: PetscCall(MatGetSize(X, NULL, &N));
1233: PetscCall(MatDenseGetLDA(X, &lda));
1234: PetscCall(MatDenseGetLDA(Y, &ldb));
1235: PetscCall(MatDenseGetArrayRead(X, &x));
1236: PetscCall(MatDenseGetArrayWrite(Y, &y));
1237: PetscCall(MatCreateDense(PetscObjectComm((PetscObject)jac->ksp[0]), m, PETSC_DECIDE, PETSC_DECIDE, N, (PetscScalar *)x, &sX));
1238: PetscCall(MatCreateDense(PetscObjectComm((PetscObject)jac->ksp[0]), m, PETSC_DECIDE, PETSC_DECIDE, N, y, &sY));
1239: PetscCall(MatDenseSetLDA(sX, lda));
1240: PetscCall(MatDenseSetLDA(sY, ldb));
1241: PetscCall(PetscLogEventBegin(PC_ApplyOnBlocks, jac->ksp[0], X, Y, 0));
1242: PetscCall(KSPMatSolve(jac->ksp[0], sX, sY));
1243: PetscCall(KSPCheckSolve(jac->ksp[0], pc, NULL));
1244: PetscCall(PetscLogEventEnd(PC_ApplyOnBlocks, jac->ksp[0], X, Y, 0));
1245: PetscCall(MatDestroy(&sY));
1246: PetscCall(MatDestroy(&sX));
1247: PetscCall(MatDenseRestoreArrayWrite(Y, &y));
1248: PetscCall(MatDenseRestoreArrayRead(X, &x));
1249: PetscCall(KSPGetConvergedReason(jac->ksp[0], &reason));
1250: if (reason == KSP_DIVERGED_PC_FAILED) pc->failedreason = PC_SUBPC_ERROR;
1251: PetscFunctionReturn(PETSC_SUCCESS);
1252: }
1254: static PetscErrorCode PCSetUp_BJacobi_Multiproc(PC pc)
1255: {
1256: PC_BJacobi *jac = (PC_BJacobi *)pc->data;
1257: PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc *)jac->data;
1258: PetscInt m, n;
1259: MPI_Comm comm, subcomm = 0;
1260: const char *prefix;
1261: PetscBool wasSetup = PETSC_TRUE;
1262: VecType vectype;
1264: PetscFunctionBegin;
1265: PetscCall(PetscObjectGetComm((PetscObject)pc, &comm));
1266: PetscCheck(jac->n_local <= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Only a single block in a subcommunicator is supported");
1267: jac->n_local = 1; /* currently only a single block is supported for a subcommunicator */
1268: if (!pc->setupcalled) {
1269: PetscInt nestlevel;
1271: wasSetup = PETSC_FALSE;
1272: PetscCall(PetscNew(&mpjac));
1273: jac->data = (void *)mpjac;
1275: /* initialize datastructure mpjac */
1276: if (!jac->psubcomm) {
1277: /* Create default contiguous subcommunicatiors if user does not provide them */
1278: PetscCall(PetscSubcommCreate(comm, &jac->psubcomm));
1279: PetscCall(PetscSubcommSetNumber(jac->psubcomm, jac->n));
1280: PetscCall(PetscSubcommSetType(jac->psubcomm, PETSC_SUBCOMM_CONTIGUOUS));
1281: }
1282: mpjac->psubcomm = jac->psubcomm;
1283: subcomm = PetscSubcommChild(mpjac->psubcomm);
1285: /* Get matrix blocks of pmat */
1286: PetscCall(MatGetMultiProcBlock(pc->pmat, subcomm, MAT_INITIAL_MATRIX, &mpjac->submats));
1288: /* create a new PC that processors in each subcomm have copy of */
1289: PetscCall(PetscMalloc1(1, &jac->ksp));
1290: PetscCall(KSPCreate(subcomm, &jac->ksp[0]));
1291: PetscCall(PCGetKSPNestLevel(pc, &nestlevel));
1292: PetscCall(KSPSetNestLevel(jac->ksp[0], nestlevel + 1));
1293: PetscCall(KSPSetErrorIfNotConverged(jac->ksp[0], pc->erroriffailure));
1294: PetscCall(PetscObjectIncrementTabLevel((PetscObject)jac->ksp[0], (PetscObject)pc, 1));
1295: PetscCall(KSPSetOperators(jac->ksp[0], mpjac->submats, mpjac->submats));
1296: PetscCall(KSPGetPC(jac->ksp[0], &mpjac->pc));
1298: PetscCall(PCGetOptionsPrefix(pc, &prefix));
1299: PetscCall(KSPSetOptionsPrefix(jac->ksp[0], prefix));
1300: PetscCall(KSPAppendOptionsPrefix(jac->ksp[0], "sub_"));
1301: PetscCall(KSPGetOptionsPrefix(jac->ksp[0], &prefix));
1302: PetscCall(MatSetOptionsPrefix(mpjac->submats, prefix));
1304: /* create dummy vectors xsub and ysub */
1305: PetscCall(MatGetLocalSize(mpjac->submats, &m, &n));
1306: PetscCall(VecCreateMPIWithArray(subcomm, 1, n, PETSC_DECIDE, NULL, &mpjac->xsub));
1307: PetscCall(VecCreateMPIWithArray(subcomm, 1, m, PETSC_DECIDE, NULL, &mpjac->ysub));
1308: PetscCall(MatGetVecType(mpjac->submats, &vectype));
1309: PetscCall(VecSetType(mpjac->xsub, vectype));
1310: PetscCall(VecSetType(mpjac->ysub, vectype));
1312: pc->ops->setuponblocks = PCSetUpOnBlocks_BJacobi_Multiproc;
1313: pc->ops->reset = PCReset_BJacobi_Multiproc;
1314: pc->ops->destroy = PCDestroy_BJacobi_Multiproc;
1315: pc->ops->apply = PCApply_BJacobi_Multiproc;
1316: pc->ops->matapply = PCMatApply_BJacobi_Multiproc;
1317: } else { /* pc->setupcalled */
1318: subcomm = PetscSubcommChild(mpjac->psubcomm);
1319: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
1320: /* destroy old matrix blocks, then get new matrix blocks */
1321: if (mpjac->submats) PetscCall(MatDestroy(&mpjac->submats));
1322: PetscCall(MatGetMultiProcBlock(pc->pmat, subcomm, MAT_INITIAL_MATRIX, &mpjac->submats));
1323: } else {
1324: PetscCall(MatGetMultiProcBlock(pc->pmat, subcomm, MAT_REUSE_MATRIX, &mpjac->submats));
1325: }
1326: PetscCall(KSPSetOperators(jac->ksp[0], mpjac->submats, mpjac->submats));
1327: }
1329: if (!wasSetup && pc->setfromoptionscalled) PetscCall(KSPSetFromOptions(jac->ksp[0]));
1330: PetscFunctionReturn(PETSC_SUCCESS);
1331: }