[petsc-users] using preconditioner with SLEPc

Thu Feb 18 06:00:58 CST 2021

Dear Jose,
thanks for your work. I just looked over the code, but I didn't have time
to implement our solver, yet.
If I understand the code correctly, it allows to set a precond-matrix which
should approximate A-sigma*B.

I will try to get our code running in the next few weeks. From user
perspective it would maybe simplify things if approximations for A as well
as B are given, since this would hide the internal ST transformations.

best wishes
Florian

On Tue, Feb 16, 2021 at 8:54 PM Jose E. Roman <jroman at dsic.upv.es> wrote:

> Florian: I have created a MR
> https://gitlab.com/slepc/slepc/-/merge_requests/149
> Let me know if it fits your needs.
>
> Jose
>
>
> > El 15 feb 2021, a las 18:44, Jose E. Roman <jroman at dsic.upv.es>
> escribió:
> >
> >
> >
> >> El 15 feb 2021, a las 14:53, Matthew Knepley <knepley at gmail.com>
> escribió:
> >>
> >> 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?
> >
> > Yes, you can do the following:
> >
> >  ierr = EPSSetOperators(eps,S,NULL);CHKERRQ(ierr);   // S is your shell
> matrix A^{-1}*B
> >  ierr = EPSSetProblemType(eps,EPS_HEP);CHKERRQ(ierr);  // symmetric
> problem though S is not symmetric
> >  ierr = EPSSetFromOptions(eps);CHKERRQ(ierr);
> >  ierr = EPSSetUp(eps);CHKERRQ(ierr);   // note explicitly calling setup
> here
> >  ierr = EPSGetBV(eps,&bv);CHKERRQ(ierr);
> >  ierr = BVSetMatrix(bv,B,PETSC_FALSE);CHKERRQ(ierr);    // replace
> solver's inner product
> >  ierr = EPSSolve(eps);CHKERRQ(ierr);
> >
> > I have tried this with test1.c and it works. The computed eigenvectors
> should be B-orthogonal in this case.
> >
> > Jose
> >
> >
> >>
> >>  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/
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20210218/8b7a6465/attachment-0001.html>