Reference documentation for deal.II version 9.5.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 - 2023 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/mutex.h>
24#include <deal.II/base/parallel.h>
25#include <deal.II/base/smartpointer.h>
26#include <deal.II/base/template_constraints.h>
27
28#include <deal.II/lac/affine_constraints.h>
29#include <deal.II/lac/block_vector_base.h>
30#include <deal.II/lac/diagonal_matrix.h>
31#include <deal.II/lac/identity_matrix.h>
32#include <deal.II/lac/solver_cg.h>
33#include <deal.II/lac/vector_memory.h>
34
35#include <limits>
36
37DEAL_II_NAMESPACE_OPEN
38
39// forward declarations
40#ifndef DOXYGEN
41template <typename number>
42class Vector;
43template <typename number>
44class SparseMatrix;
45namespace LinearAlgebra
46{
47 namespace distributed
48 {
49 template <typename, typename>
50 class Vector;
51 template <typename>
52 class BlockVector;
53 } // namespace distributed
54} // namespace LinearAlgebra
55#endif
56
57
83class PreconditionIdentity : public Subscriptor
84{
85public:
89 using size_type = types::global_dof_index;
90
95 struct AdditionalData
96 {
100 AdditionalData() = default;
101 };
102
106 PreconditionIdentity();
107
112 template <typename MatrixType>
113 void
114 initialize(const MatrixType & matrix,
115 const AdditionalData &additional_data = AdditionalData());
116
120 template <class VectorType>
121 void
122 vmult(VectorType &, const VectorType &) const;
123
128 template <class VectorType>
129 void
130 Tvmult(VectorType &, const VectorType &) const;
131
135 template <class VectorType>
136 void
137 vmult_add(VectorType &, const VectorType &) const;
138
143 template <class VectorType>
144 void
145 Tvmult_add(VectorType &, const VectorType &) const;
146
151 void
152 clear();
153
161 size_type
162 m() const;
163
171 size_type
172 n() const;
173
174private:
178 size_type n_rows;
179
183 size_type n_columns;
184};
185
186
187
198class PreconditionRichardson : public Subscriptor
199{
200public:
204 using size_type = types::global_dof_index;
205
209 class AdditionalData
210 {
211 public:
216 AdditionalData(const double relaxation = 1.);
217
221 double relaxation;
222 };
223
228 PreconditionRichardson();
229
233 void
234 initialize(const AdditionalData &parameters);
235
241 template <typename MatrixType>
242 void
243 initialize(const MatrixType &matrix, const AdditionalData &parameters);
244
248 template <class VectorType>
249 void
250 vmult(VectorType &, const VectorType &) const;
251
256 template <class VectorType>
257 void
258 Tvmult(VectorType &, const VectorType &) const;
262 template <class VectorType>
263 void
264 vmult_add(VectorType &, const VectorType &) const;
265
270 template <class VectorType>
271 void
272 Tvmult_add(VectorType &, const VectorType &) const;
273
278 void
279 clear()
280 {}
281
289 size_type
290 m() const;
291
299 size_type
300 n() const;
301
302private:
306 double relaxation;
307
311 size_type n_rows;
312
316 size_type n_columns;
317};
318
319
320
360template <typename MatrixType = SparseMatrix<double>,
361 class VectorType = Vector<double>>
362class PreconditionUseMatrix : public Subscriptor
363{
364public:
368 using function_ptr = void (MatrixType::*)(VectorType &,
369 const VectorType &) const;
370
376 PreconditionUseMatrix(const MatrixType &M, const function_ptr method);
377
382 void
383 vmult(VectorType &dst, const VectorType &src) const;
384
385private:
389 const MatrixType &matrix;
390
394 const function_ptr precondition;
395};
396
397
398
404template <typename MatrixType = SparseMatrix<double>,
405 typename PreconditionerType = IdentityMatrix>
406class PreconditionRelaxation : public Subscriptor
407{
408public:
412 using size_type = types::global_dof_index;
413
417 class AdditionalData
418 {
419 public:
423 AdditionalData(const double relaxation = 1.,
424 const unsigned int n_iterations = 1);
425
429 double relaxation;
430
434 unsigned int n_iterations;
435
436
437 /*
438 * Preconditioner.
439 */
440 std::shared_ptr<PreconditionerType> preconditioner;
441 };
442
448 void
449 initialize(const MatrixType & A,
450 const AdditionalData &parameters = AdditionalData());
451
455 void
456 clear();
457
462 size_type
463 m() const;
464
469 size_type
470 n() const;
471
475 template <class VectorType>
476 void
477 vmult(VectorType &, const VectorType &) const;
478
483 template <class VectorType>
484 void
485 Tvmult(VectorType &, const VectorType &) const;
486
490 template <class VectorType>
491 void
492 step(VectorType &x, const VectorType &rhs) const;
493
497 template <class VectorType>
498 void
499 Tstep(VectorType &x, const VectorType &rhs) const;
500
501protected:
505 SmartPointer<const MatrixType, PreconditionRelaxation<MatrixType>> A;
506
510 double relaxation;
511
515 unsigned int n_iterations;
516
517 /*
518 * Preconditioner.
519 */
520 std::shared_ptr<PreconditionerType> preconditioner;
521};
522
523
524
525#ifndef DOXYGEN
526
527namespace internal
528{
529 // a helper type-trait that leverage SFINAE to figure out if MatrixType has
530 // ... MatrixType::vmult(VectorType &, const VectorType&,
531 // std::function<...>, std::function<...>) const
532 template <typename MatrixType, typename VectorType>
533 using vmult_functions_t = decltype(std::declval<MatrixType const>().vmult(
534 std::declval<VectorType &>(),
535 std::declval<const VectorType &>(),
536 std::declval<
537 const std::function<void(const unsigned int, const unsigned int)> &>(),
538 std::declval<
539 const std::function<void(const unsigned int, const unsigned int)> &>()));
540
541 template <typename MatrixType,
542 typename VectorType,
543 typename PreconditionerType>
544 constexpr bool has_vmult_with_std_functions =
545 is_supported_operation<vmult_functions_t, MatrixType, VectorType> &&
546 std::is_same<PreconditionerType,
547 ::DiagonalMatrix<VectorType>>::value &&
548 (std::is_same<VectorType,
549 ::Vector<typename VectorType::value_type>>::value ||
550 std::is_same<
551 VectorType,
552 LinearAlgebra::distributed::Vector<typename VectorType::value_type,
553 MemorySpace::Host>>::value);
554
555
556 template <typename MatrixType, typename VectorType>
557 constexpr bool has_vmult_with_std_functions_for_precondition =
558 is_supported_operation<vmult_functions_t, MatrixType, VectorType>;
559
560 namespace PreconditionRelaxation
561 {
562 template <typename T, typename VectorType>
563 using Tvmult_t = decltype(
564 std::declval<T const>().Tvmult(std::declval<VectorType &>(),
565 std::declval<const VectorType &>()));
566
567 template <typename T, typename VectorType>
568 constexpr bool has_Tvmult = is_supported_operation<Tvmult_t, T, VectorType>;
569
570 template <typename T, typename VectorType>
571 using step_t = decltype(
572 std::declval<T const>().step(std::declval<VectorType &>(),
573 std::declval<const VectorType &>()));
574
575 template <typename T, typename VectorType>
576 constexpr bool has_step = is_supported_operation<step_t, T, VectorType>;
577
578 template <typename T, typename VectorType>
579 using step_omega_t =
580 decltype(std::declval<T const>().step(std::declval<VectorType &>(),
581 std::declval<const VectorType &>(),
582 std::declval<const double>()));
583
584 template <typename T, typename VectorType>
585 constexpr bool has_step_omega =
586 is_supported_operation<step_omega_t, T, VectorType>;
587
588 template <typename T, typename VectorType>
589 using Tstep_t = decltype(
590 std::declval<T const>().Tstep(std::declval<VectorType &>(),
591 std::declval<const VectorType &>()));
592
593 template <typename T, typename VectorType>
594 constexpr bool has_Tstep = is_supported_operation<Tstep_t, T, VectorType>;
595
596 template <typename T, typename VectorType>
597 using Tstep_omega_t =
598 decltype(std::declval<T const>().Tstep(std::declval<VectorType &>(),
599 std::declval<const VectorType &>(),
600 std::declval<const double>()));
601
602 template <typename T, typename VectorType>
603 constexpr bool has_Tstep_omega =
604 is_supported_operation<Tstep_omega_t, T, VectorType>;
605
606 template <typename T, typename VectorType>
607 using jacobi_step_t = decltype(
608 std::declval<T const>().Jacobi_step(std::declval<VectorType &>(),
609 std::declval<const VectorType &>(),
610 std::declval<const double>()));
611
612 template <typename T, typename VectorType>
613 constexpr bool has_jacobi_step =
614 is_supported_operation<jacobi_step_t, T, VectorType>;
615
616 template <typename T, typename VectorType>
617 using SOR_step_t = decltype(
618 std::declval<T const>().SOR_step(std::declval<VectorType &>(),
619 std::declval<const VectorType &>(),
620 std::declval<const double>()));
621
622 template <typename T, typename VectorType>
623 constexpr bool has_SOR_step =
624 is_supported_operation<SOR_step_t, T, VectorType>;
625
626 template <typename T, typename VectorType>
627 using SSOR_step_t = decltype(
628 std::declval<T const>().SSOR_step(std::declval<VectorType &>(),
629 std::declval<const VectorType &>(),
630 std::declval<const double>()));
631
632 template <typename T, typename VectorType>
633 constexpr bool has_SSOR_step =
634 is_supported_operation<SSOR_step_t, T, VectorType>;
635
636 template <typename MatrixType>
637 class PreconditionJacobiImpl
638 {
639 public:
640 PreconditionJacobiImpl(const MatrixType &A, const double relaxation)
641 : A(&A)
642 , relaxation(relaxation)
643 {}
644
645 template <typename VectorType>
646 void
647 vmult(VectorType &dst, const VectorType &src) const
648 {
649 this->A->precondition_Jacobi(dst, src, this->relaxation);
650 }
651
652 template <typename VectorType>
653 void
654 Tvmult(VectorType &dst, const VectorType &src) const
655 {
656 // call vmult, since preconditioner is symmetrical
657 this->vmult(dst, src);
658 }
659
660 template <typename VectorType,
661 std::enable_if_t<has_jacobi_step<MatrixType, VectorType>,
662 MatrixType> * = nullptr>
663 void
664 step(VectorType &dst, const VectorType &src) const
665 {
666 this->A->Jacobi_step(dst, src, this->relaxation);
667 }
668
669 template <typename VectorType,
670 std::enable_if_t<!has_jacobi_step<MatrixType, VectorType>,
671 MatrixType> * = nullptr>
672 void
673 step(VectorType &, const VectorType &) const
674 {
675 AssertThrow(false,
676 ExcMessage(
677 "Matrix A does not provide a Jacobi_step() function!"));
678 }
679
680 template <typename VectorType>
681 void
682 Tstep(VectorType &dst, const VectorType &src) const
683 {
684 // call step, since preconditioner is symmetrical
685 this->step(dst, src);
686 }
687
688 private:
689 const SmartPointer<const MatrixType> A;
690 const double relaxation;
691 };
692
693 template <typename MatrixType>
694 class PreconditionSORImpl
695 {
696 public:
697 PreconditionSORImpl(const MatrixType &A, const double relaxation)
698 : A(&A)
699 , relaxation(relaxation)
700 {}
701
702 template <typename VectorType>
703 void
704 vmult(VectorType &dst, const VectorType &src) const
705 {
706 this->A->precondition_SOR(dst, src, this->relaxation);
707 }
708
709 template <typename VectorType>
710 void
711 Tvmult(VectorType &dst, const VectorType &src) const
712 {
713 this->A->precondition_TSOR(dst, src, this->relaxation);
714 }
715
716 template <typename VectorType,
717 std::enable_if_t<has_SOR_step<MatrixType, VectorType>,
718 MatrixType> * = nullptr>
719 void
720 step(VectorType &dst, const VectorType &src) const
721 {
722 this->A->SOR_step(dst, src, this->relaxation);
723 }
724
725 template <typename VectorType,
726 std::enable_if_t<!has_SOR_step<MatrixType, VectorType>,
727 MatrixType> * = nullptr>
728 void
729 step(VectorType &, const VectorType &) const
730 {
731 AssertThrow(false,
732 ExcMessage(
733 "Matrix A does not provide a SOR_step() function!"));
734 }
735
736 template <typename VectorType,
737 std::enable_if_t<has_SOR_step<MatrixType, VectorType>,
738 MatrixType> * = nullptr>
739 void
740 Tstep(VectorType &dst, const VectorType &src) const
741 {
742 this->A->TSOR_step(dst, src, this->relaxation);
743 }
744
745 template <typename VectorType,
746 std::enable_if_t<!has_SOR_step<MatrixType, VectorType>,
747 MatrixType> * = nullptr>
748 void
749 Tstep(VectorType &, const VectorType &) const
750 {
751 AssertThrow(false,
752 ExcMessage(
753 "Matrix A does not provide a TSOR_step() function!"));
754 }
755
756 private:
757 const SmartPointer<const MatrixType> A;
758 const double relaxation;
759 };
760
761 template <typename MatrixType>
762 class PreconditionSSORImpl
763 {
764 public:
765 using size_type = typename MatrixType::size_type;
766
767 PreconditionSSORImpl(const MatrixType &A, const double relaxation)
768 : A(&A)
769 , relaxation(relaxation)
770 {
771 // in case we have a SparseMatrix class, we can extract information
772 // about the diagonal.
773 const SparseMatrix<typename MatrixType::value_type> *mat =
774 dynamic_cast<const SparseMatrix<typename MatrixType::value_type> *>(
775 &*this->A);
776
777 // calculate the positions first after the diagonal.
778 if (mat != nullptr)
779 {
780 const size_type n = this->A->n();
781 pos_right_of_diagonal.resize(n, static_cast<std::size_t>(-1));
782 for (size_type row = 0; row < n; ++row)
783 {
784 // find the first element in this line which is on the right of
785 // the diagonal. we need to precondition with the elements on
786 // the left only. note: the first entry in each line denotes the
787 // diagonal element, which we need not check.
788 typename SparseMatrix<
789 typename MatrixType::value_type>::const_iterator it =
790 mat->begin(row) + 1;
791 for (; it < mat->end(row); ++it)
792 if (it->column() > row)
793 break;
794 pos_right_of_diagonal[row] = it - mat->begin();
795 }
796 }
797 }
798
799 template <typename VectorType>
800 void
801 vmult(VectorType &dst, const VectorType &src) const
802 {
803 this->A->precondition_SSOR(dst,
804 src,
805 this->relaxation,
806 pos_right_of_diagonal);
807 }
808
809 template <typename VectorType>
810 void
811 Tvmult(VectorType &dst, const VectorType &src) const
812 {
813 this->A->precondition_SSOR(dst,
814 src,
815 this->relaxation,
816 pos_right_of_diagonal);
817 }
818
819 template <typename VectorType,
820 std::enable_if_t<has_SSOR_step<MatrixType, VectorType>,
821 MatrixType> * = nullptr>
822 void
823 step(VectorType &dst, const VectorType &src) const
824 {
825 this->A->SSOR_step(dst, src, this->relaxation);
826 }
827
828 template <typename VectorType,
829 std::enable_if_t<!has_SSOR_step<MatrixType, VectorType>,
830 MatrixType> * = nullptr>
831 void
832 step(VectorType &, const VectorType &) const
833 {
834 AssertThrow(false,
835 ExcMessage(
836 "Matrix A does not provide a SSOR_step() function!"));
837 }
838
839 template <typename VectorType>
840 void
841 Tstep(VectorType &dst, const VectorType &src) const
842 {
843 // call step, since preconditioner is symmetrical
844 this->step(dst, src);
845 }
846
847 private:
848 const SmartPointer<const MatrixType> A;
849 const double relaxation;
850
855 std::vector<std::size_t> pos_right_of_diagonal;
856 };
857
858 template <typename MatrixType>
859 class PreconditionPSORImpl
860 {
861 public:
862 using size_type = typename MatrixType::size_type;
863
864 PreconditionPSORImpl(const MatrixType & A,
865 const double relaxation,
866 const std::vector<size_type> &permutation,
867 const std::vector<size_type> &inverse_permutation)
868 : A(&A)
869 , relaxation(relaxation)
870 , permutation(permutation)
871 , inverse_permutation(inverse_permutation)
872 {}
873
874 template <typename VectorType>
875 void
876 vmult(VectorType &dst, const VectorType &src) const
877 {
878 dst = src;
879 this->A->PSOR(dst, permutation, inverse_permutation, this->relaxation);
880 }
881
882 template <typename VectorType>
883 void
884 Tvmult(VectorType &dst, const VectorType &src) const
885 {
886 dst = src;
887 this->A->TPSOR(dst, permutation, inverse_permutation, this->relaxation);
888 }
889
890 private:
891 const SmartPointer<const MatrixType> A;
892 const double relaxation;
893
894 const std::vector<size_type> &permutation;
895 const std::vector<size_type> &inverse_permutation;
896 };
897
898 template <typename MatrixType,
899 typename PreconditionerType,
900 typename VectorType,
901 std::enable_if_t<has_step_omega<PreconditionerType, VectorType>,
902 PreconditionerType> * = nullptr>
903 void
904 step(const MatrixType &,
905 const PreconditionerType &preconditioner,
906 VectorType & dst,
907 const VectorType & src,
908 const double relaxation,
909 VectorType &,
910 VectorType &)
911 {
912 preconditioner.step(dst, src, relaxation);
913 }
914
915 template <
916 typename MatrixType,
917 typename PreconditionerType,
918 typename VectorType,
919 std::enable_if_t<!has_step_omega<PreconditionerType, VectorType> &&
920 has_step<PreconditionerType, VectorType>,
921 PreconditionerType> * = nullptr>
922 void
923 step(const MatrixType &,
924 const PreconditionerType &preconditioner,
925 VectorType & dst,
926 const VectorType & src,
927 const double relaxation,
928 VectorType &,
929 VectorType &)
930 {
931 Assert(relaxation == 1.0, ExcInternalError());
932
933 (void)relaxation;
934
935 preconditioner.step(dst, src);
936 }
937
938 template <
939 typename MatrixType,
940 typename PreconditionerType,
941 typename VectorType,
942 std::enable_if_t<!has_step_omega<PreconditionerType, VectorType> &&
943 !has_step<PreconditionerType, VectorType>,
944 PreconditionerType> * = nullptr>
945 void
946 step(const MatrixType & A,
947 const PreconditionerType &preconditioner,
948 VectorType & dst,
949 const VectorType & src,
950 const double relaxation,
951 VectorType & residual,
952 VectorType & tmp)
953 {
954 residual.reinit(dst, true);
955 tmp.reinit(dst, true);
956
957 A.vmult(residual, dst);
958 residual.sadd(-1.0, 1.0, src);
959
960 preconditioner.vmult(tmp, residual);
961 dst.add(relaxation, tmp);
962 }
963
964 template <typename MatrixType,
965 typename PreconditionerType,
966 typename VectorType,
967 std::enable_if_t<has_Tstep_omega<PreconditionerType, VectorType>,
968 PreconditionerType> * = nullptr>
969 void
970 Tstep(const MatrixType &,
971 const PreconditionerType &preconditioner,
972 VectorType & dst,
973 const VectorType & src,
974 const double relaxation,
975 VectorType &,
976 VectorType &)
977 {
978 preconditioner.Tstep(dst, src, relaxation);
979 }
980
981 template <
982 typename MatrixType,
983 typename PreconditionerType,
984 typename VectorType,
985 std::enable_if_t<!has_Tstep_omega<PreconditionerType, VectorType> &&
986 has_Tstep<PreconditionerType, VectorType>,
987 PreconditionerType> * = nullptr>
988 void
989 Tstep(const MatrixType &,
990 const PreconditionerType &preconditioner,
991 VectorType & dst,
992 const VectorType & src,
993 const double relaxation,
994 VectorType &,
995 VectorType &)
996 {
997 Assert(relaxation == 1.0, ExcInternalError());
998
999 (void)relaxation;
1000
1001 preconditioner.Tstep(dst, src);
1002 }
1003
1004 template <typename MatrixType,
1005 typename VectorType,
1006 std::enable_if_t<has_Tvmult<MatrixType, VectorType>, MatrixType>
1007 * = nullptr>
1008 void
1009 Tvmult(const MatrixType &A, VectorType &dst, const VectorType &src)
1010 {
1011 A.Tvmult(dst, src);
1012 }
1013
1014 template <typename MatrixType,
1015 typename VectorType,
1016 std::enable_if_t<!has_Tvmult<MatrixType, VectorType>, MatrixType>
1017 * = nullptr>
1018 void
1019 Tvmult(const MatrixType &, VectorType &, const VectorType &)
1020 {
1021 AssertThrow(false,
1022 ExcMessage("Matrix A does not provide a Tvmult() function!"));
1023 }
1024
1025 template <
1026 typename MatrixType,
1027 typename PreconditionerType,
1028 typename VectorType,
1029 std::enable_if_t<!has_Tstep_omega<PreconditionerType, VectorType> &&
1030 !has_Tstep<PreconditionerType, VectorType>,
1031 PreconditionerType> * = nullptr>
1032 void
1033 Tstep(const MatrixType & A,
1034 const PreconditionerType &preconditioner,
1035 VectorType & dst,
1036 const VectorType & src,
1037 const double relaxation,
1038 VectorType & residual,
1039 VectorType & tmp)
1040 {
1041 residual.reinit(dst, true);
1042 tmp.reinit(dst, true);
1043
1044 Tvmult(A, residual, dst);
1045 residual.sadd(-1.0, 1.0, src);
1046
1047 Tvmult(preconditioner, tmp, residual);
1048 dst.add(relaxation, tmp);
1049 }
1050
1051 // 0) general implementation
1052 template <typename MatrixType,
1053 typename PreconditionerType,
1054 typename VectorType,
1055 std::enable_if_t<!has_vmult_with_std_functions_for_precondition<
1056 PreconditionerType,
1057 VectorType>,
1058 int> * = nullptr>
1059 void
1060 step_operations(const MatrixType & A,
1061 const PreconditionerType &preconditioner,
1062 VectorType & dst,
1063 const VectorType & src,
1064 const double relaxation,
1065 VectorType & tmp1,
1066 VectorType & tmp2,
1067 const unsigned int i,
1068 const bool transposed)
1069 {
1070 if (i == 0)
1071 {
1072 if (transposed)
1073 Tvmult(preconditioner, dst, src);
1074 else
1075 preconditioner.vmult(dst, src);
1076
1077 if (relaxation != 1.0)
1078 dst *= relaxation;
1079 }
1080 else
1081 {
1082 if (transposed)
1083 Tstep(A, preconditioner, dst, src, relaxation, tmp1, tmp2);
1084 else
1085 step(A, preconditioner, dst, src, relaxation, tmp1, tmp2);
1086 }
1087 }
1088
1089 // 1) specialized implementation with a preconditioner that accepts
1090 // ranges
1091 template <
1092 typename MatrixType,
1093 typename PreconditionerType,
1094 typename VectorType,
1095 std::enable_if_t<
1096 has_vmult_with_std_functions_for_precondition<PreconditionerType,
1097 VectorType> &&
1098 !has_vmult_with_std_functions_for_precondition<MatrixType,
1099 VectorType>,
1100 int> * = nullptr>
1101 void
1102 step_operations(const MatrixType & A,
1103 const PreconditionerType &preconditioner,
1104 VectorType & dst,
1105 const VectorType & src,
1106 const double relaxation,
1107 VectorType & tmp,
1108 VectorType &,
1109 const unsigned int i,
1110 const bool transposed)
1111 {
1112 (void)transposed;
1113 using Number = typename VectorType::value_type;
1114
1115 if (i == 0)
1116 {
1117 Number * dst_ptr = dst.begin();
1118 const Number *src_ptr = src.begin();
1119
1120 preconditioner.vmult(
1121 dst,
1122 src,
1123 [&](const unsigned int start_range, const unsigned int end_range) {
1124 // zero 'dst' before running the vmult operation
1125 if (end_range > start_range)
1126 std::memset(dst.begin() + start_range,
1127 0,
1128 sizeof(Number) * (end_range - start_range));
1129 },
1130 [&](const unsigned int start_range, const unsigned int end_range) {
1131 if (relaxation == 1.0)
1132 return; // nothing to do
1133
1134 const auto src_ptr = src.begin();
1135 const auto dst_ptr = dst.begin();
1136
1137 DEAL_II_OPENMP_SIMD_PRAGMA
1138 for (std::size_t i = start_range; i < end_range; ++i)
1139 dst_ptr[i] *= relaxation;
1140 });
1141 }
1142 else
1143 {
1144 tmp.reinit(src, true);
1145
1146 Assert(transposed == false, ExcNotImplemented());
1147
1148 A.vmult(tmp, dst);
1149
1150 preconditioner.vmult(
1151 dst,
1152 tmp,
1153 [&](const unsigned int start_range, const unsigned int end_range) {
1154 const auto src_ptr = src.begin();
1155 const auto tmp_ptr = tmp.begin();
1156
1157 if (relaxation == 1.0)
1158 {
1159 DEAL_II_OPENMP_SIMD_PRAGMA
1160 for (std::size_t i = start_range; i < end_range; ++i)
1161 tmp_ptr[i] = src_ptr[i] - tmp_ptr[i];
1162 }
1163 else
1164 {
1165 // note: we scale the residual here to be able to add into
1166 // the dst vector, which contains the solution from the last
1167 // iteration
1168 DEAL_II_OPENMP_SIMD_PRAGMA
1169 for (std::size_t i = start_range; i < end_range; ++i)
1170 tmp_ptr[i] = relaxation * (src_ptr[i] - tmp_ptr[i]);
1171 }
1172 },
1173 [&](const unsigned int, const unsigned int) {
1174 // nothing to do, since scaling by the relaxation factor
1175 // has been done in the pre operation
1176 });
1177 }
1178 }
1179
1180 // 2) specialized implementation with a preconditioner and a matrix that
1181 // accepts ranges
1182 template <
1183 typename MatrixType,
1184 typename PreconditionerType,
1185 typename VectorType,
1186 std::enable_if_t<
1187 has_vmult_with_std_functions_for_precondition<PreconditionerType,
1188 VectorType> &&
1189 has_vmult_with_std_functions_for_precondition<MatrixType, VectorType>,
1190 int> * = nullptr>
1191 void
1192 step_operations(const MatrixType & A,
1193 const PreconditionerType &preconditioner,
1194 VectorType & dst,
1195 const VectorType & src,
1196 const double relaxation,
1197 VectorType & tmp,
1198 VectorType &,
1199 const unsigned int i,
1200 const bool transposed)
1201 {
1202 (void)transposed;
1203 using Number = typename VectorType::value_type;
1204
1205 if (i == 0)
1206 {
1207 Number * dst_ptr = dst.begin();
1208 const Number *src_ptr = src.begin();
1209
1210 preconditioner.vmult(
1211 dst,
1212 src,
1213 [&](const unsigned int start_range, const unsigned int end_range) {
1214 // zero 'dst' before running the vmult operation
1215 if (end_range > start_range)
1216 std::memset(dst.begin() + start_range,
1217 0,
1218 sizeof(Number) * (end_range - start_range));
1219 },
1220 [&](const unsigned int start_range, const unsigned int end_range) {
1221 if (relaxation == 1.0)
1222 return; // nothing to do
1223
1224 const auto src_ptr = src.begin();
1225 const auto dst_ptr = dst.begin();
1226
1227 DEAL_II_OPENMP_SIMD_PRAGMA
1228 for (std::size_t i = start_range; i < end_range; ++i)
1229 dst_ptr[i] *= relaxation;
1230 });
1231 }
1232 else
1233 {
1234 tmp.reinit(src, true);
1235
1236 Assert(transposed == false, ExcNotImplemented());
1237
1238 A.vmult(
1239 tmp,
1240 dst,
1241 [&](const unsigned int start_range, const unsigned int end_range) {
1242 // zero 'tmp' before running the vmult
1243 // operation
1244 if (end_range > start_range)
1245 std::memset(tmp.begin() + start_range,
1246 0,
1247 sizeof(Number) * (end_range - start_range));
1248 },
1249 [&](const unsigned int start_range, const unsigned int end_range) {
1250 const auto src_ptr = src.begin();
1251 const auto tmp_ptr = tmp.begin();
1252
1253 if (relaxation == 1.0)
1254 {
1255 DEAL_II_OPENMP_SIMD_PRAGMA
1256 for (std::size_t i = start_range; i < end_range; ++i)
1257 tmp_ptr[i] = src_ptr[i] - tmp_ptr[i];
1258 }
1259 else
1260 {
1261 // note: we scale the residual here to be able to add into
1262 // the dst vector, which contains the solution from the last
1263 // iteration
1264 DEAL_II_OPENMP_SIMD_PRAGMA
1265 for (std::size_t i = start_range; i < end_range; ++i)
1266 tmp_ptr[i] = relaxation * (src_ptr[i] - tmp_ptr[i]);
1267 }
1268 });
1269
1270 preconditioner.vmult(dst, tmp, [](const auto, const auto) {
1271 // note: `dst` vector does not have to be zeroed out
1272 // since we add the result into it
1273 });
1274 }
1275 }
1276
1277 // 3) specialized implementation for inverse-diagonal preconditioner
1278 template <typename MatrixType,
1279 typename VectorType,
1280 std::enable_if_t<!IsBlockVector<VectorType>::value &&
1281 !has_vmult_with_std_functions<
1282 MatrixType,
1283 VectorType,
1284 ::DiagonalMatrix<VectorType>>,
1285 VectorType> * = nullptr>
1286 void
1287 step_operations(const MatrixType & A,
1288 const ::DiagonalMatrix<VectorType> &preconditioner,
1289 VectorType & dst,
1290 const VectorType & src,
1291 const double relaxation,
1292 VectorType & tmp,
1293 VectorType &,
1294 const unsigned int i,
1295 const bool transposed)
1296 {
1297 using Number = typename VectorType::value_type;
1298
1299 if (i == 0)
1300 {
1301 Number * dst_ptr = dst.begin();
1302 const Number *src_ptr = src.begin();
1303 const Number *diag_ptr = preconditioner.get_vector().begin();
1304
1305 if (relaxation == 1.0)
1306 {
1307 DEAL_II_OPENMP_SIMD_PRAGMA
1308 for (unsigned int i = 0; i < dst.locally_owned_size(); ++i)
1309 dst_ptr[i] = src_ptr[i] * diag_ptr[i];
1310 }
1311 else
1312 {
1313 DEAL_II_OPENMP_SIMD_PRAGMA
1314 for (unsigned int i = 0; i < dst.locally_owned_size(); ++i)
1315 dst_ptr[i] = relaxation * src_ptr[i] * diag_ptr[i];
1316 }
1317 }
1318 else
1319 {
1320 tmp.reinit(src, true);
1321
1322 Number * dst_ptr = dst.begin();
1323 const Number *src_ptr = src.begin();
1324 const Number *tmp_ptr = tmp.begin();
1325 const Number *diag_ptr = preconditioner.get_vector().begin();
1326
1327 if (transposed)
1328 Tvmult(A, tmp, dst);
1329 else
1330 A.vmult(tmp, dst);
1331
1332 if (relaxation == 1.0)
1333 {
1334 DEAL_II_OPENMP_SIMD_PRAGMA
1335 for (unsigned int i = 0; i < dst.locally_owned_size(); ++i)
1336 dst_ptr[i] += (src_ptr[i] - tmp_ptr[i]) * diag_ptr[i];
1337 }
1338 else
1339 {
1340 DEAL_II_OPENMP_SIMD_PRAGMA
1341 for (unsigned int i = 0; i < dst.locally_owned_size(); ++i)
1342 dst_ptr[i] +=
1343 relaxation * (src_ptr[i] - tmp_ptr[i]) * diag_ptr[i];
1344 }
1345 }
1346 }
1347
1348 // 4) specialized implementation for inverse-diagonal preconditioner and
1349 // matrix that accepts ranges
1350 template <typename MatrixType,
1351 typename VectorType,
1352 std::enable_if_t<!IsBlockVector<VectorType>::value &&
1353 has_vmult_with_std_functions<
1354 MatrixType,
1355 VectorType,
1356 ::DiagonalMatrix<VectorType>>,
1357 VectorType> * = nullptr>
1358 void
1359 step_operations(const MatrixType & A,
1360 const ::DiagonalMatrix<VectorType> &preconditioner,
1361 VectorType & dst,
1362 const VectorType & src,
1363 const double relaxation,
1364 VectorType & tmp,
1365 VectorType &,
1366 const unsigned int i,
1367 const bool transposed)
1368 {
1369 (void)transposed;
1370 using Number = typename VectorType::value_type;
1371
1372 if (i == 0)
1373 {
1374 Number * dst_ptr = dst.begin();
1375 const Number *src_ptr = src.begin();
1376 const Number *diag_ptr = preconditioner.get_vector().begin();
1377
1378 if (relaxation == 1.0)
1379 {
1380 DEAL_II_OPENMP_SIMD_PRAGMA
1381 for (unsigned int i = 0; i < dst.locally_owned_size(); ++i)
1382 dst_ptr[i] = src_ptr[i] * diag_ptr[i];
1383 }
1384 else
1385 {
1386 DEAL_II_OPENMP_SIMD_PRAGMA
1387 for (unsigned int i = 0; i < dst.locally_owned_size(); ++i)
1388 dst_ptr[i] = relaxation * src_ptr[i] * diag_ptr[i];
1389 }
1390 }
1391 else
1392 {
1393 tmp.reinit(src, true);
1394
1395 Assert(transposed == false, ExcNotImplemented());
1396
1397 A.vmult(
1398 tmp,
1399 dst,
1400 [&](const unsigned int start_range, const unsigned int end_range) {
1401 // zero 'tmp' before running the vmult operation
1402 if (end_range > start_range)
1403 std::memset(tmp.begin() + start_range,
1404 0,
1405 sizeof(Number) * (end_range - start_range));
1406 },
1407 [&](const unsigned int begin, const unsigned int end) {
1408 const Number *dst_ptr = dst.begin();
1409 const Number *src_ptr = src.begin();
1410 Number * tmp_ptr = tmp.begin();
1411 const Number *diag_ptr = preconditioner.get_vector().begin();
1412
1413 // for efficiency reason, write back to temp_vector that is
1414 // already read (avoid read-for-ownership)
1415 if (relaxation == 1.0)
1416 {
1417 DEAL_II_OPENMP_SIMD_PRAGMA
1418 for (std::size_t i = begin; i < end; ++i)
1419 tmp_ptr[i] =
1420 dst_ptr[i] + (src_ptr[i] - tmp_ptr[i]) * diag_ptr[i];
1421 }
1422 else
1423 {
1424 DEAL_II_OPENMP_SIMD_PRAGMA
1425 for (std::size_t i = begin; i < end; ++i)
1426 tmp_ptr[i] = dst_ptr[i] + relaxation *
1427 (src_ptr[i] - tmp_ptr[i]) *
1428 diag_ptr[i];
1429 }
1430 });
1431
1432 tmp.swap(dst);
1433 }
1434 }
1435
1436 } // namespace PreconditionRelaxation
1437} // namespace internal
1438
1439#endif
1440
1441
1442
1469template <typename MatrixType = SparseMatrix<double>>
1470class PreconditionJacobi
1471 : public PreconditionRelaxation<
1472 MatrixType,
1473 internal::PreconditionRelaxation::PreconditionJacobiImpl<MatrixType>>
1474{
1475 using PreconditionerType =
1476 internal::PreconditionRelaxation::PreconditionJacobiImpl<MatrixType>;
1477 using BaseClass = PreconditionRelaxation<MatrixType, PreconditionerType>;
1478
1479public:
1483 using AdditionalData = typename BaseClass::AdditionalData;
1484
1488 void
1489 initialize(const MatrixType & A,
1490 const AdditionalData &parameters = AdditionalData());
1491};
1492
1493
1539template <typename MatrixType = SparseMatrix<double>>
1540class PreconditionSOR
1541 : public PreconditionRelaxation<
1542 MatrixType,
1543 internal::PreconditionRelaxation::PreconditionSORImpl<MatrixType>>
1544{
1545 using PreconditionerType =
1546 internal::PreconditionRelaxation::PreconditionSORImpl<MatrixType>;
1547 using BaseClass = PreconditionRelaxation<MatrixType, PreconditionerType>;
1548
1549public:
1553 using AdditionalData = typename BaseClass::AdditionalData;
1554
1558 void
1559 initialize(const MatrixType & A,
1560 const AdditionalData &parameters = AdditionalData());
1561};
1562
1563
1564
1591template <typename MatrixType = SparseMatrix<double>>
1592class PreconditionSSOR
1593 : public PreconditionRelaxation<
1594 MatrixType,
1595 internal::PreconditionRelaxation::PreconditionSSORImpl<MatrixType>>
1596{
1597 using PreconditionerType =
1598 internal::PreconditionRelaxation::PreconditionSSORImpl<MatrixType>;
1599 using BaseClass = PreconditionRelaxation<MatrixType, PreconditionerType>;
1600
1601public:
1605 using AdditionalData = typename BaseClass::AdditionalData;
1606
1612 void
1613 initialize(const MatrixType & A,
1614 const AdditionalData &parameters = AdditionalData());
1615};
1616
1617
1647template <typename MatrixType = SparseMatrix<double>>
1648class PreconditionPSOR
1649 : public PreconditionRelaxation<
1650 MatrixType,
1651 internal::PreconditionRelaxation::PreconditionPSORImpl<MatrixType>>
1652{
1653 using PreconditionerType =
1654 internal::PreconditionRelaxation::PreconditionPSORImpl<MatrixType>;
1655 using BaseClass = PreconditionRelaxation<MatrixType, PreconditionerType>;
1656
1657public:
1661 using size_type = typename BaseClass::size_type;
1662
1666 class AdditionalData
1667 {
1668 public:
1679 AdditionalData(const std::vector<size_type> &permutation,
1680 const std::vector<size_type> &inverse_permutation,
1681 const typename BaseClass::AdditionalData &parameters =
1682 typename BaseClass::AdditionalData());
1683
1687 const std::vector<size_type> &permutation;
1691 const std::vector<size_type> &inverse_permutation;
1695 typename BaseClass::AdditionalData parameters;
1696 };
1697
1709 void
1710 initialize(const MatrixType & A,
1711 const std::vector<size_type> & permutation,
1712 const std::vector<size_type> & inverse_permutation,
1713 const typename BaseClass::AdditionalData &parameters =
1714 typename BaseClass::AdditionalData());
1715
1726 void
1727 initialize(const MatrixType &A, const AdditionalData &additional_data);
1728};
1729
1730
1731
1929template <typename MatrixType = SparseMatrix<double>,
1930 typename VectorType = Vector<double>,
1931 typename PreconditionerType = DiagonalMatrix<VectorType>>
1932class PreconditionChebyshev : public Subscriptor
1933{
1934public:
1938 using size_type = types::global_dof_index;
1939
1944 struct AdditionalData
1945 {
1950 enum class EigenvalueAlgorithm
1951 {
1957 lanczos,
1966 power_iteration
1967 };
1968
1972 enum class PolynomialType
1973 {
1977 first_kind,
1982 fourth_kind
1983 };
1984
1988 AdditionalData(
1989 const unsigned int degree = 1,
1990 const double smoothing_range = 0.,
1991 const unsigned int eig_cg_n_iterations = 8,
1992 const double eig_cg_residual = 1e-2,
1993 const double max_eigenvalue = 1,
1994 const EigenvalueAlgorithm eigenvalue_algorithm =
1995 EigenvalueAlgorithm::lanczos,
1996 const PolynomialType polynomial_type = PolynomialType::first_kind);
1997
2001 AdditionalData &
2002 operator=(const AdditionalData &other_data);
2003
2016 unsigned int degree;
2017
2029 double smoothing_range;
2030
2037 unsigned int eig_cg_n_iterations;
2038
2043 double eig_cg_residual;
2044
2050 double max_eigenvalue;
2051
2056 AffineConstraints<double> constraints;
2057
2061 std::shared_ptr<PreconditionerType> preconditioner;
2062
2066 EigenvalueAlgorithm eigenvalue_algorithm;
2067
2071 PolynomialType polynomial_type;
2072 };
2073
2074
2078 PreconditionChebyshev();
2079
2091 void
2092 initialize(const MatrixType & matrix,
2093 const AdditionalData &additional_data = AdditionalData());
2094
2099 void
2100 vmult(VectorType &dst, const VectorType &src) const;
2101
2106 void
2107 Tvmult(VectorType &dst, const VectorType &src) const;
2108
2112 void
2113 step(VectorType &dst, const VectorType &src) const;
2114
2118 void
2119 Tstep(VectorType &dst, const VectorType &src) const;
2120
2124 void
2125 clear();
2126
2131 size_type
2132 m() const;
2133
2138 size_type
2139 n() const;
2140
2174
2187 EigenvalueInformation
2188 estimate_eigenvalues(const VectorType &src) const;
2189
2190private:
2194 SmartPointer<
2195 const MatrixType,
2196 PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>>
2197 matrix_ptr;
2198
2202 mutable VectorType solution_old;
2203
2207 mutable VectorType temp_vector1;
2208
2212 mutable VectorType temp_vector2;
2213
2218 AdditionalData data;
2219
2223 double theta;
2224
2229 double delta;
2230
2235 bool eigenvalues_are_initialized;
2236
2241 mutable Threads::Mutex mutex;
2242};
2243
2244
2245
2247/* ---------------------------------- Inline functions ------------------- */
2248
2249#ifndef DOXYGEN
2250
2251inline PreconditionIdentity::PreconditionIdentity()
2252 : n_rows(0)
2253 , n_columns(0)
2254{}
2255
2256template <typename MatrixType>
2257inline void
2258PreconditionIdentity::initialize(const MatrixType &matrix,
2259 const PreconditionIdentity::AdditionalData &)
2260{
2261 n_rows = matrix.m();
2262 n_columns = matrix.n();
2263}
2264
2265
2266template <class VectorType>
2267inline void
2268PreconditionIdentity::vmult(VectorType &dst, const VectorType &src) const
2269{
2270 dst = src;
2271}
2272
2273
2274
2275template <class VectorType>
2276inline void
2277PreconditionIdentity::Tvmult(VectorType &dst, const VectorType &src) const
2278{
2279 dst = src;
2280}
2281
2282template <class VectorType>
2283inline void
2284PreconditionIdentity::vmult_add(VectorType &dst, const VectorType &src) const
2285{
2286 dst += src;
2287}
2288
2289
2290
2291template <class VectorType>
2292inline void
2293PreconditionIdentity::Tvmult_add(VectorType &dst, const VectorType &src) const
2294{
2295 dst += src;
2296}
2297
2298
2299
2300inline void
2301PreconditionIdentity::clear()
2302{}
2303
2304
2305
2306inline PreconditionIdentity::size_type
2307PreconditionIdentity::m() const
2308{
2309 Assert(n_rows != 0, ExcNotInitialized());
2310 return n_rows;
2311}
2312
2313inline PreconditionIdentity::size_type
2314PreconditionIdentity::n() const
2315{
2316 Assert(n_columns != 0, ExcNotInitialized());
2317 return n_columns;
2318}
2319
2320//---------------------------------------------------------------------------
2321
2322inline PreconditionRichardson::AdditionalData::AdditionalData(
2323 const double relaxation)
2324 : relaxation(relaxation)
2325{}
2326
2327
2328inline PreconditionRichardson::PreconditionRichardson()
2329 : relaxation(0)
2330 , n_rows(0)
2331 , n_columns(0)
2332{
2333 AdditionalData add_data;
2334 relaxation = add_data.relaxation;
2335}
2336
2337
2338
2339inline void
2340PreconditionRichardson::initialize(
2341 const PreconditionRichardson::AdditionalData &parameters)
2342{
2343 relaxation = parameters.relaxation;
2344}
2345
2346
2347
2348template <typename MatrixType>
2349inline void
2350PreconditionRichardson::initialize(
2351 const MatrixType & matrix,
2352 const PreconditionRichardson::AdditionalData &parameters)
2353{
2354 relaxation = parameters.relaxation;
2355 n_rows = matrix.m();
2356 n_columns = matrix.n();
2357}
2358
2359
2360
2361template <class VectorType>
2362inline void
2363PreconditionRichardson::vmult(VectorType &dst, const VectorType &src) const
2364{
2365 static_assert(
2366 std::is_same<size_type, typename VectorType::size_type>::value,
2367 "PreconditionRichardson and VectorType must have the same size_type.");
2368
2369 dst.equ(relaxation, src);
2370}
2371
2372
2373
2374template <class VectorType>
2375inline void
2376PreconditionRichardson::Tvmult(VectorType &dst, const VectorType &src) const
2377{
2378 static_assert(
2379 std::is_same<size_type, typename VectorType::size_type>::value,
2380 "PreconditionRichardson and VectorType must have the same size_type.");
2381
2382 dst.equ(relaxation, src);
2383}
2384
2385template <class VectorType>
2386inline void
2387PreconditionRichardson::vmult_add(VectorType &dst, const VectorType &src) const
2388{
2389 static_assert(
2390 std::is_same<size_type, typename VectorType::size_type>::value,
2391 "PreconditionRichardson and VectorType must have the same size_type.");
2392
2393 dst.add(relaxation, src);
2394}
2395
2396
2397
2398template <class VectorType>
2399inline void
2400PreconditionRichardson::Tvmult_add(VectorType &dst, const VectorType &src) const
2401{
2402 static_assert(
2403 std::is_same<size_type, typename VectorType::size_type>::value,
2404 "PreconditionRichardson and VectorType must have the same size_type.");
2405
2406 dst.add(relaxation, src);
2407}
2408
2409inline PreconditionRichardson::size_type
2410PreconditionRichardson::m() const
2411{
2412 Assert(n_rows != 0, ExcNotInitialized());
2413 return n_rows;
2414}
2415
2416inline PreconditionRichardson::size_type
2417PreconditionRichardson::n() const
2418{
2419 Assert(n_columns != 0, ExcNotInitialized());
2420 return n_columns;
2421}
2422
2423//---------------------------------------------------------------------------
2424
2425template <typename MatrixType, typename PreconditionerType>
2426inline void
2427PreconditionRelaxation<MatrixType, PreconditionerType>::initialize(
2428 const MatrixType & rA,
2429 const AdditionalData &parameters)
2430{
2431 A = &rA;
2432 relaxation = parameters.relaxation;
2433
2434 Assert(parameters.preconditioner, ExcNotInitialized());
2435
2436 preconditioner = parameters.preconditioner;
2437 n_iterations = parameters.n_iterations;
2438}
2439
2440
2441template <typename MatrixType, typename PreconditionerType>
2442inline void
2443PreconditionRelaxation<MatrixType, PreconditionerType>::clear()
2444{
2445 A = nullptr;
2446 preconditioner = nullptr;
2447}
2448
2449template <typename MatrixType, typename PreconditionerType>
2450inline
2451 typename PreconditionRelaxation<MatrixType, PreconditionerType>::size_type
2452 PreconditionRelaxation<MatrixType, PreconditionerType>::m() const
2453{
2454 Assert(A != nullptr, ExcNotInitialized());
2455 return A->m();
2456}
2457
2458template <typename MatrixType, typename PreconditionerType>
2459inline
2460 typename PreconditionRelaxation<MatrixType, PreconditionerType>::size_type
2461 PreconditionRelaxation<MatrixType, PreconditionerType>::n() const
2462{
2463 Assert(A != nullptr, ExcNotInitialized());
2464 return A->n();
2465}
2466
2467template <typename MatrixType, typename PreconditionerType>
2468template <class VectorType>
2469inline void
2470PreconditionRelaxation<MatrixType, PreconditionerType>::vmult(
2471 VectorType & dst,
2472 const VectorType &src) const
2473{
2474 Assert(this->A != nullptr, ExcNotInitialized());
2475 Assert(this->preconditioner != nullptr, ExcNotInitialized());
2476
2477 VectorType tmp1, tmp2;
2478
2479 for (unsigned int i = 0; i < n_iterations; ++i)
2480 internal::PreconditionRelaxation::step_operations(
2481 *A, *preconditioner, dst, src, relaxation, tmp1, tmp2, i, false);
2482}
2483
2484template <typename MatrixType, typename PreconditionerType>
2485template <class VectorType>
2486inline void
2487PreconditionRelaxation<MatrixType, PreconditionerType>::Tvmult(
2488 VectorType & dst,
2489 const VectorType &src) const
2490{
2491 Assert(this->A != nullptr, ExcNotInitialized());
2492 Assert(this->preconditioner != nullptr, ExcNotInitialized());
2493
2494 VectorType tmp1, tmp2;
2495
2496 for (unsigned int i = 0; i < n_iterations; ++i)
2497 internal::PreconditionRelaxation::step_operations(
2498 *A, *preconditioner, dst, src, relaxation, tmp1, tmp2, i, true);
2499}
2500
2501template <typename MatrixType, typename PreconditionerType>
2502template <class VectorType>
2503inline void
2504PreconditionRelaxation<MatrixType, PreconditionerType>::step(
2505 VectorType & dst,
2506 const VectorType &src) const
2507{
2508 Assert(this->A != nullptr, ExcNotInitialized());
2509 Assert(this->preconditioner != nullptr, ExcNotInitialized());
2510
2511 VectorType tmp1, tmp2;
2512
2513 for (unsigned int i = 1; i <= n_iterations; ++i)
2514 internal::PreconditionRelaxation::step_operations(
2515 *A, *preconditioner, dst, src, relaxation, tmp1, tmp2, i, false);
2516}
2517
2518template <typename MatrixType, typename PreconditionerType>
2519template <class VectorType>
2520inline void
2521PreconditionRelaxation<MatrixType, PreconditionerType>::Tstep(
2522 VectorType & dst,
2523 const VectorType &src) const
2524{
2525 Assert(this->A != nullptr, ExcNotInitialized());
2526 Assert(this->preconditioner != nullptr, ExcNotInitialized());
2527
2528 VectorType tmp1, tmp2;
2529
2530 for (unsigned int i = 1; i <= n_iterations; ++i)
2531 internal::PreconditionRelaxation::step_operations(
2532 *A, *preconditioner, dst, src, relaxation, tmp1, tmp2, i, true);
2533}
2534
2535//---------------------------------------------------------------------------
2536
2537template <typename MatrixType>
2538inline void
2539PreconditionJacobi<MatrixType>::initialize(const MatrixType & A,
2540 const AdditionalData &parameters_in)
2541{
2542 Assert(parameters_in.preconditioner == nullptr, ExcInternalError());
2543
2544 AdditionalData parameters;
2545 parameters.relaxation = 1.0;
2546 parameters.n_iterations = parameters_in.n_iterations;
2547 parameters.preconditioner =
2548 std::make_shared<PreconditionerType>(A, parameters_in.relaxation);
2549
2550 this->BaseClass::initialize(A, parameters);
2551}
2552
2553//---------------------------------------------------------------------------
2554
2555template <typename MatrixType>
2556inline void
2557PreconditionSOR<MatrixType>::initialize(const MatrixType & A,
2558 const AdditionalData &parameters_in)
2559{
2560 Assert(parameters_in.preconditioner == nullptr, ExcInternalError());
2561
2562 AdditionalData parameters;
2563 parameters.relaxation = 1.0;
2564 parameters.n_iterations = parameters_in.n_iterations;
2565 parameters.preconditioner =
2566 std::make_shared<PreconditionerType>(A, parameters_in.relaxation);
2567
2568 this->BaseClass::initialize(A, parameters);
2569}
2570
2571//---------------------------------------------------------------------------
2572
2573template <typename MatrixType>
2574inline void
2575PreconditionSSOR<MatrixType>::initialize(const MatrixType & A,
2576 const AdditionalData &parameters_in)
2577{
2578 Assert(parameters_in.preconditioner == nullptr, ExcInternalError());
2579
2580 AdditionalData parameters;
2581 parameters.relaxation = 1.0;
2582 parameters.n_iterations = parameters_in.n_iterations;
2583 parameters.preconditioner =
2584 std::make_shared<PreconditionerType>(A, parameters_in.relaxation);
2585
2586 this->BaseClass::initialize(A, parameters);
2587}
2588
2589
2590
2591//---------------------------------------------------------------------------
2592
2593template <typename MatrixType>
2594inline void
2595PreconditionPSOR<MatrixType>::initialize(
2596 const MatrixType & A,
2597 const std::vector<size_type> & p,
2598 const std::vector<size_type> & ip,
2599 const typename BaseClass::AdditionalData &parameters_in)
2600{
2601 Assert(parameters_in.preconditioner == nullptr, ExcInternalError());
2602
2603 typename BaseClass::AdditionalData parameters;
2604 parameters.relaxation = 1.0;
2605 parameters.n_iterations = parameters_in.n_iterations;
2606 parameters.preconditioner =
2607 std::make_shared<PreconditionerType>(A, parameters_in.relaxation, p, ip);
2608
2609 this->BaseClass::initialize(A, parameters);
2610}
2611
2612
2613template <typename MatrixType>
2614inline void
2615PreconditionPSOR<MatrixType>::initialize(const MatrixType & A,
2616 const AdditionalData &additional_data)
2617{
2618 initialize(A,
2619 additional_data.permutation,
2620 additional_data.inverse_permutation,
2621 additional_data.parameters);
2622}
2623
2624template <typename MatrixType>
2625PreconditionPSOR<MatrixType>::AdditionalData::AdditionalData(
2626 const std::vector<size_type> &permutation,
2627 const std::vector<size_type> &inverse_permutation,
2628 const typename PreconditionPSOR<MatrixType>::BaseClass::AdditionalData
2629 &parameters)
2630 : permutation(permutation)
2631 , inverse_permutation(inverse_permutation)
2632 , parameters(parameters)
2633{}
2634
2635
2636//---------------------------------------------------------------------------
2637
2638
2639template <typename MatrixType, class VectorType>
2640PreconditionUseMatrix<MatrixType, VectorType>::PreconditionUseMatrix(
2641 const MatrixType & M,
2642 const function_ptr method)
2643 : matrix(M)
2644 , precondition(method)
2645{}
2646
2647
2648
2649template <typename MatrixType, class VectorType>
2650void
2651PreconditionUseMatrix<MatrixType, VectorType>::vmult(
2652 VectorType & dst,
2653 const VectorType &src) const
2654{
2655 (matrix.*precondition)(dst, src);
2656}
2657
2658//---------------------------------------------------------------------------
2659
2660template <typename MatrixType, typename PreconditionerType>
2661inline PreconditionRelaxation<MatrixType, PreconditionerType>::AdditionalData::
2662 AdditionalData(const double relaxation, const unsigned int n_iterations)
2663 : relaxation(relaxation)
2664 , n_iterations(n_iterations)
2665{}
2666
2667
2668
2669//---------------------------------------------------------------------------
2670
2671namespace internal
2672{
2673 namespace PreconditionChebyshevImplementation
2674 {
2675 // for deal.II vectors, perform updates for Chebyshev preconditioner all
2676 // at once to reduce memory transfer. Here, we select between general
2677 // vectors and deal.II vectors where we expand the loop over the (local)
2678 // size of the vector
2679
2680 // generic part for non-deal.II vectors
2681 template <typename VectorType, typename PreconditionerType>
2682 inline void
2683 vector_updates(const VectorType & rhs,
2684 const PreconditionerType &preconditioner,
2685 const unsigned int iteration_index,
2686 const double factor1,
2687 const double factor2,
2688 VectorType & solution_old,
2689 VectorType & temp_vector1,
2690 VectorType & temp_vector2,
2691 VectorType & solution)
2692 {
2693 if (iteration_index == 0)
2694 {
2695 solution.equ(factor2, rhs);
2696 preconditioner.vmult(solution_old, solution);
2697 }
2698 else if (iteration_index == 1)
2699 {
2700 // compute t = P^{-1} * (b-A*x^{n})
2701 temp_vector1.sadd(-1.0, 1.0, rhs);
2702 preconditioner.vmult(solution_old, temp_vector1);
2703
2704 // compute x^{n+1} = x^{n} + f_1 * x^{n} + f_2 * t
2705 solution_old.sadd(factor2, 1 + factor1, solution);
2706 }
2707 else
2708 {
2709 // compute t = P^{-1} * (b-A*x^{n})
2710 temp_vector1.sadd(-1.0, 1.0, rhs);
2711 preconditioner.vmult(temp_vector2, temp_vector1);
2712
2713 // compute x^{n+1} = x^{n} + f_1 * (x^{n}-x^{n-1}) + f_2 * t
2714 solution_old.sadd(-factor1, factor2, temp_vector2);
2715 solution_old.add(1 + factor1, solution);
2716 }
2717
2718 solution.swap(solution_old);
2719 }
2720
2721 // generic part for deal.II vectors
2722 template <
2723 typename Number,
2724 typename PreconditionerType,
2725 std::enable_if_t<
2726 !has_vmult_with_std_functions_for_precondition<
2727 PreconditionerType,
2728 LinearAlgebra::distributed::Vector<Number, MemorySpace::Host>>,
2729 int> * = nullptr>
2730 inline void
2731 vector_updates(
2732 const LinearAlgebra::distributed::Vector<Number, MemorySpace::Host> &rhs,
2733 const PreconditionerType &preconditioner,
2734 const unsigned int iteration_index,
2735 const double factor1_,
2736 const double factor2_,
2737 LinearAlgebra::distributed::Vector<Number, MemorySpace::Host>
2738 &solution_old,
2739 LinearAlgebra::distributed::Vector<Number, MemorySpace::Host>
2740 &temp_vector1,
2741 LinearAlgebra::distributed::Vector<Number, MemorySpace::Host>
2742 &temp_vector2,
2743 LinearAlgebra::distributed::Vector<Number, MemorySpace::Host> &solution)
2744 {
2745 const Number factor1 = factor1_;
2746 const Number factor1_plus_1 = 1. + factor1_;
2747 const Number factor2 = factor2_;
2748
2749 if (iteration_index == 0)
2750 {
2751 const auto solution_old_ptr = solution_old.begin();
2752
2753 // compute t = P^{-1} * (b)
2754 preconditioner.vmult(solution_old, rhs);
2755
2756 // compute x^{n+1} = f_2 * t
2757 DEAL_II_OPENMP_SIMD_PRAGMA
2758 for (unsigned int i = 0; i < solution_old.locally_owned_size(); ++i)
2759 solution_old_ptr[i] = solution_old_ptr[i] * factor2;
2760 }
2761 else if (iteration_index == 1)
2762 {
2763 const auto solution_ptr = solution.begin();
2764 const auto solution_old_ptr = solution_old.begin();
2765
2766 // compute t = P^{-1} * (b-A*x^{n})
2767 temp_vector1.sadd(-1.0, 1.0, rhs);
2768
2769 preconditioner.vmult(solution_old, temp_vector1);
2770
2771 // compute x^{n+1} = x^{n} + f_1 * x^{n} + f_2 * t
2772 DEAL_II_OPENMP_SIMD_PRAGMA
2773 for (unsigned int i = 0; i < solution_old.locally_owned_size(); ++i)
2774 solution_old_ptr[i] =
2775 factor1_plus_1 * solution_ptr[i] + solution_old_ptr[i] * factor2;
2776 }
2777 else
2778 {
2779 const auto solution_ptr = solution.begin();
2780 const auto solution_old_ptr = solution_old.begin();
2781 const auto temp_vector2_ptr = temp_vector2.begin();
2782
2783 // compute t = P^{-1} * (b-A*x^{n})
2784 temp_vector1.sadd(-1.0, 1.0, rhs);
2785
2786 preconditioner.vmult(temp_vector2, temp_vector1);
2787
2788 // compute x^{n+1} = x^{n} + f_1 * (x^{n}-x^{n-1}) + f_2 * t
2789 DEAL_II_OPENMP_SIMD_PRAGMA
2790 for (unsigned int i = 0; i < solution_old.locally_owned_size(); ++i)
2791 solution_old_ptr[i] = factor1_plus_1 * solution_ptr[i] -
2792 factor1 * solution_old_ptr[i] +
2793 temp_vector2_ptr[i] * factor2;
2794 }
2795
2796 solution.swap(solution_old);
2797 }
2798
2799 template <
2800 typename Number,
2801 typename PreconditionerType,
2802 std::enable_if_t<
2803 has_vmult_with_std_functions_for_precondition<
2804 PreconditionerType,
2805 LinearAlgebra::distributed::Vector<Number, MemorySpace::Host>>,
2806 int> * = nullptr>
2807 inline void
2808 vector_updates(
2809 const LinearAlgebra::distributed::Vector<Number, MemorySpace::Host> &rhs,
2810 const PreconditionerType &preconditioner,
2811 const unsigned int iteration_index,
2812 const double factor1_,
2813 const double factor2_,
2814 LinearAlgebra::distributed::Vector<Number, MemorySpace::Host>
2815 &solution_old,
2816 LinearAlgebra::distributed::Vector<Number, MemorySpace::Host>
2817 &temp_vector1,
2818 LinearAlgebra::distributed::Vector<Number, MemorySpace::Host>
2819 &temp_vector2,
2820 LinearAlgebra::distributed::Vector<Number, MemorySpace::Host> &solution)
2821 {
2822 const Number factor1 = factor1_;
2823 const Number factor1_plus_1 = 1. + factor1_;
2824 const Number factor2 = factor2_;
2825
2826 const auto rhs_ptr = rhs.begin();
2827 const auto temp_vector1_ptr = temp_vector1.begin();
2828 const auto temp_vector2_ptr = temp_vector2.begin();
2829 const auto solution_ptr = solution.begin();
2830 const auto solution_old_ptr = solution_old.begin();
2831
2832 if (iteration_index == 0)
2833 {
2834 preconditioner.vmult(
2835 solution,
2836 rhs,
2837 [&](const auto start_range, const auto end_range) {
2838 if (end_range > start_range)
2839 std::memset(solution.begin() + start_range,
2840 0,
2841 sizeof(Number) * (end_range - start_range));
2842 },
2843 [&](const auto begin, const auto end) {
2844 DEAL_II_OPENMP_SIMD_PRAGMA
2845 for (std::size_t i = begin; i < end; ++i)
2846 solution_ptr[i] *= factor2;
2847 });
2848 }
2849 else
2850 {
2851 preconditioner.vmult(
2852 temp_vector2,
2853 temp_vector1,
2854 [&](const auto begin, const auto end) {
2855 if (end > begin)
2856 std::memset(temp_vector2.begin() + begin,
2857 0,
2858 sizeof(Number) * (end - begin));
2859
2860 DEAL_II_OPENMP_SIMD_PRAGMA
2861 for (std::size_t i = begin; i < end; ++i)
2862 temp_vector1_ptr[i] = rhs_ptr[i] - temp_vector1_ptr[i];
2863 },
2864 [&](const auto begin, const auto end) {
2865 if (iteration_index == 1)
2866 {
2867 DEAL_II_OPENMP_SIMD_PRAGMA
2868 for (std::size_t i = begin; i < end; ++i)
2869 temp_vector2_ptr[i] = factor1_plus_1 * solution_ptr[i] +
2870 factor2 * temp_vector2_ptr[i];
2871 }
2872 else
2873 {
2874 DEAL_II_OPENMP_SIMD_PRAGMA
2875 for (std::size_t i = begin; i < end; ++i)
2876 temp_vector2_ptr[i] = factor1_plus_1 * solution_ptr[i] -
2877 factor1 * solution_old_ptr[i] +
2878 factor2 * temp_vector2_ptr[i];
2879 }
2880 });
2881 }
2882
2883 if (iteration_index > 0)
2884 {
2885 solution_old.swap(temp_vector2);
2886 solution_old.swap(solution);
2887 }
2888 }
2889
2890 // worker routine for deal.II vectors. Because of vectorization, we need
2891 // to put the loop into an extra structure because the virtual function of
2892 // VectorUpdatesRange prevents the compiler from applying vectorization.
2893 template <typename Number>
2894 struct VectorUpdater
2895 {
2896 VectorUpdater(const Number * rhs,
2897 const Number * matrix_diagonal_inverse,
2898 const unsigned int iteration_index,
2899 const Number factor1,
2900 const Number factor2,
2901 Number * solution_old,
2902 Number * tmp_vector,
2903 Number * solution)
2904 : rhs(rhs)
2905 , matrix_diagonal_inverse(matrix_diagonal_inverse)
2906 , iteration_index(iteration_index)
2907 , factor1(factor1)
2908 , factor2(factor2)
2909 , solution_old(solution_old)
2910 , tmp_vector(tmp_vector)
2911 , solution(solution)
2912 {}
2913
2914 void
2915 apply_to_subrange(const std::size_t begin, const std::size_t end) const
2916 {
2917 // To circumvent a bug in gcc
2918 // (https://gcc.gnu.org/bugzilla/show_bug.cgi?id=63945), we create
2919 // copies of the variables factor1 and factor2 and do not check based on
2920 // factor1.
2921 const Number factor1 = this->factor1;
2922 const Number factor1_plus_1 = 1. + this->factor1;
2923 const Number factor2 = this->factor2;
2924 if (iteration_index == 0)
2925 {
2926 DEAL_II_OPENMP_SIMD_PRAGMA
2927 for (std::size_t i = begin; i < end; ++i)
2928 solution[i] = factor2 * matrix_diagonal_inverse[i] * rhs[i];
2929 }
2930 else if (iteration_index == 1)
2931 {
2932 // x^{n+1} = x^{n} + f_1 * x^{n} + f_2 * P^{-1} * (b-A*x^{n})
2933 DEAL_II_OPENMP_SIMD_PRAGMA
2934 for (std::size_t i = begin; i < end; ++i)
2935 // for efficiency reason, write back to temp_vector that is
2936 // already read (avoid read-for-ownership)
2937 tmp_vector[i] =
2938 factor1_plus_1 * solution[i] +
2939 factor2 * matrix_diagonal_inverse[i] * (rhs[i] - tmp_vector[i]);
2940 }
2941 else
2942 {
2943 // x^{n+1} = x^{n} + f_1 * (x^{n}-x^{n-1})
2944 // + f_2 * P^{-1} * (b-A*x^{n})
2945 DEAL_II_OPENMP_SIMD_PRAGMA
2946 for (std::size_t i = begin; i < end; ++i)
2947 // for efficiency reason, write back to temp_vector, which is
2948 // already modified during vmult (in best case, the modified
2949 // values are not written back to main memory yet so that
2950 // we do not have to pay additional costs for writing and
2951 // read-for-ownershop)
2952 tmp_vector[i] =
2953 factor1_plus_1 * solution[i] - factor1 * solution_old[i] +
2954 factor2 * matrix_diagonal_inverse[i] * (rhs[i] - tmp_vector[i]);
2955 }
2956 }
2957
2958 const Number * rhs;
2959 const Number * matrix_diagonal_inverse;
2960 const unsigned int iteration_index;
2961 const Number factor1;
2962 const Number factor2;
2963 mutable Number * solution_old;
2964 mutable Number * tmp_vector;
2965 mutable Number * solution;
2966 };
2967
2968 template <typename Number>
2969 struct VectorUpdatesRange : public ::parallel::ParallelForInteger
2970 {
2971 VectorUpdatesRange(const VectorUpdater<Number> &updater,
2972 const std::size_t size)
2973 : updater(updater)
2974 {
2975 if (size < internal::VectorImplementation::minimum_parallel_grain_size)
2976 VectorUpdatesRange::apply_to_subrange(0, size);
2977 else
2978 apply_parallel(
2979 0,
2980 size,
2981 internal::VectorImplementation::minimum_parallel_grain_size);
2982 }
2983
2984 ~VectorUpdatesRange() override = default;
2985
2986 virtual void
2987 apply_to_subrange(const std::size_t begin,
2988 const std::size_t end) const override
2989 {
2990 updater.apply_to_subrange(begin, end);
2991 }
2992
2993 const VectorUpdater<Number> &updater;
2994 };
2995
2996 // selection for diagonal matrix around deal.II vector
2997 template <typename Number>
2998 inline void
2999 vector_updates(
3000 const ::Vector<Number> & rhs,
3001 const ::DiagonalMatrix<::Vector<Number>> &jacobi,
3002 const unsigned int iteration_index,
3003 const double factor1,
3004 const double factor2,
3005 ::Vector<Number> & solution_old,
3006 ::Vector<Number> & temp_vector1,
3007 ::Vector<Number> &,
3008 ::Vector<Number> &solution)
3009 {
3010 VectorUpdater<Number> upd(rhs.begin(),
3011 jacobi.get_vector().begin(),
3012 iteration_index,
3013 factor1,
3014 factor2,
3015 solution_old.begin(),
3016 temp_vector1.begin(),
3017 solution.begin());
3018 VectorUpdatesRange<Number>(upd, rhs.size());
3019
3020 // swap vectors x^{n+1}->x^{n}, given the updates in the function above
3021 if (iteration_index == 0)
3022 {
3023 // nothing to do here because we can immediately write into the
3024 // solution vector without remembering any of the other vectors
3025 }
3026 else
3027 {
3028 solution.swap(temp_vector1);
3029 solution_old.swap(temp_vector1);
3030 }
3031 }
3032
3033 // selection for diagonal matrix around parallel deal.II vector
3034 template <typename Number>
3035 inline void
3036 vector_updates(
3037 const LinearAlgebra::distributed::Vector<Number, MemorySpace::Host> &rhs,
3038 const ::DiagonalMatrix<
3039 LinearAlgebra::distributed::Vector<Number, MemorySpace::Host>> &jacobi,
3040 const unsigned int iteration_index,
3041 const double factor1,
3042 const double factor2,
3043 LinearAlgebra::distributed::Vector<Number, MemorySpace::Host>
3044 &solution_old,
3045 LinearAlgebra::distributed::Vector<Number, MemorySpace::Host>
3046 &temp_vector1,
3047 LinearAlgebra::distributed::Vector<Number, MemorySpace::Host> &,
3048 LinearAlgebra::distributed::Vector<Number, MemorySpace::Host> &solution)
3049 {
3050 VectorUpdater<Number> upd(rhs.begin(),
3051 jacobi.get_vector().begin(),
3052 iteration_index,
3053 factor1,
3054 factor2,
3055 solution_old.begin(),
3056 temp_vector1.begin(),
3057 solution.begin());
3058 VectorUpdatesRange<Number>(upd, rhs.locally_owned_size());
3059
3060 // swap vectors x^{n+1}->x^{n}, given the updates in the function above
3061 if (iteration_index == 0)
3062 {
3063 // nothing to do here because we can immediately write into the
3064 // solution vector without remembering any of the other vectors
3065 }
3066 else
3067 {
3068 solution.swap(temp_vector1);
3069 solution_old.swap(temp_vector1);
3070 }
3071 }
3072
3073 // We need to have a separate declaration for static const members
3074
3075 // general case and the case that the preconditioner can work on
3076 // ranges (covered by vector_updates())
3077 template <
3078 typename MatrixType,
3079 typename VectorType,
3080 typename PreconditionerType,
3081 std::enable_if_t<
3082 !has_vmult_with_std_functions<MatrixType,
3083 VectorType,
3084 PreconditionerType> &&
3085 !(has_vmult_with_std_functions_for_precondition<PreconditionerType,
3086 VectorType> &&
3087 has_vmult_with_std_functions_for_precondition<MatrixType,
3088 VectorType>),
3089 int> * = nullptr>
3090 inline void
3091 vmult_and_update(const MatrixType & matrix,
3092 const PreconditionerType &preconditioner,
3093 const VectorType & rhs,
3094 const unsigned int iteration_index,
3095 const double factor1,
3096 const double factor2,
3097 VectorType & solution,
3098 VectorType & solution_old,
3099 VectorType & temp_vector1,
3100 VectorType & temp_vector2)
3101 {
3102 if (iteration_index > 0)
3103 matrix.vmult(temp_vector1, solution);
3104 vector_updates(rhs,
3105 preconditioner,
3106 iteration_index,
3107 factor1,
3108 factor2,
3109 solution_old,
3110 temp_vector1,
3111 temp_vector2,
3112 solution);
3113 }
3114
3115 // case that both the operator and the preconditioner can work on
3116 // subranges
3117 template <
3118 typename MatrixType,
3119 typename VectorType,
3120 typename PreconditionerType,
3121 std::enable_if_t<
3122 !has_vmult_with_std_functions<MatrixType,
3123 VectorType,
3124 PreconditionerType> &&
3125 (has_vmult_with_std_functions_for_precondition<PreconditionerType,
3126 VectorType> &&
3127 has_vmult_with_std_functions_for_precondition<MatrixType,
3128 VectorType>),
3129 int> * = nullptr>
3130 inline void
3131 vmult_and_update(const MatrixType & matrix,
3132 const PreconditionerType &preconditioner,
3133 const VectorType & rhs,
3134 const unsigned int iteration_index,
3135 const double factor1_,
3136 const double factor2_,
3137 VectorType & solution,
3138 VectorType & solution_old,
3139 VectorType & temp_vector1,
3140 VectorType & temp_vector2)
3141 {
3142 using Number = typename VectorType::value_type;
3143
3144 const Number factor1 = factor1_;
3145 const Number factor1_plus_1 = 1. + factor1_;
3146 const Number factor2 = factor2_;
3147
3148 if (iteration_index == 0)
3149 {
3150 preconditioner.vmult(
3151 solution,
3152 rhs,
3153 [&](const unsigned int start_range, const unsigned int end_range) {
3154 // zero 'solution' before running the vmult operation
3155 if (end_range > start_range)
3156 std::memset(solution.begin() + start_range,
3157 0,
3158 sizeof(Number) * (end_range - start_range));
3159 },
3160 [&](const unsigned int start_range, const unsigned int end_range) {
3161 const auto solution_ptr = solution.begin();
3162
3163 DEAL_II_OPENMP_SIMD_PRAGMA
3164 for (std::size_t i = start_range; i < end_range; ++i)
3165 solution_ptr[i] *= factor2;
3166 });
3167 }
3168 else
3169 {
3170 temp_vector1.reinit(rhs, true);
3171 temp_vector2.reinit(rhs, true);
3172
3173 // 1) compute rediduum (including operator application)
3174 matrix.vmult(
3175 temp_vector1,
3176 solution,
3177 [&](const unsigned int start_range, const unsigned int end_range) {
3178 // zero 'temp_vector1' before running the vmult
3179 // operation
3180 if (end_range > start_range)
3181 std::memset(temp_vector1.begin() + start_range,
3182 0,
3183 sizeof(Number) * (end_range - start_range));
3184 },
3185 [&](const unsigned int start_range, const unsigned int end_range) {
3186 const auto rhs_ptr = rhs.begin();
3187 const auto tmp_ptr = temp_vector1.begin();
3188
3189 DEAL_II_OPENMP_SIMD_PRAGMA
3190 for (std::size_t i = start_range; i < end_range; ++i)
3191 tmp_ptr[i] = rhs_ptr[i] - tmp_ptr[i];
3192 });
3193
3194 // 2) perform vector updates (including preconditioner application)
3195 preconditioner.vmult(
3196 temp_vector2,
3197 temp_vector1,
3198 [&](const unsigned int start_range, const unsigned int end_range) {
3199 // zero 'temp_vector2' before running the vmult
3200 // operation
3201 if (end_range > start_range)
3202 std::memset(temp_vector2.begin() + start_range,
3203 0,
3204 sizeof(Number) * (end_range - start_range));
3205 },
3206 [&](const unsigned int start_range, const unsigned int end_range) {
3207 const auto solution_ptr = solution.begin();
3208 const auto solution_old_ptr = solution_old.begin();
3209 const auto tmp_ptr = temp_vector2.begin();
3210
3211 if (iteration_index == 1)
3212 {
3213 DEAL_II_OPENMP_SIMD_PRAGMA
3214 for (std::size_t i = start_range; i < end_range; ++i)
3215 tmp_ptr[i] =
3216 factor1_plus_1 * solution_ptr[i] + factor2 * tmp_ptr[i];
3217 }
3218 else
3219 {
3220 DEAL_II_OPENMP_SIMD_PRAGMA
3221 for (std::size_t i = start_range; i < end_range; ++i)
3222 tmp_ptr[i] = factor1_plus_1 * solution_ptr[i] -
3223 factor1 * solution_old_ptr[i] +
3224 factor2 * tmp_ptr[i];
3225 }
3226 });
3227
3228 solution.swap(temp_vector2);
3229 solution_old.swap(temp_vector2);
3230 }
3231 }
3232
3233 // case that the operator can work on subranges and the preconditioner
3234 // is a diagonal
3235 template <typename MatrixType,
3236 typename VectorType,
3237 typename PreconditionerType,
3238 std::enable_if_t<has_vmult_with_std_functions<MatrixType,
3239 VectorType,
3240 PreconditionerType>,
3241 int> * = nullptr>
3242 inline void
3243 vmult_and_update(const MatrixType & matrix,
3244 const PreconditionerType &preconditioner,
3245 const VectorType & rhs,
3246 const unsigned int iteration_index,
3247 const double factor1,
3248 const double factor2,
3249 VectorType & solution,
3250 VectorType & solution_old,
3251 VectorType & temp_vector1,
3252 VectorType &)
3253 {
3254 using Number = typename VectorType::value_type;
3255 VectorUpdater<Number> updater(rhs.begin(),
3256 preconditioner.get_vector().begin(),
3257 iteration_index,
3258 factor1,
3259 factor2,
3260 solution_old.begin(),
3261 temp_vector1.begin(),
3262 solution.begin());
3263 if (iteration_index > 0)
3264 matrix.vmult(
3265 temp_vector1,
3266 solution,
3267 [&](const unsigned int start_range, const unsigned int end_range) {
3268 // zero 'temp_vector1' before running the vmult
3269 // operation
3270 if (end_range > start_range)
3271 std::memset(temp_vector1.begin() + start_range,
3272 0,
3273 sizeof(Number) * (end_range - start_range));
3274 },
3275 [&](const unsigned int start_range, const unsigned int end_range) {
3276 if (end_range > start_range)
3277 updater.apply_to_subrange(start_range, end_range);
3278 });
3279 else
3280 updater.apply_to_subrange(0U, solution.locally_owned_size());
3281
3282 // swap vectors x^{n+1}->x^{n}, given the updates in the function above
3283 if (iteration_index == 0)
3284 {
3285 // nothing to do here because we can immediately write into the
3286 // solution vector without remembering any of the other vectors
3287 }
3288 else
3289 {
3290 solution.swap(temp_vector1);
3291 solution_old.swap(temp_vector1);
3292 }
3293 }
3294
3295 template <typename MatrixType, typename PreconditionerType>
3296 inline void
3297 initialize_preconditioner(
3298 const MatrixType & matrix,
3299 std::shared_ptr<PreconditionerType> &preconditioner)
3300 {
3301 (void)matrix;
3302 (void)preconditioner;
3303 AssertThrow(preconditioner.get() != nullptr, ExcNotInitialized());
3304 }
3305
3306 template <typename MatrixType, typename VectorType>
3307 inline void
3308 initialize_preconditioner(
3309 const MatrixType & matrix,
3310 std::shared_ptr<::DiagonalMatrix<VectorType>> &preconditioner)
3311 {
3312 if (preconditioner.get() == nullptr || preconditioner->m() != matrix.m())
3313 {
3314 if (preconditioner.get() == nullptr)
3315 preconditioner =
3316 std::make_shared<::DiagonalMatrix<VectorType>>();
3317
3318 Assert(
3319 preconditioner->m() == 0,
3320 ExcMessage(
3321 "Preconditioner appears to be initialized but not sized correctly"));
3322
3323 // This part only works in serial
3324 if (preconditioner->m() != matrix.m())
3325 {
3326 preconditioner->get_vector().reinit(matrix.m());
3327 for (typename VectorType::size_type i = 0; i < matrix.m(); ++i)
3328 preconditioner->get_vector()(i) = 1. / matrix.el(i, i);
3329 }
3330 }
3331 }
3332
3333 template <typename VectorType>
3334 void
3335 set_initial_guess(VectorType &vector)
3336 {
3337 vector = 1. / std::sqrt(static_cast<double>(vector.size()));
3338 if (vector.locally_owned_elements().is_element(0))
3339 vector(0) = 0.;
3340 }
3341
3342 template <typename Number>
3343 void
3344 set_initial_guess(::Vector<Number> &vector)
3345 {
3346 // Choose a high-frequency mode consisting of numbers between 0 and 1
3347 // that is cheap to compute (cheaper than random numbers) but avoids
3348 // obviously re-occurring numbers in multi-component systems by choosing
3349 // a period of 11
3350 for (unsigned int i = 0; i < vector.size(); ++i)
3351 vector(i) = i % 11;
3352
3353 const Number mean_value = vector.mean_value();
3354 vector.add(-mean_value);
3355 }
3356
3357 template <typename Number>
3358 void
3359 set_initial_guess(
3360 ::LinearAlgebra::distributed::BlockVector<Number> &vector)
3361 {
3362 for (unsigned int block = 0; block < vector.n_blocks(); ++block)
3363 set_initial_guess(vector.block(block));
3364 }
3365
3366 template <typename Number, typename MemorySpace>
3367 void
3368 set_initial_guess(
3369 ::LinearAlgebra::distributed::Vector<Number, MemorySpace> &vector)
3370 {
3371 // Choose a high-frequency mode consisting of numbers between 0 and 1
3372 // that is cheap to compute (cheaper than random numbers) but avoids
3373 // obviously re-occurring numbers in multi-component systems by choosing
3374 // a period of 11.
3375 // Make initial guess robust with respect to number of processors
3376 // by operating on the global index.
3377 types::global_dof_index first_local_range = 0;
3378 if (!vector.locally_owned_elements().is_empty())
3379 first_local_range = vector.locally_owned_elements().nth_index_in_set(0);
3380
3381 const auto n_local_elements = vector.locally_owned_size();
3382 Number * values_ptr = vector.get_values();
3383 Kokkos::RangePolicy<typename MemorySpace::kokkos_space::execution_space,
3384 Kokkos::IndexType<types::global_dof_index>>
3385 policy(0, n_local_elements);
3386 Kokkos::parallel_for(
3387 "::PreconditionChebyshev::set_initial_guess",
3388 policy,
3389 KOKKOS_LAMBDA(types::global_dof_index i) {
3390 values_ptr[i] = (i + first_local_range) % 11;
3391 });
3392 const Number mean_value = vector.mean_value();
3393 vector.add(-mean_value);
3394 }
3395
3396 struct EigenvalueTracker
3397 {
3398 public:
3399 void
3400 slot(const std::vector<double> &eigenvalues)
3401 {
3402 values = eigenvalues;
3403 }
3404
3405 std::vector<double> values;
3406 };
3407
3408
3409
3410 template <typename MatrixType,
3411 typename VectorType,
3412 typename PreconditionerType>
3413 double
3414 power_iteration(const MatrixType & matrix,
3415 VectorType & eigenvector,
3416 const PreconditionerType &preconditioner,
3417 const unsigned int n_iterations)
3418 {
3419 double eigenvalue_estimate = 0.;
3420 eigenvector /= eigenvector.l2_norm();
3421 VectorType vector1, vector2;
3422 vector1.reinit(eigenvector, true);
3423 if (!std::is_same<PreconditionerType, PreconditionIdentity>::value)
3424 vector2.reinit(eigenvector, true);
3425 for (unsigned int i = 0; i < n_iterations; ++i)
3426 {
3427 if (!std::is_same<PreconditionerType, PreconditionIdentity>::value)
3428 {
3429 matrix.vmult(vector2, eigenvector);
3430 preconditioner.vmult(vector1, vector2);
3431 }
3432 else
3433 matrix.vmult(vector1, eigenvector);
3434
3435 eigenvalue_estimate = eigenvector * vector1;
3436
3437 vector1 /= vector1.l2_norm();
3438 eigenvector.swap(vector1);
3439 }
3440 return eigenvalue_estimate;
3441 }
3442 } // namespace PreconditionChebyshevImplementation
3443} // namespace internal
3444
3445
3446
3447template <typename MatrixType, class VectorType, typename PreconditionerType>
3448inline PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::
3449 AdditionalData::AdditionalData(const unsigned int degree,
3450 const double smoothing_range,
3451 const unsigned int eig_cg_n_iterations,
3452 const double eig_cg_residual,
3453 const double max_eigenvalue,
3454 const EigenvalueAlgorithm eigenvalue_algorithm,
3455 const PolynomialType polynomial_type)
3456 : degree(degree)
3457 , smoothing_range(smoothing_range)
3458 , eig_cg_n_iterations(eig_cg_n_iterations)
3459 , eig_cg_residual(eig_cg_residual)
3460 , max_eigenvalue(max_eigenvalue)
3461 , eigenvalue_algorithm(eigenvalue_algorithm)
3462 , polynomial_type(polynomial_type)
3463{}
3464
3465
3466
3467template <typename MatrixType, class VectorType, typename PreconditionerType>
3468inline typename PreconditionChebyshev<MatrixType,
3469 VectorType,
3470 PreconditionerType>::AdditionalData &
3471PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::
3472 AdditionalData::operator=(const AdditionalData &other_data)
3473{
3474 degree = other_data.degree;
3475 smoothing_range = other_data.smoothing_range;
3476 eig_cg_n_iterations = other_data.eig_cg_n_iterations;
3477 eig_cg_residual = other_data.eig_cg_residual;
3478 max_eigenvalue = other_data.max_eigenvalue;
3479 preconditioner = other_data.preconditioner;
3480 eigenvalue_algorithm = other_data.eigenvalue_algorithm;
3481 polynomial_type = other_data.polynomial_type;
3482 constraints.copy_from(other_data.constraints);
3483
3484 return *this;
3485}
3486
3487
3488
3489template <typename MatrixType, typename VectorType, typename PreconditionerType>
3490inline PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::
3491 PreconditionChebyshev()
3492 : theta(1.)
3493 , delta(1.)
3494 , eigenvalues_are_initialized(false)
3495{
3496 static_assert(
3497 std::is_same<size_type, typename VectorType::size_type>::value,
3498 "PreconditionChebyshev and VectorType must have the same size_type.");
3499}
3500
3501
3502
3503template <typename MatrixType, typename VectorType, typename PreconditionerType>
3504inline void
3505PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::initialize(
3506 const MatrixType & matrix,
3507 const AdditionalData &additional_data)
3508{
3509 matrix_ptr = &matrix;
3510 data = additional_data;
3511 Assert(data.degree > 0,
3512 ExcMessage("The degree of the Chebyshev method must be positive."));
3513 internal::PreconditionChebyshevImplementation::initialize_preconditioner(
3514 matrix, data.preconditioner);
3515 eigenvalues_are_initialized = false;
3516}
3517
3518
3519
3520template <typename MatrixType, typename VectorType, typename PreconditionerType>
3521inline void
3522PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::clear()
3523{
3524 eigenvalues_are_initialized = false;
3525 theta = delta = 1.0;
3526 matrix_ptr = nullptr;
3527 {
3528 VectorType empty_vector;
3529 solution_old.reinit(empty_vector);
3530 temp_vector1.reinit(empty_vector);
3531 temp_vector2.reinit(empty_vector);
3532 }
3533 data.preconditioner.reset();
3534}
3535
3536
3537
3538template <typename MatrixType, typename VectorType, typename PreconditionerType>
3539inline typename PreconditionChebyshev<MatrixType,
3540 VectorType,
3541 PreconditionerType>::EigenvalueInformation
3542PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::
3543 estimate_eigenvalues(const VectorType &src) const
3544{
3545 Assert(eigenvalues_are_initialized == false, ExcInternalError());
3546 Assert(data.preconditioner.get() != nullptr, ExcNotInitialized());
3547
3548 PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::
3549 EigenvalueInformation info{};
3550
3551 solution_old.reinit(src);
3552 temp_vector1.reinit(src, true);
3553
3554 if (data.eig_cg_n_iterations > 0)
3555 {
3556 Assert(data.eig_cg_n_iterations > 2,
3557 ExcMessage(
3558 "Need to set at least two iterations to find eigenvalues."));
3559
3560 internal::PreconditionChebyshevImplementation::EigenvalueTracker
3561 eigenvalue_tracker;
3562
3563 // set an initial guess that contains some high-frequency parts (to the
3564 // extent possible without knowing the discretization and the numbering)
3565 // to trigger high eigenvalues according to the external function
3566 internal::PreconditionChebyshevImplementation::set_initial_guess(
3567 temp_vector1);
3568 data.constraints.set_zero(temp_vector1);
3569
3570 if (data.eigenvalue_algorithm ==
3571 AdditionalData::EigenvalueAlgorithm::lanczos)
3572 {
3573 // set a very strict tolerance to force at least two iterations
3574 IterationNumberControl control(data.eig_cg_n_iterations,
3575 1e-10,
3576 false,
3577 false);
3578
3579 SolverCG<VectorType> solver(control);
3580 solver.connect_eigenvalues_slot(
3581 [&eigenvalue_tracker](const std::vector<double> &eigenvalues) {
3582 eigenvalue_tracker.slot(eigenvalues);
3583 });
3584
3585 solver.solve(*matrix_ptr,
3586 solution_old,
3587 temp_vector1,
3588 *data.preconditioner);
3589
3590 info.cg_iterations = control.last_step();
3591 }
3592 else if (data.eigenvalue_algorithm ==
3593 AdditionalData::EigenvalueAlgorithm::power_iteration)
3594 {
3595 Assert(data.degree != numbers::invalid_unsigned_int,
3596 ExcMessage("Cannot estimate the minimal eigenvalue with the "
3597 "power iteration"));
3598
3599 eigenvalue_tracker.values.push_back(
3600 internal::PreconditionChebyshevImplementation::power_iteration(
3601 *matrix_ptr,
3602 temp_vector1,
3603 *data.preconditioner,
3604 data.eig_cg_n_iterations));
3605 }
3606 else
3607 Assert(false, ExcNotImplemented());
3608
3609 // read the eigenvalues from the attached eigenvalue tracker
3610 if (eigenvalue_tracker.values.empty())
3611 info.min_eigenvalue_estimate = info.max_eigenvalue_estimate = 1.;
3612 else
3613 {
3614 info.min_eigenvalue_estimate = eigenvalue_tracker.values.front();
3615
3616 // include a safety factor since the CG method will in general not
3617 // be converged
3618 info.max_eigenvalue_estimate = 1.2 * eigenvalue_tracker.values.back();
3619 }
3620 }
3621 else
3622 {
3623 info.max_eigenvalue_estimate = data.max_eigenvalue;
3624 info.min_eigenvalue_estimate = data.max_eigenvalue / data.smoothing_range;
3625 }
3626
3627 const double alpha = (data.smoothing_range > 1. ?
3628 info.max_eigenvalue_estimate / data.smoothing_range :
3629 std::min(0.9 * info.max_eigenvalue_estimate,
3630 info.min_eigenvalue_estimate));
3631
3632 // in case the user set the degree to invalid unsigned int, we have to
3633 // determine the number of necessary iterations from the Chebyshev error
3634 // estimate, given the target tolerance specified by smoothing_range. This
3635 // estimate is based on the error formula given in section 5.1 of
3636 // R. S. Varga, Matrix iterative analysis, 2nd ed., Springer, 2009
3637 if (data.degree == numbers::invalid_unsigned_int)
3638 {
3639 const double actual_range = info.max_eigenvalue_estimate / alpha;
3640 const double sigma = (1. - std::sqrt(1. / actual_range)) /
3641 (1. + std::sqrt(1. / actual_range));
3642 const double eps = data.smoothing_range;
3643 const_cast<
3644 PreconditionChebyshev<MatrixType, VectorType, PreconditionerType> *>(
3645 this)
3646 ->data.degree =
3647 1 + static_cast<unsigned int>(
3648 std::log(1. / eps + std::sqrt(1. / eps / eps - 1.)) /
3649 std::log(1. / sigma));
3650 }
3651
3652 info.degree = data.degree;
3653
3654 const_cast<
3655 PreconditionChebyshev<MatrixType, VectorType, PreconditionerType> *>(this)
3656 ->delta =
3657 (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind) ?
3658 (info.max_eigenvalue_estimate) :
3659 ((info.max_eigenvalue_estimate - alpha) * 0.5);
3660 const_cast<
3661 PreconditionChebyshev<MatrixType, VectorType, PreconditionerType> *>(this)
3662 ->theta = (info.max_eigenvalue_estimate + alpha) * 0.5;
3663
3664 // We do not need the second temporary vector in case we have a
3665 // DiagonalMatrix as preconditioner and use deal.II's own vectors
3666 using NumberType = typename VectorType::value_type;
3667 if (std::is_same<PreconditionerType,
3668 ::DiagonalMatrix<VectorType>>::value == false ||
3669 (std::is_same<VectorType, ::Vector<NumberType>>::value == false &&
3670 ((std::is_same<VectorType,
3671 LinearAlgebra::distributed::
3672 Vector<NumberType, MemorySpace::Host>>::value ==
3673 false) ||
3674 (std::is_same<VectorType,
3675 LinearAlgebra::distributed::
3676 Vector<NumberType, MemorySpace::Default>>::value ==
3677 false))))
3678 temp_vector2.reinit(src, true);
3679 else
3680 {
3681 VectorType empty_vector;
3682 temp_vector2.reinit(empty_vector);
3683 }
3684
3685 const_cast<
3686 PreconditionChebyshev<MatrixType, VectorType, PreconditionerType> *>(this)
3687 ->eigenvalues_are_initialized = true;
3688
3689 return info;
3690}
3691
3692
3693
3694template <typename MatrixType, typename VectorType, typename PreconditionerType>
3695inline void
3696PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::vmult(
3697 VectorType & solution,
3698 const VectorType &rhs) const
3699{
3700 std::lock_guard<std::mutex> lock(mutex);
3701 if (eigenvalues_are_initialized == false)
3702 estimate_eigenvalues(rhs);
3703
3704 internal::PreconditionChebyshevImplementation::vmult_and_update(
3705 *matrix_ptr,
3706 *data.preconditioner,
3707 rhs,
3708 0,
3709 0.,
3710 (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind) ?
3711 (4. / (3. * delta)) :
3712 (1. / theta),
3713 solution,
3714 solution_old,
3715 temp_vector1,
3716 temp_vector2);
3717
3718 // if delta is zero, we do not need to iterate because the updates will be
3719 // zero
3720 if (data.degree < 2 || std::abs(delta) < 1e-40)
3721 return;
3722
3723 double rhok = delta / theta, sigma = theta / delta;
3724 for (unsigned int k = 0; k < data.degree - 1; ++k)
3725 {
3726 double factor1 = 0.0;
3727 double factor2 = 0.0;
3728
3729 if (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind)
3730 {
3731 factor1 = (2 * k + 1.) / (2 * k + 5.);
3732 factor2 = (8 * k + 12.) / (delta * (2 * k + 5.));
3733 }
3734 else
3735 {
3736 const double rhokp = 1. / (2. * sigma - rhok);
3737 factor1 = rhokp * rhok;
3738 factor2 = 2. * rhokp / delta;
3739 rhok = rhokp;
3740 }
3741
3742 internal::PreconditionChebyshevImplementation::vmult_and_update(
3743 *matrix_ptr,
3744 *data.preconditioner,
3745 rhs,
3746 k + 1,
3747 factor1,
3748 factor2,
3749 solution,
3750 solution_old,
3751 temp_vector1,
3752 temp_vector2);
3753 }
3754}
3755
3756
3757
3758template <typename MatrixType, typename VectorType, typename PreconditionerType>
3759inline void
3760PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::Tvmult(
3761 VectorType & solution,
3762 const VectorType &rhs) const
3763{
3764 std::lock_guard<std::mutex> lock(mutex);
3765 if (eigenvalues_are_initialized == false)
3766 estimate_eigenvalues(rhs);
3767
3768 internal::PreconditionChebyshevImplementation::vector_updates(
3769 rhs,
3770 *data.preconditioner,
3771 0,
3772 0.,
3773 (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind) ?
3774 (4. / (3. * delta)) :
3775 (1. / theta),
3776 solution_old,
3777 temp_vector1,
3778 temp_vector2,
3779 solution);
3780
3781 if (data.degree < 2 || std::abs(delta) < 1e-40)
3782 return;
3783
3784 double rhok = delta / theta, sigma = theta / delta;
3785 for (unsigned int k = 0; k < data.degree - 1; ++k)
3786 {
3787 double factor1 = 0.0;
3788 double factor2 = 0.0;
3789
3790 if (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind)
3791 {
3792 factor1 = (2 * k + 1.) / (2 * k + 5.);
3793 factor2 = (8 * k + 12.) / (delta * (2 * k + 5.));
3794 }
3795 else
3796 {
3797 const double rhokp = 1. / (2. * sigma - rhok);
3798 factor1 = rhokp * rhok;
3799 factor2 = 2. * rhokp / delta;
3800 rhok = rhokp;
3801 }
3802
3803 matrix_ptr->Tvmult(temp_vector1, solution);
3804 internal::PreconditionChebyshevImplementation::vector_updates(
3805 rhs,
3806 *data.preconditioner,
3807 k + 1,
3808 factor1,
3809 factor2,
3810 solution_old,
3811 temp_vector1,
3812 temp_vector2,
3813 solution);
3814 }
3815}
3816
3817
3818
3819template <typename MatrixType, typename VectorType, typename PreconditionerType>
3820inline void
3821PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::step(
3822 VectorType & solution,
3823 const VectorType &rhs) const
3824{
3825 std::lock_guard<std::mutex> lock(mutex);
3826 if (eigenvalues_are_initialized == false)
3827 estimate_eigenvalues(rhs);
3828
3829 internal::PreconditionChebyshevImplementation::vmult_and_update(
3830 *matrix_ptr,
3831 *data.preconditioner,
3832 rhs,
3833 1,
3834 0.,
3835 (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind) ?
3836 (4. / (3. * delta)) :
3837 (1. / theta),
3838 solution,
3839 solution_old,
3840 temp_vector1,
3841 temp_vector2);
3842
3843 if (data.degree < 2 || std::abs(delta) < 1e-40)
3844 return;
3845
3846 double rhok = delta / theta, sigma = theta / delta;
3847 for (unsigned int k = 0; k < data.degree - 1; ++k)
3848 {
3849 double factor1 = 0.0;
3850 double factor2 = 0.0;
3851
3852 if (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind)
3853 {
3854 factor1 = (2 * k + 1.) / (2 * k + 5.);
3855 factor2 = (8 * k + 12.) / (delta * (2 * k + 5.));
3856 }
3857 else
3858 {
3859 const double rhokp = 1. / (2. * sigma - rhok);
3860 factor1 = rhokp * rhok;
3861 factor2 = 2. * rhokp / delta;
3862 rhok = rhokp;
3863 }
3864
3865 internal::PreconditionChebyshevImplementation::vmult_and_update(
3866 *matrix_ptr,
3867 *data.preconditioner,
3868 rhs,
3869 k + 2,
3870 factor1,
3871 factor2,
3872 solution,
3873 solution_old,
3874 temp_vector1,
3875 temp_vector2);
3876 }
3877}
3878
3879
3880
3881template <typename MatrixType, typename VectorType, typename PreconditionerType>
3882inline void
3883PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::Tstep(
3884 VectorType & solution,
3885 const VectorType &rhs) const
3886{
3887 std::lock_guard<std::mutex> lock(mutex);
3888 if (eigenvalues_are_initialized == false)
3889 estimate_eigenvalues(rhs);
3890
3891 matrix_ptr->Tvmult(temp_vector1, solution);
3892 internal::PreconditionChebyshevImplementation::vector_updates(
3893 rhs,
3894 *data.preconditioner,
3895 1,
3896 0.,
3897 (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind) ?
3898 (4. / (3. * delta)) :
3899 (1. / theta),
3900 solution_old,
3901 temp_vector1,
3902 temp_vector2,
3903 solution);
3904
3905 if (data.degree < 2 || std::abs(delta) < 1e-40)
3906 return;
3907
3908 double rhok = delta / theta, sigma = theta / delta;
3909 for (unsigned int k = 0; k < data.degree - 1; ++k)
3910 {
3911 double factor1 = 0.0;
3912 double factor2 = 0.0;
3913
3914 if (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind)
3915 {
3916 factor1 = (2 * k + 1.) / (2 * k + 5.);
3917 factor2 = (8 * k + 12.) / (delta * (2 * k + 5.));
3918 }
3919 else
3920 {
3921 const double rhokp = 1. / (2. * sigma - rhok);
3922 factor1 = rhokp * rhok;
3923 factor2 = 2. * rhokp / delta;
3924 rhok = rhokp;
3925 }
3926
3927 matrix_ptr->Tvmult(temp_vector1, solution);
3928 internal::PreconditionChebyshevImplementation::vector_updates(
3929 rhs,
3930 *data.preconditioner,
3931 k + 2,
3932 factor1,
3933 factor2,
3934 solution_old,
3935 temp_vector1,
3936 temp_vector2,
3937 solution);
3938 }
3939}
3940
3941
3942
3943template <typename MatrixType, typename VectorType, typename PreconditionerType>
3944inline typename PreconditionChebyshev<MatrixType,
3945 VectorType,
3946 PreconditionerType>::size_type
3947PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::m() const
3948{
3949 Assert(matrix_ptr != nullptr, ExcNotInitialized());
3950 return matrix_ptr->m();
3951}
3952
3953
3954
3955template <typename MatrixType, typename VectorType, typename PreconditionerType>
3956inline typename PreconditionChebyshev<MatrixType,
3957 VectorType,
3958 PreconditionerType>::size_type
3959PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::n() const
3960{
3961 Assert(matrix_ptr != nullptr, ExcNotInitialized());
3962 return matrix_ptr->n();
3963}
3964
3965#endif // DOXYGEN
3966
3967DEAL_II_NAMESPACE_CLOSE
3968
3969#endif
affine_constraints.h
block_vector_base.h
ArrayView::begin
iterator begin() const
Definition array_view.h:594
ArrayView::empty
bool empty() const
Definition array_view.h:585
ArrayView::end
iterator end() const
Definition array_view.h:603
ArrayView::data
value_type * data() const noexcept
Definition array_view.h:553
ArrayView::reinit
void reinit(value_type *starting_element, const std::size_t n_elements)
Definition array_view.h:413
ArrayView::size
std::size_t size() const
Definition array_view.h:576
BlockVectorBase::n_blocks
unsigned int n_blocks() const
BlockVectorBase::block
BlockType & block(const unsigned int i)
LinearAlgebra::distributed::Vector::swap
void swap(Vector< Number, MemorySpace > &v)
LinearAlgebra::distributed::Vector::locally_owned_size
size_type locally_owned_size() const
LinearAlgebra::distributed::Vector::mean_value
virtual Number mean_value() const override
LinearAlgebra::distributed::Vector::sadd
virtual void sadd(const Number s, const Number a, const VectorSpaceVector< Number > &V) override
LinearAlgebra::distributed::Vector::add
virtual void add(const Number a) override
LinearAlgebra::distributed::Vector::locally_owned_elements
virtual ::IndexSet locally_owned_elements() const override
PreconditionChebyshev::Tvmult
void Tvmult(VectorType &dst, const VectorType &src) const
PreconditionChebyshev::m
size_type m() const
PreconditionChebyshev::n
size_type n() const
PreconditionChebyshev::step
void step(VectorType &dst, const VectorType &src) const
PreconditionChebyshev::estimate_eigenvalues
EigenvalueInformation estimate_eigenvalues(const VectorType &src) const
PreconditionChebyshev::Tstep
void Tstep(VectorType &dst, const VectorType &src) const
PreconditionChebyshev::matrix_ptr
SmartPointer< const MatrixType, PreconditionChebyshev< MatrixType, VectorType, PreconditionerType > > matrix_ptr
Definition precondition.h:2197
PreconditionChebyshev::vmult
void vmult(VectorType &dst, const VectorType &src) const
PreconditionChebyshev::initialize
void initialize(const MatrixType &matrix, const AdditionalData &additional_data=AdditionalData())
PreconditionIdentity::vmult_add
void vmult_add(VectorType &, const VectorType &) const
PreconditionIdentity::vmult
void vmult(VectorType &, const VectorType &) const
PreconditionIdentity::m
size_type m() const
PreconditionIdentity::initialize
void initialize(const MatrixType &matrix, const AdditionalData &additional_data=AdditionalData())
PreconditionIdentity::n
size_type n() const
PreconditionIdentity::Tvmult
void Tvmult(VectorType &, const VectorType &) const
PreconditionIdentity::Tvmult_add
void Tvmult_add(VectorType &, const VectorType &) const
PreconditionJacobi::AdditionalData
typename BaseClass::AdditionalData AdditionalData
Definition precondition.h:1483
PreconditionJacobi::PreconditionerType
internal::PreconditionRelaxation::PreconditionJacobiImpl< MatrixType > PreconditionerType
Definition precondition.h:1476
PreconditionJacobi::initialize
void initialize(const MatrixType &A, const AdditionalData &parameters=AdditionalData())
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())
PreconditionPSOR::AdditionalData::parameters
BaseClass::AdditionalData parameters
Definition precondition.h:1695
PreconditionPSOR::AdditionalData::inverse_permutation
const std::vector< size_type > & inverse_permutation
Definition precondition.h:1691
PreconditionPSOR::AdditionalData::permutation
const std::vector< size_type > & permutation
Definition precondition.h:1687
PreconditionPSOR::size_type
typename BaseClass::size_type size_type
Definition precondition.h:1661
PreconditionPSOR::initialize
void initialize(const MatrixType &A, const AdditionalData &additional_data)
PreconditionPSOR::PreconditionerType
internal::PreconditionRelaxation::PreconditionPSORImpl< MatrixType > PreconditionerType
Definition precondition.h:1654
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())
PreconditionRelaxation::AdditionalData::preconditioner
std::shared_ptr< PreconditionerType > preconditioner
Definition precondition.h:440
PreconditionRelaxation::AdditionalData::AdditionalData
AdditionalData(const double relaxation=1., const unsigned int n_iterations=1)
PreconditionRelaxation::A
SmartPointer< const MatrixType, PreconditionRelaxation< MatrixType > > A
Definition precondition.h:505
PreconditionRelaxation::Tvmult
void Tvmult(VectorType &, const VectorType &) const
PreconditionRelaxation::preconditioner
std::shared_ptr< PreconditionerType > preconditioner
Definition precondition.h:520
PreconditionRelaxation::n
size_type n() const
PreconditionRelaxation::m
size_type m() const
PreconditionRelaxation::step
void step(VectorType &x, const VectorType &rhs) const
PreconditionRelaxation::n_iterations
unsigned int n_iterations
Definition precondition.h:515
PreconditionRelaxation::initialize
void initialize(const MatrixType &A, const AdditionalData &parameters=AdditionalData())
PreconditionRelaxation::Tstep
void Tstep(VectorType &x, const VectorType &rhs) const
PreconditionRelaxation::vmult
void vmult(VectorType &, const VectorType &) const
PreconditionRichardson::AdditionalData::AdditionalData
AdditionalData(const double relaxation=1.)
PreconditionRichardson::size_type
types::global_dof_index size_type
Definition precondition.h:204
PreconditionRichardson::vmult_add
void vmult_add(VectorType &, const VectorType &) const
PreconditionRichardson::n
size_type n() const
PreconditionRichardson::initialize
void initialize(const AdditionalData &parameters)
PreconditionRichardson::initialize
void initialize(const MatrixType &matrix, const AdditionalData &parameters)
PreconditionRichardson::m
size_type m() const
PreconditionRichardson::vmult
void vmult(VectorType &, const VectorType &) const
PreconditionRichardson::Tvmult
void Tvmult(VectorType &, const VectorType &) const
PreconditionRichardson::Tvmult_add
void Tvmult_add(VectorType &, const VectorType &) const
PreconditionSOR::PreconditionerType
internal::PreconditionRelaxation::PreconditionSORImpl< MatrixType > PreconditionerType
Definition precondition.h:1546
PreconditionSOR::AdditionalData
typename BaseClass::AdditionalData AdditionalData
Definition precondition.h:1553
PreconditionSOR::initialize
void initialize(const MatrixType &A, const AdditionalData &parameters=AdditionalData())
PreconditionSSOR::initialize
void initialize(const MatrixType &A, const AdditionalData &parameters=AdditionalData())
PreconditionSSOR::AdditionalData
typename BaseClass::AdditionalData AdditionalData
Definition precondition.h:1605
PreconditionSSOR::PreconditionerType
internal::PreconditionRelaxation::PreconditionSSORImpl< MatrixType > PreconditionerType
Definition precondition.h:1598
PreconditionUseMatrix::function_ptr
void(MatrixType::*)(VectorType &, const VectorType &) const function_ptr
Definition precondition.h:369
PreconditionUseMatrix::precondition
const function_ptr precondition
Definition precondition.h:394
PreconditionUseMatrix::matrix
const MatrixType & matrix
Definition precondition.h:389
PreconditionUseMatrix::vmult
void vmult(VectorType &dst, const VectorType &src) const
PreconditionUseMatrix::PreconditionUseMatrix
PreconditionUseMatrix(const MatrixType &M, const function_ptr method)
DEAL_II_OPENMP_SIMD_PRAGMA
#define DEAL_II_OPENMP_SIMD_PRAGMA
Definition config.h:138
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
cuda_size.h
diagonal_matrix.h
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
Assert
#define Assert(cond, exc)
Definition exceptions.h:1614
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
StandardExceptions::ExcNotInitialized
static ::ExceptionBase & ExcNotInitialized()
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
AssertThrow
#define AssertThrow(cond, exc)
Definition exceptions.h:1703
identity_matrix.h
vector_memory.h
memory_space.h
LAPACKSupport::matrix
@ matrix
Contents is actually a matrix.
Definition lapack_support.h:58
TrilinosWrappers::internal::end
VectorType::value_type * end(VectorType &V)
Definition trilinos_sparse_matrix.cc:65
TrilinosWrappers::internal::begin
VectorType::value_type * begin(VectorType &V)
Definition trilinos_sparse_matrix.cc:51
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:213
std
STL namespace.
std::log
::VectorizedArray< Number, width > log(const ::VectorizedArray< Number, width > &)
Definition vectorization.h:5924
std::min
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6041
std::sqrt
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
Definition vectorization.h:5946
std::abs
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6008
types::global_dof_index
unsigned int global_dof_index
Definition types.h:82
solver_cg.h
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, const PolynomialType polynomial_type=PolynomialType::first_kind)
PreconditionChebyshev::AdditionalData::preconditioner
std::shared_ptr< PreconditionerType > preconditioner
Definition precondition.h:2061
PreconditionChebyshev::AdditionalData::operator=
AdditionalData & operator=(const AdditionalData &other_data)
PreconditionChebyshev::AdditionalData::constraints
AffineConstraints< double > constraints
Definition precondition.h:2056
eigenvalues
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)
SymmetricTensorEigenvectorMethod::jacobi
@ jacobi
template_constraints.h
Generated by doxygen 1.9.8