16 #ifndef dealii_solver_gmres_h 17 #define dealii_solver_gmres_h 48 namespace SolverGMRESImplementation
58 template <
typename VectorType>
105 std::vector<typename VectorMemory<VectorType>::Pointer>
data;
174 template <
class VectorType = Vector<
double>>
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);
246 template <
typename MatrixType,
typename PreconditionerType>
248 solve(
const MatrixType &
A,
251 const PreconditionerType &preconditioner);
259 boost::signals2::connection
260 connect_condition_number_slot(
const std::function<
void(
double)> &slot,
261 const bool every_iteration =
false);
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);
281 boost::signals2::connection
282 connect_hessenberg_slot(
284 const bool every_iteration =
true);
292 boost::signals2::connection
293 connect_krylov_space_slot(
303 boost::signals2::connection
304 connect_re_orthogonalization_slot(
const std::function<
void(
int)> &slot);
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.");
335 boost::signals2::signal<void(const std::vector<std::complex<double>> &)>
342 boost::signals2::signal<void(const std::vector<std::complex<double>> &)>
343 all_eigenvalues_signal;
349 boost::signals2::signal<void(const FullMatrix<double> &)> hessenberg_signal;
355 boost::signals2::signal<void(const FullMatrix<double> &)>
356 all_hessenberg_signal;
362 boost::signals2::signal<void(
400 modified_gram_schmidt(
402 & orthogonal_vectors,
403 const unsigned int dim,
404 const unsigned int accumulated_iterations,
407 bool & re_orthogonalize,
408 const boost::signals2::signal<
void(
int)> &re_orthogonalize_signal =
409 boost::signals2::signal<
void(
int)>());
418 compute_eigs_and_cond(
420 const unsigned int dim,
421 const boost::signals2::signal<
422 void(
const std::vector<std::complex<double>> &)> &eigenvalues_signal,
425 const boost::signals2::signal<
void(
double)> &cond_signal);
458 template <
class VectorType = Vector<
double>>
471 : max_basis_size(max_basis_size)
497 template <
typename MatrixType,
typename PreconditionerType>
499 solve(
const MatrixType &
A,
502 const PreconditionerType &preconditioner);
528 namespace SolverGMRESImplementation
530 template <
class VectorType>
539 template <
class VectorType>
551 template <
class VectorType>
557 if (
data[i] ==
nullptr)
560 data[i]->reinit(temp);
567 template <
class VectorType>
571 return (
data.size() > 0 ?
data.size() - 1 : 0);
578 complex_less_pred(
const std::complex<double> &x,
579 const std::complex<double> &y)
581 return x.real() < y.real() ||
582 (x.real() == y.real() && x.imag() < y.imag());
589 template <
class VectorType>
591 const unsigned int max_n_tmp_vectors,
592 const bool right_preconditioning,
593 const bool use_default_residual,
594 const bool force_re_orthogonalization)
595 : max_n_tmp_vectors(max_n_tmp_vectors)
596 , right_preconditioning(right_preconditioning)
597 , use_default_residual(use_default_residual)
598 , force_re_orthogonalization(force_re_orthogonalization)
600 Assert(3 <= max_n_tmp_vectors,
601 ExcMessage(
"SolverGMRES needs at least three " 602 "temporary vectors."));
607 template <
class VectorType>
612 , additional_data(data)
617 template <
class VectorType>
621 , additional_data(data)
626 template <
class VectorType>
634 for (
int i = 0; i < col; i++)
636 const double s = si(i);
637 const double c = ci(i);
638 const double dummy = h(i);
639 h(i) = c * dummy + s * h(i + 1);
640 h(i + 1) = -s * dummy + c * h(i + 1);
643 const double r = 1. / std::sqrt(h(col) * h(col) + h(col + 1) * h(col + 1));
644 si(col) = h(col + 1) * r;
645 ci(col) = h(col) * r;
646 h(col) = ci(col) * h(col) + si(col) * h(col + 1);
647 b(col + 1) = -si(col) *
b(col);
653 template <
class VectorType>
657 & orthogonal_vectors,
658 const unsigned int dim,
659 const unsigned int accumulated_iterations,
662 bool & reorthogonalize,
663 const boost::signals2::signal<
void(
int)> &reorthogonalize_signal)
666 const unsigned int inner_iteration = dim - 1;
669 double norm_vv_start = 0;
670 const bool consider_reorthogonalize =
671 (reorthogonalize ==
false) && (inner_iteration % 5 == 4);
672 if (consider_reorthogonalize)
673 norm_vv_start = vv.l2_norm();
676 h(0) = vv * orthogonal_vectors[0];
677 for (
unsigned int i = 1; i < dim; ++i)
678 h(i) = vv.add_and_dot(-h(i - 1),
679 orthogonal_vectors[i - 1],
680 orthogonal_vectors[i]);
682 std::sqrt(vv.add_and_dot(-h(dim - 1), orthogonal_vectors[dim - 1], vv));
691 if (consider_reorthogonalize)
694 10. * norm_vv_start *
701 reorthogonalize =
true;
702 if (!reorthogonalize_signal.empty())
703 reorthogonalize_signal(accumulated_iterations);
707 if (reorthogonalize ==
true)
709 double htmp = vv * orthogonal_vectors[0];
711 for (
unsigned int i = 1; i < dim; ++i)
713 htmp = vv.add_and_dot(-htmp,
714 orthogonal_vectors[i - 1],
715 orthogonal_vectors[i]);
719 std::sqrt(vv.add_and_dot(-htmp, orthogonal_vectors[dim - 1], vv));
727 template <
class VectorType>
731 const unsigned int dim,
732 const boost::signals2::signal<
void(
const std::vector<std::complex<double>> &)>
736 const boost::signals2::signal<
void(
double)> &cond_signal)
739 if ((!eigenvalues_signal.empty() || !hessenberg_signal.empty() ||
740 !cond_signal.empty()) &&
744 for (
unsigned int i = 0; i < dim; ++i)
745 for (
unsigned int j = 0; j < dim; ++j)
746 mat(i, j) = H_orig(i, j);
747 hessenberg_signal(H_orig);
749 if (!eigenvalues_signal.empty())
756 std::vector<std::complex<double>>
eigenvalues(dim);
757 for (
unsigned int i = 0; i < mat_eig.
n(); ++i)
760 std::sort(eigenvalues.begin(),
762 internal::SolverGMRESImplementation::complex_less_pred);
763 eigenvalues_signal(eigenvalues);
767 if (!cond_signal.empty() && (mat.
n() > 1))
770 double condition_number =
772 cond_signal(condition_number);
779 template <
class VectorType>
780 template <
typename MatrixType,
typename PreconditionerType>
785 const PreconditionerType &preconditioner)
795 const unsigned int n_tmp_vectors =
800 n_tmp_vectors, this->memory);
805 unsigned int accumulated_iterations = 0;
807 const bool do_eigenvalues =
808 !condition_number_signal.empty() || !all_condition_numbers_signal.empty() ||
809 !eigenvalues_signal.empty() || !all_eigenvalues_signal.empty() ||
810 !hessenberg_signal.empty() || !all_hessenberg_signal.empty();
815 H_orig.
reinit(n_tmp_vectors, n_tmp_vectors - 1);
818 H.reinit(n_tmp_vectors, n_tmp_vectors - 1);
821 ::Vector<double>
gamma(n_tmp_vectors), ci(n_tmp_vectors - 1),
822 si(n_tmp_vectors - 1), h(n_tmp_vectors - 1);
825 unsigned int dim = 0;
841 VectorType &p = tmp_vectors(n_tmp_vectors - 1, x);
847 std::unique_ptr<::Vector<double>> gamma_;
848 if (!use_default_residual)
855 gamma_ = std::make_unique<::Vector<double>>(
gamma.size());
867 h.reinit(n_tmp_vectors - 1);
869 if (left_precondition)
873 preconditioner.vmult(v, p);
881 double rho = v.l2_norm();
886 if (use_default_residual)
890 this->iteration_status(accumulated_iterations, rho, x);
897 deallog <<
"default_res=" << rho << std::endl;
899 if (left_precondition)
905 preconditioner.vmult(*r, v);
907 double res = r->l2_norm();
910 this->iteration_status(accumulated_iterations, res, x);
923 for (
unsigned int inner_iteration = 0;
924 ((inner_iteration < n_tmp_vectors - 2) &&
928 ++accumulated_iterations;
930 VectorType &vv = tmp_vectors(inner_iteration + 1, x);
932 if (left_precondition)
934 A.vmult(p, tmp_vectors[inner_iteration]);
935 preconditioner.vmult(vv, p);
939 preconditioner.vmult(p, tmp_vectors[inner_iteration]);
943 dim = inner_iteration + 1;
945 const double s = modified_gram_schmidt(tmp_vectors,
947 accumulated_iterations,
951 re_orthogonalize_signal);
952 h(inner_iteration + 1) = s;
962 for (
unsigned int i = 0; i < dim + 1; ++i)
963 H_orig(i, inner_iteration) = h(i);
969 for (
unsigned int i = 0; i < dim; ++i)
970 H(i, inner_iteration) = h(i);
975 if (use_default_residual)
979 this->iteration_status(accumulated_iterations, rho, x);
983 deallog <<
"default_res=" << rho << std::endl;
985 ::Vector<double> h_(dim);
988 H1.reinit(dim + 1, dim);
990 for (
unsigned int i = 0; i < dim + 1; ++i)
991 for (
unsigned int j = 0; j < dim; ++j)
994 H1.backward(h_, *gamma_);
996 if (left_precondition)
997 for (
unsigned int i = 0; i < dim; ++i)
998 x_->add(h_(i), tmp_vectors[i]);
1002 for (
unsigned int i = 0; i < dim; ++i)
1003 p.add(h_(i), tmp_vectors[i]);
1004 preconditioner.vmult(*r, p);
1008 r->sadd(-1., 1., b);
1010 if (left_precondition)
1012 const double res = r->l2_norm();
1016 this->iteration_status(accumulated_iterations, res, x);
1020 preconditioner.vmult(*x_, *r);
1021 const double preconditioned_res = x_->l2_norm();
1022 last_res = preconditioned_res;
1025 this->iteration_status(accumulated_iterations,
1034 H1.reinit(dim + 1, dim);
1036 for (
unsigned int i = 0; i < dim + 1; ++i)
1037 for (
unsigned int j = 0; j < dim; ++j)
1040 compute_eigs_and_cond(H_orig,
1042 all_eigenvalues_signal,
1043 all_hessenberg_signal,
1044 condition_number_signal);
1046 H1.backward(h,
gamma);
1048 if (left_precondition)
1049 for (
unsigned int i = 0; i < dim; ++i)
1050 x.add(h(i), tmp_vectors[i]);
1054 for (
unsigned int i = 0; i < dim; ++i)
1055 p.add(h(i), tmp_vectors[i]);
1056 preconditioner.vmult(v, p);
1064 compute_eigs_and_cond(H_orig,
1068 condition_number_signal);
1070 if (!krylov_space_signal.empty())
1071 krylov_space_signal(tmp_vectors);
1080 template <
class VectorType>
1081 boost::signals2::connection
1083 const std::function<
void(
double)> &slot,
1084 const bool every_iteration)
1086 if (every_iteration)
1088 return all_condition_numbers_signal.connect(slot);
1092 return condition_number_signal.connect(slot);
1098 template <
class VectorType>
1099 boost::signals2::connection
1101 const std::function<
void(
const std::vector<std::complex<double>> &)> &slot,
1102 const bool every_iteration)
1104 if (every_iteration)
1106 return all_eigenvalues_signal.connect(slot);
1110 return eigenvalues_signal.connect(slot);
1116 template <
class VectorType>
1117 boost::signals2::connection
1120 const bool every_iteration)
1122 if (every_iteration)
1124 return all_hessenberg_signal.connect(slot);
1128 return hessenberg_signal.connect(slot);
1134 template <
class VectorType>
1135 boost::signals2::connection
1137 const std::function<
void(
1140 return krylov_space_signal.connect(slot);
1145 template <
class VectorType>
1146 boost::signals2::connection
1148 const std::function<
void(
int)> &slot)
1150 return re_orthogonalize_signal.connect(slot);
1155 template <
class VectorType>
1168 template <
class VectorType>
1173 , additional_data(data)
1178 template <
class VectorType>
1182 , additional_data(data)
1187 template <
class VectorType>
1188 template <
typename MatrixType,
typename PreconditionerType>
1193 const PreconditionerType &preconditioner)
1199 const unsigned int basis_size = additional_data.max_basis_size;
1203 basis_size, this->memory);
1205 basis_size, this->memory);
1209 unsigned int accumulated_iterations = 0;
1212 H.reinit(basis_size + 1, basis_size);
1215 Vector<double> projected_rhs;
1226 aux->sadd(-1., 1., b);
1228 double beta = aux->l2_norm();
1230 iteration_state = this->iteration_status(accumulated_iterations, res, x);
1234 H.reinit(basis_size + 1, basis_size);
1237 for (
unsigned int j = 0; j < basis_size; ++j)
1240 v(j, x).equ(1. / a, *aux);
1245 preconditioner.vmult(z(j, x), v[j]);
1246 A.vmult(*aux, z[j]);
1249 H(0, j) = *aux * v[0];
1250 for (
unsigned int i = 1; i <= j; ++i)
1251 H(i, j) = aux->add_and_dot(-H(i - 1, j), v[i - 1], v[i]);
1252 H(j + 1, j) = a = std::sqrt(aux->add_and_dot(-H(j, j), v[j], *aux));
1258 H1.reinit(j + 1, j);
1259 projected_rhs.reinit(j + 1);
1261 projected_rhs(0) = beta;
1271 this->iteration_status(++accumulated_iterations, res, x);
1278 for (
unsigned int j = 0; j < y.size(); ++j)
PETScWrappers::SolverGMRES SolverGMRES
SolverFGMRES(SolverControl &cn, VectorMemory< VectorType > &mem, const AdditionalData &data=AdditionalData())
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)
boost::signals2::signal< void(const internal::SolverGMRESImplementation::TmpVectors< VectorType > &)> krylov_space_signal
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)
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)
AdditionalData additional_data
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)
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)
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
#define DEAL_II_NAMESPACE_CLOSE
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
Expression fabs(const Expression &x)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
number singular_value(const size_type i) const
VectorMemory< VectorType > & mem
virtual double criterion()
unsigned int size() const
std::array< NumberType, 3 > givens_rotation(const NumberType &x, const NumberType &y)
#define DEAL_II_NAMESPACE_OPEN
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)
bool use_default_residual
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)
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_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
T max(const T &t, const MPI_Comm &mpi_communicator)
unsigned int max_n_tmp_vectors
static ::ExceptionBase & ExcInternalError()