Make UMFPack available as a dune-solvers solver

This class simply forwards all calls to the ISTL UMFPack solver.
As an additional feature, this new implementation honours the ignoreNodes_
field.  Currently, the price for this is an additional copy of the matrix.
I'll fix that in a later release.

dune/solvers/solvers/umfpacksolver.hh

+// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=8 sw=4 sts=4:
+#ifndef DUNE_SOLVERS_SOLVERS_UMFPACKSOLVER_HH
+#define DUNE_SOLVERS_SOLVERS_UMFPACKSOLVER_HH
+
+/** \file
+    \brief A wrapper for the UMFPack sparse direct solver
+*/
+#include <algorithm>
+#include <vector>
+
+#include <dune/common/exceptions.hh>
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/fmatrix.hh>
+
+#include <dune/istl/solver.hh>
+#include <dune/istl/umfpack.hh>
+#include <dune/istl/io.hh>
+
+#include <dune/solvers/common/canignore.hh>
+#include <dune/solvers/solvers/solver.hh>
+
+namespace Dune
+{
+
+namespace Solvers
+{
+
+/** \brief Wraps the UMFPack sparse direct solver */
+template <class MatrixType, class VectorType>
+class UMFPackSolver
+    : public Solver, public CanIgnore<Dune::BitSetVector<VectorType::block_type::dimension> >
+{
+  enum {blocksize = VectorType::block_type::dimension};
+
+  typedef typename VectorType::field_type field_type;
+
+  typedef Dune::FieldMatrix<field_type, blocksize, blocksize> MatrixBlock;
+  typedef Dune::FieldVector<field_type, blocksize>            VectorBlock;
+
+public:
+
+  /** \brief Default constructor */
+  UMFPackSolver ()
+  : Solver(NumProc::FULL)
+  {}
+
+  /** \brief Constructor for a linear problem */
+  UMFPackSolver (const MatrixType& matrix,
+                 VectorType& x,
+                 const VectorType& rhs,
+                 NumProc::VerbosityMode verbosity=NumProc::FULL)
+  : Solver(verbosity),
+    matrix_(&matrix),
+    x_(&x),
+    rhs_(&rhs)
+  {}
+
+  void setProblem(const MatrixType& matrix,
+                  VectorType& x,
+                  const VectorType& rhs)
+  {
+    matrix_ = &matrix;
+    x_ = &x;
+    rhs_ = &rhs;
+  }
+
+  virtual void solve()
+  {
+    // We may use the original rhs, but ISTL modifies it, so we need a non-const type here
+    VectorType modifiedRhs = *rhs_;
+
+    if (not this->ignoreNodes_)
+    {
+      /////////////////////////////////////////////////////////////////
+      //  Solve the system
+      /////////////////////////////////////////////////////////////////
+      Dune::InverseOperatorResult statistics;
+      Dune::UMFPack<MatrixType> solver(*matrix_);
+      solver.setOption(UMFPACK_PRL, 0);   // no output at all
+      solver.apply(*x_, modifiedRhs, statistics);
+
+    }
+    else
+    {
+      MatrixType modifiedStiffness = *matrix_;
+      //Dune::printmatrix(std::cout, modifiedStiffness, "mod stiffness", "--");
+
+      ///////////////////////////////////////////
+      //   Modify Dirichlet rows
+      ///////////////////////////////////////////
+
+      for (size_t j=0; j<modifiedStiffness.N(); j++)
+        for (int k=0; k<blocksize; k++)
+          modifiedStiffness[j][j][k][k] += 0.0;
+
+      for (size_t j=0; j<modifiedStiffness.N(); j++)
+      {
+        auto cIt    = modifiedStiffness[j].begin();
+        auto cEndIt = modifiedStiffness[j].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];
+      }
+
+      //Dune::printmatrix(std::cout, modifiedStiffness, "mod stiffness", "--");
+
+      /////////////////////////////////////////////////////////////////
+      //  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);
+    }
+  }
+
+  ///////////////////////////////////////////////////////
+  // Data
+  ///////////////////////////////////////////////////////
+
+  //! The quadratic term in the quadratic energy
+  const MatrixType* matrix_;
+
+  //! Vector to store the solution
+  VectorType* x_;
+
+  //! The linear term in the quadratic energy
+  const VectorType* rhs_;
+};
+
+}   // namespace Solvers
+
+}   // namespace Dune
+
+#endif