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