00001
00002
00003
00004
00005
00006
00007
00008
00009
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
00962 using typename PreconditionBlockSOR<MATRIX,inverse_type>::AdditionalData;
00963
00964
00965
00966
00967
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
01082
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