Spectra

https://github.com/yixuan/spectra
Admin
Spectra is a header-only C++ library for large-scale eigenvalue problems, offering a redesigned...
Tokens:5,602
Snippets:54
Trust Score:9.7
License:MPL-2.0
Update:4 months ago
Show doc for...
Tokens:
Raw
### Build and Run Example

Source: https://github.com/yixuan/spectra/blob/master/examples/DavidsonSymEigs_example.md

Steps to compile and execute the Davidson Symmetric Eigenvalue Solver example.

```bash
mkdir build && cd build && cmake ../
make DavidsonSymEigs_example
./example/DavidsonSymEigs_example
```

--------------------------------

### Install Spectra using CMake

Source: https://github.com/yixuan/spectra/blob/master/README.md

Use these commands to build and install Spectra with CMake. Ensure CMake v3.10 or higher is installed. You can specify the installation directory and enable tests.

```bash
mkdir build && cd build
cmake .. -DCMAKE_INSTALL_PREFIX='intended installation directory' -DBUILD_TESTS=TRUE
make all && make test && make install
```

--------------------------------

### Configure Eigen Installation with CMake

Source: https://github.com/yixuan/spectra/blob/master/README.md

When Eigen is already installed, configure the CMake build to find it by setting CMAKE_PREFIX_PATH or Eigen3_ROOT. This example shows how to specify the installation directory and enable tests.

```bash
cmake .. -DCMAKE_INSTALL_PREFIX='intended installation directory' -DCMAKE_PREFIX_PATH='path where the installation of Eigen3 can be found' -DBUILD_TESTS=TRUE
```

--------------------------------

### CMakeLists.txt Configuration

Source: https://github.com/yixuan/spectra/blob/master/examples/CMakeLists.txt

This snippet shows how to set up a CMake build for C++ examples that use the Spectra library. It defines source files, creates executables, and links the necessary library.

```cmake
set(example_target_sources)

list(APPEND example_target_sources
        DavidsonSymEigs_example.cpp
        )

foreach(EXAMPLE_SOURCE ${example_target_sources})
    # Extract the filename without extension (NAME_WE) as a name for our executable
    get_filename_component(EXAMPLE_NAME ${EXAMPLE_SOURCE} NAME_WE)

    # Add an executable based on the source
    add_executable(${EXAMPLE_NAME} ${EXAMPLE_SOURCE})

    # Configure (include headers and link libraries) the example
    target_link_libraries(${EXAMPLE_NAME} PRIVATE Spectra)

endforeach()
```

--------------------------------

### Eigen Solver Migration: 0.9.0 to 1.0.0

Source: https://github.com/yixuan/spectra/blob/master/MIGRATION.md

Illustrates the conversion of eigen solver code from Spectra 0.9.0 to 1.0.0, reflecting changes in template parameters, runtime selection rules, and constructor arguments.

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

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

// Initialize, and compute with at most 1000 iterations
eigs.init();
int nconv = eigs.compute(1000);

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

```

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

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

// Initialize, and compute with at most 1000 iterations
eigs.init();
int nconv = eigs.compute(SortRule::LargestAlge, 1000);

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

```

--------------------------------

### Initializing the Davidson Symmetric Eigenvalue Solver

Source: https://github.com/yixuan/spectra/blob/master/examples/DavidsonSymEigs_example.md

Constructs the DavidsonSymEigsSolver, specifying the matrix operation and the number of eigenpairs to find. Internal parameters like search space size can be tuned.

```cpp
Eigen::Index num_of_eigenvalues = 5;
DavidsonSymEigsSolver<DenseSymMatProd<double>> solver(op_dense, num_of_eigenvalues);  //Create Solver

```

--------------------------------

### Creating and Configuring Test Executables in CMake

Source: https://github.com/yixuan/spectra/blob/master/test/CMakeLists.txt

This CMake loop iterates through a list of test source files, creating an executable for each. It configures each test executable by linking it to the Spectra library and a main test runner, and sets the working directory for running the tests.

```cmake
foreach(TEST_SOURCE ${test_target_sources})
    # Extract the filename without extension (NAME_WE) as a name for our executable
    get_filename_component(TEST_NAME ${TEST_SOURCE} NAME_WE)

    # Add an executable based on the source
    add_executable(${TEST_NAME} ${TEST_SOURCE})

    # Configure (include headers and link libraries) the test
    target_link_libraries(${TEST_NAME} PRIVATE Spectra tests-main)

    add_test(NAME ${TEST_NAME} COMMAND ${TEST_NAME}
            WORKING_DIRECTORY ${CMAKE_BINARY_DIR})  # the working directory is the out-of-source build directory
endforeach()
```

--------------------------------

### Constructing a Diagonally Dominant Symmetric Matrix

Source: https://github.com/yixuan/spectra/blob/master/examples/DavidsonSymEigs_example.md

Initializes a large random symmetric matrix with diagonal dominance for eigenvalue computation. Ensure the matrix is diagonally dominant for convergence.

```cpp
#include <Eigen/Dense>
#include <Spectra/DavidsonEigsSolver.h>
#include <Spectra/MatOp/DenseSymMatProd.h> 
#include <iostream>

using namespace Spectra;

int main()
{
Eigen::Index n = 1000;
    Eigen::MatrixXd mat = 0.03 * Eigen::MatrixXd::Random(n, n);
    Eigen::MatrixXd mat1 = mat + mat.transpose();
    for (Eigen::Index i=0; i<n; i++) {
        mat1(i,i) += i+1;
    }

```

--------------------------------

### Custom Matrix Operation Eigen Solver

Source: https://github.com/yixuan/spectra/blob/master/doxygen/Overview.md

Demonstrates using a user-supplied matrix operation class with SymEigsSolver. The custom class must define a 'Scalar' typedef and 'rows()', 'cols()', and 'perform_op()' methods.

```cpp
#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() const { return 10; }
    int cols() const { 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;
}

```

--------------------------------

### Complex-Valued Matrix Eigenvalue Solvers

Source: https://github.com/yixuan/spectra/blob/master/doxygen/Overview.md

Demonstrates the usage of Spectra's HermEigsSolver for complex Hermitian matrices and GenEigsSolver for general complex matrices. Ensure necessary headers are included and matrix operation objects are correctly constructed.

```cpp
#include <Eigen/Core>
#include <Spectra/HermEigsSolver.h>
#include <Spectra/GenEigsSolver.h>
#include <iostream>

using namespace Spectra;

int main()
{
    std::srand(0);

    // We are going to calculate the eigenvalues of H and G
    Eigen::MatrixXcd G = Eigen::MatrixXcd::Random(10, 10);
    // H is Hermitian
    Eigen::MatrixXcd H = G + G.adjoint();

    // Construct matrix operation objects using the wrapper
    // classes DenseHermMatProd and DenseGenMatProd
    using OpHType = DenseHermMatProd<std::complex<double>>;
    using OpGType = DenseGenMatProd<std::complex<double>>;
    OpHType opH(H);
    OpGType opG(G);

    // Construct solver object for H, requesting the largest three eigenvalues
    HermEigsSolver<OpHType> eigsH(opH, 3, 6);

    // Initialize and compute
    eigsH.init();
    int nconvH = eigsH.compute(SortRule::LargestAlge);

    // Retrieve results
    // Eigenvalues are real-valued, and eigenvectors are complex-valued
    if (eigsH.info() == CompInfo::Successful)
    {
        Eigen::VectorXd evaluesH = eigsH.eigenvalues();
        std::cout << "Eigenvalues of H found:\n" << evaluesH << std::endl;
        Eigen::MatrixXcd evecsH = eigsH.eigenvectors();
        std::cout << "Eigenvectors of H:\n" << evecsH << std::endl;
    }

    // Similar procedure for matrix G
    GenEigsSolver<OpGType> eigsG(opG, 3, 6);
    eigsG.init();
    int nconvG = eigsG.compute(SortRule::LargestMagn);
    if (eigsG.info() == CompInfo::Successful)
    {
        Eigen::VectorXcd evaluesG = eigsG.eigenvalues();
        std::cout << "Eigenvalues of G found:\n" << evaluesG << std::endl;
        Eigen::MatrixXcd evecsG = eigsG.eigenvectors();
        std::cout << "Eigenvectors of G:\n" << evecsG << std::endl;
    }

    return 0;
}

```

--------------------------------

### Creating a Dense Symmetric Matrix Product Operation

Source: https://github.com/yixuan/spectra/blob/master/examples/DavidsonSymEigs_example.md

Defines the matrix operation required by the solver using Spectra's DenseSymMatProd. This allows using custom matrix operations without needing the underlying matrix explicitly.

```cpp
DenseSymMatProd<double> op(mat); // Create the Matrix Product operation

```

--------------------------------

### Computing Eigenpairs and Retrieving Results

Source: https://github.com/yixuan/spectra/blob/master/examples/DavidsonSymEigs_example.md

Executes the solver using the compute method, specifying the sorting rule, maximum iterations, and tolerance. It retrieves and prints the converged eigenvalues or reports failure.

```cpp
    Eigen::Index iterations = 100;
    double tolerance = 1e-3;
    int nconv = solver.compute(SortRule::LargestAlge, iterations, tolerance);

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

    std::cout <<nconv<< " Eigenvalues found:\n"
              << evalues << std::endl;
    }else{
        std::cout <<"Calculation failed"<<std::endl;
    }
    return 0;
}

```

--------------------------------

### Dense Symmetric Eigen Solver

Source: https://github.com/yixuan/spectra/blob/master/doxygen/Overview.md

Calculates eigenvalues for a dense symmetric matrix using SymEigsSolver. Requires Eigen and Spectra libraries.

```cpp
#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;
}

```

--------------------------------

### Defining Test Sources in CMake

Source: https://github.com/yixuan/spectra/blob/master/test/CMakeLists.txt

This snippet defines a list of C++ source files used for testing various linear algebra and eigen-solver functionalities within the Spectra project. It categorizes sources for clarity.

```cmake
set(test_target_sources)

add_library (tests-main tests-main.cpp)
list(APPEND test_target_sources
        # Linear algebra
        Givens.cpp
        QR.cpp
        Eigen.cpp
        Schur.cpp
        BKLDLT.cpp
        Arnoldi.cpp

        # JD linear algebra
        Orthogonalization.cpp
        RitzPairs.cpp
        SearchSpace.cpp

        # Matrix operators
        DenseGenMatProd.cpp
        DenseSymMatProd.cpp
        SparseGenMatProd.cpp
        SparseSymMatProd.cpp

        # Symmetric eigen solver
        SymEigs.cpp
        SymEigsShift.cpp

        # Hermitian eigen solver
        HermEigs.cpp

        # General eigen solver
        GenEigs.cpp
        GenEigsRealShift.cpp
        GenEigsComplexShift.cpp

        # Complex eigen solver
        ComplexEigs.cpp

        # Symmetric generalized eigen solver
        SymGEigsCholesky.cpp
        SymGEigsRegInv.cpp
        SymGEigsShift.cpp

        # SVD
        SVD.cpp

        # JD eigen solver
        JDSymEigsBase.cpp
        JDSymEigsDPRConstructor.cpp
        DavidsonSymEigs.cpp

        # Examples from bug reports
        Example1.cpp
        Example2.cpp
        Example3.cpp
        Example4.cpp
        )
```

--------------------------------

### Sparse General Eigen Solver

Source: https://github.com/yixuan/spectra/blob/master/doxygen/Overview.md

Calculates eigenvalues for a sparse matrix using GenEigsSolver. Supports sparse matrices via SparseGenMatProd.

```cpp
#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;
}

```

--------------------------------

### Add Scalar Type to User-Defined Matrix Class

Source: https://github.com/yixuan/spectra/blob/master/MIGRATION.md

When defining custom matrix operation classes, ensure a public type definition named 'Scalar' is included to specify the element type of the matrix.

```cpp
// A user-defined matrix operation class
// representing the matrix A=diag(1, 2, ..., 10)
class MyDiagonalTen
{
public:
    // The line below is new
    using Scalar = double;

    int rows() const { return 10; }
    int cols() const { 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);
        }
    }
};

```
Tokens:
Raw