| Issue |
207883
|
| Summary |
[AMDGPU] Miscompiled post-barrier LDS reload in a large-workgroup HIP kernel on gfx11 (silent wrong results)
|
| Labels |
|
| Assignees |
|
| Reporter |
zpengmei
|
# [AMDGPU] Miscompiled post-barrier LDS reload in a large-workgroup HIP kernel on gfx11 (silent wrong results)
## Summary
On `gfx1100` (RDNA3), a value written to LDS (`__shared__`) by one lane, published
with `__syncthreads()`, and then reloaded by the other lanes of the workgroup can
be read back **inconsistently across lanes** when the workgroup exceeds ~512
threads. The reload appears to be miscompiled: different lanes proceed with
different values of a quantity that must be uniform, producing silently wrong
numerical results (no fault, no assertion).
This was found in rocSOLVER's small-matrix Jacobi eigensolver (`run_syevj`), where
a Jacobi rotation `(c, s)` is computed once by lane `tiy==0`, stored in LDS, and
reloaded by every row-parallel lane after a barrier to apply the rotation. Above
~512 threads the reloaded `(c, s)` differ between lanes, so different matrix rows
receive different rotations and the accumulated eigenvectors lose orthogonality.
## Environment
- GPU: AMD Radeon RX 7900 XT, `gfx1100` (RDNA3, wave32)
- clang / LLVM: `clang version 21.1.8`
- HIP: `7.1.52801`
- Compile: `hipcc -O3 --offload-arch=gfx1100 --gpu-max-threads-per-block=1024`
- Reproduces at both `-O2` and `-O3`.
## Reproducer
Self-contained, ~180 lines, depends only on `hip_runtime` (attached:
`llvm_repro_syevj.cpp`). It runs a faithful transcription of the rocSOLVER kernel
on a batch of ill-conditioned symmetric matrices and checks that the returned
eigenvector matrix `A` is orthogonal (`||A^T A - I|| ~ 1e-13`).
A single `-D` toggles only how the rotation `(c, s)` reaches the applying lanes:
- `BUG=1` (default): compute in lane `tiy==0`, store to `cosines_res/sines_diag`
in LDS, `__syncthreads()`, reload in every lane. ← upstream form
- `BUG=0`: compute `(c, s)` in every lane (a pure function of the same, unmodified
`Acpy` entries, so bit-identical across lanes). No LDS round-trip.
Both are numerically identical; the apply code that follows is byte-for-byte the
same. Only the LDS reload differs.
```
$ hipcc -O3 --offload-arch=gfx1100 --gpu-max-threads-per-block=1024 -DBUG=1 llvm_repro_syevj.cpp -o r_bug && ./r_bug
device=AMD Radeon RX 7900 XT arch=gfx1100 BUG=1
n=39 threads= 400 [small] ||A^T A - I|| worst=6.311e-14 corrupt=0/12 ok
n=48 threads= 576 [small] ||A^T A - I|| worst=4.389e+00 corrupt=12/12 <== MISCOMPILE
n=56 threads= 784 [small] ||A^T A - I|| worst=3.589e+00 corrupt=9/12 <== MISCOMPILE
n=58 threads= 841 [small] ||A^T A - I|| worst=1.047e-13 corrupt=0/12 ok
REPRODUCED (non-orthogonal eigenvectors)
$ hipcc -O3 --offload-arch=gfx1100 --gpu-max-threads-per-block=1024 -DBUG=0 llvm_repro_syevj.cpp -o r_fix && ./r_fix
... all rows: ||A^T A - I|| ~ 1e-13, clean
```
Notes:
- `n=39` → 400 threads (< 512): clean. The corruption only appears once the
workgroup crosses ~512 threads (16 wave32 waves), independent of the matrix size
(the algorithm itself is correct — `BUG=0` is clean at the same sizes).
- The failure is data dependent and, for thread counts that do not pack into whole
wavefronts (e.g. n=58 → 841 threads), non-deterministic run to run. `n=48`
(576 threads = 18 whole waves) reproduces deterministically at 12/12; a given
seed may or may not trip `n=58`. Run a few times / vary the seed to see the
intermittent sizes.
## What was ruled out (rocSOLVER-side investigation)
- Not a memory-model / visibility race: inserting `__threadfence()` before every
barrier is a no-op (bit-identical output).
- Not compiler load/store reordering across the barrier at source level: an
`asm volatile("" ::: "memory")` clobber before the apply has no effect.
- Not a broken `s_barrier`: launching extra no-op lanes past 512 threads (which
hit the barriers but do no writes) stays clean.
- Not register spilling: the kernel does not spill (VGPRs 42, ScratchSize 0).
- Not wave-size selection: rebuilding the kernel `-mwavefrontsize64` does not fix
it (if anything, more failures).
- Not `-O3`-specific: reproduces at `-O2` as well.
- Not inlining scope: moving the apply into a `__noinline__` device function does
not fix it.
- The identical row-parallel global read-modify-write pattern, in isolation
(a small standalone kernel without the surrounding convergence/norm code), is
bit-exact at 841 threads — so the write pattern and the hardware are fine. The
defect only appears in the full kernel context, and is removed by eliminating
the post-barrier LDS reload.
Taken together this points to a codegen defect around the LDS reload of a
barrier-published value under this kernel's control flow / occupancy on gfx11,
rather than an application bug (the two forms are numerically equivalent).
## Ask
- Confirmation of the miscompile and, if possible, identification of the pass /
condition (the ~512-thread threshold and whole-wave-vs-partial-wave determinism
split suggest a scheduling/occupancy interaction).
- Happy to provide LLVM IR / `-save-temps`, an `llvm-reduce`-minimized case, or
ISA on request.
## Downstream
Fixed downstream in rocSOLVER by computing the rotation per-lane instead of
broadcasting it through LDS (equivalent to `BUG=0` above): ROCm/rocm-libraries#9135. That restores correct, deterministic results at no
measurable cost.
<details>
<summary>Reproducer — <code>llvm_repro_syevj.cpp</code></summary>
```cpp
// Self-contained reproducer for an LLVM AMDGPU miscompile on gfx11 (RDNA3).
//
// hipcc -O3 --offload-arch=gfx1100 --gpu-max-threads-per-block=1024 llvm_repro_syevj.cpp -o r && ./r
//
// This is a faithful, dependency-free extraction of rocSOLVER's small-matrix Jacobi
// eigensolver device function (run_syevj / syevj_small_kernel). It diagonalizes a batch
// of symmetric matrices; the accumulated eigenvector matrix A must come out orthogonal
// (A^T A = I). Build with -DBUG=1 (default) to use the upstream form, where the rotation
// (c,s) is produced by the tiy==0 lane, stored in LDS, and reloaded by every row-parallel
// lane after __syncthreads(). Build with -DBUG=0 to compute (c,s) per lane instead.
//
// On gfx1100 the -DBUG=1 build returns non-orthogonal eigenvectors (||A^T A - I|| = O(1),
// silently) once the 2D thread array exceeds ~512 threads (n >= ~40); -DBUG=0 is clean.
// The two forms are numerically identical (same (c,s), a pure function of unmodified Acpy
// entries), so the only difference is the post-barrier LDS reload -> a codegen defect.
#include <hip/hip_runtime.h>
#include <vector>
#include <random>
#include <cmath>
#include <cstdio>
#ifndef BUG
#define BUG 1
#endif
static __device__ void givens(double f, double g, double& c, double& s, double& r){
if(g==0){ c=1; s=0; r=f; } else if(f==0){ c=0; s=1; r=g; }
else { r=std::hypot(f,g); double d=1.0/r; c=std::fabs(f)*d;
s=(f>=0? g*d : -g*d); if(f<0) r=-r; }
}
// Faithful transcription of run_syevj (evect=vector, uplo=upper, esort=none).
__global__ void __launch_bounds__(1024)
run_syevj_k(int n, double* AA, double* AcpyA, double* WW, int lda, double abstol,
int max_sweeps, int forced_ddy)
{
int tid=threadIdx.x, bid=blockIdx.z;
int even_n=n+n%2, half_n=even_n/2;
double* A = AA + (size_t)bid*lda*n;
double* Acpy= AcpyA + (size_t)bid*n*n;
double* W = WW + (size_t)bid*n;
// syevj_get_dims(n, 1024): ddy=min(256,half_n), ddx=min(1024/ddy,half_n); forced_ddy overrides ddy.
int ddy = forced_ddy>0 ? forced_ddy : (256<half_n?256:half_n);
int ddx = (1024/ddy)<half_n ? (1024/ddy) : half_n;
extern __shared__ double lmem[];
double* cosines_res=lmem; double* sines_diag=cosines_res+ddx;
int* top=(int*)(sines_diag+ddx); int* bottom=top+half_n;
int tix=tid/ddy, tiy=tid%ddy;
int dimx=ddx, dimy=ddy;
double c,s,r,f,g,mag,aij,s1,s2,temp1,temp2,local_res=0,local_diag=0;
int i,j,sweeps=0;
if(tiy==0){
for(i=tix;i<n;i+=dimx){
aij=A[i+i*lda]; local_diag+=aij*aij; Acpy[i+i*n]=aij; A[i+i*lda]=1;
for(j=n-1;j>i;j--){ aij=A[i+j*lda]; local_res+=2*aij*aij; Acpy[i+j*n]=aij; Acpy[j+i*n]=aij;
A[i+j*lda]=0; A[j+i*lda]=0; }
}
cosines_res[tix]=local_res; sines_diag[tix]=local_diag;
for(i=tix;i<half_n;i+=dimx){ top[i]=i*2; bottom[i]=i*2+1; }
}
__syncthreads();
local_res=0; local_diag=0;
for(i=0;i<dimx;i++){ local_res+=cosines_res[i]; local_diag+=sines_diag[i]; }
double tolerance=(local_res+local_diag)*abstol*abstol;
double small_num=1e-300;
int count=(half_n-1)/dimx+1;
while(sweeps<max_sweeps && local_res>tolerance){
for(int k=0;k<even_n-1;++k){
for(int cc=0;cc<count;++cc){
int kx=tix+cc*dimx;
i=kx<half_n?top[kx]:n; j=kx<half_n?bottom[kx]:n;
#if BUG
if(tiy==0 && i<n && j<n){
aij=Acpy[i+j*n]; mag=std::fabs(aij);
if(mag*mag<small_num){ c=1; s1=0; }
else { g=2*mag; f=Acpy[j+j*n]-Acpy[i+i*n];
f += (f<0)? -std::hypot(f,g): std::hypot(f,g); givens(f,g,c,s,r); s1=s*aij/mag; }
cosines_res[tix]=c; sines_diag[tix]=s1;
}
__syncthreads();
if(i<n && j<n){ c=cosines_res[tix]; s1=sines_diag[tix]; s2=s1;
#else
if(i<n && j<n){
aij=Acpy[i+j*n]; mag=std::fabs(aij);
if(mag*mag<small_num){ c=1; s1=0; }
else { g=2*mag; f=Acpy[j+j*n]-Acpy[i+i*n];
f += (f<0)? -std::hypot(f,g): std::hypot(f,g); givens(f,g,c,s,r); s1=s*aij/mag; }
}
__syncthreads();
if(i<n && j<n){ s2=s1;
#endif
for(int ky=tiy;ky<half_n;ky+=dimy){ int y1=ky*2,y2=y1+1;
temp1=Acpy[y1+i*n]; temp2=Acpy[y1+j*n];
Acpy[y1+i*n]=c*temp1+s2*temp2; Acpy[y1+j*n]=-s1*temp1+c*temp2;
if(y2<n){ temp1=Acpy[y2+i*n]; temp2=Acpy[y2+j*n];
Acpy[y2+i*n]=c*temp1+s2*temp2; Acpy[y2+j*n]=-s1*temp1+c*temp2; }
temp1=A[y1+i*lda]; temp2=A[y1+j*lda];
A[y1+i*lda]=c*temp1+s2*temp2; A[y1+j*lda]=-s1*temp1+c*temp2;
if(y2<n){ temp1=A[y2+i*lda]; temp2=A[y2+j*lda];
A[y2+i*lda]=c*temp1+s2*temp2; A[y2+j*lda]=-s1*temp1+c*temp2; }
}
}
__syncthreads();
if(i<n && j<n){
for(int ky=tiy;ky<half_n;ky+=dimy){ int y1=ky*2,y2=y1+1;
temp1=Acpy[i+y1*n]; temp2=Acpy[j+y1*n];
Acpy[i+y1*n]=c*temp1+s1*temp2; Acpy[j+y1*n]=-s2*temp1+c*temp2;
if(y2<n){ temp1=Acpy[i+y2*n]; temp2=Acpy[j+y2*n];
Acpy[i+y2*n]=c*temp1+s1*temp2; Acpy[j+y2*n]=-s2*temp1+c*temp2; }
}
}
__syncthreads();
if(tiy==0 && i<n && j<n){ Acpy[i+j*n]=0; Acpy[j+i*n]=0; }
__syncthreads();
if(tiy==0 && kx<half_n){
if(i>0){ if(i==2||i==even_n-1) top[kx]=i-1; else top[kx]=i+((i%2==0)?-2:2); }
if(j==2||j==even_n-1) bottom[kx]=j-1; else bottom[kx]=j+((j%2==0)?-2:2);
}
__syncthreads();
}
}
if(tiy==0){ local_res=0; for(i=tix;i<n;i+=dimx) for(j=0;j<i;j++) local_res+=2*Acpy[i+j*n]*Acpy[i+j*n];
cosines_res[tix]=local_res; }
__syncthreads();
local_res=0; for(i=0;i<dimx;i++) local_res+=cosines_res[i];
sweeps++;
}
if(tiy==0) for(i=tix;i<n;i+=dimx) W[i]=Acpy[i+i*n];
}
int main(){
hipDeviceProp_t p{}; hipGetDeviceProperties(&p,0);
printf("device=%s arch=%s BUG=%d\n\n",p.name,p.gcnArchName,BUG);
std::mt19937_64 rng(7);
int worst_fail=0;
for(int n:{39,48,56,58}){
int even_n=n+n%2, half_n=even_n/2, B=12;
std::normal_distribution<double> N(0,1);
// near-null SPD batch: A = Q diag(lam) Q^T, lam={1..1e-12, then 1e-16}
std::vector<double> A((size_t)B*n*n), q(n*n), lam(n);
for(int k=0;k<n;k++) lam[k]= k<8? std::pow(1e-12,double(k)/7):1e-16;
for(int b=0;b<B;b++){ for(double&x:q)x=N(rng);
for(int k=0;k<n;k++){ for(int pp=0;pp<k;pp++){ double d=0; for(int t=0;t<n;t++)d+=q[t+k*n]*q[t+pp*n];
for(int t=0;t<n;t++)q[t+k*n]-=d*q[t+pp*n]; }
double nr=0; for(int t=0;t<n;t++)nr+=q[t+k*n]*q[t+k*n]; nr=1/std::sqrt(nr);
for(int t=0;t<n;t++)q[t+k*n]*=nr; }
double*Ab=A.data()+(size_t)b*n*n;
for(int jj=0;jj<n;jj++)for(int ii=0;ii<n;ii++){ double a=0; for(int k=0;k<n;k++)a+=q[ii+k*n]*lam[k]*q[jj+k*n]; Ab[ii+jj*n]=a; }
}
double *dA,*dAcpy,*dW;
hipMalloc(&dA,sizeof(double)*A.size()); hipMalloc(&dAcpy,sizeof(double)*A.size()); hipMalloc(&dW,sizeof(double)*B*n);
hipMemcpy(dA,A.data(),sizeof(double)*A.size(),hipMemcpyHostToDevice);
int ddy = 256<half_n?256:half_n; int ddx=(1024/ddy)<half_n?(1024/ddy):half_n;
int threads=ddx*ddy, sh=(2*ddx)*sizeof(double)+(2*half_n)*sizeof(int);
run_syevj_k<<<dim3(1,1,B),dim3(threads),sh>>>(n,dA,dAcpy,dW,n,1e-12,100,0);
hipDeviceSynchronize();
std::vector<double> V(A.size()); hipMemcpy(V.data(),dA,sizeof(double)*A.size(),hipMemcpyDeviceToHost);
double worst=0; int nfail=0;
for(int b=0;b<B;b++){ const double*Vb=V.data()+(size_t)b*n*n; double o=0;
for(int a=0;a<n;a++)for(int cc=0;cc<n;cc++){ double d=0; for(int k=0;k<n;k++)d+=Vb[k+a*n]*Vb[k+cc*n]; d-=(a==cc); o+=d*d; }
o=std::sqrt(o); if(o>worst)worst=o; if(o>1e-8)nfail++; }
printf("n=%2d threads=%4d [%s] ||A^T A - I|| worst=%.3e corrupt=%d/%d %s\n",
n,threads, n<=58?"small":"blk", worst,nfail,B, nfail?" <== MISCOMPILE":" ok");
worst_fail+=nfail;
hipFree(dA);hipFree(dAcpy);hipFree(dW);
}
printf("\n%s\n", worst_fail? "REPRODUCED (non-orthogonal eigenvectors)":"clean");
return worst_fail?1:0;
}
```
</details>
![]()
_______________________________________________
llvm-bugs mailing list
[email protected]
https://lists.llvm.org/cgi-bin/mailman/listinfo/llvm-bugs
