Reduce floating-point division in operator/

diff --git a/include/deal.II/base/point.h b/include/deal.II/base/point.h
index f2836107ec8a66d6fd8838d0c5b33a82132d88d1..8a62844aba0c7604b0ac94b243db5a39aede7176 100644 (file)
--- 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 <int dim,
+              typename Number,
+              typename OtherNumber,
+              typename std::enable_if<
+                !std::is_integral<
+                  typename ProductType<Number, OtherNumber>::type>::value,
+                int>::type = 0>
+    inline DEAL_II_CUDA_HOST_DEV
+      Point<dim, typename ProductType<Number, OtherNumber>::type>
+      division_operator(const Point<dim, Number> &p, const OtherNumber factor)
+    {
+      using common_type = typename ProductType<Number, OtherNumber>::type;
+      Point<dim, common_type> tmp;
+      const auto              inverse_factor =
+        static_cast<common_type>(1.) / static_cast<common_type>(factor);
+      for (unsigned int i = 0; i < dim; ++i)
+        tmp[i] = p[i] * inverse_factor;
+      return tmp;
+    }
+
+    template <int dim,
+              typename Number,
+              typename OtherNumber,
+              typename std::enable_if<
+                std::is_integral<
+                  typename ProductType<Number, OtherNumber>::type>::value,
+                int>::type = 0>
+    inline DEAL_II_CUDA_HOST_DEV
+      Point<dim, typename ProductType<Number, OtherNumber>::type>
+      division_operator(const Point<dim, Number> &p, const OtherNumber factor)
+    {
+      Point<dim, typename ProductType<Number, OtherNumber>::type> tmp;
+      for (unsigned int i = 0; i < dim; ++i)
+        tmp[i] = p[i] / factor;
+      return tmp;
+    }
+  } // namespace PointImplementation
+} // namespace internal
+
+
 template <int dim, typename Number>
 template <typename OtherNumber>
 inline DEAL_II_CUDA_HOST_DEV
@@ -568,10 +612,7 @@ inline DEAL_II_CUDA_HOST_DEV
                              typename EnableIfScalar<OtherNumber>::type>::type>
   Point<dim, Number>::operator/(const OtherNumber factor) const
 {
-  Point<dim, typename ProductType<Number, OtherNumber>::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 c2fd168fcb7b78035303459adbd73f24b72675de..e7326c427eea6d7a4747b02f6001cc77a6d47301 100644 (file)
--- 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 <int rank,
+              int dim,
+              typename Number,
+              typename OtherNumber,
+              typename std::enable_if<
+                !std::is_integral<
+                  typename ProductType<Number, OtherNumber>::type>::value,
+                int>::type = 0>
+    DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV inline DEAL_II_ALWAYS_INLINE void
+                      division_operator(Tensor<rank, dim, Number> (&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 <int rank,
+              int dim,
+              typename Number,
+              typename OtherNumber,
+              typename std::enable_if<
+                std::is_integral<
+                  typename ProductType<Number, OtherNumber>::type>::value,
+                int>::type = 0>
+    DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV inline DEAL_II_ALWAYS_INLINE void
+                      division_operator(dealii::Tensor<rank, dim, Number> (&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 <int rank_, int dim, typename Number>
 template <typename OtherNumber>
 DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
   DEAL_II_CUDA_HOST_DEV Tensor<rank_, dim, Number> &
   Tensor<rank_, dim, Number>::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 <int rank,
+              int dim,
+              typename Number,
+              typename OtherNumber,
+              typename std::enable_if<
+                !std::is_integral<
+                  typename ProductType<Number, OtherNumber>::type>::value,
+                int>::type = 0>
+    DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV inline DEAL_II_ALWAYS_INLINE
+                      Tensor<rank, dim, typename ProductType<Number, OtherNumber>::type>
+                      division_operator(const Tensor<rank, dim, Number> &t,
+                                        const OtherNumber &              factor)
+    {
+      Tensor<rank, dim, typename ProductType<Number, OtherNumber>::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 <int rank,
+              int dim,
+              typename Number,
+              typename OtherNumber,
+              typename std::enable_if<
+                std::is_integral<
+                  typename ProductType<Number, OtherNumber>::type>::value,
+                int>::type = 0>
+    DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV inline DEAL_II_ALWAYS_INLINE
+                      Tensor<rank, dim, typename ProductType<Number, OtherNumber>::type>
+                      division_operator(const Tensor<rank, dim, Number> &t,
+                                        const OtherNumber &              factor)
+    {
+      Tensor<rank, dim, typename ProductType<Number, OtherNumber>::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<OtherNumber>::type>::type>
                   operator/(const Tensor<rank, dim, Number> &t, const OtherNumber &factor)
 {
-  // recurse over the base objects
-  Tensor<rank, dim, typename ProductType<Number, OtherNumber>::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 047fd2115c9b038c4474d214e37ab2cbb2471ff2..1b7bb4b23b848ccadf29e26f44ddd356690190c2 100644 (file)
--- a/include/deal.II/lac/block_vector_base.h
+++ b/include/deal.II/lac/block_vector_base.h
@@ -2100,8 +2100,10 @@ BlockVectorBase<VectorType>::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 5f5f6fde269300ec5b201b5830c94cc933d89864..5641ef185bdb1380a2e3cf17b64df99c99604954 100644 (file)
--- 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 @@
 <Piece NumberOfPoints="384" NumberOfCells="24" >
   <Points>
     <DataArray type="Float32" NumberOfComponents="3" format="binary">
-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=
     </DataArray>
   </Points>
 

diff --git a/tests/tensors/tensor_11.cc b/tests/tensors/tensor_11.cc
new file mode 100644 (file)
index 0000000..da68a93
--- /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 <deal.II/base/tensor.h>
+
+#include "../tests.h"
+
+template <typename Number>
+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<double>();
+  deallog.pop();
+  deallog.push("Test int");
+  test<int>();
+  deallog.pop();
+}

diff --git a/tests/tensors/tensor_11.output b/tests/tensors/tensor_11.output
new file mode 100644 (file)
index 0000000..c33b1a1
--- /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