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