Reference documentation for deal.II version Git 4fbb5374f0 2021-01-22 15:02:13 -0500
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
solver_gmres.h
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1998 - 2020 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #ifndef dealii_solver_gmres_h
17 #define dealii_solver_gmres_h
18 
19 
20 
21 #include <deal.II/base/config.h>
22 
23 #include <deal.II/base/logstream.h>
25 
29 #include <deal.II/lac/solver.h>
31 #include <deal.II/lac/vector.h>
32 
33 #include <algorithm>
34 #include <cmath>
35 #include <memory>
36 #include <vector>
37 
39 
42 
43 namespace internal
44 {
48  namespace SolverGMRESImplementation
49  {
58  template <typename VectorType>
59  class TmpVectors
60  {
61  public:
66  TmpVectors(const unsigned int max_size, VectorMemory<VectorType> &vmem);
67 
71  ~TmpVectors() = default;
72 
77  VectorType &operator[](const unsigned int i) const;
78 
85  VectorType &
86  operator()(const unsigned int i, const VectorType &temp);
87 
92  unsigned int
93  size() const;
94 
95 
96  private:
101 
105  std::vector<typename VectorMemory<VectorType>::Pointer> data;
106  };
107  } // namespace SolverGMRESImplementation
108 } // namespace internal
109 
174 template <class VectorType = Vector<double>>
175 class SolverGMRES : public SolverBase<VectorType>
176 {
177 public:
182  {
189  explicit AdditionalData(const unsigned int max_n_tmp_vectors = 30,
190  const bool right_preconditioning = false,
191  const bool use_default_residual = true,
192  const bool force_re_orthogonalization = false);
193 
200  unsigned int max_n_tmp_vectors;
201 
210 
215 
223  };
224 
230  const AdditionalData & data = AdditionalData());
231 
237 
241  SolverGMRES(const SolverGMRES<VectorType> &) = delete;
242 
246  template <typename MatrixType, typename PreconditionerType>
247  void
248  solve(const MatrixType & A,
249  VectorType & x,
250  const VectorType & b,
251  const PreconditionerType &preconditioner);
252 
259  boost::signals2::connection
260  connect_condition_number_slot(const std::function<void(double)> &slot,
261  const bool every_iteration = false);
262 
269  boost::signals2::connection
270  connect_eigenvalues_slot(
271  const std::function<void(const std::vector<std::complex<double>> &)> &slot,
272  const bool every_iteration = false);
273 
281  boost::signals2::connection
282  connect_hessenberg_slot(
283  const std::function<void(const FullMatrix<double> &)> &slot,
284  const bool every_iteration = true);
285 
292  boost::signals2::connection
293  connect_krylov_space_slot(
294  const std::function<
296  &slot);
297 
298 
303  boost::signals2::connection
304  connect_re_orthogonalization_slot(const std::function<void(int)> &slot);
305 
306 
307  DeclException1(ExcTooFewTmpVectors,
308  int,
309  << "The number of temporary vectors you gave (" << arg1
310  << ") is too small. It should be at least 10 for "
311  << "any results, and much more for reasonable ones.");
312 
313 protected:
318 
323  boost::signals2::signal<void(double)> condition_number_signal;
324 
329  boost::signals2::signal<void(double)> all_condition_numbers_signal;
330 
335  boost::signals2::signal<void(const std::vector<std::complex<double>> &)>
336  eigenvalues_signal;
337 
342  boost::signals2::signal<void(const std::vector<std::complex<double>> &)>
343  all_eigenvalues_signal;
344 
349  boost::signals2::signal<void(const FullMatrix<double> &)> hessenberg_signal;
350 
355  boost::signals2::signal<void(const FullMatrix<double> &)>
356  all_hessenberg_signal;
357 
362  boost::signals2::signal<void(
365 
370  boost::signals2::signal<void(int)> re_orthogonalize_signal;
371 
375  virtual double
376  criterion();
377 
382  void
383  givens_rotation(Vector<double> &h,
384  Vector<double> &b,
385  Vector<double> &ci,
386  Vector<double> &si,
387  int col) const;
388 
399  static double
400  modified_gram_schmidt(
402  & orthogonal_vectors,
403  const unsigned int dim,
404  const unsigned int accumulated_iterations,
405  VectorType & vv,
406  Vector<double> & h,
407  bool & re_orthogonalize,
408  const boost::signals2::signal<void(int)> &re_orthogonalize_signal =
409  boost::signals2::signal<void(int)>());
410 
417  static void
418  compute_eigs_and_cond(
419  const FullMatrix<double> &H_orig,
420  const unsigned int dim,
421  const boost::signals2::signal<
422  void(const std::vector<std::complex<double>> &)> &eigenvalues_signal,
423  const boost::signals2::signal<void(const FullMatrix<double> &)>
424  & hessenberg_signal,
425  const boost::signals2::signal<void(double)> &cond_signal);
426 
431 
436 };
437 
458 template <class VectorType = Vector<double>>
459 class SolverFGMRES : public SolverBase<VectorType>
460 {
461 public:
466  {
470  explicit AdditionalData(const unsigned int max_basis_size = 30)
471  : max_basis_size(max_basis_size)
472  {}
473 
477  unsigned int max_basis_size;
478  };
479 
485  const AdditionalData & data = AdditionalData());
486 
492  const AdditionalData &data = AdditionalData());
493 
497  template <typename MatrixType, typename PreconditionerType>
498  void
499  solve(const MatrixType & A,
500  VectorType & x,
501  const VectorType & b,
502  const PreconditionerType &preconditioner);
503 
504 private:
509 
514 
519 };
520 
522 /* --------------------- Inline and template functions ------------------- */
523 
524 
525 #ifndef DOXYGEN
526 namespace internal
527 {
528  namespace SolverGMRESImplementation
529  {
530  template <class VectorType>
531  inline TmpVectors<VectorType>::TmpVectors(const unsigned int max_size,
533  : mem(vmem)
534  , data(max_size)
535  {}
536 
537 
538 
539  template <class VectorType>
541  operator[](const unsigned int i) const
542  {
543  AssertIndexRange(i, data.size());
544 
545  Assert(data[i] != nullptr, ExcNotInitialized());
546  return *data[i];
547  }
548 
549 
550 
551  template <class VectorType>
552  inline VectorType &
553  TmpVectors<VectorType>::operator()(const unsigned int i,
554  const VectorType & temp)
555  {
556  AssertIndexRange(i, data.size());
557  if (data[i] == nullptr)
558  {
559  data[i] = std::move(typename VectorMemory<VectorType>::Pointer(mem));
560  data[i]->reinit(temp);
561  }
562  return *data[i];
563  }
564 
565 
566 
567  template <class VectorType>
568  unsigned int
570  {
571  return (data.size() > 0 ? data.size() - 1 : 0);
572  }
573 
574 
575 
576  // A comparator for better printing eigenvalues
577  inline bool
578  complex_less_pred(const std::complex<double> &x,
579  const std::complex<double> &y)
580  {
581  return x.real() < y.real() ||
582  (x.real() == y.real() && x.imag() < y.imag());
583  }
584  } // namespace SolverGMRESImplementation
585 } // namespace internal
586 
587 
588 
589 template <class VectorType>
591  const unsigned int max_n_tmp_vectors,
592  const bool right_preconditioning,
593  const bool use_default_residual,
594  const bool force_re_orthogonalization)
595  : max_n_tmp_vectors(max_n_tmp_vectors)
596  , right_preconditioning(right_preconditioning)
597  , use_default_residual(use_default_residual)
598  , force_re_orthogonalization(force_re_orthogonalization)
599 {
600  Assert(3 <= max_n_tmp_vectors,
601  ExcMessage("SolverGMRES needs at least three "
602  "temporary vectors."));
603 }
604 
605 
606 
607 template <class VectorType>
610  const AdditionalData & data)
612  , additional_data(data)
613 {}
614 
615 
616 
617 template <class VectorType>
619  const AdditionalData &data)
621  , additional_data(data)
622 {}
623 
624 
625 
626 template <class VectorType>
627 inline void
629  Vector<double> &b,
630  Vector<double> &ci,
631  Vector<double> &si,
632  int col) const
633 {
634  for (int i = 0; i < col; i++)
635  {
636  const double s = si(i);
637  const double c = ci(i);
638  const double dummy = h(i);
639  h(i) = c * dummy + s * h(i + 1);
640  h(i + 1) = -s * dummy + c * h(i + 1);
641  };
642 
643  const double r = 1. / std::sqrt(h(col) * h(col) + h(col + 1) * h(col + 1));
644  si(col) = h(col + 1) * r;
645  ci(col) = h(col) * r;
646  h(col) = ci(col) * h(col) + si(col) * h(col + 1);
647  b(col + 1) = -si(col) * b(col);
648  b(col) *= ci(col);
649 }
650 
651 
652 
653 template <class VectorType>
654 inline double
657  & orthogonal_vectors,
658  const unsigned int dim,
659  const unsigned int accumulated_iterations,
660  VectorType & vv,
661  Vector<double> & h,
662  bool & reorthogonalize,
663  const boost::signals2::signal<void(int)> &reorthogonalize_signal)
664 {
665  Assert(dim > 0, ExcInternalError());
666  const unsigned int inner_iteration = dim - 1;
667 
668  // need initial norm for detection of re-orthogonalization, see below
669  double norm_vv_start = 0;
670  const bool consider_reorthogonalize =
671  (reorthogonalize == false) && (inner_iteration % 5 == 4);
672  if (consider_reorthogonalize)
673  norm_vv_start = vv.l2_norm();
674 
675  // Orthogonalization
676  h(0) = vv * orthogonal_vectors[0];
677  for (unsigned int i = 1; i < dim; ++i)
678  h(i) = vv.add_and_dot(-h(i - 1),
679  orthogonal_vectors[i - 1],
680  orthogonal_vectors[i]);
681  double norm_vv =
682  std::sqrt(vv.add_and_dot(-h(dim - 1), orthogonal_vectors[dim - 1], vv));
683 
684  // Re-orthogonalization if loss of orthogonality detected. For the test, use
685  // a strategy discussed in C. T. Kelley, Iterative Methods for Linear and
686  // Nonlinear Equations, SIAM, Philadelphia, 1995: Compare the norm of vv
687  // after orthogonalization with its norm when starting the
688  // orthogonalization. If vv became very small (here: less than the square
689  // root of the machine precision times 10), it is almost in the span of the
690  // previous vectors, which indicates loss of precision.
691  if (consider_reorthogonalize)
692  {
693  if (norm_vv >
694  10. * norm_vv_start *
695  std::sqrt(
697  return norm_vv;
698 
699  else
700  {
701  reorthogonalize = true;
702  if (!reorthogonalize_signal.empty())
703  reorthogonalize_signal(accumulated_iterations);
704  }
705  }
706 
707  if (reorthogonalize == true)
708  {
709  double htmp = vv * orthogonal_vectors[0];
710  h(0) += htmp;
711  for (unsigned int i = 1; i < dim; ++i)
712  {
713  htmp = vv.add_and_dot(-htmp,
714  orthogonal_vectors[i - 1],
715  orthogonal_vectors[i]);
716  h(i) += htmp;
717  }
718  norm_vv =
719  std::sqrt(vv.add_and_dot(-htmp, orthogonal_vectors[dim - 1], vv));
720  }
721 
722  return norm_vv;
723 }
724 
725 
726 
727 template <class VectorType>
728 inline void
730  const FullMatrix<double> &H_orig,
731  const unsigned int dim,
732  const boost::signals2::signal<void(const std::vector<std::complex<double>> &)>
733  &eigenvalues_signal,
734  const boost::signals2::signal<void(const FullMatrix<double> &)>
735  & hessenberg_signal,
736  const boost::signals2::signal<void(double)> &cond_signal)
737 {
738  // Avoid copying the Hessenberg matrix if it isn't needed.
739  if ((!eigenvalues_signal.empty() || !hessenberg_signal.empty() ||
740  !cond_signal.empty()) &&
741  dim > 0)
742  {
743  LAPACKFullMatrix<double> mat(dim, dim);
744  for (unsigned int i = 0; i < dim; ++i)
745  for (unsigned int j = 0; j < dim; ++j)
746  mat(i, j) = H_orig(i, j);
747  hessenberg_signal(H_orig);
748  // Avoid computing eigenvalues if they are not needed.
749  if (!eigenvalues_signal.empty())
750  {
751  // Copy mat so that we can compute svd below. Necessary since
752  // compute_eigenvalues will leave mat in state
753  // LAPACKSupport::unusable.
754  LAPACKFullMatrix<double> mat_eig(mat);
755  mat_eig.compute_eigenvalues();
756  std::vector<std::complex<double>> eigenvalues(dim);
757  for (unsigned int i = 0; i < mat_eig.n(); ++i)
758  eigenvalues[i] = mat_eig.eigenvalue(i);
759  // Sort eigenvalues for nicer output.
760  std::sort(eigenvalues.begin(),
761  eigenvalues.end(),
762  internal::SolverGMRESImplementation::complex_less_pred);
763  eigenvalues_signal(eigenvalues);
764  }
765  // Calculate condition number, avoid calculating the svd if a slot
766  // isn't connected. Need at least a 2-by-2 matrix to do the estimate.
767  if (!cond_signal.empty() && (mat.n() > 1))
768  {
769  mat.compute_svd();
770  double condition_number =
771  mat.singular_value(0) / mat.singular_value(mat.n() - 1);
772  cond_signal(condition_number);
773  }
774  }
775 }
776 
777 
778 
779 template <class VectorType>
780 template <typename MatrixType, typename PreconditionerType>
781 void
782 SolverGMRES<VectorType>::solve(const MatrixType & A,
783  VectorType & x,
784  const VectorType & b,
785  const PreconditionerType &preconditioner)
786 {
787  // TODO:[?] Check, why there are two different start residuals.
788  // TODO:[GK] Make sure the parameter in the constructor means maximum basis
789  // size
790 
791  LogStream::Prefix prefix("GMRES");
792 
793  // extra call to std::max to placate static analyzers: coverity rightfully
794  // complains that data.max_n_tmp_vectors - 2 may overflow
795  const unsigned int n_tmp_vectors =
796  std::max(additional_data.max_n_tmp_vectors, 3u);
797 
798  // Generate an object where basis vectors are stored.
800  n_tmp_vectors, this->memory);
801 
802  // number of the present iteration; this
803  // number is not reset to zero upon a
804  // restart
805  unsigned int accumulated_iterations = 0;
806 
807  const bool do_eigenvalues =
808  !condition_number_signal.empty() || !all_condition_numbers_signal.empty() ||
809  !eigenvalues_signal.empty() || !all_eigenvalues_signal.empty() ||
810  !hessenberg_signal.empty() || !all_hessenberg_signal.empty();
811  // for eigenvalue computation, need to collect the Hessenberg matrix (before
812  // applying Givens rotations)
813  FullMatrix<double> H_orig;
814  if (do_eigenvalues)
815  H_orig.reinit(n_tmp_vectors, n_tmp_vectors - 1);
816 
817  // matrix used for the orthogonalization process later
818  H.reinit(n_tmp_vectors, n_tmp_vectors - 1);
819 
820  // some additional vectors, also used in the orthogonalization
821  ::Vector<double> gamma(n_tmp_vectors), ci(n_tmp_vectors - 1),
822  si(n_tmp_vectors - 1), h(n_tmp_vectors - 1);
823 
824 
825  unsigned int dim = 0;
826 
828  double last_res = -std::numeric_limits<double>::max();
829 
830  // switch to determine whether we want a left or a right preconditioner. at
831  // present, left is default, but both ways are implemented
832  const bool left_precondition = !additional_data.right_preconditioning;
833 
834  // Per default the left preconditioned GMRes uses the preconditioned
835  // residual and the right preconditioned GMRes uses the unpreconditioned
836  // residual as stopping criterion.
837  const bool use_default_residual = additional_data.use_default_residual;
838 
839  // define two aliases
840  VectorType &v = tmp_vectors(0, x);
841  VectorType &p = tmp_vectors(n_tmp_vectors - 1, x);
842 
843  // Following vectors are needed when we are not using the default residuals
844  // as stopping criterion
847  std::unique_ptr<::Vector<double>> gamma_;
848  if (!use_default_residual)
849  {
850  r = std::move(typename VectorMemory<VectorType>::Pointer(this->memory));
851  x_ = std::move(typename VectorMemory<VectorType>::Pointer(this->memory));
852  r->reinit(x);
853  x_->reinit(x);
854 
855  gamma_ = std::make_unique<::Vector<double>>(gamma.size());
856  }
857 
858  bool re_orthogonalize = additional_data.force_re_orthogonalization;
859 
861  // outer iteration: loop until we either reach convergence or the maximum
862  // number of iterations is exceeded. each cycle of this loop amounts to one
863  // restart
864  do
865  {
866  // reset this vector to the right size
867  h.reinit(n_tmp_vectors - 1);
868 
869  if (left_precondition)
870  {
871  A.vmult(p, x);
872  p.sadd(-1., 1., b);
873  preconditioner.vmult(v, p);
874  }
875  else
876  {
877  A.vmult(v, x);
878  v.sadd(-1., 1., b);
879  };
880 
881  double rho = v.l2_norm();
882 
883  // check the residual here as well since it may be that we got the exact
884  // (or an almost exact) solution vector at the outset. if we wouldn't
885  // check here, the next scaling operation would produce garbage
886  if (use_default_residual)
887  {
888  last_res = rho;
889  iteration_state =
890  this->iteration_status(accumulated_iterations, rho, x);
891 
892  if (iteration_state != SolverControl::iterate)
893  break;
894  }
895  else
896  {
897  deallog << "default_res=" << rho << std::endl;
898 
899  if (left_precondition)
900  {
901  A.vmult(*r, x);
902  r->sadd(-1., 1., b);
903  }
904  else
905  preconditioner.vmult(*r, v);
906 
907  double res = r->l2_norm();
908  last_res = res;
909  iteration_state =
910  this->iteration_status(accumulated_iterations, res, x);
911 
912  if (iteration_state != SolverControl::iterate)
913  break;
914  }
915 
916  gamma(0) = rho;
917 
918  v *= 1. / rho;
919 
920  // inner iteration doing at most as many steps as there are temporary
921  // vectors. the number of steps actually been done is propagated outside
922  // through the @p dim variable
923  for (unsigned int inner_iteration = 0;
924  ((inner_iteration < n_tmp_vectors - 2) &&
925  (iteration_state == SolverControl::iterate));
926  ++inner_iteration)
927  {
928  ++accumulated_iterations;
929  // yet another alias
930  VectorType &vv = tmp_vectors(inner_iteration + 1, x);
931 
932  if (left_precondition)
933  {
934  A.vmult(p, tmp_vectors[inner_iteration]);
935  preconditioner.vmult(vv, p);
936  }
937  else
938  {
939  preconditioner.vmult(p, tmp_vectors[inner_iteration]);
940  A.vmult(vv, p);
941  }
942 
943  dim = inner_iteration + 1;
944 
945  const double s = modified_gram_schmidt(tmp_vectors,
946  dim,
947  accumulated_iterations,
948  vv,
949  h,
950  re_orthogonalize,
951  re_orthogonalize_signal);
952  h(inner_iteration + 1) = s;
953 
954  // s=0 is a lucky breakdown, the solver will reach convergence,
955  // but we must not divide by zero here.
956  if (s != 0)
957  vv *= 1. / s;
958 
959  // for eigenvalues, get the resulting coefficients from the
960  // orthogonalization process
961  if (do_eigenvalues)
962  for (unsigned int i = 0; i < dim + 1; ++i)
963  H_orig(i, inner_iteration) = h(i);
964 
965  // Transformation into tridiagonal structure
966  givens_rotation(h, gamma, ci, si, inner_iteration);
967 
968  // append vector on matrix
969  for (unsigned int i = 0; i < dim; ++i)
970  H(i, inner_iteration) = h(i);
971 
972  // default residual
973  rho = std::fabs(gamma(dim));
974 
975  if (use_default_residual)
976  {
977  last_res = rho;
978  iteration_state =
979  this->iteration_status(accumulated_iterations, rho, x);
980  }
981  else
982  {
983  deallog << "default_res=" << rho << std::endl;
984 
985  ::Vector<double> h_(dim);
986  *x_ = x;
987  *gamma_ = gamma;
988  H1.reinit(dim + 1, dim);
989 
990  for (unsigned int i = 0; i < dim + 1; ++i)
991  for (unsigned int j = 0; j < dim; ++j)
992  H1(i, j) = H(i, j);
993 
994  H1.backward(h_, *gamma_);
995 
996  if (left_precondition)
997  for (unsigned int i = 0; i < dim; ++i)
998  x_->add(h_(i), tmp_vectors[i]);
999  else
1000  {
1001  p = 0.;
1002  for (unsigned int i = 0; i < dim; ++i)
1003  p.add(h_(i), tmp_vectors[i]);
1004  preconditioner.vmult(*r, p);
1005  x_->add(1., *r);
1006  };
1007  A.vmult(*r, *x_);
1008  r->sadd(-1., 1., b);
1009  // Now *r contains the unpreconditioned residual!!
1010  if (left_precondition)
1011  {
1012  const double res = r->l2_norm();
1013  last_res = res;
1014 
1015  iteration_state =
1016  this->iteration_status(accumulated_iterations, res, x);
1017  }
1018  else
1019  {
1020  preconditioner.vmult(*x_, *r);
1021  const double preconditioned_res = x_->l2_norm();
1022  last_res = preconditioned_res;
1023 
1024  iteration_state =
1025  this->iteration_status(accumulated_iterations,
1026  preconditioned_res,
1027  x);
1028  }
1029  }
1030  };
1031  // end of inner iteration. now calculate the solution from the temporary
1032  // vectors
1033  h.reinit(dim);
1034  H1.reinit(dim + 1, dim);
1035 
1036  for (unsigned int i = 0; i < dim + 1; ++i)
1037  for (unsigned int j = 0; j < dim; ++j)
1038  H1(i, j) = H(i, j);
1039 
1040  compute_eigs_and_cond(H_orig,
1041  dim,
1042  all_eigenvalues_signal,
1043  all_hessenberg_signal,
1044  condition_number_signal);
1045 
1046  H1.backward(h, gamma);
1047 
1048  if (left_precondition)
1049  for (unsigned int i = 0; i < dim; ++i)
1050  x.add(h(i), tmp_vectors[i]);
1051  else
1052  {
1053  p = 0.;
1054  for (unsigned int i = 0; i < dim; ++i)
1055  p.add(h(i), tmp_vectors[i]);
1056  preconditioner.vmult(v, p);
1057  x.add(1., v);
1058  };
1059  // end of outer iteration. restart if no convergence and the number of
1060  // iterations is not exceeded
1061  }
1062  while (iteration_state == SolverControl::iterate);
1063 
1064  compute_eigs_and_cond(H_orig,
1065  dim,
1066  eigenvalues_signal,
1067  hessenberg_signal,
1068  condition_number_signal);
1069 
1070  if (!krylov_space_signal.empty())
1071  krylov_space_signal(tmp_vectors);
1072 
1073  // in case of failure: throw exception
1074  AssertThrow(iteration_state == SolverControl::success,
1075  SolverControl::NoConvergence(accumulated_iterations, last_res));
1076 }
1077 
1078 
1079 
1080 template <class VectorType>
1081 boost::signals2::connection
1083  const std::function<void(double)> &slot,
1084  const bool every_iteration)
1085 {
1086  if (every_iteration)
1087  {
1088  return all_condition_numbers_signal.connect(slot);
1089  }
1090  else
1091  {
1092  return condition_number_signal.connect(slot);
1093  }
1094 }
1095 
1096 
1097 
1098 template <class VectorType>
1099 boost::signals2::connection
1101  const std::function<void(const std::vector<std::complex<double>> &)> &slot,
1102  const bool every_iteration)
1103 {
1104  if (every_iteration)
1105  {
1106  return all_eigenvalues_signal.connect(slot);
1107  }
1108  else
1109  {
1110  return eigenvalues_signal.connect(slot);
1111  }
1112 }
1113 
1114 
1115 
1116 template <class VectorType>
1117 boost::signals2::connection
1119  const std::function<void(const FullMatrix<double> &)> &slot,
1120  const bool every_iteration)
1121 {
1122  if (every_iteration)
1123  {
1124  return all_hessenberg_signal.connect(slot);
1125  }
1126  else
1127  {
1128  return hessenberg_signal.connect(slot);
1129  }
1130 }
1131 
1132 
1133 
1134 template <class VectorType>
1135 boost::signals2::connection
1137  const std::function<void(
1139 {
1140  return krylov_space_signal.connect(slot);
1141 }
1142 
1143 
1144 
1145 template <class VectorType>
1146 boost::signals2::connection
1148  const std::function<void(int)> &slot)
1149 {
1150  return re_orthogonalize_signal.connect(slot);
1151 }
1152 
1153 
1154 
1155 template <class VectorType>
1156 double
1158 {
1159  // dummy implementation. this function is not needed for the present
1160  // implementation of gmres
1161  Assert(false, ExcInternalError());
1162  return 0;
1163 }
1164 
1165 
1166 //----------------------------------------------------------------------//
1167 
1168 template <class VectorType>
1171  const AdditionalData & data)
1173  , additional_data(data)
1174 {}
1175 
1176 
1177 
1178 template <class VectorType>
1180  const AdditionalData &data)
1182  , additional_data(data)
1183 {}
1184 
1185 
1186 
1187 template <class VectorType>
1188 template <typename MatrixType, typename PreconditionerType>
1189 void
1190 SolverFGMRES<VectorType>::solve(const MatrixType & A,
1191  VectorType & x,
1192  const VectorType & b,
1193  const PreconditionerType &preconditioner)
1194 {
1195  LogStream::Prefix prefix("FGMRES");
1196 
1197  SolverControl::State iteration_state = SolverControl::iterate;
1198 
1199  const unsigned int basis_size = additional_data.max_basis_size;
1200 
1201  // Generate an object where basis vectors are stored.
1203  basis_size, this->memory);
1205  basis_size, this->memory);
1206 
1207  // number of the present iteration; this number is not reset to zero upon a
1208  // restart
1209  unsigned int accumulated_iterations = 0;
1210 
1211  // matrix used for the orthogonalization process later
1212  H.reinit(basis_size + 1, basis_size);
1213 
1214  // Vectors for projected system
1215  Vector<double> projected_rhs;
1216  Vector<double> y;
1217 
1218  // Iteration starts here
1219  double res = -std::numeric_limits<double>::max();
1220 
1221  typename VectorMemory<VectorType>::Pointer aux(this->memory);
1222  aux->reinit(x);
1223  do
1224  {
1225  A.vmult(*aux, x);
1226  aux->sadd(-1., 1., b);
1227 
1228  double beta = aux->l2_norm();
1229  res = beta;
1230  iteration_state = this->iteration_status(accumulated_iterations, res, x);
1231  if (iteration_state == SolverControl::success)
1232  break;
1233 
1234  H.reinit(basis_size + 1, basis_size);
1235  double a = beta;
1236 
1237  for (unsigned int j = 0; j < basis_size; ++j)
1238  {
1239  if (a != 0) // treat lucky breakdown
1240  v(j, x).equ(1. / a, *aux);
1241  else
1242  v(j, x) = 0.;
1243 
1244 
1245  preconditioner.vmult(z(j, x), v[j]);
1246  A.vmult(*aux, z[j]);
1247 
1248  // Gram-Schmidt
1249  H(0, j) = *aux * v[0];
1250  for (unsigned int i = 1; i <= j; ++i)
1251  H(i, j) = aux->add_and_dot(-H(i - 1, j), v[i - 1], v[i]);
1252  H(j + 1, j) = a = std::sqrt(aux->add_and_dot(-H(j, j), v[j], *aux));
1253 
1254  // Compute projected solution
1255 
1256  if (j > 0)
1257  {
1258  H1.reinit(j + 1, j);
1259  projected_rhs.reinit(j + 1);
1260  y.reinit(j);
1261  projected_rhs(0) = beta;
1262  H1.fill(H);
1263 
1264  // check convergence. note that the vector 'x' we pass to the
1265  // criterion is not the final solution we compute if we
1266  // decide to jump out of the iteration (we update 'x' again
1267  // right after the current loop)
1268  Householder<double> house(H1);
1269  res = house.least_squares(y, projected_rhs);
1270  iteration_state =
1271  this->iteration_status(++accumulated_iterations, res, x);
1272  if (iteration_state != SolverControl::iterate)
1273  break;
1274  }
1275  }
1276 
1277  // Update solution vector
1278  for (unsigned int j = 0; j < y.size(); ++j)
1279  x.add(y(j), z[j]);
1280  }
1281  while (iteration_state == SolverControl::iterate);
1282 
1283  // in case of failure: throw exception
1284  if (iteration_state != SolverControl::success)
1285  AssertThrow(false,
1286  SolverControl::NoConvergence(accumulated_iterations, res));
1287 }
1288 
1289 #endif // DOXYGEN
1290 
1292 
1293 #endif
PETScWrappers::SolverGMRES SolverGMRES
SolverFGMRES(SolverControl &cn, VectorMemory< VectorType > &mem, const AdditionalData &data=AdditionalData())
LogStream deallog
Definition: logstream.cc:37
TmpVectors(const unsigned int max_size, VectorMemory< VectorType > &vmem)
boost::signals2::connection connect_condition_number_slot(const std::function< void(double)> &slot, const bool every_iteration=false)
Continue iteration.
size_type n() const
boost::signals2::signal< void(const internal::SolverGMRESImplementation::TmpVectors< VectorType > &)> krylov_space_signal
Definition: solver_gmres.h:364
std::complex< number > eigenvalue(const size_type i) const
void solve(const MatrixType &A, VectorType &x, const VectorType &b, const PreconditionerType &preconditioner)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1691
boost::signals2::connection connect_re_orthogonalization_slot(const std::function< void(int)> &slot)
FullMatrix< double > H1
Definition: solver_gmres.h:518
boost::signals2::signal< void(double)> condition_number_signal
Definition: solver_gmres.h:323
AdditionalData additional_data
Definition: solver_gmres.h:317
static ::ExceptionBase & ExcNotInitialized()
std::vector< typename VectorMemory< VectorType >::Pointer > data
Definition: solver_gmres.h:105
#define AssertThrow(cond, exc)
Definition: exceptions.h:1576
AdditionalData additional_data
Definition: solver_gmres.h:508
FullMatrix< double > H1
Definition: solver_gmres.h:435
VectorType & operator()(const unsigned int i, const VectorType &temp)
SymmetricTensor< 2, dim, Number > epsilon(const Tensor< 2, dim, Number > &Grad_u)
AdditionalData(const unsigned int max_basis_size=30)
Definition: solver_gmres.h:470
boost::signals2::signal< void(int)> re_orthogonalize_signal
Definition: solver_gmres.h:370
FullMatrix< double > H
Definition: solver_gmres.h:430
void solve(const MatrixType &A, VectorType &x, const VectorType &b, const PreconditionerType &preconditioner)
static ::ExceptionBase & ExcMessage(std::string arg1)
AdditionalData(const unsigned int max_n_tmp_vectors=30, const bool right_preconditioning=false, const bool use_default_residual=true, const bool force_re_orthogonalization=false)
#define DeclException1(Exception1, type1, outsequence)
Definition: exceptions.h:515
Stop iteration, goal reached.
void reinit(const TableIndices< N > &new_size, const bool omit_default_initialization=false)
#define Assert(cond, exc)
Definition: exceptions.h:1466
boost::signals2::connection connect_krylov_space_slot(const std::function< void(const internal::SolverGMRESImplementation::TmpVectors< VectorType > &)> &slot)
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:380
static double modified_gram_schmidt(const internal::SolverGMRESImplementation::TmpVectors< VectorType > &orthogonal_vectors, const unsigned int dim, const unsigned int accumulated_iterations, VectorType &vv, Vector< double > &h, bool &re_orthogonalize, const boost::signals2::signal< void(int)> &re_orthogonalize_signal=boost::signals2::signal< void(int)>())
Expression fabs(const Expression &x)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
number singular_value(const size_type i) const
static const char A
virtual double criterion()
std::array< NumberType, 3 > givens_rotation(const NumberType &x, const NumberType &y)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:379
double least_squares(Vector< number2 > &dst, const Vector< number2 > &src) const
void compute_eigenvalues(const bool right_eigenvectors=false, const bool left_eigenvectors=false)
long double gamma(const unsigned int n)
Eigenvalue vector is filled.
VectorType & operator[](const unsigned int i) const
boost::signals2::connection connect_hessenberg_slot(const std::function< void(const FullMatrix< double > &)> &slot, const bool every_iteration=true)
FullMatrix< double > H
Definition: solver_gmres.h:513
SolverGMRES(SolverControl &cn, VectorMemory< VectorType > &mem, const AdditionalData &data=AdditionalData())
void givens_rotation(Vector< double > &h, Vector< double > &b, Vector< double > &ci, Vector< double > &si, int col) const
boost::signals2::connection connect_eigenvalues_slot(const std::function< void(const std::vector< std::complex< double >> &)> &slot, const bool every_iteration=false)
static void compute_eigs_and_cond(const FullMatrix< double > &H_orig, const unsigned int dim, const boost::signals2::signal< void(const std::vector< std::complex< double >> &)> &eigenvalues_signal, const boost::signals2::signal< void(const FullMatrix< double > &)> &hessenberg_signal, const boost::signals2::signal< void(double)> &cond_signal)
boost::signals2::signal< void(double)> all_condition_numbers_signal
Definition: solver_gmres.h:329
T max(const T &t, const MPI_Comm &mpi_communicator)
static ::ExceptionBase & ExcInternalError()