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/solver.h"
#include "amici/sundials_matrix_wrapper.h"
#include "amici/vector.h"

#include <memory>

namespace amici {

class Model;
class Solver;
class AmiVector;
struct SimulationState;

class NewtonSolver {

  public:
    explicit NewtonSolver(Model const& model);

    static std::unique_ptr<NewtonSolver>
    getSolver(Solver const& simulationSolver, Model const& model);

    void getStep(AmiVector& delta, Model& model, SimulationState const& state);

    void computeNewtonSensis(
        AmiVectorArray& sx, Model& model, SimulationState const& state
    );

    virtual void prepareLinearSystem(Model& model, SimulationState const& state)
        = 0;

    virtual void
    prepareLinearSystemB(Model& model, SimulationState const& state)
        = 0;

    virtual void solveLinearSystem(AmiVector& rhs) = 0;

    virtual void reinitialize() = 0;

    virtual bool is_singular(Model& model, SimulationState const& state) const
        = 0;

    virtual ~NewtonSolver() = default;

  protected:
    AmiVector xdot_;
    AmiVector x_;
    AmiVector xB_;
    AmiVector dxB_;
};

class NewtonSolverDense : public NewtonSolver {

  public:
    explicit NewtonSolverDense(Model const& model);

    NewtonSolverDense(NewtonSolverDense const&) = delete;

    NewtonSolverDense& operator=(NewtonSolverDense const& other) = delete;

    ~NewtonSolverDense() override;

    void solveLinearSystem(AmiVector& rhs) override;

    void
    prepareLinearSystem(Model& model, SimulationState const& state) override;

    void
    prepareLinearSystemB(Model& model, SimulationState const& state) override;

    void reinitialize() override;

    bool is_singular(Model& model, SimulationState const& state) const override;

  private:
    SUNMatrixWrapper Jtmp_;

    SUNLinearSolver linsol_{nullptr};
};

class NewtonSolverSparse : public NewtonSolver {

  public:
    explicit NewtonSolverSparse(Model const& model);

    NewtonSolverSparse(NewtonSolverSparse const&) = delete;

    NewtonSolverSparse& operator=(NewtonSolverSparse const& other) = delete;

    ~NewtonSolverSparse() override;

    void solveLinearSystem(AmiVector& rhs) override;

    void
    prepareLinearSystem(Model& model, SimulationState const& state) override;

    void
    prepareLinearSystemB(Model& model, SimulationState const& state) override;

    bool is_singular(Model& model, SimulationState const& state) const override;

    void reinitialize() override;

  private:
    SUNMatrixWrapper Jtmp_;

    SUNLinearSolver linsol_{nullptr};
};

} // namespace amici

#endif // NEWTON_SOLVER