block_vector_base.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2004 - 2018 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE at
// the top level of the deal.II distribution.
//
// ---------------------------------------------------------------------

#ifndef dealii_block_vector_base_h
#define dealii_block_vector_base_h


#include <deal.II/base/config.h>
#include <deal.II/base/numbers.h>
#include <deal.II/base/exceptions.h>
#include <deal.II/base/subscriptor.h>
#include <deal.II/base/memory_consumption.h>
#include <deal.II/lac/exceptions.h>
#include <deal.II/lac/block_indices.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/vector_operation.h>

#include <vector>
#include <iterator>
#include <cmath>
#include <cstddef>

DEAL_II_NAMESPACE_OPEN


template <typename> class BlockVectorBase;


template <typename VectorType>
struct IsBlockVector
{
private:
  struct yes_type
  {
    char c[1];
  };
  struct no_type
  {
    char c[2];
  };

  template <typename T>
  static yes_type check_for_block_vector (const BlockVectorBase<T> *);

  static no_type check_for_block_vector (...);

public:
  static const bool value = (sizeof(check_for_block_vector
                                    ((VectorType *)nullptr))
                             ==
                             sizeof(yes_type));
};


// instantiation of the static member
template <typename VectorType>
const bool IsBlockVector<VectorType>::value;




namespace internal
{

  namespace BlockVectorIterators
  {
    template <class BlockVectorType, bool Constness>
    struct Types
    {
    };



    template <class BlockVectorType>
    struct Types<BlockVectorType,false>
    {
      typedef typename BlockVectorType::BlockType Vector;

      typedef BlockVectorType BlockVector;

      typedef typename BlockVector::value_type value_type;

      typedef typename Vector::reference dereference_type;
    };



    template <class BlockVectorType>
    struct Types<BlockVectorType,true>
    {
      typedef const typename BlockVectorType::BlockType Vector;

      typedef const BlockVectorType BlockVector;

      typedef const typename BlockVector::value_type value_type;

      typedef value_type dereference_type;
    };


    template <class BlockVectorType, bool Constness>
    class Iterator
    {
    public:
      typedef types::global_dof_index size_type;

      typedef
      typename Types<BlockVectorType,Constness>::value_type
      value_type;

      typedef std::random_access_iterator_tag               iterator_category;
      typedef std::ptrdiff_t                                difference_type;
      typedef typename BlockVectorType::reference           reference;
      typedef value_type                                   *pointer;

      typedef
      typename Types<BlockVectorType,Constness>::dereference_type
      dereference_type;

      typedef
      typename Types<BlockVectorType,Constness>::BlockVector
      BlockVector;

      Iterator (BlockVector     &parent,
                const size_type  global_index);

      Iterator (const Iterator<BlockVectorType,!Constness> &c);


      Iterator (const Iterator &c);

    private:
      Iterator (BlockVector     &parent,
                const size_type  global_index,
                const size_type  current_block,
                const size_type  index_within_block,
                const size_type  next_break_forward,
                const size_type  next_break_backward);

    public:

      Iterator &operator = (const Iterator &c);

      dereference_type operator * () const;

      dereference_type operator [] (const difference_type d) const;

      Iterator &operator ++ ();

      Iterator operator ++ (int);

      Iterator &operator -- ();

      Iterator operator -- (int);

      template <bool OtherConstness>
      bool operator == (const Iterator<BlockVectorType, OtherConstness> &i) const;

      template <bool OtherConstness>
      bool operator != (const Iterator<BlockVectorType, OtherConstness> &i) const;

      template <bool OtherConstness>
      bool operator < (const Iterator<BlockVectorType, OtherConstness> &i) const;

      template <bool OtherConstness>
      bool operator <= (const Iterator<BlockVectorType, OtherConstness> &i) const;

      template <bool OtherConstness>
      bool operator > (const Iterator<BlockVectorType, OtherConstness> &i) const;

      template <bool OtherConstness>
      bool operator >= (const Iterator<BlockVectorType, OtherConstness> &i) const;

      template <bool OtherConstness>
      difference_type operator - (const Iterator<BlockVectorType, OtherConstness> &i) const;

      Iterator operator + (const difference_type &d) const;

      Iterator operator - (const difference_type &d) const;

      Iterator &operator += (const difference_type &d);

      Iterator &operator -= (const difference_type &d);

      DeclExceptionMsg (ExcPointerToDifferentVectors,
                        "Your program tried to compare iterators pointing to "
                        "different block vectors. There is no reasonable way "
                        "to do this.");

    private:
      BlockVector *parent;

      size_type     global_index;

      unsigned int current_block;
      size_type index_within_block;

      size_type next_break_forward;
      size_type next_break_backward;

      void move_forward ();

      void move_backward ();


      template <typename, bool>
      friend class Iterator;
    };
  }  // namespace BlockVectorIterators
}  // namespace internal


template <class VectorType>
class BlockVectorBase : public Subscriptor
{
public:
  typedef VectorType BlockType;

  /*
   * Declare standard types used in
   * all containers. These types
   * parallel those in the
   * <tt>C++</tt> standard
   * libraries
   * <tt>std::vector<...></tt>
   * class. This includes iterator
   * types.
   */
  typedef typename BlockType::value_type  value_type;
  typedef value_type                     *pointer;
  typedef const value_type               *const_pointer;
  typedef ::internal::BlockVectorIterators::Iterator<BlockVectorBase,false> iterator;
  typedef ::internal::BlockVectorIterators::Iterator<BlockVectorBase,true>  const_iterator;
  typedef typename BlockType::reference       reference;
  typedef typename BlockType::const_reference const_reference;
  typedef types::global_dof_index             size_type;

  typedef typename BlockType::real_type real_type;

  BlockVectorBase () = default;

  BlockVectorBase (const BlockVectorBase &/*V*/) = default;

  BlockVectorBase (BlockVectorBase &&/*V*/) noexcept = default;

  void collect_sizes ();

  void compress (::VectorOperation::values operation);

  BlockType &
  block (const unsigned int i);

  const BlockType &
  block (const unsigned int i) const;

  const BlockIndices &
  get_block_indices () const;

  unsigned int n_blocks () const;

  std::size_t size () const;

  IndexSet locally_owned_elements () const;

  iterator begin ();

  const_iterator begin () const;

  iterator end ();

  const_iterator end () const;

  value_type operator() (const size_type i) const;

  reference operator() (const size_type i);

  value_type operator[] (const size_type i) const;

  reference operator[] (const size_type i);

  template <typename OtherNumber>
  void extract_subvector_to (const std::vector<size_type> &indices,
                             std::vector<OtherNumber> &values) const;

  template <typename ForwardIterator, typename OutputIterator>
  void extract_subvector_to (ForwardIterator          indices_begin,
                             const ForwardIterator    indices_end,
                             OutputIterator           values_begin) const;

  BlockVectorBase &operator = (const value_type s);

  BlockVectorBase &
  operator= (const BlockVectorBase &V);

  BlockVectorBase &operator= (BlockVectorBase &&/*V*/) = default; // NOLINT

  template <class VectorType2>
  BlockVectorBase &
  operator= (const BlockVectorBase<VectorType2> &V);

  BlockVectorBase &
  operator = (const VectorType &v);

  template <class VectorType2>
  bool
  operator == (const BlockVectorBase<VectorType2> &v) const;

  value_type operator* (const BlockVectorBase &V) const;

  real_type norm_sqr () const;

  value_type mean_value () const;

  real_type l1_norm () const;

  real_type l2_norm () const;

  real_type linfty_norm () const;

  value_type add_and_dot (const value_type       a,
                          const BlockVectorBase &V,
                          const BlockVectorBase &W);

  bool in_local_range (const size_type global_index) const;

  bool all_zero () const;

  bool is_non_negative () const;

  BlockVectorBase &
  operator += (const BlockVectorBase &V);

  BlockVectorBase &
  operator -= (const BlockVectorBase &V);


  template <typename Number>
  void add (const std::vector<size_type> &indices,
            const std::vector<Number>    &values);

  template <typename Number>
  void add (const std::vector<size_type> &indices,
            const Vector<Number>         &values);

  template <typename Number>
  void add (const size_type  n_elements,
            const size_type *indices,
            const Number    *values);

  void add (const value_type s);

  void add (const value_type a, const BlockVectorBase &V);

  void add (const value_type a, const BlockVectorBase &V,
            const value_type b, const BlockVectorBase &W);

  void sadd (const value_type s, const BlockVectorBase &V);

  void sadd (const value_type s, const value_type a, const BlockVectorBase &V);

  void sadd (const value_type s, const value_type a,
             const BlockVectorBase &V,
             const value_type b, const BlockVectorBase &W);

  void sadd (const value_type s, const value_type a,
             const BlockVectorBase &V,
             const value_type b, const BlockVectorBase &W,
             const value_type c, const BlockVectorBase &X);

  BlockVectorBase &operator *= (const value_type factor);

  BlockVectorBase &operator /= (const value_type factor);

  template <class BlockVector2>
  void scale (const BlockVector2 &v);

  template <class BlockVector2>
  void equ (const value_type a, const BlockVector2 &V);

  void equ (const value_type a, const BlockVectorBase &V,
            const value_type b, const BlockVectorBase &W);

  void update_ghost_values () const;

  std::size_t memory_consumption () const;

protected:
  std::vector<VectorType> components;

  BlockIndices block_indices;

  template <typename N, bool C>
  friend class ::internal::BlockVectorIterators::Iterator;

  template <typename> friend class BlockVectorBase;
};


/*----------------------- Inline functions ----------------------------------*/


#ifndef DOXYGEN
namespace internal
{
  namespace BlockVectorIterators
  {
    template <class BlockVectorType, bool Constness>
    inline
    Iterator<BlockVectorType,Constness>::
    Iterator (const Iterator<BlockVectorType,Constness> &c)
      :
      parent (c.parent),
      global_index (c.global_index),
      current_block (c.current_block),
      index_within_block (c.index_within_block),
      next_break_forward (c.next_break_forward),
      next_break_backward (c.next_break_backward)
    {}



    template <class BlockVectorType, bool Constness>
    inline
    Iterator<BlockVectorType,Constness>::
    Iterator (const Iterator<BlockVectorType,!Constness> &c)
      :
      parent (c.parent),
      global_index (c.global_index),
      current_block (c.current_block),
      index_within_block (c.index_within_block),
      next_break_forward (c.next_break_forward),
      next_break_backward (c.next_break_backward)
    {
      // Only permit copy-constructing const iterators from non-const
      // iterators, and not vice versa (i.e., Constness must always be
      // true).
      static_assert(Constness == true,
                    "Constructing a non-const iterator from a const iterator "
                    "does not make sense.");
    }



    template <class BlockVectorType, bool Constness>
    inline
    Iterator<BlockVectorType,Constness>::
    Iterator (BlockVector        &parent,
              const size_type  global_index,
              const size_type  current_block,
              const size_type  index_within_block,
              const size_type  next_break_forward,
              const size_type  next_break_backward)
      :
      parent (&parent),
      global_index (global_index),
      current_block (current_block),
      index_within_block (index_within_block),
      next_break_forward (next_break_forward),
      next_break_backward (next_break_backward)
    {
    }



    template <class BlockVectorType, bool Constness>
    inline
    Iterator<BlockVectorType,Constness> &
    Iterator<BlockVectorType,Constness>::
    operator = (const Iterator &c)
    {
      parent              = c.parent;
      global_index        = c.global_index;
      index_within_block  = c.index_within_block;
      current_block       = c.current_block;
      next_break_forward  = c.next_break_forward;
      next_break_backward = c.next_break_backward;

      return *this;
    }



    template <class BlockVectorType, bool Constness>
    inline
    typename Iterator<BlockVectorType,Constness>::dereference_type
    Iterator<BlockVectorType,Constness>::operator * () const
    {
      return parent->block(current_block)(index_within_block);
    }



    template <class BlockVectorType, bool Constness>
    inline
    typename Iterator<BlockVectorType,Constness>::dereference_type
    Iterator<BlockVectorType,Constness>::operator [] (const difference_type d) const
    {
      // if the index pointed to is
      // still within the block we
      // currently point into, then we
      // can save the computation of
      // the block
      if ((global_index+d >= next_break_backward) &&
          (global_index+d <= next_break_forward))
        return parent->block(current_block)(index_within_block + d);

      // if the index is not within the
      // block of the block vector into
      // which we presently point, then
      // there is no way: we have to
      // search for the block. this can
      // be done through the parent
      // class as well.
      return (*parent)(global_index+d);
    }



    template <class BlockVectorType, bool Constness>
    inline
    Iterator<BlockVectorType,Constness> &
    Iterator<BlockVectorType,Constness>::operator ++ ()
    {
      move_forward ();
      return *this;
    }



    template <class BlockVectorType, bool Constness>
    inline
    Iterator<BlockVectorType,Constness>
    Iterator<BlockVectorType,Constness>::operator ++ (int)
    {
      const Iterator old_value = *this;
      move_forward ();
      return old_value;
    }



    template <class BlockVectorType, bool Constness>
    inline
    Iterator<BlockVectorType,Constness> &
    Iterator<BlockVectorType,Constness>::operator -- ()
    {
      move_backward ();
      return *this;
    }



    template <class BlockVectorType, bool Constness>
    inline
    Iterator<BlockVectorType,Constness>
    Iterator<BlockVectorType,Constness>::operator -- (int)
    {
      const Iterator old_value = *this;
      move_backward ();
      return old_value;
    }



    template <class BlockVectorType, bool Constness>
    template <bool OtherConstness>
    inline
    bool
    Iterator<BlockVectorType,Constness>::
    operator == (const Iterator<BlockVectorType, OtherConstness> &i) const
    {
      Assert (parent == i.parent, ExcPointerToDifferentVectors());

      return (global_index == i.global_index);
    }



    template <class BlockVectorType, bool Constness>
    template <bool OtherConstness>
    inline
    bool
    Iterator<BlockVectorType,Constness>::
    operator != (const Iterator<BlockVectorType, OtherConstness> &i) const
    {
      Assert (parent == i.parent, ExcPointerToDifferentVectors());

      return (global_index != i.global_index);
    }



    template <class BlockVectorType, bool Constness>
    template <bool OtherConstness>
    inline
    bool
    Iterator<BlockVectorType,Constness>::
    operator < (const Iterator<BlockVectorType, OtherConstness> &i) const
    {
      Assert (parent == i.parent, ExcPointerToDifferentVectors());

      return (global_index < i.global_index);
    }



    template <class BlockVectorType, bool Constness>
    template <bool OtherConstness>
    inline
    bool
    Iterator<BlockVectorType,Constness>::
    operator <= (const Iterator<BlockVectorType, OtherConstness> &i) const
    {
      Assert (parent == i.parent, ExcPointerToDifferentVectors());

      return (global_index <= i.global_index);
    }



    template <class BlockVectorType, bool Constness>
    template <bool OtherConstness>
    inline
    bool
    Iterator<BlockVectorType,Constness>::
    operator > (const Iterator<BlockVectorType, OtherConstness> &i) const
    {
      Assert (parent == i.parent, ExcPointerToDifferentVectors());

      return (global_index > i.global_index);
    }



    template <class BlockVectorType, bool Constness>
    template <bool OtherConstness>
    inline
    bool
    Iterator<BlockVectorType,Constness>::
    operator >= (const Iterator<BlockVectorType, OtherConstness> &i) const
    {
      Assert (parent == i.parent, ExcPointerToDifferentVectors());

      return (global_index >= i.global_index);
    }



    template <class BlockVectorType, bool Constness>
    template <bool OtherConstness>
    inline
    typename Iterator<BlockVectorType,Constness>::difference_type
    Iterator<BlockVectorType,Constness>::
    operator - (const Iterator<BlockVectorType, OtherConstness> &i) const
    {
      Assert (parent == i.parent, ExcPointerToDifferentVectors());

      return (static_cast<signed int>(global_index) -
              static_cast<signed int>(i.global_index));
    }



    template <class BlockVectorType, bool Constness>
    inline
    Iterator<BlockVectorType,Constness>
    Iterator<BlockVectorType,Constness>::
    operator + (const difference_type &d) const
    {
      // if the index pointed to is
      // still within the block we
      // currently point into, then we
      // can save the computation of
      // the block
      if ((global_index+d >= next_break_backward) &&
          (global_index+d <= next_break_forward))
        return Iterator (*parent, global_index+d, current_block,
                         index_within_block+d,
                         next_break_forward, next_break_backward);
      else
        // outside present block, so
        // have to seek new block
        // anyway
        return Iterator (*parent, global_index+d);
    }



    template <class BlockVectorType, bool Constness>
    inline
    Iterator<BlockVectorType,Constness>
    Iterator<BlockVectorType,Constness>::
    operator - (const difference_type &d) const
    {
      // if the index pointed to is
      // still within the block we
      // currently point into, then we
      // can save the computation of
      // the block
      if ((global_index-d >= next_break_backward) &&
          (global_index-d <= next_break_forward))
        return Iterator (*parent, global_index-d, current_block,
                         index_within_block-d,
                         next_break_forward, next_break_backward);
      else
        // outside present block, so
        // have to seek new block
        // anyway
        return Iterator (*parent, global_index-d);
    }



    template <class BlockVectorType, bool Constness>
    inline
    Iterator<BlockVectorType,Constness> &
    Iterator<BlockVectorType,Constness>::
    operator += (const difference_type &d)
    {
      // if the index pointed to is
      // still within the block we
      // currently point into, then we
      // can save the computation of
      // the block
      if ((global_index+d >= next_break_backward) &&
          (global_index+d <= next_break_forward))
        {
          global_index       += d;
          index_within_block += d;
        }
      else
        // outside present block, so
        // have to seek new block
        // anyway
        *this = Iterator (*parent, global_index+d);

      return *this;
    }



    template <class BlockVectorType, bool Constness>
    inline
    Iterator<BlockVectorType,Constness> &
    Iterator<BlockVectorType,Constness>::
    operator -= (const difference_type &d)
    {
      // if the index pointed to is
      // still within the block we
      // currently point into, then we
      // can save the computation of
      // the block
      if ((global_index-d >= next_break_backward) &&
          (global_index-d <= next_break_forward))
        {
          global_index       -= d;
          index_within_block -= d;
        }
      else
        // outside present block, so
        // have to seek new block
        // anyway
        *this = Iterator (*parent, global_index-d);

      return *this;
    }


    template <class BlockVectorType, bool Constness>
    Iterator<BlockVectorType,Constness>::
    Iterator (BlockVector    &parent,
              const size_type global_index)
      :
      parent (&parent),
      global_index (global_index)
    {
      // find which block we are
      // in. for this, take into
      // account that it happens at
      // times that people want to
      // initialize iterators
      // past-the-end
      if (global_index < parent.size())
        {
          const std::pair<size_type, size_type>
          indices = parent.block_indices.global_to_local(global_index);
          current_block      = indices.first;
          index_within_block = indices.second;

          next_break_backward
            = parent.block_indices.local_to_global (current_block, 0);
          next_break_forward
            = (parent.block_indices.local_to_global (current_block, 0)
               +parent.block_indices.block_size(current_block)-1);
        }
      else
        // past the end. only have one
        // value for this
        {
          this->global_index  = parent.size ();
          current_block       = parent.n_blocks();
          index_within_block  = 0;
          next_break_backward = global_index;
          next_break_forward  = numbers::invalid_size_type;
        };
    }



    template <class BlockVectorType, bool Constness>
    void
    Iterator<BlockVectorType,Constness>::move_forward ()
    {
      if (global_index != next_break_forward)
        ++index_within_block;
      else
        {
          // ok, we traverse a boundary
          // between blocks:
          index_within_block = 0;
          ++current_block;

          // break backwards is now old
          // break forward
          next_break_backward = next_break_forward+1;

          // compute new break forward
          if (current_block < parent->block_indices.size())
            next_break_forward
            += parent->block_indices.block_size(current_block);
          else
            // if we are beyond the end,
            // then move the next
            // boundary arbitrarily far
            // away
            next_break_forward = numbers::invalid_size_type;
        };

      ++global_index;
    }



    template <class BlockVectorType, bool Constness>
    void
    Iterator<BlockVectorType,Constness>::move_backward ()
    {
      if (global_index != next_break_backward)
        --index_within_block;
      else if (current_block != 0)
        {
          // ok, we traverse a boundary
          // between blocks:
          --current_block;
          index_within_block = parent->block_indices.block_size(current_block)-1;

          // break forwards is now old
          // break backward
          next_break_forward = next_break_backward-1;

          // compute new break forward
          next_break_backward
          -= parent->block_indices.block_size (current_block);
        }
      else
        // current block was 0, we now
        // get into unspecified terrain
        {
          --current_block;
          index_within_block = numbers::invalid_size_type;
          next_break_forward = 0;
          next_break_backward = 0;
        };

      --global_index;
    }


  } // namespace BlockVectorIterators

} //namespace internal



template <class VectorType>
inline
std::size_t
BlockVectorBase<VectorType>::size () const
{
  return block_indices.total_size();
}



template <class VectorType>
inline
IndexSet
BlockVectorBase<VectorType>::locally_owned_elements () const
{
  IndexSet is (size());

  // copy index sets from blocks into the global one, shifted
  // by the appropriate amount for each block
  for (unsigned int b=0; b<n_blocks(); ++b)
    {
      IndexSet x = block(b).locally_owned_elements();
      is.add_indices(x, block_indices.block_start(b));
    }

  is.compress();

  return is;
}



template <class VectorType>
inline
unsigned int
BlockVectorBase<VectorType>::n_blocks () const
{
  return block_indices.size();
}


template <class VectorType>
inline
typename BlockVectorBase<VectorType>::BlockType &
BlockVectorBase<VectorType>::block (const unsigned int i)
{
  Assert(i<n_blocks(), ExcIndexRange(i,0,n_blocks()));

  return components[i];
}



template <class VectorType>
inline
const typename BlockVectorBase<VectorType>::BlockType &
BlockVectorBase<VectorType>::block (const unsigned int i) const
{
  Assert(i<n_blocks(), ExcIndexRange(i,0,n_blocks()));

  return components[i];
}



template <class VectorType>
inline
const BlockIndices &
BlockVectorBase<VectorType>::get_block_indices () const
{
  return block_indices;
}


template <class VectorType>
inline
void
BlockVectorBase<VectorType>::collect_sizes ()
{
  std::vector<size_type> sizes (n_blocks());

  for (size_type i=0; i<n_blocks(); ++i)
    sizes[i] = block(i).size();

  block_indices.reinit(sizes);
}



template <class VectorType>
inline
void
BlockVectorBase<VectorType>::compress (::VectorOperation::values operation)
{
  for (unsigned int i=0; i<n_blocks(); ++i)
    block(i).compress (operation);
}



template <class VectorType>
inline
typename BlockVectorBase<VectorType>::iterator
BlockVectorBase<VectorType>::begin()
{
  return iterator(*this, 0U);
}



template <class VectorType>
inline
typename BlockVectorBase<VectorType>::const_iterator
BlockVectorBase<VectorType>::begin() const
{
  return const_iterator(*this, 0U);
}


template <class VectorType>
inline
typename BlockVectorBase<VectorType>::iterator
BlockVectorBase<VectorType>::end()
{
  return iterator(*this, size());
}



template <class VectorType>
inline
typename BlockVectorBase<VectorType>::const_iterator
BlockVectorBase<VectorType>::end() const
{
  return const_iterator(*this, size());
}


template <class VectorType>
inline
bool
BlockVectorBase<VectorType>::in_local_range
(const size_type global_index) const
{
  const std::pair<size_type,size_type> local_index
    = block_indices.global_to_local (global_index);

  return components[local_index.first].in_local_range (global_index);
}


template <class VectorType>
bool
BlockVectorBase<VectorType>::all_zero () const
{
  for (size_type i=0; i<n_blocks(); ++i)
    if (components[i].all_zero() == false)
      return false;

  return true;
}



template <class VectorType>
bool
BlockVectorBase<VectorType>::is_non_negative () const
{
  for (size_type i=0; i<n_blocks(); ++i)
    if (components[i].is_non_negative() == false)
      return false;

  return true;
}



template <class VectorType>
typename BlockVectorBase<VectorType>::value_type
BlockVectorBase<VectorType>::
operator * (const BlockVectorBase<VectorType> &v) const
{
  Assert (n_blocks() == v.n_blocks(),
          ExcDimensionMismatch(n_blocks(), v.n_blocks()));

  value_type sum = 0.;
  for (size_type i=0; i<n_blocks(); ++i)
    sum += components[i]*v.components[i];

  return sum;
}


template <class VectorType>
typename BlockVectorBase<VectorType>::real_type
BlockVectorBase<VectorType>::norm_sqr () const
{
  real_type sum = 0.;
  for (size_type i=0; i<n_blocks(); ++i)
    sum += components[i].norm_sqr();

  return sum;
}



template <class VectorType>
typename BlockVectorBase<VectorType>::value_type
BlockVectorBase<VectorType>::mean_value () const
{
  value_type sum = 0.;
  // need to do static_cast as otherwise it won't work with value_type=complex<T>
  for (size_type i=0; i<n_blocks(); ++i)
    sum += components[i].mean_value() * (typename numbers::NumberTraits<value_type>::real_type(components[i].size()));

  return sum/(typename numbers::NumberTraits<value_type>::real_type(size()));
}



template <class VectorType>
typename BlockVectorBase<VectorType>::real_type
BlockVectorBase<VectorType>::l1_norm () const
{
  real_type sum = 0.;
  for (size_type i=0; i<n_blocks(); ++i)
    sum += components[i].l1_norm();

  return sum;
}



template <class VectorType>
typename BlockVectorBase<VectorType>::real_type
BlockVectorBase<VectorType>::l2_norm () const
{
  return std::sqrt(norm_sqr());
}



template <class VectorType>
typename BlockVectorBase<VectorType>::real_type
BlockVectorBase<VectorType>::linfty_norm () const
{
  real_type sum = 0.;
  for (size_type i=0; i<n_blocks(); ++i)
    {
      value_type newval = components[i].linfty_norm();
      if (sum<newval)
        sum = newval;
    }
  return sum;
}



template <class VectorType>
typename BlockVectorBase<VectorType>::value_type
BlockVectorBase<VectorType>::
add_and_dot (const typename BlockVectorBase<VectorType>::value_type a,
             const BlockVectorBase<VectorType> &V,
             const BlockVectorBase<VectorType> &W)
{
  AssertDimension (n_blocks(), V.n_blocks());
  AssertDimension (n_blocks(), W.n_blocks());

  value_type sum = 0.;
  for (size_type i=0; i<n_blocks(); ++i)
    sum += components[i].add_and_dot(a, V.components[i], W.components[i]);

  return sum;
}



template <class VectorType>
BlockVectorBase<VectorType> &
BlockVectorBase<VectorType>::operator += (const BlockVectorBase<VectorType> &v)
{
  Assert (n_blocks() == v.n_blocks(),
          ExcDimensionMismatch(n_blocks(), v.n_blocks()));

  for (size_type i=0; i<n_blocks(); ++i)
    {
      components[i] += v.components[i];
    }

  return *this;
}



template <class VectorType>
BlockVectorBase<VectorType> &
BlockVectorBase<VectorType>::operator -= (const BlockVectorBase<VectorType> &v)
{
  Assert (n_blocks() == v.n_blocks(),
          ExcDimensionMismatch(n_blocks(), v.n_blocks()));

  for (size_type i=0; i<n_blocks(); ++i)
    {
      components[i] -= v.components[i];
    }
  return *this;
}



template <class VectorType>
template <typename Number>
inline
void
BlockVectorBase<VectorType>::add (const std::vector<size_type> &indices,
                                  const std::vector<Number>    &values)
{
  Assert (indices.size() == values.size(),
          ExcDimensionMismatch(indices.size(), values.size()));
  add (indices.size(), indices.data(), values.data());
}



template <class VectorType>
template <typename Number>
inline
void
BlockVectorBase<VectorType>::add (const std::vector<size_type> &indices,
                                  const Vector<Number>         &values)
{
  Assert (indices.size() == values.size(),
          ExcDimensionMismatch(indices.size(), values.size()));
  const size_type n_indices = indices.size();
  for (size_type i=0; i<n_indices; ++i)
    (*this)(indices[i]) += values(i);
}



template <class VectorType>
template <typename Number>
inline
void
BlockVectorBase<VectorType>::add (const size_type  n_indices,
                                  const size_type *indices,
                                  const Number    *values)
{
  for (size_type i=0; i<n_indices; ++i)
    (*this)(indices[i]) += values[i];
}



template <class VectorType>
void BlockVectorBase<VectorType>::add (const value_type a)
{
  AssertIsFinite(a);

  for (size_type i=0; i<n_blocks(); ++i)
    {
      components[i].add(a);
    }
}



template <class VectorType>
void BlockVectorBase<VectorType>::add (const value_type a,
                                       const BlockVectorBase<VectorType> &v)
{

  AssertIsFinite(a);

  Assert (n_blocks() == v.n_blocks(),
          ExcDimensionMismatch(n_blocks(), v.n_blocks()));

  for (size_type i=0; i<n_blocks(); ++i)
    {
      components[i].add(a, v.components[i]);
    }
}



template <class VectorType>
void BlockVectorBase<VectorType>::add (const value_type a,
                                       const BlockVectorBase<VectorType> &v,
                                       const value_type b,
                                       const BlockVectorBase<VectorType> &w)
{

  AssertIsFinite(a);
  AssertIsFinite(b);

  Assert (n_blocks() == v.n_blocks(),
          ExcDimensionMismatch(n_blocks(), v.n_blocks()));
  Assert (n_blocks() == w.n_blocks(),
          ExcDimensionMismatch(n_blocks(), w.n_blocks()));


  for (size_type i=0; i<n_blocks(); ++i)
    {
      components[i].add(a, v.components[i], b, w.components[i]);
    }
}



template <class VectorType>
void BlockVectorBase<VectorType>::sadd (const value_type x,
                                        const BlockVectorBase<VectorType> &v)
{

  AssertIsFinite(x);

  Assert (n_blocks() == v.n_blocks(),
          ExcDimensionMismatch(n_blocks(), v.n_blocks()));

  for (size_type i=0; i<n_blocks(); ++i)
    {
      components[i].sadd(x, v.components[i]);
    }
}



template <class VectorType>
void BlockVectorBase<VectorType>::sadd (const value_type x, const value_type a,
                                        const BlockVectorBase<VectorType> &v)
{

  AssertIsFinite(x);
  AssertIsFinite(a);

  Assert (n_blocks() == v.n_blocks(),
          ExcDimensionMismatch(n_blocks(), v.n_blocks()));

  for (size_type i=0; i<n_blocks(); ++i)
    {
      components[i].sadd(x, a, v.components[i]);
    }
}



template <class VectorType>
void BlockVectorBase<VectorType>::sadd (const value_type x, const value_type a,
                                        const BlockVectorBase<VectorType> &v,
                                        const value_type b,
                                        const BlockVectorBase<VectorType> &w)
{

  AssertIsFinite(x);
  AssertIsFinite(a);
  AssertIsFinite(b);

  Assert (n_blocks() == v.n_blocks(),
          ExcDimensionMismatch(n_blocks(), v.n_blocks()));
  Assert (n_blocks() == w.n_blocks(),
          ExcDimensionMismatch(n_blocks(), w.n_blocks()));

  for (size_type i=0; i<n_blocks(); ++i)
    {
      components[i].sadd(x, a, v.components[i], b, w.components[i]);
    }
}



template <class VectorType>
void BlockVectorBase<VectorType>::sadd (const value_type x, const value_type a,
                                        const BlockVectorBase<VectorType> &v,
                                        const value_type b,
                                        const BlockVectorBase<VectorType> &w,
                                        const value_type c,
                                        const BlockVectorBase<VectorType> &y)
{

  AssertIsFinite(x);
  AssertIsFinite(a);
  AssertIsFinite(b);
  AssertIsFinite(c);

  Assert (n_blocks() == v.n_blocks(),
          ExcDimensionMismatch(n_blocks(), v.n_blocks()));
  Assert (n_blocks() == w.n_blocks(),
          ExcDimensionMismatch(n_blocks(), w.n_blocks()));
  Assert (n_blocks() == y.n_blocks(),
          ExcDimensionMismatch(n_blocks(), y.n_blocks()));

  for (size_type i=0; i<n_blocks(); ++i)
    {
      components[i].sadd(x, a, v.components[i],
                         b, w.components[i], c, y.components[i]);
    }
}



template <class VectorType>
template <class BlockVector2>
void BlockVectorBase<VectorType>::scale (const BlockVector2 &v)
{
  Assert (n_blocks() == v.n_blocks(),
          ExcDimensionMismatch(n_blocks(), v.n_blocks()));
  for (size_type i=0; i<n_blocks(); ++i)
    components[i].scale(v.block(i));
}



template <class VectorType>
void BlockVectorBase<VectorType>::equ (const value_type a,
                                       const BlockVectorBase<VectorType> &v,
                                       const value_type b,
                                       const BlockVectorBase<VectorType> &w)
{

  AssertIsFinite(a);
  AssertIsFinite(b);

  Assert (n_blocks() == v.n_blocks(),
          ExcDimensionMismatch(n_blocks(), v.n_blocks()));
  Assert (n_blocks() == w.n_blocks(),
          ExcDimensionMismatch(n_blocks(), w.n_blocks()));

  for (size_type i=0; i<n_blocks(); ++i)
    {
      components[i].equ( a, v.components[i], b, w.components[i]);
    }
}



template <class VectorType>
std::size_t
BlockVectorBase<VectorType>::memory_consumption () const
{
  return (MemoryConsumption::memory_consumption (this->block_indices)
          +
          MemoryConsumption::memory_consumption (this->components));
}



template <class VectorType>
template <class BlockVector2>
void BlockVectorBase<VectorType>::equ (const value_type    a,
                                       const BlockVector2 &v)
{

  AssertIsFinite(a);

  Assert (n_blocks() == v.n_blocks(),
          ExcDimensionMismatch(n_blocks(), v.n_blocks()));

  for (size_type i=0; i<n_blocks(); ++i)
    components[i].equ( a, v.components[i]);
}



template <class VectorType>
void BlockVectorBase<VectorType>::update_ghost_values () const
{
  for (size_type i=0; i<n_blocks(); ++i)
    block(i).update_ghost_values ();
}



template <class VectorType>
BlockVectorBase<VectorType> &
BlockVectorBase<VectorType>::operator = (const value_type s)
{

  AssertIsFinite(s);

  for (size_type i=0; i<n_blocks(); ++i)
    components[i] = s;

  return *this;
}


template <class VectorType>
BlockVectorBase<VectorType> &
BlockVectorBase<VectorType>::operator = (const BlockVectorBase<VectorType> &v)
{
  AssertDimension(n_blocks(), v.n_blocks());

  for (size_type i=0; i<n_blocks(); ++i)
    components[i] = v.components[i];

  return *this;
}


template <class VectorType>
template <class VectorType2>
BlockVectorBase<VectorType> &
BlockVectorBase<VectorType>::operator = (const BlockVectorBase<VectorType2> &v)
{
  AssertDimension(n_blocks(), v.n_blocks());

  for (size_type i=0; i<n_blocks(); ++i)
    components[i] = v.components[i];

  return *this;
}



template <class VectorType>
BlockVectorBase<VectorType> &
BlockVectorBase<VectorType>::operator = (const VectorType &v)
{
  Assert (size() == v.size(),
          ExcDimensionMismatch(size(), v.size()));

  size_type index_v = 0;
  for (size_type b=0; b<n_blocks(); ++b)
    for (size_type i=0; i<block(b).size(); ++i, ++index_v)
      block(b)(i) = v(index_v);

  return *this;
}



template <class VectorType>
template <class VectorType2>
inline
bool
BlockVectorBase<VectorType>::
operator == (const BlockVectorBase<VectorType2> &v) const
{
  Assert (block_indices == v.block_indices, ExcDifferentBlockIndices());

  for (size_type i=0; i<n_blocks(); ++i)
    if ( ! (components[i] == v.components[i]))
      return false;

  return true;
}



template <class VectorType>
inline
BlockVectorBase<VectorType> &
BlockVectorBase<VectorType>::operator *= (const value_type factor)
{

  AssertIsFinite(factor);

  for (size_type i=0; i<n_blocks(); ++i)
    components[i] *= factor;

  return *this;
}



template <class VectorType>
inline
BlockVectorBase<VectorType> &
BlockVectorBase<VectorType>::operator /= (const value_type factor)
{

  AssertIsFinite(factor);
  Assert (factor != 0., ExcDivideByZero() );

  for (size_type i=0; i<n_blocks(); ++i)
    components[i] /= factor;

  return *this;
}


template <class VectorType>
inline
typename BlockVectorBase<VectorType>::value_type
BlockVectorBase<VectorType>::operator() (const size_type i) const
{
  const std::pair<unsigned int,size_type> local_index
    = block_indices.global_to_local (i);
  return components[local_index.first](local_index.second);
}



template <class VectorType>
inline
typename BlockVectorBase<VectorType>::reference
BlockVectorBase<VectorType>::operator() (const size_type i)
{
  const std::pair<unsigned int,size_type> local_index
    = block_indices.global_to_local (i);
  return components[local_index.first](local_index.second);
}



template <class VectorType>
inline
typename BlockVectorBase<VectorType>::value_type
BlockVectorBase<VectorType>::operator[] (const size_type i) const
{
  return operator()(i);
}



template <class VectorType>
inline
typename BlockVectorBase<VectorType>::reference
BlockVectorBase<VectorType>::operator[] (const size_type i)
{
  return operator()(i);
}



template <typename VectorType>
template <typename OtherNumber>
inline
void BlockVectorBase<VectorType>::extract_subvector_to (const std::vector<size_type> &indices,
                                                        std::vector<OtherNumber> &values) const
{
  for (size_type i = 0; i < indices.size(); ++i)
    values[i] = operator()(indices[i]);
}



template <typename VectorType>
template <typename ForwardIterator, typename OutputIterator>
inline
void BlockVectorBase<VectorType>::extract_subvector_to (ForwardIterator          indices_begin,
                                                        const ForwardIterator    indices_end,
                                                        OutputIterator           values_begin) const
{
  while (indices_begin != indices_end)
    {
      *values_begin = operator()(*indices_begin);
      indices_begin++;
      values_begin++;
    }
}

#endif // DOXYGEN

DEAL_II_NAMESPACE_CLOSE

#endif