Actual source code: bjacobi.c
petsc-3.14.6 2021-03-30
2: /*
3: Defines a block Jacobi preconditioner.
4: */
6: #include <../src/ksp/pc/impls/bjacobi/bjacobi.h>
8: static PetscErrorCode PCSetUp_BJacobi_Singleblock(PC,Mat,Mat);
9: static PetscErrorCode PCSetUp_BJacobi_Multiblock(PC,Mat,Mat);
10: static PetscErrorCode PCSetUp_BJacobi_Multiproc(PC);
12: static PetscErrorCode PCSetUp_BJacobi(PC pc)
13: {
14: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
15: Mat mat = pc->mat,pmat = pc->pmat;
17: PetscBool hasop;
18: PetscInt N,M,start,i,sum,end;
19: PetscInt bs,i_start=-1,i_end=-1;
20: PetscMPIInt rank,size;
21: const char *pprefix,*mprefix;
24: MPI_Comm_rank(PetscObjectComm((PetscObject)pc),&rank);
25: MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);
26: MatGetLocalSize(pc->pmat,&M,&N);
27: MatGetBlockSize(pc->pmat,&bs);
29: if (jac->n > 0 && jac->n < size) {
30: PCSetUp_BJacobi_Multiproc(pc);
31: return(0);
32: }
34: /* --------------------------------------------------------------------------
35: Determines the number of blocks assigned to each processor
36: -----------------------------------------------------------------------------*/
38: /* local block count given */
39: if (jac->n_local > 0 && jac->n < 0) {
40: MPIU_Allreduce(&jac->n_local,&jac->n,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)pc));
41: if (jac->l_lens) { /* check that user set these correctly */
42: sum = 0;
43: for (i=0; i<jac->n_local; i++) {
44: if (jac->l_lens[i]/bs*bs !=jac->l_lens[i]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Mat blocksize doesn't match block Jacobi layout");
45: sum += jac->l_lens[i];
46: }
47: if (sum != M) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local lens set incorrectly");
48: } else {
49: PetscMalloc1(jac->n_local,&jac->l_lens);
50: for (i=0; i<jac->n_local; i++) jac->l_lens[i] = bs*((M/bs)/jac->n_local + (((M/bs) % jac->n_local) > i));
51: }
52: } else if (jac->n > 0 && jac->n_local < 0) { /* global block count given */
53: /* global blocks given: determine which ones are local */
54: if (jac->g_lens) {
55: /* check if the g_lens is has valid entries */
56: for (i=0; i<jac->n; i++) {
57: if (!jac->g_lens[i]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Zero block not allowed");
58: if (jac->g_lens[i]/bs*bs != jac->g_lens[i]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Mat blocksize doesn't match block Jacobi layout");
59: }
60: if (size == 1) {
61: jac->n_local = jac->n;
62: PetscMalloc1(jac->n_local,&jac->l_lens);
63: PetscArraycpy(jac->l_lens,jac->g_lens,jac->n_local);
64: /* check that user set these correctly */
65: sum = 0;
66: for (i=0; i<jac->n_local; i++) sum += jac->l_lens[i];
67: if (sum != M) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Global lens set incorrectly");
68: } else {
69: MatGetOwnershipRange(pc->pmat,&start,&end);
70: /* loop over blocks determing first one owned by me */
71: sum = 0;
72: for (i=0; i<jac->n+1; i++) {
73: if (sum == start) { i_start = i; goto start_1;}
74: if (i < jac->n) sum += jac->g_lens[i];
75: }
76: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Block sizes used in PCBJacobiSetTotalBlocks()\nare not compatible with parallel matrix layout");
77: start_1:
78: for (i=i_start; i<jac->n+1; i++) {
79: if (sum == end) { i_end = i; goto end_1; }
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: PetscMalloc1(jac->n_local,&jac->l_lens);
86: 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: 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: if (!jac->l_lens[i]) SETERRQ(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: 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: 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: if (jac->n_local < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Number of blocks is less than number of processors");
109: /* -------------------------
110: Determines mat and pmat
111: ---------------------------*/
112: MatHasOperation(pc->mat,MATOP_GET_DIAGONAL_BLOCK,&hasop);
113: if (!hasop && size == 1) {
114: mat = pc->mat;
115: pmat = pc->pmat;
116: } else {
117: if (pc->useAmat) {
118: /* use block from Amat matrix, not Pmat for local MatMult() */
119: MatGetDiagonalBlock(pc->mat,&mat);
120: /* make submatrix have same prefix as entire matrix */
121: PetscObjectGetOptionsPrefix((PetscObject)pc->mat,&mprefix);
122: PetscObjectSetOptionsPrefix((PetscObject)mat,mprefix);
123: }
124: if (pc->pmat != pc->mat || !pc->useAmat) {
125: MatGetDiagonalBlock(pc->pmat,&pmat);
126: /* make submatrix have same prefix as entire matrix */
127: PetscObjectGetOptionsPrefix((PetscObject)pc->pmat,&pprefix);
128: PetscObjectSetOptionsPrefix((PetscObject)pmat,pprefix);
129: } else pmat = mat;
130: }
132: /* ------
133: Setup code depends on the number of blocks
134: */
135: if (jac->n_local == 1) {
136: PCSetUp_BJacobi_Singleblock(pc,mat,pmat);
137: } else {
138: PCSetUp_BJacobi_Multiblock(pc,mat,pmat);
139: }
140: return(0);
141: }
143: /* Default destroy, if it has never been setup */
144: static PetscErrorCode PCDestroy_BJacobi(PC pc)
145: {
146: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
150: PetscFree(jac->g_lens);
151: PetscFree(jac->l_lens);
152: PetscFree(pc->data);
153: return(0);
154: }
157: static PetscErrorCode PCSetFromOptions_BJacobi(PetscOptionItems *PetscOptionsObject,PC pc)
158: {
159: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
161: PetscInt blocks,i;
162: PetscBool flg;
165: PetscOptionsHead(PetscOptionsObject,"Block Jacobi options");
166: PetscOptionsInt("-pc_bjacobi_blocks","Total number of blocks","PCBJacobiSetTotalBlocks",jac->n,&blocks,&flg);
167: if (flg) {PCBJacobiSetTotalBlocks(pc,blocks,NULL);}
168: PetscOptionsInt("-pc_bjacobi_local_blocks","Local number of blocks","PCBJacobiSetLocalBlocks",jac->n_local,&blocks,&flg);
169: if (flg) {PCBJacobiSetLocalBlocks(pc,blocks,NULL);}
170: if (jac->ksp) {
171: /* The sub-KSP has already been set up (e.g., PCSetUp_BJacobi_Singleblock), but KSPSetFromOptions was not called
172: * unless we had already been called. */
173: for (i=0; i<jac->n_local; i++) {
174: KSPSetFromOptions(jac->ksp[i]);
175: }
176: }
177: PetscOptionsTail();
178: return(0);
179: }
181: #include <petscdraw.h>
182: static PetscErrorCode PCView_BJacobi(PC pc,PetscViewer viewer)
183: {
184: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
185: PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc*)jac->data;
186: PetscErrorCode ierr;
187: PetscMPIInt rank;
188: PetscInt i;
189: PetscBool iascii,isstring,isdraw;
190: PetscViewer sviewer;
193: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
194: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
195: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
196: if (iascii) {
197: if (pc->useAmat) {
198: PetscViewerASCIIPrintf(viewer," using Amat local matrix, number of blocks = %D\n",jac->n);
199: }
200: PetscViewerASCIIPrintf(viewer," number of blocks = %D\n",jac->n);
201: MPI_Comm_rank(PetscObjectComm((PetscObject)pc),&rank);
202: if (jac->same_local_solves) {
203: PetscViewerASCIIPrintf(viewer," Local solver is the same for all blocks, as in the following KSP and PC objects on rank 0:\n");
204: if (jac->ksp && !jac->psubcomm) {
205: PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
206: if (!rank) {
207: PetscViewerASCIIPushTab(viewer);
208: KSPView(jac->ksp[0],sviewer);
209: PetscViewerASCIIPopTab(viewer);
210: }
211: PetscViewerFlush(sviewer);
212: PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
213: PetscViewerFlush(viewer);
214: /* extra call needed because of the two calls to PetscViewerASCIIPushSynchronized() in PetscViewerGetSubViewer() */
215: PetscViewerASCIIPopSynchronized(viewer);
216: } else if (mpjac && jac->ksp && mpjac->psubcomm) {
217: PetscViewerGetSubViewer(viewer,mpjac->psubcomm->child,&sviewer);
218: if (!mpjac->psubcomm->color) {
219: PetscViewerASCIIPushTab(viewer);
220: KSPView(*(jac->ksp),sviewer);
221: PetscViewerASCIIPopTab(viewer);
222: }
223: PetscViewerFlush(sviewer);
224: PetscViewerRestoreSubViewer(viewer,mpjac->psubcomm->child,&sviewer);
225: PetscViewerFlush(viewer);
226: /* extra call needed because of the two calls to PetscViewerASCIIPushSynchronized() in PetscViewerGetSubViewer() */
227: PetscViewerASCIIPopSynchronized(viewer);
228: } else {
229: PetscViewerFlush(viewer);
230: }
231: } else {
232: PetscInt n_global;
233: MPIU_Allreduce(&jac->n_local,&n_global,1,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)pc));
234: PetscViewerASCIIPushSynchronized(viewer);
235: PetscViewerASCIIPrintf(viewer," Local solve info for each block is in the following KSP and PC objects:\n");
236: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] number of local blocks = %D, first local block number = %D\n",
237: rank,jac->n_local,jac->first_local);
238: PetscViewerASCIIPushTab(viewer);
239: PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
240: for (i=0; i<jac->n_local; i++) {
241: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] local block number %D\n",rank,i);
242: KSPView(jac->ksp[i],sviewer);
243: PetscViewerASCIISynchronizedPrintf(viewer,"- - - - - - - - - - - - - - - - - -\n");
244: }
245: PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
246: PetscViewerASCIIPopTab(viewer);
247: PetscViewerFlush(viewer);
248: PetscViewerASCIIPopSynchronized(viewer);
249: }
250: } else if (isstring) {
251: PetscViewerStringSPrintf(viewer," blks=%D",jac->n);
252: PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
253: if (jac->ksp) {KSPView(jac->ksp[0],sviewer);}
254: PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
255: } else if (isdraw) {
256: PetscDraw draw;
257: char str[25];
258: PetscReal x,y,bottom,h;
260: PetscViewerDrawGetDraw(viewer,0,&draw);
261: PetscDrawGetCurrentPoint(draw,&x,&y);
262: PetscSNPrintf(str,25,"Number blocks %D",jac->n);
263: PetscDrawStringBoxed(draw,x,y,PETSC_DRAW_RED,PETSC_DRAW_BLACK,str,NULL,&h);
264: bottom = y - h;
265: PetscDrawPushCurrentPoint(draw,x,bottom);
266: /* warning the communicator on viewer is different then on ksp in parallel */
267: if (jac->ksp) {KSPView(jac->ksp[0],viewer);}
268: PetscDrawPopCurrentPoint(draw);
269: }
270: return(0);
271: }
273: /* -------------------------------------------------------------------------------------*/
275: static PetscErrorCode PCBJacobiGetSubKSP_BJacobi(PC pc,PetscInt *n_local,PetscInt *first_local,KSP **ksp)
276: {
277: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
280: if (!pc->setupcalled) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must call KSPSetUp() or PCSetUp() first");
282: if (n_local) *n_local = jac->n_local;
283: if (first_local) *first_local = jac->first_local;
284: *ksp = jac->ksp;
285: jac->same_local_solves = PETSC_FALSE; /* Assume that local solves are now different;
286: not necessarily true though! This flag is
287: used only for PCView_BJacobi() */
288: return(0);
289: }
291: static PetscErrorCode PCBJacobiSetTotalBlocks_BJacobi(PC pc,PetscInt blocks,PetscInt *lens)
292: {
293: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
297: if (pc->setupcalled > 0 && jac->n!=blocks) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ORDER,"Cannot alter number of blocks after PCSetUp()/KSPSetUp() has been called");
298: jac->n = blocks;
299: if (!lens) jac->g_lens = NULL;
300: else {
301: PetscMalloc1(blocks,&jac->g_lens);
302: PetscLogObjectMemory((PetscObject)pc,blocks*sizeof(PetscInt));
303: PetscArraycpy(jac->g_lens,lens,blocks);
304: }
305: return(0);
306: }
308: static PetscErrorCode PCBJacobiGetTotalBlocks_BJacobi(PC pc, PetscInt *blocks, const PetscInt *lens[])
309: {
310: PC_BJacobi *jac = (PC_BJacobi*) pc->data;
313: *blocks = jac->n;
314: if (lens) *lens = jac->g_lens;
315: return(0);
316: }
318: static PetscErrorCode PCBJacobiSetLocalBlocks_BJacobi(PC pc,PetscInt blocks,const PetscInt lens[])
319: {
320: PC_BJacobi *jac;
324: jac = (PC_BJacobi*)pc->data;
326: jac->n_local = blocks;
327: if (!lens) jac->l_lens = NULL;
328: else {
329: PetscMalloc1(blocks,&jac->l_lens);
330: PetscLogObjectMemory((PetscObject)pc,blocks*sizeof(PetscInt));
331: PetscArraycpy(jac->l_lens,lens,blocks);
332: }
333: return(0);
334: }
336: static PetscErrorCode PCBJacobiGetLocalBlocks_BJacobi(PC pc, PetscInt *blocks, const PetscInt *lens[])
337: {
338: PC_BJacobi *jac = (PC_BJacobi*) pc->data;
341: *blocks = jac->n_local;
342: if (lens) *lens = jac->l_lens;
343: return(0);
344: }
346: /* -------------------------------------------------------------------------------------*/
348: /*@C
349: PCBJacobiGetSubKSP - Gets the local KSP contexts for all blocks on
350: this processor.
352: Not Collective
354: Input Parameter:
355: . pc - the preconditioner context
357: Output Parameters:
358: + n_local - the number of blocks on this processor, or NULL
359: . first_local - the global number of the first block on this processor, or NULL
360: - ksp - the array of KSP contexts
362: Notes:
363: After PCBJacobiGetSubKSP() the array of KSP contexts is not to be freed.
365: Currently for some matrix implementations only 1 block per processor
366: is supported.
368: You must call KSPSetUp() or PCSetUp() before calling PCBJacobiGetSubKSP().
370: Fortran Usage: You must pass in a KSP array that is large enough to contain all the local KSPs.
371: You can call PCBJacobiGetSubKSP(pc,nlocal,firstlocal,PETSC_NULL_KSP,ierr) to determine how large the
372: KSP array must be.
374: Level: advanced
376: .seealso: PCBJacobiGetSubKSP()
377: @*/
378: PetscErrorCode PCBJacobiGetSubKSP(PC pc,PetscInt *n_local,PetscInt *first_local,KSP *ksp[])
379: {
384: PetscUseMethod(pc,"PCBJacobiGetSubKSP_C",(PC,PetscInt*,PetscInt*,KSP **),(pc,n_local,first_local,ksp));
385: return(0);
386: }
388: /*@
389: PCBJacobiSetTotalBlocks - Sets the global number of blocks for the block
390: Jacobi preconditioner.
392: Collective on PC
394: Input Parameters:
395: + pc - the preconditioner context
396: . blocks - the number of blocks
397: - lens - [optional] integer array containing the size of each block
399: Options Database Key:
400: . -pc_bjacobi_blocks <blocks> - Sets the number of global blocks
402: Notes:
403: Currently only a limited number of blocking configurations are supported.
404: All processors sharing the PC must call this routine with the same data.
406: Level: intermediate
408: .seealso: PCSetUseAmat(), PCBJacobiSetLocalBlocks()
409: @*/
410: PetscErrorCode PCBJacobiSetTotalBlocks(PC pc,PetscInt blocks,const PetscInt lens[])
411: {
416: if (blocks <= 0) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Must have positive blocks");
417: PetscTryMethod(pc,"PCBJacobiSetTotalBlocks_C",(PC,PetscInt,const PetscInt[]),(pc,blocks,lens));
418: return(0);
419: }
421: /*@C
422: PCBJacobiGetTotalBlocks - Gets the global number of blocks for the block
423: Jacobi preconditioner.
425: Not Collective
427: Input Parameter:
428: . pc - the preconditioner context
430: Output parameters:
431: + blocks - the number of blocks
432: - lens - integer array containing the size of each block
434: Level: intermediate
436: .seealso: PCSetUseAmat(), PCBJacobiGetLocalBlocks()
437: @*/
438: PetscErrorCode PCBJacobiGetTotalBlocks(PC pc, PetscInt *blocks, const PetscInt *lens[])
439: {
445: PetscUseMethod(pc,"PCBJacobiGetTotalBlocks_C",(PC,PetscInt*, const PetscInt *[]),(pc,blocks,lens));
446: return(0);
447: }
449: /*@
450: PCBJacobiSetLocalBlocks - Sets the local number of blocks for the block
451: Jacobi preconditioner.
453: Not Collective
455: Input Parameters:
456: + pc - the preconditioner context
457: . blocks - the number of blocks
458: - lens - [optional] integer array containing size of each block
460: Options Database Key:
461: . -pc_bjacobi_local_blocks <blocks> - Sets the number of local blocks
463: Note:
464: Currently only a limited number of blocking configurations are supported.
466: Level: intermediate
468: .seealso: PCSetUseAmat(), PCBJacobiSetTotalBlocks()
469: @*/
470: PetscErrorCode PCBJacobiSetLocalBlocks(PC pc,PetscInt blocks,const PetscInt lens[])
471: {
476: if (blocks < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must have nonegative blocks");
477: PetscTryMethod(pc,"PCBJacobiSetLocalBlocks_C",(PC,PetscInt,const PetscInt []),(pc,blocks,lens));
478: return(0);
479: }
481: /*@C
482: PCBJacobiGetLocalBlocks - Gets the local number of blocks for the block
483: Jacobi preconditioner.
485: Not Collective
487: Input Parameters:
488: + pc - the preconditioner context
489: . blocks - the number of blocks
490: - lens - [optional] integer array containing size of each block
492: Note:
493: Currently only a limited number of blocking configurations are supported.
495: Level: intermediate
497: .seealso: PCSetUseAmat(), PCBJacobiGetTotalBlocks()
498: @*/
499: PetscErrorCode PCBJacobiGetLocalBlocks(PC pc, PetscInt *blocks, const PetscInt *lens[])
500: {
506: PetscUseMethod(pc,"PCBJacobiGetLocalBlocks_C",(PC,PetscInt*, const PetscInt *[]),(pc,blocks,lens));
507: return(0);
508: }
510: /* -----------------------------------------------------------------------------------*/
512: /*MC
513: PCBJACOBI - Use block Jacobi preconditioning, each block is (approximately) solved with
514: its own KSP object.
516: Options Database Keys:
517: + -pc_use_amat - use Amat to apply block of operator in inner Krylov method
518: - -pc_bjacobi_blocks <n> - use n total blocks
520: Notes:
521: Each processor can have one or more blocks, or a single block can be shared by several processes. Defaults to one block per processor.
523: To set options on the solvers for each block append -sub_ to all the KSP, KSP, and PC
524: options database keys. For example, -sub_pc_type ilu -sub_pc_factor_levels 1 -sub_ksp_type preonly
526: To set the options on the solvers separate for each block call PCBJacobiGetSubKSP()
527: and set the options directly on the resulting KSP object (you can access its PC
528: KSPGetPC())
530: For GPU-based vectors (CUDA, ViennaCL) it is recommended to use exactly one block per MPI process for best
531: performance. Different block partitioning may lead to additional data transfers
532: between host and GPU that lead to degraded performance.
534: The options prefix for each block is sub_, for example -sub_pc_type lu.
536: When multiple processes share a single block, each block encompasses exactly all the unknowns owned its set of processes.
538: Level: beginner
540: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
541: PCASM, PCSetUseAmat(), PCGetUseAmat(), PCBJacobiGetSubKSP(), PCBJacobiSetTotalBlocks(),
542: PCBJacobiSetLocalBlocks(), PCSetModifySubMatrices()
543: M*/
545: PETSC_EXTERN PetscErrorCode PCCreate_BJacobi(PC pc)
546: {
548: PetscMPIInt rank;
549: PC_BJacobi *jac;
552: PetscNewLog(pc,&jac);
553: MPI_Comm_rank(PetscObjectComm((PetscObject)pc),&rank);
555: pc->ops->apply = NULL;
556: pc->ops->matapply = NULL;
557: pc->ops->applytranspose = NULL;
558: pc->ops->setup = PCSetUp_BJacobi;
559: pc->ops->destroy = PCDestroy_BJacobi;
560: pc->ops->setfromoptions = PCSetFromOptions_BJacobi;
561: pc->ops->view = PCView_BJacobi;
562: pc->ops->applyrichardson = NULL;
564: pc->data = (void*)jac;
565: jac->n = -1;
566: jac->n_local = -1;
567: jac->first_local = rank;
568: jac->ksp = NULL;
569: jac->same_local_solves = PETSC_TRUE;
570: jac->g_lens = NULL;
571: jac->l_lens = NULL;
572: jac->psubcomm = NULL;
574: PetscObjectComposeFunction((PetscObject)pc,"PCBJacobiGetSubKSP_C",PCBJacobiGetSubKSP_BJacobi);
575: PetscObjectComposeFunction((PetscObject)pc,"PCBJacobiSetTotalBlocks_C",PCBJacobiSetTotalBlocks_BJacobi);
576: PetscObjectComposeFunction((PetscObject)pc,"PCBJacobiGetTotalBlocks_C",PCBJacobiGetTotalBlocks_BJacobi);
577: PetscObjectComposeFunction((PetscObject)pc,"PCBJacobiSetLocalBlocks_C",PCBJacobiSetLocalBlocks_BJacobi);
578: PetscObjectComposeFunction((PetscObject)pc,"PCBJacobiGetLocalBlocks_C",PCBJacobiGetLocalBlocks_BJacobi);
579: return(0);
580: }
582: /* --------------------------------------------------------------------------------------------*/
583: /*
584: These are for a single block per processor; works for AIJ, BAIJ; Seq and MPI
585: */
586: static PetscErrorCode PCReset_BJacobi_Singleblock(PC pc)
587: {
588: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
589: PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock*)jac->data;
590: PetscErrorCode ierr;
593: KSPReset(jac->ksp[0]);
594: VecDestroy(&bjac->x);
595: VecDestroy(&bjac->y);
596: return(0);
597: }
599: static PetscErrorCode PCDestroy_BJacobi_Singleblock(PC pc)
600: {
601: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
602: PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock*)jac->data;
603: PetscErrorCode ierr;
606: PCReset_BJacobi_Singleblock(pc);
607: KSPDestroy(&jac->ksp[0]);
608: PetscFree(jac->ksp);
609: PetscFree(jac->l_lens);
610: PetscFree(jac->g_lens);
611: PetscFree(bjac);
612: PetscFree(pc->data);
613: return(0);
614: }
616: static PetscErrorCode PCSetUpOnBlocks_BJacobi_Singleblock(PC pc)
617: {
618: PetscErrorCode ierr;
619: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
620: KSP subksp = jac->ksp[0];
621: KSPConvergedReason reason;
624: KSPSetUp(subksp);
625: KSPGetConvergedReason(subksp,&reason);
626: if (reason == KSP_DIVERGED_PC_FAILED) {
627: pc->failedreason = PC_SUBPC_ERROR;
628: }
629: return(0);
630: }
632: static PetscErrorCode PCApply_BJacobi_Singleblock(PC pc,Vec x,Vec y)
633: {
634: PetscErrorCode ierr;
635: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
636: PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock*)jac->data;
639: VecGetLocalVectorRead(x, bjac->x);
640: VecGetLocalVector(y, bjac->y);
641: /* Since the inner KSP matrix may point directly to the diagonal block of an MPI matrix the inner
642: matrix may change even if the outer KSP/PC has not updated the preconditioner, this will trigger a rebuild
643: of the inner preconditioner automatically unless we pass down the outer preconditioners reuse flag.*/
644: KSPSetReusePreconditioner(jac->ksp[0],pc->reusepreconditioner);
645: KSPSolve(jac->ksp[0],bjac->x,bjac->y);
646: KSPCheckSolve(jac->ksp[0],pc,bjac->y);
647: VecRestoreLocalVectorRead(x, bjac->x);
648: VecRestoreLocalVector(y, bjac->y);
649: return(0);
650: }
652: static PetscErrorCode PCMatApply_BJacobi_Singleblock(PC pc,Mat X,Mat Y)
653: {
654: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
655: Mat sX,sY;
659: /* Since the inner KSP matrix may point directly to the diagonal block of an MPI matrix the inner
660: matrix may change even if the outer KSP/PC has not updated the preconditioner, this will trigger a rebuild
661: of the inner preconditioner automatically unless we pass down the outer preconditioners reuse flag.*/
662: KSPSetReusePreconditioner(jac->ksp[0],pc->reusepreconditioner);
663: MatDenseGetLocalMatrix(X,&sX);
664: MatDenseGetLocalMatrix(Y,&sY);
665: KSPMatSolve(jac->ksp[0],sX,sY);
666: return(0);
667: }
669: static PetscErrorCode PCApplySymmetricLeft_BJacobi_Singleblock(PC pc,Vec x,Vec y)
670: {
671: PetscErrorCode ierr;
672: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
673: PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock*)jac->data;
674: PetscScalar *y_array;
675: const PetscScalar *x_array;
676: PC subpc;
679: /*
680: The VecPlaceArray() is to avoid having to copy the
681: y vector into the bjac->x vector. The reason for
682: the bjac->x vector is that we need a sequential vector
683: for the sequential solve.
684: */
685: VecGetArrayRead(x,&x_array);
686: VecGetArray(y,&y_array);
687: VecPlaceArray(bjac->x,x_array);
688: VecPlaceArray(bjac->y,y_array);
689: /* apply the symmetric left portion of the inner PC operator */
690: /* note this by-passes the inner KSP and its options completely */
691: KSPGetPC(jac->ksp[0],&subpc);
692: PCApplySymmetricLeft(subpc,bjac->x,bjac->y);
693: VecResetArray(bjac->x);
694: VecResetArray(bjac->y);
695: VecRestoreArrayRead(x,&x_array);
696: VecRestoreArray(y,&y_array);
697: return(0);
698: }
700: static PetscErrorCode PCApplySymmetricRight_BJacobi_Singleblock(PC pc,Vec x,Vec y)
701: {
702: PetscErrorCode ierr;
703: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
704: PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock*)jac->data;
705: PetscScalar *y_array;
706: const PetscScalar *x_array;
707: PC subpc;
710: /*
711: The VecPlaceArray() is to avoid having to copy the
712: y vector into the bjac->x vector. The reason for
713: the bjac->x vector is that we need a sequential vector
714: for the sequential solve.
715: */
716: VecGetArrayRead(x,&x_array);
717: VecGetArray(y,&y_array);
718: VecPlaceArray(bjac->x,x_array);
719: VecPlaceArray(bjac->y,y_array);
721: /* apply the symmetric right portion of the inner PC operator */
722: /* note this by-passes the inner KSP and its options completely */
724: KSPGetPC(jac->ksp[0],&subpc);
725: PCApplySymmetricRight(subpc,bjac->x,bjac->y);
727: VecRestoreArrayRead(x,&x_array);
728: VecRestoreArray(y,&y_array);
729: return(0);
730: }
732: static PetscErrorCode PCApplyTranspose_BJacobi_Singleblock(PC pc,Vec x,Vec y)
733: {
734: PetscErrorCode ierr;
735: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
736: PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock*)jac->data;
737: PetscScalar *y_array;
738: const PetscScalar *x_array;
741: /*
742: The VecPlaceArray() is to avoid having to copy the
743: y vector into the bjac->x vector. The reason for
744: the bjac->x vector is that we need a sequential vector
745: for the sequential solve.
746: */
747: VecGetArrayRead(x,&x_array);
748: VecGetArray(y,&y_array);
749: VecPlaceArray(bjac->x,x_array);
750: VecPlaceArray(bjac->y,y_array);
751: KSPSolveTranspose(jac->ksp[0],bjac->x,bjac->y);
752: KSPCheckSolve(jac->ksp[0],pc,bjac->y);
753: VecResetArray(bjac->x);
754: VecResetArray(bjac->y);
755: VecRestoreArrayRead(x,&x_array);
756: VecRestoreArray(y,&y_array);
757: return(0);
758: }
760: static PetscErrorCode PCSetUp_BJacobi_Singleblock(PC pc,Mat mat,Mat pmat)
761: {
762: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
763: PetscErrorCode ierr;
764: PetscInt m;
765: KSP ksp;
766: PC_BJacobi_Singleblock *bjac;
767: PetscBool wasSetup = PETSC_TRUE;
768: VecType vectype;
771: if (!pc->setupcalled) {
772: const char *prefix;
774: if (!jac->ksp) {
775: wasSetup = PETSC_FALSE;
777: KSPCreate(PETSC_COMM_SELF,&ksp);
778: KSPSetErrorIfNotConverged(ksp,pc->erroriffailure);
779: PetscObjectIncrementTabLevel((PetscObject)ksp,(PetscObject)pc,1);
780: PetscLogObjectParent((PetscObject)pc,(PetscObject)ksp);
781: KSPSetType(ksp,KSPPREONLY);
782: PCGetOptionsPrefix(pc,&prefix);
783: KSPSetOptionsPrefix(ksp,prefix);
784: KSPAppendOptionsPrefix(ksp,"sub_");
786: pc->ops->reset = PCReset_BJacobi_Singleblock;
787: pc->ops->destroy = PCDestroy_BJacobi_Singleblock;
788: pc->ops->apply = PCApply_BJacobi_Singleblock;
789: pc->ops->matapply = PCMatApply_BJacobi_Singleblock;
790: pc->ops->applysymmetricleft = PCApplySymmetricLeft_BJacobi_Singleblock;
791: pc->ops->applysymmetricright = PCApplySymmetricRight_BJacobi_Singleblock;
792: pc->ops->applytranspose = PCApplyTranspose_BJacobi_Singleblock;
793: pc->ops->setuponblocks = PCSetUpOnBlocks_BJacobi_Singleblock;
795: PetscMalloc1(1,&jac->ksp);
796: jac->ksp[0] = ksp;
798: PetscNewLog(pc,&bjac);
799: jac->data = (void*)bjac;
800: } else {
801: ksp = jac->ksp[0];
802: bjac = (PC_BJacobi_Singleblock*)jac->data;
803: }
805: /*
806: The reason we need to generate these vectors is to serve
807: as the right-hand side and solution vector for the solve on the
808: block. We do not need to allocate space for the vectors since
809: that is provided via VecPlaceArray() just before the call to
810: KSPSolve() on the block.
811: */
812: MatGetSize(pmat,&m,&m);
813: VecCreateSeqWithArray(PETSC_COMM_SELF,1,m,NULL,&bjac->x);
814: VecCreateSeqWithArray(PETSC_COMM_SELF,1,m,NULL,&bjac->y);
815: MatGetVecType(pmat,&vectype);
816: VecSetType(bjac->x,vectype);
817: VecSetType(bjac->y,vectype);
818: PetscLogObjectParent((PetscObject)pc,(PetscObject)bjac->x);
819: PetscLogObjectParent((PetscObject)pc,(PetscObject)bjac->y);
820: } else {
821: ksp = jac->ksp[0];
822: bjac = (PC_BJacobi_Singleblock*)jac->data;
823: }
824: if (pc->useAmat) {
825: KSPSetOperators(ksp,mat,pmat);
826: } else {
827: KSPSetOperators(ksp,pmat,pmat);
828: }
829: if (!wasSetup && pc->setfromoptionscalled) {
830: /* If PCSetFromOptions_BJacobi is called later, KSPSetFromOptions will be called at that time. */
831: KSPSetFromOptions(ksp);
832: }
833: return(0);
834: }
836: /* ---------------------------------------------------------------------------------------------*/
837: static PetscErrorCode PCReset_BJacobi_Multiblock(PC pc)
838: {
839: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
840: PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock*)jac->data;
841: PetscErrorCode ierr;
842: PetscInt i;
845: if (bjac && bjac->pmat) {
846: MatDestroyMatrices(jac->n_local,&bjac->pmat);
847: if (pc->useAmat) {
848: MatDestroyMatrices(jac->n_local,&bjac->mat);
849: }
850: }
852: for (i=0; i<jac->n_local; i++) {
853: KSPReset(jac->ksp[i]);
854: if (bjac && bjac->x) {
855: VecDestroy(&bjac->x[i]);
856: VecDestroy(&bjac->y[i]);
857: ISDestroy(&bjac->is[i]);
858: }
859: }
860: PetscFree(jac->l_lens);
861: PetscFree(jac->g_lens);
862: return(0);
863: }
865: static PetscErrorCode PCDestroy_BJacobi_Multiblock(PC pc)
866: {
867: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
868: PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock*)jac->data;
869: PetscErrorCode ierr;
870: PetscInt i;
873: PCReset_BJacobi_Multiblock(pc);
874: if (bjac) {
875: PetscFree2(bjac->x,bjac->y);
876: PetscFree(bjac->starts);
877: PetscFree(bjac->is);
878: }
879: PetscFree(jac->data);
880: for (i=0; i<jac->n_local; i++) {
881: KSPDestroy(&jac->ksp[i]);
882: }
883: PetscFree(jac->ksp);
884: PetscFree(pc->data);
885: return(0);
886: }
888: static PetscErrorCode PCSetUpOnBlocks_BJacobi_Multiblock(PC pc)
889: {
890: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
891: PetscErrorCode ierr;
892: PetscInt i,n_local = jac->n_local;
893: KSPConvergedReason reason;
896: for (i=0; i<n_local; i++) {
897: KSPSetUp(jac->ksp[i]);
898: KSPGetConvergedReason(jac->ksp[i],&reason);
899: if (reason == KSP_DIVERGED_PC_FAILED) {
900: pc->failedreason = PC_SUBPC_ERROR;
901: }
902: }
903: return(0);
904: }
906: /*
907: Preconditioner for block Jacobi
908: */
909: static PetscErrorCode PCApply_BJacobi_Multiblock(PC pc,Vec x,Vec y)
910: {
911: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
912: PetscErrorCode ierr;
913: PetscInt i,n_local = jac->n_local;
914: PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock*)jac->data;
915: PetscScalar *yin;
916: const PetscScalar *xin;
919: VecGetArrayRead(x,&xin);
920: VecGetArray(y,&yin);
921: for (i=0; i<n_local; i++) {
922: /*
923: To avoid copying the subvector from x into a workspace we instead
924: make the workspace vector array point to the subpart of the array of
925: the global vector.
926: */
927: VecPlaceArray(bjac->x[i],xin+bjac->starts[i]);
928: VecPlaceArray(bjac->y[i],yin+bjac->starts[i]);
930: PetscLogEventBegin(PC_ApplyOnBlocks,jac->ksp[i],bjac->x[i],bjac->y[i],0);
931: KSPSolve(jac->ksp[i],bjac->x[i],bjac->y[i]);
932: KSPCheckSolve(jac->ksp[i],pc,bjac->y[i]);
933: PetscLogEventEnd(PC_ApplyOnBlocks,jac->ksp[i],bjac->x[i],bjac->y[i],0);
935: VecResetArray(bjac->x[i]);
936: VecResetArray(bjac->y[i]);
937: }
938: VecRestoreArrayRead(x,&xin);
939: VecRestoreArray(y,&yin);
940: return(0);
941: }
943: /*
944: Preconditioner for block Jacobi
945: */
946: static PetscErrorCode PCApplyTranspose_BJacobi_Multiblock(PC pc,Vec x,Vec y)
947: {
948: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
949: PetscErrorCode ierr;
950: PetscInt i,n_local = jac->n_local;
951: PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock*)jac->data;
952: PetscScalar *yin;
953: const PetscScalar *xin;
956: VecGetArrayRead(x,&xin);
957: VecGetArray(y,&yin);
958: for (i=0; i<n_local; i++) {
959: /*
960: To avoid copying the subvector from x into a workspace we instead
961: make the workspace vector array point to the subpart of the array of
962: the global vector.
963: */
964: VecPlaceArray(bjac->x[i],xin+bjac->starts[i]);
965: VecPlaceArray(bjac->y[i],yin+bjac->starts[i]);
967: PetscLogEventBegin(PC_ApplyTransposeOnBlocks,jac->ksp[i],bjac->x[i],bjac->y[i],0);
968: KSPSolveTranspose(jac->ksp[i],bjac->x[i],bjac->y[i]);
969: KSPCheckSolve(jac->ksp[i],pc,bjac->y[i]);
970: PetscLogEventEnd(PC_ApplyTransposeOnBlocks,jac->ksp[i],bjac->x[i],bjac->y[i],0);
972: VecResetArray(bjac->x[i]);
973: VecResetArray(bjac->y[i]);
974: }
975: VecRestoreArrayRead(x,&xin);
976: VecRestoreArray(y,&yin);
977: return(0);
978: }
980: static PetscErrorCode PCSetUp_BJacobi_Multiblock(PC pc,Mat mat,Mat pmat)
981: {
982: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
983: PetscErrorCode ierr;
984: PetscInt m,n_local,N,M,start,i;
985: const char *prefix,*pprefix,*mprefix;
986: KSP ksp;
987: Vec x,y;
988: PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock*)jac->data;
989: PC subpc;
990: IS is;
991: MatReuse scall;
992: VecType vectype;
995: MatGetLocalSize(pc->pmat,&M,&N);
997: n_local = jac->n_local;
999: if (pc->useAmat) {
1000: PetscBool same;
1001: PetscObjectTypeCompare((PetscObject)mat,((PetscObject)pmat)->type_name,&same);
1002: if (!same) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"Matrices not of same type");
1003: }
1005: if (!pc->setupcalled) {
1006: scall = MAT_INITIAL_MATRIX;
1008: if (!jac->ksp) {
1009: pc->ops->reset = PCReset_BJacobi_Multiblock;
1010: pc->ops->destroy = PCDestroy_BJacobi_Multiblock;
1011: pc->ops->apply = PCApply_BJacobi_Multiblock;
1012: pc->ops->matapply = NULL;
1013: pc->ops->applytranspose= PCApplyTranspose_BJacobi_Multiblock;
1014: pc->ops->setuponblocks = PCSetUpOnBlocks_BJacobi_Multiblock;
1016: PetscNewLog(pc,&bjac);
1017: PetscMalloc1(n_local,&jac->ksp);
1018: PetscLogObjectMemory((PetscObject)pc,sizeof(n_local*sizeof(KSP)));
1019: PetscMalloc2(n_local,&bjac->x,n_local,&bjac->y);
1020: PetscMalloc1(n_local,&bjac->starts);
1021: PetscLogObjectMemory((PetscObject)pc,sizeof(n_local*sizeof(PetscScalar)));
1023: jac->data = (void*)bjac;
1024: PetscMalloc1(n_local,&bjac->is);
1025: PetscLogObjectMemory((PetscObject)pc,sizeof(n_local*sizeof(IS)));
1027: for (i=0; i<n_local; i++) {
1028: KSPCreate(PETSC_COMM_SELF,&ksp);
1029: KSPSetErrorIfNotConverged(ksp,pc->erroriffailure);
1030: PetscObjectIncrementTabLevel((PetscObject)ksp,(PetscObject)pc,1);
1031: PetscLogObjectParent((PetscObject)pc,(PetscObject)ksp);
1032: KSPSetType(ksp,KSPPREONLY);
1033: KSPGetPC(ksp,&subpc);
1034: PCGetOptionsPrefix(pc,&prefix);
1035: KSPSetOptionsPrefix(ksp,prefix);
1036: KSPAppendOptionsPrefix(ksp,"sub_");
1038: jac->ksp[i] = ksp;
1039: }
1040: } else {
1041: bjac = (PC_BJacobi_Multiblock*)jac->data;
1042: }
1044: start = 0;
1045: MatGetVecType(pmat,&vectype);
1046: for (i=0; i<n_local; i++) {
1047: m = jac->l_lens[i];
1048: /*
1049: The reason we need to generate these vectors is to serve
1050: as the right-hand side and solution vector for the solve on the
1051: block. We do not need to allocate space for the vectors since
1052: that is provided via VecPlaceArray() just before the call to
1053: KSPSolve() on the block.
1055: */
1056: VecCreateSeq(PETSC_COMM_SELF,m,&x);
1057: VecCreateSeqWithArray(PETSC_COMM_SELF,1,m,NULL,&y);
1058: VecSetType(x,vectype);
1059: VecSetType(y,vectype);
1060: PetscLogObjectParent((PetscObject)pc,(PetscObject)x);
1061: PetscLogObjectParent((PetscObject)pc,(PetscObject)y);
1063: bjac->x[i] = x;
1064: bjac->y[i] = y;
1065: bjac->starts[i] = start;
1067: ISCreateStride(PETSC_COMM_SELF,m,start,1,&is);
1068: bjac->is[i] = is;
1069: PetscLogObjectParent((PetscObject)pc,(PetscObject)is);
1071: start += m;
1072: }
1073: } else {
1074: bjac = (PC_BJacobi_Multiblock*)jac->data;
1075: /*
1076: Destroy the blocks from the previous iteration
1077: */
1078: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
1079: MatDestroyMatrices(n_local,&bjac->pmat);
1080: if (pc->useAmat) {
1081: MatDestroyMatrices(n_local,&bjac->mat);
1082: }
1083: scall = MAT_INITIAL_MATRIX;
1084: } else scall = MAT_REUSE_MATRIX;
1085: }
1087: MatCreateSubMatrices(pmat,n_local,bjac->is,bjac->is,scall,&bjac->pmat);
1088: if (pc->useAmat) {
1089: MatCreateSubMatrices(mat,n_local,bjac->is,bjac->is,scall,&bjac->mat);
1090: }
1091: /* Return control to the user so that the submatrices can be modified (e.g., to apply
1092: different boundary conditions for the submatrices than for the global problem) */
1093: PCModifySubMatrices(pc,n_local,bjac->is,bjac->is,bjac->pmat,pc->modifysubmatricesP);
1095: PetscObjectGetOptionsPrefix((PetscObject)pmat,&pprefix);
1096: for (i=0; i<n_local; i++) {
1097: PetscLogObjectParent((PetscObject)pc,(PetscObject)bjac->pmat[i]);
1098: PetscObjectSetOptionsPrefix((PetscObject)bjac->pmat[i],pprefix);
1099: if (pc->useAmat) {
1100: PetscLogObjectParent((PetscObject)pc,(PetscObject)bjac->mat[i]);
1101: PetscObjectGetOptionsPrefix((PetscObject)mat,&mprefix);
1102: PetscObjectSetOptionsPrefix((PetscObject)bjac->mat[i],mprefix);
1103: KSPSetOperators(jac->ksp[i],bjac->mat[i],bjac->pmat[i]);
1104: } else {
1105: KSPSetOperators(jac->ksp[i],bjac->pmat[i],bjac->pmat[i]);
1106: }
1107: if (pc->setfromoptionscalled) {
1108: KSPSetFromOptions(jac->ksp[i]);
1109: }
1110: }
1111: return(0);
1112: }
1114: /* ---------------------------------------------------------------------------------------------*/
1115: /*
1116: These are for a single block with multiple processes;
1117: */
1118: static PetscErrorCode PCReset_BJacobi_Multiproc(PC pc)
1119: {
1120: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
1121: PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc*)jac->data;
1122: PetscErrorCode ierr;
1125: VecDestroy(&mpjac->ysub);
1126: VecDestroy(&mpjac->xsub);
1127: MatDestroy(&mpjac->submats);
1128: if (jac->ksp) {KSPReset(jac->ksp[0]);}
1129: return(0);
1130: }
1132: static PetscErrorCode PCDestroy_BJacobi_Multiproc(PC pc)
1133: {
1134: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
1135: PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc*)jac->data;
1136: PetscErrorCode ierr;
1139: PCReset_BJacobi_Multiproc(pc);
1140: KSPDestroy(&jac->ksp[0]);
1141: PetscFree(jac->ksp);
1142: PetscSubcommDestroy(&mpjac->psubcomm);
1144: PetscFree(mpjac);
1145: PetscFree(pc->data);
1146: return(0);
1147: }
1149: static PetscErrorCode PCApply_BJacobi_Multiproc(PC pc,Vec x,Vec y)
1150: {
1151: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
1152: PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc*)jac->data;
1153: PetscErrorCode ierr;
1154: PetscScalar *yarray;
1155: const PetscScalar *xarray;
1156: KSPConvergedReason reason;
1159: /* place x's and y's local arrays into xsub and ysub */
1160: VecGetArrayRead(x,&xarray);
1161: VecGetArray(y,&yarray);
1162: VecPlaceArray(mpjac->xsub,xarray);
1163: VecPlaceArray(mpjac->ysub,yarray);
1165: /* apply preconditioner on each matrix block */
1166: PetscLogEventBegin(PC_ApplyOnBlocks,jac->ksp[0],mpjac->xsub,mpjac->ysub,0);
1167: KSPSolve(jac->ksp[0],mpjac->xsub,mpjac->ysub);
1168: KSPCheckSolve(jac->ksp[0],pc,mpjac->ysub);
1169: PetscLogEventEnd(PC_ApplyOnBlocks,jac->ksp[0],mpjac->xsub,mpjac->ysub,0);
1170: KSPGetConvergedReason(jac->ksp[0],&reason);
1171: if (reason == KSP_DIVERGED_PC_FAILED) {
1172: pc->failedreason = PC_SUBPC_ERROR;
1173: }
1175: VecResetArray(mpjac->xsub);
1176: VecResetArray(mpjac->ysub);
1177: VecRestoreArrayRead(x,&xarray);
1178: VecRestoreArray(y,&yarray);
1179: return(0);
1180: }
1182: static PetscErrorCode PCMatApply_BJacobi_Multiproc(PC pc,Mat X,Mat Y)
1183: {
1184: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
1185: KSPConvergedReason reason;
1186: Mat sX,sY;
1187: const PetscScalar *x;
1188: PetscScalar *y;
1189: PetscInt m,N,lda,ldb;
1190: PetscErrorCode ierr;
1193: /* apply preconditioner on each matrix block */
1194: MatGetLocalSize(X,&m,NULL);
1195: MatGetSize(X,NULL,&N);
1196: MatDenseGetLDA(X,&lda);
1197: MatDenseGetLDA(Y,&ldb);
1198: MatDenseGetArrayRead(X,&x);
1199: MatDenseGetArrayWrite(Y,&y);
1200: MatCreateDense(PetscObjectComm((PetscObject)jac->ksp[0]),m,PETSC_DECIDE,PETSC_DECIDE,N,(PetscScalar*)x,&sX);
1201: MatCreateDense(PetscObjectComm((PetscObject)jac->ksp[0]),m,PETSC_DECIDE,PETSC_DECIDE,N,y,&sY);
1202: MatDenseSetLDA(sX,lda);
1203: MatDenseSetLDA(sY,ldb);
1204: PetscLogEventBegin(PC_ApplyOnBlocks,jac->ksp[0],X,Y,0);
1205: KSPMatSolve(jac->ksp[0],sX,sY);
1206: KSPCheckSolve(jac->ksp[0],pc,NULL);
1207: PetscLogEventEnd(PC_ApplyOnBlocks,jac->ksp[0],X,Y,0);
1208: MatDestroy(&sY);
1209: MatDestroy(&sX);
1210: MatDenseRestoreArrayWrite(Y,&y);
1211: MatDenseRestoreArrayRead(X,&x);
1212: KSPGetConvergedReason(jac->ksp[0],&reason);
1213: if (reason == KSP_DIVERGED_PC_FAILED) {
1214: pc->failedreason = PC_SUBPC_ERROR;
1215: }
1216: return(0);
1217: }
1219: static PetscErrorCode PCSetUp_BJacobi_Multiproc(PC pc)
1220: {
1221: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
1222: PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc*)jac->data;
1223: PetscErrorCode ierr;
1224: PetscInt m,n;
1225: MPI_Comm comm,subcomm=0;
1226: const char *prefix;
1227: PetscBool wasSetup = PETSC_TRUE;
1228: VecType vectype;
1232: PetscObjectGetComm((PetscObject)pc,&comm);
1233: if (jac->n_local > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Only a single block in a subcommunicator is supported");
1234: jac->n_local = 1; /* currently only a single block is supported for a subcommunicator */
1235: if (!pc->setupcalled) {
1236: wasSetup = PETSC_FALSE;
1237: PetscNewLog(pc,&mpjac);
1238: jac->data = (void*)mpjac;
1240: /* initialize datastructure mpjac */
1241: if (!jac->psubcomm) {
1242: /* Create default contiguous subcommunicatiors if user does not provide them */
1243: PetscSubcommCreate(comm,&jac->psubcomm);
1244: PetscSubcommSetNumber(jac->psubcomm,jac->n);
1245: PetscSubcommSetType(jac->psubcomm,PETSC_SUBCOMM_CONTIGUOUS);
1246: PetscLogObjectMemory((PetscObject)pc,sizeof(PetscSubcomm));
1247: }
1248: mpjac->psubcomm = jac->psubcomm;
1249: subcomm = PetscSubcommChild(mpjac->psubcomm);
1251: /* Get matrix blocks of pmat */
1252: MatGetMultiProcBlock(pc->pmat,subcomm,MAT_INITIAL_MATRIX,&mpjac->submats);
1254: /* create a new PC that processors in each subcomm have copy of */
1255: PetscMalloc1(1,&jac->ksp);
1256: KSPCreate(subcomm,&jac->ksp[0]);
1257: KSPSetErrorIfNotConverged(jac->ksp[0],pc->erroriffailure);
1258: PetscObjectIncrementTabLevel((PetscObject)jac->ksp[0],(PetscObject)pc,1);
1259: PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->ksp[0]);
1260: KSPSetOperators(jac->ksp[0],mpjac->submats,mpjac->submats);
1261: KSPGetPC(jac->ksp[0],&mpjac->pc);
1263: PCGetOptionsPrefix(pc,&prefix);
1264: KSPSetOptionsPrefix(jac->ksp[0],prefix);
1265: KSPAppendOptionsPrefix(jac->ksp[0],"sub_");
1267: /* create dummy vectors xsub and ysub */
1268: MatGetLocalSize(mpjac->submats,&m,&n);
1269: VecCreateMPIWithArray(subcomm,1,n,PETSC_DECIDE,NULL,&mpjac->xsub);
1270: VecCreateMPIWithArray(subcomm,1,m,PETSC_DECIDE,NULL,&mpjac->ysub);
1271: MatGetVecType(mpjac->submats,&vectype);
1272: VecSetType(mpjac->xsub,vectype);
1273: VecSetType(mpjac->ysub,vectype);
1274: PetscLogObjectParent((PetscObject)pc,(PetscObject)mpjac->xsub);
1275: PetscLogObjectParent((PetscObject)pc,(PetscObject)mpjac->ysub);
1277: pc->ops->reset = PCReset_BJacobi_Multiproc;
1278: pc->ops->destroy = PCDestroy_BJacobi_Multiproc;
1279: pc->ops->apply = PCApply_BJacobi_Multiproc;
1280: pc->ops->matapply= PCMatApply_BJacobi_Multiproc;
1281: } else { /* pc->setupcalled */
1282: subcomm = PetscSubcommChild(mpjac->psubcomm);
1283: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
1284: /* destroy old matrix blocks, then get new matrix blocks */
1285: if (mpjac->submats) {MatDestroy(&mpjac->submats);}
1286: MatGetMultiProcBlock(pc->pmat,subcomm,MAT_INITIAL_MATRIX,&mpjac->submats);
1287: } else {
1288: MatGetMultiProcBlock(pc->pmat,subcomm,MAT_REUSE_MATRIX,&mpjac->submats);
1289: }
1290: KSPSetOperators(jac->ksp[0],mpjac->submats,mpjac->submats);
1291: }
1293: if (!wasSetup && pc->setfromoptionscalled) {
1294: KSPSetFromOptions(jac->ksp[0]);
1295: }
1296: return(0);
1297: }