diff --git a/CMakeLists.txt b/CMakeLists.txt index 07bda8ac..be199682 100755 --- a/CMakeLists.txt +++ b/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 diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 5bf4620f..2aaec2a6 100644 --- a/src/CMakeLists.txt +++ b/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) diff --git a/src/popsift/common/assist.h b/src/popsift/common/assist.h index 1c52c93d..23017f02 100644 --- a/src/popsift/common/assist.h +++ b/src/popsift/common/assist.h @@ -8,7 +8,9 @@ #pragma once #include +#include +#include #include #include #include @@ -42,6 +44,10 @@ template __device__ inline T shuffle_up ( T variable, int delta, in template __device__ inline T shuffle_down( T variable, int delta, int ws ) { return __shfl_down_sync( 0xffffffff, variable, delta, ws ); } template __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 __device__ inline T shuffle ( T variable, int src ) { return __shfl ( variable, src ); } template __device__ inline T shuffle_up ( T variable, int delta ) { return __shfl_up ( variable, delta ); } template __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( tex, fx + 1.5f, fy + 0.5f, z ); const float t01 = tex2DLayered( tex, fx + 0.5f, fy + 1.5f, z ); const float t11 = tex2DLayered( 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( 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 ) { diff --git a/src/popsift/common/excl_blk_prefix_sum.h b/src/popsift/common/excl_blk_prefix_sum.h index 7e37b76d..1bb24c50 100644 --- a/src/popsift/common/excl_blk_prefix_sum.h +++ b/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< +#include + +/* 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 diff --git a/src/popsift/common/warp_bitonic_sort.h b/src/popsift/common/warp_bitonic_sort.h index 23157a6b..6b13c9b7 100644 --- a/src/popsift/common/warp_bitonic_sort.h +++ b/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 diff --git a/src/popsift/cuda_to_hip.h b/src/popsift/cuda_to_hip.h index b4a24894..cc42b64a 100644 --- a/src/popsift/cuda_to_hip.h +++ b/src/popsift/cuda_to_hip.h @@ -17,24 +17,6 @@ #if defined(USE_HIP) || defined(__HIP_PLATFORM_AMD__) #include -// Backward-compat note (ROCm 7.2.x -> newer ROCm): newer ROCm's -// 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() -# include -# 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 diff --git a/src/popsift/features.cu b/src/popsift/features.cu index 8f6db5e2..ffbc47e0 100755 --- a/src/popsift/features.cu +++ b/src/popsift/features.cu @@ -175,22 +175,14 @@ l2_in_t0( const float4* lptr, const float4* rptr ) + mval.y * mval.y + mval.z * mval.z + mval.w * mval.w; - // 32-lane reduction over threadIdx.x (compute_distance is a 32-thread block). - // Force a width-32 sub-group on a 64-lane wavefront so the inactive upper - // half is never sampled. CUDA: width 32 is the whole warp, unchanged. -#if defined(USE_HIP) || defined(__HIP_PLATFORM_AMD__) + // Reduction over the 32 threads of compute_distance's block. The shuffle + // width is that block width, not the hardware warp size, so a 64-lane + // wavefront does not sample its inactive upper half. res += shuffle_down( res, 16, 32 ); res += shuffle_down( res, 8, 32 ); res += shuffle_down( res, 4, 32 ); res += shuffle_down( res, 2, 32 ); res += shuffle_down( res, 1, 32 ); -#else - res += shuffle_down( res, 16 ); - res += shuffle_down( res, 8 ); - res += shuffle_down( res, 4 ); - res += shuffle_down( res, 2 ); - res += shuffle_down( res, 1 ); -#endif return res; } diff --git a/src/popsift/s_desc_grid.h b/src/popsift/s_desc_grid.h index 63e3a991..bef5dcdd 100644 --- a/src/popsift/s_desc_grid.h +++ b/src/popsift/s_desc_grid.h @@ -40,10 +40,7 @@ inline static bool start_ext_desc_grid( const int octave, Octave& oct_obj ) ext_desc_grid <<>> ( octave, - POPSIFT_LAYERED_SRC( oct_obj.getDataTexPoint( ), - oct_obj.getDataSurface( ), - oct_obj.getWidth(), - oct_obj.getHeight() ) ); + oct_obj.getDataReadTexPoint( ) ); POP_SYNC_CHK; diff --git a/src/popsift/s_desc_igrid.h b/src/popsift/s_desc_igrid.h index 7bed27df..f8d4064f 100644 --- a/src/popsift/s_desc_igrid.h +++ b/src/popsift/s_desc_igrid.h @@ -41,10 +41,7 @@ inline static bool start_ext_desc_igrid( const int octave, Octave& oct_obj ) ext_desc_igrid <<>> ( octave, - POPSIFT_LAYERED_SRC( oct_obj.getDataTexLinear( ).tex, - oct_obj.getDataSurface( ), - oct_obj.getWidth(), - oct_obj.getHeight() ) ); + oct_obj.getDataReadTexLinear( ) ); POP_SYNC_CHK; diff --git a/src/popsift/s_desc_iloop.cu b/src/popsift/s_desc_iloop.cu index 0d8ed08c..4913d1ac 100644 --- a/src/popsift/s_desc_iloop.cu +++ b/src/popsift/s_desc_iloop.cu @@ -114,9 +114,8 @@ void ext_desc_iloop_sub( const float ang, dpt[0] += dpt[8]; /* reduction here */ - // 32-lane reduction over threadIdx.x; confine the shuffles to a width-32 - // sub-group on a 64-lane wavefront (see s_desc_loop.cu). CUDA unchanged. -#if defined(USE_HIP) || defined(__HIP_PLATFORM_AMD__) + // Reduction over the 32 threads of one descriptor tile; the shuffle width is + // that tile width, not the hardware warp size (see s_desc_loop.cu). for (int i = 0; i < 8; i++) { dpt[i] += popsift::shuffle_down( dpt[i], 16, 32 ); dpt[i] += popsift::shuffle_down( dpt[i], 8, 32 ); @@ -125,16 +124,6 @@ void ext_desc_iloop_sub( const float ang, dpt[i] += popsift::shuffle_down( dpt[i], 1, 32 ); dpt[i] = popsift::shuffle ( dpt[i], 0, 32 ); } -#else - for (int i = 0; i < 8; i++) { - dpt[i] += popsift::shuffle_down( dpt[i], 16 ); - dpt[i] += popsift::shuffle_down( dpt[i], 8 ); - dpt[i] += popsift::shuffle_down( dpt[i], 4 ); - dpt[i] += popsift::shuffle_down( dpt[i], 2 ); - dpt[i] += popsift::shuffle_down( dpt[i], 1 ); - dpt[i] = popsift::shuffle ( dpt[i], 0 ); - } -#endif if( threadIdx.x < 8 ) { features[tile+threadIdx.x] = dpt[threadIdx.x]; diff --git a/src/popsift/s_desc_iloop.h b/src/popsift/s_desc_iloop.h index 03a020ab..5af35c9e 100644 --- a/src/popsift/s_desc_iloop.h +++ b/src/popsift/s_desc_iloop.h @@ -35,10 +35,7 @@ inline static bool start_ext_desc_iloop( const int octave, Octave& oct_obj ) ext_desc_iloop <<>> ( octave, - POPSIFT_LAYERED_SRC( oct_obj.getDataTexLinear( ).tex, - oct_obj.getDataSurface( ), - oct_obj.getWidth(), - oct_obj.getHeight() ), + oct_obj.getDataReadTexLinear( ), oct_obj.getWidth(), oct_obj.getHeight() ); diff --git a/src/popsift/s_desc_loop.cu b/src/popsift/s_desc_loop.cu index d58a46f9..dab6352c 100644 --- a/src/popsift/s_desc_loop.cu +++ b/src/popsift/s_desc_loop.cu @@ -124,12 +124,10 @@ void ext_desc_loop_sub( const float ang, dpt[0] += dpt[8]; /* reduction here */ - // 32-lane reduction over threadIdx.x. The block is (32,4,4): each (y,z) pair - // is its own 32-thread group, so on a 64-lane wavefront (two groups per - // wavefront) the shuffles must be confined to a width-32 sub-group, else the - // reduction and the lane-0 broadcast leak across the group boundary and - // corrupt half the descriptors (NaN). On CUDA width 32 is the whole warp. -#if defined(USE_HIP) || defined(__HIP_PLATFORM_AMD__) + // The block is (32,4,4) and each (y,z) pair reduces its own 32 threads. The + // shuffle width is that group width, not the hardware warp size: on a 64-lane + // wavefront two groups share a wavefront, and an unrestricted reduction and + // lane-0 broadcast would leak across the group boundary. for (int i = 0; i < 8; i++) { dpt[i] += popsift::shuffle_down( dpt[i], 16, 32 ); dpt[i] += popsift::shuffle_down( dpt[i], 8, 32 ); @@ -138,16 +136,6 @@ void ext_desc_loop_sub( const float ang, dpt[i] += popsift::shuffle_down( dpt[i], 1, 32 ); dpt[i] = popsift::shuffle ( dpt[i], 0, 32 ); } -#else - for (int i = 0; i < 8; i++) { - dpt[i] += popsift::shuffle_down( dpt[i], 16 ); - dpt[i] += popsift::shuffle_down( dpt[i], 8 ); - dpt[i] += popsift::shuffle_down( dpt[i], 4 ); - dpt[i] += popsift::shuffle_down( dpt[i], 2 ); - dpt[i] += popsift::shuffle_down( dpt[i], 1 ); - dpt[i] = popsift::shuffle ( dpt[i], 0 ); - } -#endif if( threadIdx.x < 8 ) { features[tile+threadIdx.x] = dpt[threadIdx.x]; diff --git a/src/popsift/s_desc_loop.h b/src/popsift/s_desc_loop.h index b7e6084b..f9bf3764 100644 --- a/src/popsift/s_desc_loop.h +++ b/src/popsift/s_desc_loop.h @@ -43,10 +43,7 @@ inline static bool start_ext_desc_loop( const int octave, Octave& oct_obj ) ext_desc_loop <<>> ( octave, - POPSIFT_LAYERED_SRC( oct_obj.getDataTexPoint( ), - oct_obj.getDataSurface( ), - oct_obj.getWidth(), - oct_obj.getHeight() ), + oct_obj.getDataReadTexPoint( ), oct_obj.getWidth(), oct_obj.getHeight() ); diff --git a/src/popsift/s_desc_norm_l2.h b/src/popsift/s_desc_norm_l2.h index d1f42b9f..a8b34f43 100644 --- a/src/popsift/s_desc_norm_l2.h +++ b/src/popsift/s_desc_norm_l2.h @@ -58,41 +58,25 @@ void NormalizeL2::normalize( const float* src_desc, float* dst_desc, const bool + descr.y * descr.y + descr.z * descr.z + descr.w * descr.w; - // 32-lane reduction over threadIdx.x. The normalize block is (32,32): each - // threadIdx.y row is one descriptor. On a 64-lane wavefront two rows share a - // wavefront, so the reduction and the lane-0 broadcast must stay inside a - // width-32 sub-group, else rows cross-contaminate and the descriptor norm is - // wrong (NaN). CUDA: width 32 is the whole warp, unchanged. -#if defined(USE_HIP) || defined(__HIP_PLATFORM_AMD__) + // The normalize block is (32,32) and each threadIdx.y row holds one + // descriptor. 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 reduction and lane-0 broadcast would mix the two descriptors. norm += popsift::shuffle_down( norm, 16, 32 ); norm += popsift::shuffle_down( norm, 8, 32 ); norm += popsift::shuffle_down( norm, 4, 32 ); norm += popsift::shuffle_down( norm, 2, 32 ); norm += popsift::shuffle_down( norm, 1, 32 ); -#else - norm += popsift::shuffle_down( norm, 16 ); - norm += popsift::shuffle_down( norm, 8 ); - norm += popsift::shuffle_down( norm, 4 ); - norm += popsift::shuffle_down( norm, 2 ); - norm += popsift::shuffle_down( norm, 1 ); -#endif if( threadIdx.x == 0 ) { - // compute 1 / sqrt(sum) in round-to-nearest even mode in thread 0 - norm = __frsqrt_rn( norm ); -#if defined(USE_HIP) || defined(__HIP_PLATFORM_AMD__) - // A fully flat gradient window yields an all-zero descriptor; __frsqrt_rn(0) - // is +inf, and 0*inf below is NaN. Keep such degenerate descriptors all-zero. - if( ! isfinite( norm ) ) norm = 0.0f; -#endif + // compute 1 / sqrt(sum) in round-to-nearest even mode in thread 0. + // A sum of squares cannot be negative, and it is only zero for an + // all-zero descriptor, which stays all-zero instead of scaling by inf. + norm = ( norm > 0.0f ) ? __frsqrt_rn( norm ) : 0.0f; } // spread the inverted norm from thread 0 to all threads in the warp -#if defined(USE_HIP) || defined(__HIP_PLATFORM_AMD__) norm = popsift::shuffle( norm, 0, 32 ); -#else - norm = popsift::shuffle( norm, 0 ); -#endif // quasi-normalize all 128 floats descr.x = min( descr.x*norm, 0.2f ); @@ -108,33 +92,18 @@ void NormalizeL2::normalize( const float* src_desc, float* dst_desc, const bool + descr.y * descr.y + descr.z * descr.z + descr.w * descr.w; -#if defined(USE_HIP) || defined(__HIP_PLATFORM_AMD__) norm += popsift::shuffle_down( norm, 16, 32 ); norm += popsift::shuffle_down( norm, 8, 32 ); norm += popsift::shuffle_down( norm, 4, 32 ); norm += popsift::shuffle_down( norm, 2, 32 ); norm += popsift::shuffle_down( norm, 1, 32 ); -#else - norm += popsift::shuffle_down( norm, 16 ); - norm += popsift::shuffle_down( norm, 8 ); - norm += popsift::shuffle_down( norm, 4 ); - norm += popsift::shuffle_down( norm, 2 ); - norm += popsift::shuffle_down( norm, 1 ); -#endif if( threadIdx.x == 0 ) { - norm = __frsqrt_rn( norm ); // inverse square root -#if defined(USE_HIP) || defined(__HIP_PLATFORM_AMD__) - if( ! isfinite( norm ) ) norm = 0.0f; // see above: avoid 0*inf NaN -#endif + norm = ( norm > 0.0f ) ? __frsqrt_rn( norm ) : 0.0f; // inverse square root norm = scalbnf( norm, d_consts.norm_multi ); } -#if defined(USE_HIP) || defined(__HIP_PLATFORM_AMD__) norm = popsift::shuffle( norm, 0, 32 ); -#else - norm = popsift::shuffle( norm, 0 ); -#endif descr.x = descr.x * norm; descr.y = descr.y * norm; diff --git a/src/popsift/s_desc_norm_rs.h b/src/popsift/s_desc_norm_rs.h index fd3f7100..df8f6e22 100644 --- a/src/popsift/s_desc_norm_rs.h +++ b/src/popsift/s_desc_norm_rs.h @@ -48,10 +48,10 @@ void NormalizeRootSift::normalize( const float* src_desc, float* dst_desc, bool float sum = descr.x + descr.y + descr.z + descr.w; - // 32-lane reduction over threadIdx.x; confine to a width-32 sub-group on a - // 64-lane wavefront (normalize block is (32,32), one descriptor per row). - // CUDA unchanged. -#if defined(USE_HIP) || defined(__HIP_PLATFORM_AMD__) + // The normalize block is (32,32) and each threadIdx.y row holds one + // descriptor. 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 reduction and lane-0 broadcast would mix the two descriptors. sum += popsift::shuffle_down( sum, 16, 32 ); sum += popsift::shuffle_down( sum, 8, 32 ); sum += popsift::shuffle_down( sum, 4, 32 ); @@ -59,50 +59,28 @@ void NormalizeRootSift::normalize( const float* src_desc, float* dst_desc, bool sum += popsift::shuffle_down( sum, 1, 32 ); sum = popsift::shuffle( sum, 0, 32 ); -#else - sum += popsift::shuffle_down( sum, 16 ); - sum += popsift::shuffle_down( sum, 8 ); - sum += popsift::shuffle_down( sum, 4 ); - sum += popsift::shuffle_down( sum, 2 ); - sum += popsift::shuffle_down( sum, 1 ); - sum = popsift::shuffle( sum, 0 ); -#endif - -#if defined(USE_HIP) || defined(__HIP_PLATFORM_AMD__) - // RootSift takes sqrt(bin/sum). The fmaxf(.,0) below clamps a bin that came out - // slightly negative from the descriptor accumulation's round-toward-+inf - // intrinsics (mapped to round-to-nearest in cuda_to_hip.h; a negative bin is - // unphysical). The divisor is gated at a small threshold so a degenerate - // near-zero sum is treated as an all-zero descriptor: an all-flat window - // (sum==0) normalizes to 0, and a tiny subnormal sum can no longer make 1/sum - // overflow to +inf (which would normalize a positive bin to +inf). CUDA, where - // the directed-rounding intrinsics exist and bins stay non-negative, is unchanged. + // RootSift takes sqrt(bin/sum). A descriptor bin cannot be negative, so a + // sum that is not clearly positive means a degenerate descriptor, which + // stays all-zero instead of dividing. The test is against a small threshold + // rather than against zero because a subnormal sum makes the reciprocal + // overflow to infinity, which would normalize a positive bin to infinity. + // The per-bin test below keeps a bin that came out slightly negative + // (round-to-nearest weight accumulation on platforms without the + // round-toward-+inf intrinsics) out of the square root. const float inv = ( sum > 1e-20f ) ? __fdividef( 1.0f, sum ) : 0.0f; - float val; - val = scalbnf( __fsqrt_rn( fmaxf( descr.x * inv, 0.0f ) ), d_consts.norm_multi ); - descr.x = val; - val = scalbnf( __fsqrt_rn( fmaxf( descr.y * inv, 0.0f ) ), d_consts.norm_multi ); - descr.y = val; - val = scalbnf( __fsqrt_rn( fmaxf( descr.z * inv, 0.0f ) ), d_consts.norm_multi ); - descr.z = val; - val = scalbnf( __fsqrt_rn( fmaxf( descr.w * inv, 0.0f ) ), d_consts.norm_multi ); - descr.w = val; -#else - float val; - val = scalbnf( __fsqrt_rn( __fdividef( descr.x, sum ) ), - d_consts.norm_multi ); - descr.x = val; - val = scalbnf( __fsqrt_rn( __fdividef( descr.y, sum ) ), - d_consts.norm_multi ); - descr.y = val; - val = scalbnf( __fsqrt_rn( __fdividef( descr.z, sum ) ), - d_consts.norm_multi ); - descr.z = val; - val = scalbnf( __fsqrt_rn( __fdividef( descr.w, sum ) ), - d_consts.norm_multi ); - descr.w = val; -#endif + + if( inv <= 0.0f ) + { + descr.x = descr.y = descr.z = descr.w = 0.0f; + } + else + { + descr.x = descr.x <= 0.0f ? 0.0f : scalbnf( __fsqrt_rn( descr.x * inv ), d_consts.norm_multi ); + descr.y = descr.y <= 0.0f ? 0.0f : scalbnf( __fsqrt_rn( descr.y * inv ), d_consts.norm_multi ); + descr.z = descr.z <= 0.0f ? 0.0f : scalbnf( __fsqrt_rn( descr.z * inv ), d_consts.norm_multi ); + descr.w = descr.w <= 0.0f ? 0.0f : scalbnf( __fsqrt_rn( descr.w * inv ), d_consts.norm_multi ); + } if( ! ignoreme ) { float4* out4 = (float4*)dst_desc; diff --git a/src/popsift/s_desc_notile.cu b/src/popsift/s_desc_notile.cu index f6fe7789..2d3e2ed2 100644 --- a/src/popsift/s_desc_notile.cu +++ b/src/popsift/s_desc_notile.cu @@ -149,10 +149,7 @@ bool start_ext_desc_notile( int octave, Octave& oct_obj ) ext_desc_notile <<>> ( octave, - POPSIFT_LAYERED_SRC( oct_obj.getDataTexLinear( ).tex, - oct_obj.getDataSurface( ), - oct_obj.getWidth(), - oct_obj.getHeight() ) ); + oct_obj.getDataReadTexLinear( ) ); cudaDeviceSynchronize(); cudaError_t err = cudaGetLastError( ); POP_CUDA_FATAL_TEST(err, "cudaGetLastError failed: "); diff --git a/src/popsift/s_desc_vlfeat.cu b/src/popsift/s_desc_vlfeat.cu index 6509d875..a3b46914 100644 --- a/src/popsift/s_desc_vlfeat.cu +++ b/src/popsift/s_desc_vlfeat.cu @@ -201,10 +201,7 @@ bool start_ext_desc_vlfeat( const int octave, Octave& oct_obj ) ext_desc_vlfeat <<>> ( octave, - POPSIFT_LAYERED_SRC( oct_obj.getDataTexPoint( ), - oct_obj.getDataSurface( ), - oct_obj.getWidth(), - oct_obj.getHeight() ), + oct_obj.getDataReadTexPoint( ), oct_obj.getWidth(), oct_obj.getHeight() ); diff --git a/src/popsift/s_extrema.cu b/src/popsift/s_extrema.cu index 43ef2bbb..c276c58c 100644 --- a/src/popsift/s_extrema.cu +++ b/src/popsift/s_extrema.cu @@ -605,10 +605,7 @@ void Pyramid::find_extrema( const Config& conf ) case Config::VLFeat : find_extrema_in_dog <<>> - ( POPSIFT_LAYERED_SRC( oct_obj.getDogTexturePoint( ), - oct_obj.getDogSurface( ), - oct_obj.getWidth(), - oct_obj.getHeight() ), + ( oct_obj.getDogReadTexPoint( ), octave, cols, rows, @@ -623,10 +620,7 @@ void Pyramid::find_extrema( const Config& conf ) case Config::OpenCV : find_extrema_in_dog <<>> - ( POPSIFT_LAYERED_SRC( oct_obj.getDogTexturePoint( ), - oct_obj.getDogSurface( ), - oct_obj.getWidth(), - oct_obj.getHeight() ), + ( oct_obj.getDogReadTexPoint( ), octave, cols, rows, @@ -641,10 +635,7 @@ void Pyramid::find_extrema( const Config& conf ) default : find_extrema_in_dog <<>> - ( POPSIFT_LAYERED_SRC( oct_obj.getDogTexturePoint( ), - oct_obj.getDogSurface( ), - oct_obj.getWidth(), - oct_obj.getHeight() ), + ( oct_obj.getDogReadTexPoint( ), octave, cols, rows, diff --git a/src/popsift/s_filtergrid.cu b/src/popsift/s_filtergrid.cu index 00d1b9eb..cd726f77 100644 --- a/src/popsift/s_filtergrid.cu +++ b/src/popsift/s_filtergrid.cu @@ -11,10 +11,11 @@ #if ! POPSIFT_IS_DEFINED(POPSIFT_DISABLE_GRID_FILTER) +#include "common/thrust_setup.h" + #include #include #include -#include #include #include #include @@ -22,16 +23,6 @@ #include #include #include -#include - -// The stream-bound parallel execution policy lives in thrust::cuda on NVIDIA -// and thrust::hip on rocThrust. Fully qualified from the global namespace to -// avoid colliding with popsift::cuda (debug_macros.h) inside namespace popsift. -#if defined(USE_HIP) || defined(__HIP_PLATFORM_AMD__) -#define POPSIFT_THRUST_PAR ::thrust::hip::par -#else -#define POPSIFT_THRUST_PAR ::thrust::cuda::par -#endif namespace popsift { diff --git a/src/popsift/s_orientation.cu b/src/popsift/s_orientation.cu index 30b9b53f..1daedd0b 100644 --- a/src/popsift/s_orientation.cu +++ b/src/popsift/s_orientation.cu @@ -227,7 +227,7 @@ void ori_par( const int octave, // All threads retrieve the yval of thread 0, the largest // of all yvals. const float best_val = yval[best_index.x]; - const float yval_ref = 0.8f * popsift::shuffle( best_val, 0 ); + const float yval_ref = 0.8f * popsift::shuffle( best_val, 0, 32 ); const bool valid = ( best_val >= yval_ref ); bool written = false; @@ -411,10 +411,7 @@ void Pyramid::orientation( const Config& conf ) <<>> ( octave, hct.ext_ps[octave], - POPSIFT_LAYERED_SRC( oct_obj.getDataTexPoint( ), - oct_obj.getDataSurface( ), - oct_obj.getWidth( ), - oct_obj.getHeight( ) ), + oct_obj.getDataReadTexPoint( ), oct_obj.getWidth( ), oct_obj.getHeight( ) ); POP_SYNC_CHK; diff --git a/src/popsift/s_pyramid_build.cu b/src/popsift/s_pyramid_build.cu index 629f633d..26b9f13d 100755 --- a/src/popsift/s_pyramid_build.cu +++ b/src/popsift/s_pyramid_build.cu @@ -198,10 +198,7 @@ inline void Pyramid::downscale_from_prev_octave( int octave, cudaStream_t stream gauss::get_by_2_pick_every_second <<>> - ( POPSIFT_LAYERED_SRC( prev_oct_obj.getDataTexPoint( ), - prev_oct_obj.getDataSurface( ), - prev_oct_obj.getWidth(), - prev_oct_obj.getHeight() ), + ( prev_oct_obj.getDataReadTexPoint( ), prev_oct_obj.getWidth(), prev_oct_obj.getHeight(), _levels-PREV_LEVEL, @@ -231,10 +228,7 @@ inline void Pyramid::horiz_from_prev_level( int octave, int level, cudaStream_t gauss::absoluteSourceInterpolated::horiz <<>> - ( POPSIFT_LAYERED_SRC( oct_obj.getDataTexLinear( ).tex, - oct_obj.getDataSurface( ), - oct_obj.getWidth(), - oct_obj.getHeight() ), + ( oct_obj.getDataReadTexLinear( ), oct_obj.getIntermediateSurface( ), level ); } @@ -250,10 +244,7 @@ inline void Pyramid::horiz_from_prev_level( int octave, int level, cudaStream_t gauss::absoluteSource::horiz <<>> - ( POPSIFT_LAYERED_SRC( oct_obj.getDataTexPoint( ), - oct_obj.getDataSurface( ), - oct_obj.getWidth(), - oct_obj.getHeight() ), + ( oct_obj.getDataReadTexPoint( ), oct_obj.getIntermediateSurface( ), level ); } @@ -291,9 +282,7 @@ inline void Pyramid::vert_from_interm( int octave, int level, cudaStream_t strea gauss::absoluteSourceInterpolated::vert <<>> - ( POPSIFT_LAYERED_SRC( oct_obj.getIntermDataTexLinear( ).tex, - oct_obj.getIntermediateSurface( ), - width, height ), + ( oct_obj.getIntermReadTexLinear( ), oct_obj.getDataSurface( ), level ); } @@ -307,9 +296,7 @@ inline void Pyramid::vert_from_interm( int octave, int level, cudaStream_t strea gauss::absoluteSourceInterpolated::vert_abs0 <<>> - ( POPSIFT_LAYERED_SRC( oct_obj.getIntermDataTexLinear( ).tex, - oct_obj.getIntermediateSurface( ), - width, height ), + ( oct_obj.getIntermReadTexLinear( ), oct_obj.getDataSurface( ), level ); } @@ -323,9 +310,7 @@ inline void Pyramid::vert_from_interm( int octave, int level, cudaStream_t strea gauss::absoluteSource::vert <<>> - ( POPSIFT_LAYERED_SRC( oct_obj.getIntermDataTexPoint( ), - oct_obj.getIntermediateSurface( ), - width, height ), + ( oct_obj.getIntermReadTexPoint( ), oct_obj.getDataSurface( ), level ); } @@ -339,9 +324,7 @@ inline void Pyramid::vert_from_interm( int octave, int level, cudaStream_t strea gauss::absoluteSource::vert_abs0 <<>> - ( POPSIFT_LAYERED_SRC( oct_obj.getIntermDataTexPoint( ), - oct_obj.getIntermediateSurface( ), - width, height ), + ( oct_obj.getIntermReadTexPoint( ), oct_obj.getDataSurface( ), level ); } @@ -377,9 +360,7 @@ inline void Pyramid::vert_all_from_interm( int octave, int start_level, int max_ gauss::absoluteSourceInterpolated::vert_all_abs0 <<>> - ( POPSIFT_LAYERED_SRC( oct_obj.getIntermDataTexLinear( ).tex, - oct_obj.getIntermediateSurface( ), - width, height ), + ( oct_obj.getIntermReadTexLinear( ), oct_obj.getDataSurface( ), start_level, max_level ); @@ -394,9 +375,7 @@ inline void Pyramid::vert_all_from_interm( int octave, int start_level, int max_ gauss::absoluteSource::vert_all_abs0 <<>> - ( POPSIFT_LAYERED_SRC( oct_obj.getIntermDataTexPoint( ), - oct_obj.getIntermediateSurface( ), - width, height ), + ( oct_obj.getIntermReadTexPoint( ), oct_obj.getDataSurface( ), start_level, max_level ); @@ -429,10 +408,7 @@ inline void Pyramid::dogs_from_blurred( int octave, int max_level, cudaStream_t gauss::make_dog <<>> - ( POPSIFT_LAYERED_SRC( oct_obj.getDataTexPoint( ), - oct_obj.getDataSurface( ), - oct_obj.getWidth(), - oct_obj.getHeight() ), + ( oct_obj.getDataReadTexPoint( ), oct_obj.getDogSurface( ), oct_obj.getWidth(), oct_obj.getHeight(), diff --git a/src/popsift/s_pyramid_build_aa.cu b/src/popsift/s_pyramid_build_aa.cu index 1a6edff3..aadf964d 100755 --- a/src/popsift/s_pyramid_build_aa.cu +++ b/src/popsift/s_pyramid_build_aa.cu @@ -33,13 +33,11 @@ __global__ void horiz(LayeredReadTex src_point_texture, cudaSurfaceObject_t dst_ g = filter[span]; out += ( A + B ) * g; - // Horizontal Gauss tap exchange via warp shuffles. The block is (32,blockDim.y): - // each threadIdx.y row is an independent 32-lane group convolving one image - // row. On a 64-lane wavefront two rows share a wavefront, so the shuffles must - // be confined to a width-32 sub-group (threadIdx.x is already the in-row lane - // id since blockDim.x==32) or a lane would pull a neighbour from the wrong - // image row and corrupt the pyramid. CUDA: width 32 is the whole warp. -#if defined(USE_HIP) || defined(__HIP_PLATFORM_AMD__) + // Horizontal Gauss tap exchange via warp shuffles. The block is + // (32,blockDim.y) and each threadIdx.y row convolves one image row. The + // shuffle width is that row width, not the hardware warp size: on a 64-lane + // wavefront two rows share a wavefront and a lane would otherwise pull a + // neighbour from the wrong image row. int shiftval = 0; for( int offset=span-1; offset>0; offset-- ) { shiftval += 1; @@ -52,20 +50,6 @@ __global__ void horiz(LayeredReadTex src_point_texture, cudaSurfaceObject_t dst_ g = filter[offset]; out += ( D + E ) * g; } -#else - int shiftval = 0; - for( int offset=span-1; offset>0; offset-- ) { - shiftval += 1; - const float D1 = popsift::shuffle_down( A, shiftval ); - const float D2 = popsift::shuffle_up ( C, span - shiftval ); - const float D = threadIdx.x < (32 - shiftval) ? D1 : D2; - const float E1 = popsift::shuffle_up ( B, shiftval ); - const float E2 = popsift::shuffle_down( C, span - shiftval ); - const float E = threadIdx.x > shiftval ? E1 : E2; - g = filter[offset]; - out += ( D + E ) * g; - } -#endif surf2DLayeredwrite( out, dst_data, off_x*4, off_y, dst_level, cudaBoundaryModeZero ); } diff --git a/src/popsift/s_pyramid_fixed.cu b/src/popsift/s_pyramid_fixed.cu index c9af8be2..605ab6b9 100755 --- a/src/popsift/s_pyramid_fixed.cu +++ b/src/popsift/s_pyramid_fixed.cu @@ -30,12 +30,11 @@ inline float octave_fixed_horiz( float fval, const float* filter ) * input fval of thread N is extracted from image index N-4 * output fval of thread N should be filtered sum from N-4 to N+4 */ - // Horizontal fixed-span Gauss via warp shuffles. block.x==32 and the block - // packs multiple rows (threadIdx.y/z); on a 64-lane wavefront confine the - // shuffles to a width-32 sub-group so a lane does not pull a neighbour from - // another row. CUDA: width 32 is the whole warp, unchanged. + // Horizontal fixed-span Gauss via warp shuffles. block.x is 32 and the block + // packs several rows in threadIdx.y/z. The shuffle width is that row width, + // not the hardware warp size: on a 64-lane wavefront two rows share a + // wavefront and a lane would otherwise pull a neighbour from another row. float out = fval * filter[0]; -#if defined(USE_HIP) || defined(__HIP_PLATFORM_AMD__) #pragma unroll for( int i=1; i<=SHIFT; i++ ) { float val = popsift::shuffle_up( fval, i, 32 ) + popsift::shuffle_down( fval, i, 32 ); @@ -43,15 +42,6 @@ inline float octave_fixed_horiz( float fval, const float* filter ) } fval = popsift::shuffle_down( out, SHIFT, 32 ); -#else - #pragma unroll - for( int i=1; i<=SHIFT; i++ ) { - float val = popsift::shuffle_up( fval, i ) + popsift::shuffle_down( fval, i ); - out += val * filter[i]; - } - - fval = popsift::shuffle_down( out, SHIFT ); -#endif return fval; } @@ -272,10 +262,7 @@ inline void make_octave_sub( const Config& conf, ImageBase* base, Octave& oct_ob gauss::fixedSpan::absoluteTexAddress::octave_fixed <<>> - ( POPSIFT_LAYERED_SRC( oct_obj.getDataTexPoint( ), - oct_obj.getDataSurface( ), - oct_obj.getWidth(), - oct_obj.getHeight() ), + ( oct_obj.getDataReadTexPoint( ), oct_obj.getDataSurface( ), oct_obj.getWidth(), oct_obj.getHeight(), diff --git a/src/popsift/sift_desc.cu b/src/popsift/sift_desc.cu index c03aee5f..fbc753de 100644 --- a/src/popsift/sift_desc.cu +++ b/src/popsift/sift_desc.cu @@ -39,7 +39,7 @@ using namespace std; * This implies that a keypoint that has changed octave * in subpixelic refinement is going to be sampled from * the wrong level of the octave. - * Unfortunately, we cannot implement getDataTexPoint() + * Unfortunately, we cannot implement getDataReadTexPoint() * as a layered 2D texture to fix this issue, because that * would require to store blur levels in cudaArrays, which * are hard to write. Alternatively, we could keep a diff --git a/src/popsift/sift_octave.cu b/src/popsift/sift_octave.cu index 4ec5c908..d7879983 100755 --- a/src/popsift/sift_octave.cu +++ b/src/popsift/sift_octave.cu @@ -221,13 +221,14 @@ void Octave::alloc_data_planes() _data_ext.height = _h; _data_ext.depth = _levels; - // Observed on gfx90a (ROCm 7.2.1): layered images are incoherent across kernel - // launches (the layer dimension collapses to a single layer on read); filed as - // ROCm/clr#275 (the partial fix ROCm/rocm-systems#6683 covers only surf2DLayered). - // A non-layered 3D array with surf3D/tex3D access is coherent, so drop - // cudaArrayLayered on HIP. The blur levels are addressed by the z coordinate - // instead of the layer index. See cuda_to_hip.h and common/assist.h. CUDA keeps - // a real layered array. + // Observed on gfx90a and on gfx1100 (ROCm 7.2.1): a layered array written + // layer by layer through surf2DLayeredwrite reads back as a single layer for + // every layer index, because the write passed the layer index in the mipmap + // level slot. Filed as ROCm/clr#275; ROCm/rocm-systems#6683 corrects it, but + // that is not in ROCm 7.2.x. A non-layered 3D array with surf3D/tex3D access + // is coherent, so drop cudaArrayLayered on HIP. The blur levels are addressed + // by the z coordinate instead of the layer index. See cuda_to_hip.h, + // sift_textures.h and common/assist.h. CUDA keeps a real layered array. #if defined(USE_HIP) || defined(__HIP_PLATFORM_AMD__) err = cudaMalloc3DArray( &_data, &_data_desc, _data_ext, cudaArraySurfaceLoadStore ); #else @@ -279,15 +280,16 @@ void Octave::alloc_data_tex() tex_desc.readMode = cudaReadModeElementType; // read as float tex_desc.filterMode = cudaFilterModeLinear; // hardware bilinear (CUDA) #if defined(USE_HIP) || defined(__HIP_PLATFORM_AMD__) - // Observed on gfx90a (ROCm 7.2.1): hipCreateTextureObject rejects hardware - // linear filtering on element-read float arrays ("operation not supported"). - // Create the texture with point filtering and do bilinear interpolation in - // software in readTex(). Empirical on this device/ROCm; re-verify on RDNA. + // Observed on gfx90a (ROCm 7.2.1): creating a texture with hardware linear + // filtering over an element-read float array fails with "operation not + // supported"; gfx1100 accepts it, so this is per device. Create the texture + // with point filtering and interpolate in software in readTex(), which keeps + // one build correct on either device. tex_desc.filterMode = cudaFilterModePoint; #endif err = cudaCreateTextureObject( &_data_tex_linear.tex, &res_desc, &tex_desc, nullptr ); - POP_CUDA_FATAL_TEST(err, "Could not create Blur data point texture: "); + POP_CUDA_FATAL_TEST(err, "Could not create Blur data linear texture: "); } void Octave::free_data_tex() @@ -370,7 +372,7 @@ void Octave::alloc_interm_tex() #endif err = cudaCreateTextureObject( &_intm_tex_linear.tex, &res_desc, &tex_desc, nullptr ); - POP_CUDA_FATAL_TEST(err, "Could not create Blur intermediate point texture: "); + POP_CUDA_FATAL_TEST(err, "Could not create Blur intermediate linear texture: "); } void Octave::free_interm_tex() diff --git a/src/popsift/sift_octave.h b/src/popsift/sift_octave.h index ead94c02..e985ba6f 100755 --- a/src/popsift/sift_octave.h +++ b/src/popsift/sift_octave.h @@ -11,21 +11,13 @@ #include "sift_conf.h" #include "sift_constants.h" #include "sift_extremum.h" +#include "sift_textures.h" #include #include namespace popsift { -struct LinearTexture -{ - // Holds a texture object (assigned from cudaCreateTextureObject and consumed - // by tex2DLayered). On CUDA texture and surface handles are both unsigned - // long long so the original cudaSurfaceObject_t typing compiled; HIP uses - // distinct pointer types, so the correct texture type must be used here. - cudaTextureObject_t tex; -}; - class Octave { int _w{}; @@ -105,18 +97,6 @@ class Octave return _desc_done; } - inline LinearTexture getIntermDataTexLinear( ) { - return _intm_tex_linear; - } - inline cudaTextureObject_t getIntermDataTexPoint( ) const { - return _intm_tex_point; - } - inline LinearTexture getDataTexLinear( ) { - return _data_tex_linear; - } - inline cudaTextureObject_t getDataTexPoint( ) const { - return _data_tex_point; - } inline cudaSurfaceObject_t getDataSurface( ) const { return _data_surf; } @@ -127,8 +107,25 @@ class Octave inline cudaSurfaceObject_t& getDogSurface( ) { return _dog_3d_surf; } - inline cudaTextureObject_t& getDogTexturePoint( ) { - return _dog_3d_tex_point; + + /* The read handles the kernels take. A read handle names the array to read + * and the filtering to read it with, and it carries whatever else the + * platform needs to address that array (see sift_textures.h). + */ + inline LayeredReadTex getDataReadTexPoint( ) const { + return makeLayeredReadTex( _data_tex_point, _data_surf, _w, _h ); + } + inline LayeredReadTex getDataReadTexLinear( ) const { + return makeLayeredReadTex( _data_tex_linear.tex, _data_surf, _w, _h ); + } + inline LayeredReadTex getIntermReadTexPoint( ) const { + return makeLayeredReadTex( _intm_tex_point, _intm_surf, _w, _h ); + } + inline LayeredReadTex getIntermReadTexLinear( ) const { + return makeLayeredReadTex( _intm_tex_linear.tex, _intm_surf, _w, _h ); + } + inline LayeredReadTex getDogReadTexPoint( ) const { + return makeLayeredReadTex( _dog_3d_tex_point, _dog_3d_surf, _w, _h ); } /** diff --git a/src/popsift/sift_textures.h b/src/popsift/sift_textures.h new file mode 100644 index 00000000..26e7e3fe --- /dev/null +++ b/src/popsift/sift_textures.h @@ -0,0 +1,71 @@ +/* + * Copyright 2016, 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 + +/* The texture and surface handle types the pyramid hands to the kernels. + */ + +#include + +namespace popsift { + +/* A texture object created with linear filtering, so that a fetch through it is + * distinguishable from a fetch through a point-filtered one at the call site. + */ +struct LinearTexture +{ + cudaTextureObject_t tex; +}; + +#if defined(USE_HIP) || defined(__HIP_PLATFORM_AMD__) +/* The read handle for a pyramid array, passed to the readTex overloads in + * common/assist.h. + * + * On HIP the pyramid arrays are not layered. Observed on gfx90a and on gfx1100 + * (ROCm 7.2.1): once a layered array has been written layer by layer through + * surf2DLayeredwrite, a read in a later kernel launch returns the last written + * layer for every layer index, through tex2DLayered, through surf2DLayeredread + * and through hipMemcpy3D alike. Filed as ROCm/clr#275. The defect is on the + * write side: 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 all three read paths return correct per-layer data, but it + * is not in ROCm 7.2.x. A plain 3D array is coherent across launches, so the + * arrays are allocated non-layered (sift_octave.cu), the layer index becomes the + * z coordinate, and the read path samples the surface. That needs the surface + * and the level dimensions alongside the texture, which is what this type + * carries. + * + * On CUDA the handle is the texture object itself and the read path is + * unchanged. + */ +struct LayeredReadTex +{ + cudaTextureObject_t tex; + cudaSurfaceObject_t surf; + int width; + int height; +}; + +inline LayeredReadTex makeLayeredReadTex( cudaTextureObject_t tex, + cudaSurfaceObject_t surf, + int width, int height ) +{ + return LayeredReadTex{ tex, surf, width, height }; +} +#else +using LayeredReadTex = cudaTextureObject_t; + +inline LayeredReadTex makeLayeredReadTex( cudaTextureObject_t tex, + cudaSurfaceObject_t, + int, int ) +{ + return tex; +} +#endif + +} // namespace popsift