]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Add the Chunk version of the SparseMatrix class. Pretty much all preconditioners...
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 6 Aug 2008 05:00:56 +0000 (05:00 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 6 Aug 2008 05:00:56 +0000 (05:00 +0000)
git-svn-id: https://svn.dealii.org/trunk@16496 0785d39b-7218-0410-832d-ea1e28bc413d

49 files changed:
deal.II/lac/include/lac/chunk_sparse_matrix.h [new file with mode: 0644]
deal.II/lac/include/lac/chunk_sparse_matrix.templates.h [new file with mode: 0644]
deal.II/lac/source/chunk_sparse_matrix.cc [new file with mode: 0644]
deal.II/lac/source/chunk_sparse_matrix.inst.in [new file with mode: 0644]
tests/bits/Makefile
tests/bits/chunk_sparse_matrix_00.cc [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_00/cmp/generic [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_01.cc [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_01/cmp/generic [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_01a.cc [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_01a/cmp/generic [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_02.cc [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_02/cmp/generic [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_03.cc [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_03/cmp/generic [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_03a.cc [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_03a/cmp/generic [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_03b.cc [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_03b/cmp/generic [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_04.cc [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_04/cmp/generic [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_05.cc [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_05/cmp/generic [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_05a.cc [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_05a/cmp/generic [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_06.cc [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_06/cmp/generic [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_07.cc [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_07/cmp/generic [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_08.cc [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_08/cmp/generic [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_09.cc [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_09/cmp/generic [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_10.cc [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_10/cmp/generic [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_vector_01.cc [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_vector_01/cmp/generic [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_vector_02.cc [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_vector_02/cmp/generic [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_vector_03.cc [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_vector_03/cmp/generic [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_vector_04.cc [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_vector_04/cmp/generic [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_vector_05.cc [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_vector_05/cmp/generic [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_vector_06.cc [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_vector_06/cmp/generic [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_vector_07.cc [new file with mode: 0644]
tests/bits/chunk_sparse_matrix_vector_07/cmp/generic [new file with mode: 0644]

diff --git a/deal.II/lac/include/lac/chunk_sparse_matrix.h b/deal.II/lac/include/lac/chunk_sparse_matrix.h
new file mode 100644 (file)
index 0000000..986c53b
--- /dev/null
@@ -0,0 +1,1467 @@
+//---------------------------------------------------------------------------
+//    $Id$
+//    Version: $Name:  $
+//
+//    Copyright (C) 2008 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.
+//
+//---------------------------------------------------------------------------
+#ifndef __deal2__chunk_sparse_matrix_h
+#define __deal2__chunk_sparse_matrix_h
+
+
+#include <base/config.h>
+#include <base/subscriptor.h>
+#include <base/smartpointer.h>
+#include <lac/chunk_sparsity_pattern.h>
+#include <lac/identity_matrix.h>
+#include <lac/exceptions.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+template<typename number> class Vector;
+template<typename number> class FullMatrix;
+
+/*! @addtogroup Matrix1
+ *@{
+ */
+
+
+/**
+ * Sparse matrix. This class implements the function to store values
+ * in the locations of a sparse matrix denoted by a
+ * SparsityPattern. The separation of sparsity pattern and values is
+ * done since one can store data elements of different type in these
+ * locations without the SparsityPattern having to know this, and more
+ * importantly one can associate more than one matrix with the same
+ * sparsity pattern.
+ *
+ * @note Instantiations for this template are provided for <tt>@<float@> and
+ * @<double@></tt>; others can be generated in application programs (see the
+ * section on @ref Instantiations in the manual).
+ *
+ * @author Wolfgang Bangerth, 2008
+ */
+template <typename number>
+class ChunkSparseMatrix : public virtual Subscriptor
+{
+  public:
+                                    /**
+                                     * Type of matrix entries. In analogy to
+                                     * the STL container classes.
+                                     */
+    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.
+                                     *
+                                     * 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.
+                                      */
+    struct Traits
+    {
+                                         /**
+                                          * It is safe to elide additions of
+                                          * zeros to individual elements of
+                                          * this matrix.
+                                          */
+        static const bool zero_addition_can_be_elided = true;
+    };
+    
+/**
+ * @name Constructors and initalization.
+ */
+//@{
+                                    /**
+                                     * 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&).
+                                     */
+    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.
+                                     *
+                                     * 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.
+                                     *
+                                     * 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.
+                                     */
+    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.
+                                     */
+    ChunkSparseMatrix (const ChunkSparsityPattern &sparsity,
+                      const IdentityMatrix  &id);
+    
+                                    /**
+                                     * 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.
+                                     */
+    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.
+                                     */
+    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
+                                      * 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
+                                     * reserved.
+                                     *
+                                     * 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.
+                                     *
+                                     * 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
+                                     * previously tied to.
+                                     */
+    virtual void clear ();
+//@}    
+/**
+ * @name Information on the matrix
+ */
+//@{
+                                    /**
+                                     * 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$.
+                                     */
+    unsigned int m () const;
+    
+                                    /**
+                                     * 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.
+                                     */
+    unsigned int n_nonzero_elements () const;
+
+                                    /**
+                                     * 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.
+                                     */
+    unsigned int n_actually_nonzero_elements () const;
+    
+                                    /**
+                                     * 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.
+                                     */
+    const ChunkSparsityPattern & get_sparsity_pattern () const;
+
+                                    /**
+                                     * Determine an estimate for the
+                                     * memory consumption (in bytes)
+                                     * of this object. See
+                                     * MemoryConsumption.
+                                     */
+    unsigned int memory_consumption () const;
+    
+//@}
+/**
+ * @name Modifying entries
+ */
+//@{
+                                    /**
+                                     * 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.
+                                     */
+    void add (const unsigned int i,
+              const unsigned int j,
+             const number value);
+
+                                    /**
+                                     * Multiply the entire matrix by a
+                                     * fixed factor.
+                                     */
+    ChunkSparseMatrix & operator *= (const number 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)$.
+                                     *
+                                     * 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.
+                                     *
+                                     * 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>.
+                                     */
+    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.
+                                     *
+                                     * 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.
+                                     */
+    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.
+                                     */
+    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.
+                                     *
+                                     * 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>
+    void add (const number factor,
+             const ChunkSparseMatrix<somenumber> &matrix);    
+    
+//@}
+/**
+ * @name Entry Access
+ */
+//@{
+    
+                                    /**
+                                     * 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.
+                                     *
+                                     * 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.
+                                     *
+                                     * If you are looping over all elements,
+                                     * consider using one of the iterator
+                                     * classes instead, since they are
+                                     * tailored better to a sparse matrix
+                                     * structure.
+                                     */
+    number el (const unsigned int i,
+              const unsigned int j) const;
+
+                                    /**
+                                     * Return the main diagonal
+                                     * element in the <i>i</i>th
+                                     * row. This function throws an
+                                     * error if the matrix is not
+                                     * quadratic (see
+                                     * ChunkSparsityPattern::optimize_diagonal()).
+                                     *
+                                     * 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?
+                                     */
+    number & diag_element (const unsigned int i);
+
+//@}
+/**
+ * @name Matrix vector multiplications
+ */
+//@{
+                                    /**
+                                     * 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.
+                                     * 
+                                      * 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.
+                                      *
+                                     * 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.
+                                     */
+    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.
+                                      *
+                                     * 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.
+                                     */
+    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.
+                                      *
+                                     * 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.
+                                     */
+    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.
+                                     *
+                                     * 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.
+                                     */
+    template <typename somenumber>
+    somenumber matrix_norm_square (const Vector<somenumber> &v) const;
+
+                                    /**
+                                     * 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.
+                                      *
+                                      * Source <i>x</i> and destination
+                                      * <i>dst</i> must not be the same
+                                      * vector.
+                                     */
+    template <typename somenumber>
+    somenumber residual (Vector<somenumber>       &dst,
+                        const Vector<somenumber> &x,
+                        const Vector<somenumber> &b) const;
+    
+//@}
+/**
+ * @name Matrix norms
+ */
+//@{
+    
+                                    /**
+                                     * 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)
+                                     */
+    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.
+                                      */
+    real_type frobenius_norm () const;
+//@}
+/**
+ * @name Preconditioning methods
+ */
+//@{
+    
+                                    /**
+                                     * 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,
+                             const Vector<somenumber> &src,
+                             const number              omega = 1.) const;
+
+                                    /**
+                                     * Apply SSOR preconditioning to
+                                     * <tt>src</tt>.
+                                     */
+    template <typename somenumber>
+    void precondition_SSOR (Vector<somenumber>       &dst,
+                           const Vector<somenumber> &src,
+                           const number              om = 1.) const;
+
+                                    /**
+                                     * Apply SOR preconditioning
+                                     * matrix to <tt>src</tt>.
+                                     */
+    template <typename somenumber>
+    void precondition_SOR (Vector<somenumber>       &dst,
+                          const Vector<somenumber> &src,
+                          const number              om = 1.) const;
+    
+                                    /**
+                                     * Apply transpose SOR
+                                     * preconditioning matrix to
+                                     * <tt>src</tt>.
+                                     */
+    template <typename somenumber>
+    void precondition_TSOR (Vector<somenumber>       &dst,
+                           const Vector<somenumber> &src,
+                           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
+                                     * parameter.
+                                     */
+    template <typename somenumber>
+    void SSOR (Vector<somenumber> &v,
+              const number        omega = 1.) const;
+
+                                    /**
+                                     * 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
+                                     * relaxation parameter.
+                                     */
+    template <typename somenumber>
+    void TSOR (Vector<somenumber> &v,
+              const number        om = 1.) const;
+
+                                    /**
+                                     * 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.
+                                     *
+                                     * <tt>omega</tt> is the
+                                     * relaxation parameter.
+                                     */
+    template <typename somenumber>
+    void PSOR (Vector<somenumber> &v,
+              const std::vector<unsigned int>& permutation,
+              const std::vector<unsigned int>& inverse_permutation,
+              const number        om = 1.) const;
+
+                                    /**
+                                     * 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.
+                                     *
+                                     * <tt>omega</tt> is the
+                                     * relaxation parameter.
+                                     */
+    template <typename somenumber>
+    void TPSOR (Vector<somenumber> &v,
+               const std::vector<unsigned int>& permutation,
+               const std::vector<unsigned int>& inverse_permutation,
+               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>.
+                                     */
+    template <typename somenumber>
+    void SOR_step (Vector<somenumber> &v,
+                  const Vector<somenumber> &b,
+                  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>.
+                                     */
+    template <typename somenumber>
+    void TSOR_step (Vector<somenumber> &v,
+                   const Vector<somenumber> &b,
+                   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.
+                                     */
+    template <typename somenumber>
+    void SSOR_step (Vector<somenumber> &v,
+                   const Vector<somenumber> &b,
+                   const number        om = 1.) const;
+//@}
+/**
+ * @name Input/Output
+ */
+//@{
+
+                                    /**
+                                     * 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.
+                                     *
+                                     * 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.
+                                     *
+                                     * 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!
+                                     */
+    void print_formatted (std::ostream       &out,
+                         const unsigned int  precision   = 3,
+                         const bool          scientific  = true,
+                         const unsigned int  width       = 0,
+                         const char         *zero_string = " ",
+                         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.
+                                     */
+    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.
+                                     *
+                                     * 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. 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.
+                                     */
+    void block_read (std::istream &in);
+//@}
+                                    /** @addtogroup Exceptions
+                                     * @{ */
+
+                                    /**
+                                     * Exception
+                                     */
+    DeclException2 (ExcInvalidIndex,
+                   int, int,
+                   << "The entry with index <" << arg1 << ',' << arg2
+                   << "> does not exist.");
+                                    /**
+                                     * Exception
+                                     */
+    DeclException1 (ExcInvalidIndex1,
+                   int,
+                   << "The index " << arg1 << " is not in the allowed range.");
+                                    /**
+                                     * Exception
+                                     */
+    DeclException0 (ExcDifferentChunkSparsityPatterns);
+                                    /**
+                                     * Exception
+                                     */
+    DeclException2 (ExcIteratorRange,
+                   int, int,
+                   << "The iterators denote a range of " << arg1
+                   << " elements, but the given number of rows was " << arg2);
+                                     /**
+                                      * Exception
+                                      */
+    DeclException0 (ExcSourceEqualsDestination);
+                                    //@}    
+  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.
+                                     */
+    SmartPointer<const ChunkSparsityPattern> 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.
+                                     */
+    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()
+                                     * function.
+                                     */
+    unsigned int max_len;
+    
+                                    /**
+                                     * 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
+    template <typename somenumber> friend class ChunkSparseMatrix;
+};
+
+/*@}*/
+
+#ifndef DOXYGEN
+/*---------------------- Inline functions -----------------------------------*/
+
+
+
+template <typename number>
+inline
+unsigned int ChunkSparseMatrix<number>::m () const
+{
+  Assert (cols != 0, ExcNotInitialized());
+  return cols->rows;
+}
+
+
+template <typename number>
+inline
+unsigned int ChunkSparseMatrix<number>::n () const
+{
+  Assert (cols != 0, ExcNotInitialized());
+  return cols->cols;
+}
+
+
+
+template <typename number>
+inline
+unsigned int
+ChunkSparseMatrix<number>::compute_location (const unsigned int i,
+                                            const unsigned int j) const
+{
+  const unsigned int chunk_size = cols->get_chunk_size();
+  const unsigned int chunk_index
+    = cols->sparsity_pattern(i/chunk_size, j/chunk_size);
+
+  if (chunk_index == ChunkSparsityPattern::invalid_entry)
+    return ChunkSparsityPattern::invalid_entry;
+  else
+    {
+      return (chunk_index * chunk_size * chunk_size
+             +
+             (i % chunk_size) * chunk_size
+             +
+             (j % chunk_size));
+    }
+}
+
+
+template <typename number>
+inline
+void ChunkSparseMatrix<number>::set (const unsigned int i,
+                                    const unsigned int j,
+                                    const number value)
+{
+
+  Assert (numbers::is_finite(value), 
+          ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+
+  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
+  const unsigned int index = compute_location(i,j);
+  Assert ((index != SparsityPattern::invalid_entry) ||
+         (value == 0.),
+         ExcInvalidIndex(i,j));
+
+  if (index != SparsityPattern::invalid_entry)
+    val[index] = value;
+}
+
+
+
+template <typename number>
+inline
+void ChunkSparseMatrix<number>::add (const unsigned int i,
+                                    const unsigned int j,
+                                    const number value)
+{
+
+  Assert (numbers::is_finite(value), 
+          ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)"));
+
+  Assert (cols != 0, ExcNotInitialized());
+
+  const unsigned int index = compute_location(i,j);
+  Assert ((index != ChunkSparsityPattern::invalid_entry) ||
+         (value == 0.),
+         ExcInvalidIndex(i,j));
+
+  if (value != 0.)
+    val[index] += value;
+}
+
+
+
+template <typename number>
+inline
+ChunkSparseMatrix<number> &
+ChunkSparseMatrix<number>::operator *= (const number factor)
+{
+  Assert (cols != 0, ExcNotInitialized());
+  Assert (val != 0, ExcNotInitialized());
+
+  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
+  number             *val_ptr    = val;
+  const number *const end_ptr    = val +
+                                  cols->sparsity_pattern.n_nonzero_elements()
+                                  *
+                                  chunk_size * chunk_size;
+  while (val_ptr != end_ptr)
+    *val_ptr++ *= factor;
+
+  return *this;
+}
+
+
+
+template <typename number>
+inline
+ChunkSparseMatrix<number> &
+ChunkSparseMatrix<number>::operator /= (const number factor)
+{
+  Assert (cols != 0, ExcNotInitialized());
+  Assert (val != 0, ExcNotInitialized());
+  Assert (factor !=0, ExcDivideByZero());
+
+  const number factor_inv = 1. / 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
+  number             *val_ptr    = val;
+  const number *const end_ptr    = val +
+                                  cols->sparsity_pattern.n_nonzero_elements()
+                                  *
+                                  chunk_size * chunk_size;
+
+  while (val_ptr != end_ptr)
+    *val_ptr++ *= factor_inv;
+
+  return *this;
+}
+
+
+
+template <typename number>
+inline
+number ChunkSparseMatrix<number>::operator () (const unsigned int i,
+                                              const unsigned int j) const
+{
+  Assert (cols != 0, ExcNotInitialized());
+  AssertThrow (compute_location(i,j) != SparsityPattern::invalid_entry,
+               ExcInvalidIndex(i,j));
+  return val[compute_location(i,j)];
+}
+
+
+
+template <typename number>
+inline
+number ChunkSparseMatrix<number>::el (const unsigned int i,
+                                     const unsigned int j) const
+{
+  Assert (cols != 0, ExcNotInitialized());
+  const unsigned int index = compute_location(i,j);
+
+  if (index != ChunkSparsityPattern::invalid_entry)
+    return val[index];
+  else
+    return 0;
+}
+
+
+
+template <typename number>
+inline
+number ChunkSparseMatrix<number>::diag_element (const unsigned int i) const
+{
+  Assert (cols != 0, ExcNotInitialized());
+  Assert (cols->optimize_diagonal(),  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 (cols->optimize_diagonal(),  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>
+template <typename ForwardIterator>
+void
+ChunkSparseMatrix<number>::copy_from (const ForwardIterator begin,
+                                     const ForwardIterator end)
+{
+  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
+  typedef typename std::iterator_traits<ForwardIterator>::value_type::const_iterator inner_iterator;
+  unsigned int row=0;
+  for (ForwardIterator i=begin; i!=end; ++i, ++row)
+    {
+      const inner_iterator end_of_row = i->end();
+      for (inner_iterator j=i->begin(); j!=end_of_row; ++j)
+                                        // write entries
+       set (row, j->first, j->second);
+    }
+}
+
+
+#endif // DOXYGEN
+
+
+/*----------------------------   chunk_sparse_matrix.h     ---------------------------*/
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
+/*----------------------------   chunk_sparse_matrix.h     ---------------------------*/
diff --git a/deal.II/lac/include/lac/chunk_sparse_matrix.templates.h b/deal.II/lac/include/lac/chunk_sparse_matrix.templates.h
new file mode 100644 (file)
index 0000000..b0d63b6
--- /dev/null
@@ -0,0 +1,1591 @@
+//---------------------------------------------------------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2008 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.
+//
+//---------------------------------------------------------------------------
+#ifndef __deal2__chunk_sparse_matrix_templates_h
+#define __deal2__chunk_sparse_matrix_templates_h
+
+
+#include <base/template_constraints.h>
+#include <lac/chunk_sparse_matrix.h>
+#include <lac/vector.h>
+#include <lac/full_matrix.h>
+
+
+// we only need output streams, but older compilers did not provide
+// them in a separate include file
+#ifdef HAVE_STD_OSTREAM_HEADER
+#  include <ostream>
+#else
+#  include <iostream>
+#endif
+
+#include <iomanip>
+#include <algorithm>
+#include <functional>
+#include <cmath>
+
+#include <vector>
+#include <numeric>
+
+#include <base/thread_management.h>
+#include <base/multithread_info.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+
+namespace internal
+{
+//TODO: the goal of the ChunkSparseMatrix class is to stream data and use
+// the vectorization features of modern processors. to make this happen,
+// we will have to vectorize the functions in the following namespace, either
+// by hand or by using, for example, optimized BLAS versions for them.  
+  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.
+                                     */
+    template <typename MatrixIterator,
+             typename SrcIterator,
+             typename DstIterator>
+    inline
+    void
+    chunk_vmult_add (const unsigned int   chunk_size,
+                    const MatrixIterator matrix,
+                    const SrcIterator    src,
+                    DstIterator          dst)
+    {
+      MatrixIterator matrix_row = matrix;
+      
+      for (unsigned int i=0; i<chunk_size;
+          ++i,  matrix_row += chunk_size)
+       {
+         typename std::iterator_traits<DstIterator>::value_type
+           sum = 0;
+         
+         for (unsigned int j=0; j<chunk_size; ++j)
+           sum += matrix_row[j] * src[j];
+
+         dst[i] += sum;
+       }
+    }
+
+
+
+                                    /**
+                                     * Like the previous function, but
+                                     * subtract. We need this for computing
+                                     * the residual.
+                                     */
+    template <typename MatrixIterator,
+             typename SrcIterator,
+             typename DstIterator>
+    inline
+    void
+    chunk_vmult_subtract (const unsigned int   chunk_size,
+                         const MatrixIterator matrix,
+                         const SrcIterator    src,
+                         DstIterator          dst)
+    {
+      MatrixIterator matrix_row = matrix;
+      
+      for (unsigned int i=0; i<chunk_size;
+          ++i,  matrix_row += chunk_size)
+       {
+         typename std::iterator_traits<DstIterator>::value_type
+           sum = 0;
+         
+         for (unsigned int j=0; j<chunk_size; ++j)
+           sum += matrix_row[j] * src[j];
+
+         dst[i] -= sum;
+       }
+    }
+
+
+                                    /**
+                                     * 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,
+             typename DstIterator>
+    inline
+    void
+    chunk_Tvmult_add (const unsigned int   chunk_size,
+                     const MatrixIterator matrix,
+                     const SrcIterator    src,
+                     DstIterator          dst)
+    {
+      for (unsigned int i=0; i<chunk_size; ++i)
+       {
+         typename std::iterator_traits<DstIterator>::value_type
+           sum = 0;
+         
+         for (unsigned int j=0; j<chunk_size; ++j)
+           sum += matrix[j*chunk_size+i] * src[j];
+
+         dst[i] += sum;
+       }
+    }
+
+
+                                    /**
+                                     * Produce the result of the matrix
+                                     * scalar product $u^TMv$ for an
+                                     * individual chunk.
+                                     */
+    template <typename result_type,
+             typename MatrixIterator,
+             typename SrcIterator1,
+             typename SrcIterator2>
+    inline
+    result_type
+    chunk_matrix_scalar_product (const unsigned int   chunk_size,
+                                const MatrixIterator matrix,
+                                const SrcIterator1   u,
+                                const SrcIterator2   v)
+    {
+      result_type result = 0;
+      
+      MatrixIterator matrix_row = matrix;
+      
+      for (unsigned int i=0; i<chunk_size;
+          ++i,  matrix_row += chunk_size)
+       {
+         typename std::iterator_traits<SrcIterator2>::value_type
+           sum = 0;
+         
+         for (unsigned int j=0; j<chunk_size; ++j)
+           sum += matrix_row[j] * v[j];
+
+         result += u[i] * sum;
+       }
+
+      return result;
+    }
+}
+}
+
+
+
+template <typename number>
+ChunkSparseMatrix<number>::ChunkSparseMatrix ()
+                :
+               cols(0, "ChunkSparseMatrix"),
+               val(0),
+               max_len(0)
+{}
+
+
+
+template <typename number>
+ChunkSparseMatrix<number>::ChunkSparseMatrix (const ChunkSparseMatrix &m)
+                :
+               Subscriptor (m),
+               cols(0, "ChunkSparseMatrix"),
+               val(0),
+               max_len(0)
+{
+  Assert (m.cols==0, ExcInvalidConstructorCall());
+  Assert (m.val==0, ExcInvalidConstructorCall());
+  Assert (m.max_len==0, ExcInvalidConstructorCall());
+}
+
+
+
+template <typename number>
+ChunkSparseMatrix<number>&
+ChunkSparseMatrix<number>::operator = (const ChunkSparseMatrix<number> &m)
+{
+  Assert (m.cols==0, ExcInvalidConstructorCall());
+  Assert (m.val==0, ExcInvalidConstructorCall());
+  Assert (m.max_len==0, ExcInvalidConstructorCall());
+
+  return *this;
+}
+
+
+
+template <typename number>
+ChunkSparseMatrix<number>::ChunkSparseMatrix (const ChunkSparsityPattern &c)
+                :
+               cols(0, "ChunkSparseMatrix"),
+               val(0),
+               max_len(0)
+{
+  reinit (c);
+}
+
+
+
+template <typename number>
+ChunkSparseMatrix<number>::ChunkSparseMatrix (const ChunkSparsityPattern &c,
+                                             const IdentityMatrix  &id)
+                :
+               cols(0, "ChunkSparseMatrix"),
+               val(0),
+               max_len(0)
+{
+  Assert (c.n_rows() == id.m(), ExcDimensionMismatch (c.n_rows(), id.m()));
+  Assert (c.n_cols() == id.n(), ExcDimensionMismatch (c.n_cols(), id.n()));
+
+  reinit (c);
+  for (unsigned int i=0; i<n(); ++i)
+    this->set(i,i,1.);
+}
+
+
+
+template <typename number>
+ChunkSparseMatrix<number>::~ChunkSparseMatrix ()
+{
+  cols = 0;
+  
+  if (val != 0)
+    delete[] val;
+}
+
+
+
+template <typename number>
+ChunkSparseMatrix<number> &
+ChunkSparseMatrix<number>::operator = (const double d)
+{
+  Assert (d==0, ExcScalarAssignmentOnlyForZeroValue());
+  
+  Assert (cols != 0, ExcNotInitialized());
+  Assert (cols->sparsity_pattern.compressed || cols->empty(),
+         ChunkSparsityPattern::ExcNotCompressed());
+
+  if (val)
+    {
+      const unsigned int chunk_size = cols->get_chunk_size();
+      std::fill_n (val,
+                  val +
+                  cols->sparsity_pattern.n_nonzero_elements() *
+                  chunk_size * chunk_size,
+                  0.);
+    }
+
+  return *this;
+}
+
+
+
+template <typename number>
+ChunkSparseMatrix<number> &
+ChunkSparseMatrix<number>::operator= (const IdentityMatrix  &id)
+{
+  Assert (cols->n_rows() == id.m(),
+         ExcDimensionMismatch (cols->n_rows(), id.m()));
+  Assert (cols->n_cols() == id.n(),
+         ExcDimensionMismatch (cols->n_cols(), id.n()));
+
+  *this = 0;
+  for (unsigned int i=0; i<n(); ++i)
+    this->set(i,i,1.);
+
+  return *this;
+}
+
+
+
+template <typename number>
+void
+ChunkSparseMatrix<number>::reinit (const ChunkSparsityPattern &sparsity)
+{
+  cols = &sparsity;
+
+  if (cols->empty())
+    {
+      if (val != 0)
+        delete[] val;
+      val = 0;
+      max_len = 0;
+      return;
+    }
+
+                                  // 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;
+  if (N > max_len || max_len == 0)
+    {
+      if (val != 0)
+        delete[] val;
+      val = new number[N];
+      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);
+}
+
+
+
+template <typename number>
+void
+ChunkSparseMatrix<number>::clear ()
+{
+  cols = 0;
+  if (val) delete[] val;
+  val = 0;
+  max_len = 0;
+}
+
+
+
+template <typename number>
+bool
+ChunkSparseMatrix<number>::empty () const
+{
+  if (cols == 0)
+    return true;
+  else
+    return cols->empty();
+}
+
+
+
+template <typename number>
+unsigned int
+ChunkSparseMatrix<number>::n_nonzero_elements () const
+{
+  Assert (cols != 0, ExcNotInitialized());
+  return cols->n_nonzero_elements ();
+}
+
+
+
+template <typename number>
+unsigned int
+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
+  const unsigned int chunk_size = cols->get_chunk_size();
+  return std::count_if(&val[0],
+                      &val[cols->sparsity_pattern.n_nonzero_elements () *
+                           chunk_size * chunk_size],
+                      std::bind2nd(std::not_equal_to<double>(), 0));
+}
+
+
+
+template <typename number>
+void
+ChunkSparseMatrix<number>::symmetrize ()
+{
+  Assert (cols != 0, ExcNotInitialized());
+  Assert (cols->rows == cols->cols, ExcNotQuadratic());
+
+  Assert (false, ExcNotImplemented());
+}
+
+
+
+template <typename number>
+template <typename somenumber>
+ChunkSparseMatrix<number> &
+ChunkSparseMatrix<number>::copy_from (const ChunkSparseMatrix<somenumber> &matrix)
+{
+  Assert (cols != 0, ExcNotInitialized());
+  Assert (val != 0, ExcNotInitialized());
+  Assert (cols == matrix.cols, ExcDifferentChunkSparsityPatterns());
+
+                                  // 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()
+                        * chunk_size * chunk_size],
+            &val[0]);
+  
+  return *this;
+}
+
+
+
+template <typename number>
+template <typename somenumber>
+void
+ChunkSparseMatrix<number>::copy_from (const FullMatrix<somenumber> &matrix)
+{
+                                  // first delete previous content
+  *this = 0;
+
+                                  // then copy old matrix
+  for (unsigned int row=0; row<matrix.m(); ++row)
+    for (unsigned int col=0; col<matrix.n(); ++col)
+      if (matrix(row,col) != 0)
+       set (row, col, matrix(row,col));
+}
+
+
+
+template <typename number>
+template <typename somenumber>
+void
+ChunkSparseMatrix<number>::add (const number factor,
+                               const ChunkSparseMatrix<somenumber> &matrix)
+{
+  Assert (cols != 0, ExcNotInitialized());
+  Assert (val != 0, ExcNotInitialized());
+  Assert (cols == matrix.cols, ExcDifferentChunkSparsityPatterns());
+
+                                  // 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];
+  const number *const end_ptr    = &val[cols->sparsity_pattern.n_nonzero_elements()
+                                       * chunk_size * chunk_size];
+
+  while (val_ptr != end_ptr)
+    *val_ptr++ += factor * *matrix_ptr++;
+}
+
+
+template <typename number>
+template <class OutVector, class InVector>
+void
+ChunkSparseMatrix<number>::vmult (OutVector& dst,
+                                 const InVector& src) const
+{
+  Assert (cols != 0, ExcNotInitialized());
+  Assert (val != 0, ExcNotInitialized());
+  Assert(m() == dst.size(), ExcDimensionMismatch(m(),dst.size()));
+  Assert(n() == src.size(), ExcDimensionMismatch(n(),src.size()));
+
+  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
+  dst = 0;
+  vmult_add (dst, src);
+}
+
+
+
+template <typename number>
+template <class OutVector, class InVector>
+void
+ChunkSparseMatrix<number>::Tvmult (OutVector& dst,
+                                  const InVector& src) const
+{
+  Assert (val != 0, ExcNotInitialized());
+  Assert (cols != 0, ExcNotInitialized());
+  Assert(n() == dst.size(), ExcDimensionMismatch(n(),dst.size()));
+  Assert(m() == src.size(), ExcDimensionMismatch(m(),src.size()));
+
+  Assert (!PointerComparison::equal(&src, &dst), ExcSourceEqualsDestination());
+
+  Assert (cols != 0, ExcNotInitialized());
+  Assert (val != 0, ExcNotInitialized());
+  Assert(m() == dst.size(), ExcDimensionMismatch(m(),dst.size()));
+  Assert(n() == src.size(), ExcDimensionMismatch(n(),src.size()));
+
+  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
+  dst = 0;
+  Tvmult_add (dst, src);
+}
+
+
+
+template <typename number>
+template <class OutVector, class InVector>
+void
+ChunkSparseMatrix<number>::vmult_add (OutVector& dst,
+                                     const InVector& src) const
+{
+  Assert (cols != 0, ExcNotInitialized());
+  Assert (val != 0, ExcNotInitialized());
+  Assert(m() == dst.size(), ExcDimensionMismatch(m(),dst.size()));
+  Assert(n() == src.size(), ExcDimensionMismatch(n(),src.size()));
+
+  Assert (!PointerComparison::equal(&src, &dst), ExcSourceEqualsDestination());
+  
+  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;
+       }
+    }
+}
+
+
+template <typename number>
+template <class OutVector, class InVector>
+void
+ChunkSparseMatrix<number>::Tvmult_add (OutVector& dst,
+                                      const InVector& src) const
+{
+  Assert (cols != 0, ExcNotInitialized());
+  Assert (val != 0, ExcNotInitialized());
+  Assert(m() == dst.size(), ExcDimensionMismatch(m(),dst.size()));
+  Assert(n() == src.size(), ExcDimensionMismatch(n(),src.size()));
+
+  Assert (!PointerComparison::equal(&src, &dst), ExcSourceEqualsDestination());
+  
+  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);
+
+                                  // 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;
+          
+  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_Tvmult_add
+             (cols->chunk_size,
+              val_ptr,
+              src.begin() + chunk_row * cols->chunk_size,
+              dst.begin() + *colnum_ptr * cols->chunk_size);
+         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(*colnum_ptr * cols->chunk_size + c)
+                 += (val_ptr[r*cols->chunk_size + c] *
+                     src(chunk_row * cols->chunk_size + r));
+
+         ++colnum_ptr;
+         val_ptr += cols->chunk_size * 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(*colnum_ptr * cols->chunk_size + c)
+                   += (val_ptr[r*cols->chunk_size + c] *
+                       src(chunk_row * cols->chunk_size + r));
+           }
+         else
+                                            // we're at a chunk row and
+                                            // column that has padding
+           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)
+                 += (val_ptr[r*cols->chunk_size + c] *
+                     src(chunk_row * cols->chunk_size + r));
+
+         ++colnum_ptr;
+         val_ptr += cols->chunk_size * cols->chunk_size;
+       }
+    }
+}
+
+
+template <typename number>
+template <typename somenumber>
+somenumber
+ChunkSparseMatrix<number>::matrix_norm_square (const Vector<somenumber>& v) const
+{
+  Assert (cols != 0, ExcNotInitialized());
+  Assert (val != 0, ExcNotInitialized());
+  Assert(m() == v.size(), ExcDimensionMismatch(m(),v.size()));
+  Assert(n() == v.size(), ExcDimensionMismatch(n(),v.size()));
+
+  somenumber result = 0;
+
+                                  ////////////////
+                                  // 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
+  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 Vector<somenumber>::const_iterator v_ptr = v.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))
+           result +=
+             internal::ChunkSparseMatrix::
+             chunk_matrix_scalar_product<somenumber>
+             (cols->chunk_size,
+              val_ptr,
+              v_ptr,
+              v.begin() + *colnum_ptr * cols->chunk_size);
+         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)
+               result
+                 +=
+                 v(chunk_row * cols->chunk_size + r)
+                 * (val_ptr[r*cols->chunk_size + c] *
+                    v(*colnum_ptr * cols->chunk_size + c));
+
+         ++colnum_ptr;
+         val_ptr += cols->chunk_size * cols->chunk_size;
+       }
+      
+
+      v_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)
+                 result
+                   +=
+                   v(chunk_row * cols->chunk_size + r)
+                   * (val_ptr[r*cols->chunk_size + c] *
+                      v(*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)
+               result
+                 +=
+                 v(chunk_row * cols->chunk_size + r)
+                 * (val_ptr[r*cols->chunk_size + c] *
+                    v(*colnum_ptr * cols->chunk_size + c));
+
+         ++colnum_ptr;
+         val_ptr += cols->chunk_size * cols->chunk_size;
+       }
+    }
+
+  return result;
+}
+
+
+
+template <typename number>
+template <typename somenumber>
+somenumber
+ChunkSparseMatrix<number>::matrix_scalar_product (const Vector<somenumber>& u,
+                                                 const Vector<somenumber>& v) const
+{
+  Assert (cols != 0, ExcNotInitialized());
+  Assert (val != 0, ExcNotInitialized());
+  Assert(m() == u.size(), ExcDimensionMismatch(m(),u.size()));
+  Assert(n() == v.size(), ExcDimensionMismatch(n(),v.size()));
+
+                                  // 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
+  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 Vector<somenumber>::const_iterator u_ptr = u.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))
+           result +=
+             internal::ChunkSparseMatrix::
+             chunk_matrix_scalar_product<somenumber>
+             (cols->chunk_size,
+              val_ptr,
+              u_ptr,
+              v.begin() + *colnum_ptr * cols->chunk_size);
+         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)
+               result
+                 +=
+                 u(chunk_row * cols->chunk_size + r)
+                 * (val_ptr[r*cols->chunk_size + c] *
+                    v(*colnum_ptr * cols->chunk_size + c));
+
+         ++colnum_ptr;
+         val_ptr += cols->chunk_size * cols->chunk_size;
+       }
+      
+
+      u_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)
+                 result
+                   +=
+                   u(chunk_row * cols->chunk_size + r)
+                   * (val_ptr[r*cols->chunk_size + c] *
+                      v(*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)
+               result
+                 +=
+                 u(chunk_row * cols->chunk_size + r)
+                 * (val_ptr[r*cols->chunk_size + c] *
+                    v(*colnum_ptr * cols->chunk_size + c));
+
+         ++colnum_ptr;
+         val_ptr += cols->chunk_size * cols->chunk_size;
+       }
+    }
+
+  return result;
+}
+
+
+
+template <typename number>
+typename ChunkSparseMatrix<number>::real_type
+ChunkSparseMatrix<number>::l1_norm () const
+{
+  Assert (cols != 0, ExcNotInitialized());
+  Assert (val != 0, ExcNotInitialized());
+
+  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
+  Vector<real_type> column_sums(cols->sparsity_pattern.n_cols() *
+                               cols->chunk_size);
+
+  for (unsigned int chunk_row=0; chunk_row<n_chunk_rows; ++chunk_row)
+    for (unsigned int j=cols->sparsity_pattern.rowstart[chunk_row];
+        j<cols->sparsity_pattern.rowstart[chunk_row+1] ; ++j)
+      for (unsigned int r=0; r<cols->chunk_size; ++r)
+       for (unsigned int s=0; s<cols->chunk_size; ++s)
+         column_sums(cols->sparsity_pattern.colnums[j] *
+                     cols->chunk_size + s) +=
+           numbers::NumberTraits<number>::abs(val[j * cols->chunk_size *
+                                                  cols->chunk_size +
+                                                  r * cols->chunk_size +
+                                                  s]);
+
+  return column_sums.linfty_norm();
+}
+
+
+
+template <typename number>
+typename ChunkSparseMatrix<number>::real_type
+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
+  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
+  Vector<real_type> row_sums(cols->sparsity_pattern.n_rows() *
+                            cols->chunk_size);
+
+  for (unsigned int chunk_row=0; chunk_row<n_chunk_rows; ++chunk_row)
+    for (unsigned int j=cols->sparsity_pattern.rowstart[chunk_row];
+        j<cols->sparsity_pattern.rowstart[chunk_row+1] ; ++j)
+      for (unsigned int r=0; r<cols->chunk_size; ++r)
+       for (unsigned int s=0; s<cols->chunk_size; ++s)
+         row_sums(chunk_row * cols->chunk_size + r) +=
+           numbers::NumberTraits<number>::abs(val[j * cols->chunk_size *
+                                                  cols->chunk_size +
+                                                  r * cols->chunk_size +
+                                                  s]);
+
+  return row_sums.linfty_norm();
+}
+
+
+
+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
+                                   // reference to rows or columns
+                                  //
+                                  // padding elements are zero, so we can add
+                                  // them up as well
+  real_type norm_sqr = 0;
+  const unsigned int n_rows = m();
+  for (const number *ptr = &val[0]; ptr != &val[max_len]; ++ptr)
+    norm_sqr +=  numbers::NumberTraits<number>::abs_square(*ptr);
+
+  return std::sqrt (norm_sqr);
+}
+
+
+
+template <typename number>
+template <typename somenumber>
+somenumber
+ChunkSparseMatrix<number>::residual (Vector<somenumber>       &dst,
+                                    const Vector<somenumber> &u,
+                                    const Vector<somenumber> &b) const
+{
+  Assert (cols != 0, ExcNotInitialized());
+  Assert (val != 0, ExcNotInitialized());
+  Assert(m() == dst.size(), ExcDimensionMismatch(m(),dst.size()));
+  Assert(m() == b.size(), ExcDimensionMismatch(m(),b.size()));
+  Assert(n() == u.size(), ExcDimensionMismatch(n(),u.size()));
+
+  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
+  dst = b;
+
+                                  /////////
+                                  // 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
+  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 Vector<somenumber>::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_subtract
+             (cols->chunk_size,
+              val_ptr,
+              u.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] *
+                     u(*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] *
+                       u(*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] *
+                     u(*colnum_ptr * cols->chunk_size + c));
+
+         ++colnum_ptr;
+         val_ptr += cols->chunk_size * cols->chunk_size;
+       }
+      
+
+      dst_ptr += cols->chunk_size;
+    }
+
+                                  // finally compute the norm
+  return dst.l2_norm();
+}
+
+
+
+template <typename number>
+template <typename somenumber>
+void
+ChunkSparseMatrix<number>::precondition_Jacobi (Vector<somenumber>       &dst,
+                                               const Vector<somenumber> &src,
+                                               const number              /*om*/) const
+{
+  Assert (cols != 0, ExcNotInitialized());
+  Assert (val != 0, ExcNotInitialized());
+  Assert (cols->optimize_diagonal(),
+         typename ChunkSparsityPattern::ExcDiagonalNotOptimized());
+  
+  Assert (dst.size() == n(), ExcDimensionMismatch (dst.size(), n()));
+  Assert (src.size() == n(), ExcDimensionMismatch (src.size(), n()));
+
+  Assert (false, ExcNotImplemented());
+}
+
+
+
+template <typename number>
+template <typename somenumber>
+void
+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...
+  Assert (cols != 0, ExcNotInitialized());
+  Assert (val != 0, ExcNotInitialized());
+  Assert (cols->optimize_diagonal(),
+         typename ChunkSparsityPattern::ExcDiagonalNotOptimized());
+  
+  Assert (dst.size() == n(), ExcDimensionMismatch (dst.size(), n()));
+  Assert (src.size() == n(), ExcDimensionMismatch (src.size(), n()));
+
+  Assert (false, ExcNotImplemented());
+}
+
+
+template <typename number>
+template <typename somenumber>
+void
+ChunkSparseMatrix<number>::precondition_SOR (Vector<somenumber>& dst,
+                                            const Vector<somenumber>& src,
+                                            const number om) const
+{
+  Assert (cols != 0, ExcNotInitialized());
+  Assert (val != 0, ExcNotInitialized());
+  Assert (cols->optimize_diagonal(),
+         typename ChunkSparsityPattern::ExcDiagonalNotOptimized());
+  
+
+  dst = src;
+  SOR(dst,om);
+}
+
+
+template <typename number>
+template <typename somenumber>
+void
+ChunkSparseMatrix<number>::precondition_TSOR (Vector<somenumber>& dst,
+                                             const Vector<somenumber>& src,
+                                             const number om) const
+{
+  Assert (cols != 0, ExcNotInitialized());
+  Assert (val != 0, ExcNotInitialized());
+  Assert (cols->optimize_diagonal(),
+         typename ChunkSparsityPattern::ExcDiagonalNotOptimized());
+  
+
+  dst = src;
+  TSOR(dst,om);
+}
+
+
+template <typename number>
+template <typename somenumber>
+void
+ChunkSparseMatrix<number>::SOR (Vector<somenumber>& dst,
+                               const number /*om*/) const
+{
+  Assert (cols != 0, ExcNotInitialized());
+  Assert (val != 0, ExcNotInitialized());
+  Assert (cols->optimize_diagonal(),
+         typename ChunkSparsityPattern::ExcDiagonalNotOptimized());
+  
+  Assert (m() == dst.size(), ExcDimensionMismatch(m(),dst.size()));
+
+  Assert (false, ExcNotImplemented());
+}
+
+
+template <typename number>
+template <typename somenumber>
+void
+ChunkSparseMatrix<number>::TSOR (Vector<somenumber>& dst,
+                                const number /*om*/) const
+{
+  Assert (cols != 0, ExcNotInitialized());
+  Assert (val != 0, ExcNotInitialized());
+  Assert (cols->optimize_diagonal(),
+         typename ChunkSparsityPattern::ExcDiagonalNotOptimized());
+  
+  Assert (m() == dst.size(), ExcDimensionMismatch(m(),dst.size()));
+
+  Assert (false, ExcNotImplemented());
+}
+
+
+template <typename number>
+template <typename somenumber>
+void
+ChunkSparseMatrix<number>::PSOR (Vector<somenumber>& dst,
+                                const std::vector<unsigned int>& permutation,
+                                const std::vector<unsigned int>& inverse_permutation,
+                                const number /*om*/) const
+{
+  Assert (cols != 0, ExcNotInitialized());
+  Assert (val != 0, ExcNotInitialized());
+  Assert (cols->optimize_diagonal(),
+         typename ChunkSparsityPattern::ExcDiagonalNotOptimized());
+  
+  Assert (m() == dst.size(), ExcDimensionMismatch(m(),dst.size()));
+  Assert (m() == permutation.size(),
+         ExcDimensionMismatch(m(), permutation.size()));
+  Assert (m() == inverse_permutation.size(),
+         ExcDimensionMismatch(m(), inverse_permutation.size()));
+
+  Assert (false, ExcNotImplemented());
+}
+
+
+template <typename number>
+template <typename somenumber>
+void
+ChunkSparseMatrix<number>::TPSOR (Vector<somenumber>& dst,
+                                 const std::vector<unsigned int>& permutation,
+                                 const std::vector<unsigned int>& inverse_permutation,
+                                 const number /*om*/) const
+{
+  Assert (cols != 0, ExcNotInitialized());
+  Assert (val != 0, ExcNotInitialized());
+  Assert (cols->optimize_diagonal(),
+         typename ChunkSparsityPattern::ExcDiagonalNotOptimized());
+  
+  Assert (m() == dst.size(), ExcDimensionMismatch(m(),dst.size()));
+  Assert (m() == permutation.size(),
+         ExcDimensionMismatch(m(), permutation.size()));
+  Assert (m() == inverse_permutation.size(),
+         ExcDimensionMismatch(m(), inverse_permutation.size()));
+
+  Assert (false, ExcNotImplemented());
+}
+
+
+
+template <typename number>
+template <typename somenumber>
+void
+ChunkSparseMatrix<number>::SOR_step (Vector<somenumber> &v,
+                                    const Vector<somenumber> &b,
+                                    const number        /*om*/) const
+{
+  Assert (cols != 0, ExcNotInitialized());
+  Assert (val != 0, ExcNotInitialized());
+  Assert (cols->optimize_diagonal(),
+         typename ChunkSparsityPattern::ExcDiagonalNotOptimized());
+  
+  Assert (m() == v.size(), ExcDimensionMismatch(m(),v.size()));
+  Assert (m() == b.size(), ExcDimensionMismatch(m(),b.size()));
+
+  Assert (false, ExcNotImplemented());
+}
+
+
+
+template <typename number>
+template <typename somenumber>
+void
+ChunkSparseMatrix<number>::TSOR_step (Vector<somenumber> &v,
+                                     const Vector<somenumber> &b,
+                                     const number        /*om*/) const
+{
+  Assert (cols != 0, ExcNotInitialized());
+  Assert (val != 0, ExcNotInitialized());
+  Assert (cols->optimize_diagonal(),
+         typename ChunkSparsityPattern::ExcDiagonalNotOptimized());
+  
+  Assert (m() == v.size(), ExcDimensionMismatch(m(),v.size()));
+  Assert (m() == b.size(), ExcDimensionMismatch(m(),b.size()));
+
+  Assert (false, ExcNotImplemented());
+}
+
+
+
+template <typename number>
+template <typename somenumber>
+void
+ChunkSparseMatrix<number>::SSOR_step (Vector<somenumber> &v,
+                                     const Vector<somenumber> &b,
+                                     const number        om) const
+{
+  SOR_step(v,b,om);
+  TSOR_step(v,b,om);
+}
+
+
+
+template <typename number>
+template <typename somenumber>
+void
+ChunkSparseMatrix<number>::SSOR (Vector<somenumber>& dst,
+                                const number /*om*/) const
+{
+  Assert (cols != 0, ExcNotInitialized());
+  Assert (val != 0, ExcNotInitialized());
+  Assert (cols->optimize_diagonal(),
+         typename ChunkSparsityPattern::ExcDiagonalNotOptimized());
+  
+  Assert (m() == dst.size(), ExcDimensionMismatch(m(),dst.size()));
+
+  Assert (false, ExcNotImplemented());
+}
+
+
+
+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
+{
+  AssertThrow (out, ExcIO());
+
+  Assert (cols != 0, ExcNotInitialized());
+  Assert (val != 0, ExcNotInitialized());
+
+  Assert (false, ExcNotImplemented());
+
+  AssertThrow (out, ExcIO());
+}
+
+
+template <typename number>
+void ChunkSparseMatrix<number>::print_formatted (std::ostream &out,
+                                                const unsigned int precision,
+                                                const bool scientific,
+                                                const unsigned int width_,
+                                                const char* zero_string,
+                                                const double denominator) const
+{
+  AssertThrow (out, ExcIO());
+
+  Assert (cols != 0, ExcNotInitialized());
+  Assert (val != 0, ExcNotInitialized());
+
+  unsigned int width = width_;
+  
+  Assert (false, ExcNotImplemented());
+
+  std::ios::fmtflags old_flags = out.flags();
+  unsigned int old_precision = out.precision (precision);
+
+  if (scientific)
+    {
+      out.setf (std::ios::scientific, std::ios::floatfield);
+      if (!width)
+       width = precision+7;
+    } else {
+      out.setf (std::ios::fixed, std::ios::floatfield);
+      if (!width)
+       width = precision+2;
+    }
+
+  for (unsigned int i=0; i<m(); ++i)
+    {
+      for (unsigned int j=0; j<n(); ++j)
+       if (cols->sparsity_pattern(i,j) != SparsityPattern::invalid_entry)
+         out << std::setw(width)
+             << val[cols->sparsity_pattern(i,j)] * denominator << ' ';
+       else
+         out << std::setw(width) << zero_string << ' ';
+      out << std::endl;
+    };
+  AssertThrow (out, ExcIO());
+
+                                  // reset output format
+  out.precision(old_precision);
+  out.flags (old_flags);
+}
+
+
+
+template <typename number>
+void ChunkSparseMatrix<number>::print_pattern (std::ostream &out,
+                                              const double threshold) const
+{
+  AssertThrow (out, ExcIO());
+
+  Assert (cols != 0, ExcNotInitialized());
+  Assert (val != 0, ExcNotInitialized());
+
+  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
+  for (unsigned int i=0; i<cols->sparsity_pattern.n_rows(); ++i)
+    {
+      for (unsigned int d=0; d<chunk_size; ++d)
+       for (unsigned int j=0; j<cols->sparsity_pattern.n_cols(); ++j)
+         if (cols->sparsity_pattern(i,j) == SparsityPattern::invalid_entry)
+           {
+             for (unsigned int e=0; e<chunk_size; ++e)
+               out << '.';
+           }
+         else
+           if (std::fabs(val[cols->sparsity_pattern(i,j)]) > threshold)
+             {
+               for (unsigned int e=0; e<chunk_size; ++e)
+                 out << '*';
+             }
+           else
+             {
+               for (unsigned int e=0; e<chunk_size; ++e)
+                 out << ':';
+             }
+      out << std::endl;
+    }
+
+  AssertThrow (out, ExcIO());
+}
+
+
+
+template <typename number>
+void
+ChunkSparseMatrix<number>::block_write (std::ostream &out) const 
+{
+  AssertThrow (out, ExcIO());
+
+                                   // first the simple objects,
+                                   // bracketed in [...]
+  out << '[' << max_len << "][";
+                                   // then write out real data
+  out.write (reinterpret_cast<const char*>(&val[0]),
+            reinterpret_cast<const char*>(&val[max_len])
+            - reinterpret_cast<const char*>(&val[0]));
+  out << ']';
+  
+  AssertThrow (out, ExcIO());
+}
+
+
+
+template <typename number>
+void
+ChunkSparseMatrix<number>::block_read (std::istream &in)
+{
+  AssertThrow (in, ExcIO());
+
+  char c;
+
+                                   // first read in simple data
+  in >> c;
+  AssertThrow (c == '[', ExcIO());
+  in >> max_len;
+
+  in >> c;
+  AssertThrow (c == ']', ExcIO());
+  in >> c;
+  AssertThrow (c == '[', ExcIO());
+
+                                   // reallocate space
+  delete[] val;
+  val  = new number[max_len];
+  
+                                   // then read data
+  in.read (reinterpret_cast<char*>(&val[0]),
+           reinterpret_cast<char*>(&val[max_len])
+           - reinterpret_cast<char*>(&val[0]));
+  in >> c;
+  AssertThrow (c == ']', ExcIO());
+}
+
+
+
+template <typename number>
+unsigned int
+ChunkSparseMatrix<number>::memory_consumption () const
+{
+  return sizeof(*this) + max_len*sizeof(number);
+}
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
diff --git a/deal.II/lac/source/chunk_sparse_matrix.cc b/deal.II/lac/source/chunk_sparse_matrix.cc
new file mode 100644 (file)
index 0000000..cc28259
--- /dev/null
@@ -0,0 +1,20 @@
+//---------------------------------------------------------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2008 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.
+//
+//---------------------------------------------------------------------------
+
+
+#include <lac/chunk_sparse_matrix.templates.h>
+#include <lac/block_vector.h>
+
+DEAL_II_NAMESPACE_OPEN
+#include "chunk_sparse_matrix.inst"
+DEAL_II_NAMESPACE_CLOSE
diff --git a/deal.II/lac/source/chunk_sparse_matrix.inst.in b/deal.II/lac/source/chunk_sparse_matrix.inst.in
new file mode 100644 (file)
index 0000000..f592b88
--- /dev/null
@@ -0,0 +1,225 @@
+//----------------------------  sparse_matrix.inst.in  ---------------------------
+//    $Id: sparse_matrix_matrix.in.h 15011 2007-08-22 16:59:41Z kanschat $
+//    Version: $Name$
+//
+//    Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2006, 2007, 2008 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.
+//
+//----------------------------  sparse_matrix.inst.in  ---------------------------
+
+
+// real instantiations
+
+for (S : REAL_SCALARS)
+  {
+    template class ChunkSparseMatrix<S>;
+  }
+
+
+
+for (S1, S2 : REAL_SCALARS)
+  {
+    template ChunkSparseMatrix<S1> &
+      ChunkSparseMatrix<S1>::copy_from<S2> (const ChunkSparseMatrix<S2> &);
+
+    template 
+      void ChunkSparseMatrix<S1>::copy_from<S2> (const FullMatrix<S2> &);
+
+    template void ChunkSparseMatrix<S1>::add<S2> (const S1,
+                                            const ChunkSparseMatrix<S2> &);
+  }
+
+
+for (S1, S2 : REAL_SCALARS)
+  {
+    template S2
+      ChunkSparseMatrix<S1>::
+      matrix_norm_square<S2> (const Vector<S2> &) const;
+
+    template S2
+      ChunkSparseMatrix<S1>::
+      matrix_scalar_product<S2> (const Vector<S2> &,
+                                const Vector<S2> &) const;
+
+    template S2 ChunkSparseMatrix<S1>::
+      residual<S2> (Vector<S2> &,
+                   const Vector<S2> &,
+                   const Vector<S2> &) const;
+
+    template void ChunkSparseMatrix<S1>::
+      precondition_SSOR<S2> (Vector<S2> &,
+                            const Vector<S2> &,
+                            const S1) const;
+
+    template void ChunkSparseMatrix<S1>::
+      precondition_SOR<S2> (Vector<S2> &,
+                           const Vector<S2> &,
+                           const S1) const;
+
+    template void ChunkSparseMatrix<S1>::
+      precondition_TSOR<S2> (Vector<S2> &,
+                            const Vector<S2> &,
+                            const S1) const;
+
+    template void ChunkSparseMatrix<S1>::
+      precondition_Jacobi<S2> (Vector<S2> &,
+                              const Vector<S2> &,
+                              const S1) const;
+
+    template void ChunkSparseMatrix<S1>::
+      SOR<S2> (Vector<S2> &,
+              const S1) const;
+    template void ChunkSparseMatrix<S1>::
+      TSOR<S2> (Vector<S2> &,
+               const S1) const;
+    template void ChunkSparseMatrix<S1>::
+      SSOR<S2> (Vector<S2> &,
+               const S1) const;
+    template void ChunkSparseMatrix<S1>::
+      PSOR<S2> (Vector<S2> &,
+               const std::vector<unsigned int>&,
+               const std::vector<unsigned int>&,
+               const S1) const;
+    template void ChunkSparseMatrix<S1>::
+      TPSOR<S2> (Vector<S2> &,
+                const std::vector<unsigned int>&,
+                const std::vector<unsigned int>&,
+                const S1) const;
+    template void ChunkSparseMatrix<S1>::
+      SOR_step<S2> (Vector<S2> &,
+                   const Vector<S2> &,
+                   const S1) const;
+    template void ChunkSparseMatrix<S1>::
+      TSOR_step<S2> (Vector<S2> &,
+                    const Vector<S2> &,
+                    const S1) const;
+    template void ChunkSparseMatrix<S1>::
+      SSOR_step<S2> (Vector<S2> &,
+                    const Vector<S2> &, 
+                    const S1) const;
+  }
+
+
+for (S1, S2, S3 : REAL_SCALARS;
+     V1, V2     : DEAL_II_VEC_TEMPLATES)
+  {
+    template void ChunkSparseMatrix<S1>::
+      vmult (V1<S2> &, const V2<S3> &) const;
+    template void ChunkSparseMatrix<S1>::
+      Tvmult (V1<S2> &, const V2<S3> &) const;
+    template void ChunkSparseMatrix<S1>::
+      vmult_add (V1<S2> &, const V2<S3> &) const;
+    template void ChunkSparseMatrix<S1>::
+      Tvmult_add (V1<S2> &, const V2<S3> &) const;
+  }
+
+
+
+// complex instantiations
+
+// for (S : COMPLEX_SCALARS)
+//   {
+//     template class ChunkSparseMatrix<S>;
+//   }
+
+
+
+// for (S1, S2 : COMPLEX_SCALARS)
+//   {
+//     template ChunkSparseMatrix<S1> &
+//       ChunkSparseMatrix<S1>::copy_from<S2> (const ChunkSparseMatrix<S2> &);
+
+//     template 
+//       void ChunkSparseMatrix<S1>::copy_from<S2> (const FullMatrix<S2> &);
+
+//     template void ChunkSparseMatrix<S1>::add<S2> (const S1,
+//                                          const ChunkSparseMatrix<S2> &);
+//   }
+
+
+// for (S1, S2 : COMPLEX_SCALARS)
+//   {
+//     template S2
+//       ChunkSparseMatrix<S1>::
+//       matrix_norm_square<S2> (const Vector<S2> &) const;
+
+//     template S2
+//       ChunkSparseMatrix<S1>::
+//       matrix_scalar_product<S2> (const Vector<S2> &,
+//                              const Vector<S2> &) const;
+
+//     template S2 ChunkSparseMatrix<S1>::
+//       residual<S2> (Vector<S2> &,
+//                 const Vector<S2> &,
+//                 const Vector<S2> &) const;
+
+//     template void ChunkSparseMatrix<S1>::
+//       precondition_SSOR<S2> (Vector<S2> &,
+//                          const Vector<S2> &,
+//                          const S1) const;
+
+//     template void ChunkSparseMatrix<S1>::
+//       precondition_SOR<S2> (Vector<S2> &,
+//                         const Vector<S2> &,
+//                         const S1) const;
+
+//     template void ChunkSparseMatrix<S1>::
+//       precondition_TSOR<S2> (Vector<S2> &,
+//                          const Vector<S2> &,
+//                          const S1) const;
+
+//     template void ChunkSparseMatrix<S1>::
+//       precondition_Jacobi<S2> (Vector<S2> &,
+//                            const Vector<S2> &,
+//                            const S1) const;
+
+//     template void ChunkSparseMatrix<S1>::
+//       SOR<S2> (Vector<S2> &,
+//            const S1) const;
+//     template void ChunkSparseMatrix<S1>::
+//       TSOR<S2> (Vector<S2> &,
+//             const S1) const;
+//     template void ChunkSparseMatrix<S1>::
+//       SSOR<S2> (Vector<S2> &,
+//             const S1) const;
+//     template void ChunkSparseMatrix<S1>::
+//       PSOR<S2> (Vector<S2> &,
+//             const std::vector<unsigned int>&,
+//             const std::vector<unsigned int>&,
+//             const S1) const;
+//     template void ChunkSparseMatrix<S1>::
+//       TPSOR<S2> (Vector<S2> &,
+//              const std::vector<unsigned int>&,
+//              const std::vector<unsigned int>&,
+//              const S1) const;
+//     template void ChunkSparseMatrix<S1>::
+//       SOR_step<S2> (Vector<S2> &,
+//                 const Vector<S2> &,
+//                 const S1) const;
+//     template void ChunkSparseMatrix<S1>::
+//       TSOR_step<S2> (Vector<S2> &,
+//                  const Vector<S2> &,
+//                  const S1) const;
+//     template void ChunkSparseMatrix<S1>::
+//       SSOR_step<S2> (Vector<S2> &,
+//                  const Vector<S2> &, 
+//                  const S1) const;
+//   }
+
+
+// for (S1, S2, S3 : COMPLEX_SCALARS;
+//      V1, V2     : DEAL_II_VEC_TEMPLATES)
+//   {
+//     template void ChunkSparseMatrix<S1>::
+//       vmult (V1<S2> &, const V2<S3> &) const;
+//     template void ChunkSparseMatrix<S1>::
+//       Tvmult (V1<S2> &, const V2<S3> &) const;
+//     template void ChunkSparseMatrix<S1>::
+//       vmult_add (V1<S2> &, const V2<S3> &) const;
+//     template void ChunkSparseMatrix<S1>::
+//       Tvmult_add (V1<S2> &, const V2<S3> &) const;
+//   }
index 77f9e512ac13a99ebbbb7e4ff3ab9389833905e3..d35b1678639bb377fdb4fc387fe1e71dd3419725 100644 (file)
@@ -1,6 +1,6 @@
 ############################################################
 # Makefile,v 1.15 2002/06/13 12:51:13 hartmann Exp
-# Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007 by the deal.II authors
+# Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008 by the deal.II authors
 ############################################################
 
 ############################################################
@@ -63,6 +63,7 @@ tests_x = grid_generator_?? \
           find_cell_* \
           sparsity_pattern_* \
           sparse_matrix_* \
+          chunk_sparse_matrix_* \
          full_matrix_vector_* \
          solver_* \
          deal_solver_* \
diff --git a/tests/bits/chunk_sparse_matrix_00.cc b/tests/bits/chunk_sparse_matrix_00.cc
new file mode 100644 (file)
index 0000000..ade67fa
--- /dev/null
@@ -0,0 +1,95 @@
+//----------------------------  chunk_sparse_matrix_00.cc  ---------------------------
+//    $Id$
+//    Version: $Name$ 
+//
+//    Copyright (C) 2004, 2005, 2008 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_00.cc  ---------------------------
+
+
+// set a few elements in a chunk sparse matrix and output them again. should
+// yield the same result for all chunk sizes, of course
+
+#include "../tests.h"
+#include <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,false);
+  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 write them to the output stream
+  for (unsigned int i=0; i<m.m(); ++i)
+    {
+      for (unsigned int j=0; j<m.n(); ++j)
+       deallog << std::setprecision(2) << std::fixed << std::setw(4)
+               << m.el(i,j) << ' ';
+
+      deallog << std::endl;
+    }
+}
+
+
+
+int main ()
+{
+  std::ofstream logfile("chunk_sparse_matrix_00/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_00/cmp/generic b/tests/bits/chunk_sparse_matrix_00/cmp/generic
new file mode 100644 (file)
index 0000000..7b3b3ec
--- /dev/null
@@ -0,0 +1,31 @@
+
+DEAL::Chunk size = 1
+DEAL::0    0.50 0    0    0.50 
+DEAL::0    0    1.50 0    0    
+DEAL::0.50 0    0    3.50 0    
+DEAL::0    2.00 0    0    6.50 
+DEAL::0    0    4.50 0    0    
+DEAL::Chunk size = 2
+DEAL::0    0.50 0    0    0.50 
+DEAL::0    0    1.50 0    0    
+DEAL::0.50 0    0    3.50 0    
+DEAL::0    2.00 0    0    6.50 
+DEAL::0    0    4.50 0    0    
+DEAL::Chunk size = 4
+DEAL::0    0.50 0    0    0.50 
+DEAL::0    0    1.50 0    0    
+DEAL::0.50 0    0    3.50 0    
+DEAL::0    2.00 0    0    6.50 
+DEAL::0    0    4.50 0    0    
+DEAL::Chunk size = 5
+DEAL::0    0.50 0    0    0.50 
+DEAL::0    0    1.50 0    0    
+DEAL::0.50 0    0    3.50 0    
+DEAL::0    2.00 0    0    6.50 
+DEAL::0    0    4.50 0    0    
+DEAL::Chunk size = 7
+DEAL::0    0.50 0    0    0.50 
+DEAL::0    0    1.50 0    0    
+DEAL::0.50 0    0    3.50 0    
+DEAL::0    2.00 0    0    6.50 
+DEAL::0    0    4.50 0    0    
diff --git a/tests/bits/chunk_sparse_matrix_01.cc b/tests/bits/chunk_sparse_matrix_01.cc
new file mode 100644 (file)
index 0000000..8ab8713
--- /dev/null
@@ -0,0 +1,98 @@
+//----------------------------  chunk_sparse_matrix_01.cc  ---------------------------
+//    $Id$
+//    Version: $Name$ 
+//
+//    Copyright (C) 2004, 2005, 2008 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_01.cc  ---------------------------
+
+
+// check setting elements in a sparse matrix using
+// ChunkSparseMatrix::set()
+
+#include "../tests.h"
+#include <lac/chunk_sparse_matrix.h>    
+#include <fstream>
+#include <iomanip>
+
+
+void test (const unsigned int chunk_size)
+{
+  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 make sure we retrieve the same ones
+  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)
+        {
+          Assert (m(i,j) == i*j*.5+.5, ExcInternalError());
+          Assert (m.el(i,j) == i*j*.5+.5, ExcInternalError());
+        }
+      else
+        {
+          Assert (m.el(i,j) == 0, ExcInternalError());
+        }
+
+  deallog << "OK" << std::endl;
+}
+
+
+
+int main ()
+{
+  std::ofstream logfile("chunk_sparse_matrix_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_01/cmp/generic b/tests/bits/chunk_sparse_matrix_01/cmp/generic
new file mode 100644 (file)
index 0000000..c9bdf33
--- /dev/null
@@ -0,0 +1,6 @@
+
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
diff --git a/tests/bits/chunk_sparse_matrix_01a.cc b/tests/bits/chunk_sparse_matrix_01a.cc
new file mode 100644 (file)
index 0000000..b490daf
--- /dev/null
@@ -0,0 +1,122 @@
+//----------------------------  chunk_sparse_matrix_01a.cc  ---------------------------
+//    $Id$
+//    Version: $Name$ 
+//
+//    Copyright (C) 2004, 2005, 2008 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_01a.cc  ---------------------------
+
+
+// check setting elements in a sparse matrix using ChunkSparseMatrix::set(). make
+// sure they are correct, and make sure that for the nonexisting entries
+// ChunkSparseMatrix::el() returns zero and operator() throws an exception
+
+#include "../tests.h"
+#include <lac/chunk_sparse_matrix.h>    
+#include <fstream>
+#include <iomanip>
+
+
+void test (const unsigned int chunk_size)
+{
+  ChunkSparsityPattern sp (5,5,3,chunk_size,false);
+  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 make sure we retrieve the same ones
+  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)
+        {
+          Assert (m(i,j) == i*j*.5+.5, ExcInternalError());
+          Assert (m.el(i,j) == i*j*.5+.5, ExcInternalError());
+        }
+      else
+        {
+                                          // reading elements not in the
+                                          // sparsity pattern should return
+                                          // zero
+         const double x = m.el(i,j);
+         Assert (x == 0, ExcInternalError());
+
+                                          // if this is a sparsity_pattern
+                                          // with chunk_size==1, then we need
+                                          // to get an exception if we access
+                                          // any other element. if
+                                          // chunk_size>1, then this isn't
+                                          // necessarily true
+          bool exc_thrown = false;
+          double d;
+          try
+            {
+              d = m(i,j);
+            }
+          catch (const std::exception &)
+            {
+              exc_thrown = true;
+            }
+          Assert ((exc_thrown == true) || (chunk_size > 1),
+                 ExcInternalError());
+        }
+
+  deallog << "OK" << std::endl;
+}
+
+
+
+int main ()
+{
+  std::ofstream logfile("chunk_sparse_matrix_01a/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_01a/cmp/generic b/tests/bits/chunk_sparse_matrix_01a/cmp/generic
new file mode 100644 (file)
index 0000000..c9bdf33
--- /dev/null
@@ -0,0 +1,6 @@
+
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
diff --git a/tests/bits/chunk_sparse_matrix_02.cc b/tests/bits/chunk_sparse_matrix_02.cc
new file mode 100644 (file)
index 0000000..6a13460
--- /dev/null
@@ -0,0 +1,97 @@
+//----------------------------  chunk_sparse_matrix_02.cc  ---------------------------
+//    $Id$
+//    Version: $Name$ 
+//
+//    Copyright (C) 2004, 2005, 2008 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_02.cc  ---------------------------
+
+
+// check setting elements in a matrix using
+// ChunkSparseMatrix::add()
+
+#include "../tests.h"
+#include <lac/chunk_sparse_matrix.h>    
+#include <fstream>
+
+
+void test (const unsigned int chunk_size)
+{
+  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.add (i,j, i*j*.5+.5);
+
+                                   // then make sure we retrieve the same ones
+  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)
+        {
+          Assert (m(i,j) == i*j*.5+.5, ExcInternalError());
+          Assert (m.el(i,j) == i*j*.5+.5, ExcInternalError());
+        }
+      else
+        {
+          Assert (m.el(i,j) == 0, ExcInternalError());
+        }
+
+  deallog << "OK" << std::endl;
+}
+
+
+
+int main ()
+{
+  std::ofstream logfile("chunk_sparse_matrix_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_02/cmp/generic b/tests/bits/chunk_sparse_matrix_02/cmp/generic
new file mode 100644 (file)
index 0000000..c9bdf33
--- /dev/null
@@ -0,0 +1,6 @@
+
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
diff --git a/tests/bits/chunk_sparse_matrix_03.cc b/tests/bits/chunk_sparse_matrix_03.cc
new file mode 100644 (file)
index 0000000..54c477d
--- /dev/null
@@ -0,0 +1,104 @@
+//----------------------------  chunk_sparse_matrix_03.cc  ---------------------------
+//    $Id$
+//    Version: $Name$ 
+//
+//    Copyright (C) 2004, 2005, 2008 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_03.cc  ---------------------------
+
+
+// check setting elements in a sparse matrix using set() and add()
+// intermixed. this poses PETSc some problems, since one has to flush some
+// buffer in between these two types of operations, but it shouldn't be a
+// problem with the deal.II matrices. worth checking anyway
+
+#include "../tests.h"
+#include <lac/chunk_sparse_matrix.h>    
+#include <fstream>
+
+
+void test (const unsigned int chunk_size)
+{
+  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)
+        {
+          if (i*j % 2 == 0)
+            m.set (i,j, i*j*.5+.5);
+          else
+            m.add (i,j, i*j*.5+.5);
+        }
+  
+                                   // then make sure we retrieve the same ones
+  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)
+        {
+          Assert (m(i,j) == i*j*.5+.5, ExcInternalError());
+          Assert (m.el(i,j) == i*j*.5+.5, ExcInternalError());
+        }
+      else
+        {
+          Assert (m.el(i,j) == 0, ExcInternalError());
+        }
+
+  deallog << "OK" << std::endl;
+}
+
+
+
+int main ()
+{
+  std::ofstream logfile("chunk_sparse_matrix_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_03/cmp/generic b/tests/bits/chunk_sparse_matrix_03/cmp/generic
new file mode 100644 (file)
index 0000000..c9bdf33
--- /dev/null
@@ -0,0 +1,6 @@
+
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
diff --git a/tests/bits/chunk_sparse_matrix_03a.cc b/tests/bits/chunk_sparse_matrix_03a.cc
new file mode 100644 (file)
index 0000000..afc8289
--- /dev/null
@@ -0,0 +1,108 @@
+//----------------------------  chunk_sparse_matrix_03.cc  ---------------------------
+//    $Id$
+//    Version: $Name$ 
+//
+//    Copyright (C) 2004, 2005, 2008 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_03.cc  ---------------------------
+
+
+// check setting elements in a sparse matrix using set() and add()
+// intermixed. this poses PETSc some problems, since one has to flush some
+// buffer in between these two types of operations, but it shouldn't be a
+// problem with the deal.II matrices. worth checking anyway
+//
+// in contrast to petsc_03, we set and add the same elements here twice, to
+// get double the original value
+
+#include "../tests.h"
+#include <lac/chunk_sparse_matrix.h>    
+#include <fstream>
+
+
+void test (const unsigned int chunk_size)
+{
+  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 add the same elements again
+  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.add (i,j, i*j*.5+.5);
+  
+                                   // then make sure we retrieve the correct
+                                   // ones
+  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)
+        {
+          Assert (m(i,j) == 2*(i*j*.5+.5), ExcInternalError());
+          Assert (m.el(i,j) == 2*(i*j*.5+.5), ExcInternalError());
+        }
+      else
+        {
+          Assert (m.el(i,j) == 0, ExcInternalError());
+        }
+
+  deallog << "OK" << std::endl;
+}
+
+
+
+int main ()
+{
+  std::ofstream logfile("chunk_sparse_matrix_03a/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_03a/cmp/generic b/tests/bits/chunk_sparse_matrix_03a/cmp/generic
new file mode 100644 (file)
index 0000000..c9bdf33
--- /dev/null
@@ -0,0 +1,6 @@
+
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
diff --git a/tests/bits/chunk_sparse_matrix_03b.cc b/tests/bits/chunk_sparse_matrix_03b.cc
new file mode 100644 (file)
index 0000000..a2c8b24
--- /dev/null
@@ -0,0 +1,108 @@
+//----------------------------  chunk_sparse_matrix_03.cc  ---------------------------
+//    $Id$
+//    Version: $Name$ 
+//
+//    Copyright (C) 2004, 2005, 2008 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_03.cc  ---------------------------
+
+
+// check setting elements in a sparse matrix using set() and add()
+// intermixed. this poses PETSc some problems, since one has to flush some
+// buffer in between these two types of operations, but it shouldn't be a
+// problem with the deal.II matrices. worth checking anyway
+//
+// in contrast to petsc_03, we set and add the same elements here twice, to
+// get double the original value
+
+#include "../tests.h"
+#include <lac/chunk_sparse_matrix.h>    
+#include <fstream>
+
+
+void test (const unsigned int chunk_size)
+{
+  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 add the same elements again
+  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.add (i,j, i*j*.5+.5);
+  
+                                   // then make sure we retrieve the correct
+                                   // ones
+  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)
+        {
+          Assert (m(i,j) == 2*(i*j*.5+.5), ExcInternalError());
+          Assert (m.el(i,j) == 2*(i*j*.5+.5), ExcInternalError());
+        }
+      else
+        {
+          Assert (m.el(i,j) == 0, ExcInternalError());
+        }
+
+  deallog << "OK" << std::endl;
+}
+
+
+
+int main ()
+{
+  std::ofstream logfile("chunk_sparse_matrix_03b/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_03b/cmp/generic b/tests/bits/chunk_sparse_matrix_03b/cmp/generic
new file mode 100644 (file)
index 0000000..c9bdf33
--- /dev/null
@@ -0,0 +1,6 @@
+
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
diff --git a/tests/bits/chunk_sparse_matrix_04.cc b/tests/bits/chunk_sparse_matrix_04.cc
new file mode 100644 (file)
index 0000000..bb65fad
--- /dev/null
@@ -0,0 +1,80 @@
+//----------------------------  chunk_sparse_matrix_04.cc  ---------------------------
+//    $Id$
+//    Version: $Name$ 
+//
+//    Copyright (C) 2004, 2005, 2008 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_04.cc  ---------------------------
+
+
+// check querying matrix sizes
+
+#include "../tests.h"
+#include <lac/chunk_sparse_matrix.h>    
+#include <fstream>
+
+
+void test (const unsigned int chunk_size)
+{
+  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);
+  
+  Assert (m.m() == 5, ExcInternalError());
+  Assert (m.n() == 5, ExcInternalError());
+
+  deallog << "OK" << std::endl;
+}
+
+
+
+int main ()
+{
+  std::ofstream logfile("chunk_sparse_matrix_04/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_04/cmp/generic b/tests/bits/chunk_sparse_matrix_04/cmp/generic
new file mode 100644 (file)
index 0000000..c9bdf33
--- /dev/null
@@ -0,0 +1,6 @@
+
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
diff --git a/tests/bits/chunk_sparse_matrix_05.cc b/tests/bits/chunk_sparse_matrix_05.cc
new file mode 100644 (file)
index 0000000..e2dcbf5
--- /dev/null
@@ -0,0 +1,105 @@
+//----------------------------  chunk_sparse_matrix_05.cc  ---------------------------
+//    $Id$
+//    Version: $Name$ 
+//
+//    Copyright (C) 2004, 2005, 2008 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_05.cc  ---------------------------
+
+
+// check querying the number of nonzero elements in
+// ChunkSparseMatrix
+
+#include "../tests.h"
+#include <lac/chunk_sparse_matrix.h>    
+#include <fstream>
+
+
+void test (const unsigned int chunk_size)
+{
+  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. count how many
+                                   // entries we have. note that for square
+                                   // matrices we also always store the
+                                   // diagonal element, so add one per row,
+                                   // but don't count it when traversing the
+                                   // row
+  unsigned int counter = 0;
+  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);
+            if (i!=j)
+              ++counter;
+          }
+      ++counter;
+    }
+
+  deallog << m.n_nonzero_elements() << std::endl;
+  if (chunk_size == 1)
+    Assert (m.n_nonzero_elements() == counter,
+           ExcInternalError())
+  else
+    Assert (m.n_nonzero_elements() >= counter,
+           ExcInternalError());
+  
+  deallog << "OK" << std::endl;
+}
+
+
+
+int main ()
+{
+  std::ofstream logfile("chunk_sparse_matrix_05/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_05/cmp/generic b/tests/bits/chunk_sparse_matrix_05/cmp/generic
new file mode 100644 (file)
index 0000000..1af8241
--- /dev/null
@@ -0,0 +1,11 @@
+
+DEAL::13
+DEAL::OK
+DEAL::23
+DEAL::OK
+DEAL::25
+DEAL::OK
+DEAL::25
+DEAL::OK
+DEAL::25
+DEAL::OK
diff --git a/tests/bits/chunk_sparse_matrix_05a.cc b/tests/bits/chunk_sparse_matrix_05a.cc
new file mode 100644 (file)
index 0000000..77422f7
--- /dev/null
@@ -0,0 +1,102 @@
+//----------------------------  chunk_sparse_matrix_05a.cc  ---------------------------
+//    $Id$
+//    Version: $Name$ 
+//
+//    Copyright (C) 2004, 2005, 2008 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_05a.cc  ---------------------------
+
+
+// check querying the number of nonzero elements in
+// ChunkSparseMatrix when we don't store the diagonal elements explicitly
+
+#include "../tests.h"
+#include <lac/chunk_sparse_matrix.h>    
+#include <fstream>
+
+
+void test (const unsigned int chunk_size)
+{
+  ChunkSparsityPattern sp (5,5,3,chunk_size,false);
+  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. count how many
+                                   // entries we have. note that for square
+                                   // matrices we also always store the
+                                   // diagonal element, except when as above
+                                   // we set the special flag for the matrix
+                                   // sparsity pattern
+  unsigned int counter = 0;
+  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);
+          ++counter;
+        }
+
+  deallog << m.n_nonzero_elements() << std::endl;
+
+  if (chunk_size == 1)
+    Assert (m.n_nonzero_elements() == counter,
+           ExcInternalError())
+  else
+    Assert (m.n_nonzero_elements() >= counter,
+           ExcInternalError());
+  
+  deallog << "OK" << std::endl;
+}
+
+
+
+int main ()
+{
+  std::ofstream logfile("chunk_sparse_matrix_05a/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_05a/cmp/generic b/tests/bits/chunk_sparse_matrix_05a/cmp/generic
new file mode 100644 (file)
index 0000000..5f8ca8f
--- /dev/null
@@ -0,0 +1,11 @@
+
+DEAL::8
+DEAL::OK
+DEAL::22
+DEAL::OK
+DEAL::24
+DEAL::OK
+DEAL::25
+DEAL::OK
+DEAL::25
+DEAL::OK
diff --git a/tests/bits/chunk_sparse_matrix_06.cc b/tests/bits/chunk_sparse_matrix_06.cc
new file mode 100644 (file)
index 0000000..8e2c9ac
--- /dev/null
@@ -0,0 +1,89 @@
+//----------------------------  chunk_sparse_matrix_06.cc  ---------------------------
+//    $Id$
+//    Version: $Name$ 
+//
+//    Copyright (C) 2004, 2005, 2008 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_06.cc  ---------------------------
+
+
+// check ChunkSparseMatrix::l1_norm
+
+#include "../tests.h"
+#include <lac/chunk_sparse_matrix.h>    
+#include <fstream>
+
+
+void test (const unsigned int chunk_size)
+{
+  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. count how many
+                                   // entries we have
+  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);  
+
+                                   // compare against the exact value of the
+                                   // l1-norm (max col-sum)
+  deallog << m.l1_norm() << std::endl;
+  Assert (m.l1_norm() == 7, ExcInternalError());
+  
+  deallog << "OK" << std::endl;
+}
+
+
+
+int main ()
+{
+  std::ofstream logfile("chunk_sparse_matrix_06/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_06/cmp/generic b/tests/bits/chunk_sparse_matrix_06/cmp/generic
new file mode 100644 (file)
index 0000000..fdea211
--- /dev/null
@@ -0,0 +1,11 @@
+
+DEAL::7.00000
+DEAL::OK
+DEAL::7.00000
+DEAL::OK
+DEAL::7.00000
+DEAL::OK
+DEAL::7.00000
+DEAL::OK
+DEAL::7.00000
+DEAL::OK
diff --git a/tests/bits/chunk_sparse_matrix_07.cc b/tests/bits/chunk_sparse_matrix_07.cc
new file mode 100644 (file)
index 0000000..8946908
--- /dev/null
@@ -0,0 +1,89 @@
+//----------------------------  chunk_sparse_matrix_07.cc  ---------------------------
+//    $Id$
+//    Version: $Name$ 
+//
+//    Copyright (C) 2004, 2005, 2008 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_07.cc  ---------------------------
+
+
+// check ChunkSparseMatrix::linfty_norm
+
+#include "../tests.h"
+#include <lac/chunk_sparse_matrix.h>    
+#include <fstream>
+
+
+void test (const unsigned int chunk_size)
+{
+  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. count how many
+                                   // entries we have
+  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);  
+
+                                   // compare against the exact value of the
+                                   // linfty-norm (max row-sum)
+  deallog << m.linfty_norm() << std::endl;
+  Assert (m.linfty_norm() == 8.5, ExcInternalError());
+  
+  deallog << "OK" << std::endl;
+}
+
+
+
+int main ()
+{
+  std::ofstream logfile("chunk_sparse_matrix_07/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_07/cmp/generic b/tests/bits/chunk_sparse_matrix_07/cmp/generic
new file mode 100644 (file)
index 0000000..4514d74
--- /dev/null
@@ -0,0 +1,11 @@
+
+DEAL::8.50000
+DEAL::OK
+DEAL::8.50000
+DEAL::OK
+DEAL::8.50000
+DEAL::OK
+DEAL::8.50000
+DEAL::OK
+DEAL::8.50000
+DEAL::OK
diff --git a/tests/bits/chunk_sparse_matrix_08.cc b/tests/bits/chunk_sparse_matrix_08.cc
new file mode 100644 (file)
index 0000000..10ede05
--- /dev/null
@@ -0,0 +1,94 @@
+//----------------------------  chunk_sparse_matrix_08.cc  ---------------------------
+//    chunk_sparse_matrix_08.cc,v 1.4 2004/02/26 17:25:44 wolf Exp
+//    Version:  
+//
+//    Copyright (C) 2004, 2005, 2008 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_08.cc  ---------------------------
+
+
+// check ChunkSparseMatrix::frobenius_norm
+
+#include "../tests.h"
+#include <lac/chunk_sparse_matrix.h>    
+#include <fstream>
+
+
+void test (const unsigned int chunk_size)
+{
+  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. count how many
+                                   // entries we have
+  double norm = 0;
+  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);
+          norm += (i*j*.5+.5)*(i*j*.5+.5);
+        }
+  norm = std::sqrt(norm);
+
+                                   // compare against the exact value of the
+                                   // l2-norm (max row-sum)
+  deallog << m.frobenius_norm() << std::endl;
+  Assert (std::fabs((m.frobenius_norm() - norm)/norm) < 1e-14, ExcInternalError());
+  
+  deallog << "OK" << std::endl;
+}
+
+
+
+int main ()
+{
+  std::ofstream logfile("chunk_sparse_matrix_08/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_08/cmp/generic b/tests/bits/chunk_sparse_matrix_08/cmp/generic
new file mode 100644 (file)
index 0000000..4eb71e4
--- /dev/null
@@ -0,0 +1,11 @@
+
+DEAL::9.04157
+DEAL::OK
+DEAL::9.04157
+DEAL::OK
+DEAL::9.04157
+DEAL::OK
+DEAL::9.04157
+DEAL::OK
+DEAL::9.04157
+DEAL::OK
diff --git a/tests/bits/chunk_sparse_matrix_09.cc b/tests/bits/chunk_sparse_matrix_09.cc
new file mode 100644 (file)
index 0000000..dbc906b
--- /dev/null
@@ -0,0 +1,100 @@
+//----------------------------  chunk_sparse_matrix_09.cc  ---------------------------
+//    $Id$
+//    Version: $Name$ 
+//
+//    Copyright (C) 2004, 2005, 2008 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_09.cc  ---------------------------
+
+
+// check ChunkSparseMatrix::operator *=
+
+#include "../tests.h"
+#include <lac/chunk_sparse_matrix.h>    
+#include <fstream>
+
+
+void test (const unsigned int chunk_size)
+{
+  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 multiply everything by 1.25 and
+                                   // make sure we retrieve the values we
+                                   // expect
+  m *= 1.25;
+  
+  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)
+        {
+          Assert (m(i,j) == (i*j*.5+.5)*1.25, ExcInternalError());
+          Assert (m.el(i,j) == (i*j*.5+.5)*1.25, ExcInternalError());
+        }
+      else
+        {
+          Assert (m.el(i,j) == 0, ExcInternalError());
+        }
+
+  deallog << "OK" << std::endl;
+}
+
+
+
+int main ()
+{
+  std::ofstream logfile("chunk_sparse_matrix_09/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_09/cmp/generic b/tests/bits/chunk_sparse_matrix_09/cmp/generic
new file mode 100644 (file)
index 0000000..c9bdf33
--- /dev/null
@@ -0,0 +1,6 @@
+
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
diff --git a/tests/bits/chunk_sparse_matrix_10.cc b/tests/bits/chunk_sparse_matrix_10.cc
new file mode 100644 (file)
index 0000000..57e57c8
--- /dev/null
@@ -0,0 +1,100 @@
+//----------------------------  chunk_sparse_matrix_10.cc  ---------------------------
+//    $Id$
+//    Version: $Name$ 
+//
+//    Copyright (C) 2004, 2005, 2008 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_10.cc  ---------------------------
+
+
+// check ChunkSparseMatrix::operator /=
+
+#include "../tests.h"
+#include <lac/chunk_sparse_matrix.h>    
+#include <fstream>
+
+
+void test (const unsigned int chunk_size)
+{
+  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 divide everything by 4/3 and
+                                   // make sure we retrieve the values we
+                                   // expect
+  m /= 4./3.;
+  
+  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)
+        {
+          Assert (m(i,j) == (i*j*.5+.5)/4*3, ExcInternalError());
+          Assert (m.el(i,j) == (i*j*.5+.5)/4*3, ExcInternalError());
+        }
+      else
+        {
+          Assert (m.el(i,j) == 0, ExcInternalError());
+        }
+
+  deallog << "OK" << std::endl;
+}
+
+
+
+int main ()
+{
+  std::ofstream logfile("chunk_sparse_matrix_10/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_10/cmp/generic b/tests/bits/chunk_sparse_matrix_10/cmp/generic
new file mode 100644 (file)
index 0000000..c9bdf33
--- /dev/null
@@ -0,0 +1,6 @@
+
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
diff --git a/tests/bits/chunk_sparse_matrix_vector_01.cc b/tests/bits/chunk_sparse_matrix_vector_01.cc
new file mode 100644 (file)
index 0000000..8bd552f
--- /dev/null
@@ -0,0 +1,110 @@
+//----------------------------  chunk_sparse_matrix_vector_01.cc  ---------------------------
+//    $Id$
+//    Version: $Name$ 
+//
+//    Copyright (C) 2004, 2005, 2008 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_vector_01.cc  ---------------------------
+
+
+// check ChunkSparseMatrix::vmult
+
+#include "../tests.h"
+#include <lac/vector.h>
+#include <lac/chunk_sparse_matrix.h>
+#include <fstream>
+#include <iomanip>
+#include <vector>
+
+
+void test (const unsigned int chunk_size,
+          Vector<double> &v,
+           Vector<double> &w)
+{
+                                   // set some entries in the
+                                   // matrix. actually, set them all
+  ChunkSparsityPattern sp (v.size(),v.size(),v.size(), chunk_size);
+  for (unsigned int i=0; i<v.size(); ++i)
+    for (unsigned int j=0; j<v.size(); ++j)
+      sp.add (i,j);
+  sp.compress ();
+
+                                   // then create a matrix from that
+  ChunkSparseMatrix<double> m(sp);
+  for (unsigned int i=0; i<m.m(); ++i)
+    for (unsigned int j=0; j<m.n(); ++j)
+        m.set (i,j, i+2*j);
+
+  for (unsigned int i=0; i<v.size(); ++i)
+    v(i) = i;
+  
+  v.compress ();
+  w.compress ();
+
+                                   // w:=Mv
+  m.vmult (w,v);
+
+                                   // make sure we get the expected result
+  for (unsigned int i=0; i<v.size(); ++i)
+    {
+      Assert (v(i) == i, ExcInternalError());
+
+      double result = 0;
+      for (unsigned int j=0; j<m.n(); ++j)
+        result += (i+2*j)*j;
+      Assert (w(i) == result, ExcInternalError());
+    }
+
+  deallog << "OK" << std::endl;
+}
+
+
+
+int main () 
+{
+  std::ofstream logfile("chunk_sparse_matrix_vector_01/output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  try
+    {
+      const unsigned int chunk_sizes[] = { 1, 2, 4, 7, 11 };
+      for (unsigned int i=0; i<sizeof(chunk_sizes)/sizeof(chunk_sizes[0]);
+          ++i)
+       {
+         Vector<double> v (100);
+         Vector<double> w (100);
+         test (chunk_sizes[i], v,w);
+       }
+    }
+  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_vector_01/cmp/generic b/tests/bits/chunk_sparse_matrix_vector_01/cmp/generic
new file mode 100644 (file)
index 0000000..c9bdf33
--- /dev/null
@@ -0,0 +1,6 @@
+
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
diff --git a/tests/bits/chunk_sparse_matrix_vector_02.cc b/tests/bits/chunk_sparse_matrix_vector_02.cc
new file mode 100644 (file)
index 0000000..318b33d
--- /dev/null
@@ -0,0 +1,110 @@
+//----------------------------  chunk_sparse_matrix_vector_02.cc  ---------------------------
+//    $Id$
+//    Version: $Name$ 
+//
+//    Copyright (C) 2004, 2005, 2008 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_vector_02.cc  ---------------------------
+
+
+// check ChunkSparseMatrix::Tvmult
+
+#include "../tests.h"
+#include <lac/vector.h>
+#include <lac/chunk_sparse_matrix.h>
+#include <fstream>
+#include <iomanip>
+#include <vector>
+
+
+void test (const unsigned int chunk_size,
+          Vector<double> &v,
+           Vector<double> &w)
+{
+                                   // set some entries in the
+                                   // matrix. actually, set them all
+  ChunkSparsityPattern sp (v.size(),v.size(),v.size(), chunk_size);
+  for (unsigned int i=0; i<v.size(); ++i)
+    for (unsigned int j=0; j<v.size(); ++j)
+      sp.add (i,j);
+  sp.compress ();
+
+                                   // then create a matrix from that
+  ChunkSparseMatrix<double> m(sp);
+  for (unsigned int i=0; i<m.m(); ++i)
+    for (unsigned int j=0; j<m.n(); ++j)
+        m.set (i,j, i+2*j);
+
+  for (unsigned int i=0; i<v.size(); ++i)
+    v(i) = i;
+  
+  v.compress ();
+  w.compress ();
+
+                                   // w:=Mv
+  m.Tvmult (w,v);
+
+                                   // make sure we get the expected result
+  for (unsigned int i=0; i<v.size(); ++i)
+    {
+      Assert (v(i) == i, ExcInternalError());
+
+      double result = 0;
+      for (unsigned int j=0; j<m.n(); ++j)
+        result += (j+2*i)*j;
+      Assert (w(i) == result, ExcInternalError());
+    }
+
+  deallog << "OK" << std::endl;
+}
+
+
+
+int main () 
+{
+  std::ofstream logfile("chunk_sparse_matrix_vector_02/output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  try
+    {
+      const unsigned int chunk_sizes[] = { 1, 2, 4, 7, 11 };
+      for (unsigned int i=0; i<sizeof(chunk_sizes)/sizeof(chunk_sizes[0]);
+          ++i)
+       {
+         Vector<double> v (100);
+         Vector<double> w (100);
+         test (chunk_sizes[i], v,w);
+       }
+    }
+  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_vector_02/cmp/generic b/tests/bits/chunk_sparse_matrix_vector_02/cmp/generic
new file mode 100644 (file)
index 0000000..c9bdf33
--- /dev/null
@@ -0,0 +1,6 @@
+
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
diff --git a/tests/bits/chunk_sparse_matrix_vector_03.cc b/tests/bits/chunk_sparse_matrix_vector_03.cc
new file mode 100644 (file)
index 0000000..59e19be
--- /dev/null
@@ -0,0 +1,113 @@
+//----------------------------  chunk_sparse_matrix_vector_03.cc  ---------------------------
+//    $Id$
+//    Version: $Name$ 
+//
+//    Copyright (C) 2004, 2005, 2008 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_vector_03.cc  ---------------------------
+
+
+// check ChunkSparseMatrix::vmult_add
+
+#include "../tests.h"
+#include <lac/vector.h>
+#include <lac/chunk_sparse_matrix.h>
+#include <fstream>
+#include <iomanip>
+#include <vector>
+
+
+void test (const unsigned int chunk_size,
+          Vector<double> &v,
+           Vector<double> &w)
+{
+                                   // set some entries in the
+                                   // matrix. actually, set them all
+  ChunkSparsityPattern sp (v.size(),v.size(),v.size(), chunk_size);
+  for (unsigned int i=0; i<v.size(); ++i)
+    for (unsigned int j=0; j<v.size(); ++j)
+      sp.add (i,j);
+  sp.compress ();
+
+                                   // then create a matrix from that
+  ChunkSparseMatrix<double> m(sp);
+  for (unsigned int i=0; i<m.m(); ++i)
+    for (unsigned int j=0; j<m.n(); ++j)
+        m.set (i,j, i+2*j);
+
+  for (unsigned int i=0; i<v.size(); ++i)
+    {
+      v(i) = i;
+      w(i) = i;
+    }
+  
+  v.compress ();
+  w.compress ();
+
+                                   // w+=Mv
+  m.vmult_add (w,v);
+
+                                   // make sure we get the expected result
+  for (unsigned int i=0; i<v.size(); ++i)
+    {
+      Assert (v(i) == i, ExcInternalError());
+
+      double result = 0;
+      for (unsigned int j=0; j<m.n(); ++j)
+        result += (i+2*j)*j;
+      Assert (w(i) == i+result, ExcInternalError());
+    }
+
+  deallog << "OK" << std::endl;
+}
+
+
+
+int main () 
+{
+  std::ofstream logfile("chunk_sparse_matrix_vector_03/output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  try
+    {
+      const unsigned int chunk_sizes[] = { 1, 2, 4, 7, 11 };
+      for (unsigned int i=0; i<sizeof(chunk_sizes)/sizeof(chunk_sizes[0]);
+          ++i)
+       {
+         Vector<double> v (100);
+         Vector<double> w (100);
+         test (chunk_sizes[i], v,w);
+       }
+    }
+  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_vector_03/cmp/generic b/tests/bits/chunk_sparse_matrix_vector_03/cmp/generic
new file mode 100644 (file)
index 0000000..c9bdf33
--- /dev/null
@@ -0,0 +1,6 @@
+
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
diff --git a/tests/bits/chunk_sparse_matrix_vector_04.cc b/tests/bits/chunk_sparse_matrix_vector_04.cc
new file mode 100644 (file)
index 0000000..9f85cb0
--- /dev/null
@@ -0,0 +1,113 @@
+//----------------------------  chunk_sparse_matrix_vector_04.cc  ---------------------------
+//    $Id$
+//    Version: $Name$ 
+//
+//    Copyright (C) 2004, 2005, 2008 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_vector_04.cc  ---------------------------
+
+
+// check ChunkSparseMatrix::Tvmult_add
+
+#include "../tests.h"
+#include <lac/vector.h>
+#include <lac/chunk_sparse_matrix.h>
+#include <fstream>
+#include <iomanip>
+#include <vector>
+
+
+void test (const unsigned int chunk_size,
+          Vector<double> &v,
+           Vector<double> &w)
+{
+                                   // set some entries in the
+                                   // matrix. actually, set them all
+  ChunkSparsityPattern sp (v.size(),v.size(),v.size(), chunk_size);
+  for (unsigned int i=0; i<v.size(); ++i)
+    for (unsigned int j=0; j<v.size(); ++j)
+      sp.add (i,j);
+  sp.compress ();
+
+                                   // then create a matrix from that
+  ChunkSparseMatrix<double> m(sp);
+  for (unsigned int i=0; i<m.m(); ++i)
+    for (unsigned int j=0; j<m.n(); ++j)
+        m.set (i,j, i+2*j);
+
+  for (unsigned int i=0; i<v.size(); ++i)
+    {
+      v(i) = i;
+      w(i) = i;
+    }
+  
+  v.compress ();
+  w.compress ();
+
+                                   // w:=Mv
+  m.Tvmult_add (w,v);
+
+                                   // make sure we get the expected result
+  for (unsigned int i=0; i<v.size(); ++i)
+    {
+      Assert (v(i) == i, ExcInternalError());
+
+      double result = 0;
+      for (unsigned int j=0; j<m.n(); ++j)
+        result += (j+2*i)*j;
+      Assert (w(i) == i+result, ExcInternalError());
+    }
+
+  deallog << "OK" << std::endl;
+}
+
+
+
+int main () 
+{
+  std::ofstream logfile("chunk_sparse_matrix_vector_04/output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  try
+    {
+      const unsigned int chunk_sizes[] = { 1, 2, 4, 7, 11 };
+      for (unsigned int i=0; i<sizeof(chunk_sizes)/sizeof(chunk_sizes[0]);
+          ++i)
+       {
+         Vector<double> v (100);
+         Vector<double> w (100);
+         test (chunk_sizes[i], v,w);
+       }
+    }
+  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_vector_04/cmp/generic b/tests/bits/chunk_sparse_matrix_vector_04/cmp/generic
new file mode 100644 (file)
index 0000000..c9bdf33
--- /dev/null
@@ -0,0 +1,6 @@
+
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
diff --git a/tests/bits/chunk_sparse_matrix_vector_05.cc b/tests/bits/chunk_sparse_matrix_vector_05.cc
new file mode 100644 (file)
index 0000000..4f2efc4
--- /dev/null
@@ -0,0 +1,116 @@
+//----------------------------  chunk_sparse_matrix_vector_05.cc  ---------------------------
+//    $Id$
+//    Version: $Name$ 
+//
+//    Copyright (C) 2004, 2005, 2008 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_vector_05.cc  ---------------------------
+
+
+// check ChunkSparseMatrix::matrix_scalar_product
+
+#include "../tests.h"
+#include <lac/vector.h>
+#include <lac/chunk_sparse_matrix.h>
+#include <fstream>
+#include <iomanip>
+#include <vector>
+
+
+void test (const unsigned int chunk_size,
+          Vector<double> &v,
+           Vector<double> &w)
+{
+                                   // set some entries in the
+                                   // matrix. actually, set them all
+  ChunkSparsityPattern sp (v.size(),v.size(),v.size(), chunk_size);
+  for (unsigned int i=0; i<v.size(); ++i)
+    for (unsigned int j=0; j<v.size(); ++j)
+      sp.add (i,j);
+  sp.compress ();
+
+                                   // then create a matrix from that
+  ChunkSparseMatrix<double> m(sp);
+  for (unsigned int i=0; i<m.m(); ++i)
+    for (unsigned int j=0; j<m.n(); ++j)
+        m.set (i,j, i+2*j);
+
+  for (unsigned int i=0; i<v.size(); ++i)
+    {
+      v(i) = i;
+      w(i) = i+1;
+    }
+  
+  v.compress ();
+  w.compress ();
+
+                                   // <w,Mv>
+  const double s = m.matrix_scalar_product (w,v);
+
+                                   // make sure we get the expected result
+  for (unsigned int i=0; i<v.size(); ++i)
+    {
+      Assert (v(i) == i, ExcInternalError());
+      Assert (w(i) == i+1, ExcInternalError());
+    }
+
+  double result = 0;
+  for (unsigned int i=0; i<m.m(); ++i)
+    for (unsigned int j=0; j<m.n(); ++j)
+      result += (i+2*j)*j*(i+1);
+
+  Assert (s == result, ExcInternalError());
+
+  deallog << "OK" << std::endl;
+}
+
+
+
+int main () 
+{
+  std::ofstream logfile("chunk_sparse_matrix_vector_05/output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  try
+    {
+      const unsigned int chunk_sizes[] = { 1, 2, 4, 7, 11 };
+      for (unsigned int i=0; i<sizeof(chunk_sizes)/sizeof(chunk_sizes[0]);
+          ++i)
+       {
+         Vector<double> v (100);
+         Vector<double> w (100);
+         test (chunk_sizes[i],v,w);
+       }
+    }
+  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_vector_05/cmp/generic b/tests/bits/chunk_sparse_matrix_vector_05/cmp/generic
new file mode 100644 (file)
index 0000000..c9bdf33
--- /dev/null
@@ -0,0 +1,6 @@
+
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
diff --git a/tests/bits/chunk_sparse_matrix_vector_06.cc b/tests/bits/chunk_sparse_matrix_vector_06.cc
new file mode 100644 (file)
index 0000000..493dd45
--- /dev/null
@@ -0,0 +1,107 @@
+//----------------------------  chunk_sparse_matrix_vector_06.cc  ---------------------------
+//    $Id$
+//    Version: $Name$ 
+//
+//    Copyright (C) 2004, 2005, 2008 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_vector_06.cc  ---------------------------
+
+
+// check ChunkSparseMatrix::matrix_norm_square
+
+#include "../tests.h"
+#include <lac/vector.h>
+#include <lac/chunk_sparse_matrix.h>
+#include <fstream>
+#include <iomanip>
+#include <vector>
+
+
+void test (const unsigned int chunk_size,
+          Vector<double> &v)
+{
+                                   // set some entries in the
+                                   // matrix. actually, set them all
+  ChunkSparsityPattern sp (v.size(),v.size(),v.size(), chunk_size);
+  for (unsigned int i=0; i<v.size(); ++i)
+    for (unsigned int j=0; j<v.size(); ++j)
+      sp.add (i,j);
+  sp.compress ();
+
+                                   // then create a matrix from that
+  ChunkSparseMatrix<double> m(sp);
+  for (unsigned int i=0; i<m.m(); ++i)
+    for (unsigned int j=0; j<m.n(); ++j)
+        m.set (i,j, i+2*j);
+
+  for (unsigned int i=0; i<v.size(); ++i)
+    v(i) = i;
+  
+  v.compress ();
+
+                                   // <w,Mv>
+  const double s = m.matrix_norm_square (v);
+
+                                   // make sure we get the expected result
+  for (unsigned int i=0; i<v.size(); ++i)
+    Assert (v(i) == i, ExcInternalError());
+
+  double result = 0;
+  for (unsigned int i=0; i<m.m(); ++i)
+    for (unsigned int j=0; j<m.n(); ++j)
+      result += (i+2*j)*j*i;
+
+  Assert (s == result, ExcInternalError());
+
+  deallog << "OK" << std::endl;
+}
+
+
+
+int main () 
+{
+  std::ofstream logfile("chunk_sparse_matrix_vector_06/output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  try
+    {
+      const unsigned int chunk_sizes[] = { 1, 2, 4, 7, 11 };
+      for (unsigned int i=0; i<sizeof(chunk_sizes)/sizeof(chunk_sizes[0]);
+          ++i)
+       {
+         Vector<double> v (100);
+         test (chunk_sizes[i], v);
+       }
+    }
+  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_vector_06/cmp/generic b/tests/bits/chunk_sparse_matrix_vector_06/cmp/generic
new file mode 100644 (file)
index 0000000..c9bdf33
--- /dev/null
@@ -0,0 +1,6 @@
+
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
diff --git a/tests/bits/chunk_sparse_matrix_vector_07.cc b/tests/bits/chunk_sparse_matrix_vector_07.cc
new file mode 100644 (file)
index 0000000..a6f536d
--- /dev/null
@@ -0,0 +1,119 @@
+//----------------------------  chunk_sparse_matrix_vector_07.cc  ---------------------------
+//    chunk_sparse_matrix_vector_07.cc,v 1.5 2004/02/26 17:25:45 wolf Exp
+//    Version:  
+//
+//    Copyright (C) 2004, 2005, 2008 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_vector_07.cc  ---------------------------
+
+
+// check ChunkSparseMatrix::residual
+
+#include "../tests.h"
+#include <lac/vector.h>
+#include <lac/chunk_sparse_matrix.h>
+#include <fstream>
+#include <iomanip>
+#include <vector>
+
+
+void test (const unsigned int chunk_size,
+          Vector<double> &v,
+           Vector<double> &w,
+           Vector<double> &x)
+{
+                                   // set some entries in the
+                                   // matrix. actually, set them all
+  ChunkSparsityPattern sp (v.size(),v.size(),v.size(), chunk_size);
+  for (unsigned int i=0; i<v.size(); ++i)
+    for (unsigned int j=0; j<v.size(); ++j)
+      sp.add (i,j);
+  sp.compress ();
+
+                                   // then create a matrix from that
+  ChunkSparseMatrix<double> m(sp);
+  for (unsigned int i=0; i<m.m(); ++i)
+    for (unsigned int j=0; j<m.n(); ++j)
+        m.set (i,j, i+2*j);
+
+  for (unsigned int i=0; i<v.size(); ++i)
+    {
+      v(i) = i;
+      w(i) = i+1;
+    }
+      
+  v.compress ();
+  w.compress ();
+
+                                   // x=w-Mv
+  const double s = m.residual (x, v, w);
+
+                                   // make sure we get the expected result
+  for (unsigned int i=0; i<v.size(); ++i)
+    {
+      Assert (v(i) == i, ExcInternalError());
+      Assert (w(i) == i+1, ExcInternalError());
+
+      double result = i+1;
+      for (unsigned int j=0; j<m.n(); ++j)
+        result -= (i+2*j)*j;
+
+      Assert (x(i) == result, ExcInternalError());
+    }
+
+  Assert (std::fabs((s - x.l2_norm())/s) < 1e-14, ExcInternalError());
+
+  deallog << "OK" << std::endl;
+}
+
+
+
+int main () 
+{
+  std::ofstream logfile("chunk_sparse_matrix_vector_07/output");
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  try
+    {
+      const unsigned int chunk_sizes[] = { 1, 2, 4, 7, 11 };
+      for (unsigned int i=0; i<sizeof(chunk_sizes)/sizeof(chunk_sizes[0]);
+          ++i)
+       {
+         Vector<double> v (100);
+         Vector<double> w (100);
+         Vector<double> x (100);
+         test (chunk_sizes[i],v,w,x);
+       }
+    }
+  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_vector_07/cmp/generic b/tests/bits/chunk_sparse_matrix_vector_07/cmp/generic
new file mode 100644 (file)
index 0000000..c9bdf33
--- /dev/null
@@ -0,0 +1,6 @@
+
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
+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.