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