block_matrix_base.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2004 - 2017 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_matrix_base_h
#define dealii_block_matrix_base_h


#include <deal.II/base/config.h>
#include <deal.II/base/table.h>
#include <deal.II/base/thread_management.h>
#include <deal.II/base/utilities.h>
#include <deal.II/base/smartpointer.h>
#include <deal.II/base/memory_consumption.h>
#include <deal.II/lac/block_indices.h>
#include <deal.II/lac/exceptions.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/matrix_iterator.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/vector_operation.h>

#include <cmath>

DEAL_II_NAMESPACE_OPEN


template <typename> class MatrixIterator;



namespace BlockMatrixIterators
{
  template <class BlockMatrixType>
  class AccessorBase
  {
  public:
    typedef types::global_dof_index size_type;

    typedef typename BlockMatrixType::value_type value_type;

    AccessorBase ();

    unsigned int block_row() const;

    unsigned int block_column() const;

  protected:
    unsigned int row_block;

    unsigned int col_block;

    template <typename>
    friend class MatrixIterator;
  };



  template <class BlockMatrixType, bool Constness>
  class Accessor;


  template <class BlockMatrixType>
  class Accessor<BlockMatrixType, false>
    :
    public AccessorBase<BlockMatrixType>
  {
  public:
    typedef types::global_dof_index size_type;

    typedef BlockMatrixType MatrixType;

    typedef typename BlockMatrixType::value_type value_type;

    Accessor (BlockMatrixType *m,
              const size_type row,
              const size_type col);

    size_type row() const;

    size_type column() const;

    value_type value() const;

    void set_value(value_type newval) const;

  protected:
    BlockMatrixType *matrix;

    typename BlockMatrixType::BlockType::iterator base_iterator;

    void advance ();

    bool operator == (const Accessor &a) const;

    template <typename> friend class MatrixIterator;
    friend class Accessor<BlockMatrixType, true>;
  };

  template <class BlockMatrixType>
  class Accessor<BlockMatrixType, true>
    :
    public AccessorBase<BlockMatrixType>
  {
  public:
    typedef types::global_dof_index size_type;

    typedef const BlockMatrixType MatrixType;

    typedef typename BlockMatrixType::value_type value_type;

    Accessor (const BlockMatrixType *m,
              const size_type row,
              const size_type col);

    Accessor(const Accessor<BlockMatrixType, false> &);

    size_type row() const;

    size_type column() const;

    value_type value() const;
  protected:
    const BlockMatrixType *matrix;

    typename BlockMatrixType::BlockType::const_iterator base_iterator;

    void advance ();

    bool operator == (const Accessor &a) const;

    template <typename>
    friend class ::MatrixIterator;
  };
}



template <typename MatrixType>
class BlockMatrixBase : public Subscriptor
{
public:
  typedef MatrixType BlockType;

  typedef typename BlockType::value_type value_type;
  typedef value_type             *pointer;
  typedef const value_type       *const_pointer;
  typedef value_type             &reference;
  typedef const value_type       &const_reference;
  typedef types::global_dof_index size_type;

  typedef
  MatrixIterator<BlockMatrixIterators::Accessor<BlockMatrixBase, false> >
  iterator;

  typedef
  MatrixIterator<BlockMatrixIterators::Accessor<BlockMatrixBase, true> >
  const_iterator;


  BlockMatrixBase () = default;

  ~BlockMatrixBase ();

  template <class BlockMatrixType>
  BlockMatrixBase &
  copy_from (const BlockMatrixType &source);

  BlockType &
  block (const unsigned int row,
         const unsigned int column);


  const BlockType &
  block (const unsigned int row,
         const unsigned int column) const;

  size_type m () const;

  size_type n () const;


  unsigned int n_block_rows () const;

  unsigned int n_block_cols () const;

  void set (const size_type i,
            const size_type j,
            const value_type value);

  template <typename number>
  void set (const std::vector<size_type> &indices,
            const FullMatrix<number>     &full_matrix,
            const bool                    elide_zero_values = false);

  template <typename number>
  void set (const std::vector<size_type> &row_indices,
            const std::vector<size_type> &col_indices,
            const FullMatrix<number>     &full_matrix,
            const bool                    elide_zero_values = false);

  template <typename number>
  void set (const size_type row,
            const std::vector<size_type> &col_indices,
            const std::vector<number>    &values,
            const bool                    elide_zero_values = false);

  template <typename number>
  void set (const size_type  row,
            const size_type  n_cols,
            const size_type *col_indices,
            const number    *values,
            const bool       elide_zero_values = false);

  void add (const size_type i,
            const size_type j,
            const value_type value);

  template <typename number>
  void add (const std::vector<size_type> &indices,
            const FullMatrix<number>     &full_matrix,
            const bool                    elide_zero_values = true);

  template <typename number>
  void add (const std::vector<size_type> &row_indices,
            const std::vector<size_type> &col_indices,
            const FullMatrix<number>     &full_matrix,
            const bool                    elide_zero_values = true);

  template <typename number>
  void add (const size_type row,
            const std::vector<size_type> &col_indices,
            const std::vector<number>    &values,
            const bool                    elide_zero_values = true);

  template <typename number>
  void add (const size_type  row,
            const size_type  n_cols,
            const size_type *col_indices,
            const number    *values,
            const bool       elide_zero_values = true,
            const bool       col_indices_are_sorted = false);

  void add (const value_type                   factor,
            const BlockMatrixBase<MatrixType> &matrix);

  value_type operator () (const size_type i,
                          const size_type j) const;

  value_type el (const size_type i,
                 const size_type j) const;

  value_type diag_element (const size_type i) const;

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

  BlockMatrixBase &operator *= (const value_type factor);

  BlockMatrixBase &operator /= (const value_type factor);

  template <class BlockVectorType>
  void vmult_add (BlockVectorType       &dst,
                  const BlockVectorType &src) const;

  template <class BlockVectorType>
  void Tvmult_add (BlockVectorType       &dst,
                   const BlockVectorType &src) const;

  template <class BlockVectorType>
  value_type
  matrix_norm_square (const BlockVectorType &v) const;

  template <class BlockVectorType>
  value_type
  matrix_scalar_product (const BlockVectorType &u,
                         const BlockVectorType &v) const;

  template <class BlockVectorType>
  value_type residual (BlockVectorType       &dst,
                       const BlockVectorType &x,
                       const BlockVectorType &b) const;

  void print (std::ostream &out,
              const bool    alternative_output = false) const;

  iterator begin ();

  iterator end ();

  iterator begin (const size_type r);

  iterator end (const size_type r);
  const_iterator begin () const;

  const_iterator end () const;

  const_iterator begin (const size_type r) const;

  const_iterator end (const size_type r) const;

  const BlockIndices &get_row_indices () const;

  const BlockIndices &get_column_indices () const;

  std::size_t memory_consumption () const;

  DeclException4 (ExcIncompatibleRowNumbers,
                  int, int, int, int,
                  << "The blocks [" << arg1 << ',' << arg2 << "] and ["
                  << arg3 << ',' << arg4 << "] have differing row numbers.");
  DeclException4 (ExcIncompatibleColNumbers,
                  int, int, int, int,
                  << "The blocks [" << arg1 << ',' << arg2 << "] and ["
                  << arg3 << ',' << arg4 << "] have differing column numbers.");
protected:
  void clear ();

  BlockIndices row_block_indices;
  BlockIndices column_block_indices;

  Table<2,SmartPointer<BlockType, BlockMatrixBase<MatrixType> > > sub_objects;

  void collect_sizes ();

  template <class BlockVectorType>
  void vmult_block_block (BlockVectorType       &dst,
                          const BlockVectorType &src) const;

  template <class BlockVectorType,
            class VectorType>
  void vmult_block_nonblock (BlockVectorType          &dst,
                             const VectorType &src) const;

  template <class BlockVectorType,
            class VectorType>
  void vmult_nonblock_block (VectorType    &dst,
                             const BlockVectorType &src) const;

  template <class VectorType>
  void vmult_nonblock_nonblock (VectorType       &dst,
                                const VectorType &src) const;

  template <class BlockVectorType>
  void Tvmult_block_block (BlockVectorType       &dst,
                           const BlockVectorType &src) const;

  template <class BlockVectorType,
            class VectorType>
  void Tvmult_block_nonblock (BlockVectorType  &dst,
                              const VectorType &src) const;

  template <class BlockVectorType,
            class VectorType>
  void Tvmult_nonblock_block (VectorType    &dst,
                              const BlockVectorType &src) const;

  template <class VectorType>
  void Tvmult_nonblock_nonblock (VectorType       &dst,
                                 const VectorType &src) const;


protected:

  void prepare_add_operation();

  void prepare_set_operation();


private:

  struct TemporaryData
  {
    std::vector<size_type> counter_within_block;

    std::vector<std::vector<size_type> > column_indices;

    std::vector<std::vector<value_type> > column_values;

    Threads::Mutex mutex;

    TemporaryData &operator = (const TemporaryData &)
    {
      return *this;
    }
  };

  TemporaryData temporary_data;

  template <typename, bool>
  friend class BlockMatrixIterators::Accessor;

  template <typename>
  friend class MatrixIterator;
};


#ifndef DOXYGEN
/* ------------------------- Template functions ---------------------- */


namespace BlockMatrixIterators
{
  template <class BlockMatrixType>
  inline
  AccessorBase<BlockMatrixType>::AccessorBase()
    :
    row_block(0),
    col_block(0)
  {}


  template <class BlockMatrixType>
  inline
  unsigned int
  AccessorBase<BlockMatrixType>::block_row() const
  {
    Assert (row_block != numbers::invalid_unsigned_int,
            ExcIteratorPastEnd());

    return row_block;
  }


  template <class BlockMatrixType>
  inline
  unsigned int
  AccessorBase<BlockMatrixType>::block_column() const
  {
    Assert (col_block != numbers::invalid_unsigned_int,
            ExcIteratorPastEnd());

    return col_block;
  }


  template <class BlockMatrixType>
  inline
  Accessor<BlockMatrixType, true>::Accessor (
    const BlockMatrixType  *matrix,
    const size_type        row,
    const size_type        col)
    :
    matrix(matrix),
    base_iterator(matrix->block(0,0).begin())
  {
    (void)col;
    Assert(col==0, ExcNotImplemented());

    // check if this is a regular row or
    // the end of the matrix
    if (row < matrix->m())
      {
        const std::pair<unsigned int,size_type> indices
          = matrix->row_block_indices.global_to_local(row);

        // find the first block that does
        // have an entry in this row
        for (unsigned int bc=0; bc<matrix->n_block_cols(); ++bc)
          {
            base_iterator
              = matrix->block(indices.first, bc).begin(indices.second);
            if (base_iterator !=
                matrix->block(indices.first, bc).end(indices.second))
              {
                this->row_block = indices.first;
                this->col_block = bc;
                return;
              }
          }

        // hm, there is no block that has
        // an entry in this column. we need
        // to take the next entry then,
        // which may be the first entry of
        // the next row, or recursively the
        // next row, or so on
        *this = Accessor (matrix, row+1, 0);
      }
    else
      {
        // we were asked to create the end
        // iterator for this matrix
        this->row_block = numbers::invalid_unsigned_int;
        this->col_block = numbers::invalid_unsigned_int;
      }
  }


//   template <class BlockMatrixType>
//   inline
//   Accessor<BlockMatrixType, true>::Accessor (const Accessor<BlockMatrixType, true>& other)
//                :
//                matrix(other.matrix),
//                base_iterator(other.base_iterator)
//   {
//     this->row_block = other.row_block;
//     this->col_block = other.col_block;
//   }


  template <class BlockMatrixType>
  inline
  Accessor<BlockMatrixType, true>::Accessor (const Accessor<BlockMatrixType, false> &other)
    :
    matrix(other.matrix),
    base_iterator(other.base_iterator)
  {
    this->row_block = other.row_block;
    this->col_block = other.col_block;
  }


  template <class BlockMatrixType>
  inline
  typename Accessor<BlockMatrixType, true>::size_type
  Accessor<BlockMatrixType, true>::row() const
  {
    Assert (this->row_block != numbers::invalid_unsigned_int,
            ExcIteratorPastEnd());

    return (matrix->row_block_indices.local_to_global(this->row_block, 0) +
            base_iterator->row());
  }


  template <class BlockMatrixType>
  inline
  typename Accessor<BlockMatrixType, true>::size_type
  Accessor<BlockMatrixType, true>::column() const
  {
    Assert (this->col_block != numbers::invalid_unsigned_int,
            ExcIteratorPastEnd());

    return (matrix->column_block_indices.local_to_global(this->col_block,0) +
            base_iterator->column());
  }


  template <class BlockMatrixType>
  inline
  typename Accessor<BlockMatrixType, true>::value_type
  Accessor<BlockMatrixType, true>::value () const
  {
    Assert (this->row_block != numbers::invalid_unsigned_int,
            ExcIteratorPastEnd());
    Assert (this->col_block != numbers::invalid_unsigned_int,
            ExcIteratorPastEnd());

    return base_iterator->value();
  }



  template <class BlockMatrixType>
  inline
  void
  Accessor<BlockMatrixType, true>::advance ()
  {
    Assert (this->row_block != numbers::invalid_unsigned_int,
            ExcIteratorPastEnd());
    Assert (this->col_block != numbers::invalid_unsigned_int,
            ExcIteratorPastEnd());

    // Remember current row inside block
    size_type local_row = base_iterator->row();

    // Advance one element inside the
    // current block
    ++base_iterator;

    // while we hit the end of the row of a
    // block (which may happen multiple
    // times if rows inside a block are
    // empty), we have to jump to the next
    // block and take the
    while (base_iterator ==
           matrix->block(this->row_block, this->col_block).end(local_row))
      {
        // jump to next block in this block
        // row, if possible, otherwise go
        // to next row
        if (this->col_block < matrix->n_block_cols()-1)
          {
            ++this->col_block;
            base_iterator
              = matrix->block(this->row_block, this->col_block).begin(local_row);
          }
        else
          {
            // jump back to next row in
            // first block column
            this->col_block = 0;
            ++local_row;

            // see if this has brought us
            // past the number of rows in
            // this block. if so see
            // whether we've just fallen
            // off the end of the whole
            // matrix
            if (local_row == matrix->block(this->row_block, this->col_block).m())
              {
                local_row = 0;
                ++this->row_block;
                if (this->row_block == matrix->n_block_rows())
                  {
                    this->row_block = numbers::invalid_unsigned_int;
                    this->col_block = numbers::invalid_unsigned_int;
                    return;
                  }
              }

            base_iterator
              = matrix->block(this->row_block, this->col_block).begin(local_row);
          }
      }
  }


  template <class BlockMatrixType>
  inline
  bool
  Accessor<BlockMatrixType, true>::operator == (const Accessor &a) const
  {
    if (matrix != a.matrix)
      return false;

    if (this->row_block == a.row_block
        && this->col_block == a.col_block)
      // end iterators do not necessarily
      // have to have the same
      // base_iterator representation, but
      // valid iterators have to
      return (((this->row_block == numbers::invalid_unsigned_int)
               &&
               (this->col_block == numbers::invalid_unsigned_int))
              ||
              (base_iterator == a.base_iterator));

    return false;
  }

//----------------------------------------------------------------------//


  template <class BlockMatrixType>
  inline
  Accessor<BlockMatrixType, false>::Accessor (
    BlockMatrixType  *matrix,
    const size_type  row,
    const size_type  col)
    :
    matrix(matrix),
    base_iterator(matrix->block(0,0).begin())
  {
    (void)col;
    Assert(col==0, ExcNotImplemented());
    // check if this is a regular row or
    // the end of the matrix
    if (row < matrix->m())
      {
        const std::pair<unsigned int,size_type> indices
          = matrix->row_block_indices.global_to_local(row);

        // find the first block that does
        // have an entry in this row
        for (size_type bc=0; bc<matrix->n_block_cols(); ++bc)
          {
            base_iterator
              = matrix->block(indices.first, bc).begin(indices.second);
            if (base_iterator !=
                matrix->block(indices.first, bc).end(indices.second))
              {
                this->row_block = indices.first;
                this->col_block = bc;
                return;
              }
          }

        // hm, there is no block that has
        // an entry in this column. we need
        // to take the next entry then,
        // which may be the first entry of
        // the next row, or recursively the
        // next row, or so on
        *this = Accessor (matrix, row+1, 0);
      }
    else
      {
        // we were asked to create the end
        // iterator for this matrix
        this->row_block = numbers::invalid_size_type;
        this->col_block = numbers::invalid_size_type;
      }
  }


  template <class BlockMatrixType>
  inline
  typename Accessor<BlockMatrixType, false>::size_type
  Accessor<BlockMatrixType, false>::row() const
  {
    Assert (this->row_block != numbers::invalid_size_type,
            ExcIteratorPastEnd());

    return (matrix->row_block_indices.local_to_global(this->row_block, 0) +
            base_iterator->row());
  }


  template <class BlockMatrixType>
  inline
  typename Accessor<BlockMatrixType, false>::size_type
  Accessor<BlockMatrixType, false>::column() const
  {
    Assert (this->col_block != numbers::invalid_size_type,
            ExcIteratorPastEnd());

    return (matrix->column_block_indices.local_to_global(this->col_block,0) +
            base_iterator->column());
  }


  template <class BlockMatrixType>
  inline
  typename Accessor<BlockMatrixType, false>::value_type
  Accessor<BlockMatrixType, false>::value () const
  {
    Assert (this->row_block != numbers::invalid_size_type,
            ExcIteratorPastEnd());
    Assert (this->col_block != numbers::invalid_size_type,
            ExcIteratorPastEnd());

    return base_iterator->value();
  }



  template <class BlockMatrixType>
  inline
  void
  Accessor<BlockMatrixType, false>::set_value (typename Accessor<BlockMatrixType, false>::value_type newval) const
  {
    Assert (this->row_block != numbers::invalid_size_type,
            ExcIteratorPastEnd());
    Assert (this->col_block != numbers::invalid_size_type,
            ExcIteratorPastEnd());

    base_iterator->value() = newval;
  }



  template <class BlockMatrixType>
  inline
  void
  Accessor<BlockMatrixType, false>::advance ()
  {
    Assert (this->row_block != numbers::invalid_size_type,
            ExcIteratorPastEnd());
    Assert (this->col_block != numbers::invalid_size_type,
            ExcIteratorPastEnd());

    // Remember current row inside block
    size_type local_row = base_iterator->row();

    // Advance one element inside the
    // current block
    ++base_iterator;

    // while we hit the end of the row of a
    // block (which may happen multiple
    // times if rows inside a block are
    // empty), we have to jump to the next
    // block and take the
    while (base_iterator ==
           matrix->block(this->row_block, this->col_block).end(local_row))
      {
        // jump to next block in this block
        // row, if possible, otherwise go
        // to next row
        if (this->col_block < matrix->n_block_cols()-1)
          {
            ++this->col_block;
            base_iterator
              = matrix->block(this->row_block, this->col_block).begin(local_row);
          }
        else
          {
            // jump back to next row in
            // first block column
            this->col_block = 0;
            ++local_row;

            // see if this has brought us
            // past the number of rows in
            // this block. if so see
            // whether we've just fallen
            // off the end of the whole
            // matrix
            if (local_row == matrix->block(this->row_block, this->col_block).m())
              {
                local_row = 0;
                ++this->row_block;
                if (this->row_block == matrix->n_block_rows())
                  {
                    this->row_block = numbers::invalid_size_type;
                    this->col_block = numbers::invalid_size_type;
                    return;
                  }
              }

            base_iterator
              = matrix->block(this->row_block, this->col_block).begin(local_row);
          }
      }
  }



  template <class BlockMatrixType>
  inline
  bool
  Accessor<BlockMatrixType, false>::operator == (const Accessor &a) const
  {
    if (matrix != a.matrix)
      return false;

    if (this->row_block == a.row_block
        && this->col_block == a.col_block)
      // end iterators do not necessarily
      // have to have the same
      // base_iterator representation, but
      // valid iterators have to
      return (((this->row_block == numbers::invalid_size_type)
               &&
               (this->col_block == numbers::invalid_size_type))
              ||
              (base_iterator == a.base_iterator));

    return false;
  }
}


//---------------------------------------------------------------------------

template <typename MatrixType>
inline
BlockMatrixBase<MatrixType>::~BlockMatrixBase ()
{
  try
    {
      clear ();
    }
  catch (...)
    {}
}


template <class MatrixType>
template <class BlockMatrixType>
inline
BlockMatrixBase<MatrixType> &
BlockMatrixBase<MatrixType>::
copy_from (const BlockMatrixType &source)
{
  for (unsigned int r=0; r<n_block_rows(); ++r)
    for (unsigned int c=0; c<n_block_cols(); ++c)
      block(r,c).copy_from (source.block(r,c));

  return *this;
}


template <class MatrixType>
std::size_t
BlockMatrixBase<MatrixType>::memory_consumption () const
{
  std::size_t mem =
    MemoryConsumption::memory_consumption(row_block_indices)+
    MemoryConsumption::memory_consumption(column_block_indices)+
    MemoryConsumption::memory_consumption(sub_objects)+
    MemoryConsumption::memory_consumption(temporary_data.counter_within_block)+
    MemoryConsumption::memory_consumption(temporary_data.column_indices)+
    MemoryConsumption::memory_consumption(temporary_data.column_values)+
    sizeof(temporary_data.mutex);

  for (unsigned int r=0; r<n_block_rows(); ++r)
    for (unsigned int c=0; c<n_block_cols(); ++c)
      {
        MatrixType *p = this->sub_objects[r][c];
        mem += MemoryConsumption::memory_consumption(*p);
      }

  return mem;
}



template <class MatrixType>
inline
void
BlockMatrixBase<MatrixType>::clear ()
{
  for (unsigned int r=0; r<n_block_rows(); ++r)
    for (unsigned int c=0; c<n_block_cols(); ++c)
      {
        MatrixType *p = this->sub_objects[r][c];
        this->sub_objects[r][c] = nullptr;
        delete p;
      }
  sub_objects.reinit (0,0);

  // reset block indices to empty
  row_block_indices = column_block_indices = BlockIndices ();
}



template <class MatrixType>
inline
typename BlockMatrixBase<MatrixType>::BlockType &
BlockMatrixBase<MatrixType>::block (const unsigned int row,
                                    const unsigned int column)
{
  Assert (row<n_block_rows(),
          ExcIndexRange (row, 0, n_block_rows()));
  Assert (column<n_block_cols(),
          ExcIndexRange (column, 0, n_block_cols()));

  return *sub_objects[row][column];
}



template <class MatrixType>
inline
const typename BlockMatrixBase<MatrixType>::BlockType &
BlockMatrixBase<MatrixType>::block (const unsigned int row,
                                    const unsigned int column) const
{
  Assert (row<n_block_rows(),
          ExcIndexRange (row, 0, n_block_rows()));
  Assert (column<n_block_cols(),
          ExcIndexRange (column, 0, n_block_cols()));

  return *sub_objects[row][column];
}


template <class MatrixType>
inline
typename BlockMatrixBase<MatrixType>::size_type
BlockMatrixBase<MatrixType>::m () const
{
  return row_block_indices.total_size();
}



template <class MatrixType>
inline
typename BlockMatrixBase<MatrixType>::size_type
BlockMatrixBase<MatrixType>::n () const
{
  return column_block_indices.total_size();
}



template <class MatrixType>
inline
unsigned int
BlockMatrixBase<MatrixType>::n_block_cols () const
{
  return column_block_indices.size();
}



template <class MatrixType>
inline
unsigned int
BlockMatrixBase<MatrixType>::n_block_rows () const
{
  return row_block_indices.size();
}



// Write the single set manually,
// since the other function has a lot
// of overhead in that case.
template <class MatrixType>
inline
void
BlockMatrixBase<MatrixType>::set (const size_type i,
                                  const size_type j,
                                  const value_type value)
{
  prepare_set_operation();

  AssertIsFinite(value);

  const std::pair<unsigned int,size_type>
  row_index = row_block_indices.global_to_local (i),
  col_index = column_block_indices.global_to_local (j);
  block(row_index.first,col_index.first).set (row_index.second,
                                              col_index.second,
                                              value);
}



template <class MatrixType>
template <typename number>
inline
void
BlockMatrixBase<MatrixType>::set (const std::vector<size_type> &row_indices,
                                  const std::vector<size_type> &col_indices,
                                  const FullMatrix<number>     &values,
                                  const bool                    elide_zero_values)
{
  Assert (row_indices.size() == values.m(),
          ExcDimensionMismatch(row_indices.size(), values.m()));
  Assert (col_indices.size() == values.n(),
          ExcDimensionMismatch(col_indices.size(), values.n()));

  for (size_type i=0; i<row_indices.size(); ++i)
    set (row_indices[i], col_indices.size(), col_indices.data(), &values(i,0),
         elide_zero_values);
}



template <class MatrixType>
template <typename number>
inline
void
BlockMatrixBase<MatrixType>::set (const std::vector<size_type> &indices,
                                  const FullMatrix<number>     &values,
                                  const bool                    elide_zero_values)
{
  Assert (indices.size() == values.m(),
          ExcDimensionMismatch(indices.size(), values.m()));
  Assert (values.n() == values.m(), ExcNotQuadratic());

  for (size_type i=0; i<indices.size(); ++i)
    set (indices[i], indices.size(), indices.data(), &values(i,0),
         elide_zero_values);
}



template <class MatrixType>
template <typename number>
inline
void
BlockMatrixBase<MatrixType>::set (const size_type               row,
                                  const std::vector<size_type> &col_indices,
                                  const std::vector<number>    &values,
                                  const bool                    elide_zero_values)
{
  Assert (col_indices.size() == values.size(),
          ExcDimensionMismatch(col_indices.size(), values.size()));

  set (row, col_indices.size(), col_indices.data(), values.data(),
       elide_zero_values);
}



// This is a very messy function, since
// we need to calculate to each position
// the location in the global array.
template <class MatrixType>
template <typename number>
inline
void
BlockMatrixBase<MatrixType>::set (const size_type  row,
                                  const size_type  n_cols,
                                  const size_type *col_indices,
                                  const number    *values,
                                  const bool       elide_zero_values)
{
  prepare_set_operation();

  // lock access to the temporary data structure to
  // allow multiple threads to call this function concurrently
  Threads::Mutex::ScopedLock lock (temporary_data.mutex);

  // Resize scratch arrays
  if (temporary_data.column_indices.size() < this->n_block_cols())
    {
      temporary_data.column_indices.resize (this->n_block_cols());
      temporary_data.column_values.resize (this->n_block_cols());
      temporary_data.counter_within_block.resize (this->n_block_cols());
    }

  // Resize sub-arrays to n_cols. This
  // is a bit wasteful, but we resize
  // only a few times (then the maximum
  // row length won't increase that
  // much any more). At least we know
  // that all arrays are going to be of
  // the same size, so we can check
  // whether the size of one is large
  // enough before actually going
  // through all of them.
  if (temporary_data.column_indices[0].size() < n_cols)
    {
      for (unsigned int i=0; i<this->n_block_cols(); ++i)
        {
          temporary_data.column_indices[i].resize(n_cols);
          temporary_data.column_values[i].resize(n_cols);
        }
    }

  // Reset the number of added elements
  // in each block to zero.
  for (unsigned int i=0; i<this->n_block_cols(); ++i)
    temporary_data.counter_within_block[i] = 0;

  // Go through the column indices to
  // find out which portions of the
  // values should be set in which
  // block of the matrix. We need to
  // touch all the data, since we can't
  // be sure that the data of one block
  // is stored contiguously (in fact,
  // indices will be intermixed when it
  // comes from an element matrix).
  for (size_type j=0; j<n_cols; ++j)
    {
      number value = values[j];

      if (value == number() && elide_zero_values == true)
        continue;

      const std::pair<unsigned int, size_type>
      col_index = this->column_block_indices.global_to_local(col_indices[j]);

      const size_type local_index = temporary_data.counter_within_block[col_index.first]++;

      temporary_data.column_indices[col_index.first][local_index] = col_index.second;
      temporary_data.column_values[col_index.first][local_index] = value;
    }

#ifdef DEBUG
  // If in debug mode, do a check whether
  // the right length has been obtained.
  size_type length = 0;
  for (unsigned int i=0; i<this->n_block_cols(); ++i)
    length += temporary_data.counter_within_block[i];
  Assert (length <= n_cols, ExcInternalError());
#endif

  // Now we found out about where the
  // individual columns should start and
  // where we should start reading out
  // data. Now let's write the data into
  // the individual blocks!
  const std::pair<unsigned int,size_type>
  row_index = this->row_block_indices.global_to_local (row);
  for (unsigned int block_col=0; block_col<n_block_cols(); ++block_col)
    {
      if (temporary_data.counter_within_block[block_col] == 0)
        continue;

      block(row_index.first, block_col).set
      (row_index.second,
       temporary_data.counter_within_block[block_col],
       &temporary_data.column_indices[block_col][0],
       &temporary_data.column_values[block_col][0],
       false);
    }
}



template <class MatrixType>
inline
void
BlockMatrixBase<MatrixType>::add (const size_type  i,
                                  const size_type  j,
                                  const value_type value)
{

  AssertIsFinite(value);

  prepare_add_operation();

  // save some cycles for zero additions, but
  // only if it is safe for the matrix we are
  // working with
  typedef typename MatrixType::Traits MatrixTraits;
  if ((MatrixTraits::zero_addition_can_be_elided == true)
      &&
      (value == value_type()))
    return;

  const std::pair<unsigned int,size_type>
  row_index = row_block_indices.global_to_local (i),
  col_index = column_block_indices.global_to_local (j);
  block(row_index.first,col_index.first).add (row_index.second,
                                              col_index.second,
                                              value);
}



template <class MatrixType>
template <typename number>
inline
void
BlockMatrixBase<MatrixType>::add (const std::vector<size_type> &row_indices,
                                  const std::vector<size_type> &col_indices,
                                  const FullMatrix<number>     &values,
                                  const bool                    elide_zero_values)
{
  Assert (row_indices.size() == values.m(),
          ExcDimensionMismatch(row_indices.size(), values.m()));
  Assert (col_indices.size() == values.n(),
          ExcDimensionMismatch(col_indices.size(), values.n()));

  for (size_type i=0; i<row_indices.size(); ++i)
    add (row_indices[i], col_indices.size(), col_indices.data(), &values(i,0),
         elide_zero_values);
}



template <class MatrixType>
template <typename number>
inline
void
BlockMatrixBase<MatrixType>::add (const std::vector<size_type> &indices,
                                  const FullMatrix<number>     &values,
                                  const bool                    elide_zero_values)
{
  Assert (indices.size() == values.m(),
          ExcDimensionMismatch(indices.size(), values.m()));
  Assert (values.n() == values.m(), ExcNotQuadratic());

  for (size_type i=0; i<indices.size(); ++i)
    add (indices[i], indices.size(), indices.data(), &values(i,0),
         elide_zero_values);
}



template <class MatrixType>
template <typename number>
inline
void
BlockMatrixBase<MatrixType>::add (const size_type               row,
                                  const std::vector<size_type> &col_indices,
                                  const std::vector<number>    &values,
                                  const bool                    elide_zero_values)
{
  Assert (col_indices.size() == values.size(),
          ExcDimensionMismatch(col_indices.size(), values.size()));

  add (row, col_indices.size(), col_indices.data(), values.data(),
       elide_zero_values);
}



// This is a very messy function, since
// we need to calculate to each position
// the location in the global array.
template <class MatrixType>
template <typename number>
inline
void
BlockMatrixBase<MatrixType>::add (const size_type  row,
                                  const size_type  n_cols,
                                  const size_type *col_indices,
                                  const number    *values,
                                  const bool       elide_zero_values,
                                  const bool       col_indices_are_sorted)
{
  prepare_add_operation();

  // TODO: Look over this to find out
  // whether we can do that more
  // efficiently.
  if (col_indices_are_sorted == true)
    {
#ifdef DEBUG
      // check whether indices really are
      // sorted.
      size_type before = col_indices[0];
      for (size_type i=1; i<n_cols; ++i)
        if (col_indices[i] <= before)
          Assert (false, ExcMessage ("Flag col_indices_are_sorted is set, but "
                                     "indices appear to not be sorted."))
          else
            before = col_indices[i];
#endif
      const std::pair<unsigned int,size_type>
      row_index = this->row_block_indices.global_to_local (row);

      if (this->n_block_cols() > 1)
        {
          const size_type *first_block = Utilities::lower_bound (col_indices,
                                                                 col_indices+n_cols,
                                                                 this->column_block_indices.block_start(1));

          const size_type n_zero_block_indices = first_block - col_indices;
          block(row_index.first, 0).add (row_index.second,
                                         n_zero_block_indices,
                                         col_indices,
                                         values,
                                         elide_zero_values,
                                         col_indices_are_sorted);

          if (n_zero_block_indices < n_cols)
            this->add(row, n_cols - n_zero_block_indices, first_block,
                      values + n_zero_block_indices, elide_zero_values,
                      false);
        }
      else
        {
          block(row_index.first, 0). add (row_index.second,
                                          n_cols,
                                          col_indices,
                                          values,
                                          elide_zero_values,
                                          col_indices_are_sorted);
        }

      return;
    }

  // Lock scratch arrays, then resize them
  Threads::Mutex::ScopedLock lock (temporary_data.mutex);

  if (temporary_data.column_indices.size() < this->n_block_cols())
    {
      temporary_data.column_indices.resize (this->n_block_cols());
      temporary_data.column_values.resize (this->n_block_cols());
      temporary_data.counter_within_block.resize (this->n_block_cols());
    }

  // Resize sub-arrays to n_cols. This
  // is a bit wasteful, but we resize
  // only a few times (then the maximum
  // row length won't increase that
  // much any more). At least we know
  // that all arrays are going to be of
  // the same size, so we can check
  // whether the size of one is large
  // enough before actually going
  // through all of them.
  if (temporary_data.column_indices[0].size() < n_cols)
    {
      for (unsigned int i=0; i<this->n_block_cols(); ++i)
        {
          temporary_data.column_indices[i].resize(n_cols);
          temporary_data.column_values[i].resize(n_cols);
        }
    }

  // Reset the number of added elements
  // in each block to zero.
  for (unsigned int i=0; i<this->n_block_cols(); ++i)
    temporary_data.counter_within_block[i] = 0;

  // Go through the column indices to
  // find out which portions of the
  // values should be written into
  // which block of the matrix. We need
  // to touch all the data, since we
  // can't be sure that the data of one
  // block is stored contiguously (in
  // fact, data will be intermixed when
  // it comes from an element matrix).
  for (size_type j=0; j<n_cols; ++j)
    {
      number value = values[j];

      if (value == number() && elide_zero_values == true)
        continue;

      const std::pair<unsigned int, size_type>
      col_index = this->column_block_indices.global_to_local(col_indices[j]);

      const size_type local_index = temporary_data.counter_within_block[col_index.first]++;

      temporary_data.column_indices[col_index.first][local_index] = col_index.second;
      temporary_data.column_values[col_index.first][local_index] = value;
    }

#ifdef DEBUG
  // If in debug mode, do a check whether
  // the right length has been obtained.
  size_type length = 0;
  for (unsigned int i=0; i<this->n_block_cols(); ++i)
    length += temporary_data.counter_within_block[i];
  Assert (length <= n_cols, ExcInternalError());
#endif

  // Now we found out about where the
  // individual columns should start and
  // where we should start reading out
  // data. Now let's write the data into
  // the individual blocks!
  const std::pair<unsigned int,size_type>
  row_index = this->row_block_indices.global_to_local (row);
  for (unsigned int block_col=0; block_col<n_block_cols(); ++block_col)
    {
      if (temporary_data.counter_within_block[block_col] == 0)
        continue;

      block(row_index.first, block_col).add
      (row_index.second,
       temporary_data.counter_within_block[block_col],
       &temporary_data.column_indices[block_col][0],
       &temporary_data.column_values[block_col][0],
       false,
       col_indices_are_sorted);
    }
}



template <class MatrixType>
inline
void
BlockMatrixBase<MatrixType>::add (const value_type                   factor,
                                  const BlockMatrixBase<MatrixType> &matrix)
{
  AssertIsFinite(factor);

  prepare_add_operation();

  // save some cycles for zero additions, but
  // only if it is safe for the matrix we are
  // working with
  typedef typename MatrixType::Traits MatrixTraits;
  if ((MatrixTraits::zero_addition_can_be_elided == true)
      &&
      (factor == 0))
    return;

  for (unsigned int row=0; row<n_block_rows(); ++row)
    for (unsigned int col=0; col<n_block_cols(); ++col)
      // This function should throw if the sparsity
      // patterns of the two blocks differ
      block(row, col).add(factor, matrix.block(row,col));
}



template <class MatrixType>
inline
typename BlockMatrixBase<MatrixType>::value_type
BlockMatrixBase<MatrixType>::operator () (const size_type i,
                                          const size_type j) const
{
  const std::pair<unsigned int,size_type>
  row_index = row_block_indices.global_to_local (i),
  col_index = column_block_indices.global_to_local (j);
  return block(row_index.first,col_index.first) (row_index.second,
                                                 col_index.second);
}



template <class MatrixType>
inline
typename BlockMatrixBase<MatrixType>::value_type
BlockMatrixBase<MatrixType>::el (const size_type i,
                                 const size_type j) const
{
  const std::pair<unsigned int,size_type>
  row_index = row_block_indices.global_to_local (i),
  col_index = column_block_indices.global_to_local (j);
  return block(row_index.first,col_index.first).el (row_index.second,
                                                    col_index.second);
}



template <class MatrixType>
inline
typename BlockMatrixBase<MatrixType>::value_type
BlockMatrixBase<MatrixType>::diag_element (const size_type i) const
{
  Assert (n_block_rows() == n_block_cols(),
          ExcNotQuadratic());

  const std::pair<unsigned int,size_type>
  index = row_block_indices.global_to_local (i);
  return block(index.first,index.first).diag_element(index.second);
}



template <class MatrixType>
inline
void
BlockMatrixBase<MatrixType>::compress (::VectorOperation::values operation)
{
  for (unsigned int r=0; r<n_block_rows(); ++r)
    for (unsigned int c=0; c<n_block_cols(); ++c)
      block(r,c).compress (operation);
}



template <class MatrixType>
inline
BlockMatrixBase<MatrixType> &
BlockMatrixBase<MatrixType>::operator *= (const value_type factor)
{
  Assert (n_block_cols() != 0, ExcNotInitialized());
  Assert (n_block_rows() != 0, ExcNotInitialized());

  for (unsigned int r=0; r<n_block_rows(); ++r)
    for (unsigned int c=0; c<n_block_cols(); ++c)
      block(r,c) *= factor;

  return *this;
}



template <class MatrixType>
inline
BlockMatrixBase<MatrixType> &
BlockMatrixBase<MatrixType>::operator /= (const value_type factor)
{
  Assert (n_block_cols() != 0, ExcNotInitialized());
  Assert (n_block_rows() != 0, ExcNotInitialized());
  Assert (factor !=0, ExcDivideByZero());

  const value_type factor_inv = 1. / factor;

  for (unsigned int r=0; r<n_block_rows(); ++r)
    for (unsigned int c=0; c<n_block_cols(); ++c)
      block(r,c) *= factor_inv;

  return *this;
}



template <class MatrixType>
const BlockIndices &
BlockMatrixBase<MatrixType>::get_row_indices () const
{
  return this->row_block_indices;
}



template <class MatrixType>
const BlockIndices &
BlockMatrixBase<MatrixType>::get_column_indices () const
{
  return this->column_block_indices;
}



template <class MatrixType>
template <class BlockVectorType>
void
BlockMatrixBase<MatrixType>::
vmult_block_block (BlockVectorType       &dst,
                   const BlockVectorType &src) const
{
  Assert (dst.n_blocks() == n_block_rows(),
          ExcDimensionMismatch(dst.n_blocks(), n_block_rows()));
  Assert (src.n_blocks() == n_block_cols(),
          ExcDimensionMismatch(src.n_blocks(), n_block_cols()));

  for (size_type row=0; row<n_block_rows(); ++row)
    {
      block(row,0).vmult (dst.block(row),
                          src.block(0));
      for (size_type col=1; col<n_block_cols(); ++col)
        block(row,col).vmult_add (dst.block(row),
                                  src.block(col));
    };
}



template <class MatrixType>
template <class BlockVectorType,
          class VectorType>
void
BlockMatrixBase<MatrixType>::
vmult_nonblock_block (VectorType    &dst,
                      const BlockVectorType &src) const
{
  Assert (n_block_rows() == 1,
          ExcDimensionMismatch(1, n_block_rows()));
  Assert (src.n_blocks() == n_block_cols(),
          ExcDimensionMismatch(src.n_blocks(), n_block_cols()));

  block(0,0).vmult (dst, src.block(0));
  for (size_type col=1; col<n_block_cols(); ++col)
    block(0,col).vmult_add (dst, src.block(col));
}



template <class MatrixType>
template <class BlockVectorType,
          class VectorType>
void
BlockMatrixBase<MatrixType>::
vmult_block_nonblock (BlockVectorType  &dst,
                      const VectorType &src) const
{
  Assert (dst.n_blocks() == n_block_rows(),
          ExcDimensionMismatch(dst.n_blocks(), n_block_rows()));
  Assert (1 == n_block_cols(),
          ExcDimensionMismatch(1, n_block_cols()));

  for (size_type row=0; row<n_block_rows(); ++row)
    block(row,0).vmult (dst.block(row),
                        src);
}



template <class MatrixType>
template <class VectorType>
void
BlockMatrixBase<MatrixType>::
vmult_nonblock_nonblock (VectorType       &dst,
                         const VectorType &src) const
{
  Assert (1 == n_block_rows(),
          ExcDimensionMismatch(1, n_block_rows()));
  Assert (1 == n_block_cols(),
          ExcDimensionMismatch(1, n_block_cols()));

  block(0,0).vmult (dst, src);
}



template <class MatrixType>
template <class BlockVectorType>
void
BlockMatrixBase<MatrixType>::vmult_add (BlockVectorType       &dst,
                                        const BlockVectorType &src) const
{
  Assert (dst.n_blocks() == n_block_rows(),
          ExcDimensionMismatch(dst.n_blocks(), n_block_rows()));
  Assert (src.n_blocks() == n_block_cols(),
          ExcDimensionMismatch(src.n_blocks(), n_block_cols()));

  for (unsigned int row=0; row<n_block_rows(); ++row)
    for (unsigned int col=0; col<n_block_cols(); ++col)
      block(row,col).vmult_add (dst.block(row),
                                src.block(col));
}




template <class MatrixType>
template <class BlockVectorType>
void
BlockMatrixBase<MatrixType>::
Tvmult_block_block (BlockVectorType       &dst,
                    const BlockVectorType &src) const
{
  Assert (dst.n_blocks() == n_block_cols(),
          ExcDimensionMismatch(dst.n_blocks(), n_block_cols()));
  Assert (src.n_blocks() == n_block_rows(),
          ExcDimensionMismatch(src.n_blocks(), n_block_rows()));

  dst = 0.;

  for (unsigned int  row=0; row<n_block_rows(); ++row)
    {
      for (unsigned int col=0; col<n_block_cols(); ++col)
        block(row,col).Tvmult_add (dst.block(col),
                                   src.block(row));
    };
}



template <class MatrixType>
template <class BlockVectorType,
          class VectorType>
void
BlockMatrixBase<MatrixType>::
Tvmult_block_nonblock (BlockVectorType  &dst,
                       const VectorType &src) const
{
  Assert (dst.n_blocks() == n_block_cols(),
          ExcDimensionMismatch(dst.n_blocks(), n_block_cols()));
  Assert (1 == n_block_rows(),
          ExcDimensionMismatch(1, n_block_rows()));

  dst = 0.;

  for (unsigned int col=0; col<n_block_cols(); ++col)
    block(0,col).Tvmult_add (dst.block(col), src);
}



template <class MatrixType>
template <class BlockVectorType,
          class VectorType>
void
BlockMatrixBase<MatrixType>::
Tvmult_nonblock_block (VectorType    &dst,
                       const BlockVectorType &src) const
{
  Assert (1 == n_block_cols(),
          ExcDimensionMismatch(1, n_block_cols()));
  Assert (src.n_blocks() == n_block_rows(),
          ExcDimensionMismatch(src.n_blocks(), n_block_rows()));

  block(0,0).Tvmult (dst, src.block(0));

  for (size_type row=1; row<n_block_rows(); ++row)
    block(row,0).Tvmult_add (dst, src.block(row));
}



template <class MatrixType>
template <class VectorType>
void
BlockMatrixBase<MatrixType>::
Tvmult_nonblock_nonblock (VectorType       &dst,
                          const VectorType &src) const
{
  Assert (1 == n_block_cols(),
          ExcDimensionMismatch(1, n_block_cols()));
  Assert (1 == n_block_rows(),
          ExcDimensionMismatch(1, n_block_rows()));

  block(0,0).Tvmult (dst, src);
}



template <class MatrixType>
template <class BlockVectorType>
void
BlockMatrixBase<MatrixType>::Tvmult_add (BlockVectorType &dst,
                                         const BlockVectorType &src) const
{
  Assert (dst.n_blocks() == n_block_cols(),
          ExcDimensionMismatch(dst.n_blocks(), n_block_cols()));
  Assert (src.n_blocks() == n_block_rows(),
          ExcDimensionMismatch(src.n_blocks(), n_block_rows()));

  for (unsigned int row=0; row<n_block_rows(); ++row)
    for (unsigned int col=0; col<n_block_cols(); ++col)
      block(row,col).Tvmult_add (dst.block(col),
                                 src.block(row));
}



template <class MatrixType>
template <class BlockVectorType>
typename BlockMatrixBase<MatrixType>::value_type
BlockMatrixBase<MatrixType>::matrix_norm_square (const BlockVectorType &v) const
{
  Assert (n_block_rows() == n_block_cols(), ExcNotQuadratic());
  Assert (v.n_blocks() == n_block_rows(),
          ExcDimensionMismatch(v.n_blocks(), n_block_rows()));

  value_type norm_sqr = 0;
  for (unsigned int row=0; row<n_block_rows(); ++row)
    for (unsigned int col=0; col<n_block_cols(); ++col)
      if (row==col)
        norm_sqr += block(row,col).matrix_norm_square (v.block(row));
      else
        norm_sqr += block(row,col).matrix_scalar_product (v.block(row),
                                                          v.block(col));
  return norm_sqr;
}



template <class MatrixType>
template <class BlockVectorType>
typename BlockMatrixBase<MatrixType>::value_type
BlockMatrixBase<MatrixType>::
matrix_scalar_product (const BlockVectorType    &u,
                       const BlockVectorType &v) const
{
  Assert (u.n_blocks() == n_block_rows(),
          ExcDimensionMismatch(u.n_blocks(), n_block_rows()));
  Assert (v.n_blocks() == n_block_cols(),
          ExcDimensionMismatch(v.n_blocks(), n_block_cols()));

  value_type result = 0;
  for (unsigned int row=0; row<n_block_rows(); ++row)
    for (unsigned int col=0; col<n_block_cols(); ++col)
      result += block(row,col).matrix_scalar_product (u.block(row),
                                                      v.block(col));
  return result;
}



template <class MatrixType>
template <class BlockVectorType>
typename BlockMatrixBase<MatrixType>::value_type
BlockMatrixBase<MatrixType>::
residual (BlockVectorType          &dst,
          const BlockVectorType &x,
          const BlockVectorType    &b) const
{
  Assert (dst.n_blocks() == n_block_rows(),
          ExcDimensionMismatch(dst.n_blocks(), n_block_rows()));
  Assert (b.n_blocks() == n_block_rows(),
          ExcDimensionMismatch(b.n_blocks(), n_block_rows()));
  Assert (x.n_blocks() == n_block_cols(),
          ExcDimensionMismatch(x.n_blocks(), n_block_cols()));
  // in block notation, the residual is
  // r_i = b_i - \sum_j A_ij x_j.
  // this can be written as
  // r_i = b_i - A_i0 x_0 - \sum_{j>0} A_ij x_j.
  //
  // for the first two terms, we can
  // call the residual function of
  // A_i0. for the other terms, we
  // use vmult_add. however, we want
  // to subtract, so in order to
  // avoid a temporary vector, we
  // perform a sign change of the
  // first two term before, and after
  // adding up
  for (unsigned int row=0; row<n_block_rows(); ++row)
    {
      block(row,0).residual (dst.block(row),
                             x.block(0),
                             b.block(row));

      for (size_type i=0; i<dst.block(row).size(); ++i)
        dst.block(row)(i) = -dst.block(row)(i);

      for (unsigned int col=1; col<n_block_cols(); ++col)
        block(row,col).vmult_add (dst.block(row),
                                  x.block(col));

      for (size_type i=0; i<dst.block(row).size(); ++i)
        dst.block(row)(i) = -dst.block(row)(i);
    };

  value_type res = 0;
  for (size_type row=0; row<n_block_rows(); ++row)
    res += dst.block(row).norm_sqr ();
  return std::sqrt(res);
}



template <class MatrixType>
inline
void
BlockMatrixBase<MatrixType>::print (std::ostream &out,
                                    const bool    alternative_output) const
{
  for (unsigned int row=0; row<n_block_rows(); ++row)
    for (unsigned int col=0; col<n_block_cols(); ++col)
      {
        if (!alternative_output)
          out << "Block (" << row << ", " << col << ")" << std::endl;

        block(row, col).print(out, alternative_output);
      }
}



template <class MatrixType>
inline
typename BlockMatrixBase<MatrixType>::const_iterator
BlockMatrixBase<MatrixType>::begin () const
{
  return const_iterator(this, 0);
}



template <class MatrixType>
inline
typename BlockMatrixBase<MatrixType>::const_iterator
BlockMatrixBase<MatrixType>::end () const
{
  return const_iterator(this, m());
}



template <class MatrixType>
inline
typename BlockMatrixBase<MatrixType>::const_iterator
BlockMatrixBase<MatrixType>::begin (const size_type r) const
{
  Assert (r<m(), ExcIndexRange(r,0,m()));
  return const_iterator(this, r);
}



template <class MatrixType>
inline
typename BlockMatrixBase<MatrixType>::const_iterator
BlockMatrixBase<MatrixType>::end (const size_type r) const
{
  Assert (r<m(), ExcIndexRange(r,0,m()));
  return const_iterator(this, r+1);
}



template <class MatrixType>
inline
typename BlockMatrixBase<MatrixType>::iterator
BlockMatrixBase<MatrixType>::begin ()
{
  return iterator(this, 0);
}



template <class MatrixType>
inline
typename BlockMatrixBase<MatrixType>::iterator
BlockMatrixBase<MatrixType>::end ()
{
  return iterator(this, m());
}



template <class MatrixType>
inline
typename BlockMatrixBase<MatrixType>::iterator
BlockMatrixBase<MatrixType>::begin (const size_type r)
{
  Assert (r<m(), ExcIndexRange(r,0,m()));
  return iterator(this, r);
}



template <class MatrixType>
inline
typename BlockMatrixBase<MatrixType>::iterator
BlockMatrixBase<MatrixType>::end (const size_type r)
{
  Assert (r<m(), ExcIndexRange(r,0,m()));
  return iterator(this, r+1);
}



template <class MatrixType>
void
BlockMatrixBase<MatrixType>::collect_sizes ()
{
  std::vector<size_type> row_sizes (this->n_block_rows());
  std::vector<size_type> col_sizes (this->n_block_cols());

  // first find out the row sizes
  // from the first block column
  for (unsigned int r=0; r<this->n_block_rows(); ++r)
    row_sizes[r] = sub_objects[r][0]->m();
  // then check that the following
  // block columns have the same
  // sizes
  for (unsigned int c=1; c<this->n_block_cols(); ++c)
    for (unsigned int r=0; r<this->n_block_rows(); ++r)
      Assert (row_sizes[r] == sub_objects[r][c]->m(),
              ExcIncompatibleRowNumbers (r,0,r,c));

  // finally initialize the row
  // indices with this array
  this->row_block_indices.reinit (row_sizes);


  // then do the same with the columns
  for (unsigned int c=0; c<this->n_block_cols(); ++c)
    col_sizes[c] = sub_objects[0][c]->n();
  for (unsigned int r=1; r<this->n_block_rows(); ++r)
    for (unsigned int c=0; c<this->n_block_cols(); ++c)
      Assert (col_sizes[c] == sub_objects[r][c]->n(),
              ExcIncompatibleRowNumbers (0,c,r,c));

  // finally initialize the row
  // indices with this array
  this->column_block_indices.reinit (col_sizes);
}



template <class MatrixType>
void
BlockMatrixBase<MatrixType>::prepare_add_operation ()
{
  for (unsigned int row=0; row<n_block_rows(); ++row)
    for (unsigned int col=0; col<n_block_cols(); ++col)
      block(row, col).prepare_add();
}



template <class MatrixType>
void
BlockMatrixBase<MatrixType>::prepare_set_operation ()
{
  for (unsigned int row=0; row<n_block_rows(); ++row)
    for (unsigned int col=0; col<n_block_cols(); ++col)
      block(row, col).prepare_set();
}

#endif // DOXYGEN


DEAL_II_NAMESPACE_CLOSE

#endif    // dealii_block_matrix_base_h