From: bangerth Date: Sat, 2 Feb 2013 21:35:22 +0000 (+0000) Subject: Always use the storage pattern where the diagonal entry comes first in each row,... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=c8a2db377fc0c4c3f56ca2744c1e812399b4953d;p=dealii-svn.git Always use the storage pattern where the diagonal entry comes first in each row, as long as the matrix is square. git-svn-id: https://svn.dealii.org/trunk@28216 0785d39b-7218-0410-832d-ea1e28bc413d --- 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 8a41854df8..24115a9238 100644 --- a/deal.II/include/deal.II/lac/chunk_sparsity_pattern.h +++ b/deal.II/include/deal.II/lac/chunk_sparsity_pattern.h @@ -884,7 +884,7 @@ inline bool ChunkSparsityPattern::optimize_diagonal () const { - return sparsity_pattern.diagonal_optimized; + return sparsity_pattern.store_diagonal_first_in_row; } diff --git a/deal.II/include/deal.II/lac/sparse_matrix.h b/deal.II/include/deal.II/lac/sparse_matrix.h index f2498fb6ec..9c8527715a 100644 --- a/deal.II/include/deal.II/lac/sparse_matrix.h +++ b/deal.II/include/deal.II/lac/sparse_matrix.h @@ -278,6 +278,13 @@ namespace SparseMatrixIterators /** * STL conforming iterator for constant and non-constant matrices. * + * The typical use for these iterators is to iterate over the elements of + * a sparse matrix + * or over the elements of individual rows. Note that there is no guarantee + * that the elements of a row are actually traversed in an order in which + * columns monotonically increase. See the documentation of the + * SparsityPattern class for more information. + * * The first template argument denotes the underlying numeric type, the * second the constness of the matrix. * @@ -401,6 +408,18 @@ namespace SparseMatrixIterators * importantly one can associate more than one matrix with the same * sparsity pattern. * + * The elements of a SparseMatrix are stored in the same order in which the + * SparsityPattern class stores its entries. + * Within each row, elements are generally stored left-to-right in increasing + * column index order; the exception to this rule is that if the matrix + * is square (n_rows() == n_columns()), then the diagonal entry is stored + * as the first element in each row to make operations like applying a + * Jacobi or SSOR preconditioner faster. As a consequence, if you traverse + * the elements of a row of a SparseMatrix with the help of iterators into + * this object (using SparseMatrix::begin and SparseMatrix::end) you + * will find that the elements are not sorted by column index within each row + * whenever the matrix is square. + * * @note Instantiations for this template are provided for @ and * @; others can be generated in application programs (see the * section on @ref Instantiations in the manual). @@ -1299,6 +1318,9 @@ public: /** * STL-like iterator with the first entry of the matrix. This is the version * for constant matrices. + * + * Note the discussion in the general documentation of this class about + * the order in which elements are accessed. */ const_iterator begin () const; @@ -1310,6 +1332,9 @@ public: /** * STL-like iterator with the first entry of the matrix. This is the version * for non-constant matrices. + * + * Note the discussion in the general documentation of this class about + * the order in which elements are accessed. */ iterator begin (); @@ -1326,6 +1351,9 @@ public: * entries, then the iterator returned by this function equals * end(r). Note also that the iterator may not be dereferencable in * that case. + * + * Note also the discussion in the general documentation of this class about + * the order in which elements are accessed. */ const_iterator begin (const unsigned int r) const; @@ -1348,6 +1376,9 @@ public: * entries, then the iterator returned by this function equals * end(r). Note also that the iterator may not be dereferencable in * that case. + * + * Note the discussion in the general documentation of this class about + * the order in which elements are accessed. */ iterator begin (const unsigned int r); diff --git a/deal.II/include/deal.II/lac/sparsity_pattern.h b/deal.II/include/deal.II/lac/sparsity_pattern.h index f0f1a1490e..16f30df171 100644 --- a/deal.II/include/deal.II/lac/sparsity_pattern.h +++ b/deal.II/include/deal.II/lac/sparsity_pattern.h @@ -184,6 +184,14 @@ namespace SparsityPatternIterators /** * STL conforming iterator walking over the elements of a sparsity pattern. + * + * The typical use for these iterators is to iterate over the elements of + * a sparsity pattern (or, since they also serve as the basis for iterating + * over the elements of an associated matrix, over the elements of a sparse matrix) + * or over the elements of individual rows. Note that there is no guarantee + * that the elements of a row are actually traversed in an order in which + * columns monotonically increase. See the documentation of the + * SparsityPattern class for more information. */ class Iterator { @@ -247,11 +255,28 @@ namespace SparsityPatternIterators /** * Structure representing the sparsity pattern of a sparse matrix. - * * This class is an example of the "static" type of @ref Sparsity. - * * It uses the compressed row storage (CSR) format to store data. * + * The elements of a SparsityPattern, corresponding to the places where + * SparseMatrix objects can store nonzero entries, are stored row-by-row. + * Within each row, elements are generally stored left-to-right in increasing + * column index order; the exception to this rule is that if the matrix + * is square (n_rows() == n_columns()), then the diagonal entry is stored + * as the first element in each row to make operations like applying a + * Jacobi or SSOR preconditioner faster. As a consequence, if you traverse + * the elements of a row of a SparsityPattern with the help of iterators into + * this object (using SparsityPattern::begin and SparsityPattern::end) you + * will find that the elements are not sorted by column index within each row + * whenever the matrix is square. + * + * @note While this class forms the basis upon which SparseMatrix + * objects base their storage format, and thus plays a central role in + * setting up linear systems, it is rarely set up directly due to the + * way it stores its information. Rather, one typically goes through + * an intermediate format first, see for example the step-2 tutorial + * program as well as the documentation module @ref Sparsity. + * * @author Wolfgang Bangerth, Guido Kanschat and others */ class SparsityPattern : public Subscriptor @@ -268,6 +293,8 @@ public: /** * Typedef an iterator class that allows to walk over the nonzero elements * of a row of a sparsity pattern. + * + * @deprecated This typedef is deprecated. Use proper iterators instead. */ typedef const unsigned int *row_iterator; @@ -335,13 +362,26 @@ public: * @arg m number of rows * @arg n number of columns * @arg max_per_row maximum number of nonzero entries per row - * @arg optimize_diagonal store diagonal entries first in row; see - * optimize_diagonal(). This takes effect for quadratic matrices only. + * @arg optimize_diagonal This parameter is ignored. + * + * @deprecated Use the constructor without the last argument since + * it is ignored. */ SparsityPattern (const unsigned int m, const unsigned int n, const unsigned int max_per_row, - const bool optimize_diagonal = true); + const bool optimize_diagonal) DEAL_II_DEPRECATED; + + /** + * 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 + */ + SparsityPattern (const unsigned int m, + const unsigned int n, + const unsigned int max_per_row); /** * Initialize a rectangular matrix. @@ -350,13 +390,27 @@ public: * @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 optimize_diagonal store diagonal entries first in row; see - * optimize_diagonal(). This takes effect for quadratic matrices only. + * @arg optimize_diagonal This argument is ignored. + * + * @deprecated Use the constructor without the last argument since + * it is ignored. */ SparsityPattern (const unsigned int m, const unsigned int n, const std::vector &row_lengths, - const bool optimize_diagonal = true); + const bool optimize_diagonal) DEAL_II_DEPRECATED; + + /** + * 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. + */ + SparsityPattern (const unsigned int m, + const unsigned int n, + const std::vector &row_lengths); /** * Initialize a quadratic matrix of dimension n with at most @@ -375,12 +429,24 @@ public: * @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 optimize_diagonal store diagonal entries first in row; see - * optimize_diagonal(). + * @arg optimize_diagonal This argument is ignored. + * + * @deprecated Use the constructor without the last argument since + * it is ignored. */ SparsityPattern (const unsigned int m, const std::vector &row_lengths, - const bool optimize_diagonal = true); + const bool optimize_diagonal) DEAL_II_DEPRECATED; + + /** + * 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. + */ + SparsityPattern (const unsigned int m, + const std::vector &row_lengths); /** * Make a copy with extra off-diagonals. @@ -427,11 +493,26 @@ public: * * This function simply maps its operations to the other reinit * function. + * + * @deprecated The last argument of this function is ignored. Use the + * version of this function without the last argument. */ void reinit (const unsigned int m, const unsigned int n, const unsigned int max_per_row, - const bool optimize_diagonal = true); + const bool optimize_diagonal) DEAL_II_DEPRECATED; + + /** + * 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. + */ + void reinit (const unsigned int m, + const unsigned int n, + const unsigned int max_per_row); /** * Reallocate memory for a matrix of size m x n. The number of @@ -446,19 +527,50 @@ public: * If the number of rows equals the number of columns and the last parameter * is true, diagonal elements are stored first in each row to allow * optimized access in relaxation methods of SparseMatrix. + * + * @deprecated The last argument of this function is ignored. Use the + * version of this function without the last argument. */ void reinit (const unsigned int m, const unsigned int n, const std::vector &row_lengths, - const bool optimize_diagonal = true); + 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. + * + * 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 and the last parameter + * is true, diagonal elements are stored first in each row to allow + * optimized access in relaxation methods of SparseMatrix. + */ + void reinit (const unsigned int m, + const unsigned int n, + const std::vector &row_lengths); /** * Same as above, but with a VectorSlice argument instead. + * + * @deprecated The last argument of this function is ignored. Use the + * version of this function without the last argument. */ void reinit (const unsigned int m, const unsigned int n, const VectorSlice > &row_lengths, - const bool optimize_diagonal = true); + const bool optimize_diagonal) DEAL_II_DEPRECATED; + + /** + * Same as above, but with a VectorSlice argument instead. + */ + void reinit (const unsigned int m, + const unsigned int n, + const VectorSlice > &row_lengths); /** * This function compresses the sparsity structure that this object @@ -549,23 +661,120 @@ public: * int,double> could be replaced by std::vector >, or a list or set of such pairs, as they all return * iterators that point to such pairs. + * + * @deprecated The last argument of this function is ignored. Use the + * version of this function without the last argument. */ template void copy_from (const unsigned int n_rows, const unsigned int n_cols, const ForwardIterator begin, const ForwardIterator end, - const bool optimize_diagonal = true); + const bool optimize_diagonal) DEAL_II_DEPRECATED; + + /** + * 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: + * @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. + */ + template + void copy_from (const unsigned int n_rows, + const unsigned int n_cols, + const ForwardIterator begin, + const ForwardIterator end); /** * 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. + * + * @deprecated The last argument of this function is ignored. Use the + * version of this function without the last argument. */ template void copy_from (const CompressedSparsityType &csp, - const bool optimize_diagonal = true); + 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 + * compressed mode afterwards. + */ + template + void copy_from (const CompressedSparsityType &csp); /** * Take a full matrix and use its nonzero entries to generate a sparse @@ -573,10 +782,23 @@ public: * * Previous content of this object is lost, and the sparsity pattern is in * compressed mode afterwards. + * + * @deprecated The last argument of this function is ignored. Use the + * version of this function without the last argument. */ template void copy_from (const FullMatrix &matrix, - const bool optimize_diagonal = true); + 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. + * + * Previous content of this object is lost, and the sparsity pattern is in + * compressed mode afterwards. + */ + template + void copy_from (const FullMatrix &matrix); /** * Make the sparsity pattern symmetric by adding the sparsity pattern of the @@ -618,6 +840,9 @@ public: * 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. + * + * Note the discussion in the general documentation of this class about + * the order in which elements are accessed. */ iterator begin () const; @@ -633,6 +858,9 @@ public: * entries, then the iterator returned by this function equals * end(r). Note also that the iterator may not be dereferencable in * that case. + * + * Note also the discussion in the general documentation of this class about + * the order in which elements are accessed. */ iterator begin (const unsigned int r) const; @@ -653,8 +881,11 @@ public: * entries, then the iterator returned by this function equals * end(r). Note also that the iterator may not be dereferencable in * that case. + * + * @deprecated Use the iterators provided by the begin() and end() + * functions instead. */ - row_iterator row_begin (const unsigned int r) const; + row_iterator row_begin (const unsigned int r) const DEAL_II_DEPRECATED; /** * Final iterator of row r. It points to the first element past the @@ -663,8 +894,11 @@ public: * 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. + * + * @deprecated Use the iterators provided by the begin() and end() + * functions instead. */ - row_iterator row_end (const unsigned int r) const; + row_iterator row_end (const unsigned int r) const DEAL_II_DEPRECATED; // @} /** @@ -730,34 +964,20 @@ public: unsigned int row_length (const unsigned int row) const; /** - * Determine whether the matrix uses special convention for quadratic - * matrices. - * - * A return value true means that diagonal elements are stored - * first in each row. A number of functions in this class and the library in - * general, for example relaxation methods like Jacobi() and SOR(), require - * this to make their operations more efficient, since they need to quickly - * access the diagonal elements and do not have to search for them if they - * are the first element of each row. A side effect of this scheme is that - * each row contains at least one element, even if the row is empty - * (i.e. the diagonal element exists, but has value zero). + * Determine whether the matrix uses the special convention for quadratic + * matrices that the diagonal entries are stored first in each row. * - * A return value false means that diagonal elements are stored - * anywhere in the row, or not at all. In particular, a row or even the - * whole matrix may be empty. This can be used if you have block matrices - * where the off-diagonal blocks are quadratic but are never used for - * operations like the ones mentioned above. In this case, some memory can - * be saved by not using the diagonal storage optimization. + * @deprecated The function always returns true if the matrix is + * square and false if it is not. */ - bool optimize_diagonal () const; + bool optimize_diagonal () const DEAL_II_DEPRECATED; /** * Return whether this object stores only those entries that have been added * explicitly, or if the sparsity pattern contains elements that have been * added through other means (implicitly) while building it. For the current - * class, the result is true iff optimize_diag in the constructor or - * reinit() calls has been set to false, or if the represented matrix is not - * square. + * class, the result is true if and only if it is square because it then + * unconditionally stores the diagonal entries whether they have been added explicitly or not. * * This function mainly serves the purpose of describing the current class * in cases where several kinds of sparsity patterns can be passed as @@ -1128,7 +1348,7 @@ private: /** * Is special treatment of diagonals enabled? */ - bool diagonal_optimized; + bool store_diagonal_first_in_row; /** * Make all sparse matrices friends of this class. @@ -1442,7 +1662,7 @@ inline bool SparsityPattern::optimize_diagonal () const { - return diagonal_optimized; + return store_diagonal_first_in_row; } @@ -1450,12 +1670,7 @@ inline bool SparsityPattern::stores_only_added_elements () const { - if ((diagonal_optimized == true) - && - (n_cols() == n_rows())) - return false; - else - return true; + return store_diagonal_first_in_row; } @@ -1517,7 +1732,7 @@ SparsityPattern::save (Archive &ar, const unsigned int) const // forward to serialization function in the base class. ar &static_cast(*this); - ar &max_dim &rows &cols &max_vec_len &max_row_length &compressed &diagonal_optimized; + ar &max_dim &rows &cols &max_vec_len &max_row_length &compressed &store_diagonal_first_in_row; ar &boost::serialization::make_array(rowstart, max_dim + 1); ar &boost::serialization::make_array(colnums, max_vec_len); @@ -1533,7 +1748,7 @@ SparsityPattern::load (Archive &ar, const unsigned int) // forward to serialization function in the base class. ar &static_cast(*this); - ar &max_dim &rows &cols &max_vec_len &max_row_length &compressed &diagonal_optimized; + ar &max_dim &rows &cols &max_vec_len &max_row_length &compressed &store_diagonal_first_in_row; rowstart = new std::size_t [max_dim + 1]; colnums = new unsigned int [max_vec_len]; @@ -1555,7 +1770,7 @@ SparsityPattern::operator == (const SparsityPattern &sp2) const if (rows != sp2.rows || cols != sp2.cols || compressed != sp2.compressed || - diagonal_optimized != sp2.diagonal_optimized) + store_diagonal_first_in_row != sp2.store_diagonal_first_in_row) return false; for (unsigned int i = 0; i < rows+1; ++i) @@ -1612,7 +1827,19 @@ SparsityPattern::copy_from (const unsigned int n_rows, const unsigned int n_cols, const ForwardIterator begin, const ForwardIterator end, - const bool optimize_diag) + const bool) +{ + copy_from (n_rows, n_cols, begin, end); +} + + + +template +void +SparsityPattern::copy_from (const unsigned int n_rows, + const unsigned int n_cols, + const ForwardIterator begin, + const ForwardIterator end) { Assert (static_cast(std::distance (begin, end)) == n_rows, ExcIteratorRange (std::distance (begin, end), n_rows)); @@ -1621,14 +1848,14 @@ SparsityPattern::copy_from (const unsigned int n_rows, // 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 = optimize_diag && (n_rows == n_cols); + const bool is_square = (n_rows == n_cols); std::vector row_lengths; row_lengths.reserve(n_rows); for (ForwardIterator i=begin; i!=end; ++i) row_lengths.push_back (std::distance (i->begin(), i->end()) + (is_square ? 1 : 0)); - reinit (n_rows, n_cols, row_lengths, is_square); + reinit (n_rows, n_cols, row_lengths); // now enter all the elements into the matrix. note that if the matrix is // quadratic, then we already have the diagonal element preallocated @@ -1649,8 +1876,8 @@ SparsityPattern::copy_from (const unsigned int n_rows, if ((col!=row) || !is_square) *cols++ = col; - }; - }; + } + } // finally compress everything. this also sorts the entries within each row compress (); diff --git a/deal.II/source/lac/sparsity_pattern.cc b/deal.II/source/lac/sparsity_pattern.cc index c10384a8ad..3b2c76a383 100644 --- a/deal.II/source/lac/sparsity_pattern.cc +++ b/deal.II/source/lac/sparsity_pattern.cc @@ -43,7 +43,7 @@ SparsityPattern::SparsityPattern () rowstart(0), colnums(0), compressed(false), - diagonal_optimized(false) + store_diagonal_first_in_row(false) { reinit (0,0,0); } @@ -58,14 +58,14 @@ SparsityPattern::SparsityPattern (const SparsityPattern &s) rowstart(0), colnums(0), compressed(false), - diagonal_optimized(false) + store_diagonal_first_in_row(false) { Assert (s.rowstart == 0, ExcInvalidConstructorCall()); Assert (s.colnums == 0, ExcInvalidConstructorCall()); Assert (s.rows == 0, ExcInvalidConstructorCall()); Assert (s.cols == 0, ExcInvalidConstructorCall()); - reinit (0,0,0, false); + reinit (0,0,0); } @@ -73,32 +73,62 @@ SparsityPattern::SparsityPattern (const SparsityPattern &s) SparsityPattern::SparsityPattern (const unsigned int m, const unsigned int n, const unsigned int max_per_row, - const bool optimize_diag) + const bool) : max_dim(0), max_vec_len(0), rowstart(0), colnums(0), compressed(false), - diagonal_optimized(false) + store_diagonal_first_in_row(m == n) { - reinit (m,n,max_per_row, optimize_diag); + reinit (m,n,max_per_row); } -SparsityPattern::SparsityPattern ( - const unsigned int m, - const unsigned int n, - const std::vector &row_lengths, - const bool optimize_diag) +SparsityPattern::SparsityPattern (const unsigned int m, + const unsigned int n, + const unsigned int max_per_row) : max_dim(0), max_vec_len(0), rowstart(0), - colnums(0) + colnums(0), + compressed(false), + store_diagonal_first_in_row(m == n) { - reinit (m, n, row_lengths, optimize_diag); + reinit (m,n,max_per_row); +} + + + +SparsityPattern::SparsityPattern (const unsigned int m, + const unsigned int n, + const std::vector &row_lengths, + const bool) + : + max_dim(0), + max_vec_len(0), + rowstart(0), + colnums(0), + store_diagonal_first_in_row(m == n) +{ + reinit (m, n, row_lengths); +} + + +SparsityPattern::SparsityPattern (const unsigned int m, + const unsigned int n, + const std::vector &row_lengths) + : + max_dim(0), + max_vec_len(0), + rowstart(0), + colnums(0), + store_diagonal_first_in_row(m == n) +{ + reinit (m, n, row_lengths); } @@ -111,22 +141,33 @@ SparsityPattern::SparsityPattern (const unsigned int n, rowstart(0), colnums(0) { - reinit (n, n, max_per_row, true); + reinit (n, n, max_per_row); } -SparsityPattern::SparsityPattern ( - const unsigned int m, - const std::vector &row_lengths, - const bool optimize_diag) +SparsityPattern::SparsityPattern (const unsigned int m, + const std::vector &row_lengths, + const bool) : max_dim(0), max_vec_len(0), rowstart(0), colnums(0) { - reinit (m, m, row_lengths, optimize_diag); + reinit (m, m, row_lengths); +} + + +SparsityPattern::SparsityPattern (const unsigned int m, + const std::vector &row_lengths) + : + max_dim(0), + max_vec_len(0), + rowstart(0), + colnums(0) +{ + reinit (m, m, row_lengths); } @@ -141,11 +182,9 @@ SparsityPattern::SparsityPattern (const SparsityPattern &original, colnums(0) { Assert (original.rows==original.cols, ExcNotQuadratic()); - Assert (original.optimize_diagonal(), ExcNotQuadratic()); Assert (original.is_compressed(), ExcNotCompressed()); - reinit (original.rows, original.cols, max_per_row, - original.optimize_diagonal()); + reinit (original.rows, original.cols, max_per_row); // now copy the entries from // the other object @@ -265,22 +304,44 @@ void SparsityPattern::reinit (const unsigned int m, const unsigned int n, const unsigned int max_per_row, - const bool optimize_diag) + const bool) { // simply map this function to the // other @p{reinit} function const std::vector row_lengths (m, max_per_row); - reinit (m, n, row_lengths, optimize_diag); + reinit (m, n, row_lengths); } void -SparsityPattern::reinit ( - const unsigned int m, - const unsigned int n, - const VectorSlice > &row_lengths, - const bool optimize_diag) +SparsityPattern::reinit (const unsigned int m, + const unsigned int n, + const unsigned int max_per_row) +{ + // simply map this function to the + // other @p{reinit} function + const std::vector row_lengths (m, max_per_row); + reinit (m, n, row_lengths); +} + + + +void +SparsityPattern::reinit (const unsigned int m, + const unsigned int n, + const VectorSlice > &row_lengths, + const bool) +{ + reinit (m, n, row_lengths); +} + + + +void +SparsityPattern::reinit (const unsigned int m, + const unsigned int n, + const VectorSlice > &row_lengths) { AssertDimension (row_lengths.size(), m); @@ -303,14 +364,13 @@ SparsityPattern::reinit ( } // first, if the matrix is - // quadratic and special treatment - // of diagonals is requested, we + // quadratic, we // will have to make sure that each // row has at least one entry for // the diagonal element. make this // more obvious by having a // variable which we can query - diagonal_optimized = (m == n ? true : false) && optimize_diag; + store_diagonal_first_in_row = (m == n); // find out how many entries we // need in the @p{colnums} array. if @@ -323,7 +383,7 @@ SparsityPattern::reinit ( // of columns std::size_t vec_len = 0; for (unsigned int i=0; i 1) - std::sort ((diagonal_optimized) + std::sort ((store_diagonal_first_in_row) ? tmp_entries.begin()+1 : tmp_entries.begin(), tmp_entries.begin()+row_length); @@ -500,7 +560,7 @@ SparsityPattern::compress () // the diagonal element // (i.e. with column // index==line number) - Assert ((!diagonal_optimized) || + Assert ((!store_diagonal_first_in_row) || (new_colnums[rowstart[line]] == line), ExcInternalError()); // assert that the first entry @@ -555,7 +615,16 @@ SparsityPattern::compress () template void SparsityPattern::copy_from (const CSP &csp, - const bool optimize_diag) + const bool) +{ + copy_from (csp); +} + + + +template +void +SparsityPattern::copy_from (const CSP &csp) { // first determine row lengths for // each row. if the matrix is @@ -567,7 +636,7 @@ SparsityPattern::copy_from (const CSP &csp, // bother to check whether that // diagonal entry is in a certain // row or not - const bool do_diag_optimize = optimize_diag && (csp.n_rows() == csp.n_cols()); + const bool do_diag_optimize = (csp.n_rows() == csp.n_cols()); std::vector row_lengths (csp.n_rows()); for (unsigned int i=0; i void SparsityPattern::copy_from (const FullMatrix &matrix, const bool optimize_diag) +{ + copy_from (matrix); +} + + + +template +void SparsityPattern::copy_from (const FullMatrix &matrix) { // first init with the number of entries - // per row. if optimize_diag is set then we + // per row. if this matrix is square then we // also have to allocate memory for the // diagonal entry, unless we have already - // counted it or the matrix isn't square + // counted it std::vector entries_per_row (matrix.m(), 0); for (unsigned int row=0; row &matrix, ++entries_per_row[row]; } - reinit (matrix.m(), matrix.n(), entries_per_row, optimize_diag); + reinit (matrix.m(), matrix.n(), entries_per_row); // now set entries for (unsigned int row=0; row &matrix, void -SparsityPattern::reinit ( - const unsigned int m, - const unsigned int n, - const std::vector &row_lengths, - const bool optimize_diag) +SparsityPattern::reinit (const unsigned int m, + const unsigned int n, + const std::vector &row_lengths, + const bool) { - reinit(m, n, make_slice(row_lengths), optimize_diag); + reinit(m, n, make_slice(row_lengths)); } +void +SparsityPattern::reinit (const unsigned int m, + const unsigned int n, + const std::vector &row_lengths) +{ + reinit(m, n, make_slice(row_lengths)); +} + + + + bool SparsityPattern::empty () const { @@ -722,7 +809,7 @@ SparsityPattern::operator () (const unsigned int i, // was requested, we can get the // diagonal element faster by this // query. - if (diagonal_optimized && (i==j)) + if (store_diagonal_first_in_row && (i==j)) return rowstart[i]; // all other entries are sorted, so @@ -737,7 +824,7 @@ SparsityPattern::operator () (const unsigned int i, // top of this function, so it may // not be called for noncompressed // structures. - const unsigned int *sorted_region_start = (diagonal_optimized ? + const unsigned int *sorted_region_start = (store_diagonal_first_in_row ? &colnums[rowstart[i]+1] : &colnums[rowstart[i]]); const unsigned int *const p @@ -981,7 +1068,7 @@ SparsityPattern::block_write (std::ostream &out) const << max_vec_len << ' ' << max_row_length << ' ' << compressed << ' ' - << diagonal_optimized << "]["; + << store_diagonal_first_in_row << "]["; // then write out real data out.write (reinterpret_cast(&rowstart[0]), reinterpret_cast(&rowstart[max_dim+1]) @@ -1013,7 +1100,7 @@ SparsityPattern::block_read (std::istream &in) >> max_vec_len >> max_row_length >> compressed - >> diagonal_optimized; + >> store_diagonal_first_in_row; in >> c; AssertThrow (c == ']', ExcIO()); @@ -1072,6 +1159,12 @@ template void SparsityPattern::copy_from (const template void SparsityPattern::copy_from (const FullMatrix &, bool); template void SparsityPattern::copy_from (const FullMatrix &, bool); +template void SparsityPattern::copy_from (const CompressedSparsityPattern &); +template void SparsityPattern::copy_from (const CompressedSetSparsityPattern &); +template void SparsityPattern::copy_from (const CompressedSimpleSparsityPattern &); +template void SparsityPattern::copy_from (const FullMatrix &); +template void SparsityPattern::copy_from (const FullMatrix &); + template void SparsityPattern::add_entries (const unsigned int, const unsigned int *, const unsigned int *,