include/deal.II/lac/precondition_block.h

00001 //---------------------------------------------------------------------------
00002 //    @f$Id: precondition_block.h 25378 2012-04-03 20:24:22Z kanschat @f$
00003 //
00004 //    Copyright (C) 1999, 2000, 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__precondition_block_h
00013 #define __deal2__precondition_block_h
00014 
00015 
00016 #include <deal.II/base/config.h>
00017 #include <deal.II/base/exceptions.h>
00018 #include <deal.II/base/subscriptor.h>
00019 #include <deal.II/base/smartpointer.h>
00020 #include <deal.II/lac/precondition_block_base.h>
00021 
00022 #include <vector>
00023 
00024 DEAL_II_NAMESPACE_OPEN
00025 
00026 template<class MATRIX, typename inverse_type>
00027 class PreconditionBlockJacobi;
00028 
00074 template<class MATRIX, typename inverse_type = typename MATRIX::value_type>
00075 class PreconditionBlock
00076   : public virtual Subscriptor,
00077     protected PreconditionBlockBase<inverse_type>
00078 {
00079   private:
00083     typedef typename MATRIX::value_type number;
00084 
00088     typedef inverse_type value_type;
00089 
00090   public:
00094     class AdditionalData
00095     {
00096       public:
00103         AdditionalData (const unsigned int block_size,
00104                         const double relaxation = 1.,
00105                         const bool invert_diagonal = true,
00106                         const bool same_diagonal = false);
00107 
00111         double relaxation;
00112 
00116         unsigned int block_size;
00117 
00121         bool invert_diagonal;
00122 
00127         bool same_diagonal;
00132         typename PreconditionBlockBase<inverse_type>::Inversion inversion;
00133 
00143         double threshold;
00144     };
00145 
00146 
00150     PreconditionBlock(bool store_diagonals = false);
00151 
00155     ~PreconditionBlock();
00156 
00170     void initialize (const MATRIX& A,
00171                      const AdditionalData parameters);
00172   protected:
00195     void initialize (const MATRIX& A,
00196                      const std::vector<unsigned int>& permutation,
00197                      const std::vector<unsigned int>& inverse_permutation,
00198                      const AdditionalData parameters);
00199 
00246     void set_permutation(const std::vector<unsigned int>& permutation,
00247                          const std::vector<unsigned int>& inverse_permutation);
00248 
00254     void invert_permuted_diagblocks(
00255       const std::vector<unsigned int>& permutation,
00256       const std::vector<unsigned int>& inverse_permutation);
00257   public:
00266     void clear();
00267 
00271     bool empty () const;
00272 
00279     value_type el(unsigned int i,
00280                   unsigned int j) const;
00281 
00309     void invert_diagblocks();
00310 
00328     template <typename number2>
00329     void forward_step (
00330       Vector<number2>       &dst,
00331       const Vector<number2> &prev,
00332       const Vector<number2> &src,
00333       const bool transpose_diagonal) const;
00334 
00352     template <typename number2>
00353     void backward_step (
00354       Vector<number2>       &dst,
00355       const Vector<number2> &prev,
00356       const Vector<number2> &src,
00357       const bool transpose_diagonal) const;
00358 
00359 
00363     unsigned int block_size () const;
00364 
00372     unsigned int n_blocks() const;
00373 
00379     std::size_t memory_consumption () const;
00380 
00391     DeclException2 (ExcWrongBlockSize,
00392                     int, int,
00393                     << "The blocksize " << arg1
00394                     << " and the size of the matrix " << arg2
00395                     << " do not match.");
00396 
00400     DeclException0 (ExcInverseMatricesAlreadyExist);
00401 
00403 
00404   protected:
00410     unsigned int blocksize;
00411 
00424     SmartPointer<const MATRIX,PreconditionBlock<MATRIX,inverse_type> > A;
00429     double relaxation;
00430 
00434     std::vector<unsigned int> permutation;
00435 
00439     std::vector<unsigned int> inverse_permutation;
00440 
00445 };
00446 
00447 
00448 
00459 template<class MATRIX, typename inverse_type = typename MATRIX::value_type>
00460 class PreconditionBlockJacobi : public virtual Subscriptor,
00461                                 private PreconditionBlock<MATRIX, inverse_type>
00462 {
00463   private:
00467     typedef typename MATRIX::value_type number;
00468 
00469   public:
00473     class const_iterator
00474     {
00475       private:
00479         class Accessor
00480         {
00481           public:
00488             Accessor (const PreconditionBlockJacobi<MATRIX, inverse_type> *matrix,
00489                       const unsigned int row);
00490 
00496             unsigned int row() const;
00497 
00503             unsigned int column() const;
00504 
00508             inverse_type value() const;
00509 
00510           protected:
00514             const PreconditionBlockJacobi<MATRIX, inverse_type>* matrix;
00515 
00520             unsigned int bs;
00521 
00525             unsigned int a_block;
00526 
00530             typename FullMatrix<inverse_type>::const_iterator b_iterator;
00531 
00535             typename FullMatrix<inverse_type>::const_iterator b_end;
00536 
00541             friend class const_iterator;
00542         };
00543 
00544       public:
00548         const_iterator(const PreconditionBlockJacobi<MATRIX, inverse_type>* matrix,
00549                        const unsigned int row);
00550 
00554         const_iterator& operator++ ();
00555 
00559         const_iterator& operator++ (int);
00560 
00564         const Accessor& operator* () const;
00565 
00569         const Accessor* operator-> () const;
00570 
00577         bool operator == (const const_iterator&) const;
00581         bool operator != (const const_iterator&) const;
00582 
00592         bool operator < (const const_iterator&) const;
00593 
00594       private:
00599         Accessor accessor;
00600     };
00601 
00605     using PreconditionBlock<MATRIX, inverse_type>::initialize;
00606     using PreconditionBlock<MATRIX, inverse_type>::clear;
00607     using PreconditionBlock<MATRIX, inverse_type>::empty;
00608     using PreconditionBlock<MATRIX, inverse_type>::el;
00609     using PreconditionBlock<MATRIX, inverse_type>::set_same_diagonal;
00610     using PreconditionBlock<MATRIX, inverse_type>::invert_diagblocks;
00611     using PreconditionBlock<MATRIX, inverse_type>::block_size;
00612     using PreconditionBlockBase<inverse_type>::size;
00613     using PreconditionBlockBase<inverse_type>::inverse;
00614     using PreconditionBlockBase<inverse_type>::inverse_householder;
00615     using PreconditionBlockBase<inverse_type>::inverse_svd;
00616     using PreconditionBlockBase<inverse_type>::log_statistics;
00620     using PreconditionBlock<MATRIX, inverse_type>::n_blocks;
00621     using PreconditionBlock<MATRIX, inverse_type>::set_permutation;
00622 
00635     template <typename number2>
00636     void vmult (Vector<number2>&, const Vector<number2>&) const;
00637 
00641     template <typename number2>
00642     void Tvmult (Vector<number2>&, const Vector<number2>&) const;
00655     template <typename number2>
00656     void vmult_add (Vector<number2>&, const Vector<number2>&) const;
00657 
00662     template <typename number2>
00663     void Tvmult_add (Vector<number2>&, const Vector<number2>&) const;
00664 
00669     template <typename number2>
00670     void step (Vector<number2>& dst, const Vector<number2>& rhs) const;
00671 
00676     template <typename number2>
00677     void Tstep (Vector<number2>& dst, const Vector<number2>& rhs) const;
00678 
00683     const_iterator begin () const;
00684 
00688     const_iterator end () const;
00689 
00694     const_iterator begin (const unsigned int r) const;
00695 
00699     const_iterator end (const unsigned int r) const;
00700 
00701 
00702   private:
00711     template <typename number2>
00712     void do_vmult (Vector<number2>&,
00713                    const Vector<number2>&,
00714                    bool adding) const;
00715 
00716     friend class Accessor;
00717     friend class const_iterator;
00718 };
00719 
00720 
00721 
00758 template<class MATRIX, typename inverse_type = typename MATRIX::value_type>
00759 class PreconditionBlockSOR : public virtual Subscriptor,
00760                              protected PreconditionBlock<MATRIX, inverse_type>
00761 {
00762   public:
00766     PreconditionBlockSOR();
00767 
00771     typedef typename MATRIX::value_type number;
00772 
00776     using typename PreconditionBlock<MATRIX,inverse_type>::AdditionalData;
00777     using PreconditionBlock<MATRIX, inverse_type>::initialize;
00778     using PreconditionBlock<MATRIX, inverse_type>::clear;
00779     using PreconditionBlock<MATRIX, inverse_type>::empty;
00780     using PreconditionBlockBase<inverse_type>::size;
00781     using PreconditionBlockBase<inverse_type>::inverse;
00782     using PreconditionBlockBase<inverse_type>::inverse_householder;
00783     using PreconditionBlockBase<inverse_type>::inverse_svd;
00784     using PreconditionBlock<MATRIX, inverse_type>::el;
00785     using PreconditionBlock<MATRIX, inverse_type>::set_same_diagonal;
00786     using PreconditionBlock<MATRIX, inverse_type>::invert_diagblocks;
00787     using PreconditionBlock<MATRIX, inverse_type>::set_permutation;
00788     using PreconditionBlockBase<inverse_type>::log_statistics;
00789 
00806     template <typename number2>
00807     void vmult (Vector<number2>&, const Vector<number2>&) const;
00808 
00826     template <typename number2>
00827     void vmult_add (Vector<number2>&, const Vector<number2>&) const;
00828 
00843     template <typename number2>
00844     void Tvmult (Vector<number2>&, const Vector<number2>&) const;
00845 
00863     template <typename number2>
00864     void Tvmult_add (Vector<number2>&, const Vector<number2>&) const;
00865 
00870     template <typename number2>
00871     void step (Vector<number2>& dst, const Vector<number2>& rhs) const;
00872 
00877     template <typename number2>
00878     void Tstep (Vector<number2>& dst, const Vector<number2>& rhs) const;
00879 
00880   protected:
00885     PreconditionBlockSOR(bool store);
00886 
00900     template <typename number2>
00901     void forward (Vector<number2>&,
00902                   const Vector<number2>&,
00903                   const bool transpose_diagonal,
00904                   const bool adding) const;
00905 
00919     template <typename number2>
00920     void backward (Vector<number2>&,
00921                    const Vector<number2>&,
00922                    const bool transpose_diagonal,
00923                    const bool adding) const;
00924 };
00925 
00926 
00946 template<class MATRIX, typename inverse_type = typename MATRIX::value_type>
00947 class PreconditionBlockSSOR : public virtual Subscriptor,
00948                               private PreconditionBlockSOR<MATRIX, inverse_type>
00949 {
00950   public:
00954     typedef typename MATRIX::value_type number;
00955 
00959     PreconditionBlockSSOR ();
00960 
00961                                      // Keep AdditionalData accessible
00962     using typename PreconditionBlockSOR<MATRIX,inverse_type>::AdditionalData;
00963 
00964                                      // The following are the
00965                                      // functions of the base classes
00966                                      // which we want to keep
00967                                      // accessible.
00972     using PreconditionBlockSOR<MATRIX,inverse_type>::initialize;
00973     using PreconditionBlockSOR<MATRIX,inverse_type>::clear;
00974     using PreconditionBlockBase<inverse_type>::size;
00975     using PreconditionBlockBase<inverse_type>::inverse;
00976     using PreconditionBlockBase<inverse_type>::inverse_householder;
00977     using PreconditionBlockBase<inverse_type>::inverse_svd;
00978     using PreconditionBlockBase<inverse_type>::log_statistics;
00979     using PreconditionBlockSOR<MATRIX,inverse_type>::set_permutation;
00980     using PreconditionBlockSOR<MATRIX, inverse_type>::empty;
00981     using PreconditionBlockSOR<MATRIX, inverse_type>::el;
00982     using PreconditionBlockSOR<MATRIX,inverse_type>::set_same_diagonal;
00983     using PreconditionBlockSOR<MATRIX,inverse_type>::invert_diagblocks;
00984 
00997     template <typename number2>
00998     void vmult (Vector<number2>&, const Vector<number2>&) const;
00999 
01003     template <typename number2>
01004     void Tvmult (Vector<number2>&, const Vector<number2>&) const;
01005 
01010     template <typename number2>
01011     void step (Vector<number2>& dst, const Vector<number2>& rhs) const;
01012 
01017     template <typename number2>
01018     void Tstep (Vector<number2>& dst, const Vector<number2>& rhs) const;
01019 };
01020 
01022 //---------------------------------------------------------------------------
01023 
01024 #ifndef DOXYGEN
01025 
01026 template<class MATRIX, typename inverse_type>
01027 inline bool
01028 PreconditionBlock<MATRIX, inverse_type>::empty () const
01029 {
01030   if (A == 0)
01031     return true;
01032   return A->empty();
01033 }
01034 
01035 
01036 template<class MATRIX, typename inverse_type>
01037 inline unsigned int
01038 PreconditionBlock<MATRIX, inverse_type>::n_blocks () const
01039 {
01040   return this->size();
01041 }
01042 
01043 
01044 template<class MATRIX, typename inverse_type>
01045 inline inverse_type
01046 PreconditionBlock<MATRIX, inverse_type>::el (
01047   unsigned int i,
01048   unsigned int j) const
01049 {
01050   const unsigned int bs = blocksize;
01051   const unsigned int nb = i/bs;
01052 
01053   const FullMatrix<inverse_type>& B = this->inverse(nb);
01054 
01055   const unsigned int ib = i % bs;
01056   const unsigned int jb = j % bs;
01057 
01058   if (jb + nb*bs != j)
01059     {
01060       return 0.;
01061     }
01062 
01063   return B(ib, jb);
01064 }
01065 
01066 //---------------------------------------------------------------------------
01067 
01068 template<class MATRIX, typename inverse_type>
01069 inline
01070 PreconditionBlockJacobi<MATRIX, inverse_type>::const_iterator::Accessor::
01071 Accessor (const PreconditionBlockJacobi<MATRIX, inverse_type>* matrix,
01072           const unsigned int row)
01073                 :
01074                 matrix(matrix),
01075                 b_iterator(&matrix->inverse(0), 0, 0),
01076                 b_end(&matrix->inverse(0), 0, 0)
01077 {
01078   bs = matrix->block_size();
01079   a_block = row / bs;
01080 
01081                                    // This is the end accessor, which
01082                                    // does not hava a valid block.
01083   if (a_block == matrix->n_blocks())
01084     return;
01085 
01086   const unsigned int r = row % bs;
01087 
01088   b_iterator = matrix->inverse(a_block).begin(r);
01089   b_end = matrix->inverse(a_block).end();
01090 
01091   Assert (a_block < matrix->n_blocks(),
01092           ExcIndexRange(a_block, 0, matrix->n_blocks()));
01093 }
01094 
01095 
01096 template<class MATRIX, typename inverse_type>
01097 inline
01098 unsigned int
01099 PreconditionBlockJacobi<MATRIX, inverse_type>::const_iterator::Accessor::row() const
01100 {
01101   Assert (a_block < matrix->n_blocks(),
01102           ExcIteratorPastEnd());
01103 
01104   return bs * a_block + b_iterator->row();
01105 }
01106 
01107 
01108 template<class MATRIX, typename inverse_type>
01109 inline
01110 unsigned int
01111 PreconditionBlockJacobi<MATRIX, inverse_type>::const_iterator::Accessor::column() const
01112 {
01113   Assert (a_block < matrix->n_blocks(),
01114           ExcIteratorPastEnd());
01115 
01116   return bs * a_block + b_iterator->column();
01117 }
01118 
01119 
01120 template<class MATRIX, typename inverse_type>
01121 inline
01122 inverse_type
01123 PreconditionBlockJacobi<MATRIX, inverse_type>::const_iterator::Accessor::value() const
01124 {
01125   Assert (a_block < matrix->n_blocks(),
01126           ExcIteratorPastEnd());
01127 
01128   return b_iterator->value();
01129 }
01130 
01131 
01132 template<class MATRIX, typename inverse_type>
01133 inline
01134 PreconditionBlockJacobi<MATRIX, inverse_type>::const_iterator::
01135 const_iterator(const PreconditionBlockJacobi<MATRIX, inverse_type> *matrix,
01136                const unsigned int row)
01137                 :
01138                 accessor(matrix, row)
01139 {}
01140 
01141 
01142 template<class MATRIX, typename inverse_type>
01143 inline
01144 typename PreconditionBlockJacobi<MATRIX, inverse_type>::const_iterator &
01145 PreconditionBlockJacobi<MATRIX, inverse_type>::const_iterator::operator++ ()
01146 {
01147   Assert (*this != accessor.matrix->end(), ExcIteratorPastEnd());
01148 
01149   ++accessor.b_iterator;
01150   if (accessor.b_iterator == accessor.b_end)
01151     {
01152       ++accessor.a_block;
01153 
01154       if (accessor.a_block < accessor.matrix->n_blocks())
01155         {
01156           accessor.b_iterator = accessor.matrix->inverse(accessor.a_block).begin();
01157           accessor.b_end = accessor.matrix->inverse(accessor.a_block).end();
01158         }
01159     }
01160   return *this;
01161 }
01162 
01163 
01164 template<class MATRIX, typename inverse_type>
01165 inline
01166 const typename PreconditionBlockJacobi<MATRIX, inverse_type>::const_iterator::Accessor &
01167 PreconditionBlockJacobi<MATRIX, inverse_type>::const_iterator::operator* () const
01168 {
01169   return accessor;
01170 }
01171 
01172 
01173 template<class MATRIX, typename inverse_type>
01174 inline
01175 const typename PreconditionBlockJacobi<MATRIX, inverse_type>::const_iterator::Accessor *
01176 PreconditionBlockJacobi<MATRIX, inverse_type>::const_iterator::operator-> () const
01177 {
01178   return &accessor;
01179 }
01180 
01181 
01182 template<class MATRIX, typename inverse_type>
01183 inline
01184 bool
01185 PreconditionBlockJacobi<MATRIX, inverse_type>::const_iterator::
01186 operator == (const const_iterator& other) const
01187 {
01188   if (accessor.a_block == accessor.matrix->n_blocks() &&
01189       accessor.a_block == other.accessor.a_block)
01190     return true;
01191 
01192   if (accessor.a_block != other.accessor.a_block)
01193     return false;
01194 
01195   return (accessor.row() == other.accessor.row() &&
01196           accessor.column() == other.accessor.column());
01197 }
01198 
01199 
01200 template<class MATRIX, typename inverse_type>
01201 inline
01202 bool
01203 PreconditionBlockJacobi<MATRIX, inverse_type>::const_iterator::
01204 operator != (const const_iterator& other) const
01205 {
01206   return ! (*this == other);
01207 }
01208 
01209 
01210 template<class MATRIX, typename inverse_type>
01211 inline
01212 bool
01213 PreconditionBlockJacobi<MATRIX, inverse_type>::const_iterator::
01214 operator < (const const_iterator& other) const
01215 {
01216   return (accessor.row() < other.accessor.row() ||
01217           (accessor.row() == other.accessor.row() &&
01218            accessor.column() < other.accessor.column()));
01219 }
01220 
01221 
01222 template<class MATRIX, typename inverse_type>
01223 inline
01224 typename PreconditionBlockJacobi<MATRIX, inverse_type>::const_iterator
01225 PreconditionBlockJacobi<MATRIX, inverse_type>::begin () const
01226 {
01227   return const_iterator(this, 0);
01228 }
01229 
01230 
01231 template<class MATRIX, typename inverse_type>
01232 inline
01233 typename PreconditionBlockJacobi<MATRIX, inverse_type>::const_iterator
01234 PreconditionBlockJacobi<MATRIX, inverse_type>::end () const
01235 {
01236   return const_iterator(this, this->n_blocks() * this->block_size());
01237 }
01238 
01239 
01240 template<class MATRIX, typename inverse_type>
01241 inline
01242 typename PreconditionBlockJacobi<MATRIX, inverse_type>::const_iterator
01243 PreconditionBlockJacobi<MATRIX, inverse_type>::begin (
01244   const unsigned int r) const
01245 {
01246   Assert (r < this->A->m(), ExcIndexRange(r, 0, this->A->m()));
01247   return const_iterator(this, r);
01248 }
01249 
01250 
01251 
01252 template<class MATRIX, typename inverse_type>
01253 inline
01254 typename PreconditionBlockJacobi<MATRIX, inverse_type>::const_iterator
01255 PreconditionBlockJacobi<MATRIX, inverse_type>::end (
01256   const unsigned int r) const
01257 {
01258   Assert (r < this->A->m(), ExcIndexRange(r, 0, this->A->m()));
01259   return const_iterator(this, r+1);
01260 }
01261 
01262 #endif // DOXYGEN
01263 
01264 DEAL_II_NAMESPACE_CLOSE
01265 
01266 #endif
01267 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

deal.II documentation generated on Wed May 23 2012 06:07:39 by doxygen 1.7.3