From: kronbichler Date: Tue, 7 May 2013 17:14:12 +0000 (+0000) Subject: Implement iterator for ChunkSparseMatrix. Parallelize matrix-vector product for Chunk... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=87028b3d0bade771727316fc2ec9e019dd47cee5;p=dealii-svn.git Implement iterator for ChunkSparseMatrix. Parallelize matrix-vector product for ChunkSparseMatrix. git-svn-id: https://svn.dealii.org/trunk@29476 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/include/deal.II/lac/chunk_sparse_matrix.h b/deal.II/include/deal.II/lac/chunk_sparse_matrix.h index 9c20d69fac..3d7a5c2fc3 100644 --- a/deal.II/include/deal.II/lac/chunk_sparse_matrix.h +++ b/deal.II/include/deal.II/lac/chunk_sparse_matrix.h @@ -29,6 +29,359 @@ template class FullMatrix; *@{ */ +/** + * A namespace in which we declare iterators over the elements of sparse + * matrices. + */ +namespace ChunkSparseMatrixIterators +{ + // forward declaration + template + 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 + 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 &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 + class Accessor : public ChunkSparsityPatternIterators::Accessor + { + public: + /** + * Typedef for the type (including constness) of the matrix to be used + * here. + */ + typedef const ChunkSparseMatrix 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 &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 + 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 + class Accessor : 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 number (i.e., most commonly a + * double) or to another object of type Reference. 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 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 + 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 + * Constness=false, this class is for iterators to constant + * matrices. + */ + template + class Iterator + { + public: + /** + * Typedef for the matrix type (including constness) we are to operate on. + */ + typedef + typename Accessor::MatrixType + MatrixType; + + /** + * A typedef for the type you get when you dereference an iterator + * of the current kind. + */ + typedef + const Accessor & 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 &i); + + /** + * 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 ==. + */ + 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 accessor; + }; + +} + + /** * Sparse matrix. This class implements the function to store values @@ -56,42 +409,46 @@ public: typedef number value_type; /** - * Declare a type that has holds - * real-valued numbers with the - * same precision as the template - * argument to this class. If the - * template argument of this - * class is a real data type, - * then real_type equals the - * template argument. If the - * template argument is a - * std::complex type then - * real_type equals the type - * underlying the complex - * numbers. + * Declare a type that has holds real-valued numbers with the same precision + * as the template argument to this class. If the template argument of this + * class is a real data type, then real_type equals the template + * argument. If the template argument is a std::complex type then real_type + * equals the type underlying the complex numbers. * - * This typedef is used to - * represent the return type of - * norms. + * This typedef is used to represent the return type of norms. */ typedef typename numbers::NumberTraits::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 + 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 + iterator; + + /** + * A structure that describes some of the traits of this class in terms of + * its run-time behavior. Some other classes (such as the block matrix + * classes) that take one or other of the matrix classes as its template + * parameters can tune their behavior based on the variables in this class. */ struct Traits { /** - * It is safe to elide additions of - * zeros to individual elements of - * this matrix. + * It is safe to elide additions of zeros to individual elements of this + * matrix. */ static const bool zero_addition_can_be_elided = true; }; @@ -101,172 +458,106 @@ public: */ //@{ /** - * Constructor; initializes the matrix to - * be empty, without any structure, i.e. - * the matrix is not usable at all. This - * constructor is therefore only useful - * for matrices which are members of a - * class. All other matrices should be - * created at a point in the data flow - * where all necessary information is - * available. + * Constructor; initializes the matrix to be empty, without any structure, + * i.e. the matrix is not usable at all. This constructor is therefore only + * useful for matrices which are members of a class. All other matrices + * should be created at a point in the data flow where all necessary + * information is available. * - * You have to initialize - * the matrix before usage with - * reinit(const ChunkSparsityPattern&). + * You have to initialize the matrix before usage with reinit(const + * ChunkSparsityPattern&). */ ChunkSparseMatrix (); /** - * Copy constructor. This constructor is - * only allowed to be called if the matrix - * to be copied is empty. This is for the - * same reason as for the - * ChunkSparsityPattern, see there for the - * details. + * Copy constructor. This constructor is only allowed to be called if the + * matrix to be copied is empty. This is for the same reason as for the + * ChunkSparsityPattern, see there for the details. * - * If you really want to copy a whole - * matrix, you can do so by using the + * If you really want to copy a whole matrix, you can do so by using the * copy_from() function. */ ChunkSparseMatrix (const ChunkSparseMatrix &); /** - * Constructor. Takes the given - * matrix sparsity structure to - * represent the sparsity pattern - * of this matrix. You can change - * the sparsity pattern later on - * by calling the reinit(const - * ChunkSparsityPattern&) function. + * Constructor. Takes the given matrix sparsity structure to represent the + * sparsity pattern of this matrix. You can change the sparsity pattern + * later on by calling the reinit(const ChunkSparsityPattern&) function. * - * You have to make sure that the - * lifetime of the sparsity - * structure is at least as long - * as that of this matrix or as - * long as reinit(const - * ChunkSparsityPattern&) is not - * called with a new sparsity - * pattern. + * You have to make sure that the lifetime of the sparsity structure is at + * least as long as that of this matrix or as long as reinit(const + * ChunkSparsityPattern&) is not called with a new sparsity pattern. * - * The constructor is marked - * explicit so as to disallow - * that someone passes a sparsity - * pattern in place of a sparse - * matrix to some function, where - * an empty matrix would be - * generated then. + * The constructor is marked explicit so as to disallow that someone passes + * a sparsity pattern in place of a sparse matrix to some function, where an + * empty matrix would be generated then. */ explicit ChunkSparseMatrix (const ChunkSparsityPattern &sparsity); /** - * Copy constructor: initialize - * the matrix with the identity - * matrix. This constructor will - * throw an exception if the - * sizes of the sparsity pattern - * and the identity matrix do not - * coincide, or if the sparsity - * pattern does not provide for - * nonzero entries on the entire - * diagonal. + * Copy constructor: initialize the matrix with the identity matrix. This + * constructor will throw an exception if the sizes of the sparsity pattern + * and the identity matrix do not coincide, or if the sparsity pattern does + * not provide for nonzero entries on the entire diagonal. */ ChunkSparseMatrix (const ChunkSparsityPattern &sparsity, const IdentityMatrix &id); /** - * Destructor. Free all memory, but do not - * release the memory of the sparsity - * structure. + * Destructor. Free all memory, but do not release the memory of the + * sparsity structure. */ virtual ~ChunkSparseMatrix (); /** - * Copy operator. Since copying - * entire sparse matrices is a - * very expensive operation, we - * disallow doing so except for - * the special case of empty - * matrices of size zero. This - * doesn't seem particularly - * useful, but is exactly what - * one needs if one wanted to - * have a - * std::vector@ - * @>: 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 + * std::vector@ @>: 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 &operator = (const ChunkSparseMatrix &); /** - * 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 & 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 - * matrix=0, 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 matrix=0, which sets + * all elements of the matrix to zero, but keep the sparsity pattern * previously used. */ ChunkSparseMatrix &operator = (const double d); /** - * Reinitialize the sparse matrix - * with the given sparsity - * pattern. The latter tells the - * matrix how many nonzero - * elements there need to be + * Reinitialize the sparse matrix with the given sparsity pattern. The + * latter tells the matrix how many nonzero elements there need to be * reserved. * - * Regarding memory allocation, - * the same applies as said - * above. + * Regarding memory allocation, the same applies as said above. * - * You have to make sure that the - * lifetime of the sparsity - * structure is at least as long - * as that of this matrix or as - * long as reinit(const - * ChunkSparsityPattern &) is not - * called with a new sparsity - * structure. + * You have to make sure that the lifetime of the sparsity structure is at + * least as long as that of this matrix or as long as reinit(const + * ChunkSparsityPattern &) is not called with a new sparsity structure. * - * The elements of the matrix are - * set to zero by this function. + * The elements of the matrix are set to zero by this function. */ virtual void reinit (const ChunkSparsityPattern &sparsity); /** - * Release all memory and return - * to a state just like after - * having called the default - * constructor. It also forgets - * the sparsity pattern it was + * Release all memory and return to a state just like after having called + * the default constructor. It also forgets the sparsity pattern it was * previously tied to. */ virtual void clear (); @@ -276,74 +567,52 @@ public: */ //@{ /** - * Return whether the object is - * empty. It is empty if either - * both dimensions are zero or no - * ChunkSparsityPattern is - * associated. + * Return whether the object is empty. It is empty if either both dimensions + * are zero or no ChunkSparsityPattern is associated. */ bool empty () const; /** - * Return the dimension of the - * image space. To remember: the - * matrix is of dimension - * $m \times n$. + * Return the dimension of the image space. To remember: the matrix is of + * dimension $m \times n$. */ unsigned int m () const; /** - * Return the dimension of the - * range space. To remember: the - * matrix is of dimension - * $m \times n$. + * Return the dimension of the range space. To remember: the matrix is of + * dimension $m \times n$. */ unsigned int n () const; /** - * Return the number of nonzero - * elements of this - * matrix. Actually, it returns - * the number of entries in the - * sparsity pattern; if any of - * the entries should happen to - * be zero, it is counted anyway. + * Return the number of nonzero elements of this matrix. Actually, it + * returns the number of entries in the sparsity pattern; if any of the + * entries should happen to be zero, it is counted anyway. */ unsigned int n_nonzero_elements () const; /** - * Return the number of actually - * nonzero elements of this - * matrix. + * Return the number of actually nonzero elements of this matrix. * - * Note, that this function does - * (in contrary to - * n_nonzero_elements()) not - * count all entries of the - * sparsity pattern but only the - * ones that are nonzero. + * Note, that this function does (in contrary to n_nonzero_elements()) not + * count all entries of the sparsity pattern but only the ones that are + * nonzero. */ unsigned int n_actually_nonzero_elements () const; /** - * Return a (constant) reference - * to the underlying sparsity - * pattern of this matrix. + * Return a (constant) reference to the underlying sparsity pattern of this + * matrix. * - * Though the return value is - * declared const, 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 const, you should be aware + * that it may change if you call any nonconstant function of objects which + * operate on it. */ const ChunkSparsityPattern &get_sparsity_pattern () const; /** - * Determine an estimate for the - * memory consumption (in bytes) - * of this object. See - * MemoryConsumption. + * Determine an estimate for the memory consumption (in bytes) of this + * object. See MemoryConsumption. */ std::size_t memory_consumption () const; @@ -353,47 +622,31 @@ public: */ //@{ /** - * Set the element (i,j) - * to value. Throws an - * error if the entry does not - * exist or if value is - * not a finite number. Still, it - * is allowed to store zero - * values in non-existent fields. + * Set the element (i,j) to value. Throws an error if the + * entry does not exist or if value 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 value to the - * element (i,j). Throws - * an error if the entry does not - * exist or if value is - * not a finite number. Still, it - * is allowed to store zero - * values in non-existent fields. + * Add value to the element (i,j). Throws an error if the + * entry does not exist or if value 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 - * values in the given - * global matrix row at columns - * specified by col_indices in the - * sparse matrix. + * Add an array of values given by values in the given global + * matrix row at columns specified by col_indices in the sparse matrix. * - * The optional parameter - * elide_zero_values 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 true, - * i.e., zero values won't be added - * into the matrix. + * The optional parameter elide_zero_values 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 true, + * i.e., zero values won't be added into the matrix. */ template void add (const unsigned int row, @@ -404,140 +657,84 @@ public: const bool col_indices_are_sorted = false); /** - * Multiply the entire matrix by a - * fixed factor. + * Multiply the entire matrix by a fixed factor. */ ChunkSparseMatrix &operator *= (const number factor); /** - * Divide the entire matrix by a - * fixed factor. + * Divide the entire matrix by a fixed factor. */ ChunkSparseMatrix &operator /= (const number factor); /** - * Symmetrize the matrix by - * forming the mean value between - * the existing matrix and its - * transpose, $A = \frac 12(A+A^T)$. + * Symmetrize the matrix by forming the mean value between the existing + * matrix and its transpose, $A = \frac 12(A+A^T)$. * - * This operation assumes that - * the underlying sparsity - * pattern represents a symmetric - * object. If this is not the - * case, then the result of this - * operation will not be a - * symmetric matrix, since it - * only explicitly symmetrizes - * by looping over the lower left - * triangular part for efficiency - * reasons; if there are entries - * in the upper right triangle, - * then these elements are missed - * in the - * symmetrization. Symmetrization - * of the sparsity pattern can be - * obtain by - * ChunkSparsityPattern::symmetrize(). + * This operation assumes that the underlying sparsity pattern represents a + * symmetric object. If this is not the case, then the result of this + * operation will not be a symmetric matrix, since it only explicitly + * symmetrizes by looping over the lower left triangular part for efficiency + * reasons; if there are entries in the upper right triangle, then these + * elements are missed in the symmetrization. Symmetrization of the sparsity + * pattern can be obtain by ChunkSparsityPattern::symmetrize(). */ void symmetrize (); /** - * Copy the given matrix to this - * one. The operation throws an - * error if the sparsity patterns - * of the two involved matrices - * do not point to the same - * object, since in this case the - * copy operation is - * cheaper. Since this operation - * is notheless not for free, we - * do not make it available - * through operator =, - * 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 + * operator =, 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 - * *this. + * The function returns a reference to *this. */ template ChunkSparseMatrix & copy_from (const ChunkSparseMatrix &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 std::pair, where - * value needs to be - * convertible to the element - * type of this class, as - * specified by the - * number template + * The only difference to the cited function is that the objects which the + * inner iterator points to need to be of type std::pair, where value needs to be convertible to the element + * type of this class, as specified by the number 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 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 void copy_from (const FullMatrix &matrix); /** - * Add matrix scaled by - * factor to this matrix, - * i.e. the matrix factor*matrix - * is added to this. 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 matrix scaled by factor to this matrix, i.e. the + * matrix factor*matrix is added to this. 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 @@ -551,54 +748,31 @@ public: //@{ /** - * Return the value of the entry - * (i,j). 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,j). 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,j). 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,j). The only difference is that if this entry + * does not exist in the sparsity pattern, then instead of raising an + * exception, zero is returned. While this may be convenient in some cases, + * note that it is simple to write algorithms that are slow compared to an + * optimal solution, since the sparsity of the matrix is not used. * - * If you are looping over all elements, - * consider using one of the iterator - * classes instead, since they are - * tailored better to a sparse matrix + * If you are looping over all elements, consider using one of the iterator + * classes instead, since they are tailored better to a sparse matrix * structure. */ number el (const unsigned int i, @@ -611,21 +785,16 @@ public: * error if the matrix is not * quadratic. * - * This function is considerably - * faster than the operator()(), - * since for quadratic matrices, the - * diagonal entry may be the - * first to be stored in each row - * and access therefore does not - * involve searching for the - * right column number. + * This function is considerably faster than the operator()(), since for + * quadratic matrices, the diagonal entry may be the first to be stored in + * each row and access therefore does not involve searching for the right + * column number. */ number diag_element (const unsigned int i) const; /** - * Same as above, but return a - * writeable reference. You're - * sure you know what you do? + * Same as above, but return a writeable reference. You're sure you know + * what you do? */ number &diag_element (const unsigned int i); @@ -635,166 +804,105 @@ public: */ //@{ /** - * Matrix-vector multiplication: - * let dst = M*src with - * M being this matrix. + * Matrix-vector multiplication: let dst = M*src with M 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 void vmult (OutVector &dst, const InVector &src) const; /** - * Matrix-vector multiplication: - * let dst = MT*src with - * M being this - * matrix. This function does the - * same as vmult() but takes - * the transposed matrix. + * Matrix-vector multiplication: let dst = MT*src with + * M 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 void Tvmult (OutVector &dst, const InVector &src) const; /** - * Adding Matrix-vector - * multiplication. Add - * M*src on dst - * with M being this - * matrix. + * Adding Matrix-vector multiplication. Add M*src on dst with + * M 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 void vmult_add (OutVector &dst, const InVector &src) const; /** - * Adding Matrix-vector - * multiplication. Add - * MT*src to - * dst with M being - * this matrix. This function - * does the same as vmult_add() - * but takes the transposed - * matrix. + * Adding Matrix-vector multiplication. Add MT*src to + * dst with M 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 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 somenumber matrix_norm_square (const Vector &v) const; /** - * Compute the matrix scalar - * product $\left(u,Mv\right)$. + * Compute the matrix scalar product $\left(u,Mv\right)$. */ template somenumber matrix_scalar_product (const Vector &u, const Vector &v) const; /** - * Compute the residual of an - * equation Mx=b, where - * the residual is defined to be - * r=b-Mx. Write the - * residual into - * dst. The - * l2 norm of - * the residual vector is - * returned. + * Compute the residual of an equation Mx=b, where the residual is + * defined to be r=b-Mx. Write the residual into dst. The + * l2 norm of the residual vector is returned. * - * Source x and destination - * dst must not be the same - * vector. + * Source x and destination dst must not be the same vector. */ template somenumber residual (Vector &dst, @@ -808,37 +916,26 @@ public: //@{ /** - * Return the l1-norm of the matrix, that is - * $|M|_1=max_{all columns j}\sum_{all - * rows i} |M_ij|$, - * (max. sum of columns). - * This is the - * natural matrix norm that is compatible - * to the l1-norm for vectors, i.e. - * $|Mv|_1\leq |M|_1 |v|_1$. - * (cf. Haemmerlin-Hoffmann : Numerische Mathematik) + * Return the l1-norm of the matrix, that is $|M|_1=max_{all columns + * j}\sum_{all rows i} |M_ij|$, (max. sum of columns). This is the natural + * matrix norm that is compatible to the l1-norm for vectors, i.e. + * $|Mv|_1\leq |M|_1 |v|_1$. (cf. Haemmerlin-Hoffmann : Numerische + * Mathematik) */ real_type l1_norm () const; /** - * Return the linfty-norm of the - * matrix, that is - * $|M|_infty=max_{all rows i}\sum_{all - * columns j} |M_ij|$, - * (max. sum of rows). - * This is the - * natural matrix norm that is compatible - * to the linfty-norm of vectors, i.e. - * $|Mv|_infty \leq |M|_infty |v|_infty$. - * (cf. Haemmerlin-Hoffmann : Numerische Mathematik) + * Return the linfty-norm of the matrix, that is $|M|_infty=max_{all rows + * i}\sum_{all columns j} |M_ij|$, (max. sum of rows). This is the natural + * matrix norm that is compatible to the linfty-norm of vectors, i.e. + * $|Mv|_infty \leq |M|_infty |v|_infty$. (cf. Haemmerlin-Hoffmann : + * Numerische Mathematik) */ real_type linfty_norm () const; /** - * Return the frobenius norm of the - * matrix, i.e. the square root of the - * sum of squares of all entries in the - * matrix. + * Return the frobenius norm of the matrix, i.e. the square root of the sum + * of squares of all entries in the matrix. */ real_type frobenius_norm () const; //@} @@ -848,14 +945,9 @@ public: //@{ /** - * Apply the Jacobi - * preconditioner, which - * multiplies every element of - * the src vector by the - * inverse of the respective - * diagonal element and - * multiplies the result with the - * relaxation factor omega. + * Apply the Jacobi preconditioner, which multiplies every element of the + * src vector by the inverse of the respective diagonal element and + * multiplies the result with the relaxation factor omega. */ template void precondition_Jacobi (Vector &dst, @@ -863,8 +955,7 @@ public: const number omega = 1.) const; /** - * Apply SSOR preconditioning to - * src. + * Apply SSOR preconditioning to src. */ template void precondition_SSOR (Vector &dst, @@ -872,8 +963,7 @@ public: const number om = 1.) const; /** - * Apply SOR preconditioning - * matrix to src. + * Apply SOR preconditioning matrix to src. */ template void precondition_SOR (Vector &dst, @@ -881,9 +971,7 @@ public: const number om = 1.) const; /** - * Apply transpose SOR - * preconditioning matrix to - * src. + * Apply transpose SOR preconditioning matrix to src. */ template void precondition_TSOR (Vector &dst, @@ -891,11 +979,8 @@ public: const number om = 1.) const; /** - * Perform SSOR preconditioning - * in-place. Apply the - * preconditioner matrix without - * copying to a second vector. - * omega is the relaxation + * Perform SSOR preconditioning in-place. Apply the preconditioner matrix + * without copying to a second vector. omega is the relaxation * parameter. */ template @@ -903,18 +988,15 @@ public: const number omega = 1.) const; /** - * Perform an SOR preconditioning - * in-place. omega is - * the relaxation parameter. + * Perform an SOR preconditioning in-place. omega is the + * relaxation parameter. */ template void SOR (Vector &v, const number om = 1.) const; /** - * Perform a transpose SOR - * preconditioning in-place. - * omega is the + * Perform a transpose SOR preconditioning in-place. omega is the * relaxation parameter. */ template @@ -922,21 +1004,14 @@ public: const number om = 1.) const; /** - * Perform a permuted SOR - * preconditioning in-place. + * Perform a permuted SOR preconditioning in-place. * - * The standard SOR method is - * applied in the order - * prescribed by permutation, - * that is, first the row - * permutation[0], then - * permutation[1] 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 + * permutation, that is, first the row permutation[0], + * then permutation[1] and so on. For efficiency reasons, the + * permutation as well as its inverse are required. * - * omega is the - * relaxation parameter. + * omega is the relaxation parameter. */ template void PSOR (Vector &v, @@ -945,24 +1020,14 @@ public: const number om = 1.) const; /** - * Perform a transposed permuted SOR - * preconditioning in-place. + * Perform a transposed permuted SOR preconditioning in-place. * - * The transposed SOR method is - * applied in the order - * prescribed by - * permutation, that is, - * first the row - * permutation[m()-1], - * then - * permutation[m()-2] - * 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 + * permutation, that is, first the row permutation[m()-1], + * then permutation[m()-2] and so on. For efficiency reasons, the + * permutation as well as its inverse are required. * - * omega is the - * relaxation parameter. + * omega is the relaxation parameter. */ template void TPSOR (Vector &v, @@ -971,10 +1036,8 @@ public: const number om = 1.) const; /** - * Do one SOR step on v. - * Performs a direct SOR step - * with right hand side - * b. + * Do one SOR step on v. Performs a direct SOR step with right + * hand side b. */ template void SOR_step (Vector &v, @@ -982,10 +1045,8 @@ public: const number om = 1.) const; /** - * Do one adjoint SOR step on - * v. Performs a direct - * TSOR step with right hand side - * b. + * Do one adjoint SOR step on v. Performs a direct TSOR step with + * right hand side b. */ template void TSOR_step (Vector &v, @@ -993,16 +1054,124 @@ public: const number om = 1.) const; /** - * Do one SSOR step on - * v. Performs a direct - * SSOR step with right hand side - * b by performing TSOR - * after SOR. + * Do one SSOR step on v. Performs a direct SSOR step with right + * hand side b by performing TSOR after SOR. */ template void SSOR_step (Vector &v, const Vector &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 r. 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 + * end(r). 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 r. 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 r. 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 + * end(r). 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 r. It points to the first element past the + * end of line @p r, or past the end of the entire sparsity pattern. This is + * the version for non-constant matrices. + * + * Note that the end iterator is not necessarily dereferencable. This is in + * particular the case if it is the end iterator for the last row of a + * matrix. + * + * Note that due to the layout in ChunkSparseMatrix, iterating over matrix + * entries is considerably slower than for a sparse matrix, as the iterator + * is travels row-by-row, whereas data is stored in chunks of several rows + * and columns. + */ + iterator end (const unsigned int r); //@} /** * @name Input/Output @@ -1010,54 +1179,30 @@ public: //@{ /** - * Print the matrix to the given - * stream, using the format - * (line,col) value, - * i.e. one nonzero entry of the - * matrix per line. + * Print the matrix to the given stream, using the format (line,col) + * value, 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: precision and - * scientific are used - * to determine the number - * format, where scientific = - * false means fixed point - * notation. A zero entry for - * width 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: + * precision and scientific are used to determine the + * number format, where scientific = false means fixed point + * notation. A zero entry for width 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 large amounts - * of output if applied to a - * large matrix! + * @attention This function may produce large amounts of output if + * applied to a large matrix! */ void print_formatted (std::ostream &out, const unsigned int precision = 3, @@ -1067,69 +1212,40 @@ public: const double denominator = 1.) const; /** - * Print the actual pattern of - * the matrix. For each entry - * with an absolute value larger - * than threshold, a '*' is - * printed, a ':' for every value - * smaller and a '.' for every - * entry not allocated. + * Print the actual pattern of the matrix. For each entry with an absolute + * value larger than threshold, a '*' is printed, a ':' for every value + * smaller and a '.' for every entry not allocated. */ void print_pattern(std::ostream &out, const double threshold = 0.) const; /** - * Write the data of this object - * en bloc to a file. This is - * done in a binary mode, so the - * output is neither readable by - * humans nor (probably) by other - * computers using a different - * operating system of number - * format. + * Write the data of this object en bloc to a file. This is done in a binary + * mode, so the output is neither readable by humans nor (probably) by other + * computers using a different operating system or number format. * - * The purpose of this function - * is that you can swap out - * matrices and sparsity pattern - * if you are short of memory, - * want to communicate between - * different programs, or allow - * objects to be persistent - * across different runs of the - * program. + * The purpose of this function is that you can swap out matrices and + * sparsity pattern if you are short of memory, want to communicate between + * different programs, or allow objects to be persistent across different + * runs of the program. */ void block_write (std::ostream &out) const; /** - * Read data that has previously - * been written by block_write() - * from a file. This is done - * using the inverse operations - * to the above function, so it - * is reasonably fast because the - * bitstream is not interpreted - * except for a few numbers up - * front. + * Read data that has previously been written by block_write() from a + * file. This is done using the inverse operations to the above function, so + * it is reasonably fast because the bitstream is not interpreted except for + * a few numbers up front. * - * The object is resized on this - * operation, and all previous - * contents are lost. Note, - * however, that no checks are - * performed whether new data and - * the underlying ChunkSparsityPattern - * object fit together. It is - * your responsibility to make - * sure that the sparsity pattern - * and the data to be read match. + * The object is resized on this operation, and all previous contents are + * lost. Note, however, that no checks are performed whether new data and + * the underlying ChunkSparsityPattern object fit together. It is your + * responsibility to make sure that the sparsity pattern and the data to be + * read match. * - * A primitive form of error - * checking is performed which - * will recognize the bluntest - * attempts to interpret some - * data as a matrix stored - * bitwise to a file that wasn't - * actually created that way, but - * not more. + * A primitive form of error checking is performed which will recognize the + * bluntest attempts to interpret some data as a matrix stored bitwise to a + * file that wasn't actually created that way, but not more. */ void block_read (std::istream &in); //@} @@ -1167,53 +1283,43 @@ public: //@} private: /** - * Pointer to the sparsity - * pattern used for this - * matrix. In order to guarantee - * that it is not deleted while - * still in use, we subscribe to - * it using the SmartPointer - * class. + * Pointer to the sparsity pattern used for this matrix. In order to + * guarantee that it is not deleted while still in use, we subscribe to it + * using the SmartPointer class. */ SmartPointer > 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 friend class ChunkSparseMatrix; + + /** + * Also give access to internal details to the iterator/accessor + * classes. + */ + template friend class ChunkSparseMatrixIterators::Iterator; + template friend class ChunkSparseMatrixIterators::Accessor; }; /*@}*/ @@ -1242,6 +1348,17 @@ unsigned int ChunkSparseMatrix::n () const +template +inline +const ChunkSparsityPattern & +ChunkSparseMatrix::get_sparsity_pattern () const +{ + Assert (cols != 0, ExcNotInitialized()); + return *cols; +} + + + template inline unsigned int @@ -1250,7 +1367,7 @@ ChunkSparseMatrix::compute_location (const unsigned int i, { const unsigned int chunk_size = cols->get_chunk_size(); const unsigned int chunk_index - = cols->sparsity_pattern(i/chunk_size, j/chunk_size); + = cols->sparsity_pattern(i/chunk_size, j/chunk_size); if (chunk_index == ChunkSparsityPattern::invalid_entry) return ChunkSparsityPattern::invalid_entry; @@ -1275,10 +1392,8 @@ void ChunkSparseMatrix::set (const unsigned int i, Assert (numbers::is_finite(value), ExcNumberNotFinite()); Assert (cols != 0, ExcNotInitialized()); - // it is allowed to set elements of - // the matrix that are not part of - // the sparsity pattern, if the - // value to which we set it is zero + // it is allowed to set elements of the matrix that are not part of the + // sparsity pattern, if the value to which we set it is zero const unsigned int index = compute_location(i,j); Assert ((index != SparsityPattern::invalid_entry) || (value == 0.), @@ -1340,14 +1455,10 @@ ChunkSparseMatrix::operator *= (const number factor) const unsigned int chunk_size = cols->get_chunk_size(); - // multiply all elements of the matrix with - // the given factor. this includes the - // padding elements in chunks that overlap - // the boundaries of the actual matrix -- - // but since multiplication with a number - // does not violate the invariant of - // keeping these elements at zero nothing - // can happen + // multiply all elements of the matrix with the given factor. this includes + // the padding elements in chunks that overlap the boundaries of the actual + // matrix -- but since multiplication with a number does not violate the + // invariant of keeping these elements at zero nothing can happen number *val_ptr = val; const number *const end_ptr = val + cols->sparsity_pattern.n_nonzero_elements() @@ -1374,14 +1485,10 @@ ChunkSparseMatrix::operator /= (const number factor) const unsigned int chunk_size = cols->get_chunk_size(); - // multiply all elements of the matrix with - // the given factor. this includes the - // padding elements in chunks that overlap - // the boundaries of the actual matrix -- - // but since multiplication with a number - // does not violate the invariant of - // keeping these elements at zero nothing - // can happen + // multiply all elements of the matrix with the given factor. this includes + // the padding elements in chunks that overlap the boundaries of the actual + // matrix -- but since multiplication with a number does not violate the + // invariant of keeping these elements at zero nothing can happen number *val_ptr = val; const number *const end_ptr = val + cols->sparsity_pattern.n_nonzero_elements() @@ -1433,31 +1540,7 @@ number ChunkSparseMatrix::diag_element (const unsigned int i) const Assert (m() == n(), ExcNotQuadratic()); Assert (iget_chunk_size(); - return val[cols->sparsity_pattern.rowstart[i/chunk_size] - * - chunk_size * chunk_size - + - (i % chunk_size) * chunk_size - + - (i % chunk_size)]; -} - - - -template -inline -number &ChunkSparseMatrix::diag_element (const unsigned int i) -{ - Assert (cols != 0, ExcNotInitialized()); - Assert (m() == n(), ExcNotQuadratic()); - Assert (iget_chunk_size(); return val[cols->sparsity_pattern.rowstart[i/chunk_size] @@ -1473,6 +1556,7 @@ number &ChunkSparseMatrix::diag_element (const unsigned int i) template template +inline void ChunkSparseMatrix::copy_from (const ForwardIterator begin, const ForwardIterator end) @@ -1480,9 +1564,8 @@ ChunkSparseMatrix::copy_from (const ForwardIterator begin, Assert (static_cast(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::value_type::const_iterator inner_iterator; unsigned int row=0; for (ForwardIterator i=begin; i!=end; ++i, ++row) @@ -1495,6 +1578,466 @@ ChunkSparseMatrix::copy_from (const ForwardIterator begin, } + +//--------------------------------------------------------------------------- + + +namespace ChunkSparseMatrixIterators +{ + template + inline + Accessor:: + Accessor (const MatrixType *matrix, + const unsigned int row) + : + ChunkSparsityPatternIterators::Accessor (&matrix->get_sparsity_pattern(), + row), + matrix (matrix) + {} + + + + template + inline + Accessor:: + Accessor (const MatrixType *matrix) + : + ChunkSparsityPatternIterators::Accessor (&matrix->get_sparsity_pattern()), + matrix (matrix) + {} + + + + template + inline + Accessor:: + Accessor (const ChunkSparseMatrixIterators::Accessor &a) + : + ChunkSparsityPatternIterators::Accessor (a), + matrix (&a.get_matrix()) + {} + + + + template + inline + number + Accessor::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 + inline + typename Accessor::MatrixType & + Accessor::get_matrix () const + { + return *matrix; + } + + + + template + inline + Accessor::Reference::Reference ( + const Accessor *accessor, + const bool) + : + accessor (accessor) + {} + + + template + inline + Accessor::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 + inline + const typename Accessor::Reference & + Accessor::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 + inline + const typename Accessor::Reference & + Accessor::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 + inline + const typename Accessor::Reference & + Accessor::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 + inline + const typename Accessor::Reference & + Accessor::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 + inline + const typename Accessor::Reference & + Accessor::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 + inline + Accessor:: + Accessor (MatrixType *matrix, + const unsigned int row) + : + ChunkSparsityPatternIterators::Accessor (&matrix->get_sparsity_pattern(), + row), + matrix (matrix) + {} + + + + template + inline + Accessor:: + Accessor (MatrixType *matrix) + : + ChunkSparsityPatternIterators::Accessor (&matrix->get_sparsity_pattern()), + matrix (matrix) + {} + + + + template + inline + typename Accessor::Reference + Accessor::value() const + { + return Reference(this,true); + } + + + + + template + inline + typename Accessor::MatrixType & + Accessor::get_matrix () const + { + return *matrix; + } + + + + template + inline + Iterator:: + Iterator (MatrixType *matrix, + const unsigned int row) + : + accessor(matrix, row) + {} + + + + template + inline + Iterator:: + Iterator (MatrixType *matrix) + : + accessor(matrix) + {} + + + + template + inline + Iterator:: + Iterator (const ChunkSparseMatrixIterators::Iterator &i) + : + accessor(*i) + {} + + + + template + inline + Iterator & + Iterator::operator++ () + { + accessor.advance (); + return *this; + } + + + template + inline + Iterator + Iterator::operator++ (int) + { + const Iterator iter = *this; + accessor.advance (); + return iter; + } + + + template + inline + const Accessor & + Iterator::operator* () const + { + return accessor; + } + + + template + inline + const Accessor * + Iterator::operator-> () const + { + return &accessor; + } + + + template + inline + bool + Iterator:: + operator == (const Iterator &other) const + { + return (accessor == other.accessor); + } + + + template + inline + bool + Iterator:: + operator != (const Iterator &other) const + { + return ! (*this == other); + } + + + template + inline + bool + Iterator:: + operator < (const Iterator &other) const + { + Assert (&accessor.get_matrix() == &other.accessor.get_matrix(), + ExcInternalError()); + + return (accessor < other.accessor); + } + + + template + inline + bool + Iterator:: + operator > (const Iterator &other) const + { + return (other < *this); + } + + + template + inline + int + Iterator:: + 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 + inline + Iterator + Iterator:: + operator + (const unsigned int n) const + { + Iterator x = *this; + for (unsigned int i=0; i +inline +typename ChunkSparseMatrix::const_iterator +ChunkSparseMatrix::begin () const +{ + return const_iterator(this, 0); +} + + +template +inline +typename ChunkSparseMatrix::const_iterator +ChunkSparseMatrix::end () const +{ + return const_iterator(this); +} + + +template +inline +typename ChunkSparseMatrix::iterator +ChunkSparseMatrix::begin () +{ + return iterator(this, 0); +} + + +template +inline +typename ChunkSparseMatrix::iterator +ChunkSparseMatrix::end () +{ + return iterator(this); +} + + +template +inline +typename ChunkSparseMatrix::const_iterator +ChunkSparseMatrix::begin (const unsigned int r) const +{ + Assert (r +inline +typename ChunkSparseMatrix::const_iterator +ChunkSparseMatrix::end (const unsigned int r) const +{ + Assert (r +inline +typename ChunkSparseMatrix::iterator +ChunkSparseMatrix::begin (const unsigned int r) +{ + Assert (r +inline +typename ChunkSparseMatrix::iterator +ChunkSparseMatrix::end (const unsigned int r) +{ + Assert (r +#include #include #include #include @@ -43,11 +44,9 @@ namespace internal namespace ChunkSparseMatrix { /** - * Add the result of multiplying a chunk - * of size chunk_size times chunk_size by - * a source vector fragment of size - * chunk_size to the destination vector - * fragment. + * Add the result of multiplying a chunk of size chunk_size times + * chunk_size by a source vector fragment of size chunk_size to the + * destination vector fragment. */ template + 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 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::~ChunkSparseMatrix () +namespace internal +{ + namespace ChunkSparseMatrix + { + template + void zero_subrange (const unsigned int begin, + const unsigned int end, + T *dst) + { + std::memset (dst+begin,0,(end-begin)*sizeof(T)); + } + } +} + + + template ChunkSparseMatrix & ChunkSparseMatrix::operator = (const double d) @@ -266,14 +377,27 @@ ChunkSparseMatrix::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, + std_cxx1x::_1, std_cxx1x::_2, + val), + grain_size); + else if (matrix_size > 0) + std::memset (&val[0], 0, matrix_size*sizeof(number)); return *this; } @@ -313,10 +437,8 @@ ChunkSparseMatrix::reinit (const ChunkSparsityPattern &sparsity) return; } - // allocate not just m() * n() elements but - // enough so that we can store full - // chunks. this entails some padding - // elements + // allocate not just m() * n() elements but enough so that we can store full + // chunks. this entails some padding elements const unsigned int chunk_size = cols->get_chunk_size(); const unsigned int N = cols->sparsity_pattern.n_nonzero_elements() * chunk_size * chunk_size; @@ -328,12 +450,10 @@ ChunkSparseMatrix::reinit (const ChunkSparsityPattern &sparsity) max_len = N; } - // fill with zeros. do not just fill N - // elements but all that we allocated to - // ensure that also the padding elements - // are zero and not left at previous values - if (val != 0) - std::fill_n (&val[0], max_len, 0); + // fill with zeros. do not just fill N elements but all that we allocated to + // ensure that also the padding elements are zero and not left at previous + // values + this->operator=(0.); } @@ -378,11 +498,9 @@ ChunkSparseMatrix::n_actually_nonzero_elements () const { Assert (cols != 0, ExcNotInitialized()); - // count those elements that are nonzero, - // even if they lie in the padding around - // the matrix. since we have the invariant - // that padding elements are zero, nothing - // bad can happen here + // count those elements that are nonzero, even if they lie in the padding + // around the matrix. since we have the invariant that padding elements are + // zero, nothing bad can happen here const unsigned int chunk_size = cols->get_chunk_size(); return std::count_if(&val[0], &val[cols->sparsity_pattern.n_nonzero_elements () * @@ -413,8 +531,7 @@ ChunkSparseMatrix::copy_from (const ChunkSparseMatrix &matri Assert (val != 0, ExcNotInitialized()); Assert (cols == matrix.cols, ExcDifferentChunkSparsityPatterns()); - // copy everything, including padding - // elements + // copy everything, including padding elements const unsigned int chunk_size = cols->get_chunk_size(); std::copy (&matrix.val[0], &matrix.val[cols->sparsity_pattern.n_nonzero_elements() @@ -453,8 +570,7 @@ ChunkSparseMatrix::add (const number factor, Assert (val != 0, ExcNotInitialized()); Assert (cols == matrix.cols, ExcDifferentChunkSparsityPatterns()); - // add everything, including padding - // elements + // add everything, including padding elements const unsigned int chunk_size = cols->get_chunk_size(); number *val_ptr = &val[0]; const somenumber *matrix_ptr = &matrix.val[0]; @@ -479,10 +595,8 @@ ChunkSparseMatrix::vmult (OutVector &dst, Assert (!PointerComparison::equal(&src, &dst), ExcSourceEqualsDestination()); - // set the output vector to zero and then - // add to it the contributions of vmults - // from individual chunks. this is what - // vmult_add does + // set the output vector to zero and then add to it the contributions of + // vmults from individual chunks. this is what vmult_add does dst = 0; vmult_add (dst, src); } @@ -509,10 +623,8 @@ ChunkSparseMatrix::Tvmult (OutVector &dst, Assert (!PointerComparison::equal(&src, &dst), ExcSourceEqualsDestination()); - // set the output vector to zero and then - // add to it the contributions of vmults - // from individual chunks. this is what - // vmult_add does + // set the output vector to zero and then add to it the contributions of + // vmults from individual chunks. this is what vmult_add does dst = 0; Tvmult_add (dst, src); } @@ -531,93 +643,18 @@ ChunkSparseMatrix::vmult_add (OutVector &dst, Assert(n() == src.size(), ExcDimensionMismatch(n(),src.size())); Assert (!PointerComparison::equal(&src, &dst), ExcSourceEqualsDestination()); + parallel::apply_to_subranges (0U, cols->sparsity_pattern.n_rows(), + std_cxx1x::bind (&internal::ChunkSparseMatrix::vmult_add_on_subrange + , + 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_rowsparsity_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; rchunk_size; ++r) - for (unsigned int c=0; cchunk_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; rchunk_size; ++r) - for (unsigned int c=0; cchunk_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; rchunk_size; ++r) - for (unsigned int c=0; cchunk_size; ++c) - dst(chunk_row * cols->chunk_size + r) - += (val_ptr[r*cols->chunk_size + c] * - src(*colnum_ptr * cols->chunk_size + c)); - - ++colnum_ptr; - val_ptr += cols->chunk_size * cols->chunk_size; - } - } } @@ -636,10 +673,8 @@ ChunkSparseMatrix::Tvmult_add (OutVector &dst, const unsigned int n_chunk_rows = cols->sparsity_pattern.n_rows(); - // loop over all chunks. note that we need - // to treat the last chunk row and column - // differently if they have padding - // elements + // loop over all chunks. note that we need to treat the last chunk row and + // column differently if they have padding elements const bool rows_have_padding = (m() % cols->chunk_size != 0), cols_have_padding = (n() % cols->chunk_size != 0); @@ -648,9 +683,8 @@ ChunkSparseMatrix::Tvmult_add (OutVector &dst, n_chunk_rows-1 : n_chunk_rows); - // like in vmult_add, but don't keep an - // iterator into dst around since we're not - // traversing it sequentially this time + // like in vmult_add, but don't keep an iterator into dst around since we're + // not traversing it sequentially this time const number *val_ptr = val; const unsigned int *colnum_ptr = cols->sparsity_pattern.colnums; @@ -670,8 +704,7 @@ ChunkSparseMatrix::Tvmult_add (OutVector &dst, src.begin() + chunk_row * cols->chunk_size, dst.begin() + *colnum_ptr * cols->chunk_size); else - // we're at a chunk column that - // has padding + // we're at a chunk column that has padding for (unsigned int r=0; rchunk_size; ++r) for (unsigned int c=0; cchunk_size; ++c) dst(*colnum_ptr * cols->chunk_size + c) @@ -683,8 +716,7 @@ ChunkSparseMatrix::Tvmult_add (OutVector &dst, } } - // now deal with last chunk row if - // necessary + // now deal with last chunk row if necessary if (rows_have_padding) { const unsigned int chunk_row = n_chunk_rows - 1; @@ -698,8 +730,7 @@ ChunkSparseMatrix::Tvmult_add (OutVector &dst, || (*colnum_ptr != cols->sparsity_pattern.n_cols()-1)) { - // we're at a chunk row but not - // column that has padding + // we're at a chunk row but not column that has padding for (unsigned int r=0; rchunk_size; ++r) for (unsigned int c=0; cchunk_size; ++c) dst(*colnum_ptr * cols->chunk_size + c) @@ -707,8 +738,7 @@ ChunkSparseMatrix::Tvmult_add (OutVector &dst, src(chunk_row * cols->chunk_size + r)); } else - // we're at a chunk row and - // column that has padding + // we're at a chunk row and column that has padding for (unsigned int r=0; rchunk_size; ++r) for (unsigned int c=0; cchunk_size; ++c) dst(*colnum_ptr * cols->chunk_size + c) @@ -735,15 +765,12 @@ ChunkSparseMatrix::matrix_norm_square (const Vector &v) cons somenumber result = 0; //////////////// - // like matrix_scalar_product, except that - // the two vectors are now the same + // like matrix_scalar_product, except that the two vectors are now the same const unsigned int n_chunk_rows = cols->sparsity_pattern.n_rows(); - // loop over all chunks. note that we need - // to treat the last chunk row and column - // differently if they have padding - // elements + // loop over all chunks. note that we need to treat the last chunk row and + // column differently if they have padding elements const bool rows_have_padding = (m() % cols->chunk_size != 0), cols_have_padding = (n() % cols->chunk_size != 0); @@ -774,8 +801,7 @@ ChunkSparseMatrix::matrix_norm_square (const Vector &v) cons v_ptr, v.begin() + *colnum_ptr * cols->chunk_size); else - // we're at a chunk column that - // has padding + // we're at a chunk column that has padding for (unsigned int r=0; rchunk_size; ++r) for (unsigned int c=0; cchunk_size; ++c) result @@ -792,8 +818,7 @@ ChunkSparseMatrix::matrix_norm_square (const Vector &v) cons v_ptr += cols->chunk_size; } - // now deal with last chunk row if - // necessary + // now deal with last chunk row if necessary if (rows_have_padding) { const unsigned int chunk_row = n_chunk_rows - 1; @@ -807,8 +832,7 @@ ChunkSparseMatrix::matrix_norm_square (const Vector &v) cons || (*colnum_ptr != cols->sparsity_pattern.n_cols()-1)) { - // we're at a chunk row but not - // column that has padding + // we're at a chunk row but not column that has padding for (unsigned int r=0; rchunk_size; ++r) for (unsigned int c=0; cchunk_size; ++c) result @@ -818,8 +842,7 @@ ChunkSparseMatrix::matrix_norm_square (const Vector &v) cons v(*colnum_ptr * cols->chunk_size + c)); } else - // we're at a chunk row and - // column that has padding + // we're at a chunk row and column that has padding for (unsigned int r=0; rchunk_size; ++r) for (unsigned int c=0; cchunk_size; ++c) result @@ -849,16 +872,13 @@ ChunkSparseMatrix::matrix_scalar_product (const Vector &u, Assert(m() == u.size(), ExcDimensionMismatch(m(),u.size())); Assert(n() == v.size(), ExcDimensionMismatch(n(),v.size())); - // the following works like the vmult_add - // function + // the following works like the vmult_add function somenumber result = 0; const unsigned int n_chunk_rows = cols->sparsity_pattern.n_rows(); - // loop over all chunks. note that we need - // to treat the last chunk row and column - // differently if they have padding - // elements + // loop over all chunks. note that we need to treat the last chunk row and + // column differently if they have padding elements const bool rows_have_padding = (m() % cols->chunk_size != 0), cols_have_padding = (n() % cols->chunk_size != 0); @@ -889,8 +909,7 @@ ChunkSparseMatrix::matrix_scalar_product (const Vector &u, u_ptr, v.begin() + *colnum_ptr * cols->chunk_size); else - // we're at a chunk column that - // has padding + // we're at a chunk column that has padding for (unsigned int r=0; rchunk_size; ++r) for (unsigned int c=0; cchunk_size; ++c) result @@ -907,8 +926,7 @@ ChunkSparseMatrix::matrix_scalar_product (const Vector &u, u_ptr += cols->chunk_size; } - // now deal with last chunk row if - // necessary + // now deal with last chunk row if necessary if (rows_have_padding) { const unsigned int chunk_row = n_chunk_rows - 1; @@ -922,8 +940,7 @@ ChunkSparseMatrix::matrix_scalar_product (const Vector &u, || (*colnum_ptr != cols->sparsity_pattern.n_cols()-1)) { - // we're at a chunk row but not - // column that has padding + // we're at a chunk row but not column that has padding for (unsigned int r=0; rchunk_size; ++r) for (unsigned int c=0; cchunk_size; ++c) result @@ -933,8 +950,7 @@ ChunkSparseMatrix::matrix_scalar_product (const Vector &u, v(*colnum_ptr * cols->chunk_size + c)); } else - // we're at a chunk row and - // column that has padding + // we're at a chunk row and column that has padding for (unsigned int r=0; rchunk_size; ++r) for (unsigned int c=0; cchunk_size; ++c) result @@ -962,11 +978,9 @@ ChunkSparseMatrix::l1_norm () const const unsigned int n_chunk_rows = cols->sparsity_pattern.n_rows(); - // loop over all rows and columns; it is - // safe to also loop over the padding - // elements (they are zero) if we make sure - // that the vector into which we sum column - // sums is large enough + // loop over all rows and columns; it is safe to also loop over the padding + // elements (they are zero) if we make sure that the vector into which we + // sum column sums is large enough Vector column_sums(cols->sparsity_pattern.n_cols() * cols->chunk_size); @@ -994,19 +1008,15 @@ ChunkSparseMatrix::linfty_norm () const Assert (cols != 0, ExcNotInitialized()); Assert (val != 0, ExcNotInitialized()); - // this function works like l1_norm(). it - // can be made more efficient (without - // allocating a temporary vector) as is - // done in the SparseMatrix class but since - // it is rarely called in time critical - // places it is probably not worth it + // this function works like l1_norm(). it can be made more efficient + // (without allocating a temporary vector) as is done in the SparseMatrix + // class but since it is rarely called in time critical places it is + // probably not worth it const unsigned int n_chunk_rows = cols->sparsity_pattern.n_rows(); - // loop over all rows and columns; it is - // safe to also loop over the padding - // elements (they are zero) if we make sure - // that the vector into which we sum column - // sums is large enough + // loop over all rows and columns; it is safe to also loop over the padding + // elements (they are zero) if we make sure that the vector into which we + // sum column sums is large enough Vector row_sums(cols->sparsity_pattern.n_rows() * cols->chunk_size); @@ -1030,12 +1040,10 @@ template typename ChunkSparseMatrix::real_type ChunkSparseMatrix::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::abs_square(*ptr); @@ -1060,25 +1068,19 @@ ChunkSparseMatrix::residual (Vector &dst, Assert (&u != &dst, ExcSourceEqualsDestination()); - // set dst=b, then subtract the result of - // A*u from it. since the purpose of the - // current class is to promote streaming of - // data rather than more random access - // patterns, breaking things up into two - // loops may be reasonable + // set dst=b, then subtract the result of A*u from it. since the purpose of + // the current class is to promote streaming of data rather than more random + // access patterns, breaking things up into two loops may be reasonable dst = b; ///////// - // the rest of this function is like - // vmult_add, except that we subtract + // the rest of this function is like vmult_add, except that we subtract // rather than add A*u ///////// const unsigned int n_chunk_rows = cols->sparsity_pattern.n_rows(); - // loop over all chunks. note that we need - // to treat the last chunk row and column - // differently if they have padding - // elements + // loop over all chunks. note that we need to treat the last chunk row and + // column differently if they have padding elements const bool rows_have_padding = (m() % cols->chunk_size != 0), cols_have_padding = (n() % cols->chunk_size != 0); @@ -1107,8 +1109,7 @@ ChunkSparseMatrix::residual (Vector &dst, u.begin() + *colnum_ptr * cols->chunk_size, dst_ptr); else - // we're at a chunk column that - // has padding + // we're at a chunk column that has padding for (unsigned int r=0; rchunk_size; ++r) for (unsigned int c=0; cchunk_size; ++c) dst(chunk_row * cols->chunk_size + r) @@ -1123,8 +1124,7 @@ ChunkSparseMatrix::residual (Vector &dst, dst_ptr += cols->chunk_size; } - // now deal with last chunk row if - // necessary + // now deal with last chunk row if necessary if (rows_have_padding) { const unsigned int chunk_row = n_chunk_rows - 1; @@ -1138,8 +1138,7 @@ ChunkSparseMatrix::residual (Vector &dst, || (*colnum_ptr != cols->sparsity_pattern.n_cols()-1)) { - // we're at a chunk row but not - // column that has padding + // we're at a chunk row but not column that has padding for (unsigned int r=0; rchunk_size; ++r) for (unsigned int c=0; cchunk_size; ++c) dst(chunk_row * cols->chunk_size + r) @@ -1147,8 +1146,7 @@ ChunkSparseMatrix::residual (Vector &dst, u(*colnum_ptr * cols->chunk_size + c)); } else - // we're at a chunk row and - // column that has padding + // we're at a chunk row and column that has padding for (unsigned int r=0; rchunk_size; ++r) for (unsigned int c=0; cchunk_size; ++c) dst(chunk_row * cols->chunk_size + r) @@ -1195,10 +1193,8 @@ ChunkSparseMatrix::precondition_SSOR (Vector &dst, const Vector &src, const number /*om*/) const { - // to understand how this function works - // you may want to take a look at the CVS - // archives to see the original version - // which is much clearer... + // to understand how this function works you may want to take a look at the + // CVS archives to see the original version which is much clearer... Assert (cols != 0, ExcNotInitialized()); Assert (val != 0, ExcNotInitialized()); Assert (m() == n(), ExcMessage("This operation is only valid on square matrices.")); @@ -1387,16 +1383,6 @@ ChunkSparseMatrix::SSOR (Vector &dst, -template -const ChunkSparsityPattern & -ChunkSparseMatrix::get_sparsity_pattern () const -{ - Assert (cols != 0, ExcNotInitialized()); - return *cols; -} - - - template void ChunkSparseMatrix::print (std::ostream &out) const { @@ -1474,9 +1460,8 @@ void ChunkSparseMatrix::print_pattern (std::ostream &out, const unsigned int chunk_size = cols->get_chunk_size(); - // loop over all chunk rows and columns, - // and each time we find something repeat - // it chunk_size times in both directions + // loop over all chunk rows and columns, and each time we find something + // repeat it chunk_size times in both directions for (unsigned int i=0; isparsity_pattern.n_rows(); ++i) { for (unsigned int d=0; d::block_write (std::ostream &out) const { AssertThrow (out, ExcIO()); - // first the simple objects, - // bracketed in [...] + // first the simple objects, bracketed in [...] out << '[' << max_len << "]["; // then write out real data out.write (reinterpret_cast(&val[0]), diff --git a/deal.II/include/deal.II/lac/chunk_sparsity_pattern.h b/deal.II/include/deal.II/lac/chunk_sparsity_pattern.h index 948c8d017e..655c0ba406 100644 --- a/deal.II/include/deal.II/lac/chunk_sparsity_pattern.h +++ b/deal.II/include/deal.II/lac/chunk_sparsity_pattern.h @@ -34,6 +34,187 @@ template class ChunkSparseMatrix; */ + +/** + * Iterators on sparsity patterns + */ +namespace ChunkSparsityPatternIterators +{ + // forward declaration + class Iterator; + + /** + * Accessor class for iterators into sparsity patterns. This class is + * also the base class for both const and non-const accessor classes + * into sparse matrices. + * + * Note that this class only allows read access to elements, providing + * their row and column number. It does not allow modifying the + * sparsity pattern itself. + * + * @author Martin Kronbichler + * @date 2013 + */ + class Accessor + { + public: + /** + * Constructor. + */ + Accessor (const ChunkSparsityPattern *matrix, + const unsigned int row); + + /** + * Constructor. Construct the end accessor for the given sparsity pattern. + */ + Accessor (const ChunkSparsityPattern *matrix); + + /** + * Row number of the element represented by this object. This function can + * only be called for entries for which is_valid_entry() is true. + */ + unsigned int row () const; + + /** + * Index in row of the element represented by this object. This function + * can only be called for entries for which is_valid_entry() is true. + */ + unsigned int index () const; + + /** + * Returns the global index from the reduced sparsity pattern. + */ + std::size_t reduced_index() const; + + /** + * Column number of the element represented by this object. This function + * can only be called for entries for which is_valid_entry() is true. + */ + unsigned int column () const; + + /** + * Return whether the sparsity pattern entry pointed to by this iterator + * is valid or not. Note that after compressing the sparsity pattern, all + * entries are valid. However, before compression, the sparsity pattern + * allocated some memory to be used while still adding new nonzero + * entries; if you create iterators in this phase of the sparsity + * pattern's lifetime, you will iterate over elements that are not + * valid. If this is so, then this function will return false. + */ + bool is_valid_entry () const; + + + /** + * Comparison. True, if both iterators point to the same matrix position. + */ + bool operator == (const Accessor &) const; + + + /** + * Comparison operator. Result is true if either the first row number is + * smaller or if the row numbers are equal and the first index is smaller. + * + * This function is only valid if both iterators point into the same + * sparsity pattern. + */ + bool operator < (const Accessor &) const; + + protected: + /** + * The sparsity pattern we operate on accessed. + */ + const ChunkSparsityPattern *sparsity_pattern; + + /** + * The accessor of the (reduced) sparsity pattern. + */ + SparsityPatternIterators::Accessor reduced_accessor; + + /** + * Current row number. + */ + unsigned int chunk_row; + + /** + * Current index in row. + */ + unsigned int chunk_col; + + /** + * Move the accessor to the next nonzero entry in the matrix. + */ + void advance (); + + /** + * Grant access to iterator class. + */ + friend class Iterator; + }; + + + + /** + * STL conforming iterator walking over the elements of a sparsity pattern. + */ + class Iterator + { + public: + /** + * Constructor. Create an iterator into the sparsity pattern @p sp for the + * given row and the index within it. + */ + Iterator (const ChunkSparsityPattern *sp, + const unsigned int row); + + /** + * Prefix increment. + */ + Iterator &operator++ (); + + /** + * Postfix increment. + */ + Iterator operator++ (int); + + /** + * Dereferencing operator. + */ + const Accessor &operator* () const; + + /** + * Dereferencing operator. + */ + const Accessor *operator-> () const; + + /** + * Comparison. True, if both iterators point to the same matrix position. + */ + bool operator == (const Iterator &) const; + + /** + * Inverse of ==. + */ + bool operator != (const Iterator &) const; + + /** + * Comparison operator. Result is true if either the first row number is + * smaller or if the row numbers are equal and the first index is smaller. + * + * This function is only valid if both iterators point into the same + * matrix. + */ + bool operator < (const Iterator &) const; + + private: + /** + * Store an object of the accessor class. + */ + Accessor accessor; + }; +} + + + /** * Structure representing the sparsity pattern of a sparse matrix. * @@ -46,86 +227,67 @@ template class ChunkSparseMatrix; class ChunkSparsityPattern : public Subscriptor { public: + /** + * Typedef an iterator class that allows to walk over all nonzero elements + * of a sparsity pattern. + */ + typedef ChunkSparsityPatternIterators::Iterator const_iterator; + + /** + * Typedef an iterator class that allows to walk over all nonzero elements + * of a sparsity pattern. + * + * Since the iterator does not allow to modify the sparsity pattern, this + * type is the same as that for @p const_iterator. + */ + typedef ChunkSparsityPatternIterators::Iterator iterator; /** - * Define a value which is used - * to indicate that a certain - * value in the colnums array - * is unused, i.e. does not - * represent a certain column - * number index. + * Define a value which is used to indicate that a certain value in the + * colnums array is unused, i.e. does not represent a certain column number + * index. * - * Indices with this invalid - * value are used to insert new - * entries to the sparsity - * pattern using the add() member - * function, and are removed when + * Indices with this invalid value are used to insert new entries to the + * sparsity pattern using the add() member function, and are removed when * calling compress(). * - * You should not assume that the - * variable declared here has a - * certain value. The - * initialization is given here - * only to enable the compiler to - * perform some optimizations, - * but the actual value of the - * variable may change over time. + * You should not assume that the variable declared here has a certain + * value. The initialization is given here only to enable the compiler to + * perform some optimizations, but the actual value of the variable may + * change over time. */ static const unsigned int invalid_entry = SparsityPattern::invalid_entry; /** - * Initialize the matrix empty, - * that is with no memory - * allocated. This is useful if - * you want such objects as - * member variables in other - * classes. You can make the - * structure usable by calling - * the reinit() function. + * Initialize the matrix empty, that is with no memory allocated. This is + * useful if you want such objects as member variables in other classes. You + * can make the structure usable by calling the reinit() function. */ ChunkSparsityPattern (); /** - * Copy constructor. This - * constructor is only allowed to - * be called if the matrix - * structure to be copied is - * empty. This is so in order to - * prevent involuntary copies of - * objects for temporaries, which - * can use large amounts of - * computing time. However, copy - * constructors are needed if yo - * want to use the STL data types - * on classes like this, e.g. to - * write such statements like - * v.push_back - * (ChunkSparsityPattern());, - * with v 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 v.push_back (ChunkSparsityPattern());, with + * v a vector of ChunkSparsityPattern objects. + * + * Usually, it is sufficient to use the explicit keyword to disallow + * unwanted temporaries, but for the STL vectors, this does not work. Since + * copying a structure like this is not useful anyway because multiple + * matrices can use the same sparsity structure, copies are only allowed for + * empty objects, as described above. */ ChunkSparsityPattern (const ChunkSparsityPattern &); /** - * Initialize a rectangular - * matrix. + * Initialize a rectangular matrix. * * @arg m number of rows * @arg n number of columns - * @arg max_per_row maximum - * number of nonzero entries per row + * @arg max_per_row maximum number of nonzero entries per row */ ChunkSparsityPattern (const unsigned int m, const unsigned int n, @@ -143,16 +305,12 @@ public: const bool optimize_diagonal) DEAL_II_DEPRECATED; /** - * Initialize a rectangular - * matrix. + * Initialize a rectangular matrix. * * @arg m number of rows * @arg n number of columns - * - * @arg row_lengths possible - * number of nonzero entries for - * each row. This vector must - * have one entry for each row. + * @arg row_lengths possible number of nonzero entries for each row. This + * vector must have one entry for each row. */ ChunkSparsityPattern (const unsigned int m, const unsigned int n, @@ -170,17 +328,12 @@ public: const bool optimize_diagonal) DEAL_II_DEPRECATED; /** - * Initialize a quadratic matrix - * of dimension n with - * at most max_per_row - * nonzero entries per row. + * Initialize a quadratic matrix of dimension n with at most + * max_per_row nonzero entries per row. * - * This constructor automatically - * enables optimized storage of - * diagonal elements. To avoid - * this, use the constructor - * taking row and column numbers - * separately. + * This constructor automatically enables optimized storage of diagonal + * elements. To avoid this, use the constructor taking row and column + * numbers separately. */ ChunkSparsityPattern (const unsigned int n, const unsigned int max_per_row, @@ -190,11 +343,8 @@ public: * Initialize a quadratic matrix. * * @arg m number of rows and columns - * - * @arg row_lengths possible - * number of nonzero entries for - * each row. This vector must - * have one entry for each row. + * @arg row_lengths possible number of nonzero entries for each row. This + * vector must have one entry for each row. */ ChunkSparsityPattern (const unsigned int m, const std::vector &row_lengths, @@ -215,25 +365,19 @@ public: ~ChunkSparsityPattern (); /** - * Copy operator. For this the - * same holds as for the copy - * constructor: it is declared, - * defined and fine to be called, - * but the latter only for empty + * Copy operator. For this the same holds as for the copy constructor: it is + * declared, defined and fine to be called, but the latter only for empty * objects. */ ChunkSparsityPattern &operator = (const ChunkSparsityPattern &); /** - * Reallocate memory and set up data - * structures for a new matrix with - * m rows and n columns, - * with at most max_per_row + * Reallocate memory and set up data structures for a new matrix with m + * rows and n columns, with at most max_per_row * nonzero entries per row. * - * This function simply maps its - * operations to the other - * reinit function. + * This function simply maps its operations to the other reinit + * function. */ void reinit (const unsigned int m, const unsigned int n, @@ -251,21 +395,14 @@ public: const bool optimize_diagonal) DEAL_II_DEPRECATED; /** - * Reallocate memory for a matrix - * of size m x n. The - * number of entries for each row - * is taken from the array - * row_lengths which has to - * give this number of each row - * i=1...m. + * Reallocate memory for a matrix of size m x n. The number of + * entries for each row is taken from the array row_lengths which + * has to give this number of each row i=1...m. * - * If m*n==0 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 m*n==0 all memory is freed, resulting in a total + * reinitialization of the object. If it is nonzero, new memory is only + * allocated if the new size extends the old one. This is done to save time + * and to avoid fragmentation of the heap. * * If the number of rows equals * the number of columns then @@ -290,8 +427,7 @@ public: const bool optimize_diagonal) DEAL_II_DEPRECATED; /** - * Same as above, but with a - * VectorSlice argument instead. + * Same as above, but with a VectorSlice argument instead. */ void reinit (const unsigned int m, const unsigned int n, @@ -309,103 +445,52 @@ public: const bool optimize_diagonal) DEAL_II_DEPRECATED; /** - * This function compresses the sparsity - * structure that this object represents. - * It does so by eliminating unused - * entries and sorting the remaining ones - * to allow faster access by usage of - * binary search algorithms. A special - * sorting scheme is used for the - * diagonal entry of quadratic matrices, - * which is always the first entry of - * each row. + * This function compresses the sparsity structure that this object + * represents. It does so by eliminating unused entries and sorting the + * remaining ones to allow faster access by usage of binary search + * algorithms. A special sorting scheme is used for the diagonal entry of + * quadratic matrices, which is always the first entry of each row. * - * The memory which is no more - * needed is released. + * The memory which is no more needed is released. * - * SparseMatrix objects require the - * ChunkSparsityPattern objects they are - * initialized with to be compressed, to - * reduce memory requirements. + * SparseMatrix objects require the ChunkSparsityPattern objects they are + * initialized with to be compressed, to reduce memory requirements. */ void compress (); /** - * This function can be used as a - * replacement for reinit(), - * subsequent calls to add() and - * a final call to close() if you - * know exactly in advance the - * entries that will form the - * matrix sparsity pattern. - * - * The first two parameters - * determine the size of the - * matrix. For the two last ones, - * note that a sparse matrix can - * be described by a sequence of - * rows, each of which is - * represented by a sequence of - * pairs of column indices and - * values. In the present - * context, the begin() and - * end() parameters designate - * iterators (of forward iterator - * type) into a container, one - * representing one row. The - * distance between begin() - * and end() should therefore - * be equal to - * n_rows(). These iterators - * may be iterators of - * std::vector, - * std::list, 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 operator* or - * operator-> to one of these - * iterators) must be a container - * itself that provides functions - * begin and end - * 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 std::vector, std::list, 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 operator* or operator-> to one of + * these iterators) must be a container itself that provides functions + * begin and end 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 > column_indices (n_rows); * for (unsigned int row=0; rowbegin and - * end (namely - * std::vectors), and the - * inner iterators dereferenced - * yield unsigned integers as - * column indices. Note that we - * could have replaced each of - * the two std::vector - * occurrences by std::list, - * and the inner one by - * std::set 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 begin and end (namely + * std::vectors), and the inner iterators dereferenced yield + * unsigned integers as column indices. Note that we could have replaced + * each of the two std::vector occurrences by std::list, + * and the inner one by std::set as well. + * + * Another example would be as follows, where we initialize a whole matrix, + * not only a sparsity pattern: * @code * std::vector > entries (n_rows); * for (unsigned int row=0; rowstd::vector - * could be replaced by - * std::list, and the inner - * std::map - * could be replaced by - * std::vector >, - * 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 std::vector could + * be replaced by std::list, and the inner std::map could be replaced by std::vector >, or a list or set of such pairs, as they all return + * iterators that point to such pairs. */ template void copy_from (const unsigned int n_rows, @@ -488,12 +554,9 @@ public: const bool optimize_diagonal) DEAL_II_DEPRECATED; /** - * Copy data from an object of type - * CompressedSparsityPattern, - * CompressedSetSparsityPattern or - * CompressedSimpleSparsityPattern. - * Previous content of this object is - * lost, and the sparsity pattern is in + * Copy data from an object of type CompressedSparsityPattern, + * CompressedSetSparsityPattern or CompressedSimpleSparsityPattern. + * Previous content of this object is lost, and the sparsity pattern is in * compressed mode afterwards. */ template @@ -510,14 +573,10 @@ public: const bool optimize_diagonal) DEAL_II_DEPRECATED; /** - * Take a full matrix and use its - * nonzero entries to generate a - * sparse matrix entry pattern - * for this object. + * Take a full matrix and use its nonzero entries to generate a sparse + * matrix entry pattern for this object. * - * Previous content of this - * object is lost, and the - * sparsity pattern is in + * Previous content of this object is lost, and the sparsity pattern is in * compressed mode afterwards. */ template @@ -534,71 +593,79 @@ public: const bool optimize_diagonal) DEAL_II_DEPRECATED; /** - * Return whether the object is empty. It - * is empty if no memory is allocated, - * which is the same as that both - * dimensions are zero. + * Set the sparsity pattern of the chunk sparsity pattern to be given by + * chunk_size*chunksize blocks of the sparsity pattern for chunks + * specified. Note that the final number of rows m of the sparsity + * pattern will be approximately sparsity_pattern_for_chunks.n_rows() * + * chunk_size (modulo padding elements in the last chunk) and similarly + * for the number of columns n. + * + * 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 make_sparsity_pattern 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 + void create_from (const unsigned int m, + const unsigned int n, + const Sparsity &sparsity_pattern_for_chunks, + const unsigned int chunk_size, + const bool optimize_diagonal = true); + + /** + * Return whether the object is empty. It is empty if no memory is + * allocated, which is the same as that both dimensions are zero. */ bool empty () const; /** - * Return the chunk size given as - * argument when constructing this - * object. + * Return the chunk size given as argument when constructing this object. */ unsigned int get_chunk_size () const; /** - * Return the maximum number of entries per - * row. Before compression, this equals the - * number given to the constructor, while - * after compression, it equals the maximum - * number of entries actually allocated by - * the user. + * Return the maximum number of entries per row. Before compression, this + * equals the number given to the constructor, while after compression, it + * equals the maximum number of entries actually allocated by the user. */ unsigned int max_entries_per_row () const; /** - * Add a nonzero entry to the matrix. - * This function may only be called - * for non-compressed sparsity patterns. + * Add a nonzero entry to the matrix. This function may only be called for + * non-compressed sparsity patterns. * - * If the entry already exists, nothing - * bad happens. + * If the entry already exists, nothing bad happens. */ void add (const unsigned int i, const unsigned int j); /** - * Make the sparsity pattern - * symmetric by adding the - * sparsity pattern of the + * Make the sparsity pattern symmetric by adding the sparsity pattern of the * transpose object. * - * This function throws an - * exception if the sparsity - * pattern does not represent a - * quadratic matrix. + * This function throws an exception if the sparsity pattern does not + * represent a quadratic matrix. */ void symmetrize (); /** - * Return number of rows of this - * matrix, which equals the dimension - * of the image space. + * Return number of rows of this matrix, which equals the dimension of the + * image space. */ - inline unsigned int n_rows () const; + unsigned int n_rows () const; /** - * Return number of columns of this - * matrix, which equals the dimension - * of the range space. + * Return number of columns of this matrix, which equals the dimension of + * the range space. */ - inline unsigned int n_cols () const; + unsigned int n_cols () const; /** - * Check if a value at a certain - * position may be non-zero. + * Check if a value at a certain position may be non-zero. */ bool exists (const unsigned int i, const unsigned int j) const; @@ -609,34 +676,25 @@ public: unsigned int row_length (const unsigned int row) const; /** - * Compute the bandwidth of the matrix - * represented by this structure. The - * bandwidth is the maximum of $|i-j|$ - * for which the index pair $(i,j)$ - * represents a nonzero entry of the - * matrix. Consequently, the maximum - * bandwidth a $n\times m$ matrix can - * have is $\max\{n-1,m-1\}$. + * Compute the bandwidth of the matrix represented by this structure. The + * bandwidth is the maximum of $|i-j|$ for which the index pair $(i,j)$ + * represents a nonzero entry of the matrix. Consequently, the maximum + * bandwidth a $n\times m$ matrix can have is $\max\{n-1,m-1\}$. */ unsigned int bandwidth () const; /** - * Return the number of nonzero elements of - * this matrix. Actually, it returns the - * number of entries in the sparsity - * pattern; if any of the entries should - * happen to be zero, it is counted - * anyway. + * Return the number of nonzero elements of this matrix. Actually, it + * returns the number of entries in the sparsity pattern; if any of the + * entries should happen to be zero, it is counted anyway. * - * This function may only be called if the - * matrix struct is compressed. It does not - * make too much sense otherwise anyway. + * This function may only be called if the matrix struct is compressed. It + * does not make too much sense otherwise anyway. */ unsigned int n_nonzero_elements () const; /** - * Return whether the structure is - * compressed or not. + * Return whether the structure is compressed or not. */ bool is_compressed () const; @@ -663,94 +721,89 @@ public: bool stores_only_added_elements () const; /** - * Write the data of this object - * en bloc to a file. This is - * done in a binary mode, so the - * output is neither readable by - * humans nor (probably) by other - * computers using a different - * operating system of number - * format. - * - * The purpose of this function - * is that you can swap out - * matrices and sparsity pattern - * if you are short of memory, - * want to communicate between - * different programs, or allow - * objects to be persistent - * across different runs of the - * program. + * STL-like iterator with the first entry of the matrix. The resulting + * iterator can be used to walk over all nonzero entries of the sparsity + * pattern. + */ + iterator begin () const; + + /** + * Final iterator. + */ + iterator end () const; + + /** + * STL-like iterator with the first entry of row r. + * + * Note that if the given row is empty, i.e. does not contain any nonzero + * entries, then the iterator returned by this function equals + * end(r). Note also that the iterator may not be dereferencable in + * that case. + */ + iterator begin (const unsigned int r) const; + + /** + * Final iterator of row r. 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 - * [i,j1,j2,j3,...]. i - * is the row number and - * jn are the allocated - * columns in this row. + * Print the sparsity of the matrix. The output consists of one line per row + * of the format [i,j1,j2,j3,...]. i is the row number and + * jn are the allocated columns in this row. */ void print (std::ostream &out) const; /** - * Print the sparsity of the matrix - * in a format that gnuplot understands - * and which can be used to plot the - * sparsity pattern in a graphical - * way. The format consists of pairs - * i j 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 - * (0,0) 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 - * plot command. + * Print the sparsity of the matrix in a format that gnuplot + * understands and which can be used to plot the sparsity pattern in a + * graphical way. The format consists of pairs i j 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 (0,0) + * 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 plot command. */ void print_gnuplot (std::ostream &out) const; /** - * Determine an estimate for the - * memory consumption (in bytes) - * of this object. See - * MemoryConsumption. + * Determine an estimate for the memory consumption (in bytes) of this + * object. See MemoryConsumption. */ std::size_t memory_consumption () const; @@ -822,14 +875,12 @@ public: //@} private: /** - * Number of rows that this sparsity - * structure shall represent. + * Number of rows that this sparsity structure shall represent. */ unsigned int rows; /** - * Number of columns that this sparsity - * structure shall represent. + * Number of columns that this sparsity structure shall represent. */ unsigned int cols; @@ -839,18 +890,20 @@ private: unsigned int chunk_size; /** - * The reduced sparsity pattern. We store - * only which chunks exist, with each - * chunk a block in the matrix of size - * chunk_size by chunk_size. + * The reduced sparsity pattern. We store only which chunks exist, with each + * chunk a block in the matrix of size chunk_size by chunk_size. */ SparsityPattern sparsity_pattern; /** - * Make all the chunk sparse matrix kinds - * friends. + * Make all the chunk sparse matrix kinds friends. */ template friend class ChunkSparseMatrix; + + /** + * Make the accessor class a friend. + */ + friend class ChunkSparsityPatternIterators::Accessor; }; @@ -859,6 +912,330 @@ private: #ifndef DOXYGEN +namespace ChunkSparsityPatternIterators +{ + inline + Accessor:: + Accessor (const ChunkSparsityPattern *sparsity_pattern, + const unsigned int row) + : + sparsity_pattern(sparsity_pattern), + reduced_accessor(row==sparsity_pattern->n_rows() ? + *sparsity_pattern->sparsity_pattern.end() : + *sparsity_pattern->sparsity_pattern. + begin(row/sparsity_pattern->get_chunk_size())), + chunk_row (row==sparsity_pattern->n_rows() ? 0 : + row%sparsity_pattern->get_chunk_size()), + chunk_col (0) + {} + + + + inline + Accessor:: + Accessor (const ChunkSparsityPattern *sparsity_pattern) + : + sparsity_pattern(sparsity_pattern), + reduced_accessor(*sparsity_pattern->sparsity_pattern.end()), + chunk_row (0), + chunk_col (0) + {} + + + + inline + bool + Accessor::is_valid_entry () const + { + return reduced_accessor.is_valid_entry() + && + sparsity_pattern->get_chunk_size()*reduced_accessor.row()+chunk_row < + sparsity_pattern->n_rows() + && + sparsity_pattern->get_chunk_size()*reduced_accessor.column()+chunk_col < + sparsity_pattern->n_cols(); + } + + + + inline + unsigned int + Accessor::row() const + { + Assert (is_valid_entry() == true, ExcInvalidIterator()); + + return sparsity_pattern->get_chunk_size()*reduced_accessor.row()+chunk_row; + } + + + + inline + unsigned int + Accessor::column() const + { + Assert (is_valid_entry() == true, ExcInvalidIterator()); + + return sparsity_pattern->get_chunk_size()*reduced_accessor.column() + + chunk_col; + } + + + + inline + unsigned int + Accessor::index() const + { + Assert (is_valid_entry() == true, ExcInvalidIterator()); + + return sparsity_pattern->get_chunk_size()*reduced_accessor.index() + + chunk_col; + } + + + + inline + std::size_t + Accessor::reduced_index() const + { + Assert (is_valid_entry() == true, ExcInvalidIterator()); + + return reduced_accessor.index_within_sparsity; + } + + + + + inline + bool + Accessor::operator == (const Accessor &other) const + { + // no need to check for equality of sparsity patterns as this is done in + // the reduced case already and every ChunkSparsityPattern has its own + // reduced sparsity pattern + return (reduced_accessor == other.reduced_accessor && + chunk_row == other.chunk_row && + chunk_col == other.chunk_col); + } + + + + inline + bool + Accessor::operator < (const Accessor &other) const + { + Assert (sparsity_pattern == other.sparsity_pattern, + ExcInternalError()); + + // comparison is a bit messy because of the way ChunkSparsityPattern + // stores entry: chunk rows run faster than the indices of the reduced + // sparsity pattern, but the accessors should of course compare less based + // on the actual row, not the reduced one. + if (chunk_row == other.chunk_row) + return (reduced_accessor.index_within_sparsity < + other.reduced_accessor.index_within_sparsity || + (reduced_accessor.index_within_sparsity == + other.reduced_accessor.index_within_sparsity && + chunk_col < other.chunk_col)); + else + { + if (reduced_accessor.index_within_sparsity == + reduced_accessor.sparsity_pattern->n_nonzero_elements()) + return false; + if (other.reduced_accessor.index_within_sparsity == + reduced_accessor.sparsity_pattern->n_nonzero_elements()) + return true; + const unsigned int row = reduced_accessor.row(), + other_row = other.reduced_accessor.row(); + return (row < other_row + || + (row == other_row + && + (chunk_row < other.chunk_row + || + (chunk_row == other.chunk_row + && + (reduced_accessor.index_within_sparsity < + other.reduced_accessor.index_within_sparsity + || + (reduced_accessor.index_within_sparsity == + other.reduced_accessor.index_within_sparsity + && + chunk_col < other.chunk_col)))))); + } + } + + + inline + void + Accessor::advance () + { + const unsigned int chunk_size = sparsity_pattern->get_chunk_size(); + Assert (chunk_row < chunk_size && chunk_col < chunk_size, + ExcIteratorPastEnd()); + Assert (reduced_accessor.row() * chunk_size + chunk_row < + sparsity_pattern->n_rows() + && + reduced_accessor.column() * chunk_size + chunk_col < + sparsity_pattern->n_cols(), + ExcIteratorPastEnd()); + + ++chunk_col; + + // end of chunk + if (chunk_col == chunk_size + || + reduced_accessor.column() * chunk_size + chunk_col == + sparsity_pattern->n_cols()) + { + // end of row + if (reduced_accessor.index() + 1 == + reduced_accessor.sparsity_pattern->row_length(reduced_accessor.row())) + { + ++chunk_row; + chunk_col = 0; + const unsigned int old_reduced_row = reduced_accessor.row(); + + // end of matrix + if (old_reduced_row * chunk_size + chunk_row == + sparsity_pattern->n_rows()) + { + chunk_row = 0; + reduced_accessor = + SparsityPatternIterators::Accessor(&sparsity_pattern->sparsity_pattern, + sparsity_pattern->sparsity_pattern.n_nonzero_elements()); + } + // end of chunk rows + else if (chunk_row == chunk_size) + { + reduced_accessor = + *sparsity_pattern->sparsity_pattern.begin(old_reduced_row+1); + chunk_row = 0; + } + else + reduced_accessor = + *sparsity_pattern->sparsity_pattern.begin(old_reduced_row); + } + else + { + reduced_accessor.advance(); + chunk_col = 0; + } + } + } + + + + inline + Iterator::Iterator (const ChunkSparsityPattern *sparsity_pattern, + const unsigned int row) + : + accessor(sparsity_pattern, row) + {} + + + + inline + Iterator & + Iterator::operator++ () + { + accessor.advance (); + return *this; + } + + + + inline + Iterator + Iterator::operator++ (int) + { + const Iterator iter = *this; + accessor.advance (); + return iter; + } + + + + inline + const Accessor & + Iterator::operator* () const + { + return accessor; + } + + + + inline + const Accessor * + Iterator::operator-> () const + { + return &accessor; + } + + + inline + bool + Iterator::operator == (const Iterator &other) const + { + return (accessor == other.accessor); + } + + + + inline + bool + Iterator::operator != (const Iterator &other) const + { + return ! (*this == other); + } + + + inline + bool + Iterator::operator < (const Iterator &other) const + { + return accessor < other.accessor; + } + +} + + + +inline +ChunkSparsityPattern::iterator +ChunkSparsityPattern::begin () const +{ + return iterator(this, 0); +} + + +inline +ChunkSparsityPattern::iterator +ChunkSparsityPattern::end () const +{ + return iterator(this, n_rows()); +} + + + +inline +ChunkSparsityPattern::iterator +ChunkSparsityPattern::begin (const unsigned int r) const +{ + Assert (r +inline void ChunkSparsityPattern::copy_from (const unsigned int n_rows, const unsigned int n_cols, @@ -928,16 +1307,10 @@ ChunkSparsityPattern::copy_from (const unsigned int n_rows, Assert (static_cast(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 row_lengths; row_lengths.reserve(n_rows); @@ -947,8 +1320,7 @@ ChunkSparsityPattern::copy_from (const unsigned int n_rows, (is_square ? 1 : 0)); reinit (n_rows, n_cols, row_lengths, chunk_size); - // now enter all the elements into - // the matrix + // now enter all the elements into the matrix unsigned int row = 0; typedef typename std::iterator_traits::value_type::const_iterator inner_iterator; for (ForwardIterator i=begin; i!=end; ++i, ++row) @@ -957,16 +1329,14 @@ ChunkSparsityPattern::copy_from (const unsigned int n_rows, for (inner_iterator j=i->begin(); j!=end_of_row; ++j) { const unsigned int col - = internal::SparsityPatternTools::get_column_index_from_iterator(*j); + = internal::SparsityPatternTools::get_column_index_from_iterator(*j); Assert (col < n_cols, ExcInvalidIndex(col,n_cols)); add (row, col); } } - // finally compress - // everything. this also sorts the - // entries within each row + // finally compress everything. this also sorts the entries within each row compress (); } diff --git a/deal.II/source/lac/chunk_sparsity_pattern.cc b/deal.II/source/lac/chunk_sparsity_pattern.cc index 788ed489b9..5bcb7c0470 100644 --- a/deal.II/source/lac/chunk_sparsity_pattern.cc +++ b/deal.II/source/lac/chunk_sparsity_pattern.cc @@ -138,8 +138,7 @@ ChunkSparsityPattern::operator = (const ChunkSparsityPattern &s) Assert (s.rows == 0, ExcInvalidConstructorCall()); Assert (s.cols == 0, ExcInvalidConstructorCall()); - // perform the checks in the underlying - // object as well + // perform the checks in the underlying object as well sparsity_pattern = s.sparsity_pattern; return *this; @@ -166,8 +165,7 @@ ChunkSparsityPattern::reinit (const unsigned int m, { Assert (chunk_size > 0, ExcInvalidNumber (chunk_size)); - // simply map this function to the - // other @p{reinit} function + // simply map this function to the other @p{reinit} function const std::vector row_lengths (m, max_per_row); reinit (m, n, row_lengths, chunk_size); } @@ -201,31 +199,30 @@ ChunkSparsityPattern::reinit ( this->chunk_size = chunk_size; - // pass down to the necessary information - // to the underlying object. we need to - // calculate how many chunks we need: we - // need to round up (m/chunk_size) and - // (n/chunk_size). rounding up in integer - // arithmetic equals + // pass down to the necessary information to the underlying object. we need + // to calculate how many chunks we need: we need to round up (m/chunk_size) + // and (n/chunk_size). rounding up in integer arithmetic equals // ((m+chunk_size-1)/chunk_size): const unsigned int m_chunks = (m+chunk_size-1) / chunk_size, n_chunks = (n+chunk_size-1) / chunk_size; - // compute the maximum number of chunks in - // each row. the passed array denotes the - // number of entries in each row of the big - // matrix -- in the worst case, these are - // all in independent chunks, so we have to - // calculate it as follows (as an example: - // let chunk_size==2, - // row_lengths={2,2,...}, and entries in - // row zero at columns {0,2} and for row - // one at {4,6} --> we'll need 4 chunks for - // the first chunk row!) : + // compute the maximum number of chunks in each row. the passed array + // denotes the number of entries in each row of the big matrix -- in the + // worst case, these are all in independent chunks, so we have to calculate + // it as follows (as an example: let chunk_size==2, row_lengths={2,2,...}, + // and entries in row zero at columns {0,2} and for row one at {4,6} --> + // we'll need 4 chunks for the first chunk row!) : std::vector chunk_row_lengths (m_chunks, 0); for (unsigned int i=0; i 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 entries_per_row (csp.n_rows(), 0); for (unsigned int row = 0; row + void copy_sparsity (const SP &src, + SparsityPattern &dst) + { + dst.copy_from(src); + } + + void copy_sparsity (const SparsityPattern &src, + SparsityPattern &dst) + { + dst = src; + } + } +} + + + +template +void +ChunkSparsityPattern::create_from +(const unsigned int m, + const unsigned int n, + const Sparsity &sparsity_pattern_for_chunks, + const unsigned int chunk_size_in, + const bool) +{ + Assert (m > (sparsity_pattern_for_chunks.n_rows()-1) * chunk_size_in && + m <= sparsity_pattern_for_chunks.n_rows() * chunk_size_in, + ExcMessage("Number of rows m is not compatible with chunk size " + "and number of rows in sparsity pattern for the chunks.")); + Assert (n > (sparsity_pattern_for_chunks.n_cols()-1) * chunk_size_in && + n <= sparsity_pattern_for_chunks.n_cols() * chunk_size_in, + ExcMessage("Number of columns m is not compatible with chunk size " + "and number of columns in sparsity pattern for the chunks.")); + + internal::copy_sparsity(sparsity_pattern_for_chunks, sparsity_pattern); + chunk_size = chunk_size_in; + rows = m; + cols = n; +} + + + bool ChunkSparsityPattern::empty () const { @@ -422,27 +465,41 @@ ChunkSparsityPattern::exists (const unsigned int i, -unsigned int -ChunkSparsityPattern::row_length (const unsigned int i) const +void +ChunkSparsityPattern::symmetrize () { - Assert (icolumn() != sparsity_pattern.n_cols() - 1) + n += chunk_size; + else + n += (n_cols() % chunk_size); + return n; + } } @@ -457,18 +514,15 @@ ChunkSparsityPattern::n_nonzero_elements () const chunk_size * chunk_size); else - // some of the chunks reach beyond the - // extent of this matrix. this requires a - // somewhat more complicated - // computations, in particular if the + // some of the chunks reach beyond the extent of this matrix. this + // requires a somewhat more complicated computations, in particular if the // columns don't align { if ((n_rows() % chunk_size != 0) && (n_cols() % chunk_size == 0)) { - // columns align with chunks, but - // not rows + // columns align with chunks, but not rows unsigned int n = sparsity_pattern.n_nonzero_elements() * chunk_size * chunk_size; @@ -480,10 +534,9 @@ ChunkSparsityPattern::n_nonzero_elements () const else { - // if columns don't align, then - // just iterate over all chunks and - // see what this leads to. follow the advice in the documentation of - // the sparsity pattern iterators to do the loop over individual rows, + // if columns don't align, then just iterate over all chunks and see + // what this leads to. follow the advice in the documentation of the + // sparsity pattern iterators to do the loop over individual rows, // rather than all elements unsigned int n = 0; @@ -498,17 +551,13 @@ ChunkSparsityPattern::n_nonzero_elements () const else if ((row == sparsity_pattern.n_rows() - 1) && (p->column() != sparsity_pattern.n_cols() - 1)) - // last chunk row, but not - // last chunk column. only a - // smaller number (n_rows % - // chunk_size) of rows - // actually exist + // last chunk row, but not last chunk column. only a smaller + // number (n_rows % chunk_size) of rows actually exist n += (n_rows() % chunk_size) * chunk_size; else if ((row != sparsity_pattern.n_rows() - 1) && (p->column() == sparsity_pattern.n_cols() - 1)) - // last chunk column, but - // not row + // last chunk column, but not row n += (n_cols() % chunk_size) * chunk_size; else // bottom right chunk @@ -561,8 +610,7 @@ ChunkSparsityPattern::print_gnuplot (std::ostream &out) const AssertThrow (out, ExcIO()); - // for each entry in the underlying - // sparsity pattern, repeat everything + // for each entry in the underlying sparsity pattern, repeat everything // chunk_size x chunk_size times for (unsigned int i=0; i(i*chunk_size+e) << std::endl; @@ -593,13 +638,10 @@ ChunkSparsityPattern::print_gnuplot (std::ostream &out) const unsigned int ChunkSparsityPattern::bandwidth () const { - // calculate the bandwidth from that of the - // underlying sparsity pattern. note that - // even if the bandwidth of that is zero, - // then the bandwidth of the chunky pattern - // is chunk_size-1, if it is 1 then the - // chunky pattern has - // chunk_size+(chunk_size-1), etc + // calculate the bandwidth from that of the underlying sparsity + // pattern. note that even if the bandwidth of that is zero, then the + // bandwidth of the chunky pattern is chunk_size-1, if it is 1 then the + // chunky pattern has chunk_size+(chunk_size-1), etc // // we'll cut it off at max(n(),m()) return std::min (sparsity_pattern.bandwidth()*chunk_size @@ -625,8 +667,7 @@ ChunkSparsityPattern::block_write (std::ostream &out) const { AssertThrow (out, ExcIO()); - // first the simple objects, - // bracketed in [...] + // first the simple objects, bracketed in [...] out << '[' << rows << ' ' << cols << ' ' @@ -660,8 +701,7 @@ ChunkSparsityPattern::block_read (std::istream &in) in >> c; AssertThrow (c == '[', ExcIO()); - // then read the underlying sparsity - // pattern + // then read the underlying sparsity pattern sparsity_pattern.block_read (in); in >> c; @@ -681,6 +721,10 @@ ChunkSparsityPattern::memory_consumption () const // explicit instantiations template +void ChunkSparsityPattern::copy_from (const SparsityPattern &, + const unsigned int, + const bool); +template void ChunkSparsityPattern::copy_from (const CompressedSparsityPattern &, const unsigned int, const bool); @@ -693,6 +737,34 @@ void ChunkSparsityPattern::copy_from (const Com const unsigned int, const bool); template +void ChunkSparsityPattern::create_from +(const unsigned int, + const unsigned int, + const SparsityPattern &, + const unsigned int, + const bool); +template +void ChunkSparsityPattern::create_from +(const unsigned int, + const unsigned int, + const CompressedSparsityPattern &, + const unsigned int, + const bool); +template +void ChunkSparsityPattern::create_from +(const unsigned int, + const unsigned int, + const CompressedSetSparsityPattern &, + const unsigned int, + const bool); +template +void ChunkSparsityPattern::create_from +(const unsigned int, + const unsigned int, + const CompressedSimpleSparsityPattern &, + const unsigned int, + const bool); +template void ChunkSparsityPattern::copy_from (const FullMatrix &, const unsigned int, const bool); diff --git a/tests/bits/chunk_sparse_matrix_11.cc b/tests/bits/chunk_sparse_matrix_11.cc new file mode 100644 index 0000000000..3018ad1fa8 --- /dev/null +++ b/tests/bits/chunk_sparse_matrix_11.cc @@ -0,0 +1,101 @@ +//----------------------- chunk_sparse_matrix_11.cc ------------------------ +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2004, 2005, 2008, 2013 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//----------------------- chunk_sparse_matrix_11.cc ------------------------ + + +// set a few elements in a chunk sparse matrix, read them via its iterator +// class and print them to the log file + +#include "../tests.h" +#include +#include +#include + + +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 m(sp); + + // first set a few entries + for (unsigned int i=0; i::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 +#include +#include + + +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 m(sp); + + // first set a few entries + for (unsigned int i=0; i::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 +#include +#include + + +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 m(sp); + + // first set a few entries + for (unsigned int i=0; i::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 +#include +#include + + +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 m(sp); + + // first set a few entries + for (unsigned int i=0; i::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 +#include +#include +#include + + +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 m(sp); + + // first set a few entries + for (unsigned int i=0; i src(m.n()), dst(m.m()), dst_ref(m.m()); + for (unsigned int i=0; i::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 +#include +#include + + +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 m(sp); + ChunkSparseMatrix::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 +#include +#include + + +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 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::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 +#include +#include + + +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 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::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 +#include +#include + + +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 m(sp); + + ChunkSparseMatrix::iterator i = m.begin(); + for (; i!=m.end(); ++i) + i->value() = i->row()*i->column(); + + for (unsigned int i=0; i<5; ++i) + for (unsigned int j=0; j<5; ++j) + if (((i+2*j+1) % 3 == 0) + || + (i==j)) + { + deallog << i << ' ' << j << ' ' << m.el(i,j) + << std::endl; + Assert (std::fabs(m.el(i,j)-i*j) < 1e-14, + ExcInternalError()); + } + + deallog << "OK" << std::endl; +} + + + +int main () +{ + std::ofstream logfile("chunk_sparse_matrix_iterator_04/output"); + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + try + { + test (1); + test (2); + test (4); + test (5); + test (7); + } + catch (std::exception &exc) + { + deallog << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + deallog << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; +} diff --git a/tests/bits/chunk_sparse_matrix_iterator_04/cmp/generic b/tests/bits/chunk_sparse_matrix_iterator_04/cmp/generic new file mode 100644 index 0000000000..5ec97e7639 --- /dev/null +++ b/tests/bits/chunk_sparse_matrix_iterator_04/cmp/generic @@ -0,0 +1,76 @@ + +DEAL::Chunk size: 1 +DEAL::0 0 0 +DEAL::0 1 0 +DEAL::0 4 0 +DEAL::1 1 1.00000 +DEAL::1 2 2.00000 +DEAL::2 0 0 +DEAL::2 2 4.00000 +DEAL::2 3 6.00000 +DEAL::3 1 3.00000 +DEAL::3 3 9.00000 +DEAL::3 4 12.0000 +DEAL::4 2 8.00000 +DEAL::4 4 16.0000 +DEAL::OK +DEAL::Chunk size: 2 +DEAL::0 0 0 +DEAL::0 1 0 +DEAL::0 4 0 +DEAL::1 1 1.00000 +DEAL::1 2 2.00000 +DEAL::2 0 0 +DEAL::2 2 4.00000 +DEAL::2 3 6.00000 +DEAL::3 1 3.00000 +DEAL::3 3 9.00000 +DEAL::3 4 12.0000 +DEAL::4 2 8.00000 +DEAL::4 4 16.0000 +DEAL::OK +DEAL::Chunk size: 4 +DEAL::0 0 0 +DEAL::0 1 0 +DEAL::0 4 0 +DEAL::1 1 1.00000 +DEAL::1 2 2.00000 +DEAL::2 0 0 +DEAL::2 2 4.00000 +DEAL::2 3 6.00000 +DEAL::3 1 3.00000 +DEAL::3 3 9.00000 +DEAL::3 4 12.0000 +DEAL::4 2 8.00000 +DEAL::4 4 16.0000 +DEAL::OK +DEAL::Chunk size: 5 +DEAL::0 0 0 +DEAL::0 1 0 +DEAL::0 4 0 +DEAL::1 1 1.00000 +DEAL::1 2 2.00000 +DEAL::2 0 0 +DEAL::2 2 4.00000 +DEAL::2 3 6.00000 +DEAL::3 1 3.00000 +DEAL::3 3 9.00000 +DEAL::3 4 12.0000 +DEAL::4 2 8.00000 +DEAL::4 4 16.0000 +DEAL::OK +DEAL::Chunk size: 7 +DEAL::0 0 0 +DEAL::0 1 0 +DEAL::0 4 0 +DEAL::1 1 1.00000 +DEAL::1 2 2.00000 +DEAL::2 0 0 +DEAL::2 2 4.00000 +DEAL::2 3 6.00000 +DEAL::3 1 3.00000 +DEAL::3 3 9.00000 +DEAL::3 4 12.0000 +DEAL::4 2 8.00000 +DEAL::4 4 16.0000 +DEAL::OK diff --git a/tests/bits/chunk_sparse_matrix_iterator_05.cc b/tests/bits/chunk_sparse_matrix_iterator_05.cc new file mode 100644 index 0000000000..011df54eac --- /dev/null +++ b/tests/bits/chunk_sparse_matrix_iterator_05.cc @@ -0,0 +1,103 @@ +//---------------------------- chunk_sparse_matrix_iterator_05.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2004, 2005, 2013 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------- chunk_sparse_matrix_iterator_05.cc --------------------------- + +// test setting some elements using a non-const chunk matrix iterator and +// operator+=, and reading them back through the matrix itself + +#include "../tests.h" +#include +#include +#include + + +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 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::iterator i = m.begin(); + for (; i!=m.end(); ++i) + i->value() += i->row()*i->column(); + + for (unsigned int i=0; i<5; ++i) + for (unsigned int j=0; j<5; ++j) + if (((i+2*j+1) % 3 == 0) + || + (i==j)) + { + deallog << i << ' ' << j << ' ' << m.el(i,j) + << std::endl; + Assert (std::fabs(m.el(i,j)-(1.+i*j)) < 1e-14, + ExcInternalError()); + } + + deallog << "OK" << std::endl; +} + + + +int main () +{ + std::ofstream logfile("chunk_sparse_matrix_iterator_05/output"); + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + try + { + test (1); + test (2); + test (4); + test (5); + test (7); + } + catch (std::exception &exc) + { + deallog << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + deallog << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; +} diff --git a/tests/bits/chunk_sparse_matrix_iterator_05/cmp/generic b/tests/bits/chunk_sparse_matrix_iterator_05/cmp/generic new file mode 100644 index 0000000000..94546302bc --- /dev/null +++ b/tests/bits/chunk_sparse_matrix_iterator_05/cmp/generic @@ -0,0 +1,76 @@ + +DEAL::Chunk size: 1 +DEAL::0 0 1.00000 +DEAL::0 1 1.00000 +DEAL::0 4 1.00000 +DEAL::1 1 2.00000 +DEAL::1 2 3.00000 +DEAL::2 0 1.00000 +DEAL::2 2 5.00000 +DEAL::2 3 7.00000 +DEAL::3 1 4.00000 +DEAL::3 3 10.0000 +DEAL::3 4 13.0000 +DEAL::4 2 9.00000 +DEAL::4 4 17.0000 +DEAL::OK +DEAL::Chunk size: 2 +DEAL::0 0 1.00000 +DEAL::0 1 1.00000 +DEAL::0 4 1.00000 +DEAL::1 1 2.00000 +DEAL::1 2 3.00000 +DEAL::2 0 1.00000 +DEAL::2 2 5.00000 +DEAL::2 3 7.00000 +DEAL::3 1 4.00000 +DEAL::3 3 10.0000 +DEAL::3 4 13.0000 +DEAL::4 2 9.00000 +DEAL::4 4 17.0000 +DEAL::OK +DEAL::Chunk size: 4 +DEAL::0 0 1.00000 +DEAL::0 1 1.00000 +DEAL::0 4 1.00000 +DEAL::1 1 2.00000 +DEAL::1 2 3.00000 +DEAL::2 0 1.00000 +DEAL::2 2 5.00000 +DEAL::2 3 7.00000 +DEAL::3 1 4.00000 +DEAL::3 3 10.0000 +DEAL::3 4 13.0000 +DEAL::4 2 9.00000 +DEAL::4 4 17.0000 +DEAL::OK +DEAL::Chunk size: 5 +DEAL::0 0 1.00000 +DEAL::0 1 1.00000 +DEAL::0 4 1.00000 +DEAL::1 1 2.00000 +DEAL::1 2 3.00000 +DEAL::2 0 1.00000 +DEAL::2 2 5.00000 +DEAL::2 3 7.00000 +DEAL::3 1 4.00000 +DEAL::3 3 10.0000 +DEAL::3 4 13.0000 +DEAL::4 2 9.00000 +DEAL::4 4 17.0000 +DEAL::OK +DEAL::Chunk size: 7 +DEAL::0 0 1.00000 +DEAL::0 1 1.00000 +DEAL::0 4 1.00000 +DEAL::1 1 2.00000 +DEAL::1 2 3.00000 +DEAL::2 0 1.00000 +DEAL::2 2 5.00000 +DEAL::2 3 7.00000 +DEAL::3 1 4.00000 +DEAL::3 3 10.0000 +DEAL::3 4 13.0000 +DEAL::4 2 9.00000 +DEAL::4 4 17.0000 +DEAL::OK diff --git a/tests/bits/chunk_sparse_matrix_iterator_06.cc b/tests/bits/chunk_sparse_matrix_iterator_06.cc new file mode 100644 index 0000000000..cb9cd38ab8 --- /dev/null +++ b/tests/bits/chunk_sparse_matrix_iterator_06.cc @@ -0,0 +1,110 @@ +//---------------------------- chunk_sparse_matrix_iterator_06.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2004, 2005, 2013 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------- chunk_sparse_matrix_iterator_06.cc --------------------------- + +// test setting some elements using a non-const chunk matrix iterator and +// operator-=, and reading them back through the matrix itself + +#include "../tests.h" +#include +#include +#include + + +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 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::iterator i = m.begin(); + for (; i!=m.end(); ++i) + i->value() -= i->row()*i->column(); + + for (unsigned int i=0; i<5; ++i) + for (unsigned int j=0; j<5; ++j) + if (((i+2*j+1) % 3 == 0) + || + (i==j)) + { + deallog << i << ' ' << j << ' ' << m.el(i,j) + << std::endl; + Assert (std::fabs(m.el(i,j)-(1.-i*j)) < 1e-14, + ExcInternalError()); + } + else if (sp.exists(i,j)) + { + deallog << i << ' ' << j << ' ' << m.el(i,j) + << std::endl; + Assert (std::fabs(m.el(i,j)+i*j) < 1e-14, + ExcInternalError()); + } + + deallog << "OK" << std::endl; +} + + + +int main () +{ + std::ofstream logfile("chunk_sparse_matrix_iterator_06/output"); + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + try + { + test (1); + test (2); + test (4); + test (5); + test (7); + } + catch (std::exception &exc) + { + deallog << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + deallog << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; +} diff --git a/tests/bits/chunk_sparse_matrix_iterator_06/cmp/generic b/tests/bits/chunk_sparse_matrix_iterator_06/cmp/generic new file mode 100644 index 0000000000..5654a87bac --- /dev/null +++ b/tests/bits/chunk_sparse_matrix_iterator_06/cmp/generic @@ -0,0 +1,122 @@ + +DEAL::Chunk size: 1 +DEAL::0 0 1.00000 +DEAL::0 1 1.00000 +DEAL::0 4 1.00000 +DEAL::1 1 0 +DEAL::1 2 -1.00000 +DEAL::2 0 1.00000 +DEAL::2 2 -3.00000 +DEAL::2 3 -5.00000 +DEAL::3 1 -2.00000 +DEAL::3 3 -8.00000 +DEAL::3 4 -11.0000 +DEAL::4 2 -7.00000 +DEAL::4 4 -15.0000 +DEAL::OK +DEAL::Chunk size: 2 +DEAL::0 0 1.00000 +DEAL::0 1 1.00000 +DEAL::0 2 0 +DEAL::0 3 0 +DEAL::0 4 1.00000 +DEAL::1 0 0 +DEAL::1 1 0 +DEAL::1 2 -1.00000 +DEAL::1 3 -3.00000 +DEAL::1 4 -4.00000 +DEAL::2 0 1.00000 +DEAL::2 1 -2.00000 +DEAL::2 2 -3.00000 +DEAL::2 3 -5.00000 +DEAL::2 4 -8.00000 +DEAL::3 0 0 +DEAL::3 1 -2.00000 +DEAL::3 2 -6.00000 +DEAL::3 3 -8.00000 +DEAL::3 4 -11.0000 +DEAL::4 2 -7.00000 +DEAL::4 3 -12.0000 +DEAL::4 4 -15.0000 +DEAL::OK +DEAL::Chunk size: 4 +DEAL::0 0 1.00000 +DEAL::0 1 1.00000 +DEAL::0 2 0 +DEAL::0 3 0 +DEAL::0 4 1.00000 +DEAL::1 0 0 +DEAL::1 1 0 +DEAL::1 2 -1.00000 +DEAL::1 3 -3.00000 +DEAL::1 4 -4.00000 +DEAL::2 0 1.00000 +DEAL::2 1 -2.00000 +DEAL::2 2 -3.00000 +DEAL::2 3 -5.00000 +DEAL::2 4 -8.00000 +DEAL::3 0 0 +DEAL::3 1 -2.00000 +DEAL::3 2 -6.00000 +DEAL::3 3 -8.00000 +DEAL::3 4 -11.0000 +DEAL::4 0 0 +DEAL::4 1 -4.00000 +DEAL::4 2 -7.00000 +DEAL::4 3 -12.0000 +DEAL::4 4 -15.0000 +DEAL::OK +DEAL::Chunk size: 5 +DEAL::0 0 1.00000 +DEAL::0 1 1.00000 +DEAL::0 2 0 +DEAL::0 3 0 +DEAL::0 4 1.00000 +DEAL::1 0 0 +DEAL::1 1 0 +DEAL::1 2 -1.00000 +DEAL::1 3 -3.00000 +DEAL::1 4 -4.00000 +DEAL::2 0 1.00000 +DEAL::2 1 -2.00000 +DEAL::2 2 -3.00000 +DEAL::2 3 -5.00000 +DEAL::2 4 -8.00000 +DEAL::3 0 0 +DEAL::3 1 -2.00000 +DEAL::3 2 -6.00000 +DEAL::3 3 -8.00000 +DEAL::3 4 -11.0000 +DEAL::4 0 0 +DEAL::4 1 -4.00000 +DEAL::4 2 -7.00000 +DEAL::4 3 -12.0000 +DEAL::4 4 -15.0000 +DEAL::OK +DEAL::Chunk size: 7 +DEAL::0 0 1.00000 +DEAL::0 1 1.00000 +DEAL::0 2 0 +DEAL::0 3 0 +DEAL::0 4 1.00000 +DEAL::1 0 0 +DEAL::1 1 0 +DEAL::1 2 -1.00000 +DEAL::1 3 -3.00000 +DEAL::1 4 -4.00000 +DEAL::2 0 1.00000 +DEAL::2 1 -2.00000 +DEAL::2 2 -3.00000 +DEAL::2 3 -5.00000 +DEAL::2 4 -8.00000 +DEAL::3 0 0 +DEAL::3 1 -2.00000 +DEAL::3 2 -6.00000 +DEAL::3 3 -8.00000 +DEAL::3 4 -11.0000 +DEAL::4 0 0 +DEAL::4 1 -4.00000 +DEAL::4 2 -7.00000 +DEAL::4 3 -12.0000 +DEAL::4 4 -15.0000 +DEAL::OK diff --git a/tests/bits/chunk_sparse_matrix_iterator_07.cc b/tests/bits/chunk_sparse_matrix_iterator_07.cc new file mode 100644 index 0000000000..c380b20ea7 --- /dev/null +++ b/tests/bits/chunk_sparse_matrix_iterator_07.cc @@ -0,0 +1,110 @@ +//---------------------------- chunk_sparse_matrix_iterator_07.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2004, 2005, 2013 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------- chunk_sparse_matrix_iterator_07.cc --------------------------- + +// test setting some elements using a non-const chunk matrix iterator and +// operator*=, and reading them back through the matrix itself + +#include "../tests.h" +#include +#include +#include + + +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 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::iterator i = m.begin(); + for (; i!=m.end(); ++i) + i->value() *= 2; + + for (unsigned int i=0; i<5; ++i) + for (unsigned int j=0; j<5; ++j) + if (((i+2*j+1) % 3 == 0) + || + (i==j)) + { + deallog << i << ' ' << j << ' ' << m.el(i,j) + << std::endl; + Assert (std::fabs(m.el(i,j)-2*i*j) < 1e-14, + ExcInternalError()); + } + else if (sp.exists(i,j)) + { + deallog << i << ' ' << j << ' ' << m.el(i,j) + << std::endl; + Assert (std::fabs(m.el(i,j)) < 1e-14, + ExcInternalError()); + } + + deallog << "OK" << std::endl; +} + + + +int main () +{ + std::ofstream logfile("chunk_sparse_matrix_iterator_07/output"); + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + try + { + test (1); + test (2); + test (4); + test (5); + test (7); + } + catch (std::exception &exc) + { + deallog << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + deallog << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; +} diff --git a/tests/bits/chunk_sparse_matrix_iterator_07/cmp/generic b/tests/bits/chunk_sparse_matrix_iterator_07/cmp/generic new file mode 100644 index 0000000000..3b469fb493 --- /dev/null +++ b/tests/bits/chunk_sparse_matrix_iterator_07/cmp/generic @@ -0,0 +1,122 @@ + +DEAL::Chunk size: 1 +DEAL::0 0 0 +DEAL::0 1 0 +DEAL::0 4 0 +DEAL::1 1 2.00000 +DEAL::1 2 4.00000 +DEAL::2 0 0 +DEAL::2 2 8.00000 +DEAL::2 3 12.0000 +DEAL::3 1 6.00000 +DEAL::3 3 18.0000 +DEAL::3 4 24.0000 +DEAL::4 2 16.0000 +DEAL::4 4 32.0000 +DEAL::OK +DEAL::Chunk size: 2 +DEAL::0 0 0 +DEAL::0 1 0 +DEAL::0 2 0 +DEAL::0 3 0 +DEAL::0 4 0 +DEAL::1 0 0 +DEAL::1 1 2.00000 +DEAL::1 2 4.00000 +DEAL::1 3 0 +DEAL::1 4 0 +DEAL::2 0 0 +DEAL::2 1 0 +DEAL::2 2 8.00000 +DEAL::2 3 12.0000 +DEAL::2 4 0 +DEAL::3 0 0 +DEAL::3 1 6.00000 +DEAL::3 2 0 +DEAL::3 3 18.0000 +DEAL::3 4 24.0000 +DEAL::4 2 16.0000 +DEAL::4 3 0 +DEAL::4 4 32.0000 +DEAL::OK +DEAL::Chunk size: 4 +DEAL::0 0 0 +DEAL::0 1 0 +DEAL::0 2 0 +DEAL::0 3 0 +DEAL::0 4 0 +DEAL::1 0 0 +DEAL::1 1 2.00000 +DEAL::1 2 4.00000 +DEAL::1 3 0 +DEAL::1 4 0 +DEAL::2 0 0 +DEAL::2 1 0 +DEAL::2 2 8.00000 +DEAL::2 3 12.0000 +DEAL::2 4 0 +DEAL::3 0 0 +DEAL::3 1 6.00000 +DEAL::3 2 0 +DEAL::3 3 18.0000 +DEAL::3 4 24.0000 +DEAL::4 0 0 +DEAL::4 1 0 +DEAL::4 2 16.0000 +DEAL::4 3 0 +DEAL::4 4 32.0000 +DEAL::OK +DEAL::Chunk size: 5 +DEAL::0 0 0 +DEAL::0 1 0 +DEAL::0 2 0 +DEAL::0 3 0 +DEAL::0 4 0 +DEAL::1 0 0 +DEAL::1 1 2.00000 +DEAL::1 2 4.00000 +DEAL::1 3 0 +DEAL::1 4 0 +DEAL::2 0 0 +DEAL::2 1 0 +DEAL::2 2 8.00000 +DEAL::2 3 12.0000 +DEAL::2 4 0 +DEAL::3 0 0 +DEAL::3 1 6.00000 +DEAL::3 2 0 +DEAL::3 3 18.0000 +DEAL::3 4 24.0000 +DEAL::4 0 0 +DEAL::4 1 0 +DEAL::4 2 16.0000 +DEAL::4 3 0 +DEAL::4 4 32.0000 +DEAL::OK +DEAL::Chunk size: 7 +DEAL::0 0 0 +DEAL::0 1 0 +DEAL::0 2 0 +DEAL::0 3 0 +DEAL::0 4 0 +DEAL::1 0 0 +DEAL::1 1 2.00000 +DEAL::1 2 4.00000 +DEAL::1 3 0 +DEAL::1 4 0 +DEAL::2 0 0 +DEAL::2 1 0 +DEAL::2 2 8.00000 +DEAL::2 3 12.0000 +DEAL::2 4 0 +DEAL::3 0 0 +DEAL::3 1 6.00000 +DEAL::3 2 0 +DEAL::3 3 18.0000 +DEAL::3 4 24.0000 +DEAL::4 0 0 +DEAL::4 1 0 +DEAL::4 2 16.0000 +DEAL::4 3 0 +DEAL::4 4 32.0000 +DEAL::OK diff --git a/tests/bits/chunk_sparse_matrix_iterator_08.cc b/tests/bits/chunk_sparse_matrix_iterator_08.cc new file mode 100644 index 0000000000..5c59cd3dc1 --- /dev/null +++ b/tests/bits/chunk_sparse_matrix_iterator_08.cc @@ -0,0 +1,110 @@ +//---------------------------- chunk_sparse_matrix_iterator_08.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2004, 2005, 2013 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------- chunk_sparse_matrix_iterator_08.cc --------------------------- + +// test setting some elements using a non-const chunk matrix iterator and +// operator/=, and reading them back through the matrix itself + +#include "../tests.h" +#include +#include +#include + + +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 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::iterator i = m.begin(); + for (; i!=m.end(); ++i) + i->value() /= 2; + + for (unsigned int i=0; i<5; ++i) + for (unsigned int j=0; j<5; ++j) + if (((i+2*j+1) % 3 == 0) + || + (i==j)) + { + deallog << i << ' ' << j << ' ' << m.el(i,j) + << std::endl; + Assert (std::fabs(m.el(i,j)-(i*j/2.)) < 1e-14, + ExcInternalError()); + } + else if (sp.exists(i,j)) + { + deallog << i << ' ' << j << ' ' << m.el(i,j) + << std::endl; + Assert (std::fabs(m.el(i,j)) < 1e-14, + ExcInternalError()); + } + + deallog << "OK" << std::endl; +} + + + +int main () +{ + std::ofstream logfile("chunk_sparse_matrix_iterator_08/output"); + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + try + { + test (1); + test (2); + test (4); + test (5); + test (7); + } + catch (std::exception &exc) + { + deallog << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + deallog << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; +} diff --git a/tests/bits/chunk_sparse_matrix_iterator_08/cmp/generic b/tests/bits/chunk_sparse_matrix_iterator_08/cmp/generic new file mode 100644 index 0000000000..38a9f8ac2e --- /dev/null +++ b/tests/bits/chunk_sparse_matrix_iterator_08/cmp/generic @@ -0,0 +1,122 @@ + +DEAL::Chunk size: 1 +DEAL::0 0 0 +DEAL::0 1 0 +DEAL::0 4 0 +DEAL::1 1 0.500000 +DEAL::1 2 1.00000 +DEAL::2 0 0 +DEAL::2 2 2.00000 +DEAL::2 3 3.00000 +DEAL::3 1 1.50000 +DEAL::3 3 4.50000 +DEAL::3 4 6.00000 +DEAL::4 2 4.00000 +DEAL::4 4 8.00000 +DEAL::OK +DEAL::Chunk size: 2 +DEAL::0 0 0 +DEAL::0 1 0 +DEAL::0 2 0 +DEAL::0 3 0 +DEAL::0 4 0 +DEAL::1 0 0 +DEAL::1 1 0.500000 +DEAL::1 2 1.00000 +DEAL::1 3 0 +DEAL::1 4 0 +DEAL::2 0 0 +DEAL::2 1 0 +DEAL::2 2 2.00000 +DEAL::2 3 3.00000 +DEAL::2 4 0 +DEAL::3 0 0 +DEAL::3 1 1.50000 +DEAL::3 2 0 +DEAL::3 3 4.50000 +DEAL::3 4 6.00000 +DEAL::4 2 4.00000 +DEAL::4 3 0 +DEAL::4 4 8.00000 +DEAL::OK +DEAL::Chunk size: 4 +DEAL::0 0 0 +DEAL::0 1 0 +DEAL::0 2 0 +DEAL::0 3 0 +DEAL::0 4 0 +DEAL::1 0 0 +DEAL::1 1 0.500000 +DEAL::1 2 1.00000 +DEAL::1 3 0 +DEAL::1 4 0 +DEAL::2 0 0 +DEAL::2 1 0 +DEAL::2 2 2.00000 +DEAL::2 3 3.00000 +DEAL::2 4 0 +DEAL::3 0 0 +DEAL::3 1 1.50000 +DEAL::3 2 0 +DEAL::3 3 4.50000 +DEAL::3 4 6.00000 +DEAL::4 0 0 +DEAL::4 1 0 +DEAL::4 2 4.00000 +DEAL::4 3 0 +DEAL::4 4 8.00000 +DEAL::OK +DEAL::Chunk size: 5 +DEAL::0 0 0 +DEAL::0 1 0 +DEAL::0 2 0 +DEAL::0 3 0 +DEAL::0 4 0 +DEAL::1 0 0 +DEAL::1 1 0.500000 +DEAL::1 2 1.00000 +DEAL::1 3 0 +DEAL::1 4 0 +DEAL::2 0 0 +DEAL::2 1 0 +DEAL::2 2 2.00000 +DEAL::2 3 3.00000 +DEAL::2 4 0 +DEAL::3 0 0 +DEAL::3 1 1.50000 +DEAL::3 2 0 +DEAL::3 3 4.50000 +DEAL::3 4 6.00000 +DEAL::4 0 0 +DEAL::4 1 0 +DEAL::4 2 4.00000 +DEAL::4 3 0 +DEAL::4 4 8.00000 +DEAL::OK +DEAL::Chunk size: 7 +DEAL::0 0 0 +DEAL::0 1 0 +DEAL::0 2 0 +DEAL::0 3 0 +DEAL::0 4 0 +DEAL::1 0 0 +DEAL::1 1 0.500000 +DEAL::1 2 1.00000 +DEAL::1 3 0 +DEAL::1 4 0 +DEAL::2 0 0 +DEAL::2 1 0 +DEAL::2 2 2.00000 +DEAL::2 3 3.00000 +DEAL::2 4 0 +DEAL::3 0 0 +DEAL::3 1 1.50000 +DEAL::3 2 0 +DEAL::3 3 4.50000 +DEAL::3 4 6.00000 +DEAL::4 0 0 +DEAL::4 1 0 +DEAL::4 2 4.00000 +DEAL::4 3 0 +DEAL::4 4 8.00000 +DEAL::OK diff --git a/tests/bits/chunk_sparse_matrix_iterator_09.cc b/tests/bits/chunk_sparse_matrix_iterator_09.cc new file mode 100644 index 0000000000..6edddd9e6b --- /dev/null +++ b/tests/bits/chunk_sparse_matrix_iterator_09.cc @@ -0,0 +1,83 @@ +//---------------------------- chunk_sparse_matrix_iterator_09.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2004, 2005, 2013 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------- chunk_sparse_matrix_iterator_09.cc --------------------------- + + +// this test is sparse_matrix_iterator_09 for a ChunkSparseMatrix + +#include "../tests.h" +#include +#include +#include + + +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 A(sparsity); + + // and loop over the elements of it + for (ChunkSparseMatrix::const_iterator k=A.begin(); + k!=A.end(); ++k) + deallog << k->row() << ' ' << k->column() << ' ' << k->value() + << std::endl; +} + + + +int main () +{ + std::ofstream logfile("chunk_sparse_matrix_iterator_09/output"); + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + try + { + test (1); + test (3); + } + catch (std::exception &exc) + { + deallog << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + deallog << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; +} diff --git a/tests/bits/chunk_sparse_matrix_iterator_09/cmp/generic b/tests/bits/chunk_sparse_matrix_iterator_09/cmp/generic new file mode 100644 index 0000000000..3d7d6ec0c4 --- /dev/null +++ b/tests/bits/chunk_sparse_matrix_iterator_09/cmp/generic @@ -0,0 +1,19 @@ + +DEAL::Chunk size: 1 +DEAL::1 1 0 +DEAL::3 1 0 +DEAL::Chunk size: 3 +DEAL::0 0 0 +DEAL::0 1 0 +DEAL::0 2 0 +DEAL::1 0 0 +DEAL::1 1 0 +DEAL::1 2 0 +DEAL::2 0 0 +DEAL::2 1 0 +DEAL::2 2 0 +DEAL::3 3 0 +DEAL::3 4 0 +DEAL::3 0 0 +DEAL::3 1 0 +DEAL::3 2 0 diff --git a/tests/bits/chunk_sparse_matrix_iterator_10.cc b/tests/bits/chunk_sparse_matrix_iterator_10.cc new file mode 100644 index 0000000000..8f788d1f58 --- /dev/null +++ b/tests/bits/chunk_sparse_matrix_iterator_10.cc @@ -0,0 +1,85 @@ +//---------------------------- chunk_sparse_matrix_iterator_10.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2004, 2005, 2013 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------- chunk_sparse_matrix_iterator_10.cc --------------------------- + + +// this test is sparse_matrix_iterator_10 for a ChunkSparseMatrix and the same +// test as chunk_sparse_matrix_iterator_09 with postfix operator++ instead of +// prefix + +#include "../tests.h" +#include +#include +#include + + +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 A(sparsity); + + // and loop over the elements of it + for (ChunkSparseMatrix::const_iterator k=A.begin(); + k!=A.end(); ++k) + deallog << k->row() << ' ' << k->column() << ' ' << k->value() + << std::endl; +} + + + +int main () +{ + std::ofstream logfile("chunk_sparse_matrix_iterator_10/output"); + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + try + { + test (1); + test (3); + } + catch (std::exception &exc) + { + deallog << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + deallog << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; +} diff --git a/tests/bits/chunk_sparse_matrix_iterator_10/cmp/generic b/tests/bits/chunk_sparse_matrix_iterator_10/cmp/generic new file mode 100644 index 0000000000..3d7d6ec0c4 --- /dev/null +++ b/tests/bits/chunk_sparse_matrix_iterator_10/cmp/generic @@ -0,0 +1,19 @@ + +DEAL::Chunk size: 1 +DEAL::1 1 0 +DEAL::3 1 0 +DEAL::Chunk size: 3 +DEAL::0 0 0 +DEAL::0 1 0 +DEAL::0 2 0 +DEAL::1 0 0 +DEAL::1 1 0 +DEAL::1 2 0 +DEAL::2 0 0 +DEAL::2 1 0 +DEAL::2 2 0 +DEAL::3 3 0 +DEAL::3 4 0 +DEAL::3 0 0 +DEAL::3 1 0 +DEAL::3 2 0 diff --git a/tests/bits/chunk_sparse_matrix_iterator_11.cc b/tests/bits/chunk_sparse_matrix_iterator_11.cc new file mode 100644 index 0000000000..d5d6560a11 --- /dev/null +++ b/tests/bits/chunk_sparse_matrix_iterator_11.cc @@ -0,0 +1,95 @@ +//---------------------------- chunk_sparse_matrix_iterator_11.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2004, 2005, 2013 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------- chunk_sparse_matrix_iterator_11.cc --------------------------- + + +// comparisons between chunk sparse matrix iterators, same as +// sparse_matrix_iterators_11 otherwise. + +#include "../tests.h" +#include +#include +#include + + +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 A(sparsity); + + ChunkSparseMatrix::iterator k = A.begin(), + j = ++A.begin(); + + Assert (k < j, ExcInternalError()); + Assert (j > k, ExcInternalError()); + + Assert (!(j < k), ExcInternalError()); + Assert (!(k > j), ExcInternalError()); + + Assert (k != j, ExcInternalError()); + Assert (!(k == j), ExcInternalError()); + + Assert (k == k, ExcInternalError()); + Assert (!(k != k), ExcInternalError()); + + deallog << "OK" << std::endl; +} + + + +int main () +{ + std::ofstream logfile("chunk_sparse_matrix_iterator_11/output"); + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + try + { + test (1); + test (3); + } + catch (std::exception &exc) + { + deallog << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + deallog << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; +} diff --git a/tests/bits/chunk_sparse_matrix_iterator_11/cmp/generic b/tests/bits/chunk_sparse_matrix_iterator_11/cmp/generic new file mode 100644 index 0000000000..764c6e2ec6 --- /dev/null +++ b/tests/bits/chunk_sparse_matrix_iterator_11/cmp/generic @@ -0,0 +1,5 @@ + +DEAL::Chunk size: 1 +DEAL::OK +DEAL::Chunk size: 3 +DEAL::OK diff --git a/tests/bits/chunk_sparse_matrix_iterator_12.cc b/tests/bits/chunk_sparse_matrix_iterator_12.cc new file mode 100644 index 0000000000..63032cd162 --- /dev/null +++ b/tests/bits/chunk_sparse_matrix_iterator_12.cc @@ -0,0 +1,94 @@ +//---------------------------- chunk_sparse_matrix_iterator_12.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2004, 2005, 2013 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------- chunk_sparse_matrix_iterator_12.cc --------------------------- + + +// like chunk_sparse_matrix_iterator_11 but for const_iterator + +#include "../tests.h" +#include +#include +#include + + +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 A(sparsity); + + ChunkSparseMatrix::const_iterator k = A.begin(), + j = ++A.begin(); + + Assert (k < j, ExcInternalError()); + Assert (j > k, ExcInternalError()); + + Assert (!(j < k), ExcInternalError()); + Assert (!(k > j), ExcInternalError()); + + Assert (k != j, ExcInternalError()); + Assert (!(k == j), ExcInternalError()); + + Assert (k == k, ExcInternalError()); + Assert (!(k != k), ExcInternalError()); + + deallog << "OK" << std::endl; +} + + + +int main () +{ + std::ofstream logfile("chunk_sparse_matrix_iterator_12/output"); + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + try + { + test (1); + test (3); + } + catch (std::exception &exc) + { + deallog << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + deallog << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + deallog << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; +} diff --git a/tests/bits/chunk_sparse_matrix_iterator_12/cmp/generic b/tests/bits/chunk_sparse_matrix_iterator_12/cmp/generic new file mode 100644 index 0000000000..764c6e2ec6 --- /dev/null +++ b/tests/bits/chunk_sparse_matrix_iterator_12/cmp/generic @@ -0,0 +1,5 @@ + +DEAL::Chunk size: 1 +DEAL::OK +DEAL::Chunk size: 3 +DEAL::OK diff --git a/tests/bits/chunk_sparse_matrix_iterator_13.cc b/tests/bits/chunk_sparse_matrix_iterator_13.cc new file mode 100644 index 0000000000..8c38ca9d63 --- /dev/null +++ b/tests/bits/chunk_sparse_matrix_iterator_13.cc @@ -0,0 +1,117 @@ +//---------------------------- chunk_sparse_matrix_iterator_13.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2004, 2005, 2013 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------- chunk_sparse_matrix_iterator_13.cc --------------------------- + + +// test ChunkSparseMatrix::iterator::operator- + +#include "../tests.h" +#include +#include +#include + + +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 m(sp); + + for (unsigned int row=0; row