From 58fbca3ac185a5d57d0044755aab19ef58f9cfa4 Mon Sep 17 00:00:00 2001 From: Daniel Arndt Date: Mon, 1 Jul 2019 17:30:01 -0400 Subject: [PATCH] Annotate Tensor functions that can be used in CUDA device code --- include/deal.II/base/tensor.h | 242 ++++++++++++++++++----------- tests/cuda/cuda_tensor_01.cu | 40 ++--- tests/cuda/cuda_tensor_01.output | 24 +-- tests/cuda/cuda_tensor_02.cu | 252 +++++++++++++++++++++++++++++++ tests/cuda/cuda_tensor_02.output | 13 ++ 5 files changed, 458 insertions(+), 113 deletions(-) create mode 100644 tests/cuda/cuda_tensor_02.cu create mode 100644 tests/cuda/cuda_tensor_02.output diff --git a/include/deal.II/base/tensor.h b/include/deal.II/base/tensor.h index d9f8f4c676..961ceb7a18 100644 --- 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 - 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 - 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 - 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 - 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 - 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 - 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 factor. @@ -273,16 +289,20 @@ public: /** * Divide the scalar by factor. + * + * @note This function can also be used in CUDA device code. */ template - 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 l2 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 - constexpr Tensor(const Tensor &initializer); + constexpr DEAL_II_CUDA_HOST_DEV + Tensor(const Tensor &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 - DEAL_II_CONSTEXPR Tensor & - operator=(const Tensor &rhs); + DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV Tensor & + operator=(const Tensor &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 - DEAL_II_CONSTEXPR Tensor & - operator+=(const Tensor &); + DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV Tensor & + operator+=(const Tensor &); /** * Subtract another tensor. + * + * @note This function can also be used in CUDA device code. */ template - DEAL_II_CONSTEXPR Tensor & - operator-=(const Tensor &); + DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV Tensor & + operator-=(const Tensor &); /** * Scale the tensor by factor, i.e. multiply all components by @@ -578,16 +612,20 @@ public: /** * Scale the vector by 1/factor. + * + * @note This function can also be used in CUDA device code. */ template - 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 l2 norm of the vector. + * + * @note This function can also be used in CUDA device code. */ - DEAL_II_CONSTEXPR typename numbers::NumberTraits::real_type - norm() const; + DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV + typename numbers::NumberTraits::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 + * to create constexpr constructors for Tensor + * + * @note This function can also be used in CUDA device code. */ template - constexpr Tensor(const ArrayLike &initializer, - std_cxx14::index_sequence); + constexpr DEAL_II_CUDA_HOST_DEV + Tensor(const ArrayLike &initializer, std_cxx14::index_sequence); /** * Allow an arbitrary Tensor to access the underlying values. @@ -783,7 +826,8 @@ Tensor<0, dim, Number>::Tensor() template template -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::value(initializer)) {} @@ -791,7 +835,8 @@ constexpr Tensor<0, dim, Number>::Tensor(const OtherNumber &initializer) template template -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 template -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::value(p); @@ -871,7 +916,7 @@ Tensor<0, dim, Number>::operator=(const Tensor<0, dim, OtherNumber> &p) #ifdef __INTEL_COMPILER template -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 template -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::value(d); @@ -895,7 +940,8 @@ template 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::value || std::is_same::value), ExcMessage( @@ -918,7 +964,7 @@ Tensor<0, dim, Number>::operator!=(const Tensor<0, dim, OtherNumber> &p) const template template -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 template -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 template -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 -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 template -DEAL_II_ALWAYS_INLINE constexpr Tensor::Tensor( - const ArrayLike &initializer, - std_cxx14::index_sequence) +DEAL_II_ALWAYS_INLINE constexpr DEAL_II_CUDA_HOST_DEV +Tensor::Tensor(const ArrayLike &initializer, + std_cxx14::index_sequence) : values{Tensor(initializer[indices])...} { static_assert(sizeof...(indices) == dim, @@ -1060,15 +1106,16 @@ DEAL_II_ALWAYS_INLINE constexpr Tensor::Tensor( template -DEAL_II_ALWAYS_INLINE constexpr Tensor::Tensor( - const array_type &initializer) +DEAL_II_ALWAYS_INLINE constexpr DEAL_II_CUDA_HOST_DEV +Tensor::Tensor(const array_type &initializer) : Tensor(initializer, std_cxx14::make_index_sequence{}) {} template template -DEAL_II_ALWAYS_INLINE constexpr Tensor::Tensor( +DEAL_II_ALWAYS_INLINE constexpr DEAL_II_CUDA_HOST_DEV +Tensor::Tensor( const Tensor &initializer) : Tensor(initializer, std_cxx14::make_index_sequence{}) {} @@ -1154,8 +1201,8 @@ DEAL_II_ALWAYS_INLINE constexpr DEAL_II_CUDA_HOST_DEV const typename Tensor< Number>::value_type &Tensor:: operator[](const unsigned int i) const { - return dealii::internal::TensorSubscriptor::subscript( - values, i, std::integral_constant()); + return values[i]; /*dealii::internal::TensorSubscriptor::subscript( + values, i, std::integral_constant());*/ } @@ -1285,7 +1332,7 @@ operator!=(const Tensor &p) const template template -DEAL_II_CONSTEXPR inline Tensor & +DEAL_II_CONSTEXPR inline DEAL_II_CUDA_HOST_DEV Tensor & Tensor::operator+=(const Tensor &p) { for (unsigned int i = 0; i < dim; ++i) @@ -1296,7 +1343,7 @@ Tensor::operator+=(const Tensor &p) template template -DEAL_II_CONSTEXPR inline Tensor & +DEAL_II_CONSTEXPR inline DEAL_II_CUDA_HOST_DEV Tensor & Tensor::operator-=(const Tensor &p) { for (unsigned int i = 0; i < dim; ++i) @@ -1318,7 +1365,7 @@ Tensor::operator*=(const OtherNumber &s) template template -DEAL_II_CONSTEXPR inline Tensor & +DEAL_II_CONSTEXPR inline DEAL_II_CUDA_HOST_DEV Tensor & Tensor::operator/=(const OtherNumber &s) { for (unsigned int i = 0; i < dim; ++i) @@ -1328,7 +1375,7 @@ Tensor::operator/=(const OtherNumber &s) template -DEAL_II_CONSTEXPR inline Tensor +DEAL_II_CONSTEXPR inline DEAL_II_CUDA_HOST_DEV Tensor Tensor::operator-() const { Tensor 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 -constexpr DEAL_II_ALWAYS_INLINE typename ProductType::type -operator*(const Other &object, const Tensor<0, dim, Number> &t) +constexpr DEAL_II_CUDA_HOST_DEV DEAL_II_ALWAYS_INLINE + typename ProductType::type + operator*(const Other &object, const Tensor<0, dim, Number> &t) { return object * static_cast(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 -constexpr DEAL_II_ALWAYS_INLINE typename ProductType::type -operator*(const Tensor<0, dim, Number> &t, const Other &object) +constexpr DEAL_II_CUDA_HOST_DEV DEAL_II_ALWAYS_INLINE + typename ProductType::type + operator*(const Tensor<0, dim, Number> &t, const Other &object) { return static_cast(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 -constexpr DEAL_II_ALWAYS_INLINE typename ProductType::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::type + operator*(const Tensor<0, dim, Number> & src1, + const Tensor<0, dim, OtherNumber> &src2) { return static_cast(src1) * static_cast(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 -constexpr DEAL_II_ALWAYS_INLINE +DEAL_II_CUDA_HOST_DEV constexpr DEAL_II_ALWAYS_INLINE Tensor<0, dim, typename ProductType */ template -constexpr DEAL_II_ALWAYS_INLINE - Tensor<0, dim, typename ProductType::type> - operator+(const Tensor<0, dim, Number> & p, +constexpr DEAL_II_ALWAYS_INLINE DEAL_II_CUDA_HOST_DEV + Tensor<0, dim, typename ProductType::type> + operator+(const Tensor<0, dim, Number> & p, const Tensor<0, dim, OtherNumber> &q) { return static_cast(p) + static_cast(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 -constexpr DEAL_II_ALWAYS_INLINE - Tensor<0, dim, typename ProductType::type> - operator-(const Tensor<0, dim, Number> & p, +constexpr DEAL_II_ALWAYS_INLINE DEAL_II_CUDA_HOST_DEV + Tensor<0, dim, typename ProductType::type> + operator-(const Tensor<0, dim, Number> & p, const Tensor<0, dim, OtherNumber> &q) { return static_cast(p) - static_cast(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 -DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE - Tensor::type>::type> - operator*(const Tensor &t, const OtherNumber &factor) + operator*(const Tensor &t, const OtherNumber &factor) { // recurse over the base objects Tensor::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 -constexpr DEAL_II_ALWAYS_INLINE +DEAL_II_CUDA_HOST_DEV constexpr DEAL_II_ALWAYS_INLINE Tensor::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 -DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE - Tensor::type>::type> - operator/(const Tensor &t, const OtherNumber &factor) + operator/(const Tensor &t, const OtherNumber &factor) { // recurse over the base objects Tensor::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 -DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE - Tensor::type> - operator+(const Tensor & p, +DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV inline DEAL_II_ALWAYS_INLINE + Tensor::type> + operator+(const Tensor & p, const Tensor &q) { Tensor::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 -DEAL_II_CONSTEXPR inline DEAL_II_ALWAYS_INLINE - Tensor::type> - operator-(const Tensor & p, +DEAL_II_CONSTEXPR DEAL_II_CUDA_HOST_DEV inline DEAL_II_ALWAYS_INLINE + Tensor::type> + operator-(const Tensor & p, const Tensor &q) { Tensor::type> tmp(p); diff --git a/tests/cuda/cuda_tensor_01.cu b/tests/cuda/cuda_tensor_01.cu index d530992a4d..cec4c58bd8 100644 --- 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 #include -#include -#include - #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 6b9b206f7d..1a80647040 100644 --- 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 index 0000000000..83864d20f1 --- /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 + +#include "../tests.h" + +template +__global__ void +miscellaneous_kernel() +{ + // constructors + typename Tensor::array_type array{}; + Tensor dummy_1(array); + Tensor dummy_2; + Tensor dummy_3 = dummy_2; + + // access + Tensor initializer_1; + const auto dummy_5 = initializer_1[0]; + + // assignment + dummy_2 = dummy_3; +} + +template +__global__ void +summation_kernel(Tensor *t, + Tensor *t1, + Tensor *t2) +{ + *t2 += *t; + *t1 = *t1 + *t; +} + +template +__global__ void +subtraction_kernel(Tensor *t, + Tensor *t1, + Tensor *t2) +{ + *t2 -= *t; + *t1 = *t1 - *t; +} + +template +__global__ void +multiplication_kernel(Tensor *t, + Tensor *t1, + Tensor *t2) +{ + *t1 = *t * Number(2.); + *t2 = Number(2.) * *t; + *t *= 2.; +} + +template +__global__ void +division_kernel(Tensor *t, + Tensor *t1, + Tensor *t2) +{ + *t1 = *t / Number(2.); + *t /= 2.; + *t2 = *t1; +} + +template +__global__ void init_kernel(Tensor<0, dim, Number> *t) +{ + if (threadIdx.x == 0) + *t = 1.; +} + +template +__global__ void init_kernel(Tensor<1, dim, Number> *t) +{ + const unsigned int i = threadIdx.x; + if (i < dim) + (*t)[i] = i + 1.; +} + +template +__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 +void +test_gpu() +{ + const double tolerance = 1.e-8; + + Tensor *t_dev; + Tensor *t1_dev; + Tensor *t2_dev; + + Tensor t_host; + Tensor t1_host; + Tensor t2_host; + + Tensor reference_host; + + // Allocate objects on the device + cudaError_t cuda_error = + cudaMalloc(&t_dev, sizeof(Tensor)); + AssertCuda(cuda_error); + cuda_error = cudaMalloc(&t1_dev, sizeof(Tensor)); + AssertCuda(cuda_error); + cuda_error = cudaMalloc(&t2_dev, sizeof(Tensor)); + 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), + 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), + cudaMemcpyDeviceToHost); + AssertCuda(cuda_error); + cuda_error = cudaMemcpy(&t1_host, + t1_dev, + sizeof(Tensor), + cudaMemcpyDeviceToHost); + AssertCuda(cuda_error); + cuda_error = cudaMemcpy(&t2_host, + t2_dev, + sizeof(Tensor), + 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), + cudaMemcpyDeviceToHost); + AssertCuda(cuda_error); + cuda_error = cudaMemcpy(&t1_host, + t1_dev, + sizeof(Tensor), + 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), + cudaMemcpyDeviceToHost); + AssertCuda(cuda_error); + cuda_error = cudaMemcpy(&t2_host, + t2_dev, + sizeof(Tensor), + 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), + cudaMemcpyDeviceToHost); + AssertCuda(cuda_error); + cuda_error = cudaMemcpy(&t2_host, + t2_dev, + sizeof(Tensor), + cudaMemcpyDeviceToHost); + + reference_host /= 2.; + AssertThrow((t1_host - reference_host).norm() < tolerance, + ExcInternalError()); + AssertThrow((t2_host - reference_host).norm() < tolerance, + ExcInternalError()); + + // Miscellaneous + miscellaneous_kernel<<<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 index 0000000000..bd1f4f4e12 --- /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 -- 2.39.5