Reference documentation for deal.II version GIT 9c182271f7 2023-03-28 14:30: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\}}\)
precondition.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1999 - 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_precondition_h
17 #define dealii_precondition_h
18 
19 #include <deal.II/base/config.h>
20 
21 #include <deal.II/base/cuda_size.h>
23 #include <deal.II/base/mutex.h>
24 #include <deal.II/base/parallel.h>
27 
32 #include <deal.II/lac/solver_cg.h>
34 
35 #include <limits>
36 
38 
39 // forward declarations
40 #ifndef DOXYGEN
41 template <typename number>
42 class Vector;
43 template <typename number>
44 class SparseMatrix;
45 namespace LinearAlgebra
46 {
47  namespace distributed
48  {
49  template <typename, typename>
50  class Vector;
51  template <typename>
52  class BlockVector;
53  } // namespace distributed
54 } // namespace LinearAlgebra
55 #endif
56 
57 
84 {
85 public:
90 
96  {
100  AdditionalData() = default;
101  };
102 
107 
112  template <typename MatrixType>
113  void
114  initialize(const MatrixType & matrix,
115  const AdditionalData &additional_data = AdditionalData());
116 
120  template <class VectorType>
121  void
122  vmult(VectorType &, const VectorType &) const;
123 
128  template <class VectorType>
129  void
130  Tvmult(VectorType &, const VectorType &) const;
131 
135  template <class VectorType>
136  void
137  vmult_add(VectorType &, const VectorType &) const;
138 
143  template <class VectorType>
144  void
145  Tvmult_add(VectorType &, const VectorType &) const;
146 
151  void
152  clear();
153 
161  size_type
162  m() const;
163 
171  size_type
172  n() const;
173 
174 private:
179 
184 };
185 
186 
187 
199 {
200 public:
205 
210  {
211  public:
216  AdditionalData(const double relaxation = 1.);
217 
221  double relaxation;
222  };
223 
229 
233  void
234  initialize(const AdditionalData &parameters);
235 
241  template <typename MatrixType>
242  void
243  initialize(const MatrixType &matrix, const AdditionalData &parameters);
244 
248  template <class VectorType>
249  void
250  vmult(VectorType &, const VectorType &) const;
251 
256  template <class VectorType>
257  void
258  Tvmult(VectorType &, const VectorType &) const;
262  template <class VectorType>
263  void
264  vmult_add(VectorType &, const VectorType &) const;
265 
270  template <class VectorType>
271  void
272  Tvmult_add(VectorType &, const VectorType &) const;
273 
278  void
280  {}
281 
289  size_type
290  m() const;
291 
299  size_type
300  n() const;
301 
302 private:
306  double relaxation;
307 
312 
317 };
318 
319 
320 
360 template <typename MatrixType = SparseMatrix<double>,
361  class VectorType = Vector<double>>
363 {
364 public:
368  using function_ptr = void (MatrixType::*)(VectorType &,
369  const VectorType &) const;
370 
376  PreconditionUseMatrix(const MatrixType &M, const function_ptr method);
377 
382  void
383  vmult(VectorType &dst, const VectorType &src) const;
384 
385 private:
389  const MatrixType &matrix;
390 
395 };
396 
397 
398 
404 template <typename MatrixType = SparseMatrix<double>,
405  typename PreconditionerType = IdentityMatrix>
407 {
408 public:
413 
418  {
419  public:
423  AdditionalData(const double relaxation = 1.,
424  const unsigned int n_iterations = 1);
425 
429  double relaxation;
430 
434  unsigned int n_iterations;
435 
436 
437  /*
438  * Preconditioner.
439  */
440  std::shared_ptr<PreconditionerType> preconditioner;
441  };
442 
448  void
449  initialize(const MatrixType & A,
450  const AdditionalData &parameters = AdditionalData());
451 
455  void
456  clear();
457 
462  size_type
463  m() const;
464 
469  size_type
470  n() const;
471 
475  template <class VectorType>
476  void
477  vmult(VectorType &, const VectorType &) const;
478 
483  template <class VectorType>
484  void
485  Tvmult(VectorType &, const VectorType &) const;
486 
490  template <class VectorType>
491  void
492  step(VectorType &x, const VectorType &rhs) const;
493 
497  template <class VectorType>
498  void
499  Tstep(VectorType &x, const VectorType &rhs) const;
500 
501 protected:
506 
510  double relaxation;
511 
515  unsigned int n_iterations;
516 
517  /*
518  * Preconditioner.
519  */
520  std::shared_ptr<PreconditionerType> preconditioner;
521 };
522 
523 
524 
525 #ifndef DOXYGEN
526 
527 namespace internal
528 {
529  // a helper type-trait that leverage SFINAE to figure out if MatrixType has
530  // ... MatrixType::vmult(VectorType &, const VectorType&,
531  // std::function<...>, std::function<...>) const
532  template <typename MatrixType, typename VectorType>
533  using vmult_functions_t = decltype(std::declval<MatrixType const>().vmult(
534  std::declval<VectorType &>(),
535  std::declval<const VectorType &>(),
536  std::declval<
537  const std::function<void(const unsigned int, const unsigned int)> &>(),
538  std::declval<
539  const std::function<void(const unsigned int, const unsigned int)> &>()));
540 
541  template <typename MatrixType,
542  typename VectorType,
543  typename PreconditionerType>
544  constexpr bool has_vmult_with_std_functions =
545  is_supported_operation<vmult_functions_t, MatrixType, VectorType> &&
546  std::is_same<PreconditionerType,
548  (std::is_same<VectorType,
550  std::is_same<
551  VectorType,
552  LinearAlgebra::distributed::Vector<typename VectorType::value_type,
553  MemorySpace::Host>>::value);
554 
555 
556  template <typename MatrixType, typename VectorType>
557  constexpr bool has_vmult_with_std_functions_for_precondition =
558  is_supported_operation<vmult_functions_t, MatrixType, VectorType>;
559 
560  namespace PreconditionRelaxation
561  {
562  template <typename T, typename VectorType>
563  using Tvmult_t = decltype(
564  std::declval<T const>().Tvmult(std::declval<VectorType &>(),
565  std::declval<const VectorType &>()));
566 
567  template <typename T, typename VectorType>
568  constexpr bool has_Tvmult = is_supported_operation<Tvmult_t, T, VectorType>;
569 
570  template <typename T, typename VectorType>
571  using step_t = decltype(
572  std::declval<T const>().step(std::declval<VectorType &>(),
573  std::declval<const VectorType &>()));
574 
575  template <typename T, typename VectorType>
576  constexpr bool has_step = is_supported_operation<step_t, T, VectorType>;
577 
578  template <typename T, typename VectorType>
579  using step_omega_t =
580  decltype(std::declval<T const>().step(std::declval<VectorType &>(),
581  std::declval<const VectorType &>(),
582  std::declval<const double>()));
583 
584  template <typename T, typename VectorType>
585  constexpr bool has_step_omega =
586  is_supported_operation<step_omega_t, T, VectorType>;
587 
588  template <typename T, typename VectorType>
589  using Tstep_t = decltype(
590  std::declval<T const>().Tstep(std::declval<VectorType &>(),
591  std::declval<const VectorType &>()));
592 
593  template <typename T, typename VectorType>
594  constexpr bool has_Tstep = is_supported_operation<Tstep_t, T, VectorType>;
595 
596  template <typename T, typename VectorType>
597  using Tstep_omega_t =
598  decltype(std::declval<T const>().Tstep(std::declval<VectorType &>(),
599  std::declval<const VectorType &>(),
600  std::declval<const double>()));
601 
602  template <typename T, typename VectorType>
603  constexpr bool has_Tstep_omega =
604  is_supported_operation<Tstep_omega_t, T, VectorType>;
605 
606  template <typename T, typename VectorType>
607  using jacobi_step_t = decltype(
608  std::declval<T const>().Jacobi_step(std::declval<VectorType &>(),
609  std::declval<const VectorType &>(),
610  std::declval<const double>()));
611 
612  template <typename T, typename VectorType>
613  constexpr bool has_jacobi_step =
614  is_supported_operation<jacobi_step_t, T, VectorType>;
615 
616  template <typename T, typename VectorType>
617  using SOR_step_t = decltype(
618  std::declval<T const>().SOR_step(std::declval<VectorType &>(),
619  std::declval<const VectorType &>(),
620  std::declval<const double>()));
621 
622  template <typename T, typename VectorType>
623  constexpr bool has_SOR_step =
624  is_supported_operation<SOR_step_t, T, VectorType>;
625 
626  template <typename T, typename VectorType>
627  using SSOR_step_t = decltype(
628  std::declval<T const>().SSOR_step(std::declval<VectorType &>(),
629  std::declval<const VectorType &>(),
630  std::declval<const double>()));
631 
632  template <typename T, typename VectorType>
633  constexpr bool has_SSOR_step =
634  is_supported_operation<SSOR_step_t, T, VectorType>;
635 
636  template <typename MatrixType>
637  class PreconditionJacobiImpl
638  {
639  public:
640  PreconditionJacobiImpl(const MatrixType &A, const double relaxation)
641  : A(&A)
642  , relaxation(relaxation)
643  {}
644 
645  template <typename VectorType>
646  void
647  vmult(VectorType &dst, const VectorType &src) const
648  {
649  this->A->precondition_Jacobi(dst, src, this->relaxation);
650  }
651 
652  template <typename VectorType>
653  void
654  Tvmult(VectorType &dst, const VectorType &src) const
655  {
656  // call vmult, since preconditioner is symmetrical
657  this->vmult(dst, src);
658  }
659 
660  template <typename VectorType,
661  std::enable_if_t<has_jacobi_step<MatrixType, VectorType>,
662  MatrixType> * = nullptr>
663  void
664  step(VectorType &dst, const VectorType &src) const
665  {
666  this->A->Jacobi_step(dst, src, this->relaxation);
667  }
668 
669  template <typename VectorType,
670  std::enable_if_t<!has_jacobi_step<MatrixType, VectorType>,
671  MatrixType> * = nullptr>
672  void
673  step(VectorType &, const VectorType &) const
674  {
675  AssertThrow(false,
676  ExcMessage(
677  "Matrix A does not provide a Jacobi_step() function!"));
678  }
679 
680  template <typename VectorType>
681  void
682  Tstep(VectorType &dst, const VectorType &src) const
683  {
684  // call step, since preconditioner is symmetrical
685  this->step(dst, src);
686  }
687 
688  private:
690  const double relaxation;
691  };
692 
693  template <typename MatrixType>
694  class PreconditionSORImpl
695  {
696  public:
697  PreconditionSORImpl(const MatrixType &A, const double relaxation)
698  : A(&A)
699  , relaxation(relaxation)
700  {}
701 
702  template <typename VectorType>
703  void
704  vmult(VectorType &dst, const VectorType &src) const
705  {
706  this->A->precondition_SOR(dst, src, this->relaxation);
707  }
708 
709  template <typename VectorType>
710  void
711  Tvmult(VectorType &dst, const VectorType &src) const
712  {
713  this->A->precondition_TSOR(dst, src, this->relaxation);
714  }
715 
716  template <typename VectorType,
717  std::enable_if_t<has_SOR_step<MatrixType, VectorType>,
718  MatrixType> * = nullptr>
719  void
720  step(VectorType &dst, const VectorType &src) const
721  {
722  this->A->SOR_step(dst, src, this->relaxation);
723  }
724 
725  template <typename VectorType,
726  std::enable_if_t<!has_SOR_step<MatrixType, VectorType>,
727  MatrixType> * = nullptr>
728  void
729  step(VectorType &, const VectorType &) const
730  {
731  AssertThrow(false,
732  ExcMessage(
733  "Matrix A does not provide a SOR_step() function!"));
734  }
735 
736  template <typename VectorType,
737  std::enable_if_t<has_SOR_step<MatrixType, VectorType>,
738  MatrixType> * = nullptr>
739  void
740  Tstep(VectorType &dst, const VectorType &src) const
741  {
742  this->A->TSOR_step(dst, src, this->relaxation);
743  }
744 
745  template <typename VectorType,
746  std::enable_if_t<!has_SOR_step<MatrixType, VectorType>,
747  MatrixType> * = nullptr>
748  void
749  Tstep(VectorType &, const VectorType &) const
750  {
751  AssertThrow(false,
752  ExcMessage(
753  "Matrix A does not provide a TSOR_step() function!"));
754  }
755 
756  private:
758  const double relaxation;
759  };
760 
761  template <typename MatrixType>
762  class PreconditionSSORImpl
763  {
764  public:
765  using size_type = typename MatrixType::size_type;
766 
767  PreconditionSSORImpl(const MatrixType &A, const double relaxation)
768  : A(&A)
769  , relaxation(relaxation)
770  {
771  // in case we have a SparseMatrix class, we can extract information
772  // about the diagonal.
774  dynamic_cast<const SparseMatrix<typename MatrixType::value_type> *>(
775  &*this->A);
776 
777  // calculate the positions first after the diagonal.
778  if (mat != nullptr)
779  {
780  const size_type n = this->A->n();
781  pos_right_of_diagonal.resize(n, static_cast<std::size_t>(-1));
782  for (size_type row = 0; row < n; ++row)
783  {
784  // find the first element in this line which is on the right of
785  // the diagonal. we need to precondition with the elements on
786  // the left only. note: the first entry in each line denotes the
787  // diagonal element, which we need not check.
788  typename SparseMatrix<
789  typename MatrixType::value_type>::const_iterator it =
790  mat->begin(row) + 1;
791  for (; it < mat->end(row); ++it)
792  if (it->column() > row)
793  break;
794  pos_right_of_diagonal[row] = it - mat->begin();
795  }
796  }
797  }
798 
799  template <typename VectorType>
800  void
801  vmult(VectorType &dst, const VectorType &src) const
802  {
803  this->A->precondition_SSOR(dst,
804  src,
805  this->relaxation,
806  pos_right_of_diagonal);
807  }
808 
809  template <typename VectorType>
810  void
811  Tvmult(VectorType &dst, const VectorType &src) const
812  {
813  this->A->precondition_SSOR(dst,
814  src,
815  this->relaxation,
816  pos_right_of_diagonal);
817  }
818 
819  template <typename VectorType,
820  std::enable_if_t<has_SSOR_step<MatrixType, VectorType>,
821  MatrixType> * = nullptr>
822  void
823  step(VectorType &dst, const VectorType &src) const
824  {
825  this->A->SSOR_step(dst, src, this->relaxation);
826  }
827 
828  template <typename VectorType,
829  std::enable_if_t<!has_SSOR_step<MatrixType, VectorType>,
830  MatrixType> * = nullptr>
831  void
832  step(VectorType &, const VectorType &) const
833  {
834  AssertThrow(false,
835  ExcMessage(
836  "Matrix A does not provide a SSOR_step() function!"));
837  }
838 
839  template <typename VectorType>
840  void
841  Tstep(VectorType &dst, const VectorType &src) const
842  {
843  // call step, since preconditioner is symmetrical
844  this->step(dst, src);
845  }
846 
847  private:
849  const double relaxation;
850 
855  std::vector<std::size_t> pos_right_of_diagonal;
856  };
857 
858  template <typename MatrixType>
859  class PreconditionPSORImpl
860  {
861  public:
862  using size_type = typename MatrixType::size_type;
863 
864  PreconditionPSORImpl(const MatrixType & A,
865  const double relaxation,
866  const std::vector<size_type> &permutation,
867  const std::vector<size_type> &inverse_permutation)
868  : A(&A)
869  , relaxation(relaxation)
870  , permutation(permutation)
871  , inverse_permutation(inverse_permutation)
872  {}
873 
874  template <typename VectorType>
875  void
876  vmult(VectorType &dst, const VectorType &src) const
877  {
878  dst = src;
879  this->A->PSOR(dst, permutation, inverse_permutation, this->relaxation);
880  }
881 
882  template <typename VectorType>
883  void
884  Tvmult(VectorType &dst, const VectorType &src) const
885  {
886  dst = src;
887  this->A->TPSOR(dst, permutation, inverse_permutation, this->relaxation);
888  }
889 
890  private:
892  const double relaxation;
893 
894  const std::vector<size_type> &permutation;
895  const std::vector<size_type> &inverse_permutation;
896  };
897 
898  template <typename MatrixType,
899  typename PreconditionerType,
900  typename VectorType,
901  std::enable_if_t<has_step_omega<PreconditionerType, VectorType>,
902  PreconditionerType> * = nullptr>
903  void
904  step(const MatrixType &,
905  const PreconditionerType &preconditioner,
906  VectorType & dst,
907  const VectorType & src,
908  const double relaxation,
909  VectorType &,
910  VectorType &)
911  {
912  preconditioner.step(dst, src, relaxation);
913  }
914 
915  template <
916  typename MatrixType,
917  typename PreconditionerType,
918  typename VectorType,
919  std::enable_if_t<!has_step_omega<PreconditionerType, VectorType> &&
920  has_step<PreconditionerType, VectorType>,
921  PreconditionerType> * = nullptr>
922  void
923  step(const MatrixType &,
924  const PreconditionerType &preconditioner,
925  VectorType & dst,
926  const VectorType & src,
927  const double relaxation,
928  VectorType &,
929  VectorType &)
930  {
931  Assert(relaxation == 1.0, ExcInternalError());
932 
933  (void)relaxation;
934 
935  preconditioner.step(dst, src);
936  }
937 
938  template <
939  typename MatrixType,
940  typename PreconditionerType,
941  typename VectorType,
942  std::enable_if_t<!has_step_omega<PreconditionerType, VectorType> &&
943  !has_step<PreconditionerType, VectorType>,
944  PreconditionerType> * = nullptr>
945  void
946  step(const MatrixType & A,
947  const PreconditionerType &preconditioner,
948  VectorType & dst,
949  const VectorType & src,
950  const double relaxation,
951  VectorType & residual,
952  VectorType & tmp)
953  {
954  residual.reinit(dst, true);
955  tmp.reinit(dst, true);
956 
957  A.vmult(residual, dst);
958  residual.sadd(-1.0, 1.0, src);
959 
960  preconditioner.vmult(tmp, residual);
961  dst.add(relaxation, tmp);
962  }
963 
964  template <typename MatrixType,
965  typename PreconditionerType,
966  typename VectorType,
967  std::enable_if_t<has_Tstep_omega<PreconditionerType, VectorType>,
968  PreconditionerType> * = nullptr>
969  void
970  Tstep(const MatrixType &,
971  const PreconditionerType &preconditioner,
972  VectorType & dst,
973  const VectorType & src,
974  const double relaxation,
975  VectorType &,
976  VectorType &)
977  {
978  preconditioner.Tstep(dst, src, relaxation);
979  }
980 
981  template <
982  typename MatrixType,
983  typename PreconditionerType,
984  typename VectorType,
985  std::enable_if_t<!has_Tstep_omega<PreconditionerType, VectorType> &&
986  has_Tstep<PreconditionerType, VectorType>,
987  PreconditionerType> * = nullptr>
988  void
989  Tstep(const MatrixType &,
990  const PreconditionerType &preconditioner,
991  VectorType & dst,
992  const VectorType & src,
993  const double relaxation,
994  VectorType &,
995  VectorType &)
996  {
997  Assert(relaxation == 1.0, ExcInternalError());
998 
999  (void)relaxation;
1000 
1001  preconditioner.Tstep(dst, src);
1002  }
1003 
1004  template <typename MatrixType,
1005  typename VectorType,
1006  std::enable_if_t<has_Tvmult<MatrixType, VectorType>, MatrixType>
1007  * = nullptr>
1008  void
1009  Tvmult(const MatrixType &A, VectorType &dst, const VectorType &src)
1010  {
1011  A.Tvmult(dst, src);
1012  }
1013 
1014  template <typename MatrixType,
1015  typename VectorType,
1016  std::enable_if_t<!has_Tvmult<MatrixType, VectorType>, MatrixType>
1017  * = nullptr>
1018  void
1019  Tvmult(const MatrixType &, VectorType &, const VectorType &)
1020  {
1021  AssertThrow(false,
1022  ExcMessage("Matrix A does not provide a Tvmult() function!"));
1023  }
1024 
1025  template <
1026  typename MatrixType,
1027  typename PreconditionerType,
1028  typename VectorType,
1029  std::enable_if_t<!has_Tstep_omega<PreconditionerType, VectorType> &&
1030  !has_Tstep<PreconditionerType, VectorType>,
1031  PreconditionerType> * = nullptr>
1032  void
1033  Tstep(const MatrixType & A,
1034  const PreconditionerType &preconditioner,
1035  VectorType & dst,
1036  const VectorType & src,
1037  const double relaxation,
1038  VectorType & residual,
1039  VectorType & tmp)
1040  {
1041  residual.reinit(dst, true);
1042  tmp.reinit(dst, true);
1043 
1044  Tvmult(A, residual, dst);
1045  residual.sadd(-1.0, 1.0, src);
1046 
1047  Tvmult(preconditioner, tmp, residual);
1048  dst.add(relaxation, tmp);
1049  }
1050 
1051  // 0) general implementation
1052  template <typename MatrixType,
1053  typename PreconditionerType,
1054  typename VectorType,
1055  std::enable_if_t<!has_vmult_with_std_functions_for_precondition<
1056  PreconditionerType,
1057  VectorType>,
1058  int> * = nullptr>
1059  void
1060  step_operations(const MatrixType & A,
1061  const PreconditionerType &preconditioner,
1062  VectorType & dst,
1063  const VectorType & src,
1064  const double relaxation,
1065  VectorType & tmp1,
1066  VectorType & tmp2,
1067  const unsigned int i,
1068  const bool transposed)
1069  {
1070  if (i == 0)
1071  {
1072  if (transposed)
1073  Tvmult(preconditioner, dst, src);
1074  else
1075  preconditioner.vmult(dst, src);
1076 
1077  if (relaxation != 1.0)
1078  dst *= relaxation;
1079  }
1080  else
1081  {
1082  if (transposed)
1083  Tstep(A, preconditioner, dst, src, relaxation, tmp1, tmp2);
1084  else
1085  step(A, preconditioner, dst, src, relaxation, tmp1, tmp2);
1086  }
1087  }
1088 
1089  // 1) specialized implementation with a preconditioner that accepts
1090  // ranges
1091  template <
1092  typename MatrixType,
1093  typename PreconditionerType,
1094  typename VectorType,
1095  std::enable_if_t<
1096  has_vmult_with_std_functions_for_precondition<PreconditionerType,
1097  VectorType> &&
1098  !has_vmult_with_std_functions_for_precondition<MatrixType,
1099  VectorType>,
1100  int> * = nullptr>
1101  void
1102  step_operations(const MatrixType & A,
1103  const PreconditionerType &preconditioner,
1104  VectorType & dst,
1105  const VectorType & src,
1106  const double relaxation,
1107  VectorType & tmp,
1108  VectorType &,
1109  const unsigned int i,
1110  const bool transposed)
1111  {
1112  (void)transposed;
1113  using Number = typename VectorType::value_type;
1114 
1115  if (i == 0)
1116  {
1117  Number * dst_ptr = dst.begin();
1118  const Number *src_ptr = src.begin();
1119 
1120  preconditioner.vmult(
1121  dst,
1122  src,
1123  [&](const unsigned int start_range, const unsigned int end_range) {
1124  // zero 'dst' before running the vmult operation
1125  if (end_range > start_range)
1126  std::memset(dst.begin() + start_range,
1127  0,
1128  sizeof(Number) * (end_range - start_range));
1129  },
1130  [&](const unsigned int start_range, const unsigned int end_range) {
1131  if (relaxation == 1.0)
1132  return; // nothing to do
1133 
1134  const auto src_ptr = src.begin();
1135  const auto dst_ptr = dst.begin();
1136 
1138  for (std::size_t i = start_range; i < end_range; ++i)
1139  dst_ptr[i] *= relaxation;
1140  });
1141  }
1142  else
1143  {
1144  tmp.reinit(src, true);
1145 
1146  Assert(transposed == false, ExcNotImplemented());
1147 
1148  A.vmult(tmp, dst);
1149 
1150  preconditioner.vmult(
1151  dst,
1152  tmp,
1153  [&](const unsigned int start_range, const unsigned int end_range) {
1154  const auto src_ptr = src.begin();
1155  const auto tmp_ptr = tmp.begin();
1156 
1157  if (relaxation == 1.0)
1158  {
1160  for (std::size_t i = start_range; i < end_range; ++i)
1161  tmp_ptr[i] = src_ptr[i] - tmp_ptr[i];
1162  }
1163  else
1164  {
1165  // note: we scale the residual here to be able to add into
1166  // the dst vector, which contains the solution from the last
1167  // iteration
1169  for (std::size_t i = start_range; i < end_range; ++i)
1170  tmp_ptr[i] = relaxation * (src_ptr[i] - tmp_ptr[i]);
1171  }
1172  },
1173  [&](const unsigned int, const unsigned int) {
1174  // nothing to do, since scaling by the relaxation factor
1175  // has been done in the pre operation
1176  });
1177  }
1178  }
1179 
1180  // 2) specialized implementation with a preconditioner and a matrix that
1181  // accepts ranges
1182  template <
1183  typename MatrixType,
1184  typename PreconditionerType,
1185  typename VectorType,
1186  std::enable_if_t<
1187  has_vmult_with_std_functions_for_precondition<PreconditionerType,
1188  VectorType> &&
1189  has_vmult_with_std_functions_for_precondition<MatrixType, VectorType>,
1190  int> * = nullptr>
1191  void
1192  step_operations(const MatrixType & A,
1193  const PreconditionerType &preconditioner,
1194  VectorType & dst,
1195  const VectorType & src,
1196  const double relaxation,
1197  VectorType & tmp,
1198  VectorType &,
1199  const unsigned int i,
1200  const bool transposed)
1201  {
1202  (void)transposed;
1203  using Number = typename VectorType::value_type;
1204 
1205  if (i == 0)
1206  {
1207  Number * dst_ptr = dst.begin();
1208  const Number *src_ptr = src.begin();
1209 
1210  preconditioner.vmult(
1211  dst,
1212  src,
1213  [&](const unsigned int start_range, const unsigned int end_range) {
1214  // zero 'dst' before running the vmult operation
1215  if (end_range > start_range)
1216  std::memset(dst.begin() + start_range,
1217  0,
1218  sizeof(Number) * (end_range - start_range));
1219  },
1220  [&](const unsigned int start_range, const unsigned int end_range) {
1221  if (relaxation == 1.0)
1222  return; // nothing to do
1223 
1224  const auto src_ptr = src.begin();
1225  const auto dst_ptr = dst.begin();
1226 
1228  for (std::size_t i = start_range; i < end_range; ++i)
1229  dst_ptr[i] *= relaxation;
1230  });
1231  }
1232  else
1233  {
1234  tmp.reinit(src, true);
1235 
1236  Assert(transposed == false, ExcNotImplemented());
1237 
1238  A.vmult(
1239  tmp,
1240  dst,
1241  [&](const unsigned int start_range, const unsigned int end_range) {
1242  // zero 'tmp' before running the vmult
1243  // operation
1244  if (end_range > start_range)
1245  std::memset(tmp.begin() + start_range,
1246  0,
1247  sizeof(Number) * (end_range - start_range));
1248  },
1249  [&](const unsigned int start_range, const unsigned int end_range) {
1250  const auto src_ptr = src.begin();
1251  const auto tmp_ptr = tmp.begin();
1252 
1253  if (relaxation == 1.0)
1254  {
1256  for (std::size_t i = start_range; i < end_range; ++i)
1257  tmp_ptr[i] = src_ptr[i] - tmp_ptr[i];
1258  }
1259  else
1260  {
1261  // note: we scale the residual here to be able to add into
1262  // the dst vector, which contains the solution from the last
1263  // iteration
1265  for (std::size_t i = start_range; i < end_range; ++i)
1266  tmp_ptr[i] = relaxation * (src_ptr[i] - tmp_ptr[i]);
1267  }
1268  });
1269 
1270  preconditioner.vmult(dst, tmp, [](const auto, const auto) {
1271  // note: `dst` vector does not have to be zeroed out
1272  // since we add the result into it
1273  });
1274  }
1275  }
1276 
1277  // 3) specialized implementation for inverse-diagonal preconditioner
1278  template <typename MatrixType,
1279  typename VectorType,
1280  std::enable_if_t<!IsBlockVector<VectorType>::value &&
1281  !has_vmult_with_std_functions<
1282  MatrixType,
1283  VectorType,
1285  VectorType> * = nullptr>
1286  void
1287  step_operations(const MatrixType & A,
1288  const ::DiagonalMatrix<VectorType> &preconditioner,
1289  VectorType & dst,
1290  const VectorType & src,
1291  const double relaxation,
1292  VectorType & tmp,
1293  VectorType &,
1294  const unsigned int i,
1295  const bool transposed)
1296  {
1297  using Number = typename VectorType::value_type;
1298 
1299  if (i == 0)
1300  {
1301  Number * dst_ptr = dst.begin();
1302  const Number *src_ptr = src.begin();
1303  const Number *diag_ptr = preconditioner.get_vector().begin();
1304 
1305  if (relaxation == 1.0)
1306  {
1308  for (unsigned int i = 0; i < dst.locally_owned_size(); ++i)
1309  dst_ptr[i] = src_ptr[i] * diag_ptr[i];
1310  }
1311  else
1312  {
1314  for (unsigned int i = 0; i < dst.locally_owned_size(); ++i)
1315  dst_ptr[i] = relaxation * src_ptr[i] * diag_ptr[i];
1316  }
1317  }
1318  else
1319  {
1320  tmp.reinit(src, true);
1321 
1322  Number * dst_ptr = dst.begin();
1323  const Number *src_ptr = src.begin();
1324  const Number *tmp_ptr = tmp.begin();
1325  const Number *diag_ptr = preconditioner.get_vector().begin();
1326 
1327  if (transposed)
1328  Tvmult(A, tmp, dst);
1329  else
1330  A.vmult(tmp, dst);
1331 
1332  if (relaxation == 1.0)
1333  {
1335  for (unsigned int i = 0; i < dst.locally_owned_size(); ++i)
1336  dst_ptr[i] += (src_ptr[i] - tmp_ptr[i]) * diag_ptr[i];
1337  }
1338  else
1339  {
1341  for (unsigned int i = 0; i < dst.locally_owned_size(); ++i)
1342  dst_ptr[i] +=
1343  relaxation * (src_ptr[i] - tmp_ptr[i]) * diag_ptr[i];
1344  }
1345  }
1346  }
1347 
1348  // 4) specialized implementation for inverse-diagonal preconditioner and
1349  // matrix that accepts ranges
1350  template <typename MatrixType,
1351  typename VectorType,
1352  std::enable_if_t<!IsBlockVector<VectorType>::value &&
1353  has_vmult_with_std_functions<
1354  MatrixType,
1355  VectorType,
1357  VectorType> * = nullptr>
1358  void
1359  step_operations(const MatrixType & A,
1360  const ::DiagonalMatrix<VectorType> &preconditioner,
1361  VectorType & dst,
1362  const VectorType & src,
1363  const double relaxation,
1364  VectorType & tmp,
1365  VectorType &,
1366  const unsigned int i,
1367  const bool transposed)
1368  {
1369  (void)transposed;
1370  using Number = typename VectorType::value_type;
1371 
1372  if (i == 0)
1373  {
1374  Number * dst_ptr = dst.begin();
1375  const Number *src_ptr = src.begin();
1376  const Number *diag_ptr = preconditioner.get_vector().begin();
1377 
1378  if (relaxation == 1.0)
1379  {
1381  for (unsigned int i = 0; i < dst.locally_owned_size(); ++i)
1382  dst_ptr[i] = src_ptr[i] * diag_ptr[i];
1383  }
1384  else
1385  {
1387  for (unsigned int i = 0; i < dst.locally_owned_size(); ++i)
1388  dst_ptr[i] = relaxation * src_ptr[i] * diag_ptr[i];
1389  }
1390  }
1391  else
1392  {
1393  tmp.reinit(src, true);
1394 
1395  Assert(transposed == false, ExcNotImplemented());
1396 
1397  A.vmult(
1398  tmp,
1399  dst,
1400  [&](const unsigned int start_range, const unsigned int end_range) {
1401  // zero 'tmp' before running the vmult operation
1402  if (end_range > start_range)
1403  std::memset(tmp.begin() + start_range,
1404  0,
1405  sizeof(Number) * (end_range - start_range));
1406  },
1407  [&](const unsigned int begin, const unsigned int end) {
1408  const Number *dst_ptr = dst.begin();
1409  const Number *src_ptr = src.begin();
1410  Number * tmp_ptr = tmp.begin();
1411  const Number *diag_ptr = preconditioner.get_vector().begin();
1412 
1413  // for efficiency reason, write back to temp_vector that is
1414  // already read (avoid read-for-ownership)
1415  if (relaxation == 1.0)
1416  {
1418  for (std::size_t i = begin; i < end; ++i)
1419  tmp_ptr[i] =
1420  dst_ptr[i] + (src_ptr[i] - tmp_ptr[i]) * diag_ptr[i];
1421  }
1422  else
1423  {
1425  for (std::size_t i = begin; i < end; ++i)
1426  tmp_ptr[i] = dst_ptr[i] + relaxation *
1427  (src_ptr[i] - tmp_ptr[i]) *
1428  diag_ptr[i];
1429  }
1430  });
1431 
1432  tmp.swap(dst);
1433  }
1434  }
1435 
1436  } // namespace PreconditionRelaxation
1437 } // namespace internal
1438 
1439 #endif
1440 
1441 
1442 
1469 template <typename MatrixType = SparseMatrix<double>>
1471  : public PreconditionRelaxation<
1472  MatrixType,
1473  internal::PreconditionRelaxation::PreconditionJacobiImpl<MatrixType>>
1474 {
1476  internal::PreconditionRelaxation::PreconditionJacobiImpl<MatrixType>;
1478 
1479 public:
1484 
1488  void
1489  initialize(const MatrixType & A,
1490  const AdditionalData &parameters = AdditionalData());
1491 };
1492 
1493 
1539 template <typename MatrixType = SparseMatrix<double>>
1541  : public PreconditionRelaxation<
1542  MatrixType,
1543  internal::PreconditionRelaxation::PreconditionSORImpl<MatrixType>>
1544 {
1546  internal::PreconditionRelaxation::PreconditionSORImpl<MatrixType>;
1548 
1549 public:
1554 
1558  void
1559  initialize(const MatrixType & A,
1560  const AdditionalData &parameters = AdditionalData());
1561 };
1562 
1563 
1564 
1591 template <typename MatrixType = SparseMatrix<double>>
1593  : public PreconditionRelaxation<
1594  MatrixType,
1595  internal::PreconditionRelaxation::PreconditionSSORImpl<MatrixType>>
1596 {
1598  internal::PreconditionRelaxation::PreconditionSSORImpl<MatrixType>;
1600 
1601 public:
1606 
1612  void
1613  initialize(const MatrixType & A,
1614  const AdditionalData &parameters = AdditionalData());
1615 };
1616 
1617 
1647 template <typename MatrixType = SparseMatrix<double>>
1649  : public PreconditionRelaxation<
1650  MatrixType,
1651  internal::PreconditionRelaxation::PreconditionPSORImpl<MatrixType>>
1652 {
1654  internal::PreconditionRelaxation::PreconditionPSORImpl<MatrixType>;
1656 
1657 public:
1662 
1667  {
1668  public:
1679  AdditionalData(const std::vector<size_type> &permutation,
1680  const std::vector<size_type> &inverse_permutation,
1681  const typename BaseClass::AdditionalData &parameters =
1682  typename BaseClass::AdditionalData());
1683 
1687  const std::vector<size_type> &permutation;
1691  const std::vector<size_type> &inverse_permutation;
1696  };
1697 
1709  void
1710  initialize(const MatrixType & A,
1711  const std::vector<size_type> & permutation,
1712  const std::vector<size_type> & inverse_permutation,
1713  const typename BaseClass::AdditionalData &parameters =
1714  typename BaseClass::AdditionalData());
1715 
1726  void
1727  initialize(const MatrixType &A, const AdditionalData &additional_data);
1728 };
1729 
1730 
1731 
1929 template <typename MatrixType = SparseMatrix<double>,
1930  typename VectorType = Vector<double>,
1931  typename PreconditionerType = DiagonalMatrix<VectorType>>
1933 {
1934 public:
1939 
1945  {
1951  {
1957  lanczos,
1967  };
1968 
1972  enum class PolynomialType
1973  {
1977  first_kind,
1982  fourth_kind
1983  };
1984 
1989  const unsigned int degree = 1,
1990  const double smoothing_range = 0.,
1991  const unsigned int eig_cg_n_iterations = 8,
1992  const double eig_cg_residual = 1e-2,
1993  const double max_eigenvalue = 1,
1997 
2001  AdditionalData &
2002  operator=(const AdditionalData &other_data);
2003 
2016  unsigned int degree;
2017 
2030 
2037  unsigned int eig_cg_n_iterations;
2038 
2044 
2051 
2057 
2061  std::shared_ptr<PreconditionerType> preconditioner;
2062 
2067 
2072  };
2073 
2074 
2079 
2091  void
2092  initialize(const MatrixType & matrix,
2093  const AdditionalData &additional_data = AdditionalData());
2094 
2099  void
2100  vmult(VectorType &dst, const VectorType &src) const;
2101 
2106  void
2107  Tvmult(VectorType &dst, const VectorType &src) const;
2108 
2112  void
2113  step(VectorType &dst, const VectorType &src) const;
2114 
2118  void
2119  Tstep(VectorType &dst, const VectorType &src) const;
2120 
2124  void
2126 
2131  size_type
2132  m() const;
2133 
2138  size_type
2139  n() const;
2140 
2146  {
2158  unsigned int cg_iterations;
2163  unsigned int degree;
2168  : min_eigenvalue_estimate{std::numeric_limits<double>::max()}
2169  , max_eigenvalue_estimate{std::numeric_limits<double>::lowest()}
2170  , cg_iterations{0}
2171  , degree{0}
2172  {}
2173  };
2174 
2187  EigenvalueInformation
2188  estimate_eigenvalues(const VectorType &src) const;
2189 
2190 private:
2194  SmartPointer<
2195  const MatrixType,
2198 
2202  mutable VectorType solution_old;
2203 
2207  mutable VectorType temp_vector1;
2208 
2212  mutable VectorType temp_vector2;
2213 
2219 
2223  double theta;
2224 
2229  double delta;
2230 
2236 
2242 };
2243 
2244 
2245 
2247 /* ---------------------------------- Inline functions ------------------- */
2248 
2249 #ifndef DOXYGEN
2250 
2252  : n_rows(0)
2253  , n_columns(0)
2254 {}
2255 
2256 template <typename MatrixType>
2257 inline void
2258 PreconditionIdentity::initialize(const MatrixType &matrix,
2260 {
2261  n_rows = matrix.m();
2262  n_columns = matrix.n();
2263 }
2264 
2265 
2266 template <class VectorType>
2267 inline void
2268 PreconditionIdentity::vmult(VectorType &dst, const VectorType &src) const
2269 {
2270  dst = src;
2271 }
2272 
2273 
2274 
2275 template <class VectorType>
2276 inline void
2277 PreconditionIdentity::Tvmult(VectorType &dst, const VectorType &src) const
2278 {
2279  dst = src;
2280 }
2281 
2282 template <class VectorType>
2283 inline void
2284 PreconditionIdentity::vmult_add(VectorType &dst, const VectorType &src) const
2285 {
2286  dst += src;
2287 }
2288 
2289 
2290 
2291 template <class VectorType>
2292 inline void
2293 PreconditionIdentity::Tvmult_add(VectorType &dst, const VectorType &src) const
2294 {
2295  dst += src;
2296 }
2297 
2298 
2299 
2300 inline void
2302 {}
2303 
2304 
2305 
2308 {
2309  Assert(n_rows != 0, ExcNotInitialized());
2310  return n_rows;
2311 }
2312 
2315 {
2317  return n_columns;
2318 }
2319 
2320 //---------------------------------------------------------------------------
2321 
2323  const double relaxation)
2324  : relaxation(relaxation)
2325 {}
2326 
2327 
2329  : relaxation(0)
2330  , n_rows(0)
2331  , n_columns(0)
2332 {
2333  AdditionalData add_data;
2334  relaxation = add_data.relaxation;
2335 }
2336 
2337 
2338 
2339 inline void
2341  const PreconditionRichardson::AdditionalData &parameters)
2342 {
2343  relaxation = parameters.relaxation;
2344 }
2345 
2346 
2347 
2348 template <typename MatrixType>
2349 inline void
2351  const MatrixType & matrix,
2352  const PreconditionRichardson::AdditionalData &parameters)
2353 {
2354  relaxation = parameters.relaxation;
2355  n_rows = matrix.m();
2356  n_columns = matrix.n();
2357 }
2358 
2359 
2360 
2361 template <class VectorType>
2362 inline void
2363 PreconditionRichardson::vmult(VectorType &dst, const VectorType &src) const
2364 {
2365  static_assert(
2366  std::is_same<size_type, typename VectorType::size_type>::value,
2367  "PreconditionRichardson and VectorType must have the same size_type.");
2368 
2369  dst.equ(relaxation, src);
2370 }
2371 
2372 
2373 
2374 template <class VectorType>
2375 inline void
2376 PreconditionRichardson::Tvmult(VectorType &dst, const VectorType &src) const
2377 {
2378  static_assert(
2379  std::is_same<size_type, typename VectorType::size_type>::value,
2380  "PreconditionRichardson and VectorType must have the same size_type.");
2381 
2382  dst.equ(relaxation, src);
2383 }
2384 
2385 template <class VectorType>
2386 inline void
2387 PreconditionRichardson::vmult_add(VectorType &dst, const VectorType &src) const
2388 {
2389  static_assert(
2390  std::is_same<size_type, typename VectorType::size_type>::value,
2391  "PreconditionRichardson and VectorType must have the same size_type.");
2392 
2393  dst.add(relaxation, src);
2394 }
2395 
2396 
2397 
2398 template <class VectorType>
2399 inline void
2400 PreconditionRichardson::Tvmult_add(VectorType &dst, const VectorType &src) const
2401 {
2402  static_assert(
2403  std::is_same<size_type, typename VectorType::size_type>::value,
2404  "PreconditionRichardson and VectorType must have the same size_type.");
2405 
2406  dst.add(relaxation, src);
2407 }
2408 
2411 {
2412  Assert(n_rows != 0, ExcNotInitialized());
2413  return n_rows;
2414 }
2415 
2418 {
2420  return n_columns;
2421 }
2422 
2423 //---------------------------------------------------------------------------
2424 
2425 template <typename MatrixType, typename PreconditionerType>
2426 inline void
2428  const MatrixType & rA,
2429  const AdditionalData &parameters)
2430 {
2431  A = &rA;
2432  relaxation = parameters.relaxation;
2433 
2434  Assert(parameters.preconditioner, ExcNotInitialized());
2435 
2436  preconditioner = parameters.preconditioner;
2437  n_iterations = parameters.n_iterations;
2438 }
2439 
2440 
2441 template <typename MatrixType, typename PreconditionerType>
2442 inline void
2444 {
2445  A = nullptr;
2446  preconditioner = nullptr;
2447 }
2448 
2449 template <typename MatrixType, typename PreconditionerType>
2450 inline
2453 {
2454  Assert(A != nullptr, ExcNotInitialized());
2455  return A->m();
2456 }
2457 
2458 template <typename MatrixType, typename PreconditionerType>
2459 inline
2462 {
2463  Assert(A != nullptr, ExcNotInitialized());
2464  return A->n();
2465 }
2466 
2467 template <typename MatrixType, typename PreconditionerType>
2468 template <class VectorType>
2469 inline void
2471  VectorType & dst,
2472  const VectorType &src) const
2473 {
2474  Assert(this->A != nullptr, ExcNotInitialized());
2475  Assert(this->preconditioner != nullptr, ExcNotInitialized());
2476 
2477  VectorType tmp1, tmp2;
2478 
2479  for (unsigned int i = 0; i < n_iterations; ++i)
2480  internal::PreconditionRelaxation::step_operations(
2481  *A, *preconditioner, dst, src, relaxation, tmp1, tmp2, i, false);
2482 }
2483 
2484 template <typename MatrixType, typename PreconditionerType>
2485 template <class VectorType>
2486 inline void
2488  VectorType & dst,
2489  const VectorType &src) const
2490 {
2491  Assert(this->A != nullptr, ExcNotInitialized());
2492  Assert(this->preconditioner != nullptr, ExcNotInitialized());
2493 
2494  VectorType tmp1, tmp2;
2495 
2496  for (unsigned int i = 0; i < n_iterations; ++i)
2497  internal::PreconditionRelaxation::step_operations(
2498  *A, *preconditioner, dst, src, relaxation, tmp1, tmp2, i, true);
2499 }
2500 
2501 template <typename MatrixType, typename PreconditionerType>
2502 template <class VectorType>
2503 inline void
2505  VectorType & dst,
2506  const VectorType &src) const
2507 {
2508  Assert(this->A != nullptr, ExcNotInitialized());
2509  Assert(this->preconditioner != nullptr, ExcNotInitialized());
2510 
2511  VectorType tmp1, tmp2;
2512 
2513  for (unsigned int i = 1; i <= n_iterations; ++i)
2514  internal::PreconditionRelaxation::step_operations(
2515  *A, *preconditioner, dst, src, relaxation, tmp1, tmp2, i, false);
2516 }
2517 
2518 template <typename MatrixType, typename PreconditionerType>
2519 template <class VectorType>
2520 inline void
2522  VectorType & dst,
2523  const VectorType &src) const
2524 {
2525  Assert(this->A != nullptr, ExcNotInitialized());
2526  Assert(this->preconditioner != nullptr, ExcNotInitialized());
2527 
2528  VectorType tmp1, tmp2;
2529 
2530  for (unsigned int i = 1; i <= n_iterations; ++i)
2531  internal::PreconditionRelaxation::step_operations(
2532  *A, *preconditioner, dst, src, relaxation, tmp1, tmp2, i, true);
2533 }
2534 
2535 //---------------------------------------------------------------------------
2536 
2537 template <typename MatrixType>
2538 inline void
2540  const AdditionalData &parameters_in)
2541 {
2542  Assert(parameters_in.preconditioner == nullptr, ExcInternalError());
2543 
2544  AdditionalData parameters;
2545  parameters.relaxation = 1.0;
2546  parameters.n_iterations = parameters_in.n_iterations;
2547  parameters.preconditioner =
2548  std::make_shared<PreconditionerType>(A, parameters_in.relaxation);
2549 
2550  this->BaseClass::initialize(A, parameters);
2551 }
2552 
2553 //---------------------------------------------------------------------------
2554 
2555 template <typename MatrixType>
2556 inline void
2557 PreconditionSOR<MatrixType>::initialize(const MatrixType & A,
2558  const AdditionalData &parameters_in)
2559 {
2560  Assert(parameters_in.preconditioner == nullptr, ExcInternalError());
2561 
2562  AdditionalData parameters;
2563  parameters.relaxation = 1.0;
2564  parameters.n_iterations = parameters_in.n_iterations;
2565  parameters.preconditioner =
2566  std::make_shared<PreconditionerType>(A, parameters_in.relaxation);
2567 
2568  this->BaseClass::initialize(A, parameters);
2569 }
2570 
2571 //---------------------------------------------------------------------------
2572 
2573 template <typename MatrixType>
2574 inline void
2575 PreconditionSSOR<MatrixType>::initialize(const MatrixType & A,
2576  const AdditionalData &parameters_in)
2577 {
2578  Assert(parameters_in.preconditioner == nullptr, ExcInternalError());
2579 
2580  AdditionalData parameters;
2581  parameters.relaxation = 1.0;
2582  parameters.n_iterations = parameters_in.n_iterations;
2583  parameters.preconditioner =
2584  std::make_shared<PreconditionerType>(A, parameters_in.relaxation);
2585 
2586  this->BaseClass::initialize(A, parameters);
2587 }
2588 
2589 
2590 
2591 //---------------------------------------------------------------------------
2592 
2593 template <typename MatrixType>
2594 inline void
2596  const MatrixType & A,
2597  const std::vector<size_type> & p,
2598  const std::vector<size_type> & ip,
2599  const typename BaseClass::AdditionalData &parameters_in)
2600 {
2601  Assert(parameters_in.preconditioner == nullptr, ExcInternalError());
2602 
2603  typename BaseClass::AdditionalData parameters;
2604  parameters.relaxation = 1.0;
2605  parameters.n_iterations = parameters_in.n_iterations;
2606  parameters.preconditioner =
2607  std::make_shared<PreconditionerType>(A, parameters_in.relaxation, p, ip);
2608 
2609  this->BaseClass::initialize(A, parameters);
2610 }
2611 
2612 
2613 template <typename MatrixType>
2614 inline void
2615 PreconditionPSOR<MatrixType>::initialize(const MatrixType & A,
2616  const AdditionalData &additional_data)
2617 {
2618  initialize(A,
2619  additional_data.permutation,
2620  additional_data.inverse_permutation,
2621  additional_data.parameters);
2622 }
2623 
2624 template <typename MatrixType>
2626  const std::vector<size_type> &permutation,
2627  const std::vector<size_type> &inverse_permutation,
2629  &parameters)
2630  : permutation(permutation)
2631  , inverse_permutation(inverse_permutation)
2632  , parameters(parameters)
2633 {}
2634 
2635 
2636 //---------------------------------------------------------------------------
2637 
2638 
2639 template <typename MatrixType, class VectorType>
2641  const MatrixType & M,
2642  const function_ptr method)
2643  : matrix(M)
2644  , precondition(method)
2645 {}
2646 
2647 
2648 
2649 template <typename MatrixType, class VectorType>
2650 void
2652  VectorType & dst,
2653  const VectorType &src) const
2654 {
2655  (matrix.*precondition)(dst, src);
2656 }
2657 
2658 //---------------------------------------------------------------------------
2659 
2660 template <typename MatrixType, typename PreconditionerType>
2662  AdditionalData(const double relaxation, const unsigned int n_iterations)
2663  : relaxation(relaxation)
2664  , n_iterations(n_iterations)
2665 {}
2666 
2667 
2668 
2669 //---------------------------------------------------------------------------
2670 
2671 namespace internal
2672 {
2673  namespace PreconditionChebyshevImplementation
2674  {
2675  // for deal.II vectors, perform updates for Chebyshev preconditioner all
2676  // at once to reduce memory transfer. Here, we select between general
2677  // vectors and deal.II vectors where we expand the loop over the (local)
2678  // size of the vector
2679 
2680  // generic part for non-deal.II vectors
2681  template <typename VectorType, typename PreconditionerType>
2682  inline void
2683  vector_updates(const VectorType & rhs,
2684  const PreconditionerType &preconditioner,
2685  const unsigned int iteration_index,
2686  const double factor1,
2687  const double factor2,
2688  VectorType & solution_old,
2689  VectorType & temp_vector1,
2690  VectorType & temp_vector2,
2691  VectorType & solution)
2692  {
2693  if (iteration_index == 0)
2694  {
2695  solution.equ(factor2, rhs);
2696  preconditioner.vmult(solution_old, solution);
2697  }
2698  else if (iteration_index == 1)
2699  {
2700  // compute t = P^{-1} * (b-A*x^{n})
2701  temp_vector1.sadd(-1.0, 1.0, rhs);
2702  preconditioner.vmult(solution_old, temp_vector1);
2703 
2704  // compute x^{n+1} = x^{n} + f_1 * x^{n} + f_2 * t
2705  solution_old.sadd(factor2, 1 + factor1, solution);
2706  }
2707  else
2708  {
2709  // compute t = P^{-1} * (b-A*x^{n})
2710  temp_vector1.sadd(-1.0, 1.0, rhs);
2711  preconditioner.vmult(temp_vector2, temp_vector1);
2712 
2713  // compute x^{n+1} = x^{n} + f_1 * (x^{n}-x^{n-1}) + f_2 * t
2714  solution_old.sadd(-factor1, factor2, temp_vector2);
2715  solution_old.add(1 + factor1, solution);
2716  }
2717 
2718  solution.swap(solution_old);
2719  }
2720 
2721  // generic part for deal.II vectors
2722  template <
2723  typename Number,
2724  typename PreconditionerType,
2725  std::enable_if_t<
2726  !has_vmult_with_std_functions_for_precondition<
2727  PreconditionerType,
2729  int> * = nullptr>
2730  inline void
2731  vector_updates(
2733  const PreconditionerType &preconditioner,
2734  const unsigned int iteration_index,
2735  const double factor1_,
2736  const double factor2_,
2738  &solution_old,
2740  &temp_vector1,
2742  &temp_vector2,
2744  {
2745  const Number factor1 = factor1_;
2746  const Number factor1_plus_1 = 1. + factor1_;
2747  const Number factor2 = factor2_;
2748 
2749  if (iteration_index == 0)
2750  {
2751  const auto solution_old_ptr = solution_old.begin();
2752 
2753  // compute t = P^{-1} * (b)
2754  preconditioner.vmult(solution_old, rhs);
2755 
2756  // compute x^{n+1} = f_2 * t
2758  for (unsigned int i = 0; i < solution_old.locally_owned_size(); ++i)
2759  solution_old_ptr[i] = solution_old_ptr[i] * factor2;
2760  }
2761  else if (iteration_index == 1)
2762  {
2763  const auto solution_ptr = solution.begin();
2764  const auto solution_old_ptr = solution_old.begin();
2765 
2766  // compute t = P^{-1} * (b-A*x^{n})
2767  temp_vector1.sadd(-1.0, 1.0, rhs);
2768 
2769  preconditioner.vmult(solution_old, temp_vector1);
2770 
2771  // compute x^{n+1} = x^{n} + f_1 * x^{n} + f_2 * t
2773  for (unsigned int i = 0; i < solution_old.locally_owned_size(); ++i)
2774  solution_old_ptr[i] =
2775  factor1_plus_1 * solution_ptr[i] + solution_old_ptr[i] * factor2;
2776  }
2777  else
2778  {
2779  const auto solution_ptr = solution.begin();
2780  const auto solution_old_ptr = solution_old.begin();
2781  const auto temp_vector2_ptr = temp_vector2.begin();
2782 
2783  // compute t = P^{-1} * (b-A*x^{n})
2784  temp_vector1.sadd(-1.0, 1.0, rhs);
2785 
2786  preconditioner.vmult(temp_vector2, temp_vector1);
2787 
2788  // compute x^{n+1} = x^{n} + f_1 * (x^{n}-x^{n-1}) + f_2 * t
2790  for (unsigned int i = 0; i < solution_old.locally_owned_size(); ++i)
2791  solution_old_ptr[i] = factor1_plus_1 * solution_ptr[i] -
2792  factor1 * solution_old_ptr[i] +
2793  temp_vector2_ptr[i] * factor2;
2794  }
2795 
2796  solution.swap(solution_old);
2797  }
2798 
2799  template <
2800  typename Number,
2801  typename PreconditionerType,
2802  std::enable_if_t<
2803  has_vmult_with_std_functions_for_precondition<
2804  PreconditionerType,
2806  int> * = nullptr>
2807  inline void
2808  vector_updates(
2810  const PreconditionerType &preconditioner,
2811  const unsigned int iteration_index,
2812  const double factor1_,
2813  const double factor2_,
2815  &solution_old,
2817  &temp_vector1,
2819  &temp_vector2,
2821  {
2822  const Number factor1 = factor1_;
2823  const Number factor1_plus_1 = 1. + factor1_;
2824  const Number factor2 = factor2_;
2825 
2826  const auto rhs_ptr = rhs.begin();
2827  const auto temp_vector1_ptr = temp_vector1.begin();
2828  const auto temp_vector2_ptr = temp_vector2.begin();
2829  const auto solution_ptr = solution.begin();
2830  const auto solution_old_ptr = solution_old.begin();
2831 
2832  if (iteration_index == 0)
2833  {
2834  preconditioner.vmult(
2835  solution,
2836  rhs,
2837  [&](const auto start_range, const auto end_range) {
2838  if (end_range > start_range)
2839  std::memset(solution.begin() + start_range,
2840  0,
2841  sizeof(Number) * (end_range - start_range));
2842  },
2843  [&](const auto begin, const auto end) {
2845  for (std::size_t i = begin; i < end; ++i)
2846  solution_ptr[i] *= factor2;
2847  });
2848  }
2849  else
2850  {
2851  preconditioner.vmult(
2852  temp_vector2,
2853  temp_vector1,
2854  [&](const auto begin, const auto end) {
2855  if (end > begin)
2856  std::memset(temp_vector2.begin() + begin,
2857  0,
2858  sizeof(Number) * (end - begin));
2859 
2861  for (std::size_t i = begin; i < end; ++i)
2862  temp_vector1_ptr[i] = rhs_ptr[i] - temp_vector1_ptr[i];
2863  },
2864  [&](const auto begin, const auto end) {
2865  if (iteration_index == 1)
2866  {
2868  for (std::size_t i = begin; i < end; ++i)
2869  temp_vector2_ptr[i] = factor1_plus_1 * solution_ptr[i] +
2870  factor2 * temp_vector2_ptr[i];
2871  }
2872  else
2873  {
2875  for (std::size_t i = begin; i < end; ++i)
2876  temp_vector2_ptr[i] = factor1_plus_1 * solution_ptr[i] -
2877  factor1 * solution_old_ptr[i] +
2878  factor2 * temp_vector2_ptr[i];
2879  }
2880  });
2881  }
2882 
2883  if (iteration_index > 0)
2884  {
2885  solution_old.swap(temp_vector2);
2886  solution_old.swap(solution);
2887  }
2888  }
2889 
2890  // worker routine for deal.II vectors. Because of vectorization, we need
2891  // to put the loop into an extra structure because the virtual function of
2892  // VectorUpdatesRange prevents the compiler from applying vectorization.
2893  template <typename Number>
2894  struct VectorUpdater
2895  {
2896  VectorUpdater(const Number * rhs,
2897  const Number * matrix_diagonal_inverse,
2898  const unsigned int iteration_index,
2899  const Number factor1,
2900  const Number factor2,
2901  Number * solution_old,
2902  Number * tmp_vector,
2903  Number * solution)
2904  : rhs(rhs)
2905  , matrix_diagonal_inverse(matrix_diagonal_inverse)
2906  , iteration_index(iteration_index)
2907  , factor1(factor1)
2908  , factor2(factor2)
2909  , solution_old(solution_old)
2910  , tmp_vector(tmp_vector)
2911  , solution(solution)
2912  {}
2913 
2914  void
2915  apply_to_subrange(const std::size_t begin, const std::size_t end) const
2916  {
2917  // To circumvent a bug in gcc
2918  // (https://gcc.gnu.org/bugzilla/show_bug.cgi?id=63945), we create
2919  // copies of the variables factor1 and factor2 and do not check based on
2920  // factor1.
2921  const Number factor1 = this->factor1;
2922  const Number factor1_plus_1 = 1. + this->factor1;
2923  const Number factor2 = this->factor2;
2924  if (iteration_index == 0)
2925  {
2927  for (std::size_t i = begin; i < end; ++i)
2928  solution[i] = factor2 * matrix_diagonal_inverse[i] * rhs[i];
2929  }
2930  else if (iteration_index == 1)
2931  {
2932  // x^{n+1} = x^{n} + f_1 * x^{n} + f_2 * P^{-1} * (b-A*x^{n})
2934  for (std::size_t i = begin; i < end; ++i)
2935  // for efficiency reason, write back to temp_vector that is
2936  // already read (avoid read-for-ownership)
2937  tmp_vector[i] =
2938  factor1_plus_1 * solution[i] +
2939  factor2 * matrix_diagonal_inverse[i] * (rhs[i] - tmp_vector[i]);
2940  }
2941  else
2942  {
2943  // x^{n+1} = x^{n} + f_1 * (x^{n}-x^{n-1})
2944  // + f_2 * P^{-1} * (b-A*x^{n})
2946  for (std::size_t i = begin; i < end; ++i)
2947  // for efficiency reason, write back to temp_vector, which is
2948  // already modified during vmult (in best case, the modified
2949  // values are not written back to main memory yet so that
2950  // we do not have to pay additional costs for writing and
2951  // read-for-ownershop)
2952  tmp_vector[i] =
2953  factor1_plus_1 * solution[i] - factor1 * solution_old[i] +
2954  factor2 * matrix_diagonal_inverse[i] * (rhs[i] - tmp_vector[i]);
2955  }
2956  }
2957 
2958  const Number * rhs;
2959  const Number * matrix_diagonal_inverse;
2960  const unsigned int iteration_index;
2961  const Number factor1;
2962  const Number factor2;
2963  mutable Number * solution_old;
2964  mutable Number * tmp_vector;
2965  mutable Number * solution;
2966  };
2967 
2968  template <typename Number>
2969  struct VectorUpdatesRange : public ::parallel::ParallelForInteger
2970  {
2971  VectorUpdatesRange(const VectorUpdater<Number> &updater,
2972  const std::size_t size)
2973  : updater(updater)
2974  {
2976  VectorUpdatesRange::apply_to_subrange(0, size);
2977  else
2978  apply_parallel(
2979  0,
2980  size,
2982  }
2983 
2984  ~VectorUpdatesRange() override = default;
2985 
2986  virtual void
2987  apply_to_subrange(const std::size_t begin,
2988  const std::size_t end) const override
2989  {
2990  updater.apply_to_subrange(begin, end);
2991  }
2992 
2993  const VectorUpdater<Number> &updater;
2994  };
2995 
2996  // selection for diagonal matrix around deal.II vector
2997  template <typename Number>
2998  inline void
2999  vector_updates(
3000  const ::Vector<Number> & rhs,
3001  const ::DiagonalMatrix<::Vector<Number>> &jacobi,
3002  const unsigned int iteration_index,
3003  const double factor1,
3004  const double factor2,
3005  ::Vector<Number> & solution_old,
3006  ::Vector<Number> & temp_vector1,
3007  ::Vector<Number> &,
3008  ::Vector<Number> &solution)
3009  {
3010  VectorUpdater<Number> upd(rhs.begin(),
3011  jacobi.get_vector().begin(),
3012  iteration_index,
3013  factor1,
3014  factor2,
3015  solution_old.begin(),
3016  temp_vector1.begin(),
3017  solution.begin());
3018  VectorUpdatesRange<Number>(upd, rhs.size());
3019 
3020  // swap vectors x^{n+1}->x^{n}, given the updates in the function above
3021  if (iteration_index == 0)
3022  {
3023  // nothing to do here because we can immediately write into the
3024  // solution vector without remembering any of the other vectors
3025  }
3026  else
3027  {
3028  solution.swap(temp_vector1);
3029  solution_old.swap(temp_vector1);
3030  }
3031  }
3032 
3033  // selection for diagonal matrix around parallel deal.II vector
3034  template <typename Number>
3035  inline void
3036  vector_updates(
3038  const ::DiagonalMatrix<
3040  const unsigned int iteration_index,
3041  const double factor1,
3042  const double factor2,
3044  &solution_old,
3046  &temp_vector1,
3049  {
3050  VectorUpdater<Number> upd(rhs.begin(),
3051  jacobi.get_vector().begin(),
3052  iteration_index,
3053  factor1,
3054  factor2,
3055  solution_old.begin(),
3056  temp_vector1.begin(),
3057  solution.begin());
3058  VectorUpdatesRange<Number>(upd, rhs.locally_owned_size());
3059 
3060  // swap vectors x^{n+1}->x^{n}, given the updates in the function above
3061  if (iteration_index == 0)
3062  {
3063  // nothing to do here because we can immediately write into the
3064  // solution vector without remembering any of the other vectors
3065  }
3066  else
3067  {
3068  solution.swap(temp_vector1);
3069  solution_old.swap(temp_vector1);
3070  }
3071  }
3072 
3073  // We need to have a separate declaration for static const members
3074 
3075  // general case and the case that the preconditioner can work on
3076  // ranges (covered by vector_updates())
3077  template <
3078  typename MatrixType,
3079  typename VectorType,
3080  typename PreconditionerType,
3081  std::enable_if_t<
3082  !has_vmult_with_std_functions<MatrixType,
3083  VectorType,
3084  PreconditionerType> &&
3085  !(has_vmult_with_std_functions_for_precondition<PreconditionerType,
3086  VectorType> &&
3087  has_vmult_with_std_functions_for_precondition<MatrixType,
3088  VectorType>),
3089  int> * = nullptr>
3090  inline void
3091  vmult_and_update(const MatrixType & matrix,
3092  const PreconditionerType &preconditioner,
3093  const VectorType & rhs,
3094  const unsigned int iteration_index,
3095  const double factor1,
3096  const double factor2,
3097  VectorType & solution,
3098  VectorType & solution_old,
3099  VectorType & temp_vector1,
3100  VectorType & temp_vector2)
3101  {
3102  if (iteration_index > 0)
3103  matrix.vmult(temp_vector1, solution);
3104  vector_updates(rhs,
3105  preconditioner,
3106  iteration_index,
3107  factor1,
3108  factor2,
3109  solution_old,
3110  temp_vector1,
3111  temp_vector2,
3112  solution);
3113  }
3114 
3115  // case that both the operator and the preconditioner can work on
3116  // subranges
3117  template <
3118  typename MatrixType,
3119  typename VectorType,
3120  typename PreconditionerType,
3121  std::enable_if_t<
3122  !has_vmult_with_std_functions<MatrixType,
3123  VectorType,
3124  PreconditionerType> &&
3125  (has_vmult_with_std_functions_for_precondition<PreconditionerType,
3126  VectorType> &&
3127  has_vmult_with_std_functions_for_precondition<MatrixType,
3128  VectorType>),
3129  int> * = nullptr>
3130  inline void
3131  vmult_and_update(const MatrixType & matrix,
3132  const PreconditionerType &preconditioner,
3133  const VectorType & rhs,
3134  const unsigned int iteration_index,
3135  const double factor1_,
3136  const double factor2_,
3137  VectorType & solution,
3138  VectorType & solution_old,
3139  VectorType & temp_vector1,
3140  VectorType & temp_vector2)
3141  {
3142  using Number = typename VectorType::value_type;
3143 
3144  const Number factor1 = factor1_;
3145  const Number factor1_plus_1 = 1. + factor1_;
3146  const Number factor2 = factor2_;
3147 
3148  if (iteration_index == 0)
3149  {
3150  preconditioner.vmult(
3151  solution,
3152  rhs,
3153  [&](const unsigned int start_range, const unsigned int end_range) {
3154  // zero 'solution' before running the vmult operation
3155  if (end_range > start_range)
3156  std::memset(solution.begin() + start_range,
3157  0,
3158  sizeof(Number) * (end_range - start_range));
3159  },
3160  [&](const unsigned int start_range, const unsigned int end_range) {
3161  const auto solution_ptr = solution.begin();
3162 
3164  for (std::size_t i = start_range; i < end_range; ++i)
3165  solution_ptr[i] *= factor2;
3166  });
3167  }
3168  else
3169  {
3170  temp_vector1.reinit(rhs, true);
3171  temp_vector2.reinit(rhs, true);
3172 
3173  // 1) compute rediduum (including operator application)
3174  matrix.vmult(
3175  temp_vector1,
3176  solution,
3177  [&](const unsigned int start_range, const unsigned int end_range) {
3178  // zero 'temp_vector1' before running the vmult
3179  // operation
3180  if (end_range > start_range)
3181  std::memset(temp_vector1.begin() + start_range,
3182  0,
3183  sizeof(Number) * (end_range - start_range));
3184  },
3185  [&](const unsigned int start_range, const unsigned int end_range) {
3186  const auto rhs_ptr = rhs.begin();
3187  const auto tmp_ptr = temp_vector1.begin();
3188 
3190  for (std::size_t i = start_range; i < end_range; ++i)
3191  tmp_ptr[i] = rhs_ptr[i] - tmp_ptr[i];
3192  });
3193 
3194  // 2) perform vector updates (including preconditioner application)
3195  preconditioner.vmult(
3196  temp_vector2,
3197  temp_vector1,
3198  [&](const unsigned int start_range, const unsigned int end_range) {
3199  // zero 'temp_vector2' before running the vmult
3200  // operation
3201  if (end_range > start_range)
3202  std::memset(temp_vector2.begin() + start_range,
3203  0,
3204  sizeof(Number) * (end_range - start_range));
3205  },
3206  [&](const unsigned int start_range, const unsigned int end_range) {
3207  const auto solution_ptr = solution.begin();
3208  const auto solution_old_ptr = solution_old.begin();
3209  const auto tmp_ptr = temp_vector2.begin();
3210 
3211  if (iteration_index == 1)
3212  {
3214  for (std::size_t i = start_range; i < end_range; ++i)
3215  tmp_ptr[i] =
3216  factor1_plus_1 * solution_ptr[i] + factor2 * tmp_ptr[i];
3217  }
3218  else
3219  {
3221  for (std::size_t i = start_range; i < end_range; ++i)
3222  tmp_ptr[i] = factor1_plus_1 * solution_ptr[i] -
3223  factor1 * solution_old_ptr[i] +
3224  factor2 * tmp_ptr[i];
3225  }
3226  });
3227 
3228  solution.swap(temp_vector2);
3229  solution_old.swap(temp_vector2);
3230  }
3231  }
3232 
3233  // case that the operator can work on subranges and the preconditioner
3234  // is a diagonal
3235  template <typename MatrixType,
3236  typename VectorType,
3237  typename PreconditionerType,
3238  std::enable_if_t<has_vmult_with_std_functions<MatrixType,
3239  VectorType,
3240  PreconditionerType>,
3241  int> * = nullptr>
3242  inline void
3243  vmult_and_update(const MatrixType & matrix,
3244  const PreconditionerType &preconditioner,
3245  const VectorType & rhs,
3246  const unsigned int iteration_index,
3247  const double factor1,
3248  const double factor2,
3249  VectorType & solution,
3250  VectorType & solution_old,
3251  VectorType & temp_vector1,
3252  VectorType &)
3253  {
3254  using Number = typename VectorType::value_type;
3255  VectorUpdater<Number> updater(rhs.begin(),
3256  preconditioner.get_vector().begin(),
3257  iteration_index,
3258  factor1,
3259  factor2,
3260  solution_old.begin(),
3261  temp_vector1.begin(),
3262  solution.begin());
3263  if (iteration_index > 0)
3264  matrix.vmult(
3265  temp_vector1,
3266  solution,
3267  [&](const unsigned int start_range, const unsigned int end_range) {
3268  // zero 'temp_vector1' before running the vmult
3269  // operation
3270  if (end_range > start_range)
3271  std::memset(temp_vector1.begin() + start_range,
3272  0,
3273  sizeof(Number) * (end_range - start_range));
3274  },
3275  [&](const unsigned int start_range, const unsigned int end_range) {
3276  if (end_range > start_range)
3277  updater.apply_to_subrange(start_range, end_range);
3278  });
3279  else
3280  updater.apply_to_subrange(0U, solution.locally_owned_size());
3281 
3282  // swap vectors x^{n+1}->x^{n}, given the updates in the function above
3283  if (iteration_index == 0)
3284  {
3285  // nothing to do here because we can immediately write into the
3286  // solution vector without remembering any of the other vectors
3287  }
3288  else
3289  {
3290  solution.swap(temp_vector1);
3291  solution_old.swap(temp_vector1);
3292  }
3293  }
3294 
3295  template <typename MatrixType, typename PreconditionerType>
3296  inline void
3297  initialize_preconditioner(
3298  const MatrixType & matrix,
3299  std::shared_ptr<PreconditionerType> &preconditioner)
3300  {
3301  (void)matrix;
3302  (void)preconditioner;
3303  AssertThrow(preconditioner.get() != nullptr, ExcNotInitialized());
3304  }
3305 
3306  template <typename MatrixType, typename VectorType>
3307  inline void
3308  initialize_preconditioner(
3309  const MatrixType & matrix,
3310  std::shared_ptr<::DiagonalMatrix<VectorType>> &preconditioner)
3311  {
3312  if (preconditioner.get() == nullptr || preconditioner->m() != matrix.m())
3313  {
3314  if (preconditioner.get() == nullptr)
3315  preconditioner =
3316  std::make_shared<::DiagonalMatrix<VectorType>>();
3317 
3318  Assert(
3319  preconditioner->m() == 0,
3320  ExcMessage(
3321  "Preconditioner appears to be initialized but not sized correctly"));
3322 
3323  // This part only works in serial
3324  if (preconditioner->m() != matrix.m())
3325  {
3326  preconditioner->get_vector().reinit(matrix.m());
3327  for (typename VectorType::size_type i = 0; i < matrix.m(); ++i)
3328  preconditioner->get_vector()(i) = 1. / matrix.el(i, i);
3329  }
3330  }
3331  }
3332 
3333  template <typename VectorType>
3334  void
3335  set_initial_guess(VectorType &vector)
3336  {
3337  vector = 1. / std::sqrt(static_cast<double>(vector.size()));
3338  if (vector.locally_owned_elements().is_element(0))
3339  vector(0) = 0.;
3340  }
3341 
3342  template <typename Number>
3343  void
3344  set_initial_guess(::Vector<Number> &vector)
3345  {
3346  // Choose a high-frequency mode consisting of numbers between 0 and 1
3347  // that is cheap to compute (cheaper than random numbers) but avoids
3348  // obviously re-occurring numbers in multi-component systems by choosing
3349  // a period of 11
3350  for (unsigned int i = 0; i < vector.size(); ++i)
3351  vector(i) = i % 11;
3352 
3353  const Number mean_value = vector.mean_value();
3354  vector.add(-mean_value);
3355  }
3356 
3357  template <typename Number>
3358  void
3359  set_initial_guess(
3361  &vector)
3362  {
3363  // Choose a high-frequency mode consisting of numbers between 0 and 1
3364  // that is cheap to compute (cheaper than random numbers) but avoids
3365  // obviously re-occurring numbers in multi-component systems by choosing
3366  // a period of 11.
3367  // Make initial guess robust with respect to number of processors
3368  // by operating on the global index.
3369  types::global_dof_index first_local_range = 0;
3370  if (!vector.locally_owned_elements().is_empty())
3371  first_local_range = vector.locally_owned_elements().nth_index_in_set(0);
3372  for (unsigned int i = 0; i < vector.locally_owned_size(); ++i)
3373  vector.local_element(i) = (i + first_local_range) % 11;
3374 
3375  const Number mean_value = vector.mean_value();
3376  vector.add(-mean_value);
3377  }
3378 
3379  template <typename Number>
3380  void
3381  set_initial_guess(
3383  {
3384  for (unsigned int block = 0; block < vector.n_blocks(); ++block)
3385  set_initial_guess(vector.block(block));
3386  }
3387 
3388 
3389 # ifdef DEAL_II_WITH_CUDA
3390  template <typename Number>
3391  __global__ void
3392  set_initial_guess_kernel(const types::global_dof_index offset,
3393  const unsigned int locally_owned_size,
3394  Number * values)
3395 
3396  {
3397  const unsigned int index = threadIdx.x + blockDim.x * blockIdx.x;
3398  if (index < locally_owned_size)
3399  values[index] = (index + offset) % 11;
3400  }
3401 
3402  template <typename Number>
3403  void
3404  set_initial_guess(
3406  &vector)
3407  {
3408  // Choose a high-frequency mode consisting of numbers between 0 and 1
3409  // that is cheap to compute (cheaper than random numbers) but avoids
3410  // obviously re-occurring numbers in multi-component systems by choosing
3411  // a period of 11.
3412  // Make initial guess robust with respect to number of processors
3413  // by operating on the global index.
3414  types::global_dof_index first_local_range = 0;
3415  if (!vector.locally_owned_elements().is_empty())
3416  first_local_range = vector.locally_owned_elements().nth_index_in_set(0);
3417 
3418  const auto n_local_elements = vector.locally_owned_size();
3419  const int n_blocks =
3420  1 + (n_local_elements - 1) / CUDAWrappers::block_size;
3421  set_initial_guess_kernel<<<n_blocks, CUDAWrappers::block_size>>>(
3422  first_local_range, n_local_elements, vector.get_values());
3423  AssertCudaKernel();
3424 
3425  const Number mean_value = vector.mean_value();
3426  vector.add(-mean_value);
3427  }
3428 # endif // DEAL_II_WITH_CUDA
3429 
3430  struct EigenvalueTracker
3431  {
3432  public:
3433  void
3434  slot(const std::vector<double> &eigenvalues)
3435  {
3436  values = eigenvalues;
3437  }
3438 
3439  std::vector<double> values;
3440  };
3441 
3442 
3443 
3444  template <typename MatrixType,
3445  typename VectorType,
3446  typename PreconditionerType>
3447  double
3448  power_iteration(const MatrixType & matrix,
3449  VectorType & eigenvector,
3450  const PreconditionerType &preconditioner,
3451  const unsigned int n_iterations)
3452  {
3453  double eigenvalue_estimate = 0.;
3454  eigenvector /= eigenvector.l2_norm();
3455  VectorType vector1, vector2;
3456  vector1.reinit(eigenvector, true);
3457  if (!std::is_same<PreconditionerType, PreconditionIdentity>::value)
3458  vector2.reinit(eigenvector, true);
3459  for (unsigned int i = 0; i < n_iterations; ++i)
3460  {
3461  if (!std::is_same<PreconditionerType, PreconditionIdentity>::value)
3462  {
3463  matrix.vmult(vector2, eigenvector);
3464  preconditioner.vmult(vector1, vector2);
3465  }
3466  else
3467  matrix.vmult(vector1, eigenvector);
3468 
3469  eigenvalue_estimate = eigenvector * vector1;
3470 
3471  vector1 /= vector1.l2_norm();
3472  eigenvector.swap(vector1);
3473  }
3474  return eigenvalue_estimate;
3475  }
3476  } // namespace PreconditionChebyshevImplementation
3477 } // namespace internal
3478 
3479 
3480 
3481 template <typename MatrixType, class VectorType, typename PreconditionerType>
3483  AdditionalData::AdditionalData(const unsigned int degree,
3484  const double smoothing_range,
3485  const unsigned int eig_cg_n_iterations,
3486  const double eig_cg_residual,
3487  const double max_eigenvalue,
3488  const EigenvalueAlgorithm eigenvalue_algorithm,
3489  const PolynomialType polynomial_type)
3490  : degree(degree)
3491  , smoothing_range(smoothing_range)
3492  , eig_cg_n_iterations(eig_cg_n_iterations)
3493  , eig_cg_residual(eig_cg_residual)
3494  , max_eigenvalue(max_eigenvalue)
3495  , eigenvalue_algorithm(eigenvalue_algorithm)
3496  , polynomial_type(polynomial_type)
3497 {}
3498 
3499 
3500 
3501 template <typename MatrixType, class VectorType, typename PreconditionerType>
3502 inline typename PreconditionChebyshev<MatrixType,
3503  VectorType,
3504  PreconditionerType>::AdditionalData &
3506  AdditionalData::operator=(const AdditionalData &other_data)
3507 {
3508  degree = other_data.degree;
3509  smoothing_range = other_data.smoothing_range;
3510  eig_cg_n_iterations = other_data.eig_cg_n_iterations;
3511  eig_cg_residual = other_data.eig_cg_residual;
3512  max_eigenvalue = other_data.max_eigenvalue;
3513  preconditioner = other_data.preconditioner;
3514  eigenvalue_algorithm = other_data.eigenvalue_algorithm;
3515  polynomial_type = other_data.polynomial_type;
3516  constraints.copy_from(other_data.constraints);
3517 
3518  return *this;
3519 }
3520 
3521 
3522 
3523 template <typename MatrixType, typename VectorType, typename PreconditionerType>
3526  : theta(1.)
3527  , delta(1.)
3528  , eigenvalues_are_initialized(false)
3529 {
3530  static_assert(
3531  std::is_same<size_type, typename VectorType::size_type>::value,
3532  "PreconditionChebyshev and VectorType must have the same size_type.");
3533 }
3534 
3535 
3536 
3537 template <typename MatrixType, typename VectorType, typename PreconditionerType>
3538 inline void
3540  const MatrixType & matrix,
3541  const AdditionalData &additional_data)
3542 {
3543  matrix_ptr = &matrix;
3544  data = additional_data;
3545  Assert(data.degree > 0,
3546  ExcMessage("The degree of the Chebyshev method must be positive."));
3547  internal::PreconditionChebyshevImplementation::initialize_preconditioner(
3548  matrix, data.preconditioner);
3549  eigenvalues_are_initialized = false;
3550 }
3551 
3552 
3553 
3554 template <typename MatrixType, typename VectorType, typename PreconditionerType>
3555 inline void
3557 {
3558  eigenvalues_are_initialized = false;
3559  theta = delta = 1.0;
3560  matrix_ptr = nullptr;
3561  {
3562  VectorType empty_vector;
3563  solution_old.reinit(empty_vector);
3564  temp_vector1.reinit(empty_vector);
3565  temp_vector2.reinit(empty_vector);
3566  }
3567  data.preconditioner.reset();
3568 }
3569 
3570 
3571 
3572 template <typename MatrixType, typename VectorType, typename PreconditionerType>
3573 inline typename PreconditionChebyshev<MatrixType,
3574  VectorType,
3575  PreconditionerType>::EigenvalueInformation
3577  estimate_eigenvalues(const VectorType &src) const
3578 {
3579  Assert(eigenvalues_are_initialized == false, ExcInternalError());
3580  Assert(data.preconditioner.get() != nullptr, ExcNotInitialized());
3581 
3583  EigenvalueInformation info{};
3584 
3585  solution_old.reinit(src);
3586  temp_vector1.reinit(src, true);
3587 
3588  if (data.eig_cg_n_iterations > 0)
3589  {
3590  Assert(data.eig_cg_n_iterations > 2,
3591  ExcMessage(
3592  "Need to set at least two iterations to find eigenvalues."));
3593 
3594  internal::PreconditionChebyshevImplementation::EigenvalueTracker
3595  eigenvalue_tracker;
3596 
3597  // set an initial guess that contains some high-frequency parts (to the
3598  // extent possible without knowing the discretization and the numbering)
3599  // to trigger high eigenvalues according to the external function
3600  internal::PreconditionChebyshevImplementation::set_initial_guess(
3601  temp_vector1);
3602  data.constraints.set_zero(temp_vector1);
3603 
3604  if (data.eigenvalue_algorithm ==
3605  AdditionalData::EigenvalueAlgorithm::lanczos)
3606  {
3607  // set a very strict tolerance to force at least two iterations
3608  IterationNumberControl control(data.eig_cg_n_iterations,
3609  1e-10,
3610  false,
3611  false);
3612 
3613  SolverCG<VectorType> solver(control);
3614  solver.connect_eigenvalues_slot(
3615  [&eigenvalue_tracker](const std::vector<double> &eigenvalues) {
3616  eigenvalue_tracker.slot(eigenvalues);
3617  });
3618 
3619  solver.solve(*matrix_ptr,
3620  solution_old,
3621  temp_vector1,
3622  *data.preconditioner);
3623 
3624  info.cg_iterations = control.last_step();
3625  }
3626  else if (data.eigenvalue_algorithm ==
3627  AdditionalData::EigenvalueAlgorithm::power_iteration)
3628  {
3629  Assert(data.degree != numbers::invalid_unsigned_int,
3630  ExcMessage("Cannot estimate the minimal eigenvalue with the "
3631  "power iteration"));
3632 
3633  eigenvalue_tracker.values.push_back(
3634  internal::PreconditionChebyshevImplementation::power_iteration(
3635  *matrix_ptr,
3636  temp_vector1,
3637  *data.preconditioner,
3638  data.eig_cg_n_iterations));
3639  }
3640  else
3641  Assert(false, ExcNotImplemented());
3642 
3643  // read the eigenvalues from the attached eigenvalue tracker
3644  if (eigenvalue_tracker.values.empty())
3645  info.min_eigenvalue_estimate = info.max_eigenvalue_estimate = 1.;
3646  else
3647  {
3648  info.min_eigenvalue_estimate = eigenvalue_tracker.values.front();
3649 
3650  // include a safety factor since the CG method will in general not
3651  // be converged
3652  info.max_eigenvalue_estimate = 1.2 * eigenvalue_tracker.values.back();
3653  }
3654  }
3655  else
3656  {
3657  info.max_eigenvalue_estimate = data.max_eigenvalue;
3658  info.min_eigenvalue_estimate = data.max_eigenvalue / data.smoothing_range;
3659  }
3660 
3661  const double alpha = (data.smoothing_range > 1. ?
3662  info.max_eigenvalue_estimate / data.smoothing_range :
3663  std::min(0.9 * info.max_eigenvalue_estimate,
3664  info.min_eigenvalue_estimate));
3665 
3666  // in case the user set the degree to invalid unsigned int, we have to
3667  // determine the number of necessary iterations from the Chebyshev error
3668  // estimate, given the target tolerance specified by smoothing_range. This
3669  // estimate is based on the error formula given in section 5.1 of
3670  // R. S. Varga, Matrix iterative analysis, 2nd ed., Springer, 2009
3671  if (data.degree == numbers::invalid_unsigned_int)
3672  {
3673  const double actual_range = info.max_eigenvalue_estimate / alpha;
3674  const double sigma = (1. - std::sqrt(1. / actual_range)) /
3675  (1. + std::sqrt(1. / actual_range));
3676  const double eps = data.smoothing_range;
3677  const_cast<
3679  this)
3680  ->data.degree =
3681  1 + static_cast<unsigned int>(
3682  std::log(1. / eps + std::sqrt(1. / eps / eps - 1.)) /
3683  std::log(1. / sigma));
3684  }
3685 
3686  info.degree = data.degree;
3687 
3688  const_cast<
3690  ->delta =
3691  (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind) ?
3692  (info.max_eigenvalue_estimate) :
3693  ((info.max_eigenvalue_estimate - alpha) * 0.5);
3694  const_cast<
3696  ->theta = (info.max_eigenvalue_estimate + alpha) * 0.5;
3697 
3698  // We do not need the second temporary vector in case we have a
3699  // DiagonalMatrix as preconditioner and use deal.II's own vectors
3700  using NumberType = typename VectorType::value_type;
3701  if (std::is_same<PreconditionerType,
3702  ::DiagonalMatrix<VectorType>>::value == false ||
3703  (std::is_same<VectorType, ::Vector<NumberType>>::value == false &&
3704  ((std::is_same<VectorType,
3707  false) ||
3708  (std::is_same<VectorType,
3709  LinearAlgebra::distributed::
3710  Vector<NumberType, MemorySpace::CUDA>>::value ==
3711  false))))
3712  temp_vector2.reinit(src, true);
3713  else
3714  {
3715  VectorType empty_vector;
3716  temp_vector2.reinit(empty_vector);
3717  }
3718 
3719  const_cast<
3721  ->eigenvalues_are_initialized = true;
3722 
3723  return info;
3724 }
3725 
3726 
3727 
3728 template <typename MatrixType, typename VectorType, typename PreconditionerType>
3729 inline void
3731  VectorType & solution,
3732  const VectorType &rhs) const
3733 {
3734  std::lock_guard<std::mutex> lock(mutex);
3735  if (eigenvalues_are_initialized == false)
3736  estimate_eigenvalues(rhs);
3737 
3738  internal::PreconditionChebyshevImplementation::vmult_and_update(
3739  *matrix_ptr,
3740  *data.preconditioner,
3741  rhs,
3742  0,
3743  0.,
3744  (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind) ?
3745  (4. / (3. * delta)) :
3746  (1. / theta),
3747  solution,
3748  solution_old,
3749  temp_vector1,
3750  temp_vector2);
3751 
3752  // if delta is zero, we do not need to iterate because the updates will be
3753  // zero
3754  if (data.degree < 2 || std::abs(delta) < 1e-40)
3755  return;
3756 
3757  double rhok = delta / theta, sigma = theta / delta;
3758  for (unsigned int k = 0; k < data.degree - 1; ++k)
3759  {
3760  double factor1 = 0.0;
3761  double factor2 = 0.0;
3762 
3763  if (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind)
3764  {
3765  factor1 = (2 * k + 1.) / (2 * k + 5.);
3766  factor2 = (8 * k + 12.) / (delta * (2 * k + 5.));
3767  }
3768  else
3769  {
3770  const double rhokp = 1. / (2. * sigma - rhok);
3771  factor1 = rhokp * rhok;
3772  factor2 = 2. * rhokp / delta;
3773  rhok = rhokp;
3774  }
3775 
3776  internal::PreconditionChebyshevImplementation::vmult_and_update(
3777  *matrix_ptr,
3778  *data.preconditioner,
3779  rhs,
3780  k + 1,
3781  factor1,
3782  factor2,
3783  solution,
3784  solution_old,
3785  temp_vector1,
3786  temp_vector2);
3787  }
3788 }
3789 
3790 
3791 
3792 template <typename MatrixType, typename VectorType, typename PreconditionerType>
3793 inline void
3795  VectorType & solution,
3796  const VectorType &rhs) const
3797 {
3798  std::lock_guard<std::mutex> lock(mutex);
3799  if (eigenvalues_are_initialized == false)
3800  estimate_eigenvalues(rhs);
3801 
3802  internal::PreconditionChebyshevImplementation::vector_updates(
3803  rhs,
3804  *data.preconditioner,
3805  0,
3806  0.,
3807  (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind) ?
3808  (4. / (3. * delta)) :
3809  (1. / theta),
3810  solution_old,
3811  temp_vector1,
3812  temp_vector2,
3813  solution);
3814 
3815  if (data.degree < 2 || std::abs(delta) < 1e-40)
3816  return;
3817 
3818  double rhok = delta / theta, sigma = theta / delta;
3819  for (unsigned int k = 0; k < data.degree - 1; ++k)
3820  {
3821  double factor1 = 0.0;
3822  double factor2 = 0.0;
3823 
3824  if (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind)
3825  {
3826  factor1 = (2 * k + 1.) / (2 * k + 5.);
3827  factor2 = (8 * k + 12.) / (delta * (2 * k + 5.));
3828  }
3829  else
3830  {
3831  const double rhokp = 1. / (2. * sigma - rhok);
3832  factor1 = rhokp * rhok;
3833  factor2 = 2. * rhokp / delta;
3834  rhok = rhokp;
3835  }
3836 
3837  matrix_ptr->Tvmult(temp_vector1, solution);
3838  internal::PreconditionChebyshevImplementation::vector_updates(
3839  rhs,
3840  *data.preconditioner,
3841  k + 1,
3842  factor1,
3843  factor2,
3844  solution_old,
3845  temp_vector1,
3846  temp_vector2,
3847  solution);
3848  }
3849 }
3850 
3851 
3852 
3853 template <typename MatrixType, typename VectorType, typename PreconditionerType>
3854 inline void
3856  VectorType & solution,
3857  const VectorType &rhs) const
3858 {
3859  std::lock_guard<std::mutex> lock(mutex);
3860  if (eigenvalues_are_initialized == false)
3861  estimate_eigenvalues(rhs);
3862 
3863  internal::PreconditionChebyshevImplementation::vmult_and_update(
3864  *matrix_ptr,
3865  *data.preconditioner,
3866  rhs,
3867  1,
3868  0.,
3869  (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind) ?
3870  (4. / (3. * delta)) :
3871  (1. / theta),
3872  solution,
3873  solution_old,
3874  temp_vector1,
3875  temp_vector2);
3876 
3877  if (data.degree < 2 || std::abs(delta) < 1e-40)
3878  return;
3879 
3880  double rhok = delta / theta, sigma = theta / delta;
3881  for (unsigned int k = 0; k < data.degree - 1; ++k)
3882  {
3883  double factor1 = 0.0;
3884  double factor2 = 0.0;
3885 
3886  if (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind)
3887  {
3888  factor1 = (2 * k + 1.) / (2 * k + 5.);
3889  factor2 = (8 * k + 12.) / (delta * (2 * k + 5.));
3890  }
3891  else
3892  {
3893  const double rhokp = 1. / (2. * sigma - rhok);
3894  factor1 = rhokp * rhok;
3895  factor2 = 2. * rhokp / delta;
3896  rhok = rhokp;
3897  }
3898 
3899  internal::PreconditionChebyshevImplementation::vmult_and_update(
3900  *matrix_ptr,
3901  *data.preconditioner,
3902  rhs,
3903  k + 2,
3904  factor1,
3905  factor2,
3906  solution,
3907  solution_old,
3908  temp_vector1,
3909  temp_vector2);
3910  }
3911 }
3912 
3913 
3914 
3915 template <typename MatrixType, typename VectorType, typename PreconditionerType>
3916 inline void
3918  VectorType & solution,
3919  const VectorType &rhs) const
3920 {
3921  std::lock_guard<std::mutex> lock(mutex);
3922  if (eigenvalues_are_initialized == false)
3923  estimate_eigenvalues(rhs);
3924 
3925  matrix_ptr->Tvmult(temp_vector1, solution);
3926  internal::PreconditionChebyshevImplementation::vector_updates(
3927  rhs,
3928  *data.preconditioner,
3929  1,
3930  0.,
3931  (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind) ?
3932  (4. / (3. * delta)) :
3933  (1. / theta),
3934  solution_old,
3935  temp_vector1,
3936  temp_vector2,
3937  solution);
3938 
3939  if (data.degree < 2 || std::abs(delta) < 1e-40)
3940  return;
3941 
3942  double rhok = delta / theta, sigma = theta / delta;
3943  for (unsigned int k = 0; k < data.degree - 1; ++k)
3944  {
3945  double factor1 = 0.0;
3946  double factor2 = 0.0;
3947 
3948  if (data.polynomial_type == AdditionalData::PolynomialType::fourth_kind)
3949  {
3950  factor1 = (2 * k + 1.) / (2 * k + 5.);
3951  factor2 = (8 * k + 12.) / (delta * (2 * k + 5.));
3952  }
3953  else
3954  {
3955  const double rhokp = 1. / (2. * sigma - rhok);
3956  factor1 = rhokp * rhok;
3957  factor2 = 2. * rhokp / delta;
3958  rhok = rhokp;
3959  }
3960 
3961  matrix_ptr->Tvmult(temp_vector1, solution);
3962  internal::PreconditionChebyshevImplementation::vector_updates(
3963  rhs,
3964  *data.preconditioner,
3965  k + 2,
3966  factor1,
3967  factor2,
3968  solution_old,
3969  temp_vector1,
3970  temp_vector2,
3971  solution);
3972  }
3973 }
3974 
3975 
3976 
3977 template <typename MatrixType, typename VectorType, typename PreconditionerType>
3978 inline typename PreconditionChebyshev<MatrixType,
3979  VectorType,
3980  PreconditionerType>::size_type
3982 {
3983  Assert(matrix_ptr != nullptr, ExcNotInitialized());
3984  return matrix_ptr->m();
3985 }
3986 
3987 
3988 
3989 template <typename MatrixType, typename VectorType, typename PreconditionerType>
3990 inline typename PreconditionChebyshev<MatrixType,
3991  VectorType,
3992  PreconditionerType>::size_type
3994 {
3995  Assert(matrix_ptr != nullptr, ExcNotInitialized());
3996  return matrix_ptr->n();
3997 }
3998 
3999 #endif // DOXYGEN
4000 
4002 
4003 #endif
unsigned int n_blocks() const
BlockType & block(const unsigned int i)
VectorType & get_vector()
bool is_empty() const
Definition: index_set.h:1788
size_type nth_index_in_set(const size_type local_index) const
Definition: index_set.h:1844
Number local_element(const size_type local_index) const
void swap(Vector< Number, MemorySpace > &v)
size_type locally_owned_size() const
virtual Number mean_value() const override
virtual void sadd(const Number s, const Number a, const VectorSpaceVector< Number > &V) override
virtual void add(const Number a) override
virtual ::IndexSet locally_owned_elements() const override
void Tvmult(VectorType &dst, const VectorType &src) const
size_type m() const
Threads::Mutex mutex
size_type n() const
void step(VectorType &dst, const VectorType &src) const
AdditionalData data
EigenvalueInformation estimate_eigenvalues(const VectorType &src) const
void Tstep(VectorType &dst, const VectorType &src) const
SmartPointer< const MatrixType, PreconditionChebyshev< MatrixType, VectorType, PreconditionerType > > matrix_ptr
types::global_dof_index size_type
void vmult(VectorType &dst, const VectorType &src) const
void initialize(const MatrixType &matrix, const AdditionalData &additional_data=AdditionalData())
void vmult_add(VectorType &, const VectorType &) const
void vmult(VectorType &, const VectorType &) const
size_type m() const
void initialize(const MatrixType &matrix, const AdditionalData &additional_data=AdditionalData())
size_type n() const
void Tvmult(VectorType &, const VectorType &) const
types::global_dof_index size_type
Definition: precondition.h:89
void Tvmult_add(VectorType &, const VectorType &) const
typename BaseClass::AdditionalData AdditionalData
internal::PreconditionRelaxation::PreconditionJacobiImpl< MatrixType > PreconditionerType
void initialize(const MatrixType &A, const AdditionalData &parameters=AdditionalData())
AdditionalData(const std::vector< size_type > &permutation, const std::vector< size_type > &inverse_permutation, const typename BaseClass::AdditionalData &parameters=typename BaseClass::AdditionalData())
BaseClass::AdditionalData parameters
const std::vector< size_type > & inverse_permutation
const std::vector< size_type > & permutation
typename BaseClass::size_type size_type
void initialize(const MatrixType &A, const AdditionalData &additional_data)
internal::PreconditionRelaxation::PreconditionPSORImpl< MatrixType > PreconditionerType
void initialize(const MatrixType &A, const std::vector< size_type > &permutation, const std::vector< size_type > &inverse_permutation, const typename BaseClass::AdditionalData &parameters=typename BaseClass::AdditionalData())
std::shared_ptr< PreconditionerType > preconditioner
Definition: precondition.h:440
AdditionalData(const double relaxation=1., const unsigned int n_iterations=1)
SmartPointer< const MatrixType, PreconditionRelaxation< MatrixType > > A
Definition: precondition.h:505
void Tvmult(VectorType &, const VectorType &) const
std::shared_ptr< PreconditionerType > preconditioner
Definition: precondition.h:520
size_type n() const
size_type m() const
void step(VectorType &x, const VectorType &rhs) const
unsigned int n_iterations
Definition: precondition.h:515
void initialize(const MatrixType &A, const AdditionalData &parameters=AdditionalData())
void Tstep(VectorType &x, const VectorType &rhs) const
void vmult(VectorType &, const VectorType &) const
types::global_dof_index size_type
Definition: precondition.h:412
AdditionalData(const double relaxation=1.)
types::global_dof_index size_type
Definition: precondition.h:204
void vmult_add(VectorType &, const VectorType &) const
size_type n() const
void initialize(const AdditionalData &parameters)
void initialize(const MatrixType &matrix, const AdditionalData &parameters)
size_type m() const
void vmult(VectorType &, const VectorType &) const
void Tvmult(VectorType &, const VectorType &) const
void Tvmult_add(VectorType &, const VectorType &) const
internal::PreconditionRelaxation::PreconditionSORImpl< MatrixType > PreconditionerType
typename BaseClass::AdditionalData AdditionalData
void initialize(const MatrixType &A, const AdditionalData &parameters=AdditionalData())
void initialize(const MatrixType &A, const AdditionalData &parameters=AdditionalData())
typename BaseClass::AdditionalData AdditionalData
internal::PreconditionRelaxation::PreconditionSSORImpl< MatrixType > PreconditionerType
void(MatrixType::*)(VectorType &, const VectorType &) const function_ptr
Definition: precondition.h:369
const function_ptr precondition
Definition: precondition.h:394
const MatrixType & matrix
Definition: precondition.h:389
void vmult(VectorType &dst, const VectorType &src) const
PreconditionUseMatrix(const MatrixType &M, const function_ptr method)
const_iterator end() const
const_iterator begin() const
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)
Definition: vector.h:109
void add(const std::vector< size_type > &indices, const std::vector< OtherNumber > &values)
Number mean_value() const
size_type size() const
virtual void swap(Vector< Number > &v)
iterator begin()
#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 & ExcInternalError()
#define AssertCudaKernel()
Definition: exceptions.h:1967
static ::ExceptionBase & ExcNotInitialized()
#define Assert(cond, exc)
Definition: exceptions.h:1586
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1675
constexpr int block_size
Definition: cuda_size.h:29
@ matrix
Contents is actually a matrix.
@ eigenvalues
Eigenvalue vector is filled.
static const char U
static const char A
types::global_dof_index size_type
Definition: cuda_kernels.h:45
std::enable_if_t< IsBlockVector< VectorType >::value, unsigned int > n_blocks(const VectorType &vector)
Definition: operators.h:50
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
VectorType::value_type * begin(VectorType &V)
VectorType::value_type * end(VectorType &V)
std::array< std::pair< Number, Tensor< 1, dim, Number > >, dim > jacobi(::SymmetricTensor< 2, dim, Number > A)
unsigned int minimum_parallel_grain_size
Definition: parallel.cc:34
static const unsigned int invalid_unsigned_int
Definition: types.h:213
unsigned int global_dof_index
Definition: types.h:82
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
AdditionalData(const unsigned int degree=1, const double smoothing_range=0., const unsigned int eig_cg_n_iterations=8, const double eig_cg_residual=1e-2, const double max_eigenvalue=1, const EigenvalueAlgorithm eigenvalue_algorithm=EigenvalueAlgorithm::lanczos, const PolynomialType polynomial_type=PolynomialType::first_kind)
std::shared_ptr< PreconditionerType > preconditioner
AffineConstraints< double > constraints
EigenvalueAlgorithm eigenvalue_algorithm
AdditionalData & operator=(const AdditionalData &other_data)