*@{
*/
+/**
+ * A namespace in which we declare iterators over the elements of sparse
+ * matrices.
+ */
+namespace ChunkSparseMatrixIterators
+{
+ // forward declaration
+ template <typename number, bool Constness>
+ class Iterator;
+
+ /**
+ * General template for sparse matrix accessors. The first template argument
+ * denotes the underlying numeric type, the second the constness of the
+ * matrix.
+ *
+ * The general template is not implemented, only the specializations for the
+ * two possible values of the second template argument. Therefore, the
+ * interface listed here only serves as a template provided since doxygen
+ * does not link the specializations.
+ */
+ template <typename number, bool Constness>
+ class Accessor : public ChunkSparsityPatternIterators::Accessor
+ {
+ public:
+ /**
+ * Value of this matrix entry.
+ */
+ number value() const;
+
+ /**
+ * Value of this matrix entry.
+ */
+ number &value();
+
+ /**
+ * Return a reference to the matrix into which this accessor points. Note
+ * that in the present case, this is a constant reference.
+ */
+ const ChunkSparseMatrix<number> &get_matrix () const;
+ };
+
+
+
+ /**
+ * Accessor class for constant matrices, used in the const_iterators. This
+ * class builds on the accessor classes used for sparsity patterns to loop
+ * over all nonzero entries, and only adds the accessor functions to gain
+ * access to the actual value stored at a certain location.
+ */
+ template <typename number>
+ class Accessor<number,true> : public ChunkSparsityPatternIterators::Accessor
+ {
+ public:
+ /**
+ * Typedef for the type (including constness) of the matrix to be used
+ * here.
+ */
+ typedef const ChunkSparseMatrix<number> MatrixType;
+
+ /**
+ * Constructor.
+ */
+ Accessor (MatrixType *matrix,
+ const unsigned int row);
+
+ /**
+ * Constructor. Construct the end accessor for the given matrix.
+ */
+ Accessor (MatrixType *matrix);
+
+ /**
+ * Copy constructor to get from a non-const accessor to a const accessor.
+ */
+ Accessor (const ChunkSparseMatrixIterators::Accessor<number,false> &a);
+
+ /**
+ * Value of this matrix entry.
+ */
+ number value() const;
+
+ /**
+ * Return a reference to the matrix into which this accessor points. Note
+ * that in the present case, this is a constant reference.
+ */
+ MatrixType &get_matrix () const;
+
+ private:
+ /**
+ * Pointer to the matrix we use.
+ */
+ MatrixType *matrix;
+
+ /**
+ * Make the advance function of the base class available.
+ */
+ using ChunkSparsityPatternIterators::Accessor::advance;
+
+ /**
+ * Make iterator class a friend.
+ */
+ template <typename, bool>
+ friend class Iterator;
+ };
+
+
+ /**
+ * Accessor class for non-constant matrices, used in the iterators. This
+ * class builds on the accessor classes used for sparsity patterns to loop
+ * over all nonzero entries, and only adds the accessor functions to gain
+ * access to the actual value stored at a certain location.
+ */
+ template <typename number>
+ class Accessor<number,false> : public ChunkSparsityPatternIterators::Accessor
+ {
+ private:
+ /**
+ * Reference class. This is what the accessor class returns when you call
+ * the value() function. The reference acts just as if it were a reference
+ * to the actual value of a matrix entry, i.e. you can read and write it,
+ * you can add and multiply to it, etc, but since the matrix does not give
+ * away the address of this matrix entry, we have to go through functions
+ * to do all this.
+ *
+ * The constructor takes a pointer to an accessor object that describes
+ * which element of the matrix it points to. This creates an ambiguity
+ * when one writes code like iterator->value()=0 (instead of
+ * iterator->value()=0.0), since the right hand side is an integer that
+ * can both be converted to a <tt>number</tt> (i.e., most commonly a
+ * double) or to another object of type <tt>Reference</tt>. The compiler
+ * then complains about not knowing which conversion to take.
+ *
+ * For some reason, adding another overload operator=(int) doesn't seem to
+ * cure the problem. We avoid it, however, by adding a second, dummy
+ * argument to the Reference constructor, that is unused, but makes sure
+ * there is no second matching conversion sequence using a one-argument
+ * right hand side.
+ */
+ class Reference
+ {
+ public:
+ /**
+ * Constructor. For the second argument, see the general class
+ * documentation.
+ */
+ Reference (const Accessor *accessor,
+ const bool dummy);
+
+ /**
+ * Conversion operator to the data type of the matrix.
+ */
+ operator number () const;
+
+ /**
+ * Set the element of the matrix we presently point to to @p n.
+ */
+ const Reference &operator = (const number n) const;
+
+ /**
+ * Add @p n to the element of the matrix we presently point to.
+ */
+ const Reference &operator += (const number n) const;
+
+ /**
+ * Subtract @p n from the element of the matrix we presently point to.
+ */
+ const Reference &operator -= (const number n) const;
+
+ /**
+ * Multiply the element of the matrix we presently point to by @p n.
+ */
+ const Reference &operator *= (const number n) const;
+
+ /**
+ * Divide the element of the matrix we presently point to by @p n.
+ */
+ const Reference &operator /= (const number n) const;
+
+ private:
+ /**
+ * Pointer to the accessor that denotes which element we presently point
+ * to.
+ */
+ const Accessor *accessor;
+ };
+
+ public:
+ /**
+ * Typedef for the type (including constness) of the matrix to be used
+ * here.
+ */
+ typedef ChunkSparseMatrix<number> MatrixType;
+
+ /**
+ * Constructor.
+ */
+ Accessor (MatrixType *matrix,
+ const unsigned int row);
+
+ /**
+ * Constructor. Construct the end accessor for the given matrix.
+ */
+ Accessor (MatrixType *matrix);
+
+ /**
+ * Value of this matrix entry, returned as a read- and writable reference.
+ */
+ Reference value() const;
+
+ /**
+ * Return a reference to the matrix into which this accessor points. Note
+ * that in the present case, this is a non-constant reference.
+ */
+ MatrixType &get_matrix () const;
+
+ private:
+ /**
+ * Pointer to the matrix we use.
+ */
+ MatrixType *matrix;
+
+ /**
+ * Make the advance function of the base class available.
+ */
+ using ChunkSparsityPatternIterators::Accessor::advance;
+
+ /**
+ * Make iterator class a friend.
+ */
+ template <typename, bool>
+ friend class Iterator;
+
+ /**
+ * Make the inner reference class a friend if the compiler has a bug and
+ * requires this.
+ */
+ };
+
+
+
+ /**
+ * STL conforming iterator for constant and non-constant matrices.
+ *
+ * The first template argument denotes the underlying numeric type, the
+ * second the constness of the matrix.
+ *
+ * Since there is a specialization of this class for
+ * <tt>Constness=false</tt>, this class is for iterators to constant
+ * matrices.
+ */
+ template <typename number, bool Constness>
+ class Iterator
+ {
+ public:
+ /**
+ * Typedef for the matrix type (including constness) we are to operate on.
+ */
+ typedef
+ typename Accessor<number,Constness>::MatrixType
+ MatrixType;
+
+ /**
+ * A typedef for the type you get when you dereference an iterator
+ * of the current kind.
+ */
+ typedef
+ const Accessor<number,Constness> & value_type;
+
+ /**
+ * Constructor. Create an iterator into the matrix @p matrix for the given
+ * row and the index within it.
+ */
+ Iterator (MatrixType *matrix,
+ const unsigned int row);
+
+ /**
+ * Constructor. Create the end iterator for the given matrix.
+ */
+ Iterator (MatrixType *matrix);
+
+ /**
+ * Conversion constructor to get from a non-const iterator to a const
+ * iterator.
+ */
+ Iterator (const ChunkSparseMatrixIterators::Iterator<number,false> &i);
+
+ /**
+ * Prefix increment.
+ */
+ Iterator &operator++ ();
+
+ /**
+ * Postfix increment.
+ */
+ Iterator operator++ (int);
+
+ /**
+ * Dereferencing operator.
+ */
+ const Accessor<number,Constness> &operator* () const;
+
+ /**
+ * Dereferencing operator.
+ */
+ const Accessor<number,Constness> *operator-> () const;
+
+ /**
+ * Comparison. True, if both iterators point to the same matrix position.
+ */
+ bool operator == (const Iterator &) const;
+
+ /**
+ * Inverse of <tt>==</tt>.
+ */
+ bool operator != (const Iterator &) const;
+
+ /**
+ * Comparison operator. Result is true if either the first row number is
+ * smaller or if the row numbers are equal and the first index is smaller.
+ *
+ * This function is only valid if both iterators point into the same
+ * matrix.
+ */
+ bool operator < (const Iterator &) const;
+
+ /**
+ * Comparison operator. Works in the same way as above operator, just the
+ * other way round.
+ */
+ bool operator > (const Iterator &) const;
+
+ /**
+ * Return the distance between the current iterator and the argument.
+ * The distance is given by how many times one has to apply operator++
+ * to the current iterator to get the argument (for a positive return
+ * value), or operator-- (for a negative return value).
+ */
+ int operator - (const Iterator &p) const;
+
+ /**
+ * Return an iterator that is @p n ahead of the current one.
+ */
+ Iterator operator + (const unsigned int n) const;
+
+ private:
+ /**
+ * Store an object of the accessor class.
+ */
+ Accessor<number,Constness> accessor;
+ };
+
+}
+
+
/**
* Sparse matrix. This class implements the function to store values
typedef number value_type;
/**
- * Declare a type that has holds
- * real-valued numbers with the
- * same precision as the template
- * argument to this class. If the
- * template argument of this
- * class is a real data type,
- * then real_type equals the
- * template argument. If the
- * template argument is a
- * std::complex type then
- * real_type equals the type
- * underlying the complex
- * numbers.
+ * Declare a type that has holds real-valued numbers with the same precision
+ * as the template argument to this class. If the template argument of this
+ * class is a real data type, then real_type equals the template
+ * argument. If the template argument is a std::complex type then real_type
+ * equals the type underlying the complex numbers.
*
- * This typedef is used to
- * represent the return type of
- * norms.
+ * This typedef is used to represent the return type of norms.
*/
typedef typename numbers::NumberTraits<number>::real_type real_type;
/**
- * A structure that describes some of the
- * traits of this class in terms of its
- * run-time behavior. Some other classes
- * (such as the block matrix classes)
- * that take one or other of the matrix
- * classes as its template parameters can
- * tune their behavior based on the
- * variables in this class.
+ * Typedef of an STL conforming iterator class walking over all the nonzero
+ * entries of this matrix. This iterator cannot change the values of the
+ * matrix.
+ */
+ typedef
+ ChunkSparseMatrixIterators::Iterator<number,true>
+ const_iterator;
+
+ /**
+ * Typedef of an STL conforming iterator class walking over all the nonzero
+ * entries of this matrix. This iterator @em can change the values of the
+ * matrix, but of course can't change the sparsity pattern as this is fixed
+ * once a sparse matrix is attached to it.
+ */
+ typedef
+ ChunkSparseMatrixIterators::Iterator<number,false>
+ iterator;
+
+ /**
+ * A structure that describes some of the traits of this class in terms of
+ * its run-time behavior. Some other classes (such as the block matrix
+ * classes) that take one or other of the matrix classes as its template
+ * parameters can tune their behavior based on the variables in this class.
*/
struct Traits
{
/**
- * It is safe to elide additions of
- * zeros to individual elements of
- * this matrix.
+ * It is safe to elide additions of zeros to individual elements of this
+ * matrix.
*/
static const bool zero_addition_can_be_elided = true;
};
*/
//@{
/**
- * Constructor; initializes the matrix to
- * be empty, without any structure, i.e.
- * the matrix is not usable at all. This
- * constructor is therefore only useful
- * for matrices which are members of a
- * class. All other matrices should be
- * created at a point in the data flow
- * where all necessary information is
- * available.
+ * Constructor; initializes the matrix to be empty, without any structure,
+ * i.e. the matrix is not usable at all. This constructor is therefore only
+ * useful for matrices which are members of a class. All other matrices
+ * should be created at a point in the data flow where all necessary
+ * information is available.
*
- * You have to initialize
- * the matrix before usage with
- * reinit(const ChunkSparsityPattern&).
+ * You have to initialize the matrix before usage with reinit(const
+ * ChunkSparsityPattern&).
*/
ChunkSparseMatrix ();
/**
- * Copy constructor. This constructor is
- * only allowed to be called if the matrix
- * to be copied is empty. This is for the
- * same reason as for the
- * ChunkSparsityPattern, see there for the
- * details.
+ * Copy constructor. This constructor is only allowed to be called if the
+ * matrix to be copied is empty. This is for the same reason as for the
+ * ChunkSparsityPattern, see there for the details.
*
- * If you really want to copy a whole
- * matrix, you can do so by using the
+ * If you really want to copy a whole matrix, you can do so by using the
* copy_from() function.
*/
ChunkSparseMatrix (const ChunkSparseMatrix &);
/**
- * Constructor. Takes the given
- * matrix sparsity structure to
- * represent the sparsity pattern
- * of this matrix. You can change
- * the sparsity pattern later on
- * by calling the reinit(const
- * ChunkSparsityPattern&) function.
+ * Constructor. Takes the given matrix sparsity structure to represent the
+ * sparsity pattern of this matrix. You can change the sparsity pattern
+ * later on by calling the reinit(const ChunkSparsityPattern&) function.
*
- * You have to make sure that the
- * lifetime of the sparsity
- * structure is at least as long
- * as that of this matrix or as
- * long as reinit(const
- * ChunkSparsityPattern&) is not
- * called with a new sparsity
- * pattern.
+ * You have to make sure that the lifetime of the sparsity structure is at
+ * least as long as that of this matrix or as long as reinit(const
+ * ChunkSparsityPattern&) is not called with a new sparsity pattern.
*
- * The constructor is marked
- * explicit so as to disallow
- * that someone passes a sparsity
- * pattern in place of a sparse
- * matrix to some function, where
- * an empty matrix would be
- * generated then.
+ * The constructor is marked explicit so as to disallow that someone passes
+ * a sparsity pattern in place of a sparse matrix to some function, where an
+ * empty matrix would be generated then.
*/
explicit ChunkSparseMatrix (const ChunkSparsityPattern &sparsity);
/**
- * Copy constructor: initialize
- * the matrix with the identity
- * matrix. This constructor will
- * throw an exception if the
- * sizes of the sparsity pattern
- * and the identity matrix do not
- * coincide, or if the sparsity
- * pattern does not provide for
- * nonzero entries on the entire
- * diagonal.
+ * Copy constructor: initialize the matrix with the identity matrix. This
+ * constructor will throw an exception if the sizes of the sparsity pattern
+ * and the identity matrix do not coincide, or if the sparsity pattern does
+ * not provide for nonzero entries on the entire diagonal.
*/
ChunkSparseMatrix (const ChunkSparsityPattern &sparsity,
const IdentityMatrix &id);
/**
- * Destructor. Free all memory, but do not
- * release the memory of the sparsity
- * structure.
+ * Destructor. Free all memory, but do not release the memory of the
+ * sparsity structure.
*/
virtual ~ChunkSparseMatrix ();
/**
- * Copy operator. Since copying
- * entire sparse matrices is a
- * very expensive operation, we
- * disallow doing so except for
- * the special case of empty
- * matrices of size zero. This
- * doesn't seem particularly
- * useful, but is exactly what
- * one needs if one wanted to
- * have a
- * <code>std::vector@<ChunkSparseMatrix@<double@>
- * @></code>: in that case, one
- * can create a vector (which
- * needs the ability to copy
- * objects) of empty matrices
- * that are then later filled
- * with something useful.
+ * Copy operator. Since copying entire sparse matrices is a very expensive
+ * operation, we disallow doing so except for the special case of empty
+ * matrices of size zero. This doesn't seem particularly useful, but is
+ * exactly what one needs if one wanted to have a
+ * <code>std::vector@<ChunkSparseMatrix@<double@> @></code>: in that case,
+ * one can create a vector (which needs the ability to copy objects) of
+ * empty matrices that are then later filled with something useful.
*/
ChunkSparseMatrix<number> &operator = (const ChunkSparseMatrix<number> &);
/**
- * Copy operator: initialize
- * the matrix with the identity
- * matrix. This operator will
- * throw an exception if the
- * sizes of the sparsity pattern
- * and the identity matrix do not
- * coincide, or if the sparsity
- * pattern does not provide for
- * nonzero entries on the entire
- * diagonal.
+ * Copy operator: initialize the matrix with the identity matrix. This
+ * operator will throw an exception if the sizes of the sparsity pattern and
+ * the identity matrix do not coincide, or if the sparsity pattern does not
+ * provide for nonzero entries on the entire diagonal.
*/
ChunkSparseMatrix<number> &
operator= (const IdentityMatrix &id);
/**
- * This operator assigns a scalar to
- * a matrix. Since this does usually
- * not make much sense (should we set
- * all matrix entries to this value?
- * Only the nonzero entries of the
- * sparsity pattern?), this operation
- * is only allowed if the actual
- * value to be assigned is zero. This
- * operator only exists to allow for
- * the obvious notation
- * <tt>matrix=0</tt>, which sets all
- * elements of the matrix to zero,
- * but keep the sparsity pattern
+ * This operator assigns a scalar to a matrix. Since this does usually not
+ * make much sense (should we set all matrix entries to this value? Only
+ * the nonzero entries of the sparsity pattern?), this operation is only
+ * allowed if the actual value to be assigned is zero. This operator only
+ * exists to allow for the obvious notation <tt>matrix=0</tt>, which sets
+ * all elements of the matrix to zero, but keep the sparsity pattern
* previously used.
*/
ChunkSparseMatrix &operator = (const double d);
/**
- * Reinitialize the sparse matrix
- * with the given sparsity
- * pattern. The latter tells the
- * matrix how many nonzero
- * elements there need to be
+ * Reinitialize the sparse matrix with the given sparsity pattern. The
+ * latter tells the matrix how many nonzero elements there need to be
* reserved.
*
- * Regarding memory allocation,
- * the same applies as said
- * above.
+ * Regarding memory allocation, the same applies as said above.
*
- * You have to make sure that the
- * lifetime of the sparsity
- * structure is at least as long
- * as that of this matrix or as
- * long as reinit(const
- * ChunkSparsityPattern &) is not
- * called with a new sparsity
- * structure.
+ * You have to make sure that the lifetime of the sparsity structure is at
+ * least as long as that of this matrix or as long as reinit(const
+ * ChunkSparsityPattern &) is not called with a new sparsity structure.
*
- * The elements of the matrix are
- * set to zero by this function.
+ * The elements of the matrix are set to zero by this function.
*/
virtual void reinit (const ChunkSparsityPattern &sparsity);
/**
- * Release all memory and return
- * to a state just like after
- * having called the default
- * constructor. It also forgets
- * the sparsity pattern it was
+ * Release all memory and return to a state just like after having called
+ * the default constructor. It also forgets the sparsity pattern it was
* previously tied to.
*/
virtual void clear ();
*/
//@{
/**
- * Return whether the object is
- * empty. It is empty if either
- * both dimensions are zero or no
- * ChunkSparsityPattern is
- * associated.
+ * Return whether the object is empty. It is empty if either both dimensions
+ * are zero or no ChunkSparsityPattern is associated.
*/
bool empty () const;
/**
- * Return the dimension of the
- * image space. To remember: the
- * matrix is of dimension
- * $m \times n$.
+ * Return the dimension of the image space. To remember: the matrix is of
+ * dimension $m \times n$.
*/
unsigned int m () const;
/**
- * Return the dimension of the
- * range space. To remember: the
- * matrix is of dimension
- * $m \times n$.
+ * Return the dimension of the range space. To remember: the matrix is of
+ * dimension $m \times n$.
*/
unsigned int n () const;
/**
- * Return the number of nonzero
- * elements of this
- * matrix. Actually, it returns
- * the number of entries in the
- * sparsity pattern; if any of
- * the entries should happen to
- * be zero, it is counted anyway.
+ * Return the number of nonzero elements of this matrix. Actually, it
+ * returns the number of entries in the sparsity pattern; if any of the
+ * entries should happen to be zero, it is counted anyway.
*/
unsigned int n_nonzero_elements () const;
/**
- * Return the number of actually
- * nonzero elements of this
- * matrix.
+ * Return the number of actually nonzero elements of this matrix.
*
- * Note, that this function does
- * (in contrary to
- * n_nonzero_elements()) not
- * count all entries of the
- * sparsity pattern but only the
- * ones that are nonzero.
+ * Note, that this function does (in contrary to n_nonzero_elements()) not
+ * count all entries of the sparsity pattern but only the ones that are
+ * nonzero.
*/
unsigned int n_actually_nonzero_elements () const;
/**
- * Return a (constant) reference
- * to the underlying sparsity
- * pattern of this matrix.
+ * Return a (constant) reference to the underlying sparsity pattern of this
+ * matrix.
*
- * Though the return value is
- * declared <tt>const</tt>, you
- * should be aware that it may
- * change if you call any
- * nonconstant function of
- * objects which operate on it.
+ * Though the return value is declared <tt>const</tt>, you should be aware
+ * that it may change if you call any nonconstant function of objects which
+ * operate on it.
*/
const ChunkSparsityPattern &get_sparsity_pattern () const;
/**
- * Determine an estimate for the
- * memory consumption (in bytes)
- * of this object. See
- * MemoryConsumption.
+ * Determine an estimate for the memory consumption (in bytes) of this
+ * object. See MemoryConsumption.
*/
std::size_t memory_consumption () const;
*/
//@{
/**
- * Set the element (<i>i,j</i>)
- * to <tt>value</tt>. Throws an
- * error if the entry does not
- * exist or if <tt>value</tt> is
- * not a finite number. Still, it
- * is allowed to store zero
- * values in non-existent fields.
+ * Set the element (<i>i,j</i>) to <tt>value</tt>. Throws an error if the
+ * entry does not exist or if <tt>value</tt> is not a finite number. Still,
+ * it is allowed to store zero values in non-existent fields.
*/
void set (const unsigned int i,
const unsigned int j,
const number value);
/**
- * Add <tt>value</tt> to the
- * element (<i>i,j</i>). Throws
- * an error if the entry does not
- * exist or if <tt>value</tt> is
- * not a finite number. Still, it
- * is allowed to store zero
- * values in non-existent fields.
+ * Add <tt>value</tt> to the element (<i>i,j</i>). Throws an error if the
+ * entry does not exist or if <tt>value</tt> is not a finite number. Still,
+ * it is allowed to store zero values in non-existent fields.
*/
void add (const unsigned int i,
const unsigned int j,
const number value);
/**
- * Add an array of values given by
- * <tt>values</tt> in the given
- * global matrix row at columns
- * specified by col_indices in the
- * sparse matrix.
+ * Add an array of values given by <tt>values</tt> in the given global
+ * matrix row at columns specified by col_indices in the sparse matrix.
*
- * The optional parameter
- * <tt>elide_zero_values</tt> can be
- * used to specify whether zero
- * values should be added anyway or
- * these should be filtered away and
- * only non-zero data is added. The
- * default value is <tt>true</tt>,
- * i.e., zero values won't be added
- * into the matrix.
+ * The optional parameter <tt>elide_zero_values</tt> can be used to specify
+ * whether zero values should be added anyway or these should be filtered
+ * away and only non-zero data is added. The default value is <tt>true</tt>,
+ * i.e., zero values won't be added into the matrix.
*/
template <typename number2>
void add (const unsigned int row,
const bool col_indices_are_sorted = false);
/**
- * Multiply the entire matrix by a
- * fixed factor.
+ * Multiply the entire matrix by a fixed factor.
*/
ChunkSparseMatrix &operator *= (const number factor);
/**
- * Divide the entire matrix by a
- * fixed factor.
+ * Divide the entire matrix by a fixed factor.
*/
ChunkSparseMatrix &operator /= (const number factor);
/**
- * Symmetrize the matrix by
- * forming the mean value between
- * the existing matrix and its
- * transpose, $A = \frac 12(A+A^T)$.
+ * Symmetrize the matrix by forming the mean value between the existing
+ * matrix and its transpose, $A = \frac 12(A+A^T)$.
*
- * This operation assumes that
- * the underlying sparsity
- * pattern represents a symmetric
- * object. If this is not the
- * case, then the result of this
- * operation will not be a
- * symmetric matrix, since it
- * only explicitly symmetrizes
- * by looping over the lower left
- * triangular part for efficiency
- * reasons; if there are entries
- * in the upper right triangle,
- * then these elements are missed
- * in the
- * symmetrization. Symmetrization
- * of the sparsity pattern can be
- * obtain by
- * ChunkSparsityPattern::symmetrize().
+ * This operation assumes that the underlying sparsity pattern represents a
+ * symmetric object. If this is not the case, then the result of this
+ * operation will not be a symmetric matrix, since it only explicitly
+ * symmetrizes by looping over the lower left triangular part for efficiency
+ * reasons; if there are entries in the upper right triangle, then these
+ * elements are missed in the symmetrization. Symmetrization of the sparsity
+ * pattern can be obtain by ChunkSparsityPattern::symmetrize().
*/
void symmetrize ();
/**
- * Copy the given matrix to this
- * one. The operation throws an
- * error if the sparsity patterns
- * of the two involved matrices
- * do not point to the same
- * object, since in this case the
- * copy operation is
- * cheaper. Since this operation
- * is notheless not for free, we
- * do not make it available
- * through <tt>operator =</tt>,
- * since this may lead to
- * unwanted usage, e.g. in copy
- * arguments to functions, which
- * should really be arguments by
- * reference.
+ * Copy the given matrix to this one. The operation throws an error if the
+ * sparsity patterns of the two involved matrices do not point to the same
+ * object, since in this case the copy operation is cheaper. Since this
+ * operation is notheless not for free, we do not make it available through
+ * <tt>operator =</tt>, since this may lead to unwanted usage, e.g. in copy
+ * arguments to functions, which should really be arguments by reference.
*
- * The source matrix may be a matrix
- * of arbitrary type, as long as its
- * data type is convertible to the
- * data type of this matrix.
+ * The source matrix may be a matrix of arbitrary type, as long as its data
+ * type is convertible to the data type of this matrix.
*
- * The function returns a reference to
- * <tt>*this</tt>.
+ * The function returns a reference to <tt>*this</tt>.
*/
template <typename somenumber>
ChunkSparseMatrix<number> &
copy_from (const ChunkSparseMatrix<somenumber> &source);
/**
- * This function is complete
- * analogous to the
- * ChunkSparsityPattern::copy_from()
- * function in that it allows to
- * initialize a whole matrix in
- * one step. See there for more
- * information on argument types
- * and their meaning. You can
- * also find a small example on
- * how to use this function
- * there.
+ * This function is complete analogous to the
+ * ChunkSparsityPattern::copy_from() function in that it allows to
+ * initialize a whole matrix in one step. See there for more information on
+ * argument types and their meaning. You can also find a small example on
+ * how to use this function there.
*
- * The only difference to the
- * cited function is that the
- * objects which the inner
- * iterator points to need to be
- * of type <tt>std::pair<unsigned
- * int, value</tt>, where
- * <tt>value</tt> needs to be
- * convertible to the element
- * type of this class, as
- * specified by the
- * <tt>number</tt> template
+ * The only difference to the cited function is that the objects which the
+ * inner iterator points to need to be of type <tt>std::pair<unsigned int,
+ * value</tt>, where <tt>value</tt> needs to be convertible to the element
+ * type of this class, as specified by the <tt>number</tt> template
* argument.
*
- * Previous content of the matrix
- * is overwritten. Note that the
- * entries specified by the input
- * parameters need not
- * necessarily cover all elements
- * of the matrix. Elements not
- * covered remain untouched.
+ * Previous content of the matrix is overwritten. Note that the entries
+ * specified by the input parameters need not necessarily cover all elements
+ * of the matrix. Elements not covered remain untouched.
*/
template <typename ForwardIterator>
void copy_from (const ForwardIterator begin,
const ForwardIterator end);
/**
- * Copy the nonzero entries of a
- * full matrix into this
- * object. Previous content is
- * deleted. Note that the
- * underlying sparsity pattern
- * must be appropriate to hold
- * the nonzero entries of the
- * full matrix.
+ * Copy the nonzero entries of a full matrix into this object. Previous
+ * content is deleted. Note that the underlying sparsity pattern must be
+ * appropriate to hold the nonzero entries of the full matrix.
*/
template <typename somenumber>
void copy_from (const FullMatrix<somenumber> &matrix);
/**
- * Add <tt>matrix</tt> scaled by
- * <tt>factor</tt> to this matrix,
- * i.e. the matrix <tt>factor*matrix</tt>
- * is added to <tt>this</tt>. This
- * function throws an error if the
- * sparsity patterns of the two involved
- * matrices do not point to the same
- * object, since in this case the
- * operation is cheaper.
+ * Add <tt>matrix</tt> scaled by <tt>factor</tt> to this matrix, i.e. the
+ * matrix <tt>factor*matrix</tt> is added to <tt>this</tt>. This function
+ * throws an error if the sparsity patterns of the two involved matrices do
+ * not point to the same object, since in this case the operation is
+ * cheaper.
*
- * The source matrix may be a sparse
- * matrix over an arbitrary underlying
- * scalar type, as long as its data type
- * is convertible to the data type of
+ * The source matrix may be a sparse matrix over an arbitrary underlying
+ * scalar type, as long as its data type is convertible to the data type of
* this matrix.
*/
template <typename somenumber>
//@{
/**
- * Return the value of the entry
- * (<i>i,j</i>). This may be an
- * expensive operation and you
- * should always take care where
- * to call this function. In
- * order to avoid abuse, this
- * function throws an exception
- * if the required element does
- * not exist in the matrix.
+ * Return the value of the entry (<i>i,j</i>). This may be an expensive
+ * operation and you should always take care where to call this function.
+ * In order to avoid abuse, this function throws an exception if the
+ * required element does not exist in the matrix.
*
- * In case you want a function
- * that returns zero instead (for
- * entries that are not in the
- * sparsity pattern of the
- * matrix), use the el()
- * function.
+ * In case you want a function that returns zero instead (for entries that
+ * are not in the sparsity pattern of the matrix), use the el() function.
*
- * If you are looping over all elements,
- * consider using one of the iterator
- * classes instead, since they are
- * tailored better to a sparse matrix
+ * If you are looping over all elements, consider using one of the iterator
+ * classes instead, since they are tailored better to a sparse matrix
* structure.
*/
number operator () (const unsigned int i,
const unsigned int j) const;
/**
- * This function is mostly like
- * operator()() in that it
- * returns the value of the
- * matrix entry (<i>i,j</i>). The
- * only difference is that if
- * this entry does not exist in
- * the sparsity pattern, then
- * instead of raising an
- * exception, zero is
- * returned. While this may be
- * convenient in some cases, note
- * that it is simple to write
- * algorithms that are slow
- * compared to an optimal
- * solution, since the sparsity
- * of the matrix is not used.
+ * This function is mostly like operator()() in that it returns the value of
+ * the matrix entry (<i>i,j</i>). The only difference is that if this entry
+ * does not exist in the sparsity pattern, then instead of raising an
+ * exception, zero is returned. While this may be convenient in some cases,
+ * note that it is simple to write algorithms that are slow compared to an
+ * optimal solution, since the sparsity of the matrix is not used.
*
- * If you are looping over all elements,
- * consider using one of the iterator
- * classes instead, since they are
- * tailored better to a sparse matrix
+ * If you are looping over all elements, consider using one of the iterator
+ * classes instead, since they are tailored better to a sparse matrix
* structure.
*/
number el (const unsigned int i,
* error if the matrix is not
* quadratic.
*
- * This function is considerably
- * faster than the operator()(),
- * since for quadratic matrices, the
- * diagonal entry may be the
- * first to be stored in each row
- * and access therefore does not
- * involve searching for the
- * right column number.
+ * This function is considerably faster than the operator()(), since for
+ * quadratic matrices, the diagonal entry may be the first to be stored in
+ * each row and access therefore does not involve searching for the right
+ * column number.
*/
number diag_element (const unsigned int i) const;
/**
- * Same as above, but return a
- * writeable reference. You're
- * sure you know what you do?
+ * Same as above, but return a writeable reference. You're sure you know
+ * what you do?
*/
number &diag_element (const unsigned int i);
*/
//@{
/**
- * Matrix-vector multiplication:
- * let <i>dst = M*src</i> with
- * <i>M</i> being this matrix.
+ * Matrix-vector multiplication: let <i>dst = M*src</i> with <i>M</i> being
+ * this matrix.
*
- * Note that while this function can
- * operate on all vectors that offer
- * iterator classes, it is only really
- * effective for objects of type @ref
- * Vector. For all classes for which
- * iterating over elements, or random
- * member access is expensive, this
- * function is not efficient. In
- * particular, if you want to multiply
- * with BlockVector objects, you should
- * consider using a BlockChunkSparseMatrix as
- * well.
+ * Note that while this function can operate on all vectors that offer
+ * iterator classes, it is only really effective for objects of type @ref
+ * Vector. For all classes for which iterating over elements, or random
+ * member access is expensive, this function is not efficient. In
+ * particular, if you want to multiply with BlockVector objects, you should
+ * consider using a BlockChunkSparseMatrix as well.
*
- * Source and destination must
- * not be the same vector.
+ * Source and destination must not be the same vector.
*/
template <class OutVector, class InVector>
void vmult (OutVector &dst,
const InVector &src) const;
/**
- * Matrix-vector multiplication:
- * let <i>dst = M<sup>T</sup>*src</i> with
- * <i>M</i> being this
- * matrix. This function does the
- * same as vmult() but takes
- * the transposed matrix.
+ * Matrix-vector multiplication: let <i>dst = M<sup>T</sup>*src</i> with
+ * <i>M</i> being this matrix. This function does the same as vmult() but
+ * takes the transposed matrix.
*
- * Note that while this function can
- * operate on all vectors that offer
- * iterator classes, it is only really
- * effective for objects of type @ref
- * Vector. For all classes for which
- * iterating over elements, or random
- * member access is expensive, this
- * function is not efficient. In
- * particular, if you want to multiply
- * with BlockVector objects, you should
- * consider using a BlockChunkSparseMatrix as
- * well.
+ * Note that while this function can operate on all vectors that offer
+ * iterator classes, it is only really effective for objects of type @ref
+ * Vector. For all classes for which iterating over elements, or random
+ * member access is expensive, this function is not efficient. In
+ * particular, if you want to multiply with BlockVector objects, you should
+ * consider using a BlockChunkSparseMatrix as well.
*
- * Source and destination must
- * not be the same vector.
+ * Source and destination must not be the same vector.
*/
template <class OutVector, class InVector>
void Tvmult (OutVector &dst,
const InVector &src) const;
/**
- * Adding Matrix-vector
- * multiplication. Add
- * <i>M*src</i> on <i>dst</i>
- * with <i>M</i> being this
- * matrix.
+ * Adding Matrix-vector multiplication. Add <i>M*src</i> on <i>dst</i> with
+ * <i>M</i> being this matrix.
*
- * Note that while this function can
- * operate on all vectors that offer
- * iterator classes, it is only really
- * effective for objects of type @ref
- * Vector. For all classes for which
- * iterating over elements, or random
- * member access is expensive, this
- * function is not efficient. In
- * particular, if you want to multiply
- * with BlockVector objects, you should
- * consider using a BlockChunkSparseMatrix as
- * well.
+ * Note that while this function can operate on all vectors that offer
+ * iterator classes, it is only really effective for objects of type @ref
+ * Vector. For all classes for which iterating over elements, or random
+ * member access is expensive, this function is not efficient. In
+ * particular, if you want to multiply with BlockVector objects, you should
+ * consider using a BlockChunkSparseMatrix as well.
*
- * Source and destination must
- * not be the same vector.
+ * Source and destination must not be the same vector.
*/
template <class OutVector, class InVector>
void vmult_add (OutVector &dst,
const InVector &src) const;
/**
- * Adding Matrix-vector
- * multiplication. Add
- * <i>M<sup>T</sup>*src</i> to
- * <i>dst</i> with <i>M</i> being
- * this matrix. This function
- * does the same as vmult_add()
- * but takes the transposed
- * matrix.
+ * Adding Matrix-vector multiplication. Add <i>M<sup>T</sup>*src</i> to
+ * <i>dst</i> with <i>M</i> being this matrix. This function does the same
+ * as vmult_add() but takes the transposed matrix.
*
- * Note that while this function can
- * operate on all vectors that offer
- * iterator classes, it is only really
- * effective for objects of type @ref
- * Vector. For all classes for which
- * iterating over elements, or random
- * member access is expensive, this
- * function is not efficient. In
- * particular, if you want to multiply
- * with BlockVector objects, you should
- * consider using a BlockChunkSparseMatrix as
- * well.
+ * Note that while this function can operate on all vectors that offer
+ * iterator classes, it is only really effective for objects of type @ref
+ * Vector. For all classes for which iterating over elements, or random
+ * member access is expensive, this function is not efficient. In
+ * particular, if you want to multiply with BlockVector objects, you should
+ * consider using a BlockChunkSparseMatrix as well.
*
- * Source and destination must
- * not be the same vector.
+ * Source and destination must not be the same vector.
*/
template <class OutVector, class InVector>
void Tvmult_add (OutVector &dst,
const InVector &src) const;
/**
- * Return the square of the norm
- * of the vector $v$ with respect
- * to the norm induced by this
- * matrix,
- * i.e. $\left(v,Mv\right)$. This
- * is useful, e.g. in the finite
- * element context, where the
- * $L_2$ norm of a function
- * equals the matrix norm with
- * respect to the mass matrix of
- * the vector representing the
- * nodal values of the finite
- * element function.
+ * Return the square of the norm of the vector $v$ with respect to the norm
+ * induced by this matrix, i.e. $\left(v,Mv\right)$. This is useful, e.g. in
+ * the finite element context, where the $L_2$ norm of a function equals the
+ * matrix norm with respect to the mass matrix of the vector representing
+ * the nodal values of the finite element function.
*
- * Obviously, the matrix needs to be
- * quadratic for this operation, and for
- * the result to actually be a norm it
- * also needs to be either real symmetric
- * or complex hermitian.
+ * Obviously, the matrix needs to be quadratic for this operation, and for
+ * the result to actually be a norm it also needs to be either real
+ * symmetric or complex hermitian.
*
- * The underlying template types of both
- * this matrix and the given vector
- * should either both be real or
- * complex-valued, but not mixed, for
- * this function to make sense.
+ * The underlying template types of both this matrix and the given vector
+ * should either both be real or complex-valued, but not mixed, for this
+ * function to make sense.
*/
template <typename somenumber>
somenumber matrix_norm_square (const Vector<somenumber> &v) const;
/**
- * Compute the matrix scalar
- * product $\left(u,Mv\right)$.
+ * Compute the matrix scalar product $\left(u,Mv\right)$.
*/
template <typename somenumber>
somenumber matrix_scalar_product (const Vector<somenumber> &u,
const Vector<somenumber> &v) const;
/**
- * Compute the residual of an
- * equation <i>Mx=b</i>, where
- * the residual is defined to be
- * <i>r=b-Mx</i>. Write the
- * residual into
- * <tt>dst</tt>. The
- * <i>l<sub>2</sub></i> norm of
- * the residual vector is
- * returned.
+ * Compute the residual of an equation <i>Mx=b</i>, where the residual is
+ * defined to be <i>r=b-Mx</i>. Write the residual into <tt>dst</tt>. The
+ * <i>l<sub>2</sub></i> norm of the residual vector is returned.
*
- * Source <i>x</i> and destination
- * <i>dst</i> must not be the same
- * vector.
+ * Source <i>x</i> and destination <i>dst</i> must not be the same vector.
*/
template <typename somenumber>
somenumber residual (Vector<somenumber> &dst,
//@{
/**
- * Return the l1-norm of the matrix, that is
- * $|M|_1=max_{all columns j}\sum_{all
- * rows i} |M_ij|$,
- * (max. sum of columns).
- * This is the
- * natural matrix norm that is compatible
- * to the l1-norm for vectors, i.e.
- * $|Mv|_1\leq |M|_1 |v|_1$.
- * (cf. Haemmerlin-Hoffmann : Numerische Mathematik)
+ * Return the l1-norm of the matrix, that is $|M|_1=max_{all columns
+ * j}\sum_{all rows i} |M_ij|$, (max. sum of columns). This is the natural
+ * matrix norm that is compatible to the l1-norm for vectors, i.e.
+ * $|Mv|_1\leq |M|_1 |v|_1$. (cf. Haemmerlin-Hoffmann : Numerische
+ * Mathematik)
*/
real_type l1_norm () const;
/**
- * Return the linfty-norm of the
- * matrix, that is
- * $|M|_infty=max_{all rows i}\sum_{all
- * columns j} |M_ij|$,
- * (max. sum of rows).
- * This is the
- * natural matrix norm that is compatible
- * to the linfty-norm of vectors, i.e.
- * $|Mv|_infty \leq |M|_infty |v|_infty$.
- * (cf. Haemmerlin-Hoffmann : Numerische Mathematik)
+ * Return the linfty-norm of the matrix, that is $|M|_infty=max_{all rows
+ * i}\sum_{all columns j} |M_ij|$, (max. sum of rows). This is the natural
+ * matrix norm that is compatible to the linfty-norm of vectors, i.e.
+ * $|Mv|_infty \leq |M|_infty |v|_infty$. (cf. Haemmerlin-Hoffmann :
+ * Numerische Mathematik)
*/
real_type linfty_norm () const;
/**
- * Return the frobenius norm of the
- * matrix, i.e. the square root of the
- * sum of squares of all entries in the
- * matrix.
+ * Return the frobenius norm of the matrix, i.e. the square root of the sum
+ * of squares of all entries in the matrix.
*/
real_type frobenius_norm () const;
//@}
//@{
/**
- * Apply the Jacobi
- * preconditioner, which
- * multiplies every element of
- * the <tt>src</tt> vector by the
- * inverse of the respective
- * diagonal element and
- * multiplies the result with the
- * relaxation factor <tt>omega</tt>.
+ * Apply the Jacobi preconditioner, which multiplies every element of the
+ * <tt>src</tt> vector by the inverse of the respective diagonal element and
+ * multiplies the result with the relaxation factor <tt>omega</tt>.
*/
template <typename somenumber>
void precondition_Jacobi (Vector<somenumber> &dst,
const number omega = 1.) const;
/**
- * Apply SSOR preconditioning to
- * <tt>src</tt>.
+ * Apply SSOR preconditioning to <tt>src</tt>.
*/
template <typename somenumber>
void precondition_SSOR (Vector<somenumber> &dst,
const number om = 1.) const;
/**
- * Apply SOR preconditioning
- * matrix to <tt>src</tt>.
+ * Apply SOR preconditioning matrix to <tt>src</tt>.
*/
template <typename somenumber>
void precondition_SOR (Vector<somenumber> &dst,
const number om = 1.) const;
/**
- * Apply transpose SOR
- * preconditioning matrix to
- * <tt>src</tt>.
+ * Apply transpose SOR preconditioning matrix to <tt>src</tt>.
*/
template <typename somenumber>
void precondition_TSOR (Vector<somenumber> &dst,
const number om = 1.) const;
/**
- * Perform SSOR preconditioning
- * in-place. Apply the
- * preconditioner matrix without
- * copying to a second vector.
- * <tt>omega</tt> is the relaxation
+ * Perform SSOR preconditioning in-place. Apply the preconditioner matrix
+ * without copying to a second vector. <tt>omega</tt> is the relaxation
* parameter.
*/
template <typename somenumber>
const number omega = 1.) const;
/**
- * Perform an SOR preconditioning
- * in-place. <tt>omega</tt> is
- * the relaxation parameter.
+ * Perform an SOR preconditioning in-place. <tt>omega</tt> is the
+ * relaxation parameter.
*/
template <typename somenumber>
void SOR (Vector<somenumber> &v,
const number om = 1.) const;
/**
- * Perform a transpose SOR
- * preconditioning in-place.
- * <tt>omega</tt> is the
+ * Perform a transpose SOR preconditioning in-place. <tt>omega</tt> is the
* relaxation parameter.
*/
template <typename somenumber>
const number om = 1.) const;
/**
- * Perform a permuted SOR
- * preconditioning in-place.
+ * Perform a permuted SOR preconditioning in-place.
*
- * The standard SOR method is
- * applied in the order
- * prescribed by <tt>permutation</tt>,
- * that is, first the row
- * <tt>permutation[0]</tt>, then
- * <tt>permutation[1]</tt> and so
- * on. For efficiency reasons,
- * the permutation as well as its
- * inverse are required.
+ * The standard SOR method is applied in the order prescribed by
+ * <tt>permutation</tt>, that is, first the row <tt>permutation[0]</tt>,
+ * then <tt>permutation[1]</tt> and so on. For efficiency reasons, the
+ * permutation as well as its inverse are required.
*
- * <tt>omega</tt> is the
- * relaxation parameter.
+ * <tt>omega</tt> is the relaxation parameter.
*/
template <typename somenumber>
void PSOR (Vector<somenumber> &v,
const number om = 1.) const;
/**
- * Perform a transposed permuted SOR
- * preconditioning in-place.
+ * Perform a transposed permuted SOR preconditioning in-place.
*
- * The transposed SOR method is
- * applied in the order
- * prescribed by
- * <tt>permutation</tt>, that is,
- * first the row
- * <tt>permutation[m()-1]</tt>,
- * then
- * <tt>permutation[m()-2]</tt>
- * and so on. For efficiency
- * reasons, the permutation as
- * well as its inverse are
- * required.
+ * The transposed SOR method is applied in the order prescribed by
+ * <tt>permutation</tt>, that is, first the row <tt>permutation[m()-1]</tt>,
+ * then <tt>permutation[m()-2]</tt> and so on. For efficiency reasons, the
+ * permutation as well as its inverse are required.
*
- * <tt>omega</tt> is the
- * relaxation parameter.
+ * <tt>omega</tt> is the relaxation parameter.
*/
template <typename somenumber>
void TPSOR (Vector<somenumber> &v,
const number om = 1.) const;
/**
- * Do one SOR step on <tt>v</tt>.
- * Performs a direct SOR step
- * with right hand side
- * <tt>b</tt>.
+ * Do one SOR step on <tt>v</tt>. Performs a direct SOR step with right
+ * hand side <tt>b</tt>.
*/
template <typename somenumber>
void SOR_step (Vector<somenumber> &v,
const number om = 1.) const;
/**
- * Do one adjoint SOR step on
- * <tt>v</tt>. Performs a direct
- * TSOR step with right hand side
- * <tt>b</tt>.
+ * Do one adjoint SOR step on <tt>v</tt>. Performs a direct TSOR step with
+ * right hand side <tt>b</tt>.
*/
template <typename somenumber>
void TSOR_step (Vector<somenumber> &v,
const number om = 1.) const;
/**
- * Do one SSOR step on
- * <tt>v</tt>. Performs a direct
- * SSOR step with right hand side
- * <tt>b</tt> by performing TSOR
- * after SOR.
+ * Do one SSOR step on <tt>v</tt>. Performs a direct SSOR step with right
+ * hand side <tt>b</tt> by performing TSOR after SOR.
*/
template <typename somenumber>
void SSOR_step (Vector<somenumber> &v,
const Vector<somenumber> &b,
const number om = 1.) const;
+//@}
+ /**
+ * @name Iterators
+ */
+//@{
+
+ /**
+ * STL-like iterator with the first entry of the matrix. This is the version
+ * for constant matrices.
+ *
+ * Note that due to the layout in ChunkSparseMatrix, iterating over matrix
+ * entries is considerably slower than for a sparse matrix, as the iterator
+ * is travels row-by-row, whereas data is stored in chunks of several rows
+ * and columns.
+ */
+ const_iterator begin () const;
+
+ /**
+ * Final iterator. This is the version for constant matrices.
+ *
+ * Note that due to the layout in ChunkSparseMatrix, iterating over matrix
+ * entries is considerably slower than for a sparse matrix, as the iterator
+ * is travels row-by-row, whereas data is stored in chunks of several rows
+ * and columns.
+ */
+ const_iterator end () const;
+
+ /**
+ * STL-like iterator with the first entry of the matrix. This is the version
+ * for non-constant matrices.
+ *
+ * Note that due to the layout in ChunkSparseMatrix, iterating over matrix
+ * entries is considerably slower than for a sparse matrix, as the iterator
+ * is travels row-by-row, whereas data is stored in chunks of several rows
+ * and columns.
+ */
+ iterator begin ();
+
+ /**
+ * Final iterator. This is the version for non-constant matrices.
+ *
+ * Note that due to the layout in ChunkSparseMatrix, iterating over matrix
+ * entries is considerably slower than for a sparse matrix, as the iterator
+ * is travels row-by-row, whereas data is stored in chunks of several rows
+ * and columns.
+ */
+ iterator end ();
+
+ /**
+ * STL-like iterator with the first entry of row <tt>r</tt>. This is the
+ * version for constant matrices.
+ *
+ * Note that if the given row is empty, i.e. does not contain any nonzero
+ * entries, then the iterator returned by this function equals
+ * <tt>end(r)</tt>. Note also that the iterator may not be dereferencable in
+ * that case.
+ *
+ * Note that due to the layout in ChunkSparseMatrix, iterating over matrix
+ * entries is considerably slower than for a sparse matrix, as the iterator
+ * is travels row-by-row, whereas data is stored in chunks of several rows
+ * and columns.
+ */
+ const_iterator begin (const unsigned int r) const;
+
+ /**
+ * Final iterator of row <tt>r</tt>. It points to the first element past the
+ * end of line @p r, or past the end of the entire sparsity pattern. This is
+ * the version for constant matrices.
+ *
+ * Note that the end iterator is not necessarily dereferencable. This is in
+ * particular the case if it is the end iterator for the last row of a
+ * matrix.
+ *
+ * Note that due to the layout in ChunkSparseMatrix, iterating over matrix
+ * entries is considerably slower than for a sparse matrix, as the iterator
+ * is travels row-by-row, whereas data is stored in chunks of several rows
+ * and columns.
+ */
+ const_iterator end (const unsigned int r) const;
+
+ /**
+ * STL-like iterator with the first entry of row <tt>r</tt>. This is the
+ * version for non-constant matrices.
+ *
+ * Note that if the given row is empty, i.e. does not contain any nonzero
+ * entries, then the iterator returned by this function equals
+ * <tt>end(r)</tt>. Note also that the iterator may not be dereferencable in
+ * that case.
+ *
+ * Note that due to the layout in ChunkSparseMatrix, iterating over matrix
+ * entries is considerably slower than for a sparse matrix, as the iterator
+ * is travels row-by-row, whereas data is stored in chunks of several rows
+ * and columns.
+ */
+ iterator begin (const unsigned int r);
+
+ /**
+ * Final iterator of row <tt>r</tt>. It points to the first element past the
+ * end of line @p r, or past the end of the entire sparsity pattern. This is
+ * the version for non-constant matrices.
+ *
+ * Note that the end iterator is not necessarily dereferencable. This is in
+ * particular the case if it is the end iterator for the last row of a
+ * matrix.
+ *
+ * Note that due to the layout in ChunkSparseMatrix, iterating over matrix
+ * entries is considerably slower than for a sparse matrix, as the iterator
+ * is travels row-by-row, whereas data is stored in chunks of several rows
+ * and columns.
+ */
+ iterator end (const unsigned int r);
//@}
/**
* @name Input/Output
//@{
/**
- * Print the matrix to the given
- * stream, using the format
- * <tt>(line,col) value</tt>,
- * i.e. one nonzero entry of the
- * matrix per line.
+ * Print the matrix to the given stream, using the format <tt>(line,col)
+ * value</tt>, i.e. one nonzero entry of the matrix per line.
*/
void print (std::ostream &out) const;
/**
- * Print the matrix in the usual
- * format, i.e. as a matrix and
- * not as a list of nonzero
- * elements. For better
- * readability, elements not in
- * the matrix are displayed as
- * empty space, while matrix
- * elements which are explicitly
- * set to zero are displayed as
- * such.
+ * Print the matrix in the usual format, i.e. as a matrix and not as a list
+ * of nonzero elements. For better readability, elements not in the matrix
+ * are displayed as empty space, while matrix elements which are explicitly
+ * set to zero are displayed as such.
*
- * The parameters allow for a
- * flexible setting of the output
- * format: <tt>precision</tt> and
- * <tt>scientific</tt> are used
- * to determine the number
- * format, where <tt>scientific =
- * false</tt> means fixed point
- * notation. A zero entry for
- * <tt>width</tt> makes the
- * function compute a width, but
- * it may be changed to a
- * positive value, if output is
- * crude.
+ * The parameters allow for a flexible setting of the output format:
+ * <tt>precision</tt> and <tt>scientific</tt> are used to determine the
+ * number format, where <tt>scientific = false</tt> means fixed point
+ * notation. A zero entry for <tt>width</tt> makes the function compute a
+ * width, but it may be changed to a positive value, if output is crude.
*
- * Additionally, a character for
- * an empty value may be
- * specified.
+ * Additionally, a character for an empty value may be specified.
*
- * Finally, the whole matrix can
- * be multiplied with a common
- * denominator to produce more
- * readable output, even
- * integers.
+ * Finally, the whole matrix can be multiplied with a common denominator to
+ * produce more readable output, even integers.
*
- * @attention This function may
- * produce <b>large</b> amounts
- * of output if applied to a
- * large matrix!
+ * @attention This function may produce <b>large</b> amounts of output if
+ * applied to a large matrix!
*/
void print_formatted (std::ostream &out,
const unsigned int precision = 3,
const double denominator = 1.) const;
/**
- * Print the actual pattern of
- * the matrix. For each entry
- * with an absolute value larger
- * than threshold, a '*' is
- * printed, a ':' for every value
- * smaller and a '.' for every
- * entry not allocated.
+ * Print the actual pattern of the matrix. For each entry with an absolute
+ * value larger than threshold, a '*' is printed, a ':' for every value
+ * smaller and a '.' for every entry not allocated.
*/
void print_pattern(std::ostream &out,
const double threshold = 0.) const;
/**
- * Write the data of this object
- * en bloc to a file. This is
- * done in a binary mode, so the
- * output is neither readable by
- * humans nor (probably) by other
- * computers using a different
- * operating system of number
- * format.
+ * Write the data of this object en bloc to a file. This is done in a binary
+ * mode, so the output is neither readable by humans nor (probably) by other
+ * computers using a different operating system or number format.
*
- * The purpose of this function
- * is that you can swap out
- * matrices and sparsity pattern
- * if you are short of memory,
- * want to communicate between
- * different programs, or allow
- * objects to be persistent
- * across different runs of the
- * program.
+ * The purpose of this function is that you can swap out matrices and
+ * sparsity pattern if you are short of memory, want to communicate between
+ * different programs, or allow objects to be persistent across different
+ * runs of the program.
*/
void block_write (std::ostream &out) const;
/**
- * Read data that has previously
- * been written by block_write()
- * from a file. This is done
- * using the inverse operations
- * to the above function, so it
- * is reasonably fast because the
- * bitstream is not interpreted
- * except for a few numbers up
- * front.
+ * Read data that has previously been written by block_write() from a
+ * file. This is done using the inverse operations to the above function, so
+ * it is reasonably fast because the bitstream is not interpreted except for
+ * a few numbers up front.
*
- * The object is resized on this
- * operation, and all previous
- * contents are lost. Note,
- * however, that no checks are
- * performed whether new data and
- * the underlying ChunkSparsityPattern
- * object fit together. It is
- * your responsibility to make
- * sure that the sparsity pattern
- * and the data to be read match.
+ * The object is resized on this operation, and all previous contents are
+ * lost. Note, however, that no checks are performed whether new data and
+ * the underlying ChunkSparsityPattern object fit together. It is your
+ * responsibility to make sure that the sparsity pattern and the data to be
+ * read match.
*
- * A primitive form of error
- * checking is performed which
- * will recognize the bluntest
- * attempts to interpret some
- * data as a matrix stored
- * bitwise to a file that wasn't
- * actually created that way, but
- * not more.
+ * A primitive form of error checking is performed which will recognize the
+ * bluntest attempts to interpret some data as a matrix stored bitwise to a
+ * file that wasn't actually created that way, but not more.
*/
void block_read (std::istream &in);
//@}
//@}
private:
/**
- * Pointer to the sparsity
- * pattern used for this
- * matrix. In order to guarantee
- * that it is not deleted while
- * still in use, we subscribe to
- * it using the SmartPointer
- * class.
+ * Pointer to the sparsity pattern used for this matrix. In order to
+ * guarantee that it is not deleted while still in use, we subscribe to it
+ * using the SmartPointer class.
*/
SmartPointer<const ChunkSparsityPattern,ChunkSparseMatrix<number> > cols;
/**
- * Array of values for all the
- * nonzero entries. The position
- * within the matrix, i.e. the
- * row and column number for a
- * given entry can only be
- * deduced using the sparsity
- * pattern. The same holds for
- * the more common operation of
- * finding an entry by its
- * coordinates.
+ * Array of values for all the nonzero entries. The position within the
+ * matrix, i.e. the row and column number for a given entry can only be
+ * deduced using the sparsity pattern. The same holds for the more common
+ * operation of finding an entry by its coordinates.
*/
number *val;
/**
- * Allocated size of #val. This
- * can be larger than the
- * actually used part if the size
- * of the matrix was reduced
- * somewhen in the past by
- * associating a sparsity pattern
- * with a smaller size to this
- * object, using the reinit()
+ * Allocated size of #val. This can be larger than the actually used part if
+ * the size of the matrix was reduced somewhen in the past by associating a
+ * sparsity pattern with a smaller size to this object, using the reinit()
* function.
*/
unsigned int max_len;
/**
- * Return the location of entry
- * $(i,j)$ within the val array.
+ * Return the location of entry $(i,j)$ within the val array.
*/
unsigned int compute_location (const unsigned int i,
const unsigned int j) const;
- // make all other sparse matrices
- // friends
+ // make all other sparse matrices friends
template <typename somenumber> friend class ChunkSparseMatrix;
+
+ /**
+ * Also give access to internal details to the iterator/accessor
+ * classes.
+ */
+ template <typename,bool> friend class ChunkSparseMatrixIterators::Iterator;
+ template <typename,bool> friend class ChunkSparseMatrixIterators::Accessor;
};
/*@}*/
+template <typename number>
+inline
+const ChunkSparsityPattern &
+ChunkSparseMatrix<number>::get_sparsity_pattern () const
+{
+ Assert (cols != 0, ExcNotInitialized());
+ return *cols;
+}
+
+
+
template <typename number>
inline
unsigned int
{
const unsigned int chunk_size = cols->get_chunk_size();
const unsigned int chunk_index
- = cols->sparsity_pattern(i/chunk_size, j/chunk_size);
+ = cols->sparsity_pattern(i/chunk_size, j/chunk_size);
if (chunk_index == ChunkSparsityPattern::invalid_entry)
return ChunkSparsityPattern::invalid_entry;
Assert (numbers::is_finite(value), ExcNumberNotFinite());
Assert (cols != 0, ExcNotInitialized());
- // it is allowed to set elements of
- // the matrix that are not part of
- // the sparsity pattern, if the
- // value to which we set it is zero
+ // it is allowed to set elements of the matrix that are not part of the
+ // sparsity pattern, if the value to which we set it is zero
const unsigned int index = compute_location(i,j);
Assert ((index != SparsityPattern::invalid_entry) ||
(value == 0.),
const unsigned int chunk_size = cols->get_chunk_size();
- // multiply all elements of the matrix with
- // the given factor. this includes the
- // padding elements in chunks that overlap
- // the boundaries of the actual matrix --
- // but since multiplication with a number
- // does not violate the invariant of
- // keeping these elements at zero nothing
- // can happen
+ // multiply all elements of the matrix with the given factor. this includes
+ // the padding elements in chunks that overlap the boundaries of the actual
+ // matrix -- but since multiplication with a number does not violate the
+ // invariant of keeping these elements at zero nothing can happen
number *val_ptr = val;
const number *const end_ptr = val +
cols->sparsity_pattern.n_nonzero_elements()
const unsigned int chunk_size = cols->get_chunk_size();
- // multiply all elements of the matrix with
- // the given factor. this includes the
- // padding elements in chunks that overlap
- // the boundaries of the actual matrix --
- // but since multiplication with a number
- // does not violate the invariant of
- // keeping these elements at zero nothing
- // can happen
+ // multiply all elements of the matrix with the given factor. this includes
+ // the padding elements in chunks that overlap the boundaries of the actual
+ // matrix -- but since multiplication with a number does not violate the
+ // invariant of keeping these elements at zero nothing can happen
number *val_ptr = val;
const number *const end_ptr = val +
cols->sparsity_pattern.n_nonzero_elements()
Assert (m() == n(), ExcNotQuadratic());
Assert (i<m(), ExcInvalidIndex1(i));
- // Use that the first element in each row
- // of a quadratic matrix is the main
- // diagonal of the chunk sparsity pattern
- const unsigned int chunk_size = cols->get_chunk_size();
- return val[cols->sparsity_pattern.rowstart[i/chunk_size]
- *
- chunk_size * chunk_size
- +
- (i % chunk_size) * chunk_size
- +
- (i % chunk_size)];
-}
-
-
-
-template <typename number>
-inline
-number &ChunkSparseMatrix<number>::diag_element (const unsigned int i)
-{
- Assert (cols != 0, ExcNotInitialized());
- Assert (m() == n(), ExcNotQuadratic());
- Assert (i<m(), ExcInvalidIndex1(i));
-
- // Use that the first element in each row
- // of a quadratic matrix is the main
+ // Use that the first element in each row of a quadratic matrix is the main
// diagonal of the chunk sparsity pattern
const unsigned int chunk_size = cols->get_chunk_size();
return val[cols->sparsity_pattern.rowstart[i/chunk_size]
template <typename number>
template <typename ForwardIterator>
+inline
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
+ // 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)
}
+
+//---------------------------------------------------------------------------
+
+
+namespace ChunkSparseMatrixIterators
+{
+ template <typename number>
+ inline
+ Accessor<number,true>::
+ Accessor (const MatrixType *matrix,
+ const unsigned int row)
+ :
+ ChunkSparsityPatternIterators::Accessor (&matrix->get_sparsity_pattern(),
+ row),
+ matrix (matrix)
+ {}
+
+
+
+ template <typename number>
+ inline
+ Accessor<number,true>::
+ Accessor (const MatrixType *matrix)
+ :
+ ChunkSparsityPatternIterators::Accessor (&matrix->get_sparsity_pattern()),
+ matrix (matrix)
+ {}
+
+
+
+ template <typename number>
+ inline
+ Accessor<number,true>::
+ Accessor (const ChunkSparseMatrixIterators::Accessor<number,false> &a)
+ :
+ ChunkSparsityPatternIterators::Accessor (a),
+ matrix (&a.get_matrix())
+ {}
+
+
+
+ template <typename number>
+ inline
+ number
+ Accessor<number, true>::value () const
+ {
+ const unsigned int chunk_size = matrix->get_sparsity_pattern().get_chunk_size();
+ return matrix->val[reduced_index() * chunk_size * chunk_size
+ +
+ chunk_row * chunk_size
+ +
+ chunk_col];
+ }
+
+
+
+ template <typename number>
+ inline
+ typename Accessor<number, true>::MatrixType &
+ Accessor<number, true>::get_matrix () const
+ {
+ return *matrix;
+ }
+
+
+
+ template <typename number>
+ inline
+ Accessor<number, false>::Reference::Reference (
+ const Accessor *accessor,
+ const bool)
+ :
+ accessor (accessor)
+ {}
+
+
+ template <typename number>
+ inline
+ Accessor<number, false>::Reference::operator number() const
+ {
+ const unsigned int chunk_size = accessor->matrix->get_sparsity_pattern().get_chunk_size();
+ return accessor->matrix->val[accessor->reduced_index() * chunk_size * chunk_size
+ +
+ accessor->chunk_row * chunk_size
+ +
+ accessor->chunk_col];
+ }
+
+
+
+ template <typename number>
+ inline
+ const typename Accessor<number, false>::Reference &
+ Accessor<number, false>::Reference::operator = (const number n) const
+ {
+ const unsigned int chunk_size = accessor->matrix->get_sparsity_pattern().get_chunk_size();
+ accessor->matrix->val[accessor->reduced_index() * chunk_size * chunk_size
+ +
+ accessor->chunk_row * chunk_size
+ +
+ accessor->chunk_col] = n;
+ return *this;
+ }
+
+
+
+ template <typename number>
+ inline
+ const typename Accessor<number, false>::Reference &
+ Accessor<number, false>::Reference::operator += (const number n) const
+ {
+ const unsigned int chunk_size = accessor->matrix->get_sparsity_pattern().get_chunk_size();
+ accessor->matrix->val[accessor->reduced_index() * chunk_size * chunk_size
+ +
+ accessor->chunk_row * chunk_size
+ +
+ accessor->chunk_col] += n;
+ return *this;
+ }
+
+
+
+ template <typename number>
+ inline
+ const typename Accessor<number, false>::Reference &
+ Accessor<number, false>::Reference::operator -= (const number n) const
+ {
+ const unsigned int chunk_size = accessor->matrix->get_sparsity_pattern().get_chunk_size();
+ accessor->matrix->val[accessor->reduced_index() * chunk_size * chunk_size
+ +
+ accessor->chunk_row * chunk_size
+ +
+ accessor->chunk_col] -= n;
+ return *this;
+ }
+
+
+
+ template <typename number>
+ inline
+ const typename Accessor<number, false>::Reference &
+ Accessor<number, false>::Reference::operator *= (const number n) const
+ {
+ const unsigned int chunk_size = accessor->matrix->get_sparsity_pattern().get_chunk_size();
+ accessor->matrix->val[accessor->reduced_index() * chunk_size * chunk_size
+ +
+ accessor->chunk_row * chunk_size
+ +
+ accessor->chunk_col] *= n;
+ return *this;
+ }
+
+
+
+ template <typename number>
+ inline
+ const typename Accessor<number, false>::Reference &
+ Accessor<number, false>::Reference::operator /= (const number n) const
+ {
+ const unsigned int chunk_size = accessor->matrix->get_sparsity_pattern().get_chunk_size();
+ accessor->matrix->val[accessor->reduced_index() * chunk_size * chunk_size
+ +
+ accessor->chunk_row * chunk_size
+ +
+ accessor->chunk_col] /= n;
+ return *this;
+ }
+
+
+
+ template <typename number>
+ inline
+ Accessor<number,false>::
+ Accessor (MatrixType *matrix,
+ const unsigned int row)
+ :
+ ChunkSparsityPatternIterators::Accessor (&matrix->get_sparsity_pattern(),
+ row),
+ matrix (matrix)
+ {}
+
+
+
+ template <typename number>
+ inline
+ Accessor<number,false>::
+ Accessor (MatrixType *matrix)
+ :
+ ChunkSparsityPatternIterators::Accessor (&matrix->get_sparsity_pattern()),
+ matrix (matrix)
+ {}
+
+
+
+ template <typename number>
+ inline
+ typename Accessor<number, false>::Reference
+ Accessor<number, false>::value() const
+ {
+ return Reference(this,true);
+ }
+
+
+
+
+ template <typename number>
+ inline
+ typename Accessor<number, false>::MatrixType &
+ Accessor<number, false>::get_matrix () const
+ {
+ return *matrix;
+ }
+
+
+
+ template <typename number, bool Constness>
+ inline
+ Iterator<number, Constness>::
+ Iterator (MatrixType *matrix,
+ const unsigned int row)
+ :
+ accessor(matrix, row)
+ {}
+
+
+
+ template <typename number, bool Constness>
+ inline
+ Iterator<number, Constness>::
+ Iterator (MatrixType *matrix)
+ :
+ accessor(matrix)
+ {}
+
+
+
+ template <typename number, bool Constness>
+ inline
+ Iterator<number, Constness>::
+ Iterator (const ChunkSparseMatrixIterators::Iterator<number,false> &i)
+ :
+ accessor(*i)
+ {}
+
+
+
+ template <typename number, bool Constness>
+ inline
+ Iterator<number, Constness> &
+ Iterator<number,Constness>::operator++ ()
+ {
+ accessor.advance ();
+ return *this;
+ }
+
+
+ template <typename number, bool Constness>
+ inline
+ Iterator<number,Constness>
+ Iterator<number,Constness>::operator++ (int)
+ {
+ const Iterator iter = *this;
+ accessor.advance ();
+ return iter;
+ }
+
+
+ template <typename number, bool Constness>
+ inline
+ const Accessor<number,Constness> &
+ Iterator<number,Constness>::operator* () const
+ {
+ return accessor;
+ }
+
+
+ template <typename number, bool Constness>
+ inline
+ const Accessor<number,Constness> *
+ Iterator<number,Constness>::operator-> () const
+ {
+ return &accessor;
+ }
+
+
+ template <typename number, bool Constness>
+ inline
+ bool
+ Iterator<number,Constness>::
+ operator == (const Iterator &other) const
+ {
+ return (accessor == other.accessor);
+ }
+
+
+ template <typename number, bool Constness>
+ inline
+ bool
+ Iterator<number,Constness>::
+ operator != (const Iterator &other) const
+ {
+ return ! (*this == other);
+ }
+
+
+ template <typename number, bool Constness>
+ inline
+ bool
+ Iterator<number,Constness>::
+ operator < (const Iterator &other) const
+ {
+ Assert (&accessor.get_matrix() == &other.accessor.get_matrix(),
+ ExcInternalError());
+
+ return (accessor < other.accessor);
+ }
+
+
+ template <typename number, bool Constness>
+ inline
+ bool
+ Iterator<number,Constness>::
+ operator > (const Iterator &other) const
+ {
+ return (other < *this);
+ }
+
+
+ template <typename number, bool Constness>
+ inline
+ int
+ Iterator<number,Constness>::
+ operator - (const Iterator &other) const
+ {
+ Assert (&accessor.get_matrix() == &other.accessor.get_matrix(),
+ ExcInternalError());
+
+ // TODO: can be optimized
+ int difference = 0;
+ if (*this < other)
+ {
+ Iterator copy = *this;
+ while (copy < other)
+ {
+ ++copy;
+ --difference;
+ }
+ }
+ else
+ {
+ Iterator copy = other;
+ while (copy < *this)
+ {
+ ++copy;
+ ++difference;
+ }
+ }
+ return difference;
+ }
+
+
+
+ template <typename number, bool Constness>
+ inline
+ Iterator<number,Constness>
+ Iterator<number,Constness>::
+ operator + (const unsigned int n) const
+ {
+ Iterator x = *this;
+ for (unsigned int i=0; i<n; ++i)
+ ++x;
+
+ return x;
+ }
+
+}
+
+
+
+template <typename number>
+inline
+typename ChunkSparseMatrix<number>::const_iterator
+ChunkSparseMatrix<number>::begin () const
+{
+ return const_iterator(this, 0);
+}
+
+
+template <typename number>
+inline
+typename ChunkSparseMatrix<number>::const_iterator
+ChunkSparseMatrix<number>::end () const
+{
+ return const_iterator(this);
+}
+
+
+template <typename number>
+inline
+typename ChunkSparseMatrix<number>::iterator
+ChunkSparseMatrix<number>::begin ()
+{
+ return iterator(this, 0);
+}
+
+
+template <typename number>
+inline
+typename ChunkSparseMatrix<number>::iterator
+ChunkSparseMatrix<number>::end ()
+{
+ return iterator(this);
+}
+
+
+template <typename number>
+inline
+typename ChunkSparseMatrix<number>::const_iterator
+ChunkSparseMatrix<number>::begin (const unsigned int r) const
+{
+ Assert (r<m(), ExcIndexRange(r,0,m()));
+ return const_iterator(this, r);
+}
+
+
+
+template <typename number>
+inline
+typename ChunkSparseMatrix<number>::const_iterator
+ChunkSparseMatrix<number>::end (const unsigned int r) const
+{
+ Assert (r<m(), ExcIndexRange(r,0,m()));
+ return const_iterator(this, r+1);
+}
+
+
+
+template <typename number>
+inline
+typename ChunkSparseMatrix<number>::iterator
+ChunkSparseMatrix<number>::begin (const unsigned int r)
+{
+ Assert (r<m(), ExcIndexRange(r,0,m()));
+ return iterator(this, r);
+}
+
+
+
+template <typename number>
+inline
+typename ChunkSparseMatrix<number>::iterator
+ChunkSparseMatrix<number>::end (const unsigned int r)
+{
+ Assert (r<m(), ExcIndexRange(r,0,m()));
+ return iterator(this, r+1);
+}
+
+
+
+
#endif // DOXYGEN
#include <deal.II/base/template_constraints.h>
+#include <deal.II/base/parallel.h>
#include <deal.II/lac/chunk_sparse_matrix.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/full_matrix.h>
namespace ChunkSparseMatrix
{
/**
- * Add the result of multiplying a chunk
- * of size chunk_size times chunk_size by
- * a source vector fragment of size
- * chunk_size to the destination vector
- * fragment.
+ * Add the result of multiplying a chunk of size chunk_size times
+ * chunk_size by a source vector fragment of size chunk_size to the
+ * destination vector fragment.
*/
template <typename MatrixIterator,
typename SrcIterator,
/**
- * Like the previous function, but
- * subtract. We need this for computing
+ * Like the previous function, but subtract. We need this for computing
* the residual.
*/
template <typename MatrixIterator,
/**
- * Add the result of multiplying the
- * transpose of a chunk of size
- * chunk_size times chunk_size by a
- * source vector fragment of size
- * chunk_size to the destination vector
- * fragment.
+ * Add the result of multiplying the transpose of a chunk of size
+ * chunk_size times chunk_size by a source vector fragment of size
+ * chunk_size to the destination vector fragment.
*/
template <typename MatrixIterator,
typename SrcIterator,
/**
- * Produce the result of the matrix
- * scalar product $u^TMv$ for an
+ * Produce the result of the matrix scalar product $u^TMv$ for an
* individual chunk.
*/
template <typename result_type,
return result;
}
+
+
+
+ /**
+ * Perform a vmult_add using the ChunkSparseMatrix data structures, but
+ * only using a subinterval of the matrix rows.
+ *
+ * In the sequential case, this function is called on all rows, in the
+ * parallel case it may be called on a subrange, at the discretion of the
+ * task scheduler.
+ */
+ template <typename number,
+ typename InVector,
+ typename OutVector>
+ void vmult_add_on_subrange (const ChunkSparsityPattern &cols,
+ const unsigned int begin_row,
+ const unsigned int end_row,
+ const number *values,
+ const std::size_t *rowstart,
+ const unsigned int *colnums,
+ const InVector &src,
+ OutVector &dst)
+ {
+ const unsigned int m = cols.n_rows();
+ const unsigned int n = cols.n_cols();
+ const unsigned int chunk_size = cols.get_chunk_size();
+
+ // loop over all chunks. note that we need to treat the last chunk row
+ // and column differently if they have padding elements
+ const unsigned int n_filled_last_rows = m % chunk_size;
+ const unsigned int n_filled_last_cols = n % chunk_size;
+
+ const unsigned int last_regular_row = n_filled_last_rows > 0 ?
+ std::min(m/chunk_size, end_row) : end_row;
+ const unsigned int irregular_col = n/chunk_size;
+
+ typename OutVector::iterator dst_ptr = dst.begin()+chunk_size*begin_row;
+ const number *val_ptr= &values[rowstart[begin_row]*chunk_size*chunk_size];
+ const unsigned int *colnum_ptr = &colnums[rowstart[begin_row]];
+ for (unsigned int chunk_row=begin_row; chunk_row<last_regular_row;
+ ++chunk_row)
+ {
+ const number *const val_end_of_row = &values[rowstart[chunk_row+1] *
+ chunk_size * chunk_size];
+ while (val_ptr != val_end_of_row)
+ {
+ if (*colnum_ptr != irregular_col)
+ chunk_vmult_add (chunk_size,
+ val_ptr,
+ src.begin() + *colnum_ptr * chunk_size,
+ dst_ptr);
+ else
+ // we're at a chunk column that has padding
+ for (unsigned int r=0; r<chunk_size; ++r)
+ for (unsigned int c=0; c<n_filled_last_cols; ++c)
+ dst_ptr[r] += (val_ptr[r*chunk_size + c] *
+ src(*colnum_ptr * chunk_size + c));
+
+ ++colnum_ptr;
+ val_ptr += chunk_size * chunk_size;
+ }
+
+ dst_ptr += chunk_size;
+ }
+
+ // now deal with last chunk row if necessary
+ if (n_filled_last_rows > 0 && end_row == (m/chunk_size+1))
+ {
+ const unsigned int chunk_row = last_regular_row;
+
+ const number *const val_end_of_row = &values[rowstart[chunk_row+1] *
+ chunk_size * chunk_size];
+ while (val_ptr != val_end_of_row)
+ {
+ if (*colnum_ptr != irregular_col)
+ {
+ // we're at a chunk row but not column that has padding
+ for (unsigned int r=0; r<n_filled_last_rows; ++r)
+ for (unsigned int c=0; c<chunk_size; ++c)
+ dst_ptr[r]
+ += (val_ptr[r*chunk_size + c] *
+ src(*colnum_ptr * chunk_size + c));
+ }
+ else
+ // we're at a chunk row and column that has padding
+ for (unsigned int r=0; r<n_filled_last_rows; ++r)
+ for (unsigned int c=0; c<n_filled_last_cols; ++c)
+ dst_ptr[r]
+ += (val_ptr[r*chunk_size + c] *
+ src(*colnum_ptr * chunk_size + c));
+
+ ++colnum_ptr;
+ val_ptr += chunk_size * chunk_size;
+ }
+ }
+ Assert(std::size_t(colnum_ptr-&colnums[0]) == rowstart[end_row],
+ ExcInternalError());
+ Assert(std::size_t(val_ptr-&values[0]) ==
+ rowstart[end_row] * chunk_size * chunk_size,
+ ExcInternalError());
+ }
}
}
+namespace internal
+{
+ namespace ChunkSparseMatrix
+ {
+ template<typename T>
+ void zero_subrange (const unsigned int begin,
+ const unsigned int end,
+ T *dst)
+ {
+ std::memset (dst+begin,0,(end-begin)*sizeof(T));
+ }
+ }
+}
+
+
+
template <typename number>
ChunkSparseMatrix<number> &
ChunkSparseMatrix<number>::operator = (const double d)
Assert (cols->sparsity_pattern.compressed || cols->empty(),
ChunkSparsityPattern::ExcNotCompressed());
- if (val)
- {
- const unsigned int chunk_size = cols->get_chunk_size();
- std::fill_n (val,
- cols->sparsity_pattern.n_nonzero_elements() *
- chunk_size * chunk_size,
- 0.);
- }
+ // do initial zeroing of elements in parallel. Try to achieve a similar
+ // layout as when doing matrix-vector products, as on some NUMA systems, a
+ // memory block is assigned to memory banks where the first access is
+ // generated. For sparse matrices, the first operations is usually the
+ // operator=. The grain size is chosen to reflect the number of rows in
+ // minimum_parallel_grain_size, weighted by the number of nonzero entries
+ // per row on average.
+ const unsigned int matrix_size = cols->sparsity_pattern.n_nonzero_elements()
+ * cols->chunk_size * cols->chunk_size;
+ const unsigned int grain_size =
+ internal::SparseMatrix::minimum_parallel_grain_size *
+ (matrix_size+m()) / m();
+ if (matrix_size>grain_size)
+ parallel::apply_to_subranges (0U, matrix_size,
+ std_cxx1x::bind(&internal::ChunkSparseMatrix::template
+ zero_subrange<number>,
+ std_cxx1x::_1, std_cxx1x::_2,
+ val),
+ grain_size);
+ else if (matrix_size > 0)
+ std::memset (&val[0], 0, matrix_size*sizeof(number));
return *this;
}
return;
}
- // allocate not just m() * n() elements but
- // enough so that we can store full
- // chunks. this entails some padding
- // elements
+ // allocate not just m() * n() elements but enough so that we can store full
+ // chunks. this entails some padding elements
const unsigned int chunk_size = cols->get_chunk_size();
const unsigned int N = cols->sparsity_pattern.n_nonzero_elements() *
chunk_size * chunk_size;
max_len = N;
}
- // fill with zeros. do not just fill N
- // elements but all that we allocated to
- // ensure that also the padding elements
- // are zero and not left at previous values
- if (val != 0)
- std::fill_n (&val[0], max_len, 0);
+ // fill with zeros. do not just fill N elements but all that we allocated to
+ // ensure that also the padding elements are zero and not left at previous
+ // values
+ this->operator=(0.);
}
{
Assert (cols != 0, ExcNotInitialized());
- // count those elements that are nonzero,
- // even if they lie in the padding around
- // the matrix. since we have the invariant
- // that padding elements are zero, nothing
- // bad can happen here
+ // count those elements that are nonzero, even if they lie in the padding
+ // around the matrix. since we have the invariant that padding elements are
+ // zero, nothing bad can happen here
const unsigned int chunk_size = cols->get_chunk_size();
return std::count_if(&val[0],
&val[cols->sparsity_pattern.n_nonzero_elements () *
Assert (val != 0, ExcNotInitialized());
Assert (cols == matrix.cols, ExcDifferentChunkSparsityPatterns());
- // copy everything, including padding
- // elements
+ // copy everything, including padding elements
const unsigned int chunk_size = cols->get_chunk_size();
std::copy (&matrix.val[0],
&matrix.val[cols->sparsity_pattern.n_nonzero_elements()
Assert (val != 0, ExcNotInitialized());
Assert (cols == matrix.cols, ExcDifferentChunkSparsityPatterns());
- // add everything, including padding
- // elements
+ // add everything, including padding elements
const unsigned int chunk_size = cols->get_chunk_size();
number *val_ptr = &val[0];
const somenumber *matrix_ptr = &matrix.val[0];
Assert (!PointerComparison::equal(&src, &dst), ExcSourceEqualsDestination());
- // set the output vector to zero and then
- // add to it the contributions of vmults
- // from individual chunks. this is what
- // vmult_add does
+ // set the output vector to zero and then add to it the contributions of
+ // vmults from individual chunks. this is what vmult_add does
dst = 0;
vmult_add (dst, src);
}
Assert (!PointerComparison::equal(&src, &dst), ExcSourceEqualsDestination());
- // set the output vector to zero and then
- // add to it the contributions of vmults
- // from individual chunks. this is what
- // vmult_add does
+ // set the output vector to zero and then add to it the contributions of
+ // vmults from individual chunks. this is what vmult_add does
dst = 0;
Tvmult_add (dst, src);
}
Assert(n() == src.size(), ExcDimensionMismatch(n(),src.size()));
Assert (!PointerComparison::equal(&src, &dst), ExcSourceEqualsDestination());
+ parallel::apply_to_subranges (0U, cols->sparsity_pattern.n_rows(),
+ std_cxx1x::bind (&internal::ChunkSparseMatrix::vmult_add_on_subrange
+ <number,InVector,OutVector>,
+ std_cxx1x::cref(*cols),
+ std_cxx1x::_1, std_cxx1x::_2,
+ val,
+ cols->sparsity_pattern.rowstart,
+ cols->sparsity_pattern.colnums,
+ std_cxx1x::cref(src),
+ std_cxx1x::ref(dst)),
+ internal::SparseMatrix::minimum_parallel_grain_size/cols->chunk_size+1);
- const unsigned int n_chunk_rows = cols->sparsity_pattern.n_rows();
-
- // loop over all chunks. note that we need
- // to treat the last chunk row and column
- // differently if they have padding
- // elements
- const bool rows_have_padding = (m() % cols->chunk_size != 0),
- cols_have_padding = (n() % cols->chunk_size != 0);
-
- const unsigned int n_regular_chunk_rows
- = (rows_have_padding ?
- n_chunk_rows-1 :
- n_chunk_rows);
-
- const number *val_ptr = val;
- const unsigned int *colnum_ptr = cols->sparsity_pattern.colnums;
- typename OutVector::iterator dst_ptr = dst.begin();
-
- for (unsigned int chunk_row=0; chunk_row<n_regular_chunk_rows; ++chunk_row)
- {
- const number *const val_end_of_row = &val[cols->sparsity_pattern.rowstart[chunk_row+1]
- * cols->chunk_size
- * cols->chunk_size];
- while (val_ptr != val_end_of_row)
- {
- if ((cols_have_padding == false)
- ||
- (*colnum_ptr != cols->sparsity_pattern.n_cols()-1))
- internal::ChunkSparseMatrix::chunk_vmult_add
- (cols->chunk_size,
- val_ptr,
- src.begin() + *colnum_ptr * cols->chunk_size,
- dst_ptr);
- else
- // we're at a chunk column that
- // has padding
- for (unsigned int r=0; r<cols->chunk_size; ++r)
- for (unsigned int c=0; c<n() % cols->chunk_size; ++c)
- dst(chunk_row * cols->chunk_size + r)
- += (val_ptr[r*cols->chunk_size + c] *
- src(*colnum_ptr * cols->chunk_size + c));
-
- ++colnum_ptr;
- val_ptr += cols->chunk_size * cols->chunk_size;
- }
-
-
- dst_ptr += cols->chunk_size;
- }
-
- // now deal with last chunk row if
- // necessary
- if (rows_have_padding)
- {
- const unsigned int chunk_row = n_chunk_rows - 1;
-
- const number *const val_end_of_row = &val[cols->sparsity_pattern.rowstart[chunk_row+1]
- * cols->chunk_size
- * cols->chunk_size];
- while (val_ptr != val_end_of_row)
- {
- if ((cols_have_padding == false)
- ||
- (*colnum_ptr != cols->sparsity_pattern.n_cols()-1))
- {
- // we're at a chunk row but not
- // column that has padding
- for (unsigned int r=0; r<m() % cols->chunk_size; ++r)
- for (unsigned int c=0; c<cols->chunk_size; ++c)
- dst(chunk_row * cols->chunk_size + r)
- += (val_ptr[r*cols->chunk_size + c] *
- src(*colnum_ptr * cols->chunk_size + c));
- }
- else
- // we're at a chunk row and
- // column that has padding
- for (unsigned int r=0; r<m() % cols->chunk_size; ++r)
- for (unsigned int c=0; c<n() % cols->chunk_size; ++c)
- dst(chunk_row * cols->chunk_size + r)
- += (val_ptr[r*cols->chunk_size + c] *
- src(*colnum_ptr * cols->chunk_size + c));
-
- ++colnum_ptr;
- val_ptr += cols->chunk_size * cols->chunk_size;
- }
- }
}
const unsigned int n_chunk_rows = cols->sparsity_pattern.n_rows();
- // loop over all chunks. note that we need
- // to treat the last chunk row and column
- // differently if they have padding
- // elements
+ // loop over all chunks. note that we need to treat the last chunk row and
+ // column differently if they have padding elements
const bool rows_have_padding = (m() % cols->chunk_size != 0),
cols_have_padding = (n() % cols->chunk_size != 0);
n_chunk_rows-1 :
n_chunk_rows);
- // like in vmult_add, but don't keep an
- // iterator into dst around since we're not
- // traversing it sequentially this time
+ // like in vmult_add, but don't keep an iterator into dst around since we're
+ // not traversing it sequentially this time
const number *val_ptr = val;
const unsigned int *colnum_ptr = cols->sparsity_pattern.colnums;
src.begin() + chunk_row * cols->chunk_size,
dst.begin() + *colnum_ptr * cols->chunk_size);
else
- // we're at a chunk column that
- // has padding
+ // we're at a chunk column that has padding
for (unsigned int r=0; r<cols->chunk_size; ++r)
for (unsigned int c=0; c<n() % cols->chunk_size; ++c)
dst(*colnum_ptr * cols->chunk_size + c)
}
}
- // now deal with last chunk row if
- // necessary
+ // now deal with last chunk row if necessary
if (rows_have_padding)
{
const unsigned int chunk_row = n_chunk_rows - 1;
||
(*colnum_ptr != cols->sparsity_pattern.n_cols()-1))
{
- // we're at a chunk row but not
- // column that has padding
+ // we're at a chunk row but not column that has padding
for (unsigned int r=0; r<m() % cols->chunk_size; ++r)
for (unsigned int c=0; c<cols->chunk_size; ++c)
dst(*colnum_ptr * cols->chunk_size + c)
src(chunk_row * cols->chunk_size + r));
}
else
- // we're at a chunk row and
- // column that has padding
+ // we're at a chunk row and column that has padding
for (unsigned int r=0; r<m() % cols->chunk_size; ++r)
for (unsigned int c=0; c<n() % cols->chunk_size; ++c)
dst(*colnum_ptr * cols->chunk_size + c)
somenumber result = 0;
////////////////
- // like matrix_scalar_product, except that
- // the two vectors are now the same
+ // like matrix_scalar_product, except that the two vectors are now the same
const unsigned int n_chunk_rows = cols->sparsity_pattern.n_rows();
- // loop over all chunks. note that we need
- // to treat the last chunk row and column
- // differently if they have padding
- // elements
+ // loop over all chunks. note that we need to treat the last chunk row and
+ // column differently if they have padding elements
const bool rows_have_padding = (m() % cols->chunk_size != 0),
cols_have_padding = (n() % cols->chunk_size != 0);
v_ptr,
v.begin() + *colnum_ptr * cols->chunk_size);
else
- // we're at a chunk column that
- // has padding
+ // we're at a chunk column that has padding
for (unsigned int r=0; r<cols->chunk_size; ++r)
for (unsigned int c=0; c<n() % cols->chunk_size; ++c)
result
v_ptr += cols->chunk_size;
}
- // now deal with last chunk row if
- // necessary
+ // now deal with last chunk row if necessary
if (rows_have_padding)
{
const unsigned int chunk_row = n_chunk_rows - 1;
||
(*colnum_ptr != cols->sparsity_pattern.n_cols()-1))
{
- // we're at a chunk row but not
- // column that has padding
+ // we're at a chunk row but not column that has padding
for (unsigned int r=0; r<m() % cols->chunk_size; ++r)
for (unsigned int c=0; c<cols->chunk_size; ++c)
result
v(*colnum_ptr * cols->chunk_size + c));
}
else
- // we're at a chunk row and
- // column that has padding
+ // we're at a chunk row and column that has padding
for (unsigned int r=0; r<m() % cols->chunk_size; ++r)
for (unsigned int c=0; c<n() % cols->chunk_size; ++c)
result
Assert(m() == u.size(), ExcDimensionMismatch(m(),u.size()));
Assert(n() == v.size(), ExcDimensionMismatch(n(),v.size()));
- // the following works like the vmult_add
- // function
+ // the following works like the vmult_add function
somenumber result = 0;
const unsigned int n_chunk_rows = cols->sparsity_pattern.n_rows();
- // loop over all chunks. note that we need
- // to treat the last chunk row and column
- // differently if they have padding
- // elements
+ // loop over all chunks. note that we need to treat the last chunk row and
+ // column differently if they have padding elements
const bool rows_have_padding = (m() % cols->chunk_size != 0),
cols_have_padding = (n() % cols->chunk_size != 0);
u_ptr,
v.begin() + *colnum_ptr * cols->chunk_size);
else
- // we're at a chunk column that
- // has padding
+ // we're at a chunk column that has padding
for (unsigned int r=0; r<cols->chunk_size; ++r)
for (unsigned int c=0; c<n() % cols->chunk_size; ++c)
result
u_ptr += cols->chunk_size;
}
- // now deal with last chunk row if
- // necessary
+ // now deal with last chunk row if necessary
if (rows_have_padding)
{
const unsigned int chunk_row = n_chunk_rows - 1;
||
(*colnum_ptr != cols->sparsity_pattern.n_cols()-1))
{
- // we're at a chunk row but not
- // column that has padding
+ // we're at a chunk row but not column that has padding
for (unsigned int r=0; r<m() % cols->chunk_size; ++r)
for (unsigned int c=0; c<cols->chunk_size; ++c)
result
v(*colnum_ptr * cols->chunk_size + c));
}
else
- // we're at a chunk row and
- // column that has padding
+ // we're at a chunk row and column that has padding
for (unsigned int r=0; r<m() % cols->chunk_size; ++r)
for (unsigned int c=0; c<n() % cols->chunk_size; ++c)
result
const unsigned int n_chunk_rows = cols->sparsity_pattern.n_rows();
- // loop over all rows and columns; it is
- // safe to also loop over the padding
- // elements (they are zero) if we make sure
- // that the vector into which we sum column
- // sums is large enough
+ // loop over all rows and columns; it is safe to also loop over the padding
+ // elements (they are zero) if we make sure that the vector into which we
+ // sum column sums is large enough
Vector<real_type> column_sums(cols->sparsity_pattern.n_cols() *
cols->chunk_size);
Assert (cols != 0, ExcNotInitialized());
Assert (val != 0, ExcNotInitialized());
- // this function works like l1_norm(). it
- // can be made more efficient (without
- // allocating a temporary vector) as is
- // done in the SparseMatrix class but since
- // it is rarely called in time critical
- // places it is probably not worth it
+ // this function works like l1_norm(). it can be made more efficient
+ // (without allocating a temporary vector) as is done in the SparseMatrix
+ // class but since it is rarely called in time critical places it is
+ // probably not worth it
const unsigned int n_chunk_rows = cols->sparsity_pattern.n_rows();
- // loop over all rows and columns; it is
- // safe to also loop over the padding
- // elements (they are zero) if we make sure
- // that the vector into which we sum column
- // sums is large enough
+ // loop over all rows and columns; it is safe to also loop over the padding
+ // elements (they are zero) if we make sure that the vector into which we
+ // sum column sums is large enough
Vector<real_type> row_sums(cols->sparsity_pattern.n_rows() *
cols->chunk_size);
typename ChunkSparseMatrix<number>::real_type
ChunkSparseMatrix<number>::frobenius_norm () const
{
- // simply add up all entries in the
- // sparsity pattern, without taking any
+ // simply add up all entries in the sparsity pattern, without taking any
// reference to rows or columns
//
- // padding elements are zero, so we can add
- // them up as well
+ // padding elements are zero, so we can add them up as well
real_type norm_sqr = 0;
for (const number *ptr = &val[0]; ptr != &val[max_len]; ++ptr)
norm_sqr += numbers::NumberTraits<number>::abs_square(*ptr);
Assert (&u != &dst, ExcSourceEqualsDestination());
- // set dst=b, then subtract the result of
- // A*u from it. since the purpose of the
- // current class is to promote streaming of
- // data rather than more random access
- // patterns, breaking things up into two
- // loops may be reasonable
+ // set dst=b, then subtract the result of A*u from it. since the purpose of
+ // the current class is to promote streaming of data rather than more random
+ // access patterns, breaking things up into two loops may be reasonable
dst = b;
/////////
- // the rest of this function is like
- // vmult_add, except that we subtract
+ // the rest of this function is like vmult_add, except that we subtract
// rather than add A*u
/////////
const unsigned int n_chunk_rows = cols->sparsity_pattern.n_rows();
- // loop over all chunks. note that we need
- // to treat the last chunk row and column
- // differently if they have padding
- // elements
+ // loop over all chunks. note that we need to treat the last chunk row and
+ // column differently if they have padding elements
const bool rows_have_padding = (m() % cols->chunk_size != 0),
cols_have_padding = (n() % cols->chunk_size != 0);
u.begin() + *colnum_ptr * cols->chunk_size,
dst_ptr);
else
- // we're at a chunk column that
- // has padding
+ // we're at a chunk column that has padding
for (unsigned int r=0; r<cols->chunk_size; ++r)
for (unsigned int c=0; c<n() % cols->chunk_size; ++c)
dst(chunk_row * cols->chunk_size + r)
dst_ptr += cols->chunk_size;
}
- // now deal with last chunk row if
- // necessary
+ // now deal with last chunk row if necessary
if (rows_have_padding)
{
const unsigned int chunk_row = n_chunk_rows - 1;
||
(*colnum_ptr != cols->sparsity_pattern.n_cols()-1))
{
- // we're at a chunk row but not
- // column that has padding
+ // we're at a chunk row but not column that has padding
for (unsigned int r=0; r<m() % cols->chunk_size; ++r)
for (unsigned int c=0; c<cols->chunk_size; ++c)
dst(chunk_row * cols->chunk_size + r)
u(*colnum_ptr * cols->chunk_size + c));
}
else
- // we're at a chunk row and
- // column that has padding
+ // we're at a chunk row and column that has padding
for (unsigned int r=0; r<m() % cols->chunk_size; ++r)
for (unsigned int c=0; c<n() % cols->chunk_size; ++c)
dst(chunk_row * cols->chunk_size + r)
const Vector<somenumber> &src,
const number /*om*/) const
{
- // to understand how this function works
- // you may want to take a look at the CVS
- // archives to see the original version
- // which is much clearer...
+ // to understand how this function works you may want to take a look at the
+ // CVS archives to see the original version which is much clearer...
Assert (cols != 0, ExcNotInitialized());
Assert (val != 0, ExcNotInitialized());
Assert (m() == n(), ExcMessage("This operation is only valid on square matrices."));
-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
{
const unsigned int chunk_size = cols->get_chunk_size();
- // loop over all chunk rows and columns,
- // and each time we find something repeat
- // it chunk_size times in both directions
+ // loop over all chunk rows and columns, and each time we find something
+ // repeat it chunk_size times in both directions
for (unsigned int i=0; i<cols->sparsity_pattern.n_rows(); ++i)
{
for (unsigned int d=0; d<chunk_size; ++d)
{
AssertThrow (out, ExcIO());
- // first the simple objects,
- // bracketed in [...]
+ // first the simple objects, bracketed in [...]
out << '[' << max_len << "][";
// then write out real data
out.write (reinterpret_cast<const char *>(&val[0]),
*/
+
+/**
+ * Iterators on sparsity patterns
+ */
+namespace ChunkSparsityPatternIterators
+{
+ // forward declaration
+ class Iterator;
+
+ /**
+ * Accessor class for iterators into sparsity patterns. This class is
+ * also the base class for both const and non-const accessor classes
+ * into sparse matrices.
+ *
+ * Note that this class only allows read access to elements, providing
+ * their row and column number. It does not allow modifying the
+ * sparsity pattern itself.
+ *
+ * @author Martin Kronbichler
+ * @date 2013
+ */
+ class Accessor
+ {
+ public:
+ /**
+ * Constructor.
+ */
+ Accessor (const ChunkSparsityPattern *matrix,
+ const unsigned int row);
+
+ /**
+ * Constructor. Construct the end accessor for the given sparsity pattern.
+ */
+ Accessor (const ChunkSparsityPattern *matrix);
+
+ /**
+ * Row number of the element represented by this object. This function can
+ * only be called for entries for which is_valid_entry() is true.
+ */
+ unsigned int row () const;
+
+ /**
+ * Index in row of the element represented by this object. This function
+ * can only be called for entries for which is_valid_entry() is true.
+ */
+ unsigned int index () const;
+
+ /**
+ * Returns the global index from the reduced sparsity pattern.
+ */
+ std::size_t reduced_index() const;
+
+ /**
+ * Column number of the element represented by this object. This function
+ * can only be called for entries for which is_valid_entry() is true.
+ */
+ unsigned int column () const;
+
+ /**
+ * Return whether the sparsity pattern entry pointed to by this iterator
+ * is valid or not. Note that after compressing the sparsity pattern, all
+ * entries are valid. However, before compression, the sparsity pattern
+ * allocated some memory to be used while still adding new nonzero
+ * entries; if you create iterators in this phase of the sparsity
+ * pattern's lifetime, you will iterate over elements that are not
+ * valid. If this is so, then this function will return false.
+ */
+ bool is_valid_entry () const;
+
+
+ /**
+ * Comparison. True, if both iterators point to the same matrix position.
+ */
+ bool operator == (const Accessor &) const;
+
+
+ /**
+ * Comparison operator. Result is true if either the first row number is
+ * smaller or if the row numbers are equal and the first index is smaller.
+ *
+ * This function is only valid if both iterators point into the same
+ * sparsity pattern.
+ */
+ bool operator < (const Accessor &) const;
+
+ protected:
+ /**
+ * The sparsity pattern we operate on accessed.
+ */
+ const ChunkSparsityPattern *sparsity_pattern;
+
+ /**
+ * The accessor of the (reduced) sparsity pattern.
+ */
+ SparsityPatternIterators::Accessor reduced_accessor;
+
+ /**
+ * Current row number.
+ */
+ unsigned int chunk_row;
+
+ /**
+ * Current index in row.
+ */
+ unsigned int chunk_col;
+
+ /**
+ * Move the accessor to the next nonzero entry in the matrix.
+ */
+ void advance ();
+
+ /**
+ * Grant access to iterator class.
+ */
+ friend class Iterator;
+ };
+
+
+
+ /**
+ * STL conforming iterator walking over the elements of a sparsity pattern.
+ */
+ class Iterator
+ {
+ public:
+ /**
+ * Constructor. Create an iterator into the sparsity pattern @p sp for the
+ * given row and the index within it.
+ */
+ Iterator (const ChunkSparsityPattern *sp,
+ const unsigned int row);
+
+ /**
+ * Prefix increment.
+ */
+ Iterator &operator++ ();
+
+ /**
+ * Postfix increment.
+ */
+ Iterator operator++ (int);
+
+ /**
+ * Dereferencing operator.
+ */
+ const Accessor &operator* () const;
+
+ /**
+ * Dereferencing operator.
+ */
+ const Accessor *operator-> () const;
+
+ /**
+ * Comparison. True, if both iterators point to the same matrix position.
+ */
+ bool operator == (const Iterator &) const;
+
+ /**
+ * Inverse of <tt>==</tt>.
+ */
+ bool operator != (const Iterator &) const;
+
+ /**
+ * Comparison operator. Result is true if either the first row number is
+ * smaller or if the row numbers are equal and the first index is smaller.
+ *
+ * This function is only valid if both iterators point into the same
+ * matrix.
+ */
+ bool operator < (const Iterator &) const;
+
+ private:
+ /**
+ * Store an object of the accessor class.
+ */
+ Accessor accessor;
+ };
+}
+
+
+
/**
* Structure representing the sparsity pattern of a sparse matrix.
*
class ChunkSparsityPattern : public Subscriptor
{
public:
+ /**
+ * Typedef an iterator class that allows to walk over all nonzero elements
+ * of a sparsity pattern.
+ */
+ typedef ChunkSparsityPatternIterators::Iterator const_iterator;
+
+ /**
+ * Typedef an iterator class that allows to walk over all nonzero elements
+ * of a sparsity pattern.
+ *
+ * Since the iterator does not allow to modify the sparsity pattern, this
+ * type is the same as that for @p const_iterator.
+ */
+ typedef ChunkSparsityPatternIterators::Iterator iterator;
/**
- * Define a value which is used
- * to indicate that a certain
- * value in the colnums array
- * is unused, i.e. does not
- * represent a certain column
- * number index.
+ * Define a value which is used to indicate that a certain value in the
+ * colnums array is unused, i.e. does not represent a certain column number
+ * index.
*
- * Indices with this invalid
- * value are used to insert new
- * entries to the sparsity
- * pattern using the add() member
- * function, and are removed when
+ * Indices with this invalid value are used to insert new entries to the
+ * sparsity pattern using the add() member function, and are removed when
* calling compress().
*
- * You should not assume that the
- * variable declared here has a
- * certain value. The
- * initialization is given here
- * only to enable the compiler to
- * perform some optimizations,
- * but the actual value of the
- * variable may change over time.
+ * You should not assume that the variable declared here has a certain
+ * value. The initialization is given here only to enable the compiler to
+ * perform some optimizations, but the actual value of the variable may
+ * change over time.
*/
static const unsigned int invalid_entry = SparsityPattern::invalid_entry;
/**
- * Initialize the matrix empty,
- * that is with no memory
- * allocated. This is useful if
- * you want such objects as
- * member variables in other
- * classes. You can make the
- * structure usable by calling
- * the reinit() function.
+ * Initialize the matrix empty, that is with no memory allocated. This is
+ * useful if you want such objects as member variables in other classes. You
+ * can make the structure usable by calling the reinit() function.
*/
ChunkSparsityPattern ();
/**
- * Copy constructor. This
- * constructor is only allowed to
- * be called if the matrix
- * structure to be copied is
- * empty. This is so in order to
- * prevent involuntary copies of
- * objects for temporaries, which
- * can use large amounts of
- * computing time. However, copy
- * constructors are needed if yo
- * want to use the STL data types
- * on classes like this, e.g. to
- * write such statements like
- * <tt>v.push_back
- * (ChunkSparsityPattern());</tt>,
- * with <tt>v</tt> a vector of
- * ChunkSparsityPattern objects.
- *
- * Usually, it is sufficient to
- * use the explicit keyword to
- * disallow unwanted temporaries,
- * but for the STL vectors, this
- * does not work. Since copying a
- * structure like this is not
- * useful anyway because multiple
- * matrices can use the same
- * sparsity structure, copies are
- * only allowed for empty
- * objects, as described above.
+ * Copy constructor. This constructor is only allowed to be called if the
+ * matrix structure to be copied is empty. This is so in order to prevent
+ * involuntary copies of objects for temporaries, which can use large
+ * amounts of computing time. However, copy constructors are needed if yo
+ * want to use the STL data types on classes like this, e.g. to write such
+ * statements like <tt>v.push_back (ChunkSparsityPattern());</tt>, with
+ * <tt>v</tt> a vector of ChunkSparsityPattern objects.
+ *
+ * Usually, it is sufficient to use the explicit keyword to disallow
+ * unwanted temporaries, but for the STL vectors, this does not work. Since
+ * copying a structure like this is not useful anyway because multiple
+ * matrices can use the same sparsity structure, copies are only allowed for
+ * empty objects, as described above.
*/
ChunkSparsityPattern (const ChunkSparsityPattern &);
/**
- * Initialize a rectangular
- * matrix.
+ * Initialize a rectangular matrix.
*
* @arg m number of rows
* @arg n number of columns
- * @arg max_per_row maximum
- * number of nonzero entries per row
+ * @arg max_per_row maximum number of nonzero entries per row
*/
ChunkSparsityPattern (const unsigned int m,
const unsigned int n,
const bool optimize_diagonal) DEAL_II_DEPRECATED;
/**
- * Initialize a rectangular
- * matrix.
+ * Initialize a rectangular matrix.
*
* @arg m number of rows
* @arg n number of columns
- *
- * @arg row_lengths possible
- * number of nonzero entries for
- * each row. This vector must
- * have one entry for each row.
+ * @arg row_lengths possible number of nonzero entries for each row. This
+ * vector must have one entry for each row.
*/
ChunkSparsityPattern (const unsigned int m,
const unsigned int n,
const bool optimize_diagonal) DEAL_II_DEPRECATED;
/**
- * Initialize a quadratic matrix
- * of dimension <tt>n</tt> with
- * at most <tt>max_per_row</tt>
- * nonzero entries per row.
+ * Initialize a quadratic matrix of dimension <tt>n</tt> with at most
+ * <tt>max_per_row</tt> nonzero entries per row.
*
- * This constructor automatically
- * enables optimized storage of
- * diagonal elements. To avoid
- * this, use the constructor
- * taking row and column numbers
- * separately.
+ * This constructor automatically enables optimized storage of diagonal
+ * elements. To avoid this, use the constructor taking row and column
+ * numbers separately.
*/
ChunkSparsityPattern (const unsigned int n,
const unsigned int max_per_row,
* Initialize a quadratic matrix.
*
* @arg m number of rows and columns
- *
- * @arg row_lengths possible
- * number of nonzero entries for
- * each row. This vector must
- * have one entry for each row.
+ * @arg row_lengths possible number of nonzero entries for each row. This
+ * vector must have one entry for each row.
*/
ChunkSparsityPattern (const unsigned int m,
const std::vector<unsigned int> &row_lengths,
~ChunkSparsityPattern ();
/**
- * Copy operator. For this the
- * same holds as for the copy
- * constructor: it is declared,
- * defined and fine to be called,
- * but the latter only for empty
+ * Copy operator. For this the same holds as for the copy constructor: it is
+ * declared, defined and fine to be called, but the latter only for empty
* objects.
*/
ChunkSparsityPattern &operator = (const ChunkSparsityPattern &);
/**
- * Reallocate memory and set up data
- * structures for a new matrix with
- * <tt>m </tt>rows and <tt>n</tt> columns,
- * with at most <tt>max_per_row</tt>
+ * Reallocate memory and set up data structures for a new matrix with <tt>m
+ * </tt>rows and <tt>n</tt> columns, with at most <tt>max_per_row</tt>
* nonzero entries per row.
*
- * This function simply maps its
- * operations to the other
- * <tt>reinit</tt> function.
+ * This function simply maps its operations to the other <tt>reinit</tt>
+ * function.
*/
void reinit (const unsigned int m,
const unsigned int n,
const bool optimize_diagonal) DEAL_II_DEPRECATED;
/**
- * Reallocate memory for a matrix
- * of size <tt>m x n</tt>. The
- * number of entries for each row
- * is taken from the array
- * <tt>row_lengths</tt> which has to
- * give this number of each row
- * <tt>i=1...m</tt>.
+ * Reallocate memory for a matrix of size <tt>m x n</tt>. The number of
+ * entries for each row is taken from the array <tt>row_lengths</tt> which
+ * has to give this number of each row <tt>i=1...m</tt>.
*
- * If <tt>m*n==0</tt> all memory is freed,
- * resulting in a total reinitialization
- * of the object. If it is nonzero, new
- * memory is only allocated if the new
- * size extends the old one. This is done
- * to save time and to avoid fragmentation
- * of the heap.
+ * If <tt>m*n==0</tt> all memory is freed, resulting in a total
+ * reinitialization of the object. If it is nonzero, new memory is only
+ * allocated if the new size extends the old one. This is done to save time
+ * and to avoid fragmentation of the heap.
*
* If the number of rows equals
* the number of columns then
const bool optimize_diagonal) DEAL_II_DEPRECATED;
/**
- * Same as above, but with a
- * VectorSlice argument instead.
+ * Same as above, but with a VectorSlice argument instead.
*/
void reinit (const unsigned int m,
const unsigned int n,
const bool optimize_diagonal) DEAL_II_DEPRECATED;
/**
- * This function compresses the sparsity
- * structure that this object represents.
- * It does so by eliminating unused
- * entries and sorting the remaining ones
- * to allow faster access by usage of
- * binary search algorithms. A special
- * sorting scheme is used for the
- * diagonal entry of quadratic matrices,
- * which is always the first entry of
- * each row.
+ * This function compresses the sparsity structure that this object
+ * represents. It does so by eliminating unused entries and sorting the
+ * remaining ones to allow faster access by usage of binary search
+ * algorithms. A special sorting scheme is used for the diagonal entry of
+ * quadratic matrices, which is always the first entry of each row.
*
- * The memory which is no more
- * needed is released.
+ * The memory which is no more needed is released.
*
- * SparseMatrix objects require the
- * ChunkSparsityPattern objects they are
- * initialized with to be compressed, to
- * reduce memory requirements.
+ * SparseMatrix objects require the ChunkSparsityPattern objects they are
+ * initialized with to be compressed, to reduce memory requirements.
*/
void compress ();
/**
- * This function can be used as a
- * replacement for reinit(),
- * subsequent calls to add() and
- * a final call to close() if you
- * know exactly in advance the
- * entries that will form the
- * matrix sparsity pattern.
- *
- * The first two parameters
- * determine the size of the
- * matrix. For the two last ones,
- * note that a sparse matrix can
- * be described by a sequence of
- * rows, each of which is
- * represented by a sequence of
- * pairs of column indices and
- * values. In the present
- * context, the begin() and
- * end() parameters designate
- * iterators (of forward iterator
- * type) into a container, one
- * representing one row. The
- * distance between begin()
- * and end() should therefore
- * be equal to
- * n_rows(). These iterators
- * may be iterators of
- * <tt>std::vector</tt>,
- * <tt>std::list</tt>, pointers into a
- * C-style array, or any other
- * iterator satisfying the
- * requirements of a forward
- * iterator. The objects pointed
- * to by these iterators
- * (i.e. what we get after
- * applying <tt>operator*</tt> or
- * <tt>operator-></tt> to one of these
- * iterators) must be a container
- * itself that provides functions
- * <tt>begin</tt> and <tt>end</tt>
- * designating a range of
- * iterators that describe the
- * contents of one
- * line. Dereferencing these
- * inner iterators must either
- * yield a pair of an unsigned
- * integer as column index and a
- * value of arbitrary type (such
- * a type would be used if we
- * wanted to describe a sparse
- * matrix with one such object),
- * or simply an unsigned integer
- * (of we only wanted to describe
- * a sparsity pattern). The
- * function is able to determine
- * itself whether an unsigned
- * integer or a pair is what we
- * get after dereferencing the
- * inner iterators, through some
- * template magic.
- *
- * While the order of the outer
- * iterators denotes the
- * different rows of the matrix,
- * the order of the inner
- * iterator denoting the columns
- * does not matter, as they are
- * sorted internal to this
- * function anyway.
- *
- * Since that all sounds very
- * complicated, consider the
- * following example code, which
- * may be used to fill a sparsity
- * pattern:
+ * This function can be used as a replacement for reinit(), subsequent calls
+ * to add() and a final call to close() if you know exactly in advance the
+ * entries that will form the matrix sparsity pattern.
+ *
+ * The first two parameters determine the size of the matrix. For the two
+ * last ones, note that a sparse matrix can be described by a sequence of
+ * rows, each of which is represented by a sequence of pairs of column
+ * indices and values. In the present context, the begin() and end()
+ * parameters designate iterators (of forward iterator type) into a
+ * container, one representing one row. The distance between begin() and
+ * end() should therefore be equal to n_rows(). These iterators may be
+ * iterators of <tt>std::vector</tt>, <tt>std::list</tt>, pointers into a
+ * C-style array, or any other iterator satisfying the requirements of a
+ * forward iterator. The objects pointed to by these iterators (i.e. what we
+ * get after applying <tt>operator*</tt> or <tt>operator-></tt> to one of
+ * these iterators) must be a container itself that provides functions
+ * <tt>begin</tt> and <tt>end</tt> designating a range of iterators that
+ * describe the contents of one line. Dereferencing these inner iterators
+ * must either yield a pair of an unsigned integer as column index and a
+ * value of arbitrary type (such a type would be used if we wanted to
+ * describe a sparse matrix with one such object), or simply an unsigned
+ * integer (of we only wanted to describe a sparsity pattern). The function
+ * is able to determine itself whether an unsigned integer or a pair is what
+ * we get after dereferencing the inner iterators, through some template
+ * magic.
+ *
+ * While the order of the outer iterators denotes the different rows of the
+ * matrix, the order of the inner iterator denoting the columns does not
+ * matter, as they are sorted internal to this function anyway.
+ *
+ * Since that all sounds very complicated, consider the following example
+ * code, which may be used to fill a sparsity pattern:
* @code
* std::vector<std::vector<unsigned int> > column_indices (n_rows);
* for (unsigned int row=0; row<n_rows; ++row)
* column_indices.end());
* @endcode
*
- * Note that this example works
- * since the iterators
- * dereferenced yield containers
- * with functions <tt>begin</tt> and
- * <tt>end</tt> (namely
- * <tt>std::vector</tt>s), and the
- * inner iterators dereferenced
- * yield unsigned integers as
- * column indices. Note that we
- * could have replaced each of
- * the two <tt>std::vector</tt>
- * occurrences by <tt>std::list</tt>,
- * and the inner one by
- * <tt>std::set</tt> as well.
- *
- * Another example would be as
- * follows, where we initialize a
- * whole matrix, not only a
- * sparsity pattern:
+ * Note that this example works since the iterators dereferenced yield
+ * containers with functions <tt>begin</tt> and <tt>end</tt> (namely
+ * <tt>std::vector</tt>s), and the inner iterators dereferenced yield
+ * unsigned integers as column indices. Note that we could have replaced
+ * each of the two <tt>std::vector</tt> occurrences by <tt>std::list</tt>,
+ * and the inner one by <tt>std::set</tt> as well.
+ *
+ * Another example would be as follows, where we initialize a whole matrix,
+ * not only a sparsity pattern:
* @code
* std::vector<std::map<unsigned int,double> > entries (n_rows);
* for (unsigned int row=0; row<n_rows; ++row)
* column_indices.end());
* @endcode
*
- * This example works because
- * dereferencing iterators of the
- * inner type yields a pair of
- * unsigned integers and a value,
- * the first of which we take as
- * column index. As previously,
- * the outer <tt>std::vector</tt>
- * could be replaced by
- * <tt>std::list</tt>, and the inner
- * <tt>std::map<unsigned int,double></tt>
- * could be replaced by
- * <tt>std::vector<std::pair<unsigned int,double> ></tt>,
- * or a list or set of such
- * pairs, as they all return
- * iterators that point to such
- * pairs.
+ * This example works because dereferencing iterators of the inner type
+ * yields a pair of unsigned integers and a value, the first of which we
+ * take as column index. As previously, the outer <tt>std::vector</tt> could
+ * be replaced by <tt>std::list</tt>, and the inner <tt>std::map<unsigned
+ * int,double></tt> could be replaced by <tt>std::vector<std::pair<unsigned
+ * int,double> ></tt>, or a list or set of such pairs, as they all return
+ * iterators that point to such pairs.
*/
template <typename ForwardIterator>
void copy_from (const unsigned int n_rows,
const bool optimize_diagonal) DEAL_II_DEPRECATED;
/**
- * Copy data from an object of type
- * CompressedSparsityPattern,
- * CompressedSetSparsityPattern or
- * CompressedSimpleSparsityPattern.
- * Previous content of this object is
- * lost, and the sparsity pattern is in
+ * Copy data from an object of type CompressedSparsityPattern,
+ * CompressedSetSparsityPattern or CompressedSimpleSparsityPattern.
+ * Previous content of this object is lost, and the sparsity pattern is in
* compressed mode afterwards.
*/
template <typename SparsityType>
const bool optimize_diagonal) DEAL_II_DEPRECATED;
/**
- * Take a full matrix and use its
- * nonzero entries to generate a
- * sparse matrix entry pattern
- * for this object.
+ * Take a full matrix and use its nonzero entries to generate a sparse
+ * matrix entry pattern for this object.
*
- * Previous content of this
- * object is lost, and the
- * sparsity pattern is in
+ * Previous content of this object is lost, and the sparsity pattern is in
* compressed mode afterwards.
*/
template <typename number>
const bool optimize_diagonal) DEAL_II_DEPRECATED;
/**
- * Return whether the object is empty. It
- * is empty if no memory is allocated,
- * which is the same as that both
- * dimensions are zero.
+ * Set the sparsity pattern of the chunk sparsity pattern to be given by
+ * <tt>chunk_size*chunksize</tt> blocks of the sparsity pattern for chunks
+ * specified. Note that the final number of rows <tt>m</tt> of the sparsity
+ * pattern will be approximately <tt>sparsity_pattern_for_chunks.n_rows() *
+ * chunk_size</tt> (modulo padding elements in the last chunk) and similarly
+ * for the number of columns <tt>n</tt>.
+ *
+ * This is a special initialization option in case you can tell the position
+ * of the chunk already from the beginning without generating the sparsity
+ * pattern using <tt>make_sparsity_pattern</tt> calls. This bypasses the
+ * search for chunks but of course needs to be handled with care in order to
+ * give a correct sparsity pattern.
+ *
+ * Previous content of this object is lost, and the sparsity pattern is in
+ * compressed mode afterwards.
+ */
+ template <typename Sparsity>
+ void create_from (const unsigned int m,
+ const unsigned int n,
+ const Sparsity &sparsity_pattern_for_chunks,
+ const unsigned int chunk_size,
+ const bool optimize_diagonal = true);
+
+ /**
+ * Return whether the object is empty. It is empty if no memory is
+ * allocated, which is the same as that both dimensions are zero.
*/
bool empty () const;
/**
- * Return the chunk size given as
- * argument when constructing this
- * object.
+ * Return the chunk size given as argument when constructing this object.
*/
unsigned int get_chunk_size () const;
/**
- * Return the maximum number of entries per
- * row. Before compression, this equals the
- * number given to the constructor, while
- * after compression, it equals the maximum
- * number of entries actually allocated by
- * the user.
+ * Return the maximum number of entries per row. Before compression, this
+ * equals the number given to the constructor, while after compression, it
+ * equals the maximum number of entries actually allocated by the user.
*/
unsigned int max_entries_per_row () const;
/**
- * Add a nonzero entry to the matrix.
- * This function may only be called
- * for non-compressed sparsity patterns.
+ * Add a nonzero entry to the matrix. This function may only be called for
+ * non-compressed sparsity patterns.
*
- * If the entry already exists, nothing
- * bad happens.
+ * If the entry already exists, nothing bad happens.
*/
void add (const unsigned int i,
const unsigned int j);
/**
- * Make the sparsity pattern
- * symmetric by adding the
- * sparsity pattern of the
+ * Make the sparsity pattern symmetric by adding the sparsity pattern of the
* transpose object.
*
- * This function throws an
- * exception if the sparsity
- * pattern does not represent a
- * quadratic matrix.
+ * This function throws an exception if the sparsity pattern does not
+ * represent a quadratic matrix.
*/
void symmetrize ();
/**
- * Return number of rows of this
- * matrix, which equals the dimension
- * of the image space.
+ * Return number of rows of this matrix, which equals the dimension of the
+ * image space.
*/
- inline unsigned int n_rows () const;
+ unsigned int n_rows () const;
/**
- * Return number of columns of this
- * matrix, which equals the dimension
- * of the range space.
+ * Return number of columns of this matrix, which equals the dimension of
+ * the range space.
*/
- inline unsigned int n_cols () const;
+ unsigned int n_cols () const;
/**
- * Check if a value at a certain
- * position may be non-zero.
+ * Check if a value at a certain position may be non-zero.
*/
bool exists (const unsigned int i,
const unsigned int j) const;
unsigned int row_length (const unsigned int row) const;
/**
- * Compute the bandwidth of the matrix
- * represented by this structure. The
- * bandwidth is the maximum of $|i-j|$
- * for which the index pair $(i,j)$
- * represents a nonzero entry of the
- * matrix. Consequently, the maximum
- * bandwidth a $n\times m$ matrix can
- * have is $\max\{n-1,m-1\}$.
+ * Compute the bandwidth of the matrix represented by this structure. The
+ * bandwidth is the maximum of $|i-j|$ for which the index pair $(i,j)$
+ * represents a nonzero entry of the matrix. Consequently, the maximum
+ * bandwidth a $n\times m$ matrix can have is $\max\{n-1,m-1\}$.
*/
unsigned int bandwidth () const;
/**
- * Return the number of nonzero elements of
- * this matrix. Actually, it returns the
- * number of entries in the sparsity
- * pattern; if any of the entries should
- * happen to be zero, it is counted
- * anyway.
+ * Return the number of nonzero elements of this matrix. Actually, it
+ * returns the number of entries in the sparsity pattern; if any of the
+ * entries should happen to be zero, it is counted anyway.
*
- * This function may only be called if the
- * matrix struct is compressed. It does not
- * make too much sense otherwise anyway.
+ * This function may only be called if the matrix struct is compressed. It
+ * does not make too much sense otherwise anyway.
*/
unsigned int n_nonzero_elements () const;
/**
- * Return whether the structure is
- * compressed or not.
+ * Return whether the structure is compressed or not.
*/
bool is_compressed () const;
bool stores_only_added_elements () const;
/**
- * Write the data of this object
- * en bloc to a file. This is
- * done in a binary mode, so the
- * output is neither readable by
- * humans nor (probably) by other
- * computers using a different
- * operating system of number
- * format.
- *
- * The purpose of this function
- * is that you can swap out
- * matrices and sparsity pattern
- * if you are short of memory,
- * want to communicate between
- * different programs, or allow
- * objects to be persistent
- * across different runs of the
- * program.
+ * STL-like iterator with the first entry of the matrix. The resulting
+ * iterator can be used to walk over all nonzero entries of the sparsity
+ * pattern.
+ */
+ iterator begin () const;
+
+ /**
+ * Final iterator.
+ */
+ iterator end () const;
+
+ /**
+ * STL-like iterator with the first entry of row <tt>r</tt>.
+ *
+ * Note that if the given row is empty, i.e. does not contain any nonzero
+ * entries, then the iterator returned by this function equals
+ * <tt>end(r)</tt>. Note also that the iterator may not be dereferencable in
+ * that case.
+ */
+ iterator begin (const unsigned int r) const;
+
+ /**
+ * Final iterator of row <tt>r</tt>. It points to the first element past the
+ * end of line @p r, or past the end of the entire sparsity pattern.
+ *
+ * Note that the end iterator is not necessarily dereferencable. This is in
+ * particular the case if it is the end iterator for the last row of a
+ * matrix.
+ */
+ iterator end (const unsigned int r) const;
+
+ /**
+ * Write the data of this object en bloc to a file. This is done in a binary
+ * mode, so the output is neither readable by humans nor (probably) by other
+ * computers using a different operating system of number format.
+ *
+ * The purpose of this function is that you can swap out matrices and
+ * sparsity pattern if you are short of memory, want to communicate between
+ * different programs, or allow objects to be persistent across different
+ * runs of the program.
*/
void block_write (std::ostream &out) const;
/**
- * Read data that has previously
- * been written by block_write()
- * from a file. This is done
- * using the inverse operations
- * to the above function, so it
- * is reasonably fast because the
- * bitstream is not interpreted
- * except for a few numbers up
- * front.
- *
- * The object is resized on this
- * operation, and all previous
- * contents are lost.
- *
- * A primitive form of error
- * checking is performed which
- * will recognize the bluntest
- * attempts to interpret some
- * data as a vector stored
- * bitwise to a file, but not
- * more.
+ * Read data that has previously been written by block_write() from a
+ * file. This is done using the inverse operations to the above function, so
+ * it is reasonably fast because the bitstream is not interpreted except for
+ * a few numbers up front.
+ *
+ * The object is resized on this operation, and all previous contents are
+ * lost.
+ *
+ * A primitive form of error checking is performed which will recognize the
+ * bluntest attempts to interpret some data as a vector stored bitwise to a
+ * file, but not more.
*/
void block_read (std::istream &in);
/**
- * Print the sparsity of the
- * matrix. The output consists of
- * one line per row of the format
- * <tt>[i,j1,j2,j3,...]</tt>. <i>i</i>
- * is the row number and
- * <i>jn</i> are the allocated
- * columns in this row.
+ * Print the sparsity of the matrix. The output consists of one line per row
+ * of the format <tt>[i,j1,j2,j3,...]</tt>. <i>i</i> is the row number and
+ * <i>jn</i> are the allocated columns in this row.
*/
void print (std::ostream &out) const;
/**
- * Print the sparsity of the matrix
- * in a format that <tt>gnuplot</tt> understands
- * and which can be used to plot the
- * sparsity pattern in a graphical
- * way. The format consists of pairs
- * <tt>i j</tt> of nonzero elements, each
- * representing one entry of this
- * matrix, one per line of the output
- * file. Indices are counted from
- * zero on, as usual. Since sparsity
- * patterns are printed in the same
- * way as matrices are displayed, we
- * print the negative of the column
- * index, which means that the
- * <tt>(0,0)</tt> element is in the top left
- * rather than in the bottom left
- * corner.
- *
- * Print the sparsity pattern in
- * gnuplot by setting the data style
- * to dots or points and use the
- * <tt>plot</tt> command.
+ * Print the sparsity of the matrix in a format that <tt>gnuplot</tt>
+ * understands and which can be used to plot the sparsity pattern in a
+ * graphical way. The format consists of pairs <tt>i j</tt> of nonzero
+ * elements, each representing one entry of this matrix, one per line of the
+ * output file. Indices are counted from zero on, as usual. Since sparsity
+ * patterns are printed in the same way as matrices are displayed, we print
+ * the negative of the column index, which means that the <tt>(0,0)</tt>
+ * element is in the top left rather than in the bottom left corner.
+ *
+ * Print the sparsity pattern in gnuplot by setting the data style to dots
+ * or points and use the <tt>plot</tt> command.
*/
void print_gnuplot (std::ostream &out) const;
/**
- * Determine an estimate for the
- * memory consumption (in bytes)
- * of this object. See
- * MemoryConsumption.
+ * Determine an estimate for the memory consumption (in bytes) of this
+ * object. See MemoryConsumption.
*/
std::size_t memory_consumption () const;
//@}
private:
/**
- * Number of rows that this sparsity
- * structure shall represent.
+ * Number of rows that this sparsity structure shall represent.
*/
unsigned int rows;
/**
- * Number of columns that this sparsity
- * structure shall represent.
+ * Number of columns that this sparsity structure shall represent.
*/
unsigned int cols;
unsigned int chunk_size;
/**
- * The reduced sparsity pattern. We store
- * only which chunks exist, with each
- * chunk a block in the matrix of size
- * chunk_size by chunk_size.
+ * The reduced sparsity pattern. We store only which chunks exist, with each
+ * chunk a block in the matrix of size chunk_size by chunk_size.
*/
SparsityPattern sparsity_pattern;
/**
- * Make all the chunk sparse matrix kinds
- * friends.
+ * Make all the chunk sparse matrix kinds friends.
*/
template <typename> friend class ChunkSparseMatrix;
+
+ /**
+ * Make the accessor class a friend.
+ */
+ friend class ChunkSparsityPatternIterators::Accessor;
};
#ifndef DOXYGEN
+namespace ChunkSparsityPatternIterators
+{
+ inline
+ Accessor::
+ Accessor (const ChunkSparsityPattern *sparsity_pattern,
+ const unsigned int row)
+ :
+ sparsity_pattern(sparsity_pattern),
+ reduced_accessor(row==sparsity_pattern->n_rows() ?
+ *sparsity_pattern->sparsity_pattern.end() :
+ *sparsity_pattern->sparsity_pattern.
+ begin(row/sparsity_pattern->get_chunk_size())),
+ chunk_row (row==sparsity_pattern->n_rows() ? 0 :
+ row%sparsity_pattern->get_chunk_size()),
+ chunk_col (0)
+ {}
+
+
+
+ inline
+ Accessor::
+ Accessor (const ChunkSparsityPattern *sparsity_pattern)
+ :
+ sparsity_pattern(sparsity_pattern),
+ reduced_accessor(*sparsity_pattern->sparsity_pattern.end()),
+ chunk_row (0),
+ chunk_col (0)
+ {}
+
+
+
+ inline
+ bool
+ Accessor::is_valid_entry () const
+ {
+ return reduced_accessor.is_valid_entry()
+ &&
+ sparsity_pattern->get_chunk_size()*reduced_accessor.row()+chunk_row <
+ sparsity_pattern->n_rows()
+ &&
+ sparsity_pattern->get_chunk_size()*reduced_accessor.column()+chunk_col <
+ sparsity_pattern->n_cols();
+ }
+
+
+
+ inline
+ unsigned int
+ Accessor::row() const
+ {
+ Assert (is_valid_entry() == true, ExcInvalidIterator());
+
+ return sparsity_pattern->get_chunk_size()*reduced_accessor.row()+chunk_row;
+ }
+
+
+
+ inline
+ unsigned int
+ Accessor::column() const
+ {
+ Assert (is_valid_entry() == true, ExcInvalidIterator());
+
+ return sparsity_pattern->get_chunk_size()*reduced_accessor.column() +
+ chunk_col;
+ }
+
+
+
+ inline
+ unsigned int
+ Accessor::index() const
+ {
+ Assert (is_valid_entry() == true, ExcInvalidIterator());
+
+ return sparsity_pattern->get_chunk_size()*reduced_accessor.index() +
+ chunk_col;
+ }
+
+
+
+ inline
+ std::size_t
+ Accessor::reduced_index() const
+ {
+ Assert (is_valid_entry() == true, ExcInvalidIterator());
+
+ return reduced_accessor.index_within_sparsity;
+ }
+
+
+
+
+ inline
+ bool
+ Accessor::operator == (const Accessor &other) const
+ {
+ // no need to check for equality of sparsity patterns as this is done in
+ // the reduced case already and every ChunkSparsityPattern has its own
+ // reduced sparsity pattern
+ return (reduced_accessor == other.reduced_accessor &&
+ chunk_row == other.chunk_row &&
+ chunk_col == other.chunk_col);
+ }
+
+
+
+ inline
+ bool
+ Accessor::operator < (const Accessor &other) const
+ {
+ Assert (sparsity_pattern == other.sparsity_pattern,
+ ExcInternalError());
+
+ // comparison is a bit messy because of the way ChunkSparsityPattern
+ // stores entry: chunk rows run faster than the indices of the reduced
+ // sparsity pattern, but the accessors should of course compare less based
+ // on the actual row, not the reduced one.
+ if (chunk_row == other.chunk_row)
+ return (reduced_accessor.index_within_sparsity <
+ other.reduced_accessor.index_within_sparsity ||
+ (reduced_accessor.index_within_sparsity ==
+ other.reduced_accessor.index_within_sparsity &&
+ chunk_col < other.chunk_col));
+ else
+ {
+ if (reduced_accessor.index_within_sparsity ==
+ reduced_accessor.sparsity_pattern->n_nonzero_elements())
+ return false;
+ if (other.reduced_accessor.index_within_sparsity ==
+ reduced_accessor.sparsity_pattern->n_nonzero_elements())
+ return true;
+ const unsigned int row = reduced_accessor.row(),
+ other_row = other.reduced_accessor.row();
+ return (row < other_row
+ ||
+ (row == other_row
+ &&
+ (chunk_row < other.chunk_row
+ ||
+ (chunk_row == other.chunk_row
+ &&
+ (reduced_accessor.index_within_sparsity <
+ other.reduced_accessor.index_within_sparsity
+ ||
+ (reduced_accessor.index_within_sparsity ==
+ other.reduced_accessor.index_within_sparsity
+ &&
+ chunk_col < other.chunk_col))))));
+ }
+ }
+
+
+ inline
+ void
+ Accessor::advance ()
+ {
+ const unsigned int chunk_size = sparsity_pattern->get_chunk_size();
+ Assert (chunk_row < chunk_size && chunk_col < chunk_size,
+ ExcIteratorPastEnd());
+ Assert (reduced_accessor.row() * chunk_size + chunk_row <
+ sparsity_pattern->n_rows()
+ &&
+ reduced_accessor.column() * chunk_size + chunk_col <
+ sparsity_pattern->n_cols(),
+ ExcIteratorPastEnd());
+
+ ++chunk_col;
+
+ // end of chunk
+ if (chunk_col == chunk_size
+ ||
+ reduced_accessor.column() * chunk_size + chunk_col ==
+ sparsity_pattern->n_cols())
+ {
+ // end of row
+ if (reduced_accessor.index() + 1 ==
+ reduced_accessor.sparsity_pattern->row_length(reduced_accessor.row()))
+ {
+ ++chunk_row;
+ chunk_col = 0;
+ const unsigned int old_reduced_row = reduced_accessor.row();
+
+ // end of matrix
+ if (old_reduced_row * chunk_size + chunk_row ==
+ sparsity_pattern->n_rows())
+ {
+ chunk_row = 0;
+ reduced_accessor =
+ SparsityPatternIterators::Accessor(&sparsity_pattern->sparsity_pattern,
+ sparsity_pattern->sparsity_pattern.n_nonzero_elements());
+ }
+ // end of chunk rows
+ else if (chunk_row == chunk_size)
+ {
+ reduced_accessor =
+ *sparsity_pattern->sparsity_pattern.begin(old_reduced_row+1);
+ chunk_row = 0;
+ }
+ else
+ reduced_accessor =
+ *sparsity_pattern->sparsity_pattern.begin(old_reduced_row);
+ }
+ else
+ {
+ reduced_accessor.advance();
+ chunk_col = 0;
+ }
+ }
+ }
+
+
+
+ inline
+ Iterator::Iterator (const ChunkSparsityPattern *sparsity_pattern,
+ const unsigned int row)
+ :
+ accessor(sparsity_pattern, row)
+ {}
+
+
+
+ inline
+ Iterator &
+ Iterator::operator++ ()
+ {
+ accessor.advance ();
+ return *this;
+ }
+
+
+
+ inline
+ Iterator
+ Iterator::operator++ (int)
+ {
+ const Iterator iter = *this;
+ accessor.advance ();
+ return iter;
+ }
+
+
+
+ inline
+ const Accessor &
+ Iterator::operator* () const
+ {
+ return accessor;
+ }
+
+
+
+ inline
+ const Accessor *
+ Iterator::operator-> () const
+ {
+ return &accessor;
+ }
+
+
+ inline
+ bool
+ Iterator::operator == (const Iterator &other) const
+ {
+ return (accessor == other.accessor);
+ }
+
+
+
+ inline
+ bool
+ Iterator::operator != (const Iterator &other) const
+ {
+ return ! (*this == other);
+ }
+
+
+ inline
+ bool
+ Iterator::operator < (const Iterator &other) const
+ {
+ return accessor < other.accessor;
+ }
+
+}
+
+
+
+inline
+ChunkSparsityPattern::iterator
+ChunkSparsityPattern::begin () const
+{
+ return iterator(this, 0);
+}
+
+
+inline
+ChunkSparsityPattern::iterator
+ChunkSparsityPattern::end () const
+{
+ return iterator(this, n_rows());
+}
+
+
+
+inline
+ChunkSparsityPattern::iterator
+ChunkSparsityPattern::begin (const unsigned int r) const
+{
+ Assert (r<n_rows(), ExcIndexRange(r,0,n_rows()));
+ return iterator(this, r);
+}
+
+
+
+inline
+ChunkSparsityPattern::iterator
+ChunkSparsityPattern::end (const unsigned int r) const
+{
+ Assert (r<n_rows(), ExcIndexRange(r,0,n_rows()))
+ return iterator(this, r+1);
+}
+
+
inline
unsigned int
}
+
inline
bool
ChunkSparsityPattern::optimize_diagonal () const
template <typename ForwardIterator>
+inline
void
ChunkSparsityPattern::copy_from (const unsigned int n_rows,
const unsigned int n_cols,
Assert (static_cast<unsigned int>(std::distance (begin, end)) == n_rows,
ExcIteratorRange (std::distance (begin, end), n_rows));
- // first determine row lengths for
- // each row. if the matrix is
- // quadratic, then we might have to
- // add an additional entry for the
- // diagonal, if that is not yet
- // present. as we have to call
- // compress anyway later on, don't
- // bother to check whether that
- // diagonal entry is in a certain
- // row or not
+ // first determine row lengths for each row. if the matrix is quadratic,
+ // then we might have to add an additional entry for the diagonal, if that
+ // is not yet present. as we have to call compress anyway later on, don't
+ // bother to check whether that diagonal entry is in a certain row or not
const bool is_square = (n_rows == n_cols);
std::vector<unsigned int> row_lengths;
row_lengths.reserve(n_rows);
(is_square ? 1 : 0));
reinit (n_rows, n_cols, row_lengths, chunk_size);
- // now enter all the elements into
- // the matrix
+ // now enter all the elements into the matrix
unsigned int row = 0;
typedef typename std::iterator_traits<ForwardIterator>::value_type::const_iterator inner_iterator;
for (ForwardIterator i=begin; i!=end; ++i, ++row)
for (inner_iterator j=i->begin(); j!=end_of_row; ++j)
{
const unsigned int col
- = internal::SparsityPatternTools::get_column_index_from_iterator(*j);
+ = internal::SparsityPatternTools::get_column_index_from_iterator(*j);
Assert (col < n_cols, ExcInvalidIndex(col,n_cols));
add (row, col);
}
}
- // finally compress
- // everything. this also sorts the
- // entries within each row
+ // finally compress everything. this also sorts the entries within each row
compress ();
}
Assert (s.rows == 0, ExcInvalidConstructorCall());
Assert (s.cols == 0, ExcInvalidConstructorCall());
- // perform the checks in the underlying
- // object as well
+ // perform the checks in the underlying object as well
sparsity_pattern = s.sparsity_pattern;
return *this;
{
Assert (chunk_size > 0, ExcInvalidNumber (chunk_size));
- // simply map this function to the
- // other @p{reinit} function
+ // simply map this function to the other @p{reinit} function
const std::vector<unsigned int> row_lengths (m, max_per_row);
reinit (m, n, row_lengths, chunk_size);
}
this->chunk_size = chunk_size;
- // pass down to the necessary information
- // to the underlying object. we need to
- // calculate how many chunks we need: we
- // need to round up (m/chunk_size) and
- // (n/chunk_size). rounding up in integer
- // arithmetic equals
+ // pass down to the necessary information to the underlying object. we need
+ // to calculate how many chunks we need: we need to round up (m/chunk_size)
+ // and (n/chunk_size). rounding up in integer arithmetic equals
// ((m+chunk_size-1)/chunk_size):
const unsigned int m_chunks = (m+chunk_size-1) / chunk_size,
n_chunks = (n+chunk_size-1) / chunk_size;
- // compute the maximum number of chunks in
- // each row. the passed array denotes the
- // number of entries in each row of the big
- // matrix -- in the worst case, these are
- // all in independent chunks, so we have to
- // calculate it as follows (as an example:
- // let chunk_size==2,
- // row_lengths={2,2,...}, and entries in
- // row zero at columns {0,2} and for row
- // one at {4,6} --> we'll need 4 chunks for
- // the first chunk row!) :
+ // compute the maximum number of chunks in each row. the passed array
+ // denotes the number of entries in each row of the big matrix -- in the
+ // worst case, these are all in independent chunks, so we have to calculate
+ // it as follows (as an example: let chunk_size==2, row_lengths={2,2,...},
+ // and entries in row zero at columns {0,2} and for row one at {4,6} -->
+ // we'll need 4 chunks for the first chunk row!) :
std::vector<unsigned int> chunk_row_lengths (m_chunks, 0);
for (unsigned int i=0; i<m; ++i)
chunk_row_lengths[i/chunk_size] += row_lengths[i];
+ // for the case that the reduced sparsity pattern optimizes the diagonal but
+ // the actual sparsity pattern does not, need to take one more entry in the
+ // row to fit the user-required entry
+ if (m != n && m_chunks == n_chunks)
+ for (unsigned int i=0; i<m_chunks; ++i)
+ ++chunk_row_lengths[i];
+
sparsity_pattern.reinit (m_chunks,
n_chunks,
chunk_row_lengths);
{
Assert (chunk_size > 0, ExcInvalidNumber (chunk_size));
- // count number of entries per row, then
- // initialize the underlying sparsity
+ // count number of entries per row, then initialize the underlying sparsity
// pattern
std::vector<unsigned int> entries_per_row (csp.n_rows(), 0);
for (unsigned int row = 0; row<csp.n_rows(); ++row)
+namespace internal
+{
+ namespace
+ {
+ template <typename SP>
+ void copy_sparsity (const SP &src,
+ SparsityPattern &dst)
+ {
+ dst.copy_from(src);
+ }
+
+ void copy_sparsity (const SparsityPattern &src,
+ SparsityPattern &dst)
+ {
+ dst = src;
+ }
+ }
+}
+
+
+
+template <typename Sparsity>
+void
+ChunkSparsityPattern::create_from
+(const unsigned int m,
+ const unsigned int n,
+ const Sparsity &sparsity_pattern_for_chunks,
+ const unsigned int chunk_size_in,
+ const bool)
+{
+ Assert (m > (sparsity_pattern_for_chunks.n_rows()-1) * chunk_size_in &&
+ m <= sparsity_pattern_for_chunks.n_rows() * chunk_size_in,
+ ExcMessage("Number of rows m is not compatible with chunk size "
+ "and number of rows in sparsity pattern for the chunks."));
+ Assert (n > (sparsity_pattern_for_chunks.n_cols()-1) * chunk_size_in &&
+ n <= sparsity_pattern_for_chunks.n_cols() * chunk_size_in,
+ ExcMessage("Number of columns m is not compatible with chunk size "
+ "and number of columns in sparsity pattern for the chunks."));
+
+ internal::copy_sparsity(sparsity_pattern_for_chunks, sparsity_pattern);
+ chunk_size = chunk_size_in;
+ rows = m;
+ cols = n;
+}
+
+
+
bool
ChunkSparsityPattern::empty () const
{
-unsigned int
-ChunkSparsityPattern::row_length (const unsigned int i) const
+void
+ChunkSparsityPattern::symmetrize ()
{
- Assert (i<rows, ExcIndexRange(i,0,rows));
+ // matrix must be square. note that the for some matrix sizes, the current
+ // sparsity pattern may not be square even if the underlying sparsity
+ // pattern is (e.g. a 10x11 matrix with chunk_size 4)
+ Assert (rows==cols, ExcNotQuadratic());
- return sparsity_pattern.row_length (i/chunk_size) * chunk_size;
+ sparsity_pattern.symmetrize ();
}
-void
-ChunkSparsityPattern::symmetrize ()
+unsigned int
+ChunkSparsityPattern::row_length (const unsigned int i) const
{
- // matrix must be square. note that the for
- // some matrix sizes, the current sparsity
- // pattern may not be square even if the
- // underlying sparsity pattern is (e.g. a
- // 10x11 matrix with chunk_size 4)
- Assert (rows==cols, ExcNotQuadratic());
+ Assert (i<rows, ExcIndexRange(i,0,rows));
- sparsity_pattern.symmetrize ();
+ // find out if we did padding and if this row is affected by it
+ if (n_cols() % chunk_size == 0)
+ return sparsity_pattern.row_length (i/chunk_size) * chunk_size;
+ else
+ // if columns don't align, then just iterate over all chunks and see
+ // what this leads to
+ {
+ SparsityPattern::const_iterator p = sparsity_pattern.begin(i/chunk_size),
+ end = sparsity_pattern.end(i/chunk_size);
+ unsigned int n = 0;
+ for ( ; p != end; ++p)
+ if (p->column() != sparsity_pattern.n_cols() - 1)
+ n += chunk_size;
+ else
+ n += (n_cols() % chunk_size);
+ return n;
+ }
}
chunk_size *
chunk_size);
else
- // some of the chunks reach beyond the
- // extent of this matrix. this requires a
- // somewhat more complicated
- // computations, in particular if the
+ // some of the chunks reach beyond the extent of this matrix. this
+ // requires a somewhat more complicated computations, in particular if the
// columns don't align
{
if ((n_rows() % chunk_size != 0)
&&
(n_cols() % chunk_size == 0))
{
- // columns align with chunks, but
- // not rows
+ // columns align with chunks, but not rows
unsigned int n = sparsity_pattern.n_nonzero_elements() *
chunk_size *
chunk_size;
else
{
- // if columns don't align, then
- // just iterate over all chunks and
- // see what this leads to. follow the advice in the documentation of
- // the sparsity pattern iterators to do the loop over individual rows,
+ // if columns don't align, then just iterate over all chunks and see
+ // what this leads to. follow the advice in the documentation of the
+ // sparsity pattern iterators to do the loop over individual rows,
// rather than all elements
unsigned int n = 0;
else if ((row == sparsity_pattern.n_rows() - 1)
&&
(p->column() != sparsity_pattern.n_cols() - 1))
- // last chunk row, but not
- // last chunk column. only a
- // smaller number (n_rows %
- // chunk_size) of rows
- // actually exist
+ // last chunk row, but not last chunk column. only a smaller
+ // number (n_rows % chunk_size) of rows actually exist
n += (n_rows() % chunk_size) * chunk_size;
else if ((row != sparsity_pattern.n_rows() - 1)
&&
(p->column() == sparsity_pattern.n_cols() - 1))
- // last chunk column, but
- // not row
+ // last chunk column, but not row
n += (n_cols() % chunk_size) * chunk_size;
else
// bottom right chunk
AssertThrow (out, ExcIO());
- // for each entry in the underlying
- // sparsity pattern, repeat everything
+ // for each entry in the underlying sparsity pattern, repeat everything
// chunk_size x chunk_size times
for (unsigned int i=0; i<sparsity_pattern.rows; ++i)
for (unsigned int j=sparsity_pattern.rowstart[i];
for (unsigned int e=0;
(e<chunk_size) && (i*chunk_size + e < n_rows());
++e)
- // while matrix entries are
- // usually written (i,j), with i
- // vertical and j horizontal,
- // gnuplot output is x-y, that is
- // we have to exchange the order
- // of output
+ // while matrix entries are usually written (i,j), with i vertical
+ // and j horizontal, gnuplot output is x-y, that is we have to
+ // exchange the order of output
out << sparsity_pattern.colnums[j]*chunk_size+d << " "
<< -static_cast<signed int>(i*chunk_size+e)
<< std::endl;
unsigned int
ChunkSparsityPattern::bandwidth () const
{
- // calculate the bandwidth from that of the
- // underlying sparsity pattern. note that
- // even if the bandwidth of that is zero,
- // then the bandwidth of the chunky pattern
- // is chunk_size-1, if it is 1 then the
- // chunky pattern has
- // chunk_size+(chunk_size-1), etc
+ // calculate the bandwidth from that of the underlying sparsity
+ // pattern. note that even if the bandwidth of that is zero, then the
+ // bandwidth of the chunky pattern is chunk_size-1, if it is 1 then the
+ // chunky pattern has chunk_size+(chunk_size-1), etc
//
// we'll cut it off at max(n(),m())
return std::min (sparsity_pattern.bandwidth()*chunk_size
{
AssertThrow (out, ExcIO());
- // first the simple objects,
- // bracketed in [...]
+ // first the simple objects, bracketed in [...]
out << '['
<< rows << ' '
<< cols << ' '
in >> c;
AssertThrow (c == '[', ExcIO());
- // then read the underlying sparsity
- // pattern
+ // then read the underlying sparsity pattern
sparsity_pattern.block_read (in);
in >> c;
// explicit instantiations
template
+void ChunkSparsityPattern::copy_from<SparsityPattern> (const SparsityPattern &,
+ const unsigned int,
+ const bool);
+template
void ChunkSparsityPattern::copy_from<CompressedSparsityPattern> (const CompressedSparsityPattern &,
const unsigned int,
const bool);
const unsigned int,
const bool);
template
+void ChunkSparsityPattern::create_from<SparsityPattern>
+(const unsigned int,
+ const unsigned int,
+ const SparsityPattern &,
+ const unsigned int,
+ const bool);
+template
+void ChunkSparsityPattern::create_from<CompressedSparsityPattern>
+(const unsigned int,
+ const unsigned int,
+ const CompressedSparsityPattern &,
+ const unsigned int,
+ const bool);
+template
+void ChunkSparsityPattern::create_from<CompressedSetSparsityPattern>
+(const unsigned int,
+ const unsigned int,
+ const CompressedSetSparsityPattern &,
+ const unsigned int,
+ const bool);
+template
+void ChunkSparsityPattern::create_from<CompressedSimpleSparsityPattern>
+(const unsigned int,
+ const unsigned int,
+ const CompressedSimpleSparsityPattern &,
+ const unsigned int,
+ const bool);
+template
void ChunkSparsityPattern::copy_from<float> (const FullMatrix<float> &,
const unsigned int,
const bool);
--- /dev/null
+//----------------------- chunk_sparse_matrix_11.cc ------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 2004, 2005, 2008, 2013 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//----------------------- chunk_sparse_matrix_11.cc ------------------------
+
+
+// set a few elements in a chunk sparse matrix, read them via its iterator
+// class and print them to the log file
+
+#include "../tests.h"
+#include <deal.II/lac/chunk_sparse_matrix.h>
+#include <fstream>
+#include <iomanip>
+
+
+void test (const unsigned int chunk_size)
+{
+ deallog << "Chunk size = " << chunk_size << std::endl;
+
+ ChunkSparsityPattern sp (5,5,3,chunk_size);
+ for (unsigned int i=0; i<5; ++i)
+ for (unsigned int j=0; j<5; ++j)
+ if ((i+2*j+1) % 3 == 0)
+ sp.add (i,j);
+ sp.compress ();
+
+ ChunkSparseMatrix<double> m(sp);
+
+ // first set a few entries
+ for (unsigned int i=0; i<m.m(); ++i)
+ for (unsigned int j=0; j<m.n(); ++j)
+ if ((i+2*j+1) % 3 == 0)
+ m.set (i,j, i*j*.5+.5);
+
+ // then extract the elements (note that
+ // some may be zero or even outside the
+ // matrix
+ for (unsigned int i=0; i<m.m(); ++i)
+ {
+ deallog << "row " << i << ": ";
+ for (ChunkSparseMatrix<double>::const_iterator it = m.begin(i);
+ it != m.end(i); ++it)
+ {
+ deallog << "[" << it->column() << ","
+ << std::setprecision(2) << std::fixed << std::setw (4)
+ << it->value() << "] ";
+ }
+ deallog << std::endl;
+ }
+}
+
+
+
+int main ()
+{
+ std::ofstream logfile("chunk_sparse_matrix_11/output");
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ try
+ {
+ const unsigned int chunk_sizes[] = { 1, 2, 4, 5, 7 };
+ for (unsigned int i=0;
+ i<sizeof(chunk_sizes)/sizeof(chunk_sizes[0]);
+ ++i)
+ test (chunk_sizes[i]);
+ }
+ catch (std::exception &exc)
+ {
+ deallog << std::endl << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ deallog << "Exception on processing: " << std::endl
+ << exc.what() << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+
+ return 1;
+ }
+ catch (...)
+ {
+ deallog << std::endl << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ deallog << "Unknown exception!" << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ return 1;
+ };
+}
--- /dev/null
+
+DEAL::Chunk size = 1
+DEAL::row 0: [0,0 ] [1,0.50] [4,0.50]
+DEAL::row 1: [1,0 ] [2,1.50]
+DEAL::row 2: [2,0 ] [0,0.50] [3,3.50]
+DEAL::row 3: [3,0 ] [1,2.00] [4,6.50]
+DEAL::row 4: [4,0 ] [2,4.50]
+DEAL::Chunk size = 2
+DEAL::row 0: [0,0 ] [1,0.50] [2,0 ] [3,0 ] [4,0.50]
+DEAL::row 1: [0,0 ] [1,0 ] [2,1.50] [3,0 ] [4,0 ]
+DEAL::row 2: [2,0 ] [3,3.50] [0,0.50] [1,0 ] [4,0 ]
+DEAL::row 3: [2,0 ] [3,0 ] [0,0 ] [1,2.00] [4,6.50]
+DEAL::row 4: [4,0 ] [2,4.50] [3,0 ]
+DEAL::Chunk size = 4
+DEAL::row 0: [0,0 ] [1,0.50] [2,0 ] [3,0 ] [4,0.50]
+DEAL::row 1: [0,0 ] [1,0 ] [2,1.50] [3,0 ] [4,0 ]
+DEAL::row 2: [0,0.50] [1,0 ] [2,0 ] [3,3.50] [4,0 ]
+DEAL::row 3: [0,0 ] [1,2.00] [2,0 ] [3,0 ] [4,6.50]
+DEAL::row 4: [4,0 ] [0,0 ] [1,0 ] [2,4.50] [3,0 ]
+DEAL::Chunk size = 5
+DEAL::row 0: [0,0 ] [1,0.50] [2,0 ] [3,0 ] [4,0.50]
+DEAL::row 1: [0,0 ] [1,0 ] [2,1.50] [3,0 ] [4,0 ]
+DEAL::row 2: [0,0.50] [1,0 ] [2,0 ] [3,3.50] [4,0 ]
+DEAL::row 3: [0,0 ] [1,2.00] [2,0 ] [3,0 ] [4,6.50]
+DEAL::row 4: [0,0 ] [1,0 ] [2,4.50] [3,0 ] [4,0 ]
+DEAL::Chunk size = 7
+DEAL::row 0: [0,0 ] [1,0.50] [2,0 ] [3,0 ] [4,0.50]
+DEAL::row 1: [0,0 ] [1,0 ] [2,1.50] [3,0 ] [4,0 ]
+DEAL::row 2: [0,0.50] [1,0 ] [2,0 ] [3,3.50] [4,0 ]
+DEAL::row 3: [0,0 ] [1,2.00] [2,0 ] [3,0 ] [4,6.50]
+DEAL::row 4: [0,0 ] [1,0 ] [2,4.50] [3,0 ] [4,0 ]
--- /dev/null
+//----------------------- chunk_sparse_matrix_12.cc ------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 2004, 2005, 2008, 2013 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//----------------------- chunk_sparse_matrix_12.cc ------------------------
+
+
+// set a few elements in a non-square chunk sparse matrix, read them via its
+// iterator class and print them to the log file
+
+#include "../tests.h"
+#include <deal.II/lac/chunk_sparse_matrix.h>
+#include <fstream>
+#include <iomanip>
+
+
+void test (const unsigned int chunk_size)
+{
+ deallog << "Chunk size = " << chunk_size << std::endl;
+
+ ChunkSparsityPattern sp (5,6,3,chunk_size);
+ for (unsigned int i=0; i<5; ++i)
+ for (unsigned int j=0; j<6; ++j)
+ if ((i+2*j+1) % 3 == 0)
+ sp.add (i,j);
+ sp.compress ();
+
+ ChunkSparseMatrix<double> m(sp);
+
+ // first set a few entries
+ for (unsigned int i=0; i<m.m(); ++i)
+ for (unsigned int j=0; j<m.n(); ++j)
+ if ((i+2*j+1) % 3 == 0)
+ m.set (i,j, i*j*.5+.5);
+
+ // then extract the elements (note that
+ // some may be zero or even outside the
+ // matrix
+ for (unsigned int i=0; i<m.m(); ++i)
+ {
+ deallog << "row " << i << ": ";
+ for (ChunkSparseMatrix<double>::const_iterator it = m.begin(i);
+ it != m.end(i); ++it)
+ {
+ deallog << "[" << it->column() << ","
+ << std::setprecision(2) << std::fixed << std::setw (4)
+ << it->value() << "] ";
+ }
+ deallog << std::endl;
+ }
+}
+
+
+
+int main ()
+{
+ std::ofstream logfile("chunk_sparse_matrix_12/output");
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ try
+ {
+ const unsigned int chunk_sizes[] = { 1, 2, 4, 5, 7 };
+ for (unsigned int i=0;
+ i<sizeof(chunk_sizes)/sizeof(chunk_sizes[0]);
+ ++i)
+ test (chunk_sizes[i]);
+ }
+ catch (std::exception &exc)
+ {
+ deallog << std::endl << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ deallog << "Exception on processing: " << std::endl
+ << exc.what() << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+
+ return 1;
+ }
+ catch (...)
+ {
+ deallog << std::endl << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ deallog << "Unknown exception!" << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ return 1;
+ };
+}
--- /dev/null
+
+DEAL::Chunk size = 1
+DEAL::row 0: [1,0.50] [4,0.50]
+DEAL::row 1: [2,1.50] [5,3.00]
+DEAL::row 2: [0,0.50] [3,3.50]
+DEAL::row 3: [1,2.00] [4,6.50]
+DEAL::row 4: [2,4.50] [5,10.50]
+DEAL::Chunk size = 2
+DEAL::row 0: [0,0 ] [1,0.50] [2,0 ] [3,0 ] [4,0.50] [5,0 ]
+DEAL::row 1: [0,0 ] [1,0 ] [2,1.50] [3,0 ] [4,0 ] [5,3.00]
+DEAL::row 2: [2,0 ] [3,3.50] [0,0.50] [1,0 ] [4,0 ] [5,0 ]
+DEAL::row 3: [2,0 ] [3,0 ] [0,0 ] [1,2.00] [4,6.50] [5,0 ]
+DEAL::row 4: [4,0 ] [5,10.50] [2,4.50] [3,0 ]
+DEAL::Chunk size = 4
+DEAL::row 0: [0,0 ] [1,0.50] [2,0 ] [3,0 ] [4,0.50] [5,0 ]
+DEAL::row 1: [0,0 ] [1,0 ] [2,1.50] [3,0 ] [4,0 ] [5,3.00]
+DEAL::row 2: [0,0.50] [1,0 ] [2,0 ] [3,3.50] [4,0 ] [5,0 ]
+DEAL::row 3: [0,0 ] [1,2.00] [2,0 ] [3,0 ] [4,6.50] [5,0 ]
+DEAL::row 4: [4,0 ] [5,10.50] [0,0 ] [1,0 ] [2,4.50] [3,0 ]
+DEAL::Chunk size = 5
+DEAL::row 0: [0,0 ] [1,0.50] [2,0 ] [3,0 ] [4,0.50] [5,0 ]
+DEAL::row 1: [0,0 ] [1,0 ] [2,1.50] [3,0 ] [4,0 ] [5,3.00]
+DEAL::row 2: [0,0.50] [1,0 ] [2,0 ] [3,3.50] [4,0 ] [5,0 ]
+DEAL::row 3: [0,0 ] [1,2.00] [2,0 ] [3,0 ] [4,6.50] [5,0 ]
+DEAL::row 4: [0,0 ] [1,0 ] [2,4.50] [3,0 ] [4,0 ] [5,10.50]
+DEAL::Chunk size = 7
+DEAL::row 0: [0,0 ] [1,0.50] [2,0 ] [3,0 ] [4,0.50] [5,0 ]
+DEAL::row 1: [0,0 ] [1,0 ] [2,1.50] [3,0 ] [4,0 ] [5,3.00]
+DEAL::row 2: [0,0.50] [1,0 ] [2,0 ] [3,3.50] [4,0 ] [5,0 ]
+DEAL::row 3: [0,0 ] [1,2.00] [2,0 ] [3,0 ] [4,6.50] [5,0 ]
+DEAL::row 4: [0,0 ] [1,0 ] [2,4.50] [3,0 ] [4,0 ] [5,10.50]
--- /dev/null
+//----------------------- chunk_sparse_matrix_13.cc ------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 2004, 2005, 2008, 2013 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//----------------------- chunk_sparse_matrix_13.cc ------------------------
+
+
+// set a few elements in a non-square chunk sparse matrix, read them via its
+// iterator class and print them to the log file
+
+#include "../tests.h"
+#include <deal.II/lac/chunk_sparse_matrix.h>
+#include <fstream>
+#include <iomanip>
+
+
+void test (const unsigned int chunk_size)
+{
+ deallog << "Chunk size = " << chunk_size << std::endl;
+
+ ChunkSparsityPattern sp (6,5,3,chunk_size);
+ for (unsigned int i=0; i<6; ++i)
+ for (unsigned int j=0; j<5; ++j)
+ if ((i+2*j+1) % 3 == 0)
+ sp.add (i,j);
+ sp.compress ();
+
+ ChunkSparseMatrix<double> m(sp);
+
+ // first set a few entries
+ for (unsigned int i=0; i<m.m(); ++i)
+ for (unsigned int j=0; j<m.n(); ++j)
+ if ((i+2*j+1) % 3 == 0)
+ m.set (i,j, i*j*.5+.5);
+
+ // then extract the elements (note that
+ // some may be zero or even outside the
+ // matrix
+ for (unsigned int i=0; i<m.m(); ++i)
+ {
+ deallog << "row " << i << ": ";
+ for (ChunkSparseMatrix<double>::const_iterator it = m.begin(i);
+ it != m.end(i); ++it)
+ {
+ deallog << "[" << it->column() << ","
+ << std::setprecision(2) << std::fixed << std::setw (4)
+ << it->value() << "] ";
+ }
+ deallog << std::endl;
+ }
+}
+
+
+
+int main ()
+{
+ std::ofstream logfile("chunk_sparse_matrix_13/output");
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ try
+ {
+ const unsigned int chunk_sizes[] = { 1, 2, 4, 5, 7 };
+ for (unsigned int i=0;
+ i<sizeof(chunk_sizes)/sizeof(chunk_sizes[0]);
+ ++i)
+ test (chunk_sizes[i]);
+ }
+ catch (std::exception &exc)
+ {
+ deallog << std::endl << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ deallog << "Exception on processing: " << std::endl
+ << exc.what() << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+
+ return 1;
+ }
+ catch (...)
+ {
+ deallog << std::endl << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ deallog << "Unknown exception!" << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ return 1;
+ };
+}
--- /dev/null
+
+DEAL::Chunk size = 1
+DEAL::row 0: [1,0.50] [4,0.50]
+DEAL::row 1: [2,1.50]
+DEAL::row 2: [0,0.50] [3,3.50]
+DEAL::row 3: [1,2.00] [4,6.50]
+DEAL::row 4: [2,4.50]
+DEAL::row 5: [0,0.50] [3,8.00]
+DEAL::Chunk size = 2
+DEAL::row 0: [0,0 ] [1,0.50] [2,0 ] [3,0 ] [4,0.50]
+DEAL::row 1: [0,0 ] [1,0 ] [2,1.50] [3,0 ] [4,0 ]
+DEAL::row 2: [2,0 ] [3,3.50] [0,0.50] [1,0 ] [4,0 ]
+DEAL::row 3: [2,0 ] [3,0 ] [0,0 ] [1,2.00] [4,6.50]
+DEAL::row 4: [4,0 ] [0,0 ] [1,0 ] [2,4.50] [3,0 ]
+DEAL::row 5: [4,0 ] [0,0.50] [1,0 ] [2,0 ] [3,8.00]
+DEAL::Chunk size = 4
+DEAL::row 0: [0,0 ] [1,0.50] [2,0 ] [3,0 ] [4,0.50]
+DEAL::row 1: [0,0 ] [1,0 ] [2,1.50] [3,0 ] [4,0 ]
+DEAL::row 2: [0,0.50] [1,0 ] [2,0 ] [3,3.50] [4,0 ]
+DEAL::row 3: [0,0 ] [1,2.00] [2,0 ] [3,0 ] [4,6.50]
+DEAL::row 4: [4,0 ] [0,0 ] [1,0 ] [2,4.50] [3,0 ]
+DEAL::row 5: [4,0 ] [0,0.50] [1,0 ] [2,0 ] [3,8.00]
+DEAL::Chunk size = 5
+DEAL::row 0: [0,0 ] [1,0.50] [2,0 ] [3,0 ] [4,0.50]
+DEAL::row 1: [0,0 ] [1,0 ] [2,1.50] [3,0 ] [4,0 ]
+DEAL::row 2: [0,0.50] [1,0 ] [2,0 ] [3,3.50] [4,0 ]
+DEAL::row 3: [0,0 ] [1,2.00] [2,0 ] [3,0 ] [4,6.50]
+DEAL::row 4: [0,0 ] [1,0 ] [2,4.50] [3,0 ] [4,0 ]
+DEAL::row 5: [0,0.50] [1,0 ] [2,0 ] [3,8.00] [4,0 ]
+DEAL::Chunk size = 7
+DEAL::row 0: [0,0 ] [1,0.50] [2,0 ] [3,0 ] [4,0.50]
+DEAL::row 1: [0,0 ] [1,0 ] [2,1.50] [3,0 ] [4,0 ]
+DEAL::row 2: [0,0.50] [1,0 ] [2,0 ] [3,3.50] [4,0 ]
+DEAL::row 3: [0,0 ] [1,2.00] [2,0 ] [3,0 ] [4,6.50]
+DEAL::row 4: [0,0 ] [1,0 ] [2,4.50] [3,0 ] [4,0 ]
+DEAL::row 5: [0,0.50] [1,0 ] [2,0 ] [3,8.00] [4,0 ]
--- /dev/null
+//----------------------- chunk_sparse_matrix_14.cc ------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 2004, 2005, 2008, 2013 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//----------------------- chunk_sparse_matrix_14.cc ------------------------
+
+
+// set a few elements in a chunk sparse matrix and test for iterator
+// inequality
+
+#include "../tests.h"
+#include <deal.II/lac/chunk_sparse_matrix.h>
+#include <fstream>
+#include <iomanip>
+
+
+void test (const unsigned int chunk_size)
+{
+ deallog << "Chunk size = " << chunk_size << std::endl;
+
+ ChunkSparsityPattern sp (5,5,3,chunk_size);
+ for (unsigned int i=0; i<5; ++i)
+ for (unsigned int j=0; j<5; ++j)
+ if ((i+2*j+1) % 3 == 0)
+ sp.add (i,j);
+ sp.compress ();
+
+ ChunkSparseMatrix<double> m(sp);
+
+ // first set a few entries
+ for (unsigned int i=0; i<m.m(); ++i)
+ for (unsigned int j=0; j<m.n(); ++j)
+ if ((i+2*j+1) % 3 == 0)
+ m.set (i,j, i*j*.5+.5);
+
+ // then extract the elements (note that
+ // some may be zero or even outside the
+ // matrix
+ AssertDimension(m.end()-m.begin(), m.n_nonzero_elements());
+ for (unsigned int i=0; i<m.m(); ++i)
+ {
+ deallog << "row " << i << ": ";
+ AssertDimension(m.end(i)-m.begin(i),
+ m.get_sparsity_pattern().row_length(i));
+ for (ChunkSparseMatrix<double>::const_iterator it = m.begin(i);
+ it != m.end(i); ++it)
+ {
+ deallog << " done " << (it-m.begin(i)) << ", left " << (it-m.end(i));
+ }
+ deallog << std::endl;
+ }
+}
+
+
+
+int main ()
+{
+ std::ofstream logfile("chunk_sparse_matrix_14/output");
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ try
+ {
+ const unsigned int chunk_sizes[] = { 1, 2, 4, 5, 7 };
+ for (unsigned int i=0;
+ i<sizeof(chunk_sizes)/sizeof(chunk_sizes[0]);
+ ++i)
+ test (chunk_sizes[i]);
+ }
+ catch (std::exception &exc)
+ {
+ deallog << std::endl << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ deallog << "Exception on processing: " << std::endl
+ << exc.what() << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+
+ return 1;
+ }
+ catch (...)
+ {
+ deallog << std::endl << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ deallog << "Unknown exception!" << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ return 1;
+ };
+}
--- /dev/null
+
+DEAL::Chunk size = 1
+DEAL::row 0: done 0, left -3 done 1, left -2 done 2, left -1
+DEAL::row 1: done 0, left -2 done 1, left -1
+DEAL::row 2: done 0, left -3 done 1, left -2 done 2, left -1
+DEAL::row 3: done 0, left -3 done 1, left -2 done 2, left -1
+DEAL::row 4: done 0, left -2 done 1, left -1
+DEAL::Chunk size = 2
+DEAL::row 0: done 0, left -5 done 1, left -4 done 2, left -3 done 3, left -2 done 4, left -1
+DEAL::row 1: done 0, left -5 done 1, left -4 done 2, left -3 done 3, left -2 done 4, left -1
+DEAL::row 2: done 0, left -5 done 1, left -4 done 2, left -3 done 3, left -2 done 4, left -1
+DEAL::row 3: done 0, left -5 done 1, left -4 done 2, left -3 done 3, left -2 done 4, left -1
+DEAL::row 4: done 0, left -3 done 1, left -2 done 2, left -1
+DEAL::Chunk size = 4
+DEAL::row 0: done 0, left -5 done 1, left -4 done 2, left -3 done 3, left -2 done 4, left -1
+DEAL::row 1: done 0, left -5 done 1, left -4 done 2, left -3 done 3, left -2 done 4, left -1
+DEAL::row 2: done 0, left -5 done 1, left -4 done 2, left -3 done 3, left -2 done 4, left -1
+DEAL::row 3: done 0, left -5 done 1, left -4 done 2, left -3 done 3, left -2 done 4, left -1
+DEAL::row 4: done 0, left -5 done 1, left -4 done 2, left -3 done 3, left -2 done 4, left -1
+DEAL::Chunk size = 5
+DEAL::row 0: done 0, left -5 done 1, left -4 done 2, left -3 done 3, left -2 done 4, left -1
+DEAL::row 1: done 0, left -5 done 1, left -4 done 2, left -3 done 3, left -2 done 4, left -1
+DEAL::row 2: done 0, left -5 done 1, left -4 done 2, left -3 done 3, left -2 done 4, left -1
+DEAL::row 3: done 0, left -5 done 1, left -4 done 2, left -3 done 3, left -2 done 4, left -1
+DEAL::row 4: done 0, left -5 done 1, left -4 done 2, left -3 done 3, left -2 done 4, left -1
+DEAL::Chunk size = 7
+DEAL::row 0: done 0, left -5 done 1, left -4 done 2, left -3 done 3, left -2 done 4, left -1
+DEAL::row 1: done 0, left -5 done 1, left -4 done 2, left -3 done 3, left -2 done 4, left -1
+DEAL::row 2: done 0, left -5 done 1, left -4 done 2, left -3 done 3, left -2 done 4, left -1
+DEAL::row 3: done 0, left -5 done 1, left -4 done 2, left -3 done 3, left -2 done 4, left -1
+DEAL::row 4: done 0, left -5 done 1, left -4 done 2, left -3 done 3, left -2 done 4, left -1
--- /dev/null
+//----------------------- chunk_sparse_matrix_15.cc ------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 2004, 2005, 2008, 2013 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//----------------------- chunk_sparse_matrix_15.cc ------------------------
+
+
+// set a few elements in a chunk sparse matrix, perform a matrix-vector
+// product through the iterator and compare with a matrix-vector product
+
+#include "../tests.h"
+#include <deal.II/lac/chunk_sparse_matrix.h>
+#include <deal.II/lac/vector.h>
+#include <fstream>
+#include <iomanip>
+
+
+void test (const unsigned int chunk_size)
+{
+ deallog << "Chunk size = " << chunk_size << std::endl;
+
+ for (unsigned int n_cols = 4; n_cols<7; ++n_cols)
+ {
+ deallog << "n_cols = " << n_cols << std::endl;
+ ChunkSparsityPattern sp (5,n_cols,3,chunk_size);
+ for (unsigned int i=0; i<5; ++i)
+ for (unsigned int j=0; j<n_cols; ++j)
+ if ((i+2*j+1) % 3 == 0)
+ sp.add (i,j);
+ sp.compress ();
+
+ ChunkSparseMatrix<double> m(sp);
+
+ // first set a few entries
+ for (unsigned int i=0; i<m.m(); ++i)
+ for (unsigned int j=0; j<m.n(); ++j)
+ if ((i+2*j+1) % 3 == 0)
+ m.set (i,j, i*j*.5+.5);
+
+ // next perform a matrix-vector product using the entries as given by
+ // the iterator and compare it with the exact value
+ Vector<double> src(m.n()), dst(m.m()), dst_ref(m.m());
+ for (unsigned int i=0; i<src.size(); ++i)
+ src(i) = (double)rand()/RAND_MAX;
+ for (unsigned int i=0; i<m.m(); ++i)
+ {
+ double sum = 0;
+ for (ChunkSparseMatrix<double>::const_iterator it = m.begin(i);
+ it != m.end(i); ++it)
+ sum += it->value() * src(it->column());
+ dst(i) = sum;
+ }
+ m.vmult(dst_ref, src);
+ dst -= dst_ref;
+ deallog << "Error in matrix-vector product done via iterator: "
+ << dst.linfty_norm() << std::endl;
+ }
+}
+
+
+
+int main ()
+{
+ std::ofstream logfile("chunk_sparse_matrix_15/output");
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ try
+ {
+ const unsigned int chunk_sizes[] = { 1, 2, 4, 5, 7 };
+ for (unsigned int i=0;
+ i<sizeof(chunk_sizes)/sizeof(chunk_sizes[0]);
+ ++i)
+ test (chunk_sizes[i]);
+ }
+ catch (std::exception &exc)
+ {
+ deallog << std::endl << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ deallog << "Exception on processing: " << std::endl
+ << exc.what() << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+
+ return 1;
+ }
+ catch (...)
+ {
+ deallog << std::endl << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ deallog << "Unknown exception!" << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ return 1;
+ };
+}
--- /dev/null
+
+DEAL::Chunk size = 1
+DEAL::n_cols = 4
+DEAL::Error in matrix-vector product done via iterator: 0
+DEAL::n_cols = 5
+DEAL::Error in matrix-vector product done via iterator: 0
+DEAL::n_cols = 6
+DEAL::Error in matrix-vector product done via iterator: 0
+DEAL::Chunk size = 2
+DEAL::n_cols = 4
+DEAL::Error in matrix-vector product done via iterator: 0
+DEAL::n_cols = 5
+DEAL::Error in matrix-vector product done via iterator: 0
+DEAL::n_cols = 6
+DEAL::Error in matrix-vector product done via iterator: 0
+DEAL::Chunk size = 4
+DEAL::n_cols = 4
+DEAL::Error in matrix-vector product done via iterator: 0
+DEAL::n_cols = 5
+DEAL::Error in matrix-vector product done via iterator: 0
+DEAL::n_cols = 6
+DEAL::Error in matrix-vector product done via iterator: 0
+DEAL::Chunk size = 5
+DEAL::n_cols = 4
+DEAL::Error in matrix-vector product done via iterator: 0
+DEAL::n_cols = 5
+DEAL::Error in matrix-vector product done via iterator: 0
+DEAL::n_cols = 6
+DEAL::Error in matrix-vector product done via iterator: 0
+DEAL::Chunk size = 7
+DEAL::n_cols = 4
+DEAL::Error in matrix-vector product done via iterator: 0
+DEAL::n_cols = 5
+DEAL::Error in matrix-vector product done via iterator: 0
+DEAL::n_cols = 6
+DEAL::Error in matrix-vector product done via iterator: 0
--- /dev/null
+//---------------------------- chunk_sparse_matrix_iterator_01.cc ---------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 2004, 2005, 2013 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//---------------------------- chunk_sparse_matrix_iterator_01.cc ---------------------------
+
+
+// ChunkSparseMatrix::const_iterator::operator++(int)
+
+#include "../tests.h"
+#include <deal.II/lac/chunk_sparse_matrix.h>
+#include <fstream>
+#include <iomanip>
+
+
+void test (const unsigned int chunk_size)
+{
+ deallog << "Chunk size: " << chunk_size << std::endl;
+ ChunkSparsityPattern sp (5,5,3,chunk_size);
+ for (unsigned int i=0; i<5; ++i)
+ for (unsigned int j=0; j<5; ++j)
+ if ((i+2*j+1) % 3 == 0)
+ sp.add (i,j);
+ sp.compress ();
+
+ ChunkSparseMatrix<double> m(sp);
+ ChunkSparseMatrix<double>::const_iterator i = m.begin();
+ deallog << i->value() << std::endl;
+ ++i;
+ deallog << i->value() << std::endl;
+ i++;
+ deallog << i->value() << std::endl;
+
+ deallog << "OK" << std::endl;
+}
+
+
+
+int main ()
+{
+ std::ofstream logfile("chunk_sparse_matrix_iterator_01/output");
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ try
+ {
+ const unsigned int chunk_sizes[] = { 1, 2, 4, 5, 7 };
+ for (unsigned int i=0;
+ i<sizeof(chunk_sizes)/sizeof(chunk_sizes[0]);
+ ++i)
+ test (chunk_sizes[i]);
+ }
+ catch (std::exception &exc)
+ {
+ deallog << std::endl << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ deallog << "Exception on processing: " << std::endl
+ << exc.what() << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+
+ return 1;
+ }
+ catch (...)
+ {
+ deallog << std::endl << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ deallog << "Unknown exception!" << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ return 1;
+ };
+}
--- /dev/null
+
+DEAL::Chunk size: 1
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::OK
+DEAL::Chunk size: 2
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::OK
+DEAL::Chunk size: 4
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::OK
+DEAL::Chunk size: 5
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::OK
+DEAL::Chunk size: 7
+DEAL::0
+DEAL::0
+DEAL::0
+DEAL::OK
--- /dev/null
+//---------------------------- chunk_sparse_matrix_iterator_02.cc ---------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 2004, 2005, 2013 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//---------------------------- chunk_sparse_matrix_iterator_02.cc ---------------------------
+
+
+// test setting some elements and reading them back from a const matrix
+// iterator
+
+#include "../tests.h"
+#include <deal.II/lac/chunk_sparse_matrix.h>
+#include <fstream>
+#include <iomanip>
+
+
+void test (const unsigned int chunk_size)
+{
+ deallog << "Chunk size: " << chunk_size << std::endl;
+ ChunkSparsityPattern sp (5,5,3,chunk_size);
+ for (unsigned int i=0; i<5; ++i)
+ for (unsigned int j=0; j<5; ++j)
+ if (((i+2*j+1) % 3 == 0)
+ ||
+ (i==j))
+ sp.add (i,j);
+ sp.compress ();
+
+ ChunkSparseMatrix<double> m(sp);
+ for (unsigned int i=0; i<5; ++i)
+ for (unsigned int j=0; j<5; ++j)
+ if (((i+2*j+1) % 3 == 0)
+ ||
+ (i==j))
+ m.set (i,j, i*j);
+
+ ChunkSparseMatrix<double>::const_iterator i = m.begin();
+ for (; i!=m.end(); ++i)
+ {
+ deallog << i->row() << ' ' << i->column() << ' '
+ << i->value() << std::endl;
+ if (((i->row()+2*i->column()+1) % 3 == 0)
+ ||
+ (i->row()==i->column()))
+ {
+ Assert (std::fabs(i->value() - i->row()*i->column()) < 1e-14,
+ ExcInternalError());
+ }
+ else
+ Assert (i->value() == 0, ExcInternalError());
+ }
+
+ deallog << "OK" << std::endl;
+}
+
+
+
+int main ()
+{
+ std::ofstream logfile("chunk_sparse_matrix_iterator_02/output");
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ try
+ {
+ const unsigned int chunk_sizes[] = { 1, 2, 4, 5, 7 };
+ for (unsigned int i=0;
+ i<sizeof(chunk_sizes)/sizeof(chunk_sizes[0]);
+ ++i)
+ test (chunk_sizes[i]);
+ }
+ catch (std::exception &exc)
+ {
+ deallog << std::endl << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ deallog << "Exception on processing: " << std::endl
+ << exc.what() << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+
+ return 1;
+ }
+ catch (...)
+ {
+ deallog << std::endl << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ deallog << "Unknown exception!" << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ return 1;
+ };
+}
--- /dev/null
+
+DEAL::Chunk size: 1
+DEAL::0 0 0
+DEAL::0 1 0
+DEAL::0 4 0
+DEAL::1 1 1.00000
+DEAL::1 2 2.00000
+DEAL::2 2 4.00000
+DEAL::2 0 0
+DEAL::2 3 6.00000
+DEAL::3 3 9.00000
+DEAL::3 1 3.00000
+DEAL::3 4 12.0000
+DEAL::4 4 16.0000
+DEAL::4 2 8.00000
+DEAL::OK
+DEAL::Chunk size: 2
+DEAL::0 0 0
+DEAL::0 1 0
+DEAL::0 2 0
+DEAL::0 3 0
+DEAL::0 4 0
+DEAL::1 0 0
+DEAL::1 1 1.00000
+DEAL::1 2 2.00000
+DEAL::1 3 0
+DEAL::1 4 0
+DEAL::2 2 4.00000
+DEAL::2 3 6.00000
+DEAL::2 0 0
+DEAL::2 1 0
+DEAL::2 4 0
+DEAL::3 2 0
+DEAL::3 3 9.00000
+DEAL::3 0 0
+DEAL::3 1 3.00000
+DEAL::3 4 12.0000
+DEAL::4 4 16.0000
+DEAL::4 2 8.00000
+DEAL::4 3 0
+DEAL::OK
+DEAL::Chunk size: 4
+DEAL::0 0 0
+DEAL::0 1 0
+DEAL::0 2 0
+DEAL::0 3 0
+DEAL::0 4 0
+DEAL::1 0 0
+DEAL::1 1 1.00000
+DEAL::1 2 2.00000
+DEAL::1 3 0
+DEAL::1 4 0
+DEAL::2 0 0
+DEAL::2 1 0
+DEAL::2 2 4.00000
+DEAL::2 3 6.00000
+DEAL::2 4 0
+DEAL::3 0 0
+DEAL::3 1 3.00000
+DEAL::3 2 0
+DEAL::3 3 9.00000
+DEAL::3 4 12.0000
+DEAL::4 4 16.0000
+DEAL::4 0 0
+DEAL::4 1 0
+DEAL::4 2 8.00000
+DEAL::4 3 0
+DEAL::OK
+DEAL::Chunk size: 5
+DEAL::0 0 0
+DEAL::0 1 0
+DEAL::0 2 0
+DEAL::0 3 0
+DEAL::0 4 0
+DEAL::1 0 0
+DEAL::1 1 1.00000
+DEAL::1 2 2.00000
+DEAL::1 3 0
+DEAL::1 4 0
+DEAL::2 0 0
+DEAL::2 1 0
+DEAL::2 2 4.00000
+DEAL::2 3 6.00000
+DEAL::2 4 0
+DEAL::3 0 0
+DEAL::3 1 3.00000
+DEAL::3 2 0
+DEAL::3 3 9.00000
+DEAL::3 4 12.0000
+DEAL::4 0 0
+DEAL::4 1 0
+DEAL::4 2 8.00000
+DEAL::4 3 0
+DEAL::4 4 16.0000
+DEAL::OK
+DEAL::Chunk size: 7
+DEAL::0 0 0
+DEAL::0 1 0
+DEAL::0 2 0
+DEAL::0 3 0
+DEAL::0 4 0
+DEAL::1 0 0
+DEAL::1 1 1.00000
+DEAL::1 2 2.00000
+DEAL::1 3 0
+DEAL::1 4 0
+DEAL::2 0 0
+DEAL::2 1 0
+DEAL::2 2 4.00000
+DEAL::2 3 6.00000
+DEAL::2 4 0
+DEAL::3 0 0
+DEAL::3 1 3.00000
+DEAL::3 2 0
+DEAL::3 3 9.00000
+DEAL::3 4 12.0000
+DEAL::4 0 0
+DEAL::4 1 0
+DEAL::4 2 8.00000
+DEAL::4 3 0
+DEAL::4 4 16.0000
+DEAL::OK
--- /dev/null
+//---------------------------- chunk_sparse_matrix_iterator_03.cc ---------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 2004, 2005, 2013 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//---------------------------- chunk_sparse_matrix_iterator_03.cc ---------------------------
+
+
+// test setting some elements and reading them back from a non-const chunk
+// matrix iterator
+
+#include "../tests.h"
+#include <deal.II/lac/chunk_sparse_matrix.h>
+#include <fstream>
+#include <iomanip>
+
+
+void test (const unsigned int chunk_size)
+{
+ deallog << "Chunk size: " << chunk_size << std::endl;
+ ChunkSparsityPattern sp (5,5,3,chunk_size);
+ for (unsigned int i=0; i<5; ++i)
+ for (unsigned int j=0; j<5; ++j)
+ if (((i+2*j+1) % 3 == 0)
+ ||
+ (i==j))
+ sp.add (i,j);
+ sp.compress ();
+
+ ChunkSparseMatrix<double> m(sp);
+ for (unsigned int i=0; i<5; ++i)
+ for (unsigned int j=0; j<5; ++j)
+ if (((i+2*j+1) % 3 == 0)
+ ||
+ (i==j))
+ m.set (i,j, i*j);
+
+ ChunkSparseMatrix<double>::iterator i = m.begin();
+ for (; i!=m.end(); ++i)
+ {
+ deallog << i->row() << ' ' << i->column() << ' '
+ << i->value() << std::endl;
+ if (((i->row()+2*i->column()+1) % 3 == 0)
+ ||
+ (i->row()==i->column()))
+ {
+ Assert (std::fabs(i->value() - i->row()*i->column()) < 1e-14,
+ ExcInternalError());
+ }
+ else
+ Assert (i->value() == 0, ExcInternalError());
+ }
+
+ deallog << "OK" << std::endl;
+}
+
+
+
+int main ()
+{
+ std::ofstream logfile("chunk_sparse_matrix_iterator_03/output");
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ try
+ {
+ const unsigned int chunk_sizes[] = { 1, 2, 4, 5, 7 };
+ for (unsigned int i=0;
+ i<sizeof(chunk_sizes)/sizeof(chunk_sizes[0]);
+ ++i)
+ test (chunk_sizes[i]);
+ }
+ catch (std::exception &exc)
+ {
+ deallog << std::endl << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ deallog << "Exception on processing: " << std::endl
+ << exc.what() << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+
+ return 1;
+ }
+ catch (...)
+ {
+ deallog << std::endl << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ deallog << "Unknown exception!" << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ return 1;
+ };
+}
--- /dev/null
+
+DEAL::Chunk size: 1
+DEAL::0 0 0.00000
+DEAL::0 1 0.00000
+DEAL::0 4 0.00000
+DEAL::1 1 1.00000
+DEAL::1 2 2.00000
+DEAL::2 2 4.00000
+DEAL::2 0 0.00000
+DEAL::2 3 6.00000
+DEAL::3 3 9.00000
+DEAL::3 1 3.00000
+DEAL::3 4 12.0000
+DEAL::4 4 16.0000
+DEAL::4 2 8.00000
+DEAL::OK
+DEAL::Chunk size: 2
+DEAL::0 0 0.00000
+DEAL::0 1 0.00000
+DEAL::0 2 0.00000
+DEAL::0 3 0.00000
+DEAL::0 4 0.00000
+DEAL::1 0 0.00000
+DEAL::1 1 1.00000
+DEAL::1 2 2.00000
+DEAL::1 3 0.00000
+DEAL::1 4 0.00000
+DEAL::2 2 4.00000
+DEAL::2 3 6.00000
+DEAL::2 0 0.00000
+DEAL::2 1 0.00000
+DEAL::2 4 0.00000
+DEAL::3 2 0.00000
+DEAL::3 3 9.00000
+DEAL::3 0 0.00000
+DEAL::3 1 3.00000
+DEAL::3 4 12.0000
+DEAL::4 4 16.0000
+DEAL::4 2 8.00000
+DEAL::4 3 0.00000
+DEAL::OK
+DEAL::Chunk size: 4
+DEAL::0 0 0.00000
+DEAL::0 1 0.00000
+DEAL::0 2 0.00000
+DEAL::0 3 0.00000
+DEAL::0 4 0.00000
+DEAL::1 0 0.00000
+DEAL::1 1 1.00000
+DEAL::1 2 2.00000
+DEAL::1 3 0.00000
+DEAL::1 4 0.00000
+DEAL::2 0 0.00000
+DEAL::2 1 0.00000
+DEAL::2 2 4.00000
+DEAL::2 3 6.00000
+DEAL::2 4 0.00000
+DEAL::3 0 0.00000
+DEAL::3 1 3.00000
+DEAL::3 2 0.00000
+DEAL::3 3 9.00000
+DEAL::3 4 12.0000
+DEAL::4 4 16.0000
+DEAL::4 0 0.00000
+DEAL::4 1 0.00000
+DEAL::4 2 8.00000
+DEAL::4 3 0.00000
+DEAL::OK
+DEAL::Chunk size: 5
+DEAL::0 0 0.00000
+DEAL::0 1 0.00000
+DEAL::0 2 0.00000
+DEAL::0 3 0.00000
+DEAL::0 4 0.00000
+DEAL::1 0 0.00000
+DEAL::1 1 1.00000
+DEAL::1 2 2.00000
+DEAL::1 3 0.00000
+DEAL::1 4 0.00000
+DEAL::2 0 0.00000
+DEAL::2 1 0.00000
+DEAL::2 2 4.00000
+DEAL::2 3 6.00000
+DEAL::2 4 0.00000
+DEAL::3 0 0.00000
+DEAL::3 1 3.00000
+DEAL::3 2 0.00000
+DEAL::3 3 9.00000
+DEAL::3 4 12.0000
+DEAL::4 0 0.00000
+DEAL::4 1 0.00000
+DEAL::4 2 8.00000
+DEAL::4 3 0.00000
+DEAL::4 4 16.0000
+DEAL::OK
+DEAL::Chunk size: 7
+DEAL::0 0 0.00000
+DEAL::0 1 0.00000
+DEAL::0 2 0.00000
+DEAL::0 3 0.00000
+DEAL::0 4 0.00000
+DEAL::1 0 0.00000
+DEAL::1 1 1.00000
+DEAL::1 2 2.00000
+DEAL::1 3 0.00000
+DEAL::1 4 0.00000
+DEAL::2 0 0.00000
+DEAL::2 1 0.00000
+DEAL::2 2 4.00000
+DEAL::2 3 6.00000
+DEAL::2 4 0.00000
+DEAL::3 0 0.00000
+DEAL::3 1 3.00000
+DEAL::3 2 0.00000
+DEAL::3 3 9.00000
+DEAL::3 4 12.0000
+DEAL::4 0 0.00000
+DEAL::4 1 0.00000
+DEAL::4 2 8.00000
+DEAL::4 3 0.00000
+DEAL::4 4 16.0000
+DEAL::OK
--- /dev/null
+//---------------------------- chunk_sparse_matrix_iterator_04.cc ---------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 2004, 2005, 2013 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//---------------------------- chunk_sparse_matrix_iterator_04.cc ---------------------------
+
+// test setting some elements using a non-const chunk matrix iterator and
+// operator=, and reading them back through the matrix itself
+
+#include "../tests.h"
+#include <deal.II/lac/chunk_sparse_matrix.h>
+#include <fstream>
+#include <iomanip>
+
+
+void test (const unsigned int chunk_size)
+{
+ deallog << "Chunk size: " << chunk_size << std::endl;
+ ChunkSparsityPattern sp (5,5,3,chunk_size);
+ for (unsigned int i=0; i<5; ++i)
+ for (unsigned int j=0; j<5; ++j)
+ if (((i+2*j+1) % 3 == 0)
+ ||
+ (i==j))
+ sp.add (i,j);
+ sp.compress ();
+ ChunkSparseMatrix<double> m(sp);
+
+ ChunkSparseMatrix<double>::iterator i = m.begin();
+ for (; i!=m.end(); ++i)
+ i->value() = i->row()*i->column();
+
+ for (unsigned int i=0; i<5; ++i)
+ for (unsigned int j=0; j<5; ++j)
+ if (((i+2*j+1) % 3 == 0)
+ ||
+ (i==j))
+ {
+ deallog << i << ' ' << j << ' ' << m.el(i,j)
+ << std::endl;
+ Assert (std::fabs(m.el(i,j)-i*j) < 1e-14,
+ ExcInternalError());
+ }
+
+ deallog << "OK" << std::endl;
+}
+
+
+
+int main ()
+{
+ std::ofstream logfile("chunk_sparse_matrix_iterator_04/output");
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ try
+ {
+ test (1);
+ test (2);
+ test (4);
+ test (5);
+ test (7);
+ }
+ catch (std::exception &exc)
+ {
+ deallog << std::endl << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ deallog << "Exception on processing: " << std::endl
+ << exc.what() << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+
+ return 1;
+ }
+ catch (...)
+ {
+ deallog << std::endl << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ deallog << "Unknown exception!" << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ return 1;
+ };
+}
--- /dev/null
+
+DEAL::Chunk size: 1
+DEAL::0 0 0
+DEAL::0 1 0
+DEAL::0 4 0
+DEAL::1 1 1.00000
+DEAL::1 2 2.00000
+DEAL::2 0 0
+DEAL::2 2 4.00000
+DEAL::2 3 6.00000
+DEAL::3 1 3.00000
+DEAL::3 3 9.00000
+DEAL::3 4 12.0000
+DEAL::4 2 8.00000
+DEAL::4 4 16.0000
+DEAL::OK
+DEAL::Chunk size: 2
+DEAL::0 0 0
+DEAL::0 1 0
+DEAL::0 4 0
+DEAL::1 1 1.00000
+DEAL::1 2 2.00000
+DEAL::2 0 0
+DEAL::2 2 4.00000
+DEAL::2 3 6.00000
+DEAL::3 1 3.00000
+DEAL::3 3 9.00000
+DEAL::3 4 12.0000
+DEAL::4 2 8.00000
+DEAL::4 4 16.0000
+DEAL::OK
+DEAL::Chunk size: 4
+DEAL::0 0 0
+DEAL::0 1 0
+DEAL::0 4 0
+DEAL::1 1 1.00000
+DEAL::1 2 2.00000
+DEAL::2 0 0
+DEAL::2 2 4.00000
+DEAL::2 3 6.00000
+DEAL::3 1 3.00000
+DEAL::3 3 9.00000
+DEAL::3 4 12.0000
+DEAL::4 2 8.00000
+DEAL::4 4 16.0000
+DEAL::OK
+DEAL::Chunk size: 5
+DEAL::0 0 0
+DEAL::0 1 0
+DEAL::0 4 0
+DEAL::1 1 1.00000
+DEAL::1 2 2.00000
+DEAL::2 0 0
+DEAL::2 2 4.00000
+DEAL::2 3 6.00000
+DEAL::3 1 3.00000
+DEAL::3 3 9.00000
+DEAL::3 4 12.0000
+DEAL::4 2 8.00000
+DEAL::4 4 16.0000
+DEAL::OK
+DEAL::Chunk size: 7
+DEAL::0 0 0
+DEAL::0 1 0
+DEAL::0 4 0
+DEAL::1 1 1.00000
+DEAL::1 2 2.00000
+DEAL::2 0 0
+DEAL::2 2 4.00000
+DEAL::2 3 6.00000
+DEAL::3 1 3.00000
+DEAL::3 3 9.00000
+DEAL::3 4 12.0000
+DEAL::4 2 8.00000
+DEAL::4 4 16.0000
+DEAL::OK
--- /dev/null
+//---------------------------- chunk_sparse_matrix_iterator_05.cc ---------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 2004, 2005, 2013 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//---------------------------- chunk_sparse_matrix_iterator_05.cc ---------------------------
+
+// test setting some elements using a non-const chunk matrix iterator and
+// operator+=, and reading them back through the matrix itself
+
+#include "../tests.h"
+#include <deal.II/lac/chunk_sparse_matrix.h>
+#include <fstream>
+#include <iomanip>
+
+
+void test (const unsigned int chunk_size)
+{
+ deallog << "Chunk size: " << chunk_size << std::endl;
+ ChunkSparsityPattern sp (5,5,3,chunk_size);
+ for (unsigned int i=0; i<5; ++i)
+ for (unsigned int j=0; j<5; ++j)
+ if (((i+2*j+1) % 3 == 0)
+ ||
+ (i==j))
+ sp.add (i,j);
+ sp.compress ();
+
+ ChunkSparseMatrix<double> m(sp);
+ for (unsigned int i=0; i<5; ++i)
+ for (unsigned int j=0; j<5; ++j)
+ if (((i+2*j+1) % 3 == 0)
+ ||
+ (i==j))
+ m.set(i,j,1.);
+
+ ChunkSparseMatrix<double>::iterator i = m.begin();
+ for (; i!=m.end(); ++i)
+ i->value() += i->row()*i->column();
+
+ for (unsigned int i=0; i<5; ++i)
+ for (unsigned int j=0; j<5; ++j)
+ if (((i+2*j+1) % 3 == 0)
+ ||
+ (i==j))
+ {
+ deallog << i << ' ' << j << ' ' << m.el(i,j)
+ << std::endl;
+ Assert (std::fabs(m.el(i,j)-(1.+i*j)) < 1e-14,
+ ExcInternalError());
+ }
+
+ deallog << "OK" << std::endl;
+}
+
+
+
+int main ()
+{
+ std::ofstream logfile("chunk_sparse_matrix_iterator_05/output");
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ try
+ {
+ test (1);
+ test (2);
+ test (4);
+ test (5);
+ test (7);
+ }
+ catch (std::exception &exc)
+ {
+ deallog << std::endl << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ deallog << "Exception on processing: " << std::endl
+ << exc.what() << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+
+ return 1;
+ }
+ catch (...)
+ {
+ deallog << std::endl << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ deallog << "Unknown exception!" << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ return 1;
+ };
+}
--- /dev/null
+
+DEAL::Chunk size: 1
+DEAL::0 0 1.00000
+DEAL::0 1 1.00000
+DEAL::0 4 1.00000
+DEAL::1 1 2.00000
+DEAL::1 2 3.00000
+DEAL::2 0 1.00000
+DEAL::2 2 5.00000
+DEAL::2 3 7.00000
+DEAL::3 1 4.00000
+DEAL::3 3 10.0000
+DEAL::3 4 13.0000
+DEAL::4 2 9.00000
+DEAL::4 4 17.0000
+DEAL::OK
+DEAL::Chunk size: 2
+DEAL::0 0 1.00000
+DEAL::0 1 1.00000
+DEAL::0 4 1.00000
+DEAL::1 1 2.00000
+DEAL::1 2 3.00000
+DEAL::2 0 1.00000
+DEAL::2 2 5.00000
+DEAL::2 3 7.00000
+DEAL::3 1 4.00000
+DEAL::3 3 10.0000
+DEAL::3 4 13.0000
+DEAL::4 2 9.00000
+DEAL::4 4 17.0000
+DEAL::OK
+DEAL::Chunk size: 4
+DEAL::0 0 1.00000
+DEAL::0 1 1.00000
+DEAL::0 4 1.00000
+DEAL::1 1 2.00000
+DEAL::1 2 3.00000
+DEAL::2 0 1.00000
+DEAL::2 2 5.00000
+DEAL::2 3 7.00000
+DEAL::3 1 4.00000
+DEAL::3 3 10.0000
+DEAL::3 4 13.0000
+DEAL::4 2 9.00000
+DEAL::4 4 17.0000
+DEAL::OK
+DEAL::Chunk size: 5
+DEAL::0 0 1.00000
+DEAL::0 1 1.00000
+DEAL::0 4 1.00000
+DEAL::1 1 2.00000
+DEAL::1 2 3.00000
+DEAL::2 0 1.00000
+DEAL::2 2 5.00000
+DEAL::2 3 7.00000
+DEAL::3 1 4.00000
+DEAL::3 3 10.0000
+DEAL::3 4 13.0000
+DEAL::4 2 9.00000
+DEAL::4 4 17.0000
+DEAL::OK
+DEAL::Chunk size: 7
+DEAL::0 0 1.00000
+DEAL::0 1 1.00000
+DEAL::0 4 1.00000
+DEAL::1 1 2.00000
+DEAL::1 2 3.00000
+DEAL::2 0 1.00000
+DEAL::2 2 5.00000
+DEAL::2 3 7.00000
+DEAL::3 1 4.00000
+DEAL::3 3 10.0000
+DEAL::3 4 13.0000
+DEAL::4 2 9.00000
+DEAL::4 4 17.0000
+DEAL::OK
--- /dev/null
+//---------------------------- chunk_sparse_matrix_iterator_06.cc ---------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 2004, 2005, 2013 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//---------------------------- chunk_sparse_matrix_iterator_06.cc ---------------------------
+
+// test setting some elements using a non-const chunk matrix iterator and
+// operator-=, and reading them back through the matrix itself
+
+#include "../tests.h"
+#include <deal.II/lac/chunk_sparse_matrix.h>
+#include <fstream>
+#include <iomanip>
+
+
+void test (const unsigned int chunk_size)
+{
+ deallog << "Chunk size: " << chunk_size << std::endl;
+ ChunkSparsityPattern sp (5,5,3,chunk_size);
+ for (unsigned int i=0; i<5; ++i)
+ for (unsigned int j=0; j<5; ++j)
+ if (((i+2*j+1) % 3 == 0)
+ ||
+ (i==j))
+ sp.add (i,j);
+ sp.compress ();
+
+ ChunkSparseMatrix<double> m(sp);
+ for (unsigned int i=0; i<5; ++i)
+ for (unsigned int j=0; j<5; ++j)
+ if (((i+2*j+1) % 3 == 0)
+ ||
+ (i==j))
+ m.set(i,j,1.);
+
+ ChunkSparseMatrix<double>::iterator i = m.begin();
+ for (; i!=m.end(); ++i)
+ i->value() -= i->row()*i->column();
+
+ for (unsigned int i=0; i<5; ++i)
+ for (unsigned int j=0; j<5; ++j)
+ if (((i+2*j+1) % 3 == 0)
+ ||
+ (i==j))
+ {
+ deallog << i << ' ' << j << ' ' << m.el(i,j)
+ << std::endl;
+ Assert (std::fabs(m.el(i,j)-(1.-i*j)) < 1e-14,
+ ExcInternalError());
+ }
+ else if (sp.exists(i,j))
+ {
+ deallog << i << ' ' << j << ' ' << m.el(i,j)
+ << std::endl;
+ Assert (std::fabs(m.el(i,j)+i*j) < 1e-14,
+ ExcInternalError());
+ }
+
+ deallog << "OK" << std::endl;
+}
+
+
+
+int main ()
+{
+ std::ofstream logfile("chunk_sparse_matrix_iterator_06/output");
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ try
+ {
+ test (1);
+ test (2);
+ test (4);
+ test (5);
+ test (7);
+ }
+ catch (std::exception &exc)
+ {
+ deallog << std::endl << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ deallog << "Exception on processing: " << std::endl
+ << exc.what() << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+
+ return 1;
+ }
+ catch (...)
+ {
+ deallog << std::endl << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ deallog << "Unknown exception!" << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ return 1;
+ };
+}
--- /dev/null
+
+DEAL::Chunk size: 1
+DEAL::0 0 1.00000
+DEAL::0 1 1.00000
+DEAL::0 4 1.00000
+DEAL::1 1 0
+DEAL::1 2 -1.00000
+DEAL::2 0 1.00000
+DEAL::2 2 -3.00000
+DEAL::2 3 -5.00000
+DEAL::3 1 -2.00000
+DEAL::3 3 -8.00000
+DEAL::3 4 -11.0000
+DEAL::4 2 -7.00000
+DEAL::4 4 -15.0000
+DEAL::OK
+DEAL::Chunk size: 2
+DEAL::0 0 1.00000
+DEAL::0 1 1.00000
+DEAL::0 2 0
+DEAL::0 3 0
+DEAL::0 4 1.00000
+DEAL::1 0 0
+DEAL::1 1 0
+DEAL::1 2 -1.00000
+DEAL::1 3 -3.00000
+DEAL::1 4 -4.00000
+DEAL::2 0 1.00000
+DEAL::2 1 -2.00000
+DEAL::2 2 -3.00000
+DEAL::2 3 -5.00000
+DEAL::2 4 -8.00000
+DEAL::3 0 0
+DEAL::3 1 -2.00000
+DEAL::3 2 -6.00000
+DEAL::3 3 -8.00000
+DEAL::3 4 -11.0000
+DEAL::4 2 -7.00000
+DEAL::4 3 -12.0000
+DEAL::4 4 -15.0000
+DEAL::OK
+DEAL::Chunk size: 4
+DEAL::0 0 1.00000
+DEAL::0 1 1.00000
+DEAL::0 2 0
+DEAL::0 3 0
+DEAL::0 4 1.00000
+DEAL::1 0 0
+DEAL::1 1 0
+DEAL::1 2 -1.00000
+DEAL::1 3 -3.00000
+DEAL::1 4 -4.00000
+DEAL::2 0 1.00000
+DEAL::2 1 -2.00000
+DEAL::2 2 -3.00000
+DEAL::2 3 -5.00000
+DEAL::2 4 -8.00000
+DEAL::3 0 0
+DEAL::3 1 -2.00000
+DEAL::3 2 -6.00000
+DEAL::3 3 -8.00000
+DEAL::3 4 -11.0000
+DEAL::4 0 0
+DEAL::4 1 -4.00000
+DEAL::4 2 -7.00000
+DEAL::4 3 -12.0000
+DEAL::4 4 -15.0000
+DEAL::OK
+DEAL::Chunk size: 5
+DEAL::0 0 1.00000
+DEAL::0 1 1.00000
+DEAL::0 2 0
+DEAL::0 3 0
+DEAL::0 4 1.00000
+DEAL::1 0 0
+DEAL::1 1 0
+DEAL::1 2 -1.00000
+DEAL::1 3 -3.00000
+DEAL::1 4 -4.00000
+DEAL::2 0 1.00000
+DEAL::2 1 -2.00000
+DEAL::2 2 -3.00000
+DEAL::2 3 -5.00000
+DEAL::2 4 -8.00000
+DEAL::3 0 0
+DEAL::3 1 -2.00000
+DEAL::3 2 -6.00000
+DEAL::3 3 -8.00000
+DEAL::3 4 -11.0000
+DEAL::4 0 0
+DEAL::4 1 -4.00000
+DEAL::4 2 -7.00000
+DEAL::4 3 -12.0000
+DEAL::4 4 -15.0000
+DEAL::OK
+DEAL::Chunk size: 7
+DEAL::0 0 1.00000
+DEAL::0 1 1.00000
+DEAL::0 2 0
+DEAL::0 3 0
+DEAL::0 4 1.00000
+DEAL::1 0 0
+DEAL::1 1 0
+DEAL::1 2 -1.00000
+DEAL::1 3 -3.00000
+DEAL::1 4 -4.00000
+DEAL::2 0 1.00000
+DEAL::2 1 -2.00000
+DEAL::2 2 -3.00000
+DEAL::2 3 -5.00000
+DEAL::2 4 -8.00000
+DEAL::3 0 0
+DEAL::3 1 -2.00000
+DEAL::3 2 -6.00000
+DEAL::3 3 -8.00000
+DEAL::3 4 -11.0000
+DEAL::4 0 0
+DEAL::4 1 -4.00000
+DEAL::4 2 -7.00000
+DEAL::4 3 -12.0000
+DEAL::4 4 -15.0000
+DEAL::OK
--- /dev/null
+//---------------------------- chunk_sparse_matrix_iterator_07.cc ---------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 2004, 2005, 2013 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//---------------------------- chunk_sparse_matrix_iterator_07.cc ---------------------------
+
+// test setting some elements using a non-const chunk matrix iterator and
+// operator*=, and reading them back through the matrix itself
+
+#include "../tests.h"
+#include <deal.II/lac/chunk_sparse_matrix.h>
+#include <fstream>
+#include <iomanip>
+
+
+void test (const unsigned int chunk_size)
+{
+ deallog << "Chunk size: " << chunk_size << std::endl;
+ ChunkSparsityPattern sp (5,5,3,chunk_size);
+ for (unsigned int i=0; i<5; ++i)
+ for (unsigned int j=0; j<5; ++j)
+ if (((i+2*j+1) % 3 == 0)
+ ||
+ (i==j))
+ sp.add (i,j);
+ sp.compress ();
+
+ ChunkSparseMatrix<double> m(sp);
+ for (unsigned int i=0; i<5; ++i)
+ for (unsigned int j=0; j<5; ++j)
+ if (((i+2*j+1) % 3 == 0)
+ ||
+ (i==j))
+ m.set(i,j,i*j);
+
+ ChunkSparseMatrix<double>::iterator i = m.begin();
+ for (; i!=m.end(); ++i)
+ i->value() *= 2;
+
+ for (unsigned int i=0; i<5; ++i)
+ for (unsigned int j=0; j<5; ++j)
+ if (((i+2*j+1) % 3 == 0)
+ ||
+ (i==j))
+ {
+ deallog << i << ' ' << j << ' ' << m.el(i,j)
+ << std::endl;
+ Assert (std::fabs(m.el(i,j)-2*i*j) < 1e-14,
+ ExcInternalError());
+ }
+ else if (sp.exists(i,j))
+ {
+ deallog << i << ' ' << j << ' ' << m.el(i,j)
+ << std::endl;
+ Assert (std::fabs(m.el(i,j)) < 1e-14,
+ ExcInternalError());
+ }
+
+ deallog << "OK" << std::endl;
+}
+
+
+
+int main ()
+{
+ std::ofstream logfile("chunk_sparse_matrix_iterator_07/output");
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ try
+ {
+ test (1);
+ test (2);
+ test (4);
+ test (5);
+ test (7);
+ }
+ catch (std::exception &exc)
+ {
+ deallog << std::endl << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ deallog << "Exception on processing: " << std::endl
+ << exc.what() << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+
+ return 1;
+ }
+ catch (...)
+ {
+ deallog << std::endl << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ deallog << "Unknown exception!" << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ return 1;
+ };
+}
--- /dev/null
+
+DEAL::Chunk size: 1
+DEAL::0 0 0
+DEAL::0 1 0
+DEAL::0 4 0
+DEAL::1 1 2.00000
+DEAL::1 2 4.00000
+DEAL::2 0 0
+DEAL::2 2 8.00000
+DEAL::2 3 12.0000
+DEAL::3 1 6.00000
+DEAL::3 3 18.0000
+DEAL::3 4 24.0000
+DEAL::4 2 16.0000
+DEAL::4 4 32.0000
+DEAL::OK
+DEAL::Chunk size: 2
+DEAL::0 0 0
+DEAL::0 1 0
+DEAL::0 2 0
+DEAL::0 3 0
+DEAL::0 4 0
+DEAL::1 0 0
+DEAL::1 1 2.00000
+DEAL::1 2 4.00000
+DEAL::1 3 0
+DEAL::1 4 0
+DEAL::2 0 0
+DEAL::2 1 0
+DEAL::2 2 8.00000
+DEAL::2 3 12.0000
+DEAL::2 4 0
+DEAL::3 0 0
+DEAL::3 1 6.00000
+DEAL::3 2 0
+DEAL::3 3 18.0000
+DEAL::3 4 24.0000
+DEAL::4 2 16.0000
+DEAL::4 3 0
+DEAL::4 4 32.0000
+DEAL::OK
+DEAL::Chunk size: 4
+DEAL::0 0 0
+DEAL::0 1 0
+DEAL::0 2 0
+DEAL::0 3 0
+DEAL::0 4 0
+DEAL::1 0 0
+DEAL::1 1 2.00000
+DEAL::1 2 4.00000
+DEAL::1 3 0
+DEAL::1 4 0
+DEAL::2 0 0
+DEAL::2 1 0
+DEAL::2 2 8.00000
+DEAL::2 3 12.0000
+DEAL::2 4 0
+DEAL::3 0 0
+DEAL::3 1 6.00000
+DEAL::3 2 0
+DEAL::3 3 18.0000
+DEAL::3 4 24.0000
+DEAL::4 0 0
+DEAL::4 1 0
+DEAL::4 2 16.0000
+DEAL::4 3 0
+DEAL::4 4 32.0000
+DEAL::OK
+DEAL::Chunk size: 5
+DEAL::0 0 0
+DEAL::0 1 0
+DEAL::0 2 0
+DEAL::0 3 0
+DEAL::0 4 0
+DEAL::1 0 0
+DEAL::1 1 2.00000
+DEAL::1 2 4.00000
+DEAL::1 3 0
+DEAL::1 4 0
+DEAL::2 0 0
+DEAL::2 1 0
+DEAL::2 2 8.00000
+DEAL::2 3 12.0000
+DEAL::2 4 0
+DEAL::3 0 0
+DEAL::3 1 6.00000
+DEAL::3 2 0
+DEAL::3 3 18.0000
+DEAL::3 4 24.0000
+DEAL::4 0 0
+DEAL::4 1 0
+DEAL::4 2 16.0000
+DEAL::4 3 0
+DEAL::4 4 32.0000
+DEAL::OK
+DEAL::Chunk size: 7
+DEAL::0 0 0
+DEAL::0 1 0
+DEAL::0 2 0
+DEAL::0 3 0
+DEAL::0 4 0
+DEAL::1 0 0
+DEAL::1 1 2.00000
+DEAL::1 2 4.00000
+DEAL::1 3 0
+DEAL::1 4 0
+DEAL::2 0 0
+DEAL::2 1 0
+DEAL::2 2 8.00000
+DEAL::2 3 12.0000
+DEAL::2 4 0
+DEAL::3 0 0
+DEAL::3 1 6.00000
+DEAL::3 2 0
+DEAL::3 3 18.0000
+DEAL::3 4 24.0000
+DEAL::4 0 0
+DEAL::4 1 0
+DEAL::4 2 16.0000
+DEAL::4 3 0
+DEAL::4 4 32.0000
+DEAL::OK
--- /dev/null
+//---------------------------- chunk_sparse_matrix_iterator_08.cc ---------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 2004, 2005, 2013 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//---------------------------- chunk_sparse_matrix_iterator_08.cc ---------------------------
+
+// test setting some elements using a non-const chunk matrix iterator and
+// operator/=, and reading them back through the matrix itself
+
+#include "../tests.h"
+#include <deal.II/lac/chunk_sparse_matrix.h>
+#include <fstream>
+#include <iomanip>
+
+
+void test (const unsigned int chunk_size)
+{
+ deallog << "Chunk size: " << chunk_size << std::endl;
+ ChunkSparsityPattern sp (5,5,3,chunk_size);
+ for (unsigned int i=0; i<5; ++i)
+ for (unsigned int j=0; j<5; ++j)
+ if (((i+2*j+1) % 3 == 0)
+ ||
+ (i==j))
+ sp.add (i,j);
+ sp.compress ();
+
+ ChunkSparseMatrix<double> m(sp);
+ for (unsigned int i=0; i<5; ++i)
+ for (unsigned int j=0; j<5; ++j)
+ if (((i+2*j+1) % 3 == 0)
+ ||
+ (i==j))
+ m.set(i,j,i*j);
+
+ ChunkSparseMatrix<double>::iterator i = m.begin();
+ for (; i!=m.end(); ++i)
+ i->value() /= 2;
+
+ for (unsigned int i=0; i<5; ++i)
+ for (unsigned int j=0; j<5; ++j)
+ if (((i+2*j+1) % 3 == 0)
+ ||
+ (i==j))
+ {
+ deallog << i << ' ' << j << ' ' << m.el(i,j)
+ << std::endl;
+ Assert (std::fabs(m.el(i,j)-(i*j/2.)) < 1e-14,
+ ExcInternalError());
+ }
+ else if (sp.exists(i,j))
+ {
+ deallog << i << ' ' << j << ' ' << m.el(i,j)
+ << std::endl;
+ Assert (std::fabs(m.el(i,j)) < 1e-14,
+ ExcInternalError());
+ }
+
+ deallog << "OK" << std::endl;
+}
+
+
+
+int main ()
+{
+ std::ofstream logfile("chunk_sparse_matrix_iterator_08/output");
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ try
+ {
+ test (1);
+ test (2);
+ test (4);
+ test (5);
+ test (7);
+ }
+ catch (std::exception &exc)
+ {
+ deallog << std::endl << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ deallog << "Exception on processing: " << std::endl
+ << exc.what() << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+
+ return 1;
+ }
+ catch (...)
+ {
+ deallog << std::endl << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ deallog << "Unknown exception!" << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ return 1;
+ };
+}
--- /dev/null
+
+DEAL::Chunk size: 1
+DEAL::0 0 0
+DEAL::0 1 0
+DEAL::0 4 0
+DEAL::1 1 0.500000
+DEAL::1 2 1.00000
+DEAL::2 0 0
+DEAL::2 2 2.00000
+DEAL::2 3 3.00000
+DEAL::3 1 1.50000
+DEAL::3 3 4.50000
+DEAL::3 4 6.00000
+DEAL::4 2 4.00000
+DEAL::4 4 8.00000
+DEAL::OK
+DEAL::Chunk size: 2
+DEAL::0 0 0
+DEAL::0 1 0
+DEAL::0 2 0
+DEAL::0 3 0
+DEAL::0 4 0
+DEAL::1 0 0
+DEAL::1 1 0.500000
+DEAL::1 2 1.00000
+DEAL::1 3 0
+DEAL::1 4 0
+DEAL::2 0 0
+DEAL::2 1 0
+DEAL::2 2 2.00000
+DEAL::2 3 3.00000
+DEAL::2 4 0
+DEAL::3 0 0
+DEAL::3 1 1.50000
+DEAL::3 2 0
+DEAL::3 3 4.50000
+DEAL::3 4 6.00000
+DEAL::4 2 4.00000
+DEAL::4 3 0
+DEAL::4 4 8.00000
+DEAL::OK
+DEAL::Chunk size: 4
+DEAL::0 0 0
+DEAL::0 1 0
+DEAL::0 2 0
+DEAL::0 3 0
+DEAL::0 4 0
+DEAL::1 0 0
+DEAL::1 1 0.500000
+DEAL::1 2 1.00000
+DEAL::1 3 0
+DEAL::1 4 0
+DEAL::2 0 0
+DEAL::2 1 0
+DEAL::2 2 2.00000
+DEAL::2 3 3.00000
+DEAL::2 4 0
+DEAL::3 0 0
+DEAL::3 1 1.50000
+DEAL::3 2 0
+DEAL::3 3 4.50000
+DEAL::3 4 6.00000
+DEAL::4 0 0
+DEAL::4 1 0
+DEAL::4 2 4.00000
+DEAL::4 3 0
+DEAL::4 4 8.00000
+DEAL::OK
+DEAL::Chunk size: 5
+DEAL::0 0 0
+DEAL::0 1 0
+DEAL::0 2 0
+DEAL::0 3 0
+DEAL::0 4 0
+DEAL::1 0 0
+DEAL::1 1 0.500000
+DEAL::1 2 1.00000
+DEAL::1 3 0
+DEAL::1 4 0
+DEAL::2 0 0
+DEAL::2 1 0
+DEAL::2 2 2.00000
+DEAL::2 3 3.00000
+DEAL::2 4 0
+DEAL::3 0 0
+DEAL::3 1 1.50000
+DEAL::3 2 0
+DEAL::3 3 4.50000
+DEAL::3 4 6.00000
+DEAL::4 0 0
+DEAL::4 1 0
+DEAL::4 2 4.00000
+DEAL::4 3 0
+DEAL::4 4 8.00000
+DEAL::OK
+DEAL::Chunk size: 7
+DEAL::0 0 0
+DEAL::0 1 0
+DEAL::0 2 0
+DEAL::0 3 0
+DEAL::0 4 0
+DEAL::1 0 0
+DEAL::1 1 0.500000
+DEAL::1 2 1.00000
+DEAL::1 3 0
+DEAL::1 4 0
+DEAL::2 0 0
+DEAL::2 1 0
+DEAL::2 2 2.00000
+DEAL::2 3 3.00000
+DEAL::2 4 0
+DEAL::3 0 0
+DEAL::3 1 1.50000
+DEAL::3 2 0
+DEAL::3 3 4.50000
+DEAL::3 4 6.00000
+DEAL::4 0 0
+DEAL::4 1 0
+DEAL::4 2 4.00000
+DEAL::4 3 0
+DEAL::4 4 8.00000
+DEAL::OK
--- /dev/null
+//---------------------------- chunk_sparse_matrix_iterator_09.cc ---------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 2004, 2005, 2013 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//---------------------------- chunk_sparse_matrix_iterator_09.cc ---------------------------
+
+
+// this test is sparse_matrix_iterator_09 for a ChunkSparseMatrix
+
+#include "../tests.h"
+#include <deal.II/lac/chunk_sparse_matrix.h>
+#include <fstream>
+#include <iomanip>
+
+
+void test (const unsigned int chunk_size)
+{
+ deallog << "Chunk size: " << chunk_size << std::endl;
+
+ // create a sparsity pattern with totally
+ // empty lines (not even diagonals, since
+ // not quadratic)
+ ChunkSparsityPattern sparsity(4,5,1,chunk_size);
+ sparsity.add (1,1);
+ sparsity.add (3,1);
+ sparsity.compress ();
+
+ // attach a sparse matrix to it
+ ChunkSparseMatrix<double> A(sparsity);
+
+ // and loop over the elements of it
+ for (ChunkSparseMatrix<double>::const_iterator k=A.begin();
+ k!=A.end(); ++k)
+ deallog << k->row() << ' ' << k->column() << ' ' << k->value()
+ << std::endl;
+}
+
+
+
+int main ()
+{
+ std::ofstream logfile("chunk_sparse_matrix_iterator_09/output");
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ try
+ {
+ test (1);
+ test (3);
+ }
+ catch (std::exception &exc)
+ {
+ deallog << std::endl << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ deallog << "Exception on processing: " << std::endl
+ << exc.what() << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+
+ return 1;
+ }
+ catch (...)
+ {
+ deallog << std::endl << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ deallog << "Unknown exception!" << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ return 1;
+ };
+}
--- /dev/null
+
+DEAL::Chunk size: 1
+DEAL::1 1 0
+DEAL::3 1 0
+DEAL::Chunk size: 3
+DEAL::0 0 0
+DEAL::0 1 0
+DEAL::0 2 0
+DEAL::1 0 0
+DEAL::1 1 0
+DEAL::1 2 0
+DEAL::2 0 0
+DEAL::2 1 0
+DEAL::2 2 0
+DEAL::3 3 0
+DEAL::3 4 0
+DEAL::3 0 0
+DEAL::3 1 0
+DEAL::3 2 0
--- /dev/null
+//---------------------------- chunk_sparse_matrix_iterator_10.cc ---------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 2004, 2005, 2013 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//---------------------------- chunk_sparse_matrix_iterator_10.cc ---------------------------
+
+
+// this test is sparse_matrix_iterator_10 for a ChunkSparseMatrix and the same
+// test as chunk_sparse_matrix_iterator_09 with postfix operator++ instead of
+// prefix
+
+#include "../tests.h"
+#include <deal.II/lac/chunk_sparse_matrix.h>
+#include <fstream>
+#include <iomanip>
+
+
+void test (const unsigned int chunk_size)
+{
+ deallog << "Chunk size: " << chunk_size << std::endl;
+
+ // create a sparsity pattern with totally
+ // empty lines (not even diagonals, since
+ // not quadratic)
+ ChunkSparsityPattern sparsity(4,5,1,chunk_size);
+ sparsity.add (1,1);
+ sparsity.add (3,1);
+ sparsity.compress ();
+
+ // attach a sparse matrix to it
+ ChunkSparseMatrix<double> A(sparsity);
+
+ // and loop over the elements of it
+ for (ChunkSparseMatrix<double>::const_iterator k=A.begin();
+ k!=A.end(); ++k)
+ deallog << k->row() << ' ' << k->column() << ' ' << k->value()
+ << std::endl;
+}
+
+
+
+int main ()
+{
+ std::ofstream logfile("chunk_sparse_matrix_iterator_10/output");
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ try
+ {
+ test (1);
+ test (3);
+ }
+ catch (std::exception &exc)
+ {
+ deallog << std::endl << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ deallog << "Exception on processing: " << std::endl
+ << exc.what() << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+
+ return 1;
+ }
+ catch (...)
+ {
+ deallog << std::endl << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ deallog << "Unknown exception!" << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ return 1;
+ };
+}
--- /dev/null
+
+DEAL::Chunk size: 1
+DEAL::1 1 0
+DEAL::3 1 0
+DEAL::Chunk size: 3
+DEAL::0 0 0
+DEAL::0 1 0
+DEAL::0 2 0
+DEAL::1 0 0
+DEAL::1 1 0
+DEAL::1 2 0
+DEAL::2 0 0
+DEAL::2 1 0
+DEAL::2 2 0
+DEAL::3 3 0
+DEAL::3 4 0
+DEAL::3 0 0
+DEAL::3 1 0
+DEAL::3 2 0
--- /dev/null
+//---------------------------- chunk_sparse_matrix_iterator_11.cc ---------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 2004, 2005, 2013 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//---------------------------- chunk_sparse_matrix_iterator_11.cc ---------------------------
+
+
+// comparisons between chunk sparse matrix iterators, same as
+// sparse_matrix_iterators_11 otherwise.
+
+#include "../tests.h"
+#include <deal.II/lac/chunk_sparse_matrix.h>
+#include <fstream>
+#include <iomanip>
+
+
+void test (const unsigned int chunk_size)
+{
+ deallog << "Chunk size: " << chunk_size << std::endl;
+
+ // create a sparsity pattern with totally
+ // empty lines (not even diagonals, since
+ // not quadratic)
+ ChunkSparsityPattern sparsity(4,5,1,chunk_size);
+ sparsity.add (1,1);
+ sparsity.add (3,1);
+ sparsity.compress ();
+
+ // attach a sparse matrix to it
+ ChunkSparseMatrix<double> A(sparsity);
+
+ ChunkSparseMatrix<double>::iterator k = A.begin(),
+ j = ++A.begin();
+
+ Assert (k < j, ExcInternalError());
+ Assert (j > k, ExcInternalError());
+
+ Assert (!(j < k), ExcInternalError());
+ Assert (!(k > j), ExcInternalError());
+
+ Assert (k != j, ExcInternalError());
+ Assert (!(k == j), ExcInternalError());
+
+ Assert (k == k, ExcInternalError());
+ Assert (!(k != k), ExcInternalError());
+
+ deallog << "OK" << std::endl;
+}
+
+
+
+int main ()
+{
+ std::ofstream logfile("chunk_sparse_matrix_iterator_11/output");
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ try
+ {
+ test (1);
+ test (3);
+ }
+ catch (std::exception &exc)
+ {
+ deallog << std::endl << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ deallog << "Exception on processing: " << std::endl
+ << exc.what() << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+
+ return 1;
+ }
+ catch (...)
+ {
+ deallog << std::endl << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ deallog << "Unknown exception!" << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ return 1;
+ };
+}
--- /dev/null
+
+DEAL::Chunk size: 1
+DEAL::OK
+DEAL::Chunk size: 3
+DEAL::OK
--- /dev/null
+//---------------------------- chunk_sparse_matrix_iterator_12.cc ---------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 2004, 2005, 2013 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//---------------------------- chunk_sparse_matrix_iterator_12.cc ---------------------------
+
+
+// like chunk_sparse_matrix_iterator_11 but for const_iterator
+
+#include "../tests.h"
+#include <deal.II/lac/chunk_sparse_matrix.h>
+#include <fstream>
+#include <iomanip>
+
+
+void test (const unsigned int chunk_size)
+{
+ deallog << "Chunk size: " << chunk_size << std::endl;
+
+ // create a sparsity pattern with totally
+ // empty lines (not even diagonals, since
+ // not quadratic)
+ ChunkSparsityPattern sparsity(4,5,1,chunk_size);
+ sparsity.add (1,1);
+ sparsity.add (3,1);
+ sparsity.compress ();
+
+ // attach a sparse matrix to it
+ ChunkSparseMatrix<double> A(sparsity);
+
+ ChunkSparseMatrix<double>::const_iterator k = A.begin(),
+ j = ++A.begin();
+
+ Assert (k < j, ExcInternalError());
+ Assert (j > k, ExcInternalError());
+
+ Assert (!(j < k), ExcInternalError());
+ Assert (!(k > j), ExcInternalError());
+
+ Assert (k != j, ExcInternalError());
+ Assert (!(k == j), ExcInternalError());
+
+ Assert (k == k, ExcInternalError());
+ Assert (!(k != k), ExcInternalError());
+
+ deallog << "OK" << std::endl;
+}
+
+
+
+int main ()
+{
+ std::ofstream logfile("chunk_sparse_matrix_iterator_12/output");
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ try
+ {
+ test (1);
+ test (3);
+ }
+ catch (std::exception &exc)
+ {
+ deallog << std::endl << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ deallog << "Exception on processing: " << std::endl
+ << exc.what() << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+
+ return 1;
+ }
+ catch (...)
+ {
+ deallog << std::endl << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ deallog << "Unknown exception!" << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ return 1;
+ };
+}
--- /dev/null
+
+DEAL::Chunk size: 1
+DEAL::OK
+DEAL::Chunk size: 3
+DEAL::OK
--- /dev/null
+//---------------------------- chunk_sparse_matrix_iterator_13.cc ---------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 2004, 2005, 2013 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//---------------------------- chunk_sparse_matrix_iterator_13.cc ---------------------------
+
+
+// test ChunkSparseMatrix::iterator::operator-
+
+#include "../tests.h"
+#include <deal.II/lac/chunk_sparse_matrix.h>
+#include <fstream>
+#include <iomanip>
+
+
+void test (const unsigned int chunk_size)
+{
+ deallog << "Chunk size: " << chunk_size << std::endl;
+
+ // create a sparsity pattern with totally
+ // empty lines (not even diagonals, since
+ // not quadratic)
+ ChunkSparsityPattern sp(5,5,3,chunk_size);
+ for (unsigned int i=0; i<5; ++i)
+ for (unsigned int j=0; j<5; ++j)
+ if (((i+2*j+1) % 3 == 0)
+ ||
+ (i==j))
+ sp.add (i,j);
+ sp.compress ();
+
+ // attach a sparse matrix to it
+ ChunkSparseMatrix<double> m(sp);
+
+ for (unsigned int row=0; row<sp.n_rows(); ++row)
+ Assert (m.begin(row)-m.begin(row) == 0,
+ ExcInternalError());
+
+ for (unsigned int row=0; row<sp.n_rows(); ++row)
+ Assert (m.end(row)-m.begin(row) == (int)sp.row_length(row),
+ ExcInternalError());
+ for (unsigned int row=0; row<sp.n_rows(); ++row)
+ Assert (m.begin(row)-m.end(row) == -(int)sp.row_length(row),
+ ExcInternalError());
+
+ {
+ unsigned int counter = 0;
+ for (unsigned int row=0; row<sp.n_rows(); ++row)
+ {
+ Assert (m.begin(row)-m.begin(0) == (int)counter,
+ ExcInternalError());
+ Assert (m.begin(0)-m.begin(row) == -(int)counter,
+ ExcInternalError());
+ counter += sp.row_length(row);
+ }
+ }
+
+ Assert (m.begin() - m.begin(0) == 0, ExcInternalError());
+ Assert (m.begin(0) - m.begin() == 0, ExcInternalError());
+ Assert (m.end(sp.n_rows()-1) - m.end() == 0, ExcInternalError());
+ Assert (m.end() - m.end(sp.n_rows()-1) == 0, ExcInternalError());
+ Assert (m.end() - m.begin() == (int)sp.n_nonzero_elements(),
+ ExcInternalError());
+ Assert (m.begin() - m.end() == -(int)sp.n_nonzero_elements(),
+ ExcInternalError());
+
+ deallog << "OK" << std::endl;
+}
+
+
+
+int main ()
+{
+ std::ofstream logfile("chunk_sparse_matrix_iterator_13/output");
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ try
+ {
+ test (1);
+ test (2);
+ test (4);
+ test (7);
+ }
+ catch (std::exception &exc)
+ {
+ deallog << std::endl << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ deallog << "Exception on processing: " << std::endl
+ << exc.what() << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+
+ return 1;
+ }
+ catch (...)
+ {
+ deallog << std::endl << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ deallog << "Unknown exception!" << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ return 1;
+ };
+}
--- /dev/null
+
+DEAL::Chunk size: 1
+DEAL::OK
+DEAL::Chunk size: 2
+DEAL::OK
+DEAL::Chunk size: 4
+DEAL::OK
+DEAL::Chunk size: 7
+DEAL::OK
projects / dealii-svn.git / commitdiff