Add CUDA support for tensor.

diff --git a/include/deal.II/base/numbers.h b/include/deal.II/base/numbers.h
index d480e3b963407059dca71cc2201596eacb278d9d..6e3782926b14ecd0a51143f360a07c4cc76af9e6 100644 (file)
--- a/include/deal.II/base/numbers.h
+++ b/include/deal.II/base/numbers.h
@@ -24,6 +24,13 @@
 #include <cstdlib>
 #include <complex>
 
+#ifdef DEAL_II_WITH_CUDA
+#  include <cuda_runtime_api.h>
+#  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 <typename number>
+  DEAL_II_CUDA_HOST_DEV
   typename NumberTraits<number>::real_type
   NumberTraits<number>::abs_square (const number &x)
   {

diff --git a/include/deal.II/base/tensor.h b/include/deal.II/base/tensor.h
index 39ba02107ec56c92ce26940c8bad421dec2c317d..58748bae7284289da30a488e3f45cd31aa4ab981 100644 (file)
--- a/include/deal.II/base/tensor.h
+++ b/include/deal.II/base/tensor.h
@@ -27,6 +27,7 @@
 #include <ostream>
 #include <vector>
 
+
 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 <tt>factor</tt>.
+   *
+   * @ingroup CUDAWrappers
    */
   template<typename OtherNumber>
-  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 <tt>factor</tt>.
@@ -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 <tt>indices</tt>
@@ -449,9 +466,11 @@ public:
   /**
    * Scale the tensor by <tt>factor</tt>, i.e. multiply all components by
    * <tt>factor</tt>.
+   *
+   * @ingroup CUDAWrappers
    */
   template <typename OtherNumber>
-  Tensor<rank_,dim,Number> &operator *= (const OtherNumber factor);
+  DEAL_II_CUDA_HOST_DEV Tensor<rank_,dim,Number> &operator *= (const OtherNumber factor);
 
   /**
    * Scale the vector by <tt>1/factor</tt>.
@@ -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<Number>::real_type norm_square() const;
+  DEAL_II_CUDA_HOST_DEV typename numbers::NumberTraits<Number>::real_type norm_square() const;
 
   /**
    * Fill a vector with all tensor elements.
@@ -569,7 +590,7 @@ private:
 
 template <int dim,typename Number>
 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 <int dim, typename Number>
 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 <int dim, typename Number>
 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 <int dim, typename Number>
 template <typename OtherNumber>
 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 <int dim, typename Number>
 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<Number>::abs_square (value);
@@ -742,7 +763,7 @@ void Tensor<0,dim,Number>::serialize(Archive &ar, const unsigned int)
 
 template <int rank_, int dim, typename Number>
 inline
-Tensor<rank_,dim,Number>::Tensor ()
+DEAL_II_CUDA_HOST_DEV Tensor<rank_,dim,Number>::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 <typename ArrayElementType, int dim>
     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 <typename ArrayElementType>
+    DEAL_II_CUDA_HOST_DEV
     ArrayElementType &
     subscript (ArrayElementType *,
                const unsigned int,
@@ -822,6 +845,7 @@ namespace internal
 
 template <int rank_, int dim, typename Number>
 inline DEAL_II_ALWAYS_INLINE
+DEAL_II_CUDA_HOST_DEV
 typename Tensor<rank_,dim,Number>::value_type &
 Tensor<rank_,dim,Number>::operator[] (const unsigned int i)
 {
@@ -831,6 +855,7 @@ Tensor<rank_,dim,Number>::operator[] (const unsigned int i)
 
 template <int rank_, int dim, typename Number>
 inline DEAL_II_ALWAYS_INLINE
+DEAL_II_CUDA_HOST_DEV
 const typename Tensor<rank_,dim,Number>::value_type &
 Tensor<rank_,dim,Number>::operator[] (const unsigned int i) const
 {
@@ -950,6 +975,7 @@ Tensor<rank_,dim,Number>::operator -= (const Tensor<rank_,dim,OtherNumber> &p)
 template <int rank_, int dim, typename Number>
 template <typename OtherNumber>
 inline
+DEAL_II_CUDA_HOST_DEV
 Tensor<rank_,dim,Number> &
 Tensor<rank_,dim,Number>::operator *= (const OtherNumber s)
 {
@@ -996,6 +1022,7 @@ Tensor<rank_,dim,Number>::norm () const
 
 template <int rank_, int dim, typename Number>
 inline
+DEAL_II_CUDA_HOST_DEV
 typename numbers::NumberTraits<Number>::real_type
 Tensor<rank_,dim,Number>::norm_square () const
 {

diff --git a/tests/cuda/cuda_tensor_01.cu b/tests/cuda/cuda_tensor_01.cu
new file mode 100644 (file)
index 0000000..b8264e0
--- /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 <deal.II/base/tensor.h>
+#include <deal.II/base/logstream.h>
+#include <fstream>
+#include <iomanip>
+
+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<dim; ++i)
+    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();
+}
+
+__global__ void init_kernel(Tensor<2,3> *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 (file)
index 0000000..6b9b206
--- /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