You are missing a -ksp_reuse_preconditioner false when you want to compute a new factorization.
Wow, such a stupid mistake 
Now everything work perfectly, thank you!
Dear @prj ,
I would like to try out reusing LU factorization when solving a polynomial eigenvalue problem with PEPSolve in a similar way that you recommanded above. Something like first solving the problem with -st_type sinvert -st_pc_type lu, then slightly change the coefficient matrices, and use an iterative method also with reusing the previous factorization to precondition the new system. The naive way would be setting -st_ksp_type fgmres -st_pc_type lu -st_ksp_reuse_preconditioner; however, I suppose the internal KSP object created by PEPSolve (or EPSSolve) is destroyed after the call. Is there any easy way to do such a thing? If so, where should I start looking?
Thanks for the help in advance!
We could do exactly as for EPS (FreeFem-sources/SLEPc-code.hpp at develop · FreeFem/FreeFem-sources · GitHub), i.e., call PEPGetST(), and if there is a KSP attached to one of the Mat, use it via STSetKSP(). But I’m not sure of what you want to do exactly, so if you could share a small example (possibly in private), it would be easier for me to help you out.
Thanks for the quick reply. Based on the info you provided, I modified the Navier-Stokes SLEPc example, which hopefully illustrated what I have in mind. I seem to be able to access the KSP object. However, running the script with the -log_view option shows MatLUFactorNum 2, which indicates that the LU factorization is computed once again (or I also managed to crash the program).
I would really appreciate if you could take a look at the code.
navier-stokes-2d-SLEPc-complex.edp (3.4 KB)
You are not attaching any KSP information (set(J, sparams = "...");) to J prior to the first call to EPSSolve(), so it is expected that there are two numerical factorizations.
Interestingly, setting the same option for both the matrix J and EPSSolve before the first call does not change the fact that LU factorization is performed twice. Furthermore, I find it odd that even if I do not do set(J, sparams = "..."); before the first EPSSolve call, set(J, sparams = " -ksp_monitor"); changes the behavior of the program.
navier-stokes-2d-SLEPc-complex.edp (2.5 KB)
I don’t understand why the MUMPS instance is being destroyed at the end of the first EPSSolve(), I’ll need to investigate this in the upcoming days (busy schedule). I’ll keep you posted ASAP.
1 Like
Thank you very much!
This should be fixed in Fixes reuse of KSP in ST. · FreeFem/FreeFem-sources@6c1db49 · GitHub, could you please give it a go? Then, we will move to PEPSolve(), if that’s OK for you.
Thanks for the quick fix. I will try it out in the evening and let you know about the results.
I compiled FreeFem with the fix, it works perfectly. Not sure why, but setting only set(J, sparams = " -ksp_type fgmres -pc_type lu -ksp_reuse_preconditioner true "); achieves the desired effect, while changing only -st_ksp_type fgmres -st_pc_type lu -st_ksp_reuse_preconditioner true does not do the job (but I guess this does not matter - it works 
).
I would like it very much if we could try this out with PEPSolve. I think I could rewrite my polynomial eigenvalue problem with the companion matrix method into a normal eigenvalue problem, assemble the system matrix and use EPSSolve; however, I prefer SLEPc to handle all that.
Furthermore, I have some other questions which I think are important if I want to try out iterative solvers.
-
The following code suggests me that if I give
EPSSolve(andPEPSolve?) nonzero vectors, those are used as a guess for the followingEPSSolvecall? Because the SLEPc documentation suggests me that this might not be a good idea for the Krylov-Schur method. -
EPSSetDeflationSpacecan be accessed with passing the inputdeflation=...to the eigenvalue solver? How does this work - I am struggling to figure this out on my own.
Thank you very much for the help!
Not sure why, but setting only […]
The st_ prefix is added when the KSP is created by SLEPc. If you create it yourself with set(...), the option should not be prefixed, unless you add the flag set(..., prefix = "st_").
[…] however, I prefer SLEPc to handle all that.
Would it be possible to ask you for a piece of code that calls PEPSolve() and where I could try the -ksp_reuse_preconditioner option? You can share that in private if it’s best for you.
[…] nonzero vectors, those are used as a guess […]
I’ve never seen tremendous improvements using this feature.
How does this work - I am struggling to figure this out on my own.
Please see Slepc ESPSetDeflationSpace function andd let me know if you need more details.
Thanks for the detailed answers.
I slightly modified the blasius-stability-1d-SLEPc-complex.edp example with a second PEPSolve call - this is very similar to what I am trying to achieve, yet it is close to an MWE.
blasius-stability-1d-SLEPc-complex.edp (5.6 KB)
I think this works with Gets a ST from a PEP (stored in a Mat). · FreeFem/FreeFem-sources@cc9ba13 · GitHub and the updated
blasius-stability-1d-SLEPc-complex.edp (5.7 KB). I don’t like the code very much, so it’s currently in its own branch. If this works OK for you, I’ll try to tidy this up. Please let me know.
Thank you very much. I will test this feature on my problem tomorrow and let you know about the results I got.
Dear @prj ,
I tested reusing the LU factorization of the polynomial eigenvalue problem for my main problem. It works well; unfortunately, only a small speedup could be achieved and only in the case of using a small krylov subspace. Therefore it is unlikely that I will use this feature for my research. Nevertheless, thank you very much for the implementation!
Thanks for letting me know. I’ll babysit the pep branch for a while until we find out if it’s really not worth it.