#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 <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);

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

  template <typename MatrixType>
  void
  computeForSymmetricMatrix(const MatrixType& A, SmallArray<double>& eigenvalues)
  {
    Assert((eigenvalues.size() - A.numberOfRows() == 0) and A.isSquare());

    switch (m_options.library()) {
    case ESLibrary::eigen3: {
      throw NotImplementedError("Eigen3");
    }
    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((eigenvalues.size() - A.numberOfRows() == 0) and (A.numberOfRows() - eigenvectors.size() == 0) and
           A.isSquare());

    switch (m_options.library()) {
    case ESLibrary::eigen3: {
      throw NotImplementedError("Eigen3");
    }
    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((eigenvalues.size() - A.numberOfRows() == 0) and (A.numberOfRows() - P.numberOfRows() == 0) and
           A.isSquare() and P.isSquare());

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

  EigenvalueSolver(const EigenvalueSolverOptions& options = EigenvalueSolverOptions::default_options);
  ~EigenvalueSolver() = default;
};

#endif   // EIGENVALUE_SOLVER_HPP