[petsc-users] using preconditioner with SLEPc

Matthew Knepley knepley at gmail.com
Mon Feb 15 07:53:49 CST 2021 

On Mon, Feb 15, 2021 at 7:27 AM Jose E. Roman <jroman at dsic.upv.es> wrote:

> I will think about the viability of adding an interface function to pass
> the preconditioner matrix.
>
> Regarding the question about the B-orthogonality of computed vectors, in
> the symmetric solver the B-orthogonality is enforced during the
> computation, so you have guarantee that the computed vectors satisfy it.
> But if solved as non-symetric, the computed vectors may depart from
> B-orthogonality, unless the tolerance is very small.
>

Yes, the vectors I generate are not B-orthogonal.

Jose, do you think there is a way to reformulate what I am doing to use the
symmetric solver, even if we only have the action of B?

  Thanks,

     Matt


> Jose
>
>
> > El 14 feb 2021, a las 21:41, Barry Smith <bsmith at petsc.dev> escribió:
> >
> >
> >   Florian,
> >
> >    I'm sorry I don't know the answers; I can only speculate. There is a
> STGetShift().
> >
> >    All I was saying is theoretically there could/should be such support
> in SLEPc.
> >
> >   Barry
> >
> >
> >> On Feb 13, 2021, at 6:43 PM, Florian Bruckner <e0425375 at gmail.com>
> wrote:
> >>
> >> Dear Barry,
> >> thank you for your clarification. What I wanted to say is that even if
> I could reset the KSP operators directly I would require to know which
> transformation ST applies in order to provide the preconditioning matrix
> for the correct operator.
> >> The more general solution would be that SLEPc provides the interface to
> pass the preconditioning matrix for A0 and ST applies the same
> transformations as for the operator.
> >>
> >> If you write "SLEPc could provide an interface", do you mean someone
> should implement it, or should it already be possible and I am not using it
> correctly?
> >> I wrote a small standalone example based on ex9.py from slepc4py, where
> i tried to use an operator.
> >>
> >> best wishes
> >> Florian
> >>
> >> On Sat, Feb 13, 2021 at 7:15 PM Barry Smith <bsmith at petsc.dev> wrote:
> >>
> >>
> >>> On Feb 13, 2021, at 2:47 AM, Pierre Jolivet <pierre at joliv.et> wrote:
> >>>
> >>>
> >>>
> >>>> On 13 Feb 2021, at 7:25 AM, Florian Bruckner <e0425375 at gmail.com>
> wrote:
> >>>>
> >>>> Dear Jose, Dear Barry,
> >>>> thanks again for your reply. One final question about the B0
> orthogonality. Do you mean that eigenvectors are not B0 orthogonal, but
> they are i*B0 orthogonal? or is there an issue with Matt's approach?
> >>>> For my problem I can show that eigenvalues fulfill an orthogonality
> relation (phi_i, A0 phi_j ) = omega_i (phi_i, B0 phi_j) = delta_ij. This
> should be independent of the solving method, right?
> >>>>
> >>>> Regarding Barry's advice this is what I first tried:
> >>>> es = SLEPc.EPS().create(comm=fd.COMM_WORLD)
> >>>> st = es.getST()
> >>>> ksp = st.getKSP()
> >>>> ksp.setOperators(self.A0, self.P0)
> >>>>
> >>>> But it seems that the provided P0 is not used. Furthermore the
> interface is maybe a bit confusing if ST performs some transformation. In
> this case P0 needs to approximate A0^{-1}*B0 and not A0, right?
> >>>
> >>> No, you need to approximate (A0-sigma B0)^-1. If you have a null
> shift, which looks like it is the case, you end up with A0^-1.
> >>
> >>   Just trying to provide more clarity with the terms.
> >>
> >> If ST transforms the operator in the KSP to (A0-sigma B0) and you are
> providing the "sparse matrix from which the preconditioner is to be built"
> then you need to provide something that approximates (A0-sigma B0). Since
> the PC will use your matrix to construct a preconditioner that approximates
> the inverse of  (A0-sigma B0), you don't need to directly provide something
> that approximates (A0-sigma B0)^-1
> >>
> >>  Yes, I would think SLEPc could provide an interface where it manages
> "the matrix from which to construct the preconditioner" and transforms that
> matrix just like the true matrix. To do it by hand you simply need to know
> what A0 and B0 are and which sigma ST has selected and then you can
> construct your modA0 - sigma modB0 and pass it to the KSP. Where modA0 and
> modB0 are your "sparser approximations".
> >>
> >>   Barry
> >>
> >>
> >>>
> >>>> Nevertheless I think it would be the best solution if one could
> provide P0 (approx A0) and SLEPc derives the preconditioner from this.
> Would this be hard to implement?
> >>>
> >>> This is what Barry’s suggestion is implementing. Don’t know why it
> doesn’t work with your Python operator though.
> >>>
> >>> Thanks,
> >>> Pierre
> >>>
> >>>> best wishes
> >>>> Florian
> >>>>
> >>>>
> >>>> On Sat, Feb 13, 2021 at 4:19 AM Barry Smith <bsmith at petsc.dev> wrote:
> >>>>
> >>>>
> >>>>> On Feb 12, 2021, at 2:32 AM, Florian Bruckner <e0425375 at gmail.com>
> wrote:
> >>>>>
> >>>>> Dear Jose, Dear Matt,
> >>>>>
> >>>>> I needed some time to think about your answers.
> >>>>> If I understand correctly, the eigenmode solver internally uses
> A0^{-1}*B0, which is normally handled by the ST object, which creates a KSP
> solver and a corresponding preconditioner.
> >>>>> What I would need is an interface to provide not only the system
> Matrix A0 (which is an operator), but also a preconditioning matrix (sparse
> approximation of the operator).
> >>>>> Unfortunately this interface is not available, right?
> >>>>
> >>>>    If SLEPc does not provide this directly it is still intended to be
> trivial to provide the "preconditioner matrix" (that is matrix from which
> the preconditioner is built). Just get the KSP from the ST object and use
> KSPSetOperators() to provide the "preconditioner matrix" .
> >>>>
> >>>>   Barry
> >>>>
> >>>>>
> >>>>> Matt directly creates A0^{-1}*B0 as a matshell operator. The
> operator uses a KSP with a proper PC internally. SLEPc would directly get
> A0^{-1}*B0 and solve a standard eigenvalue problem with this modified
> operator. Did I understand this correctly?
> >>>>>
> >>>>> I have two further points, which I did not mention yet: the matrix
> B0 is Hermitian, but it is (purely) imaginary (B0.real=0). Right now, I am
> using Firedrake to set up the PETSc system matrices A0, i*B0 (which is
> real). Then I convert them into ScipyLinearOperators and use
> scipy.sparse.eigsh(B0, b=A0, Minv=Minv) to calculate the eigenvalues.
> Minv=A0^-1 is also solving within scipy using a preconditioned gmres.
> Advantage of this setup is that the imaginary B0 can be handled efficiently
> and also the post-processing of the eigenvectors (which requires complex
> arithmetics) is simplified.
> >>>>>
> >>>>> Nevertheless I think that the mixing of PETSc and Scipy looks too
> complicated and is not very flexible.
> >>>>> If I would use Matt's approach, could I then simply switch between
> multiple standard eigenvalue methods (e.g. LOBPCG)? or is it limited due to
> the use of matshell?
> >>>>> Is there a solution for the imaginary B0, or do I have to use the
> non-hermitian methods? Is this a large performance drawback?
> >>>>>
> >>>>> thanks again,
> >>>>> and best wishes
> >>>>> Florian
> >>>>>
> >>>>> On Mon, Feb 8, 2021 at 3:37 PM Jose E. Roman <jroman at dsic.upv.es>
> wrote:
> >>>>> The problem can be written as A0*v=omega*B0*v and you want the
> eigenvalues omega closest to zero. If the matrices were explicitly
> available, you would do shift-and-invert with target=0, that is
> >>>>>
> >>>>>   (A0-sigma*B0)^{-1}*B0*v=theta*v    for sigma=0, that is
> >>>>>
> >>>>>   A0^{-1}*B0*v=theta*v
> >>>>>
> >>>>> and you compute EPS_LARGEST_MAGNITUDE eigenvalues theta=1/omega.
> >>>>>
> >>>>> Matt: I guess you should have EPS_LARGEST_MAGNITUDE instead of
> EPS_SMALLEST_REAL in your code. Are you getting the eigenvalues you need?
> EPS_SMALLEST_REAL will give slow convergence.
> >>>>>
> >>>>> Florian: I would not recommend setting the KSP matrices directly, it
> may produce strange side-effects. We should have an interface function to
> pass this matrix. Currently there is STPrecondSetMatForPC() but it has two
> problems: (1) it is intended for STPRECOND, so cannot be used with
> Krylov-Schur, and (2) it is not currently available in the python interface.
> >>>>>
> >>>>> The approach used by Matt is a workaround that does not use ST, so
> you can handle linear solves with a KSP of your own.
> >>>>>
> >>>>> As an alternative, since your problem is symmetric, you could try
> LOBPCG, assuming that the leftmost eigenvalues are those that you want
> (e.g. if all eigenvalues are non-negative). In that case you could use
> STPrecondSetMatForPC(), but the remaining issue is calling it from python.
> >>>>>
> >>>>> If you are using the git repo, I could add the relevant code.
> >>>>>
> >>>>> Jose
> >>>>>
> >>>>>
> >>>>>
> >>>>> > El 8 feb 2021, a las 14:22, Matthew Knepley <knepley at gmail.com>
> escribió:
> >>>>> >
> >>>>> > On Mon, Feb 8, 2021 at 7:04 AM Florian Bruckner <
> e0425375 at gmail.com> wrote:
> >>>>> > Dear PETSc / SLEPc Users,
> >>>>> >
> >>>>> > my question is very similar to the one posted here:
> >>>>> >
> https://lists.mcs.anl.gov/pipermail/petsc-users/2018-August/035878.html
> >>>>> >
> >>>>> > The eigensystem I would like to solve looks like:
> >>>>> > B0 v = 1/omega A0 v
> >>>>> > B0 and A0 are both hermitian, A0 is positive definite, but only
> given as a linear operator (matshell). I am looking for the largest
> eigenvalues (=smallest omega).
> >>>>> >
> >>>>> > I also have a sparse approximation P0 of the A0 operator, which i
> would like to use as precondtioner, using something like this:
> >>>>> >
> >>>>> >         es = SLEPc.EPS().create(comm=fd.COMM_WORLD)
> >>>>> >         st = es.getST()
> >>>>> >         ksp = st.getKSP()
> >>>>> >         ksp.setOperators(self.A0, self.P0)
> >>>>> >
> >>>>> > Unfortunately PETSc still complains that it cannot create a
> preconditioner for a type 'python' matrix although P0.type == 'seqaij' (but
> A0.type == 'python').
> >>>>> > By the way, should P0 be an approximation of A0 or does it have to
> include B0?
> >>>>> >
> >>>>> > Right now I am using the krylov-schur method. Are there any
> alternatives if A0 is only given as an operator?
> >>>>> >
> >>>>> > Jose can correct me if I say something wrong.
> >>>>> >
> >>>>> > When I did this, I made a shell operator for the action of A0^{-1}
> B0 which has a KSPSolve() in it, so you can use your P0 preconditioning
> matrix, and
> >>>>> > then handed that to EPS. You can see me do it here:
> >>>>> >
> >>>>> >
> https://gitlab.com/knepley/bamg/-/blob/master/src/coarse/bamgCoarseSpace.c#L123
> >>>>> >
> >>>>> > I had a hard time getting the embedded solver to work the way I
> wanted, but maybe that is the better way.
> >>>>> >
> >>>>> >   Thanks,
> >>>>> >
> >>>>> >      Matt
> >>>>> >
> >>>>> > thanks for any advice
> >>>>> > best wishes
> >>>>> > Florian
> >>>>> >
> >>>>> >
> >>>>> > --
> >>>>> > What most experimenters take for granted before they begin their
> experiments is infinitely more interesting than any results to which their
> experiments lead.
> >>>>> > -- Norbert Wiener
> >>>>> >
> >>>>> > https://www.cse.buffalo.edu/~knepley/
> >>>>>
> >>>>
> >>>
> >>
> >> <test.py>
> >
>
>

-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener

https://www.cse.buffalo.edu/~knepley/ <http://www.cse.buffalo.edu/~knepley/>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20210215/03201146/attachment.html>

More information about the petsc-users mailing list