Thank you both, I am getting the correct results now with these changes.
On 8/5/22 00:32, Jose E. Roman wrote:
You are setting the matrix values in the lower triangular part only. You are
assuming that SLEPc uses the same storage as LAPACK, but this is not true.
EPSSolve() does not work correctly because you are specifying a symmetric
problem EPS_GHEP but your matrices are not symmetric. You should set the upper
triangular part as well, or use the special storage SBAIJ. Anyway, your
matrices are dense, so it is better that you use a DENSE storage format. Note
that SLEPc is more appropriate for working with sparse matrices, not dense.
With both upper and lower triangular parts:
$ ./mwe -eps_nev 10 -eps_target 0 -st_type sinvert -mat_type dense
Number of requested eigenvalues: 10
Linear eigensolve converged (14 eigenpairs) due to CONVERGED_TOL; iterations 6
---------------------- --------------------
k ||Ax-kBx||/||kx||
---------------------- --------------------
1.309567 5.65098e-14
1.648674 1.67452e-09
1.709879 5.74072e-10
1.709879 7.52919e-09
2.142210 1.1093e-10
2.142210 3.24847e-09
2.145726 1.68857e-10
2.216526 9.25784e-11
2.216526 2.22353e-09
2.695837 2.30037e-09
2.712654 6.93094e-10
2.720031 1.63047e-09
2.734190 8.88298e-12
2.805326 1.31304e-09
---------------------- --------------------
As Pierre said, you want to avoid -eps_smallest_magnitude, see
https://slepc.upv.es/documentation/faq.htm#faq12
Jose
El 5 ago 2022, a las 7:47, Pierre Jolivet <[email protected]> escribió:
On 5 Aug 2022, at 5:04 AM, Patrick Alken <[email protected]> wrote:
I am trying to solve a generalized symmetric eigensystem problem using SLEPc,
and I am getting much different eigenvalues compared with the dense LAPACK
solver DSYGVD. I have attached a minimal working example (MWE). The code relies
on the external libraries GSL, Lapack and Lapacke.
The matrices (A,B) come from discretizing the Laplacian in spherical
coordinates using a radial basis function method. Here are the eigenvalues
computed by the dense LAPACK solver:
=== LAPACK eigenvalues ===
eval(0) = 1.309567289750e+00
eval(1) = 1.648674331144e+00
eval(2) = 1.709879393690e+00
eval(3) = 1.709879393690e+00
eval(4) = 2.142210251638e+00
eval(5) = 2.142210251638e+00
eval(6) = 2.145726457896e+00
eval(7) = 2.216526311324e+00
eval(8) = 2.216526311324e+00
eval(9) = 2.695837030236e+00
Here are the eigenvalues computed by SLEPc, using the command:
$ ./mwe -eps_nev 10 -eps_smallest_magnitude -eps_tol 1e-8
Number of requested eigenvalues: 10
Linear eigensolve converged (10 eigenpairs) due to CONVERGED_TOL; iterations 196
---------------------- --------------------
k ||Ax-kBx||/||kx||
---------------------- --------------------
0.094520 906.031
0.177268 363.291
1.324501 61.9577
1.786472 35.7082
2.187857 32.8364
-2.886905 27.7314
2.899206 24.2224
4.195222 20.3007
4.787192 12.8346
7.221589 9.58513
---------------------- --------------------
As we can see, the results are quite different than LAPACK and the error norms
in the right column are quite large. Using a stricter tolerance (eps_tol) does
not seem to help. Can anyone suggest how to diagnose and fix the problem?
You probably want to use a spectral transformation to find the smallest
magnitude eigenpairs.
I cannot try your code because of its dependencies, but I bet that if you
replace -eps_smallest_magnitude by -eps_target 0.0 -st_type sinvert, you’ll get
a match with the LAPACK eigenvalues.
Thanks,
Pierre
<laplace3d.c><le_laplace3d.h><Makefile.txt><mwe.c>
