#ifndef EIGENVALUE_SOLVER_HPP
#define EIGENVALUE_SOLVER_HPP

#include <algebra/EigenvalueSolverOptions.hpp>

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

class EigenvalueSolver
{
 private:
  struct Internals;

  const EigenvalueSolverOptions m_options;

  void _slepscComputeForSymmetricMatrix(const CRSMatrix<double>& A, SmallArray<double>& eigenvalues);

  void _slepscComputeForSymmetricMatrix(const DenseMatrixWrapper<double>& A, SmallArray<double>& eigenvalues);

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

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

  void _slepscComputeForSymmetricMatrix(const CRSMatrix<double>& A,
                                        SmallArray<double>& eigenvalues,
                                        SmallMatrix<double>& P);

  void _slepscComputeForSymmetricMatrix(const DenseMatrixWrapper<double>& A,
                                        SmallArray<double>& eigenvalues,
                                        SmallMatrix<double>& P);

  void _eigen3ComputeForSymmetricMatrix(const CRSMatrix<double>& A, SmallArray<double>& eigenvalues);

  void _eigen3ComputeForSymmetricMatrix(const DenseMatrixWrapper<double>& A, SmallArray<double>& eigenvalues);

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

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

  void _eigen3ComputeForSymmetricMatrix(const CRSMatrix<double>& A,
                                        SmallArray<double>& eigenvalues,
                                        SmallMatrix<double>& P);

  void _eigen3ComputeForSymmetricMatrix(const DenseMatrixWrapper<double>& A,
                                        SmallArray<double>& eigenvalues,
                                        SmallMatrix<double>& P);

 public:
  static bool hasLibrary(ESLibrary library);
  static void checkHasLibrary(const ESLibrary library);

  template <typename MatrixType>
  void
  computeForSymmetricMatrix(const MatrixType& A, SmallArray<double>& eigenvalues)
  {
    Assert(A.isSquare());

    switch (m_options.library()) {
    case ESLibrary::eigen3: {
      this->_eigen3ComputeForSymmetricMatrix(A, eigenvalues);
      break;
    }
    case ESLibrary::slepsc: {
      this->_slepscComputeForSymmetricMatrix(A, eigenvalues);
      break;
    }
    default: {
      throw UnexpectedError(::name(m_options.library()) + " cannot compute for symmetric matrices");
    }
    }
  }

  template <typename MatrixType>
  void
  computeForSymmetricMatrix(const MatrixType& A,
                            SmallArray<double>& eigenvalues,
                            std::vector<SmallVector<double>>& eigenvectors)
  {
    Assert(A.isSquare());

    switch (m_options.library()) {
    case ESLibrary::eigen3: {
      this->_eigen3ComputeForSymmetricMatrix(A, eigenvalues, eigenvectors);
      break;
    }
    case ESLibrary::slepsc: {
      this->_slepscComputeForSymmetricMatrix(A, eigenvalues, eigenvectors);
      break;
    }
    default: {
      throw UnexpectedError(::name(m_options.library()) + " cannot compute for symmetric matrices");
    }
    }
  }

  template <typename MatrixType>
  void
  computeForSymmetricMatrix(const MatrixType& A, SmallArray<double>& eigenvalues, SmallMatrix<double>& P)
  {
    Assert(A.isSquare());

    switch (m_options.library()) {
    case ESLibrary::eigen3: {
      this->_eigen3ComputeForSymmetricMatrix(A, eigenvalues, P);
      break;
    }
    case ESLibrary::slepsc: {
      this->_slepscComputeForSymmetricMatrix(A, eigenvalues, P);
      break;
    }
    default: {
      throw UnexpectedError(::name(m_options.library()) + " cannot compute for symmetric matrices");
    }
    }
  }

  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(const EigenvalueSolverOptions& options = EigenvalueSolverOptions::default_options);
  ~EigenvalueSolver() = default;
};

#endif   // EIGENVALUE_SOLVER_HPP