--- /dev/null
+//---------------------------------------------------------------------------
+// $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 ---------------------------*/
--- /dev/null
+//---------------------------------------------------------------------------
+// $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
--- /dev/null
+//---------------------------------------------------------------------------
+// $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
--- /dev/null
+//---------------------------- 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;
+// }
############################################################
# 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
############################################################
############################################################
find_cell_* \
sparsity_pattern_* \
sparse_matrix_* \
+ chunk_sparse_matrix_* \
full_matrix_vector_* \
solver_* \
deal_solver_* \
--- /dev/null
+//---------------------------- 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;
+ };
+}
--- /dev/null
+
+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
--- /dev/null
+//---------------------------- 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;
+ };
+}
--- /dev/null
+
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
--- /dev/null
+//---------------------------- 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;
+ };
+}
--- /dev/null
+
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
--- /dev/null
+//---------------------------- 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;
+ };
+}
--- /dev/null
+
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
--- /dev/null
+//---------------------------- 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;
+ };
+}
--- /dev/null
+
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
--- /dev/null
+//---------------------------- 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;
+ };
+}
--- /dev/null
+
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
--- /dev/null
+//---------------------------- 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;
+ };
+}
--- /dev/null
+
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
--- /dev/null
+//---------------------------- 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;
+ };
+}
--- /dev/null
+
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
--- /dev/null
+//---------------------------- 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;
+ };
+}
--- /dev/null
+
+DEAL::13
+DEAL::OK
+DEAL::23
+DEAL::OK
+DEAL::25
+DEAL::OK
+DEAL::25
+DEAL::OK
+DEAL::25
+DEAL::OK
--- /dev/null
+//---------------------------- 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;
+ };
+}
--- /dev/null
+
+DEAL::8
+DEAL::OK
+DEAL::22
+DEAL::OK
+DEAL::24
+DEAL::OK
+DEAL::25
+DEAL::OK
+DEAL::25
+DEAL::OK
--- /dev/null
+//---------------------------- 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;
+ };
+}
--- /dev/null
+
+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
--- /dev/null
+//---------------------------- 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;
+ };
+}
--- /dev/null
+
+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
--- /dev/null
+//---------------------------- 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;
+ };
+}
--- /dev/null
+
+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
--- /dev/null
+//---------------------------- 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;
+ };
+}
--- /dev/null
+
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
--- /dev/null
+//---------------------------- 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;
+ };
+}
--- /dev/null
+
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
--- /dev/null
+//---------------------------- 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;
+ };
+}
--- /dev/null
+
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
--- /dev/null
+//---------------------------- 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;
+ };
+}
--- /dev/null
+
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
--- /dev/null
+//---------------------------- 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;
+ };
+}
--- /dev/null
+
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
--- /dev/null
+//---------------------------- 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;
+ };
+}
--- /dev/null
+
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
--- /dev/null
+//---------------------------- 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;
+ };
+}
--- /dev/null
+
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
--- /dev/null
+//---------------------------- 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;
+ };
+}
--- /dev/null
+
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
--- /dev/null
+//---------------------------- 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;
+ };
+}
--- /dev/null
+
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK