[petsc-users] using preconditioner with SLEPc

[Jose E. Roman]

> El 12 feb 2021, a las 9:32, Florian Bruckner <e0425375 at gmail.com> escribió:
> 
> 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?

Yes, when using shift-and-invert with target=0 the solver internally uses A0^{-1}*B0. It also uses the B0-inner product to preserve symmetry.

> 
> 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?

Yes, the difference here is that you have to solve the problem as non-symmetric. This is not going to be slower, maybe just a bit less accurate. And the computed eigenvectors will not be B0-orthogonal.

> 
> 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?

LOBPCG can be used only in Hermitian problems, see Table 2.4 in the users manual.

The problem A0*v=omega*(i*B0)*v can be solved as A0*v=(i*omega)*B0*v, no problem with that. Just solve it as non-Hermitian. As I said above, the computational cost should not be too different. But solving as non-Hermitian restricts the available solvers.

Jose

> 
> 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/
>