cmake/modules/FindSuperLU.cmake

@@ -125,8 +125,20 @@ int main(void)
 }"
 SUPERLU_MIN_VERSION_5)
 
+if(SUPERLU_MIN_VERSION_4_3)
+  include(CheckIncludeFiles)
+  set(HAVE_SLU_DDEFS_H 1)
+  check_include_files(slu_sdefs.h HAVE_SLU_SDEFS_H)
+  check_include_files(slu_cdefs.h HAVE_SLU_CDEFS_H)
+  check_include_files(slu_zdefs.h HAVE_SLU_ZDEFS_H)
+endif(SUPERLU_MIN_VERSION_4_3)
+
 cmake_pop_check_state()
 
+set(SUPERLU_INT_TYPE "int" CACHE STRING
+  "The integer version that SuperLU was compiled for (Default is int.
+  Should be the same as int_t define in e.g. slu_sdefs.h")
+
 if(NOT SUPERLU_MIN_VERSION_4)
   set(SUPERLU_WITH_VERSION "SuperLU < 4.0" CACHE STRING
     "Human readable string containing SuperLU version information.")

config.h.cmake

@@ -31,6 +31,22 @@
 /* Define to ENABLE_SUPERLU if the SuperLU library is available */
 #cmakedefine HAVE_SUPERLU ENABLE_SUPERLU
 
+/* Define to the integer type that SuperLU was compiled for
+   See e.g. what int_t is defined to in slu_sdefs.h */
+#cmakedefine SUPERLU_INT_TYPE @SUPERLU_INT_TYPE@
+
+/* Define to 1 if header slu_sdefs.h is there. */
+#cmakedefine01 HAVE_SLU_SDEFS_H
+
+/* Define to 1 if header slu_ddefs.h is there. */
+#cmakedefine01 HAVE_SLU_DDEFS_H
+
+/* Define to 1 if header slu_cdefs.h is there. */
+#cmakedefine01 HAVE_SLU_CDEFS_H
+
+/* Define to 1 if header slu_zdefs.h is there. */
+#cmakedefine01 HAVE_SLU_ZDEFS_H
+
 /* Define to ENABLE_ARPACKPP if the ARPACK++ library is available */
 #cmakedefine HAVE_ARPACKPP ENABLE_ARPACKPP
 

doc/doxygen/modules.txt

 /**
   @defgroup ISTL Iterative Solvers Template Library (ISTL)
-    
+
    @brief Iterative Solvers supporting block recursive matrix and
    vector classes at compile time.
 
-   The Iterative Solver Template Library applies generic programming 
-   in C++ to the 
+   The Iterative Solver Template Library applies generic programming
+   in C++ to the
    domain of iterative solvers of linear systems stemming
    from finite element discretizations. Those discretizations exhibit a
    lot of structure, e.g:
@@ -13,7 +13,7 @@
    <ol>
    <li>Certain discretizations for systems of PDEs or higher order
   methods result in matrices where individual entries are replaced by
-  small blocks, say of size \f$2\times 2\f$ or \f$4\times 4\f$. 
+  small blocks, say of size \f$2\times 2\f$ or \f$4\times 4\f$.
   Dense blocks of different sizes
   e.g. arise in \f$hp\f$ Discontinuous Galerkin discretization methods.
   It is straightforward and
@@ -30,13 +30,22 @@
 
    Our matrix
    and vector interface supports a block recursive structure. Each
-   sparse matrix entry can itself be either a sparse or a small 
+   sparse matrix entry can itself be either a sparse or a small
    dense matrix.
-   
-   The solvers use this recursive block structure via template meta 
+
+   The solvers use this recursive block structure via template meta
    programming at compile time.
 
-   ISTL consists of the \ref ISTL_SPMV "matrix and vector API" and 
+   ISTL consists of the \ref ISTL_SPMV "matrix and vector API" and
    the \ref ISTL_Solvers "solvers" which use the \ref ISTL_Prec preconditioners.
  */
 
+/** @defgroup ISTL_Solvers Iterative Solvers
+    @ingroup ISTL
+*/
+/** @defgroup ISTL_Eigenvalue Eigenvalue Solvers
+    @ingroup ISTL_Solvers
+*/
+/** @defgroup ISTL_Comm Communication Interface
+    @ingroup ISTL
+*/

doc/istl.tex

@@ -760,7 +760,7 @@ might be needed on new platforms.
 
 
 \section{Performance}
-We evaluated the performance of our implementation on a Petium 4 Mobile
+We evaluated the performance of our implementation on a Pentium 4 Mobile
 2.4 GHz with a measured memory bandwidth of 1084 MB/s for the daypy
 operation ($x = y + \alpha z$) in Tables
 \ref{tab:istl_performance}. 

dune/istl/CMakeLists.txt

@@ -42,6 +42,7 @@ install(FILES
    solvertype.hh
    spqr.hh
    superlu.hh
+   superlufunctions.hh
    supermatrix.hh
    umfpack.hh
    vbvector.hh

dune/istl/basearray.hh

@@ -18,6 +18,9 @@
 
 namespace Dune {
 
+/** \brief Everything in this namespace is internal to dune-istl, and may change without warning */
+namespace Imp {
+
   /**  \brief A simple array container for objects of type B
 
      Implement.
@@ -748,6 +751,8 @@ namespace Dune {
     size_type* j;     // the index set
   };
 
+} // end namespace Imp
+
 } // end namespace
 
 #endif

dune/istl/bcrsmatrix.hh

@@ -451,7 +451,7 @@ namespace Dune {
     typedef A allocator_type;
 
     //! implement row_type with compressed vector
-    typedef CompressedBlockVectorWindow<B,A> row_type;
+    typedef Imp::CompressedBlockVectorWindow<B,A> row_type;
 
     //! The type for the index access and the size
     typedef typename A::size_type size_type;

dune/istl/bvector.hh

@@ -28,6 +28,9 @@ namespace Dune {
   template<class B, class A=std::allocator<B> >
   class BlockVectorWindow;
 
+/** \brief Everything in this namespace is internal to dune-istl, and may change without warning */
+namespace Imp {
+
   /**
       \brief An unmanaged vector of blocks.
 
@@ -294,6 +297,7 @@ namespace Dune {
     {       }
   };
 
+} // end namespace Imp
   /**
      @addtogroup ISTL_SPMV
      @{
@@ -309,7 +313,7 @@ namespace Dune {
           enables error checking.
    */
   template<class B, class A=std::allocator<B> >
-  class BlockVector : public block_vector_unmanaged<B,A>
+  class BlockVector : public Imp::block_vector_unmanaged<B,A>
   {
   public:
 
@@ -334,15 +338,15 @@ namespace Dune {
     };
 
     //! make iterators available as types
-    typedef typename block_vector_unmanaged<B,A>::Iterator Iterator;
+    typedef typename Imp::block_vector_unmanaged<B,A>::Iterator Iterator;
 
     //! make iterators available as types
-    typedef typename block_vector_unmanaged<B,A>::ConstIterator ConstIterator;
+    typedef typename Imp::block_vector_unmanaged<B,A>::ConstIterator ConstIterator;
 
     //===== constructors and such
 
     //! makes empty vector
-    BlockVector () : block_vector_unmanaged<B,A>(),
+    BlockVector () : Imp::block_vector_unmanaged<B,A>(),
                      capacity_(0)
     {}
 
@@ -436,7 +440,7 @@ namespace Dune {
      */
     void reserve(size_type capacity, bool copyOldValues=true)
     {
-      if(capacity >= block_vector_unmanaged<B,A>::N() && capacity != capacity_) {
+      if(capacity >= Imp::block_vector_unmanaged<B,A>::N() && capacity != capacity_) {
         // save the old data
         B* pold = this->p;
 
@@ -450,7 +454,7 @@ namespace Dune {
             B* to = this->p;
             B* from = pold;
 
-            for(size_type i=0; i < block_vector_unmanaged<B,A>::N(); ++i, ++from, ++to)
+            for(size_type i=0; i < Imp::block_vector_unmanaged<B,A>::N(); ++i, ++from, ++to)
               *to = *from;
           }
           if(capacity_ > 0) {
@@ -498,7 +502,7 @@ namespace Dune {
      */
     void resize(size_type size, bool copyOldValues=true)
     {
-      if(size > block_vector_unmanaged<B,A>::N())
+      if (size > Imp::block_vector_unmanaged<B,A>::N())
         if(capacity_ < size)
           this->reserve(size, copyOldValues);
       this->n = size;
@@ -509,7 +513,7 @@ namespace Dune {
 
     //! copy constructor
     BlockVector (const BlockVector& a) :
-      block_vector_unmanaged<B,A>(a)
+      Imp::block_vector_unmanaged<B,A>(a)
     {
       // allocate memory with same size as a
       this->n = a.n;
@@ -575,7 +579,7 @@ namespace Dune {
     BlockVector& operator= (const field_type& k)
     {
       // forward to operator= in base class
-      (static_cast<block_vector_unmanaged<B,A>&>(*this)) = k;
+      (static_cast<Imp::block_vector_unmanaged<B,A>&>(*this)) = k;
       return *this;
     }
 
@@ -623,6 +627,9 @@ namespace Dune {
     return s;
   }
 
+/** \brief Everything in this namespace is internal to dune-istl, and may change without warning */
+namespace Imp {
+
   /** BlockVectorWindow adds window manipulation functions
           to the block_vector_unmanaged template.
 
@@ -646,7 +653,7 @@ namespace Dune {
 #else
   template<class B, class A=std::allocator<B> >
 #endif
-  class BlockVectorWindow : public block_vector_unmanaged<B,A>
+  class BlockVectorWindow : public Imp::block_vector_unmanaged<B,A>
   {
   public:
 
@@ -671,15 +678,15 @@ namespace Dune {
     };
 
     //! make iterators available as types
-    typedef typename block_vector_unmanaged<B,A>::Iterator Iterator;
+    typedef typename Imp::block_vector_unmanaged<B,A>::Iterator Iterator;
 
     //! make iterators available as types
-    typedef typename block_vector_unmanaged<B,A>::ConstIterator ConstIterator;
+    typedef typename Imp::block_vector_unmanaged<B,A>::ConstIterator ConstIterator;
 
 
     //===== constructors and such
     //! makes empty array
-    BlockVectorWindow () : block_vector_unmanaged<B,A>()
+    BlockVectorWindow () : Imp::block_vector_unmanaged<B,A>()
     {       }
 
     //! make array from given pointer and size
@@ -715,7 +722,7 @@ namespace Dune {
     //! assign from scalar
     BlockVectorWindow& operator= (const field_type& k)
     {
-      (static_cast<block_vector_unmanaged<B,A>&>(*this)) = k;
+      (static_cast<Imp::block_vector_unmanaged<B,A>&>(*this)) = k;
       return *this;
     }
 
@@ -1156,6 +1163,8 @@ namespace Dune {
     }
   };
 
-} // end namespace
+} // end namespace 'Imp'
+
+} // end namespace 'Dune'
 
 #endif

dune/istl/eigenvalue/arpackpp.hh

@@ -3,7 +3,7 @@
 #ifndef DUNE_ISTL_EIGENVALUE_ARPACKPP_HH
 #define DUNE_ISTL_EIGENVALUE_ARPACKPP_HH
 
-#if HAVE_ARPACKPP
+#if HAVE_ARPACKPP || defined DOXYGEN
 
 #include <cmath>  // provides std::abs, std::pow, std::sqrt
 
@@ -26,176 +26,182 @@
 #endif
 #include "arssym.h"  // provides ARSymStdEig
 
-
 namespace Dune
 {
 
-  /**
-   * \brief Wrapper for a DUNE-ISTL BCRSMatrix which can be used
-   *        together with those algorithms of the ARPACK++ library
-   *        which solely perform the products A*v and/or A^T*A*v
-   *        and/or A*A^T*v.
-   *
-   * \todo The current implementation is limited to DUNE-ISTL
-   *       BCRSMatrix types with blocklevel 2. An extension to
-   *       blocklevel >= 2 might be provided in a future version.
-   *
-   * \tparam BCRSMatrix Type of a DUNE-ISTL BCRSMatrix;
-   *                    is assumed to have blocklevel 2.
-   *
-   * \author Sebastian Westerheide.
-   */
-  template <class BCRSMatrix>
-  class ArPackPlusPlus_BCRSMatrixWrapper
-  {
-  public:
-    //! Type of the underlying field of the matrix
-    typedef typename BCRSMatrix::field_type Real;
-
-  public:
-    //! Construct from BCRSMatrix A
-    ArPackPlusPlus_BCRSMatrixWrapper (const BCRSMatrix& A)
-      : A_(A),
-        m_(A_.M() * mBlock), n_(A_.N() * nBlock)
-    {
-      // assert that BCRSMatrix type has blocklevel 2
-      static_assert
-        (BCRSMatrix::blocklevel == 2,
-         "Only BCRSMatrices with blocklevel 2 are supported.");
-
-      // allocate memory for auxiliary block vector objects
-      // which are compatible to matrix rows / columns
-      domainBlockVector.resize(A_.N(),false);
-      rangeBlockVector.resize(A_.M(),false);
-    }
-
-    //! Perform matrix-vector product w = A*v
-    inline void multMv (Real* v, Real* w)
-    {
-      // get vector v as an object of appropriate type
-      arrayToDomainBlockVector(v,domainBlockVector);
+  /** @addtogroup ISTL_Eigenvalue
+    @{
+  */
 
-      // perform matrix-vector product
-      A_.mv(domainBlockVector,rangeBlockVector);
-
-      // get vector w from object of appropriate type
-      rangeBlockVectorToArray(rangeBlockVector,w);
-    };
-
-    //! Perform matrix-vector product w = A^T*A*v
-    inline void multMtMv (Real* v, Real* w)
+  namespace Impl {
+    /**
+     *        Wrapper for a DUNE-ISTL BCRSMatrix which can be used
+     *        together with those algorithms of the ARPACK++ library
+     *        which solely perform the products A*v and/or A^T*A*v
+     *        and/or A*A^T*v.
+     *
+     * \todo The current implementation is limited to DUNE-ISTL
+     *       BCRSMatrix types with blocklevel 2. An extension to
+     *       blocklevel >= 2 might be provided in a future version.
+     *
+     * \tparam BCRSMatrix Type of a DUNE-ISTL BCRSMatrix;
+     *                    is assumed to have blocklevel 2.
+     *
+     * \author Sebastian Westerheide.
+     */
+    template <class BCRSMatrix>
+    class ArPackPlusPlus_BCRSMatrixWrapper
     {
-      // get vector v as an object of appropriate type
-      arrayToDomainBlockVector(v,domainBlockVector);
-
-      // perform matrix-vector product
-      A_.mv(domainBlockVector,rangeBlockVector);
-      A_.mtv(rangeBlockVector,domainBlockVector);
+    public:
+      //! Type of the underlying field of the matrix
+      typedef typename BCRSMatrix::field_type Real;
+
+    public:
+      //! Construct from BCRSMatrix A
+      ArPackPlusPlus_BCRSMatrixWrapper (const BCRSMatrix& A)
+        : A_(A),
+          m_(A_.M() * mBlock), n_(A_.N() * nBlock)
+      {
+        // assert that BCRSMatrix type has blocklevel 2
+        static_assert
+          (BCRSMatrix::blocklevel == 2,
+            "Only BCRSMatrices with blocklevel 2 are supported.");
+
+        // allocate memory for auxiliary block vector objects
+        // which are compatible to matrix rows / columns
+        domainBlockVector.resize(A_.N(),false);
+        rangeBlockVector.resize(A_.M(),false);
+      }
 
-      // get vector w from object of appropriate type
-      domainBlockVectorToArray(domainBlockVector,w);
-    };
+      //! Perform matrix-vector product w = A*v
+      inline void multMv (Real* v, Real* w)
+      {
+        // get vector v as an object of appropriate type
+        arrayToDomainBlockVector(v,domainBlockVector);
 
-    //! Perform matrix-vector product w = A*A^T*v
-    inline void multMMtv (Real* v, Real* w)
-    {
-      // get vector v as an object of appropriate type
-      arrayToRangeBlockVector(v,rangeBlockVector);
+        // perform matrix-vector product
+        A_.mv(domainBlockVector,rangeBlockVector);
 
-      // perform matrix-vector product
-      A_.mtv(rangeBlockVector,domainBlockVector);
-      A_.mv(domainBlockVector,rangeBlockVector);
+        // get vector w from object of appropriate type
+        rangeBlockVectorToArray(rangeBlockVector,w);
+      };
 
-      // get vector w from object of appropriate type
-      rangeBlockVectorToArray(rangeBlockVector,w);
-    };
+      //! Perform matrix-vector product w = A^T*A*v
+      inline void multMtMv (Real* v, Real* w)
+      {
+        // get vector v as an object of appropriate type
+        arrayToDomainBlockVector(v,domainBlockVector);
 
-    //! Return number of rows in the matrix
-    inline int nrows () const { return m_; }
+        // perform matrix-vector product
+        A_.mv(domainBlockVector,rangeBlockVector);
+        A_.mtv(rangeBlockVector,domainBlockVector);
 
-    //! Return number of columns in the matrix
-    inline int ncols () const { return n_; }
+        // get vector w from object of appropriate type
+        domainBlockVectorToArray(domainBlockVector,w);
+      };
 
-  protected:
-    // Number of rows and columns in each block of the matrix
-    constexpr static int mBlock = BCRSMatrix::block_type::rows;
-    constexpr static int nBlock = BCRSMatrix::block_type::cols;
-
-    // Type of vectors in the domain of the linear map associated with
-    // the matrix, i.e. block vectors compatible to matrix rows
-    constexpr static int dbvBlockSize = nBlock;
-    typedef Dune::FieldVector<Real,dbvBlockSize> DomainBlockVectorBlock;
-    typedef Dune::BlockVector<DomainBlockVectorBlock> DomainBlockVector;
-
-    // Type of vectors in the range of the linear map associated with
-    // the matrix, i.e. block vectors compatible to matrix columns
-    constexpr static int rbvBlockSize = mBlock;
-    typedef Dune::FieldVector<Real,rbvBlockSize> RangeBlockVectorBlock;
-    typedef Dune::BlockVector<RangeBlockVectorBlock> RangeBlockVector;
-
-    // Types for vector index access
-    typedef typename DomainBlockVector::size_type dbv_size_type;
-    typedef typename RangeBlockVector::size_type rbv_size_type;
-    typedef typename DomainBlockVectorBlock::size_type dbvb_size_type;
-    typedef typename RangeBlockVectorBlock::size_type rbvb_size_type;
-
-    // Get vector v from a block vector object which is compatible to
-    // matrix rows
-    static inline void
+      //! Perform matrix-vector product w = A*A^T*v
+      inline void multMMtv (Real* v, Real* w)
+      {
+        // get vector v as an object of appropriate type
+        arrayToRangeBlockVector(v,rangeBlockVector);
+
+        // perform matrix-vector product
+        A_.mtv(rangeBlockVector,domainBlockVector);
+        A_.mv(domainBlockVector,rangeBlockVector);
+
+        // get vector w from object of appropriate type
+        rangeBlockVectorToArray(rangeBlockVector,w);
+      };
+
+      //! Return number of rows in the matrix
+      inline int nrows () const { return m_; }
+
+      //! Return number of columns in the matrix
+      inline int ncols () const { return n_; }
+
+    protected:
+      // Number of rows and columns in each block of the matrix
+      constexpr static int mBlock = BCRSMatrix::block_type::rows;
+      constexpr static int nBlock = BCRSMatrix::block_type::cols;
+
+      // Type of vectors in the domain of the linear map associated with
+      // the matrix, i.e. block vectors compatible to matrix rows
+      constexpr static int dbvBlockSize = nBlock;
+      typedef Dune::FieldVector<Real,dbvBlockSize> DomainBlockVectorBlock;
+      typedef Dune::BlockVector<DomainBlockVectorBlock> DomainBlockVector;
+
+      // Type of vectors in the range of the linear map associated with
+      // the matrix, i.e. block vectors compatible to matrix columns
+      constexpr static int rbvBlockSize = mBlock;
+      typedef Dune::FieldVector<Real,rbvBlockSize> RangeBlockVectorBlock;
+      typedef Dune::BlockVector<RangeBlockVectorBlock> RangeBlockVector;
+
+      // Types for vector index access
+      typedef typename DomainBlockVector::size_type dbv_size_type;
+      typedef typename RangeBlockVector::size_type rbv_size_type;
+      typedef typename DomainBlockVectorBlock::size_type dbvb_size_type;
+      typedef typename RangeBlockVectorBlock::size_type rbvb_size_type;
+
+      // Get vector v from a block vector object which is compatible to
+      // matrix rows
+      static inline void
       domainBlockVectorToArray (const DomainBlockVector& dbv, Real* v)
-    {
-      for (dbv_size_type block = 0; block < dbv.N(); ++block)
-        for (dbvb_size_type iBlock = 0; iBlock < dbvBlockSize; ++iBlock)
-          v[block*dbvBlockSize + iBlock] = dbv[block][iBlock];
-    }
+      {
+        for (dbv_size_type block = 0; block < dbv.N(); ++block)
+          for (dbvb_size_type iBlock = 0; iBlock < dbvBlockSize; ++iBlock)
+            v[block*dbvBlockSize + iBlock] = dbv[block][iBlock];
+      }
 
-    // Get vector v from a block vector object which is compatible to
-    // matrix columns
-    static inline void
+      // Get vector v from a block vector object which is compatible to
+      // matrix columns
+      static inline void
       rangeBlockVectorToArray (const RangeBlockVector& rbv, Real* v)
-    {
-      for (rbv_size_type block = 0; block < rbv.N(); ++block)
-        for (rbvb_size_type iBlock = 0; iBlock < rbvBlockSize; ++iBlock)
-          v[block*rbvBlockSize + iBlock] = rbv[block][iBlock];
-    }
-
-  public:
-    //! Get vector v as a block vector object which is compatible to
-    //! matrix rows
-    static inline void arrayToDomainBlockVector (const Real* v,
-                                                 DomainBlockVector& dbv)
-    {
-      for (dbv_size_type block = 0; block < dbv.N(); ++block)
-        for (dbvb_size_type iBlock = 0; iBlock < dbvBlockSize; ++iBlock)
-          dbv[block][iBlock] = v[block*dbvBlockSize + iBlock];
-    }
+      {
+        for (rbv_size_type block = 0; block < rbv.N(); ++block)
+          for (rbvb_size_type iBlock = 0; iBlock < rbvBlockSize; ++iBlock)
+            v[block*rbvBlockSize + iBlock] = rbv[block][iBlock];
+      }
 
-    //! Get vector v as a block vector object which is compatible to
-    //! matrix columns
-    static inline void arrayToRangeBlockVector (const Real* v,
-                                                RangeBlockVector& rbv)
-    {
-      for (rbv_size_type block = 0; block < rbv.N(); ++block)
-        for (rbvb_size_type iBlock = 0; iBlock < rbvBlockSize; ++iBlock)
-          rbv[block][iBlock] = v[block*rbvBlockSize + iBlock];
-    }
+    public:
+      //! Get vector v as a block vector object which is compatible to
+      //! matrix rows
+      static inline void arrayToDomainBlockVector (const Real* v,
+        DomainBlockVector& dbv)
+      {
+        for (dbv_size_type block = 0; block < dbv.N(); ++block)
+          for (dbvb_size_type iBlock = 0; iBlock < dbvBlockSize; ++iBlock)
+            dbv[block][iBlock] = v[block*dbvBlockSize + iBlock];
+      }
 
-  protected:
-    // The DUNE-ISTL BCRSMatrix
-    const BCRSMatrix& A_;
+      //! Get vector v as a block vector object which is compatible to
+      //! matrix columns
+      static inline void arrayToRangeBlockVector (const Real* v,
+        RangeBlockVector& rbv)
+      {
+        for (rbv_size_type block = 0; block < rbv.N(); ++block)
+          for (rbvb_size_type iBlock = 0; iBlock < rbvBlockSize; ++iBlock)
+            rbv[block][iBlock] = v[block*rbvBlockSize + iBlock];
+      }
 
-    // Number of rows and columns in the matrix
-    const int m_, n_;
+    protected:
+      // The DUNE-ISTL BCRSMatrix
+      const BCRSMatrix& A_;
 
-    // Auxiliary block vector objects which are
-    // compatible to matrix rows / columns
-    mutable DomainBlockVector domainBlockVector;
-    mutable RangeBlockVector rangeBlockVector;
-  };
+      // Number of rows and columns in the matrix
+      const int m_, n_;
 
+      // Auxiliary block vector objects which are
+      // compatible to matrix rows / columns
+      mutable DomainBlockVector domainBlockVector;
+      mutable RangeBlockVector rangeBlockVector;
+    };
+  } // end namespace Impl
 
   /**
-   * \brief A class template for performing some eigenvalue algorithms
+   * \brief Wrapper to use a range of ARPACK++ eigenvalue solvers
+   *
+   *        A class template for performing some eigenvalue algorithms
    *        provided by the ARPACK++ library which is based on the implicitly
    *        restarted Arnoldi/Lanczos method (IRAM/IRLM), a synthesis of the
    *        Arnoldi/Lanczos process with the implicitily shifted QR technique.
@@ -288,7 +294,7 @@ namespace Dune
 
       // use type ArPackPlusPlus_BCRSMatrixWrapper to store matrix information
       // and to perform the product A*v (LU decomposition is not used)
-      typedef ArPackPlusPlus_BCRSMatrixWrapper<BCRSMatrix> WrappedMatrix;
+      typedef Impl::ArPackPlusPlus_BCRSMatrixWrapper<BCRSMatrix> WrappedMatrix;
       WrappedMatrix A(m_);
 
       // get number of rows and columns in A
@@ -390,7 +396,7 @@ namespace Dune
 
       // use type ArPackPlusPlus_BCRSMatrixWrapper to store matrix information
       // and to perform the product A*v (LU decomposition is not used)
-      typedef ArPackPlusPlus_BCRSMatrixWrapper<BCRSMatrix> WrappedMatrix;
+      typedef Impl::ArPackPlusPlus_BCRSMatrixWrapper<BCRSMatrix> WrappedMatrix;
       WrappedMatrix A(m_);
 
       // get number of rows and columns in A
@@ -491,7 +497,7 @@ namespace Dune
 
       // use type ArPackPlusPlus_BCRSMatrixWrapper to store matrix information
       // and to perform the product A*v (LU decomposition is not used)
-      typedef ArPackPlusPlus_BCRSMatrixWrapper<BCRSMatrix> WrappedMatrix;
+      typedef Impl::ArPackPlusPlus_BCRSMatrixWrapper<BCRSMatrix> WrappedMatrix;
       WrappedMatrix A(m_);
 
       // get number of rows and columns in A
@@ -608,7 +614,7 @@ namespace Dune
 
       // use type ArPackPlusPlus_BCRSMatrixWrapper to store matrix information
       // and to perform the product A^T*A*v (LU decomposition is not used)
-      typedef ArPackPlusPlus_BCRSMatrixWrapper<BCRSMatrix> WrappedMatrix;
+      typedef Impl::ArPackPlusPlus_BCRSMatrixWrapper<BCRSMatrix> WrappedMatrix;
       WrappedMatrix A(m_);
 
       // get number of rows and columns in A
@@ -720,7 +726,7 @@ namespace Dune
 
       // use type ArPackPlusPlus_BCRSMatrixWrapper to store matrix information
       // and to perform the product A^T*A*v (LU decomposition is not used)
-      typedef ArPackPlusPlus_BCRSMatrixWrapper<BCRSMatrix> WrappedMatrix;
+      typedef Impl::ArPackPlusPlus_BCRSMatrixWrapper<BCRSMatrix> WrappedMatrix;
       WrappedMatrix A(m_);
 
       // get number of rows and columns in A
@@ -828,7 +834,7 @@ namespace Dune
 
       // use type ArPackPlusPlus_BCRSMatrixWrapper to store matrix information
       // and to perform the product A^T*A*v (LU decomposition is not used)
-      typedef ArPackPlusPlus_BCRSMatrixWrapper<BCRSMatrix> WrappedMatrix;
+      typedef Impl::ArPackPlusPlus_BCRSMatrixWrapper<BCRSMatrix> WrappedMatrix;
       WrappedMatrix A(m_);
 
       // get number of rows and columns in A
@@ -958,8 +964,9 @@ namespace Dune
     const std::string blank_;
   };
 
-}  // namespace Dune
+  /** @} */
 
+}  // namespace Dune
 
 #endif  // HAVE_ARPACKPP
 

dune/istl/eigenvalue/poweriteration.hh

@@ -24,154 +24,108 @@
 #include <dune/istl/io.hh>              // provides Dune::printvector(...)
 #include <dune/istl/solvers.hh>         // provides Dune::InverseOperatorResult
 
-
 namespace Dune
 {
 
-  /**
-   * \brief A linear operator scaling vectors by a scalar value.
-   *        The scalar value can be changed as it is given in a
-   *        form decomposed into an immutable and a mutable part.
-   *
-   * \author Sebastian Westerheide.
-   */
-  template <class X, class Y = X>
-  class ScalingLinearOperator : public Dune::LinearOperator<X,Y>
-  {
-  public:
-    typedef X domain_type;
-    typedef Y range_type;
-    typedef typename X::field_type field_type;
-
-    enum {category = Dune::SolverCategory::sequential};
-
-    ScalingLinearOperator (field_type immutable_scaling,
-                           const field_type& mutable_scaling)
-      : immutable_scaling_(immutable_scaling),
-        mutable_scaling_(mutable_scaling)
-    {}
+  /** @addtogroup ISTL_Eigenvalue
+    @{
+  */
 
-    virtual void apply (const X& x, Y& y) const
-    {
-      y = x;
-      y *= immutable_scaling_*mutable_scaling_;
-    }
-
-    virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const
+  namespace Impl {
+    /**
+     * \brief A linear operator scaling vectors by a scalar value.
+     *        The scalar value can be changed as it is given in a
+     *        form decomposed into an immutable and a mutable part.
+     *
+     * \author Sebastian Westerheide.
+     */
+    template <class X, class Y = X>
+    class ScalingLinearOperator : public Dune::LinearOperator<X,Y>
     {
-      X temp(x);
-      temp *= immutable_scaling_*mutable_scaling_;
-      y.axpy(alpha,temp);
-    }
+    public:
+      typedef X domain_type;
+      typedef Y range_type;
+      typedef typename X::field_type field_type;
 
-  protected:
-    const field_type immutable_scaling_;
-    const field_type& mutable_scaling_;
-  };
+      enum {category = Dune::SolverCategory::sequential};
 
+      ScalingLinearOperator (field_type immutable_scaling,
+        const field_type& mutable_scaling)
+        : immutable_scaling_(immutable_scaling),
+          mutable_scaling_(mutable_scaling)
+      {}
 
-  /**
-   * \brief A linear operator representing the sum of two linear operators.
-   *
-   * \tparam OP1 Type of the first linear operator.
-   * \tparam OP2 Type of the second linear operator.
-   *
-   * \author Sebastian Westerheide.
-   */
-  template <class OP1, class OP2>
-  class LinearOperatorSum
-    : public Dune::LinearOperator<typename OP1::domain_type,
-                                  typename OP1::range_type>
-  {
-  public:
-    typedef typename OP1::domain_type domain_type;
-    typedef typename OP1::range_type range_type;
-    typedef typename domain_type::field_type field_type;
+      virtual void apply (const X& x, Y& y) const
+      {
+        y = x;
+        y *= immutable_scaling_*mutable_scaling_;
+      }
 
-    enum {category = Dune::SolverCategory::sequential};
+      virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const
+      {
+        X temp(x);
+        temp *= immutable_scaling_*mutable_scaling_;
+        y.axpy(alpha,temp);
+      }
 
-    LinearOperatorSum (const OP1& op1, const OP2& op2)
-      : op1_(op1), op2_(op2)
-    {
-      static_assert(std::is_same<typename OP2::domain_type,domain_type>::value,
-                    "Domain type of both operators doesn't match!");
-      static_assert(std::is_same<typename OP2::range_type,range_type>::value,
-                    "Range type of both operators doesn't match!");
-    }
+    protected:
+      const field_type immutable_scaling_;
+      const field_type& mutable_scaling_;
+    };
 
-    virtual void apply (const domain_type& x, range_type& y) const
-    {
-      op1_.apply(x,y);
-      op2_.applyscaleadd(1.0,x,y);
-    }
 
-    virtual void applyscaleadd (field_type alpha,
-                                const domain_type& x, range_type& y) const
+    /**
+     * \brief A linear operator representing the sum of two linear operators.
+     *
+     * \tparam OP1 Type of the first linear operator.
+     * \tparam OP2 Type of the second linear operator.
+     *
+     * \author Sebastian Westerheide.
+     */
+    template <class OP1, class OP2>
+    class LinearOperatorSum
+      : public Dune::LinearOperator<typename OP1::domain_type,
+                                    typename OP1::range_type>
     {
-      range_type temp(y);
-      op1_.apply(x,temp);
-      op2_.applyscaleadd(1.0,x,temp);
-      y.axpy(alpha,temp);
-    }
+    public:
+      typedef typename OP1::domain_type domain_type;
+      typedef typename OP1::range_type range_type;
+      typedef typename domain_type::field_type field_type;
 
-  protected:
-    const OP1& op1_;
-    const OP2& op2_;
-  };
+      enum {category = Dune::SolverCategory::sequential};
 
+      LinearOperatorSum (const OP1& op1, const OP2& op2)
+        : op1_(op1), op2_(op2)
+      {
+        static_assert(std::is_same<typename OP2::domain_type,domain_type>::value,
+          "Domain type of both operators doesn't match!");
+        static_assert(std::is_same<typename OP2::range_type,range_type>::value,
+          "Range type of both operators doesn't match!");
+      }
 
-  /**
-   * \brief Helper class for notifying a DUNE-ISTL linear solver about
-   *        a change of the iteration matrix object in a unified way,
-   *        i.e. independent from the solver's type (direct/iterative).
-   *
-   * \author Sebastian Westerheide.
-   */
-  template <typename ISTLLinearSolver, typename BCRSMatrix>
-  class SolverHelper
-  {
-  public:
-    static void setMatrix (ISTLLinearSolver& solver,
-                           const BCRSMatrix& matrix)
-    {
-      static const bool is_direct_solver
-        = Dune::IsDirectSolver<ISTLLinearSolver>::value;
-      SolverHelper<ISTLLinearSolver,BCRSMatrix>::
-        Implementation<is_direct_solver>::setMatrix(solver,matrix);
-    }
-
-  protected:
-    /**
-     * \brief Implementation that works together with iterative ISTL
-     *        solvers, e.g. Dune::CGSolver or Dune::BiCGSTABSolver.
-     */
-    template <bool is_direct_solver, typename Dummy = void>
-    struct Implementation
-    {
-      static void setMatrix (ISTLLinearSolver&,
-                             const BCRSMatrix&)
-      {}
-    };
+      virtual void apply (const domain_type& x, range_type& y) const
+      {
+        op1_.apply(x,y);
+        op2_.applyscaleadd(1.0,x,y);
+      }
 
-    /**
-     * \brief Implementation that works together with direct ISTL
-     *        solvers, e.g. Dune::SuperLU or Dune::UMFPack.
-     */
-    template <typename Dummy>
-    struct Implementation<true,Dummy>
-    {
-      static void setMatrix (ISTLLinearSolver& solver,
-                             const BCRSMatrix& matrix)
+      virtual void applyscaleadd (field_type alpha,
+        const domain_type& x, range_type& y) const
       {
-        solver.setMatrix(matrix);
+        range_type temp(y);
+        op1_.apply(x,temp);
+        op2_.applyscaleadd(1.0,x,temp);
+        y.axpy(alpha,temp);
       }
-    };
-  };
 
+    protected:
+      const OP1& op1_;
+      const OP2& op2_;
+    };
+  } // end namespace Impl
 
   /**
-   * \brief A class template for performing some iterative eigenvalue algorithms
-   *        based on power iteration.
+   * \brief Iterative eigenvalue algorithms based on power iteration.
    *
    * Given a square matrix whose eigenvalues shall be considered, this class
    * template provides methods for performing the power iteration algorithm,
@@ -213,8 +167,8 @@ namespace Dune
     // Type definitions for type of iteration operator (m_ - mu_*I)
     typedef typename Dune::MatrixAdapter<BCRSMatrix,BlockVector,BlockVector>
       MatrixOperator;
-    typedef ScalingLinearOperator<BlockVector> ScalingOperator;
-    typedef LinearOperatorSum<MatrixOperator,ScalingOperator> OperatorSum;
+    typedef Impl::ScalingLinearOperator<BlockVector> ScalingOperator;
+    typedef Impl::LinearOperatorSum<MatrixOperator,ScalingOperator> OperatorSum;
 
   public:
     //! Type of underlying field
@@ -1078,7 +1032,8 @@ namespace Dune
     const std::string blank_;
   };
 
-}  // namespace Dune
+  /** @} */
 
+}  // namespace Dune
 
 #endif  // DUNE_ISTL_EIGENVALUE_POWERITERATION_HH

dune/istl/matrix.hh

@@ -31,7 +31,7 @@ namespace MatrixImp
        elegant solution.
    */
   template<class B, class A=std::allocator<B> >
-  class DenseMatrixBase : public block_vector_unmanaged<B,A>
+  class DenseMatrixBase : public Imp::block_vector_unmanaged<B,A>
                               // this derivation gives us all the blas level 1 and norms
                               // on the large array. However, access operators have to be
                               // overwritten.
@@ -61,7 +61,7 @@ namespace MatrixImp
     typedef BlockVector<B,A> block_type;
 
     // just a shorthand
-    typedef BlockVectorWindow<B,A> window_type;
+    typedef Imp::BlockVectorWindow<B,A> window_type;
 
     typedef window_type reference;
 
@@ -73,7 +73,7 @@ namespace MatrixImp
     /** constructor without arguments makes empty vector,
             object cannot be used yet
      */
-    DenseMatrixBase () : block_vector_unmanaged<B,A>()
+    DenseMatrixBase () : Imp::block_vector_unmanaged<B,A>()
     {
       // nothing is known ...
       rows_ = 0;
@@ -86,7 +86,7 @@ namespace MatrixImp
             \param _nblocks Number of blocks
             \param m Number of elements in each block
      */
-    DenseMatrixBase (size_type rows, size_type columns) : block_vector_unmanaged<B,A>()
+    DenseMatrixBase (size_type rows, size_type columns) : Imp::block_vector_unmanaged<B,A>()
     {
       // and we can allocate the big array in the base class
       this->n = rows*columns;
@@ -222,7 +222,7 @@ namespace MatrixImp
     //! assign from scalar
     DenseMatrixBase& operator= (const field_type& k)
     {
-      (static_cast<block_vector_unmanaged<B,A>&>(*this)) = k;
+      (static_cast<Imp::block_vector_unmanaged<B,A>&>(*this)) = k;
       return *this;
     }
 

dune/istl/paamg/amg.hh

@@ -370,6 +370,69 @@ namespace Dune
       rhs_=nullptr;
     }
 
+    template <class Matrix,
+              class Vector>
+    struct DirectSolverSelector
+    {
+      typedef typename Matrix :: field_type field_type;
+      enum SolverType { umfpack, superlu, none };
+
+      static constexpr SolverType solver =
+#if HAVE_SUITESPARSE_UMFPACK
+        UMFPackMethodChooser< field_type > :: valid ? umfpack : none ;
+#elif HAVE_SUPERLU
+        superlu ;
+#else
+        none;
+#endif
+
+      template <class M, SolverType>
+      struct Solver
+      {
+        typedef InverseOperator<Vector,Vector> type;
+        static type* create(const M& mat, bool verbose, bool reusevector )
+        {
+          DUNE_THROW(NotImplemented,"DirectSolver not selected");
+          return nullptr;
+        }
+        static std::string name () { return "None"; }
+      };
+#if HAVE_SUITESPARSE_UMFPACK
+      template <class M>
+      struct Solver< M, umfpack >
+      {
+        typedef UMFPack< M > type;
+        static type* create(const M& mat, bool verbose, bool reusevector )
+        {
+          return new type(mat, verbose, reusevector );
+        }
+        static std::string name () { return "UMFPack"; }
+      };
+#endif
+#if HAVE_SUPERLU
+      template <class M>
+      struct Solver< M, superlu >
+      {
+        typedef SuperLU< M > type;
+        static type* create(const M& mat, bool verbose, bool reusevector )
+        {
+          return new type(mat, verbose, reusevector );
+        }
+        static std::string name () { return "SuperLU"; }
+      };
+#endif
+
+      // define direct solver type to be used
+      typedef Solver< Matrix, solver > SelectedSolver ;
+      typedef typename SelectedSolver :: type   DirectSolver;
+      static constexpr bool isDirectSolver = solver != none;
+      static std::string name() { return SelectedSolver :: name (); }
+      static DirectSolver* create(const Matrix& mat, bool verbose, bool reusevector )
+      {
+        return SelectedSolver :: create( mat, verbose, reusevector );
+      }
+    };
+
     template<class M, class X, class S, class PI, class A>
     template<class C>
     void AMG<M,X,S,PI,A>::createHierarchies(C& criterion, Operator& matrix,
@@ -383,7 +446,8 @@ namespace Dune
       // build the necessary smoother hierarchies
       matrices_->coarsenSmoother(*smoothers_, smootherArgs_);
 
-            if(buildHierarchy_ && matrices_->levels()==matrices_->maxlevels()) {
+      if(buildHierarchy_ && matrices_->levels()==matrices_->maxlevels())
+      {
         // We have the carsest level. Create the coarse Solver
         SmootherArgs sargs(smootherArgs_);
         sargs.iterations = 1;
@@ -402,32 +466,35 @@ namespace Dune
         coarseSmoother_.reset(ConstructionTraits<Smoother>::construct(cargs));
         scalarProduct_.reset(ScalarProductChooser::construct(cargs.getComm()));
 
-#if HAVE_SUPERLU || HAVE_SUITESPARSE_UMFPACK
-#if HAVE_SUITESPARSE_UMFPACK
-#define DIRECTSOLVER UMFPack
-#else
-#define DIRECTSOLVER SuperLU
-#endif
+        typedef DirectSolverSelector< typename M::matrix_type, X > SolverSelector;
+
         // Use superlu if we are purely sequential or with only one processor on the coarsest level.
-        if(std::is_same<ParallelInformation,SequentialInformation>::value // sequential mode
+        if( SolverSelector::isDirectSolver &&
+            (std::is_same<ParallelInformation,SequentialInformation>::value // sequential mode
            || matrices_->parallelInformation().coarsest()->communicator().size()==1 //parallel mode and only one processor
            || (matrices_->parallelInformation().coarsest().isRedistributed()
                && matrices_->parallelInformation().coarsest().getRedistributed().communicator().size()==1
-               && matrices_->parallelInformation().coarsest().getRedistributed().communicator().size()>0)) { // redistribute and 1 proc
+               && matrices_->parallelInformation().coarsest().getRedistributed().communicator().size()>0) )
+          )
+        { // redistribute and 1 proc
           if(matrices_->parallelInformation().coarsest().isRedistributed())
           {
             if(matrices_->matrices().coarsest().getRedistributed().getmat().N()>0)
+            {
               // We are still participating on this level
-              solver_.reset(new DIRECTSOLVER<typename M::matrix_type>(matrices_->matrices().coarsest().getRedistributed().getmat(), false, false));
+              solver_.reset(SolverSelector::create(matrices_->matrices().coarsest().getRedistributed().getmat(), false, false));
+            }
             else
               solver_.reset();
-          }else
-            solver_.reset(new DIRECTSOLVER<typename M::matrix_type>(matrices_->matrices().coarsest()->getmat(), false, false));
+          }
+          else
+          {
+            solver_.reset(SolverSelector::create(matrices_->matrices().coarsest()->getmat(), false, false));
+          }
           if(verbosity_>0 && matrices_->parallelInformation().coarsest()->communicator().rank()==0)
-            std::cout<< "Using a direct coarse solver (" << static_cast< DIRECTSOLVER<typename M::matrix_type>* >(solver_.get())->name() << ")" << std::endl;
-        }else
-#undef DIRECTSOLVER
-#endif // HAVE_SUPERLU|| HAVE_SUITESPARSE_UMFPACK
+            std::cout<< "Using a direct coarse solver (" << SolverSelector::name() << ")" << std::endl;
+        }
+        else
         {
           if(matrices_->parallelInformation().coarsest().isRedistributed())
           {

dune/istl/solver.hh

@@ -6,6 +6,7 @@
 
 #include <iomanip>
 #include <ostream>
+#include "solvertype.hh"
 
 namespace Dune
 {
@@ -155,6 +156,54 @@ namespace Dune
     }
   };
 
+  /**
+   * \brief Helper class for notifying a DUNE-ISTL linear solver about
+   *        a change of the iteration matrix object in a unified way,
+   *        i.e. independent from the solver's type (direct/iterative).
+   *
+   * \author Sebastian Westerheide.
+   */
+  template <typename ISTLLinearSolver, typename BCRSMatrix>
+  class SolverHelper
+  {
+  public:
+    static void setMatrix (ISTLLinearSolver& solver,
+                           const BCRSMatrix& matrix)
+    {
+      static const bool is_direct_solver
+        = Dune::IsDirectSolver<ISTLLinearSolver>::value;
+      SolverHelper<ISTLLinearSolver,BCRSMatrix>::
+        Implementation<is_direct_solver>::setMatrix(solver,matrix);
+    }
+
+  protected:
+    /**
+     * \brief Implementation that works together with iterative ISTL
+     *        solvers, e.g. Dune::CGSolver or Dune::BiCGSTABSolver.
+     */
+    template <bool is_direct_solver, typename Dummy = void>
+    struct Implementation
+    {
+      static void setMatrix (ISTLLinearSolver&,
+                             const BCRSMatrix&)
+      {}
+    };
+
+    /**
+     * \brief Implementation that works together with direct ISTL
+     *        solvers, e.g. Dune::SuperLU or Dune::UMFPack.
+     */
+    template <typename Dummy>
+    struct Implementation<true,Dummy>
+    {
+      static void setMatrix (ISTLLinearSolver& solver,
+                             const BCRSMatrix& matrix)
+      {
+        solver.setMatrix(matrix);
+      }
+    };
+  };
+
 /**
  * @}
  */

dune/istl/solvers.hh

@@ -25,9 +25,6 @@
 #include <dune/common/typetraits.hh>
 
 namespace Dune {
-  /** @defgroup ISTL_Solvers Iterative Solvers
-      @ingroup ISTL
-   */
   /** @addtogroup ISTL_Solvers
       @{
    */

dune/istl/superlu.hh

@@ -4,32 +4,8 @@
 #define DUNE_ISTL_SUPERLU_HH
 
 #if HAVE_SUPERLU
-#ifdef TRUE
-#undef TRUE
-#endif
-#ifdef FALSE
-#undef FALSE
-#endif
-#ifndef SUPERLU_NTYPE
-#define  SUPERLU_NTYPE 1
-#endif
-
-#if SUPERLU_NTYPE==0
-#include "slu_sdefs.h"
-#endif
-
-#if SUPERLU_NTYPE==1
-#include "slu_ddefs.h"
-#endif
-
-#if SUPERLU_NTYPE==2
-#include "slu_cdefs.h"
-#endif
-
-#if SUPERLU_NTYPE>=3
-#include "slu_zdefs.h"
-#endif
 
+#include "superlufunctions.hh"
 #include "solvers.hh"
 #include "supermatrix.hh"
 #include <algorithm>
@@ -81,7 +57,7 @@ namespace Dune
   struct QuerySpaceChooser
   {};
 
-#if SUPERLU_NTYPE==0
+#if HAVE_SLU_SDEFS_H
   template<>
   struct SuperLUDenseMatChooser<float>
   {
@@ -129,7 +105,7 @@ namespace Dune
 
 #endif
 
-#if SUPERLU_NTYPE==1
+#if HAVE_SLU_DDEFS_H
 
   template<>
   struct SuperLUDenseMatChooser<double>
@@ -176,7 +152,7 @@ namespace Dune
   };
 #endif
 
-#if SUPERLU_NTYPE>=3
+#if HAVE_SLU_ZDEFS_H
   template<>
   struct SuperLUDenseMatChooser<std::complex<double> >
   {
@@ -223,7 +199,7 @@ namespace Dune
   };
 #endif
 
-#if SUPERLU_NTYPE==2
+#if HAVE_SLU_CDEFS_H
   template<>
   struct SuperLUDenseMatChooser<std::complex<float> >
   {

dune/istl/superlufunctions.hh

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_ISTL_SUPERLUFUNCTIONS_HH
+#define DUNE_ISTL_SUPERLUFUNCTIONS_HH
+#if HAVE_SUPERLU
+
+
+#define int_t SUPERLU_INT_TYPE
+#include "supermatrix.h"
+#include "slu_util.h"
+#undef int_t
+
+#if HAVE_SLU_SDEFS_H
+extern "C" {
+  extern void
+  sgssvx(superlu_options_t *, SuperMatrix *, int *, int *, int *,
+         char *, float *, float *, SuperMatrix *, SuperMatrix *,
+         void *, int, SuperMatrix *, SuperMatrix *,
+         float *, float *, float *, float *,
+#if SUPERLU_MIN_VERSION_5
+         GlobalLU_t*,
+#endif
+         mem_usage_t *, SuperLUStat_t *, int *);
+
+  extern void
+  sCreate_Dense_Matrix(SuperMatrix *, int, int, float *, int,
+                       Stype_t, Dtype_t, Mtype_t);
+  extern void
+  sCreate_CompCol_Matrix(SuperMatrix *, int, int, int, float *,
+                         int *, int *, Stype_t, Dtype_t, Mtype_t);
+  extern int     sQuerySpace (SuperMatrix *, SuperMatrix *, mem_usage_t *);
+
+  extern void    sPrint_CompCol_Matrix(char *, SuperMatrix *);
+}
+#endif
+
+#if HAVE_SLU_DDEFS_H
+extern "C" {
+  extern void
+  dgssvx(superlu_options_t *, SuperMatrix *, int *, int *, int *,
+         char *, double *, double *, SuperMatrix *, SuperMatrix *,
+         void *, int, SuperMatrix *, SuperMatrix *,
+         double *, double *, double *, double *,
+#if SUPERLU_MIN_VERSION_5
+         GlobalLU_t*,
+#endif
+         mem_usage_t *, SuperLUStat_t *, int *);
+
+  extern void
+  dCreate_CompCol_Matrix(SuperMatrix *, int, int, int, double *,
+                         int *, int *, Stype_t, Dtype_t, Mtype_t);
+
+  extern void
+  dCreate_Dense_Matrix(SuperMatrix *, int, int, double *, int,
+                       Stype_t, Dtype_t, Mtype_t);
+
+  extern int     dQuerySpace (SuperMatrix *, SuperMatrix *, mem_usage_t *);
+
+  extern void    dPrint_CompCol_Matrix(char *, SuperMatrix *);
+}
+#endif
+
+#if HAVE_SLU_CDEFS_H
+#include "slu_scomplex.h"
+
+extern "C" {
+  extern void
+  cgssvx(superlu_options_t *, SuperMatrix *, int *, int *, int *,
+         char *, float *, float *, SuperMatrix *, SuperMatrix *,
+         void *, int, SuperMatrix *, SuperMatrix *,
+         float *, float *, float *, float *,
+#if SUPERLU_MIN_VERSION_5
+         GlobalLU_t*,
+#endif
+         mem_usage_t *, SuperLUStat_t *, int *);
+
+
+  extern void
+  cCreate_Dense_Matrix(SuperMatrix *, int, int, ::complex *, int,
+                       Stype_t, Dtype_t, Mtype_t);
+
+
+  extern void
+  cCreate_CompCol_Matrix(SuperMatrix *, int, int, int, ::complex *,
+                         int *, int *, Stype_t, Dtype_t, Mtype_t);
+
+  extern int     cQuerySpace (SuperMatrix *, SuperMatrix *, mem_usage_t *);
+
+  extern void    cPrint_CompCol_Matrix(char *, SuperMatrix *);
+}
+#endif
+
+#if HAVE_SLU_ZDEFS_H
+#include "slu_dcomplex.h"
+extern "C" {
+  extern void
+  zgssvx(superlu_options_t *, SuperMatrix *, int *, int *, int *,
+         char *, double *, double *, SuperMatrix *, SuperMatrix *,
+         void *, int, SuperMatrix *, SuperMatrix *,
+         double *, double *, double *, double *,
+#if SUPERLU_MIN_VERSION_5
+         GlobalLU_t*,
+#endif
+         mem_usage_t *, SuperLUStat_t *, int *);
+
+
+  extern void
+  zCreate_CompCol_Matrix(SuperMatrix *, int, int, int, doublecomplex *,
+                         int *, int *, Stype_t, Dtype_t, Mtype_t);
+
+  extern void
+  zCreate_Dense_Matrix(SuperMatrix *, int, int, doublecomplex *, int,
+                       Stype_t, Dtype_t, Mtype_t);
+
+  extern int     zQuerySpace (SuperMatrix *, SuperMatrix *, mem_usage_t *);
+
+  extern void    zPrint_CompCol_Matrix(char *, SuperMatrix *);
+}
+#endif
+
+
+#endif
+#endif

dune/istl/supermatrix.hh

@@ -5,26 +5,6 @@
 
 #if HAVE_SUPERLU
 
-#ifndef SUPERLU_NTYPE
-#define SUPERLU_NTYPE 1
-#endif
-
-#if SUPERLU_NTYPE==0
-#include "slu_sdefs.h"
-#endif
-
-#if SUPERLU_NTYPE==1
-#include "slu_ddefs.h"
-#endif
-
-#if SUPERLU_NTYPE==2
-#include "slu_cdefs.h"
-#endif
-
-#if SUPERLU_NTYPE>=3
-#include "slu_zdefs.h"
-#endif
-
 #include "bcrsmatrix.hh"
 #include "bvector.hh"
 #include <dune/common/fmatrix.hh>
@@ -34,6 +14,8 @@
 
 #include"colcompmatrix.hh"
 
+#include "superlufunctions.hh"
+
 namespace Dune
 {
 
@@ -45,7 +27,7 @@ namespace Dune
   struct SuperMatrixPrinter
   {};
 
-#if SUPERLU_NTYPE==0
+#if HAVE_SLU_SDEFS_H
   template<>
   struct SuperMatrixCreateSparseChooser<float>
   {
@@ -68,7 +50,7 @@ namespace Dune
   };
 #endif
 
-#if SUPERLU_NTYPE==1
+#if HAVE_SLU_DDEFS_H
   template<>
   struct SuperMatrixCreateSparseChooser<double>
   {
@@ -91,7 +73,7 @@ namespace Dune
   };
 #endif
 
-#if SUPERLU_NTYPE==2
+#if HAVE_SLU_CDEFS_H
   template<>
   struct SuperMatrixCreateSparseChooser<std::complex<float> >
   {
@@ -114,7 +96,7 @@ namespace Dune
   };
 #endif
 
-#if SUPERLU_NTYPE>=3
+#if HAVE_SLU_ZDEFS_H
   template<>
   struct SuperMatrixCreateSparseChooser<std::complex<double> >
   {

dune/istl/test/CMakeLists.txt

@@ -45,6 +45,8 @@ dune_add_test(SOURCES inverseoperator2prectest.cc)
 
 dune_add_test(SOURCES scaledidmatrixtest.cc)
 
+dune_add_test(SOURCES solvertest.cc)
+
 dune_add_test(SOURCES solveraborttest.cc)
 
 # Pardiso tests

dune/istl/test/basearraytest.cc

@@ -8,8 +8,8 @@ using namespace Dune;
 
 int main()
 {
-  base_array<double> v1(10);
-  base_array<double> v2 = v1;
+  Imp::base_array<double> v1(10);
+  Imp::base_array<double> v2 = v1;
 
   v1.resize(20);
 

dune/istl/test/superlutest.cc

@@ -14,36 +14,11 @@
 
 #include "laplacian.hh"
 
-#ifndef SUPERLU_NTYPE
-#define SUPERLU_NTYPE 1
-#endif
-
-#if SUPERLU_NTYPE==1
-typedef double FIELD_TYPE;
-#endif
-
-#if SUPERLU_NTYPE==0
-typedef float FIELD_TYPE;
-#endif
-
-#if SUPERLU_NTYPE==2
-typedef std::complex<float> FIELD_TYPE;
-#endif
-
-#if SUPERLU_NTYPE>=3
-typedef std::complex<double> FIELD_TYPE;
-#endif
-
-int main(int argc, char** argv)
-try
-{
 #if HAVE_SUPERLU
+template<typename FIELD_TYPE>
+void runSuperLU(std::size_t N)
+{
   const int BS=1;
-  std::size_t N=100;
-
-  if(argc>1)
-    N = atoi(argv[1]);
-  std::cout<<"testing for N="<<N<<" BS="<<1<<std::endl;
 
 
   typedef Dune::FieldMatrix<FIELD_TYPE,BS,BS> MatrixBlock;
@@ -83,8 +58,34 @@ try
 
   solver1.apply(x1,b1, res);
   solver1.apply(reinterpret_cast<FIELD_TYPE*>(&x1[0]), reinterpret_cast<FIELD_TYPE*>(&b1[0]));
+}
+#endif
+
+int main(int argc, char** argv)
+try
+{
+
+  std::size_t N=100;
+
+  if(argc>1)
+    N = atoi(argv[1]);
+  std::cout<<"testing for N="<<N<<" BS="<<1<<std::endl;
+#if HAVE_SUPERLU
+#if HAVE_SLU_SDEFS_H
+  runSuperLU<float>(N);
+#endif
+#if HAVE_SLU_DDEFS_H
+  runSuperLU<double>(N);
+#endif
+#if HAVE_SLU_CDEFS_H
+  runSuperLU<std::complex<float> >(N);
+#endif
+#if HAVE_SLU_ZDEFS_H
+  runSuperLU<std::complex<double> >(N);
+#endif
 
   return 0;
+
 #else // HAVE_SUPERLU
   std::cerr << "You need SuperLU to run this test." << std::endl;
   return 77;