From: turcksin <turcksin@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Thu, 9 May 2013 19:11:46 +0000 (+0000)
Subject: Merge from mainline.
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=d8b4faa90f6e52f907930fa1cf03736b0f36f7af;p=dealii-svn.git

Merge from mainline.

git-svn-id: https://svn.dealii.org/branches/branch_bigger_global_dof_indices_4@29488 0785d39b-7218-0410-832d-ea1e28bc413d
---

d8b4faa90f6e52f907930fa1cf03736b0f36f7af
diff --cc deal.II/include/deal.II/lac/chunk_sparse_matrix.h
index bb3b4d8ce8,8ceb359cfa..5ecc459342
--- a/deal.II/include/deal.II/lac/chunk_sparse_matrix.h
+++ b/deal.II/include/deal.II/lac/chunk_sparse_matrix.h
@@@ -290,57 -573,40 +578,40 @@@ public
    bool empty () const;
  
    /**
-    * Return the dimension of the
-    * image space.  To remember: the
-    * matrix is of dimension
-    * $m \times n$.
+    * Return the dimension of the image space.  To remember: the matrix is of
+    * dimension $m \times n$.
     */
 -  unsigned int m () const;
 +  size_type m () const;
  
    /**
-    * Return the dimension of the
-    * range space.  To remember: the
-    * matrix is of dimension
-    * $m \times n$.
+    * Return the dimension of the range space.  To remember: the matrix is of
+    * dimension $m \times n$.
     */
 -  unsigned int n () const;
 +  size_type n () const;
  
    /**
-    * Return the number of nonzero
-    * elements of this
-    * matrix. Actually, it returns
-    * the number of entries in the
-    * sparsity pattern; if any of
-    * the entries should happen to
-    * be zero, it is counted anyway.
+    * Return the number of nonzero elements of this matrix. Actually, it
+    * returns the number of entries in the sparsity pattern; if any of the
+    * entries should happen to be zero, it is counted anyway.
     */
 -  unsigned int n_nonzero_elements () const;
 +  size_type n_nonzero_elements () const;
  
    /**
-    * Return the number of actually
-    * nonzero elements of this
-    * matrix.
+    * Return the number of actually nonzero elements of this matrix.
     *
-    * Note, that this function does
-    * (in contrary to
-    * n_nonzero_elements()) not
-    * count all entries of the
-    * sparsity pattern but only the
-    * ones that are nonzero.
+    * Note, that this function does (in contrary to n_nonzero_elements()) not
+    * count all entries of the sparsity pattern but only the ones that are
+    * nonzero.
     */
 -  unsigned int n_actually_nonzero_elements () const;
 +  size_type n_actually_nonzero_elements () const;
  
    /**
-    * Return a (constant) reference
-    * to the underlying sparsity
-    * pattern of this matrix.
+    * Return a (constant) reference to the underlying sparsity pattern of this
+    * matrix.
     *
-    * Though the return value is
-    * declared <tt>const</tt>, you
-    * should be aware that it may
-    * change if you call any
-    * nonconstant function of
-    * objects which operate on it.
+    * Though the return value is declared <tt>const</tt>, you should be aware
+    * that it may change if you call any nonconstant function of objects which
+    * operate on it.
     */
    const ChunkSparsityPattern &get_sparsity_pattern () const;
  
@@@ -358,59 -622,42 +627,42 @@@
     */
  //@{
    /**
-    * Set the element (<i>i,j</i>)
-    * to <tt>value</tt>. Throws an
-    * error if the entry does not
-    * exist or if <tt>value</tt> is
-    * not a finite number. Still, it
-    * is allowed to store zero
-    * values in non-existent fields.
+    * Set the element (<i>i,j</i>) to <tt>value</tt>. Throws an error if the
+    * entry does not exist or if <tt>value</tt> is not a finite number. Still,
+    * it is allowed to store zero values in non-existent fields.
     */
 -  void set (const unsigned int i,
 -            const unsigned int j,
 +  void set (const size_type i,
 +            const size_type j,
              const number value);
  
    /**
-    * Add <tt>value</tt> to the
-    * element (<i>i,j</i>).  Throws
-    * an error if the entry does not
-    * exist or if <tt>value</tt> is
-    * not a finite number. Still, it
-    * is allowed to store zero
-    * values in non-existent fields.
+    * Add <tt>value</tt> to the element (<i>i,j</i>).  Throws an error if the
+    * entry does not exist or if <tt>value</tt> is not a finite number. Still,
+    * it is allowed to store zero values in non-existent fields.
     */
 -  void add (const unsigned int i,
 -            const unsigned int j,
 +  void add (const size_type i,
 +            const size_type j,
              const number value);
  
    /**
-    * Add an array of values given by
-    * <tt>values</tt> in the given
-    * global matrix row at columns
-    * specified by col_indices in the
-    * sparse matrix.
+    * Add an array of values given by <tt>values</tt> in the given global
+    * matrix row at columns specified by col_indices in the sparse matrix.
     *
-    * The optional parameter
-    * <tt>elide_zero_values</tt> can be
-    * used to specify whether zero
-    * values should be added anyway or
-    * these should be filtered away and
-    * only non-zero data is added. The
-    * default value is <tt>true</tt>,
-    * i.e., zero values won't be added
-    * into the matrix.
+    * The optional parameter <tt>elide_zero_values</tt> can be used to specify
+    * whether zero values should be added anyway or these should be filtered
+    * away and only non-zero data is added. The default value is <tt>true</tt>,
+    * i.e., zero values won't be added into the matrix.
     */
    template <typename number2>
 -  void add (const unsigned int  row,
 -            const unsigned int  n_cols,
 -            const unsigned int *col_indices,
 -            const number2      *values,
 -            const bool          elide_zero_values = true,
 -            const bool          col_indices_are_sorted = false);
 +  void add (const size_type  row,
 +            const size_type  n_cols,
 +            const size_type *col_indices,
 +            const number2   *values,
 +            const bool       elide_zero_values = true,
 +            const bool       col_indices_are_sorted = false);
  
    /**
-    * Multiply the entire matrix by a
-    * fixed factor.
+    * Multiply the entire matrix by a fixed factor.
     */
    ChunkSparseMatrix &operator *= (const number factor);
  
@@@ -556,58 -748,35 +753,35 @@@
  //@{
  
    /**
-    * Return the value of the entry
-    * (<i>i,j</i>).  This may be an
-    * expensive operation and you
-    * should always take care where
-    * to call this function.  In
-    * order to avoid abuse, this
-    * function throws an exception
-    * if the required element does
-    * not exist in the matrix.
+    * Return the value of the entry (<i>i,j</i>).  This may be an expensive
+    * operation and you should always take care where to call this function.
+    * In order to avoid abuse, this function throws an exception if the
+    * required element does not exist in the matrix.
     *
-    * In case you want a function
-    * that returns zero instead (for
-    * entries that are not in the
-    * sparsity pattern of the
-    * matrix), use the el()
-    * function.
+    * In case you want a function that returns zero instead (for entries that
+    * are not in the sparsity pattern of the matrix), use the el() function.
     *
-    * If you are looping over all elements,
-    * consider using one of the iterator
-    * classes instead, since they are
-    * tailored better to a sparse matrix
+    * If you are looping over all elements, consider using one of the iterator
+    * classes instead, since they are tailored better to a sparse matrix
     * structure.
     */
 -  number operator () (const unsigned int i,
 -                      const unsigned int j) const;
 +  number operator () (const size_type i,
 +                      const size_type j) const;
  
    /**
-    * This function is mostly like
-    * operator()() in that it
-    * returns the value of the
-    * matrix entry (<i>i,j</i>). The
-    * only difference is that if
-    * this entry does not exist in
-    * the sparsity pattern, then
-    * instead of raising an
-    * exception, zero is
-    * returned. While this may be
-    * convenient in some cases, note
-    * that it is simple to write
-    * algorithms that are slow
-    * compared to an optimal
-    * solution, since the sparsity
-    * of the matrix is not used.
+    * This function is mostly like operator()() in that it returns the value of
+    * the matrix entry (<i>i,j</i>). The only difference is that if this entry
+    * does not exist in the sparsity pattern, then instead of raising an
+    * exception, zero is returned. While this may be convenient in some cases,
+    * note that it is simple to write algorithms that are slow compared to an
+    * optimal solution, since the sparsity of the matrix is not used.
     *
-    * If you are looping over all elements,
-    * consider using one of the iterator
-    * classes instead, since they are
-    * tailored better to a sparse matrix
+    * If you are looping over all elements, consider using one of the iterator
+    * classes instead, since they are tailored better to a sparse matrix
     * structure.
     */
 -  number el (const unsigned int i,
 -             const unsigned int j) const;
 +  number el (const size_type i,
 +             const size_type j) const;
  
    /**
     * Return the main diagonal
@@@ -616,23 -785,18 +790,18 @@@
     * error if the matrix is not
     * quadratic.
     *
-    * This function is considerably
-    * faster than the operator()(),
-    * since for quadratic matrices, the
-    * diagonal entry may be the
-    * first to be stored in each row
-    * and access therefore does not
-    * involve searching for the
-    * right column number.
+    * This function is considerably faster than the operator()(), since for
+    * quadratic matrices, the diagonal entry may be the first to be stored in
+    * each row and access therefore does not involve searching for the right
+    * column number.
     */
 -  number diag_element (const unsigned int i) const;
 +  number diag_element (const size_type i) const;
  
    /**
-    * Same as above, but return a
-    * writeable reference. You're
-    * sure you know what you do?
+    * Same as above, but return a writeable reference. You're sure you know
+    * what you do?
     */
 -  number &diag_element (const unsigned int i);
 +  number &diag_element (const size_type i);
  
  //@}
    /**
@@@ -1197,28 -1298,28 +1303,28 @@@ private
    number *val;
  
    /**
-    * Allocated size of #val. This
-    * can be larger than the
-    * actually used part if the size
-    * of the matrix was reduced
-    * somewhen in the past by
-    * associating a sparsity pattern
-    * with a smaller size to this
-    * object, using the reinit()
+    * Allocated size of #val. This can be larger than the actually used part if
+    * the size of the matrix was reduced somewhen in the past by associating a
+    * sparsity pattern with a smaller size to this object, using the reinit()
     * function.
     */
 -  unsigned int max_len;
 +  size_type max_len;
  
    /**
-    * Return the location of entry
-    * $(i,j)$ within the val array.
+    * Return the location of entry $(i,j)$ within the val array.
     */
 -  unsigned int compute_location (const unsigned int i,
 -                                 const unsigned int j) const;
 +  size_type compute_location (const size_type i,
 +                              const size_type j) const;
  
-   // make all other sparse matrices
-   // friends
+   // make all other sparse matrices friends
    template <typename somenumber> friend class ChunkSparseMatrix;
+ 
+   /**
+    * Also give access to internal details to the iterator/accessor
+    * classes.
+    */
+   template <typename,bool> friend class ChunkSparseMatrixIterators::Iterator;
+   template <typename,bool> friend class ChunkSparseMatrixIterators::Accessor;
  };
  
  /*@}*/
@@@ -1249,15 -1348,26 +1355,26 @@@ ChunkSparseMatrix<number>::n () cons
  
  
  
+ template <typename number>
+ inline
+ const ChunkSparsityPattern &
+ ChunkSparseMatrix<number>::get_sparsity_pattern () const
+ {
+   Assert (cols != 0, ExcNotInitialized());
+   return *cols;
+ }
+ 
+ 
+ 
  template <typename number>
  inline
 -unsigned int
 -ChunkSparseMatrix<number>::compute_location (const unsigned int i,
 -                                             const unsigned int j) const
 +typename ChunkSparseMatrix<number>::size_type 
 +ChunkSparseMatrix<number>::compute_location (const size_type i,
 +                                             const size_type j) const
  {
 -  const unsigned int chunk_size = cols->get_chunk_size();
 -  const unsigned int chunk_index
 +  const size_type chunk_size = cols->get_chunk_size();
 +  const size_type chunk_index
-     = cols->sparsity_pattern(i/chunk_size, j/chunk_size);
+   = cols->sparsity_pattern(i/chunk_size, j/chunk_size);
  
    if (chunk_index == ChunkSparsityPattern::invalid_entry)
      return ChunkSparsityPattern::invalid_entry;
@@@ -1282,11 -1392,9 +1399,9 @@@ void ChunkSparseMatrix<number>::set (co
    Assert (numbers::is_finite(value), ExcNumberNotFinite());
  
    Assert (cols != 0, ExcNotInitialized());
-   // it is allowed to set elements of
-   // the matrix that are not part of
-   // the sparsity pattern, if the
-   // value to which we set it is zero
+   // it is allowed to set elements of the matrix that are not part of the
+   // sparsity pattern, if the value to which we set it is zero
 -  const unsigned int index = compute_location(i,j);
 +  const size_type index = compute_location(i,j);
    Assert ((index != SparsityPattern::invalid_entry) ||
            (value == 0.),
            ExcInvalidIndex(i,j));
@@@ -1345,16 -1453,12 +1460,12 @@@ ChunkSparseMatrix<number>::operator *= 
    Assert (cols != 0, ExcNotInitialized());
    Assert (val != 0, ExcNotInitialized());
  
 -  const unsigned int chunk_size = cols->get_chunk_size();
 +  const size_type chunk_size = cols->get_chunk_size();
  
-   // multiply all elements of the matrix with
-   // the given factor. this includes the
-   // padding elements in chunks that overlap
-   // the boundaries of the actual matrix --
-   // but since multiplication with a number
-   // does not violate the invariant of
-   // keeping these elements at zero nothing
-   // can happen
+   // multiply all elements of the matrix with the given factor. this includes
+   // the padding elements in chunks that overlap the boundaries of the actual
+   // matrix -- but since multiplication with a number does not violate the
+   // invariant of keeping these elements at zero nothing can happen
    number             *val_ptr    = val;
    const number *const end_ptr    = val +
                                     cols->sparsity_pattern.n_nonzero_elements()
@@@ -1379,16 -1483,12 +1490,12 @@@ ChunkSparseMatrix<number>::operator /= 
  
    const number factor_inv = 1. / factor;
  
 -  const unsigned int chunk_size = cols->get_chunk_size();
 +  const size_type chunk_size = cols->get_chunk_size();
  
-   // multiply all elements of the matrix with
-   // the given factor. this includes the
-   // padding elements in chunks that overlap
-   // the boundaries of the actual matrix --
-   // but since multiplication with a number
-   // does not violate the invariant of
-   // keeping these elements at zero nothing
-   // can happen
+   // multiply all elements of the matrix with the given factor. this includes
+   // the padding elements in chunks that overlap the boundaries of the actual
+   // matrix -- but since multiplication with a number does not violate the
+   // invariant of keeping these elements at zero nothing can happen
    number             *val_ptr    = val;
    const number *const end_ptr    = val +
                                     cols->sparsity_pattern.n_nonzero_elements()
@@@ -1440,33 -1540,9 +1547,9 @@@ number ChunkSparseMatrix<number>::diag_
    Assert (m() == n(),  ExcNotQuadratic());
    Assert (i<m(), ExcInvalidIndex1(i));
  
-   // Use that the first element in each row
-   // of a quadratic matrix is the main
-   // diagonal of the chunk sparsity pattern
-   const size_type chunk_size = cols->get_chunk_size();
-   return val[cols->sparsity_pattern.rowstart[i/chunk_size]
-              *
-              chunk_size * chunk_size
-              +
-              (i % chunk_size) * chunk_size
-              +
-              (i % chunk_size)];
- }
- 
- 
- 
- template <typename number>
- inline
- number &ChunkSparseMatrix<number>::diag_element (const size_type i)
- {
-   Assert (cols != 0, ExcNotInitialized());
-   Assert (m() == n(),  ExcNotQuadratic());
-   Assert (i<m(), ExcInvalidIndex1(i));
- 
-   // Use that the first element in each row
-   // of a quadratic matrix is the main
+   // Use that the first element in each row of a quadratic matrix is the main
    // diagonal of the chunk sparsity pattern
 -  const unsigned int chunk_size = cols->get_chunk_size();
 +  const size_type chunk_size = cols->get_chunk_size();
    return val[cols->sparsity_pattern.rowstart[i/chunk_size]
               *
               chunk_size * chunk_size
@@@ -1484,14 -1561,13 +1568,13 @@@ voi
  ChunkSparseMatrix<number>::copy_from (const ForwardIterator begin,
                                        const ForwardIterator end)
  {
 -  Assert (static_cast<unsigned int>(std::distance (begin, end)) == m(),
 +  Assert (static_cast<size_type >(std::distance (begin, end)) == m(),
            ExcIteratorRange (std::distance (begin, end), m()));
  
-   // for use in the inner loop, we
-   // define a typedef to the type of
-   // the inner iterators
+   // for use in the inner loop, we define a typedef to the type of the inner
+   // iterators
    typedef typename std::iterator_traits<ForwardIterator>::value_type::const_iterator inner_iterator;
 -  unsigned int row=0;
 +  size_type row=0;
    for (ForwardIterator i=begin; i!=end; ++i, ++row)
      {
        const inner_iterator end_of_row = i->end();
diff --cc deal.II/include/deal.II/lac/chunk_sparse_matrix.templates.h
index 19ada3ca0f,bbc1fb9944..d6fad7a4d6
--- a/deal.II/include/deal.II/lac/chunk_sparse_matrix.templates.h
+++ b/deal.II/include/deal.II/lac/chunk_sparse_matrix.templates.h
@@@ -318,12 -437,10 +444,10 @@@ ChunkSparseMatrix<number>::reinit (cons
        return;
      }
  
-   // allocate not just m() * n() elements but
-   // enough so that we can store full
-   // chunks. this entails some padding
-   // elements
+   // allocate not just m() * n() elements but enough so that we can store full
+   // chunks. this entails some padding elements
 -  const unsigned int chunk_size = cols->get_chunk_size();
 -  const unsigned int N = cols->sparsity_pattern.n_nonzero_elements() *
 +  const size_type chunk_size = cols->get_chunk_size();
 +  const size_type N = cols->sparsity_pattern.n_nonzero_elements() *
                           chunk_size * chunk_size;
    if (N > max_len || max_len == 0)
      {
@@@ -383,12 -498,10 +505,10 @@@ ChunkSparseMatrix<number>::n_actually_n
  {
    Assert (cols != 0, ExcNotInitialized());
  
-   // count those elements that are nonzero,
-   // even if they lie in the padding around
-   // the matrix. since we have the invariant
-   // that padding elements are zero, nothing
-   // bad can happen here
+   // count those elements that are nonzero, even if they lie in the padding
+   // around the matrix. since we have the invariant that padding elements are
+   // zero, nothing bad can happen here
 -  const unsigned int chunk_size = cols->get_chunk_size();
 +  const size_type chunk_size = cols->get_chunk_size();
    return std::count_if(&val[0],
                         &val[cols->sparsity_pattern.n_nonzero_elements () *
                              chunk_size * chunk_size],
@@@ -418,9 -531,8 +538,8 @@@ ChunkSparseMatrix<number>::copy_from (c
    Assert (val != 0, ExcNotInitialized());
    Assert (cols == matrix.cols, ExcDifferentChunkSparsityPatterns());
  
-   // copy everything, including padding
-   // elements
+   // copy everything, including padding elements
 -  const unsigned int chunk_size = cols->get_chunk_size();
 +  const size_type chunk_size = cols->get_chunk_size();
    std::copy (&matrix.val[0],
               &matrix.val[cols->sparsity_pattern.n_nonzero_elements()
                           * chunk_size * chunk_size],
@@@ -458,9 -570,8 +577,8 @@@ ChunkSparseMatrix<number>::add (const n
    Assert (val != 0, ExcNotInitialized());
    Assert (cols == matrix.cols, ExcDifferentChunkSparsityPatterns());
  
-   // add everything, including padding
-   // elements
+   // add everything, including padding elements
 -  const unsigned int chunk_size = cols->get_chunk_size();
 +  const size_type     chunk_size = cols->get_chunk_size();
    number             *val_ptr    = &val[0];
    const somenumber   *matrix_ptr = &matrix.val[0];
    const number *const end_ptr    = &val[cols->sparsity_pattern.n_nonzero_elements()
@@@ -639,12 -671,10 +678,10 @@@ ChunkSparseMatrix<number>::Tvmult_add (
  
    Assert (!PointerComparison::equal(&src, &dst), ExcSourceEqualsDestination());
  
 -  const unsigned int n_chunk_rows = cols->sparsity_pattern.n_rows();
 +  const size_type n_chunk_rows = cols->sparsity_pattern.n_rows();
  
-   // loop over all chunks. note that we need
-   // to treat the last chunk row and column
-   // differently if they have padding
-   // elements
+   // loop over all chunks. note that we need to treat the last chunk row and
+   // column differently if they have padding elements
    const bool rows_have_padding = (m() % cols->chunk_size != 0),
               cols_have_padding = (n() % cols->chunk_size != 0);
  
@@@ -653,13 -683,12 +690,12 @@@
         n_chunk_rows-1 :
         n_chunk_rows);
  
-   // like in vmult_add, but don't keep an
-   // iterator into dst around since we're not
-   // traversing it sequentially this time
+   // like in vmult_add, but don't keep an iterator into dst around since we're
+   // not traversing it sequentially this time
 -  const number       *val_ptr    = val;
 -  const unsigned int *colnum_ptr = cols->sparsity_pattern.colnums;
 +  const number    *val_ptr    = val;
 +  const size_type *colnum_ptr = cols->sparsity_pattern.colnums;
  
 -  for (unsigned int chunk_row=0; chunk_row<n_regular_chunk_rows; ++chunk_row)
 +  for (size_type chunk_row=0; chunk_row<n_regular_chunk_rows; ++chunk_row)
      {
        const number *const val_end_of_row = &val[cols->sparsity_pattern.rowstart[chunk_row+1]
                                                  * cols->chunk_size
@@@ -675,10 -704,9 +711,9 @@@
               src.begin() + chunk_row * cols->chunk_size,
               dst.begin() + *colnum_ptr * cols->chunk_size);
            else
-             // we're at a chunk column that
-             // has padding
+             // we're at a chunk column that has padding
 -            for (unsigned int r=0; r<cols->chunk_size; ++r)
 -              for (unsigned int c=0; c<n() % cols->chunk_size; ++c)
 +            for (size_type r=0; r<cols->chunk_size; ++r)
 +              for (size_type c=0; c<n() % cols->chunk_size; ++c)
                  dst(*colnum_ptr * cols->chunk_size + c)
                  += (val_ptr[r*cols->chunk_size + c] *
                      src(chunk_row * cols->chunk_size + r));
@@@ -688,11 -716,10 +723,10 @@@
          }
      }
  
-   // now deal with last chunk row if
-   // necessary
+   // now deal with last chunk row if necessary
    if (rows_have_padding)
      {
 -      const unsigned int chunk_row = n_chunk_rows - 1;
 +      const size_type chunk_row = n_chunk_rows - 1;
  
        const number *const val_end_of_row = &val[cols->sparsity_pattern.rowstart[chunk_row+1]
                                                  * cols->chunk_size
@@@ -703,19 -730,17 +737,17 @@@
                ||
                (*colnum_ptr != cols->sparsity_pattern.n_cols()-1))
              {
-               // we're at a chunk row but not
-               // column that has padding
+               // we're at a chunk row but not column that has padding
 -              for (unsigned int r=0; r<m() % cols->chunk_size; ++r)
 -                for (unsigned int c=0; c<cols->chunk_size; ++c)
 +              for (size_type r=0; r<m() % cols->chunk_size; ++r)
 +                for (size_type c=0; c<cols->chunk_size; ++c)
                    dst(*colnum_ptr * cols->chunk_size + c)
                    += (val_ptr[r*cols->chunk_size + c] *
                        src(chunk_row * cols->chunk_size + r));
              }
            else
-             // we're at a chunk row and
-             // column that has padding
+             // we're at a chunk row and column that has padding
 -            for (unsigned int r=0; r<m() % cols->chunk_size; ++r)
 -              for (unsigned int c=0; c<n() % cols->chunk_size; ++c)
 +            for (size_type r=0; r<m() % cols->chunk_size; ++r)
 +              for (size_type c=0; c<n() % cols->chunk_size; ++c)
                  dst(*colnum_ptr * cols->chunk_size + c)
                  += (val_ptr[r*cols->chunk_size + c] *
                      src(chunk_row * cols->chunk_size + r));
@@@ -740,15 -765,12 +772,12 @@@ ChunkSparseMatrix<number>::matrix_norm_
    somenumber result = 0;
  
    ////////////////
-   // like matrix_scalar_product, except that
-   // the two vectors are now the same
+   // like matrix_scalar_product, except that the two vectors are now the same
  
 -  const unsigned int n_chunk_rows = cols->sparsity_pattern.n_rows();
 +  const size_type n_chunk_rows = cols->sparsity_pattern.n_rows();
  
-   // loop over all chunks. note that we need
-   // to treat the last chunk row and column
-   // differently if they have padding
-   // elements
+   // loop over all chunks. note that we need to treat the last chunk row and
+   // column differently if they have padding elements
    const bool rows_have_padding = (m() % cols->chunk_size != 0),
               cols_have_padding = (n() % cols->chunk_size != 0);
  
@@@ -779,10 -801,9 +808,9 @@@
                 v_ptr,
                 v.begin() + *colnum_ptr * cols->chunk_size);
            else
-             // we're at a chunk column that
-             // has padding
+             // we're at a chunk column that has padding
 -            for (unsigned int r=0; r<cols->chunk_size; ++r)
 -              for (unsigned int c=0; c<n() % cols->chunk_size; ++c)
 +            for (size_type r=0; r<cols->chunk_size; ++r)
 +              for (size_type c=0; c<n() % cols->chunk_size; ++c)
                  result
                  +=
                    v(chunk_row * cols->chunk_size + r)
@@@ -797,11 -818,10 +825,10 @@@
        v_ptr += cols->chunk_size;
      }
  
-   // now deal with last chunk row if
-   // necessary
+   // now deal with last chunk row if necessary
    if (rows_have_padding)
      {
 -      const unsigned int chunk_row = n_chunk_rows - 1;
 +      const size_type chunk_row = n_chunk_rows - 1;
  
        const number *const val_end_of_row = &val[cols->sparsity_pattern.rowstart[chunk_row+1]
                                                  * cols->chunk_size
@@@ -812,10 -832,9 +839,9 @@@
                ||
                (*colnum_ptr != cols->sparsity_pattern.n_cols()-1))
              {
-               // we're at a chunk row but not
-               // column that has padding
+               // we're at a chunk row but not column that has padding
 -              for (unsigned int r=0; r<m() % cols->chunk_size; ++r)
 -                for (unsigned int c=0; c<cols->chunk_size; ++c)
 +              for (size_type r=0; r<m() % cols->chunk_size; ++r)
 +                for (size_type c=0; c<cols->chunk_size; ++c)
                    result
                    +=
                      v(chunk_row * cols->chunk_size + r)
@@@ -823,10 -842,9 +849,9 @@@
                         v(*colnum_ptr * cols->chunk_size + c));
              }
            else
-             // we're at a chunk row and
-             // column that has padding
+             // we're at a chunk row and column that has padding
 -            for (unsigned int r=0; r<m() % cols->chunk_size; ++r)
 -              for (unsigned int c=0; c<n() % cols->chunk_size; ++c)
 +            for (size_type r=0; r<m() % cols->chunk_size; ++r)
 +              for (size_type c=0; c<n() % cols->chunk_size; ++c)
                  result
                  +=
                    v(chunk_row * cols->chunk_size + r)
@@@ -854,16 -872,13 +879,13 @@@ ChunkSparseMatrix<number>::matrix_scala
    Assert(m() == u.size(), ExcDimensionMismatch(m(),u.size()));
    Assert(n() == v.size(), ExcDimensionMismatch(n(),v.size()));
  
-   // the following works like the vmult_add
-   // function
+   // the following works like the vmult_add function
    somenumber result = 0;
  
 -  const unsigned int n_chunk_rows = cols->sparsity_pattern.n_rows();
 +  const size_type n_chunk_rows = cols->sparsity_pattern.n_rows();
  
-   // loop over all chunks. note that we need
-   // to treat the last chunk row and column
-   // differently if they have padding
-   // elements
+   // loop over all chunks. note that we need to treat the last chunk row and
+   // column differently if they have padding elements
    const bool rows_have_padding = (m() % cols->chunk_size != 0),
               cols_have_padding = (n() % cols->chunk_size != 0);
  
@@@ -894,10 -909,9 +916,9 @@@
                 u_ptr,
                 v.begin() + *colnum_ptr * cols->chunk_size);
            else
-             // we're at a chunk column that
-             // has padding
+             // we're at a chunk column that has padding
 -            for (unsigned int r=0; r<cols->chunk_size; ++r)
 -              for (unsigned int c=0; c<n() % cols->chunk_size; ++c)
 +            for (size_type r=0; r<cols->chunk_size; ++r)
 +              for (size_type c=0; c<n() % cols->chunk_size; ++c)
                  result
                  +=
                    u(chunk_row * cols->chunk_size + r)
@@@ -912,11 -926,10 +933,10 @@@
        u_ptr += cols->chunk_size;
      }
  
-   // now deal with last chunk row if
-   // necessary
+   // now deal with last chunk row if necessary
    if (rows_have_padding)
      {
 -      const unsigned int chunk_row = n_chunk_rows - 1;
 +      const size_type chunk_row = n_chunk_rows - 1;
  
        const number *const val_end_of_row = &val[cols->sparsity_pattern.rowstart[chunk_row+1]
                                                  * cols->chunk_size
@@@ -927,10 -940,9 +947,9 @@@
                ||
                (*colnum_ptr != cols->sparsity_pattern.n_cols()-1))
              {
-               // we're at a chunk row but not
-               // column that has padding
+               // we're at a chunk row but not column that has padding
 -              for (unsigned int r=0; r<m() % cols->chunk_size; ++r)
 -                for (unsigned int c=0; c<cols->chunk_size; ++c)
 +              for (size_type r=0; r<m() % cols->chunk_size; ++r)
 +                for (size_type c=0; c<cols->chunk_size; ++c)
                    result
                    +=
                      u(chunk_row * cols->chunk_size + r)
@@@ -938,10 -950,9 +957,9 @@@
                         v(*colnum_ptr * cols->chunk_size + c));
              }
            else
-             // we're at a chunk row and
-             // column that has padding
+             // we're at a chunk row and column that has padding
 -            for (unsigned int r=0; r<m() % cols->chunk_size; ++r)
 -              for (unsigned int c=0; c<n() % cols->chunk_size; ++c)
 +            for (size_type r=0; r<m() % cols->chunk_size; ++r)
 +              for (size_type c=0; c<n() % cols->chunk_size; ++c)
                  result
                  +=
                    u(chunk_row * cols->chunk_size + r)
@@@ -965,13 -976,11 +983,11 @@@ ChunkSparseMatrix<number>::l1_norm () c
    Assert (cols != 0, ExcNotInitialized());
    Assert (val != 0, ExcNotInitialized());
  
 -  const unsigned int n_chunk_rows = cols->sparsity_pattern.n_rows();
 +  const size_type n_chunk_rows = cols->sparsity_pattern.n_rows();
  
-   // loop over all rows and columns; it is
-   // safe to also loop over the padding
-   // elements (they are zero) if we make sure
-   // that the vector into which we sum column
-   // sums is large enough
+   // loop over all rows and columns; it is safe to also loop over the padding
+   // elements (they are zero) if we make sure that the vector into which we
+   // sum column sums is large enough
    Vector<real_type> column_sums(cols->sparsity_pattern.n_cols() *
                                  cols->chunk_size);
  
@@@ -999,19 -1008,15 +1015,15 @@@ ChunkSparseMatrix<number>::linfty_norm 
    Assert (cols != 0, ExcNotInitialized());
    Assert (val != 0, ExcNotInitialized());
  
-   // this function works like l1_norm(). it
-   // can be made more efficient (without
-   // allocating a temporary vector) as is
-   // done in the SparseMatrix class but since
-   // it is rarely called in time critical
-   // places it is probably not worth it
+   // this function works like l1_norm(). it can be made more efficient
+   // (without allocating a temporary vector) as is done in the SparseMatrix
+   // class but since it is rarely called in time critical places it is
+   // probably not worth it
 -  const unsigned int n_chunk_rows = cols->sparsity_pattern.n_rows();
 +  const size_type n_chunk_rows = cols->sparsity_pattern.n_rows();
  
-   // loop over all rows and columns; it is
-   // safe to also loop over the padding
-   // elements (they are zero) if we make sure
-   // that the vector into which we sum column
-   // sums is large enough
+   // loop over all rows and columns; it is safe to also loop over the padding
+   // elements (they are zero) if we make sure that the vector into which we
+   // sum column sums is large enough
    Vector<real_type> row_sums(cols->sparsity_pattern.n_rows() *
                               cols->chunk_size);
  
@@@ -1074,16 -1074,13 +1081,13 @@@ ChunkSparseMatrix<number>::residual (Ve
    dst = b;
  
    /////////
-   // the rest of this function is like
-   // vmult_add, except that we subtract
+   // the rest of this function is like vmult_add, except that we subtract
    // rather than add A*u
    /////////
 -  const unsigned int n_chunk_rows = cols->sparsity_pattern.n_rows();
 +  const size_type n_chunk_rows = cols->sparsity_pattern.n_rows();
  
-   // loop over all chunks. note that we need
-   // to treat the last chunk row and column
-   // differently if they have padding
-   // elements
+   // loop over all chunks. note that we need to treat the last chunk row and
+   // column differently if they have padding elements
    const bool rows_have_padding = (m() % cols->chunk_size != 0),
               cols_have_padding = (n() % cols->chunk_size != 0);
  
@@@ -1112,10 -1109,9 +1116,9 @@@
               u.begin() + *colnum_ptr * cols->chunk_size,
               dst_ptr);
            else
-             // we're at a chunk column that
-             // has padding
+             // we're at a chunk column that has padding
 -            for (unsigned int r=0; r<cols->chunk_size; ++r)
 -              for (unsigned int c=0; c<n() % cols->chunk_size; ++c)
 +            for (size_type r=0; r<cols->chunk_size; ++r)
 +              for (size_type c=0; c<n() % cols->chunk_size; ++c)
                  dst(chunk_row * cols->chunk_size + r)
                  -= (val_ptr[r*cols->chunk_size + c] *
                      u(*colnum_ptr * cols->chunk_size + c));
@@@ -1128,11 -1124,10 +1131,10 @@@
        dst_ptr += cols->chunk_size;
      }
  
-   // now deal with last chunk row if
-   // necessary
+   // now deal with last chunk row if necessary
    if (rows_have_padding)
      {
 -      const unsigned int chunk_row = n_chunk_rows - 1;
 +      const size_type chunk_row = n_chunk_rows - 1;
  
        const number *const val_end_of_row = &val[cols->sparsity_pattern.rowstart[chunk_row+1]
                                                  * cols->chunk_size
@@@ -1143,19 -1138,17 +1145,17 @@@
                ||
                (*colnum_ptr != cols->sparsity_pattern.n_cols()-1))
              {
-               // we're at a chunk row but not
-               // column that has padding
+               // we're at a chunk row but not column that has padding
 -              for (unsigned int r=0; r<m() % cols->chunk_size; ++r)
 -                for (unsigned int c=0; c<cols->chunk_size; ++c)
 +              for (size_type r=0; r<m() % cols->chunk_size; ++r)
 +                for (size_type c=0; c<cols->chunk_size; ++c)
                    dst(chunk_row * cols->chunk_size + r)
                    -= (val_ptr[r*cols->chunk_size + c] *
                        u(*colnum_ptr * cols->chunk_size + c));
              }
            else
-             // we're at a chunk row and
-             // column that has padding
+             // we're at a chunk row and column that has padding
 -            for (unsigned int r=0; r<m() % cols->chunk_size; ++r)
 -              for (unsigned int c=0; c<n() % cols->chunk_size; ++c)
 +            for (size_type r=0; r<m() % cols->chunk_size; ++r)
 +              for (size_type c=0; c<n() % cols->chunk_size; ++c)
                  dst(chunk_row * cols->chunk_size + r)
                  -= (val_ptr[r*cols->chunk_size + c] *
                      u(*colnum_ptr * cols->chunk_size + c));
@@@ -1477,18 -1458,17 +1465,17 @@@ void ChunkSparseMatrix<number>::print_p
    Assert (cols != 0, ExcNotInitialized());
    Assert (val != 0, ExcNotInitialized());
  
 -  const unsigned int chunk_size = cols->get_chunk_size();
 +  const size_type chunk_size = cols->get_chunk_size();
  
-   // loop over all chunk rows and columns,
-   // and each time we find something repeat
-   // it chunk_size times in both directions
+   // loop over all chunk rows and columns, and each time we find something
+   // repeat it chunk_size times in both directions
 -  for (unsigned int i=0; i<cols->sparsity_pattern.n_rows(); ++i)
 +  for (size_type i=0; i<cols->sparsity_pattern.n_rows(); ++i)
      {
 -      for (unsigned int d=0; d<chunk_size; ++d)
 -        for (unsigned int j=0; j<cols->sparsity_pattern.n_cols(); ++j)
 +      for (size_type d=0; d<chunk_size; ++d)
 +        for (size_type j=0; j<cols->sparsity_pattern.n_cols(); ++j)
            if (cols->sparsity_pattern(i,j) == SparsityPattern::invalid_entry)
              {
 -              for (unsigned int e=0; e<chunk_size; ++e)
 +              for (size_type e=0; e<chunk_size; ++e)
                  out << '.';
              }
            else if (std::fabs(val[cols->sparsity_pattern(i,j)]) > threshold)
diff --cc deal.II/include/deal.II/lac/chunk_sparsity_pattern.h
index 8c6618ba6f,200c80303d..73ba343023
--- a/deal.II/include/deal.II/lac/chunk_sparsity_pattern.h
+++ b/deal.II/include/deal.II/lac/chunk_sparsity_pattern.h
@@@ -46,46 -221,41 +221,45 @@@ namespace ChunkSparsityPatternIterator
  class ChunkSparsityPattern : public Subscriptor
  {
  public:
 +  /**
 +   * Declare the type for container size.
 +   */
 +  typedef types::global_dof_index size_type;
+   /**
+    * Typedef an iterator class that allows to walk over all nonzero elements
+    * of a sparsity pattern.
+    */
+   typedef ChunkSparsityPatternIterators::Iterator const_iterator;
+ 
+   /**
+    * Typedef an iterator class that allows to walk over all nonzero elements
+    * of a sparsity pattern.
+    *
+    * Since the iterator does not allow to modify the sparsity pattern, this
+    * type is the same as that for @p const_iterator.
+    */
+   typedef ChunkSparsityPatternIterators::Iterator iterator;
  
    /**
-    * Define a value which is used
-    * to indicate that a certain
-    * value in the colnums array
-    * is unused, i.e. does not
-    * represent a certain column
-    * number index.
+    * Define a value which is used to indicate that a certain value in the
+    * colnums array is unused, i.e. does not represent a certain column number
+    * index.
     *
-    * Indices with this invalid
-    * value are used to insert new
-    * entries to the sparsity
-    * pattern using the add() member
-    * function, and are removed when
+    * Indices with this invalid value are used to insert new entries to the
+    * sparsity pattern using the add() member function, and are removed when
     * calling compress().
     *
-    * You should not assume that the
-    * variable declared here has a
-    * certain value. The
-    * initialization is given here
-    * only to enable the compiler to
-    * perform some optimizations,
-    * but the actual value of the
-    * variable may change over time.
+    * You should not assume that the variable declared here has a certain
+    * value. The initialization is given here only to enable the compiler to
+    * perform some optimizations, but the actual value of the variable may
+    * change over time.
     */
 -  static const unsigned int invalid_entry = SparsityPattern::invalid_entry;
 +  static const size_type invalid_entry = SparsityPattern::invalid_entry;
  
    /**
-    * Initialize the matrix empty,
-    * that is with no memory
-    * allocated. This is useful if
-    * you want such objects as
-    * member variables in other
-    * classes. You can make the
-    * structure usable by calling
-    * the reinit() function.
+    * Initialize the matrix empty, that is with no memory allocated. This is
+    * useful if you want such objects as member variables in other classes. You
+    * can make the structure usable by calling the reinit() function.
     */
    ChunkSparsityPattern ();
  
@@@ -128,13 -281,12 +285,12 @@@
     *
     * @arg m number of rows
     * @arg n number of columns
-    * @arg max_per_row maximum
-    * number of nonzero entries per row
+    * @arg max_per_row maximum number of nonzero entries per row
     */
 -  ChunkSparsityPattern (const unsigned int m,
 -                        const unsigned int n,
 -                        const unsigned int max_chunks_per_row,
 -                        const unsigned int chunk_size);
 +  ChunkSparsityPattern (const size_type m,
 +                        const size_type n,
 +                        const size_type max_chunks_per_row,
 +                        const size_type chunk_size);
  
    /**
     * @deprecated This constructor is deprecated. Use the version
@@@ -152,16 -303,13 +307,13 @@@
     *
     * @arg m number of rows
     * @arg n number of columns
-    *
-    * @arg row_lengths possible
-    * number of nonzero entries for
-    * each row.  This vector must
-    * have one entry for each row.
+    * @arg row_lengths possible number of nonzero entries for each row.  This
+    * vector must have one entry for each row.
     */
 -  ChunkSparsityPattern (const unsigned int               m,
 -                        const unsigned int               n,
 -                        const std::vector<unsigned int> &row_lengths,
 -                        const unsigned int chunk_size);
 +  ChunkSparsityPattern (const size_type m,
 +                        const size_type n,
 +                        const std::vector<size_type> &row_lengths,
 +                        const size_type chunk_size);
  
    /**
     * @deprecated This constructor is deprecated. Use the version
@@@ -174,35 -322,27 +326,27 @@@
                          const bool optimize_diagonal) DEAL_II_DEPRECATED;
  
    /**
-    * Initialize a quadratic matrix
-    * of dimension <tt>n</tt> with
-    * at most <tt>max_per_row</tt>
-    * nonzero entries per row.
+    * Initialize a quadratic matrix of dimension <tt>n</tt> with at most
+    * <tt>max_per_row</tt> nonzero entries per row.
     *
-    * This constructor automatically
-    * enables optimized storage of
-    * diagonal elements. To avoid
-    * this, use the constructor
-    * taking row and column numbers
-    * separately.
+    * This constructor automatically enables optimized storage of diagonal
+    * elements. To avoid this, use the constructor taking row and column
+    * numbers separately.
     */
 -  ChunkSparsityPattern (const unsigned int n,
 -                        const unsigned int max_per_row,
 -                        const unsigned int chunk_size);
 +  ChunkSparsityPattern (const size_type n,
 +                        const size_type max_per_row,
 +                        const size_type chunk_size);
  
    /**
     * Initialize a quadratic matrix.
     *
     * @arg m number of rows and columns
-    *
-    * @arg row_lengths possible
-    * number of nonzero entries for
-    * each row.  This vector must
-    * have one entry for each row.
+    * @arg row_lengths possible number of nonzero entries for each row.  This
+    * vector must have one entry for each row.
     */
 -  ChunkSparsityPattern (const unsigned int               m,
 -                        const std::vector<unsigned int> &row_lengths,
 -                        const unsigned int               chunk_size);
 +  ChunkSparsityPattern (const size_type                m,
 +                        const std::vector<size_type>  &row_lengths,
 +                        const size_type                chunk_size);
  
    /**
     * @deprecated This constructor is deprecated. Use the version
@@@ -229,20 -366,17 +370,17 @@@
    ChunkSparsityPattern &operator = (const ChunkSparsityPattern &);
  
    /**
-    * Reallocate memory and set up data
-    * structures for a new matrix with
-    * <tt>m </tt>rows and <tt>n</tt> columns,
-    * with at most <tt>max_per_row</tt>
+    * Reallocate memory and set up data structures for a new matrix with <tt>m
+    * </tt>rows and <tt>n</tt> columns, with at most <tt>max_per_row</tt>
     * nonzero entries per row.
     *
-    * This function simply maps its
-    * operations to the other
-    * <tt>reinit</tt> function.
+    * This function simply maps its operations to the other <tt>reinit</tt>
+    * function.
     */
 -  void reinit (const unsigned int m,
 -               const unsigned int n,
 -               const unsigned int max_per_row,
 -               const unsigned int chunk_size);
 +  void reinit (const size_type m,
 +               const size_type n,
 +               const size_type max_per_row,
 +               const size_type chunk_size);
  
    /**
     * @deprecated This function is deprecated. Use the function
@@@ -287,20 -414,19 +418,19 @@@
     * @deprecated This function is deprecated. Use the function
     * without the last argument
     */
 -  void reinit (const unsigned int               m,
 -               const unsigned int               n,
 -               const std::vector<unsigned int> &row_lengths,
 -               const unsigned int chunk_size,
 -               const bool optimize_diagonal) DEAL_II_DEPRECATED;
 +  void reinit (const size_type                m,
 +               const size_type                n,
 +               const std::vector<size_type > &row_lengths,
 +               const size_type                chunk_size,
 +               const bool                     optimize_diagonal) DEAL_II_DEPRECATED;
  
    /**
-    * Same as above, but with a
-    * VectorSlice argument instead.
+    * Same as above, but with a VectorSlice argument instead.
     */
 -  void reinit (const unsigned int               m,
 -               const unsigned int               n,
 -               const VectorSlice<const std::vector<unsigned int> > &row_lengths,
 -               const unsigned int chunk_size);
 +  void reinit (const size_type m,
 +               const size_type n,
 +               const VectorSlice<const std::vector<size_type> > &row_lengths,
 +               const size_type chunk_size);
  
    /**
     * @deprecated This function is deprecated. Use the function
@@@ -335,84 -453,41 +457,41 @@@
    void compress ();
  
    /**
-    * This function can be used as a
-    * replacement for reinit(),
-    * subsequent calls to add() and
-    * a final call to close() if you
-    * know exactly in advance the
-    * entries that will form the
-    * matrix sparsity pattern.
-    *
-    * The first two parameters
-    * determine the size of the
-    * matrix. For the two last ones,
-    * note that a sparse matrix can
-    * be described by a sequence of
-    * rows, each of which is
-    * represented by a sequence of
-    * pairs of column indices and
-    * values. In the present
-    * context, the begin() and
-    * end() parameters designate
-    * iterators (of forward iterator
-    * type) into a container, one
-    * representing one row. The
-    * distance between begin()
-    * and end() should therefore
-    * be equal to
-    * n_rows(). These iterators
-    * may be iterators of
-    * <tt>std::vector</tt>,
-    * <tt>std::list</tt>, pointers into a
-    * C-style array, or any other
-    * iterator satisfying the
-    * requirements of a forward
-    * iterator. The objects pointed
-    * to by these iterators
-    * (i.e. what we get after
-    * applying <tt>operator*</tt> or
-    * <tt>operator-></tt> to one of these
-    * iterators) must be a container
-    * itself that provides functions
-    * <tt>begin</tt> and <tt>end</tt>
-    * designating a range of
-    * iterators that describe the
-    * contents of one
-    * line. Dereferencing these
-    * inner iterators must either
-    * yield a pair of a size_type 
-    * as column index and a
-    * value of arbitrary type (such
-    * a type would be used if we
-    * wanted to describe a sparse
-    * matrix with one such object),
-    * or simply a size_type
-    * (of we only wanted to describe
-    * a sparsity pattern). The
-    * function is able to determine
-    * itself whether an unsigned
-    * integer or a pair is what we
-    * get after dereferencing the
-    * inner iterators, through some
-    * template magic.
-    *
-    * While the order of the outer
-    * iterators denotes the
-    * different rows of the matrix,
-    * the order of the inner
-    * iterator denoting the columns
-    * does not matter, as they are
-    * sorted internal to this
-    * function anyway.
-    *
-    * Since that all sounds very
-    * complicated, consider the
-    * following example code, which
-    * may be used to fill a sparsity
-    * pattern:
+    * This function can be used as a replacement for reinit(), subsequent calls
+    * to add() and a final call to close() if you know exactly in advance the
+    * entries that will form the matrix sparsity pattern.
+    *
+    * The first two parameters determine the size of the matrix. For the two
+    * last ones, note that a sparse matrix can be described by a sequence of
+    * rows, each of which is represented by a sequence of pairs of column
+    * indices and values. In the present context, the begin() and end()
+    * parameters designate iterators (of forward iterator type) into a
+    * container, one representing one row. The distance between begin() and
+    * end() should therefore be equal to n_rows(). These iterators may be
+    * iterators of <tt>std::vector</tt>, <tt>std::list</tt>, pointers into a
+    * C-style array, or any other iterator satisfying the requirements of a
+    * forward iterator. The objects pointed to by these iterators (i.e. what we
+    * get after applying <tt>operator*</tt> or <tt>operator-></tt> to one of
+    * these iterators) must be a container itself that provides functions
+    * <tt>begin</tt> and <tt>end</tt> designating a range of iterators that
+    * describe the contents of one line. Dereferencing these inner iterators
+    * must either yield a pair of an unsigned integer as column index and a
+    * value of arbitrary type (such a type would be used if we wanted to
+    * describe a sparse matrix with one such object), or simply an unsigned
+    * integer (of we only wanted to describe a sparsity pattern). The function
+    * is able to determine itself whether an unsigned integer or a pair is what
+    * we get after dereferencing the inner iterators, through some template
+    * magic.
+    *
+    * While the order of the outer iterators denotes the different rows of the
+    * matrix, the order of the inner iterator denoting the columns does not
+    * matter, as they are sorted internal to this function anyway.
+    *
+    * Since that all sounds very complicated, consider the following example
+    * code, which may be used to fill a sparsity pattern:
     * @code
 -   * std::vector<std::vector<unsigned int> > column_indices (n_rows);
 -   * for (unsigned int row=0; row<n_rows; ++row)
 +   * std::vector<std::vector<size_type> > column_indices (n_rows);
 +   * for (size_type row=0; row<n_rows; ++row)
     *         // generate necessary columns in this row
     *   fill_row (column_indices[row]);
     *
@@@ -421,28 -496,18 +500,18 @@@
     *                     column_indices.end());
     * @endcode
     *
-    * Note that this example works
-    * since the iterators
-    * dereferenced yield containers
-    * with functions <tt>begin</tt> and
-    * <tt>end</tt> (namely
-    * <tt>std::vector</tt>s), and the
-    * inner iterators dereferenced
-    * yield size_type as
-    * column indices. Note that we
-    * could have replaced each of
-    * the two <tt>std::vector</tt>
-    * occurrences by <tt>std::list</tt>,
-    * and the inner one by
-    * <tt>std::set</tt> as well.
-    *
-    * Another example would be as
-    * follows, where we initialize a
-    * whole matrix, not only a
-    * sparsity pattern:
+    * Note that this example works since the iterators dereferenced yield
+    * containers with functions <tt>begin</tt> and <tt>end</tt> (namely
+    * <tt>std::vector</tt>s), and the inner iterators dereferenced yield
+    * unsigned integers as column indices. Note that we could have replaced
+    * each of the two <tt>std::vector</tt> occurrences by <tt>std::list</tt>,
+    * and the inner one by <tt>std::set</tt> as well.
+    *
+    * Another example would be as follows, where we initialize a whole matrix,
+    * not only a sparsity pattern:
     * @code
 -   * std::vector<std::map<unsigned int,double> > entries (n_rows);
 -   * for (unsigned int row=0; row<n_rows; ++row)
 +   * std::vector<std::map<size_type,double> > entries (n_rows);
 +   * for (size_type row=0; row<n_rows; ++row)
     *         // generate necessary pairs of columns
     *         // and corresponding values in this row
     *   fill_row (entries[row]);
@@@ -455,29 -520,20 +524,20 @@@
     *                   column_indices.end());
     * @endcode
     *
-    * This example works because
-    * dereferencing iterators of the
-    * inner type yields a pair of
-    * size_type and a value,
-    * the first of which we take as
-    * column index. As previously,
-    * the outer <tt>std::vector</tt>
-    * could be replaced by
-    * <tt>std::list</tt>, and the inner
-    * <tt>std::map<size_type,double></tt>
-    * could be replaced by
-    * <tt>std::vector<std::pair<size_type,double> ></tt>,
-    * or a list or set of such
-    * pairs, as they all return
-    * iterators that point to such
-    * pairs.
+    * This example works because dereferencing iterators of the inner type
+    * yields a pair of unsigned integers and a value, the first of which we
+    * take as column index. As previously, the outer <tt>std::vector</tt> could
+    * be replaced by <tt>std::list</tt>, and the inner <tt>std::map<unsigned
+    * int,double></tt> could be replaced by <tt>std::vector<std::pair<unsigned
+    * int,double> ></tt>, or a list or set of such pairs, as they all return
+    * iterators that point to such pairs.
     */
    template <typename ForwardIterator>
 -  void copy_from (const unsigned int    n_rows,
 -                  const unsigned int    n_cols,
 +  void copy_from (const size_type n_rows,
 +                  const size_type n_cols,
                    const ForwardIterator begin,
                    const ForwardIterator end,
 -                  const unsigned int chunk_size);
 +                  const size_type chunk_size);
  
     /**
      * @deprecated This function is deprecated. Use the function
@@@ -546,101 -617,78 +621,78 @@@
    bool empty () const;
  
    /**
-    * Return the chunk size given as
-    * argument when constructing this
-    * object.
+    * Return the chunk size given as argument when constructing this object.
     */
 -  unsigned int get_chunk_size () const;
 +  size_type get_chunk_size () const;
  
    /**
-    * Return the maximum number of entries per
-    * row. Before compression, this equals the
-    * number given to the constructor, while
-    * after compression, it equals the maximum
-    * number of entries actually allocated by
-    * the user.
+    * Return the maximum number of entries per row. Before compression, this
+    * equals the number given to the constructor, while after compression, it
+    * equals the maximum number of entries actually allocated by the user.
     */
 -  unsigned int max_entries_per_row () const;
 +  size_type max_entries_per_row () const;
  
    /**
-    * Add a nonzero entry to the matrix.
-    * This function may only be called
-    * for non-compressed sparsity patterns.
+    * Add a nonzero entry to the matrix. This function may only be called for
+    * non-compressed sparsity patterns.
     *
-    * If the entry already exists, nothing
-    * bad happens.
+    * If the entry already exists, nothing bad happens.
     */
 -  void add (const unsigned int i,
 -            const unsigned int j);
 +  void add (const size_type i,
 +            const size_type j);
  
    /**
-    * Make the sparsity pattern
-    * symmetric by adding the
-    * sparsity pattern of the
+    * Make the sparsity pattern symmetric by adding the sparsity pattern of the
     * transpose object.
     *
-    * This function throws an
-    * exception if the sparsity
-    * pattern does not represent a
-    * quadratic matrix.
+    * This function throws an exception if the sparsity pattern does not
+    * represent a quadratic matrix.
     */
    void symmetrize ();
  
    /**
-    * Return number of rows of this
-    * matrix, which equals the dimension
-    * of the image space.
+    * Return number of rows of this matrix, which equals the dimension of the
+    * image space.
     */
 -  unsigned int n_rows () const;
 +  inline size_type n_rows () const;
  
    /**
-    * Return number of columns of this
-    * matrix, which equals the dimension
-    * of the range space.
+    * Return number of columns of this matrix, which equals the dimension of
+    * the range space.
     */
 -  unsigned int n_cols () const;
 +  inline size_type n_cols () const;
  
    /**
-    * Check if a value at a certain
-    * position may be non-zero.
+    * Check if a value at a certain position may be non-zero.
     */
 -  bool exists (const unsigned int i,
 -               const unsigned int j) const;
 +  bool exists (const size_type i,
 +               const size_type j) const;
  
    /**
     * Number of entries in a specific row.
     */
 -  unsigned int row_length (const unsigned int row) const;
 +  size_type row_length (const size_type row) const;
  
    /**
-    * Compute the bandwidth of the matrix
-    * represented by this structure. The
-    * bandwidth is the maximum of $|i-j|$
-    * for which the index pair $(i,j)$
-    * represents a nonzero entry of the
-    * matrix. Consequently, the maximum
-    * bandwidth a $n\times m$ matrix can
-    * have is $\max\{n-1,m-1\}$.
+    * Compute the bandwidth of the matrix represented by this structure. The
+    * bandwidth is the maximum of $|i-j|$ for which the index pair $(i,j)$
+    * represents a nonzero entry of the matrix. Consequently, the maximum
+    * bandwidth a $n\times m$ matrix can have is $\max\{n-1,m-1\}$.
     */
 -  unsigned int bandwidth () const;
 +  size_type bandwidth () const;
  
    /**
-    * Return the number of nonzero elements of
-    * this matrix. Actually, it returns the
-    * number of entries in the sparsity
-    * pattern; if any of the entries should
-    * happen to be zero, it is counted
-    * anyway.
+    * Return the number of nonzero elements of this matrix. Actually, it
+    * returns the number of entries in the sparsity pattern; if any of the
+    * entries should happen to be zero, it is counted anyway.
     *
-    * This function may only be called if the
-    * matrix struct is compressed. It does not
-    * make too much sense otherwise anyway.
+    * This function may only be called if the matrix struct is compressed. It
+    * does not make too much sense otherwise anyway.
     */
 -  unsigned int n_nonzero_elements () const;
 +  size_type n_nonzero_elements () const;
  
    /**
-    * Return whether the structure is
-    * compressed or not.
+    * Return whether the structure is compressed or not.
     */
    bool is_compressed () const;
  
@@@ -826,27 -869,23 +873,23 @@@
    //@}
  private:
    /**
-    * Number of rows that this sparsity
-    * structure shall represent.
+    * Number of rows that this sparsity structure shall represent.
     */
 -  unsigned int rows;
 +  size_type rows;
  
    /**
-    * Number of columns that this sparsity
-    * structure shall represent.
+    * Number of columns that this sparsity structure shall represent.
     */
 -  unsigned int cols;
 +  size_type cols;
  
    /**
     * The size of chunks.
     */
 -  unsigned int chunk_size;
 +  size_type chunk_size;
  
    /**
-    * The reduced sparsity pattern. We store
-    * only which chunks exist, with each
-    * chunk a block in the matrix of size
-    * chunk_size by chunk_size.
+    * The reduced sparsity pattern. We store only which chunks exist, with each
+    * chunk a block in the matrix of size chunk_size by chunk_size.
     */
    SparsityPattern sparsity_pattern;
  
@@@ -863,9 -906,312 +910,312 @@@
  
  #ifndef DOXYGEN
  
+ namespace ChunkSparsityPatternIterators
+ {
+   inline
+   Accessor::
+   Accessor (const ChunkSparsityPattern *sparsity_pattern,
+             const unsigned int          row)
+     :
+     sparsity_pattern(sparsity_pattern),
+     reduced_accessor(row==sparsity_pattern->n_rows() ? 
+                      *sparsity_pattern->sparsity_pattern.end() :
+                      *sparsity_pattern->sparsity_pattern.
+                      begin(row/sparsity_pattern->get_chunk_size())),
+     chunk_row (row==sparsity_pattern->n_rows() ? 0 :
+                row%sparsity_pattern->get_chunk_size()),
+     chunk_col (0)
+   {}
+ 
+ 
+ 
+   inline
+   Accessor::
+   Accessor (const ChunkSparsityPattern *sparsity_pattern)
+     :
+     sparsity_pattern(sparsity_pattern),
+     reduced_accessor(*sparsity_pattern->sparsity_pattern.end()),
+     chunk_row (0),
+     chunk_col (0)
+   {}
+ 
+ 
+ 
+   inline
+   bool
+   Accessor::is_valid_entry () const
+   {
+     return reduced_accessor.is_valid_entry()
+       &&
+       sparsity_pattern->get_chunk_size()*reduced_accessor.row()+chunk_row <
+       sparsity_pattern->n_rows()
+       &&
+       sparsity_pattern->get_chunk_size()*reduced_accessor.column()+chunk_col <
+       sparsity_pattern->n_cols();
+   }
+ 
+ 
+ 
+   inline
+   unsigned int
+   Accessor::row() const
+   {
+     Assert (is_valid_entry() == true, ExcInvalidIterator());
+ 
+     return sparsity_pattern->get_chunk_size()*reduced_accessor.row() +
+       chunk_row;
+   }
+ 
+ 
+ 
+   inline
+   unsigned int
+   Accessor::column() const
+   {
+     Assert (is_valid_entry() == true, ExcInvalidIterator());
+ 
+     return sparsity_pattern->get_chunk_size()*reduced_accessor.column() +
+       chunk_col;
+   }
+ 
+ 
+ 
+   inline
+   std::size_t
+   Accessor::reduced_index() const
+   {
+     Assert (is_valid_entry() == true, ExcInvalidIterator());
+ 
+     return reduced_accessor.index_within_sparsity;
+   }
+ 
+ 
+ 
+ 
+   inline
+   bool
+   Accessor::operator == (const Accessor &other) const
+   {
+     // no need to check for equality of sparsity patterns as this is done in
+     // the reduced case already and every ChunkSparsityPattern has its own
+     // reduced sparsity pattern
+     return (reduced_accessor == other.reduced_accessor &&
+             chunk_row == other.chunk_row &&
+             chunk_col == other.chunk_col);
+   }
+ 
+ 
+ 
+   inline
+   bool
+   Accessor::operator < (const Accessor &other) const
+   {
+     Assert (sparsity_pattern == other.sparsity_pattern,
+             ExcInternalError());
+ 
+     if (chunk_row != other.chunk_row)
+       {
+         if (reduced_accessor.index_within_sparsity ==
+             reduced_accessor.sparsity_pattern->n_nonzero_elements())
+           return false;
+         if (other.reduced_accessor.index_within_sparsity ==
+             reduced_accessor.sparsity_pattern->n_nonzero_elements())
+           return true;
+ 
+         const unsigned int
+           global_row = sparsity_pattern->get_chunk_size()*reduced_accessor.row()
+           +chunk_row,
+           other_global_row = sparsity_pattern->get_chunk_size()*
+           other.reduced_accessor.row()+other.chunk_row;
+         if (global_row < other_global_row)
+           return true;
+         else if (global_row > other_global_row)
+           return false;
+       }
+ 
+     return (reduced_accessor.index_within_sparsity <
+             other.reduced_accessor.index_within_sparsity ||
+             (reduced_accessor.index_within_sparsity ==
+              other.reduced_accessor.index_within_sparsity &&
+              chunk_col < other.chunk_col));
+   }
+ 
+ 
+   inline
+   void
+   Accessor::advance ()
+   {
+     const unsigned int chunk_size = sparsity_pattern->get_chunk_size();
+     Assert (chunk_row < chunk_size && chunk_col < chunk_size,
+             ExcIteratorPastEnd());
+     Assert (reduced_accessor.row() * chunk_size + chunk_row <
+             sparsity_pattern->n_rows()
+             &&
+             reduced_accessor.column() * chunk_size + chunk_col <
+             sparsity_pattern->n_cols(),
+             ExcIteratorPastEnd());
+     if (chunk_size == 1)
+       {
+         reduced_accessor.advance();
+         return;
+       }
+ 
+     ++chunk_col;
+ 
+     // end of chunk
+     if (chunk_col == chunk_size
+         ||
+         reduced_accessor.column() * chunk_size + chunk_col ==
+         sparsity_pattern->n_cols())
+       {
+         const unsigned int reduced_row = reduced_accessor.row();
+         // end of row
+         if (reduced_accessor.index_within_sparsity + 1 ==
+             reduced_accessor.sparsity_pattern->rowstart[reduced_row+1])
+           {
+             ++chunk_row;
+ 
+             chunk_col = 0;
+ 
+             // end of chunk rows or end of matrix
+             if (chunk_row == chunk_size ||
+                 (reduced_row * chunk_size + chunk_row ==
+                  sparsity_pattern->n_rows()))
+               {
+                 chunk_row = 0;
+                 reduced_accessor.advance();
+               }
+             // go back to the beginning of the same reduced row but with
+             // chunk_row increased by one
+             else
+               reduced_accessor.index_within_sparsity =
+                 reduced_accessor.sparsity_pattern->rowstart[reduced_row];
+           }
+         // advance within chunk
+         else
+           {
+             reduced_accessor.advance();
+             chunk_col = 0;
+           }
+       }
+   }
+ 
+ 
+ 
+   inline
+   Iterator::Iterator (const ChunkSparsityPattern *sparsity_pattern,
+                       const unsigned int          row)
+     :
+     accessor(sparsity_pattern, row)
+   {}
+ 
+ 
+ 
+   inline
+   Iterator &
+   Iterator::operator++ ()
+   {
+     accessor.advance ();
+     return *this;
+   }
+ 
+ 
+ 
+   inline
+   Iterator
+   Iterator::operator++ (int)
+   {
+     const Iterator iter = *this;
+     accessor.advance ();
+     return iter;
+   }
+ 
+ 
+ 
+   inline
+   const Accessor &
+   Iterator::operator* () const
+   {
+     return accessor;
+   }
+ 
+ 
+ 
+   inline
+   const Accessor *
+   Iterator::operator-> () const
+   {
+     return &accessor;
+   }
+ 
+ 
+   inline
+   bool
+   Iterator::operator == (const Iterator &other) const
+   {
+     return (accessor == other.accessor);
+   }
+ 
+ 
+ 
+   inline
+   bool
+   Iterator::operator != (const Iterator &other) const
+   {
+     return ! (accessor == other.accessor);
+   }
+ 
+ 
+   inline
+   bool
+   Iterator::operator < (const Iterator &other) const
+   {
+     return accessor < other.accessor;
+   }
+ 
+ }
+ 
+ 
+ 
+ inline
+ ChunkSparsityPattern::iterator
+ ChunkSparsityPattern::begin () const
+ {
+   return iterator(this, 0);
+ }
+ 
+ 
+ inline
+ ChunkSparsityPattern::iterator
+ ChunkSparsityPattern::end () const
+ {
+   return iterator(this, n_rows());
+ }
+ 
+ 
+ 
+ inline
+ ChunkSparsityPattern::iterator
+ ChunkSparsityPattern::begin (const unsigned int r) const
+ {
+   Assert (r<n_rows(), ExcIndexRange(r,0,n_rows()));
+   return iterator(this, r);
+ }
+ 
+ 
+ 
+ inline
+ ChunkSparsityPattern::iterator
+ ChunkSparsityPattern::end (const unsigned int r) const
+ {
+   Assert (r<n_rows(), ExcIndexRange(r,0,n_rows()))
+     return iterator(this, r+1);
+ }
+ 
+ 
  
  inline
 -unsigned int
 +ChunkSparsityPattern::size_type
  ChunkSparsityPattern::n_rows () const
  {
    return rows;
@@@ -908,12 -1255,13 +1259,13 @@@ ChunkSparsityPattern::optimize_diagona
  
  
  template <typename ForwardIterator>
+ inline
  void
 -ChunkSparsityPattern::copy_from (const unsigned int    n_rows,
 -                                 const unsigned int    n_cols,
 +ChunkSparsityPattern::copy_from (const size_type n_rows,
 +                                 const size_type n_cols,
                                   const ForwardIterator begin,
                                   const ForwardIterator end,
 -                                 const unsigned int    chunk_size,
 +                                 const size_type chunk_size,
                                   const bool)
  {
    copy_from (n_rows, n_cols, begin, end, chunk_size);
@@@ -923,27 -1271,21 +1275,21 @@@
  
  template <typename ForwardIterator>
  void
 -ChunkSparsityPattern::copy_from (const unsigned int    n_rows,
 -                                 const unsigned int    n_cols,
 +ChunkSparsityPattern::copy_from (const size_type       n_rows,
 +                                 const size_type       n_cols,
                                   const ForwardIterator begin,
                                   const ForwardIterator end,
 -                                 const unsigned int    chunk_size)
 +                                 const size_type       chunk_size)
  {
 -  Assert (static_cast<unsigned int>(std::distance (begin, end)) == n_rows,
 +  Assert (static_cast<size_type>(std::distance (begin, end)) == n_rows,
            ExcIteratorRange (std::distance (begin, end), n_rows));
  
-   // first determine row lengths for
-   // each row. if the matrix is
-   // quadratic, then we might have to
-   // add an additional entry for the
-   // diagonal, if that is not yet
-   // present. as we have to call
-   // compress anyway later on, don't
-   // bother to check whether that
-   // diagonal entry is in a certain
-   // row or not
+   // first determine row lengths for each row. if the matrix is quadratic,
+   // then we might have to add an additional entry for the diagonal, if that
+   // is not yet present. as we have to call compress anyway later on, don't
+   // bother to check whether that diagonal entry is in a certain row or not
    const bool is_square = (n_rows == n_cols);
 -  std::vector<unsigned int> row_lengths;
 +  std::vector<size_type> row_lengths;
    row_lengths.reserve(n_rows);
    for (ForwardIterator i=begin; i!=end; ++i)
      row_lengths.push_back (std::distance (i->begin(), i->end())
@@@ -951,17 -1293,16 +1297,16 @@@
                             (is_square ? 1 : 0));
    reinit (n_rows, n_cols, row_lengths, chunk_size);
  
-   // now enter all the elements into
-   // the matrix
+   // now enter all the elements into the matrix
 -  unsigned int row = 0;
 +  size_type row = 0;
    typedef typename std::iterator_traits<ForwardIterator>::value_type::const_iterator inner_iterator;
    for (ForwardIterator i=begin; i!=end; ++i, ++row)
      {
        const inner_iterator end_of_row = i->end();
        for (inner_iterator j=i->begin(); j!=end_of_row; ++j)
          {
 -          const unsigned int col
 +          const size_type col
-             = internal::SparsityPatternTools::get_column_index_from_iterator(*j);
+           = internal::SparsityPatternTools::get_column_index_from_iterator(*j);
            Assert (col < n_cols, ExcInvalidIndex(col,n_cols));
  
            add (row, col);
diff --cc deal.II/include/deal.II/lac/parallel_vector.h
index 4c2dde075a,64cade9c6f..ac46606f6d
--- a/deal.II/include/deal.II/lac/parallel_vector.h
+++ b/deal.II/include/deal.II/lac/parallel_vector.h
@@@ -523,14 -523,14 +523,14 @@@ namespace paralle
         * are owned locally and for indices not
         * present at all.
         */
 -      bool is_ghost_entry (const types::global_dof_index global_index) const;
 +      bool is_ghost_entry (const size_type global_index) const;
  
        /**
-        * Make the @p Vector class a bit like
-        * the <tt>vector<></tt> class of the C++
-        * standard library by returning
-        * iterators to the start and end of the
-        * locally owned elements of this vector.
+        * Make the @p Vector class a bit like the <tt>vector<></tt> class of
+        * the C++ standard library by returning iterators to the start and end
+        * of the <i>locally owned</i> elements of this vector.
+        *
+        * It holds that end() - begin() == local_size().
         */
        iterator begin ();
  
diff --cc deal.II/include/deal.II/lac/trilinos_block_sparse_matrix.h
index 5f295a977a,1550c8a246..f9e8eb4f04
--- a/deal.II/include/deal.II/lac/trilinos_block_sparse_matrix.h
+++ b/deal.II/include/deal.II/lac/trilinos_block_sparse_matrix.h
@@@ -289,164 -290,26 +290,26 @@@ namespace TrilinosWrapper
       * elements of this
       * matrix.
       */
 -    unsigned int n_nonzero_elements () const;
 +    size_type n_nonzero_elements () const;
  
      /**
-      * Matrix-vector multiplication:
-      * let $dst = M*src$ with $M$
-      * being this matrix.
+      * Matrix-vector multiplication: let $dst = M*src$ with $M$ being this
+      * matrix. The vector types can be block vectors or non-block vectors
+      * (only if the matrix has only one row or column, respectively), and need
+      * to define TrilinosWrappers::SparseMatrix::vmult.
       */
-     void vmult (MPI::BlockVector       &dst,
-                 const MPI::BlockVector &src) const;
- 
- 
-     /**
-      * Matrix-vector multiplication:
-      * let $dst = M*src$ with $M$
-      * being this matrix, now applied
-      * to localized block vectors
-      * (works only when run on one
-      * processor).
-      */
-     void vmult (BlockVector       &dst,
-                 const BlockVector &src) const;
- 
-     /**
-      * Matrix-vector
-      * multiplication. Just like the
-      * previous function, but only
-      * applicable if the matrix has
-      * only one block column.
-      */
-     void vmult (MPI::BlockVector  &dst,
-                 const MPI::Vector &src) const;
- 
-     /**
-      * Matrix-vector
-      * multiplication. Just like the
-      * previous function, but only
-      * applicable if the matrix has
-      * only one block column, now
-      * applied to localized vectors
-      * (works only when run on one
-      * processor).
-      */
-     void vmult (BlockVector  &dst,
-                 const Vector &src) const;
- 
-     /**
-      * Matrix-vector
-      * multiplication. Just like the
-      * previous function, but only
-      * applicable if the matrix has
-      * only one block row.
-      */
-     void vmult (MPI::Vector            &dst,
-                 const MPI::BlockVector &src) const;
- 
-     /**
-      * Matrix-vector
-      * multiplication. Just like the
-      * previous function, but only
-      * applicable if the matrix has
-      * only one block row, now
-      * applied to localized vectors
-      * (works only when run on one
-      * processor).
-      */
-     void vmult (Vector            &dst,
-                 const BlockVector &src) const;
- 
-     /**
-      * Matrix-vector
-      * multiplication. Just like the
-      * previous function, but only
-      * applicable if the matrix has
-      * only one block.
-      */
-     void vmult (VectorBase       &dst,
-                 const VectorBase &src) const;
- 
-     /**
-      * Matrix-vector multiplication:
-      * let $dst = M^T*src$ with $M$
-      * being this matrix. This
-      * function does the same as
-      * vmult() but takes the
-      * transposed matrix.
-      */
-     void Tvmult (MPI::BlockVector       &dst,
-                  const MPI::BlockVector &src) const;
- 
-     /**
-      * Matrix-vector multiplication:
-      * let $dst = M^T*src$ with $M$
-      * being this matrix. This
-      * function does the same as
-      * vmult() but takes the
-      * transposed matrix, now applied
-      * to localized Trilinos vectors
-      * (works only when run on one
-      * processor).
-      */
-     void Tvmult (BlockVector       &dst,
-                  const BlockVector &src) const;
- 
-     /**
-      * Matrix-vector
-      * multiplication. Just like the
-      * previous function, but only
-      * applicable if the matrix has
-      * only one block row.
-      */
-     void Tvmult (MPI::BlockVector  &dst,
-                  const MPI::Vector &src) const;
+     template <typename VectorType1, typename VectorType2>
+     void vmult (VectorType1       &dst,
+                 const VectorType2 &src) const;
  
      /**
-      * Matrix-vector
-      * multiplication. Just like the
-      * previous function, but only
-      * applicable if the matrix has
-      * only one block row, now
-      * applied to localized Trilinos
-      * vectors (works only when run
-      * on one processor).
-      */
-     void Tvmult (BlockVector  &dst,
-                  const Vector &src) const;
- 
-     /**
-      * Matrix-vector
-      * multiplication. Just like the
-      * previous function, but only
-      * applicable if the matrix has
-      * only one block column.
-      */
-     void Tvmult (MPI::Vector    &dst,
-                  const MPI::BlockVector &src) const;
- 
-     /**
-      * Matrix-vector
-      * multiplication. Just like the
-      * previous function, but only
-      * applicable if the matrix has
-      * only one block column, now
-      * applied to localized Trilinos
-      * vectors (works only when run
-      * on one processor).
-      */
-     void Tvmult (Vector    &dst,
-                  const BlockVector &src) const;
- 
-     /**
-      * Matrix-vector
-      * multiplication. Just like the
-      * previous function, but only
-      * applicable if the matrix has
-      * only one block.
+      * Matrix-vector multiplication: let $dst = M^T*src$ with $M$ being this
+      * matrix. This function does the same as vmult() but takes the transposed
+      * matrix.
       */
-     void Tvmult (VectorBase       &dst,
-                  const VectorBase &src) const;
+     template <typename VectorType1, typename VectorType2>
+     void Tvmult (VectorType1       &dst,
+                  const VectorType2 &src) const;
  
      /**
       * Compute the residual of an
diff --cc deal.II/include/deal.II/lac/trilinos_sparse_matrix.h
index 648d4b3cbb,ff56ed065c..9cd050798d
--- a/deal.II/include/deal.II/lac/trilinos_sparse_matrix.h
+++ b/deal.II/include/deal.II/lac/trilinos_sparse_matrix.h
@@@ -3206,16 -3096,23 +3116,23 @@@ namespace TrilinosWrapper
  
      last_action = Insert;
  
 -    int *col_index_ptr;
 +    TrilinosWrappers::types::int_type *col_index_ptr;
-     TrilinosScalar const *col_value_ptr;
+     TrilinosScalar *col_value_ptr;
 -    int n_columns;
 +    TrilinosWrappers::types::int_type n_columns;
  
-     // If we don't elide zeros, the pointers
-     // are already available...
+     TrilinosScalar short_val_array[100];
+     int short_index_array[100];
+     std::vector<TrilinosScalar> long_val_array;
+     std::vector<int> long_index_array;
+ 
+ 
+     // If we don't elide zeros, the pointers are already available... need to
+     // cast to non-const pointers as that is the format taken by Trilinos (but
+     // we will not modify const data)
      if (elide_zero_values == false)
        {
 -        col_index_ptr = (int *)col_indices;
 +        col_index_ptr = (TrilinosWrappers::types::int_type *)col_indices;
-         col_value_ptr = values;
+         col_value_ptr = const_cast<TrilinosScalar*>(values);
          n_columns = n_cols;
        }
      else
@@@ -3242,41 -3145,29 +3165,30 @@@
                }
            }
  
 -        Assert(n_columns <= (int)n_cols, ExcInternalError());
 +        Assert(n_columns <= (TrilinosWrappers::types::int_type)n_cols, ExcInternalError());
- 
-         col_index_ptr = (TrilinosWrappers::types::int_type *)&column_indices[0];
-         col_value_ptr = &column_values[0];
        }
  
  
-     // If the calling matrix owns the row to
-     // which we want to insert values, we
-     // can directly call the Epetra_CrsMatrix
-     // input function, which is much faster
-     // than the Epetra_FECrsMatrix
-     // function. We distinguish between two
-     // cases: the first one is when the matrix
-     // is not filled (i.e., it is possible to
-     // add new elements to the sparsity pattern),
-     // and the second one is when the pattern is
-     // already fixed. In the former case, we
-     // add the possibility to insert new values,
-     // and in the second we just replace
+     // If the calling matrix owns the row to which we want to insert values,
+     // we can directly call the Epetra_CrsMatrix input function, which is much
+     // faster than the Epetra_FECrsMatrix function. We distinguish between two
+     // cases: the first one is when the matrix is not filled (i.e., it is
+     // possible to add new elements to the sparsity pattern), and the second
+     // one is when the pattern is already fixed. In the former case, we add
+     // the possibility to insert new values, and in the second we just replace
      // data.
 -    if (row_partitioner().MyGID(static_cast<int>(row)) == true)
 +    if (row_partitioner().MyGID(static_cast<TrilinosWrappers::types::int_type>(row)) == true)
        {
          if (matrix->Filled() == false)
            {
 -            ierr = matrix->Epetra_CrsMatrix::InsertGlobalValues(row, n_columns,
 -                                                                col_value_ptr,
 -                                                                col_index_ptr);
 +            ierr = matrix->Epetra_CrsMatrix::InsertGlobalValues(
 +                static_cast<TrilinosWrappers::types::int_type>(row), 
 +                static_cast<int>(n_columns),const_cast<double *>(col_value_ptr), 
 +                col_index_ptr);
  
-             // When inserting elements, we do
-             // not want to create exceptions in
-             // the case when inserting non-local
-             // data (since that's what we want
-             // to do right now).
+             // When inserting elements, we do not want to create exceptions in
+             // the case when inserting non-local data (since that's what we
+             // want to do right now).
              if (ierr > 0)
                ierr = 0;
            }
@@@ -3438,19 -3319,25 +3342,25 @@@
  
      last_action = Add;
  
 -    int *col_index_ptr;
 +    TrilinosWrappers::types::int_type *col_index_ptr;
-     TrilinosScalar const *col_value_ptr;
+     TrilinosScalar *col_value_ptr;
 -    int n_columns;
 +    TrilinosWrappers::types::int_type n_columns;
  
-     // If we don't elide zeros, the pointers
-     // are already available...
+     double short_val_array[100];
+     int short_index_array[100];
+     std::vector<TrilinosScalar> long_val_array;
+     std::vector<int> long_index_array;
+ 
+     // If we don't elide zeros, the pointers are already available... need to
+     // cast to non-const pointers as that is the format taken by Trilinos (but
+     // we will not modify const data)
      if (elide_zero_values == false)
        {
 -        col_index_ptr = (int *)col_indices;
 +        col_index_ptr = (TrilinosWrappers::types::int_type *)col_indices;
-         col_value_ptr = values;
+         col_value_ptr = const_cast<TrilinosScalar*>(values);
          n_columns = n_cols;
  #ifdef DEBUG
 -        for (unsigned int j=0; j<n_cols; ++j)
 +        for (size_type j=0; j<n_cols; ++j)
            Assert (numbers::is_finite(values[j]), ExcNumberNotFinite());
  #endif
        }
@@@ -3465,9 -3359,10 +3382,10 @@@
            }
  
          n_columns = 0;
 -        for (unsigned int j=0; j<n_cols; ++j)
 +        for (size_type j=0; j<n_cols; ++j)
            {
              const double value = values[j];
+ 
              Assert (numbers::is_finite(value), ExcNumberNotFinite());
              if (value != 0)
                {
@@@ -3477,39 -3372,30 +3395,31 @@@
                }
            }
  
 -        Assert(n_columns <= (int)n_cols, ExcInternalError());
 +        Assert(n_columns <= (TrilinosWrappers::types::int_type)n_cols, ExcInternalError());
  
-         col_index_ptr = (TrilinosWrappers::types::int_type *)&column_indices[0];
-         col_value_ptr = &column_values[0];
        }
  
-     // If the calling matrix owns the row to
-     // which we want to add values, we
-     // can directly call the Epetra_CrsMatrix
-     // input function, which is much faster
-     // than the Epetra_FECrsMatrix function.
+     // If the calling matrix owns the row to which we want to add values, we
+     // can directly call the Epetra_CrsMatrix input function, which is much
+     // faster than the Epetra_FECrsMatrix function.
 -    if (row_partitioner().MyGID(static_cast<int>(row)) == true)
 +    if (row_partitioner().MyGID(static_cast<TrilinosWrappers::types::int_type>(row)) == true)
        {
          ierr = matrix->Epetra_CrsMatrix::SumIntoGlobalValues(row, n_columns,
-                                                              const_cast<double *>(col_value_ptr),
+                                                              col_value_ptr,
                                                               col_index_ptr);
        }
      else
        {
-         // When we're at off-processor data, we
-         // have to stick with the standard
-         // SumIntoGlobalValues
-         // function. Nevertheless, the way we
-         // call it is the fastest one (any other
-         // will lead to repeated allocation and
-         // deallocation of memory in order to
-         // call the function we already use,
-         // which is very unefficient if writing
-         // one element at a time).
+         // When we're at off-processor data, we have to stick with the
+         // standard SumIntoGlobalValues function. Nevertheless, the way we
+         // call it is the fastest one (any other will lead to repeated
+         // allocation and deallocation of memory in order to call the function
+         // we already use, which is very inefficient if writing one element at
+         // a time).
          compressed = false;
  
 -        ierr = matrix->SumIntoGlobalValues (1, (int *)&row, n_columns,
 +        ierr = matrix->SumIntoGlobalValues (1,
 +                                            (TrilinosWrappers::types::int_type *)&row, n_columns,
                                              col_index_ptr,
                                              &col_value_ptr,
                                              Epetra_FECrsMatrix::ROW_MAJOR);
@@@ -3789,16 -3642,11 +3683,11 @@@
      Assert (&src != &dst, ExcSourceEqualsDestination());
      Assert (matrix->Filled(), ExcMatrixNotCompressed());
  
- #ifndef DEAL_II_USE_LARGE_INDEX_TYPE
-     const int n_global_elements = matrix->DomainMap().NumGlobalElements();
- #else
-     const long long int n_global_elements = matrix->DomainMap().NumGlobalElements64();
- #endif
- 
-     AssertDimension (static_cast<size_type>(matrix->DomainMap().NumMyElements()),
-                      static_cast<size_type>(n_global_elements));
-     AssertDimension (dst.size(), static_cast<size_type>(matrix->RangeMap().NumMyElements()));
-     AssertDimension (src.size(), static_cast<size_type>(matrix->DomainMap().NumMyElements()));
+     internal::SparseMatrix::check_vector_map_equality(*matrix, src, dst);
 -    const int dst_local_size = dst.end() - dst.begin();
 -    AssertDimension (dst_local_size, matrix->RangeMap().NumMyElements());
 -    const int src_local_size = src.end() - src.begin();
 -    AssertDimension (src_local_size, matrix->DomainMap().NumMyElements());
++    const size_type dst_local_size = dst.end() - dst.begin();
++    AssertDimension (dst_local_size, static_cast<size_type>(matrix->RangeMap().NumMyElements()));
++    const size_type src_local_size = src.end() - src.begin();
++    AssertDimension (src_local_size, static_cast<size_type>(matrix->DomainMap().NumMyElements()));
  
      Epetra_Vector tril_dst (View, matrix->RangeMap(), dst.begin());
      Epetra_Vector tril_src (View, matrix->DomainMap(),
@@@ -3862,16 -3668,11 +3709,11 @@@
      Assert (&src != &dst, ExcSourceEqualsDestination());
      Assert (matrix->Filled(), ExcMatrixNotCompressed());
  
- #ifndef DEAL_II_USE_LARGE_INDEX_TYPE
-     const int n_global_elements = matrix->DomainMap().NumGlobalElements();
- #else
-     const long long int n_global_elements = matrix->DomainMap().NumGlobalElements64();
- #endif
- 
-     AssertDimension (static_cast<size_type>(matrix->DomainMap().NumMyElements()),
-                      static_cast<size_type>(n_global_elements));
-     AssertDimension (dst.size(), static_cast<size_type>(matrix->DomainMap().NumMyElements()));
-     AssertDimension (src.size(), static_cast<size_type>(matrix->RangeMap().NumMyElements()));
+     internal::SparseMatrix::check_vector_map_equality(*matrix, dst, src);
 -    const int dst_local_size = dst.end() - dst.begin();
 -    AssertDimension (dst_local_size, matrix->DomainMap().NumMyElements());
 -    const int src_local_size = src.end() - src.begin();
 -    AssertDimension (src_local_size, matrix->RangeMap().NumMyElements());
++    const size_type dst_local_size = dst.end() - dst.begin();
++    AssertDimension (dst_local_size, static_cast<size_type>(matrix->DomainMap().NumMyElements()));
++    const size_type src_local_size = src.end() - src.begin();
++    AssertDimension (src_local_size, static_cast<size_type>(matrix->RangeMap().NumMyElements()));
  
      Epetra_Vector tril_dst (View, matrix->DomainMap(), dst.begin());
      Epetra_Vector tril_src (View, matrix->RangeMap(),
diff --cc deal.II/include/deal.II/lac/trilinos_vector_base.h
index 2740b3c509,f9c2145a90..aed6e15110
--- a/deal.II/include/deal.II/lac/trilinos_vector_base.h
+++ b/deal.II/include/deal.II/lac/trilinos_vector_base.h
@@@ -272,10 -267,12 +272,12 @@@ namespace TrilinosWrapper
       * <tt>C</tt> standard libraries
       * <tt>vector<...></tt> class.
       */
-     typedef TrilinosScalar            value_type;
-     typedef TrilinosScalar            real_type;
-     typedef dealii::types::global_dof_index   size_type;
-     typedef internal::VectorReference reference;
+     typedef TrilinosScalar                  value_type;
+     typedef TrilinosScalar                  real_type;
 -    typedef std::size_t                     size_type;
++    typedef dealii::types::global_dof_index size_type;
+     typedef value_type                     *iterator;
+     typedef const value_type               *const_iterator;
+     typedef internal::VectorReference       reference;
      typedef const internal::VectorReference const_reference;
  
      /**
@@@ -662,8 -659,35 +664,35 @@@
       * elements sits on another
       * process.
       */
 -    TrilinosScalar el (const unsigned int index) const;
 +    TrilinosScalar el (const size_type index) const;
  
+     /**
+      * Make the Vector class a bit like the <tt>vector<></tt> class of
+      * the C++ standard library by returning iterators to the start and end
+      * of the locally owned elements of this vector. The ordering of local elements corresponds to the one given 
+      *
+      * It holds that end() - begin() == local_size().
+      */
+     iterator begin ();
+ 
+     /**
+      * Return constant iterator to the start of the locally owned elements
+      * of the vector.
+      */
+     const_iterator begin () const;
+ 
+     /**
+      * Return an iterator pointing to the element past the end of the array
+      * of locally owned entries.
+      */
+     iterator end ();
+ 
+     /**
+      * Return a constant iterator pointing to the element past the end of
+      * the array of the locally owned entries.
+      */
+     const_iterator end () const;
+ 
      /**
       * A collective set operation:
       * instead of setting individual
diff --cc deal.II/include/deal.II/numerics/vector_tools.templates.h
index 29d8821b42,3dec6628b3..575a416a2e
--- a/deal.II/include/deal.II/numerics/vector_tools.templates.h
+++ b/deal.II/include/deal.II/numerics/vector_tools.templates.h
@@@ -4267,11 -4278,10 +4278,10 @@@ namespace VectorTool
                          "to imposing Dirichlet values on the vector-valued "
                          "quantity."));
  
 -    std::vector<unsigned int> face_dofs;
 +    std::vector<types::global_dof_index> face_dofs;
  
-     // create FE and mapping
-     // collections for all elements in
-     // use by this DoFHandler
+     // create FE and mapping collections for all elements in use by this
+     // DoFHandler
      hp::FECollection<dim,spacedim>      fe_collection (dof_handler.get_fe());
      hp::MappingCollection<dim,spacedim> mapping_collection;
      for (unsigned int i=0; i<fe_collection.size(); ++i)
diff --cc deal.II/source/lac/chunk_sparsity_pattern.cc
index 2c5a9442be,5bcb7c0470..498568f4b1
--- a/deal.II/source/lac/chunk_sparsity_pattern.cc
+++ b/deal.II/source/lac/chunk_sparsity_pattern.cc
@@@ -166,9 -165,8 +165,8 @@@ ChunkSparsityPattern::reinit (const siz
  {
    Assert (chunk_size > 0, ExcInvalidNumber (chunk_size));
  
-   // simply map this function to the
-   // other @p{reinit} function
+   // simply map this function to the other @p{reinit} function
 -  const std::vector<unsigned int> row_lengths (m, max_per_row);
 +  const std::vector<size_type> row_lengths (m, max_per_row);
    reinit (m, n, row_lengths, chunk_size);
  }
  
@@@ -201,31 -199,30 +199,30 @@@ ChunkSparsityPattern::reinit 
  
    this->chunk_size = chunk_size;
  
-   // pass down to the necessary information
-   // to the underlying object. we need to
-   // calculate how many chunks we need: we
-   // need to round up (m/chunk_size) and
-   // (n/chunk_size). rounding up in integer
-   // arithmetic equals
+   // pass down to the necessary information to the underlying object. we need
+   // to calculate how many chunks we need: we need to round up (m/chunk_size)
+   // and (n/chunk_size). rounding up in integer arithmetic equals
    // ((m+chunk_size-1)/chunk_size):
 -  const unsigned int m_chunks = (m+chunk_size-1) / chunk_size,
 -                     n_chunks = (n+chunk_size-1) / chunk_size;
 +  const size_type m_chunks = (m+chunk_size-1) / chunk_size,
 +                  n_chunks = (n+chunk_size-1) / chunk_size;
  
-   // compute the maximum number of chunks in
-   // each row. the passed array denotes the
-   // number of entries in each row of the big
-   // matrix -- in the worst case, these are
-   // all in independent chunks, so we have to
-   // calculate it as follows (as an example:
-   // let chunk_size==2,
-   // row_lengths={2,2,...}, and entries in
-   // row zero at columns {0,2} and for row
-   // one at {4,6} --> we'll need 4 chunks for
-   // the first chunk row!) :
+   // compute the maximum number of chunks in each row. the passed array
+   // denotes the number of entries in each row of the big matrix -- in the
+   // worst case, these are all in independent chunks, so we have to calculate
+   // it as follows (as an example: let chunk_size==2, row_lengths={2,2,...},
+   // and entries in row zero at columns {0,2} and for row one at {4,6} -->
+   // we'll need 4 chunks for the first chunk row!) :
    std::vector<unsigned int> chunk_row_lengths (m_chunks, 0);
 -  for (unsigned int i=0; i<m; ++i)
 +  for (size_type i=0; i<m; ++i)
      chunk_row_lengths[i/chunk_size] += row_lengths[i];
  
+   // for the case that the reduced sparsity pattern optimizes the diagonal but
+   // the actual sparsity pattern does not, need to take one more entry in the
+   // row to fit the user-required entry
+   if (m != n && m_chunks == n_chunks)
+     for (unsigned int i=0; i<m_chunks; ++i)
+       ++chunk_row_lengths[i];
+ 
    sparsity_pattern.reinit (m_chunks,
                             n_chunks,
                             chunk_row_lengths);
@@@ -292,11 -284,10 +289,10 @@@ ChunkSparsityPattern::copy_from (const 
  {
    Assert (chunk_size > 0, ExcInvalidNumber (chunk_size));
  
-   // count number of entries per row, then
-   // initialize the underlying sparsity
+   // count number of entries per row, then initialize the underlying sparsity
    // pattern
 -  std::vector<unsigned int> entries_per_row (csp.n_rows(), 0);
 -  for (unsigned int row = 0; row<csp.n_rows(); ++row)
 +  std::vector<size_type> entries_per_row (csp.n_rows(), 0);
 +  for (size_type row = 0; row<csp.n_rows(); ++row)
      entries_per_row[row] = csp.row_length(row);
  
    reinit (csp.n_rows(), csp.n_cols(),
@@@ -437,17 -478,28 +483,28 @@@ ChunkSparsityPattern::symmetrize (
  
  
  
- void
- ChunkSparsityPattern::symmetrize ()
 -unsigned int
 -ChunkSparsityPattern::row_length (const unsigned int i) const
++ChunkSparsityPattern::size_type 
++ChunkSparsityPattern::row_length (const size_type i) const
  {
-   // matrix must be square. note that the for
-   // some matrix sizes, the current sparsity
-   // pattern may not be square even if the
-   // underlying sparsity pattern is (e.g. a
-   // 10x11 matrix with chunk_size 4)
-   Assert (rows==cols, ExcNotQuadratic());
+   Assert (i<rows, ExcIndexRange(i,0,rows));
  
-   sparsity_pattern.symmetrize ();
+   // find out if we did padding and if this row is affected by it
+   if (n_cols() % chunk_size == 0)
+     return sparsity_pattern.row_length (i/chunk_size) * chunk_size;
+   else
+     // if columns don't align, then just iterate over all chunks and see
+     // what this leads to
+     {
+       SparsityPattern::const_iterator p = sparsity_pattern.begin(i/chunk_size),
+                                       end = sparsity_pattern.end(i/chunk_size);
+       unsigned int n = 0;
+       for ( ; p != end; ++p)
+         if (p->column() != sparsity_pattern.n_cols() - 1)
+           n += chunk_size;
+         else
+           n += (n_cols() % chunk_size);
+       return n;
+     }
  }
  
  
@@@ -472,11 -522,10 +527,10 @@@ ChunkSparsityPattern::n_nonzero_element
            &&
            (n_cols() % chunk_size == 0))
          {
-           // columns align with chunks, but
-           // not rows
+           // columns align with chunks, but not rows
 -          unsigned int n = sparsity_pattern.n_nonzero_elements() *
 -                           chunk_size *
 -                           chunk_size;
 +          size_type n = sparsity_pattern.n_nonzero_elements() *
 +                        chunk_size *
 +                        chunk_size;
            n -= (sparsity_pattern.n_rows() * chunk_size - n_rows()) *
                 sparsity_pattern.row_length(sparsity_pattern.n_rows()-1) *
                 chunk_size;
@@@ -485,14 -534,13 +539,13 @@@
  
        else
          {
-     	  // if columns don't align, then
-     	  // just iterate over all chunks and
-     	  // see what this leads to. follow the advice in the documentation of
-     	  // the sparsity pattern iterators to do the loop over individual rows,
+     	  // if columns don't align, then just iterate over all chunks and see
+     	  // what this leads to. follow the advice in the documentation of the
+     	  // sparsity pattern iterators to do the loop over individual rows,
      	  // rather than all elements
 -    	  unsigned int n = 0;
 +    	  size_type n = 0;
  
 -    	  for (unsigned int row = 0; row < sparsity_pattern.n_rows(); ++row)
 +    	  for (size_type row = 0; row < sparsity_pattern.n_rows(); ++row)
      	    {
      	      SparsityPattern::const_iterator p = sparsity_pattern.begin(row);
      	      for (; p!=sparsity_pattern.end(row); ++p)
@@@ -566,26 -610,22 +615,22 @@@ ChunkSparsityPattern::print_gnuplot (st
  
    AssertThrow (out, ExcIO());
  
-   // for each entry in the underlying
-   // sparsity pattern, repeat everything
+   // for each entry in the underlying sparsity pattern, repeat everything
    // chunk_size x chunk_size times
 -  for (unsigned int i=0; i<sparsity_pattern.rows; ++i)
 -    for (unsigned int j=sparsity_pattern.rowstart[i];
 +  for (size_type i=0; i<sparsity_pattern.rows; ++i)
 +    for (size_type j=sparsity_pattern.rowstart[i];
           j<sparsity_pattern.rowstart[i+1]; ++j)
        if (sparsity_pattern.colnums[j] != sparsity_pattern.invalid_entry)
 -        for (unsigned int d=0;
 +        for (size_type d=0;
               ((d<chunk_size) &&
                (sparsity_pattern.colnums[j]*chunk_size+d < n_cols()));
               ++d)
 -          for (unsigned int e=0;
 +          for (size_type e=0;
                 (e<chunk_size) && (i*chunk_size + e < n_rows());
                 ++e)
-             // while matrix entries are
-             // usually written (i,j), with i
-             // vertical and j horizontal,
-             // gnuplot output is x-y, that is
-             // we have to exchange the order
-             // of output
+             // while matrix entries are usually written (i,j), with i vertical
+             // and j horizontal, gnuplot output is x-y, that is we have to
+             // exchange the order of output
              out << sparsity_pattern.colnums[j]*chunk_size+d << " "
                  << -static_cast<signed int>(i*chunk_size+e)
                  << std::endl;
@@@ -595,16 -635,13 +640,13 @@@
  
  
  
 -unsigned int
 +ChunkSparsityPattern::size_type 
  ChunkSparsityPattern::bandwidth () const
  {
-   // calculate the bandwidth from that of the
-   // underlying sparsity pattern. note that
-   // even if the bandwidth of that is zero,
-   // then the bandwidth of the chunky pattern
-   // is chunk_size-1, if it is 1 then the
-   // chunky pattern has
-   // chunk_size+(chunk_size-1), etc
+   // calculate the bandwidth from that of the underlying sparsity
+   // pattern. note that even if the bandwidth of that is zero, then the
+   // bandwidth of the chunky pattern is chunk_size-1, if it is 1 then the
+   // chunky pattern has chunk_size+(chunk_size-1), etc
    //
    // we'll cut it off at max(n(),m())
    return std::min (sparsity_pattern.bandwidth()*chunk_size
@@@ -686,8 -721,12 +726,12 @@@ ChunkSparsityPattern::memory_consumptio
  
  // explicit instantiations
  template
+ void ChunkSparsityPattern::copy_from<SparsityPattern> (const SparsityPattern &,
+                                                        const unsigned int,
+                                                        const bool);
+ template
  void ChunkSparsityPattern::copy_from<CompressedSparsityPattern> (const CompressedSparsityPattern &,
 -    const unsigned int,
 +    const size_type ,
      const bool);
  template
  void ChunkSparsityPattern::copy_from<CompressedSetSparsityPattern> (const CompressedSetSparsityPattern &,
@@@ -695,11 -734,39 +739,39 @@@
      const bool);
  template
  void ChunkSparsityPattern::copy_from<CompressedSimpleSparsityPattern> (const CompressedSimpleSparsityPattern &,
 -    const unsigned int,
 +    const size_type ,
      const bool);
  template
+ void ChunkSparsityPattern::create_from<SparsityPattern>
+ (const unsigned int,
+  const unsigned int,
+  const SparsityPattern &,
+  const unsigned int,
+  const bool);
+ template
+ void ChunkSparsityPattern::create_from<CompressedSparsityPattern>
+ (const unsigned int,
+  const unsigned int,
+  const CompressedSparsityPattern &,
+  const unsigned int,
+  const bool);
+ template
+ void ChunkSparsityPattern::create_from<CompressedSetSparsityPattern>
+ (const unsigned int,
+  const unsigned int,
+  const CompressedSetSparsityPattern &,
+  const unsigned int,
+  const bool);
+ template
+ void ChunkSparsityPattern::create_from<CompressedSimpleSparsityPattern>
+ (const unsigned int,
+  const unsigned int,
+  const CompressedSimpleSparsityPattern &,
+  const unsigned int,
+  const bool);
+ template
  void ChunkSparsityPattern::copy_from<float> (const FullMatrix<float> &,
 -                                             const unsigned int,
 +                                             const size_type ,
                                               const bool);
  template
  void ChunkSparsityPattern::copy_from<double> (const FullMatrix<double> &,
diff --cc deal.II/source/numerics/matrix_tools.cc
index e20a1334c4,77fa2af43d..e0bd90bbd2
--- a/deal.II/source/numerics/matrix_tools.cc
+++ b/deal.II/source/numerics/matrix_tools.cc
@@@ -206,9 -206,10 +206,10 @@@ namespace MatrixCreato
  
        struct CopyData
        {
 -        std::vector<unsigned int> dof_indices;
 +        std::vector<types::global_dof_index> dof_indices;
          FullMatrix<double>        cell_matrix;
          dealii::Vector<double>    cell_rhs;
+         const ConstraintMatrix   *constraints;
        };
      }
  
diff --cc tests/trilinos/slowness_04.cc
index 8dc2febce6,4048a0e734..9c3b472da3
--- a/tests/trilinos/slowness_04.cc
+++ b/tests/trilinos/slowness_04.cc
@@@ -93,8 -92,8 +93,8 @@@ void test (
              matrix.set(global, global+N, rand());
            }
        }
-   matrix.compress ();
+   matrix.compress (VectorOperation::insert);
 -  
 +
                                     // then do a single matrix-vector
                                     // multiplication with subsequent formation
                                     // of the matrix norm