sparse_gmres.hpp
No OneTemporary
Actions

Subscribers

None

File Metadata

	#ifndef SPECMICP_SPARSESOLVERS_SPARSEGMRES_HPP
	#define SPECMICP_SPARSESOLVERS_SPARSEGMRES_HPP

	#include "sparse_solver_structs.hpp"
	#include "sparse_solver_base.hpp"

	// needed for eigen 3.2.2
	#include <iostream>

	#include <Eigen/IterativeSolvers>

	namespace specmicp {
	namespace sparse_solvers {

	template <class MatrixT, class DerivedR, class DerivedS>
	class SparseSolverGMRES: public SparseSolverBase<MatrixT, DerivedR, DerivedS>
	{
	using SolverT = Eigen::GMRES<MatrixT, Eigen::IncompleteLUT<typename MatrixT::Scalar>> ;

	public:
	//! \brief Decompose the jacboian
	SparseSolverReturnCode decompose(MatrixT& jacobian)
	{
	m_solver.compute(jacobian);
	if (m_solver.info() != Eigen::Success)
	{
	return SparseSolverReturnCode::FailedDecomposition;
	}
	return SparseSolverReturnCode::Success;
	}

	//! \brief Solve the problem
	SparseSolverReturnCode solve(
	const DerivedR& residuals,
	DerivedS& solution
	)
	{
	solution = m_solver.solve(-residuals);
	if (m_solver.info() != Eigen::Success)
	{
	return SparseSolverReturnCode::FailedSystemSolving;
	}
	return SparseSolverReturnCode::Success;
	}

	SparseSolverReturnCode solve_scaling(
	const DerivedR& residuals,
	const DerivedR& scaling,
	DerivedS& solution
	)
	{
	solution = m_solver.solve(scaling.asDiagonal()*(-residuals));
	if (m_solver.info() != Eigen::Success)
	{
	return SparseSolverReturnCode::FailedSystemSolving;
	}
	return SparseSolverReturnCode::Success;
	}

	private:
	SolverT m_solver;
	};

	} // end namespace sparse_solvers
	} // end namespace specmicp

	#endif //SPECMICP_SPARSESOLVERS_SPARSEGMRES_HPP