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