Refactor NewtonSolver (#2521)

Use SUNLinSolWrapper. Get rid of subclasses. Simplify.

Closes #1164.

We should also be able to get rid of the remaining data members,
but that's for another time...

include/amici/newton_solver.h

@@ -2,11 +2,8 @@
 #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;
@@ -26,35 +23,30 @@ class NewtonSolver {
  * @brief Initializes solver according to the dimensions in the provided
  * model
  *
- * @param model pointer to the model object
+ * @param model the model object
+ * @param linsol_type type of linear solver to use
  */
- explicit NewtonSolver(Model const& model);
+ explicit NewtonSolver(Model const& model, LinearSolver linsol_type);
 
- /**
- * @brief Factory method to create a NewtonSolver based on linsolType
- *
- * @param simulationSolver solver with settings
- * @param model pointer to the model instance
- * @return solver NewtonSolver according to the specified linsolType
- */
- static std::unique_ptr<NewtonSolver>
- getSolver(Solver const& simulationSolver, Model const& model);
+ NewtonSolver(NewtonSolver const&) = delete;
+
+ NewtonSolver& operator=(NewtonSolver const& other) = delete;
 
  /**
  * @brief Computes the solution of one Newton iteration
  *
  * @param delta containing the RHS of the linear system, will be
  * overwritten by solution to the linear system
- * @param model pointer to the model instance
+ * @param model the model instance
  * @param state current simulation state
  */
  void getStep(AmiVector& delta, Model& model, SimulationState const& state);
 
  /**
  * @brief Computes steady state sensitivities
  *
- * @param sx pointer to state variable sensitivities
- * @param model pointer to the model instance
+ * @param sx state variable sensitivities
+ * @param model the model instance
  * @param state current simulation state
  */
  void computeNewtonSensis(
@@ -65,51 +57,45 @@ class NewtonSolver {
  * @brief Writes the Jacobian for the Newton iteration and passes it to the
  * linear solver
  *
- * @param model pointer to the model instance
+ * @param model the model instance
  * @param state current simulation state
  */
- virtual void prepareLinearSystem(Model& model, SimulationState const& state)
- = 0;
+ void prepareLinearSystem(Model& model, SimulationState const& state);
 
  /**
  * Writes the Jacobian (JB) for the Newton iteration and passes it to the
  * linear solver
  *
- * @param model pointer to the model instance
+ * @param model the model instance
  * @param state current simulation state
  */
- virtual void
- prepareLinearSystemB(Model& model, SimulationState const& state)
- = 0;
+ void prepareLinearSystemB(Model& model, SimulationState const& state);
 
  /**
  * @brief Solves the linear system for the Newton step
  *
  * @param rhs containing the RHS of the linear system, will be
  * overwritten by solution to the linear system
  */
- virtual void solveLinearSystem(AmiVector& rhs) = 0;
+ void solveLinearSystem(AmiVector& rhs);
 
  /**
  * @brief Reinitialize the linear solver
  *
  */
- virtual void reinitialize() = 0;
+ void reinitialize();
 
  /**
- * @brief Checks whether linear system is singular
+ * @brief Checks whether the linear system is singular
  *
- * @param model pointer to the model instance
+ * @param model the model instance
  * @param state current simulation state
  * @return boolean indicating whether the linear system is singular
  * (condition number < 1/machine precision)
  */
- virtual bool is_singular(Model& model, SimulationState const& state) const
- = 0;
-
- virtual ~NewtonSolver() = default;
+ bool is_singular(Model& model, SimulationState const& state) const;
 
- protected:
+ private:
  /** dummy rhs, used as dummy argument when computing J and JB */
  AmiVector xdot_;
  /** dummy state, attached to linear solver */
@@ -119,88 +105,9 @@ class NewtonSolver {
  /** dummy differential adjoint state, used as dummy argument when computing
  * JB */
  AmiVector dxB_;
-};
-
-/**
- * @brief The NewtonSolverDense provides access to the dense linear solver for
- * the Newton method.
- */
-
-class NewtonSolverDense : public NewtonSolver {
-
- public:
- /**
- * @brief constructor for sparse solver
- *
- * @param model model instance that provides problem dimensions
- */
- 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:
- /** temporary storage of Jacobian */
- SUNMatrixWrapper Jtmp_;
-
- /** dense linear solver */
- SUNLinearSolver linsol_{nullptr};
-};
-
-/**
- * @brief The NewtonSolverSparse provides access to the sparse linear solver for
- * the Newton method.
- */
-
-class NewtonSolverSparse : public NewtonSolver {
-
- public:
- /**
- * @brief constructor for dense solver
- *
- * @param model model instance that provides problem dimensions
- */
- 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:
- /** temporary storage of Jacobian */
- SUNMatrixWrapper Jtmp_;
 
- /** sparse linear solver */
- SUNLinearSolver linsol_{nullptr};
+ /** linear solver */
+ std::unique_ptr<SUNLinSolWrapper> linsol_;
 };
 
 } // namespace amici

include/amici/steadystateproblem.h

@@ -452,7 +452,7 @@ class SteadystateProblem {
  realtype rtol_quad_{NAN};
 
  /** newton solver */
- std::unique_ptr<NewtonSolver> newton_solver_{nullptr};
+ NewtonSolver newton_solver_;
 
  /** damping factor flag */
  NewtonDampingFactorMode damping_factor_mode_{NewtonDampingFactorMode::on};

include/amici/sundials_linsol_wrapper.h

@@ -195,7 +195,8 @@ class SUNLinSolKLU : public SUNLinSolWrapper {
  * @param ordering
  */
  SUNLinSolKLU(
- AmiVector const& x, int nnz, int sparsetype, StateOrdering ordering = StateOrdering::COLAMD
+ AmiVector const& x, int nnz, int sparsetype,
+ StateOrdering ordering = StateOrdering::COLAMD
  );
 
  /**
@@ -214,6 +215,14 @@ class SUNLinSolKLU : public SUNLinSolWrapper {
  * @param ordering
  */
  void setOrdering(StateOrdering ordering);
+
+ /**
+ * @brief Checks whether the linear system is singular
+ *
+ * @return boolean indicating whether the linear system is singular
+ * (condition number < 1/machine precision)
+ */
+ bool is_singular() const;
 };
 
 #ifdef SUNDIALS_SUPERLUMT