Program Listing for File newton_solver.h¶
↰ Return to documentation for file (include/amici/newton_solver.h
)
#ifndef amici_newton_solver_h
#define amici_newton_solver_h
#include "amici/vector.h"
#include "amici/defines.h"
#include "amici/sundials_matrix_wrapper.h"
#include "amici/sundials_linsol_wrapper.h"
#include <memory>
namespace amici {
class Model;
class Solver;
class AmiVector;
class NewtonSolver {
public:
NewtonSolver(realtype *t, AmiVector *x, Model *model);
static std::unique_ptr<NewtonSolver> getSolver(
realtype *t, AmiVector *x, Solver &simulationSolver, Model *model);
void getStep(int ntry, int nnewt, AmiVector &delta);
void computeNewtonSensis(AmiVectorArray &sx);
const std::vector<int> &getNumLinSteps() const {
return num_lin_steps_;
}
virtual void prepareLinearSystem(int ntry, int nnewt) = 0;
virtual void prepareLinearSystemB(int ntry, int nnewt) = 0;
virtual void solveLinearSystem(AmiVector &rhs) = 0;
virtual ~NewtonSolver() = default;
int max_lin_steps_ {0};
int max_steps {0};
realtype atol_ {1e-16};
realtype rtol_ {1e-8};
NewtonDampingFactorMode damping_factor_mode_ {NewtonDampingFactorMode::on};
realtype damping_factor_lower_bound {1e-8};
protected:
realtype *t_;
Model *model_;
AmiVector xdot_;
AmiVector *x_;
AmiVector dx_;
std::vector<int> num_lin_steps_;
AmiVector xB_;
AmiVector dxB_;
};
class NewtonSolverDense : public NewtonSolver {
public:
NewtonSolverDense(realtype *t, AmiVector *x, Model *model);
NewtonSolverDense(const NewtonSolverDense&) = delete;
NewtonSolverDense& operator=(const NewtonSolverDense& other) = delete;
~NewtonSolverDense() override;
void solveLinearSystem(AmiVector &rhs) override;
void prepareLinearSystem(int ntry, int nnewt) override;
void prepareLinearSystemB(int ntry, int nnewt) override;
private:
SUNMatrixWrapper Jtmp_;
SUNLinearSolver linsol_ {nullptr};
};
class NewtonSolverSparse : public NewtonSolver {
public:
NewtonSolverSparse(realtype *t, AmiVector *x, Model *model);
NewtonSolverSparse(const NewtonSolverSparse&) = delete;
NewtonSolverSparse& operator=(const NewtonSolverSparse& other) = delete;
~NewtonSolverSparse() override;
void solveLinearSystem(AmiVector &rhs) override;
void prepareLinearSystem(int ntry, int nnewt) override;
void prepareLinearSystemB(int ntry, int nnewt) override;
private:
SUNMatrixWrapper Jtmp_;
SUNLinearSolver linsol_ {nullptr};
};
class NewtonSolverIterative : public NewtonSolver {
public:
NewtonSolverIterative(realtype *t, AmiVector *x, Model *model);
~NewtonSolverIterative() override = default;
void solveLinearSystem(AmiVector &rhs) override;
void prepareLinearSystem(int ntry, int nnewt) override;
void prepareLinearSystemB(int ntry, int nnewt) override;
void linsolveSPBCG(int ntry, int nnewt, AmiVector &ns_delta);
private:
int newton_try_ {0};
int i_newton_ {0};
AmiVector ns_p_;
AmiVector ns_h_;
AmiVector ns_t_;
AmiVector ns_s_;
AmiVector ns_r_;
AmiVector ns_rt_;
AmiVector ns_v_;
AmiVector ns_Jv_;
AmiVector ns_tmp_;
AmiVector ns_Jdiag_;
SUNMatrixWrapper ns_J_;
};
} // namespace amici
#endif // NEWTON_SOLVER