Reference documentation for deal.II version 9.5.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\}}\)
Loading...
Searching...
No Matches
solver_cg.h
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 1998 - 2023 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
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 std::enable_if_t<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>>
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 std::enable_if_t<has_apply<PreconditionerType>, U>
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 && beta != Number())
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 if (iteration_index % 2 == 0 && beta == Number())
856 {
857 // Case where beta is zero: we cannot reconstruct p_{j-1} in the
858 // next iteration, and must hence compute x here. This can happen
859 // before termination
860 for (unsigned int j = start_range; j < end_regular; j += n_lanes)
861 {
862 VectorizedArray<Number> rj, xj, vj, pj, prec_rj;
863 rj.load(r + j);
864 vj.load(v + j);
865 xj.load(x + j);
866 pj.load(p + j);
867 xj += alpha * pj;
868 xj.store(x + j);
869 rj -= alpha * vj;
870 rj.store(r + j);
872 for (unsigned int l = 0; l < n_lanes; ++l)
873 prec_rj[l] = this->preconditioner.apply(j + l, rj[l]);
874 prec_rj.store(p + j);
876 rj.store(v + j);
877 }
878 for (unsigned int j = end_regular; j < end_range; ++j)
879 {
880 r[j] -= alpha * v[j];
881 x[j] += alpha * p[j];
882 p[j] = this->preconditioner.apply(j, r[j]);
883 v[j] = Number();
884 }
885 }
886 else
887 {
888 const Number alpha_plus_previous_alpha_over_beta =
889 alpha + this->previous_alpha / this->previous_beta;
890 const Number previous_alpha_over_beta =
891 this->previous_alpha / this->previous_beta;
892 for (unsigned int j = start_range; j < end_regular; j += n_lanes)
893 {
894 VectorizedArray<Number> rj, vj, pj, xj, prec_rj, prec_vj;
895 xj.load(x + j);
896 pj.load(p + j);
897 xj += alpha_plus_previous_alpha_over_beta * pj;
898 rj.load(r + j);
899 vj.load(v + j);
901 for (unsigned int l = 0; l < n_lanes; ++l)
902 {
903 prec_rj[l] = this->preconditioner.apply(j + l, rj[l]);
904 prec_vj[l] = this->preconditioner.apply(j + l, vj[l]);
905 }
906 xj -= previous_alpha_over_beta * prec_rj;
907 xj.store(x + j);
908 rj -= alpha * vj;
909 rj.store(r + j);
910 prec_rj -= alpha * prec_vj;
911 pj = beta * pj + prec_rj;
912 pj.store(p + j);
914 rj.store(v + j);
915 }
916 for (unsigned int j = end_regular; j < end_range; ++j)
917 {
918 x[j] += alpha_plus_previous_alpha_over_beta * p[j];
919 x[j] -= previous_alpha_over_beta *
920 this->preconditioner.apply(j, r[j]);
921 r[j] -= alpha * v[j];
922 p[j] = beta * p[j] + this->preconditioner.apply(j, r[j]);
923 v[j] = Number();
924 }
925 }
926 }
927
928 // Function that we use if the PreconditionerType implements an apply()
929 // function
930 template <typename U = void>
931 std::enable_if_t<has_apply<PreconditionerType>, U>
932 operation_after_loop(
933 const unsigned int start_range,
934 const unsigned int end_range,
935 std::array<VectorizedArray<Number>, 7> &vectorized_sums) const
936 {
937 const Number * r = this->r.begin();
938 const Number * p = this->p.begin();
939 const Number * v = this->v.begin();
940 std::array<VectorizedArray<Number>, 7> my_sums = {};
941 constexpr unsigned int n_lanes = VectorizedArray<Number>::size();
942 const unsigned int end_regular =
943 start_range + (end_range - start_range) / n_lanes * n_lanes;
944 for (unsigned int j = start_range; j < end_regular; j += n_lanes)
945 {
946 VectorizedArray<Number> pj, vj, rj, prec_vj, prec_rj;
947 pj.load(p + j);
948 vj.load(v + j);
949 rj.load(r + j);
951 for (unsigned int l = 0; l < n_lanes; ++l)
952 {
953 prec_vj[l] = this->preconditioner.apply(j + l, vj[l]);
954 prec_rj[l] = this->preconditioner.apply(j + l, rj[l]);
955 }
956 my_sums[0] += pj * vj;
957 my_sums[1] += vj * vj;
958 my_sums[2] += rj * vj;
959 my_sums[3] += rj * rj;
960 my_sums[4] += rj * prec_vj;
961 my_sums[5] += vj * prec_vj;
962 my_sums[6] += rj * prec_rj;
963 }
964 for (unsigned int j = end_regular; j < end_range; ++j)
965 {
966 const Number prec_v = this->preconditioner.apply(j, v[j]);
967 const Number prec_r = this->preconditioner.apply(j, r[j]);
968 my_sums[0][0] += p[j] * v[j];
969 my_sums[1][0] += v[j] * v[j];
970 my_sums[2][0] += r[j] * v[j];
971 my_sums[3][0] += r[j] * r[j];
972 my_sums[4][0] += r[j] * prec_v;
973 my_sums[5][0] += v[j] * prec_v;
974 my_sums[6][0] += r[j] * prec_r;
975 }
976 for (unsigned int i = 0; i < vectorized_sums.size(); ++i)
977 vectorized_sums[i] += my_sums[i];
978 }
979
980 // Function that we use if the PreconditionerType implements an apply()
981 // function
982 template <typename U = void>
983 std::enable_if_t<has_apply<PreconditionerType>, U>
984 finalize_after_convergence(const unsigned int iteration_index)
985 {
986 if (iteration_index % 2 == 1 || this->beta == Number())
987 this->x.add(this->alpha, this->p);
988 else
989 {
990 using Number = typename VectorType::value_type;
991 const unsigned int end_range = this->x.locally_owned_size();
992
993 Number *x = this->x.begin();
994 Number *r = this->r.begin();
995 Number *p = this->p.begin();
996
997 // Note that we use 'beta' here rather than 'previous_beta' in the
998 // formula above, which is because the shift in beta ->
999 // previous_beta has not been applied at this stage, allowing us
1000 // to recover the previous search direction
1001 const Number alpha_plus_previous_alpha_over_beta =
1002 this->alpha + this->previous_alpha / this->beta;
1003 const Number previous_alpha_over_beta =
1004 this->previous_alpha / this->beta;
1005
1007 for (unsigned int j = 0; j < end_range; ++j)
1008 {
1009 x[j] += alpha_plus_previous_alpha_over_beta * p[j] -
1010 previous_alpha_over_beta *
1011 this->preconditioner.apply(j, r[j]);
1012 }
1013 }
1014 }
1015
1016 // Function that we use if the PreconditionerType does not implement an
1017 // apply() function, where we instead need to choose the
1018 // apply_to_subrange function
1019 template <typename U = void>
1020 std::enable_if_t<!has_apply<PreconditionerType>, U>
1021 operation_before_loop(const unsigned int iteration_index,
1022 const unsigned int start_range,
1023 const unsigned int end_range) const
1024 {
1025 Number * x = this->x.begin() + start_range;
1026 Number * r = this->r.begin() + start_range;
1027 Number * p = this->p.begin() + start_range;
1028 Number * v = this->v.begin() + start_range;
1029 const Number alpha = this->alpha;
1030 const Number beta = this->beta;
1031 constexpr unsigned int grain_size = 128;
1032 std::array<Number, grain_size> prec_r;
1033 if (iteration_index == 1)
1034 {
1035 for (unsigned int j = start_range; j < end_range; j += grain_size)
1036 {
1037 const unsigned int length = std::min(grain_size, end_range - j);
1038 this->preconditioner.apply_to_subrange(j,
1039 j + length,
1040 r,
1041 prec_r.data());
1043 for (unsigned int i = 0; i < length; ++i)
1044 {
1045 p[i] = prec_r[i];
1046 v[i] = Number();
1047 }
1048 p += length;
1049 r += length;
1050 v += length;
1051 }
1052 }
1053 else if (iteration_index % 2 == 0 && beta != Number())
1054 {
1055 for (unsigned int j = start_range; j < end_range; j += grain_size)
1056 {
1057 const unsigned int length = std::min(grain_size, end_range - j);
1059 for (unsigned int i = 0; i < length; ++i)
1060 r[i] -= this->alpha * v[i];
1061 this->preconditioner.apply_to_subrange(j,
1062 j + length,
1063 r,
1064 prec_r.data());
1066 for (unsigned int i = 0; i < length; ++i)
1067 {
1068 p[i] = this->beta * p[i] + prec_r[i];
1069 v[i] = Number();
1070 }
1071 p += length;
1072 r += length;
1073 v += length;
1074 }
1075 }
1076 else if (iteration_index % 2 == 0 && beta == Number())
1077 {
1078 // Case where beta is zero: We cannot reconstruct p_{j-1} in the
1079 // next iteration, and must hence compute x here. This can happen
1080 // before termination
1081 for (unsigned int j = start_range; j < end_range; j += grain_size)
1082 {
1083 const unsigned int length = std::min(grain_size, end_range - j);
1085 for (unsigned int i = 0; i < length; ++i)
1086 r[i] -= this->alpha * v[i];
1087 this->preconditioner.apply_to_subrange(j,
1088 j + length,
1089 r,
1090 prec_r.data());
1092 for (unsigned int i = 0; i < length; ++i)
1093 {
1094 x[i] += this->alpha * p[i];
1095 p[i] = prec_r[i];
1096 v[i] = Number();
1097 }
1098 p += length;
1099 r += length;
1100 v += length;
1101 }
1102 }
1103 else
1104 {
1105 const Number alpha_plus_previous_alpha_over_beta =
1106 this->alpha + this->previous_alpha / this->previous_beta;
1107 const Number previous_alpha_over_beta =
1108 this->previous_alpha / this->previous_beta;
1109 for (unsigned int j = start_range; j < end_range; j += grain_size)
1110 {
1111 const unsigned int length = std::min(grain_size, end_range - j);
1112 this->preconditioner.apply_to_subrange(j,
1113 j + length,
1114 r,
1115 prec_r.data());
1117 for (unsigned int i = 0; i < length; ++i)
1118 {
1119 x[i] += alpha_plus_previous_alpha_over_beta * p[i] -
1120 previous_alpha_over_beta * prec_r[i];
1121 r[i] -= this->alpha * v[i];
1122 }
1123 this->preconditioner.apply_to_subrange(j,
1124 j + length,
1125 r,
1126 prec_r.data());
1128 for (unsigned int i = 0; i < length; ++i)
1129 {
1130 p[i] = this->beta * p[i] + prec_r[i];
1131 v[i] = Number();
1132 }
1133 p += length;
1134 r += length;
1135 v += length;
1136 x += length;
1137 }
1138 }
1139 }
1140
1141 // Function that we use if the PreconditionerType does not implement an
1142 // apply() function and where we instead need to use the
1143 // apply_to_subrange function
1144 template <typename U = void>
1145 std::enable_if_t<!has_apply<PreconditionerType>, U>
1146 operation_after_loop(
1147 const unsigned int start_range,
1148 const unsigned int end_range,
1149 std::array<VectorizedArray<Number>, 7> &vectorized_sums) const
1150 {
1151 const Number * r = this->r.begin();
1152 const Number * p = this->p.begin();
1153 const Number * v = this->v.begin();
1154 std::array<VectorizedArray<Number>, 7> my_sums = {};
1155 constexpr unsigned int grain_size = 128;
1156 Assert(grain_size % VectorizedArray<Number>::size() == 0,
1158 const unsigned int end_regular =
1159 start_range + (end_range - start_range) / grain_size * grain_size;
1160 std::array<Number, grain_size> prec_r;
1161 std::array<Number, grain_size> prec_v;
1162 for (unsigned int j = start_range; j < end_regular; j += grain_size)
1163 {
1164 this->preconditioner.apply_to_subrange(j,
1165 j + grain_size,
1166 r + j,
1167 prec_r.data());
1168 this->preconditioner.apply_to_subrange(j,
1169 j + grain_size,
1170 v + j,
1171 prec_v.data());
1172 VectorizedArray<Number> pj, vj, rj, prec_vj, prec_rj;
1173 for (unsigned int i = 0; i < grain_size;
1175 {
1176 pj.load(p + j + i);
1177 vj.load(v + j + i);
1178 rj.load(r + j + i);
1179 prec_rj.load(prec_r.data() + i);
1180 prec_vj.load(prec_v.data() + i);
1181
1182 my_sums[0] += pj * vj;
1183 my_sums[1] += vj * vj;
1184 my_sums[2] += rj * vj;
1185 my_sums[3] += rj * rj;
1186 my_sums[4] += rj * prec_vj;
1187 my_sums[5] += vj * prec_vj;
1188 my_sums[6] += rj * prec_rj;
1189 }
1190 }
1191 const unsigned int length = end_range - end_regular;
1192 AssertIndexRange(length, grain_size);
1193 this->preconditioner.apply_to_subrange(end_regular,
1194 end_regular + length,
1195 r + end_regular,
1196 prec_r.data());
1197 this->preconditioner.apply_to_subrange(end_regular,
1198 end_regular + length,
1199 v + end_regular,
1200 prec_v.data());
1201 for (unsigned int j = end_regular; j < end_range; ++j)
1202 {
1203 my_sums[0][0] += p[j] * v[j];
1204 my_sums[1][0] += v[j] * v[j];
1205 my_sums[2][0] += r[j] * v[j];
1206 my_sums[3][0] += r[j] * r[j];
1207 my_sums[4][0] += r[j] * prec_v[j - end_regular];
1208 my_sums[5][0] += v[j] * prec_v[j - end_regular];
1209 my_sums[6][0] += r[j] * prec_r[j - end_regular];
1210 }
1211 for (unsigned int i = 0; i < vectorized_sums.size(); ++i)
1212 vectorized_sums[i] += my_sums[i];
1213 }
1214
1215 // Function that we use if the PreconditionerType does not implement an
1216 // apply() function, where we instead need to choose the
1217 // apply_to_subrange function
1218 template <typename U = void>
1219 std::enable_if_t<!has_apply<PreconditionerType>, U>
1220 finalize_after_convergence(const unsigned int iteration_index)
1221 {
1222 if (iteration_index % 2 == 1 || this->beta == Number())
1223 this->x.add(this->alpha, this->p);
1224 else
1225 {
1226 const unsigned int end_range = this->x.locally_owned_size();
1227
1228 Number * x = this->x.begin();
1229 Number * r = this->r.begin();
1230 Number * p = this->p.begin();
1231 const Number alpha_plus_previous_alpha_over_beta =
1232 this->alpha + this->previous_alpha / this->beta;
1233 const Number previous_alpha_over_beta =
1234 this->previous_alpha / this->beta;
1235
1236 constexpr unsigned int grain_size = 128;
1237 std::array<Number, grain_size> prec_r;
1238 for (unsigned int j = 0; j < end_range; j += grain_size)
1239 {
1240 const unsigned int length = std::min(grain_size, end_range - j);
1241 this->preconditioner.apply_to_subrange(j,
1242 j + length,
1243 r,
1244 prec_r.data());
1246 for (unsigned int i = 0; i < length; ++i)
1247 x[i] += alpha_plus_previous_alpha_over_beta * p[i] -
1248 previous_alpha_over_beta * prec_r[i];
1249
1250 x += length;
1251 r += length;
1252 p += length;
1253 }
1254 }
1255 }
1256 };
1257 } // namespace SolverCG
1258} // namespace internal
1259
1260
1261
1262template <typename VectorType>
1263template <typename MatrixType, typename PreconditionerType>
1264void
1265SolverCG<VectorType>::solve(const MatrixType & A,
1266 VectorType & x,
1267 const VectorType & b,
1268 const PreconditionerType &preconditioner)
1269{
1270 using number = typename VectorType::value_type;
1271
1273
1274 LogStream::Prefix prefix("cg");
1275
1276 // Should we build the matrix for eigenvalue computations?
1277 const bool do_eigenvalues =
1278 !condition_number_signal.empty() || !all_condition_numbers_signal.empty() ||
1279 !eigenvalues_signal.empty() || !all_eigenvalues_signal.empty();
1280
1281 // vectors used for eigenvalue computations
1282 std::vector<typename VectorType::value_type> diagonal;
1283 std::vector<typename VectorType::value_type> offdiagonal;
1284
1285 typename VectorType::value_type eigen_beta_alpha = 0;
1286
1287 int it = 0;
1288
1289 internal::SolverCG::
1290 IterationWorker<VectorType, MatrixType, PreconditionerType>
1291 worker(
1292 A, preconditioner, determine_beta_by_flexible_formula, this->memory, x);
1293
1294 worker.startup(b);
1295
1296 solver_state = this->iteration_status(0, worker.residual_norm, x);
1297 if (solver_state != SolverControl::iterate)
1298 return;
1299
1300 while (solver_state == SolverControl::iterate)
1301 {
1302 it++;
1303
1304 worker.do_iteration(it);
1305
1306 print_vectors(it, x, worker.r, worker.p);
1307
1308 if (it > 1)
1309 {
1310 this->coefficients_signal(worker.previous_alpha, worker.beta);
1311 // set up the vectors containing the diagonal and the off diagonal
1312 // of the projected matrix.
1313 if (do_eigenvalues)
1314 {
1315 diagonal.push_back(number(1.) / worker.previous_alpha +
1316 eigen_beta_alpha);
1317 eigen_beta_alpha = worker.beta / worker.previous_alpha;
1318 offdiagonal.push_back(std::sqrt(worker.beta) /
1319 worker.previous_alpha);
1320 }
1321 compute_eigs_and_cond(diagonal,
1322 offdiagonal,
1323 all_eigenvalues_signal,
1324 all_condition_numbers_signal);
1325 }
1326
1327 solver_state = this->iteration_status(it, worker.residual_norm, x);
1328 }
1329
1330 worker.finalize_after_convergence(it);
1331
1332 compute_eigs_and_cond(diagonal,
1333 offdiagonal,
1334 eigenvalues_signal,
1335 condition_number_signal);
1336
1337 AssertThrow(solver_state == SolverControl::success,
1338 SolverControl::NoConvergence(it, worker.residual_norm));
1339}
1340
1341
1342
1343template <typename VectorType>
1344boost::signals2::connection
1346 const std::function<void(typename VectorType::value_type,
1347 typename VectorType::value_type)> &slot)
1348{
1349 return coefficients_signal.connect(slot);
1350}
1351
1352
1353
1354template <typename VectorType>
1355boost::signals2::connection
1357 const std::function<void(double)> &slot,
1358 const bool every_iteration)
1359{
1360 if (every_iteration)
1361 {
1362 return all_condition_numbers_signal.connect(slot);
1363 }
1364 else
1365 {
1366 return condition_number_signal.connect(slot);
1367 }
1368}
1369
1370
1371
1372template <typename VectorType>
1373boost::signals2::connection
1375 const std::function<void(const std::vector<double> &)> &slot,
1376 const bool every_iteration)
1377{
1378 if (every_iteration)
1379 {
1380 return all_eigenvalues_signal.connect(slot);
1381 }
1382 else
1383 {
1384 return eigenvalues_signal.connect(slot);
1385 }
1386}
1387
1388
1389
1390template <typename VectorType>
1393 const AdditionalData &)
1394 : SolverCG<VectorType>(cn, mem)
1395{
1397}
1398
1399
1400
1401template <typename VectorType>
1403 const AdditionalData &)
1404 : SolverCG<VectorType>(cn)
1405{
1407}
1408
1409
1410
1411#endif // DOXYGEN
1412
1414
1415#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())
VectorizedArrayIterator< T > begin()
void store(OtherNumber *ptr) const
void load(const OtherNumber *ptr)
#define DEAL_II_OPENMP_SIMD_PRAGMA
Definition config.h:138
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcDivideByZero()
#define AssertThrow(cond, exc)
@ diagonal
Matrix is diagonal.
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:82
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)