Actual source code: bjacobi.c
petsc-3.4.5 2014-06-29
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: MPI_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: PetscMalloc(jac->n_local*sizeof(PetscInt),&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: PetscMalloc(jac->n_local*sizeof(PetscInt),&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: PetscMalloc(jac->n_local*sizeof(PetscInt),&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: PetscMalloc(jac->n_local*sizeof(PetscInt),&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: PetscMalloc(sizeof(PetscInt),&jac->l_lens);
101: jac->l_lens[0] = M;
102: } else { /* jac->n > 0 && jac->n_local > 0 */
103: if (!jac->l_lens) {
104: PetscMalloc(jac->n_local*sizeof(PetscInt),&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(PC pc)
163: {
164: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
166: PetscInt blocks,i;
167: PetscBool flg;
170: PetscOptionsHead("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: PetscViewerGetSingleton(viewer,&sviewer);
213: if (!rank) {
214: PetscViewerASCIIPushTab(viewer);
215: KSPView(jac->ksp[0],sviewer);
216: PetscViewerASCIIPopTab(viewer);
217: }
218: PetscViewerRestoreSingleton(viewer,&sviewer);
219: } else if (jac->psubcomm && !jac->psubcomm->color) {
220: PetscViewerASCIIGetStdout(mpjac->psubcomm->comm,&sviewer);
221: PetscViewerASCIIPushTab(viewer);
222: KSPView(*(jac->ksp),sviewer);
223: PetscViewerASCIIPopTab(viewer);
224: }
225: } else {
226: PetscInt n_global;
227: MPI_Allreduce(&jac->n_local,&n_global,1,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)pc));
228: PetscViewerASCIISynchronizedAllow(viewer,PETSC_TRUE);
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: for (i=0; i<n_global; i++) {
234: PetscViewerGetSingleton(viewer,&sviewer);
235: if (i < jac->n_local) {
236: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] local block number %D\n",rank,i);
237: KSPView(jac->ksp[i],sviewer);
238: PetscViewerASCIISynchronizedPrintf(viewer,"- - - - - - - - - - - - - - - - - -\n");
239: }
240: PetscViewerRestoreSingleton(viewer,&sviewer);
241: }
242: PetscViewerASCIIPopTab(viewer);
243: PetscViewerFlush(viewer);
244: PetscViewerASCIISynchronizedAllow(viewer,PETSC_FALSE);
245: }
246: } else if (isstring) {
247: PetscViewerStringSPrintf(viewer," blks=%D",jac->n);
248: PetscViewerGetSingleton(viewer,&sviewer);
249: if (jac->ksp) {KSPView(jac->ksp[0],sviewer);}
250: PetscViewerRestoreSingleton(viewer,&sviewer);
251: } else if (isdraw) {
252: PetscDraw draw;
253: char str[25];
254: PetscReal x,y,bottom,h;
256: PetscViewerDrawGetDraw(viewer,0,&draw);
257: PetscDrawGetCurrentPoint(draw,&x,&y);
258: PetscSNPrintf(str,25,"Number blocks %D",jac->n);
259: PetscDrawBoxedString(draw,x,y,PETSC_DRAW_RED,PETSC_DRAW_BLACK,str,NULL,&h);
260: bottom = y - h;
261: PetscDrawPushCurrentPoint(draw,x,bottom);
262: /* warning the communicator on viewer is different then on ksp in parallel */
263: if (jac->ksp) {KSPView(jac->ksp[0],viewer);}
264: PetscDrawPopCurrentPoint(draw);
265: }
266: return(0);
267: }
269: /* -------------------------------------------------------------------------------------*/
273: static PetscErrorCode PCBJacobiGetSubKSP_BJacobi(PC pc,PetscInt *n_local,PetscInt *first_local,KSP **ksp)
274: {
275: PC_BJacobi *jac = (PC_BJacobi*)pc->data;;
278: if (!pc->setupcalled) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must call KSPSetUp() or PCSetUp() first");
280: if (n_local) *n_local = jac->n_local;
281: if (first_local) *first_local = jac->first_local;
282: *ksp = jac->ksp;
283: jac->same_local_solves = PETSC_FALSE; /* Assume that local solves are now different;
284: not necessarily true though! This flag is
285: used only for PCView_BJacobi() */
286: return(0);
287: }
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 = 0;
300: else {
301: PetscMalloc(blocks*sizeof(PetscInt),&jac->g_lens);
302: PetscLogObjectMemory(pc,blocks*sizeof(PetscInt));
303: PetscMemcpy(jac->g_lens,lens,blocks*sizeof(PetscInt));
304: }
305: return(0);
306: }
310: static PetscErrorCode PCBJacobiGetTotalBlocks_BJacobi(PC pc, PetscInt *blocks, const PetscInt *lens[])
311: {
312: PC_BJacobi *jac = (PC_BJacobi*) pc->data;
315: *blocks = jac->n;
316: if (lens) *lens = jac->g_lens;
317: return(0);
318: }
322: static PetscErrorCode PCBJacobiSetLocalBlocks_BJacobi(PC pc,PetscInt blocks,const PetscInt lens[])
323: {
324: PC_BJacobi *jac;
328: jac = (PC_BJacobi*)pc->data;
330: jac->n_local = blocks;
331: if (!lens) jac->l_lens = 0;
332: else {
333: PetscMalloc(blocks*sizeof(PetscInt),&jac->l_lens);
334: PetscLogObjectMemory(pc,blocks*sizeof(PetscInt));
335: PetscMemcpy(jac->l_lens,lens,blocks*sizeof(PetscInt));
336: }
337: return(0);
338: }
342: static PetscErrorCode PCBJacobiGetLocalBlocks_BJacobi(PC pc, PetscInt *blocks, const PetscInt *lens[])
343: {
344: PC_BJacobi *jac = (PC_BJacobi*) pc->data;
347: *blocks = jac->n_local;
348: if (lens) *lens = jac->l_lens;
349: return(0);
350: }
352: /* -------------------------------------------------------------------------------------*/
356: /*@C
357: PCBJacobiGetSubKSP - Gets the local KSP contexts for all blocks on
358: this processor.
360: Note Collective
362: Input Parameter:
363: . pc - the preconditioner context
365: Output Parameters:
366: + n_local - the number of blocks on this processor, or NULL
367: . first_local - the global number of the first block on this processor, or NULL
368: - ksp - the array of KSP contexts
370: Notes:
371: After PCBJacobiGetSubKSP() the array of KSP contexts is not to be freed.
373: Currently for some matrix implementations only 1 block per processor
374: is supported.
376: You must call KSPSetUp() or PCSetUp() before calling PCBJacobiGetSubKSP().
378: Fortran Usage: You must pass in a KSP array that is large enough to contain all the local KSPs.
379: You can call PCBJacobiGetSubKSP(pc,nlocal,firstlocal,NULL_OBJECT,ierr) to determine how large the
380: KSP array must be.
382: Level: advanced
384: .keywords: block, Jacobi, get, sub, KSP, context
386: .seealso: PCBJacobiGetSubKSP()
387: @*/
388: PetscErrorCode PCBJacobiGetSubKSP(PC pc,PetscInt *n_local,PetscInt *first_local,KSP *ksp[])
389: {
394: PetscUseMethod(pc,"PCBJacobiGetSubKSP_C",(PC,PetscInt*,PetscInt*,KSP **),(pc,n_local,first_local,ksp));
395: return(0);
396: }
400: /*@
401: PCBJacobiSetTotalBlocks - Sets the global number of blocks for the block
402: Jacobi preconditioner.
404: Collective on PC
406: Input Parameters:
407: + pc - the preconditioner context
408: . blocks - the number of blocks
409: - lens - [optional] integer array containing the size of each block
411: Options Database Key:
412: . -pc_bjacobi_blocks <blocks> - Sets the number of global blocks
414: Notes:
415: Currently only a limited number of blocking configurations are supported.
416: All processors sharing the PC must call this routine with the same data.
418: Level: intermediate
420: .keywords: set, number, Jacobi, global, total, blocks
422: .seealso: PCSetUseAmat(), PCBJacobiSetLocalBlocks()
423: @*/
424: PetscErrorCode PCBJacobiSetTotalBlocks(PC pc,PetscInt blocks,const PetscInt lens[])
425: {
430: if (blocks <= 0) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Must have positive blocks");
431: PetscTryMethod(pc,"PCBJacobiSetTotalBlocks_C",(PC,PetscInt,const PetscInt[]),(pc,blocks,lens));
432: return(0);
433: }
437: /*@C
438: PCBJacobiGetTotalBlocks - Gets the global number of blocks for the block
439: Jacobi preconditioner.
441: Not Collective
443: Input Parameter:
444: . pc - the preconditioner context
446: Output parameters:
447: + blocks - the number of blocks
448: - lens - integer array containing the size of each block
450: Level: intermediate
452: .keywords: get, number, Jacobi, global, total, blocks
454: .seealso: PCSetUseAmat(), PCBJacobiGetLocalBlocks()
455: @*/
456: PetscErrorCode PCBJacobiGetTotalBlocks(PC pc, PetscInt *blocks, const PetscInt *lens[])
457: {
463: PetscUseMethod(pc,"PCBJacobiGetTotalBlocks_C",(PC,PetscInt*, const PetscInt *[]),(pc,blocks,lens));
464: return(0);
465: }
469: /*@
470: PCBJacobiSetLocalBlocks - Sets the local number of blocks for the block
471: Jacobi preconditioner.
473: Not Collective
475: Input Parameters:
476: + pc - the preconditioner context
477: . blocks - the number of blocks
478: - lens - [optional] integer array containing size of each block
480: Note:
481: Currently only a limited number of blocking configurations are supported.
483: Level: intermediate
485: .keywords: PC, set, number, Jacobi, local, blocks
487: .seealso: PCSetUseAmat(), PCBJacobiSetTotalBlocks()
488: @*/
489: PetscErrorCode PCBJacobiSetLocalBlocks(PC pc,PetscInt blocks,const PetscInt lens[])
490: {
495: if (blocks < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must have nonegative blocks");
496: PetscTryMethod(pc,"PCBJacobiSetLocalBlocks_C",(PC,PetscInt,const PetscInt []),(pc,blocks,lens));
497: return(0);
498: }
502: /*@C
503: PCBJacobiGetLocalBlocks - Gets the local number of blocks for the block
504: Jacobi preconditioner.
506: Not Collective
508: Input Parameters:
509: + pc - the preconditioner context
510: . blocks - the number of blocks
511: - lens - [optional] integer array containing size of each block
513: Note:
514: Currently only a limited number of blocking configurations are supported.
516: Level: intermediate
518: .keywords: PC, get, number, Jacobi, local, blocks
520: .seealso: PCSetUseAmat(), PCBJacobiGetTotalBlocks()
521: @*/
522: PetscErrorCode PCBJacobiGetLocalBlocks(PC pc, PetscInt *blocks, const PetscInt *lens[])
523: {
529: PetscUseMethod(pc,"PCBJacobiGetLocalBlocks_C",(PC,PetscInt*, const PetscInt *[]),(pc,blocks,lens));
530: return(0);
531: }
533: /* -----------------------------------------------------------------------------------*/
535: /*MC
536: PCBJACOBI - Use block Jacobi preconditioning, each block is (approximately) solved with
537: its own KSP object.
539: Options Database Keys:
540: . -pc_use_amat - use Amat to apply block of operator in inner Krylov method
542: Notes: Each processor can have one or more blocks, but a block cannot be shared by more
543: than one processor. Defaults to one block per processor.
545: To set options on the solvers for each block append -sub_ to all the KSP, KSP, and PC
546: options database keys. For example, -sub_pc_type ilu -sub_pc_factor_levels 1 -sub_ksp_type preonly
548: To set the options on the solvers separate for each block call PCBJacobiGetSubKSP()
549: and set the options directly on the resulting KSP object (you can access its PC
550: KSPGetPC())
552: Level: beginner
554: Concepts: block Jacobi
556: Developer Notes: This preconditioner does not currently work with CUDA/CUSP for a couple of reasons.
557: (1) It creates seq vectors as work vectors that should be cusp
558: (2) The use of VecPlaceArray() is not handled properly by CUSP (that is it will not know where
559: the ownership of the vector is so may use wrong values) even if it did know the ownership
560: it may induce extra copy ups and downs. Satish suggests a VecTransplantArray() to handle two
561: vectors sharing the same pointer and handling the CUSP side as well instead of VecGetArray()/VecPlaceArray().
564: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
565: PCASM, PCSetUseAmat(), PCGetUseAmat(), PCBJacobiGetSubKSP(), PCBJacobiSetTotalBlocks(),
566: PCBJacobiSetLocalBlocks(), PCSetModifySubmatrices()
567: M*/
571: PETSC_EXTERN PetscErrorCode PCCreate_BJacobi(PC pc)
572: {
574: PetscMPIInt rank;
575: PC_BJacobi *jac;
578: PetscNewLog(pc,PC_BJacobi,&jac);
579: MPI_Comm_rank(PetscObjectComm((PetscObject)pc),&rank);
581: pc->ops->apply = 0;
582: pc->ops->applytranspose = 0;
583: pc->ops->setup = PCSetUp_BJacobi;
584: pc->ops->destroy = PCDestroy_BJacobi;
585: pc->ops->setfromoptions = PCSetFromOptions_BJacobi;
586: pc->ops->view = PCView_BJacobi;
587: pc->ops->applyrichardson = 0;
589: pc->data = (void*)jac;
590: jac->n = -1;
591: jac->n_local = -1;
592: jac->first_local = rank;
593: jac->ksp = 0;
594: jac->same_local_solves = PETSC_TRUE;
595: jac->g_lens = 0;
596: jac->l_lens = 0;
597: jac->psubcomm = 0;
599: PetscObjectComposeFunction((PetscObject)pc,"PCBJacobiGetSubKSP_C",PCBJacobiGetSubKSP_BJacobi);
600: PetscObjectComposeFunction((PetscObject)pc,"PCBJacobiSetTotalBlocks_C",PCBJacobiSetTotalBlocks_BJacobi);
601: PetscObjectComposeFunction((PetscObject)pc,"PCBJacobiGetTotalBlocks_C",PCBJacobiGetTotalBlocks_BJacobi);
602: PetscObjectComposeFunction((PetscObject)pc,"PCBJacobiSetLocalBlocks_C",PCBJacobiSetLocalBlocks_BJacobi);
603: PetscObjectComposeFunction((PetscObject)pc,"PCBJacobiGetLocalBlocks_C",PCBJacobiGetLocalBlocks_BJacobi);
604: return(0);
605: }
607: /* --------------------------------------------------------------------------------------------*/
608: /*
609: These are for a single block per processor; works for AIJ, BAIJ; Seq and MPI
610: */
613: PetscErrorCode PCReset_BJacobi_Singleblock(PC pc)
614: {
615: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
616: PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock*)jac->data;
617: PetscErrorCode ierr;
620: KSPReset(jac->ksp[0]);
621: VecDestroy(&bjac->x);
622: VecDestroy(&bjac->y);
623: return(0);
624: }
628: PetscErrorCode PCDestroy_BJacobi_Singleblock(PC pc)
629: {
630: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
631: PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock*)jac->data;
632: PetscErrorCode ierr;
635: PCReset_BJacobi_Singleblock(pc);
636: KSPDestroy(&jac->ksp[0]);
637: PetscFree(jac->ksp);
638: PetscFree(jac->l_lens);
639: PetscFree(jac->g_lens);
640: PetscFree(bjac);
641: PetscFree(pc->data);
642: return(0);
643: }
647: PetscErrorCode PCSetUpOnBlocks_BJacobi_Singleblock(PC pc)
648: {
650: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
653: KSPSetUp(jac->ksp[0]);
654: return(0);
655: }
659: 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;
664: PetscScalar *x_array,*y_array;
667: /*
668: The VecPlaceArray() is to avoid having to copy the
669: y vector into the bjac->x vector. The reason for
670: the bjac->x vector is that we need a sequential vector
671: for the sequential solve.
672: */
673: VecGetArray(x,&x_array);
674: VecGetArray(y,&y_array);
675: VecPlaceArray(bjac->x,x_array);
676: VecPlaceArray(bjac->y,y_array);
677: KSPSolve(jac->ksp[0],bjac->x,bjac->y);
678: VecResetArray(bjac->x);
679: VecResetArray(bjac->y);
680: VecRestoreArray(x,&x_array);
681: VecRestoreArray(y,&y_array);
682: return(0);
683: }
687: PetscErrorCode PCApplySymmetricLeft_BJacobi_Singleblock(PC pc,Vec x,Vec y)
688: {
689: PetscErrorCode ierr;
690: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
691: PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock*)jac->data;
692: PetscScalar *x_array,*y_array;
693: PC subpc;
696: /*
697: The VecPlaceArray() is to avoid having to copy the
698: y vector into the bjac->x vector. The reason for
699: the bjac->x vector is that we need a sequential vector
700: for the sequential solve.
701: */
702: VecGetArray(x,&x_array);
703: VecGetArray(y,&y_array);
704: VecPlaceArray(bjac->x,x_array);
705: VecPlaceArray(bjac->y,y_array);
706: /* apply the symmetric left portion of the inner PC operator */
707: /* note this by-passes the inner KSP and its options completely */
708: KSPGetPC(jac->ksp[0],&subpc);
709: PCApplySymmetricLeft(subpc,bjac->x,bjac->y);
710: VecResetArray(bjac->x);
711: VecResetArray(bjac->y);
712: VecRestoreArray(x,&x_array);
713: VecRestoreArray(y,&y_array);
714: return(0);
715: }
719: PetscErrorCode PCApplySymmetricRight_BJacobi_Singleblock(PC pc,Vec x,Vec y)
720: {
721: PetscErrorCode ierr;
722: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
723: PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock*)jac->data;
724: PetscScalar *x_array,*y_array;
725: PC subpc;
728: /*
729: The VecPlaceArray() is to avoid having to copy the
730: y vector into the bjac->x vector. The reason for
731: the bjac->x vector is that we need a sequential vector
732: for the sequential solve.
733: */
734: VecGetArray(x,&x_array);
735: VecGetArray(y,&y_array);
736: VecPlaceArray(bjac->x,x_array);
737: VecPlaceArray(bjac->y,y_array);
739: /* apply the symmetric right portion of the inner PC operator */
740: /* note this by-passes the inner KSP and its options completely */
742: KSPGetPC(jac->ksp[0],&subpc);
743: PCApplySymmetricRight(subpc,bjac->x,bjac->y);
745: VecRestoreArray(x,&x_array);
746: VecRestoreArray(y,&y_array);
747: return(0);
748: }
752: PetscErrorCode PCApplyTranspose_BJacobi_Singleblock(PC pc,Vec x,Vec y)
753: {
754: PetscErrorCode ierr;
755: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
756: PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock*)jac->data;
757: PetscScalar *x_array,*y_array;
760: /*
761: The VecPlaceArray() is to avoid having to copy the
762: y vector into the bjac->x vector. The reason for
763: the bjac->x vector is that we need a sequential vector
764: for the sequential solve.
765: */
766: VecGetArray(x,&x_array);
767: VecGetArray(y,&y_array);
768: VecPlaceArray(bjac->x,x_array);
769: VecPlaceArray(bjac->y,y_array);
770: KSPSolveTranspose(jac->ksp[0],bjac->x,bjac->y);
771: VecResetArray(bjac->x);
772: VecResetArray(bjac->y);
773: VecRestoreArray(x,&x_array);
774: VecRestoreArray(y,&y_array);
775: return(0);
776: }
780: static PetscErrorCode PCSetUp_BJacobi_Singleblock(PC pc,Mat mat,Mat pmat)
781: {
782: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
783: PetscErrorCode ierr;
784: PetscInt m;
785: KSP ksp;
786: PC_BJacobi_Singleblock *bjac;
787: PetscBool wasSetup = PETSC_TRUE;
790: if (!pc->setupcalled) {
791: const char *prefix;
793: if (!jac->ksp) {
794: wasSetup = PETSC_FALSE;
796: KSPCreate(PETSC_COMM_SELF,&ksp);
797: PetscObjectIncrementTabLevel((PetscObject)ksp,(PetscObject)pc,1);
798: PetscLogObjectParent(pc,ksp);
799: KSPSetType(ksp,KSPPREONLY);
800: PCGetOptionsPrefix(pc,&prefix);
801: KSPSetOptionsPrefix(ksp,prefix);
802: KSPAppendOptionsPrefix(ksp,"sub_");
804: pc->ops->reset = PCReset_BJacobi_Singleblock;
805: pc->ops->destroy = PCDestroy_BJacobi_Singleblock;
806: pc->ops->apply = PCApply_BJacobi_Singleblock;
807: pc->ops->applysymmetricleft = PCApplySymmetricLeft_BJacobi_Singleblock;
808: pc->ops->applysymmetricright = PCApplySymmetricRight_BJacobi_Singleblock;
809: pc->ops->applytranspose = PCApplyTranspose_BJacobi_Singleblock;
810: pc->ops->setuponblocks = PCSetUpOnBlocks_BJacobi_Singleblock;
812: PetscMalloc(sizeof(KSP),&jac->ksp);
813: jac->ksp[0] = ksp;
815: PetscNewLog(pc,PC_BJacobi_Singleblock,&bjac);
816: jac->data = (void*)bjac;
817: } else {
818: ksp = jac->ksp[0];
819: bjac = (PC_BJacobi_Singleblock*)jac->data;
820: }
822: /*
823: The reason we need to generate these vectors is to serve
824: as the right-hand side and solution vector for the solve on the
825: block. We do not need to allocate space for the vectors since
826: that is provided via VecPlaceArray() just before the call to
827: KSPSolve() on the block.
828: */
829: MatGetSize(pmat,&m,&m);
830: VecCreateSeqWithArray(PETSC_COMM_SELF,1,m,NULL,&bjac->x);
831: VecCreateSeqWithArray(PETSC_COMM_SELF,1,m,NULL,&bjac->y);
832: PetscLogObjectParent(pc,bjac->x);
833: PetscLogObjectParent(pc,bjac->y);
834: } else {
835: ksp = jac->ksp[0];
836: bjac = (PC_BJacobi_Singleblock*)jac->data;
837: }
838: if (pc->useAmat) {
839: KSPSetOperators(ksp,mat,pmat,pc->flag);
840: } else {
841: KSPSetOperators(ksp,pmat,pmat,pc->flag);
842: }
843: if (!wasSetup && pc->setfromoptionscalled) {
844: /* If PCSetFromOptions_BJacobi is called later, KSPSetFromOptions will be called at that time. */
845: KSPSetFromOptions(ksp);
846: }
847: return(0);
848: }
850: /* ---------------------------------------------------------------------------------------------*/
853: PetscErrorCode PCReset_BJacobi_Multiblock(PC pc)
854: {
855: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
856: PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock*)jac->data;
857: PetscErrorCode ierr;
858: PetscInt i;
861: if (bjac && bjac->pmat) {
862: MatDestroyMatrices(jac->n_local,&bjac->pmat);
863: if (pc->useAmat) {
864: MatDestroyMatrices(jac->n_local,&bjac->mat);
865: }
866: }
868: for (i=0; i<jac->n_local; i++) {
869: KSPReset(jac->ksp[i]);
870: if (bjac && bjac->x) {
871: VecDestroy(&bjac->x[i]);
872: VecDestroy(&bjac->y[i]);
873: ISDestroy(&bjac->is[i]);
874: }
875: }
876: PetscFree(jac->l_lens);
877: PetscFree(jac->g_lens);
878: return(0);
879: }
883: PetscErrorCode PCDestroy_BJacobi_Multiblock(PC pc)
884: {
885: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
886: PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock*)jac->data;
887: PetscErrorCode ierr;
888: PetscInt i;
891: PCReset_BJacobi_Multiblock(pc);
892: if (bjac) {
893: PetscFree2(bjac->x,bjac->y);
894: PetscFree(bjac->starts);
895: PetscFree(bjac->is);
896: }
897: PetscFree(jac->data);
898: for (i=0; i<jac->n_local; i++) {
899: KSPDestroy(&jac->ksp[i]);
900: }
901: PetscFree(jac->ksp);
902: PetscFree(pc->data);
903: return(0);
904: }
908: PetscErrorCode PCSetUpOnBlocks_BJacobi_Multiblock(PC pc)
909: {
910: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
912: PetscInt i,n_local = jac->n_local;
915: for (i=0; i<n_local; i++) {
916: KSPSetUp(jac->ksp[i]);
917: }
918: return(0);
919: }
921: /*
922: Preconditioner for block Jacobi
923: */
926: PetscErrorCode PCApply_BJacobi_Multiblock(PC pc,Vec x,Vec y)
927: {
928: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
929: PetscErrorCode ierr;
930: PetscInt i,n_local = jac->n_local;
931: PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock*)jac->data;
932: PetscScalar *xin,*yin;
935: VecGetArray(x,&xin);
936: VecGetArray(y,&yin);
937: for (i=0; i<n_local; i++) {
938: /*
939: To avoid copying the subvector from x into a workspace we instead
940: make the workspace vector array point to the subpart of the array of
941: the global vector.
942: */
943: VecPlaceArray(bjac->x[i],xin+bjac->starts[i]);
944: VecPlaceArray(bjac->y[i],yin+bjac->starts[i]);
946: PetscLogEventBegin(PC_ApplyOnBlocks,jac->ksp[i],bjac->x[i],bjac->y[i],0);
947: KSPSolve(jac->ksp[i],bjac->x[i],bjac->y[i]);
948: PetscLogEventEnd(PC_ApplyOnBlocks,jac->ksp[i],bjac->x[i],bjac->y[i],0);
950: VecResetArray(bjac->x[i]);
951: VecResetArray(bjac->y[i]);
952: }
953: VecRestoreArray(x,&xin);
954: VecRestoreArray(y,&yin);
955: return(0);
956: }
958: /*
959: Preconditioner for block Jacobi
960: */
963: PetscErrorCode PCApplyTranspose_BJacobi_Multiblock(PC pc,Vec x,Vec y)
964: {
965: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
966: PetscErrorCode ierr;
967: PetscInt i,n_local = jac->n_local;
968: PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock*)jac->data;
969: PetscScalar *xin,*yin;
972: VecGetArray(x,&xin);
973: VecGetArray(y,&yin);
974: for (i=0; i<n_local; i++) {
975: /*
976: To avoid copying the subvector from x into a workspace we instead
977: make the workspace vector array point to the subpart of the array of
978: the global vector.
979: */
980: VecPlaceArray(bjac->x[i],xin+bjac->starts[i]);
981: VecPlaceArray(bjac->y[i],yin+bjac->starts[i]);
983: PetscLogEventBegin(PC_ApplyTransposeOnBlocks,jac->ksp[i],bjac->x[i],bjac->y[i],0);
984: KSPSolveTranspose(jac->ksp[i],bjac->x[i],bjac->y[i]);
985: PetscLogEventEnd(PC_ApplyTransposeOnBlocks,jac->ksp[i],bjac->x[i],bjac->y[i],0);
987: VecResetArray(bjac->x[i]);
988: VecResetArray(bjac->y[i]);
989: }
990: VecRestoreArray(x,&xin);
991: VecRestoreArray(y,&yin);
992: return(0);
993: }
997: static PetscErrorCode PCSetUp_BJacobi_Multiblock(PC pc,Mat mat,Mat pmat)
998: {
999: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
1000: PetscErrorCode ierr;
1001: PetscInt m,n_local,N,M,start,i;
1002: const char *prefix,*pprefix,*mprefix;
1003: KSP ksp;
1004: Vec x,y;
1005: PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock*)jac->data;
1006: PC subpc;
1007: IS is;
1008: MatReuse scall;
1011: MatGetLocalSize(pc->pmat,&M,&N);
1013: n_local = jac->n_local;
1015: if (pc->useAmat) {
1016: PetscBool same;
1017: PetscObjectTypeCompare((PetscObject)mat,((PetscObject)pmat)->type_name,&same);
1018: if (!same) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"Matrices not of same type");
1019: }
1021: if (!pc->setupcalled) {
1022: scall = MAT_INITIAL_MATRIX;
1024: if (!jac->ksp) {
1025: pc->ops->reset = PCReset_BJacobi_Multiblock;
1026: pc->ops->destroy = PCDestroy_BJacobi_Multiblock;
1027: pc->ops->apply = PCApply_BJacobi_Multiblock;
1028: pc->ops->applytranspose= PCApplyTranspose_BJacobi_Multiblock;
1029: pc->ops->setuponblocks = PCSetUpOnBlocks_BJacobi_Multiblock;
1031: PetscNewLog(pc,PC_BJacobi_Multiblock,&bjac);
1032: PetscMalloc(n_local*sizeof(KSP),&jac->ksp);
1033: PetscLogObjectMemory(pc,sizeof(n_local*sizeof(KSP)));
1034: PetscMalloc2(n_local,Vec,&bjac->x,n_local,Vec,&bjac->y);
1035: PetscMalloc(n_local*sizeof(PetscScalar),&bjac->starts);
1036: PetscLogObjectMemory(pc,sizeof(n_local*sizeof(PetscScalar)));
1038: jac->data = (void*)bjac;
1039: PetscMalloc(n_local*sizeof(IS),&bjac->is);
1040: PetscLogObjectMemory(pc,sizeof(n_local*sizeof(IS)));
1042: for (i=0; i<n_local; i++) {
1043: KSPCreate(PETSC_COMM_SELF,&ksp);
1044: PetscObjectIncrementTabLevel((PetscObject)ksp,(PetscObject)pc,1);
1045: PetscLogObjectParent(pc,ksp);
1046: KSPSetType(ksp,KSPPREONLY);
1047: KSPGetPC(ksp,&subpc);
1048: PCGetOptionsPrefix(pc,&prefix);
1049: KSPSetOptionsPrefix(ksp,prefix);
1050: KSPAppendOptionsPrefix(ksp,"sub_");
1052: jac->ksp[i] = ksp;
1053: }
1054: } else {
1055: bjac = (PC_BJacobi_Multiblock*)jac->data;
1056: }
1058: start = 0;
1059: for (i=0; i<n_local; i++) {
1060: m = jac->l_lens[i];
1061: /*
1062: The reason we need to generate these vectors is to serve
1063: as the right-hand side and solution vector for the solve on the
1064: block. We do not need to allocate space for the vectors since
1065: that is provided via VecPlaceArray() just before the call to
1066: KSPSolve() on the block.
1068: */
1069: VecCreateSeq(PETSC_COMM_SELF,m,&x);
1070: VecCreateSeqWithArray(PETSC_COMM_SELF,1,m,NULL,&y);
1071: PetscLogObjectParent(pc,x);
1072: PetscLogObjectParent(pc,y);
1074: bjac->x[i] = x;
1075: bjac->y[i] = y;
1076: bjac->starts[i] = start;
1078: ISCreateStride(PETSC_COMM_SELF,m,start,1,&is);
1079: bjac->is[i] = is;
1080: PetscLogObjectParent(pc,is);
1082: start += m;
1083: }
1084: } else {
1085: bjac = (PC_BJacobi_Multiblock*)jac->data;
1086: /*
1087: Destroy the blocks from the previous iteration
1088: */
1089: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
1090: MatDestroyMatrices(n_local,&bjac->pmat);
1091: if (pc->useAmat) {
1092: MatDestroyMatrices(n_local,&bjac->mat);
1093: }
1094: scall = MAT_INITIAL_MATRIX;
1095: } else scall = MAT_REUSE_MATRIX;
1096: }
1098: MatGetSubMatrices(pmat,n_local,bjac->is,bjac->is,scall,&bjac->pmat);
1099: if (pc->useAmat) {
1100: PetscObjectGetOptionsPrefix((PetscObject)mat,&mprefix);
1101: MatGetSubMatrices(mat,n_local,bjac->is,bjac->is,scall,&bjac->mat);
1102: }
1103: /* Return control to the user so that the submatrices can be modified (e.g., to apply
1104: different boundary conditions for the submatrices than for the global problem) */
1105: PCModifySubMatrices(pc,n_local,bjac->is,bjac->is,bjac->pmat,pc->modifysubmatricesP);
1107: PetscObjectGetOptionsPrefix((PetscObject)pmat,&pprefix);
1108: for (i=0; i<n_local; i++) {
1109: PetscLogObjectParent(pc,bjac->pmat[i]);
1110: PetscObjectSetOptionsPrefix((PetscObject)bjac->pmat[i],pprefix);
1111: if (pc->useAmat) {
1112: PetscLogObjectParent(pc,bjac->mat[i]);
1113: PetscObjectSetOptionsPrefix((PetscObject)bjac->mat[i],mprefix);
1114: KSPSetOperators(jac->ksp[i],bjac->mat[i],bjac->pmat[i],pc->flag);
1115: } else {
1116: KSPSetOperators(jac->ksp[i],bjac->pmat[i],bjac->pmat[i],pc->flag);
1117: }
1118: if (pc->setfromoptionscalled) {
1119: KSPSetFromOptions(jac->ksp[i]);
1120: }
1121: }
1122: return(0);
1123: }
1125: /* ---------------------------------------------------------------------------------------------*/
1126: /*
1127: These are for a single block with multiple processes;
1128: */
1131: static PetscErrorCode PCReset_BJacobi_Multiproc(PC pc)
1132: {
1133: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
1134: PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc*)jac->data;
1135: PetscErrorCode ierr;
1138: VecDestroy(&mpjac->ysub);
1139: VecDestroy(&mpjac->xsub);
1140: MatDestroy(&mpjac->submats);
1141: if (jac->ksp) {KSPReset(jac->ksp[0]);}
1142: return(0);
1143: }
1147: static PetscErrorCode PCDestroy_BJacobi_Multiproc(PC pc)
1148: {
1149: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
1150: PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc*)jac->data;
1151: PetscErrorCode ierr;
1154: PCReset_BJacobi_Multiproc(pc);
1155: KSPDestroy(&jac->ksp[0]);
1156: PetscFree(jac->ksp);
1157: PetscSubcommDestroy(&mpjac->psubcomm);
1159: PetscFree(mpjac);
1160: PetscFree(pc->data);
1161: return(0);
1162: }
1166: static PetscErrorCode PCApply_BJacobi_Multiproc(PC pc,Vec x,Vec y)
1167: {
1168: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
1169: PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc*)jac->data;
1170: PetscErrorCode ierr;
1171: PetscScalar *xarray,*yarray;
1174: /* place x's and y's local arrays into xsub and ysub */
1175: VecGetArray(x,&xarray);
1176: VecGetArray(y,&yarray);
1177: VecPlaceArray(mpjac->xsub,xarray);
1178: VecPlaceArray(mpjac->ysub,yarray);
1180: /* apply preconditioner on each matrix block */
1181: PetscLogEventBegin(PC_ApplyOnMproc,jac->ksp[0],mpjac->xsub,mpjac->ysub,0);
1182: KSPSolve(jac->ksp[0],mpjac->xsub,mpjac->ysub);
1183: PetscLogEventEnd(PC_ApplyOnMproc,jac->ksp[0],mpjac->xsub,mpjac->ysub,0);
1185: VecResetArray(mpjac->xsub);
1186: VecResetArray(mpjac->ysub);
1187: VecRestoreArray(x,&xarray);
1188: VecRestoreArray(y,&yarray);
1189: return(0);
1190: }
1192: #include <petsc-private/matimpl.h>
1195: static PetscErrorCode PCSetUp_BJacobi_Multiproc(PC pc)
1196: {
1197: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
1198: PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc*)jac->data;
1199: PetscErrorCode ierr;
1200: PetscInt m,n;
1201: MPI_Comm comm,subcomm=0;
1202: const char *prefix;
1203: PetscBool wasSetup = PETSC_TRUE;
1206: PetscObjectGetComm((PetscObject)pc,&comm);
1207: if (jac->n_local > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Only a single block in a subcommunicator is supported");
1208: jac->n_local = 1; /* currently only a single block is supported for a subcommunicator */
1209: if (!pc->setupcalled) {
1210: wasSetup = PETSC_FALSE;
1211: PetscNewLog(pc,PC_BJacobi_Multiproc,&mpjac);
1212: jac->data = (void*)mpjac;
1214: /* initialize datastructure mpjac */
1215: if (!jac->psubcomm) {
1216: /* Create default contiguous subcommunicatiors if user does not provide them */
1217: PetscSubcommCreate(comm,&jac->psubcomm);
1218: PetscSubcommSetNumber(jac->psubcomm,jac->n);
1219: PetscSubcommSetType(jac->psubcomm,PETSC_SUBCOMM_CONTIGUOUS);
1220: PetscLogObjectMemory(pc,sizeof(PetscSubcomm));
1221: }
1222: mpjac->psubcomm = jac->psubcomm;
1223: subcomm = mpjac->psubcomm->comm;
1225: /* Get matrix blocks of pmat */
1226: if (!pc->pmat->ops->getmultiprocblock) SETERRQ(PetscObjectComm((PetscObject)pc->pmat),PETSC_ERR_SUP,"No support for the requested operation");
1227: (*pc->pmat->ops->getmultiprocblock)(pc->pmat,subcomm,MAT_INITIAL_MATRIX,&mpjac->submats);
1229: /* create a new PC that processors in each subcomm have copy of */
1230: PetscMalloc(sizeof(KSP),&jac->ksp);
1231: KSPCreate(subcomm,&jac->ksp[0]);
1232: PetscObjectIncrementTabLevel((PetscObject)jac->ksp[0],(PetscObject)pc,1);
1233: PetscLogObjectParent(pc,jac->ksp[0]);
1234: KSPSetOperators(jac->ksp[0],mpjac->submats,mpjac->submats,pc->flag);
1235: KSPGetPC(jac->ksp[0],&mpjac->pc);
1237: PCGetOptionsPrefix(pc,&prefix);
1238: KSPSetOptionsPrefix(jac->ksp[0],prefix);
1239: KSPAppendOptionsPrefix(jac->ksp[0],"sub_");
1240: /*
1241: PetscMPIInt rank,subsize,subrank;
1242: MPI_Comm_rank(comm,&rank);
1243: MPI_Comm_size(subcomm,&subsize);
1244: MPI_Comm_rank(subcomm,&subrank);
1246: MatGetLocalSize(mpjac->submats,&m,NULL);
1247: MatGetSize(mpjac->submats,&n,NULL);
1248: PetscSynchronizedPrintf(comm,"[%d], sub-size %d,sub-rank %d\n",rank,subsize,subrank);
1249: PetscSynchronizedFlush(comm);
1250: */
1252: /* create dummy vectors xsub and ysub */
1253: MatGetLocalSize(mpjac->submats,&m,&n);
1254: VecCreateMPIWithArray(subcomm,1,n,PETSC_DECIDE,NULL,&mpjac->xsub);
1255: VecCreateMPIWithArray(subcomm,1,m,PETSC_DECIDE,NULL,&mpjac->ysub);
1256: PetscLogObjectParent(pc,mpjac->xsub);
1257: PetscLogObjectParent(pc,mpjac->ysub);
1259: pc->ops->reset = PCReset_BJacobi_Multiproc;
1260: pc->ops->destroy = PCDestroy_BJacobi_Multiproc;
1261: pc->ops->apply = PCApply_BJacobi_Multiproc;
1262: } else { /* pc->setupcalled */
1263: subcomm = mpjac->psubcomm->comm;
1264: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
1265: /* destroy old matrix blocks, then get new matrix blocks */
1266: if (mpjac->submats) {MatDestroy(&mpjac->submats);}
1267: (*pc->pmat->ops->getmultiprocblock)(pc->pmat,subcomm,MAT_INITIAL_MATRIX,&mpjac->submats);
1268: } else {
1269: (*pc->pmat->ops->getmultiprocblock)(pc->pmat,subcomm,MAT_REUSE_MATRIX,&mpjac->submats);
1270: }
1271: KSPSetOperators(jac->ksp[0],mpjac->submats,mpjac->submats,pc->flag);
1272: }
1274: if (!wasSetup && pc->setfromoptionscalled) {
1275: KSPSetFromOptions(jac->ksp[0]);
1276: }
1277: return(0);
1278: }