* <h3>Template arguments</h3>
*
* 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 <tt>operator()</tt>,
* 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 MATRIX, class VECTOR=Vector<typename MATRIX::value_type> >
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
- * <tt>std::map<unsigned,value_type></tt>,
- * 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 <class ConstraintList>
void add_constraints (const ConstraintList &new_constraints);
*/
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<IndexValuePair> &column_entries,
- const bool matrix_is_symmetric) const;
};
/*@}*/
#include <lac/vector.h>
#include <lac/block_vector.h>
-
template <class MATRIX, class VECTOR>
FilteredMatrix<MATRIX,VECTOR>::FilteredMatrix ()
{}
-
-
template <class MATRIX, class VECTOR>
FilteredMatrix<MATRIX,VECTOR>::
FilteredMatrix (const FilteredMatrix &fm)
set_referenced_matrix (const MATRIX &m)
{
matrix = &m;
- allocate_tmp_vector ();
+ allocate_tmp_vector();
}
void
FilteredMatrix<MATRIX,VECTOR>::
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<IndexValuePair> 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;
}
void
FilteredMatrix<MATRIX,VECTOR>::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;
}
template <class MATRIX, class VECTOR>
void
FilteredMatrix<MATRIX,VECTOR>::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);
}
template <class MATRIX, class VECTOR>
void
FilteredMatrix<MATRIX,VECTOR>::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);
}
template <class MATRIX, class VECTOR>
typename FilteredMatrix<MATRIX,VECTOR>::value_type
FilteredMatrix<MATRIX,VECTOR>::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);
}
template <class MATRIX, class VECTOR>
void
FilteredMatrix<MATRIX,VECTOR>::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 <class MATRIX, class VECTOR>
typename FilteredMatrix<MATRIX,VECTOR>::value_type
FilteredMatrix<MATRIX,VECTOR>::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;
}
void
FilteredMatrix<MATRIX,VECTOR>::
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);
}
unsigned int
FilteredMatrix<MATRIX,VECTOR>::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));
}
#include <lac/filtered_matrix.templates.h>
+#define FILT(MM,VV) \
+template <> \
+void FilteredMatrix<MM , VV >::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<SparseMatrix<TYPE>,Vector<TYPE> >:: \
-get_column_entries (const unsigned int index, \
- std::vector<IndexValuePair> &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; i<row_length; ++i) \
- { \
- const unsigned int c = *(col_nums+i); \
- \
- /* if not diagonal entry,*/ \
- /* add to list */ \
- if (c != index) \
- column_entries.push_back (std::make_pair(c, (*matrix)(c,index))); \
- }; \
- } \
- else \
- { \
- /* otherwise check each row for */ \
- /* occurrence of an entry in */ \
- /* this column */ \
- for (unsigned int row=0; row<n(); ++row) \
- if (row != index) \
- { \
- const unsigned int \
- global_index = matrix->get_sparsity_pattern()(row,index); \
- if (global_index != SparsityPattern::invalid_entry) \
- column_entries.push_back (std::make_pair(row, \
- (*matrix)(row,index))); \
- }; \
- }; \
-} \
- \
- \
- \
-template <> \
-void \
-FilteredMatrix<BlockSparseMatrix<TYPE>,BlockVector<TYPE> >:: \
-get_column_entries (const unsigned int /*index*/, \
- std::vector<IndexValuePair> &/*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<SparseMatrix<TYPE>,Vector<TYPE> >:: \
-allocate_tmp_vector () \
-{ \
- Threads::ThreadMutex::ScopedLock lock (tmp_mutex); \
- tmp_vector.reinit (matrix->n(), true); \
-} \
- \
- \
- \
-template <> \
-void \
-FilteredMatrix<BlockSparseMatrix<TYPE>,BlockVector<TYPE> >:: \
+FilteredMatrix<MM ,VV >:: \
allocate_tmp_vector () \
{ \
std::vector<unsigned int> block_sizes (matrix->n_block_rows()); \
tmp_vector.reinit (block_sizes, true); \
}
-
-// now instantiate
-DEFINE_SPECIALIZATIONS(double)
+FILT(SparseMatrix<double>, Vector<double>)
+FILT(SparseMatrix<float>, Vector<float>)
template class FilteredMatrix<SparseMatrix<double>,Vector<double> >;
template class FilteredMatrix<BlockSparseMatrix<double>,BlockVector<double> >;
-DEFINE_SPECIALIZATIONS(float)
+BFILT(BlockSparseMatrix<double>, BlockVector<double>)
+BFILT(BlockSparseMatrix<float>, BlockVector<float>)
template class FilteredMatrix<SparseMatrix<float>,Vector<float> >;
template class FilteredMatrix<BlockSparseMatrix<float>,BlockVector<float> >;
-