Reference documentation for deal.II version GIT 35969cdc9b 2023-12-09 01:10:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
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 
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_v<PreconditionerType, PreconditionIdentity> == false)
607  {
608  preconditioner.vmult(v, r);
609  r_dot_preconditioner_dot_r = r * v;
610  }
611  else
612  r_dot_preconditioner_dot_r = residual_norm * residual_norm;
613 
614  const VectorType &direction =
615  std::is_same_v<PreconditionerType, PreconditionIdentity> ? r : v;
616 
617  if (iteration_index > 1)
618  {
619  Assert(std::abs(previous_r_dot_preconditioner_dot_r) != 0.,
620  ExcDivideByZero());
621  beta =
622  r_dot_preconditioner_dot_r / previous_r_dot_preconditioner_dot_r;
623  if (this->flexible)
624  beta -= (r * z) / previous_r_dot_preconditioner_dot_r;
625  p.sadd(beta, 1., direction);
626  }
627  else
628  p.equ(1., direction);
629 
630  if (this->flexible)
631  z.swap(v);
632 
633  A.vmult(v, p);
634 
635  const Number p_dot_A_dot_p = p * v;
636  Assert(std::abs(p_dot_A_dot_p) != 0., ExcDivideByZero());
637 
638  this->previous_alpha = alpha;
639  alpha = r_dot_preconditioner_dot_r / p_dot_A_dot_p;
640 
641  x.add(alpha, p);
642  residual_norm = std::sqrt(std::abs(r.add_and_dot(-alpha, v, r)));
643  }
644 
645  void
646  finalize_after_convergence(const unsigned int)
647  {}
648  };
649 
650 
651  // In the following, we provide a specialization of the above
652  // IterationWorker class that picks up particular features in the matrix
653  // and preconditioners.
654 
655  // a helper type-trait that leverage SFINAE to figure out if MatrixType has
656  // ... MatrixType::vmult(VectorType &, const VectorType&,
657  // std::function<...>, std::function<...>) const
658  template <typename MatrixType, typename VectorType>
659  using vmult_functions_t = decltype(std::declval<const MatrixType>().vmult(
660  std::declval<VectorType &>(),
661  std::declval<const VectorType &>(),
662  std::declval<
663  const std::function<void(const unsigned int, const unsigned int)> &>(),
664  std::declval<const std::function<void(const unsigned int,
665  const unsigned int)> &>()));
666 
667  template <typename MatrixType, typename VectorType>
668  constexpr bool has_vmult_functions =
669  is_supported_operation<vmult_functions_t, MatrixType, VectorType>;
670 
671  // a helper type-trait that leverage SFINAE to figure out if
672  // PreconditionerType has ... PreconditionerType::apply_to_subrange(const
673  // unsigned int, const unsigned int, const Number*, Number*) const
674  template <typename PreconditionerType>
675  using apply_to_subrange_t =
676  decltype(std::declval<const PreconditionerType>()
677  .apply_to_subrange(0U, 0U, nullptr, nullptr));
678 
679  template <typename PreconditionerType>
680  constexpr bool has_apply_to_subrange =
681  is_supported_operation<apply_to_subrange_t, PreconditionerType>;
682 
683  // a helper type-trait that leverage SFINAE to figure out if
684  // PreconditionerType has ... PreconditionerType::apply(const
685  // unsigned int, const Number) const
686  template <typename PreconditionerType>
687  using apply_t =
688  decltype(std::declval<const PreconditionerType>().apply(0U, 0.0));
689 
690  template <typename PreconditionerType>
691  constexpr bool has_apply =
692  is_supported_operation<apply_t, PreconditionerType>;
693 
694 
695  // Internal function to run one iteration of the conjugate gradient solver
696  // for matrices and preconditioners that support interleaving the vector
697  // updates with the matrix-vector product.
698  template <typename VectorType,
699  typename MatrixType,
700  typename PreconditionerType>
701  struct IterationWorker<
702  VectorType,
703  MatrixType,
704  PreconditionerType,
705  std::enable_if_t<has_vmult_functions<MatrixType, VectorType> &&
706  (has_apply_to_subrange<PreconditionerType> ||
707  has_apply<PreconditionerType>)&&std::
708  is_same_v<VectorType,
709  LinearAlgebra::distributed::Vector<
710  typename VectorType::value_type,
711  MemorySpace::Host>>,
712  int>>
713  : public IterationWorkerBase<VectorType, MatrixType, PreconditionerType>
714  {
715  using Number = typename VectorType::value_type;
716 
717  Number next_r_dot_preconditioner_dot_r;
718  Number previous_beta;
719 
720  IterationWorker(const MatrixType &A,
721  const PreconditionerType &preconditioner,
722  const bool flexible,
723  VectorMemory<VectorType> &memory,
724  VectorType &x)
725  : IterationWorkerBase<VectorType, MatrixType, PreconditionerType>(
726  A,
727  preconditioner,
728  flexible,
729  memory,
730  x)
731  , next_r_dot_preconditioner_dot_r(0.)
732  , previous_beta(0.)
733  {}
734 
735  // This is the main iteration function, that will use some of the
736  // specialized functions below
737  void
738  do_iteration(const unsigned int iteration_index)
739  {
740  if (iteration_index > 1)
741  {
742  previous_beta = this->beta;
743  this->beta = next_r_dot_preconditioner_dot_r /
744  this->r_dot_preconditioner_dot_r;
745  }
746 
747  std::array<VectorizedArray<Number>, 7> vectorized_sums = {};
748 
749  this->A.vmult(
750  this->v,
751  this->p,
752  [&](const unsigned int begin, const unsigned int end) {
753  operation_before_loop(iteration_index, begin, end);
754  },
755  [&](const unsigned int begin, const unsigned int end) {
756  operation_after_loop(begin, end, vectorized_sums);
757  });
758 
759  std::array<Number, 7> scalar_sums;
760  for (unsigned int i = 0; i < 7; ++i)
761  scalar_sums[i] = vectorized_sums[i][0];
762  for (unsigned int l = 1; l < VectorizedArray<Number>::size(); ++l)
763  for (unsigned int i = 0; i < 7; ++i)
764  scalar_sums[i] += vectorized_sums[i][l];
765 
766  Utilities::MPI::sum(::ArrayView<const Number>(scalar_sums.data(),
767  7),
768  this->r.get_mpi_communicator(),
769  ::ArrayView<Number>(scalar_sums.data(), 7));
770 
771  this->r_dot_preconditioner_dot_r = scalar_sums[6];
772 
773  const Number p_dot_A_dot_p = scalar_sums[0];
774  Assert(std::abs(p_dot_A_dot_p) != 0., ExcDivideByZero());
775 
776  this->previous_alpha = this->alpha;
777  this->alpha = this->r_dot_preconditioner_dot_r / p_dot_A_dot_p;
778 
779  // Round-off errors near zero might yield negative values, so take
780  // the absolute value in the next two formulas
781  this->residual_norm = std::sqrt(std::abs(
782  scalar_sums[3] +
783  this->alpha * (-2. * scalar_sums[2] + this->alpha * scalar_sums[1])));
784 
785  next_r_dot_preconditioner_dot_r = std::abs(
786  this->r_dot_preconditioner_dot_r +
787  this->alpha * (-2. * scalar_sums[4] + this->alpha * scalar_sums[5]));
788  }
789 
790  // Function that we use if the PreconditionerType implements an apply()
791  // function
792  template <typename U = void>
793  std::enable_if_t<has_apply<PreconditionerType>, U>
794  operation_before_loop(const unsigned int iteration_index,
795  const unsigned int start_range,
796  const unsigned int end_range) const
797  {
798  Number *x = this->x.begin();
799  Number *r = this->r.begin();
800  Number *p = this->p.begin();
801  Number *v = this->v.begin();
802  const Number alpha = this->alpha;
803  const Number beta = this->beta;
804  constexpr unsigned int n_lanes = VectorizedArray<Number>::size();
805  const unsigned int end_regular =
806  start_range + (end_range - start_range) / n_lanes * n_lanes;
807  if (iteration_index == 1)
808  {
809  // Vectorize by hand since compilers are often pretty bad at
810  // doing these steps automatically even with
811  // DEAL_II_OPENMP_SIMD_PRAGMA
812  for (unsigned int j = start_range; j < end_regular; j += n_lanes)
813  {
815  rj.load(r + j);
817  for (unsigned int l = 0; l < n_lanes; ++l)
818  pj[l] = this->preconditioner.apply(j + l, rj[l]);
819  pj.store(p + j);
821  rj.store(v + j);
822  }
823  for (unsigned int j = end_regular; j < end_range; ++j)
824  {
825  p[j] = this->preconditioner.apply(j, r[j]);
826  v[j] = Number();
827  }
828  }
829  else if (iteration_index % 2 == 0 && beta != Number())
830  {
831  for (unsigned int j = start_range; j < end_regular; j += n_lanes)
832  {
833  VectorizedArray<Number> rj, vj, pj, prec_rj;
834  rj.load(r + j);
835  vj.load(v + j);
836  rj -= alpha * vj;
837  rj.store(r + j);
839  for (unsigned int l = 0; l < n_lanes; ++l)
840  prec_rj[l] = this->preconditioner.apply(j + l, rj[l]);
841  pj.load(p + j);
842  pj = beta * pj + prec_rj;
843  pj.store(p + j);
845  rj.store(v + j);
846  }
847  for (unsigned int j = end_regular; j < end_range; ++j)
848  {
849  r[j] -= alpha * v[j];
850  p[j] = beta * p[j] + this->preconditioner.apply(j, r[j]);
851  v[j] = Number();
852  }
853  }
854  else if (iteration_index % 2 == 0 && beta == Number())
855  {
856  // Case where beta is zero: we cannot reconstruct p_{j-1} in the
857  // next iteration, and must hence compute x here. This can happen
858  // before termination
859  for (unsigned int j = start_range; j < end_regular; j += n_lanes)
860  {
861  VectorizedArray<Number> rj, xj, vj, pj, prec_rj;
862  rj.load(r + j);
863  vj.load(v + j);
864  xj.load(x + j);
865  pj.load(p + j);
866  xj += alpha * pj;
867  xj.store(x + j);
868  rj -= alpha * vj;
869  rj.store(r + j);
871  for (unsigned int l = 0; l < n_lanes; ++l)
872  prec_rj[l] = this->preconditioner.apply(j + l, rj[l]);
873  prec_rj.store(p + j);
875  rj.store(v + j);
876  }
877  for (unsigned int j = end_regular; j < end_range; ++j)
878  {
879  r[j] -= alpha * v[j];
880  x[j] += alpha * p[j];
881  p[j] = this->preconditioner.apply(j, r[j]);
882  v[j] = Number();
883  }
884  }
885  else
886  {
887  const Number alpha_plus_previous_alpha_over_beta =
888  alpha + this->previous_alpha / this->previous_beta;
889  const Number previous_alpha_over_beta =
890  this->previous_alpha / this->previous_beta;
891  for (unsigned int j = start_range; j < end_regular; j += n_lanes)
892  {
893  VectorizedArray<Number> rj, vj, pj, xj, prec_rj, prec_vj;
894  xj.load(x + j);
895  pj.load(p + j);
896  xj += alpha_plus_previous_alpha_over_beta * pj;
897  rj.load(r + j);
898  vj.load(v + j);
900  for (unsigned int l = 0; l < n_lanes; ++l)
901  {
902  prec_rj[l] = this->preconditioner.apply(j + l, rj[l]);
903  prec_vj[l] = this->preconditioner.apply(j + l, vj[l]);
904  }
905  xj -= previous_alpha_over_beta * prec_rj;
906  xj.store(x + j);
907  rj -= alpha * vj;
908  rj.store(r + j);
909  prec_rj -= alpha * prec_vj;
910  pj = beta * pj + prec_rj;
911  pj.store(p + j);
913  rj.store(v + j);
914  }
915  for (unsigned int j = end_regular; j < end_range; ++j)
916  {
917  x[j] += alpha_plus_previous_alpha_over_beta * p[j];
918  x[j] -= previous_alpha_over_beta *
919  this->preconditioner.apply(j, r[j]);
920  r[j] -= alpha * v[j];
921  p[j] = beta * p[j] + this->preconditioner.apply(j, r[j]);
922  v[j] = Number();
923  }
924  }
925  }
926 
927  // Function that we use if the PreconditionerType implements an apply()
928  // function
929  template <typename U = void>
930  std::enable_if_t<has_apply<PreconditionerType>, U>
931  operation_after_loop(
932  const unsigned int start_range,
933  const unsigned int end_range,
934  std::array<VectorizedArray<Number>, 7> &vectorized_sums) const
935  {
936  const Number *r = this->r.begin();
937  const Number *p = this->p.begin();
938  const Number *v = this->v.begin();
939  std::array<VectorizedArray<Number>, 7> my_sums = {};
940  constexpr unsigned int n_lanes = VectorizedArray<Number>::size();
941  const unsigned int end_regular =
942  start_range + (end_range - start_range) / n_lanes * n_lanes;
943  for (unsigned int j = start_range; j < end_regular; j += n_lanes)
944  {
945  VectorizedArray<Number> pj, vj, rj, prec_vj, prec_rj;
946  pj.load(p + j);
947  vj.load(v + j);
948  rj.load(r + j);
950  for (unsigned int l = 0; l < n_lanes; ++l)
951  {
952  prec_vj[l] = this->preconditioner.apply(j + l, vj[l]);
953  prec_rj[l] = this->preconditioner.apply(j + l, rj[l]);
954  }
955  my_sums[0] += pj * vj;
956  my_sums[1] += vj * vj;
957  my_sums[2] += rj * vj;
958  my_sums[3] += rj * rj;
959  my_sums[4] += rj * prec_vj;
960  my_sums[5] += vj * prec_vj;
961  my_sums[6] += rj * prec_rj;
962  }
963  for (unsigned int j = end_regular; j < end_range; ++j)
964  {
965  const Number prec_v = this->preconditioner.apply(j, v[j]);
966  const Number prec_r = this->preconditioner.apply(j, r[j]);
967  my_sums[0][0] += p[j] * v[j];
968  my_sums[1][0] += v[j] * v[j];
969  my_sums[2][0] += r[j] * v[j];
970  my_sums[3][0] += r[j] * r[j];
971  my_sums[4][0] += r[j] * prec_v;
972  my_sums[5][0] += v[j] * prec_v;
973  my_sums[6][0] += r[j] * prec_r;
974  }
975  for (unsigned int i = 0; i < vectorized_sums.size(); ++i)
976  vectorized_sums[i] += my_sums[i];
977  }
978 
979  // Function that we use if the PreconditionerType implements an apply()
980  // function
981  template <typename U = void>
982  std::enable_if_t<has_apply<PreconditionerType>, U>
983  finalize_after_convergence(const unsigned int iteration_index)
984  {
985  if (iteration_index % 2 == 1 || this->beta == Number())
986  this->x.add(this->alpha, this->p);
987  else
988  {
989  using Number = typename VectorType::value_type;
990  const unsigned int end_range = this->x.locally_owned_size();
991 
992  Number *x = this->x.begin();
993  Number *r = this->r.begin();
994  Number *p = this->p.begin();
995 
996  // Note that we use 'beta' here rather than 'previous_beta' in the
997  // formula above, which is because the shift in beta ->
998  // previous_beta has not been applied at this stage, allowing us
999  // to recover the previous search direction
1000  const Number alpha_plus_previous_alpha_over_beta =
1001  this->alpha + this->previous_alpha / this->beta;
1002  const Number previous_alpha_over_beta =
1003  this->previous_alpha / this->beta;
1004 
1006  for (unsigned int j = 0; j < end_range; ++j)
1007  {
1008  x[j] += alpha_plus_previous_alpha_over_beta * p[j] -
1009  previous_alpha_over_beta *
1010  this->preconditioner.apply(j, r[j]);
1011  }
1012  }
1013  }
1014 
1015  // Function that we use if the PreconditionerType does not implement an
1016  // apply() function, where we instead need to choose the
1017  // apply_to_subrange function
1018  template <typename U = void>
1019  std::enable_if_t<!has_apply<PreconditionerType>, U>
1020  operation_before_loop(const unsigned int iteration_index,
1021  const unsigned int start_range,
1022  const unsigned int end_range) const
1023  {
1024  Number *x = this->x.begin() + start_range;
1025  Number *r = this->r.begin() + start_range;
1026  Number *p = this->p.begin() + start_range;
1027  Number *v = this->v.begin() + start_range;
1028  const Number alpha = this->alpha;
1029  const Number beta = this->beta;
1030  constexpr unsigned int grain_size = 128;
1031  std::array<Number, grain_size> prec_r;
1032  if (iteration_index == 1)
1033  {
1034  for (unsigned int j = start_range; j < end_range; j += grain_size)
1035  {
1036  const unsigned int length = std::min(grain_size, end_range - j);
1037  this->preconditioner.apply_to_subrange(j,
1038  j + length,
1039  r,
1040  prec_r.data());
1042  for (unsigned int i = 0; i < length; ++i)
1043  {
1044  p[i] = prec_r[i];
1045  v[i] = Number();
1046  }
1047  p += length;
1048  r += length;
1049  v += length;
1050  }
1051  }
1052  else if (iteration_index % 2 == 0 && beta != Number())
1053  {
1054  for (unsigned int j = start_range; j < end_range; j += grain_size)
1055  {
1056  const unsigned int length = std::min(grain_size, end_range - j);
1058  for (unsigned int i = 0; i < length; ++i)
1059  r[i] -= this->alpha * v[i];
1060  this->preconditioner.apply_to_subrange(j,
1061  j + length,
1062  r,
1063  prec_r.data());
1065  for (unsigned int i = 0; i < length; ++i)
1066  {
1067  p[i] = this->beta * p[i] + prec_r[i];
1068  v[i] = Number();
1069  }
1070  p += length;
1071  r += length;
1072  v += length;
1073  }
1074  }
1075  else if (iteration_index % 2 == 0 && beta == Number())
1076  {
1077  // Case where beta is zero: We cannot reconstruct p_{j-1} in the
1078  // next iteration, and must hence compute x here. This can happen
1079  // before termination
1080  for (unsigned int j = start_range; j < end_range; j += grain_size)
1081  {
1082  const unsigned int length = std::min(grain_size, end_range - j);
1084  for (unsigned int i = 0; i < length; ++i)
1085  r[i] -= this->alpha * v[i];
1086  this->preconditioner.apply_to_subrange(j,
1087  j + length,
1088  r,
1089  prec_r.data());
1091  for (unsigned int i = 0; i < length; ++i)
1092  {
1093  x[i] += this->alpha * p[i];
1094  p[i] = prec_r[i];
1095  v[i] = Number();
1096  }
1097  p += length;
1098  r += length;
1099  v += length;
1100  }
1101  }
1102  else
1103  {
1104  const Number alpha_plus_previous_alpha_over_beta =
1105  this->alpha + this->previous_alpha / this->previous_beta;
1106  const Number previous_alpha_over_beta =
1107  this->previous_alpha / this->previous_beta;
1108  for (unsigned int j = start_range; j < end_range; j += grain_size)
1109  {
1110  const unsigned int length = std::min(grain_size, end_range - j);
1111  this->preconditioner.apply_to_subrange(j,
1112  j + length,
1113  r,
1114  prec_r.data());
1116  for (unsigned int i = 0; i < length; ++i)
1117  {
1118  x[i] += alpha_plus_previous_alpha_over_beta * p[i] -
1119  previous_alpha_over_beta * prec_r[i];
1120  r[i] -= this->alpha * v[i];
1121  }
1122  this->preconditioner.apply_to_subrange(j,
1123  j + length,
1124  r,
1125  prec_r.data());
1127  for (unsigned int i = 0; i < length; ++i)
1128  {
1129  p[i] = this->beta * p[i] + prec_r[i];
1130  v[i] = Number();
1131  }
1132  p += length;
1133  r += length;
1134  v += length;
1135  x += length;
1136  }
1137  }
1138  }
1139 
1140  // Function that we use if the PreconditionerType does not implement an
1141  // apply() function and where we instead need to use the
1142  // apply_to_subrange function
1143  template <typename U = void>
1144  std::enable_if_t<!has_apply<PreconditionerType>, U>
1145  operation_after_loop(
1146  const unsigned int start_range,
1147  const unsigned int end_range,
1148  std::array<VectorizedArray<Number>, 7> &vectorized_sums) const
1149  {
1150  const Number *r = this->r.begin();
1151  const Number *p = this->p.begin();
1152  const Number *v = this->v.begin();
1153  std::array<VectorizedArray<Number>, 7> my_sums = {};
1154  constexpr unsigned int grain_size = 128;
1155  Assert(grain_size % VectorizedArray<Number>::size() == 0,
1156  ExcNotImplemented());
1157  const unsigned int end_regular =
1158  start_range + (end_range - start_range) / grain_size * grain_size;
1159  std::array<Number, grain_size> prec_r;
1160  std::array<Number, grain_size> prec_v;
1161  for (unsigned int j = start_range; j < end_regular; j += grain_size)
1162  {
1163  this->preconditioner.apply_to_subrange(j,
1164  j + grain_size,
1165  r + j,
1166  prec_r.data());
1167  this->preconditioner.apply_to_subrange(j,
1168  j + grain_size,
1169  v + j,
1170  prec_v.data());
1171  VectorizedArray<Number> pj, vj, rj, prec_vj, prec_rj;
1172  for (unsigned int i = 0; i < grain_size;
1174  {
1175  pj.load(p + j + i);
1176  vj.load(v + j + i);
1177  rj.load(r + j + i);
1178  prec_rj.load(prec_r.data() + i);
1179  prec_vj.load(prec_v.data() + i);
1180 
1181  my_sums[0] += pj * vj;
1182  my_sums[1] += vj * vj;
1183  my_sums[2] += rj * vj;
1184  my_sums[3] += rj * rj;
1185  my_sums[4] += rj * prec_vj;
1186  my_sums[5] += vj * prec_vj;
1187  my_sums[6] += rj * prec_rj;
1188  }
1189  }
1190  const unsigned int length = end_range - end_regular;
1191  AssertIndexRange(length, grain_size);
1192  this->preconditioner.apply_to_subrange(end_regular,
1193  end_regular + length,
1194  r + end_regular,
1195  prec_r.data());
1196  this->preconditioner.apply_to_subrange(end_regular,
1197  end_regular + length,
1198  v + end_regular,
1199  prec_v.data());
1200  for (unsigned int j = end_regular; j < end_range; ++j)
1201  {
1202  my_sums[0][0] += p[j] * v[j];
1203  my_sums[1][0] += v[j] * v[j];
1204  my_sums[2][0] += r[j] * v[j];
1205  my_sums[3][0] += r[j] * r[j];
1206  my_sums[4][0] += r[j] * prec_v[j - end_regular];
1207  my_sums[5][0] += v[j] * prec_v[j - end_regular];
1208  my_sums[6][0] += r[j] * prec_r[j - end_regular];
1209  }
1210  for (unsigned int i = 0; i < vectorized_sums.size(); ++i)
1211  vectorized_sums[i] += my_sums[i];
1212  }
1213 
1214  // Function that we use if the PreconditionerType does not implement an
1215  // apply() function, where we instead need to choose the
1216  // apply_to_subrange function
1217  template <typename U = void>
1218  std::enable_if_t<!has_apply<PreconditionerType>, U>
1219  finalize_after_convergence(const unsigned int iteration_index)
1220  {
1221  if (iteration_index % 2 == 1 || this->beta == Number())
1222  this->x.add(this->alpha, this->p);
1223  else
1224  {
1225  const unsigned int end_range = this->x.locally_owned_size();
1226 
1227  Number *x = this->x.begin();
1228  Number *r = this->r.begin();
1229  Number *p = this->p.begin();
1230  const Number alpha_plus_previous_alpha_over_beta =
1231  this->alpha + this->previous_alpha / this->beta;
1232  const Number previous_alpha_over_beta =
1233  this->previous_alpha / this->beta;
1234 
1235  constexpr unsigned int grain_size = 128;
1236  std::array<Number, grain_size> prec_r;
1237  for (unsigned int j = 0; j < end_range; j += grain_size)
1238  {
1239  const unsigned int length = std::min(grain_size, end_range - j);
1240  this->preconditioner.apply_to_subrange(j,
1241  j + length,
1242  r,
1243  prec_r.data());
1245  for (unsigned int i = 0; i < length; ++i)
1246  x[i] += alpha_plus_previous_alpha_over_beta * p[i] -
1247  previous_alpha_over_beta * prec_r[i];
1248 
1249  x += length;
1250  r += length;
1251  p += length;
1252  }
1253  }
1254  }
1255  };
1256  } // namespace SolverCG
1257 } // namespace internal
1258 
1259 
1260 
1261 template <typename VectorType>
1262 template <typename MatrixType, typename PreconditionerType>
1263 void
1264 SolverCG<VectorType>::solve(const MatrixType &A,
1265  VectorType &x,
1266  const VectorType &b,
1267  const PreconditionerType &preconditioner)
1268 {
1269  using number = typename VectorType::value_type;
1270 
1272 
1273  LogStream::Prefix prefix("cg");
1274 
1275  // Should we build the matrix for eigenvalue computations?
1276  const bool do_eigenvalues =
1277  !condition_number_signal.empty() || !all_condition_numbers_signal.empty() ||
1278  !eigenvalues_signal.empty() || !all_eigenvalues_signal.empty();
1279 
1280  // vectors used for eigenvalue computations
1281  std::vector<typename VectorType::value_type> diagonal;
1282  std::vector<typename VectorType::value_type> offdiagonal;
1283 
1284  typename VectorType::value_type eigen_beta_alpha = 0;
1285 
1286  int it = 0;
1287 
1288  internal::SolverCG::
1289  IterationWorker<VectorType, MatrixType, PreconditionerType>
1290  worker(
1291  A, preconditioner, determine_beta_by_flexible_formula, this->memory, x);
1292 
1293  worker.startup(b);
1294 
1295  solver_state = this->iteration_status(0, worker.residual_norm, x);
1296  if (solver_state != SolverControl::iterate)
1297  return;
1298 
1299  while (solver_state == SolverControl::iterate)
1300  {
1301  it++;
1302 
1303  worker.do_iteration(it);
1304 
1305  print_vectors(it, x, worker.r, worker.p);
1306 
1307  if (it > 1)
1308  {
1309  this->coefficients_signal(worker.previous_alpha, worker.beta);
1310  // set up the vectors containing the diagonal and the off diagonal
1311  // of the projected matrix.
1312  if (do_eigenvalues)
1313  {
1314  diagonal.push_back(number(1.) / worker.previous_alpha +
1315  eigen_beta_alpha);
1316  eigen_beta_alpha = worker.beta / worker.previous_alpha;
1317  offdiagonal.push_back(std::sqrt(worker.beta) /
1318  worker.previous_alpha);
1319  }
1320  compute_eigs_and_cond(diagonal,
1321  offdiagonal,
1322  all_eigenvalues_signal,
1323  all_condition_numbers_signal);
1324  }
1325 
1326  solver_state = this->iteration_status(it, worker.residual_norm, x);
1327  }
1328 
1329  worker.finalize_after_convergence(it);
1330 
1331  compute_eigs_and_cond(diagonal,
1332  offdiagonal,
1333  eigenvalues_signal,
1334  condition_number_signal);
1335 
1336  AssertThrow(solver_state == SolverControl::success,
1337  SolverControl::NoConvergence(it, worker.residual_norm));
1338 }
1339 
1340 
1341 
1342 template <typename VectorType>
1343 boost::signals2::connection
1345  const std::function<void(typename VectorType::value_type,
1346  typename VectorType::value_type)> &slot)
1347 {
1348  return coefficients_signal.connect(slot);
1349 }
1350 
1351 
1352 
1353 template <typename VectorType>
1354 boost::signals2::connection
1356  const std::function<void(double)> &slot,
1357  const bool every_iteration)
1358 {
1359  if (every_iteration)
1360  {
1361  return all_condition_numbers_signal.connect(slot);
1362  }
1363  else
1364  {
1365  return condition_number_signal.connect(slot);
1366  }
1367 }
1368 
1369 
1370 
1371 template <typename VectorType>
1372 boost::signals2::connection
1374  const std::function<void(const std::vector<double> &)> &slot,
1375  const bool every_iteration)
1376 {
1377  if (every_iteration)
1378  {
1379  return all_eigenvalues_signal.connect(slot);
1380  }
1381  else
1382  {
1383  return eigenvalues_signal.connect(slot);
1384  }
1385 }
1386 
1387 
1388 
1389 template <typename VectorType>
1392  const AdditionalData &)
1393  : SolverCG<VectorType>(cn, mem)
1394 {
1396 }
1397 
1398 
1399 
1400 template <typename VectorType>
1402  const AdditionalData &)
1403  : SolverCG<VectorType>(cn)
1404 {
1406 }
1407 
1408 
1409 
1410 #endif // DOXYGEN
1411 
1413 
1414 #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:110
void store(OtherNumber *ptr) const
void load(const OtherNumber *ptr)
#define DEAL_II_OPENMP_SIMD_PRAGMA
Definition: config.h:144
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
static ::ExceptionBase & ExcDivideByZero()
#define Assert(cond, exc)
Definition: exceptions.h:1631
static ::ExceptionBase & ExcNotImplemented()
#define AssertIndexRange(index, range)
Definition: exceptions.h:1888
#define AssertThrow(cond, exc)
Definition: exceptions.h:1732
void apply(const Kokkos::TeamPolicy< MemorySpace::Default::kokkos_space::execution_space >::member_type &team_member, const Kokkos::View< Number *, MemorySpace::Default::kokkos_space > shape_data, const ViewTypeIn in, ViewTypeOut out)
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)
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 > &)