update 5/30

.github/workflows/continuous.yml

@@ -42,18 +42,19 @@ jobs:
           libglu1-mesa-dev \
           xorg-dev \
           mpi \
-          ccache
+          ccache \
+          libsuitesparse-dev
           echo 'CACHE_PATH=~/.ccache' >> "$GITHUB_ENV"
 
       - name: Dependencies (macOS)
         if: runner.os == 'macOS'
         run: |
-          brew install ccache open-mpi
+          brew install ccache open-mpi suitesparse
           echo 'CACHE_PATH=~/Library/Caches/ccache' >> "$GITHUB_ENV"
 
       - name: Cache Build
         id: cache-build
-        uses: actions/cache@v4.0.2
+        uses: actions/cache@v4
         with:
           path: ${{ env.CACHE_PATH }}
           key: ${{ runner.os }}-${{ matrix.config }}-cache-${{ github.sha }}
@@ -105,7 +106,7 @@ jobs:
 
       - name: Cache build
         id: cache-build
-        uses: actions/cache@v4.0.2
+        uses: actions/cache@v4
         with:
           path: ${{ env.appdata }}\Mozilla\sccache
           key: ${{ runner.os }}-${{ matrix.config }}-cache-${{ github.sha }}

.github/workflows/coverage.yml

@@ -41,7 +41,7 @@ jobs:
 
       - name: Cache Build
         id: cache-build
-        uses: actions/cache@v4.0.2
+        uses: actions/cache@v4
         with:
           path: ${{ env.CACHE_PATH }}
           key: ${{ runner.os }}-Release-${{ matrix.threading }}-cache-${{ github.sha }}

CMakeLists.txt

@@ -10,6 +10,12 @@ endif()
 set(REQUIRED_CMAKE_VERSION "3.18.0")
 if(POLYSOLVE_TOPLEVEL_PROJECT)
     cmake_minimum_required(VERSION ${REQUIRED_CMAKE_VERSION})
+    SET(CMAKE_POLICY_VERSION_MINIMUM ${REQUIRED_CMAKE_VERSION})
+    if(CMAKE_VERSION VERSION_GREATER_EQUAL "4.0.0")
+        message(WARNING "CMake version is ${CMAKE_VERSION}, setting stuff for dependencies.")
+        SET(CMAKE_POLICY_VERSION_MINIMUM ${REQUIRED_CMAKE_VERSION})
+        set(POLYSOLVE_WITH_AMGCL OFF CACHE BOOL "Use AMGCL for solving linear systems")
+    endif()
 else()
     # Don't use cmake_minimum_required here to avoid implicitly overriding parent policies
     if(${CMAKE_VERSION} VERSION_LESS ${REQUIRED_CMAKE_VERSION})
@@ -74,6 +80,7 @@ option(POLYSOLVE_WITH_ACCELERATE    "Enable Apple Accelerate" ${POLYSOLVE_ON_APP
 option(POLYSOLVE_WITH_CHOLMOD       "Enable Cholmod library"                             ON)
 option(POLYSOLVE_WITH_UMFPACK       "Enable UmfPack library"                             ON)
 option(POLYSOLVE_WITH_SUPERLU       "Enable SuperLU library"                             ON)
+option(POLYSOLVE_WITH_SPQR          "Enable SPQR library"                                ON)
 option(POLYSOLVE_WITH_MKL           "Enable MKL library"  ${POLYSOLVE_NOT_ON_APPLE_SILICON})
 option(POLYSOLVE_WITH_CUSOLVER      "Enable cuSOLVER library"                           OFF)
 option(POLYSOLVE_WITH_PARDISO       "Enable Pardiso library"                            OFF)
@@ -269,6 +276,17 @@ if(POLYSOLVE_WITH_SUPERLU)
     endif()
 endif()
 
+# SuperLU solver
+if(POLYSOLVE_WITH_SPQR)
+    include(spqr)
+    if(TARGET SuiteSparse::SPQR)
+        target_link_libraries(polysolve_linear PRIVATE SuiteSparse::SPQR)
+        target_compile_definitions(polysolve_linear PUBLIC POLYSOLVE_WITH_SPQR)
+    else()
+        message(WARNING "SPQR Not found, solver will not be available.")
+    endif()
+endif()
+
 # AMGCL solver
 if(POLYSOLVE_WITH_AMGCL)
     include(amgcl)

cmake/recipes/boost.cmake

@@ -20,7 +20,9 @@ set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fPIC")
 set(OLD_CMAKE_POSITION_INDEPENDENT_CODE ${CMAKE_POSITION_INDEPENDENT_CODE})
 set(CMAKE_POSITION_INDEPENDENT_CODE ON)
 
-set(BOOST_URL "https://boostorg.jfrog.io/artifactory/main/release/1.83.0/source/boost_1_83_0.tar.bz2" CACHE STRING "Boost download URL")
+set(BOOST_URL 
+    "https://archives.boost.io/release/1.83.0/source/boost_1_83_0.tar.bz2"
+    CACHE STRING "Boost download URL")
 set(BOOST_URL_SHA256 "6478edfe2f3305127cffe8caf73ea0176c53769f4bf1585be237eb30798c3b8e" CACHE STRING "Boost download URL SHA256 checksum")
 
 include(CPM)

cmake/recipes/spqr.cmake

@@ -0,0 +1,11 @@
+# SPQR solver
+
+if(TARGET SparseSuite::SPQR)
+    return()
+endif()
+
+message(STATUS "Third-party: creating targets 'SuiteSparse::SPQR'")
+
+# We do not have a build recipe for this, so find it as a system installed library.
+find_package(SPQR)
+

src/polysolve/Utils.cpp

@@ -1,6 +1,10 @@
 #include "Utils.hpp"
 
+#if defined(SPDLOG_FMT_EXTERNAL)
+#include <fmt/color.h>
+#else
 #include <spdlog/fmt/bundled/color.h>
+#endif
 
 namespace polysolve
 {
@@ -49,7 +53,7 @@ namespace polysolve
 
     void StopWatch::log_msg()
     {
-        const static std::string log_fmt_text =
+        const static auto log_fmt_text =
             fmt::format("[{}] {{}} {{:.3g}}s", fmt::format(fmt::fg(fmt::terminal_color::magenta), "timing"));
 
         if (!m_name.empty())
@@ -60,7 +64,7 @@ namespace polysolve
 
     void log_and_throw_error(spdlog::logger &logger, const std::string &msg)
     {
-        logger.error(msg);
+        logger.error("{}", msg);
         throw std::runtime_error(msg);
     }
 
@@ -79,4 +83,4 @@ namespace polysolve
         return json[name];
     }
 
-} // namespace polysolve
+} // namespace polysolve

src/polysolve/Utils.hpp

@@ -55,7 +55,7 @@ namespace polysolve
     template <typename... Args>
     [[noreturn]] void log_and_throw_error(spdlog::logger &logger, const std::string &msg, const Args &...args)
     {
-        log_and_throw_error(logger, fmt::format(msg, args...));
+        log_and_throw_error(logger, fmt::format(fmt::runtime(msg), args...));
     }
 
     Eigen::SparseMatrix<double> sparse_identity(int rows, int cols);

src/polysolve/linear/Solver.cpp

@@ -11,6 +11,12 @@
 #include <fstream>
 
 // -----------------------------------------------------------------------------
+// 
+// Subsequent macros assume a single template parameter and SparseQR fails due to requiring 2 parameters. this is because the OrderingType is not filled in.
+// SparseLU has a default declaration of _OrderingType to use COLAMDOrdering but SparseQR doesn't - so this just mimics that behavior. If Eigen adds such a default in the future this line will need to be guarded to avoid multiple defaults
+namespace Eigen {
+template <typename _MatrixType, typename _OrderingType = COLAMDOrdering<typename _MatrixType::StorageIndex> > class SparseQR;
+}
 #include <Eigen/Sparse>
 #ifdef POLYSOLVE_WITH_ACCELERATE
 #include <Eigen/AccelerateSupport>
@@ -21,6 +27,24 @@
 #ifdef POLYSOLVE_WITH_UMFPACK
 #include <Eigen/UmfPackSupport>
 #endif
+#ifdef POLYSOLVE_WITH_SPQR
+#include <Eigen/SPQRSupport>
+namespace polysolve::linear {
+    template <>
+    void EigenDirect<Eigen::SPQR<StiffnessMatrix>>::analyze_pattern(const StiffnessMatrix& A, const int precond_num) {
+        m_Solver.compute(A);
+    }
+    template <>
+    void EigenDirect<Eigen::SPQR<StiffnessMatrix>>::factorize(const StiffnessMatrix &A)
+    {
+        m_Solver.compute(A);
+        if (m_Solver.info() == Eigen::NumericalIssue)
+        {
+            throw std::runtime_error("[EigenDirect] NumericalIssue encountered.");
+        }
+    }
+}
+#endif
 #ifdef POLYSOLVE_WITH_SUPERLU
 #include <Eigen/SuperLUSupport>
 #endif
@@ -293,6 +317,10 @@ namespace polysolve::linear
         else if (solver == "Eigen::SparseLU")
         {
             RETURN_DIRECT_SOLVER_PTR(SparseLU, "Eigen::SparseLU");
+        }
+        else if (solver == "Eigen::SparseQR")
+        {
+            RETURN_DIRECT_SOLVER_PTR(SparseQR, "Eigen::SparseQR");
 #ifdef POLYSOLVE_WITH_ACCELERATE
         }
         else if (solver == "Eigen::AccelerateLLT")
@@ -335,6 +363,12 @@ namespace polysolve::linear
         {
             RETURN_DIRECT_SOLVER_PTR(SuperLU, "Eigen::SuperLU");
 #endif
+#ifdef POLYSOLVE_WITH_SPQR
+        }
+        else if (solver == "Eigen::SPQR")
+        {
+            RETURN_DIRECT_SOLVER_PTR(SPQR, "Eigen::SPQR");
+#endif
 #ifdef POLYSOLVE_WITH_MKL
         }
         else if (solver == "Eigen::PardisoLLT")
@@ -465,6 +499,7 @@ namespace polysolve::linear
         return {{
             "Eigen::SimplicialLDLT",
             "Eigen::SparseLU",
+            "Eigen::SparseQR",
 #ifdef POLYSOLVE_WITH_ACCELERATE
             "Eigen::AccelerateLLT",
             "Eigen::AccelerateLDLT",
@@ -481,6 +516,9 @@ namespace polysolve::linear
 #ifdef POLYSOLVE_WITH_SUPERLU
             "Eigen::SuperLU",
 #endif
+#ifdef POLYSOLVE_WITH_SPQR
+            "Eigen::SPQR",
+#endif
 #ifdef POLYSOLVE_WITH_MKL
             "Eigen::PardisoLLT",
             "Eigen::PardisoLDLT",

src/polysolve/nonlinear/Criteria.cpp

@@ -29,7 +29,10 @@ namespace polysolve::nonlinear
 
     void Criteria::print(std::ostream &os) const
     {
-        os << fmt::format(
+        os << print_message();
+    }
+    std::string Criteria::print_message() const {
+        return fmt::format(
             "iters={:d} Δf={:g} ‖∇f‖={:g} ‖Δx‖={:g} Δx⋅∇f(x)={:g}",
             iterations, fDelta, gradNorm, xDelta, xDeltaDotGrad);
     }
@@ -61,47 +64,39 @@ namespace polysolve::nonlinear
         return Status::Continue;
     }
 
-    std::ostream &operator<<(std::ostream &os, const Status &s)
-    {
+    std::string_view status_message(Status s) {
         switch (s)
         {
         case Status::NotStarted:
-            os << "Solver not started";
-            break;
+            return "Solver not started";
         case Status::Continue:
-            os << "Convergence criteria not reached";
-            break;
+            return "Convergence criteria not reached";
         case Status::IterationLimit:
-            os << "Iteration limit reached";
-            break;
+            return "Iteration limit reached";
         case Status::XDeltaTolerance:
-            os << "Change in parameter vector too small";
-            break;
+            return "Change in parameter vector too small";
         case Status::FDeltaTolerance:
-            os << "Change in cost function value too small";
-            break;
+            return "Change in cost function value too small";
         case Status::GradNormTolerance:
-            os << "Gradient vector norm too small";
-            break;
+            return "Gradient vector norm too small";
         case Status::ObjectiveCustomStop:
-            os << "Objective function specified to stop";
-            break;
+            return "Objective function specified to stop";
         case Status::NanEncountered:
-            os << "Objective or gradient function returned NaN";
-            break;
+            return "Objective or gradient function returned NaN";
         case Status::NotDescentDirection:
-            os << "Search direction not a descent direction";
-            break;
+            return "Search direction not a descent direction";
         case Status::LineSearchFailed:
-            os << "Line search failed";
-            break;
+            return "Line search failed";
         case Status::UpdateDirectionFailed:
-            os << "Update direction could not be computed";
-            break;
+            return "Update direction could not be computed";
         default:
-            os << "Unknown status";
-            break;
+            return "Unknown status";
         }
+    }
+
+    std::ostream &operator<<(std::ostream &os, const Status &s)
+    {
+        os << status_message(s);
         return os;
     }
 
@@ -110,4 +105,4 @@ namespace polysolve::nonlinear
         c.print(os);
         return os;
     }
-} // namespace polysolve::nonlinear
+} // namespace polysolve::nonlinear

src/polysolve/nonlinear/Criteria.hpp

@@ -2,6 +2,7 @@
 
 #include <cstddef>
 #include <iostream>
+#include <string_view>
 
 namespace polysolve::nonlinear
 {
@@ -43,11 +44,15 @@ namespace polysolve::nonlinear
         void reset();
 
         void print(std::ostream &os) const;
+        std::string print_message() const;
     };
 
     Status checkConvergence(const Criteria &stop, const Criteria &current);
 
+    std::string_view status_message(Status s);
+    std::string criteria_message(const Criteria& s);
+
     std::ostream &operator<<(std::ostream &os, const Status &s);
 
     std::ostream &operator<<(std::ostream &os, const Criteria &c);
-} // namespace polysolve::nonlinear
+} // namespace polysolve::nonlinear