15#ifndef dealii_solver_cg_h
16#define dealii_solver_cg_h
42 template <
typename,
typename>
194template <
typename VectorType = Vector<
double>>
254 template <
typename MatrixType,
typename PreconditionerType>
261 const PreconditionerType &preconditioner);
270 connect_coefficients_slot(
291 connect_eigenvalues_slot(
302 print_vectors(
const unsigned int step,
313 compute_eigs_and_cond(
355 all_eigenvalues_signal;
365 bool determine_beta_by_flexible_formula;
463template <
typename VectorType>
466 const bool use_default_residual)
467 : use_default_residual(use_default_residual)
470template <
typename VectorType>
474 const AdditionalData &
data)
476 , additional_data(
data)
477 , determine_beta_by_flexible_formula(
false)
482template <
typename VectorType>
486 , additional_data(
data)
487 , determine_beta_by_flexible_formula(
false)
492template <
typename VectorType>
497 const VectorType &)
const
502template <
typename VectorType>
505 const std::vector<typename VectorType::value_type> &diagonal,
506 const std::vector<typename VectorType::value_type> &
offdiagonal,
507 const boost::signals2::signal<
void(
const std::vector<double> &)>
509 const boost::signals2::signal<
void(
double)> &
cond_signal)
512 if (!
cond_signal.empty() || !eigenvalues_signal.empty())
516 for (size_type i = 0; i <
diagonal.size(); ++i)
522 T.compute_eigenvalues();
533 if (!eigenvalues_signal.empty())
536 for (
unsigned int j = 0;
j < T.n(); ++
j)
562 typename PreconditionerType>
565 using Number =
typename VectorType::value_type;
568 const PreconditionerType &preconditioner;
596 const bool use_default_residual;
600 const PreconditionerType &preconditioner,
605 const bool &use_default_residual)
607 , preconditioner(preconditioner)
626 , use_default_residual(use_default_residual)
639 if (!use_default_residual)
662 typename PreconditionerType,
672 const PreconditionerType &preconditioner,
677 const bool &use_default_residual)
684 use_default_residual)
688 using BaseClass::alpha;
690 using BaseClass::beta;
691 using BaseClass::explicit_r;
693 using BaseClass::preconditioner;
695 using BaseClass::r_dot_preconditioner_dot_r;
696 using BaseClass::residual_norm;
697 using BaseClass::use_default_residual;
705 using Number =
typename VectorType::value_type;
710 if (std::is_same_v<PreconditionerType, PreconditionIdentity> ==
false)
712 preconditioner.vmult(v,
r);
719 std::is_same_v<PreconditionerType, PreconditionIdentity> ?
r : v;
729 p.sadd(
beta, 1., direction);
732 p.equ(1., direction);
742 this->previous_alpha =
alpha;
748 if (use_default_residual)
778 template <
typename MatrixType,
typename VectorType>
780 std::declval<VectorType &>(),
781 std::declval<const VectorType &>(),
783 const std::function<
void(
const unsigned int,
const unsigned int)> &>(),
784 std::declval<
const std::function<
void(
const unsigned int,
785 const unsigned int)> &>()));
787 template <
typename MatrixType,
typename VectorType>
794 template <
typename PreconditionerType>
796 decltype(std::declval<const PreconditionerType>()
797 .apply_to_subrange(0U, 0U,
nullptr,
nullptr));
799 template <
typename PreconditionerType>
806 template <
typename PreconditionerType>
808 decltype(std::declval<const PreconditionerType>().apply(0U, 0.0));
810 template <
typename PreconditionerType>
820 typename PreconditionerType>
826 (has_apply_to_subrange<PreconditionerType> ||
827 has_apply<PreconditionerType>)&&std::
828 is_same_v<VectorType,
829 LinearAlgebra::distributed::Vector<
830 typename VectorType::value_type,
835 using Number =
typename VectorType::value_type;
841 const PreconditionerType &preconditioner,
846 const bool &use_default_residual)
854 use_default_residual)
868 this->r_dot_preconditioner_dot_r;
876 [&](
const unsigned int begin,
const unsigned int end) {
879 [&](
const unsigned int begin,
const unsigned int end) {
884 for (
unsigned int i = 0; i < 7; ++i)
887 for (
unsigned int i = 0; i < 7; ++i)
892 this->r.get_mpi_communicator(),
900 this->previous_alpha = this->alpha;
901 this->alpha = this->r_dot_preconditioner_dot_r /
p_dot_A_dot_p;
910 this->r_dot_preconditioner_dot_r +
916 template <
typename U =
void>
917 std::enable_if_t<has_apply<PreconditionerType>, U>
922 Number *
x = this->x.begin();
923 Number *
r = this->r.begin();
924 Number *p = this->p.begin();
925 Number *v = this->v.begin();
926 const Number
alpha = this->alpha;
927 const Number
beta = this->beta;
941 for (
unsigned int l = 0;
l < n_lanes; ++
l)
942 pj[l] = this->preconditioner.apply(
j + l,
rj[l]);
949 p[
j] = this->preconditioner.apply(
j,
r[
j]);
963 for (
unsigned int l = 0;
l < n_lanes; ++
l)
964 prec_rj[l] = this->preconditioner.apply(
j + l,
rj[l]);
974 p[
j] =
beta * p[
j] + this->preconditioner.apply(
j,
r[
j]);
995 for (
unsigned int l = 0;
l < n_lanes; ++
l)
996 prec_rj[l] = this->preconditioner.apply(
j + l,
rj[l]);
1005 p[
j] = this->preconditioner.apply(
j,
r[
j]);
1012 alpha + this->previous_alpha / this->previous_beta;
1014 this->previous_alpha / this->previous_beta;
1024 for (
unsigned int l = 0;
l < n_lanes; ++
l)
1026 prec_rj[
l] = this->preconditioner.apply(
j + l,
rj[l]);
1027 prec_vj[
l] = this->preconditioner.apply(
j + l,
vj[l]);
1043 this->preconditioner.apply(
j,
r[
j]);
1045 p[
j] =
beta * p[
j] + this->preconditioner.apply(
j,
r[
j]);
1053 template <
typename U =
void>
1054 std::enable_if_t<has_apply<PreconditionerType>, U>
1060 const Number *
r = this->r.begin();
1061 const Number *p = this->p.begin();
1062 const Number *v = this->v.begin();
1063 std::array<VectorizedArray<Number>, 7>
my_sums = {};
1074 for (
unsigned int l = 0;
l < n_lanes; ++
l)
1076 prec_vj[
l] = this->preconditioner.apply(
j + l,
vj[l]);
1077 prec_rj[
l] = this->preconditioner.apply(
j + l,
rj[l]);
1089 const Number
prec_v = this->preconditioner.apply(
j, v[
j]);
1090 const Number
prec_r = this->preconditioner.apply(
j,
r[
j]);
1105 template <
typename U =
void>
1106 std::enable_if_t<has_apply<PreconditionerType>, U>
1110 this->x.add(this->alpha, this->p);
1113 using Number =
typename VectorType::value_type;
1114 const unsigned int end_range = this->x.locally_owned_size();
1116 Number *
x = this->x.begin();
1117 Number *
r = this->r.begin();
1118 Number *p = this->p.begin();
1125 this->alpha + this->previous_alpha / this->beta;
1127 this->previous_alpha / this->beta;
1134 this->preconditioner.apply(
j,
r[
j]);
1142 template <
typename U =
void>
1143 std::enable_if_t<!has_apply<PreconditionerType>, U>
1152 const Number
alpha = this->alpha;
1153 const Number
beta = this->beta;
1155 std::array<Number, grain_size>
prec_r;
1161 this->preconditioner.apply_to_subrange(
j,
1166 for (
unsigned int i = 0; i < length; ++i)
1182 for (
unsigned int i = 0; i < length; ++i)
1183 r[i] -= this->alpha * v[i];
1184 this->preconditioner.apply_to_subrange(
j,
1189 for (
unsigned int i = 0; i < length; ++i)
1191 p[i] = this->beta * p[i] +
prec_r[i];
1208 for (
unsigned int i = 0; i < length; ++i)
1209 r[i] -= this->alpha * v[i];
1210 this->preconditioner.apply_to_subrange(
j,
1215 for (
unsigned int i = 0; i < length; ++i)
1217 x[i] += this->alpha * p[i];
1229 this->alpha + this->previous_alpha / this->previous_beta;
1231 this->previous_alpha / this->previous_beta;
1235 this->preconditioner.apply_to_subrange(
j,
1240 for (
unsigned int i = 0; i < length; ++i)
1244 r[i] -= this->alpha * v[i];
1246 this->preconditioner.apply_to_subrange(
j,
1251 for (
unsigned int i = 0; i < length; ++i)
1253 p[i] = this->beta * p[i] +
prec_r[i];
1267 template <
typename U =
void>
1268 std::enable_if_t<!has_apply<PreconditionerType>, U>
1274 const Number *
r = this->r.begin();
1275 const Number *p = this->p.begin();
1276 const Number *v = this->v.begin();
1277 std::array<VectorizedArray<Number>, 7>
my_sums = {};
1283 std::array<Number, grain_size>
prec_r;
1284 std::array<Number, grain_size>
prec_v;
1287 this->preconditioner.apply_to_subrange(
j,
1291 this->preconditioner.apply_to_subrange(
j,
1316 this->preconditioner.apply_to_subrange(
end_regular,
1320 this->preconditioner.apply_to_subrange(
end_regular,
1341 template <
typename U =
void>
1342 std::enable_if_t<!has_apply<PreconditionerType>, U>
1346 this->x.add(this->alpha, this->p);
1349 const unsigned int end_range = this->x.locally_owned_size();
1351 Number *
x = this->x.begin();
1352 Number *
r = this->r.begin();
1353 Number *p = this->p.begin();
1355 this->alpha + this->previous_alpha / this->beta;
1357 this->previous_alpha / this->beta;
1360 std::array<Number, grain_size>
prec_r;
1364 this->preconditioner.apply_to_subrange(
j,
1369 for (
unsigned int i = 0; i < length; ++i)
1385template <
typename VectorType>
1387template <
typename MatrixType,
typename PreconditionerType>
1393 const VectorType &b,
1394 const PreconditionerType &preconditioner)
1396 using number =
typename VectorType::value_type;
1404 !condition_number_signal.empty() || !all_condition_numbers_signal.empty() ||
1405 !eigenvalues_signal.empty() || !all_eigenvalues_signal.empty();
1408 std::vector<typename VectorType::value_type>
diagonal;
1409 std::vector<typename VectorType::value_type>
offdiagonal;
1415 internal::SolverCG::
1416 IterationWorker<VectorType, MatrixType, PreconditionerType>
1419 determine_beta_by_flexible_formula,
1427 solver_state = this->iteration_status(0, worker.residual_norm,
x);
1435 worker.do_iteration(
it);
1437 print_vectors(
it,
x, worker.r, worker.p);
1441 this->coefficients_signal(worker.previous_alpha, worker.beta);
1446 diagonal.push_back(number(1.) / worker.previous_alpha +
1450 worker.previous_alpha);
1452 compute_eigs_and_cond(diagonal,
1454 all_eigenvalues_signal,
1455 all_condition_numbers_signal);
1461 worker.finalize_after_convergence(
it);
1463 compute_eigs_and_cond(diagonal,
1466 condition_number_signal);
1474template <
typename VectorType>
1477 const std::function<
void(
typename VectorType::value_type,
1478 typename VectorType::value_type)> &
slot)
1480 return coefficients_signal.connect(
slot);
1485template <
typename VectorType>
1488 const std::function<
void(
double)> &
slot,
1493 return all_condition_numbers_signal.connect(
slot);
1497 return condition_number_signal.connect(
slot);
1503template <
typename VectorType>
1506 const std::function<
void(
const std::vector<double> &)> &
slot,
1511 return all_eigenvalues_signal.connect(
slot);
1515 return eigenvalues_signal.connect(
slot);
1521template <
typename VectorType>
1524 const bool use_default_residual)
1525 : use_default_residual(use_default_residual)
1529template <
typename VectorType>
1533 const AdditionalData &
data)
1536 this->determine_beta_by_flexible_formula =
true;
1537 this->additional_data =
data;
1542template <
typename VectorType>
1545 const AdditionalData &
data)
1548 this->determine_beta_by_flexible_formula =
true;
1549 this->additional_data =
data;
SolverCG(SolverControl &cn, VectorMemory< VectorType > &mem, const AdditionalData &data=AdditionalData())
virtual void print_vectors(const unsigned int step, const VectorType &x, const VectorType &r, const VectorType &d) const
boost::signals2::connection connect_condition_number_slot(const std::function< void(double)> &slot, const bool every_iteration=false)
boost::signals2::connection connect_eigenvalues_slot(const std::function< void(const std::vector< double > &)> &slot, const bool every_iteration=false)
boost::signals2::connection connect_coefficients_slot(const std::function< void(typename VectorType::value_type, typename VectorType::value_type)> &slot)
virtual ~SolverCG() override=default
void solve(const MatrixType &A, VectorType &x, const VectorType &b, const PreconditionerType &preconditioner)
SolverCG(SolverControl &cn, const AdditionalData &data=AdditionalData())
static void compute_eigs_and_cond(const std::vector< typename VectorType::value_type > &diagonal, const std::vector< typename VectorType::value_type > &offdiagonal, const boost::signals2::signal< void(const std::vector< double > &)> &eigenvalues_signal, const boost::signals2::signal< void(double)> &cond_signal)
@ iterate
Continue iteration.
@ success
Stop iteration, goal reached.
AdditionalData additional_data
SolverFlexibleCG(SolverControl &cn, const AdditionalData &data=AdditionalData())
SolverFlexibleCG(SolverControl &cn, VectorMemory< VectorType > &mem, const AdditionalData &data=AdditionalData())
static constexpr std::size_t size()
#define DEAL_II_OPENMP_SIMD_PRAGMA
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_CXX20_REQUIRES(condition)
#define DEAL_II_NAMESPACE_CLOSE
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcDivideByZero()
#define AssertThrow(cond, exc)
std::vector< index_type > data
@ diagonal
Matrix is diagonal.
Tpetra::Vector< Number, LO, GO, NodeType< MemorySpace > > VectorType
Tpetra::CrsMatrix< Number, LO, GO, NodeType< MemorySpace > > MatrixType
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
VectorType::value_type * end(VectorType &V)
VectorType::value_type * begin(VectorType &V)
T sum(const T &t, const MPI_Comm mpi_communicator)
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
unsigned int global_dof_index
AdditionalData(const bool use_default_residual=true)
bool use_default_residual
bool use_default_residual
AdditionalData(const bool use_default_residual=true)
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)