From 3d989cd549294af2defc2c53ded691b161ba16ff Mon Sep 17 00:00:00 2001 From: bangerth Date: Sat, 9 Feb 2013 05:15:09 +0000 Subject: [PATCH] Push the now generally true assumption that sparsity pattern for square matrices *always* store the diagonal element first in each row, into the ChunkSparsityPattern as well. git-svn-id: https://svn.dealii.org/trunk@28293 0785d39b-7218-0410-832d-ea1e28bc413d --- .../include/deal.II/lac/chunk_sparse_matrix.h | 7 +- .../lac/chunk_sparse_matrix.templates.h | 35 +-- .../deal.II/lac/chunk_sparsity_pattern.h | 217 ++++++++++-------- deal.II/source/lac/chunk_sparsity_pattern.cc | 116 ++++++++-- 4 files changed, 236 insertions(+), 139 deletions(-) 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 59cdb14d5b..1179da45fa 100644 --- a/deal.II/include/deal.II/lac/chunk_sparse_matrix.h +++ b/deal.II/include/deal.II/lac/chunk_sparse_matrix.h @@ -609,8 +609,7 @@ public: * element in the ith * row. This function throws an * error if the matrix is not - * quadratic (see - * ChunkSparsityPattern::optimize_diagonal()). + * quadratic. * * This function is considerably * faster than the operator()(), @@ -1431,7 +1430,7 @@ inline number ChunkSparseMatrix::diag_element (const unsigned int i) const { Assert (cols != 0, ExcNotInitialized()); - Assert (cols->optimize_diagonal(), ExcNotQuadratic()); + Assert (m() == n(), ExcNotQuadratic()); Assert (i::diag_element (const unsigned int i) { Assert (cols != 0, ExcNotInitialized()); - Assert (cols->optimize_diagonal(), ExcNotQuadratic()); + Assert (m() == n(), ExcNotQuadratic()); Assert (i::precondition_Jacobi (Vector &dst, { Assert (cols != 0, ExcNotInitialized()); Assert (val != 0, ExcNotInitialized()); - Assert (cols->optimize_diagonal(), - typename ChunkSparsityPattern::ExcDiagonalNotOptimized()); + Assert (m() == n(), ExcMessage("This operation is only valid on square matrices.")); Assert (dst.size() == n(), ExcDimensionMismatch (dst.size(), n())); Assert (src.size() == n(), ExcDimensionMismatch (src.size(), n())); @@ -1202,8 +1201,7 @@ ChunkSparseMatrix::precondition_SSOR (Vector &dst, // which is much clearer... Assert (cols != 0, ExcNotInitialized()); Assert (val != 0, ExcNotInitialized()); - Assert (cols->optimize_diagonal(), - typename ChunkSparsityPattern::ExcDiagonalNotOptimized()); + Assert (m() == n(), ExcMessage("This operation is only valid on square matrices.")); Assert (dst.size() == n(), ExcDimensionMismatch (dst.size(), n())); Assert (src.size() == n(), ExcDimensionMismatch (src.size(), n())); @@ -1221,8 +1219,7 @@ ChunkSparseMatrix::precondition_SOR (Vector &dst, { Assert (cols != 0, ExcNotInitialized()); Assert (val != 0, ExcNotInitialized()); - Assert (cols->optimize_diagonal(), - typename ChunkSparsityPattern::ExcDiagonalNotOptimized()); + Assert (m() == n(), ExcMessage("This operation is only valid on square matrices.")); dst = src; @@ -1239,8 +1236,7 @@ ChunkSparseMatrix::precondition_TSOR (Vector &dst, { Assert (cols != 0, ExcNotInitialized()); Assert (val != 0, ExcNotInitialized()); - Assert (cols->optimize_diagonal(), - typename ChunkSparsityPattern::ExcDiagonalNotOptimized()); + Assert (m() == n(), ExcMessage("This operation is only valid on square matrices.")); dst = src; @@ -1256,9 +1252,7 @@ ChunkSparseMatrix::SOR (Vector &dst, { Assert (cols != 0, ExcNotInitialized()); Assert (val != 0, ExcNotInitialized()); - Assert (cols->optimize_diagonal(), - typename ChunkSparsityPattern::ExcDiagonalNotOptimized()); - + Assert (m() == n(), ExcMessage("This operation is only valid on square matrices.")); Assert (m() == dst.size(), ExcDimensionMismatch(m(),dst.size())); Assert (false, ExcNotImplemented()); @@ -1273,9 +1267,7 @@ ChunkSparseMatrix::TSOR (Vector &dst, { Assert (cols != 0, ExcNotInitialized()); Assert (val != 0, ExcNotInitialized()); - Assert (cols->optimize_diagonal(), - typename ChunkSparsityPattern::ExcDiagonalNotOptimized()); - + Assert (m() == n(), ExcMessage("This operation is only valid on square matrices.")); Assert (m() == dst.size(), ExcDimensionMismatch(m(),dst.size())); Assert (false, ExcNotImplemented()); @@ -1292,8 +1284,7 @@ ChunkSparseMatrix::PSOR (Vector &dst, { Assert (cols != 0, ExcNotInitialized()); Assert (val != 0, ExcNotInitialized()); - Assert (cols->optimize_diagonal(), - typename ChunkSparsityPattern::ExcDiagonalNotOptimized()); + Assert (m() == n(), ExcMessage("This operation is only valid on square matrices.")); Assert (m() == dst.size(), ExcDimensionMismatch(m(),dst.size())); Assert (m() == permutation.size(), @@ -1315,8 +1306,7 @@ ChunkSparseMatrix::TPSOR (Vector &dst, { Assert (cols != 0, ExcNotInitialized()); Assert (val != 0, ExcNotInitialized()); - Assert (cols->optimize_diagonal(), - typename ChunkSparsityPattern::ExcDiagonalNotOptimized()); + Assert (m() == n(), ExcMessage("This operation is only valid on square matrices.")); Assert (m() == dst.size(), ExcDimensionMismatch(m(),dst.size())); Assert (m() == permutation.size(), @@ -1338,8 +1328,7 @@ ChunkSparseMatrix::SOR_step (Vector &v, { Assert (cols != 0, ExcNotInitialized()); Assert (val != 0, ExcNotInitialized()); - Assert (cols->optimize_diagonal(), - typename ChunkSparsityPattern::ExcDiagonalNotOptimized()); + Assert (m() == n(), ExcMessage("This operation is only valid on square matrices.")); Assert (m() == v.size(), ExcDimensionMismatch(m(),v.size())); Assert (m() == b.size(), ExcDimensionMismatch(m(),b.size())); @@ -1358,8 +1347,7 @@ ChunkSparseMatrix::TSOR_step (Vector &v, { Assert (cols != 0, ExcNotInitialized()); Assert (val != 0, ExcNotInitialized()); - Assert (cols->optimize_diagonal(), - typename ChunkSparsityPattern::ExcDiagonalNotOptimized()); + Assert (m() == n(), ExcMessage("This operation is only valid on square matrices.")); Assert (m() == v.size(), ExcDimensionMismatch(m(),v.size())); Assert (m() == b.size(), ExcDimensionMismatch(m(),b.size())); @@ -1390,8 +1378,7 @@ ChunkSparseMatrix::SSOR (Vector &dst, { Assert (cols != 0, ExcNotInitialized()); Assert (val != 0, ExcNotInitialized()); - Assert (cols->optimize_diagonal(), - typename ChunkSparsityPattern::ExcDiagonalNotOptimized()); + Assert (m() == n(), ExcMessage("This operation is only valid on square matrices.")); Assert (m() == dst.size(), ExcDimensionMismatch(m(),dst.size())); 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 24115a9238..8f7dc5965b 100644 --- a/deal.II/include/deal.II/lac/chunk_sparsity_pattern.h +++ b/deal.II/include/deal.II/lac/chunk_sparsity_pattern.h @@ -126,18 +126,21 @@ public: * @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. + */ + ChunkSparsityPattern (const unsigned int m, + const unsigned int n, + const unsigned int max_chunks_per_row, + const unsigned int chunk_size); + + /** + * @deprecated This constructor is deprecated. Use the version + * without the last argument */ ChunkSparsityPattern (const unsigned int m, const unsigned int n, const unsigned int max_chunks_per_row, const unsigned int chunk_size, - const bool optimize_diagonal = true); + const bool optimize_diagonal) DEAL_II_DEPRECATED; /** * Initialize a rectangular @@ -150,18 +153,21 @@ public: * 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. + */ + ChunkSparsityPattern (const unsigned int m, + const unsigned int n, + const std::vector &row_lengths, + const unsigned int chunk_size); + + /** + * @deprecated This constructor is deprecated. Use the version + * without the last argument */ ChunkSparsityPattern (const unsigned int m, const unsigned int n, const std::vector &row_lengths, const unsigned int chunk_size, - const bool optimize_diagonal = true); + const bool optimize_diagonal) DEAL_II_DEPRECATED; /** * Initialize a quadratic matrix @@ -189,15 +195,19 @@ public: * 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(). + */ + ChunkSparsityPattern (const unsigned int m, + const std::vector &row_lengths, + const unsigned int chunk_size); + + /** + * @deprecated This constructor is deprecated. Use the version + * without the last argument */ ChunkSparsityPattern (const unsigned int m, const std::vector &row_lengths, const unsigned int chunk_size, - const bool optimize_diagonal = true); + const bool optimize_diagonal) DEAL_II_DEPRECATED; /** * Destructor. @@ -225,11 +235,20 @@ public: * operations to the other * reinit function. */ + void reinit (const unsigned int m, + const unsigned int n, + const unsigned int max_per_row, + const unsigned int chunk_size); + + /** + * @deprecated This function is deprecated. Use the function + * without the last argument + */ void reinit (const unsigned int m, const unsigned int n, const unsigned int max_per_row, const unsigned int chunk_size, - const bool optimize_diagonal = true); + const bool optimize_diagonal) DEAL_II_DEPRECATED; /** * Reallocate memory for a matrix @@ -249,28 +268,45 @@ public: * of the heap. * * If the number of rows equals - * the number of columns and the - * last parameter is true, + * the number of columns then * 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, + const unsigned int chunk_size); + + /** + * @deprecated This function is deprecated. Use the function + * without the last argument + */ void reinit (const unsigned int m, const unsigned int n, const std::vector &row_lengths, const unsigned int chunk_size, - 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, + const unsigned int chunk_size); + + /** + * @deprecated This function is deprecated. Use the function + * without the last argument + */ void reinit (const unsigned int m, const unsigned int n, const VectorSlice > &row_lengths, const unsigned int chunk_size, - const bool optimize_diagonal = true); + const bool optimize_diagonal) DEAL_II_DEPRECATED; /** * This function compresses the sparsity @@ -433,12 +469,23 @@ public: * pairs. */ template + void copy_from (const unsigned int n_rows, + const unsigned int n_cols, + const ForwardIterator begin, + const ForwardIterator end, + const unsigned int chunk_size); + + /** + * @deprecated This function is deprecated. Use the 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 unsigned int chunk_size, - const bool optimize_diagonal = true); + const bool optimize_diagonal) DEAL_II_DEPRECATED; /** * Copy data from an object of type @@ -450,9 +497,17 @@ public: * compressed mode afterwards. */ template + void copy_from (const SparsityType &csp, + const unsigned int chunk_size); + + /** + * @deprecated This function is deprecated. Use the function + * without the last argument + */ + template void copy_from (const SparsityType &csp, const unsigned int chunk_size, - const bool optimize_diagonal = true); + const bool optimize_diagonal) DEAL_II_DEPRECATED; /** * Take a full matrix and use its @@ -466,9 +521,17 @@ public: * compressed mode afterwards. */ template + void copy_from (const FullMatrix &matrix, + const unsigned int chunk_size); + + /** + * @deprecated This function is deprecated. Use the function + * without the last argument + */ + template void copy_from (const FullMatrix &matrix, const unsigned int chunk_size, - const bool optimize_diagonal = true); + const bool optimize_diagonal) DEAL_II_DEPRECATED; /** * Return whether the object is empty. It @@ -578,62 +641,23 @@ public: bool is_compressed () 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). - * - * 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. - */ - bool optimize_diagonal () const; - - /** - * 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 false because we store - * entire chunks, not individual - * elements, and adding one entry to the - * sparsity pattern requires also adding - * all the other elements of a chunk. The - * only exception is if - * chunk_size==1, the - * sparsity pattern is nonsymmetric or - * optimize_diag has been set to false. - * - * This function mainly serves the - * purpose of describing the current - * class in cases where several kinds of - * sparsity patterns can be passed as + * Determine whether the matrix uses the special convention for quadratic + * matrices that the diagonal entries are stored first in each row. + * + * @deprecated The function always returns true if the matrix is + * square and false if it is not. + */ + 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 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 * template arguments. */ bool stores_only_added_elements () const; @@ -770,16 +794,6 @@ public: * Exception */ DeclException0 (ExcInvalidConstructorCall); - /** - * This exception is thrown if - * the matrix does not follow the - * convention of storing diagonal - * elements first in row. Refer - * to - * SparityPattern::optimize_diagonal() - * for more information. - */ - DeclException0 (ExcDiagonalNotOptimized); /** * Exception */ @@ -896,7 +910,20 @@ ChunkSparsityPattern::copy_from (const unsigned int n_rows, const ForwardIterator begin, const ForwardIterator end, const unsigned int chunk_size, - const bool optimize_diag) + const bool) +{ + copy_from (n_rows, n_cols, begin, end, chunk_size); +} + + + +template +void +ChunkSparsityPattern::copy_from (const unsigned int n_rows, + const unsigned int n_cols, + const ForwardIterator begin, + const ForwardIterator end, + const unsigned int chunk_size) { Assert (static_cast(std::distance (begin, end)) == n_rows, ExcIteratorRange (std::distance (begin, end), n_rows)); @@ -911,7 +938,7 @@ ChunkSparsityPattern::copy_from (const unsigned int n_rows, // 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) diff --git a/deal.II/source/lac/chunk_sparsity_pattern.cc b/deal.II/source/lac/chunk_sparsity_pattern.cc index 5c15733a95..2983d2b747 100644 --- a/deal.II/source/lac/chunk_sparsity_pattern.cc +++ b/deal.II/source/lac/chunk_sparsity_pattern.cc @@ -37,7 +37,7 @@ ChunkSparsityPattern::ChunkSparsityPattern (const ChunkSparsityPattern &s) Assert (s.rows == 0, ExcInvalidConstructorCall()); Assert (s.cols == 0, ExcInvalidConstructorCall()); - reinit (0,0,0,0, false); + reinit (0,0,0,0); } @@ -46,11 +46,22 @@ ChunkSparsityPattern::ChunkSparsityPattern (const unsigned int m, const unsigned int n, const unsigned int max_per_row, const unsigned int chunk_size, - const bool optimize_diag) + const bool) { Assert (chunk_size > 0, ExcInvalidNumber (chunk_size)); - reinit (m,n,max_per_row, chunk_size, optimize_diag); + reinit (m,n,max_per_row, chunk_size); +} + + +ChunkSparsityPattern::ChunkSparsityPattern (const unsigned int m, + const unsigned int n, + const unsigned int max_per_row, + const unsigned int chunk_size) +{ + Assert (chunk_size > 0, ExcInvalidNumber (chunk_size)); + + reinit (m,n,max_per_row, chunk_size); } @@ -64,7 +75,20 @@ ChunkSparsityPattern::ChunkSparsityPattern ( { Assert (chunk_size > 0, ExcInvalidNumber (chunk_size)); - reinit (m, n, row_lengths, chunk_size, optimize_diag); + reinit (m, n, row_lengths, chunk_size); +} + + + +ChunkSparsityPattern::ChunkSparsityPattern ( + const unsigned int m, + const unsigned int n, + const std::vector &row_lengths, + const unsigned int chunk_size) +{ + Assert (chunk_size > 0, ExcInvalidNumber (chunk_size)); + + reinit (m, n, row_lengths, chunk_size); } @@ -73,7 +97,7 @@ ChunkSparsityPattern::ChunkSparsityPattern (const unsigned int n, const unsigned int max_per_row, const unsigned int chunk_size) { - reinit (n, n, max_per_row, chunk_size, true); + reinit (n, n, max_per_row, chunk_size); } @@ -86,7 +110,19 @@ ChunkSparsityPattern::ChunkSparsityPattern ( { Assert (chunk_size > 0, ExcInvalidNumber (chunk_size)); - reinit (m, m, row_lengths, chunk_size, optimize_diag); + reinit (m, m, row_lengths, chunk_size); +} + + + +ChunkSparsityPattern::ChunkSparsityPattern ( + const unsigned int m, + const std::vector &row_lengths, + const unsigned int chunk_size) +{ + Assert (chunk_size > 0, ExcInvalidNumber (chunk_size)); + + reinit (m, m, row_lengths, chunk_size); } @@ -116,14 +152,24 @@ ChunkSparsityPattern::reinit (const unsigned int m, const unsigned int n, const unsigned int max_per_row, const unsigned int chunk_size, - const bool optimize_diag) + const bool) +{ + reinit (m, n, max_per_row, chunk_size); +} + + +void +ChunkSparsityPattern::reinit (const unsigned int m, + const unsigned int n, + const unsigned int max_per_row, + const unsigned int chunk_size) { Assert (chunk_size > 0, ExcInvalidNumber (chunk_size)); // 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, optimize_diag); + reinit (m, n, row_lengths, chunk_size); } @@ -134,7 +180,18 @@ ChunkSparsityPattern::reinit ( const unsigned int n, const VectorSlice > &row_lengths, const unsigned int chunk_size, - const bool optimize_diag) + const bool) +{ + reinit (m, n, row_lengths, chunk_size); +} + + +void +ChunkSparsityPattern::reinit ( + const unsigned int m, + const unsigned int n, + const VectorSlice > &row_lengths, + const unsigned int chunk_size) { Assert (row_lengths.size() == m, ExcInvalidNumber (m)); Assert (chunk_size > 0, ExcInvalidNumber (chunk_size)); @@ -171,8 +228,7 @@ ChunkSparsityPattern::reinit ( sparsity_pattern.reinit (m_chunks, n_chunks, - chunk_row_lengths, - optimize_diag); + chunk_row_lengths); } @@ -189,7 +245,16 @@ template void ChunkSparsityPattern::copy_from (const SparsityType &csp, const unsigned int chunk_size, - const bool optimize_diag) + const bool) +{ + copy_from (csp, chunk_size); +} + + +template +void +ChunkSparsityPattern::copy_from (const SparsityType &csp, + const unsigned int chunk_size) { Assert (chunk_size > 0, ExcInvalidNumber (chunk_size)); @@ -202,7 +267,7 @@ ChunkSparsityPattern::copy_from (const SparsityType &csp, reinit (csp.n_rows(), csp.n_cols(), entries_per_row, - chunk_size, optimize_diag); + chunk_size); // then actually fill it for (unsigned int row = 0; row void ChunkSparsityPattern::copy_from (const FullMatrix &matrix, const unsigned int chunk_size, const bool optimize_diag) +{ + copy_from (matrix, chunk_size); +} + + +template +void ChunkSparsityPattern::copy_from (const FullMatrix &matrix, + const unsigned int chunk_size) { Assert (chunk_size > 0, ExcInvalidNumber (chunk_size)); @@ -238,7 +311,7 @@ void ChunkSparsityPattern::copy_from (const FullMatrix &matrix, reinit (matrix.m(), matrix.n(), entries_per_row, - chunk_size, optimize_diag); + chunk_size); // then actually fill it for (unsigned int row=0; row &row_lengths, const unsigned int chunk_size, - const bool optimize_diag) + const bool) +{ + reinit (m, n, row_lengths, chunk_size); +} + + +void +ChunkSparsityPattern::reinit ( + const unsigned int m, + const unsigned int n, + const std::vector &row_lengths, + const unsigned int chunk_size) { Assert (chunk_size > 0, ExcInvalidNumber (chunk_size)); - reinit(m, n, make_slice(row_lengths), chunk_size, optimize_diag); + reinit(m, n, make_slice(row_lengths), chunk_size); } -- 2.39.5