Reference documentation for deal.II version 9.4.1
\(\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.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_h
17#define dealii_precondition_h
18
19#include <deal.II/base/config.h>
20
21#include <deal.II/base/cuda_size.h>
22#include <deal.II/base/memory_space.h>
23#include <deal.II/base/parallel.h>
24#include <deal.II/base/smartpointer.h>
25#include <deal.II/base/template_constraints.h>
26#include <deal.II/base/thread_management.h>
27#include <deal.II/base/utilities.h>
28
29#include <deal.II/lac/affine_constraints.h>
30#include <deal.II/lac/block_vector_base.h>
31#include <deal.II/lac/diagonal_matrix.h>
32#include <deal.II/lac/identity_matrix.h>
33#include <deal.II/lac/solver_cg.h>
34#include <deal.II/lac/vector_memory.h>
35
36DEAL_II_NAMESPACE_OPEN
37
38// forward declarations
39#ifndef DOXYGEN
40template <typename number>
41class Vector;
42template <typename number>
43class SparseMatrix;
44namespace LinearAlgebra
45{
46 namespace distributed
47 {
48 template <typename, typename>
49 class Vector;
50 template <typename>
51 class BlockVector;
52 } // namespace distributed
53} // namespace LinearAlgebra
54#endif
55
56
81class PreconditionIdentity : public Subscriptor
82{
83public:
87 using size_type = types::global_dof_index;
88
93 struct AdditionalData
94 {
98 AdditionalData() = default;
99 };
100
104 PreconditionIdentity();
105
110 template <typename MatrixType>
111 void
112 initialize(const MatrixType & matrix,
113 const AdditionalData &additional_data = AdditionalData());
114
118 template <class VectorType>
119 void
120 vmult(VectorType &, const VectorType &) const;
121
126 template <class VectorType>
127 void
128 Tvmult(VectorType &, const VectorType &) const;
129
133 template <class VectorType>
134 void
135 vmult_add(VectorType &, const VectorType &) const;
136
141 template <class VectorType>
142 void
143 Tvmult_add(VectorType &, const VectorType &) const;
144
149 void
150 clear();
151
159 size_type
160 m() const;
161
169 size_type
170 n() const;
171
172private:
176 size_type n_rows;
177
181 size_type n_columns;
182};
183
184
185
196class PreconditionRichardson : public Subscriptor
197{
198public:
202 using size_type = types::global_dof_index;
203
207 class AdditionalData
208 {
209 public:
214 AdditionalData(const double relaxation = 1.);
215
219 double relaxation;
220 };
221
226 PreconditionRichardson();
227
231 void
232 initialize(const AdditionalData &parameters);
233
239 template <typename MatrixType>
240 void
241 initialize(const MatrixType &matrix, const AdditionalData &parameters);
242
246 template <class VectorType>
247 void
248 vmult(VectorType &, const VectorType &) const;
249
254 template <class VectorType>
255 void
256 Tvmult(VectorType &, const VectorType &) const;
260 template <class VectorType>
261 void
262 vmult_add(VectorType &, const VectorType &) const;
263
268 template <class VectorType>
269 void
270 Tvmult_add(VectorType &, const VectorType &) const;
271
276 void
277 clear()
278 {}
279
287 size_type
288 m() const;
289
297 size_type
298 n() const;
299
300private:
304 double relaxation;
305
309 size_type n_rows;
310
314 size_type n_columns;
315};
316
317
318
358template <typename MatrixType = SparseMatrix<double>,
359 class VectorType = Vector<double>>
360class PreconditionUseMatrix : public Subscriptor
361{
362public:
366 using function_ptr = void (MatrixType::*)(VectorType &,
367 const VectorType &) const;
368
374 PreconditionUseMatrix(const MatrixType &M, const function_ptr method);
375
380 void
381 vmult(VectorType &dst, const VectorType &src) const;
382
383private:
387 const MatrixType &matrix;
388
392 const function_ptr precondition;
393};
394
395
396
402template <typename MatrixType = SparseMatrix<double>,
403 typename PreconditionerType = IdentityMatrix>
404class PreconditionRelaxation : public Subscriptor
405{
406public:
410 using size_type = types::global_dof_index;
411
415 class AdditionalData
416 {
417 public:
421 AdditionalData(const double relaxation = 1.,
422 const unsigned int n_iterations = 1);
423
427 double relaxation;
428
432 unsigned int n_iterations;
433
434
435 /*
436 * Preconditioner.
437 */
438 std::shared_ptr<PreconditionerType> preconditioner;
439 };
440
446 void
447 initialize(const MatrixType & A,
448 const AdditionalData &parameters = AdditionalData());
449
453 void
454 clear();
455
460 size_type
461 m() const;
462
467 size_type
468 n() const;
469
473 template <class VectorType>
474 void
475 vmult(VectorType &, const VectorType &) const;
476
481 template <class VectorType>
482 void
483 Tvmult(VectorType &, const VectorType &) const;
484
488 template <class VectorType>
489 void
490 step(VectorType &x, const VectorType &rhs) const;
491
495 template <class VectorType>
496 void
497 Tstep(VectorType &x, const VectorType &rhs) const;
498
499protected:
503 SmartPointer<const MatrixType, PreconditionRelaxation<MatrixType>> A;
504
508 double relaxation;
509
513 unsigned int n_iterations;
514
515 /*
516 * Preconditioner.
517 */
518 std::shared_ptr<PreconditionerType> preconditioner;
519};
520
521
522
523#ifndef DOXYGEN
524
525namespace internal
526{
527 namespace PreconditionRelaxation
528 {
529 template <typename T, typename VectorType>
530 using Tvmult_t = decltype(
531 std::declval<T const>().Tvmult(std::declval<VectorType &>(),
532 std::declval<const VectorType &>()));
533
534 template <typename T, typename VectorType>
535 constexpr bool has_Tvmult = is_supported_operation<Tvmult_t, T, VectorType>;
536
537 template <typename T, typename VectorType>
538 using step_t = decltype(
539 std::declval<T const>().step(std::declval<VectorType &>(),
540 std::declval<const VectorType &>()));
541
542 template <typename T, typename VectorType>
543 constexpr bool has_step = is_supported_operation<step_t, T, VectorType>;
544
545 template <typename T, typename VectorType>
546 using step_omega_t =
547 decltype(std::declval<T const>().step(std::declval<VectorType &>(),
548 std::declval<const VectorType &>(),
549 std::declval<const double>()));
550
551 template <typename T, typename VectorType>
552 constexpr bool has_step_omega =
553 is_supported_operation<step_omega_t, T, VectorType>;
554
555 template <typename T, typename VectorType>
556 using Tstep_t = decltype(
557 std::declval<T const>().Tstep(std::declval<VectorType &>(),
558 std::declval<const VectorType &>()));
559
560 template <typename T, typename VectorType>
561 constexpr bool has_Tstep = is_supported_operation<Tstep_t, T, VectorType>;
562
563 template <typename T, typename VectorType>
564 using Tstep_omega_t =
565 decltype(std::declval<T const>().Tstep(std::declval<VectorType &>(),
566 std::declval<const VectorType &>(),
567 std::declval<const double>()));
568
569 template <typename T, typename VectorType>
570 constexpr bool has_Tstep_omega =
571 is_supported_operation<Tstep_omega_t, T, VectorType>;
572
573 template <typename T, typename VectorType>
574 using jacobi_step_t = decltype(
575 std::declval<T const>().Jacobi_step(std::declval<VectorType &>(),
576 std::declval<const VectorType &>(),
577 std::declval<const double>()));
578
579 template <typename T, typename VectorType>
580 constexpr bool has_jacobi_step =
581 is_supported_operation<jacobi_step_t, T, VectorType>;
582
583 template <typename T, typename VectorType>
584 using SOR_step_t = decltype(
585 std::declval<T const>().SOR_step(std::declval<VectorType &>(),
586 std::declval<const VectorType &>(),
587 std::declval<const double>()));
588
589 template <typename T, typename VectorType>
590 constexpr bool has_SOR_step =
591 is_supported_operation<SOR_step_t, T, VectorType>;
592
593 template <typename T, typename VectorType>
594 using SSOR_step_t = decltype(
595 std::declval<T const>().SSOR_step(std::declval<VectorType &>(),
596 std::declval<const VectorType &>(),
597 std::declval<const double>()));
598
599 template <typename T, typename VectorType>
600 constexpr bool has_SSOR_step =
601 is_supported_operation<SSOR_step_t, T, VectorType>;
602
603 template <typename MatrixType>
604 class PreconditionJacobiImpl
605 {
606 public:
607 PreconditionJacobiImpl(const MatrixType &A, const double relaxation)
608 : A(&A)
609 , relaxation(relaxation)
610 {}
611
612 template <typename VectorType>
613 void
614 vmult(VectorType &dst, const VectorType &src) const
615 {
616 this->A->precondition_Jacobi(dst, src, this->relaxation);
617 }
618
619 template <typename VectorType>
620 void
621 Tvmult(VectorType &dst, const VectorType &src) const
622 {
623 // call vmult, since preconditioner is symmetrical
624 this->vmult(dst, src);
625 }
626
627 template <typename VectorType,
628 typename std::enable_if<has_jacobi_step<MatrixType, VectorType>,
629 MatrixType>::type * = nullptr>
630 void
631 step(VectorType &dst, const VectorType &src) const
632 {
633 this->A->Jacobi_step(dst, src, this->relaxation);
634 }
635
636 template <
637 typename VectorType,
638 typename std::enable_if<!has_jacobi_step<MatrixType, VectorType>,
639 MatrixType>::type * = nullptr>
640 void
641 step(VectorType &, const VectorType &) const
642 {
643 AssertThrow(false,
644 ExcMessage(
645 "Matrix A does not provide a Jacobi_step() function!"));
646 }
647
648 template <typename VectorType>
649 void
650 Tstep(VectorType &dst, const VectorType &src) const
651 {
652 // call step, since preconditioner is symmetrical
653 this->step(dst, src);
654 }
655
656 private:
657 const SmartPointer<const MatrixType> A;
658 const double relaxation;
659 };
660
661 template <typename MatrixType>
662 class PreconditionSORImpl
663 {
664 public:
665 PreconditionSORImpl(const MatrixType &A, const double relaxation)
666 : A(&A)
667 , relaxation(relaxation)
668 {}
669
670 template <typename VectorType>
671 void
672 vmult(VectorType &dst, const VectorType &src) const
673 {
674 this->A->precondition_SOR(dst, src, this->relaxation);
675 }
676
677 template <typename VectorType>
678 void
679 Tvmult(VectorType &dst, const VectorType &src) const
680 {
681 this->A->precondition_TSOR(dst, src, this->relaxation);
682 }
683
684 template <typename VectorType,
685 typename std::enable_if<has_SOR_step<MatrixType, VectorType>,
686 MatrixType>::type * = nullptr>
687 void
688 step(VectorType &dst, const VectorType &src) const
689 {
690 this->A->SOR_step(dst, src, this->relaxation);
691 }
692
693 template <typename VectorType,
694 typename std::enable_if<!has_SOR_step<MatrixType, VectorType>,
695 MatrixType>::type * = nullptr>
696 void
697 step(VectorType &, const VectorType &) const
698 {
699 AssertThrow(false,
700 ExcMessage(
701 "Matrix A does not provide a SOR_step() function!"));
702 }
703
704 template <typename VectorType,
705 typename std::enable_if<has_SOR_step<MatrixType, VectorType>,
706 MatrixType>::type * = nullptr>
707 void
708 Tstep(VectorType &dst, const VectorType &src) const
709 {
710 this->A->TSOR_step(dst, src, this->relaxation);
711 }
712
713 template <typename VectorType,
714 typename std::enable_if<!has_SOR_step<MatrixType, VectorType>,
715 MatrixType>::type * = nullptr>
716 void
717 Tstep(VectorType &, const VectorType &) const
718 {
719 AssertThrow(false,
720 ExcMessage(
721 "Matrix A does not provide a TSOR_step() function!"));
722 }
723
724 private:
725 const SmartPointer<const MatrixType> A;
726 const double relaxation;
727 };
728
729 template <typename MatrixType>
730 class PreconditionSSORImpl
731 {
732 public:
733 using size_type = typename MatrixType::size_type;
734
735 PreconditionSSORImpl(const MatrixType &A, const double relaxation)
736 : A(&A)
737 , relaxation(relaxation)
738 {
739 // in case we have a SparseMatrix class, we can extract information
740 // about the diagonal.
741 const SparseMatrix<typename MatrixType::value_type> *mat =
742 dynamic_cast<const SparseMatrix<typename MatrixType::value_type> *>(
743 &*this->A);
744
745 // calculate the positions first after the diagonal.
746 if (mat != nullptr)
747 {
748 const size_type n = this->A->n();
749 pos_right_of_diagonal.resize(n, static_cast<std::size_t>(-1));
750 for (size_type row = 0; row < n; ++row)
751 {
752 // find the first element in this line which is on the right of
753 // the diagonal. we need to precondition with the elements on
754 // the left only. note: the first entry in each line denotes the
755 // diagonal element, which we need not check.
756 typename SparseMatrix<
757 typename MatrixType::value_type>::const_iterator it =
758 mat->begin(row) + 1;
759 for (; it < mat->end(row); ++it)
760 if (it->column() > row)
761 break;
762 pos_right_of_diagonal[row] = it - mat->begin();
763 }
764 }
765 }
766
767 template <typename VectorType>
768 void
769 vmult(VectorType &dst, const VectorType &src) const
770 {
771 this->A->precondition_SSOR(dst,
772 src,
773 this->relaxation,
774 pos_right_of_diagonal);
775 }
776
777 template <typename VectorType>
778 void
779 Tvmult(VectorType &dst, const VectorType &src) const
780 {
781 this->A->precondition_SSOR(dst,
782 src,
783 this->relaxation,
784 pos_right_of_diagonal);
785 }
786
787 template <typename VectorType,
788 typename std::enable_if<has_SSOR_step<MatrixType, VectorType>,
789 MatrixType>::type * = nullptr>
790 void
791 step(VectorType &dst, const VectorType &src) const
792 {
793 this->A->SSOR_step(dst, src, this->relaxation);
794 }
795
796 template <typename VectorType,
797 typename std::enable_if<!has_SSOR_step<MatrixType, VectorType>,
798 MatrixType>::type * = nullptr>
799 void
800 step(VectorType &, const VectorType &) const
801 {
802 AssertThrow(false,
803 ExcMessage(
804 "Matrix A does not provide a SSOR_step() function!"));
805 }
806
807 template <typename VectorType>
808 void
809 Tstep(VectorType &dst, const VectorType &src) const
810 {
811 // call step, since preconditioner is symmetrical
812 this->step(dst, src);
813 }
814
815 private:
816 const SmartPointer<const MatrixType> A;
817 const double relaxation;
818
823 std::vector<std::size_t> pos_right_of_diagonal;
824 };
825
826 template <typename MatrixType>
827 class PreconditionPSORImpl
828 {
829 public:
830 using size_type = typename MatrixType::size_type;
831
832 PreconditionPSORImpl(const MatrixType & A,
833 const double relaxation,
834 const std::vector<size_type> &permutation,
835 const std::vector<size_type> &inverse_permutation)
836 : A(&A)
837 , relaxation(relaxation)
838 , permutation(permutation)
839 , inverse_permutation(inverse_permutation)
840 {}
841
842 template <typename VectorType>
843 void
844 vmult(VectorType &dst, const VectorType &src) const
845 {
846 dst = src;
847 this->A->PSOR(dst, permutation, inverse_permutation, this->relaxation);
848 }
849
850 template <typename VectorType>
851 void
852 Tvmult(VectorType &dst, const VectorType &src) const
853 {
854 dst = src;
855 this->A->TPSOR(dst, permutation, inverse_permutation, this->relaxation);
856 }
857
858 private:
859 const SmartPointer<const MatrixType> A;
860 const double relaxation;
861
862 const std::vector<size_type> &permutation;
863 const std::vector<size_type> &inverse_permutation;
864 };
865
866 template <
867 typename MatrixType,
868 typename PreconditionerType,
869 typename VectorType,
870 typename std::enable_if<has_step_omega<PreconditionerType, VectorType>,
871 PreconditionerType>::type * = nullptr>
872 void
873 step(const MatrixType &,
874 const PreconditionerType &preconditioner,
875 VectorType & dst,
876 const VectorType & src,
877 const double relaxation,
878 VectorType &,
879 VectorType &)
880 {
881 preconditioner.step(dst, src, relaxation);
882 }
883
884 template <
885 typename MatrixType,
886 typename PreconditionerType,
887 typename VectorType,
888 typename std::enable_if<!has_step_omega<PreconditionerType, VectorType> &&
889 has_step<PreconditionerType, VectorType>,
890 PreconditionerType>::type * = nullptr>
891 void
892 step(const MatrixType &,
893 const PreconditionerType &preconditioner,
894 VectorType & dst,
895 const VectorType & src,
896 const double relaxation,
897 VectorType &,
898 VectorType &)
899 {
900 Assert(relaxation == 1.0, ExcInternalError());
901
902 (void)relaxation;
903
904 preconditioner.step(dst, src);
905 }
906
907 template <
908 typename MatrixType,
909 typename PreconditionerType,
910 typename VectorType,
911 typename std::enable_if<!has_step_omega<PreconditionerType, VectorType> &&
912 !has_step<PreconditionerType, VectorType>,
913 PreconditionerType>::type * = nullptr>
914 void
915 step(const MatrixType & A,
916 const PreconditionerType &preconditioner,
917 VectorType & dst,
918 const VectorType & src,
919 const double relaxation,
920 VectorType & residual,
921 VectorType & tmp)
922 {
923 residual.reinit(dst, true);
924 tmp.reinit(dst, true);
925
926 A.vmult(residual, dst);
927 residual.sadd(-1.0, 1.0, src);
928
929 preconditioner.vmult(tmp, residual);
930 dst.add(relaxation, tmp);
931 }
932
933 template <
934 typename MatrixType,
935 typename PreconditionerType,
936 typename VectorType,
937 typename std::enable_if<has_Tstep_omega<PreconditionerType, VectorType>,
938 PreconditionerType>::type * = nullptr>
939 void
940 Tstep(const MatrixType &,
941 const PreconditionerType &preconditioner,
942 VectorType & dst,
943 const VectorType & src,
944 const double relaxation,
945 VectorType &,
946 VectorType &)
947 {
948 preconditioner.Tstep(dst, src, relaxation);
949 }
950
951 template <typename MatrixType,
952 typename PreconditionerType,
953 typename VectorType,
954 typename std::enable_if<
955 !has_Tstep_omega<PreconditionerType, VectorType> &&
956 has_Tstep<PreconditionerType, VectorType>,
957 PreconditionerType>::type * = nullptr>
958 void
959 Tstep(const MatrixType &,
960 const PreconditionerType &preconditioner,
961 VectorType & dst,
962 const VectorType & src,
963 const double relaxation,
964 VectorType &,
965 VectorType &)
966 {
967 Assert(relaxation == 1.0, ExcInternalError());
968
969 (void)relaxation;
970
971 preconditioner.Tstep(dst, src);
972 }
973
974 template <typename MatrixType,
975 typename VectorType,
976 typename std::enable_if<has_Tvmult<MatrixType, VectorType>,
977 MatrixType>::type * = nullptr>
978 void
979 Tvmult(const MatrixType &A, VectorType &dst, const VectorType &src)
980 {
981 A.Tvmult(dst, src);
982 }
983
984 template <typename MatrixType,
985 typename VectorType,
986 typename std::enable_if<!has_Tvmult<MatrixType, VectorType>,
987 MatrixType>::type * = nullptr>
988 void
989 Tvmult(const MatrixType &, VectorType &, const VectorType &)
990 {
991 AssertThrow(false,
992 ExcMessage("Matrix A does not provide a Tvmult() function!"));
993 }
994
995 template <typename MatrixType,
996 typename PreconditionerType,
997 typename VectorType,
998 typename std::enable_if<
999 !has_Tstep_omega<PreconditionerType, VectorType> &&
1000 !has_Tstep<PreconditionerType, VectorType>,
1001 PreconditionerType>::type * = nullptr>
1002 void
1003 Tstep(const MatrixType & A,
1004 const PreconditionerType &preconditioner,
1005 VectorType & dst,
1006 const VectorType & src,
1007 const double relaxation,
1008 VectorType & residual,
1009 VectorType & tmp)
1010 {
1011 residual.reinit(dst, true);
1012 tmp.reinit(dst, true);
1013
1014 Tvmult(A, residual, dst);
1015 residual.sadd(-1.0, 1.0, src);
1016
1017 Tvmult(preconditioner, tmp, residual);
1018 dst.add(relaxation, tmp);
1019 }
1020
1021 template <typename MatrixType,
1022 typename PreconditionerType,
1023 typename VectorType>
1024 void
1025 step_operations(const MatrixType & A,
1026 const PreconditionerType &preconditioner,
1027 VectorType & dst,
1028 const VectorType & src,
1029 const double relaxation,
1030 VectorType & tmp1,
1031 VectorType & tmp2,
1032 const unsigned int i,
1033 const bool transposed)
1034 {
1035 if (i == 0)
1036 {
1037 if (transposed)
1038 Tvmult(preconditioner, dst, src);
1039 else
1040 preconditioner.vmult(dst, src);
1041
1042 if (relaxation != 1.0)
1043 dst *= relaxation;
1044 }
1045 else
1046 {
1047 if (transposed)
1048 Tstep(A, preconditioner, dst, src, relaxation, tmp1, tmp2);
1049 else
1050 step(A, preconditioner, dst, src, relaxation, tmp1, tmp2);
1051 }
1052 }
1053
1054 template <typename MatrixType,
1055 typename VectorType,
1056 typename std::enable_if<!IsBlockVector<VectorType>::value,
1057 VectorType>::type * = nullptr>
1058 void
1059 step_operations(const MatrixType & A,
1060 const DiagonalMatrix<VectorType> &preconditioner,
1061 VectorType & dst,
1062 const VectorType & src,
1063 const double relaxation,
1064 VectorType & residual,
1065 VectorType &,
1066 const unsigned int i,
1067 const bool transposed)
1068 {
1069 if (i == 0)
1070 {
1071 const auto dst_ptr = dst.begin();
1072 const auto src_ptr = src.begin();
1073 const auto diag_ptr = preconditioner.get_vector().begin();
1074
1075 for (unsigned int i = 0; i < dst.locally_owned_size(); ++i)
1076 dst_ptr[i] = relaxation * src_ptr[i] * diag_ptr[i];
1077 }
1078 else
1079 {
1080 residual.reinit(src, true);
1081
1082 const auto dst_ptr = dst.begin();
1083 const auto src_ptr = src.begin();
1084 const auto residual_ptr = residual.begin();
1085 const auto diag_ptr = preconditioner.get_vector().begin();
1086
1087 if (transposed)
1088 A.Tvmult(residual, dst);
1089 else
1090 A.vmult(residual, dst);
1091
1092 for (unsigned int i = 0; i < dst.locally_owned_size(); ++i)
1093 dst_ptr[i] +=
1094 relaxation * (src_ptr[i] - residual_ptr[i]) * diag_ptr[i];
1095 }
1096 }
1097
1098 } // namespace PreconditionRelaxation
1099} // namespace internal
1100
1101#endif
1102
1103
1104
1131template <typename MatrixType = SparseMatrix<double>>
1132class PreconditionJacobi
1133 : public PreconditionRelaxation<
1134 MatrixType,
1135 internal::PreconditionRelaxation::PreconditionJacobiImpl<MatrixType>>
1136{
1137 using PreconditionerType =
1138 internal::PreconditionRelaxation::PreconditionJacobiImpl<MatrixType>;
1139 using BaseClass = PreconditionRelaxation<MatrixType, PreconditionerType>;
1140
1141public:
1145 using AdditionalData = typename BaseClass::AdditionalData;
1146
1150 void
1151 initialize(const MatrixType & A,
1152 const AdditionalData &parameters = AdditionalData());
1153};
1154
1155
1201template <typename MatrixType = SparseMatrix<double>>
1202class PreconditionSOR
1203 : public PreconditionRelaxation<
1204 MatrixType,
1205 internal::PreconditionRelaxation::PreconditionSORImpl<MatrixType>>
1206{
1207 using PreconditionerType =
1208 internal::PreconditionRelaxation::PreconditionSORImpl<MatrixType>;
1209 using BaseClass = PreconditionRelaxation<MatrixType, PreconditionerType>;
1210
1211public:
1215 using AdditionalData = typename BaseClass::AdditionalData;
1216
1220 void
1221 initialize(const MatrixType & A,
1222 const AdditionalData &parameters = AdditionalData());
1223};
1224
1225
1226
1253template <typename MatrixType = SparseMatrix<double>>
1254class PreconditionSSOR
1255 : public PreconditionRelaxation<
1256 MatrixType,
1257 internal::PreconditionRelaxation::PreconditionSSORImpl<MatrixType>>
1258{
1259 using PreconditionerType =
1260 internal::PreconditionRelaxation::PreconditionSSORImpl<MatrixType>;
1261 using BaseClass = PreconditionRelaxation<MatrixType, PreconditionerType>;
1262
1263public:
1267 using AdditionalData = typename BaseClass::AdditionalData;
1268
1274 void
1275 initialize(const MatrixType & A,
1276 const AdditionalData &parameters = AdditionalData());
1277};
1278
1279
1309template <typename MatrixType = SparseMatrix<double>>
1310class PreconditionPSOR
1311 : public PreconditionRelaxation<
1312 MatrixType,
1313 internal::PreconditionRelaxation::PreconditionPSORImpl<MatrixType>>
1314{
1315 using PreconditionerType =
1316 internal::PreconditionRelaxation::PreconditionPSORImpl<MatrixType>;
1317 using BaseClass = PreconditionRelaxation<MatrixType, PreconditionerType>;
1318
1319public:
1323 using size_type = typename BaseClass::size_type;
1324
1328 class AdditionalData
1329 {
1330 public:
1341 AdditionalData(const std::vector<size_type> &permutation,
1342 const std::vector<size_type> &inverse_permutation,
1343 const typename BaseClass::AdditionalData &parameters =
1344 typename BaseClass::AdditionalData());
1345
1349 const std::vector<size_type> &permutation;
1353 const std::vector<size_type> &inverse_permutation;
1357 typename BaseClass::AdditionalData parameters;
1358 };
1359
1371 void
1372 initialize(const MatrixType & A,
1373 const std::vector<size_type> & permutation,
1374 const std::vector<size_type> & inverse_permutation,
1375 const typename BaseClass::AdditionalData &parameters =
1376 typename BaseClass::AdditionalData());
1377
1388 void
1389 initialize(const MatrixType &A, const AdditionalData &additional_data);
1390};
1391
1392
1393
1591template <typename MatrixType = SparseMatrix<double>,
1592 typename VectorType = Vector<double>,
1593 typename PreconditionerType = DiagonalMatrix<VectorType>>
1594class PreconditionChebyshev : public Subscriptor
1595{
1596public:
1600 using size_type = types::global_dof_index;
1601
1606 struct AdditionalData
1607 {
1612 enum class EigenvalueAlgorithm
1613 {
1619 lanczos,
1628 power_iteration
1629 };
1630
1634 AdditionalData(const unsigned int degree = 1,
1635 const double smoothing_range = 0.,
1636 const unsigned int eig_cg_n_iterations = 8,
1637 const double eig_cg_residual = 1e-2,
1638 const double max_eigenvalue = 1,
1639 const EigenvalueAlgorithm eigenvalue_algorithm =
1640 EigenvalueAlgorithm::lanczos);
1641
1645 AdditionalData &
1646 operator=(const AdditionalData &other_data);
1647
1660 unsigned int degree;
1661
1673 double smoothing_range;
1674
1681 unsigned int eig_cg_n_iterations;
1682
1687 double eig_cg_residual;
1688
1694 double max_eigenvalue;
1695
1700 AffineConstraints<double> constraints;
1701
1705 std::shared_ptr<PreconditionerType> preconditioner;
1706
1710 EigenvalueAlgorithm eigenvalue_algorithm;
1711 };
1712
1713
1717 PreconditionChebyshev();
1718
1730 void
1731 initialize(const MatrixType & matrix,
1732 const AdditionalData &additional_data = AdditionalData());
1733
1738 void
1739 vmult(VectorType &dst, const VectorType &src) const;
1740
1745 void
1746 Tvmult(VectorType &dst, const VectorType &src) const;
1747
1751 void
1752 step(VectorType &dst, const VectorType &src) const;
1753
1757 void
1758 Tstep(VectorType &dst, const VectorType &src) const;
1759
1763 void
1764 clear();
1765
1770 size_type
1771 m() const;
1772
1777 size_type
1778 n() const;
1779
1784 struct EigenvalueInformation
1785 {
1789 double min_eigenvalue_estimate;
1793 double max_eigenvalue_estimate;
1797 unsigned int cg_iterations;
1802 unsigned int degree;
1806 EigenvalueInformation()
1807 : min_eigenvalue_estimate{std::numeric_limits<double>::max()}
1808 , max_eigenvalue_estimate{std::numeric_limits<double>::lowest()}
1809 , cg_iterations{0}
1810 , degree{0}
1811 {}
1812 };
1813
1826 EigenvalueInformation
1827 estimate_eigenvalues(const VectorType &src) const;
1828
1829private:
1833 SmartPointer<
1834 const MatrixType,
1835 PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>>
1836 matrix_ptr;
1837
1841 mutable VectorType solution_old;
1842
1846 mutable VectorType temp_vector1;
1847
1851 mutable VectorType temp_vector2;
1852
1857 AdditionalData data;
1858
1862 double theta;
1863
1868 double delta;
1869
1874 bool eigenvalues_are_initialized;
1875
1880 mutable Threads::Mutex mutex;
1881};
1882
1883
1884
1886/* ---------------------------------- Inline functions ------------------- */
1887
1888#ifndef DOXYGEN
1889
1890inline PreconditionIdentity::PreconditionIdentity()
1891 : n_rows(0)
1892 , n_columns(0)
1893{}
1894
1895template <typename MatrixType>
1896inline void
1897PreconditionIdentity::initialize(const MatrixType &matrix,
1898 const PreconditionIdentity::AdditionalData &)
1899{
1900 n_rows = matrix.m();
1901 n_columns = matrix.n();
1902}
1903
1904
1905template <class VectorType>
1906inline void
1907PreconditionIdentity::vmult(VectorType &dst, const VectorType &src) const
1908{
1909 dst = src;
1910}
1911
1912
1913
1914template <class VectorType>
1915inline void
1916PreconditionIdentity::Tvmult(VectorType &dst, const VectorType &src) const
1917{
1918 dst = src;
1919}
1920
1921template <class VectorType>
1922inline void
1923PreconditionIdentity::vmult_add(VectorType &dst, const VectorType &src) const
1924{
1925 dst += src;
1926}
1927
1928
1929
1930template <class VectorType>
1931inline void
1932PreconditionIdentity::Tvmult_add(VectorType &dst, const VectorType &src) const
1933{
1934 dst += src;
1935}
1936
1937
1938
1939inline void
1940PreconditionIdentity::clear()
1941{}
1942
1943
1944
1945inline PreconditionIdentity::size_type
1946PreconditionIdentity::m() const
1947{
1948 Assert(n_rows != 0, ExcNotInitialized());
1949 return n_rows;
1950}
1951
1952inline PreconditionIdentity::size_type
1953PreconditionIdentity::n() const
1954{
1955 Assert(n_columns != 0, ExcNotInitialized());
1956 return n_columns;
1957}
1958
1959//---------------------------------------------------------------------------
1960
1961inline PreconditionRichardson::AdditionalData::AdditionalData(
1962 const double relaxation)
1963 : relaxation(relaxation)
1964{}
1965
1966
1967inline PreconditionRichardson::PreconditionRichardson()
1968 : relaxation(0)
1969 , n_rows(0)
1970 , n_columns(0)
1971{
1972 AdditionalData add_data;
1973 relaxation = add_data.relaxation;
1974}
1975
1976
1977
1978inline void
1979PreconditionRichardson::initialize(
1980 const PreconditionRichardson::AdditionalData &parameters)
1981{
1982 relaxation = parameters.relaxation;
1983}
1984
1985
1986
1987template <typename MatrixType>
1988inline void
1989PreconditionRichardson::initialize(
1990 const MatrixType & matrix,
1991 const PreconditionRichardson::AdditionalData &parameters)
1992{
1993 relaxation = parameters.relaxation;
1994 n_rows = matrix.m();
1995 n_columns = matrix.n();
1996}
1997
1998
1999
2000template <class VectorType>
2001inline void
2002PreconditionRichardson::vmult(VectorType &dst, const VectorType &src) const
2003{
2004 static_assert(
2005 std::is_same<size_type, typename VectorType::size_type>::value,
2006 "PreconditionRichardson and VectorType must have the same size_type.");
2007
2008 dst.equ(relaxation, src);
2009}
2010
2011
2012
2013template <class VectorType>
2014inline void
2015PreconditionRichardson::Tvmult(VectorType &dst, const VectorType &src) const
2016{
2017 static_assert(
2018 std::is_same<size_type, typename VectorType::size_type>::value,
2019 "PreconditionRichardson and VectorType must have the same size_type.");
2020
2021 dst.equ(relaxation, src);
2022}
2023
2024template <class VectorType>
2025inline void
2026PreconditionRichardson::vmult_add(VectorType &dst, const VectorType &src) const
2027{
2028 static_assert(
2029 std::is_same<size_type, typename VectorType::size_type>::value,
2030 "PreconditionRichardson and VectorType must have the same size_type.");
2031
2032 dst.add(relaxation, src);
2033}
2034
2035
2036
2037template <class VectorType>
2038inline void
2039PreconditionRichardson::Tvmult_add(VectorType &dst, const VectorType &src) const
2040{
2041 static_assert(
2042 std::is_same<size_type, typename VectorType::size_type>::value,
2043 "PreconditionRichardson and VectorType must have the same size_type.");
2044
2045 dst.add(relaxation, src);
2046}
2047
2048inline PreconditionRichardson::size_type
2049PreconditionRichardson::m() const
2050{
2051 Assert(n_rows != 0, ExcNotInitialized());
2052 return n_rows;
2053}
2054
2055inline PreconditionRichardson::size_type
2056PreconditionRichardson::n() const
2057{
2058 Assert(n_columns != 0, ExcNotInitialized());
2059 return n_columns;
2060}
2061
2062//---------------------------------------------------------------------------
2063
2064template <typename MatrixType, typename PreconditionerType>
2065inline void
2066PreconditionRelaxation<MatrixType, PreconditionerType>::initialize(
2067 const MatrixType & rA,
2068 const AdditionalData &parameters)
2069{
2070 A = &rA;
2071 relaxation = parameters.relaxation;
2072
2073 Assert(parameters.preconditioner, ExcNotInitialized());
2074
2075 preconditioner = parameters.preconditioner;
2076 n_iterations = parameters.n_iterations;
2077}
2078
2079
2080template <typename MatrixType, typename PreconditionerType>
2081inline void
2082PreconditionRelaxation<MatrixType, PreconditionerType>::clear()
2083{
2084 A = nullptr;
2085 preconditioner = nullptr;
2086}
2087
2088template <typename MatrixType, typename PreconditionerType>
2089inline
2090 typename PreconditionRelaxation<MatrixType, PreconditionerType>::size_type
2091 PreconditionRelaxation<MatrixType, PreconditionerType>::m() const
2092{
2093 Assert(A != nullptr, ExcNotInitialized());
2094 return A->m();
2095}
2096
2097template <typename MatrixType, typename PreconditionerType>
2098inline
2099 typename PreconditionRelaxation<MatrixType, PreconditionerType>::size_type
2100 PreconditionRelaxation<MatrixType, PreconditionerType>::n() const
2101{
2102 Assert(A != nullptr, ExcNotInitialized());
2103 return A->n();
2104}
2105
2106template <typename MatrixType, typename PreconditionerType>
2107template <class VectorType>
2108inline void
2109PreconditionRelaxation<MatrixType, PreconditionerType>::vmult(
2110 VectorType & dst,
2111 const VectorType &src) const
2112{
2113 Assert(this->A != nullptr, ExcNotInitialized());
2114 Assert(this->preconditioner != nullptr, ExcNotInitialized());
2115
2116 VectorType tmp1, tmp2;
2117
2118 for (unsigned int i = 0; i < n_iterations; ++i)
2119 internal::PreconditionRelaxation::step_operations(
2120 *A, *preconditioner, dst, src, relaxation, tmp1, tmp2, i, false);
2121}
2122
2123template <typename MatrixType, typename PreconditionerType>
2124template <class VectorType>
2125inline void
2126PreconditionRelaxation<MatrixType, PreconditionerType>::Tvmult(
2127 VectorType & dst,
2128 const VectorType &src) const
2129{
2130 Assert(this->A != nullptr, ExcNotInitialized());
2131 Assert(this->preconditioner != nullptr, ExcNotInitialized());
2132
2133 VectorType tmp1, tmp2;
2134
2135 for (unsigned int i = 0; i < n_iterations; ++i)
2136 internal::PreconditionRelaxation::step_operations(
2137 *A, *preconditioner, dst, src, relaxation, tmp1, tmp2, i, true);
2138}
2139
2140template <typename MatrixType, typename PreconditionerType>
2141template <class VectorType>
2142inline void
2143PreconditionRelaxation<MatrixType, PreconditionerType>::step(
2144 VectorType & dst,
2145 const VectorType &src) const
2146{
2147 Assert(this->A != nullptr, ExcNotInitialized());
2148 Assert(this->preconditioner != nullptr, ExcNotInitialized());
2149
2150 VectorType tmp1, tmp2;
2151
2152 for (unsigned int i = 1; i <= n_iterations; ++i)
2153 internal::PreconditionRelaxation::step_operations(
2154 *A, *preconditioner, dst, src, relaxation, tmp1, tmp2, i, false);
2155}
2156
2157template <typename MatrixType, typename PreconditionerType>
2158template <class VectorType>
2159inline void
2160PreconditionRelaxation<MatrixType, PreconditionerType>::Tstep(
2161 VectorType & dst,
2162 const VectorType &src) const
2163{
2164 Assert(this->A != nullptr, ExcNotInitialized());
2165 Assert(this->preconditioner != nullptr, ExcNotInitialized());
2166
2167 VectorType tmp1, tmp2;
2168
2169 for (unsigned int i = 1; i <= n_iterations; ++i)
2170 internal::PreconditionRelaxation::step_operations(
2171 *A, *preconditioner, dst, src, relaxation, tmp1, tmp2, i, true);
2172}
2173
2174//---------------------------------------------------------------------------
2175
2176template <typename MatrixType>
2177inline void
2178PreconditionJacobi<MatrixType>::initialize(const MatrixType & A,
2179 const AdditionalData &parameters_in)
2180{
2181 Assert(parameters_in.preconditioner == nullptr, ExcInternalError());
2182
2183 AdditionalData parameters;
2184 parameters.relaxation = 1.0;
2185 parameters.n_iterations = parameters_in.n_iterations;
2186 parameters.preconditioner =
2187 std::make_shared<PreconditionerType>(A, parameters_in.relaxation);
2188
2189 this->BaseClass::initialize(A, parameters);
2190}
2191
2192//---------------------------------------------------------------------------
2193
2194template <typename MatrixType>
2195inline void
2196PreconditionSOR<MatrixType>::initialize(const MatrixType & A,
2197 const AdditionalData &parameters_in)
2198{
2199 Assert(parameters_in.preconditioner == nullptr, ExcInternalError());
2200
2201 AdditionalData parameters;
2202 parameters.relaxation = 1.0;
2203 parameters.n_iterations = parameters_in.n_iterations;
2204 parameters.preconditioner =
2205 std::make_shared<PreconditionerType>(A, parameters_in.relaxation);
2206
2207 this->BaseClass::initialize(A, parameters);
2208}
2209
2210//---------------------------------------------------------------------------
2211
2212template <typename MatrixType>
2213inline void
2214PreconditionSSOR<MatrixType>::initialize(const MatrixType & A,
2215 const AdditionalData &parameters_in)
2216{
2217 Assert(parameters_in.preconditioner == nullptr, ExcInternalError());
2218
2219 AdditionalData parameters;
2220 parameters.relaxation = 1.0;
2221 parameters.n_iterations = parameters_in.n_iterations;
2222 parameters.preconditioner =
2223 std::make_shared<PreconditionerType>(A, parameters_in.relaxation);
2224
2225 this->BaseClass::initialize(A, parameters);
2226}
2227
2228
2229
2230//---------------------------------------------------------------------------
2231
2232template <typename MatrixType>
2233inline void
2234PreconditionPSOR<MatrixType>::initialize(
2235 const MatrixType & A,
2236 const std::vector<size_type> & p,
2237 const std::vector<size_type> & ip,
2238 const typename BaseClass::AdditionalData &parameters_in)
2239{
2240 Assert(parameters_in.preconditioner == nullptr, ExcInternalError());
2241
2242 typename BaseClass::AdditionalData parameters;
2243 parameters.relaxation = 1.0;
2244 parameters.n_iterations = parameters_in.n_iterations;
2245 parameters.preconditioner =
2246 std::make_shared<PreconditionerType>(A, parameters_in.relaxation, p, ip);
2247
2248 this->BaseClass::initialize(A, parameters);
2249}
2250
2251
2252template <typename MatrixType>
2253inline void
2254PreconditionPSOR<MatrixType>::initialize(const MatrixType & A,
2255 const AdditionalData &additional_data)
2256{
2257 initialize(A,
2258 additional_data.permutation,
2259 additional_data.inverse_permutation,
2260 additional_data.parameters);
2261}
2262
2263template <typename MatrixType>
2264PreconditionPSOR<MatrixType>::AdditionalData::AdditionalData(
2265 const std::vector<size_type> &permutation,
2266 const std::vector<size_type> &inverse_permutation,
2267 const typename PreconditionPSOR<MatrixType>::BaseClass::AdditionalData
2268 &parameters)
2269 : permutation(permutation)
2270 , inverse_permutation(inverse_permutation)
2271 , parameters(parameters)
2272{}
2273
2274
2275//---------------------------------------------------------------------------
2276
2277
2278template <typename MatrixType, class VectorType>
2279PreconditionUseMatrix<MatrixType, VectorType>::PreconditionUseMatrix(
2280 const MatrixType & M,
2281 const function_ptr method)
2282 : matrix(M)
2283 , precondition(method)
2284{}
2285
2286
2287
2288template <typename MatrixType, class VectorType>
2289void
2290PreconditionUseMatrix<MatrixType, VectorType>::vmult(
2291 VectorType & dst,
2292 const VectorType &src) const
2293{
2294 (matrix.*precondition)(dst, src);
2295}
2296
2297//---------------------------------------------------------------------------
2298
2299template <typename MatrixType, typename PreconditionerType>
2300inline PreconditionRelaxation<MatrixType, PreconditionerType>::AdditionalData::
2301 AdditionalData(const double relaxation, const unsigned int n_iterations)
2302 : relaxation(relaxation)
2303 , n_iterations(n_iterations)
2304{}
2305
2306
2307
2308//---------------------------------------------------------------------------
2309
2310namespace internal
2311{
2312 namespace PreconditionChebyshevImplementation
2313 {
2314 // for deal.II vectors, perform updates for Chebyshev preconditioner all
2315 // at once to reduce memory transfer. Here, we select between general
2316 // vectors and deal.II vectors where we expand the loop over the (local)
2317 // size of the vector
2318
2319 // generic part for non-deal.II vectors
2320 template <typename VectorType, typename PreconditionerType>
2321 inline void
2322 vector_updates(const VectorType & rhs,
2323 const PreconditionerType &preconditioner,
2324 const unsigned int iteration_index,
2325 const double factor1,
2326 const double factor2,
2327 VectorType & solution_old,
2328 VectorType & temp_vector1,
2329 VectorType & temp_vector2,
2330 VectorType & solution)
2331 {
2332 if (iteration_index == 0)
2333 {
2334 solution.equ(factor2, rhs);
2335 preconditioner.vmult(solution_old, solution);
2336 }
2337 else if (iteration_index == 1)
2338 {
2339 // compute t = P^{-1} * (b-A*x^{n})
2340 temp_vector1.sadd(-1.0, 1.0, rhs);
2341 preconditioner.vmult(solution_old, temp_vector1);
2342
2343 // compute x^{n+1} = x^{n} + f_1 * x^{n} + f_2 * t
2344 solution_old.sadd(factor2, 1 + factor1, solution);
2345 }
2346 else
2347 {
2348 // compute t = P^{-1} * (b-A*x^{n})
2349 temp_vector1.sadd(-1.0, 1.0, rhs);
2350 preconditioner.vmult(temp_vector2, temp_vector1);
2351
2352 // compute x^{n+1} = x^{n} + f_1 * (x^{n}-x^{n-1}) + f_2 * t
2353 solution_old.sadd(-factor1, factor2, temp_vector2);
2354 solution_old.add(1 + factor1, solution);
2355 }
2356
2357 solution.swap(solution_old);
2358 }
2359
2360 // worker routine for deal.II vectors. Because of vectorization, we need
2361 // to put the loop into an extra structure because the virtual function of
2362 // VectorUpdatesRange prevents the compiler from applying vectorization.
2363 template <typename Number>
2364 struct VectorUpdater
2365 {
2366 VectorUpdater(const Number * rhs,
2367 const Number * matrix_diagonal_inverse,
2368 const unsigned int iteration_index,
2369 const Number factor1,
2370 const Number factor2,
2371 Number * solution_old,
2372 Number * tmp_vector,
2373 Number * solution)
2374 : rhs(rhs)
2375 , matrix_diagonal_inverse(matrix_diagonal_inverse)
2376 , iteration_index(iteration_index)
2377 , factor1(factor1)
2378 , factor2(factor2)
2379 , solution_old(solution_old)
2380 , tmp_vector(tmp_vector)
2381 , solution(solution)
2382 {}
2383
2384 void
2385 apply_to_subrange(const std::size_t begin, const std::size_t end) const
2386 {
2387 // To circumvent a bug in gcc
2388 // (https://gcc.gnu.org/bugzilla/show_bug.cgi?id=63945), we create
2389 // copies of the variables factor1 and factor2 and do not check based on
2390 // factor1.
2391 const Number factor1 = this->factor1;
2392 const Number factor1_plus_1 = 1. + this->factor1;
2393 const Number factor2 = this->factor2;
2394 if (iteration_index == 0)
2395 {
2396 DEAL_II_OPENMP_SIMD_PRAGMA
2397 for (std::size_t i = begin; i < end; ++i)
2398 solution[i] = factor2 * matrix_diagonal_inverse[i] * rhs[i];
2399 }
2400 else if (iteration_index == 1)
2401 {
2402 // x^{n+1} = x^{n} + f_1 * x^{n} + f_2 * P^{-1} * (b-A*x^{n})
2403 DEAL_II_OPENMP_SIMD_PRAGMA
2404 for (std::size_t i = begin; i < end; ++i)
2405 // for efficiency reason, write back to temp_vector that is
2406 // already read (avoid read-for-ownership)
2407 tmp_vector[i] =
2408 factor1_plus_1 * solution[i] +
2409 factor2 * matrix_diagonal_inverse[i] * (rhs[i] - tmp_vector[i]);
2410 }
2411 else
2412 {
2413 // x^{n+1} = x^{n} + f_1 * (x^{n}-x^{n-1})
2414 // + f_2 * P^{-1} * (b-A*x^{n})
2415 DEAL_II_OPENMP_SIMD_PRAGMA
2416 for (std::size_t i = begin; i < end; ++i)
2417 solution_old[i] =
2418 factor1_plus_1 * solution[i] - factor1 * solution_old[i] +
2419 factor2 * matrix_diagonal_inverse[i] * (rhs[i] - tmp_vector[i]);
2420 }
2421 }
2422
2423 const Number * rhs;
2424 const Number * matrix_diagonal_inverse;
2425 const unsigned int iteration_index;
2426 const Number factor1;
2427 const Number factor2;
2428 mutable Number * solution_old;
2429 mutable Number * tmp_vector;
2430 mutable Number * solution;
2431 };
2432
2433 template <typename Number>
2434 struct VectorUpdatesRange : public ::parallel::ParallelForInteger
2435 {
2436 VectorUpdatesRange(const VectorUpdater<Number> &updater,
2437 const std::size_t size)
2438 : updater(updater)
2439 {
2440 if (size < internal::VectorImplementation::minimum_parallel_grain_size)
2441 VectorUpdatesRange::apply_to_subrange(0, size);
2442 else
2443 apply_parallel(
2444 0,
2445 size,
2446 internal::VectorImplementation::minimum_parallel_grain_size);
2447 }
2448
2449 ~VectorUpdatesRange() override = default;
2450
2451 virtual void
2452 apply_to_subrange(const std::size_t begin,
2453 const std::size_t end) const override
2454 {
2455 updater.apply_to_subrange(begin, end);
2456 }
2457
2458 const VectorUpdater<Number> &updater;
2459 };
2460
2461 // selection for diagonal matrix around deal.II vector
2462 template <typename Number>
2463 inline void
2464 vector_updates(const ::Vector<Number> & rhs,
2465 const DiagonalMatrix<::Vector<Number>> &jacobi,
2466 const unsigned int iteration_index,
2467 const double factor1,
2468 const double factor2,
2469 ::Vector<Number> &solution_old,
2470 ::Vector<Number> &temp_vector1,
2471 ::Vector<Number> &,
2472 ::Vector<Number> &solution)
2473 {
2474 VectorUpdater<Number> upd(rhs.begin(),
2475 jacobi.get_vector().begin(),
2476 iteration_index,
2477 factor1,
2478 factor2,
2479 solution_old.begin(),
2480 temp_vector1.begin(),
2481 solution.begin());
2482 VectorUpdatesRange<Number>(upd, rhs.size());
2483
2484 // swap vectors x^{n+1}->x^{n}, given the updates in the function above
2485 if (iteration_index == 0)
2486 {
2487 // nothing to do here because we can immediately write into the
2488 // solution vector without remembering any of the other vectors
2489 }
2490 else if (iteration_index == 1)
2491 {
2492 solution.swap(temp_vector1);
2493 solution_old.swap(temp_vector1);
2494 }
2495 else
2496 solution.swap(solution_old);
2497 }
2498
2499 // selection for diagonal matrix around parallel deal.II vector
2500 template <typename Number>
2501 inline void
2502 vector_updates(
2503 const LinearAlgebra::distributed::Vector<Number, MemorySpace::Host> &rhs,
2504 const DiagonalMatrix<
2505 LinearAlgebra::distributed::Vector<Number, MemorySpace::Host>> &jacobi,
2506 const unsigned int iteration_index,
2507 const double factor1,
2508 const double factor2,
2509 LinearAlgebra::distributed::Vector<Number, MemorySpace::Host>
2510 &solution_old,
2511 LinearAlgebra::distributed::Vector<Number, MemorySpace::Host>
2512 &temp_vector1,
2513 LinearAlgebra::distributed::Vector<Number, MemorySpace::Host> &,
2514 LinearAlgebra::distributed::Vector<Number, MemorySpace::Host> &solution)
2515 {
2516 VectorUpdater<Number> upd(rhs.begin(),
2517 jacobi.get_vector().begin(),
2518 iteration_index,
2519 factor1,
2520 factor2,
2521 solution_old.begin(),
2522 temp_vector1.begin(),
2523 solution.begin());
2524 VectorUpdatesRange<Number>(upd, rhs.locally_owned_size());
2525
2526 // swap vectors x^{n+1}->x^{n}, given the updates in the function above
2527 if (iteration_index == 0)
2528 {
2529 // nothing to do here because we can immediately write into the
2530 // solution vector without remembering any of the other vectors
2531 }
2532 else if (iteration_index == 1)
2533 {
2534 solution.swap(temp_vector1);
2535 solution_old.swap(temp_vector1);
2536 }
2537 else
2538 solution.swap(solution_old);
2539 }
2540
2541 // a helper type-trait that leverage SFINAE to figure out if MatrixType has
2542 // ... MatrixType::vmult(VectorType &, const VectorType&,
2543 // std::function<...>, std::function<...>) const
2544 template <typename MatrixType, typename VectorType>
2545 using vmult_functions_t = decltype(std::declval<MatrixType const>().vmult(
2546 std::declval<VectorType &>(),
2547 std::declval<const VectorType &>(),
2548 std::declval<
2549 const std::function<void(const unsigned int, const unsigned int)> &>(),
2550 std::declval<const std::function<void(const unsigned int,
2551 const unsigned int)> &>()));
2552
2553 template <typename MatrixType,
2554 typename VectorType,
2555 typename PreconditionerType>
2556 constexpr bool has_vmult_with_std_functions =
2557 is_supported_operation<vmult_functions_t, MatrixType, VectorType>
2558 && std::is_same<PreconditionerType, DiagonalMatrix<VectorType>>::value
2559 &&std::is_same<
2560 VectorType,
2561 LinearAlgebra::distributed::Vector<typename VectorType::value_type,
2562 MemorySpace::Host>>::value;
2563
2564 // We need to have a separate declaration for static const members
2565
2566 template <
2567 typename MatrixType,
2568 typename VectorType,
2569 typename PreconditionerType,
2570 typename std::enable_if<!has_vmult_with_std_functions<MatrixType,
2571 VectorType,
2572 PreconditionerType>,
2573 int>::type * = nullptr>
2574 inline void
2575 vmult_and_update(const MatrixType & matrix,
2576 const PreconditionerType &preconditioner,
2577 const VectorType & rhs,
2578 const unsigned int iteration_index,
2579 const double factor1,
2580 const double factor2,
2581 VectorType & solution,
2582 VectorType & solution_old,
2583 VectorType & temp_vector1,
2584 VectorType & temp_vector2)
2585 {
2586 if (iteration_index > 0)
2587 matrix.vmult(temp_vector1, solution);
2588 vector_updates(rhs,
2589 preconditioner,
2590 iteration_index,
2591 factor1,
2592 factor2,
2593 solution_old,
2594 temp_vector1,
2595 temp_vector2,
2596 solution);
2597 }
2598
2599 template <
2600 typename MatrixType,
2601 typename VectorType,
2602 typename PreconditionerType,
2603 typename std::enable_if<has_vmult_with_std_functions<MatrixType,
2604 VectorType,
2605 PreconditionerType>,
2606 int>::type * = nullptr>
2607 inline void
2608 vmult_and_update(const MatrixType & matrix,
2609 const PreconditionerType &preconditioner,
2610 const VectorType & rhs,
2611 const unsigned int iteration_index,
2612 const double factor1,
2613 const double factor2,
2614 VectorType & solution,
2615 VectorType & solution_old,
2616 VectorType & temp_vector1,
2617 VectorType &)
2618 {
2619 using Number = typename VectorType::value_type;
2620 VectorUpdater<Number> updater(rhs.begin(),
2621 preconditioner.get_vector().begin(),
2622 iteration_index,
2623 factor1,
2624 factor2,
2625 solution_old.begin(),
2626 temp_vector1.begin(),
2627 solution.begin());
2628 if (iteration_index > 0)
2629 matrix.vmult(
2630 temp_vector1,
2631 solution,
2632 [&](const unsigned int start_range, const unsigned int end_range) {
2633 // zero 'temp_vector1' before running the vmult
2634 // operation
2635 if (end_range > start_range)
2636 std::memset(temp_vector1.begin() + start_range,
2637 0,
2638 sizeof(Number) * (end_range - start_range));
2639 },
2640 [&](const unsigned int start_range, const unsigned int end_range) {
2641 if (end_range > start_range)
2642 updater.apply_to_subrange(start_range, end_range);
2643 });
2644 else
2645 updater.apply_to_subrange(0U, solution.locally_owned_size());
2646
2647 // swap vectors x^{n+1}->x^{n}, given the updates in the function above
2648 if (iteration_index == 0)
2649 {
2650 // nothing to do here because we can immediately write into the
2651 // solution vector without remembering any of the other vectors
2652 }
2653 else if (iteration_index == 1)
2654 {
2655 solution.swap(temp_vector1);
2656 solution_old.swap(temp_vector1);
2657 }
2658 else
2659 solution.swap(solution_old);
2660 }
2661
2662 template <typename MatrixType, typename PreconditionerType>
2663 inline void
2664 initialize_preconditioner(
2665 const MatrixType & matrix,
2666 std::shared_ptr<PreconditionerType> &preconditioner)
2667 {
2668 (void)matrix;
2669 (void)preconditioner;
2670 AssertThrow(preconditioner.get() != nullptr, ExcNotInitialized());
2671 }
2672
2673 template <typename MatrixType, typename VectorType>
2674 inline void
2675 initialize_preconditioner(
2676 const MatrixType & matrix,
2677 std::shared_ptr<DiagonalMatrix<VectorType>> &preconditioner)
2678 {
2679 if (preconditioner.get() == nullptr || preconditioner->m() != matrix.m())
2680 {
2681 if (preconditioner.get() == nullptr)
2682 preconditioner = std::make_shared<DiagonalMatrix<VectorType>>();
2683
2684 Assert(
2685 preconditioner->m() == 0,
2686 ExcMessage(
2687 "Preconditioner appears to be initialized but not sized correctly"));
2688
2689 // This part only works in serial
2690 if (preconditioner->m() != matrix.m())
2691 {
2692 preconditioner->get_vector().reinit(matrix.m());
2693 for (typename VectorType::size_type i = 0; i < matrix.m(); ++i)
2694 preconditioner->get_vector()(i) = 1. / matrix.el(i, i);
2695 }
2696 }
2697 }
2698
2699 template <typename VectorType>
2700 void
2701 set_initial_guess(VectorType &vector)
2702 {
2703 vector = 1. / std::sqrt(static_cast<double>(vector.size()));
2704 if (vector.locally_owned_elements().is_element(0))
2705 vector(0) = 0.;
2706 }
2707
2708 template <typename Number>
2709 void
2710 set_initial_guess(::Vector<Number> &vector)
2711 {
2712 // Choose a high-frequency mode consisting of numbers between 0 and 1
2713 // that is cheap to compute (cheaper than random numbers) but avoids
2714 // obviously re-occurring numbers in multi-component systems by choosing
2715 // a period of 11
2716 for (unsigned int i = 0; i < vector.size(); ++i)
2717 vector(i) = i % 11;
2718
2719 const Number mean_value = vector.mean_value();
2720 vector.add(-mean_value);
2721 }
2722
2723 template <typename Number>
2724 void
2725 set_initial_guess(
2726 ::LinearAlgebra::distributed::Vector<Number, MemorySpace::Host>
2727 &vector)
2728 {
2729 // Choose a high-frequency mode consisting of numbers between 0 and 1
2730 // that is cheap to compute (cheaper than random numbers) but avoids
2731 // obviously re-occurring numbers in multi-component systems by choosing
2732 // a period of 11.
2733 // Make initial guess robust with respect to number of processors
2734 // by operating on the global index.
2735 types::global_dof_index first_local_range = 0;
2736 if (!vector.locally_owned_elements().is_empty())
2737 first_local_range = vector.locally_owned_elements().nth_index_in_set(0);
2738 for (unsigned int i = 0; i < vector.locally_owned_size(); ++i)
2739 vector.local_element(i) = (i + first_local_range) % 11;
2740
2741 const Number mean_value = vector.mean_value();
2742 vector.add(-mean_value);
2743 }
2744
2745 template <typename Number>
2746 void
2747 set_initial_guess(
2748 ::LinearAlgebra::distributed::BlockVector<Number> &vector)
2749 {
2750 for (unsigned int block = 0; block < vector.n_blocks(); ++block)
2751 set_initial_guess(vector.block(block));
2752 }
2753
2754
2755# ifdef DEAL_II_COMPILER_CUDA_AWARE
2756 template <typename Number>
2757 __global__ void
2758 set_initial_guess_kernel(const types::global_dof_index offset,
2759 const unsigned int locally_owned_size,
2760 Number * values)
2761
2762 {
2763 const unsigned int index = threadIdx.x + blockDim.x * blockIdx.x;
2764 if (index < locally_owned_size)
2765 values[index] = (index + offset) % 11;
2766 }
2767
2768 template <typename Number>
2769 void
2770 set_initial_guess(
2771 ::LinearAlgebra::distributed::Vector<Number, MemorySpace::CUDA>
2772 &vector)
2773 {
2774 // Choose a high-frequency mode consisting of numbers between 0 and 1
2775 // that is cheap to compute (cheaper than random numbers) but avoids
2776 // obviously re-occurring numbers in multi-component systems by choosing
2777 // a period of 11.
2778 // Make initial guess robust with respect to number of processors
2779 // by operating on the global index.
2780 types::global_dof_index first_local_range = 0;
2781 if (!vector.locally_owned_elements().is_empty())
2782 first_local_range = vector.locally_owned_elements().nth_index_in_set(0);
2783
2784 const auto n_local_elements = vector.locally_owned_size();
2785 const int n_blocks =
2786 1 + (n_local_elements - 1) / CUDAWrappers::block_size;
2787 set_initial_guess_kernel<<<n_blocks, CUDAWrappers::block_size>>>(
2788 first_local_range, n_local_elements, vector.get_values());
2789 AssertCudaKernel();
2790
2791 const Number mean_value = vector.mean_value();
2792 vector.add(-mean_value);
2793 }
2794# endif // DEAL_II_COMPILER_CUDA_AWARE
2795
2796 struct EigenvalueTracker
2797 {
2798 public:
2799 void
2800 slot(const std::vector<double> &eigenvalues)
2801 {
2802 values = eigenvalues;
2803 }
2804
2805 std::vector<double> values;
2806 };
2807
2808
2809
2810 template <typename MatrixType,
2811 typename VectorType,
2812 typename PreconditionerType>
2813 double
2814 power_iteration(const MatrixType & matrix,
2815 VectorType & eigenvector,
2816 const PreconditionerType &preconditioner,
2817 const unsigned int n_iterations)
2818 {
2819 double eigenvalue_estimate = 0.;
2820 eigenvector /= eigenvector.l2_norm();
2821 VectorType vector1, vector2;
2822 vector1.reinit(eigenvector, true);
2823 if (!std::is_same<PreconditionerType, PreconditionIdentity>::value)
2824 vector2.reinit(eigenvector, true);
2825 for (unsigned int i = 0; i < n_iterations; ++i)
2826 {
2827 if (!std::is_same<PreconditionerType, PreconditionIdentity>::value)
2828 {
2829 matrix.vmult(vector2, eigenvector);
2830 preconditioner.vmult(vector1, vector2);
2831 }
2832 else
2833 matrix.vmult(vector1, eigenvector);
2834
2835 eigenvalue_estimate = eigenvector * vector1;
2836
2837 vector1 /= vector1.l2_norm();
2838 eigenvector.swap(vector1);
2839 }
2840 return eigenvalue_estimate;
2841 }
2842 } // namespace PreconditionChebyshevImplementation
2843} // namespace internal
2844
2845
2846
2847template <typename MatrixType, class VectorType, typename PreconditionerType>
2848inline PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::
2849 AdditionalData::AdditionalData(const unsigned int degree,
2850 const double smoothing_range,
2851 const unsigned int eig_cg_n_iterations,
2852 const double eig_cg_residual,
2853 const double max_eigenvalue,
2854 const EigenvalueAlgorithm eigenvalue_algorithm)
2855 : degree(degree)
2856 , smoothing_range(smoothing_range)
2857 , eig_cg_n_iterations(eig_cg_n_iterations)
2858 , eig_cg_residual(eig_cg_residual)
2859 , max_eigenvalue(max_eigenvalue)
2860 , eigenvalue_algorithm(eigenvalue_algorithm)
2861{}
2862
2863
2864
2865template <typename MatrixType, class VectorType, typename PreconditionerType>
2866inline typename PreconditionChebyshev<MatrixType,
2867 VectorType,
2868 PreconditionerType>::AdditionalData &
2869PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::
2870 AdditionalData::operator=(const AdditionalData &other_data)
2871{
2872 degree = other_data.degree;
2873 smoothing_range = other_data.smoothing_range;
2874 eig_cg_n_iterations = other_data.eig_cg_n_iterations;
2875 eig_cg_residual = other_data.eig_cg_residual;
2876 max_eigenvalue = other_data.max_eigenvalue;
2877 preconditioner = other_data.preconditioner;
2878 eigenvalue_algorithm = other_data.eigenvalue_algorithm;
2879 constraints.copy_from(other_data.constraints);
2880
2881 return *this;
2882}
2883
2884
2885
2886template <typename MatrixType, typename VectorType, typename PreconditionerType>
2887inline PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::
2888 PreconditionChebyshev()
2889 : theta(1.)
2890 , delta(1.)
2891 , eigenvalues_are_initialized(false)
2892{
2893 static_assert(
2894 std::is_same<size_type, typename VectorType::size_type>::value,
2895 "PreconditionChebyshev and VectorType must have the same size_type.");
2896}
2897
2898
2899
2900template <typename MatrixType, typename VectorType, typename PreconditionerType>
2901inline void
2902PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::initialize(
2903 const MatrixType & matrix,
2904 const AdditionalData &additional_data)
2905{
2906 matrix_ptr = &matrix;
2907 data = additional_data;
2908 Assert(data.degree > 0,
2909 ExcMessage("The degree of the Chebyshev method must be positive."));
2910 internal::PreconditionChebyshevImplementation::initialize_preconditioner(
2911 matrix, data.preconditioner);
2912 eigenvalues_are_initialized = false;
2913}
2914
2915
2916
2917template <typename MatrixType, typename VectorType, typename PreconditionerType>
2918inline void
2919PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::clear()
2920{
2921 eigenvalues_are_initialized = false;
2922 theta = delta = 1.0;
2923 matrix_ptr = nullptr;
2924 {
2925 VectorType empty_vector;
2926 solution_old.reinit(empty_vector);
2927 temp_vector1.reinit(empty_vector);
2928 temp_vector2.reinit(empty_vector);
2929 }
2930 data.preconditioner.reset();
2931}
2932
2933
2934
2935template <typename MatrixType, typename VectorType, typename PreconditionerType>
2936inline typename PreconditionChebyshev<MatrixType,
2937 VectorType,
2938 PreconditionerType>::EigenvalueInformation
2939PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::
2940 estimate_eigenvalues(const VectorType &src) const
2941{
2942 Assert(eigenvalues_are_initialized == false, ExcInternalError());
2943 Assert(data.preconditioner.get() != nullptr, ExcNotInitialized());
2944
2945 PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::
2946 EigenvalueInformation info{};
2947
2948 solution_old.reinit(src);
2949 temp_vector1.reinit(src, true);
2950
2951 if (data.eig_cg_n_iterations > 0)
2952 {
2953 Assert(data.eig_cg_n_iterations > 2,
2954 ExcMessage(
2955 "Need to set at least two iterations to find eigenvalues."));
2956
2957 internal::PreconditionChebyshevImplementation::EigenvalueTracker
2958 eigenvalue_tracker;
2959
2960 // set an initial guess that contains some high-frequency parts (to the
2961 // extent possible without knowing the discretization and the numbering)
2962 // to trigger high eigenvalues according to the external function
2963 internal::PreconditionChebyshevImplementation::set_initial_guess(
2964 temp_vector1);
2965 data.constraints.set_zero(temp_vector1);
2966
2967 if (data.eigenvalue_algorithm ==
2968 AdditionalData::EigenvalueAlgorithm::lanczos)
2969 {
2970 // set a very strict tolerance to force at least two iterations
2971 IterationNumberControl control(data.eig_cg_n_iterations,
2972 1e-10,
2973 false,
2974 false);
2975
2976 SolverCG<VectorType> solver(control);
2977 solver.connect_eigenvalues_slot(
2978 [&eigenvalue_tracker](const std::vector<double> &eigenvalues) {
2979 eigenvalue_tracker.slot(eigenvalues);
2980 });
2981
2982 solver.solve(*matrix_ptr,
2983 solution_old,
2984 temp_vector1,
2985 *data.preconditioner);
2986
2987 info.cg_iterations = control.last_step();
2988 }
2989 else if (data.eigenvalue_algorithm ==
2990 AdditionalData::EigenvalueAlgorithm::power_iteration)
2991 {
2992 Assert(data.degree != numbers::invalid_unsigned_int,
2993 ExcMessage("Cannot estimate the minimal eigenvalue with the "
2994 "power iteration"));
2995
2996 eigenvalue_tracker.values.push_back(
2997 internal::PreconditionChebyshevImplementation::power_iteration(
2998 *matrix_ptr,
2999 temp_vector1,
3000 *data.preconditioner,
3001 data.eig_cg_n_iterations));
3002 }
3003 else
3004 Assert(false, ExcNotImplemented());
3005
3006 // read the eigenvalues from the attached eigenvalue tracker
3007 if (eigenvalue_tracker.values.empty())
3008 info.min_eigenvalue_estimate = info.max_eigenvalue_estimate = 1.;
3009 else
3010 {
3011 info.min_eigenvalue_estimate = eigenvalue_tracker.values.front();
3012
3013 // include a safety factor since the CG method will in general not
3014 // be converged
3015 info.max_eigenvalue_estimate = 1.2 * eigenvalue_tracker.values.back();
3016 }
3017 }
3018 else
3019 {
3020 info.max_eigenvalue_estimate = data.max_eigenvalue;
3021 info.min_eigenvalue_estimate = data.max_eigenvalue / data.smoothing_range;
3022 }
3023
3024 const double alpha = (data.smoothing_range > 1. ?
3025 info.max_eigenvalue_estimate / data.smoothing_range :
3026 std::min(0.9 * info.max_eigenvalue_estimate,
3027 info.min_eigenvalue_estimate));
3028
3029 // in case the user set the degree to invalid unsigned int, we have to
3030 // determine the number of necessary iterations from the Chebyshev error
3031 // estimate, given the target tolerance specified by smoothing_range. This
3032 // estimate is based on the error formula given in section 5.1 of
3033 // R. S. Varga, Matrix iterative analysis, 2nd ed., Springer, 2009
3034 if (data.degree == numbers::invalid_unsigned_int)
3035 {
3036 const double actual_range = info.max_eigenvalue_estimate / alpha;
3037 const double sigma = (1. - std::sqrt(1. / actual_range)) /
3038 (1. + std::sqrt(1. / actual_range));
3039 const double eps = data.smoothing_range;
3040 const_cast<
3041 PreconditionChebyshev<MatrixType, VectorType, PreconditionerType> *>(
3042 this)
3043 ->data.degree =
3044 1 + static_cast<unsigned int>(
3045 std::log(1. / eps + std::sqrt(1. / eps / eps - 1.)) /
3046 std::log(1. / sigma));
3047 }
3048
3049 info.degree = data.degree;
3050
3051 const_cast<
3052 PreconditionChebyshev<MatrixType, VectorType, PreconditionerType> *>(this)
3053 ->delta = (info.max_eigenvalue_estimate - alpha) * 0.5;
3054 const_cast<
3055 PreconditionChebyshev<MatrixType, VectorType, PreconditionerType> *>(this)
3056 ->theta = (info.max_eigenvalue_estimate + alpha) * 0.5;
3057
3058 // We do not need the second temporary vector in case we have a
3059 // DiagonalMatrix as preconditioner and use deal.II's own vectors
3060 using NumberType = typename VectorType::value_type;
3061 if (std::is_same<PreconditionerType, DiagonalMatrix<VectorType>>::value ==
3062 false ||
3063 (std::is_same<VectorType, ::Vector<NumberType>>::value == false &&
3064 ((std::is_same<VectorType,
3065 LinearAlgebra::distributed::
3066 Vector<NumberType, MemorySpace::Host>>::value ==
3067 false) ||
3068 (std::is_same<VectorType,
3069 LinearAlgebra::distributed::
3070 Vector<NumberType, MemorySpace::CUDA>>::value ==
3071 false))))
3072 temp_vector2.reinit(src, true);
3073 else
3074 {
3075 VectorType empty_vector;
3076 temp_vector2.reinit(empty_vector);
3077 }
3078
3079 const_cast<
3080 PreconditionChebyshev<MatrixType, VectorType, PreconditionerType> *>(this)
3081 ->eigenvalues_are_initialized = true;
3082
3083 return info;
3084}
3085
3086
3087
3088template <typename MatrixType, typename VectorType, typename PreconditionerType>
3089inline void
3090PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::vmult(
3091 VectorType & solution,
3092 const VectorType &rhs) const
3093{
3094 std::lock_guard<std::mutex> lock(mutex);
3095 if (eigenvalues_are_initialized == false)
3096 estimate_eigenvalues(rhs);
3097
3098 internal::PreconditionChebyshevImplementation::vmult_and_update(
3099 *matrix_ptr,
3100 *data.preconditioner,
3101 rhs,
3102 0,
3103 0.,
3104 1. / theta,
3105 solution,
3106 solution_old,
3107 temp_vector1,
3108 temp_vector2);
3109
3110 // if delta is zero, we do not need to iterate because the updates will be
3111 // zero
3112 if (data.degree < 2 || std::abs(delta) < 1e-40)
3113 return;
3114
3115 double rhok = delta / theta, sigma = theta / delta;
3116 for (unsigned int k = 0; k < data.degree - 1; ++k)
3117 {
3118 const double rhokp = 1. / (2. * sigma - rhok);
3119 const double factor1 = rhokp * rhok, factor2 = 2. * rhokp / delta;
3120 rhok = rhokp;
3121 internal::PreconditionChebyshevImplementation::vmult_and_update(
3122 *matrix_ptr,
3123 *data.preconditioner,
3124 rhs,
3125 k + 1,
3126 factor1,
3127 factor2,
3128 solution,
3129 solution_old,
3130 temp_vector1,
3131 temp_vector2);
3132 }
3133}
3134
3135
3136
3137template <typename MatrixType, typename VectorType, typename PreconditionerType>
3138inline void
3139PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::Tvmult(
3140 VectorType & solution,
3141 const VectorType &rhs) const
3142{
3143 std::lock_guard<std::mutex> lock(mutex);
3144 if (eigenvalues_are_initialized == false)
3145 estimate_eigenvalues(rhs);
3146
3147 internal::PreconditionChebyshevImplementation::vector_updates(
3148 rhs,
3149 *data.preconditioner,
3150 0,
3151 0.,
3152 1. / theta,
3153 solution_old,
3154 temp_vector1,
3155 temp_vector2,
3156 solution);
3157
3158 if (data.degree < 2 || std::abs(delta) < 1e-40)
3159 return;
3160
3161 double rhok = delta / theta, sigma = theta / delta;
3162 for (unsigned int k = 0; k < data.degree - 1; ++k)
3163 {
3164 const double rhokp = 1. / (2. * sigma - rhok);
3165 const double factor1 = rhokp * rhok, factor2 = 2. * rhokp / delta;
3166 rhok = rhokp;
3167 matrix_ptr->Tvmult(temp_vector1, solution);
3168 internal::PreconditionChebyshevImplementation::vector_updates(
3169 rhs,
3170 *data.preconditioner,
3171 k + 1,
3172 factor1,
3173 factor2,
3174 solution_old,
3175 temp_vector1,
3176 temp_vector2,
3177 solution);
3178 }
3179}
3180
3181
3182
3183template <typename MatrixType, typename VectorType, typename PreconditionerType>
3184inline void
3185PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::step(
3186 VectorType & solution,
3187 const VectorType &rhs) const
3188{
3189 std::lock_guard<std::mutex> lock(mutex);
3190 if (eigenvalues_are_initialized == false)
3191 estimate_eigenvalues(rhs);
3192
3193 internal::PreconditionChebyshevImplementation::vmult_and_update(
3194 *matrix_ptr,
3195 *data.preconditioner,
3196 rhs,
3197 1,
3198 0.,
3199 1. / theta,
3200 solution,
3201 solution_old,
3202 temp_vector1,
3203 temp_vector2);
3204
3205 if (data.degree < 2 || std::abs(delta) < 1e-40)
3206 return;
3207
3208 double rhok = delta / theta, sigma = theta / delta;
3209 for (unsigned int k = 0; k < data.degree - 1; ++k)
3210 {
3211 const double rhokp = 1. / (2. * sigma - rhok);
3212 const double factor1 = rhokp * rhok, factor2 = 2. * rhokp / delta;
3213 rhok = rhokp;
3214 internal::PreconditionChebyshevImplementation::vmult_and_update(
3215 *matrix_ptr,
3216 *data.preconditioner,
3217 rhs,
3218 k + 2,
3219 factor1,
3220 factor2,
3221 solution,
3222 solution_old,
3223 temp_vector1,
3224 temp_vector2);
3225 }
3226}
3227
3228
3229
3230template <typename MatrixType, typename VectorType, typename PreconditionerType>
3231inline void
3232PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::Tstep(
3233 VectorType & solution,
3234 const VectorType &rhs) const
3235{
3236 std::lock_guard<std::mutex> lock(mutex);
3237 if (eigenvalues_are_initialized == false)
3238 estimate_eigenvalues(rhs);
3239
3240 matrix_ptr->Tvmult(temp_vector1, solution);
3241 internal::PreconditionChebyshevImplementation::vector_updates(
3242 rhs,
3243 *data.preconditioner,
3244 1,
3245 0.,
3246 1. / theta,
3247 solution_old,
3248 temp_vector1,
3249 temp_vector2,
3250 solution);
3251
3252 if (data.degree < 2 || std::abs(delta) < 1e-40)
3253 return;
3254
3255 double rhok = delta / theta, sigma = theta / delta;
3256 for (unsigned int k = 0; k < data.degree - 1; ++k)
3257 {
3258 const double rhokp = 1. / (2. * sigma - rhok);
3259 const double factor1 = rhokp * rhok, factor2 = 2. * rhokp / delta;
3260 rhok = rhokp;
3261 matrix_ptr->Tvmult(temp_vector1, solution);
3262 internal::PreconditionChebyshevImplementation::vector_updates(
3263 rhs,
3264 *data.preconditioner,
3265 k + 2,
3266 factor1,
3267 factor2,
3268 solution_old,
3269 temp_vector1,
3270 temp_vector2,
3271 solution);
3272 }
3273}
3274
3275
3276
3277template <typename MatrixType, typename VectorType, typename PreconditionerType>
3278inline typename PreconditionChebyshev<MatrixType,
3279 VectorType,
3280 PreconditionerType>::size_type
3281PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::m() const
3282{
3283 Assert(matrix_ptr != nullptr, ExcNotInitialized());
3284 return matrix_ptr->m();
3285}
3286
3287
3288
3289template <typename MatrixType, typename VectorType, typename PreconditionerType>
3290inline typename PreconditionChebyshev<MatrixType,
3291 VectorType,
3292 PreconditionerType>::size_type
3293PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::n() const
3294{
3295 Assert(matrix_ptr != nullptr, ExcNotInitialized());
3296 return matrix_ptr->n();
3297}
3298
3299#endif // DOXYGEN
3300
3301DEAL_II_NAMESPACE_CLOSE
3302
3303#endif
affine_constraints.h
block_vector_base.h
IndexSet::is_empty
bool is_empty() const
Definition: index_set.h:1826
IndexSet::nth_index_in_set
size_type nth_index_in_set(const size_type local_index) const
Definition: index_set.h:1882
SparseMatrix::end
const_iterator end() const
SparseMatrix::begin
const_iterator begin() const
Vector
Definition: vector.h:109
DEAL_II_OPENMP_SIMD_PRAGMA
#define DEAL_II_OPENMP_SIMD_PRAGMA
Definition: config.h:142
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
cuda_size.h
diagonal_matrix.h
AssertCudaKernel
#define AssertCudaKernel()
Definition: exceptions.h:1871
BlockVectorBase::n_blocks
unsigned int n_blocks() const
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1473
DataOutBase::eps
@ eps
Definition: data_out_base.h:1621
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
BlockVectorBase::block
BlockType & block(const unsigned int i)
StandardExceptions::ExcNotInitialized
static ::ExceptionBase & ExcNotInitialized()
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
AssertThrow
#define AssertThrow(cond, exc)
Definition: exceptions.h:1583
PreconditionRichardson::size_type
types::global_dof_index size_type
Definition: precondition.h:202
PreconditionRichardson::vmult_add
void vmult_add(VectorType &, const VectorType &) const
PreconditionSOR::PreconditionerType
internal::PreconditionRelaxation::PreconditionSORImpl< MatrixType > PreconditionerType
Definition: precondition.h:1208
PreconditionJacobi::AdditionalData
typename BaseClass::AdditionalData AdditionalData
Definition: precondition.h:1145
PreconditionRelaxation::A
SmartPointer< const MatrixType, PreconditionRelaxation< MatrixType > > A
Definition: precondition.h:503
PreconditionChebyshev::Tvmult
void Tvmult(VectorType &dst, const VectorType &src) const
PreconditionPSOR::AdditionalData::AdditionalData
AdditionalData(const std::vector< size_type > &permutation, const std::vector< size_type > &inverse_permutation, const typename BaseClass::AdditionalData &parameters=typename BaseClass::AdditionalData())
PreconditionRelaxation::clear
void clear()
PreconditionRichardson::n_columns
size_type n_columns
Definition: precondition.h:314
PreconditionPSOR::size_type
typename BaseClass::size_type size_type
Definition: precondition.h:1323
PreconditionRelaxation::AdditionalData::relaxation
double relaxation
Definition: precondition.h:427
PreconditionIdentity::PreconditionIdentity
PreconditionIdentity()
PreconditionRelaxation::Tvmult
void Tvmult(VectorType &, const VectorType &) const
PreconditionRelaxation::preconditioner
std::shared_ptr< PreconditionerType > preconditioner
Definition: precondition.h:518
PreconditionChebyshev::m
size_type m() const
PreconditionChebyshev::AdditionalData::eig_cg_residual
double eig_cg_residual
Definition: precondition.h:1687
PreconditionChebyshev::mutex
Threads::Mutex mutex
Definition: precondition.h:1880
PreconditionRichardson::n
size_type n() const
PreconditionRelaxation::AdditionalData::n_iterations
unsigned int n_iterations
Definition: precondition.h:432
PreconditionChebyshev::AdditionalData::preconditioner
std::shared_ptr< PreconditionerType > preconditioner
Definition: precondition.h:1705
PreconditionIdentity::n_columns
size_type n_columns
Definition: precondition.h:181
PreconditionUseMatrix::function_ptr
void(MatrixType::*)(VectorType &, const VectorType &) const function_ptr
Definition: precondition.h:367
PreconditionIdentity::vmult_add
void vmult_add(VectorType &, const VectorType &) const
PreconditionRichardson::relaxation
double relaxation
Definition: precondition.h:304
PreconditionRelaxation::n
size_type n() const
PreconditionRichardson::n_rows
size_type n_rows
Definition: precondition.h:309
PreconditionChebyshev::AdditionalData::smoothing_range
double smoothing_range
Definition: precondition.h:1673
PreconditionChebyshev::n
size_type n() const
PreconditionUseMatrix::precondition
const function_ptr precondition
Definition: precondition.h:392
PreconditionChebyshev::step
void step(VectorType &dst, const VectorType &src) const
PreconditionChebyshev::EigenvalueInformation::min_eigenvalue_estimate
double min_eigenvalue_estimate
Definition: precondition.h:1789
PreconditionChebyshev::data
AdditionalData data
Definition: precondition.h:1857
PreconditionChebyshev::solution_old
VectorType solution_old
Definition: precondition.h:1841
PreconditionChebyshev::EigenvalueInformation::max_eigenvalue_estimate
double max_eigenvalue_estimate
Definition: precondition.h:1793
PreconditionSSOR::initialize
void initialize(const MatrixType &A, const AdditionalData &parameters=AdditionalData())
PreconditionPSOR::AdditionalData::parameters
BaseClass::AdditionalData parameters
Definition: precondition.h:1357
PreconditionRichardson::initialize
void initialize(const AdditionalData &parameters)
PreconditionIdentity::vmult
void vmult(VectorType &, const VectorType &) const
PreconditionChebyshev::AdditionalData::operator=
AdditionalData & operator=(const AdditionalData &other_data)
PreconditionIdentity::m
size_type m() const
PreconditionRichardson::AdditionalData::relaxation
double relaxation
Definition: precondition.h:219
PreconditionPSOR::initialize
void initialize(const MatrixType &A, const AdditionalData &additional_data)
PreconditionChebyshev::estimate_eigenvalues
EigenvalueInformation estimate_eigenvalues(const VectorType &src) const
PreconditionChebyshev::AdditionalData::constraints
AffineConstraints< double > constraints
Definition: precondition.h:1700
PreconditionChebyshev::Tstep
void Tstep(VectorType &dst, const VectorType &src) const
PreconditionRelaxation::AdditionalData::preconditioner
std::shared_ptr< PreconditionerType > preconditioner
Definition: precondition.h:438
PreconditionPSOR::AdditionalData::inverse_permutation
const std::vector< size_type > & inverse_permutation
Definition: precondition.h:1353
PreconditionIdentity::AdditionalData::AdditionalData
AdditionalData()=default
PreconditionUseMatrix::matrix
const MatrixType & matrix
Definition: precondition.h:387
PreconditionSOR::AdditionalData
typename BaseClass::AdditionalData AdditionalData
Definition: precondition.h:1215
PreconditionIdentity::initialize
void initialize(const MatrixType &matrix, const AdditionalData &additional_data=AdditionalData())
PreconditionChebyshev::clear
void clear()
PreconditionRichardson::initialize
void initialize(const MatrixType &matrix, const AdditionalData &parameters)
PreconditionChebyshev::matrix_ptr
SmartPointer< const MatrixType, PreconditionChebyshev< MatrixType, VectorType, PreconditionerType > > matrix_ptr
Definition: precondition.h:1836
PreconditionChebyshev::AdditionalData::eig_cg_n_iterations
unsigned int eig_cg_n_iterations
Definition: precondition.h:1681
PreconditionSOR::initialize
void initialize(const MatrixType &A, const AdditionalData &parameters=AdditionalData())
PreconditionChebyshev::EigenvalueInformation::cg_iterations
unsigned int cg_iterations
Definition: precondition.h:1797
PreconditionChebyshev::delta
double delta
Definition: precondition.h:1868
PreconditionRelaxation::m
size_type m() const
PreconditionRichardson::m
size_type m() const
PreconditionIdentity::clear
void clear()
PreconditionIdentity::n_rows
size_type n_rows
Definition: precondition.h:176
PreconditionRichardson::AdditionalData::AdditionalData
AdditionalData(const double relaxation=1.)
PreconditionJacobi::PreconditionerType
internal::PreconditionRelaxation::PreconditionJacobiImpl< MatrixType > PreconditionerType
Definition: precondition.h:1138
PreconditionChebyshev::AdditionalData::degree
unsigned int degree
Definition: precondition.h:1660
PreconditionJacobi::initialize
void initialize(const MatrixType &A, const AdditionalData &parameters=AdditionalData())
PreconditionRelaxation::step
void step(VectorType &x, const VectorType &rhs) const
PreconditionRelaxation::AdditionalData::AdditionalData
AdditionalData(const double relaxation=1., const unsigned int n_iterations=1)
PreconditionRelaxation::n_iterations
unsigned int n_iterations
Definition: precondition.h:513
PreconditionIdentity::n
size_type n() const
PreconditionChebyshev::AdditionalData::EigenvalueAlgorithm
EigenvalueAlgorithm
Definition: precondition.h:1613
PreconditionPSOR::AdditionalData::permutation
const std::vector< size_type > & permutation
Definition: precondition.h:1349
PreconditionRichardson::vmult
void vmult(VectorType &, const VectorType &) const
PreconditionChebyshev::AdditionalData::AdditionalData
AdditionalData(const unsigned int degree=1, const double smoothing_range=0., const unsigned int eig_cg_n_iterations=8, const double eig_cg_residual=1e-2, const double max_eigenvalue=1, const EigenvalueAlgorithm eigenvalue_algorithm=EigenvalueAlgorithm::lanczos)
PreconditionRichardson::Tvmult
void Tvmult(VectorType &, const VectorType &) const
PreconditionRichardson::Tvmult_add
void Tvmult_add(VectorType &, const VectorType &) const
PreconditionChebyshev::temp_vector2
VectorType temp_vector2
Definition: precondition.h:1851
PreconditionChebyshev::AdditionalData::max_eigenvalue
double max_eigenvalue
Definition: precondition.h:1694
PreconditionRelaxation::initialize
void initialize(const MatrixType &A, const AdditionalData &parameters=AdditionalData())
PreconditionChebyshev::EigenvalueInformation::EigenvalueInformation
EigenvalueInformation()
Definition: precondition.h:1806
PreconditionChebyshev::vmult
void vmult(VectorType &dst, const VectorType &src) const
PreconditionChebyshev::temp_vector1
VectorType temp_vector1
Definition: precondition.h:1846
PreconditionPSOR::PreconditionerType
internal::PreconditionRelaxation::PreconditionPSORImpl< MatrixType > PreconditionerType
Definition: precondition.h:1316
PreconditionSSOR::AdditionalData
typename BaseClass::AdditionalData AdditionalData
Definition: precondition.h:1267
PreconditionSSOR::PreconditionerType
internal::PreconditionRelaxation::PreconditionSSORImpl< MatrixType > PreconditionerType
Definition: precondition.h:1260
PreconditionUseMatrix::vmult
void vmult(VectorType &dst, const VectorType &src) const
PreconditionRelaxation::Tstep
void Tstep(VectorType &x, const VectorType &rhs) const
PreconditionChebyshev::theta
double theta
Definition: precondition.h:1862
PreconditionRelaxation::vmult
void vmult(VectorType &, const VectorType &) const
PreconditionChebyshev::initialize
void initialize(const MatrixType &matrix, const AdditionalData &additional_data=AdditionalData())
PreconditionUseMatrix::PreconditionUseMatrix
PreconditionUseMatrix(const MatrixType &M, const function_ptr method)
PreconditionPSOR::initialize
void initialize(const MatrixType &A, const std::vector< size_type > &permutation, const std::vector< size_type > &inverse_permutation, const typename BaseClass::AdditionalData &parameters=typename BaseClass::AdditionalData())
PreconditionChebyshev::eigenvalues_are_initialized
bool eigenvalues_are_initialized
Definition: precondition.h:1874
PreconditionIdentity::Tvmult
void Tvmult(VectorType &, const VectorType &) const
PreconditionRichardson::clear
void clear()
Definition: precondition.h:277
PreconditionChebyshev::AdditionalData::eigenvalue_algorithm
EigenvalueAlgorithm eigenvalue_algorithm
Definition: precondition.h:1710
PreconditionIdentity::Tvmult_add
void Tvmult_add(VectorType &, const VectorType &) const
PreconditionChebyshev::PreconditionChebyshev
PreconditionChebyshev()
PreconditionRelaxation::relaxation
double relaxation
Definition: precondition.h:508
PreconditionRichardson::PreconditionRichardson
PreconditionRichardson()
PreconditionChebyshev::EigenvalueInformation::degree
unsigned int degree
Definition: precondition.h:1802
PreconditionRelaxation::size_type
types::global_dof_index size_type
Definition: precondition.h:410
PreconditionChebyshev::AdditionalData::EigenvalueAlgorithm::power_iteration
@ power_iteration
PreconditionChebyshev::AdditionalData::EigenvalueAlgorithm::lanczos
@ lanczos
Vector::add
void add(const std::vector< size_type > &indices, const std::vector< OtherNumber > &values)
LinearAlgebra::distributed::Vector::local_element
Number local_element(const size_type local_index) const
LinearAlgebra::distributed::Vector::get_values
Number * get_values() const
LinearAlgebra::distributed::Vector::swap
void swap(Vector< Number, MemorySpace > &v)
Vector::mean_value
Number mean_value() const
LinearAlgebra::distributed::Vector::locally_owned_size
size_type locally_owned_size() const
LinearAlgebra::distributed::Vector::mean_value
virtual Number mean_value() const override
Vector::size
size_type size() const
Vector::swap
virtual void swap(Vector< Number > &v)
LinearAlgebra::distributed::Vector::add
virtual void add(const Number a) override
Vector::begin
iterator begin()
LinearAlgebra::distributed::Vector::begin
iterator begin()
LinearAlgebra::distributed::Vector::locally_owned_elements
virtual ::IndexSet locally_owned_elements() const override
identity_matrix.h
vector_memory.h
memory_space.h
CUDAWrappers::block_size
constexpr int block_size
Definition: cuda_size.h:29
LAPACKSupport::matrix
@ matrix
Contents is actually a matrix.
Definition: lapack_support.h:58
MatrixFreeOperators::BlockHelper::n_blocks
std::enable_if< IsBlockVector< VectorType >::value, unsignedint >::type n_blocks(const VectorType &vector)
Definition: operators.h:50
TrilinosWrappers::internal::end
VectorType::value_type * end(VectorType &V)
Definition: trilinos_sparse_matrix.cc:67
internal::VectorImplementation::minimum_parallel_grain_size
unsigned int minimum_parallel_grain_size
Definition: parallel.cc:34
numbers::invalid_unsigned_int
static const unsigned int invalid_unsigned_int
Definition: types.h:201
std
STL namespace.
std::log
::VectorizedArray< Number, width > log(const ::VectorizedArray< Number, width > &)
Definition: vectorization.h:5488
std::min
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
Definition: vectorization.h:5605
std::sqrt
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
Definition: vectorization.h:5510
std::abs
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
Definition: vectorization.h:5572
types::global_dof_index
unsigned int global_dof_index
Definition: types.h:76
solver_cg.h
eigenvalues
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)
SymmetricTensorEigenvectorMethod::jacobi
@ jacobi
template_constraints.h
thread_management.h
Generated by doxygen 1.9.5