Reference documentation for deal.II version 8.1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
precondition_block.h
1 // ---------------------------------------------------------------------
2 // @f$Id: precondition_block.h 30040 2013-07-18 17:06:48Z maier @f$
3 //
4 // Copyright (C) 1999 - 2013 by the deal.II authors
5 //
6 // This file is part of the deal.II library.
7 //
8 // The deal.II library is free software; you can use it, redistribute
9 // it, and/or modify it under the terms of the GNU Lesser General
10 // Public License as published by the Free Software Foundation; either
11 // version 2.1 of the License, or (at your option) any later version.
12 // The full text of the license can be found in the file LICENSE at
13 // the top level of the deal.II distribution.
14 //
15 // ---------------------------------------------------------------------
16 
17 #ifndef __deal2__precondition_block_h
18 #define __deal2__precondition_block_h
19 
20 
21 #include <deal.II/base/config.h>
23 #include <deal.II/base/subscriptor.h>
24 #include <deal.II/base/smartpointer.h>
25 #include <deal.II/lac/precondition_block_base.h>
26 
27 #include <vector>
28 
29 DEAL_II_NAMESPACE_OPEN
30 
31 template<class MATRIX, typename inverse_type>
33 
86 template<class MATRIX, typename inverse_type = typename MATRIX::value_type>
88  : public virtual Subscriptor,
89  protected PreconditionBlockBase<inverse_type>
90 {
91 private:
95  typedef typename MATRIX::value_type number;
96 
101 
102 public:
107 
112  {
113  public:
121  const double relaxation = 1.,
122  const bool invert_diagonal = true,
123  const bool same_diagonal = false);
124 
128  double relaxation;
129 
134 
139 
150 
160  double threshold;
161  };
162 
163 
167  PreconditionBlock(bool store_diagonals = false);
168 
173 
187  void initialize (const MATRIX &A,
188  const AdditionalData parameters);
189 protected:
212  void initialize (const MATRIX &A,
213  const std::vector<size_type> &permutation,
214  const std::vector<size_type> &inverse_permutation,
215  const AdditionalData parameters);
216 
263  void set_permutation(const std::vector<size_type> &permutation,
264  const std::vector<size_type> &inverse_permutation);
265 
272  const std::vector<size_type> &permutation,
273  const std::vector<size_type> &inverse_permutation);
274 public:
283  void clear();
284 
288  bool empty () const;
289 
297  size_type j) const;
298 
326  void invert_diagblocks();
327 
345  template <typename number2>
346  void forward_step (
347  Vector<number2> &dst,
348  const Vector<number2> &prev,
349  const Vector<number2> &src,
350  const bool transpose_diagonal) const;
351 
369  template <typename number2>
370  void backward_step (
371  Vector<number2> &dst,
372  const Vector<number2> &prev,
373  const Vector<number2> &src,
374  const bool transpose_diagonal) const;
375 
376 
380  size_type block_size () const;
381 
389  unsigned int n_blocks() const DEAL_II_DEPRECATED;
390 
396  std::size_t memory_consumption () const;
397 
408  DeclException2 (ExcWrongBlockSize,
409  int, int,
410  << "The blocksize " << arg1
411  << " and the size of the matrix " << arg2
412  << " do not match.");
413 
417  DeclException0 (ExcInverseMatricesAlreadyExist);
418 
420 
421 protected:
428 
446  double relaxation;
447 
451  std::vector<size_type> permutation;
452 
457 
462 };
463 
464 
465 
476 template<class MATRIX, typename inverse_type = typename MATRIX::value_type>
477 class PreconditionBlockJacobi : public virtual Subscriptor,
478  private PreconditionBlock<MATRIX, inverse_type>
479 {
480 private:
484  typedef typename MATRIX::value_type number;
485 
486 public:
491 
496  {
497  private:
501  class Accessor
502  {
503  public:
511  const size_type row);
512 
518  size_type row() const;
519 
525  size_type column() const;
526 
530  inverse_type value() const;
531 
532  protected:
537 
542  size_type bs;
543 
547  size_type a_block;
548 
553 
558 
563  friend class const_iterator;
564  };
565 
566  public:
571  const size_type row);
572 
576  const_iterator &operator++ ();
577 
581  const_iterator &operator++ (int);
582 
586  const Accessor &operator* () const;
587 
591  const Accessor *operator-> () const;
592 
599  bool operator == (const const_iterator &) const;
603  bool operator != (const const_iterator &) const;
604 
614  bool operator < (const const_iterator &) const;
615 
616  private:
622  };
623 
644 
657  template <typename number2>
658  void vmult (Vector<number2> &, const Vector<number2> &) const;
659 
663  template <typename number2>
664  void Tvmult (Vector<number2> &, const Vector<number2> &) const;
677  template <typename number2>
678  void vmult_add (Vector<number2> &, const Vector<number2> &) const;
679 
684  template <typename number2>
685  void Tvmult_add (Vector<number2> &, const Vector<number2> &) const;
686 
691  template <typename number2>
692  void step (Vector<number2> &dst, const Vector<number2> &rhs) const;
693 
698  template <typename number2>
699  void Tstep (Vector<number2> &dst, const Vector<number2> &rhs) const;
700 
705  const_iterator begin () const;
706 
710  const_iterator end () const;
711 
716  const_iterator begin (const size_type r) const;
717 
721  const_iterator end (const size_type r) const;
722 
723 
724 private:
733  template <typename number2>
734  void do_vmult (Vector<number2> &,
735  const Vector<number2> &,
736  bool adding) const;
737 
738  friend class Accessor;
739  friend class const_iterator;
740 };
741 
742 
743 
780 template<class MATRIX, typename inverse_type = typename MATRIX::value_type>
781 class PreconditionBlockSOR : public virtual Subscriptor,
782  protected PreconditionBlock<MATRIX, inverse_type>
783 {
784 public:
789 
794 
798  typedef typename MATRIX::value_type number;
799 
816 
833  template <typename number2>
834  void vmult (Vector<number2> &, const Vector<number2> &) const;
835 
853  template <typename number2>
854  void vmult_add (Vector<number2> &, const Vector<number2> &) const;
855 
870  template <typename number2>
871  void Tvmult (Vector<number2> &, const Vector<number2> &) const;
872 
890  template <typename number2>
891  void Tvmult_add (Vector<number2> &, const Vector<number2> &) const;
892 
897  template <typename number2>
898  void step (Vector<number2> &dst, const Vector<number2> &rhs) const;
899 
904  template <typename number2>
905  void Tstep (Vector<number2> &dst, const Vector<number2> &rhs) const;
906 
907 protected:
912  PreconditionBlockSOR(bool store);
913 
927  template <typename number2>
928  void forward (Vector<number2> &,
929  const Vector<number2> &,
930  const bool transpose_diagonal,
931  const bool adding) const;
932 
946  template <typename number2>
947  void backward (Vector<number2> &,
948  const Vector<number2> &,
949  const bool transpose_diagonal,
950  const bool adding) const;
951 };
952 
953 
973 template<class MATRIX, typename inverse_type = typename MATRIX::value_type>
974 class PreconditionBlockSSOR : public virtual Subscriptor,
975  private PreconditionBlockSOR<MATRIX, inverse_type>
976 {
977 public:
982 
986  typedef typename MATRIX::value_type number;
987 
992 
993  // Keep AdditionalData accessible
995 
996  // The following are the
997  // functions of the base classes
998  // which we want to keep
999  // accessible.
1016 
1029  template <typename number2>
1030  void vmult (Vector<number2> &, const Vector<number2> &) const;
1031 
1035  template <typename number2>
1036  void Tvmult (Vector<number2> &, const Vector<number2> &) const;
1037 
1042  template <typename number2>
1043  void step (Vector<number2> &dst, const Vector<number2> &rhs) const;
1044 
1049  template <typename number2>
1050  void Tstep (Vector<number2> &dst, const Vector<number2> &rhs) const;
1051 };
1052 
1054 //---------------------------------------------------------------------------
1055 
1056 #ifndef DOXYGEN
1057 
1058 template<class MATRIX, typename inverse_type>
1059 inline bool
1061 {
1062  if (A == 0)
1063  return true;
1064  return A->empty();
1065 }
1066 
1067 
1068 template<class MATRIX, typename inverse_type>
1069 inline unsigned int
1071 {
1072  return this->size();
1073 }
1074 
1075 
1076 template<class MATRIX, typename inverse_type>
1077 inline inverse_type
1079  size_type i,
1080  size_type j) const
1081 {
1082  const size_type bs = blocksize;
1083  const unsigned int nb = i/bs;
1084 
1085  const FullMatrix<inverse_type> &B = this->inverse(nb);
1086 
1087  const size_type ib = i % bs;
1088  const size_type jb = j % bs;
1089 
1090  if (jb + nb*bs != j)
1091  {
1092  return 0.;
1093  }
1094 
1095  return B(ib, jb);
1096 }
1097 
1098 //---------------------------------------------------------------------------
1099 
1100 template<class MATRIX, typename inverse_type>
1101 inline
1104  const size_type row)
1105  :
1106  matrix(matrix),
1107  b_iterator(&matrix->inverse(0), 0, 0),
1108  b_end(&matrix->inverse(0), 0, 0)
1109 {
1110  bs = matrix->block_size();
1111  a_block = row / bs;
1112 
1113  // This is the end accessor, which
1114  // does not hava a valid block.
1115  if (a_block == matrix->size())
1116  return;
1117 
1118  const size_type r = row % bs;
1119 
1120  b_iterator = matrix->inverse(a_block).begin(r);
1121  b_end = matrix->inverse(a_block).end();
1122 
1123  Assert (a_block < matrix->size(),
1124  ExcIndexRange(a_block, 0, matrix->size()));
1125 }
1126 
1127 
1128 template<class MATRIX, typename inverse_type>
1129 inline
1132 {
1133  Assert (a_block < matrix->size(),
1134  ExcIteratorPastEnd());
1135 
1136  return bs * a_block + b_iterator->row();
1137 }
1138 
1139 
1140 template<class MATRIX, typename inverse_type>
1141 inline
1144 {
1145  Assert (a_block < matrix->size(),
1146  ExcIteratorPastEnd());
1147 
1148  return bs * a_block + b_iterator->column();
1149 }
1150 
1151 
1152 template<class MATRIX, typename inverse_type>
1153 inline
1156 {
1157  Assert (a_block < matrix->size(),
1158  ExcIteratorPastEnd());
1159 
1160  return b_iterator->value();
1161 }
1162 
1163 
1164 template<class MATRIX, typename inverse_type>
1165 inline
1168  const size_type row)
1169  :
1170  accessor(matrix, row)
1171 {}
1172 
1173 
1174 template<class MATRIX, typename inverse_type>
1175 inline
1178 {
1179  Assert (*this != accessor.matrix->end(), ExcIteratorPastEnd());
1180 
1181  ++accessor.b_iterator;
1183  {
1184  ++accessor.a_block;
1185 
1186  if (accessor.a_block < accessor.matrix->size())
1187  {
1188  accessor.b_iterator = accessor.matrix->inverse(accessor.a_block).begin();
1189  accessor.b_end = accessor.matrix->inverse(accessor.a_block).end();
1190  }
1191  }
1192  return *this;
1193 }
1194 
1195 
1196 template<class MATRIX, typename inverse_type>
1197 inline
1200 {
1201  return accessor;
1202 }
1203 
1204 
1205 template<class MATRIX, typename inverse_type>
1206 inline
1209 {
1210  return &accessor;
1211 }
1212 
1213 
1214 template<class MATRIX, typename inverse_type>
1215 inline
1216 bool
1218 operator == (const const_iterator &other) const
1219 {
1220  if (accessor.a_block == accessor.matrix->size() &&
1221  accessor.a_block == other.accessor.a_block)
1222  return true;
1223 
1224  if (accessor.a_block != other.accessor.a_block)
1225  return false;
1226 
1227  return (accessor.row() == other.accessor.row() &&
1228  accessor.column() == other.accessor.column());
1229 }
1230 
1231 
1232 template<class MATRIX, typename inverse_type>
1233 inline
1234 bool
1236 operator != (const const_iterator &other) const
1237 {
1238  return ! (*this == other);
1239 }
1240 
1241 
1242 template<class MATRIX, typename inverse_type>
1243 inline
1244 bool
1246 operator < (const const_iterator &other) const
1247 {
1248  return (accessor.row() < other.accessor.row() ||
1249  (accessor.row() == other.accessor.row() &&
1250  accessor.column() < other.accessor.column()));
1251 }
1252 
1253 
1254 template<class MATRIX, typename inverse_type>
1255 inline
1258 {
1259  return const_iterator(this, 0);
1260 }
1261 
1262 
1263 template<class MATRIX, typename inverse_type>
1264 inline
1267 {
1268  return const_iterator(this, this->size() * this->block_size());
1269 }
1270 
1271 
1272 template<class MATRIX, typename inverse_type>
1273 inline
1276  const size_type r) const
1277 {
1278  Assert (r < this->A->m(), ExcIndexRange(r, 0, this->A->m()));
1279  return const_iterator(this, r);
1280 }
1281 
1282 
1283 
1284 template<class MATRIX, typename inverse_type>
1285 inline
1288  const size_type r) const
1289 {
1290  Assert (r < this->A->m(), ExcIndexRange(r, 0, this->A->m()));
1291  return const_iterator(this, r+1);
1292 }
1293 
1294 #endif // DOXYGEN
1295 
1296 DEAL_II_NAMESPACE_CLOSE
1297 
1298 #endif
types::global_dof_index size_type
FullMatrix< inverse_type >::const_iterator b_end
types::global_dof_index size_type
unsigned int n_blocks() const DEAL_II_DEPRECATED
void backward_step(Vector< number2 > &dst, const Vector< number2 > &prev, const Vector< number2 > &src, const bool transpose_diagonal) const
bool operator<(const const_iterator &) const
bool operator==(const const_iterator &) const
MATRIX::value_type number
AdditionalData(const size_type block_size, const double relaxation=1., const bool invert_diagonal=true, const bool same_diagonal=false)
std::vector< size_type > inverse_permutation
void initialize(const MATRIX &A, const AdditionalData parameters)
MATRIX::value_type number
void forward_step(Vector< number2 > &dst, const Vector< number2 > &prev, const Vector< number2 > &src, const bool transpose_diagonal) const
std::size_t memory_consumption() const
FullMatrix< inverse_type > & inverse(size_type i)
const_iterator end() const
unsigned int global_dof_index
Definition: types.h:100
const Accessor * operator->() const
#define Assert(cond, exc)
Definition: exceptions.h:299
const PreconditionBlockJacobi< MATRIX, inverse_type > * matrix
DeclException0(ExcInverseMatricesAlreadyExist)
Accessor(const PreconditionBlockJacobi< MATRIX, inverse_type > *matrix, const size_type row)
void invert_permuted_diagblocks(const std::vector< size_type > &permutation, const std::vector< size_type > &inverse_permutation)
::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
const Accessor & operator*() const
const_iterator(const PreconditionBlockJacobi< MATRIX, inverse_type > *matrix, const size_type row)
BlockCompressedSparsityPattern CompressedBlockSparsityPattern DEAL_II_DEPRECATED
const_iterator begin() const
MATRIX::value_type number
const_iterator begin() const
bool empty() const
void set_permutation(const std::vector< size_type > &permutation, const std::vector< size_type > &inverse_permutation)
PreconditionBlockBase< inverse_type >::Inversion inversion
::ExceptionBase & ExcIteratorPastEnd()
DeclException2(ExcWrongBlockSize, int, int,<< "The blocksize "<< arg1<< " and the size of the matrix "<< arg2<< " do not match.")
std::vector< size_type > permutation
PreconditionBlock(bool store_diagonals=false)
value_type el(size_type i, size_type j) const
MATRIX::value_type number
FullMatrix< inverse_type >::const_iterator b_iterator
bool operator!=(const const_iterator &) const
SmartPointer< const MATRIX, PreconditionBlock< MATRIX, inverse_type > > A
const_iterator end() const
types::global_dof_index size_type
types::global_dof_index size_type