A header-only C++ library for large scale eigenvalue problems

Overview

Build Status Basic CI codecov

NOTE: Spectra 1.0.0 is released, with a lot of API-breaking changes. Please see the migration guide for a smooth transition to the new version.

NOTE: If you are interested in the future development of Spectra, please join this thread to share your comments and suggestions.

Spectra stands for Sparse Eigenvalue Computation Toolkit as a Redesigned ARPACK. It is a C++ library for large scale eigenvalue problems, built on top of Eigen, an open source linear algebra library.

Spectra is implemented as a header-only C++ library, whose only dependency, Eigen, is also header-only. Hence Spectra can be easily embedded in C++ projects that require calculating eigenvalues of large matrices.

Relation to ARPACK

ARPACK is a software package written in FORTRAN for solving large scale eigenvalue problems. The development of Spectra is much inspired by ARPACK, and as the full name indicates, Spectra is a redesign of the ARPACK library using the C++ language.

In fact, Spectra is based on the algorithm described in the ARPACK Users' Guide, the implicitly restarted Arnoldi/Lanczos method. However, it does not use the ARPACK code, and it is NOT a clone of ARPACK for C++. In short, Spectra implements the major algorithms in ARPACK, but Spectra provides a completely different interface, and it does not depend on ARPACK.

Common Usage

Spectra is designed to calculate a specified number (k) of eigenvalues of a large square matrix (A). Usually k is much smaller than the size of the matrix (n), so that only a few eigenvalues and eigenvectors are computed, which in general is more efficient than calculating the whole spectral decomposition. Users can choose eigenvalue selection rules to pick the eigenvalues of interest, such as the largest k eigenvalues, or eigenvalues with largest real parts, etc.

To use the eigen solvers in this library, the user does not need to directly provide the whole matrix, but instead, the algorithm only requires certain operations defined on A. In the basic setting, it is simply the matrix-vector multiplication. Therefore, if the matrix-vector product A * x can be computed efficiently, which is the case when A is sparse, Spectra will be very powerful for large scale eigenvalue problems.

There are two major steps to use the Spectra library:

  1. Define a class that implements a certain matrix operation, for example the matrix-vector multiplication y = A * x or the shift-solve operation y = inv(A - σ * I) * x. Spectra has defined a number of helper classes to quickly create such operations from a matrix object. See the documentation of DenseGenMatProd, DenseSymShiftSolve, etc.
  2. Create an object of one of the eigen solver classes, for example SymEigsSolver for symmetric matrices, and GenEigsSolver for general matrices. Member functions of this object can then be called to conduct the computation and retrieve the eigenvalues and/or eigenvectors.

Below is a list of the available eigen solvers in Spectra:

Examples

Below is an example that demonstrates the use of the eigen solver for symmetric matrices.

#include <Eigen/Core>
#include <Spectra/SymEigsSolver.h>
// <Spectra/MatOp/DenseSymMatProd.h> is implicitly included
#include <iostream>

using namespace Spectra;

int main()
{
    // We are going to calculate the eigenvalues of M
    Eigen::MatrixXd A = Eigen::MatrixXd::Random(10, 10);
    Eigen::MatrixXd M = A + A.transpose();

    // Construct matrix operation object using the wrapper class DenseSymMatProd
    DenseSymMatProd<double> op(M);

    // Construct eigen solver object, requesting the largest three eigenvalues
    SymEigsSolver<DenseSymMatProd<double>> eigs(op, 3, 6);

    // Initialize and compute
    eigs.init();
    int nconv = eigs.compute(SortRule::LargestAlge);

    // Retrieve results
    Eigen::VectorXd evalues;
    if(eigs.info() == CompInfo::Successful)
        evalues = eigs.eigenvalues();

    std::cout << "Eigenvalues found:\n" << evalues << std::endl;

    return 0;
}

Sparse matrix is supported via classes such as SparseGenMatProd and SparseSymMatProd.

#include <Eigen/Core>
#include <Eigen/SparseCore>
#include <Spectra/GenEigsSolver.h>
#include <Spectra/MatOp/SparseGenMatProd.h>
#include <iostream>

using namespace Spectra;

int main()
{
    // A band matrix with 1 on the main diagonal, 2 on the below-main subdiagonal,
    // and 3 on the above-main subdiagonal
    const int n = 10;
    Eigen::SparseMatrix<double> M(n, n);
    M.reserve(Eigen::VectorXi::Constant(n, 3));
    for(int i = 0; i < n; i++)
    {
        M.insert(i, i) = 1.0;
        if(i > 0)
            M.insert(i - 1, i) = 3.0;
        if(i < n - 1)
            M.insert(i + 1, i) = 2.0;
    }

    // Construct matrix operation object using the wrapper class SparseGenMatProd
    SparseGenMatProd<double> op(M);

    // Construct eigen solver object, requesting the largest three eigenvalues
    GenEigsSolver<SparseGenMatProd<double>> eigs(op, 3, 6);

    // Initialize and compute
    eigs.init();
    int nconv = eigs.compute(SortRule::LargestMagn);

    // Retrieve results
    Eigen::VectorXcd evalues;
    if(eigs.info() == CompInfo::Successful)
        evalues = eigs.eigenvalues();

    std::cout << "Eigenvalues found:\n" << evalues << std::endl;

    return 0;
}

And here is an example for user-supplied matrix operation class.

#include <Eigen/Core>
#include <Spectra/SymEigsSolver.h>
#include <iostream>

using namespace Spectra;

// M = diag(1, 2, ..., 10)
class MyDiagonalTen
{
public:
    using Scalar = double;  // A typedef named "Scalar" is required
    int rows() { return 10; }
    int cols() { return 10; }
    // y_out = M * x_in
    void perform_op(const double *x_in, double *y_out) const
    {
        for(int i = 0; i < rows(); i++)
        {
            y_out[i] = x_in[i] * (i + 1);
        }
    }
};

int main()
{
    MyDiagonalTen op;
    SymEigsSolver<MyDiagonalTen> eigs(op, 3, 6);
    eigs.init();
    eigs.compute(SortRule::LargestAlge);
    if(eigs.info() == CompInfo::Successful)
    {
        Eigen::VectorXd evalues = eigs.eigenvalues();
        std::cout << "Eigenvalues found:\n" << evalues << std::endl;
    }

    return 0;
}

Shift-and-invert Mode

When it is needed to find eigenvalues that are closest to a number σ, for example to find the smallest eigenvalues of a positive definite matrix (in which case σ = 0), it is advised to use the shift-and-invert mode of eigen solvers.

In the shift-and-invert mode, selection rules are applied to 1/(λ - σ) rather than λ, where λ are eigenvalues of A. To use this mode, users need to define the shift-solve matrix operation. See the documentation of SymEigsShiftSolver for details.

Documentation

The API reference page contains the documentation of Spectra generated by Doxygen, including all the background knowledge, example code and class APIs.

More information can be found in the project page https://spectralib.org.

Installation

An optional CMake installation is supported, if you have CMake with at least v3.10 installed. You can install the headers using the following commands:

mkdir build && cd build
cmake .. -DCMAKE_INSTALL_PREFIX='intended installation directory' -DCMAKE_PREFIX_PATH='path where the installation of Eigen3 can be found' -DBUILD_TESTS=TRUE
make all && make tests && make install

By installing Spectra in this way, you also create a CMake target Spectra::Spectra that can be used in subsequent build procedures for other programs.

License

Spectra is an open source project licensed under MPL2, the same license used by Eigen.

Comments
  • Build issues because of Eigen and warnings with Clang-format

    Build issues because of Eigen and warnings with Clang-format

    Hi,

    I'm trying to add Spectra to a libigl project I have. My CMakeLists.txt is this one:

    cmake_minimum_required(VERSION 3.1 FATAL_ERROR)
    project(TestProject)
    set(CMAKE_CXX_STANDARD 14)
    
    set(SRC_DIR ${CMAKE_CURRENT_SOURCE_DIR}/src)
    set(INCLUDE_DIR ${CMAKE_CURRENT_SOURCE_DIR}/include)
    
    add_subdirectory(${LIBIGL_DIR} libigl)
    add_subdirectory(${SPECTRA_DIR} spectra)
    
    # Add your project files
    add_executable(test ${SRC_DIR}/main.cpp )
    target_link_libraries(test igl::opengl_glfw)
    target_include_directories(test PUBLIC ${INCLUDE_DIR})
    

    However when I run this I get:

    /usr/local/bin/cmake -DCMAKE_BUILD_TYPE=Debug -DLIBIGL_DIR=../third_party/libigl -DSPECTRA_DIR=../third_party/spectra/ -DCMAKE_MODULE_PATH=../third_party/spectra/cmake -G "CodeBlocks - Unix Makefiles" /home/TestProject
    -- Libigl top level project: OFF
    -- Creating target: igl::core (igl)
    -- Creating target: igl::comiso (igl_comiso)
    CMake Warning (dev) at /home/third_party/libigl/external/embree/common/cmake/test.cmake:31 (SET):
      implicitly converting 'INT' to 'STRING' type.
    Call Stack (most recent call first):
      /home/third_party/libigl/external/embree/CMakeLists.txt:101 (INCLUDE)
    This warning is for project developers.  Use -Wno-dev to suppress it.
    
    -- Creating target: igl::embree (igl_embree)
    -- Creating target: igl::opengl (igl_opengl)
    -- Creating target: igl::opengl_glfw (igl_opengl_glfw)
    -- Using X11 for window creation
    -- Creating target: igl::opengl_glfw_imgui (igl_opengl_glfw_imgui)
    -- Creating target: igl::png (igl_png)
    -- Creating target: igl::tetgen (igl_tetgen)
    -- Creating target: igl::triangle (igl_triangle)
    -- Creating target: igl::predicates (igl_predicates)
    -- Creating target: igl::xml (igl_xml)
    -- Found Eigen3 Version:  Path: /usr/lib/cmake/eigen3
    CMake Error: File /home/TestProject/cmake/spectra-config.cmake.in does not exist.
    CMake Error at /usr/local/share/cmake-3.14/Modules/CMakePackageConfigHelpers.cmake:330 (configure_file):
      configure_file Problem configuring file
    Call Stack (most recent call first):
      /home/third_party/spectra/CMakeLists.txt:69 (configure_package_config_file)
    
    
    -- Could NOT find CLANG_FORMAT: Found unsuitable version "3.8.0", but required is at least "7.0.1" (found /usr/bin/clang-format)
    -- The following OPTIONAL packages have been found:
    
     * Git
     * Freetype
     * PkgConfig
     * Fontconfig
    
    -- The following REQUIRED packages have been found:
    
     * BLAS
     * OpenGL
     * Threads
     * X11
     * Eigen3
       C++ vector data structures
    
    -- The following OPTIONAL packages have not been found:
    
     * CLANG_FORMAT (required version >= 7.0.1)
    
    -- Configuring incomplete, errors occurred!
    See also "/home/TestProject/cmake-build-debug/CMakeFiles/CMakeOutput.log".
    See also "/home/TestProject/cmake-build-debug/CMakeFiles/CMakeError.log".
    
    [Failed to reload]
    

    However I have clang-format-7 installed (I did it using apt get). Can you help?

    (I've also tried to follow: https://github.com/yixuan/spectra#installation, but I get the same error).

    opened by lukkio88 36
  • the sign of the eigenvectors returned by SymGEigsSolver is not correct

    the sign of the eigenvectors returned by SymGEigsSolver is not correct

    Hi Yixuan, Recently I am using spectra to solve a general sparse Ax=lamdaBx problem. I compared the results with Matlab. I found the eigenvectors had the same value but the sign is messed up. Totally half value is different. I suspect there is something wrong with the code.

    The code I used is from your example code. typedef SparseSymMatProd OpType; typedef SparseRegularInverse BOpType; SymGEigsSolver<double, Spectra::SMALLEST_MAGN, OpType, BOpType, GEIGS_REGULAR_INVERSE> eigs( &op, &Bop, nev, ncv ); eigs.init( ); int nconv = eigs.compute( n*3); // maxit = 100 to reduce running time for failed cases int niter = eigs.num_iterations( ); int nops = eigs.num_operations( ); Vector evals = eigs.eigenvalues( ); Matrix evecs = eigs.eigenvectors( ); could you help have a check? Thanks. data.zip

    opened by YuchuanQiao 16
  • SymGEigsSolver settings for large sparse problem

    SymGEigsSolver settings for large sparse problem

    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!

    opened by armingeiser 14
  • SymGEigsSolver: Strange results since the update from

    SymGEigsSolver: Strange results since the update from "update badge link" (30. may 20) to "merge pull request" (12. jan 21)

    Thanks to your migration guide the update was easy for me.

    I am (still) using SymGEigsSolver. But after the update, I get too many eigenmodes/-frequencies returned for one of my test cases. The test case has only a few degrees of freedom (ndof = 6) but by default I cause a certain number of eigenmodes (e.g. nef = 10) to be calculated. Previously this was not a problem and I received few (Due to the nature of my matrices 3) but fully reasonable eigenmodes. After the update, I now receive in total min(nef, ndof-1) (5) eigenmodes/-frequencies. Some (1) are now missing and some (3) of them are not reasonable.

    In all my other test cases ndof is (usually much) larger than nef and among them, there is only one case with a very very small numerical deviation of no significance. All others are exactly the same up to the very last digit.

    I will try to put together a test code for you. But at the moment I can't estimate if and until when I will succeed.

    Motivation: My test system is mechanically a simple cantilever made from a simple bar element with a point mass at the end. I then feed SymGEigsSolver the following problem: ([M] + (1/w^2) [C+sM]) x = 0 -> f = sqrt(w^2 - s)/2pi det([M]) >= 0 det([C+sM]) > 0

    This mechanical system has, if you want to see it that way, 3 infinitely large natural frequencies, which did not disturb so far. Previous, meaningful natural frequencies [Hz] 1: .13891786E+01 2: .13891786E+01 3: .54689064E+02

    Now, additional, non-sense eigenfrequencies [Hz] 1: .00000000E+00 2: .13891786E+01 3: .20267187E+02 4: .25771982E+02 5: .54689064E+02

    opened by Spammed 13
  • fail to find eigenvalue for some cases

    fail to find eigenvalue for some cases

    I have a matrix of size 30000x30000, very sparse, denote as A. We know A's rank should be 29998. I need to compute serveral singular values (ones close to 0) of A and corresponding right singular vectors. This equivalents to compute eigenvalues of AA=A'*A and corresponding right eigenvectors. My code is as following

    Eigen::SparseMatrix<double> AA = A.transpose()*A;
    Spectra::SparseSymMatProd<double> op(AA);
    Spectra::SymEigsSolver< double, Spectra::SMALLEST_ALGE, Spectra::SparseSymMatProd<double>> eigs(&op, kd, max(3 * kd, 20)); // kd = 2
    eigs.init();
    int nconv = eigs.compute();
    

    I got no eigenvalue converged and takes more then 10 seconds.

    We can use matlab to do this:

    [V,D] = eigs(A'*A,2,'SM')
    

    it takes about 0.3 second and successfully obtained 2 eigenvalues 4.2590e-18 and -5.7298e-17 Matlab makes a call to arpack internally to compute eigenvalues.

    The matrix is here: https://1drv.ms/u/s!AoRwrBYa5axzleIlJ2KzlxF1Fged8w, it stores in COO format.

    In another case with a matrix of 2298x2304, spectra did succeed to compute eigenvalues. So I wonder if size or matrix structures matters (two matrix have similar structure)

    opened by cfwen 13
  • Eigenvalues are different for direct EIGEN implementation vs Spectra based EIGEN implementation.

    Eigenvalues are different for direct EIGEN implementation vs Spectra based EIGEN implementation.

    Hi, Thanks for the fastest tool for EIGEN implementation, it's so useful for very large matrices. I am working with a FE project which is related to solving eigenvalues of a system.

    I am getting different eigen values with Spectra vs direct EIGEN/MATLAB implementation. For the same matrices check.txt as A and Mcheck.txt as MM, I solved using Spectra and Eigen as follows, The both matrices are symmetric and positive definite.

    //-----------------in -----Spectra ---------- MatrixXd A(n,n); //some file reading code here and memcpy into A and M_lin matrices here. DenseSymMatProd op(A); SparseMatrix MM(n,n); MM.reserve(Eigen::VectorXi::Constant(n, 3));

    SparseCholesky Bop(MM); int num_eigval_want = 1398; int converge_speed = 1400; SymGEigsSolver<double, WHICH_SM, DenseSymMatProd, SparseCholesky, GEIGS_CHOLESKY>geigs(&op, &Bop, num_eigval_want, converge_speed); geigs.init(); int nconv = geigs.compute(); VectorXd evalues; MatrixXd evecs; if(geigs.info() == SUCCESSFUL){ evalues = geigs.eigenvalues(); evecs = geigs.eigenvectors();} std::cout << "Generalized eigenvalues found:\n" << evalues << std::endl; //std::cout << "Generalized eigenvectors found:\n" << evecs.topRows(num_eigval_want) << std::endl;

    //-----------------in -----direct Eigen ---------- MatrixXd A(n,n); MatrixXd MM(n,n); //file copying code here into A and MM GeneralizedEigenSolver ges; ges.compute(A, MM); cout << "The (complex) generalzied eigenvalues are (alphas./beta): " << ges.eigenvalues().transpose() << endl;

    Can you please correct me, if I am doing something wrong here. The attached eig.txt (contains (freq)^2) is the values obtained by direct EIGEN/MATLAB implementation. Seig.txt (contains (freq)^2 in first column) is Spectra implementation, which is having different eigenvalues than direct EIGEN/MATLAB.

    Thanks in advance.

    check.txt eig.txt Mcheck.txt Seig.txt

    opened by Maharshi-kintada 12
  • rank deficient matrices

    rank deficient matrices

    Hi @yixuan , Can I ask you what is the expected behavior of SymEigsShiftSolver in presence of matrices with some zero eigenvalues? It seems that SparseSymShiftSolve relies on Eigen::SparseLU , which can have a problem on line 83 of SparseGenRealShiftDolver.h : y.noalias() = m_solver.solve(x); for non-invertible matrices...

    opened by jdbancal 11
  • performance problem with SymGEigsSolver

    performance problem with SymGEigsSolver

    Sir:

    I have a performance problem with SymGEigsSolver.

    I have a generalized eigenvalue problem, the system is sparse, and with size 960x960.

    I want find 5 smallest eigenvalue and corresponding eigenvector.

    I using spectra, I found it is hard to convergence, at I tune the parameter to let solver convergence, it takes around 14s on my computer.

    I using matlab to solve the same problem, it takes 0.07s to obtain smallest 5 eigenvalue and takes 0.16s to obtain all eigen value.

    I do not know why there exist a so large gap, I already enable the -O3 optimization. Maybe I misused the library.

    My code using spectra looks like

    #include <Eigen/Core>
    #include <Eigen/SparseCore>
    #include <Eigen/Eigenvalues>
    #include <Spectra/SymGEigsSolver.h>
    #include <Spectra/MatOp/SparseGenMatProd.h>
    #include <Spectra/MatOp/SparseCholesky.h>
    #include <iostream>
    #include <fstream>
    std::vector<Eigen::Triplet<double>> readTripletFromFile(const std::string &file_name, const int nonzero)
    {
        std::vector<Eigen::Triplet<double>> coefficients(nonzero);
        std::ifstream fin;
        fin.open(file_name);
        int row, col;
        double value;
        for (int i = 0; i < nonzero; i++)
        {
            fin >> row;
            fin >> col;
            fin >> value;
            coefficients[i] = Eigen::Triplet<double>(row, col, value);
        }
        fin.close();
        return coefficients;
    }
    
    int main()
    {
        const int n = 960;
        const int nonzero = 184320;
        std::vector<Eigen::Triplet<double>> coefficients_A = readTripletFromFile("data_A.txt", nonzero);
        std::vector<Eigen::Triplet<double>> coefficients_B = readTripletFromFile("data_B.txt", nonzero);
        Eigen::SparseMatrix<double> A(n, n);
        A.setFromTriplets(coefficients_A.begin(), coefficients_A.end());
        Eigen::SparseMatrix<double> B(n, n);
        B.setFromTriplets(coefficients_B.begin(), coefficients_B.end());
    
        Spectra::SparseSymMatProd<double> Aop(A);
        Spectra::SparseCholesky<double> Bop(B);
    
        const int num_solution = 5;
        const int convergence_parameter = 100;
    
        Spectra::SymGEigsSolver<double, Spectra::SMALLEST_ALGE,
                                Spectra::SparseSymMatProd<double>,
                                Spectra::SparseCholesky<double>,
                                Spectra::GEIGS_CHOLESKY>
            eigen_solver(&Aop, &Bop, num_solution, convergence_parameter);
    
        eigen_solver.init();
        int num_convergence = eigen_solver.compute();
    
        if (eigen_solver.info() == Spectra::SUCCESSFUL)
        {
            std::cout << "success." << std::endl;
            std::cout << eigen_solver.eigenvalues()
                      << std::endl;
        }
        else
        {
            std::cout << "failed." << std::endl;
        }
        return 0;
    }
    

    The full code can be found in https://github.com/qzxuhui/MyDiscussion/tree/master/spectra/1

    Could anyone help me? Thanks for your time.

    opened by qzxuhui 10
  • Issues with Sparse Matrices Loaded from File

    Issues with Sparse Matrices Loaded from File

    Hi, I've been trying to use the Spectra library for analyzing some sparse matrices that I have. My basic function for getting eigen pairs looks like this:

    pair<Eigen::MatrixXcd, Eigen::VectorXcd> getEigPairs(SMatrixXd &matrix, string &fname, int k){
    
         bool status = Eigen::loadMarket(matrix, fname);
         if (status)
         {
             printf("Matrix succesfully loaded\n");
             printf("Stats:\n");
             printf("Rows: %ld Cols: %ld, NNZ: %ld\n",  matrix.rows(), 
                matrix.cols(), matrix.nonZeros());
             // printSparseData(matrix);
         } else{
            printf("Couldn't load matrix\n");
            return {};
         } 
    
         Eigen::VectorXcd eigvalues; 
         Eigen::MatrixXcd eigvectors;
    
         Spectra::SparseGenMatProd<double> op(matrix); 
         // // last param is the convergence speed 
         cout << "Creating eig solver now...\n";
         Spectra::SymEigsSolver<double, Spectra::LARGEST_MAGN, Spectra::SparseGenMatProd<double> > eigs(&op, k, k*2);
         // Spectra::GenEigsSolver<double, Spectra::LARGEST_MAGN, Spectra::SparseGenMatProd<double> > eigs(&op, k, k*2);  
    
         cout << eigs.info();
    
    
         int nconv = eigs.compute(); 
         cout << nconv << std::endl;
         cout << eigs.eigenvalues() << "\n";
         if (eigs.info() == Spectra::SUCCESSFUL)
         { 
            eigvalues = eigs.eigenvalues();
            eigvectors = eigs.eigenvectors(k);
         } else{
            cout << "Error in computing eigenpairs";
         }
    
         return {eigvectors, eigvalues};
    }
    
    

    Here SMatrixXd is a typedef for a rowmajor column double matrix. I have tried both GenEigsSolver and SymEigsSolver. In either of them I get an error that says:

    image

    And for real symmetric matrices, I get the same error but with Lanczos number. Is there anything wrong I'm doing? Or is it something else?

    opened by aicaffeinelife 10
  • The core dumped can happen when using row major Eigen::SparseMatrix

    The core dumped can happen when using row major Eigen::SparseMatrix

    The core dumped can happen when using row major EigenSparseMatrix.

    I am not sure it is a my code bug or spectra's bug.

    The follow code can used to reproduce problem.

    #include <Eigen/Core>
    #include <Eigen/SparseCore>
    #include <Eigen/Eigenvalues>
    #include <Spectra/SymGEigsSolver.h>
    #include <Spectra/MatOp/SparseGenMatProd.h>
    #include <Spectra/MatOp/SparseCholesky.h>
    #include <iostream>
    int main()
    {
        // We are going to solve the generalized eigenvalue problem A * x = lambda * B * x
        const int n = 10;
    
        // Define the A matrix, a band matrix with 3 on the diagonal and 1 on the subdiagonals
        // Define the B matrix, a band matrix with 2 on the diagonal and 1 on the subdiagonals
        Eigen::SparseMatrix<double, Eigen::RowMajor> A(n, n);
        Eigen::SparseMatrix<double, Eigen::RowMajor> B(n, n);
        std::vector<Eigen::Triplet<double>> coefficients_A;
        std::vector<Eigen::Triplet<double>> coefficients_B;
        for (int i = 0; i < n; i++)
        {
            coefficients_A.push_back(Eigen::Triplet<double>(i, i, 3));
            coefficients_B.push_back(Eigen::Triplet<double>(i, i, 2));
            if (i > 0)
            {
                coefficients_A.push_back(Eigen::Triplet<double>(i - 1, i, 1));
                coefficients_B.push_back(Eigen::Triplet<double>(i - 1, i, 1));
            }
    
            if (i < n - 1)
            {
                coefficients_A.push_back(Eigen::Triplet<double>(i + 1, i, 1));
                coefficients_B.push_back(Eigen::Triplet<double>(i + 1, i, 1));
            }
        }
        A.setFromTriplets(coefficients_A.begin(), coefficients_A.end());
        B.setFromTriplets(coefficients_B.begin(), coefficients_B.end());
    
        // Construct matrix operation object using the wrapper classes
        Spectra::SparseSymMatProd<double> Aop(A);
        Spectra::SparseCholesky<double> Bop(B);
    
        // Construct generalized eigen solver object, requesting the smallest three generalized eigenvalues
        const int num_solution = 5;
        const int convergence_parameter = 8;
    
        Spectra::SymGEigsSolver<double, Spectra::SMALLEST_MAGN,
                                Spectra::SparseSymMatProd<double>,
                                Spectra::SparseCholesky<double>,
                                Spectra::GEIGS_CHOLESKY>
            eigen_solver(&Aop, &Bop, num_solution, convergence_parameter);
    
        // Initialize and compute
        eigen_solver.init();
        int num_convergence = eigen_solver.compute();
        std::cout << "num convergence = " << num_convergence << std::endl;
    
        // Retrieve results
        Eigen::VectorXd eigen_value;
        Eigen::MatrixXd eigen_vector;
        if (eigen_solver.info() == Spectra::SUCCESSFUL)
        {
            eigen_value = eigen_solver.eigenvalues();
            eigen_vector = eigen_solver.eigenvectors();
        }
        std::cout << "Generalized eigenvalues found:\n"
                  << eigen_value << std::endl;
        std::cout << "Generalized eigenvectors found:\n"
                  << eigen_vector << std::endl;
    
        // Verify results using the generalized eigen solver in Eigen
        Eigen::MatrixXd Adense = A;
        Eigen::MatrixXd Bdense = B;
        Eigen::GeneralizedSelfAdjointEigenSolver<Eigen::MatrixXd> standard_eigen_solver(Adense, Bdense);
        std::cout << "Generalized eigenvalues:\n"
                  << standard_eigen_solver.eigenvalues() << std::endl;
        std::cout << "Generalized eigenvectors:\n"
                  << standard_eigen_solver.eigenvectors() << std::endl;
    
        return 0;
    }
    

    The follow script can used to run the project on my computer

    eigen_path=/opt/Eigen/eigen-3.3.7
    spectra_path=/opt/spectra/spectra-0.8.1/include/
    
    rm EigenSolverExample.out
    
    g++ -I ${eigen_path} -I ${spectra_path} EigenSolverExample.cpp -o EigenSolverExample.out
    
    ./EigenSolverExample.out
    

    One can see in code's top I defined two matrix

    Eigen::SparseMatrix<double, Eigen::RowMajor> A(n, n);
    Eigen::SparseMatrix<double, Eigen::RowMajor> B(n, n);
    

    If I modified to

    Eigen::SparseMatrix<double, Eigen::ColMajor> A(n, n);
    Eigen::SparseMatrix<double, Eigen::ColMajor> B(n, n);
    

    there will be no core dumped happen.

    Could any one help me? Thanks for your time.

    opened by qzxuhui 9
  • SymGEigsShiftSolver problem

    SymGEigsShiftSolver problem

    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

    opened by wo80 7
  • SymEigsSolver misses some eigenvalues

    SymEigsSolver misses some eigenvalues

    I am adapting the first example from the Quick Start page to compute eigenvalues of a graph Laplacian matrix. I'm using the code given at the bottom of the issue, which takes arguments n and k where n is the number of vertices in your graph (dimensions of the Laplacian matrix) and k is the number of eigenvalues to compute.

    The code:

    • constructs the normalised Laplacian of the cycle graph on n vertices
    • uses the example code from the Quick Start page to compute k eigenvalues

    When I call the program with different values of k, I get inconsistent results.

    >> ./demo 20 5
    Eigenvalues found:
          2
    1.95106
    1.80902
    1.58779
    1.30902
    
    >> ./demo 20 6
    Eigenvalues found:
          2
    1.95106
    1.95106
    1.80902
    1.80902
    1.58779
    

    Notice that the multiplicities of the eigenvalues seem to be ignored when asking for 5 eigenvalues, but the multiple eigenvalues are returned correctly when we ask for 6.

    Smallest Magnitude

    If I modify the code to look for the smallest magnitude eigenvalues, then the results are also inconsistent.

    >> ./demo 20 5
    Eigenvalues found:
            1
     0.690983
     0.412215
     0.190983
    0.0489435
    
    >> ./demo 20 6
    Eigenvalues found:
        0.690983
        0.412215
        0.190983
       0.0489435
       0.0489435
    -9.62807e-23
    

    When asking for 5 eigenvalues, the smallest eigenvalue (which is 0) is not found, and multiplicities are ignored. When asking for 6 eigenvalues, the smallest eigenvalue is found although the multiplicities are still not correct - every eigenvalue other than 0 should have multiplicity 2.

    Code

    #include <Eigen/Core>
    #include <Spectra/SymEigsSolver.h>
    #include <iostream>
    #include <cstdlib>
    
    // Function declarations
    Eigen::MatrixXd construct_cycle_laplacian(int n);
    
    int main(int argc, char *argv[])
    {
      // Assume that we have been given two arguments: n and k
      //   n is the dimension of the matrix to construct
      //   k is the number of eigenvalues to compute
      int n = atoi(argv[1]);
      int k = atoi(argv[2]);
    
      // We are going to calculate the eigenvalues of M
      Eigen::MatrixXd M = construct_cycle_laplacian(n);
    
      // Construct matrix operation object using the wrapper class DenseSymMatProd
      Spectra::DenseSymMatProd<double> op(M);
    
      // Construct eigen solver object, requesting the largest k eigenvalues
      Spectra::SymEigsSolver<Spectra::DenseSymMatProd<double>> eigs(op, k, fmin(2*k, n));
    
      // Initialize and compute
      eigs.init();
      eigs.compute(Spectra::SortRule::LargestMagn);
    
      // Retrieve results
      Eigen::VectorXd eigenvalues;
      if (eigs.info() == Spectra::CompInfo::Successful) eigenvalues = eigs.eigenvalues();
    
      std::cout << "Eigenvalues found:\n" << eigenvalues << std::endl;
    
      return 0;
    }
    
    Eigen::MatrixXd construct_cycle_laplacian(int n) {
      // Initialise the laplacian matrix
      Eigen::MatrixXd L(n, n);
    
      // Add the matrix entries, iterating over the rows
      for (int i = 0; i < n; i++) {
        L(i, i) = 1;
        L(i, (i+n-1) % n) = -0.5;
        L(i, (i+1) % n) = -0.5;
      }
    
      return L;
    }
    
    opened by pmacg 0
  • Minimum error tolerance - GenEigsSolver

    Minimum error tolerance - GenEigsSolver

    Hi,

    Using GenEigsSolver, I calculated the eigenvalues of a large sparse matrix [6840 x 6840]. Theoretically, the real part of all eigenvalues of this matrix should be zero, but I got values mostly around 1.e-12. I think this might be related to an error tolerance set in GenEigsSolver, but I could not find any information regarding the tolerance.

    Q1. What is the minimum (or initial setting) error tolerance when using GenEigsSolver? Q2. Is there any way I can control the error tolerance?

    Sincerely,

    opened by yaeji43 0
  • More iterations are needed for some test in test/GenEigsRealShift.cpp

    More iterations are needed for some test in test/GenEigsRealShift.cpp

    Hello,

    I am working on packaging spectra 1.0.1 for Debian GNU/Linux.

    On my computer, when I run test/GenEigsRealShift.cpp through CMake, I get:


    Eigensolver of general real matrix [10x10] Smallest Real Part

    ./test/GenEigsRealShift.cpp:101 ...............................................................................

    ./test/GenEigsRealShift.cpp:61: FAILED: REQUIRE( eigs.info() == CompInfo::Successful ) with expansion: 2 == 0 with messages: nconv = 1 niter = 201 nops = 754

    I rose the number of iterations from 200 to 500 on line 44 of test/GenEigsRealShift.cpp, and then the test passed. I did not try with other values between 200 and 500. I am using Debian testing (Bookworm) on amd64.

    Best wishes,

    Pierre

    opened by pgrt 1
  • -Wuse-after-free with gcc 12.1.0

    -Wuse-after-free with gcc 12.1.0

    with gcc 12.1.0 and CMAKE_CXX_FLAGS="-Wall -DNDEBUG" spectra issues a -Wuse-after-free warning:

    [ 75%] Building CXX object test/CMakeFiles/GenEigs.dir/GenEigs.cpp.o
    In file included from /usr/include/eigen3/Eigen/Core:166,
                     from /home/schueller/projects/spectra/test/GenEigs.cpp:1:
    In function 'void Eigen::internal::aligned_free(void*)',
        inlined from 'void Eigen::internal::conditional_aligned_free(void*) [with bool Align = true]' at /usr/include/eigen3/Eigen/src/Core/util/Memory.h:259:15,
        inlined from 'void Eigen::internal::conditional_aligned_delete_auto(T*, std::size_t) [with T = double; bool Align = true]' at /usr/include/eigen3/Eigen/src/Core/util/Memory.h:446:34,
        inlined from 'Eigen::DenseStorage<T, -1, -1, _Cols, _Options>::~DenseStorage() [with T = double; int _Cols = 1; int _Options = 0]' at /usr/include/eigen3/Eigen/src/Core/DenseStorage.h:621:109,
        inlined from 'Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1> >::~PlainObjectBase()' at /usr/include/eigen3/Eigen/src/Core/PlainObjectBase.h:98:7,
        inlined from 'Eigen::Matrix<double, -1, 1>::~Matrix()' at /usr/include/eigen3/Eigen/src/Core/Matrix.h:178:7,
        inlined from 'void Spectra::UpperHessenbergEigen<Scalar>::doComputeEigenvectors() [with Scalar = double]' at /home/schueller/projects/spectra/include/Spectra/LinAlg/UpperHessenbergEigen.h:39:5:
    /usr/include/eigen3/Eigen/src/Core/util/Memory.h:203:9: warning: pointer may be used after 'void free(void*)' [-Wuse-after-free]
      203 |     free(ptr);
          |     ~~~~^~~~~
    In function 'void Eigen::internal::aligned_free(void*)',
        inlined from 'void Eigen::internal::conditional_aligned_free(void*) [with bool Align = true]' at /usr/include/eigen3/Eigen/src/Core/util/Memory.h:259:15,
        inlined from 'void Eigen::internal::conditional_aligned_delete_auto(T*, std::size_t) [with T = double; bool Align = true]' at /usr/include/eigen3/Eigen/src/Core/util/Memory.h:446:34,
        inlined from 'void Eigen::DenseStorage<T, -1, -1, _Cols, _Options>::resize(Eigen::Index, Eigen::Index, Eigen::Index) [with T = double; int _Cols = 1; int _Options = 0]' at /usr/include/eigen3/Eigen/src/Core/DenseStorage.h:637:77,
        inlined from 'void Eigen::PlainObjectBase<Derived>::resize(Eigen::Index, Eigen::Index) [with Derived = Eigen::Matrix<double, -1, 1>]' at /usr/include/eigen3/Eigen/src/Core/PlainObjectBase.h:285:25,
        inlined from 'static void Eigen::internal::Assignment<DstXprType, Eigen::Product<Lhs, Rhs, Options>, Eigen::internal::assign_op<Scalar, Scalar>, Eigen::internal::Dense2Dense, typename Eigen::internal::enable_if<((Options == Eigen::DefaultProduct) || (Options == Eigen::AliasFreeProduct))>::type>::run(DstXprType&, const SrcXprType&, const Eigen::internal::assign_op<Scalar, Scalar>&) [with DstXprType = Eigen::Matrix<double, -1, 1>; Lhs = Eigen::Block<Eigen::Matrix<double, -1, -1>, -1, -1, true>; Rhs = Eigen::Block<Eigen::Block<Eigen::Matrix<double, -1, -1>, -1, 1, true>, -1, 1, false>; int Options = 0; Scalar = double]' at /usr/include/eigen3/Eigen/src/Core/ProductEvaluators.h:146:17,
        inlined from 'void Eigen::internal::call_assignment_no_alias(Dst&, const Src&, const Func&) [with Dst = Eigen::Matrix<double, -1, 1>; Src = Eigen::Product<Eigen::Block<Eigen::Matrix<double, -1, -1>, -1, -1, true>, Eigen::Block<Eigen::Block<Eigen::Matrix<double, -1, -1>, -1, 1, true>, -1, 1, false>, 0>; Func = assign_op<double, double>]' at /usr/include/eigen3/Eigen/src/Core/AssignEvaluator.h:890:49,
        inlined from 'ExpressionType& Eigen::NoAlias<ExpressionType, StorageBase>::operator=(const StorageBase<OtherDerived>&) [with OtherDerived = Eigen::Product<Eigen::Block<Eigen::Matrix<double, -1, -1>, -1, -1, true>, Eigen::Block<Eigen::Block<Eigen::Matrix<double, -1, -1>, -1, 1, true>, -1, 1, false>, 0>; ExpressionType = Eigen::Matrix<double, -1, 1>; StorageBase = Eigen::MatrixBase]' at /usr/include/eigen3/Eigen/src/Core/NoAlias.h:43:31,
        inlined from 'void Spectra::UpperHessenbergEigen<Scalar>::doComputeEigenvectors() [with Scalar = double]' at /home/schueller/projects/spectra/include/Spectra/LinAlg/UpperHessenbergEigen.h:36:29:
    /usr/include/eigen3/Eigen/src/Core/util/Memory.h:203:9: note: call to 'void free(void*)' here
      203 |     free(ptr);
          |     ~~~~^~~~~
    [100%] Linking CXX executable GenEigs
    

    here is a minimal example to reproduce:

    #include <Spectra/GenEigsSolver.h>
    
    class KLMatProd
    {
    public:
      using Scalar = double;
    
      KLMatProd()
      {
      }
      
      virtual ~KLMatProd() {}
      
      virtual int rows() const
      {
        return 0;
      }
    
      virtual int cols() const
      {
        return 0;
      }
    
      virtual void perform_op(const Scalar * x_in, Scalar * y_out) const
      {
        (void)x_in;
        (void)y_out;
      }
    };
    
    int main()
    {
      KLMatProd op;
      Spectra::GenEigsSolver<KLMatProd> solver(op, 8,3);
      solver.init();
      solver.compute(Spectra::SortRule::LargestMagn);
      return 0;
    }
    
    

    in UpperHessenbergEigen::doComputeEigenvectors, if I remove the ".noalias()" below the warning goes away:

           // Back transformation to get eigenvectors of original matrix
            Vector m_tmp(size);
            for (Index j = size - 1; j >= 0; j--)
            {
                m_tmp.noalias() = m_eivec.leftCols(j + 1) * m_matT.col(j).segment(0, j + 1);
                m_eivec.col(j) = m_tmp;
            }
    

    I dont know if it"s an actual spectra issue though

    opened by jschueller 0
  • Eigenvalues test is different between macos and linux for 1.0.1 release

    Eigenvalues test is different between macos and linux for 1.0.1 release

    It looks like the Eigenvalues behavior is different between linux and macos for 1.0.1 release.

    Full build log ref, https://github.com/Homebrew/homebrew-core/runs/6248946918 relates to https://github.com/Homebrew/homebrew-core/pull/98697

    opened by chenrui333 1
  • Include Spectra as a subdirectory in a CMake build

    Include Spectra as a subdirectory in a CMake build

    Summary

    This PR adds support to include Spectra as a subdirectory project via the add_subdirectory call in CMake. For example, the following should add the Spectra::Spectra target to the build tree so it can be linked to further down the line without the need for installation:

    # ...
    add_subdirectory(third-party/spectra)
    find_package (Spectra 1.0.1 REQUIRED)
    add_library(MyProject main.cpp)
    target_link_libraries(MyProject PUBLIC Spectra::Spectra)
    

    assuming spectra is in the directory third-party/. This should enable Spectra to be utilized as a git submodule.

    Notes

    The installation behaviour should have remained unchanged. I also added some stuff to the .gitignore.

    opened by cbritopacheco 0
Releases(v1.0.1)
  • v1.0.1(Apr 6, 2022)

    Added

    • Added SIMD support for UpperHessenbergSchur. This should accelerate general eigen solvers such as GenEigsSolver
    • Added test code for UpperHessenbergSchur

    Changed

    • Fixed several bugs in the examples caused by the const keyword, reported by @alexpghayes (#135, #137)
    • Updated the included Catch2 to v2.13.8
    Source code(tar.gz)
    Source code(zip)
  • v1.0.0(Jul 1, 2021)

    Finally, Spectra v1.0.0 is released!

    This new version is the accumulation of dozens of new features and bug fixes, with the help of many users and contributors. The most visible changes compared to v0.9.0 include fully refactored code base using C++11, a redesigned API, and the newly added Davidson eigen solver by Jens Wehner and his group. More details can be found below and in the migration guide.

    Added

    • Added version macros SPECTRA_MAJOR_VERSION, SPECTRA_MINOR_VERSION, SPECTRA_PATCH_VERSION, and SPECTRA_VERSION that are included by all eigen solvers
    • Added the wrapper class SparseGenComplexShiftSolve for eigen solver with complex shifts
    • Added the SymGEigsShiftSolver class for symmetric generalized eigen solver with real shifts
    • Added the wrapper class SymShiftInvert that can be used with SymGEigsShiftSolver
    • Added test code for symmetric generalized eigen solver with real shifts
    • Added an internal class UpperHessenbergSchur to compute the Schur decomposition of upper Hessenberg matrices more efficiently
    • Added a Flags template parameter to every matrix operation class (e.g. DenseCholesky and DenseSymMatProd), whose possible values are Eigen::ColMajor and Eigen::RowMajor. This parameter allows these wrapper classes to handle row-major matrices. If the input matrix is inconsistent with the Flags parameter (e.g., if Flags is Eigen::ColMajor but the input matrix is row-major), a compiler error will occur
    • Added the member function info() and convergence tests to SparseRegularInverse, suggested by @Spammed (#111)
    • Added symmetric Davidson eigen solver DavidsonSymEigsSolver, written by Felipe Zapata, Nicolas Renaud, Victor Azizi, Pablo Lopez-Tarifa, and Jens Wehner from the Netherlands eScience Center
    • Extended matrix operations in DenseGenMatProd, DenseSymMatProd, SparseGenMatProd, and SparseSymMatProd to handle matrix-matrix products and coefficient-wise accessors

    Changed

    • API change: Spectra now requires C++11
    • API change: All enumerations have been converted to enum classes (e.g. LARGEST_MAGN is now SortRule::LargestMagn)
    • API change: Selection rules are no longer template parameters. They are now specified in the compute() member function as arguments
    • API change: The Scalar template parameter has been removed from eigen solvers. Instead, matrix operation classes now need to define a public type named Scalar
    • API change: Constructors of solvers now request references of matrix operators instead of pointers
    • Clang-Format now uses the C++11 standard to format code
    • Updated documentation to reflect the new API
    • Many internal changes to make use of C++11 features
    • Added a SPECTRA_ prefix to each header guard to prevent potential name clash
    • Changed the default value of the Flags template parameter that exists in various class templates from 0 to the more readable constant Eigen::ColMajor
    • Renamed the function mat_prod to perform_op in the SparseRegularInverse wrapper class. This makes the API more consistent when implementing new generalized eigen solvers
    • Improved the precision of UpperHessenbergQR and TridiagQR by computing the Givens rotations in a more stable way
    • Added a deflation test to TridiagQR to accelerate the convergence of eigen solvers
    • Improved the precision of TridiagQR::matrix_QtHQ() by directly applying rotations to the original input matrix
    • Improved the precision of DoubleShiftQR by computing the Householder reflectors in a more stable way
    • Improved the deflation test in DoubleShiftQR
    • More careful computation of residual vectors in the Lanczos process
    • Initial vectors in the Lanczos and Arnoldi processes are now forced to be in the range of the A matrix
    • More sensible test for orthogonality in generating new random vectors in the Lanczos and Arnoldi processes
    • In symmetric eigen solvers large shifts are applied first to increase precision
    • Updated the included Catch2 to v2.13.6
    Source code(tar.gz)
    Source code(zip)
  • v0.9.0(May 19, 2020)

  • v0.8.1(Jun 7, 2019)

    Changed

    • Fixed a bug in BKLDLT in which a wrong type was used, thanks to @jdbancal for the issue #64
    • Fixed a bug in BKLDLT that caused segmentation fault in some edge cases, also reported by @jdbancal in issue #66
    • The type Eigen::Index is now globally used for indices and sizes, in order to handle potentially large matrices. This was suggested by Yuan Yao in issue #19
    Source code(tar.gz)
    Source code(zip)
  • v0.8.0(Apr 4, 2019)

    Added

    • Added a BKLDLT class that implements the Bunch-Kaufman LDLT decomposition for symmetric indefinite matrices. According to the Eigen documentation, currently Eigen::LDLT cannot handle some special indefinite matrices such as [0, 1; 1, 0], but BKLDLT is applicable to any symmetric matrices as long as it is not singular. LDLT decomposition is used in shift-and-invert solvers (see below)
    • Added a unit test for BKLDLT

    Changed

    • DenseSymShiftSolve now uses the newly added BKLDLT class to do the decomposition. This change broadens the class of matrices that DenseSymShiftSolve can handle, reduces memory use, and should also improve the numerical stability of the solver
    • Replaced Eigen::SimplicialLDLT with Eigen::SparseLU in the SparseSymShiftSolve class, as some edge-case indefinite matrices may break Eigen::SimplicialLDLT
    • SparseSymShiftSolve and SparseGenRealShiftSolve will throw an error if the factorization failed, for example, on singular matrices
    • Fixed a missing #include in DenseCholesky.h, thanks to Lennart Trunk for the issue #59
    • Fixed errors in examples (#60), thanks to @linuxfreebird
    • Updated the included Catch2 to v2.7.0
    Source code(tar.gz)
    Source code(zip)
  • v0.7.0(Jan 10, 2019)

    Added

    • Added a directory contrib to include code contributed by users. It is not formally a part of the Spectra library, but it may contain useful solvers and applications based on Spectra. Code in contrib may not be fully tested, so please use with caution. Feedback and report of issues are always welcome
    • Added an eigen solver LOBPCGSolver in the contrib directory using the LOBPCG algorithm, contributed by Anna Araslanova
    • Added a partial SVD solver PartialSVDSolver in the contrib directory
    • Added two internal classes Arnoldi and Lanczos to compute the Arnoldi/Lanczos factorization in eigen solvers
    • Added a few other internal classes to refactor the eigen solver classes (see below)

    Changed

    • API change: Spectra now requires Eigen >= 3.3
    • API change: The library header files are moved into a directory named Spectra. Hence the recommended include directive would look like #include <Spectra/SymEigsSolver.h>
    • All eigen solvers have been refactored using a cleaner class hierarchy. It may potentially make the implementation of new eigen solvers easier, especially for generalized eigen problems
    • The matrix operation classes (e.g. DenseSymMatProd and SparseSymMatProd) are now internally using an Eigen::Ref object to wrap the user matrices, thanks to Dario Mangoni who raised this issue in #16
    • Fixed inappropriate range of random numbers in the tests
    Source code(tar.gz)
    Source code(zip)
  • v0.6.2(May 22, 2018)

    Changed

    • Fixed regressions in v0.6.0 on some edge cases
    • Improved the accuracy of restarting processes in SymEigsSolver and GenEigsSolver
    • Updated the included Catch2 to v2.2.2
    • Code and documentation cleanup
    Source code(tar.gz)
    Source code(zip)
  • v0.6.1(Mar 4, 2018)

  • v0.6.0(Mar 3, 2018)

    Added

    • Added virtual destructors to the SymEigsSolver and UpperHessenbergQR classes to fix compiler warnings, by Julian Kent
    • Added a NUMERICAL_ISSUE entry to the COMPUTATION_INFO enumeration to indicate the status of Cholesky decomposition
    • Added the info() member function to DenseCholesky and SparseCholesky to report the status of the decomposition
    • Added a missing #include item in SparseCholesky.h, thanks to Maxim Torgonsky
    • Added a TypeTraits class to retrieve additional numeric limits of scalar value types

    Changed

    • Documentation updates
    • Updated the project URL to https://spectralib.org
    • Some internal improvements, such as pre-allocating vectors in loops, and changing return type to reference, thanks to Angelos Mantzaflaris
    • Improved the accuracy of symmetric and general eigen solvers
    • Reduced the memory use of UpperHessenbergQR and TridiagQR decompositions
    • Updated the included Catch2 to v2.0.1
    • Updated the testing code using the new API of Catch2
    • Updated Travis CI script
    Source code(tar.gz)
    Source code(zip)
  • v0.5.0(Feb 5, 2017)

    Added

    • Added the generalized eigen solver SymGEigsSolver in the regular inverse mode
    • Added the wrapper class SparseRegularInverse that can be used with SymGEigsSolver in the regular inverse mode
    • Added test code for generalized eigen solver in the regular inverse mode

    Changed

    • Improved the numerical precision and stability of some internal linear algebra classes, including TridiagEigen, UpperHessenbergEigen, and DoubleShiftQR
    • API change: The x_in argument in matrix operation functions, e.g. perform_op(), is now labelled to be constant
    • Fixed a bug that GenEigsComplexShiftSolver gave wrong results when transforming back the eigenvalues, discovered by @jdbancal
    • Updated included Catch to v1.7.0
    • Documentation improvement
    Source code(tar.gz)
    Source code(zip)
  • v0.4.0(Nov 15, 2016)

    Added

    • Added an Uplo template parameter to the DenseSymShiftSolve class
    • Added the generalized eigen solver SymGEigsSolver in Cholesky mode
    • Added the wrapper classes DenseCholesky and SparseCholesky that can be used in SymGEigsSolver
    • Added test code for generalized eigen solver

    Changed

    • Allowing basic math functions such as abs() and sqrt() to be overloaded (avoid using std::abs and std::sqrt directly), thanks to @jdbancal. This makes it possible to use user-defined float number types with Spectra
    • Replaced other std functions by their Eigen counterparts, for example using Eigen::NumTraits<Scalar>::epsilon() to substitute std::numeric_limits<Scalar>::epsilon()
    • Fixed an out-of-bound bug detected by @jdbancal
    • Improved numerical stability, e.g. the function hypot(x, y) is used to compute sqrt(x^2 + y^2)
    • More careful in using "approximate zero" constants
    • Updated included Catch to v1.5.7
    • Improved documentation
    • Updated Travis CI script
    Source code(tar.gz)
    Source code(zip)
  • v0.3.0(Jul 3, 2016)

    Added

    • Added the wrapper classes SparseSymMatProd and SparseSymShiftSolve for sparse symmetric matrices
    • Added the wrapper class SparseGenRealShiftSolve for general sparse matrices
    • Added tests for sparse matrices
    • Using Travis CI for automatic unit test

    Changed

    • Updated included Catch to v1.5.6
    • API change: Each eigen solver was moved to its own header file. For example to use SymEigsShiftSolver one needs to include <SymEigsShiftSolver.h>
    • Header files for internal use were relocated
    Source code(tar.gz)
    Source code(zip)
  • v0.2.0(Feb 28, 2016)

    Added

    • Benchmark script now outputs number of matrix operations
    • Added a simple built-in random number generator, so that the algorithm was made to be deterministic
    • Added the wrapper class DenseSymMatProd for symmetric matrices

    Changed

    • Improved Arnoldi factorization
      • Iteratively corrects orthogonality
      • Creates new residual vector when invariant subspace is found
      • Stability for matrices with repeated eigenvalues is greatly improved
    • Adjusted deflation tolerance in double shift QR
    • Updated result analyzer
    • Updated copyright information
    • API change: Default operator of SymEigsSolver was changed from DenseGenMatProd to DenseSymMatProd
    Source code(tar.gz)
    Source code(zip)
  • v0.1.0(Dec 19, 2015)

C++ library for solving large sparse linear systems with algebraic multigrid method

AMGCL AMGCL is a header-only C++ library for solving large sparse linear systems with algebraic multigrid (AMG) method. AMG is one of the most effecti

Denis Demidov 552 Sep 30, 2022
A C++ header-only library of statistical distribution functions.

StatsLib StatsLib is a templated C++ library of statistical distribution functions, featuring unique compile-time computing capabilities and seamless

Keith O'Hara 404 Sep 28, 2022
RcppFastFloat: Rcpp Bindings for the fastfloat C++ Header-Only Library

Converting ascii text into (floating-point) numeric values is a very common problem. The fast_float header-only C++ library by Daniel Lemire does this very well, and very fast at up to or over to 1 gigabyte per second as described in more detail in a recent arXiv paper.

Dirk Eddelbuettel 18 May 2, 2022
Header only FFT library

dj_fft: Header-only FFT library Details This repository provides a header-only library to compute fourier transforms in 1D, 2D, and 3D. Its goal is to

Jonathan Dupuy 121 Sep 1, 2022
C++ header-only fixed-point math library

fpm A C++ header-only fixed-point math library. "fpm" stands for "fixed-point math". It is designed to serve as a drop-in replacement for floating-poi

Mike Lankamp 348 Oct 2, 2022
C++ header-only library with methods to efficiently encode/decode Morton codes in/from 2D/3D coordinates

Libmorton v0.2.7 Libmorton is a C++ header-only library with methods to efficiently encode/decode 64, 32 and 16-bit Morton codes and coordinates, in 2

Jeroen Baert 473 Sep 27, 2022
Extremely simple yet powerful header-only C++ plotting library built on the popular matplotlib

matplotlib-cpp Welcome to matplotlib-cpp, possibly the simplest C++ plotting library. It is built to resemble the plotting API used by Matlab and matp

Benno Evers 3.4k Oct 4, 2022
A modern, C++20-native, single-file header-only dense 2D matrix library.

A modern, C++20-native, single-file header-only dense 2D matrix library. Contents Example usage creating matrices basic operations row, col, size, sha

feng wang 52 Sep 26, 2022
libmpc++ is a C++ header-only library to solve linear and non-linear MPC

libmpc++ libmpc++ is a C++ library to solve linear and non-linear MPC. The library is written in modern C++17 and it is tested to work on Linux, macOS

Nicola Piccinelli 39 Jul 8, 2022
Header-only C++11 library to handle physical measures

cpp-measures Header-only C++11 library to handle physical measures License: This project is released under the Mozilla Public License 2.0. Purpose Thi

Carlo Milanesi 20 Jun 28, 2018
Header only, single file, simple and efficient C++ library to compute the signed distance function to a triangle mesh

TriangleMeshDistance Header only, single file, simple and efficient C++11 library to compute the signed distance function to a triangle mesh. The dist

Interactive Computer Graphics 95 Sep 23, 2022
A work-in-progress C++20/23 header-only maths library for game development, embedded, kernel and general-purpose that works in constant context.

kMath /kmæθ/ A work-in-progress general-purpose C++20/23 header-only maths library that works in constant context Abstract The kMath Project aims to p

The λ Project 13 Sep 5, 2022
A C++ header only library for decomposition of spectra into a sum of response functions whose weights are positive definite.

DecompLib A C++ header only library for decomposition of spectra into a sum of response functions whose weights are positive definite. Introduction Bu

James T. Matta 2 Jan 22, 2022
A matrix header-only library, uses graphs internally, helpful when your matrix is part of a simulation where it needs to grow many times (or auto expand)

GraphMat Header-only Library Matrix implemented as a graph, specially for the use case when it should be auto expanding at custom rate, specially in s

Aditya Gupta 3 Oct 25, 2021
linalg.h is a single header, public domain, short vector math library for C++

linalg.h linalg.h is a single header, public domain, short vector math library for C++. It is inspired by the syntax of popular shading and compute la

Sterling Orsten 746 Oct 1, 2022
MIRACL Cryptographic SDK: Multiprecision Integer and Rational Arithmetic Cryptographic Library is a C software library that is widely regarded by developers as the gold standard open source SDK for elliptic curve cryptography (ECC).

MIRACL What is MIRACL? Multiprecision Integer and Rational Arithmetic Cryptographic Library – the MIRACL Crypto SDK – is a C software library that is

MIRACL 498 Sep 28, 2022
A C library for statistical and scientific computing

Apophenia is an open statistical library for working with data sets and statistical or simulation models. It provides functions on the same level as t

null 186 Sep 11, 2022
P(R*_{3, 0, 1}) specialized SIMD Geometric Algebra Library

Klein ?? ?? Project Site ?? ?? Description Do you need to do any of the following? Quickly? Really quickly even? Projecting points onto lines, lines t

Jeremy Ong 621 Sep 30, 2022
LibTomMath is a free open source portable number theoretic multiple-precision integer library written entirely in C.

libtommath This is the git repository for LibTomMath, a free open source portable number theoretic multiple-precision integer (MPI) library written en

libtom 531 Sep 29, 2022