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

#include <memory>

namespace amici {

class Model;
class Solver;
class AmiVector;

class NewtonSolver {

  public:
    explicit NewtonSolver(const Model &model);

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

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

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

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

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

    virtual void solveLinearSystem(AmiVector &rhs) = 0;

    virtual void reinitialize() = 0;

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

    virtual ~NewtonSolver() = default;

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

class NewtonSolverDense : public NewtonSolver {

  public:
    explicit NewtonSolverDense(const Model &model);

    NewtonSolverDense(const NewtonSolverDense &) = delete;

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

    ~NewtonSolverDense() override;

    void solveLinearSystem(AmiVector &rhs) override;

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

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

    void reinitialize() override;

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

  private:
    SUNMatrixWrapper Jtmp_;

    SUNLinearSolver linsol_{nullptr};
};

class NewtonSolverSparse : public NewtonSolver {

  public:
    explicit NewtonSolverSparse(const Model &model);

    NewtonSolverSparse(const NewtonSolverSparse &) = delete;

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

    ~NewtonSolverSparse() override;

    void solveLinearSystem(AmiVector &rhs) override;

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

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

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

    void reinitialize() override;

  private:
    SUNMatrixWrapper Jtmp_;

    SUNLinearSolver linsol_{nullptr};
};

} // namespace amici

#endif // NEWTON_SOLVER