[ROCm] HIP port for AMD GPUs

[jeffdaily]

Ports PopSift (GPU SIFT) to ROCm/HIP. A single compatibility header src/popsift/cuda_to_hip.h aliases the CUDA spellings to HIP under USE_HIP and is a no-op include of the CUDA headers otherwise; the .cu sources build as HIP, and the CUDA build path is byte-for-byte unchanged (every change is wrapped in USE_HIP / if(NOT USE_HIP)). README.md documents the USE_HIP build option, the build recipe, and the supported AMD targets.

Getting correct, deterministic SIFT output required addressing five distinct HIP/ROCm fault classes, each behind a USE_HIP guard.

1. Hardware linear-filter textures. HIP rejects a cudaFilterModeLinear + element-read float texture at creation, so the linear fetches do manual bilinear interpolation in software (point-sample the 4 neighbors and lerp, CUDA's -0.5 texel-center convention).

2. Layered cudaArray dimension collapse. A multi-layer cudaArrayLayered array written one layer per kernel launch collapses to a single, last-written layer on any later-launch read on this ROCm (tex2DLayered, surf2DLayeredread, and even host hipMemcpy3D all return the last layer for every index; see Layered cudaArray collapses to one layer: surf2DLayeredread/tex2DLayered return the last-written layer for every layer index (gfx90a) ROCm/clr#275, which PopSift motivated, and the partial fix [HIP] Fix incorrect ockl binding for layered surface functions JIRA ID : SWDEV-597321 ROCm/rocm-systems#6683 covering only surf2DLayered). The Gaussian pyramid and DoG arrays are converted to non-layered 3D arrays (surf3Dwrite/surf3Dread/tex3D with the level as a real z coordinate, done once in the compat header's write wrapper plus the texture-fetch helper), which is fully coherent across launches.

3. Wavefront-width-generic warp collectives. The bitonic sort, exclusive prefix sums, descriptor reductions and orientation ballots operate per warpSize-lane group rather than assuming a 32-lane warp. The extrema counter in s_extrema.cu is the critical case: the in-wavefront lane is the block's linear thread index modulo the actual warpSize, the ballot/popc/shuffle widths follow warpSize, and the leader is lane 0 of the real wavefront. On wave64 this reduces exactly to a 64-bit ballot with an exclusive-prefix slot assignment; on wave32 it is correct where a fixed 64-lane geometry would have lost atomicAdds and read phantom lanes.

4. RootSift and L2 NaN. The default RootSift normalization sqrt(bin/L1) and the L2 path are guarded against sqrt of a slightly-negative bin and the 0*inf of an all-zero descriptor.

5. Example linkage. The PopSift library links roc::rocthrust PRIVATE and exposes hip::host PUBLIC, so the plain-C++ examples (main.cpp/match.cpp/pgmread.cpp, which have no device code) compile as CXX and just link the HIP runtime, rather than being forced through the HIP compiler by a leaked hip::device interface. The compat header's forward-compat <hip/hip_bf16.h> include is gated on clang so a gcc host consumer never parses that clang-only header; it is pulled in before the _shfl_sync compat macros so newer ROCm's real _shfl_sync<...> functions are defined before the macros (guarded with __has_include so older ROCm, where the header has no such functions, is a harmless no-op).

Tested on three AMD targets. PopSift ships no CPU reference, so the gates are run-to-run determinism, value sanity (finite, correctly normalized descriptors; in-bounds keypoints), and cross-architecture agreement between gfx90a and gfx1100 on the same images.

• gfx90a (AMD Instinct MI250X, Linux, ROCm 7.2.1): popsift-demo gives 895 features / 1494 descriptors, identical across 5/5 runs; 0 NaN/Inf in all 191232 descriptor values; per-descriptor L2 in [0.9993, 1.0009], mean 1.0000; a rotated image gives 896/1484, also deterministic. popsift-match produces real finite distances and a sane ratio test.
• gfx1100 (AMD Radeon Pro W7800, Linux, RDNA3 wave32, ROCm 7.2.1): the Oxford feature-detection images reproduce the gfx90a reference feature counts across all six images, deterministic across runs, 0 NaN/Inf across all descriptor modes.
• gfx1151 (AMD Radeon 8060S APU, Windows, TheRock ROCm): SIFT on a real image gives non-zero, deterministic features (187/221 across runs), 0 NaN/Inf, RootSift L2 ~1.0, keypoints in bounds.

Authored with the assistance of Claude (Anthropic).

[gemini-code-assist]

Code Review

This pull request adds support for AMD GPUs via ROCm/HIP to the PopSift library, enabling compilation and execution on AMD hardware while preserving the existing CUDA path. It introduces a CUDA-to-HIP compatibility shim, addresses warp-size differences (32-lane vs. 64-lane) in reduction and shuffle operations, and works around ROCm-specific limitations regarding layered image coherency and hardware bilinear filtering. The review feedback highlights two key areas for improvement: addressing a potential NaN corruption in s_desc_norm_rs.h when sum is extremely small, and initializing write_index in s_extrema.cu to prevent compiler or static analysis warnings regarding uninitialized variables.

Important

The consumer version of Gemini Code Assist on GitHub is being sunset. Starting June 18, 2026, new organization installations will be blocked, and all code review activity will officially cease on July 17, 2026.
For more details on the timeline and next steps, please review the Help Documentation.

[jeffdaily]

On s_extrema.cu: agreed, write_index will be initialized; __shfl reads only lane 0 so it was harmless, but the initializer silences the warning.

On s_desc_norm_rs.h: the 0 * inf = NaN case is already handled -- the result feeds fmaxf(descr.x * inv, 0.0f), and fmaxf returns 0.0f when one argument is NaN, so a zero bin normalizes to 0 rather than NaN (validation over 191232 descriptor values showed no NaN or Inf). The real residual is that a subnormal sum makes 1.0f/sum overflow to +inf, and a positive bin then yields +inf, which fmaxf(., 0) does not clamp. Guarding the divisor (treating a near-zero sum as an all-zero descriptor) closes that, so I've applied the threshold.

[griwodz]

Hi Jeff,

I'm really thankful that you have made the HIP port!

I'm working abroad this year and I have only a Mac with me. I'll try to get my fingers of an NVidia and an AMD machine from my hosts, but I don't know if I can get it.

The develop branch doesn't have a CPU-only version as you write, but there is the "sycl" branch, which can support CPU only including Mac via SYCL. However, that version is stripped down to one or at most two algorithms in every step of the pipeline.

But in any case, our reference is always vlfeat. For a very long time, I've been annoyed that the vlfeat command line has stopped worked. But a few months ago, I finally found the bug in vlfeat. The README for fixing the bug and the python/opencv scripts for running the tests are 3 commits in the sycl branch.
Commit 80bb8ab gives PopSift a compatible output, but this may not fit the develop branch 100%
Commit 1a7b3f9 adds the Python scripts
Commit 7e3fccc adds a README for fixing the vlfeat command line.

Another note on SYCL: Intel's SYCL compiler doesn't do very well with the PopSift port, while U.Stuttgart's compiler doesn't support hardware interpolation via sycl::image. CUDA and HIP still win, although the SYCL port shows that PopSift's DoG and normalize code could be better.

[jeffdaily]

Thanks Carsten, and thanks for the vlfeat pointers.

I ran your vlfeat comparison on an AMD GPU (gfx90a) using your own setup: vlfeat built from source with the vl/pgm.c fix from your README, popsift-demo with --root-sift --write-as-uchar --write-as-orientation --norm-multi 9 --threshold 0, and your testScripts/openCVMatches.py. On Oxford boat img1: vlfeat produces 14397 descriptors, popsift on AMD produces 14509, and the cross-checked L2 BFMatcher finds 13914 mutual matches. 99.6% of those matched pairs land within 0.5px of the same keypoint (median offset 0.0px), so it is genuine agreement with the vlfeat reference rather than a descriptor-space coincidence. I applied your 80bb8ab output-format change by hand since on develop the descriptor printing lives in features.cu rather than features.cc, but it was otherwise a clean port of your flag.

On the AMD side this is already running on gfx90a, gfx1100 (RDNA3), and gfx1201 (RDNA4) through the same SIFT pipeline. The AMD support is entirely behind USE_HIP and the CUDA path is unchanged; I also confirmed the NVIDIA build still compiles clean under CUDA 12.6 (the full library, sm_86), so existing CUDA users are unaffected. You don't need to source an AMD machine to take this. If it's useful I can also send the side-by-side keypoint image that openCVMatches.py produces.

One unrelated note while I was in there: develop does not build under CUDA 13, because the newer Thrust no longer transitively includes <thrust/tuple.h> and s_filtergrid.cu then fails on thrust::get/make_tuple. It is pre-existing (the unmodified file fails the same way) and trivial to fix with an explicit include.

[griwodz]

I still don't have access to an AMD machine, but at least I can accress my NVidia machine again.
These are my current comments.

[griwodz]

I've proposed in PR #190 to remove support for CUDA 10 and earlier. C++17 can be the default for CUDA as well.

[griwodz]

I've always found that the CUDA docs are missing a good explanation for this +0.5f offset. I believe that the weird CUDA addressing exists because normalized addressing requires that the left-/topmost corner or the left-/topmost pixel have index 0 (instead of the pixel center). If index 0 were the center of the pixel, normalized addressing would have to work on negative coordinates.

[griwodz]

Since x is always positive, x - floorf(x) is always the same as x-static_cast(x)
For some reason, CUDA doesn't have a floor intrinsic and the cast is faster. Maybe HIP is better?

[griwodz]

Have you considered using fmaf() ?
It was a real time-saver in the s_desc* functions.
Or perhaps the equivalent of cuda::std::lerp() exists for HIP.

[griwodz]

That is very unfortunate. The main reason for going through the hassle with textures and surfaces was to benefit from the hardware interpolation. If you have to do that manually, wouldn't it be more efficient to drop textures altogether and use "normal" cudaMalloc?
Also getting rid of the half-pixel offset.

[griwodz]

I'm not clear how this happens.
The CUDA code uses one layer of one texture backed by an array for reading, and it uses a different layer of the same surface backed by the same array for writing.
Layers are never touched by one another. Since you don't have hardware interpolation, you could use a simple 3D point texture.
Does the problem appear because the reading from a texture and writing to a surface backed by the same Array is undefined in ROCm even when the cache lines never touch? That is inconvenient.

But I was wondering: I want a cudaArray for CUDA because I can only interpolate in X/Y while remaining discrete in Z when I used a Layered2D instead of 3D. And I can only define Layered2D over a cudaArray. If linear interpolation is not working anyway, you could use plain 3D memory, and the caching would happen in L1. Although that would cost read performance unless you transpose the Intermediate layers.

[griwodz]

PS: It could be a really good idea to create an alternative Octave class that works on flat memory (3D or Pitched). That alternative would have transposed Intermediate planes and need no texture at all. It might make the HIP code easier to read, and it would overcome the CUDA limit for large images, which are today constrained by the maximum Texture dimensions.

[griwodz]

oops

[griwodz]

My original code was buggy. The final commit shouldn't keep the comment because it would confuse future readers.

[griwodz]

The need for rdc / CUDA_SEPARABLE_COMPILATION is an artifact of a mistake I made a long time ago: the host code that calls a global function should always be in the same compile unit but it isn't.
There is a branch that removes a multitude of alternatives that don't work well, where this mistake is also fixed.
I'm afraid it has to stay in the CUDA case for now.

[griwodz]

My comment string for the Tegra case is wrong. CMake is unable to choose the right CCs for Tegra. If we don't force it, it will configure for discrete-GPU CCs also on ARM hosts. But Tegra-family machines cannot host discrete NVidia GPUs.

[jeffdaily]

Thanks for the careful review. All comments are addressed on the branch; here is a summary.

Every shuffle now passes its width explicitly, on both platforms, so the #ifdef USE_HIP pairs around the shuffles are gone. The widths are the algorithmic group widths (a descriptor row, the bitonic network, block.x), which is why they are literal 32s rather than warpSize.

The shuffle compatibility now goes through PopSift_HAVE_SHFL_DOWN_SYNC as you suggested: it is 0 on the HIP build, so assist.h takes its existing pre-CUDA-9 branch, which is exactly the spelling HIP provides. The seven __*_sync macros are gone from the compatibility header, and with them the bf16 include-order workaround they had caused.

The texture and surface types now live in one header, src/popsift/sift_textures.h (LayeredReadTex and LinearTexture), and the layered-source macro is replaced by a makeLayeredReadTex() helper there.

Octave now returns filled read handles (getDataReadTexPoint/Linear, getIntermReadTexPoint/Linear, getDogReadTexPoint) built from its own members, and all 22 launch sites use them. Since you wrote "instead of Octave::getDataTexPoint()", the five superseded raw accessors are removed; say the word if you want them kept and I will restore them.

The RootSift normalization has the structure you wrote: the early if (inv <= 0.0f) zeroing and the per-bin ternary with scalbnf. One note: the small-sum threshold (sum > 1e-20f) is kept, and now applies on both platforms. It was added earlier in this PR because a subnormal sum makes __fdividef(1.0f, sum) overflow to infinity; your inv <= 0.0f test is unaffected by it. The NVIDIA path previously divided each bin by the sum directly, so this is a tiny behavioral change on that path, called out in the commit message.

The L2 normalization guard is one line on both paths: norm = (norm > 0.0f) ? __frsqrt_rn(norm) : 0.0f.

Also done: the copy-pasted "point texture" error string on the linear-texture creation, the stale comments (including recording that gfx1100 accepts hardware linear filtering, so the software bilinear is a per-device workaround kept so one binary serves both), a shared Thrust setup header (common/thrust_setup.h), and fmaf in the manual bilinear lerp (the compiler had already contracted it, so the binary is unchanged).

One suggestion I did not take: x - static_cast<int>(x) for x - floorf(x). The x there is not always positive: the horizontal filter reads readTex(src, off_x - span, ...) and the gradient reads x - 1.0f, so x goes a few texels negative at the image border, where truncation and floor differ. floorf keeps the clamped-border behavior.

On the other points: agreed that rdc stays; the CUDA 10 / C++17 CMake modernization makes sense in #190 rather than here; and replacing the texture path with flat 3D memory and a leaner Octave class reads like the right long-term direction, which this change deliberately does not start.

Since your review the branch was validated again end to end on an MI250X (gfx90a) and a Radeon 8060S (gfx1151, Windows), with earlier passes of the same changes on a Radeon Pro W7800 (gfx1100) and two more Windows machines (gfx1201 and gfx1101); the CUDA arm compile-checks clean with nvcc 12.8. On gfx90a one of the six benchmark images yields 9451 descriptors where the other architectures give 9452, stable across runs -- a single boundary keypoint on the 64-wide wavefront, noted here for completeness.

[jeffdaily]

Note: this reply was drafted by an AI assistant.

Merged develop to clear the conflict on src/popsift/s_filtergrid.cu. The overlap was your CUDA 13.3 include fix in #188 landing on the same include block this branch had reworked to route through common/thrust_setup.h. Kept both: your two headers in their sorted positions, our shared-header include. thrust_setup.h only pulls in execution_policy.h and version.h, so nothing is duplicated.

Rebuilt and re-ran the correctness checks after the merge. The Oxford boat feature and descriptor counts are unchanged on every image, and the CUDA build is still clean with nvcc 12.8 (sm_86, 0 errors, 0 warnings).