[!271] Implement umfpack for scalar matrices

Merge branch 'implement-umfpack-for-scalar-matrices' into 'master'

See merge request [!271]

  [!271]: Nonecore/dune-istl/merge_requests/271

dune/istl/test/umfpacktest.cc

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
 #include <config.h>
 
-#include <complex>
 #include <iostream>
 
 #include <dune/common/fmatrix.hh>
 #include <dune/common/fvector.hh>
 #include <dune/common/timer.hh>
 #include <dune/istl/bvector.hh>
-#include <dune/istl/colcompmatrix.hh>
-#include <dune/istl/io.hh>
-#include <dune/istl/operators.hh>
 #include <dune/istl/umfpack.hh>
 
 #include "laplacian.hh"
 
 
-int main(int argc, char** argv)
+template<typename Matrix, typename Vector>
+void runUMFPack(std::size_t N)
 {
-#if HAVE_SUITESPARSE_UMFPACK
-  try
-  {
-    typedef double FIELD_TYPE;
-    //typedef std::complex<double> FIELD_TYPE;
+  typedef Dune::MatrixAdapter<Matrix,Vector,Vector> Operator;
 
-    const int BS=1;
-    std::size_t N=100;
+  Matrix mat;
+  Operator fop(mat);
+  Vector b(N*N), x(N*N), b1(N/2), x1(N/2);
 
-    if(argc>1)
-      N = atoi(argv[1]);
-    std::cout<<"testing for N="<<N<<" BS="<<1<<std::endl;
+  setupLaplacian(mat,N);
+  b=1;
+  b1=1;
+  x=0;
+  x1=0;
 
-    typedef Dune::FieldMatrix<FIELD_TYPE,BS,BS> MatrixBlock;
-    typedef Dune::BCRSMatrix<MatrixBlock> BCRSMat;
-    typedef Dune::FieldVector<FIELD_TYPE,BS> VectorBlock;
-    typedef Dune::BlockVector<VectorBlock> Vector;
-    typedef Dune::MatrixAdapter<BCRSMat,Vector,Vector> Operator;
+  Dune::Timer watch;
 
-    BCRSMat mat;
-    Operator fop(mat);
-    Vector b(N*N), x(N*N), b1(N/2), x1(N/2);
+  watch.reset();
 
-    setupLaplacian(mat,N);
-    b=1;
-    b1=1;
-    x=0;
-    x1=0;
+  Dune::UMFPack<Matrix> solver(mat,1);
 
-    Dune::Timer watch;
+  Dune::InverseOperatorResult res;
 
-    watch.reset();
+  solver.apply(x, b, res);
+  solver.free();
 
-    Dune::UMFPack<BCRSMat> solver(mat,1);
+  Dune::UMFPack<Matrix> solver1;
 
-    Dune::InverseOperatorResult res;
+  std::set<std::size_t> mrs;
+  for(std::size_t s=0; s < N/2; ++s)
+    mrs.insert(s);
 
-    solver.apply(x, b, res);
-    solver.free();
+  solver1.setSubMatrix(mat,mrs);
+  solver1.setVerbosity(true);
 
-    Dune::UMFPack<BCRSMat> solver1;
+  solver1.apply(x1,b1, res);
+  solver1.apply(reinterpret_cast<typename Matrix::field_type*>(&x1[0]),
+                reinterpret_cast<typename Matrix::field_type*>(&b1[0]));
 
-    std::set<std::size_t> mrs;
-    for(std::size_t s=0; s < N/2; ++s)
-      mrs.insert(s);
+  Dune::UMFPack<Matrix> save_solver(mat,"umfpack_decomp",0);
+  Dune::UMFPack<Matrix> load_solver(mat,"umfpack_decomp",0);
+}
 
-    solver1.setSubMatrix(mat,mrs);
-    solver1.setVerbosity(true);
+int main(int argc, char** argv) try
+{
+#if HAVE_SUITESPARSE_UMFPACK
+  std::size_t N=100;
 
-    solver1.apply(x1,b1, res);
-    solver1.apply(reinterpret_cast<FIELD_TYPE*>(&x1[0]), reinterpret_cast<FIELD_TYPE*>(&b1[0]));
+  if(argc>1)
+    N = atoi(argv[1]);
 
-    Dune::UMFPack<BCRSMat> save_solver(mat,"umfpack_decomp",0);
-    Dune::UMFPack<BCRSMat> load_solver(mat,"umfpack_decomp",0);
-    return 0;
+  // ------------------------------------------------------------------------------
+  std::cout<<"testing for N="<<N<<" BCRSMatrix<double>"<<std::endl;
+  {
+    using Matrix = Dune::BCRSMatrix<double>;
+    using Vector = Dune::BlockVector<double>;
+    runUMFPack<Matrix,Vector>(N);
   }
-  catch (std::exception &e)
+
+  // ------------------------------------------------------------------------------
+  std::cout<<"testing for N="<<N<<" BCRSMatrix<FieldMatrix<double,1,1> >"<<std::endl;
   {
-    std::cout << "ERROR: " << e.what() << std::endl;
-    return 1;
+    using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1> >;
+    using Vector = Dune::BlockVector<Dune::FieldVector<double,1> >;
+    runUMFPack<Matrix,Vector>(N);
   }
-  catch (...)
+
+  // ------------------------------------------------------------------------------
+  std::cout<<"testing for N="<<N<<" BCRSMatrix<FieldMatrix<double,2,2> >"<<std::endl;
   {
-    std::cerr << "Dune reported an unknown error." << std::endl;
-    exit(1);
+    using Matrix = Dune::BCRSMatrix<Dune::FieldMatrix<double,2,2> >;
+    using Vector = Dune::BlockVector<Dune::FieldVector<double,2> >;
+    runUMFPack<Matrix,Vector>(N);
   }
+
+  return 0;
 #else // HAVE_SUITESPARSE_UMFPACK
   std::cerr << "You need SuiteSparse's UMFPack to run this test." << std::endl;
   return 77;
 #endif // HAVE_SUITESPARSE_UMFPACK
 }
+catch (std::exception &e)
+{
+  std::cout << "ERROR: " << e.what() << std::endl;
+  return 1;
+}

dune/istl/umfpack.hh

@@ -40,15 +40,6 @@ namespace Dune {
   template<class T, bool tag>
   struct SeqOverlappingSchwarzAssemblerHelper;
 
-  /** @brief Use the %UMFPack package to directly solve linear systems -- empty default class
-   * @tparam Matrix the matrix type defining the system
-   * Details on UMFPack can be found on
-   * http://www.cise.ufl.edu/research/sparse/umfpack/
-   */
-  template<class Matrix>
-  class UMFPack
-  {};
-
   // wrapper class for C-Function Calls in the backend. Choose the right function namespace
   // depending on the template parameter used.
   template<typename T>
@@ -172,43 +163,68 @@ namespace Dune {
     }
   };
 
-  /** @brief The %UMFPack direct sparse solver for matrices of type BCRSMatrix
+  namespace Impl
+  {
+    template<class M>
+    struct UMFPackVectorChooser
+    {};
+
+    template<typename T, typename A, int n, int m>
+    struct UMFPackVectorChooser<BCRSMatrix<FieldMatrix<T,n,m>,A > >
+    {
+      /** @brief The type of the domain of the solver */
+      using domain_type = BlockVector<
+                              FieldVector<T,m>,
+                              typename A::template rebind<FieldVector<T,m> >::other>;
+      /** @brief The type of the range of the solver */
+      using range_type  = BlockVector<
+                              FieldVector<T,n>,
+                              typename A::template rebind<FieldVector<T,n> >::other>;
+    };
+
+    template<typename T, typename A>
+    struct UMFPackVectorChooser<BCRSMatrix<T,A> >
+    {
+      /** @brief The type of the domain of the solver */
+      using domain_type = BlockVector<T, A>;
+      /** @brief The type of the range of the solver */
+      using range_type  = BlockVector<T, A>;
+    };
+  }
+
+  /** @brief The %UMFPack direct sparse solver
    *
-   * Specialization for the Dune::BCRSMatrix. %UMFPack will always go double
-   * precision and supports complex numbers
-   * too (use std::complex<double> for that).
+   * Details on UMFPack can be found on
+   * http://www.cise.ufl.edu/research/sparse/umfpack/
+   *
+   * %UMFPack will always use double precision.
+   * For complex matrices use a matrix type with std::complex<double>
+   * as the underlying number type.
    *
-   * \tparam T Number type.  Only double and std::complex<double> is supported
-   * \tparam A STL-compatible allocator type
-   * \tparam n Number of rows in a matrix block
-   * \tparam m Number of columns in a matrix block
+   * \tparam Matrix the matrix type defining the system
    *
    * \note This will only work if dune-istl has been configured to use UMFPack
    */
-  template<typename T, typename A, int n, int m>
-  class UMFPack<BCRSMatrix<FieldMatrix<T,n,m>,A > >
+  template<typename M>
+  class UMFPack
       : public InverseOperator<
-          BlockVector<FieldVector<T,m>,
-              typename A::template rebind<FieldVector<T,m> >::other>,
-          BlockVector<FieldVector<T,n>,
-              typename A::template rebind<FieldVector<T,n> >::other> >
+          typename Impl::UMFPackVectorChooser<M>::domain_type,
+          typename Impl::UMFPackVectorChooser<M>::range_type >
   {
+    using T = typename M::field_type;
+
     public:
     /** @brief The matrix type. */
-    typedef Dune::BCRSMatrix<FieldMatrix<T,n,m>,A> Matrix;
-    typedef Dune::BCRSMatrix<FieldMatrix<T,n,m>,A> matrix_type;
-    /** @brief The corresponding SuperLU Matrix type.*/
+    using Matrix = M;
+    using matrix_type = M;
+    /** @brief The corresponding UMFPack matrix type.*/
     typedef Dune::ColCompMatrix<Matrix> UMFPackMatrix;
     /** @brief Type of an associated initializer class. */
-    typedef ColCompMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> > MatrixInitializer;
+    typedef ColCompMatrixInitializer<M> MatrixInitializer;
     /** @brief The type of the domain of the solver. */
-    typedef Dune::BlockVector<
-        FieldVector<T,m>,
-        typename A::template rebind<FieldVector<T,m> >::other> domain_type;
+    using domain_type = typename Impl::UMFPackVectorChooser<M>::domain_type;
     /** @brief The type of the range of the solver. */
-    typedef Dune::BlockVector<
-        FieldVector<T,n>,
-        typename A::template rebind<FieldVector<T,n> >::other> range_type;
+    using range_type = typename Impl::UMFPackVectorChooser<M>::range_type;
 
     //! Category of the solver (see SolverCategory::Category)
     virtual SolverCategory::Category category() const
@@ -216,7 +232,7 @@ namespace Dune {
       return SolverCategory::Category::sequential;
     }
 
-    /** @brief Construct a solver object from a BCRSMatrix
+    /** @brief Construct a solver object from a matrix
      *
      * This computes the matrix decomposition, and may take a long time
      * (and use a lot of memory).
@@ -483,7 +499,7 @@ namespace Dune {
     private:
     typedef typename Dune::UMFPackMethodChooser<T> Caller;
 
-    template<class M,class X, class TM, class TD, class T1>
+    template<class Mat,class X, class TM, class TD, class T1>
     friend class SeqOverlappingSchwarz;
     friend struct SeqOverlappingSchwarzAssemblerHelper<UMFPack<Matrix>,true>;