Homec4science
Diffusion SpecMiCP / ReactMiCP (cython) Browse

specmicp/src/utils/sparse_solver.hpp
29fe36a29e5acython

Last Change
Raw File

sparse_solver.hpp
View Options

#ifndef SPECMICP_UTILS_SPARSESOLVER_HPP
#define SPECMICP_UTILS_SPARSESOLVER_HPP
#include <Eigen/Sparse>
#include "sparse_solver_structs.hpp"
namespace specmicp {
// header
// ======
//! \brief Solve a sparse linear system
template <class DerivedV1, class DerivedV2, class DerivedM>
SparseSolverReturnCode solve_sparse_linear_system(
DerivedV1& residual,
DerivedM& jacobian,
DerivedV2& update,
SparseSolver solver
);
//! \brief Solve a sparse linear system using the QR decomposition
template <class DerivedV1, class DerivedV2, class DerivedM>
SparseSolverReturnCode solve_linear_system_sparseqr(
DerivedV1& residual,
DerivedM& jacobian,
DerivedV2& update
);
//! \brief Solve a sparse linear system using the LU decomposition
template <class DerivedV1, class DerivedV2, class DerivedM>
SparseSolverReturnCode solve_linear_system_sparselu(
DerivedV1& residual,
DerivedM& jacobian,
DerivedV2& update
);
//! \brief Solve a sparse linear system using the BiCGSTAB algorithm
template <class DerivedV1, class DerivedV2, class DerivedM>
SparseSolverReturnCode solve_linear_system_bicgstab(
DerivedV1& residual,
DerivedM& jacobian,
DerivedV2& update
);
//! \brief A sparse system solver
template <class SolverType>
class SparseSystemSolver
{
public:
//! \brief Type of the matrix;
using MatrixType = typename SolverType::MatrixType;
//! \brief Analyse the sparsity matrix of the jacobian
void analyse_sparsity() {
m_solver.analysePattern(m_jacobian);
m_structure_is_known = true;
}
//! \brief Compute the decomposition of the jacobian (if available)
SparseSolverReturnCode compute();
//! \brief Solve the problem, using 'rhs' and storing solution in 'solution'
template <class RHSType, class UpdateType>
SparseSolverReturnCode solve(const RHSType& rhs, UpdateType& solution);
//! \brief Return a reference to the jacobian, so it can be updated
MatrixType& jacobian() {
m_jacobian_has_changed=true;
return m_jacobian;
}
//! \brief Const reference to the jacobian
const MatrixType& jacobian() const {return m_jacobian;}
private:
bool m_structure_is_known; //!< True if the structure of the jacobian is known
bool m_jacobian_has_changed; //!< True if the jacobian has changed
MatrixType m_jacobian; //!< The jacobian
SolverType m_solver; //!< The solver
};
// Implementation
// ==============
template <class SolverType>
SparseSolverReturnCode SparseSystemSolver<SolverType>::compute()
{
// first compress the matrix to avoid copy
if (not m_jacobian.isCompressed()) m_jacobian.makeCompressed();
if (m_structure_is_known)
m_solver.factorize(m_jacobian);
else
m_solver.compute(m_jacobian);
if (m_solver.info() != Eigen::Success)
return SparseSolverReturnCode::FailedDecomposition;
return SparseSolverReturnCode::Success;
}
template <class SolverType>
template <class RHSType, class UpdateType>
SparseSolverReturnCode SparseSystemSolver<SolverType>::solve(const RHSType& rhs, UpdateType& solution)
{
if (m_jacobian_has_changed) compute();
solution = m_solver.solve(rhs);
if (m_solver.info() != Eigen::Success)
return SparseSolverReturnCode::FailedSystemSolving;
return SparseSolverReturnCode::Success;
}
// Implementation
// ==============
template <class DerivedV1, class DerivedV2, class DerivedM>
SparseSolverReturnCode solve_sparse_linear_system(
DerivedV1& residual,
DerivedM& jacobian,
DerivedV2& update,
SparseSolver solver)
{
if (not jacobian.isCompressed()) jacobian.makeCompressed();
if (solver == SparseSolver::SparseQR) {
return solve_linear_system_sparseqr(residual, jacobian, update);
} else if (solver == SparseSolver::SparseLU) {
return solve_linear_system_sparselu(residual, jacobian, update);
}
else {
return solve_linear_system_bicgstab(residual, jacobian, update);
}
}
template <class DerivedV1, class DerivedV2, class DerivedM>
SparseSolverReturnCode solve_linear_system_sparseqr(DerivedV1& residual, DerivedM& jacobian, DerivedV2& update)
{
Eigen::SparseQR<DerivedM, Eigen::COLAMDOrdering<int>> solver(jacobian);
if (solver.info() != Eigen::Success)
{
return SparseSolverReturnCode::FailedDecomposition;
}
update = solver.solve(-residual);
if (solver.info() != Eigen::Success)
{
return SparseSolverReturnCode::FailedSystemSolving;
}
return SparseSolverReturnCode::Success;
}
template <class DerivedV1, class DerivedV2, class DerivedM>
SparseSolverReturnCode solve_linear_system_sparselu(DerivedV1& residual, DerivedM& jacobian, DerivedV2& update)
{
Eigen::SparseLU<DerivedM> solver(jacobian);
if (solver.info() != Eigen::Success)
{
return SparseSolverReturnCode::FailedDecomposition;;
}
update = solver.solve(-residual);
if (solver.info() != Eigen::Success)
{
return SparseSolverReturnCode::FailedSystemSolving;
}
return SparseSolverReturnCode::Success;
}
template <class DerivedV1, class DerivedV2, class DerivedM>
SparseSolverReturnCode solve_linear_system_bicgstab(DerivedV1& residual, DerivedM& jacobian, DerivedV2& update)
{
Eigen::BiCGSTAB<DerivedM, Eigen::IncompleteLUT<typename DerivedM::Scalar>> solver(jacobian);
update = solver.solve(-residual);
if (solver.info() != Eigen::Success)
{
return SparseSolverReturnCode::FailedSystemSolving;
}
return SparseSolverReturnCode::Success;
}
} // end namespace specmicp
#endif //SPECMICP_UTILS_SPARSESOLVER_HPP
c4science · Help