16 #ifndef dealii_solver_gmres_h
17 #define dealii_solver_gmres_h
50 namespace SolverGMRESImplementation
60 template <
typename VectorType>
89 operator()(
const unsigned int i,
const VectorType &temp);
108 std::vector<typename VectorMemory<VectorType>::Pointer>
data;
177 template <
class VectorType = Vector<
double>>
249 template <
typename MatrixType,
typename PreconditionerType>
253 const VectorType &
b,
254 const PreconditionerType &preconditioner);
262 boost::signals2::connection
264 const bool every_iteration =
false);
272 boost::signals2::connection
274 const std::function<
void(
const std::vector<std::complex<double>> &)> &slot,
275 const bool every_iteration =
false);
284 boost::signals2::connection
287 const bool every_iteration =
true);
295 boost::signals2::connection
306 boost::signals2::connection
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.");
338 boost::signals2::signal<void(
const std::vector<std::complex<double>> &)>
345 boost::signals2::signal<void(
const std::vector<std::complex<double>> &)>
365 boost::signals2::signal<void(
405 & orthogonal_vectors,
406 const unsigned int dim,
407 const unsigned int accumulated_iterations,
410 bool & re_orthogonalize,
412 boost::signals2::signal<
void(
int)>());
423 const unsigned int dim,
424 const boost::signals2::signal<
428 const boost::signals2::signal<
void(
double)> &cond_signal);
461 template <
class VectorType = Vector<
double>>
500 template <
typename MatrixType,
typename PreconditionerType>
504 const VectorType &
b,
505 const PreconditionerType &preconditioner);
531 namespace SolverGMRESImplementation
533 template <
class VectorType>
542 template <
class VectorType>
554 template <
class VectorType>
557 const VectorType & temp)
560 if (data[i] ==
nullptr)
563 data[i]->reinit(temp);
570 template <
class VectorType>
574 return (data.size() > 0 ? data.size() - 1 : 0);
581 complex_less_pred(
const std::complex<double> &x,
582 const std::complex<double> &y)
584 return x.real() < y.real() ||
585 (x.real() == y.real() && x.imag() < y.imag());
592 template <
class VectorType>
594 const unsigned int max_n_tmp_vectors,
595 const bool right_preconditioning,
596 const bool use_default_residual,
597 const bool force_re_orthogonalization)
598 : max_n_tmp_vectors(max_n_tmp_vectors)
599 , right_preconditioning(right_preconditioning)
600 , use_default_residual(use_default_residual)
601 , force_re_orthogonalization(force_re_orthogonalization)
603 Assert(3 <= max_n_tmp_vectors,
604 ExcMessage(
"SolverGMRES needs at least three "
605 "temporary vectors."));
610 template <
class VectorType>
613 const AdditionalData & data)
615 , additional_data(data)
620 template <
class VectorType>
622 const AdditionalData &data)
624 , additional_data(data)
629 template <
class VectorType>
637 for (
int i = 0; i < col; ++i)
639 const double s = si(i);
640 const double c = ci(i);
641 const double dummy = h(i);
642 h(i) = c * dummy + s * h(i + 1);
643 h(i + 1) = -s * dummy + c * h(i + 1);
646 const double r = 1. /
std::sqrt(h(col) * h(col) + h(col + 1) * h(col + 1));
647 si(col) = h(col + 1) * r;
648 ci(col) = h(col) * r;
649 h(col) = ci(col) * h(col) + si(col) * h(col + 1);
650 b(col + 1) = -si(col) *
b(col);
656 template <
class VectorType>
660 & orthogonal_vectors,
661 const unsigned int dim,
662 const unsigned int accumulated_iterations,
665 bool & reorthogonalize,
666 const boost::signals2::signal<
void(
int)> &reorthogonalize_signal)
669 const unsigned int inner_iteration = dim - 1;
672 double norm_vv_start = 0;
673 const bool consider_reorthogonalize =
674 (reorthogonalize ==
false) && (inner_iteration % 5 == 4);
675 if (consider_reorthogonalize)
676 norm_vv_start = vv.l2_norm();
679 h(0) = vv * orthogonal_vectors[0];
680 for (
unsigned int i = 1; i < dim; ++i)
682 orthogonal_vectors[i - 1],
683 orthogonal_vectors[i]);
685 std::sqrt(vv.add_and_dot(-h(dim - 1), orthogonal_vectors[dim - 1], vv));
694 if (consider_reorthogonalize)
697 10. * norm_vv_start *
704 reorthogonalize =
true;
705 if (!reorthogonalize_signal.empty())
706 reorthogonalize_signal(accumulated_iterations);
710 if (reorthogonalize ==
true)
712 double htmp = vv * orthogonal_vectors[0];
714 for (
unsigned int i = 1; i < dim; ++i)
716 htmp = vv.add_and_dot(-htmp,
717 orthogonal_vectors[i - 1],
718 orthogonal_vectors[i]);
722 std::sqrt(vv.add_and_dot(-htmp, orthogonal_vectors[dim - 1], vv));
730 template <
class VectorType>
734 const unsigned int dim,
735 const boost::signals2::signal<
void(
const std::vector<std::complex<double>> &)>
739 const boost::signals2::signal<
void(
double)> &cond_signal)
742 if ((!eigenvalues_signal.empty() || !hessenberg_signal.empty() ||
743 !cond_signal.empty()) &&
747 for (
unsigned int i = 0; i < dim; ++i)
748 for (
unsigned int j = 0; j < dim; ++j)
749 mat(i, j) = H_orig(i, j);
750 hessenberg_signal(H_orig);
752 if (!eigenvalues_signal.empty())
758 mat_eig.compute_eigenvalues();
759 std::vector<std::complex<double>>
eigenvalues(dim);
760 for (
unsigned int i = 0; i < mat_eig.n(); ++i)
765 internal::SolverGMRESImplementation::complex_less_pred);
770 if (!cond_signal.empty() && (mat.n() > 1))
773 double condition_number =
774 mat.singular_value(0) / mat.singular_value(mat.n() - 1);
775 cond_signal(condition_number);
782 template <
class VectorType>
783 template <
typename MatrixType,
typename PreconditionerType>
787 const VectorType &
b,
788 const PreconditionerType &preconditioner)
798 const unsigned int n_tmp_vectors =
803 n_tmp_vectors, this->memory);
808 unsigned int accumulated_iterations = 0;
810 const bool do_eigenvalues =
811 !condition_number_signal.empty() || !all_condition_numbers_signal.empty() ||
812 !eigenvalues_signal.empty() || !all_eigenvalues_signal.empty() ||
813 !hessenberg_signal.empty() || !all_hessenberg_signal.empty();
818 H_orig.reinit(n_tmp_vectors, n_tmp_vectors - 1);
821 H.reinit(n_tmp_vectors, n_tmp_vectors - 1);
825 si(n_tmp_vectors - 1), h(n_tmp_vectors - 1);
828 unsigned int dim = 0;
831 double last_res = std::numeric_limits<double>::lowest();
843 VectorType &v = tmp_vectors(0, x);
844 VectorType &p = tmp_vectors(n_tmp_vectors - 1, x);
850 std::unique_ptr<::Vector<double>> gamma_;
851 if (!use_default_residual)
858 gamma_ = std::make_unique<::Vector<double>>(
gamma.size());
870 h.
reinit(n_tmp_vectors - 1);
872 if (left_precondition)
876 preconditioner.vmult(v, p);
884 double rho = v.l2_norm();
889 if (use_default_residual)
893 this->iteration_status(accumulated_iterations, rho, x);
900 deallog <<
"default_res=" << rho << std::endl;
902 if (left_precondition)
908 preconditioner.vmult(*r, v);
910 double res = r->l2_norm();
913 this->iteration_status(accumulated_iterations, res, x);
926 for (
unsigned int inner_iteration = 0;
927 ((inner_iteration < n_tmp_vectors - 2) &&
931 ++accumulated_iterations;
933 VectorType &vv = tmp_vectors(inner_iteration + 1, x);
935 if (left_precondition)
937 A.vmult(p, tmp_vectors[inner_iteration]);
938 preconditioner.vmult(vv, p);
942 preconditioner.vmult(p, tmp_vectors[inner_iteration]);
946 dim = inner_iteration + 1;
948 const double s = modified_gram_schmidt(tmp_vectors,
950 accumulated_iterations,
954 re_orthogonalize_signal);
955 h(inner_iteration + 1) = s;
965 for (
unsigned int i = 0; i < dim + 1; ++i)
966 H_orig(i, inner_iteration) = h(i);
972 for (
unsigned int i = 0; i < dim; ++i)
973 H(i, inner_iteration) = h(i);
978 if (use_default_residual)
982 this->iteration_status(accumulated_iterations, rho, x);
986 deallog <<
"default_res=" << rho << std::endl;
991 H1.reinit(dim + 1, dim);
993 for (
unsigned int i = 0; i < dim + 1; ++i)
994 for (
unsigned int j = 0; j < dim; ++j)
997 H1.backward(h_, *gamma_);
999 if (left_precondition)
1000 for (
unsigned int i = 0; i < dim; ++i)
1001 x_->add(h_(i), tmp_vectors[i]);
1005 for (
unsigned int i = 0; i < dim; ++i)
1006 p.add(h_(i), tmp_vectors[i]);
1007 preconditioner.vmult(*r, p);
1011 r->sadd(-1., 1.,
b);
1013 if (left_precondition)
1015 const double res = r->l2_norm();
1019 this->iteration_status(accumulated_iterations, res, x);
1023 preconditioner.vmult(*x_, *r);
1024 const double preconditioned_res = x_->l2_norm();
1025 last_res = preconditioned_res;
1028 this->iteration_status(accumulated_iterations,
1037 H1.reinit(dim + 1, dim);
1039 for (
unsigned int i = 0; i < dim + 1; ++i)
1040 for (
unsigned int j = 0; j < dim; ++j)
1043 compute_eigs_and_cond(H_orig,
1045 all_eigenvalues_signal,
1046 all_hessenberg_signal,
1047 condition_number_signal);
1049 H1.backward(h,
gamma);
1051 if (left_precondition)
1052 for (
unsigned int i = 0; i < dim; ++i)
1053 x.add(h(i), tmp_vectors[i]);
1057 for (
unsigned int i = 0; i < dim; ++i)
1058 p.add(h(i), tmp_vectors[i]);
1059 preconditioner.vmult(v, p);
1067 compute_eigs_and_cond(H_orig,
1071 condition_number_signal);
1073 if (!krylov_space_signal.empty())
1074 krylov_space_signal(tmp_vectors);
1083 template <
class VectorType>
1084 boost::signals2::connection
1086 const std::function<
void(
double)> &slot,
1087 const bool every_iteration)
1089 if (every_iteration)
1091 return all_condition_numbers_signal.connect(slot);
1095 return condition_number_signal.connect(slot);
1101 template <
class VectorType>
1102 boost::signals2::connection
1104 const std::function<
void(
const std::vector<std::complex<double>> &)> &slot,
1105 const bool every_iteration)
1107 if (every_iteration)
1109 return all_eigenvalues_signal.connect(slot);
1113 return eigenvalues_signal.connect(slot);
1119 template <
class VectorType>
1120 boost::signals2::connection
1123 const bool every_iteration)
1125 if (every_iteration)
1127 return all_hessenberg_signal.connect(slot);
1131 return hessenberg_signal.connect(slot);
1137 template <
class VectorType>
1138 boost::signals2::connection
1140 const std::function<
void(
1143 return krylov_space_signal.connect(slot);
1148 template <
class VectorType>
1149 boost::signals2::connection
1151 const std::function<
void(
int)> &slot)
1153 return re_orthogonalize_signal.connect(slot);
1158 template <
class VectorType>
1171 template <
class VectorType>
1174 const AdditionalData & data)
1176 , additional_data(data)
1181 template <
class VectorType>
1183 const AdditionalData &data)
1185 , additional_data(data)
1190 template <
class VectorType>
1191 template <
typename MatrixType,
typename PreconditionerType>
1195 const VectorType &
b,
1196 const PreconditionerType &preconditioner)
1202 const unsigned int basis_size = additional_data.max_basis_size;
1206 basis_size, this->memory);
1208 basis_size, this->memory);
1212 unsigned int accumulated_iterations = 0;
1215 H.reinit(basis_size + 1, basis_size);
1222 double res = std::numeric_limits<double>::lowest();
1229 aux->sadd(-1., 1.,
b);
1231 double beta = aux->l2_norm();
1233 iteration_state = this->iteration_status(accumulated_iterations, res, x);
1237 H.reinit(basis_size + 1, basis_size);
1240 for (
unsigned int j = 0; j < basis_size; ++j)
1243 v(j, x).equ(1. / a, *aux);
1248 preconditioner.vmult(z(j, x), v[j]);
1249 A.vmult(*aux, z[j]);
1252 H(0, j) = *aux * v[0];
1253 for (
unsigned int i = 1; i <= j; ++i)
1254 H(i, j) = aux->add_and_dot(-H(i - 1, j), v[i - 1], v[i]);
1255 H(j + 1, j) = a =
std::sqrt(aux->add_and_dot(-H(j, j), v[j], *aux));
1261 H1.reinit(j + 1, j);
1262 projected_rhs.
reinit(j + 1);
1264 projected_rhs(0) = beta;
1272 res = house.least_squares(y, projected_rhs);
1274 this->iteration_status(++accumulated_iterations, res, x);
1281 for (
unsigned int j = 0; j < y.
size(); ++j)
@ 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)
SolverFGMRES(SolverControl &cn, const AdditionalData &data=AdditionalData())
AdditionalData additional_data
boost::signals2::signal< void(const FullMatrix< double > &)> hessenberg_signal
AdditionalData additional_data
boost::signals2::signal< void(double)> condition_number_signal
SolverGMRES(const SolverGMRES< VectorType > &)=delete
boost::signals2::signal< void(const std::vector< std::complex< double >> &)> eigenvalues_signal
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
boost::signals2::signal< void(const FullMatrix< double > &)> all_hessenberg_signal
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::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
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())
boost::signals2::connection connect_eigenvalues_slot(const std::function< void(const std::vector< std::complex< double >> &)> &slot, const bool every_iteration=false)
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
boost::signals2::signal< void(int)> re_orthogonalize_signal
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)
Number add_and_dot(const Number a, const Vector< Number > &V, const Vector< Number > &W)
virtual void reinit(const size_type N, const bool omit_zeroing_entries=false)
VectorMemory< VectorType > & mem
std::vector< typename VectorMemory< VectorType >::Pointer > data
VectorType & operator()(const unsigned int i, const VectorType &temp)
TmpVectors(const unsigned int max_size, VectorMemory< VectorType > &vmem)
unsigned int size() const
VectorType & operator[](const unsigned int i) const
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcNotInitialized()
#define Assert(cond, exc)
static ::ExceptionBase & ExcTooFewTmpVectors(int arg1)
#define AssertIndexRange(index, range)
#define DeclException1(Exception1, type1, outsequence)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Expression fabs(const Expression &x)
SymmetricTensor< 2, dim, Number > epsilon(const Tensor< 2, dim, Number > &Grad_u)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
std::array< NumberType, 3 > givens_rotation(const NumberType &x, const NumberType &y)
long double gamma(const unsigned int n)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
AdditionalData(const unsigned int max_basis_size=30)
unsigned int max_basis_size
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)
bool right_preconditioning
bool force_re_orthogonalization
bool use_default_residual
unsigned int max_n_tmp_vectors