Tensor: Use native atomic addition for double reductions

benchmarks/GPU/bench_gpu_reduction.cu

diff --git a/benchmarks/GPU/bench_gpu_reduction.cu b/benchmarks/GPU/bench_gpu_reduction.cu
index a59535c..2ea8b26 100644
--- a/benchmarks/GPU/bench_gpu_reduction.cu
+++ b/benchmarks/GPU/bench_gpu_reduction.cu
@@ -223,8 +223,8 @@
 
 // Explicit shapes rather than a product: the product of the two ranges asks for tensors of up to 2^40 elements.
 // Args is (rows, columns), so for the outer reduction it reads as (outputs, reduced extent) and the last three
-// shapes straddle the band OuterReducer::run admits doubles to OuterReductionKernel in: a reduced extent below
-// its floor, a shape inside the band, and an output count above its ceiling on a 36-multiprocessor device. For
+// shapes straddle the former CUDA double-sum band: a reduced extent below its floor, a shape inside the band,
+// and an output count above its ceiling on a 36-multiprocessor device. HIP retains that band. For
 // the inner reduction and the CUB baseline the two extents swap roles.
 #define EIGEN_GPU_PARTIAL_REDUCTION_SHAPES(NAME, TYPE) \
   BENCHMARK_TEMPLATE(NAME, TYPE)                       \

contrib/Eigen/src/Tensor/TensorReductionGpu.h

diff --git a/contrib/Eigen/src/Tensor/TensorReductionGpu.h b/contrib/Eigen/src/Tensor/TensorReductionGpu.h
index 1dccb6c..ef17207 100644
--- a/contrib/Eigen/src/Tensor/TensorReductionGpu.h
+++ b/contrib/Eigen/src/Tensor/TensorReductionGpu.h
@@ -156,6 +156,11 @@
   atomicAdd(output, accum);
 }
 
+template <>
+__device__ inline void atomicReduce(double* output, double accum, SumReducer<double>&) {
+  atomicAdd(output, accum);
+}
+
 template <typename CoeffType, typename Index>
 __global__ EIGEN_HIP_LAUNCH_BOUNDS_1024 void ReductionInitKernel(const CoeffType val, Index num_preserved_coeffs,
                                                                  CoeffType* output) {
@@ -889,9 +894,14 @@
       return true;
     }
 
-    // Double reductions use compare-and-swap atomics. Restrict their parallel reduction to the shape band
-    // measured to amortize that contention on sm_89; outside it, use one thread per output.
-    if (std::is_same<OutputType, double>::value) {
+    // CUDA has native FP64 atomicAdd on every supported architecture (sm_60+). Keep the measured shape band
+    // for the other double reducers and for HIP, where atomicAdd may still use compare-and-swap.
+#if defined(EIGEN_CUDACC)
+    constexpr bool has_native_double_sum = std::is_same<Op, SumReducer<double>>::value;
+#else
+    constexpr bool has_native_double_sum = false;
+#endif
+    if (std::is_same<OutputType, double>::value && !has_native_double_sum) {
       const Index multi_processors = device.getNumGpuMultiProcessors();
       if (num_coeffs_to_reduce < 64 || num_preserved_vals < 8 * multi_processors ||
           num_preserved_vals > 64 * multi_processors) {

contrib/test/tensor_reduction_gpu.cu

diff --git a/contrib/test/tensor_reduction_gpu.cu b/contrib/test/tensor_reduction_gpu.cu
index 0e91784..e4f9ddf 100644
--- a/contrib/test/tensor_reduction_gpu.cu
+++ b/contrib/test/tensor_reduction_gpu.cu
@@ -53,11 +53,7 @@
 }
 
 template <typename Type, int DataLayout>
-static void test_first_dim_reductions() {
-  int dim_x = 33;
-  int dim_y = 1;
-  int dim_z = 128;
-
+static void test_first_dim_reductions(int dim_x = 33, int dim_y = 1, int dim_z = 128) {
   Tensor<Type, 3, DataLayout> in(dim_x, dim_y, dim_z);
   in.setRandom();
 
@@ -144,13 +140,16 @@
   CALL_SUBTEST_6((test_last_dim_reductions<float, RowMajor>()));
   CALL_SUBTEST_6((test_last_dim_reductions<double, RowMajor>()));
 
-  // A double outer reduction reaches OuterReductionKernel only inside the band of output counts OuterReducer::run
-  // gates on, and that band scales with the device's multiprocessor count; every shape above lies outside it and
-  // takes the generic path. The float cases are not gated and reach the kernel already.
+  // CUDA double sums also reach OuterReductionKernel below and above the former shape band. HIP keeps that band.
   {
     Eigen::GpuStreamDevice stream;
     Eigen::GpuDevice device(&stream);
-    const int outputs_in_band = 16 * device.getNumGpuMultiProcessors();
-    CALL_SUBTEST_7((test_last_dim_reductions<double, ColMajor>(outputs_in_band, 1, 128)));
+    const int multi_processors = device.getNumGpuMultiProcessors();
+    for (int outputs : {4 * multi_processors, 16 * multi_processors, 65 * multi_processors}) {
+      CALL_SUBTEST_7((test_last_dim_reductions<double, ColMajor>(outputs, 1, 128)));
+      CALL_SUBTEST_7((test_first_dim_reductions<double, RowMajor>(128, 1, outputs)));
+    }
+    CALL_SUBTEST_7((test_last_dim_reductions<double, ColMajor>(16 * multi_processors, 1, 33)));
+    CALL_SUBTEST_7((test_first_dim_reductions<double, RowMajor>(33, 1, 16 * multi_processors)));
   }
 }