Cholmod: Change template specialization to vector type and extend the handling of ignored entries

We interpret Cholmod as an InverseOperator between BlockVectors. Therefore
we dump the BCRSMatrix template.

To clarify the usage of ignore nodes the apply() method now ignores the corresponding entries directly.
The information is stored in a mapping from all indices to the not ignored indices.

Also, the unit tests for cholmod are extended.

dune/istl/cholmod.hh

@@ -6,7 +6,7 @@
 #include <dune/istl/bvector.hh>
 #include<dune/istl/solver.hh>
 
-
+#include <vector>
 #include <memory>
 
 #include <cholmod.h>
@@ -35,25 +35,20 @@ namespace Impl{
 template<class T>
 class Cholmod;
 
-
+/** @brief Dune wrapper for SuiteSparse/CHOLMOD solver
+  *
+  * This class implements an InverseOperator between BlockVector types
+  */
 template<class T, class A, int k>
-class Cholmod<BCRSMatrix<FieldMatrix<T,k,k>, A>>
-  : public InverseOperator<
-      BlockVector<FieldVector<T,k>,
-        typename A::template rebind<FieldVector<T,k>>::other>,
-      BlockVector<FieldVector<T,k>,
-        typename A::template rebind<FieldVector<T,k>>::other>>
+class Cholmod<BlockVector<FieldVector<T,k>, A>>
+  : public InverseOperator<BlockVector<FieldVector<T,k>, A>, BlockVector<FieldVector<T,k>, A>>
 {
 public:
 
   // type of unknown
-  using X = BlockVector<FieldVector<T,k>,
-    typename A::template rebind<FieldVector<T,k>>::other>;
+  using X = BlockVector<FieldVector<T,k>, A>;
   // type of rhs
-  using B = BlockVector<FieldVector<T,k>,
-    typename A::template rebind<FieldVector<T,k>>::other>;
-  // type of external matrix
-  using Matrix = BCRSMatrix<FieldMatrix<T,k,k>, A>;
+  using B = BlockVector<FieldVector<T,k>, A>;
 
 
   /** @brief Default constructor
@@ -78,13 +73,12 @@ public:
     cholmod_finish(&c_);
   }
 
-  // forbit this to avoid freeing memory twice
+  // forbid copying to avoid freeing memory twice
   Cholmod(const Cholmod&) = delete;
   Cholmod& operator=(const Cholmod&) = delete;
 
 
-  /**
-    *  \copydoc InverseOperator::apply(X&,Y&,double,InverseOperatorResult&)
+  /** @brief simple forward to apply(X&, Y&, InverseOperatorResult&)
     */
   void apply (X& x, B& b, double reduction, InverseOperatorResult& res)
   {
@@ -92,20 +86,37 @@ public:
     apply(x,b,res);
   }
 
-  /**
-  *  \copydoc InverseOperator::apply(X&, Y&, InverseOperatorResult&)
+  /** @brief solve the linear system Ax=b (possibly with respect to some ignore field)
+   *
+   * The method assumes that setMatrix() was called before
+   * In the case of a given ignore field the corresponding entries of both in x and b will stay untouched in this method.
   */
   void apply(X& x, B& b, InverseOperatorResult& res)
   {
+    if (x.size() != b.size())
+      DUNE_THROW(Exception, "Error in apply(): sizes of x and b do not match!");
+
+    const auto& blocksize = k;
+
     // cast to double array
     auto b2 = std::make_unique<double[]>(L_->n);
     auto x2 = std::make_unique<double[]>(L_->n);
 
     // copy to cholmod
     auto bp = b2.get();
-    for (const auto& bi : b)
-      for (const auto& bii : bi)
-        *bp++ = bii;
+    if (inverseSubIndices_.empty()) // no ignore field given
+    {
+      // simply copy all values
+      for (const auto& bi : b)
+        for (const auto& bii : bi)
+          *bp++ = bii;
+    }
+    else // use the mapping from not ignored entries
+    {
+      // iterate over inverseSubIndices and resolve the block indices
+      for (const auto& idx : inverseSubIndices_)
+        *bp++ = b[ idx / blocksize ][ idx % blocksize ];
+    }
 
     // create a cholmod dense object
     auto b3 = make_cholmod_dense(cholmod_allocate_dense(L_->n, 1, L_->n, CHOLMOD_REAL, &c_), &c_);
@@ -119,10 +130,19 @@ public:
     auto xp = static_cast<double*>(x3->x);
 
     // copy into x
-    x.resize(L_->n);
-    for (auto& xi : x)
-      for (auto& xii : xi)
-        xii = *xp++;
+    if (inverseSubIndices_.empty()) // no ignore field given
+    {
+      // simply copy all values
+      for (auto i=0; i<x.size(); i++)
+        for (auto ii=0; ii<x[i].size(); ii++)
+          x[i][ii] = *xp++;
+    }
+    else // use the mapping from not ignored entries
+    {
+      // iterate over inverseSubIndices and resolve the block indices
+      for (const auto& idx : inverseSubIndices_)
+        x[ idx / blocksize ][ idx % blocksize ] = *xp++;
+    }
 
     // statistics for a direct solver
     res.iterations = 1;
@@ -135,7 +155,7 @@ public:
    * This method forwards a nullptr to the setMatrix method
    * with ignore nodes
    */
-  void setMatrix(const Matrix& matrix)
+  void setMatrix(const BCRSMatrix<FieldMatrix<T,k,k>>& matrix)
   {
     const Impl::NoIgnore* noIgnore = nullptr;
     setMatrix(matrix, noIgnore);
@@ -144,27 +164,34 @@ public:
   /** @brief Set matrix and ignore nodes
    *
    * The input matrix is copied to CHOLMOD compatible form.
-   * The ignore argument consists of a compatible BitVector,
-   * indicating the dof's which has to be deleted from the matrix
+   * It is possible to ignore some degrees of freedom, provided an ignore field is given with same block structure
+   * like the BlockVector template of the class.
+   *
+   * The ignore field causes the method to ignore both rows and cols of the matrix and therefore operates only
+   * on the reduced quadratic matrix. In case of, e.g., Dirichlet values the user has to take care of proper
+   * adjusting of the rhs before calling apply().
+   *
    * Decomposing the matrix at the end takes a lot of time
    * \param [in] matrix Matrix to be decomposed. In BCRS compatible form
    * \param [in] ignore Pointer to a compatible BitVector
    */
   template<class Ignore>
-  void setMatrix(const Matrix& matrix, const Ignore* ignore)
+  void setMatrix(const BCRSMatrix<FieldMatrix<T,k,k>>& matrix, const Ignore* ignore)
   {
 
+    const auto blocksize = k;
+
     // Total number of rows
-    int N = k * matrix.N();
+    int N = blocksize * matrix.N();
     if ( ignore )
       N -= ignore->count();
 
     // number of nonzeroes
-    const int nnz = k * k * matrix.nonzeroes();
+    const int nnz = blocksize * blocksize * matrix.nonzeroes();
     // number of diagonal entries
-    const int nDiag = k * k * matrix.N();
+    const int nDiag = blocksize * blocksize * matrix.N();
     // number of nonzeroes in the dialgonal
-    const int nnzDiag = (k * (k+1)) / 2 * matrix.N();
+    const int nnzDiag = (blocksize * (blocksize+1)) / 2 * matrix.N();
     // number of nonzeroes in triangular submatrix (including diagonal)
     const int nnzTri = (nnz - nDiag) / 2 + nnzDiag;
 
@@ -192,14 +219,14 @@ public:
     int* Ai = static_cast<int*>(M->i);
     double* Ax = static_cast<double*>(M->x);
 
-    std::size_t blocksize = k;
-
-    // create a vector that maps each remaining matrix index to it's number in the condensed matrix
+    // vector mapping all indices in flat order to the not ignored indices
     std::vector<std::size_t> subIndices;
 
     if ( ignore )
     {
-      subIndices.resize(matrix.M()*blocksize);
+      // init the mappings
+      inverseSubIndices_.resize(N);            // size = number of not ignored entries
+      subIndices.resize(matrix.M()*blocksize); // size = number of all entries
 
       std::size_t j=0;
       for( std::size_t block=0; block<matrix.N(); block++ )
@@ -208,7 +235,9 @@ public:
         {
           if( not (*ignore)[block][i] )
           {
-            subIndices[ block*blocksize + i ] = j++;
+            subIndices[ block*blocksize + i ] = j;
+            inverseSubIndices_[j] = block*blocksize + i;
+            j++;
           }
         }
       }
@@ -279,7 +308,11 @@ private:
 
   cholmod_common c_;
   cholmod_factor* L_ = nullptr;
-};
 
+  // mapping from the not ignored indices in flat order to all indices in flat order
+  // it also holds the info about ignore nodes: if it is empty => no ignore field
+  std::vector<std::size_t> inverseSubIndices_;
+
+};
 
 } /* namespace Dune */

dune/istl/test/cholmodtest.cc

@@ -9,6 +9,9 @@
 #include <dune/istl/bvector.hh>
 #include <dune/istl/io.hh>
 #include <dune/istl/operators.hh>
+#include<array>
+#include<vector>
+
 #include <dune/istl/cholmod.hh>
 
 
@@ -16,17 +19,59 @@
 
 using namespace Dune;
 
-// dummy for empty ignore block
+// One ignore block
+template<int k>
 struct IgnoreBlock
 {
-  bool operator[](std::size_t) const {return false;}
+  std::array<bool,k> data;
+
+  bool operator[](std::size_t i) const
+  {
+    return data[i];
+
+  }
+
+  bool& operator[](std::size_t i)
+  {
+    return data[i];
+
+  }
+
+  std::size_t count() const
+  {
+    size_t res = 0;
+    for(const auto& d : data)
+      if (d)
+        res++;
+
+    return res;
+  }
 };
 
-// dummy for empty ignore block vector
+// Block ignore vector
+template<int blocksize>
 struct Ignore
 {
-  IgnoreBlock operator[](std::size_t) const { return IgnoreBlock(); }
-  std::size_t count() const { return 0; }
+  std::vector<IgnoreBlock<blocksize>> data;
+
+  IgnoreBlock<blocksize> operator[](std::size_t i) const
+  {
+    return data[i];
+  }
+
+  IgnoreBlock<blocksize>& operator[](std::size_t i)
+  {
+    return data[i];
+  }
+
+  std::size_t count() const
+  {
+    size_t res = 0;
+    for(const auto& d : data)
+      res += d.count();
+
+    return res;
+  }
 };
 
 
@@ -36,8 +81,8 @@ int main(int argc, char** argv)
   try
   {
 
-    int N = 300; // number of nodes
-    const int bs = 1; // block size
+    int N = 30; // number of nodes
+    const int bs = 2; // block size
 
     // fill matrix with external method
     BCRSMatrix<FieldMatrix<double,bs,bs>> A;
@@ -45,12 +90,13 @@ int main(int argc, char** argv)
 
     BlockVector<FieldVector<double,bs>> b,x;
     b.resize(A.N());
+    x.resize(A.N());
     b = 1;
 
     InverseOperatorResult res;
 
     // test without ignore nodes
-    Cholmod<BCRSMatrix<FieldMatrix<double,bs,bs>>> cholmod;
+    Cholmod<BlockVector<FieldVector<double,bs>>> cholmod;
     cholmod.setMatrix(A);
     cholmod.apply(x,b,res);
 
@@ -64,16 +110,31 @@ int main(int argc, char** argv)
     b = 1;
 
     // test with ignore nodes
-    Ignore ignore;
-    Cholmod<BCRSMatrix<FieldMatrix<double,bs,bs>>> cholmod2;
+    Ignore<bs> ignore;
+    ignore.data.resize(A.N());
+    // ignore one random entry in x and b
+    ignore[12][0] = true;
+    b[12][0] = 666;
+    x[12][0] = 123;
+
+
+    Cholmod<BlockVector<FieldVector<double,bs>>> cholmod2;
     cholmod2.setMatrix(A,&ignore);
     cholmod2.apply(x,b,res);
 
-    // test
+    // 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;
+
+    // reset the x value
+    x[12][0] = 0;
+    // test -> this should result in zero in every line except entry [12][1]
     A.mmv(x,b);
+    auto b_12_0 = b[12][0];
 
-    if ( b.two_norm() > 1e-9 )
-      std::cerr << " Error in CHOLMOD, residual is too large: " << b.two_norm() << std::endl;
+    // 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;
 
   }
   catch (std::exception &e)