May 2012 - petsc-dev - mailman3.cels.anl.gov

A PETSc question on scicomp I can't completely answer (behavior of PETScSetxxx)
by Aron Ahmadia 23 May '12

23 May '12

http://scicomp.stackexchange.com/questions/2303/petscs-xxxsetxxx-methods-ow… My understanding is that Set for the most part acts like a pass-by-reference, and my answer reflects that ( http://scicomp.stackexchange.com/a/2311/9) If there's a better answer than that, Jed or Matt should post it or fix mine :) A

1 0

Re: [petsc-dev] Multigrid is confusing
by Barry Smith 22 May '12

22 May '12

On May 22, 2012, at 9:58 PM, Jed Brown wrote: > On Tue, May 22, 2012 at 9:52 PM, Barry Smith <bsmith(a)mcs.anl.gov> wrote: > How about a default nonlinear Krylov accelerator? fgrmes is not much more expensive than gmres under causal circumstances. > > I could handle that. The extra storage concerns me more than time. > > I'd eventually like to make PCMG linear by default, especially considering that the convergence with GMRES(1) smoothers tends to be worse in my experience (except when coarse level operators are poorly scaled/unstable) and adds a ton of synchronization points. > > Oh, FGMRES by default would necessarily switch from left to right preconditioning. I'm not sure that's a good thing since it requires much tighter tolerances for ill-conditioned systems to reach a given error (especially penalty boundaries). There are advantages both ways and I keep waffling. I expect that changing the default will cause some surprises. So does it make sense to use ILU(0) smoothing with Chebyshev(3) acceleration (with eigenvalues precomputed by GMRES?). I'm fine with switching now, I just want something reasonably robust so for example ex19 converges well with the defaults. Barry

1 0

Re: [petsc-dev] Multigrid is confusing
by Jed Brown 22 May '12

22 May '12

On Tue, May 22, 2012 at 9:52 PM, Barry Smith <bsmith(a)mcs.anl.gov> wrote: > How about a default nonlinear Krylov accelerator? fgrmes is not much more > expensive than gmres under causal circumstances. I could handle that. The extra storage concerns me more than time. I'd eventually like to make PCMG linear by default, especially considering that the convergence with GMRES(1) smoothers tends to be worse in my experience (except when coarse level operators are poorly scaled/unstable) and adds a ton of synchronization points. Oh, FGMRES by default would necessarily switch from left to right preconditioning. I'm not sure that's a good thing since it requires much tighter tolerances for ill-conditioned systems to reach a given error (especially penalty boundaries). There are advantages both ways and I keep waffling. I expect that changing the default will cause some surprises.

1 0

Re: [petsc-dev] Multigrid is confusing
by Barry Smith 22 May '12

22 May '12

How about a default nonlinear Krylov accelerator? fgrmes is not much more expensive than gmres under causal circumstances. Barry On May 22, 2012, at 7:54 PM, Jed Brown wrote: > ~/petsc/src/snes/examples/tutorials$ ./ex5 -da_refine 4 -ksp_monitor_true_residual -pc_type mg -snes_max_it 1 -ksp_rtol 1e-8 > 0 KSP preconditioned resid norm 5.024049081009e+00 true resid norm 1.141253746897e+00 ||r(i)||/||b|| 1.000000000000e+00 > 1 KSP preconditioned resid norm 2.430101772479e-01 true resid norm 3.124735997267e-02 ||r(i)||/||b|| 2.737985312875e-02 > 2 KSP preconditioned resid norm 1.925127800529e-03 true resid norm 1.276301220139e-03 ||r(i)||/||b|| 1.118332556287e-03 > 3 KSP preconditioned resid norm 3.150352958424e-05 true resid norm 2.212937904665e-03 ||r(i)||/||b|| 1.939041085896e-03 > 4 KSP preconditioned resid norm 5.059562742413e-07 true resid norm 2.180293054518e-03 ||r(i)||/||b|| 1.910436710895e-03 > 5 KSP preconditioned resid norm 1.095269606918e-08 true resid norm 2.180508594793e-03 ||r(i)||/||b|| 1.910625573605e-03 > > I think we need to make default -pc_type mg linear. We talked about this before, but now we are getting ready for the release and that's bound to cause some turmoil. > > Thoughts?

1 0

Re: [petsc-dev] anyone feel like fixing a bug?
by Jed Brown 22 May '12

22 May '12

Well, here's the problem: MatSetBlockSize() used to call mat->ops->setblocksize() which the MPI matrices forwarded down to their seq submatrices. But then this patch removed that behavior so MatSetBlockSize() is processed immediately and not forwarded down. http://petsc.cs.iit.edu/petsc/petsc-dev/rev/80ec81eb9db9 Since the Seq inner matrices have different PetscLayouts, they don't inherit the block sizes set on their parent. Barry, why did you remove this MatOp? How do you want the MPI matrices to forward block sizes now? On Sat, May 19, 2012 at 10:57 AM, Mark F. Adams <mark.adams(a)columbia.edu>wrote: > I'm going to be tied up for a few days, should be able to get to this > Wed., if not sooner, but if anyone wants to take a crack at it, here it is. > This is run with 8 procs. > > Mark > > ==51037== Command: ./ex56 -ne 19 -alpha 1.e-3 -ksp_type cg -ksp_monitor > -ksp_max_it 15 -pc_type gamg -pc_gamg_type agg -pc_gamg_agg_nsmooths 1 > -pc_gamg_verbose 2 -ksp_converged_reason -options_left -theta -30.0 > -epsilon .01 -blob_center 0.,0. -pc_gamg_threshold .0 > -pc_gamg_coarse_eq_limit 100 -use_coordinates -out_verbose 2 -aggmg_smooths > 1 -mg_levels_ksp_type chebyshev > -mg_levels_ksp_chebyshev_estimate_eigenvalues 0,0.2,0,1.1 > -mg_levels_pc_type pbjacobi -mg_levels_ksp_max_it 1 > > ==51039== Invalid read of size 8 > ==51039== at 0x91E38B: PCApply_PBJacobi_3 (pbjacobi.c:81) > ==51039== by 0xA0FE23: PCApply (precon.c:384) > ==51039== by 0xADA9B2: KSPInitialResidual (itres.c:64) > ==51039== by 0xA7EA39: KSPSolve_GMRES (gmres.c:230) > ==51039== by 0xAB3CA2: KSPSolve (itfunc.c:446) > ==51039== by 0xAA6FF8: KSPSolve_Chebyshev (cheby.c:269) > ==51039== by 0xAB3CA2: KSPSolve (itfunc.c:446) > ==51039== by 0x9ECD2F: PCMGMCycle_Private (mg.c:20) > ==51039== by 0x9F0BEF: PCApply_MG (mg.c:321) > ==51039== by 0xA0FE23: PCApply (precon.c:384) > ==51039== by 0xA43EBF: KSPSolve_CG (cg.c:139) > ==51039== by 0xAB3CA2: KSPSolve (itfunc.c:446) > ==51039== by 0x4805: main (in ./ex56) > ==51039== Address 0x27f4c90 is 0 bytes after a block of size 24,032 > alloc'd > ==51039== at 0x1B2AF16: malloc (vg_replace_malloc.c:236) > ==51039== by 0x8A8C7: PetscMallocAlign (mal.c:37) > ==51039== by 0x61E602: MatInvertBlockDiagonal_SeqAIJ (aij.c:2816) > ==51039== by 0x2F0AC3: MatInvertBlockDiagonal (matrix.c:9280) > ==51039== by 0x6D1494: MatInvertBlockDiagonal_MPIAIJ (mpiaij.c:2981) > ==51039== by 0x2F0AC3: MatInvertBlockDiagonal (matrix.c:9280) > ==51039== by 0x920A85: PCSetUp_PBJacobi (pbjacobi.c:213) > ==51039== by 0xA1783B: PCSetUp (precon.c:832) > ==51039== by 0xAB2125: KSPSetUp (itfunc.c:278) > ==51039== by 0x9F536C: PCSetUp_MG (mg.c:678) > ==51039== by 0x98C26C: PCSetUp_GAMG (gamg.c:984) > ==51039== by 0xA1783B: PCSetUp (precon.c:832) > ==51039== by 0xAB2125: KSPSetUp (itfunc.c:278) > ==51039== by 0x46AD: main (in ./ex56) > ==51039== > >

1 0

Multigrid is confusing
by Jed Brown 22 May '12

22 May '12

~/petsc/src/snes/examples/tutorials$ ./ex5 -da_refine 4 -ksp_monitor_true_residual -pc_type mg -snes_max_it 1 -ksp_rtol 1e-8 0 KSP preconditioned resid norm 5.024049081009e+00 true resid norm 1.141253746897e+00 ||r(i)||/||b|| 1.000000000000e+00 1 KSP preconditioned resid norm 2.430101772479e-01 true resid norm 3.124735997267e-02 ||r(i)||/||b|| 2.737985312875e-02 2 KSP preconditioned resid norm 1.925127800529e-03 true resid norm 1.276301220139e-03 ||r(i)||/||b|| 1.118332556287e-03 3 KSP preconditioned resid norm 3.150352958424e-05 true resid norm 2.212937904665e-03 ||r(i)||/||b|| 1.939041085896e-03 4 KSP preconditioned resid norm 5.059562742413e-07 true resid norm 2.180293054518e-03 ||r(i)||/||b|| 1.910436710895e-03 5 KSP preconditioned resid norm 1.095269606918e-08 true resid norm 2.180508594793e-03 ||r(i)||/||b|| 1.910625573605e-03 I think we need to make default -pc_type mg linear. We talked about this before, but now we are getting ready for the release and that's bound to cause some turmoil. Thoughts?

1 0

Re: [petsc-dev] API changes in MatIS
by Jed Brown 22 May '12

22 May '12

On Mon, May 21, 2012 at 3:45 PM, Stefano Zampini <stefano.zampini(a)gmail.com>wrote: > Yes. In the rectangular case you will need both global_to_local and > local_to global scatters. The square case is a special case for which they > are the same (except that they are performed in SCATTER_FORWARD or > SCATTER_REVERSE mode). With rectangular matrices, you can perform MatMult > using the SCATTER_FORWARD mode for the two different scatters; > MatMultTranspose using SCATTER REVERSEs. > Yes, scatter is fine. An ISGlobalToLocalMapping (allowing you to translate global indices into local indices) is not. > > The Dohrmann's approach for solving the local saddle point problems is > strictly related on how the are currently solved in PCBDDC. I'm planning to > add a new feature in PCBDDC which can give the possibility to factor and > solve the local saddle point problems directly (factoring the whole local > saddle point problems). I used the change of basis approach in the past, > but I'm more inclined to adopt the saddle point approach, because you can > (almost) easily switch inexact subdomain solvers. However, I think the > change of basis can also be added to the current method with minor code > changes. > Okay. Presumably the change of basis can be done with a MatPtAP if desirable. > Actually in PCBDDC the local coarse matrix is computed using a theoretical > equivalence of the PtAP operation of coarse basis and the unassembled MATIS > matrix (coarse basis are continuous only at vertex and constraints dofs). > PtAP (where P is dense) is just avoided for its computational costs. > Are you basically just doing a local PtAP or do you use the equivalence K \Psi = C^T \Lambda (notation of Dohrmann's Eq 2) or something else? > I'm not familiar with multigrid, but if you can implement an MG framework > for the application of a nonoverlapping (or strip) preconditioner it would > be great. I can help you with the Neumann-Neumann and BDDC cases. > Cool, we can definitely do this. The software reuse is somewhat delicate when asked to also support iterating in the reduced space, but I don't think that's as big of a deal in 3D (the reduced space isn't massively smaller than the full space unless you have huge subdomains, which cause their own set of problems including limited suitability of direct solvers). > I think it is a problem. All MATIS code is inherenlty written for > uniprocessor subdomains and one subdomain per process. Maybe it would be > better to think at a brand new class of matrices (your PA I think) for > which construct the new class and subclasses of preconditioners. Maybe such > a new class can be implemented with SF instead of VecScatter. > That sounds fine to me, hopefully we can unify later. > PCISCreateRestrictionAndProlongation_BDDC: (in case of exact solvers for >>> the Dirichlet problems) default P will be of size >>> n_coarse_dofs*sum^N_{i=1}pcis->n_B with local matrices of P the actual >>> pcbddc->coarse_phi_B >>> >>> PCISCreateRestrictionAndProlongation_BDDC: (in case of inexact solvers >>> for the Dirichlet problems) default P will be of size >>> n_coarse_dofs*sum^N_{i=1}pcis->n with local matrices of P the actual >>> pcbddc->coarse_phi_B concatenated with pcbddc->coarse_phi_D >>> (pcis->n=pcis->n_B+pcis->n_D) >>> >> >> This sounds okay. The reduced space iteration with exact Dirichlet >> solvers bothers me somewhat. We should be able to implement as running >> PCFieldSplit to restrict the inner iteration to the interface problem, but >> with our current data structures, that may have thrown away the interior >> information that we need. >> > > I'm not getting the point here. Can you clarify? > So when exact subdomain solves are used, people ask for the interior dofs to be eliminated and the Krylov space is constructed in the interface space. This means that each application of the operator involves solving a subdomain Dirichlet problem. This is similar to standard FETI-DP where the iteration takes place in the space of Lagrange multipliers and the operator solves the local pinned Neumann problems and coarse grid correction. Working in the reduced space requires exact solves, but reduces the size of the vectors in the Krylov space (just the length of the vectors, it's irrelevant for convergence rate).

1 0

Re: [petsc-dev] cheby.c
by Hong Zhang 22 May '12

22 May '12

Barry, Do you recall the paper or detailed Chebyshev algorithm that you implemented the Chebyshev in petsc? I added 'hybrid' to it, which is not correct. Jie is helping to fix it. Having an original algorithm would help. On Tue, May 22, 2012 at 1:21 PM, Jie Chen <jiechen(a)mcs.anl.gov> wrote: > >> Then do you know who wrote the code? Hong said it was you. >> > > It was written in the 90s, before version control history. > Sorry, my mistake. I recall that Barry said that he wrote it when Jed was in dia... ;-( Hong > > >> >> Hong is trying to modify the code but there is something wrong. I am >> helping to see what is the problem. The code itself is fine. > > > Hard to speculate without more than "something wrong". > > >> >> Jie >> >> >> >> ----- Original Message ----- >> From: "Jed Brown" <jedbrown(a)mcs.anl.gov> >> To: "Jie Chen" <jiechen(a)mcs.anl.gov>, "For users of the development >> version of PETSc" <petsc-dev(a)mcs.anl.gov> >> Sent: Tuesday, May 22, 2012 1:18:07 PM >> Subject: Re: cheby.c >> >> >> On Tue, May 22, 2012 at 12:37 PM, Jie Chen < jiechen(a)mcs.anl.gov > wrote: >> >> >> Jed, I am reading your implementation of the chebyshev code. The >> algorithm does not seem to be the same as the version I am familiar with. >> Did you implement it according to some paper or book? >> >> >> >> I didn't write it. >> >> >> What do you find strange? >> > >

1 0

Re: [petsc-dev] cheby.c
by Jed Brown 22 May '12

22 May '12

On Tue, May 22, 2012 at 1:21 PM, Jie Chen <jiechen(a)mcs.anl.gov> wrote: > Then do you know who wrote the code? Hong said it was you. > It was written in the 90s, before version control history. > > Hong is trying to modify the code but there is something wrong. I am > helping to see what is the problem. The code itself is fine. > Hard to speculate without more than "something wrong". > > Jie > > > > ----- Original Message ----- > From: "Jed Brown" <jedbrown(a)mcs.anl.gov> > To: "Jie Chen" <jiechen(a)mcs.anl.gov>, "For users of the development > version of PETSc" <petsc-dev(a)mcs.anl.gov> > Sent: Tuesday, May 22, 2012 1:18:07 PM > Subject: Re: cheby.c > > > On Tue, May 22, 2012 at 12:37 PM, Jie Chen < jiechen(a)mcs.anl.gov > wrote: > > > Jed, I am reading your implementation of the chebyshev code. The algorithm > does not seem to be the same as the version I am familiar with. Did you > implement it according to some paper or book? > > > > I didn't write it. > > > What do you find strange? >

1 0

Re: [petsc-dev] cheby.c
by Jie Chen 22 May '12

22 May '12

Then do you know who wrote the code? Hong said it was you. Hong is trying to modify the code but there is something wrong. I am helping to see what is the problem. The code itself is fine. Jie ----- Original Message ----- From: "Jed Brown" <jedbrown(a)mcs.anl.gov> To: "Jie Chen" <jiechen(a)mcs.anl.gov>, "For users of the development version of PETSc" <petsc-dev(a)mcs.anl.gov> Sent: Tuesday, May 22, 2012 1:18:07 PM Subject: Re: cheby.c On Tue, May 22, 2012 at 12:37 PM, Jie Chen < jiechen(a)mcs.anl.gov > wrote: Jed, I am reading your implementation of the chebyshev code. The algorithm does not seem to be the same as the version I am familiar with. Did you implement it according to some paper or book? I didn't write it. What do you find strange?

1 0