Annotate Tensor functions that can be used in CUDA device code

diff --git a/include/deal.II/base/tensor.h b/include/deal.II/base/tensor.h
index d9f8f4c676f012007c63e9570241376f94cc4999..961ceb7a18c228a717995e5bdedb2192b7fe895a 100644 (file)
--- a/include/deal.II/base/tensor.h
+++ b/include/deal.II/base/tensor.h
@@ -150,15 +150,21 @@ public:
    * Constructor from tensors with different underlying scalar type. This
    * obviously requires that the @p OtherNumber type is convertible to @p
    * Number.
+   *
+   * @note This function can also be used in CUDA device code.
    */
   template <typename OtherNumber>
-  constexpr Tensor(const Tensor<0, dim, OtherNumber> &initializer);
+  constexpr DEAL_II_CUDA_HOST_DEV
+  Tensor(const Tensor<0, dim, OtherNumber> &initializer);
 
   /**
    * Constructor, where the data is copied from a C-style array.
+   *
+   * @note This function can also be used in CUDA device code.
    */
   template <typename OtherNumber>
-  constexpr Tensor(const OtherNumber &initializer);
+  constexpr DEAL_II_CUDA_HOST_DEV
+  Tensor(const OtherNumber &initializer);
 
   /**
    * Return a pointer to the first element of the underlying storage.
@@ -210,10 +216,12 @@ public:
    * Assignment from tensors with different underlying scalar type. This
    * obviously requires that the @p OtherNumber type is convertible to @p
    * Number.
+   *
+   * @note This function can also be used in CUDA device code.
    */
   template <typename OtherNumber>
-  DEAL_II_CONSTEXPR Tensor &
-                    operator=(const Tensor<0, dim, OtherNumber> &rhs);
+  DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV Tensor &
+                                          operator=(const Tensor<0, dim, OtherNumber> &rhs);
 
 #ifdef __INTEL_COMPILER
   /**
@@ -221,18 +229,22 @@ public:
    * This is needed for ICC15 because it can't generate a suitable
    * copy constructor for Sacado::Rad::ADvar types automatically.
    * See https://github.com/dealii/dealii/pull/5865.
+   *
+   * @note This function can also be used in CUDA device code.
    */
-  DEAL_II_CONSTEXPR Tensor &
-                    operator=(const Tensor<0, dim, Number> &rhs);
+  DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV Tensor &
+                                          operator=(const Tensor<0, dim, Number> &rhs);
 #endif
 
   /**
    * This operator assigns a scalar to a tensor. This obviously requires
    * that the @p OtherNumber type is convertible to @p Number.
+   *
+   * @note This function can also be used in CUDA device code.
    */
   template <typename OtherNumber>
-  DEAL_II_CONSTEXPR Tensor &
-                    operator=(const OtherNumber &d);
+  DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV Tensor &
+                                          operator=(const OtherNumber &d);
 
   /**
    * Test for equality of two tensors.
@@ -249,18 +261,22 @@ public:
   operator!=(const Tensor<0, dim, OtherNumber> &rhs) const;
 
   /**
-   * Add another scalar
+   * Add another scalar.
+   *
+   * @note This function can also be used in CUDA device code.
    */
   template <typename OtherNumber>
-  DEAL_II_CONSTEXPR Tensor &
-                    operator+=(const Tensor<0, dim, OtherNumber> &rhs);
+  DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV Tensor &
+                                          operator+=(const Tensor<0, dim, OtherNumber> &rhs);
 
   /**
    * Subtract another scalar.
+   *
+   * @note This function can also be used in CUDA device code.
    */
   template <typename OtherNumber>
-  DEAL_II_CONSTEXPR Tensor &
-                    operator-=(const Tensor<0, dim, OtherNumber> &rhs);
+  DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV Tensor &
+                                          operator-=(const Tensor<0, dim, OtherNumber> &rhs);
 
   /**
    * Multiply the scalar with a <tt>factor</tt>.
@@ -273,16 +289,20 @@ public:
 
   /**
    * Divide the scalar by <tt>factor</tt>.
+   *
+   * @note This function can also be used in CUDA device code.
    */
   template <typename OtherNumber>
-  DEAL_II_CONSTEXPR Tensor &
-                    operator/=(const OtherNumber &factor);
+  DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV Tensor &
+                                          operator/=(const OtherNumber &factor);
 
   /**
    * Tensor with inverted entries.
+   *
+   * @note This function can also be used in CUDA device code.
    */
-  constexpr Tensor
-  operator-() const;
+  constexpr DEAL_II_CUDA_HOST_DEV Tensor
+                                  operator-() const;
 
   /**
    * Reset all values to zero.
@@ -303,9 +323,11 @@ public:
    * Return the Frobenius-norm of a tensor, i.e. the square root of the sum of
    * the absolute squares of all entries. For the present case of rank-1
    * tensors, this equals the usual <tt>l<sub>2</sub></tt> norm of the vector.
+   *
+   * @note This function can also be used in CUDA device code.
    */
-  DEAL_II_CONSTEXPR real_type
-                    norm() const;
+  DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV real_type
+                                          norm() const;
 
   /**
    * Return the square of the Frobenius-norm of a tensor, i.e. the sum of the
@@ -444,16 +466,22 @@ public:
 
   /**
    * Constructor, where the data is copied from a C-style array.
+   *
+   * @note This function can also be used in CUDA device code.
    */
-  constexpr explicit Tensor(const array_type &initializer);
+  constexpr DEAL_II_CUDA_HOST_DEV explicit Tensor(
+    const array_type &initializer);
 
   /**
    * Constructor from tensors with different underlying scalar type. This
    * obviously requires that the @p OtherNumber type is convertible to @p
    * Number.
+   *
+   * @note This function can also be used in CUDA device code.
    */
   template <typename OtherNumber>
-  constexpr Tensor(const Tensor<rank_, dim, OtherNumber> &initializer);
+  constexpr DEAL_II_CUDA_HOST_DEV
+  Tensor(const Tensor<rank_, dim, OtherNumber> &initializer);
 
   /**
    * Constructor that converts from a "tensor of tensors".
@@ -524,10 +552,12 @@ public:
    * Assignment operator from tensors with different underlying scalar type.
    * This obviously requires that the @p OtherNumber type is convertible to @p
    * Number.
+   *
+   * @note This function can also be used in CUDA device code.
    */
   template <typename OtherNumber>
-  DEAL_II_CONSTEXPR Tensor &
-                    operator=(const Tensor<rank_, dim, OtherNumber> &rhs);
+  DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV Tensor &
+                                          operator=(const Tensor<rank_, dim, OtherNumber> &rhs);
 
   /**
    * This operator assigns a scalar to a tensor. To avoid confusion with what
@@ -554,17 +584,21 @@ public:
 
   /**
    * Add another tensor.
+   *
+   * @note This function can also be used in CUDA device code.
    */
   template <typename OtherNumber>
-  DEAL_II_CONSTEXPR Tensor &
-                    operator+=(const Tensor<rank_, dim, OtherNumber> &);
+  DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV Tensor &
+                                          operator+=(const Tensor<rank_, dim, OtherNumber> &);
 
   /**
    * Subtract another tensor.
+   *
+   * @note This function can also be used in CUDA device code.
    */
   template <typename OtherNumber>
-  DEAL_II_CONSTEXPR Tensor &
-                    operator-=(const Tensor<rank_, dim, OtherNumber> &);
+  DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV Tensor &
+                                          operator-=(const Tensor<rank_, dim, OtherNumber> &);
 
   /**
    * Scale the tensor by <tt>factor</tt>, i.e. multiply all components by
@@ -578,16 +612,20 @@ public:
 
   /**
    * Scale the vector by <tt>1/factor</tt>.
+   *
+   * @note This function can also be used in CUDA device code.
    */
   template <typename OtherNumber>
-  DEAL_II_CONSTEXPR Tensor &
-                    operator/=(const OtherNumber &factor);
+  DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV Tensor &
+                                          operator/=(const OtherNumber &factor);
 
   /**
    * Unary minus operator. Negate all entries of a tensor.
+   *
+   * @note This function can also be used in CUDA device code.
    */
-  DEAL_II_CONSTEXPR Tensor
-                    operator-() const;
+  DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV Tensor
+                                          operator-() const;
 
   /**
    * Reset all values to zero.
@@ -608,10 +646,13 @@ public:
    * Return the Frobenius-norm of a tensor, i.e. the square root of the sum of
    * the absolute squares of all entries. For the present case of rank-1
    * tensors, this equals the usual <tt>l<sub>2</sub></tt> norm of the vector.
+   *
+   * @note This function can also be used in CUDA device code.
    */
 
-  DEAL_II_CONSTEXPR typename numbers::NumberTraits<Number>::real_type
-  norm() const;
+  DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV
+    typename numbers::NumberTraits<Number>::real_type
+    norm() const;
 
   /**
    * Return the square of the Frobenius-norm of a tensor, i.e. the sum of the
@@ -687,11 +728,13 @@ private:
 
   /**
    * This constructor is for internal use. It provides a way
-   *  to create constexpr constructors for Tensor<rank, dim, Number>
+   * to create constexpr constructors for Tensor<rank, dim, Number>
+   *
+   * @note This function can also be used in CUDA device code.
    */
   template <typename ArrayLike, std::size_t... Indices>
-  constexpr Tensor(const ArrayLike &initializer,
-                   std_cxx14::index_sequence<Indices...>);
+  constexpr DEAL_II_CUDA_HOST_DEV
+  Tensor(const ArrayLike &initializer, std_cxx14::index_sequence<Indices...>);
 
   /**
    * Allow an arbitrary Tensor to access the underlying values.
@@ -783,7 +826,8 @@ Tensor<0, dim, Number>::Tensor()
 
 template <int dim, typename Number>
 template <typename OtherNumber>
-constexpr Tensor<0, dim, Number>::Tensor(const OtherNumber &initializer)
+constexpr DEAL_II_CUDA_HOST_DEV
+Tensor<0, dim, Number>::Tensor(const OtherNumber &initializer)
   : value(internal::NumberType<Number>::value(initializer))
 {}
 
@@ -791,7 +835,8 @@ constexpr Tensor<0, dim, Number>::Tensor(const OtherNumber &initializer)
 
 template <int dim, typename Number>
 template <typename OtherNumber>
-constexpr Tensor<0, dim, Number>::Tensor(const Tensor<0, dim, OtherNumber> &p)
+constexpr DEAL_II_CUDA_HOST_DEV
+Tensor<0, dim, Number>::Tensor(const Tensor<0, dim, OtherNumber> &p)
   : Tensor{p.value}
 {}
 
@@ -861,7 +906,7 @@ DEAL_II_CONSTEXPR inline DEAL_II_CUDA_HOST_DEV Tensor<0, dim, Number>::
 
 template <int dim, typename Number>
 template <typename OtherNumber>
-DEAL_II_CONSTEXPR inline Tensor<0, dim, Number> &
+DEAL_II_CONSTEXPR inline DEAL_II_CUDA_HOST_DEV Tensor<0, dim, Number> &
 Tensor<0, dim, Number>::operator=(const Tensor<0, dim, OtherNumber> &p)
 {
   value = internal::NumberType<Number>::value(p);
@@ -871,7 +916,7 @@ Tensor<0, dim, Number>::operator=(const Tensor<0, dim, OtherNumber> &p)
 
 #ifdef __INTEL_COMPILER
 template <int dim, typename Number>
-DEAL_II_CONSTEXPR inline Tensor<0, dim, Number> &
+DEAL_II_CONSTEXPR inline DEAL_II_CUDA_HOS_DEV Tensor<0, dim, Number> &
 Tensor<0, dim, Number>::operator=(const Tensor<0, dim, Number> &p)
 {
   value = p.value;
@@ -882,7 +927,7 @@ Tensor<0, dim, Number>::operator=(const Tensor<0, dim, Number> &p)
 
 template <int dim, typename Number>
 template <typename OtherNumber>
-DEAL_II_CONSTEXPR inline Tensor<0, dim, Number> &
+DEAL_II_CONSTEXPR inline DEAL_II_CUDA_HOST_DEV Tensor<0, dim, Number> &
 Tensor<0, dim, Number>::operator=(const OtherNumber &d)
 {
   value = internal::NumberType<Number>::value(d);
@@ -895,7 +940,8 @@ template <typename OtherNumber>
 DEAL_II_CONSTEXPR inline bool
 Tensor<0, dim, Number>::operator==(const Tensor<0, dim, OtherNumber> &p) const
 {
-#ifdef DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING
+#if defined(DEAL_II_ADOLC_WITH_ADVANCED_BRANCHING) && \
+  !defined(DEAL_II_COMPILER_CUDA_AWARE)
   Assert(!(std::is_same<Number, adouble>::value ||
            std::is_same<OtherNumber, adouble>::value),
          ExcMessage(
@@ -918,7 +964,7 @@ Tensor<0, dim, Number>::operator!=(const Tensor<0, dim, OtherNumber> &p) const
 
 template <int dim, typename Number>
 template <typename OtherNumber>
-DEAL_II_CONSTEXPR inline Tensor<0, dim, Number> &
+DEAL_II_CONSTEXPR inline DEAL_II_CUDA_HOST_DEV Tensor<0, dim, Number> &
 Tensor<0, dim, Number>::operator+=(const Tensor<0, dim, OtherNumber> &p)
 {
   value += p.value;
@@ -928,7 +974,7 @@ Tensor<0, dim, Number>::operator+=(const Tensor<0, dim, OtherNumber> &p)
 
 template <int dim, typename Number>
 template <typename OtherNumber>
-DEAL_II_CONSTEXPR inline Tensor<0, dim, Number> &
+DEAL_II_CONSTEXPR inline DEAL_II_CUDA_HOST_DEV Tensor<0, dim, Number> &
 Tensor<0, dim, Number>::operator-=(const Tensor<0, dim, OtherNumber> &p)
 {
   value -= p.value;
@@ -974,7 +1020,7 @@ Tensor<0, dim, Number>::operator*=(const OtherNumber &s)
 
 template <int dim, typename Number>
 template <typename OtherNumber>
-DEAL_II_CONSTEXPR inline Tensor<0, dim, Number> &
+DEAL_II_CONSTEXPR inline DEAL_II_CUDA_HOST_DEV Tensor<0, dim, Number> &
 Tensor<0, dim, Number>::operator/=(const OtherNumber &s)
 {
   value /= s;
@@ -983,7 +1029,7 @@ Tensor<0, dim, Number>::operator/=(const OtherNumber &s)
 
 
 template <int dim, typename Number>
-constexpr Tensor<0, dim, Number>
+constexpr DEAL_II_CUDA_HOST_DEV Tensor<0, dim, Number>
 Tensor<0, dim, Number>::operator-() const
 {
   return -value;
@@ -1049,9 +1095,9 @@ Tensor<0, dim, Number>::serialize(Archive &ar, const unsigned int)
 
 template <int rank_, int dim, typename Number>
 template <typename ArrayLike, std::size_t... indices>
-DEAL_II_ALWAYS_INLINE constexpr Tensor<rank_, dim, Number>::Tensor(
-  const ArrayLike &initializer,
-  std_cxx14::index_sequence<indices...>)
+DEAL_II_ALWAYS_INLINE constexpr DEAL_II_CUDA_HOST_DEV
+Tensor<rank_, dim, Number>::Tensor(const ArrayLike &initializer,
+                                   std_cxx14::index_sequence<indices...>)
   : values{Tensor<rank_ - 1, dim, Number>(initializer[indices])...}
 {
   static_assert(sizeof...(indices) == dim,
@@ -1060,15 +1106,16 @@ DEAL_II_ALWAYS_INLINE constexpr Tensor<rank_, dim, Number>::Tensor(
 
 
 template <int rank_, int dim, typename Number>
-DEAL_II_ALWAYS_INLINE constexpr Tensor<rank_, dim, Number>::Tensor(
-  const array_type &initializer)
+DEAL_II_ALWAYS_INLINE constexpr DEAL_II_CUDA_HOST_DEV
+Tensor<rank_, dim, Number>::Tensor(const array_type &initializer)
   : Tensor(initializer, std_cxx14::make_index_sequence<dim>{})
 {}
 
 
 template <int rank_, int dim, typename Number>
 template <typename OtherNumber>
-DEAL_II_ALWAYS_INLINE constexpr Tensor<rank_, dim, Number>::Tensor(
+DEAL_II_ALWAYS_INLINE constexpr DEAL_II_CUDA_HOST_DEV
+Tensor<rank_, dim, Number>::Tensor(
   const Tensor<rank_, dim, OtherNumber> &initializer)
   : Tensor(initializer, std_cxx14::make_index_sequence<dim>{})
 {}
@@ -1154,8 +1201,8 @@ DEAL_II_ALWAYS_INLINE constexpr DEAL_II_CUDA_HOST_DEV const typename Tensor<
   Number>::value_type &Tensor<rank_, dim, Number>::
                        operator[](const unsigned int i) const
 {
-  return dealii::internal::TensorSubscriptor::subscript(
-    values, i, std::integral_constant<int, dim>());
+  return values[i]; /*dealii::internal::TensorSubscriptor::subscript(
+     values, i, std::integral_constant<int, dim>());*/
 }
 
 
@@ -1285,7 +1332,7 @@ operator!=(const Tensor<rank_, dim, OtherNumber> &p) const
 
 template <int rank_, int dim, typename Number>
 template <typename OtherNumber>
-DEAL_II_CONSTEXPR inline Tensor<rank_, dim, Number> &
+DEAL_II_CONSTEXPR inline DEAL_II_CUDA_HOST_DEV Tensor<rank_, dim, Number> &
 Tensor<rank_, dim, Number>::operator+=(const Tensor<rank_, dim, OtherNumber> &p)
 {
   for (unsigned int i = 0; i < dim; ++i)
@@ -1296,7 +1343,7 @@ Tensor<rank_, dim, Number>::operator+=(const Tensor<rank_, dim, OtherNumber> &p)
 
 template <int rank_, int dim, typename Number>
 template <typename OtherNumber>
-DEAL_II_CONSTEXPR inline Tensor<rank_, dim, Number> &
+DEAL_II_CONSTEXPR inline DEAL_II_CUDA_HOST_DEV Tensor<rank_, dim, Number> &
 Tensor<rank_, dim, Number>::operator-=(const Tensor<rank_, dim, OtherNumber> &p)
 {
   for (unsigned int i = 0; i < dim; ++i)
@@ -1318,7 +1365,7 @@ Tensor<rank_, dim, Number>::operator*=(const OtherNumber &s)
 
 template <int rank_, int dim, typename Number>
 template <typename OtherNumber>
-DEAL_II_CONSTEXPR inline Tensor<rank_, dim, Number> &
+DEAL_II_CONSTEXPR 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)
@@ -1328,7 +1375,7 @@ Tensor<rank_, dim, Number>::operator/=(const OtherNumber &s)
 
 
 template <int rank_, int dim, typename Number>
-DEAL_II_CONSTEXPR inline Tensor<rank_, dim, Number>
+DEAL_II_CONSTEXPR inline DEAL_II_CUDA_HOST_DEV Tensor<rank_, dim, Number>
 Tensor<rank_, dim, Number>::operator-() const
 {
   Tensor<rank_, dim, Number> tmp;
@@ -1543,11 +1590,14 @@ operator<<(std::ostream &out, const Tensor<0, dim, Number> &p)
  * This function unwraps the underlying @p Number stored in the Tensor and
  * multiplies @p object with it.
  *
+ * @note This function can also be used in CUDA device code.
+ *
  * @relatesalso Tensor<0,dim,Number>
  */
 template <int dim, typename Number, typename Other>
-constexpr DEAL_II_ALWAYS_INLINE typename ProductType<Other, Number>::type
-operator*(const Other &object, const Tensor<0, dim, Number> &t)
+constexpr DEAL_II_CUDA_HOST_DEV DEAL_II_ALWAYS_INLINE
+  typename ProductType<Other, Number>::type
+  operator*(const Other &object, const Tensor<0, dim, Number> &t)
 {
   return object * static_cast<const Number &>(t);
 }
@@ -1560,11 +1610,14 @@ operator*(const Other &object, const Tensor<0, dim, Number> &t)
  * This function unwraps the underlying @p Number stored in the Tensor and
  * multiplies @p object with it.
  *
+ * @note This function can also be used in CUDA device code.
+ *
  * @relatesalso Tensor<0,dim,Number>
  */
 template <int dim, typename Number, typename Other>
-constexpr DEAL_II_ALWAYS_INLINE typename ProductType<Number, Other>::type
-operator*(const Tensor<0, dim, Number> &t, const Other &object)
+constexpr DEAL_II_CUDA_HOST_DEV DEAL_II_ALWAYS_INLINE
+  typename ProductType<Number, Other>::type
+  operator*(const Tensor<0, dim, Number> &t, const Other &object)
 {
   return static_cast<const Number &>(t) * object;
 }
@@ -1577,12 +1630,15 @@ operator*(const Tensor<0, dim, Number> &t, const Other &object)
  * OtherNumber that are stored within the Tensor and multiplies them. It
  * returns an unwrapped number of product type.
  *
+ * @note This function can also be used in CUDA device code.
+ *
  * @relatesalso Tensor<0,dim,Number>
  */
 template <int dim, typename Number, typename OtherNumber>
-constexpr DEAL_II_ALWAYS_INLINE typename ProductType<Number, OtherNumber>::type
-operator*(const Tensor<0, dim, Number> &     src1,
-          const Tensor<0, dim, OtherNumber> &src2)
+DEAL_II_CUDA_HOST_DEV constexpr DEAL_II_ALWAYS_INLINE
+  typename ProductType<Number, OtherNumber>::type
+  operator*(const Tensor<0, dim, Number> &     src1,
+            const Tensor<0, dim, OtherNumber> &src2)
 {
   return static_cast<const Number &>(src1) *
          static_cast<const OtherNumber &>(src2);
@@ -1592,10 +1648,12 @@ operator*(const Tensor<0, dim, Number> &     src1,
 /**
  * Division of a tensor of rank 0 by a scalar number.
  *
+ * @note This function can also be used in CUDA device code.
+ *
  * @relatesalso Tensor<0,dim,Number>
  */
 template <int dim, typename Number, typename OtherNumber>
-constexpr DEAL_II_ALWAYS_INLINE
+DEAL_II_CUDA_HOST_DEV constexpr DEAL_II_ALWAYS_INLINE
   Tensor<0,
          dim,
          typename ProductType<Number,
@@ -1609,12 +1667,14 @@ constexpr DEAL_II_ALWAYS_INLINE
 /**
  * Add two tensors of rank 0.
  *
+ * @note This function can also be used in CUDA device code.
+ *
  * @relatesalso Tensor<0,dim,Number>
  */
 template <int dim, typename Number, typename OtherNumber>
-constexpr DEAL_II_ALWAYS_INLINE
-  Tensor<0, dim, typename ProductType<Number, OtherNumber>::type>
-  operator+(const Tensor<0, dim, Number> &     p,
+constexpr DEAL_II_ALWAYS_INLINE DEAL_II_CUDA_HOST_DEV
+                                Tensor<0, dim, typename ProductType<Number, OtherNumber>::type>
+                                operator+(const Tensor<0, dim, Number> &     p,
             const Tensor<0, dim, OtherNumber> &q)
 {
   return static_cast<const Number &>(p) + static_cast<const OtherNumber &>(q);
@@ -1624,12 +1684,14 @@ constexpr DEAL_II_ALWAYS_INLINE
 /**
  * Subtract two tensors of rank 0.
  *
+ * @note This function can also be used in CUDA device code.
+ *
  * @relatesalso Tensor<0,dim,Number>
  */
 template <int dim, typename Number, typename OtherNumber>
-constexpr DEAL_II_ALWAYS_INLINE
-  Tensor<0, dim, typename ProductType<Number, OtherNumber>::type>
-  operator-(const Tensor<0, dim, Number> &     p,
+constexpr DEAL_II_ALWAYS_INLINE DEAL_II_CUDA_HOST_DEV
+                                Tensor<0, dim, typename ProductType<Number, OtherNumber>::type>
+                                operator-(const Tensor<0, dim, Number> &     p,
             const Tensor<0, dim, OtherNumber> &q)
 {
   return static_cast<const Number &>(p) - static_cast<const OtherNumber &>(q);
@@ -1644,15 +1706,17 @@ constexpr DEAL_II_ALWAYS_INLINE
  * number, a complex floating point number, etc.) is allowed, see the
  * documentation of EnableIfScalar for details.
  *
+ * @note This function can also be used in CUDA device code.
+ *
  * @relatesalso Tensor
  */
 template <int rank, int dim, typename Number, typename OtherNumber>
-DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
-  Tensor<rank,
+DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV inline DEAL_II_ALWAYS_INLINE
+                  Tensor<rank,
          dim,
          typename ProductType<Number,
                               typename EnableIfScalar<OtherNumber>::type>::type>
-  operator*(const Tensor<rank, dim, Number> &t, const OtherNumber &factor)
+                  operator*(const Tensor<rank, dim, Number> &t, const OtherNumber &factor)
 {
   // recurse over the base objects
   Tensor<rank, dim, typename ProductType<Number, OtherNumber>::type> tt;
@@ -1670,10 +1734,12 @@ DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
  * number, a complex floating point number, etc.) is allowed, see the
  * documentation of EnableIfScalar for details.
  *
+ * @note This function can also be used in CUDA device code.
+ *
  * @relatesalso Tensor
  */
 template <int rank, int dim, typename Number, typename OtherNumber>
-constexpr DEAL_II_ALWAYS_INLINE
+DEAL_II_CUDA_HOST_DEV constexpr DEAL_II_ALWAYS_INLINE
   Tensor<rank,
          dim,
          typename ProductType<typename EnableIfScalar<Number>::type,
@@ -1690,15 +1756,17 @@ constexpr DEAL_II_ALWAYS_INLINE
  * discussion on operator*() above for more information about template
  * arguments and the return type.
  *
+ * @note This function can also be used in CUDA device code.
+ *
  * @relatesalso Tensor
  */
 template <int rank, int dim, typename Number, typename OtherNumber>
-DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
-  Tensor<rank,
+DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV inline DEAL_II_ALWAYS_INLINE
+                  Tensor<rank,
          dim,
          typename ProductType<Number,
                               typename EnableIfScalar<OtherNumber>::type>::type>
-  operator/(const Tensor<rank, dim, Number> &t, const OtherNumber &factor)
+                  operator/(const Tensor<rank, dim, Number> &t, const OtherNumber &factor)
 {
   // recurse over the base objects
   Tensor<rank, dim, typename ProductType<Number, OtherNumber>::type> tt;
@@ -1713,12 +1781,14 @@ DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
  *
  * @tparam rank The rank of both tensors.
  *
+ * @note This function can also be used in CUDA device code.
+ *
  * @relatesalso Tensor
  */
 template <int rank, int dim, typename Number, typename OtherNumber>
-DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
-  Tensor<rank, dim, typename ProductType<Number, OtherNumber>::type>
-  operator+(const Tensor<rank, dim, Number> &     p,
+DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV inline DEAL_II_ALWAYS_INLINE
+                  Tensor<rank, dim, typename ProductType<Number, OtherNumber>::type>
+                  operator+(const Tensor<rank, dim, Number> &     p,
             const Tensor<rank, dim, OtherNumber> &q)
 {
   Tensor<rank, dim, typename ProductType<Number, OtherNumber>::type> tmp(p);
@@ -1735,12 +1805,14 @@ DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
  *
  * @tparam rank The rank of both tensors.
  *
+ * @note This function can also be used in CUDA device code.
+ *
  * @relatesalso Tensor
  */
 template <int rank, int dim, typename Number, typename OtherNumber>
-DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE
-  Tensor<rank, dim, typename ProductType<Number, OtherNumber>::type>
-  operator-(const Tensor<rank, dim, Number> &     p,
+DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV inline DEAL_II_ALWAYS_INLINE
+                  Tensor<rank, dim, typename ProductType<Number, OtherNumber>::type>
+                  operator-(const Tensor<rank, dim, Number> &     p,
             const Tensor<rank, dim, OtherNumber> &q)
 {
   Tensor<rank, dim, typename ProductType<Number, OtherNumber>::type> tmp(p);

diff --git a/tests/cuda/cuda_tensor_01.cu b/tests/cuda/cuda_tensor_01.cu
index d530992a4d3c6a022e024e66d2f01e276266376a..cec4c58bd8ce1e7f5f0c22489dc7664bcd5aeda7 100644 (file)
--- a/tests/cuda/cuda_tensor_01.cu
+++ b/tests/cuda/cuda_tensor_01.cu
@@ -13,14 +13,10 @@
 //
 // ---------------------------------------------------------------------
 
-// Test operator[] and norm_square of cuda_tensor.
+// Test operator[], norm and norm_square of cuda_tensor.
 
-#include <deal.II/base/logstream.h>
 #include <deal.II/base/tensor.h>
 
-#include <fstream>
-#include <iomanip>
-
 #include "../tests.h"
 
 void
@@ -33,16 +29,14 @@ test_cpu()
     for (unsigned int j = 0; j < dim; ++j)
       t[i][j] = a[i][j];
 
-
   deallog.push("values");
   for (unsigned int i = 0; i < dim; ++i)
     for (unsigned int j = 0; j < dim; ++j)
       deallog << t[i][j] << std::endl;
   deallog.pop();
 
-  deallog.push("norm_square");
-  deallog << t.norm_square() << std::endl;
-  deallog.pop();
+  deallog << "norm: " << t.norm() << std::endl;
+  deallog << "norm_square: " << t.norm_square() << std::endl;
 }
 
 __global__ void init_kernel(Tensor<2, 3> *t, const unsigned int N)
@@ -53,10 +47,13 @@ __global__ void init_kernel(Tensor<2, 3> *t, const unsigned int N)
     (*t)[i][j] = j + i * N + 1.;
 }
 
-__global__ void norm_kernel(Tensor<2, 3> *t, double *norm)
+__global__ void norm_kernel(Tensor<2, 3> *t, double *norm, double *norm_square)
 {
   if (threadIdx.x == 0)
-    *norm = t->norm_square();
+    {
+      *norm        = t->norm();
+      *norm_square = t->norm_square();
+    }
 }
 
 void
@@ -65,6 +62,8 @@ test_gpu()
   const unsigned int dim = 3;
   double *           norm_dev;
   double             norm_host;
+  double *           norm_square_dev;
+  double             norm_square_host;
   Tensor<2, dim> *   t_dev;
 
   // Allocate objects on the device
@@ -72,34 +71,41 @@ test_gpu()
   AssertCuda(cuda_error);
   cuda_error = cudaMalloc(&norm_dev, sizeof(double));
   AssertCuda(cuda_error);
+  cuda_error = cudaMalloc(&norm_square_dev, sizeof(double));
+  AssertCuda(cuda_error);
 
   // Launch the kernels.
   dim3 block_dim(dim, dim);
   init_kernel<<<1, block_dim>>>(t_dev, dim);
-  norm_kernel<<<1, 1>>>(t_dev, norm_dev);
+  norm_kernel<<<1, 1>>>(t_dev, norm_dev, norm_square_dev);
 
   // Copy the result to the device
   cuda_error =
     cudaMemcpy(&norm_host, norm_dev, sizeof(double), cudaMemcpyDeviceToHost);
   AssertCuda(cuda_error);
+  cuda_error = cudaMemcpy(&norm_square_host,
+                          norm_square_dev,
+                          sizeof(double),
+                          cudaMemcpyDeviceToHost);
+  AssertCuda(cuda_error);
 
   // Free memory
   cuda_error = cudaFree(t_dev);
   AssertCuda(cuda_error);
   cuda_error = cudaFree(norm_dev);
   AssertCuda(cuda_error);
+  cuda_error = cudaFree(norm_square_dev);
+  AssertCuda(cuda_error);
 
   // Output result
-  deallog.push("norm_square GPU");
-  deallog << norm_host << std::endl;
+  deallog << "norm GPU: " << norm_host << std::endl;
+  deallog << "norm_square GPU: " << norm_square_host << std::endl;
 }
 
 int
 main()
 {
-  std::ofstream logfile("output");
-  deallog << std::setprecision(5);
-  deallog.attach(logfile);
+  initlog();
 
   init_cuda();
 

diff --git a/tests/cuda/cuda_tensor_01.output b/tests/cuda/cuda_tensor_01.output
index 6b9b206f7d6a76c09152e18c5cb2a9e5ec0c93b8..1a8064704064f8837e806a70447b617ab4dbd6de 100644 (file)
--- a/tests/cuda/cuda_tensor_01.output
+++ b/tests/cuda/cuda_tensor_01.output
@@ -1,12 +1,14 @@
 
-DEAL:values::1.0000
-DEAL:values::2.0000
-DEAL:values::3.0000
-DEAL:values::4.0000
-DEAL:values::5.0000
-DEAL:values::6.0000
-DEAL:values::7.0000
-DEAL:values::8.0000
-DEAL:values::9.0000
-DEAL:norm_square::285.00
-DEAL:norm_square GPU::285.00
+DEAL:values::1.00000
+DEAL:values::2.00000
+DEAL:values::3.00000
+DEAL:values::4.00000
+DEAL:values::5.00000
+DEAL:values::6.00000
+DEAL:values::7.00000
+DEAL:values::8.00000
+DEAL:values::9.00000
+DEAL::norm: 16.8819
+DEAL::norm_square: 285.000
+DEAL::norm GPU: 16.8819
+DEAL::norm_square GPU: 285.000

diff --git a/tests/cuda/cuda_tensor_02.cu b/tests/cuda/cuda_tensor_02.cu
new file mode 100644 (file)
index 0000000..83864d2
--- /dev/null
+++ b/tests/cuda/cuda_tensor_02.cu
@@ -0,0 +1,252 @@
+// ---------------------------------------------------------------------
+//
+// 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 operator[] and norm_square of cuda_tensor.
+
+#include <deal.II/base/tensor.h>
+
+#include "../tests.h"
+
+template <int rank, int dim, typename Number>
+__global__ void
+miscellaneous_kernel()
+{
+  // constructors
+  typename Tensor<rank, dim, Number>::array_type array{};
+  Tensor<rank, dim, Number>                      dummy_1(array);
+  Tensor<rank, dim, Number>                      dummy_2;
+  Tensor<rank, dim, Number>                      dummy_3 = dummy_2;
+
+  // access
+  Tensor<rank + 1, dim, Number> initializer_1;
+  const auto                    dummy_5 = initializer_1[0];
+
+  // assignment
+  dummy_2 = dummy_3;
+}
+
+template <int rank, int dim, typename Number>
+__global__ void
+summation_kernel(Tensor<rank, dim, Number> *t,
+                 Tensor<rank, dim, Number> *t1,
+                 Tensor<rank, dim, Number> *t2)
+{
+  *t2 += *t;
+  *t1 = *t1 + *t;
+}
+
+template <int rank, int dim, typename Number>
+__global__ void
+subtraction_kernel(Tensor<rank, dim, Number> *t,
+                   Tensor<rank, dim, Number> *t1,
+                   Tensor<rank, dim, Number> *t2)
+{
+  *t2 -= *t;
+  *t1 = *t1 - *t;
+}
+
+template <int rank, int dim, typename Number>
+__global__ void
+multiplication_kernel(Tensor<rank, dim, Number> *t,
+                      Tensor<rank, dim, Number> *t1,
+                      Tensor<rank, dim, Number> *t2)
+{
+  *t1 = *t * Number(2.);
+  *t2 = Number(2.) * *t;
+  *t *= 2.;
+}
+
+template <int rank, int dim, typename Number>
+__global__ void
+division_kernel(Tensor<rank, dim, Number> *t,
+                Tensor<rank, dim, Number> *t1,
+                Tensor<rank, dim, Number> *t2)
+{
+  *t1 = *t / Number(2.);
+  *t /= 2.;
+  *t2 = *t1;
+}
+
+template <int dim, typename Number>
+__global__ void init_kernel(Tensor<0, dim, Number> *t)
+{
+  if (threadIdx.x == 0)
+    *t = 1.;
+}
+
+template <int dim, typename Number>
+__global__ void init_kernel(Tensor<1, dim, Number> *t)
+{
+  const unsigned int i = threadIdx.x;
+  if (i < dim)
+    (*t)[i] = i + 1.;
+}
+
+template <int dim, typename Number>
+__global__ void init_kernel(Tensor<2, dim, Number> *t)
+{
+  const unsigned int i = threadIdx.y;
+  const unsigned int j = threadIdx.x;
+  if ((i < dim) && (j < dim))
+    (*t)[i][j] = j + i * dim + 1.;
+}
+
+
+template <int rank, int dim, typename Number>
+void
+test_gpu()
+{
+  const double tolerance = 1.e-8;
+
+  Tensor<rank, dim, Number> *t_dev;
+  Tensor<rank, dim, Number> *t1_dev;
+  Tensor<rank, dim, Number> *t2_dev;
+
+  Tensor<rank, dim, Number> t_host;
+  Tensor<rank, dim, Number> t1_host;
+  Tensor<rank, dim, Number> t2_host;
+
+  Tensor<rank, dim, Number> reference_host;
+
+  // Allocate objects on the device
+  cudaError_t cuda_error =
+    cudaMalloc(&t_dev, sizeof(Tensor<rank, dim, Number>));
+  AssertCuda(cuda_error);
+  cuda_error = cudaMalloc(&t1_dev, sizeof(Tensor<rank, dim, Number>));
+  AssertCuda(cuda_error);
+  cuda_error = cudaMalloc(&t2_dev, sizeof(Tensor<rank, dim, Number>));
+  AssertCuda(cuda_error);
+
+  // Initialize
+  dim3 block_dim(dim, dim);
+  init_kernel<<<1, block_dim>>>(t_dev);
+  cuda_error = cudaMemcpy(&reference_host,
+                          t_dev,
+                          sizeof(Tensor<rank, dim, Number>),
+                          cudaMemcpyDeviceToHost);
+  AssertCuda(cuda_error);
+
+  // Test multiplication.
+  multiplication_kernel<<<1, 1>>>(t_dev, t1_dev, t2_dev);
+
+  cuda_error = cudaMemcpy(&t_host,
+                          t_dev,
+                          sizeof(Tensor<rank, dim, Number>),
+                          cudaMemcpyDeviceToHost);
+  AssertCuda(cuda_error);
+  cuda_error = cudaMemcpy(&t1_host,
+                          t1_dev,
+                          sizeof(Tensor<rank, dim, Number>),
+                          cudaMemcpyDeviceToHost);
+  AssertCuda(cuda_error);
+  cuda_error = cudaMemcpy(&t2_host,
+                          t2_dev,
+                          sizeof(Tensor<rank, dim, Number>),
+                          cudaMemcpyDeviceToHost);
+  AssertCuda(cuda_error);
+
+  reference_host *= 2;
+  AssertThrow((t_host - reference_host).norm() < tolerance, ExcInternalError());
+  AssertThrow((t1_host - reference_host).norm() < tolerance,
+              ExcInternalError());
+  AssertThrow((t2_host - reference_host).norm() < tolerance,
+              ExcInternalError());
+
+  deallog << "multiplication OK" << std::endl;
+
+  // Test division.
+  division_kernel<<<1, 1>>>(t_dev, t1_dev, t2_dev);
+  cuda_error = cudaMemcpy(&t_host,
+                          t_dev,
+                          sizeof(Tensor<rank, dim, Number>),
+                          cudaMemcpyDeviceToHost);
+  AssertCuda(cuda_error);
+  cuda_error = cudaMemcpy(&t1_host,
+                          t1_dev,
+                          sizeof(Tensor<rank, dim, Number>),
+                          cudaMemcpyDeviceToHost);
+  AssertCuda(cuda_error);
+
+  reference_host /= 2.;
+  AssertThrow((t_host - reference_host).norm() < tolerance, ExcInternalError());
+  AssertThrow((t1_host - reference_host).norm() < tolerance,
+              ExcInternalError());
+
+  deallog << "division OK" << std::endl;
+
+  // Test summation
+  summation_kernel<<<1, 1>>>(t_dev, t1_dev, t2_dev);
+  cuda_error = cudaMemcpy(&t1_host,
+                          t1_dev,
+                          sizeof(Tensor<rank, dim, Number>),
+                          cudaMemcpyDeviceToHost);
+  AssertCuda(cuda_error);
+  cuda_error = cudaMemcpy(&t2_host,
+                          t2_dev,
+                          sizeof(Tensor<rank, dim, Number>),
+                          cudaMemcpyDeviceToHost);
+  AssertCuda(cuda_error);
+
+  reference_host *= 2.;
+  AssertThrow((t1_host - reference_host).norm() < tolerance,
+              ExcInternalError());
+  AssertThrow((t2_host - reference_host).norm() < tolerance,
+              ExcInternalError());
+
+
+  // Test subtraction
+  subtraction_kernel<<<1, 1>>>(t_dev, t1_dev, t2_dev);
+  cuda_error = cudaMemcpy(&t1_host,
+                          t1_dev,
+                          sizeof(Tensor<rank, dim, Number>),
+                          cudaMemcpyDeviceToHost);
+  AssertCuda(cuda_error);
+  cuda_error = cudaMemcpy(&t2_host,
+                          t2_dev,
+                          sizeof(Tensor<rank, dim, Number>),
+                          cudaMemcpyDeviceToHost);
+
+  reference_host /= 2.;
+  AssertThrow((t1_host - reference_host).norm() < tolerance,
+              ExcInternalError());
+  AssertThrow((t2_host - reference_host).norm() < tolerance,
+              ExcInternalError());
+
+  // Miscellaneous
+  miscellaneous_kernel<rank, dim, Number><<<1, 1>>>();
+
+  // Free memory
+  cuda_error = cudaFree(t_dev);
+  AssertCuda(cuda_error);
+  cuda_error = cudaFree(t1_dev);
+  AssertCuda(cuda_error);
+  cuda_error = cudaFree(t2_dev);
+  AssertCuda(cuda_error);
+}
+
+int
+main()
+{
+  initlog();
+
+  init_cuda();
+
+  test_gpu<0, 3, double>();
+  test_gpu<1, 3, double>();
+  test_gpu<2, 3, double>();
+  test_gpu<0, 3, float>();
+  test_gpu<1, 3, float>();
+  test_gpu<2, 3, float>();
+}

diff --git a/tests/cuda/cuda_tensor_02.output b/tests/cuda/cuda_tensor_02.output
new file mode 100644 (file)
index 0000000..bd1f4f4
--- /dev/null
+++ b/tests/cuda/cuda_tensor_02.output
@@ -0,0 +1,13 @@
+
+DEAL::multiplication OK
+DEAL::division OK
+DEAL::multiplication OK
+DEAL::division OK
+DEAL::multiplication OK
+DEAL::division OK
+DEAL::multiplication OK
+DEAL::division OK
+DEAL::multiplication OK
+DEAL::division OK
+DEAL::multiplication OK
+DEAL::division OK