EigenvalueSolver.hpp

#ifndef EIGENVALUE_SOLVER_HPP
#define EIGENVALUE_SOLVER_HPP

#include <algebra/PETScUtils.hpp>

#include <algebra/SmallMatrix.hpp>
#include <algebra/SmallVector.hpp>
#include <analysis/Polynomial.hpp>
#include <utils/Exceptions.hpp>
#include <utils/SmallArray.hpp>

class EigenvalueSolver
{
 private:
  struct Internals;
  TinyMatrix<3> eigenvectors0{1, 0, 0, 0, 1, 0, 0, 0, 1};
  TinyVector<3> eigenvalues0{0, 0, 0};

#ifdef PUGS_HAS_SLEPC
  void computeForSymmetricMatrix(const PETScAijMatrixEmbedder& A, SmallArray<double>& eigenvalues);

  void computeForSymmetricMatrix(const PETScAijMatrixEmbedder& A,
                                 SmallArray<double>& eigenvalues,
                                 std::vector<SmallVector<double>>& eigenvectors);

  void computeForSymmetricMatrix(const PETScAijMatrixEmbedder& A,
                                 SmallArray<double>& eigenvalues,
                                 SmallMatrix<double>& P);
  void computeExpForSymmetricMatrix(const PETScAijMatrixEmbedder& A, SmallMatrix<double>& expA);
#endif   // PUGS_HAS_SLEPC

 public:
  template <typename MatrixType>
  void
  computeForSymmetricMatrix([[maybe_unused]] const MatrixType& A, [[maybe_unused]] SmallArray<double>& eigenvalues)
  {
#ifdef PUGS_HAS_SLEPC
    this->computeForSymmetricMatrix(PETScAijMatrixEmbedder{A}, eigenvalues);
#else    // PUGS_HAS_SLEPC
    throw NotImplementedError("SLEPc is required to solve eigenvalue problems");
#endif   // PUGS_HAS_SLEPC
  }

  template <typename MatrixType>
  void
  computeForSymmetricMatrix([[maybe_unused]] const MatrixType& A,
                            [[maybe_unused]] SmallArray<double>& eigenvalues,
                            [[maybe_unused]] std::vector<SmallVector<double>>& eigenvectors)
  {
#ifdef PUGS_HAS_SLEPC
    this->computeForSymmetricMatrix(PETScAijMatrixEmbedder{A}, eigenvalues, eigenvectors);
#else    // PUGS_HAS_SLEPC
    throw NotImplementedError("SLEPc is required to solve eigenvalue problems");
#endif   // PUGS_HAS_SLEPC
  }

  template <typename MatrixType>
  void
  computeForSymmetricMatrix([[maybe_unused]] const MatrixType& A,
                            [[maybe_unused]] SmallArray<double>& eigenvalues,
                            [[maybe_unused]] SmallMatrix<double>& P)
  {
#ifdef PUGS_HAS_SLEPC
    this->computeForSymmetricMatrix(PETScAijMatrixEmbedder{A}, eigenvalues, P);
#else    // PUGS_HAS_SLEPC
    throw NotImplementedError("SLEPc is required to solve eigenvalue problems");
#endif   // PUGS_HAS_SLEPC
  }

  template <typename MatrixType>
  void
  computeExpForSymmetricMatrix([[maybe_unused]] const MatrixType& A, [[maybe_unused]] SmallMatrix<double>& expA)
  {
#ifdef PUGS_HAS_SLEPC
    this->computeExpForSymmetricMatrix(PETScAijMatrixEmbedder{A}, expA);
#else    // PUGS_HAS_SLEPC
    throw NotImplementedError("SLEPc is required to solve eigenvalue problems");
#endif   // PUGS_HAS_SLEPC
  }

  template <typename T>
  PUGS_INLINE TinyMatrix<3, 3, T> swap(TinyMatrix<3, 3, T>& matrix, size_t i, size_t j) const;

  template <typename T>
  PUGS_INLINE constexpr TinyMatrix<3, 3, T> rowReduce(const TinyMatrix<3, 3, T>& matrix, const double epsilon) const;

  TinyMatrix<3, 3> findEigenvector(const TinyMatrix<3, 3>& A,
                                   const double eigenvalue,
                                   const size_t space_size,
                                   const double epsilon) const;

  bool isDiagonal(const TinyMatrix<3, 3>& A, const double epsilon);

  TinyMatrix<3, 3> findEigenvectors(const TinyMatrix<3, 3>& A, TinyVector<3> eigenvalues, const double epsilon) const;

  std::tuple<TinyVector<3>, TinyMatrix<3>> findEigen(const TinyMatrix<3> A);
  TinyVector<3> _findEigenValues(const TinyMatrix<3>& A, const double epsilon) const;

  double _findFirstRoot(Polynomial<3> P) const;

  TinyVector<2> _findLastRoots(Polynomial<2> P, const double epsilon) const;

  TinyMatrix<3> computeExpForTinyMatrix(const TinyMatrix<3>& A);

  EigenvalueSolver();
  ~EigenvalueSolver() = default;
};

#endif   // EIGENVALUE_SOLVER_HPP