Reference documentation for deal.II version 9.2.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
precondition_block.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1999 - 2020 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 
25 
27 
28 #include <vector>
29 
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 
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 
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  , bs(matrix->block_size())
943  , a_block(row / bs)
944  , b_iterator(&matrix->inverse(0), 0, 0)
945  , b_end(&matrix->inverse(0), 0, 0)
946 {
947  // This is the end accessor, which
948  // does not have a valid block.
949  if (a_block == matrix->size())
950  return;
951 
952  const size_type r = row % bs;
953 
954  b_iterator = matrix->inverse(a_block).begin(r);
955  b_end = matrix->inverse(a_block).end();
956 
957  AssertIndexRange(a_block, matrix->size());
958 }
959 
960 
961 template <typename MatrixType, typename inverse_type>
963 PreconditionBlockJacobi<MatrixType,
964  inverse_type>::const_iterator::Accessor::row() const
965 {
966  Assert(a_block < matrix->size(), ExcIteratorPastEnd());
967 
968  return bs * a_block + b_iterator->row();
969 }
970 
971 
972 template <typename MatrixType, typename inverse_type>
974 PreconditionBlockJacobi<MatrixType,
975  inverse_type>::const_iterator::Accessor::column() const
976 {
977  Assert(a_block < matrix->size(), ExcIteratorPastEnd());
978 
979  return bs * a_block + b_iterator->column();
980 }
981 
982 
983 template <typename MatrixType, typename inverse_type>
984 inline inverse_type
985 PreconditionBlockJacobi<MatrixType,
987 {
988  Assert(a_block < matrix->size(), ExcIteratorPastEnd());
989 
990  return b_iterator->value();
991 }
992 
993 
994 template <typename MatrixType, typename inverse_type>
998  const size_type row)
999  : accessor(matrix, row)
1000 {}
1001 
1002 
1003 template <typename MatrixType, typename inverse_type>
1004 inline
1007  operator++()
1008 {
1009  Assert(*this != accessor.matrix->end(), ExcIteratorPastEnd());
1010 
1011  ++accessor.b_iterator;
1012  if (accessor.b_iterator == accessor.b_end)
1013  {
1014  ++accessor.a_block;
1015 
1016  if (accessor.a_block < accessor.matrix->size())
1017  {
1018  accessor.b_iterator =
1019  accessor.matrix->inverse(accessor.a_block).begin();
1020  accessor.b_end = accessor.matrix->inverse(accessor.a_block).end();
1021  }
1022  }
1023  return *this;
1024 }
1025 
1026 
1027 template <typename MatrixType, typename inverse_type>
1031  operator*() const
1032 {
1033  return accessor;
1034 }
1035 
1036 
1037 template <typename MatrixType, typename inverse_type>
1041  operator->() const
1042 {
1043  return &accessor;
1044 }
1045 
1046 
1047 template <typename MatrixType, typename inverse_type>
1048 inline bool
1050 operator==(const const_iterator &other) const
1051 {
1052  if (accessor.a_block == accessor.matrix->size() &&
1053  accessor.a_block == other.accessor.a_block)
1054  return true;
1055 
1056  if (accessor.a_block != other.accessor.a_block)
1057  return false;
1058 
1059  return (accessor.row() == other.accessor.row() &&
1060  accessor.column() == other.accessor.column());
1061 }
1062 
1063 
1064 template <typename MatrixType, typename inverse_type>
1065 inline bool
1067 operator!=(const const_iterator &other) const
1068 {
1069  return !(*this == other);
1070 }
1071 
1072 
1073 template <typename MatrixType, typename inverse_type>
1074 inline bool
1076 operator<(const const_iterator &other) const
1077 {
1078  return (accessor.row() < other.accessor.row() ||
1079  (accessor.row() == other.accessor.row() &&
1080  accessor.column() < other.accessor.column()));
1081 }
1082 
1083 
1084 template <typename MatrixType, typename inverse_type>
1085 inline
1088 {
1089  return const_iterator(this, 0);
1090 }
1091 
1092 
1093 template <typename MatrixType, typename inverse_type>
1094 inline
1097 {
1098  return const_iterator(this, this->size() * this->block_size());
1099 }
1100 
1101 
1102 template <typename MatrixType, typename inverse_type>
1103 inline
1106  const size_type r) const
1107 {
1108  AssertIndexRange(r, this->A->m());
1109  return const_iterator(this, r);
1110 }
1111 
1112 
1113 
1114 template <typename MatrixType, typename inverse_type>
1115 inline
1118  const size_type r) const
1119 {
1120  AssertIndexRange(r, this->A->m());
1121  return const_iterator(this, r + 1);
1122 }
1123 
1124 #endif // DOXYGEN
1125 
1127 
1128 #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
Contents is actually a matrix.
FullMatrix< inverse_type >::const_iterator b_iterator
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
bool operator!=(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1649
const_iterator end() const
Accessor(const PreconditionBlockJacobi< MatrixType, inverse_type > *matrix, const size_type row)
bool operator<(const SynchronousIterators< Iterators > &a, const SynchronousIterators< Iterators > &b)
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)
std::vector< size_type > inverse_permutation
bool operator==(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
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
#define Assert(cond, exc)
Definition: exceptions.h:1419
value_type el(size_type i, size_type j) const
#define DeclException0(Exception0)
Definition: exceptions.h:473
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:359
VectorType::value_type * end(VectorType &V)
PreconditionBlockBase< inverse_type >::Inversion inversion
iterator end(const size_type r)
SynchronousIterators< Iterators > operator++(SynchronousIterators< Iterators > &a)
const_iterator begin() const
PreconditionBlock(bool store_diagonals=false)
MatrixTableIterators::Accessor< TransposeTable< T >, Constness, MatrixTableIterators::Storage::column_major > Accessor
Definition: table.h:1907
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:76
void invert_diagblocks()
std::vector< size_type > permutation
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:358
VectorType::value_type * begin(VectorType &V)
bool empty() const
const Accessor & operator*() const
void invert_permuted_diagblocks()
std::enable_if< std::is_floating_point< T >::value &&std::is_floating_point< U >::value, typename ProductType< std::complex< T >, std::complex< U > >::type >::type operator*(const std::complex< T > &left, const std::complex< U > &right)
const PreconditionBlockJacobi< MatrixType, inverse_type > * matrix
iterator begin(const size_type r)
static const bool value
bool operator<(const const_iterator &) const
~PreconditionBlock() override=default