remove SPBCG solver (#1729)

* remove SPBCG solver

* remove test

* fixup

* retry deleting from hdf5

* fixup remove options

* remove python test

include/amici/newton_solver.h

@@ -254,96 +254,6 @@ class NewtonSolverSparse : public NewtonSolver {
     SUNLinearSolver linsol_ {nullptr};
 };
 
-/**
- * @brief The NewtonSolverIterative provides access to the iterative linear
- * solver for the Newton method.
- */
-
-class NewtonSolverIterative : public NewtonSolver {
-
-  public:
-    /**
-     * @brief Constructor, initializes all members with the provided objects
-     * @param t pointer to time variable
-     * @param x pointer to state variables
-     * @param model pointer to the model object
-     */
-    NewtonSolverIterative(realtype *t, AmiVector *x, Model *model);
-
-    ~NewtonSolverIterative() override = default;
-
-    /**
-     * @brief Solves the linear system for the Newton step by passing it to
-     * linsolveSPBCG
-     *
-     * @param rhs containing the RHS of the linear system, will be
-     * overwritten by solution to the linear system
-     */
-    void solveLinearSystem(AmiVector &rhs) override;
-
-    /**
-     * Writes the Jacobian (J) for the Newton iteration and passes it to the linear
-     * solver.
-     * Also wraps around getSensis for iterative linear solver.
-     *
-     * @param ntry integer newton_try integer start number of Newton solver
-     * (1 or 2)
-     * @param nnewt integer number of current Newton step
-     */
-    void prepareLinearSystem(int ntry, int nnewt) override;
-
-    /**
-     * Writes the Jacobian (JB) for the Newton iteration and passes it to the linear
-     * solver.
-     * Also wraps around getSensis for iterative linear solver.
-     *
-     * @param ntry integer newton_try integer start number of Newton solver
-     * (1 or 2)
-     * @param nnewt integer number of current Newton step
-     */
-    void prepareLinearSystemB(int ntry, int nnewt) override;
-
-    /**
-     * Iterative linear solver created from SPILS BiCG-Stab.
-     * Solves the linear system within each Newton step if iterative solver is
-     * chosen.
-     *
-     * @param ntry integer newton_try integer start number of Newton solver
-     * (1 or 2)
-     * @param nnewt integer number of current Newton step
-     * @param ns_delta Newton step
-     */
-    void linsolveSPBCG(int ntry, int nnewt, AmiVector &ns_delta);
-
-  private:
-    /** number of tries  */
-    int newton_try_ {0};
-    /** number of iterations  */
-    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_;
-    /** temporary storage of Jacobian */
-    SUNMatrixWrapper ns_J_;
-};
-
 
 } // namespace amici
 

python/tests/test_preequilibration.py

@@ -359,19 +359,3 @@ def test_newton_solver_equilibration(preeq_fixture):
             rdatas[settings[1]][variable],
             1e-6, 1e-6
         ).all(), variable
-
-    # test failure for iterative linear solver with sensitivities
-    edata.fixedParametersPreequilibration = ()
-    edata.t_presim = 0.0
-    edata.fixedParametersPresimulation = ()
-
-    solver.setLinearSolver(amici.LinearSolver.SPBCG)
-    solver.setSensitivityMethod(amici.SensitivityMethod.adjoint)
-    solver.setSensitivityOrder(amici.SensitivityOrder.first)
-    model.setSteadyStateSensitivityMode(
-        amici.SteadyStateSensitivityMode.newtonOnly)
-    solver.setNewtonMaxSteps(10)
-    solver.setNewtonMaxLinearSteps(10)
-    rdata_spbcg = amici.runAmiciSimulation(model, solver, edata)
-
-    assert rdata_spbcg['status'] == amici.AMICI_ERROR

src/newton_solver.cpp

@@ -49,8 +49,7 @@ std::unique_ptr<NewtonSolver> NewtonSolver::getSolver(realtype *t, AmiVector *x,
         throw NewtonFailure(AMICI_NOT_IMPLEMENTED, "getSolver");
 
     case LinearSolver::SPBCG:
-        solver.reset(new NewtonSolverIterative(t, x, model));
-        break;
+        throw NewtonFailure(AMICI_NOT_IMPLEMENTED, "getSolver");
 
     case LinearSolver::SPTFQMR:
         throw NewtonFailure(AMICI_NOT_IMPLEMENTED, "getSolver");
@@ -71,8 +70,6 @@ std::unique_ptr<NewtonSolver> NewtonSolver::getSolver(realtype *t, AmiVector *x,
     solver->damping_factor_mode_ = simulationSolver.getNewtonDampingFactorMode();
     solver->damping_factor_lower_bound =
         simulationSolver.getNewtonDampingFactorLowerBound();
-    if (simulationSolver.getLinearSolver() == LinearSolver::SPBCG)
-        solver->num_lin_steps_.resize(simulationSolver.getNewtonMaxSteps(), 0);
 
     return solver;
 }
@@ -220,152 +217,5 @@ NewtonSolverSparse::~NewtonSolverSparse() {
         SUNLinSolFree_KLU(linsol_);
 }
 
-/* ------------------------------------------------------------------------- */
-/* - Iterative linear solver------------------------------------------------ */
-/* ------------------------------------------------------------------------- */
-
-NewtonSolverIterative::NewtonSolverIterative(realtype *t, AmiVector *x,
-                                             Model *model)
-    : NewtonSolver(t, x, model), ns_p_(model->nx_solver),
-    ns_h_(model->nx_solver), ns_t_(model->nx_solver), ns_s_(model->nx_solver),
-    ns_r_(model->nx_solver), ns_rt_(model->nx_solver), ns_v_(model->nx_solver),
-    ns_Jv_(model->nx_solver), ns_tmp_(model->nx_solver),
-    ns_Jdiag_(model->nx_solver), ns_J_(model->nx_solver, model->nx_solver)
-    {
-}
-
-/* ------------------------------------------------------------------------- */
-
-void NewtonSolverIterative::prepareLinearSystem(int ntry, int nnewt) {
-    newton_try_ = ntry;
-    i_newton_ = nnewt;
-    if (nnewt == -1) {
-        throw NewtonFailure(AMICI_NOT_IMPLEMENTED,
-                            "Linear solver SPBCG does not support sensitivity "
-                            "computation for steady state problems.");
-    }
-
-    // Get the Jacobian and its diagonal for preconditioning
-    model_->fJ(*t_, 0.0, *x_, dx_, xdot_, ns_J_.get());
-    ns_J_.refresh();
-    model_->fJDiag(*t_, ns_Jdiag_, 0.0, *x_, dx_);
-
-    // Ensure positivity of entries in ns_Jdiag
-    ns_p_.set(1.0);
-    ns_Jdiag_.abs();
-    N_VCompare(1e-15, ns_Jdiag_.getNVector(), ns_tmp_.getNVector());
-    linearSum(-1.0, ns_tmp_, 1.0, ns_p_, ns_tmp_);
-    linearSum(1.0, ns_Jdiag_, 1.0, ns_tmp_, ns_Jdiag_);
-}
-
-/* ------------------------------------------------------------------------- */
-
-void NewtonSolverIterative::prepareLinearSystemB(int ntry, int nnewt) {
-    newton_try_ = ntry;
-    i_newton_ = nnewt;
-    if (nnewt == -1) {
-        throw AmiException("Linear solver SPBCG does not support sensitivity "
-                           "computation for steady state problems.");
-    }
-
-    // Get the Jacobian and its diagonal for preconditioning
-    model_->fJB(*t_, 0.0, *x_, dx_, xB_, dxB_, xdot_, ns_J_.get());
-    ns_J_.refresh();
-    // Get the diagonal and ensure negativity of entries is ns_J. Note that diag(JB) = -diag(J).
-    model_->fJDiag(*t_, ns_Jdiag_, 0.0, *x_, dx_);
-
-    ns_p_.set(1.0);
-    ns_Jdiag_.abs();
-    N_VCompare(1e-15, ns_Jdiag_.getNVector(), ns_tmp_.getNVector());
-    linearSum(-1.0, ns_tmp_, 1.0, ns_p_, ns_tmp_);
-    linearSum(1.0, ns_Jdiag_, 1.0, ns_tmp_, ns_Jdiag_);
-    ns_Jdiag_.minus();
-}
-
-/* ------------------------------------------------------------------------- */
-
-void NewtonSolverIterative::solveLinearSystem(AmiVector &rhs) {
-    linsolveSPBCG(newton_try_, i_newton_, rhs);
-    rhs.minus();
-}
-
-/* ------------------------------------------------------------------------- */
-
-void NewtonSolverIterative::linsolveSPBCG(int /*ntry*/, int nnewt,
-                                          AmiVector &ns_delta) {
-    xdot_ = ns_delta;
-    xdot_.minus();
-
-    // Initialize for linear solve
-    ns_p_.zero();
-    ns_v_.zero();
-    ns_delta.zero();
-    ns_tmp_.zero();
-    double rho = 1.0;
-    double omega = 1.0;
-    double alpha = 1.0;
-
-    ns_J_.multiply(ns_Jv_, ns_delta);
-
-    // ns_r = xdot - ns_Jv;
-    linearSum(-1.0, ns_Jv_, 1.0, xdot_, ns_r_);
-    ns_r_ /= ns_Jdiag_;
-    ns_rt_ = ns_r_;
-
-    for (int i_linstep = 0; i_linstep < max_lin_steps_; i_linstep++) {
-        // Compute factors
-        double rho1 = rho;
-        rho = dotProd(ns_rt_, ns_r_);
-        double beta = rho * alpha / (rho1 * omega);
-
-        // ns_p = ns_r + beta * (ns_p - omega * ns_v);
-        linearSum(1.0, ns_p_, -omega, ns_v_, ns_p_);
-        linearSum(1.0, ns_r_, beta, ns_p_, ns_p_);
-
-        // ns_v = J * ns_p
-        ns_v_.zero();
-        ns_J_.multiply(ns_v_, ns_p_);
-        ns_v_ /= ns_Jdiag_;
-
-        // Compute factor
-        alpha = rho / dotProd(ns_rt_, ns_v_);
-
-        // ns_h = ns_delta + alpha * ns_p;
-        linearSum(1.0, ns_delta, alpha, ns_p_, ns_h_);
-        // ns_s = ns_r - alpha * ns_v;
-        linearSum(1.0, ns_r_, -alpha, ns_v_, ns_s_);
-
-        // ns_t = J * ns_s
-        ns_t_.zero();
-        ns_J_.multiply(ns_t_, ns_s_);
-        ns_t_ /= ns_Jdiag_;
-
-        // Compute factor
-        omega = dotProd(ns_t_, ns_s_) / dotProd(ns_t_, ns_t_);
-
-        // ns_delta = ns_h + omega * ns_s;
-        linearSum(1.0, ns_h_, omega, ns_s_, ns_delta);
-        // ns_r = ns_s - omega * ns_t;
-        linearSum(1.0, ns_s_, -omega, ns_t_, ns_r_);
-
-        // Compute the (unscaled) residual
-        ns_r_ *= ns_Jdiag_;
-        double res = sqrt(dotProd(ns_r_, ns_r_));
-
-        // Test convergence
-        if (res < atol_) {
-            // Write number of steps needed
-            num_lin_steps_.at(nnewt) = i_linstep + 1;
-
-            // Return success
-            return;
-        }
-
-        // Scale back
-        ns_r_ /= ns_Jdiag_;
-    }
-    throw NewtonFailure(AMICI_CONV_FAILURE, "linsolveSPBCG");
-}
-
 
 } // namespace amici

src/steadystateproblem.cpp

@@ -26,11 +26,6 @@ SteadystateProblem::SteadystateProblem(const Solver &solver, const Model &model)
       xB_(model.nJ * model.nx_solver), xQ_(model.nJ * model.nx_solver),
       xQB_(model.nplist()), xQBdot_(model.nplist()),
       dJydx_(model.nJ * model.nx_solver * model.nt(), 0.0) {
-          /* maxSteps must be adapted if iterative linear solvers are used */
-          if (solver.getLinearSolver() == LinearSolver::SPBCG) {
-              max_steps_ = solver.getNewtonMaxSteps();
-              numlinsteps_.resize(2 * max_steps_, 0);
-          }
           /* Check for compatibility of options */
           if (solver.getSensitivityMethod() == SensitivityMethod::forward &&
               solver.getSensitivityMethodPreequilibration() == SensitivityMethod::adjoint &&

tests/cpp/steadystate/tests1.cpp

@@ -174,12 +174,6 @@ TEST(ExampleSteadystate, SensitivityForwardDense)
     amici::simulateVerifyWrite("/model_steadystate/sensiforwarddense/");
 }
 
-TEST(ExampleSteadystate, SensitivityForwardSPBCG)
-{
-    amici::simulateVerifyWrite(
-      "/model_steadystate/nosensiSPBCG/", 10 * TEST_ATOL, 10 * TEST_RTOL);
-}
-
 TEST(ExampleSteadystate, SensiFwdNewtonPreeq)
 {
     amici::simulateVerifyWrite("/model_steadystate/sensifwdnewtonpreeq/");

tests/generateTestConfig/example_steadystate.py

@@ -54,16 +54,6 @@ def writeSensiForwardDense(filename):
     ex.writeToFile(filename, '/model_steadystate/sensiforwarddense/')
 
 
-def writeNosensiSPBCG(filename):
-    ex = ExampleSteadystate()
-
-    ex.modelOptions['ts'] = np.append(np.linspace(0, 100, 50), np.inf)
-    ex.solverOptions['sensi'] = 0
-    ex.solverOptions['linsol'] = 7
-
-    ex.writeToFile(filename, '/model_steadystate/nosensiSPBCG/')
-
-
 def writeSensiForwardErrorInt(filename):
     ex = ExampleSteadystate()
 
@@ -379,7 +369,6 @@ def main():
     writeSensiForward(filename)
     writeSensiForwardPlist(filename)
     writeSensiForwardDense(filename)
-    writeNosensiSPBCG(filename)
     writeSensiForwardErrorInt(filename)
     writeSensiForwardErrorNewt(filename)
     writeSensiFwdNewtonPreeq(filename)