]> https://gitweb.dealii.org/ - dealii.git/commitdiff
New Version of Filtered Matrix, which requires only vector matrix mult.
authorLuca Heltai <luca.heltai@sissa.it>
Tue, 10 Jan 2006 00:34:22 +0000 (00:34 +0000)
committerLuca Heltai <luca.heltai@sissa.it>
Tue, 10 Jan 2006 00:34:22 +0000 (00:34 +0000)
git-svn-id: https://svn.dealii.org/trunk@11979 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/include/lac/filtered_matrix.h
deal.II/lac/include/lac/filtered_matrix.templates.h
deal.II/lac/source/filtered_matrix.cc

index 85085f8699b1eabb82e9fcd066b4be116711b7ff..52464ec18650282fda5cad360651080bff1fa736 100644 (file)
@@ -169,7 +169,7 @@ template <typename number> class Vector;
  * <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>,
@@ -190,7 +190,7 @@ template <typename number> 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 MATRIX, class VECTOR=Vector<typename MATRIX::value_type> >
 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
-                                     * <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);
@@ -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<IndexValuePair> &column_entries,
-                            const bool                   matrix_is_symmetric) const;
 };
 
 /*@}*/
index f5947de0b220f5dcab3d07b1c2b41f8183ce6b42..2a3131f22e1e7c67504efe0cee9b9796cb9f288e 100644 (file)
 #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)
@@ -67,7 +64,7 @@ FilteredMatrix<MATRIX,VECTOR>::
 set_referenced_matrix (const MATRIX &m)
 {
   matrix = &m;
-  allocate_tmp_vector ();
+  allocate_tmp_vector();
 }
 
 
@@ -87,55 +84,20 @@ template <class MATRIX, class 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;
 }
 
 
@@ -144,12 +106,12 @@ template <class MATRIX, class VECTOR>
 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;
 }
 
 
@@ -157,14 +119,14 @@ FilteredMatrix<MATRIX,VECTOR>::pre_filter (VECTOR &v) const
 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);
 }
 
 
@@ -172,19 +134,19 @@ FilteredMatrix<MATRIX,VECTOR>::post_filter (const VECTOR &in,
 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);
 }
 
 
@@ -192,34 +154,34 @@ FilteredMatrix<MATRIX,VECTOR>::vmult (VECTOR       &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);
 }
 
 
@@ -227,41 +189,41 @@ FilteredMatrix<MATRIX,VECTOR>::residual (VECTOR       &dst,
 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;
 }
 
 
@@ -270,23 +232,23 @@ template <class MATRIX, class VECTOR>
 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);
 }
 
 
@@ -295,9 +257,9 @@ template <class MATRIX, class VECTOR>
 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));
 }
 
 
index 3bd0d036f6211d15b3cc5ca1dae508bd5a1d89cc..bc6421cde0ac3b46c451b922738719d0fd1693fe 100644 (file)
 
 #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());   \
@@ -116,13 +36,12 @@ allocate_tmp_vector ()    \
   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> >;
-

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.