[petsc-users] using preconditioner with SLEPc

Fri Feb 12 02:32:26 CST 2021

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?

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/
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20210212/3041f3cf/attachment.html>