Feature request: CUDA backend for geocodeSlc using existing CUDA Geocode infrastructure

[s-sasaki-earthsea-wizard]

Summary

isce3::geocode::geocodeSlc (the kernel behind Sentinel-1 CSLC generation in COMPASS, and SLC-domain geocoding generally) has no CUDA backend. Profiling on a real Sentinel-1 IW burst (~16 s wall per burst on a 16-core host) shows its time is concentrated in the per-pixel geo2rdr Newton solve — functionality that isce3::cuda::geocode::Geocode already implements on-device and
that the production workflow code (nisar/workflows/geocode_insar.py) already drives.

This issue proposes building a CUDA geocodeSlc backend on the existing cuda::geocode::Geocode infrastructure rather than writing a parallel CUDA implementation from scratch. Whether that takes the form of an extension of Geocode, composition, or a shared lower-level backend is a design question I'd like maintainer input on (below).
An end-to-end trial of the existing CUDA class on the same real burst supports feasibility: it reproduces CPU geometry and amplitude to within quantified tolerances and runs the shared portion ~5× faster than the CPU path on a consumer GPU. I am prepared to implement this and want to align on scope and design constraints first.

Key results

measurement	result
geo2rdr share of the geocodeSlc call (py-spy ×2 + perf)	69–77 %
isce3::core::Orbit::interpolate share of whole-run user cycles	~54 %
Existing CUDA Geocode vs CPU parity, RTX 5080, 1 burst (n=5)	5.1×
Same trial, A100, 9-burst IW2 aggregate (n=3/burst; CPU host-shared)	6.7× warm / 5.7× batch
Full frame, A100 GPU sequential batch (27 bursts, 2.04 Gpx)	62.0 s wall
Valid-pixel masks, RTX 1 burst / A100 9 bursts	identical (42,841,555 px / all 9 bursts)
Amplitude agreement, RTX / A100	r ≈ 1.0; max rel. err 3.0e-2 / worst-burst 4.3e-2

All numbers are hardware- and workload-specific; measurement protocols, fairness notes, and caveats are in the collapsed sections
below. Full provenance: RTX 5080 evidence bundle ·
A100 evidence bundle.

Motivation: where the CPU time actually goes (profiling methodology and attribution)

Workload: one Sentinel-1 IW3 burst (Boso, Japan), geogrid 3225 × 20939 (~67.5 Mpx), biquintic DEM interpolation, SINC SLC interpolation.
Host: Intel Core Ultra 9 285H (16 cores), 96 GB RAM, NVIDIA RTX 5080 16 GiB (measurement provenance and scripts:
RTX evidence bundle).
The geocodeSlc call takes ~16 s wall per burst (array mode, OpenMP, warm caches).

Three measurements agree on the attribution (the parallel py-spy run had worker-0 sampling bias and lag warnings and serves as a
consistency check; the OMP=1 and perf runs are the primary evidence):

method	geo2rdr share of the geocodeSlc call
py-spy --native, parallel run (5079 samples)	76.5 %
py-spy --native, OMP_NUM_THREADS=1 (13884 samples)	69 %
perf record (cycles:u, all threads, 36562 samples)	~77 %

Sharper still: isce3::core::Orbit::interpolate alone is ~54 % of all user cycles of the whole run — the Newton solve's cost is
dominated by per-iteration Hermite orbit interpolation (detail::geo2rdr, Geo2Rdr.icc). The remainder of the kernel is
SINC interpolation (~20 % of single-thread kernel CPU) and DEM interpolation (~10 %). Note the different denominators: the 69-77 % figures are shares of the geocodeSlc call; the 54 % figure is a share of whole-run user cycles. perf stat: IPC ≈ 4.0, branch-miss < 0.7 % — consistent with compute-bound behavior (no cache/bandwidth counters were collected, so stated as "consistent with", not proven).

This is why a CUDA path is worth having: the hot loop is exactly the part cuda::geocode::Geocode already runs on-device (geoToRdrIndices kernel calls the device geo2rdr, all-double).

Feasibility evidence: existing CUDA Geocode on the same burst (trial protocol, timings, fp64, agreement)

To bound what reuse buys without writing any new CUDA code, I ran the existing isce3.cuda.geocode.Geocode (SINC mode) against the CPU isce3.geocode.geocode_slc on the same burst, same geogrid, same DEM, same geo2rdr parameters (delta_range matched to the CPU's hardcoded 1e-8).
Because the CUDA class lacks the SLC-specific phase steps (gap list below), the primary CPU comparator is a "parity" variant with
carrier/Doppler/flattening neutralized, so both sides do the same work; the production-semantics CPU run is also shown.

RTX 5080 (16 GiB, consumer, 1:64 fp64), n = 5, median (IQR):

variant	wall
CUDA Geocode, warm	2.95 s (0.05)
CUDA Geocode, first call	3.45 s
CPU geocode_slc, parity semantics	15.04 s (0.27)
CPU geocode_slc, production semantics	16.27 s (0.58)

≈ 5.1× on the semantic-parity public-API comparison. Notes on fairness: the timed GPU region includes raster read, H2D, kernels,
D2H, and the raster API writes (final close/flush falls outside the timed region), while the CPU region is pure compute on preloaded
arrays — the asymmetry biases against the GPU. This is a production-shaped comparison, not a kernel microbenchmark.
Notably the device geo2rdr is all-double, yet the RTX 5080's weak nominal fp64 ratio (1:64) did not prevent the observed 5.1×
public-API speedup (GPU kernel counters were not collected, so no stronger claim is made). The SLC phase steps to be ported (carrier deramp/reramp, Doppler demodulation, flattening) also need double precision: a standalone kernel study measured ~10 rad unwrapped phase error in fp32 flattening at this slant-range scale, so fp64 is mandatory there. In that study (synthetic 1046×645 grid, flattening kernel only) the fp64 variant ran at 5.5× the fp32 one — far from the 64× nominal throughput gap — but the study covers neither the full phase chain nor burst-scale arrays, so per-step cost on real bursts remains to be measured during the port.

Agreement, measured through the actual code paths (a coordinate-encoded cf32 raster, so indices are quantized by float32):

• Valid-pixel masks identical: 42,841,555 valid pixels on both sides (≈63 % of the geogrid; the remainder lies outside the burst's
  subswath validity region), zero mismatches.
• Radar-index agreement: range max |Δ| 0.0176 px (RMS 1.0e-4), azimuth max |Δ| 0.0011 px (RMS 7.1e-6); p99 = 0 on both axes, but that zero is quantization-dominated (cf32 encoding cannot represent sub-quantum differences), so max/RMS are the informative figures.
• Amplitude (carrier-invariant): Pearson r = 0.99999999999, mean relative error 2.0e-7, p99 = 0; max 3.0e-2 (tail localization
  remains to be characterized).
• Deterministic: bit-exact across repeated runs at fixed lines_per_block.
• Block-size sensitivity: changing lines_per_block 200→100 alters 4 px of 42.8 M (0.093 ppm, max |Δ| 0.0204); the coordinate rasters differ at the same 4 px, locating the origin in the per-block geometry stage — consistent with per-block DEM interpolation or margin handling — rather than in the SLC interpolator. Left as an open question for the shared infrastructure.

Scale-up: datacenter GPU, full subswath and full frame (A100 9-burst / 27-burst results, VRAM scan)

To test whether the single-burst consumer-GPU result holds at production scale and on datacenter hardware, I repeated the trial on
an NVIDIA A100-SXM4-80GB (cloud instance; CPU = 2× AMD EPYC 7742, 128 cores, host-shared — CPU figures below carry that caveat; the RTX numbers above remain the controlled comparison).
Same isce3 source SHA, same trial semantics, lines_per_block = 200 for comparability. Full artifacts:
A100 evidence bundle.

All 9 IW2 VV bursts of the same Sentinel-1 acquisition (geogrids ~3480 × 23300 each, 734 Mpx total), n = 3 per burst:

metric	value
GPU, sum of per-burst warm medians	23.8 s
GPU, sequential batch wall (cold start incl.)	27.7 s
CPU parity semantics, sum of medians	159.2 s
CPU production semantics, sum of medians	166.1 s
speedup (parity, warm)	6.7×
speedup (parity vs batch wall)	5.7×

The geometry/amplitude agreement replicates on every burst of the subswath: valid masks identical on all 9 (Jaccard 1.0), range index max |Δ| ≤ 0.0195 px, azimuth ≤ 0.0013 px, amplitude r ≈ 1.0 with worst per-burst max relative amplitude error 4.3e-2 (per-burst maxima 0.4-4.3e-2; same tail character as the single-burst run).
A full frame (27 VV bursts across IW1/IW2/IW3, 2.04 Gpx) geocodes in a 62.0 s sequential batch wall (~33 Mpx/s sustained), demonstrating the workload at frame scale.

The 80 GB card also allows a lines_per_block scan that doubles as memory evidence for the chip-buffer issue (collapsed section further below).
Wall time and peak device memory (out-of-process busy-poll NVML sampling, ~58k queries/s; coarser 20-50 ms samplers can miss or severely understate this transient — its >4 GiB dwell is only ~40-215 ms, and sampler cadence and process scheduling both contribute):

lines_per_block	wall (s)	measured peak (MiB)	chip model (GiB)
200	2.37	4,849	2.84
1000	1.59	16,951	14.19
2000	1.49	32,049	28.39

Two conclusions: (a) baseline-subtracted peak scaling is consistent with the chip-buffer model (n_elem_out × 81 × 8 B, SINC_ONE² = 81; baseline 1,565 MiB) — at lines_per_block = 1000 (the runconfig-level default in adjacent workflows) the measured allocation footprint of 16,951 MiB exceeds a 16 GiB consumer card's total VRAM, so the OOM hazard on such cards is evidence-backed rather than estimated (no OOM was provoked on the RTX itself); (b) larger blocks do pay on large-memory cards (2.37 s → 1.49 s), with most of the gain already at 1000 — supporting memory-budget-based/adaptive block sizing rather than any fixed default. Outputs were bit-exact across all scanned block sizes on this burst; the 4-px block-size sensitivity observed on the RTX burst did not reproduce here.

Proposal: build on the cuda::geocode::Geocode infrastructure, don't fork it

The shared machinery is already there and driven by the production workflow code (geocode_insar.py): per-pixel
device geo2rdr with azimuth/range correction LUTs, projection + DEM interpolation, masked block machinery with the same interpolation margin as the CPU path, and a SINC interpolator with the same SINC_LEN/SINC_SUB and thrust::complex<float> support.

What geocodeSlc needs that the CUDA class lacks today (from a line-by-line diff of geocodeSlc.cpp vs Geocode.cu); this is the implementation scope:

1. Carrier phase deramp before interpolation (per-radar-pixel az/rg carrier evaluation).
2. Doppler demodulation inside the chip load + add-back after interpolation (no-op for zero-Doppler NISAR, required for S1
   native-Doppler data).
3. Carrier reramp + flattening, and the carrier/flattening phase output rasters (+ uncorrectedSRange option).
4. slicedRadarGrid masking.
5. Subswath validity: CPU checks all chip corners, CUDA checks the chip center only.
6. Geo-grid mask export (subswath labels) — currently internal-only.
7. Carrier API: the CUDA class currently has no carrier arguments at all (its LUT2d surface covers native Doppler and correction LUTs) — a new carrier API is required. CPU geocodeSlc accepts both LUT2d and Poly2d carriers; the existing device LUT2d plumbing may be reusable for the LUT2d case.

Semantic differences that should become explicit porting requirements with acceptance tests (they currently differ silently):

• Correction-LUT out-of-bounds: CPU skips the correction and keeps the pixel; CUDA masks the pixel invalid.
• Non-converged geo2rdr pixels: CPU leaves NaN indices; CUDA writes 0.0 + mask. (The porting requirement should target public-API mask/output semantics, not internal representation.)
• geo2rdr tolerance: CPU hardcodes delta_range = 1e-8; CUDA takes it as a parameter.

Pre-existing sizing issue the port should address (SINC chip-buffer transient; OOM hazard on 16 GiB cards)

The SINC path preallocates n_elem_out × SINC_ONE² chips per block (self-described "VERY inefficient" in the source). At this geogrid width, lines_per_block = 1000 implies a ~13-15 GB device transient for cf32 — OOM territory on 16 GiB consumer cards. High-rate NVML sampling on an A100 (scale-up section above) recorded a 16,951 MiB total allocation peak (baseline 1,565 MiB) at lines_per_block = 1000 on a 23.5k-wide geogrid — scaling consistent with the model, so the hazard is evidence-backed rather than estimated; the transient's dwell is ~tens-to-hundreds of ms, so coarse samplers understate it. Memory-budget-based or adaptive block sizing — or a chip-free interpolation path — would remove a real footgun.

What this would mean for users (bounded expectations: 4-5× hypothesis, workflow-level gzip floor)

Bounded, honest framing:

• Phase-level: the existing shared path measured 5.1× (consumer GPU). This is feasibility evidence for the shared infrastructure, not a forecast for a completed SLC backend. As an Amdahl-style hypothesis: the SLC-specific steps account for ~7.6 % of CPU kernel wall here (production 16.27 s vs parity 15.04 s), so an effective 4-5× is plausible if they port with comparable efficiency — to be verified, not promised.
• Workflow-level: for COMPASS-style CSLC runs the geocoding phase also contains single-threaded gzip HDF5 output (~14 s of a ~30 s phase on my host, with the current serial workflow and gzip-4 settings). That is an output-configuration concern outside isce3 and is reported separately; with those settings a CUDA kernel alone cannot shrink the phase below that ~14 s floor (overlapping compression with GPU compute could change this in the future). Stated so this proposal doesn't overclaim.

Maintenance posture and design questions

• Builds on the maintained CUDA class rather than adding a parallel implementation; CPU/CUDA share semantics via the acceptance tests above.
• Tests follow the existing ctest GPU-label separation so CI without a GPU is unaffected.
• Design questions I'd like maintainer input on: (a) whether SLC semantics belong in Geocode itself or a derived/composed API;
  (b) whether the Poly2d carrier path is required for the first increment; (c) for each semantic difference listed above, align-to-CPU or document-and-keep — e.g. should the extension default delta_range to the CPU path's 1e-8 (keeping the parameter), and should correction-LUT out-of-bounds skip the correction (CPU behavior) or mask the pixel (current CUDA behavior)?