SymGEigsShiftSolver problem

[wo80]

I started testing your latest code from the 1.y.z branch and found a problem. I'm solving a generalized problem using shift-invert and buckling mode with shift 1.0. While ARPACK gives the correct result for both modes

Testing ARPACK++ class ARluSymGenEig
Real symmetric generalized eigenvalue problem: A*x - lambda*B*x

Dimension of the system            : 100
Number of 'requested' eigenvalues  : 4
Number of 'converged' eigenvalues  : 4
Number of Arnoldi vectors generated: 9
Number of iterations taken         : 2

Eigenvalues:
  lambda[1]: 9,870400174642766
  lambda[2]: 39,49115121244263
  lambda[3]: 88,89091388108714
  lambda[4]: 158,11748682937144

||A*x(0) - lambda(0)*B*x(0)||: 7,256312682381566E-14
||A*x(1) - lambda(1)*B*x(1)||: 2,0601765361521015E-14
||A*x(2) - lambda(2)*B*x(2)||: 2,398614911913575E-09
||A*x(3) - lambda(3)*B*x(3)||: 2,15870050563603E-08

Spectra fails using the shift-invert mode, while the buckling mode succeeds:

Eigenvalues:
  lambda[1]: 1,993675588829463
  lambda[2]: 1,98875025628225
  lambda[3]: 1,9746778716421958
  lambda[4]: 1,8986869850963182

||A*x(0) - lambda(0)*B*x(0)||: 0,6188761805550974
||A*x(1) - lambda(1)*B*x(1)||: 0,4608353958197358
||A*x(2) - lambda(2)*B*x(2)||: 0,3007288554098257
||A*x(3) - lambda(3)*B*x(3)||: 0,1329643481854819

Here's the relevant code:

using OpType = SymShiftInvert<double, Eigen::Sparse, Eigen::Sparse>;
using BOpType = SparseSymMatProd<double>;

OpType op(A, B);
BOpType Bop(A);

SymGEigsShiftSolver<double, OpType, BOpType, GEigsMode::ShiftInvert> eigs(op, Bop, 4, 12, 1.0); // FAILS
//SymGEigsShiftSolver<double, OpType, BOpType, GEigsMode::Buckling> eigs(op, Bop, 4, 12, 1.0); // OK
eigs.init();
eigs.compute(SortRule::LargestMagn, 100, 1e-5); // SmallestMagn also fails

The test matrices are stored in MatrixMarket format: A.txt B.txt

[wo80]

To clarify: the main problem I see here is that Spectra returns after 33 iterations with compinfo SUCCESSFUL.

[yixuan]

I guess this is an easy fix: The roles of A and B are different in the two solvers.

ShiftInvert:

// A * x = lambda * B * x, A is symmetric, B is positive definite
// A is dense, B is sparse
using OpType = SymShiftInvert<double, Eigen::Dense, Eigen::Sparse>;
using BOpType = SparseSymMatProd<double>;
OpType op(A, B);
BOpType Bop(B);

Buckling:

// K * x = lambda * KG * x, K is positive definite, KG is symmetric
// K is sparse, KG is dense
using OpType = SymShiftInvert<double, Eigen::Sparse, Eigen::Dense>;
using BOpType = SparseSymMatProd<double>;
OpType op(K, KG);
BOpType Bop(K);

[wo80]

In this case, both matrices are spd, so this shouldn't be the problem.

[wo80]

Btw, I'm comparing results with the tests from ARPACK++:

Matrix A from examples/matrices/sym/lsmatrxc.h Matrix B from examples/matrices/sym/lsmatrxd.h

Shift-invert mode test from examples/superlu/sym/lsymgshf.cc Buckling mode test from examples/superlu/sym/lsymgbkl.cc

[yixuan]

I mean, BOpType is B in the shift-and-invert mode, and is K in the buckling mode.

[wo80]

Ok, I see, this does make a difference. But now the shift-invert mode does not converge at all.

I'm (still) working on a C wrapper for Spectra. Maybe you can have a short look at the relevant code: https://github.com/wo80/vs-spectra/blob/master/src/Visual%20Studio/shared/libspectra/spectra_di_sg.cpp#L80

Buckling and Caley mode work fine, but the regular shift-invert mode just doesn't ...

[wo80]

I had selection set to SortRule::SmallestMagn. Now it works as expected. Thanks!