16 #ifndef dealii_linear_operator_h
17 #define dealii_linear_operator_h
27 #include <type_traits>
35 namespace LinearOperatorImplementation
41 template <
typename Number>
46 template <
typename Range = Vector<
double>,
47 typename Domain = Range,
49 internal::LinearOperatorImplementation::EmptyPayload>
55 typename Domain = Range,
57 typename OperatorExemplar,
64 typename Domain = Range,
72 typename Domain = Range,
77 template <
typename Range,
typename Domain,
typename Payload>
200 template <
typename Range,
typename Domain,
typename Payload>
216 vmult = [](Range &,
const Domain &) {
218 ExcMessage(
"Uninitialized LinearOperator<Range, "
219 "Domain>::vmult called"));
222 vmult_add = [](Range &,
const Domain &) {
224 ExcMessage(
"Uninitialized LinearOperator<Range, "
225 "Domain>::vmult_add called"));
228 Tvmult = [](Domain &,
const Range &) {
230 ExcMessage(
"Uninitialized LinearOperator<Range, "
231 "Domain>::Tvmult called"));
236 ExcMessage(
"Uninitialized LinearOperator<Range, "
237 "Domain>::Tvmult_add called"));
242 ExcMessage(
"Uninitialized LinearOperator<Range, "
243 "Domain>::reinit_range_vector method called"));
248 ExcMessage(
"Uninitialized LinearOperator<Range, "
249 "Domain>::reinit_domain_vector method called"));
263 template <
typename Op,
264 typename = std::enable_if_t<
265 !std::is_base_of_v<LinearOperator<Range, Domain, Payload>, Op>>>
268 *
this = linear_operator<Range, Domain, Payload, Op>(op);
281 template <
typename Op,
282 typename = std::enable_if_t<
283 !std::is_base_of_v<LinearOperator<Range, Domain, Payload>, Op>>>
287 *
this = linear_operator<Range, Domain, Payload, Op>(op);
295 std::function<void(Range &v,
const Domain &u)>
vmult;
301 std::function<void(Range &v,
const Domain &u)>
vmult_add;
307 std::function<void(Domain &v,
const Range &u)>
Tvmult;
331 std::function<void(Domain &v,
bool omit_zeroing_entries)>
346 *
this = *
this + second_op;
357 *
this = *
this - second_op;
368 *
this = *
this * second_op;
379 *
this = *
this * number;
408 template <
typename Range,
typename Domain,
typename Payload>
424 static_cast<const Payload &
>(first_op) +
425 static_cast<const Payload &
>(second_op)};
433 return_op.vmult = [first_op, second_op](Range &v,
const Domain &u) {
434 first_op.
vmult(v, u);
438 return_op.vmult_add = [first_op, second_op](Range &v,
const Domain &u) {
443 return_op.Tvmult = [first_op, second_op](Domain &v,
const Range &u) {
448 return_op.Tvmult_add = [first_op, second_op](Domain &v,
const Range &u) {
467 template <
typename Range,
typename Domain,
typename Payload>
474 return -1. * second_op;
483 return first_op + (-1. * second_op);
505 template <
typename Range,
typename Domain,
typename Payload>
511 std::is_convertible<
typename Range::value_type,
512 typename Domain::value_type>::value,
513 "Range and Domain must have implicitly convertible 'value_type's");
519 else if (number == 0.)
530 return_op.
vmult = [number, op](Range &v,
const Domain &u) {
535 return_op.
vmult_add = [number, op](Range &v,
const Domain &u) {
541 return_op.
Tvmult = [number, op](Domain &v,
const Range &u) {
546 return_op.
Tvmult_add = [number, op](Domain &v,
const Range &u) {
572 template <
typename Range,
typename Domain,
typename Payload>
575 typename Domain::value_type number)
578 std::is_convertible<
typename Range::value_type,
579 typename Domain::value_type>::value,
580 "Range and Domain must have implicitly convertible 'value_type's");
602 template <
typename Range,
603 typename Intermediate,
620 static_cast<const Payload &
>(first_op) *
621 static_cast<const Payload &
>(second_op)};
629 return_op.vmult = [first_op, second_op](Range &v,
const Domain &u) {
634 second_op.
vmult(*i, u);
635 first_op.
vmult(v, *i);
638 return_op.vmult_add = [first_op, second_op](Range &v,
const Domain &u) {
643 second_op.
vmult(*i, u);
647 return_op.Tvmult = [first_op, second_op](Domain &v,
const Range &u) {
656 return_op.Tvmult_add = [first_op, second_op](Domain &v,
const Range &u) {
678 template <
typename Range,
typename Domain,
typename Payload>
687 return_op.vmult = op.
Tvmult;
689 return_op.Tvmult = op.
vmult;
717 typename Preconditioner,
718 typename Range =
typename Solver::vector_type,
719 typename Domain = Range>
723 const Preconditioner &preconditioner)
726 op.inverse_payload(solver, preconditioner)};
731 return_op.vmult = [op, &solver, &preconditioner](Range &v,
const Domain &u) {
733 solver.solve(op, v, u, preconditioner);
736 return_op.vmult_add = [op, &solver, &preconditioner](Range &v,
742 solver.solve(op, *v2, u, preconditioner);
746 return_op.Tvmult = [op, &solver, &preconditioner](Range &v,
const Domain &u) {
751 return_op.Tvmult_add = [op, &solver, &preconditioner](Range &v,
775 typename Range =
typename Solver::vector_type,
776 typename Domain = Range>
783 op.inverse_payload(solver, preconditioner)};
788 return_op.vmult = [op, &solver, preconditioner](Range &v,
const Domain &u) {
790 solver.solve(op, v, u, preconditioner);
793 return_op.vmult_add = [op, &solver, preconditioner](Range &v,
799 solver.solve(op, *v2, u, preconditioner);
803 return_op.Tvmult = [op, &solver, preconditioner](Range &v,
const Domain &u) {
808 return_op.Tvmult_add = [op, &solver, preconditioner](Range &v,
833 typename Range =
typename Solver::vector_type,
834 typename Domain = Range>
853 typename Range =
typename Solver::vector_type,
854 typename Domain = Range>
890 return_op.reinit_range_vector = reinit_vector;
891 return_op.reinit_domain_vector = reinit_vector;
893 return_op.vmult = [](Range &v,
const Range &u) { v = u; };
895 return_op.vmult_add = [](Range &v,
const Range &u) { v += u; };
897 return_op.Tvmult = [](Range &v,
const Range &u) { v = u; };
899 return_op.Tvmult_add = [](Range &v,
const Range &u) { v += u; };
917 template <
typename Range,
typename Domain,
typename Payload>
922 static_cast<Payload &
>(return_op) = op.identity_payload();
937 template <
typename Range,
typename Domain,
typename Payload>
948 return_op.vmult = [](Range &v,
const Domain &) { v = 0.; };
950 return_op.vmult_add = [](Range &,
const Domain &) {};
952 return_op.Tvmult = [](Domain &v,
const Range &) { v = 0.; };
954 return_op.Tvmult_add = [](Domain &,
const Range &) {};
980 return_op.reinit_range_vector = reinit_vector;
981 return_op.reinit_domain_vector = reinit_vector;
983 return_op.vmult = [](Range &v,
const Range &u) {
984 const auto mean = u.mean_value();
990 return_op.vmult_add = [](Range &v,
const Range &u) {
991 const auto mean = u.mean_value();
997 return_op.Tvmult = return_op.vmult_add;
998 return_op.Tvmult_add = return_op.vmult_add;
1016 template <
typename Range,
typename Domain,
typename Payload>
1021 static_cast<Payload &
>(return_op) = op.identity_payload();
1029 namespace LinearOperatorImplementation
1042 template <
typename Vector>
1057 template <
typename Matrix>
1061 bool omit_zeroing_entries)
1077 template <
typename Matrix>
1081 bool omit_zeroing_entries)
1109 template <
typename... Args>
1147 template <
typename Solver,
typename Preconditioner>
1179 template <
typename Range,
typename Domain,
typename T>
1182 template <
typename C>
1183 static std::false_type
1186 template <
typename C>
1189 -> decltype(std::declval<C>().vmult_add(*r, *
d),
1190 std::declval<C>().Tvmult_add(*
d, *r),
1197 using type = decltype(test<T>(
nullptr,
nullptr));
1203 template <
typename Function,
typename Range,
typename Domain>
1226 template <
typename Range,
typename Domain,
typename Payload>
1230 template <
typename Matrix>
1291 template <
typename Range,
typename Domain,
typename Payload>
1295 template <
typename Matrix>
1398 template <
typename Range,
typename Domain,
typename Payload,
typename Matrix>
1403 return linear_operator<Range, Domain, Payload, Matrix, Matrix>(
matrix,
1422 template <
typename Range,
1425 typename OperatorExemplar,
1440 return_op.reinit_range_vector =
1441 [&operator_exemplar](Range &v,
bool omit_zeroing_entries) {
1443 Range>::reinit_range_vector(operator_exemplar, v, omit_zeroing_entries);
1446 return_op.reinit_domain_vector = [&operator_exemplar](
1447 Domain &v,
bool omit_zeroing_entries) {
1449 Domain>::reinit_domain_vector(operator_exemplar, v, omit_zeroing_entries);
1457 operator()(return_op,
matrix);
1480 template <
typename Range,
typename Domain,
typename Payload,
typename Matrix>
1487 auto return_op = operator_exemplar;
1494 operator()(return_op,
matrix);
1511 typename Domain = Range,
1513 typename OperatorExemplar,
1515 typename = std::enable_if_t<!std::is_lvalue_reference_v<Matrix>>>
1521 typename Domain = Range,
1523 typename OperatorExemplar,
1525 typename = std::enable_if_t<!std::is_lvalue_reference_v<OperatorExemplar>>,
1526 typename = std::enable_if_t<
1527 !std::is_same_v<OperatorExemplar, LinearOperator<Range, Domain, Payload>>>>
1533 typename Domain = Range,
1535 typename OperatorExemplar,
1537 typename = std::enable_if_t<!std::is_lvalue_reference_v<Matrix>>,
1538 typename = std::enable_if_t<!std::is_lvalue_reference_v<OperatorExemplar>>,
1539 typename = std::enable_if_t<
1540 !std::is_same_v<OperatorExemplar, LinearOperator<Range, Domain, Payload>>>>
1546 typename Domain = Range,
1549 typename = std::enable_if_t<!std::is_lvalue_reference_v<Matrix>>>
1552 Matrix &&) =
delete;
1556 typename Domain = Range,
1559 typename = std::enable_if_t<!std::is_lvalue_reference_v<Matrix>>>
1566 typename Preconditioner,
1567 typename Range =
typename Solver::vector_type,
1568 typename Domain = Range,
1569 typename = std::enable_if_t<!std::is_lvalue_reference_v<Preconditioner>>,
1571 std::enable_if_t<!std::is_same_v<Preconditioner, PreconditionIdentity>>,
1572 typename = std::enable_if_t<
1573 !std::is_same_v<Preconditioner, LinearOperator<Range, Domain, Payload>>>>
1577 Preconditioner &&) =
delete;
LinearOperator< Range, Domain, Payload > & operator*=(const LinearOperator< Domain, Domain, Payload > &second_op)
LinearOperator< Range, Domain, Payload > operator*=(typename Domain::value_type number)
LinearOperator< Range, Domain, Payload > & operator-=(const LinearOperator< Range, Domain, Payload > &second_op)
std::function< void(Range &v, const Domain &u)> vmult_add
std::function< void(Domain &v, const Range &u)> Tvmult
std::function< void(Domain &v, bool omit_zeroing_entries)> reinit_domain_vector
LinearOperator(const Payload &payload=Payload())
std::function< void(Range &v, const Domain &u)> vmult
LinearOperator< Range, Domain, Payload > & operator+=(const LinearOperator< Range, Domain, Payload > &second_op)
std::function< void(Range &v, bool omit_zeroing_entries)> reinit_range_vector
std::function< void(Domain &v, const Range &u)> Tvmult_add
LinearOperator(const Op &op)
LinearOperator< Range, Domain, Payload > & operator=(const Op &op)
LinearOperator< Range, Domain, Payload > & operator=(const LinearOperator< Range, Domain, Payload > &)=default
LinearOperator(const LinearOperator< Range, Domain, Payload > &)=default
virtual void reinit(const size_type N, const bool omit_zeroing_entries=false)
EmptyPayload(const Args &...)
EmptyPayload transpose_payload() const
EmptyPayload identity_payload() const
EmptyPayload null_payload() const
EmptyPayload inverse_payload(Solver &, const Preconditioner &) const
void operator()(LinearOperator< Range, Domain, Payload > &op, const Matrix &matrix)
void operator()(LinearOperator< Range, Domain, Payload > &op, const Matrix &matrix)
static void reinit_domain_vector(const Matrix &matrix, Vector &v, bool omit_zeroing_entries)
static void reinit_range_vector(const Matrix &matrix, Vector &v, bool omit_zeroing_entries)
static std::false_type test(...)
static auto test(Range *r, Domain *d) -> decltype(std::declval< C >().vmult_add(*r, *d), std::declval< C >().Tvmult_add(*d, *r), std::true_type())
decltype(test< T >(nullptr, nullptr)) type
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
#define Assert(cond, exc)
static ::ExceptionBase & ExcMessage(std::string arg1)
LinearOperator< Domain, Range, Payload > inverse_operator(const LinearOperator< Range, Domain, Payload > &op, Solver &solver)
LinearOperator< Range, Domain, Payload > operator*(const LinearOperator< Range, Intermediate, Payload > &first_op, const LinearOperator< Intermediate, Domain, Payload > &second_op)
LinearOperator< Range, Domain, Payload > linear_operator(const OperatorExemplar &operator_exemplar, const Matrix &matrix)
LinearOperator< Range, Domain, Payload > operator-(const LinearOperator< Range, Domain, Payload > &first_op, const LinearOperator< Range, Domain, Payload > &second_op)
LinearOperator< Range, Range, Payload > identity_operator(const std::function< void(Range &, bool)> &reinit_vector)
LinearOperator< Range, Domain, Payload > null_operator(const LinearOperator< Range, Domain, Payload > &op)
LinearOperator< Range, Domain, Payload > identity_operator(const LinearOperator< Range, Domain, Payload > &op)
LinearOperator< Domain, Range, Payload > inverse_operator(const LinearOperator< Range, Domain, Payload > &op, Solver &solver, const PreconditionIdentity &)
LinearOperator< Domain, Range, Payload > inverse_operator(const LinearOperator< Range, Domain, Payload > &op, Solver &solver, const LinearOperator< Range, Domain, Payload > &preconditioner)
LinearOperator< Range, Domain, Payload > linear_operator(const Matrix &matrix)
LinearOperator< Domain, Range, Payload > transpose_operator(const LinearOperator< Range, Domain, Payload > &op)
LinearOperator< Range, Domain, Payload > mean_value_filter(const LinearOperator< Range, Domain, Payload > &op)
LinearOperator< Range, Domain, Payload > operator*(typename Range::value_type number, const LinearOperator< Range, Domain, Payload > &op)
LinearOperator< Domain, Range, Payload > inverse_operator(const LinearOperator< Range, Domain, Payload > &op, Solver &solver, const Preconditioner &preconditioner)
LinearOperator< Range, Domain, Payload > operator*(const LinearOperator< Range, Domain, Payload > &op, typename Domain::value_type number)
LinearOperator< Range, Domain, Payload > linear_operator(const LinearOperator< Range, Domain, Payload > &operator_exemplar, const Matrix &matrix)
LinearOperator< Range, Range, Payload > mean_value_filter(const std::function< void(Range &, bool)> &reinit_vector)
LinearOperator< Range, Domain, Payload > operator+(const LinearOperator< Range, Domain, Payload > &first_op, const LinearOperator< Range, Domain, Payload > &second_op)
@ matrix
Contents is actually a matrix.
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
EmptyPayload operator+(const EmptyPayload &, const EmptyPayload &)
void apply_with_intermediate_storage(Function function, Range &v, const Domain &u, bool add)
EmptyPayload operator*(const EmptyPayload &, const EmptyPayload &)
static bool equal(const T *p1, const T *p2)