.gitlab-ci.yml

@@ -2,26 +2,39 @@
 before_script:
   - duneci-install-module https://gitlab.dune-project.org/core/dune-common.git
 
+variables:
+# Suitesparse, as installed with Debian, is thread-parallel using OpenMP.
+# OpenMP silently assumes, it can spawn as many threads as there are cores.
+# In a worst case scenario, this leads to a number of threads quadratic in
+# the core count, if you also do parallel test execution with the maximum
+# number of cores. We solve the issue by disallowing OpenMP to allocate more
+# than one thread.
+  OMP_NUM_THREADS: 1
+
+debian:10  gcc:c++17:
+  image: duneci/base:10
+  script: duneci-standard-test
+  variables: {DUNECI_OPTS: /duneci/opts.gcc.c++17}
+  tags: [duneci]
+
 debian:9--gcc:
   image: duneci/base:9
   script: duneci-standard-test
+  tags: [duneci]
 
 debian:9--clang:
   image: duneci/base:9
-  script: duneci-standard-test --opts=/duneci/opts.clang
-
-debian:8--gcc:
-  image: duneci/base:8
   script: duneci-standard-test
-
-debian:8-backports--clang:
-  image: duneci/base:8-backports
-  script: duneci-standard-test --opts=/duneci/opts.clang
+  variables: {DUNECI_OPTS: /duneci/opts.clang}
+  tags: [duneci]
 
 ubuntu:16.04--gcc:
   image: duneci/base:16.04
   script: duneci-standard-test
+  tags: [duneci]
 
 ubuntu:16.04--clang:
   image: duneci/base:16.04
-  script: duneci-standard-test --opts=/duneci/opts.clang
+  script: duneci-standard-test
+  variables: {DUNECI_OPTS: /duneci/opts.clang}
+  tags: [duneci]

CHANGELOG.md

+# Master (will become release 2.7)
+
+- Deprecated the preconditioner implementations `SeqILU0` and `SeqILUn`.
+  Use `SeqILU` instead, which implements incomplete LU decomposition
+  of any order.
+
+# Release 2.6
+
+- `BDMatrix` objects can now be constructed and assigned from `std::initializer_list`.
+
+- `BDMatrix` and `BTDMatrix` now implement the `setSize` method, which allows to
+  resize existing matrix objects.
+
+- The solver infrastructure was updated to support SIMD data types (see
+  current changes in `dune-common`). This allows to solve multiple systems
+  simultaniously using vectorization.

CMakeLists.txt

@@ -2,7 +2,7 @@
 project(dune-istl C CXX)
 
 # general stuff
-cmake_minimum_required(VERSION 2.8.12)
+cmake_minimum_required(VERSION 3.1)
 
 # guess build tree of dune-common
 if(NOT (dune-common_DIR OR dune-common_ROOT OR

cmake/modules/AddARPACKPPFlags.cmake

@@ -14,10 +14,8 @@ function(add_dune_arpackpp_flags _targets)
   if(ARPACKPP_FOUND)
     foreach(_target ${_targets})
       target_link_libraries(${_target} ${ARPACKPP_DUNE_LIBRARIES})
-      get_target_property(_props ${_target} COMPILE_FLAGS)
-      string(REPLACE "_props-NOTFOUND" "" _props "${_props}")
-      set_target_properties(${_target} PROPERTIES COMPILE_FLAGS
-        "${_props} ${ARPACKPP_DUNE_COMPILE_FLAGS} -DENABLE_ARPACKPP=1")
+      target_compile_definitions(${_target} PUBLIC ENABLE_ARPACKPP=1)
+      target_compile_options(${_target} PUBLIC ${ARPACKPP_DUNE_COMPILE_FLAGS})
     endforeach()
   endif()
 endfunction(add_dune_arpackpp_flags)

cmake/modules/DuneIstlMacros.cmake

@@ -12,3 +12,7 @@ find_package(ARPACKPP)
 include(AddARPACKPPFlags)
 find_package(SuiteSparse OPTIONAL_COMPONENTS LDL SPQR UMFPACK)
 include(AddSuiteSparseFlags)
+
+# enable / disable backwards compatibility w.r.t. category
+set(DUNE_ISTL_SUPPORT_OLD_CATEGORY_INTERFACE 1
+  CACHE BOOL "Enable/Disable the backwards compatibility of the category enum/method in dune-istl solvers, preconditioner, etc. '1'")

cmake/modules/FindARPACK.cmake

@@ -88,3 +88,8 @@ else()
     "Determing location of ARPACK failed:\n"
     "Libraries to link against: ${ARPACK_LIBRARIES}\n\n")
 endif()
+
+# text for feature summary
+set_package_properties("ARPACK" PROPERTIES
+  DESCRIPTION "ARnoldi PACKage"
+  PURPOSE "Solve large scale eigenvalue problems")

cmake/modules/FindARPACKPP.cmake

@@ -100,7 +100,15 @@ if(ARPACKPP_FOUND)
     "Determing location of ARPACK++ succeeded:\n"
     "Include directory: ${ARPACKPP_INCLUDE_DIRS}\n"
     "Libraries to link against: ${ARPACKPP_LIBRARIES}\n\n")
-  set(ARPACKPP_DUNE_COMPILE_FLAGS "-I${ARPACKPP_INCLUDE_DIRS}"
+
+  # the following is a pretty roundabout way of setting include directories, but it's the
+  # only way we can force -isystem. And we want the compiler to treat ARPACK++ as a system
+  # library to avoid scaring our users with the horrible warnings triggered by the bitrotted
+  # ARPACK++ sources.
+  #
+  # For this to work, only set the COMPILE_OPTIONS (those replaced COMPILE_FLAGS a while ago), never
+  # the INCLUDE_DIRECTORIES!
+  set(ARPACKPP_DUNE_COMPILE_FLAGS "$<$<BOOL:${ARPACKPP_INCLUDE_DIRS}>:-isystem$<JOIN:${ARPACKPP_INCLUDE_DIRS}, -isystem>>"
     CACHE STRING "Compile flags used by DUNE when compiling ARPACK++ programs")
   set(ARPACKPP_DUNE_LIBRARIES ${ARPACKPP_LIBRARIES}
     CACHE STRING "Libraries used by DUNE when linking ARPACK++ programs")
@@ -119,5 +127,10 @@ set(HAVE_ARPACKPP ${ARPACKPP_FOUND})
 if(ARPACKPP_FOUND)
   dune_register_package_flags(COMPILE_DEFINITIONS "ENABLE_ARPACKPP=1"
                               LIBRARIES "${ARPACKPP_LIBRARIES}"
-                              INCLUDE_DIRS "${ARPACKPP_INCLUDE_DIRS}")
+                              COMPILE_OPTIONS "${ARPACKPP_DUNE_COMPILE_FLAGS}")
 endif()
+
+# text for feature summary
+set_package_properties("ARPACKPP" PROPERTIES
+  DESCRIPTION "ARPACK++"
+  PURPOSE "C++ interface for ARPACK")

cmake/modules/FindSuperLU.cmake

@@ -168,8 +168,6 @@ find_package_handle_standard_args(
 
 mark_as_advanced(SUPERLU_INCLUDE_DIR SUPERLU_LIBRARY)
 
-set_package_info("SuperLU" "Direct linear solver library")
-
 # if both headers and library are found, store results
 if(SUPERLU_FOUND)
   set(SUPERLU_INCLUDE_DIRS ${SUPERLU_INCLUDE_DIR})
@@ -201,3 +199,8 @@ if(SUPERLU_FOUND)
                               LIBRARIES "${SUPERLU_DUNE_LIBRARIES}"
                               INCLUDE_DIRS "${SUPERLU_INCLUDE_DIRS}")
 endif()
+
+# text for feature summary
+set_package_properties("SuperLU" PROPERTIES
+  DESCRIPTION "Supernodal LU"
+  PURPOSE "Direct solver for linear system, based on LU decomposition")

config.h.cmake

@@ -71,6 +71,9 @@
 /* Define to the revision of dune-istl */
 #define DUNE_ISTL_VERSION_REVISION ${DUNE_ISTL_VERSION_REVISION}
 
+/* Enable/Disable the backwards compatibility of the category enum/method in dune-istl solvers, preconditioner, etc. */
+#cmakedefine DUNE_ISTL_SUPPORT_OLD_CATEGORY_INTERFACE @DUNE_ISTL_SUPPORT_OLD_CATEGORY_INTERFACE@
+
 /* end dune-istl
    Everything below here will be overwritten
 */

dune.module

 Module: dune-istl
-Version: 2.6-git
-Maintainer: dune-devel@dune-project.org
+Version: 2.7-git
+Maintainer: dune-devel@lists.dune-project.org
 Depends: dune-common (>= 2.6)
 Whitespace-Hook: Yes

dune/istl/CMakeLists.txt

 add_subdirectory("eigenvalue")
 add_subdirectory("paamg")
-add_subdirectory("tutorial")
 add_subdirectory(test)
 
 #install headers
 install(FILES
+   allocator.hh
    basearray.hh
    bcrsmatrix.hh
    bdmatrix.hh
@@ -12,6 +12,7 @@ install(FILES
    bvector.hh
    colcompmatrix.hh
    gsetc.hh
+   ildl.hh
    ilu.hh
    ilusubdomainsolver.hh
    io.hh

dune/istl/allocator.hh

+#ifndef DUNE_ISTL_ALLOCATOR_HH
+#define DUNE_ISTL_ALLOCATOR_HH
+
+#include <dune/common/typetraits.hh>
+#include <memory>
+
+namespace Dune {
+
+    template<typename T>
+    struct exists{
+        static const bool value = true;
+    };
+
+    template<typename T, typename = void>
+    struct DefaultAllocatorTraits
+    {
+        using type = std::allocator<T>;
+    };
+
+    template<typename T>
+    struct DefaultAllocatorTraits<T, void_t<typename T::allocator_type> >
+    {
+        using type = typename T::allocator_type;
+    };
+
+    template<typename T>
+    struct AllocatorTraits : public DefaultAllocatorTraits<T> {};
+
+    template<typename T>
+    using AllocatorType = typename AllocatorTraits<T>::type;
+
+    template<typename T, typename X>
+    using ReboundAllocatorType = typename AllocatorTraits<T>::type::template rebind<X>::other;
+
+} // end namespace Dune
+
+#endif // DUNE_ISTL_ALLOCATOR_HH

dune/istl/basearray.hh

@@ -61,11 +61,16 @@ namespace Imp {
     //! the type for the index access
     typedef typename A::size_type size_type;
 
+    //! the type used for references
+    using reference = B&;
+
+    //! the type used for const references
+    using const_reference = const B&;
 
     //===== access to components
 
     //! random access to blocks
-    B& operator[] (size_type i)
+    reference operator[] (size_type i)
     {
 #ifdef DUNE_ISTL_WITH_CHECKING
       if (i>=n) DUNE_THROW(ISTLError,"index out of range");
@@ -74,7 +79,7 @@ namespace Imp {
     }
 
     //! same for read only access
-    const B& operator[] (size_type i) const
+    const_reference operator[] (size_type i) const
     {
 #ifdef DUNE_ISTL_WITH_CHECKING
       if (i>=n) DUNE_THROW(ISTLError,"index out of range");
@@ -147,13 +152,13 @@ namespace Imp {
       }
 
       // Needed for operator[] of the iterator
-      B& elementAt (std::ptrdiff_t offset) const
+      reference elementAt (std::ptrdiff_t offset) const
       {
         return *(i+offset);
       }
 
       //! dereferencing
-      B& dereference () const
+      reference dereference () const
       {
         return *i;
       }
@@ -247,6 +252,18 @@ namespace Imp {
       return n;
     }
 
+    //! Returns pointer to the underlying array
+    const B* data() const
+    {
+      return p;
+    }
+
+    //! Returns pointer to the underlying array
+    B* data()
+    {
+      return p;
+    }
+
   protected:
     //! makes empty array
     base_array_unmanaged ()
@@ -262,248 +279,6 @@ namespace Imp {
 
 
 
-  /**  \brief Extend base_array_unmanaged by functions to manipulate
-
-           This container has *NO* memory management at all,
-           copy constuctor, assignment and destructor are all default.
-           A container can be constructed as empty or from a given pointer
-       and size. This class is used to implement a view into a larger
-       array.
-
-           You can copy such an object to a base_array to make a real copy.
-
-           Error checking: no error checking is provided normally.
-           Setting the compile time switch DUNE_ISTL_WITH_CHECKING
-           enables error checking.
-   \internal This class is an implementation detail, and should not be used outside of dune-istl.
-   */
-  template<class B, class A=std::allocator<B> >
-  class base_array_window : public base_array_unmanaged<B,A>
-  {
-  public:
-
-    //===== type definitions and constants
-
-    //! export the type representing the components
-    typedef B member_type;
-
-    //! export the allocator type
-    typedef A allocator_type;
-
-    //! make iterators available as types
-    typedef typename base_array_unmanaged<B,A>::iterator iterator;
-
-    //! make iterators available as types
-    typedef typename base_array_unmanaged<B,A>::const_iterator const_iterator;
-
-    //! The type used for the index access
-    typedef typename base_array_unmanaged<B,A>::size_type size_type;
-
-    //! The type used for the difference between two iterator positions
-    typedef typename A::difference_type difference_type;
-
-    //===== constructors and such
-
-    //! makes empty array
-    base_array_window ()
-      : base_array_unmanaged<B,A>()
-    {   }
-
-    //! make array from given pointer and size
-    base_array_window (B* _p, size_type _n)
-      : base_array_unmanaged<B,A>(_n ,_p)
-    {}
-
-    //===== window manipulation methods
-
-    //! set pointer and length
-    void set (size_type _n, B* _p)
-    {
-      this->n = _n;
-      this->p = _p;
-    }
-
-    //! advance pointer by newsize elements and then set size to new size
-    void advance (difference_type newsize)
-    {
-      this->p += this->n;
-      this->n = newsize;
-    }
-
-    //! increment pointer by offset and set size
-    void move (difference_type offset, size_type newsize)
-    {
-      this->p += offset;
-      this->n = newsize;
-    }
-
-    //! increment pointer by offset, leave size
-    void move (difference_type offset)
-    {
-      this->p += offset;
-    }
-
-    //! return the pointer
-    B* getptr ()
-    {
-      return this->p;
-    }
-  };
-
-
-
-  /**  \brief  This container extends base_array_unmanaged by memory management
-        with the usual copy semantics providing the full range of
-        copy constructor, destructor and assignment operators.
-
-            You can make
-
-        - empty array
-        - array with n components dynamically allocated
-        - resize an array with complete loss of data
-        - assign/construct from a base_array_window to make a copy of the data
-
-            Error checking: no error checking is provided normally.
-            Setting the compile time switch DUNE_ISTL_WITH_CHECKING
-            enables error checking.
-   \internal This class is an implementation detail, and should not be used outside of dune-istl.
-   */
-  template<class B, class A=std::allocator<B> >
-  class base_array : public base_array_unmanaged<B,A>
-  {
-  public:
-
-    //===== type definitions and constants
-
-    //! export the type representing the components
-    typedef B member_type;
-
-    //! export the allocator type
-    typedef A allocator_type;
-
-    //! make iterators available as types
-    typedef typename base_array_unmanaged<B,A>::iterator iterator;
-
-    //! make iterators available as types
-    typedef typename base_array_unmanaged<B,A>::const_iterator const_iterator;
-
-    //! The type used for the index access
-    typedef typename base_array_unmanaged<B,A>::size_type size_type;
-
-    //! The type used for the difference between two iterator positions
-    typedef typename A::difference_type difference_type;
-
-    //===== constructors and such
-
-    //! makes empty array
-    base_array ()
-      : base_array_unmanaged<B,A>()
-    {}
-
-    //! make array with _n components
-    base_array (size_type _n)
-      : base_array_unmanaged<B,A>(_n, 0)
-    {
-      if (this->n>0) {
-        this->p = allocator_.allocate(this->n);
-        new (this->p)B[this->n];
-      } else
-      {
-        this->n = 0;
-        this->p = 0;
-      }
-    }
-
-    //! copy constructor
-    base_array (const base_array& a)
-    {
-      // allocate memory with same size as a
-      this->n = a.n;
-
-      if (this->n>0) {
-        this->p = allocator_.allocate(this->n);
-        new (this->p)B[this->n];
-      } else
-      {
-        this->n = 0;
-        this->p = 0;
-      }
-
-      // and copy elements
-      for (size_type i=0; i<this->n; i++) this->p[i]=a.p[i];
-    }
-
-    //! free dynamic memory
-    ~base_array ()
-    {
-      if (this->n>0) {
-        int i=this->n;
-        while (i)
-          this->p[--i].~B();
-        allocator_.deallocate(this->p,this->n);
-      }
-    }
-
-    //! reallocate array to given size, any data is lost
-    void resize (size_type _n)
-    {
-      if (this->n==_n) return;
-
-      if (this->n>0) {
-        int i=this->n;
-        while (i)
-          this->p[--i].~B();
-        allocator_.deallocate(this->p,this->n);
-      }
-      this->n = _n;
-      if (this->n>0) {
-        this->p = allocator_.allocate(this->n);
-        new (this->p)B[this->n];
-      } else
-      {
-        this->n = 0;
-        this->p = 0;
-      }
-    }
-
-    //! assignment
-    base_array& operator= (const base_array& a)
-    {
-      if (&a!=this)     // check if this and a are different objects
-      {
-        // adjust size of array
-        if (this->n!=a.n)           // check if size is different
-        {
-          if (this->n>0) {
-            int i=this->n;
-            while (i)
-              this->p[--i].~B();
-            allocator_.deallocate(this->p,this->n);                     // delete old memory
-          }
-          this->n = a.n;
-          if (this->n>0) {
-            this->p = allocator_.allocate(this->n);
-            new (this->p)B[this->n];
-          } else
-          {
-            this->n = 0;
-            this->p = 0;
-          }
-        }
-        // copy data
-        for (size_type i=0; i<this->n; i++) this->p[i]=a.p[i];
-      }
-      return *this;
-    }
-
-  protected:
-
-    A allocator_;
-  };
-
-
-
-
   /** \brief A simple array container with non-consecutive index set.
 
        Elements in the array are of type B. This class provides
@@ -541,10 +316,16 @@ namespace Imp {
     //! The type used for the index access
     typedef typename A::size_type size_type;
 
+    //! the type used for references
+    using reference = B&;
+
+    //! the type used for const references
+    using const_reference = const B&;
+
     //===== access to components
 
     //! random access to blocks, assumes ascending ordering
-    B& operator[] (size_type i)
+    reference operator[] (size_type i)
     {
       const size_type* lb = std::lower_bound(j, j+n, i);
       if (lb == j+n || *lb != i)
@@ -553,7 +334,7 @@ namespace Imp {
     }
 
     //! same for read only access, assumes ascending ordering
-    const B& operator[] (size_type i) const
+    const_reference operator[] (size_type i) const
     {
       const size_type* lb = std::lower_bound(j, j+n, i);
       if (lb == j+n || *lb != i)
@@ -646,7 +427,7 @@ namespace Imp {
       }
 
       //! dereferencing
-      B& dereference () const
+      reference dereference () const
       {
         return p[i];
       }

dune/istl/bcrsmatrix.hh

@@ -770,15 +770,7 @@ namespace Dune {
         DUNE_THROW(InvalidStateException,"BCRSMatrix can only be copy-constructed when source matrix is completely empty (size not set) or fully built)");
 
       // deep copy in global array
-      size_type _nnz = Mat.nnz_;
-
-      // in case of row-wise allocation
-      if (_nnz<=0)
-      {
-        _nnz = 0;
-        for (size_type i=0; i<Mat.n; i++)
-          _nnz += Mat.r[i].getsize();
-      }
+      size_type _nnz = Mat.nonzeroes();
 
       j_ = Mat.j_; // enable column index sharing, release array in case of row-wise allocation
       allocate(Mat.n, Mat.m, _nnz, true, true);
@@ -895,12 +887,7 @@ namespace Dune {
         rowAllocator_.deallocate(r,n);
       }
 
-      nnz_ = Mat.nnz_;
-      if (nnz_ <= 0)
-      {
-        for (size_type i=0; i<Mat.n; i++)
-          nnz_ += Mat.r[i].getsize();
-      }
+      nnz_ = Mat.nonzeroes();
 
       // allocate a, share j_
       j_ = Mat.j_;
@@ -940,6 +927,9 @@ namespace Dune {
           DUNE_THROW(BCRSMatrixError,"creation only allowed for uninitialized matrix");
         if(Mat.build_mode!=row_wise)
           DUNE_THROW(BCRSMatrixError,"creation only allowed if row wise allocation was requested in the constructor");
+        if(i==0 && _Mat.N()==0)
+          // empty Matrix is always built.
+           Mat.ready = built;
       }
 
       //! prefix increment
@@ -1234,6 +1224,7 @@ namespace Dune {
           if (j.index() >= m) {
             dwarn << "WARNING: size of row "<< i.index()<<" is "<<j.offset()<<". But was specified as being "<< (*i).end().offset()
                   <<". This means you are wasting valuable space and creating additional cache misses!"<<std::endl;
+            nnz_ -= ((*i).end().offset() - j.offset());
             r[i.index()].setsize(j.offset());
             break;
           }
@@ -1813,7 +1804,7 @@ namespace Dune {
 
     //! infinity norm (row sum norm, how to generalize for blocks?)
     template <typename ft = field_type,
-              typename std::enable_if<!has_nan<ft>::value, int>::type = 0>
+              typename std::enable_if<!HasNaN<ft>::value, int>::type = 0>
     typename FieldTraits<ft>::real_type infinity_norm() const {
       if (ready != built)
         DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
@@ -1833,7 +1824,7 @@ namespace Dune {
 
     //! simplified infinity norm (uses Manhattan norm for complex values)
     template <typename ft = field_type,
-              typename std::enable_if<!has_nan<ft>::value, int>::type = 0>
+              typename std::enable_if<!HasNaN<ft>::value, int>::type = 0>
     typename FieldTraits<ft>::real_type infinity_norm_real() const {
       if (ready != built)
         DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
@@ -1853,7 +1844,7 @@ namespace Dune {
 
     //! infinity norm (row sum norm, how to generalize for blocks?)
     template <typename ft = field_type,
-              typename std::enable_if<has_nan<ft>::value, int>::type = 0>
+              typename std::enable_if<HasNaN<ft>::value, int>::type = 0>
     typename FieldTraits<ft>::real_type infinity_norm() const {
       if (ready != built)
         DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
@@ -1876,7 +1867,7 @@ namespace Dune {
 
     //! simplified infinity norm (uses Manhattan norm for complex values)
     template <typename ft = field_type,
-              typename std::enable_if<has_nan<ft>::value, int>::type = 0>
+              typename std::enable_if<HasNaN<ft>::value, int>::type = 0>
     typename FieldTraits<ft>::real_type infinity_norm_real() const {
       if (ready != built)
         DUNE_THROW(BCRSMatrixError,"You can only call arithmetic operations on fully built BCRSMatrix instances");
@@ -1914,6 +1905,9 @@ namespace Dune {
     //! number of blocks that are stored (the number of blocks that possibly are nonzero)
     size_type nonzeroes () const
     {
+      // in case of row-wise allocation
+      if( nnz_ <= 0 )
+        nnz_ = std::accumulate( begin(), end(), size_type( 0 ), [] ( size_type s, const row_type &row ) { return s+row.getsize(); } );
       return nnz_;
     }
 
@@ -1957,7 +1951,7 @@ namespace Dune {
     // size of the matrix
     size_type n;       // number of rows
     size_type m;       // number of columns
-    size_type nnz_;     // number of nonzeroes contained in the matrix
+    mutable size_type nnz_;     // number of nonzeroes contained in the matrix
     size_type allocationSize_; //allocated size of a and j arrays, except in implicit mode: nnz_==allocationsSize_
     // zero means that memory is allocated separately for each row.
 
@@ -2158,6 +2152,9 @@ namespace Dune {
           if (r)
             DUNE_THROW(InvalidStateException,"Rows have already been allocated, cannot allocate a second time");
           r = rowAllocator_.allocate(rows);
+          // initialize row entries
+          for(row_type* ri=r; ri!=r+rows; ++ri)
+            rowAllocator_.construct(ri, row_type());
         }else{
           r = 0;
         }
@@ -2172,8 +2169,6 @@ namespace Dune {
           j_.reset(sizeAllocator_.allocate(allocationSize_),Deallocator(sizeAllocator_));
       }else{
         j_.reset();
-        for(row_type* ri=r; ri!=r+rows; ++ri)
-          rowAllocator_.construct(ri, row_type());
       }
 
       // Mark the matrix as not built.

dune/istl/bdmatrix.hh

@@ -5,6 +5,8 @@
 
 #include <memory>
 
+#include <dune/common/rangeutilities.hh>
+
 #include <dune/istl/bcrsmatrix.hh>
 
 /** \file
@@ -65,6 +67,33 @@ namespace Dune {
 
     }
 
+    /** \brief Construct from a std::initializer_list */
+    BDMatrix (std::initializer_list<B> const &list)
+    : BDMatrix(list.size())
+    {
+      size_t i=0;
+      for (auto it = list.begin(); it != list.end(); ++it, ++i)
+        (*this)[i][i] = *it;
+    }
+
+    /** \brief Resize the matrix.  Invalidates the content! */
+    void setSize(size_type size)
+    {
+      this->BCRSMatrix<B,A>::setSize(size,   // rows
+                                     size,   // columns
+                                     size);  // nonzeros
+
+      for (auto i : range(size))
+        this->BCRSMatrix<B,A>::setrowsize(i, 1);
+
+      this->BCRSMatrix<B,A>::endrowsizes();
+
+      for (auto i : range(size))
+        this->BCRSMatrix<B,A>::addindex(i, i);
+
+      this->BCRSMatrix<B,A>::endindices();
+    }
+
     //! assignment
     BDMatrix& operator= (const BDMatrix& other) {
       this->BCRSMatrix<B,A>::operator=(other);

dune/istl/btdmatrix.hh

@@ -49,46 +49,53 @@ namespace Dune {
     /** \brief Default constructor */
     BTDMatrix() : BCRSMatrix<B,A>() {}
 
-    explicit BTDMatrix(int size)
+    explicit BTDMatrix(size_type size)
       : BCRSMatrix<B,A>(size, size, BCRSMatrix<B,A>::random)
     {
-      // special handling for 1x1 matrices
-      if (size==1) {
-
-        this->BCRSMatrix<B,A>::setrowsize(0, 1);
-        this->BCRSMatrix<B,A>::endrowsizes();
+      // Set number of entries for each row
+      // All rows get three entries, except for the first and the last one
+      for (size_t i=0; i<size; i++)
+        this->BCRSMatrix<B,A>::setrowsize(i, 3 - (i==0) - (i==(size-1)));
 
-        this->BCRSMatrix<B,A>::addindex(0, 0);
-        this->BCRSMatrix<B,A>::endindices();
+      this->BCRSMatrix<B,A>::endrowsizes();
 
-        return;
+      // The actual entries for each row
+      for (size_t i=0; i<size; i++) {
+        if (i>0)
+          this->BCRSMatrix<B,A>::addindex(i, i-1);
+        this->BCRSMatrix<B,A>::addindex(i, i  );
+        if (i<size-1)
+          this->BCRSMatrix<B,A>::addindex(i, i+1);
       }
 
-      // Set number of entries for each row
-      this->BCRSMatrix<B,A>::setrowsize(0, 2);
+      this->BCRSMatrix<B,A>::endindices();
+    }
 
-      for (int i=1; i<size-1; i++)
-        this->BCRSMatrix<B,A>::setrowsize(i, 3);
+    /** \brief Resize the matrix.  Invalidates the content! */
+    void setSize(size_type size)
+    {
+      auto nonZeros = (size==0) ? 0 : size + 2*(size-1);
+      this->BCRSMatrix<B,A>::setSize(size,   // rows
+                                     size,   // columns
+                                     nonZeros);
 
-      this->BCRSMatrix<B,A>::setrowsize(size-1, 2);
+      // Set number of entries for each row
+      // All rows get three entries, except for the first and the last one
+      for (size_t i=0; i<size; i++)
+        this->BCRSMatrix<B,A>::setrowsize(i, 3 - (i==0) - (i==(size-1)));
 
       this->BCRSMatrix<B,A>::endrowsizes();
 
       // The actual entries for each row
-      this->BCRSMatrix<B,A>::addindex(0, 0);
-      this->BCRSMatrix<B,A>::addindex(0, 1);
-
-      for (int i=1; i<size-1; i++) {
-        this->BCRSMatrix<B,A>::addindex(i, i-1);
+      for (size_t i=0; i<size; i++) {
+        if (i>0)
+          this->BCRSMatrix<B,A>::addindex(i, i-1);
         this->BCRSMatrix<B,A>::addindex(i, i  );
-        this->BCRSMatrix<B,A>::addindex(i, i+1);
+        if (i<size-1)
+          this->BCRSMatrix<B,A>::addindex(i, i+1);
       }
 
-      this->BCRSMatrix<B,A>::addindex(size-1, size-2);
-      this->BCRSMatrix<B,A>::addindex(size-1, size-1);
-
       this->BCRSMatrix<B,A>::endindices();
-
     }
 
     //! assignment

dune/istl/bvector.hh

@@ -4,17 +4,24 @@
 #ifndef DUNE_ISTL_BVECTOR_HH
 #define DUNE_ISTL_BVECTOR_HH
 
+#include <algorithm>
 #include <cmath>
 #include <complex>
-#include <memory>
+#include <initializer_list>
 #include <limits>
+#include <memory>
+#include <utility>
+#include <vector>
 
-#include <dune/common/promotiontraits.hh>
+#include <dune/common/deprecated.hh>
 #include <dune/common/dotproduct.hh>
 #include <dune/common/ftraits.hh>
+#include <dune/common/promotiontraits.hh>
+#include <dune/common/typetraits.hh>
+#include <dune/common/unused.hh>
 
-#include "istlexception.hh"
 #include "basearray.hh"
+#include "istlexception.hh"
 
 /*! \file
 
@@ -25,9 +32,6 @@
 
 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 {
 
@@ -195,6 +199,7 @@ namespace Imp {
     //! two norm sqrt(sum over squared values of entries)
     typename FieldTraits<field_type>::real_type two_norm () const
     {
+      using std::sqrt;
       typename FieldTraits<field_type>::real_type sum=0;
       for (size_type i=0; i<this->n; ++i) sum += (*this)[i].two_norm2();
       return sqrt(sum);
@@ -210,7 +215,7 @@ namespace Imp {
 
     //! infinity norm (maximum of absolute values of entries)
     template <typename ft = field_type,
-              typename std::enable_if<!has_nan<ft>::value, int>::type = 0>
+              typename std::enable_if<!HasNaN<ft>::value, int>::type = 0>
     typename FieldTraits<ft>::real_type infinity_norm() const {
       using real_type = typename FieldTraits<ft>::real_type;
       using std::max;
@@ -225,7 +230,7 @@ namespace Imp {
 
     //! simplified infinity norm (uses Manhattan norm for complex values)
     template <typename ft = field_type,
-              typename std::enable_if<!has_nan<ft>::value, int>::type = 0>
+              typename std::enable_if<!HasNaN<ft>::value, int>::type = 0>
     typename FieldTraits<ft>::real_type infinity_norm_real() const {
       using real_type = typename FieldTraits<ft>::real_type;
       using std::max;
@@ -240,7 +245,7 @@ namespace Imp {
 
     //! infinity norm (maximum of absolute values of entries)
     template <typename ft = field_type,
-              typename std::enable_if<has_nan<ft>::value, int>::type = 0>
+              typename std::enable_if<HasNaN<ft>::value, int>::type = 0>
     typename FieldTraits<ft>::real_type infinity_norm() const {
       using real_type = typename FieldTraits<ft>::real_type;
       using std::max;
@@ -258,7 +263,7 @@ namespace Imp {
 
     //! simplified infinity norm (uses Manhattan norm for complex values)
     template <typename ft = field_type,
-              typename std::enable_if<has_nan<ft>::value, int>::type = 0>
+              typename std::enable_if<HasNaN<ft>::value, int>::type = 0>
     typename FieldTraits<ft>::real_type infinity_norm_real() const {
       using real_type = typename FieldTraits<ft>::real_type;
       using std::max;
@@ -297,6 +302,39 @@ namespace Imp {
     {       }
   };
 
+  //! simple scope guard, execute the provided functor on scope exit
+  /**
+   * The guard may not be copied or moved.  This avoids executing the cleanup
+   * function twice.  The cleanup function should not throw, as it may be
+   * called during stack unwinding.
+   */
+  template<class F>
+  class ScopeGuard {
+    F cleanupFunc_;
+  public:
+    ScopeGuard(F cleanupFunc) : cleanupFunc_(std::move(cleanupFunc)) {}
+    ScopeGuard(const ScopeGuard &) = delete;
+    ScopeGuard(ScopeGuard &&) = delete;
+    ScopeGuard &operator=(ScopeGuard) = delete;
+    ~ScopeGuard() { cleanupFunc_(); }
+  };
+
+  //! create a scope guard
+  /**
+   * Use like
+   * ```c++
+   * {
+   *   const auto &guard = makeScopeGuard([&]{ doSomething(); });
+   *   doSomethingThatMightThrow();
+   * }
+   * ```
+   */
+  template<class F>
+  ScopeGuard<F> makeScopeGuard(F cleanupFunc)
+  {
+    return { std::move(cleanupFunc) };
+  }
+
 } // end namespace Imp
   /**
      @addtogroup ISTL_SPMV
@@ -346,44 +384,21 @@ namespace Imp {
     //===== constructors and such
 
     //! makes empty vector
-    BlockVector () : Imp::block_vector_unmanaged<B,A>(),
-                     capacity_(0)
-    {}
+    BlockVector ()
+    {
+      syncBaseArray();
+    }
 
     //! make vector with _n components
-    explicit BlockVector (size_type _n)
+    explicit BlockVector (size_type _n) : storage_(_n)
     {
-      this->n = _n;
-      capacity_ = _n;
-      if (capacity_>0) {
-        this->p = this->allocator_.allocate(capacity_);
-        // actually construct the objects
-        new(this->p)B[capacity_];
-      } else
-      {
-        this->p = 0;
-        this->n = 0;
-        capacity_ = 0;
-      }
+      syncBaseArray();
     }
 
     /** \brief Construct from a std::initializer_list */
-    BlockVector (std::initializer_list<B> const &l)
+    BlockVector (std::initializer_list<B> const &l) : storage_(l)
     {
-      this->n = l.size();
-      capacity_ = l.size();
-      if (capacity_>0) {
-        this->p = this->allocator_.allocate(capacity_);
-        // actually construct the objects
-        new(this->p)B[capacity_];
-
-        std::copy_n(l.begin(), l.size(), this->p);
-      } else
-      {
-        this->p = 0;
-        this->n = 0;
-        capacity_ = 0;
-      }
+      syncBaseArray();
     }
 
 
@@ -403,25 +418,29 @@ namespace Imp {
     {
       static_assert(std::numeric_limits<S>::is_integer,
         "capacity must be an unsigned integral type (be aware, that this constructor does not set the default value!)" );
-      size_type capacity = _capacity;
-      this->n = _n;
-      if(this->n > capacity)
-        capacity_ = _n;
-      else
-        capacity_ = capacity;
-
-      if (capacity_>0) {
-        this->p = this->allocator_.allocate(capacity_);
-        new (this->p)B[capacity_];
-      } else
-      {
-        this->p = 0;
-        this->n = 0;
-        capacity_ = 0;
-      }
+      if(_capacity > _n)
+        storage_.reserve(_capacity);
+      storage_.resize(_n);
+      syncBaseArray();
     }
 
 
+    /**
+     * @brief Reserve space.
+     *
+     * Allocate storage for up to `capacity` blocks.  Resizing won't cause
+     * reallocation until the size exceeds the `capacity`
+     *
+     * @param capacity The maximum number of elements the vector
+     *        needs to hold.
+     */
+    void reserve(size_type capacity)
+    {
+      DUNE_UNUSED const auto &guard =
+        Imp::makeScopeGuard([this]{ syncBaseArray(); });
+      storage_.reserve(capacity);
+    }
+
     /**
      * @brief Reserve space.
      *
@@ -435,44 +454,17 @@ namespace Imp {
      *
      * @param capacity The maximum number of elements the vector
      * needs to hold.
-     * @param copyOldValues If false no object will be copied and the data might be
-     * lost. Default value is true.
+     * @param copyOldValues Ignored, values are always copied.
+     *
+     * \deprecated This method is deprecated, since the extra copyOldValues
+     *             parameter is unusual, and the previous implementation did
+     *             not always reduce memory anyway.
      */
-    void reserve(size_type capacity, bool copyOldValues=true)
-    {
-      if(capacity >= Imp::block_vector_unmanaged<B,A>::N() && capacity != capacity_) {
-        // save the old data
-        B* pold = this->p;
-
-        if(capacity>0) {
-          // create new array with capacity
-          this->p = this->allocator_.allocate(capacity);
-          new (this->p)B[capacity];
-
-          if(copyOldValues) {
-            // copy the old values
-            B* to = this->p;
-            B* from = pold;
-
-            for(size_type i=0; i < Imp::block_vector_unmanaged<B,A>::N(); ++i, ++from, ++to)
-              *to = *from;
-          }
-          if(capacity_ > 0) {
-            // Destruct old objects and free memory
-            int i=capacity_;
-            while (i)
-              pold[--i].~B();
-            this->allocator_.deallocate(pold,capacity_);
-          }
-        }else{
-          if(capacity_ > 0)
-            // free old data
-            this->p = 0;
-          capacity_ = 0;
-        }
-
-        capacity_ = capacity;
-      }
+    DUNE_DEPRECATED_MSG("Use the overload without the second parameter, "
+                        "values are always copied")
+    void reserve(size_type capacity, bool copyOldValues)
+    {
+      reserve(capacity);
     }
 
     /**
@@ -483,7 +475,24 @@ namespace Imp {
      */
     size_type capacity() const
     {
-      return capacity_;
+      return storage_.capacity();
+    }
+
+    /**
+     * @brief Resize the vector.
+     *
+     * Resize the vector to the given number of blocks.  Blocks below the
+     * given size are copied (moved if possible).  Old blocks above the given
+     * size are destructed, new blocks above the given size are
+     * value-initialized.
+     *
+     * @param size The new number of blocks of the vector.
+     */
+    void resize(size_type size)
+    {
+      DUNE_UNUSED const auto &guard =
+        Imp::makeScopeGuard([this]{ syncBaseArray(); });
+      storage_.resize(size);
     }
 
     /**
@@ -497,84 +506,67 @@ namespace Imp {
      * will be copied if the capacity changes.  If it is false
      * the old values are lost.
      * @param size The new size of the vector.
-     * @param copyOldValues If false no object will be copied and the data might be
-     * lost.
+     * @param copyOldValues Ignored, values are always copied.
+     *
+     * \deprecated This method is deprecated, since the extra copyOldValues
+     *             parameter is unusual and conflicts with the usual meaning
+     *             (default value for newly created elements).
      */
-    void resize(size_type size, bool copyOldValues=true)
+    DUNE_DEPRECATED_MSG("Use the overload without the second parameter, "
+                        "values are always copied")
+    void resize(size_type size, bool copyOldValues)
     {
-      if (size > Imp::block_vector_unmanaged<B,A>::N())
-        if(capacity_ < size)
-          this->reserve(size, copyOldValues);
-      this->n = size;
+      resize(size);
     }
 
 
 
 
     //! copy constructor
-    BlockVector (const BlockVector& a) :
-      Imp::block_vector_unmanaged<B,A>(a)
+    BlockVector(const BlockVector &a)
+      noexcept(noexcept(std::declval<BlockVector>().storage_ = a.storage_))
     {
-      // allocate memory with same size as a
-      this->n = a.n;
-      capacity_ = a.capacity_;
-
-      if (capacity_>0) {
-        this->p = this->allocator_.allocate(capacity_);
-        new (this->p)B[capacity_];
-      } else
-      {
-        this->n = 0;
-        this->p = 0;
-      }
-
-      // and copy elements
-      for (size_type i=0; i<this->n; i++) this->p[i]=a.p[i];
+      storage_ = a.storage_;
+      syncBaseArray();
     }
 
-    //! free dynamic memory
-    ~BlockVector ()
+    //! move constructor
+    BlockVector(BlockVector &&a)
+      noexcept(noexcept(std::declval<BlockVector>().swap(a)))
     {
-      if (capacity_>0) {
-        int i=capacity_;
-        while (i)
-          this->p[--i].~B();
-        this->allocator_.deallocate(this->p,capacity_);
-      }
+      swap(a);
     }
 
     //! assignment
     BlockVector& operator= (const BlockVector& a)
+      noexcept(noexcept(std::declval<BlockVector>().storage_ = a.storage_))
     {
-      if (&a!=this)     // check if this and a are different objects
-      {
-        // adjust size of vector
-        if (capacity_!=a.capacity_)           // check if size is different
-        {
-          if (capacity_>0) {
-            int i=capacity_;
-            while (i)
-              this->p[--i].~B();
-            this->allocator_.deallocate(this->p,capacity_);                     // free old memory
-          }
-          capacity_ = a.capacity_;
-          if (capacity_>0) {
-            this->p = this->allocator_.allocate(capacity_);
-            new (this->p)B[capacity_];
-          } else
-          {
-            this->p = 0;
-            capacity_ = 0;
-          }
-        }
-        this->n = a.n;
-        // copy data
-        for (size_type i=0; i<this->n; i++)
-          this->p[i]=a.p[i];
-      }
+      DUNE_UNUSED const auto &guard =
+        Imp::makeScopeGuard([this]{ syncBaseArray(); });
+      storage_ = a.storage_;
       return *this;
     }
 
+    //! move assignment
+    BlockVector& operator= (BlockVector&& a)
+      noexcept(noexcept(std::declval<BlockVector>().swap(a)))
+    {
+      swap(a);
+      return *this;
+    }
+
+    //! swap operation
+    void swap(BlockVector &other)
+      noexcept(noexcept(
+            std::declval<BlockVector&>().storage_.swap(other.storage_)))
+    {
+      DUNE_UNUSED const auto &guard = Imp::makeScopeGuard([&]{
+          syncBaseArray();
+          other.syncBaseArray();
+        });
+      storage_.swap(other.storage_);
+    }
+
     //! assign from scalar
     BlockVector& operator= (const field_type& k)
     {
@@ -583,21 +575,14 @@ namespace Imp {
       return *this;
     }
 
-    //! Assignment from BlockVectorWindow
-    template<class OtherAlloc>
-    BlockVector& operator= (const BlockVectorWindow<B,OtherAlloc>& other)
+  private:
+    void syncBaseArray() noexcept
     {
-      resize(other.size());
-      for(std::size_t i=0; i<other.size(); ++i)
-        (*this)[i] = other[i];
-      return *this;
+      this->p = storage_.data();
+      this->n = storage_.size();
     }
 
-  protected:
-    size_type capacity_;
-
-    A allocator_;
-
+    std::vector<B, A> storage_;
   };
 
   /** @} */
@@ -726,6 +711,14 @@ namespace Imp {
       return *this;
     }
 
+    //! copy into an independent BlockVector object
+    operator BlockVector<B, A>() const {
+      auto bv = BlockVector<B, A>(this->n);
+
+      std::copy(this->begin(), this->end(), bv.begin());
+
+      return bv;
+    }
 
     //===== window manipulation methods
 
@@ -755,7 +748,7 @@ namespace Imp {
     }
 
     //! get size
-    size_type getsize ()
+    size_type getsize () const
     {
       return this->n;
     }
@@ -895,6 +888,7 @@ namespace Imp {
     //! two norm sqrt(sum over squared values of entries)
     typename FieldTraits<field_type>::real_type two_norm () const
     {
+      using std::sqrt;
       typename FieldTraits<field_type>::real_type sum=0;
       for (size_type i=0; i<this->n; ++i) sum += (this->p)[i].two_norm2();
       return sqrt(sum);
@@ -910,7 +904,7 @@ namespace Imp {
 
     //! infinity norm (maximum of absolute values of entries)
     template <typename ft = field_type,
-              typename std::enable_if<!has_nan<ft>::value, int>::type = 0>
+              typename std::enable_if<!HasNaN<ft>::value, int>::type = 0>
     typename FieldTraits<ft>::real_type infinity_norm() const {
       using real_type = typename FieldTraits<ft>::real_type;
       using std::max;
@@ -925,7 +919,7 @@ namespace Imp {
 
     //! simplified infinity norm (uses Manhattan norm for complex values)
     template <typename ft = field_type,
-              typename std::enable_if<!has_nan<ft>::value, int>::type = 0>
+              typename std::enable_if<!HasNaN<ft>::value, int>::type = 0>
     typename FieldTraits<ft>::real_type infinity_norm_real() const {
       using real_type = typename FieldTraits<ft>::real_type;
       using std::max;
@@ -940,7 +934,7 @@ namespace Imp {
 
     //! infinity norm (maximum of absolute values of entries)
     template <typename ft = field_type,
-              typename std::enable_if<has_nan<ft>::value, int>::type = 0>
+              typename std::enable_if<HasNaN<ft>::value, int>::type = 0>
     typename FieldTraits<ft>::real_type infinity_norm() const {
       using real_type = typename FieldTraits<ft>::real_type;
       using std::max;
@@ -958,7 +952,7 @@ namespace Imp {
 
     //! simplified infinity norm (uses Manhattan norm for complex values)
     template <typename ft = field_type,
-              typename std::enable_if<has_nan<ft>::value, int>::type = 0>
+              typename std::enable_if<HasNaN<ft>::value, int>::type = 0>
     typename FieldTraits<ft>::real_type infinity_norm_real() const {
       using real_type = typename FieldTraits<ft>::real_type;
       using std::max;
@@ -1165,6 +1159,15 @@ namespace Imp {
 
 } // end namespace 'Imp'
 
+
+  //! Specialization for the proxies of `BlockVectorWindow`
+  template<typename B, typename A>
+  struct AutonomousValueType<Imp::BlockVectorWindow<B,A>>
+  {
+    using type = BlockVector<B, A>;
+  };
+
+
 } // end namespace 'Dune'
 
 #endif

dune/istl/eigenvalue/arpackpp.hh

@@ -69,8 +69,8 @@ namespace Dune
 
         // allocate memory for auxiliary block vector objects
         // which are compatible to matrix rows / columns
-        domainBlockVector.resize(A_.N(),false);
-        rangeBlockVector.resize(A_.M(),false);
+        domainBlockVector.resize(A_.N());
+        rangeBlockVector.resize(A_.M());
       }
 
       //! Perform matrix-vector product w = A*v
@@ -308,7 +308,7 @@ namespace Dune
 
       // allocate memory for variables, set parameters
       const int nev = 1;                     // Number of eigenvalues to compute
-      const int ncv = 20;                    // Number of Arnoldi vectors generated at each iteration (0 == auto)
+      int ncv = std::min(20, nrows);         // Number of Arnoldi vectors generated at each iteration (0 == auto)
       const Real tol = epsilon;              // Stopping tolerance (relative accuracy of Ritz values) (0 == machine precision)
       const int maxit = nIterationsMax_*nev; // Maximum number of Arnoldi update iterations allowed   (0 == 100*nev)
       Real* ev = new Real[nev];              // Computed eigenvalues of A
@@ -410,7 +410,7 @@ namespace Dune
 
       // allocate memory for variables, set parameters
       const int nev = 1;                     // Number of eigenvalues to compute
-      const int ncv = 20;                    // Number of Arnoldi vectors generated at each iteration (0 == auto)
+      int ncv = std::min(20, nrows);         // Number of Arnoldi vectors generated at each iteration (0 == auto)
       const Real tol = epsilon;              // Stopping tolerance (relative accuracy of Ritz values) (0 == machine precision)
       const int maxit = nIterationsMax_*nev; // Maximum number of Arnoldi update iterations allowed   (0 == 100*nev)
       Real* ev = new Real[nev];              // Computed eigenvalues of A
@@ -511,7 +511,7 @@ namespace Dune
 
       // allocate memory for variables, set parameters
       const int nev = 2;                     // Number of eigenvalues to compute
-      const int ncv = 20;                    // Number of Arnoldi vectors generated at each iteration (0 == auto)
+      int ncv = std::min(20, nrows);         // Number of Arnoldi vectors generated at each iteration (0 == auto)
       const Real tol = epsilon;              // Stopping tolerance (relative accuracy of Ritz values) (0 == machine precision)
       const int maxit = nIterationsMax_*nev; // Maximum number of Arnoldi update iterations allowed   (0 == 100*nev)
       Real* ev = new Real[nev];              // Computed eigenvalues of A
@@ -629,7 +629,7 @@ namespace Dune
 
       // allocate memory for variables, set parameters
       const int nev = 1;                     // Number of eigenvalues to compute
-      const int ncv = 20;                    // Number of Arnoldi vectors generated at each iteration (0 == auto)
+      int ncv = std::min(20, nrows);         // Number of Arnoldi vectors generated at each iteration (0 == auto)
       const Real tol = epsilon;              // Stopping tolerance (relative accuracy of Ritz values) (0 == machine precision)
       const int maxit = nIterationsMax_*nev; // Maximum number of Arnoldi update iterations allowed   (0 == 100*nev)
       Real* ev = new Real[nev];              // Computed eigenvalues of A^T*A
@@ -741,7 +741,7 @@ namespace Dune
 
       // allocate memory for variables, set parameters
       const int nev = 1;                     // Number of eigenvalues to compute
-      const int ncv = 20;                    // Number of Arnoldi vectors generated at each iteration (0 == auto)
+      int ncv = std::min(20, nrows);         // Number of Arnoldi vectors generated at each iteration (0 == auto)
       const Real tol = epsilon;              // Stopping tolerance (relative accuracy of Ritz values) (0 == machine precision)
       const int maxit = nIterationsMax_*nev; // Maximum number of Arnoldi update iterations allowed   (0 == 100*nev)
       Real* ev = new Real[nev];              // Computed eigenvalues of A^T*A
@@ -849,7 +849,7 @@ namespace Dune
 
       // allocate memory for variables, set parameters
       const int nev = 2;                     // Number of eigenvalues to compute
-      const int ncv = 20;                    // Number of Arnoldi vectors generated at each iteration (0 == auto)
+      int ncv = std::min(20, nrows);         // Number of Arnoldi vectors generated at each iteration (0 == auto)
       const Real tol = epsilon;              // Stopping tolerance (relative accuracy of Ritz values) (0 == machine precision)
       const int maxit = nIterationsMax_*nev; // Maximum number of Arnoldi update iterations allowed   (0 == 100*nev)
       Real* ev = new Real[nev];              // Computed eigenvalues of A^T*A

dune/istl/eigenvalue/poweriteration.hh

@@ -47,8 +47,6 @@ namespace Dune
       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),
@@ -68,6 +66,12 @@ namespace Dune
         y.axpy(alpha,temp);
       }
 
+      //! Category of the linear operator (see SolverCategory::Category)
+      virtual SolverCategory::Category category() const
+      {
+        return SolverCategory::sequential;
+      }
+
     protected:
       const field_type immutable_scaling_;
       const field_type& mutable_scaling_;
@@ -92,8 +96,6 @@ namespace Dune
       typedef typename OP1::range_type range_type;
       typedef typename domain_type::field_type field_type;
 
-      enum {category = Dune::SolverCategory::sequential};
-
       LinearOperatorSum (const OP1& op1, const OP2& op2)
         : op1_(op1), op2_(op2)
       {
@@ -118,6 +120,12 @@ namespace Dune
         y.axpy(alpha,temp);
       }
 
+      //! Category of the linear operator (see SolverCategory::Category)
+      virtual SolverCategory::Category category() const
+      {
+        return SolverCategory::sequential;
+      }
+
     protected:
       const OP1& op1_;
       const OP2& op2_;
@@ -935,7 +943,7 @@ namespace Dune
     {
       // create iteration matrix on demand
       if (!itMatrix_)
-        itMatrix_ = std::unique_ptr<BCRSMatrix>(new BCRSMatrix(m_));
+        itMatrix_ = std::make_unique<BCRSMatrix>(m_);
 
       // return iteration matrix
       return *itMatrix_;

dune/istl/eigenvalue/test/matrixinfo.hh

@@ -5,7 +5,7 @@
 
 #include <cmath>    // provides std::abs and std::sqrt
 #include <cassert>  // provides assert
-
+#include <limits>
 #include <iostream>  // provides std::cout, std::endl
 
 #include <dune/common/exceptions.hh>  // provides DUNE_THROW(...), Dune::Exception
@@ -191,7 +191,7 @@ protected:
     // 7) get smallest magnitude eigenvalue via TLIME iteration
     //    (assume that m_ is symmetric)
     {
-      const Real epsilon = 1e-11;
+      const Real epsilon = std::sqrt(std::numeric_limits<Real>::epsilon());
       x = 1.0;
       // 7.1) perform TLIME iteration for smallest magnitude
       //      eigenvalue
@@ -322,7 +322,7 @@ protected:
     // 9) get smallest singular value as square root of the smallest
     //    magnitude eigenvalue of m^t*m via TLIME iteration
     {
-      const Real epsilon = 1e-11;
+      const Real epsilon = std::sqrt(std::numeric_limits<Real>::epsilon());
       x = 1.0;
       // 9.1) perform TLIME iteration for smallest magnitude
       //      eigenvalue of m^t*m