From 8e3efa78882a542c772ef96f5a3ee774d1774dc3 Mon Sep 17 00:00:00 2001 From: Daniel Arndt Date: Fri, 5 Jul 2019 14:26:34 -0400 Subject: [PATCH] Annotate Point usable for CUDA kernels --- doc/news/changes/minor/20190705DanielArndt | 4 + include/deal.II/base/point.h | 166 +++++++++++++-------- tests/cuda/cuda_point.cu | 80 ++++++++++ tests/cuda/cuda_point.output | 7 + 4 files changed, 198 insertions(+), 59 deletions(-) create mode 100644 doc/news/changes/minor/20190705DanielArndt create mode 100644 tests/cuda/cuda_point.cu create mode 100644 tests/cuda/cuda_point.output 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. +
+(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 dim==1 since the usage is considered unsafe for points with * dim!=1 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: * dim!=2 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 dim!=3 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 i, i.e., a vector * that is zero in all coordinates except for a single 1 in the ith * coordinate. + * + * @note This function can also be used in CUDA device code. */ - static Point - unit_vector(const unsigned int i); + static DEAL_II_CUDA_HOST_DEV Point + unit_vector(const unsigned int i); /** * Read access to the indexth 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 indexth 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 - operator+(const Tensor<1, dim, Number> &) const; + DEAL_II_CUDA_HOST_DEV Point + 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@. + * + * @note This function can also be used in CUDA device code. */ - Tensor<1, dim, Number> - operator-(const Point &) const; + DEAL_II_CUDA_HOST_DEV Tensor<1, dim, Number> + operator-(const Point &) 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@ object rather than as a Tensor@<1,dim@>. + * + * @note This function can also be used in CUDA device code. */ - Point - operator-(const Tensor<1, dim, Number> &) const; + DEAL_II_CUDA_HOST_DEV Point + operator-(const Tensor<1, dim, Number> &) const; /** * The opposite vector. + * + * @note This function can also be used in CUDA device code. */ - Point - operator-() const; + DEAL_II_CUDA_HOST_DEV Point + 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 - Point::type>::type> + DEAL_II_CUDA_HOST_DEV Point< + dim, + typename ProductType::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 - Point::type>::type> + DEAL_II_CUDA_HOST_DEV Point< + dim, + typename ProductType::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::norm_square() which returns the square of the * Frobenius norm. + * + * @note This function can also be used in CUDA device code. */ - typename numbers::NumberTraits::real_type + DEAL_II_CUDA_HOST_DEV typename numbers::NumberTraits::real_type square() const; /** * Return the Euclidean distance of this point to the point * p, i.e. the l_2 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::real_type + DEAL_II_CUDA_HOST_DEV typename numbers::NumberTraits::real_type distance(const Point &p) const; /** * Return the squared Euclidean distance of this point to the point * p. + * + * @note This function can also be used in CUDA device code. */ - typename numbers::NumberTraits::real_type + DEAL_II_CUDA_HOST_DEV typename numbers::NumberTraits::real_type distance_square(const Point &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::Point () = default' here. template -inline Point::Point() // NOLINT +inline DEAL_II_CUDA_HOST_DEV +Point::Point() // NOLINT {} template -inline Point::Point(const Tensor<1, dim, Number> &t) +inline DEAL_II_CUDA_HOST_DEV +Point::Point(const Tensor<1, dim, Number> &t) : Tensor<1, dim, Number>(t) {} template -inline Point::Point(const Number x) +inline DEAL_II_CUDA_HOST_DEV +Point::Point(const Number x) { Assert(dim == 1, ExcMessage( @@ -332,7 +376,8 @@ inline Point::Point(const Number x) template -inline Point::Point(const Number x, const Number y) +inline DEAL_II_CUDA_HOST_DEV +Point::Point(const Number x, const Number y) { Assert(dim == 2, ExcMessage( @@ -351,7 +396,8 @@ inline Point::Point(const Number x, const Number y) template -inline Point::Point(const Number x, const Number y, const Number z) +inline DEAL_II_CUDA_HOST_DEV +Point::Point(const Number x, const Number y, const Number z) { Assert(dim == 3, ExcMessage( @@ -395,8 +441,8 @@ inline Point::Point( template -inline Point -Point::unit_vector(unsigned int i) +inline DEAL_II_CUDA_HOST_DEV Point + Point::unit_vector(unsigned int i) { Point p; p[i] = 1.; @@ -405,7 +451,7 @@ Point::unit_vector(unsigned int i) template -inline Number +inline DEAL_II_CUDA_HOST_DEV Number Point::operator()(const unsigned int index) const { AssertIndexRange(index, dim); @@ -415,7 +461,7 @@ Point::operator()(const unsigned int index) const template -inline Number & +inline DEAL_II_CUDA_HOST_DEV Number & Point::operator()(const unsigned int index) { AssertIndexRange(index, dim); @@ -425,7 +471,7 @@ Point::operator()(const unsigned int index) template -inline Point +inline DEAL_II_CUDA_HOST_DEV Point Point::operator+(const Tensor<1, dim, Number> &p) const { Point tmp = *this; @@ -436,7 +482,7 @@ Point::operator+(const Tensor<1, dim, Number> &p) const template -inline Tensor<1, dim, Number> +inline DEAL_II_CUDA_HOST_DEV Tensor<1, dim, Number> Point::operator-(const Point &p) const { return (Tensor<1, dim, Number>(*this) -= p); @@ -445,7 +491,7 @@ Point::operator-(const Point &p) const template -inline Point +inline DEAL_II_CUDA_HOST_DEV Point Point::operator-(const Tensor<1, dim, Number> &p) const { Point tmp = *this; @@ -456,7 +502,7 @@ Point::operator-(const Tensor<1, dim, Number> &p) const template -inline Point +inline DEAL_II_CUDA_HOST_DEV Point Point::operator-() const { Point result; @@ -469,11 +515,11 @@ Point::operator-() const template template -inline Point< - dim, - typename ProductType::type>::type> - Point::operator*(const OtherNumber factor) const +inline DEAL_II_CUDA_HOST_DEV + Point::type>::type> + Point::operator*(const OtherNumber factor) const { Point::type> tmp; for (unsigned int i = 0; i < dim; ++i) @@ -485,11 +531,11 @@ inline Point< template template -inline Point< - dim, - typename ProductType::type>::type> -Point::operator/(const OtherNumber factor) const +inline DEAL_II_CUDA_HOST_DEV + Point::type>::type> + Point::operator/(const OtherNumber factor) const { Point::type> tmp; for (unsigned int i = 0; i < dim; ++i) @@ -500,8 +546,8 @@ Point::operator/(const OtherNumber factor) const template -inline Number Point:: - operator*(const Tensor<1, dim, Number> &p) const +inline DEAL_II_CUDA_HOST_DEV Number Point:: + 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:: template -inline typename numbers::NumberTraits::real_type +inline DEAL_II_CUDA_HOST_DEV typename numbers::NumberTraits::real_type Point::square() const { return this->norm_square(); @@ -520,7 +566,7 @@ Point::square() const template -inline typename numbers::NumberTraits::real_type +inline DEAL_II_CUDA_HOST_DEV typename numbers::NumberTraits::real_type Point::distance(const Point &p) const { return std::sqrt(distance_square(p)); @@ -529,7 +575,7 @@ Point::distance(const Point &p) const template -inline typename numbers::NumberTraits::real_type +inline DEAL_II_CUDA_HOST_DEV typename numbers::NumberTraits::real_type Point::distance_square(const Point &p) const { Number sum = internal::NumberType::value(0.0); @@ -563,15 +609,17 @@ Point::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 -inline Point< - dim, - typename ProductType::type>::type> -operator*(const OtherNumber factor, const Point &p) +inline DEAL_II_CUDA_HOST_DEV + Point::type>::type> + operator*(const OtherNumber factor, const Point &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 + +#include "../tests.h" + +template +__global__ void +miscellaneous_kernel() +{ + Point p_1; + Point p_2(Tensor<1, dim, Number>{}); + if (dim == 1) + Point p(1.); + if (dim == 2) + Point p(1., 2.); + if (dim == 3) + Point p(1., 2., 3.); + + auto p_3 = Point::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 +void +test_gpu() +{ + // Miscellaneous + miscellaneous_kernel<<<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 -- 2.39.5