--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+#ifndef dealii__cuda_atomic_cuh
+#define dealii__cuda_atomic_cuh
+
+#include <deal.II/base/config.h>
+
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace LinearAlgebra
+{
+ namespace CUDAWrappers
+ {
+ /**
+ * Provide atomicAdd for floats.
+ */
+ inline __device__ float atomicAdd_wrapper(float *address, float val)
+ {
+ return atomicAdd(address,val);
+ }
+
+
+
+ /**
+ * Provide atomicAdd for doubles.
+ */
+ inline __device__ double atomicAdd_wrapper(double *address, double val)
+ {
+ // Use native instruction for CUDA 8 on Pascal or newer architecture
+#if __CUDACC_VER_MAJOR__ >= 8 && ( !defined(__CUDA_ARCH__) || __CUDA_ARCH__ >= 600 )
+ return atomicAdd(address,val);
+#else
+
+ unsigned long long int *address_as_ull =
+ reinterpret_cast<unsigned long long int *>(address);
+ unsigned long long int old = *address_as_ull, assumed;
+ do
+ {
+ assumed = old;
+ old = atomicCAS(address_as_ull, assumed,
+ __double_as_longlong(val + __longlong_as_double(assumed)));
+ }
+ while (assumed != old);
+
+ return __longlong_as_double(old);
+#endif
+ }
+
+
+
+ /**
+ * Provide atomicMax for floats.
+ */
+ inline __device__ float atomicMax_wrapper(float *address, float val)
+ {
+ int *address_as_int = reinterpret_cast<int *>(address);
+ int old = *address_as_int, assumed;
+ do
+ {
+ assumed = old;
+ old = atomicCAS(address_as_int, assumed,
+ atomicMax(address_as_int, __float_as_int(val)));
+ }
+ while (assumed != old);
+
+ return __longlong_as_double(old);
+ }
+
+
+
+ /**
+ * Provide atomicMax for doubles.
+ */
+ inline __device__ double atomicMax_wrapper(double *address, double val)
+ {
+ unsigned long long int *address_as_ull =
+ reinterpret_cast<unsigned long long int *>(address);
+ unsigned long long int old = *address_as_ull, assumed;
+ do
+ {
+ assumed = old;
+ old = atomicCAS(address_as_ull, assumed,
+ atomicMax(address_as_ull,
+ static_cast<unsigned long long int>(__double_as_longlong(val))));
+ }
+ while (assumed != old);
+
+ return __longlong_as_double(old);
+ }
+ }
+}
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
// ---------------------------------------------------------------------
#include <deal.II/lac/cuda_vector.h>
+#include <deal.II/lac/cuda_atomic.cuh>
#include <deal.II/lac/read_write_vector.h>
#include <deal.II/base/exceptions.h>
#include <cmath>
{
__device__ static Number reduction_op(const Number a, const Number b)
{
- return std::abs(a) + std::abs(b);
+ return (a + b);
}
- __device__ static void atomic_op(Number *dst, const Number a)
+ __device__ static Number atomic_op(Number *dst, const Number a)
{
- *dst = std::abs(*dst) + std::abs(a);
+ return atomicAdd_wrapper(dst, a);
+ }
+
+ __device__ static Number element_wise_op(const Number a)
+ {
+ return std::fabs(a);
}
__device__ static Number null_value()
{
__device__ static Number reduction_op(const Number a, const Number b)
{
- if (std::abs(a) > std::abs(b))
- return std::abs(a);
+ if (a > b)
+ return a;
else
- return std::abs(b);
+ return b;
}
- __device__ static void atomic_op(Number *dst, const Number a)
+ __device__ static Number atomic_op(Number *dst, const Number a)
{
- if (std::abs(*dst) < std::abs(a))
- *dst = std::abs(a);
- else
- *dst = std::abs(*dst);
+ return atomicMax_wrapper(dst, a);
+ }
+
+ __device__ static Number element_wise_op(const Number a)
+ {
+ return std::fabs(a);
}
__device__ static Number null_value()
const typename Vector<Number>::size_type local_idx = threadIdx.x;
if (global_idx<N)
- result_buffer[local_idx] = v[global_idx];
+ result_buffer[local_idx] = Operation::element_wise_op(v[global_idx]);
else
result_buffer[local_idx] = Operation::null_value();
return a+b;
}
- __device__ static void atomic_op(Number *dst, const Number a)
+ __device__ static Number atomic_op(Number *dst, const Number a)
{
- *dst += a;
+ return atomicAdd_wrapper(dst, a);
}
__device__ static Number null_value()