00001
00002
00003
00004
00005
00006
00007
00008
00009
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
00261
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
00650
00651
00652
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
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
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
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
00884
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
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 ) 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
00935
00936
00937 matrix->vmult_add(v, *tmp_vector);
00938
00939
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
00955
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
00971
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
00993
00994 tmp_vector->reinit(src, true);
00995 *tmp_vector = src;
00996 pre_filter (*tmp_vector);
00997
00998 matrix->vmult (dst, *tmp_vector);
00999 mem.free(tmp_vector);
01000 }
01001 else
01002 {
01003 matrix->vmult (dst, src);
01004 }
01005
01006
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
01022
01023 tmp_vector->reinit(src, true);
01024 *tmp_vector = src;
01025 pre_filter (*tmp_vector);
01026
01027 matrix->Tvmult (dst, *tmp_vector);
01028 mem.free(tmp_vector);
01029 }
01030 else
01031 {
01032 matrix->Tvmult (dst, src);
01033 }
01034
01035
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
01051
01052 tmp_vector->reinit(src, true);
01053 *tmp_vector = src;
01054 pre_filter (*tmp_vector);
01055
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
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
01080
01081 tmp_vector->reinit(src, true);
01082 *tmp_vector = src;
01083 pre_filter (*tmp_vector);
01084
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
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
01114