include/deal.II/multigrid/mg_smoother.h

00001 //---------------------------------------------------------------------------
00002 //    @f$Id: mg_smoother.h 25460 2012-04-27 09:19:06Z bangerth @f$
00003 //
00004 //    Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 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__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  * MGSmootherBase is defined in mg_base.h
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                                       * Use the transpose of the
00145                                       * relaxation method instead of
00146                                       * the method itself. This has no
00147                                       * effect if #symmetric smoothing
00148                                       * is chosen.
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 //     template <class MATRIX2>
00292 //     void initialize (const MGLevelObject<MatrixBlock<MATRIX2> >& matrices,
00293 //                   const typename RELAX::AdditionalData & additional_data = typename RELAX::AdditionalData());
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 /* ------------------------------- Inline functions -------------------------- */
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 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

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