Homec4science
Diffusion SpecMiCP / ReactMiCP (fixpoint) Browse

specmicp/src/dfpmsolver/sparse_solver.hpp
eebb86eba7d7fixpoint

Last Change
Raw File

sparse_solver.hpp
View Options

#ifndef SPECMICP_DFPMSOLVER_SPARSESOLVER_HPP
#define SPECMICP_DFPMSOLVER_SPARSESOLVER_HPP
#include <Eigen/Dense>
#include <Eigen/Sparse>
#include "sparse_solver_structs.hpp"
namespace specmicp {
namespace dfpmsolver {
// header
// ======
//! \brief Solve a sparse linear system
template <class Vector, class Matrix>
int solve_sparse_linear_system(Vector& residual, Matrix& jacobian, Vector& update, SparseSolver solver);
//! \brief Solve a sparse linear system using the QR decomposition
template <class Vector, class Matrix>
int solve_linear_system_sparseqr(Vector& residual, Matrix& jacobian, Vector& update);
//! \brief Solve a sparse linear system using the LU decomposition
template <class Vector, class Matrix>
int solve_linear_system_sparselu(Vector& residual, Matrix& jacobian, Vector& update);
//! \brief Solve a sparse linear system using the BiCGSTAB algorithm
template <class Vector, class Matrix>
int solve_linear_system_bicgstab(Vector& residual, Matrix& jacobian, Vector& update);
// Implementation
// ==============
template <class Vector, class Matrix>
int solve_sparse_linear_system(Vector& residual, Matrix& jacobian, Vector& 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 Vector, class Matrix>
int solve_linear_system_sparseqr(Vector& residual, Matrix& jacobian, Vector& update)
{
Eigen::SparseQR<Matrix, 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 Vector, class Matrix>
int solve_linear_system_sparselu(Vector& residual, Matrix& jacobian, Vector& update)
{
Eigen::SparseLU<Matrix> solver(jacobian);
if (solver.info() != Eigen::Success)
{
return -2;
}
update = solver.solve(-residual);
if (solver.info() != Eigen::Success)
{
return -1;
}
return 1;
}
template <class Vector, class Matrix>
int solve_linear_system_bicgstab(Vector& residual, Matrix& jacobian, Vector& update)
{
Eigen::BiCGSTAB<Matrix, Eigen::IncompleteLUT<typename Matrix::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