Actual source code: bjacobi.c
2: /*
3: Defines a block Jacobi preconditioner.
4: */
6: #include <../src/ksp/pc/impls/bjacobi/bjacobi.h>
8: static PetscErrorCode PCSetUp_BJacobi_Singleblock(PC,Mat,Mat);
9: static PetscErrorCode PCSetUp_BJacobi_Multiblock(PC,Mat,Mat);
10: static PetscErrorCode PCSetUp_BJacobi_Multiproc(PC);
12: static PetscErrorCode PCSetUp_BJacobi(PC pc)
13: {
14: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
15: Mat mat = pc->mat,pmat = pc->pmat;
17: PetscBool hasop;
18: PetscInt N,M,start,i,sum,end;
19: PetscInt bs,i_start=-1,i_end=-1;
20: PetscMPIInt rank,size;
23: MPI_Comm_rank(PetscObjectComm((PetscObject)pc),&rank);
24: MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);
25: MatGetLocalSize(pc->pmat,&M,&N);
26: MatGetBlockSize(pc->pmat,&bs);
28: if (jac->n > 0 && jac->n < size) {
29: PCSetUp_BJacobi_Multiproc(pc);
30: return(0);
31: }
33: /* --------------------------------------------------------------------------
34: Determines the number of blocks assigned to each processor
35: -----------------------------------------------------------------------------*/
37: /* local block count given */
38: if (jac->n_local > 0 && jac->n < 0) {
39: MPIU_Allreduce(&jac->n_local,&jac->n,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)pc));
40: if (jac->l_lens) { /* check that user set these correctly */
41: sum = 0;
42: for (i=0; i<jac->n_local; i++) {
43: 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");
44: sum += jac->l_lens[i];
45: }
46: if (sum != M) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Local lens set incorrectly");
47: } else {
48: PetscMalloc1(jac->n_local,&jac->l_lens);
49: for (i=0; i<jac->n_local; i++) jac->l_lens[i] = bs*((M/bs)/jac->n_local + (((M/bs) % jac->n_local) > i));
50: }
51: } else if (jac->n > 0 && jac->n_local < 0) { /* global block count given */
52: /* global blocks given: determine which ones are local */
53: if (jac->g_lens) {
54: /* check if the g_lens is has valid entries */
55: for (i=0; i<jac->n; i++) {
56: if (!jac->g_lens[i]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Zero block not allowed");
57: 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");
58: }
59: if (size == 1) {
60: jac->n_local = jac->n;
61: PetscMalloc1(jac->n_local,&jac->l_lens);
62: PetscArraycpy(jac->l_lens,jac->g_lens,jac->n_local);
63: /* check that user set these correctly */
64: sum = 0;
65: for (i=0; i<jac->n_local; i++) sum += jac->l_lens[i];
66: if (sum != M) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Global lens set incorrectly");
67: } else {
68: MatGetOwnershipRange(pc->pmat,&start,&end);
69: /* loop over blocks determing first one owned by me */
70: sum = 0;
71: for (i=0; i<jac->n+1; i++) {
72: if (sum == start) { i_start = i; goto start_1;}
73: if (i < jac->n) sum += jac->g_lens[i];
74: }
75: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Block sizes used in PCBJacobiSetTotalBlocks()\nare not compatible with parallel matrix layout");
76: start_1:
77: for (i=i_start; i<jac->n+1; i++) {
78: if (sum == end) { i_end = i; goto end_1; }
79: if (i < jac->n) sum += jac->g_lens[i];
80: }
81: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Block sizes used in PCBJacobiSetTotalBlocks()\nare not compatible with parallel matrix layout");
82: end_1:
83: jac->n_local = i_end - i_start;
84: PetscMalloc1(jac->n_local,&jac->l_lens);
85: PetscArraycpy(jac->l_lens,jac->g_lens+i_start,jac->n_local);
86: }
87: } else { /* no global blocks given, determine then using default layout */
88: jac->n_local = jac->n/size + ((jac->n % size) > rank);
89: PetscMalloc1(jac->n_local,&jac->l_lens);
90: for (i=0; i<jac->n_local; i++) {
91: jac->l_lens[i] = ((M/bs)/jac->n_local + (((M/bs) % jac->n_local) > i))*bs;
92: if (!jac->l_lens[i]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Too many blocks given");
93: }
94: }
95: } else if (jac->n < 0 && jac->n_local < 0) { /* no blocks given */
96: jac->n = size;
97: jac->n_local = 1;
98: PetscMalloc1(1,&jac->l_lens);
99: jac->l_lens[0] = M;
100: } else { /* jac->n > 0 && jac->n_local > 0 */
101: if (!jac->l_lens) {
102: PetscMalloc1(jac->n_local,&jac->l_lens);
103: for (i=0; i<jac->n_local; i++) jac->l_lens[i] = bs*((M/bs)/jac->n_local + (((M/bs) % jac->n_local) > i));
104: }
105: }
106: if (jac->n_local < 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Number of blocks is less than number of processors");
108: /* -------------------------
109: Determines mat and pmat
110: ---------------------------*/
111: MatHasOperation(pc->mat,MATOP_GET_DIAGONAL_BLOCK,&hasop);
112: if (!hasop && size == 1) {
113: mat = pc->mat;
114: pmat = pc->pmat;
115: } else {
116: if (pc->useAmat) {
117: /* use block from Amat matrix, not Pmat for local MatMult() */
118: MatGetDiagonalBlock(pc->mat,&mat);
119: }
120: if (pc->pmat != pc->mat || !pc->useAmat) {
121: MatGetDiagonalBlock(pc->pmat,&pmat);
122: } else pmat = mat;
123: }
125: /* ------
126: Setup code depends on the number of blocks
127: */
128: if (jac->n_local == 1) {
129: PCSetUp_BJacobi_Singleblock(pc,mat,pmat);
130: } else {
131: PCSetUp_BJacobi_Multiblock(pc,mat,pmat);
132: }
133: return(0);
134: }
136: /* Default destroy, if it has never been setup */
137: static PetscErrorCode PCDestroy_BJacobi(PC pc)
138: {
139: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
143: PetscFree(jac->g_lens);
144: PetscFree(jac->l_lens);
145: PetscFree(pc->data);
146: return(0);
147: }
150: static PetscErrorCode PCSetFromOptions_BJacobi(PetscOptionItems *PetscOptionsObject,PC pc)
151: {
152: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
154: PetscInt blocks,i;
155: PetscBool flg;
158: PetscOptionsHead(PetscOptionsObject,"Block Jacobi options");
159: PetscOptionsInt("-pc_bjacobi_blocks","Total number of blocks","PCBJacobiSetTotalBlocks",jac->n,&blocks,&flg);
160: if (flg) {PCBJacobiSetTotalBlocks(pc,blocks,NULL);}
161: PetscOptionsInt("-pc_bjacobi_local_blocks","Local number of blocks","PCBJacobiSetLocalBlocks",jac->n_local,&blocks,&flg);
162: if (flg) {PCBJacobiSetLocalBlocks(pc,blocks,NULL);}
163: if (jac->ksp) {
164: /* The sub-KSP has already been set up (e.g., PCSetUp_BJacobi_Singleblock), but KSPSetFromOptions was not called
165: * unless we had already been called. */
166: for (i=0; i<jac->n_local; i++) {
167: KSPSetFromOptions(jac->ksp[i]);
168: }
169: }
170: PetscOptionsTail();
171: return(0);
172: }
174: #include <petscdraw.h>
175: static PetscErrorCode PCView_BJacobi(PC pc,PetscViewer viewer)
176: {
177: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
178: PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc*)jac->data;
179: PetscErrorCode ierr;
180: PetscMPIInt rank;
181: PetscInt i;
182: PetscBool iascii,isstring,isdraw;
183: PetscViewer sviewer;
184: PetscViewerFormat format;
185: const char *prefix;
188: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
189: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
190: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
191: if (iascii) {
192: if (pc->useAmat) {
193: PetscViewerASCIIPrintf(viewer," using Amat local matrix, number of blocks = %D\n",jac->n);
194: }
195: PetscViewerASCIIPrintf(viewer," number of blocks = %D\n",jac->n);
196: MPI_Comm_rank(PetscObjectComm((PetscObject)pc),&rank);
197: PetscViewerGetFormat(viewer,&format);
198: if (format != PETSC_VIEWER_ASCII_INFO_DETAIL) {
199: PetscViewerASCIIPrintf(viewer," Local solver information for first block is in the following KSP and PC objects on rank 0:\n");
200: PCGetOptionsPrefix(pc,&prefix);
201: PetscViewerASCIIPrintf(viewer," Use -%sksp_view ::ascii_info_detail to display information for all blocks\n",prefix?prefix:"");
202: if (jac->ksp && !jac->psubcomm) {
203: PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
204: if (!rank) {
205: PetscViewerASCIIPushTab(viewer);
206: KSPView(jac->ksp[0],sviewer);
207: PetscViewerASCIIPopTab(viewer);
208: }
209: PetscViewerFlush(sviewer);
210: PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
211: PetscViewerFlush(viewer);
212: /* extra call needed because of the two calls to PetscViewerASCIIPushSynchronized() in PetscViewerGetSubViewer() */
213: PetscViewerASCIIPopSynchronized(viewer);
214: } else if (mpjac && jac->ksp && mpjac->psubcomm) {
215: PetscViewerGetSubViewer(viewer,mpjac->psubcomm->child,&sviewer);
216: if (!mpjac->psubcomm->color) {
217: PetscViewerASCIIPushTab(viewer);
218: KSPView(*(jac->ksp),sviewer);
219: PetscViewerASCIIPopTab(viewer);
220: }
221: PetscViewerFlush(sviewer);
222: PetscViewerRestoreSubViewer(viewer,mpjac->psubcomm->child,&sviewer);
223: PetscViewerFlush(viewer);
224: /* extra call needed because of the two calls to PetscViewerASCIIPushSynchronized() in PetscViewerGetSubViewer() */
225: PetscViewerASCIIPopSynchronized(viewer);
226: } else {
227: PetscViewerFlush(viewer);
228: }
229: } else {
230: PetscInt n_global;
231: MPIU_Allreduce(&jac->n_local,&n_global,1,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)pc));
232: PetscViewerASCIIPushSynchronized(viewer);
233: PetscViewerASCIIPrintf(viewer," Local solver information for each block is in the following KSP and PC objects:\n");
234: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] number of local blocks = %D, first local block number = %D\n",
235: rank,jac->n_local,jac->first_local);
236: PetscViewerASCIIPushTab(viewer);
237: PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
238: for (i=0; i<jac->n_local; i++) {
239: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] local block number %D\n",rank,i);
240: KSPView(jac->ksp[i],sviewer);
241: PetscViewerASCIISynchronizedPrintf(viewer,"- - - - - - - - - - - - - - - - - -\n");
242: }
243: PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
244: PetscViewerASCIIPopTab(viewer);
245: PetscViewerFlush(viewer);
246: PetscViewerASCIIPopSynchronized(viewer);
247: }
248: } else if (isstring) {
249: PetscViewerStringSPrintf(viewer," blks=%D",jac->n);
250: PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
251: if (jac->ksp) {KSPView(jac->ksp[0],sviewer);}
252: PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
253: } else if (isdraw) {
254: PetscDraw draw;
255: char str[25];
256: PetscReal x,y,bottom,h;
258: PetscViewerDrawGetDraw(viewer,0,&draw);
259: PetscDrawGetCurrentPoint(draw,&x,&y);
260: PetscSNPrintf(str,25,"Number blocks %D",jac->n);
261: PetscDrawStringBoxed(draw,x,y,PETSC_DRAW_RED,PETSC_DRAW_BLACK,str,NULL,&h);
262: bottom = y - h;
263: PetscDrawPushCurrentPoint(draw,x,bottom);
264: /* warning the communicator on viewer is different then on ksp in parallel */
265: if (jac->ksp) {KSPView(jac->ksp[0],viewer);}
266: PetscDrawPopCurrentPoint(draw);
267: }
268: return(0);
269: }
271: /* -------------------------------------------------------------------------------------*/
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: if (ksp) *ksp = jac->ksp;
283: return(0);
284: }
286: static PetscErrorCode PCBJacobiSetTotalBlocks_BJacobi(PC pc,PetscInt blocks,PetscInt *lens)
287: {
288: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
292: 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");
293: jac->n = blocks;
294: if (!lens) jac->g_lens = NULL;
295: else {
296: PetscMalloc1(blocks,&jac->g_lens);
297: PetscLogObjectMemory((PetscObject)pc,blocks*sizeof(PetscInt));
298: PetscArraycpy(jac->g_lens,lens,blocks);
299: }
300: return(0);
301: }
303: static PetscErrorCode PCBJacobiGetTotalBlocks_BJacobi(PC pc, PetscInt *blocks, const PetscInt *lens[])
304: {
305: PC_BJacobi *jac = (PC_BJacobi*) pc->data;
308: *blocks = jac->n;
309: if (lens) *lens = jac->g_lens;
310: return(0);
311: }
313: static PetscErrorCode PCBJacobiSetLocalBlocks_BJacobi(PC pc,PetscInt blocks,const PetscInt lens[])
314: {
315: PC_BJacobi *jac;
319: jac = (PC_BJacobi*)pc->data;
321: jac->n_local = blocks;
322: if (!lens) jac->l_lens = NULL;
323: else {
324: PetscMalloc1(blocks,&jac->l_lens);
325: PetscLogObjectMemory((PetscObject)pc,blocks*sizeof(PetscInt));
326: PetscArraycpy(jac->l_lens,lens,blocks);
327: }
328: return(0);
329: }
331: static PetscErrorCode PCBJacobiGetLocalBlocks_BJacobi(PC pc, PetscInt *blocks, const PetscInt *lens[])
332: {
333: PC_BJacobi *jac = (PC_BJacobi*) pc->data;
336: *blocks = jac->n_local;
337: if (lens) *lens = jac->l_lens;
338: return(0);
339: }
341: /* -------------------------------------------------------------------------------------*/
343: /*@C
344: PCBJacobiGetSubKSP - Gets the local KSP contexts for all blocks on
345: this processor.
347: Not Collective
349: Input Parameter:
350: . pc - the preconditioner context
352: Output Parameters:
353: + n_local - the number of blocks on this processor, or NULL
354: . first_local - the global number of the first block on this processor, or NULL
355: - ksp - the array of KSP contexts
357: Notes:
358: After PCBJacobiGetSubKSP() the array of KSP contexts is not to be freed.
360: Currently for some matrix implementations only 1 block per processor
361: is supported.
363: You must call KSPSetUp() or PCSetUp() before calling PCBJacobiGetSubKSP().
365: Fortran Usage: You must pass in a KSP array that is large enough to contain all the local KSPs.
366: You can call PCBJacobiGetSubKSP(pc,nlocal,firstlocal,PETSC_NULL_KSP,ierr) to determine how large the
367: KSP array must be.
369: Level: advanced
371: .seealso: PCASMGetSubKSP()
372: @*/
373: PetscErrorCode PCBJacobiGetSubKSP(PC pc,PetscInt *n_local,PetscInt *first_local,KSP *ksp[])
374: {
379: PetscUseMethod(pc,"PCBJacobiGetSubKSP_C",(PC,PetscInt*,PetscInt*,KSP **),(pc,n_local,first_local,ksp));
380: return(0);
381: }
383: /*@
384: PCBJacobiSetTotalBlocks - Sets the global number of blocks for the block
385: Jacobi preconditioner.
387: Collective on PC
389: Input Parameters:
390: + pc - the preconditioner context
391: . blocks - the number of blocks
392: - lens - [optional] integer array containing the size of each block
394: Options Database Key:
395: . -pc_bjacobi_blocks <blocks> - Sets the number of global blocks
397: Notes:
398: Currently only a limited number of blocking configurations are supported.
399: All processors sharing the PC must call this routine with the same data.
401: Level: intermediate
403: .seealso: PCSetUseAmat(), PCBJacobiSetLocalBlocks()
404: @*/
405: PetscErrorCode PCBJacobiSetTotalBlocks(PC pc,PetscInt blocks,const PetscInt lens[])
406: {
411: if (blocks <= 0) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Must have positive blocks");
412: PetscTryMethod(pc,"PCBJacobiSetTotalBlocks_C",(PC,PetscInt,const PetscInt[]),(pc,blocks,lens));
413: return(0);
414: }
416: /*@C
417: PCBJacobiGetTotalBlocks - Gets the global number of blocks for the block
418: Jacobi preconditioner.
420: Not Collective
422: Input Parameter:
423: . pc - the preconditioner context
425: Output parameters:
426: + blocks - the number of blocks
427: - lens - integer array containing the size of each block
429: Level: intermediate
431: .seealso: PCSetUseAmat(), PCBJacobiGetLocalBlocks()
432: @*/
433: PetscErrorCode PCBJacobiGetTotalBlocks(PC pc, PetscInt *blocks, const PetscInt *lens[])
434: {
440: PetscUseMethod(pc,"PCBJacobiGetTotalBlocks_C",(PC,PetscInt*, const PetscInt *[]),(pc,blocks,lens));
441: return(0);
442: }
444: /*@
445: PCBJacobiSetLocalBlocks - Sets the local number of blocks for the block
446: Jacobi preconditioner.
448: Not Collective
450: Input Parameters:
451: + pc - the preconditioner context
452: . blocks - the number of blocks
453: - lens - [optional] integer array containing size of each block
455: Options Database Key:
456: . -pc_bjacobi_local_blocks <blocks> - Sets the number of local blocks
458: Note:
459: Currently only a limited number of blocking configurations are supported.
461: Level: intermediate
463: .seealso: PCSetUseAmat(), PCBJacobiSetTotalBlocks()
464: @*/
465: PetscErrorCode PCBJacobiSetLocalBlocks(PC pc,PetscInt blocks,const PetscInt lens[])
466: {
471: if (blocks < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must have nonegative blocks");
472: PetscTryMethod(pc,"PCBJacobiSetLocalBlocks_C",(PC,PetscInt,const PetscInt []),(pc,blocks,lens));
473: return(0);
474: }
476: /*@C
477: PCBJacobiGetLocalBlocks - Gets the local number of blocks for the block
478: Jacobi preconditioner.
480: Not Collective
482: Input Parameters:
483: + pc - the preconditioner context
484: . blocks - the number of blocks
485: - lens - [optional] integer array containing size of each block
487: Note:
488: Currently only a limited number of blocking configurations are supported.
490: Level: intermediate
492: .seealso: PCSetUseAmat(), PCBJacobiGetTotalBlocks()
493: @*/
494: PetscErrorCode PCBJacobiGetLocalBlocks(PC pc, PetscInt *blocks, const PetscInt *lens[])
495: {
501: PetscUseMethod(pc,"PCBJacobiGetLocalBlocks_C",(PC,PetscInt*, const PetscInt *[]),(pc,blocks,lens));
502: return(0);
503: }
505: /* -----------------------------------------------------------------------------------*/
507: /*MC
508: PCBJACOBI - Use block Jacobi preconditioning, each block is (approximately) solved with
509: its own KSP object.
511: Options Database Keys:
512: + -pc_use_amat - use Amat to apply block of operator in inner Krylov method
513: - -pc_bjacobi_blocks <n> - use n total blocks
515: Notes:
516: Each processor can have one or more blocks, or a single block can be shared by several processes. Defaults to one block per processor.
518: To set options on the solvers for each block append -sub_ to all the KSP, KSP, and PC
519: options database keys. For example, -sub_pc_type ilu -sub_pc_factor_levels 1 -sub_ksp_type preonly
521: To set the options on the solvers separate for each block call PCBJacobiGetSubKSP()
522: and set the options directly on the resulting KSP object (you can access its PC
523: KSPGetPC())
525: For GPU-based vectors (CUDA, ViennaCL) it is recommended to use exactly one block per MPI process for best
526: performance. Different block partitioning may lead to additional data transfers
527: between host and GPU that lead to degraded performance.
529: The options prefix for each block is sub_, for example -sub_pc_type lu.
531: When multiple processes share a single block, each block encompasses exactly all the unknowns owned its set of processes.
533: Level: beginner
535: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
536: PCASM, PCSetUseAmat(), PCGetUseAmat(), PCBJacobiGetSubKSP(), PCBJacobiSetTotalBlocks(),
537: PCBJacobiSetLocalBlocks(), PCSetModifySubMatrices()
538: M*/
540: PETSC_EXTERN PetscErrorCode PCCreate_BJacobi(PC pc)
541: {
543: PetscMPIInt rank;
544: PC_BJacobi *jac;
547: PetscNewLog(pc,&jac);
548: MPI_Comm_rank(PetscObjectComm((PetscObject)pc),&rank);
550: pc->ops->apply = NULL;
551: pc->ops->matapply = NULL;
552: pc->ops->applytranspose = NULL;
553: pc->ops->setup = PCSetUp_BJacobi;
554: pc->ops->destroy = PCDestroy_BJacobi;
555: pc->ops->setfromoptions = PCSetFromOptions_BJacobi;
556: pc->ops->view = PCView_BJacobi;
557: pc->ops->applyrichardson = NULL;
559: pc->data = (void*)jac;
560: jac->n = -1;
561: jac->n_local = -1;
562: jac->first_local = rank;
563: jac->ksp = NULL;
564: jac->g_lens = NULL;
565: jac->l_lens = NULL;
566: jac->psubcomm = NULL;
568: PetscObjectComposeFunction((PetscObject)pc,"PCBJacobiGetSubKSP_C",PCBJacobiGetSubKSP_BJacobi);
569: PetscObjectComposeFunction((PetscObject)pc,"PCBJacobiSetTotalBlocks_C",PCBJacobiSetTotalBlocks_BJacobi);
570: PetscObjectComposeFunction((PetscObject)pc,"PCBJacobiGetTotalBlocks_C",PCBJacobiGetTotalBlocks_BJacobi);
571: PetscObjectComposeFunction((PetscObject)pc,"PCBJacobiSetLocalBlocks_C",PCBJacobiSetLocalBlocks_BJacobi);
572: PetscObjectComposeFunction((PetscObject)pc,"PCBJacobiGetLocalBlocks_C",PCBJacobiGetLocalBlocks_BJacobi);
573: return(0);
574: }
576: /* --------------------------------------------------------------------------------------------*/
577: /*
578: These are for a single block per processor; works for AIJ, BAIJ; Seq and MPI
579: */
580: static PetscErrorCode PCReset_BJacobi_Singleblock(PC pc)
581: {
582: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
583: PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock*)jac->data;
584: PetscErrorCode ierr;
587: KSPReset(jac->ksp[0]);
588: VecDestroy(&bjac->x);
589: VecDestroy(&bjac->y);
590: return(0);
591: }
593: static PetscErrorCode PCDestroy_BJacobi_Singleblock(PC pc)
594: {
595: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
596: PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock*)jac->data;
597: PetscErrorCode ierr;
600: PCReset_BJacobi_Singleblock(pc);
601: KSPDestroy(&jac->ksp[0]);
602: PetscFree(jac->ksp);
603: PetscFree(jac->l_lens);
604: PetscFree(jac->g_lens);
605: PetscFree(bjac);
606: PetscFree(pc->data);
607: return(0);
608: }
610: static PetscErrorCode PCSetUpOnBlocks_BJacobi_Singleblock(PC pc)
611: {
612: PetscErrorCode ierr;
613: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
614: KSP subksp = jac->ksp[0];
615: KSPConvergedReason reason;
618: KSPSetUp(subksp);
619: KSPGetConvergedReason(subksp,&reason);
620: if (reason == KSP_DIVERGED_PC_FAILED) {
621: pc->failedreason = PC_SUBPC_ERROR;
622: }
623: return(0);
624: }
626: static PetscErrorCode PCApply_BJacobi_Singleblock(PC pc,Vec x,Vec y)
627: {
628: PetscErrorCode ierr;
629: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
630: PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock*)jac->data;
633: VecGetLocalVectorRead(x, bjac->x);
634: VecGetLocalVector(y, bjac->y);
635: /* Since the inner KSP matrix may point directly to the diagonal block of an MPI matrix the inner
636: matrix may change even if the outer KSP/PC has not updated the preconditioner, this will trigger a rebuild
637: of the inner preconditioner automatically unless we pass down the outer preconditioners reuse flag.*/
638: KSPSetReusePreconditioner(jac->ksp[0],pc->reusepreconditioner);
639: KSPSolve(jac->ksp[0],bjac->x,bjac->y);
640: KSPCheckSolve(jac->ksp[0],pc,bjac->y);
641: VecRestoreLocalVectorRead(x, bjac->x);
642: VecRestoreLocalVector(y, bjac->y);
643: return(0);
644: }
646: static PetscErrorCode PCMatApply_BJacobi_Singleblock(PC pc,Mat X,Mat Y)
647: {
648: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
649: Mat sX,sY;
653: /* Since the inner KSP matrix may point directly to the diagonal block of an MPI matrix the inner
654: matrix may change even if the outer KSP/PC has not updated the preconditioner, this will trigger a rebuild
655: of the inner preconditioner automatically unless we pass down the outer preconditioners reuse flag.*/
656: KSPSetReusePreconditioner(jac->ksp[0],pc->reusepreconditioner);
657: MatDenseGetLocalMatrix(X,&sX);
658: MatDenseGetLocalMatrix(Y,&sY);
659: KSPMatSolve(jac->ksp[0],sX,sY);
660: return(0);
661: }
663: static PetscErrorCode PCApplySymmetricLeft_BJacobi_Singleblock(PC pc,Vec x,Vec y)
664: {
665: PetscErrorCode ierr;
666: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
667: PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock*)jac->data;
668: PetscScalar *y_array;
669: const PetscScalar *x_array;
670: PC subpc;
673: /*
674: The VecPlaceArray() is to avoid having to copy the
675: y vector into the bjac->x vector. The reason for
676: the bjac->x vector is that we need a sequential vector
677: for the sequential solve.
678: */
679: VecGetArrayRead(x,&x_array);
680: VecGetArray(y,&y_array);
681: VecPlaceArray(bjac->x,x_array);
682: VecPlaceArray(bjac->y,y_array);
683: /* apply the symmetric left portion of the inner PC operator */
684: /* note this by-passes the inner KSP and its options completely */
685: KSPGetPC(jac->ksp[0],&subpc);
686: PCApplySymmetricLeft(subpc,bjac->x,bjac->y);
687: VecResetArray(bjac->x);
688: VecResetArray(bjac->y);
689: VecRestoreArrayRead(x,&x_array);
690: VecRestoreArray(y,&y_array);
691: return(0);
692: }
694: static PetscErrorCode PCApplySymmetricRight_BJacobi_Singleblock(PC pc,Vec x,Vec y)
695: {
696: PetscErrorCode ierr;
697: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
698: PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock*)jac->data;
699: PetscScalar *y_array;
700: const PetscScalar *x_array;
701: PC subpc;
704: /*
705: The VecPlaceArray() is to avoid having to copy the
706: y vector into the bjac->x vector. The reason for
707: the bjac->x vector is that we need a sequential vector
708: for the sequential solve.
709: */
710: VecGetArrayRead(x,&x_array);
711: VecGetArray(y,&y_array);
712: VecPlaceArray(bjac->x,x_array);
713: VecPlaceArray(bjac->y,y_array);
715: /* apply the symmetric right portion of the inner PC operator */
716: /* note this by-passes the inner KSP and its options completely */
718: KSPGetPC(jac->ksp[0],&subpc);
719: PCApplySymmetricRight(subpc,bjac->x,bjac->y);
721: VecRestoreArrayRead(x,&x_array);
722: VecRestoreArray(y,&y_array);
723: return(0);
724: }
726: static PetscErrorCode PCApplyTranspose_BJacobi_Singleblock(PC pc,Vec x,Vec y)
727: {
728: PetscErrorCode ierr;
729: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
730: PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock*)jac->data;
731: PetscScalar *y_array;
732: const PetscScalar *x_array;
735: /*
736: The VecPlaceArray() is to avoid having to copy the
737: y vector into the bjac->x vector. The reason for
738: the bjac->x vector is that we need a sequential vector
739: for the sequential solve.
740: */
741: VecGetArrayRead(x,&x_array);
742: VecGetArray(y,&y_array);
743: VecPlaceArray(bjac->x,x_array);
744: VecPlaceArray(bjac->y,y_array);
745: KSPSolveTranspose(jac->ksp[0],bjac->x,bjac->y);
746: KSPCheckSolve(jac->ksp[0],pc,bjac->y);
747: VecResetArray(bjac->x);
748: VecResetArray(bjac->y);
749: VecRestoreArrayRead(x,&x_array);
750: VecRestoreArray(y,&y_array);
751: return(0);
752: }
754: static PetscErrorCode PCSetUp_BJacobi_Singleblock(PC pc,Mat mat,Mat pmat)
755: {
756: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
757: PetscErrorCode ierr;
758: PetscInt m;
759: KSP ksp;
760: PC_BJacobi_Singleblock *bjac;
761: PetscBool wasSetup = PETSC_TRUE;
762: VecType vectype;
763: const char *prefix;
766: if (!pc->setupcalled) {
768: if (!jac->ksp) {
769: wasSetup = PETSC_FALSE;
771: KSPCreate(PETSC_COMM_SELF,&ksp);
772: KSPSetErrorIfNotConverged(ksp,pc->erroriffailure);
773: PetscObjectIncrementTabLevel((PetscObject)ksp,(PetscObject)pc,1);
774: PetscLogObjectParent((PetscObject)pc,(PetscObject)ksp);
775: KSPSetType(ksp,KSPPREONLY);
776: PCGetOptionsPrefix(pc,&prefix);
777: KSPSetOptionsPrefix(ksp,prefix);
778: KSPAppendOptionsPrefix(ksp,"sub_");
780: pc->ops->reset = PCReset_BJacobi_Singleblock;
781: pc->ops->destroy = PCDestroy_BJacobi_Singleblock;
782: pc->ops->apply = PCApply_BJacobi_Singleblock;
783: pc->ops->matapply = PCMatApply_BJacobi_Singleblock;
784: pc->ops->applysymmetricleft = PCApplySymmetricLeft_BJacobi_Singleblock;
785: pc->ops->applysymmetricright = PCApplySymmetricRight_BJacobi_Singleblock;
786: pc->ops->applytranspose = PCApplyTranspose_BJacobi_Singleblock;
787: pc->ops->setuponblocks = PCSetUpOnBlocks_BJacobi_Singleblock;
789: PetscMalloc1(1,&jac->ksp);
790: jac->ksp[0] = ksp;
792: PetscNewLog(pc,&bjac);
793: jac->data = (void*)bjac;
794: } else {
795: ksp = jac->ksp[0];
796: bjac = (PC_BJacobi_Singleblock*)jac->data;
797: }
799: /*
800: The reason we need to generate these vectors is to serve
801: as the right-hand side and solution vector for the solve on the
802: block. We do not need to allocate space for the vectors since
803: that is provided via VecPlaceArray() just before the call to
804: KSPSolve() on the block.
805: */
806: MatGetSize(pmat,&m,&m);
807: VecCreateSeqWithArray(PETSC_COMM_SELF,1,m,NULL,&bjac->x);
808: VecCreateSeqWithArray(PETSC_COMM_SELF,1,m,NULL,&bjac->y);
809: MatGetVecType(pmat,&vectype);
810: VecSetType(bjac->x,vectype);
811: VecSetType(bjac->y,vectype);
812: PetscLogObjectParent((PetscObject)pc,(PetscObject)bjac->x);
813: PetscLogObjectParent((PetscObject)pc,(PetscObject)bjac->y);
814: } else {
815: ksp = jac->ksp[0];
816: bjac = (PC_BJacobi_Singleblock*)jac->data;
817: }
818: KSPGetOptionsPrefix(ksp,&prefix);
819: if (pc->useAmat) {
820: KSPSetOperators(ksp,mat,pmat);
821: MatSetOptionsPrefix(mat,prefix);
822: } else {
823: KSPSetOperators(ksp,pmat,pmat);
824: }
825: MatSetOptionsPrefix(pmat,prefix);
826: if (!wasSetup && pc->setfromoptionscalled) {
827: /* If PCSetFromOptions_BJacobi is called later, KSPSetFromOptions will be called at that time. */
828: KSPSetFromOptions(ksp);
829: }
830: return(0);
831: }
833: /* ---------------------------------------------------------------------------------------------*/
834: static PetscErrorCode PCReset_BJacobi_Multiblock(PC pc)
835: {
836: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
837: PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock*)jac->data;
838: PetscErrorCode ierr;
839: PetscInt i;
842: if (bjac && bjac->pmat) {
843: MatDestroyMatrices(jac->n_local,&bjac->pmat);
844: if (pc->useAmat) {
845: MatDestroyMatrices(jac->n_local,&bjac->mat);
846: }
847: }
849: for (i=0; i<jac->n_local; i++) {
850: KSPReset(jac->ksp[i]);
851: if (bjac && bjac->x) {
852: VecDestroy(&bjac->x[i]);
853: VecDestroy(&bjac->y[i]);
854: ISDestroy(&bjac->is[i]);
855: }
856: }
857: PetscFree(jac->l_lens);
858: PetscFree(jac->g_lens);
859: return(0);
860: }
862: static PetscErrorCode PCDestroy_BJacobi_Multiblock(PC pc)
863: {
864: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
865: PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock*)jac->data;
866: PetscErrorCode ierr;
867: PetscInt i;
870: PCReset_BJacobi_Multiblock(pc);
871: if (bjac) {
872: PetscFree2(bjac->x,bjac->y);
873: PetscFree(bjac->starts);
874: PetscFree(bjac->is);
875: }
876: PetscFree(jac->data);
877: for (i=0; i<jac->n_local; i++) {
878: KSPDestroy(&jac->ksp[i]);
879: }
880: PetscFree(jac->ksp);
881: PetscFree(pc->data);
882: return(0);
883: }
885: static PetscErrorCode PCSetUpOnBlocks_BJacobi_Multiblock(PC pc)
886: {
887: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
888: PetscErrorCode ierr;
889: PetscInt i,n_local = jac->n_local;
890: KSPConvergedReason reason;
893: for (i=0; i<n_local; i++) {
894: KSPSetUp(jac->ksp[i]);
895: KSPGetConvergedReason(jac->ksp[i],&reason);
896: if (reason == KSP_DIVERGED_PC_FAILED) {
897: pc->failedreason = PC_SUBPC_ERROR;
898: }
899: }
900: return(0);
901: }
903: /*
904: Preconditioner for block Jacobi
905: */
906: static PetscErrorCode PCApply_BJacobi_Multiblock(PC pc,Vec x,Vec y)
907: {
908: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
909: PetscErrorCode ierr;
910: PetscInt i,n_local = jac->n_local;
911: PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock*)jac->data;
912: PetscScalar *yin;
913: const PetscScalar *xin;
916: VecGetArrayRead(x,&xin);
917: VecGetArray(y,&yin);
918: for (i=0; i<n_local; i++) {
919: /*
920: To avoid copying the subvector from x into a workspace we instead
921: make the workspace vector array point to the subpart of the array of
922: the global vector.
923: */
924: VecPlaceArray(bjac->x[i],xin+bjac->starts[i]);
925: VecPlaceArray(bjac->y[i],yin+bjac->starts[i]);
927: PetscLogEventBegin(PC_ApplyOnBlocks,jac->ksp[i],bjac->x[i],bjac->y[i],0);
928: KSPSolve(jac->ksp[i],bjac->x[i],bjac->y[i]);
929: KSPCheckSolve(jac->ksp[i],pc,bjac->y[i]);
930: PetscLogEventEnd(PC_ApplyOnBlocks,jac->ksp[i],bjac->x[i],bjac->y[i],0);
932: VecResetArray(bjac->x[i]);
933: VecResetArray(bjac->y[i]);
934: }
935: VecRestoreArrayRead(x,&xin);
936: VecRestoreArray(y,&yin);
937: return(0);
938: }
940: /*
941: Preconditioner for block Jacobi
942: */
943: static PetscErrorCode PCApplyTranspose_BJacobi_Multiblock(PC pc,Vec x,Vec y)
944: {
945: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
946: PetscErrorCode ierr;
947: PetscInt i,n_local = jac->n_local;
948: PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock*)jac->data;
949: PetscScalar *yin;
950: const PetscScalar *xin;
953: VecGetArrayRead(x,&xin);
954: VecGetArray(y,&yin);
955: for (i=0; i<n_local; i++) {
956: /*
957: To avoid copying the subvector from x into a workspace we instead
958: make the workspace vector array point to the subpart of the array of
959: the global vector.
960: */
961: VecPlaceArray(bjac->x[i],xin+bjac->starts[i]);
962: VecPlaceArray(bjac->y[i],yin+bjac->starts[i]);
964: PetscLogEventBegin(PC_ApplyTransposeOnBlocks,jac->ksp[i],bjac->x[i],bjac->y[i],0);
965: KSPSolveTranspose(jac->ksp[i],bjac->x[i],bjac->y[i]);
966: KSPCheckSolve(jac->ksp[i],pc,bjac->y[i]);
967: PetscLogEventEnd(PC_ApplyTransposeOnBlocks,jac->ksp[i],bjac->x[i],bjac->y[i],0);
969: VecResetArray(bjac->x[i]);
970: VecResetArray(bjac->y[i]);
971: }
972: VecRestoreArrayRead(x,&xin);
973: VecRestoreArray(y,&yin);
974: return(0);
975: }
977: static PetscErrorCode PCSetUp_BJacobi_Multiblock(PC pc,Mat mat,Mat pmat)
978: {
979: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
980: PetscErrorCode ierr;
981: PetscInt m,n_local,N,M,start,i;
982: const char *prefix;
983: KSP ksp;
984: Vec x,y;
985: PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock*)jac->data;
986: PC subpc;
987: IS is;
988: MatReuse scall;
989: VecType vectype;
992: MatGetLocalSize(pc->pmat,&M,&N);
994: n_local = jac->n_local;
996: if (pc->useAmat) {
997: PetscBool same;
998: PetscObjectTypeCompare((PetscObject)mat,((PetscObject)pmat)->type_name,&same);
999: if (!same) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"Matrices not of same type");
1000: }
1002: if (!pc->setupcalled) {
1003: scall = MAT_INITIAL_MATRIX;
1005: if (!jac->ksp) {
1006: pc->ops->reset = PCReset_BJacobi_Multiblock;
1007: pc->ops->destroy = PCDestroy_BJacobi_Multiblock;
1008: pc->ops->apply = PCApply_BJacobi_Multiblock;
1009: pc->ops->matapply = NULL;
1010: pc->ops->applytranspose= PCApplyTranspose_BJacobi_Multiblock;
1011: pc->ops->setuponblocks = PCSetUpOnBlocks_BJacobi_Multiblock;
1013: PetscNewLog(pc,&bjac);
1014: PetscMalloc1(n_local,&jac->ksp);
1015: PetscLogObjectMemory((PetscObject)pc,sizeof(n_local*sizeof(KSP)));
1016: PetscMalloc2(n_local,&bjac->x,n_local,&bjac->y);
1017: PetscMalloc1(n_local,&bjac->starts);
1018: PetscLogObjectMemory((PetscObject)pc,sizeof(n_local*sizeof(PetscScalar)));
1020: jac->data = (void*)bjac;
1021: PetscMalloc1(n_local,&bjac->is);
1022: PetscLogObjectMemory((PetscObject)pc,sizeof(n_local*sizeof(IS)));
1024: for (i=0; i<n_local; i++) {
1025: KSPCreate(PETSC_COMM_SELF,&ksp);
1026: KSPSetErrorIfNotConverged(ksp,pc->erroriffailure);
1027: PetscObjectIncrementTabLevel((PetscObject)ksp,(PetscObject)pc,1);
1028: PetscLogObjectParent((PetscObject)pc,(PetscObject)ksp);
1029: KSPSetType(ksp,KSPPREONLY);
1030: KSPGetPC(ksp,&subpc);
1031: PCGetOptionsPrefix(pc,&prefix);
1032: KSPSetOptionsPrefix(ksp,prefix);
1033: KSPAppendOptionsPrefix(ksp,"sub_");
1035: jac->ksp[i] = ksp;
1036: }
1037: } else {
1038: bjac = (PC_BJacobi_Multiblock*)jac->data;
1039: }
1041: start = 0;
1042: MatGetVecType(pmat,&vectype);
1043: for (i=0; i<n_local; i++) {
1044: m = jac->l_lens[i];
1045: /*
1046: The reason we need to generate these vectors is to serve
1047: as the right-hand side and solution vector for the solve on the
1048: block. We do not need to allocate space for the vectors since
1049: that is provided via VecPlaceArray() just before the call to
1050: KSPSolve() on the block.
1052: */
1053: VecCreateSeq(PETSC_COMM_SELF,m,&x);
1054: VecCreateSeqWithArray(PETSC_COMM_SELF,1,m,NULL,&y);
1055: VecSetType(x,vectype);
1056: VecSetType(y,vectype);
1057: PetscLogObjectParent((PetscObject)pc,(PetscObject)x);
1058: PetscLogObjectParent((PetscObject)pc,(PetscObject)y);
1060: bjac->x[i] = x;
1061: bjac->y[i] = y;
1062: bjac->starts[i] = start;
1064: ISCreateStride(PETSC_COMM_SELF,m,start,1,&is);
1065: bjac->is[i] = is;
1066: PetscLogObjectParent((PetscObject)pc,(PetscObject)is);
1068: start += m;
1069: }
1070: } else {
1071: bjac = (PC_BJacobi_Multiblock*)jac->data;
1072: /*
1073: Destroy the blocks from the previous iteration
1074: */
1075: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
1076: MatDestroyMatrices(n_local,&bjac->pmat);
1077: if (pc->useAmat) {
1078: MatDestroyMatrices(n_local,&bjac->mat);
1079: }
1080: scall = MAT_INITIAL_MATRIX;
1081: } else scall = MAT_REUSE_MATRIX;
1082: }
1084: MatCreateSubMatrices(pmat,n_local,bjac->is,bjac->is,scall,&bjac->pmat);
1085: if (pc->useAmat) {
1086: MatCreateSubMatrices(mat,n_local,bjac->is,bjac->is,scall,&bjac->mat);
1087: }
1088: /* Return control to the user so that the submatrices can be modified (e.g., to apply
1089: different boundary conditions for the submatrices than for the global problem) */
1090: PCModifySubMatrices(pc,n_local,bjac->is,bjac->is,bjac->pmat,pc->modifysubmatricesP);
1092: for (i=0; i<n_local; i++) {
1093: PetscLogObjectParent((PetscObject)pc,(PetscObject)bjac->pmat[i]);
1094: KSPGetOptionsPrefix(jac->ksp[i],&prefix);
1095: if (pc->useAmat) {
1096: PetscLogObjectParent((PetscObject)pc,(PetscObject)bjac->mat[i]);
1097: KSPSetOperators(jac->ksp[i],bjac->mat[i],bjac->pmat[i]);
1098: MatSetOptionsPrefix(bjac->mat[i],prefix);
1099: } else {
1100: KSPSetOperators(jac->ksp[i],bjac->pmat[i],bjac->pmat[i]);
1101: }
1102: MatSetOptionsPrefix(bjac->pmat[i],prefix);
1103: if (pc->setfromoptionscalled) {
1104: KSPSetFromOptions(jac->ksp[i]);
1105: }
1106: }
1107: return(0);
1108: }
1110: /* ---------------------------------------------------------------------------------------------*/
1111: /*
1112: These are for a single block with multiple processes;
1113: */
1114: static PetscErrorCode PCReset_BJacobi_Multiproc(PC pc)
1115: {
1116: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
1117: PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc*)jac->data;
1118: PetscErrorCode ierr;
1121: VecDestroy(&mpjac->ysub);
1122: VecDestroy(&mpjac->xsub);
1123: MatDestroy(&mpjac->submats);
1124: if (jac->ksp) {KSPReset(jac->ksp[0]);}
1125: return(0);
1126: }
1128: static PetscErrorCode PCDestroy_BJacobi_Multiproc(PC pc)
1129: {
1130: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
1131: PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc*)jac->data;
1132: PetscErrorCode ierr;
1135: PCReset_BJacobi_Multiproc(pc);
1136: KSPDestroy(&jac->ksp[0]);
1137: PetscFree(jac->ksp);
1138: PetscSubcommDestroy(&mpjac->psubcomm);
1140: PetscFree(mpjac);
1141: PetscFree(pc->data);
1142: return(0);
1143: }
1145: static PetscErrorCode PCApply_BJacobi_Multiproc(PC pc,Vec x,Vec y)
1146: {
1147: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
1148: PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc*)jac->data;
1149: PetscErrorCode ierr;
1150: PetscScalar *yarray;
1151: const PetscScalar *xarray;
1152: KSPConvergedReason reason;
1155: /* place x's and y's local arrays into xsub and ysub */
1156: VecGetArrayRead(x,&xarray);
1157: VecGetArray(y,&yarray);
1158: VecPlaceArray(mpjac->xsub,xarray);
1159: VecPlaceArray(mpjac->ysub,yarray);
1161: /* apply preconditioner on each matrix block */
1162: PetscLogEventBegin(PC_ApplyOnBlocks,jac->ksp[0],mpjac->xsub,mpjac->ysub,0);
1163: KSPSolve(jac->ksp[0],mpjac->xsub,mpjac->ysub);
1164: KSPCheckSolve(jac->ksp[0],pc,mpjac->ysub);
1165: PetscLogEventEnd(PC_ApplyOnBlocks,jac->ksp[0],mpjac->xsub,mpjac->ysub,0);
1166: KSPGetConvergedReason(jac->ksp[0],&reason);
1167: if (reason == KSP_DIVERGED_PC_FAILED) {
1168: pc->failedreason = PC_SUBPC_ERROR;
1169: }
1171: VecResetArray(mpjac->xsub);
1172: VecResetArray(mpjac->ysub);
1173: VecRestoreArrayRead(x,&xarray);
1174: VecRestoreArray(y,&yarray);
1175: return(0);
1176: }
1178: static PetscErrorCode PCMatApply_BJacobi_Multiproc(PC pc,Mat X,Mat Y)
1179: {
1180: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
1181: KSPConvergedReason reason;
1182: Mat sX,sY;
1183: const PetscScalar *x;
1184: PetscScalar *y;
1185: PetscInt m,N,lda,ldb;
1186: PetscErrorCode ierr;
1189: /* apply preconditioner on each matrix block */
1190: MatGetLocalSize(X,&m,NULL);
1191: MatGetSize(X,NULL,&N);
1192: MatDenseGetLDA(X,&lda);
1193: MatDenseGetLDA(Y,&ldb);
1194: MatDenseGetArrayRead(X,&x);
1195: MatDenseGetArrayWrite(Y,&y);
1196: MatCreateDense(PetscObjectComm((PetscObject)jac->ksp[0]),m,PETSC_DECIDE,PETSC_DECIDE,N,(PetscScalar*)x,&sX);
1197: MatCreateDense(PetscObjectComm((PetscObject)jac->ksp[0]),m,PETSC_DECIDE,PETSC_DECIDE,N,y,&sY);
1198: MatDenseSetLDA(sX,lda);
1199: MatDenseSetLDA(sY,ldb);
1200: PetscLogEventBegin(PC_ApplyOnBlocks,jac->ksp[0],X,Y,0);
1201: KSPMatSolve(jac->ksp[0],sX,sY);
1202: KSPCheckSolve(jac->ksp[0],pc,NULL);
1203: PetscLogEventEnd(PC_ApplyOnBlocks,jac->ksp[0],X,Y,0);
1204: MatDestroy(&sY);
1205: MatDestroy(&sX);
1206: MatDenseRestoreArrayWrite(Y,&y);
1207: MatDenseRestoreArrayRead(X,&x);
1208: KSPGetConvergedReason(jac->ksp[0],&reason);
1209: if (reason == KSP_DIVERGED_PC_FAILED) {
1210: pc->failedreason = PC_SUBPC_ERROR;
1211: }
1212: return(0);
1213: }
1215: static PetscErrorCode PCSetUp_BJacobi_Multiproc(PC pc)
1216: {
1217: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
1218: PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc*)jac->data;
1219: PetscErrorCode ierr;
1220: PetscInt m,n;
1221: MPI_Comm comm,subcomm=0;
1222: const char *prefix;
1223: PetscBool wasSetup = PETSC_TRUE;
1224: VecType vectype;
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: MatGetMultiProcBlock(pc->pmat,subcomm,MAT_INITIAL_MATRIX,&mpjac->submats);
1250: /* create a new PC that processors in each subcomm have copy of */
1251: PetscMalloc1(1,&jac->ksp);
1252: KSPCreate(subcomm,&jac->ksp[0]);
1253: KSPSetErrorIfNotConverged(jac->ksp[0],pc->erroriffailure);
1254: PetscObjectIncrementTabLevel((PetscObject)jac->ksp[0],(PetscObject)pc,1);
1255: PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->ksp[0]);
1256: KSPSetOperators(jac->ksp[0],mpjac->submats,mpjac->submats);
1257: KSPGetPC(jac->ksp[0],&mpjac->pc);
1259: PCGetOptionsPrefix(pc,&prefix);
1260: KSPSetOptionsPrefix(jac->ksp[0],prefix);
1261: KSPAppendOptionsPrefix(jac->ksp[0],"sub_");
1262: KSPGetOptionsPrefix(jac->ksp[0],&prefix);
1263: MatSetOptionsPrefix(mpjac->submats,prefix);
1265: /* create dummy vectors xsub and ysub */
1266: MatGetLocalSize(mpjac->submats,&m,&n);
1267: VecCreateMPIWithArray(subcomm,1,n,PETSC_DECIDE,NULL,&mpjac->xsub);
1268: VecCreateMPIWithArray(subcomm,1,m,PETSC_DECIDE,NULL,&mpjac->ysub);
1269: MatGetVecType(mpjac->submats,&vectype);
1270: VecSetType(mpjac->xsub,vectype);
1271: VecSetType(mpjac->ysub,vectype);
1272: PetscLogObjectParent((PetscObject)pc,(PetscObject)mpjac->xsub);
1273: PetscLogObjectParent((PetscObject)pc,(PetscObject)mpjac->ysub);
1275: pc->ops->reset = PCReset_BJacobi_Multiproc;
1276: pc->ops->destroy = PCDestroy_BJacobi_Multiproc;
1277: pc->ops->apply = PCApply_BJacobi_Multiproc;
1278: pc->ops->matapply= PCMatApply_BJacobi_Multiproc;
1279: } else { /* pc->setupcalled */
1280: subcomm = PetscSubcommChild(mpjac->psubcomm);
1281: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
1282: /* destroy old matrix blocks, then get new matrix blocks */
1283: if (mpjac->submats) {MatDestroy(&mpjac->submats);}
1284: MatGetMultiProcBlock(pc->pmat,subcomm,MAT_INITIAL_MATRIX,&mpjac->submats);
1285: } else {
1286: MatGetMultiProcBlock(pc->pmat,subcomm,MAT_REUSE_MATRIX,&mpjac->submats);
1287: }
1288: KSPSetOperators(jac->ksp[0],mpjac->submats,mpjac->submats);
1289: }
1291: if (!wasSetup && pc->setfromoptionscalled) {
1292: KSPSetFromOptions(jac->ksp[0]);
1293: }
1294: return(0);
1295: }