[petsc-users] using preconditioner with SLEPc

Sat Feb 13 18:43:47 CST 2021

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/
>>>
>>>
>>
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20210214/a9715c96/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: test.py
Type: text/x-python
Size: 2509 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20210214/a9715c96/attachment-0001.py>