Do not create a temporary copy of the matrix even if the ignoreNodes bitset is not empty

This save some memory when the matrix is large.  The dune-istl UMFPack solver
(which we use internally here) actually has support for sending only submatrices
to UMFPack proper.  Unfortunately, this feature is not well-documented...
Also, it only works if all dofs in a block-vector block are either completely
ignored or completely un-ignored.

dune/solvers/solvers/umfpacksolver.hh

@@ -68,7 +68,7 @@ public:
   virtual void solve()
   {
     // We may use the original rhs, but ISTL modifies it, so we need a non-const type here
-    VectorType modifiedRhs = *rhs_;
+    VectorType mutableRhs = *rhs_;
 
     if (not this->ignoreNodes_)
     {
@@ -78,49 +78,67 @@ public:
       Dune::InverseOperatorResult statistics;
       Dune::UMFPack<MatrixType> solver(*matrix_);
       solver.setOption(UMFPACK_PRL, 0);   // no output at all
-      solver.apply(*x_, modifiedRhs, statistics);
+      solver.apply(*x_, mutableRhs, statistics);
 
     }
     else
     {
-      MatrixType modifiedStiffness = *matrix_;
-      //Dune::printmatrix(std::cout, modifiedStiffness, "mod stiffness", "--");
-
-      ///////////////////////////////////////////
-      //   Modify Dirichlet rows
-      ///////////////////////////////////////////
+      ///////////////////////////////////////////////////////////////////////////////////////////
+      //  Extract the set of matrix rows that do not correspond to ignored degrees of freedom.
+      //  Unfortunately, not all cases are handled by the ISTL UMFPack solver.  Currently,
+      //  you can only remove complete block rows.  If a block is only partially ignored,
+      //  the dune-istl UMFPack solver cannot do it, and the code here will throw an exception.
+      //  All this can be fixed, but it needs going into the istl UMFPack code.
+      ///////////////////////////////////////////////////////////////////////////////////////////
+      std::set<std::size_t> nonIgnoreRows;
+      for (size_t i=0; i<matrix_->N(); i++)
+      {
+        auto ignoreCount = (*this->ignoreNodes_)[i].count();
+        if (ignoreCount==0)
+          nonIgnoreRows.insert(i);
+        else if (ignoreCount!=blocksize)
+          DUNE_THROW(Dune::NotImplemented, "Individual blocks must be either ignored completely, or not at all");
+      }
 
-      for (size_t j=0; j<modifiedStiffness.N(); j++)
+      // Construct the solver
+      Dune::InverseOperatorResult statistics;
+      Dune::UMFPack<MatrixType> solver;
+      solver.setOption(UMFPACK_PRL, 0);   // no output at all
+      // We eliminate all rows and columns(!) from the matrix that correspond to ignored degrees of freedom.
+      // Here is where the sparse LU decomposition is happenening.
+      solver.setSubMatrix(*matrix_,nonIgnoreRows);
+
+      // We need to modify the rhs vector by static condensation.
+      VectorType shortRhs(nonIgnoreRows.size());
+      int shortRowCount=0;
+      for (auto it = nonIgnoreRows.begin(); it!=nonIgnoreRows.end(); ++shortRowCount, ++it)
+        shortRhs[shortRowCount] = mutableRhs[*it];
+
+      shortRowCount = 0;
+      for (size_t i=0; i<matrix_->N(); i++)
       {
-        auto cIt    = modifiedStiffness[j].begin();
-        auto cEndIt = modifiedStiffness[j].end();
+        if ((*this->ignoreNodes_)[i][0])
+          continue;
+
+        auto cIt    = (*matrix_)[i].begin();
+        auto cEndIt = (*matrix_)[i].end();
 
         for (; cIt!=cEndIt; ++cIt)
-        {
-          for (int k=0; k<blocksize; k++)
-            if ((*this->ignoreNodes_)[j][k])
-            {
-              (*cIt)[k] = 0;
-              if (cIt.index()==j)
-                (*cIt)[k][k] = 1.0;
-            }
-
-        }
-
-        for (int k=0; k<blocksize; k++)
-          if ((*this->ignoreNodes_)[j][k])
-            modifiedRhs[j][k] = (*x_)[j][k];
+          if ((*this->ignoreNodes_)[cIt.index()][0])
+            cIt->mmv((*x_)[cIt.index()], shortRhs[shortRowCount]);
+
+        shortRowCount++;
       }
 
-      //Dune::printmatrix(std::cout, modifiedStiffness, "mod stiffness", "--");
+      // Solve the reduced system
+      VectorType shortX(nonIgnoreRows.size());
+      solver.apply(shortX, shortRhs, statistics);
+
+      // Blow up the solution vector
+      shortRowCount=0;
+      for (auto it = nonIgnoreRows.begin(); it!=nonIgnoreRows.end(); ++shortRowCount, ++it)
+        (*x_)[*it] = shortX[shortRowCount];
 
-      /////////////////////////////////////////////////////////////////
-      //  Solve the system
-      /////////////////////////////////////////////////////////////////
-      Dune::InverseOperatorResult statistics;
-      Dune::UMFPack<MatrixType> solver(modifiedStiffness);
-      solver.setOption(UMFPACK_PRL, 0);   // no output at all
-      solver.apply(*x_, modifiedRhs, statistics);
     }
   }