]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Implement iterator for ChunkSparseMatrix. Parallelize matrix-vector product for Chunk...
authorkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 7 May 2013 17:14:12 +0000 (17:14 +0000)
committerkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 7 May 2013 17:14:12 +0000 (17:14 +0000)
git-svn-id: https://svn.dealii.org/trunk@29476 0785d39b-7218-0410-832d-ea1e28bc413d

40 files changed:
deal.II/include/deal.II/lac/chunk_sparse_matrix.h
deal.II/include/deal.II/lac/chunk_sparse_matrix.templates.h
deal.II/include/deal.II/lac/chunk_sparsity_pattern.h
deal.II/source/lac/chunk_sparsity_pattern.cc
tests/bits/chunk_sparse_matrix_11.cc [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_11/cmp/generic [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_12.cc [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_12/cmp/generic [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_13.cc [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_13/cmp/generic [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_14.cc [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_14/cmp/generic [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_15.cc [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_15/cmp/generic [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_iterator_01.cc [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_iterator_01/cmp/generic [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_iterator_02.cc [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_iterator_02/cmp/generic [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_iterator_03.cc [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_iterator_03/cmp/generic [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_iterator_04.cc [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_iterator_04/cmp/generic [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_iterator_05.cc [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_iterator_05/cmp/generic [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_iterator_06.cc [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_iterator_06/cmp/generic [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_iterator_07.cc [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_iterator_07/cmp/generic [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_iterator_08.cc [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_iterator_08/cmp/generic [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_iterator_09.cc [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_iterator_09/cmp/generic [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_iterator_10.cc [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_iterator_10/cmp/generic [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_iterator_11.cc [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_iterator_11/cmp/generic [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_iterator_12.cc [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_iterator_12/cmp/generic [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_iterator_13.cc [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_iterator_13/cmp/generic [new file with mode: 0644]

index 9c20d69face932dc5943dd287eb5f9fe692845c2..3d7a5c2fc3518606c90d4676b3d1df72ded2e8ea 100644 (file)
@@ -29,6 +29,359 @@ template<typename number> class FullMatrix;
  *@{
  */
 
+/**
+ * A namespace in which we declare iterators over the elements of sparse
+ * matrices.
+ */
+namespace ChunkSparseMatrixIterators
+{
+  // forward declaration
+  template <typename number, bool Constness>
+  class Iterator;
+
+  /**
+   * General template for sparse matrix accessors. The first template argument
+   * denotes the underlying numeric type, the second the constness of the
+   * matrix.
+   *
+   * The general template is not implemented, only the specializations for the
+   * two possible values of the second template argument. Therefore, the
+   * interface listed here only serves as a template provided since doxygen
+   * does not link the specializations.
+   */
+  template <typename number, bool Constness>
+  class Accessor : public ChunkSparsityPatternIterators::Accessor
+  {
+  public:
+    /**
+     * Value of this matrix entry.
+     */
+    number value() const;
+
+    /**
+     * Value of this matrix entry.
+     */
+    number &value();
+
+    /**
+     * Return a reference to the matrix into which this accessor points. Note
+     * that in the present case, this is a constant reference.
+     */
+    const ChunkSparseMatrix<number> &get_matrix () const;
+  };
+
+
+
+  /**
+   * Accessor class for constant matrices, used in the const_iterators. This
+   * class builds on the accessor classes used for sparsity patterns to loop
+   * over all nonzero entries, and only adds the accessor functions to gain
+   * access to the actual value stored at a certain location.
+   */
+  template <typename number>
+  class Accessor<number,true> : public ChunkSparsityPatternIterators::Accessor
+  {
+  public:
+    /**
+     * Typedef for the type (including constness) of the matrix to be used
+     * here.
+     */
+    typedef const ChunkSparseMatrix<number> MatrixType;
+
+    /**
+     * Constructor.
+     */
+    Accessor (MatrixType         *matrix,
+              const unsigned int  row);
+
+    /**
+     * Constructor. Construct the end accessor for the given matrix.
+     */
+    Accessor (MatrixType         *matrix);
+
+    /**
+     * Copy constructor to get from a non-const accessor to a const accessor.
+     */
+    Accessor (const ChunkSparseMatrixIterators::Accessor<number,false> &a);
+
+    /**
+     * Value of this matrix entry.
+     */
+    number value() const;
+
+    /**
+     * Return a reference to the matrix into which this accessor points. Note
+     * that in the present case, this is a constant reference.
+     */
+    MatrixType &get_matrix () const;
+
+  private:
+    /**
+     * Pointer to the matrix we use.
+     */
+    MatrixType *matrix;
+
+    /**
+     * Make the advance function of the base class available.
+     */
+    using ChunkSparsityPatternIterators::Accessor::advance;
+
+    /**
+     * Make iterator class a friend.
+     */
+    template <typename, bool>
+    friend class Iterator;
+  };
+
+
+  /**
+   * Accessor class for non-constant matrices, used in the iterators. This
+   * class builds on the accessor classes used for sparsity patterns to loop
+   * over all nonzero entries, and only adds the accessor functions to gain
+   * access to the actual value stored at a certain location.
+   */
+  template <typename number>
+  class Accessor<number,false> : public ChunkSparsityPatternIterators::Accessor
+  {
+  private:
+    /**
+     * Reference class. This is what the accessor class returns when you call
+     * the value() function. The reference acts just as if it were a reference
+     * to the actual value of a matrix entry, i.e. you can read and write it,
+     * you can add and multiply to it, etc, but since the matrix does not give
+     * away the address of this matrix entry, we have to go through functions
+     * to do all this.
+     *
+     * The constructor takes a pointer to an accessor object that describes
+     * which element of the matrix it points to. This creates an ambiguity
+     * when one writes code like iterator->value()=0 (instead of
+     * iterator->value()=0.0), since the right hand side is an integer that
+     * can both be converted to a <tt>number</tt> (i.e., most commonly a
+     * double) or to another object of type <tt>Reference</tt>. The compiler
+     * then complains about not knowing which conversion to take.
+     *
+     * For some reason, adding another overload operator=(int) doesn't seem to
+     * cure the problem. We avoid it, however, by adding a second, dummy
+     * argument to the Reference constructor, that is unused, but makes sure
+     * there is no second matching conversion sequence using a one-argument
+     * right hand side.
+     */
+    class Reference
+    {
+    public:
+      /**
+       * Constructor. For the second argument, see the general class
+       * documentation.
+       */
+      Reference (const Accessor *accessor,
+                 const bool dummy);
+
+      /**
+       * Conversion operator to the data type of the matrix.
+       */
+      operator number () const;
+
+      /**
+       * Set the element of the matrix we presently point to to @p n.
+       */
+      const Reference &operator = (const number n) const;
+
+      /**
+       * Add @p n to the element of the matrix we presently point to.
+       */
+      const Reference &operator += (const number n) const;
+
+      /**
+       * Subtract @p n from the element of the matrix we presently point to.
+       */
+      const Reference &operator -= (const number n) const;
+
+      /**
+       * Multiply the element of the matrix we presently point to by @p n.
+       */
+      const Reference &operator *= (const number n) const;
+
+      /**
+       * Divide the element of the matrix we presently point to by @p n.
+       */
+      const Reference &operator /= (const number n) const;
+
+    private:
+      /**
+       * Pointer to the accessor that denotes which element we presently point
+       * to.
+       */
+      const Accessor *accessor;
+    };
+
+  public:
+    /**
+     * Typedef for the type (including constness) of the matrix to be used
+     * here.
+     */
+    typedef ChunkSparseMatrix<number> MatrixType;
+
+    /**
+     * Constructor.
+     */
+    Accessor (MatrixType         *matrix,
+              const unsigned int  row);
+
+    /**
+     * Constructor. Construct the end accessor for the given matrix.
+     */
+    Accessor (MatrixType         *matrix);
+
+    /**
+     * Value of this matrix entry, returned as a read- and writable reference.
+     */
+    Reference value() const;
+
+    /**
+     * Return a reference to the matrix into which this accessor points. Note
+     * that in the present case, this is a non-constant reference.
+     */
+    MatrixType &get_matrix () const;
+
+  private:
+    /**
+     * Pointer to the matrix we use.
+     */
+    MatrixType *matrix;
+
+    /**
+     * Make the advance function of the base class available.
+     */
+    using ChunkSparsityPatternIterators::Accessor::advance;
+
+    /**
+     * Make iterator class a friend.
+     */
+    template <typename, bool>
+    friend class Iterator;
+
+    /**
+     * Make the inner reference class a friend if the compiler has a bug and
+     * requires this.
+     */
+  };
+
+
+
+  /**
+   * STL conforming iterator for constant and non-constant matrices.
+   *
+   * The first template argument denotes the underlying numeric type, the
+   * second the constness of the matrix.
+   *
+   * Since there is a specialization of this class for
+   * <tt>Constness=false</tt>, this class is for iterators to constant
+   * matrices.
+   */
+  template <typename number, bool Constness>
+  class Iterator
+  {
+  public:
+    /**
+     * Typedef for the matrix type (including constness) we are to operate on.
+     */
+    typedef
+    typename Accessor<number,Constness>::MatrixType
+    MatrixType;
+
+    /**
+     * A typedef for the type you get when you dereference an iterator
+     * of the current kind.
+     */
+    typedef
+    const Accessor<number,Constness> & value_type;
+
+    /**
+     * Constructor. Create an iterator into the matrix @p matrix for the given
+     * row and the index within it.
+     */
+    Iterator (MatrixType        *matrix,
+              const unsigned int row);
+
+    /**
+     * Constructor. Create the end iterator for the given matrix.
+     */
+    Iterator (MatrixType *matrix);
+
+    /**
+     * Conversion constructor to get from a non-const iterator to a const
+     * iterator.
+     */
+    Iterator (const ChunkSparseMatrixIterators::Iterator<number,false> &i);
+
+    /**
+     * Prefix increment.
+     */
+    Iterator &operator++ ();
+
+    /**
+     * Postfix increment.
+     */
+    Iterator operator++ (int);
+
+    /**
+     * Dereferencing operator.
+     */
+    const Accessor<number,Constness> &operator* () const;
+
+    /**
+     * Dereferencing operator.
+     */
+    const Accessor<number,Constness> *operator-> () const;
+
+    /**
+     * Comparison. True, if both iterators point to the same matrix position.
+     */
+    bool operator == (const Iterator &) const;
+
+    /**
+     * Inverse of <tt>==</tt>.
+     */
+    bool operator != (const Iterator &) const;
+
+    /**
+     * Comparison operator. Result is true if either the first row number is
+     * smaller or if the row numbers are equal and the first index is smaller.
+     *
+     * This function is only valid if both iterators point into the same
+     * matrix.
+     */
+    bool operator < (const Iterator &) const;
+
+    /**
+     * Comparison operator. Works in the same way as above operator, just the
+     * other way round.
+     */
+    bool operator > (const Iterator &) const;
+
+    /**
+     * Return the distance between the current iterator and the argument.
+     * The distance is given by how many times one has to apply operator++
+     * to the current iterator to get the argument (for a positive return
+     * value), or operator-- (for a negative return value).
+     */
+    int operator - (const Iterator &p) const;
+
+    /**
+     * Return an iterator that is @p n ahead of the current one.
+     */
+    Iterator operator + (const unsigned int n) const;
+
+  private:
+    /**
+     * Store an object of the accessor class.
+     */
+    Accessor<number,Constness> accessor;
+  };
+
+}
+
+
 
 /**
  * Sparse matrix. This class implements the function to store values
@@ -56,42 +409,46 @@ public:
   typedef number value_type;
 
   /**
-   * Declare a type that has holds
-   * real-valued numbers with the
-   * same precision as the template
-   * argument to this class. If the
-   * template argument of this
-   * class is a real data type,
-   * then real_type equals the
-   * template argument. If the
-   * template argument is a
-   * std::complex type then
-   * real_type equals the type
-   * underlying the complex
-   * numbers.
+   * Declare a type that has holds real-valued numbers with the same precision
+   * as the template argument to this class. If the template argument of this
+   * class is a real data type, then real_type equals the template
+   * argument. If the template argument is a std::complex type then real_type
+   * equals the type underlying the complex numbers.
    *
-   * This typedef is used to
-   * represent the return type of
-   * norms.
+   * This typedef is used to represent the return type of norms.
    */
   typedef typename numbers::NumberTraits<number>::real_type real_type;
 
   /**
-   * A structure that describes some of the
-   * traits of this class in terms of its
-   * run-time behavior. Some other classes
-   * (such as the block matrix classes)
-   * that take one or other of the matrix
-   * classes as its template parameters can
-   * tune their behavior based on the
-   * variables in this class.
+   * Typedef of an STL conforming iterator class walking over all the nonzero
+   * entries of this matrix. This iterator cannot change the values of the
+   * matrix.
+   */
+  typedef
+  ChunkSparseMatrixIterators::Iterator<number,true>
+  const_iterator;
+
+  /**
+   * Typedef of an STL conforming iterator class walking over all the nonzero
+   * entries of this matrix. This iterator @em can change the values of the
+   * matrix, but of course can't change the sparsity pattern as this is fixed
+   * once a sparse matrix is attached to it.
+   */
+  typedef
+  ChunkSparseMatrixIterators::Iterator<number,false>
+  iterator;
+
+  /**
+   * A structure that describes some of the traits of this class in terms of
+   * its run-time behavior. Some other classes (such as the block matrix
+   * classes) that take one or other of the matrix classes as its template
+   * parameters can tune their behavior based on the variables in this class.
    */
   struct Traits
   {
     /**
-     * It is safe to elide additions of
-     * zeros to individual elements of
-     * this matrix.
+     * It is safe to elide additions of zeros to individual elements of this
+     * matrix.
      */
     static const bool zero_addition_can_be_elided = true;
   };
@@ -101,172 +458,106 @@ public:
    */
 //@{
   /**
-   * Constructor; initializes the matrix to
-   * be empty, without any structure, i.e.
-   * the matrix is not usable at all. This
-   * constructor is therefore only useful
-   * for matrices which are members of a
-   * class. All other matrices should be
-   * created at a point in the data flow
-   * where all necessary information is
-   * available.
+   * Constructor; initializes the matrix to be empty, without any structure,
+   * i.e.  the matrix is not usable at all. This constructor is therefore only
+   * useful for matrices which are members of a class. All other matrices
+   * should be created at a point in the data flow where all necessary
+   * information is available.
    *
-   * You have to initialize
-   * the matrix before usage with
-   * reinit(const ChunkSparsityPattern&).
+   * You have to initialize the matrix before usage with reinit(const
+   * ChunkSparsityPattern&).
    */
   ChunkSparseMatrix ();
 
   /**
-   * Copy constructor. This constructor is
-   * only allowed to be called if the matrix
-   * to be copied is empty. This is for the
-   * same reason as for the
-   * ChunkSparsityPattern, see there for the
-   * details.
+   * Copy constructor. This constructor is only allowed to be called if the
+   * matrix to be copied is empty. This is for the same reason as for the
+   * ChunkSparsityPattern, see there for the details.
    *
-   * If you really want to copy a whole
-   * matrix, you can do so by using the
+   * If you really want to copy a whole matrix, you can do so by using the
    * copy_from() function.
    */
   ChunkSparseMatrix (const ChunkSparseMatrix &);
 
   /**
-   * Constructor. Takes the given
-   * matrix sparsity structure to
-   * represent the sparsity pattern
-   * of this matrix. You can change
-   * the sparsity pattern later on
-   * by calling the reinit(const
-   * ChunkSparsityPattern&) function.
+   * Constructor. Takes the given matrix sparsity structure to represent the
+   * sparsity pattern of this matrix. You can change the sparsity pattern
+   * later on by calling the reinit(const ChunkSparsityPattern&) function.
    *
-   * You have to make sure that the
-   * lifetime of the sparsity
-   * structure is at least as long
-   * as that of this matrix or as
-   * long as reinit(const
-   * ChunkSparsityPattern&) is not
-   * called with a new sparsity
-   * pattern.
+   * You have to make sure that the lifetime of the sparsity structure is at
+   * least as long as that of this matrix or as long as reinit(const
+   * ChunkSparsityPattern&) is not called with a new sparsity pattern.
    *
-   * The constructor is marked
-   * explicit so as to disallow
-   * that someone passes a sparsity
-   * pattern in place of a sparse
-   * matrix to some function, where
-   * an empty matrix would be
-   * generated then.
+   * The constructor is marked explicit so as to disallow that someone passes
+   * a sparsity pattern in place of a sparse matrix to some function, where an
+   * empty matrix would be generated then.
    */
   explicit ChunkSparseMatrix (const ChunkSparsityPattern &sparsity);
 
   /**
-   * Copy constructor: initialize
-   * the matrix with the identity
-   * matrix. This constructor will
-   * throw an exception if the
-   * sizes of the sparsity pattern
-   * and the identity matrix do not
-   * coincide, or if the sparsity
-   * pattern does not provide for
-   * nonzero entries on the entire
-   * diagonal.
+   * Copy constructor: initialize the matrix with the identity matrix. This
+   * constructor will throw an exception if the sizes of the sparsity pattern
+   * and the identity matrix do not coincide, or if the sparsity pattern does
+   * not provide for nonzero entries on the entire diagonal.
    */
   ChunkSparseMatrix (const ChunkSparsityPattern &sparsity,
                      const IdentityMatrix  &id);
 
   /**
-   * Destructor. Free all memory, but do not
-   * release the memory of the sparsity
-   * structure.
+   * Destructor. Free all memory, but do not release the memory of the
+   * sparsity structure.
    */
   virtual ~ChunkSparseMatrix ();
 
   /**
-   * Copy operator. Since copying
-   * entire sparse matrices is a
-   * very expensive operation, we
-   * disallow doing so except for
-   * the special case of empty
-   * matrices of size zero. This
-   * doesn't seem particularly
-   * useful, but is exactly what
-   * one needs if one wanted to
-   * have a
-   * <code>std::vector@<ChunkSparseMatrix@<double@>
-   * @></code>: in that case, one
-   * can create a vector (which
-   * needs the ability to copy
-   * objects) of empty matrices
-   * that are then later filled
-   * with something useful.
+   * Copy operator. Since copying entire sparse matrices is a very expensive
+   * operation, we disallow doing so except for the special case of empty
+   * matrices of size zero. This doesn't seem particularly useful, but is
+   * exactly what one needs if one wanted to have a
+   * <code>std::vector@<ChunkSparseMatrix@<double@> @></code>: in that case,
+   * one can create a vector (which needs the ability to copy objects) of
+   * empty matrices that are then later filled with something useful.
    */
   ChunkSparseMatrix<number> &operator = (const ChunkSparseMatrix<number> &);
 
   /**
-   * Copy operator: initialize
-   * the matrix with the identity
-   * matrix. This operator will
-   * throw an exception if the
-   * sizes of the sparsity pattern
-   * and the identity matrix do not
-   * coincide, or if the sparsity
-   * pattern does not provide for
-   * nonzero entries on the entire
-   * diagonal.
+   * Copy operator: initialize the matrix with the identity matrix. This
+   * operator will throw an exception if the sizes of the sparsity pattern and
+   * the identity matrix do not coincide, or if the sparsity pattern does not
+   * provide for nonzero entries on the entire diagonal.
    */
   ChunkSparseMatrix<number> &
   operator= (const IdentityMatrix  &id);
 
   /**
-   * This operator assigns a scalar to
-   * a matrix. Since this does usually
-   * not make much sense (should we set
-   * all matrix entries to this value?
-   * Only the nonzero entries of the
-   * sparsity pattern?), this operation
-   * is only allowed if the actual
-   * value to be assigned is zero. This
-   * operator only exists to allow for
-   * the obvious notation
-   * <tt>matrix=0</tt>, which sets all
-   * elements of the matrix to zero,
-   * but keep the sparsity pattern
+   * This operator assigns a scalar to a matrix. Since this does usually not
+   * make much sense (should we set all matrix entries to this value?  Only
+   * the nonzero entries of the sparsity pattern?), this operation is only
+   * allowed if the actual value to be assigned is zero. This operator only
+   * exists to allow for the obvious notation <tt>matrix=0</tt>, which sets
+   * all elements of the matrix to zero, but keep the sparsity pattern
    * previously used.
    */
   ChunkSparseMatrix &operator = (const double d);
 
   /**
-   * Reinitialize the sparse matrix
-   * with the given sparsity
-   * pattern. The latter tells the
-   * matrix how many nonzero
-   * elements there need to be
+   * Reinitialize the sparse matrix with the given sparsity pattern. The
+   * latter tells the matrix how many nonzero elements there need to be
    * reserved.
    *
-   * Regarding memory allocation,
-   * the same applies as said
-   * above.
+   * Regarding memory allocation, the same applies as said above.
    *
-   * You have to make sure that the
-   * lifetime of the sparsity
-   * structure is at least as long
-   * as that of this matrix or as
-   * long as reinit(const
-   * ChunkSparsityPattern &) is not
-   * called with a new sparsity
-   * structure.
+   * You have to make sure that the lifetime of the sparsity structure is at
+   * least as long as that of this matrix or as long as reinit(const
+   * ChunkSparsityPattern &) is not called with a new sparsity structure.
    *
-   * The elements of the matrix are
-   * set to zero by this function.
+   * The elements of the matrix are set to zero by this function.
    */
   virtual void reinit (const ChunkSparsityPattern &sparsity);
 
   /**
-   * Release all memory and return
-   * to a state just like after
-   * having called the default
-   * constructor. It also forgets
-   * the sparsity pattern it was
+   * Release all memory and return to a state just like after having called
+   * the default constructor. It also forgets the sparsity pattern it was
    * previously tied to.
    */
   virtual void clear ();
@@ -276,74 +567,52 @@ public:
    */
 //@{
   /**
-   * Return whether the object is
-   * empty. It is empty if either
-   * both dimensions are zero or no
-   * ChunkSparsityPattern is
-   * associated.
+   * Return whether the object is empty. It is empty if either both dimensions
+   * are zero or no ChunkSparsityPattern is associated.
    */
   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;
 
   /**
-   * 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;
 
   /**
-   * 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;
 
   /**
-   * 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;
 
   /**
-   * 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;
 
   /**
-   * Determine an estimate for the
-   * memory consumption (in bytes)
-   * of this object. See
-   * MemoryConsumption.
+   * Determine an estimate for the memory consumption (in bytes) of this
+   * object. See MemoryConsumption.
    */
   std::size_t memory_consumption () const;
 
@@ -353,47 +622,31 @@ public:
    */
 //@{
   /**
-   * 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,
             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,
             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,
@@ -404,140 +657,84 @@ public:
             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);
 
   /**
-   * Divide the entire matrix by a
-   * fixed factor.
+   * Divide the entire matrix by a fixed factor.
    */
   ChunkSparseMatrix &operator /= (const number factor);
 
   /**
-   * Symmetrize the matrix by
-   * forming the mean value between
-   * the existing matrix and its
-   * transpose, $A = \frac 12(A+A^T)$.
+   * Symmetrize the matrix by forming the mean value between the existing
+   * matrix and its transpose, $A = \frac 12(A+A^T)$.
    *
-   * This operation assumes that
-   * the underlying sparsity
-   * pattern represents a symmetric
-   * object. If this is not the
-   * case, then the result of this
-   * operation will not be a
-   * symmetric matrix, since it
-   * only explicitly symmetrizes
-   * by looping over the lower left
-   * triangular part for efficiency
-   * reasons; if there are entries
-   * in the upper right triangle,
-   * then these elements are missed
-   * in the
-   * symmetrization. Symmetrization
-   * of the sparsity pattern can be
-   * obtain by
-   * ChunkSparsityPattern::symmetrize().
+   * This operation assumes that the underlying sparsity pattern represents a
+   * symmetric object. If this is not the case, then the result of this
+   * operation will not be a symmetric matrix, since it only explicitly
+   * symmetrizes by looping over the lower left triangular part for efficiency
+   * reasons; if there are entries in the upper right triangle, then these
+   * elements are missed in the symmetrization. Symmetrization of the sparsity
+   * pattern can be obtain by ChunkSparsityPattern::symmetrize().
    */
   void symmetrize ();
 
   /**
-   * Copy the given matrix to this
-   * one.  The operation throws an
-   * error if the sparsity patterns
-   * of the two involved matrices
-   * do not point to the same
-   * object, since in this case the
-   * copy operation is
-   * cheaper. Since this operation
-   * is notheless not for free, we
-   * do not make it available
-   * through <tt>operator =</tt>,
-   * since this may lead to
-   * unwanted usage, e.g. in copy
-   * arguments to functions, which
-   * should really be arguments by
-   * reference.
+   * Copy the given matrix to this one.  The operation throws an error if the
+   * sparsity patterns of the two involved matrices do not point to the same
+   * object, since in this case the copy operation is cheaper. Since this
+   * operation is notheless not for free, we do not make it available through
+   * <tt>operator =</tt>, since this may lead to unwanted usage, e.g. in copy
+   * arguments to functions, which should really be arguments by reference.
    *
-   * The source matrix may be a matrix
-   * of arbitrary type, as long as its
-   * data type is convertible to the
-   * data type of this matrix.
+   * The source matrix may be a matrix of arbitrary type, as long as its data
+   * type is convertible to the data type of this matrix.
    *
-   * The function returns a reference to
-   * <tt>*this</tt>.
+   * The function returns a reference to <tt>*this</tt>.
    */
   template <typename somenumber>
   ChunkSparseMatrix<number> &
   copy_from (const ChunkSparseMatrix<somenumber> &source);
 
   /**
-   * This function is complete
-   * analogous to the
-   * ChunkSparsityPattern::copy_from()
-   * function in that it allows to
-   * initialize a whole matrix in
-   * one step. See there for more
-   * information on argument types
-   * and their meaning. You can
-   * also find a small example on
-   * how to use this function
-   * there.
+   * This function is complete analogous to the
+   * ChunkSparsityPattern::copy_from() function in that it allows to
+   * initialize a whole matrix in one step. See there for more information on
+   * argument types and their meaning. You can also find a small example on
+   * how to use this function there.
    *
-   * The only difference to the
-   * cited function is that the
-   * objects which the inner
-   * iterator points to need to be
-   * of type <tt>std::pair<unsigned
-   * int, value</tt>, where
-   * <tt>value</tt> needs to be
-   * convertible to the element
-   * type of this class, as
-   * specified by the
-   * <tt>number</tt> template
+   * The only difference to the cited function is that the objects which the
+   * inner iterator points to need to be of type <tt>std::pair<unsigned int,
+   * value</tt>, where <tt>value</tt> needs to be convertible to the element
+   * type of this class, as specified by the <tt>number</tt> template
    * argument.
    *
-   * Previous content of the matrix
-   * is overwritten. Note that the
-   * entries specified by the input
-   * parameters need not
-   * necessarily cover all elements
-   * of the matrix. Elements not
-   * covered remain untouched.
+   * Previous content of the matrix is overwritten. Note that the entries
+   * specified by the input parameters need not necessarily cover all elements
+   * of the matrix. Elements not covered remain untouched.
    */
   template <typename ForwardIterator>
   void copy_from (const ForwardIterator begin,
                   const ForwardIterator end);
 
   /**
-   * Copy the nonzero entries of a
-   * full matrix into this
-   * object. Previous content is
-   * deleted. Note that the
-   * underlying sparsity pattern
-   * must be appropriate to hold
-   * the nonzero entries of the
-   * full matrix.
+   * Copy the nonzero entries of a full matrix into this object. Previous
+   * content is deleted. Note that the underlying sparsity pattern must be
+   * appropriate to hold the nonzero entries of the full matrix.
    */
   template <typename somenumber>
   void copy_from (const FullMatrix<somenumber> &matrix);
 
   /**
-   * Add <tt>matrix</tt> scaled by
-   * <tt>factor</tt> to this matrix,
-   * i.e. the matrix <tt>factor*matrix</tt>
-   * is added to <tt>this</tt>. This
-   * function throws an error if the
-   * sparsity patterns of the two involved
-   * matrices do not point to the same
-   * object, since in this case the
-   * operation is cheaper.
+   * Add <tt>matrix</tt> scaled by <tt>factor</tt> to this matrix, i.e. the
+   * matrix <tt>factor*matrix</tt> is added to <tt>this</tt>. This function
+   * throws an error if the sparsity patterns of the two involved matrices do
+   * not point to the same object, since in this case the operation is
+   * cheaper.
    *
-   * The source matrix may be a sparse
-   * matrix over an arbitrary underlying
-   * scalar type, as long as its data type
-   * is convertible to the data type of
+   * The source matrix may be a sparse matrix over an arbitrary underlying
+   * scalar type, as long as its data type is convertible to the data type of
    * this matrix.
    */
   template <typename somenumber>
@@ -551,54 +748,31 @@ public:
 //@{
 
   /**
-   * 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;
 
   /**
-   * 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,
@@ -611,21 +785,16 @@ public:
    * 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;
 
   /**
-   * 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);
 
@@ -635,166 +804,105 @@ public:
    */
 //@{
   /**
-   * Matrix-vector multiplication:
-   * let <i>dst = M*src</i> with
-   * <i>M</i> being this matrix.
+   * Matrix-vector multiplication: let <i>dst = M*src</i> with <i>M</i> being
+   * this matrix.
    *
-   * Note that while this function can
-   * operate on all vectors that offer
-   * iterator classes, it is only really
-   * effective for objects of type @ref
-   * Vector. For all classes for which
-   * iterating over elements, or random
-   * member access is expensive, this
-   * function is not efficient. In
-   * particular, if you want to multiply
-   * with BlockVector objects, you should
-   * consider using a BlockChunkSparseMatrix as
-   * well.
+   * Note that while this function can operate on all vectors that offer
+   * iterator classes, it is only really effective for objects of type @ref
+   * Vector. For all classes for which iterating over elements, or random
+   * member access is expensive, this function is not efficient. In
+   * particular, if you want to multiply with BlockVector objects, you should
+   * consider using a BlockChunkSparseMatrix as well.
    *
-   * Source and destination must
-   * not be the same vector.
+   * Source and destination must not be the same vector.
    */
   template <class OutVector, class InVector>
   void vmult (OutVector &dst,
               const InVector &src) const;
 
   /**
-   * Matrix-vector multiplication:
-   * let <i>dst = M<sup>T</sup>*src</i> with
-   * <i>M</i> being this
-   * matrix. This function does the
-   * same as vmult() but takes
-   * the transposed matrix.
+   * Matrix-vector multiplication: let <i>dst = M<sup>T</sup>*src</i> with
+   * <i>M</i> being this matrix. This function does the same as vmult() but
+   * takes the transposed matrix.
    *
-   * Note that while this function can
-   * operate on all vectors that offer
-   * iterator classes, it is only really
-   * effective for objects of type @ref
-   * Vector. For all classes for which
-   * iterating over elements, or random
-   * member access is expensive, this
-   * function is not efficient. In
-   * particular, if you want to multiply
-   * with BlockVector objects, you should
-   * consider using a BlockChunkSparseMatrix as
-   * well.
+   * Note that while this function can operate on all vectors that offer
+   * iterator classes, it is only really effective for objects of type @ref
+   * Vector. For all classes for which iterating over elements, or random
+   * member access is expensive, this function is not efficient. In
+   * particular, if you want to multiply with BlockVector objects, you should
+   * consider using a BlockChunkSparseMatrix as well.
    *
-   * Source and destination must
-   * not be the same vector.
+   * Source and destination must not be the same vector.
    */
   template <class OutVector, class InVector>
   void Tvmult (OutVector &dst,
                const InVector &src) const;
 
   /**
-   * Adding Matrix-vector
-   * multiplication. Add
-   * <i>M*src</i> on <i>dst</i>
-   * with <i>M</i> being this
-   * matrix.
+   * Adding Matrix-vector multiplication. Add <i>M*src</i> on <i>dst</i> with
+   * <i>M</i> being this matrix.
    *
-   * Note that while this function can
-   * operate on all vectors that offer
-   * iterator classes, it is only really
-   * effective for objects of type @ref
-   * Vector. For all classes for which
-   * iterating over elements, or random
-   * member access is expensive, this
-   * function is not efficient. In
-   * particular, if you want to multiply
-   * with BlockVector objects, you should
-   * consider using a BlockChunkSparseMatrix as
-   * well.
+   * Note that while this function can operate on all vectors that offer
+   * iterator classes, it is only really effective for objects of type @ref
+   * Vector. For all classes for which iterating over elements, or random
+   * member access is expensive, this function is not efficient. In
+   * particular, if you want to multiply with BlockVector objects, you should
+   * consider using a BlockChunkSparseMatrix as well.
    *
-   * Source and destination must
-   * not be the same vector.
+   * Source and destination must not be the same vector.
    */
   template <class OutVector, class InVector>
   void vmult_add (OutVector &dst,
                   const InVector &src) const;
 
   /**
-   * Adding Matrix-vector
-   * multiplication. Add
-   * <i>M<sup>T</sup>*src</i> to
-   * <i>dst</i> with <i>M</i> being
-   * this matrix. This function
-   * does the same as vmult_add()
-   * but takes the transposed
-   * matrix.
+   * Adding Matrix-vector multiplication. Add <i>M<sup>T</sup>*src</i> to
+   * <i>dst</i> with <i>M</i> being this matrix. This function does the same
+   * as vmult_add() but takes the transposed matrix.
    *
-   * Note that while this function can
-   * operate on all vectors that offer
-   * iterator classes, it is only really
-   * effective for objects of type @ref
-   * Vector. For all classes for which
-   * iterating over elements, or random
-   * member access is expensive, this
-   * function is not efficient. In
-   * particular, if you want to multiply
-   * with BlockVector objects, you should
-   * consider using a BlockChunkSparseMatrix as
-   * well.
+   * Note that while this function can operate on all vectors that offer
+   * iterator classes, it is only really effective for objects of type @ref
+   * Vector. For all classes for which iterating over elements, or random
+   * member access is expensive, this function is not efficient. In
+   * particular, if you want to multiply with BlockVector objects, you should
+   * consider using a BlockChunkSparseMatrix as well.
    *
-   * Source and destination must
-   * not be the same vector.
+   * Source and destination must not be the same vector.
    */
   template <class OutVector, class InVector>
   void Tvmult_add (OutVector &dst,
                    const InVector &src) const;
 
   /**
-   * Return the square of the norm
-   * of the vector $v$ with respect
-   * to the norm induced by this
-   * matrix,
-   * i.e. $\left(v,Mv\right)$. This
-   * is useful, e.g. in the finite
-   * element context, where the
-   * $L_2$ norm of a function
-   * equals the matrix norm with
-   * respect to the mass matrix of
-   * the vector representing the
-   * nodal values of the finite
-   * element function.
+   * Return the square of the norm of the vector $v$ with respect to the norm
+   * induced by this matrix, i.e. $\left(v,Mv\right)$. This is useful, e.g. in
+   * the finite element context, where the $L_2$ norm of a function equals the
+   * matrix norm with respect to the mass matrix of the vector representing
+   * the nodal values of the finite element function.
    *
-   * Obviously, the matrix needs to be
-   * quadratic for this operation, and for
-   * the result to actually be a norm it
-   * also needs to be either real symmetric
-   * or complex hermitian.
+   * Obviously, the matrix needs to be quadratic for this operation, and for
+   * the result to actually be a norm it also needs to be either real
+   * symmetric or complex hermitian.
    *
-   * The underlying template types of both
-   * this matrix and the given vector
-   * should either both be real or
-   * complex-valued, but not mixed, for
-   * this function to make sense.
+   * The underlying template types of both this matrix and the given vector
+   * should either both be real or complex-valued, but not mixed, for this
+   * function to make sense.
    */
   template <typename somenumber>
   somenumber matrix_norm_square (const Vector<somenumber> &v) const;
 
   /**
-   * Compute the matrix scalar
-   * product $\left(u,Mv\right)$.
+   * Compute the matrix scalar product $\left(u,Mv\right)$.
    */
   template <typename somenumber>
   somenumber matrix_scalar_product (const Vector<somenumber> &u,
                                     const Vector<somenumber> &v) const;
   /**
-   * Compute the residual of an
-   * equation <i>Mx=b</i>, where
-   * the residual is defined to be
-   * <i>r=b-Mx</i>. Write the
-   * residual into
-   * <tt>dst</tt>. The
-   * <i>l<sub>2</sub></i> norm of
-   * the residual vector is
-   * returned.
+   * Compute the residual of an equation <i>Mx=b</i>, where the residual is
+   * defined to be <i>r=b-Mx</i>. Write the residual into <tt>dst</tt>. The
+   * <i>l<sub>2</sub></i> norm of the residual vector is returned.
    *
-   * Source <i>x</i> and destination
-   * <i>dst</i> must not be the same
-   * vector.
+   * Source <i>x</i> and destination <i>dst</i> must not be the same vector.
    */
   template <typename somenumber>
   somenumber residual (Vector<somenumber>       &dst,
@@ -808,37 +916,26 @@ public:
 //@{
 
   /**
-   * Return the l1-norm of the matrix, that is
-   * $|M|_1=max_{all columns j}\sum_{all
-   * rows i} |M_ij|$,
-   * (max. sum of columns).
-   * This is the
-   * natural matrix norm that is compatible
-   * to the l1-norm for vectors, i.e.
-   * $|Mv|_1\leq |M|_1 |v|_1$.
-   * (cf. Haemmerlin-Hoffmann : Numerische Mathematik)
+   * Return the l1-norm of the matrix, that is $|M|_1=max_{all columns
+   * j}\sum_{all rows i} |M_ij|$, (max. sum of columns).  This is the natural
+   * matrix norm that is compatible to the l1-norm for vectors, i.e.
+   * $|Mv|_1\leq |M|_1 |v|_1$.  (cf. Haemmerlin-Hoffmann : Numerische
+   * Mathematik)
    */
   real_type l1_norm () const;
 
   /**
-   * Return the linfty-norm of the
-   * matrix, that is
-   * $|M|_infty=max_{all rows i}\sum_{all
-   * columns j} |M_ij|$,
-   * (max. sum of rows).
-   * This is the
-   * natural matrix norm that is compatible
-   * to the linfty-norm of vectors, i.e.
-   * $|Mv|_infty \leq |M|_infty |v|_infty$.
-   * (cf. Haemmerlin-Hoffmann : Numerische Mathematik)
+   * Return the linfty-norm of the matrix, that is $|M|_infty=max_{all rows
+   * i}\sum_{all columns j} |M_ij|$, (max. sum of rows).  This is the natural
+   * matrix norm that is compatible to the linfty-norm of vectors, i.e.
+   * $|Mv|_infty \leq |M|_infty |v|_infty$.  (cf. Haemmerlin-Hoffmann :
+   * Numerische Mathematik)
    */
   real_type linfty_norm () const;
 
   /**
-   * Return the frobenius norm of the
-   * matrix, i.e. the square root of the
-   * sum of squares of all entries in the
-   * matrix.
+   * Return the frobenius norm of the matrix, i.e. the square root of the sum
+   * of squares of all entries in the matrix.
    */
   real_type frobenius_norm () const;
 //@}
@@ -848,14 +945,9 @@ public:
 //@{
 
   /**
-   * Apply the Jacobi
-   * preconditioner, which
-   * multiplies every element of
-   * the <tt>src</tt> vector by the
-   * inverse of the respective
-   * diagonal element and
-   * multiplies the result with the
-   * relaxation factor <tt>omega</tt>.
+   * Apply the Jacobi preconditioner, which multiplies every element of the
+   * <tt>src</tt> vector by the inverse of the respective diagonal element and
+   * multiplies the result with the relaxation factor <tt>omega</tt>.
    */
   template <typename somenumber>
   void precondition_Jacobi (Vector<somenumber>       &dst,
@@ -863,8 +955,7 @@ public:
                             const number              omega = 1.) const;
 
   /**
-   * Apply SSOR preconditioning to
-   * <tt>src</tt>.
+   * Apply SSOR preconditioning to <tt>src</tt>.
    */
   template <typename somenumber>
   void precondition_SSOR (Vector<somenumber>       &dst,
@@ -872,8 +963,7 @@ public:
                           const number              om = 1.) const;
 
   /**
-   * Apply SOR preconditioning
-   * matrix to <tt>src</tt>.
+   * Apply SOR preconditioning matrix to <tt>src</tt>.
    */
   template <typename somenumber>
   void precondition_SOR (Vector<somenumber>       &dst,
@@ -881,9 +971,7 @@ public:
                          const number              om = 1.) const;
 
   /**
-   * Apply transpose SOR
-   * preconditioning matrix to
-   * <tt>src</tt>.
+   * Apply transpose SOR preconditioning matrix to <tt>src</tt>.
    */
   template <typename somenumber>
   void precondition_TSOR (Vector<somenumber>       &dst,
@@ -891,11 +979,8 @@ public:
                           const number              om = 1.) const;
 
   /**
-   * Perform SSOR preconditioning
-   * in-place.  Apply the
-   * preconditioner matrix without
-   * copying to a second vector.
-   * <tt>omega</tt> is the relaxation
+   * Perform SSOR preconditioning in-place.  Apply the preconditioner matrix
+   * without copying to a second vector.  <tt>omega</tt> is the relaxation
    * parameter.
    */
   template <typename somenumber>
@@ -903,18 +988,15 @@ public:
              const number        omega = 1.) const;
 
   /**
-   * Perform an SOR preconditioning
-   * in-place.  <tt>omega</tt> is
-   * the relaxation parameter.
+   * Perform an SOR preconditioning in-place.  <tt>omega</tt> is the
+   * relaxation parameter.
    */
   template <typename somenumber>
   void SOR (Vector<somenumber> &v,
             const number        om = 1.) const;
 
   /**
-   * Perform a transpose SOR
-   * preconditioning in-place.
-   * <tt>omega</tt> is the
+   * Perform a transpose SOR preconditioning in-place.  <tt>omega</tt> is the
    * relaxation parameter.
    */
   template <typename somenumber>
@@ -922,21 +1004,14 @@ public:
              const number        om = 1.) const;
 
   /**
-   * Perform a permuted SOR
-   * preconditioning in-place.
+   * Perform a permuted SOR preconditioning in-place.
    *
-   * The standard SOR method is
-   * applied in the order
-   * prescribed by <tt>permutation</tt>,
-   * that is, first the row
-   * <tt>permutation[0]</tt>, then
-   * <tt>permutation[1]</tt> and so
-   * on. For efficiency reasons,
-   * the permutation as well as its
-   * inverse are required.
+   * The standard SOR method is applied in the order prescribed by
+   * <tt>permutation</tt>, that is, first the row <tt>permutation[0]</tt>,
+   * then <tt>permutation[1]</tt> and so on. For efficiency reasons, the
+   * permutation as well as its inverse are required.
    *
-   * <tt>omega</tt> is the
-   * relaxation parameter.
+   * <tt>omega</tt> is the relaxation parameter.
    */
   template <typename somenumber>
   void PSOR (Vector<somenumber> &v,
@@ -945,24 +1020,14 @@ public:
              const number        om = 1.) const;
 
   /**
-   * Perform a transposed permuted SOR
-   * preconditioning in-place.
+   * Perform a transposed permuted SOR preconditioning in-place.
    *
-   * The transposed SOR method is
-   * applied in the order
-   * prescribed by
-   * <tt>permutation</tt>, that is,
-   * first the row
-   * <tt>permutation[m()-1]</tt>,
-   * then
-   * <tt>permutation[m()-2]</tt>
-   * and so on. For efficiency
-   * reasons, the permutation as
-   * well as its inverse are
-   * required.
+   * The transposed SOR method is applied in the order prescribed by
+   * <tt>permutation</tt>, that is, first the row <tt>permutation[m()-1]</tt>,
+   * then <tt>permutation[m()-2]</tt> and so on. For efficiency reasons, the
+   * permutation as well as its inverse are required.
    *
-   * <tt>omega</tt> is the
-   * relaxation parameter.
+   * <tt>omega</tt> is the relaxation parameter.
    */
   template <typename somenumber>
   void TPSOR (Vector<somenumber> &v,
@@ -971,10 +1036,8 @@ public:
               const number        om = 1.) const;
 
   /**
-   * Do one SOR step on <tt>v</tt>.
-   * Performs a direct SOR step
-   * with right hand side
-   * <tt>b</tt>.
+   * Do one SOR step on <tt>v</tt>.  Performs a direct SOR step with right
+   * hand side <tt>b</tt>.
    */
   template <typename somenumber>
   void SOR_step (Vector<somenumber> &v,
@@ -982,10 +1045,8 @@ public:
                  const number        om = 1.) const;
 
   /**
-   * Do one adjoint SOR step on
-   * <tt>v</tt>.  Performs a direct
-   * TSOR step with right hand side
-   * <tt>b</tt>.
+   * Do one adjoint SOR step on <tt>v</tt>.  Performs a direct TSOR step with
+   * right hand side <tt>b</tt>.
    */
   template <typename somenumber>
   void TSOR_step (Vector<somenumber> &v,
@@ -993,16 +1054,124 @@ public:
                   const number        om = 1.) const;
 
   /**
-   * Do one SSOR step on
-   * <tt>v</tt>.  Performs a direct
-   * SSOR step with right hand side
-   * <tt>b</tt> by performing TSOR
-   * after SOR.
+   * Do one SSOR step on <tt>v</tt>.  Performs a direct SSOR step with right
+   * hand side <tt>b</tt> by performing TSOR after SOR.
    */
   template <typename somenumber>
   void SSOR_step (Vector<somenumber> &v,
                   const Vector<somenumber> &b,
                   const number        om = 1.) const;
+//@}
+  /**
+   * @name Iterators
+   */
+//@{
+
+  /**
+   * STL-like iterator with the first entry of the matrix. This is the version
+   * for constant matrices.
+   *
+   * Note that due to the layout in ChunkSparseMatrix, iterating over matrix
+   * entries is considerably slower than for a sparse matrix, as the iterator
+   * is travels row-by-row, whereas data is stored in chunks of several rows
+   * and columns.
+   */
+  const_iterator begin () const;
+
+  /**
+   * Final iterator. This is the version for constant matrices.
+   *
+   * Note that due to the layout in ChunkSparseMatrix, iterating over matrix
+   * entries is considerably slower than for a sparse matrix, as the iterator
+   * is travels row-by-row, whereas data is stored in chunks of several rows
+   * and columns.
+   */
+  const_iterator end () const;
+
+  /**
+   * STL-like iterator with the first entry of the matrix. This is the version
+   * for non-constant matrices.
+   *
+   * Note that due to the layout in ChunkSparseMatrix, iterating over matrix
+   * entries is considerably slower than for a sparse matrix, as the iterator
+   * is travels row-by-row, whereas data is stored in chunks of several rows
+   * and columns.
+   */
+  iterator begin ();
+
+  /**
+   * Final iterator. This is the version for non-constant matrices.
+   *
+   * Note that due to the layout in ChunkSparseMatrix, iterating over matrix
+   * entries is considerably slower than for a sparse matrix, as the iterator
+   * is travels row-by-row, whereas data is stored in chunks of several rows
+   * and columns.
+   */
+  iterator end ();
+
+  /**
+   * STL-like iterator with the first entry of row <tt>r</tt>. This is the
+   * version for constant matrices.
+   *
+   * Note that if the given row is empty, i.e. does not contain any nonzero
+   * entries, then the iterator returned by this function equals
+   * <tt>end(r)</tt>. Note also that the iterator may not be dereferencable in
+   * that case.
+   *
+   * Note that due to the layout in ChunkSparseMatrix, iterating over matrix
+   * entries is considerably slower than for a sparse matrix, as the iterator
+   * is travels row-by-row, whereas data is stored in chunks of several rows
+   * and columns.
+   */
+  const_iterator begin (const unsigned int r) const;
+
+  /**
+   * Final iterator of row <tt>r</tt>. It points to the first element past the
+   * end of line @p r, or past the end of the entire sparsity pattern. This is
+   * the version for constant matrices.
+   *
+   * Note that the end iterator is not necessarily dereferencable. This is in
+   * particular the case if it is the end iterator for the last row of a
+   * matrix.
+   *
+   * Note that due to the layout in ChunkSparseMatrix, iterating over matrix
+   * entries is considerably slower than for a sparse matrix, as the iterator
+   * is travels row-by-row, whereas data is stored in chunks of several rows
+   * and columns.
+   */
+  const_iterator end (const unsigned int r) const;
+
+  /**
+   * STL-like iterator with the first entry of row <tt>r</tt>. This is the
+   * version for non-constant matrices.
+   *
+   * Note that if the given row is empty, i.e. does not contain any nonzero
+   * entries, then the iterator returned by this function equals
+   * <tt>end(r)</tt>. Note also that the iterator may not be dereferencable in
+   * that case.
+   *
+   * Note that due to the layout in ChunkSparseMatrix, iterating over matrix
+   * entries is considerably slower than for a sparse matrix, as the iterator
+   * is travels row-by-row, whereas data is stored in chunks of several rows
+   * and columns.
+   */
+  iterator begin (const unsigned int r);
+
+  /**
+   * Final iterator of row <tt>r</tt>. It points to the first element past the
+   * end of line @p r, or past the end of the entire sparsity pattern. This is
+   * the version for non-constant matrices.
+   *
+   * Note that the end iterator is not necessarily dereferencable. This is in
+   * particular the case if it is the end iterator for the last row of a
+   * matrix.
+   *
+   * Note that due to the layout in ChunkSparseMatrix, iterating over matrix
+   * entries is considerably slower than for a sparse matrix, as the iterator
+   * is travels row-by-row, whereas data is stored in chunks of several rows
+   * and columns.
+   */
+  iterator end (const unsigned int r);
 //@}
   /**
    * @name Input/Output
@@ -1010,54 +1179,30 @@ public:
 //@{
 
   /**
-   * Print the matrix to the given
-   * stream, using the format
-   * <tt>(line,col) value</tt>,
-   * i.e. one nonzero entry of the
-   * matrix per line.
+   * Print the matrix to the given stream, using the format <tt>(line,col)
+   * value</tt>, i.e. one nonzero entry of the matrix per line.
    */
   void print (std::ostream &out) const;
 
   /**
-   * Print the matrix in the usual
-   * format, i.e. as a matrix and
-   * not as a list of nonzero
-   * elements. For better
-   * readability, elements not in
-   * the matrix are displayed as
-   * empty space, while matrix
-   * elements which are explicitly
-   * set to zero are displayed as
-   * such.
+   * Print the matrix in the usual format, i.e. as a matrix and not as a list
+   * of nonzero elements. For better readability, elements not in the matrix
+   * are displayed as empty space, while matrix elements which are explicitly
+   * set to zero are displayed as such.
    *
-   * The parameters allow for a
-   * flexible setting of the output
-   * format: <tt>precision</tt> and
-   * <tt>scientific</tt> are used
-   * to determine the number
-   * format, where <tt>scientific =
-   * false</tt> means fixed point
-   * notation.  A zero entry for
-   * <tt>width</tt> makes the
-   * function compute a width, but
-   * it may be changed to a
-   * positive value, if output is
-   * crude.
+   * The parameters allow for a flexible setting of the output format:
+   * <tt>precision</tt> and <tt>scientific</tt> are used to determine the
+   * number format, where <tt>scientific = false</tt> means fixed point
+   * notation.  A zero entry for <tt>width</tt> makes the function compute a
+   * width, but it may be changed to a positive value, if output is crude.
    *
-   * Additionally, a character for
-   * an empty value may be
-   * specified.
+   * Additionally, a character for an empty value may be specified.
    *
-   * Finally, the whole matrix can
-   * be multiplied with a common
-   * denominator to produce more
-   * readable output, even
-   * integers.
+   * Finally, the whole matrix can be multiplied with a common denominator to
+   * produce more readable output, even integers.
    *
-   * @attention This function may
-   * produce <b>large</b> amounts
-   * of output if applied to a
-   * large matrix!
+   * @attention This function may produce <b>large</b> amounts of output if
+   * applied to a large matrix!
    */
   void print_formatted (std::ostream       &out,
                         const unsigned int  precision   = 3,
@@ -1067,69 +1212,40 @@ public:
                         const double        denominator = 1.) const;
 
   /**
-   * Print the actual pattern of
-   * the matrix. For each entry
-   * with an absolute value larger
-   * than threshold, a '*' is
-   * printed, a ':' for every value
-   * smaller and a '.' for every
-   * entry not allocated.
+   * Print the actual pattern of the matrix. For each entry with an absolute
+   * value larger than threshold, a '*' is printed, a ':' for every value
+   * smaller and a '.' for every entry not allocated.
    */
   void print_pattern(std::ostream &out,
                      const double threshold = 0.) const;
 
   /**
-   * Write the data of this object
-   * en bloc to a file. This is
-   * done in a binary mode, so the
-   * output is neither readable by
-   * humans nor (probably) by other
-   * computers using a different
-   * operating system of number
-   * format.
+   * Write the data of this object en bloc to a file. This is done in a binary
+   * mode, so the output is neither readable by humans nor (probably) by other
+   * computers using a different operating system or number format.
    *
-   * The purpose of this function
-   * is that you can swap out
-   * matrices and sparsity pattern
-   * if you are short of memory,
-   * want to communicate between
-   * different programs, or allow
-   * objects to be persistent
-   * across different runs of the
-   * program.
+   * The purpose of this function is that you can swap out matrices and
+   * sparsity pattern if you are short of memory, want to communicate between
+   * different programs, or allow objects to be persistent across different
+   * runs of the program.
    */
   void block_write (std::ostream &out) const;
 
   /**
-   * Read data that has previously
-   * been written by block_write()
-   * from a file. This is done
-   * using the inverse operations
-   * to the above function, so it
-   * is reasonably fast because the
-   * bitstream is not interpreted
-   * except for a few numbers up
-   * front.
+   * Read data that has previously been written by block_write() from a
+   * file. This is done using the inverse operations to the above function, so
+   * it is reasonably fast because the bitstream is not interpreted except for
+   * a few numbers up front.
    *
-   * The object is resized on this
-   * operation, and all previous
-   * contents are lost. Note,
-   * however, that no checks are
-   * performed whether new data and
-   * the underlying ChunkSparsityPattern
-   * object fit together. It is
-   * your responsibility to make
-   * sure that the sparsity pattern
-   * and the data to be read match.
+   * The object is resized on this operation, and all previous contents are
+   * lost. Note, however, that no checks are performed whether new data and
+   * the underlying ChunkSparsityPattern object fit together. It is your
+   * responsibility to make sure that the sparsity pattern and the data to be
+   * read match.
    *
-   * A primitive form of error
-   * checking is performed which
-   * will recognize the bluntest
-   * attempts to interpret some
-   * data as a matrix stored
-   * bitwise to a file that wasn't
-   * actually created that way, but
-   * not more.
+   * A primitive form of error checking is performed which will recognize the
+   * bluntest attempts to interpret some data as a matrix stored bitwise to a
+   * file that wasn't actually created that way, but not more.
    */
   void block_read (std::istream &in);
 //@}
@@ -1167,53 +1283,43 @@ public:
   //@}
 private:
   /**
-   * Pointer to the sparsity
-   * pattern used for this
-   * matrix. In order to guarantee
-   * that it is not deleted while
-   * still in use, we subscribe to
-   * it using the SmartPointer
-   * class.
+   * Pointer to the sparsity pattern used for this matrix. In order to
+   * guarantee that it is not deleted while still in use, we subscribe to it
+   * using the SmartPointer class.
    */
   SmartPointer<const ChunkSparsityPattern,ChunkSparseMatrix<number> > cols;
 
   /**
-   * Array of values for all the
-   * nonzero entries. The position
-   * within the matrix, i.e.  the
-   * row and column number for a
-   * given entry can only be
-   * deduced using the sparsity
-   * pattern. The same holds for
-   * the more common operation of
-   * finding an entry by its
-   * coordinates.
+   * Array of values for all the nonzero entries. The position within the
+   * matrix, i.e.  the row and column number for a given entry can only be
+   * deduced using the sparsity pattern. The same holds for the more common
+   * operation of finding an entry by its coordinates.
    */
   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;
 
   /**
-   * 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;
 
-  // 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;
 };
 
 /*@}*/
@@ -1242,6 +1348,17 @@ unsigned int ChunkSparseMatrix<number>::n () const
 
 
 
+template <typename number>
+inline
+const ChunkSparsityPattern &
+ChunkSparseMatrix<number>::get_sparsity_pattern () const
+{
+  Assert (cols != 0, ExcNotInitialized());
+  return *cols;
+}
+
+
+
 template <typename number>
 inline
 unsigned int
@@ -1250,7 +1367,7 @@ ChunkSparseMatrix<number>::compute_location (const unsigned int i,
 {
   const unsigned int chunk_size = cols->get_chunk_size();
   const unsigned int 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;
@@ -1275,10 +1392,8 @@ void ChunkSparseMatrix<number>::set (const unsigned int i,
   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);
   Assert ((index != SparsityPattern::invalid_entry) ||
           (value == 0.),
@@ -1340,14 +1455,10 @@ ChunkSparseMatrix<number>::operator *= (const number factor)
 
   const unsigned int 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()
@@ -1374,14 +1485,10 @@ ChunkSparseMatrix<number>::operator /= (const number factor)
 
   const unsigned int 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()
@@ -1433,31 +1540,7 @@ number ChunkSparseMatrix<number>::diag_element (const unsigned int i) const
   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 unsigned int 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 unsigned int 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();
   return val[cols->sparsity_pattern.rowstart[i/chunk_size]
@@ -1473,6 +1556,7 @@ number &ChunkSparseMatrix<number>::diag_element (const unsigned int i)
 
 template <typename number>
 template <typename ForwardIterator>
+inline
 void
 ChunkSparseMatrix<number>::copy_from (const ForwardIterator begin,
                                       const ForwardIterator end)
@@ -1480,9 +1564,8 @@ ChunkSparseMatrix<number>::copy_from (const ForwardIterator begin,
   Assert (static_cast<unsigned int>(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;
   for (ForwardIterator i=begin; i!=end; ++i, ++row)
@@ -1495,6 +1578,466 @@ ChunkSparseMatrix<number>::copy_from (const ForwardIterator begin,
 }
 
 
+
+//---------------------------------------------------------------------------
+
+
+namespace ChunkSparseMatrixIterators
+{
+  template <typename number>
+  inline
+  Accessor<number,true>::
+  Accessor (const MatrixType   *matrix,
+            const unsigned int  row)
+    :
+    ChunkSparsityPatternIterators::Accessor (&matrix->get_sparsity_pattern(),
+                                             row),
+    matrix (matrix)
+  {}
+
+
+
+  template <typename number>
+  inline
+  Accessor<number,true>::
+  Accessor (const MatrixType *matrix)
+    :
+    ChunkSparsityPatternIterators::Accessor (&matrix->get_sparsity_pattern()),
+    matrix (matrix)
+  {}
+
+
+
+  template <typename number>
+  inline
+  Accessor<number,true>::
+  Accessor (const ChunkSparseMatrixIterators::Accessor<number,false> &a)
+    :
+    ChunkSparsityPatternIterators::Accessor (a),
+    matrix (&a.get_matrix())
+  {}
+
+
+
+  template <typename number>
+  inline
+  number
+  Accessor<number, true>::value () const
+  {
+    const unsigned int chunk_size = matrix->get_sparsity_pattern().get_chunk_size();
+    return matrix->val[reduced_index() * chunk_size * chunk_size
+                       +
+                       chunk_row * chunk_size
+                       +
+                       chunk_col];
+  }
+
+
+
+  template <typename number>
+  inline
+  typename Accessor<number, true>::MatrixType &
+  Accessor<number, true>::get_matrix () const
+  {
+    return *matrix;
+  }
+
+
+
+  template <typename number>
+  inline
+  Accessor<number, false>::Reference::Reference (
+    const Accessor *accessor,
+    const bool)
+    :
+    accessor (accessor)
+  {}
+
+
+  template <typename number>
+  inline
+  Accessor<number, false>::Reference::operator number() const
+  {
+    const unsigned int chunk_size = accessor->matrix->get_sparsity_pattern().get_chunk_size();
+    return accessor->matrix->val[accessor->reduced_index() * chunk_size * chunk_size
+                                 +
+                                 accessor->chunk_row * chunk_size
+                                 +
+                                 accessor->chunk_col];
+  }
+
+
+
+  template <typename number>
+  inline
+  const typename Accessor<number, false>::Reference &
+  Accessor<number, false>::Reference::operator = (const number n) const
+  {
+    const unsigned int chunk_size = accessor->matrix->get_sparsity_pattern().get_chunk_size();
+    accessor->matrix->val[accessor->reduced_index() * chunk_size * chunk_size
+                          +
+                          accessor->chunk_row * chunk_size
+                          +
+                          accessor->chunk_col] = n;
+    return *this;
+  }
+
+
+
+  template <typename number>
+  inline
+  const typename Accessor<number, false>::Reference &
+  Accessor<number, false>::Reference::operator += (const number n) const
+  {
+    const unsigned int chunk_size = accessor->matrix->get_sparsity_pattern().get_chunk_size();
+    accessor->matrix->val[accessor->reduced_index() * chunk_size * chunk_size
+                          +
+                          accessor->chunk_row * chunk_size
+                          +
+                          accessor->chunk_col] += n;
+    return *this;
+  }
+
+
+
+  template <typename number>
+  inline
+  const typename Accessor<number, false>::Reference &
+  Accessor<number, false>::Reference::operator -= (const number n) const
+  {
+    const unsigned int chunk_size = accessor->matrix->get_sparsity_pattern().get_chunk_size();
+    accessor->matrix->val[accessor->reduced_index() * chunk_size * chunk_size
+                          +
+                          accessor->chunk_row * chunk_size
+                          +
+                          accessor->chunk_col] -= n;
+    return *this;
+  }
+
+
+
+  template <typename number>
+  inline
+  const typename Accessor<number, false>::Reference &
+  Accessor<number, false>::Reference::operator *= (const number n) const
+  {
+    const unsigned int chunk_size = accessor->matrix->get_sparsity_pattern().get_chunk_size();
+    accessor->matrix->val[accessor->reduced_index() * chunk_size * chunk_size
+                          +
+                          accessor->chunk_row * chunk_size
+                          +
+                          accessor->chunk_col] *= n;
+    return *this;
+  }
+
+
+
+  template <typename number>
+  inline
+  const typename Accessor<number, false>::Reference &
+  Accessor<number, false>::Reference::operator /= (const number n) const
+  {
+    const unsigned int chunk_size = accessor->matrix->get_sparsity_pattern().get_chunk_size();
+    accessor->matrix->val[accessor->reduced_index() * chunk_size * chunk_size
+                          +
+                          accessor->chunk_row * chunk_size
+                          +
+                          accessor->chunk_col] /= n;
+    return *this;
+  }
+
+
+
+  template <typename number>
+  inline
+  Accessor<number,false>::
+  Accessor (MatrixType         *matrix,
+            const unsigned int  row)
+    :
+    ChunkSparsityPatternIterators::Accessor (&matrix->get_sparsity_pattern(),
+                                             row),
+    matrix (matrix)
+  {}
+
+
+
+  template <typename number>
+  inline
+  Accessor<number,false>::
+  Accessor (MatrixType         *matrix)
+    :
+    ChunkSparsityPatternIterators::Accessor (&matrix->get_sparsity_pattern()),
+    matrix (matrix)
+  {}
+
+
+
+  template <typename number>
+  inline
+  typename Accessor<number, false>::Reference
+  Accessor<number, false>::value() const
+  {
+    return Reference(this,true);
+  }
+
+
+
+
+  template <typename number>
+  inline
+  typename Accessor<number, false>::MatrixType &
+  Accessor<number, false>::get_matrix () const
+  {
+    return *matrix;
+  }
+
+
+
+  template <typename number, bool Constness>
+  inline
+  Iterator<number, Constness>::
+  Iterator (MatrixType        *matrix,
+            const unsigned int row)
+    :
+    accessor(matrix, row)
+  {}
+
+
+
+  template <typename number, bool Constness>
+  inline
+  Iterator<number, Constness>::
+  Iterator (MatrixType *matrix)
+    :
+    accessor(matrix)
+  {}
+
+
+
+  template <typename number, bool Constness>
+  inline
+  Iterator<number, Constness>::
+  Iterator (const ChunkSparseMatrixIterators::Iterator<number,false> &i)
+    :
+    accessor(*i)
+  {}
+
+
+
+  template <typename number, bool Constness>
+  inline
+  Iterator<number, Constness> &
+  Iterator<number,Constness>::operator++ ()
+  {
+    accessor.advance ();
+    return *this;
+  }
+
+
+  template <typename number, bool Constness>
+  inline
+  Iterator<number,Constness>
+  Iterator<number,Constness>::operator++ (int)
+  {
+    const Iterator iter = *this;
+    accessor.advance ();
+    return iter;
+  }
+
+
+  template <typename number, bool Constness>
+  inline
+  const Accessor<number,Constness> &
+  Iterator<number,Constness>::operator* () const
+  {
+    return accessor;
+  }
+
+
+  template <typename number, bool Constness>
+  inline
+  const Accessor<number,Constness> *
+  Iterator<number,Constness>::operator-> () const
+  {
+    return &accessor;
+  }
+
+
+  template <typename number, bool Constness>
+  inline
+  bool
+  Iterator<number,Constness>::
+  operator == (const Iterator &other) const
+  {
+    return (accessor == other.accessor);
+  }
+
+
+  template <typename number, bool Constness>
+  inline
+  bool
+  Iterator<number,Constness>::
+  operator != (const Iterator &other) const
+  {
+    return ! (*this == other);
+  }
+
+
+  template <typename number, bool Constness>
+  inline
+  bool
+  Iterator<number,Constness>::
+  operator < (const Iterator &other) const
+  {
+    Assert (&accessor.get_matrix() == &other.accessor.get_matrix(),
+            ExcInternalError());
+
+    return (accessor < other.accessor);
+  }
+
+
+  template <typename number, bool Constness>
+  inline
+  bool
+  Iterator<number,Constness>::
+  operator > (const Iterator &other) const
+  {
+    return (other < *this);
+  }
+
+
+  template <typename number, bool Constness>
+  inline
+  int
+  Iterator<number,Constness>::
+  operator - (const Iterator &other) const
+  {
+    Assert (&accessor.get_matrix() == &other.accessor.get_matrix(),
+            ExcInternalError());
+
+    // TODO: can be optimized
+    int difference = 0;
+    if (*this < other)
+      {
+        Iterator copy = *this;
+        while (copy < other)
+          {
+            ++copy;
+            --difference;
+          }
+      }
+    else
+      {
+        Iterator copy = other;
+        while (copy < *this)
+          {
+            ++copy;
+            ++difference;
+          }
+      }
+    return difference;
+  }
+
+
+
+  template <typename number, bool Constness>
+  inline
+  Iterator<number,Constness>
+  Iterator<number,Constness>::
+  operator + (const unsigned int n) const
+  {
+    Iterator x = *this;
+    for (unsigned int i=0; i<n; ++i)
+      ++x;
+
+    return x;
+  }
+
+}
+
+
+
+template <typename number>
+inline
+typename ChunkSparseMatrix<number>::const_iterator
+ChunkSparseMatrix<number>::begin () const
+{
+  return const_iterator(this, 0);
+}
+
+
+template <typename number>
+inline
+typename ChunkSparseMatrix<number>::const_iterator
+ChunkSparseMatrix<number>::end () const
+{
+  return const_iterator(this);
+}
+
+
+template <typename number>
+inline
+typename ChunkSparseMatrix<number>::iterator
+ChunkSparseMatrix<number>::begin ()
+{
+  return iterator(this, 0);
+}
+
+
+template <typename number>
+inline
+typename ChunkSparseMatrix<number>::iterator
+ChunkSparseMatrix<number>::end ()
+{
+  return iterator(this);
+}
+
+
+template <typename number>
+inline
+typename ChunkSparseMatrix<number>::const_iterator
+ChunkSparseMatrix<number>::begin (const unsigned int r) const
+{
+  Assert (r<m(), ExcIndexRange(r,0,m()));
+  return const_iterator(this, r);
+}
+
+
+
+template <typename number>
+inline
+typename ChunkSparseMatrix<number>::const_iterator
+ChunkSparseMatrix<number>::end (const unsigned int r) const
+{
+  Assert (r<m(), ExcIndexRange(r,0,m()));
+  return const_iterator(this, r+1);
+}
+
+
+
+template <typename number>
+inline
+typename ChunkSparseMatrix<number>::iterator
+ChunkSparseMatrix<number>::begin (const unsigned int r)
+{
+  Assert (r<m(), ExcIndexRange(r,0,m()));
+  return iterator(this, r);
+}
+
+
+
+template <typename number>
+inline
+typename ChunkSparseMatrix<number>::iterator
+ChunkSparseMatrix<number>::end (const unsigned int r)
+{
+  Assert (r<m(), ExcIndexRange(r,0,m()));
+  return iterator(this, r+1);
+}
+
+
+
+
 #endif // DOXYGEN
 
 
index 3dccb9ae6d212761385daf4dbe2273738ff023f4..bbc1fb99442846f5bbdc8f97869959b9e319d3a7 100644 (file)
@@ -14,6 +14,7 @@
 
 
 #include <deal.II/base/template_constraints.h>
+#include <deal.II/base/parallel.h>
 #include <deal.II/lac/chunk_sparse_matrix.h>
 #include <deal.II/lac/vector.h>
 #include <deal.II/lac/full_matrix.h>
@@ -43,11 +44,9 @@ namespace internal
   namespace ChunkSparseMatrix
   {
     /**
-     * Add the result of multiplying a chunk
-     * of size chunk_size times chunk_size by
-     * a source vector fragment of size
-     * chunk_size to the destination vector
-     * fragment.
+     * Add the result of multiplying a chunk of size chunk_size times
+     * chunk_size by a source vector fragment of size chunk_size to the
+     * destination vector fragment.
      */
     template <typename MatrixIterator,
              typename SrcIterator,
@@ -77,8 +76,7 @@ namespace internal
 
 
     /**
-     * Like the previous function, but
-     * subtract. We need this for computing
+     * Like the previous function, but subtract. We need this for computing
      * the residual.
      */
     template <typename MatrixIterator,
@@ -108,12 +106,9 @@ namespace internal
 
 
     /**
-     * Add the result of multiplying the
-     * transpose of a chunk of size
-     * chunk_size times chunk_size by a
-     * source vector fragment of size
-     * chunk_size to the destination vector
-     * fragment.
+     * Add the result of multiplying the transpose of a chunk of size
+     * chunk_size times chunk_size by a source vector fragment of size
+     * chunk_size to the destination vector fragment.
      */
     template <typename MatrixIterator,
              typename SrcIterator,
@@ -139,8 +134,7 @@ namespace internal
 
 
     /**
-     * Produce the result of the matrix
-     * scalar product $u^TMv$ for an
+     * Produce the result of the matrix scalar product $u^TMv$ for an
      * individual chunk.
      */
     template <typename result_type,
@@ -172,6 +166,107 @@ namespace internal
 
       return result;
     }
+
+
+
+    /**
+     * Perform a vmult_add using the ChunkSparseMatrix data structures, but
+     * only using a subinterval of the matrix rows.
+     *
+     * In the sequential case, this function is called on all rows, in the
+     * parallel case it may be called on a subrange, at the discretion of the
+     * task scheduler.
+     */
+    template <typename number,
+             typename InVector,
+             typename OutVector>
+    void vmult_add_on_subrange (const ChunkSparsityPattern &cols,
+                                const unsigned int  begin_row,
+                                const unsigned int  end_row,
+                                const number       *values,
+                                const std::size_t  *rowstart,
+                                const unsigned int *colnums,
+                                const InVector     &src,
+                                OutVector          &dst)
+    {
+      const unsigned int m = cols.n_rows();
+      const unsigned int n = cols.n_cols();
+      const unsigned int chunk_size = cols.get_chunk_size();
+
+      // loop over all chunks. note that we need to treat the last chunk row
+      // and column differently if they have padding elements
+      const unsigned int n_filled_last_rows = m % chunk_size;
+      const unsigned int n_filled_last_cols = n % chunk_size;
+
+      const unsigned int last_regular_row = n_filled_last_rows > 0 ?
+                                            std::min(m/chunk_size, end_row) : end_row;
+      const unsigned int irregular_col = n/chunk_size;
+
+      typename OutVector::iterator dst_ptr = dst.begin()+chunk_size*begin_row;
+      const number *val_ptr= &values[rowstart[begin_row]*chunk_size*chunk_size];
+      const unsigned int *colnum_ptr = &colnums[rowstart[begin_row]];
+      for (unsigned int chunk_row=begin_row; chunk_row<last_regular_row;
+           ++chunk_row)
+        {
+          const number *const val_end_of_row = &values[rowstart[chunk_row+1] *
+                                                       chunk_size * chunk_size];
+          while (val_ptr != val_end_of_row)
+            {
+              if (*colnum_ptr != irregular_col)
+                chunk_vmult_add (chunk_size,
+                                 val_ptr,
+                                 src.begin() + *colnum_ptr * chunk_size,
+                                 dst_ptr);
+              else
+                // we're at a chunk column that has padding
+                for (unsigned int r=0; r<chunk_size; ++r)
+                  for (unsigned int c=0; c<n_filled_last_cols; ++c)
+                    dst_ptr[r] += (val_ptr[r*chunk_size + c] *
+                                   src(*colnum_ptr * chunk_size + c));
+
+              ++colnum_ptr;
+              val_ptr += chunk_size * chunk_size;
+            }
+
+          dst_ptr += chunk_size;
+        }
+
+      // now deal with last chunk row if necessary
+      if (n_filled_last_rows > 0 && end_row == (m/chunk_size+1))
+        {
+          const unsigned int chunk_row = last_regular_row;
+
+          const number *const val_end_of_row = &values[rowstart[chunk_row+1] *
+                                                       chunk_size * chunk_size];
+          while (val_ptr != val_end_of_row)
+            {
+              if (*colnum_ptr != irregular_col)
+                {
+                  // we're at a chunk row but not column that has padding
+                  for (unsigned int r=0; r<n_filled_last_rows; ++r)
+                    for (unsigned int c=0; c<chunk_size; ++c)
+                      dst_ptr[r]
+                      += (val_ptr[r*chunk_size + c] *
+                          src(*colnum_ptr * chunk_size + c));
+                }
+              else
+                // we're at a chunk row and column that has padding
+                for (unsigned int r=0; r<n_filled_last_rows; ++r)
+                  for (unsigned int c=0; c<n_filled_last_cols; ++c)
+                    dst_ptr[r]
+                    += (val_ptr[r*chunk_size + c] *
+                        src(*colnum_ptr * chunk_size + c));
+
+              ++colnum_ptr;
+              val_ptr += chunk_size * chunk_size;
+            }
+        }
+      Assert(std::size_t(colnum_ptr-&colnums[0]) == rowstart[end_row],
+             ExcInternalError());
+      Assert(std::size_t(val_ptr-&values[0]) ==
+             rowstart[end_row] * chunk_size * chunk_size,
+             ExcInternalError());
+    }
   }
 }
 
@@ -256,6 +351,22 @@ ChunkSparseMatrix<number>::~ChunkSparseMatrix ()
 
 
 
+namespace internal
+{
+  namespace ChunkSparseMatrix
+  {
+    template<typename T>
+    void zero_subrange (const unsigned int begin,
+                        const unsigned int end,
+                        T *dst)
+    {
+      std::memset (dst+begin,0,(end-begin)*sizeof(T));
+    }
+  }
+}
+
+
+
 template <typename number>
 ChunkSparseMatrix<number> &
 ChunkSparseMatrix<number>::operator = (const double d)
@@ -266,14 +377,27 @@ ChunkSparseMatrix<number>::operator = (const double d)
   Assert (cols->sparsity_pattern.compressed || cols->empty(),
           ChunkSparsityPattern::ExcNotCompressed());
 
-  if (val)
-    {
-      const unsigned int chunk_size = cols->get_chunk_size();
-      std::fill_n (val,
-                   cols->sparsity_pattern.n_nonzero_elements() *
-                   chunk_size * chunk_size,
-                   0.);
-    }
+  // do initial zeroing of elements in parallel. Try to achieve a similar
+  // layout as when doing matrix-vector products, as on some NUMA systems, a
+  // memory block is assigned to memory banks where the first access is
+  // generated. For sparse matrices, the first operations is usually the
+  // operator=. The grain size is chosen to reflect the number of rows in
+  // minimum_parallel_grain_size, weighted by the number of nonzero entries
+  // per row on average.
+  const unsigned int matrix_size = cols->sparsity_pattern.n_nonzero_elements()
+                                   * cols->chunk_size * cols->chunk_size;
+  const unsigned int grain_size =
+    internal::SparseMatrix::minimum_parallel_grain_size *
+    (matrix_size+m()) / m();
+  if (matrix_size>grain_size)
+    parallel::apply_to_subranges (0U, matrix_size,
+                                  std_cxx1x::bind(&internal::ChunkSparseMatrix::template
+                                                  zero_subrange<number>,
+                                                  std_cxx1x::_1, std_cxx1x::_2,
+                                                  val),
+                                  grain_size);
+  else if (matrix_size > 0)
+    std::memset (&val[0], 0, matrix_size*sizeof(number));
 
   return *this;
 }
@@ -313,10 +437,8 @@ ChunkSparseMatrix<number>::reinit (const ChunkSparsityPattern &sparsity)
       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() *
                          chunk_size * chunk_size;
@@ -328,12 +450,10 @@ ChunkSparseMatrix<number>::reinit (const ChunkSparsityPattern &sparsity)
       max_len = N;
     }
 
-  // fill with zeros. do not just fill N
-  // elements but all that we allocated to
-  // ensure that also the padding elements
-  // are zero and not left at previous values
-  if (val != 0)
-    std::fill_n (&val[0], max_len, 0);
+  // fill with zeros. do not just fill N elements but all that we allocated to
+  // ensure that also the padding elements are zero and not left at previous
+  // values
+  this->operator=(0.);
 }
 
 
@@ -378,11 +498,9 @@ ChunkSparseMatrix<number>::n_actually_nonzero_elements () const
 {
   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();
   return std::count_if(&val[0],
                        &val[cols->sparsity_pattern.n_nonzero_elements () *
@@ -413,8 +531,7 @@ ChunkSparseMatrix<number>::copy_from (const ChunkSparseMatrix<somenumber> &matri
   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();
   std::copy (&matrix.val[0],
              &matrix.val[cols->sparsity_pattern.n_nonzero_elements()
@@ -453,8 +570,7 @@ ChunkSparseMatrix<number>::add (const number factor,
   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();
   number             *val_ptr    = &val[0];
   const somenumber   *matrix_ptr = &matrix.val[0];
@@ -479,10 +595,8 @@ ChunkSparseMatrix<number>::vmult (OutVector &dst,
 
   Assert (!PointerComparison::equal(&src, &dst), ExcSourceEqualsDestination());
 
-  // set the output vector to zero and then
-  // add to it the contributions of vmults
-  // from individual chunks. this is what
-  // vmult_add does
+  // set the output vector to zero and then add to it the contributions of
+  // vmults from individual chunks. this is what vmult_add does
   dst = 0;
   vmult_add (dst, src);
 }
@@ -509,10 +623,8 @@ ChunkSparseMatrix<number>::Tvmult (OutVector &dst,
 
   Assert (!PointerComparison::equal(&src, &dst), ExcSourceEqualsDestination());
 
-  // set the output vector to zero and then
-  // add to it the contributions of vmults
-  // from individual chunks. this is what
-  // vmult_add does
+  // set the output vector to zero and then add to it the contributions of
+  // vmults from individual chunks. this is what vmult_add does
   dst = 0;
   Tvmult_add (dst, src);
 }
@@ -531,93 +643,18 @@ ChunkSparseMatrix<number>::vmult_add (OutVector &dst,
   Assert(n() == src.size(), ExcDimensionMismatch(n(),src.size()));
 
   Assert (!PointerComparison::equal(&src, &dst), ExcSourceEqualsDestination());
+  parallel::apply_to_subranges (0U, cols->sparsity_pattern.n_rows(),
+                                std_cxx1x::bind (&internal::ChunkSparseMatrix::vmult_add_on_subrange
+                                                 <number,InVector,OutVector>,
+                                                 std_cxx1x::cref(*cols),
+                                                 std_cxx1x::_1, std_cxx1x::_2,
+                                                 val,
+                                                 cols->sparsity_pattern.rowstart,
+                                                 cols->sparsity_pattern.colnums,
+                                                 std_cxx1x::cref(src),
+                                                 std_cxx1x::ref(dst)),
+                                internal::SparseMatrix::minimum_parallel_grain_size/cols->chunk_size+1);
 
-  const unsigned int 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
-  const bool rows_have_padding = (m() % cols->chunk_size != 0),
-             cols_have_padding = (n() % cols->chunk_size != 0);
-
-  const unsigned int n_regular_chunk_rows
-    = (rows_have_padding ?
-       n_chunk_rows-1 :
-       n_chunk_rows);
-
-  const number       *val_ptr    = val;
-  const unsigned int *colnum_ptr = cols->sparsity_pattern.colnums;
-  typename OutVector::iterator dst_ptr = dst.begin();
-
-  for (unsigned int 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
-                                                * cols->chunk_size];
-      while (val_ptr != val_end_of_row)
-        {
-          if ((cols_have_padding == false)
-              ||
-              (*colnum_ptr != cols->sparsity_pattern.n_cols()-1))
-            internal::ChunkSparseMatrix::chunk_vmult_add
-            (cols->chunk_size,
-             val_ptr,
-             src.begin() + *colnum_ptr * cols->chunk_size,
-             dst_ptr);
-          else
-            // 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)
-                dst(chunk_row * cols->chunk_size + r)
-                += (val_ptr[r*cols->chunk_size + c] *
-                    src(*colnum_ptr * cols->chunk_size + c));
-
-          ++colnum_ptr;
-          val_ptr += cols->chunk_size * cols->chunk_size;
-        }
-
-
-      dst_ptr += cols->chunk_size;
-    }
-
-  // now deal with last chunk row if
-  // necessary
-  if (rows_have_padding)
-    {
-      const unsigned int chunk_row = n_chunk_rows - 1;
-
-      const number *const val_end_of_row = &val[cols->sparsity_pattern.rowstart[chunk_row+1]
-                                                * cols->chunk_size
-                                                * cols->chunk_size];
-      while (val_ptr != val_end_of_row)
-        {
-          if ((cols_have_padding == false)
-              ||
-              (*colnum_ptr != cols->sparsity_pattern.n_cols()-1))
-            {
-              // 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)
-                  dst(chunk_row * cols->chunk_size + r)
-                  += (val_ptr[r*cols->chunk_size + c] *
-                      src(*colnum_ptr * cols->chunk_size + c));
-            }
-          else
-            // 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)
-                dst(chunk_row * cols->chunk_size + r)
-                += (val_ptr[r*cols->chunk_size + c] *
-                    src(*colnum_ptr * cols->chunk_size + c));
-
-          ++colnum_ptr;
-          val_ptr += cols->chunk_size * cols->chunk_size;
-        }
-    }
 }
 
 
@@ -636,10 +673,8 @@ ChunkSparseMatrix<number>::Tvmult_add (OutVector &dst,
 
   const unsigned int 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);
 
@@ -648,9 +683,8 @@ ChunkSparseMatrix<number>::Tvmult_add (OutVector &dst,
        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;
 
@@ -670,8 +704,7 @@ ChunkSparseMatrix<number>::Tvmult_add (OutVector &dst,
              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)
                 dst(*colnum_ptr * cols->chunk_size + c)
@@ -683,8 +716,7 @@ ChunkSparseMatrix<number>::Tvmult_add (OutVector &dst,
         }
     }
 
-  // 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;
@@ -698,8 +730,7 @@ ChunkSparseMatrix<number>::Tvmult_add (OutVector &dst,
               ||
               (*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)
                   dst(*colnum_ptr * cols->chunk_size + c)
@@ -707,8 +738,7 @@ ChunkSparseMatrix<number>::Tvmult_add (OutVector &dst,
                       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)
                 dst(*colnum_ptr * cols->chunk_size + c)
@@ -735,15 +765,12 @@ ChunkSparseMatrix<number>::matrix_norm_square (const Vector<somenumber> &v) cons
   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();
 
-  // 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);
 
@@ -774,8 +801,7 @@ ChunkSparseMatrix<number>::matrix_norm_square (const Vector<somenumber> &v) cons
                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)
                 result
@@ -792,8 +818,7 @@ ChunkSparseMatrix<number>::matrix_norm_square (const Vector<somenumber> &v) cons
       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;
@@ -807,8 +832,7 @@ ChunkSparseMatrix<number>::matrix_norm_square (const Vector<somenumber> &v) cons
               ||
               (*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)
                   result
@@ -818,8 +842,7 @@ ChunkSparseMatrix<number>::matrix_norm_square (const Vector<somenumber> &v) cons
                        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)
                 result
@@ -849,16 +872,13 @@ ChunkSparseMatrix<number>::matrix_scalar_product (const Vector<somenumber> &u,
   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();
 
-  // 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);
 
@@ -889,8 +909,7 @@ ChunkSparseMatrix<number>::matrix_scalar_product (const Vector<somenumber> &u,
                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)
                 result
@@ -907,8 +926,7 @@ ChunkSparseMatrix<number>::matrix_scalar_product (const Vector<somenumber> &u,
       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;
@@ -922,8 +940,7 @@ ChunkSparseMatrix<number>::matrix_scalar_product (const Vector<somenumber> &u,
               ||
               (*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)
                   result
@@ -933,8 +950,7 @@ ChunkSparseMatrix<number>::matrix_scalar_product (const Vector<somenumber> &u,
                        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)
                 result
@@ -962,11 +978,9 @@ ChunkSparseMatrix<number>::l1_norm () const
 
   const unsigned int 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);
 
@@ -994,19 +1008,15 @@ ChunkSparseMatrix<number>::linfty_norm () const
   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();
 
-  // 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);
 
@@ -1030,12 +1040,10 @@ template <typename number>
 typename ChunkSparseMatrix<number>::real_type
 ChunkSparseMatrix<number>::frobenius_norm () const
 {
-  // simply add up all entries in the
-  // sparsity pattern, without taking any
+  // simply add up all entries in the sparsity pattern, without taking any
   // reference to rows or columns
   //
-  // padding elements are zero, so we can add
-  // them up as well
+  // padding elements are zero, so we can add them up as well
   real_type norm_sqr = 0;
   for (const number *ptr = &val[0]; ptr != &val[max_len]; ++ptr)
     norm_sqr +=  numbers::NumberTraits<number>::abs_square(*ptr);
@@ -1060,25 +1068,19 @@ ChunkSparseMatrix<number>::residual (Vector<somenumber>       &dst,
 
   Assert (&u != &dst, ExcSourceEqualsDestination());
 
-  // set dst=b, then subtract the result of
-  // A*u from it. since the purpose of the
-  // current class is to promote streaming of
-  // data rather than more random access
-  // patterns, breaking things up into two
-  // loops may be reasonable
+  // set dst=b, then subtract the result of A*u from it. since the purpose of
+  // the current class is to promote streaming of data rather than more random
+  // access patterns, breaking things up into two loops may be reasonable
   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();
 
-  // 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);
 
@@ -1107,8 +1109,7 @@ ChunkSparseMatrix<number>::residual (Vector<somenumber>       &dst,
              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)
                 dst(chunk_row * cols->chunk_size + r)
@@ -1123,8 +1124,7 @@ ChunkSparseMatrix<number>::residual (Vector<somenumber>       &dst,
       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;
@@ -1138,8 +1138,7 @@ ChunkSparseMatrix<number>::residual (Vector<somenumber>       &dst,
               ||
               (*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)
                   dst(chunk_row * cols->chunk_size + r)
@@ -1147,8 +1146,7 @@ ChunkSparseMatrix<number>::residual (Vector<somenumber>       &dst,
                       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)
                 dst(chunk_row * cols->chunk_size + r)
@@ -1195,10 +1193,8 @@ ChunkSparseMatrix<number>::precondition_SSOR (Vector<somenumber>       &dst,
                                               const Vector<somenumber> &src,
                                               const number              /*om*/) const
 {
-  // to understand how this function works
-  // you may want to take a look at the CVS
-  // archives to see the original version
-  // which is much clearer...
+  // to understand how this function works you may want to take a look at the
+  // CVS archives to see the original version which is much clearer...
   Assert (cols != 0, ExcNotInitialized());
   Assert (val != 0, ExcNotInitialized());
   Assert (m() == n(), ExcMessage("This operation is only valid on square matrices."));
@@ -1387,16 +1383,6 @@ ChunkSparseMatrix<number>::SSOR (Vector<somenumber> &dst,
 
 
 
-template <typename number>
-const ChunkSparsityPattern &
-ChunkSparseMatrix<number>::get_sparsity_pattern () const
-{
-  Assert (cols != 0, ExcNotInitialized());
-  return *cols;
-}
-
-
-
 template <typename number>
 void ChunkSparseMatrix<number>::print (std::ostream &out) const
 {
@@ -1474,9 +1460,8 @@ void ChunkSparseMatrix<number>::print_pattern (std::ostream &out,
 
   const unsigned int 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 (unsigned int d=0; d<chunk_size; ++d)
@@ -1510,8 +1495,7 @@ ChunkSparseMatrix<number>::block_write (std::ostream &out) const
 {
   AssertThrow (out, ExcIO());
 
-  // first the simple objects,
-  // bracketed in [...]
+  // first the simple objects, bracketed in [...]
   out << '[' << max_len << "][";
   // then write out real data
   out.write (reinterpret_cast<const char *>(&val[0]),
index 948c8d017e989727add68074e4d172f4d8e91d8a..655c0ba406483a2e7304e1745cc19db39f20844b 100644 (file)
@@ -34,6 +34,187 @@ template <typename> class ChunkSparseMatrix;
  */
 
 
+
+/**
+ * Iterators on sparsity patterns
+ */
+namespace ChunkSparsityPatternIterators
+{
+  // forward declaration
+  class Iterator;
+
+  /**
+   * Accessor class for iterators into sparsity patterns. This class is
+   * also the base class for both const and non-const accessor classes
+   * into sparse matrices.
+   *
+   * Note that this class only allows read access to elements, providing
+   * their row and column number. It does not allow modifying the
+   * sparsity pattern itself.
+   *
+   * @author Martin Kronbichler
+   * @date 2013
+   */
+  class Accessor
+  {
+  public:
+    /**
+     * Constructor.
+     */
+    Accessor (const ChunkSparsityPattern *matrix,
+              const unsigned int          row);
+
+    /**
+     * Constructor. Construct the end accessor for the given sparsity pattern.
+     */
+    Accessor (const ChunkSparsityPattern *matrix);
+
+    /**
+     * Row number of the element represented by this object. This function can
+     * only be called for entries for which is_valid_entry() is true.
+     */
+    unsigned int row () const;
+
+    /**
+     * Index in row of the element represented by this object. This function
+     * can only be called for entries for which is_valid_entry() is true.
+     */
+    unsigned int index () const;
+
+    /**
+     * Returns the global index from the reduced sparsity pattern.
+     */
+    std::size_t reduced_index() const;
+
+    /**
+     * Column number of the element represented by this object. This function
+     * can only be called for entries for which is_valid_entry() is true.
+     */
+    unsigned int column () const;
+
+    /**
+     * Return whether the sparsity pattern entry pointed to by this iterator
+     * is valid or not. Note that after compressing the sparsity pattern, all
+     * entries are valid. However, before compression, the sparsity pattern
+     * allocated some memory to be used while still adding new nonzero
+     * entries; if you create iterators in this phase of the sparsity
+     * pattern's lifetime, you will iterate over elements that are not
+     * valid. If this is so, then this function will return false.
+     */
+    bool is_valid_entry () const;
+
+
+    /**
+     * Comparison. True, if both iterators point to the same matrix position.
+     */
+    bool operator == (const Accessor &) const;
+
+
+    /**
+     * Comparison operator. Result is true if either the first row number is
+     * smaller or if the row numbers are equal and the first index is smaller.
+     *
+     * This function is only valid if both iterators point into the same
+     * sparsity pattern.
+     */
+    bool operator < (const Accessor &) const;
+
+  protected:
+    /**
+     * The sparsity pattern we operate on accessed.
+     */
+    const ChunkSparsityPattern *sparsity_pattern;
+
+    /**
+     * The accessor of the (reduced) sparsity pattern.
+     */
+    SparsityPatternIterators::Accessor reduced_accessor;
+
+    /**
+     * Current row number.
+     */
+    unsigned int chunk_row;
+
+    /**
+     * Current index in row.
+     */
+    unsigned int chunk_col;
+
+    /**
+     * Move the accessor to the next nonzero entry in the matrix.
+     */
+    void advance ();
+
+    /**
+     * Grant access to iterator class.
+     */
+    friend class Iterator;
+  };
+
+
+
+  /**
+   * STL conforming iterator walking over the elements of a sparsity pattern.
+   */
+  class Iterator
+  {
+  public:
+    /**
+     * Constructor. Create an iterator into the sparsity pattern @p sp for the
+     * given row and the index within it.
+     */
+    Iterator (const ChunkSparsityPattern *sp,
+              const unsigned int          row);
+
+    /**
+     * Prefix increment.
+     */
+    Iterator &operator++ ();
+
+    /**
+     * Postfix increment.
+     */
+    Iterator operator++ (int);
+
+    /**
+     * Dereferencing operator.
+     */
+    const Accessor &operator* () const;
+
+    /**
+     * Dereferencing operator.
+     */
+    const Accessor *operator-> () const;
+
+    /**
+     * Comparison. True, if both iterators point to the same matrix position.
+     */
+    bool operator == (const Iterator &) const;
+
+    /**
+     * Inverse of <tt>==</tt>.
+     */
+    bool operator != (const Iterator &) const;
+
+    /**
+     * Comparison operator. Result is true if either the first row number is
+     * smaller or if the row numbers are equal and the first index is smaller.
+     *
+     * This function is only valid if both iterators point into the same
+     * matrix.
+     */
+    bool operator < (const Iterator &) const;
+
+  private:
+    /**
+     * Store an object of the accessor class.
+     */
+    Accessor accessor;
+  };
+}
+
+
+
 /**
  * Structure representing the sparsity pattern of a sparse matrix.
  *
@@ -46,86 +227,67 @@ template <typename> class ChunkSparseMatrix;
 class ChunkSparsityPattern : public Subscriptor
 {
 public:
+  /**
+   * 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;
 
   /**
-   * 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 ();
 
   /**
-   * Copy constructor. This
-   * constructor is only allowed to
-   * be called if the matrix
-   * structure to be copied is
-   * empty. This is so in order to
-   * prevent involuntary copies of
-   * objects for temporaries, which
-   * can use large amounts of
-   * computing time.  However, copy
-   * constructors are needed if yo
-   * want to use the STL data types
-   * on classes like this, e.g. to
-   * write such statements like
-   * <tt>v.push_back
-   * (ChunkSparsityPattern());</tt>,
-   * with <tt>v</tt> a vector of
-   * ChunkSparsityPattern objects.
-   *
-   * Usually, it is sufficient to
-   * use the explicit keyword to
-   * disallow unwanted temporaries,
-   * but for the STL vectors, this
-   * does not work. Since copying a
-   * structure like this is not
-   * useful anyway because multiple
-   * matrices can use the same
-   * sparsity structure, copies are
-   * only allowed for empty
-   * objects, as described above.
+   * Copy constructor. This constructor is only allowed to be called if the
+   * matrix structure to be copied is empty. This is so in order to prevent
+   * involuntary copies of objects for temporaries, which can use large
+   * amounts of computing time.  However, copy constructors are needed if yo
+   * want to use the STL data types on classes like this, e.g. to write such
+   * statements like <tt>v.push_back (ChunkSparsityPattern());</tt>, with
+   * <tt>v</tt> a vector of ChunkSparsityPattern objects.
+   *
+   * Usually, it is sufficient to use the explicit keyword to disallow
+   * unwanted temporaries, but for the STL vectors, this does not work. Since
+   * copying a structure like this is not useful anyway because multiple
+   * matrices can use the same sparsity structure, copies are only allowed for
+   * empty objects, as described above.
    */
   ChunkSparsityPattern (const ChunkSparsityPattern &);
 
   /**
-   * Initialize a rectangular
-   * matrix.
+   * Initialize a rectangular matrix.
    *
    * @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,
@@ -143,16 +305,12 @@ public:
                         const bool optimize_diagonal) DEAL_II_DEPRECATED;
 
   /**
-   * Initialize a rectangular
-   * matrix.
+   * Initialize a rectangular matrix.
    *
    * @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,
@@ -170,17 +328,12 @@ public:
                         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,
@@ -190,11 +343,8 @@ public:
    * 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,
@@ -215,25 +365,19 @@ public:
   ~ChunkSparsityPattern ();
 
   /**
-   * Copy operator. For this the
-   * same holds as for the copy
-   * constructor: it is declared,
-   * defined and fine to be called,
-   * but the latter only for empty
+   * Copy operator. For this the same holds as for the copy constructor: it is
+   * declared, defined and fine to be called, but the latter only for empty
    * objects.
    */
   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,
@@ -251,21 +395,14 @@ public:
                const bool optimize_diagonal) DEAL_II_DEPRECATED;
 
   /**
-   * Reallocate memory for a matrix
-   * of size <tt>m x n</tt>. The
-   * number of entries for each row
-   * is taken from the array
-   * <tt>row_lengths</tt> which has to
-   * give this number of each row
-   * <tt>i=1...m</tt>.
+   * Reallocate memory for a matrix of size <tt>m x n</tt>. The number of
+   * entries for each row is taken from the array <tt>row_lengths</tt> which
+   * has to give this number of each row <tt>i=1...m</tt>.
    *
-   * If <tt>m*n==0</tt> all memory is freed,
-   * resulting in a total reinitialization
-   * of the object. If it is nonzero, new
-   * memory is only allocated if the new
-   * size extends the old one. This is done
-   * to save time and to avoid fragmentation
-   * of the heap.
+   * If <tt>m*n==0</tt> all memory is freed, resulting in a total
+   * reinitialization of the object. If it is nonzero, new memory is only
+   * allocated if the new size extends the old one. This is done to save time
+   * and to avoid fragmentation of the heap.
    *
    * If the number of rows equals
    * the number of columns then
@@ -290,8 +427,7 @@ public:
                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,
@@ -309,103 +445,52 @@ public:
                const bool optimize_diagonal) DEAL_II_DEPRECATED;
 
   /**
-   * This function compresses the sparsity
-   * structure that this object represents.
-   * It does so by eliminating unused
-   * entries and sorting the remaining ones
-   * to allow faster access by usage of
-   * binary search algorithms. A special
-   * sorting scheme is used for the
-   * diagonal entry of quadratic matrices,
-   * which is always the first entry of
-   * each row.
+   * This function compresses the sparsity structure that this object
+   * represents.  It does so by eliminating unused entries and sorting the
+   * remaining ones to allow faster access by usage of binary search
+   * algorithms. A special sorting scheme is used for the diagonal entry of
+   * quadratic matrices, which is always the first entry of each row.
    *
-   * The memory which is no more
-   * needed is released.
+   * The memory which is no more needed is released.
    *
-   * SparseMatrix objects require the
-   * ChunkSparsityPattern objects they are
-   * initialized with to be compressed, to
-   * reduce memory requirements.
+   * SparseMatrix objects require the ChunkSparsityPattern objects they are
+   * initialized with to be compressed, to reduce memory requirements.
    */
   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 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:
+   * 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)
@@ -417,25 +502,15 @@ public:
    *                     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 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:
+   * 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)
@@ -451,22 +526,13 @@ public:
    *                   column_indices.end());
    * @endcode
    *
-   * 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.
+   * 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,
@@ -488,12 +554,9 @@ public:
                   const bool optimize_diagonal) DEAL_II_DEPRECATED;
 
   /**
-   * Copy data from an object of type
-   * CompressedSparsityPattern,
-   * CompressedSetSparsityPattern or
-   * CompressedSimpleSparsityPattern.
-   * Previous content of this object is
-   * lost, and the sparsity pattern is in
+   * Copy data from an object of type CompressedSparsityPattern,
+   * CompressedSetSparsityPattern or CompressedSimpleSparsityPattern.
+   * Previous content of this object is lost, and the sparsity pattern is in
    * compressed mode afterwards.
    */
   template <typename SparsityType>
@@ -510,14 +573,10 @@ public:
                   const bool          optimize_diagonal) DEAL_II_DEPRECATED;
 
   /**
-   * Take a full matrix and use its
-   * nonzero entries to generate a
-   * sparse matrix entry pattern
-   * for this object.
+   * Take a full matrix and use its nonzero entries to generate a sparse
+   * matrix entry pattern for this object.
    *
-   * Previous content of this
-   * object is lost, and the
-   * sparsity pattern is in
+   * Previous content of this object is lost, and the sparsity pattern is in
    * compressed mode afterwards.
    */
   template <typename number>
@@ -534,71 +593,79 @@ public:
                   const bool optimize_diagonal) DEAL_II_DEPRECATED;
 
   /**
-   * Return whether the object is empty. It
-   * is empty if no memory is allocated,
-   * which is the same as that both
-   * dimensions are zero.
+   * Set the sparsity pattern of the chunk sparsity pattern to be given by
+   * <tt>chunk_size*chunksize</tt> blocks of the sparsity pattern for chunks
+   * specified. Note that the final number of rows <tt>m</tt> of the sparsity
+   * pattern will be approximately <tt>sparsity_pattern_for_chunks.n_rows() *
+   * chunk_size</tt> (modulo padding elements in the last chunk) and similarly
+   * for the number of columns <tt>n</tt>.
+   *
+   * This is a special initialization option in case you can tell the position
+   * of the chunk already from the beginning without generating the sparsity
+   * pattern using <tt>make_sparsity_pattern</tt> calls. This bypasses the
+   * search for chunks but of course needs to be handled with care in order to
+   * give a correct sparsity pattern.
+   *
+   * Previous content of this object is lost, and the sparsity pattern is in
+   * compressed mode afterwards.
+   */
+  template <typename Sparsity>
+  void create_from (const unsigned int  m,
+                    const unsigned int  n,
+                    const Sparsity     &sparsity_pattern_for_chunks,
+                    const unsigned int  chunk_size,
+                    const bool          optimize_diagonal = true);
+
+  /**
+   * Return whether the object is empty. It is empty if no memory is
+   * allocated, which is the same as that both dimensions are zero.
    */
   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;
 
   /**
-   * 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;
 
   /**
-   * 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);
 
   /**
-   * 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.
    */
-  inline unsigned int n_rows () const;
+  unsigned int 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.
    */
-  inline unsigned int n_cols () const;
+  unsigned int 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;
@@ -609,34 +676,25 @@ public:
   unsigned int row_length (const unsigned int 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;
 
   /**
-   * 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;
 
   /**
-   * Return whether the structure is
-   * compressed or not.
+   * Return whether the structure is compressed or not.
    */
   bool is_compressed () const;
 
@@ -663,94 +721,89 @@ public:
   bool stores_only_added_elements () const;
 
   /**
-   * Write the data of this object
-   * en bloc to a file. This is
-   * done in a binary mode, so the
-   * output is neither readable by
-   * humans nor (probably) by other
-   * computers using a different
-   * operating system of number
-   * format.
-   *
-   * The purpose of this function
-   * is that you can swap out
-   * matrices and sparsity pattern
-   * if you are short of memory,
-   * want to communicate between
-   * different programs, or allow
-   * objects to be persistent
-   * across different runs of the
-   * program.
+   * STL-like iterator with the first entry of the matrix. The resulting
+   * iterator can be used to walk over all nonzero entries of the sparsity
+   * pattern.
+   */
+  iterator begin () const;
+
+  /**
+   * Final iterator.
+   */
+  iterator end () const;
+
+  /**
+   * STL-like iterator with the first entry of row <tt>r</tt>.
+   *
+   * Note that if the given row is empty, i.e. does not contain any nonzero
+   * entries, then the iterator returned by this function equals
+   * <tt>end(r)</tt>. Note also that the iterator may not be dereferencable in
+   * that case.
+   */
+  iterator begin (const unsigned int r) const;
+
+  /**
+   * Final iterator of row <tt>r</tt>. It points to the first element past the
+   * end of line @p r, or past the end of the entire sparsity pattern.
+   *
+   * Note that the end iterator is not necessarily dereferencable. This is in
+   * particular the case if it is the end iterator for the last row of a
+   * matrix.
+   */
+  iterator end (const unsigned int r) const;
+
+  /**
+   * Write the data of this object en bloc to a file. This is done in a binary
+   * mode, so the output is neither readable by humans nor (probably) by other
+   * computers using a different operating system of number format.
+   *
+   * The purpose of this function is that you can swap out matrices and
+   * sparsity pattern if you are short of memory, want to communicate between
+   * different programs, or allow objects to be persistent across different
+   * runs of the program.
    */
   void block_write (std::ostream &out) const;
 
   /**
-   * Read data that has previously
-   * been written by block_write()
-   * from a file. This is done
-   * using the inverse operations
-   * to the above function, so it
-   * is reasonably fast because the
-   * bitstream is not interpreted
-   * except for a few numbers up
-   * front.
-   *
-   * The object is resized on this
-   * operation, and all previous
-   * contents are lost.
-   *
-   * A primitive form of error
-   * checking is performed which
-   * will recognize the bluntest
-   * attempts to interpret some
-   * data as a vector stored
-   * bitwise to a file, but not
-   * more.
+   * Read data that has previously been written by block_write() from a
+   * file. This is done using the inverse operations to the above function, so
+   * it is reasonably fast because the bitstream is not interpreted except for
+   * a few numbers up front.
+   *
+   * The object is resized on this operation, and all previous contents are
+   * lost.
+   *
+   * A primitive form of error checking is performed which will recognize the
+   * bluntest attempts to interpret some data as a vector stored bitwise to a
+   * file, but not more.
    */
   void block_read (std::istream &in);
 
   /**
-   * Print the sparsity of the
-   * matrix. The output consists of
-   * one line per row of the format
-   * <tt>[i,j1,j2,j3,...]</tt>. <i>i</i>
-   * is the row number and
-   * <i>jn</i> are the allocated
-   * columns in this row.
+   * Print the sparsity of the matrix. The output consists of one line per row
+   * of the format <tt>[i,j1,j2,j3,...]</tt>. <i>i</i> is the row number and
+   * <i>jn</i> are the allocated columns in this row.
    */
   void print (std::ostream &out) const;
 
   /**
-   * Print the sparsity of the matrix
-   * in a format that <tt>gnuplot</tt> understands
-   * and which can be used to plot the
-   * sparsity pattern in a graphical
-   * way. The format consists of pairs
-   * <tt>i j</tt> of nonzero elements, each
-   * representing one entry of this
-   * matrix, one per line of the output
-   * file. Indices are counted from
-   * zero on, as usual. Since sparsity
-   * patterns are printed in the same
-   * way as matrices are displayed, we
-   * print the negative of the column
-   * index, which means that the
-   * <tt>(0,0)</tt> element is in the top left
-   * rather than in the bottom left
-   * corner.
-   *
-   * Print the sparsity pattern in
-   * gnuplot by setting the data style
-   * to dots or points and use the
-   * <tt>plot</tt> command.
+   * Print the sparsity of the matrix in a format that <tt>gnuplot</tt>
+   * understands and which can be used to plot the sparsity pattern in a
+   * graphical way. The format consists of pairs <tt>i j</tt> of nonzero
+   * elements, each representing one entry of this matrix, one per line of the
+   * output file. Indices are counted from zero on, as usual. Since sparsity
+   * patterns are printed in the same way as matrices are displayed, we print
+   * the negative of the column index, which means that the <tt>(0,0)</tt>
+   * element is in the top left rather than in the bottom left corner.
+   *
+   * Print the sparsity pattern in gnuplot by setting the data style to dots
+   * or points and use the <tt>plot</tt> command.
    */
   void print_gnuplot (std::ostream &out) const;
 
   /**
-   * Determine an estimate for the
-   * memory consumption (in bytes)
-   * of this object. See
-   * MemoryConsumption.
+   * Determine an estimate for the memory consumption (in bytes) of this
+   * object. See MemoryConsumption.
    */
   std::size_t memory_consumption () const;
 
@@ -822,14 +875,12 @@ public:
   //@}
 private:
   /**
-   * Number of rows that this sparsity
-   * structure shall represent.
+   * Number of rows that this sparsity structure shall represent.
    */
   unsigned int rows;
 
   /**
-   * Number of columns that this sparsity
-   * structure shall represent.
+   * Number of columns that this sparsity structure shall represent.
    */
   unsigned int cols;
 
@@ -839,18 +890,20 @@ private:
   unsigned int 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;
 
   /**
-   * Make all the chunk sparse matrix kinds
-   * friends.
+   * Make all the chunk sparse matrix kinds friends.
    */
   template <typename> friend class ChunkSparseMatrix;
+
+  /**
+   * Make the accessor class a friend.
+   */
+  friend class ChunkSparsityPatternIterators::Accessor;
 };
 
 
@@ -859,6 +912,330 @@ private:
 
 #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
+  unsigned int
+  Accessor::index() const
+  {
+    Assert (is_valid_entry() == true, ExcInvalidIterator());
+
+    return sparsity_pattern->get_chunk_size()*reduced_accessor.index() +
+      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());
+
+    // comparison is a bit messy because of the way ChunkSparsityPattern
+    // stores entry: chunk rows run faster than the indices of the reduced
+    // sparsity pattern, but the accessors should of course compare less based
+    // on the actual row, not the reduced one.
+    if (chunk_row == other.chunk_row)
+      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));
+    else
+      {
+        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 row = reduced_accessor.row(),
+          other_row = other.reduced_accessor.row();
+        return (row < other_row
+                ||
+                (row == other_row
+                 &&
+                 (chunk_row < other.chunk_row
+                  ||
+                  (chunk_row == other.chunk_row
+                   &&
+                   (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());
+
+    ++chunk_col;
+
+    // end of chunk
+    if (chunk_col == chunk_size
+        ||
+        reduced_accessor.column() * chunk_size + chunk_col ==
+        sparsity_pattern->n_cols())
+      {
+        // end of row
+        if (reduced_accessor.index() + 1 ==
+            reduced_accessor.sparsity_pattern->row_length(reduced_accessor.row()))
+          {
+            ++chunk_row;
+            chunk_col = 0;
+            const unsigned int old_reduced_row = reduced_accessor.row();
+
+            // end of matrix
+            if (old_reduced_row * chunk_size + chunk_row ==
+                sparsity_pattern->n_rows())
+              {
+                chunk_row = 0;
+                reduced_accessor =
+                  SparsityPatternIterators::Accessor(&sparsity_pattern->sparsity_pattern,
+                                                     sparsity_pattern->sparsity_pattern.n_nonzero_elements());
+              }
+            // end of chunk rows
+            else if (chunk_row == chunk_size)
+              {
+                reduced_accessor =
+                  *sparsity_pattern->sparsity_pattern.begin(old_reduced_row+1);
+                chunk_row = 0;
+              }
+            else
+              reduced_accessor =
+                *sparsity_pattern->sparsity_pattern.begin(old_reduced_row);
+          }
+        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 ! (*this == other);
+  }
+
+
+  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
@@ -894,6 +1271,7 @@ ChunkSparsityPattern::is_compressed () const
 }
 
 
+
 inline
 bool
 ChunkSparsityPattern::optimize_diagonal () const
@@ -904,6 +1282,7 @@ ChunkSparsityPattern::optimize_diagonal () const
 
 
 template <typename ForwardIterator>
+inline
 void
 ChunkSparsityPattern::copy_from (const unsigned int    n_rows,
                                  const unsigned int    n_cols,
@@ -928,16 +1307,10 @@ ChunkSparsityPattern::copy_from (const unsigned int    n_rows,
   Assert (static_cast<unsigned int>(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;
   row_lengths.reserve(n_rows);
@@ -947,8 +1320,7 @@ ChunkSparsityPattern::copy_from (const unsigned int    n_rows,
                            (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;
   typedef typename std::iterator_traits<ForwardIterator>::value_type::const_iterator inner_iterator;
   for (ForwardIterator i=begin; i!=end; ++i, ++row)
@@ -957,16 +1329,14 @@ ChunkSparsityPattern::copy_from (const unsigned int    n_rows,
       for (inner_iterator j=i->begin(); j!=end_of_row; ++j)
         {
           const unsigned int 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);
         }
     }
 
-  // finally compress
-  // everything. this also sorts the
-  // entries within each row
+  // finally compress everything. this also sorts the entries within each row
   compress ();
 }
 
index 788ed489b96bf15fd0ff88a379d4989593023ab5..5bcb7c047087db0dc7bdb76dbc5eea6271adcc25 100644 (file)
@@ -138,8 +138,7 @@ ChunkSparsityPattern::operator = (const ChunkSparsityPattern &s)
   Assert (s.rows == 0, ExcInvalidConstructorCall());
   Assert (s.cols == 0, ExcInvalidConstructorCall());
 
-  // perform the checks in the underlying
-  // object as well
+  // perform the checks in the underlying object as well
   sparsity_pattern = s.sparsity_pattern;
 
   return *this;
@@ -166,8 +165,7 @@ ChunkSparsityPattern::reinit (const unsigned int m,
 {
   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);
   reinit (m, n, row_lengths, chunk_size);
 }
@@ -201,31 +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;
 
-  // 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)
     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);
@@ -287,8 +284,7 @@ ChunkSparsityPattern::copy_from (const SparsityType &csp,
 {
   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)
@@ -382,6 +378,53 @@ ChunkSparsityPattern::reinit (
 
 
 
+namespace internal
+{
+  namespace
+  {
+    template <typename SP>
+    void copy_sparsity (const SP        &src,
+                        SparsityPattern &dst)
+    {
+      dst.copy_from(src);
+    }
+
+    void copy_sparsity (const SparsityPattern &src,
+                        SparsityPattern       &dst)
+    {
+      dst = src;
+    }
+  }
+}
+
+
+
+template <typename Sparsity>
+void
+ChunkSparsityPattern::create_from
+(const unsigned int m,
+ const unsigned int n,
+ const Sparsity    &sparsity_pattern_for_chunks,
+ const unsigned int chunk_size_in,
+ const bool)
+{
+  Assert (m > (sparsity_pattern_for_chunks.n_rows()-1) * chunk_size_in &&
+          m <= sparsity_pattern_for_chunks.n_rows() * chunk_size_in,
+          ExcMessage("Number of rows m is not compatible with chunk size "
+                     "and number of rows in sparsity pattern for the chunks."));
+  Assert (n > (sparsity_pattern_for_chunks.n_cols()-1) * chunk_size_in &&
+          n <= sparsity_pattern_for_chunks.n_cols() * chunk_size_in,
+          ExcMessage("Number of columns m is not compatible with chunk size "
+                     "and number of columns in sparsity pattern for the chunks."));
+
+  internal::copy_sparsity(sparsity_pattern_for_chunks, sparsity_pattern);
+  chunk_size = chunk_size_in;
+  rows = m;
+  cols = n;
+}
+
+
+
 bool
 ChunkSparsityPattern::empty () const
 {
@@ -422,27 +465,41 @@ ChunkSparsityPattern::exists (const unsigned int i,
 
 
 
-unsigned int
-ChunkSparsityPattern::row_length (const unsigned int i) const
+void
+ChunkSparsityPattern::symmetrize ()
 {
-  Assert (i<rows, ExcIndexRange(i,0,rows));
+  // 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());
 
-  return sparsity_pattern.row_length (i/chunk_size) * chunk_size;
+  sparsity_pattern.symmetrize ();
 }
 
 
 
-void
-ChunkSparsityPattern::symmetrize ()
+unsigned int
+ChunkSparsityPattern::row_length (const unsigned int 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;
+    }
 }
 
 
@@ -457,18 +514,15 @@ ChunkSparsityPattern::n_nonzero_elements () const
             chunk_size *
             chunk_size);
   else
-    // some of the chunks reach beyond the
-    // extent of this matrix. this requires a
-    // somewhat more complicated
-    // computations, in particular if the
+    // some of the chunks reach beyond the extent of this matrix. this
+    // requires a somewhat more complicated computations, in particular if the
     // columns don't align
     {
       if ((n_rows() % chunk_size != 0)
           &&
           (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;
@@ -480,10 +534,9 @@ ChunkSparsityPattern::n_nonzero_elements () const
 
       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;
 
@@ -498,17 +551,13 @@ ChunkSparsityPattern::n_nonzero_elements () const
                else if ((row == sparsity_pattern.n_rows() - 1)
                    &&
                    (p->column() != sparsity_pattern.n_cols() - 1))
-                 // last chunk row, but not
-                 // last chunk column. only a
-                 // smaller number (n_rows %
-                 // chunk_size) of rows
-                 // actually exist
+                 // last chunk row, but not last chunk column. only a smaller
+                 // number (n_rows % chunk_size) of rows actually exist
                  n += (n_rows() % chunk_size) * chunk_size;
                else if ((row != sparsity_pattern.n_rows() - 1)
                    &&
                    (p->column() == sparsity_pattern.n_cols() - 1))
-                 // last chunk column, but
-                 // not row
+                 // last chunk column, but not row
                  n += (n_cols() % chunk_size) * chunk_size;
                else
                  // bottom right chunk
@@ -561,8 +610,7 @@ ChunkSparsityPattern::print_gnuplot (std::ostream &out) const
 
   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];
@@ -575,12 +623,9 @@ ChunkSparsityPattern::print_gnuplot (std::ostream &out) const
           for (unsigned int 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;
@@ -593,13 +638,10 @@ ChunkSparsityPattern::print_gnuplot (std::ostream &out) const
 unsigned int
 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
@@ -625,8 +667,7 @@ ChunkSparsityPattern::block_write (std::ostream &out) const
 {
   AssertThrow (out, ExcIO());
 
-  // first the simple objects,
-  // bracketed in [...]
+  // first the simple objects, bracketed in [...]
   out << '['
       << rows << ' '
       << cols << ' '
@@ -660,8 +701,7 @@ ChunkSparsityPattern::block_read (std::istream &in)
   in >> c;
   AssertThrow (c == '[', ExcIO());
 
-  // then read the underlying sparsity
-  // pattern
+  // then read the underlying sparsity pattern
   sparsity_pattern.block_read (in);
 
   in >> c;
@@ -681,6 +721,10 @@ ChunkSparsityPattern::memory_consumption () const
 
 // 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 bool);
@@ -693,6 +737,34 @@ void ChunkSparsityPattern::copy_from<CompressedSimpleSparsityPattern> (const Com
     const unsigned int,
     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 bool);
diff --git a/tests/bits/chunk_sparse_matrix_11.cc b/tests/bits/chunk_sparse_matrix_11.cc
new file mode 100644 (file)
index 0000000..3018ad1
--- /dev/null
@@ -0,0 +1,101 @@
+//-----------------------  chunk_sparse_matrix_11.cc  ------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2004, 2005, 2008, 2013 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//-----------------------  chunk_sparse_matrix_11.cc  ------------------------
+
+
+// set a few elements in a chunk sparse matrix, read them via its iterator
+// class and print them to the log file
+
+#include "../tests.h"
+#include <deal.II/lac/chunk_sparse_matrix.h>
+#include <fstream>
+#include <iomanip>
+
+
+void test (const unsigned int chunk_size)
+{
+  deallog << "Chunk size = " << chunk_size << std::endl;
+
+  ChunkSparsityPattern sp (5,5,3,chunk_size);
+  for (unsigned int i=0; i<5; ++i)
+    for (unsigned int j=0; j<5; ++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);
+
+                                   // then extract the elements (note that
+                                   // some may be zero or even outside the
+                                   // matrix
+  for (unsigned int i=0; i<m.m(); ++i)
+    {
+      deallog << "row " << i << ": ";
+      for (ChunkSparseMatrix<double>::const_iterator it = m.begin(i);
+           it != m.end(i); ++it)
+        {
+          deallog << "[" << it->column() << ","
+                  << std::setprecision(2) << std::fixed << std::setw (4)
+                  << it->value() << "] ";
+        }
+      deallog << std::endl;
+    }
+}
+
+
+
+int main ()
+{
+  std::ofstream logfile("chunk_sparse_matrix_11/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_11/cmp/generic b/tests/bits/chunk_sparse_matrix_11/cmp/generic
new file mode 100644 (file)
index 0000000..a4f1ebb
--- /dev/null
@@ -0,0 +1,31 @@
+
+DEAL::Chunk size = 1
+DEAL::row 0: [0,0   ] [1,0.50] [4,0.50] 
+DEAL::row 1: [1,0   ] [2,1.50] 
+DEAL::row 2: [2,0   ] [0,0.50] [3,3.50] 
+DEAL::row 3: [3,0   ] [1,2.00] [4,6.50] 
+DEAL::row 4: [4,0   ] [2,4.50] 
+DEAL::Chunk size = 2
+DEAL::row 0: [0,0   ] [1,0.50] [2,0   ] [3,0   ] [4,0.50] 
+DEAL::row 1: [0,0   ] [1,0   ] [2,1.50] [3,0   ] [4,0   ] 
+DEAL::row 2: [2,0   ] [3,3.50] [0,0.50] [1,0   ] [4,0   ] 
+DEAL::row 3: [2,0   ] [3,0   ] [0,0   ] [1,2.00] [4,6.50] 
+DEAL::row 4: [4,0   ] [2,4.50] [3,0   ] 
+DEAL::Chunk size = 4
+DEAL::row 0: [0,0   ] [1,0.50] [2,0   ] [3,0   ] [4,0.50] 
+DEAL::row 1: [0,0   ] [1,0   ] [2,1.50] [3,0   ] [4,0   ] 
+DEAL::row 2: [0,0.50] [1,0   ] [2,0   ] [3,3.50] [4,0   ] 
+DEAL::row 3: [0,0   ] [1,2.00] [2,0   ] [3,0   ] [4,6.50] 
+DEAL::row 4: [4,0   ] [0,0   ] [1,0   ] [2,4.50] [3,0   ] 
+DEAL::Chunk size = 5
+DEAL::row 0: [0,0   ] [1,0.50] [2,0   ] [3,0   ] [4,0.50] 
+DEAL::row 1: [0,0   ] [1,0   ] [2,1.50] [3,0   ] [4,0   ] 
+DEAL::row 2: [0,0.50] [1,0   ] [2,0   ] [3,3.50] [4,0   ] 
+DEAL::row 3: [0,0   ] [1,2.00] [2,0   ] [3,0   ] [4,6.50] 
+DEAL::row 4: [0,0   ] [1,0   ] [2,4.50] [3,0   ] [4,0   ] 
+DEAL::Chunk size = 7
+DEAL::row 0: [0,0   ] [1,0.50] [2,0   ] [3,0   ] [4,0.50] 
+DEAL::row 1: [0,0   ] [1,0   ] [2,1.50] [3,0   ] [4,0   ] 
+DEAL::row 2: [0,0.50] [1,0   ] [2,0   ] [3,3.50] [4,0   ] 
+DEAL::row 3: [0,0   ] [1,2.00] [2,0   ] [3,0   ] [4,6.50] 
+DEAL::row 4: [0,0   ] [1,0   ] [2,4.50] [3,0   ] [4,0   ] 
diff --git a/tests/bits/chunk_sparse_matrix_12.cc b/tests/bits/chunk_sparse_matrix_12.cc
new file mode 100644 (file)
index 0000000..6c8e7bf
--- /dev/null
@@ -0,0 +1,101 @@
+//-----------------------  chunk_sparse_matrix_12.cc  ------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2004, 2005, 2008, 2013 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//-----------------------  chunk_sparse_matrix_12.cc  ------------------------
+
+
+// set a few elements in a non-square chunk sparse matrix, read them via its
+// iterator class and print them to the log file
+
+#include "../tests.h"
+#include <deal.II/lac/chunk_sparse_matrix.h>
+#include <fstream>
+#include <iomanip>
+
+
+void test (const unsigned int chunk_size)
+{
+  deallog << "Chunk size = " << chunk_size << std::endl;
+
+  ChunkSparsityPattern sp (5,6,3,chunk_size);
+  for (unsigned int i=0; i<5; ++i)
+    for (unsigned int j=0; j<6; ++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);
+
+                                   // then extract the elements (note that
+                                   // some may be zero or even outside the
+                                   // matrix
+  for (unsigned int i=0; i<m.m(); ++i)
+    {
+      deallog << "row " << i << ": ";
+      for (ChunkSparseMatrix<double>::const_iterator it = m.begin(i);
+           it != m.end(i); ++it)
+        {
+          deallog << "[" << it->column() << ","
+                  << std::setprecision(2) << std::fixed << std::setw (4)
+                  << it->value() << "] ";
+        }
+      deallog << std::endl;
+    }
+}
+
+
+
+int main ()
+{
+  std::ofstream logfile("chunk_sparse_matrix_12/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_12/cmp/generic b/tests/bits/chunk_sparse_matrix_12/cmp/generic
new file mode 100644 (file)
index 0000000..ca0af65
--- /dev/null
@@ -0,0 +1,31 @@
+
+DEAL::Chunk size = 1
+DEAL::row 0: [1,0.50] [4,0.50] 
+DEAL::row 1: [2,1.50] [5,3.00] 
+DEAL::row 2: [0,0.50] [3,3.50] 
+DEAL::row 3: [1,2.00] [4,6.50] 
+DEAL::row 4: [2,4.50] [5,10.50] 
+DEAL::Chunk size = 2
+DEAL::row 0: [0,0   ] [1,0.50] [2,0   ] [3,0   ] [4,0.50] [5,0   ] 
+DEAL::row 1: [0,0   ] [1,0   ] [2,1.50] [3,0   ] [4,0   ] [5,3.00] 
+DEAL::row 2: [2,0   ] [3,3.50] [0,0.50] [1,0   ] [4,0   ] [5,0   ] 
+DEAL::row 3: [2,0   ] [3,0   ] [0,0   ] [1,2.00] [4,6.50] [5,0   ] 
+DEAL::row 4: [4,0   ] [5,10.50] [2,4.50] [3,0   ] 
+DEAL::Chunk size = 4
+DEAL::row 0: [0,0   ] [1,0.50] [2,0   ] [3,0   ] [4,0.50] [5,0   ] 
+DEAL::row 1: [0,0   ] [1,0   ] [2,1.50] [3,0   ] [4,0   ] [5,3.00] 
+DEAL::row 2: [0,0.50] [1,0   ] [2,0   ] [3,3.50] [4,0   ] [5,0   ] 
+DEAL::row 3: [0,0   ] [1,2.00] [2,0   ] [3,0   ] [4,6.50] [5,0   ] 
+DEAL::row 4: [4,0   ] [5,10.50] [0,0   ] [1,0   ] [2,4.50] [3,0   ] 
+DEAL::Chunk size = 5
+DEAL::row 0: [0,0   ] [1,0.50] [2,0   ] [3,0   ] [4,0.50] [5,0   ] 
+DEAL::row 1: [0,0   ] [1,0   ] [2,1.50] [3,0   ] [4,0   ] [5,3.00] 
+DEAL::row 2: [0,0.50] [1,0   ] [2,0   ] [3,3.50] [4,0   ] [5,0   ] 
+DEAL::row 3: [0,0   ] [1,2.00] [2,0   ] [3,0   ] [4,6.50] [5,0   ] 
+DEAL::row 4: [0,0   ] [1,0   ] [2,4.50] [3,0   ] [4,0   ] [5,10.50] 
+DEAL::Chunk size = 7
+DEAL::row 0: [0,0   ] [1,0.50] [2,0   ] [3,0   ] [4,0.50] [5,0   ] 
+DEAL::row 1: [0,0   ] [1,0   ] [2,1.50] [3,0   ] [4,0   ] [5,3.00] 
+DEAL::row 2: [0,0.50] [1,0   ] [2,0   ] [3,3.50] [4,0   ] [5,0   ] 
+DEAL::row 3: [0,0   ] [1,2.00] [2,0   ] [3,0   ] [4,6.50] [5,0   ] 
+DEAL::row 4: [0,0   ] [1,0   ] [2,4.50] [3,0   ] [4,0   ] [5,10.50] 
diff --git a/tests/bits/chunk_sparse_matrix_13.cc b/tests/bits/chunk_sparse_matrix_13.cc
new file mode 100644 (file)
index 0000000..2b3e9fa
--- /dev/null
@@ -0,0 +1,101 @@
+//-----------------------  chunk_sparse_matrix_13.cc  ------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2004, 2005, 2008, 2013 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//-----------------------  chunk_sparse_matrix_13.cc  ------------------------
+
+
+// set a few elements in a non-square chunk sparse matrix, read them via its
+// iterator class and print them to the log file
+
+#include "../tests.h"
+#include <deal.II/lac/chunk_sparse_matrix.h>
+#include <fstream>
+#include <iomanip>
+
+
+void test (const unsigned int chunk_size)
+{
+  deallog << "Chunk size = " << chunk_size << std::endl;
+
+  ChunkSparsityPattern sp (6,5,3,chunk_size);
+  for (unsigned int i=0; i<6; ++i)
+    for (unsigned int j=0; j<5; ++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);
+
+                                   // then extract the elements (note that
+                                   // some may be zero or even outside the
+                                   // matrix
+  for (unsigned int i=0; i<m.m(); ++i)
+    {
+      deallog << "row " << i << ": ";
+      for (ChunkSparseMatrix<double>::const_iterator it = m.begin(i);
+           it != m.end(i); ++it)
+        {
+          deallog << "[" << it->column() << ","
+                  << std::setprecision(2) << std::fixed << std::setw (4)
+                  << it->value() << "] ";
+        }
+      deallog << std::endl;
+    }
+}
+
+
+
+int main ()
+{
+  std::ofstream logfile("chunk_sparse_matrix_13/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_13/cmp/generic b/tests/bits/chunk_sparse_matrix_13/cmp/generic
new file mode 100644 (file)
index 0000000..1f0bdfe
--- /dev/null
@@ -0,0 +1,36 @@
+
+DEAL::Chunk size = 1
+DEAL::row 0: [1,0.50] [4,0.50] 
+DEAL::row 1: [2,1.50] 
+DEAL::row 2: [0,0.50] [3,3.50] 
+DEAL::row 3: [1,2.00] [4,6.50] 
+DEAL::row 4: [2,4.50] 
+DEAL::row 5: [0,0.50] [3,8.00] 
+DEAL::Chunk size = 2
+DEAL::row 0: [0,0   ] [1,0.50] [2,0   ] [3,0   ] [4,0.50] 
+DEAL::row 1: [0,0   ] [1,0   ] [2,1.50] [3,0   ] [4,0   ] 
+DEAL::row 2: [2,0   ] [3,3.50] [0,0.50] [1,0   ] [4,0   ] 
+DEAL::row 3: [2,0   ] [3,0   ] [0,0   ] [1,2.00] [4,6.50] 
+DEAL::row 4: [4,0   ] [0,0   ] [1,0   ] [2,4.50] [3,0   ] 
+DEAL::row 5: [4,0   ] [0,0.50] [1,0   ] [2,0   ] [3,8.00] 
+DEAL::Chunk size = 4
+DEAL::row 0: [0,0   ] [1,0.50] [2,0   ] [3,0   ] [4,0.50] 
+DEAL::row 1: [0,0   ] [1,0   ] [2,1.50] [3,0   ] [4,0   ] 
+DEAL::row 2: [0,0.50] [1,0   ] [2,0   ] [3,3.50] [4,0   ] 
+DEAL::row 3: [0,0   ] [1,2.00] [2,0   ] [3,0   ] [4,6.50] 
+DEAL::row 4: [4,0   ] [0,0   ] [1,0   ] [2,4.50] [3,0   ] 
+DEAL::row 5: [4,0   ] [0,0.50] [1,0   ] [2,0   ] [3,8.00] 
+DEAL::Chunk size = 5
+DEAL::row 0: [0,0   ] [1,0.50] [2,0   ] [3,0   ] [4,0.50] 
+DEAL::row 1: [0,0   ] [1,0   ] [2,1.50] [3,0   ] [4,0   ] 
+DEAL::row 2: [0,0.50] [1,0   ] [2,0   ] [3,3.50] [4,0   ] 
+DEAL::row 3: [0,0   ] [1,2.00] [2,0   ] [3,0   ] [4,6.50] 
+DEAL::row 4: [0,0   ] [1,0   ] [2,4.50] [3,0   ] [4,0   ] 
+DEAL::row 5: [0,0.50] [1,0   ] [2,0   ] [3,8.00] [4,0   ] 
+DEAL::Chunk size = 7
+DEAL::row 0: [0,0   ] [1,0.50] [2,0   ] [3,0   ] [4,0.50] 
+DEAL::row 1: [0,0   ] [1,0   ] [2,1.50] [3,0   ] [4,0   ] 
+DEAL::row 2: [0,0.50] [1,0   ] [2,0   ] [3,3.50] [4,0   ] 
+DEAL::row 3: [0,0   ] [1,2.00] [2,0   ] [3,0   ] [4,6.50] 
+DEAL::row 4: [0,0   ] [1,0   ] [2,4.50] [3,0   ] [4,0   ] 
+DEAL::row 5: [0,0.50] [1,0   ] [2,0   ] [3,8.00] [4,0   ] 
diff --git a/tests/bits/chunk_sparse_matrix_14.cc b/tests/bits/chunk_sparse_matrix_14.cc
new file mode 100644 (file)
index 0000000..0c9f262
--- /dev/null
@@ -0,0 +1,102 @@
+//-----------------------  chunk_sparse_matrix_14.cc  ------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2004, 2005, 2008, 2013 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//-----------------------  chunk_sparse_matrix_14.cc  ------------------------
+
+
+// set a few elements in a chunk sparse matrix and test for iterator
+// inequality
+
+#include "../tests.h"
+#include <deal.II/lac/chunk_sparse_matrix.h>
+#include <fstream>
+#include <iomanip>
+
+
+void test (const unsigned int chunk_size)
+{
+  deallog << "Chunk size = " << chunk_size << std::endl;
+
+  ChunkSparsityPattern sp (5,5,3,chunk_size);
+  for (unsigned int i=0; i<5; ++i)
+    for (unsigned int j=0; j<5; ++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);
+
+                                   // then extract the elements (note that
+                                   // some may be zero or even outside the
+                                   // matrix
+  AssertDimension(m.end()-m.begin(), m.n_nonzero_elements());
+  for (unsigned int i=0; i<m.m(); ++i)
+    {
+      deallog << "row " << i << ": ";
+      AssertDimension(m.end(i)-m.begin(i),
+                      m.get_sparsity_pattern().row_length(i));
+      for (ChunkSparseMatrix<double>::const_iterator it = m.begin(i);
+           it != m.end(i); ++it)
+        {
+          deallog << " done " << (it-m.begin(i)) << ", left " << (it-m.end(i));
+        }
+      deallog << std::endl;
+    }
+}
+
+
+
+int main ()
+{
+  std::ofstream logfile("chunk_sparse_matrix_14/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_14/cmp/generic b/tests/bits/chunk_sparse_matrix_14/cmp/generic
new file mode 100644 (file)
index 0000000..e6cf488
--- /dev/null
@@ -0,0 +1,31 @@
+
+DEAL::Chunk size = 1
+DEAL::row 0:  done 0, left -3 done 1, left -2 done 2, left -1
+DEAL::row 1:  done 0, left -2 done 1, left -1
+DEAL::row 2:  done 0, left -3 done 1, left -2 done 2, left -1
+DEAL::row 3:  done 0, left -3 done 1, left -2 done 2, left -1
+DEAL::row 4:  done 0, left -2 done 1, left -1
+DEAL::Chunk size = 2
+DEAL::row 0:  done 0, left -5 done 1, left -4 done 2, left -3 done 3, left -2 done 4, left -1
+DEAL::row 1:  done 0, left -5 done 1, left -4 done 2, left -3 done 3, left -2 done 4, left -1
+DEAL::row 2:  done 0, left -5 done 1, left -4 done 2, left -3 done 3, left -2 done 4, left -1
+DEAL::row 3:  done 0, left -5 done 1, left -4 done 2, left -3 done 3, left -2 done 4, left -1
+DEAL::row 4:  done 0, left -3 done 1, left -2 done 2, left -1
+DEAL::Chunk size = 4
+DEAL::row 0:  done 0, left -5 done 1, left -4 done 2, left -3 done 3, left -2 done 4, left -1
+DEAL::row 1:  done 0, left -5 done 1, left -4 done 2, left -3 done 3, left -2 done 4, left -1
+DEAL::row 2:  done 0, left -5 done 1, left -4 done 2, left -3 done 3, left -2 done 4, left -1
+DEAL::row 3:  done 0, left -5 done 1, left -4 done 2, left -3 done 3, left -2 done 4, left -1
+DEAL::row 4:  done 0, left -5 done 1, left -4 done 2, left -3 done 3, left -2 done 4, left -1
+DEAL::Chunk size = 5
+DEAL::row 0:  done 0, left -5 done 1, left -4 done 2, left -3 done 3, left -2 done 4, left -1
+DEAL::row 1:  done 0, left -5 done 1, left -4 done 2, left -3 done 3, left -2 done 4, left -1
+DEAL::row 2:  done 0, left -5 done 1, left -4 done 2, left -3 done 3, left -2 done 4, left -1
+DEAL::row 3:  done 0, left -5 done 1, left -4 done 2, left -3 done 3, left -2 done 4, left -1
+DEAL::row 4:  done 0, left -5 done 1, left -4 done 2, left -3 done 3, left -2 done 4, left -1
+DEAL::Chunk size = 7
+DEAL::row 0:  done 0, left -5 done 1, left -4 done 2, left -3 done 3, left -2 done 4, left -1
+DEAL::row 1:  done 0, left -5 done 1, left -4 done 2, left -3 done 3, left -2 done 4, left -1
+DEAL::row 2:  done 0, left -5 done 1, left -4 done 2, left -3 done 3, left -2 done 4, left -1
+DEAL::row 3:  done 0, left -5 done 1, left -4 done 2, left -3 done 3, left -2 done 4, left -1
+DEAL::row 4:  done 0, left -5 done 1, left -4 done 2, left -3 done 3, left -2 done 4, left -1
diff --git a/tests/bits/chunk_sparse_matrix_15.cc b/tests/bits/chunk_sparse_matrix_15.cc
new file mode 100644 (file)
index 0000000..d07f993
--- /dev/null
@@ -0,0 +1,108 @@
+//-----------------------  chunk_sparse_matrix_15.cc  ------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2004, 2005, 2008, 2013 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//-----------------------  chunk_sparse_matrix_15.cc  ------------------------
+
+
+// set a few elements in a chunk sparse matrix, perform a matrix-vector
+// product through the iterator and compare with a matrix-vector product
+
+#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 iterator 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)rand()/RAND_MAX;
+      for (unsigned int i=0; i<m.m(); ++i)
+        {
+          double sum = 0;
+          for (ChunkSparseMatrix<double>::const_iterator it = m.begin(i);
+               it != m.end(i); ++it)
+            sum += it->value() * src(it->column());
+          dst(i) = sum;
+        }
+      m.vmult(dst_ref, src);
+      dst -= dst_ref;
+      deallog << "Error in matrix-vector product done via iterator: "
+              << dst.linfty_norm() << std::endl;
+    }
+}
+
+
+
+int main ()
+{
+  std::ofstream logfile("chunk_sparse_matrix_15/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_15/cmp/generic b/tests/bits/chunk_sparse_matrix_15/cmp/generic
new file mode 100644 (file)
index 0000000..e8ef5ab
--- /dev/null
@@ -0,0 +1,36 @@
+
+DEAL::Chunk size = 1
+DEAL::n_cols = 4
+DEAL::Error in matrix-vector product done via iterator: 0
+DEAL::n_cols = 5
+DEAL::Error in matrix-vector product done via iterator: 0
+DEAL::n_cols = 6
+DEAL::Error in matrix-vector product done via iterator: 0
+DEAL::Chunk size = 2
+DEAL::n_cols = 4
+DEAL::Error in matrix-vector product done via iterator: 0
+DEAL::n_cols = 5
+DEAL::Error in matrix-vector product done via iterator: 0
+DEAL::n_cols = 6
+DEAL::Error in matrix-vector product done via iterator: 0
+DEAL::Chunk size = 4
+DEAL::n_cols = 4
+DEAL::Error in matrix-vector product done via iterator: 0
+DEAL::n_cols = 5
+DEAL::Error in matrix-vector product done via iterator: 0
+DEAL::n_cols = 6
+DEAL::Error in matrix-vector product done via iterator: 0
+DEAL::Chunk size = 5
+DEAL::n_cols = 4
+DEAL::Error in matrix-vector product done via iterator: 0
+DEAL::n_cols = 5
+DEAL::Error in matrix-vector product done via iterator: 0
+DEAL::n_cols = 6
+DEAL::Error in matrix-vector product done via iterator: 0
+DEAL::Chunk size = 7
+DEAL::n_cols = 4
+DEAL::Error in matrix-vector product done via iterator: 0
+DEAL::n_cols = 5
+DEAL::Error in matrix-vector product done via iterator: 0
+DEAL::n_cols = 6
+DEAL::Error in matrix-vector product done via iterator: 0
diff --git a/tests/bits/chunk_sparse_matrix_iterator_01.cc b/tests/bits/chunk_sparse_matrix_iterator_01.cc
new file mode 100644 (file)
index 0000000..eb2261c
--- /dev/null
@@ -0,0 +1,85 @@
+//----------------------------  chunk_sparse_matrix_iterator_01.cc  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2004, 2005, 2013 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  chunk_sparse_matrix_iterator_01.cc  ---------------------------
+
+
+// ChunkSparseMatrix::const_iterator::operator++(int)
+
+#include "../tests.h"
+#include <deal.II/lac/chunk_sparse_matrix.h>
+#include <fstream>
+#include <iomanip>
+
+
+void test (const unsigned int chunk_size)
+{
+  deallog << "Chunk size: " << chunk_size << std::endl;
+  ChunkSparsityPattern sp (5,5,3,chunk_size);
+  for (unsigned int i=0; i<5; ++i)
+    for (unsigned int j=0; j<5; ++j)
+      if ((i+2*j+1) % 3 == 0)
+        sp.add (i,j);
+  sp.compress ();
+
+  ChunkSparseMatrix<double> m(sp);
+  ChunkSparseMatrix<double>::const_iterator i = m.begin();
+  deallog << i->value() << std::endl;
+  ++i;
+  deallog << i->value() << std::endl;
+  i++;
+  deallog << i->value() << std::endl;
+
+  deallog << "OK" << std::endl;
+}
+
+
+
+int main ()
+{
+  std::ofstream logfile("chunk_sparse_matrix_iterator_01/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_iterator_01/cmp/generic b/tests/bits/chunk_sparse_matrix_iterator_01/cmp/generic
new file mode 100644 (file)
index 0000000..58c65b8
--- /dev/null
@@ -0,0 +1,26 @@
+
+DEAL::Chunk size: 1
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::OK
+DEAL::Chunk size: 2
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::OK
+DEAL::Chunk size: 4
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::OK
+DEAL::Chunk size: 5
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::OK
+DEAL::Chunk size: 7
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::OK
diff --git a/tests/bits/chunk_sparse_matrix_iterator_02.cc b/tests/bits/chunk_sparse_matrix_iterator_02.cc
new file mode 100644 (file)
index 0000000..4a6d5d3
--- /dev/null
@@ -0,0 +1,104 @@
+//----------------------------  chunk_sparse_matrix_iterator_02.cc  ---------------------------
+//    $Id$
+//    Version: $Name$ 
+//
+//    Copyright (C) 2004, 2005, 2013 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  chunk_sparse_matrix_iterator_02.cc  ---------------------------
+
+
+// test setting some elements and reading them back from a const matrix
+// iterator
+
+#include "../tests.h"
+#include <deal.II/lac/chunk_sparse_matrix.h>    
+#include <fstream>
+#include <iomanip>
+
+
+void test (const unsigned int chunk_size)
+{
+  deallog << "Chunk size: " << chunk_size << std::endl;
+  ChunkSparsityPattern sp (5,5,3,chunk_size);
+  for (unsigned int i=0; i<5; ++i)
+    for (unsigned int j=0; j<5; ++j)
+      if (((i+2*j+1) % 3 == 0)
+          ||
+          (i==j))
+        sp.add (i,j);
+  sp.compress ();
+
+  ChunkSparseMatrix<double> m(sp);
+  for (unsigned int i=0; i<5; ++i)
+    for (unsigned int j=0; j<5; ++j)
+      if (((i+2*j+1) % 3 == 0)
+          ||
+          (i==j))
+        m.set (i,j, i*j);
+  
+  ChunkSparseMatrix<double>::const_iterator i = m.begin();
+  for (; i!=m.end(); ++i)
+    {
+      deallog << i->row() << ' ' << i->column() << ' '
+              << i->value() << std::endl;
+      if (((i->row()+2*i->column()+1) % 3 == 0)
+          ||
+          (i->row()==i->column()))
+        {
+          Assert (std::fabs(i->value() - i->row()*i->column()) < 1e-14,
+                  ExcInternalError());
+        }
+      else
+        Assert (i->value() == 0, ExcInternalError());
+    }
+
+  deallog << "OK" << std::endl;
+}
+
+
+
+int main ()
+{
+  std::ofstream logfile("chunk_sparse_matrix_iterator_02/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_iterator_02/cmp/generic b/tests/bits/chunk_sparse_matrix_iterator_02/cmp/generic
new file mode 100644 (file)
index 0000000..65368c9
--- /dev/null
@@ -0,0 +1,122 @@
+
+DEAL::Chunk size: 1
+DEAL::0 0 0
+DEAL::0 1 0
+DEAL::0 4 0
+DEAL::1 1 1.00000
+DEAL::1 2 2.00000
+DEAL::2 2 4.00000
+DEAL::2 0 0
+DEAL::2 3 6.00000
+DEAL::3 3 9.00000
+DEAL::3 1 3.00000
+DEAL::3 4 12.0000
+DEAL::4 4 16.0000
+DEAL::4 2 8.00000
+DEAL::OK
+DEAL::Chunk size: 2
+DEAL::0 0 0
+DEAL::0 1 0
+DEAL::0 2 0
+DEAL::0 3 0
+DEAL::0 4 0
+DEAL::1 0 0
+DEAL::1 1 1.00000
+DEAL::1 2 2.00000
+DEAL::1 3 0
+DEAL::1 4 0
+DEAL::2 2 4.00000
+DEAL::2 3 6.00000
+DEAL::2 0 0
+DEAL::2 1 0
+DEAL::2 4 0
+DEAL::3 2 0
+DEAL::3 3 9.00000
+DEAL::3 0 0
+DEAL::3 1 3.00000
+DEAL::3 4 12.0000
+DEAL::4 4 16.0000
+DEAL::4 2 8.00000
+DEAL::4 3 0
+DEAL::OK
+DEAL::Chunk size: 4
+DEAL::0 0 0
+DEAL::0 1 0
+DEAL::0 2 0
+DEAL::0 3 0
+DEAL::0 4 0
+DEAL::1 0 0
+DEAL::1 1 1.00000
+DEAL::1 2 2.00000
+DEAL::1 3 0
+DEAL::1 4 0
+DEAL::2 0 0
+DEAL::2 1 0
+DEAL::2 2 4.00000
+DEAL::2 3 6.00000
+DEAL::2 4 0
+DEAL::3 0 0
+DEAL::3 1 3.00000
+DEAL::3 2 0
+DEAL::3 3 9.00000
+DEAL::3 4 12.0000
+DEAL::4 4 16.0000
+DEAL::4 0 0
+DEAL::4 1 0
+DEAL::4 2 8.00000
+DEAL::4 3 0
+DEAL::OK
+DEAL::Chunk size: 5
+DEAL::0 0 0
+DEAL::0 1 0
+DEAL::0 2 0
+DEAL::0 3 0
+DEAL::0 4 0
+DEAL::1 0 0
+DEAL::1 1 1.00000
+DEAL::1 2 2.00000
+DEAL::1 3 0
+DEAL::1 4 0
+DEAL::2 0 0
+DEAL::2 1 0
+DEAL::2 2 4.00000
+DEAL::2 3 6.00000
+DEAL::2 4 0
+DEAL::3 0 0
+DEAL::3 1 3.00000
+DEAL::3 2 0
+DEAL::3 3 9.00000
+DEAL::3 4 12.0000
+DEAL::4 0 0
+DEAL::4 1 0
+DEAL::4 2 8.00000
+DEAL::4 3 0
+DEAL::4 4 16.0000
+DEAL::OK
+DEAL::Chunk size: 7
+DEAL::0 0 0
+DEAL::0 1 0
+DEAL::0 2 0
+DEAL::0 3 0
+DEAL::0 4 0
+DEAL::1 0 0
+DEAL::1 1 1.00000
+DEAL::1 2 2.00000
+DEAL::1 3 0
+DEAL::1 4 0
+DEAL::2 0 0
+DEAL::2 1 0
+DEAL::2 2 4.00000
+DEAL::2 3 6.00000
+DEAL::2 4 0
+DEAL::3 0 0
+DEAL::3 1 3.00000
+DEAL::3 2 0
+DEAL::3 3 9.00000
+DEAL::3 4 12.0000
+DEAL::4 0 0
+DEAL::4 1 0
+DEAL::4 2 8.00000
+DEAL::4 3 0
+DEAL::4 4 16.0000
+DEAL::OK
diff --git a/tests/bits/chunk_sparse_matrix_iterator_03.cc b/tests/bits/chunk_sparse_matrix_iterator_03.cc
new file mode 100644 (file)
index 0000000..95a5150
--- /dev/null
@@ -0,0 +1,104 @@
+//----------------------------  chunk_sparse_matrix_iterator_03.cc  ---------------------------
+//    $Id$
+//    Version: $Name$ 
+//
+//    Copyright (C) 2004, 2005, 2013 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  chunk_sparse_matrix_iterator_03.cc  ---------------------------
+
+
+// test setting some elements and reading them back from a non-const chunk
+// matrix iterator
+
+#include "../tests.h"
+#include <deal.II/lac/chunk_sparse_matrix.h>    
+#include <fstream>
+#include <iomanip>
+
+
+void test (const unsigned int chunk_size)
+{
+  deallog << "Chunk size: " << chunk_size << std::endl;
+  ChunkSparsityPattern sp (5,5,3,chunk_size);
+  for (unsigned int i=0; i<5; ++i)
+    for (unsigned int j=0; j<5; ++j)
+      if (((i+2*j+1) % 3 == 0)
+          ||
+          (i==j))
+        sp.add (i,j);
+  sp.compress ();
+
+  ChunkSparseMatrix<double> m(sp);
+  for (unsigned int i=0; i<5; ++i)
+    for (unsigned int j=0; j<5; ++j)
+      if (((i+2*j+1) % 3 == 0)
+          ||
+          (i==j))
+        m.set (i,j, i*j);
+  
+  ChunkSparseMatrix<double>::iterator i = m.begin();
+  for (; i!=m.end(); ++i)
+    {
+      deallog << i->row() << ' ' << i->column() << ' '
+              << i->value() << std::endl;
+      if (((i->row()+2*i->column()+1) % 3 == 0)
+          ||
+          (i->row()==i->column()))
+        {
+          Assert (std::fabs(i->value() - i->row()*i->column()) < 1e-14,
+                  ExcInternalError());
+        }
+      else
+        Assert (i->value() == 0, ExcInternalError());
+    }
+
+  deallog << "OK" << std::endl;
+}
+
+
+
+int main ()
+{
+  std::ofstream logfile("chunk_sparse_matrix_iterator_03/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_iterator_03/cmp/generic b/tests/bits/chunk_sparse_matrix_iterator_03/cmp/generic
new file mode 100644 (file)
index 0000000..401253e
--- /dev/null
@@ -0,0 +1,122 @@
+
+DEAL::Chunk size: 1
+DEAL::0 0 0.00000
+DEAL::0 1 0.00000
+DEAL::0 4 0.00000
+DEAL::1 1 1.00000
+DEAL::1 2 2.00000
+DEAL::2 2 4.00000
+DEAL::2 0 0.00000
+DEAL::2 3 6.00000
+DEAL::3 3 9.00000
+DEAL::3 1 3.00000
+DEAL::3 4 12.0000
+DEAL::4 4 16.0000
+DEAL::4 2 8.00000
+DEAL::OK
+DEAL::Chunk size: 2
+DEAL::0 0 0.00000
+DEAL::0 1 0.00000
+DEAL::0 2 0.00000
+DEAL::0 3 0.00000
+DEAL::0 4 0.00000
+DEAL::1 0 0.00000
+DEAL::1 1 1.00000
+DEAL::1 2 2.00000
+DEAL::1 3 0.00000
+DEAL::1 4 0.00000
+DEAL::2 2 4.00000
+DEAL::2 3 6.00000
+DEAL::2 0 0.00000
+DEAL::2 1 0.00000
+DEAL::2 4 0.00000
+DEAL::3 2 0.00000
+DEAL::3 3 9.00000
+DEAL::3 0 0.00000
+DEAL::3 1 3.00000
+DEAL::3 4 12.0000
+DEAL::4 4 16.0000
+DEAL::4 2 8.00000
+DEAL::4 3 0.00000
+DEAL::OK
+DEAL::Chunk size: 4
+DEAL::0 0 0.00000
+DEAL::0 1 0.00000
+DEAL::0 2 0.00000
+DEAL::0 3 0.00000
+DEAL::0 4 0.00000
+DEAL::1 0 0.00000
+DEAL::1 1 1.00000
+DEAL::1 2 2.00000
+DEAL::1 3 0.00000
+DEAL::1 4 0.00000
+DEAL::2 0 0.00000
+DEAL::2 1 0.00000
+DEAL::2 2 4.00000
+DEAL::2 3 6.00000
+DEAL::2 4 0.00000
+DEAL::3 0 0.00000
+DEAL::3 1 3.00000
+DEAL::3 2 0.00000
+DEAL::3 3 9.00000
+DEAL::3 4 12.0000
+DEAL::4 4 16.0000
+DEAL::4 0 0.00000
+DEAL::4 1 0.00000
+DEAL::4 2 8.00000
+DEAL::4 3 0.00000
+DEAL::OK
+DEAL::Chunk size: 5
+DEAL::0 0 0.00000
+DEAL::0 1 0.00000
+DEAL::0 2 0.00000
+DEAL::0 3 0.00000
+DEAL::0 4 0.00000
+DEAL::1 0 0.00000
+DEAL::1 1 1.00000
+DEAL::1 2 2.00000
+DEAL::1 3 0.00000
+DEAL::1 4 0.00000
+DEAL::2 0 0.00000
+DEAL::2 1 0.00000
+DEAL::2 2 4.00000
+DEAL::2 3 6.00000
+DEAL::2 4 0.00000
+DEAL::3 0 0.00000
+DEAL::3 1 3.00000
+DEAL::3 2 0.00000
+DEAL::3 3 9.00000
+DEAL::3 4 12.0000
+DEAL::4 0 0.00000
+DEAL::4 1 0.00000
+DEAL::4 2 8.00000
+DEAL::4 3 0.00000
+DEAL::4 4 16.0000
+DEAL::OK
+DEAL::Chunk size: 7
+DEAL::0 0 0.00000
+DEAL::0 1 0.00000
+DEAL::0 2 0.00000
+DEAL::0 3 0.00000
+DEAL::0 4 0.00000
+DEAL::1 0 0.00000
+DEAL::1 1 1.00000
+DEAL::1 2 2.00000
+DEAL::1 3 0.00000
+DEAL::1 4 0.00000
+DEAL::2 0 0.00000
+DEAL::2 1 0.00000
+DEAL::2 2 4.00000
+DEAL::2 3 6.00000
+DEAL::2 4 0.00000
+DEAL::3 0 0.00000
+DEAL::3 1 3.00000
+DEAL::3 2 0.00000
+DEAL::3 3 9.00000
+DEAL::3 4 12.0000
+DEAL::4 0 0.00000
+DEAL::4 1 0.00000
+DEAL::4 2 8.00000
+DEAL::4 3 0.00000
+DEAL::4 4 16.0000
+DEAL::OK
diff --git a/tests/bits/chunk_sparse_matrix_iterator_04.cc b/tests/bits/chunk_sparse_matrix_iterator_04.cc
new file mode 100644 (file)
index 0000000..59c87d1
--- /dev/null
@@ -0,0 +1,96 @@
+//----------------------------  chunk_sparse_matrix_iterator_04.cc  ---------------------------
+//    $Id$
+//    Version: $Name$ 
+//
+//    Copyright (C) 2004, 2005, 2013 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  chunk_sparse_matrix_iterator_04.cc  ---------------------------
+
+// test setting some elements using a non-const chunk matrix iterator and
+// operator=, and reading them back through the matrix itself
+
+#include "../tests.h"
+#include <deal.II/lac/chunk_sparse_matrix.h>    
+#include <fstream>
+#include <iomanip>
+
+
+void test (const unsigned int chunk_size)
+{
+  deallog << "Chunk size: " << chunk_size << std::endl;
+  ChunkSparsityPattern sp (5,5,3,chunk_size);
+  for (unsigned int i=0; i<5; ++i)
+    for (unsigned int j=0; j<5; ++j)
+      if (((i+2*j+1) % 3 == 0)
+          ||
+          (i==j))
+        sp.add (i,j);
+  sp.compress ();
+  ChunkSparseMatrix<double> m(sp);
+
+  ChunkSparseMatrix<double>::iterator i = m.begin();
+  for (; i!=m.end(); ++i)
+    i->value() = i->row()*i->column();
+  
+  for (unsigned int i=0; i<5; ++i)
+    for (unsigned int j=0; j<5; ++j)
+      if (((i+2*j+1) % 3 == 0)
+          ||
+          (i==j))
+        {
+          deallog << i << ' ' << j << ' ' << m.el(i,j)
+                  << std::endl;
+          Assert (std::fabs(m.el(i,j)-i*j) < 1e-14,
+                  ExcInternalError());
+        }
+  
+  deallog << "OK" << std::endl;
+}
+
+
+
+int main ()
+{
+  std::ofstream logfile("chunk_sparse_matrix_iterator_04/output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  try
+    {
+      test (1);
+      test (2);
+      test (4);
+      test (5);
+      test (7);
+    }
+  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_iterator_04/cmp/generic b/tests/bits/chunk_sparse_matrix_iterator_04/cmp/generic
new file mode 100644 (file)
index 0000000..5ec97e7
--- /dev/null
@@ -0,0 +1,76 @@
+
+DEAL::Chunk size: 1
+DEAL::0 0 0
+DEAL::0 1 0
+DEAL::0 4 0
+DEAL::1 1 1.00000
+DEAL::1 2 2.00000
+DEAL::2 0 0
+DEAL::2 2 4.00000
+DEAL::2 3 6.00000
+DEAL::3 1 3.00000
+DEAL::3 3 9.00000
+DEAL::3 4 12.0000
+DEAL::4 2 8.00000
+DEAL::4 4 16.0000
+DEAL::OK
+DEAL::Chunk size: 2
+DEAL::0 0 0
+DEAL::0 1 0
+DEAL::0 4 0
+DEAL::1 1 1.00000
+DEAL::1 2 2.00000
+DEAL::2 0 0
+DEAL::2 2 4.00000
+DEAL::2 3 6.00000
+DEAL::3 1 3.00000
+DEAL::3 3 9.00000
+DEAL::3 4 12.0000
+DEAL::4 2 8.00000
+DEAL::4 4 16.0000
+DEAL::OK
+DEAL::Chunk size: 4
+DEAL::0 0 0
+DEAL::0 1 0
+DEAL::0 4 0
+DEAL::1 1 1.00000
+DEAL::1 2 2.00000
+DEAL::2 0 0
+DEAL::2 2 4.00000
+DEAL::2 3 6.00000
+DEAL::3 1 3.00000
+DEAL::3 3 9.00000
+DEAL::3 4 12.0000
+DEAL::4 2 8.00000
+DEAL::4 4 16.0000
+DEAL::OK
+DEAL::Chunk size: 5
+DEAL::0 0 0
+DEAL::0 1 0
+DEAL::0 4 0
+DEAL::1 1 1.00000
+DEAL::1 2 2.00000
+DEAL::2 0 0
+DEAL::2 2 4.00000
+DEAL::2 3 6.00000
+DEAL::3 1 3.00000
+DEAL::3 3 9.00000
+DEAL::3 4 12.0000
+DEAL::4 2 8.00000
+DEAL::4 4 16.0000
+DEAL::OK
+DEAL::Chunk size: 7
+DEAL::0 0 0
+DEAL::0 1 0
+DEAL::0 4 0
+DEAL::1 1 1.00000
+DEAL::1 2 2.00000
+DEAL::2 0 0
+DEAL::2 2 4.00000
+DEAL::2 3 6.00000
+DEAL::3 1 3.00000
+DEAL::3 3 9.00000
+DEAL::3 4 12.0000
+DEAL::4 2 8.00000
+DEAL::4 4 16.0000
+DEAL::OK
diff --git a/tests/bits/chunk_sparse_matrix_iterator_05.cc b/tests/bits/chunk_sparse_matrix_iterator_05.cc
new file mode 100644 (file)
index 0000000..011df54
--- /dev/null
@@ -0,0 +1,103 @@
+//----------------------------  chunk_sparse_matrix_iterator_05.cc  ---------------------------
+//    $Id$
+//    Version: $Name$ 
+//
+//    Copyright (C) 2004, 2005, 2013 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  chunk_sparse_matrix_iterator_05.cc  ---------------------------
+
+// test setting some elements using a non-const chunk matrix iterator and
+// operator+=, and reading them back through the matrix itself
+
+#include "../tests.h"
+#include <deal.II/lac/chunk_sparse_matrix.h>    
+#include <fstream>
+#include <iomanip>
+
+
+void test (const unsigned int chunk_size)
+{
+  deallog << "Chunk size: " << chunk_size << std::endl;
+  ChunkSparsityPattern sp (5,5,3,chunk_size);
+  for (unsigned int i=0; i<5; ++i)
+    for (unsigned int j=0; j<5; ++j)
+      if (((i+2*j+1) % 3 == 0)
+          ||
+          (i==j))
+        sp.add (i,j);
+  sp.compress ();
+
+  ChunkSparseMatrix<double> m(sp);
+  for (unsigned int i=0; i<5; ++i)
+    for (unsigned int j=0; j<5; ++j)
+      if (((i+2*j+1) % 3 == 0)
+          ||
+          (i==j))
+        m.set(i,j,1.);
+
+  ChunkSparseMatrix<double>::iterator i = m.begin();
+  for (; i!=m.end(); ++i)
+    i->value() += i->row()*i->column();
+  
+  for (unsigned int i=0; i<5; ++i)
+    for (unsigned int j=0; j<5; ++j)
+      if (((i+2*j+1) % 3 == 0)
+          ||
+          (i==j))
+        {
+          deallog << i << ' ' << j << ' ' << m.el(i,j)
+                  << std::endl;
+          Assert (std::fabs(m.el(i,j)-(1.+i*j)) < 1e-14,
+                  ExcInternalError());
+        }
+  
+  deallog << "OK" << std::endl;
+}
+
+
+
+int main ()
+{
+  std::ofstream logfile("chunk_sparse_matrix_iterator_05/output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  try
+    {
+      test (1);
+      test (2);
+      test (4);
+      test (5);
+      test (7);
+    }
+  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_iterator_05/cmp/generic b/tests/bits/chunk_sparse_matrix_iterator_05/cmp/generic
new file mode 100644 (file)
index 0000000..9454630
--- /dev/null
@@ -0,0 +1,76 @@
+
+DEAL::Chunk size: 1
+DEAL::0 0 1.00000
+DEAL::0 1 1.00000
+DEAL::0 4 1.00000
+DEAL::1 1 2.00000
+DEAL::1 2 3.00000
+DEAL::2 0 1.00000
+DEAL::2 2 5.00000
+DEAL::2 3 7.00000
+DEAL::3 1 4.00000
+DEAL::3 3 10.0000
+DEAL::3 4 13.0000
+DEAL::4 2 9.00000
+DEAL::4 4 17.0000
+DEAL::OK
+DEAL::Chunk size: 2
+DEAL::0 0 1.00000
+DEAL::0 1 1.00000
+DEAL::0 4 1.00000
+DEAL::1 1 2.00000
+DEAL::1 2 3.00000
+DEAL::2 0 1.00000
+DEAL::2 2 5.00000
+DEAL::2 3 7.00000
+DEAL::3 1 4.00000
+DEAL::3 3 10.0000
+DEAL::3 4 13.0000
+DEAL::4 2 9.00000
+DEAL::4 4 17.0000
+DEAL::OK
+DEAL::Chunk size: 4
+DEAL::0 0 1.00000
+DEAL::0 1 1.00000
+DEAL::0 4 1.00000
+DEAL::1 1 2.00000
+DEAL::1 2 3.00000
+DEAL::2 0 1.00000
+DEAL::2 2 5.00000
+DEAL::2 3 7.00000
+DEAL::3 1 4.00000
+DEAL::3 3 10.0000
+DEAL::3 4 13.0000
+DEAL::4 2 9.00000
+DEAL::4 4 17.0000
+DEAL::OK
+DEAL::Chunk size: 5
+DEAL::0 0 1.00000
+DEAL::0 1 1.00000
+DEAL::0 4 1.00000
+DEAL::1 1 2.00000
+DEAL::1 2 3.00000
+DEAL::2 0 1.00000
+DEAL::2 2 5.00000
+DEAL::2 3 7.00000
+DEAL::3 1 4.00000
+DEAL::3 3 10.0000
+DEAL::3 4 13.0000
+DEAL::4 2 9.00000
+DEAL::4 4 17.0000
+DEAL::OK
+DEAL::Chunk size: 7
+DEAL::0 0 1.00000
+DEAL::0 1 1.00000
+DEAL::0 4 1.00000
+DEAL::1 1 2.00000
+DEAL::1 2 3.00000
+DEAL::2 0 1.00000
+DEAL::2 2 5.00000
+DEAL::2 3 7.00000
+DEAL::3 1 4.00000
+DEAL::3 3 10.0000
+DEAL::3 4 13.0000
+DEAL::4 2 9.00000
+DEAL::4 4 17.0000
+DEAL::OK
diff --git a/tests/bits/chunk_sparse_matrix_iterator_06.cc b/tests/bits/chunk_sparse_matrix_iterator_06.cc
new file mode 100644 (file)
index 0000000..cb9cd38
--- /dev/null
@@ -0,0 +1,110 @@
+//----------------------------  chunk_sparse_matrix_iterator_06.cc  ---------------------------
+//    $Id$
+//    Version: $Name$ 
+//
+//    Copyright (C) 2004, 2005, 2013 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  chunk_sparse_matrix_iterator_06.cc  ---------------------------
+
+// test setting some elements using a non-const chunk matrix iterator and
+// operator-=, and reading them back through the matrix itself
+
+#include "../tests.h"
+#include <deal.II/lac/chunk_sparse_matrix.h>    
+#include <fstream>
+#include <iomanip>
+
+
+void test (const unsigned int chunk_size)
+{
+  deallog << "Chunk size: " << chunk_size << std::endl;
+  ChunkSparsityPattern sp (5,5,3,chunk_size);
+  for (unsigned int i=0; i<5; ++i)
+    for (unsigned int j=0; j<5; ++j)
+      if (((i+2*j+1) % 3 == 0)
+          ||
+          (i==j))
+        sp.add (i,j);
+  sp.compress ();
+
+  ChunkSparseMatrix<double> m(sp);
+  for (unsigned int i=0; i<5; ++i)
+    for (unsigned int j=0; j<5; ++j)
+      if (((i+2*j+1) % 3 == 0)
+          ||
+          (i==j))
+        m.set(i,j,1.);
+
+  ChunkSparseMatrix<double>::iterator i = m.begin();
+  for (; i!=m.end(); ++i)
+    i->value() -= i->row()*i->column();
+  
+  for (unsigned int i=0; i<5; ++i)
+    for (unsigned int j=0; j<5; ++j)
+      if (((i+2*j+1) % 3 == 0)
+          ||
+          (i==j))
+        {
+          deallog << i << ' ' << j << ' ' << m.el(i,j)
+                  << std::endl;
+          Assert (std::fabs(m.el(i,j)-(1.-i*j)) < 1e-14,
+                  ExcInternalError());
+        }
+      else if (sp.exists(i,j))
+        {
+          deallog << i << ' ' << j << ' ' << m.el(i,j)
+                  << std::endl;
+          Assert (std::fabs(m.el(i,j)+i*j) < 1e-14,
+                  ExcInternalError());
+        }
+  
+  deallog << "OK" << std::endl;
+}
+
+
+
+int main ()
+{
+  std::ofstream logfile("chunk_sparse_matrix_iterator_06/output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  try
+    {
+      test (1);
+      test (2);
+      test (4);
+      test (5);
+      test (7);
+    }
+  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_iterator_06/cmp/generic b/tests/bits/chunk_sparse_matrix_iterator_06/cmp/generic
new file mode 100644 (file)
index 0000000..5654a87
--- /dev/null
@@ -0,0 +1,122 @@
+
+DEAL::Chunk size: 1
+DEAL::0 0 1.00000
+DEAL::0 1 1.00000
+DEAL::0 4 1.00000
+DEAL::1 1 0
+DEAL::1 2 -1.00000
+DEAL::2 0 1.00000
+DEAL::2 2 -3.00000
+DEAL::2 3 -5.00000
+DEAL::3 1 -2.00000
+DEAL::3 3 -8.00000
+DEAL::3 4 -11.0000
+DEAL::4 2 -7.00000
+DEAL::4 4 -15.0000
+DEAL::OK
+DEAL::Chunk size: 2
+DEAL::0 0 1.00000
+DEAL::0 1 1.00000
+DEAL::0 2 0
+DEAL::0 3 0
+DEAL::0 4 1.00000
+DEAL::1 0 0
+DEAL::1 1 0
+DEAL::1 2 -1.00000
+DEAL::1 3 -3.00000
+DEAL::1 4 -4.00000
+DEAL::2 0 1.00000
+DEAL::2 1 -2.00000
+DEAL::2 2 -3.00000
+DEAL::2 3 -5.00000
+DEAL::2 4 -8.00000
+DEAL::3 0 0
+DEAL::3 1 -2.00000
+DEAL::3 2 -6.00000
+DEAL::3 3 -8.00000
+DEAL::3 4 -11.0000
+DEAL::4 2 -7.00000
+DEAL::4 3 -12.0000
+DEAL::4 4 -15.0000
+DEAL::OK
+DEAL::Chunk size: 4
+DEAL::0 0 1.00000
+DEAL::0 1 1.00000
+DEAL::0 2 0
+DEAL::0 3 0
+DEAL::0 4 1.00000
+DEAL::1 0 0
+DEAL::1 1 0
+DEAL::1 2 -1.00000
+DEAL::1 3 -3.00000
+DEAL::1 4 -4.00000
+DEAL::2 0 1.00000
+DEAL::2 1 -2.00000
+DEAL::2 2 -3.00000
+DEAL::2 3 -5.00000
+DEAL::2 4 -8.00000
+DEAL::3 0 0
+DEAL::3 1 -2.00000
+DEAL::3 2 -6.00000
+DEAL::3 3 -8.00000
+DEAL::3 4 -11.0000
+DEAL::4 0 0
+DEAL::4 1 -4.00000
+DEAL::4 2 -7.00000
+DEAL::4 3 -12.0000
+DEAL::4 4 -15.0000
+DEAL::OK
+DEAL::Chunk size: 5
+DEAL::0 0 1.00000
+DEAL::0 1 1.00000
+DEAL::0 2 0
+DEAL::0 3 0
+DEAL::0 4 1.00000
+DEAL::1 0 0
+DEAL::1 1 0
+DEAL::1 2 -1.00000
+DEAL::1 3 -3.00000
+DEAL::1 4 -4.00000
+DEAL::2 0 1.00000
+DEAL::2 1 -2.00000
+DEAL::2 2 -3.00000
+DEAL::2 3 -5.00000
+DEAL::2 4 -8.00000
+DEAL::3 0 0
+DEAL::3 1 -2.00000
+DEAL::3 2 -6.00000
+DEAL::3 3 -8.00000
+DEAL::3 4 -11.0000
+DEAL::4 0 0
+DEAL::4 1 -4.00000
+DEAL::4 2 -7.00000
+DEAL::4 3 -12.0000
+DEAL::4 4 -15.0000
+DEAL::OK
+DEAL::Chunk size: 7
+DEAL::0 0 1.00000
+DEAL::0 1 1.00000
+DEAL::0 2 0
+DEAL::0 3 0
+DEAL::0 4 1.00000
+DEAL::1 0 0
+DEAL::1 1 0
+DEAL::1 2 -1.00000
+DEAL::1 3 -3.00000
+DEAL::1 4 -4.00000
+DEAL::2 0 1.00000
+DEAL::2 1 -2.00000
+DEAL::2 2 -3.00000
+DEAL::2 3 -5.00000
+DEAL::2 4 -8.00000
+DEAL::3 0 0
+DEAL::3 1 -2.00000
+DEAL::3 2 -6.00000
+DEAL::3 3 -8.00000
+DEAL::3 4 -11.0000
+DEAL::4 0 0
+DEAL::4 1 -4.00000
+DEAL::4 2 -7.00000
+DEAL::4 3 -12.0000
+DEAL::4 4 -15.0000
+DEAL::OK
diff --git a/tests/bits/chunk_sparse_matrix_iterator_07.cc b/tests/bits/chunk_sparse_matrix_iterator_07.cc
new file mode 100644 (file)
index 0000000..c380b20
--- /dev/null
@@ -0,0 +1,110 @@
+//----------------------------  chunk_sparse_matrix_iterator_07.cc  ---------------------------
+//    $Id$
+//    Version: $Name$ 
+//
+//    Copyright (C) 2004, 2005, 2013 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  chunk_sparse_matrix_iterator_07.cc  ---------------------------
+
+// test setting some elements using a non-const chunk matrix iterator and
+// operator*=, and reading them back through the matrix itself
+
+#include "../tests.h"
+#include <deal.II/lac/chunk_sparse_matrix.h>    
+#include <fstream>
+#include <iomanip>
+
+
+void test (const unsigned int chunk_size)
+{
+  deallog << "Chunk size: " << chunk_size << std::endl;
+  ChunkSparsityPattern sp (5,5,3,chunk_size);
+  for (unsigned int i=0; i<5; ++i)
+    for (unsigned int j=0; j<5; ++j)
+      if (((i+2*j+1) % 3 == 0)
+          ||
+          (i==j))
+        sp.add (i,j);
+  sp.compress ();
+
+  ChunkSparseMatrix<double> m(sp);
+  for (unsigned int i=0; i<5; ++i)
+    for (unsigned int j=0; j<5; ++j)
+      if (((i+2*j+1) % 3 == 0)
+          ||
+          (i==j))
+        m.set(i,j,i*j);
+
+  ChunkSparseMatrix<double>::iterator i = m.begin();
+  for (; i!=m.end(); ++i)
+    i->value() *= 2;
+  
+  for (unsigned int i=0; i<5; ++i)
+    for (unsigned int j=0; j<5; ++j)
+      if (((i+2*j+1) % 3 == 0)
+          ||
+          (i==j))
+        {
+          deallog << i << ' ' << j << ' ' << m.el(i,j)
+                  << std::endl;
+          Assert (std::fabs(m.el(i,j)-2*i*j) < 1e-14,
+                  ExcInternalError());
+        }
+      else if (sp.exists(i,j))
+        {
+          deallog << i << ' ' << j << ' ' << m.el(i,j)
+                  << std::endl;
+          Assert (std::fabs(m.el(i,j)) < 1e-14,
+                  ExcInternalError());
+        }
+  
+  deallog << "OK" << std::endl;
+}
+
+
+
+int main ()
+{
+  std::ofstream logfile("chunk_sparse_matrix_iterator_07/output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  try
+    {
+      test (1);
+      test (2);
+      test (4);
+      test (5);
+      test (7);
+    }
+  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_iterator_07/cmp/generic b/tests/bits/chunk_sparse_matrix_iterator_07/cmp/generic
new file mode 100644 (file)
index 0000000..3b469fb
--- /dev/null
@@ -0,0 +1,122 @@
+
+DEAL::Chunk size: 1
+DEAL::0 0 0
+DEAL::0 1 0
+DEAL::0 4 0
+DEAL::1 1 2.00000
+DEAL::1 2 4.00000
+DEAL::2 0 0
+DEAL::2 2 8.00000
+DEAL::2 3 12.0000
+DEAL::3 1 6.00000
+DEAL::3 3 18.0000
+DEAL::3 4 24.0000
+DEAL::4 2 16.0000
+DEAL::4 4 32.0000
+DEAL::OK
+DEAL::Chunk size: 2
+DEAL::0 0 0
+DEAL::0 1 0
+DEAL::0 2 0
+DEAL::0 3 0
+DEAL::0 4 0
+DEAL::1 0 0
+DEAL::1 1 2.00000
+DEAL::1 2 4.00000
+DEAL::1 3 0
+DEAL::1 4 0
+DEAL::2 0 0
+DEAL::2 1 0
+DEAL::2 2 8.00000
+DEAL::2 3 12.0000
+DEAL::2 4 0
+DEAL::3 0 0
+DEAL::3 1 6.00000
+DEAL::3 2 0
+DEAL::3 3 18.0000
+DEAL::3 4 24.0000
+DEAL::4 2 16.0000
+DEAL::4 3 0
+DEAL::4 4 32.0000
+DEAL::OK
+DEAL::Chunk size: 4
+DEAL::0 0 0
+DEAL::0 1 0
+DEAL::0 2 0
+DEAL::0 3 0
+DEAL::0 4 0
+DEAL::1 0 0
+DEAL::1 1 2.00000
+DEAL::1 2 4.00000
+DEAL::1 3 0
+DEAL::1 4 0
+DEAL::2 0 0
+DEAL::2 1 0
+DEAL::2 2 8.00000
+DEAL::2 3 12.0000
+DEAL::2 4 0
+DEAL::3 0 0
+DEAL::3 1 6.00000
+DEAL::3 2 0
+DEAL::3 3 18.0000
+DEAL::3 4 24.0000
+DEAL::4 0 0
+DEAL::4 1 0
+DEAL::4 2 16.0000
+DEAL::4 3 0
+DEAL::4 4 32.0000
+DEAL::OK
+DEAL::Chunk size: 5
+DEAL::0 0 0
+DEAL::0 1 0
+DEAL::0 2 0
+DEAL::0 3 0
+DEAL::0 4 0
+DEAL::1 0 0
+DEAL::1 1 2.00000
+DEAL::1 2 4.00000
+DEAL::1 3 0
+DEAL::1 4 0
+DEAL::2 0 0
+DEAL::2 1 0
+DEAL::2 2 8.00000
+DEAL::2 3 12.0000
+DEAL::2 4 0
+DEAL::3 0 0
+DEAL::3 1 6.00000
+DEAL::3 2 0
+DEAL::3 3 18.0000
+DEAL::3 4 24.0000
+DEAL::4 0 0
+DEAL::4 1 0
+DEAL::4 2 16.0000
+DEAL::4 3 0
+DEAL::4 4 32.0000
+DEAL::OK
+DEAL::Chunk size: 7
+DEAL::0 0 0
+DEAL::0 1 0
+DEAL::0 2 0
+DEAL::0 3 0
+DEAL::0 4 0
+DEAL::1 0 0
+DEAL::1 1 2.00000
+DEAL::1 2 4.00000
+DEAL::1 3 0
+DEAL::1 4 0
+DEAL::2 0 0
+DEAL::2 1 0
+DEAL::2 2 8.00000
+DEAL::2 3 12.0000
+DEAL::2 4 0
+DEAL::3 0 0
+DEAL::3 1 6.00000
+DEAL::3 2 0
+DEAL::3 3 18.0000
+DEAL::3 4 24.0000
+DEAL::4 0 0
+DEAL::4 1 0
+DEAL::4 2 16.0000
+DEAL::4 3 0
+DEAL::4 4 32.0000
+DEAL::OK
diff --git a/tests/bits/chunk_sparse_matrix_iterator_08.cc b/tests/bits/chunk_sparse_matrix_iterator_08.cc
new file mode 100644 (file)
index 0000000..5c59cd3
--- /dev/null
@@ -0,0 +1,110 @@
+//----------------------------  chunk_sparse_matrix_iterator_08.cc  ---------------------------
+//    $Id$
+//    Version: $Name$ 
+//
+//    Copyright (C) 2004, 2005, 2013 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  chunk_sparse_matrix_iterator_08.cc  ---------------------------
+
+// test setting some elements using a non-const chunk matrix iterator and
+// operator/=, and reading them back through the matrix itself
+
+#include "../tests.h"
+#include <deal.II/lac/chunk_sparse_matrix.h>    
+#include <fstream>
+#include <iomanip>
+
+
+void test (const unsigned int chunk_size)
+{
+  deallog << "Chunk size: " << chunk_size << std::endl;
+  ChunkSparsityPattern sp (5,5,3,chunk_size);
+  for (unsigned int i=0; i<5; ++i)
+    for (unsigned int j=0; j<5; ++j)
+      if (((i+2*j+1) % 3 == 0)
+          ||
+          (i==j))
+        sp.add (i,j);
+  sp.compress ();
+
+  ChunkSparseMatrix<double> m(sp);
+  for (unsigned int i=0; i<5; ++i)
+    for (unsigned int j=0; j<5; ++j)
+      if (((i+2*j+1) % 3 == 0)
+          ||
+          (i==j))
+        m.set(i,j,i*j);
+
+  ChunkSparseMatrix<double>::iterator i = m.begin();
+  for (; i!=m.end(); ++i)
+    i->value() /= 2;
+  
+  for (unsigned int i=0; i<5; ++i)
+    for (unsigned int j=0; j<5; ++j)
+      if (((i+2*j+1) % 3 == 0)
+          ||
+          (i==j))
+        {
+          deallog << i << ' ' << j << ' ' << m.el(i,j)
+                  << std::endl;
+          Assert (std::fabs(m.el(i,j)-(i*j/2.)) < 1e-14,
+                  ExcInternalError());
+        }
+      else if (sp.exists(i,j))
+        {
+          deallog << i << ' ' << j << ' ' << m.el(i,j)
+                  << std::endl;
+          Assert (std::fabs(m.el(i,j)) < 1e-14,
+                  ExcInternalError());
+        }
+  
+  deallog << "OK" << std::endl;
+}
+
+
+
+int main ()
+{
+  std::ofstream logfile("chunk_sparse_matrix_iterator_08/output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  try
+    {
+      test (1);
+      test (2);
+      test (4);
+      test (5);
+      test (7);
+    }
+  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_iterator_08/cmp/generic b/tests/bits/chunk_sparse_matrix_iterator_08/cmp/generic
new file mode 100644 (file)
index 0000000..38a9f8a
--- /dev/null
@@ -0,0 +1,122 @@
+
+DEAL::Chunk size: 1
+DEAL::0 0 0
+DEAL::0 1 0
+DEAL::0 4 0
+DEAL::1 1 0.500000
+DEAL::1 2 1.00000
+DEAL::2 0 0
+DEAL::2 2 2.00000
+DEAL::2 3 3.00000
+DEAL::3 1 1.50000
+DEAL::3 3 4.50000
+DEAL::3 4 6.00000
+DEAL::4 2 4.00000
+DEAL::4 4 8.00000
+DEAL::OK
+DEAL::Chunk size: 2
+DEAL::0 0 0
+DEAL::0 1 0
+DEAL::0 2 0
+DEAL::0 3 0
+DEAL::0 4 0
+DEAL::1 0 0
+DEAL::1 1 0.500000
+DEAL::1 2 1.00000
+DEAL::1 3 0
+DEAL::1 4 0
+DEAL::2 0 0
+DEAL::2 1 0
+DEAL::2 2 2.00000
+DEAL::2 3 3.00000
+DEAL::2 4 0
+DEAL::3 0 0
+DEAL::3 1 1.50000
+DEAL::3 2 0
+DEAL::3 3 4.50000
+DEAL::3 4 6.00000
+DEAL::4 2 4.00000
+DEAL::4 3 0
+DEAL::4 4 8.00000
+DEAL::OK
+DEAL::Chunk size: 4
+DEAL::0 0 0
+DEAL::0 1 0
+DEAL::0 2 0
+DEAL::0 3 0
+DEAL::0 4 0
+DEAL::1 0 0
+DEAL::1 1 0.500000
+DEAL::1 2 1.00000
+DEAL::1 3 0
+DEAL::1 4 0
+DEAL::2 0 0
+DEAL::2 1 0
+DEAL::2 2 2.00000
+DEAL::2 3 3.00000
+DEAL::2 4 0
+DEAL::3 0 0
+DEAL::3 1 1.50000
+DEAL::3 2 0
+DEAL::3 3 4.50000
+DEAL::3 4 6.00000
+DEAL::4 0 0
+DEAL::4 1 0
+DEAL::4 2 4.00000
+DEAL::4 3 0
+DEAL::4 4 8.00000
+DEAL::OK
+DEAL::Chunk size: 5
+DEAL::0 0 0
+DEAL::0 1 0
+DEAL::0 2 0
+DEAL::0 3 0
+DEAL::0 4 0
+DEAL::1 0 0
+DEAL::1 1 0.500000
+DEAL::1 2 1.00000
+DEAL::1 3 0
+DEAL::1 4 0
+DEAL::2 0 0
+DEAL::2 1 0
+DEAL::2 2 2.00000
+DEAL::2 3 3.00000
+DEAL::2 4 0
+DEAL::3 0 0
+DEAL::3 1 1.50000
+DEAL::3 2 0
+DEAL::3 3 4.50000
+DEAL::3 4 6.00000
+DEAL::4 0 0
+DEAL::4 1 0
+DEAL::4 2 4.00000
+DEAL::4 3 0
+DEAL::4 4 8.00000
+DEAL::OK
+DEAL::Chunk size: 7
+DEAL::0 0 0
+DEAL::0 1 0
+DEAL::0 2 0
+DEAL::0 3 0
+DEAL::0 4 0
+DEAL::1 0 0
+DEAL::1 1 0.500000
+DEAL::1 2 1.00000
+DEAL::1 3 0
+DEAL::1 4 0
+DEAL::2 0 0
+DEAL::2 1 0
+DEAL::2 2 2.00000
+DEAL::2 3 3.00000
+DEAL::2 4 0
+DEAL::3 0 0
+DEAL::3 1 1.50000
+DEAL::3 2 0
+DEAL::3 3 4.50000
+DEAL::3 4 6.00000
+DEAL::4 0 0
+DEAL::4 1 0
+DEAL::4 2 4.00000
+DEAL::4 3 0
+DEAL::4 4 8.00000
+DEAL::OK
diff --git a/tests/bits/chunk_sparse_matrix_iterator_09.cc b/tests/bits/chunk_sparse_matrix_iterator_09.cc
new file mode 100644 (file)
index 0000000..6edddd9
--- /dev/null
@@ -0,0 +1,83 @@
+//----------------------------  chunk_sparse_matrix_iterator_09.cc  ---------------------------
+//    $Id$
+//    Version: $Name$ 
+//
+//    Copyright (C) 2004, 2005, 2013 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  chunk_sparse_matrix_iterator_09.cc  ---------------------------
+
+
+// this test is sparse_matrix_iterator_09 for a ChunkSparseMatrix
+
+#include "../tests.h"
+#include <deal.II/lac/chunk_sparse_matrix.h>
+#include <fstream>
+#include <iomanip>
+
+
+void test (const unsigned int chunk_size)
+{
+  deallog << "Chunk size: " << chunk_size << std::endl;
+
+                                   // create a sparsity pattern with totally
+                                   // empty lines (not even diagonals, since
+                                   // not quadratic)
+  ChunkSparsityPattern sparsity(4,5,1,chunk_size);
+  sparsity.add (1,1);
+  sparsity.add (3,1);
+  sparsity.compress ();
+
+                                   // attach a sparse matrix to it
+  ChunkSparseMatrix<double> A(sparsity);
+
+                                   // and loop over the elements of it
+  for (ChunkSparseMatrix<double>::const_iterator k=A.begin();
+       k!=A.end(); ++k)
+    deallog << k->row() << ' ' << k->column() << ' ' << k->value()
+            << std::endl;
+}
+
+
+
+int main ()
+{
+  std::ofstream logfile("chunk_sparse_matrix_iterator_09/output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  try
+    {
+      test (1);
+      test (3);
+    }
+  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_iterator_09/cmp/generic b/tests/bits/chunk_sparse_matrix_iterator_09/cmp/generic
new file mode 100644 (file)
index 0000000..3d7d6ec
--- /dev/null
@@ -0,0 +1,19 @@
+
+DEAL::Chunk size: 1
+DEAL::1 1 0
+DEAL::3 1 0
+DEAL::Chunk size: 3
+DEAL::0 0 0
+DEAL::0 1 0
+DEAL::0 2 0
+DEAL::1 0 0
+DEAL::1 1 0
+DEAL::1 2 0
+DEAL::2 0 0
+DEAL::2 1 0
+DEAL::2 2 0
+DEAL::3 3 0
+DEAL::3 4 0
+DEAL::3 0 0
+DEAL::3 1 0
+DEAL::3 2 0
diff --git a/tests/bits/chunk_sparse_matrix_iterator_10.cc b/tests/bits/chunk_sparse_matrix_iterator_10.cc
new file mode 100644 (file)
index 0000000..8f788d1
--- /dev/null
@@ -0,0 +1,85 @@
+//----------------------------  chunk_sparse_matrix_iterator_10.cc  ---------------------------
+//    $Id$
+//    Version: $Name$ 
+//
+//    Copyright (C) 2004, 2005, 2013 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  chunk_sparse_matrix_iterator_10.cc  ---------------------------
+
+
+// this test is sparse_matrix_iterator_10 for a ChunkSparseMatrix and the same
+// test as chunk_sparse_matrix_iterator_09 with postfix operator++ instead of
+// prefix
+
+#include "../tests.h"
+#include <deal.II/lac/chunk_sparse_matrix.h>
+#include <fstream>
+#include <iomanip>
+
+
+void test (const unsigned int chunk_size)
+{
+  deallog << "Chunk size: " << chunk_size << std::endl;
+
+                                   // create a sparsity pattern with totally
+                                   // empty lines (not even diagonals, since
+                                   // not quadratic)
+  ChunkSparsityPattern sparsity(4,5,1,chunk_size);
+  sparsity.add (1,1);
+  sparsity.add (3,1);
+  sparsity.compress ();
+
+                                   // attach a sparse matrix to it
+  ChunkSparseMatrix<double> A(sparsity);
+
+                                   // and loop over the elements of it
+  for (ChunkSparseMatrix<double>::const_iterator k=A.begin();
+       k!=A.end(); ++k)
+    deallog << k->row() << ' ' << k->column() << ' ' << k->value()
+            << std::endl;
+}
+
+
+
+int main ()
+{
+  std::ofstream logfile("chunk_sparse_matrix_iterator_10/output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  try
+    {
+      test (1);
+      test (3);
+    }
+  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_iterator_10/cmp/generic b/tests/bits/chunk_sparse_matrix_iterator_10/cmp/generic
new file mode 100644 (file)
index 0000000..3d7d6ec
--- /dev/null
@@ -0,0 +1,19 @@
+
+DEAL::Chunk size: 1
+DEAL::1 1 0
+DEAL::3 1 0
+DEAL::Chunk size: 3
+DEAL::0 0 0
+DEAL::0 1 0
+DEAL::0 2 0
+DEAL::1 0 0
+DEAL::1 1 0
+DEAL::1 2 0
+DEAL::2 0 0
+DEAL::2 1 0
+DEAL::2 2 0
+DEAL::3 3 0
+DEAL::3 4 0
+DEAL::3 0 0
+DEAL::3 1 0
+DEAL::3 2 0
diff --git a/tests/bits/chunk_sparse_matrix_iterator_11.cc b/tests/bits/chunk_sparse_matrix_iterator_11.cc
new file mode 100644 (file)
index 0000000..d5d6560
--- /dev/null
@@ -0,0 +1,95 @@
+//----------------------------  chunk_sparse_matrix_iterator_11.cc  ---------------------------
+//    $Id$
+//    Version: $Name$ 
+//
+//    Copyright (C) 2004, 2005, 2013 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  chunk_sparse_matrix_iterator_11.cc  ---------------------------
+
+
+// comparisons between chunk sparse matrix iterators, same as
+// sparse_matrix_iterators_11 otherwise.
+
+#include "../tests.h"
+#include <deal.II/lac/chunk_sparse_matrix.h>
+#include <fstream>
+#include <iomanip>
+
+
+void test (const unsigned int chunk_size)
+{
+  deallog << "Chunk size: " << chunk_size << std::endl;
+
+                                   // create a sparsity pattern with totally
+                                   // empty lines (not even diagonals, since
+                                   // not quadratic)
+  ChunkSparsityPattern sparsity(4,5,1,chunk_size);
+  sparsity.add (1,1);
+  sparsity.add (3,1);
+  sparsity.compress ();
+
+                                   // attach a sparse matrix to it
+  ChunkSparseMatrix<double> A(sparsity);
+
+  ChunkSparseMatrix<double>::iterator k = A.begin(),
+                                      j = ++A.begin();
+
+  Assert (k < j, ExcInternalError());
+  Assert (j > k, ExcInternalError());
+
+  Assert (!(j < k), ExcInternalError());
+  Assert (!(k > j), ExcInternalError());
+
+  Assert (k != j, ExcInternalError());
+  Assert (!(k == j), ExcInternalError());
+
+  Assert (k == k, ExcInternalError());
+  Assert (!(k != k), ExcInternalError());
+  
+  deallog << "OK" << std::endl;
+}
+
+
+
+int main ()
+{
+  std::ofstream logfile("chunk_sparse_matrix_iterator_11/output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  try
+    {
+      test (1);
+      test (3);
+    }
+  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_iterator_11/cmp/generic b/tests/bits/chunk_sparse_matrix_iterator_11/cmp/generic
new file mode 100644 (file)
index 0000000..764c6e2
--- /dev/null
@@ -0,0 +1,5 @@
+
+DEAL::Chunk size: 1
+DEAL::OK
+DEAL::Chunk size: 3
+DEAL::OK
diff --git a/tests/bits/chunk_sparse_matrix_iterator_12.cc b/tests/bits/chunk_sparse_matrix_iterator_12.cc
new file mode 100644 (file)
index 0000000..63032cd
--- /dev/null
@@ -0,0 +1,94 @@
+//----------------------------  chunk_sparse_matrix_iterator_12.cc  ---------------------------
+//    $Id$
+//    Version: $Name$ 
+//
+//    Copyright (C) 2004, 2005, 2013 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  chunk_sparse_matrix_iterator_12.cc  ---------------------------
+
+
+// like chunk_sparse_matrix_iterator_11 but for const_iterator
+
+#include "../tests.h"
+#include <deal.II/lac/chunk_sparse_matrix.h>
+#include <fstream>
+#include <iomanip>
+
+
+void test (const unsigned int chunk_size)
+{
+  deallog << "Chunk size: " << chunk_size << std::endl;
+
+                                   // create a sparsity pattern with totally
+                                   // empty lines (not even diagonals, since
+                                   // not quadratic)
+  ChunkSparsityPattern sparsity(4,5,1,chunk_size);
+  sparsity.add (1,1);
+  sparsity.add (3,1);
+  sparsity.compress ();
+
+                                   // attach a sparse matrix to it
+  ChunkSparseMatrix<double> A(sparsity);
+
+  ChunkSparseMatrix<double>::const_iterator k = A.begin(),
+                                            j = ++A.begin();
+
+  Assert (k < j, ExcInternalError());
+  Assert (j > k, ExcInternalError());
+
+  Assert (!(j < k), ExcInternalError());
+  Assert (!(k > j), ExcInternalError());
+
+  Assert (k != j, ExcInternalError());
+  Assert (!(k == j), ExcInternalError());
+
+  Assert (k == k, ExcInternalError());
+  Assert (!(k != k), ExcInternalError());
+  
+  deallog << "OK" << std::endl;
+}
+
+
+
+int main ()
+{
+  std::ofstream logfile("chunk_sparse_matrix_iterator_12/output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  try
+    {
+      test (1);
+      test (3);
+    }
+  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_iterator_12/cmp/generic b/tests/bits/chunk_sparse_matrix_iterator_12/cmp/generic
new file mode 100644 (file)
index 0000000..764c6e2
--- /dev/null
@@ -0,0 +1,5 @@
+
+DEAL::Chunk size: 1
+DEAL::OK
+DEAL::Chunk size: 3
+DEAL::OK
diff --git a/tests/bits/chunk_sparse_matrix_iterator_13.cc b/tests/bits/chunk_sparse_matrix_iterator_13.cc
new file mode 100644 (file)
index 0000000..8c38ca9
--- /dev/null
@@ -0,0 +1,117 @@
+//----------------------------  chunk_sparse_matrix_iterator_13.cc  ---------------------------
+//    $Id$
+//    Version: $Name$ 
+//
+//    Copyright (C) 2004, 2005, 2013 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  chunk_sparse_matrix_iterator_13.cc  ---------------------------
+
+
+// test ChunkSparseMatrix::iterator::operator-
+
+#include "../tests.h"
+#include <deal.II/lac/chunk_sparse_matrix.h>
+#include <fstream>
+#include <iomanip>
+
+
+void test (const unsigned int chunk_size)
+{
+  deallog << "Chunk size: " << chunk_size << std::endl;
+
+                                   // create a sparsity pattern with totally
+                                   // empty lines (not even diagonals, since
+                                   // not quadratic)
+  ChunkSparsityPattern sp(5,5,3,chunk_size);
+  for (unsigned int i=0; i<5; ++i)
+    for (unsigned int j=0; j<5; ++j)
+      if (((i+2*j+1) % 3 == 0)
+          ||
+          (i==j))
+        sp.add (i,j);
+  sp.compress ();
+
+                                   // attach a sparse matrix to it
+  ChunkSparseMatrix<double> m(sp);
+
+  for (unsigned int row=0; row<sp.n_rows(); ++row)
+    Assert (m.begin(row)-m.begin(row) == 0,
+           ExcInternalError());
+
+  for (unsigned int row=0; row<sp.n_rows(); ++row)
+    Assert (m.end(row)-m.begin(row) == (int)sp.row_length(row),
+           ExcInternalError());
+  for (unsigned int row=0; row<sp.n_rows(); ++row)
+    Assert (m.begin(row)-m.end(row) == -(int)sp.row_length(row),
+           ExcInternalError());
+
+  {
+    unsigned int counter = 0;
+    for (unsigned int row=0; row<sp.n_rows(); ++row)
+      {
+       Assert (m.begin(row)-m.begin(0) == (int)counter,
+               ExcInternalError());
+       Assert (m.begin(0)-m.begin(row) == -(int)counter,
+               ExcInternalError());
+       counter += sp.row_length(row);
+      }
+  }
+
+  Assert (m.begin() - m.begin(0) == 0, ExcInternalError());
+  Assert (m.begin(0) - m.begin() == 0, ExcInternalError());
+  Assert (m.end(sp.n_rows()-1) - m.end() == 0, ExcInternalError());
+  Assert (m.end() - m.end(sp.n_rows()-1) == 0, ExcInternalError());
+  Assert (m.end() - m.begin() == (int)sp.n_nonzero_elements(),
+         ExcInternalError());
+  Assert (m.begin() - m.end() == -(int)sp.n_nonzero_elements(),
+         ExcInternalError());
+
+  deallog << "OK" << std::endl;
+}
+
+
+
+int main ()
+{
+  std::ofstream logfile("chunk_sparse_matrix_iterator_13/output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  try
+    {
+      test (1);
+      test (2);
+      test (4);
+      test (7);
+    }
+  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_iterator_13/cmp/generic b/tests/bits/chunk_sparse_matrix_iterator_13/cmp/generic
new file mode 100644 (file)
index 0000000..b4056d7
--- /dev/null
@@ -0,0 +1,9 @@
+
+DEAL::Chunk size: 1
+DEAL::OK
+DEAL::Chunk size: 2
+DEAL::OK
+DEAL::Chunk size: 4
+DEAL::OK
+DEAL::Chunk size: 7
+DEAL::OK

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.