Reference documentation for deal.II version GIT 5dcc62ab57 2022-07-04 21:05: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 - 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 
51 
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  typename std::enable_if<has_vmult_functions<MatrixType, VectorType> &&
707  (has_apply_to_subrange<PreconditionerType> ||
708  has_apply<PreconditionerType>)&&std::
709  is_same<VectorType,
710  LinearAlgebra::distributed::Vector<
711  typename VectorType::value_type,
712  MemorySpace::Host>>::value,
713  int>::type>
714  : public IterationWorkerBase<VectorType, MatrixType, PreconditionerType>
715  {
716  using Number = typename VectorType::value_type;
717 
718  Number next_r_dot_preconditioner_dot_r;
719  Number previous_beta;
720 
721  IterationWorker(const MatrixType & A,
722  const PreconditionerType &preconditioner,
723  const bool flexible,
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  typename std::enable_if<has_apply<PreconditionerType>, U>::type
795  operation_before_loop(const unsigned int iteration_index,
796  const unsigned int start_range,
797  const unsigned int end_range) const
798  {
799  Number * x = this->x.begin();
800  Number * r = this->r.begin();
801  Number * p = this->p.begin();
802  Number * v = this->v.begin();
803  const Number alpha = this->alpha;
804  const Number beta = this->beta;
805  constexpr unsigned int n_lanes = VectorizedArray<Number>::size();
806  const unsigned int end_regular =
807  start_range + (end_range - start_range) / n_lanes * n_lanes;
808  if (iteration_index == 1)
809  {
810  // Vectorize by hand since compilers are often pretty bad at
811  // doing these steps automatically even with
812  // DEAL_II_OPENMP_SIMD_PRAGMA
813  for (unsigned int j = start_range; j < end_regular; j += n_lanes)
814  {
816  rj.load(r + j);
818  for (unsigned int l = 0; l < n_lanes; ++l)
819  pj[l] = this->preconditioner.apply(j + l, rj[l]);
820  pj.store(p + j);
822  rj.store(v + j);
823  }
824  for (unsigned int j = end_regular; j < end_range; ++j)
825  {
826  p[j] = this->preconditioner.apply(j, r[j]);
827  v[j] = Number();
828  }
829  }
830  else if (iteration_index % 2 == 0)
831  {
832  for (unsigned int j = start_range; j < end_regular; j += n_lanes)
833  {
834  VectorizedArray<Number> rj, vj, pj, prec_rj;
835  rj.load(r + j);
836  vj.load(v + j);
837  rj -= alpha * vj;
838  rj.store(r + j);
840  for (unsigned int l = 0; l < n_lanes; ++l)
841  prec_rj[l] = this->preconditioner.apply(j + l, rj[l]);
842  pj.load(p + j);
843  pj = beta * pj + prec_rj;
844  pj.store(p + j);
846  rj.store(v + j);
847  }
848  for (unsigned int j = end_regular; j < end_range; ++j)
849  {
850  r[j] -= alpha * v[j];
851  p[j] = beta * p[j] + this->preconditioner.apply(j, r[j]);
852  v[j] = Number();
853  }
854  }
855  else
856  {
857  const Number alpha_plus_previous_alpha_over_beta =
858  alpha + this->previous_alpha / this->previous_beta;
859  const Number previous_alpha_over_beta =
860  this->previous_alpha / this->previous_beta;
861  for (unsigned int j = start_range; j < end_regular; j += n_lanes)
862  {
863  VectorizedArray<Number> rj, vj, pj, xj, prec_rj, prec_vj;
864  xj.load(x + j);
865  pj.load(p + j);
866  xj += alpha_plus_previous_alpha_over_beta * pj;
867  rj.load(r + j);
868  vj.load(v + j);
870  for (unsigned int l = 0; l < n_lanes; ++l)
871  {
872  prec_rj[l] = this->preconditioner.apply(j + l, rj[l]);
873  prec_vj[l] = this->preconditioner.apply(j + l, vj[l]);
874  }
875  xj -= previous_alpha_over_beta * prec_rj;
876  xj.store(x + j);
877  rj -= alpha * vj;
878  rj.store(r + j);
879  prec_rj -= alpha * prec_vj;
880  pj = beta * pj + prec_rj;
881  pj.store(p + j);
883  rj.store(v + j);
884  }
885  for (unsigned int j = end_regular; j < end_range; ++j)
886  {
887  x[j] += alpha_plus_previous_alpha_over_beta * p[j];
888  x[j] -= previous_alpha_over_beta *
889  this->preconditioner.apply(j, r[j]);
890  r[j] -= alpha * v[j];
891  p[j] = beta * p[j] + this->preconditioner.apply(j, r[j]);
892  v[j] = Number();
893  }
894  }
895  }
896 
897  // Function that we use if the PreconditionerType implements an apply()
898  // function
899  template <typename U = void>
900  typename std::enable_if<has_apply<PreconditionerType>, U>::type
901  operation_after_loop(
902  const unsigned int start_range,
903  const unsigned int end_range,
904  std::array<VectorizedArray<Number>, 7> &vectorized_sums) const
905  {
906  const Number * r = this->r.begin();
907  const Number * p = this->p.begin();
908  const Number * v = this->v.begin();
909  std::array<VectorizedArray<Number>, 7> my_sums = {};
910  constexpr unsigned int n_lanes = VectorizedArray<Number>::size();
911  const unsigned int end_regular =
912  start_range + (end_range - start_range) / n_lanes * n_lanes;
913  for (unsigned int j = start_range; j < end_regular; j += n_lanes)
914  {
915  VectorizedArray<Number> pj, vj, rj, prec_vj, prec_rj;
916  pj.load(p + j);
917  vj.load(v + j);
918  rj.load(r + j);
920  for (unsigned int l = 0; l < n_lanes; ++l)
921  {
922  prec_vj[l] = this->preconditioner.apply(j + l, vj[l]);
923  prec_rj[l] = this->preconditioner.apply(j + l, rj[l]);
924  }
925  my_sums[0] += pj * vj;
926  my_sums[1] += vj * vj;
927  my_sums[2] += rj * vj;
928  my_sums[3] += rj * rj;
929  my_sums[4] += rj * prec_vj;
930  my_sums[5] += vj * prec_vj;
931  my_sums[6] += rj * prec_rj;
932  }
933  for (unsigned int j = end_regular; j < end_range; ++j)
934  {
935  const Number prec_v = this->preconditioner.apply(j, v[j]);
936  const Number prec_r = this->preconditioner.apply(j, r[j]);
937  my_sums[0][0] += p[j] * v[j];
938  my_sums[1][0] += v[j] * v[j];
939  my_sums[2][0] += r[j] * v[j];
940  my_sums[3][0] += r[j] * r[j];
941  my_sums[4][0] += r[j] * prec_v;
942  my_sums[5][0] += v[j] * prec_v;
943  my_sums[6][0] += r[j] * prec_r;
944  }
945  for (unsigned int i = 0; i < vectorized_sums.size(); ++i)
946  vectorized_sums[i] += my_sums[i];
947  }
948 
949  // Function that we use if the PreconditionerType implements an apply()
950  // function
951  template <typename U = void>
952  typename std::enable_if<has_apply<PreconditionerType>, U>::type
953  finalize_after_convergence(const unsigned int iteration_index)
954  {
955  if (iteration_index % 2 == 1 || iteration_index == 2)
956  this->x.add(this->alpha, this->p);
957  else
958  {
959  using Number = typename VectorType::value_type;
960  const unsigned int end_range = this->x.locally_owned_size();
961 
962  Number *x = this->x.begin();
963  Number *r = this->r.begin();
964  Number *p = this->p.begin();
965 
966  // Note that we use 'beta' here rather than 'previous_beta' in the
967  // formula above, which is because the shift in beta ->
968  // previous_beta has not been applied at this stage, allowing us
969  // to recover the previous search direction
970  const Number alpha_plus_previous_alpha_over_beta =
971  this->alpha + this->previous_alpha / this->previous_beta;
972  const Number previous_alpha_over_beta =
973  this->previous_alpha / this->previous_beta;
974 
976  for (unsigned int j = 0; j < end_range; ++j)
977  {
978  x[j] += alpha_plus_previous_alpha_over_beta * p[j] -
979  previous_alpha_over_beta *
980  this->preconditioner.apply(j, r[j]);
981  }
982  }
983  }
984 
985  // Function that we use if the PreconditionerType does not implement an
986  // apply() function, where we instead need to choose the
987  // apply_to_subrange function
988  template <typename U = void>
989  typename std::enable_if<!has_apply<PreconditionerType>, U>::type
990  operation_before_loop(const unsigned int iteration_index,
991  const unsigned int start_range,
992  const unsigned int end_range) const
993  {
994  Number * x = this->x.begin() + start_range;
995  Number * r = this->r.begin() + start_range;
996  Number * p = this->p.begin() + start_range;
997  Number * v = this->v.begin() + start_range;
998  const Number alpha = this->alpha;
999  const Number beta = this->beta;
1000  constexpr unsigned int grain_size = 128;
1001  std::array<Number, grain_size> prec_r;
1002  if (iteration_index == 1)
1003  {
1004  for (unsigned int j = start_range; j < end_range; j += grain_size)
1005  {
1006  const unsigned int length = std::min(grain_size, end_range - j);
1007  this->preconditioner.apply_to_subrange(j,
1008  j + length,
1009  r,
1010  prec_r.data());
1012  for (unsigned int i = 0; i < length; ++i)
1013  {
1014  p[i] = prec_r[i];
1015  v[i] = Number();
1016  }
1017  p += length;
1018  r += length;
1019  v += length;
1020  }
1021  }
1022  else if (iteration_index % 2 == 0)
1023  {
1024  for (unsigned int j = start_range; j < end_range; j += grain_size)
1025  {
1026  const unsigned int length = std::min(grain_size, end_range - j);
1028  for (unsigned int i = 0; i < length; ++i)
1029  r[i] -= this->alpha * v[i];
1030  this->preconditioner.apply_to_subrange(j,
1031  j + length,
1032  r,
1033  prec_r.data());
1035  for (unsigned int i = 0; i < length; ++i)
1036  {
1037  p[i] = this->beta * p[i] + prec_r[i];
1038  v[i] = Number();
1039  }
1040  p += length;
1041  r += length;
1042  v += length;
1043  }
1044  }
1045  else
1046  {
1047  const Number alpha_plus_previous_alpha_over_beta =
1048  this->alpha + this->previous_alpha / this->previous_beta;
1049  const Number previous_alpha_over_beta =
1050  this->previous_alpha / this->previous_beta;
1051  for (unsigned int j = start_range; j < end_range; j += grain_size)
1052  {
1053  const unsigned int length = std::min(grain_size, end_range - j);
1054  this->preconditioner.apply_to_subrange(j,
1055  j + length,
1056  r,
1057  prec_r.data());
1059  for (unsigned int i = 0; i < length; ++i)
1060  {
1061  x[i] += alpha_plus_previous_alpha_over_beta * p[i] -
1062  previous_alpha_over_beta * prec_r[i];
1063  r[i] -= this->alpha * v[i];
1064  }
1065  this->preconditioner.apply_to_subrange(j,
1066  j + length,
1067  r,
1068  prec_r.data());
1070  for (unsigned int i = 0; i < length; ++i)
1071  {
1072  p[i] = this->beta * p[i] + prec_r[i];
1073  v[i] = Number();
1074  }
1075  p += length;
1076  r += length;
1077  v += length;
1078  x += length;
1079  }
1080  }
1081  }
1082 
1083  // Function that we use if the PreconditionerType does not implement an
1084  // apply() function and where we instead need to use the
1085  // apply_to_subrange function
1086  template <typename U = void>
1087  typename std::enable_if<!has_apply<PreconditionerType>, U>::type
1088  operation_after_loop(
1089  const unsigned int start_range,
1090  const unsigned int end_range,
1091  std::array<VectorizedArray<Number>, 7> &vectorized_sums) const
1092  {
1093  const Number * r = this->r.begin();
1094  const Number * p = this->p.begin();
1095  const Number * v = this->v.begin();
1096  std::array<VectorizedArray<Number>, 7> my_sums = {};
1097  constexpr unsigned int grain_size = 128;
1098  Assert(grain_size % VectorizedArray<Number>::size() == 0,
1099  ExcNotImplemented());
1100  const unsigned int end_regular =
1101  start_range + (end_range - start_range) / grain_size * grain_size;
1102  std::array<Number, grain_size> prec_r;
1103  std::array<Number, grain_size> prec_v;
1104  for (unsigned int j = start_range; j < end_regular; j += grain_size)
1105  {
1106  this->preconditioner.apply_to_subrange(j,
1107  j + grain_size,
1108  r + j,
1109  prec_r.data());
1110  this->preconditioner.apply_to_subrange(j,
1111  j + grain_size,
1112  v + j,
1113  prec_v.data());
1114  VectorizedArray<Number> pj, vj, rj, prec_vj, prec_rj;
1115  for (unsigned int i = 0; i < grain_size;
1117  {
1118  pj.load(p + j + i);
1119  vj.load(v + j + i);
1120  rj.load(r + j + i);
1121  prec_rj.load(prec_r.data() + i);
1122  prec_vj.load(prec_v.data() + i);
1123 
1124  my_sums[0] += pj * vj;
1125  my_sums[1] += vj * vj;
1126  my_sums[2] += rj * vj;
1127  my_sums[3] += rj * rj;
1128  my_sums[4] += rj * prec_vj;
1129  my_sums[5] += vj * prec_vj;
1130  my_sums[6] += rj * prec_rj;
1131  }
1132  }
1133  const unsigned int length = end_range - end_regular;
1134  AssertIndexRange(length, grain_size);
1135  this->preconditioner.apply_to_subrange(end_regular,
1136  end_regular + length,
1137  r + end_regular,
1138  prec_r.data());
1139  this->preconditioner.apply_to_subrange(end_regular,
1140  end_regular + length,
1141  v + end_regular,
1142  prec_v.data());
1143  for (unsigned int j = end_regular; j < end_range; ++j)
1144  {
1145  my_sums[0][0] += p[j] * v[j];
1146  my_sums[1][0] += v[j] * v[j];
1147  my_sums[2][0] += r[j] * v[j];
1148  my_sums[3][0] += r[j] * r[j];
1149  my_sums[4][0] += r[j] * prec_v[j - end_regular];
1150  my_sums[5][0] += v[j] * prec_v[j - end_regular];
1151  my_sums[6][0] += r[j] * prec_r[j - end_regular];
1152  }
1153  for (unsigned int i = 0; i < vectorized_sums.size(); ++i)
1154  vectorized_sums[i] += my_sums[i];
1155  }
1156 
1157  // Function that we use if the PreconditionerType does not implement an
1158  // apply() function, where we instead need to choose the
1159  // apply_to_subrange function
1160  template <typename U = void>
1161  typename std::enable_if<!has_apply<PreconditionerType>, U>::type
1162  finalize_after_convergence(const unsigned int iteration_index)
1163  {
1164  if (iteration_index % 2 == 1 || iteration_index == 2)
1165  this->x.add(this->alpha, this->p);
1166  else
1167  {
1168  const unsigned int end_range = this->x.locally_owned_size();
1169 
1170  Number * x = this->x.begin();
1171  Number * r = this->r.begin();
1172  Number * p = this->p.begin();
1173  const Number alpha_plus_previous_alpha_over_beta =
1174  this->alpha + this->previous_alpha / this->previous_beta;
1175  const Number previous_alpha_over_beta =
1176  this->previous_alpha / this->previous_beta;
1177 
1178  constexpr unsigned int grain_size = 128;
1179  std::array<Number, grain_size> prec_r;
1180  for (unsigned int j = 0; j < end_range; j += grain_size)
1181  {
1182  const unsigned int length = std::min(grain_size, end_range - j);
1183  this->preconditioner.apply_to_subrange(j,
1184  j + length,
1185  r,
1186  prec_r.data());
1188  for (unsigned int i = 0; i < length; ++i)
1189  x[i] += alpha_plus_previous_alpha_over_beta * p[i] -
1190  previous_alpha_over_beta * prec_r[i];
1191 
1192  x += length;
1193  r += length;
1194  p += length;
1195  }
1196  }
1197  }
1198  };
1199  } // namespace SolverCG
1200 } // namespace internal
1201 
1202 
1203 
1204 template <typename VectorType>
1205 template <typename MatrixType, typename PreconditionerType>
1206 void
1207 SolverCG<VectorType>::solve(const MatrixType & A,
1208  VectorType & x,
1209  const VectorType & b,
1210  const PreconditionerType &preconditioner)
1211 {
1212  using number = typename VectorType::value_type;
1213 
1215 
1216  LogStream::Prefix prefix("cg");
1217 
1218  // Should we build the matrix for eigenvalue computations?
1219  const bool do_eigenvalues =
1220  !condition_number_signal.empty() || !all_condition_numbers_signal.empty() ||
1221  !eigenvalues_signal.empty() || !all_eigenvalues_signal.empty();
1222 
1223  // vectors used for eigenvalue computations
1224  std::vector<typename VectorType::value_type> diagonal;
1225  std::vector<typename VectorType::value_type> offdiagonal;
1226 
1227  typename VectorType::value_type eigen_beta_alpha = 0;
1228 
1229  int it = 0;
1230 
1231  internal::SolverCG::
1232  IterationWorker<VectorType, MatrixType, PreconditionerType>
1233  worker(
1234  A, preconditioner, determine_beta_by_flexible_formula, this->memory, x);
1235 
1236  worker.startup(b);
1237 
1238  solver_state = this->iteration_status(0, worker.residual_norm, x);
1239  if (solver_state != SolverControl::iterate)
1240  return;
1241 
1242  while (solver_state == SolverControl::iterate)
1243  {
1244  it++;
1245 
1246  worker.do_iteration(it);
1247 
1248  print_vectors(it, x, worker.r, worker.p);
1249 
1250  if (it > 1)
1251  {
1252  this->coefficients_signal(worker.previous_alpha, worker.beta);
1253  // set up the vectors containing the diagonal and the off diagonal
1254  // of the projected matrix.
1255  if (do_eigenvalues)
1256  {
1257  diagonal.push_back(number(1.) / worker.previous_alpha +
1258  eigen_beta_alpha);
1259  eigen_beta_alpha = worker.beta / worker.previous_alpha;
1260  offdiagonal.push_back(std::sqrt(worker.beta) /
1261  worker.previous_alpha);
1262  }
1263  compute_eigs_and_cond(diagonal,
1264  offdiagonal,
1265  all_eigenvalues_signal,
1266  all_condition_numbers_signal);
1267  }
1268 
1269  solver_state = this->iteration_status(it, worker.residual_norm, x);
1270  }
1271 
1272  worker.finalize_after_convergence(it);
1273 
1274  compute_eigs_and_cond(diagonal,
1275  offdiagonal,
1276  eigenvalues_signal,
1277  condition_number_signal);
1278 
1279  AssertThrow(solver_state == SolverControl::success,
1280  SolverControl::NoConvergence(it, worker.residual_norm));
1281 }
1282 
1283 
1284 
1285 template <typename VectorType>
1286 boost::signals2::connection
1288  const std::function<void(typename VectorType::value_type,
1289  typename VectorType::value_type)> &slot)
1290 {
1291  return coefficients_signal.connect(slot);
1292 }
1293 
1294 
1295 
1296 template <typename VectorType>
1297 boost::signals2::connection
1299  const std::function<void(double)> &slot,
1300  const bool every_iteration)
1301 {
1302  if (every_iteration)
1303  {
1304  return all_condition_numbers_signal.connect(slot);
1305  }
1306  else
1307  {
1308  return condition_number_signal.connect(slot);
1309  }
1310 }
1311 
1312 
1313 
1314 template <typename VectorType>
1315 boost::signals2::connection
1317  const std::function<void(const std::vector<double> &)> &slot,
1318  const bool every_iteration)
1319 {
1320  if (every_iteration)
1321  {
1322  return all_eigenvalues_signal.connect(slot);
1323  }
1324  else
1325  {
1326  return eigenvalues_signal.connect(slot);
1327  }
1328 }
1329 
1330 
1331 
1332 template <typename VectorType>
1335  const AdditionalData &)
1336  : SolverCG<VectorType>(cn, mem)
1337 {
1339 }
1340 
1341 
1342 
1343 template <typename VectorType>
1345  const AdditionalData &)
1346  : SolverCG<VectorType>(cn)
1347 {
1349 }
1350 
1351 
1352 
1353 #endif // DOXYGEN
1354 
1356 
1357 #endif
SolverCG(SolverControl &cn, VectorMemory< VectorType > &mem, const AdditionalData &data=AdditionalData())
boost::signals2::signal< void(const std::vector< double > &)> eigenvalues_signal
Definition: solver_cg.h:306
bool determine_beta_by_flexible_formula
Definition: solver_cg.h:323
virtual void print_vectors(const unsigned int step, const VectorType &x, const VectorType &r, const VectorType &d) const
boost::signals2::signal< void(typename VectorType::value_type, typename VectorType::value_type)> coefficients_signal
Definition: solver_cg.h:288
boost::signals2::connection connect_condition_number_slot(const std::function< void(double)> &slot, const bool every_iteration=false)
boost::signals2::connection connect_eigenvalues_slot(const std::function< void(const std::vector< double > &)> &slot, const bool every_iteration=false)
boost::signals2::signal< void(double)> condition_number_signal
Definition: solver_cg.h:294
boost::signals2::connection connect_coefficients_slot(const std::function< void(typename VectorType::value_type, typename VectorType::value_type)> &slot)
virtual ~SolverCG() override=default
void solve(const MatrixType &A, VectorType &x, const VectorType &b, const PreconditionerType &preconditioner)
boost::signals2::signal< void(double)> all_condition_numbers_signal
Definition: solver_cg.h:300
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(Number *ptr) const
void load(const Number *ptr)
#define DEAL_II_OPENMP_SIMD_PRAGMA
Definition: config.h:142
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
static ::ExceptionBase & ExcDivideByZero()
#define Assert(cond, exc)
Definition: exceptions.h:1473
static ::ExceptionBase & ExcNotImplemented()
#define AssertIndexRange(index, range)
Definition: exceptions.h:1732
#define AssertThrow(cond, exc)
Definition: exceptions.h:1583
static const char U
static const char A
@ diagonal
Matrix is diagonal.
static const char T
types::global_dof_index size_type
Definition: cuda_kernels.h:45
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:76
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)