Homec4science
Diffusion SpecMiCP / ReactMiCP (reactmicp/unsaturated) Browse

specmicp/src/dfpmsolver/sparse_solver.hpp
849e20a62100reactmicp/unsaturated

Last Change
Raw File

sparse_solver.hpp
View Options

#ifndef SPECMICP_DFPMSOLVER_SPARSESOLVER_HPP
#define SPECMICP_DFPMSOLVER_SPARSESOLVER_HPP
#include <Eigen/Sparse>
#include "sparse_solver_structs.hpp"
namespace specmicp {
namespace dfpmsolver {
// header
// ======
//! \brief Solve a sparse linear system
template <class DerivedV1, class DerivedV2, class DerivedM>
int 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>
int 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>
int 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>
int solve_linear_system_bicgstab(DerivedV1& residual, DerivedM& jacobian, DerivedV2& update);
// Implementation
// ==============
template <class DerivedV1, class DerivedV2, class DerivedM>
int solve_sparse_linear_system(DerivedV1& residual, DerivedM& jacobian, DerivedV2& update, SparseSolver solver)
{
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>
int solve_linear_system_sparseqr(DerivedV1& residual, DerivedM& jacobian, DerivedV2& update)
{
Eigen::SparseQR<DerivedM, Eigen::COLAMDOrdering<int>> solver(jacobian);
if (solver.info() != Eigen::Success)
{
return -2;
}
update = solver.solve(-residual);
if (solver.info() != Eigen::Success)
{
return -1;
}
return 1;
}
template <class DerivedV1, class DerivedV2, class DerivedM>
int solve_linear_system_sparselu(DerivedV1& residual, DerivedM& jacobian, DerivedV2& update)
{
Eigen::SparseLU<DerivedM> solver(jacobian);
if (solver.info() != Eigen::Success)
{
return -2;
}
update = solver.solve(-residual);
if (solver.info() != Eigen::Success)
{
return -1;
}
return 1;
}
template <class DerivedV1, class DerivedV2, class DerivedM>
int 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 -1;
}
return 1;
}
} // end namespace dfpmsolver
} // end namespace specmicp
#endif //SPECMICP_DFPMSOLVER_SPARSESOLVER_HPP
c4science · Help