From d8d164c69c5a06cca504cb4c75892583c89da490 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fabian=20Fr=C3=B6hlich?= Date: Wed, 16 Mar 2022 16:47:30 -0400 Subject: [PATCH] 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 | 90 ----------- python/tests/test_preequilibration.py | 16 -- src/newton_solver.cpp | 152 +----------------- src/steadystateproblem.cpp | 5 - tests/cpp/expectedResults.h5 | Bin 4198368 -> 4198368 bytes tests/cpp/steadystate/tests1.cpp | 6 - tests/cpp/testOptions.h5 | Bin 345728 -> 345728 bytes .../generateTestConfig/example_steadystate.py | 11 -- 8 files changed, 1 insertion(+), 279 deletions(-) diff --git a/include/amici/newton_solver.h b/include/amici/newton_solver.h index 01c033aeb6..9dd3c8a92f 100644 --- a/include/amici/newton_solver.h +++ b/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 diff --git a/python/tests/test_preequilibration.py b/python/tests/test_preequilibration.py index 2af6da888b..5cd67a39d3 100644 --- a/python/tests/test_preequilibration.py +++ b/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 diff --git a/src/newton_solver.cpp b/src/newton_solver.cpp index ca032d982a..f10006ff92 100644 --- a/src/newton_solver.cpp +++ b/src/newton_solver.cpp @@ -49,8 +49,7 @@ std::unique_ptr 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::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 diff --git a/src/steadystateproblem.cpp b/src/steadystateproblem.cpp index d7c3850e64..1bdae812e8 100644 --- a/src/steadystateproblem.cpp +++ b/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 && diff --git a/tests/cpp/expectedResults.h5 b/tests/cpp/expectedResults.h5 index a24be851566e794a6dffd0f0564bed2660e06996..a68de6ea7eda08275a442ed7b639639b8429aacd 100644 GIT binary patch delta 314 zcmW;GJ5B;o06kY#mOxMG0+eJzE7g=NfZ9xJsAvfr zUZ7I*UlFNa%b~&qoEV+fWe}}SQzm4SI{3ky% zObM@Jm_GkBFcPjqx;@_XA7+!S4}Js?#1TR`#tFjE5y2_W5QV}yVz@vYm$*U#NnGOw yw@4w446?XG4)@5TfFk~%N_fN*$}mtt6*bh+z%!a?;RS7U&_xe@XR>{L@9qKhI))hl delta 345 zcmXZPJ5B;o0EJ=hz=aV(6u3T7L*>1a$`fVCNG zD|;GCc3?&CAK)Zka{k$SUz=eZ7LI|kl(3X$DW|1epGtl-G*j5$4ViWd_vEiBQ}G|i za^sx7IAYabtNzXFlaMLrBcmnZwvCoYANsnly-v6H(l)El^-*o=fsQqJu?`>n*gycA z*g_B?grTsF9qeKc`#3-Zhd4qM$B5wsahxK7B+hV-6fSUyG_G)s8yL979q#df46?}K S5qT8wgd$2P+qG%st@sZ$&4B6v diff --git a/tests/cpp/steadystate/tests1.cpp b/tests/cpp/steadystate/tests1.cpp index 505f8f15ed..93f48ceb27 100644 --- a/tests/cpp/steadystate/tests1.cpp +++ b/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/"); diff --git a/tests/cpp/testOptions.h5 b/tests/cpp/testOptions.h5 index e2bf3f658e83688bab6aa03e5fe7a117d1360861..25280e9cec208ed2aaecb54b27a2106c61723607 100644 GIT binary patch delta 73 zcmZqZ6>R{b7RDB)EzA=`85y=u3T4(|W@MkfKZ03vxR{b7RDB)EzA=`85y@v3T4(|X5^T@KZ03PDS!b2PV8m|GZ>-#4^Y0qzUhx5 jn7?g{Vm`5+CojJ^HLo}`IKau-z5NO&%l0drtZQrl-;o