Actual source code: bjacobi.c
petsc-3.7.7 2017-09-25
2: /*
3: Defines a block Jacobi preconditioner.
4: */
5: #include <petsc/private/pcimpl.h> /*I "petscpc.h" I*/
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);
14: static PetscErrorCode PCSetUp_BJacobi(PC pc)
15: {
16: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
17: Mat mat = pc->mat,pmat = pc->pmat;
18: PetscErrorCode ierr,(*f)(Mat,Mat*);
19: PetscInt N,M,start,i,sum,end;
20: PetscInt bs,i_start=-1,i_end=-1;
21: PetscMPIInt rank,size;
22: const char *pprefix,*mprefix;
25: MPI_Comm_rank(PetscObjectComm((PetscObject)pc),&rank);
26: MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);
27: MatGetLocalSize(pc->pmat,&M,&N);
28: MatGetBlockSize(pc->pmat,&bs);
30: if (jac->n > 0 && jac->n < size) {
31: PCSetUp_BJacobi_Multiproc(pc);
32: return(0);
33: }
35: /* --------------------------------------------------------------------------
36: Determines the number of blocks assigned to each processor
37: -----------------------------------------------------------------------------*/
39: /* local block count given */
40: if (jac->n_local > 0 && jac->n < 0) {
41: MPIU_Allreduce(&jac->n_local,&jac->n,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)pc));
42: if (jac->l_lens) { /* check that user set these correctly */
43: sum = 0;
44: for (i=0; i<jac->n_local; i++) {
45: 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");
46: sum += jac->l_lens[i];
47: }
48: if (sum != M) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local lens set incorrectly");
49: } else {
50: PetscMalloc1(jac->n_local,&jac->l_lens);
51: for (i=0; i<jac->n_local; i++) jac->l_lens[i] = bs*((M/bs)/jac->n_local + (((M/bs) % jac->n_local) > i));
52: }
53: } else if (jac->n > 0 && jac->n_local < 0) { /* global block count given */
54: /* global blocks given: determine which ones are local */
55: if (jac->g_lens) {
56: /* check if the g_lens is has valid entries */
57: for (i=0; i<jac->n; i++) {
58: if (!jac->g_lens[i]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Zero block not allowed");
59: 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");
60: }
61: if (size == 1) {
62: jac->n_local = jac->n;
63: PetscMalloc1(jac->n_local,&jac->l_lens);
64: PetscMemcpy(jac->l_lens,jac->g_lens,jac->n_local*sizeof(PetscInt));
65: /* check that user set these correctly */
66: sum = 0;
67: for (i=0; i<jac->n_local; i++) sum += jac->l_lens[i];
68: if (sum != M) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Global lens set incorrectly");
69: } else {
70: MatGetOwnershipRange(pc->pmat,&start,&end);
71: /* loop over blocks determing first one owned by me */
72: sum = 0;
73: for (i=0; i<jac->n+1; i++) {
74: if (sum == start) { i_start = i; goto start_1;}
75: if (i < jac->n) sum += jac->g_lens[i];
76: }
77: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Block sizes used in PCBJacobiSetTotalBlocks()\nare not compatible with parallel matrix layout");
78: start_1:
79: for (i=i_start; i<jac->n+1; i++) {
80: if (sum == end) { i_end = i; goto end_1; }
81: if (i < jac->n) sum += jac->g_lens[i];
82: }
83: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Block sizes used in PCBJacobiSetTotalBlocks()\nare not compatible with parallel matrix layout");
84: end_1:
85: jac->n_local = i_end - i_start;
86: PetscMalloc1(jac->n_local,&jac->l_lens);
87: PetscMemcpy(jac->l_lens,jac->g_lens+i_start,jac->n_local*sizeof(PetscInt));
88: }
89: } else { /* no global blocks given, determine then using default layout */
90: jac->n_local = jac->n/size + ((jac->n % size) > rank);
91: PetscMalloc1(jac->n_local,&jac->l_lens);
92: for (i=0; i<jac->n_local; i++) {
93: jac->l_lens[i] = ((M/bs)/jac->n_local + (((M/bs) % jac->n_local) > i))*bs;
94: if (!jac->l_lens[i]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Too many blocks given");
95: }
96: }
97: } else if (jac->n < 0 && jac->n_local < 0) { /* no blocks given */
98: jac->n = size;
99: jac->n_local = 1;
100: PetscMalloc1(1,&jac->l_lens);
101: jac->l_lens[0] = M;
102: } else { /* jac->n > 0 && jac->n_local > 0 */
103: if (!jac->l_lens) {
104: PetscMalloc1(jac->n_local,&jac->l_lens);
105: for (i=0; i<jac->n_local; i++) jac->l_lens[i] = bs*((M/bs)/jac->n_local + (((M/bs) % jac->n_local) > i));
106: }
107: }
108: if (jac->n_local < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Number of blocks is less than number of processors");
110: /* -------------------------
111: Determines mat and pmat
112: ---------------------------*/
113: PetscObjectQueryFunction((PetscObject)pc->mat,"MatGetDiagonalBlock_C",&f);
114: if (!f && size == 1) {
115: mat = pc->mat;
116: pmat = pc->pmat;
117: } else {
118: if (pc->useAmat) {
119: /* use block from Amat matrix, not Pmat for local MatMult() */
120: MatGetDiagonalBlock(pc->mat,&mat);
121: /* make submatrix have same prefix as entire matrix */
122: PetscObjectGetOptionsPrefix((PetscObject)pc->mat,&mprefix);
123: PetscObjectSetOptionsPrefix((PetscObject)mat,mprefix);
124: }
125: if (pc->pmat != pc->mat || !pc->useAmat) {
126: MatGetDiagonalBlock(pc->pmat,&pmat);
127: /* make submatrix have same prefix as entire matrix */
128: PetscObjectGetOptionsPrefix((PetscObject)pc->pmat,&pprefix);
129: PetscObjectSetOptionsPrefix((PetscObject)pmat,pprefix);
130: } else pmat = mat;
131: }
133: /* ------
134: Setup code depends on the number of blocks
135: */
136: if (jac->n_local == 1) {
137: PCSetUp_BJacobi_Singleblock(pc,mat,pmat);
138: } else {
139: PCSetUp_BJacobi_Multiblock(pc,mat,pmat);
140: }
141: return(0);
142: }
144: /* Default destroy, if it has never been setup */
147: static PetscErrorCode PCDestroy_BJacobi(PC pc)
148: {
149: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
153: PetscFree(jac->g_lens);
154: PetscFree(jac->l_lens);
155: PetscFree(pc->data);
156: return(0);
157: }
162: static PetscErrorCode PCSetFromOptions_BJacobi(PetscOptionItems *PetscOptionsObject,PC pc)
163: {
164: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
166: PetscInt blocks,i;
167: PetscBool flg;
170: PetscOptionsHead(PetscOptionsObject,"Block Jacobi options");
171: PetscOptionsInt("-pc_bjacobi_blocks","Total number of blocks","PCBJacobiSetTotalBlocks",jac->n,&blocks,&flg);
172: if (flg) {
173: PCBJacobiSetTotalBlocks(pc,blocks,NULL);
174: }
175: if (jac->ksp) {
176: /* The sub-KSP has already been set up (e.g., PCSetUp_BJacobi_Singleblock), but KSPSetFromOptions was not called
177: * unless we had already been called. */
178: for (i=0; i<jac->n_local; i++) {
179: KSPSetFromOptions(jac->ksp[i]);
180: }
181: }
182: PetscOptionsTail();
183: return(0);
184: }
186: #include <petscdraw.h>
189: static PetscErrorCode PCView_BJacobi(PC pc,PetscViewer viewer)
190: {
191: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
192: PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc*)jac->data;
193: PetscErrorCode ierr;
194: PetscMPIInt rank;
195: PetscInt i;
196: PetscBool iascii,isstring,isdraw;
197: PetscViewer sviewer;
200: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
201: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
202: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
203: if (iascii) {
204: if (pc->useAmat) {
205: PetscViewerASCIIPrintf(viewer," block Jacobi: using Amat local matrix, number of blocks = %D\n",jac->n);
206: }
207: PetscViewerASCIIPrintf(viewer," block Jacobi: number of blocks = %D\n",jac->n);
208: MPI_Comm_rank(PetscObjectComm((PetscObject)pc),&rank);
209: if (jac->same_local_solves) {
210: PetscViewerASCIIPrintf(viewer," Local solve is same for all blocks, in the following KSP and PC objects:\n");
211: if (jac->ksp && !jac->psubcomm) {
212: PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
213: if (!rank) {
214: PetscViewerASCIIPushTab(viewer);
215: KSPView(jac->ksp[0],sviewer);
216: PetscViewerASCIIPopTab(viewer);
217: }
218: PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
219: } else if (jac->psubcomm && !jac->psubcomm->color) {
220: PetscViewerASCIIGetStdout(PetscSubcommChild(mpjac->psubcomm),&sviewer);
221: PetscViewerASCIIPushTab(viewer);
222: KSPView(*(jac->ksp),sviewer);
223: PetscViewerASCIIPopTab(viewer);
224: }
225: } else {
226: PetscInt n_global;
227: MPIU_Allreduce(&jac->n_local,&n_global,1,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)pc));
228: PetscViewerASCIIPushSynchronized(viewer);
229: PetscViewerASCIIPrintf(viewer," Local solve info for each block is in the following KSP and PC objects:\n");
230: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] number of local blocks = %D, first local block number = %D\n",
231: rank,jac->n_local,jac->first_local);
232: PetscViewerASCIIPushTab(viewer);
233: PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
234: for (i=0; i<jac->n_local; i++) {
235: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] local block number %D\n",rank,i);
236: KSPView(jac->ksp[i],sviewer);
237: PetscViewerASCIISynchronizedPrintf(viewer,"- - - - - - - - - - - - - - - - - -\n");
238: }
239: PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
240: PetscViewerASCIIPopTab(viewer);
241: PetscViewerFlush(viewer);
242: PetscViewerASCIIPopSynchronized(viewer);
243: }
244: } else if (isstring) {
245: PetscViewerStringSPrintf(viewer," blks=%D",jac->n);
246: PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
247: if (jac->ksp) {KSPView(jac->ksp[0],sviewer);}
248: PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
249: } else if (isdraw) {
250: PetscDraw draw;
251: char str[25];
252: PetscReal x,y,bottom,h;
254: PetscViewerDrawGetDraw(viewer,0,&draw);
255: PetscDrawGetCurrentPoint(draw,&x,&y);
256: PetscSNPrintf(str,25,"Number blocks %D",jac->n);
257: PetscDrawStringBoxed(draw,x,y,PETSC_DRAW_RED,PETSC_DRAW_BLACK,str,NULL,&h);
258: bottom = y - h;
259: PetscDrawPushCurrentPoint(draw,x,bottom);
260: /* warning the communicator on viewer is different then on ksp in parallel */
261: if (jac->ksp) {KSPView(jac->ksp[0],viewer);}
262: PetscDrawPopCurrentPoint(draw);
263: }
264: return(0);
265: }
267: /* -------------------------------------------------------------------------------------*/
271: static PetscErrorCode PCBJacobiGetSubKSP_BJacobi(PC pc,PetscInt *n_local,PetscInt *first_local,KSP **ksp)
272: {
273: PC_BJacobi *jac = (PC_BJacobi*)pc->data;;
276: if (!pc->setupcalled) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must call KSPSetUp() or PCSetUp() first");
278: if (n_local) *n_local = jac->n_local;
279: if (first_local) *first_local = jac->first_local;
280: *ksp = jac->ksp;
281: jac->same_local_solves = PETSC_FALSE; /* Assume that local solves are now different;
282: not necessarily true though! This flag is
283: used only for PCView_BJacobi() */
284: return(0);
285: }
289: static PetscErrorCode PCBJacobiSetTotalBlocks_BJacobi(PC pc,PetscInt blocks,PetscInt *lens)
290: {
291: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
295: 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");
296: jac->n = blocks;
297: if (!lens) jac->g_lens = 0;
298: else {
299: PetscMalloc1(blocks,&jac->g_lens);
300: PetscLogObjectMemory((PetscObject)pc,blocks*sizeof(PetscInt));
301: PetscMemcpy(jac->g_lens,lens,blocks*sizeof(PetscInt));
302: }
303: return(0);
304: }
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: }
320: static PetscErrorCode PCBJacobiSetLocalBlocks_BJacobi(PC pc,PetscInt blocks,const PetscInt lens[])
321: {
322: PC_BJacobi *jac;
326: jac = (PC_BJacobi*)pc->data;
328: jac->n_local = blocks;
329: if (!lens) jac->l_lens = 0;
330: else {
331: PetscMalloc1(blocks,&jac->l_lens);
332: PetscLogObjectMemory((PetscObject)pc,blocks*sizeof(PetscInt));
333: PetscMemcpy(jac->l_lens,lens,blocks*sizeof(PetscInt));
334: }
335: return(0);
336: }
340: static PetscErrorCode PCBJacobiGetLocalBlocks_BJacobi(PC pc, PetscInt *blocks, const PetscInt *lens[])
341: {
342: PC_BJacobi *jac = (PC_BJacobi*) pc->data;
345: *blocks = jac->n_local;
346: if (lens) *lens = jac->l_lens;
347: return(0);
348: }
350: /* -------------------------------------------------------------------------------------*/
354: /*@C
355: PCBJacobiGetSubKSP - Gets the local KSP contexts for all blocks on
356: this processor.
358: Note Collective
360: Input Parameter:
361: . pc - the preconditioner context
363: Output Parameters:
364: + n_local - the number of blocks on this processor, or NULL
365: . first_local - the global number of the first block on this processor, or NULL
366: - ksp - the array of KSP contexts
368: Notes:
369: After PCBJacobiGetSubKSP() the array of KSP contexts is not to be freed.
371: Currently for some matrix implementations only 1 block per processor
372: is supported.
374: You must call KSPSetUp() or PCSetUp() before calling PCBJacobiGetSubKSP().
376: Fortran Usage: You must pass in a KSP array that is large enough to contain all the local KSPs.
377: You can call PCBJacobiGetSubKSP(pc,nlocal,firstlocal,NULL_OBJECT,ierr) to determine how large the
378: KSP array must be.
380: Level: advanced
382: .keywords: block, Jacobi, get, sub, KSP, context
384: .seealso: PCBJacobiGetSubKSP()
385: @*/
386: PetscErrorCode PCBJacobiGetSubKSP(PC pc,PetscInt *n_local,PetscInt *first_local,KSP *ksp[])
387: {
392: PetscUseMethod(pc,"PCBJacobiGetSubKSP_C",(PC,PetscInt*,PetscInt*,KSP **),(pc,n_local,first_local,ksp));
393: return(0);
394: }
398: /*@
399: PCBJacobiSetTotalBlocks - Sets the global number of blocks for the block
400: Jacobi preconditioner.
402: Collective on PC
404: Input Parameters:
405: + pc - the preconditioner context
406: . blocks - the number of blocks
407: - lens - [optional] integer array containing the size of each block
409: Options Database Key:
410: . -pc_bjacobi_blocks <blocks> - Sets the number of global blocks
412: Notes:
413: Currently only a limited number of blocking configurations are supported.
414: All processors sharing the PC must call this routine with the same data.
416: Level: intermediate
418: .keywords: set, number, Jacobi, global, total, blocks
420: .seealso: PCSetUseAmat(), PCBJacobiSetLocalBlocks()
421: @*/
422: PetscErrorCode PCBJacobiSetTotalBlocks(PC pc,PetscInt blocks,const PetscInt lens[])
423: {
428: if (blocks <= 0) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Must have positive blocks");
429: PetscTryMethod(pc,"PCBJacobiSetTotalBlocks_C",(PC,PetscInt,const PetscInt[]),(pc,blocks,lens));
430: return(0);
431: }
435: /*@C
436: PCBJacobiGetTotalBlocks - Gets the global number of blocks for the block
437: Jacobi preconditioner.
439: Not Collective
441: Input Parameter:
442: . pc - the preconditioner context
444: Output parameters:
445: + blocks - the number of blocks
446: - lens - integer array containing the size of each block
448: Level: intermediate
450: .keywords: get, number, Jacobi, global, total, blocks
452: .seealso: PCSetUseAmat(), PCBJacobiGetLocalBlocks()
453: @*/
454: PetscErrorCode PCBJacobiGetTotalBlocks(PC pc, PetscInt *blocks, const PetscInt *lens[])
455: {
461: PetscUseMethod(pc,"PCBJacobiGetTotalBlocks_C",(PC,PetscInt*, const PetscInt *[]),(pc,blocks,lens));
462: return(0);
463: }
467: /*@
468: PCBJacobiSetLocalBlocks - Sets the local number of blocks for the block
469: Jacobi preconditioner.
471: Not Collective
473: Input Parameters:
474: + pc - the preconditioner context
475: . blocks - the number of blocks
476: - lens - [optional] integer array containing size of each block
478: Note:
479: Currently only a limited number of blocking configurations are supported.
481: Level: intermediate
483: .keywords: PC, set, number, Jacobi, local, blocks
485: .seealso: PCSetUseAmat(), PCBJacobiSetTotalBlocks()
486: @*/
487: PetscErrorCode PCBJacobiSetLocalBlocks(PC pc,PetscInt blocks,const PetscInt lens[])
488: {
493: if (blocks < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must have nonegative blocks");
494: PetscTryMethod(pc,"PCBJacobiSetLocalBlocks_C",(PC,PetscInt,const PetscInt []),(pc,blocks,lens));
495: return(0);
496: }
500: /*@C
501: PCBJacobiGetLocalBlocks - Gets the local number of blocks for the block
502: Jacobi preconditioner.
504: Not Collective
506: Input Parameters:
507: + pc - the preconditioner context
508: . blocks - the number of blocks
509: - lens - [optional] integer array containing size of each block
511: Note:
512: Currently only a limited number of blocking configurations are supported.
514: Level: intermediate
516: .keywords: PC, get, number, Jacobi, local, blocks
518: .seealso: PCSetUseAmat(), PCBJacobiGetTotalBlocks()
519: @*/
520: PetscErrorCode PCBJacobiGetLocalBlocks(PC pc, PetscInt *blocks, const PetscInt *lens[])
521: {
527: PetscUseMethod(pc,"PCBJacobiGetLocalBlocks_C",(PC,PetscInt*, const PetscInt *[]),(pc,blocks,lens));
528: return(0);
529: }
531: /* -----------------------------------------------------------------------------------*/
533: /*MC
534: PCBJACOBI - Use block Jacobi preconditioning, each block is (approximately) solved with
535: its own KSP object.
537: Options Database Keys:
538: . -pc_use_amat - use Amat to apply block of operator in inner Krylov method
540: Notes: Each processor can have one or more blocks, but a block cannot be shared by more
541: than one processor. Defaults to one block per processor.
543: To set options on the solvers for each block append -sub_ to all the KSP, KSP, and PC
544: options database keys. For example, -sub_pc_type ilu -sub_pc_factor_levels 1 -sub_ksp_type preonly
546: To set the options on the solvers separate for each block call PCBJacobiGetSubKSP()
547: and set the options directly on the resulting KSP object (you can access its PC
548: KSPGetPC())
550: For CUSP vectors it is recommended to use exactly one block per MPI process for best
551: performance. Different block partitioning may lead to additional data transfers
552: between host and GPU that lead to degraded performance.
554: Level: beginner
556: Concepts: block Jacobi
559: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
560: PCASM, PCSetUseAmat(), PCGetUseAmat(), PCBJacobiGetSubKSP(), PCBJacobiSetTotalBlocks(),
561: PCBJacobiSetLocalBlocks(), PCSetModifySubmatrices()
562: M*/
566: PETSC_EXTERN PetscErrorCode PCCreate_BJacobi(PC pc)
567: {
569: PetscMPIInt rank;
570: PC_BJacobi *jac;
573: PetscNewLog(pc,&jac);
574: MPI_Comm_rank(PetscObjectComm((PetscObject)pc),&rank);
576: pc->ops->apply = 0;
577: pc->ops->applytranspose = 0;
578: pc->ops->setup = PCSetUp_BJacobi;
579: pc->ops->destroy = PCDestroy_BJacobi;
580: pc->ops->setfromoptions = PCSetFromOptions_BJacobi;
581: pc->ops->view = PCView_BJacobi;
582: pc->ops->applyrichardson = 0;
584: pc->data = (void*)jac;
585: jac->n = -1;
586: jac->n_local = -1;
587: jac->first_local = rank;
588: jac->ksp = 0;
589: jac->same_local_solves = PETSC_TRUE;
590: jac->g_lens = 0;
591: jac->l_lens = 0;
592: jac->psubcomm = 0;
594: PetscObjectComposeFunction((PetscObject)pc,"PCBJacobiGetSubKSP_C",PCBJacobiGetSubKSP_BJacobi);
595: PetscObjectComposeFunction((PetscObject)pc,"PCBJacobiSetTotalBlocks_C",PCBJacobiSetTotalBlocks_BJacobi);
596: PetscObjectComposeFunction((PetscObject)pc,"PCBJacobiGetTotalBlocks_C",PCBJacobiGetTotalBlocks_BJacobi);
597: PetscObjectComposeFunction((PetscObject)pc,"PCBJacobiSetLocalBlocks_C",PCBJacobiSetLocalBlocks_BJacobi);
598: PetscObjectComposeFunction((PetscObject)pc,"PCBJacobiGetLocalBlocks_C",PCBJacobiGetLocalBlocks_BJacobi);
599: return(0);
600: }
602: /* --------------------------------------------------------------------------------------------*/
603: /*
604: These are for a single block per processor; works for AIJ, BAIJ; Seq and MPI
605: */
608: static PetscErrorCode PCReset_BJacobi_Singleblock(PC pc)
609: {
610: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
611: PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock*)jac->data;
612: PetscErrorCode ierr;
615: KSPReset(jac->ksp[0]);
616: VecDestroy(&bjac->x);
617: VecDestroy(&bjac->y);
618: return(0);
619: }
623: static PetscErrorCode PCDestroy_BJacobi_Singleblock(PC pc)
624: {
625: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
626: PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock*)jac->data;
627: PetscErrorCode ierr;
630: PCReset_BJacobi_Singleblock(pc);
631: KSPDestroy(&jac->ksp[0]);
632: PetscFree(jac->ksp);
633: PetscFree(jac->l_lens);
634: PetscFree(jac->g_lens);
635: PetscFree(bjac);
636: PetscFree(pc->data);
637: return(0);
638: }
640: #include <petsc/private/kspimpl.h>
643: static PetscErrorCode PCSetUpOnBlocks_BJacobi_Singleblock(PC pc)
644: {
646: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
647: KSP subksp = jac->ksp[0];
650: KSPSetUp(subksp);
651: if (subksp->reason == KSP_DIVERGED_PCSETUP_FAILED) {
652: pc->failedreason = PC_SUBPC_ERROR;
653: }
654: return(0);
655: }
659: static PetscErrorCode PCApply_BJacobi_Singleblock(PC pc,Vec x,Vec y)
660: {
661: PetscErrorCode ierr;
662: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
663: PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock*)jac->data;
666: VecGetLocalVectorRead(x, bjac->x);
667: VecGetLocalVector(y, bjac->y);
668: /* Since the inner KSP matrix may point directly to the diagonal block of an MPI matrix the inner
669: matrix may change even if the outter KSP/PC has not updated the preconditioner, this will trigger a rebuild
670: of the inner preconditioner automatically unless we pass down the outter preconditioners reuse flag.*/
671: KSPSetReusePreconditioner(jac->ksp[0],pc->reusepreconditioner);
672: KSPSolve(jac->ksp[0],bjac->x,bjac->y);
673: VecRestoreLocalVectorRead(x, bjac->x);
674: VecRestoreLocalVector(y, bjac->y);
675: return(0);
676: }
680: static PetscErrorCode PCApplySymmetricLeft_BJacobi_Singleblock(PC pc,Vec x,Vec y)
681: {
682: PetscErrorCode ierr;
683: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
684: PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock*)jac->data;
685: PetscScalar *y_array;
686: const PetscScalar *x_array;
687: PC subpc;
690: /*
691: The VecPlaceArray() is to avoid having to copy the
692: y vector into the bjac->x vector. The reason for
693: the bjac->x vector is that we need a sequential vector
694: for the sequential solve.
695: */
696: VecGetArrayRead(x,&x_array);
697: VecGetArray(y,&y_array);
698: VecPlaceArray(bjac->x,x_array);
699: VecPlaceArray(bjac->y,y_array);
700: /* apply the symmetric left portion of the inner PC operator */
701: /* note this by-passes the inner KSP and its options completely */
702: KSPGetPC(jac->ksp[0],&subpc);
703: PCApplySymmetricLeft(subpc,bjac->x,bjac->y);
704: VecResetArray(bjac->x);
705: VecResetArray(bjac->y);
706: VecRestoreArrayRead(x,&x_array);
707: VecRestoreArray(y,&y_array);
708: return(0);
709: }
713: static PetscErrorCode PCApplySymmetricRight_BJacobi_Singleblock(PC pc,Vec x,Vec y)
714: {
715: PetscErrorCode ierr;
716: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
717: PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock*)jac->data;
718: PetscScalar *y_array;
719: const PetscScalar *x_array;
720: PC subpc;
723: /*
724: The VecPlaceArray() is to avoid having to copy the
725: y vector into the bjac->x vector. The reason for
726: the bjac->x vector is that we need a sequential vector
727: for the sequential solve.
728: */
729: VecGetArrayRead(x,&x_array);
730: VecGetArray(y,&y_array);
731: VecPlaceArray(bjac->x,x_array);
732: VecPlaceArray(bjac->y,y_array);
734: /* apply the symmetric right portion of the inner PC operator */
735: /* note this by-passes the inner KSP and its options completely */
737: KSPGetPC(jac->ksp[0],&subpc);
738: PCApplySymmetricRight(subpc,bjac->x,bjac->y);
740: VecRestoreArrayRead(x,&x_array);
741: VecRestoreArray(y,&y_array);
742: return(0);
743: }
747: static PetscErrorCode PCApplyTranspose_BJacobi_Singleblock(PC pc,Vec x,Vec y)
748: {
749: PetscErrorCode ierr;
750: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
751: PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock*)jac->data;
752: PetscScalar *y_array;
753: const PetscScalar *x_array;
756: /*
757: The VecPlaceArray() is to avoid having to copy the
758: y vector into the bjac->x vector. The reason for
759: the bjac->x vector is that we need a sequential vector
760: for the sequential solve.
761: */
762: VecGetArrayRead(x,&x_array);
763: VecGetArray(y,&y_array);
764: VecPlaceArray(bjac->x,x_array);
765: VecPlaceArray(bjac->y,y_array);
766: KSPSolveTranspose(jac->ksp[0],bjac->x,bjac->y);
767: VecResetArray(bjac->x);
768: VecResetArray(bjac->y);
769: VecRestoreArrayRead(x,&x_array);
770: VecRestoreArray(y,&y_array);
771: return(0);
772: }
776: static PetscErrorCode PCSetUp_BJacobi_Singleblock(PC pc,Mat mat,Mat pmat)
777: {
778: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
779: PetscErrorCode ierr;
780: PetscInt m;
781: KSP ksp;
782: PC_BJacobi_Singleblock *bjac;
783: PetscBool wasSetup = PETSC_TRUE;
786: if (!pc->setupcalled) {
787: const char *prefix;
789: if (!jac->ksp) {
790: wasSetup = PETSC_FALSE;
792: KSPCreate(PETSC_COMM_SELF,&ksp);
793: KSPSetErrorIfNotConverged(ksp,pc->erroriffailure);
794: PetscObjectIncrementTabLevel((PetscObject)ksp,(PetscObject)pc,1);
795: PetscLogObjectParent((PetscObject)pc,(PetscObject)ksp);
796: KSPSetType(ksp,KSPPREONLY);
797: PCGetOptionsPrefix(pc,&prefix);
798: KSPSetOptionsPrefix(ksp,prefix);
799: KSPAppendOptionsPrefix(ksp,"sub_");
801: pc->ops->reset = PCReset_BJacobi_Singleblock;
802: pc->ops->destroy = PCDestroy_BJacobi_Singleblock;
803: pc->ops->apply = PCApply_BJacobi_Singleblock;
804: pc->ops->applysymmetricleft = PCApplySymmetricLeft_BJacobi_Singleblock;
805: pc->ops->applysymmetricright = PCApplySymmetricRight_BJacobi_Singleblock;
806: pc->ops->applytranspose = PCApplyTranspose_BJacobi_Singleblock;
807: pc->ops->setuponblocks = PCSetUpOnBlocks_BJacobi_Singleblock;
809: PetscMalloc1(1,&jac->ksp);
810: jac->ksp[0] = ksp;
812: PetscNewLog(pc,&bjac);
813: jac->data = (void*)bjac;
814: } else {
815: ksp = jac->ksp[0];
816: bjac = (PC_BJacobi_Singleblock*)jac->data;
817: }
819: /*
820: The reason we need to generate these vectors is to serve
821: as the right-hand side and solution vector for the solve on the
822: block. We do not need to allocate space for the vectors since
823: that is provided via VecPlaceArray() just before the call to
824: KSPSolve() on the block.
825: */
826: MatGetSize(pmat,&m,&m);
827: VecCreateSeqWithArray(PETSC_COMM_SELF,1,m,NULL,&bjac->x);
828: VecCreateSeqWithArray(PETSC_COMM_SELF,1,m,NULL,&bjac->y);
829: #if defined(PETSC_HAVE_CUSP)
830: VecSetType(bjac->x,VECCUSP);
831: VecSetType(bjac->y,VECCUSP);
832: #elif defined(PETSC_HAVE_VECCUDA)
833: VecSetType(bjac->x,VECCUDA);
834: VecSetType(bjac->y,VECCUDA);
835: #endif
836: PetscLogObjectParent((PetscObject)pc,(PetscObject)bjac->x);
837: PetscLogObjectParent((PetscObject)pc,(PetscObject)bjac->y);
838: } else {
839: ksp = jac->ksp[0];
840: bjac = (PC_BJacobi_Singleblock*)jac->data;
841: }
842: if (pc->useAmat) {
843: KSPSetOperators(ksp,mat,pmat);
844: } else {
845: KSPSetOperators(ksp,pmat,pmat);
846: }
847: if (!wasSetup && pc->setfromoptionscalled) {
848: /* If PCSetFromOptions_BJacobi is called later, KSPSetFromOptions will be called at that time. */
849: KSPSetFromOptions(ksp);
850: }
851: return(0);
852: }
854: /* ---------------------------------------------------------------------------------------------*/
857: static PetscErrorCode PCReset_BJacobi_Multiblock(PC pc)
858: {
859: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
860: PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock*)jac->data;
861: PetscErrorCode ierr;
862: PetscInt i;
865: if (bjac && bjac->pmat) {
866: MatDestroyMatrices(jac->n_local,&bjac->pmat);
867: if (pc->useAmat) {
868: MatDestroyMatrices(jac->n_local,&bjac->mat);
869: }
870: }
872: for (i=0; i<jac->n_local; i++) {
873: KSPReset(jac->ksp[i]);
874: if (bjac && bjac->x) {
875: VecDestroy(&bjac->x[i]);
876: VecDestroy(&bjac->y[i]);
877: ISDestroy(&bjac->is[i]);
878: }
879: }
880: PetscFree(jac->l_lens);
881: PetscFree(jac->g_lens);
882: return(0);
883: }
887: static PetscErrorCode PCDestroy_BJacobi_Multiblock(PC pc)
888: {
889: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
890: PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock*)jac->data;
891: PetscErrorCode ierr;
892: PetscInt i;
895: PCReset_BJacobi_Multiblock(pc);
896: if (bjac) {
897: PetscFree2(bjac->x,bjac->y);
898: PetscFree(bjac->starts);
899: PetscFree(bjac->is);
900: }
901: PetscFree(jac->data);
902: for (i=0; i<jac->n_local; i++) {
903: KSPDestroy(&jac->ksp[i]);
904: }
905: PetscFree(jac->ksp);
906: PetscFree(pc->data);
907: return(0);
908: }
912: static PetscErrorCode PCSetUpOnBlocks_BJacobi_Multiblock(PC pc)
913: {
914: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
916: PetscInt i,n_local = jac->n_local;
919: for (i=0; i<n_local; i++) {
920: KSPSetUp(jac->ksp[i]);
921: if (jac->ksp[i]->reason == KSP_DIVERGED_PCSETUP_FAILED) {
922: pc->failedreason = PC_SUBPC_ERROR;
923: }
924: }
925: return(0);
926: }
928: /*
929: Preconditioner for block Jacobi
930: */
933: static PetscErrorCode PCApply_BJacobi_Multiblock(PC pc,Vec x,Vec y)
934: {
935: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
936: PetscErrorCode ierr;
937: PetscInt i,n_local = jac->n_local;
938: PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock*)jac->data;
939: PetscScalar *yin;
940: const PetscScalar *xin;
943: VecGetArrayRead(x,&xin);
944: VecGetArray(y,&yin);
945: for (i=0; i<n_local; i++) {
946: /*
947: To avoid copying the subvector from x into a workspace we instead
948: make the workspace vector array point to the subpart of the array of
949: the global vector.
950: */
951: VecPlaceArray(bjac->x[i],xin+bjac->starts[i]);
952: VecPlaceArray(bjac->y[i],yin+bjac->starts[i]);
954: PetscLogEventBegin(PC_ApplyOnBlocks,jac->ksp[i],bjac->x[i],bjac->y[i],0);
955: KSPSolve(jac->ksp[i],bjac->x[i],bjac->y[i]);
956: PetscLogEventEnd(PC_ApplyOnBlocks,jac->ksp[i],bjac->x[i],bjac->y[i],0);
958: VecResetArray(bjac->x[i]);
959: VecResetArray(bjac->y[i]);
960: }
961: VecRestoreArrayRead(x,&xin);
962: VecRestoreArray(y,&yin);
963: return(0);
964: }
966: /*
967: Preconditioner for block Jacobi
968: */
971: static PetscErrorCode PCApplyTranspose_BJacobi_Multiblock(PC pc,Vec x,Vec y)
972: {
973: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
974: PetscErrorCode ierr;
975: PetscInt i,n_local = jac->n_local;
976: PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock*)jac->data;
977: PetscScalar *yin;
978: const PetscScalar *xin;
981: VecGetArrayRead(x,&xin);
982: VecGetArray(y,&yin);
983: for (i=0; i<n_local; i++) {
984: /*
985: To avoid copying the subvector from x into a workspace we instead
986: make the workspace vector array point to the subpart of the array of
987: the global vector.
988: */
989: VecPlaceArray(bjac->x[i],xin+bjac->starts[i]);
990: VecPlaceArray(bjac->y[i],yin+bjac->starts[i]);
992: PetscLogEventBegin(PC_ApplyTransposeOnBlocks,jac->ksp[i],bjac->x[i],bjac->y[i],0);
993: KSPSolveTranspose(jac->ksp[i],bjac->x[i],bjac->y[i]);
994: PetscLogEventEnd(PC_ApplyTransposeOnBlocks,jac->ksp[i],bjac->x[i],bjac->y[i],0);
996: VecResetArray(bjac->x[i]);
997: VecResetArray(bjac->y[i]);
998: }
999: VecRestoreArrayRead(x,&xin);
1000: VecRestoreArray(y,&yin);
1001: return(0);
1002: }
1006: static PetscErrorCode PCSetUp_BJacobi_Multiblock(PC pc,Mat mat,Mat pmat)
1007: {
1008: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
1009: PetscErrorCode ierr;
1010: PetscInt m,n_local,N,M,start,i;
1011: const char *prefix,*pprefix,*mprefix;
1012: KSP ksp;
1013: Vec x,y;
1014: PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock*)jac->data;
1015: PC subpc;
1016: IS is;
1017: MatReuse scall;
1020: MatGetLocalSize(pc->pmat,&M,&N);
1022: n_local = jac->n_local;
1024: if (pc->useAmat) {
1025: PetscBool same;
1026: PetscObjectTypeCompare((PetscObject)mat,((PetscObject)pmat)->type_name,&same);
1027: if (!same) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"Matrices not of same type");
1028: }
1030: if (!pc->setupcalled) {
1031: scall = MAT_INITIAL_MATRIX;
1033: if (!jac->ksp) {
1034: pc->ops->reset = PCReset_BJacobi_Multiblock;
1035: pc->ops->destroy = PCDestroy_BJacobi_Multiblock;
1036: pc->ops->apply = PCApply_BJacobi_Multiblock;
1037: pc->ops->applytranspose= PCApplyTranspose_BJacobi_Multiblock;
1038: pc->ops->setuponblocks = PCSetUpOnBlocks_BJacobi_Multiblock;
1040: PetscNewLog(pc,&bjac);
1041: PetscMalloc1(n_local,&jac->ksp);
1042: PetscLogObjectMemory((PetscObject)pc,sizeof(n_local*sizeof(KSP)));
1043: PetscMalloc2(n_local,&bjac->x,n_local,&bjac->y);
1044: PetscMalloc1(n_local,&bjac->starts);
1045: PetscLogObjectMemory((PetscObject)pc,sizeof(n_local*sizeof(PetscScalar)));
1047: jac->data = (void*)bjac;
1048: PetscMalloc1(n_local,&bjac->is);
1049: PetscLogObjectMemory((PetscObject)pc,sizeof(n_local*sizeof(IS)));
1051: for (i=0; i<n_local; i++) {
1052: KSPCreate(PETSC_COMM_SELF,&ksp);
1053: KSPSetErrorIfNotConverged(ksp,pc->erroriffailure);
1054: PetscObjectIncrementTabLevel((PetscObject)ksp,(PetscObject)pc,1);
1055: PetscLogObjectParent((PetscObject)pc,(PetscObject)ksp);
1056: KSPSetType(ksp,KSPPREONLY);
1057: KSPGetPC(ksp,&subpc);
1058: PCGetOptionsPrefix(pc,&prefix);
1059: KSPSetOptionsPrefix(ksp,prefix);
1060: KSPAppendOptionsPrefix(ksp,"sub_");
1062: jac->ksp[i] = ksp;
1063: }
1064: } else {
1065: bjac = (PC_BJacobi_Multiblock*)jac->data;
1066: }
1068: start = 0;
1069: for (i=0; i<n_local; i++) {
1070: m = jac->l_lens[i];
1071: /*
1072: The reason we need to generate these vectors is to serve
1073: as the right-hand side and solution vector for the solve on the
1074: block. We do not need to allocate space for the vectors since
1075: that is provided via VecPlaceArray() just before the call to
1076: KSPSolve() on the block.
1078: */
1079: VecCreateSeq(PETSC_COMM_SELF,m,&x);
1080: VecCreateSeqWithArray(PETSC_COMM_SELF,1,m,NULL,&y);
1081: #if defined(PETSC_HAVE_CUSP)
1082: VecSetType(x,VECCUSP);
1083: VecSetType(y,VECCUSP);
1084: #elif defined(PETSC_HAVE_VECCUDA)
1085: VecSetType(x,VECCUDA);
1086: VecSetType(y,VECCUDA);
1087: #endif
1088: PetscLogObjectParent((PetscObject)pc,(PetscObject)x);
1089: PetscLogObjectParent((PetscObject)pc,(PetscObject)y);
1091: bjac->x[i] = x;
1092: bjac->y[i] = y;
1093: bjac->starts[i] = start;
1095: ISCreateStride(PETSC_COMM_SELF,m,start,1,&is);
1096: bjac->is[i] = is;
1097: PetscLogObjectParent((PetscObject)pc,(PetscObject)is);
1099: start += m;
1100: }
1101: } else {
1102: bjac = (PC_BJacobi_Multiblock*)jac->data;
1103: /*
1104: Destroy the blocks from the previous iteration
1105: */
1106: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
1107: MatDestroyMatrices(n_local,&bjac->pmat);
1108: if (pc->useAmat) {
1109: MatDestroyMatrices(n_local,&bjac->mat);
1110: }
1111: scall = MAT_INITIAL_MATRIX;
1112: } else scall = MAT_REUSE_MATRIX;
1113: }
1115: MatGetSubMatrices(pmat,n_local,bjac->is,bjac->is,scall,&bjac->pmat);
1116: if (pc->useAmat) {
1117: PetscObjectGetOptionsPrefix((PetscObject)mat,&mprefix);
1118: MatGetSubMatrices(mat,n_local,bjac->is,bjac->is,scall,&bjac->mat);
1119: }
1120: /* Return control to the user so that the submatrices can be modified (e.g., to apply
1121: different boundary conditions for the submatrices than for the global problem) */
1122: PCModifySubMatrices(pc,n_local,bjac->is,bjac->is,bjac->pmat,pc->modifysubmatricesP);
1124: PetscObjectGetOptionsPrefix((PetscObject)pmat,&pprefix);
1125: for (i=0; i<n_local; i++) {
1126: PetscLogObjectParent((PetscObject)pc,(PetscObject)bjac->pmat[i]);
1127: PetscObjectSetOptionsPrefix((PetscObject)bjac->pmat[i],pprefix);
1128: if (pc->useAmat) {
1129: PetscLogObjectParent((PetscObject)pc,(PetscObject)bjac->mat[i]);
1130: PetscObjectSetOptionsPrefix((PetscObject)bjac->mat[i],mprefix);
1131: KSPSetOperators(jac->ksp[i],bjac->mat[i],bjac->pmat[i]);
1132: } else {
1133: KSPSetOperators(jac->ksp[i],bjac->pmat[i],bjac->pmat[i]);
1134: }
1135: if (pc->setfromoptionscalled) {
1136: KSPSetFromOptions(jac->ksp[i]);
1137: }
1138: }
1139: return(0);
1140: }
1142: /* ---------------------------------------------------------------------------------------------*/
1143: /*
1144: These are for a single block with multiple processes;
1145: */
1148: static PetscErrorCode PCReset_BJacobi_Multiproc(PC pc)
1149: {
1150: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
1151: PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc*)jac->data;
1152: PetscErrorCode ierr;
1155: VecDestroy(&mpjac->ysub);
1156: VecDestroy(&mpjac->xsub);
1157: MatDestroy(&mpjac->submats);
1158: if (jac->ksp) {KSPReset(jac->ksp[0]);}
1159: return(0);
1160: }
1164: static PetscErrorCode PCDestroy_BJacobi_Multiproc(PC pc)
1165: {
1166: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
1167: PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc*)jac->data;
1168: PetscErrorCode ierr;
1171: PCReset_BJacobi_Multiproc(pc);
1172: KSPDestroy(&jac->ksp[0]);
1173: PetscFree(jac->ksp);
1174: PetscSubcommDestroy(&mpjac->psubcomm);
1176: PetscFree(mpjac);
1177: PetscFree(pc->data);
1178: return(0);
1179: }
1183: static PetscErrorCode PCApply_BJacobi_Multiproc(PC pc,Vec x,Vec y)
1184: {
1185: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
1186: PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc*)jac->data;
1187: PetscErrorCode ierr;
1188: PetscScalar *yarray;
1189: const PetscScalar *xarray;
1192: /* place x's and y's local arrays into xsub and ysub */
1193: VecGetArrayRead(x,&xarray);
1194: VecGetArray(y,&yarray);
1195: VecPlaceArray(mpjac->xsub,xarray);
1196: VecPlaceArray(mpjac->ysub,yarray);
1198: /* apply preconditioner on each matrix block */
1199: PetscLogEventBegin(PC_ApplyOnMproc,jac->ksp[0],mpjac->xsub,mpjac->ysub,0);
1200: KSPSolve(jac->ksp[0],mpjac->xsub,mpjac->ysub);
1201: PetscLogEventEnd(PC_ApplyOnMproc,jac->ksp[0],mpjac->xsub,mpjac->ysub,0);
1203: if (jac->ksp[0]->reason == KSP_DIVERGED_PCSETUP_FAILED) {
1204: pc->failedreason = PC_SUBPC_ERROR;
1205: }
1207: VecResetArray(mpjac->xsub);
1208: VecResetArray(mpjac->ysub);
1209: VecRestoreArrayRead(x,&xarray);
1210: VecRestoreArray(y,&yarray);
1211: return(0);
1212: }
1214: #include <petsc/private/matimpl.h>
1217: static PetscErrorCode PCSetUp_BJacobi_Multiproc(PC pc)
1218: {
1219: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
1220: PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc*)jac->data;
1221: PetscErrorCode ierr;
1222: PetscInt m,n;
1223: MPI_Comm comm,subcomm=0;
1224: const char *prefix;
1225: PetscBool wasSetup = PETSC_TRUE;
1228: PetscObjectGetComm((PetscObject)pc,&comm);
1229: if (jac->n_local > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Only a single block in a subcommunicator is supported");
1230: jac->n_local = 1; /* currently only a single block is supported for a subcommunicator */
1231: if (!pc->setupcalled) {
1232: wasSetup = PETSC_FALSE;
1233: PetscNewLog(pc,&mpjac);
1234: jac->data = (void*)mpjac;
1236: /* initialize datastructure mpjac */
1237: if (!jac->psubcomm) {
1238: /* Create default contiguous subcommunicatiors if user does not provide them */
1239: PetscSubcommCreate(comm,&jac->psubcomm);
1240: PetscSubcommSetNumber(jac->psubcomm,jac->n);
1241: PetscSubcommSetType(jac->psubcomm,PETSC_SUBCOMM_CONTIGUOUS);
1242: PetscLogObjectMemory((PetscObject)pc,sizeof(PetscSubcomm));
1243: }
1244: mpjac->psubcomm = jac->psubcomm;
1245: subcomm = PetscSubcommChild(mpjac->psubcomm);
1247: /* Get matrix blocks of pmat */
1248: if (!pc->pmat->ops->getmultiprocblock) SETERRQ(PetscObjectComm((PetscObject)pc->pmat),PETSC_ERR_SUP,"No support for the requested operation");
1249: (*pc->pmat->ops->getmultiprocblock)(pc->pmat,subcomm,MAT_INITIAL_MATRIX,&mpjac->submats);
1251: /* create a new PC that processors in each subcomm have copy of */
1252: PetscMalloc1(1,&jac->ksp);
1253: KSPCreate(subcomm,&jac->ksp[0]);
1254: KSPSetErrorIfNotConverged(jac->ksp[0],pc->erroriffailure);
1255: PetscObjectIncrementTabLevel((PetscObject)jac->ksp[0],(PetscObject)pc,1);
1256: PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->ksp[0]);
1257: KSPSetOperators(jac->ksp[0],mpjac->submats,mpjac->submats);
1258: KSPGetPC(jac->ksp[0],&mpjac->pc);
1260: PCGetOptionsPrefix(pc,&prefix);
1261: KSPSetOptionsPrefix(jac->ksp[0],prefix);
1262: KSPAppendOptionsPrefix(jac->ksp[0],"sub_");
1263: /*
1264: PetscMPIInt rank,subsize,subrank;
1265: MPI_Comm_rank(comm,&rank);
1266: MPI_Comm_size(subcomm,&subsize);
1267: MPI_Comm_rank(subcomm,&subrank);
1269: MatGetLocalSize(mpjac->submats,&m,NULL);
1270: MatGetSize(mpjac->submats,&n,NULL);
1271: PetscSynchronizedPrintf(comm,"[%d], sub-size %d,sub-rank %d\n",rank,subsize,subrank);
1272: PetscSynchronizedFlush(comm,PETSC_STDOUT);
1273: */
1275: /* create dummy vectors xsub and ysub */
1276: MatGetLocalSize(mpjac->submats,&m,&n);
1277: VecCreateMPIWithArray(subcomm,1,n,PETSC_DECIDE,NULL,&mpjac->xsub);
1278: VecCreateMPIWithArray(subcomm,1,m,PETSC_DECIDE,NULL,&mpjac->ysub);
1279: #if defined(PETSC_HAVE_CUSP)
1280: VecSetType(mpjac->xsub,VECMPICUSP);
1281: VecSetType(mpjac->ysub,VECMPICUSP);
1282: #elif defined(PETSC_HAVE_VECCUDA)
1283: VecSetType(mpjac->xsub,VECMPICUDA);
1284: VecSetType(mpjac->ysub,VECMPICUDA);
1285: #endif
1286: PetscLogObjectParent((PetscObject)pc,(PetscObject)mpjac->xsub);
1287: PetscLogObjectParent((PetscObject)pc,(PetscObject)mpjac->ysub);
1289: pc->ops->reset = PCReset_BJacobi_Multiproc;
1290: pc->ops->destroy = PCDestroy_BJacobi_Multiproc;
1291: pc->ops->apply = PCApply_BJacobi_Multiproc;
1292: } else { /* pc->setupcalled */
1293: subcomm = PetscSubcommChild(mpjac->psubcomm);
1294: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
1295: /* destroy old matrix blocks, then get new matrix blocks */
1296: if (mpjac->submats) {MatDestroy(&mpjac->submats);}
1297: (*pc->pmat->ops->getmultiprocblock)(pc->pmat,subcomm,MAT_INITIAL_MATRIX,&mpjac->submats);
1298: } else {
1299: (*pc->pmat->ops->getmultiprocblock)(pc->pmat,subcomm,MAT_REUSE_MATRIX,&mpjac->submats);
1300: }
1301: KSPSetOperators(jac->ksp[0],mpjac->submats,mpjac->submats);
1302: }
1304: if (!wasSetup && pc->setfromoptionscalled) {
1305: KSPSetFromOptions(jac->ksp[0]);
1306: }
1307: return(0);
1308: }