From ae7e7f5c2966673ca6df8d1d75cf10e552d030e7 Mon Sep 17 00:00:00 2001 From: Luca Heltai Date: Tue, 10 Jan 2006 00:34:22 +0000 Subject: [PATCH] New Version of Filtered Matrix, which requires only vector matrix mult. git-svn-id: https://svn.dealii.org/trunk@11979 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/lac/include/lac/filtered_matrix.h | 34 +-- .../include/lac/filtered_matrix.templates.h | 258 ++++++++---------- deal.II/lac/source/filtered_matrix.cc | 107 +------- 3 files changed, 126 insertions(+), 273 deletions(-) diff --git a/deal.II/lac/include/lac/filtered_matrix.h b/deal.II/lac/include/lac/filtered_matrix.h index 85085f8699..52464ec186 100644 --- a/deal.II/lac/include/lac/filtered_matrix.h +++ b/deal.II/lac/include/lac/filtered_matrix.h @@ -169,7 +169,7 @@ template class Vector; *

Template arguments

* * This class takes as template arguments a matrix and a vector - * class. The former must provide @p vmult, @p Tvmult, and + * class. The former must provide @p vmult, @p vmult_add, @p Tvmult, and * @p residual member function that operate on the vector type (the * second template argument). The latter template parameter must * provide access to indivual elements through operator(), @@ -190,7 +190,7 @@ template class Vector; * bottleneck. If you don't want this serialization of operations, you * have to use several objects of this type. * - * @author Wolfgang Bangerth 2001 + * @author Wolfgang Bangerth 2001, Luca Heltai 2006 */ template > class FilteredMatrix : public Subscriptor @@ -278,19 +278,7 @@ class FilteredMatrix : public Subscriptor * contains an entry for a degree * of freedom that has already * been constrained - * previously. Furthermore, it is - * assumed that the list of - * constraints is sorted. If the - * input argument is a - * std::map, - * this is automatically the - * case. - * - * Note that we do not check - * whether the input range is - * sorted, as this would be too - * expensive. You have to ensure - * this yourself. + * previously. */ template void add_constraints (const ConstraintList &new_constraints); @@ -552,22 +540,6 @@ class FilteredMatrix : public Subscriptor */ void allocate_tmp_vector (); - /** - * Determine all entries in the - * given column of the matrix - * except for the diagonal entry - * and return their index/value - * pairs. If the matrix is - * symmetric, use a faster - * algorithm. - * - * This function needs to be - * specialized for the different - * matrix types. - */ - void get_column_entries (const unsigned int index, - std::vector &column_entries, - const bool matrix_is_symmetric) const; }; /*@}*/ diff --git a/deal.II/lac/include/lac/filtered_matrix.templates.h b/deal.II/lac/include/lac/filtered_matrix.templates.h index f5947de0b2..2a3131f22e 100644 --- a/deal.II/lac/include/lac/filtered_matrix.templates.h +++ b/deal.II/lac/include/lac/filtered_matrix.templates.h @@ -22,13 +22,10 @@ #include #include - template FilteredMatrix::FilteredMatrix () {} - - template FilteredMatrix:: FilteredMatrix (const FilteredMatrix &fm) @@ -67,7 +64,7 @@ FilteredMatrix:: set_referenced_matrix (const MATRIX &m) { matrix = &m; - allocate_tmp_vector (); + allocate_tmp_vector(); } @@ -87,55 +84,20 @@ template void FilteredMatrix:: apply_constraints (VECTOR &v, - const bool matrix_is_symmetric) const + const bool /* matrix_is_symmetric */) const { - // array that will hold the pairs - // of index/value of all nonzero - // entries in a given column - std::vector column_entries; - - // iterate over all constraints and - // treat them one after the other - const_index_value_iterator i = constraints.begin(); - const const_index_value_iterator e = constraints.end(); - for (; i!=e; ++i) - { - // define abbreviations - const unsigned index = i->first; - const value_type value = i->second; - - // check whether the value is - // zero, since in that case we do - // not have to modify other nodes - if (value != 0) - { - // first clear array of - // previous content - column_entries.clear (); - - // then get all entries in - // the present column - get_column_entries (index, column_entries, matrix_is_symmetric); - - // modify rhs for each entry - const_index_value_iterator col = column_entries.begin(); - const const_index_value_iterator col_end = column_entries.end(); - for (; col!=col_end; ++col) - v(col->first) -= col->second * value; - }; - }; - - - // finally set constrained - // entries themselves. we can't - // do it in the above loop - // since we might end up - // modifying an entry that we - // have already set if - // constrained dofs couple to - // each other - for (i=constraints.begin(); i!=e; ++i) - v(i->first) = i->second; + tmp_vector = 0; + const_index_value_iterator i = constraints.begin(); + const const_index_value_iterator e = constraints.end(); + for (; i!=e; ++i) + tmp_vector(i->first) = -i->second; + + // This vmult is without bc, to get the rhs correction in a correct way. + matrix->vmult_add(v, tmp_vector); + + // finally set constrained entries themselves + for (i=constraints.begin(); i!=e; ++i) + v(i->first) = i->second; } @@ -144,12 +106,12 @@ template void FilteredMatrix::pre_filter (VECTOR &v) const { - // iterate over all constraints and - // zero out value - const_index_value_iterator i = constraints.begin(); - const const_index_value_iterator e = constraints.end(); - for (; i!=e; ++i) - v(i->first) = 0; + // iterate over all constraints and + // zero out value + const_index_value_iterator i = constraints.begin(); + const const_index_value_iterator e = constraints.end(); + for (; i!=e; ++i) + v(i->first) = 0; } @@ -157,14 +119,14 @@ FilteredMatrix::pre_filter (VECTOR &v) const template void FilteredMatrix::post_filter (const VECTOR &in, - VECTOR &out) const + VECTOR &out) const { - // iterate over all constraints and - // set value correctly - const_index_value_iterator i = constraints.begin(); - const const_index_value_iterator e = constraints.end(); - for (; i!=e; ++i) - out(i->first) = in(i->first); + // iterate over all constraints and + // set value correctly + const_index_value_iterator i = constraints.begin(); + const const_index_value_iterator e = constraints.end(); + for (; i!=e; ++i) + out(i->first) = in(i->first); } @@ -172,19 +134,19 @@ FilteredMatrix::post_filter (const VECTOR &in, template void FilteredMatrix::vmult (VECTOR &dst, - const VECTOR &src) const + const VECTOR &src) const { - tmp_mutex.acquire (); - // first copy over src vector and - // pre-filter - tmp_vector = src; - pre_filter (tmp_vector); - // then let matrix do its work - matrix->vmult (dst, tmp_vector); - // tmp_vector now no more needed - tmp_mutex.release (); - // finally do post-filtering - post_filter (src, dst); + tmp_mutex.acquire (); + // first copy over src vector and + // pre-filter + tmp_vector = src; + pre_filter (tmp_vector); + // then let matrix do its work + matrix->vmult (dst, tmp_vector); + // tmp_vector now no more needed + tmp_mutex.release (); + // finally do post-filtering + post_filter (src, dst); } @@ -192,34 +154,34 @@ FilteredMatrix::vmult (VECTOR &dst, template typename FilteredMatrix::value_type FilteredMatrix::residual (VECTOR &dst, - const VECTOR &x, - const VECTOR &b) const + const VECTOR &x, + const VECTOR &b) const { - tmp_mutex.acquire (); - // first copy over x vector and - // pre-filter - tmp_vector = x; - pre_filter (tmp_vector); - // then let matrix do its work - value_type res = matrix->residual (dst, tmp_vector, b); - value_type res2 = res*res; - // tmp_vector now no more needed - tmp_mutex.release (); - // finally do post-filtering. here, - // we set constrained indices to - // zero, but have to subtract their - // contributions to the residual - const_index_value_iterator i = constraints.begin(); - const const_index_value_iterator e = constraints.end(); - for (; i!=e; ++i) + tmp_mutex.acquire (); + // first copy over x vector and + // pre-filter + tmp_vector = x; + pre_filter (tmp_vector); + // then let matrix do its work + value_type res = matrix->residual (dst, tmp_vector, b); + value_type res2 = res*res; + // tmp_vector now no more needed + tmp_mutex.release (); + // finally do post-filtering. here, + // we set constrained indices to + // zero, but have to subtract their + // contributions to the residual + const_index_value_iterator i = constraints.begin(); + const const_index_value_iterator e = constraints.end(); + for (; i!=e; ++i) { - const value_type v = dst(i->first); - res2 -= v*v; - dst(i->first) = 0; + const value_type v = dst(i->first); + res2 -= v*v; + dst(i->first) = 0; }; - - Assert (res2>=0, ExcInternalError()); - return std::sqrt (res2); + + Assert (res2>=0, ExcInternalError()); + return std::sqrt (res2); } @@ -227,41 +189,41 @@ FilteredMatrix::residual (VECTOR &dst, template void FilteredMatrix::Tvmult (VECTOR &dst, - const VECTOR &src) const + const VECTOR &src) const { - tmp_mutex.acquire (); - // first copy over src vector and - // pre-filter - tmp_vector = src; - pre_filter (tmp_vector); - // then let matrix do its work - matrix->Tvmult (dst, tmp_vector); - // tmp_vector now no more needed - tmp_mutex.release (); - // finally do post-filtering - post_filter (src, dst); + tmp_mutex.acquire (); + // first copy over src vector and + // pre-filter + tmp_vector = src; + pre_filter (tmp_vector); + // then let matrix do its work + matrix->Tvmult (dst, tmp_vector); + // tmp_vector now no more needed + tmp_mutex.release (); + // finally do post-filtering + post_filter (src, dst); } - + template typename FilteredMatrix::value_type FilteredMatrix::matrix_norm_square (const VECTOR &v) const { - tmp_mutex.acquire (); - tmp_vector = v; - - // zero out constrained entries and - // form matrix norm with original - // matrix. this is equivalent to - // forming the matrix norm of the - // original vector with the matrix - // where we have zeroed out rows - // and columns - pre_filter (tmp_vector); - const value_type ret = matrix->matrix_norm_square (tmp_vector); - tmp_mutex.release (); - return ret; + tmp_mutex.acquire (); + tmp_vector = v; + + // zero out constrained entries and + // form matrix norm with original + // matrix. this is equivalent to + // forming the matrix norm of the + // original vector with the matrix + // where we have zeroed out rows + // and columns + pre_filter (tmp_vector); + const value_type ret = matrix->matrix_norm_square (tmp_vector); + tmp_mutex.release (); + return ret; } @@ -270,23 +232,23 @@ template void FilteredMatrix:: precondition_Jacobi (VECTOR &dst, - const VECTOR &src, - const value_type omega) const + const VECTOR &src, + const value_type omega) const { - // first precondition as usual, - // using the fast algorithms of the - // matrix class - matrix->precondition_Jacobi (dst, src, omega); - - // then modify the constrained - // degree of freedom. as the - // diagonal entries of the filtered - // matrix would be 1.0, simply copy - // over old and new values - const_index_value_iterator i = constraints.begin(); - const const_index_value_iterator e = constraints.end(); - for (; i!=e; ++i) - dst(i->first) = src(i->first); + // first precondition as usual, + // using the fast algorithms of the + // matrix class + matrix->precondition_Jacobi (dst, src, omega); + + // then modify the constrained + // degree of freedom. as the + // diagonal entries of the filtered + // matrix would be 1.0, simply copy + // over old and new values + const_index_value_iterator i = constraints.begin(); + const const_index_value_iterator e = constraints.end(); + for (; i!=e; ++i) + dst(i->first) = src(i->first); } @@ -295,9 +257,9 @@ template unsigned int FilteredMatrix::memory_consumption () const { - return (MemoryConsumption::memory_consumption (matrix) + - MemoryConsumption::memory_consumption (constraints) + - MemoryConsumption::memory_consumption (tmp_vector)); + return (MemoryConsumption::memory_consumption (matrix) + + MemoryConsumption::memory_consumption (constraints) + + MemoryConsumption::memory_consumption (tmp_vector)); } diff --git a/deal.II/lac/source/filtered_matrix.cc b/deal.II/lac/source/filtered_matrix.cc index 3bd0d036f6..bc6421cde0 100644 --- a/deal.II/lac/source/filtered_matrix.cc +++ b/deal.II/lac/source/filtered_matrix.cc @@ -14,98 +14,18 @@ #include +#define FILT(MM,VV) \ +template <> \ +void FilteredMatrix::allocate_tmp_vector() \ +{\ + Threads::ThreadMutex::ScopedLock lock (tmp_mutex); \ + tmp_vector.reinit (matrix->n(), true); \ +} -// Following, we define some partial specializations. since C++ does not allow -// this properly, we have to work around it using some macro ugliness - -#define DEFINE_SPECIALIZATIONS(TYPE) \ -template <> \ -void \ -FilteredMatrix,Vector >:: \ -get_column_entries (const unsigned int index, \ - std::vector &column_entries, \ - const bool matrix_is_symmetric) const \ -{ \ - /* depending on whether the matrix */ \ - /* can be assumed symmetric or not,*/ \ - /* either use a fast or a slow */ \ - /* algorithm */ \ - if (matrix_is_symmetric == true) \ - /* ok, matrix is symmetric. we */ \ - /* may determine the matrix */ \ - /* entries in this column by */ \ - /* looking at the matrix entries */ \ - /* in this row which is */ \ - /* significantly faster since we */ \ - /* can traverse them linearly and */ \ - /* do not have to check each row */ \ - /* for the possible existence of */ \ - /* a matrix entry */ \ - { \ - const unsigned int * \ - col_nums = &(matrix->get_sparsity_pattern().get_column_numbers() \ - [matrix->get_sparsity_pattern().get_rowstart_indices()[index]]); \ - const unsigned int \ - row_length = matrix->get_sparsity_pattern().row_length(index); \ - \ - for (unsigned int i=0; iget_sparsity_pattern()(row,index); \ - if (global_index != SparsityPattern::invalid_entry) \ - column_entries.push_back (std::make_pair(row, \ - (*matrix)(row,index))); \ - }; \ - }; \ -} \ - \ - \ - \ -template <> \ -void \ -FilteredMatrix,BlockVector >:: \ -get_column_entries (const unsigned int /*index*/, \ - std::vector &/*column_entries*/, \ - const bool /*matrix_is_symmetric*/) const \ -{ \ - /* presently not implemented, but */ \ - /* should be fairly simple to do */ \ - Assert (false, ExcNotImplemented()); \ -} \ - \ - \ - \ - \ +#define BFILT(MM,VV) \ template <> \ void \ -FilteredMatrix,Vector >:: \ -allocate_tmp_vector () \ -{ \ - Threads::ThreadMutex::ScopedLock lock (tmp_mutex); \ - tmp_vector.reinit (matrix->n(), true); \ -} \ - \ - \ - \ -template <> \ -void \ -FilteredMatrix,BlockVector >:: \ +FilteredMatrix:: \ allocate_tmp_vector () \ { \ std::vector block_sizes (matrix->n_block_rows()); \ @@ -116,13 +36,12 @@ allocate_tmp_vector () \ tmp_vector.reinit (block_sizes, true); \ } - -// now instantiate -DEFINE_SPECIALIZATIONS(double) +FILT(SparseMatrix, Vector) +FILT(SparseMatrix, Vector) template class FilteredMatrix,Vector >; template class FilteredMatrix,BlockVector >; -DEFINE_SPECIALIZATIONS(float) +BFILT(BlockSparseMatrix, BlockVector) +BFILT(BlockSparseMatrix, BlockVector) template class FilteredMatrix,Vector >; template class FilteredMatrix,BlockVector >; - -- 2.39.5