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