sparse_lu.hpp
14f16a9d1764reactivetransport

sparse_lu.hpp
View Options

			#ifndef SPECMICP_SPARSESOLVERS_SPARSELU_HPP
			#define SPECMICP_SPARSESOLVERS_SPARSELU_HPP

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

			#include <Eigen/SparseLU>

			namespace specmicp {
			namespace sparse_solvers {

			template <class MatrixT, class DerivedR, class DerivedS>
			class SparseSolverLU: public SparseSolverBase<MatrixT, DerivedR, DerivedS>
			{
			using SolverT = Eigen::SparseLU<MatrixT, Eigen::COLAMDOrdering<int>>;

			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_SPARSELU_HPP