From 66956aa05d9cd5d75ff934ad1bcc31dc51cf6ead Mon Sep 17 00:00:00 2001 From: Daniel Arndt Date: Fri, 2 Aug 2019 14:02:45 -0400 Subject: [PATCH] Reduce floating-point division in operator/ --- include/deal.II/base/point.h | 49 ++++++++- include/deal.II/base/tensor.h | 101 ++++++++++++++++-- include/deal.II/lac/block_vector_base.h | 4 +- .../data_out_base_vtu_04.with_zlib=on.output | 2 +- tests/tensors/tensor_11.cc | 49 +++++++++ tests/tensors/tensor_11.output | 7 ++ 6 files changed, 199 insertions(+), 13 deletions(-) create mode 100644 tests/tensors/tensor_11.cc create mode 100644 tests/tensors/tensor_11.output diff --git a/include/deal.II/base/point.h b/include/deal.II/base/point.h index f2836107ec..8a62844aba 100644 --- a/include/deal.II/base/point.h +++ b/include/deal.II/base/point.h @@ -560,6 +560,50 @@ inline DEAL_II_CUDA_HOST_DEV +namespace internal +{ + namespace PointImplementation + { + template ::type>::value, + int>::type = 0> + inline DEAL_II_CUDA_HOST_DEV + Point::type> + division_operator(const Point &p, const OtherNumber factor) + { + using common_type = typename ProductType::type; + Point tmp; + const auto inverse_factor = + static_cast(1.) / static_cast(factor); + for (unsigned int i = 0; i < dim; ++i) + tmp[i] = p[i] * inverse_factor; + return tmp; + } + + template ::type>::value, + int>::type = 0> + inline DEAL_II_CUDA_HOST_DEV + Point::type> + division_operator(const Point &p, const OtherNumber factor) + { + Point::type> tmp; + for (unsigned int i = 0; i < dim; ++i) + tmp[i] = p[i] / factor; + return tmp; + } + } // namespace PointImplementation +} // namespace internal + + template template inline DEAL_II_CUDA_HOST_DEV @@ -568,10 +612,7 @@ inline DEAL_II_CUDA_HOST_DEV typename EnableIfScalar::type>::type> Point::operator/(const OtherNumber factor) const { - Point::type> tmp; - for (unsigned int i = 0; i < dim; ++i) - tmp[i] = this->operator[](i) / factor; - return tmp; + return internal::PointImplementation::division_operator(*this, factor); } diff --git a/include/deal.II/base/tensor.h b/include/deal.II/base/tensor.h index c2fd168fcb..e7326c427e 100644 --- a/include/deal.II/base/tensor.h +++ b/include/deal.II/base/tensor.h @@ -1368,14 +1368,56 @@ DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE } +namespace internal +{ + namespace TensorImplementation + { + template ::type>::value, + int>::type = 0> + DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV inline DEAL_II_ALWAYS_INLINE void + division_operator(Tensor (&t)[dim], + const OtherNumber &factor) + { + const Number inverse_factor = Number(1.) / factor; + // recurse over the base objects + for (unsigned int d = 0; d < dim; ++d) + t[d] *= inverse_factor; + } + + + template ::type>::value, + int>::type = 0> + DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV inline DEAL_II_ALWAYS_INLINE void + division_operator(dealii::Tensor (&t)[dim], + const OtherNumber &factor) + { + // recurse over the base objects + for (unsigned int d = 0; d < dim; ++d) + t[d] /= factor; + } + } // namespace TensorImplementation +} // namespace internal + + template template DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE DEAL_II_CUDA_HOST_DEV Tensor & Tensor::operator/=(const OtherNumber &s) { - for (unsigned int i = 0; i < dim; ++i) - values[i] /= s; + internal::TensorImplementation::division_operator(values, s); return *this; } @@ -1758,6 +1800,55 @@ DEAL_II_CUDA_HOST_DEV constexpr DEAL_II_ALWAYS_INLINE } +namespace internal +{ + namespace TensorImplementation + { + template ::type>::value, + int>::type = 0> + DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV inline DEAL_II_ALWAYS_INLINE + Tensor::type> + division_operator(const Tensor &t, + const OtherNumber & factor) + { + Tensor::type> tt; + const Number inverse_factor = Number(1.) / factor; + // recurse over the base objects + for (unsigned int d = 0; d < dim; ++d) + tt[d] = t[d] * inverse_factor; + return tt; + } + + + template ::type>::value, + int>::type = 0> + DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV inline DEAL_II_ALWAYS_INLINE + Tensor::type> + division_operator(const Tensor &t, + const OtherNumber & factor) + { + Tensor::type> tt; + // recurse over the base objects + for (unsigned int d = 0; d < dim; ++d) + tt[d] = t[d] / factor; + return tt; + } + } // namespace TensorImplementation +} // namespace internal + + /** * Division of a tensor of general rank with a scalar number. See the * discussion on operator*() above for more information about template @@ -1775,11 +1866,7 @@ DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV inline DEAL_II_ALWAYS_INLINE typename EnableIfScalar::type>::type> operator/(const Tensor &t, const OtherNumber &factor) { - // recurse over the base objects - Tensor::type> tt; - for (unsigned int d = 0; d < dim; ++d) - tt[d] = t[d] / factor; - return tt; + return internal::TensorImplementation::division_operator(t, factor); } diff --git a/include/deal.II/lac/block_vector_base.h b/include/deal.II/lac/block_vector_base.h index 047fd2115c..1b7bb4b23b 100644 --- a/include/deal.II/lac/block_vector_base.h +++ b/include/deal.II/lac/block_vector_base.h @@ -2100,8 +2100,10 @@ BlockVectorBase::operator/=(const value_type factor) AssertIsFinite(factor); Assert(factor != 0., ExcDivideByZero()); + const value_type inverse_factor = value_type(1.) / factor; + for (size_type i = 0; i < n_blocks(); ++i) - components[i] /= factor; + components[i] *= inverse_factor; return *this; } diff --git a/tests/data_out/data_out_base_vtu_04.with_zlib=on.output b/tests/data_out/data_out_base_vtu_04.with_zlib=on.output index 5f5f6fde26..5641ef185b 100644 --- a/tests/data_out/data_out_base_vtu_04.with_zlib=on.output +++ b/tests/data_out/data_out_base_vtu_04.with_zlib=on.output @@ -9,7 +9,7 @@ -AQAAAAASAAAAEgAARQUAAA==eNqdV01onFUUnVo0lApCkxbaZhFqJqigrbQFW5jv2VJFsCBYZqVWurFoN+LCjVJCK2jFFES6K9JIsnCykKIFBb97AxHEhcQisdCfhdRVIVBcCFER733fve+debOamc1733yXyTkv59x3T6s1G1r2ObX7/XBz9elK999dfDtcm/wq7m98/ZrWxLql3lupfufO0+HD//bEmrVzr4Z/PrgS9/sOvxiWevfivtvtpvrrO14OB5+ciN/ff/eFcHT907i/sPxM6HbXqqbq2VR/YaIT7r7yYPx+65kD4aPu+bj//fkZqVkewLay61gYm5uOe1312X4qWF0ftv0r+8Pl33bE73XV54bjvUrrSmxfnHk07Jva0vyOrPrccFyrtK7Etm1iPHz8619xr6s+N3CWq6ZueL5LvSdS/ZFtj4U929c7up95fSpcfu/NWHP453Hh2LPz357qNx4ZD2ffWY31R048HB4/fjzWfHZik3C8VPlZZT38W91cvdpp9HC/uja51878jqyz1Sh8Edu57x8I37xxK+511eeGY6/SuhLb3Z/+qO5scNzrqs8Nx0uV1pXYVnb9UI3NLVSNHhYqfW7wzFZNXb82VO+qe92rD9QPVt9umU9QG6p31b3u1QfqB91/O7PRdp+gNlTvqnvdqw/UD7r//KlP2tkn+az0/686iNxFF6qPpuaXdtZNxiaYyTELF3Iu8mHniNgEMzlm4ULORTiSc0Rsgpkcs3Ah5yIcyTkiNsFMjlm4kHMRjpT1Pxxf1IbqXXWve/WB+kH3L41x232C2lC9q+51rz5QP0QB3P6ynX2StaF6V903Z3u1o36w82xnnwzHF7EJ5toxCxdyLsKRnCNiE8y1YxYu5FyEIzlHxCaYa8csXMi5CEcC/SdtSM9k75mysvdSec+5/2dtSM9k75mysvdSec+5/2dtSM9k75mysvdSec+5/+ezkh5C3kNkZe8t8p6h/ydsonf2O0t8wH6XCV+2c+nDJnpnv7PEB+x3mfgj1SM20Tv7nSU+YL/LxB+pHrHJ/5+9h4su2Hu76CXVD8sXtSE9k7xnysreS+U95/6ftSE9k7xnykreS+U95/6ftSE9k7xnykreS+U9Qf8fii9iE72z31niA/a7TPzBf8//OV1iE72z31niA/a7TPzBD52cmy6xid7J7yzxAfldppqXmulSG6YZAi0RaIxLbZhmCLREoDEqtWGaIdASgcaoPCs7Q4KzJTjzAWzmWQYvM3icS2zmWQYvM3icS2zmWQYvM3icS2ymYQJtM2ieR+GL2jDN1KAlAo1RqQ3TTA1aItAYldowzdSgJQKN0Sh8EZt5lsDLDB7nEpt5lsDLBB7nEpt5lsDLBB6nUhs2MzDMEoZ/dtF80qcNmxkYZom43zvfWTSf9GnDZgaGWSLur588umg+6Tsru0MZ7lbX86LpRj5TCZvNbAyzHMOMF/fz9bGEzWY2hlmOYcaL+y3rzyVsNrMxzHIMM94ANpthGGYbhplnJL6oDZsZGGaJVK8+KbVhMwPDLJHq1SelNmxmIJglKONpjcR38sdNCZvNbAyzHMOMF/eHLm5O2GxmY5jlGGY8+1tXFhybzWwEsxzBjEelNiwzMmRJhow5oA3LjAxZkiFjUqkNy4wMWZIhY1J5VpahGLIVQeaiEptldoIs73xTpkNsltkJsjyVGR+xWWYnyPJUZnzEZhmWINtSmXmH5YvasMzIkCUJMiaV2rDMSJAlCTLmgDYsMxJkSYKMSaPwRWyW2QmyfF1mfMRmmZ0gy9dlxkdsltkJsnxdZvz/AR5U9Is= +AQAAAAASAAAAEgAAVwUAAA==eNqdV0FoXVUQ/VW0FhWhSQutWYSaX1SQVqrQFv67tlQR7EJa/kLUSjeKdlNcuFEkVMFGTEGkuy4ayXeRgBK0gtI3E4ggLiQWiUKkglQ3hUDQUqiKOHPfzL3n36z+z+bOyxvSc17PmTun1ZoMLfs5ef9bYXX50Urrr86dDpfGPov1z5+/qD2xb37utdS/Y8cr4b3/dsWelTMvhH/evRjrvQefCfNz12Pd7XZT/5Xtx8Ljj4zG36+/8XQ4vPZhrKcW94dud6Vqup5I/VOjnXDt+Tvi7+8+9Vg4230n1r89tVt6FjdgW9p5JGyenoi1nvpsfypYXx+2fUv7woWftsff66nPDcfrlfaV2D4+9UDYO76l+Tty6nPDcaXSvhLb1tGR8P6PN2Otpz43cBarpm9wvvNzD6f+Q1sfDLu2rXW03v3SeLjw5qux5+D3I8Jxzr7/ttR/676R8Pbry7H/0PF7wkNHj8aej45vEo7nK/9WWQ//VqvLC51GD+vVpbE99s2vyjlZDcMXsZ25fFv44uVfYq2nPjcc5yrtK7Fd++736uotjrWe+txwPF9pX4ltaec31ebp2arRw2ylzw2eyarp69eG6l11r7X6QP3QtPzQbplPUBuqd9W91uoD9YPWnx74q+0+QW2o3lX3WqsP1A/RR/un2tkn+Vvp/7/qIHIXXag+Mh7XzXrCJpjJMQsXci7yw87x19X1hE0wk2MWLuRchCM5x4W7ziZsgpkcs3Ah5yIcyTkiNsFMjlm4kHMRjpT1Pxhf1IbqXXWvtfpA/aD1vR/UbfcJakP1rrrXWn2gftD66z8+aWefZG2o3lX3zbdd6Kgfmp5n29kng/F97sblhE0w145ZuJBzEY7kHI992UvYBHPtmIULORfhSM6x1TqdsAnm2jELF3IuwpFA/0kbMjPZZ6ac7LNU3nOe/1kbMjPZZ6ac7LNU3nOe/1kbMjPZZ6ac7LNU3nOe//lbyQwhnyFyss8Wec8w/xM20Tv7nSU+YL/LhK/2tEtsonf2O0t8wH6XiT9SP2ITvbPfWeID9rtM/JH6EZv8/7PPcNEF+2wXvaT+QfmiNmRmks9MOdlnqbznPP+zNmRmks9MOclnqbznPP+zNmRmks9MOclnqbwnmP8D8UVsonf2O0t8wH6XiT/475k/J0psonf2O0t8wH6XiT/4zhPTEyU20Tv5nSU+IL/LVPPSM1FqwzRDoCUCjXGpDdMMgZYINEalNkwzBFoi0BiV38q+IcG3JfjmG7CZZxm8zOBxLrGZZxm8zOBxLrGZZxm8zOBxLrGZhgm0zaB5HoYvasM0U4OWCDRGpTZMMzVoiUBjVGrDNFODlgg0RsPwRWzmWQIvM3icS2zmWQIvE3icS2zmWQIvE3icSm3YzsCwSxj+yZ75pE8btjMw7BKx3jPT6ZlP+rRhOwPDLhHrKycO98wnfd/K7lCGu9X13DPdyM94wmY7G8Mux7DjxXqmPpKw2c7GsMsx7Hix3rL2ZMJmOxvDLsew423AZjsMw27DsPMMxRe1YTsDwy6R+tUnpTZsZ2DYJVK/+qTUhu0MBLsEZTytofiOfbspYbOdjWGXY9jxYn3g3O0Jm+1sDLscw45n/9bFWcdmOxvBLkew41GpDcuMDFmSIWNu0IZlRoYsyZAxqdSGZUaGLMmQMan8VpahGLIVQeaiEptldoIs73xTpkNsltkJsjyVGR+xWWYnyPJUZnzEZhmWINtSmXkH5YvasMzIkCUJMiaV2rDMSJAlCTLmBm1YZiTIkgQZk4bhi9gssxNk+brM+IjNMjtBlq/LjI/YLLMTZPm6zPj/A1hg9+k= diff --git a/tests/tensors/tensor_11.cc b/tests/tensors/tensor_11.cc new file mode 100644 index 0000000000..da68a93fa5 --- /dev/null +++ b/tests/tensors/tensor_11.cc @@ -0,0 +1,49 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2019 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + +// check Tensor::operator/=(int) and Tensor::operator/(int) + +#include + +#include "../tests.h" + +template +void +test() +{ + Number a[3][3] = {{1, 2, 3}, {3, 4, 5}, {6, 7, 8}}; + + const unsigned int dim = 3; + Tensor<2, dim, Number> t(a); + + deallog << "Before: " << t << std::endl; + Tensor<2, dim, Number> u = t / 2; + deallog << "operator/: " << u << std::endl; + t /= 2; + deallog << "operator/=: " << t << std::endl; +} + +int +main() +{ + initlog(); + deallog.push("Test double"); + test(); + deallog.pop(); + deallog.push("Test int"); + test(); + deallog.pop(); +} diff --git a/tests/tensors/tensor_11.output b/tests/tensors/tensor_11.output new file mode 100644 index 0000000000..c33b1a1f0a --- /dev/null +++ b/tests/tensors/tensor_11.output @@ -0,0 +1,7 @@ + +DEAL:Test double::Before: 1.00000 2.00000 3.00000 3.00000 4.00000 5.00000 6.00000 7.00000 8.00000 +DEAL:Test double::operator/: 0.500000 1.00000 1.50000 1.50000 2.00000 2.50000 3.00000 3.50000 4.00000 +DEAL:Test double::operator/=: 0.500000 1.00000 1.50000 1.50000 2.00000 2.50000 3.00000 3.50000 4.00000 +DEAL:Test int::Before: 1 2 3 3 4 5 6 7 8 +DEAL:Test int::operator/: 0 1 1 1 2 2 3 3 4 +DEAL:Test int::operator/=: 0 1 1 1 2 2 3 3 4 -- 2.39.5