Reference documentation for deal.II version GIT 85919f3e70 2023-05-28 07:10:01+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\}}\)
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 
23 #include <deal.II/base/logstream.h>
26 
27 #include <deal.II/lac/solver.h>
30 
31 #include <cmath>
32 
34 
35 // forward declaration
36 #ifndef DOXYGEN
38 namespace LinearAlgebra
39 {
40  namespace distributed
41  {
42  template <typename, typename>
43  class Vector;
44  }
45 } // namespace LinearAlgebra
46 #endif
47 
48 
176 template <typename VectorType = Vector<double>>
177 class SolverCG : public SolverBase<VectorType>
178 {
179 public:
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 
253 protected:
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 
353 template <typename VectorType = Vector<double>>
354 class SolverFlexibleCG : public SolverCG<VectorType>
355 {
356 public:
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 
394 template <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 
405 template <typename VectorType>
406 SolverCG<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 
414 template <typename VectorType>
415 void
416 SolverCG<VectorType>::print_vectors(const unsigned int,
417  const VectorType &,
418  const VectorType &,
419  const VectorType &) const
420 {}
421 
422 
423 
424 template <typename VectorType>
425 inline 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 
471 namespace 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,
518  VectorMemory<VectorType> &memory,
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,
581  VectorMemory<VectorType> &memory,
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.,
621  ExcDivideByZero());
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,
724  VectorMemory<VectorType> &memory,
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 
767  Utilities::MPI::sum(::ArrayView<const Number>(scalar_sums.data(),
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,
1157  ExcNotImplemented());
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 
1262 template <typename VectorType>
1263 template <typename MatrixType, typename PreconditionerType>
1264 void
1265 SolverCG<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 
1343 template <typename VectorType>
1344 boost::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 
1354 template <typename VectorType>
1355 boost::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 
1372 template <typename VectorType>
1373 boost::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 
1390 template <typename VectorType>
1393  const AdditionalData &)
1394  : SolverCG<VectorType>(cn, mem)
1395 {
1397 }
1398 
1399 
1400 
1401 template <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
types::global_dof_index size_type
Definition: solver_cg.h:183
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.
types::global_dof_index size_type
Definition: solver_cg.h:360
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(OtherNumber *ptr) const
void load(const OtherNumber *ptr)
#define DEAL_II_OPENMP_SIMD_PRAGMA
Definition: config.h:140
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:474
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:475
static ::ExceptionBase & ExcDivideByZero()
#define Assert(cond, exc)
Definition: exceptions.h:1614
static ::ExceptionBase & ExcNotImplemented()
#define AssertIndexRange(index, range)
Definition: exceptions.h:1855
#define AssertThrow(cond, exc)
Definition: exceptions.h:1703
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
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
VectorType::value_type * begin(VectorType &V)
VectorType::value_type * end(VectorType &V)
T sum(const T &t, const MPI_Comm mpi_communicator)
auto apply(F &&fn, Tuple &&t) -> decltype(apply_impl(std::forward< F >(fn), std::forward< Tuple >(t), std::make_index_sequence< std::tuple_size< typename std::remove_reference< Tuple >::type >::value >()))
Definition: tuple.h:36
unsigned int global_dof_index
Definition: types.h:82
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)