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.md at
// the top level directory of deal.II.
//
// ---------------------------------------------------------------------

#ifndef dealii_block_vector_base_h
#define dealii_block_vector_base_h


#include <deal.II/base/config.h>

#include <deal.II/base/exceptions.h>
#include <deal.II/base/memory_consumption.h>
#include <deal.II/base/numbers.h>
#include <deal.II/base/subscriptor.h>

#include <deal.II/lac/block_indices.h>
#include <deal.II/lac/exceptions.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/vector_operation.h>

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

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>
    {
      using Vector = typename BlockVectorType::BlockType;

      using BlockVector = BlockVectorType;

      using value_type = typename BlockVector::value_type;

      using dereference_type = typename Vector::reference;
    };



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

      using BlockVector = const BlockVectorType;

      using value_type = const typename BlockVector::value_type;

      using dereference_type = value_type;
    };


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

      using value_type = typename Types<BlockVectorType, Constness>::value_type;

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

      using dereference_type =
        typename Types<BlockVectorType, Constness>::dereference_type;

      using BlockVector =
        typename Types<BlockVectorType, Constness>::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:
  using BlockType = VectorType;

  /*
   * 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.
   */
  using value_type    = typename BlockType::value_type;
  using pointer       = value_type *;
  using const_pointer = const value_type *;
  using iterator =
    ::internal::BlockVectorIterators::Iterator<BlockVectorBase, false>;
  using const_iterator =
    ::internal::BlockVectorIterators::Iterator<BlockVectorBase, true>;
  using reference       = typename BlockType::reference;
  using const_reference = typename BlockType::const_reference;
  using size_type       = types::global_dof_index;

  using real_type = typename BlockType::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