Performance: per-pixel dynamic scheduling and unneeded atomics in RTC `_normalizeRtcArea` cost ~11% of GCOV geocode time at NISAR scale

Summary

_normalizeRtcArea in cxx/isce3/geometry/RTC.cpp runs its gamma-naught normalization loop with omp parallel for schedule(dynamic) collapse(2) (no chunk size, i.e. one iteration per dispatch) and performs a per-pixel omp atomic write/update. Each (i, j) is only touched by its own iteration, so the loop is race-free with or without the atomics — they are pure overhead — and the pixel-granularity dynamic scheduling turns a memory-bound pass over two float arrays into ~1.2e9 contended acquisitions of the shared iteration counter per call pair.

On a real NISAR L1 RSLC (frequency A, 29240 x 21232 radar grid) processed through the GCOV workflow on a 16-core CPU host, replacing this with a plain row-wise omp parallel for and a conditional expression:

• cuts the RTC-AP journal timer by 16.5% (median 966.2 s -> 806.5 s),
• cuts the enclosing GEO-AP timer (the single C++ geocodeAreaProj call, ~89% of workflow wall) by 11.1% (1506.1 s -> 1338.8 s, i.e. -167 s),
• reduces total workflow wall by 6.1% (the remainder of the wall is dominated by unrelated I/O),
• leaves all four GTiff output products bit-identical.

Problem

The loop as of current develop (RTC.cpp#L268-L286 @ bdf1f6f):

// cxx/isce3/geometry/RTC.cpp, _normalizeRtcArea
_Pragma("omp parallel for schedule(dynamic) collapse(2)")
    for (int i = 0; i < numerator_array.length(); ++i)
        for (int j = 0; j < numerator_array.width(); ++j) {
            const float denominator_value = denominator_array(i, j);
            if (denominator_value == 0) {
                _Pragma("omp atomic write")
                    numerator_array(i, j) =
                        std::numeric_limits<float>::quiet_NaN();
                continue;
            }
            _Pragma("omp atomic update")
                numerator_array(i, j) /= denominator_value;
        }

Two independent issues compound here:

1. Scheduling granularity. schedule(dynamic) without a chunk size dispatches one iteration at a time; with collapse(2) that is one pixel at a time. Every dispatch is an atomic fetch-and-add on a shared counter contended by all threads. The per-iteration work (one load, one compare, one divide, one store) is far smaller than the dispatch cost. The GCOV workflow calls this function twice per frequency (gamma-naught normalization and the gamma-to-sigma factor), so the frequency-A grid above pays ~1.24e9 dispatches.
2. Unnecessary atomics. The collapse(2) iteration space maps 1:1 onto pixels; no two iterations touch the same (i, j). The atomics do not guard any actual race — the loop is already deterministic without them — and disassembly of the compiled object shows they lower to lock cmpxchg CAS loops on every pixel.

The dispatch and atomics together account for ~20% of the samples in a whole-run CPU profile (see below), attributed to gomp_iter_dynamic_next plus the _normalizeRtcArea outlined function, for a pass that is otherwise memory-bandwidth bound (~7.4 GB of traffic against a >2 GB/s-per-core memory system — sub-second territory if scheduling overhead were removed).

Possible fix

    // Each (i, j) is read and written by exactly one iteration, so no
    // atomics are required. Row-wise static scheduling avoids per-pixel
    // dynamic dispatch overhead.
    _Pragma("omp parallel for")
        for (int i = 0; i < numerator_array.length(); ++i)
            for (int j = 0; j < numerator_array.width(); ++j) {
                const float denominator_value = denominator_array(i, j);
                numerator_array(i, j) = denominator_value == 0
                        ? std::numeric_limits<float>::quiet_NaN()
                        : numerator_array(i, j) / denominator_value;
            }

Row count (tens of thousands) vastly exceeds thread count and per-pixel cost is uniform, so static row-wise scheduling has no load-balance downside here. All measurements below were taken with exactly this change applied — happy to open a PR if this direction looks reasonable to you.

Measurement

Setup: 16-core x86-64 host, GCC 13.3, RelWithDebInfo, isce3 built from develop @ bdf1f6f; NISAR L1 RSLC frequency-A granule (L-band, 29240 x 21232), GCOV CPU workflow (gpu_enabled: False), GTiff output mode, OMP_NUM_THREADS=16. BEFORE = clean develop build; AFTER = same base plus the one-function patch above. Runs were interleaved A -> B -> A (3 timing runs per phase, so 6 BEFORE / 3 AFTER) to control for environment drift; the two BEFORE phases agree within 1.9% (median), and each phase also produced one whole-run perf record -e cpu-clock:u -F 99 profile (BEFORE/AFTER).

Journal timers, medians (all 9 runs exited 0):

timer	BEFORE (n=6, pooled)	AFTER (n=3)	delta
RTC-AP	966.2 s	806.5 s	-16.5%
GEO-AP (contains RTC-AP)	1506.1 s	1338.8 s	-11.1%
workflow wall	28:28	26:43	-6.1%

Per-run RTC-AP timings (all 9 runs)

run	A1 BEFORE	B AFTER	A2 BEFORE
run1	975.2 s	780.5 s	921.2 s
run2	991.5 s	806.5 s	957.2 s
run3	902.2 s	834.5 s	987.6 s
median	975.2 s	806.5 s	957.2 s

Phase medians A1 vs A2 agree within 1.9% (no drift over the ~10 h campaign). Pooled BEFORE (n=6): RTC-AP 966.2 s, GEO-AP 1506.1 s.

The BEFORE and AFTER distributions do not overlap (slowest AFTER RTC-AP 834.5 s vs fastest BEFORE 902.2 s). Total workflow wall improves less in relative terms because ~11% of wall is Python-side I/O wait that the patch does not touch. (The GEO-AP delta, -167.3 s, nominally exceeds the RTC-AP delta, -159.7 s; the patch only touches code inside RTC-AP, and the difference is within the run-to-run spread of either timer.)

Whole-run profile (cpu-clock:u, F=99), share of user CPU:

symbol	BEFORE	AFTER
gomp_iter_dynamic_next (libgomp)	10.91%	absent (<0.5%)
_normalizeRtcArea [omp_fn]	9.86%	absent (<0.5%)
isce3::core::Orbit::interpolate	24.39%	31.12%

The two overhead symbols vanish; other symbols' shares rise only because total CPU time shrank (15021 -> 12018 thread-seconds; Orbit::interpolate in absolute terms is 3664 vs 3740 thread-seconds, unchanged within noise). In BEFORE, dispatch + normalize together are ~3120 thread-seconds, i.e. ~195 s of wall on 16 threads — consistent with the timer deltas.

Output equality: the four GTiff products (HHHH, HVHV, numberOfLooks, rtcGammaToSigmaFactor) are bit-identical across all 9 runs (single md5 per product). This is expected by construction, not just observed: the schedule clause changes which thread executes which iteration and in what order, never which iteration owns which pixel; collapse(2) flattens the nest 1:1 with no tiling or reduction; and each pixel's value is produced by a single IEEE division (or a constant qNaN store). The same argument applies to the unpatched loop — which is why the atomics were pure overhead rather than a correctness device.

Three scope notes for fair reading:

• Measured on one granule / one 16-thread host. The overhead is contention on a shared counter, so the recoverable fraction should grow with thread count, and shrink on smaller grids in absolute terms (though in a frequency-B profile the dispatch share reached ~39% of user CPU — single run, share not timing, run aborted during HDF5 output — because the normalize pass scales with radar grid pixels while the area-projection integration scales with geogrid pixels x upsampling^2).
• The patched loop remains scalar with GCC 13 at -O2/-O3 (the NaN guard defeats if-conversion), so the measured gain is attributable entirely to removing the dynamic dispatch and the atomics. Vectorizing the division (unconditional divide + separate NaN-mask pass) is possible future work, independent of this report.
• This change removes scheduling/synchronization overhead, not compute. The remaining ~800 s of RTC-AP is genuine computation (orbit interpolation, per-pixel radar geometry, the area scatter-add), so the case for accelerating the area-projection RTC by other means — algorithmic or hardware — is essentially unchanged by this fix.

Out-of-scope siblings

One footnote on siblings, deliberately left out of scope: _applyRtcMinValueDb (RTC.cpp:255) uses the identical pixel-granularity schedule(dynamic) collapse(2) + atomic-write pattern and would take the same fix, but it did not execute in the measured configuration. The sigma-naught-ellipsoid loop near RTC.cpp:1032 is not the same case: its per-iteration work is heavy and variable (orbit interpolation + Newton iterations), so dynamic scheduling is plausibly justified there, and the atomics in the adjacent area-accumulation are a genuine scatter-add that must stay.

Read cxx/isce3/geometry/RTC.cpp at _normalizeRtcArea and compare the current loop with the proposed row-wise OpenMP version. Run the GCOV CPU workflow with the stated NISAR configuration, then verify that RTC-AP improves and the four GTiff products (HHHH, HVHV, numberOfLooks, rtcGammaToSigmaFactor) remain bit-identical.