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;
16: PetscBool hasop;
17: PetscInt N,M,start,i,sum,end;
18: PetscInt bs,i_start=-1,i_end=-1;
19: PetscMPIInt rank,size;
21: MPI_Comm_rank(PetscObjectComm((PetscObject)pc),&rank);
22: MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);
23: MatGetLocalSize(pc->pmat,&M,&N);
24: MatGetBlockSize(pc->pmat,&bs);
26: if (jac->n > 0 && jac->n < size) {
27: PCSetUp_BJacobi_Multiproc(pc);
28: return 0;
29: }
31: /* --------------------------------------------------------------------------
32: Determines the number of blocks assigned to each processor
33: -----------------------------------------------------------------------------*/
35: /* local block count given */
36: if (jac->n_local > 0 && jac->n < 0) {
37: MPIU_Allreduce(&jac->n_local,&jac->n,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)pc));
38: if (jac->l_lens) { /* check that user set these correctly */
39: sum = 0;
40: for (i=0; i<jac->n_local; i++) {
42: sum += jac->l_lens[i];
43: }
45: } else {
46: PetscMalloc1(jac->n_local,&jac->l_lens);
47: for (i=0; i<jac->n_local; i++) jac->l_lens[i] = bs*((M/bs)/jac->n_local + (((M/bs) % jac->n_local) > i));
48: }
49: } else if (jac->n > 0 && jac->n_local < 0) { /* global block count given */
50: /* global blocks given: determine which ones are local */
51: if (jac->g_lens) {
52: /* check if the g_lens is has valid entries */
53: for (i=0; i<jac->n; i++) {
56: }
57: if (size == 1) {
58: jac->n_local = jac->n;
59: PetscMalloc1(jac->n_local,&jac->l_lens);
60: PetscArraycpy(jac->l_lens,jac->g_lens,jac->n_local);
61: /* check that user set these correctly */
62: sum = 0;
63: for (i=0; i<jac->n_local; i++) sum += jac->l_lens[i];
65: } else {
66: MatGetOwnershipRange(pc->pmat,&start,&end);
67: /* loop over blocks determing first one owned by me */
68: sum = 0;
69: for (i=0; i<jac->n+1; i++) {
70: if (sum == start) { i_start = i; goto start_1;}
71: if (i < jac->n) sum += jac->g_lens[i];
72: }
73: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Block sizes used in PCBJacobiSetTotalBlocks()\nare not compatible with parallel matrix layout");
74: start_1:
75: for (i=i_start; i<jac->n+1; i++) {
76: if (sum == end) { i_end = i; goto end_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: end_1:
81: jac->n_local = i_end - i_start;
82: PetscMalloc1(jac->n_local,&jac->l_lens);
83: PetscArraycpy(jac->l_lens,jac->g_lens+i_start,jac->n_local);
84: }
85: } else { /* no global blocks given, determine then using default layout */
86: jac->n_local = jac->n/size + ((jac->n % size) > rank);
87: PetscMalloc1(jac->n_local,&jac->l_lens);
88: for (i=0; i<jac->n_local; i++) {
89: jac->l_lens[i] = ((M/bs)/jac->n_local + (((M/bs) % jac->n_local) > i))*bs;
91: }
92: }
93: } else if (jac->n < 0 && jac->n_local < 0) { /* no blocks given */
94: jac->n = size;
95: jac->n_local = 1;
96: PetscMalloc1(1,&jac->l_lens);
97: jac->l_lens[0] = M;
98: } else { /* jac->n > 0 && jac->n_local > 0 */
99: if (!jac->l_lens) {
100: PetscMalloc1(jac->n_local,&jac->l_lens);
101: for (i=0; i<jac->n_local; i++) jac->l_lens[i] = bs*((M/bs)/jac->n_local + (((M/bs) % jac->n_local) > i));
102: }
103: }
106: /* -------------------------
107: Determines mat and pmat
108: ---------------------------*/
109: MatHasOperation(pc->mat,MATOP_GET_DIAGONAL_BLOCK,&hasop);
110: if (!hasop && size == 1) {
111: mat = pc->mat;
112: pmat = pc->pmat;
113: } else {
114: if (pc->useAmat) {
115: /* use block from Amat matrix, not Pmat for local MatMult() */
116: MatGetDiagonalBlock(pc->mat,&mat);
117: }
118: if (pc->pmat != pc->mat || !pc->useAmat) {
119: MatGetDiagonalBlock(pc->pmat,&pmat);
120: } else pmat = mat;
121: }
123: /* ------
124: Setup code depends on the number of blocks
125: */
126: if (jac->n_local == 1) {
127: PCSetUp_BJacobi_Singleblock(pc,mat,pmat);
128: } else {
129: PCSetUp_BJacobi_Multiblock(pc,mat,pmat);
130: }
131: return 0;
132: }
134: /* Default destroy, if it has never been setup */
135: static PetscErrorCode PCDestroy_BJacobi(PC pc)
136: {
137: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
139: PetscFree(jac->g_lens);
140: PetscFree(jac->l_lens);
141: PetscFree(pc->data);
142: return 0;
143: }
145: static PetscErrorCode PCSetFromOptions_BJacobi(PetscOptionItems *PetscOptionsObject,PC pc)
146: {
147: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
148: PetscInt blocks,i;
149: PetscBool flg;
151: PetscOptionsHead(PetscOptionsObject,"Block Jacobi options");
152: PetscOptionsInt("-pc_bjacobi_blocks","Total number of blocks","PCBJacobiSetTotalBlocks",jac->n,&blocks,&flg);
153: if (flg) PCBJacobiSetTotalBlocks(pc,blocks,NULL);
154: PetscOptionsInt("-pc_bjacobi_local_blocks","Local number of blocks","PCBJacobiSetLocalBlocks",jac->n_local,&blocks,&flg);
155: if (flg) PCBJacobiSetLocalBlocks(pc,blocks,NULL);
156: if (jac->ksp) {
157: /* The sub-KSP has already been set up (e.g., PCSetUp_BJacobi_Singleblock), but KSPSetFromOptions was not called
158: * unless we had already been called. */
159: for (i=0; i<jac->n_local; i++) {
160: KSPSetFromOptions(jac->ksp[i]);
161: }
162: }
163: PetscOptionsTail();
164: return 0;
165: }
167: #include <petscdraw.h>
168: static PetscErrorCode PCView_BJacobi(PC pc,PetscViewer viewer)
169: {
170: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
171: PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc*)jac->data;
172: PetscErrorCode ierr;
173: PetscMPIInt rank;
174: PetscInt i;
175: PetscBool iascii,isstring,isdraw;
176: PetscViewer sviewer;
177: PetscViewerFormat format;
178: const char *prefix;
180: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
181: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
182: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
183: if (iascii) {
184: if (pc->useAmat) {
185: PetscViewerASCIIPrintf(viewer," using Amat local matrix, number of blocks = %D\n",jac->n);
186: }
187: PetscViewerASCIIPrintf(viewer," number of blocks = %D\n",jac->n);
188: MPI_Comm_rank(PetscObjectComm((PetscObject)pc),&rank);
189: PetscViewerGetFormat(viewer,&format);
190: if (format != PETSC_VIEWER_ASCII_INFO_DETAIL) {
191: PetscViewerASCIIPrintf(viewer," Local solver information for first block is in the following KSP and PC objects on rank 0:\n");
192: PCGetOptionsPrefix(pc,&prefix);
193: PetscViewerASCIIPrintf(viewer," Use -%sksp_view ::ascii_info_detail to display information for all blocks\n",prefix?prefix:"");
194: if (jac->ksp && !jac->psubcomm) {
195: PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
196: if (rank == 0) {
197: PetscViewerASCIIPushTab(viewer);
198: KSPView(jac->ksp[0],sviewer);
199: PetscViewerASCIIPopTab(viewer);
200: }
201: PetscViewerFlush(sviewer);
202: PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
203: PetscViewerFlush(viewer);
204: /* extra call needed because of the two calls to PetscViewerASCIIPushSynchronized() in PetscViewerGetSubViewer() */
205: PetscViewerASCIIPopSynchronized(viewer);
206: } else if (mpjac && jac->ksp && mpjac->psubcomm) {
207: PetscViewerGetSubViewer(viewer,mpjac->psubcomm->child,&sviewer);
208: if (!mpjac->psubcomm->color) {
209: PetscViewerASCIIPushTab(viewer);
210: KSPView(*(jac->ksp),sviewer);
211: PetscViewerASCIIPopTab(viewer);
212: }
213: PetscViewerFlush(sviewer);
214: PetscViewerRestoreSubViewer(viewer,mpjac->psubcomm->child,&sviewer);
215: PetscViewerFlush(viewer);
216: /* extra call needed because of the two calls to PetscViewerASCIIPushSynchronized() in PetscViewerGetSubViewer() */
217: PetscViewerASCIIPopSynchronized(viewer);
218: } else {
219: PetscViewerFlush(viewer);
220: }
221: } else {
222: PetscInt n_global;
223: MPIU_Allreduce(&jac->n_local,&n_global,1,MPIU_INT,MPI_MAX,PetscObjectComm((PetscObject)pc));
224: PetscViewerASCIIPushSynchronized(viewer);
225: PetscViewerASCIIPrintf(viewer," Local solver information for each block is in the following KSP and PC objects:\n");
226: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] number of local blocks = %D, first local block number = %D\n",
227: rank,jac->n_local,jac->first_local);
228: PetscViewerASCIIPushTab(viewer);
229: PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
230: for (i=0; i<jac->n_local; i++) {
231: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] local block number %D\n",rank,i);
232: KSPView(jac->ksp[i],sviewer);
233: PetscViewerASCIISynchronizedPrintf(viewer,"- - - - - - - - - - - - - - - - - -\n");
234: }
235: PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
236: PetscViewerASCIIPopTab(viewer);
237: PetscViewerFlush(viewer);
238: PetscViewerASCIIPopSynchronized(viewer);
239: }
240: } else if (isstring) {
241: PetscViewerStringSPrintf(viewer," blks=%D",jac->n);
242: PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
243: if (jac->ksp) KSPView(jac->ksp[0],sviewer);
244: PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
245: } else if (isdraw) {
246: PetscDraw draw;
247: char str[25];
248: PetscReal x,y,bottom,h;
250: PetscViewerDrawGetDraw(viewer,0,&draw);
251: PetscDrawGetCurrentPoint(draw,&x,&y);
252: PetscSNPrintf(str,25,"Number blocks %D",jac->n);
253: PetscDrawStringBoxed(draw,x,y,PETSC_DRAW_RED,PETSC_DRAW_BLACK,str,NULL,&h);
254: bottom = y - h;
255: PetscDrawPushCurrentPoint(draw,x,bottom);
256: /* warning the communicator on viewer is different then on ksp in parallel */
257: if (jac->ksp) KSPView(jac->ksp[0],viewer);
258: PetscDrawPopCurrentPoint(draw);
259: }
260: return 0;
261: }
263: /* -------------------------------------------------------------------------------------*/
265: static PetscErrorCode PCBJacobiGetSubKSP_BJacobi(PC pc,PetscInt *n_local,PetscInt *first_local,KSP **ksp)
266: {
267: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
271: if (n_local) *n_local = jac->n_local;
272: if (first_local) *first_local = jac->first_local;
273: if (ksp) *ksp = jac->ksp;
274: return 0;
275: }
277: static PetscErrorCode PCBJacobiSetTotalBlocks_BJacobi(PC pc,PetscInt blocks,PetscInt *lens)
278: {
279: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
282: jac->n = blocks;
283: if (!lens) jac->g_lens = NULL;
284: else {
285: PetscMalloc1(blocks,&jac->g_lens);
286: PetscLogObjectMemory((PetscObject)pc,blocks*sizeof(PetscInt));
287: PetscArraycpy(jac->g_lens,lens,blocks);
288: }
289: return 0;
290: }
292: static PetscErrorCode PCBJacobiGetTotalBlocks_BJacobi(PC pc, PetscInt *blocks, const PetscInt *lens[])
293: {
294: PC_BJacobi *jac = (PC_BJacobi*) pc->data;
296: *blocks = jac->n;
297: if (lens) *lens = jac->g_lens;
298: return 0;
299: }
301: static PetscErrorCode PCBJacobiSetLocalBlocks_BJacobi(PC pc,PetscInt blocks,const PetscInt lens[])
302: {
303: PC_BJacobi *jac;
305: jac = (PC_BJacobi*)pc->data;
307: jac->n_local = blocks;
308: if (!lens) jac->l_lens = NULL;
309: else {
310: PetscMalloc1(blocks,&jac->l_lens);
311: PetscLogObjectMemory((PetscObject)pc,blocks*sizeof(PetscInt));
312: PetscArraycpy(jac->l_lens,lens,blocks);
313: }
314: return 0;
315: }
317: static PetscErrorCode PCBJacobiGetLocalBlocks_BJacobi(PC pc, PetscInt *blocks, const PetscInt *lens[])
318: {
319: PC_BJacobi *jac = (PC_BJacobi*) pc->data;
321: *blocks = jac->n_local;
322: if (lens) *lens = jac->l_lens;
323: return 0;
324: }
326: /* -------------------------------------------------------------------------------------*/
328: /*@C
329: PCBJacobiGetSubKSP - Gets the local KSP contexts for all blocks on
330: this processor.
332: Not Collective
334: Input Parameter:
335: . pc - the preconditioner context
337: Output Parameters:
338: + n_local - the number of blocks on this processor, or NULL
339: . first_local - the global number of the first block on this processor, or NULL
340: - ksp - the array of KSP contexts
342: Notes:
343: After PCBJacobiGetSubKSP() the array of KSP contexts is not to be freed.
345: Currently for some matrix implementations only 1 block per processor
346: is supported.
348: You must call KSPSetUp() or PCSetUp() before calling PCBJacobiGetSubKSP().
350: Fortran Usage: You must pass in a KSP array that is large enough to contain all the local KSPs.
351: You can call PCBJacobiGetSubKSP(pc,nlocal,firstlocal,PETSC_NULL_KSP,ierr) to determine how large the
352: KSP array must be.
354: Level: advanced
356: .seealso: PCASMGetSubKSP()
357: @*/
358: PetscErrorCode PCBJacobiGetSubKSP(PC pc,PetscInt *n_local,PetscInt *first_local,KSP *ksp[])
359: {
361: PetscUseMethod(pc,"PCBJacobiGetSubKSP_C",(PC,PetscInt*,PetscInt*,KSP **),(pc,n_local,first_local,ksp));
362: return 0;
363: }
365: /*@
366: PCBJacobiSetTotalBlocks - Sets the global number of blocks for the block
367: Jacobi preconditioner.
369: Collective on PC
371: Input Parameters:
372: + pc - the preconditioner context
373: . blocks - the number of blocks
374: - lens - [optional] integer array containing the size of each block
376: Options Database Key:
377: . -pc_bjacobi_blocks <blocks> - Sets the number of global blocks
379: Notes:
380: Currently only a limited number of blocking configurations are supported.
381: All processors sharing the PC must call this routine with the same data.
383: Level: intermediate
385: .seealso: PCSetUseAmat(), PCBJacobiSetLocalBlocks()
386: @*/
387: PetscErrorCode PCBJacobiSetTotalBlocks(PC pc,PetscInt blocks,const PetscInt lens[])
388: {
391: PetscTryMethod(pc,"PCBJacobiSetTotalBlocks_C",(PC,PetscInt,const PetscInt[]),(pc,blocks,lens));
392: return 0;
393: }
395: /*@C
396: PCBJacobiGetTotalBlocks - Gets the global number of blocks for the block
397: Jacobi preconditioner.
399: Not Collective
401: Input Parameter:
402: . pc - the preconditioner context
404: Output parameters:
405: + blocks - the number of blocks
406: - lens - integer array containing the size of each block
408: Level: intermediate
410: .seealso: PCSetUseAmat(), PCBJacobiGetLocalBlocks()
411: @*/
412: PetscErrorCode PCBJacobiGetTotalBlocks(PC pc, PetscInt *blocks, const PetscInt *lens[])
413: {
416: PetscUseMethod(pc,"PCBJacobiGetTotalBlocks_C",(PC,PetscInt*, const PetscInt *[]),(pc,blocks,lens));
417: return 0;
418: }
420: /*@
421: PCBJacobiSetLocalBlocks - Sets the local number of blocks for the block
422: Jacobi preconditioner.
424: Not Collective
426: Input Parameters:
427: + pc - the preconditioner context
428: . blocks - the number of blocks
429: - lens - [optional] integer array containing size of each block
431: Options Database Key:
432: . -pc_bjacobi_local_blocks <blocks> - Sets the number of local blocks
434: Note:
435: Currently only a limited number of blocking configurations are supported.
437: Level: intermediate
439: .seealso: PCSetUseAmat(), PCBJacobiSetTotalBlocks()
440: @*/
441: PetscErrorCode PCBJacobiSetLocalBlocks(PC pc,PetscInt blocks,const PetscInt lens[])
442: {
445: PetscTryMethod(pc,"PCBJacobiSetLocalBlocks_C",(PC,PetscInt,const PetscInt []),(pc,blocks,lens));
446: return 0;
447: }
449: /*@C
450: PCBJacobiGetLocalBlocks - Gets the local number of blocks for the block
451: Jacobi preconditioner.
453: Not Collective
455: Input Parameters:
456: + pc - the preconditioner context
457: . blocks - the number of blocks
458: - lens - [optional] integer array containing size of each block
460: Note:
461: Currently only a limited number of blocking configurations are supported.
463: Level: intermediate
465: .seealso: PCSetUseAmat(), PCBJacobiGetTotalBlocks()
466: @*/
467: PetscErrorCode PCBJacobiGetLocalBlocks(PC pc, PetscInt *blocks, const PetscInt *lens[])
468: {
471: PetscUseMethod(pc,"PCBJacobiGetLocalBlocks_C",(PC,PetscInt*, const PetscInt *[]),(pc,blocks,lens));
472: return 0;
473: }
475: /* -----------------------------------------------------------------------------------*/
477: /*MC
478: PCBJACOBI - Use block Jacobi preconditioning, each block is (approximately) solved with
479: its own KSP object.
481: Options Database Keys:
482: + -pc_use_amat - use Amat to apply block of operator in inner Krylov method
483: - -pc_bjacobi_blocks <n> - use n total blocks
485: Notes:
486: Each processor can have one or more blocks, or a single block can be shared by several processes. Defaults to one block per processor.
488: To set options on the solvers for each block append -sub_ to all the KSP, KSP, and PC
489: options database keys. For example, -sub_pc_type ilu -sub_pc_factor_levels 1 -sub_ksp_type preonly
491: To set the options on the solvers separate for each block call PCBJacobiGetSubKSP()
492: and set the options directly on the resulting KSP object (you can access its PC
493: KSPGetPC())
495: For GPU-based vectors (CUDA, ViennaCL) it is recommended to use exactly one block per MPI process for best
496: performance. Different block partitioning may lead to additional data transfers
497: between host and GPU that lead to degraded performance.
499: The options prefix for each block is sub_, for example -sub_pc_type lu.
501: When multiple processes share a single block, each block encompasses exactly all the unknowns owned its set of processes.
503: See PCJACOBI for point Jacobi preconditioning, PCVPBJACOBI for variable size point block Jacobi and PCPBJACOBI for large blocks
505: Level: beginner
507: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
508: PCASM, PCSetUseAmat(), PCGetUseAmat(), PCBJacobiGetSubKSP(), PCBJacobiSetTotalBlocks(),
509: PCBJacobiSetLocalBlocks(), PCSetModifySubMatrices(), PCJACOBI, PCVPBJACOBI, PCPBJACOBI
510: M*/
512: PETSC_EXTERN PetscErrorCode PCCreate_BJacobi(PC pc)
513: {
514: PetscMPIInt rank;
515: PC_BJacobi *jac;
517: PetscNewLog(pc,&jac);
518: MPI_Comm_rank(PetscObjectComm((PetscObject)pc),&rank);
520: pc->ops->apply = NULL;
521: pc->ops->matapply = NULL;
522: pc->ops->applytranspose = NULL;
523: pc->ops->setup = PCSetUp_BJacobi;
524: pc->ops->destroy = PCDestroy_BJacobi;
525: pc->ops->setfromoptions = PCSetFromOptions_BJacobi;
526: pc->ops->view = PCView_BJacobi;
527: pc->ops->applyrichardson = NULL;
529: pc->data = (void*)jac;
530: jac->n = -1;
531: jac->n_local = -1;
532: jac->first_local = rank;
533: jac->ksp = NULL;
534: jac->g_lens = NULL;
535: jac->l_lens = NULL;
536: jac->psubcomm = NULL;
538: PetscObjectComposeFunction((PetscObject)pc,"PCBJacobiGetSubKSP_C",PCBJacobiGetSubKSP_BJacobi);
539: PetscObjectComposeFunction((PetscObject)pc,"PCBJacobiSetTotalBlocks_C",PCBJacobiSetTotalBlocks_BJacobi);
540: PetscObjectComposeFunction((PetscObject)pc,"PCBJacobiGetTotalBlocks_C",PCBJacobiGetTotalBlocks_BJacobi);
541: PetscObjectComposeFunction((PetscObject)pc,"PCBJacobiSetLocalBlocks_C",PCBJacobiSetLocalBlocks_BJacobi);
542: PetscObjectComposeFunction((PetscObject)pc,"PCBJacobiGetLocalBlocks_C",PCBJacobiGetLocalBlocks_BJacobi);
543: return 0;
544: }
546: /* --------------------------------------------------------------------------------------------*/
547: /*
548: These are for a single block per processor; works for AIJ, BAIJ; Seq and MPI
549: */
550: static PetscErrorCode PCReset_BJacobi_Singleblock(PC pc)
551: {
552: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
553: PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock*)jac->data;
555: KSPReset(jac->ksp[0]);
556: VecDestroy(&bjac->x);
557: VecDestroy(&bjac->y);
558: return 0;
559: }
561: static PetscErrorCode PCDestroy_BJacobi_Singleblock(PC pc)
562: {
563: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
564: PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock*)jac->data;
566: PCReset_BJacobi_Singleblock(pc);
567: KSPDestroy(&jac->ksp[0]);
568: PetscFree(jac->ksp);
569: PetscFree(jac->l_lens);
570: PetscFree(jac->g_lens);
571: PetscFree(bjac);
572: PetscFree(pc->data);
573: return 0;
574: }
576: static PetscErrorCode PCSetUpOnBlocks_BJacobi_Singleblock(PC pc)
577: {
578: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
579: KSP subksp = jac->ksp[0];
580: KSPConvergedReason reason;
582: KSPSetUp(subksp);
583: KSPGetConvergedReason(subksp,&reason);
584: if (reason == KSP_DIVERGED_PC_FAILED) {
585: pc->failedreason = PC_SUBPC_ERROR;
586: }
587: return 0;
588: }
590: static PetscErrorCode PCApply_BJacobi_Singleblock(PC pc,Vec x,Vec y)
591: {
592: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
593: PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock*)jac->data;
595: VecGetLocalVectorRead(x, bjac->x);
596: VecGetLocalVector(y, bjac->y);
597: /* Since the inner KSP matrix may point directly to the diagonal block of an MPI matrix the inner
598: matrix may change even if the outer KSP/PC has not updated the preconditioner, this will trigger a rebuild
599: of the inner preconditioner automatically unless we pass down the outer preconditioners reuse flag.*/
600: KSPSetReusePreconditioner(jac->ksp[0],pc->reusepreconditioner);
601: KSPSolve(jac->ksp[0],bjac->x,bjac->y);
602: KSPCheckSolve(jac->ksp[0],pc,bjac->y);
603: VecRestoreLocalVectorRead(x, bjac->x);
604: VecRestoreLocalVector(y, bjac->y);
605: return 0;
606: }
608: static PetscErrorCode PCMatApply_BJacobi_Singleblock(PC pc,Mat X,Mat Y)
609: {
610: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
611: Mat sX,sY;
613: /* Since the inner KSP matrix may point directly to the diagonal block of an MPI matrix the inner
614: matrix may change even if the outer KSP/PC has not updated the preconditioner, this will trigger a rebuild
615: of the inner preconditioner automatically unless we pass down the outer preconditioners reuse flag.*/
616: KSPSetReusePreconditioner(jac->ksp[0],pc->reusepreconditioner);
617: MatDenseGetLocalMatrix(X,&sX);
618: MatDenseGetLocalMatrix(Y,&sY);
619: KSPMatSolve(jac->ksp[0],sX,sY);
620: return 0;
621: }
623: static PetscErrorCode PCApplySymmetricLeft_BJacobi_Singleblock(PC pc,Vec x,Vec y)
624: {
625: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
626: PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock*)jac->data;
627: PetscScalar *y_array;
628: const PetscScalar *x_array;
629: PC subpc;
631: /*
632: The VecPlaceArray() is to avoid having to copy the
633: y vector into the bjac->x vector. The reason for
634: the bjac->x vector is that we need a sequential vector
635: for the sequential solve.
636: */
637: VecGetArrayRead(x,&x_array);
638: VecGetArray(y,&y_array);
639: VecPlaceArray(bjac->x,x_array);
640: VecPlaceArray(bjac->y,y_array);
641: /* apply the symmetric left portion of the inner PC operator */
642: /* note this by-passes the inner KSP and its options completely */
643: KSPGetPC(jac->ksp[0],&subpc);
644: PCApplySymmetricLeft(subpc,bjac->x,bjac->y);
645: VecResetArray(bjac->x);
646: VecResetArray(bjac->y);
647: VecRestoreArrayRead(x,&x_array);
648: VecRestoreArray(y,&y_array);
649: return 0;
650: }
652: static PetscErrorCode PCApplySymmetricRight_BJacobi_Singleblock(PC pc,Vec x,Vec y)
653: {
654: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
655: PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock*)jac->data;
656: PetscScalar *y_array;
657: const PetscScalar *x_array;
658: PC subpc;
660: /*
661: The VecPlaceArray() is to avoid having to copy the
662: y vector into the bjac->x vector. The reason for
663: the bjac->x vector is that we need a sequential vector
664: for the sequential solve.
665: */
666: VecGetArrayRead(x,&x_array);
667: VecGetArray(y,&y_array);
668: VecPlaceArray(bjac->x,x_array);
669: VecPlaceArray(bjac->y,y_array);
671: /* apply the symmetric right portion of the inner PC operator */
672: /* note this by-passes the inner KSP and its options completely */
674: KSPGetPC(jac->ksp[0],&subpc);
675: PCApplySymmetricRight(subpc,bjac->x,bjac->y);
677: VecRestoreArrayRead(x,&x_array);
678: VecRestoreArray(y,&y_array);
679: return 0;
680: }
682: static PetscErrorCode PCApplyTranspose_BJacobi_Singleblock(PC pc,Vec x,Vec y)
683: {
684: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
685: PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock*)jac->data;
686: PetscScalar *y_array;
687: const PetscScalar *x_array;
689: /*
690: The VecPlaceArray() is to avoid having to copy the
691: y vector into the bjac->x vector. The reason for
692: the bjac->x vector is that we need a sequential vector
693: for the sequential solve.
694: */
695: VecGetArrayRead(x,&x_array);
696: VecGetArray(y,&y_array);
697: VecPlaceArray(bjac->x,x_array);
698: VecPlaceArray(bjac->y,y_array);
699: KSPSolveTranspose(jac->ksp[0],bjac->x,bjac->y);
700: KSPCheckSolve(jac->ksp[0],pc,bjac->y);
701: VecResetArray(bjac->x);
702: VecResetArray(bjac->y);
703: VecRestoreArrayRead(x,&x_array);
704: VecRestoreArray(y,&y_array);
705: return 0;
706: }
708: static PetscErrorCode PCSetUp_BJacobi_Singleblock(PC pc,Mat mat,Mat pmat)
709: {
710: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
711: PetscInt m;
712: KSP ksp;
713: PC_BJacobi_Singleblock *bjac;
714: PetscBool wasSetup = PETSC_TRUE;
715: VecType vectype;
716: const char *prefix;
718: if (!pc->setupcalled) {
719: if (!jac->ksp) {
720: wasSetup = PETSC_FALSE;
722: KSPCreate(PETSC_COMM_SELF,&ksp);
723: KSPSetErrorIfNotConverged(ksp,pc->erroriffailure);
724: PetscObjectIncrementTabLevel((PetscObject)ksp,(PetscObject)pc,1);
725: PetscLogObjectParent((PetscObject)pc,(PetscObject)ksp);
726: KSPSetType(ksp,KSPPREONLY);
727: PCGetOptionsPrefix(pc,&prefix);
728: KSPSetOptionsPrefix(ksp,prefix);
729: KSPAppendOptionsPrefix(ksp,"sub_");
731: pc->ops->reset = PCReset_BJacobi_Singleblock;
732: pc->ops->destroy = PCDestroy_BJacobi_Singleblock;
733: pc->ops->apply = PCApply_BJacobi_Singleblock;
734: pc->ops->matapply = PCMatApply_BJacobi_Singleblock;
735: pc->ops->applysymmetricleft = PCApplySymmetricLeft_BJacobi_Singleblock;
736: pc->ops->applysymmetricright = PCApplySymmetricRight_BJacobi_Singleblock;
737: pc->ops->applytranspose = PCApplyTranspose_BJacobi_Singleblock;
738: pc->ops->setuponblocks = PCSetUpOnBlocks_BJacobi_Singleblock;
740: PetscMalloc1(1,&jac->ksp);
741: jac->ksp[0] = ksp;
743: PetscNewLog(pc,&bjac);
744: jac->data = (void*)bjac;
745: } else {
746: ksp = jac->ksp[0];
747: bjac = (PC_BJacobi_Singleblock*)jac->data;
748: }
750: /*
751: The reason we need to generate these vectors is to serve
752: as the right-hand side and solution vector for the solve on the
753: block. We do not need to allocate space for the vectors since
754: that is provided via VecPlaceArray() just before the call to
755: KSPSolve() on the block.
756: */
757: MatGetSize(pmat,&m,&m);
758: VecCreateSeqWithArray(PETSC_COMM_SELF,1,m,NULL,&bjac->x);
759: VecCreateSeqWithArray(PETSC_COMM_SELF,1,m,NULL,&bjac->y);
760: MatGetVecType(pmat,&vectype);
761: VecSetType(bjac->x,vectype);
762: VecSetType(bjac->y,vectype);
763: PetscLogObjectParent((PetscObject)pc,(PetscObject)bjac->x);
764: PetscLogObjectParent((PetscObject)pc,(PetscObject)bjac->y);
765: } else {
766: ksp = jac->ksp[0];
767: bjac = (PC_BJacobi_Singleblock*)jac->data;
768: }
769: KSPGetOptionsPrefix(ksp,&prefix);
770: if (pc->useAmat) {
771: KSPSetOperators(ksp,mat,pmat);
772: MatSetOptionsPrefix(mat,prefix);
773: } else {
774: KSPSetOperators(ksp,pmat,pmat);
775: }
776: MatSetOptionsPrefix(pmat,prefix);
777: if (!wasSetup && pc->setfromoptionscalled) {
778: /* If PCSetFromOptions_BJacobi is called later, KSPSetFromOptions will be called at that time. */
779: KSPSetFromOptions(ksp);
780: }
781: return 0;
782: }
784: /* ---------------------------------------------------------------------------------------------*/
785: static PetscErrorCode PCReset_BJacobi_Multiblock(PC pc)
786: {
787: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
788: PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock*)jac->data;
789: PetscInt i;
791: if (bjac && bjac->pmat) {
792: MatDestroyMatrices(jac->n_local,&bjac->pmat);
793: if (pc->useAmat) {
794: MatDestroyMatrices(jac->n_local,&bjac->mat);
795: }
796: }
798: for (i=0; i<jac->n_local; i++) {
799: KSPReset(jac->ksp[i]);
800: if (bjac && bjac->x) {
801: VecDestroy(&bjac->x[i]);
802: VecDestroy(&bjac->y[i]);
803: ISDestroy(&bjac->is[i]);
804: }
805: }
806: PetscFree(jac->l_lens);
807: PetscFree(jac->g_lens);
808: return 0;
809: }
811: static PetscErrorCode PCDestroy_BJacobi_Multiblock(PC pc)
812: {
813: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
814: PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock*)jac->data;
815: PetscInt i;
817: PCReset_BJacobi_Multiblock(pc);
818: if (bjac) {
819: PetscFree2(bjac->x,bjac->y);
820: PetscFree(bjac->starts);
821: PetscFree(bjac->is);
822: }
823: PetscFree(jac->data);
824: for (i=0; i<jac->n_local; i++) {
825: KSPDestroy(&jac->ksp[i]);
826: }
827: PetscFree(jac->ksp);
828: PetscFree(pc->data);
829: return 0;
830: }
832: static PetscErrorCode PCSetUpOnBlocks_BJacobi_Multiblock(PC pc)
833: {
834: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
835: PetscInt i,n_local = jac->n_local;
836: KSPConvergedReason reason;
838: for (i=0; i<n_local; i++) {
839: KSPSetUp(jac->ksp[i]);
840: KSPGetConvergedReason(jac->ksp[i],&reason);
841: if (reason == KSP_DIVERGED_PC_FAILED) {
842: pc->failedreason = PC_SUBPC_ERROR;
843: }
844: }
845: return 0;
846: }
848: /*
849: Preconditioner for block Jacobi
850: */
851: static PetscErrorCode PCApply_BJacobi_Multiblock(PC pc,Vec x,Vec y)
852: {
853: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
854: PetscInt i,n_local = jac->n_local;
855: PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock*)jac->data;
856: PetscScalar *yin;
857: const PetscScalar *xin;
859: VecGetArrayRead(x,&xin);
860: VecGetArray(y,&yin);
861: for (i=0; i<n_local; i++) {
862: /*
863: To avoid copying the subvector from x into a workspace we instead
864: make the workspace vector array point to the subpart of the array of
865: the global vector.
866: */
867: VecPlaceArray(bjac->x[i],xin+bjac->starts[i]);
868: VecPlaceArray(bjac->y[i],yin+bjac->starts[i]);
870: PetscLogEventBegin(PC_ApplyOnBlocks,jac->ksp[i],bjac->x[i],bjac->y[i],0);
871: KSPSolve(jac->ksp[i],bjac->x[i],bjac->y[i]);
872: KSPCheckSolve(jac->ksp[i],pc,bjac->y[i]);
873: PetscLogEventEnd(PC_ApplyOnBlocks,jac->ksp[i],bjac->x[i],bjac->y[i],0);
875: VecResetArray(bjac->x[i]);
876: VecResetArray(bjac->y[i]);
877: }
878: VecRestoreArrayRead(x,&xin);
879: VecRestoreArray(y,&yin);
880: return 0;
881: }
883: /*
884: Preconditioner for block Jacobi
885: */
886: static PetscErrorCode PCApplyTranspose_BJacobi_Multiblock(PC pc,Vec x,Vec y)
887: {
888: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
889: PetscInt i,n_local = jac->n_local;
890: PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock*)jac->data;
891: PetscScalar *yin;
892: const PetscScalar *xin;
894: VecGetArrayRead(x,&xin);
895: VecGetArray(y,&yin);
896: for (i=0; i<n_local; i++) {
897: /*
898: To avoid copying the subvector from x into a workspace we instead
899: make the workspace vector array point to the subpart of the array of
900: the global vector.
901: */
902: VecPlaceArray(bjac->x[i],xin+bjac->starts[i]);
903: VecPlaceArray(bjac->y[i],yin+bjac->starts[i]);
905: PetscLogEventBegin(PC_ApplyTransposeOnBlocks,jac->ksp[i],bjac->x[i],bjac->y[i],0);
906: KSPSolveTranspose(jac->ksp[i],bjac->x[i],bjac->y[i]);
907: KSPCheckSolve(jac->ksp[i],pc,bjac->y[i]);
908: PetscLogEventEnd(PC_ApplyTransposeOnBlocks,jac->ksp[i],bjac->x[i],bjac->y[i],0);
910: VecResetArray(bjac->x[i]);
911: VecResetArray(bjac->y[i]);
912: }
913: VecRestoreArrayRead(x,&xin);
914: VecRestoreArray(y,&yin);
915: return 0;
916: }
918: static PetscErrorCode PCSetUp_BJacobi_Multiblock(PC pc,Mat mat,Mat pmat)
919: {
920: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
921: PetscInt m,n_local,N,M,start,i;
922: const char *prefix;
923: KSP ksp;
924: Vec x,y;
925: PC_BJacobi_Multiblock *bjac = (PC_BJacobi_Multiblock*)jac->data;
926: PC subpc;
927: IS is;
928: MatReuse scall;
929: VecType vectype;
931: MatGetLocalSize(pc->pmat,&M,&N);
933: n_local = jac->n_local;
935: if (pc->useAmat) {
936: PetscBool same;
937: PetscObjectTypeCompare((PetscObject)mat,((PetscObject)pmat)->type_name,&same);
939: }
941: if (!pc->setupcalled) {
942: scall = MAT_INITIAL_MATRIX;
944: if (!jac->ksp) {
945: pc->ops->reset = PCReset_BJacobi_Multiblock;
946: pc->ops->destroy = PCDestroy_BJacobi_Multiblock;
947: pc->ops->apply = PCApply_BJacobi_Multiblock;
948: pc->ops->matapply = NULL;
949: pc->ops->applytranspose= PCApplyTranspose_BJacobi_Multiblock;
950: pc->ops->setuponblocks = PCSetUpOnBlocks_BJacobi_Multiblock;
952: PetscNewLog(pc,&bjac);
953: PetscMalloc1(n_local,&jac->ksp);
954: PetscLogObjectMemory((PetscObject)pc,sizeof(n_local*sizeof(KSP)));
955: PetscMalloc2(n_local,&bjac->x,n_local,&bjac->y);
956: PetscMalloc1(n_local,&bjac->starts);
957: PetscLogObjectMemory((PetscObject)pc,sizeof(n_local*sizeof(PetscScalar)));
959: jac->data = (void*)bjac;
960: PetscMalloc1(n_local,&bjac->is);
961: PetscLogObjectMemory((PetscObject)pc,sizeof(n_local*sizeof(IS)));
963: for (i=0; i<n_local; i++) {
964: KSPCreate(PETSC_COMM_SELF,&ksp);
965: KSPSetErrorIfNotConverged(ksp,pc->erroriffailure);
966: PetscObjectIncrementTabLevel((PetscObject)ksp,(PetscObject)pc,1);
967: PetscLogObjectParent((PetscObject)pc,(PetscObject)ksp);
968: KSPSetType(ksp,KSPPREONLY);
969: KSPGetPC(ksp,&subpc);
970: PCGetOptionsPrefix(pc,&prefix);
971: KSPSetOptionsPrefix(ksp,prefix);
972: KSPAppendOptionsPrefix(ksp,"sub_");
974: jac->ksp[i] = ksp;
975: }
976: } else {
977: bjac = (PC_BJacobi_Multiblock*)jac->data;
978: }
980: start = 0;
981: MatGetVecType(pmat,&vectype);
982: for (i=0; i<n_local; i++) {
983: m = jac->l_lens[i];
984: /*
985: The reason we need to generate these vectors is to serve
986: as the right-hand side and solution vector for the solve on the
987: block. We do not need to allocate space for the vectors since
988: that is provided via VecPlaceArray() just before the call to
989: KSPSolve() on the block.
991: */
992: VecCreateSeq(PETSC_COMM_SELF,m,&x);
993: VecCreateSeqWithArray(PETSC_COMM_SELF,1,m,NULL,&y);
994: VecSetType(x,vectype);
995: VecSetType(y,vectype);
996: PetscLogObjectParent((PetscObject)pc,(PetscObject)x);
997: PetscLogObjectParent((PetscObject)pc,(PetscObject)y);
999: bjac->x[i] = x;
1000: bjac->y[i] = y;
1001: bjac->starts[i] = start;
1003: ISCreateStride(PETSC_COMM_SELF,m,start,1,&is);
1004: bjac->is[i] = is;
1005: PetscLogObjectParent((PetscObject)pc,(PetscObject)is);
1007: start += m;
1008: }
1009: } else {
1010: bjac = (PC_BJacobi_Multiblock*)jac->data;
1011: /*
1012: Destroy the blocks from the previous iteration
1013: */
1014: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
1015: MatDestroyMatrices(n_local,&bjac->pmat);
1016: if (pc->useAmat) {
1017: MatDestroyMatrices(n_local,&bjac->mat);
1018: }
1019: scall = MAT_INITIAL_MATRIX;
1020: } else scall = MAT_REUSE_MATRIX;
1021: }
1023: MatCreateSubMatrices(pmat,n_local,bjac->is,bjac->is,scall,&bjac->pmat);
1024: if (pc->useAmat) {
1025: MatCreateSubMatrices(mat,n_local,bjac->is,bjac->is,scall,&bjac->mat);
1026: }
1027: /* Return control to the user so that the submatrices can be modified (e.g., to apply
1028: different boundary conditions for the submatrices than for the global problem) */
1029: PCModifySubMatrices(pc,n_local,bjac->is,bjac->is,bjac->pmat,pc->modifysubmatricesP);
1031: for (i=0; i<n_local; i++) {
1032: PetscLogObjectParent((PetscObject)pc,(PetscObject)bjac->pmat[i]);
1033: KSPGetOptionsPrefix(jac->ksp[i],&prefix);
1034: if (pc->useAmat) {
1035: PetscLogObjectParent((PetscObject)pc,(PetscObject)bjac->mat[i]);
1036: KSPSetOperators(jac->ksp[i],bjac->mat[i],bjac->pmat[i]);
1037: MatSetOptionsPrefix(bjac->mat[i],prefix);
1038: } else {
1039: KSPSetOperators(jac->ksp[i],bjac->pmat[i],bjac->pmat[i]);
1040: }
1041: MatSetOptionsPrefix(bjac->pmat[i],prefix);
1042: if (pc->setfromoptionscalled) {
1043: KSPSetFromOptions(jac->ksp[i]);
1044: }
1045: }
1046: return 0;
1047: }
1049: /* ---------------------------------------------------------------------------------------------*/
1050: /*
1051: These are for a single block with multiple processes
1052: */
1053: static PetscErrorCode PCSetUpOnBlocks_BJacobi_Multiproc(PC pc)
1054: {
1055: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
1056: KSP subksp = jac->ksp[0];
1057: KSPConvergedReason reason;
1059: KSPSetUp(subksp);
1060: KSPGetConvergedReason(subksp,&reason);
1061: if (reason == KSP_DIVERGED_PC_FAILED) {
1062: pc->failedreason = PC_SUBPC_ERROR;
1063: }
1064: return 0;
1065: }
1067: static PetscErrorCode PCReset_BJacobi_Multiproc(PC pc)
1068: {
1069: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
1070: PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc*)jac->data;
1072: VecDestroy(&mpjac->ysub);
1073: VecDestroy(&mpjac->xsub);
1074: MatDestroy(&mpjac->submats);
1075: if (jac->ksp) KSPReset(jac->ksp[0]);
1076: return 0;
1077: }
1079: static PetscErrorCode PCDestroy_BJacobi_Multiproc(PC pc)
1080: {
1081: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
1082: PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc*)jac->data;
1084: PCReset_BJacobi_Multiproc(pc);
1085: KSPDestroy(&jac->ksp[0]);
1086: PetscFree(jac->ksp);
1087: PetscSubcommDestroy(&mpjac->psubcomm);
1089: PetscFree(mpjac);
1090: PetscFree(pc->data);
1091: return 0;
1092: }
1094: static PetscErrorCode PCApply_BJacobi_Multiproc(PC pc,Vec x,Vec y)
1095: {
1096: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
1097: PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc*)jac->data;
1098: PetscScalar *yarray;
1099: const PetscScalar *xarray;
1100: KSPConvergedReason reason;
1102: /* place x's and y's local arrays into xsub and ysub */
1103: VecGetArrayRead(x,&xarray);
1104: VecGetArray(y,&yarray);
1105: VecPlaceArray(mpjac->xsub,xarray);
1106: VecPlaceArray(mpjac->ysub,yarray);
1108: /* apply preconditioner on each matrix block */
1109: PetscLogEventBegin(PC_ApplyOnBlocks,jac->ksp[0],mpjac->xsub,mpjac->ysub,0);
1110: KSPSolve(jac->ksp[0],mpjac->xsub,mpjac->ysub);
1111: KSPCheckSolve(jac->ksp[0],pc,mpjac->ysub);
1112: PetscLogEventEnd(PC_ApplyOnBlocks,jac->ksp[0],mpjac->xsub,mpjac->ysub,0);
1113: KSPGetConvergedReason(jac->ksp[0],&reason);
1114: if (reason == KSP_DIVERGED_PC_FAILED) {
1115: pc->failedreason = PC_SUBPC_ERROR;
1116: }
1118: VecResetArray(mpjac->xsub);
1119: VecResetArray(mpjac->ysub);
1120: VecRestoreArrayRead(x,&xarray);
1121: VecRestoreArray(y,&yarray);
1122: return 0;
1123: }
1125: static PetscErrorCode PCMatApply_BJacobi_Multiproc(PC pc,Mat X,Mat Y)
1126: {
1127: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
1128: KSPConvergedReason reason;
1129: Mat sX,sY;
1130: const PetscScalar *x;
1131: PetscScalar *y;
1132: PetscInt m,N,lda,ldb;
1134: /* apply preconditioner on each matrix block */
1135: MatGetLocalSize(X,&m,NULL);
1136: MatGetSize(X,NULL,&N);
1137: MatDenseGetLDA(X,&lda);
1138: MatDenseGetLDA(Y,&ldb);
1139: MatDenseGetArrayRead(X,&x);
1140: MatDenseGetArrayWrite(Y,&y);
1141: MatCreateDense(PetscObjectComm((PetscObject)jac->ksp[0]),m,PETSC_DECIDE,PETSC_DECIDE,N,(PetscScalar*)x,&sX);
1142: MatCreateDense(PetscObjectComm((PetscObject)jac->ksp[0]),m,PETSC_DECIDE,PETSC_DECIDE,N,y,&sY);
1143: MatDenseSetLDA(sX,lda);
1144: MatDenseSetLDA(sY,ldb);
1145: PetscLogEventBegin(PC_ApplyOnBlocks,jac->ksp[0],X,Y,0);
1146: KSPMatSolve(jac->ksp[0],sX,sY);
1147: KSPCheckSolve(jac->ksp[0],pc,NULL);
1148: PetscLogEventEnd(PC_ApplyOnBlocks,jac->ksp[0],X,Y,0);
1149: MatDestroy(&sY);
1150: MatDestroy(&sX);
1151: MatDenseRestoreArrayWrite(Y,&y);
1152: MatDenseRestoreArrayRead(X,&x);
1153: KSPGetConvergedReason(jac->ksp[0],&reason);
1154: if (reason == KSP_DIVERGED_PC_FAILED) {
1155: pc->failedreason = PC_SUBPC_ERROR;
1156: }
1157: return 0;
1158: }
1160: static PetscErrorCode PCSetUp_BJacobi_Multiproc(PC pc)
1161: {
1162: PC_BJacobi *jac = (PC_BJacobi*)pc->data;
1163: PC_BJacobi_Multiproc *mpjac = (PC_BJacobi_Multiproc*)jac->data;
1164: PetscInt m,n;
1165: MPI_Comm comm,subcomm=0;
1166: const char *prefix;
1167: PetscBool wasSetup = PETSC_TRUE;
1168: VecType vectype;
1170: PetscObjectGetComm((PetscObject)pc,&comm);
1172: jac->n_local = 1; /* currently only a single block is supported for a subcommunicator */
1173: if (!pc->setupcalled) {
1174: wasSetup = PETSC_FALSE;
1175: PetscNewLog(pc,&mpjac);
1176: jac->data = (void*)mpjac;
1178: /* initialize datastructure mpjac */
1179: if (!jac->psubcomm) {
1180: /* Create default contiguous subcommunicatiors if user does not provide them */
1181: PetscSubcommCreate(comm,&jac->psubcomm);
1182: PetscSubcommSetNumber(jac->psubcomm,jac->n);
1183: PetscSubcommSetType(jac->psubcomm,PETSC_SUBCOMM_CONTIGUOUS);
1184: PetscLogObjectMemory((PetscObject)pc,sizeof(PetscSubcomm));
1185: }
1186: mpjac->psubcomm = jac->psubcomm;
1187: subcomm = PetscSubcommChild(mpjac->psubcomm);
1189: /* Get matrix blocks of pmat */
1190: MatGetMultiProcBlock(pc->pmat,subcomm,MAT_INITIAL_MATRIX,&mpjac->submats);
1192: /* create a new PC that processors in each subcomm have copy of */
1193: PetscMalloc1(1,&jac->ksp);
1194: KSPCreate(subcomm,&jac->ksp[0]);
1195: KSPSetErrorIfNotConverged(jac->ksp[0],pc->erroriffailure);
1196: PetscObjectIncrementTabLevel((PetscObject)jac->ksp[0],(PetscObject)pc,1);
1197: PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->ksp[0]);
1198: KSPSetOperators(jac->ksp[0],mpjac->submats,mpjac->submats);
1199: KSPGetPC(jac->ksp[0],&mpjac->pc);
1201: PCGetOptionsPrefix(pc,&prefix);
1202: KSPSetOptionsPrefix(jac->ksp[0],prefix);
1203: KSPAppendOptionsPrefix(jac->ksp[0],"sub_");
1204: KSPGetOptionsPrefix(jac->ksp[0],&prefix);
1205: MatSetOptionsPrefix(mpjac->submats,prefix);
1207: /* create dummy vectors xsub and ysub */
1208: MatGetLocalSize(mpjac->submats,&m,&n);
1209: VecCreateMPIWithArray(subcomm,1,n,PETSC_DECIDE,NULL,&mpjac->xsub);
1210: VecCreateMPIWithArray(subcomm,1,m,PETSC_DECIDE,NULL,&mpjac->ysub);
1211: MatGetVecType(mpjac->submats,&vectype);
1212: VecSetType(mpjac->xsub,vectype);
1213: VecSetType(mpjac->ysub,vectype);
1214: PetscLogObjectParent((PetscObject)pc,(PetscObject)mpjac->xsub);
1215: PetscLogObjectParent((PetscObject)pc,(PetscObject)mpjac->ysub);
1217: pc->ops->setuponblocks = PCSetUpOnBlocks_BJacobi_Multiproc;
1218: pc->ops->reset = PCReset_BJacobi_Multiproc;
1219: pc->ops->destroy = PCDestroy_BJacobi_Multiproc;
1220: pc->ops->apply = PCApply_BJacobi_Multiproc;
1221: pc->ops->matapply = PCMatApply_BJacobi_Multiproc;
1222: } else { /* pc->setupcalled */
1223: subcomm = PetscSubcommChild(mpjac->psubcomm);
1224: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
1225: /* destroy old matrix blocks, then get new matrix blocks */
1226: if (mpjac->submats) MatDestroy(&mpjac->submats);
1227: MatGetMultiProcBlock(pc->pmat,subcomm,MAT_INITIAL_MATRIX,&mpjac->submats);
1228: } else {
1229: MatGetMultiProcBlock(pc->pmat,subcomm,MAT_REUSE_MATRIX,&mpjac->submats);
1230: }
1231: KSPSetOperators(jac->ksp[0],mpjac->submats,mpjac->submats);
1232: }
1234: if (!wasSetup && pc->setfromoptionscalled) {
1235: KSPSetFromOptions(jac->ksp[0]);
1236: }
1237: return 0;
1238: }