sparse_lu.hpp
No OneTemporary
Actions

Subscribers

None

File Metadata

	#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<typename MatrixT::Index>>;

	public:
	void analyse_pattern(MatrixT& jacobian)
	{
	m_solver.analyzePattern(jacobian);
	}

	//! \brief Decompose the jacboian
	SparseSolverReturnCode decompose(MatrixT& jacobian)
	{
	m_solver.factorize(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