sparse_gmres.hpp
14f16a9d1764reactivetransport

sparse_gmres.hpp
View Options

			#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