Reference documentation for deal.II version 9.5.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\}}\)
Loading...
Searching...
No Matches
precondition_block.h
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 1999 - 2022 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
81template <typename MatrixType,
82 typename inverse_type = typename MatrixType::value_type>
83class PreconditionBlock : public virtual Subscriptor,
84 protected PreconditionBlockBase<inverse_type>
85{
86private:
90 using number = typename MatrixType::value_type;
91
95 using value_type = inverse_type;
96
97public:
102
107 {
108 public:
114 const double relaxation = 1.,
115 const bool invert_diagonal = true,
116 const bool same_diagonal = false);
117
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
171protected:
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
222public:
228 void
230
234 bool
235 empty() const;
236
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
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
334protected:
340
352
356 std::vector<size_type> permutation;
357
361 std::vector<size_type> inverse_permutation;
362};
363
364
365
377template <typename MatrixType,
378 typename inverse_type = typename MatrixType::value_type>
380 : public virtual Subscriptor,
381 private PreconditionBlock<MatrixType, inverse_type>
382{
383private:
387 using number = typename MatrixType::value_type;
388
389public:
394
399 {
400 private:
405 {
406 public:
412 const size_type row);
413
418 row() const;
419
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
496 bool
498
503 bool
504 operator<(const const_iterator &) const;
505
506 private:
511 };
512
517 using PreconditionBlock<MatrixType, inverse_type>::initialize;
518 using PreconditionBlock<MatrixType, inverse_type>::clear;
519 using PreconditionBlock<MatrixType, inverse_type>::empty;
520 using PreconditionBlock<MatrixType, inverse_type>::el;
521 using PreconditionBlock<MatrixType, inverse_type>::invert_diagblocks;
522 using PreconditionBlock<MatrixType, inverse_type>::block_size;
523 using PreconditionBlockBase<inverse_type>::size;
524 using PreconditionBlockBase<inverse_type>::inverse;
526 using PreconditionBlockBase<inverse_type>::inverse_svd;
527 using PreconditionBlockBase<inverse_type>::log_statistics;
528 using PreconditionBlock<MatrixType, inverse_type>::set_permutation;
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
604private:
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
654template <typename MatrixType,
655 typename inverse_type = typename MatrixType::value_type>
657 : public virtual Subscriptor,
658 protected PreconditionBlock<MatrixType, inverse_type>
659{
660public:
665
670
674 using number = typename MatrixType::value_type;
675
680 using PreconditionBlock<MatrixType, inverse_type>::initialize;
681 using PreconditionBlock<MatrixType, inverse_type>::clear;
682 using PreconditionBlockBase<inverse_type>::size;
683 using PreconditionBlockBase<inverse_type>::inverse;
685 using PreconditionBlockBase<inverse_type>::inverse_svd;
686 using PreconditionBlock<MatrixType, inverse_type>::invert_diagblocks;
687 using PreconditionBlock<MatrixType, inverse_type>::set_permutation;
688 using PreconditionBlockBase<inverse_type>::log_statistics;
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
758protected:
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
817template <typename MatrixType,
818 typename inverse_type = typename MatrixType::value_type>
820 : public virtual Subscriptor,
821 private PreconditionBlockSOR<MatrixType, inverse_type>
822{
823public:
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.
849 using PreconditionBlockSOR<MatrixType, inverse_type>::initialize;
850 using PreconditionBlockSOR<MatrixType, inverse_type>::clear;
851 using PreconditionBlockBase<inverse_type>::size;
852 using PreconditionBlockBase<inverse_type>::inverse;
854 using PreconditionBlockBase<inverse_type>::inverse_svd;
855 using PreconditionBlockBase<inverse_type>::log_statistics;
856 using PreconditionBlockSOR<MatrixType, inverse_type>::set_permutation;
857 using PreconditionBlockSOR<MatrixType, inverse_type>::empty;
858 using PreconditionBlockSOR<MatrixType, inverse_type>::el;
859 using PreconditionBlockSOR<MatrixType, inverse_type>::invert_diagblocks;
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
899template <typename MatrixType, typename inverse_type>
900inline bool
902{
903 if (A == nullptr)
904 return true;
905 return A->empty();
906}
907
908
909template <typename MatrixType, typename inverse_type>
910inline inverse_type
911PreconditionBlock<MatrixType, inverse_type>::el(size_type i, size_type j) const
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
931template <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
956template <typename MatrixType, typename inverse_type>
958PreconditionBlockJacobi<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
967template <typename MatrixType, typename inverse_type>
969PreconditionBlockJacobi<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
978template <typename MatrixType, typename inverse_type>
979inline inverse_type
980PreconditionBlockJacobi<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
989template <typename MatrixType, typename inverse_type>
993 const size_type row)
994 : accessor(matrix, row)
995{}
996
997
998template <typename MatrixType, typename inverse_type>
999inline typename PreconditionBlockJacobi<MatrixType,
1000 inverse_type>::const_iterator &
1002{
1004
1007 {
1008 ++accessor.a_block;
1009
1011 {
1015 }
1016 }
1017 return *this;
1018}
1019
1020
1021template <typename MatrixType, typename inverse_type>
1023 const_iterator::Accessor &
1025 const
1026{
1027 return accessor;
1028}
1029
1030
1031template <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
1041template <typename MatrixType, typename inverse_type>
1042inline 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
1058template <typename MatrixType, typename inverse_type>
1059inline bool
1061 const const_iterator &other) const
1062{
1063 return !(*this == other);
1064}
1065
1066
1067template <typename MatrixType, typename inverse_type>
1068inline 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
1078template <typename MatrixType, typename inverse_type>
1079inline
1082{
1083 return const_iterator(this, 0);
1084}
1085
1086
1087template <typename MatrixType, typename inverse_type>
1088inline
1091{
1092 return const_iterator(this, this->size() * this->block_size());
1093}
1094
1095
1096template <typename MatrixType, typename inverse_type>
1097inline
1100 const size_type r) const
1101{
1102 AssertIndexRange(r, this->A->m());
1103 return const_iterator(this, r);
1104}
1105
1106
1107
1108template <typename MatrixType, typename inverse_type>
1109inline
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
typename Table< 2, number >::const_iterator const_iterator
iterator begin(const size_type r)
iterator end(const size_type r)
Householder< number > & inverse_householder(size_type i)
FullMatrix< number > & inverse(size_type i)
LAPACKFullMatrix< number > & inverse_svd(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)
bool operator!=(const const_iterator &) const
const_iterator(const PreconditionBlockJacobi< MatrixType, inverse_type > *matrix, const size_type row)
bool operator==(const const_iterator &) const
bool operator<(const const_iterator &) const
const Accessor & operator*() const
const Accessor * operator->() 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
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()
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
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)
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
#define DeclException0(Exception0)
Definition exceptions.h:465
static ::ExceptionBase & ExcInverseMatricesAlreadyExist()
#define Assert(cond, exc)
static ::ExceptionBase & ExcIteratorPastEnd()
#define DeclException2(Exception2, type1, type2, outsequence)
Definition exceptions.h:533
static ::ExceptionBase & ExcWrongBlockSize(int arg1, int arg2)
#define AssertIndexRange(index, range)
constexpr int block_size
Definition cuda_size.h:29
@ matrix
Contents is actually a matrix.
types::global_dof_index size_type
unsigned int global_dof_index
Definition types.h:82