Reference documentation for deal.II version GIT 85a786e179 2022-11-27 23:15: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  AdditionalData(const unsigned int degree = 1,
1973  const double smoothing_range = 0.,
1974  const unsigned int eig_cg_n_iterations = 8,
1975  const double eig_cg_residual = 1e-2,
1976  const double max_eigenvalue = 1,
1979 
1983  AdditionalData &
1984  operator=(const AdditionalData &other_data);
1985 
1998  unsigned int degree;
1999 
2012 
2019  unsigned int eig_cg_n_iterations;
2020 
2026 
2033 
2039 
2043  std::shared_ptr<PreconditionerType> preconditioner;
2044 
2049  };
2050 
2051 
2056 
2068  void
2069  initialize(const MatrixType & matrix,
2070  const AdditionalData &additional_data = AdditionalData());
2071 
2076  void
2077  vmult(VectorType &dst, const VectorType &src) const;
2078 
2083  void
2084  Tvmult(VectorType &dst, const VectorType &src) const;
2085 
2089  void
2090  step(VectorType &dst, const VectorType &src) const;
2091 
2095  void
2096  Tstep(VectorType &dst, const VectorType &src) const;
2097 
2101  void
2103 
2108  size_type
2109  m() const;
2110 
2115  size_type
2116  n() const;
2117 
2123  {
2135  unsigned int cg_iterations;
2140  unsigned int degree;
2145  : min_eigenvalue_estimate{std::numeric_limits<double>::max()}
2146  , max_eigenvalue_estimate{std::numeric_limits<double>::lowest()}
2147  , cg_iterations{0}
2148  , degree{0}
2149  {}
2150  };
2151 
2164  EigenvalueInformation
2165  estimate_eigenvalues(const VectorType &src) const;
2166 
2167 private:
2171  SmartPointer<
2172  const MatrixType,
2175 
2179  mutable VectorType solution_old;
2180 
2184  mutable VectorType temp_vector1;
2185 
2189  mutable VectorType temp_vector2;
2190 
2196 
2200  double theta;
2201 
2206  double delta;
2207 
2213 
2219 };
2220 
2221 
2222 
2224 /* ---------------------------------- Inline functions ------------------- */
2225 
2226 #ifndef DOXYGEN
2227 
2229  : n_rows(0)
2230  , n_columns(0)
2231 {}
2232 
2233 template <typename MatrixType>
2234 inline void
2235 PreconditionIdentity::initialize(const MatrixType &matrix,
2237 {
2238  n_rows = matrix.m();
2239  n_columns = matrix.n();
2240 }
2241 
2242 
2243 template <class VectorType>
2244 inline void
2245 PreconditionIdentity::vmult(VectorType &dst, const VectorType &src) const
2246 {
2247  dst = src;
2248 }
2249 
2250 
2251 
2252 template <class VectorType>
2253 inline void
2254 PreconditionIdentity::Tvmult(VectorType &dst, const VectorType &src) const
2255 {
2256  dst = src;
2257 }
2258 
2259 template <class VectorType>
2260 inline void
2261 PreconditionIdentity::vmult_add(VectorType &dst, const VectorType &src) const
2262 {
2263  dst += src;
2264 }
2265 
2266 
2267 
2268 template <class VectorType>
2269 inline void
2270 PreconditionIdentity::Tvmult_add(VectorType &dst, const VectorType &src) const
2271 {
2272  dst += src;
2273 }
2274 
2275 
2276 
2277 inline void
2279 {}
2280 
2281 
2282 
2285 {
2286  Assert(n_rows != 0, ExcNotInitialized());
2287  return n_rows;
2288 }
2289 
2292 {
2294  return n_columns;
2295 }
2296 
2297 //---------------------------------------------------------------------------
2298 
2300  const double relaxation)
2301  : relaxation(relaxation)
2302 {}
2303 
2304 
2306  : relaxation(0)
2307  , n_rows(0)
2308  , n_columns(0)
2309 {
2310  AdditionalData add_data;
2311  relaxation = add_data.relaxation;
2312 }
2313 
2314 
2315 
2316 inline void
2318  const PreconditionRichardson::AdditionalData &parameters)
2319 {
2320  relaxation = parameters.relaxation;
2321 }
2322 
2323 
2324 
2325 template <typename MatrixType>
2326 inline void
2328  const MatrixType & matrix,
2329  const PreconditionRichardson::AdditionalData &parameters)
2330 {
2331  relaxation = parameters.relaxation;
2332  n_rows = matrix.m();
2333  n_columns = matrix.n();
2334 }
2335 
2336 
2337 
2338 template <class VectorType>
2339 inline void
2340 PreconditionRichardson::vmult(VectorType &dst, const VectorType &src) const
2341 {
2342  static_assert(
2343  std::is_same<size_type, typename VectorType::size_type>::value,
2344  "PreconditionRichardson and VectorType must have the same size_type.");
2345 
2346  dst.equ(relaxation, src);
2347 }
2348 
2349 
2350 
2351 template <class VectorType>
2352 inline void
2353 PreconditionRichardson::Tvmult(VectorType &dst, const VectorType &src) const
2354 {
2355  static_assert(
2356  std::is_same<size_type, typename VectorType::size_type>::value,
2357  "PreconditionRichardson and VectorType must have the same size_type.");
2358 
2359  dst.equ(relaxation, src);
2360 }
2361 
2362 template <class VectorType>
2363 inline void
2364 PreconditionRichardson::vmult_add(VectorType &dst, const VectorType &src) const
2365 {
2366  static_assert(
2367  std::is_same<size_type, typename VectorType::size_type>::value,
2368  "PreconditionRichardson and VectorType must have the same size_type.");
2369 
2370  dst.add(relaxation, src);
2371 }
2372 
2373 
2374 
2375 template <class VectorType>
2376 inline void
2377 PreconditionRichardson::Tvmult_add(VectorType &dst, const VectorType &src) const
2378 {
2379  static_assert(
2380  std::is_same<size_type, typename VectorType::size_type>::value,
2381  "PreconditionRichardson and VectorType must have the same size_type.");
2382 
2383  dst.add(relaxation, src);
2384 }
2385 
2388 {
2389  Assert(n_rows != 0, ExcNotInitialized());
2390  return n_rows;
2391 }
2392 
2395 {
2397  return n_columns;
2398 }
2399 
2400 //---------------------------------------------------------------------------
2401 
2402 template <typename MatrixType, typename PreconditionerType>
2403 inline void
2405  const MatrixType & rA,
2406  const AdditionalData &parameters)
2407 {
2408  A = &rA;
2409  relaxation = parameters.relaxation;
2410 
2411  Assert(parameters.preconditioner, ExcNotInitialized());
2412 
2413  preconditioner = parameters.preconditioner;
2414  n_iterations = parameters.n_iterations;
2415 }
2416 
2417 
2418 template <typename MatrixType, typename PreconditionerType>
2419 inline void
2421 {
2422  A = nullptr;
2423  preconditioner = nullptr;
2424 }
2425 
2426 template <typename MatrixType, typename PreconditionerType>
2427 inline
2430 {
2431  Assert(A != nullptr, ExcNotInitialized());
2432  return A->m();
2433 }
2434 
2435 template <typename MatrixType, typename PreconditionerType>
2436 inline
2439 {
2440  Assert(A != nullptr, ExcNotInitialized());
2441  return A->n();
2442 }
2443 
2444 template <typename MatrixType, typename PreconditionerType>
2445 template <class VectorType>
2446 inline void
2448  VectorType & dst,
2449  const VectorType &src) const
2450 {
2451  Assert(this->A != nullptr, ExcNotInitialized());
2452  Assert(this->preconditioner != nullptr, ExcNotInitialized());
2453 
2454  VectorType tmp1, tmp2;
2455 
2456  for (unsigned int i = 0; i < n_iterations; ++i)
2457  internal::PreconditionRelaxation::step_operations(
2458  *A, *preconditioner, dst, src, relaxation, tmp1, tmp2, i, false);
2459 }
2460 
2461 template <typename MatrixType, typename PreconditionerType>
2462 template <class VectorType>
2463 inline void
2465  VectorType & dst,
2466  const VectorType &src) const
2467 {
2468  Assert(this->A != nullptr, ExcNotInitialized());
2469  Assert(this->preconditioner != nullptr, ExcNotInitialized());
2470 
2471  VectorType tmp1, tmp2;
2472 
2473  for (unsigned int i = 0; i < n_iterations; ++i)
2474  internal::PreconditionRelaxation::step_operations(
2475  *A, *preconditioner, dst, src, relaxation, tmp1, tmp2, i, true);
2476 }
2477 
2478 template <typename MatrixType, typename PreconditionerType>
2479 template <class VectorType>
2480 inline void
2482  VectorType & dst,
2483  const VectorType &src) const
2484 {
2485  Assert(this->A != nullptr, ExcNotInitialized());
2486  Assert(this->preconditioner != nullptr, ExcNotInitialized());
2487 
2488  VectorType tmp1, tmp2;
2489 
2490  for (unsigned int i = 1; i <= n_iterations; ++i)
2491  internal::PreconditionRelaxation::step_operations(
2492  *A, *preconditioner, dst, src, relaxation, tmp1, tmp2, i, false);
2493 }
2494 
2495 template <typename MatrixType, typename PreconditionerType>
2496 template <class VectorType>
2497 inline void
2499  VectorType & dst,
2500  const VectorType &src) const
2501 {
2502  Assert(this->A != nullptr, ExcNotInitialized());
2503  Assert(this->preconditioner != nullptr, ExcNotInitialized());
2504 
2505  VectorType tmp1, tmp2;
2506 
2507  for (unsigned int i = 1; i <= n_iterations; ++i)
2508  internal::PreconditionRelaxation::step_operations(
2509  *A, *preconditioner, dst, src, relaxation, tmp1, tmp2, i, true);
2510 }
2511 
2512 //---------------------------------------------------------------------------
2513 
2514 template <typename MatrixType>
2515 inline void
2517  const AdditionalData &parameters_in)
2518 {
2519  Assert(parameters_in.preconditioner == nullptr, ExcInternalError());
2520 
2521  AdditionalData parameters;
2522  parameters.relaxation = 1.0;
2523  parameters.n_iterations = parameters_in.n_iterations;
2524  parameters.preconditioner =
2525  std::make_shared<PreconditionerType>(A, parameters_in.relaxation);
2526 
2527  this->BaseClass::initialize(A, parameters);
2528 }
2529 
2530 //---------------------------------------------------------------------------
2531 
2532 template <typename MatrixType>
2533 inline void
2534 PreconditionSOR<MatrixType>::initialize(const MatrixType & A,
2535  const AdditionalData &parameters_in)
2536 {
2537  Assert(parameters_in.preconditioner == nullptr, ExcInternalError());
2538 
2539  AdditionalData parameters;
2540  parameters.relaxation = 1.0;
2541  parameters.n_iterations = parameters_in.n_iterations;
2542  parameters.preconditioner =
2543  std::make_shared<PreconditionerType>(A, parameters_in.relaxation);
2544 
2545  this->BaseClass::initialize(A, parameters);
2546 }
2547 
2548 //---------------------------------------------------------------------------
2549 
2550 template <typename MatrixType>
2551 inline void
2552 PreconditionSSOR<MatrixType>::initialize(const MatrixType & A,
2553  const AdditionalData &parameters_in)
2554 {
2555  Assert(parameters_in.preconditioner == nullptr, ExcInternalError());
2556 
2557  AdditionalData parameters;
2558  parameters.relaxation = 1.0;
2559  parameters.n_iterations = parameters_in.n_iterations;
2560  parameters.preconditioner =
2561  std::make_shared<PreconditionerType>(A, parameters_in.relaxation);
2562 
2563  this->BaseClass::initialize(A, parameters);
2564 }
2565 
2566 
2567 
2568 //---------------------------------------------------------------------------
2569 
2570 template <typename MatrixType>
2571 inline void
2573  const MatrixType & A,
2574  const std::vector<size_type> & p,
2575  const std::vector<size_type> & ip,
2576  const typename BaseClass::AdditionalData &parameters_in)
2577 {
2578  Assert(parameters_in.preconditioner == nullptr, ExcInternalError());
2579 
2580  typename BaseClass::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, p, ip);
2585 
2586  this->BaseClass::initialize(A, parameters);
2587 }
2588 
2589 
2590 template <typename MatrixType>
2591 inline void
2592 PreconditionPSOR<MatrixType>::initialize(const MatrixType & A,
2593  const AdditionalData &additional_data)
2594 {
2595  initialize(A,
2596  additional_data.permutation,
2597  additional_data.inverse_permutation,
2598  additional_data.parameters);
2599 }
2600 
2601 template <typename MatrixType>
2603  const std::vector<size_type> &permutation,
2604  const std::vector<size_type> &inverse_permutation,
2606  &parameters)
2607  : permutation(permutation)
2608  , inverse_permutation(inverse_permutation)
2609  , parameters(parameters)
2610 {}
2611 
2612 
2613 //---------------------------------------------------------------------------
2614 
2615 
2616 template <typename MatrixType, class VectorType>
2618  const MatrixType & M,
2619  const function_ptr method)
2620  : matrix(M)
2621  , precondition(method)
2622 {}
2623 
2624 
2625 
2626 template <typename MatrixType, class VectorType>
2627 void
2629  VectorType & dst,
2630  const VectorType &src) const
2631 {
2632  (matrix.*precondition)(dst, src);
2633 }
2634 
2635 //---------------------------------------------------------------------------
2636 
2637 template <typename MatrixType, typename PreconditionerType>
2639  AdditionalData(const double relaxation, const unsigned int n_iterations)
2640  : relaxation(relaxation)
2641  , n_iterations(n_iterations)
2642 {}
2643 
2644 
2645 
2646 //---------------------------------------------------------------------------
2647 
2648 namespace internal
2649 {
2650  namespace PreconditionChebyshevImplementation
2651  {
2652  // for deal.II vectors, perform updates for Chebyshev preconditioner all
2653  // at once to reduce memory transfer. Here, we select between general
2654  // vectors and deal.II vectors where we expand the loop over the (local)
2655  // size of the vector
2656 
2657  // generic part for non-deal.II vectors
2658  template <typename VectorType, typename PreconditionerType>
2659  inline void
2660  vector_updates(const VectorType & rhs,
2661  const PreconditionerType &preconditioner,
2662  const unsigned int iteration_index,
2663  const double factor1,
2664  const double factor2,
2665  VectorType & solution_old,
2666  VectorType & temp_vector1,
2667  VectorType & temp_vector2,
2668  VectorType & solution)
2669  {
2670  if (iteration_index == 0)
2671  {
2672  solution.equ(factor2, rhs);
2673  preconditioner.vmult(solution_old, solution);
2674  }
2675  else if (iteration_index == 1)
2676  {
2677  // compute t = P^{-1} * (b-A*x^{n})
2678  temp_vector1.sadd(-1.0, 1.0, rhs);
2679  preconditioner.vmult(solution_old, temp_vector1);
2680 
2681  // compute x^{n+1} = x^{n} + f_1 * x^{n} + f_2 * t
2682  solution_old.sadd(factor2, 1 + factor1, solution);
2683  }
2684  else
2685  {
2686  // compute t = P^{-1} * (b-A*x^{n})
2687  temp_vector1.sadd(-1.0, 1.0, rhs);
2688  preconditioner.vmult(temp_vector2, temp_vector1);
2689 
2690  // compute x^{n+1} = x^{n} + f_1 * (x^{n}-x^{n-1}) + f_2 * t
2691  solution_old.sadd(-factor1, factor2, temp_vector2);
2692  solution_old.add(1 + factor1, solution);
2693  }
2694 
2695  solution.swap(solution_old);
2696  }
2697 
2698  // generic part for deal.II vectors
2699  template <
2700  typename Number,
2701  typename PreconditionerType,
2702  std::enable_if_t<
2703  !has_vmult_with_std_functions_for_precondition<
2704  PreconditionerType,
2706  int> * = nullptr>
2707  inline void
2708  vector_updates(
2710  const PreconditionerType &preconditioner,
2711  const unsigned int iteration_index,
2712  const double factor1_,
2713  const double factor2_,
2715  &solution_old,
2717  &temp_vector1,
2719  &temp_vector2,
2721  {
2722  const Number factor1 = factor1_;
2723  const Number factor1_plus_1 = 1. + factor1_;
2724  const Number factor2 = factor2_;
2725 
2726  if (iteration_index == 0)
2727  {
2728  const auto solution_old_ptr = solution_old.begin();
2729 
2730  // compute t = P^{-1} * (b)
2731  preconditioner.vmult(solution_old, rhs);
2732 
2733  // compute x^{n+1} = f_2 * t
2735  for (unsigned int i = 0; i < solution_old.locally_owned_size(); ++i)
2736  solution_old_ptr[i] = solution_old_ptr[i] * factor2;
2737  }
2738  else if (iteration_index == 1)
2739  {
2740  const auto solution_ptr = solution.begin();
2741  const auto solution_old_ptr = solution_old.begin();
2742 
2743  // compute t = P^{-1} * (b-A*x^{n})
2744  temp_vector1.sadd(-1.0, 1.0, rhs);
2745 
2746  preconditioner.vmult(solution_old, temp_vector1);
2747 
2748  // compute x^{n+1} = x^{n} + f_1 * x^{n} + f_2 * t
2750  for (unsigned int i = 0; i < solution_old.locally_owned_size(); ++i)
2751  solution_old_ptr[i] =
2752  factor1_plus_1 * solution_ptr[i] + solution_old_ptr[i] * factor2;
2753  }
2754  else
2755  {
2756  const auto solution_ptr = solution.begin();
2757  const auto solution_old_ptr = solution_old.begin();
2758  const auto temp_vector2_ptr = temp_vector2.begin();
2759 
2760  // compute t = P^{-1} * (b-A*x^{n})
2761  temp_vector1.sadd(-1.0, 1.0, rhs);
2762 
2763  preconditioner.vmult(temp_vector2, temp_vector1);
2764 
2765  // compute x^{n+1} = x^{n} + f_1 * (x^{n}-x^{n-1}) + f_2 * t
2767  for (unsigned int i = 0; i < solution_old.locally_owned_size(); ++i)
2768  solution_old_ptr[i] = factor1_plus_1 * solution_ptr[i] -
2769  factor1 * solution_old_ptr[i] +
2770  temp_vector2_ptr[i] * factor2;
2771  }
2772 
2773  solution.swap(solution_old);
2774  }
2775 
2776  template <
2777  typename Number,
2778  typename PreconditionerType,
2779  std::enable_if_t<
2780  has_vmult_with_std_functions_for_precondition<
2781  PreconditionerType,
2783  int> * = nullptr>
2784  inline void
2785  vector_updates(
2787  const PreconditionerType &preconditioner,
2788  const unsigned int iteration_index,
2789  const double factor1_,
2790  const double factor2_,
2792  &solution_old,
2794  &temp_vector1,
2796  &temp_vector2,
2798  {
2799  const Number factor1 = factor1_;
2800  const Number factor1_plus_1 = 1. + factor1_;
2801  const Number factor2 = factor2_;
2802 
2803  const auto rhs_ptr = rhs.begin();
2804  const auto temp_vector1_ptr = temp_vector1.begin();
2805  const auto temp_vector2_ptr = temp_vector2.begin();
2806  const auto solution_ptr = solution.begin();
2807  const auto solution_old_ptr = solution_old.begin();
2808 
2809  if (iteration_index == 0)
2810  {
2811  preconditioner.vmult(
2812  solution,
2813  rhs,
2814  [&](const auto start_range, const auto end_range) {
2815  if (end_range > start_range)
2816  std::memset(solution.begin() + start_range,
2817  0,
2818  sizeof(Number) * (end_range - start_range));
2819  },
2820  [&](const auto begin, const auto end) {
2822  for (std::size_t i = begin; i < end; ++i)
2823  solution_ptr[i] *= factor2;
2824  });
2825  }
2826  else
2827  {
2828  preconditioner.vmult(
2829  temp_vector2,
2830  temp_vector1,
2831  [&](const auto begin, const auto end) {
2832  if (end > begin)
2833  std::memset(temp_vector2.begin() + begin,
2834  0,
2835  sizeof(Number) * (end - begin));
2836 
2838  for (std::size_t i = begin; i < end; ++i)
2839  temp_vector1_ptr[i] = rhs_ptr[i] - temp_vector1_ptr[i];
2840  },
2841  [&](const auto begin, const auto end) {
2842  if (iteration_index == 1)
2843  {
2845  for (std::size_t i = begin; i < end; ++i)
2846  temp_vector2_ptr[i] = factor1_plus_1 * solution_ptr[i] +
2847  factor2 * temp_vector2_ptr[i];
2848  }
2849  else
2850  {
2852  for (std::size_t i = begin; i < end; ++i)
2853  temp_vector2_ptr[i] = factor1_plus_1 * solution_ptr[i] -
2854  factor1 * solution_old_ptr[i] +
2855  factor2 * temp_vector2_ptr[i];
2856  }
2857  });
2858  }
2859 
2860  if (iteration_index > 0)
2861  {
2862  solution_old.swap(temp_vector2);
2863  solution_old.swap(solution);
2864  }
2865  }
2866 
2867  // worker routine for deal.II vectors. Because of vectorization, we need
2868  // to put the loop into an extra structure because the virtual function of
2869  // VectorUpdatesRange prevents the compiler from applying vectorization.
2870  template <typename Number>
2871  struct VectorUpdater
2872  {
2873  VectorUpdater(const Number * rhs,
2874  const Number * matrix_diagonal_inverse,
2875  const unsigned int iteration_index,
2876  const Number factor1,
2877  const Number factor2,
2878  Number * solution_old,
2879  Number * tmp_vector,
2880  Number * solution)
2881  : rhs(rhs)
2882  , matrix_diagonal_inverse(matrix_diagonal_inverse)
2883  , iteration_index(iteration_index)
2884  , factor1(factor1)
2885  , factor2(factor2)
2886  , solution_old(solution_old)
2887  , tmp_vector(tmp_vector)
2888  , solution(solution)
2889  {}
2890 
2891  void
2892  apply_to_subrange(const std::size_t begin, const std::size_t end) const
2893  {
2894  // To circumvent a bug in gcc
2895  // (https://gcc.gnu.org/bugzilla/show_bug.cgi?id=63945), we create
2896  // copies of the variables factor1 and factor2 and do not check based on
2897  // factor1.
2898  const Number factor1 = this->factor1;
2899  const Number factor1_plus_1 = 1. + this->factor1;
2900  const Number factor2 = this->factor2;
2901  if (iteration_index == 0)
2902  {
2904  for (std::size_t i = begin; i < end; ++i)
2905  solution[i] = factor2 * matrix_diagonal_inverse[i] * rhs[i];
2906  }
2907  else if (iteration_index == 1)
2908  {
2909  // x^{n+1} = x^{n} + f_1 * x^{n} + f_2 * P^{-1} * (b-A*x^{n})
2911  for (std::size_t i = begin; i < end; ++i)
2912  // for efficiency reason, write back to temp_vector that is
2913  // already read (avoid read-for-ownership)
2914  tmp_vector[i] =
2915  factor1_plus_1 * solution[i] +
2916  factor2 * matrix_diagonal_inverse[i] * (rhs[i] - tmp_vector[i]);
2917  }
2918  else
2919  {
2920  // x^{n+1} = x^{n} + f_1 * (x^{n}-x^{n-1})
2921  // + f_2 * P^{-1} * (b-A*x^{n})
2923  for (std::size_t i = begin; i < end; ++i)
2924  // for efficiency reason, write back to temp_vector, which is
2925  // already modified during vmult (in best case, the modified
2926  // values are not written back to main memory yet so that
2927  // we do not have to pay additional costs for writing and
2928  // read-for-ownershop)
2929  tmp_vector[i] =
2930  factor1_plus_1 * solution[i] - factor1 * solution_old[i] +
2931  factor2 * matrix_diagonal_inverse[i] * (rhs[i] - tmp_vector[i]);
2932  }
2933  }
2934 
2935  const Number * rhs;
2936  const Number * matrix_diagonal_inverse;
2937  const unsigned int iteration_index;
2938  const Number factor1;
2939  const Number factor2;
2940  mutable Number * solution_old;
2941  mutable Number * tmp_vector;
2942  mutable Number * solution;
2943  };
2944 
2945  template <typename Number>
2946  struct VectorUpdatesRange : public ::parallel::ParallelForInteger
2947  {
2948  VectorUpdatesRange(const VectorUpdater<Number> &updater,
2949  const std::size_t size)
2950  : updater(updater)
2951  {
2953  VectorUpdatesRange::apply_to_subrange(0, size);
2954  else
2955  apply_parallel(
2956  0,
2957  size,
2959  }
2960 
2961  ~VectorUpdatesRange() override = default;
2962 
2963  virtual void
2964  apply_to_subrange(const std::size_t begin,
2965  const std::size_t end) const override
2966  {
2967  updater.apply_to_subrange(begin, end);
2968  }
2969 
2970  const VectorUpdater<Number> &updater;
2971  };
2972 
2973  // selection for diagonal matrix around deal.II vector
2974  template <typename Number>
2975  inline void
2976  vector_updates(
2977  const ::Vector<Number> & rhs,
2978  const ::DiagonalMatrix<::Vector<Number>> &jacobi,
2979  const unsigned int iteration_index,
2980  const double factor1,
2981  const double factor2,
2982  ::Vector<Number> & solution_old,
2983  ::Vector<Number> & temp_vector1,
2984  ::Vector<Number> &,
2985  ::Vector<Number> &solution)
2986  {
2987  VectorUpdater<Number> upd(rhs.begin(),
2988  jacobi.get_vector().begin(),
2989  iteration_index,
2990  factor1,
2991  factor2,
2992  solution_old.begin(),
2993  temp_vector1.begin(),
2994  solution.begin());
2995  VectorUpdatesRange<Number>(upd, rhs.size());
2996 
2997  // swap vectors x^{n+1}->x^{n}, given the updates in the function above
2998  if (iteration_index == 0)
2999  {
3000  // nothing to do here because we can immediately write into the
3001  // solution vector without remembering any of the other vectors
3002  }
3003  else
3004  {
3005  solution.swap(temp_vector1);
3006  solution_old.swap(temp_vector1);
3007  }
3008  }
3009 
3010  // selection for diagonal matrix around parallel deal.II vector
3011  template <typename Number>
3012  inline void
3013  vector_updates(
3015  const ::DiagonalMatrix<
3017  const unsigned int iteration_index,
3018  const double factor1,
3019  const double factor2,
3021  &solution_old,
3023  &temp_vector1,
3026  {
3027  VectorUpdater<Number> upd(rhs.begin(),
3028  jacobi.get_vector().begin(),
3029  iteration_index,
3030  factor1,
3031  factor2,
3032  solution_old.begin(),
3033  temp_vector1.begin(),
3034  solution.begin());
3035  VectorUpdatesRange<Number>(upd, rhs.locally_owned_size());
3036 
3037  // swap vectors x^{n+1}->x^{n}, given the updates in the function above
3038  if (iteration_index == 0)
3039  {
3040  // nothing to do here because we can immediately write into the
3041  // solution vector without remembering any of the other vectors
3042  }
3043  else
3044  {
3045  solution.swap(temp_vector1);
3046  solution_old.swap(temp_vector1);
3047  }
3048  }
3049 
3050  // We need to have a separate declaration for static const members
3051 
3052  // general case and the case that the preconditioner can work on
3053  // ranges (covered by vector_updates())
3054  template <
3055  typename MatrixType,
3056  typename VectorType,
3057  typename PreconditionerType,
3058  std::enable_if_t<
3059  !has_vmult_with_std_functions<MatrixType,
3060  VectorType,
3061  PreconditionerType> &&
3062  !(has_vmult_with_std_functions_for_precondition<PreconditionerType,
3063  VectorType> &&
3064  has_vmult_with_std_functions_for_precondition<MatrixType,
3065  VectorType>),
3066  int> * = nullptr>
3067  inline void
3068  vmult_and_update(const MatrixType & matrix,
3069  const PreconditionerType &preconditioner,
3070  const VectorType & rhs,
3071  const unsigned int iteration_index,
3072  const double factor1,
3073  const double factor2,
3074  VectorType & solution,
3075  VectorType & solution_old,
3076  VectorType & temp_vector1,
3077  VectorType & temp_vector2)
3078  {
3079  if (iteration_index > 0)
3080  matrix.vmult(temp_vector1, solution);
3081  vector_updates(rhs,
3082  preconditioner,
3083  iteration_index,
3084  factor1,
3085  factor2,
3086  solution_old,
3087  temp_vector1,
3088  temp_vector2,
3089  solution);
3090  }
3091 
3092  // case that both the operator and the preconditioner can work on
3093  // subranges
3094  template <
3095  typename MatrixType,
3096  typename VectorType,
3097  typename PreconditionerType,
3098  std::enable_if_t<
3099  !has_vmult_with_std_functions<MatrixType,
3100  VectorType,
3101  PreconditionerType> &&
3102  (has_vmult_with_std_functions_for_precondition<PreconditionerType,
3103  VectorType> &&
3104  has_vmult_with_std_functions_for_precondition<MatrixType,
3105  VectorType>),
3106  int> * = nullptr>
3107  inline void
3108  vmult_and_update(const MatrixType & matrix,
3109  const PreconditionerType &preconditioner,
3110  const VectorType & rhs,
3111  const unsigned int iteration_index,
3112  const double factor1_,
3113  const double factor2_,
3114  VectorType & solution,
3115  VectorType & solution_old,
3116  VectorType & temp_vector1,
3117  VectorType & temp_vector2)
3118  {
3119  using Number = typename VectorType::value_type;
3120 
3121  const Number factor1 = factor1_;
3122  const Number factor1_plus_1 = 1. + factor1_;
3123  const Number factor2 = factor2_;
3124 
3125  if (iteration_index == 0)
3126  {
3127  preconditioner.vmult(
3128  solution,
3129  rhs,
3130  [&](const unsigned int start_range, const unsigned int end_range) {
3131  // zero 'solution' before running the vmult operation
3132  if (end_range > start_range)
3133  std::memset(solution.begin() + start_range,
3134  0,
3135  sizeof(Number) * (end_range - start_range));
3136  },
3137  [&](const unsigned int start_range, const unsigned int end_range) {
3138  const auto solution_ptr = solution.begin();
3139 
3141  for (std::size_t i = start_range; i < end_range; ++i)
3142  solution_ptr[i] *= factor2;
3143  });
3144  }
3145  else
3146  {
3147  temp_vector1.reinit(rhs, true);
3148  temp_vector2.reinit(rhs, true);
3149 
3150  // 1) compute rediduum (including operator application)
3151  matrix.vmult(
3152  temp_vector1,
3153  solution,
3154  [&](const unsigned int start_range, const unsigned int end_range) {
3155  // zero 'temp_vector1' before running the vmult
3156  // operation
3157  if (end_range > start_range)
3158  std::memset(temp_vector1.begin() + start_range,
3159  0,
3160  sizeof(Number) * (end_range - start_range));
3161  },
3162  [&](const unsigned int start_range, const unsigned int end_range) {
3163  const auto rhs_ptr = rhs.begin();
3164  const auto tmp_ptr = temp_vector1.begin();
3165 
3167  for (std::size_t i = start_range; i < end_range; ++i)
3168  tmp_ptr[i] = rhs_ptr[i] - tmp_ptr[i];
3169  });
3170 
3171  // 2) perform vector updates (including preconditioner application)
3172  preconditioner.vmult(
3173  temp_vector2,
3174  temp_vector1,
3175  [&](const unsigned int start_range, const unsigned int end_range) {
3176  // zero 'temp_vector2' before running the vmult
3177  // operation
3178  if (end_range > start_range)
3179  std::memset(temp_vector2.begin() + start_range,
3180  0,
3181  sizeof(Number) * (end_range - start_range));
3182  },
3183  [&](const unsigned int start_range, const unsigned int end_range) {
3184  const auto solution_ptr = solution.begin();
3185  const auto solution_old_ptr = solution_old.begin();
3186  const auto tmp_ptr = temp_vector2.begin();
3187 
3188  if (iteration_index == 1)
3189  {
3191  for (std::size_t i = start_range; i < end_range; ++i)
3192  tmp_ptr[i] =
3193  factor1_plus_1 * solution_ptr[i] + factor2 * tmp_ptr[i];
3194  }
3195  else
3196  {
3198  for (std::size_t i = start_range; i < end_range; ++i)
3199  tmp_ptr[i] = factor1_plus_1 * solution_ptr[i] -
3200  factor1 * solution_old_ptr[i] +
3201  factor2 * tmp_ptr[i];
3202  }
3203  });
3204 
3205  solution.swap(temp_vector2);
3206  solution_old.swap(temp_vector2);
3207  }
3208  }
3209 
3210  // case that the operator can work on subranges and the preconditioner
3211  // is a diagonal
3212  template <typename MatrixType,
3213  typename VectorType,
3214  typename PreconditionerType,
3215  std::enable_if_t<has_vmult_with_std_functions<MatrixType,
3216  VectorType,
3217  PreconditionerType>,
3218  int> * = nullptr>
3219  inline void
3220  vmult_and_update(const MatrixType & matrix,
3221  const PreconditionerType &preconditioner,
3222  const VectorType & rhs,
3223  const unsigned int iteration_index,
3224  const double factor1,
3225  const double factor2,
3226  VectorType & solution,
3227  VectorType & solution_old,
3228  VectorType & temp_vector1,
3229  VectorType &)
3230  {
3231  using Number = typename VectorType::value_type;
3232  VectorUpdater<Number> updater(rhs.begin(),
3233  preconditioner.get_vector().begin(),
3234  iteration_index,
3235  factor1,
3236  factor2,
3237  solution_old.begin(),
3238  temp_vector1.begin(),
3239  solution.begin());
3240  if (iteration_index > 0)
3241  matrix.vmult(
3242  temp_vector1,
3243  solution,
3244  [&](const unsigned int start_range, const unsigned int end_range) {
3245  // zero 'temp_vector1' before running the vmult
3246  // operation
3247  if (end_range > start_range)
3248  std::memset(temp_vector1.begin() + start_range,
3249  0,
3250  sizeof(Number) * (end_range - start_range));
3251  },
3252  [&](const unsigned int start_range, const unsigned int end_range) {
3253  if (end_range > start_range)
3254  updater.apply_to_subrange(start_range, end_range);
3255  });
3256  else
3257  updater.apply_to_subrange(0U, solution.locally_owned_size());
3258 
3259  // swap vectors x^{n+1}->x^{n}, given the updates in the function above
3260  if (iteration_index == 0)
3261  {
3262  // nothing to do here because we can immediately write into the
3263  // solution vector without remembering any of the other vectors
3264  }
3265  else
3266  {
3267  solution.swap(temp_vector1);
3268  solution_old.swap(temp_vector1);
3269  }
3270  }
3271 
3272  template <typename MatrixType, typename PreconditionerType>
3273  inline void
3274  initialize_preconditioner(
3275  const MatrixType & matrix,
3276  std::shared_ptr<PreconditionerType> &preconditioner)
3277  {
3278  (void)matrix;
3279  (void)preconditioner;
3280  AssertThrow(preconditioner.get() != nullptr, ExcNotInitialized());
3281  }
3282 
3283  template <typename MatrixType, typename VectorType>
3284  inline void
3285  initialize_preconditioner(
3286  const MatrixType & matrix,
3287  std::shared_ptr<::DiagonalMatrix<VectorType>> &preconditioner)
3288  {
3289  if (preconditioner.get() == nullptr || preconditioner->m() != matrix.m())
3290  {
3291  if (preconditioner.get() == nullptr)
3292  preconditioner =
3293  std::make_shared<::DiagonalMatrix<VectorType>>();
3294 
3295  Assert(
3296  preconditioner->m() == 0,
3297  ExcMessage(
3298  "Preconditioner appears to be initialized but not sized correctly"));
3299 
3300  // This part only works in serial
3301  if (preconditioner->m() != matrix.m())
3302  {
3303  preconditioner->get_vector().reinit(matrix.m());
3304  for (typename VectorType::size_type i = 0; i < matrix.m(); ++i)
3305  preconditioner->get_vector()(i) = 1. / matrix.el(i, i);
3306  }
3307  }
3308  }
3309 
3310  template <typename VectorType>
3311  void
3312  set_initial_guess(VectorType &vector)
3313  {
3314  vector = 1. / std::sqrt(static_cast<double>(vector.size()));
3315  if (vector.locally_owned_elements().is_element(0))
3316  vector(0) = 0.;
3317  }
3318 
3319  template <typename Number>
3320  void
3321  set_initial_guess(::Vector<Number> &vector)
3322  {
3323  // Choose a high-frequency mode consisting of numbers between 0 and 1
3324  // that is cheap to compute (cheaper than random numbers) but avoids
3325  // obviously re-occurring numbers in multi-component systems by choosing
3326  // a period of 11
3327  for (unsigned int i = 0; i < vector.size(); ++i)
3328  vector(i) = i % 11;
3329 
3330  const Number mean_value = vector.mean_value();
3331  vector.add(-mean_value);
3332  }
3333 
3334  template <typename Number>
3335  void
3336  set_initial_guess(
3338  &vector)
3339  {
3340  // Choose a high-frequency mode consisting of numbers between 0 and 1
3341  // that is cheap to compute (cheaper than random numbers) but avoids
3342  // obviously re-occurring numbers in multi-component systems by choosing
3343  // a period of 11.
3344  // Make initial guess robust with respect to number of processors
3345  // by operating on the global index.
3346  types::global_dof_index first_local_range = 0;
3347  if (!vector.locally_owned_elements().is_empty())
3348  first_local_range = vector.locally_owned_elements().nth_index_in_set(0);
3349  for (unsigned int i = 0; i < vector.locally_owned_size(); ++i)
3350  vector.local_element(i) = (i + first_local_range) % 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 
3366 # ifdef DEAL_II_COMPILER_CUDA_AWARE
3367  template <typename Number>
3368  __global__ void
3369  set_initial_guess_kernel(const types::global_dof_index offset,
3370  const unsigned int locally_owned_size,
3371  Number * values)
3372 
3373  {
3374  const unsigned int index = threadIdx.x + blockDim.x * blockIdx.x;
3375  if (index < locally_owned_size)
3376  values[index] = (index + offset) % 11;
3377  }
3378 
3379  template <typename Number>
3380  void
3381  set_initial_guess(
3383  &vector)
3384  {
3385  // Choose a high-frequency mode consisting of numbers between 0 and 1
3386  // that is cheap to compute (cheaper than random numbers) but avoids
3387  // obviously re-occurring numbers in multi-component systems by choosing
3388  // a period of 11.
3389  // Make initial guess robust with respect to number of processors
3390  // by operating on the global index.
3391  types::global_dof_index first_local_range = 0;
3392  if (!vector.locally_owned_elements().is_empty())
3393  first_local_range = vector.locally_owned_elements().nth_index_in_set(0);
3394 
3395  const auto n_local_elements = vector.locally_owned_size();
3396  const int n_blocks =
3397  1 + (n_local_elements - 1) / CUDAWrappers::block_size;
3398  set_initial_guess_kernel<<<n_blocks, CUDAWrappers::block_size>>>(
3399  first_local_range, n_local_elements, vector.get_values());
3400  AssertCudaKernel();
3401 
3402  const Number mean_value = vector.mean_value();
3403  vector.add(-mean_value);
3404  }
3405 # endif // DEAL_II_COMPILER_CUDA_AWARE
3406 
3407  struct EigenvalueTracker
3408  {
3409  public:
3410  void
3411  slot(const std::vector<double> &eigenvalues)
3412  {
3413  values = eigenvalues;
3414  }
3415 
3416  std::vector<double> values;
3417  };
3418 
3419 
3420 
3421  template <typename MatrixType,
3422  typename VectorType,
3423  typename PreconditionerType>
3424  double
3425  power_iteration(const MatrixType & matrix,
3426  VectorType & eigenvector,
3427  const PreconditionerType &preconditioner,
3428  const unsigned int n_iterations)
3429  {
3430  double eigenvalue_estimate = 0.;
3431  eigenvector /= eigenvector.l2_norm();
3432  VectorType vector1, vector2;
3433  vector1.reinit(eigenvector, true);
3434  if (!std::is_same<PreconditionerType, PreconditionIdentity>::value)
3435  vector2.reinit(eigenvector, true);
3436  for (unsigned int i = 0; i < n_iterations; ++i)
3437  {
3438  if (!std::is_same<PreconditionerType, PreconditionIdentity>::value)
3439  {
3440  matrix.vmult(vector2, eigenvector);
3441  preconditioner.vmult(vector1, vector2);
3442  }
3443  else
3444  matrix.vmult(vector1, eigenvector);
3445 
3446  eigenvalue_estimate = eigenvector * vector1;
3447 
3448  vector1 /= vector1.l2_norm();
3449  eigenvector.swap(vector1);
3450  }
3451  return eigenvalue_estimate;
3452  }
3453  } // namespace PreconditionChebyshevImplementation
3454 } // namespace internal
3455 
3456 
3457 
3458 template <typename MatrixType, class VectorType, typename PreconditionerType>
3460  AdditionalData::AdditionalData(const unsigned int degree,
3461  const double smoothing_range,
3462  const unsigned int eig_cg_n_iterations,
3463  const double eig_cg_residual,
3464  const double max_eigenvalue,
3465  const EigenvalueAlgorithm eigenvalue_algorithm)
3466  : degree(degree)
3467  , smoothing_range(smoothing_range)
3468  , eig_cg_n_iterations(eig_cg_n_iterations)
3469  , eig_cg_residual(eig_cg_residual)
3470  , max_eigenvalue(max_eigenvalue)
3471  , eigenvalue_algorithm(eigenvalue_algorithm)
3472 {}
3473 
3474 
3475 
3476 template <typename MatrixType, class VectorType, typename PreconditionerType>
3477 inline typename PreconditionChebyshev<MatrixType,
3478  VectorType,
3479  PreconditionerType>::AdditionalData &
3481  AdditionalData::operator=(const AdditionalData &other_data)
3482 {
3483  degree = other_data.degree;
3484  smoothing_range = other_data.smoothing_range;
3485  eig_cg_n_iterations = other_data.eig_cg_n_iterations;
3486  eig_cg_residual = other_data.eig_cg_residual;
3487  max_eigenvalue = other_data.max_eigenvalue;
3488  preconditioner = other_data.preconditioner;
3489  eigenvalue_algorithm = other_data.eigenvalue_algorithm;
3490  constraints.copy_from(other_data.constraints);
3491 
3492  return *this;
3493 }
3494 
3495 
3496 
3497 template <typename MatrixType, typename VectorType, typename PreconditionerType>
3500  : theta(1.)
3501  , delta(1.)
3502  , eigenvalues_are_initialized(false)
3503 {
3504  static_assert(
3505  std::is_same<size_type, typename VectorType::size_type>::value,
3506  "PreconditionChebyshev and VectorType must have the same size_type.");
3507 }
3508 
3509 
3510 
3511 template <typename MatrixType, typename VectorType, typename PreconditionerType>
3512 inline void
3514  const MatrixType & matrix,
3515  const AdditionalData &additional_data)
3516 {
3517  matrix_ptr = &matrix;
3518  data = additional_data;
3519  Assert(data.degree > 0,
3520  ExcMessage("The degree of the Chebyshev method must be positive."));
3521  internal::PreconditionChebyshevImplementation::initialize_preconditioner(
3522  matrix, data.preconditioner);
3523  eigenvalues_are_initialized = false;
3524 }
3525 
3526 
3527 
3528 template <typename MatrixType, typename VectorType, typename PreconditionerType>
3529 inline void
3531 {
3532  eigenvalues_are_initialized = false;
3533  theta = delta = 1.0;
3534  matrix_ptr = nullptr;
3535  {
3536  VectorType empty_vector;
3537  solution_old.reinit(empty_vector);
3538  temp_vector1.reinit(empty_vector);
3539  temp_vector2.reinit(empty_vector);
3540  }
3541  data.preconditioner.reset();
3542 }
3543 
3544 
3545 
3546 template <typename MatrixType, typename VectorType, typename PreconditionerType>
3547 inline typename PreconditionChebyshev<MatrixType,
3548  VectorType,
3549  PreconditionerType>::EigenvalueInformation
3551  estimate_eigenvalues(const VectorType &src) const
3552 {
3553  Assert(eigenvalues_are_initialized == false, ExcInternalError());
3554  Assert(data.preconditioner.get() != nullptr, ExcNotInitialized());
3555 
3557  EigenvalueInformation info{};
3558 
3559  solution_old.reinit(src);
3560  temp_vector1.reinit(src, true);
3561 
3562  if (data.eig_cg_n_iterations > 0)
3563  {
3564  Assert(data.eig_cg_n_iterations > 2,
3565  ExcMessage(
3566  "Need to set at least two iterations to find eigenvalues."));
3567 
3568  internal::PreconditionChebyshevImplementation::EigenvalueTracker
3569  eigenvalue_tracker;
3570 
3571  // set an initial guess that contains some high-frequency parts (to the
3572  // extent possible without knowing the discretization and the numbering)
3573  // to trigger high eigenvalues according to the external function
3574  internal::PreconditionChebyshevImplementation::set_initial_guess(
3575  temp_vector1);
3576  data.constraints.set_zero(temp_vector1);
3577 
3578  if (data.eigenvalue_algorithm ==
3579  AdditionalData::EigenvalueAlgorithm::lanczos)
3580  {
3581  // set a very strict tolerance to force at least two iterations
3582  IterationNumberControl control(data.eig_cg_n_iterations,
3583  1e-10,
3584  false,
3585  false);
3586 
3587  SolverCG<VectorType> solver(control);
3588  solver.connect_eigenvalues_slot(
3589  [&eigenvalue_tracker](const std::vector<double> &eigenvalues) {
3590  eigenvalue_tracker.slot(eigenvalues);
3591  });
3592 
3593  solver.solve(*matrix_ptr,
3594  solution_old,
3595  temp_vector1,
3596  *data.preconditioner);
3597 
3598  info.cg_iterations = control.last_step();
3599  }
3600  else if (data.eigenvalue_algorithm ==
3601  AdditionalData::EigenvalueAlgorithm::power_iteration)
3602  {
3603  Assert(data.degree != numbers::invalid_unsigned_int,
3604  ExcMessage("Cannot estimate the minimal eigenvalue with the "
3605  "power iteration"));
3606 
3607  eigenvalue_tracker.values.push_back(
3608  internal::PreconditionChebyshevImplementation::power_iteration(
3609  *matrix_ptr,
3610  temp_vector1,
3611  *data.preconditioner,
3612  data.eig_cg_n_iterations));
3613  }
3614  else
3615  Assert(false, ExcNotImplemented());
3616 
3617  // read the eigenvalues from the attached eigenvalue tracker
3618  if (eigenvalue_tracker.values.empty())
3619  info.min_eigenvalue_estimate = info.max_eigenvalue_estimate = 1.;
3620  else
3621  {
3622  info.min_eigenvalue_estimate = eigenvalue_tracker.values.front();
3623 
3624  // include a safety factor since the CG method will in general not
3625  // be converged
3626  info.max_eigenvalue_estimate = 1.2 * eigenvalue_tracker.values.back();
3627  }
3628  }
3629  else
3630  {
3631  info.max_eigenvalue_estimate = data.max_eigenvalue;
3632  info.min_eigenvalue_estimate = data.max_eigenvalue / data.smoothing_range;
3633  }
3634 
3635  const double alpha = (data.smoothing_range > 1. ?
3636  info.max_eigenvalue_estimate / data.smoothing_range :
3637  std::min(0.9 * info.max_eigenvalue_estimate,
3638  info.min_eigenvalue_estimate));
3639 
3640  // in case the user set the degree to invalid unsigned int, we have to
3641  // determine the number of necessary iterations from the Chebyshev error
3642  // estimate, given the target tolerance specified by smoothing_range. This
3643  // estimate is based on the error formula given in section 5.1 of
3644  // R. S. Varga, Matrix iterative analysis, 2nd ed., Springer, 2009
3645  if (data.degree == numbers::invalid_unsigned_int)
3646  {
3647  const double actual_range = info.max_eigenvalue_estimate / alpha;
3648  const double sigma = (1. - std::sqrt(1. / actual_range)) /
3649  (1. + std::sqrt(1. / actual_range));
3650  const double eps = data.smoothing_range;
3651  const_cast<
3653  this)
3654  ->data.degree =
3655  1 + static_cast<unsigned int>(
3656  std::log(1. / eps + std::sqrt(1. / eps / eps - 1.)) /
3657  std::log(1. / sigma));
3658  }
3659 
3660  info.degree = data.degree;
3661 
3662  const_cast<
3664  ->delta = (info.max_eigenvalue_estimate - alpha) * 0.5;
3665  const_cast<
3667  ->theta = (info.max_eigenvalue_estimate + alpha) * 0.5;
3668 
3669  // We do not need the second temporary vector in case we have a
3670  // DiagonalMatrix as preconditioner and use deal.II's own vectors
3671  using NumberType = typename VectorType::value_type;
3672  if (std::is_same<PreconditionerType,
3673  ::DiagonalMatrix<VectorType>>::value == false ||
3674  (std::is_same<VectorType, ::Vector<NumberType>>::value == false &&
3675  ((std::is_same<VectorType,
3678  false) ||
3679  (std::is_same<VectorType,
3680  LinearAlgebra::distributed::
3681  Vector<NumberType, MemorySpace::CUDA>>::value ==
3682  false))))
3683  temp_vector2.reinit(src, true);
3684  else
3685  {
3686  VectorType empty_vector;
3687  temp_vector2.reinit(empty_vector);
3688  }
3689 
3690  const_cast<
3692  ->eigenvalues_are_initialized = true;
3693 
3694  return info;
3695 }
3696 
3697 
3698 
3699 template <typename MatrixType, typename VectorType, typename PreconditionerType>
3700 inline void
3702  VectorType & solution,
3703  const VectorType &rhs) const
3704 {
3705  std::lock_guard<std::mutex> lock(mutex);
3706  if (eigenvalues_are_initialized == false)
3707  estimate_eigenvalues(rhs);
3708 
3709  internal::PreconditionChebyshevImplementation::vmult_and_update(
3710  *matrix_ptr,
3711  *data.preconditioner,
3712  rhs,
3713  0,
3714  0.,
3715  1. / theta,
3716  solution,
3717  solution_old,
3718  temp_vector1,
3719  temp_vector2);
3720 
3721  // if delta is zero, we do not need to iterate because the updates will be
3722  // zero
3723  if (data.degree < 2 || std::abs(delta) < 1e-40)
3724  return;
3725 
3726  double rhok = delta / theta, sigma = theta / delta;
3727  for (unsigned int k = 0; k < data.degree - 1; ++k)
3728  {
3729  const double rhokp = 1. / (2. * sigma - rhok);
3730  const double factor1 = rhokp * rhok, factor2 = 2. * rhokp / delta;
3731  rhok = rhokp;
3732  internal::PreconditionChebyshevImplementation::vmult_and_update(
3733  *matrix_ptr,
3734  *data.preconditioner,
3735  rhs,
3736  k + 1,
3737  factor1,
3738  factor2,
3739  solution,
3740  solution_old,
3741  temp_vector1,
3742  temp_vector2);
3743  }
3744 }
3745 
3746 
3747 
3748 template <typename MatrixType, typename VectorType, typename PreconditionerType>
3749 inline void
3751  VectorType & solution,
3752  const VectorType &rhs) const
3753 {
3754  std::lock_guard<std::mutex> lock(mutex);
3755  if (eigenvalues_are_initialized == false)
3756  estimate_eigenvalues(rhs);
3757 
3758  internal::PreconditionChebyshevImplementation::vector_updates(
3759  rhs,
3760  *data.preconditioner,
3761  0,
3762  0.,
3763  1. / theta,
3764  solution_old,
3765  temp_vector1,
3766  temp_vector2,
3767  solution);
3768 
3769  if (data.degree < 2 || std::abs(delta) < 1e-40)
3770  return;
3771 
3772  double rhok = delta / theta, sigma = theta / delta;
3773  for (unsigned int k = 0; k < data.degree - 1; ++k)
3774  {
3775  const double rhokp = 1. / (2. * sigma - rhok);
3776  const double factor1 = rhokp * rhok, factor2 = 2. * rhokp / delta;
3777  rhok = rhokp;
3778  matrix_ptr->Tvmult(temp_vector1, solution);
3779  internal::PreconditionChebyshevImplementation::vector_updates(
3780  rhs,
3781  *data.preconditioner,
3782  k + 1,
3783  factor1,
3784  factor2,
3785  solution_old,
3786  temp_vector1,
3787  temp_vector2,
3788  solution);
3789  }
3790 }
3791 
3792 
3793 
3794 template <typename MatrixType, typename VectorType, typename PreconditionerType>
3795 inline void
3797  VectorType & solution,
3798  const VectorType &rhs) const
3799 {
3800  std::lock_guard<std::mutex> lock(mutex);
3801  if (eigenvalues_are_initialized == false)
3802  estimate_eigenvalues(rhs);
3803 
3804  internal::PreconditionChebyshevImplementation::vmult_and_update(
3805  *matrix_ptr,
3806  *data.preconditioner,
3807  rhs,
3808  1,
3809  0.,
3810  1. / theta,
3811  solution,
3812  solution_old,
3813  temp_vector1,
3814  temp_vector2);
3815 
3816  if (data.degree < 2 || std::abs(delta) < 1e-40)
3817  return;
3818 
3819  double rhok = delta / theta, sigma = theta / delta;
3820  for (unsigned int k = 0; k < data.degree - 1; ++k)
3821  {
3822  const double rhokp = 1. / (2. * sigma - rhok);
3823  const double factor1 = rhokp * rhok, factor2 = 2. * rhokp / delta;
3824  rhok = rhokp;
3825  internal::PreconditionChebyshevImplementation::vmult_and_update(
3826  *matrix_ptr,
3827  *data.preconditioner,
3828  rhs,
3829  k + 2,
3830  factor1,
3831  factor2,
3832  solution,
3833  solution_old,
3834  temp_vector1,
3835  temp_vector2);
3836  }
3837 }
3838 
3839 
3840 
3841 template <typename MatrixType, typename VectorType, typename PreconditionerType>
3842 inline void
3844  VectorType & solution,
3845  const VectorType &rhs) const
3846 {
3847  std::lock_guard<std::mutex> lock(mutex);
3848  if (eigenvalues_are_initialized == false)
3849  estimate_eigenvalues(rhs);
3850 
3851  matrix_ptr->Tvmult(temp_vector1, solution);
3852  internal::PreconditionChebyshevImplementation::vector_updates(
3853  rhs,
3854  *data.preconditioner,
3855  1,
3856  0.,
3857  1. / theta,
3858  solution_old,
3859  temp_vector1,
3860  temp_vector2,
3861  solution);
3862 
3863  if (data.degree < 2 || std::abs(delta) < 1e-40)
3864  return;
3865 
3866  double rhok = delta / theta, sigma = theta / delta;
3867  for (unsigned int k = 0; k < data.degree - 1; ++k)
3868  {
3869  const double rhokp = 1. / (2. * sigma - rhok);
3870  const double factor1 = rhokp * rhok, factor2 = 2. * rhokp / delta;
3871  rhok = rhokp;
3872  matrix_ptr->Tvmult(temp_vector1, solution);
3873  internal::PreconditionChebyshevImplementation::vector_updates(
3874  rhs,
3875  *data.preconditioner,
3876  k + 2,
3877  factor1,
3878  factor2,
3879  solution_old,
3880  temp_vector1,
3881  temp_vector2,
3882  solution);
3883  }
3884 }
3885 
3886 
3887 
3888 template <typename MatrixType, typename VectorType, typename PreconditionerType>
3889 inline typename PreconditionChebyshev<MatrixType,
3890  VectorType,
3891  PreconditionerType>::size_type
3893 {
3894  Assert(matrix_ptr != nullptr, ExcNotInitialized());
3895  return matrix_ptr->m();
3896 }
3897 
3898 
3899 
3900 template <typename MatrixType, typename VectorType, typename PreconditionerType>
3901 inline typename PreconditionChebyshev<MatrixType,
3902  VectorType,
3903  PreconditionerType>::size_type
3905 {
3906  Assert(matrix_ptr != nullptr, ExcNotInitialized());
3907  return matrix_ptr->n();
3908 }
3909 
3910 #endif // DOXYGEN
3911 
3913 
3914 #endif
unsigned int n_blocks() const
BlockType & block(const unsigned int i)
VectorType & get_vector()
bool is_empty() const
Definition: index_set.h:1766
size_type nth_index_in_set(const size_type local_index) const
Definition: index_set.h:1822
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:142
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:458
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:459
static ::ExceptionBase & ExcInternalError()
#define AssertCudaKernel()
Definition: exceptions.h:1899
static ::ExceptionBase & ExcNotInitialized()
#define Assert(cond, exc)
Definition: exceptions.h:1501
static ::ExceptionBase & ExcNotImplemented()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1611
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:206
unsigned int global_dof_index
Definition: types.h:81
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
std::shared_ptr< PreconditionerType > preconditioner
AffineConstraints< double > constraints
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)
EigenvalueAlgorithm eigenvalue_algorithm
AdditionalData & operator=(const AdditionalData &other_data)