Hi Jose, I have now gotten the quadratic problem working decently using the PEP package with appropriate scaling and preconditioning, so thanks for all the suggestions! For the case where K is a shell matrix, I used a scaling based on an approximation of K, and that seems to be working well.
So now that both linear and quadratic problems are working, I wanted to get your suggestions on solving a non-linear problem. In some of our cases, we have a non-linear source term S(lambda) on the right-hand side of the equation as follows: (K + lambda*C + lambda^2*M)*x = S(lambda)*x, where the source can sometimes be simplified as S(lambda) = exp(lambda*t)*A, where A is a constant matrix. I am currently solving this non-linear problem iteratively. For each eigenvalue, I compute the source and add it into the K matrix, and then iterate until convergence. For this reason, I end up solving the system multiple times which makes it very slow. I saw some examples of non-linear problems included in the NEP package. I just wanted to get your thoughts if I would benefit from using the NEP package for this particular problem? Will I be able to use preconditioning and scaling as with the PEP package to speed up the computation for the case where K is a shell matrix? Thanks for your help. Regards, Varun On Thu, Sep 30, 2021 at 10:12 PM Varun Hiremath <varunhirem...@gmail.com> wrote: > Hi Jose, > > Thanks again for your valuable suggestions. I am still working on this but > wanted to give you a quick update. > > For the linear problem, I tried different KSP solvers, and finally, I'm > getting good convergence using CGS with LU (using MUMPS) inexact inverse. > So thank you very much for your help! > > But for the quadratic problem, I'm still struggling. As you suggested, I > have now started using the PEP solver. For the simple case where the K > matrix is explicitly known, everything works fine. But for the case where K > is a shell matrix, it struggles to converge. I am yet to try the scaling > option and some other preconditioning options. I will get back to you on > this if I have any questions. Appreciate your help! > > Thanks, > Varun > > On Tue, Sep 28, 2021 at 8:09 AM Jose E. Roman <jro...@dsic.upv.es> wrote: > >> >> >> > El 28 sept 2021, a las 7:50, Varun Hiremath <varunhirem...@gmail.com> >> escribió: >> > >> > Hi Jose, >> > >> > I implemented the LU factorized preconditioner and tested it using >> PREONLY + LU, but that actually is converging to the wrong eigenvalues, >> compared to just using BICGS + BJACOBI, or simply computing >> EPS_SMALLEST_MAGNITUDE without any preconditioning. My preconditioning >> matrix is only a 1st order approximation, and the off-diagonal terms are >> not very accurate, so I'm guessing this is why the LU factorization doesn't >> help much? Nonetheless, using BICGS + BJACOBI with slightly relaxed >> tolerances seems to be working fine. >> >> If your PCMAT is not an exact inverse, then you have to iterate, i.e. not >> use KSPPREONLY but KSPBCGS or another. >> >> > >> > I now want to test the same preconditioning idea for a quadratic >> problem. I am solving a quadratic equation similar to Eqn.(5.1) in the >> SLEPc manual: >> > (K + lambda*C + lambda^2*M)*x = 0, >> > I don't use the PEP package directly, but solve this by linearizing >> similar to Eqn.(5.3) and calling EPS. Without explicitly forming the full >> matrix, I just use the block matrix structure as explained in the below >> example and that works nicely for my case: >> > https://slepc.upv.es/documentation/current/src/eps/tutorials/ex9.c.html >> >> Using PEP is generally recommended. The default solver TOAR is >> memory-efficient and performs less computation than a trivial >> linearization. In addition, PEP allows you to do scaling, which is often >> very important to get accurate results in some problems, depending on >> conditioning. >> >> In your case K is a shell matrix, so things may not be trivial. If I am >> not wrong, you should be able to use STSetPreconditionerMat() for a PEP, >> where the preconditioner in this case should be built to approximate >> Q(sigma), where Q(.) is the quadratic polynomial and sigma is the target. >> >> > >> > In my case, K is not explicitly known, and for linear problems, where C >> = 0, I am using a 1st order approximation of K as the preconditioner. Now >> could you please tell me if there is a way to conveniently set the >> preconditioner for the quadratic problem, which will be of the form [-K 0; >> 0 I]? Note that K is constructed in parallel (the rows are distributed), so >> I wasn't sure how to construct this preconditioner matrix which will be >> compatible with the shell matrix structure that I'm using to define the >> MatMult function as in ex9. >> >> The shell matrix of ex9.c interleaves the local parts of the first block >> and the second block. In other words, a process' local part consists of the >> local rows of the first block followed by the local rows of the second >> block. In your case, the local rows of K followed by the local rows of the >> identity (appropriately padded with zeros). >> >> Jose >> >> >> > >> > Thanks, >> > Varun >> > >> > On Fri, Sep 24, 2021 at 11:50 PM Varun Hiremath < >> varunhirem...@gmail.com> wrote: >> > Ok, great! I will give that a try, thanks for your help! >> > >> > On Fri, Sep 24, 2021 at 11:12 PM Jose E. Roman <jro...@dsic.upv.es> >> wrote: >> > Yes, you can use PCMAT >> https://petsc.org/release/docs/manualpages/PC/PCMAT.html then pass a >> preconditioner matrix that performs the inverse via a shell matrix. >> > >> > > El 25 sept 2021, a las 8:07, Varun Hiremath <varunhirem...@gmail.com> >> escribió: >> > > >> > > Hi Jose, >> > > >> > > Thanks for checking my code and providing suggestions. >> > > >> > > In my particular case, I don't know the matrix A explicitly, I >> compute A*x in a matrix-free way within a shell matrix, so I can't use any >> of the direct factorization methods. But just a question regarding your >> suggestion to compute a (parallel) LU factorization. In our work, we do use >> MUMPS to compute the parallel factorization. For solving the generalized >> problem, A*x = lambda*B*x, we are computing inv(B)*A*x within a shell >> matrix, where factorization of B is computed using MUMPS. (We don't call >> MUMPS through SLEPc as we have our own MPI wrapper and other user settings >> to handle.) >> > > >> > > So for the preconditioning, instead of using the iterative solvers, >> can I provide a shell matrix that computes inv(P)*x corrections (where P is >> the preconditioner matrix) using MUMPS direct solver? >> > > >> > > And yes, thanks, #define PETSC_USE_COMPLEX 1 is not needed, it works >> without it. >> > > >> > > Regards, >> > > Varun >> > > >> > > On Fri, Sep 24, 2021 at 9:14 AM Jose E. Roman <jro...@dsic.upv.es> >> wrote: >> > > If you do >> > > $ ./acoustic_matrix_test.o -shell 0 -st_type sinvert -deflate 1 >> > > then it is using an LU factorization (the default), which is fast. >> > > >> > > Use -eps_view to see which solver settings are you using. >> > > >> > > BiCGStab with block Jacobi does not work for you matrix, it exceeds >> the maximum 10000 iterations. So this is not viable unless you can find a >> better preconditioner for your problem. If not, just using >> EPS_SMALLEST_MAGNITUDE will be faster. >> > > >> > > Computing smallest magnitude eigenvalues is a difficult task. The >> most robust way is to compute a (parallel) LU factorization if you can >> afford it. >> > > >> > > >> > > A side note: don't add this to your source code >> > > #define PETSC_USE_COMPLEX 1 >> > > This define is taken from PETSc's include files, you should not mess >> with it. Instead, you probably want to add something like this AFTER >> #include <slepceps.h>: >> > > #if !defined(PETSC_USE_COMPLEX) >> > > #error "Requires complex scalars" >> > > #endif >> > > >> > > Jose >> > > >> > > >> > > > El 22 sept 2021, a las 19:38, Varun Hiremath < >> varunhirem...@gmail.com> escribió: >> > > > >> > > > Hi Jose, >> > > > >> > > > Thank you, that explains it and my example code works now without >> specifying "-eps_target 0" in the command line. >> > > > >> > > > However, both the Krylov inexact shift-invert and JD solvers are >> struggling to converge for some of my actual problems. The issue seems to >> be related to non-symmetric general matrices. I have extracted one such >> matrix attached here as MatA.gz (size 100k), and have also included a short >> program that loads this matrix and then computes the smallest eigenvalues >> as I described earlier. >> > > > >> > > > For this matrix, if I compute the eigenvalues directly (without >> using the shell matrix) using shift-and-invert (as below) then it converges >> in less than a minute. >> > > > $ ./acoustic_matrix_test.o -shell 0 -st_type sinvert -deflate 1 >> > > > >> > > > However, if I use the shell matrix and use any of the >> preconditioned solvers JD or Krylov shift-invert (as shown below) with the >> same matrix as the preconditioner, then they struggle to converge. >> > > > $ ./acoustic_matrix_test.o -usejd 1 -deflate 1 >> > > > $ ./acoustic_matrix_test.o -sinvert 1 -deflate 1 >> > > > >> > > > Could you please check the attached code and suggest any changes in >> settings that might help with convergence for these kinds of matrices? I >> appreciate your help! >> > > > >> > > > Thanks, >> > > > Varun >> > > > >> > > > On Tue, Sep 21, 2021 at 11:14 AM Jose E. Roman <jro...@dsic.upv.es> >> wrote: >> > > > I will have a look at your code when I have more time. Meanwhile, I >> am answering 3) below... >> > > > >> > > > > El 21 sept 2021, a las 0:23, Varun Hiremath < >> varunhirem...@gmail.com> escribió: >> > > > > >> > > > > Hi Jose, >> > > > > >> > > > > Sorry, it took me a while to test these settings in the new >> builds. I am getting good improvement in performance using the >> preconditioned solvers, so thanks for the suggestions! But I have some >> questions related to the usage. >> > > > > >> > > > > We are using SLEPc to solve the acoustic modal eigenvalue >> problem. Attached is a simple standalone program that computes acoustic >> modes in a simple rectangular box. This program illustrates the general >> setup I am using, though here the shell matrix and the preconditioner >> matrix are the same, while in my actual program the shell matrix computes >> A*x without explicitly forming A, and the preconditioner is a 0th order >> approximation of A. >> > > > > >> > > > > In the attached program I have tested both >> > > > > 1) the Krylov-Schur with inexact shift-and-invert (implemented >> under the option sinvert); >> > > > > 2) the JD solver with preconditioner (implemented under the >> option usejd) >> > > > > >> > > > > Both the solvers seem to work decently, compared to no >> preconditioning. This is how I run the two solvers (for a mesh size of >> 1600x400): >> > > > > $ ./acoustic_box_test.o -nx 1600 -ny 400 -usejd 1 -deflate 1 >> -eps_target 0 >> > > > > $ ./acoustic_box_test.o -nx 1600 -ny 400 -sinvert 1 -deflate 1 >> -eps_target 0 >> > > > > Both finish in about ~10 minutes on my system in serial. JD seems >> to be slightly faster and more accurate (for the imaginary part of >> eigenvalue). >> > > > > The program also runs in parallel using mpiexec. I use complex >> builds, as in my main program the matrix can be complex. >> > > > > >> > > > > Now here are my questions: >> > > > > 1) For this particular problem type, could you please check if >> these are the best settings that one could use? I have tried different >> combinations of KSP/PC types e.g. GMRES, GAMG, etc, but BCGSL + BJACOBI >> seems to work the best in serial and parallel. >> > > > > >> > > > > 2) When I tested these settings in my main program, for some >> reason the JD solver was not converging. After further testing, I found the >> issue was related to the setting of "-eps_target 0". I have included >> "EPSSetTarget(eps,0.0);" in the program and I assumed this is equivalent to >> passing "-eps_target 0" from the command line, but that doesn't seem to be >> the case. For instance, if I run the attached program without "-eps_target >> 0" in the command line then it doesn't converge. >> > > > > $ ./acoustic_box_test.o -nx 1600 -ny 400 -usejd 1 -deflate 1 >> -eps_target 0 >> > > > > the above finishes in about 10 minutes >> > > > > $ ./acoustic_box_test.o -nx 1600 -ny 400 -usejd 1 -deflate 1 >> > > > > the above doesn't converge even though "EPSSetTarget(eps,0.0);" >> is included in the code >> > > > > >> > > > > This only seems to affect the JD solver, not the Krylov >> shift-and-invert (-sinvert 1) option. So is there any difference between >> passing "-eps_target 0" from the command line vs using >> "EPSSetTarget(eps,0.0);" in the code? I cannot pass any command line >> arguments in my actual program, so need to set everything internally. >> > > > > >> > > > > 3) Also, another minor related issue. While using the inexact >> shift-and-invert option, I was running into the following error: >> > > > > >> > > > > "" >> > > > > Missing or incorrect user input >> > > > > Shift-and-invert requires a target 'which' (see >> EPSSetWhichEigenpairs), for instance -st_type sinvert -eps_target 0 >> -eps_target_magnitude >> > > > > "" >> > > > > >> > > > > I already have the below two lines in the code: >> > > > > EPSSetWhichEigenpairs(eps,EPS_SMALLEST_MAGNITUDE); >> > > > > EPSSetTarget(eps,0.0); >> > > > > >> > > > > so shouldn't these be enough? If I comment out the first line >> "EPSSetWhichEigenpairs", then the code works fine. >> > > > >> > > > You should either do >> > > > >> > > > EPSSetWhichEigenpairs(eps,EPS_SMALLEST_MAGNITUDE); >> > > > >> > > > without shift-and-invert or >> > > > >> > > > EPSSetWhichEigenpairs(eps,EPS_TARGET_MAGNITUDE); >> > > > EPSSetTarget(eps,0.0); >> > > > >> > > > with shift-and-invert. The latter can also be used without >> shift-and-invert (e.g. in JD). >> > > > >> > > > I have to check, but a possible explanation why in your comment >> above (2) the command-line option -eps_target 0 works differently is that >> it also sets -eps_target_magnitude if omitted, so to be equivalent in >> source code you have to call both >> > > > EPSSetWhichEigenpairs(eps,EPS_TARGET_MAGNITUDE); >> > > > EPSSetTarget(eps,0.0); >> > > > >> > > > Jose >> > > > >> > > > > I have some more questions regarding setting the preconditioner >> for a quadratic eigenvalue problem, which I will ask in a follow-up email. >> > > > > >> > > > > Thanks for your help! >> > > > > >> > > > > -Varun >> > > > > >> > > > > >> > > > > On Thu, Jul 1, 2021 at 5:01 AM Varun Hiremath < >> varunhirem...@gmail.com> wrote: >> > > > > Thank you very much for these suggestions! We are currently using >> version 3.12, so I'll try to update to the latest version and try your >> suggestions. Let me get back to you, thanks! >> > > > > >> > > > > On Thu, Jul 1, 2021, 4:45 AM Jose E. Roman <jro...@dsic.upv.es> >> wrote: >> > > > > Then I would try Davidson methods https://doi.org/10.1145/2543696 >> > > > > You can also try Krylov-Schur with "inexact" shift-and-invert, >> for instance, with preconditioned BiCGStab or GMRES, see section 3.4.1 of >> the users manual. >> > > > > >> > > > > In both cases, you have to pass matrix A in the call to >> EPSSetOperators() and the preconditioner matrix via >> STSetPreconditionerMat() - note this function was introduced in version >> 3.15. >> > > > > >> > > > > Jose >> > > > > >> > > > > >> > > > > >> > > > > > El 1 jul 2021, a las 13:36, Varun Hiremath < >> varunhirem...@gmail.com> escribió: >> > > > > > >> > > > > > Thanks. I actually do have a 1st order approximation of matrix >> A, that I can explicitly compute and also invert. Can I use that matrix as >> preconditioner to speed things up? Is there some example that explains how >> to setup and call SLEPc for this scenario? >> > > > > > >> > > > > > On Thu, Jul 1, 2021, 4:29 AM Jose E. Roman <jro...@dsic.upv.es> >> wrote: >> > > > > > For smallest real parts one could adapt ex34.c, but it is going >> to be costly >> https://slepc.upv.es/documentation/current/src/eps/tutorials/ex36.c.html >> > > > > > Also, if eigenvalues are clustered around the origin, >> convergence may still be very slow. >> > > > > > >> > > > > > It is a tough problem, unless you are able to compute a good >> preconditioner of A (no need to compute the exact inverse). >> > > > > > >> > > > > > Jose >> > > > > > >> > > > > > >> > > > > > > El 1 jul 2021, a las 13:23, Varun Hiremath < >> varunhirem...@gmail.com> escribió: >> > > > > > > >> > > > > > > I'm solving for the smallest eigenvalues in magnitude. Though >> is it cheaper to solve smallest in real part, as that might also work in my >> case? Thanks for your help. >> > > > > > > >> > > > > > > On Thu, Jul 1, 2021, 4:08 AM Jose E. Roman < >> jro...@dsic.upv.es> wrote: >> > > > > > > Smallest eigenvalue in magnitude or real part? >> > > > > > > >> > > > > > > >> > > > > > > > El 1 jul 2021, a las 11:58, Varun Hiremath < >> varunhirem...@gmail.com> escribió: >> > > > > > > > >> > > > > > > > Sorry, no both A and B are general sparse matrices >> (non-hermitian). So is there anything else I could try? >> > > > > > > > >> > > > > > > > On Thu, Jul 1, 2021 at 2:43 AM Jose E. Roman < >> jro...@dsic.upv.es> wrote: >> > > > > > > > Is the problem symmetric (GHEP)? In that case, you can try >> LOBPCG on the pair (A,B). But this will likely be slow as well, unless you >> can provide a good preconditioner. >> > > > > > > > >> > > > > > > > Jose >> > > > > > > > >> > > > > > > > >> > > > > > > > > El 1 jul 2021, a las 11:37, Varun Hiremath < >> varunhirem...@gmail.com> escribió: >> > > > > > > > > >> > > > > > > > > Hi All, >> > > > > > > > > >> > > > > > > > > I am trying to compute the smallest eigenvalues of a >> generalized system A*x= lambda*B*x. I don't explicitly know the matrix A >> (so I am using a shell matrix with a custom matmult function) however, the >> matrix B is explicitly known so I compute inv(B)*A within the shell matrix >> and solve inv(B)*A*x = lambda*x. >> > > > > > > > > >> > > > > > > > > To compute the smallest eigenvalues it is recommended to >> solve the inverted system, but since matrix A is not explicitly known I >> can't invert the system. Moreover, the size of the system can be really >> big, and with the default Krylov solver, it is extremely slow. So is there >> a better way for me to compute the smallest eigenvalues of this system? >> > > > > > > > > >> > > > > > > > > Thanks, >> > > > > > > > > Varun >> > > > > > > > >> > > > > > > >> > > > > > >> > > > > >> > > > > <acoustic_box_test.cpp> >> > > > >> > > > <acoustic_matrix_test.cpp><MatA.gz> >> > > >> > >> >>