Reference documentation for deal.II version 9.1.1
\(\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 
464  friend class const_iterator;
465  };
466 
467  public:
473  const size_type row);
474 
479  operator++();
480 
484  const Accessor &operator*() const;
485 
489  const Accessor *operator->() const;
490 
494  bool
495  operator==(const const_iterator &) const;
499  bool
500  operator!=(const const_iterator &) const;
501 
506  bool
507  operator<(const const_iterator &) const;
508 
509  private:
514  };
515 
532 
540  template <typename number2>
541  void
542  vmult(Vector<number2> &, const Vector<number2> &) const;
543 
547  template <typename number2>
548  void
549  Tvmult(Vector<number2> &, const Vector<number2> &) const;
557  template <typename number2>
558  void
559  vmult_add(Vector<number2> &, const Vector<number2> &) const;
560 
564  template <typename number2>
565  void
566  Tvmult_add(Vector<number2> &, const Vector<number2> &) const;
567 
571  template <typename number2>
572  void
573  step(Vector<number2> &dst, const Vector<number2> &rhs) const;
574 
578  template <typename number2>
579  void
580  Tstep(Vector<number2> &dst, const Vector<number2> &rhs) const;
581 
586  begin() const;
587 
592  end() const;
593 
598  begin(const size_type r) const;
599 
604  end(const size_type r) const;
605 
606 
607 private:
614  template <typename number2>
615  void
616  do_vmult(Vector<number2> &, const Vector<number2> &, bool adding) const;
617 
618  friend class Accessor;
619  friend class const_iterator;
620 };
621 
622 
623 
659 template <typename MatrixType,
660  typename inverse_type = typename MatrixType::value_type>
662  : public virtual Subscriptor,
663  protected PreconditionBlock<MatrixType, inverse_type>
664 {
665 public:
670 
675 
679  using number = typename MatrixType::value_type;
680 
694 
705  template <typename number2>
706  void
707  vmult(Vector<number2> &, const Vector<number2> &) const;
708 
719  template <typename number2>
720  void
721  vmult_add(Vector<number2> &, const Vector<number2> &) const;
722 
731  template <typename number2>
732  void
733  Tvmult(Vector<number2> &, const Vector<number2> &) const;
734 
745  template <typename number2>
746  void
747  Tvmult_add(Vector<number2> &, const Vector<number2> &) const;
748 
752  template <typename number2>
753  void
754  step(Vector<number2> &dst, const Vector<number2> &rhs) const;
755 
759  template <typename number2>
760  void
761  Tstep(Vector<number2> &dst, const Vector<number2> &rhs) const;
762 
763 protected:
767  PreconditionBlockSOR(bool store);
768 
778  template <typename number2>
779  void
780  forward(Vector<number2> &,
781  const Vector<number2> &,
782  const bool transpose_diagonal,
783  const bool adding) const;
784 
794  template <typename number2>
795  void
796  backward(Vector<number2> &,
797  const Vector<number2> &,
798  const bool transpose_diagonal,
799  const bool adding) const;
800 };
801 
802 
824 template <typename MatrixType,
825  typename inverse_type = typename MatrixType::value_type>
827  : public virtual Subscriptor,
828  private PreconditionBlockSOR<MatrixType, inverse_type>
829 {
830 public:
835 
839  using number = typename MatrixType::value_type;
840 
845 
846  // Keep AdditionalData accessible
848 
849  // The following are the
850  // functions of the base classes
851  // which we want to keep
852  // accessible.
867 
875  template <typename number2>
876  void
877  vmult(Vector<number2> &, const Vector<number2> &) const;
878 
882  template <typename number2>
883  void
884  Tvmult(Vector<number2> &, const Vector<number2> &) const;
885 
889  template <typename number2>
890  void
891  step(Vector<number2> &dst, const Vector<number2> &rhs) const;
892 
896  template <typename number2>
897  void
898  Tstep(Vector<number2> &dst, const Vector<number2> &rhs) const;
899 };
900 
902 //---------------------------------------------------------------------------
903 
904 #ifndef DOXYGEN
905 
906 template <typename MatrixType, typename inverse_type>
907 inline bool
909 {
910  if (A == nullptr)
911  return true;
912  return A->empty();
913 }
914 
915 
916 template <typename MatrixType, typename inverse_type>
917 inline inverse_type
918 PreconditionBlock<MatrixType, inverse_type>::el(size_type i, size_type j) const
919 {
920  const size_type bs = blocksize;
921  const unsigned int nb = i / bs;
922 
923  const FullMatrix<inverse_type> &B = this->inverse(nb);
924 
925  const size_type ib = i % bs;
926  const size_type jb = j % bs;
927 
928  if (jb + nb * bs != j)
929  {
930  return 0.;
931  }
932 
933  return B(ib, jb);
934 }
935 
936 //---------------------------------------------------------------------------
937 
938 template <typename MatrixType, typename inverse_type>
942  const size_type row)
943  : matrix(matrix)
944  , b_iterator(&matrix->inverse(0), 0, 0)
945  , b_end(&matrix->inverse(0), 0, 0)
946 {
947  bs = matrix->block_size();
948  a_block = row / bs;
949 
950  // This is the end accessor, which
951  // does not have a valid block.
952  if (a_block == matrix->size())
953  return;
954 
955  const size_type r = row % bs;
956 
957  b_iterator = matrix->inverse(a_block).begin(r);
958  b_end = matrix->inverse(a_block).end();
959 
960  Assert(a_block < matrix->size(), ExcIndexRange(a_block, 0, matrix->size()));
961 }
962 
963 
964 template <typename MatrixType, typename inverse_type>
966 PreconditionBlockJacobi<MatrixType,
967  inverse_type>::const_iterator::Accessor::row() const
968 {
969  Assert(a_block < matrix->size(), ExcIteratorPastEnd());
970 
971  return bs * a_block + b_iterator->row();
972 }
973 
974 
975 template <typename MatrixType, typename inverse_type>
977 PreconditionBlockJacobi<MatrixType,
978  inverse_type>::const_iterator::Accessor::column() const
979 {
980  Assert(a_block < matrix->size(), ExcIteratorPastEnd());
981 
982  return bs * a_block + b_iterator->column();
983 }
984 
985 
986 template <typename MatrixType, typename inverse_type>
987 inline inverse_type
988 PreconditionBlockJacobi<MatrixType,
989  inverse_type>::const_iterator::Accessor::value() const
990 {
991  Assert(a_block < matrix->size(), ExcIteratorPastEnd());
992 
993  return b_iterator->value();
994 }
995 
996 
997 template <typename MatrixType, typename inverse_type>
1001  const size_type row)
1002  : accessor(matrix, row)
1003 {}
1004 
1005 
1006 template <typename MatrixType, typename inverse_type>
1007 inline
1010  operator++()
1011 {
1012  Assert(*this != accessor.matrix->end(), ExcIteratorPastEnd());
1013 
1014  ++accessor.b_iterator;
1016  {
1017  ++accessor.a_block;
1018 
1019  if (accessor.a_block < accessor.matrix->size())
1020  {
1022  accessor.matrix->inverse(accessor.a_block).begin();
1023  accessor.b_end = accessor.matrix->inverse(accessor.a_block).end();
1024  }
1025  }
1026  return *this;
1027 }
1028 
1029 
1030 template <typename MatrixType, typename inverse_type>
1032  const_iterator::Accessor &
1034  operator*() const
1035 {
1036  return accessor;
1037 }
1038 
1039 
1040 template <typename MatrixType, typename inverse_type>
1042  const_iterator::Accessor *
1044  operator->() const
1045 {
1046  return &accessor;
1047 }
1048 
1049 
1050 template <typename MatrixType, typename inverse_type>
1051 inline bool
1053 operator==(const const_iterator &other) const
1054 {
1055  if (accessor.a_block == accessor.matrix->size() &&
1056  accessor.a_block == other.accessor.a_block)
1057  return true;
1058 
1059  if (accessor.a_block != other.accessor.a_block)
1060  return false;
1061 
1062  return (accessor.row() == other.accessor.row() &&
1063  accessor.column() == other.accessor.column());
1064 }
1065 
1066 
1067 template <typename MatrixType, typename inverse_type>
1068 inline bool
1070 operator!=(const const_iterator &other) const
1071 {
1072  return !(*this == other);
1073 }
1074 
1075 
1076 template <typename MatrixType, typename inverse_type>
1077 inline bool
1079 operator<(const const_iterator &other) const
1080 {
1081  return (accessor.row() < other.accessor.row() ||
1082  (accessor.row() == other.accessor.row() &&
1083  accessor.column() < other.accessor.column()));
1084 }
1085 
1086 
1087 template <typename MatrixType, typename inverse_type>
1088 inline
1091 {
1092  return const_iterator(this, 0);
1093 }
1094 
1095 
1096 template <typename MatrixType, typename inverse_type>
1097 inline
1100 {
1101  return const_iterator(this, this->size() * this->block_size());
1102 }
1103 
1104 
1105 template <typename MatrixType, typename inverse_type>
1106 inline
1109  const size_type r) const
1110 {
1111  Assert(r < this->A->m(), ExcIndexRange(r, 0, this->A->m()));
1112  return const_iterator(this, r);
1113 }
1114 
1115 
1116 
1117 template <typename MatrixType, typename inverse_type>
1118 inline
1121  const size_type r) const
1122 {
1123  Assert(r < this->A->m(), ExcIndexRange(r, 0, this->A->m()));
1124  return const_iterator(this, r + 1);
1125 }
1126 
1127 #endif // DOXYGEN
1128 
1129 DEAL_II_NAMESPACE_CLOSE
1130 
1131 #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
void Tvmult_add(Vector< number2 > &, const Vector< number2 > &) const
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:541
Contents is actually a matrix.
FullMatrix< inverse_type >::const_iterator b_iterator
void Tstep(Vector< number2 > &dst, const Vector< number2 > &rhs) const
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)
void step(Vector< number2 > &dst, const Vector< number2 > &rhs) const
const Accessor * operator->() const
void Tvmult(Vector< number2 > &, const Vector< number2 > &) 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)
void Tvmult_add(Vector< number2 > &, const Vector< number2 > &) const
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()
void step(Vector< number2 > &dst, const Vector< number2 > &rhs) const
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
void vmult(Vector< number2 > &, const Vector< number2 > &) const
void Tvmult(Vector< number2 > &, const Vector< number2 > &) const
void vmult_add(Vector< number2 > &, const Vector< number2 > &) const
#define DeclException0(Exception0)
Definition: exceptions.h:473
types::global_dof_index size_type
PreconditionBlockBase< inverse_type >::Inversion inversion
void vmult_add(Vector< number2 > &, const Vector< number2 > &) const
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
void Tvmult(Vector< number2 > &, const Vector< number2 > &) const
static ::ExceptionBase & ExcIteratorPastEnd()
void do_vmult(Vector< number2 > &, const Vector< number2 > &, bool adding) const
void backward(Vector< number2 > &, const Vector< number2 > &, const bool transpose_diagonal, const bool adding) const
unsigned int global_dof_index
Definition: types.h:89
void vmult(Vector< number2 > &, const Vector< number2 > &) const
types::global_dof_index size_type
void Tstep(Vector< number2 > &dst, const Vector< number2 > &rhs) const
void invert_diagblocks()
std::vector< size_type > permutation
void vmult(Vector< number2 > &, const Vector< number2 > &) const
bool empty() const
void Tstep(Vector< number2 > &dst, const Vector< number2 > &rhs) const
void step(Vector< number2 > &dst, const Vector< number2 > &rhs) const
const Accessor & operator*() const
void invert_permuted_diagblocks()
const PreconditionBlockJacobi< MatrixType, inverse_type > * matrix
void forward(Vector< number2 > &, const Vector< number2 > &, const bool transpose_diagonal, const bool adding) const
bool operator<(const const_iterator &) const
~PreconditionBlock() override=default
inverse_type value_type