From a115d2ac2e970779422c612678be08b126724b59 Mon Sep 17 00:00:00 2001 From: Matthias Maier Date: Mon, 4 May 2015 17:26:09 +0200 Subject: [PATCH] use intermediate storage if vmult is called with same dest and src --- include/deal.II/lac/linear_operator.h | 168 ++++++++++++++++++-------- 1 file changed, 115 insertions(+), 53 deletions(-) diff --git a/include/deal.II/lac/linear_operator.h b/include/deal.II/lac/linear_operator.h index f6fc3f2f57..f0bc25b8f3 100644 --- a/include/deal.II/lac/linear_operator.h +++ b/include/deal.II/lac/linear_operator.h @@ -781,10 +781,32 @@ namespace }; - // A helper class to add the full matrix interface to a LinearOperator + // A helper function to apply a given vmult, or Tvmult to a vector with + // intermediate storage + template + void apply_with_intermediate_storage(Function function, Range &v, + const Domain &u, bool add) { + static GrowingVectorMemory vector_memory; + + Range *i = vector_memory.alloc(); + i->reinit(v, /*bool fast =*/true); + + function(*i, u); + + if (add) + v += *i; + else + v = *i; + + vector_memory.free(i); + } + + + // A helper class to add a reduced matrix interface to a LinearOperator + // (typically provided by Preconditioner classes) template - class MatrixInterfaceWithVmultAdd + class MatrixInterfaceWithoutVmultAdd { public: template @@ -792,67 +814,107 @@ namespace { op.vmult = [&matrix](Range &v, const Domain &u) { - matrix.vmult(v,u); + if (static_cast(&v) == static_cast(&u)) + { + // If v and u are the same memory location use intermediate storage + apply_with_intermediate_storage([&matrix](Range &b, const Domain &a) + { + matrix.vmult(b, a); + }, + v, u, /*bool add =*/ false); + } + else + { + matrix.vmult(v,u); + } }; op.vmult_add = [&matrix](Range &v, const Domain &u) { - matrix.vmult_add(v,u); + // use intermediate storage to implement vmult_add with vmult + apply_with_intermediate_storage([&matrix](Range &b, const Domain &a) + { + matrix.vmult(b, a); + }, + v, u, /*bool add =*/ true); }; op.Tvmult = [&matrix](Domain &v, const Range &u) { - matrix.Tvmult(v,u); + if (static_cast(&v) == static_cast(&u)) + { + // If v and u are the same memory location use intermediate storage + apply_with_intermediate_storage([&matrix](Domain &b, const Range &a) + { + matrix.Tvmult(b, a); + }, + v, u, /*bool add =*/ false); + } + else + { + matrix.Tvmult(v,u); + } + }; op.Tvmult_add = [&matrix](Domain &v, const Range &u) { - matrix.Tvmult_add(v,u); + // use intermediate storage to implement Tvmult_add with Tvmult + apply_with_intermediate_storage([&matrix](Domain &b, const Range &a) + { + matrix.Tvmult(b, a); + }, + v, u, /*bool add =*/ true); }; } }; - // A helper class to add a reduced matrix interface to a LinearOperator - // (typically provided by Preconditioner classes) + // A helper class to add the full matrix interface to a LinearOperator template - class MatrixInterfaceWithoutVmultAdd + class MatrixInterfaceWithVmultAdd { public: template void operator()(LinearOperator &op, const Matrix &matrix) { - op.vmult = [&matrix](Range &v, const Domain &u) - { - matrix.vmult(v,u); - }; + // As above ... - op.vmult_add = [op](Range &v, const Domain &u) - { - static GrowingVectorMemory vector_memory; + MatrixInterfaceWithoutVmultAdd().operator()(op, matrix); - Range *i = vector_memory.alloc(); - op.reinit_range_vector(*i, /*bool fast =*/true); - op.vmult(*i, u); - v += *i; - vector_memory.free(i); - }; + // ... but add native vmult_add and Tvmult_add variants: - op.Tvmult = [&matrix](Domain &v, const Range &u) + op.vmult_add = [&matrix](Range &v, const Domain &u) { - matrix.Tvmult(v,u); + if (static_cast(&v) == static_cast(&u)) + { + apply_with_intermediate_storage([&matrix](Range &b, const Domain &a) + { + matrix.vmult(b, a); + }, + v, u, /*bool add =*/ false); + } + else + { + matrix.vmult_add(v,u); + } }; - op.Tvmult_add = [op](Domain &v, const Range &u) + op.Tvmult_add = [&matrix](Domain &v, const Range &u) { - static GrowingVectorMemory vector_memory; - - Domain *i = vector_memory.alloc(); - op.reinit_domain_vector(*i, /*bool fast =*/true); - op.Tvmult(*i, u); - v += *i; - vector_memory.free(i); + if (static_cast(&v) == static_cast(&u)) + { + apply_with_intermediate_storage([&matrix](Domain &b, const Range &a) + { + matrix.Tvmult(b, a); + }, + v, u, /*bool add =*/ true); + } + else + { + matrix.Tvmult_add(v,u); + } }; } }; @@ -2091,11 +2153,11 @@ operator*(const PackagedOperation &comp, * Create a PackagedOperation object from a reference to a matrix @p m and * a reference to a vector @p u of the Domain space. The object stores the * PackagedOperation $\text{m} \,u$ (in matrix notation). If a different - * result type than Vector is desired, the matrix should be + * result type than Domain is desired, the matrix should be * explicitly wrapper with linear_operator instead. * * return (return_add) are implemented with - * vmult(__1,u) (vmult_add(__1,u)). + * vmult(_1,u) (vmult_add(_1,u)). * * The PackagedOperation object that is created stores references to @p u * and @p m. Thus, both must remain valid references for the whole lifetime @@ -2108,13 +2170,13 @@ operator*(const PackagedOperation &comp, template ::type::value>::type, - typename = typename std::enable_if, Domain>, Matrix>::value>::type, - typename = typename std::enable_if, Domain, Matrix>::type::value>::type> -PackagedOperation > + typename = typename std::enable_if, Matrix>::value>::type, + typename = typename std::enable_if::type::value>::type> +PackagedOperation operator*(const Matrix &m, const Domain &u) { - return linear_operator, Domain>(m) * u; + return linear_operator(m) * u; } @@ -2128,13 +2190,13 @@ operator*(const Matrix &m, template ::type::value>::type, - typename = typename std::enable_if, Domain>, Matrix>::value>::type, - typename = typename std::enable_if, Domain, Matrix>::type::value>::type> -PackagedOperation > + typename = typename std::enable_if, Matrix>::value>::type, + typename = typename std::enable_if::type::value>::type> +PackagedOperation operator*(const Matrix &m, const PackagedOperation &u) { - return linear_operator, Domain>(m) * u; + return linear_operator(m) * u; } @@ -2142,13 +2204,13 @@ operator*(const Matrix &m, * @relates PackagedOperation * * Create a PackagedOperation object from a reference to a matrix @p m and - * a reference to a vector @p u of the Domain space. The object stores the + * a reference to a vector @p u of the Range space. The object stores the * PackagedOperation $\text{m}^T \,u$ (in matrix notation). If a different - * result type than Vector is desired, the matrix should be + * result type than Range is desired, the matrix should be * explicitly wrapper with linear_operator instead. * * return (return_add) are implemented with - * Tvmult(__1,u) (Tvmult_add(__1,u)). + * Tvmult(_1,u) (Tvmult_add(_1,u)). * * The PackagedOperation object that is created stores references to @p u * and @p m. Thus, both must remain valid references for the whole lifetime @@ -2161,13 +2223,13 @@ operator*(const Matrix &m, template ::type::value>::type, - typename = typename std::enable_if >, Matrix>::value>::type, - typename = typename std::enable_if, Matrix>::type::value>::type> -PackagedOperation > + typename = typename std::enable_if, Matrix>::value>::type, + typename = typename std::enable_if::type::value>::type> +PackagedOperation operator*(const Range &u, const Matrix &m) { - return u * linear_operator>(m); + return u * linear_operator(m); } @@ -2181,13 +2243,13 @@ operator*(const Range &u, template ::type::value>::type, - typename = typename std::enable_if >, Matrix>::value>::type, - typename = typename std::enable_if, Matrix>::type::value>::type> -PackagedOperation > + typename = typename std::enable_if, Matrix>::value>::type, + typename = typename std::enable_if::type::value>::type> +PackagedOperation operator*(const PackagedOperation &u, const Matrix &m) { - return u * linear_operator>(m); + return u * linear_operator(m); } //@} -- 2.39.5