<div dir="ltr"><div dir="ltr">On Tue, Mar 9, 2021 at 8:23 AM Florian Bruckner <<a href="mailto:e0425375@gmail.com">e0425375@gmail.com</a>> wrote:<br></div><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div>Dear Jose, <br></div><div><br></div><div>good news: something is running. I have to admit that I don't understand the difference of what you suggested, but if I run</div><div> es = SLEPc.EPS().create(comm=fd.COMM_WORLD) <br> es.setDimensions(k)<br> es.setType(SLEPc.EPS.Type.KRYLOVSCHUR) <br> es.setProblemType(SLEPc.EPS.ProblemType.GNHEP) # Generalized Non-Hermitian eigenproblem with positive definite B<br> es.setTarget(0.0)<br> es.setWhichEigenpairs(SLEPc.EPS.Which.TARGET_MAGNITUDE) <br> es.setTolerances(1e-10)<br> es.setOperators(self.A0, self.B0) <br> es.setFromOptions()<br> st = es.getST()<br> st.setType(SLEPc.ST.Type.SINVERT)<br> st.setPreconditionerMat(self.P0) <br></div><div> es.solve() <br> es.view()</div><div>python run_slepc.py -st_ksp_type gmres -st_pc_type bjacobi<br></div><div><br></div><div>i get the correct eigenvalues if as -1j*es.getEigenvalue(i) (the -1j is because B0 should be purely imaginary).</div><div>Omitting the preconditioning matrix gives the same results, but as expected is significantly slower.</div><div><br></div><div>If running it with the -st_ksp_type preonly -st_pc_type lu it still converges against totally wrong values.</div><div>But this could be due to the preonly option. If only the PC is applied only, it is clear that A0 is not used at all, right?<br></div><div><br></div><div></div><div>I didn't set LU manually. Perhaps firedrake changes the defaults somehow? But this should not be a problem.</div></div></blockquote><div><br></div><div>Hi Florian,</div><div><br></div><div>This is a mismatch of assumptions by you and Firedrake. The solver (preonly, LU) factors the _preconditioning_ matrix and</div><div>solves it directly. Usually this is very robust, but it assumes that the system matrix and preconditioning matrix are the same,</div><div>so you get the solution to your actual problem. Here you have an approximate preconditioning matrix, so (preonly, LU) solves</div><div>only that approximate problem. The solver (gmres, LU) will use GMRES to solve the system matrix, preconditioned by the LU</div><div>solve of the preconditioning matrix, which gives you the right result.</div><div><br></div><div>Does that make sense?</div><div><br></div><div> Thanks,</div><div><br></div><div> Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div>So, the code seems to work. Just to be sure, solving with (A0, B0) for SMALLEST_MAGNITUDE should be similar to target=0 and TARGET_MAGNITUDE, or is there any difference?</div><div>Finally, also solving (B0, A0) with LARGEST_MAGNITUDE should be similar, but one gets 1/omega as eigenvalues. In all cases the provided precond matrix should approximate A0, right?</div><div>Is this correct, or are the differences in the implementation when using the different formulations of the problem?</div><div><br></div><div>again, many many thanks.</div><div>best wishes <br></div><div>Florian <br></div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Tue, Mar 9, 2021 at 1:20 PM Jose E. Roman <<a href="mailto:jroman@dsic.upv.es" target="_blank">jroman@dsic.upv.es</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><br>
<br>
> El 9 mar 2021, a las 13:07, Florian Bruckner <<a href="mailto:e0425375@gmail.com" target="_blank">e0425375@gmail.com</a>> escribió:<br>
> <br>
> Dear Jose, <br>
> I appended the output of eps-view for the original method and for the method you proposed. <br>
> Unfortunately the new method does not converge. This I what i did:<br>
> <br>
> es = SLEPc.EPS().create(comm=fd.COMM_WORLD)<br>
> es.setDimensions(k)<br>
> es.setType(SLEPc.EPS.Type.KRYLOVSCHUR)<br>
> es.setProblemType(SLEPc.EPS.ProblemType.PGNHEP) # Generalized Non-Hermitian eigenproblem with positive definite B<br>
> es.setTarget(0.0)<br>
> es.setWhichEigenpairs(SLEPc.EPS.Which.TARGET_MAGNITUDE)<br>
> es.setTolerances(1e-10)<br>
> es.setOperators(self.B0, self.A0)<br>
> es.setFromOptions()<br>
> st = es.getST()<br>
> st.setPreconditionerMat(self.P0)<br>
<br>
No, you did not set SINVERT. Also, you should set (A,B) and not (B,A). And forget about PGNHEP for the moment.<br>
<br>
> <br>
> Is TARGET_MAGNITUDE correct? If I change it back to LARGEST_MAGNITUDE i can reproduce the (wrong) results from before. <br>
> Why do I need the shift-invert mode at all? I thought this is only necessary if I would like to solve for the smallest eigenmodes?<br>
> <br>
> Without st.setPreconditionerMat(self.P0) the code does not work, because A0 is a matshell and the preconditioner cannot be set up. <br>
> If I use pc_type = None the method converges, but results are totally wrong (1e-12 GHz instead of 7GHz).<br>
> <br>
> What confuses me most, is that the slepc results (without target=0 and with the precond matrix) produces nearly correct results, <br>
> which perfectly fit the results from scipy when using P0 instead of A0. This is a strange coincidence, and looks like P0 is used somewhere instead of A0.<br>
<br>
As I said, the problem is that it is using PREONLY+LU when it should use GMRES+BJACOBI by default. Are you setting PREONLY+LU somehow? Try running with -st_ksp_type gmres -st_pc_type bjacobi<br>
<br>
Jose<br>
<br>
> <br>
> thanks for your help<br>
> Florian <br>
> <br>
> On Tue, Mar 9, 2021 at 10:48 AM Jose E. Roman <<a href="mailto:jroman@dsic.upv.es" target="_blank">jroman@dsic.upv.es</a>> wrote:<br>
> The reason may be that it is using a direct solver instead of an iterative solver. What do you get for -eps_view ?<br>
> <br>
> Does the code work correctly if you comment out st.setPreconditionerMat(self.P0) ?<br>
> <br>
> Your approach should work, but I would first try as is done in the example <a href="https://slepc.upv.es/slepc-main/src/eps/tutorials/ex46.c.html" rel="noreferrer" target="_blank">https://slepc.upv.es/slepc-main/src/eps/tutorials/ex46.c.html</a><br>
> that is, shift-and-invert with target=0 and target_magnitude.<br>
> <br>
> Jose<br>
> <br>
> <br>
> > El 9 mar 2021, a las 9:07, Florian Bruckner <<a href="mailto:e0425375@gmail.com" target="_blank">e0425375@gmail.com</a>> escribió:<br>
> > <br>
> > Dear Jose, <br>
> > I asked Lawrence Mitchell from the firedrake people to help me with the slepc update (I think they are applying some modifications for petsc, which is why simply updating petsc within my docker container did not work).<br>
> > Now the latest slepc version runs and I already get some results of the eigenmode solver. The good thing is that the solver runs significantly faster. The bad thing is that the results are still wrong :-)<br>
> > <br>
> > Could you have a short look at the code:<br>
> > es = SLEPc.EPS().create(comm=fd.COMM_WORLD)<br>
> > es.setDimensions(k)<br>
> > es.setType(SLEPc.EPS.Type.KRYLOVSCHUR)<br>
> > es.setProblemType(SLEPc.EPS.ProblemType.PGNHEP) # Generalized Non-Hermitian eigenproblem with positive definite B<br>
> > es.setWhichEigenpairs(SLEPc.EPS.Which.LARGEST_MAGNITUDE)<br>
> > #es.setTrueResidual(True)<br>
> > es.setTolerances(1e-10)<br>
> > es.setOperators(self.B0, self.A0)<br>
> > es.setFromOptions()<br>
> > <br>
> > st = es.getST()<br>
> > st.setPreconditionerMat(self.P0)<br>
> > <br>
> > You wrote that when using shift-and-invert with target=0 the solver internally uses A0^{-1}*B0. <br>
> > Nevertheless I think the precond P0 mat should be an approximation of A0, right? <br>
> > This is because the solver uses the B0-inner product to preserve symmetry. <br>
> > Or is the B0 missing in my code?<br>
> > <br>
> > As I mentioned before, convergence of the method is extremely fast. I thought that maybe the tolerance is set too low, but increasing it did not change the situation. <br>
> > With using setTrueResidual, there is no convergence at all.<br>
> > <br>
> > Figures show the different results for the original scipy method (which has been compared to published results) as well as the new slepc method.<br>
> > For some strange reason I get nearly the same (wrong) results if i replace A0 with P0 in the original scipy code. <br>
> > In my case A0 is a non-local field operator and P0 only contains local and next-neighbour interaction. <br>
> > Is it possible that the wrong operator (P0 instead of A0) is used internally?<br>
> > <br>
> > best wishes<br>
> > Florian<br>
> > <br>
> > On Thu, Feb 18, 2021 at 1:00 PM Florian Bruckner <<a href="mailto:e0425375@gmail.com" target="_blank">e0425375@gmail.com</a>> wrote:<br>
> > Dear Jose,<br>
> > thanks for your work. I just looked over the code, but I didn't have time to implement our solver, yet. <br>
> > If I understand the code correctly, it allows to set a precond-matrix which should approximate A-sigma*B.<br>
> > <br>
> > 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.<br>
> > <br>
> > best wishes<br>
> > Florian <br>
> > <br>
> > On Tue, Feb 16, 2021 at 8:54 PM Jose E. Roman <<a href="mailto:jroman@dsic.upv.es" target="_blank">jroman@dsic.upv.es</a>> wrote:<br>
> > Florian: I have created a MR <a href="https://gitlab.com/slepc/slepc/-/merge_requests/149" rel="noreferrer" target="_blank">https://gitlab.com/slepc/slepc/-/merge_requests/149</a><br>
> > Let me know if it fits your needs.<br>
> > <br>
> > Jose<br>
> > <br>
> > <br>
> > > El 15 feb 2021, a las 18:44, Jose E. Roman <<a href="mailto:jroman@dsic.upv.es" target="_blank">jroman@dsic.upv.es</a>> escribió:<br>
> > > <br>
> > > <br>
> > > <br>
> > >> El 15 feb 2021, a las 14:53, Matthew Knepley <<a href="mailto:knepley@gmail.com" target="_blank">knepley@gmail.com</a>> escribió:<br>
> > >> <br>
> > >> On Mon, Feb 15, 2021 at 7:27 AM Jose E. Roman <<a href="mailto:jroman@dsic.upv.es" target="_blank">jroman@dsic.upv.es</a>> wrote:<br>
> > >> I will think about the viability of adding an interface function to pass the preconditioner matrix.<br>
> > >> <br>
> > >> 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.<br>
> > >> <br>
> > >> Yes, the vectors I generate are not B-orthogonal.<br>
> > >> <br>
> > >> 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?<br>
> > > <br>
> > > Yes, you can do the following:<br>
> > > <br>
> > > ierr = EPSSetOperators(eps,S,NULL);CHKERRQ(ierr); // S is your shell matrix A^{-1}*B<br>
> > > ierr = EPSSetProblemType(eps,EPS_HEP);CHKERRQ(ierr); // symmetric problem though S is not symmetric<br>
> > > ierr = EPSSetFromOptions(eps);CHKERRQ(ierr);<br>
> > > ierr = EPSSetUp(eps);CHKERRQ(ierr); // note explicitly calling setup here<br>
> > > ierr = EPSGetBV(eps,&bv);CHKERRQ(ierr);<br>
> > > ierr = BVSetMatrix(bv,B,PETSC_FALSE);CHKERRQ(ierr); // replace solver's inner product<br>
> > > ierr = EPSSolve(eps);CHKERRQ(ierr);<br>
> > > <br>
> > > I have tried this with test1.c and it works. The computed eigenvectors should be B-orthogonal in this case.<br>
> > > <br>
> > > Jose<br>
> > > <br>
> > > <br>
> > >> <br>
> > >> Thanks,<br>
> > >> <br>
> > >> Matt<br>
> > >> <br>
> > >> Jose<br>
> > >> <br>
> > >> <br>
> > >>> El 14 feb 2021, a las 21:41, Barry Smith <<a href="mailto:bsmith@petsc.dev" target="_blank">bsmith@petsc.dev</a>> escribió:<br>
> > >>> <br>
> > >>> <br>
> > >>> Florian,<br>
> > >>> <br>
> > >>> I'm sorry I don't know the answers; I can only speculate. There is a STGetShift(). <br>
> > >>> <br>
> > >>> All I was saying is theoretically there could/should be such support in SLEPc.<br>
> > >>> <br>
> > >>> Barry<br>
> > >>> <br>
> > >>> <br>
> > >>>> On Feb 13, 2021, at 6:43 PM, Florian Bruckner <<a href="mailto:e0425375@gmail.com" target="_blank">e0425375@gmail.com</a>> wrote:<br>
> > >>>> <br>
> > >>>> Dear Barry,<br>
> > >>>> 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.<br>
> > >>>> 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.<br>
> > >>>> <br>
> > >>>> 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?<br>
> > >>>> I wrote a small standalone example based on ex9.py from slepc4py, where i tried to use an operator.<br>
> > >>>> <br>
> > >>>> best wishes<br>
> > >>>> Florian<br>
> > >>>> <br>
> > >>>> On Sat, Feb 13, 2021 at 7:15 PM Barry Smith <<a href="mailto:bsmith@petsc.dev" target="_blank">bsmith@petsc.dev</a>> wrote:<br>
> > >>>> <br>
> > >>>> <br>
> > >>>>> On Feb 13, 2021, at 2:47 AM, Pierre Jolivet <<a href="mailto:pierre@joliv.et" target="_blank">pierre@joliv.et</a>> wrote:<br>
> > >>>>> <br>
> > >>>>> <br>
> > >>>>> <br>
> > >>>>>> On 13 Feb 2021, at 7:25 AM, Florian Bruckner <<a href="mailto:e0425375@gmail.com" target="_blank">e0425375@gmail.com</a>> wrote:<br>
> > >>>>>> <br>
> > >>>>>> Dear Jose, Dear Barry,<br>
> > >>>>>> 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?<br>
> > >>>>>> 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?<br>
> > >>>>>> <br>
> > >>>>>> Regarding Barry's advice this is what I first tried:<br>
> > >>>>>> es = SLEPc.EPS().create(comm=fd.COMM_WORLD)<br>
> > >>>>>> st = es.getST()<br>
> > >>>>>> ksp = st.getKSP()<br>
> > >>>>>> ksp.setOperators(self.A0, self.P0)<br>
> > >>>>>> <br>
> > >>>>>> 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?<br>
> > >>>>> <br>
> > >>>>> 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.<br>
> > >>>> <br>
> > >>>> Just trying to provide more clarity with the terms.<br>
> > >>>> <br>
> > >>>> 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<br>
> > >>>> <br>
> > >>>> 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".<br>
> > >>>> <br>
> > >>>> Barry<br>
> > >>>> <br>
> > >>>> <br>
> > >>>>> <br>
> > >>>>>> 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?<br>
> > >>>>> <br>
> > >>>>> This is what Barry’s suggestion is implementing. Don’t know why it doesn’t work with your Python operator though.<br>
> > >>>>> <br>
> > >>>>> Thanks,<br>
> > >>>>> Pierre<br>
> > >>>>> <br>
> > >>>>>> best wishes<br>
> > >>>>>> Florian<br>
> > >>>>>> <br>
> > >>>>>> <br>
> > >>>>>> On Sat, Feb 13, 2021 at 4:19 AM Barry Smith <<a href="mailto:bsmith@petsc.dev" target="_blank">bsmith@petsc.dev</a>> wrote:<br>
> > >>>>>> <br>
> > >>>>>> <br>
> > >>>>>>> On Feb 12, 2021, at 2:32 AM, Florian Bruckner <<a href="mailto:e0425375@gmail.com" target="_blank">e0425375@gmail.com</a>> wrote:<br>
> > >>>>>>> <br>
> > >>>>>>> Dear Jose, Dear Matt, <br>
> > >>>>>>> <br>
> > >>>>>>> I needed some time to think about your answers. <br>
> > >>>>>>> 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.<br>
> > >>>>>>> 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).<br>
> > >>>>>>> Unfortunately this interface is not available, right?<br>
> > >>>>>> <br>
> > >>>>>> 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" . <br>
> > >>>>>> <br>
> > >>>>>> Barry<br>
> > >>>>>> <br>
> > >>>>>>> <br>
> > >>>>>>> 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?<br>
> > >>>>>>> <br>
> > >>>>>>> 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. <br>
> > >>>>>>> <br>
> > >>>>>>> Nevertheless I think that the mixing of PETSc and Scipy looks too complicated and is not very flexible. <br>
> > >>>>>>> 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?<br>
> > >>>>>>> Is there a solution for the imaginary B0, or do I have to use the non-hermitian methods? Is this a large performance drawback?<br>
> > >>>>>>> <br>
> > >>>>>>> thanks again,<br>
> > >>>>>>> and best wishes<br>
> > >>>>>>> Florian<br>
> > >>>>>>> <br>
> > >>>>>>> On Mon, Feb 8, 2021 at 3:37 PM Jose E. Roman <<a href="mailto:jroman@dsic.upv.es" target="_blank">jroman@dsic.upv.es</a>> wrote:<br>
> > >>>>>>> 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<br>
> > >>>>>>> <br>
> > >>>>>>> (A0-sigma*B0)^{-1}*B0*v=theta*v for sigma=0, that is<br>
> > >>>>>>> <br>
> > >>>>>>> A0^{-1}*B0*v=theta*v<br>
> > >>>>>>> <br>
> > >>>>>>> and you compute EPS_LARGEST_MAGNITUDE eigenvalues theta=1/omega.<br>
> > >>>>>>> <br>
> > >>>>>>> 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.<br>
> > >>>>>>> <br>
> > >>>>>>> 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.<br>
> > >>>>>>> <br>
> > >>>>>>> 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.<br>
> > >>>>>>> <br>
> > >>>>>>> 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.<br>
> > >>>>>>> <br>
> > >>>>>>> If you are using the git repo, I could add the relevant code.<br>
> > >>>>>>> <br>
> > >>>>>>> Jose<br>
> > >>>>>>> <br>
> > >>>>>>> <br>
> > >>>>>>> <br>
> > >>>>>>>> El 8 feb 2021, a las 14:22, Matthew Knepley <<a href="mailto:knepley@gmail.com" target="_blank">knepley@gmail.com</a>> escribió:<br>
> > >>>>>>>> <br>
> > >>>>>>>> On Mon, Feb 8, 2021 at 7:04 AM Florian Bruckner <<a href="mailto:e0425375@gmail.com" target="_blank">e0425375@gmail.com</a>> wrote:<br>
> > >>>>>>>> Dear PETSc / SLEPc Users,<br>
> > >>>>>>>> <br>
> > >>>>>>>> my question is very similar to the one posted here: <br>
> > >>>>>>>> <a href="https://lists.mcs.anl.gov/pipermail/petsc-users/2018-August/035878.html" rel="noreferrer" target="_blank">https://lists.mcs.anl.gov/pipermail/petsc-users/2018-August/035878.html</a><br>
> > >>>>>>>> <br>
> > >>>>>>>> The eigensystem I would like to solve looks like:<br>
> > >>>>>>>> B0 v = 1/omega A0 v<br>
> > >>>>>>>> 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). <br>
> > >>>>>>>> <br>
> > >>>>>>>> I also have a sparse approximation P0 of the A0 operator, which i would like to use as precondtioner, using something like this:<br>
> > >>>>>>>> <br>
> > >>>>>>>> es = SLEPc.EPS().create(comm=fd.COMM_WORLD)<br>
> > >>>>>>>> st = es.getST()<br>
> > >>>>>>>> ksp = st.getKSP()<br>
> > >>>>>>>> ksp.setOperators(self.A0, self.P0)<br>
> > >>>>>>>> <br>
> > >>>>>>>> Unfortunately PETSc still complains that it cannot create a preconditioner for a type 'python' matrix although P0.type == 'seqaij' (but A0.type == 'python'). <br>
> > >>>>>>>> By the way, should P0 be an approximation of A0 or does it have to include B0?<br>
> > >>>>>>>> <br>
> > >>>>>>>> Right now I am using the krylov-schur method. Are there any alternatives if A0 is only given as an operator?<br>
> > >>>>>>>> <br>
> > >>>>>>>> Jose can correct me if I say something wrong.<br>
> > >>>>>>>> <br>
> > >>>>>>>> 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<br>
> > >>>>>>>> then handed that to EPS. You can see me do it here:<br>
> > >>>>>>>> <br>
> > >>>>>>>> <a href="https://gitlab.com/knepley/bamg/-/blob/master/src/coarse/bamgCoarseSpace.c#L123" rel="noreferrer" target="_blank">https://gitlab.com/knepley/bamg/-/blob/master/src/coarse/bamgCoarseSpace.c#L123</a><br>
> > >>>>>>>> <br>
> > >>>>>>>> I had a hard time getting the embedded solver to work the way I wanted, but maybe that is the better way.<br>
> > >>>>>>>> <br>
> > >>>>>>>> Thanks,<br>
> > >>>>>>>> <br>
> > >>>>>>>> Matt<br>
> > >>>>>>>> <br>
> > >>>>>>>> thanks for any advice<br>
> > >>>>>>>> best wishes<br>
> > >>>>>>>> Florian<br>
> > >>>>>>>> <br>
> > >>>>>>>> <br>
> > >>>>>>>> -- <br>
> > >>>>>>>> What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>
> > >>>>>>>> -- Norbert Wiener<br>
> > >>>>>>>> <br>
> > >>>>>>>> <a href="https://www.cse.buffalo.edu/~knepley/" rel="noreferrer" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br>
> > >>>>>>> <br>
> > >>>>>> <br>
> > >>>>> <br>
> > >>>> <br>
> > >>>> <test.py><br>
> > >>> <br>
> > >> <br>
> > >> <br>
> > >> <br>
> > >> -- <br>
> > >> What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>
> > >> -- Norbert Wiener<br>
> > >> <br>
> > >> <a href="https://www.cse.buffalo.edu/~knepley/" rel="noreferrer" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br>
> > <br>
> > <frequencies_scipy_P0_instead_A0.png><frequencies_slepc.png><frequencies_scipy.png><br>
> <br>
> <epsview_original.txt><epsview_target0.txt><br>
<br>
</blockquote></div>
</blockquote></div><br clear="all"><div><br></div>-- <br><div dir="ltr" class="gmail_signature"><div dir="ltr"><div><div dir="ltr"><div><div dir="ltr"><div>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>-- Norbert Wiener</div><div><br></div><div><a href="http://www.cse.buffalo.edu/~knepley/" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br></div></div></div></div></div></div></div></div>