Reference documentation for deal.II version GIT 6da2e5d553 2022-07-01 18:55: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 - 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/parallel.h>
27 #include <deal.II/base/utilities.h>
28 
33 #include <deal.II/lac/solver_cg.h>
35 
37 
38 // forward declarations
39 #ifndef DOXYGEN
40 template <typename number>
41 class Vector;
42 template <typename number>
43 class SparseMatrix;
44 namespace LinearAlgebra
45 {
46  namespace distributed
47  {
48  template <typename, typename>
49  class Vector;
50  template <typename>
51  class BlockVector;
52  } // namespace distributed
53 } // namespace LinearAlgebra
54 #endif
55 
56 
82 {
83 public:
88 
94  {
98  AdditionalData() = default;
99  };
100 
105 
110  template <typename MatrixType>
111  void
112  initialize(const MatrixType & matrix,
113  const AdditionalData &additional_data = AdditionalData());
114 
118  template <class VectorType>
119  void
120  vmult(VectorType &, const VectorType &) const;
121 
126  template <class VectorType>
127  void
128  Tvmult(VectorType &, const VectorType &) const;
129 
133  template <class VectorType>
134  void
135  vmult_add(VectorType &, const VectorType &) const;
136 
141  template <class VectorType>
142  void
143  Tvmult_add(VectorType &, const VectorType &) const;
144 
149  void
150  clear();
151 
159  size_type
160  m() const;
161 
169  size_type
170  n() const;
171 
172 private:
177 
182 };
183 
184 
185 
197 {
198 public:
203 
208  {
209  public:
214  AdditionalData(const double relaxation = 1.);
215 
219  double relaxation;
220  };
221 
227 
231  void
232  initialize(const AdditionalData &parameters);
233 
239  template <typename MatrixType>
240  void
241  initialize(const MatrixType &matrix, const AdditionalData &parameters);
242 
246  template <class VectorType>
247  void
248  vmult(VectorType &, const VectorType &) const;
249 
254  template <class VectorType>
255  void
256  Tvmult(VectorType &, const VectorType &) const;
260  template <class VectorType>
261  void
262  vmult_add(VectorType &, const VectorType &) const;
263 
268  template <class VectorType>
269  void
270  Tvmult_add(VectorType &, const VectorType &) const;
271 
276  void
278  {}
279 
287  size_type
288  m() const;
289 
297  size_type
298  n() const;
299 
300 private:
304  double relaxation;
305 
310 
315 };
316 
317 
318 
358 template <typename MatrixType = SparseMatrix<double>,
359  class VectorType = Vector<double>>
361 {
362 public:
366  using function_ptr = void (MatrixType::*)(VectorType &,
367  const VectorType &) const;
368 
374  PreconditionUseMatrix(const MatrixType &M, const function_ptr method);
375 
380  void
381  vmult(VectorType &dst, const VectorType &src) const;
382 
383 private:
387  const MatrixType &matrix;
388 
393 };
394 
395 
396 
402 template <typename MatrixType = SparseMatrix<double>,
403  typename PreconditionerType = IdentityMatrix>
405 {
406 public:
411 
416  {
417  public:
421  AdditionalData(const double relaxation = 1.,
422  const unsigned int n_iterations = 1);
423 
427  double relaxation;
428 
432  unsigned int n_iterations;
433 
434 
435  /*
436  * Preconditioner.
437  */
438  std::shared_ptr<PreconditionerType> preconditioner;
439  };
440 
446  void
447  initialize(const MatrixType & A,
448  const AdditionalData &parameters = AdditionalData());
449 
453  void
454  clear();
455 
460  size_type
461  m() const;
462 
467  size_type
468  n() const;
469 
473  template <class VectorType>
474  void
475  vmult(VectorType &, const VectorType &) const;
476 
481  template <class VectorType>
482  void
483  Tvmult(VectorType &, const VectorType &) const;
484 
488  template <class VectorType>
489  void
490  step(VectorType &x, const VectorType &rhs) const;
491 
495  template <class VectorType>
496  void
497  Tstep(VectorType &x, const VectorType &rhs) const;
498 
499 protected:
504 
508  double relaxation;
509 
513  unsigned int n_iterations;
514 
515  /*
516  * Preconditioner.
517  */
518  std::shared_ptr<PreconditionerType> preconditioner;
519 };
520 
521 
522 
523 #ifndef DOXYGEN
524 
525 namespace internal
526 {
527  namespace PreconditionRelaxation
528  {
529  template <typename T, typename VectorType>
530  using Tvmult_t = decltype(
531  std::declval<T const>().Tvmult(std::declval<VectorType &>(),
532  std::declval<const VectorType &>()));
533 
534  template <typename T, typename VectorType>
535  constexpr bool has_Tvmult = is_supported_operation<Tvmult_t, T, VectorType>;
536 
537  template <typename T, typename VectorType>
538  using step_t = decltype(
539  std::declval<T const>().step(std::declval<VectorType &>(),
540  std::declval<const VectorType &>()));
541 
542  template <typename T, typename VectorType>
543  constexpr bool has_step = is_supported_operation<step_t, T, VectorType>;
544 
545  template <typename T, typename VectorType>
546  using step_omega_t =
547  decltype(std::declval<T const>().step(std::declval<VectorType &>(),
548  std::declval<const VectorType &>(),
549  std::declval<const double>()));
550 
551  template <typename T, typename VectorType>
552  constexpr bool has_step_omega =
553  is_supported_operation<step_omega_t, T, VectorType>;
554 
555  template <typename T, typename VectorType>
556  using Tstep_t = decltype(
557  std::declval<T const>().Tstep(std::declval<VectorType &>(),
558  std::declval<const VectorType &>()));
559 
560  template <typename T, typename VectorType>
561  constexpr bool has_Tstep = is_supported_operation<Tstep_t, T, VectorType>;
562 
563  template <typename T, typename VectorType>
564  using Tstep_omega_t =
565  decltype(std::declval<T const>().Tstep(std::declval<VectorType &>(),
566  std::declval<const VectorType &>(),
567  std::declval<const double>()));
568 
569  template <typename T, typename VectorType>
570  constexpr bool has_Tstep_omega =
571  is_supported_operation<Tstep_omega_t, T, VectorType>;
572 
573  template <typename T, typename VectorType>
574  using jacobi_step_t = decltype(
575  std::declval<T const>().Jacobi_step(std::declval<VectorType &>(),
576  std::declval<const VectorType &>(),
577  std::declval<const double>()));
578 
579  template <typename T, typename VectorType>
580  constexpr bool has_jacobi_step =
581  is_supported_operation<jacobi_step_t, T, VectorType>;
582 
583  template <typename T, typename VectorType>
584  using SOR_step_t = decltype(
585  std::declval<T const>().SOR_step(std::declval<VectorType &>(),
586  std::declval<const VectorType &>(),
587  std::declval<const double>()));
588 
589  template <typename T, typename VectorType>
590  constexpr bool has_SOR_step =
591  is_supported_operation<SOR_step_t, T, VectorType>;
592 
593  template <typename T, typename VectorType>
594  using SSOR_step_t = decltype(
595  std::declval<T const>().SSOR_step(std::declval<VectorType &>(),
596  std::declval<const VectorType &>(),
597  std::declval<const double>()));
598 
599  template <typename T, typename VectorType>
600  constexpr bool has_SSOR_step =
601  is_supported_operation<SSOR_step_t, T, VectorType>;
602 
603  template <typename MatrixType>
604  class PreconditionJacobiImpl
605  {
606  public:
607  PreconditionJacobiImpl(const MatrixType &A, const double relaxation)
608  : A(&A)
609  , relaxation(relaxation)
610  {}
611 
612  template <typename VectorType>
613  void
614  vmult(VectorType &dst, const VectorType &src) const
615  {
616  this->A->precondition_Jacobi(dst, src, this->relaxation);
617  }
618 
619  template <typename VectorType>
620  void
621  Tvmult(VectorType &dst, const VectorType &src) const
622  {
623  // call vmult, since preconditioner is symmetrical
624  this->vmult(dst, src);
625  }
626 
627  template <typename VectorType,
628  typename std::enable_if<has_jacobi_step<MatrixType, VectorType>,
629  MatrixType>::type * = nullptr>
630  void
631  step(VectorType &dst, const VectorType &src) const
632  {
633  this->A->Jacobi_step(dst, src, this->relaxation);
634  }
635 
636  template <
637  typename VectorType,
638  typename std::enable_if<!has_jacobi_step<MatrixType, VectorType>,
639  MatrixType>::type * = nullptr>
640  void
641  step(VectorType &, const VectorType &) const
642  {
643  AssertThrow(false,
644  ExcMessage(
645  "Matrix A does not provide a Jacobi_step() function!"));
646  }
647 
648  template <typename VectorType>
649  void
650  Tstep(VectorType &dst, const VectorType &src) const
651  {
652  // call step, since preconditioner is symmetrical
653  this->step(dst, src);
654  }
655 
656  private:
658  const double relaxation;
659  };
660 
661  template <typename MatrixType>
662  class PreconditionSORImpl
663  {
664  public:
665  PreconditionSORImpl(const MatrixType &A, const double relaxation)
666  : A(&A)
667  , relaxation(relaxation)
668  {}
669 
670  template <typename VectorType>
671  void
672  vmult(VectorType &dst, const VectorType &src) const
673  {
674  this->A->precondition_SOR(dst, src, this->relaxation);
675  }
676 
677  template <typename VectorType>
678  void
679  Tvmult(VectorType &dst, const VectorType &src) const
680  {
681  this->A->precondition_TSOR(dst, src, this->relaxation);
682  }
683 
684  template <typename VectorType,
685  typename std::enable_if<has_SOR_step<MatrixType, VectorType>,
686  MatrixType>::type * = nullptr>
687  void
688  step(VectorType &dst, const VectorType &src) const
689  {
690  this->A->SOR_step(dst, src, this->relaxation);
691  }
692 
693  template <typename VectorType,
694  typename std::enable_if<!has_SOR_step<MatrixType, VectorType>,
695  MatrixType>::type * = nullptr>
696  void
697  step(VectorType &, const VectorType &) const
698  {
699  AssertThrow(false,
700  ExcMessage(
701  "Matrix A does not provide a SOR_step() function!"));
702  }
703 
704  template <typename VectorType,
705  typename std::enable_if<has_SOR_step<MatrixType, VectorType>,
706  MatrixType>::type * = nullptr>
707  void
708  Tstep(VectorType &dst, const VectorType &src) const
709  {
710  this->A->TSOR_step(dst, src, this->relaxation);
711  }
712 
713  template <typename VectorType,
714  typename std::enable_if<!has_SOR_step<MatrixType, VectorType>,
715  MatrixType>::type * = nullptr>
716  void
717  Tstep(VectorType &, const VectorType &) const
718  {
719  AssertThrow(false,
720  ExcMessage(
721  "Matrix A does not provide a TSOR_step() function!"));
722  }
723 
724  private:
726  const double relaxation;
727  };
728 
729  template <typename MatrixType>
730  class PreconditionSSORImpl
731  {
732  public:
733  using size_type = typename MatrixType::size_type;
734 
735  PreconditionSSORImpl(const MatrixType &A, const double relaxation)
736  : A(&A)
737  , relaxation(relaxation)
738  {
739  // in case we have a SparseMatrix class, we can extract information
740  // about the diagonal.
742  dynamic_cast<const SparseMatrix<typename MatrixType::value_type> *>(
743  &*this->A);
744 
745  // calculate the positions first after the diagonal.
746  if (mat != nullptr)
747  {
748  const size_type n = this->A->n();
749  pos_right_of_diagonal.resize(n, static_cast<std::size_t>(-1));
750  for (size_type row = 0; row < n; ++row)
751  {
752  // find the first element in this line which is on the right of
753  // the diagonal. we need to precondition with the elements on
754  // the left only. note: the first entry in each line denotes the
755  // diagonal element, which we need not check.
756  typename SparseMatrix<
757  typename MatrixType::value_type>::const_iterator it =
758  mat->begin(row) + 1;
759  for (; it < mat->end(row); ++it)
760  if (it->column() > row)
761  break;
762  pos_right_of_diagonal[row] = it - mat->begin();
763  }
764  }
765  }
766 
767  template <typename VectorType>
768  void
769  vmult(VectorType &dst, const VectorType &src) const
770  {
771  this->A->precondition_SSOR(dst,
772  src,
773  this->relaxation,
774  pos_right_of_diagonal);
775  }
776 
777  template <typename VectorType>
778  void
779  Tvmult(VectorType &dst, const VectorType &src) const
780  {
781  this->A->precondition_SSOR(dst,
782  src,
783  this->relaxation,
784  pos_right_of_diagonal);
785  }
786 
787  template <typename VectorType,
788  typename std::enable_if<has_SSOR_step<MatrixType, VectorType>,
789  MatrixType>::type * = nullptr>
790  void
791  step(VectorType &dst, const VectorType &src) const
792  {
793  this->A->SSOR_step(dst, src, this->relaxation);
794  }
795 
796  template <typename VectorType,
797  typename std::enable_if<!has_SSOR_step<MatrixType, VectorType>,
798  MatrixType>::type * = nullptr>
799  void
800  step(VectorType &, const VectorType &) const
801  {
802  AssertThrow(false,
803  ExcMessage(
804  "Matrix A does not provide a SSOR_step() function!"));
805  }
806 
807  template <typename VectorType>
808  void
809  Tstep(VectorType &dst, const VectorType &src) const
810  {
811  // call step, since preconditioner is symmetrical
812  this->step(dst, src);
813  }
814 
815  private:
817  const double relaxation;
818 
823  std::vector<std::size_t> pos_right_of_diagonal;
824  };
825 
826  template <typename MatrixType>
827  class PreconditionPSORImpl
828  {
829  public:
830  using size_type = typename MatrixType::size_type;
831 
832  PreconditionPSORImpl(const MatrixType & A,
833  const double relaxation,
834  const std::vector<size_type> &permutation,
835  const std::vector<size_type> &inverse_permutation)
836  : A(&A)
837  , relaxation(relaxation)
838  , permutation(permutation)
839  , inverse_permutation(inverse_permutation)
840  {}
841 
842  template <typename VectorType>
843  void
844  vmult(VectorType &dst, const VectorType &src) const
845  {
846  dst = src;
847  this->A->PSOR(dst, permutation, inverse_permutation, this->relaxation);
848  }
849 
850  template <typename VectorType>
851  void
852  Tvmult(VectorType &dst, const VectorType &src) const
853  {
854  dst = src;
855  this->A->TPSOR(dst, permutation, inverse_permutation, this->relaxation);
856  }
857 
858  private:
860  const double relaxation;
861 
862  const std::vector<size_type> &permutation;
863  const std::vector<size_type> &inverse_permutation;
864  };
865 
866  template <
867  typename MatrixType,
868  typename PreconditionerType,
869  typename VectorType,
870  typename std::enable_if<has_step_omega<PreconditionerType, VectorType>,
871  PreconditionerType>::type * = nullptr>
872  void
873  step(const MatrixType &,
874  const PreconditionerType &preconditioner,
875  VectorType & dst,
876  const VectorType & src,
877  const double relaxation,
878  VectorType &,
879  VectorType &)
880  {
881  preconditioner.step(dst, src, relaxation);
882  }
883 
884  template <
885  typename MatrixType,
886  typename PreconditionerType,
887  typename VectorType,
888  typename std::enable_if<!has_step_omega<PreconditionerType, VectorType> &&
889  has_step<PreconditionerType, VectorType>,
890  PreconditionerType>::type * = nullptr>
891  void
892  step(const MatrixType &,
893  const PreconditionerType &preconditioner,
894  VectorType & dst,
895  const VectorType & src,
896  const double relaxation,
897  VectorType &,
898  VectorType &)
899  {
900  Assert(relaxation == 1.0, ExcInternalError());
901 
902  (void)relaxation;
903 
904  preconditioner.step(dst, src);
905  }
906 
907  template <
908  typename MatrixType,
909  typename PreconditionerType,
910  typename VectorType,
911  typename std::enable_if<!has_step_omega<PreconditionerType, VectorType> &&
912  !has_step<PreconditionerType, VectorType>,
913  PreconditionerType>::type * = nullptr>
914  void
915  step(const MatrixType & A,
916  const PreconditionerType &preconditioner,
917  VectorType & dst,
918  const VectorType & src,
919  const double relaxation,
920  VectorType & residual,
921  VectorType & tmp)
922  {
923  residual.reinit(dst, true);
924  tmp.reinit(dst, true);
925 
926  A.vmult(residual, dst);
927  residual.sadd(-1.0, 1.0, src);
928 
929  preconditioner.vmult(tmp, residual);
930  dst.add(relaxation, tmp);
931  }
932 
933  template <
934  typename MatrixType,
935  typename PreconditionerType,
936  typename VectorType,
937  typename std::enable_if<has_Tstep_omega<PreconditionerType, VectorType>,
938  PreconditionerType>::type * = nullptr>
939  void
940  Tstep(const MatrixType &,
941  const PreconditionerType &preconditioner,
942  VectorType & dst,
943  const VectorType & src,
944  const double relaxation,
945  VectorType &,
946  VectorType &)
947  {
948  preconditioner.Tstep(dst, src, relaxation);
949  }
950 
951  template <typename MatrixType,
952  typename PreconditionerType,
953  typename VectorType,
954  typename std::enable_if<
955  !has_Tstep_omega<PreconditionerType, VectorType> &&
956  has_Tstep<PreconditionerType, VectorType>,
957  PreconditionerType>::type * = nullptr>
958  void
959  Tstep(const MatrixType &,
960  const PreconditionerType &preconditioner,
961  VectorType & dst,
962  const VectorType & src,
963  const double relaxation,
964  VectorType &,
965  VectorType &)
966  {
967  Assert(relaxation == 1.0, ExcInternalError());
968 
969  (void)relaxation;
970 
971  preconditioner.Tstep(dst, src);
972  }
973 
974  template <typename MatrixType,
975  typename VectorType,
976  typename std::enable_if<has_Tvmult<MatrixType, VectorType>,
977  MatrixType>::type * = nullptr>
978  void
979  Tvmult(const MatrixType &A, VectorType &dst, const VectorType &src)
980  {
981  A.Tvmult(dst, src);
982  }
983 
984  template <typename MatrixType,
985  typename VectorType,
986  typename std::enable_if<!has_Tvmult<MatrixType, VectorType>,
987  MatrixType>::type * = nullptr>
988  void
989  Tvmult(const MatrixType &, VectorType &, const VectorType &)
990  {
991  AssertThrow(false,
992  ExcMessage("Matrix A does not provide a Tvmult() function!"));
993  }
994 
995  template <typename MatrixType,
996  typename PreconditionerType,
997  typename VectorType,
998  typename std::enable_if<
999  !has_Tstep_omega<PreconditionerType, VectorType> &&
1000  !has_Tstep<PreconditionerType, VectorType>,
1001  PreconditionerType>::type * = nullptr>
1002  void
1003  Tstep(const MatrixType & A,
1004  const PreconditionerType &preconditioner,
1005  VectorType & dst,
1006  const VectorType & src,
1007  const double relaxation,
1008  VectorType & residual,
1009  VectorType & tmp)
1010  {
1011  residual.reinit(dst, true);
1012  tmp.reinit(dst, true);
1013 
1014  Tvmult(A, residual, dst);
1015  residual.sadd(-1.0, 1.0, src);
1016 
1017  Tvmult(preconditioner, tmp, residual);
1018  dst.add(relaxation, tmp);
1019  }
1020 
1021  template <typename MatrixType,
1022  typename PreconditionerType,
1023  typename VectorType>
1024  void
1025  step_operations(const MatrixType & A,
1026  const PreconditionerType &preconditioner,
1027  VectorType & dst,
1028  const VectorType & src,
1029  const double relaxation,
1030  VectorType & tmp1,
1031  VectorType & tmp2,
1032  const unsigned int i,
1033  const bool transposed)
1034  {
1035  if (i == 0)
1036  {
1037  if (transposed)
1038  Tvmult(preconditioner, dst, src);
1039  else
1040  preconditioner.vmult(dst, src);
1041 
1042  if (relaxation != 1.0)
1043  dst *= relaxation;
1044  }
1045  else
1046  {
1047  if (transposed)
1048  Tstep(A, preconditioner, dst, src, relaxation, tmp1, tmp2);
1049  else
1050  step(A, preconditioner, dst, src, relaxation, tmp1, tmp2);
1051  }
1052  }
1053 
1054  template <typename MatrixType,
1055  typename VectorType,
1056  typename std::enable_if<!IsBlockVector<VectorType>::value,
1057  VectorType>::type * = nullptr>
1058  void
1059  step_operations(const MatrixType & A,
1060  const DiagonalMatrix<VectorType> &preconditioner,
1061  VectorType & dst,
1062  const VectorType & src,
1063  const double relaxation,
1064  VectorType & residual,
1065  VectorType &,
1066  const unsigned int i,
1067  const bool transposed)
1068  {
1069  if (i == 0)
1070  {
1071  const auto dst_ptr = dst.begin();
1072  const auto src_ptr = src.begin();
1073  const auto diag_ptr = preconditioner.get_vector().begin();
1074 
1075  for (unsigned int i = 0; i < dst.locally_owned_size(); ++i)
1076  dst_ptr[i] = relaxation * src_ptr[i] * diag_ptr[i];
1077  }
1078  else
1079  {
1080  residual.reinit(src, true);
1081 
1082  const auto dst_ptr = dst.begin();
1083  const auto src_ptr = src.begin();
1084  const auto residual_ptr = residual.begin();
1085  const auto diag_ptr = preconditioner.get_vector().begin();
1086 
1087  if (transposed)
1088  A.Tvmult(residual, dst);
1089  else
1090  A.vmult(residual, dst);
1091 
1092  for (unsigned int i = 0; i < dst.locally_owned_size(); ++i)
1093  dst_ptr[i] +=
1094  relaxation * (src_ptr[i] - residual_ptr[i]) * diag_ptr[i];
1095  }
1096  }
1097 
1098  } // namespace PreconditionRelaxation
1099 } // namespace internal
1100 
1101 #endif
1102 
1103 
1104 
1131 template <typename MatrixType = SparseMatrix<double>>
1133  : public PreconditionRelaxation<
1134  MatrixType,
1135  internal::PreconditionRelaxation::PreconditionJacobiImpl<MatrixType>>
1136 {
1138  internal::PreconditionRelaxation::PreconditionJacobiImpl<MatrixType>;
1140 
1141 public:
1146 
1150  void
1151  initialize(const MatrixType & A,
1152  const AdditionalData &parameters = AdditionalData());
1153 };
1154 
1155 
1201 template <typename MatrixType = SparseMatrix<double>>
1203  : public PreconditionRelaxation<
1204  MatrixType,
1205  internal::PreconditionRelaxation::PreconditionSORImpl<MatrixType>>
1206 {
1208  internal::PreconditionRelaxation::PreconditionSORImpl<MatrixType>;
1210 
1211 public:
1216 
1220  void
1221  initialize(const MatrixType & A,
1222  const AdditionalData &parameters = AdditionalData());
1223 };
1224 
1225 
1226 
1253 template <typename MatrixType = SparseMatrix<double>>
1255  : public PreconditionRelaxation<
1256  MatrixType,
1257  internal::PreconditionRelaxation::PreconditionSSORImpl<MatrixType>>
1258 {
1260  internal::PreconditionRelaxation::PreconditionSSORImpl<MatrixType>;
1262 
1263 public:
1268 
1274  void
1275  initialize(const MatrixType & A,
1276  const AdditionalData &parameters = AdditionalData());
1277 };
1278 
1279 
1309 template <typename MatrixType = SparseMatrix<double>>
1311  : public PreconditionRelaxation<
1312  MatrixType,
1313  internal::PreconditionRelaxation::PreconditionPSORImpl<MatrixType>>
1314 {
1316  internal::PreconditionRelaxation::PreconditionPSORImpl<MatrixType>;
1318 
1319 public:
1324 
1329  {
1330  public:
1341  AdditionalData(const std::vector<size_type> &permutation,
1342  const std::vector<size_type> &inverse_permutation,
1343  const typename BaseClass::AdditionalData &parameters =
1344  typename BaseClass::AdditionalData());
1345 
1349  const std::vector<size_type> &permutation;
1353  const std::vector<size_type> &inverse_permutation;
1358  };
1359 
1371  void
1372  initialize(const MatrixType & A,
1373  const std::vector<size_type> & permutation,
1374  const std::vector<size_type> & inverse_permutation,
1375  const typename BaseClass::AdditionalData &parameters =
1376  typename BaseClass::AdditionalData());
1377 
1388  void
1389  initialize(const MatrixType &A, const AdditionalData &additional_data);
1390 };
1391 
1392 
1393 
1591 template <typename MatrixType = SparseMatrix<double>,
1592  typename VectorType = Vector<double>,
1593  typename PreconditionerType = DiagonalMatrix<VectorType>>
1595 {
1596 public:
1601 
1607  {
1613  {
1619  lanczos,
1629  };
1630 
1634  AdditionalData(const unsigned int degree = 1,
1635  const double smoothing_range = 0.,
1636  const unsigned int eig_cg_n_iterations = 8,
1637  const double eig_cg_residual = 1e-2,
1638  const double max_eigenvalue = 1,
1641 
1645  AdditionalData &
1646  operator=(const AdditionalData &other_data);
1647 
1660  unsigned int degree;
1661 
1674 
1681  unsigned int eig_cg_n_iterations;
1682 
1688 
1695 
1701 
1705  std::shared_ptr<PreconditionerType> preconditioner;
1706 
1711  };
1712 
1713 
1718 
1730  void
1731  initialize(const MatrixType & matrix,
1732  const AdditionalData &additional_data = AdditionalData());
1733 
1738  void
1739  vmult(VectorType &dst, const VectorType &src) const;
1740 
1745  void
1746  Tvmult(VectorType &dst, const VectorType &src) const;
1747 
1751  void
1752  step(VectorType &dst, const VectorType &src) const;
1753 
1757  void
1758  Tstep(VectorType &dst, const VectorType &src) const;
1759 
1763  void
1765 
1770  size_type
1771  m() const;
1772 
1777  size_type
1778  n() const;
1779 
1785  {
1797  unsigned int cg_iterations;
1802  unsigned int degree;
1807  : min_eigenvalue_estimate{std::numeric_limits<double>::max()}
1808  , max_eigenvalue_estimate{std::numeric_limits<double>::lowest()}
1809  , cg_iterations{0}
1810  , degree{0}
1811  {}
1812  };
1813 
1826  EigenvalueInformation
1827  estimate_eigenvalues(const VectorType &src) const;
1828 
1829 private:
1833  SmartPointer<
1834  const MatrixType,
1837 
1841  mutable VectorType solution_old;
1842 
1846  mutable VectorType temp_vector1;
1847 
1851  mutable VectorType temp_vector2;
1852 
1858 
1862  double theta;
1863 
1868  double delta;
1869 
1875 
1881 };
1882 
1883 
1884 
1886 /* ---------------------------------- Inline functions ------------------- */
1887 
1888 #ifndef DOXYGEN
1889 
1891  : n_rows(0)
1892  , n_columns(0)
1893 {}
1894 
1895 template <typename MatrixType>
1896 inline void
1897 PreconditionIdentity::initialize(const MatrixType &matrix,
1899 {
1900  n_rows = matrix.m();
1901  n_columns = matrix.n();
1902 }
1903 
1904 
1905 template <class VectorType>
1906 inline void
1907 PreconditionIdentity::vmult(VectorType &dst, const VectorType &src) const
1908 {
1909  dst = src;
1910 }
1911 
1912 
1913 
1914 template <class VectorType>
1915 inline void
1916 PreconditionIdentity::Tvmult(VectorType &dst, const VectorType &src) const
1917 {
1918  dst = src;
1919 }
1920 
1921 template <class VectorType>
1922 inline void
1923 PreconditionIdentity::vmult_add(VectorType &dst, const VectorType &src) const
1924 {
1925  dst += src;
1926 }
1927 
1928 
1929 
1930 template <class VectorType>
1931 inline void
1932 PreconditionIdentity::Tvmult_add(VectorType &dst, const VectorType &src) const
1933 {
1934  dst += src;
1935 }
1936 
1937 
1938 
1939 inline void
1941 {}
1942 
1943 
1944 
1947 {
1948  Assert(n_rows != 0, ExcNotInitialized());
1949  return n_rows;
1950 }
1951 
1954 {
1956  return n_columns;
1957 }
1958 
1959 //---------------------------------------------------------------------------
1960 
1962  const double relaxation)
1963  : relaxation(relaxation)
1964 {}
1965 
1966 
1968  : relaxation(0)
1969  , n_rows(0)
1970  , n_columns(0)
1971 {
1972  AdditionalData add_data;
1973  relaxation = add_data.relaxation;
1974 }
1975 
1976 
1977 
1978 inline void
1980  const PreconditionRichardson::AdditionalData &parameters)
1981 {
1982  relaxation = parameters.relaxation;
1983 }
1984 
1985 
1986 
1987 template <typename MatrixType>
1988 inline void
1990  const MatrixType & matrix,
1991  const PreconditionRichardson::AdditionalData &parameters)
1992 {
1993  relaxation = parameters.relaxation;
1994  n_rows = matrix.m();
1995  n_columns = matrix.n();
1996 }
1997 
1998 
1999 
2000 template <class VectorType>
2001 inline void
2002 PreconditionRichardson::vmult(VectorType &dst, const VectorType &src) const
2003 {
2004  static_assert(
2005  std::is_same<size_type, typename VectorType::size_type>::value,
2006  "PreconditionRichardson and VectorType must have the same size_type.");
2007 
2008  dst.equ(relaxation, src);
2009 }
2010 
2011 
2012 
2013 template <class VectorType>
2014 inline void
2015 PreconditionRichardson::Tvmult(VectorType &dst, const VectorType &src) const
2016 {
2017  static_assert(
2018  std::is_same<size_type, typename VectorType::size_type>::value,
2019  "PreconditionRichardson and VectorType must have the same size_type.");
2020 
2021  dst.equ(relaxation, src);
2022 }
2023 
2024 template <class VectorType>
2025 inline void
2026 PreconditionRichardson::vmult_add(VectorType &dst, const VectorType &src) const
2027 {
2028  static_assert(
2029  std::is_same<size_type, typename VectorType::size_type>::value,
2030  "PreconditionRichardson and VectorType must have the same size_type.");
2031 
2032  dst.add(relaxation, src);
2033 }
2034 
2035 
2036 
2037 template <class VectorType>
2038 inline void
2039 PreconditionRichardson::Tvmult_add(VectorType &dst, const VectorType &src) const
2040 {
2041  static_assert(
2042  std::is_same<size_type, typename VectorType::size_type>::value,
2043  "PreconditionRichardson and VectorType must have the same size_type.");
2044 
2045  dst.add(relaxation, src);
2046 }
2047 
2050 {
2051  Assert(n_rows != 0, ExcNotInitialized());
2052  return n_rows;
2053 }
2054 
2057 {
2059  return n_columns;
2060 }
2061 
2062 //---------------------------------------------------------------------------
2063 
2064 template <typename MatrixType, typename PreconditionerType>
2065 inline void
2067  const MatrixType & rA,
2068  const AdditionalData &parameters)
2069 {
2070  A = &rA;
2071  relaxation = parameters.relaxation;
2072 
2073  Assert(parameters.preconditioner, ExcNotInitialized());
2074 
2075  preconditioner = parameters.preconditioner;
2076  n_iterations = parameters.n_iterations;
2077 }
2078 
2079 
2080 template <typename MatrixType, typename PreconditionerType>
2081 inline void
2083 {
2084  A = nullptr;
2085  preconditioner = nullptr;
2086 }
2087 
2088 template <typename MatrixType, typename PreconditionerType>
2089 inline
2092 {
2093  Assert(A != nullptr, ExcNotInitialized());
2094  return A->m();
2095 }
2096 
2097 template <typename MatrixType, typename PreconditionerType>
2098 inline
2101 {
2102  Assert(A != nullptr, ExcNotInitialized());
2103  return A->n();
2104 }
2105 
2106 template <typename MatrixType, typename PreconditionerType>
2107 template <class VectorType>
2108 inline void
2110  VectorType & dst,
2111  const VectorType &src) const
2112 {
2113  Assert(this->A != nullptr, ExcNotInitialized());
2114  Assert(this->preconditioner != nullptr, ExcNotInitialized());
2115 
2116  VectorType tmp1, tmp2;
2117 
2118  for (unsigned int i = 0; i < n_iterations; ++i)
2119  internal::PreconditionRelaxation::step_operations(
2120  *A, *preconditioner, dst, src, relaxation, tmp1, tmp2, i, false);
2121 }
2122 
2123 template <typename MatrixType, typename PreconditionerType>
2124 template <class VectorType>
2125 inline void
2127  VectorType & dst,
2128  const VectorType &src) const
2129 {
2130  Assert(this->A != nullptr, ExcNotInitialized());
2131  Assert(this->preconditioner != nullptr, ExcNotInitialized());
2132 
2133  VectorType tmp1, tmp2;
2134 
2135  for (unsigned int i = 0; i < n_iterations; ++i)
2136  internal::PreconditionRelaxation::step_operations(
2137  *A, *preconditioner, dst, src, relaxation, tmp1, tmp2, i, true);
2138 }
2139 
2140 template <typename MatrixType, typename PreconditionerType>
2141 template <class VectorType>
2142 inline void
2144  VectorType & dst,
2145  const VectorType &src) const
2146 {
2147  Assert(this->A != nullptr, ExcNotInitialized());
2148  Assert(this->preconditioner != nullptr, ExcNotInitialized());
2149 
2150  VectorType tmp1, tmp2;
2151 
2152  for (unsigned int i = 1; i <= n_iterations; ++i)
2153  internal::PreconditionRelaxation::step_operations(
2154  *A, *preconditioner, dst, src, relaxation, tmp1, tmp2, i, false);
2155 }
2156 
2157 template <typename MatrixType, typename PreconditionerType>
2158 template <class VectorType>
2159 inline void
2161  VectorType & dst,
2162  const VectorType &src) const
2163 {
2164  Assert(this->A != nullptr, ExcNotInitialized());
2165  Assert(this->preconditioner != nullptr, ExcNotInitialized());
2166 
2167  VectorType tmp1, tmp2;
2168 
2169  for (unsigned int i = 1; i <= n_iterations; ++i)
2170  internal::PreconditionRelaxation::step_operations(
2171  *A, *preconditioner, dst, src, relaxation, tmp1, tmp2, i, true);
2172 }
2173 
2174 //---------------------------------------------------------------------------
2175 
2176 template <typename MatrixType>
2177 inline void
2179  const AdditionalData &parameters_in)
2180 {
2181  Assert(parameters_in.preconditioner == nullptr, ExcInternalError());
2182 
2183  AdditionalData parameters;
2184  parameters.relaxation = 1.0;
2185  parameters.n_iterations = parameters_in.n_iterations;
2186  parameters.preconditioner =
2187  std::make_shared<PreconditionerType>(A, parameters_in.relaxation);
2188 
2189  this->BaseClass::initialize(A, parameters);
2190 }
2191 
2192 //---------------------------------------------------------------------------
2193 
2194 template <typename MatrixType>
2195 inline void
2196 PreconditionSOR<MatrixType>::initialize(const MatrixType & A,
2197  const AdditionalData &parameters_in)
2198 {
2199  Assert(parameters_in.preconditioner == nullptr, ExcInternalError());
2200 
2201  AdditionalData parameters;
2202  parameters.relaxation = 1.0;
2203  parameters.n_iterations = parameters_in.n_iterations;
2204  parameters.preconditioner =
2205  std::make_shared<PreconditionerType>(A, parameters_in.relaxation);
2206 
2207  this->BaseClass::initialize(A, parameters);
2208 }
2209 
2210 //---------------------------------------------------------------------------
2211 
2212 template <typename MatrixType>
2213 inline void
2214 PreconditionSSOR<MatrixType>::initialize(const MatrixType & A,
2215  const AdditionalData &parameters_in)
2216 {
2217  Assert(parameters_in.preconditioner == nullptr, ExcInternalError());
2218 
2219  AdditionalData parameters;
2220  parameters.relaxation = 1.0;
2221  parameters.n_iterations = parameters_in.n_iterations;
2222  parameters.preconditioner =
2223  std::make_shared<PreconditionerType>(A, parameters_in.relaxation);
2224 
2225  this->BaseClass::initialize(A, parameters);
2226 }
2227 
2228 
2229 
2230 //---------------------------------------------------------------------------
2231 
2232 template <typename MatrixType>
2233 inline void
2235  const MatrixType & A,
2236  const std::vector<size_type> & p,
2237  const std::vector<size_type> & ip,
2238  const typename BaseClass::AdditionalData &parameters_in)
2239 {
2240  Assert(parameters_in.preconditioner == nullptr, ExcInternalError());
2241 
2242  typename BaseClass::AdditionalData parameters;
2243  parameters.relaxation = 1.0;
2244  parameters.n_iterations = parameters_in.n_iterations;
2245  parameters.preconditioner =
2246  std::make_shared<PreconditionerType>(A, parameters_in.relaxation, p, ip);
2247 
2248  this->BaseClass::initialize(A, parameters);
2249 }
2250 
2251 
2252 template <typename MatrixType>
2253 inline void
2254 PreconditionPSOR<MatrixType>::initialize(const MatrixType & A,
2255  const AdditionalData &additional_data)
2256 {
2257  initialize(A,
2258  additional_data.permutation,
2259  additional_data.inverse_permutation,
2260  additional_data.parameters);
2261 }
2262 
2263 template <typename MatrixType>
2265  const std::vector<size_type> &permutation,
2266  const std::vector<size_type> &inverse_permutation,
2268  &parameters)
2269  : permutation(permutation)
2270  , inverse_permutation(inverse_permutation)
2271  , parameters(parameters)
2272 {}
2273 
2274 
2275 //---------------------------------------------------------------------------
2276 
2277 
2278 template <typename MatrixType, class VectorType>
2280  const MatrixType & M,
2281  const function_ptr method)
2282  : matrix(M)
2283  , precondition(method)
2284 {}
2285 
2286 
2287 
2288 template <typename MatrixType, class VectorType>
2289 void
2291  VectorType & dst,
2292  const VectorType &src) const
2293 {
2294  (matrix.*precondition)(dst, src);
2295 }
2296 
2297 //---------------------------------------------------------------------------
2298 
2299 template <typename MatrixType, typename PreconditionerType>
2301  AdditionalData(const double relaxation, const unsigned int n_iterations)
2302  : relaxation(relaxation)
2303  , n_iterations(n_iterations)
2304 {}
2305 
2306 
2307 
2308 //---------------------------------------------------------------------------
2309 
2310 namespace internal
2311 {
2312  namespace PreconditionChebyshevImplementation
2313  {
2314  // for deal.II vectors, perform updates for Chebyshev preconditioner all
2315  // at once to reduce memory transfer. Here, we select between general
2316  // vectors and deal.II vectors where we expand the loop over the (local)
2317  // size of the vector
2318 
2319  // generic part for non-deal.II vectors
2320  template <typename VectorType, typename PreconditionerType>
2321  inline void
2322  vector_updates(const VectorType & rhs,
2323  const PreconditionerType &preconditioner,
2324  const unsigned int iteration_index,
2325  const double factor1,
2326  const double factor2,
2327  VectorType & solution_old,
2328  VectorType & temp_vector1,
2329  VectorType & temp_vector2,
2330  VectorType & solution)
2331  {
2332  if (iteration_index == 0)
2333  {
2334  solution.equ(factor2, rhs);
2335  preconditioner.vmult(solution_old, solution);
2336  }
2337  else if (iteration_index == 1)
2338  {
2339  // compute t = P^{-1} * (b-A*x^{n})
2340  temp_vector1.sadd(-1.0, 1.0, rhs);
2341  preconditioner.vmult(solution_old, temp_vector1);
2342 
2343  // compute x^{n+1} = x^{n} + f_1 * x^{n} + f_2 * t
2344  solution_old.sadd(factor2, 1 + factor1, solution);
2345  }
2346  else
2347  {
2348  // compute t = P^{-1} * (b-A*x^{n})
2349  temp_vector1.sadd(-1.0, 1.0, rhs);
2350  preconditioner.vmult(temp_vector2, temp_vector1);
2351 
2352  // compute x^{n+1} = x^{n} + f_1 * (x^{n}-x^{n-1}) + f_2 * t
2353  solution_old.sadd(-factor1, factor2, temp_vector2);
2354  solution_old.add(1 + factor1, solution);
2355  }
2356 
2357  solution.swap(solution_old);
2358  }
2359 
2360  // worker routine for deal.II vectors. Because of vectorization, we need
2361  // to put the loop into an extra structure because the virtual function of
2362  // VectorUpdatesRange prevents the compiler from applying vectorization.
2363  template <typename Number>
2364  struct VectorUpdater
2365  {
2366  VectorUpdater(const Number * rhs,
2367  const Number * matrix_diagonal_inverse,
2368  const unsigned int iteration_index,
2369  const Number factor1,
2370  const Number factor2,
2371  Number * solution_old,
2372  Number * tmp_vector,
2373  Number * solution)
2374  : rhs(rhs)
2375  , matrix_diagonal_inverse(matrix_diagonal_inverse)
2376  , iteration_index(iteration_index)
2377  , factor1(factor1)
2378  , factor2(factor2)
2379  , solution_old(solution_old)
2380  , tmp_vector(tmp_vector)
2381  , solution(solution)
2382  {}
2383 
2384  void
2385  apply_to_subrange(const std::size_t begin, const std::size_t end) const
2386  {
2387  // To circumvent a bug in gcc
2388  // (https://gcc.gnu.org/bugzilla/show_bug.cgi?id=63945), we create
2389  // copies of the variables factor1 and factor2 and do not check based on
2390  // factor1.
2391  const Number factor1 = this->factor1;
2392  const Number factor1_plus_1 = 1. + this->factor1;
2393  const Number factor2 = this->factor2;
2394  if (iteration_index == 0)
2395  {
2397  for (std::size_t i = begin; i < end; ++i)
2398  solution[i] = factor2 * matrix_diagonal_inverse[i] * rhs[i];
2399  }
2400  else if (iteration_index == 1)
2401  {
2402  // x^{n+1} = x^{n} + f_1 * x^{n} + f_2 * P^{-1} * (b-A*x^{n})
2404  for (std::size_t i = begin; i < end; ++i)
2405  // for efficiency reason, write back to temp_vector that is
2406  // already read (avoid read-for-ownership)
2407  tmp_vector[i] =
2408  factor1_plus_1 * solution[i] +
2409  factor2 * matrix_diagonal_inverse[i] * (rhs[i] - tmp_vector[i]);
2410  }
2411  else
2412  {
2413  // x^{n+1} = x^{n} + f_1 * (x^{n}-x^{n-1})
2414  // + f_2 * P^{-1} * (b-A*x^{n})
2416  for (std::size_t i = begin; i < end; ++i)
2417  solution_old[i] =
2418  factor1_plus_1 * solution[i] - factor1 * solution_old[i] +
2419  factor2 * matrix_diagonal_inverse[i] * (rhs[i] - tmp_vector[i]);
2420  }
2421  }
2422 
2423  const Number * rhs;
2424  const Number * matrix_diagonal_inverse;
2425  const unsigned int iteration_index;
2426  const Number factor1;
2427  const Number factor2;
2428  mutable Number * solution_old;
2429  mutable Number * tmp_vector;
2430  mutable Number * solution;
2431  };
2432 
2433  template <typename Number>
2434  struct VectorUpdatesRange : public ::parallel::ParallelForInteger
2435  {
2436  VectorUpdatesRange(const VectorUpdater<Number> &updater,
2437  const std::size_t size)
2438  : updater(updater)
2439  {
2441  VectorUpdatesRange::apply_to_subrange(0, size);
2442  else
2443  apply_parallel(
2444  0,
2445  size,
2447  }
2448 
2449  ~VectorUpdatesRange() override = default;
2450 
2451  virtual void
2452  apply_to_subrange(const std::size_t begin,
2453  const std::size_t end) const override
2454  {
2455  updater.apply_to_subrange(begin, end);
2456  }
2457 
2458  const VectorUpdater<Number> &updater;
2459  };
2460 
2461  // selection for diagonal matrix around deal.II vector
2462  template <typename Number>
2463  inline void
2464  vector_updates(const ::Vector<Number> & rhs,
2465  const DiagonalMatrix<::Vector<Number>> &jacobi,
2466  const unsigned int iteration_index,
2467  const double factor1,
2468  const double factor2,
2469  ::Vector<Number> &solution_old,
2470  ::Vector<Number> &temp_vector1,
2471  ::Vector<Number> &,
2472  ::Vector<Number> &solution)
2473  {
2474  VectorUpdater<Number> upd(rhs.begin(),
2475  jacobi.get_vector().begin(),
2476  iteration_index,
2477  factor1,
2478  factor2,
2479  solution_old.begin(),
2480  temp_vector1.begin(),
2481  solution.begin());
2482  VectorUpdatesRange<Number>(upd, rhs.size());
2483 
2484  // swap vectors x^{n+1}->x^{n}, given the updates in the function above
2485  if (iteration_index == 0)
2486  {
2487  // nothing to do here because we can immediately write into the
2488  // solution vector without remembering any of the other vectors
2489  }
2490  else if (iteration_index == 1)
2491  {
2492  solution.swap(temp_vector1);
2493  solution_old.swap(temp_vector1);
2494  }
2495  else
2496  solution.swap(solution_old);
2497  }
2498 
2499  // selection for diagonal matrix around parallel deal.II vector
2500  template <typename Number>
2501  inline void
2502  vector_updates(
2504  const DiagonalMatrix<
2506  const unsigned int iteration_index,
2507  const double factor1,
2508  const double factor2,
2510  &solution_old,
2512  &temp_vector1,
2515  {
2516  VectorUpdater<Number> upd(rhs.begin(),
2517  jacobi.get_vector().begin(),
2518  iteration_index,
2519  factor1,
2520  factor2,
2521  solution_old.begin(),
2522  temp_vector1.begin(),
2523  solution.begin());
2524  VectorUpdatesRange<Number>(upd, rhs.locally_owned_size());
2525 
2526  // swap vectors x^{n+1}->x^{n}, given the updates in the function above
2527  if (iteration_index == 0)
2528  {
2529  // nothing to do here because we can immediately write into the
2530  // solution vector without remembering any of the other vectors
2531  }
2532  else if (iteration_index == 1)
2533  {
2534  solution.swap(temp_vector1);
2535  solution_old.swap(temp_vector1);
2536  }
2537  else
2538  solution.swap(solution_old);
2539  }
2540 
2541  // a helper type-trait that leverage SFINAE to figure out if MatrixType has
2542  // ... MatrixType::vmult(VectorType &, const VectorType&,
2543  // std::function<...>, std::function<...>) const
2544  template <typename MatrixType, typename VectorType>
2545  using vmult_functions_t = decltype(std::declval<MatrixType const>().vmult(
2546  std::declval<VectorType &>(),
2547  std::declval<const VectorType &>(),
2548  std::declval<
2549  const std::function<void(const unsigned int, const unsigned int)> &>(),
2550  std::declval<const std::function<void(const unsigned int,
2551  const unsigned int)> &>()));
2552 
2553  template <typename MatrixType,
2554  typename VectorType,
2555  typename PreconditionerType>
2556  constexpr bool has_vmult_with_std_functions =
2557  is_supported_operation<vmult_functions_t, MatrixType, VectorType>
2558  && std::is_same<PreconditionerType, DiagonalMatrix<VectorType>>::value
2559  &&std::is_same<
2560  VectorType,
2561  LinearAlgebra::distributed::Vector<typename VectorType::value_type,
2562  MemorySpace::Host>>::value;
2563 
2564  // We need to have a separate declaration for static const members
2565 
2566  template <
2567  typename MatrixType,
2568  typename VectorType,
2569  typename PreconditionerType,
2570  typename std::enable_if<!has_vmult_with_std_functions<MatrixType,
2571  VectorType,
2572  PreconditionerType>,
2573  int>::type * = nullptr>
2574  inline void
2575  vmult_and_update(const MatrixType & matrix,
2576  const PreconditionerType &preconditioner,
2577  const VectorType & rhs,
2578  const unsigned int iteration_index,
2579  const double factor1,
2580  const double factor2,
2581  VectorType & solution,
2582  VectorType & solution_old,
2583  VectorType & temp_vector1,
2584  VectorType & temp_vector2)
2585  {
2586  if (iteration_index > 0)
2587  matrix.vmult(temp_vector1, solution);
2588  vector_updates(rhs,
2589  preconditioner,
2590  iteration_index,
2591  factor1,
2592  factor2,
2593  solution_old,
2594  temp_vector1,
2595  temp_vector2,
2596  solution);
2597  }
2598 
2599  template <
2600  typename MatrixType,
2601  typename VectorType,
2602  typename PreconditionerType,
2603  typename std::enable_if<has_vmult_with_std_functions<MatrixType,
2604  VectorType,
2605  PreconditionerType>,
2606  int>::type * = nullptr>
2607  inline void
2608  vmult_and_update(const MatrixType & matrix,
2609  const PreconditionerType &preconditioner,
2610  const VectorType & rhs,
2611  const unsigned int iteration_index,
2612  const double factor1,
2613  const double factor2,
2614  VectorType & solution,
2615  VectorType & solution_old,
2616  VectorType & temp_vector1,
2617  VectorType &)
2618  {
2619  using Number = typename VectorType::value_type;
2620  VectorUpdater<Number> updater(rhs.begin(),
2621  preconditioner.get_vector().begin(),
2622  iteration_index,
2623  factor1,
2624  factor2,
2625  solution_old.begin(),
2626  temp_vector1.begin(),
2627  solution.begin());
2628  if (iteration_index > 0)
2629  matrix.vmult(
2630  temp_vector1,
2631  solution,
2632  [&](const unsigned int start_range, const unsigned int end_range) {
2633  // zero 'temp_vector1' before running the vmult
2634  // operation
2635  if (end_range > start_range)
2636  std::memset(temp_vector1.begin() + start_range,
2637  0,
2638  sizeof(Number) * (end_range - start_range));
2639  },
2640  [&](const unsigned int start_range, const unsigned int end_range) {
2641  if (end_range > start_range)
2642  updater.apply_to_subrange(start_range, end_range);
2643  });
2644  else
2645  updater.apply_to_subrange(0U, solution.locally_owned_size());
2646 
2647  // swap vectors x^{n+1}->x^{n}, given the updates in the function above
2648  if (iteration_index == 0)
2649  {
2650  // nothing to do here because we can immediately write into the
2651  // solution vector without remembering any of the other vectors
2652  }
2653  else if (iteration_index == 1)
2654  {
2655  solution.swap(temp_vector1);
2656  solution_old.swap(temp_vector1);
2657  }
2658  else
2659  solution.swap(solution_old);
2660  }
2661 
2662  template <typename MatrixType, typename PreconditionerType>
2663  inline void
2664  initialize_preconditioner(
2665  const MatrixType & matrix,
2666  std::shared_ptr<PreconditionerType> &preconditioner)
2667  {
2668  (void)matrix;
2669  (void)preconditioner;
2670  AssertThrow(preconditioner.get() != nullptr, ExcNotInitialized());
2671  }
2672 
2673  template <typename MatrixType, typename VectorType>
2674  inline void
2675  initialize_preconditioner(
2676  const MatrixType & matrix,
2677  std::shared_ptr<DiagonalMatrix<VectorType>> &preconditioner)
2678  {
2679  if (preconditioner.get() == nullptr || preconditioner->m() != matrix.m())
2680  {
2681  if (preconditioner.get() == nullptr)
2682  preconditioner = std::make_shared<DiagonalMatrix<VectorType>>();
2683 
2684  Assert(
2685  preconditioner->m() == 0,
2686  ExcMessage(
2687  "Preconditioner appears to be initialized but not sized correctly"));
2688 
2689  // This part only works in serial
2690  if (preconditioner->m() != matrix.m())
2691  {
2692  preconditioner->get_vector().reinit(matrix.m());
2693  for (typename VectorType::size_type i = 0; i < matrix.m(); ++i)
2694  preconditioner->get_vector()(i) = 1. / matrix.el(i, i);
2695  }
2696  }
2697  }
2698 
2699  template <typename VectorType>
2700  void
2701  set_initial_guess(VectorType &vector)
2702  {
2703  vector = 1. / std::sqrt(static_cast<double>(vector.size()));
2704  if (vector.locally_owned_elements().is_element(0))
2705  vector(0) = 0.;
2706  }
2707 
2708  template <typename Number>
2709  void
2710  set_initial_guess(::Vector<Number> &vector)
2711  {
2712  // Choose a high-frequency mode consisting of numbers between 0 and 1
2713  // that is cheap to compute (cheaper than random numbers) but avoids
2714  // obviously re-occurring numbers in multi-component systems by choosing
2715  // a period of 11
2716  for (unsigned int i = 0; i < vector.size(); ++i)
2717  vector(i) = i % 11;
2718 
2719  const Number mean_value = vector.mean_value();
2720  vector.add(-mean_value);
2721  }
2722 
2723  template <typename Number>
2724  void
2725  set_initial_guess(
2727  &vector)
2728  {
2729  // Choose a high-frequency mode consisting of numbers between 0 and 1
2730  // that is cheap to compute (cheaper than random numbers) but avoids
2731  // obviously re-occurring numbers in multi-component systems by choosing
2732  // a period of 11.
2733  // Make initial guess robust with respect to number of processors
2734  // by operating on the global index.
2735  types::global_dof_index first_local_range = 0;
2736  if (!vector.locally_owned_elements().is_empty())
2737  first_local_range = vector.locally_owned_elements().nth_index_in_set(0);
2738  for (unsigned int i = 0; i < vector.locally_owned_size(); ++i)
2739  vector.local_element(i) = (i + first_local_range) % 11;
2740 
2741  const Number mean_value = vector.mean_value();
2742  vector.add(-mean_value);
2743  }
2744 
2745  template <typename Number>
2746  void
2747  set_initial_guess(
2749  {
2750  for (unsigned int block = 0; block < vector.n_blocks(); ++block)
2751  set_initial_guess(vector.block(block));
2752  }
2753 
2754 
2755 # ifdef DEAL_II_COMPILER_CUDA_AWARE
2756  template <typename Number>
2757  __global__ void
2758  set_initial_guess_kernel(const types::global_dof_index offset,
2759  const unsigned int locally_owned_size,
2760  Number * values)
2761 
2762  {
2763  const unsigned int index = threadIdx.x + blockDim.x * blockIdx.x;
2764  if (index < locally_owned_size)
2765  values[index] = (index + offset) % 11;
2766  }
2767 
2768  template <typename Number>
2769  void
2770  set_initial_guess(
2772  &vector)
2773  {
2774  // Choose a high-frequency mode consisting of numbers between 0 and 1
2775  // that is cheap to compute (cheaper than random numbers) but avoids
2776  // obviously re-occurring numbers in multi-component systems by choosing
2777  // a period of 11.
2778  // Make initial guess robust with respect to number of processors
2779  // by operating on the global index.
2780  types::global_dof_index first_local_range = 0;
2781  if (!vector.locally_owned_elements().is_empty())
2782  first_local_range = vector.locally_owned_elements().nth_index_in_set(0);
2783 
2784  const auto n_local_elements = vector.locally_owned_size();
2785  const int n_blocks =
2786  1 + (n_local_elements - 1) / CUDAWrappers::block_size;
2787  set_initial_guess_kernel<<<n_blocks, CUDAWrappers::block_size>>>(
2788  first_local_range, n_local_elements, vector.get_values());
2789  AssertCudaKernel();
2790 
2791  const Number mean_value = vector.mean_value();
2792  vector.add(-mean_value);
2793  }
2794 # endif // DEAL_II_COMPILER_CUDA_AWARE
2795 
2796  struct EigenvalueTracker
2797  {
2798  public:
2799  void
2800  slot(const std::vector<double> &eigenvalues)
2801  {
2802  values = eigenvalues;
2803  }
2804 
2805  std::vector<double> values;
2806  };
2807 
2808 
2809 
2810  template <typename MatrixType,
2811  typename VectorType,
2812  typename PreconditionerType>
2813  double
2814  power_iteration(const MatrixType & matrix,
2815  VectorType & eigenvector,
2816  const PreconditionerType &preconditioner,
2817  const unsigned int n_iterations)
2818  {
2819  double eigenvalue_estimate = 0.;
2820  eigenvector /= eigenvector.l2_norm();
2821  VectorType vector1, vector2;
2822  vector1.reinit(eigenvector, true);
2823  if (!std::is_same<PreconditionerType, PreconditionIdentity>::value)
2824  vector2.reinit(eigenvector, true);
2825  for (unsigned int i = 0; i < n_iterations; ++i)
2826  {
2827  if (!std::is_same<PreconditionerType, PreconditionIdentity>::value)
2828  {
2829  matrix.vmult(vector2, eigenvector);
2830  preconditioner.vmult(vector1, vector2);
2831  }
2832  else
2833  matrix.vmult(vector1, eigenvector);
2834 
2835  eigenvalue_estimate = eigenvector * vector1;
2836 
2837  vector1 /= vector1.l2_norm();
2838  eigenvector.swap(vector1);
2839  }
2840  return eigenvalue_estimate;
2841  }
2842  } // namespace PreconditionChebyshevImplementation
2843 } // namespace internal
2844 
2845 
2846 
2847 template <typename MatrixType, class VectorType, typename PreconditionerType>
2849  AdditionalData::AdditionalData(const unsigned int degree,
2850  const double smoothing_range,
2851  const unsigned int eig_cg_n_iterations,
2852  const double eig_cg_residual,
2853  const double max_eigenvalue,
2854  const EigenvalueAlgorithm eigenvalue_algorithm)
2855  : degree(degree)
2856  , smoothing_range(smoothing_range)
2857  , eig_cg_n_iterations(eig_cg_n_iterations)
2858  , eig_cg_residual(eig_cg_residual)
2859  , max_eigenvalue(max_eigenvalue)
2860  , eigenvalue_algorithm(eigenvalue_algorithm)
2861 {}
2862 
2863 
2864 
2865 template <typename MatrixType, class VectorType, typename PreconditionerType>
2866 inline typename PreconditionChebyshev<MatrixType,
2867  VectorType,
2868  PreconditionerType>::AdditionalData &
2870  AdditionalData::operator=(const AdditionalData &other_data)
2871 {
2872  degree = other_data.degree;
2873  smoothing_range = other_data.smoothing_range;
2874  eig_cg_n_iterations = other_data.eig_cg_n_iterations;
2875  eig_cg_residual = other_data.eig_cg_residual;
2876  max_eigenvalue = other_data.max_eigenvalue;
2877  preconditioner = other_data.preconditioner;
2878  eigenvalue_algorithm = other_data.eigenvalue_algorithm;
2879  constraints.copy_from(other_data.constraints);
2880 
2881  return *this;
2882 }
2883 
2884 
2885 
2886 template <typename MatrixType, typename VectorType, typename PreconditionerType>
2889  : theta(1.)
2890  , delta(1.)
2891  , eigenvalues_are_initialized(false)
2892 {
2893  static_assert(
2894  std::is_same<size_type, typename VectorType::size_type>::value,
2895  "PreconditionChebyshev and VectorType must have the same size_type.");
2896 }
2897 
2898 
2899 
2900 template <typename MatrixType, typename VectorType, typename PreconditionerType>
2901 inline void
2903  const MatrixType & matrix,
2904  const AdditionalData &additional_data)
2905 {
2906  matrix_ptr = &matrix;
2907  data = additional_data;
2908  Assert(data.degree > 0,
2909  ExcMessage("The degree of the Chebyshev method must be positive."));
2910  internal::PreconditionChebyshevImplementation::initialize_preconditioner(
2911  matrix, data.preconditioner);
2912  eigenvalues_are_initialized = false;
2913 }
2914 
2915 
2916 
2917 template <typename MatrixType, typename VectorType, typename PreconditionerType>
2918 inline void
2920 {
2921  eigenvalues_are_initialized = false;
2922  theta = delta = 1.0;
2923  matrix_ptr = nullptr;
2924  {
2925  VectorType empty_vector;
2926  solution_old.reinit(empty_vector);
2927  temp_vector1.reinit(empty_vector);
2928  temp_vector2.reinit(empty_vector);
2929  }
2930  data.preconditioner.reset();
2931 }
2932 
2933 
2934 
2935 template <typename MatrixType, typename VectorType, typename PreconditionerType>
2936 inline typename PreconditionChebyshev<MatrixType,
2937  VectorType,
2938  PreconditionerType>::EigenvalueInformation
2940  estimate_eigenvalues(const VectorType &src) const
2941 {
2942  Assert(eigenvalues_are_initialized == false, ExcInternalError());
2943  Assert(data.preconditioner.get() != nullptr, ExcNotInitialized());
2944 
2946  EigenvalueInformation info{};
2947 
2948  solution_old.reinit(src);
2949  temp_vector1.reinit(src, true);
2950 
2951  if (data.eig_cg_n_iterations > 0)
2952  {
2953  Assert(data.eig_cg_n_iterations > 2,
2954  ExcMessage(
2955  "Need to set at least two iterations to find eigenvalues."));
2956 
2957  internal::PreconditionChebyshevImplementation::EigenvalueTracker
2958  eigenvalue_tracker;
2959 
2960  // set an initial guess that contains some high-frequency parts (to the
2961  // extent possible without knowing the discretization and the numbering)
2962  // to trigger high eigenvalues according to the external function
2963  internal::PreconditionChebyshevImplementation::set_initial_guess(
2964  temp_vector1);
2965  data.constraints.set_zero(temp_vector1);
2966 
2967  if (data.eigenvalue_algorithm ==
2968  AdditionalData::EigenvalueAlgorithm::lanczos)
2969  {
2970  // set a very strict tolerance to force at least two iterations
2971  IterationNumberControl control(data.eig_cg_n_iterations,
2972  1e-10,
2973  false,
2974  false);
2975 
2976  SolverCG<VectorType> solver(control);
2977  solver.connect_eigenvalues_slot(
2978  [&eigenvalue_tracker](const std::vector<double> &eigenvalues) {
2979  eigenvalue_tracker.slot(eigenvalues);
2980  });
2981 
2982  solver.solve(*matrix_ptr,
2983  solution_old,
2984  temp_vector1,
2985  *data.preconditioner);
2986 
2987  info.cg_iterations = control.last_step();
2988  }
2989  else if (data.eigenvalue_algorithm ==
2990  AdditionalData::EigenvalueAlgorithm::power_iteration)
2991  {
2992  Assert(data.degree != numbers::invalid_unsigned_int,
2993  ExcMessage("Cannot estimate the minimal eigenvalue with the "
2994  "power iteration"));
2995 
2996  eigenvalue_tracker.values.push_back(
2997  internal::PreconditionChebyshevImplementation::power_iteration(
2998  *matrix_ptr,
2999  temp_vector1,
3000  *data.preconditioner,
3001  data.eig_cg_n_iterations));
3002  }
3003  else
3004  Assert(false, ExcNotImplemented());
3005 
3006  // read the eigenvalues from the attached eigenvalue tracker
3007  if (eigenvalue_tracker.values.empty())
3008  info.min_eigenvalue_estimate = info.max_eigenvalue_estimate = 1.;
3009  else
3010  {
3011  info.min_eigenvalue_estimate = eigenvalue_tracker.values.front();
3012 
3013  // include a safety factor since the CG method will in general not
3014  // be converged
3015  info.max_eigenvalue_estimate = 1.2 * eigenvalue_tracker.values.back();
3016  }
3017  }
3018  else
3019  {
3020  info.max_eigenvalue_estimate = data.max_eigenvalue;
3021  info.min_eigenvalue_estimate = data.max_eigenvalue / data.smoothing_range;
3022  }
3023 
3024  const double alpha = (data.smoothing_range > 1. ?
3025  info.max_eigenvalue_estimate / data.smoothing_range :
3026  std::min(0.9 * info.max_eigenvalue_estimate,
3027  info.min_eigenvalue_estimate));
3028 
3029  // in case the user set the degree to invalid unsigned int, we have to
3030  // determine the number of necessary iterations from the Chebyshev error
3031  // estimate, given the target tolerance specified by smoothing_range. This
3032  // estimate is based on the error formula given in section 5.1 of
3033  // R. S. Varga, Matrix iterative analysis, 2nd ed., Springer, 2009
3034  if (data.degree == numbers::invalid_unsigned_int)
3035  {
3036  const double actual_range = info.max_eigenvalue_estimate / alpha;
3037  const double sigma = (1. - std::sqrt(1. / actual_range)) /
3038  (1. + std::sqrt(1. / actual_range));
3039  const double eps = data.smoothing_range;
3040  const_cast<
3042  this)
3043  ->data.degree =
3044  1 + static_cast<unsigned int>(
3045  std::log(1. / eps + std::sqrt(1. / eps / eps - 1.)) /
3046  std::log(1. / sigma));
3047  }
3048 
3049  info.degree = data.degree;
3050 
3051  const_cast<
3053  ->delta = (info.max_eigenvalue_estimate - alpha) * 0.5;
3054  const_cast<
3056  ->theta = (info.max_eigenvalue_estimate + alpha) * 0.5;
3057 
3058  // We do not need the second temporary vector in case we have a
3059  // DiagonalMatrix as preconditioner and use deal.II's own vectors
3060  using NumberType = typename VectorType::value_type;
3061  if (std::is_same<PreconditionerType, DiagonalMatrix<VectorType>>::value ==
3062  false ||
3063  (std::is_same<VectorType, ::Vector<NumberType>>::value == false &&
3064  ((std::is_same<VectorType,
3067  false) ||
3068  (std::is_same<VectorType,
3069  LinearAlgebra::distributed::
3070  Vector<NumberType, MemorySpace::CUDA>>::value ==
3071  false))))
3072  temp_vector2.reinit(src, true);
3073  else
3074  {
3075  VectorType empty_vector;
3076  temp_vector2.reinit(empty_vector);
3077  }
3078 
3079  const_cast<
3081  ->eigenvalues_are_initialized = true;
3082 
3083  return info;
3084 }
3085 
3086 
3087 
3088 template <typename MatrixType, typename VectorType, typename PreconditionerType>
3089 inline void
3091  VectorType & solution,
3092  const VectorType &rhs) const
3093 {
3094  std::lock_guard<std::mutex> lock(mutex);
3095  if (eigenvalues_are_initialized == false)
3096  estimate_eigenvalues(rhs);
3097 
3098  internal::PreconditionChebyshevImplementation::vmult_and_update(
3099  *matrix_ptr,
3100  *data.preconditioner,
3101  rhs,
3102  0,
3103  0.,
3104  1. / theta,
3105  solution,
3106  solution_old,
3107  temp_vector1,
3108  temp_vector2);
3109 
3110  // if delta is zero, we do not need to iterate because the updates will be
3111  // zero
3112  if (data.degree < 2 || std::abs(delta) < 1e-40)
3113  return;
3114 
3115  double rhok = delta / theta, sigma = theta / delta;
3116  for (unsigned int k = 0; k < data.degree - 1; ++k)
3117  {
3118  const double rhokp = 1. / (2. * sigma - rhok);
3119  const double factor1 = rhokp * rhok, factor2 = 2. * rhokp / delta;
3120  rhok = rhokp;
3121  internal::PreconditionChebyshevImplementation::vmult_and_update(
3122  *matrix_ptr,
3123  *data.preconditioner,
3124  rhs,
3125  k + 1,
3126  factor1,
3127  factor2,
3128  solution,
3129  solution_old,
3130  temp_vector1,
3131  temp_vector2);
3132  }
3133 }
3134 
3135 
3136 
3137 template <typename MatrixType, typename VectorType, typename PreconditionerType>
3138 inline void
3140  VectorType & solution,
3141  const VectorType &rhs) const
3142 {
3143  std::lock_guard<std::mutex> lock(mutex);
3144  if (eigenvalues_are_initialized == false)
3145  estimate_eigenvalues(rhs);
3146 
3147  internal::PreconditionChebyshevImplementation::vector_updates(
3148  rhs,
3149  *data.preconditioner,
3150  0,
3151  0.,
3152  1. / theta,
3153  solution_old,
3154  temp_vector1,
3155  temp_vector2,
3156  solution);
3157 
3158  if (data.degree < 2 || std::abs(delta) < 1e-40)
3159  return;
3160 
3161  double rhok = delta / theta, sigma = theta / delta;
3162  for (unsigned int k = 0; k < data.degree - 1; ++k)
3163  {
3164  const double rhokp = 1. / (2. * sigma - rhok);
3165  const double factor1 = rhokp * rhok, factor2 = 2. * rhokp / delta;
3166  rhok = rhokp;
3167  matrix_ptr->Tvmult(temp_vector1, solution);
3168  internal::PreconditionChebyshevImplementation::vector_updates(
3169  rhs,
3170  *data.preconditioner,
3171  k + 1,
3172  factor1,
3173  factor2,
3174  solution_old,
3175  temp_vector1,
3176  temp_vector2,
3177  solution);
3178  }
3179 }
3180 
3181 
3182 
3183 template <typename MatrixType, typename VectorType, typename PreconditionerType>
3184 inline void
3186  VectorType & solution,
3187  const VectorType &rhs) const
3188 {
3189  std::lock_guard<std::mutex> lock(mutex);
3190  if (eigenvalues_are_initialized == false)
3191  estimate_eigenvalues(rhs);
3192 
3193  internal::PreconditionChebyshevImplementation::vmult_and_update(
3194  *matrix_ptr,
3195  *data.preconditioner,
3196  rhs,
3197  1,
3198  0.,
3199  1. / theta,
3200  solution,
3201  solution_old,
3202  temp_vector1,
3203  temp_vector2);
3204 
3205  if (data.degree < 2 || std::abs(delta) < 1e-40)
3206  return;
3207 
3208  double rhok = delta / theta, sigma = theta / delta;
3209  for (unsigned int k = 0; k < data.degree - 1; ++k)
3210  {
3211  const double rhokp = 1. / (2. * sigma - rhok);
3212  const double factor1 = rhokp * rhok, factor2 = 2. * rhokp / delta;
3213  rhok = rhokp;
3214  internal::PreconditionChebyshevImplementation::vmult_and_update(
3215  *matrix_ptr,
3216  *data.preconditioner,
3217  rhs,
3218  k + 2,
3219  factor1,
3220  factor2,
3221  solution,
3222  solution_old,
3223  temp_vector1,
3224  temp_vector2);
3225  }
3226 }
3227 
3228 
3229 
3230 template <typename MatrixType, typename VectorType, typename PreconditionerType>
3231 inline void
3233  VectorType & solution,
3234  const VectorType &rhs) const
3235 {
3236  std::lock_guard<std::mutex> lock(mutex);
3237  if (eigenvalues_are_initialized == false)
3238  estimate_eigenvalues(rhs);
3239 
3240  matrix_ptr->Tvmult(temp_vector1, solution);
3241  internal::PreconditionChebyshevImplementation::vector_updates(
3242  rhs,
3243  *data.preconditioner,
3244  1,
3245  0.,
3246  1. / theta,
3247  solution_old,
3248  temp_vector1,
3249  temp_vector2,
3250  solution);
3251 
3252  if (data.degree < 2 || std::abs(delta) < 1e-40)
3253  return;
3254 
3255  double rhok = delta / theta, sigma = theta / delta;
3256  for (unsigned int k = 0; k < data.degree - 1; ++k)
3257  {
3258  const double rhokp = 1. / (2. * sigma - rhok);
3259  const double factor1 = rhokp * rhok, factor2 = 2. * rhokp / delta;
3260  rhok = rhokp;
3261  matrix_ptr->Tvmult(temp_vector1, solution);
3262  internal::PreconditionChebyshevImplementation::vector_updates(
3263  rhs,
3264  *data.preconditioner,
3265  k + 2,
3266  factor1,
3267  factor2,
3268  solution_old,
3269  temp_vector1,
3270  temp_vector2,
3271  solution);
3272  }
3273 }
3274 
3275 
3276 
3277 template <typename MatrixType, typename VectorType, typename PreconditionerType>
3278 inline typename PreconditionChebyshev<MatrixType,
3279  VectorType,
3280  PreconditionerType>::size_type
3282 {
3283  Assert(matrix_ptr != nullptr, ExcNotInitialized());
3284  return matrix_ptr->m();
3285 }
3286 
3287 
3288 
3289 template <typename MatrixType, typename VectorType, typename PreconditionerType>
3290 inline typename PreconditionChebyshev<MatrixType,
3291  VectorType,
3292  PreconditionerType>::size_type
3294 {
3295  Assert(matrix_ptr != nullptr, ExcNotInitialized());
3296  return matrix_ptr->n();
3297 }
3298 
3299 #endif // DOXYGEN
3300 
3302 
3303 #endif
bool is_empty() const
Definition: index_set.h:1775
size_type nth_index_in_set(const size_type local_index) const
Definition: index_set.h:1831
const_iterator end() const
const_iterator begin() const
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)
Definition: vector.h:109
#define DEAL_II_OPENMP_SIMD_PRAGMA
Definition: config.h:142
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
static ::ExceptionBase & ExcInternalError()
#define AssertCudaKernel()
Definition: exceptions.h:1871
unsigned int n_blocks() const
static ::ExceptionBase & ExcNotInitialized()
#define Assert(cond, exc)
Definition: exceptions.h:1473
static ::ExceptionBase & ExcNotImplemented()
BlockType & block(const unsigned int i)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1583
types::global_dof_index size_type
Definition: precondition.h:202
void vmult_add(VectorType &, const VectorType &) const
internal::PreconditionRelaxation::PreconditionSORImpl< MatrixType > PreconditionerType
typename BaseClass::AdditionalData AdditionalData
SmartPointer< const MatrixType, PreconditionRelaxation< MatrixType > > A
Definition: precondition.h:503
void Tvmult(VectorType &dst, const VectorType &src) const
AdditionalData(const std::vector< size_type > &permutation, const std::vector< size_type > &inverse_permutation, const typename BaseClass::AdditionalData &parameters=typename BaseClass::AdditionalData())
typename BaseClass::size_type size_type
void Tvmult(VectorType &, const VectorType &) const
std::shared_ptr< PreconditionerType > preconditioner
Definition: precondition.h:518
size_type m() const
Threads::Mutex mutex
size_type n() const
std::shared_ptr< PreconditionerType > preconditioner
void(MatrixType::*)(VectorType &, const VectorType &) const function_ptr
Definition: precondition.h:367
void vmult_add(VectorType &, const VectorType &) const
size_type n() const
size_type n() const
const function_ptr precondition
Definition: precondition.h:392
void step(VectorType &dst, const VectorType &src) const
AdditionalData data
void initialize(const MatrixType &A, const AdditionalData &parameters=AdditionalData())
BaseClass::AdditionalData parameters
void initialize(const AdditionalData &parameters)
void vmult(VectorType &, const VectorType &) const
size_type m() const
void initialize(const MatrixType &A, const AdditionalData &additional_data)
EigenvalueInformation estimate_eigenvalues(const VectorType &src) const
AffineConstraints< double > constraints
void Tstep(VectorType &dst, const VectorType &src) const
std::shared_ptr< PreconditionerType > preconditioner
Definition: precondition.h:438
const std::vector< size_type > & inverse_permutation
const MatrixType & matrix
Definition: precondition.h:387
typename BaseClass::AdditionalData AdditionalData
void initialize(const MatrixType &matrix, const AdditionalData &additional_data=AdditionalData())
void initialize(const MatrixType &matrix, const AdditionalData &parameters)
SmartPointer< const MatrixType, PreconditionChebyshev< MatrixType, VectorType, PreconditionerType > > matrix_ptr
void initialize(const MatrixType &A, const AdditionalData &parameters=AdditionalData())
size_type m() const
size_type m() const
AdditionalData(const double relaxation=1.)
internal::PreconditionRelaxation::PreconditionJacobiImpl< MatrixType > PreconditionerType
void initialize(const MatrixType &A, const AdditionalData &parameters=AdditionalData())
void step(VectorType &x, const VectorType &rhs) const
AdditionalData(const double relaxation=1., const unsigned int n_iterations=1)
unsigned int n_iterations
Definition: precondition.h:513
size_type n() const
types::global_dof_index size_type
const std::vector< size_type > & permutation
void vmult(VectorType &, const VectorType &) const
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)
void Tvmult(VectorType &, const VectorType &) const
void Tvmult_add(VectorType &, const VectorType &) const
void initialize(const MatrixType &A, const AdditionalData &parameters=AdditionalData())
void vmult(VectorType &dst, const VectorType &src) const
internal::PreconditionRelaxation::PreconditionPSORImpl< MatrixType > PreconditionerType
typename BaseClass::AdditionalData AdditionalData
internal::PreconditionRelaxation::PreconditionSSORImpl< MatrixType > PreconditionerType
void vmult(VectorType &dst, const VectorType &src) const
void Tstep(VectorType &x, const VectorType &rhs) const
void vmult(VectorType &, const VectorType &) const
void initialize(const MatrixType &matrix, const AdditionalData &additional_data=AdditionalData())
PreconditionUseMatrix(const MatrixType &M, const function_ptr method)
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())
void Tvmult(VectorType &, const VectorType &) const
EigenvalueAlgorithm eigenvalue_algorithm
types::global_dof_index size_type
Definition: precondition.h:87
void Tvmult_add(VectorType &, const VectorType &) const
AdditionalData & operator=(const AdditionalData &other_data)
types::global_dof_index size_type
Definition: precondition.h:410
void add(const std::vector< size_type > &indices, const std::vector< OtherNumber > &values)
Number local_element(const size_type local_index) const
void swap(Vector< Number, MemorySpace > &v)
Number mean_value() const
size_type locally_owned_size() const
virtual Number mean_value() const override
size_type size() const
virtual void swap(Vector< Number > &v)
virtual void add(const Number a) override
iterator begin()
virtual ::IndexSet locally_owned_elements() const override
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< IsBlockVector< VectorType >::value, unsigned int >::type 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:201
unsigned int global_dof_index
Definition: types.h:76
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)