Reference documentation for deal.II version Git 423cd11810 2020-05-28 01:07:25 +0200
\(\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 
177 template <class VectorType = Vector<double>>
178 class SolverGMRES : public SolverBase<VectorType>
179 {
180 public:
185  {
192  explicit AdditionalData(const unsigned int max_n_tmp_vectors = 30,
193  const bool right_preconditioning = false,
194  const bool use_default_residual = true,
195  const bool force_re_orthogonalization = false);
196 
203  unsigned int max_n_tmp_vectors;
204 
213 
218 
226  };
227 
233  const AdditionalData & data = AdditionalData());
234 
240 
244  SolverGMRES(const SolverGMRES<VectorType> &) = delete;
245 
249  template <typename MatrixType, typename PreconditionerType>
250  void
251  solve(const MatrixType & A,
252  VectorType & x,
253  const VectorType & b,
254  const PreconditionerType &preconditioner);
255 
262  boost::signals2::connection
263  connect_condition_number_slot(const std::function<void(double)> &slot,
264  const bool every_iteration = false);
265 
272  boost::signals2::connection
273  connect_eigenvalues_slot(
274  const std::function<void(const std::vector<std::complex<double>> &)> &slot,
275  const bool every_iteration = false);
276 
284  boost::signals2::connection
285  connect_hessenberg_slot(
286  const std::function<void(const FullMatrix<double> &)> &slot,
287  const bool every_iteration = true);
288 
295  boost::signals2::connection
296  connect_krylov_space_slot(
297  const std::function<
299  &slot);
300 
301 
306  boost::signals2::connection
307  connect_re_orthogonalization_slot(const std::function<void(int)> &slot);
308 
309 
310  DeclException1(ExcTooFewTmpVectors,
311  int,
312  << "The number of temporary vectors you gave (" << arg1
313  << ") is too small. It should be at least 10 for "
314  << "any results, and much more for reasonable ones.");
315 
316 protected:
321 
326  boost::signals2::signal<void(double)> condition_number_signal;
327 
332  boost::signals2::signal<void(double)> all_condition_numbers_signal;
333 
338  boost::signals2::signal<void(const std::vector<std::complex<double>> &)>
339  eigenvalues_signal;
340 
345  boost::signals2::signal<void(const std::vector<std::complex<double>> &)>
346  all_eigenvalues_signal;
347 
352  boost::signals2::signal<void(const FullMatrix<double> &)> hessenberg_signal;
353 
358  boost::signals2::signal<void(const FullMatrix<double> &)>
359  all_hessenberg_signal;
360 
365  boost::signals2::signal<void(
368 
373  boost::signals2::signal<void(int)> re_orthogonalize_signal;
374 
378  virtual double
379  criterion();
380 
385  void
387  Vector<double> &b,
388  Vector<double> &ci,
389  Vector<double> &si,
390  int col) const;
391 
402  static double
403  modified_gram_schmidt(
405  & orthogonal_vectors,
406  const unsigned int dim,
407  const unsigned int accumulated_iterations,
408  VectorType & vv,
409  Vector<double> & h,
410  bool & re_orthogonalize,
411  const boost::signals2::signal<void(int)> &re_orthogonalize_signal =
412  boost::signals2::signal<void(int)>());
413 
420  static void
421  compute_eigs_and_cond(
422  const FullMatrix<double> &H_orig,
423  const unsigned int dim,
424  const boost::signals2::signal<
425  void(const std::vector<std::complex<double>> &)> &eigenvalues_signal,
426  const boost::signals2::signal<void(const FullMatrix<double> &)>
427  & hessenberg_signal,
428  const boost::signals2::signal<void(double)> &cond_signal);
429 
434 
439 };
440 
463 template <class VectorType = Vector<double>>
464 class SolverFGMRES : public SolverBase<VectorType>
465 {
466 public:
471  {
475  explicit AdditionalData(const unsigned int max_basis_size = 30)
476  : max_basis_size(max_basis_size)
477  {}
478 
484  AdditionalData(const unsigned int max_basis_size,
485  const bool use_default_residual)
486  : max_basis_size(max_basis_size)
487  {
488  (void)use_default_residual;
489  }
490 
494  unsigned int max_basis_size;
495  };
496 
502  const AdditionalData & data = AdditionalData());
503 
509  const AdditionalData &data = AdditionalData());
510 
514  template <typename MatrixType, typename PreconditionerType>
515  void
516  solve(const MatrixType & A,
517  VectorType & x,
518  const VectorType & b,
519  const PreconditionerType &preconditioner);
520 
521 private:
526 
531 
536 };
537 
539 /* --------------------- Inline and template functions ------------------- */
540 
541 
542 #ifndef DOXYGEN
543 namespace internal
544 {
545  namespace SolverGMRESImplementation
546  {
547  template <class VectorType>
548  inline TmpVectors<VectorType>::TmpVectors(const unsigned int max_size,
550  : mem(vmem)
551  , data(max_size)
552  {}
553 
554 
555 
556  template <class VectorType>
558  operator[](const unsigned int i) const
559  {
560  AssertIndexRange(i, data.size());
561 
562  Assert(data[i] != nullptr, ExcNotInitialized());
563  return *data[i];
564  }
565 
566 
567 
568  template <class VectorType>
569  inline VectorType &
570  TmpVectors<VectorType>::operator()(const unsigned int i,
571  const VectorType & temp)
572  {
573  AssertIndexRange(i, data.size());
574  if (data[i] == nullptr)
575  {
576  data[i] = std::move(typename VectorMemory<VectorType>::Pointer(mem));
577  data[i]->reinit(temp);
578  }
579  return *data[i];
580  }
581 
582 
583 
584  template <class VectorType>
585  unsigned int
587  {
588  return (data.size() > 0 ? data.size() - 1 : 0);
589  }
590 
591 
592 
593  // A comparator for better printing eigenvalues
594  inline bool
595  complex_less_pred(const std::complex<double> &x,
596  const std::complex<double> &y)
597  {
598  return x.real() < y.real() ||
599  (x.real() == y.real() && x.imag() < y.imag());
600  }
601  } // namespace SolverGMRESImplementation
602 } // namespace internal
603 
604 
605 
606 template <class VectorType>
608  const unsigned int max_n_tmp_vectors,
609  const bool right_preconditioning,
610  const bool use_default_residual,
611  const bool force_re_orthogonalization)
612  : max_n_tmp_vectors(max_n_tmp_vectors)
613  , right_preconditioning(right_preconditioning)
614  , use_default_residual(use_default_residual)
615  , force_re_orthogonalization(force_re_orthogonalization)
616 {
617  Assert(3 <= max_n_tmp_vectors,
618  ExcMessage("SolverGMRES needs at least three "
619  "temporary vectors."));
620 }
621 
622 
623 
624 template <class VectorType>
627  const AdditionalData & data)
629  , additional_data(data)
630 {}
631 
632 
633 
634 template <class VectorType>
636  const AdditionalData &data)
638  , additional_data(data)
639 {}
640 
641 
642 
643 template <class VectorType>
644 inline void
646  Vector<double> &b,
647  Vector<double> &ci,
648  Vector<double> &si,
649  int col) const
650 {
651  for (int i = 0; i < col; i++)
652  {
653  const double s = si(i);
654  const double c = ci(i);
655  const double dummy = h(i);
656  h(i) = c * dummy + s * h(i + 1);
657  h(i + 1) = -s * dummy + c * h(i + 1);
658  };
659 
660  const double r = 1. / std::sqrt(h(col) * h(col) + h(col + 1) * h(col + 1));
661  si(col) = h(col + 1) * r;
662  ci(col) = h(col) * r;
663  h(col) = ci(col) * h(col) + si(col) * h(col + 1);
664  b(col + 1) = -si(col) * b(col);
665  b(col) *= ci(col);
666 }
667 
668 
669 
670 template <class VectorType>
671 inline double
674  & orthogonal_vectors,
675  const unsigned int dim,
676  const unsigned int accumulated_iterations,
677  VectorType & vv,
678  Vector<double> & h,
679  bool & reorthogonalize,
680  const boost::signals2::signal<void(int)> &reorthogonalize_signal)
681 {
682  Assert(dim > 0, ExcInternalError());
683  const unsigned int inner_iteration = dim - 1;
684 
685  // need initial norm for detection of re-orthogonalization, see below
686  double norm_vv_start = 0;
687  const bool consider_reorthogonalize =
688  (reorthogonalize == false) && (inner_iteration % 5 == 4);
689  if (consider_reorthogonalize)
690  norm_vv_start = vv.l2_norm();
691 
692  // Orthogonalization
693  h(0) = vv * orthogonal_vectors[0];
694  for (unsigned int i = 1; i < dim; ++i)
695  h(i) = vv.add_and_dot(-h(i - 1),
696  orthogonal_vectors[i - 1],
697  orthogonal_vectors[i]);
698  double norm_vv =
699  std::sqrt(vv.add_and_dot(-h(dim - 1), orthogonal_vectors[dim - 1], vv));
700 
701  // Re-orthogonalization if loss of orthogonality detected. For the test, use
702  // a strategy discussed in C. T. Kelley, Iterative Methods for Linear and
703  // Nonlinear Equations, SIAM, Philadelphia, 1995: Compare the norm of vv
704  // after orthogonalization with its norm when starting the
705  // orthogonalization. If vv became very small (here: less than the square
706  // root of the machine precision times 10), it is almost in the span of the
707  // previous vectors, which indicates loss of precision.
708  if (consider_reorthogonalize)
709  {
710  if (norm_vv >
711  10. * norm_vv_start *
712  std::sqrt(
714  return norm_vv;
715 
716  else
717  {
718  reorthogonalize = true;
719  if (!reorthogonalize_signal.empty())
720  reorthogonalize_signal(accumulated_iterations);
721  }
722  }
723 
724  if (reorthogonalize == true)
725  {
726  double htmp = vv * orthogonal_vectors[0];
727  h(0) += htmp;
728  for (unsigned int i = 1; i < dim; ++i)
729  {
730  htmp = vv.add_and_dot(-htmp,
731  orthogonal_vectors[i - 1],
732  orthogonal_vectors[i]);
733  h(i) += htmp;
734  }
735  norm_vv =
736  std::sqrt(vv.add_and_dot(-htmp, orthogonal_vectors[dim - 1], vv));
737  }
738 
739  return norm_vv;
740 }
741 
742 
743 
744 template <class VectorType>
745 inline void
747  const FullMatrix<double> &H_orig,
748  const unsigned int dim,
749  const boost::signals2::signal<void(const std::vector<std::complex<double>> &)>
750  &eigenvalues_signal,
751  const boost::signals2::signal<void(const FullMatrix<double> &)>
752  & hessenberg_signal,
753  const boost::signals2::signal<void(double)> &cond_signal)
754 {
755  // Avoid copying the Hessenberg matrix if it isn't needed.
756  if ((!eigenvalues_signal.empty() || !hessenberg_signal.empty() ||
757  !cond_signal.empty()) &&
758  dim > 0)
759  {
760  LAPACKFullMatrix<double> mat(dim, dim);
761  for (unsigned int i = 0; i < dim; ++i)
762  for (unsigned int j = 0; j < dim; ++j)
763  mat(i, j) = H_orig(i, j);
764  hessenberg_signal(H_orig);
765  // Avoid computing eigenvalues if they are not needed.
766  if (!eigenvalues_signal.empty())
767  {
768  // Copy mat so that we can compute svd below. Necessary since
769  // compute_eigenvalues will leave mat in state
770  // LAPACKSupport::unusable.
771  LAPACKFullMatrix<double> mat_eig(mat);
772  mat_eig.compute_eigenvalues();
773  std::vector<std::complex<double>> eigenvalues(dim);
774  for (unsigned int i = 0; i < mat_eig.n(); ++i)
775  eigenvalues[i] = mat_eig.eigenvalue(i);
776  // Sort eigenvalues for nicer output.
777  std::sort(eigenvalues.begin(),
778  eigenvalues.end(),
779  internal::SolverGMRESImplementation::complex_less_pred);
780  eigenvalues_signal(eigenvalues);
781  }
782  // Calculate condition number, avoid calculating the svd if a slot
783  // isn't connected. Need at least a 2-by-2 matrix to do the estimate.
784  if (!cond_signal.empty() && (mat.n() > 1))
785  {
786  mat.compute_svd();
787  double condition_number =
788  mat.singular_value(0) / mat.singular_value(mat.n() - 1);
789  cond_signal(condition_number);
790  }
791  }
792 }
793 
794 
795 
796 template <class VectorType>
797 template <typename MatrixType, typename PreconditionerType>
798 void
799 SolverGMRES<VectorType>::solve(const MatrixType & A,
800  VectorType & x,
801  const VectorType & b,
802  const PreconditionerType &preconditioner)
803 {
804  // TODO:[?] Check, why there are two different start residuals.
805  // TODO:[GK] Make sure the parameter in the constructor means maximum basis
806  // size
807 
808  LogStream::Prefix prefix("GMRES");
809 
810  // extra call to std::max to placate static analyzers: coverity rightfully
811  // complains that data.max_n_tmp_vectors - 2 may overflow
812  const unsigned int n_tmp_vectors =
813  std::max(additional_data.max_n_tmp_vectors, 3u);
814 
815  // Generate an object where basis vectors are stored.
817  n_tmp_vectors, this->memory);
818 
819  // number of the present iteration; this
820  // number is not reset to zero upon a
821  // restart
822  unsigned int accumulated_iterations = 0;
823 
824  const bool do_eigenvalues =
825  !condition_number_signal.empty() || !all_condition_numbers_signal.empty() ||
826  !eigenvalues_signal.empty() || !all_eigenvalues_signal.empty() ||
827  !hessenberg_signal.empty() || !all_hessenberg_signal.empty();
828  // for eigenvalue computation, need to collect the Hessenberg matrix (before
829  // applying Givens rotations)
830  FullMatrix<double> H_orig;
831  if (do_eigenvalues)
832  H_orig.reinit(n_tmp_vectors, n_tmp_vectors - 1);
833 
834  // matrix used for the orthogonalization process later
835  H.reinit(n_tmp_vectors, n_tmp_vectors - 1);
836 
837  // some additional vectors, also used in the orthogonalization
838  ::Vector<double> gamma(n_tmp_vectors), ci(n_tmp_vectors - 1),
839  si(n_tmp_vectors - 1), h(n_tmp_vectors - 1);
840 
841 
842  unsigned int dim = 0;
843 
845  double last_res = -std::numeric_limits<double>::max();
846 
847  // switch to determine whether we want a left or a right preconditioner. at
848  // present, left is default, but both ways are implemented
849  const bool left_precondition = !additional_data.right_preconditioning;
850 
851  // Per default the left preconditioned GMRes uses the preconditioned
852  // residual and the right preconditioned GMRes uses the unpreconditioned
853  // residual as stopping criterion.
854  const bool use_default_residual = additional_data.use_default_residual;
855 
856  // define two aliases
857  VectorType &v = tmp_vectors(0, x);
858  VectorType &p = tmp_vectors(n_tmp_vectors - 1, x);
859 
860  // Following vectors are needed when we are not using the default residuals
861  // as stopping criterion
864  std::unique_ptr<::Vector<double>> gamma_;
865  if (!use_default_residual)
866  {
867  r = std::move(typename VectorMemory<VectorType>::Pointer(this->memory));
868  x_ = std::move(typename VectorMemory<VectorType>::Pointer(this->memory));
869  r->reinit(x);
870  x_->reinit(x);
871 
872  gamma_ = std::make_unique<::Vector<double>>(gamma.size());
873  }
874 
875  bool re_orthogonalize = additional_data.force_re_orthogonalization;
876 
878  // outer iteration: loop until we either reach convergence or the maximum
879  // number of iterations is exceeded. each cycle of this loop amounts to one
880  // restart
881  do
882  {
883  // reset this vector to the right size
884  h.reinit(n_tmp_vectors - 1);
885 
886  if (left_precondition)
887  {
888  A.vmult(p, x);
889  p.sadd(-1., 1., b);
890  preconditioner.vmult(v, p);
891  }
892  else
893  {
894  A.vmult(v, x);
895  v.sadd(-1., 1., b);
896  };
897 
898  double rho = v.l2_norm();
899 
900  // check the residual here as well since it may be that we got the exact
901  // (or an almost exact) solution vector at the outset. if we wouldn't
902  // check here, the next scaling operation would produce garbage
903  if (use_default_residual)
904  {
905  last_res = rho;
906  iteration_state =
907  this->iteration_status(accumulated_iterations, rho, x);
908 
909  if (iteration_state != SolverControl::iterate)
910  break;
911  }
912  else
913  {
914  deallog << "default_res=" << rho << std::endl;
915 
916  if (left_precondition)
917  {
918  A.vmult(*r, x);
919  r->sadd(-1., 1., b);
920  }
921  else
922  preconditioner.vmult(*r, v);
923 
924  double res = r->l2_norm();
925  last_res = res;
926  iteration_state =
927  this->iteration_status(accumulated_iterations, res, x);
928 
929  if (iteration_state != SolverControl::iterate)
930  break;
931  }
932 
933  gamma(0) = rho;
934 
935  v *= 1. / rho;
936 
937  // inner iteration doing at most as many steps as there are temporary
938  // vectors. the number of steps actually been done is propagated outside
939  // through the @p dim variable
940  for (unsigned int inner_iteration = 0;
941  ((inner_iteration < n_tmp_vectors - 2) &&
942  (iteration_state == SolverControl::iterate));
943  ++inner_iteration)
944  {
945  ++accumulated_iterations;
946  // yet another alias
947  VectorType &vv = tmp_vectors(inner_iteration + 1, x);
948 
949  if (left_precondition)
950  {
951  A.vmult(p, tmp_vectors[inner_iteration]);
952  preconditioner.vmult(vv, p);
953  }
954  else
955  {
956  preconditioner.vmult(p, tmp_vectors[inner_iteration]);
957  A.vmult(vv, p);
958  }
959 
960  dim = inner_iteration + 1;
961 
962  const double s = modified_gram_schmidt(tmp_vectors,
963  dim,
964  accumulated_iterations,
965  vv,
966  h,
967  re_orthogonalize,
968  re_orthogonalize_signal);
969  h(inner_iteration + 1) = s;
970 
971  // s=0 is a lucky breakdown, the solver will reach convergence,
972  // but we must not divide by zero here.
973  if (s != 0)
974  vv *= 1. / s;
975 
976  // for eigenvalues, get the resulting coefficients from the
977  // orthogonalization process
978  if (do_eigenvalues)
979  for (unsigned int i = 0; i < dim + 1; ++i)
980  H_orig(i, inner_iteration) = h(i);
981 
982  // Transformation into tridiagonal structure
983  givens_rotation(h, gamma, ci, si, inner_iteration);
984 
985  // append vector on matrix
986  for (unsigned int i = 0; i < dim; ++i)
987  H(i, inner_iteration) = h(i);
988 
989  // default residual
990  rho = std::fabs(gamma(dim));
991 
992  if (use_default_residual)
993  {
994  last_res = rho;
995  iteration_state =
996  this->iteration_status(accumulated_iterations, rho, x);
997  }
998  else
999  {
1000  deallog << "default_res=" << rho << std::endl;
1001 
1002  ::Vector<double> h_(dim);
1003  *x_ = x;
1004  *gamma_ = gamma;
1005  H1.reinit(dim + 1, dim);
1006 
1007  for (unsigned int i = 0; i < dim + 1; ++i)
1008  for (unsigned int j = 0; j < dim; ++j)
1009  H1(i, j) = H(i, j);
1010 
1011  H1.backward(h_, *gamma_);
1012 
1013  if (left_precondition)
1014  for (unsigned int i = 0; i < dim; ++i)
1015  x_->add(h_(i), tmp_vectors[i]);
1016  else
1017  {
1018  p = 0.;
1019  for (unsigned int i = 0; i < dim; ++i)
1020  p.add(h_(i), tmp_vectors[i]);
1021  preconditioner.vmult(*r, p);
1022  x_->add(1., *r);
1023  };
1024  A.vmult(*r, *x_);
1025  r->sadd(-1., 1., b);
1026  // Now *r contains the unpreconditioned residual!!
1027  if (left_precondition)
1028  {
1029  const double res = r->l2_norm();
1030  last_res = res;
1031 
1032  iteration_state =
1033  this->iteration_status(accumulated_iterations, res, x);
1034  }
1035  else
1036  {
1037  preconditioner.vmult(*x_, *r);
1038  const double preconditioned_res = x_->l2_norm();
1039  last_res = preconditioned_res;
1040 
1041  iteration_state =
1042  this->iteration_status(accumulated_iterations,
1043  preconditioned_res,
1044  x);
1045  }
1046  }
1047  };
1048  // end of inner iteration. now calculate the solution from the temporary
1049  // vectors
1050  h.reinit(dim);
1051  H1.reinit(dim + 1, dim);
1052 
1053  for (unsigned int i = 0; i < dim + 1; ++i)
1054  for (unsigned int j = 0; j < dim; ++j)
1055  H1(i, j) = H(i, j);
1056 
1057  compute_eigs_and_cond(H_orig,
1058  dim,
1059  all_eigenvalues_signal,
1060  all_hessenberg_signal,
1061  condition_number_signal);
1062 
1063  H1.backward(h, gamma);
1064 
1065  if (left_precondition)
1066  for (unsigned int i = 0; i < dim; ++i)
1067  x.add(h(i), tmp_vectors[i]);
1068  else
1069  {
1070  p = 0.;
1071  for (unsigned int i = 0; i < dim; ++i)
1072  p.add(h(i), tmp_vectors[i]);
1073  preconditioner.vmult(v, p);
1074  x.add(1., v);
1075  };
1076  // end of outer iteration. restart if no convergence and the number of
1077  // iterations is not exceeded
1078  }
1079  while (iteration_state == SolverControl::iterate);
1080 
1081  compute_eigs_and_cond(H_orig,
1082  dim,
1083  eigenvalues_signal,
1084  hessenberg_signal,
1085  condition_number_signal);
1086 
1087  if (!krylov_space_signal.empty())
1088  krylov_space_signal(tmp_vectors);
1089 
1090  // in case of failure: throw exception
1091  AssertThrow(iteration_state == SolverControl::success,
1092  SolverControl::NoConvergence(accumulated_iterations, last_res));
1093 }
1094 
1095 
1096 
1097 template <class VectorType>
1098 boost::signals2::connection
1100  const std::function<void(double)> &slot,
1101  const bool every_iteration)
1102 {
1103  if (every_iteration)
1104  {
1105  return all_condition_numbers_signal.connect(slot);
1106  }
1107  else
1108  {
1109  return condition_number_signal.connect(slot);
1110  }
1111 }
1112 
1113 
1114 
1115 template <class VectorType>
1116 boost::signals2::connection
1118  const std::function<void(const std::vector<std::complex<double>> &)> &slot,
1119  const bool every_iteration)
1120 {
1121  if (every_iteration)
1122  {
1123  return all_eigenvalues_signal.connect(slot);
1124  }
1125  else
1126  {
1127  return eigenvalues_signal.connect(slot);
1128  }
1129 }
1130 
1131 
1132 
1133 template <class VectorType>
1134 boost::signals2::connection
1136  const std::function<void(const FullMatrix<double> &)> &slot,
1137  const bool every_iteration)
1138 {
1139  if (every_iteration)
1140  {
1141  return all_hessenberg_signal.connect(slot);
1142  }
1143  else
1144  {
1145  return hessenberg_signal.connect(slot);
1146  }
1147 }
1148 
1149 
1150 
1151 template <class VectorType>
1152 boost::signals2::connection
1154  const std::function<void(
1156 {
1157  return krylov_space_signal.connect(slot);
1158 }
1159 
1160 
1161 
1162 template <class VectorType>
1163 boost::signals2::connection
1165  const std::function<void(int)> &slot)
1166 {
1167  return re_orthogonalize_signal.connect(slot);
1168 }
1169 
1170 
1171 
1172 template <class VectorType>
1173 double
1175 {
1176  // dummy implementation. this function is not needed for the present
1177  // implementation of gmres
1178  Assert(false, ExcInternalError());
1179  return 0;
1180 }
1181 
1182 
1183 //----------------------------------------------------------------------//
1184 
1185 template <class VectorType>
1188  const AdditionalData & data)
1190  , additional_data(data)
1191 {}
1192 
1193 
1194 
1195 template <class VectorType>
1197  const AdditionalData &data)
1199  , additional_data(data)
1200 {}
1201 
1202 
1203 
1204 template <class VectorType>
1205 template <typename MatrixType, typename PreconditionerType>
1206 void
1207 SolverFGMRES<VectorType>::solve(const MatrixType & A,
1208  VectorType & x,
1209  const VectorType & b,
1210  const PreconditionerType &preconditioner)
1211 {
1212  LogStream::Prefix prefix("FGMRES");
1213 
1214  SolverControl::State iteration_state = SolverControl::iterate;
1215 
1216  const unsigned int basis_size = additional_data.max_basis_size;
1217 
1218  // Generate an object where basis vectors are stored.
1220  basis_size, this->memory);
1222  basis_size, this->memory);
1223 
1224  // number of the present iteration; this number is not reset to zero upon a
1225  // restart
1226  unsigned int accumulated_iterations = 0;
1227 
1228  // matrix used for the orthogonalization process later
1229  H.reinit(basis_size + 1, basis_size);
1230 
1231  // Vectors for projected system
1232  Vector<double> projected_rhs;
1233  Vector<double> y;
1234 
1235  // Iteration starts here
1236  double res = -std::numeric_limits<double>::max();
1237 
1238  typename VectorMemory<VectorType>::Pointer aux(this->memory);
1239  aux->reinit(x);
1240  do
1241  {
1242  A.vmult(*aux, x);
1243  aux->sadd(-1., 1., b);
1244 
1245  double beta = aux->l2_norm();
1246  res = beta;
1247  iteration_state = this->iteration_status(accumulated_iterations, res, x);
1248  if (iteration_state == SolverControl::success)
1249  break;
1250 
1251  H.reinit(basis_size + 1, basis_size);
1252  double a = beta;
1253 
1254  for (unsigned int j = 0; j < basis_size; ++j)
1255  {
1256  if (a != 0) // treat lucky breakdown
1257  v(j, x).equ(1. / a, *aux);
1258  else
1259  v(j, x) = 0.;
1260 
1261 
1262  preconditioner.vmult(z(j, x), v[j]);
1263  A.vmult(*aux, z[j]);
1264 
1265  // Gram-Schmidt
1266  H(0, j) = *aux * v[0];
1267  for (unsigned int i = 1; i <= j; ++i)
1268  H(i, j) = aux->add_and_dot(-H(i - 1, j), v[i - 1], v[i]);
1269  H(j + 1, j) = a = std::sqrt(aux->add_and_dot(-H(j, j), v[j], *aux));
1270 
1271  // Compute projected solution
1272 
1273  if (j > 0)
1274  {
1275  H1.reinit(j + 1, j);
1276  projected_rhs.reinit(j + 1);
1277  y.reinit(j);
1278  projected_rhs(0) = beta;
1279  H1.fill(H);
1280 
1281  // check convergence. note that the vector 'x' we pass to the
1282  // criterion is not the final solution we compute if we
1283  // decide to jump out of the iteration (we update 'x' again
1284  // right after the current loop)
1285  Householder<double> house(H1);
1286  res = house.least_squares(y, projected_rhs);
1287  iteration_state =
1288  this->iteration_status(++accumulated_iterations, res, x);
1289  if (iteration_state != SolverControl::iterate)
1290  break;
1291  }
1292  }
1293 
1294  // Update solution vector
1295  for (unsigned int j = 0; j < y.size(); ++j)
1296  x.add(y(j), z[j]);
1297  }
1298  while (iteration_state == SolverControl::iterate);
1299 
1300  // in case of failure: throw exception
1301  if (iteration_state != SolverControl::success)
1302  AssertThrow(false,
1303  SolverControl::NoConvergence(accumulated_iterations, res));
1304 }
1305 
1306 #endif // DOXYGEN
1307 
1309 
1310 #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:367
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:1649
boost::signals2::connection connect_re_orthogonalization_slot(const std::function< void(int)> &slot)
FullMatrix< double > H1
Definition: solver_gmres.h:535
boost::signals2::signal< void(double)> condition_number_signal
Definition: solver_gmres.h:326
AdditionalData additional_data
Definition: solver_gmres.h:320
static ::ExceptionBase & ExcNotInitialized()
std::vector< typename VectorMemory< VectorType >::Pointer > data
Definition: solver_gmres.h:105
#define AssertThrow(cond, exc)
Definition: exceptions.h:1531
AdditionalData additional_data
Definition: solver_gmres.h:525
FullMatrix< double > H1
Definition: solver_gmres.h:438
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:475
boost::signals2::signal< void(int)> re_orthogonalize_signal
Definition: solver_gmres.h:373
FullMatrix< double > H
Definition: solver_gmres.h:433
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:518
Stop iteration, goal reached.
void reinit(const TableIndices< N > &new_size, const bool omit_default_initialization=false)
#define Assert(cond, exc)
Definition: exceptions.h:1419
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:362
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)
AdditionalData(const unsigned int max_basis_size, const bool use_default_residual)
Definition: solver_gmres.h:484
const types::global_dof_index * dummy()
Definition: dof_handler.cc:59
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:361
size_type size() const
virtual void reinit(const size_type N, const bool omit_zeroing_entries=false)
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:530
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
#define DEAL_II_DEPRECATED
Definition: config.h:151
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:332
T max(const T &t, const MPI_Comm &mpi_communicator)
static ::ExceptionBase & ExcInternalError()