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