Actual source code: bjacobi.c
petsc-3.3-p7 2013-05-11
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(((PetscObject)pc)->comm,&rank);
26: MPI_Comm_size(((PetscObject)pc)->comm,&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,((PetscObject)pc)->comm);
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++) {
52: jac->l_lens[i] = bs*((M/bs)/jac->n_local + (((M/bs) % jac->n_local) > i));
53: }
54: }
55: } else if (jac->n > 0 && jac->n_local < 0) { /* global block count given */
56: /* global blocks given: determine which ones are local */
57: if (jac->g_lens) {
58: /* check if the g_lens is has valid entries */
59: for (i=0; i<jac->n; i++) {
60: if (!jac->g_lens[i]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Zero block not allowed");
61: 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");
62: }
63: if (size == 1) {
64: jac->n_local = jac->n;
65: PetscMalloc(jac->n_local*sizeof(PetscInt),&jac->l_lens);
66: PetscMemcpy(jac->l_lens,jac->g_lens,jac->n_local*sizeof(PetscInt));
67: /* check that user set these correctly */
68: sum = 0;
69: for (i=0; i<jac->n_local; i++) sum += jac->l_lens[i];
70: if (sum != M) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Global lens set incorrectly");
71: } else {
72: MatGetOwnershipRange(pc->pmat,&start,&end);
73: /* loop over blocks determing first one owned by me */
74: sum = 0;
75: for (i=0; i<jac->n+1; i++) {
76: if (sum == start) { i_start = i; goto start_1;}
77: if (i < jac->n) sum += jac->g_lens[i];
78: }
79: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Block sizes used in PCBJacobiSetTotalBlocks()\nare not compatible with parallel matrix layout");
80: start_1:
81: for (i=i_start; i<jac->n+1; i++) {
82: if (sum == end) { i_end = i; goto end_1; }
83: if (i < jac->n) sum += jac->g_lens[i];
84: }
85: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Block sizes used in PCBJacobiSetTotalBlocks()\nare not compatible with parallel matrix layout");
86: end_1:
87: jac->n_local = i_end - i_start;
88: PetscMalloc(jac->n_local*sizeof(PetscInt),&jac->l_lens);
89: PetscMemcpy(jac->l_lens,jac->g_lens+i_start,jac->n_local*sizeof(PetscInt));
90: }
91: } else { /* no global blocks given, determine then using default layout */
92: jac->n_local = jac->n/size + ((jac->n % size) > rank);
93: PetscMalloc(jac->n_local*sizeof(PetscInt),&jac->l_lens);
94: for (i=0; i<jac->n_local; i++) {
95: jac->l_lens[i] = ((M/bs)/jac->n_local + (((M/bs) % jac->n_local) > i))*bs;
96: if (!jac->l_lens[i]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Too many blocks given");
97: }
98: }
99: } else if (jac->n < 0 && jac->n_local < 0) { /* no blocks given */
100: jac->n = size;
101: jac->n_local = 1;
102: PetscMalloc(sizeof(PetscInt),&jac->l_lens);
103: jac->l_lens[0] = M;
104: } else { /* jac->n > 0 && jac->n_local > 0 */
105: if (!jac->l_lens) {
106: PetscMalloc(jac->n_local*sizeof(PetscInt),&jac->l_lens);
107: for (i=0; i<jac->n_local; i++) {
108: jac->l_lens[i] = bs*((M/bs)/jac->n_local + (((M/bs) % jac->n_local) > i));
109: }
110: }
111: }
112: if (jac->n_local < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Number of blocks is less than number of processors");
114: /* -------------------------
115: Determines mat and pmat
116: ---------------------------*/
117: PetscObjectQueryFunction((PetscObject)pc->mat,"MatGetDiagonalBlock_C",(void (**)(void))&f);
118: if (!f && size == 1) {
119: mat = pc->mat;
120: pmat = pc->pmat;
121: } else {
122: if (jac->use_true_local) {
123: /* use block from true matrix, not preconditioner matrix for local MatMult() */
124: MatGetDiagonalBlock(pc->mat,&mat);
125: /* make submatrix have same prefix as entire matrix */
126: PetscObjectGetOptionsPrefix((PetscObject)pc->mat,&mprefix);
127: PetscObjectSetOptionsPrefix((PetscObject)mat,mprefix);
128: }
129: if (pc->pmat != pc->mat || !jac->use_true_local) {
130: MatGetDiagonalBlock(pc->pmat,&pmat);
131: /* make submatrix have same prefix as entire matrix */
132: PetscObjectGetOptionsPrefix((PetscObject)pc->pmat,&pprefix);
133: PetscObjectSetOptionsPrefix((PetscObject)pmat,pprefix);
134: } else {
135: pmat = mat;
136: }
137: }
139: /* ------
140: Setup code depends on the number of blocks
141: */
142: if (jac->n_local == 1) {
143: PCSetUp_BJacobi_Singleblock(pc,mat,pmat);
144: } else {
145: PCSetUp_BJacobi_Multiblock(pc,mat,pmat);
146: }
147: return(0);
148: }
150: /* Default destroy, if it has never been setup */
153: static PetscErrorCode PCDestroy_BJacobi(PC pc)
154: {
155: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
159: PetscFree(jac->g_lens);
160: PetscFree(jac->l_lens);
161: PetscFree(pc->data);
162: return(0);
163: }
168: static PetscErrorCode PCSetFromOptions_BJacobi(PC pc)
169: {
170: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
172: PetscInt blocks;
173: PetscBool flg;
176: PetscOptionsHead("Block Jacobi options");
177: PetscOptionsInt("-pc_bjacobi_blocks","Total number of blocks","PCBJacobiSetTotalBlocks",jac->n,&blocks,&flg);
178: if (flg) {
179: PCBJacobiSetTotalBlocks(pc,blocks,PETSC_NULL);
180: }
181: flg = PETSC_FALSE;
182: PetscOptionsBool("-pc_bjacobi_truelocal","Use the true matrix, not preconditioner matrix to define matrix vector product in sub-problems","PCBJacobiSetUseTrueLocal",flg,&flg,PETSC_NULL);
183: if (flg) {
184: PCBJacobiSetUseTrueLocal(pc);
185: }
186: PetscOptionsTail();
187: return(0);
188: }
192: static PetscErrorCode PCView_BJacobi(PC pc,PetscViewer viewer)
193: {
194: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
195: PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc*)jac->data;
196: PetscErrorCode ierr;
197: PetscMPIInt rank;
198: PetscInt i;
199: PetscBool iascii,isstring;
200: PetscViewer sviewer;
203: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
204: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
205: if (iascii) {
206: if (jac->use_true_local) {
207: PetscViewerASCIIPrintf(viewer," block Jacobi: using true local matrix, number of blocks = %D\n",jac->n);
208: }
209: PetscViewerASCIIPrintf(viewer," block Jacobi: number of blocks = %D\n",jac->n);
210: MPI_Comm_rank(((PetscObject)pc)->comm,&rank);
211: if (jac->same_local_solves) {
212: PetscViewerASCIIPrintf(viewer," Local solve is same for all blocks, in the following KSP and PC objects:\n");
213: if (jac->ksp && !jac->psubcomm) {
214: PetscViewerGetSingleton(viewer,&sviewer);
215: if (!rank){
216: PetscViewerASCIIPushTab(viewer);
217: KSPView(jac->ksp[0],sviewer);
218: PetscViewerASCIIPopTab(viewer);
219: }
220: PetscViewerRestoreSingleton(viewer,&sviewer);
221: } else if (jac->psubcomm && !jac->psubcomm->color){
222: PetscViewerASCIIGetStdout(mpjac->psubcomm->comm,&sviewer);
223: PetscViewerASCIIPushTab(viewer);
224: KSPView(*(jac->ksp),sviewer);
225: PetscViewerASCIIPopTab(viewer);
226: }
227: } else {
228: PetscInt n_global;
229: MPI_Allreduce(&jac->n_local,&n_global,1,MPIU_INT,MPI_MAX,((PetscObject)pc)->comm);
230: PetscViewerASCIISynchronizedAllow(viewer,PETSC_TRUE);
231: PetscViewerASCIIPrintf(viewer," Local solve info for each block is in the following KSP and PC objects:\n");
232: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] number of local blocks = %D, first local block number = %D\n",
233: rank,jac->n_local,jac->first_local);
234: PetscViewerASCIIPushTab(viewer);
235: for (i=0; i<n_global; i++) {
236: PetscViewerGetSingleton(viewer,&sviewer);
237: if (i < jac->n_local) {
238: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] local block number %D\n",rank,i);
239: KSPView(jac->ksp[i],sviewer);
240: PetscViewerASCIISynchronizedPrintf(viewer,"- - - - - - - - - - - - - - - - - -\n");
241: }
242: PetscViewerRestoreSingleton(viewer,&sviewer);
243: }
244: PetscViewerASCIIPopTab(viewer);
245: PetscViewerFlush(viewer);
246: PetscViewerASCIISynchronizedAllow(viewer,PETSC_FALSE);
247: }
248: } else if (isstring) {
249: PetscViewerStringSPrintf(viewer," blks=%D",jac->n);
250: PetscViewerGetSingleton(viewer,&sviewer);
251: if (jac->ksp) {KSPView(jac->ksp[0],sviewer);}
252: PetscViewerRestoreSingleton(viewer,&sviewer);
253: } else {
254: SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_SUP,"Viewer type %s not supported for block Jacobi",((PetscObject)viewer)->type_name);
255: }
256: return(0);
257: }
259: /* -------------------------------------------------------------------------------------*/
261: EXTERN_C_BEGIN
264: PetscErrorCode PCBJacobiSetUseTrueLocal_BJacobi(PC pc)
265: {
266: PC_BJacobi *jac;
269: jac = (PC_BJacobi*)pc->data;
270: jac->use_true_local = PETSC_TRUE;
271: return(0);
272: }
273: EXTERN_C_END
275: EXTERN_C_BEGIN
278: PetscErrorCode PCBJacobiGetSubKSP_BJacobi(PC pc,PetscInt *n_local,PetscInt *first_local,KSP **ksp)
279: {
280: PC_BJacobi *jac = (PC_BJacobi*)pc->data;;
283: if (!pc->setupcalled) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"Must call KSPSetUp() or PCSetUp() first");
285: if (n_local) *n_local = jac->n_local;
286: if (first_local) *first_local = jac->first_local;
287: *ksp = jac->ksp;
288: jac->same_local_solves = PETSC_FALSE; /* Assume that local solves are now different;
289: not necessarily true though! This flag is
290: used only for PCView_BJacobi() */
291: return(0);
292: }
293: EXTERN_C_END
295: EXTERN_C_BEGIN
298: PetscErrorCode PCBJacobiSetTotalBlocks_BJacobi(PC pc,PetscInt blocks,PetscInt *lens)
299: {
300: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
305: if (pc->setupcalled > 0 && jac->n!=blocks) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ORDER,"Cannot alter number of blocks after PCSetUp()/KSPSetUp() has been called");
306: jac->n = blocks;
307: if (!lens) {
308: jac->g_lens = 0;
309: } else {
310: PetscMalloc(blocks*sizeof(PetscInt),&jac->g_lens);
311: PetscLogObjectMemory(pc,blocks*sizeof(PetscInt));
312: PetscMemcpy(jac->g_lens,lens,blocks*sizeof(PetscInt));
313: }
314: return(0);
315: }
316: EXTERN_C_END
318: EXTERN_C_BEGIN
321: PetscErrorCode PCBJacobiGetTotalBlocks_BJacobi(PC pc, PetscInt *blocks, const PetscInt *lens[])
322: {
323: PC_BJacobi *jac = (PC_BJacobi*) pc->data;
326: *blocks = jac->n;
327: if (lens) *lens = jac->g_lens;
328: return(0);
329: }
330: EXTERN_C_END
332: EXTERN_C_BEGIN
335: PetscErrorCode PCBJacobiSetLocalBlocks_BJacobi(PC pc,PetscInt blocks,const PetscInt lens[])
336: {
337: PC_BJacobi *jac;
341: jac = (PC_BJacobi*)pc->data;
343: jac->n_local = blocks;
344: if (!lens) {
345: jac->l_lens = 0;
346: } else {
347: PetscMalloc(blocks*sizeof(PetscInt),&jac->l_lens);
348: PetscLogObjectMemory(pc,blocks*sizeof(PetscInt));
349: PetscMemcpy(jac->l_lens,lens,blocks*sizeof(PetscInt));
350: }
351: return(0);
352: }
353: EXTERN_C_END
355: EXTERN_C_BEGIN
358: PetscErrorCode PCBJacobiGetLocalBlocks_BJacobi(PC pc, PetscInt *blocks, const PetscInt *lens[])
359: {
360: PC_BJacobi *jac = (PC_BJacobi*) pc->data;
363: *blocks = jac->n_local;
364: if (lens) *lens = jac->l_lens;
365: return(0);
366: }
367: EXTERN_C_END
369: /* -------------------------------------------------------------------------------------*/
373: /*@
374: PCBJacobiSetUseTrueLocal - Sets a flag to indicate that the block
375: problem is associated with the linear system matrix instead of the
376: default (where it is associated with the preconditioning matrix).
377: That is, if the local system is solved iteratively then it iterates
378: on the block from the matrix using the block from the preconditioner
379: as the preconditioner for the local block.
381: Logically Collective on PC
383: Input Parameters:
384: . pc - the preconditioner context
386: Options Database Key:
387: . -pc_bjacobi_truelocal - Activates PCBJacobiSetUseTrueLocal()
389: Notes:
390: For the common case in which the preconditioning and linear
391: system matrices are identical, this routine is unnecessary.
393: Level: intermediate
395: .keywords: block, Jacobi, set, true, local, flag
397: .seealso: PCSetOperators(), PCBJacobiSetLocalBlocks()
398: @*/
399: PetscErrorCode PCBJacobiSetUseTrueLocal(PC pc)
400: {
405: PetscTryMethod(pc,"PCBJacobiSetUseTrueLocal_C",(PC),(pc));
406: return(0);
407: }
411: /*@C
412: PCBJacobiGetSubKSP - Gets the local KSP contexts for all blocks on
413: this processor.
414:
415: Note Collective
417: Input Parameter:
418: . pc - the preconditioner context
420: Output Parameters:
421: + n_local - the number of blocks on this processor, or PETSC_NULL
422: . first_local - the global number of the first block on this processor, or PETSC_NULL
423: - ksp - the array of KSP contexts
425: Notes:
426: After PCBJacobiGetSubKSP() the array of KSP contexts is not to be freed.
427:
428: Currently for some matrix implementations only 1 block per processor
429: is supported.
430:
431: You must call KSPSetUp() or PCSetUp() before calling PCBJacobiGetSubKSP().
433: Level: advanced
435: .keywords: block, Jacobi, get, sub, KSP, context
437: .seealso: PCBJacobiGetSubKSP()
438: @*/
439: PetscErrorCode PCBJacobiGetSubKSP(PC pc,PetscInt *n_local,PetscInt *first_local,KSP *ksp[])
440: {
445: PetscUseMethod(pc,"PCBJacobiGetSubKSP_C",(PC,PetscInt *,PetscInt *,KSP **),(pc,n_local,first_local,ksp));
446: return(0);
447: }
451: /*@
452: PCBJacobiSetTotalBlocks - Sets the global number of blocks for the block
453: Jacobi preconditioner.
455: Collective on PC
457: Input Parameters:
458: + pc - the preconditioner context
459: . blocks - the number of blocks
460: - lens - [optional] integer array containing the size of each block
462: Options Database Key:
463: . -pc_bjacobi_blocks <blocks> - Sets the number of global blocks
465: Notes:
466: Currently only a limited number of blocking configurations are supported.
467: All processors sharing the PC must call this routine with the same data.
469: Level: intermediate
471: .keywords: set, number, Jacobi, global, total, blocks
473: .seealso: PCBJacobiSetUseTrueLocal(), PCBJacobiSetLocalBlocks()
474: @*/
475: PetscErrorCode PCBJacobiSetTotalBlocks(PC pc,PetscInt blocks,const PetscInt lens[])
476: {
481: if (blocks <= 0) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Must have positive blocks");
482: PetscTryMethod(pc,"PCBJacobiSetTotalBlocks_C",(PC,PetscInt,const PetscInt[]),(pc,blocks,lens));
483: return(0);
484: }
488: /*@C
489: PCBJacobiGetTotalBlocks - Gets the global number of blocks for the block
490: Jacobi preconditioner.
492: Not Collective
494: Input Parameter:
495: . pc - the preconditioner context
497: Output parameters:
498: + blocks - the number of blocks
499: - lens - integer array containing the size of each block
501: Level: intermediate
503: .keywords: get, number, Jacobi, global, total, blocks
505: .seealso: PCBJacobiSetUseTrueLocal(), PCBJacobiGetLocalBlocks()
506: @*/
507: PetscErrorCode PCBJacobiGetTotalBlocks(PC pc, PetscInt *blocks, const PetscInt *lens[])
508: {
514: PetscUseMethod(pc,"PCBJacobiGetTotalBlocks_C",(PC,PetscInt*, const PetscInt *[]),(pc,blocks,lens));
515: return(0);
516: }
517:
520: /*@
521: PCBJacobiSetLocalBlocks - Sets the local number of blocks for the block
522: Jacobi preconditioner.
524: Not Collective
526: Input Parameters:
527: + pc - the preconditioner context
528: . blocks - the number of blocks
529: - lens - [optional] integer array containing size of each block
531: Note:
532: Currently only a limited number of blocking configurations are supported.
534: Level: intermediate
536: .keywords: PC, set, number, Jacobi, local, blocks
538: .seealso: PCBJacobiSetUseTrueLocal(), PCBJacobiSetTotalBlocks()
539: @*/
540: PetscErrorCode PCBJacobiSetLocalBlocks(PC pc,PetscInt blocks,const PetscInt lens[])
541: {
546: if (blocks < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must have nonegative blocks");
547: PetscTryMethod(pc,"PCBJacobiSetLocalBlocks_C",(PC,PetscInt,const PetscInt []),(pc,blocks,lens));
548: return(0);
549: }
550:
553: /*@C
554: PCBJacobiGetLocalBlocks - Gets the local number of blocks for the block
555: Jacobi preconditioner.
557: Not Collective
559: Input Parameters:
560: + pc - the preconditioner context
561: . blocks - the number of blocks
562: - lens - [optional] integer array containing size of each block
564: Note:
565: Currently only a limited number of blocking configurations are supported.
567: Level: intermediate
569: .keywords: PC, get, number, Jacobi, local, blocks
571: .seealso: PCBJacobiSetUseTrueLocal(), PCBJacobiGetTotalBlocks()
572: @*/
573: PetscErrorCode PCBJacobiGetLocalBlocks(PC pc, PetscInt *blocks, const PetscInt *lens[])
574: {
580: PetscUseMethod(pc,"PCBJacobiGetLocalBlocks_C",(PC,PetscInt*, const PetscInt *[]),(pc,blocks,lens));
581: return(0);
582: }
584: /* -----------------------------------------------------------------------------------*/
586: /*MC
587: PCBJACOBI - Use block Jacobi preconditioning, each block is (approximately) solved with
588: its own KSP object.
590: Options Database Keys:
591: . -pc_bjacobi_truelocal - Activates PCBJacobiSetUseTrueLocal()
593: Notes: Each processor can have one or more blocks, but a block cannot be shared by more
594: than one processor. Defaults to one block per processor.
596: To set options on the solvers for each block append -sub_ to all the KSP, KSP, and PC
597: options database keys. For example, -sub_pc_type ilu -sub_pc_factor_levels 1 -sub_ksp_type preonly
598:
599: To set the options on the solvers separate for each block call PCBJacobiGetSubKSP()
600: and set the options directly on the resulting KSP object (you can access its PC
601: KSPGetPC())
603: Level: beginner
605: Concepts: block Jacobi
607: Developer Notes: This preconditioner does not currently work with CUDA/CUSP for a couple of reasons.
608: (1) It creates seq vectors as work vectors that should be cusp
609: (2) The use of VecPlaceArray() is not handled properly by CUSP (that is it will not know where
610: the ownership of the vector is so may use wrong values) even if it did know the ownership
611: it may induce extra copy ups and downs. Satish suggests a VecTransplantArray() to handle two
612: vectors sharing the same pointer and handling the CUSP side as well instead of VecGetArray()/VecPlaceArray().
615: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
616: PCASM, PCBJacobiSetUseTrueLocal(), PCBJacobiGetSubKSP(), PCBJacobiSetTotalBlocks(),
617: PCBJacobiSetLocalBlocks(), PCSetModifySubmatrices()
618: M*/
620: EXTERN_C_BEGIN
623: PetscErrorCode PCCreate_BJacobi(PC pc)
624: {
626: PetscMPIInt rank;
627: PC_BJacobi *jac;
630: PetscNewLog(pc,PC_BJacobi,&jac);
631: MPI_Comm_rank(((PetscObject)pc)->comm,&rank);
632: pc->ops->apply = 0;
633: pc->ops->applytranspose = 0;
634: pc->ops->setup = PCSetUp_BJacobi;
635: pc->ops->destroy = PCDestroy_BJacobi;
636: pc->ops->setfromoptions = PCSetFromOptions_BJacobi;
637: pc->ops->view = PCView_BJacobi;
638: pc->ops->applyrichardson = 0;
640: pc->data = (void*)jac;
641: jac->n = -1;
642: jac->n_local = -1;
643: jac->first_local = rank;
644: jac->ksp = 0;
645: jac->use_true_local = PETSC_FALSE;
646: jac->same_local_solves = PETSC_TRUE;
647: jac->g_lens = 0;
648: jac->l_lens = 0;
649: jac->psubcomm = 0;
651: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCBJacobiSetUseTrueLocal_C",
652: "PCBJacobiSetUseTrueLocal_BJacobi",
653: PCBJacobiSetUseTrueLocal_BJacobi);
654: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCBJacobiGetSubKSP_C","PCBJacobiGetSubKSP_BJacobi",
655: PCBJacobiGetSubKSP_BJacobi);
656: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCBJacobiSetTotalBlocks_C","PCBJacobiSetTotalBlocks_BJacobi",
657: PCBJacobiSetTotalBlocks_BJacobi);
658: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCBJacobiGetTotalBlocks_C","PCBJacobiGetTotalBlocks_BJacobi",
659: PCBJacobiGetTotalBlocks_BJacobi);
660: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCBJacobiSetLocalBlocks_C","PCBJacobiSetLocalBlocks_BJacobi",
661: PCBJacobiSetLocalBlocks_BJacobi);
662: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCBJacobiGetLocalBlocks_C","PCBJacobiGetLocalBlocks_BJacobi",
663: PCBJacobiGetLocalBlocks_BJacobi);
665: return(0);
666: }
667: EXTERN_C_END
669: /* --------------------------------------------------------------------------------------------*/
670: /*
671: These are for a single block per processor; works for AIJ, BAIJ; Seq and MPI
672: */
675: PetscErrorCode PCReset_BJacobi_Singleblock(PC pc)
676: {
677: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
678: PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock*)jac->data;
679: PetscErrorCode ierr;
682: KSPReset(jac->ksp[0]);
683: VecDestroy(&bjac->x);
684: VecDestroy(&bjac->y);
685: return(0);
686: }
690: PetscErrorCode PCDestroy_BJacobi_Singleblock(PC pc)
691: {
692: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
693: PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock*)jac->data;
694: PetscErrorCode ierr;
697: PCReset_BJacobi_Singleblock(pc);
698: KSPDestroy(&jac->ksp[0]);
699: PetscFree(jac->ksp);
700: PetscFree(jac->l_lens);
701: PetscFree(jac->g_lens);
702: PetscFree(bjac);
703: PetscFree(pc->data);
704: return(0);
705: }
709: PetscErrorCode PCSetUpOnBlocks_BJacobi_Singleblock(PC pc)
710: {
712: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
715: KSPSetUp(jac->ksp[0]);
716: return(0);
717: }
721: PetscErrorCode PCApply_BJacobi_Singleblock(PC pc,Vec x,Vec y)
722: {
723: PetscErrorCode ierr;
724: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
725: PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock*)jac->data;
726: PetscScalar *x_array,*y_array;
729: /*
730: The VecPlaceArray() is to avoid having to copy the
731: y vector into the bjac->x vector. The reason for
732: the bjac->x vector is that we need a sequential vector
733: for the sequential solve.
734: */
735: VecGetArray(x,&x_array);
736: VecGetArray(y,&y_array);
737: VecPlaceArray(bjac->x,x_array);
738: VecPlaceArray(bjac->y,y_array);
739: KSPSolve(jac->ksp[0],bjac->x,bjac->y);
740: VecResetArray(bjac->x);
741: VecResetArray(bjac->y);
742: VecRestoreArray(x,&x_array);
743: VecRestoreArray(y,&y_array);
744: return(0);
745: }
749: PetscErrorCode PCApplySymmetricLeft_BJacobi_Singleblock(PC pc,Vec x,Vec y)
750: {
751: PetscErrorCode ierr;
752: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
753: PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock*)jac->data;
754: PetscScalar *x_array,*y_array;
755: PC subpc;
758: /*
759: The VecPlaceArray() is to avoid having to copy the
760: y vector into the bjac->x vector. The reason for
761: the bjac->x vector is that we need a sequential vector
762: for the sequential solve.
763: */
764: VecGetArray(x,&x_array);
765: VecGetArray(y,&y_array);
766: VecPlaceArray(bjac->x,x_array);
767: VecPlaceArray(bjac->y,y_array);
769: /* apply the symmetric left portion of the inner PC operator */
770: /* note this by-passes the inner KSP and its options completely */
772: KSPGetPC(jac->ksp[0],&subpc);
773: PCApplySymmetricLeft(subpc,bjac->x,bjac->y);
774: VecResetArray(bjac->x);
775: VecResetArray(bjac->y);
777: VecRestoreArray(x,&x_array);
778: VecRestoreArray(y,&y_array);
779: return(0);
780: }
784: PetscErrorCode PCApplySymmetricRight_BJacobi_Singleblock(PC pc,Vec x,Vec y)
785: {
786: PetscErrorCode ierr;
787: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
788: PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock*)jac->data;
789: PetscScalar *x_array,*y_array;
790: PC subpc;
793: /*
794: The VecPlaceArray() is to avoid having to copy the
795: y vector into the bjac->x vector. The reason for
796: the bjac->x vector is that we need a sequential vector
797: for the sequential solve.
798: */
799: VecGetArray(x,&x_array);
800: VecGetArray(y,&y_array);
801: VecPlaceArray(bjac->x,x_array);
802: VecPlaceArray(bjac->y,y_array);
804: /* apply the symmetric right portion of the inner PC operator */
805: /* note this by-passes the inner KSP and its options completely */
807: KSPGetPC(jac->ksp[0],&subpc);
808: PCApplySymmetricRight(subpc,bjac->x,bjac->y);
810: VecRestoreArray(x,&x_array);
811: VecRestoreArray(y,&y_array);
812: return(0);
813: }
817: PetscErrorCode PCApplyTranspose_BJacobi_Singleblock(PC pc,Vec x,Vec y)
818: {
819: PetscErrorCode ierr;
820: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
821: PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock*)jac->data;
822: PetscScalar *x_array,*y_array;
825: /*
826: The VecPlaceArray() is to avoid having to copy the
827: y vector into the bjac->x vector. The reason for
828: the bjac->x vector is that we need a sequential vector
829: for the sequential solve.
830: */
831: VecGetArray(x,&x_array);
832: VecGetArray(y,&y_array);
833: VecPlaceArray(bjac->x,x_array);
834: VecPlaceArray(bjac->y,y_array);
835: KSPSolveTranspose(jac->ksp[0],bjac->x,bjac->y);
836: VecResetArray(bjac->x);
837: VecResetArray(bjac->y);
838: VecRestoreArray(x,&x_array);
839: VecRestoreArray(y,&y_array);
840: return(0);
841: }
845: static PetscErrorCode PCSetUp_BJacobi_Singleblock(PC pc,Mat mat,Mat pmat)
846: {
847: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
848: PetscErrorCode ierr;
849: PetscInt m;
850: KSP ksp;
851: PC_BJacobi_Singleblock *bjac;
852: PetscBool wasSetup = PETSC_TRUE;
856: if (!pc->setupcalled) {
857: const char *prefix;
859: if (!jac->ksp) {
860: wasSetup = PETSC_FALSE;
861: KSPCreate(PETSC_COMM_SELF,&ksp);
862: PetscObjectIncrementTabLevel((PetscObject)ksp,(PetscObject)pc,1);
863: PetscLogObjectParent(pc,ksp);
864: KSPSetType(ksp,KSPPREONLY);
865: PCGetOptionsPrefix(pc,&prefix);
866: KSPSetOptionsPrefix(ksp,prefix);
867: KSPAppendOptionsPrefix(ksp,"sub_");
869: pc->ops->reset = PCReset_BJacobi_Singleblock;
870: pc->ops->destroy = PCDestroy_BJacobi_Singleblock;
871: pc->ops->apply = PCApply_BJacobi_Singleblock;
872: pc->ops->applysymmetricleft = PCApplySymmetricLeft_BJacobi_Singleblock;
873: pc->ops->applysymmetricright = PCApplySymmetricRight_BJacobi_Singleblock;
874: pc->ops->applytranspose = PCApplyTranspose_BJacobi_Singleblock;
875: pc->ops->setuponblocks = PCSetUpOnBlocks_BJacobi_Singleblock;
877: PetscMalloc(sizeof(KSP),&jac->ksp);
878: jac->ksp[0] = ksp;
880: PetscNewLog(pc,PC_BJacobi_Singleblock,&bjac);
881: jac->data = (void*)bjac;
882: } else {
883: ksp = jac->ksp[0];
884: bjac = (PC_BJacobi_Singleblock *)jac->data;
885: }
887: /*
888: The reason we need to generate these vectors is to serve
889: as the right-hand side and solution vector for the solve on the
890: block. We do not need to allocate space for the vectors since
891: that is provided via VecPlaceArray() just before the call to
892: KSPSolve() on the block.
893: */
894: MatGetSize(pmat,&m,&m);
895: VecCreateSeqWithArray(PETSC_COMM_SELF,1,m,PETSC_NULL,&bjac->x);
896: VecCreateSeqWithArray(PETSC_COMM_SELF,1,m,PETSC_NULL,&bjac->y);
897: PetscLogObjectParent(pc,bjac->x);
898: PetscLogObjectParent(pc,bjac->y);
899: } else {
900: ksp = jac->ksp[0];
901: bjac = (PC_BJacobi_Singleblock *)jac->data;
902: }
903: if (jac->use_true_local) {
904: KSPSetOperators(ksp,mat,pmat,pc->flag);
905: } else {
906: KSPSetOperators(ksp,pmat,pmat,pc->flag);
907: }
908: if (!wasSetup && pc->setfromoptionscalled) {
909: KSPSetFromOptions(ksp);
910: }
911: return(0);
912: }
914: /* ---------------------------------------------------------------------------------------------*/
917: PetscErrorCode PCReset_BJacobi_Multiblock(PC pc)
918: {
919: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
920: PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock*)jac->data;
921: PetscErrorCode ierr;
922: PetscInt i;
925: if (bjac && bjac->pmat) {
926: MatDestroyMatrices(jac->n_local,&bjac->pmat);
927: if (jac->use_true_local) {
928: MatDestroyMatrices(jac->n_local,&bjac->mat);
929: }
930: }
932: for (i=0; i<jac->n_local; i++) {
933: KSPReset(jac->ksp[i]);
934: if (bjac && bjac->x) {
935: VecDestroy(&bjac->x[i]);
936: VecDestroy(&bjac->y[i]);
937: ISDestroy(&bjac->is[i]);
938: }
939: }
940: PetscFree(jac->l_lens);
941: PetscFree(jac->g_lens);
942: return(0);
943: }
947: PetscErrorCode PCDestroy_BJacobi_Multiblock(PC pc)
948: {
949: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
950: PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock*)jac->data;
951: PetscErrorCode ierr;
952: PetscInt i;
955: PCReset_BJacobi_Multiblock(pc);
956: if (bjac) {
957: PetscFree2(bjac->x,bjac->y);
958: PetscFree(bjac->starts);
959: PetscFree(bjac->is);
960: }
961: PetscFree(jac->data);
962: for (i=0; i<jac->n_local; i++) {
963: KSPDestroy(&jac->ksp[i]);
964: }
965: PetscFree(jac->ksp);
966: PetscFree(pc->data);
967: return(0);
968: }
972: PetscErrorCode PCSetUpOnBlocks_BJacobi_Multiblock(PC pc)
973: {
974: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
976: PetscInt i,n_local = jac->n_local;
979: for (i=0; i<n_local; i++) {
980: KSPSetUp(jac->ksp[i]);
981: }
982: return(0);
983: }
985: /*
986: Preconditioner for block Jacobi
987: */
990: PetscErrorCode PCApply_BJacobi_Multiblock(PC pc,Vec x,Vec y)
991: {
992: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
993: PetscErrorCode ierr;
994: PetscInt i,n_local = jac->n_local;
995: PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock*)jac->data;
996: PetscScalar *xin,*yin;
999: VecGetArray(x,&xin);
1000: VecGetArray(y,&yin);
1001: for (i=0; i<n_local; i++) {
1002: /*
1003: To avoid copying the subvector from x into a workspace we instead
1004: make the workspace vector array point to the subpart of the array of
1005: the global vector.
1006: */
1007: VecPlaceArray(bjac->x[i],xin+bjac->starts[i]);
1008: VecPlaceArray(bjac->y[i],yin+bjac->starts[i]);
1010: PetscLogEventBegin(PC_ApplyOnBlocks,jac->ksp[i],bjac->x[i],bjac->y[i],0);
1011: KSPSolve(jac->ksp[i],bjac->x[i],bjac->y[i]);
1012: PetscLogEventEnd(PC_ApplyOnBlocks,jac->ksp[i],bjac->x[i],bjac->y[i],0);
1014: VecResetArray(bjac->x[i]);
1015: VecResetArray(bjac->y[i]);
1016: }
1017: VecRestoreArray(x,&xin);
1018: VecRestoreArray(y,&yin);
1019: return(0);
1020: }
1022: /*
1023: Preconditioner for block Jacobi
1024: */
1027: PetscErrorCode PCApplyTranspose_BJacobi_Multiblock(PC pc,Vec x,Vec y)
1028: {
1029: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
1030: PetscErrorCode ierr;
1031: PetscInt i,n_local = jac->n_local;
1032: PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock*)jac->data;
1033: PetscScalar *xin,*yin;
1036: VecGetArray(x,&xin);
1037: VecGetArray(y,&yin);
1038: for (i=0; i<n_local; i++) {
1039: /*
1040: To avoid copying the subvector from x into a workspace we instead
1041: make the workspace vector array point to the subpart of the array of
1042: the global vector.
1043: */
1044: VecPlaceArray(bjac->x[i],xin+bjac->starts[i]);
1045: VecPlaceArray(bjac->y[i],yin+bjac->starts[i]);
1047: PetscLogEventBegin(PC_ApplyTransposeOnBlocks,jac->ksp[i],bjac->x[i],bjac->y[i],0);
1048: KSPSolveTranspose(jac->ksp[i],bjac->x[i],bjac->y[i]);
1049: PetscLogEventEnd(PC_ApplyTransposeOnBlocks,jac->ksp[i],bjac->x[i],bjac->y[i],0);
1051: VecResetArray(bjac->x[i]);
1052: VecResetArray(bjac->y[i]);
1053: }
1054: VecRestoreArray(x,&xin);
1055: VecRestoreArray(y,&yin);
1056: return(0);
1057: }
1061: static PetscErrorCode PCSetUp_BJacobi_Multiblock(PC pc,Mat mat,Mat pmat)
1062: {
1063: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
1064: PetscErrorCode ierr;
1065: PetscInt m,n_local,N,M,start,i;
1066: const char *prefix,*pprefix,*mprefix;
1067: KSP ksp;
1068: Vec x,y;
1069: PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock*)jac->data;
1070: PC subpc;
1071: IS is;
1072: MatReuse scall;
1075: MatGetLocalSize(pc->pmat,&M,&N);
1077: n_local = jac->n_local;
1079: if (jac->use_true_local) {
1080: PetscBool same;
1081: PetscObjectTypeCompare((PetscObject)mat,((PetscObject)pmat)->type_name,&same);
1082: if (!same) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_INCOMP,"Matrices not of same type");
1083: }
1085: if (!pc->setupcalled) {
1086: scall = MAT_INITIAL_MATRIX;
1088: if (!jac->ksp) {
1089: pc->ops->reset = PCReset_BJacobi_Multiblock;
1090: pc->ops->destroy = PCDestroy_BJacobi_Multiblock;
1091: pc->ops->apply = PCApply_BJacobi_Multiblock;
1092: pc->ops->applytranspose= PCApplyTranspose_BJacobi_Multiblock;
1093: pc->ops->setuponblocks = PCSetUpOnBlocks_BJacobi_Multiblock;
1095: PetscNewLog(pc,PC_BJacobi_Multiblock,&bjac);
1096: PetscMalloc(n_local*sizeof(KSP),&jac->ksp);
1097: PetscLogObjectMemory(pc,sizeof(n_local*sizeof(KSP)));
1098: PetscMalloc2(n_local,Vec,&bjac->x,n_local,Vec,&bjac->y);
1099: PetscMalloc(n_local*sizeof(PetscScalar),&bjac->starts);
1100: PetscLogObjectMemory(pc,sizeof(n_local*sizeof(PetscScalar)));
1101:
1102: jac->data = (void*)bjac;
1103: PetscMalloc(n_local*sizeof(IS),&bjac->is);
1104: PetscLogObjectMemory(pc,sizeof(n_local*sizeof(IS)));
1106: for (i=0; i<n_local; i++) {
1107: KSPCreate(PETSC_COMM_SELF,&ksp);
1108: PetscObjectIncrementTabLevel((PetscObject)ksp,(PetscObject)pc,1);
1109: PetscLogObjectParent(pc,ksp);
1110: KSPSetType(ksp,KSPPREONLY);
1111: KSPGetPC(ksp,&subpc);
1112: PCGetOptionsPrefix(pc,&prefix);
1113: KSPSetOptionsPrefix(ksp,prefix);
1114: KSPAppendOptionsPrefix(ksp,"sub_");
1115: jac->ksp[i] = ksp;
1116: }
1117: } else {
1118: bjac = (PC_BJacobi_Multiblock*)jac->data;
1119: }
1121: start = 0;
1122: for (i=0; i<n_local; i++) {
1123: m = jac->l_lens[i];
1124: /*
1125: The reason we need to generate these vectors is to serve
1126: as the right-hand side and solution vector for the solve on the
1127: block. We do not need to allocate space for the vectors since
1128: that is provided via VecPlaceArray() just before the call to
1129: KSPSolve() on the block.
1131: */
1132: VecCreateSeq(PETSC_COMM_SELF,m,&x);
1133: VecCreateSeqWithArray(PETSC_COMM_SELF,1,m,PETSC_NULL,&y);
1134: PetscLogObjectParent(pc,x);
1135: PetscLogObjectParent(pc,y);
1136: bjac->x[i] = x;
1137: bjac->y[i] = y;
1138: bjac->starts[i] = start;
1140: ISCreateStride(PETSC_COMM_SELF,m,start,1,&is);
1141: bjac->is[i] = is;
1142: PetscLogObjectParent(pc,is);
1144: start += m;
1145: }
1146: } else {
1147: bjac = (PC_BJacobi_Multiblock*)jac->data;
1148: /*
1149: Destroy the blocks from the previous iteration
1150: */
1151: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
1152: MatDestroyMatrices(n_local,&bjac->pmat);
1153: if (jac->use_true_local) {
1154: MatDestroyMatrices(n_local,&bjac->mat);
1155: }
1156: scall = MAT_INITIAL_MATRIX;
1157: } else {
1158: scall = MAT_REUSE_MATRIX;
1159: }
1160: }
1162: MatGetSubMatrices(pmat,n_local,bjac->is,bjac->is,scall,&bjac->pmat);
1163: if (jac->use_true_local) {
1164: PetscObjectGetOptionsPrefix((PetscObject)mat,&mprefix);
1165: MatGetSubMatrices(mat,n_local,bjac->is,bjac->is,scall,&bjac->mat);
1166: }
1167: /* Return control to the user so that the submatrices can be modified (e.g., to apply
1168: different boundary conditions for the submatrices than for the global problem) */
1169: PCModifySubMatrices(pc,n_local,bjac->is,bjac->is,bjac->pmat,pc->modifysubmatricesP);
1171: PetscObjectGetOptionsPrefix((PetscObject)pmat,&pprefix);
1172: for (i=0; i<n_local; i++) {
1173: PetscLogObjectParent(pc,bjac->pmat[i]);
1174: PetscObjectSetOptionsPrefix((PetscObject)bjac->pmat[i],pprefix);
1175: if (jac->use_true_local) {
1176: PetscLogObjectParent(pc,bjac->mat[i]);
1177: PetscObjectSetOptionsPrefix((PetscObject)bjac->mat[i],mprefix);
1178: KSPSetOperators(jac->ksp[i],bjac->mat[i],bjac->pmat[i],pc->flag);
1179: } else {
1180: KSPSetOperators(jac->ksp[i],bjac->pmat[i],bjac->pmat[i],pc->flag);
1181: }
1182: if(pc->setfromoptionscalled){
1183: KSPSetFromOptions(jac->ksp[i]);
1184: }
1185: }
1186: return(0);
1187: }
1189: /* ---------------------------------------------------------------------------------------------*/
1190: /*
1191: These are for a single block with multiple processes;
1192: */
1195: static PetscErrorCode PCReset_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;
1202: VecDestroy(&mpjac->ysub);
1203: VecDestroy(&mpjac->xsub);
1204: MatDestroy(&mpjac->submats);
1205: if (jac->ksp){KSPReset(jac->ksp[0]);}
1206: return(0);
1207: }
1211: static PetscErrorCode PCDestroy_BJacobi_Multiproc(PC pc)
1212: {
1213: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
1214: PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc*)jac->data;
1215: PetscErrorCode ierr;
1218: PCReset_BJacobi_Multiproc(pc);
1219: KSPDestroy(&jac->ksp[0]);
1220: PetscFree(jac->ksp);
1221: PetscSubcommDestroy(&mpjac->psubcomm);
1223: PetscFree(mpjac);
1224: PetscFree(pc->data);
1225: return(0);
1226: }
1230: static PetscErrorCode PCApply_BJacobi_Multiproc(PC pc,Vec x,Vec y)
1231: {
1232: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
1233: PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc*)jac->data;
1234: PetscErrorCode ierr;
1235: PetscScalar *xarray,*yarray;
1238: /* place x's and y's local arrays into xsub and ysub */
1239: VecGetArray(x,&xarray);
1240: VecGetArray(y,&yarray);
1241: VecPlaceArray(mpjac->xsub,xarray);
1242: VecPlaceArray(mpjac->ysub,yarray);
1244: /* apply preconditioner on each matrix block */
1245: PetscLogEventBegin(PC_ApplyOnMproc,jac->ksp[0],mpjac->xsub,mpjac->ysub,0);
1246: KSPSolve(jac->ksp[0],mpjac->xsub,mpjac->ysub);
1247: PetscLogEventEnd(PC_ApplyOnMproc,jac->ksp[0],mpjac->xsub,mpjac->ysub,0);
1249: VecResetArray(mpjac->xsub);
1250: VecResetArray(mpjac->ysub);
1251: VecRestoreArray(x,&xarray);
1252: VecRestoreArray(y,&yarray);
1253: return(0);
1254: }
1256: extern PetscErrorCode MatGetMultiProcBlock_MPIAIJ(Mat,MPI_Comm,MatReuse,Mat*);
1259: static PetscErrorCode PCSetUp_BJacobi_Multiproc(PC pc)
1260: {
1261: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
1262: PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc*)jac->data;
1263: PetscErrorCode ierr;
1264: PetscInt m,n;
1265: MPI_Comm comm = ((PetscObject)pc)->comm,subcomm=0;
1266: const char *prefix;
1267: PetscBool wasSetup = PETSC_TRUE;
1270: if (jac->n_local > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Only a single block in a subcommunicator is supported");
1271: jac->n_local = 1; /* currently only a single block is supported for a subcommunicator */
1272: if (!pc->setupcalled) {
1273: wasSetup = PETSC_FALSE;
1274: PetscNewLog(pc,PC_BJacobi_Multiproc,&mpjac);
1275: jac->data = (void*)mpjac;
1277: /* initialize datastructure mpjac */
1278: if (!jac->psubcomm) {
1279: /* Create default contiguous subcommunicatiors if user does not provide them */
1280: PetscSubcommCreate(comm,&jac->psubcomm);
1281: PetscSubcommSetNumber(jac->psubcomm,jac->n);
1282: PetscSubcommSetType(jac->psubcomm,PETSC_SUBCOMM_CONTIGUOUS);
1283: PetscLogObjectMemory(pc,sizeof(PetscSubcomm));
1284: }
1285: mpjac->psubcomm = jac->psubcomm;
1286: subcomm = mpjac->psubcomm->comm;
1288: /* Get matrix blocks of pmat */
1289: MatGetMultiProcBlock_MPIAIJ(pc->pmat,subcomm,MAT_INITIAL_MATRIX,&mpjac->submats);
1291: /* create a new PC that processors in each subcomm have copy of */
1292: PetscMalloc(sizeof(KSP),&jac->ksp);
1293: KSPCreate(subcomm,&jac->ksp[0]);
1294: PetscObjectIncrementTabLevel((PetscObject)jac->ksp[0],(PetscObject)pc,1);
1295: PetscLogObjectParent(pc,jac->ksp[0]);
1296: KSPSetOperators(jac->ksp[0],mpjac->submats,mpjac->submats,pc->flag);
1297: KSPGetPC(jac->ksp[0],&mpjac->pc);
1299: PCGetOptionsPrefix(pc,&prefix);
1300: KSPSetOptionsPrefix(jac->ksp[0],prefix);
1301: KSPAppendOptionsPrefix(jac->ksp[0],"sub_");
1302: /*
1303: PetscMPIInt rank,subsize,subrank;
1304: MPI_Comm_rank(comm,&rank);
1305: MPI_Comm_size(subcomm,&subsize);
1306: MPI_Comm_rank(subcomm,&subrank);
1308: MatGetLocalSize(mpjac->submats,&m,PETSC_NULL);
1309: MatGetSize(mpjac->submats,&n,PETSC_NULL);
1310: PetscSynchronizedPrintf(comm,"[%d], sub-size %d,sub-rank %d\n",rank,subsize,subrank);
1311: PetscSynchronizedFlush(comm);
1312: */
1314: /* create dummy vectors xsub and ysub */
1315: MatGetLocalSize(mpjac->submats,&m,&n);
1316: VecCreateMPIWithArray(subcomm,1,n,PETSC_DECIDE,PETSC_NULL,&mpjac->xsub);
1317: VecCreateMPIWithArray(subcomm,1,m,PETSC_DECIDE,PETSC_NULL,&mpjac->ysub);
1318: PetscLogObjectParent(pc,mpjac->xsub);
1319: PetscLogObjectParent(pc,mpjac->ysub);
1321: pc->ops->reset = PCReset_BJacobi_Multiproc;
1322: pc->ops->destroy = PCDestroy_BJacobi_Multiproc;
1323: pc->ops->apply = PCApply_BJacobi_Multiproc;
1324: } else { /* pc->setupcalled */
1325: subcomm = mpjac->psubcomm->comm;
1326: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
1327: /* destroy old matrix blocks, then get new matrix blocks */
1328: if (mpjac->submats){MatDestroy(&mpjac->submats);}
1329: MatGetMultiProcBlock_MPIAIJ(pc->pmat,subcomm,MAT_INITIAL_MATRIX,&mpjac->submats);
1330: } else {
1331: MatGetMultiProcBlock_MPIAIJ(pc->pmat,subcomm,MAT_REUSE_MATRIX,&mpjac->submats);
1332: }
1333: KSPSetOperators(jac->ksp[0],mpjac->submats,mpjac->submats,pc->flag);
1334: }
1336: if (!wasSetup && pc->setfromoptionscalled){
1337: KSPSetFromOptions(jac->ksp[0]);
1338: }
1339: return(0);
1340: }