Reference documentation for deal.II version GIT 112f7bbc69 2023-05-28 21:25:02+00:00
\(\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 - 2021 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 
81 template <typename MatrixType,
82  typename inverse_type = typename MatrixType::value_type>
83 class PreconditionBlock : public virtual Subscriptor,
84  protected PreconditionBlockBase<inverse_type>
85 {
86 private:
90  using number = typename MatrixType::value_type;
91 
95  using value_type = inverse_type;
96 
97 public:
102 
107  {
108  public:
114  const double relaxation = 1.,
115  const bool invert_diagonal = true,
116  const bool same_diagonal = false);
117 
121  double relaxation;
122 
127 
132 
141 
147  double threshold;
148  };
149 
150 
155 
159  ~PreconditionBlock() override = default;
160 
168  void
169  initialize(const MatrixType &A, const AdditionalData parameters);
170 
171 protected:
183  void
184  initialize(const MatrixType & A,
185  const std::vector<size_type> &permutation,
186  const std::vector<size_type> &inverse_permutation,
187  const AdditionalData parameters);
188 
212  void
213  set_permutation(const std::vector<size_type> &permutation,
214  const std::vector<size_type> &inverse_permutation);
215 
219  void
221 
222 public:
228  void
229  clear();
230 
234  bool
235  empty() const;
236 
241  value_type
242  el(size_type i, size_type j) const;
243 
259  void
261 
273  template <typename number2>
274  void
276  const Vector<number2> &prev,
277  const Vector<number2> &src,
278  const bool transpose_diagonal) const;
279 
291  template <typename number2>
292  void
294  const Vector<number2> &prev,
295  const Vector<number2> &src,
296  const bool transpose_diagonal) const;
297 
298 
302  size_type
303  block_size() const;
304 
309  std::size_t
311 
322  int,
323  int,
324  << "The blocksize " << arg1 << " and the size of the matrix "
325  << arg2 << " do not match.");
326 
331 
334 protected:
340 
351  double relaxation;
352 
356  std::vector<size_type> permutation;
357 
361  std::vector<size_type> inverse_permutation;
362 };
363 
364 
365 
377 template <typename MatrixType,
378  typename inverse_type = typename MatrixType::value_type>
380  : public virtual Subscriptor,
381  private PreconditionBlock<MatrixType, inverse_type>
382 {
383 private:
387  using number = typename MatrixType::value_type;
388 
389 public:
394 
399  {
400  private:
404  class Accessor
405  {
406  public:
412  const size_type row);
413 
417  size_type
418  row() const;
419 
423  size_type
424  column() const;
425 
429  inverse_type
430  value() const;
431 
432  protected:
437 
442 
447 
452 
457 
458  // Make enclosing class a friend.
459  friend class const_iterator;
460  };
461 
462  public:
468  const size_type row);
469 
475 
479  const Accessor &
480  operator*() const;
481 
485  const Accessor *
486  operator->() const;
487 
491  bool
492  operator==(const const_iterator &) const;
496  bool
497  operator!=(const const_iterator &) const;
498 
503  bool
504  operator<(const const_iterator &) const;
505 
506  private:
511  };
512 
529 
537  template <typename number2>
538  void
540 
544  template <typename number2>
545  void
554  template <typename number2>
555  void
557 
561  template <typename number2>
562  void
564 
568  template <typename number2>
569  void
570  step(Vector<number2> &dst, const Vector<number2> &rhs) const;
571 
575  template <typename number2>
576  void
577  Tstep(Vector<number2> &dst, const Vector<number2> &rhs) const;
578 
583  begin() const;
584 
589  end() const;
590 
595  begin(const size_type r) const;
596 
601  end(const size_type r) const;
602 
603 
604 private:
611  template <typename number2>
612  void
613  do_vmult(Vector<number2> &, const Vector<number2> &, bool adding) const;
614 
615  friend class Accessor;
616  friend class const_iterator;
617 };
618 
619 
620 
654 template <typename MatrixType,
655  typename inverse_type = typename MatrixType::value_type>
657  : public virtual Subscriptor,
658  protected PreconditionBlock<MatrixType, inverse_type>
659 {
660 public:
665 
670 
674  using number = typename MatrixType::value_type;
675 
689 
700  template <typename number2>
701  void
703 
714  template <typename number2>
715  void
717 
726  template <typename number2>
727  void
729 
740  template <typename number2>
741  void
743 
747  template <typename number2>
748  void
749  step(Vector<number2> &dst, const Vector<number2> &rhs) const;
750 
754  template <typename number2>
755  void
756  Tstep(Vector<number2> &dst, const Vector<number2> &rhs) const;
757 
758 protected:
762  PreconditionBlockSOR(bool store);
763 
773  template <typename number2>
774  void
776  const Vector<number2> &,
777  const bool transpose_diagonal,
778  const bool adding) const;
779 
789  template <typename number2>
790  void
792  const Vector<number2> &,
793  const bool transpose_diagonal,
794  const bool adding) const;
795 };
796 
797 
817 template <typename MatrixType,
818  typename inverse_type = typename MatrixType::value_type>
820  : public virtual Subscriptor,
821  private PreconditionBlockSOR<MatrixType, inverse_type>
822 {
823 public:
828 
832  using number = typename MatrixType::value_type;
833 
838 
839  // Keep AdditionalData accessible
841 
842  // The following are the
843  // functions of the base classes
844  // which we want to keep
845  // accessible.
860 
868  template <typename number2>
869  void
871 
875  template <typename number2>
876  void
878 
882  template <typename number2>
883  void
884  step(Vector<number2> &dst, const Vector<number2> &rhs) const;
885 
889  template <typename number2>
890  void
891  Tstep(Vector<number2> &dst, const Vector<number2> &rhs) const;
892 };
893 
895 //---------------------------------------------------------------------------
896 
897 #ifndef DOXYGEN
898 
899 template <typename MatrixType, typename inverse_type>
900 inline bool
902 {
903  if (A == nullptr)
904  return true;
905  return A->empty();
906 }
907 
908 
909 template <typename MatrixType, typename inverse_type>
910 inline inverse_type
912 {
913  const size_type bs = blocksize;
914  const unsigned int nb = i / bs;
915 
916  const FullMatrix<inverse_type> &B = this->inverse(nb);
917 
918  const size_type ib = i % bs;
919  const size_type jb = j % bs;
920 
921  if (jb + nb * bs != j)
922  {
923  return 0.;
924  }
925 
926  return B(ib, jb);
927 }
928 
929 //---------------------------------------------------------------------------
930 
931 template <typename MatrixType, typename inverse_type>
935  const size_type row)
936  : matrix(matrix)
937  , bs(matrix->block_size())
938  , a_block(row / bs)
939  , b_iterator(&matrix->inverse(0), 0, 0)
940  , b_end(&matrix->inverse(0), 0, 0)
941 {
942  // This is the end accessor, which
943  // does not have a valid block.
944  if (a_block == matrix->size())
945  return;
946 
947  const size_type r = row % bs;
948 
949  b_iterator = matrix->inverse(a_block).begin(r);
950  b_end = matrix->inverse(a_block).end();
951 
952  AssertIndexRange(a_block, matrix->size());
953 }
954 
955 
956 template <typename MatrixType, typename inverse_type>
958 PreconditionBlockJacobi<MatrixType,
959  inverse_type>::const_iterator::Accessor::row() const
960 {
961  Assert(a_block < matrix->size(), ExcIteratorPastEnd());
962 
963  return bs * a_block + b_iterator->row();
964 }
965 
966 
967 template <typename MatrixType, typename inverse_type>
969 PreconditionBlockJacobi<MatrixType,
970  inverse_type>::const_iterator::Accessor::column() const
971 {
972  Assert(a_block < matrix->size(), ExcIteratorPastEnd());
973 
974  return bs * a_block + b_iterator->column();
975 }
976 
977 
978 template <typename MatrixType, typename inverse_type>
979 inline inverse_type
980 PreconditionBlockJacobi<MatrixType,
981  inverse_type>::const_iterator::Accessor::value() const
982 {
983  Assert(a_block < matrix->size(), ExcIteratorPastEnd());
984 
985  return b_iterator->value();
986 }
987 
988 
989 template <typename MatrixType, typename inverse_type>
993  const size_type row)
994  : accessor(matrix, row)
995 {}
996 
997 
998 template <typename MatrixType, typename inverse_type>
999 inline typename PreconditionBlockJacobi<MatrixType,
1000  inverse_type>::const_iterator &
1002 {
1003  Assert(*this != accessor.matrix->end(), ExcIteratorPastEnd());
1004 
1005  ++accessor.b_iterator;
1006  if (accessor.b_iterator == accessor.b_end)
1007  {
1008  ++accessor.a_block;
1009 
1010  if (accessor.a_block < accessor.matrix->size())
1011  {
1012  accessor.b_iterator =
1013  accessor.matrix->inverse(accessor.a_block).begin();
1014  accessor.b_end = accessor.matrix->inverse(accessor.a_block).end();
1015  }
1016  }
1017  return *this;
1018 }
1019 
1020 
1021 template <typename MatrixType, typename inverse_type>
1023  const_iterator::Accessor &
1025  const
1026 {
1027  return accessor;
1028 }
1029 
1030 
1031 template <typename MatrixType, typename inverse_type>
1033  const_iterator::Accessor *
1034  PreconditionBlockJacobi<MatrixType,
1035  inverse_type>::const_iterator::operator->() const
1036 {
1037  return &accessor;
1038 }
1039 
1040 
1041 template <typename MatrixType, typename inverse_type>
1042 inline bool
1044  const const_iterator &other) const
1045 {
1046  if (accessor.a_block == accessor.matrix->size() &&
1047  accessor.a_block == other.accessor.a_block)
1048  return true;
1049 
1050  if (accessor.a_block != other.accessor.a_block)
1051  return false;
1052 
1053  return (accessor.row() == other.accessor.row() &&
1054  accessor.column() == other.accessor.column());
1055 }
1056 
1057 
1058 template <typename MatrixType, typename inverse_type>
1059 inline bool
1061  const const_iterator &other) const
1062 {
1063  return !(*this == other);
1064 }
1065 
1066 
1067 template <typename MatrixType, typename inverse_type>
1068 inline bool
1070  const const_iterator &other) const
1071 {
1072  return (accessor.row() < other.accessor.row() ||
1073  (accessor.row() == other.accessor.row() &&
1074  accessor.column() < other.accessor.column()));
1075 }
1076 
1077 
1078 template <typename MatrixType, typename inverse_type>
1079 inline
1082 {
1083  return const_iterator(this, 0);
1084 }
1085 
1086 
1087 template <typename MatrixType, typename inverse_type>
1088 inline
1091 {
1092  return const_iterator(this, this->size() * this->block_size());
1093 }
1094 
1095 
1096 template <typename MatrixType, typename inverse_type>
1097 inline
1100  const size_type r) const
1101 {
1102  AssertIndexRange(r, this->A->m());
1103  return const_iterator(this, r);
1104 }
1105 
1106 
1107 
1108 template <typename MatrixType, typename inverse_type>
1109 inline
1112  const size_type r) const
1113 {
1114  AssertIndexRange(r, this->A->m());
1115  return const_iterator(this, r + 1);
1116 }
1117 
1118 #endif // DOXYGEN
1119 
1121 
1122 #endif
FullMatrix< typename MatrixType::value_type > & inverse(size_type i)
const PreconditionBlockJacobi< MatrixType, inverse_type > * matrix
FullMatrix< inverse_type >::const_iterator b_end
FullMatrix< inverse_type >::const_iterator b_iterator
Accessor(const PreconditionBlockJacobi< MatrixType, inverse_type > *matrix, const size_type row)
const Accessor * operator->() const
bool operator!=(const const_iterator &) const
const_iterator(const PreconditionBlockJacobi< MatrixType, inverse_type > *matrix, const size_type row)
const Accessor & operator*() const
bool operator==(const const_iterator &) const
bool operator<(const const_iterator &) const
void vmult(Vector< number2 > &, const Vector< number2 > &) const
void step(Vector< number2 > &dst, const Vector< number2 > &rhs) const
void vmult_add(Vector< number2 > &, const Vector< number2 > &) const
const_iterator begin(const size_type r) const
const_iterator end() const
void Tvmult(Vector< number2 > &, const Vector< number2 > &) const
void do_vmult(Vector< number2 > &, const Vector< number2 > &, bool adding) const
void Tvmult_add(Vector< number2 > &, const Vector< number2 > &) const
typename MatrixType::value_type number
const_iterator begin() const
const_iterator end(const size_type r) const
void Tstep(Vector< number2 > &dst, const Vector< number2 > &rhs) const
types::global_dof_index size_type
types::global_dof_index size_type
PreconditionBlockSOR(bool store)
void vmult_add(Vector< number2 > &, const Vector< number2 > &) const
void vmult(Vector< number2 > &, const Vector< number2 > &) const
typename MatrixType::value_type number
void step(Vector< number2 > &dst, const Vector< number2 > &rhs) const
void Tstep(Vector< number2 > &dst, const Vector< number2 > &rhs) const
void forward(Vector< number2 > &, const Vector< number2 > &, const bool transpose_diagonal, const bool adding) const
void backward(Vector< number2 > &, const Vector< number2 > &, const bool transpose_diagonal, const bool adding) const
void Tvmult(Vector< number2 > &, const Vector< number2 > &) const
void Tvmult_add(Vector< number2 > &, const Vector< number2 > &) const
void vmult(Vector< number2 > &, const Vector< number2 > &) const
void step(Vector< number2 > &dst, const Vector< number2 > &rhs) const
typename MatrixType::value_type number
void Tstep(Vector< number2 > &dst, const Vector< number2 > &rhs) const
types::global_dof_index size_type
void Tvmult(Vector< number2 > &, const Vector< number2 > &) const
AdditionalData(const size_type block_size, const double relaxation=1., const bool invert_diagonal=true, const bool same_diagonal=false)
PreconditionBlockBase< inverse_type >::Inversion inversion
PreconditionBlock(bool store_diagonals=false)
~PreconditionBlock() override=default
value_type el(size_type i, size_type j) const
SmartPointer< const MatrixType, PreconditionBlock< MatrixType, inverse_type > > A
void forward_step(Vector< number2 > &dst, const Vector< number2 > &prev, const Vector< number2 > &src, const bool transpose_diagonal) const
bool empty() const
void initialize(const MatrixType &A, const AdditionalData parameters)
void invert_diagblocks()
types::global_dof_index size_type
void invert_permuted_diagblocks()
void backward_step(Vector< number2 > &dst, const Vector< number2 > &prev, const Vector< number2 > &src, const bool transpose_diagonal) const
std::size_t memory_consumption() const
std::vector< size_type > inverse_permutation
inverse_type value_type
void set_permutation(const std::vector< size_type > &permutation, const std::vector< size_type > &inverse_permutation)
std::vector< size_type > permutation
size_type block_size() const
typename MatrixType::value_type number
void initialize(const MatrixType &A, const std::vector< size_type > &permutation, const std::vector< size_type > &inverse_permutation, const AdditionalData parameters)
Definition: vector.h:109
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:474
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:475
static ::ExceptionBase & ExcWrongBlockSize(int arg1, int arg2)
static ::ExceptionBase & ExcInverseMatricesAlreadyExist()
#define DeclException0(Exception0)
Definition: exceptions.h:465
static ::ExceptionBase & ExcIteratorPastEnd()
#define Assert(cond, exc)
Definition: exceptions.h:1614
#define DeclException2(Exception2, type1, type2, outsequence)
Definition: exceptions.h:533
#define AssertIndexRange(index, range)
Definition: exceptions.h:1855
constexpr int block_size
Definition: cuda_size.h:29
@ matrix
Contents is actually a matrix.
types::global_dof_index size_type
Definition: cuda_kernels.h:45
unsigned int global_dof_index
Definition: types.h:82