Reference documentation for deal.II version Git 292c2606a1 2021-03-03 00:48:14 +0100
\(\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 
80 template <typename MatrixType,
81  typename inverse_type = typename MatrixType::value_type>
82 class PreconditionBlock : public virtual Subscriptor,
83  protected PreconditionBlockBase<inverse_type>
84 {
85 private:
89  using number = typename MatrixType::value_type;
90 
95 
96 public:
101 
106  {
107  public:
113  const double relaxation = 1.,
114  const bool invert_diagonal = true,
115  const bool same_diagonal = false);
116 
120  double relaxation;
121 
126 
131 
140 
146  double threshold;
147  };
148 
149 
153  PreconditionBlock(bool store_diagonals = false);
154 
158  ~PreconditionBlock() override = default;
159 
167  void
168  initialize(const MatrixType &A, const AdditionalData parameters);
169 
170 protected:
182  void
183  initialize(const MatrixType & A,
184  const std::vector<size_type> &permutation,
185  const std::vector<size_type> &inverse_permutation,
186  const AdditionalData parameters);
187 
211  void
212  set_permutation(const std::vector<size_type> &permutation,
213  const std::vector<size_type> &inverse_permutation);
214 
218  void
220 
221 public:
227  void
228  clear();
229 
233  bool
234  empty() const;
235 
240  value_type
241  el(size_type i, size_type j) const;
242 
258  void
260 
272  template <typename number2>
273  void
274  forward_step(Vector<number2> & dst,
275  const Vector<number2> &prev,
276  const Vector<number2> &src,
277  const bool transpose_diagonal) const;
278 
290  template <typename number2>
291  void
292  backward_step(Vector<number2> & dst,
293  const Vector<number2> &prev,
294  const Vector<number2> &src,
295  const bool transpose_diagonal) const;
296 
297 
301  size_type
302  block_size() const;
303 
308  std::size_t
309  memory_consumption() const;
310 
321  int,
322  int,
323  << "The blocksize " << arg1 << " and the size of the matrix "
324  << arg2 << " do not match.");
325 
330 
332 
333 protected:
339 
350  double relaxation;
351 
355  std::vector<size_type> permutation;
356 
360  std::vector<size_type> inverse_permutation;
361 };
362 
363 
364 
376 template <typename MatrixType,
377  typename inverse_type = typename MatrixType::value_type>
379  : public virtual Subscriptor,
380  private PreconditionBlock<MatrixType, inverse_type>
381 {
382 private:
386  using number = typename MatrixType::value_type;
387 
388 public:
393 
398  {
399  private:
403  class Accessor
404  {
405  public:
411  const size_type row);
412 
416  size_type
417  row() const;
418 
422  size_type
423  column() const;
424 
429  value() const;
430 
431  protected:
436 
441 
446 
451 
456 
457  // Make enclosing class a friend.
458  friend class const_iterator;
459  };
460 
461  public:
467  const size_type row);
468 
473  operator++();
474 
478  const Accessor &operator*() const;
479 
483  const Accessor *operator->() const;
484 
488  bool
489  operator==(const const_iterator &) const;
493  bool
494  operator!=(const const_iterator &) const;
495 
500  bool
501  operator<(const const_iterator &) const;
502 
503  private:
508  };
509 
526 
534  template <typename number2>
535  void
536  vmult(Vector<number2> &, const Vector<number2> &) const;
537 
541  template <typename number2>
542  void
543  Tvmult(Vector<number2> &, const Vector<number2> &) const;
551  template <typename number2>
552  void
553  vmult_add(Vector<number2> &, const Vector<number2> &) const;
554 
558  template <typename number2>
559  void
560  Tvmult_add(Vector<number2> &, const Vector<number2> &) const;
561 
565  template <typename number2>
566  void
567  step(Vector<number2> &dst, const Vector<number2> &rhs) const;
568 
572  template <typename number2>
573  void
574  Tstep(Vector<number2> &dst, const Vector<number2> &rhs) const;
575 
580  begin() const;
581 
586  end() const;
587 
592  begin(const size_type r) const;
593 
598  end(const size_type r) const;
599 
600 
601 private:
608  template <typename number2>
609  void
610  do_vmult(Vector<number2> &, const Vector<number2> &, bool adding) const;
611 
612  friend class Accessor;
613  friend class const_iterator;
614 };
615 
616 
617 
651 template <typename MatrixType,
652  typename inverse_type = typename MatrixType::value_type>
654  : public virtual Subscriptor,
655  protected PreconditionBlock<MatrixType, inverse_type>
656 {
657 public:
662 
667 
671  using number = typename MatrixType::value_type;
672 
686 
697  template <typename number2>
698  void
699  vmult(Vector<number2> &, const Vector<number2> &) const;
700 
711  template <typename number2>
712  void
713  vmult_add(Vector<number2> &, const Vector<number2> &) const;
714 
723  template <typename number2>
724  void
725  Tvmult(Vector<number2> &, const Vector<number2> &) const;
726 
737  template <typename number2>
738  void
739  Tvmult_add(Vector<number2> &, const Vector<number2> &) const;
740 
744  template <typename number2>
745  void
746  step(Vector<number2> &dst, const Vector<number2> &rhs) const;
747 
751  template <typename number2>
752  void
753  Tstep(Vector<number2> &dst, const Vector<number2> &rhs) const;
754 
755 protected:
759  PreconditionBlockSOR(bool store);
760 
770  template <typename number2>
771  void
772  forward(Vector<number2> &,
773  const Vector<number2> &,
774  const bool transpose_diagonal,
775  const bool adding) const;
776 
786  template <typename number2>
787  void
788  backward(Vector<number2> &,
789  const Vector<number2> &,
790  const bool transpose_diagonal,
791  const bool adding) const;
792 };
793 
794 
814 template <typename MatrixType,
815  typename inverse_type = typename MatrixType::value_type>
817  : public virtual Subscriptor,
818  private PreconditionBlockSOR<MatrixType, inverse_type>
819 {
820 public:
825 
829  using number = typename MatrixType::value_type;
830 
835 
836  // Keep AdditionalData accessible
838 
839  // The following are the
840  // functions of the base classes
841  // which we want to keep
842  // accessible.
857 
865  template <typename number2>
866  void
867  vmult(Vector<number2> &, const Vector<number2> &) const;
868 
872  template <typename number2>
873  void
874  Tvmult(Vector<number2> &, const Vector<number2> &) const;
875 
879  template <typename number2>
880  void
881  step(Vector<number2> &dst, const Vector<number2> &rhs) const;
882 
886  template <typename number2>
887  void
888  Tstep(Vector<number2> &dst, const Vector<number2> &rhs) const;
889 };
890 
892 //---------------------------------------------------------------------------
893 
894 #ifndef DOXYGEN
895 
896 template <typename MatrixType, typename inverse_type>
897 inline bool
899 {
900  if (A == nullptr)
901  return true;
902  return A->empty();
903 }
904 
905 
906 template <typename MatrixType, typename inverse_type>
907 inline inverse_type
909 {
910  const size_type bs = blocksize;
911  const unsigned int nb = i / bs;
912 
913  const FullMatrix<inverse_type> &B = this->inverse(nb);
914 
915  const size_type ib = i % bs;
916  const size_type jb = j % bs;
917 
918  if (jb + nb * bs != j)
919  {
920  return 0.;
921  }
922 
923  return B(ib, jb);
924 }
925 
926 //---------------------------------------------------------------------------
927 
928 template <typename MatrixType, typename inverse_type>
932  const size_type row)
933  : matrix(matrix)
934  , bs(matrix->block_size())
935  , a_block(row / bs)
936  , b_iterator(&matrix->inverse(0), 0, 0)
937  , b_end(&matrix->inverse(0), 0, 0)
938 {
939  // This is the end accessor, which
940  // does not have a valid block.
941  if (a_block == matrix->size())
942  return;
943 
944  const size_type r = row % bs;
945 
946  b_iterator = matrix->inverse(a_block).begin(r);
947  b_end = matrix->inverse(a_block).end();
948 
949  AssertIndexRange(a_block, matrix->size());
950 }
951 
952 
953 template <typename MatrixType, typename inverse_type>
955 PreconditionBlockJacobi<MatrixType,
956  inverse_type>::const_iterator::Accessor::row() const
957 {
958  Assert(a_block < matrix->size(), ExcIteratorPastEnd());
959 
960  return bs * a_block + b_iterator->row();
961 }
962 
963 
964 template <typename MatrixType, typename inverse_type>
966 PreconditionBlockJacobi<MatrixType,
967  inverse_type>::const_iterator::Accessor::column() const
968 {
969  Assert(a_block < matrix->size(), ExcIteratorPastEnd());
970 
971  return bs * a_block + b_iterator->column();
972 }
973 
974 
975 template <typename MatrixType, typename inverse_type>
976 inline inverse_type
977 PreconditionBlockJacobi<MatrixType,
978  inverse_type>::const_iterator::Accessor::value() const
979 {
980  Assert(a_block < matrix->size(), ExcIteratorPastEnd());
981 
982  return b_iterator->value();
983 }
984 
985 
986 template <typename MatrixType, typename inverse_type>
990  const size_type row)
991  : accessor(matrix, row)
992 {}
993 
994 
995 template <typename MatrixType, typename inverse_type>
996 inline
999  operator++()
1000 {
1001  Assert(*this != accessor.matrix->end(), ExcIteratorPastEnd());
1002 
1003  ++accessor.b_iterator;
1004  if (accessor.b_iterator == accessor.b_end)
1005  {
1006  ++accessor.a_block;
1007 
1008  if (accessor.a_block < accessor.matrix->size())
1009  {
1010  accessor.b_iterator =
1011  accessor.matrix->inverse(accessor.a_block).begin();
1012  accessor.b_end = accessor.matrix->inverse(accessor.a_block).end();
1013  }
1014  }
1015  return *this;
1016 }
1017 
1018 
1019 template <typename MatrixType, typename inverse_type>
1021  const_iterator::Accessor &
1023  operator*() const
1024 {
1025  return accessor;
1026 }
1027 
1028 
1029 template <typename MatrixType, typename inverse_type>
1031  const_iterator::Accessor *
1033  operator->() const
1034 {
1035  return &accessor;
1036 }
1037 
1038 
1039 template <typename MatrixType, typename inverse_type>
1040 inline bool
1042 operator==(const const_iterator &other) const
1043 {
1044  if (accessor.a_block == accessor.matrix->size() &&
1045  accessor.a_block == other.accessor.a_block)
1046  return true;
1047 
1048  if (accessor.a_block != other.accessor.a_block)
1049  return false;
1050 
1051  return (accessor.row() == other.accessor.row() &&
1052  accessor.column() == other.accessor.column());
1053 }
1054 
1055 
1056 template <typename MatrixType, typename inverse_type>
1057 inline bool
1059 operator!=(const const_iterator &other) const
1060 {
1061  return !(*this == other);
1062 }
1063 
1064 
1065 template <typename MatrixType, typename inverse_type>
1066 inline bool
1068 operator<(const const_iterator &other) const
1069 {
1070  return (accessor.row() < other.accessor.row() ||
1071  (accessor.row() == other.accessor.row() &&
1072  accessor.column() < other.accessor.column()));
1073 }
1074 
1075 
1076 template <typename MatrixType, typename inverse_type>
1077 inline
1080 {
1081  return const_iterator(this, 0);
1082 }
1083 
1084 
1085 template <typename MatrixType, typename inverse_type>
1086 inline
1089 {
1090  return const_iterator(this, this->size() * this->block_size());
1091 }
1092 
1093 
1094 template <typename MatrixType, typename inverse_type>
1095 inline
1098  const size_type r) const
1099 {
1100  AssertIndexRange(r, this->A->m());
1101  return const_iterator(this, r);
1102 }
1103 
1104 
1105 
1106 template <typename MatrixType, typename inverse_type>
1107 inline
1110  const size_type r) const
1111 {
1112  AssertIndexRange(r, this->A->m());
1113  return const_iterator(this, r + 1);
1114 }
1115 
1116 #endif // DOXYGEN
1117 
1119 
1120 #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:538
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:1691
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:1466
value_type el(size_type i, size_type j) const
#define DeclException0(Exception0)
Definition: exceptions.h:470
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:394
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)
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:393
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)
bool operator<(const const_iterator &) const
~PreconditionBlock() override=default