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

Re: [petsc-dev] PCASM uses DMCreateDecomposition which splits field, NOT expected behavior
by Dmitry Karpeev 26 May '12

26 May '12

On Sat, May 26, 2012 at 9:07 PM, Jed Brown <jedbrown(a)mcs.anl.gov> wrote: > On Sat, May 26, 2012 at 8:04 PM, Dmitry Karpeev <karpeev(a)mcs.anl.gov>wrote: > >> >> >> On Saturday, May 26, 2012, Jed Brown <jedbrown(a)mcs.anl.gov> wrote: >> > Dmitry, I have problems with this commit. >> > http://petsc.cs.iit.edu/petsc/petsc-dev/rev/440ec9350d51 >> > >> > If I try running -pc_type asm with a vector problem defined on a DM, >> PCASM calls DMCreateDecomposition which splits fields on each process. This >> is of very limited value since usually every field is coupled to every >> other, so after increasing the overlap by one, you end up solving the >> entire subdomain twice. However, at this point, you are not using overlap >> so instead, we get subdomain and field blocked Jacobi. (Equivalent to >> BJacobi with additive fieldsplit inside.) This is certainly unexpected >> behavior. >> > It also crashes with Seq(S)BAIJ because MatGetSubMatrices_SeqBAIJ can >> only extract when the block size matches. >> > And it does the wrong thing when used with MPI(S)BAIJ because >> MatGetSubMatrices_MPIBAIJ calls ISCompressIndicesGeneral which silently >> "fills in" the blocks, leading to a Mat and Vec size mismatch. Hong, how >> should we enforce that the original index set actually addresses blocks? >> Why is it not okay to sort? (Don't we rely on those indices being sorted >> somewhere else anyway?) >> >> The problem is that this API mixes field splits and (spatial) domain >> decompositions. >> Initially I thought the two kinds could be served up via the same API >> because they use the same syntax, but the semantics are different. In >> particular, domain decompositions don't change the block size, while field >> splits typically do. Then it's usually downwards -- when splitting up a >> block of fields pointwise, but sometimes upwards -- in a primal-dual split >> of a KKT system the primal subsystem might have a higher natural block size >> than the full system. >> >> I think the solution is to split up the API: DMCreateFieldSplit() and >> DMCreateDomainDecomposition(). I'll look into it later tonight, when I get >> to my computer. > > > If the arguments are definitely going to be the same, there could be an > enum to designate which form is being requested. > I think that might make it more cumbersome for impls to provide only DD or only FS decompositions. > > > What are we going to do now? Leave the code there deactivated like it is > now for the release and discuss later? > I want to fix and reactivate the code for the release so I can use it from libMesh. I'll push roughly this fix by tomorrow: DMCreateDecomposition(DM)()--> DMCreateFieldSplit(DM)()/DMCreateDomainDecomposition(DM)() used, respectively, by PCFIELDSPLIT/PC(G)ASM. That should address the problems you were encountering -- unintended forwarding of field splits to PCASM. > > >> >> >> > Our test coverage must be pretty sucky if we went two months with ASM >> this broken for systems. >> > >> > As a lower priority, -pc_asm_decomposition should provide a list of >> valid decompositions rather than be a magic string. Also, this form cannot >> be called outside of a PetscOptionsBegin/End pair, otherwise it doesn't >> print -help output. It should be moved to PCSetUp_ASM. >> > ierr = PetscOptionsString("-pc_asm_decomposition", "Name of the DM >> defining the composition", "PCSetDM", ddm_name, ddm_name,1024,&flg); >> CHKERRQ(ierr); >> >> The problem here is that the PC object doesn't really know what >> decompositions are available -- the DM does. We could have >> DMGetSubdomainDecompositionNames() advertise what'd available. Another >> problem is that what decompositions make sense is DM-impl-dependent. It is >> more stark for fieldsplits: I have a fork of libMesh that can generate >> decompositions based on the system's variables, and essentially any >> groupings of variable names are valid splits ( e.g., v1,v2;v3,v4. I think >> it makes sense to maintain that flexibility). > > > Cool, what about deeper nestings like [v1,[[v2,v3],v4]]. > Nested splits are to be handled recursively by the individual splits' DMs. And each split should also be domain-decomposable, so that geometrically-defined subdomains (e.g., Exodus II blocks) can be used with PC(G)ASM. I am testing that functionality in libMesh/Moose right now. > > >> Because of that I think the best help line would only exhibit the syntax, >> rather than enumerate all possible splits. > > > How about having a DM-specific function for that designation? > DMXXXGetFieldSplits()/DMXXXGetDomainDecompositions(), for example? I think that would work, if the implementation defined only a small number of decompositions. If instead there is a short grammar defining well-formed strings naming such decompositions, I'm not sure how to communicate that to the user, other than in documentation. I have DMLibMeshGetVars() to return all vars. FieldSplits can be defined as arbitrary groupings of those. Each split will then advertise only the vars in the groupings that defines the split. Dmitry. > > >> Then DMSetUp_Xxx should handle that. >> >> Dmitry. >> >> > >> > >> > I have disabled that code for now, but we have to decide whether it's >> conceptually right (considering the effect on overlap) and if so, how to >> fix it. >> > http://petsc.cs.iit.edu/petsc/petsc-dev/rev/d9a0e3f9dacd >> > >

1 0

Re: [petsc-dev] PCASM uses DMCreateDecomposition which splits field, NOT expected behavior
by Jed Brown 26 May '12

26 May '12

On Sat, May 26, 2012 at 8:04 PM, Dmitry Karpeev <karpeev(a)mcs.anl.gov> wrote: > > > On Saturday, May 26, 2012, Jed Brown <jedbrown(a)mcs.anl.gov> wrote: > > Dmitry, I have problems with this commit. > > http://petsc.cs.iit.edu/petsc/petsc-dev/rev/440ec9350d51 > > > > If I try running -pc_type asm with a vector problem defined on a DM, > PCASM calls DMCreateDecomposition which splits fields on each process. This > is of very limited value since usually every field is coupled to every > other, so after increasing the overlap by one, you end up solving the > entire subdomain twice. However, at this point, you are not using overlap > so instead, we get subdomain and field blocked Jacobi. (Equivalent to > BJacobi with additive fieldsplit inside.) This is certainly unexpected > behavior. > > It also crashes with Seq(S)BAIJ because MatGetSubMatrices_SeqBAIJ can > only extract when the block size matches. > > And it does the wrong thing when used with MPI(S)BAIJ because > MatGetSubMatrices_MPIBAIJ calls ISCompressIndicesGeneral which silently > "fills in" the blocks, leading to a Mat and Vec size mismatch. Hong, how > should we enforce that the original index set actually addresses blocks? > Why is it not okay to sort? (Don't we rely on those indices being sorted > somewhere else anyway?) > > The problem is that this API mixes field splits and (spatial) domain > decompositions. > Initially I thought the two kinds could be served up via the same API > because they use the same syntax, but the semantics are different. In > particular, domain decompositions don't change the block size, while field > splits typically do. Then it's usually downwards -- when splitting up a > block of fields pointwise, but sometimes upwards -- in a primal-dual split > of a KKT system the primal subsystem might have a higher natural block size > than the full system. > > I think the solution is to split up the API: DMCreateFieldSplit() and > DMCreateDomainDecomposition(). I'll look into it later tonight, when I get > to my computer. If the arguments are definitely going to be the same, there could be an enum to designate which form is being requested. What are we going to do now? Leave the code there deactivated like it is now for the release and discuss later? > > > > Our test coverage must be pretty sucky if we went two months with ASM > this broken for systems. > > > > As a lower priority, -pc_asm_decomposition should provide a list of > valid decompositions rather than be a magic string. Also, this form cannot > be called outside of a PetscOptionsBegin/End pair, otherwise it doesn't > print -help output. It should be moved to PCSetUp_ASM. > > ierr = PetscOptionsString("-pc_asm_decomposition", "Name of the DM > defining the composition", "PCSetDM", ddm_name, ddm_name,1024,&flg); > CHKERRQ(ierr); > > The problem here is that the PC object doesn't really know what > decompositions are available -- the DM does. We could have > DMGetSubdomainDecompositionNames() advertise what'd available. Another > problem is that what decompositions make sense is DM-impl-dependent. It is > more stark for fieldsplits: I have a fork of libMesh that can generate > decompositions based on the system's variables, and essentially any > groupings of variable names are valid splits ( e.g., v1,v2;v3,v4. I think > it makes sense to maintain that flexibility). Cool, what about deeper nestings like [v1,[[v2,v3],v4]]. > Because of that I think the best help line would only exhibit the syntax, > rather than enumerate all possible splits. How about having a DM-specific function for that designation? > Then DMSetUp_Xxx should handle that. > > Dmitry. > > > > > > > I have disabled that code for now, but we have to decide whether it's > conceptually right (considering the effect on overlap) and if so, how to > fix it. > > http://petsc.cs.iit.edu/petsc/petsc-dev/rev/d9a0e3f9dacd >

1 0

PCASM uses DMCreateDecomposition which splits field, NOT expected behavior
by Dmitry Karpeev 26 May '12

26 May '12

On Saturday, May 26, 2012, Jed Brown <jedbrown(a)mcs.anl.gov> wrote: > Dmitry, I have problems with this commit. > http://petsc.cs.iit.edu/petsc/petsc-dev/rev/440ec9350d51 > > If I try running -pc_type asm with a vector problem defined on a DM, PCASM calls DMCreateDecomposition which splits fields on each process. This is of very limited value since usually every field is coupled to every other, so after increasing the overlap by one, you end up solving the entire subdomain twice. However, at this point, you are not using overlap so instead, we get subdomain and field blocked Jacobi. (Equivalent to BJacobi with additive fieldsplit inside.) This is certainly unexpected behavior. > It also crashes with Seq(S)BAIJ because MatGetSubMatrices_SeqBAIJ can only extract when the block size matches. > And it does the wrong thing when used with MPI(S)BAIJ because MatGetSubMatrices_MPIBAIJ calls ISCompressIndicesGeneral which silently "fills in" the blocks, leading to a Mat and Vec size mismatch. Hong, how should we enforce that the original index set actually addresses blocks? Why is it not okay to sort? (Don't we rely on those indices being sorted somewhere else anyway?) The problem is that this API mixes field splits and (spatial) domain decompositions. Initially I thought the two kinds could be served up via the same API because they use the same syntax, but the semantics are different. In particular, domain decompositions don't change the block size, while field splits typically do. Then it's usually downwards -- when splitting up a block of fields pointwise, but sometimes upwards -- in a primal-dual split of a KKT system the primal subsystem might have a higher natural block size than the full system. I think the solution is to split up the API: DMCreateFieldSplit() and DMCreateDomainDecomposition(). I'll look into it later tonight, when I get to my computer. > Our test coverage must be pretty sucky if we went two months with ASM this broken for systems. > > As a lower priority, -pc_asm_decomposition should provide a list of valid decompositions rather than be a magic string. Also, this form cannot be called outside of a PetscOptionsBegin/End pair, otherwise it doesn't print -help output. It should be moved to PCSetUp_ASM. > ierr = PetscOptionsString("-pc_asm_decomposition", "Name of the DM defining the composition", "PCSetDM", ddm_name, ddm_name,1024,&flg); CHKERRQ(ierr); The problem here is that the PC object doesn't really know what decompositions are available -- the DM does. We could have DMGetSubdomainDecompositionNames() advertise what'd available. Another problem is that what decompositions make sense is DM-impl-dependent. It is more stark for fieldsplits: I have a fork of libMesh that can generate decompositions based on the system's variables, and essentially any groupings of variable names are valid splits ( e.g., v1,v2;v3,v4. I think it makes sense to maintain that flexibility). Because of that I think the best help line would only exhibit the syntax, rather than enumerate all possible splits. Then DMSetUp_Xxx should handle that. Dmitry. > > > I have disabled that code for now, but we have to decide whether it's conceptually right (considering the effect on overlap) and if so, how to fix it. > http://petsc.cs.iit.edu/petsc/petsc-dev/rev/d9a0e3f9dacd

1 0

PCASM uses DMCreateDecomposition which splits field, NOT expected behavior
by Jed Brown 26 May '12

26 May '12

Dmitry, I have problems with this commit. http://petsc.cs.iit.edu/petsc/petsc-dev/rev/440ec9350d51 If I try running -pc_type asm with a vector problem defined on a DM, PCASM calls DMCreateDecomposition which splits fields on each process. This is of very limited value since usually every field is coupled to every other, so after increasing the overlap by one, you end up solving the entire subdomain twice. However, at this point, you are not using overlap so instead, we get subdomain and field blocked Jacobi. (Equivalent to BJacobi with additive fieldsplit inside.) This is certainly unexpected behavior. It also crashes with Seq(S)BAIJ because MatGetSubMatrices_SeqBAIJ can only extract when the block size matches. And it does the wrong thing when used with MPI(S)BAIJ because MatGetSubMatrices_MPIBAIJ calls ISCompressIndicesGeneral which silently "fills in" the blocks, leading to a Mat and Vec size mismatch. Hong, how should we enforce that the original index set actually addresses blocks? Why is it not okay to sort? (Don't we rely on those indices being sorted somewhere else anyway?) Our test coverage must be pretty sucky if we went two months with ASM this broken for systems. As a lower priority, -pc_asm_decomposition should provide a list of valid decompositions rather than be a magic string. Also, this form cannot be called outside of a PetscOptionsBegin/End pair, otherwise it doesn't print -help output. It should be moved to PCSetUp_ASM. ierr = PetscOptionsString("-pc_asm_decomposition", "Name of the DM defining the composition", "PCSetDM", ddm_name, ddm_name,1024,&flg); CHKERRQ(ierr); I have disabled that code for now, but we have to decide whether it's conceptually right (considering the effect on overlap) and if so, how to fix it. http://petsc.cs.iit.edu/petsc/petsc-dev/rev/d9a0e3f9dacd

1 0

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

26 May '12

On Sat, May 26, 2012 at 9:39 AM, Stefano Zampini <stefano.zampini(a)gmail.com>wrote: > One we have the constraint matrix, we can easily obtain the change of > basis matrix T (as in Klawonn-Widlund papers). > Note that the change of basis approach will be very effective for exact > applications with reduced iterations. I think we should include in the new > matrix class the possibility of doing iterations on the reduced space > instead of the whole space of dofs. > Indeed. You might be aware of the "second approach" described here (CPAM 2006) https://ftp.cs.nyu.edu/web/Research/TechReports/TR2004-855/TR2004-855.pdf I think it's worth testing, but since the number of integral Lagrange multipliers is reasonably small, I would not a priori be concerned by solving the saddle point problem in which they are ordered last. If there are enough vertex constraints to keep the subdomains nonsingular, you can easily handle the problem manually (as Dohrmann does), but I wouldn't expect modern LDL^T packages to be slowed greatly by a handful of Lagrange multipliers and you want it to reorder anyway to reduce fill during the factorization. > >> >>> 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 using the equivalence. > > You said you would have the new matrix class to support either more > subdomains per core, or more cores per subdomain. In the latter case, > threaded or mpi matrices (on subcomms)? > I would plan to make it support any combination. Tell me if you think it's a bad idea, but my thought was to pre-select a few coarse vertices (ideally enough to prevent floating subdomains, but the method should be able to tolerate if not). Then the form that gets assembled due to normal MatSetValuesLocal() would have nonsingular subdomains. We can factor those subdomains and then adaptively select additional vertex or integral constraints to enrich the coarse space. As soon as we have enough constraints to make the subdomain nonsingular (i.e. we would usually start nonsingular), we can enforce all additional constraints with Lagrange multipliers so that the factorization does not need to be repeated. Once we have enriched the coarse space, we have a choice of whether to apply a change of basis and re-factor or to reuse the current factorization and keep enforcing the enrichment by Lagrange multipliers.

1 0

Re: [petsc-dev] API changes in MatIS
by Stefano Zampini 26 May '12

26 May '12

> Okay. Presumably the change of basis can be done with a MatPtAP if > desirable. > One we have the constraint matrix, we can easily obtain the change of basis matrix T (as in Klawonn-Widlund papers). Note that the change of basis approach will be very effective for exact applications with reduced iterations. I think we should include in the new matrix class the possibility of doing iterations on the reduced space instead of the whole space of dofs. > > >> 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 using the equivalence. You said you would have the new matrix class to support either more subdomains per core, or more cores per subdomain. In the latter case, threaded or mpi matrices (on subcomms)? -- Stefano

1 0

Re: [petsc-dev] What needs to be done before a release?
by Barry Smith 25 May '12

25 May '12

Tim, Thanks. But what I want is all bundled up a routine that takes in the Fortran array of character strings and allocates and copies over them to the C ptr business. Likely I could figure it out myself with but Blaise is such a Fortran wizard he can do it optimally in his sleep :-) Barry On May 25, 2012, at 8:14 PM, Tim Gallagher wrote: > For what it's worth, I answered a question about how to pass string arrays awhile ago on StackOverflow: > > http://stackoverflow.com/questions/9686532/arrays-of-strings-in-fortran-c-b… > > Tim > > ----- Original Message ----- > From: "Barry Smith" <bsmith(a)mcs.anl.gov> > To: "For users of the development version of PETSc" <petsc-dev(a)mcs.anl.gov> > Sent: Friday, May 25, 2012 8:57:52 PM > Subject: Re: [petsc-dev] What needs to be done before a release? > > Blaise > > If you provide a fortran function to convert an array of Fortran strings to an array of C strings and tell us how to delete the result then we'll provide a PetscOptionsGetEnum() for Fortran. > > Barry > > On May 2, 2012, at 12:02 PM, Blaise Bourdin wrote: > >> Barry, >> >>>> Also, fortran wrappers over PetscOptionsEnum and PetscBagRegisterEnum would be nice. I understand that it >>>> is a bit tricky since it involves passing fortran arrays of strings. >>> >>> Yup. We'd be willing to support these from FORTRAN but someone else has to do the dirty work of figuring out how to handle those array of strings (maybe F2003?). >> >> I am attaching a simple proof of concept example that demonstrates how to write a fortran interface and call petsc C functions from fortran using F2003 C-interoperability instead of going through the fortran bindings. >> >> Don't be put off by all the warnings if you use gfortran. They would be easy to avoid. Basically, right now, you detect the fortran integer and real kinds that match that of the C ones. It would be quite simple to get the fortran C-interoperable integer and real types associated with the C type of PetscInt, PetscReal etc. >> >> Basically, we (I?) could write a fortran interface and a fortran function for the tricky functions involving arrays of string, instead of creating a custom binding in C. Are you willing to let some fortran pollution into the petsc code? Where would these function have to be? >> >> Blaise >> >> PS: you need FCFLAGS='-ffixed-line-length-none -ffree-line-length-none -ffree-form' or some subset of these to compile with gfortran. I'm not sure why these are not the default with a file suffix .F90 >> >> <PetscF2003Binding.F90> >> -- >> Department of Mathematics and Center for Computation & Technology >> Louisiana State University, Baton Rouge, LA 70803, USA >> Tel. +1 (225) 578 1612, Fax +1 (225) 578 4276 http://www.math.lsu.edu/~bourdin >> >> >> >> >> >> >> >

1 0

Re: [petsc-dev] What needs to be done before a release?
by Tim Gallagher 25 May '12

25 May '12

For what it's worth, I answered a question about how to pass string arrays awhile ago on StackOverflow: http://stackoverflow.com/questions/9686532/arrays-of-strings-in-fortran-c-b… Tim ----- Original Message ----- From: "Barry Smith" <bsmith(a)mcs.anl.gov> To: "For users of the development version of PETSc" <petsc-dev(a)mcs.anl.gov> Sent: Friday, May 25, 2012 8:57:52 PM Subject: Re: [petsc-dev] What needs to be done before a release? Blaise If you provide a fortran function to convert an array of Fortran strings to an array of C strings and tell us how to delete the result then we'll provide a PetscOptionsGetEnum() for Fortran. Barry On May 2, 2012, at 12:02 PM, Blaise Bourdin wrote: > Barry, > >>> Also, fortran wrappers over PetscOptionsEnum and PetscBagRegisterEnum would be nice. I understand that it >>> is a bit tricky since it involves passing fortran arrays of strings. >> >> Yup. We'd be willing to support these from FORTRAN but someone else has to do the dirty work of figuring out how to handle those array of strings (maybe F2003?). > > I am attaching a simple proof of concept example that demonstrates how to write a fortran interface and call petsc C functions from fortran using F2003 C-interoperability instead of going through the fortran bindings. > > Don't be put off by all the warnings if you use gfortran. They would be easy to avoid. Basically, right now, you detect the fortran integer and real kinds that match that of the C ones. It would be quite simple to get the fortran C-interoperable integer and real types associated with the C type of PetscInt, PetscReal etc. > > Basically, we (I?) could write a fortran interface and a fortran function for the tricky functions involving arrays of string, instead of creating a custom binding in C. Are you willing to let some fortran pollution into the petsc code? Where would these function have to be? > > Blaise > > PS: you need FCFLAGS='-ffixed-line-length-none -ffree-line-length-none -ffree-form' or some subset of these to compile with gfortran. I'm not sure why these are not the default with a file suffix .F90 > > <PetscF2003Binding.F90> > -- > Department of Mathematics and Center for Computation & Technology > Louisiana State University, Baton Rouge, LA 70803, USA > Tel. +1 (225) 578 1612, Fax +1 (225) 578 4276 http://www.math.lsu.edu/~bourdin > > > > > > >

1 0

Re: [petsc-dev] What needs to be done before a release?
by Barry Smith 25 May '12

25 May '12

Blaise If you provide a fortran function to convert an array of Fortran strings to an array of C strings and tell us how to delete the result then we'll provide a PetscOptionsGetEnum() for Fortran. Barry On May 2, 2012, at 12:02 PM, Blaise Bourdin wrote: > Barry, > >>> Also, fortran wrappers over PetscOptionsEnum and PetscBagRegisterEnum would be nice. I understand that it >>> is a bit tricky since it involves passing fortran arrays of strings. >> >> Yup. We'd be willing to support these from FORTRAN but someone else has to do the dirty work of figuring out how to handle those array of strings (maybe F2003?). > > I am attaching a simple proof of concept example that demonstrates how to write a fortran interface and call petsc C functions from fortran using F2003 C-interoperability instead of going through the fortran bindings. > > Don't be put off by all the warnings if you use gfortran. They would be easy to avoid. Basically, right now, you detect the fortran integer and real kinds that match that of the C ones. It would be quite simple to get the fortran C-interoperable integer and real types associated with the C type of PetscInt, PetscReal etc. > > Basically, we (I?) could write a fortran interface and a fortran function for the tricky functions involving arrays of string, instead of creating a custom binding in C. Are you willing to let some fortran pollution into the petsc code? Where would these function have to be? > > Blaise > > PS: you need FCFLAGS='-ffixed-line-length-none -ffree-line-length-none -ffree-form' or some subset of these to compile with gfortran. I'm not sure why these are not the default with a file suffix .F90 > > <PetscF2003Binding.F90> > -- > Department of Mathematics and Center for Computation & Technology > Louisiana State University, Baton Rouge, LA 70803, USA > Tel. +1 (225) 578 1612, Fax +1 (225) 578 4276 http://www.math.lsu.edu/~bourdin > > > > > > >

1 0

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

25 May '12

On Fri, May 25, 2012 at 6:57 PM, Mark F. Adams <mark.adams(a)columbia.edu>wrote: > Yes, I forgot what we were arguing about there for a minute :) > > My thinking was that Cheby _could_ degrade performance over Richardson > (for G-S in particular). That is that I could construct a > smoother/operator for which Cheby would not pick the correct damping (so if > the PC/operator was damped perfectly Richardson could win). > > I'm not sure if that is true now. > So you can of course make it worse for some modes, but only by a modest factor for real eigenvalues, and Cheby generally hits a broader part of the spectrum. If the preconditioned operator has eigenvalues with large imaginary components, Cheby can become unstable. For that case (and similar in the nonlinear context), I would like to be able to estimate a few eigenvalues and optimize a polynomial smoother that damps those regions. For example, attached is an unoptimized example of a fourth order polynomial. octave> xx = linspace(0,3,30); yy = linspace(-3,3,60); [x,y] = meshgrid(xx,yy); z = x+1i*y; contourf(x,y,1-z/(1+2i)).*(1-z/(1-2i)).*(1-z/(1.5-.9i)).*(1-z/(1.5+.9i))), [.2 .4 .6 .8 1]) Everything in the dark blue region has a damping factor smaller than 0.4, everything in the light blue region is damped by 0.6. If the extreme eigenvalues get closer to the real axis, we can find polynomials that damp better (or reduce the degree of the polynomial). There is plenty of related work in the ODE literature, but they also care about accuracy, which doesn't matter at all to us.

1 0