On May 25, 2013, at 9:06 PM, Jed Brown <[email protected]> wrote:
> Barry Smith <[email protected]> writes: > >>> Instead, it >>> uses block Jacobi in which all blocks but one are empty. >> >> This seems a rather ad hoc way of handling it. Why not have a >> proper mechanism for telescoping solvers rather than an ad hoc >> block Jacobi with 0 sizes on most blocks? > > How would that "proper mechanism" look? Beats me, I didn't write this solver infrastructure but it must be possible to design a reasonable infrastructure to handle this kind of thing. If it is not possible then MPI sucks. > >>>> And why does it have different solvers on each process? >>> >>> It doesn't, but it calls PCBJacobiGetSubKSP() so as to configure the >>> inner PC to be PCLU, but PCBJacobiGetSubKSP conservatively sets this >>> flag: >>> >>> jac->same_local_solves = PETSC_FALSE; /* Assume that local solves >>> are now different; >>> not necessarily true >>> though! This flag is >>> used only for >>> PCView_BJacobi() */ >> >> Why not just fix this instead. Fix it so it doesn't set the flag in this >> situation. > > Sure, where would we set that flag? Or, how do we compare solvers (on > different subcomms) for "equivalent" configuration? { /* coarse grid */ KSP smoother,*k2; PC subpc,pc2; PetscInt ii,first; Mat Lmat = Aarr[(level=pc_gamg->Nlevels-1)]; lidx = 0; ierr = PCMGGetSmoother(pc, lidx, &smoother);CHKERRQ(ierr); ierr = KSPSetOperators(smoother, Lmat, Lmat, SAME_NONZERO_PATTERN);CHKERRQ(ierr); ierr = KSPSetNormType(smoother, KSP_NORM_NONE);CHKERRQ(ierr); ierr = KSPGetPC(smoother, &subpc);CHKERRQ(ierr); ierr = PCSetType(subpc, PCBJACOBI);CHKERRQ(ierr); ierr = PCSetUp(subpc);CHKERRQ(ierr); ierr = PCBJacobiGetSubKSP(subpc,&ii,&first,&k2);CHKERRQ(ierr); if (ii != 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"ii %D is not one",ii); ierr = KSPGetPC(k2[0],&pc2);CHKERRQ(ierr); ierr = PCSetType(pc2, PCLU);CHKERRQ(ierr); ierr = PCFactorSetShiftType(pc2,MAT_SHIFT_INBLOCKS);CHKERRQ(ierr); ierr = KSPSetTolerances(k2[0],PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,1);CHKERRQ(ierr); } > The code is already completely hardwired to use block Jacobi and LU so for maint I would just hardwire setting the flag here to same . That's right, just #include <../src/ksp/pc/impls/bjacobi/bjacobi.h> and go to town. In the future the code can be refactored to not use an ad hoc bjacobi to manage the telescoping. Barry
