From: Daniel Arndt <arndtd@ornl.gov>
Date: Fri, 5 Jul 2019 18:26:34 +0000 (-0400)
Subject: Annotate Point usable for CUDA kernels
X-Git-Tag: v9.2.0-rc1~1403^2~3
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=8e3efa78882a542c772ef96f5a3ee774d1774dc3;p=dealii.git

Annotate Point usable for CUDA kernels
---

diff --git a/doc/news/changes/minor/20190705DanielArndt b/doc/news/changes/minor/20190705DanielArndt
new file mode 100644
index 0000000000..272cdf6bd4
--- /dev/null
+++ b/doc/news/changes/minor/20190705DanielArndt
@@ -0,0 +1,4 @@
+Improved: All member funtions and free functions related to Point
+that can be used in CUDA code so far have been annotated accordingly.
+<br>
+(Daniel Arndt, 2019/07/05)
diff --git a/include/deal.II/base/point.h b/include/deal.II/base/point.h
index 11291018b5..93b14f0bb0 100644
--- a/include/deal.II/base/point.h
+++ b/include/deal.II/base/point.h
@@ -113,21 +113,28 @@ public:
   /**
    * Standard constructor. Creates an object that corresponds to the origin,
    * i.e., all coordinates are set to zero.
+   *
+   * @note This function can also be used in CUDA device code.
    */
+  DEAL_II_CUDA_HOST_DEV
   Point();
 
   /**
    * Convert a tensor to a point.
    */
-  explicit Point(const Tensor<1, dim, Number> &);
+  explicit DEAL_II_CUDA_HOST_DEV
+  Point(const Tensor<1, dim, Number> &);
 
   /**
    * Constructor for one dimensional points. This function is only implemented
    * for <tt>dim==1</tt> since the usage is considered unsafe for points with
    * <tt>dim!=1</tt> as it would leave some components of the point
    * coordinates uninitialized.
+   *
+   * @note This function can also be used in CUDA device code.
    */
-  explicit Point(const Number x);
+  explicit DEAL_II_CUDA_HOST_DEV
+  Point(const Number x);
 
   /**
    * Constructor for two dimensional points. This function is only implemented
@@ -135,7 +142,10 @@ public:
    * <tt>dim!=2</tt> as it would leave some components of the point
    * coordinates uninitialized (if dim>2) or would not use some arguments (if
    * dim<2).
+   *
+   * @note This function can also be used in CUDA device code.
    */
+  DEAL_II_CUDA_HOST_DEV
   Point(const Number x, const Number y);
 
   /**
@@ -144,7 +154,10 @@ public:
    * points with <tt>dim!=3</tt> as it would leave some components of the
    * point coordinates uninitialized (if dim>3) or would not use some
    * arguments (if dim<3).
+   *
+   * @note This function can also be used in CUDA device code.
    */
+  DEAL_II_CUDA_HOST_DEV
   Point(const Number x, const Number y, const Number z);
 
   /**
@@ -160,21 +173,27 @@ public:
    * Return a unit vector in coordinate direction <tt>i</tt>, i.e., a vector
    * that is zero in all coordinates except for a single 1 in the <tt>i</tt>th
    * coordinate.
+   *
+   * @note This function can also be used in CUDA device code.
    */
-  static Point<dim, Number>
-  unit_vector(const unsigned int i);
+  static DEAL_II_CUDA_HOST_DEV Point<dim, Number>
+                               unit_vector(const unsigned int i);
 
   /**
    * Read access to the <tt>index</tt>th coordinate.
+   *
+   * @note This function can also be used in CUDA device code.
    */
-  Number
-  operator()(const unsigned int index) const;
+  DEAL_II_CUDA_HOST_DEV Number
+                        operator()(const unsigned int index) const;
 
   /**
    * Read and write access to the <tt>index</tt>th coordinate.
+   *
+   * @note This function can also be used in CUDA device code.
    */
-  Number &
-  operator()(const unsigned int index);
+  DEAL_II_CUDA_HOST_DEV Number &
+                        operator()(const unsigned int index);
 
   /**
    * @name Addition and subtraction of points.
@@ -183,9 +202,11 @@ public:
 
   /**
    * Add an offset given as Tensor<1,dim,Number> to a point.
+   *
+   * @note This function can also be used in CUDA device code.
    */
-  Point<dim, Number>
-  operator+(const Tensor<1, dim, Number> &) const;
+  DEAL_II_CUDA_HOST_DEV Point<dim, Number>
+                        operator+(const Tensor<1, dim, Number> &) const;
 
   /**
    * Subtract two points, i.e., obtain the vector that connects the two. As
@@ -193,24 +214,30 @@ public:
    * results in a vector anchored at one of the two points (rather than at the
    * origin) and, consequently, the result is returned as a Tensor@<1,dim@>
    * rather than as a Point@<dim@>.
+   *
+   * @note This function can also be used in CUDA device code.
    */
-  Tensor<1, dim, Number>
-  operator-(const Point<dim, Number> &) const;
+  DEAL_II_CUDA_HOST_DEV Tensor<1, dim, Number>
+                        operator-(const Point<dim, Number> &) const;
 
   /**
    * Subtract a difference vector (represented by a Tensor@<1,dim@>) from the
    * current point. This results in another point and, as discussed in the
    * documentation of this class, the result is then naturally returned as a
    * Point@<dim@> object rather than as a Tensor@<1,dim@>.
+   *
+   * @note This function can also be used in CUDA device code.
    */
-  Point<dim, Number>
-  operator-(const Tensor<1, dim, Number> &) const;
+  DEAL_II_CUDA_HOST_DEV Point<dim, Number>
+                        operator-(const Tensor<1, dim, Number> &) const;
 
   /**
    * The opposite vector.
+   *
+   * @note This function can also be used in CUDA device code.
    */
-  Point<dim, Number>
-  operator-() const;
+  DEAL_II_CUDA_HOST_DEV Point<dim, Number>
+                        operator-() const;
 
   /**
    * @}
@@ -224,27 +251,35 @@ public:
   /**
    * Multiply the current point by a factor.
    *
+   * @note This function can also be used in CUDA device code.
+   *
    * @relatesalso EnableIfScalar
    */
   template <typename OtherNumber>
-  Point<dim,
-        typename ProductType<Number,
-                             typename EnableIfScalar<OtherNumber>::type>::type>
+  DEAL_II_CUDA_HOST_DEV Point<
+    dim,
+    typename ProductType<Number,
+                         typename EnableIfScalar<OtherNumber>::type>::type>
   operator*(const OtherNumber) const;
 
   /**
    * Divide the current point by a factor.
+   *
+   * @note This function can also be used in CUDA device code.
    */
   template <typename OtherNumber>
-  Point<dim,
-        typename ProductType<Number,
-                             typename EnableIfScalar<OtherNumber>::type>::type>
+  DEAL_II_CUDA_HOST_DEV Point<
+    dim,
+    typename ProductType<Number,
+                         typename EnableIfScalar<OtherNumber>::type>::type>
   operator/(const OtherNumber) const;
 
   /**
    * Return the scalar product of the vectors representing two points.
+   *
+   * @note This function can also be used in CUDA device code.
    */
-  Number operator*(const Tensor<1, dim, Number> &p) const;
+  DEAL_II_CUDA_HOST_DEV Number operator*(const Tensor<1, dim, Number> &p) const;
 
   /**
    * Return the scalar product of this point vector with itself, i.e. the
@@ -255,23 +290,29 @@ public:
    * @note This function is equivalent to
    * Tensor<rank,dim,Number>::norm_square() which returns the square of the
    * Frobenius norm.
+   *
+   * @note This function can also be used in CUDA device code.
    */
-  typename numbers::NumberTraits<Number>::real_type
+  DEAL_II_CUDA_HOST_DEV typename numbers::NumberTraits<Number>::real_type
   square() const;
 
   /**
    * Return the Euclidean distance of <tt>this</tt> point to the point
    * <tt>p</tt>, i.e. the <tt>l_2</tt> norm of the difference between the
    * vectors representing the two points.
+   *
+   * @note This function can also be used in CUDA device code.
    */
-  typename numbers::NumberTraits<Number>::real_type
+  DEAL_II_CUDA_HOST_DEV typename numbers::NumberTraits<Number>::real_type
   distance(const Point<dim, Number> &p) const;
 
   /**
    * Return the squared Euclidean distance of <tt>this</tt> point to the point
    * <tt>p</tt>.
+   *
+   * @note This function can also be used in CUDA device code.
    */
-  typename numbers::NumberTraits<Number>::real_type
+  DEAL_II_CUDA_HOST_DEV typename numbers::NumberTraits<Number>::real_type
   distance_square(const Point<dim, Number> &p) const;
 
   /**
@@ -294,20 +335,23 @@ public:
 // At least clang-3.7 requires us to have a user-defined constructor
 // and we can't use 'Point<dim,Number>::Point () = default' here.
 template <int dim, typename Number>
-inline Point<dim, Number>::Point() // NOLINT
+inline DEAL_II_CUDA_HOST_DEV
+Point<dim, Number>::Point() // NOLINT
 {}
 
 
 
 template <int dim, typename Number>
-inline Point<dim, Number>::Point(const Tensor<1, dim, Number> &t)
+inline DEAL_II_CUDA_HOST_DEV
+Point<dim, Number>::Point(const Tensor<1, dim, Number> &t)
   : Tensor<1, dim, Number>(t)
 {}
 
 
 
 template <int dim, typename Number>
-inline Point<dim, Number>::Point(const Number x)
+inline DEAL_II_CUDA_HOST_DEV
+Point<dim, Number>::Point(const Number x)
 {
   Assert(dim == 1,
          ExcMessage(
@@ -332,7 +376,8 @@ inline Point<dim, Number>::Point(const Number x)
 
 
 template <int dim, typename Number>
-inline Point<dim, Number>::Point(const Number x, const Number y)
+inline DEAL_II_CUDA_HOST_DEV
+Point<dim, Number>::Point(const Number x, const Number y)
 {
   Assert(dim == 2,
          ExcMessage(
@@ -351,7 +396,8 @@ inline Point<dim, Number>::Point(const Number x, const Number y)
 
 
 template <int dim, typename Number>
-inline Point<dim, Number>::Point(const Number x, const Number y, const Number z)
+inline DEAL_II_CUDA_HOST_DEV
+Point<dim, Number>::Point(const Number x, const Number y, const Number z)
 {
   Assert(dim == 3,
          ExcMessage(
@@ -395,8 +441,8 @@ inline Point<dim, Number>::Point(
 
 
 template <int dim, typename Number>
-inline Point<dim, Number>
-Point<dim, Number>::unit_vector(unsigned int i)
+inline DEAL_II_CUDA_HOST_DEV Point<dim, Number>
+                             Point<dim, Number>::unit_vector(unsigned int i)
 {
   Point<dim, Number> p;
   p[i] = 1.;
@@ -405,7 +451,7 @@ Point<dim, Number>::unit_vector(unsigned int i)
 
 
 template <int dim, typename Number>
-inline Number
+inline DEAL_II_CUDA_HOST_DEV Number
 Point<dim, Number>::operator()(const unsigned int index) const
 {
   AssertIndexRange(index, dim);
@@ -415,7 +461,7 @@ Point<dim, Number>::operator()(const unsigned int index) const
 
 
 template <int dim, typename Number>
-inline Number &
+inline DEAL_II_CUDA_HOST_DEV Number &
 Point<dim, Number>::operator()(const unsigned int index)
 {
   AssertIndexRange(index, dim);
@@ -425,7 +471,7 @@ Point<dim, Number>::operator()(const unsigned int index)
 
 
 template <int dim, typename Number>
-inline Point<dim, Number>
+inline DEAL_II_CUDA_HOST_DEV Point<dim, Number>
 Point<dim, Number>::operator+(const Tensor<1, dim, Number> &p) const
 {
   Point<dim, Number> tmp = *this;
@@ -436,7 +482,7 @@ Point<dim, Number>::operator+(const Tensor<1, dim, Number> &p) const
 
 
 template <int dim, typename Number>
-inline Tensor<1, dim, Number>
+inline DEAL_II_CUDA_HOST_DEV Tensor<1, dim, Number>
 Point<dim, Number>::operator-(const Point<dim, Number> &p) const
 {
   return (Tensor<1, dim, Number>(*this) -= p);
@@ -445,7 +491,7 @@ Point<dim, Number>::operator-(const Point<dim, Number> &p) const
 
 
 template <int dim, typename Number>
-inline Point<dim, Number>
+inline DEAL_II_CUDA_HOST_DEV Point<dim, Number>
 Point<dim, Number>::operator-(const Tensor<1, dim, Number> &p) const
 {
   Point<dim, Number> tmp = *this;
@@ -456,7 +502,7 @@ Point<dim, Number>::operator-(const Tensor<1, dim, Number> &p) const
 
 
 template <int dim, typename Number>
-inline Point<dim, Number>
+inline DEAL_II_CUDA_HOST_DEV Point<dim, Number>
 Point<dim, Number>::operator-() const
 {
   Point<dim, Number> result;
@@ -469,11 +515,11 @@ Point<dim, Number>::operator-() const
 
 template <int dim, typename Number>
 template <typename OtherNumber>
-inline Point<
-  dim,
-  typename ProductType<Number,
-                       typename EnableIfScalar<OtherNumber>::type>::type>
-  Point<dim, Number>::operator*(const OtherNumber factor) const
+inline DEAL_II_CUDA_HOST_DEV
+    Point<dim,
+        typename ProductType<Number,
+                             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)
@@ -485,11 +531,11 @@ inline Point<
 
 template <int dim, typename Number>
 template <typename OtherNumber>
-inline Point<
-  dim,
-  typename ProductType<Number,
-                       typename EnableIfScalar<OtherNumber>::type>::type>
-Point<dim, Number>::operator/(const OtherNumber factor) const
+inline DEAL_II_CUDA_HOST_DEV
+  Point<dim,
+        typename ProductType<Number,
+                             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)
@@ -500,8 +546,8 @@ Point<dim, Number>::operator/(const OtherNumber factor) const
 
 
 template <int dim, typename Number>
-inline Number Point<dim, Number>::
-              operator*(const Tensor<1, dim, Number> &p) const
+inline DEAL_II_CUDA_HOST_DEV Number Point<dim, Number>::
+                                    operator*(const Tensor<1, dim, Number> &p) const
 {
   Number res = Number();
   for (unsigned int i = 0; i < dim; ++i)
@@ -511,7 +557,7 @@ inline Number Point<dim, Number>::
 
 
 template <int dim, typename Number>
-inline typename numbers::NumberTraits<Number>::real_type
+inline DEAL_II_CUDA_HOST_DEV typename numbers::NumberTraits<Number>::real_type
 Point<dim, Number>::square() const
 {
   return this->norm_square();
@@ -520,7 +566,7 @@ Point<dim, Number>::square() const
 
 
 template <int dim, typename Number>
-inline typename numbers::NumberTraits<Number>::real_type
+inline DEAL_II_CUDA_HOST_DEV typename numbers::NumberTraits<Number>::real_type
 Point<dim, Number>::distance(const Point<dim, Number> &p) const
 {
   return std::sqrt(distance_square(p));
@@ -529,7 +575,7 @@ Point<dim, Number>::distance(const Point<dim, Number> &p) const
 
 
 template <int dim, typename Number>
-inline typename numbers::NumberTraits<Number>::real_type
+inline DEAL_II_CUDA_HOST_DEV typename numbers::NumberTraits<Number>::real_type
 Point<dim, Number>::distance_square(const Point<dim, Number> &p) const
 {
   Number sum = internal::NumberType<Number>::value(0.0);
@@ -563,15 +609,17 @@ Point<dim, Number>::serialize(Archive &ar, const unsigned int)
 /**
  * Global operator scaling a point vector by a scalar.
  *
+ * @note This function can also be used in CUDA device code.
+ *
  * @relatesalso Point
  * @relatesalso EnableIfScalar
  */
 template <int dim, typename Number, typename OtherNumber>
-inline Point<
-  dim,
-  typename ProductType<Number,
-                       typename EnableIfScalar<OtherNumber>::type>::type>
-operator*(const OtherNumber factor, const Point<dim, Number> &p)
+inline DEAL_II_CUDA_HOST_DEV
+  Point<dim,
+        typename ProductType<Number,
+                             typename EnableIfScalar<OtherNumber>::type>::type>
+  operator*(const OtherNumber factor, const Point<dim, Number> &p)
 {
   return p * factor;
 }
diff --git a/tests/cuda/cuda_point.cu b/tests/cuda/cuda_point.cu
new file mode 100644
index 0000000000..43a285693f
--- /dev/null
+++ b/tests/cuda/cuda_point.cu
@@ -0,0 +1,80 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2016 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.
+//
+// ---------------------------------------------------------------------
+
+// Test that Point operations on a CUDA device can be used.
+
+#include <deal.II/base/point.h>
+
+#include "../tests.h"
+
+template <int dim, typename Number>
+__global__ void
+miscellaneous_kernel()
+{
+  Point<dim, Number> p_1;
+  Point<dim, Number> p_2(Tensor<1, dim, Number>{});
+  if (dim == 1)
+    Point<dim, Number> p(1.);
+  if (dim == 2)
+    Point<dim, Number> p(1., 2.);
+  if (dim == 3)
+    Point<dim, Number> p(1., 2., 3.);
+
+  auto p_3 = Point<dim, Number>::unit_vector(0);
+
+  auto entry_1 = p_1(0);
+  p_1(0)       = Number{1.};
+
+  auto p_4 = p_1 + Tensor<1, dim, Number>{};
+  auto p_5 = p_1 - Tensor<1, dim, Number>{};
+  auto t_1 = p_1 - p_2;
+  auto p_6 = -p_3;
+  auto p_7 = p_4 / 2.;
+  auto p_8 = p_2 * 5.;
+
+  auto s_1 = p_1 * t_1;
+  auto s_2 = p_2.square();
+  auto s_3 = p_3.distance(p_5);
+  auto s_4 = p_4.distance_square(p_1);
+}
+
+template <int dim, typename Number>
+void
+test_gpu()
+{
+  // Miscellaneous
+  miscellaneous_kernel<dim, Number><<<1, 1>>>();
+  // Check that the kernel was launched correctly
+  AssertCuda(cudaGetLastError());
+  // Check that there was no problem during the execution of the kernel
+  AssertCuda(cudaDeviceSynchronize());
+
+  deallog << "OK" << std::endl;
+}
+
+int
+main()
+{
+  initlog();
+
+  init_cuda();
+
+  test_gpu<1, double>();
+  test_gpu<2, double>();
+  test_gpu<3, float>();
+  test_gpu<1, float>();
+  test_gpu<2, float>();
+  test_gpu<3, float>();
+}
diff --git a/tests/cuda/cuda_point.output b/tests/cuda/cuda_point.output
new file mode 100644
index 0000000000..0aa61ff573
--- /dev/null
+++ b/tests/cuda/cuda_point.output
@@ -0,0 +1,7 @@
+
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK
+DEAL::OK