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 8a925f6346582465d93f9bc8916497fe0a7a6478..c91751d4867f86159d1402d4747de44ad72ff0a5 100644 (file)
--- 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 4d6b9830395ff341f4aeab64ce749a6b1485714d..144315e5c0744af61edd396d6fb90334a80f95ec 100644 (file)
--- 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 491104467f46c00f7cac9d0c6e0be50ca659c996..cb8c95fdc8afda8d4e78370160099b6191910d6d 100644 (file)
--- 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 1e3b295edfc0b8f7afaee770bdf02f1bf8b6701f..c0aa183792c0c7f67c53008874ed3b921bb5ad1d 100644 (file)
--- 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 c17999c603770475345195315fd68dfa4024af25..09b317cfc075a2d3bfbc092b1b6037083c5f2aed 100644 (file)
--- 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;