diff --git a/.travis.yml b/.travis.yml
index e36de0abe3..0eacd38bde 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -41,25 +41,25 @@ jobs:
   include:
     - stage: Staticcheck
       python: 3.6.7
-      env: xTHREADING=0 xMPI=0 xSIONLIB=0 xGSL=0 xLIBNEUROSIM=0 xLTDL=0 xREADLINE=0 xPYTHON=0 xMUSIC=0 xSTATIC_ANALYSIS=1 xRUN_BUILD_AND_TESTSUITE=0 CACHE_NAME=JOB   # only static code analysis
+      env: xTHREADING=0 xMPI=0 xSIONLIB=0 xGSL=0 xLIBNEUROSIM=0 xLTDL=0 xREADLINE=0 xLIBBOOST=0 xPYTHON=0 xMUSIC=0 xSTATIC_ANALYSIS=1 xRUN_BUILD_AND_TESTSUITE=0 CACHE_NAME=JOB   # only static code analysis
     - stage: MPI-Threading-Python
       python: 3.6.7
-      env: xTHREADING=0 xMPI=1 xSIONLIB=0 xGSL=0 xLIBNEUROSIM=0 xLTDL=1 xREADLINE=1 xPYTHON=0 xMUSIC=0 xSTATIC_ANALYSIS=0 xRUN_BUILD_AND_TESTSUITE=1 CACHE_NAME=JOB   # only MPI
+      env: xTHREADING=0 xMPI=1 xSIONLIB=0 xGSL=0 xLIBNEUROSIM=0 xLTDL=1 xREADLINE=1 xLIBBOOST=1 xPYTHON=0 xMUSIC=0 xSTATIC_ANALYSIS=0 xRUN_BUILD_AND_TESTSUITE=1 CACHE_NAME=JOB   # only MPI
     - stage: MPI-Threading-Python
       python: 3.6.7
-      env: xTHREADING=1 xMPI=0 xSIONLIB=0 xGSL=0 xLIBNEUROSIM=0 xLTDL=1 xREADLINE=1 xPYTHON=0 xMUSIC=0 xSTATIC_ANALYSIS=0 xRUN_BUILD_AND_TESTSUITE=1 CACHE_NAME=JOB   # only threading
+      env: xTHREADING=1 xMPI=0 xSIONLIB=0 xGSL=0 xLIBNEUROSIM=0 xLTDL=1 xREADLINE=1 xLIBBOOST=1 xPYTHON=0 xMUSIC=0 xSTATIC_ANALYSIS=0 xRUN_BUILD_AND_TESTSUITE=1 CACHE_NAME=JOB   # only threading
     - stage: MPI-Threading-Python
       python: 3.6.7
-      env: xTHREADING=1 xMPI=0 xSIONLIB=0 xGSL=0 xLIBNEUROSIM=0 xLTDL=0 xREADLINE=0 xPYTHON=1 xMUSIC=0 xSTATIC_ANALYSIS=0 xRUN_BUILD_AND_TESTSUITE=1 CACHE_NAME=JOB   # Python & Threading
+      env: xTHREADING=1 xMPI=0 xSIONLIB=0 xGSL=0 xLIBNEUROSIM=0 xLTDL=0 xREADLINE=0 xLIBBOOST=1 xPYTHON=1 xMUSIC=0 xSTATIC_ANALYSIS=0 xRUN_BUILD_AND_TESTSUITE=1 CACHE_NAME=JOB   # Python & Threading
     - stage: MPI-Threading-Python
       python: 3.6.7
-      env: xTHREADING=0 xMPI=1 xSIONLIB=0 xGSL=0 xLIBNEUROSIM=0 xLTDL=0 xREADLINE=0 xPYTHON=1 xMUSIC=0 xSTATIC_ANALYSIS=0 xRUN_BUILD_AND_TESTSUITE=1 CACHE_NAME=JOB   # Python & MPI
+      env: xTHREADING=0 xMPI=1 xSIONLIB=0 xGSL=0 xLIBNEUROSIM=0 xLTDL=0 xREADLINE=0 xLIBBOOST=1 xPYTHON=1 xMUSIC=0 xSTATIC_ANALYSIS=0 xRUN_BUILD_AND_TESTSUITE=1 CACHE_NAME=JOB   # Python & MPI
     - stage: MPI-Threading-Python
       python: 3.6.7
-      env: xTHREADING=0 xMPI=0 xSIONLIB=0 xGSL=0 xLIBNEUROSIM=0 xLTDL=0 xREADLINE=0 xPYTHON=1 xMUSIC=0 xSTATIC_ANALYSIS=0 xRUN_BUILD_AND_TESTSUITE=1 CACHE_NAME=JOB   # only Python
+      env: xTHREADING=0 xMPI=0 xSIONLIB=0 xGSL=0 xLIBNEUROSIM=0 xLTDL=0 xREADLINE=0 xLIBBOOST=1 xPYTHON=1 xMUSIC=0 xSTATIC_ANALYSIS=0 xRUN_BUILD_AND_TESTSUITE=1 CACHE_NAME=JOB   # only Python
     - stage: Python-Full-build
       python: 3.6.7
-      env: xTHREADING=1 xMPI=1 xSIONLIB=1 xGSL=1 xLIBNEUROSIM=1 xLTDL=1 xREADLINE=1 xPYTHON=1 xMUSIC=1 xSTATIC_ANALYSIS=0 xRUN_BUILD_AND_TESTSUITE=1 CACHE_NAME=JOB   # full
+      env: xTHREADING=1 xMPI=1 xSIONLIB=1 xGSL=1 xLIBNEUROSIM=1 xLTDL=1 xREADLINE=1 xLIBBOOST=1 xPYTHON=1 xMUSIC=1 xSTATIC_ANALYSIS=0 xRUN_BUILD_AND_TESTSUITE=1 CACHE_NAME=JOB   # full
     - stage: Clang7
       language: cpp
       env: MATRIX_EVAL="CC=clang-7 && CXX=clang++-7" xRUN_BUILD_AND_TESTSUITE=1
@@ -139,6 +139,7 @@ before_install:
          cp extras/install_csa-libneurosim.sh $HOME
          cp extras/install_music.sh $HOME
          cp extras/install_sionlib.sh $HOME
+         cp extras/install_libboost.sh $HOME
          cd $HOME/build
          echo $PATH
          # Upgrade pip and setuptools
diff --git a/cmake/ProcessOptions.cmake b/cmake/ProcessOptions.cmake
index 0655a9d708..92da8c7b68 100644
--- a/cmake/ProcessOptions.cmake
+++ b/cmake/ProcessOptions.cmake
@@ -458,7 +458,7 @@ function( NEST_PROCESS_WITH_OPENMP )
   if ( NOT TARGET OpenMP::OpenMP_CXX )
     add_library(OpenMP::OpenMP_CXX INTERFACE IMPORTED)
   endif()
- 
+
 endfunction()
 
 function( NEST_PROCESS_WITH_MPI )
@@ -580,8 +580,10 @@ function( NEST_PROCESS_WITH_BOOST )
       set( BOOST_ROOT "${with-boost}" )
     endif ()
 
+    set(Boost_USE_DEBUG_LIBS OFF)  # ignore debug libs
+    set(Boost_USE_RELEASE_LIBS ON) # only find release libs
     # Needs Boost version >=1.58.0 to use Boost sorting
-    find_package( Boost 1.58.0 COMPONENTS unit_test_framework )
+    find_package( Boost 1.58.0 )
     if ( Boost_FOUND )
       # export found variables to parent scope
       set( HAVE_BOOST ON PARENT_SCOPE )
diff --git a/extras/install_libboost.sh b/extras/install_libboost.sh
new file mode 100644
index 0000000000..4c606ba483
--- /dev/null
+++ b/extras/install_libboost.sh
@@ -0,0 +1,9 @@
+#!/bin/bash
+
+# Install the boost library
+wget --no-verbose https://dl.bintray.com/boostorg/release/1.72.0/source/boost_1_72_0.tar.gz
+tar -xzf boost_1_72_0.tar.gz
+rm -fr boost_1_72_0.tar.gz
+cp -fr boost_1_72_0 $HOME/.cache/boost_1_72_0.install
+rm -fr boost_1_72_0
+
diff --git a/extras/travis_build.sh b/extras/travis_build.sh
index 4e1e479f88..d51326a7de 100755
--- a/extras/travis_build.sh
+++ b/extras/travis_build.sh
@@ -32,11 +32,6 @@
 # Exit shell if any subcommand or pipline returns a non-zero status.
 set -e
 
-mkdir -p $HOME/.matplotlib
-cat > $HOME/.matplotlib/matplotlibrc <<EOF
-    backend : svg
-EOF
-
 # Set the NEST CMake-build configuration according to the build matrix in '.travis.yml'.
 if [ "$xTHREADING" = "1" ] ; then
     CONFIGURE_THREADING="-Dwith-openmp=ON"
@@ -57,6 +52,10 @@ if [ "$xPYTHON" = "1" ] ; then
    if [[ $OSTYPE = darwin* ]]; then
       CONFIGURE_PYTHON="-DPYTHON_LIBRARY=/usr/local/Cellar/python/3.7.5/Frameworks/Python.framework/Versions/3.7/lib/libpython3.7.dylib -DPYTHON_INCLUDE_DIR=/usr/local/Cellar/python/3.7.5/Frameworks/Python.framework/Versions/3.7/include//python3.7m/"
    fi
+   mkdir -p $HOME/.matplotlib
+   cat > $HOME/.matplotlib/matplotlibrc <<EOF 
+   backend : svg
+EOF
 else
     CONFIGURE_PYTHON="-Dwith-python=OFF"
 fi
@@ -87,6 +86,14 @@ else
     CONFIGURE_READLINE="-Dwith-readline=OFF"
 fi
 
+if [ "$xLIBBOOST" = "1" ] ; then
+    CONFIGURE_BOOST="-Dwith-boost=$HOME/.cache/boost_1_72_0.install"
+    chmod +x extras/install_libboost.sh
+    ./extras/install_libboost.sh
+else
+    CONFIGURE_BOOST="-Dwith-boost=OFF"
+fi
+
 if [ "$xSIONLIB" = "1" ] ; then
     CONFIGURE_SIONLIB="-Dwith-sionlib=$HOME/.cache/sionlib.install"
     chmod +x extras/install_sionlib.sh
@@ -106,9 +113,6 @@ fi
 if [[ $OSTYPE = darwin* ]]; then
     export CC=$(ls /usr/local/bin/gcc-* | grep '^/usr/local/bin/gcc-\d$')
     export CXX=$(ls /usr/local/bin/g++-* | grep '^/usr/local/bin/g++-\d$')
-    CONFIGURE_BOOST="-Dwith-boost=OFF"
-else
-    CONFIGURE_BOOST="-Dwith-boost=ON"
 fi
  
 NEST_VPATH=build
@@ -193,7 +197,7 @@ if [ "$xSTATIC_ANALYSIS" = "1" ]; then
     # Set the command line arguments for the static code analysis script and execute it.
 
     # The names of the static code analysis tools executables.
-    VERA=vera++
+    VERA=vera++                   
     CPPCHECK=cppcheck
     CLANG_FORMAT=clang-format
     PEP8=pep8
diff --git a/libnestutil/CMakeLists.txt b/libnestutil/CMakeLists.txt
index f4a379cdea..f8dc3bf2c6 100644
--- a/libnestutil/CMakeLists.txt
+++ b/libnestutil/CMakeLists.txt
@@ -40,6 +40,7 @@ target_link_libraries( nestutil ${GSL_LIBRARIES} ${SIONLIB_LIBS} )
 
 target_include_directories( nestutil PRIVATE
     ${PROJECT_BINARY_DIR}/libnestutil
+    ${Boost_INCLUDE_DIRS}
     )
 
 install( TARGETS nestutil
diff --git a/libnestutil/block_vector.h b/libnestutil/block_vector.h
index e1517b9a7d..847bcdeb43 100644
--- a/libnestutil/block_vector.h
+++ b/libnestutil/block_vector.h
@@ -29,6 +29,8 @@
 #include <iterator>
 #include <cassert>
 
+#include "sliexceptions.h"
+
 template < typename value_type_ >
 class BlockVector;
 template < typename value_type_, typename ref_, typename ptr_ >
@@ -75,7 +77,7 @@ class bv_iterator
   using value_type = value_type_;
   using pointer = ptr_;
   using reference = ref_;
-  using difference_type = long int;
+  using difference_type = typename BlockVector< value_type >::difference_type;
 
   bv_iterator() = default;
 
@@ -107,17 +109,23 @@ class bv_iterator
   bv_iterator& operator--();
   bv_iterator& operator+=( difference_type );
   bv_iterator& operator-=( difference_type );
-  bv_iterator operator+( difference_type );
-  bv_iterator operator-( difference_type );
+  bv_iterator operator+( difference_type ) const;
+  bv_iterator operator-( difference_type ) const;
   bv_iterator operator++( int );
+  bv_iterator operator--( int );
   reference operator*() const;
   pointer operator->() const;
   difference_type operator-( const iterator& ) const;
   difference_type operator-( const const_iterator& ) const;
 
+  reference operator[]( difference_type n ) const;
+
   bool operator==( const bv_iterator& ) const;
   bool operator!=( const bv_iterator& ) const;
   bool operator<( const bv_iterator& ) const;
+  bool operator>( const bv_iterator& ) const;
+  bool operator<=( const bv_iterator& ) const;
+  bool operator>=( const bv_iterator& ) const;
 
 private:
   /**
@@ -143,8 +151,15 @@ class BlockVector
   friend class bv_iterator;
 
 public:
+  using value_type = value_type_;
+  using difference_type = typename std::vector< value_type >::difference_type;
+  using const_reference = const value_type&;
+  using const_pointer = const value_type*;
   using iterator = bv_iterator< value_type_, value_type_&, value_type_* >;
   using const_iterator = bv_iterator< value_type_, const value_type_&, const value_type_* >;
+  using reverse_iterator = std::reverse_iterator< iterator >;
+  using const_reverse_iterator = std::reverse_iterator< const_iterator >;
+  using size_type = size_t;
 
   /**
    * @brief Creates an empty BlockVector.
@@ -247,6 +262,39 @@ class BlockVector
    */
   int get_max_block_size() const;
 
+  /**
+   * @brief Returns the size() of the largest possible BlockVector.
+   */
+  size_type max_size() const;
+
+  /**
+   * @brief Returns a read/write reverse iterator that points to the last
+   * element in the BlockVector. Iteration is done in reverse element
+   * order.
+   */
+  reverse_iterator rbegin();
+
+  /**
+   * @brief Returns a read-only (constant) reverse iterator that points to
+   * the last element in the BlockVector. Iteration is done in reverse
+   * element order.
+   */
+  reverse_iterator rbegin() const;
+
+  /**
+   * @brief Returns a read/write reverse iterator that points to one
+   * before the first element in the BlockVector. Iteration is done in
+   * reverse element order.
+   */
+  reverse_iterator rend();
+
+  /**
+   * @brief Returns a read-only (constant) reverse iterator that points to
+   * one before the first element in the BlockVector. Iteration is done in
+   * reverse element order.
+   */
+  reverse_iterator rend() const;
+
 private:
   //! Vector holding blocks containing data.
   std::vector< std::vector< value_type_ > > blockmap_;
@@ -463,6 +511,41 @@ BlockVector< value_type_ >::get_max_block_size() const
   return max_block_size;
 }
 
+template < typename value_type_ >
+inline typename BlockVector< value_type_ >::size_type
+BlockVector< value_type_ >::max_size() const
+{
+  throw NotImplemented( "BlockVector max_size() is not implemented." );
+}
+
+template < typename value_type_ >
+inline typename BlockVector< value_type_ >::reverse_iterator
+BlockVector< value_type_ >::rbegin()
+{
+  throw NotImplemented( "BlockVector rbegin() is not implemented." );
+}
+
+template < typename value_type_ >
+inline typename BlockVector< value_type_ >::reverse_iterator
+BlockVector< value_type_ >::rbegin() const
+{
+  throw NotImplemented( "BlockVector rbegin() is not implemented." );
+}
+
+template < typename value_type_ >
+inline typename BlockVector< value_type_ >::reverse_iterator
+BlockVector< value_type_ >::rend()
+{
+  throw NotImplemented( "BlockVector rend() is not implemented." );
+}
+
+template < typename value_type_ >
+inline typename BlockVector< value_type_ >::reverse_iterator
+BlockVector< value_type_ >::rend() const
+{
+  throw NotImplemented( "BlockVector rend() is not implemented." );
+}
+
 /////////////////////////////////////////////////////////////
 //        BlockVector iterator method implementation       //
 /////////////////////////////////////////////////////////////
@@ -531,6 +614,10 @@ inline bv_iterator< value_type_, ref_, ptr_ >& bv_iterator< value_type_, ref_, p
 template < typename value_type_, typename ref_, typename ptr_ >
 inline bv_iterator< value_type_, ref_, ptr_ >& bv_iterator< value_type_, ref_, ptr_ >::operator+=( difference_type val )
 {
+  if ( val < 0 )
+  {
+    return operator-=( -val );
+  }
   for ( difference_type i = 0; i < val; ++i )
   {
     operator++();
@@ -541,6 +628,10 @@ inline bv_iterator< value_type_, ref_, ptr_ >& bv_iterator< value_type_, ref_, p
 template < typename value_type_, typename ref_, typename ptr_ >
 inline bv_iterator< value_type_, ref_, ptr_ >& bv_iterator< value_type_, ref_, ptr_ >::operator-=( difference_type val )
 {
+  if ( val < 0 )
+  {
+    return operator+=( -val );
+  }
   for ( difference_type i = 0; i < val; ++i )
   {
     operator--();
@@ -549,14 +640,16 @@ inline bv_iterator< value_type_, ref_, ptr_ >& bv_iterator< value_type_, ref_, p
 }
 
 template < typename value_type_, typename ref_, typename ptr_ >
-inline bv_iterator< value_type_, ref_, ptr_ > bv_iterator< value_type_, ref_, ptr_ >::operator+( difference_type val )
+inline bv_iterator< value_type_, ref_, ptr_ > bv_iterator< value_type_, ref_, ptr_ >::operator+(
+  difference_type val ) const
 {
   bv_iterator tmp = *this;
   return tmp += val;
 }
 
 template < typename value_type_, typename ref_, typename ptr_ >
-inline bv_iterator< value_type_, ref_, ptr_ > bv_iterator< value_type_, ref_, ptr_ >::operator-( difference_type val )
+inline bv_iterator< value_type_, ref_, ptr_ > bv_iterator< value_type_, ref_, ptr_ >::operator-(
+  difference_type val ) const
 {
   bv_iterator tmp = *this;
   return tmp -= val;
@@ -570,6 +663,14 @@ inline bv_iterator< value_type_, ref_, ptr_ > bv_iterator< value_type_, ref_, pt
   return old;
 }
 
+template < typename value_type_, typename ref_, typename ptr_ >
+inline bv_iterator< value_type_, ref_, ptr_ > bv_iterator< value_type_, ref_, ptr_ >::operator--( int )
+{
+  bv_iterator< value_type_, ref_, ptr_ > old( *this );
+  --( *this );
+  return old;
+}
+
 template < typename value_type_, typename ref_, typename ptr_ >
 inline typename bv_iterator< value_type_, ref_, ptr_ >::reference bv_iterator< value_type_, ref_, ptr_ >::
 operator*() const
@@ -606,6 +707,13 @@ operator-( const const_iterator& other ) const
   return ( block_index_ - other.block_index_ ) * max_block_size + ( this_element_index - other_element_index );
 }
 
+template < typename value_type_, typename ref_, typename ptr_ >
+inline typename bv_iterator< value_type_, ref_, ptr_ >::reference bv_iterator< value_type_, ref_, ptr_ >::operator[](
+  difference_type n ) const
+{
+  return *( *this + n );
+}
+
 template < typename value_type_, typename ref_, typename ptr_ >
 inline bool bv_iterator< value_type_, ref_, ptr_ >::operator==(
   const bv_iterator< value_type_, ref_, ptr_ >& rhs ) const
@@ -626,6 +734,24 @@ inline bool bv_iterator< value_type_, ref_, ptr_ >::operator<( const bv_iterator
   return ( block_index_ < rhs.block_index_ or ( block_index_ == rhs.block_index_ and block_it_ < rhs.block_it_ ) );
 }
 
+template < typename value_type_, typename ref_, typename ptr_ >
+inline bool bv_iterator< value_type_, ref_, ptr_ >::operator>( const bv_iterator& rhs ) const
+{
+  return ( block_index_ > rhs.block_index_ or ( block_index_ == rhs.block_index_ and block_it_ > rhs.block_it_ ) );
+}
+
+template < typename value_type_, typename ref_, typename ptr_ >
+inline bool bv_iterator< value_type_, ref_, ptr_ >::operator<=( const bv_iterator& rhs ) const
+{
+  return operator<( rhs ) or operator==( rhs );
+}
+
+template < typename value_type_, typename ref_, typename ptr_ >
+inline bool bv_iterator< value_type_, ref_, ptr_ >::operator>=( const bv_iterator& rhs ) const
+{
+  return operator>( rhs ) or operator==( rhs );
+}
+
 template < typename value_type_, typename ref_, typename ptr_ >
 inline typename bv_iterator< value_type_, ref_, ptr_ >::iterator
 bv_iterator< value_type_, ref_, ptr_ >::const_cast_() const
@@ -633,4 +759,12 @@ bv_iterator< value_type_, ref_, ptr_ >::const_cast_() const
   return iterator( block_vector_, block_index_, block_it_, current_block_end_ );
 }
 
+template < typename value_type_, typename ref_, typename ptr_ >
+inline bv_iterator< value_type_, ref_, ptr_ > operator+(
+  typename bv_iterator< value_type_, ref_, ptr_ >::difference_type n,
+  bv_iterator< value_type_, ref_, ptr_ >& x )
+{
+  return x + n;
+}
+
 #endif /* BLOCK_VECTOR_H_ */
diff --git a/testsuite/cpptests/run_all.cpp b/testsuite/cpptests/run_all.cpp
index 7aa23d740b..24650185d9 100644
--- a/testsuite/cpptests/run_all.cpp
+++ b/testsuite/cpptests/run_all.cpp
@@ -22,7 +22,7 @@
 
 #define BOOST_TEST_MODULE cpptests
 #define BOOST_TEST_DYN_LINK
-#include <boost/test/unit_test.hpp>
+#include <boost/test/included/unit_test.hpp>
 
 // Includes from cpptests
 #include "test_block_vector.h"
diff --git a/testsuite/cpptests/test_block_vector.h b/testsuite/cpptests/test_block_vector.h
index 439c8f9625..ba6e9cc507 100644
--- a/testsuite/cpptests/test_block_vector.h
+++ b/testsuite/cpptests/test_block_vector.h
@@ -32,7 +32,31 @@
 // Includes from libnestutil:
 #include "block_vector.h"
 
-BOOST_AUTO_TEST_SUITE( test_seque )
+/**
+ * Fixture filling a BlockVector and a vector with linearly increasing values.
+ */
+struct bv_vec_reference_fixture
+{
+  bv_vec_reference_fixture()
+    : N( block_vector.get_max_block_size() + 10 )
+  {
+    for ( int i = 0; i < N; ++i )
+    {
+      block_vector.push_back( i );
+      reference.push_back( i );
+    }
+  }
+
+  ~bv_vec_reference_fixture()
+  {
+  }
+
+  BlockVector< int > block_vector;
+  std::vector< int > reference;
+  int N;
+};
+
+BOOST_AUTO_TEST_SUITE( test_blockvector )
 
 BOOST_AUTO_TEST_CASE( test_size )
 {
@@ -98,27 +122,27 @@ BOOST_AUTO_TEST_CASE( test_erase )
     block_vector.push_back( i );
   }
 
-  auto seque_mid = block_vector;
-  seque_mid.erase( seque_mid.begin() + 2, seque_mid.begin() + 8 );
-  BOOST_REQUIRE( seque_mid.size() == 4 );
-  BOOST_REQUIRE( seque_mid[ 0 ] == 0 );
-  BOOST_REQUIRE( seque_mid[ 1 ] == 1 );
-  BOOST_REQUIRE( seque_mid[ 2 ] == 8 );
-  BOOST_REQUIRE( seque_mid[ 3 ] == 9 );
-
-  auto seque_front = block_vector;
-  seque_front.erase( seque_front.begin(), seque_front.begin() + 7 );
-  BOOST_REQUIRE( seque_front.size() == 3 );
-  BOOST_REQUIRE( seque_front[ 0 ] == 7 );
-  BOOST_REQUIRE( seque_front[ 1 ] == 8 );
-  BOOST_REQUIRE( seque_front[ 2 ] == 9 );
-
-  auto seque_back = block_vector;
-  seque_back.erase( seque_back.begin() + 3, seque_back.end() );
-  BOOST_REQUIRE( seque_back.size() == 3 );
-  BOOST_REQUIRE( seque_back[ 0 ] == 0 );
-  BOOST_REQUIRE( seque_back[ 1 ] == 1 );
-  BOOST_REQUIRE( seque_back[ 2 ] == 2 );
+  auto bv_mid = block_vector;
+  bv_mid.erase( bv_mid.begin() + 2, bv_mid.begin() + 8 );
+  BOOST_REQUIRE( bv_mid.size() == 4 );
+  BOOST_REQUIRE( bv_mid[ 0 ] == 0 );
+  BOOST_REQUIRE( bv_mid[ 1 ] == 1 );
+  BOOST_REQUIRE( bv_mid[ 2 ] == 8 );
+  BOOST_REQUIRE( bv_mid[ 3 ] == 9 );
+
+  auto bv_front = block_vector;
+  bv_front.erase( bv_front.begin(), bv_front.begin() + 7 );
+  BOOST_REQUIRE( bv_front.size() == 3 );
+  BOOST_REQUIRE( bv_front[ 0 ] == 7 );
+  BOOST_REQUIRE( bv_front[ 1 ] == 8 );
+  BOOST_REQUIRE( bv_front[ 2 ] == 9 );
+
+  auto bv_back = block_vector;
+  bv_back.erase( bv_back.begin() + 3, bv_back.end() );
+  BOOST_REQUIRE( bv_back.size() == 3 );
+  BOOST_REQUIRE( bv_back[ 0 ] == 0 );
+  BOOST_REQUIRE( bv_back[ 1 ] == 1 );
+  BOOST_REQUIRE( bv_back[ 2 ] == 2 );
 }
 
 BOOST_AUTO_TEST_CASE( test_begin )
@@ -209,10 +233,10 @@ BOOST_AUTO_TEST_CASE( test_iterator_dereference )
   BOOST_REQUIRE( *( block_vector.begin() ) == block_vector[ 0 ] );
 
   // Using operator->
-  BlockVector< std::vector< int > > nested_seque;
+  BlockVector< std::vector< int > > nested_bv;
   std::vector< int > tmp = { 42 };
-  nested_seque.push_back( tmp );
-  BOOST_REQUIRE( nested_seque.begin()->size() == 1 );
+  nested_bv.push_back( tmp );
+  BOOST_REQUIRE( nested_bv.begin()->size() == 1 );
 }
 
 BOOST_AUTO_TEST_CASE( test_iterator_assign )
@@ -250,16 +274,145 @@ BOOST_AUTO_TEST_CASE( test_iterator_compare )
   }
   BOOST_REQUIRE( block_vector.begin() < block_vector.end() );
 
-  auto it_a = block_vector.begin();
-  auto it_b = block_vector.begin();
-  BOOST_REQUIRE( it_a == it_b );
+  // Test comparison with iterator shifted one step, shifted to the end of
+  // the block, and shifted to the next block.
+  std::vector< int > it_shifts = { 1, block_vector.get_max_block_size() - 1, N - 1 };
+  for ( auto& shift : it_shifts )
+  {
+    auto it_a = block_vector.begin();
+    auto it_b = block_vector.begin();
+    BOOST_REQUIRE( it_a == it_b );
+    BOOST_REQUIRE( not( it_a != it_b ) );
+
+    it_b += shift;
+
+    BOOST_REQUIRE( it_a != it_b );
+    BOOST_REQUIRE( it_a < it_b );
+    BOOST_REQUIRE( it_a <= it_b );
+    BOOST_REQUIRE( it_b > it_a );
+    BOOST_REQUIRE( it_b >= it_a );
+
+    BOOST_REQUIRE( not( it_a == it_b ) );
+    BOOST_REQUIRE( not( it_b < it_a ) );
+    BOOST_REQUIRE( not( it_b <= it_a ) );
+    BOOST_REQUIRE( not( it_a > it_b ) );
+    BOOST_REQUIRE( not( it_a >= it_b ) );
+  }
+}
 
-  ++it_b;
-  BOOST_REQUIRE( it_a != it_b );
-  BOOST_REQUIRE( it_a < it_b );
-  BOOST_REQUIRE( not( it_b < it_a ) );
+BOOST_FIXTURE_TEST_CASE( test_operator_pp, bv_vec_reference_fixture )
+{
+  // operator++()
+  auto bvi = block_vector.begin();
+  for ( auto& ref : reference )
+  {
+    BOOST_REQUIRE( *bvi == ref );
+    ++bvi;
+  }
+  BOOST_REQUIRE( bvi == block_vector.end() );
+}
+
+BOOST_FIXTURE_TEST_CASE( test_operator_p, bv_vec_reference_fixture )
+{
+  // operator+()
+  for ( int i = 0; i < N; ++i )
+  {
+    auto bvi = block_vector.begin();
+    auto ref_it = reference.begin();
+    auto new_bvi = bvi + i;
+    auto new_ref_it = ref_it + i;
+    BOOST_REQUIRE( *new_bvi == *new_ref_it );
+  }
+}
+
+BOOST_FIXTURE_TEST_CASE( test_operator_p_eq, bv_vec_reference_fixture )
+{
+  // operator+=()
+  for ( int i = 0; i < N; ++i )
+  {
+    auto bvi = block_vector.begin();
+    auto bvi_last = block_vector.end() - 1;
+    auto ref_it = reference.begin();
+    auto ref_it_last = reference.end() - 1;
+    bvi += i;
+    ref_it += i;
+    BOOST_REQUIRE( *bvi == *ref_it );
+    bvi_last += -i;
+    ref_it_last += -i;
+    BOOST_REQUIRE( *bvi_last == *ref_it_last );
+  }
+}
+
+BOOST_FIXTURE_TEST_CASE( test_operator_m_eq, bv_vec_reference_fixture )
+{
+  // operator-=()
+  for ( int i = 1; i < N - 1; ++i )
+  {
+    auto bvi = block_vector.end();
+    auto bvi_first = block_vector.begin();
+    auto ref_it = reference.end();
+    auto ref_it_first = reference.begin();
+    bvi -= i;
+    ref_it -= i;
+    BOOST_REQUIRE( *bvi == *ref_it );
+    bvi_first -= -( i - 1 );
+    ref_it_first -= -( i - 1 );
+    BOOST_REQUIRE( *bvi_first == *ref_it_first );
+  }
+}
+
+BOOST_FIXTURE_TEST_CASE( test_operator_m, bv_vec_reference_fixture )
+{
+  // operator-()
+  for ( int i = 1; i < N - 1; ++i )
+  {
+    auto bvi = block_vector.end();
+    auto ref_it = reference.end();
+    auto new_bvi = bvi - i;
+    auto new_ref_it = ref_it - i;
+    BOOST_REQUIRE( *new_bvi == *new_ref_it );
+  }
+}
+
+BOOST_FIXTURE_TEST_CASE( test_operator_mm, bv_vec_reference_fixture )
+{
+  // operator--()
+  auto bvi = block_vector.end() - 1;
+  for ( auto ref_it = reference.end() - 1; ref_it >= reference.begin(); --ref_it, --bvi )
+  {
+    BOOST_REQUIRE( *bvi == *ref_it );
+  }
+  BOOST_REQUIRE( bvi == block_vector.begin() - 1 );
+}
+
+BOOST_FIXTURE_TEST_CASE( test_operator_eq, bv_vec_reference_fixture )
+{
+  // operator==()
+  auto bvi_pp = block_vector.begin() + 1;
+  auto bvi_copy = block_vector.begin();
+  auto bvi_mm = block_vector.begin() - 1;
+  for ( auto bvi = block_vector.begin(); bvi != block_vector.end(); ++bvi, ++bvi_copy, ++bvi_mm, ++bvi_pp )
+  {
+    BOOST_REQUIRE( bvi == bvi_copy );
+    BOOST_REQUIRE( not( bvi == bvi_pp ) );
+    BOOST_REQUIRE( not( bvi == bvi_mm ) );
+  }
+}
+
+BOOST_FIXTURE_TEST_CASE( test_operator_neq, bv_vec_reference_fixture )
+{
+  // operator!=()
+  auto bvi_pp = block_vector.begin() + 1;
+  auto bvi_copy = block_vector.begin();
+  auto bvi_mm = block_vector.begin() - 1;
+  for ( auto bvi = block_vector.begin(); bvi != block_vector.end(); ++bvi, ++bvi_copy, ++bvi_mm, ++bvi_pp )
+  {
+    BOOST_REQUIRE( not( bvi != bvi_copy ) );
+    BOOST_REQUIRE( bvi != bvi_pp );
+    BOOST_REQUIRE( bvi != bvi_mm );
+  }
 }
 
 BOOST_AUTO_TEST_SUITE_END()
 
-#endif /* TEST_SORT_H */
+#endif /* TEST_BLOCK_VECTOR_H */
diff --git a/testsuite/cpptests/test_sort.h b/testsuite/cpptests/test_sort.h
index f1f879c798..67ebcf1a65 100644
--- a/testsuite/cpptests/test_sort.h
+++ b/testsuite/cpptests/test_sort.h
@@ -27,82 +27,211 @@
 #include <boost/test/unit_test.hpp>
 
 // C++ includes:
+#include <algorithm>
 #include <vector>
 
 // Includes from libnestutil:
 #include "sort.h"
 
-namespace nest
-{
-
+/**
+ * Wrapper for quicksort3way.
+ *
+ * When calling nest::sort() directly it is impossible to sort with the
+ * built in quicksort3way from tests. This is because if NEST is compiled
+ * with Boost support, nest::sort() sorts with Boost, and if NEST is not
+ * compiled with Boost, nest::sort() calls quicksort3way, but the Boost
+ * tests are not compiled.
+ */
 void
-nest_quicksort( BlockVector< size_t >& bv0, BlockVector< size_t >& bv1 )
+nest_quicksort( BlockVector< int >& bv0, BlockVector< int >& bv1 )
 {
   nest::quicksort3way( bv0, bv1, 0, bv0.size() - 1 );
 }
 
-bool
-is_sorted( BlockVector< size_t >::const_iterator begin, BlockVector< size_t >::const_iterator end )
+/**
+ * Fixture filling two BlockVectors and a vector with lineary decreasing numbers.
+ * The vector is then sorted.
+ */
+struct fill_bv_vec_linear
 {
-  for ( BlockVector< size_t >::const_iterator it = begin; it < --end; )
+  fill_bv_vec_linear()
+    : N( 20000 )
+    , N_small( boost::sort::spreadsort::detail::min_sort_size - 10 )
+    , bv_sort( N )
+    , bv_perm( N )
+    , vec_sort( N )
+    , bv_sort_small( N_small )
+    , bv_perm_small( N_small )
+    , vec_sort_small( N_small )
   {
-    if ( *it > *( ++it ) )
+    for ( int i = 0; i < N; ++i )
     {
-      return false;
+      int element = N - i;
+      bv_sort[ i ] = element;
+      bv_perm[ i ] = element;
+      vec_sort[ i ] = element;
+
+      if ( i < N_small )
+      {
+        bv_sort_small[ i ] = element;
+        bv_perm_small[ i ] = element;
+        vec_sort_small[ i ] = element;
+      }
     }
+    std::sort( vec_sort.begin(), vec_sort.end() );
+    std::sort( vec_sort_small.begin(), vec_sort_small.end() );
   }
-  return true;
-}
+
+  const int N;
+  const int N_small;
+
+  BlockVector< int > bv_sort;
+  BlockVector< int > bv_perm;
+  std::vector< int > vec_sort;
+
+  BlockVector< int > bv_sort_small;
+  BlockVector< int > bv_perm_small;
+  std::vector< int > vec_sort_small;
+};
+
+/**
+ * Fixture filling two BlockVectors and a vector with random numbers.
+ * The vector is then sorted.
+ */
+struct fill_bv_vec_random
+{
+  fill_bv_vec_random()
+    : N( 20000 )
+    , N_small( boost::sort::spreadsort::detail::min_sort_size - 10 )
+    , bv_sort( N )
+    , bv_perm( N )
+    , vec_sort( N )
+    , bv_sort_small( N_small )
+    , bv_perm_small( N_small )
+    , vec_sort_small( N_small )
+  {
+    for ( int i = 0; i < N; ++i )
+    {
+      const size_t k = std::rand() % N;
+      bv_sort[ i ] = k;
+      bv_perm[ i ] = k;
+      vec_sort[ i ] = k;
+
+      if ( i < N_small )
+      {
+        bv_sort_small[ i ] = k;
+        bv_perm_small[ i ] = k;
+        vec_sort_small[ i ] = k;
+      }
+    }
+    std::sort( vec_sort.begin(), vec_sort.end() );
+    std::sort( vec_sort_small.begin(), vec_sort_small.end() );
+  }
+
+  const int N;
+  const int N_small;
+
+  BlockVector< int > bv_sort;
+  BlockVector< int > bv_perm;
+  std::vector< int > vec_sort;
+
+  BlockVector< int > bv_sort_small;
+  BlockVector< int > bv_perm_small;
+  std::vector< int > vec_sort_small;
+};
 
 BOOST_AUTO_TEST_SUITE( test_sort )
 
 /**
  * Tests whether two arrays with randomly generated numbers are sorted
- * correctly by a single call to sort.
+ * correctly when sorting with NEST's own quicksort.
+ */
+BOOST_FIXTURE_TEST_CASE( test_quicksort_random, fill_bv_vec_random )
+{
+  nest_quicksort( bv_sort, bv_perm );
+
+  // We test here that the BlockVectors bv_sort and bv_perm are both
+  // sorted here. However, if something goes wrong and they for example
+  // contain only zeros, is_sorted will also return true, but the test
+  // should not pass. Therefore, in addition to require the BlockVectors
+  // to be sorted, we also require them to be equal to the sorted
+  // reference vector, vec_sort.
+  BOOST_REQUIRE( std::is_sorted( bv_sort.begin(), bv_sort.end() ) );
+  BOOST_REQUIRE( std::is_sorted( bv_perm.begin(), bv_perm.end() ) );
+
+  BOOST_REQUIRE( std::equal( vec_sort.begin(), vec_sort.end(), bv_sort.begin() ) );
+  BOOST_REQUIRE( std::equal( vec_sort.begin(), vec_sort.end(), bv_perm.begin() ) );
+}
+
+/**
+ * Tests whether two arrays with linearly decreasing numbers are sorted
+ * correctly when sorting with the built-in quicksort.
+ */
+BOOST_FIXTURE_TEST_CASE( test_quicksort_linear, fill_bv_vec_linear )
+{
+  nest_quicksort( bv_sort, bv_perm );
+
+  BOOST_REQUIRE( std::is_sorted( bv_sort.begin(), bv_sort.end() ) );
+  BOOST_REQUIRE( std::is_sorted( bv_perm.begin(), bv_perm.end() ) );
+
+  BOOST_REQUIRE( std::equal( vec_sort.begin(), vec_sort.end(), bv_sort.begin() ) );
+  BOOST_REQUIRE( std::equal( vec_sort.begin(), vec_sort.end(), bv_perm.begin() ) );
+}
+
+/**
+ * Tests whether two arrays with randomly generated numbers are sorted
+ * correctly when sorting with Boost.
  */
-BOOST_AUTO_TEST_CASE( test_random )
+BOOST_FIXTURE_TEST_CASE( test_boost_random, fill_bv_vec_random )
 {
-  const size_t N = 20000;
-  BlockVector< size_t > bv0( N );
-  BlockVector< size_t > bv1( N );
+  // Making sure we are sorting with boost
+  static_assert( HAVE_BOOST, "Compiling Boost tests, but HAVE_BOOST!=1." );
 
-  for ( size_t i = 0; i < N; ++i )
-  {
-    const size_t k = std::rand() % N;
-    bv0[ i ] = k;
-    bv1[ i ] = k;
-  }
+  nest::sort( bv_sort, bv_perm );
+
+  BOOST_REQUIRE( std::is_sorted( bv_sort.begin(), bv_sort.end() ) );
+  BOOST_REQUIRE( std::is_sorted( bv_perm.begin(), bv_perm.end() ) );
 
-  nest_quicksort( bv0, bv1 );
+  BOOST_REQUIRE( std::equal( vec_sort.begin(), vec_sort.end(), bv_sort.begin() ) );
+  BOOST_REQUIRE( std::equal( vec_sort.begin(), vec_sort.end(), bv_perm.begin() ) );
 
-  BOOST_REQUIRE( is_sorted( bv0.begin(), bv0.end() ) );
-  BOOST_REQUIRE( is_sorted( bv1.begin(), bv1.end() ) );
+  // Using smaller data sets to sort with the fallback algorithm.
+  nest::sort( bv_sort_small, bv_perm_small );
+
+  BOOST_REQUIRE( std::is_sorted( bv_sort_small.begin(), bv_sort_small.end() ) );
+  BOOST_REQUIRE( std::is_sorted( bv_perm_small.begin(), bv_perm_small.end() ) );
+
+  BOOST_REQUIRE( std::equal( vec_sort_small.begin(), vec_sort_small.end(), bv_sort_small.begin() ) );
+  BOOST_REQUIRE( std::equal( vec_sort_small.begin(), vec_sort_small.end(), bv_perm_small.begin() ) );
 }
 
 /**
  * Tests whether two arrays with linearly increasing numbers are sorted
- * correctly by a single call to sort.
+ * correctly when sorting with Boost.
  */
-BOOST_AUTO_TEST_CASE( test_linear )
+BOOST_FIXTURE_TEST_CASE( test_boost_linear, fill_bv_vec_linear )
 {
-  const size_t N = 20000;
-  BlockVector< size_t > bv0( N );
-  BlockVector< size_t > bv1( N );
+  // Making sure we are sorting with boost
+  static_assert( HAVE_BOOST, "Compiling Boost tests, but HAVE_BOOST!=1." );
 
-  for ( size_t i = 0; i < N; ++i )
-  {
-    bv0[ i ] = N - i - 1;
-    bv1[ i ] = N - i - 1;
-  }
+  nest::sort( bv_sort, bv_perm );
+
+  BOOST_REQUIRE( std::is_sorted( bv_sort.begin(), bv_sort.end() ) );
+  BOOST_REQUIRE( std::is_sorted( bv_perm.begin(), bv_perm.end() ) );
+
+  BOOST_REQUIRE( std::equal( vec_sort.begin(), vec_sort.end(), bv_sort.begin() ) );
+  BOOST_REQUIRE( std::equal( vec_sort.begin(), vec_sort.end(), bv_perm.begin() ) );
 
-  nest_quicksort( bv0, bv1 );
+  // Using smaller data sets to sort with the fallback algorithm.
+  nest::sort( bv_sort_small, bv_perm_small );
 
-  BOOST_REQUIRE( is_sorted( bv0.begin(), bv0.end() ) );
-  BOOST_REQUIRE( is_sorted( bv1.begin(), bv1.end() ) );
+  BOOST_REQUIRE( std::is_sorted( bv_sort_small.begin(), bv_sort_small.end() ) );
+  BOOST_REQUIRE( std::is_sorted( bv_perm_small.begin(), bv_perm_small.end() ) );
+
+  BOOST_REQUIRE( std::equal( vec_sort_small.begin(), vec_sort_small.end(), bv_sort_small.begin() ) );
+  BOOST_REQUIRE( std::equal( vec_sort_small.begin(), vec_sort_small.end(), bv_perm_small.begin() ) );
 }
 
 BOOST_AUTO_TEST_SUITE_END()
 
-} // of namespace nest
-
 #endif /* TEST_SORT_H */