GPU PERMANOVA reads an uninitialized pseudo-F: device buffer sized n_perm, kernel launched over n_perm+1

[wasade]

On the accelerated paths, permanova_perm_fp_sW_T allocates and copies back n_perm elements of permutted_sWs, but the host array is n_perm+1 and the kernel is launched over n_perm+1 groupings. The last permutation's pseudo-F is therefore never written on the host side, and the counting loop reads it.

Line numbers below are against main at 372eb9b.

The mismatch

The host array is n_perm+1 (src/distance/permanova.cpp:328, aliased into the sW routine at :292):

TFloat *permutted_fstats = new TFloat[n_perm+1];

The device buffer is n_perm (:135, :146):

skbb_acc_nv::acc_create_buf(permutted_sWs,&permutted_sWs_device,n_perm);
skbb_acc_amd::acc_create_buf(permutted_sWs,&permutted_sWs_device,n_perm);

The chunk loop covers n_perm+1 (:152-153). With a GPU PERM_CHUNK — 200*props.multiProcessorCount, so at least 200 (permanova_dyn_impl.hpp:83, :95) — and the common n_perm = 99, step_perms is 100 and the single chunk passes max_p - tp == 100 as n_grouping_dims. pmn_f_stat_sW_cuda_one launches one block per grouping element (permanova_dyn_impl.hpp:321), takes grouping_el = blockIdx.x (:238) and writes group_sWs[grouping_el] (:309), i.e. indices 0..99 into a 99-element allocation. That is a one-element out-of-bounds device write; cudaMalloc padding is presumably why it has not faulted.

Copy-back is also n_perm (:232, :241):

skbb_acc_nv::acc_copyout_buf(permutted_sWs,permutted_sWs_device,n_perm);

so permutted_fstats[n_perm] keeps whatever new TFloat[n_perm+1] happened to pick up. permanova_all_T transforms all n_perm+1 entries including that one, and permanova_T then counts it (:336-337):

for (uint32_t i=0; i<n_perm; i++) {
  if (permutted_fstats[i+1] >= myfstat) count_larger++;   // i+1 reaches n_perm
}

Impact

Every GPU PERMANOVA p-value is one count of uninitialized heap away from correct — 1/(n_perm+1), so 0.01 at the common n_perm = 99. fstat is index 0 and is unaffected. The CPU paths write all n_perm+1 entries directly and are correct.

Run one compute at a time, the garbage tends to be stable and the result looks reproducible, which is likely why this has gone unnoticed. Run several computes concurrently in one process and the p-value moves between calls at the same seed. That is how we found it, from downstream in unifrac-binaries: four threads x 25 iterations of skbb_permanova_fp64 at a fixed seed, all against the same distance matrix, disagreed with the serial answer on a CUDA build while agreeing exactly on every CPU build.

The same sizing appears on the AMD branch. Under OMPGPU / _OPENACC, where buf_device == buf_host, the kernel would write one element past the mapped region rather than past a cudaMalloc.

Suggested fix

Size both the allocation and the copy-back to match the host array, at all four call sites:

-    skbb_acc_nv::acc_create_buf(permutted_sWs,&permutted_sWs_device,n_perm);
+    skbb_acc_nv::acc_create_buf(permutted_sWs,&permutted_sWs_device,n_perm+1);

-    skbb_acc_nv::acc_copyout_buf(permutted_sWs,permutted_sWs_device,n_perm);
+    skbb_acc_nv::acc_copyout_buf(permutted_sWs,permutted_sWs_device,n_perm+1);

and the two skbb_acc_amd equivalents.

Test gap

src/tests/test_permanova.cpp does not catch this, and neither does unifrac-binaries' CI, because every binary there that runs the skbb PERMANOVA tests ends up on the CPU path. A regression test that runs the same seed twice and requires an identical p-value would catch it, but only in a build where the accelerator is actually engaged.

Downstream workaround in the meantime

biocore/unifrac-binaries#88 keeps fstat asserted in every configuration and gates its p-value equality assertion on skbb_get_acc_mode() == SKBB_ACC_CPU, to be removed once this is fixed.