SLEPc convergence and error

Dear FreeFem users,

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.

Can you reproduce this on a minimal working example in plain C or C++? I’m not sure if FreeFEM or its interface to SLEPc is at fault to be honest.

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
there are several .sh files, but they do not seem to work. On the other hand, in the directory
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.

OK, before doing that, I’ll try to see if Jose Roman could provide a quick answer to your original question.

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.

Thanks for the quick reply.

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).

I think I figured it out: within the convergence radius, all eigenvalues need to meet the convergence criterion :person_facepalming: Sorry for the stupid question.

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().

Thanks for the reply regarding my original question.

I really wish I could do such thing without causing you extra work. However, I do not know c++ at all. Here is what I got so far:

  • I need to define a new input like here.

  • I need to define Monitor and MonitorDestroy in the namespace and this class in the eigensolver template.

  • The key part is this piece of code. I think this should be added roughly here and here with eigensolver instead of LinearSolver.

I appreciate any help either with me coding or you adding these features yourself.

Dear @prj ,

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)

I’ll have a look at this shortly. In the meantime, could you please also send a FreeFEM .edp doing what you’d like to do?

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 :slight_smile:
Blasius-stability-1d-PEP-SLEPc-complex.edp (7.7 KB)

@prj ,

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.

Nice, I’ll check this out!

Now that your feature branch has been merged, could you please let me know if this provides a good enough workaround for your initial issue?

Sure think, I will let you know regarding my findings.