Allow to use CHOLMOD with 'SuiteSparse_long' integer type

CHOLMOD supports this, it just wasn't exposed in the interface.
It is needed for large problems.

CHANGELOG.md

@@ -53,6 +53,11 @@ SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
   used by `CHOLMOD` itself to store the Cholesky factor.  This can be used to
   use the more advanced features of `CHOLMOD`.
 
+- The `Cholmod` class now has an additional template parameter `Index`, which specifies the
+  type used by `CHOLMOD` for integers.  The types `int` (the default) and `SuiteSparse_long`
+  (which is usually `long int`) are supported.  Use `SuiteSparse_long` if your problems
+  are too big for `int`.
+
 - You can now multiply objects of type `ScaledIdentityMatrix` by scalars
   using `operator*`.
 

dune/istl/cholmod.hh

@@ -66,16 +66,157 @@ namespace Impl{
     });
   }
 
+  // wrapper class for C function calls to CHOLMOD itself.
+  // The CHOLMOD API has different functions for different index types.
+  template <class Index>
+  struct CholmodMethodChooser;
+
+  // specialization using 'int' to store indices
+  template <>
+  struct CholmodMethodChooser<int>
+  {
+    [[nodiscard]]
+    static cholmod_dense* allocate_dense(size_t nrow, size_t ncol, size_t d, int xtype, cholmod_common *c)
+    {
+      return ::cholmod_allocate_dense(nrow,ncol,d,xtype,c);
+    }
+
+    [[nodiscard]]
+    static cholmod_sparse* allocate_sparse(size_t nrow, size_t ncol, size_t nzmax, int sorted, int packed, int stype, int xtype, cholmod_common *c)
+    {
+      return ::cholmod_allocate_sparse(nrow,ncol,nzmax,sorted,packed,stype,xtype,c);
+    }
+
+    [[nodiscard]]
+    static cholmod_factor* analyze(cholmod_sparse *A, cholmod_common *c)
+    {
+      return ::cholmod_analyze(A,c);
+    }
+
+    static int defaults(cholmod_common *c)
+    {
+      return ::cholmod_defaults(c);
+    }
+
+    static int factorize(cholmod_sparse *A, cholmod_factor *L, cholmod_common *c)
+    {
+      return ::cholmod_factorize(A,L,c);
+    }
+
+    static int finish(cholmod_common *c)
+    {
+      return ::cholmod_finish(c);
+    }
+
+    static int free_dense (cholmod_dense **X, cholmod_common *c)
+    {
+      return ::cholmod_free_dense(X,c);
+    }
+
+    static int free_factor(cholmod_factor **L, cholmod_common *c)
+    {
+      return ::cholmod_free_factor(L,c);
+    }
+
+    static int free_sparse(cholmod_sparse **A, cholmod_common *c)
+    {
+      return ::cholmod_free_sparse(A,c);
+    }
+
+    [[nodiscard]]
+    static cholmod_dense* solve(int sys, cholmod_factor *L, cholmod_dense *B, cholmod_common *c)
+    {
+      return ::cholmod_solve(sys,L,B,c);
+    }
+
+    static int start(cholmod_common *c)
+    {
+      return ::cholmod_start(c);
+    }
+  };
+
+  // specialization using 'SuiteSparse_long' to store indices
+  template <>
+  struct CholmodMethodChooser<SuiteSparse_long>
+  {
+    [[nodiscard]]
+    static cholmod_dense* allocate_dense(size_t nrow, size_t ncol, size_t d, int xtype, cholmod_common *c)
+    {
+      return ::cholmod_l_allocate_dense(nrow,ncol,d,xtype,c);
+    }
+
+    [[nodiscard]]
+    static cholmod_sparse* allocate_sparse(size_t nrow, size_t ncol, size_t nzmax, int sorted, int packed, int stype, int xtype, cholmod_common *c)
+    {
+      return ::cholmod_l_allocate_sparse(nrow,ncol,nzmax,sorted,packed,stype,xtype,c);
+    }
+
+    [[nodiscard]]
+    static cholmod_factor* analyze(cholmod_sparse *A, cholmod_common *c)
+    {
+      return ::cholmod_l_analyze(A,c);
+    }
+
+    static int defaults(cholmod_common *c)
+    {
+      return ::cholmod_l_defaults(c);
+    }
+
+    static int factorize(cholmod_sparse *A, cholmod_factor *L, cholmod_common *c)
+    {
+      return ::cholmod_l_factorize(A,L,c);
+    }
+
+    static int finish(cholmod_common *c)
+    {
+      return ::cholmod_l_finish(c);
+    }
+
+    static int free_dense (cholmod_dense **X, cholmod_common *c)
+    {
+      return ::cholmod_l_free_dense(X,c);
+    }
+
+    static int free_factor (cholmod_factor **L, cholmod_common *c)
+    {
+      return ::cholmod_l_free_factor(L,c);
+    }
+
+    static int free_sparse(cholmod_sparse **A, cholmod_common *c)
+    {
+      return ::cholmod_l_free_sparse(A,c);
+    }
+
+    [[nodiscard]]
+    static cholmod_dense* solve(int sys, cholmod_factor *L, cholmod_dense *B, cholmod_common *c)
+    {
+      return ::cholmod_l_solve(sys,L,B,c);
+    }
+
+    static int start(cholmod_common *c)
+    {
+      return ::cholmod_l_start(c);
+    }
+  };
 
 } //namespace Impl
 
 /** @brief Dune wrapper for SuiteSparse/CHOLMOD solver
   *
   * This class implements an InverseOperator between Vector types
+  *
+  * \tparam Vector Data type for solution and right-hand-side vectors
+  * \tparam Index Type used by CHOLMOD for indices.  Must be either 'int' or 'SuiteSparse_long'
+  *   (which is usually `long int`).
   */
-template<class Vector>
+template<class Vector, class Index=int>
 class Cholmod : public InverseOperator<Vector, Vector>
 {
+  static_assert(std::is_same_v<Index,int> || std::is_same_v<Index,SuiteSparse_long>,
+                "Index type must be either 'int' or 'SuiteSparse_long'!");
+
+  using CholmodMethod = Impl::CholmodMethodChooser<Index>;
+
 public:
 
   /** @brief Default constructor
@@ -85,7 +226,7 @@ public:
    */
   Cholmod()
   {
-    cholmod_start(&c_);
+    CholmodMethod::start(&c_);
   }
 
   /** @brief Destructor
@@ -96,8 +237,8 @@ public:
   ~Cholmod()
   {
     if (L_)
-      cholmod_free_factor(&L_, &c_);
-    cholmod_finish(&c_);
+      CholmodMethod::free_factor(&L_, &c_);
+    CholmodMethod::finish(&c_);
   }
 
   // forbid copying to avoid freeing memory twice
@@ -144,13 +285,14 @@ public:
     });
 
       // create a cholmod dense object
-    auto b3 = make_cholmod_dense(cholmod_allocate_dense(L_->n, 1, L_->n, CHOLMOD_REAL, &c_), &c_);
+    auto b3 = make_cholmod_dense(CholmodMethod::allocate_dense(L_->n, 1, L_->n, CHOLMOD_REAL, &c_), &c_);
+
     // cast because void-ptr
     auto b4 = static_cast<double*>(b3->x);
     std::copy(b2.get(), b2.get() + L_->n, b4);
 
     // solve for a cholmod x object
-    auto x3 = make_cholmod_dense(cholmod_solve(CHOLMOD_A, L_, b3.get(), &c_), &c_);
+    auto x3 = make_cholmod_dense(CholmodMethod::solve(CHOLMOD_A, L_, b3.get(), &c_), &c_);
     // cast because void-ptr
     auto xp = static_cast<double*>(x3->x);
 
@@ -232,22 +374,22 @@ public:
     * by DUNE.  So we can just store Mᵀ instead of M (as M = Mᵀ).
     */
     const auto deleter = [c = &this->c_](auto* p) {
-      cholmod_free_sparse(&p, c);
+      CholmodMethod::free_sparse(&p, c);
     };
     auto M = std::unique_ptr<cholmod_sparse, decltype(deleter)>(
-      cholmod_allocate_sparse(N,             // # rows
-                              N,             // # cols
-                              nonZeros,      // # of nonzeroes
-                              1,             // indices are sorted ( 1 = true)
-                              1,             // matrix is "packed" ( 1 = true)
-                              -1,            // stype of matrix ( -1 = consider the lower part only )
-                              CHOLMOD_REAL,  // xtype of matrix ( CHOLMOD_REAL = single array, no complex numbers)
-                              &c_            // cholmod_common ptr
-                             ), deleter);
+      CholmodMethod::allocate_sparse(N,             // # rows
+                                     N,             // # cols
+                                     nonZeros,      // # of nonzeroes
+                                     1,             // indices are sorted ( 1 = true)
+                                     1,             // matrix is "packed" ( 1 = true)
+                                     -1,            // stype of matrix ( -1 = consider the lower part only )
+                                     CHOLMOD_REAL,  // xtype of matrix ( CHOLMOD_REAL = single array, no complex numbers)
+                                     &c_            // cholmod_common ptr
+      ), deleter);
 
     // copy the data of BCRS matrix to Cholmod Sparse matrix
-    int* Ap = static_cast<int*>(M->p);
-    int* Ai = static_cast<int*>(M->i);
+    Index* Ap = static_cast<Index*>(M->p);
+    Index* Ai = static_cast<Index*>(M->i);
     double* Ax = static_cast<double*>(M->x);
 
 
@@ -318,10 +460,10 @@ public:
     });
 
     // Now analyse the pattern and optimal row order
-    L_ = cholmod_analyze(M.get(), &c_);
+    L_ = CholmodMethod::analyze(M.get(), &c_);
 
     // Do the factorization (this may take some time)
-    cholmod_factorize(M.get(), L_, &c_);
+    CholmodMethod::factorize(M.get(), L_, &c_);
   }
 
   virtual SolverCategory::Category category() const
@@ -332,7 +474,7 @@ public:
   /** \brief return a reference to the CHOLMOD common object for advanced option settings
    *
    *  The CHOLMOD common object stores all parameters and options for the solver to run
-   *  and can be modified in several ways, see CHOLMOD Userguide for further information
+   *  and can be modified in several ways, see CHOLMOD Userguide for further information.
    */
   cholmod_common& cholmodCommonObject()
   {
@@ -365,7 +507,7 @@ private:
   auto make_cholmod_dense(cholmod_dense* x, cholmod_common* c)
   {
     const auto deleter = [c](auto* p) {
-      cholmod_free_dense(&p, c);
+      CholmodMethod::free_dense(&p, c);
     };
     return std::unique_ptr<cholmod_dense, decltype(deleter)>(x, deleter);
   }

dune/istl/test/cholmodtest.cc

@@ -26,12 +26,15 @@
 
 using namespace Dune;
 
-int main(int argc, char** argv)
+/** \brief Do the actual tests
+ *
+ * \tparam Index Type used by CHOLMOD for indices (int or long int)
+ */
+template <typename Index>
+bool test()
 {
+  bool passed = true;
 #if HAVE_SUITESPARSE_UMFPACK
-  try
-  {
-
     int N = 30; // number of nodes
     const int bs = 2; // block size
 
@@ -55,7 +58,10 @@ int main(int argc, char** argv)
     A.mmv(x,b);
 
     if ( b.two_norm() > 1e-9 )
+    {
       std::cerr << " Error in CHOLMOD, residual is too large: " << b.two_norm() << std::endl;
+      passed = false;
+    }
 
     x = 0;
     b = 1;
@@ -75,7 +81,10 @@ int main(int argc, char** argv)
 
     // check that x[12][0] is untouched
     if ( std::abs(x[12][0] - 123) > 1e-15 )
+    {
       std::cerr << " Error in CHOLMOD, x was NOT ignored!"<< std::endl;
+      passed = false;
+    }
 
     // reset the x value
     x[12][0] = 0;
@@ -85,8 +94,10 @@ int main(int argc, char** argv)
 
     // check that error is completely caused by this entry
     if ( std::abs( b.two_norm() - std::abs(b_12_0) ) > 1e-15 )
+    {
       std::cerr << " Error in CHOLMOD, b was NOT ignored correctly: " << std::abs( b.two_norm() - std::abs(b_12_0) ) << std::endl;
-
+      passed = false;
+    }
 
     using BCRS = BCRSMatrix<FieldMatrix<double,bs,bs>>;
     using MTRow = MultiTypeBlockVector<BCRS,BCRS>;
@@ -125,7 +136,10 @@ int main(int argc, char** argv)
     A2.mmv(x2,b2);
 
     if ( b2.two_norm() > 1e-9 )
+    {
       std::cerr << " Error in CHOLMOD, residual is too large: " << b2.two_norm() << std::endl;
+      passed = false;
+    }
 
 
     x2 = 0;
@@ -147,7 +161,10 @@ int main(int argc, char** argv)
 
     // check that x was ignored
     if ( std::abs( x2[_0][12][0] - 123 ) > 1e-15 )
+    {
       std::cerr << " Error in CHOLMOD, x was NOT ignored correctly: " << std::abs( x2[_0][12][0] - 123 ) << std::endl;
+      passed = false;
+    }
 
     x2[_0][12][0] = 0;
 
@@ -158,22 +175,36 @@ int main(int argc, char** argv)
 
     // check that error is completely caused by this entry
     if ( std::abs( b2.two_norm() - std::abs(b_0_12_0) ) > 1e-15 )
+    {
       std::cerr << " Error in CHOLMOD, b was NOT ignored correctly: " << std::abs( b2.two_norm() - std::abs(b_0_12_0) ) << std::endl;
+      passed = false;
+    }
 
-
-  }
-  catch (std::exception &e)
-  {
-    std::cout << "ERROR: " << e.what() << std::endl;
-    return 1;
-  }
-  catch (...)
-  {
-    std::cerr << "Dune reported an unknown error." << std::endl;
-    exit(1);
-  }
 #else // HAVE_SUITESPARSE_UMFPACK
   std::cerr << "You need SuiteSparse to run the CHOLMOD test." << std::endl;
-  return 42;
+  passed = false;
 #endif // HAVE_SUITESPARSE_UMFPACK
+  return passed;
+}
+
+int main(int argc, char** argv) try
+{
+  bool passed = true;
+
+  std::cout << "testing int" << std::endl;
+  passed &= test<int>();
+  std::cout << "testing SuiteSparse_long (i.e., long int)" << std::endl;
+  passed &= test<SuiteSparse_long>();
+
+  return (passed) ? 0 : 1;
+}
+catch (std::exception &e)
+{
+  std::cerr << "ERROR: " << e.what() << std::endl;
+  return 1;
+}
+catch (...)
+{
+  std::cerr << "Dune reported an unknown error." << std::endl;
+  exit(1);
 }