16 #ifndef dealii_solver_gmres_h 17 #define dealii_solver_gmres_h 21 #include <deal.II/base/config.h> 23 #include <deal.II/base/logstream.h> 24 #include <deal.II/base/std_cxx14/memory.h> 25 #include <deal.II/base/subscriptor.h> 27 #include <deal.II/lac/full_matrix.h> 28 #include <deal.II/lac/householder.h> 29 #include <deal.II/lac/lapack_full_matrix.h> 30 #include <deal.II/lac/solver.h> 31 #include <deal.II/lac/solver_control.h> 32 #include <deal.II/lac/vector.h> 38 DEAL_II_NAMESPACE_OPEN
48 namespace SolverGMRESImplementation
58 template <
typename VectorType>
77 VectorType &
operator[](
const unsigned int i)
const;
86 operator()(
const unsigned int i,
const VectorType &temp);
105 std::vector<typename VectorMemory<VectorType>::Pointer>
data;
177 template <
class VectorType = Vector<
double>>
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);
244 template <
typename MatrixType,
typename PreconditionerType>
246 solve(
const MatrixType & A,
248 const VectorType & b,
249 const PreconditionerType &preconditioner);
257 boost::signals2::connection
258 connect_condition_number_slot(
const std::function<
void(
double)> &slot,
259 const bool every_iteration =
false);
267 boost::signals2::connection
268 connect_eigenvalues_slot(
269 const std::function<
void(
const std::vector<std::complex<double>> &)> &slot,
270 const bool every_iteration =
false);
279 boost::signals2::connection
280 connect_hessenberg_slot(
282 const bool every_iteration =
true);
290 boost::signals2::connection
291 connect_krylov_space_slot(
301 boost::signals2::connection
302 connect_re_orthogonalization_slot(
const std::function<
void(
int)> &slot);
307 <<
"The number of temporary vectors you gave (" << arg1
308 <<
") is too small. It should be at least 10 for " 309 <<
"any results, and much more for reasonable ones.");
333 boost::signals2::signal<void(const std::vector<std::complex<double>> &)>
340 boost::signals2::signal<void(const std::vector<std::complex<double>> &)>
341 all_eigenvalues_signal;
347 boost::signals2::signal<void(const FullMatrix<double> &)> hessenberg_signal;
353 boost::signals2::signal<void(const FullMatrix<double> &)>
354 all_hessenberg_signal;
360 boost::signals2::signal<void(
398 modified_gram_schmidt(
400 & orthogonal_vectors,
401 const unsigned int dim,
402 const unsigned int accumulated_iterations,
405 bool & re_orthogonalize,
406 const boost::signals2::signal<
void(
int)> &re_orthogonalize_signal =
407 boost::signals2::signal<
void(
int)>());
416 compute_eigs_and_cond(
418 const unsigned int dim,
419 const boost::signals2::signal<
420 void(
const std::vector<std::complex<double>> &)> &eigenvalues_signal,
423 const boost::signals2::signal<
void(
double)> &cond_signal);
464 template <
class VectorType = Vector<
double>>
478 : max_basis_size(max_basis_size)
504 template <
typename MatrixType,
typename PreconditionerType>
506 solve(
const MatrixType & A,
508 const VectorType & b,
509 const PreconditionerType &preconditioner);
535 namespace SolverGMRESImplementation
537 template <
class VectorType>
546 template <
class VectorType>
558 template <
class VectorType>
561 const VectorType & temp)
564 if (
data[i] ==
nullptr)
567 data[i]->reinit(temp);
574 template <
class VectorType>
578 return (
data.size() > 0 ?
data.size() - 1 : 0);
585 complex_less_pred(
const std::complex<double> &x,
586 const std::complex<double> &y)
588 return x.real() < y.real() ||
589 (x.real() == y.real() && x.imag() < y.imag());
596 template <
class VectorType>
598 const unsigned int max_n_tmp_vectors,
599 const bool right_preconditioning,
600 const bool use_default_residual,
601 const bool force_re_orthogonalization)
602 : max_n_tmp_vectors(max_n_tmp_vectors)
603 , right_preconditioning(right_preconditioning)
604 , use_default_residual(use_default_residual)
605 , force_re_orthogonalization(force_re_orthogonalization)
607 Assert(3 <= max_n_tmp_vectors,
608 ExcMessage(
"SolverGMRES needs at least three " 609 "temporary vectors."));
614 template <
class VectorType>
619 , additional_data(data)
624 template <
class VectorType>
628 , additional_data(data)
633 template <
class VectorType>
641 for (
int i = 0; i < col; i++)
643 const double s = si(i);
644 const double c = ci(i);
645 const double dummy = h(i);
646 h(i) = c * dummy + s * h(i + 1);
647 h(i + 1) = -s * dummy + c * h(i + 1);
650 const double r = 1. / std::sqrt(h(col) * h(col) + h(col + 1) * h(col + 1));
651 si(col) = h(col + 1) * r;
652 ci(col) = h(col) * r;
653 h(col) = ci(col) * h(col) + si(col) * h(col + 1);
654 b(col + 1) = -si(col) * b(col);
660 template <
class VectorType>
664 & orthogonal_vectors,
665 const unsigned int dim,
666 const unsigned int accumulated_iterations,
669 bool & reorthogonalize,
670 const boost::signals2::signal<
void(
int)> &reorthogonalize_signal)
673 const unsigned int inner_iteration = dim - 1;
676 double norm_vv_start = 0;
677 const bool consider_reorthogonalize =
678 (reorthogonalize ==
false) && (inner_iteration % 5 == 4);
679 if (consider_reorthogonalize)
680 norm_vv_start = vv.l2_norm();
683 h(0) = vv * orthogonal_vectors[0];
684 for (
unsigned int i = 1; i < dim; ++i)
686 orthogonal_vectors[i - 1],
687 orthogonal_vectors[i]);
689 std::sqrt(vv.add_and_dot(-h(dim - 1), orthogonal_vectors[dim - 1], vv));
698 if (consider_reorthogonalize)
701 10. * norm_vv_start *
703 std::numeric_limits<typename VectorType::value_type>::epsilon()))
708 reorthogonalize =
true;
709 if (!reorthogonalize_signal.empty())
710 reorthogonalize_signal(accumulated_iterations);
714 if (reorthogonalize ==
true)
716 double htmp = vv * orthogonal_vectors[0];
718 for (
unsigned int i = 1; i < dim; ++i)
721 orthogonal_vectors[i - 1],
722 orthogonal_vectors[i]);
726 std::sqrt(vv.add_and_dot(-htmp, orthogonal_vectors[dim - 1], vv));
734 template <
class VectorType>
738 const unsigned int dim,
739 const boost::signals2::signal<
void(
const std::vector<std::complex<double>> &)>
743 const boost::signals2::signal<
void(
double)> &cond_signal)
746 if ((!eigenvalues_signal.empty() || !hessenberg_signal.empty() ||
747 !cond_signal.empty()) &&
751 for (
unsigned int i = 0; i < dim; ++i)
752 for (
unsigned int j = 0; j < dim; ++j)
753 mat(i, j) = H_orig(i, j);
754 hessenberg_signal(H_orig);
756 if (!eigenvalues_signal.empty())
763 std::vector<std::complex<double>>
eigenvalues(dim);
764 for (
unsigned int i = 0; i < mat_eig.
n(); ++i)
767 std::sort(eigenvalues.begin(),
769 internal::SolverGMRESImplementation::complex_less_pred);
770 eigenvalues_signal(eigenvalues);
774 if (!cond_signal.empty() && (mat.
n() > 1))
777 double condition_number =
779 cond_signal(condition_number);
786 template <
class VectorType>
787 template <
typename MatrixType,
typename PreconditionerType>
791 const VectorType & b,
792 const PreconditionerType &preconditioner)
802 const unsigned int n_tmp_vectors =
807 n_tmp_vectors, this->memory);
812 unsigned int accumulated_iterations = 0;
814 const bool do_eigenvalues =
815 !condition_number_signal.empty() || !all_condition_numbers_signal.empty() ||
816 !eigenvalues_signal.empty() || !all_eigenvalues_signal.empty() ||
817 !hessenberg_signal.empty() || !all_hessenberg_signal.empty();
822 H_orig.
reinit(n_tmp_vectors, n_tmp_vectors - 1);
825 H.reinit(n_tmp_vectors, n_tmp_vectors - 1);
829 si(n_tmp_vectors - 1), h(n_tmp_vectors - 1);
832 unsigned int dim = 0;
835 double last_res = -std::numeric_limits<double>::max();
847 VectorType &v = tmp_vectors(0, x);
848 VectorType &p = tmp_vectors(n_tmp_vectors - 1, x);
854 std::unique_ptr<::Vector<double>> gamma_;
855 if (!use_default_residual)
862 gamma_ = std_cxx14::make_unique<::Vector<double>>(gamma.size());
874 h.
reinit(n_tmp_vectors - 1);
876 if (left_precondition)
880 preconditioner.vmult(v, p);
888 double rho = v.l2_norm();
893 if (use_default_residual)
897 this->iteration_status(accumulated_iterations, rho, x);
904 deallog <<
"default_res=" << rho << std::endl;
906 if (left_precondition)
912 preconditioner.vmult(*r, v);
914 double res = r->l2_norm();
917 this->iteration_status(accumulated_iterations, res, x);
930 for (
unsigned int inner_iteration = 0;
931 ((inner_iteration < n_tmp_vectors - 2) &&
935 ++accumulated_iterations;
937 VectorType &vv = tmp_vectors(inner_iteration + 1, x);
939 if (left_precondition)
941 A.vmult(p, tmp_vectors[inner_iteration]);
942 preconditioner.vmult(vv, p);
946 preconditioner.vmult(p, tmp_vectors[inner_iteration]);
950 dim = inner_iteration + 1;
952 const double s = modified_gram_schmidt(tmp_vectors,
954 accumulated_iterations,
958 re_orthogonalize_signal);
959 h(inner_iteration + 1) = s;
969 for (
unsigned int i = 0; i < dim + 1; ++i)
970 H_orig(i, inner_iteration) = h(i);
973 givens_rotation(h, gamma, ci, si, inner_iteration);
976 for (
unsigned int i = 0; i < dim; ++i)
977 H(i, inner_iteration) = h(i);
980 rho = std::fabs(gamma(dim));
982 if (use_default_residual)
986 this->iteration_status(accumulated_iterations, rho, x);
990 deallog <<
"default_res=" << rho << std::endl;
995 H1.reinit(dim + 1, dim);
997 for (
unsigned int i = 0; i < dim + 1; ++i)
998 for (
unsigned int j = 0; j < dim; ++j)
1001 H1.backward(h_, *gamma_);
1003 if (left_precondition)
1004 for (
unsigned int i = 0; i < dim; ++i)
1005 x_->add(h_(i), tmp_vectors[i]);
1009 for (
unsigned int i = 0; i < dim; ++i)
1010 p.add(h_(i), tmp_vectors[i]);
1011 preconditioner.vmult(*r, p);
1015 r->sadd(-1., 1., b);
1017 if (left_precondition)
1019 const double res = r->l2_norm();
1023 this->iteration_status(accumulated_iterations, res, x);
1027 preconditioner.vmult(*x_, *r);
1028 const double preconditioned_res = x_->l2_norm();
1029 last_res = preconditioned_res;
1032 this->iteration_status(accumulated_iterations,
1041 H1.reinit(dim + 1, dim);
1043 for (
unsigned int i = 0; i < dim + 1; ++i)
1044 for (
unsigned int j = 0; j < dim; ++j)
1047 compute_eigs_and_cond(H_orig,
1049 all_eigenvalues_signal,
1050 all_hessenberg_signal,
1051 condition_number_signal);
1053 H1.backward(h, gamma);
1055 if (left_precondition)
1056 for (
unsigned int i = 0; i < dim; ++i)
1057 x.add(h(i), tmp_vectors[i]);
1061 for (
unsigned int i = 0; i < dim; ++i)
1062 p.add(h(i), tmp_vectors[i]);
1063 preconditioner.vmult(v, p);
1071 compute_eigs_and_cond(H_orig,
1075 condition_number_signal);
1077 if (!krylov_space_signal.empty())
1078 krylov_space_signal(tmp_vectors);
1087 template <
class VectorType>
1088 boost::signals2::connection
1090 const std::function<
void(
double)> &slot,
1091 const bool every_iteration)
1093 if (every_iteration)
1095 return all_condition_numbers_signal.connect(slot);
1099 return condition_number_signal.connect(slot);
1105 template <
class VectorType>
1106 boost::signals2::connection
1108 const std::function<
void(
const std::vector<std::complex<double>> &)> &slot,
1109 const bool every_iteration)
1111 if (every_iteration)
1113 return all_eigenvalues_signal.connect(slot);
1117 return eigenvalues_signal.connect(slot);
1123 template <
class VectorType>
1124 boost::signals2::connection
1127 const bool every_iteration)
1129 if (every_iteration)
1131 return all_hessenberg_signal.connect(slot);
1135 return hessenberg_signal.connect(slot);
1141 template <
class VectorType>
1142 boost::signals2::connection
1144 const std::function<
void(
1147 return krylov_space_signal.connect(slot);
1152 template <
class VectorType>
1153 boost::signals2::connection
1155 const std::function<
void(
int)> &slot)
1157 return re_orthogonalize_signal.connect(slot);
1162 template <
class VectorType>
1175 template <
class VectorType>
1180 , additional_data(data)
1185 template <
class VectorType>
1189 , additional_data(data)
1194 template <
class VectorType>
1195 template <
typename MatrixType,
typename PreconditionerType>
1199 const VectorType & b,
1200 const PreconditionerType &preconditioner)
1206 const unsigned int basis_size = additional_data.max_basis_size;
1210 basis_size, this->memory);
1212 basis_size, this->memory);
1216 unsigned int accumulated_iterations = 0;
1219 H.reinit(basis_size + 1, basis_size);
1226 double res = -std::numeric_limits<double>::max();
1233 aux->sadd(-1., 1., b);
1235 double beta = aux->l2_norm();
1237 iteration_state = this->iteration_status(accumulated_iterations, res, x);
1241 H.reinit(basis_size + 1, basis_size);
1244 for (
unsigned int j = 0; j < basis_size; ++j)
1247 v(j, x).equ(1. / a, *aux);
1252 preconditioner.vmult(z(j, x), v[j]);
1253 A.vmult(*aux, z[j]);
1256 H(0, j) = *aux * v[0];
1257 for (
unsigned int i = 1; i <= j; ++i)
1258 H(i, j) = aux->add_and_dot(-H(i - 1, j), v[i - 1], v[i]);
1259 H(j + 1, j) = a = std::sqrt(aux->add_and_dot(-H(j, j), v[j], *aux));
1265 H1.reinit(j + 1, j);
1266 projected_rhs.
reinit(j + 1);
1268 projected_rhs(0) = beta;
1278 this->iteration_status(++accumulated_iterations, res, x);
1285 for (
unsigned int j = 0; j < y.
size(); ++j)
1298 DEAL_II_NAMESPACE_CLOSE
unsigned int size() const
SolverFGMRES(SolverControl &cn, VectorMemory< VectorType > &mem, const AdditionalData &data=AdditionalData())
Number add_and_dot(const Number a, const Vector< Number > &V, const Vector< Number > &W)
TmpVectors(const unsigned int max_size, VectorMemory< VectorType > &vmem)
number singular_value(const size_type i) const
boost::signals2::connection connect_condition_number_slot(const std::function< void(double)> &slot, const bool every_iteration=false)
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)
boost::signals2::signal< void(const internal::SolverGMRESImplementation::TmpVectors< VectorType > &)> krylov_space_signal
void solve(const MatrixType &A, VectorType &x, const VectorType &b, const PreconditionerType &preconditioner)
boost::signals2::connection connect_re_orthogonalization_slot(const std::function< void(int)> &slot)
boost::signals2::signal< void(double)> condition_number_signal
AdditionalData additional_data
static::ExceptionBase & ExcNotInitialized()
std::vector< typename VectorMemory< VectorType >::Pointer > data
#define AssertThrow(cond, exc)
static::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
AdditionalData additional_data
VectorType & operator()(const unsigned int i, const VectorType &temp)
boost::signals2::signal< void(int)> re_orthogonalize_signal
bool right_preconditioning
void solve(const MatrixType &A, VectorType &x, const VectorType &b, const PreconditionerType &preconditioner)
static::ExceptionBase & ExcMessage(std::string arg1)
VectorType & operator[](const unsigned int i) const
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)
Stop iteration, goal reached.
void reinit(const TableIndices< N > &new_size, const bool omit_default_initialization=false)
#define Assert(cond, exc)
boost::signals2::connection connect_krylov_space_slot(const std::function< void(const internal::SolverGMRESImplementation::TmpVectors< VectorType > &)> &slot)
unsigned int max_basis_size
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)>())
bool force_re_orthogonalization
void givens_rotation(Vector< double > &h, Vector< double > &b, Vector< double > &ci, Vector< double > &si, int col) const
VectorMemory< VectorType > & mem
virtual double criterion()
virtual void reinit(const size_type N, const bool omit_zeroing_entries=false)
void compute_eigenvalues(const bool right_eigenvectors=false, const bool left_eigenvectors=false)
std::complex< number > eigenvalue(const size_type i) const
bool use_default_residual
double least_squares(Vector< number2 > &dst, const Vector< number2 > &src) const
boost::signals2::connection connect_hessenberg_slot(const std::function< void(const FullMatrix< double > &)> &slot, const bool every_iteration=true)
AdditionalData(const unsigned int max_basis_size=30, const bool=true)
SolverGMRES(SolverControl &cn, VectorMemory< VectorType > &mem, 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)
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
unsigned int max_n_tmp_vectors
static::ExceptionBase & ExcInternalError()