Page MenuHomec4science
Files F65772510

sparse_gmres.hpp
No OneTemporary
Actions

File Metadata

Created
Thu, Jun 6, 03:13

sparse_gmres.hpp
View Options

#ifndef SPECMICP_SPARSESOLVERS_SPARSEGMRES_HPP
#define SPECMICP_SPARSESOLVERS_SPARSEGMRES_HPP
#include "sparse_solver_structs.hpp"
#include "sparse_solver_base.hpp"
// needed for eigen 3.2.2
#include <iostream>
#include <Eigen/IterativeSolvers>
namespace specmicp {
namespace sparse_solvers {
template <class MatrixT, class DerivedR, class DerivedS>
class SparseSolverGMRES: public SparseSolverBase<MatrixT, DerivedR, DerivedS>
{
using SolverT = Eigen::GMRES<MatrixT, Eigen::IncompleteLUT<typename MatrixT::Scalar>> ;
public:
void analyse_pattern(MatrixT& jacobian)
{
}
//! \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_SPARSEGMRES_HPP

Event Timeline

c4science · Help