diff --git a/cmake/recipes/mkl.cmake b/cmake/recipes/mkl.cmake index 999b096..814b351 100644 --- a/cmake/recipes/mkl.cmake +++ b/cmake/recipes/mkl.cmake @@ -34,8 +34,8 @@ set_property(CACHE MKL_LINKING PROPERTY STRINGS ${MKL_LINKINK_CHOICES}) message(STATUS "MKL linking strategy: ${MKL_LINKING}") # MKL version -set(MKL_VERSION "2021.3.0" CACHE STRING "MKL version to use (2020.4 or 2021.3.0)") -set(MKL_VERSION_CHOICES 2020.4 2021.3.0) +set(MKL_VERSION "2022.2.1" CACHE STRING "MKL version to use (2022.2.1)") +set(MKL_VERSION_CHOICES 2022.2.1) set_property(CACHE MKL_VERSION PROPERTY STRINGS ${MKL_VERSION_CHOICES}) message(STATUS "MKL version: ${MKL_VERSION}") @@ -78,66 +78,34 @@ elseif(UNIX) set(MKL_PLATFORM linux-64) endif() -if(MKL_VERSION VERSION_EQUAL 2020.4) - # To compute the md5 checksums for each lib, use the following bash script (replace the target version number): - # for f in mkl mkl-include mkl-static mkl-devel; do for os in linux osx win; do cat <(printf "$f-$os-64-md5") <(conda search --override-channel --channel intel $f=2020.4 --platform $os-64 -i | grep md5 | cut -d : -f 2); done; done - set(mkl-linux-64-md5 f85891f97a04c7b2fbf619d1b903d7f5) - set(mkl-osx-64-md5 735d7f93c7fbcffe658f1ccf67418cb3) - set(mkl-win-64-md5 37ade09cace5cd73053b16574a3ee3c3) - set(mkl-include-linux-64-md5 8b2bf0e42bd95dd700d9877add1ca6de) - set(mkl-include-osx-64-md5 26043328553952cdb064c5aab8b50c78) - set(mkl-include-win-64-md5 87e9a73a6e6757a8ed0dbc87d50d7f60) - set(mkl-static-linux-64-md5 9f589a1508fb083c3e73427db459ca4c) - set(mkl-static-osx-64-md5 2f9e1b8b6d6b0903e81a573084e4494f) - set(mkl-static-win-64-md5 5ae780c06edd0be62966c6d8ab47d5fb) - set(mkl-devel-linux-64-md5 b571698ef237c0e61abe15b7d300f157) - set(mkl-devel-osx-64-md5 ee58da0463676d910eeab9aec0470f0e) - set(mkl-devel-win-64-md5 8a7736b81b9bc2d5c044b88d6ac8af6e) - - # To compute file names, we use the following bash script: - # for f in mkl mkl-include mkl-static mkl-devel; do for os in linux osx win; do cat <(printf "$f-$os-64-file") <(conda search --override-channel --channel intel $f=2020.4 --platform $os-64 -i | grep file | cut -d : -f 2); done; done - set(mkl-linux-64-file mkl-2020.4-intel_304.tar.bz2) - set(mkl-osx-64-file mkl-2020.4-intel_301.tar.bz2) - set(mkl-win-64-file mkl-2020.4-intel_311.tar.bz2) - set(mkl-include-linux-64-file mkl-include-2020.4-intel_304.tar.bz2) - set(mkl-include-osx-64-file mkl-include-2020.4-intel_301.tar.bz2) - set(mkl-include-win-64-file mkl-include-2020.4-intel_311.tar.bz2) - set(mkl-static-linux-64-file mkl-static-2020.4-intel_304.tar.bz2) - set(mkl-static-osx-64-file mkl-static-2020.4-intel_301.tar.bz2) - set(mkl-static-win-64-file mkl-static-2020.4-intel_311.tar.bz2) - set(mkl-devel-linux-64-file mkl-devel-2020.4-intel_304.tar.bz2) - set(mkl-devel-osx-64-file mkl-devel-2020.4-intel_301.tar.bz2) - set(mkl-devel-win-64-file mkl-devel-2020.4-intel_311.tar.bz2) -elseif(MKL_VERSION VERSION_EQUAL 2021.3.0) +if(MKL_VERSION VERSION_EQUAL 2022.2.1) # To compute the md5 checksums for each lib, use the following bash script (replace the target version number): # for f in mkl mkl-include mkl-static mkl-devel; do for os in linux osx win; do cat <(printf "$f-$os-64-md5") <(conda search --override-channel --channel intel $f=2021.3.0 --platform $os-64 -i | grep md5 | cut -d : -f 2); done; done - set(mkl-linux-64-md5 2501643729c00b24fddb9530b339aea7) - set(mkl-osx-64-md5 d6129ae9dfba58671667a65c160d0776) - set(mkl-win-64-md5 264213ea4c5cb6b6d81ea97f59e757ab) - set(mkl-include-linux-64-md5 70b4f9a53401a3d11ce27d7ddb0e2511) - set(mkl-include-osx-64-md5 6da50c06992b78c4127a1881d39c1804) - set(mkl-include-win-64-md5 28d785eb22d28512d4e40e5890a817dc) - set(mkl-static-linux-64-md5 1469ad60a34269d4d0c5666bc131b82a) - set(mkl-static-osx-64-md5 4a099581ba95cc50bb538598b26389e4) - set(mkl-static-win-64-md5 69aef10428893314bc486e81397e1b25) - set(mkl-devel-linux-64-md5 2432ad963e3f7e4619ffc7f896178fbe) - set(mkl-devel-osx-64-md5 61b84a60715a3855a2097a3b619a00c8) - set(mkl-devel-win-64-md5 6128dee67d2b20ff534cf54757f623e0) + # set(mkl-linux-64-md5 2501643729c00b24fddb9530b339aea7) + # set(mkl-win-64-md5 264213ea4c5cb6b6d81ea97f59e757ab) + + # set(mkl-include-linux-64-md5 70b4f9a53401a3d11ce27d7ddb0e2511) + # set(mkl-include-win-64-md5 28d785eb22d28512d4e40e5890a817dc) + + # set(mkl-static-linux-64-md5 1469ad60a34269d4d0c5666bc131b82a) + # set(mkl-static-win-64-md5 69aef10428893314bc486e81397e1b25) + + # set(mkl-devel-linux-64-md5 2432ad963e3f7e4619ffc7f896178fbe) + # set(mkl-devel-win-64-md5 6128dee67d2b20ff534cf54757f623e0) # To compute file names, we use the following bash script: # for f in mkl mkl-include mkl-static mkl-devel; do for os in linux osx win; do cat <(printf "$f-$os-64-file") <(conda search --override-channel --channel intel $f=2021.3.0 --platform $os-64 -i | grep file | cut -d : -f 2); done; done - set(mkl-linux-64-file mkl-2021.3.0-intel_520.tar.bz2) - set(mkl-osx-64-file mkl-2021.3.0-intel_517.tar.bz2) - set(mkl-win-64-file mkl-2021.3.0-intel_524.tar.bz2) - set(mkl-include-linux-64-file mkl-include-2021.3.0-intel_520.tar.bz2) - set(mkl-include-osx-64-file mkl-include-2021.3.0-intel_517.tar.bz2) - set(mkl-include-win-64-file mkl-include-2021.3.0-intel_524.tar.bz2) - set(mkl-static-linux-64-file mkl-static-2021.3.0-intel_520.tar.bz2) - set(mkl-static-osx-64-file mkl-static-2021.3.0-intel_517.tar.bz2) - set(mkl-static-win-64-file mkl-static-2021.3.0-intel_524.tar.bz2) - set(mkl-devel-linux-64-file mkl-devel-2021.3.0-intel_520.tar.bz2) - set(mkl-devel-osx-64-file mkl-devel-2021.3.0-intel_517.tar.bz2) - set(mkl-devel-win-64-file mkl-devel-2021.3.0-intel_524.tar.bz2) + set(mkl-linux-64-file mkl-2022.2.1-h6508926_16999.tar.bz2) + set(mkl-win-64-file mkl-2022.2.1-h4060db9_19760.tar.bz2) + + set(mkl-include-linux-64-file mkl-include-2022.2.1-ha957f24_16999.tar.bz2) + set(mkl-include-win-64-file mkl-include-2022.2.1-h66d3029_19760.tar.bz2) + + set(mkl-static-linux-64-file mkl-static-2022.2.1-h6508926_16999.tar.bz2) + set(mkl-static-win-64-file mkl-static-2022.2.1-h4060db9_19760.tar.bz2) + + set(mkl-devel-linux-64-file mkl-devel-2022.2.1-ha957f24_16999.tar.bz2) + set(mkl-devel-win-64-file mkl-devel-2022.2.1-h66d3029_19760.tar.bz2) endif() # On Windows, `mkl-devel` contains the .lib files (needed at link time), @@ -153,8 +121,8 @@ include(CPM) foreach(name IN ITEMS ${MKL_REMOTES}) CPMAddPackage( NAME ${name} - URL https://anaconda.org/intel/${name}/${MKL_VERSION}/download/${MKL_PLATFORM}/${${name}-${MKL_PLATFORM}-file} - URL_MD5 ${${name}-${MKL_PLATFORM}-md5} + URL https://anaconda.org/conda-forge/${name}/${MKL_VERSION}/download/${MKL_PLATFORM}/${${name}-${MKL_PLATFORM}-file} + # URL_MD5 ${${name}-${MKL_PLATFORM}-md5} ) endforeach() diff --git a/src/polysolve/linear/Solver.hpp b/src/polysolve/linear/Solver.hpp index deaae01..58e281f 100644 --- a/src/polysolve/linear/Solver.hpp +++ b/src/polysolve/linear/Solver.hpp @@ -51,22 +51,22 @@ namespace polysolve::linear /// Selects the correct solver based on params using the fallback or the list of solvers if necessary static void select_valid_solver(json ¶ms, spdlog::logger &logger); - // Static constructor - // - // @param[in] params Parameter of the solver, including name and preconditioner - // @param[in] logger Logger used for error - // @param[in] strict_validation strict validation of the input paraams - // @param[out] a pointer to a linear solver + /// @brief Static constructor + /// + /// @param[in] params Parameter of the solver, including name and preconditioner + /// @param[in] logger Logger used for error + /// @param[in] strict_validation strict validation of the input paraams + /// @return a pointer to a linear solver // static std::unique_ptr create(const json ¶ms, spdlog::logger &logger, const bool strict_validation = true); - // Static constructor - // - // @param[in] solver Solver type - // @param[in] precond Preconditioner for iterative solvers - // + /// @brief Static constructor + /// + /// @param[in] solver Solver type + /// @param[in] precond Preconditioner for iterative solvers + /// static std::unique_ptr create(const std::string &solver, const std::string &precond); // List available solvers @@ -86,44 +86,45 @@ namespace polysolve::linear // Public interface // ////////////////////// - // Set solver parameters + /// Set solver parameters virtual void set_parameters(const json ¶ms) {} - // Get info on the last solve step + /// Get info on the last solve step virtual void get_info(json ¶ms) const {}; - // Analyze sparsity pattern + /// Analyze sparsity pattern virtual void analyze_pattern(const StiffnessMatrix &A, const int precond_num) {} - // Factorize system matrix + /// Factorize system matrix virtual void factorize(const StiffnessMatrix &A) {} - // Analyze sparsity pattern of a dense matrix + /// Analyze sparsity pattern of a dense matrix virtual void analyze_pattern_dense(const Eigen::MatrixXd &A, const int precond_num) {} - // Factorize system matrix of a dense matrix + /// Factorize system matrix of a dense matrix virtual void factorize_dense(const Eigen::MatrixXd &A) {} - // If solver uses dense matrices + /// If solver uses dense matrices virtual bool is_dense() const { return false; } - // - // @brief { Solve the linear system Ax = b } - // - // @param[in] b { Right-hand side. } - // @param[in,out] x { Unknown to compute. When using an iterative - // solver, the input unknown vector is used as an - // initial guess, and must thus be properly allocated - // and initialized. } - // + /// Set block size for multigrid solvers + virtual void set_block_size(int block_size) {} + + /// If the problem is nullspace for multigrid solvers + virtual void set_is_nullspace(const VectorXd &x) {} + + /// + /// @brief { Solve the linear system Ax = b } + /// + /// @param[in] b { Right-hand side. } + /// @param[in,out] x { Unknown to compute. When using an iterative + /// solver, the input unknown vector is used as an + /// initial guess, and must thus be properly allocated + /// and initialized. } + /// virtual void solve(const Ref b, Ref x) = 0; - public: - /////////// - // Debug // - /////////// - - // Name of the solver type (for debugging purposes) + /// @brief Name of the solver type (for debugging purposes) virtual std::string name() const { return ""; } }; diff --git a/src/polysolve/nonlinear/Criteria.cpp b/src/polysolve/nonlinear/Criteria.cpp index 37fed88..97ef602 100644 --- a/src/polysolve/nonlinear/Criteria.cpp +++ b/src/polysolve/nonlinear/Criteria.cpp @@ -2,6 +2,8 @@ // License: MIT #include "Criteria.hpp" +#include + namespace polysolve::nonlinear { bool is_converged_status(const Status s) @@ -27,13 +29,9 @@ namespace polysolve::nonlinear void Criteria::print(std::ostream &os) const { - os << "Iterations: " << iterations << std::endl; - os << "xDelta: " << xDelta << std::endl; - os << "fDelta: " << fDelta << std::endl; - os << "GradNorm: " << gradNorm << std::endl; - os << "xDeltaDotGrad: " << xDeltaDotGrad << std::endl; - os << "FirstGradNorm: " << firstGradNorm << std::endl; - os << "fDeltaCount: " << fDeltaCount << std::endl; + os << fmt::format( + "iters={:d} Δf={:g} ‖∇f‖={:g} ‖Δx‖={:g} Δx⋅∇f(x)={:g}", + iterations, fDelta, gradNorm, xDelta, xDeltaDotGrad); } Status checkConvergence(const Criteria &stop, const Criteria ¤t) @@ -42,6 +40,11 @@ namespace polysolve::nonlinear { return Status::IterationLimit; } + const double stopGradNorm = current.iterations == 0 ? stop.firstGradNorm : stop.gradNorm; + if (stopGradNorm > 0 && current.gradNorm < stopGradNorm) + { + return Status::GradNormTolerance; + } if (stop.xDelta > 0 && current.xDelta < stop.xDelta) { return Status::XDeltaTolerance; @@ -50,11 +53,6 @@ namespace polysolve::nonlinear { return Status::FDeltaTolerance; } - const double stopGradNorm = current.iterations == 0 ? stop.firstGradNorm : stop.gradNorm; - if (stopGradNorm > 0 && current.gradNorm < stopGradNorm) - { - return Status::GradNormTolerance; - } // Δx⋅∇f ≥ 0 means that the search direction is not a descent direction if (stop.xDeltaDotGrad < 0 && current.xDeltaDotGrad > stop.xDeltaDotGrad) { @@ -68,37 +66,40 @@ namespace polysolve::nonlinear switch (s) { case Status::NotStarted: - os << "Solver not started."; + os << "Solver not started"; break; case Status::Continue: - os << "Convergence criteria not reached."; + os << "Convergence criteria not reached"; break; case Status::IterationLimit: - os << "Iteration limit reached."; + os << "Iteration limit reached"; break; case Status::XDeltaTolerance: - os << "Change in parameter vector too small."; + os << "Change in parameter vector too small"; break; case Status::FDeltaTolerance: - os << "Change in cost function value too small."; + os << "Change in cost function value too small"; break; case Status::GradNormTolerance: - os << "Gradient vector norm too small."; + os << "Gradient vector norm too small"; break; case Status::ObjectiveCustomStop: - os << "Objective function specified to stop."; + os << "Objective function specified to stop"; break; case Status::NanEncountered: - os << "Objective or gradient function returned NaN."; + os << "Objective or gradient function returned NaN"; break; case Status::NotDescentDirection: - os << "Search direction not a descent direction."; + os << "Search direction not a descent direction"; break; case Status::LineSearchFailed: - os << "Line search failed."; + os << "Line search failed"; + break; + case Status::UpdateDirectionFailed: + os << "Update direction could not be computed"; break; default: - os << "Unknown status."; + os << "Unknown status"; break; } return os; diff --git a/src/polysolve/nonlinear/Criteria.hpp b/src/polysolve/nonlinear/Criteria.hpp index f318ed5..86a89e1 100644 --- a/src/polysolve/nonlinear/Criteria.hpp +++ b/src/polysolve/nonlinear/Criteria.hpp @@ -19,9 +19,10 @@ namespace polysolve::nonlinear GradNormTolerance, ///< The norm of the gradient vector is below the tolerance ObjectiveCustomStop, ///< The objective function specified to stop // Failure cases - NanEncountered, ///< The objective function returned NaN - NotDescentDirection, ///< The search direction is not a descent direction - LineSearchFailed, ///< The line search failed + NanEncountered, ///< The objective function returned NaN + NotDescentDirection, ///< The search direction is not a descent direction + LineSearchFailed, ///< The line search failed + UpdateDirectionFailed, ///< The update direction could not be computed }; bool is_converged_status(const Status s); diff --git a/src/polysolve/nonlinear/Problem.hpp b/src/polysolve/nonlinear/Problem.hpp index 8c50113..b71e212 100644 --- a/src/polysolve/nonlinear/Problem.hpp +++ b/src/polysolve/nonlinear/Problem.hpp @@ -10,6 +10,7 @@ namespace polysolve::nonlinear { + /// @brief Class defining optimization problem to be solved. To be defined by user code class Problem { public: diff --git a/src/polysolve/nonlinear/Solver.cpp b/src/polysolve/nonlinear/Solver.cpp index 976115e..990c7b9 100644 --- a/src/polysolve/nonlinear/Solver.cpp +++ b/src/polysolve/nonlinear/Solver.cpp @@ -205,6 +205,7 @@ namespace polysolve::nonlinear m_stop.xDeltaDotGrad = -solver_params["advanced"]["derivative_along_delta_x_tol"].get(); // Make these relative to the characteristic length + m_logger.trace("Using a characteristic length of {:g}", characteristic_length); m_stop.xDelta *= characteristic_length; m_stop.fDelta *= characteristic_length; m_stop.gradNorm *= characteristic_length; @@ -277,11 +278,8 @@ namespace polysolve::nonlinear stop_watch.start(); m_logger.debug( - "Starting {} with {} solve f₀={:g} " - "(stopping criteria: max_iters={:d} Δf={:g} ‖∇f‖={:g} ‖Δx‖={:g})", - descent_strategy_name(), m_line_search->name(), - objFunc(x), m_stop.iterations, - m_stop.fDelta, m_stop.gradNorm, m_stop.xDelta); + "Starting {} with {} solve f₀={:g} (stopping criteria: {})", + descent_strategy_name(), m_line_search->name(), objFunc(x), m_stop); update_solver_info(objFunc(x)); objFunc.post_step(PostStepData(m_current.iterations, solver_info, x, grad)); @@ -325,56 +323,69 @@ namespace polysolve::nonlinear log_and_throw_error(m_logger, "[{}][{}] Gradient is nan; stopping", descent_strategy_name(), m_line_search->name()); } + // Check convergence without these values to avoid impossible linear solves. + m_current.xDelta = NaN; + m_current.xDeltaDotGrad = NaN; m_status = checkConvergence(m_stop, m_current); if (m_status != Status::Continue) break; // --- Update direction -------------------------------------------- - const bool ok = compute_update_direction(objFunc, x, grad, delta_x); - m_current.xDeltaDotGrad = delta_x.dot(grad); - - if (!ok || (m_strategies[m_descent_strategy]->is_direction_descent() && m_current.gradNorm != 0 && m_current.xDeltaDotGrad >= 0)) + bool update_direction_successful; { - const std::string current_name = descent_strategy_name(); + POLYSOLVE_SCOPED_STOPWATCH("compute update direction", update_direction_time, m_logger); + update_direction_successful = compute_update_direction(objFunc, x, grad, delta_x); + } + m_current.xDelta = delta_x.norm(); + if (!update_direction_successful || std::isnan(m_current.xDelta)) + { + const auto current_name = descent_strategy_name(); if (!m_strategies[m_descent_strategy]->handle_error()) ++m_descent_strategy; if (m_descent_strategy >= m_strategies.size()) { - m_status = Status::NotDescentDirection; + m_status = Status::UpdateDirectionFailed; log_and_throw_error( - m_logger, "[{}][{}] direction is not a descent direction on last strategy (‖Δx‖={:g}; ‖g‖={:g}; Δx⋅g={:g}≥0); stopping", - current_name, m_line_search->name(), delta_x.norm(), compute_grad_norm(x, grad), m_current.xDeltaDotGrad); - } - else - { - m_status = Status::Continue; - m_logger.debug( - "[{}][{}] direction is not a descent direction (‖Δx‖={:g}; ‖g‖={:g}; Δx⋅g={:g}≥0); reverting to {}", - current_name, m_line_search->name(), delta_x.norm(), compute_grad_norm(x, grad), m_current.xDeltaDotGrad, - descent_strategy_name()); + m_logger, "[{}][{}] {} on last strategy; stopping", + current_name, m_line_search->name(), m_status); } + + m_logger.debug( + "[{}][{}] {}; reverting to {}", current_name, m_line_search->name(), + Status::UpdateDirectionFailed, descent_strategy_name()); + m_status = Status::Continue; continue; } - m_current.xDelta = delta_x.norm(); - if (std::isnan(m_current.xDelta)) + m_current.xDeltaDotGrad = delta_x.dot(grad); + + if (m_strategies[m_descent_strategy]->is_direction_descent() && m_current.gradNorm != 0 && m_current.xDeltaDotGrad >= 0) { - const auto current_name = descent_strategy_name(); + const std::string current_name = descent_strategy_name(); + if (!m_strategies[m_descent_strategy]->handle_error()) ++m_descent_strategy; if (m_descent_strategy >= m_strategies.size()) { - m_status = Status::NanEncountered; - log_and_throw_error(m_logger, "[{}][{}] Δx is nan on last strategy; stopping", - current_name, m_line_search->name()); + m_status = Status::NotDescentDirection; + log_and_throw_error( + m_logger, "[{}][{}] {} on last strategy (‖Δx‖={:g}; ‖g‖={:g}; Δx⋅g={:g}≥0); stopping", + current_name, m_line_search->name(), m_status, delta_x.norm(), compute_grad_norm(x, grad), + m_current.xDeltaDotGrad); + } + else + { + m_status = Status::Continue; + m_logger.debug( + "[{}][{}] {} (‖Δx‖={:g}; ‖g‖={:g}; Δx⋅g={:g}≥0); reverting to {}", + current_name, m_line_search->name(), Status::NotDescentDirection, + delta_x.norm(), compute_grad_norm(x, grad), m_current.xDeltaDotGrad, + descent_strategy_name()); } - - m_logger.debug("[{}][{}] Δx is nan; reverting to {}", current_name, m_line_search->name(), descent_strategy_name()); - m_status = Status::Continue; continue; } @@ -393,7 +404,12 @@ namespace polysolve::nonlinear m_current.iterations, energy, m_current.gradNorm); // Perform a line_search to compute step scale - double rate = m_line_search->line_search(x, delta_x, objFunc); + double rate; + { + POLYSOLVE_SCOPED_STOPWATCH("line search", line_search_time, m_logger); + rate = m_line_search->line_search(x, delta_x, objFunc); + } + if (std::isnan(rate)) { const auto current_name = descent_strategy_name(); @@ -421,7 +437,6 @@ namespace polysolve::nonlinear // if we did enough lower strategy, we revert back to normal if (m_descent_strategy != 0 && current_strategy_iter >= m_iter_per_strategy[m_descent_strategy]) { - const auto current_name = descent_strategy_name(); const std::string prev_strategy_name = descent_strategy_name(); @@ -442,6 +457,9 @@ namespace polysolve::nonlinear // ----------- const double step = (rate * delta_x).norm(); + // m_logger.debug("[{}][{}] rate={:g} ‖step‖={:g}", + // descent_strategy_name(), m_line_search->name(), rate, step); + update_solver_info(energy); objFunc.post_step(PostStepData(m_current.iterations, solver_info, x, grad)); @@ -454,12 +472,8 @@ namespace polysolve::nonlinear m_current.fDeltaCount = (m_current.fDelta < m_stop.fDelta) ? (m_current.fDeltaCount + 1) : 0; m_logger.debug( - "[{}][{}] iter={:d} f={:g} Δf={:g} ‖∇f‖={:g} ‖Δx‖={:g} Δx⋅∇f(x)={:g} rate={:g} ‖step‖={:g}" - " (stopping criteria: max_iters={:d} Δf={:g} ‖∇f‖={:g} ‖Δx‖={:g})", - descent_strategy_name(), m_line_search->name(), - m_current.iterations, energy, m_current.fDelta, - m_current.gradNorm, m_current.xDelta, delta_x.dot(grad), rate, step, - m_stop.iterations, m_stop.fDelta, m_stop.gradNorm, m_stop.xDelta); + "[{}][{}] {} (stopping criteria: {})", + descent_strategy_name(), m_line_search->name(), m_current, m_stop); if (++m_current.iterations >= m_stop.iterations) m_status = Status::IterationLimit; @@ -480,12 +494,9 @@ namespace polysolve::nonlinear const bool succeeded = m_status == Status::GradNormTolerance; m_logger.log( succeeded ? spdlog::level::info : spdlog::level::err, - "[{}][{}] Finished: {} Took {:g}s (niters={:d} f={:g} Δf={:g} ‖∇f‖={:g} ‖Δx‖={:g})" - " (stopping criteria: max_iters={:d} Δf={:g} ‖∇f‖={:g} ‖Δx‖={:g})", - descent_strategy_name(), m_line_search->name(), - m_status, tot_time, m_current.iterations, - old_energy, m_current.fDelta, m_current.gradNorm, m_current.xDelta, - m_stop.iterations, m_stop.fDelta, m_stop.gradNorm, m_stop.xDelta); + "[{}][{}] Finished: {} took {:g}s ({}) (stopping criteria: {})", + descent_strategy_name(), m_line_search->name(), m_status, tot_time, + m_current, m_stop); log_times(); update_solver_info(objFunc(x)); @@ -511,14 +522,13 @@ namespace polysolve::nonlinear void Solver::reset_times() { total_time = 0; + obj_fun_time = 0; grad_time = 0; + update_direction_time = 0; line_search_time = 0; - obj_fun_time = 0; constraint_set_update_time = 0; if (m_line_search) - { m_line_search->reset_times(); - } for (auto &s : m_strategies) s->reset_times(); } @@ -535,47 +545,30 @@ namespace polysolve::nonlinear double per_iteration = m_current.iterations ? m_current.iterations : 1; solver_info["total_time"] = total_time; + solver_info["time_obj_fun"] = obj_fun_time / per_iteration; solver_info["time_grad"] = grad_time / per_iteration; + // Do not save update_direction_time as it is redundant with the strategies solver_info["time_line_search"] = line_search_time / per_iteration; solver_info["time_constraint_set_update"] = constraint_set_update_time / per_iteration; - solver_info["time_obj_fun"] = obj_fun_time / per_iteration; for (auto &s : m_strategies) s->update_solver_info(solver_info, per_iteration); - if (m_line_search) - { - solver_info["line_search_iterations"] = m_line_search->iterations(); - - solver_info["time_checking_for_nan_inf"] = - m_line_search->checking_for_nan_inf_time / per_iteration; - solver_info["time_broad_phase_ccd"] = - m_line_search->broad_phase_ccd_time / per_iteration; - solver_info["time_ccd"] = m_line_search->ccd_time / per_iteration; - // Remove double counting - solver_info["time_classical_line_search"] = - (m_line_search->classical_line_search_time - - m_line_search->constraint_set_update_time) - / per_iteration; - solver_info["time_line_search_constraint_set_update"] = - m_line_search->constraint_set_update_time / per_iteration; - } + m_line_search->update_solver_info(solver_info, per_iteration); } - void Solver::log_times() + void Solver::log_times() const { m_logger.debug( - "[{}] grad {:.3g}s, " - "line_search {:.3g}s, constraint_set_update {:.3g}s, " - "obj_fun {:.3g}s, checking_for_nan_inf {:.3g}s, " - "broad_phase_ccd {:.3g}s, ccd {:.3g}s, " - "classical_line_search {:.3g}s", + "[{}] f: {:.2e}s, grad_f: {:.2e}s, update_direction: {:.2e}s, " + "line_search: {:.2e}s, constraint_set_update: {:.2e}s", fmt::format(fmt::fg(fmt::terminal_color::magenta), "timing"), - grad_time, line_search_time, - constraint_set_update_time + (m_line_search ? m_line_search->constraint_set_update_time : 0), - obj_fun_time, m_line_search ? m_line_search->checking_for_nan_inf_time : 0, - m_line_search ? m_line_search->broad_phase_ccd_time : 0, m_line_search ? m_line_search->ccd_time : 0, - m_line_search ? m_line_search->classical_line_search_time : 0); + obj_fun_time, grad_time, update_direction_time, line_search_time, + constraint_set_update_time); + for (auto &s : m_strategies) + s->log_times(); + if (m_line_search) + m_line_search->log_times(); } void Solver::verify_gradient(Problem &objFunc, const TVector &x, const TVector &grad) diff --git a/src/polysolve/nonlinear/Solver.hpp b/src/polysolve/nonlinear/Solver.hpp index c00f739..211ea98 100644 --- a/src/polysolve/nonlinear/Solver.hpp +++ b/src/polysolve/nonlinear/Solver.hpp @@ -41,7 +41,7 @@ namespace polysolve::nonlinear spdlog::logger &logger, const bool strict_validation = true); - // List available solvers + /// @brief List available solvers static std::vector available_solvers(); public: @@ -73,7 +73,7 @@ namespace polysolve::nonlinear Criteria &stop_criteria() { return m_stop; } const Criteria &stop_criteria() const { return m_stop; } const Criteria ¤t_criteria() const { return m_current; } - const Status &status() const { return m_status; } + Status status() const { return m_status; } void set_strategies_iterations(const json &solver_params); void set_line_search(const json ¶ms); @@ -82,10 +82,17 @@ namespace polysolve::nonlinear /// @brief If true the solver will not throw an error if the maximum number of iterations is reached bool allow_out_of_iterations = false; + /// @brief Get the line search object const std::shared_ptr &line_search() const { return m_line_search; }; protected: + /// @brief Compute direction in which the argument should be updated + /// @param objFunc Problem to be minimized + /// @param x Current input (n x 1) + /// @param grad Gradient at current step (n x 1) + /// @param[out] direction Current update direction (n x 1) + /// @return True if update direction was found, False otherwises virtual bool compute_update_direction( Problem &objFunc, const TVector &x, @@ -116,6 +123,10 @@ namespace polysolve::nonlinear FiniteDiffStrategy gradient_fd_strategy = FiniteDiffStrategy::NONE; double gradient_fd_eps = 1e-7; + /// @brief Check gradient versus finite difference results + /// @param objFunc Problem defining relevant objective function + /// @param x Current input (n x 1) + /// @param grad Current gradient (n x 1) virtual void verify_gradient(Problem &objFunc, const TVector &x, const TVector &grad) final; private: @@ -128,8 +139,9 @@ namespace polysolve::nonlinear // ==================================================================== // Solver state // ==================================================================== - - // Reset the solver at the start of a minimization + + /// @brief Reset the solver at the start of a minimization + /// @param ndof number of degrees of freedom void reset(const int ndof); std::string descent_strategy_name() const { return m_strategies[m_descent_strategy]->name(); }; @@ -143,18 +155,25 @@ namespace polysolve::nonlinear // Solver info // ==================================================================== + /// @brief Update solver info JSON object + /// @param energy void update_solver_info(const double energy); + + /// @brief Reset timing members to 0 void reset_times(); - void log_times(); + + /// @brief Log time taken in different phases of the solve + void log_times() const; json solver_info; // Timers double total_time; + double obj_fun_time; double grad_time; + double update_direction_time; double line_search_time; double constraint_set_update_time; - double obj_fun_time; // ==================================================================== // END diff --git a/src/polysolve/nonlinear/descent_strategies/DescentStrategy.hpp b/src/polysolve/nonlinear/descent_strategies/DescentStrategy.hpp index 6e3f857..e8d0404 100644 --- a/src/polysolve/nonlinear/descent_strategies/DescentStrategy.hpp +++ b/src/polysolve/nonlinear/descent_strategies/DescentStrategy.hpp @@ -13,6 +13,10 @@ namespace polysolve::nonlinear using TVector = Problem::TVector; using Scalar = Problem::Scalar; + /// @brief Constructor + /// @param solver_params_ JSON of solver parameters + /// @param characteristic_length + /// @param logger DescentStrategy(const json &solver_params_, const double characteristic_length, spdlog::logger &logger) @@ -24,11 +28,22 @@ namespace polysolve::nonlinear virtual void reset(const int ndof) {} virtual void reset_times() {} + + /// @brief Update solver info after finding descent direction + /// @param solver_info JSON of solver parameters + /// @param per_iteration Number of iterations (used to normalize timings) virtual void update_solver_info(json &solver_info, const double per_iteration) {} + virtual void log_times() const {} virtual bool is_direction_descent() { return true; } virtual bool handle_error() { return false; } + /// @brief Compute descent direction along which to do line search + /// @param objFunc Problem to be minimized + /// @param x Current input (n x 1) + /// @param grad Current gradient (n x 1) + /// @param direction Current descent direction (n x 1) + /// @return True if a descent direction was successfully found virtual bool compute_update_direction( Problem &objFunc, const TVector &x, diff --git a/src/polysolve/nonlinear/descent_strategies/Newton.cpp b/src/polysolve/nonlinear/descent_strategies/Newton.cpp index 2fea6d6..4dce8f5 100644 --- a/src/polysolve/nonlinear/descent_strategies/Newton.cpp +++ b/src/polysolve/nonlinear/descent_strategies/Newton.cpp @@ -2,6 +2,8 @@ #include +#include + namespace polysolve::nonlinear { @@ -57,14 +59,14 @@ namespace polysolve::nonlinear const double characteristic_length, spdlog::logger &logger) : Superclass(solver_params, characteristic_length, logger), - is_sparse(sparse), m_characteristic_length(characteristic_length), m_residual_tolerance(residual_tolerance) + is_sparse(sparse), characteristic_length(characteristic_length), residual_tolerance(residual_tolerance) { linear_solver = polysolve::linear::Solver::create(linear_solver_params, logger); if (linear_solver->is_dense() == sparse) log_and_throw_error(logger, "Newton linear solver must be {}, instead got {}", sparse ? "sparse" : "dense", linear_solver->name()); - if (m_residual_tolerance <= 0) - log_and_throw_error(logger, "Newton residual_tolerance must be > 0, instead got {}", m_residual_tolerance); + if (residual_tolerance <= 0) + log_and_throw_error(logger, "Newton residual_tolerance must be > 0, instead got {}", residual_tolerance); } Newton::Newton( @@ -139,10 +141,10 @@ namespace polysolve::nonlinear is_sparse ? solve_sparse_linear_system(objFunc, x, grad, direction) : solve_dense_linear_system(objFunc, x, grad, direction); - if (std::isnan(residual) || residual > m_residual_tolerance * m_characteristic_length) + if (std::isnan(residual) || residual > residual_tolerance * characteristic_length) { - m_logger.debug("[{}] large (or nan) linear solve residual {}>{} (||∇f||={})", - name(), residual, m_residual_tolerance * m_characteristic_length, grad.norm()); + m_logger.debug("[{}] large (or nan) linear solve residual {}>{} (‖∇f‖={})", + name(), residual, residual_tolerance * characteristic_length, grad.norm()); return false; } @@ -156,8 +158,6 @@ namespace polysolve::nonlinear // ======================================================================= - // ======================================================================= - double Newton::solve_sparse_linear_system(Problem &objFunc, const TVector &x, const TVector &grad, @@ -326,6 +326,22 @@ namespace polysolve::nonlinear solver_info["time_inverting"] = inverting_time / per_iteration; } + void Newton::reset_times() + { + assembly_time = 0; + inverting_time = 0; + } + + void Newton::log_times() const + { + if (assembly_time <= 0 && inverting_time <= 0) + return; // nothing to log + m_logger.debug( + "[{}][{}] assembly: {:.2e}s; linear_solve: {:.2e}s", + fmt::format(fmt::fg(fmt::terminal_color::magenta), "timing"), + name(), assembly_time, inverting_time); + } + // ======================================================================= } // namespace polysolve::nonlinear diff --git a/src/polysolve/nonlinear/descent_strategies/Newton.hpp b/src/polysolve/nonlinear/descent_strategies/Newton.hpp index b50718a..56cc5f1 100644 --- a/src/polysolve/nonlinear/descent_strategies/Newton.hpp +++ b/src/polysolve/nonlinear/descent_strategies/Newton.hpp @@ -47,8 +47,8 @@ namespace polysolve::nonlinear json internal_solver_info = json::array(); const bool is_sparse; - const double m_characteristic_length; - double m_residual_tolerance; + const double characteristic_length; + double residual_tolerance; std::unique_ptr linear_solver; ///< Linear solver used to solve the linear system @@ -71,12 +71,8 @@ namespace polysolve::nonlinear void reset(const int ndof) override; void update_solver_info(json &solver_info, const double per_iteration) override; - - void reset_times() override - { - assembly_time = 0; - inverting_time = 0; - } + void reset_times() override; + void log_times() const override; }; class ProjectedNewton : public Newton diff --git a/src/polysolve/nonlinear/line_search/Armijo.hpp b/src/polysolve/nonlinear/line_search/Armijo.hpp index 14ff81f..c66fc26 100644 --- a/src/polysolve/nonlinear/line_search/Armijo.hpp +++ b/src/polysolve/nonlinear/line_search/Armijo.hpp @@ -13,7 +13,7 @@ namespace polysolve::nonlinear::line_search Armijo(const json ¶ms, spdlog::logger &logger); - virtual std::string name() override { return "Armijo"; } + virtual std::string name() const override { return "Armijo"; } protected: virtual void init_compute_descent_step_size( diff --git a/src/polysolve/nonlinear/line_search/Backtracking.hpp b/src/polysolve/nonlinear/line_search/Backtracking.hpp index 1fdebd4..78663a0 100644 --- a/src/polysolve/nonlinear/line_search/Backtracking.hpp +++ b/src/polysolve/nonlinear/line_search/Backtracking.hpp @@ -13,7 +13,7 @@ namespace polysolve::nonlinear::line_search Backtracking(const json ¶ms, spdlog::logger &logger); - virtual std::string name() override { return "Backtracking"; } + virtual std::string name() const override { return "Backtracking"; } double compute_descent_step_size( const TVector &x, diff --git a/src/polysolve/nonlinear/line_search/LineSearch.cpp b/src/polysolve/nonlinear/line_search/LineSearch.cpp index 1815fbe..1102d3c 100644 --- a/src/polysolve/nonlinear/line_search/LineSearch.cpp +++ b/src/polysolve/nonlinear/line_search/LineSearch.cpp @@ -9,7 +9,7 @@ #include -#include +#include #include @@ -117,9 +117,9 @@ namespace polysolve::nonlinear::line_search } { - POLYSOLVE_SCOPED_STOPWATCH("CCD narrow-phase", ccd_time, m_logger); + POLYSOLVE_SCOPED_STOPWATCH("CCD narrow-phase", narrow_phase_ccd_time, m_logger); m_logger.trace("Performing narrow-phase CCD"); - step_size = compute_collision_free_step_size(x, delta_x, objFunc, step_size); + step_size = compute_max_step_size(x, delta_x, objFunc, step_size); if (std::isnan(step_size)) return NaN; } @@ -213,7 +213,8 @@ namespace polysolve::nonlinear::line_search return step_size; } - double LineSearch::compute_collision_free_step_size( + // change to max_step_size + double LineSearch::compute_max_step_size( const TVector &x, const TVector &delta_x, Problem &objFunc, @@ -242,4 +243,33 @@ namespace polysolve::nonlinear::line_search return step_size; } + + void LineSearch::update_solver_info(json &solver_info, const double per_iteration) + { + solver_info["line_search_iterations"] = iterations(); + solver_info["time_checking_for_nan_inf"] = checking_for_nan_inf_time / per_iteration; + solver_info["time_broad_phase_ccd"] = broad_phase_ccd_time / per_iteration; + solver_info["time_ccd"] = narrow_phase_ccd_time / per_iteration; + solver_info["time_classical_line_search"] = (classical_line_search_time - constraint_set_update_time) / per_iteration; // Remove double counting + solver_info["time_line_search_constraint_set_update"] = constraint_set_update_time / per_iteration; + } + + void LineSearch::reset_times() + { + checking_for_nan_inf_time = 0; + broad_phase_ccd_time = 0; + narrow_phase_ccd_time = 0; + constraint_set_update_time = 0; + classical_line_search_time = 0; + } + + void LineSearch::log_times() const + { + m_logger.debug( + "[{}][{}] constraint_set_update {:.2e}s, checking_for_nan_inf {:.2e}s, " + "broad_phase_ccd {:.2e}s, narrow_phase_ccd {:.2e}s, classical_line_search {:.2e}s", + fmt::format(fmt::fg(fmt::terminal_color::magenta), "timing"), + name(), constraint_set_update_time, checking_for_nan_inf_time, + broad_phase_ccd_time, narrow_phase_ccd_time, classical_line_search_time); + } } // namespace polysolve::nonlinear::line_search diff --git a/src/polysolve/nonlinear/line_search/LineSearch.hpp b/src/polysolve/nonlinear/line_search/LineSearch.hpp index e195955..52e7a4d 100644 --- a/src/polysolve/nonlinear/line_search/LineSearch.hpp +++ b/src/polysolve/nonlinear/line_search/LineSearch.hpp @@ -14,73 +14,77 @@ namespace polysolve::nonlinear::line_search public: using Scalar = typename Problem::Scalar; using TVector = typename Problem::TVector; - + + public: + /// @brief Constructor for creating new LineSearch Object + /// @param params JSON of solver parameters + /// @param m_logger LineSearch(const json ¶ms, spdlog::logger &m_logger); virtual ~LineSearch() = default; + /// @brief + /// @param x Current input vector (n x 1) + /// @param x_delta Current descent direction (n x 1) + /// @param objFunc Objective function + /// @return double line_search( const TVector &x, - const TVector &grad, + const TVector &x_delta, Problem &objFunc); + /// @brief Dispatch function for creating appropriate subclass + /// @param params JSON of solver parameters + /// @param logger + /// @return Pointer to object of the specified subclass static std::shared_ptr create( const json ¶ms, spdlog::logger &logger); + /// @brief Get list of available method names + /// @return Vector of names of available methods static std::vector available_methods(); - void reset_times() - { - checking_for_nan_inf_time = 0; - broad_phase_ccd_time = 0; - ccd_time = 0; - constraint_set_update_time = 0; - classical_line_search_time = 0; - } + virtual std::string name() const = 0; + + void update_solver_info(json &solver_info, const double per_iteration); + void reset_times(); + void log_times() const; void set_is_final_strategy(const bool val) { is_final_strategy = val; } - inline double current_min_step_size() const + double current_min_step_size() const { return is_final_strategy ? min_step_size_final : min_step_size; } - inline int current_max_step_size_iter() const + int current_max_step_size_iter() const { return is_final_strategy ? max_step_size_iter_final : max_step_size_iter; } + int iterations() const { return cur_iter; } + double checking_for_nan_inf_time; double broad_phase_ccd_time; - double ccd_time; + double narrow_phase_ccd_time; double constraint_set_update_time; double classical_line_search_time; double use_grad_norm_tol = -1; - virtual std::string name() = 0; - - inline int iterations() const { return cur_iter; } - - private: - double min_step_size; - int max_step_size_iter; - double min_step_size_final; - int max_step_size_iter_final; - - bool is_final_strategy; - - double default_init_step_size; - protected: - int cur_iter; - spdlog::logger &m_logger; - - double step_ratio; - + /// @brief Compute step size to use during line search + /// @param x Current input (n x 1) + /// @param delta_x Current step direction (n x 1) + /// @param objFunc Problem to be minimized + /// @param use_grad_norm Whether to compare grad norm or energy norm in stopping criteria + /// @param old_energy Previous energy (scalar) + /// @param old_grad Previous gradient (n x 1) + /// @param starting_step_size Initial step size + /// @return Step size to use in line search virtual double compute_descent_step_size( const TVector &x, const TVector &delta_x, @@ -90,17 +94,43 @@ namespace polysolve::nonlinear::line_search const TVector &old_grad, const double starting_step_size) = 0; + spdlog::logger &m_logger; + double step_ratio; + int cur_iter; + private: + /// @brief Compute step size that avoids nan/infinite energy + /// @param x Current input (n x 1) + /// @param delta_x Current step direction (n x 1) + /// @param objFunc Problem to be minimized + /// @param starting_step_size Initial step size + /// @param rate Rate at which to decrease step size if too large + /// @return Nan free step size double compute_nan_free_step_size( const TVector &x, const TVector &delta_x, Problem &objFunc, const double starting_step_size, const double rate); - double compute_collision_free_step_size( + /// @brief Compute maximum valid step size + /// @param x Current input (n x 1) + /// @param delta_x Current step direction (n x 1) + /// @param objFunc Problem to be minimized + /// @param starting_step_size Initial step size + /// @return Maximum valid step size (NaN if it is 0) + double compute_max_step_size( const TVector &x, const TVector &delta_x, Problem &objFunc, const double starting_step_size); + + double min_step_size; + int max_step_size_iter; + double min_step_size_final; + int max_step_size_iter_final; + + bool is_final_strategy; + + double default_init_step_size; }; } // namespace polysolve::nonlinear::line_search diff --git a/src/polysolve/nonlinear/line_search/NoLineSearch.hpp b/src/polysolve/nonlinear/line_search/NoLineSearch.hpp index f7eb101..5f257f2 100644 --- a/src/polysolve/nonlinear/line_search/NoLineSearch.hpp +++ b/src/polysolve/nonlinear/line_search/NoLineSearch.hpp @@ -12,7 +12,7 @@ namespace polysolve::nonlinear::line_search NoLineSearch(const json ¶ms, spdlog::logger &logger); - virtual std::string name() override { return "None"; } + virtual std::string name() const override { return "None"; } protected: double compute_descent_step_size( diff --git a/src/polysolve/nonlinear/line_search/RobustArmijo.hpp b/src/polysolve/nonlinear/line_search/RobustArmijo.hpp index 9b57316..97937fb 100644 --- a/src/polysolve/nonlinear/line_search/RobustArmijo.hpp +++ b/src/polysolve/nonlinear/line_search/RobustArmijo.hpp @@ -14,7 +14,7 @@ namespace polysolve::nonlinear::line_search RobustArmijo(const json ¶ms, spdlog::logger &logger); - virtual std::string name() override { return "RobustArmijo"; } + virtual std::string name() const override { return "RobustArmijo"; } protected: bool criteria(