include/deal.II/lac/filtered_matrix.h

00001 //---------------------------------------------------------------------------
00002 //    @f$Id: filtered_matrix.h 25345 2012-03-31 08:37:04Z bangerth @f$
00003 //
00004 //    Copyright (C) 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012 by the deal.II authors
00005 //
00006 //    This file is subject to QPL and may not be  distributed
00007 //    without copyright and license information. Please refer
00008 //    to the file deal.II/doc/license.html for the  text  and
00009 //    further information on this license.
00010 //
00011 //---------------------------------------------------------------------------
00012 #ifndef __deal2__filtered_matrix_h
00013 #define __deal2__filtered_matrix_h
00014 
00015 
00016 
00017 #include <deal.II/base/config.h>
00018 #include <deal.II/base/smartpointer.h>
00019 #include <deal.II/base/thread_management.h>
00020 #include <deal.II/base/memory_consumption.h>
00021 #include <deal.II/lac/pointer_matrix.h>
00022 #include <deal.II/lac/vector_memory.h>
00023 #include <vector>
00024 #include <algorithm>
00025 
00026 DEAL_II_NAMESPACE_OPEN
00027 
00028 template <typename number> class Vector;
00029 template <class VECTOR> class FilteredMatrixBlock;
00030 
00031 
00203 template <class VECTOR>
00204 class FilteredMatrix : public Subscriptor
00205 {
00206   public:
00207         class const_iterator;
00208 
00212     class Accessor
00213     {
00220         Accessor (const FilteredMatrix<VECTOR> *matrix,
00221                   const unsigned int index);
00222 
00223       public:
00229         unsigned int row() const;
00230 
00236         unsigned int column() const;
00237 
00242         double value() const;
00243 
00244       private:
00248         void advance ();
00249 
00253         const FilteredMatrix<VECTOR>* matrix;
00254 
00258         unsigned int index;
00259                                          /*
00260                                           * Make enclosing class a
00261                                           * friend.
00262                                           */
00263         friend class const_iterator;
00264     };
00265 
00269     class const_iterator
00270     {
00271       public:
00275         const_iterator(const FilteredMatrix<VECTOR> *matrix,
00276                        const unsigned int index);
00277 
00281         const_iterator& operator++ ();
00282 
00286         const_iterator& operator++ (int);
00287 
00291         const Accessor& operator* () const;
00292 
00296         const Accessor* operator-> () const;
00297 
00304         bool operator == (const const_iterator&) const;
00308         bool operator != (const const_iterator&) const;
00309 
00317         bool operator < (const const_iterator&) const;
00318 
00324         bool operator > (const const_iterator&) const;
00325 
00326       private:
00331         Accessor accessor;
00332     };
00333 
00340     typedef std::pair<unsigned int, double> IndexValuePair;
00341 
00352     FilteredMatrix ();
00353 
00360     FilteredMatrix (const FilteredMatrix &fm);
00361 
00373     template <class MATRIX>
00374     FilteredMatrix (const MATRIX &matrix,
00375                     bool expect_constrained_source = false);
00376 
00382     FilteredMatrix & operator = (const FilteredMatrix &fm);
00383 
00399     template <class MATRIX>
00400     void initialize (const MATRIX &m,
00401                      bool expect_constrained_source = false);
00402 
00407     void clear ();
00409 
00419     void add_constraint (const unsigned int i, const double v);
00420 
00451     template <class ConstraintList>
00452     void add_constraints (const ConstraintList &new_constraints);
00453 
00458     void clear_constraints ();
00460 
00476     void apply_constraints (VECTOR     &v,
00477                             const bool  matrix_is_symmetric) const;
00478 
00486     void vmult (VECTOR       &dst,
00487                 const VECTOR &src) const;
00488 
00497     void Tvmult (VECTOR       &dst,
00498                  const VECTOR &src) const;
00499 
00512     void vmult_add (VECTOR       &dst,
00513                     const VECTOR &src) const;
00514 
00527     void Tvmult_add (VECTOR       &dst,
00528                      const VECTOR &src) const;
00530 
00539     const_iterator begin () const;
00543     const_iterator end () const;
00545 
00554     std::size_t memory_consumption () const;
00555 
00556   private:
00576     bool expect_constrained_source;
00577 
00586     typedef typename std::vector<IndexValuePair>::const_iterator const_index_value_iterator;
00587 
00594     struct PairComparison
00595     {
00601         bool operator () (const IndexValuePair &i1,
00602                           const IndexValuePair &i2) const;
00603     };
00604 
00610     std_cxx1x::shared_ptr<PointerMatrixBase<VECTOR> > matrix;
00611 
00618     std::vector<IndexValuePair> constraints;
00619 
00626     void pre_filter (VECTOR &v) const;
00627 
00637     void post_filter (const VECTOR &in,
00638                       VECTOR       &out) const;
00639 
00640     friend class Accessor;
00645     friend class FilteredMatrixBlock<VECTOR>;
00646 };
00647 
00649 /*---------------------- Inline functions -----------------------------------*/
00650 
00651 
00652 //--------------------------------Iterators--------------------------------------//
00653 
00654 template<class VECTOR>
00655 inline
00656 FilteredMatrix<VECTOR>::Accessor::Accessor(
00657   const FilteredMatrix<VECTOR> *matrix,
00658   const unsigned int index)
00659                 :
00660                 matrix(matrix),
00661                 index(index)
00662 {
00663   Assert (index <= matrix->constraints.size(),
00664           ExcIndexRange(index, 0, matrix->constraints.size()));
00665 }
00666 
00667 
00668 
00669 template<class VECTOR>
00670 inline
00671 unsigned int
00672 FilteredMatrix<VECTOR>::Accessor::row() const
00673 {
00674   return matrix->constraints[index].first;
00675 }
00676 
00677 
00678 
00679 template<class VECTOR>
00680 inline
00681 unsigned int
00682 FilteredMatrix<VECTOR>::Accessor::column() const
00683 {
00684   return matrix->constraints[index].first;
00685 }
00686 
00687 
00688 
00689 template<class VECTOR>
00690 inline
00691 double
00692 FilteredMatrix<VECTOR>::Accessor::value() const
00693 {
00694   return matrix->constraints[index].second;
00695 }
00696 
00697 
00698 
00699 template<class VECTOR>
00700 inline
00701 void
00702 FilteredMatrix<VECTOR>::Accessor::advance()
00703 {
00704   Assert (index < matrix->constraints.size(), ExcIteratorPastEnd());
00705   ++index;
00706 }
00707 
00708 
00709 
00710 
00711 template<class VECTOR>
00712 inline
00713 FilteredMatrix<VECTOR>::const_iterator::const_iterator(
00714   const FilteredMatrix<VECTOR> *matrix,
00715   const unsigned int index)
00716                 :
00717                 accessor(matrix, index)
00718 {}
00719 
00720 
00721 
00722 template<class VECTOR>
00723 inline
00724 typename FilteredMatrix<VECTOR>::const_iterator&
00725 FilteredMatrix<VECTOR>::const_iterator::operator++ ()
00726 {
00727   accessor.advance();
00728   return *this;
00729 }
00730 
00731 
00732 template <typename number>
00733 inline
00734 const typename FilteredMatrix<number>::Accessor &
00735 FilteredMatrix<number>::const_iterator::operator* () const
00736 {
00737   return accessor;
00738 }
00739 
00740 
00741 template <typename number>
00742 inline
00743 const typename FilteredMatrix<number>::Accessor *
00744 FilteredMatrix<number>::const_iterator::operator-> () const
00745 {
00746   return &accessor;
00747 }
00748 
00749 
00750 template <typename number>
00751 inline
00752 bool
00753 FilteredMatrix<number>::const_iterator::
00754 operator == (const const_iterator& other) const
00755 {
00756   return (accessor.index == other.accessor.index
00757           && accessor.matrix == other.accessor.matrix);
00758 }
00759 
00760 
00761 template <typename number>
00762 inline
00763 bool
00764 FilteredMatrix<number>::const_iterator::
00765 operator != (const const_iterator& other) const
00766 {
00767   return ! (*this == other);
00768 }
00769 
00770 
00771 
00772 //------------------------------- FilteredMatrix ---------------------------------------//
00773 
00774 template <typename number>
00775 inline
00776 typename FilteredMatrix<number>::const_iterator
00777 FilteredMatrix<number>::begin () const
00778 {
00779   return const_iterator(this, 0);
00780 }
00781 
00782 
00783 template <typename number>
00784 inline
00785 typename FilteredMatrix<number>::const_iterator
00786 FilteredMatrix<number>::end () const
00787 {
00788   return const_iterator(this, constraints.size());
00789 }
00790 
00791 
00792 template <class VECTOR>
00793 inline
00794 bool
00795 FilteredMatrix<VECTOR>::PairComparison::
00796 operator () (const IndexValuePair &i1,
00797              const IndexValuePair &i2) const
00798 {
00799   return (i1.first < i2.first);
00800 }
00801 
00802 
00803 
00804 template <class VECTOR>
00805 template <class MATRIX>
00806 inline
00807 void
00808 FilteredMatrix<VECTOR>::initialize (const MATRIX &m, bool ecs)
00809 {
00810   matrix.reset (new_pointer_matrix_base(m, VECTOR()));
00811 
00812   expect_constrained_source = ecs;
00813 }
00814 
00815 
00816 
00817 template <class VECTOR>
00818 inline
00819 FilteredMatrix<VECTOR>::FilteredMatrix ()
00820 {}
00821 
00822 
00823 
00824 template <class VECTOR>
00825 inline
00826 FilteredMatrix<VECTOR>::FilteredMatrix (const FilteredMatrix &fm)
00827                 :
00828                 Subscriptor(),
00829                 expect_constrained_source(fm.expect_constrained_source),
00830                 matrix(fm.matrix),
00831                 constraints (fm.constraints)
00832 {}
00833 
00834 
00835 
00836 template <class VECTOR>
00837 template <class MATRIX>
00838 inline
00839 FilteredMatrix<VECTOR>::
00840 FilteredMatrix (const MATRIX &m, bool ecs)
00841 {
00842   initialize (m, ecs);
00843 }
00844 
00845 
00846 
00847 template <class VECTOR>
00848 inline
00849 FilteredMatrix<VECTOR> &
00850 FilteredMatrix<VECTOR>::operator = (const FilteredMatrix &fm)
00851 {
00852   matrix = fm.matrix;
00853   expect_constrained_source = fm.expect_constrained_source;
00854   constraints = fm.constraints;
00855   return *this;
00856 }
00857 
00858 
00859 
00860 template <class VECTOR>
00861 inline
00862 void
00863 FilteredMatrix<VECTOR>::add_constraint (const unsigned int index, const double value)
00864 {
00865                                    // add new constraint to end
00866   constraints.push_back(IndexValuePair(index, value));
00867 }
00868 
00869 
00870 
00871 template <class VECTOR>
00872 template <class ConstraintList>
00873 inline
00874 void
00875 FilteredMatrix<VECTOR>::add_constraints (const ConstraintList &new_constraints)
00876 {
00877                                    // add new constraints to end
00878   const unsigned int old_size = constraints.size();
00879   constraints.reserve (old_size + new_constraints.size());
00880   constraints.insert (constraints.end(),
00881                       new_constraints.begin(),
00882                       new_constraints.end());
00883                                    // then merge the two arrays to
00884                                    // form one sorted one
00885   std::inplace_merge (constraints.begin(),
00886                       constraints.begin()+old_size,
00887                       constraints.end(),
00888                       PairComparison());
00889 }
00890 
00891 
00892 
00893 template <class VECTOR>
00894 inline
00895 void
00896 FilteredMatrix<VECTOR>::clear_constraints ()
00897 {
00898                                    // swap vectors to release memory
00899   std::vector<IndexValuePair> empty;
00900   constraints.swap (empty);
00901 }
00902 
00903 
00904 
00905 template <class VECTOR>
00906 inline
00907 void
00908 FilteredMatrix<VECTOR>::clear ()
00909 {
00910   clear_constraints();
00911   matrix.reset();
00912 }
00913 
00914 
00915 
00916 template <class VECTOR>
00917 inline
00918 void
00919 FilteredMatrix<VECTOR>::apply_constraints (
00920   VECTOR     &v,
00921   const bool  /* matrix_is_symmetric */) const
00922 {
00923   GrowingVectorMemory<VECTOR> mem;
00924   typename VectorMemory<VECTOR>::Pointer tmp_vector(mem);
00925   tmp_vector->reinit(v);
00926   const_index_value_iterator       i = constraints.begin();
00927   const const_index_value_iterator e = constraints.end();
00928   for (; i!=e; ++i)
00929     {
00930       Assert(numbers::is_finite(i->second), ExcNumberNotFinite());
00931       (*tmp_vector)(i->first) = -i->second;
00932     }
00933 
00934                                    // This vmult is without bc, to get
00935                                    // the rhs correction in a correct
00936                                    // way.
00937   matrix->vmult_add(v, *tmp_vector);
00938                                    // finally set constrained
00939                                    // entries themselves
00940   for (i=constraints.begin(); i!=e; ++i)
00941     {
00942       Assert(numbers::is_finite(i->second), ExcNumberNotFinite());
00943       v(i->first) = i->second;
00944     }
00945 }
00946 
00947 
00948 
00949 template <class VECTOR>
00950 inline
00951 void
00952 FilteredMatrix<VECTOR>::pre_filter (VECTOR &v) const
00953 {
00954     // iterate over all constraints and
00955     // zero out value
00956     const_index_value_iterator       i = constraints.begin();
00957     const const_index_value_iterator e = constraints.end();
00958     for (; i!=e; ++i)
00959         v(i->first) = 0;
00960 }
00961 
00962 
00963 
00964 template <class VECTOR>
00965 inline
00966 void
00967 FilteredMatrix<VECTOR>::post_filter (const VECTOR &in,
00968         VECTOR       &out) const
00969 {
00970     // iterate over all constraints and
00971     // set value correctly
00972   const_index_value_iterator       i = constraints.begin();
00973     const const_index_value_iterator e = constraints.end();
00974     for (; i!=e; ++i)
00975       {
00976         Assert(numbers::is_finite(in(i->first)), ExcNumberNotFinite());
00977         out(i->first) = in(i->first);
00978       }
00979 }
00980 
00981 
00982 
00983 template <class VECTOR>
00984 inline
00985 void
00986 FilteredMatrix<VECTOR>::vmult (VECTOR& dst, const VECTOR& src) const
00987 {
00988   if (!expect_constrained_source)
00989     {
00990       GrowingVectorMemory<VECTOR> mem;
00991       VECTOR* tmp_vector = mem.alloc();
00992                                        // first copy over src vector and
00993                                        // pre-filter
00994       tmp_vector->reinit(src, true);
00995       *tmp_vector = src;
00996       pre_filter (*tmp_vector);
00997                                        // then let matrix do its work
00998       matrix->vmult (dst, *tmp_vector);
00999       mem.free(tmp_vector);
01000     }
01001   else
01002     {
01003       matrix->vmult (dst, src);
01004     }
01005 
01006                                    // finally do post-filtering
01007   post_filter (src, dst);
01008 }
01009 
01010 
01011 
01012 template <class VECTOR>
01013 inline
01014 void
01015 FilteredMatrix<VECTOR>::Tvmult (VECTOR& dst, const VECTOR& src) const
01016 {
01017   if (!expect_constrained_source)
01018     {
01019       GrowingVectorMemory<VECTOR> mem;
01020       VECTOR* tmp_vector = mem.alloc();
01021                                        // first copy over src vector and
01022                                        // pre-filter
01023       tmp_vector->reinit(src, true);
01024       *tmp_vector = src;
01025       pre_filter (*tmp_vector);
01026                                        // then let matrix do its work
01027       matrix->Tvmult (dst, *tmp_vector);
01028       mem.free(tmp_vector);
01029     }
01030   else
01031     {
01032       matrix->Tvmult (dst, src);
01033     }
01034 
01035                                    // finally do post-filtering
01036   post_filter (src, dst);
01037 }
01038 
01039 
01040 
01041 template <class VECTOR>
01042 inline
01043 void
01044 FilteredMatrix<VECTOR>::vmult_add (VECTOR& dst, const VECTOR& src) const
01045 {
01046   if (!expect_constrained_source)
01047     {
01048       GrowingVectorMemory<VECTOR> mem;
01049       VECTOR* tmp_vector = mem.alloc();
01050                                        // first copy over src vector and
01051                                        // pre-filter
01052       tmp_vector->reinit(src, true);
01053       *tmp_vector = src;
01054       pre_filter (*tmp_vector);
01055                                        // then let matrix do its work
01056       matrix->vmult_add (dst, *tmp_vector);
01057       mem.free(tmp_vector);
01058     }
01059   else
01060     {
01061       matrix->vmult_add (dst, src);
01062     }
01063 
01064                                    // finally do post-filtering
01065   post_filter (src, dst);
01066 }
01067 
01068 
01069 
01070 template <class VECTOR>
01071 inline
01072 void
01073 FilteredMatrix<VECTOR>::Tvmult_add (VECTOR& dst, const VECTOR& src) const
01074 {
01075   if (!expect_constrained_source)
01076     {
01077       GrowingVectorMemory<VECTOR> mem;
01078       VECTOR* tmp_vector = mem.alloc();
01079                                        // first copy over src vector and
01080                                        // pre-filter
01081       tmp_vector->reinit(src, true);
01082       *tmp_vector = src;
01083       pre_filter (*tmp_vector);
01084                                        // then let matrix do its work
01085       matrix->Tvmult_add (dst, *tmp_vector);
01086       mem.free(tmp_vector);
01087     }
01088   else
01089     {
01090       matrix->Tvmult_add (dst, src);
01091     }
01092 
01093                                    // finally do post-filtering
01094   post_filter (src, dst);
01095 }
01096 
01097 
01098 
01099 template <class VECTOR>
01100 inline
01101 std::size_t
01102 FilteredMatrix<VECTOR>::memory_consumption () const
01103 {
01104   return (MemoryConsumption::memory_consumption (matrix) +
01105           MemoryConsumption::memory_consumption (constraints));
01106 }
01107 
01108 
01109 
01110 DEAL_II_NAMESPACE_CLOSE
01111 
01112 #endif
01113 /*----------------------------   filtered_matrix.h     ---------------------------*/
01114 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

deal.II documentation generated on Tue May 22 2012 12:06:11 by doxygen 1.7.3