Any precision control for solving Generalized Eigenvalues in FEM

[xpgo]

Dear yixuan,

Thank you very much for developing the package as a c++ alternative of ARPACK. It helps me a lot. Currently, I am using SPECTRA to solve the frequency problem in Finite Element Method.
Say, K u = v M u, where K is the stiffness matrix, M is the mass matrix, v and u are eigen pairs.

In order to find some smallest eigen values for the system, I created an OP for the shift-invet mode to solve the equivalent problem: (K - s B)^-1 B u = q u, let C = K - s B, by using SparseLU:

SparseLU<SparseMatrix, Eigen::COLAMDOrdering<int>> m_decomp;
...
in Construct function:
m_decomp.analyzePattern(matC.template selfadjointView<Uplo>());
m_decomp.factorize(matC.template selfadjointView<Uplo>());
...
    void perform_op(const Scalar* x_in, Scalar* y_out)
    {
        MapConstVec x(x_in, m_n);
        MapVec      y(y_out, m_n);
        m_cache = m_matB.template selfadjointView<Uplo>() * x;
        y = m_decomp.solve(m_cache);
    }

It works and very fast. However, the results are always different from that obtained from ARPACK computaiton, for example:

 No.         ARPACK             SPECTRA
  1   0.6770787E+10  0.68465780E+10
  2   0.1473508E+11  0.14354579E+11
  3   0.2330940E+12  0.20000338E+12
  4   0.2985047E+12  0.28088213E+12
  5   0.4432748E+12  0.47451875E+12
  6   0.1048882E+13  0.15452426E+13
  7   0.1542167E+13  0.20888521E+13
  8   0.2590512E+13  0.26107036E+13
  9   0.2692186E+13  0.41441039E+13
 10   0.4887708E+13  0.54881724E+13

We can see that, the first several values are acceptable, but it get worse as the frequecies increase, also the eigen vectors become worse.

I tried to improve the results by changing the tolenrance parameter in the compute() function. However, it hardly affect the results. Did I missed something? or is there any way I can improve the results?

I appreciate if you could give some hint. Thank you.

[xpgo]

I found the answer finally. For whom encounters the same problem: since C^-1 is not symmetric any more, we can not use symmetric related solver, we need to use GenEigsSolver.

[yixuan]

Good catch. Thanks for sharing these experiences!