Reference documentation for deal.II version GIT relicensing-37-gdcd04c9418 2024-03-02 19:00:02+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
25
26#include <deal.II/lac/solver.h>
29
30#include <cmath>
31
33
34// forward declaration
35#ifndef DOXYGEN
37namespace LinearAlgebra
38{
39 namespace distributed
40 {
41 template <typename, typename>
42 class Vector;
43 }
44} // namespace LinearAlgebra
45#endif
46
47
175template <typename VectorType = Vector<double>>
176class SolverCG : public SolverBase<VectorType>
177{
178public:
183
190 {};
191
197 const AdditionalData &data = AdditionalData());
198
204
208 virtual ~SolverCG() override = default;
209
213 template <typename MatrixType, typename PreconditionerType>
214 void
215 solve(const MatrixType &A,
216 VectorType &x,
217 const VectorType &b,
218 const PreconditionerType &preconditioner);
219
226 boost::signals2::connection
228 const std::function<void(typename VectorType::value_type,
229 typename VectorType::value_type)> &slot);
230
237 boost::signals2::connection
238 connect_condition_number_slot(const std::function<void(double)> &slot,
239 const bool every_iteration = false);
240
247 boost::signals2::connection
249 const std::function<void(const std::vector<double> &)> &slot,
250 const bool every_iteration = false);
251
252protected:
258 virtual void
259 print_vectors(const unsigned int step,
260 const VectorType &x,
261 const VectorType &r,
262 const VectorType &d) const;
263
269 static void
271 const std::vector<typename VectorType::value_type> &diagonal,
272 const std::vector<typename VectorType::value_type> &offdiagonal,
273 const boost::signals2::signal<void(const std::vector<double> &)>
275 const boost::signals2::signal<void(double)> &cond_signal);
276
281
285 boost::signals2::signal<void(typename VectorType::value_type,
286 typename VectorType::value_type)>
288
293 boost::signals2::signal<void(double)> condition_number_signal;
294
299 boost::signals2::signal<void(double)> all_condition_numbers_signal;
300
305 boost::signals2::signal<void(const std::vector<double> &)> eigenvalues_signal;
306
311 boost::signals2::signal<void(const std::vector<double> &)>
313
323};
324
325
326
352template <typename VectorType = Vector<double>>
353class SolverFlexibleCG : public SolverCG<VectorType>
354{
355public:
360
367 {};
368
374 const AdditionalData &data = AdditionalData());
375
381 const AdditionalData &data = AdditionalData());
382};
383
384
387/*------------------------- Implementation ----------------------------*/
388
389#ifndef DOXYGEN
390
391
392
393template <typename VectorType>
396 const AdditionalData &data)
397 : SolverBase<VectorType>(cn, mem)
398 , additional_data(data)
399 , determine_beta_by_flexible_formula(false)
400{}
401
402
403
404template <typename VectorType>
405SolverCG<VectorType>::SolverCG(SolverControl &cn, const AdditionalData &data)
406 : SolverBase<VectorType>(cn)
407 , additional_data(data)
408 , determine_beta_by_flexible_formula(false)
409{}
410
411
412
413template <typename VectorType>
414void
415SolverCG<VectorType>::print_vectors(const unsigned int,
416 const VectorType &,
417 const VectorType &,
418 const VectorType &) const
419{}
420
421
422
423template <typename VectorType>
424inline void
426 const std::vector<typename VectorType::value_type> &diagonal,
427 const std::vector<typename VectorType::value_type> &offdiagonal,
428 const boost::signals2::signal<void(const std::vector<double> &)>
429 &eigenvalues_signal,
430 const boost::signals2::signal<void(double)> &cond_signal)
431{
432 // Avoid computing eigenvalues unless they are needed.
433 if (!cond_signal.empty() || !eigenvalues_signal.empty())
434 {
436 true);
437 for (size_type i = 0; i < diagonal.size(); ++i)
438 {
439 T(i, i) = diagonal[i];
440 if (i < diagonal.size() - 1)
441 T(i, i + 1) = offdiagonal[i];
442 }
443 T.compute_eigenvalues();
444 // Need two eigenvalues to estimate the condition number.
445 if (diagonal.size() > 1)
446 {
447 auto condition_number = T.eigenvalue(T.n() - 1) / T.eigenvalue(0);
448 // Condition number is real valued and nonnegative; simply take
449 // the absolute value:
450 cond_signal(std::abs(condition_number));
451 }
452 // Avoid copying the eigenvalues of T to a vector unless a signal is
453 // connected.
454 if (!eigenvalues_signal.empty())
455 {
456 std::vector<double> eigenvalues(T.n());
457 for (unsigned int j = 0; j < T.n(); ++j)
458 {
459 // for a hermitian matrix, all eigenvalues are real-valued
460 // and non-negative, simply return the absolute value:
461 eigenvalues[j] = std::abs(T.eigenvalue(j));
462 }
463 eigenvalues_signal(eigenvalues);
464 }
465 }
466}
467
468
469
470namespace internal
471{
472 namespace SolverCG
473 {
474 // This base class is used to select different variants of the conjugate
475 // gradient solver. The default variant is used for standard matrix and
476 // preconditioner arguments, as provided by the derived class
477 // IterationWork below, but there is also a specialized variant further
478 // down that uses SFINAE to identify whether matrices and preconditioners
479 // support special operations on sub-ranges of the vectors.
480 template <typename VectorType,
481 typename MatrixType,
482 typename PreconditionerType>
483 struct IterationWorkerBase
484 {
485 using Number = typename VectorType::value_type;
486
487 const MatrixType &A;
488 const PreconditionerType &preconditioner;
489 const bool flexible;
490 VectorType &x;
491
492 typename VectorMemory<VectorType>::Pointer r_pointer;
493 typename VectorMemory<VectorType>::Pointer p_pointer;
494 typename VectorMemory<VectorType>::Pointer v_pointer;
495 typename VectorMemory<VectorType>::Pointer z_pointer;
496
497 // Define some aliases for simpler access, using the variables 'r' for
498 // the residual b - A*x, 'p' for the search direction, and 'v' for the
499 // auxiliary vector. This naming convention is used e.g. by the
500 // description on
501 // https://en.wikipedia.org/wiki/Conjugate_gradient_method. The variable
502 // 'z' gets only used for the flexible variant of the CG method.
503 VectorType &r;
504 VectorType &p;
505 VectorType &v;
506 VectorType &z;
507
508 Number r_dot_preconditioner_dot_r;
509 Number alpha;
510 Number beta;
511 double residual_norm;
512 Number previous_alpha;
513
514 IterationWorkerBase(const MatrixType &A,
515 const PreconditionerType &preconditioner,
516 const bool flexible,
518 VectorType &x)
519 : A(A)
520 , preconditioner(preconditioner)
521 , flexible(flexible)
522 , x(x)
523 , r_pointer(memory)
524 , p_pointer(memory)
525 , v_pointer(memory)
526 , z_pointer(memory)
527 , r(*r_pointer)
528 , p(*p_pointer)
529 , v(*v_pointer)
530 , z(*z_pointer)
531 , r_dot_preconditioner_dot_r(Number())
532 , alpha(Number())
533 , beta(Number())
534 , residual_norm(0.0)
535 , previous_alpha(Number())
536 {}
537
538 void
539 startup(const VectorType &b)
540 {
541 // Initialize without setting the vector entries, as those would soon
542 // be overwritten anyway
543 r.reinit(x, true);
544 p.reinit(x, true);
545 v.reinit(x, true);
546 if (flexible)
547 z.reinit(x, true);
548
549 // compute residual. if vector is zero, then short-circuit the full
550 // computation
551 if (!x.all_zero())
552 {
553 A.vmult(r, x);
554 r.sadd(-1., 1., b);
555 }
556 else
557 r.equ(1., b);
558
559 residual_norm = r.l2_norm();
560 }
561 };
562
563
564
565 // Implementation of a conjugate gradient operation with matrices and
566 // preconditioners without special capabilities
567 template <typename VectorType,
568 typename MatrixType,
569 typename PreconditionerType,
570 typename = int>
571 struct IterationWorker
572 : public IterationWorkerBase<VectorType, MatrixType, PreconditionerType>
573 {
574 using BaseClass =
575 IterationWorkerBase<VectorType, MatrixType, PreconditionerType>;
576
577 IterationWorker(const MatrixType &A,
578 const PreconditionerType &preconditioner,
579 const bool flexible,
581 VectorType &x)
582 : BaseClass(A, preconditioner, flexible, memory, x)
583 {}
584
585 using BaseClass::A;
586 using BaseClass::alpha;
587 using BaseClass::beta;
588 using BaseClass::p;
589 using BaseClass::preconditioner;
590 using BaseClass::r;
591 using BaseClass::r_dot_preconditioner_dot_r;
592 using BaseClass::residual_norm;
593 using BaseClass::v;
594 using BaseClass::x;
595 using BaseClass::z;
596
597 void
598 do_iteration(const unsigned int iteration_index)
599 {
600 using Number = typename VectorType::value_type;
601
602 const Number previous_r_dot_preconditioner_dot_r =
603 r_dot_preconditioner_dot_r;
604
605 if (std::is_same_v<PreconditionerType, PreconditionIdentity> == false)
606 {
607 preconditioner.vmult(v, r);
608 r_dot_preconditioner_dot_r = r * v;
609 }
610 else
611 r_dot_preconditioner_dot_r = residual_norm * residual_norm;
612
613 const VectorType &direction =
614 std::is_same_v<PreconditionerType, PreconditionIdentity> ? r : v;
615
616 if (iteration_index > 1)
617 {
618 Assert(std::abs(previous_r_dot_preconditioner_dot_r) != 0.,
620 beta =
621 r_dot_preconditioner_dot_r / previous_r_dot_preconditioner_dot_r;
622 if (this->flexible)
623 beta -= (r * z) / previous_r_dot_preconditioner_dot_r;
624 p.sadd(beta, 1., direction);
625 }
626 else
627 p.equ(1., direction);
628
629 if (this->flexible)
630 z.swap(v);
631
632 A.vmult(v, p);
633
634 const Number p_dot_A_dot_p = p * v;
635 Assert(std::abs(p_dot_A_dot_p) != 0., ExcDivideByZero());
636
637 this->previous_alpha = alpha;
638 alpha = r_dot_preconditioner_dot_r / p_dot_A_dot_p;
639
640 x.add(alpha, p);
641 residual_norm = std::sqrt(std::abs(r.add_and_dot(-alpha, v, r)));
642 }
643
644 void
645 finalize_after_convergence(const unsigned int)
646 {}
647 };
648
649
650 // In the following, we provide a specialization of the above
651 // IterationWorker class that picks up particular features in the matrix
652 // and preconditioners.
653
654 // a helper type-trait that leverage SFINAE to figure out if MatrixType has
655 // ... MatrixType::vmult(VectorType &, const VectorType&,
656 // std::function<...>, std::function<...>) const
657 template <typename MatrixType, typename VectorType>
658 using vmult_functions_t = decltype(std::declval<const MatrixType>().vmult(
659 std::declval<VectorType &>(),
660 std::declval<const VectorType &>(),
661 std::declval<
662 const std::function<void(const unsigned int, const unsigned int)> &>(),
663 std::declval<const std::function<void(const unsigned int,
664 const unsigned int)> &>()));
665
666 template <typename MatrixType, typename VectorType>
667 constexpr bool has_vmult_functions =
668 is_supported_operation<vmult_functions_t, MatrixType, VectorType>;
669
670 // a helper type-trait that leverage SFINAE to figure out if
671 // PreconditionerType has ... PreconditionerType::apply_to_subrange(const
672 // unsigned int, const unsigned int, const Number*, Number*) const
673 template <typename PreconditionerType>
674 using apply_to_subrange_t =
675 decltype(std::declval<const PreconditionerType>()
676 .apply_to_subrange(0U, 0U, nullptr, nullptr));
677
678 template <typename PreconditionerType>
679 constexpr bool has_apply_to_subrange =
680 is_supported_operation<apply_to_subrange_t, PreconditionerType>;
681
682 // a helper type-trait that leverage SFINAE to figure out if
683 // PreconditionerType has ... PreconditionerType::apply(const
684 // unsigned int, const Number) const
685 template <typename PreconditionerType>
686 using apply_t =
687 decltype(std::declval<const PreconditionerType>().apply(0U, 0.0));
688
689 template <typename PreconditionerType>
690 constexpr bool has_apply =
691 is_supported_operation<apply_t, PreconditionerType>;
692
693
694 // Internal function to run one iteration of the conjugate gradient solver
695 // for matrices and preconditioners that support interleaving the vector
696 // updates with the matrix-vector product.
697 template <typename VectorType,
698 typename MatrixType,
699 typename PreconditionerType>
700 struct IterationWorker<
701 VectorType,
702 MatrixType,
703 PreconditionerType,
704 std::enable_if_t<has_vmult_functions<MatrixType, VectorType> &&
705 (has_apply_to_subrange<PreconditionerType> ||
706 has_apply<PreconditionerType>)&&std::
707 is_same_v<VectorType,
708 LinearAlgebra::distributed::Vector<
709 typename VectorType::value_type,
710 MemorySpace::Host>>,
711 int>>
712 : public IterationWorkerBase<VectorType, MatrixType, PreconditionerType>
713 {
714 using Number = typename VectorType::value_type;
715
716 Number next_r_dot_preconditioner_dot_r;
717 Number previous_beta;
718
719 IterationWorker(const MatrixType &A,
720 const PreconditionerType &preconditioner,
721 const bool flexible,
723 VectorType &x)
724 : IterationWorkerBase<VectorType, MatrixType, PreconditionerType>(
725 A,
726 preconditioner,
727 flexible,
728 memory,
729 x)
730 , next_r_dot_preconditioner_dot_r(0.)
731 , previous_beta(0.)
732 {}
733
734 // This is the main iteration function, that will use some of the
735 // specialized functions below
736 void
737 do_iteration(const unsigned int iteration_index)
738 {
739 if (iteration_index > 1)
740 {
741 previous_beta = this->beta;
742 this->beta = next_r_dot_preconditioner_dot_r /
743 this->r_dot_preconditioner_dot_r;
744 }
745
746 std::array<VectorizedArray<Number>, 7> vectorized_sums = {};
747
748 this->A.vmult(
749 this->v,
750 this->p,
751 [&](const unsigned int begin, const unsigned int end) {
752 operation_before_loop(iteration_index, begin, end);
753 },
754 [&](const unsigned int begin, const unsigned int end) {
755 operation_after_loop(begin, end, vectorized_sums);
756 });
757
758 std::array<Number, 7> scalar_sums;
759 for (unsigned int i = 0; i < 7; ++i)
760 scalar_sums[i] = vectorized_sums[i][0];
761 for (unsigned int l = 1; l < VectorizedArray<Number>::size(); ++l)
762 for (unsigned int i = 0; i < 7; ++i)
763 scalar_sums[i] += vectorized_sums[i][l];
764
766 7),
767 this->r.get_mpi_communicator(),
768 ::ArrayView<Number>(scalar_sums.data(), 7));
769
770 this->r_dot_preconditioner_dot_r = scalar_sums[6];
771
772 const Number p_dot_A_dot_p = scalar_sums[0];
773 Assert(std::abs(p_dot_A_dot_p) != 0., ExcDivideByZero());
774
775 this->previous_alpha = this->alpha;
776 this->alpha = this->r_dot_preconditioner_dot_r / p_dot_A_dot_p;
777
778 // Round-off errors near zero might yield negative values, so take
779 // the absolute value in the next two formulas
780 this->residual_norm = std::sqrt(std::abs(
781 scalar_sums[3] +
782 this->alpha * (-2. * scalar_sums[2] + this->alpha * scalar_sums[1])));
783
784 next_r_dot_preconditioner_dot_r = std::abs(
785 this->r_dot_preconditioner_dot_r +
786 this->alpha * (-2. * scalar_sums[4] + this->alpha * scalar_sums[5]));
787 }
788
789 // Function that we use if the PreconditionerType implements an apply()
790 // function
791 template <typename U = void>
792 std::enable_if_t<has_apply<PreconditionerType>, U>
793 operation_before_loop(const unsigned int iteration_index,
794 const unsigned int start_range,
795 const unsigned int end_range) const
796 {
797 Number *x = this->x.begin();
798 Number *r = this->r.begin();
799 Number *p = this->p.begin();
800 Number *v = this->v.begin();
801 const Number alpha = this->alpha;
802 const Number beta = this->beta;
803 constexpr unsigned int n_lanes = VectorizedArray<Number>::size();
804 const unsigned int end_regular =
805 start_range + (end_range - start_range) / n_lanes * n_lanes;
806 if (iteration_index == 1)
807 {
808 // Vectorize by hand since compilers are often pretty bad at
809 // doing these steps automatically even with
810 // DEAL_II_OPENMP_SIMD_PRAGMA
811 for (unsigned int j = start_range; j < end_regular; j += n_lanes)
812 {
814 rj.load(r + j);
816 for (unsigned int l = 0; l < n_lanes; ++l)
817 pj[l] = this->preconditioner.apply(j + l, rj[l]);
818 pj.store(p + j);
820 rj.store(v + j);
821 }
822 for (unsigned int j = end_regular; j < end_range; ++j)
823 {
824 p[j] = this->preconditioner.apply(j, r[j]);
825 v[j] = Number();
826 }
827 }
828 else if (iteration_index % 2 == 0 && beta != Number())
829 {
830 for (unsigned int j = start_range; j < end_regular; j += n_lanes)
831 {
832 VectorizedArray<Number> rj, vj, pj, prec_rj;
833 rj.load(r + j);
834 vj.load(v + j);
835 rj -= alpha * vj;
836 rj.store(r + j);
838 for (unsigned int l = 0; l < n_lanes; ++l)
839 prec_rj[l] = this->preconditioner.apply(j + l, rj[l]);
840 pj.load(p + j);
841 pj = beta * pj + prec_rj;
842 pj.store(p + j);
844 rj.store(v + j);
845 }
846 for (unsigned int j = end_regular; j < end_range; ++j)
847 {
848 r[j] -= alpha * v[j];
849 p[j] = beta * p[j] + this->preconditioner.apply(j, r[j]);
850 v[j] = Number();
851 }
852 }
853 else if (iteration_index % 2 == 0 && beta == Number())
854 {
855 // Case where beta is zero: we cannot reconstruct p_{j-1} in the
856 // next iteration, and must hence compute x here. This can happen
857 // before termination
858 for (unsigned int j = start_range; j < end_regular; j += n_lanes)
859 {
860 VectorizedArray<Number> rj, xj, vj, pj, prec_rj;
861 rj.load(r + j);
862 vj.load(v + j);
863 xj.load(x + j);
864 pj.load(p + j);
865 xj += alpha * pj;
866 xj.store(x + j);
867 rj -= alpha * vj;
868 rj.store(r + j);
870 for (unsigned int l = 0; l < n_lanes; ++l)
871 prec_rj[l] = this->preconditioner.apply(j + l, rj[l]);
872 prec_rj.store(p + j);
874 rj.store(v + j);
875 }
876 for (unsigned int j = end_regular; j < end_range; ++j)
877 {
878 r[j] -= alpha * v[j];
879 x[j] += alpha * p[j];
880 p[j] = this->preconditioner.apply(j, r[j]);
881 v[j] = Number();
882 }
883 }
884 else
885 {
886 const Number alpha_plus_previous_alpha_over_beta =
887 alpha + this->previous_alpha / this->previous_beta;
888 const Number previous_alpha_over_beta =
889 this->previous_alpha / this->previous_beta;
890 for (unsigned int j = start_range; j < end_regular; j += n_lanes)
891 {
892 VectorizedArray<Number> rj, vj, pj, xj, prec_rj, prec_vj;
893 xj.load(x + j);
894 pj.load(p + j);
895 xj += alpha_plus_previous_alpha_over_beta * pj;
896 rj.load(r + j);
897 vj.load(v + j);
899 for (unsigned int l = 0; l < n_lanes; ++l)
900 {
901 prec_rj[l] = this->preconditioner.apply(j + l, rj[l]);
902 prec_vj[l] = this->preconditioner.apply(j + l, vj[l]);
903 }
904 xj -= previous_alpha_over_beta * prec_rj;
905 xj.store(x + j);
906 rj -= alpha * vj;
907 rj.store(r + j);
908 prec_rj -= alpha * prec_vj;
909 pj = beta * pj + prec_rj;
910 pj.store(p + j);
912 rj.store(v + j);
913 }
914 for (unsigned int j = end_regular; j < end_range; ++j)
915 {
916 x[j] += alpha_plus_previous_alpha_over_beta * p[j];
917 x[j] -= previous_alpha_over_beta *
918 this->preconditioner.apply(j, r[j]);
919 r[j] -= alpha * v[j];
920 p[j] = beta * p[j] + this->preconditioner.apply(j, r[j]);
921 v[j] = Number();
922 }
923 }
924 }
925
926 // Function that we use if the PreconditionerType implements an apply()
927 // function
928 template <typename U = void>
929 std::enable_if_t<has_apply<PreconditionerType>, U>
930 operation_after_loop(
931 const unsigned int start_range,
932 const unsigned int end_range,
933 std::array<VectorizedArray<Number>, 7> &vectorized_sums) const
934 {
935 const Number *r = this->r.begin();
936 const Number *p = this->p.begin();
937 const Number *v = this->v.begin();
938 std::array<VectorizedArray<Number>, 7> my_sums = {};
939 constexpr unsigned int n_lanes = VectorizedArray<Number>::size();
940 const unsigned int end_regular =
941 start_range + (end_range - start_range) / n_lanes * n_lanes;
942 for (unsigned int j = start_range; j < end_regular; j += n_lanes)
943 {
944 VectorizedArray<Number> pj, vj, rj, prec_vj, prec_rj;
945 pj.load(p + j);
946 vj.load(v + j);
947 rj.load(r + j);
949 for (unsigned int l = 0; l < n_lanes; ++l)
950 {
951 prec_vj[l] = this->preconditioner.apply(j + l, vj[l]);
952 prec_rj[l] = this->preconditioner.apply(j + l, rj[l]);
953 }
954 my_sums[0] += pj * vj;
955 my_sums[1] += vj * vj;
956 my_sums[2] += rj * vj;
957 my_sums[3] += rj * rj;
958 my_sums[4] += rj * prec_vj;
959 my_sums[5] += vj * prec_vj;
960 my_sums[6] += rj * prec_rj;
961 }
962 for (unsigned int j = end_regular; j < end_range; ++j)
963 {
964 const Number prec_v = this->preconditioner.apply(j, v[j]);
965 const Number prec_r = this->preconditioner.apply(j, r[j]);
966 my_sums[0][0] += p[j] * v[j];
967 my_sums[1][0] += v[j] * v[j];
968 my_sums[2][0] += r[j] * v[j];
969 my_sums[3][0] += r[j] * r[j];
970 my_sums[4][0] += r[j] * prec_v;
971 my_sums[5][0] += v[j] * prec_v;
972 my_sums[6][0] += r[j] * prec_r;
973 }
974 for (unsigned int i = 0; i < vectorized_sums.size(); ++i)
975 vectorized_sums[i] += my_sums[i];
976 }
977
978 // Function that we use if the PreconditionerType implements an apply()
979 // function
980 template <typename U = void>
981 std::enable_if_t<has_apply<PreconditionerType>, U>
982 finalize_after_convergence(const unsigned int iteration_index)
983 {
984 if (iteration_index % 2 == 1 || this->beta == Number())
985 this->x.add(this->alpha, this->p);
986 else
987 {
988 using Number = typename VectorType::value_type;
989 const unsigned int end_range = this->x.locally_owned_size();
990
991 Number *x = this->x.begin();
992 Number *r = this->r.begin();
993 Number *p = this->p.begin();
994
995 // Note that we use 'beta' here rather than 'previous_beta' in the
996 // formula above, which is because the shift in beta ->
997 // previous_beta has not been applied at this stage, allowing us
998 // to recover the previous search direction
999 const Number alpha_plus_previous_alpha_over_beta =
1000 this->alpha + this->previous_alpha / this->beta;
1001 const Number previous_alpha_over_beta =
1002 this->previous_alpha / this->beta;
1003
1005 for (unsigned int j = 0; j < end_range; ++j)
1006 {
1007 x[j] += alpha_plus_previous_alpha_over_beta * p[j] -
1008 previous_alpha_over_beta *
1009 this->preconditioner.apply(j, r[j]);
1010 }
1011 }
1012 }
1013
1014 // Function that we use if the PreconditionerType does not implement an
1015 // apply() function, where we instead need to choose the
1016 // apply_to_subrange function
1017 template <typename U = void>
1018 std::enable_if_t<!has_apply<PreconditionerType>, U>
1019 operation_before_loop(const unsigned int iteration_index,
1020 const unsigned int start_range,
1021 const unsigned int end_range) const
1022 {
1023 Number *x = this->x.begin() + start_range;
1024 Number *r = this->r.begin() + start_range;
1025 Number *p = this->p.begin() + start_range;
1026 Number *v = this->v.begin() + start_range;
1027 const Number alpha = this->alpha;
1028 const Number beta = this->beta;
1029 constexpr unsigned int grain_size = 128;
1030 std::array<Number, grain_size> prec_r;
1031 if (iteration_index == 1)
1032 {
1033 for (unsigned int j = start_range; j < end_range; j += grain_size)
1034 {
1035 const unsigned int length = std::min(grain_size, end_range - j);
1036 this->preconditioner.apply_to_subrange(j,
1037 j + length,
1038 r,
1039 prec_r.data());
1041 for (unsigned int i = 0; i < length; ++i)
1042 {
1043 p[i] = prec_r[i];
1044 v[i] = Number();
1045 }
1046 p += length;
1047 r += length;
1048 v += length;
1049 }
1050 }
1051 else if (iteration_index % 2 == 0 && beta != Number())
1052 {
1053 for (unsigned int j = start_range; j < end_range; j += grain_size)
1054 {
1055 const unsigned int length = std::min(grain_size, end_range - j);
1057 for (unsigned int i = 0; i < length; ++i)
1058 r[i] -= this->alpha * v[i];
1059 this->preconditioner.apply_to_subrange(j,
1060 j + length,
1061 r,
1062 prec_r.data());
1064 for (unsigned int i = 0; i < length; ++i)
1065 {
1066 p[i] = this->beta * p[i] + prec_r[i];
1067 v[i] = Number();
1068 }
1069 p += length;
1070 r += length;
1071 v += length;
1072 }
1073 }
1074 else if (iteration_index % 2 == 0 && beta == Number())
1075 {
1076 // Case where beta is zero: We cannot reconstruct p_{j-1} in the
1077 // next iteration, and must hence compute x here. This can happen
1078 // before termination
1079 for (unsigned int j = start_range; j < end_range; j += grain_size)
1080 {
1081 const unsigned int length = std::min(grain_size, end_range - j);
1083 for (unsigned int i = 0; i < length; ++i)
1084 r[i] -= this->alpha * v[i];
1085 this->preconditioner.apply_to_subrange(j,
1086 j + length,
1087 r,
1088 prec_r.data());
1090 for (unsigned int i = 0; i < length; ++i)
1091 {
1092 x[i] += this->alpha * p[i];
1093 p[i] = prec_r[i];
1094 v[i] = Number();
1095 }
1096 p += length;
1097 r += length;
1098 v += length;
1099 }
1100 }
1101 else
1102 {
1103 const Number alpha_plus_previous_alpha_over_beta =
1104 this->alpha + this->previous_alpha / this->previous_beta;
1105 const Number previous_alpha_over_beta =
1106 this->previous_alpha / this->previous_beta;
1107 for (unsigned int j = start_range; j < end_range; j += grain_size)
1108 {
1109 const unsigned int length = std::min(grain_size, end_range - j);
1110 this->preconditioner.apply_to_subrange(j,
1111 j + length,
1112 r,
1113 prec_r.data());
1115 for (unsigned int i = 0; i < length; ++i)
1116 {
1117 x[i] += alpha_plus_previous_alpha_over_beta * p[i] -
1118 previous_alpha_over_beta * prec_r[i];
1119 r[i] -= this->alpha * v[i];
1120 }
1121 this->preconditioner.apply_to_subrange(j,
1122 j + length,
1123 r,
1124 prec_r.data());
1126 for (unsigned int i = 0; i < length; ++i)
1127 {
1128 p[i] = this->beta * p[i] + prec_r[i];
1129 v[i] = Number();
1130 }
1131 p += length;
1132 r += length;
1133 v += length;
1134 x += length;
1135 }
1136 }
1137 }
1138
1139 // Function that we use if the PreconditionerType does not implement an
1140 // apply() function and where we instead need to use the
1141 // apply_to_subrange function
1142 template <typename U = void>
1143 std::enable_if_t<!has_apply<PreconditionerType>, U>
1144 operation_after_loop(
1145 const unsigned int start_range,
1146 const unsigned int end_range,
1147 std::array<VectorizedArray<Number>, 7> &vectorized_sums) const
1148 {
1149 const Number *r = this->r.begin();
1150 const Number *p = this->p.begin();
1151 const Number *v = this->v.begin();
1152 std::array<VectorizedArray<Number>, 7> my_sums = {};
1153 constexpr unsigned int grain_size = 128;
1154 Assert(grain_size % VectorizedArray<Number>::size() == 0,
1156 const unsigned int end_regular =
1157 start_range + (end_range - start_range) / grain_size * grain_size;
1158 std::array<Number, grain_size> prec_r;
1159 std::array<Number, grain_size> prec_v;
1160 for (unsigned int j = start_range; j < end_regular; j += grain_size)
1161 {
1162 this->preconditioner.apply_to_subrange(j,
1163 j + grain_size,
1164 r + j,
1165 prec_r.data());
1166 this->preconditioner.apply_to_subrange(j,
1167 j + grain_size,
1168 v + j,
1169 prec_v.data());
1170 VectorizedArray<Number> pj, vj, rj, prec_vj, prec_rj;
1171 for (unsigned int i = 0; i < grain_size;
1173 {
1174 pj.load(p + j + i);
1175 vj.load(v + j + i);
1176 rj.load(r + j + i);
1177 prec_rj.load(prec_r.data() + i);
1178 prec_vj.load(prec_v.data() + i);
1179
1180 my_sums[0] += pj * vj;
1181 my_sums[1] += vj * vj;
1182 my_sums[2] += rj * vj;
1183 my_sums[3] += rj * rj;
1184 my_sums[4] += rj * prec_vj;
1185 my_sums[5] += vj * prec_vj;
1186 my_sums[6] += rj * prec_rj;
1187 }
1188 }
1189 const unsigned int length = end_range - end_regular;
1190 AssertIndexRange(length, grain_size);
1191 this->preconditioner.apply_to_subrange(end_regular,
1192 end_regular + length,
1193 r + end_regular,
1194 prec_r.data());
1195 this->preconditioner.apply_to_subrange(end_regular,
1196 end_regular + length,
1197 v + end_regular,
1198 prec_v.data());
1199 for (unsigned int j = end_regular; j < end_range; ++j)
1200 {
1201 my_sums[0][0] += p[j] * v[j];
1202 my_sums[1][0] += v[j] * v[j];
1203 my_sums[2][0] += r[j] * v[j];
1204 my_sums[3][0] += r[j] * r[j];
1205 my_sums[4][0] += r[j] * prec_v[j - end_regular];
1206 my_sums[5][0] += v[j] * prec_v[j - end_regular];
1207 my_sums[6][0] += r[j] * prec_r[j - end_regular];
1208 }
1209 for (unsigned int i = 0; i < vectorized_sums.size(); ++i)
1210 vectorized_sums[i] += my_sums[i];
1211 }
1212
1213 // Function that we use if the PreconditionerType does not implement an
1214 // apply() function, where we instead need to choose the
1215 // apply_to_subrange function
1216 template <typename U = void>
1217 std::enable_if_t<!has_apply<PreconditionerType>, U>
1218 finalize_after_convergence(const unsigned int iteration_index)
1219 {
1220 if (iteration_index % 2 == 1 || this->beta == Number())
1221 this->x.add(this->alpha, this->p);
1222 else
1223 {
1224 const unsigned int end_range = this->x.locally_owned_size();
1225
1226 Number *x = this->x.begin();
1227 Number *r = this->r.begin();
1228 Number *p = this->p.begin();
1229 const Number alpha_plus_previous_alpha_over_beta =
1230 this->alpha + this->previous_alpha / this->beta;
1231 const Number previous_alpha_over_beta =
1232 this->previous_alpha / this->beta;
1233
1234 constexpr unsigned int grain_size = 128;
1235 std::array<Number, grain_size> prec_r;
1236 for (unsigned int j = 0; j < end_range; j += grain_size)
1237 {
1238 const unsigned int length = std::min(grain_size, end_range - j);
1239 this->preconditioner.apply_to_subrange(j,
1240 j + length,
1241 r,
1242 prec_r.data());
1244 for (unsigned int i = 0; i < length; ++i)
1245 x[i] += alpha_plus_previous_alpha_over_beta * p[i] -
1246 previous_alpha_over_beta * prec_r[i];
1247
1248 x += length;
1249 r += length;
1250 p += length;
1251 }
1252 }
1253 }
1254 };
1255 } // namespace SolverCG
1256} // namespace internal
1257
1258
1259
1260template <typename VectorType>
1261template <typename MatrixType, typename PreconditionerType>
1262void
1263SolverCG<VectorType>::solve(const MatrixType &A,
1264 VectorType &x,
1265 const VectorType &b,
1266 const PreconditionerType &preconditioner)
1267{
1268 using number = typename VectorType::value_type;
1269
1271
1272 LogStream::Prefix prefix("cg");
1273
1274 // Should we build the matrix for eigenvalue computations?
1275 const bool do_eigenvalues =
1276 !condition_number_signal.empty() || !all_condition_numbers_signal.empty() ||
1277 !eigenvalues_signal.empty() || !all_eigenvalues_signal.empty();
1278
1279 // vectors used for eigenvalue computations
1280 std::vector<typename VectorType::value_type> diagonal;
1281 std::vector<typename VectorType::value_type> offdiagonal;
1282
1283 typename VectorType::value_type eigen_beta_alpha = 0;
1284
1285 int it = 0;
1286
1287 internal::SolverCG::
1288 IterationWorker<VectorType, MatrixType, PreconditionerType>
1289 worker(
1290 A, preconditioner, determine_beta_by_flexible_formula, this->memory, x);
1291
1292 worker.startup(b);
1293
1294 solver_state = this->iteration_status(0, worker.residual_norm, x);
1295 if (solver_state != SolverControl::iterate)
1296 return;
1297
1298 while (solver_state == SolverControl::iterate)
1299 {
1300 ++it;
1301
1302 worker.do_iteration(it);
1303
1304 print_vectors(it, x, worker.r, worker.p);
1305
1306 if (it > 1)
1307 {
1308 this->coefficients_signal(worker.previous_alpha, worker.beta);
1309 // set up the vectors containing the diagonal and the off diagonal
1310 // of the projected matrix.
1311 if (do_eigenvalues)
1312 {
1313 diagonal.push_back(number(1.) / worker.previous_alpha +
1314 eigen_beta_alpha);
1315 eigen_beta_alpha = worker.beta / worker.previous_alpha;
1316 offdiagonal.push_back(std::sqrt(worker.beta) /
1317 worker.previous_alpha);
1318 }
1319 compute_eigs_and_cond(diagonal,
1320 offdiagonal,
1321 all_eigenvalues_signal,
1322 all_condition_numbers_signal);
1323 }
1324
1325 solver_state = this->iteration_status(it, worker.residual_norm, x);
1326 }
1327
1328 worker.finalize_after_convergence(it);
1329
1330 compute_eigs_and_cond(diagonal,
1331 offdiagonal,
1332 eigenvalues_signal,
1333 condition_number_signal);
1334
1335 AssertThrow(solver_state == SolverControl::success,
1336 SolverControl::NoConvergence(it, worker.residual_norm));
1337}
1338
1339
1340
1341template <typename VectorType>
1342boost::signals2::connection
1344 const std::function<void(typename VectorType::value_type,
1345 typename VectorType::value_type)> &slot)
1346{
1347 return coefficients_signal.connect(slot);
1348}
1349
1350
1351
1352template <typename VectorType>
1353boost::signals2::connection
1355 const std::function<void(double)> &slot,
1356 const bool every_iteration)
1357{
1358 if (every_iteration)
1359 {
1360 return all_condition_numbers_signal.connect(slot);
1361 }
1362 else
1363 {
1364 return condition_number_signal.connect(slot);
1365 }
1366}
1367
1368
1369
1370template <typename VectorType>
1371boost::signals2::connection
1373 const std::function<void(const std::vector<double> &)> &slot,
1374 const bool every_iteration)
1375{
1376 if (every_iteration)
1377 {
1378 return all_eigenvalues_signal.connect(slot);
1379 }
1380 else
1381 {
1382 return eigenvalues_signal.connect(slot);
1383 }
1384}
1385
1386
1387
1388template <typename VectorType>
1391 const AdditionalData &)
1392 : SolverCG<VectorType>(cn, mem)
1393{
1395}
1396
1397
1398
1399template <typename VectorType>
1401 const AdditionalData &)
1402 : SolverCG<VectorType>(cn)
1403{
1405}
1406
1407
1408
1409#endif // DOXYGEN
1410
1412
1413#endif
SolverCG(SolverControl &cn, VectorMemory< VectorType > &mem, const AdditionalData &data=AdditionalData())
boost::signals2::signal< void(const std::vector< double > &)> eigenvalues_signal
Definition solver_cg.h:305
bool determine_beta_by_flexible_formula
Definition solver_cg.h:322
virtual void print_vectors(const unsigned int step, const VectorType &x, const VectorType &r, const VectorType &d) const
boost::signals2::signal< void(typename VectorType::value_type, typename VectorType::value_type)> coefficients_signal
Definition solver_cg.h:287
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::signal< void(double)> condition_number_signal
Definition solver_cg.h:293
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)
boost::signals2::signal< void(double)> all_condition_numbers_signal
Definition solver_cg.h:299
SolverCG(SolverControl &cn, const AdditionalData &data=AdditionalData())
boost::signals2::signal< void(const std::vector< double > &)> all_eigenvalues_signal
Definition solver_cg.h:312
AdditionalData additional_data
Definition solver_cg.h:280
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.
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:143
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcDivideByZero()
#define AssertThrow(cond, exc)
spacedim const Point< spacedim > & p
Definition grid_tools.h:981
@ diagonal
Matrix is diagonal.
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
T sum(const T &t, const MPI_Comm mpi_communicator)
const Iterator const std_cxx20::type_identity_t< Iterator > & end
Definition parallel.h:610
const Iterator & begin
Definition parallel.h:609
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
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)