SymGEigsSolver settings for large sparse problem

[armingeiser]

Hello,

I am trying to use the library to solve for a structural eigenfrequency probem. This gives me a sparse (60kx60k) generalized eigenvalue problem.

I am using the Spectra::SparseCholesky and Spectra::SparseSymMatProd I am looking for the 5 lowest eigenvalues, however the solution with the Spectra lib fails.

I transfered the matrices A and B to python and solved the eigenvalue problem using the scipy "eigs" function, that, as far as i know uses the ARPACK library.

from scipy.sparse.linalg import eigs
evals_small, evecs_small = eigs(A, 5, B, sigma=0, which='LM')

There i was able to obtain the correct eigenvalues.

Is there a way to use the SymGEigsSolver with shift mode?

Thanks!

[yixuan]

I haven't verified it yet, but here is what you can try: SymGEigsSolver has two operators to define, one for A and one for B. The built-in class SparseCholesky takes care of the B operator, so what's left is a class that computes inv(A) * x. You can borrow some code from the SparseSymShiftSolve class, but in your case the sigma parameter is fixed to zero.

[armingeiser]

Thank you for the hint, I will try it!

[dilevin]

Hi,

I'm also interested in finding the smallest Eigenvalues of a large sparse matrix. If this worked for the previous poster, would you mind posting the code ?

Thanks!

[yixuan]

@dilevin Are you talking about the eigenvalue of a single matrix, or the generalized eigenvalue involving two matrices?

[dilevin]

Thanks for the quick response! I’d like to solve the generalized eigenvalue problem, but find the smallest eigenvalue efficiently.

[WindingWinter]

@yixuan , I second @dilevin . It would be supremely useful if you can do a class for generalized eigenvalue problem with shift inverted mode.

[armingeiser]

I did not find the time yet to try the suggestions by @yixuan, sorry for the late response.

[WindingWinter]

@armingeiser , how do you solve this problem of finding eigenvalues/modes to generalized eigenvalue problem with shift inverted mode?

[armingeiser]

@nsoonhui as i wrote above, unfortunately i did not solve/try it yet.

[dilevin]

@armingeiser @nsoonhui

For my particular problem, modal analysis for finite element systems, I found a pretty good solution which is to solve a mass shifted generalized eigenvalue problem:

   Eigen::SparseMatrix<DataType> K = -A + shift*B;
    Eigen::SparseMatrix<DataType> M = B;

    Spectra::SparseSymMatProd<DataType> Aop(M);
    Spectra::SparseCholesky<DataType> Bop(K);

    // Construct eigen solver object, requesting the smallest three eigenvalues
    Spectra::SymGEigsSolver<DataType, Spectra::LARGEST_MAGN,      Spectra::SparseSymMatProd<DataType>, Spectra::SparseCholesky<DataType>, Spectra::GEIGS_CHOLESKY > eigs(&Aop, &Bop, numVecs, 5*numVecs);

Then correct the eigenvalues using: -1.0/(eigs.eigenvalues()[ii]) + shift);

Here A is my input stiffness matrix (symmetric -ve semidefinite), B is my mass matrix and shift is a small number used to stabilize the stiffness matrix.

You have to renormalize your eigenvectors (wrt to the mass matrix) after retrieving them from spectra, but this seems to work and is much faster than using SMALLEST_MAGN.

[yixuan]

For people who are still looking for the shift-and-invert mode of SymGEigsSolver, I have just implemented a new SymGEigsShiftSolver class in the 1.y.z branch. Example code is available in the testing program.

Note that this new branch contains several API-breaking changes that will be included in the next major release, so some (minor) code migration is needed.

[giacomo-b]

For anyone having similar problems: I followed the approach suggested by @dilevin, unfortunately without success.

I am also dealing with modal analysis, my matrices are about 20,000 x 20,000.

The following worked for me:

int n_eigs = 3, convergence_ratio = 5;

double shift = 0;

// K is the stiffness matrix, M the mass matrix
Eigen::SparseMatrix<double> K_new = K + shift * M;
Eigen::SparseMatrix<double> M_new = M;

Spectra::SparseSymMatProd<double> Aop(M_new);
Spectra::SparseCholesky<double> Bop(K_new);

Spectra::SymGEigsSolver<double, Spectra::LARGEST_MAGN, Spectra::SparseSymMatProd<double>, Spectra::SparseCholesky<double>, Spectra::GEIGS_CHOLESKY> eigs(&Aop, &Bop, n_eigs, convergence_ratio * n_eigs);

eigs.init();
int n_conv = eigs.compute();

Eigen::VectorXd evalues;
if (eigs.info() == SUCCESSFUL)
    evalues = 1.0 / eigs.eigenvalues().array() + shift;

std::cout << "Generalized eigenvalues found:\n" << evalues << std::endl;

Note that I left shift = 0 in case someone needed to shift the eigenvalues, but I could just write:

Spectra::SparseSymMatProd<double> Aop(M);
Spectra::SparseCholesky<double> Bop(K);

[armingeiser]

Thanks @yixuan!

I have been successful using the SymGEigsShiftSolver exactly following the mentioned test program. FFR: i used it to solve a modal analysis of a solid structure with ~100k dofs.

From my point of view this issue can be closed.

[yixuan]

Spectra v1.0.0 has been released, and SymGEigsShiftSolver is fully supported now. Let me close this issue.