Response to the upstream review of PR #186

CMakeLists.txt

@@ -105,8 +105,9 @@ if(USE_HIP)
   ###################
   #  HIP (AMD)
   ###################
-  # rocThrust supplies the thrust::cuda::par compatibility namespace, so the
-  # Thrust code in s_filtergrid.cu compiles unchanged on ROCm.
+  # rocThrust is the AMD build of Thrust. Its stream-bound execution policy is
+  # in a different namespace than the NVIDIA one, which is what
+  # src/popsift/common/thrust_setup.h hides from the Thrust callers.
   find_package(rocthrust REQUIRED CONFIG)
 
   set(PopSift_CXX_STANDARD 17)
@@ -116,9 +117,9 @@ if(USE_HIP)
   set(CMAKE_HIP_STANDARD_REQUIRED ON)
   set(CMAKE_HIP_SEPARABLE_COMPILATION ON)
 
-  # HIP has __shfl_down etc. unconditionally; assist.h selects the _sync
-  # spelling on this flag and cuda_to_hip.h maps _sync to the mask-free builtin.
-  set(PopSift_HAVE_SHFL_DOWN_SYNC   1)
+  # HIP provides the mask-free __shfl_down etc., i.e. the spelling assist.h
+  # selects when this flag is off.
+  set(PopSift_HAVE_SHFL_DOWN_SYNC   0)
 else()
   ###################
   #  CUDA

src/CMakeLists.txt

@@ -9,6 +9,7 @@ add_library(popsift
         popsift/s_image.cu popsift/s_image.h
         popsift/sift_pyramid.cu popsift/sift_pyramid.h
         popsift/sift_octave.cu popsift/sift_octave.h
+        popsift/sift_textures.h
         popsift/s_pyramid_build.cu
         popsift/s_pyramid_build_aa.cu popsift/s_pyramid_build_aa.h
         popsift/s_pyramid_build_ai.cu popsift/s_pyramid_build_ai.h
@@ -39,6 +40,7 @@ add_library(popsift
         popsift/common/warp_bitonic_sort.h
         popsift/common/excl_blk_prefix_sum.h
         popsift/common/vec_macros.h
+        popsift/common/thrust_setup.h
         popsift/common/clamp.h)
 
 if(USE_HIP)

src/popsift/common/assist.h

@@ -8,7 +8,9 @@
 #pragma once
 
 #include <popsift/sift_config.h>
+#include <popsift/sift_textures.h>
 
+#include <algorithm>
 #include <cuda_runtime.h>
 #include <iostream>
 #include <thread>
@@ -42,6 +44,10 @@ template<typename T> __device__ inline T shuffle_up  ( T variable, int delta, in
 template<typename T> __device__ inline T shuffle_down( T variable, int delta, int ws ) { return __shfl_down_sync( 0xffffffff, variable, delta, ws ); }
 template<typename T> __device__ inline T shuffle_xor ( T variable, int delta, int ws ) { return __shfl_xor_sync ( 0xffffffff, variable, delta, ws ); }
 #else
+/* The mask-free builtins. This is the pre-CUDA-9 spelling, and it is also the
+ * spelling HIP provides, where the width parameter may be up to the 64-lane
+ * wavefront size.
+ */
 template<typename T> __device__ inline T shuffle     ( T variable, int src   ) { return __shfl     ( variable, src   ); }
 template<typename T> __device__ inline T shuffle_up  ( T variable, int delta ) { return __shfl_up  ( variable, delta ); }
 template<typename T> __device__ inline T shuffle_down( T variable, int delta ) { return __shfl_down( variable, delta ); }
@@ -102,12 +108,11 @@ float readTex( cudaTextureObject_t tex, float x, float y, float z )
      * close by) iff we add 0.5f to X and Y coordinate.
      */
 #if defined(USE_HIP) || defined(__HIP_PLATFORM_AMD__)
-    /* Observed on gfx90a, ROCm 7.2.1: hipCreateTextureObject rejects a
-     * hardware-linear-filtered texture over an element-read float array (see
-     * sift_octave.cu), so the linear
-     * textures are created with point filtering and we reproduce CUDA's bilinear
-     * filter in software here. This is an empirical limitation on this device/ROCm;
-     * it may not hold on other arches (re-verify on RDNA).
+    /* Observed on gfx90a, ROCm 7.2.1: creating a hardware-linear-filtered
+     * texture over an element-read float array is rejected (see sift_octave.cu;
+     * gfx1100 accepts it, so this is per device). The linear textures are
+     * therefore created with point filtering and the bilinear filter is done
+     * here in software, which keeps one build correct on either device.
      * CUDA's unnormalized linear filter on tex2DLayered(c) samples at index c-0.5,
      * i.e. i0=floor(c-0.5) with weight frac=(c-0.5)-i0; readTex passes c=x+0.5,
      * so i0=floor(x), frac=x-floor(x). Point textures (used with integer x,y) hit
@@ -123,70 +128,39 @@ float readTex( cudaTextureObject_t tex, float x, float y, float z )
     const float t10 = tex2DLayered<float>( tex, fx + 1.5f, fy + 0.5f, z );
     const float t01 = tex2DLayered<float>( tex, fx + 0.5f, fy + 1.5f, z );
     const float t11 = tex2DLayered<float>( tex, fx + 1.5f, fy + 1.5f, z );
-    const float top = t00 + ax * ( t10 - t00 );
-    const float bot = t01 + ax * ( t11 - t01 );
-    return top + ay * ( bot - top );
+    const float top = fmaf( ax, t10 - t00, t00 );
+    const float bot = fmaf( ax, t11 - t01, t01 );
+    return fmaf( ay, bot - top, top );
 #else
     return tex2DLayered<float>( tex, x+0.5f, y+0.5f, z );
 #endif
 }
 
 #if defined(USE_HIP) || defined(__HIP_PLATFORM_AMD__)
-/* HIP-only: bundle a pyramid array's surface with its per-level width/height.
- *
- * Observed on gfx90a/CDNA2 (ROCm 7.2.1): LAYERED images are broken. After a layered
- * array is written one layer at a time via surf2DLayeredwrite, a read in a later
- * kernel launch (tex2DLayered, surf2DLayeredread, and even host hipMemcpy3D) returns
- * one single (last-written) layer's data for EVERY layer index -- the layer
- * dimension is collapsed. Filed as ROCm/clr#275 (popsift is the motivating case);
- * a standalone reproducer mirrors that issue's array3d_check. AMD's partial
- * fix ROCm/rocm-systems#6683 corrects only surf2DLayered; tex2DLayered and
- * hipMemcpy3D may still collapse the layer dimension independently. A standalone
- * repro confirms this AND confirms a non-layered 3D array
- * (surf3Dwrite/surf3Dread/tex3D) is fully coherent across launches. So on
- * HIP the pyramid arrays are allocated as non-layered 3D arrays (sift_octave.cu)
- * and every layered access maps to the 3D form with the layer index as the z
- * coordinate. Reads funnel through readTex below and sample the surface via
- * surf3Dread; the surface address is a 1:1 match for the surf2DLayeredwrite the
- * producer kernels issue (now mapped to surf3Dwrite in cuda_to_hip.h), so the
- * blurred value a producer wrote at (x,y,level) is exactly what the consumer reads.
- * tex is carried for parity with the CUDA signature but is unused on HIP.
- * On CUDA this struct is unused (callers pass the texture object straight to the
- * tex2DLayered readTex above) and the read path is byte-for-byte unchanged.
- */
-struct LayeredTex
-{
-    cudaTextureObject_t tex;
-    cudaSurfaceObject_t surf;
-    int                 width;
-    int                 height;
-};
-
-/* Point-fetch one texel from the 3D pyramid surface, clamping the integer coords
- * to [0,width-1]x[0,height-1] to reproduce the cudaAddressModeClamp behaviour of
- * the point/linear textures these reads replace (an out-of-range surf3Dread
- * returns 0 on HIP, which would corrupt image borders). The surface x is a byte
- * offset (sizeof(float)), matching the surf2DLayeredwrite call sites (idx*4); z is
- * the (integer) blur level.
+/* Fetch one texel of the pyramid array through its surface, emulating a
+ * point-filtered texture fetch: the integer coordinates are clamped into
+ * [0,width-1]x[0,height-1] the way cudaAddressModeClamp clamps the textures this
+ * read replaces (an out-of-range surface read returns 0, which would corrupt the
+ * image borders), and the byte offset of the surface x coordinate is hidden. The
+ * z coordinate is the blur level, i.e. the layer index of the CUDA build.
  */
 __device__ static inline
-float surfFetchClamped( const LayeredTex& s, int ix, int iy, int layer )
+float texFetchClamped( const LayeredReadTex& s, int ix, int iy, int layer )
 {
-    ix = ix < 0 ? 0 : ( ix >= s.width  ? s.width  - 1 : ix );
-    iy = iy < 0 ? 0 : ( iy >= s.height ? s.height - 1 : iy );
+    ix = std::clamp( ix, 0, s.width  - 1 );
+    iy = std::clamp( iy, 0, s.height - 1 );
     float v;
     surf3Dread( &v, s.surf, ix * 4, iy, layer );
     return v;
 }
 
 __device__ static inline
-float readTex( const LayeredTex& s, float x, float y, float z )
+float readTex( const LayeredReadTex& s, float x, float y, float z )
 {
     /* Same -0.5 texel-center bilinear convention as the texture readTex above,
-     * but sampling the coherent 3D surface instead of the broken layered texture.
-     * SIFT blurs are per-level, so z is always an integer level: interpolate in
-     * x,y at the fixed z slice (no z interpolation), matching the CUDA layered
-     * point/linear texture fetch.
+     * but sampling the 3D surface (see sift_textures.h). SIFT blurs are
+     * per-level, so z is always an integer level: interpolate in x,y at the
+     * fixed z slice, matching the CUDA layered point/linear texture fetch.
      */
     const int   layer = (int)z;
     const float fx = floorf( x );
@@ -195,32 +169,16 @@ float readTex( const LayeredTex& s, float x, float y, float z )
     const float ay = y - fy;
     const int   ix = (int)fx;
     const int   iy = (int)fy;
-    const float t00 = surfFetchClamped( s, ix,   iy,   layer );
-    const float t10 = surfFetchClamped( s, ix+1, iy,   layer );
-    const float t01 = surfFetchClamped( s, ix,   iy+1, layer );
-    const float t11 = surfFetchClamped( s, ix+1, iy+1, layer );
-    const float top = t00 + ax * ( t10 - t00 );
-    const float bot = t01 + ax * ( t11 - t01 );
-    return top + ay * ( bot - top );
+    const float t00 = texFetchClamped( s, ix,   iy,   layer );
+    const float t10 = texFetchClamped( s, ix+1, iy,   layer );
+    const float t01 = texFetchClamped( s, ix,   iy+1, layer );
+    const float t11 = texFetchClamped( s, ix+1, iy+1, layer );
+    const float top = fmaf( ax, t10 - t00, t00 );
+    const float bot = fmaf( ax, t11 - t01, t01 );
+    return fmaf( ay, bot - top, top );
 }
 #endif
 
-/* Source handle for a pyramid-array layered read funnelled through readTex.
- * On HIP it is a LayeredTex (coherent 3D surface + dims); on CUDA it is the plain
- * texture object, so kernel signatures and the read path stay identical to upstream
- * on CUDA. Build it at the launch site with POPSIFT_LAYERED_SRC, which on CUDA
- * expands to exactly the texture argument (no behavioural change) and on HIP packs
- * the array's surface and per-level width/height for surf3Dread.
- */
-#if defined(USE_HIP) || defined(__HIP_PLATFORM_AMD__)
-using LayeredReadTex = LayeredTex;
-#define POPSIFT_LAYERED_SRC( tex, surf, width, height ) \
-    ::popsift::LayeredTex{ (tex), (surf), (width), (height) }
-#else
-using LayeredReadTex = cudaTextureObject_t;
-#define POPSIFT_LAYERED_SRC( tex, surf, width, height ) (tex)
-#endif
-
 __device__ static inline
 float readTex( cudaTextureObject_t tex, float x, float y )
 {

src/popsift/common/excl_blk_prefix_sum.h

@@ -91,16 +91,12 @@ class Block
             int self = (valid) ? _reader.get(cell) : 0;
 
             // This loop is an exclusive prefix sum for one warp. The block is
-            // (32,blockDim.y), one warp per threadIdx.y row. On a 64-lane
-            // wavefront two rows share a wavefront, so the scan shuffle must be
-            // confined to a width-32 sub-group (threadIdx.x is the in-row lane id)
-            // or odd rows pull partial sums from the wrong row. CUDA: unchanged.
+            // (32,blockDim.y), one row of 32 threads per threadIdx.y. The
+            // shuffle width is that row width, not the hardware warp size: on a
+            // 64-lane wavefront two rows share a wavefront and an unrestricted
+            // shuffle would pull partial sums from the wrong row.
             for( int s=0; s<5; s++ ) {
-#if defined(USE_HIP) || defined(__HIP_PLATFORM_AMD__)
                 const int add = popsift::shuffle_up( ews+self, 1<<s, 32 );
-#else
-                const int add = popsift::shuffle_up( ews+self, 1<<s );
-#endif
                 ews += threadIdx.x < (1<<s) ? 0 : add;
             }
 
@@ -117,11 +113,7 @@ class Block
                 int self = sum[threadIdx.x];
 
                 for( int s=0; s<5; s++ ) {
-#if defined(USE_HIP) || defined(__HIP_PLATFORM_AMD__)
                     const int add = popsift::shuffle_up( ebs+self, 1<<s, 32 );
-#else
-                    const int add = popsift::shuffle_up( ebs+self, 1<<s );
-#endif
                     ebs += threadIdx.x < (1<<s) ? 0 : add;
                 }
 

src/popsift/common/thrust_setup.h

@@ -0,0 +1,27 @@
+/*
+ * Copyright 2016-2017, Simula Research Laboratory
+ *
+ * This Source Code Form is subject to the terms of the Mozilla Public
+ * License, v. 2.0. If a copy of the MPL was not distributed with this
+ * file, You can obtain one at http://mozilla.org/MPL/2.0/.
+ */
+#pragma once
+
+/* Everything a translation unit needs before it can call a Thrust algorithm on
+ * one of PopSift's streams. Grid filtering is the only user today, but Thrust is
+ * usable elsewhere.
+ */
+
+#include <thrust/execution_policy.h>
+#include <thrust/version.h>
+
+/* The stream-bound parallel execution policy. It lives in thrust::cuda on
+ * NVIDIA and in thrust::hip in the AMD build of Thrust. Fully qualified from
+ * the global namespace, because inside namespace popsift the name cuda would
+ * otherwise resolve to popsift::cuda (common/debug_macros.h).
+ */
+#if defined(USE_HIP) || defined(__HIP_PLATFORM_AMD__)
+#define POPSIFT_THRUST_PAR ::thrust::hip::par
+#else
+#define POPSIFT_THRUST_PAR ::thrust::cuda::par
+#endif

src/popsift/common/warp_bitonic_sort.h

@@ -57,32 +57,20 @@ class Warp32
     __device__ inline
     int shiftit( const int my_index, const int shift, const int direction, const bool increasing )
     {
-        // This is a 32-lane bitonic step. On a 64-lane wavefront ori_par runs as
-        // a single 32-thread block (lanes 0-31), and sort64 holds 2 elements per
-        // lane; the shuffles must stay inside a 32-lane group, so force width 32
-        // on HIP (a no-op on CUDA's 32-lane warp) and key the swap direction off
-        // the in-group lane id.
-#if defined(USE_HIP) || defined(__HIP_PLATFORM_AMD__)
-        const int  lane_id     = threadIdx.x & 31;
+        // This is a 32-lane bitonic step: the caller launches a 32-thread block
+        // and sort64 holds 2 elements per thread. The shuffle width is that
+        // network width, not the hardware warp size, so a 64-lane wavefront does
+        // not exchange elements across the two halves of the sorting network.
         const T    my_val      = _array[my_index];
         const T    other_val   = popsift::shuffle_xor( my_val, 1 << shift, 32 );
-#else
-        const int  lane_id     = threadIdx.x;
-        const T    my_val      = _array[my_index];
-        const T    other_val   = popsift::shuffle_xor( my_val, 1 << shift );
-#endif
-        const bool reverse     = ( lane_id & ( 1 << direction ) );
-        const bool id_less     = ( ( lane_id & ( 1 << shift ) ) == 0 );
+        const bool reverse     = ( threadIdx.x & ( 1 << direction ) );
+        const bool id_less     = ( ( threadIdx.x & ( 1 << shift ) ) == 0 );
         const bool my_more     = id_less ? ( my_val > other_val )
                                          : ( my_val < other_val );
         const bool must_swap   = ! ( my_more ^ reverse ^ increasing );
 
         int lane = must_swap ? ( 1 << shift ) : 0;
-#if defined(USE_HIP) || defined(__HIP_PLATFORM_AMD__)
         return popsift::shuffle_xor( my_index, lane, 32 );
-#else
-        return popsift::shuffle_xor( my_index, lane );
-#endif
     }
 
     __device__ inline

src/popsift/cuda_to_hip.h

@@ -17,24 +17,6 @@
 #if defined(USE_HIP) || defined(__HIP_PLATFORM_AMD__)
 
 #include <hip/hip_runtime.h>
-// Backward-compat note (ROCm 7.2.x -> newer ROCm): newer ROCm's
-// <hip/hip_bf16.h> defines real __shfl_*_sync<...> functions whose names collide
-// with the function-like __shfl_*_sync macros defined further below -- if a
-// rocThrust/rocprim include pulls that header in AFTER the macros, the macros
-// rewrite the header's own definitions and it fails to compile ("use of
-// undeclared identifier 'mask'"). Pull the header in FIRST so the real functions
-// are defined before the macros. On older ROCm (7.2.x) the header exists and has
-// no such functions, so this is a harmless no-op there -- guarded with
-// __has_include so it also tolerates any ROCm that ships without the header.
-// Only needed where rocThrust/rocprim is pulled in (the device .cu TUs, compiled
-// by HIP clang); gate on __clang__ so the plain-C++ host example consumers (which
-// may be built by gcc and never include rocThrust) do not try to parse hip_bf16.h,
-// whose vector intrinsics require clang builtins.
-#if defined(__clang__) && defined(__has_include)
-#  if __has_include(<hip/hip_bf16.h>)
-#    include <hip/hip_bf16.h>
-#  endif
-#endif
 
 // ---- Error handling ----
 #define cudaError_t            hipError_t
@@ -150,12 +132,15 @@
 #define cudaBoundaryModeClamp     hipBoundaryModeClamp
 #define cudaBoundaryModeTrap      hipBoundaryModeTrap
 
-// HIP layered-image coherency is broken on gfx90a/CDNA2 (observed ROCm 7.2.1),
-// filed as ROCm/clr#275 (the partial fix ROCm/rocm-systems#6683 covers only
-// surf2DLayered): after a layered array is written layer-by-layer via
-// surf2DLayeredwrite, a read in a later kernel launch (tex2DLayered OR
-// surf2DLayeredread, host hipMemcpy3D too) returns a single (last-written) layer's
-// data for EVERY layer index -- the layer dimension is effectively collapsed.
+// HIP layered images are broken on gfx90a/CDNA2 and on gfx1100/RDNA3 (observed
+// ROCm 7.2.1), filed as ROCm/clr#275: after a layered array is written
+// layer-by-layer via surf2DLayeredwrite, a read in a later kernel launch
+// (tex2DLayered OR surf2DLayeredread, host hipMemcpy3D too) returns a single
+// (last-written) layer's data for EVERY layer index -- the layer dimension is
+// effectively collapsed. The defect is in the write: surf2DLayeredwrite passed
+// the layer index in the mipmap level slot, so every layer landed in the same
+// slot. ROCm/rocm-systems#6683 corrects that, and with it every read path above
+// returns correct per-layer data, but it is not in ROCm 7.2.x.
 // A standalone reproducer confirms it and also
 // confirms a NON-layered 3D array (surf3Dwrite/surf3Dread/tex3D) is fully
 // coherent across launches. So on HIP the pyramid arrays are allocated as
@@ -186,19 +171,8 @@ __device__ __forceinline__ void popsift_surf2DLayeredwrite(
 #define __fmul_ru(a, b)     __fmul_rn((a), (b))
 
 // ---- Warp intrinsics ----
-// PopSift (assist.h) passes the CUDA 32-bit full mask 0xffffffff to the *_sync
-// builtins. On a 64-lane CDNA wavefront a uint32 mask is meaningless; HIP's
-// mask-free builtins poll the whole active wavefront, which is the faithful
-// equivalent of "operate over the (sub)warp" on the blocks PopSift launches.
-// Map the _sync forms (both the 2-arg and the explicit-width 3/4-arg overloads
-// resolve to these) to the mask-free HIP builtins.
-#define __shfl_sync(mask, ...)        __shfl(__VA_ARGS__)
-#define __shfl_up_sync(mask, ...)     __shfl_up(__VA_ARGS__)
-#define __shfl_down_sync(mask, ...)   __shfl_down(__VA_ARGS__)
-#define __shfl_xor_sync(mask, ...)    __shfl_xor(__VA_ARGS__)
-#define __ballot_sync(mask, pred)     __ballot(pred)
-#define __any_sync(mask, pred)        __any(pred)
-#define __all_sync(mask, pred)        __all(pred)
+// Nothing to map: HIP provides the mask-free __shfl/__ballot/__any/__all
+// builtins, so assist.h picks them through PopSift_HAVE_SHFL_DOWN_SYNC=0.
 
 #else // NVIDIA / CUDA