From: Bruno Turcksin Date: Sun, 14 May 2017 03:00:33 +0000 (-0400) Subject: Add CUDA support for tensor. X-Git-Tag: v9.0.0-rc1~1588^2~5 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=6127759382b092ff34368864119092bf25898aca;p=dealii.git Add CUDA support for tensor. --- diff --git a/include/deal.II/base/numbers.h b/include/deal.II/base/numbers.h index d480e3b963..6e3782926b 100644 --- a/include/deal.II/base/numbers.h +++ b/include/deal.II/base/numbers.h @@ -24,6 +24,13 @@ #include #include +#ifdef DEAL_II_WITH_CUDA +# include +# define DEAL_II_CUDA_HOST_DEV __host__ __device__ +#else +# define DEAL_II_CUDA_HOST_DEV +#endif + DEAL_II_NAMESPACE_OPEN // forward declarations to support abs or sqrt operations on VectorizedArray @@ -210,8 +217,11 @@ namespace numbers * Return the square of the absolute value of the given number. Since the * general template is chosen for types not equal to std::complex, this * function simply returns the square of the given number. + * + * @ingroup CUDAWrappers */ static + DEAL_II_CUDA_HOST_DEV real_type abs_square (const number &x); /** @@ -341,6 +351,7 @@ namespace numbers template + DEAL_II_CUDA_HOST_DEV typename NumberTraits::real_type NumberTraits::abs_square (const number &x) { diff --git a/include/deal.II/base/tensor.h b/include/deal.II/base/tensor.h index 39ba02107e..58748bae72 100644 --- a/include/deal.II/base/tensor.h +++ b/include/deal.II/base/tensor.h @@ -27,6 +27,7 @@ #include #include + DEAL_II_NAMESPACE_OPEN // Forward declarations: @@ -129,8 +130,10 @@ public: /** * Constructor. Set to zero. + * + * @ingroup CUDAWrappers */ - Tensor (); + DEAL_II_CUDA_HOST_DEV Tensor (); /** * Constructor from tensors with different underlying scalar type. This @@ -152,16 +155,20 @@ public: * * This is the non-const conversion operator that returns a writable * reference. + * + * @ingroup CUDAWrappers */ - operator Number &(); + DEAL_II_CUDA_HOST_DEV operator Number &(); /** * Return a reference to the encapsulated Number object. Since rank-0 * tensors are scalars, this is a natural operation. * * This is the const conversion operator that returns a read-only reference. + * + * @ingroup CUDAWrappers */ - operator const Number &() const; + DEAL_II_CUDA_HOST_DEV operator const Number &() const; /** * Assignment from tensors with different underlying scalar type. This @@ -197,9 +204,11 @@ public: /** * Multiply the scalar with a factor. + * + * @ingroup CUDAWrappers */ template - Tensor<0,dim,Number> &operator *= (const OtherNumber factor); + DEAL_II_CUDA_HOST_DEV Tensor<0,dim,Number> &operator *= (const OtherNumber factor); /** * Divide the scalar by factor. @@ -236,8 +245,10 @@ public: /** * Return the square of the Frobenius-norm of a tensor, i.e. the sum of the * absolute squares of all entries. + * + * @ingroup CUDAWrappers */ - real_type norm_square () const; + DEAL_II_CUDA_HOST_DEV real_type norm_square () const; /** * Read or write the data of this object to or from a stream for the purpose @@ -358,8 +369,10 @@ public: /** * Constructor. Initialize all entries to zero. + * + * @ingroup CUDAWrappers */ - Tensor (); + DEAL_II_CUDA_HOST_DEV Tensor (); /** * Constructor, where the data is copied from a C-style array. @@ -388,13 +401,17 @@ public: /** * Read-Write access operator. + * + * @ingroup CUDAWrappers */ - value_type &operator [] (const unsigned int i); + DEAL_II_CUDA_HOST_DEV value_type &operator [] (const unsigned int i); /** * Read-only access operator. + * + * @ingroup CUDAWrappers */ - const value_type &operator[](const unsigned int i) const; + DEAL_II_CUDA_HOST_DEV const value_type &operator[](const unsigned int i) const; /** * Read access using TableIndices indices @@ -449,9 +466,11 @@ public: /** * Scale the tensor by factor, i.e. multiply all components by * factor. + * + * @ingroup CUDAWrappers */ template - Tensor &operator *= (const OtherNumber factor); + DEAL_II_CUDA_HOST_DEV Tensor &operator *= (const OtherNumber factor); /** * Scale the vector by 1/factor. @@ -489,8 +508,10 @@ public: /** * Return the square of the Frobenius-norm of a tensor, i.e. the sum of the * absolute squares of all entries. + * + * @ingroup CUDAWrappers */ - typename numbers::NumberTraits::real_type norm_square() const; + DEAL_II_CUDA_HOST_DEV typename numbers::NumberTraits::real_type norm_square() const; /** * Fill a vector with all tensor elements. @@ -569,7 +590,7 @@ private: template inline -Tensor<0,dim,Number>::Tensor () +DEAL_II_CUDA_HOST_DEV Tensor<0,dim,Number>::Tensor () : value() { } @@ -595,7 +616,7 @@ Tensor<0,dim,Number>::Tensor (const Tensor<0,dim,OtherNumber> &p) template inline -Tensor<0,dim,Number>::operator Number &() +DEAL_II_CUDA_HOST_DEV Tensor<0,dim,Number>::operator Number &() { Assert(dim != 0, ExcMessage("Cannot access an object of type Tensor<0,0,Number>")); return value; @@ -604,7 +625,7 @@ Tensor<0,dim,Number>::operator Number &() template inline -Tensor<0,dim,Number>::operator const Number &() const +DEAL_II_CUDA_HOST_DEV Tensor<0,dim,Number>::operator const Number &() const { Assert(dim != 0, ExcMessage("Cannot access an object of type Tensor<0,0,Number>")); return value; @@ -662,7 +683,7 @@ Tensor<0,dim,Number> &Tensor<0,dim,Number>::operator -= (const Tensor<0,dim,Othe template template inline -Tensor<0,dim,Number> &Tensor<0,dim,Number>::operator *= (const OtherNumber s) +DEAL_II_CUDA_HOST_DEV Tensor<0,dim,Number> &Tensor<0,dim,Number>::operator *= (const OtherNumber s) { value *= s; return *this; @@ -700,7 +721,7 @@ Tensor<0,dim,Number>::norm () const template inline typename Tensor<0,dim,Number>::real_type -Tensor<0,dim,Number>::norm_square () const +DEAL_II_CUDA_HOST_DEV Tensor<0,dim,Number>::norm_square () const { Assert(dim != 0, ExcMessage("Cannot access an object of type Tensor<0,0,Number>")); return numbers::NumberTraits::abs_square (value); @@ -742,7 +763,7 @@ void Tensor<0,dim,Number>::serialize(Archive &ar, const unsigned int) template inline -Tensor::Tensor () +DEAL_II_CUDA_HOST_DEV Tensor::Tensor () { // All members of the c-style array values are already default initialized // and thus all values are already set to zero recursively. @@ -796,6 +817,7 @@ namespace internal { template inline DEAL_II_ALWAYS_INLINE + DEAL_II_CUDA_HOST_DEV ArrayElementType & subscript (ArrayElementType *values, const unsigned int i, @@ -807,6 +829,7 @@ namespace internal template + DEAL_II_CUDA_HOST_DEV ArrayElementType & subscript (ArrayElementType *, const unsigned int, @@ -822,6 +845,7 @@ namespace internal template inline DEAL_II_ALWAYS_INLINE +DEAL_II_CUDA_HOST_DEV typename Tensor::value_type & Tensor::operator[] (const unsigned int i) { @@ -831,6 +855,7 @@ Tensor::operator[] (const unsigned int i) template inline DEAL_II_ALWAYS_INLINE +DEAL_II_CUDA_HOST_DEV const typename Tensor::value_type & Tensor::operator[] (const unsigned int i) const { @@ -950,6 +975,7 @@ Tensor::operator -= (const Tensor &p) template template inline +DEAL_II_CUDA_HOST_DEV Tensor & Tensor::operator *= (const OtherNumber s) { @@ -996,6 +1022,7 @@ Tensor::norm () const template inline +DEAL_II_CUDA_HOST_DEV typename numbers::NumberTraits::real_type Tensor::norm_square () const { diff --git a/tests/cuda/cuda_tensor_01.cu b/tests/cuda/cuda_tensor_01.cu new file mode 100644 index 0000000000..b8264e0652 --- /dev/null +++ b/tests/cuda/cuda_tensor_01.cu @@ -0,0 +1,104 @@ +// --------------------------------------------------------------------- +// +// 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 at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +// Test operator[] and norm_square of cuda_tensor. + +#include "../tests.h" +#include +#include +#include +#include + +void test_cpu() +{ + double a[3][3] = {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}; + const unsigned int dim=3; + Tensor<2,dim> t; + for (unsigned int i=0; i *t, + const unsigned int N) +{ + const unsigned int i = threadIdx.y; + const unsigned int j = threadIdx.x; + if ((i < N) && (j < N)) + (*t)[i][j] = j + i*N + 1.; +} + +__global__ void norm_kernel(Tensor<2,3> *t, double *norm) +{ + if (threadIdx.x == 0) + *norm = t->norm_square(); +} + +void test_gpu() +{ + const unsigned int dim=3; + double *norm_dev; + double norm_host; + Tensor<2,dim> *t_dev; + + // Allocate objects on the device + cudaError_t cuda_error = cudaMalloc(&t_dev, sizeof(Tensor<2,dim>)); + AssertCuda(cuda_error); + cuda_error = cudaMalloc(&norm_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); + + // Copy the result to the device + cuda_error = cudaMemcpy(&norm_host, norm_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); + + // Output result + deallog.push("norm_square GPU"); + deallog << norm_host << std::endl; +} + +int main () +{ + std::ofstream logfile("output"); + deallog << std::setprecision(5); + deallog.attach(logfile); + deallog.threshold_double(1.e-10); + + test_cpu(); + + test_gpu(); +} diff --git a/tests/cuda/cuda_tensor_01.output b/tests/cuda/cuda_tensor_01.output new file mode 100644 index 0000000000..6b9b206f7d --- /dev/null +++ b/tests/cuda/cuda_tensor_01.output @@ -0,0 +1,12 @@ + +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