Reference documentation for deal.II version Git 24fafe1087 2021-11-29 17:11:51 +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 &
479  operator*() const;
480 
484  const Accessor *
485  operator->() const;
486 
490  bool
491  operator==(const const_iterator &) const;
495  bool
496  operator!=(const const_iterator &) const;
497 
502  bool
503  operator<(const const_iterator &) const;
504 
505  private:
510  };
511 
528 
536  template <typename number2>
537  void
538  vmult(Vector<number2> &, const Vector<number2> &) const;
539 
543  template <typename number2>
544  void
545  Tvmult(Vector<number2> &, const Vector<number2> &) const;
553  template <typename number2>
554  void
555  vmult_add(Vector<number2> &, const Vector<number2> &) const;
556 
560  template <typename number2>
561  void
562  Tvmult_add(Vector<number2> &, const Vector<number2> &) const;
563 
567  template <typename number2>
568  void
569  step(Vector<number2> &dst, const Vector<number2> &rhs) const;
570 
574  template <typename number2>
575  void
576  Tstep(Vector<number2> &dst, const Vector<number2> &rhs) const;
577 
582  begin() const;
583 
588  end() const;
589 
594  begin(const size_type r) const;
595 
600  end(const size_type r) const;
601 
602 
603 private:
610  template <typename number2>
611  void
612  do_vmult(Vector<number2> &, const Vector<number2> &, bool adding) const;
613 
614  friend class Accessor;
615  friend class const_iterator;
616 };
617 
618 
619 
653 template <typename MatrixType,
654  typename inverse_type = typename MatrixType::value_type>
656  : public virtual Subscriptor,
657  protected PreconditionBlock<MatrixType, inverse_type>
658 {
659 public:
664 
669 
673  using number = typename MatrixType::value_type;
674 
688 
699  template <typename number2>
700  void
701  vmult(Vector<number2> &, const Vector<number2> &) const;
702 
713  template <typename number2>
714  void
715  vmult_add(Vector<number2> &, const Vector<number2> &) const;
716 
725  template <typename number2>
726  void
727  Tvmult(Vector<number2> &, const Vector<number2> &) const;
728 
739  template <typename number2>
740  void
741  Tvmult_add(Vector<number2> &, const Vector<number2> &) const;
742 
746  template <typename number2>
747  void
748  step(Vector<number2> &dst, const Vector<number2> &rhs) const;
749 
753  template <typename number2>
754  void
755  Tstep(Vector<number2> &dst, const Vector<number2> &rhs) const;
756 
757 protected:
761  PreconditionBlockSOR(bool store);
762 
772  template <typename number2>
773  void
774  forward(Vector<number2> &,
775  const Vector<number2> &,
776  const bool transpose_diagonal,
777  const bool adding) const;
778 
788  template <typename number2>
789  void
790  backward(Vector<number2> &,
791  const Vector<number2> &,
792  const bool transpose_diagonal,
793  const bool adding) const;
794 };
795 
796 
816 template <typename MatrixType,
817  typename inverse_type = typename MatrixType::value_type>
819  : public virtual Subscriptor,
820  private PreconditionBlockSOR<MatrixType, inverse_type>
821 {
822 public:
827 
831  using number = typename MatrixType::value_type;
832 
837 
838  // Keep AdditionalData accessible
840 
841  // The following are the
842  // functions of the base classes
843  // which we want to keep
844  // accessible.
859 
867  template <typename number2>
868  void
869  vmult(Vector<number2> &, const Vector<number2> &) const;
870 
874  template <typename number2>
875  void
876  Tvmult(Vector<number2> &, const Vector<number2> &) const;
877 
881  template <typename number2>
882  void
883  step(Vector<number2> &dst, const Vector<number2> &rhs) const;
884 
888  template <typename number2>
889  void
890  Tstep(Vector<number2> &dst, const Vector<number2> &rhs) const;
891 };
892 
894 //---------------------------------------------------------------------------
895 
896 #ifndef DOXYGEN
897 
898 template <typename MatrixType, typename inverse_type>
899 inline bool
901 {
902  if (A == nullptr)
903  return true;
904  return A->empty();
905 }
906 
907 
908 template <typename MatrixType, typename inverse_type>
909 inline inverse_type
911 {
912  const size_type bs = blocksize;
913  const unsigned int nb = i / bs;
914 
915  const FullMatrix<inverse_type> &B = this->inverse(nb);
916 
917  const size_type ib = i % bs;
918  const size_type jb = j % bs;
919 
920  if (jb + nb * bs != j)
921  {
922  return 0.;
923  }
924 
925  return B(ib, jb);
926 }
927 
928 //---------------------------------------------------------------------------
929 
930 template <typename MatrixType, typename inverse_type>
934  const size_type row)
935  : matrix(matrix)
936  , bs(matrix->block_size())
937  , a_block(row / bs)
938  , b_iterator(&matrix->inverse(0), 0, 0)
939  , b_end(&matrix->inverse(0), 0, 0)
940 {
941  // This is the end accessor, which
942  // does not have a valid block.
943  if (a_block == matrix->size())
944  return;
945 
946  const size_type r = row % bs;
947 
948  b_iterator = matrix->inverse(a_block).begin(r);
949  b_end = matrix->inverse(a_block).end();
950 
951  AssertIndexRange(a_block, matrix->size());
952 }
953 
954 
955 template <typename MatrixType, typename inverse_type>
957 PreconditionBlockJacobi<MatrixType,
958  inverse_type>::const_iterator::Accessor::row() const
959 {
960  Assert(a_block < matrix->size(), ExcIteratorPastEnd());
961 
962  return bs * a_block + b_iterator->row();
963 }
964 
965 
966 template <typename MatrixType, typename inverse_type>
968 PreconditionBlockJacobi<MatrixType,
969  inverse_type>::const_iterator::Accessor::column() const
970 {
971  Assert(a_block < matrix->size(), ExcIteratorPastEnd());
972 
973  return bs * a_block + b_iterator->column();
974 }
975 
976 
977 template <typename MatrixType, typename inverse_type>
978 inline inverse_type
979 PreconditionBlockJacobi<MatrixType,
980  inverse_type>::const_iterator::Accessor::value() const
981 {
982  Assert(a_block < matrix->size(), ExcIteratorPastEnd());
983 
984  return b_iterator->value();
985 }
986 
987 
988 template <typename MatrixType, typename inverse_type>
992  const size_type row)
993  : accessor(matrix, row)
994 {}
995 
996 
997 template <typename MatrixType, typename inverse_type>
998 inline typename PreconditionBlockJacobi<MatrixType,
999  inverse_type>::const_iterator &
1001 {
1002  Assert(*this != accessor.matrix->end(), ExcIteratorPastEnd());
1003 
1004  ++accessor.b_iterator;
1005  if (accessor.b_iterator == accessor.b_end)
1006  {
1007  ++accessor.a_block;
1008 
1009  if (accessor.a_block < accessor.matrix->size())
1010  {
1011  accessor.b_iterator =
1012  accessor.matrix->inverse(accessor.a_block).begin();
1013  accessor.b_end = accessor.matrix->inverse(accessor.a_block).end();
1014  }
1015  }
1016  return *this;
1017 }
1018 
1019 
1020 template <typename MatrixType, typename inverse_type>
1022  const_iterator::Accessor &
1024  const
1025 {
1026  return accessor;
1027 }
1028 
1029 
1030 template <typename MatrixType, typename inverse_type>
1032  const_iterator::Accessor *
1033  PreconditionBlockJacobi<MatrixType,
1034  inverse_type>::const_iterator::operator->() const
1035 {
1036  return &accessor;
1037 }
1038 
1039 
1040 template <typename MatrixType, typename inverse_type>
1041 inline bool
1043  const const_iterator &other) const
1044 {
1045  if (accessor.a_block == accessor.matrix->size() &&
1046  accessor.a_block == other.accessor.a_block)
1047  return true;
1048 
1049  if (accessor.a_block != other.accessor.a_block)
1050  return false;
1051 
1052  return (accessor.row() == other.accessor.row() &&
1053  accessor.column() == other.accessor.column());
1054 }
1055 
1056 
1057 template <typename MatrixType, typename inverse_type>
1058 inline bool
1060  const const_iterator &other) const
1061 {
1062  return !(*this == other);
1063 }
1064 
1065 
1066 template <typename MatrixType, typename inverse_type>
1067 inline bool
1069  const const_iterator &other) const
1070 {
1071  return (accessor.row() < other.accessor.row() ||
1072  (accessor.row() == other.accessor.row() &&
1073  accessor.column() < other.accessor.column()));
1074 }
1075 
1076 
1077 template <typename MatrixType, typename inverse_type>
1078 inline
1081 {
1082  return const_iterator(this, 0);
1083 }
1084 
1085 
1086 template <typename MatrixType, typename inverse_type>
1087 inline
1090 {
1091  return const_iterator(this, this->size() * this->block_size());
1092 }
1093 
1094 
1095 template <typename MatrixType, typename inverse_type>
1096 inline
1099  const size_type r) const
1100 {
1101  AssertIndexRange(r, this->A->m());
1102  return const_iterator(this, r);
1103 }
1104 
1105 
1106 
1107 template <typename MatrixType, typename inverse_type>
1108 inline
1111  const size_type r) const
1112 {
1113  AssertIndexRange(r, this->A->m());
1114  return const_iterator(this, r + 1);
1115 }
1116 
1117 #endif // DOXYGEN
1118 
1120 
1121 #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:532
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)
bool operator!=(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1720
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:1461
value_type el(size_type i, size_type j) const
#define DeclException0(Exception0)
Definition: exceptions.h:464
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:401
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:400
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