Merge branch 'umfpacksolver-with-scalar-matrix-entries' into 'master'

UMFPackSolver with scalar matrix entries

See merge request !51

CHANGELOG.md

 # Master (will become release 2.8)
 
+- `UMFPackSolver` can now handle matrices and vectors with scalar entries.
+
 - CholmodSolver: Added `errorCode_` member to report errors and warnings during the matrix decomposition
 
 ## Deprecations and removals

dune/solvers/common/defaultbitvector.hh

@@ -18,11 +18,18 @@ namespace Solvers {
 
 namespace Imp {
 
+/** \brief Construct a bitvector type with a nesting that matches the one of a given vector
+ *
+ * The default implementation implements the case where the vector is a number type
+ * (to end recursion), and returns 'char' in such cases.
+ *
+ * \tparam Vector The vector types whose nesting shall be matched
+ */
 template<class Vector>
 struct DefaultBitVector
 {
-  static_assert(AlwaysFalse<Vector>::value, "No DefaultBitVector known for this type.");
-  using type = void;
+  static_assert(IsNumber<Vector>::value, "No DefaultBitVector known for this type.");
+  using type = char;
 };
 
 template<class T, int i>

dune/solvers/norms/twonorm.hh

@@ -50,7 +50,7 @@ class TwoNorm : public Norm<V>
                 typename VectorType::block_type block = f1[i];
                 block -= f2[i];
 
-                r += block.two_norm2();
+                r += Dune::Impl::asVector(block).two_norm2();
             }
 
             return std::sqrt(r);

dune/solvers/solvers/umfpacksolver.hh

@@ -88,11 +88,17 @@ public:
       for (size_t i=0; i<matrix_->N(); i++)
       {
         auto const &ignore = this->ignore()[i];
+        if constexpr (std::is_convertible<decltype(ignore), bool>::value)
+        {
+          if (!ignore)
+            nonIgnoreRows.insert(i);
+        } else {
           if (ignore.none())
             nonIgnoreRows.insert(i);
           else if (not ignore.all())
             DUNE_THROW(Dune::NotImplemented, "Individual blocks must be either ignored completely, or not at all");
         }
+      }
 
       // Construct the solver
       Dune::InverseOperatorResult statistics;
@@ -111,15 +117,25 @@ public:
       shortRowCount = 0;
       for (size_t i=0; i<matrix_->N(); i++)
       {
+        if constexpr (std::is_convertible<decltype(this->ignore()[i]), const bool>::value) {
+          if (this->ignore()[i])
+            continue;
+        } else {
           if (this->ignore()[i].all())
             continue;
+        }
 
         auto cIt    = (*matrix_)[i].begin();
         auto cEndIt = (*matrix_)[i].end();
 
         for (; cIt!=cEndIt; ++cIt)
+          if constexpr (std::is_convertible<decltype(this->ignore()[cIt.index()]), const bool>::value) {
+            if (this->ignore()[cIt.index()])
+              Dune::Impl::asMatrix(*cIt).mmv((*x_)[cIt.index()], shortRhs[shortRowCount]);
+          } else {
             if (this->ignore()[cIt.index()].all())
-            cIt->mmv((*x_)[cIt.index()], shortRhs[shortRowCount]);
+              Dune::Impl::asMatrix(*cIt).mmv((*x_)[cIt.index()], shortRhs[shortRowCount]);
+          }
 
         shortRowCount++;
       }

dune/solvers/test/common.hh

@@ -18,7 +18,7 @@
 
 #include <dune/grid/utility/structuredgridfactory.hh>
 
-#if HAVE_UG
+#if HAVE_DUNE_UGGRID
 #include <dune/grid/uggrid.hh>
 #endif
 
@@ -26,6 +26,9 @@
 #include <dune/alugrid/grid.hh>
 #endif
 
+#include <dune/matrix-vector/addtodiagonal.hh>
+
+#include <dune/solvers/common/defaultbitvector.hh>
 #include <dune/solvers/norms/energynorm.hh>
 #include <dune/solvers/norms/twonorm.hh>
 
@@ -64,13 +67,6 @@ void constructPQ1Pattern(const GridView& gridView, Matrix& matrix)
     indices.exportIdx(matrix);
 }
 
-template<class T, int n, int m, class F>
-void addToDiagonal(Dune::FieldMatrix<T, n, m>& M, const F& x)
-{
-    for(int i=0; i<n; ++i)
-        M[i][i] += x;
-}
-
 template<class GridView, class Matrix>
 void assemblePQ1Stiffness(const GridView& gridView, Matrix& matrix)
 {
@@ -131,11 +127,11 @@ void assemblePQ1Stiffness(const GridView& gridView, Matrix& matrix)
                     int jGlobal = indexSet.subIndex(element, fe.localCoefficients().localKey(j).subEntity(), dim);
 
                     double zij = (gradients[i] * gradients[j]) * z;
-                    addToDiagonal(matrix[iGlobal][jGlobal], zij);
-                    addToDiagonal(matrix[jGlobal][iGlobal], zij);
+                    Dune::MatrixVector::addToDiagonal(matrix[iGlobal][jGlobal], zij);
+                    Dune::MatrixVector::addToDiagonal(matrix[jGlobal][iGlobal], zij);
                 }
                 double zii = (gradients[i] * gradients[i]) * z;
-                addToDiagonal(matrix[iGlobal][iGlobal], zii);
+                Dune::MatrixVector::addToDiagonal(matrix[iGlobal][iGlobal], zii);
             }
         }
     }
@@ -314,15 +310,19 @@ void markBoundaryDOFs(const GridView& gridView, BitVector& isBoundary)
     }
 }
 
+/**
+ *
+ * \tparam blocksize If this is zero, then double will be used instead of FieldVector/FieldMatrix
+ */
 template <size_t blocksize, class GridView, bool trivialDirichletOnly = true,
           bool makePositiveDefinitive = not trivialDirichletOnly>
 class SymmetricSampleProblem {
 public:
-  typedef typename Dune::FieldMatrix<double, blocksize, blocksize> MatrixBlock;
-  typedef typename Dune::FieldVector<double, blocksize> VectorBlock;
+  using MatrixBlock = std::conditional_t<blocksize==0, double, Dune::FieldMatrix<double, blocksize, blocksize> >;
+  using VectorBlock = std::conditional_t<blocksize==0, double, Dune::FieldVector<double, blocksize> >;
   typedef typename Dune::BCRSMatrix<MatrixBlock> Matrix;
   typedef typename Dune::BlockVector<VectorBlock> Vector;
-  typedef typename Dune::BitSetVector<blocksize> BitVector;
+  using BitVector = Dune::Solvers::DefaultBitVector_t<Vector>;
   typedef EnergyNorm<Matrix, Vector> Norm;
 
   SymmetricSampleProblem(GridView const& gridView) {
@@ -333,15 +333,21 @@ public:
     energyNorm.setMatrix(&A);
 
     ignore.resize(A.N());
+    if constexpr (std::is_same_v<decltype(ignore),Dune::BitSetVector<blocksize> >)
       ignore.unsetAll();
+    else
+      std::fill(ignore.begin(), ignore.end(), false);
     if (trivialDirichletOnly)
         markBoundaryDOFs(gridView, ignore);
     else {
         // Mark the first component only
-        Dune::BitSetVector<blocksize> boundaryNodes(A.N(), false);
+        BitVector boundaryNodes(A.N(), false);
         markBoundaryDOFs(gridView, boundaryNodes);
         for (size_t i=0; i < boundaryNodes.size(); ++i)
+            if constexpr (std::is_same_v<decltype(ignore),Dune::BitSetVector<blocksize> >)
                 ignore[i][0] = boundaryNodes[i][0];
+            else
+                ignore[i] = boundaryNodes[i];
     }
 
     u.resize(A.N());
@@ -357,21 +363,34 @@ public:
     // Set up a random 'solution'
     u = 0;
     for (size_t i = 0; i < u.size(); ++i)
+      if constexpr (blocksize==0)
+      {
+        if (makePositiveDefinitive)
+          A[i][i] += 0.5*std::abs(A[0][0]);
+        u_ex[i] = distribution(generator);
+      } else {
         for (size_t j = 0; j < blocksize; ++j) {
           if (makePositiveDefinitive)
             A[i][i][j][j] += 0.5*std::abs(A[0][0][0][0]);
           u_ex[i][j] = distribution(generator);
         }
+      }
 
     // Construct right hand side corresponding to the 'solution'
     A.mv(u_ex, rhs);
 
     // set dirichlet values
     for (size_t i = 0; i < u.size(); ++i)
+      if constexpr (blocksize==0)
+      {
+        if (ignore[i])
+          u[i] = u_ex[i];
+      } else {
         for (size_t j = 0; j < blocksize; ++j)
           if (ignore[i][j])
             u[i][j] = u_ex[i][j];
       }
+  }
 
   Matrix A;
   Vector u;

dune/solvers/test/gssteptest.cc

@@ -250,7 +250,7 @@ bool checkWithSmallGrids(TestSuite& suite)
     passed = passed and checkWithGrid<Dune::YaspGrid<2, Dune::EquidistantOffsetCoordinates<double,2> > >(suite, 3);
     passed = passed and checkWithGrid<Dune::YaspGrid<3, Dune::EquidistantOffsetCoordinates<double,3> > >(suite, 1);
 
-#if HAVE_UG
+#if HAVE_DUNE_UGGRID
     passed = passed and checkWithGrid<Dune::UGGrid<2> >(suite, 2);
     passed = passed and checkWithGrid<Dune::UGGrid<3> >(suite, 1);
 #endif

dune/solvers/test/umfpacksolvertest.cc

@@ -68,8 +68,10 @@ int main(int argc, char** argv)
 
     bool passed = true;
 
+    UMFPackSolverTestSuite<0> testsuite0;
     UMFPackSolverTestSuite<1> testsuite1;
     UMFPackSolverTestSuite<2> testsuite2;
+    passed = checkWithStandardGrids(testsuite0);
     passed = checkWithStandardGrids(testsuite1);
     passed = checkWithStandardGrids(testsuite2);