Merge pull request #84 from mtao/mtao/add_spqr

Add QR

Author: teseoch

.github/workflows/continuous.yml

@@ -42,13 +42,14 @@ 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

CMakeLists.txt

@@ -74,6 +74,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 +270,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/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/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",