Adapt CholmodSolver to (almost) arbitrary matrix types

This is a follow-up change from the dune-istl Cholmod changes in
https://gitlab.dune-project.org/core/dune-istl/-/merge_requests/424

CHANGELOG.md

@@ -4,6 +4,8 @@
 
 - CholmodSolver: Added `errorCode_` member to report errors and warnings during the matrix decomposition
 
+- CholmodSolver: `CholmodSolver` can be used with any blocked matrix/vector type supported by `flatMatrixForEach` and `flatVectorForEach` in dune-istl.
+
 ## Deprecations and removals
 
 - The file `blockgsstep.hh` has been removed after four years of deprecation.  Use the replacement

dune/solvers/solvers/cholmodsolver.hh

@@ -13,11 +13,13 @@
 
 #include <dune/istl/solver.hh>
 #include <dune/istl/cholmod.hh>
+#include <dune/istl/foreach.hh>
 #include <dune/istl/matrixindexset.hh>
 #include <dune/istl/io.hh>
 
 #include <dune/solvers/common/canignore.hh>
 #include <dune/solvers/solvers/linearsolver.hh>
+#include <dune/solvers/common/defaultbitvector.hh>
 
 namespace Dune
 {
@@ -33,17 +35,12 @@ namespace Solvers
 *  This class wraps the dune-istl CHOLMOD solver for dune-solvers to provide
 *  a Solvers::LinearSolver interface.
 */
-template <class MatrixType, class VectorType>
+template <class MatrixType, class VectorType, class BitVectorType = DefaultBitVector_t<VectorType>>
 class CholmodSolver
-: public LinearSolver<MatrixType,VectorType>, public CanIgnore<Dune::BitSetVector<VectorType::block_type::dimension> >
+: public LinearSolver<MatrixType,VectorType>, public CanIgnore<BitVectorType>
 {
-  enum {blocksize = VectorType::block_type::dimension};
-
   using field_type = typename VectorType::field_type;
 
-  using  MatrixBlock = Dune::FieldMatrix<field_type, blocksize, blocksize> ;
-  using  VectorBlock = Dune::FieldVector<field_type, blocksize>;
-
 public:
 
   /** \brief Default constructor */
@@ -110,46 +107,45 @@ public:
       // get the error code from Cholmod in case something happened
       errorCode_ = solver.cholmodCommonObject().status;
 
-      auto modifiedRhs = *rhs_;
+      // create flat vectors
+      using T = typename VectorType::field_type;
+      std::size_t flatSize = flatVectorForEach(ignore, [](...){});
+      std::vector<bool> flatIgnore(flatSize);
+      std::vector<T> flatX(flatSize), flatRhsModifier(flatSize);
 
-      // pull all ignored columns to the rhs
-      for(auto rowIt = matrix_->begin(); rowIt != matrix_->end(); rowIt++)
+      // define a lambda that copies the entries of a blocked vector into a flat one
+      auto copyToFlat = [&](auto&& blocked, auto&& flat)
       {
-        const auto row = rowIt.index();
-
-        for(auto ii=0; ii<blocksize; ii++)
+        flatVectorForEach(blocked, [&](auto&& entry, auto&& offset)
         {
-          if( ignore[row][ii] )
-            continue;
-
-          for(auto colIt = rowIt->begin(); colIt != rowIt->end(); colIt++)
-          {
-            const auto col = colIt.index();
-            for(auto jj=0; jj<blocksize; jj++)
-            {
-              if( ignore[col][jj] )
-              {
-                // we assume values only in the upper part
-                if ( row < col )
-                {
-                  // upper part -> use canonical indices
-                  modifiedRhs[row][ii] -= (*matrix_)[row][col][ii][jj] * (*x_)[col][jj];
-                }
-                else if (row == col and ii > jj )
-                {
-                  // diagonal block but there in the lower part -> swap ii <-> jj
-                  modifiedRhs[row][ii] -= (*matrix_)[row][col][jj][ii] * (*x_)[col][jj];
-                }
-                else
-                {
-                  // lower part -> swap col <-> row and ii <-> jj
-                  modifiedRhs[row][ii] -= (*matrix_)[col][row][jj][ii] * (*x_)[col][jj];
-                }
-              }
-            }
-          }
-        }
-      }
+          flat[offset] = entry;
+        });
+      };
+
+      // copy x and the ignore field for the modification of the rhs
+      copyToFlat(ignore,flatIgnore);
+      copyToFlat(*x_,flatX);
+
+      flatMatrixForEach( *matrix_, [&]( auto&& entry, auto&& rowOffset, auto&& colOffset){
+
+        // we assume entries only in the upper part and do nothing on the diagonal
+        if ( rowOffset >= colOffset )
+          return;
+
+        // upper part: if either col or row is ignored: move the entry or the symmetric duplicate to the rhs
+        if ( flatIgnore[colOffset] and not flatIgnore[rowOffset] )
+            flatRhsModifier[rowOffset] += entry * flatX[colOffset];
+
+        else if ( not flatIgnore[colOffset] and flatIgnore[rowOffset] )
+            flatRhsModifier[colOffset] += entry * flatX[rowOffset];
+      });
+
+      // update the rhs
+      auto modifiedRhs = *rhs_;
+      flatVectorForEach(modifiedRhs, [&](auto&& entry, auto&& offset){
+        if ( not flatIgnore[offset] )
+          entry -= flatRhsModifier[offset];
+      });
 
       // Solve the modified system
       solver.apply(*x_, modifiedRhs, statistics);