[petsc-users] [SLEPc] generalized symmetric eigensystem

[Patrick Alken]

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 <pierre at joliv.et> escribió:
>>
>>
>>
>>> On 5 Aug 2022, at 5:04 AM, Patrick Alken <patrick.alken at geomag.info> 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>