Loading [MathJax]/extensions/TeX/AMSsymbols.js
 deal.II version GIT relicensing-3067-g24f3489fdc 2025-04-14 21:10:00+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages Concepts
solver_cg.h
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 1998 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15#ifndef dealii_solver_cg_h
16#define dealii_solver_cg_h
17
18
19#include <deal.II/base/config.h>
20
26
27#include <deal.II/lac/solver.h>
30
31#include <boost/signals2.hpp>
32
33#include <cmath>
34
36
37// forward declaration
38#ifndef DOXYGEN
40namespace LinearAlgebra
41{
42 namespace distributed
43 {
44 template <typename, typename>
45 class Vector;
46 }
47} // namespace LinearAlgebra
48#endif
49
50
196template <typename VectorType = Vector<double>>
198class SolverCG : public SolverBase<VectorType>
199{
200public:
205
206
211 {
226 explicit AdditionalData(const bool use_default_residual = true);
227
232 };
233
234
241
247
251 virtual ~SolverCG() override = default;
252
256 template <typename MatrixType, typename PreconditionerType>
260 void solve(const MatrixType &A,
261 VectorType &x,
262 const VectorType &b,
263 const PreconditionerType &preconditioner);
264
271 boost::signals2::connection
272 connect_coefficients_slot(
273 const std::function<void(typename VectorType::value_type,
274 typename VectorType::value_type)> &slot);
275
282 boost::signals2::connection
283 connect_condition_number_slot(const std::function<void(double)> &slot,
284 const bool every_iteration = false);
285
292 boost::signals2::connection
293 connect_eigenvalues_slot(
294 const std::function<void(const std::vector<double> &)> &slot,
295 const bool every_iteration = false);
296
297protected:
303 virtual void
304 print_vectors(const unsigned int step,
305 const VectorType &x,
306 const VectorType &r,
307 const VectorType &d) const;
308
314 static void
315 compute_eigs_and_cond(
316 const std::vector<typename VectorType::value_type> &diagonal,
317 const std::vector<typename VectorType::value_type> &offdiagonal,
318 const boost::signals2::signal<void(const std::vector<double> &)>
319 &eigenvalues_signal,
320 const boost::signals2::signal<void(double)> &cond_signal);
321
325 AdditionalData additional_data;
326
330 boost::signals2::signal<void(typename VectorType::value_type,
331 typename VectorType::value_type)>
332 coefficients_signal;
333
338 boost::signals2::signal<void(double)> condition_number_signal;
339
344 boost::signals2::signal<void(double)> all_condition_numbers_signal;
345
350 boost::signals2::signal<void(const std::vector<double> &)> eigenvalues_signal;
351
356 boost::signals2::signal<void(const std::vector<double> &)>
357 all_eigenvalues_signal;
358
367 bool determine_beta_by_flexible_formula;
368};
369
370
371
397template <typename VectorType = Vector<double>>
398DEAL_II_CXX20_REQUIRES(concepts::is_vector_space_vector<VectorType>)
399class SolverFlexibleCG : public SolverCG<VectorType>
400{
401public:
406
407
412 {
427 explicit AdditionalData(const bool use_default_residual = true);
428
433 };
434
435
442
449
450protected:
455};
456
457
460/*------------------------- Implementation ----------------------------*/
461
462#ifndef DOXYGEN
463
464
465template <typename VectorType>
468 const bool use_default_residual)
469 : use_default_residual(use_default_residual)
470{}
471
472template <typename VectorType>
476 const AdditionalData &data)
477 : SolverBase<VectorType>(cn, mem)
478 , additional_data(data)
479 , determine_beta_by_flexible_formula(false)
480{}
481
482
483
484template <typename VectorType>
486SolverCG<VectorType>::SolverCG(SolverControl &cn, const AdditionalData &data)
488 , additional_data(data)
489 , determine_beta_by_flexible_formula(false)
490{}
491
492
493
494template <typename VectorType>
496void SolverCG<VectorType>::print_vectors(const unsigned int,
497 const VectorType &,
498 const VectorType &,
499 const VectorType &) const
500{}
501
502
503
504template <typename VectorType>
507 const std::vector<typename VectorType::value_type> &diagonal,
508 const std::vector<typename VectorType::value_type> &offdiagonal,
509 const boost::signals2::signal<void(const std::vector<double> &)>
510 &eigenvalues_signal,
511 const boost::signals2::signal<void(double)> &cond_signal)
512{
513 // Avoid computing eigenvalues unless they are needed.
514 if (!cond_signal.empty() || !eigenvalues_signal.empty())
515 {
517 true);
518 for (size_type i = 0; i < diagonal.size(); ++i)
519 {
520 T(i, i) = diagonal[i];
521 if (i < diagonal.size() - 1)
522 T(i, i + 1) = offdiagonal[i];
523 }
524 T.compute_eigenvalues();
525 // Need two eigenvalues to estimate the condition number.
526 if (diagonal.size() > 1)
527 {
528 auto condition_number = T.eigenvalue(T.n() - 1) / T.eigenvalue(0);
529 // Condition number is real valued and nonnegative; simply take
530 // the absolute value:
531 cond_signal(std::abs(condition_number));
532 }
533 // Avoid copying the eigenvalues of T to a vector unless a signal is
534 // connected.
535 if (!eigenvalues_signal.empty())
536 {
537 std::vector<double> eigenvalues(T.n());
538 for (unsigned int j = 0; j < T.n(); ++j)
539 {
540 // for a hermitian matrix, all eigenvalues are real-valued
541 // and non-negative, simply return the absolute value:
542 eigenvalues[j] = std::abs(T.eigenvalue(j));
543 }
544 eigenvalues_signal(eigenvalues);
545 }
546 }
547}
548
549
550
551namespace internal
552{
553
554 namespace SolverCG
555 {
556 // This base class is used to select different variants of the conjugate
557 // gradient solver. The default variant is used for standard matrix and
558 // preconditioner arguments, as provided by the derived class
559 // IterationWork below, but there is also a specialized variant further
560 // down that uses SFINAE to identify whether matrices and preconditioners
561 // support special operations on sub-ranges of the vectors.
562 template <typename VectorType,
563 typename MatrixType,
564 typename PreconditionerType>
565 struct IterationWorkerBase
566 {
567 using Number = typename VectorType::value_type;
568
569 const MatrixType &A;
570 const PreconditionerType &preconditioner;
571 const bool flexible;
572 VectorType &x;
573 const VectorType &b;
574
575 typename VectorMemory<VectorType>::Pointer r_pointer;
576 typename VectorMemory<VectorType>::Pointer p_pointer;
577 typename VectorMemory<VectorType>::Pointer v_pointer;
578 typename VectorMemory<VectorType>::Pointer z_pointer;
579 typename VectorMemory<VectorType>::Pointer explicit_r_pointer;
580
581 // Define some aliases for simpler access, using the variables 'r' for
582 // the residual b - A*x, 'p' for the search direction, and 'v' for the
583 // auxiliary vector. This naming convention is used e.g. by the
584 // description on
585 // https://en.wikipedia.org/wiki/Conjugate_gradient_method. The variable
586 // 'z' gets only used for the flexible variant of the CG method.
587 VectorType &r;
588 VectorType &p;
589 VectorType &v;
590 VectorType &z;
591 VectorType &explicit_r;
592
593 Number r_dot_preconditioner_dot_r;
594 Number alpha;
595 Number beta;
596 double residual_norm;
597 Number previous_alpha;
598 const bool use_default_residual;
599
600
601 IterationWorkerBase(const MatrixType &A,
602 const PreconditionerType &preconditioner,
603 const bool flexible,
605 VectorType &x,
606 const VectorType &b,
607 const bool &use_default_residual)
608 : A(A)
609 , preconditioner(preconditioner)
610 , flexible(flexible)
611 , x(x)
612 , b(b)
613 , r_pointer(memory)
614 , p_pointer(memory)
615 , v_pointer(memory)
616 , z_pointer(memory)
617 , explicit_r_pointer(memory)
618 , r(*r_pointer)
619 , p(*p_pointer)
620 , v(*v_pointer)
621 , z(*z_pointer)
622 , explicit_r(*explicit_r_pointer)
623 , r_dot_preconditioner_dot_r(Number())
624 , alpha(Number())
625 , beta(Number())
626 , residual_norm(0.0)
627 , previous_alpha(Number())
628 , use_default_residual(use_default_residual)
629 {}
630
631 void
632 startup()
633 {
634 // Initialize without setting the vector entries, as those would soon
635 // be overwritten anyway
636 r.reinit(x, true);
637 p.reinit(x, true);
638 v.reinit(x, true);
639 if (flexible)
640 z.reinit(x, true);
641 if (!use_default_residual)
642 explicit_r.reinit(x, true);
643
644 // compute residual. if vector is zero, then short-circuit the full
645 // computation
646 if (!x.all_zero())
647 {
648 A.vmult(r, x);
649 r.sadd(-1., 1., b);
650 }
651 else
652 r.equ(1., b);
653
654 residual_norm = r.l2_norm();
655 }
656 };
657
658
659
660 // Implementation of a conjugate gradient operation with matrices and
661 // preconditioners without special capabilities
662 template <typename VectorType,
663 typename MatrixType,
664 typename PreconditionerType,
665 typename = int>
666 struct IterationWorker
667 : public IterationWorkerBase<VectorType, MatrixType, PreconditionerType>
668 {
669 using BaseClass =
670 IterationWorkerBase<VectorType, MatrixType, PreconditionerType>;
671
672
673 IterationWorker(const MatrixType &A,
674 const PreconditionerType &preconditioner,
675 const bool flexible,
677 VectorType &x,
678 const VectorType &b,
679 const bool &use_default_residual)
680 : BaseClass(A,
681 preconditioner,
682 flexible,
683 memory,
684 x,
685 b,
686 use_default_residual)
687 {}
688
689 using BaseClass::A;
690 using BaseClass::alpha;
691 using BaseClass::b;
692 using BaseClass::beta;
693 using BaseClass::explicit_r;
694 using BaseClass::p;
695 using BaseClass::preconditioner;
696 using BaseClass::r;
697 using BaseClass::r_dot_preconditioner_dot_r;
698 using BaseClass::residual_norm;
699 using BaseClass::use_default_residual;
700 using BaseClass::v;
701 using BaseClass::x;
702 using BaseClass::z;
703
704 void
705 do_iteration(const unsigned int iteration_index)
706 {
707 using Number = typename VectorType::value_type;
708
709 const Number previous_r_dot_preconditioner_dot_r =
710 r_dot_preconditioner_dot_r;
711
712 if (std::is_same_v<PreconditionerType, PreconditionIdentity> == false)
713 {
714 preconditioner.vmult(v, r);
715 r_dot_preconditioner_dot_r = r * v;
716 }
717 else
718 r_dot_preconditioner_dot_r = residual_norm * residual_norm;
719
720 const VectorType &direction =
721 std::is_same_v<PreconditionerType, PreconditionIdentity> ? r : v;
722
723 if (iteration_index > 1)
724 {
725 Assert(std::abs(previous_r_dot_preconditioner_dot_r) != 0.,
727 beta =
728 r_dot_preconditioner_dot_r / previous_r_dot_preconditioner_dot_r;
729 if (this->flexible)
730 beta -= (r * z) / previous_r_dot_preconditioner_dot_r;
731 p.sadd(beta, 1., direction);
732 }
733 else
734 p.equ(1., direction);
735
736 if (this->flexible)
737 z.swap(v);
738
739 A.vmult(v, p);
740
741 const Number p_dot_A_dot_p = p * v;
742 Assert(std::abs(p_dot_A_dot_p) != 0., ExcDivideByZero());
743
744 this->previous_alpha = alpha;
745 alpha = r_dot_preconditioner_dot_r / p_dot_A_dot_p;
746
747 x.add(alpha, p);
748
749 // compute the residual norm with implicit residual
750 if (use_default_residual)
751 {
752 residual_norm = std::sqrt(std::abs(r.add_and_dot(-alpha, v, r)));
753 }
754 // compute the residual norm with the explicit residual, i.e.
755 // compute l2 norm of Ax - b.
756 else
757 {
758 // compute the residual conjugate gradient update
759 r.add(-alpha, v);
760 // compute explicit residual
761 A.vmult(explicit_r, x);
762 explicit_r.add(-1, b);
763 residual_norm = explicit_r.l2_norm();
764 }
765 }
766
767 void
768 finalize_after_convergence(const unsigned int)
769 {}
770 };
771
772
773 // In the following, we provide a specialization of the above
774 // IterationWorker class that picks up particular features in the matrix
775 // and preconditioners.
776
777 // a helper type-trait that leverage SFINAE to figure out if MatrixType has
778 // ... MatrixType::vmult(VectorType &, const VectorType&,
779 // std::function<...>, std::function<...>) const
780 template <typename MatrixType, typename VectorType>
781 using vmult_functions_t = decltype(std::declval<const MatrixType>().vmult(
782 std::declval<VectorType &>(),
783 std::declval<const VectorType &>(),
784 std::declval<
785 const std::function<void(const unsigned int, const unsigned int)> &>(),
786 std::declval<const std::function<void(const unsigned int,
787 const unsigned int)> &>()));
788
789 template <typename MatrixType, typename VectorType>
790 constexpr bool has_vmult_functions =
791 is_supported_operation<vmult_functions_t, MatrixType, VectorType>;
792
793 // a helper type-trait that leverage SFINAE to figure out if
794 // PreconditionerType has ... PreconditionerType::apply_to_subrange(const
795 // unsigned int, const unsigned int, const Number*, Number*) const
796 template <typename PreconditionerType>
797 using apply_to_subrange_t =
798 decltype(std::declval<const PreconditionerType>()
799 .apply_to_subrange(0U, 0U, nullptr, nullptr));
800
801 template <typename PreconditionerType>
802 constexpr bool has_apply_to_subrange =
803 is_supported_operation<apply_to_subrange_t, PreconditionerType>;
804
805 // a helper type-trait that leverage SFINAE to figure out if
806 // PreconditionerType has ... PreconditionerType::apply(const
807 // unsigned int, const Number) const
808 template <typename PreconditionerType>
809 using apply_t =
810 decltype(std::declval<const PreconditionerType>().apply(0U, 0.0));
811
812 template <typename PreconditionerType>
813 constexpr bool has_apply =
814 is_supported_operation<apply_t, PreconditionerType>;
815
816
817 // Internal function to run one iteration of the conjugate gradient solver
818 // for matrices and preconditioners that support interleaving the vector
819 // updates with the matrix-vector product.
820 template <typename VectorType,
821 typename MatrixType,
822 typename PreconditionerType>
823 struct IterationWorker<
826 PreconditionerType,
827 std::enable_if_t<has_vmult_functions<MatrixType, VectorType> &&
828 (has_apply_to_subrange<PreconditionerType> ||
829 has_apply<PreconditionerType>)&&std::
830 is_same_v<VectorType,
831 LinearAlgebra::distributed::Vector<
832 typename VectorType::value_type,
833 MemorySpace::Host>>,
834 int>>
835 : public IterationWorkerBase<VectorType, MatrixType, PreconditionerType>
836 {
837 using Number = typename VectorType::value_type;
838
839 Number next_r_dot_preconditioner_dot_r;
840 Number previous_beta;
841
842 IterationWorker(const MatrixType &A,
843 const PreconditionerType &preconditioner,
844 const bool flexible,
846 VectorType &x,
847 const VectorType &b,
848 const bool &use_default_residual)
849 : IterationWorkerBase<VectorType, MatrixType, PreconditionerType>(
850 A,
851 preconditioner,
852 flexible,
853 memory,
854 x,
855 b,
856 use_default_residual)
857 , next_r_dot_preconditioner_dot_r(0.)
858 , previous_beta(0.)
859 {}
860
861 // This is the main iteration function, that will use some of the
862 // specialized functions below
863 void
864 do_iteration(const unsigned int iteration_index)
865 {
866 if (iteration_index > 1)
867 {
868 previous_beta = this->beta;
869 this->beta = next_r_dot_preconditioner_dot_r /
870 this->r_dot_preconditioner_dot_r;
871 }
872
873 std::array<VectorizedArray<Number>, 7> vectorized_sums = {};
874
875 this->A.vmult(
876 this->v,
877 this->p,
878 [&](const unsigned int begin, const unsigned int end) {
879 operation_before_loop(iteration_index, begin, end);
880 },
881 [&](const unsigned int begin, const unsigned int end) {
882 operation_after_loop(begin, end, vectorized_sums);
883 });
884
885 std::array<Number, 7> scalar_sums;
886 for (unsigned int i = 0; i < 7; ++i)
887 scalar_sums[i] = vectorized_sums[i][0];
888 for (unsigned int l = 1; l < VectorizedArray<Number>::size(); ++l)
889 for (unsigned int i = 0; i < 7; ++i)
890 scalar_sums[i] += vectorized_sums[i][l];
891
893 7),
894 this->r.get_mpi_communicator(),
895 ::ArrayView<Number>(scalar_sums.data(), 7));
896
897 this->r_dot_preconditioner_dot_r = scalar_sums[6];
898
899 const Number p_dot_A_dot_p = scalar_sums[0];
900 Assert(std::abs(p_dot_A_dot_p) != 0., ExcDivideByZero());
901
902 this->previous_alpha = this->alpha;
903 this->alpha = this->r_dot_preconditioner_dot_r / p_dot_A_dot_p;
904
905 // Round-off errors near zero might yield negative values, so take
906 // the absolute value in the next two formulas
907 this->residual_norm = std::sqrt(std::abs(
908 scalar_sums[3] +
909 this->alpha * (-2. * scalar_sums[2] + this->alpha * scalar_sums[1])));
910
911 next_r_dot_preconditioner_dot_r = std::abs(
912 this->r_dot_preconditioner_dot_r +
913 this->alpha * (-2. * scalar_sums[4] + this->alpha * scalar_sums[5]));
914 }
915
916 // Function that we use if the PreconditionerType implements an apply()
917 // function
918 template <typename U = void>
919 std::enable_if_t<has_apply<PreconditionerType>, U>
920 operation_before_loop(const unsigned int iteration_index,
921 const unsigned int start_range,
922 const unsigned int end_range) const
923 {
924 Number *x = this->x.begin();
925 Number *r = this->r.begin();
926 Number *p = this->p.begin();
927 Number *v = this->v.begin();
928 const Number alpha = this->alpha;
929 const Number beta = this->beta;
930 constexpr unsigned int n_lanes = VectorizedArray<Number>::size();
931 const unsigned int end_regular =
932 start_range + (end_range - start_range) / n_lanes * n_lanes;
933 if (iteration_index == 1)
934 {
935 // Vectorize by hand since compilers are often pretty bad at
936 // doing these steps automatically even with
937 // DEAL_II_OPENMP_SIMD_PRAGMA
938 for (unsigned int j = start_range; j < end_regular; j += n_lanes)
939 {
941 rj.load(r + j);
943 for (unsigned int l = 0; l < n_lanes; ++l)
944 pj[l] = this->preconditioner.apply(j + l, rj[l]);
945 pj.store(p + j);
947 rj.store(v + j);
948 }
949 for (unsigned int j = end_regular; j < end_range; ++j)
950 {
951 p[j] = this->preconditioner.apply(j, r[j]);
952 v[j] = Number();
953 }
954 }
955 else if (iteration_index % 2 == 0 && beta != Number())
956 {
957 for (unsigned int j = start_range; j < end_regular; j += n_lanes)
958 {
959 VectorizedArray<Number> rj, vj, pj, prec_rj;
960 rj.load(r + j);
961 vj.load(v + j);
962 rj -= alpha * vj;
963 rj.store(r + j);
965 for (unsigned int l = 0; l < n_lanes; ++l)
966 prec_rj[l] = this->preconditioner.apply(j + l, rj[l]);
967 pj.load(p + j);
968 pj = beta * pj + prec_rj;
969 pj.store(p + j);
971 rj.store(v + j);
972 }
973 for (unsigned int j = end_regular; j < end_range; ++j)
974 {
975 r[j] -= alpha * v[j];
976 p[j] = beta * p[j] + this->preconditioner.apply(j, r[j]);
977 v[j] = Number();
978 }
979 }
980 else if (iteration_index % 2 == 0 && beta == Number())
981 {
982 // Case where beta is zero: we cannot reconstruct p_{j-1} in the
983 // next iteration, and must hence compute x here. This can happen
984 // before termination
985 for (unsigned int j = start_range; j < end_regular; j += n_lanes)
986 {
987 VectorizedArray<Number> rj, xj, vj, pj, prec_rj;
988 rj.load(r + j);
989 vj.load(v + j);
990 xj.load(x + j);
991 pj.load(p + j);
992 xj += alpha * pj;
993 xj.store(x + j);
994 rj -= alpha * vj;
995 rj.store(r + j);
997 for (unsigned int l = 0; l < n_lanes; ++l)
998 prec_rj[l] = this->preconditioner.apply(j + l, rj[l]);
999 prec_rj.store(p + j);
1001 rj.store(v + j);
1002 }
1003 for (unsigned int j = end_regular; j < end_range; ++j)
1004 {
1005 r[j] -= alpha * v[j];
1006 x[j] += alpha * p[j];
1007 p[j] = this->preconditioner.apply(j, r[j]);
1008 v[j] = Number();
1009 }
1010 }
1011 else
1012 {
1013 const Number alpha_plus_previous_alpha_over_beta =
1014 alpha + this->previous_alpha / this->previous_beta;
1015 const Number previous_alpha_over_beta =
1016 this->previous_alpha / this->previous_beta;
1017 for (unsigned int j = start_range; j < end_regular; j += n_lanes)
1018 {
1019 VectorizedArray<Number> rj, vj, pj, xj, prec_rj, prec_vj;
1020 xj.load(x + j);
1021 pj.load(p + j);
1022 xj += alpha_plus_previous_alpha_over_beta * pj;
1023 rj.load(r + j);
1024 vj.load(v + j);
1026 for (unsigned int l = 0; l < n_lanes; ++l)
1027 {
1028 prec_rj[l] = this->preconditioner.apply(j + l, rj[l]);
1029 prec_vj[l] = this->preconditioner.apply(j + l, vj[l]);
1030 }
1031 xj -= previous_alpha_over_beta * prec_rj;
1032 xj.store(x + j);
1033 rj -= alpha * vj;
1034 rj.store(r + j);
1035 prec_rj -= alpha * prec_vj;
1036 pj = beta * pj + prec_rj;
1037 pj.store(p + j);
1039 rj.store(v + j);
1040 }
1041 for (unsigned int j = end_regular; j < end_range; ++j)
1042 {
1043 x[j] += alpha_plus_previous_alpha_over_beta * p[j];
1044 x[j] -= previous_alpha_over_beta *
1045 this->preconditioner.apply(j, r[j]);
1046 r[j] -= alpha * v[j];
1047 p[j] = beta * p[j] + this->preconditioner.apply(j, r[j]);
1048 v[j] = Number();
1049 }
1050 }
1051 }
1052
1053 // Function that we use if the PreconditionerType implements an apply()
1054 // function
1055 template <typename U = void>
1056 std::enable_if_t<has_apply<PreconditionerType>, U>
1057 operation_after_loop(
1058 const unsigned int start_range,
1059 const unsigned int end_range,
1060 std::array<VectorizedArray<Number>, 7> &vectorized_sums) const
1061 {
1062 const Number *r = this->r.begin();
1063 const Number *p = this->p.begin();
1064 const Number *v = this->v.begin();
1065 std::array<VectorizedArray<Number>, 7> my_sums = {};
1066 constexpr unsigned int n_lanes = VectorizedArray<Number>::size();
1067 const unsigned int end_regular =
1068 start_range + (end_range - start_range) / n_lanes * n_lanes;
1069 for (unsigned int j = start_range; j < end_regular; j += n_lanes)
1070 {
1071 VectorizedArray<Number> pj, vj, rj, prec_vj, prec_rj;
1072 pj.load(p + j);
1073 vj.load(v + j);
1074 rj.load(r + j);
1076 for (unsigned int l = 0; l < n_lanes; ++l)
1077 {
1078 prec_vj[l] = this->preconditioner.apply(j + l, vj[l]);
1079 prec_rj[l] = this->preconditioner.apply(j + l, rj[l]);
1080 }
1081 my_sums[0] += pj * vj;
1082 my_sums[1] += vj * vj;
1083 my_sums[2] += rj * vj;
1084 my_sums[3] += rj * rj;
1085 my_sums[4] += rj * prec_vj;
1086 my_sums[5] += vj * prec_vj;
1087 my_sums[6] += rj * prec_rj;
1088 }
1089 for (unsigned int j = end_regular; j < end_range; ++j)
1090 {
1091 const Number prec_v = this->preconditioner.apply(j, v[j]);
1092 const Number prec_r = this->preconditioner.apply(j, r[j]);
1093 my_sums[0][0] += p[j] * v[j];
1094 my_sums[1][0] += v[j] * v[j];
1095 my_sums[2][0] += r[j] * v[j];
1096 my_sums[3][0] += r[j] * r[j];
1097 my_sums[4][0] += r[j] * prec_v;
1098 my_sums[5][0] += v[j] * prec_v;
1099 my_sums[6][0] += r[j] * prec_r;
1100 }
1101 for (unsigned int i = 0; i < vectorized_sums.size(); ++i)
1102 vectorized_sums[i] += my_sums[i];
1103 }
1104
1105 // Function that we use if the PreconditionerType implements an apply()
1106 // function
1107 template <typename U = void>
1108 std::enable_if_t<has_apply<PreconditionerType>, U>
1109 finalize_after_convergence(const unsigned int iteration_index)
1110 {
1111 if (iteration_index % 2 == 1 || this->beta == Number())
1112 this->x.add(this->alpha, this->p);
1113 else
1114 {
1115 using Number = typename VectorType::value_type;
1116 const unsigned int end_range = this->x.locally_owned_size();
1117
1118 Number *x = this->x.begin();
1119 Number *r = this->r.begin();
1120 Number *p = this->p.begin();
1121
1122 // Note that we use 'beta' here rather than 'previous_beta' in the
1123 // formula above, which is because the shift in beta ->
1124 // previous_beta has not been applied at this stage, allowing us
1125 // to recover the previous search direction
1126 const Number alpha_plus_previous_alpha_over_beta =
1127 this->alpha + this->previous_alpha / this->beta;
1128 const Number previous_alpha_over_beta =
1129 this->previous_alpha / this->beta;
1130
1132 for (unsigned int j = 0; j < end_range; ++j)
1133 {
1134 x[j] += alpha_plus_previous_alpha_over_beta * p[j] -
1135 previous_alpha_over_beta *
1136 this->preconditioner.apply(j, r[j]);
1137 }
1138 }
1139 }
1140
1141 // Function that we use if the PreconditionerType does not implement an
1142 // apply() function, where we instead need to choose the
1143 // apply_to_subrange function
1144 template <typename U = void>
1145 std::enable_if_t<!has_apply<PreconditionerType>, U>
1146 operation_before_loop(const unsigned int iteration_index,
1147 const unsigned int start_range,
1148 const unsigned int end_range) const
1149 {
1150 Number *x = this->x.begin() + start_range;
1151 Number *r = this->r.begin() + start_range;
1152 Number *p = this->p.begin() + start_range;
1153 Number *v = this->v.begin() + start_range;
1154 const Number alpha = this->alpha;
1155 const Number beta = this->beta;
1156 constexpr unsigned int grain_size = 128;
1157 std::array<Number, grain_size> prec_r;
1158 if (iteration_index == 1)
1159 {
1160 for (unsigned int j = start_range; j < end_range; j += grain_size)
1161 {
1162 const unsigned int length = std::min(grain_size, end_range - j);
1163 this->preconditioner.apply_to_subrange(j,
1164 j + length,
1165 r,
1166 prec_r.data());
1168 for (unsigned int i = 0; i < length; ++i)
1169 {
1170 p[i] = prec_r[i];
1171 v[i] = Number();
1172 }
1173 p += length;
1174 r += length;
1175 v += length;
1176 }
1177 }
1178 else if (iteration_index % 2 == 0 && beta != Number())
1179 {
1180 for (unsigned int j = start_range; j < end_range; j += grain_size)
1181 {
1182 const unsigned int length = std::min(grain_size, end_range - j);
1184 for (unsigned int i = 0; i < length; ++i)
1185 r[i] -= this->alpha * v[i];
1186 this->preconditioner.apply_to_subrange(j,
1187 j + length,
1188 r,
1189 prec_r.data());
1191 for (unsigned int i = 0; i < length; ++i)
1192 {
1193 p[i] = this->beta * p[i] + prec_r[i];
1194 v[i] = Number();
1195 }
1196 p += length;
1197 r += length;
1198 v += length;
1199 }
1200 }
1201 else if (iteration_index % 2 == 0 && beta == Number())
1202 {
1203 // Case where beta is zero: We cannot reconstruct p_{j-1} in the
1204 // next iteration, and must hence compute x here. This can happen
1205 // before termination
1206 for (unsigned int j = start_range; j < end_range; j += grain_size)
1207 {
1208 const unsigned int length = std::min(grain_size, end_range - j);
1210 for (unsigned int i = 0; i < length; ++i)
1211 r[i] -= this->alpha * v[i];
1212 this->preconditioner.apply_to_subrange(j,
1213 j + length,
1214 r,
1215 prec_r.data());
1217 for (unsigned int i = 0; i < length; ++i)
1218 {
1219 x[i] += this->alpha * p[i];
1220 p[i] = prec_r[i];
1221 v[i] = Number();
1222 }
1223 p += length;
1224 r += length;
1225 v += length;
1226 }
1227 }
1228 else
1229 {
1230 const Number alpha_plus_previous_alpha_over_beta =
1231 this->alpha + this->previous_alpha / this->previous_beta;
1232 const Number previous_alpha_over_beta =
1233 this->previous_alpha / this->previous_beta;
1234 for (unsigned int j = start_range; j < end_range; j += grain_size)
1235 {
1236 const unsigned int length = std::min(grain_size, end_range - j);
1237 this->preconditioner.apply_to_subrange(j,
1238 j + length,
1239 r,
1240 prec_r.data());
1242 for (unsigned int i = 0; i < length; ++i)
1243 {
1244 x[i] += alpha_plus_previous_alpha_over_beta * p[i] -
1245 previous_alpha_over_beta * prec_r[i];
1246 r[i] -= this->alpha * v[i];
1247 }
1248 this->preconditioner.apply_to_subrange(j,
1249 j + length,
1250 r,
1251 prec_r.data());
1253 for (unsigned int i = 0; i < length; ++i)
1254 {
1255 p[i] = this->beta * p[i] + prec_r[i];
1256 v[i] = Number();
1257 }
1258 p += length;
1259 r += length;
1260 v += length;
1261 x += length;
1262 }
1263 }
1264 }
1265
1266 // Function that we use if the PreconditionerType does not implement an
1267 // apply() function and where we instead need to use the
1268 // apply_to_subrange function
1269 template <typename U = void>
1270 std::enable_if_t<!has_apply<PreconditionerType>, U>
1271 operation_after_loop(
1272 const unsigned int start_range,
1273 const unsigned int end_range,
1274 std::array<VectorizedArray<Number>, 7> &vectorized_sums) const
1275 {
1276 const Number *r = this->r.begin();
1277 const Number *p = this->p.begin();
1278 const Number *v = this->v.begin();
1279 std::array<VectorizedArray<Number>, 7> my_sums = {};
1280 constexpr unsigned int grain_size = 128;
1281 Assert(grain_size % VectorizedArray<Number>::size() == 0,
1283 const unsigned int end_regular =
1284 start_range + (end_range - start_range) / grain_size * grain_size;
1285 std::array<Number, grain_size> prec_r;
1286 std::array<Number, grain_size> prec_v;
1287 for (unsigned int j = start_range; j < end_regular; j += grain_size)
1288 {
1289 this->preconditioner.apply_to_subrange(j,
1290 j + grain_size,
1291 r + j,
1292 prec_r.data());
1293 this->preconditioner.apply_to_subrange(j,
1294 j + grain_size,
1295 v + j,
1296 prec_v.data());
1297 VectorizedArray<Number> pj, vj, rj, prec_vj, prec_rj;
1298 for (unsigned int i = 0; i < grain_size;
1300 {
1301 pj.load(p + j + i);
1302 vj.load(v + j + i);
1303 rj.load(r + j + i);
1304 prec_rj.load(prec_r.data() + i);
1305 prec_vj.load(prec_v.data() + i);
1306
1307 my_sums[0] += pj * vj;
1308 my_sums[1] += vj * vj;
1309 my_sums[2] += rj * vj;
1310 my_sums[3] += rj * rj;
1311 my_sums[4] += rj * prec_vj;
1312 my_sums[5] += vj * prec_vj;
1313 my_sums[6] += rj * prec_rj;
1314 }
1315 }
1316 const unsigned int length = end_range - end_regular;
1317 AssertIndexRange(length, grain_size);
1318 this->preconditioner.apply_to_subrange(end_regular,
1319 end_regular + length,
1320 r + end_regular,
1321 prec_r.data());
1322 this->preconditioner.apply_to_subrange(end_regular,
1323 end_regular + length,
1324 v + end_regular,
1325 prec_v.data());
1326 for (unsigned int j = end_regular; j < end_range; ++j)
1327 {
1328 my_sums[0][0] += p[j] * v[j];
1329 my_sums[1][0] += v[j] * v[j];
1330 my_sums[2][0] += r[j] * v[j];
1331 my_sums[3][0] += r[j] * r[j];
1332 my_sums[4][0] += r[j] * prec_v[j - end_regular];
1333 my_sums[5][0] += v[j] * prec_v[j - end_regular];
1334 my_sums[6][0] += r[j] * prec_r[j - end_regular];
1335 }
1336 for (unsigned int i = 0; i < vectorized_sums.size(); ++i)
1337 vectorized_sums[i] += my_sums[i];
1338 }
1339
1340 // Function that we use if the PreconditionerType does not implement an
1341 // apply() function, where we instead need to choose the
1342 // apply_to_subrange function
1343 template <typename U = void>
1344 std::enable_if_t<!has_apply<PreconditionerType>, U>
1345 finalize_after_convergence(const unsigned int iteration_index)
1346 {
1347 if (iteration_index % 2 == 1 || this->beta == Number())
1348 this->x.add(this->alpha, this->p);
1349 else
1350 {
1351 const unsigned int end_range = this->x.locally_owned_size();
1352
1353 Number *x = this->x.begin();
1354 Number *r = this->r.begin();
1355 Number *p = this->p.begin();
1356 const Number alpha_plus_previous_alpha_over_beta =
1357 this->alpha + this->previous_alpha / this->beta;
1358 const Number previous_alpha_over_beta =
1359 this->previous_alpha / this->beta;
1360
1361 constexpr unsigned int grain_size = 128;
1362 std::array<Number, grain_size> prec_r;
1363 for (unsigned int j = 0; j < end_range; j += grain_size)
1364 {
1365 const unsigned int length = std::min(grain_size, end_range - j);
1366 this->preconditioner.apply_to_subrange(j,
1367 j + length,
1368 r,
1369 prec_r.data());
1371 for (unsigned int i = 0; i < length; ++i)
1372 x[i] += alpha_plus_previous_alpha_over_beta * p[i] -
1373 previous_alpha_over_beta * prec_r[i];
1374
1375 x += length;
1376 r += length;
1377 p += length;
1378 }
1379 }
1380 }
1381 };
1382 } // namespace SolverCG
1383} // namespace internal
1384
1385
1386
1387template <typename VectorType>
1389template <typename MatrixType, typename PreconditionerType>
1393void SolverCG<VectorType>::solve(const MatrixType &A,
1394 VectorType &x,
1395 const VectorType &b,
1396 const PreconditionerType &preconditioner)
1397{
1398 using number = typename VectorType::value_type;
1399
1401
1402 LogStream::Prefix prefix("cg");
1403
1404 // Should we build the matrix for eigenvalue computations?
1405 const bool do_eigenvalues =
1406 !condition_number_signal.empty() || !all_condition_numbers_signal.empty() ||
1407 !eigenvalues_signal.empty() || !all_eigenvalues_signal.empty();
1408
1409 // vectors used for eigenvalue computations
1410 std::vector<typename VectorType::value_type> diagonal;
1411 std::vector<typename VectorType::value_type> offdiagonal;
1412
1413 typename VectorType::value_type eigen_beta_alpha = 0;
1414
1415 int it = 0;
1416
1417 internal::SolverCG::
1418 IterationWorker<VectorType, MatrixType, PreconditionerType>
1419 worker(A,
1420 preconditioner,
1421 determine_beta_by_flexible_formula,
1422 this->memory,
1423 x,
1424 b,
1425 additional_data.use_default_residual);
1426
1427 worker.startup();
1428
1429 solver_state = this->iteration_status(0, worker.residual_norm, x);
1430 if (solver_state != SolverControl::iterate)
1431 return;
1432
1433 while (solver_state == SolverControl::iterate)
1434 {
1435 ++it;
1436
1437 worker.do_iteration(it);
1438
1439 print_vectors(it, x, worker.r, worker.p);
1440
1441 if (it > 1)
1442 {
1443 this->coefficients_signal(worker.previous_alpha, worker.beta);
1444 // set up the vectors containing the diagonal and the off diagonal
1445 // of the projected matrix.
1446 if (do_eigenvalues)
1447 {
1448 diagonal.push_back(number(1.) / worker.previous_alpha +
1449 eigen_beta_alpha);
1450 eigen_beta_alpha = worker.beta / worker.previous_alpha;
1451 offdiagonal.push_back(std::sqrt(worker.beta) /
1452 worker.previous_alpha);
1453 }
1454 compute_eigs_and_cond(diagonal,
1455 offdiagonal,
1456 all_eigenvalues_signal,
1457 all_condition_numbers_signal);
1458 }
1459
1460 solver_state = this->iteration_status(it, worker.residual_norm, x);
1461 }
1462
1463 worker.finalize_after_convergence(it);
1464
1465 compute_eigs_and_cond(diagonal,
1466 offdiagonal,
1467 eigenvalues_signal,
1468 condition_number_signal);
1469
1470 AssertThrow(solver_state == SolverControl::success,
1471 SolverControl::NoConvergence(it, worker.residual_norm));
1472}
1473
1474
1475
1476template <typename VectorType>
1478boost::signals2::connection SolverCG<VectorType>::connect_coefficients_slot(
1479 const std::function<void(typename VectorType::value_type,
1480 typename VectorType::value_type)> &slot)
1481{
1482 return coefficients_signal.connect(slot);
1483}
1484
1485
1486
1487template <typename VectorType>
1489boost::signals2::connection SolverCG<VectorType>::connect_condition_number_slot(
1490 const std::function<void(double)> &slot,
1491 const bool every_iteration)
1492{
1493 if (every_iteration)
1494 {
1495 return all_condition_numbers_signal.connect(slot);
1496 }
1497 else
1498 {
1499 return condition_number_signal.connect(slot);
1500 }
1501}
1502
1503
1504
1505template <typename VectorType>
1507boost::signals2::connection SolverCG<VectorType>::connect_eigenvalues_slot(
1508 const std::function<void(const std::vector<double> &)> &slot,
1509 const bool every_iteration)
1510{
1511 if (every_iteration)
1512 {
1513 return all_eigenvalues_signal.connect(slot);
1514 }
1515 else
1516 {
1517 return eigenvalues_signal.connect(slot);
1518 }
1519}
1520
1521
1522
1523template <typename VectorType>
1526 const bool use_default_residual)
1527 : use_default_residual(use_default_residual)
1528{}
1529
1530
1531template <typename VectorType>
1535 const AdditionalData &data)
1536 : SolverCG<VectorType>(cn, mem)
1537{
1538 this->determine_beta_by_flexible_formula = true;
1539 this->additional_data = data;
1540}
1541
1542
1543
1544template <typename VectorType>
1547 const AdditionalData &data)
1548 : SolverCG<VectorType>(cn)
1549{
1550 this->determine_beta_by_flexible_formula = true;
1551 this->additional_data = data;
1552}
1553
1554
1555#endif // DOXYGEN
1556
1558
1559#endif
SolverCG(SolverControl &cn, VectorMemory< VectorType > &mem, const AdditionalData &data=AdditionalData())
virtual void print_vectors(const unsigned int step, const VectorType &x, const VectorType &r, const VectorType &d) const
boost::signals2::connection connect_condition_number_slot(const std::function< void(double)> &slot, const bool every_iteration=false)
boost::signals2::connection connect_eigenvalues_slot(const std::function< void(const std::vector< double > &)> &slot, const bool every_iteration=false)
boost::signals2::connection connect_coefficients_slot(const std::function< void(typename VectorType::value_type, typename VectorType::value_type)> &slot)
virtual ~SolverCG() override=default
void solve(const MatrixType &A, VectorType &x, const VectorType &b, const PreconditionerType &preconditioner)
SolverCG(SolverControl &cn, const AdditionalData &data=AdditionalData())
static void compute_eigs_and_cond(const std::vector< typename VectorType::value_type > &diagonal, const std::vector< typename VectorType::value_type > &offdiagonal, const boost::signals2::signal< void(const std::vector< double > &)> &eigenvalues_signal, const boost::signals2::signal< void(double)> &cond_signal)
@ iterate
Continue iteration.
@ success
Stop iteration, goal reached.
AdditionalData additional_data
Definition solver_cg.h:454
SolverFlexibleCG(SolverControl &cn, const AdditionalData &data=AdditionalData())
SolverFlexibleCG(SolverControl &cn, VectorMemory< VectorType > &mem, const AdditionalData &data=AdditionalData())
void store(OtherNumber *ptr) const
void load(const OtherNumber *ptr)
#define DEAL_II_OPENMP_SIMD_PRAGMA
Definition config.h:208
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:35
#define DEAL_II_CXX20_REQUIRES(condition)
Definition config.h:242
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:36
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcDivideByZero()
#define AssertThrow(cond, exc)
std::vector< index_type > data
Definition mpi.cc:746
@ diagonal
Matrix is diagonal.
constexpr char U
constexpr char T
constexpr char A
Tpetra::Vector< Number, LO, GO, NodeType< MemorySpace > > VectorType
Tpetra::CrsMatrix< Number, LO, GO, NodeType< MemorySpace > > MatrixType
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
VectorType::value_type * end(VectorType &V)
VectorType::value_type * begin(VectorType &V)
T sum(const T &t, const MPI_Comm mpi_communicator)
STL namespace.
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
unsigned int global_dof_index
Definition types.h:94
AdditionalData(const bool use_default_residual=true)
AdditionalData(const bool use_default_residual=true)
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)