diff --git a/libvmaf/src/feature/cuda/integer_adm/adm_cm.cu b/libvmaf/src/feature/cuda/integer_adm/adm_cm.cu
index ea3d88de..e368402d 100644
--- a/libvmaf/src/feature/cuda/integer_adm/adm_cm.cu
+++ b/libvmaf/src/feature/cuda/integer_adm/adm_cm.cu
@@ -61,39 +61,34 @@ extern "C" {
int32_t accum_thread = 0;
if (i < end_row && j < end_col) {
- int16_t offset_i[2] = {-1, 1};
- if (i == 0 && top <= 0) {
- offset_i[0] = 1;
- } else if (i == (h - 1) && bottom > (h - 1)) {
- offset_i[1] = 0;
- }
-
- int16_t offset_j[2] = {-1, 1};
- if (j == 0 && left <= 0) {
- offset_j[0] = 1;
- } else if (j == (w - 1) && right > (w - 1)) {
- offset_j[1] = 0;
- }
+ /* Border handling identical to the CPU I4_ADM_CM_THRESH_S_* macros:
+ * the top/left borders mirror {1,0,1}, the
+ * bottom/right borders replicate {n-2,n-1,n-1}; both patterns are
+ * min(abs(p), n-1). Interior positions pass through unchanged, and
+ * no read can land outside [0,h-1]x[0,w-1] (the old offset_i
+ * scheme read flt rows {1,2,3} at i==0 -- wrong mirror, with row 3
+ * one past the csf-written region). */
+ const int pos_i[3] = {min(abs(i - 1), h - 1), i, min(i + 1, h - 1)};
+ const int pos_j[3] = {min(abs(j - 1), w - 1), j, min(j + 1, w - 1)};
int32_t thr = 0;
for (int theta = 0; theta < 3; ++theta) {
int32_t sum = 0;
- int32_t src = angles[theta][src_stride * (i + offset_i[0] + 1) + j];
- int32_t *flt_ptr = flt_angles[theta];
- flt_ptr += (src_stride * (i + offset_i[0]));
- sum += flt_ptr[j + offset_j[0]];
- sum += flt_ptr[j];
- sum += flt_ptr[j + offset_j[1]];
- flt_ptr += src_stride;
- sum += flt_ptr[j + offset_j[0]];
+ int32_t src = angles[theta][src_stride * i + j];
+ int32_t *flt_ptr = flt_angles[theta] + src_stride * pos_i[0];
+ sum += flt_ptr[pos_j[0]];
+ sum += flt_ptr[pos_j[1]];
+ sum += flt_ptr[pos_j[2]];
+ flt_ptr = flt_angles[theta] + src_stride * pos_i[1];
+ sum += flt_ptr[pos_j[0]];
sum += (int32_t)((((int64_t)I4_ONE_BY_15 * abs((int32_t)src)) +
add_bef_shift_flt) >>
shift_flt);
- sum += flt_ptr[j + offset_j[1]];
- flt_ptr += src_stride * offset_i[1];
- sum += flt_ptr[j + offset_j[0]];
- sum += flt_ptr[j];
- sum += flt_ptr[j + offset_j[1]];
+ sum += flt_ptr[pos_j[2]];
+ flt_ptr = flt_angles[theta] + src_stride * pos_i[2];
+ sum += flt_ptr[pos_j[0]];
+ sum += flt_ptr[pos_j[1]];
+ sum += flt_ptr[pos_j[2]];
thr += sum;
}
int32_t x = (int32_t)((((int64_t)src_band[i * src_stride + j] *
@@ -115,15 +110,6 @@ __constant__ const int fixed_shift[3] = {4, 4, 3};
__constant__ const int32_t shift_xsq[3] = {29, 29, 30};
__constant__ const int32_t add_shift_xsq[3] = {268435456, 268435456, 536870912};
-// HACK: the 256 byte alignment is required to ensure that the struct is not moved to lmem
-struct WarpShift
-{
- uint32_t shift_cub[3];
- uint32_t add_shift_cub[3];
- uint32_t shift_sq[3];
- uint32_t add_shift_sq[3];
-};
-
template <int rows_per_thread>
__device__ __forceinline__ void adm_cm_line_kernel(AdmBufferCuda buf, int h, int w, int top,
int bottom, int left, int right,
@@ -131,15 +117,7 @@ __device__ __forceinline__ void adm_cm_line_kernel(AdmBufferCuda buf, int h, int
int end_col, int src_stride,
int csf_a_stride, int buffer_h,
int buffer_stride, int32_t *accum_per_block,
- AdmFixedParametersCuda params,
- // reduce
- int scale, int64_t* accum_global,
-
- // shift warp
- WarpShift ws,
- // shift global
- const uint32_t shift_inner_accum, const uint32_t add_shift_inner_accum
- ) {
+ AdmFixedParametersCuda params) {
const cuda_adm_dwt_band_t *src = &buf.decouple_r;
const cuda_adm_dwt_band_t *csf_f = &buf.csf_f;
const cuda_adm_dwt_band_t *csf_a = &buf.csf_a;
@@ -220,37 +198,31 @@ __device__ __forceinline__ void adm_cm_line_kernel(AdmBufferCuda buf, int h, int
accum_thread_reg[row] = max(0, sb);
}
- const int band2 = blockIdx.z;
- int64_t accum = 0;
-
- // the compiler does not assume that parameters are constant, move them to local variables to give the compiler
- // a hint that those values have to be loaded only once from constant memory.
- int32_t add_shift_cub = ws.add_shift_cub[band2];
- int32_t shift_cub = ws.shift_cub[band2];
- int32_t add_shift_sq = ws.add_shift_sq[band2];
- int32_t shift_sq = ws.shift_sq[band2];
-
- // accumulate per thread
- for (int row = 0;row < rows_per_thread;++row) {
- int32_t accum_thread = accum_thread_reg[row];
- const int32_t x_sq = (int32_t)((((int64_t)accum_thread * accum_thread) + add_shift_sq >> shift_sq));
- accum += (((int64_t)x_sq * accum_thread) + add_shift_cub) >> shift_cub;
- }
-
- // accumulate warp
- accum = warp_reduce(accum);
-
- if (threadIdx.x % 32 == 0)
- {
- accum = (accum + add_shift_inner_accum) >> shift_inner_accum;
- atomicAdd_int64(&accum_global[band2],
- accum);
+ /* Store the per-pixel masked values; adm_cm_reduce_line_kernel cubes and
+ * reduces them with ONE rounded shift_inner_accum shift per row, matching
+ * the CPU adm_cm order. The previous in-kernel tail rounded
+ * per 32x8 warp tile, which drifts the accumulator by O(1) per tile --
+ * invisible in the O(1e12) normal-path accumulators but a 2x error in the
+ * near-zero AIM accumulators at tiny frame sizes. */
+ if (cta_x < buffer_stride) {
+ for (int row = 0;row < rows_per_thread;++row) {
+ const int buffer_row = cta_y + row;
+ if (buffer_row < buffer_h)
+ accum_per_block[(blockIdx.z * buffer_h + buffer_row) *
+ buffer_stride + cta_x] = accum_thread_reg[row];
+ }
}
}
-template <int val_per_thread>
-__device__ __forceinline__ void adm_cm_reduce_line_kernel(int h, int w, int scale, int buffer_h,
- int buffer_stride,
+/* One CTA per (band, line). The CPU (i4_adm_cm) sums the whole row into
+ * accum_inner and applies `(accum_inner + add_shift_inner_accum) >>
+ * shift_inner_accum` ONCE per row; rounding is not distributive over
+ * addition, so applying that shift to warp-partials (the previous scheme)
+ * drifted the int64 accumulator by O(1) per extra warp whenever a row spans
+ * more than one warp. The whole row is therefore reduced
+ * unrounded within a single block and shifted once by thread 0. */
+extern "C" __global__ void adm_cm_reduce_line_kernel(int h, int w, int scale,
+ int buffer_h, int buffer_stride,
const int32_t *buffer,
int64_t *accum) {
const int band = blockIdx.z;
@@ -270,55 +242,46 @@ __device__ __forceinline__ void adm_cm_reduce_line_kernel(int h, int w, int scal
}
int64_t temp_value = 0;
- const int buffer_col = (blockDim.x * blockIdx.x + threadIdx.x) * val_per_thread;
- const int32_t *buffer_loc = buffer + b_off + buffer_col;
- for (int i = 0; i < val_per_thread; ++i) {
- if ((buffer_col + i) < buffer_stride) {
- const int32_t x = buffer_loc[i];
- const int32_t x_sq =
- (int32_t)((((int64_t)x * x) + add_shift_sq) >> shift_sq);
- temp_value += (((int64_t)x_sq * x) + add_shift_cub) >> shift_cub;
- }
+ for (int col = threadIdx.x; col < buffer_stride; col += blockDim.x) {
+ const int32_t x = buffer[b_off + col];
+ const int32_t x_sq =
+ (int32_t)((((int64_t)x * x) + add_shift_sq) >> shift_sq);
+ temp_value += (((int64_t)x_sq * x) + add_shift_cub) >> shift_cub;
}
temp_value = warp_reduce(temp_value);
- if ((threadIdx.x % VMAF_CUDA_THREADS_PER_WARP) == 0) {
+ __shared__ int64_t warp_sums[32]; // upper bound: 1024 threads / warp size
+ const int num_warps = blockDim.x / VMAF_CUDA_THREADS_PER_WARP;
+ if ((threadIdx.x % VMAF_CUDA_THREADS_PER_WARP) == 0)
+ warp_sums[threadIdx.x / VMAF_CUDA_THREADS_PER_WARP] = temp_value;
+ __syncthreads();
+
+ if (threadIdx.x == 0) {
+ int64_t row_accum = 0;
+ for (int i = 0; i < num_warps; ++i)
+ row_accum += warp_sums[i];
const uint32_t shift_inner_accum = __float2uint_ru(__log2f(h));
const uint32_t add_shift_inner_accum = 1 << (shift_inner_accum - 1);
atomicAdd_int64(&accum[band],
- (temp_value + add_shift_inner_accum) >> shift_inner_accum);
+ (row_accum + add_shift_inner_accum) >> shift_inner_accum);
}
}
-#define ADM_CM_REDUCE_LINE(val_per_thread) \
- __global__ void adm_cm_reduce_line_kernel_##val_per_thread ( \
- int h, int w, int scale, int buffer_h, \
- int buffer_stride, const int32_t *buffer, int64_t *accum) \
-{ \
- adm_cm_reduce_line_kernel<val_per_thread>( \
- h, w, scale, buffer_h, buffer_stride, buffer, accum); \
-}
-
#define ADM_CM_LINE(rows_per_thread) \
__global__ void adm_cm_line_kernel_##rows_per_thread ( \
AdmBufferCuda buf, int h, int w, int top, \
int bottom, int left, int right, int start_row, int end_row, int start_col, \
int end_col, int src_stride, int csf_a_stride, int buffer_h, \
- int buffer_stride, int32_t *accum_per_block, AdmFixedParametersCuda params, \
- int scale, int64_t* accum_global, WarpShift ws, \
- const uint32_t shift_inner_accum, const uint32_t add_shift_inner_accum) \
+ int buffer_stride, int32_t *accum_per_block, AdmFixedParametersCuda params) \
{ \
adm_cm_line_kernel<rows_per_thread>( \
buf, h, w, top, bottom, left, right, start_row, end_row, start_col, \
end_col, src_stride, csf_a_stride, buffer_h, buffer_stride, \
- accum_per_block, params,scale, accum_global, \
- ws, shift_inner_accum, add_shift_inner_accum); \
+ accum_per_block, params); \
}
extern "C" {
- // 128 = warps_per_thread * val_per_thread = 32 * 4 -- assuming 32 threads per warp, this might change in the future
- ADM_CM_REDUCE_LINE(4); // adm_cm_reduce_line_kernel_4
ADM_CM_LINE(8); // adm_cm_line_kernel_8
}
diff --git a/libvmaf/src/feature/cuda/integer_adm_cuda.c b/libvmaf/src/feature/cuda/integer_adm_cuda.c
index 51148057..c63dce23 100644
--- a/libvmaf/src/feature/cuda/integer_adm_cuda.c
+++ b/libvmaf/src/feature/cuda/integer_adm_cuda.c
@@ -36,13 +36,6 @@
#define RES_BUFFER_SIZE 4 * 3 * 2
-typedef struct WarpShift {
- uint32_t shift_cub[3];
- uint32_t add_shift_cub[3];
- uint32_t shift_sq[3];
- uint32_t add_shift_sq[3];
-} WarpShift;
-
typedef struct AdmStateCuda {
size_t integer_stride;
AdmBufferCuda buf;
@@ -75,7 +68,7 @@ typedef struct AdmStateCuda {
func_adm_csf_den_scale_line_kernel,
func_adm_csf_den_s123_line_kernel,
// adm_cm kernel
- func_adm_cm_reduce_line_kernel_4,
+ func_adm_cm_reduce_line_kernel,
func_adm_cm_line_kernel_8,
func_i4_adm_cm_line_kernel;
@@ -440,15 +433,17 @@ void i4_adm_cm_device(AdmStateCuda *s, AdmBufferCuda *buf, int w, int h, int src
0, c_stream, args, NULL));
}
{
- const int val_per_thread = 4;
const int warps_per_cta = 4;
const int BLOCKX = VMAF_CUDA_THREADS_PER_WARP * warps_per_cta;
void* args[] = {
&h, &w, &scale, &buffer_h, &buffer_stride,
&buf->tmp_accum->data, &buf->adm_cm[scale]};
- CHECK_CUDA(cu_f, cuLaunchKernel(s->func_adm_cm_reduce_line_kernel_4,
- DIV_ROUND_UP(buffer_stride, BLOCKX * val_per_thread), buffer_h, 3,
+ /* One block per (band, row): the row-wide sum must receive the
+ * rounded shift_inner_accum shift exactly once, as the CPU does,
+ * so the whole row reduces within a single block. */
+ CHECK_CUDA(cu_f, cuLaunchKernel(s->func_adm_cm_reduce_line_kernel,
+ 1, buffer_h, 3,
BLOCKX, 1, 1,
0, c_stream, args, NULL));
}
@@ -472,34 +467,6 @@ void adm_cm_device(AdmStateCuda *s, AdmBufferCuda *buf, int w, int h, int src_st
int buffer_stride = end_col - start_col;
int buffer_h = end_row - start_row;
- // precompute warp shift per band
- //const int32_t shift_sub[3] = {10, 10, 12};
- const int fixed_shift[3] = {4, 4, 3};
-
- // accumulation
- const int32_t shift_xsq[3] = {29, 29, 30};
- const int32_t add_shift_xsq[3] = {268435456, 268435456, 536870912};
-
- const int NUM_BANDS = 3;
- WarpShift ws;
- for (int band = 0;band < NUM_BANDS;++band) {
- ws.shift_cub[band] = (uint32_t)(ceil(log2f(w)));
- if (scale == 0) {
- ws.shift_cub[band] -= fixed_shift[band];
- ws.shift_sq[band] = shift_xsq[band];
- ws.add_shift_sq[band] = add_shift_xsq[band];
- } else {
- ws.shift_sq[band] = 30;
- ws.add_shift_sq[band] = (1 << (ws.shift_sq[band]-1));
- }
- ws.add_shift_cub[band] = 1 << (ws.shift_cub[band] - 1);
- }
-
- // precompute global shift
- uint32_t shift_inner_accum = (uint32_t)(ceil(log2f(h)));
- uint32_t add_shift_inner_accum = 1 << (shift_inner_accum - 1);
-
- // fused
{
const int rows_per_thread = 8;
const int BLOCKX = 32, BLOCKY = 4;
@@ -507,9 +474,7 @@ void adm_cm_device(AdmStateCuda *s, AdmBufferCuda *buf, int w, int h, int src_st
void* args[] = {
&*buf, &h, &w, &top, &bottom, &left, &right, &start_row, &end_row, &start_col,
&end_col, &src_stride, &csf_a_stride, &buffer_h, &buffer_stride,
- &buf->tmp_accum->data, &*p,
- &scale, &buf->adm_cm[scale], &ws,
- &shift_inner_accum, &add_shift_inner_accum
+ &buf->tmp_accum->data, &*p
};
CHECK_CUDA(cu_f, cuLaunchKernel(s->func_adm_cm_line_kernel_8,
@@ -517,6 +482,21 @@ void adm_cm_device(AdmStateCuda *s, AdmBufferCuda *buf, int w, int h, int src_st
BLOCKX, BLOCKY, 1,
0, c_stream, args, NULL));
}
+ {
+ const int warps_per_cta = 4;
+ const int BLOCKX = VMAF_CUDA_THREADS_PER_WARP * warps_per_cta;
+
+ void* args[] = {
+ &h, &w, &scale, &buffer_h, &buffer_stride,
+ &buf->tmp_accum->data, &buf->adm_cm[scale]};
+ /* One block per (band, row), one rounded shift per row -- same
+ * two-pass structure as i4_adm_cm_device; the fused in-kernel
+ * reduction rounded per warp tile instead of per row. */
+ CHECK_CUDA(cu_f, cuLaunchKernel(s->func_adm_cm_reduce_line_kernel,
+ 1, buffer_h, 3,
+ BLOCKX, 1, 1,
+ 0, c_stream, args, NULL));
+ }
}
static void conclude_adm_cm(int64_t *accum, int h,
@@ -1038,7 +1018,7 @@ static int init_fex_cuda(VmafFeatureExtractor *fex, enum VmafPixelFormat pix_fmt
CHECK_CUDA(cu_f, cuModuleGetFunction(&s->func_adm_csf_den_scale_line_kernel, adm_csf_den_module, "adm_csf_den_scale_line_kernel_8_128"));
CHECK_CUDA(cu_f, cuModuleGetFunction(&s->func_adm_csf_den_s123_line_kernel, adm_csf_den_module, "adm_csf_den_s123_line_kernel_8_128"));
- CHECK_CUDA(cu_f, cuModuleGetFunction(&s->func_adm_cm_reduce_line_kernel_4, adm_cm_module, "adm_cm_reduce_line_kernel_4"));
+ CHECK_CUDA(cu_f, cuModuleGetFunction(&s->func_adm_cm_reduce_line_kernel, adm_cm_module, "adm_cm_reduce_line_kernel"));
CHECK_CUDA(cu_f, cuModuleGetFunction(&s->func_adm_cm_line_kernel_8, adm_cm_module, "adm_cm_line_kernel_8"));
CHECK_CUDA(cu_f, cuModuleGetFunction(&s->func_i4_adm_cm_line_kernel, adm_cm_module, "i4_adm_cm_line_kernel"));
While validating a CUDA port of the integer ADM feature we found three independent bit-exactness bugs on current
master(ec3b116), all in upstream code. Tested patches for all three are inlined below (git apply-able against master); happy to open PRs if you prefer.Bug 1 — CUDA
i4_adm_cm_line_kernel: border branches use the wrong mirror and read one row past the csf regionlibvmaf/src/feature/cuda/integer_adm/adm_cm.cuThe
offset_i/offset_jscheme can only express three consecutive-ish rows. Ati == 0(active when a scale'stopclamps to<= 0, i.e. frames smaller than ~150px in a dimension) it reads flt rows{1,2,3}with the center at row 2, while the CPU macros (I4_ADM_CM_THRESH_S_0_*) read{1,0,1}with the center at row 0. Row 3 is also one row past what the csf pass wrote at that scale, so the value picked up is whatever the previous scale's csf left in the buffer — content-dependent stale data. Same for thej == 0columns.Effect: CPU-vs-CUDA
integer_adm2diverges by up to ~1e-1 on tiny frames (ABR-ladder rungs below ~150px hit the same branches at scale 3).Fix: replace the offset scheme with position clamping
pos = min(abs(p), n-1), which reproduces every CPU border pattern exactly — note the CPU is asymmetric: top/left mirror{1,0,1}, bottom/right replicate{n-2,n-1,n-1}— and can never index outside[0,h-1]x[0,w-1]. Interior positions are unchanged.Bug 2 — CUDA adm_cm reductions round per WARP-partial; the CPU rounds once per ROW
adm_cm.cu(adm_cm_reduce_line_kernel_*and the fusedadm_cm_line_kernel_8)The CPU (
adm_cm/i4_adm_cm) accumulatesaccum_innerover an entire row and applies(accum_inner + add_shift_inner_accum) >> shift_inner_accumonce per row. The CUDA kernels apply that same rounded shift to each warp-partial before the atomic add. Rounding is not distributive over addition, so whenever a row spans more than one warp the int64 accumulators drift by O(1) per extra warp — e.g. (traced against a scalar, AVX-disabled CPU run):On the normal path this is O(10) counts in O(1e9..1e12) and only surfaces as an occasional one-ulp flicker of the
%.6fjson output on rounding-boundary frames. But for near-zero accumulators the same grouping error is catastrophic in relative terms: the CPU rounds each small row-sum down to 0 while a 32x8 warp tile spanning 8 rows rounds up to 1 (we measured a clean 2x on a near-zero accumulator at 64x64).Fix:
adm_cm_reduce_line_kernelbecomes one CTA per (band, row) — block-stride loop, shared-memory combine, one rounded shift per row — and the fused scale-0 kernel stores its per-pixel masked values through its previously-unusedaccum_per_blockbuffer and goes through the same reduce kernel (the two-pass structure the i4 path already uses; the reduce kernel already carried the scale-0 shift constants from the original split design). No measurable perf change at 1080p (0.40s vs 0.41s for a 120-frame adm-only run).Bug 3 — x86 SIMD
adm_dwt2*: vector loop computes the mirrored last output column with non-mirrored tapslibvmaf/src/feature/x86/adm_avx2.c,adm_avx512.cThe horizontal-pass vector loops load their four taps at
ind_x[k][j]and assume the indices advance arithmetically across the vector (ind_x[k][j+n] == ind_x[k][j] + 2n).dwt2_src_indices_filtmirrors the taps of the last output column, breaking that assumption. The cutoffequals
half_wwheneverhalf_w ≡ 1 (mod N): the scalar tail becomes empty, the mirrored column is computed inside the vector loop with wrong taps, and the loads run a few elements past the end of thetmplorow into the adjacenttmphibuffer (an OOB-read hazard for sanitizers at exact allocation boundaries).Affected kernels (identical pattern):
adm_dwt2_s123_combined_avx2(N=4),adm_dwt2_16_avx2(N=16),adm_dwt2_s123_combined_avx512(N=8),adm_dwt2_16_avx512(N=64),adm_dwt2_8_avx512(N=32).adm_dwt2_8_avx2is unreachable in the bad alignment (dispatch requires!(w % 8)).The corruption exists at buffer level at any size where some scale's dwt input has
half_w ≡ 1 (mod N), but the wrong values sit in the last output column, which the ADM border trim excludes from every accumulation region at normal sizes — it only reaches the score when a scale is small enough for the region to include it (frames ≲150px wide). Concretely, at 66x34 (scale-1..3 dwt inputs 33/17/9 wide →half_w17/9/5, all ≡ 1 mod 4):10-bit input hits
adm_dwt2_16_avx2the same way (half_w33 ≡ 1 mod 16 at scale 0): 0.952292 vs 0.991931.Fix: cap the cutoff at
half_w - 1so the last column always goes through the (correct) scalar tail:Repro
Validation of the patches
On
master(ec3b116) + both patches, built with-Denable_cuda -Denable_nvcc(CUDA 13.3, sm_120) and compared frame-by-frame at zero tolerance on all shared adm keys (integer_adm2,integer_adm_scale0..3):--cpumask 24): bit-exact at 66x34 8-bit and 66x34 10-bit (both diverge before the patch).Patches
Both apply cleanly on current master with
git apply(save each diff block to a file first).Patch 1/2 — CUDA adm_cm: fix border indexing + round once per row (bugs 1 and 2)
Patch 2/2 — x86 SIMD dwt: leave the mirrored last output column to the scalar tail (bug 3)