From: kronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Tue, 21 Jan 2014 08:30:25 +0000 (+0000)
Subject: Test extract_row_copy, fix two bugs
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=9bfb7f475b066b102b6eba3a9703376dfc4b38d5;p=dealii-svn.git

Test extract_row_copy, fix two bugs

git-svn-id: https://svn.dealii.org/trunk@32276 0785d39b-7218-0410-832d-ea1e28bc413d
---

diff --git a/deal.II/include/deal.II/lac/chunk_sparse_matrix.templates.h b/deal.II/include/deal.II/lac/chunk_sparse_matrix.templates.h
index b493ed3f26..188a020e40 100644
--- a/deal.II/include/deal.II/lac/chunk_sparse_matrix.templates.h
+++ b/deal.II/include/deal.II/lac/chunk_sparse_matrix.templates.h
@@ -605,24 +605,46 @@ ChunkSparseMatrix<number>::extract_row_copy (const size_type row,
                                              size_type *column_indices,
                                              number *values) const
 {
-  row_length = cols->row_length(row);
+  AssertIndexRange(cols->row_length(row), array_length+1);
+  AssertIndexRange(row, m());
   const unsigned int chunk_size = cols->get_chunk_size();
   const size_type reduced_row = row/chunk_size;
-  AssertIndexRange(cols->sparsity_pattern.row_length(reduced_row)*chunk_size,
-                   array_length+1);
 
   SparsityPattern::iterator it = cols->sparsity_pattern.begin(reduced_row),
     itend = cols->sparsity_pattern.end(reduced_row);
   const number *val_ptr = &val[(it-cols->sparsity_pattern.begin(0))*chunk_size*chunk_size
                                +(row%chunk_size)*chunk_size];
-  for ( ; it < itend; ++it)
+
+  // find out if we did padding and if this row is affected by it
+  if (cols->n_cols() % chunk_size == 0)
+    {
+      for ( ; it < itend; ++it)
+        {
+          for (unsigned int c=0; c<chunk_size; ++c)
+            {
+              *values++ = val_ptr[c];
+              *column_indices++ = it->column()*chunk_size+c;
+            }
+          val_ptr += chunk_size*chunk_size;
+        }
+      row_length = chunk_size * (cols->sparsity_pattern.row_length(reduced_row));
+    }
+  else
     {
-      for (unsigned int c=0; c<chunk_size; ++c)
+      const unsigned int last_chunk_size = cols->n_cols() % chunk_size;
+      row_length = 0;
+      for ( ; it < itend; ++it)
         {
-          *values++ = val_ptr[c];
-          *column_indices++ = it->column()*chunk_size+c;
+          const unsigned int next_chunk_size =
+            (it->column()==cols->sparsity_pattern.n_cols()-1) ?
+            last_chunk_size : chunk_size;
+          for (unsigned int c=0; c<next_chunk_size; ++c, ++row_length)
+            {
+              *values++ = val_ptr[c];
+              *column_indices++ = it->column()*chunk_size+c;
+            }
+          val_ptr += chunk_size*chunk_size;
         }
-      val_ptr += chunk_size*chunk_size;
     }
 }
 
diff --git a/deal.II/source/lac/sparsity_pattern.cc b/deal.II/source/lac/sparsity_pattern.cc
index 617b820c16..63841efa32 100644
--- a/deal.II/source/lac/sparsity_pattern.cc
+++ b/deal.II/source/lac/sparsity_pattern.cc
@@ -190,42 +190,29 @@ SparsityPattern::SparsityPattern (const SparsityPattern &original,
 
   reinit (original.rows, original.cols, max_per_row);
 
-  // now copy the entries from
-  // the other object
+  // now copy the entries from the other object
   for (size_type row=0; row<original.rows; ++row)
     {
-      // copy the elements of this row
-      // of the other object
+      // copy the elements of this row of the other object
       //
-      // note that the first object actually
-      // is the main-diagonal element,
+      // note that the first object actually is the main-diagonal element,
       // which we need not copy
       //
-      // we do the copying in two steps:
-      // first we note that the elements in
-      // @p{original} are sorted, so we may
-      // first copy all the elements up to
-      // the first side-diagonal one which
-      // is to be filled in. then we insert
-      // the side-diagonals, finally copy
-      // the rest from that element onwards
-      // which is not a side-diagonal any
-      // more.
+      // we do the copying in two steps: first we note that the elements in
+      // @p{original} are sorted, so we may first copy all the elements up to
+      // the first side-diagonal one which is to be filled in. then we insert
+      // the side-diagonals, finally copy the rest from that element onwards
+      // which is not a side-diagonal any more.
       const size_type *const
       original_row_start = &original.colnums[original.rowstart[row]] + 1;
-      // the following requires that
-      // @p{original} be compressed since
+      // the following requires that @p{original} be compressed since
       // otherwise there might be invalid_entry's
       const size_type *const
       original_row_end   = &original.colnums[original.rowstart[row+1]];
 
-      // find pointers before and
-      // after extra
-      // off-diagonals. if at top or
-      // bottom of matrix, then set
-      // these pointers such that no
-      // copying is necessary (see
-      // the @p{copy} commands)
+      // find pointers before and after extra off-diagonals. if at top or
+      // bottom of matrix, then set these pointers such that no copying is
+      // necessary (see the @p{copy} commands)
       const size_type *const
       original_last_before_side_diagonals
         = (row > extra_off_diagonals ?
@@ -242,9 +229,8 @@ SparsityPattern::SparsityPattern (const SparsityPattern &original,
                              row+extra_off_diagonals) :
            original_row_end);
 
-      // find first free slot. the
-      // first slot in each row is
-      // the diagonal element
+      // find first free slot. the first slot in each row is the diagonal
+      // element
       size_type *next_free_slot = &colnums[rowstart[row]] + 1;
 
       // copy elements before side-diagonals
@@ -265,12 +251,9 @@ SparsityPattern::SparsityPattern (const SparsityPattern &original,
                                   original_row_end,
                                   next_free_slot);
 
-      // this error may happen if the
-      // sum of previous elements per row
-      // and those of the new diagonals
-      // exceeds the maximum number of
-      // elements per row given to this
-      // constructor
+      // this error may happen if the sum of previous elements per row and
+      // those of the new diagonals exceeds the maximum number of elements per
+      // row given to this constructor
       Assert (next_free_slot <= &colnums[rowstart[row+1]],
               ExcNotEnoughSpace (0,rowstart[row+1]-rowstart[row]));
     };
@@ -310,8 +293,7 @@ SparsityPattern::reinit (const size_type m,
                          const unsigned int max_per_row,
                          const bool)
 {
-  // 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);
   reinit (m, n, row_lengths);
 }
@@ -323,8 +305,7 @@ SparsityPattern::reinit (const size_type m,
                          const size_type n,
                          const unsigned int max_per_row)
 {
-  // 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);
   reinit (m, n, row_lengths);
 }
@@ -360,31 +341,23 @@ SparsityPattern::reinit (const size_type m,
       rowstart = 0;
       colnums = 0;
       max_vec_len = max_dim = rows = cols = 0;
-      // if dimension is zero: ignore
-      // max_per_row
+      // if dimension is zero: ignore max_per_row
       max_row_length = 0;
       compressed = false;
       return;
     }
 
-  // first, if the matrix is
-  // quadratic, we
-  // will have to make sure that each
-  // row has at least one entry for
-  // the diagonal element. make this
-  // more obvious by having a
-  // variable which we can query
+  // first, if the matrix is quadratic, we will have to make sure that each
+  // row has at least one entry for the diagonal element. make this more
+  // obvious by having a variable which we can query
   store_diagonal_first_in_row = (m == n);
 
-  // find out how many entries we
-  // need in the @p{colnums} array. if
-  // this number is larger than
-  // @p{max_vec_len}, then we will need
-  // to reallocate memory
+  // find out how many entries we need in the @p{colnums} array. if this
+  // number is larger than @p{max_vec_len}, then we will need to reallocate
+  // memory
   //
-  // note that the number of elements
-  // per row is bounded by the number
-  // of columns
+  // note that the number of elements per row is bounded by the number of
+  // columns
   //
   std::size_t vec_len = 0;
   for (size_type i=0; i<m; ++i)
@@ -393,13 +366,9 @@ SparsityPattern::reinit (const size_type m,
                                                row_lengths[i]),
                         n);
 
-  // sometimes, no entries are
-  // requested in the matrix (this
-  // most often happens when blocks
-  // in a block matrix are simply
-  // zero). in that case, allocate
-  // exactly one element, to have a
-  // valid pointer to some memory
+  // sometimes, no entries are requested in the matrix (this most often
+  // happens when blocks in a block matrix are simply zero). in that case,
+  // allocate exactly one element, to have a valid pointer to some memory
   if (vec_len == 0)
     {
       vec_len = 1;
@@ -422,20 +391,12 @@ SparsityPattern::reinit (const size_type m,
   if (store_diagonal_first_in_row && (max_row_length==0) && (m!=0))
     max_row_length = 1;
 
-  // allocate memory for the rowstart
-  // values, if necessary. even
-  // though we re-set the pointers
-  // again immediately after deleting
-  // their old content, set them to
-  // zero in between because the
-  // allocation might fail, in which
-  // case we get an exception and the
-  // destructor of this object will
-  // be called -- where we look at
-  // the non-nullness of the (now
-  // invalid) pointer again and try
-  // to delete the memory a second
-  // time.
+  // allocate memory for the rowstart values, if necessary. even though we
+  // re-set the pointers again immediately after deleting their old content,
+  // set them to zero in between because the allocation might fail, in which
+  // case we get an exception and the destructor of this object will be called
+  // -- where we look at the non-nullness of the (now invalid) pointer again
+  // and try to delete the memory a second time.
   if (rows > max_dim)
     {
       if (rowstart)
@@ -448,8 +409,7 @@ SparsityPattern::reinit (const size_type m,
       rowstart = new std::size_t[max_dim+1];
     }
 
-  // allocate memory for the column
-  // numbers if necessary
+  // allocate memory for the column numbers if necessary
   if (vec_len > max_vec_len)
     {
       if (colnums)
@@ -475,14 +435,11 @@ SparsityPattern::reinit (const size_type m,
           ((vec_len == 1) && (rowstart[rows] == 0)),
           ExcInternalError());
 
-  // preset the column numbers by a
-  // value indicating it is not in
-  // use
+  // preset the column numbers by a value indicating it is not in use
   std::fill_n (&colnums[0], vec_len, invalid_entry);
 
-  // if diagonal elements are
-  // special: let the first entry in
-  // each row be the diagonal value
+  // if diagonal elements are special: let the first entry in each row be the
+  // diagonal value
   if (store_diagonal_first_in_row)
     for (size_type i=0; i<rows; i++)
       colnums[rowstart[i]] = i;
@@ -505,10 +462,8 @@ SparsityPattern::compress ()
             next_row_start  = 0,
             row_length      = 0;
 
-  // first find out how many non-zero
-  // elements there are, in order to
-  // allocate the right amount of
-  // memory
+  // first find out how many non-zero elements there are, in order to allocate
+  // the right amount of memory
   const std::size_t nonzero_elements
     = std::count_if (&colnums[rowstart[0]],
                      &colnums[rowstart[rows]],
@@ -517,69 +472,50 @@ SparsityPattern::compress ()
   size_type *new_colnums = new size_type[nonzero_elements];
 
 
-  // reserve temporary storage to
-  // store the entries of one row
+  // reserve temporary storage to store the entries of one row
   std::vector<size_type> tmp_entries (max_row_length);
 
   // Traverse all rows
   for (size_type line=0; line<rows; ++line)
     {
-      // copy used entries, break if
-      // first unused entry is reached
+      // copy used entries, break if first unused entry is reached
       row_length = 0;
       for (size_type j=rowstart[line]; j<rowstart[line+1]; ++j,++row_length)
         if (colnums[j] != invalid_entry)
           tmp_entries[row_length] = colnums[j];
         else
           break;
-      // now @p{rowstart} is
-      // the number of entries in
-      // this line
+      // now @p{rowstart} is the number of entries in this line
 
-      // Sort only beginning at the
-      // second entry, if optimized
-      // storage of diagonal entries
-      // is on.
+      // Sort only beginning at the second entry, if optimized storage of
+      // diagonal entries is on.
 
-      // if this line is empty or has
-      // only one entry, don't sort
+      // if this line is empty or has only one entry, don't sort
       if (row_length > 1)
         std::sort ((store_diagonal_first_in_row)
                    ? tmp_entries.begin()+1
                    : tmp_entries.begin(),
                    tmp_entries.begin()+row_length);
 
-      // insert column numbers
-      // into the new field
+      // insert column numbers into the new field
       for (size_type j=0; j<row_length; ++j)
         new_colnums[next_free_entry++] = tmp_entries[j];
 
-      // note new start of this and
-      // the next row
+      // note new start of this and the next row
       rowstart[line] = next_row_start;
       next_row_start = next_free_entry;
 
-      // some internal checks: either
-      // the matrix is not quadratic,
-      // or if it is, then the first
-      // element of this row must be
-      // the diagonal element
-      // (i.e. with column
-      // index==line number)
+      // some internal checks: either the matrix is not quadratic, or if it
+      // is, then the first element of this row must be the diagonal element
+      // (i.e. with column index==line number)
       Assert ((!store_diagonal_first_in_row) ||
               (new_colnums[rowstart[line]] == line),
               ExcInternalError());
-      // assert that the first entry
-      // does not show up in
-      // the remaining ones and that
-      // the remaining ones are unique
-      // among themselves (this handles
-      // both cases, quadratic and
-      // rectangular matrices)
+      // assert that the first entry does not show up in the remaining ones
+      // and that the remaining ones are unique among themselves (this handles
+      // both cases, quadratic and rectangular matrices)
       //
-      // the only exception here is
-      // if the row contains no
-      // entries at all
+      // the only exception here is if the row contains no entries at all
       Assert ((rowstart[line] == next_row_start)
               ||
               (std::find (&new_colnums[rowstart[line]+1],
@@ -595,18 +531,14 @@ SparsityPattern::compress ()
               ExcInternalError());
     };
 
-  // assert that we have used all
-  // allocated space, no more and no
-  // less
+  // assert that we have used all allocated space, no more and no less
   Assert (next_free_entry == nonzero_elements,
           ExcInternalError());
 
   // set iterator-past-the-end
   rowstart[rows] = next_row_start;
 
-  // set colnums to the newly
-  // allocated array and delete the
-  // old one
+  // set colnums to the newly allocated array and delete the old one
   delete[] colnums;
   colnums = new_colnums;
 
@@ -632,16 +564,10 @@ template <typename CSP>
 void
 SparsityPattern::copy_from (const CSP &csp)
 {
-  // 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 do_diag_optimize = (csp.n_rows() == csp.n_cols());
   std::vector<unsigned int> row_lengths (csp.n_rows());
   for (size_type i=0; i<csp.n_rows(); ++i)
@@ -652,11 +578,8 @@ SparsityPattern::copy_from (const CSP &csp)
     }
   reinit (csp.n_rows(), csp.n_cols(), row_lengths);
 
-  // now enter all the elements into
-  // the matrix, if there are
-  // any. note that if the matrix is
-  // quadratic, then we already have
-  // the diagonal element
+  // now enter all the elements into the matrix, if there are any. note that
+  // if the matrix is quadratic, then we already have the diagonal element
   // preallocated
   if (n_rows() != 0 && n_cols() != 0)
     for (size_type row = 0; row<csp.n_rows(); ++row)
@@ -673,10 +596,8 @@ SparsityPattern::copy_from (const CSP &csp)
           }
       }
 
-  // do not need to compress the sparsity
-  // pattern since we already have
-  // allocated the right amount of data,
-  // and the CSP data is sorted, too.
+  // do not need to compress the sparsity pattern since we already have
+  // allocated the right amount of data, and the CSP data is sorted, too.
   compressed = true;
 }
 
@@ -694,11 +615,9 @@ void SparsityPattern::copy_from (const FullMatrix<number> &matrix,
 template <typename number>
 void SparsityPattern::copy_from (const FullMatrix<number> &matrix)
 {
-  // first init with the number of entries
-  // per row. if this matrix is square then we
-  // also have to allocate memory for the
-  // diagonal entry, unless we have already
-  // counted it
+  // first init with the number of entries per row. if this matrix is square
+  // then we also have to allocate memory for the diagonal entry, unless we
+  // have already counted it
   std::vector<unsigned int> entries_per_row (matrix.m(), 0);
   for (size_type row=0; row<matrix.m(); ++row)
     {
@@ -751,15 +670,11 @@ SparsityPattern::reinit (const size_type               m,
 bool
 SparsityPattern::empty () const
 {
-  // let's try to be on the safe side of
-  // life by using multiple possibilities in
-  // the check for emptiness... (sorry for
-  // this kludge -- emptying matrices and
-  // freeing memory was not present in the
-  // original implementation and I don't
-  // know at how many places I missed
-  // something in adding it, so I try to
-  // be cautious. wb)
+  // let's try to be on the safe side of life by using multiple possibilities
+  // in the check for emptiness... (sorry for this kludge -- emptying matrices
+  // and freeing memory was not present in the original implementation and I
+  // don't know at how many places I missed something in adding it, so I try
+  // to be cautious. wb)
   if ((rowstart==0) || (rows==0) || (cols==0))
     {
       Assert (rowstart==0, ExcInternalError());
@@ -778,18 +693,15 @@ SparsityPattern::empty () const
 SparsityPattern::size_type
 SparsityPattern::max_entries_per_row () const
 {
-  // if compress() has not yet been
-  // called, we can get the maximum
-  // number of elements per row using
-  // the stored value
+  // if compress() has not yet been called, we can get the maximum number of
+  // elements per row using the stored value
   if (!compressed)
     return max_row_length;
 
-  // if compress() was called, we
-  // use a better algorithm which
-  // gives us a sharp bound
+  // if compress() was called, we use a better algorithm which gives us a
+  // sharp bound
   size_type m = 0;
-  for (size_type i=1; i<rows; ++i)
+  for (size_type i=1; i<=rows; ++i)
     m = std::max (m, static_cast<size_type>(rowstart[i]-rowstart[i-1]));
 
   return m;
@@ -806,30 +718,21 @@ SparsityPattern::operator () (const size_type i,
   Assert (j<cols, ExcIndexRange(j,0,cols));
   Assert (compressed, ExcNotCompressed());
 
-  // let's see whether there is
-  // something in this line
+  // let's see whether there is something in this line
   if (rowstart[i] == rowstart[i+1])
     return invalid_entry;
 
-  // If special storage of diagonals
-  // was requested, we can get the
-  // diagonal element faster by this
-  // query.
+  // If special storage of diagonals was requested, we can get the diagonal
+  // element faster by this query.
   if (store_diagonal_first_in_row && (i==j))
     return rowstart[i];
 
-  // all other entries are sorted, so
-  // we can use a binary search
-  // algorithm
+  // all other entries are sorted, so we can use a binary search algorithm
   //
-  // note that the entries are only
-  // sorted upon compression, so this
-  // would fail for non-compressed
-  // sparsity patterns; however, that
-  // is why the Assertion is at the
-  // top of this function, so it may
-  // not be called for noncompressed
-  // structures.
+  // note that the entries are only sorted upon compression, so this would
+  // fail for non-compressed sparsity patterns; however, that is why the
+  // Assertion is at the top of this function, so it may not be called for
+  // noncompressed structures.
   const size_type *sorted_region_start = (store_diagonal_first_in_row ?
                                           &colnums[rowstart[i]+1] :
                                           &colnums[rowstart[i]]);
@@ -959,22 +862,16 @@ SparsityPattern::matrix_position (const size_type global_index) const
   Assert (global_index < n_nonzero_elements(),
           ExcIndexRange (global_index, 0, n_nonzero_elements()));
 
-  // first find the row in which the
-  // entry is located. for this note
-  // that the rowstart array indexes
-  // the global indices at which each
-  // row starts. since it is sorted,
-  // and since there is an element
-  // for the one-past-last row, we
-  // can simply use a bisection
-  // search on it
+  // first find the row in which the entry is located. for this note that the
+  // rowstart array indexes the global indices at which each row starts. since
+  // it is sorted, and since there is an element for the one-past-last row, we
+  // can simply use a bisection search on it
   const size_type row
     = (std::upper_bound (&rowstart[0], &rowstart[rows], global_index)
        - &rowstart[0] - 1);
 
-  // now, the column index is simple
-  // since that is what the colnums
-  // array stores:
+  // now, the column index is simple since that is what the colnums array
+  // stores:
   const size_type col = colnums[global_index];
 
   // so return the respective pair
@@ -988,37 +885,28 @@ SparsityPattern::symmetrize ()
 {
   Assert ((rowstart!=0) && (colnums!=0), ExcEmptyObject());
   Assert (compressed==false, ExcMatrixIsCompressed());
-  // Note that we only require a
-  // quadratic matrix here, no
-  // special treatment of diagonals
+  // Note that we only require a quadratic matrix here, no special treatment
+  // of diagonals
   Assert (rows==cols, ExcNotQuadratic());
 
-  // loop over all elements presently
-  // in the sparsity pattern and add
-  // the transpose element. note:
+  // loop over all elements presently in the sparsity pattern and add the
+  // transpose element. note:
   //
-  // 1. that the sparsity pattern
-  // changes which we work on, but
-  // not the present row
+  // 1. that the sparsity pattern changes which we work on, but not the
+  // present row
   //
-  // 2. that the @p{add} function can
-  // be called on elements that
-  // already exist without any harm
+  // 2. that the @p{add} function can be called on elements that already exist
+  // without any harm
   for (size_type row=0; row<rows; ++row)
     for (size_type k=rowstart[row]; k<rowstart[row+1]; k++)
       {
-        // check whether we are at
-        // the end of the entries of
-        // this row. if so, go to
-        // next row
+        // check whether we are at the end of the entries of this row. if so,
+        // go to next row
         if (colnums[k] == invalid_entry)
           break;
 
-        // otherwise add the
-        // transpose entry if this is
-        // not the diagonal (that
-        // would not harm, only take
-        // time to check up)
+        // otherwise add the transpose entry if this is not the diagonal (that
+        // would not harm, only take time to check up)
         if (colnums[k] != row)
           add (colnums[k], row);
       };
@@ -1057,13 +945,9 @@ SparsityPattern::print_gnuplot (std::ostream &out) const
   for (size_type i=0; i<rows; ++i)
     for (size_type j=rowstart[i]; j<rowstart[i+1]; ++j)
       if (colnums[j] != invalid_entry)
-        // 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 << colnums[j] << " " << -static_cast<signed int>(i) << std::endl;
 
   AssertThrow (out, ExcIO());
@@ -1084,8 +968,7 @@ SparsityPattern::bandwidth () const
             b = std::abs(static_cast<signed int>(i-colnums[j]));
         }
       else
-        // leave if at the end of
-        // the entries of this line
+        // leave if at the end of the entries of this line
         break;
   return b;
 }
@@ -1096,8 +979,7 @@ SparsityPattern::block_write (std::ostream &out) const
 {
   AssertThrow (out, ExcIO());
 
-  // first the simple objects,
-  // bracketed in [...]
+  // first the simple objects, bracketed in [...]
   out << '[' << max_dim << ' '
       << rows << ' '
       << cols << ' '
diff --git a/tests/bits/chunk_sparse_matrix_16.cc b/tests/bits/chunk_sparse_matrix_16.cc
new file mode 100644
index 0000000000..f9a99c45f7
--- /dev/null
+++ b/tests/bits/chunk_sparse_matrix_16.cc
@@ -0,0 +1,115 @@
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 2004 - 2013 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+
+// compare matrix-vector product in chunk sparse matrix with usual vmult and
+// extract_row_copy.
+
+#include "../tests.h"
+#include <deal.II/lac/chunk_sparse_matrix.h>
+#include <deal.II/lac/vector.h>
+#include <fstream>
+#include <iomanip>
+
+
+void test (const unsigned int chunk_size)
+{
+  deallog << "Chunk size = " << chunk_size << std::endl;
+
+  for (unsigned int n_cols = 4; n_cols<7; ++n_cols)
+    {
+      deallog << "n_cols = " << n_cols << std::endl;
+      ChunkSparsityPattern sp (5,n_cols,3,chunk_size);
+      for (unsigned int i=0; i<5; ++i)
+        for (unsigned int j=0; j<n_cols; ++j)
+          if ((i+2*j+1) % 3 == 0)
+            sp.add (i,j);
+      sp.compress ();
+
+      ChunkSparseMatrix<double> m(sp);
+
+      // first set a few entries
+      for (unsigned int i=0; i<m.m(); ++i)
+        for (unsigned int j=0; j<m.n(); ++j)
+          if ((i+2*j+1) % 3 == 0)
+            m.set (i,j, i*j*.5+.5);
+
+      // next perform a matrix-vector product using the entries as given by
+      // the entries in extract_row_copy and compare it with the exact value
+      Vector<double> src(m.n()), dst(m.m()), dst_ref(m.m());
+      for (unsigned int i=0; i<src.size(); ++i)
+        src(i) = (double)Testing::rand()/RAND_MAX;
+      std::vector<unsigned int> indices(sp.max_entries_per_row());
+      std::vector<double> values(sp.max_entries_per_row());
+      for (unsigned int i=0; i<m.m(); ++i)
+        {
+          unsigned int n_entries = numbers::invalid_unsigned_int;
+          m.extract_row_copy(i, values.size(), n_entries, &indices[0], &values[0]);
+          double sum = 0;
+          for (unsigned int j=0; j<n_entries; ++j)
+            sum += values[j] * src(indices[j]);
+          dst(i) = sum;
+        }
+      m.vmult(dst_ref, src);
+      dst -= dst_ref;
+      deallog << "Error in matrix-vector product done via extract_row_copy: "
+              << dst.linfty_norm() << std::endl;
+    }
+}
+
+
+
+int main ()
+{
+  std::ofstream logfile("output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  try
+    {
+      const unsigned int chunk_sizes[] = { 1, 2, 4, 5, 7 };
+      for (unsigned int i=0;
+           i<sizeof(chunk_sizes)/sizeof(chunk_sizes[0]);
+           ++i)
+        test (chunk_sizes[i]);
+    }
+  catch (std::exception &exc)
+    {
+      deallog << std::endl << std::endl
+              << "----------------------------------------------------"
+              << std::endl;
+      deallog << "Exception on processing: " << std::endl
+              << exc.what() << std::endl
+              << "Aborting!" << std::endl
+              << "----------------------------------------------------"
+              << std::endl;
+
+      return 1;
+    }
+  catch (...)
+    {
+      deallog << std::endl << std::endl
+              << "----------------------------------------------------"
+              << std::endl;
+      deallog << "Unknown exception!" << std::endl
+              << "Aborting!" << std::endl
+              << "----------------------------------------------------"
+              << std::endl;
+      return 1;
+    };
+}
diff --git a/tests/bits/chunk_sparse_matrix_16.output b/tests/bits/chunk_sparse_matrix_16.output
new file mode 100644
index 0000000000..8b12087e16
--- /dev/null
+++ b/tests/bits/chunk_sparse_matrix_16.output
@@ -0,0 +1,36 @@
+
+DEAL::Chunk size = 1
+DEAL::n_cols = 4
+DEAL::Error in matrix-vector product done via extract_row_copy: 0
+DEAL::n_cols = 5
+DEAL::Error in matrix-vector product done via extract_row_copy: 0
+DEAL::n_cols = 6
+DEAL::Error in matrix-vector product done via extract_row_copy: 0
+DEAL::Chunk size = 2
+DEAL::n_cols = 4
+DEAL::Error in matrix-vector product done via extract_row_copy: 0
+DEAL::n_cols = 5
+DEAL::Error in matrix-vector product done via extract_row_copy: 0
+DEAL::n_cols = 6
+DEAL::Error in matrix-vector product done via extract_row_copy: 0
+DEAL::Chunk size = 4
+DEAL::n_cols = 4
+DEAL::Error in matrix-vector product done via extract_row_copy: 0
+DEAL::n_cols = 5
+DEAL::Error in matrix-vector product done via extract_row_copy: 0
+DEAL::n_cols = 6
+DEAL::Error in matrix-vector product done via extract_row_copy: 0
+DEAL::Chunk size = 5
+DEAL::n_cols = 4
+DEAL::Error in matrix-vector product done via extract_row_copy: 0
+DEAL::n_cols = 5
+DEAL::Error in matrix-vector product done via extract_row_copy: 0
+DEAL::n_cols = 6
+DEAL::Error in matrix-vector product done via extract_row_copy: 0
+DEAL::Chunk size = 7
+DEAL::n_cols = 4
+DEAL::Error in matrix-vector product done via extract_row_copy: 0
+DEAL::n_cols = 5
+DEAL::Error in matrix-vector product done via extract_row_copy: 0
+DEAL::n_cols = 6
+DEAL::Error in matrix-vector product done via extract_row_copy: 0