> El 28 sept 2021, a las 7:50, Varun Hiremath <[email protected]>
> 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 <[email protected]>
> 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 <[email protected]> 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 <[email protected]>
> > 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 <[email protected]> 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 <[email protected]>
> > > 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 <[email protected]> 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 <[email protected]>
> > > > 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 <[email protected]>
> > > > 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 <[email protected]> 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 <[email protected]>
> > > > > 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 <[email protected]> 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
> > > > > > <[email protected]> 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 <[email protected]>
> > > > > > wrote:
> > > > > > Smallest eigenvalue in magnitude or real part?
> > > > > >
> > > > > >
> > > > > > > El 1 jul 2021, a las 11:58, Varun Hiremath
> > > > > > > <[email protected]> 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 <[email protected]>
> > > > > > > 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
> > > > > > > > <[email protected]> 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>
> >
>
