I am using FreeFem to solve Polynomial EVP with SLEPc PEPSolve. When I use the option -pep_monitor_all, the program displays more eigenvalues than the converged ones. As an example, when running the following example , the number of converged eigenvalues might be smaller then the ones thats residual meets the convergence criterion. Stating a bit differently, there are eigenvalues that are not listed among the converged ones but according to the monitor meet the convergence criterion. I would like to get all the eigenvalues and eigenvectors that meet the convergence criterion according to the monitor, as I can drop the unwanted ones based on physical reasoning.
Has anyone experienced similar behavior? I have a vague idea what can cause this: according the to SLEPc manual, since a spectral-transformation is used, the residual estimate might be optimistic. However, I have no idea how to get all the eigenvalues I want: this behaviour seems to persist even if I use PEPSetConvergenceTest with different options.
Well, I could give it a try although I must admit I have no experience with running SLEPc on its own.
I immediately face a question: I do not know how to run the examples. In the directory petsc/arch-FreeFem-complex/externalpackages/git.slepc/arch-FreeFem-complex/tests/pep/tutorials/nlevp
there are several .sh files, but they do not seem to work. On the other hand, in the directory petsc/arch-FreeFem-complex/share/slepc/examples/src/pep/tutorials/nlevp
there are the files which are mentioned in this tutorial. However, I do not know to what value should I set SLEPC_DIR, PETSC_DIR and PETSC_ARCH.
If you could can help me out how I can run SLEPC on its own, I can try to modify one of the examples to investigate the strange error behaviour.
Also, very trivial check, but do you use the proper -pep_nev? If you are using a value smaller than the number of converged eigenpairs, SLEPc will automatically discard the ncv-nev eigenpairs furthest from your target.
Good point regarding the -pep_nev as it can. However, in the Blasius example, whether I run the code with -pep_basis monomial -pep_general -st_type sinvert -st_pc_type lu -pep_monitor_all -pep_true_residual true -pep_target 0.4+0.1i -pep_nev 100 -pep_ncv 300 or -pep_nev 170, at the end of the solution the monitor reports nconv=170. The monitor also displays and 300 eigenvalues with errors, and actually 199 meets the convergence criterion according to the errors, not 170.
Also a quick note: in my other problem (from which I cannot send and MWE as it is very long and relies on other codes) using the more accurate but more costly -pep_conv_norm actually helps a bit, but in the Blasius example it does not make a significant difference (nconv=172 and actually 204 eigenvalues meet the criterion).
Nevertheless, this merely poses a new question: is there a way to access not only the converged eigenvectors/eigenvalues? I think I can set, based on a physical reasoning, a much more computationally economic convergence criterion (meaning, I do not have to do repeated PEP calculations over and over).
Browsing the SLEPc documentation I found that this would be possible with the -pep_conv_user or -eps_conv_user option; however, my guess is it is not yet interfaced in FreeFem. Is there a chance you could interface it in the near future?
Sorry for the late answer, but I got the following answer from Jose.
I think now I understand what the user is concerned about: the first 171 residuals are smaller than the tolerance, and the 172th one is larger, but looking at residuals after that there are also many with a residual below the tolerance. However, SLEPc solvers stop checking for residuals as soon as one of them does not satisfy the tolerance (they are checked in order). It is done in this way because accepting a subsequent eigenpair with a small residual does not guarantee that the actual residual will be of that order of magnitude. If the user wants more eigenpairs, he should increase the value of nev, and then the solver will do one iteration more.
I could (or you could!?) interface EPSSetConvergenceTestFunction() and PEPSetConvergenceTestFunction(), it’s mostly a copy/paste of what’s already in PETSc-code.hpp for, e.g., KSPMonitorSet().
I tried to figure out how to implement EPSSetConvergenceTestFunction. In the attachment you can find SLEPC-code.hpp I modified based on your suggestion. The new pieces of code are marked with newcode comments. Could you please take a look at it? I would really appreciate some feedback - as I do not know C++ at all, I have no idea whether I am working in the right direction or not. (I uploaded a .idp file as .hpp files cannot be uploaded) SLEPc-code-dev.idp (33.5 KB)
Thanks! My idea is illustrated in the attached script.
However, as I am thinking about it, it might be better approach to use a lower tolerance, interface EPSComputeError, and select the eigenvalues of interest after EPSSolve (or PEPSolve) is finished. Nevertheless, I would like to implement the convergence criterion functions and try them out. Hopefully I am not miles away from the solution, and if EPSSetConvergenceTestFunction does not work out for my problem, implementing EPSComputeError will be not that hard Blasius-stability-1d-PEP-SLEPc-complex.edp (7.7 KB)
I added EPSGetErrorEstimate (also PEPGetErrorEstimate/NEPGetErrorEstimate) as this was the easiest to implement, yet it may be enough to solve my problem as I can lower the eigenvalue tolerance, and filter the eigenvalues myself. Please check out this pull request; any comments are appreciated.
Sorry for the stupid question.
