generalizedeigensolver.hh

Go to the documentation of this file.

// SPDX-FileCopyrightText: 2021-2025 The Ikarus Developers ikarus@ibb.uni-stuttgart.de
// SPDX-License-Identifier: LGPL-3.0-or-later
 
#pragma once
 
#include <optional>
#include <Spectra/MatOp/SparseCholesky.h>
#include <Spectra/MatOp/SparseGenMatProd.h>
#include <Spectra/SymEigsSolver.h>
#include <Spectra/SymGEigsSolver.h>
#include <utility>
 
#include <Eigen/Core>
#include <Eigen/SparseCore>
 
#include <ikarus/utils/concepts.hh>
#include <ikarus/utils/makeenum.hh>
 
namespace Ikarus {
 
MAKE_ENUM(EigenValueSolverType, Spectra, Eigen);
 
template <EigenValueSolverType SolverType, Concepts::DenseOrSparseEigenMatrix MT>
struct GeneralizedSymEigenSolver
{
};
 
template <Concepts::DenseOrSparseEigenMatrix MT>
struct GeneralizedSymEigenSolver<EigenValueSolverType::Spectra, MT>
{
  using ScalarType              = typename MT::Scalar;
  static constexpr bool isDense = Concepts::EigenMatrix<MT>;
 
  using MatrixType = MT;
  using ProductType =
      std::conditional_t<isDense, Spectra::DenseSymMatProd<ScalarType>, Spectra::SparseSymMatProd<ScalarType>>;
  using CholeskyType =
      std::conditional_t<isDense, Spectra::DenseCholesky<ScalarType>, Spectra::SparseCholesky<ScalarType>>;
 
  using SolverType = Spectra::SymGEigsSolver<ProductType, CholeskyType, Spectra::GEigsMode::Cholesky>;
 
  template <typename MATA, typename MATB>
  requires(Concepts::DenseOrSparseEigenMatrix<std::remove_cvref_t<MATA>>)
  GeneralizedSymEigenSolver(MATA&& A, MATB&& B)
      : nev_(A.rows()),
        nevsPartition_(static_cast<Eigen::Index>(std::ceil(static_cast<double>(nev_) / 2)),
                       static_cast<Eigen::Index>(nev_ - std::ceil(static_cast<double>(nev_) / 2))),
        aOP_(std::forward<MATA>(A)),
        bOP_(std::forward<MATB>(B)),
        solverSmallest_(aOP_, bOP_, nevsPartition_.first, std::min(nevsPartition_.second * 2, nev_)),
        solverGreatest_(aOP_, bOP_, nevsPartition_.second, nevsPartition_.second * 2) {
    if ((A.cols() != B.cols()) or (A.rows() != B.rows()) or (A.cols() != A.rows()))
      DUNE_THROW(Dune::InvalidStateException,
                 "GeneralizedSymEigenSolver: The passed matrices should have the same size");
    eigenvalues_.resize(nev_);
    eigenvectors_.resize(A.rows(), nev_);
    assert(nevsPartition_.first + nevsPartition_.second == nev_);
  }
 
  template <Concepts::FlatAssembler AssemblerA, Concepts::FlatAssembler AssemblerB>
  GeneralizedSymEigenSolver(std::shared_ptr<AssemblerA> assemblerA, const std::shared_ptr<AssemblerB>& assemblerB)
      : GeneralizedSymEigenSolver(assemblerA->matrix(), assemblerB->matrix()) {}
 
  bool compute(Eigen::Index maxit = 1000, ScalarType tolerance = 1e-10) {
    solverSmallest_.init();
    solverSmallest_.compute(Spectra::SortRule::SmallestAlge, maxit, tolerance, Spectra::SortRule::SmallestAlge);
 
    solverGreatest_.init();
    solverGreatest_.compute(Spectra::SortRule::LargestAlge, maxit, tolerance, Spectra::SortRule::SmallestAlge);
 
    computed_ = solverSmallest_.info() == Spectra::CompInfo::Successful and
                solverGreatest_.info() == Spectra::CompInfo::Successful;
    if (not computed_)
      return false;
 
    eigenvalues_.head(nevsPartition_.first)      = solverSmallest_.eigenvalues();
    eigenvectors_.leftCols(nevsPartition_.first) = solverSmallest_.eigenvectors();
 
    eigenvalues_.tail(nevsPartition_.second)       = solverGreatest_.eigenvalues();
    eigenvectors_.rightCols(nevsPartition_.second) = solverGreatest_.eigenvectors();
 
    return computed_;
  }
 
  auto& eigenvalues() const {
    assertCompute();
    return eigenvalues_;
  }
 
  auto& eigenvectors() const {
    assertCompute();
    return eigenvectors_;
  }
 
  Eigen::Index nev() const { return nev_; }
 
private:
  Eigen::Index nev_;
  std::pair<Eigen::Index, Eigen::Index> nevsPartition_;
  ProductType aOP_;
  CholeskyType bOP_;
  SolverType solverSmallest_;
  SolverType solverGreatest_;
 
  Eigen::VectorXd eigenvalues_;
  Eigen::MatrixXd eigenvectors_;
 
  bool computed_{};
 
  void assertCompute() const {
    if (not computed_)
      DUNE_THROW(Dune::InvalidStateException, "Eigenvalues and -vectors not yet computed, please call compute() first");
  }
};
 
template <Concepts::EigenMatrix MT>
struct GeneralizedSymEigenSolver<EigenValueSolverType::Eigen, MT>
{
  using ScalarType = typename MT::Scalar;
  using MatrixType = MT;
  using SolverType = Eigen::GeneralizedSelfAdjointEigenSolver<MatrixType>;
 
  template <typename MATA, typename MATB>
  requires(Concepts::EigenMatrix<std::remove_cvref_t<MATA>>)
  GeneralizedSymEigenSolver(MATA&& A, MATB&& B)
      : matA_(std::forward<MATA>(A)),
        matB_(std::forward<MATB>(B)),
        solver_(A.size()) {
    if ((A.cols() != B.cols()) or (A.rows() != B.rows()) or (A.cols() != A.rows()))
      DUNE_THROW(Dune::InvalidStateException,
                 "GeneralizedSymEigenSolver: The passed matrices should have the same size");
  }
 
  template <Concepts::FlatAssembler AssemblerA, Concepts::FlatAssembler AssemblerB>
  GeneralizedSymEigenSolver(std::shared_ptr<AssemblerA> assemblerA, std::shared_ptr<AssemblerB> assemblerB)
      : GeneralizedSymEigenSolver(assemblerA->matrix(), assemblerB->matrix()) {}
 
  bool compute(int options = Eigen::ComputeEigenvectors) {
    solver_.compute(matA_, matB_, options);
    computed_ = true; // SelfAdjointEigenSolver will always be successful if prerequisites are met
    return computed_;
  }
 
  auto& eigenvalues() const {
    assertCompute();
    return solver_.eigenvalues();
  }
 
  auto& eigenvectors() const {
    assertCompute();
    return solver_.eigenvectors();
  }
 
  Eigen::Index nev() const { return matA_.rows(); }
 
private:
  MatrixType matA_;
  MatrixType matB_;
  SolverType solver_;
  bool computed_{};
 
  void assertCompute() const {
    if (not computed_)
      DUNE_THROW(Dune::InvalidStateException, "Eigenvalues and -vectors not yet computed, please call compute() first");
  }
};
 
namespace Impl {
  void assertNev(std::optional<Eigen::Index> nev, Eigen::Index nevMax) {
    if (nev.has_value() and nev.value() > nevMax)
      DUNE_THROW(Dune::InvalidStateException, "You can not ask for more eigenvalues or eigenvectors than calculated.");
  }
} // namespace Impl
 
template <Concepts::DenseOrSparseEigenMatrix MT>
struct PartialGeneralizedSymEigenSolver
{
  using ScalarType              = typename MT::Scalar;
  static constexpr bool isDense = Concepts::EigenMatrix<MT>;
 
  using MatrixType = MT;
  using ProductType =
      std::conditional_t<isDense, Spectra::DenseSymMatProd<ScalarType>, Spectra::SparseSymMatProd<ScalarType>>;
  using CholeskyType =
      std::conditional_t<isDense, Spectra::DenseCholesky<ScalarType>, Spectra::SparseCholesky<ScalarType>>;
 
  using SolverType = Spectra::SymGEigsSolver<ProductType, CholeskyType, Spectra::GEigsMode::Cholesky>;
 
  template <typename MATA, typename MATB>
  requires(Concepts::DenseOrSparseEigenMatrix<std::remove_cvref_t<MATA>>)
  PartialGeneralizedSymEigenSolver(MATA&& A, MATB&& B, Eigen::Index nev)
      : nev_(nev),
        aOP_(std::forward<MATA>(A)),
        bOP_(std::forward<MATB>(B)),
        solver_(aOP_, bOP_, nev, 2 * nev <= A.cols() ? 2 * nev : A.cols()) {
    if ((A.cols() != B.cols()) or (A.rows() != B.rows()) or (A.cols() != A.rows()))
      DUNE_THROW(Dune::InvalidStateException,
                 "PartialGeneralizedSymEigenSolver: The passed matrices should have the same size");
  }
 
  template <Concepts::FlatAssembler AssemblerA, Concepts::FlatAssembler AssemblerB>
  PartialGeneralizedSymEigenSolver(std::shared_ptr<AssemblerA> assemblerA, std::shared_ptr<AssemblerB> assemblerB,
                                   Eigen::Index nev)
      : PartialGeneralizedSymEigenSolver(assemblerA->matrix(), assemblerB->matrix(), nev) {}
 
  bool compute(Spectra::SortRule selection = Spectra::SortRule::SmallestAlge,
               Spectra::SortRule sortRule = Spectra::SortRule::SmallestAlge, Eigen::Index maxit = 1000,
               ScalarType tolerance = 1e-10) {
    solver_.init();
    solver_.compute(selection, maxit, tolerance, sortRule);
 
    computed_ = solver_.info() == Spectra::CompInfo::Successful;
    return computed_;
  }
 
  Eigen::VectorXd eigenvalues() const {
    assertCompute();
    return solver_.eigenvalues();
  }
 
  auto eigenvectors(std::optional<Eigen::Index> _nev = std::nullopt) const {
    Impl::assertNev(_nev, nev_);
    assertCompute();
    return solver_.eigenvectors(_nev.value_or(nev_));
  }
 
  Eigen::Index nev() const { return nev_; }
 
private:
  Eigen::Index nev_;
  ProductType aOP_;
  CholeskyType bOP_;
  SolverType solver_;
 
  bool computed_{};
 
  void assertCompute() const {
    if (not computed_)
      DUNE_THROW(Dune::InvalidStateException, "Eigenvalues and -vectors not yet computed, please call compute() first");
  }
};
 
#ifndef DOXYGEN
template <Concepts::FlatAssembler AssemblerA, Concepts::FlatAssembler AssemblerB>
requires(std::same_as<typename AssemblerA::MatrixType, typename AssemblerB::MatrixType>)
PartialGeneralizedSymEigenSolver(std::shared_ptr<AssemblerA> as1, std::shared_ptr<AssemblerB> as2,
                                 Eigen::Index nev) -> PartialGeneralizedSymEigenSolver<typename AssemblerA::MatrixType>;
 
template <typename MATA, typename MATB>
requires(Concepts::DenseOrSparseEigenMatrix<std::remove_cvref_t<MATA>> &&
         std::same_as<std::remove_cvref_t<MATA>, std::remove_cvref_t<MATB>>)
PartialGeneralizedSymEigenSolver(MATA&&, MATB&&,
                                 Eigen::Index) -> PartialGeneralizedSymEigenSolver<std::remove_cvref_t<MATA>>;
 
#endif
} // namespace Ikarus