Reference documentation for deal.II version GIT relicensing-426-g7976cfd195 2024-04-18 21:10:01+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\}}\)
Loading...
Searching...
No Matches
precondition_block.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 1999 - 2023 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#ifndef dealii_precondition_block_h
16#define dealii_precondition_block_h
17
18
19#include <deal.II/base/config.h>
20
24
26
27#include <vector>
28
30
80template <typename MatrixType,
81 typename inverse_type = typename MatrixType::value_type>
82class PreconditionBlock : public virtual Subscriptor,
83 protected PreconditionBlockBase<inverse_type>
84{
85private:
89 using number = typename MatrixType::value_type;
90
94 using value_type = inverse_type;
95
96public:
101
106 {
107 public:
113 const double relaxation = 1.,
114 const bool invert_diagonal = true,
115 const bool same_diagonal = false);
116
121
126
131
140
146 double threshold;
147 };
148
149
154
158 ~PreconditionBlock() override = default;
159
167 void
168 initialize(const MatrixType &A, const AdditionalData parameters);
169
170protected:
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
221public:
227 void
229
233 bool
234 empty() const;
235
241 el(size_type i, size_type j) const;
242
258 void
260
272 template <typename number2>
273 void
275 const Vector<number2> &prev,
276 const Vector<number2> &src,
277 const bool transpose_diagonal) const;
278
290 template <typename number2>
291 void
293 const Vector<number2> &prev,
294 const Vector<number2> &src,
295 const bool transpose_diagonal) const;
296
297
302 block_size() const;
303
308 std::size_t
310
321 int,
322 int,
323 << "The blocksize " << arg1 << " and the size of the matrix "
324 << arg2 << " do not match.");
325
330
333protected:
339
351
355 std::vector<size_type> permutation;
356
360 std::vector<size_type> inverse_permutation;
361};
362
363
364
376template <typename MatrixType,
377 typename inverse_type = typename MatrixType::value_type>
379 : public virtual Subscriptor,
380 private PreconditionBlock<MatrixType, inverse_type>
381{
382private:
386 using number = typename MatrixType::value_type;
387
388public:
393
398 {
399 private:
404 {
405 public:
411 const size_type row);
412
417 row() const;
418
423 column() const;
424
428 inverse_type
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
474
478 const Accessor &
479 operator*() const;
480
484 const Accessor *
485 operator->() const;
486
490 bool
495 bool
497
502 bool
503 operator<(const const_iterator &) const;
504
505 private:
510 };
511
516 using PreconditionBlock<MatrixType, inverse_type>::initialize;
517 using PreconditionBlock<MatrixType, inverse_type>::clear;
518 using PreconditionBlock<MatrixType, inverse_type>::empty;
519 using PreconditionBlock<MatrixType, inverse_type>::el;
520 using PreconditionBlock<MatrixType, inverse_type>::invert_diagblocks;
521 using PreconditionBlock<MatrixType, inverse_type>::block_size;
522 using PreconditionBlockBase<inverse_type>::size;
523 using PreconditionBlockBase<inverse_type>::inverse;
525 using PreconditionBlockBase<inverse_type>::inverse_svd;
526 using PreconditionBlockBase<inverse_type>::log_statistics;
527 using PreconditionBlock<MatrixType, inverse_type>::set_permutation;
528
536 template <typename number2>
537 void
539
543 template <typename number2>
544 void
553 template <typename number2>
554 void
556
560 template <typename number2>
561 void
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
603private:
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
653template <typename MatrixType,
654 typename inverse_type = typename MatrixType::value_type>
656 : public virtual Subscriptor,
657 protected PreconditionBlock<MatrixType, inverse_type>
658{
659public:
664
669
673 using number = typename MatrixType::value_type;
674
679 using PreconditionBlock<MatrixType, inverse_type>::initialize;
680 using PreconditionBlock<MatrixType, inverse_type>::clear;
681 using PreconditionBlockBase<inverse_type>::size;
682 using PreconditionBlockBase<inverse_type>::inverse;
684 using PreconditionBlockBase<inverse_type>::inverse_svd;
685 using PreconditionBlock<MatrixType, inverse_type>::invert_diagblocks;
686 using PreconditionBlock<MatrixType, inverse_type>::set_permutation;
687 using PreconditionBlockBase<inverse_type>::log_statistics;
688
699 template <typename number2>
700 void
702
713 template <typename number2>
714 void
716
725 template <typename number2>
726 void
728
739 template <typename number2>
740 void
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
757protected:
762
772 template <typename number2>
773 void
775 const Vector<number2> &,
776 const bool transpose_diagonal,
777 const bool adding) const;
778
788 template <typename number2>
789 void
791 const Vector<number2> &,
792 const bool transpose_diagonal,
793 const bool adding) const;
794};
795
796
816template <typename MatrixType,
817 typename inverse_type = typename MatrixType::value_type>
819 : public virtual Subscriptor,
820 private PreconditionBlockSOR<MatrixType, inverse_type>
821{
822public:
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.
848 using PreconditionBlockSOR<MatrixType, inverse_type>::initialize;
849 using PreconditionBlockSOR<MatrixType, inverse_type>::clear;
850 using PreconditionBlockBase<inverse_type>::size;
851 using PreconditionBlockBase<inverse_type>::inverse;
853 using PreconditionBlockBase<inverse_type>::inverse_svd;
854 using PreconditionBlockBase<inverse_type>::log_statistics;
855 using PreconditionBlockSOR<MatrixType, inverse_type>::set_permutation;
856 using PreconditionBlockSOR<MatrixType, inverse_type>::empty;
857 using PreconditionBlockSOR<MatrixType, inverse_type>::el;
858 using PreconditionBlockSOR<MatrixType, inverse_type>::invert_diagblocks;
859
867 template <typename number2>
868 void
870
874 template <typename number2>
875 void
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
898template <typename MatrixType, typename inverse_type>
899inline bool
901{
902 if (A == nullptr)
903 return true;
904 return A->empty();
905}
906
907
908template <typename MatrixType, typename inverse_type>
909inline inverse_type
910PreconditionBlock<MatrixType, inverse_type>::el(size_type i, size_type j) const
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
930template <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
955template <typename MatrixType, typename inverse_type>
957PreconditionBlockJacobi<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
966template <typename MatrixType, typename inverse_type>
968PreconditionBlockJacobi<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
977template <typename MatrixType, typename inverse_type>
978inline inverse_type
979PreconditionBlockJacobi<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
988template <typename MatrixType, typename inverse_type>
992 const size_type row)
993 : accessor(matrix, row)
994{}
995
996
997template <typename MatrixType, typename inverse_type>
998inline typename PreconditionBlockJacobi<MatrixType,
999 inverse_type>::const_iterator &
1001{
1003
1006 {
1007 ++accessor.a_block;
1008
1010 {
1014 }
1015 }
1016 return *this;
1017}
1018
1019
1020template <typename MatrixType, typename inverse_type>
1022 const_iterator::Accessor &
1024 const
1025{
1026 return accessor;
1027}
1028
1029
1030template <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
1040template <typename MatrixType, typename inverse_type>
1041inline 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
1057template <typename MatrixType, typename inverse_type>
1058inline bool
1060 const const_iterator &other) const
1061{
1062 return !(*this == other);
1063}
1064
1065
1066template <typename MatrixType, typename inverse_type>
1067inline 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
1077template <typename MatrixType, typename inverse_type>
1078inline
1081{
1082 return const_iterator(this, 0);
1083}
1084
1085
1086template <typename MatrixType, typename inverse_type>
1087inline
1090{
1091 return const_iterator(this, this->size() * this->block_size());
1092}
1093
1094
1095template <typename MatrixType, typename inverse_type>
1096inline
1099 const size_type r) const
1100{
1101 AssertIndexRange(r, this->A->m());
1102 return const_iterator(this, r);
1103}
1104
1105
1106
1107template <typename MatrixType, typename inverse_type>
1108inline
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
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:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
#define DeclException0(Exception0)
Definition exceptions.h:471
static ::ExceptionBase & ExcInverseMatricesAlreadyExist()
#define Assert(cond, exc)
static ::ExceptionBase & ExcIteratorPastEnd()
#define DeclException2(Exception2, type1, type2, outsequence)
Definition exceptions.h:539
static ::ExceptionBase & ExcWrongBlockSize(int arg1, int arg2)
#define AssertIndexRange(index, range)
constexpr int block_size
Definition cuda_size.h:28
@ matrix
Contents is actually a matrix.
types::global_dof_index size_type
unsigned int global_dof_index
Definition types.h:81