From: Bruno Turcksin <bruno.turcksin@gmail.com>
Date: Mon, 31 Oct 2016 19:15:09 +0000 (-0400)
Subject: Add mean_value to cuda_vector.
X-Git-Tag: v8.5.0-rc1~514^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=1624c99f40004ff3e5772a416d2863dbb04a4225;p=dealii.git

Add mean_value to cuda_vector.
---

diff --git a/include/deal.II/lac/cuda_vector.h b/include/deal.II/lac/cuda_vector.h
index 8a925f6346..c91751d486 100644
--- a/include/deal.II/lac/cuda_vector.h
+++ b/include/deal.II/lac/cuda_vector.h
@@ -44,6 +44,7 @@ namespace LinearAlgebra
     class Vector: public VectorSpaceVector<Number>
     {
     public:
+      typedef typename VectorSpaceVector<Number>::value_type value_type;
       typedef typename VectorSpaceVector<Number>::size_type  size_type;
       typedef typename VectorSpaceVector<Number>::real_type  real_type;
 
@@ -155,6 +156,11 @@ namespace LinearAlgebra
        */
       virtual void equ(const Number a, const VectorSpaceVector<Number> &V) override;
 
+      /**
+       * Return the mean value of all the entries of this vector.
+       */
+      virtual value_type mean_value() const override;
+
       /**
        * Return the l<sub>1</sub> norm of the vector (i.e., the sum of the
        * absolute values of all entries among all processors).
diff --git a/include/deal.II/lac/la_vector.h b/include/deal.II/lac/la_vector.h
index 4d6b983039..144315e5c0 100644
--- a/include/deal.II/lac/la_vector.h
+++ b/include/deal.II/lac/la_vector.h
@@ -179,7 +179,7 @@ namespace LinearAlgebra
     virtual void equ(const Number a, const VectorSpaceVector<Number> &V);
 
     /**
-     * Return the mean value of the element of this vector.
+     * Return the mean value of all the entries of this vector.
      */
     virtual typename VectorSpaceVector<Number>::value_type mean_value() const;
 
diff --git a/include/deal.II/lac/vector_space_vector.h b/include/deal.II/lac/vector_space_vector.h
index 491104467f..cb8c95fdc8 100644
--- a/include/deal.II/lac/vector_space_vector.h
+++ b/include/deal.II/lac/vector_space_vector.h
@@ -129,7 +129,7 @@ namespace LinearAlgebra
     virtual void equ(const Number a, const VectorSpaceVector<Number> &V) = 0;
 
     /**
-     * Return the mean value of the element of this vector.
+     * Return the mean value of all the entries of this vector.
      */
     virtual value_type mean_value() const = 0;
 
diff --git a/source/lac/cuda_vector.cu b/source/lac/cuda_vector.cu
index 1e3b295edf..c0aa183792 100644
--- a/source/lac/cuda_vector.cu
+++ b/source/lac/cuda_vector.cu
@@ -93,6 +93,33 @@ namespace LinearAlgebra
       }
 
 
+
+      template <typename Number>
+      struct ElemSum
+      {
+        __device__ static Number reduction_op(const Number a, const Number b)
+        {
+          return (a + b);
+        }
+
+        __device__ static Number atomic_op(Number *dst, const Number a)
+        {
+          return atomicAdd_wrapper(dst, a);
+        }
+
+        __device__ static Number element_wise_op(const Number a)
+        {
+          return a;
+        }
+
+        __device__ static Number null_value()
+        {
+          return Number();
+        }
+      };
+
+
+
       template <typename Number>
       struct L1Norm
       {
@@ -118,6 +145,7 @@ namespace LinearAlgebra
       };
 
 
+
       template <typename Number>
       struct LInfty
       {
@@ -864,6 +892,34 @@ namespace LinearAlgebra
 
 
 
+    template <typename Number>
+    typename Vector<Number>::value_type Vector<Number>::mean_value() const
+    {
+      Number *result_device;
+      cudaError_t error_code = cudaMalloc(&result_device, sizeof(Number));
+      AssertCuda(error_code);
+      error_code = cudaMemset(result_device, Number(), sizeof(Number));
+
+      const int n_blocks = 1 + (n_elements-1)/(CHUNK_SIZE*BLOCK_SIZE);
+      internal::reduction<Number, internal::ElemSum<Number>>
+                                                          <<<dim3(n_blocks,1),dim3(BLOCK_SIZE)>>> (
+                                                            result_device, val,
+                                                            n_elements);
+
+      // Copy the result back to the host
+      Number result;
+      error_code = cudaMemcpy(&result, result_device, sizeof(Number),
+                              cudaMemcpyDeviceToHost);
+      AssertCuda(error_code);
+      // Free the memory on the device
+      error_code = cudaFree(result_device);
+      AssertCuda(error_code);
+
+      return result/static_cast<typename Vector<Number>::value_type>(n_elements);
+    }
+
+
+
     template <typename Number>
     typename Vector<Number>::real_type Vector<Number>::l1_norm() const
     {
diff --git a/tests/cuda/cuda_vector_02.cc b/tests/cuda/cuda_vector_02.cc
index c17999c603..09b317cfc0 100644
--- a/tests/cuda/cuda_vector_02.cc
+++ b/tests/cuda/cuda_vector_02.cc
@@ -84,6 +84,8 @@ void test()
     AssertThrow(2.*read_write_2[i]==read_write_3[i],
                 ExcMessage("Problem in equ."));
 
+  AssertThrow(b.mean_value()==54.50, ExcMessage("Problem in mean_value."));
+
   AssertThrow(b.l1_norm()==5450., ExcMessage("Problem in l1_norm."));
 
   const double eps=1e-3;