Merge branch 'main' into add-comments

src/polysolve/linear/Solver.hpp

@@ -86,27 +86,33 @@ namespace polysolve::linear
         // Public interface //
         //////////////////////
 
-        // Set solver parameters
+        /// Set solver parameters
         virtual void set_parameters(const json &params) {}
 
-        // Get info on the last solve step
+        /// Get info on the last solve step
         virtual void get_info(json &params) 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; }
 
+        /// 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 }
         ///
@@ -118,11 +124,6 @@ namespace polysolve::linear
         ///
         virtual void solve(const Ref<const VectorXd> b, Ref<VectorXd> x) = 0;
 
-    public:
-        ///////////
-        // Debug //
-        ///////////
-
         /// @brief Name of the solver type (for debugging purposes)
         virtual std::string name() const { return ""; }
     };

src/polysolve/nonlinear/Criteria.cpp

@@ -2,6 +2,8 @@
 // License: MIT
 #include "Criteria.hpp"
 
+#include <spdlog/fmt/fmt.h>
+
 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 &current)
@@ -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;

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);