From e861bfce1f60e6af47ea0889c2f84ab55bbf61d5 Mon Sep 17 00:00:00 2001 From: Denis Davydov Date: Tue, 9 Oct 2018 09:59:46 +0200 Subject: [PATCH] add DynamicSparsityPattern::compute_Tmmult_pattern() by moving existing code --- .../deal.II/lac/dynamic_sparsity_pattern.h | 10 ++++- include/deal.II/lac/sparse_matrix.templates.h | 35 ++-------------- include/deal.II/lac/sparsity_pattern.h | 1 + source/lac/dynamic_sparsity_pattern.cc | 42 +++++++++++++++++++ 4 files changed, 55 insertions(+), 33 deletions(-) diff --git a/include/deal.II/lac/dynamic_sparsity_pattern.h b/include/deal.II/lac/dynamic_sparsity_pattern.h index e188579eac..c33fe8edcb 100644 --- a/include/deal.II/lac/dynamic_sparsity_pattern.h +++ b/include/deal.II/lac/dynamic_sparsity_pattern.h @@ -32,7 +32,7 @@ DEAL_II_NAMESPACE_OPEN class DynamicSparsityPattern; - +class SparsityPattern; /*! @addtogroup Sparsity *@{ @@ -446,6 +446,14 @@ public: compute_mmult_pattern(const SparsityPatternTypeLeft & left, const SparsityPatternTypeRight &right); + /** + * Construct and store in this object the sparsity pattern corresponding to + * the product of transposed @p left and and non-transpose @p right sparsity pattern. + */ + void + compute_Tmmult_pattern(const SparsityPattern &left, + const SparsityPattern &right); + /** * Print the sparsity pattern. The output consists of one line per row of * the format [i,j1,j2,j3,...]. i is the row number and diff --git a/include/deal.II/lac/sparse_matrix.templates.h b/include/deal.II/lac/sparse_matrix.templates.h index 2566bb0cb9..b7414ff150 100644 --- a/include/deal.II/lac/sparse_matrix.templates.h +++ b/include/deal.II/lac/sparse_matrix.templates.h @@ -1145,39 +1145,10 @@ SparseMatrix::Tmmult(SparseMatrix & C, C.clear(); sp_C.reinit(0, 0, 0); - // create a sparsity pattern for the matrix. we will go through all the - // rows in the matrix A, and for each column in a row we add the whole - // row of matrix B with that row number. This means that we will insert - // a lot of entries to each row, which is best handled by the - // DynamicSparsityPattern class. + // create a sparsity pattern for the matrix. { - DynamicSparsityPattern dsp(n(), B.n()); - for (size_type i = 0; i < sp_A.n_rows(); ++i) - { - const size_type * rows = &sp_A.colnums[sp_A.rowstart[i]]; - const size_type *const end_rows = - &sp_A.colnums[sp_A.rowstart[i + 1]]; - // cast away constness to conform with dsp.add_entries interface - size_type *new_cols = - const_cast(&sp_B.colnums[sp_B.rowstart[i]]); - size_type *end_new_cols = - const_cast(&sp_B.colnums[sp_B.rowstart[i + 1]]); - - if (sp_B.n_rows() == sp_B.n_cols()) - ++new_cols; - - for (; rows != end_rows; ++rows) - { - const size_type row = *rows; - - // if B has a diagonal, need to add that manually. this way, - // we maintain sortedness. - if (sp_B.n_rows() == sp_B.n_cols()) - dsp.add(row, i); - - dsp.add_entries(row, new_cols, end_new_cols, true); - } - } + DynamicSparsityPattern dsp; + dsp.compute_Tmmult_pattern(sp_A, sp_B); sp_C.copy_from(dsp); } diff --git a/include/deal.II/lac/sparsity_pattern.h b/include/deal.II/lac/sparsity_pattern.h index 37a1be647b..3caaa9d4a7 100644 --- a/include/deal.II/lac/sparsity_pattern.h +++ b/include/deal.II/lac/sparsity_pattern.h @@ -1129,6 +1129,7 @@ private: friend class ChunkSparseMatrix; friend class ChunkSparsityPattern; + friend class DynamicSparsityPattern; /** * Also give access to internal details to the iterator/accessor classes. diff --git a/source/lac/dynamic_sparsity_pattern.cc b/source/lac/dynamic_sparsity_pattern.cc index 567340242a..62ca57ff73 100644 --- a/source/lac/dynamic_sparsity_pattern.cc +++ b/source/lac/dynamic_sparsity_pattern.cc @@ -405,6 +405,48 @@ DynamicSparsityPattern::symmetrize() +void +DynamicSparsityPattern::compute_Tmmult_pattern(const SparsityPattern &sp_A, + const SparsityPattern &sp_B) +{ + Assert(sp_A.n_rows() == sp_B.n_rows(), + ExcDimensionMismatch(sp_A.n_rows(), sp_B.n_rows())); + + this->reinit(sp_A.n_cols(), sp_B.n_cols()); + // we will go through all the + // rows in the matrix A, and for each column in a row we add the whole + // row of matrix B with that row number. This means that we will insert + // a lot of entries to each row, which is best handled by the + // DynamicSparsityPattern class. + for (size_type i = 0; i < sp_A.n_rows(); ++i) + { + const size_type * rows = &sp_A.colnums[sp_A.rowstart[i]]; + const size_type *const end_rows = &sp_A.colnums[sp_A.rowstart[i + 1]]; + // cast away constness to conform with dsp.add_entries interface + size_type *new_cols = + const_cast(&sp_B.colnums[sp_B.rowstart[i]]); + size_type *end_new_cols = + const_cast(&sp_B.colnums[sp_B.rowstart[i + 1]]); + + if (sp_B.n_rows() == sp_B.n_cols()) + ++new_cols; + + for (; rows != end_rows; ++rows) + { + const size_type row = *rows; + + // if B has a diagonal, need to add that manually. this way, + // we maintain sortedness. + if (sp_B.n_rows() == sp_B.n_cols()) + this->add(row, i); + + this->add_entries(row, new_cols, end_new_cols, true); + } + } +} + + + template void DynamicSparsityPattern::compute_mmult_pattern( -- 2.39.5