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