From dd7636d15ae96911ee010a8d0eabb4121b3e7883 Mon Sep 17 00:00:00 2001 From: ESeNonFossiIo Date: Wed, 3 Jun 2015 16:09:55 +0200 Subject: [PATCH] is_null_operator --- include/deal.II/lac/linear_operator.h | 265 +++++++++++++++----------- 1 file changed, 154 insertions(+), 111 deletions(-) diff --git a/include/deal.II/lac/linear_operator.h b/include/deal.II/lac/linear_operator.h index a30e1bda1d..71b44f5db1 100644 --- a/include/deal.II/lac/linear_operator.h +++ b/include/deal.II/lac/linear_operator.h @@ -46,6 +46,16 @@ template , typename Matrix> LinearOperator linear_operator (const Matrix &); +template , + typename Range = Domain> +LinearOperator +null_operator(const std::function &, + const std::function &); + +template > +LinearOperator +null_operator(const std::function &); + /** * A class to store the abstract concept of a linear operator. @@ -110,8 +120,11 @@ public: * member objects are initialized with default variants that throw an * exception upon invocation. */ - LinearOperator () + LinearOperator (bool is_null_operator=false) + : + is_null_operator(is_null_operator) { + vmult = [](Range &, const Domain &) { Assert(false, ExcMessage("Uninitialized LinearOperator operator+(const LinearOperator &first_op, const LinearOperator &second_op) { - LinearOperator return_op; + if (first_op.is_null_operator) + { + return second_op; + } + else if (second_op.is_null_operator) + { + return first_op; + } + else + { + LinearOperator return_op; - return_op.reinit_range_vector = first_op.reinit_range_vector; - return_op.reinit_domain_vector = first_op.reinit_domain_vector; + return_op.reinit_range_vector = first_op.reinit_range_vector; + return_op.reinit_domain_vector = first_op.reinit_domain_vector; - // ensure to have valid computation objects by catching first_op and - // second_op by value + // ensure to have valid computation objects by catching first_op and + // second_op by value - return_op.vmult = [first_op, second_op](Range &v, const Domain &u) - { - first_op.vmult(v, u); - second_op.vmult_add(v, u); - }; + return_op.vmult = [first_op, second_op](Range &v, const Domain &u) + { + first_op.vmult(v, u); + second_op.vmult_add(v, u); + }; - return_op.vmult_add = [first_op, second_op](Range &v, const Domain &u) - { - first_op.vmult_add(v, u); - second_op.vmult_add(v, u); - }; + return_op.vmult_add = [first_op, second_op](Range &v, const Domain &u) + { + first_op.vmult_add(v, u); + second_op.vmult_add(v, u); + }; - return_op.Tvmult = [first_op, second_op](Domain &v, const Range &u) - { - second_op.Tvmult(v, u); - first_op.Tvmult_add(v, u); - }; + return_op.Tvmult = [first_op, second_op](Domain &v, const Range &u) + { + second_op.Tvmult(v, u); + first_op.Tvmult_add(v, u); + }; - return_op.Tvmult_add = [first_op, second_op](Domain &v, const Range &u) - { - second_op.Tvmult_add(v, u); - first_op.Tvmult_add(v, u); - }; + return_op.Tvmult_add = [first_op, second_op](Domain &v, const Range &u) + { + second_op.Tvmult_add(v, u); + first_op.Tvmult_add(v, u); + }; - return return_op; + return return_op; + } } @@ -346,8 +377,19 @@ LinearOperator operator-(const LinearOperator &first_op, const LinearOperator &second_op) { - // implement with addition and scalar multiplication - return first_op + (-1. * second_op); + if (first_op.is_null_operator) + { + return -1 * second_op; + } + else if (second_op.is_null_operator) + { + return first_op; + } + else + { + // implement with addition and scalar multiplication + return first_op + (-1. * second_op); + } } @@ -377,58 +419,51 @@ operator*(typename Range::value_type number, std::is_convertible::value, "Range and Domain must have implicitly convertible 'value_type's"); - LinearOperator return_op = op; + if (op.is_null_operator) + { + return op; + } + else if (number == 0.) + { + LinearOperator return_op = op; + return_op.is_null_operator = true; + return return_op; + } + else + { + LinearOperator return_op = op; - // ensure to have valid computation objects by catching number and op by - // value + // ensure to have valid computation objects by catching number and op by + // value - return_op.vmult = [number, op](Range &v, const Domain &u) - { - if (number == 0.) - { - v = 0.; - } - else + return_op.vmult = [number, op](Range &v, const Domain &u) { op.vmult(v,u); v *= number; - } - }; + }; - return_op.vmult_add = [number, op](Range &v, const Domain &u) - { - if (number != 0.) + return_op.vmult_add = [number, op](Range &v, const Domain &u) { v /= number; op.vmult_add(v,u); v *= number; - } - }; + }; - return_op.Tvmult = [number, op](Domain &v, const Range &u) - { - if (number == 0.) - { - v = 0.; - } - else + return_op.Tvmult = [number, op](Domain &v, const Range &u) { op.Tvmult(v,u); v *= number; - } - }; + }; - return_op.Tvmult_add = [number, op](Domain &v, const Range &u) - { - if (number != 0.) + return_op.Tvmult_add = [number, op](Domain &v, const Range &u) { v /= number; op.Tvmult_add(v,u); v *= number; - } - }; + }; - return return_op; + return return_op; + } } @@ -480,59 +515,66 @@ LinearOperator operator*(const LinearOperator &first_op, const LinearOperator &second_op) { - LinearOperator return_op; + if (first_op.is_null_operator || second_op.is_null_operator) + { + return null_operator(second_op.reinit_domain_vector, first_op.reinit_range_vector); + } + else + { + LinearOperator return_op; - return_op.reinit_domain_vector = second_op.reinit_domain_vector; - return_op.reinit_range_vector = first_op.reinit_range_vector; + return_op.reinit_domain_vector = second_op.reinit_domain_vector; + return_op.reinit_range_vector = first_op.reinit_range_vector; - // ensure to have valid computation objects by catching first_op and - // second_op by value + // ensure to have valid computation objects by catching first_op and + // second_op by value - return_op.vmult = [first_op, second_op](Range &v, const Domain &u) - { - static GrowingVectorMemory vector_memory; + return_op.vmult = [first_op, second_op](Range &v, const Domain &u) + { + static GrowingVectorMemory vector_memory; - Intermediate *i = vector_memory.alloc(); - second_op.reinit_range_vector(*i, /*bool fast =*/ true); - second_op.vmult(*i, u); - first_op.vmult(v, *i); - vector_memory.free(i); - }; + Intermediate *i = vector_memory.alloc(); + second_op.reinit_range_vector(*i, /*bool fast =*/ true); + second_op.vmult(*i, u); + first_op.vmult(v, *i); + vector_memory.free(i); + }; - return_op.vmult_add = [first_op, second_op](Range &v, const Domain &u) - { - static GrowingVectorMemory vector_memory; + return_op.vmult_add = [first_op, second_op](Range &v, const Domain &u) + { + static GrowingVectorMemory vector_memory; - Intermediate *i = vector_memory.alloc(); - second_op.reinit_range_vector(*i, /*bool fast =*/ true); - second_op.vmult(*i, u); - first_op.vmult_add(v, *i); - vector_memory.free(i); - }; + Intermediate *i = vector_memory.alloc(); + second_op.reinit_range_vector(*i, /*bool fast =*/ true); + second_op.vmult(*i, u); + first_op.vmult_add(v, *i); + vector_memory.free(i); + }; - return_op.Tvmult = [first_op, second_op](Domain &v, const Range &u) - { - static GrowingVectorMemory vector_memory; + return_op.Tvmult = [first_op, second_op](Domain &v, const Range &u) + { + static GrowingVectorMemory vector_memory; - Intermediate *i = vector_memory.alloc(); - first_op.reinit_domain_vector(*i, /*bool fast =*/ true); - first_op.Tvmult(*i, u); - second_op.Tvmult(v, *i); - vector_memory.free(i); - }; + Intermediate *i = vector_memory.alloc(); + first_op.reinit_domain_vector(*i, /*bool fast =*/ true); + first_op.Tvmult(*i, u); + second_op.Tvmult(v, *i); + vector_memory.free(i); + }; - return_op.Tvmult_add = [first_op, second_op](Domain &v, const Range &u) - { - static GrowingVectorMemory vector_memory; + return_op.Tvmult_add = [first_op, second_op](Domain &v, const Range &u) + { + static GrowingVectorMemory vector_memory; - Intermediate *i = vector_memory.alloc(); - first_op.reinit_domain_vector(*i, /*bool fast =*/ true); - first_op.Tvmult(*i, u); - second_op.Tvmult_add(v, *i); - vector_memory.free(i); - }; + Intermediate *i = vector_memory.alloc(); + first_op.reinit_domain_vector(*i, /*bool fast =*/ true); + first_op.Tvmult(*i, u); + second_op.Tvmult_add(v, *i); + vector_memory.free(i); + }; - return return_op; + return return_op; + } } /** @@ -695,31 +737,32 @@ identity_operator(const std::function &reinit_vector) * * @ingroup LAOperators */ + template + typename Range> LinearOperator null_operator(const std::function &reinit_domain_vector, - const std::function &reinit_range_vector) + const std::function &reinit_range_vector) { - LinearOperator return_op; + LinearOperator return_op(/*is_null_operator = */true); + return_op.reinit_range_vector = reinit_range_vector; return_op.reinit_domain_vector = reinit_domain_vector; - return_op.reinit_range_vector = reinit_range_vector; - return_op.vmult = [](Range &v, const Range &u) + return_op.vmult = [](Range &v, const Domain &u) { v = 0; }; - return_op.vmult_add = [](Range &v, const Range &u) + return_op.vmult_add = [](Range &v, const Domain &u) {}; - return_op.Tvmult = [](Range &v, const Range &u) + return_op.Tvmult = [](Range &v, const Domain &u) { v = 0; }; - return_op.Tvmult_add = [](Range &v, const Range &u) + return_op.Tvmult_add = [](Range &v, const Domain &u) {}; return return_op; -- 2.39.5