15#ifndef dealii_solver_cg_h
16#define dealii_solver_cg_h
41 template <
typename,
typename>
175template <
typename VectorType = Vector<
double>>
213 template <
typename MatrixType,
typename PreconditionerType>
218 const PreconditionerType &preconditioner);
226 boost::signals2::connection
228 const std::function<
void(
typename VectorType::value_type,
229 typename VectorType::value_type)> &slot);
237 boost::signals2::connection
239 const bool every_iteration =
false);
247 boost::signals2::connection
249 const std::function<
void(
const std::vector<double> &)> &slot,
250 const bool every_iteration =
false);
262 const VectorType &d)
const;
271 const std::vector<typename VectorType::value_type> &diagonal,
272 const std::vector<typename VectorType::value_type> &offdiagonal,
273 const boost::signals2::signal<
void(
const std::vector<double> &)>
275 const boost::signals2::signal<
void(
double)> &cond_signal);
285 boost::signals2::signal<void(
typename VectorType::value_type,
286 typename VectorType::value_type)>
311 boost::signals2::signal<void(
const std::vector<double> &)>
352template <
typename VectorType = Vector<
double>>
393template <
typename VectorType>
396 const AdditionalData &data)
398 , additional_data(data)
399 , determine_beta_by_flexible_formula(false)
404template <
typename VectorType>
407 , additional_data(data)
408 , determine_beta_by_flexible_formula(false)
413template <
typename VectorType>
418 const VectorType &)
const
423template <
typename VectorType>
426 const std::vector<typename VectorType::value_type> &diagonal,
427 const std::vector<typename VectorType::value_type> &offdiagonal,
428 const boost::signals2::signal<
void(
const std::vector<double> &)>
430 const boost::signals2::signal<
void(
double)> &cond_signal)
433 if (!cond_signal.empty() || !eigenvalues_signal.empty())
437 for (size_type i = 0; i <
diagonal.size(); ++i)
441 T(i, i + 1) = offdiagonal[i];
443 T.compute_eigenvalues();
447 auto condition_number = T.eigenvalue(T.n() - 1) / T.eigenvalue(0);
450 cond_signal(
std::abs(condition_number));
454 if (!eigenvalues_signal.empty())
457 for (
unsigned int j = 0; j < T.n(); ++j)
480 template <
typename VectorType,
482 typename PreconditionerType>
483 struct IterationWorkerBase
485 using Number =
typename VectorType::value_type;
488 const PreconditionerType &preconditioner;
508 Number r_dot_preconditioner_dot_r;
511 double residual_norm;
512 Number previous_alpha;
514 IterationWorkerBase(
const MatrixType &A,
515 const PreconditionerType &preconditioner,
520 , preconditioner(preconditioner)
531 , r_dot_preconditioner_dot_r(Number())
535 , previous_alpha(Number())
539 startup(
const VectorType &b)
559 residual_norm = r.l2_norm();
567 template <
typename VectorType,
569 typename PreconditionerType,
571 struct IterationWorker
572 :
public IterationWorkerBase<VectorType, MatrixType, PreconditionerType>
575 IterationWorkerBase<VectorType, MatrixType, PreconditionerType>;
577 IterationWorker(
const MatrixType &A,
578 const PreconditionerType &preconditioner,
582 : BaseClass(A, preconditioner, flexible, memory, x)
586 using BaseClass::alpha;
587 using BaseClass::beta;
589 using BaseClass::preconditioner;
591 using BaseClass::r_dot_preconditioner_dot_r;
592 using BaseClass::residual_norm;
598 do_iteration(
const unsigned int iteration_index)
600 using Number =
typename VectorType::value_type;
602 const Number previous_r_dot_preconditioner_dot_r =
603 r_dot_preconditioner_dot_r;
605 if (std::is_same_v<PreconditionerType, PreconditionIdentity> ==
false)
607 preconditioner.vmult(v, r);
608 r_dot_preconditioner_dot_r = r * v;
611 r_dot_preconditioner_dot_r = residual_norm * residual_norm;
613 const VectorType &direction =
614 std::is_same_v<PreconditionerType, PreconditionIdentity> ? r : v;
616 if (iteration_index > 1)
621 r_dot_preconditioner_dot_r / previous_r_dot_preconditioner_dot_r;
623 beta -= (r * z) / previous_r_dot_preconditioner_dot_r;
624 p.sadd(beta, 1., direction);
627 p.equ(1., direction);
634 const Number p_dot_A_dot_p =
p * v;
637 this->previous_alpha = alpha;
638 alpha = r_dot_preconditioner_dot_r / p_dot_A_dot_p;
645 finalize_after_convergence(
const unsigned int)
657 template <
typename MatrixType,
typename VectorType>
658 using vmult_functions_t =
decltype(std::declval<const MatrixType>().vmult(
659 std::declval<VectorType &>(),
660 std::declval<const VectorType &>(),
662 const std::function<
void(
const unsigned int,
const unsigned int)> &>(),
663 std::declval<
const std::function<
void(
const unsigned int,
664 const unsigned int)> &>()));
666 template <
typename MatrixType,
typename VectorType>
667 constexpr bool has_vmult_functions =
668 is_supported_operation<vmult_functions_t, MatrixType, VectorType>;
673 template <
typename PreconditionerType>
674 using apply_to_subrange_t =
675 decltype(std::declval<const PreconditionerType>()
676 .apply_to_subrange(0U, 0U,
nullptr,
nullptr));
678 template <
typename PreconditionerType>
679 constexpr bool has_apply_to_subrange =
680 is_supported_operation<apply_to_subrange_t, PreconditionerType>;
685 template <
typename PreconditionerType>
687 decltype(std::declval<const PreconditionerType>().apply(0U, 0.0));
689 template <
typename PreconditionerType>
690 constexpr bool has_apply =
691 is_supported_operation<apply_t, PreconditionerType>;
697 template <
typename VectorType,
699 typename PreconditionerType>
700 struct IterationWorker<
704 std::enable_if_t<has_vmult_functions<MatrixType, VectorType> &&
705 (has_apply_to_subrange<PreconditionerType> ||
706 has_apply<PreconditionerType>)&&std::
707 is_same_v<VectorType,
708 LinearAlgebra::distributed::Vector<
709 typename VectorType::value_type,
712 :
public IterationWorkerBase<VectorType, MatrixType, PreconditionerType>
714 using Number =
typename VectorType::value_type;
716 Number next_r_dot_preconditioner_dot_r;
717 Number previous_beta;
719 IterationWorker(
const MatrixType &A,
720 const PreconditionerType &preconditioner,
724 : IterationWorkerBase<VectorType, MatrixType, PreconditionerType>(
730 , next_r_dot_preconditioner_dot_r(0.)
737 do_iteration(
const unsigned int iteration_index)
739 if (iteration_index > 1)
741 previous_beta = this->beta;
742 this->beta = next_r_dot_preconditioner_dot_r /
743 this->r_dot_preconditioner_dot_r;
746 std::array<VectorizedArray<Number>, 7> vectorized_sums = {};
751 [&](
const unsigned int begin,
const unsigned int end) {
752 operation_before_loop(iteration_index, begin, end);
754 [&](
const unsigned int begin,
const unsigned int end) {
755 operation_after_loop(begin, end, vectorized_sums);
758 std::array<Number, 7> scalar_sums;
759 for (
unsigned int i = 0; i < 7; ++i)
760 scalar_sums[i] = vectorized_sums[i][0];
761 for (
unsigned int l = 1; l < VectorizedArray<Number>::size(); ++
l)
762 for (
unsigned int i = 0; i < 7; ++i)
763 scalar_sums[i] += vectorized_sums[i][l];
767 this->r.get_mpi_communicator(),
770 this->r_dot_preconditioner_dot_r = scalar_sums[6];
772 const Number p_dot_A_dot_p = scalar_sums[0];
775 this->previous_alpha = this->alpha;
776 this->alpha = this->r_dot_preconditioner_dot_r / p_dot_A_dot_p;
782 this->alpha * (-2. * scalar_sums[2] + this->alpha * scalar_sums[1])));
784 next_r_dot_preconditioner_dot_r =
std::abs(
785 this->r_dot_preconditioner_dot_r +
786 this->alpha * (-2. * scalar_sums[4] + this->alpha * scalar_sums[5]));
791 template <
typename U =
void>
792 std::enable_if_t<has_apply<PreconditionerType>, U>
793 operation_before_loop(
const unsigned int iteration_index,
794 const unsigned int start_range,
795 const unsigned int end_range)
const
797 Number *x = this->x.begin();
798 Number *r = this->r.begin();
799 Number *
p = this->p.begin();
800 Number *v = this->v.begin();
801 const Number alpha = this->alpha;
802 const Number beta = this->beta;
804 const unsigned int end_regular =
805 start_range + (end_range - start_range) / n_lanes * n_lanes;
806 if (iteration_index == 1)
811 for (
unsigned int j = start_range; j < end_regular; j += n_lanes)
816 for (
unsigned int l = 0;
l < n_lanes; ++
l)
817 pj[l] = this->preconditioner.apply(j + l, rj[l]);
822 for (
unsigned int j = end_regular; j < end_range; ++j)
824 p[j] = this->preconditioner.apply(j, r[j]);
828 else if (iteration_index % 2 == 0 && beta != Number())
830 for (
unsigned int j = start_range; j < end_regular; j += n_lanes)
838 for (
unsigned int l = 0;
l < n_lanes; ++
l)
839 prec_rj[l] = this->preconditioner.apply(j + l, rj[l]);
841 pj = beta * pj + prec_rj;
846 for (
unsigned int j = end_regular; j < end_range; ++j)
848 r[j] -= alpha * v[j];
849 p[j] = beta *
p[j] + this->preconditioner.apply(j, r[j]);
853 else if (iteration_index % 2 == 0 && beta == Number())
858 for (
unsigned int j = start_range; j < end_regular; j += n_lanes)
870 for (
unsigned int l = 0;
l < n_lanes; ++
l)
871 prec_rj[l] = this->preconditioner.apply(j + l, rj[l]);
872 prec_rj.
store(p + j);
876 for (
unsigned int j = end_regular; j < end_range; ++j)
878 r[j] -= alpha * v[j];
879 x[j] += alpha *
p[j];
880 p[j] = this->preconditioner.apply(j, r[j]);
886 const Number alpha_plus_previous_alpha_over_beta =
887 alpha + this->previous_alpha / this->previous_beta;
888 const Number previous_alpha_over_beta =
889 this->previous_alpha / this->previous_beta;
890 for (
unsigned int j = start_range; j < end_regular; j += n_lanes)
895 xj += alpha_plus_previous_alpha_over_beta * pj;
899 for (
unsigned int l = 0;
l < n_lanes; ++
l)
901 prec_rj[
l] = this->preconditioner.apply(j + l, rj[l]);
902 prec_vj[
l] = this->preconditioner.apply(j + l, vj[l]);
904 xj -= previous_alpha_over_beta * prec_rj;
908 prec_rj -= alpha * prec_vj;
909 pj = beta * pj + prec_rj;
914 for (
unsigned int j = end_regular; j < end_range; ++j)
916 x[j] += alpha_plus_previous_alpha_over_beta *
p[j];
917 x[j] -= previous_alpha_over_beta *
918 this->preconditioner.apply(j, r[j]);
919 r[j] -= alpha * v[j];
920 p[j] = beta *
p[j] + this->preconditioner.apply(j, r[j]);
928 template <
typename U =
void>
929 std::enable_if_t<has_apply<PreconditionerType>, U>
930 operation_after_loop(
931 const unsigned int start_range,
932 const unsigned int end_range,
935 const Number *r = this->r.begin();
936 const Number *
p = this->p.begin();
937 const Number *v = this->v.begin();
938 std::array<VectorizedArray<Number>, 7> my_sums = {};
940 const unsigned int end_regular =
941 start_range + (end_range - start_range) / n_lanes * n_lanes;
942 for (
unsigned int j = start_range; j < end_regular; j += n_lanes)
949 for (
unsigned int l = 0;
l < n_lanes; ++
l)
951 prec_vj[
l] = this->preconditioner.apply(j + l, vj[l]);
952 prec_rj[
l] = this->preconditioner.apply(j + l, rj[l]);
954 my_sums[0] += pj * vj;
955 my_sums[1] += vj * vj;
956 my_sums[2] += rj * vj;
957 my_sums[3] += rj * rj;
958 my_sums[4] += rj * prec_vj;
959 my_sums[5] += vj * prec_vj;
960 my_sums[6] += rj * prec_rj;
962 for (
unsigned int j = end_regular; j < end_range; ++j)
964 const Number prec_v = this->preconditioner.apply(j, v[j]);
965 const Number prec_r = this->preconditioner.apply(j, r[j]);
966 my_sums[0][0] +=
p[j] * v[j];
967 my_sums[1][0] += v[j] * v[j];
968 my_sums[2][0] += r[j] * v[j];
969 my_sums[3][0] += r[j] * r[j];
970 my_sums[4][0] += r[j] * prec_v;
971 my_sums[5][0] += v[j] * prec_v;
972 my_sums[6][0] += r[j] * prec_r;
974 for (
unsigned int i = 0; i < vectorized_sums.size(); ++i)
975 vectorized_sums[i] += my_sums[i];
980 template <
typename U =
void>
981 std::enable_if_t<has_apply<PreconditionerType>, U>
982 finalize_after_convergence(
const unsigned int iteration_index)
984 if (iteration_index % 2 == 1 || this->beta == Number())
985 this->x.add(this->alpha, this->p);
988 using Number =
typename VectorType::value_type;
989 const unsigned int end_range = this->x.locally_owned_size();
991 Number *x = this->x.begin();
992 Number *r = this->r.begin();
993 Number *
p = this->p.begin();
999 const Number alpha_plus_previous_alpha_over_beta =
1000 this->alpha + this->previous_alpha / this->beta;
1001 const Number previous_alpha_over_beta =
1002 this->previous_alpha / this->beta;
1005 for (
unsigned int j = 0; j < end_range; ++j)
1007 x[j] += alpha_plus_previous_alpha_over_beta *
p[j] -
1008 previous_alpha_over_beta *
1009 this->preconditioner.apply(j, r[j]);
1017 template <
typename U =
void>
1018 std::enable_if_t<!has_apply<PreconditionerType>, U>
1019 operation_before_loop(
const unsigned int iteration_index,
1020 const unsigned int start_range,
1021 const unsigned int end_range)
const
1023 Number *x = this->x.begin() + start_range;
1024 Number *r = this->r.begin() + start_range;
1025 Number *
p = this->p.begin() + start_range;
1026 Number *v = this->v.begin() + start_range;
1027 const Number alpha = this->alpha;
1028 const Number beta = this->beta;
1029 constexpr unsigned int grain_size = 128;
1030 std::array<Number, grain_size> prec_r;
1031 if (iteration_index == 1)
1033 for (
unsigned int j = start_range; j < end_range; j += grain_size)
1035 const unsigned int length =
std::min(grain_size, end_range - j);
1036 this->preconditioner.apply_to_subrange(j,
1041 for (
unsigned int i = 0; i < length; ++i)
1051 else if (iteration_index % 2 == 0 && beta != Number())
1053 for (
unsigned int j = start_range; j < end_range; j += grain_size)
1055 const unsigned int length =
std::min(grain_size, end_range - j);
1057 for (
unsigned int i = 0; i < length; ++i)
1058 r[i] -= this->alpha * v[i];
1059 this->preconditioner.apply_to_subrange(j,
1064 for (
unsigned int i = 0; i < length; ++i)
1066 p[i] = this->beta *
p[i] + prec_r[i];
1074 else if (iteration_index % 2 == 0 && beta == Number())
1079 for (
unsigned int j = start_range; j < end_range; j += grain_size)
1081 const unsigned int length =
std::min(grain_size, end_range - j);
1083 for (
unsigned int i = 0; i < length; ++i)
1084 r[i] -= this->alpha * v[i];
1085 this->preconditioner.apply_to_subrange(j,
1090 for (
unsigned int i = 0; i < length; ++i)
1092 x[i] += this->alpha *
p[i];
1103 const Number alpha_plus_previous_alpha_over_beta =
1104 this->alpha + this->previous_alpha / this->previous_beta;
1105 const Number previous_alpha_over_beta =
1106 this->previous_alpha / this->previous_beta;
1107 for (
unsigned int j = start_range; j < end_range; j += grain_size)
1109 const unsigned int length =
std::min(grain_size, end_range - j);
1110 this->preconditioner.apply_to_subrange(j,
1115 for (
unsigned int i = 0; i < length; ++i)
1117 x[i] += alpha_plus_previous_alpha_over_beta *
p[i] -
1118 previous_alpha_over_beta * prec_r[i];
1119 r[i] -= this->alpha * v[i];
1121 this->preconditioner.apply_to_subrange(j,
1126 for (
unsigned int i = 0; i < length; ++i)
1128 p[i] = this->beta *
p[i] + prec_r[i];
1142 template <
typename U =
void>
1143 std::enable_if_t<!has_apply<PreconditionerType>, U>
1144 operation_after_loop(
1145 const unsigned int start_range,
1146 const unsigned int end_range,
1149 const Number *r = this->r.begin();
1150 const Number *
p = this->p.begin();
1151 const Number *v = this->v.begin();
1152 std::array<VectorizedArray<Number>, 7> my_sums = {};
1153 constexpr unsigned int grain_size = 128;
1156 const unsigned int end_regular =
1157 start_range + (end_range - start_range) / grain_size * grain_size;
1158 std::array<Number, grain_size> prec_r;
1159 std::array<Number, grain_size> prec_v;
1160 for (
unsigned int j = start_range; j < end_regular; j += grain_size)
1162 this->preconditioner.apply_to_subrange(j,
1166 this->preconditioner.apply_to_subrange(j,
1171 for (
unsigned int i = 0; i < grain_size;
1177 prec_rj.
load(prec_r.data() + i);
1178 prec_vj.
load(prec_v.data() + i);
1180 my_sums[0] += pj * vj;
1181 my_sums[1] += vj * vj;
1182 my_sums[2] += rj * vj;
1183 my_sums[3] += rj * rj;
1184 my_sums[4] += rj * prec_vj;
1185 my_sums[5] += vj * prec_vj;
1186 my_sums[6] += rj * prec_rj;
1189 const unsigned int length = end_range - end_regular;
1191 this->preconditioner.apply_to_subrange(end_regular,
1192 end_regular + length,
1195 this->preconditioner.apply_to_subrange(end_regular,
1196 end_regular + length,
1199 for (
unsigned int j = end_regular; j < end_range; ++j)
1201 my_sums[0][0] +=
p[j] * v[j];
1202 my_sums[1][0] += v[j] * v[j];
1203 my_sums[2][0] += r[j] * v[j];
1204 my_sums[3][0] += r[j] * r[j];
1205 my_sums[4][0] += r[j] * prec_v[j - end_regular];
1206 my_sums[5][0] += v[j] * prec_v[j - end_regular];
1207 my_sums[6][0] += r[j] * prec_r[j - end_regular];
1209 for (
unsigned int i = 0; i < vectorized_sums.size(); ++i)
1210 vectorized_sums[i] += my_sums[i];
1216 template <
typename U =
void>
1217 std::enable_if_t<!has_apply<PreconditionerType>, U>
1218 finalize_after_convergence(
const unsigned int iteration_index)
1220 if (iteration_index % 2 == 1 || this->beta == Number())
1221 this->x.add(this->alpha, this->p);
1224 const unsigned int end_range = this->x.locally_owned_size();
1226 Number *x = this->x.begin();
1227 Number *r = this->r.begin();
1228 Number *
p = this->p.begin();
1229 const Number alpha_plus_previous_alpha_over_beta =
1230 this->alpha + this->previous_alpha / this->beta;
1231 const Number previous_alpha_over_beta =
1232 this->previous_alpha / this->beta;
1234 constexpr unsigned int grain_size = 128;
1235 std::array<Number, grain_size> prec_r;
1236 for (
unsigned int j = 0; j < end_range; j += grain_size)
1238 const unsigned int length =
std::min(grain_size, end_range - j);
1239 this->preconditioner.apply_to_subrange(j,
1244 for (
unsigned int i = 0; i < length; ++i)
1245 x[i] += alpha_plus_previous_alpha_over_beta * p[i] -
1246 previous_alpha_over_beta * prec_r[i];
1260template <
typename VectorType>
1261template <
typename MatrixType,
typename PreconditionerType>
1265 const VectorType &b,
1266 const PreconditionerType &preconditioner)
1268 using number =
typename VectorType::value_type;
1275 const bool do_eigenvalues =
1276 !condition_number_signal.empty() || !all_condition_numbers_signal.empty() ||
1277 !eigenvalues_signal.empty() || !all_eigenvalues_signal.empty();
1280 std::vector<typename VectorType::value_type>
diagonal;
1281 std::vector<typename VectorType::value_type> offdiagonal;
1283 typename VectorType::value_type eigen_beta_alpha = 0;
1287 internal::SolverCG::
1288 IterationWorker<VectorType, MatrixType, PreconditionerType>
1290 A, preconditioner, determine_beta_by_flexible_formula, this->memory, x);
1294 solver_state = this->iteration_status(0, worker.residual_norm, x);
1302 worker.do_iteration(it);
1304 print_vectors(it, x, worker.r, worker.p);
1308 this->coefficients_signal(worker.previous_alpha, worker.beta);
1313 diagonal.push_back(number(1.) / worker.previous_alpha +
1315 eigen_beta_alpha = worker.beta / worker.previous_alpha;
1316 offdiagonal.push_back(
std::sqrt(worker.beta) /
1317 worker.previous_alpha);
1319 compute_eigs_and_cond(diagonal,
1321 all_eigenvalues_signal,
1322 all_condition_numbers_signal);
1325 solver_state = this->iteration_status(it, worker.residual_norm, x);
1328 worker.finalize_after_convergence(it);
1330 compute_eigs_and_cond(diagonal,
1333 condition_number_signal);
1341template <
typename VectorType>
1342boost::signals2::connection
1344 const std::function<
void(
typename VectorType::value_type,
1345 typename VectorType::value_type)> &slot)
1347 return coefficients_signal.connect(slot);
1352template <
typename VectorType>
1353boost::signals2::connection
1355 const std::function<
void(
double)> &slot,
1356 const bool every_iteration)
1358 if (every_iteration)
1360 return all_condition_numbers_signal.connect(slot);
1364 return condition_number_signal.connect(slot);
1370template <
typename VectorType>
1371boost::signals2::connection
1373 const std::function<
void(
const std::vector<double> &)> &slot,
1374 const bool every_iteration)
1376 if (every_iteration)
1378 return all_eigenvalues_signal.connect(slot);
1382 return eigenvalues_signal.connect(slot);
1388template <
typename VectorType>
1391 const AdditionalData &)
1399template <
typename VectorType>
1401 const AdditionalData &)
SolverCG(SolverControl &cn, VectorMemory< VectorType > &mem, const AdditionalData &data=AdditionalData())
boost::signals2::signal< void(const std::vector< double > &)> eigenvalues_signal
bool determine_beta_by_flexible_formula
virtual void print_vectors(const unsigned int step, const VectorType &x, const VectorType &r, const VectorType &d) const
boost::signals2::signal< void(typename VectorType::value_type, typename VectorType::value_type)> coefficients_signal
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::signal< void(double)> condition_number_signal
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)
boost::signals2::signal< void(double)> all_condition_numbers_signal
SolverCG(SolverControl &cn, const AdditionalData &data=AdditionalData())
boost::signals2::signal< void(const std::vector< double > &)> all_eigenvalues_signal
AdditionalData additional_data
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.
SolverFlexibleCG(SolverControl &cn, const AdditionalData &data=AdditionalData())
SolverFlexibleCG(SolverControl &cn, VectorMemory< VectorType > &mem, const AdditionalData &data=AdditionalData())
static constexpr std::size_t size()
void store(OtherNumber *ptr) const
void load(const OtherNumber *ptr)
#define DEAL_II_OPENMP_SIMD_PRAGMA
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcDivideByZero()
#define AssertThrow(cond, exc)
@ diagonal
Matrix is diagonal.
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
T sum(const T &t, const MPI_Comm mpi_communicator)
const Iterator const std_cxx20::type_identity_t< Iterator > & end
::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
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)