Reference documentation for deal.II version Git 7026f387cc 2019-10-15 14:19:01 -0400
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
precondition_block.h
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1999 - 2018 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii_precondition_block_h
17 #define dealii_precondition_block_h
18 
19 
20 #include <deal.II/base/config.h>
21 
22 #include <deal.II/base/exceptions.h>
23 #include <deal.II/base/smartpointer.h>
24 #include <deal.II/base/subscriptor.h>
25 
26 #include <deal.II/lac/precondition_block_base.h>
27 
28 #include <vector>
29 
30 DEAL_II_NAMESPACE_OPEN
31 
82 template <typename MatrixType,
83  typename inverse_type = typename MatrixType::value_type>
84 class PreconditionBlock : public virtual Subscriptor,
85  protected PreconditionBlockBase<inverse_type>
86 {
87 private:
91  using number = typename MatrixType::value_type;
92 
96  using value_type = inverse_type;
97 
98 public:
103 
108  {
109  public:
115  const double relaxation = 1.,
116  const bool invert_diagonal = true,
117  const bool same_diagonal = false);
118 
122  double relaxation;
123 
128 
133 
142 
148  double threshold;
149  };
150 
151 
155  PreconditionBlock(bool store_diagonals = false);
156 
160  ~PreconditionBlock() override = default;
161 
169  void
170  initialize(const MatrixType &A, const AdditionalData parameters);
171 
172 protected:
184  void
185  initialize(const MatrixType & A,
186  const std::vector<size_type> &permutation,
187  const std::vector<size_type> &inverse_permutation,
188  const AdditionalData parameters);
189 
213  void
214  set_permutation(const std::vector<size_type> &permutation,
215  const std::vector<size_type> &inverse_permutation);
216 
220  void
222 
223 public:
229  void
230  clear();
231 
235  bool
236  empty() const;
237 
242  value_type
243  el(size_type i, size_type j) const;
244 
260  void
262 
274  template <typename number2>
275  void
276  forward_step(Vector<number2> & dst,
277  const Vector<number2> &prev,
278  const Vector<number2> &src,
279  const bool transpose_diagonal) const;
280 
292  template <typename number2>
293  void
294  backward_step(Vector<number2> & dst,
295  const Vector<number2> &prev,
296  const Vector<number2> &src,
297  const bool transpose_diagonal) const;
298 
299 
303  size_type
304  block_size() const;
305 
310  std::size_t
311  memory_consumption() const;
312 
323  int,
324  int,
325  << "The blocksize " << arg1 << " and the size of the matrix "
326  << arg2 << " do not match.");
327 
332 
334 
335 protected:
341 
352  double relaxation;
353 
357  std::vector<size_type> permutation;
358 
362  std::vector<size_type> inverse_permutation;
363 };
364 
365 
366 
380 template <typename MatrixType,
381  typename inverse_type = typename MatrixType::value_type>
383  : public virtual Subscriptor,
384  private PreconditionBlock<MatrixType, inverse_type>
385 {
386 private:
390  using number = typename MatrixType::value_type;
391 
392 public:
397 
402  {
403  private:
407  class Accessor
408  {
409  public:
415  const size_type row);
416 
420  size_type
421  row() const;
422 
426  size_type
427  column() const;
428 
432  inverse_type
433  value() const;
434 
435  protected:
440 
445 
450 
455 
460 
461  // Make enclosing class a friend.
462  friend class const_iterator;
463  };
464 
465  public:
471  const size_type row);
472 
477  operator++();
478 
482  const Accessor &operator*() const;
483 
487  const Accessor *operator->() const;
488 
492  bool
493  operator==(const const_iterator &) const;
497  bool
498  operator!=(const const_iterator &) const;
499 
504  bool
505  operator<(const const_iterator &) const;
506 
507  private:
512  };
513 
530 
538  template <typename number2>
539  void
540  vmult(Vector<number2> &, const Vector<number2> &) const;
541 
545  template <typename number2>
546  void
547  Tvmult(Vector<number2> &, const Vector<number2> &) const;
555  template <typename number2>
556  void
557  vmult_add(Vector<number2> &, const Vector<number2> &) const;
558 
562  template <typename number2>
563  void
564  Tvmult_add(Vector<number2> &, const Vector<number2> &) const;
565 
569  template <typename number2>
570  void
571  step(Vector<number2> &dst, const Vector<number2> &rhs) const;
572 
576  template <typename number2>
577  void
578  Tstep(Vector<number2> &dst, const Vector<number2> &rhs) const;
579 
584  begin() const;
585 
590  end() const;
591 
596  begin(const size_type r) const;
597 
602  end(const size_type r) const;
603 
604 
605 private:
612  template <typename number2>
613  void
614  do_vmult(Vector<number2> &, const Vector<number2> &, bool adding) const;
615 
616  friend class Accessor;
617  friend class const_iterator;
618 };
619 
620 
621 
657 template <typename MatrixType,
658  typename inverse_type = typename MatrixType::value_type>
660  : public virtual Subscriptor,
661  protected PreconditionBlock<MatrixType, inverse_type>
662 {
663 public:
668 
673 
677  using number = typename MatrixType::value_type;
678 
692 
703  template <typename number2>
704  void
705  vmult(Vector<number2> &, const Vector<number2> &) const;
706 
717  template <typename number2>
718  void
719  vmult_add(Vector<number2> &, const Vector<number2> &) const;
720 
729  template <typename number2>
730  void
731  Tvmult(Vector<number2> &, const Vector<number2> &) const;
732 
743  template <typename number2>
744  void
745  Tvmult_add(Vector<number2> &, const Vector<number2> &) const;
746 
750  template <typename number2>
751  void
752  step(Vector<number2> &dst, const Vector<number2> &rhs) const;
753 
757  template <typename number2>
758  void
759  Tstep(Vector<number2> &dst, const Vector<number2> &rhs) const;
760 
761 protected:
765  PreconditionBlockSOR(bool store);
766 
776  template <typename number2>
777  void
778  forward(Vector<number2> &,
779  const Vector<number2> &,
780  const bool transpose_diagonal,
781  const bool adding) const;
782 
792  template <typename number2>
793  void
794  backward(Vector<number2> &,
795  const Vector<number2> &,
796  const bool transpose_diagonal,
797  const bool adding) const;
798 };
799 
800 
822 template <typename MatrixType,
823  typename inverse_type = typename MatrixType::value_type>
825  : public virtual Subscriptor,
826  private PreconditionBlockSOR<MatrixType, inverse_type>
827 {
828 public:
833 
837  using number = typename MatrixType::value_type;
838 
843 
844  // Keep AdditionalData accessible
846 
847  // The following are the
848  // functions of the base classes
849  // which we want to keep
850  // accessible.
865 
873  template <typename number2>
874  void
875  vmult(Vector<number2> &, const Vector<number2> &) const;
876 
880  template <typename number2>
881  void
882  Tvmult(Vector<number2> &, const Vector<number2> &) const;
883 
887  template <typename number2>
888  void
889  step(Vector<number2> &dst, const Vector<number2> &rhs) const;
890 
894  template <typename number2>
895  void
896  Tstep(Vector<number2> &dst, const Vector<number2> &rhs) const;
897 };
898 
900 //---------------------------------------------------------------------------
901 
902 #ifndef DOXYGEN
903 
904 template <typename MatrixType, typename inverse_type>
905 inline bool
907 {
908  if (A == nullptr)
909  return true;
910  return A->empty();
911 }
912 
913 
914 template <typename MatrixType, typename inverse_type>
915 inline inverse_type
917 {
918  const size_type bs = blocksize;
919  const unsigned int nb = i / bs;
920 
921  const FullMatrix<inverse_type> &B = this->inverse(nb);
922 
923  const size_type ib = i % bs;
924  const size_type jb = j % bs;
925 
926  if (jb + nb * bs != j)
927  {
928  return 0.;
929  }
930 
931  return B(ib, jb);
932 }
933 
934 //---------------------------------------------------------------------------
935 
936 template <typename MatrixType, typename inverse_type>
940  const size_type row)
941  : matrix(matrix)
942  , b_iterator(&matrix->inverse(0), 0, 0)
943  , b_end(&matrix->inverse(0), 0, 0)
944 {
945  bs = matrix->block_size();
946  a_block = row / bs;
947 
948  // This is the end accessor, which
949  // does not have a valid block.
950  if (a_block == matrix->size())
951  return;
952 
953  const size_type r = row % bs;
954 
955  b_iterator = matrix->inverse(a_block).begin(r);
956  b_end = matrix->inverse(a_block).end();
957 
958  Assert(a_block < matrix->size(), ExcIndexRange(a_block, 0, matrix->size()));
959 }
960 
961 
962 template <typename MatrixType, typename inverse_type>
964 PreconditionBlockJacobi<MatrixType,
965  inverse_type>::const_iterator::Accessor::row() const
966 {
967  Assert(a_block < matrix->size(), ExcIteratorPastEnd());
968 
969  return bs * a_block + b_iterator->row();
970 }
971 
972 
973 template <typename MatrixType, typename inverse_type>
975 PreconditionBlockJacobi<MatrixType,
976  inverse_type>::const_iterator::Accessor::column() const
977 {
978  Assert(a_block < matrix->size(), ExcIteratorPastEnd());
979 
980  return bs * a_block + b_iterator->column();
981 }
982 
983 
984 template <typename MatrixType, typename inverse_type>
985 inline inverse_type
986 PreconditionBlockJacobi<MatrixType,
987  inverse_type>::const_iterator::Accessor::value() const
988 {
989  Assert(a_block < matrix->size(), ExcIteratorPastEnd());
990 
991  return b_iterator->value();
992 }
993 
994 
995 template <typename MatrixType, typename inverse_type>
999  const size_type row)
1000  : accessor(matrix, row)
1001 {}
1002 
1003 
1004 template <typename MatrixType, typename inverse_type>
1005 inline
1008  operator++()
1009 {
1010  Assert(*this != accessor.matrix->end(), ExcIteratorPastEnd());
1011 
1012  ++accessor.b_iterator;
1013  if (accessor.b_iterator == accessor.b_end)
1014  {
1015  ++accessor.a_block;
1016 
1017  if (accessor.a_block < accessor.matrix->size())
1018  {
1019  accessor.b_iterator =
1020  accessor.matrix->inverse(accessor.a_block).begin();
1021  accessor.b_end = accessor.matrix->inverse(accessor.a_block).end();
1022  }
1023  }
1024  return *this;
1025 }
1026 
1027 
1028 template <typename MatrixType, typename inverse_type>
1030  const_iterator::Accessor &
1032  operator*() const
1033 {
1034  return accessor;
1035 }
1036 
1037 
1038 template <typename MatrixType, typename inverse_type>
1040  const_iterator::Accessor *
1042  operator->() const
1043 {
1044  return &accessor;
1045 }
1046 
1047 
1048 template <typename MatrixType, typename inverse_type>
1049 inline bool
1051 operator==(const const_iterator &other) const
1052 {
1053  if (accessor.a_block == accessor.matrix->size() &&
1054  accessor.a_block == other.accessor.a_block)
1055  return true;
1056 
1057  if (accessor.a_block != other.accessor.a_block)
1058  return false;
1059 
1060  return (accessor.row() == other.accessor.row() &&
1061  accessor.column() == other.accessor.column());
1062 }
1063 
1064 
1065 template <typename MatrixType, typename inverse_type>
1066 inline bool
1068 operator!=(const const_iterator &other) const
1069 {
1070  return !(*this == other);
1071 }
1072 
1073 
1074 template <typename MatrixType, typename inverse_type>
1075 inline bool
1077 operator<(const const_iterator &other) const
1078 {
1079  return (accessor.row() < other.accessor.row() ||
1080  (accessor.row() == other.accessor.row() &&
1081  accessor.column() < other.accessor.column()));
1082 }
1083 
1084 
1085 template <typename MatrixType, typename inverse_type>
1086 inline
1089 {
1090  return const_iterator(this, 0);
1091 }
1092 
1093 
1094 template <typename MatrixType, typename inverse_type>
1095 inline
1098 {
1099  return const_iterator(this, this->size() * this->block_size());
1100 }
1101 
1102 
1103 template <typename MatrixType, typename inverse_type>
1104 inline
1107  const size_type r) const
1108 {
1109  Assert(r < this->A->m(), ExcIndexRange(r, 0, this->A->m()));
1110  return const_iterator(this, r);
1111 }
1112 
1113 
1114 
1115 template <typename MatrixType, typename inverse_type>
1116 inline
1119  const size_type r) const
1120 {
1121  Assert(r < this->A->m(), ExcIndexRange(r, 0, this->A->m()));
1122  return const_iterator(this, r + 1);
1123 }
1124 
1125 #endif // DOXYGEN
1126 
1127 DEAL_II_NAMESPACE_CLOSE
1128 
1129 #endif
FullMatrix< inverse_type >::const_iterator b_end
std::size_t memory_consumption() const
void set_permutation(const std::vector< size_type > &permutation, const std::vector< size_type > &inverse_permutation)
SmartPointer< const MatrixType, PreconditionBlock< MatrixType, inverse_type > > A
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:541
FullMatrix< inverse_type >::const_iterator b_iterator
DEAL_II_CONSTEXPR SymmetricTensor< rank_, dim, Number > operator*(const SymmetricTensor< rank_, dim, Number > &t, const Number &factor)
AdditionalData(const size_type block_size, const double relaxation=1., const bool invert_diagonal=true, const bool same_diagonal=false)
static ::ExceptionBase & ExcWrongBlockSize(int arg1, int arg2)
const Accessor * operator->() const
const_iterator end() const
Accessor(const PreconditionBlockJacobi< MatrixType, inverse_type > *matrix, const size_type row)
void backward_step(Vector< number2 > &dst, const Vector< number2 > &prev, const Vector< number2 > &src, const bool transpose_diagonal) const
void initialize(const MatrixType &A, const AdditionalData parameters)
static ::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
std::vector< size_type > inverse_permutation
const_iterator(const PreconditionBlockJacobi< MatrixType, inverse_type > *matrix, const size_type row)
size_type block_size() const
typename MatrixType::value_type number
static ::ExceptionBase & ExcInverseMatricesAlreadyExist()
FullMatrix< inverse_type > & inverse(size_type i)
bool operator!=(const const_iterator &) const
bool operator==(const const_iterator &) const
types::global_dof_index size_type
#define Assert(cond, exc)
Definition: exceptions.h:1407
value_type el(size_type i, size_type j) const
#define DeclException0(Exception0)
Definition: exceptions.h:473
types::global_dof_index size_type
PreconditionBlockBase< inverse_type >::Inversion inversion
iterator end(const size_type r)
const_iterator begin() const
PreconditionBlock(bool store_diagonals=false)
void forward_step(Vector< number2 > &dst, const Vector< number2 > &prev, const Vector< number2 > &src, const bool transpose_diagonal) const
static ::ExceptionBase & ExcIteratorPastEnd()
unsigned int global_dof_index
Definition: types.h:89
void invert_diagblocks()
std::vector< size_type > permutation
bool empty() const
const Accessor & operator*() const
void invert_permuted_diagblocks()
const PreconditionBlockJacobi< MatrixType, inverse_type > * matrix
iterator begin(const size_type r)
bool operator<(const const_iterator &) const
~PreconditionBlock() override=default
inverse_type value_type