Loading [MathJax]/extensions/TeX/newcommand.js
 deal.II version GIT relicensing-2791-g2e2a880805 2025-03-11 02:00:00+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\}}
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages Concepts
mg_smoother.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_mg_smoother_h
16#define dealii_mg_smoother_h
17
18
19#include <deal.II/base/config.h>
20
24
27
29
30#include <vector>
31
33
34/*
35 * MGSmootherBase is defined in mg_base.h
36 */
37
49template <typename VectorType>
50class MGSmoother : public MGSmootherBase<VectorType>
51{
52public:
56 MGSmoother(const unsigned int steps = 1,
57 const bool variable = false,
58 const bool symmetric = false,
59 const bool transpose = false);
60
64 void
65 set_steps(const unsigned int);
66
70 void
71 set_variable(const bool);
72
76 void
77 set_symmetric(const bool);
78
83 void
84 set_transpose(const bool);
85
90 void
91 set_debug(const unsigned int level);
92
93protected:
101
106 unsigned int steps;
107
113
119
125
129 unsigned int debug;
130};
131
132
139template <typename VectorType>
140class MGSmootherIdentity : public MGSmootherBase<VectorType>
141{
142public:
148 virtual void
149 smooth(const unsigned int level, VectorType &u, const VectorType &rhs) const;
150
151 virtual void
153};
154
155
156namespace mg
157{
184 template <typename RelaxationType, typename VectorType>
185 class SmootherRelaxation : public MGLevelObject<RelaxationType>,
186 public MGSmoother<VectorType>
187 {
188 public:
192 SmootherRelaxation(const unsigned int steps = 1,
193 const bool variable = false,
194 const bool symmetric = false,
195 const bool transpose = false);
196
205 template <typename MatrixType2>
206 void
208 const typename RelaxationType::AdditionalData &additional_data =
209 typename RelaxationType::AdditionalData());
210
218 template <typename MatrixType2, typename DataType>
219 void
221 const MGLevelObject<DataType> &additional_data);
222
226 void
227 clear() override;
228
232 virtual void
233 smooth(const unsigned int level,
234 VectorType &u,
235 const VectorType &rhs) const override;
236
253 virtual void
254 apply(const unsigned int level,
255 VectorType &u,
256 const VectorType &rhs) const override;
257
261 std::size_t
263 };
264} // namespace mg
265
295template <typename MatrixType, typename RelaxationType, typename VectorType>
296class MGSmootherRelaxation : public MGSmoother<VectorType>
297{
298public:
302 MGSmootherRelaxation(const unsigned int steps = 1,
303 const bool variable = false,
304 const bool symmetric = false,
305 const bool transpose = false);
306
315 template <typename MatrixType2>
316 void
318 const typename RelaxationType::AdditionalData &additional_data =
319 typename RelaxationType::AdditionalData());
320
329 template <typename MatrixType2, typename DataType>
330 void
332 const MGLevelObject<DataType> &additional_data);
333
343 template <typename MatrixType2, typename DataType>
344 void
346 const DataType &additional_data,
347 const unsigned int block_row,
348 const unsigned int block_col);
349
359 template <typename MatrixType2, typename DataType>
360 void
362 const MGLevelObject<DataType> &additional_data,
363 const unsigned int block_row,
364 const unsigned int block_col);
365
369 void
371
375 virtual void
376 smooth(const unsigned int level, VectorType &u, const VectorType &rhs) const;
377
393 virtual void
394 apply(const unsigned int level, VectorType &u, const VectorType &rhs) const;
395
400
404 std::size_t
406
407
408private:
413};
414
415
416
445template <typename MatrixType, typename PreconditionerType, typename VectorType>
446class MGSmootherPrecondition : public MGSmoother<VectorType>
447{
448public:
452 MGSmootherPrecondition(const unsigned int steps = 1,
453 const bool variable = false,
454 const bool symmetric = false,
455 const bool transpose = false);
456
466 template <typename MatrixType2>
467 void
469 const typename PreconditionerType::AdditionalData &
470 additional_data = typename PreconditionerType::AdditionalData());
471
479 template <typename MatrixType2>
480 void
482
492 template <typename MatrixType2, typename DataType>
493 void
495 const MGLevelObject<DataType> &additional_data);
496
507 template <typename MatrixType2, typename DataType>
508 void
510 const DataType &additional_data,
511 const unsigned int block_row,
512 const unsigned int block_col);
513
524 template <typename MatrixType2, typename DataType>
525 void
527 const MGLevelObject<DataType> &additional_data,
528 const unsigned int block_row,
529 const unsigned int block_col);
530
534 void
535 clear() override;
536
540 virtual void
541 smooth(const unsigned int level,
542 VectorType &u,
543 const VectorType &rhs) const override;
544
560 virtual void
561 apply(const unsigned int level,
562 VectorType &u,
563 const VectorType &rhs) const override;
564
569
573 std::size_t
575
576
577private:
582};
583
586/* ------------------------------- Inline functions --------------------------
587 */
588
589#ifndef DOXYGEN
590
591template <typename VectorType>
592inline void
594 VectorType &,
595 const VectorType &) const
596{}
597
598template <typename VectorType>
599inline void
601{}
602
603//---------------------------------------------------------------------------
604
605template <typename VectorType>
606inline MGSmoother<VectorType>::MGSmoother(const unsigned int steps,
607 const bool variable,
608 const bool symmetric,
609 const bool transpose)
610 : steps(steps)
611 , variable(variable)
614 , debug(0)
615{}
616
617
618template <typename VectorType>
619inline void
620MGSmoother<VectorType>::set_steps(const unsigned int s)
621{
622 steps = s;
623}
624
625
626template <typename VectorType>
627inline void
628MGSmoother<VectorType>::set_debug(const unsigned int s)
629{
630 debug = s;
631}
632
633
634template <typename VectorType>
635inline void
637{
638 variable = flag;
639}
640
641
642template <typename VectorType>
643inline void
645{
646 symmetric = flag;
647}
648
649
650template <typename VectorType>
651inline void
653{
654 transpose = flag;
655}
656
657//----------------------------------------------------------------------//
658
659namespace mg
660{
661 template <typename RelaxationType, typename VectorType>
663 const unsigned int steps,
664 const bool variable,
665 const bool symmetric,
666 const bool transpose)
667 : MGSmoother<VectorType>(steps, variable, symmetric, transpose)
668 {}
669
670
671 template <typename RelaxationType, typename VectorType>
672 inline void
673 SmootherRelaxation<RelaxationType, VectorType>::clear()
674 {
676 }
677
678
679 template <typename RelaxationType, typename VectorType>
680 template <typename MatrixType2>
681 inline void
682 SmootherRelaxation<RelaxationType, VectorType>::initialize(
684 const typename RelaxationType::AdditionalData &data)
685 {
686 const unsigned int min = m.min_level();
687 const unsigned int max = m.max_level();
688
689 this->resize(min, max);
690
691 for (unsigned int i = min; i <= max; ++i)
692 (*this)[i].initialize(Utilities::get_underlying_value(m[i]), data);
693 }
694
695
696 template <typename RelaxationType, typename VectorType>
697 template <typename MatrixType2, typename DataType>
698 inline void
699 SmootherRelaxation<RelaxationType, VectorType>::initialize(
702 {
703 const unsigned int min = std::max(m.min_level(), data.min_level());
704 const unsigned int max = std::min(m.max_level(), data.max_level());
705
706 this->resize(min, max);
707
708 for (unsigned int i = min; i <= max; ++i)
709 (*this)[i].initialize(Utilities::get_underlying_value(m[i]), data[i]);
710 }
711
712
713 template <typename RelaxationType, typename VectorType>
714 inline void
715 SmootherRelaxation<RelaxationType, VectorType>::smooth(
716 const unsigned int level,
717 VectorType &u,
718 const VectorType &rhs) const
719 {
720 unsigned int maxlevel = this->max_level();
721 unsigned int steps2 = this->steps;
722
723 if (this->variable)
724 steps2 *= (1 << (maxlevel - level));
725
726 bool T = this->transpose;
727 if (this->symmetric && (steps2 % 2 == 0))
728 T = false;
729 if (this->debug > 0)
730 deallog << 'S' << level << ' ';
731
732 for (unsigned int i = 0; i < steps2; ++i)
733 {
734 if (T)
735 (*this)[level].Tstep(u, rhs);
736 else
737 (*this)[level].step(u, rhs);
738 if (this->symmetric)
739 T = !T;
740 }
741 }
742
743
744 template <typename RelaxationType, typename VectorType>
745 inline void
746 SmootherRelaxation<RelaxationType, VectorType>::apply(
747 const unsigned int level,
748 VectorType &u,
749 const VectorType &rhs) const
750 {
751 unsigned int maxlevel = this->max_level();
752 unsigned int steps2 = this->steps;
753
754 if (this->variable)
755 steps2 *= (1 << (maxlevel - level));
756
757 bool T = this->transpose;
758 if (this->symmetric && (steps2 % 2 == 0))
759 T = false;
760 if (this->debug > 0)
761 deallog << 'S' << level << ' ';
762
763 if (T)
764 (*this)[level].Tvmult(u, rhs);
765 else
766 (*this)[level].vmult(u, rhs);
767 if (this->symmetric)
768 T = !T;
769 for (unsigned int i = 1; i < steps2; ++i)
770 {
771 if (T)
772 (*this)[level].Tstep(u, rhs);
773 else
774 (*this)[level].step(u, rhs);
775 if (this->symmetric)
776 T = !T;
777 }
778 }
779
780
781 template <typename RelaxationType, typename VectorType>
782 inline std::size_t
783 SmootherRelaxation<RelaxationType, VectorType>::memory_consumption() const
784 {
785 return sizeof(*this) - sizeof(MGLevelObject<RelaxationType>) +
787 this->vector_memory.memory_consumption();
788 }
789} // namespace mg
790
791
792//----------------------------------------------------------------------//
793
794template <typename MatrixType, typename RelaxationType, typename VectorType>
796 MGSmootherRelaxation(const unsigned int steps,
797 const bool variable,
798 const bool symmetric,
799 const bool transpose)
800 : MGSmoother<VectorType>(steps, variable, symmetric, transpose)
801{}
802
803
804
805template <typename MatrixType, typename RelaxationType, typename VectorType>
806inline void
808{
809 smoothers.clear_elements();
810
811 unsigned int i = matrices.min_level(), max_level = matrices.max_level();
812 for (; i <= max_level; ++i)
813 matrices[i] = LinearOperator<VectorType>();
814}
815
816
817template <typename MatrixType, typename RelaxationType, typename VectorType>
818template <typename MatrixType2>
819inline void
822 const typename RelaxationType::AdditionalData &data)
823{
824 const unsigned int min = m.min_level();
825 const unsigned int max = m.max_level();
826
827 matrices.resize(min, max);
828 smoothers.resize(min, max);
829
830 for (unsigned int i = min; i <= max; ++i)
831 {
832 // Workaround: Unfortunately, not every "m[i]" object has a rich
833 // enough interface to populate reinit_(domain|range)_vector. Thus,
834 // apply an empty LinearOperator exemplar.
835 matrices[i] =
836 linear_operator<VectorType>(LinearOperator<VectorType>(),
838 smoothers[i].initialize(Utilities::get_underlying_value(m[i]), data);
839 }
840}
841
842template <typename MatrixType, typename RelaxationType, typename VectorType>
843template <typename MatrixType2, typename DataType>
844inline void
848{
849 const unsigned int min = m.min_level();
850 const unsigned int max = m.max_level();
851
852 Assert(data.min_level() == min, ExcDimensionMismatch(data.min_level(), min));
853 Assert(data.max_level() == max, ExcDimensionMismatch(data.max_level(), max));
854
855 matrices.resize(min, max);
856 smoothers.resize(min, max);
857
858 for (unsigned int i = min; i <= max; ++i)
859 {
860 // Workaround: Unfortunately, not every "m[i]" object has a rich
861 // enough interface to populate reinit_(domain|range)_vector. Thus,
862 // apply an empty LinearOperator exemplar.
863 matrices[i] =
864 linear_operator<VectorType>(LinearOperator<VectorType>(),
866 smoothers[i].initialize(Utilities::get_underlying_value(m[i]), data[i]);
867 }
868}
869
870template <typename MatrixType, typename RelaxationType, typename VectorType>
871template <typename MatrixType2, typename DataType>
872inline void
875 const DataType &data,
876 const unsigned int row,
877 const unsigned int col)
878{
879 const unsigned int min = m.min_level();
880 const unsigned int max = m.max_level();
881
882 matrices.resize(min, max);
883 smoothers.resize(min, max);
884
885 for (unsigned int i = min; i <= max; ++i)
886 {
887 // Workaround: Unfortunately, not every "m[i]" object has a rich
888 // enough interface to populate reinit_(domain|range)_vector. Thus,
889 // apply an empty LinearOperator exemplar.
890 matrices[i] = linear_operator<VectorType>(LinearOperator<VectorType>(),
891 m[i].block(row, col));
892 smoothers[i].initialize(m[i].block(row, col), data);
893 }
894}
895
896template <typename MatrixType, typename RelaxationType, typename VectorType>
897template <typename MatrixType2, typename DataType>
898inline void
902 const unsigned int row,
903 const unsigned int col)
904{
905 const unsigned int min = m.min_level();
906 const unsigned int max = m.max_level();
907
908 Assert(data.min_level() == min, ExcDimensionMismatch(data.min_level(), min));
909 Assert(data.max_level() == max, ExcDimensionMismatch(data.max_level(), max));
910
911 matrices.resize(min, max);
912 smoothers.resize(min, max);
913
914 for (unsigned int i = min; i <= max; ++i)
915 {
916 // Workaround: Unfortunately, not every "m[i]" object has a rich
917 // enough interface to populate reinit_(domain|range)_vector. Thus,
918 // apply an empty LinearOperator exemplar.
919 matrices[i] = linear_operator<VectorType>(LinearOperator<VectorType>(),
920 m[i].block(row, col));
921 smoothers[i].initialize(m[i].block(row, col), data[i]);
922 }
923}
924
925
926template <typename MatrixType, typename RelaxationType, typename VectorType>
927inline void
929 const unsigned int level,
930 VectorType &u,
931 const VectorType &rhs) const
932{
933 unsigned int maxlevel = smoothers.max_level();
934 unsigned int steps2 = this->steps;
935
936 if (this->variable)
937 steps2 *= (1 << (maxlevel - level));
938
939 bool T = this->transpose;
940 if (this->symmetric && (steps2 % 2 == 0))
941 T = false;
942 if (this->debug > 0)
943 deallog << 'S' << level << ' ';
944
945 for (unsigned int i = 0; i < steps2; ++i)
946 {
947 if (T)
948 smoothers[level].Tstep(u, rhs);
949 else
950 smoothers[level].step(u, rhs);
951 if (this->symmetric)
952 T = !T;
953 }
954}
955
956
957template <typename MatrixType, typename RelaxationType, typename VectorType>
958inline void
960 const unsigned int level,
961 VectorType &u,
962 const VectorType &rhs) const
963{
964 unsigned int maxlevel = smoothers.max_level();
965 unsigned int steps2 = this->steps;
966
967 if (this->variable)
968 steps2 *= (1 << (maxlevel - level));
969
970 bool T = this->transpose;
971 if (this->symmetric && (steps2 % 2 == 0))
972 T = false;
973 if (this->debug > 0)
974 deallog << 'S' << level << ' ';
975
976 if (T)
977 smoothers[level].Tvmult(u, rhs);
978 else
979 smoothers[level].vmult(u, rhs);
980 if (this->symmetric)
981 T = !T;
982 for (unsigned int i = 1; i < steps2; ++i)
983 {
984 if (T)
985 smoothers[level].Tstep(u, rhs);
986 else
987 smoothers[level].step(u, rhs);
988 if (this->symmetric)
989 T = !T;
990 }
991}
992
993
994
995template <typename MatrixType, typename RelaxationType, typename VectorType>
996inline std::size_t
999{
1000 return sizeof(*this) + matrices.memory_consumption() +
1001 smoothers.memory_consumption() +
1002 this->vector_memory.memory_consumption();
1003}
1004
1005
1006//----------------------------------------------------------------------//
1007
1008template <typename MatrixType, typename PreconditionerType, typename VectorType>
1010 MGSmootherPrecondition(const unsigned int steps,
1011 const bool variable,
1012 const bool symmetric,
1013 const bool transpose)
1014 : MGSmoother<VectorType>(steps, variable, symmetric, transpose)
1015{}
1016
1017
1018
1019template <typename MatrixType, typename PreconditionerType, typename VectorType>
1020inline void
1022{
1023 smoothers.clear_elements();
1024
1025 unsigned int i = matrices.min_level(), max_level = matrices.max_level();
1026 for (; i <= max_level; ++i)
1027 matrices[i] = LinearOperator<VectorType>();
1028}
1029
1030
1031
1032template <typename MatrixType, typename PreconditionerType, typename VectorType>
1033template <typename MatrixType2>
1034inline void
1037 const typename PreconditionerType::AdditionalData &data)
1038{
1039 const unsigned int min = m.min_level();
1040 const unsigned int max = m.max_level();
1041
1042 matrices.resize(min, max);
1043 smoothers.resize(min, max);
1044
1045 for (unsigned int i = min; i <= max; ++i)
1046 {
1047 // Workaround: Unfortunately, not every "m[i]" object has a rich
1048 // enough interface to populate reinit_(domain|range)_vector. Thus,
1049 // apply an empty LinearOperator exemplar.
1050 matrices[i] =
1051 linear_operator<VectorType>(LinearOperator<VectorType>(),
1053 smoothers[i].initialize(Utilities::get_underlying_value(m[i]), data);
1054 }
1055}
1056
1057
1058
1059template <typename MatrixType, typename PreconditionerType, typename VectorType>
1060template <typename MatrixType2>
1061inline void
1064{
1065 const unsigned int min = m.min_level();
1066 const unsigned int max = m.max_level();
1067
1068 matrices.resize(min, max);
1069 smoothers.resize(min, max);
1070
1071 for (unsigned int i = min; i <= max; ++i)
1072 {
1073 // Workaround: Unfortunately, not every "m[i]" object has a rich
1074 // enough interface to populate reinit_(domain|range)_vector. Thus,
1075 // apply an empty LinearOperator exemplar.
1076 matrices[i] =
1077 linear_operator<VectorType>(LinearOperator<VectorType>(),
1079 }
1080}
1081
1082
1083
1084template <typename MatrixType, typename PreconditionerType, typename VectorType>
1085template <typename MatrixType2, typename DataType>
1086inline void
1090{
1091 const unsigned int min = m.min_level();
1092 const unsigned int max = m.max_level();
1093
1094 Assert(data.min_level() == min, ExcDimensionMismatch(data.min_level(), min));
1095 Assert(data.max_level() == max, ExcDimensionMismatch(data.max_level(), max));
1096
1097 matrices.resize(min, max);
1098 smoothers.resize(min, max);
1099
1100 for (unsigned int i = min; i <= max; ++i)
1101 {
1102 // Workaround: Unfortunately, not every "m[i]" object has a rich
1103 // enough interface to populate reinit_(domain|range)_vector. Thus,
1104 // apply an empty LinearOperator exemplar.
1105 matrices[i] =
1106 linear_operator<VectorType>(LinearOperator<VectorType>(),
1108 smoothers[i].initialize(Utilities::get_underlying_value(m[i]), data[i]);
1109 }
1110}
1111
1112
1113
1114template <typename MatrixType, typename PreconditionerType, typename VectorType>
1115template <typename MatrixType2, typename DataType>
1116inline void
1119 const DataType &data,
1120 const unsigned int row,
1121 const unsigned int col)
1122{
1123 const unsigned int min = m.min_level();
1124 const unsigned int max = m.max_level();
1125
1126 matrices.resize(min, max);
1127 smoothers.resize(min, max);
1128
1129 for (unsigned int i = min; i <= max; ++i)
1130 {
1131 matrices[i] = &(m[i].block(row, col));
1132 smoothers[i].initialize(m[i].block(row, col), data);
1133 }
1134}
1135
1136
1137
1138template <typename MatrixType, typename PreconditionerType, typename VectorType>
1139template <typename MatrixType2, typename DataType>
1140inline void
1144 const unsigned int row,
1145 const unsigned int col)
1146{
1147 const unsigned int min = m.min_level();
1148 const unsigned int max = m.max_level();
1149
1150 Assert(data.min_level() == min, ExcDimensionMismatch(data.min_level(), min));
1151 Assert(data.max_level() == max, ExcDimensionMismatch(data.max_level(), max));
1152
1153 matrices.resize(min, max);
1154 smoothers.resize(min, max);
1155
1156 for (unsigned int i = min; i <= max; ++i)
1157 {
1158 matrices[i] = &(m[i].block(row, col));
1159 smoothers[i].initialize(m[i].block(row, col), data[i]);
1160 }
1161}
1162
1163
1164
1165template <typename MatrixType, typename PreconditionerType, typename VectorType>
1166inline void
1168 const unsigned int level,
1169 VectorType &u,
1170 const VectorType &rhs) const
1171{
1172 unsigned int maxlevel = matrices.max_level();
1173 unsigned int steps2 = this->steps;
1174
1175 if (this->variable)
1176 steps2 *= (1 << (maxlevel - level));
1177
1178 typename VectorMemory<VectorType>::Pointer r(this->vector_memory);
1179 typename VectorMemory<VectorType>::Pointer d(this->vector_memory);
1180
1181 r->reinit(u, true);
1182 d->reinit(u, true);
1183
1184 bool T = this->transpose;
1185 if (this->symmetric && (steps2 % 2 == 0))
1186 T = false;
1187 if (this->debug > 0)
1188 deallog << 'S' << level << ' ';
1189
1190 for (unsigned int i = 0; i < steps2; ++i)
1191 {
1192 if (T)
1193 {
1194 if (this->debug > 0)
1195 deallog << 'T';
1196 matrices[level].Tvmult(*r, u);
1197 r->sadd(-1., 1., rhs);
1198 if (this->debug > 2)
1199 deallog << ' ' << r->l2_norm() << ' ';
1200 smoothers[level].Tvmult(*d, *r);
1201 if (this->debug > 1)
1202 deallog << ' ' << d->l2_norm() << ' ';
1203 }
1204 else
1205 {
1206 if (this->debug > 0)
1207 deallog << 'N';
1208 matrices[level].vmult(*r, u);
1209 r->sadd(-1., rhs);
1210 if (this->debug > 2)
1211 deallog << ' ' << r->l2_norm() << ' ';
1212 smoothers[level].vmult(*d, *r);
1213 if (this->debug > 1)
1214 deallog << ' ' << d->l2_norm() << ' ';
1215 }
1216 u += *d;
1217 if (this->symmetric)
1218 T = !T;
1219 }
1220 if (this->debug > 0)
1221 deallog << std::endl;
1222}
1223
1224
1225
1226template <typename MatrixType, typename PreconditionerType, typename VectorType>
1227inline void
1229 const unsigned int level,
1230 VectorType &u,
1231 const VectorType &rhs) const
1232{
1233 unsigned int maxlevel = matrices.max_level();
1234 unsigned int steps2 = this->steps;
1235
1236 if (this->variable)
1237 steps2 *= (1 << (maxlevel - level));
1238
1239 bool T = this->transpose;
1240 if (this->symmetric && (steps2 % 2 == 0))
1241 T = false;
1242 if (this->debug > 0)
1243 deallog << 'S' << level << ' ';
1244
1245 // first step where we overwrite the result
1246 if (this->debug > 2)
1247 deallog << ' ' << rhs.l2_norm() << ' ';
1248 if (this->debug > 0)
1249 deallog << (T ? 'T' : 'N');
1250 if (T)
1251 smoothers[level].Tvmult(u, rhs);
1252 else
1253 smoothers[level].vmult(u, rhs);
1254 if (this->debug > 1)
1255 deallog << ' ' << u.l2_norm() << ' ';
1256 if (this->symmetric)
1257 T = !T;
1258
1259 typename VectorMemory<VectorType>::Pointer r(this->vector_memory);
1260 typename VectorMemory<VectorType>::Pointer d(this->vector_memory);
1261
1262 if (steps2 > 1)
1263 {
1264 r->reinit(u, true);
1265 d->reinit(u, true);
1266 }
1267
1268 for (unsigned int i = 1; i < steps2; ++i)
1269 {
1270 if (T)
1271 {
1272 if (this->debug > 0)
1273 deallog << 'T';
1274 matrices[level].Tvmult(*r, u);
1275 r->sadd(-1., 1., rhs);
1276 if (this->debug > 2)
1277 deallog << ' ' << r->l2_norm() << ' ';
1278 smoothers[level].Tvmult(*d, *r);
1279 if (this->debug > 1)
1280 deallog << ' ' << d->l2_norm() << ' ';
1281 }
1282 else
1283 {
1284 if (this->debug > 0)
1285 deallog << 'N';
1286 matrices[level].vmult(*r, u);
1287 r->sadd(-1., rhs);
1288 if (this->debug > 2)
1289 deallog << ' ' << r->l2_norm() << ' ';
1290 smoothers[level].vmult(*d, *r);
1291 if (this->debug > 1)
1292 deallog << ' ' << d->l2_norm() << ' ';
1293 }
1294 u += *d;
1295 if (this->symmetric)
1296 T = !T;
1297 }
1298 if (this->debug > 0)
1299 deallog << std::endl;
1300}
1301
1302
1303
1304template <typename MatrixType, typename PreconditionerType, typename VectorType>
1305inline std::size_t
1307 memory_consumption() const
1308{
1309 return sizeof(*this) + matrices.memory_consumption() +
1310 smoothers.memory_consumption() +
1311 this->vector_memory.memory_consumption();
1312}
1313
1314
1315#endif // DOXYGEN
1316
1318
1319#endif
std::size_t memory_consumption() const
std::size_t memory_consumption() const
unsigned int max_level() const
unsigned int min_level() const
virtual void clear()
virtual void smooth(const unsigned int level, VectorType &u, const VectorType &rhs) const
std::size_t memory_consumption() const
MGLevelObject< PreconditionerType > smoothers
void initialize(const MGLevelObject< MatrixType2 > &matrices, const MGLevelObject< DataType > &additional_data)
void initialize_matrices(const MGLevelObject< MatrixType2 > &matrices)
virtual void smooth(const unsigned int level, VectorType &u, const VectorType &rhs) const override
void initialize(const MGLevelObject< MatrixType2 > &matrices, const typename PreconditionerType::AdditionalData &additional_data=typename PreconditionerType::AdditionalData())
virtual void apply(const unsigned int level, VectorType &u, const VectorType &rhs) const override
void clear() override
void initialize(const MGLevelObject< MatrixType2 > &matrices, const DataType &additional_data, const unsigned int block_row, const unsigned int block_col)
MGLevelObject< LinearOperator< VectorType > > matrices
MGSmootherPrecondition(const unsigned int steps=1, const bool variable=false, const bool symmetric=false, const bool transpose=false)
void initialize(const MGLevelObject< MatrixType2 > &matrices, const MGLevelObject< DataType > &additional_data, const unsigned int block_row, const unsigned int block_col)
virtual void smooth(const unsigned int level, VectorType &u, const VectorType &rhs) const
MGLevelObject< LinearOperator< VectorType > > matrices
void initialize(const MGLevelObject< MatrixType2 > &matrices, const DataType &additional_data, const unsigned int block_row, const unsigned int block_col)
void initialize(const MGLevelObject< MatrixType2 > &matrices, const MGLevelObject< DataType > &additional_data)
void initialize(const MGLevelObject< MatrixType2 > &matrices, const MGLevelObject< DataType > &additional_data, const unsigned int block_row, const unsigned int block_col)
std::size_t memory_consumption() const
virtual void apply(const unsigned int level, VectorType &u, const VectorType &rhs) const
MGLevelObject< RelaxationType > smoothers
MGSmootherRelaxation(const unsigned int steps=1, const bool variable=false, const bool symmetric=false, const bool transpose=false)
void initialize(const MGLevelObject< MatrixType2 > &matrices, const typename RelaxationType::AdditionalData &additional_data=typename RelaxationType::AdditionalData())
GrowingVectorMemory< VectorType > vector_memory
MGSmoother(const unsigned int steps=1, const bool variable=false, const bool symmetric=false, const bool transpose=false)
void set_debug(const unsigned int level)
void set_steps(const unsigned int)
unsigned int debug
void set_symmetric(const bool)
void set_transpose(const bool)
unsigned int steps
void set_variable(const bool)
std::size_t memory_consumption() const
void initialize(const MGLevelObject< MatrixType2 > &matrices, const MGLevelObject< DataType > &additional_data)
SmootherRelaxation(const unsigned int steps=1, const bool variable=false, const bool symmetric=false, const bool transpose=false)
virtual void smooth(const unsigned int level, VectorType &u, const VectorType &rhs) const override
virtual void apply(const unsigned int level, VectorType &u, const VectorType &rhs) const override
void clear() override
void initialize(const MGLevelObject< MatrixType2 > &matrices, const typename RelaxationType::AdditionalData &additional_data=typename RelaxationType::AdditionalData())
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:40
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:41
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
unsigned int level
Definition grid_out.cc:4632
#define Assert(cond, exc)
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
LogStream deallog
Definition logstream.cc:38
std::vector< index_type > data
Definition mpi.cc:744
@ symmetric
Matrix is symmetric.
constexpr char T
Tpetra::Vector< Number, LO, GO, NodeType< MemorySpace > > VectorType
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
T & get_underlying_value(T &p)
Definition utilities.h:1640
Definition mg.h:81
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)