CODEWITHSUNDEEP
c++visual-studio-2015lapackarmadillo
Aaron
Aaron

Reputation: 57

Armadillo eigs_sym: decomposition failed

I am writing C++ code within Visual Studio and using Armadillo 7.900.1

I am having no luck getting eigs_sym function working under Armadillo (I am using the lapack and blas versions that came with Armadillo). The typical error message is as below:

error: lapack::stedc(): failed to compute all eigenvalues
warning: eigs_sym(): decomposition failed

The code that is producing this output is:

sp_mat T(locations, values);
arma::eigs_sym(eigval, eigvec, T, num_eigs_wanted, "sm", tol);

where 

T is a 480,000x480,000 sparse matrix.

I can get the code working if T is small (ie. 2000x2000) and tol is high (tol = 5). But once T is about 20000x20000, then no matter what values of tol or num_eigs_wanted are used then the above error occurs. (I obviously change "values" and locations" to alter the size of T).

The matrix T is symmetric, real, positive definite.

The exact same matrix has no problems when I call the eigs function within Matlab.

Any ideas? To me it seems that eigs_sym in Armadillo is broken...any alternatives that people have used?

Cheers

Upvotes: 3

Views: 1279

Answers (2)

Aaron
Aaron

Reputation: 57

Thanks @Henri for all your help! Spectra (https://spectralib.org) is the way to go!

The eigenvalue/vestor decomposition using SymEigsSolver on a sparse matrix of size 480,000 x 480,000 takes me 3 mins 50 sec (i7/16GB RAM).

A sparse matrix of size 120,000x120,000 takes 20 seconds.

No stability problems and the programmer has direct access to the Arnoldi Algorithm, Ritz values and tolerance threshold.

Easy to use as well since Spectra is purely header based. Just #include the header files - simple. You'll need Eigen as well which is also available at the above address (and it's also header based).

Upvotes: 1

Henri Menke
Henri Menke

Reputation: 10939

Not an answer but way too long for a comment. I can reproduce with the code below. The matrix is just a diagonal matrix with ascending values on the diagonal starting at 1. The eigenvalues are obviously 1,2,3,... and the matrix is positive definite and symmetric.

#include <iostream>
#include <armadillo>

int main()
{
  constexpr int N = 20000;
  arma::umat locations(2,N);
  arma::vec values(N);

  for (int i = 0; i < N; ++i)
  {
    locations(0,i) = i;
    locations(1,i) = i;
    values(i) = i+1;
  }

  arma::sp_mat T(locations, values);

  int num_eigs_wanted = 1;
  arma::vec eigval(num_eigs_wanted);
  arma::mat eigvec(N,num_eigs_wanted);
  double tol = 1e-6;
  arma::eigs_sym(eigval, eigvec, T, num_eigs_wanted, "sm", tol);

  std::cout << eigval << '\n';
}

Compiler flags are

clang++-5.0 -std=c++11 test.cpp -larmadillo

I'm using Armadillo with MKL backend and I receive:

warning: eigs_sym(): decomposition failed
[matrix size: 0x1]

It seems to be a problem of Armadillo, because when I use the unsupported Arpack support of Eigen it works just fine.

#include <iostream>
#include <Eigen/Sparse>
#include <unsupported/Eigen/ArpackSupport>

int main()
{
  constexpr int N = 20000;

  std::vector< Eigen::Triplet<double> > triplets;
  triplets.reserve(N);
  for (int i = 0; i < N; ++i)
    triplets.push_back( {i,i,i+1} );

  Eigen::SparseMatrix<double> T(N,N);
  T.setFromTriplets(triplets.begin(), triplets.end());

  int num_eigs_wanted = 1;
  double tol = 1e-6;
  Eigen::ArpackGeneralizedSelfAdjointEigenSolver< Eigen::SparseMatrix<double> > eigs_sym;
  eigs_sym.compute(T, num_eigs_wanted, "SM", Eigen::ComputeEigenvectors, tol);
  std::cout << eigs_sym.eigenvalues() << '\n';
}

You might also want to take a look at https://github.com/yixuan/spectra/

Upvotes: 2

Related Questions