00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 #ifndef __deal2__mg_smoother_h
00013 #define __deal2__mg_smoother_h
00014
00015
00016 #include <deal.II/base/config.h>
00017 #include <deal.II/base/smartpointer.h>
00018 #include <deal.II/lac/pointer_matrix.h>
00019 #include <deal.II/lac/vector_memory.h>
00020 #include <deal.II/multigrid/mg_base.h>
00021 #include <deal.II/base/mg_level_object.h>
00022 #include <vector>
00023
00024 DEAL_II_NAMESPACE_OPEN
00025
00026 template <int dim, int spacedim> class MGDoFHandler;
00027
00028
00029
00030
00031
00032
00035
00044 template <class VECTOR>
00045 class MGSmoother : public MGSmootherBase<VECTOR>
00046 {
00047 public:
00054 MGSmoother(const unsigned int steps = 1,
00055 const bool variable = false,
00056 const bool symmetric = false,
00057 const bool transpose = false);
00058
00070 MGSmoother(VectorMemory<VECTOR>& mem,
00071 const unsigned int steps = 1,
00072 const bool variable = false,
00073 const bool symmetric = false,
00074 const bool transpose = false);
00075
00080 void set_steps (const unsigned int);
00081
00086 void set_variable (const bool);
00087
00092 void set_symmetric (const bool);
00093
00099 void set_transpose (const bool);
00100
00107 void set_debug (const unsigned int level);
00108
00109 private:
00115 GrowingVectorMemory<VECTOR> my_memory;
00116
00117 protected:
00125 unsigned int steps;
00126
00133 bool variable;
00134
00141 bool symmetric;
00142
00143
00144
00145
00146
00147
00148
00149
00150 bool transpose;
00151
00157 unsigned int debug;
00161 SmartPointer<VectorMemory<VECTOR>, MGSmoother<VECTOR> > mem;
00162
00163 };
00164
00165
00174 template <class VECTOR>
00175 class MGSmootherIdentity : public MGSmootherBase<VECTOR>
00176 {
00177 public:
00188 virtual void smooth (const unsigned int level,
00189 VECTOR& u,
00190 const VECTOR& rhs) const;
00191 virtual void clear ();
00192 };
00193
00194
00195 namespace mg
00196 {
00231 template<class RELAX, class VECTOR>
00232 class SmootherRelaxation : public MGLevelObject<RELAX>, public MGSmoother<VECTOR>
00233 {
00234 public:
00239 SmootherRelaxation(const unsigned int steps = 1,
00240 const bool variable = false,
00241 const bool symmetric = false,
00242 const bool transpose = false);
00243
00257 template <class MATRIX2>
00258 void initialize (const MGLevelObject<MATRIX2>& matrices,
00259 const typename RELAX::AdditionalData & additional_data = typename RELAX::AdditionalData());
00260
00274 template <class MATRIX2, class DATA>
00275 void initialize (const MGLevelObject<MATRIX2>& matrices,
00276 const MGLevelObject<DATA>& additional_data);
00291
00292
00293
00294
00298 void clear ();
00299
00303 virtual void smooth (const unsigned int level,
00304 VECTOR& u,
00305 const VECTOR& rhs) const;
00306
00310 std::size_t memory_consumption () const;
00311 };
00312 }
00313
00355 template<class MATRIX, class RELAX, class VECTOR>
00356 class MGSmootherRelaxation : public MGSmoother<VECTOR>
00357 {
00358 public:
00363 MGSmootherRelaxation(VectorMemory<VECTOR>& mem,
00364 const unsigned int steps = 1,
00365 const bool variable = false,
00366 const bool symmetric = false,
00367 const bool transpose = false);
00368
00384 template <class MATRIX2>
00385 void initialize (const MGLevelObject<MATRIX2>& matrices,
00386 const typename RELAX::AdditionalData & additional_data = typename RELAX::AdditionalData());
00387
00403 template <class MATRIX2, class DATA>
00404 void initialize (const MGLevelObject<MATRIX2>& matrices,
00405 const MGLevelObject<DATA>& additional_data);
00406
00427 template <class MATRIX2, class DATA>
00428 void initialize (const MGLevelObject<MATRIX2>& matrices,
00429 const DATA& additional_data,
00430 const unsigned int block_row,
00431 const unsigned int block_col);
00432
00453 template <class MATRIX2, class DATA>
00454 void initialize (const MGLevelObject<MATRIX2>& matrices,
00455 const MGLevelObject<DATA>& additional_data,
00456 const unsigned int block_row,
00457 const unsigned int block_col);
00458
00462 void clear ();
00463
00467 virtual void smooth (const unsigned int level,
00468 VECTOR& u,
00469 const VECTOR& rhs) const;
00470
00475 MGLevelObject<RELAX> smoothers;
00476
00480 std::size_t memory_consumption () const;
00481
00482
00483 private:
00487 MGLevelObject<PointerMatrix<MATRIX, VECTOR> > matrices;
00488
00489 };
00490
00491
00492
00524 template<class MATRIX, class PRECONDITIONER, class VECTOR>
00525 class MGSmootherPrecondition : public MGSmoother<VECTOR>
00526 {
00527 public:
00532 MGSmootherPrecondition(VectorMemory<VECTOR>& mem,
00533 const unsigned int steps = 1,
00534 const bool variable = false,
00535 const bool symmetric = false,
00536 const bool transpose = false);
00537
00553 template <class MATRIX2>
00554 void initialize (const MGLevelObject<MATRIX2>& matrices,
00555 const typename PRECONDITIONER::AdditionalData& additional_data = typename PRECONDITIONER::AdditionalData());
00556
00572 template <class MATRIX2, class DATA>
00573 void initialize (const MGLevelObject<MATRIX2>& matrices,
00574 const MGLevelObject<DATA>& additional_data);
00575
00596 template <class MATRIX2, class DATA>
00597 void initialize (const MGLevelObject<MATRIX2>& matrices,
00598 const DATA& additional_data,
00599 const unsigned int block_row,
00600 const unsigned int block_col);
00601
00622 template <class MATRIX2, class DATA>
00623 void initialize (const MGLevelObject<MATRIX2>& matrices,
00624 const MGLevelObject<DATA>& additional_data,
00625 const unsigned int block_row,
00626 const unsigned int block_col);
00627
00631 void clear ();
00632
00636 virtual void smooth (const unsigned int level,
00637 VECTOR& u,
00638 const VECTOR& rhs) const;
00639
00644 MGLevelObject<PRECONDITIONER> smoothers;
00645
00649 std::size_t memory_consumption () const;
00650
00651
00652 private:
00656 MGLevelObject<PointerMatrix<MATRIX, VECTOR> > matrices;
00657
00658 };
00659
00662
00663
00664 #ifndef DOXYGEN
00665
00666 template <class VECTOR>
00667 inline void
00668 MGSmootherIdentity<VECTOR>::smooth (
00669 const unsigned int, VECTOR&,
00670 const VECTOR&) const
00671 {}
00672
00673 template <class VECTOR>
00674 inline void
00675 MGSmootherIdentity<VECTOR>::clear ()
00676 {}
00677
00678
00679
00680 template <class VECTOR>
00681 inline
00682 MGSmoother<VECTOR>::MGSmoother(
00683 const unsigned int steps,
00684 const bool variable,
00685 const bool symmetric,
00686 const bool transpose)
00687 :
00688 steps(steps),
00689 variable(variable),
00690 symmetric(symmetric),
00691 transpose(transpose),
00692 debug(0),
00693 mem(&my_memory)
00694 {}
00695
00696
00697 template <class VECTOR>
00698 inline
00699 MGSmoother<VECTOR>::MGSmoother(
00700 VectorMemory<VECTOR>& mem,
00701 const unsigned int steps,
00702 const bool variable,
00703 const bool symmetric,
00704 const bool transpose)
00705 :
00706 steps(steps),
00707 variable(variable),
00708 symmetric(symmetric),
00709 transpose(transpose),
00710 debug(0),
00711 mem(&mem)
00712 {}
00713
00714
00715 template <class VECTOR>
00716 inline void
00717 MGSmoother<VECTOR>::set_steps (const unsigned int s)
00718 {
00719 steps = s;
00720 }
00721
00722
00723 template <class VECTOR>
00724 inline void
00725 MGSmoother<VECTOR>::set_debug (const unsigned int s)
00726 {
00727 debug = s;
00728 }
00729
00730
00731 template <class VECTOR>
00732 inline void
00733 MGSmoother<VECTOR>::set_variable (const bool flag)
00734 {
00735 variable = flag;
00736 }
00737
00738
00739 template <class VECTOR>
00740 inline void
00741 MGSmoother<VECTOR>::set_symmetric (const bool flag)
00742 {
00743 symmetric = flag;
00744 }
00745
00746
00747 template <class VECTOR>
00748 inline void
00749 MGSmoother<VECTOR>::set_transpose (const bool flag)
00750 {
00751 transpose = flag;
00752 }
00753
00754
00755
00756 namespace mg
00757 {
00758 template <class RELAX, class VECTOR>
00759 inline
00760 SmootherRelaxation<RELAX, VECTOR>::SmootherRelaxation(
00761 const unsigned int steps,
00762 const bool variable,
00763 const bool symmetric,
00764 const bool transpose)
00765 : MGSmoother<VECTOR>(steps, variable, symmetric, transpose)
00766 {}
00767
00768
00769 template <class RELAX, class VECTOR>
00770 inline void
00771 SmootherRelaxation<RELAX, VECTOR>::clear ()
00772 {
00773 MGLevelObject<RELAX>::clear();
00774 }
00775
00776
00777 template <class RELAX, class VECTOR>
00778 template <class MATRIX2>
00779 inline void
00780 SmootherRelaxation<RELAX, VECTOR>::initialize (
00781 const MGLevelObject<MATRIX2>& m,
00782 const typename RELAX::AdditionalData& data)
00783 {
00784 const unsigned int min = m.get_minlevel();
00785 const unsigned int max = m.get_maxlevel();
00786
00787 this->resize(min, max);
00788
00789 for (unsigned int i=min;i<=max;++i)
00790 (*this)[i].initialize(m[i], data);
00791 }
00792
00793
00794 template <class RELAX, class VECTOR>
00795 template <class MATRIX2, class DATA>
00796 inline void
00797 SmootherRelaxation<RELAX, VECTOR>::initialize (
00798 const MGLevelObject<MATRIX2>& m,
00799 const MGLevelObject<DATA>& data)
00800 {
00801 const unsigned int min = std::max(m.get_minlevel(), data.get_minlevel());
00802 const unsigned int max = std::min(m.get_maxlevel(), data.get_maxlevel());
00803
00804 this->resize(min, max);
00805
00806 for (unsigned int i=min;i<=max;++i)
00807 (*this)[i].initialize(m[i], data[i]);
00808 }
00809
00810
00811 template <class RELAX, class VECTOR>
00812 inline void
00813 SmootherRelaxation<RELAX, VECTOR>::smooth(
00814 const unsigned int level,
00815 VECTOR& u,
00816 const VECTOR& rhs) const
00817 {
00818 unsigned int maxlevel = this->get_maxlevel();
00819 unsigned int steps2 = this->steps;
00820
00821 if (this->variable)
00822 steps2 *= (1<<(maxlevel-level));
00823
00824 bool T = this->transpose;
00825 if (this->symmetric && (steps2 % 2 == 0))
00826 T = false;
00827 if (this->debug > 0)
00828 deallog << 'S' << level << ' ';
00829
00830 for (unsigned int i=0; i<steps2; ++i)
00831 {
00832 if (T)
00833 (*this)[level].Tstep(u, rhs);
00834 else
00835 (*this)[level].step(u, rhs);
00836 if (this->symmetric)
00837 T = !T;
00838 }
00839 }
00840
00841
00842 template <class RELAX, class VECTOR>
00843 inline
00844 std::size_t
00845 SmootherRelaxation<RELAX, VECTOR>::
00846 memory_consumption () const
00847 {
00848 return sizeof(*this)
00849 -sizeof(MGLevelObject<RELAX>)
00850 + MGLevelObject<RELAX>::memory_consumption()
00851 + this->mem->memory_consumption();
00852 }
00853 }
00854
00855
00856
00857
00858 template <class MATRIX, class RELAX, class VECTOR>
00859 inline
00860 MGSmootherRelaxation<MATRIX, RELAX, VECTOR>::MGSmootherRelaxation(
00861 VectorMemory<VECTOR>& mem,
00862 const unsigned int steps,
00863 const bool variable,
00864 const bool symmetric,
00865 const bool transpose)
00866 :
00867 MGSmoother<VECTOR>(mem, steps, variable, symmetric, transpose)
00868 {}
00869
00870
00871 template <class MATRIX, class RELAX, class VECTOR>
00872 inline void
00873 MGSmootherRelaxation<MATRIX, RELAX, VECTOR>::clear ()
00874 {
00875 smoothers.clear();
00876
00877 unsigned int i=matrices.get_minlevel(),
00878 max_level=matrices.get_maxlevel();
00879 for (; i<=max_level; ++i)
00880 matrices[i]=0;
00881 }
00882
00883
00884 template <class MATRIX, class RELAX, class VECTOR>
00885 template <class MATRIX2>
00886 inline void
00887 MGSmootherRelaxation<MATRIX, RELAX, VECTOR>::initialize (
00888 const MGLevelObject<MATRIX2>& m,
00889 const typename RELAX::AdditionalData& data)
00890 {
00891 const unsigned int min = m.get_minlevel();
00892 const unsigned int max = m.get_maxlevel();
00893
00894 matrices.resize(min, max);
00895 smoothers.resize(min, max);
00896
00897 for (unsigned int i=min;i<=max;++i)
00898 {
00899 matrices[i] = &m[i];
00900 smoothers[i].initialize(m[i], data);
00901 }
00902 }
00903
00904 template <class MATRIX, class RELAX, class VECTOR>
00905 template <class MATRIX2, class DATA>
00906 inline void
00907 MGSmootherRelaxation<MATRIX, RELAX, VECTOR>::initialize (
00908 const MGLevelObject<MATRIX2>& m,
00909 const MGLevelObject<DATA>& data)
00910 {
00911 const unsigned int min = m.get_minlevel();
00912 const unsigned int max = m.get_maxlevel();
00913
00914 Assert (data.get_minlevel() == min,
00915 ExcDimensionMismatch(data.get_minlevel(), min));
00916 Assert (data.get_maxlevel() == max,
00917 ExcDimensionMismatch(data.get_maxlevel(), max));
00918
00919 matrices.resize(min, max);
00920 smoothers.resize(min, max);
00921
00922 for (unsigned int i=min;i<=max;++i)
00923 {
00924 matrices[i] = &m[i];
00925 smoothers[i].initialize(m[i], data[i]);
00926 }
00927 }
00928
00929 template <class MATRIX, class RELAX, class VECTOR>
00930 template <class MATRIX2, class DATA>
00931 inline void
00932 MGSmootherRelaxation<MATRIX, RELAX, VECTOR>::initialize (
00933 const MGLevelObject<MATRIX2>& m,
00934 const DATA& data,
00935 const unsigned int row,
00936 const unsigned int col)
00937 {
00938 const unsigned int min = m.get_minlevel();
00939 const unsigned int max = m.get_maxlevel();
00940
00941 matrices.resize(min, max);
00942 smoothers.resize(min, max);
00943
00944 for (unsigned int i=min;i<=max;++i)
00945 {
00946 matrices[i] = &(m[i].block(row, col));
00947 smoothers[i].initialize(m[i].block(row, col), data);
00948 }
00949 }
00950
00951 template <class MATRIX, class RELAX, class VECTOR>
00952 template <class MATRIX2, class DATA>
00953 inline void
00954 MGSmootherRelaxation<MATRIX, RELAX, VECTOR>::initialize (
00955 const MGLevelObject<MATRIX2>& m,
00956 const MGLevelObject<DATA>& data,
00957 const unsigned int row,
00958 const unsigned int col)
00959 {
00960 const unsigned int min = m.get_minlevel();
00961 const unsigned int max = m.get_maxlevel();
00962
00963 Assert (data.get_minlevel() == min,
00964 ExcDimensionMismatch(data.get_minlevel(), min));
00965 Assert (data.get_maxlevel() == max,
00966 ExcDimensionMismatch(data.get_maxlevel(), max));
00967
00968 matrices.resize(min, max);
00969 smoothers.resize(min, max);
00970
00971 for (unsigned int i=min;i<=max;++i)
00972 {
00973 matrices[i] = &(m[i].block(row, col));
00974 smoothers[i].initialize(m[i].block(row, col), data[i]);
00975 }
00976 }
00977
00978 #ifdef DEAL_II_MULTIGRID_COMPATIBILITY
00979
00980 template <class MATRIX, class RELAX, class VECTOR>
00981 inline void
00982 MGSmootherRelaxation<MATRIX, RELAX, VECTOR>::smooth(
00983 const unsigned int level,
00984 VECTOR& u,
00985 const VECTOR& rhs) const
00986 {
00987 unsigned int maxlevel = matrices.get_maxlevel();
00988 unsigned int steps2 = this->steps;
00989
00990 if (this->variable)
00991 steps2 *= (1<<(maxlevel-level));
00992
00993 typename VectorMemory<VECTOR>::Pointer r(*this->mem);
00994 typename VectorMemory<VECTOR>::Pointer d(*this->mem);
00995 r->reinit(u);
00996 d->reinit(u);
00997
00998 bool T = this->transpose;
00999 if (this->symmetric && (steps2 % 2 == 0))
01000 T = false;
01001 if (this->debug > 0)
01002 deallog << 'S' << level << ' ';
01003
01004 for (unsigned int i=0; i<steps2; ++i)
01005 {
01006 if (T)
01007 {
01008 if (this->debug > 0)
01009 deallog << 'T';
01010 matrices[level].vmult(*r,u);
01011 r->sadd(-1.,1.,rhs);
01012 if (this->debug > 2)
01013 deallog << ' ' << r->l2_norm() << ' ';
01014 smoothers[level].Tvmult(*d, *r);
01015 if (this->debug > 1)
01016 deallog << ' ' << d->l2_norm() << ' ';
01017 } else {
01018 if (this->debug > 0)
01019 deallog << 'N';
01020 matrices[level].vmult(*r,u);
01021 r->sadd(-1.,1.,rhs);
01022 if (this->debug > 2)
01023 deallog << ' ' << r->l2_norm() << ' ';
01024 smoothers[level].vmult(*d, *r);
01025 if (this->debug > 1)
01026 deallog << ' ' << d->l2_norm() << ' ';
01027 }
01028 u += *d;
01029 if (this->symmetric)
01030 T = !T;
01031 }
01032 if (this->debug > 0)
01033 deallog << std::endl;
01034 }
01035
01036 #else
01037
01038 template <class MATRIX, class RELAX, class VECTOR>
01039 inline void
01040 MGSmootherRelaxation<MATRIX, RELAX, VECTOR>::smooth(
01041 const unsigned int level,
01042 VECTOR& u,
01043 const VECTOR& rhs) const
01044 {
01045 unsigned int maxlevel = smoothers.get_maxlevel();
01046 unsigned int steps2 = this->steps;
01047
01048 if (this->variable)
01049 steps2 *= (1<<(maxlevel-level));
01050
01051 bool T = this->transpose;
01052 if (this->symmetric && (steps2 % 2 == 0))
01053 T = false;
01054 if (this->debug > 0)
01055 deallog << 'S' << level << ' ';
01056
01057 for (unsigned int i=0; i<steps2; ++i)
01058 {
01059 if (T)
01060 smoothers[level].Tstep(u, rhs);
01061 else
01062 smoothers[level].step(u, rhs);
01063 if (this->symmetric)
01064 T = !T;
01065 }
01066 }
01067
01068 #endif
01069
01070
01071 template <class MATRIX, class RELAX, class VECTOR>
01072 inline
01073 std::size_t
01074 MGSmootherRelaxation<MATRIX, RELAX, VECTOR>::
01075 memory_consumption () const
01076 {
01077 return sizeof(*this)
01078 + matrices.memory_consumption()
01079 + smoothers.memory_consumption()
01080 + this->mem->memory_consumption();
01081 }
01082
01083
01084
01085
01086 template <class MATRIX, class PRECONDITIONER, class VECTOR>
01087 inline
01088 MGSmootherPrecondition<MATRIX, PRECONDITIONER, VECTOR>::MGSmootherPrecondition(
01089 VectorMemory<VECTOR>& mem,
01090 const unsigned int steps,
01091 const bool variable,
01092 const bool symmetric,
01093 const bool transpose)
01094 :
01095 MGSmoother<VECTOR>(mem, steps, variable, symmetric, transpose)
01096 {}
01097
01098
01099
01100 template <class MATRIX, class PRECONDITIONER, class VECTOR>
01101 inline void
01102 MGSmootherPrecondition<MATRIX, PRECONDITIONER, VECTOR>::clear ()
01103 {
01104 smoothers.clear();
01105
01106 unsigned int i=matrices.get_minlevel(),
01107 max_level=matrices.get_maxlevel();
01108 for (; i<=max_level; ++i)
01109 matrices[i]=0;
01110 }
01111
01112
01113
01114 template <class MATRIX, class PRECONDITIONER, class VECTOR>
01115 template <class MATRIX2>
01116 inline void
01117 MGSmootherPrecondition<MATRIX, PRECONDITIONER, VECTOR>::initialize (
01118 const MGLevelObject<MATRIX2>& m,
01119 const typename PRECONDITIONER::AdditionalData& data)
01120 {
01121 const unsigned int min = m.get_minlevel();
01122 const unsigned int max = m.get_maxlevel();
01123
01124 matrices.resize(min, max);
01125 smoothers.resize(min, max);
01126
01127 for (unsigned int i=min;i<=max;++i)
01128 {
01129 matrices[i] = &m[i];
01130 smoothers[i].initialize(m[i], data);
01131 }
01132 }
01133
01134
01135
01136 template <class MATRIX, class PRECONDITIONER, class VECTOR>
01137 template <class MATRIX2, class DATA>
01138 inline void
01139 MGSmootherPrecondition<MATRIX, PRECONDITIONER, VECTOR>::initialize (
01140 const MGLevelObject<MATRIX2>& m,
01141 const MGLevelObject<DATA>& data)
01142 {
01143 const unsigned int min = m.get_minlevel();
01144 const unsigned int max = m.get_maxlevel();
01145
01146 Assert (data.get_minlevel() == min,
01147 ExcDimensionMismatch(data.get_minlevel(), min));
01148 Assert (data.get_maxlevel() == max,
01149 ExcDimensionMismatch(data.get_maxlevel(), max));
01150
01151 matrices.resize(min, max);
01152 smoothers.resize(min, max);
01153
01154 for (unsigned int i=min;i<=max;++i)
01155 {
01156 matrices[i] = &m[i];
01157 smoothers[i].initialize(m[i], data[i]);
01158 }
01159 }
01160
01161
01162
01163 template <class MATRIX, class PRECONDITIONER, class VECTOR>
01164 template <class MATRIX2, class DATA>
01165 inline void
01166 MGSmootherPrecondition<MATRIX, PRECONDITIONER, VECTOR>::initialize (
01167 const MGLevelObject<MATRIX2>& m,
01168 const DATA& data,
01169 const unsigned int row,
01170 const unsigned int col)
01171 {
01172 const unsigned int min = m.get_minlevel();
01173 const unsigned int max = m.get_maxlevel();
01174
01175 matrices.resize(min, max);
01176 smoothers.resize(min, max);
01177
01178 for (unsigned int i=min;i<=max;++i)
01179 {
01180 matrices[i] = &(m[i].block(row, col));
01181 smoothers[i].initialize(m[i].block(row, col), data);
01182 }
01183 }
01184
01185
01186
01187 template <class MATRIX, class PRECONDITIONER, class VECTOR>
01188 template <class MATRIX2, class DATA>
01189 inline void
01190 MGSmootherPrecondition<MATRIX, PRECONDITIONER, VECTOR>::initialize (
01191 const MGLevelObject<MATRIX2>& m,
01192 const MGLevelObject<DATA>& data,
01193 const unsigned int row,
01194 const unsigned int col)
01195 {
01196 const unsigned int min = m.get_minlevel();
01197 const unsigned int max = m.get_maxlevel();
01198
01199 Assert (data.get_minlevel() == min,
01200 ExcDimensionMismatch(data.get_minlevel(), min));
01201 Assert (data.get_maxlevel() == max,
01202 ExcDimensionMismatch(data.get_maxlevel(), max));
01203
01204 matrices.resize(min, max);
01205 smoothers.resize(min, max);
01206
01207 for (unsigned int i=min;i<=max;++i)
01208 {
01209 matrices[i] = &(m[i].block(row, col));
01210 smoothers[i].initialize(m[i].block(row, col), data[i]);
01211 }
01212 }
01213
01214
01215
01216 template <class MATRIX, class PRECONDITIONER, class VECTOR>
01217 inline void
01218 MGSmootherPrecondition<MATRIX, PRECONDITIONER, VECTOR>::smooth(
01219 const unsigned int level,
01220 VECTOR& u,
01221 const VECTOR& rhs) const
01222 {
01223 unsigned int maxlevel = matrices.get_maxlevel();
01224 unsigned int steps2 = this->steps;
01225
01226 if (this->variable)
01227 steps2 *= (1<<(maxlevel-level));
01228
01229 typename VectorMemory<VECTOR>::Pointer r(*this->mem);
01230 typename VectorMemory<VECTOR>::Pointer d(*this->mem);
01231
01232 r->reinit(u,true);
01233 d->reinit(u,true);
01234
01235 bool T = this->transpose;
01236 if (this->symmetric && (steps2 % 2 == 0))
01237 T = false;
01238 if (this->debug > 0)
01239 deallog << 'S' << level << ' ';
01240
01241 for (unsigned int i=0; i<steps2; ++i)
01242 {
01243 if (T)
01244 {
01245 if (this->debug > 0)
01246 deallog << 'T';
01247 matrices[level].Tvmult(*r,u);
01248 r->sadd(-1.,1.,rhs);
01249 if (this->debug > 2)
01250 deallog << ' ' << r->l2_norm() << ' ';
01251 smoothers[level].Tvmult(*d, *r);
01252 if (this->debug > 1)
01253 deallog << ' ' << d->l2_norm() << ' ';
01254 } else {
01255 if (this->debug > 0)
01256 deallog << 'N';
01257 matrices[level].vmult(*r,u);
01258 r->sadd(-1.,rhs);
01259 if (this->debug > 2)
01260 deallog << ' ' << r->l2_norm() << ' ';
01261 smoothers[level].vmult(*d, *r);
01262 if (this->debug > 1)
01263 deallog << ' ' << d->l2_norm() << ' ';
01264 }
01265 u += *d;
01266 if (this->symmetric)
01267 T = !T;
01268 }
01269 if (this->debug > 0)
01270 deallog << std::endl;
01271 }
01272
01273
01274
01275 template <class MATRIX, class PRECONDITIONER, class VECTOR>
01276 inline
01277 std::size_t
01278 MGSmootherPrecondition<MATRIX, PRECONDITIONER, VECTOR>::
01279 memory_consumption () const
01280 {
01281 return sizeof(*this)
01282 + matrices.memory_consumption()
01283 + smoothers.memory_consumption()
01284 + this->mem->memory_consumption();
01285 }
01286
01287
01288 #endif // DOXYGEN
01289
01290 DEAL_II_NAMESPACE_CLOSE
01291
01292 #endif
01293