diff --git a/examples/hmctest/adjointtest4.cpp b/examples/hmctest/adjointtest4.cpp
index bc738efb..c28ee171 100644
--- a/examples/hmctest/adjointtest4.cpp
+++ b/examples/hmctest/adjointtest4.cpp
@@ -771,16 +771,10 @@ static int runTest17(int gpuDevice) {
             { 2, 7, 2, 0 },
             { 3, 6, 3, 0 }
         };
-        /* GPU writes kPaddedStateCount² values with stride kPaddedStateCount;
-         * CPU writes kStateCount² values with stride kStateCount.
-         * Allocate the padded size to prevent overflow, then re-index to S×S. */
-        int Sp = useGpu ? ((S <= 4) ? 4 : (S <= 16) ? 16 : (S <= 32) ? 32 : S) : S;
-        std::vector<double> gradBuf(Sp * Sp, 0.0);
-        beagleCalculateAdjointDerivative(inst, branchOps, 0, 0, 4, 0, 4, gradBuf.data(), NULL);
+        /* Both CPU and GPU write kStateCount² values with stride kStateCount
+         * (GPU downsamples from its internal kPaddedStateCount² buffer). */
         std::vector<double> grad(S * S, 0.0);
-        for (int ls = 0; ls < S; ls++)
-            for (int rs = 0; rs < S; rs++)
-                grad[ls*S+rs] = gradBuf[ls*Sp+rs];
+        beagleCalculateAdjointDerivative(inst, branchOps, 0, 0, 4, 0, 4, grad.data(), NULL);
 
         const char* preLabels[] = { "pre6", "pre7", "pre8", "pre9" };
         for (int b = 0; b < 4; b++) {
@@ -794,12 +788,12 @@ static int runTest17(int gpuDevice) {
         const char* bNames[4] = { "chimp", "human", "gorilla", "internal" };
         fprintf(stdout, "  Per-branch gradient G[ls,rs] (first 4x4 sub-block):\n");
         for (int b = 0; b < 4; b++) {
-            std::vector<double> gbuf(Sp * Sp, 0.0);
+            std::vector<double> gbuf(S * S, 0.0);
             beagleCalculateAdjointDerivative(inst, &branchOps[b], 0, 0, 4, 0, 1, gbuf.data(), NULL);
             fprintf(stdout, "    Branch %s (transMatrix=%d):\n", bNames[b], branchOps[b].branchTransitionMatrix);
             for (int ls = 0; ls < 4; ls++) {
                 for (int rs = 0; rs < 4; rs++)
-                    fprintf(stdout, " %12.4f", gbuf[ls*Sp+rs]);
+                    fprintf(stdout, " %12.4f", gbuf[ls*S+rs]);
                 fprintf(stdout, "\n");
             }
         }
diff --git a/libhmsbeagle/CPU/BeagleCPUPlugin.cpp b/libhmsbeagle/CPU/BeagleCPUPlugin.cpp
index 7b69eecf..44ad5351 100644
--- a/libhmsbeagle/CPU/BeagleCPUPlugin.cpp
+++ b/libhmsbeagle/CPU/BeagleCPUPlugin.cpp
@@ -50,7 +50,7 @@ Plugin("CPU", "CPU")
 
 
 extern "C" {
-void* plugin_init(void){
+__attribute__((visibility("default"))) void* plugin_init(void){
 	return new beagle::cpu::BeagleCPUPlugin();
 }
 }
diff --git a/libhmsbeagle/CPU/BeagleCPUSSEPlugin.cpp b/libhmsbeagle/CPU/BeagleCPUSSEPlugin.cpp
index c8d5aaab..32eceb58 100644
--- a/libhmsbeagle/CPU/BeagleCPUSSEPlugin.cpp
+++ b/libhmsbeagle/CPU/BeagleCPUSSEPlugin.cpp
@@ -176,7 +176,7 @@ bool check_sse2(){
 #endif
 
 
-void* plugin_init(void){
+__attribute__((visibility("default"))) void* plugin_init(void){
 	if(!check_sse2()){
 		return NULL;	// no SSE no plugin?!
 	}
diff --git a/libhmsbeagle/CPU/BeagleCPUSpectralPlugin.cpp b/libhmsbeagle/CPU/BeagleCPUSpectralPlugin.cpp
index dc4a1180..e2203870 100644
--- a/libhmsbeagle/CPU/BeagleCPUSpectralPlugin.cpp
+++ b/libhmsbeagle/CPU/BeagleCPUSpectralPlugin.cpp
@@ -48,7 +48,7 @@ Plugin("CPU-Spectral", "CPU-Sprectral")
 
 
 extern "C" {
-void* plugin_init(void){
+__attribute__((visibility("default"))) void* plugin_init(void){
 	return new beagle::cpu::BeagleCPUSpectralPlugin();
 }
 }
diff --git a/libhmsbeagle/CPU/CMakeLists.txt b/libhmsbeagle/CPU/CMakeLists.txt
index 657b6165..9512cccf 100644
--- a/libhmsbeagle/CPU/CMakeLists.txt
+++ b/libhmsbeagle/CPU/CMakeLists.txt
@@ -25,6 +25,8 @@ SET_TARGET_PROPERTIES(hmsbeagle-cpu
 	SUFFIX "${BEAGLE_PLUGIN_SUFFIX}"
     )
 
+target_compile_options(hmsbeagle-cpu PRIVATE -fvisibility=hidden -fvisibility-inlines-hidden)
+
 if(BUILD_SSE)
 add_library(hmsbeagle-cpu-sse SHARED
         BeagleCPU4StateSSEImpl.h
@@ -56,6 +58,8 @@ SET_TARGET_PROPERTIES(hmsbeagle-cpu-sse
     SOVERSION "${BEAGLE_PLUGIN_VERSION_EXTENDED}"
 	SUFFIX "${BEAGLE_PLUGIN_SUFFIX}"
     )
+
+target_compile_options(hmsbeagle-cpu-sse PRIVATE -fvisibility=hidden -fvisibility-inlines-hidden)
 endif(BUILD_SSE)
 
 if(BUILD_SPECTRAL)
@@ -93,4 +97,6 @@ SET_TARGET_PROPERTIES(hmsbeagle-cpu-spectral
     SOVERSION "${BEAGLE_PLUGIN_VERSION_EXTENDED}"
 	SUFFIX "${BEAGLE_PLUGIN_SUFFIX}"
     )
+
+target_compile_options(hmsbeagle-cpu-spectral PRIVATE -fvisibility=hidden -fvisibility-inlines-hidden)
 endif(BUILD_SPECTRAL)
diff --git a/libhmsbeagle/CPU/EigenDecompositionSpectral.hpp b/libhmsbeagle/CPU/EigenDecompositionSpectral.hpp
index b3ad816d..058668d6 100644
--- a/libhmsbeagle/CPU/EigenDecompositionSpectral.hpp
+++ b/libhmsbeagle/CPU/EigenDecompositionSpectral.hpp
@@ -87,6 +87,13 @@ void EigenDecompositionSpectral<BEAGLE_CPU_EIGEN_GENERIC>::setEigenDecomposition
     memcpy(tIvec.data(), evec.data(), kMatrixStride * kStateCount * sizeof(REALTYPE));
     transposeInPlace(tIvec.data());
 
+    if (getenv("BEAGLE_DEBUG_EIGEN")) {
+        fprintf(stderr, "[CPU setEigenDecomposition] eigenIndex=%d kStateCount=%d kMatrixStride=%d\n", eigenIndex, kStateCount, kMatrixStride);
+        fprintf(stderr, "[CPU] eval: "); for (int i=0;i<kEigenValuesSize;i++) fprintf(stderr, "%.6f ", (double)eval[i]); fprintf(stderr, "\n");
+        fprintf(stderr, "[CPU] evec:\n"); for (int i=0;i<kStateCount;i++){ for(int j=0;j<kStateCount;j++) fprintf(stderr, "%.6f ", (double)evec[i*kMatrixStride+j]); fprintf(stderr, "\n"); }
+        fprintf(stderr, "[CPU] ivec:\n"); for (int i=0;i<kStateCount;i++){ for(int j=0;j<kStateCount;j++) fprintf(stderr, "%.6f ", (double)ivec[i*kMatrixStride+j]); fprintf(stderr, "\n"); }
+        fflush(stderr);
+    }
     adjointMethodsStorage[eigenIndex] = std::make_shared<AdjointIntegralPlan<REALTYPE>>(
         eigenValuesStorage[eigenIndex].data(), kStateCount, !isComplex);
 }
diff --git a/libhmsbeagle/GPU/BeagleGPUImpl.hpp b/libhmsbeagle/GPU/BeagleGPUImpl.hpp
index 446a8b43..234a8fa6 100644
--- a/libhmsbeagle/GPU/BeagleGPUImpl.hpp
+++ b/libhmsbeagle/GPU/BeagleGPUImpl.hpp
@@ -462,6 +462,14 @@ int BeagleGPUImpl<BEAGLE_GPU_GENERIC>::createInstance(int tipCount,
         kFlags |= BEAGLE_FLAG_EIGEN_REAL;
     }
 
+    // kFlags is rebuilt from scratch above (kFlags = 0) and only specific bits are copied
+    // over from preferenceFlags/requirementFlags; BEAGLE_FLAG_SPECTRAL_REPRESENTATION must be
+    // copied too, since BeagleGPUSpectralImpl and the checks just below key off this bit in
+    // kFlags (not off the caller-supplied flags directly).
+    if (preferenceFlags & BEAGLE_FLAG_SPECTRAL_REPRESENTATION || requirementFlags & BEAGLE_FLAG_SPECTRAL_REPRESENTATION) {
+        kFlags |= BEAGLE_FLAG_SPECTRAL_REPRESENTATION;
+    }
+
     if (requirementFlags & BEAGLE_FLAG_INVEVEC_TRANSPOSED || preferenceFlags & BEAGLE_FLAG_INVEVEC_TRANSPOSED)
         kFlags |= BEAGLE_FLAG_INVEVEC_TRANSPOSED;
     else
@@ -497,7 +505,12 @@ int BeagleGPUImpl<BEAGLE_GPU_GENERIC>::createInstance(int tipCount,
     kPartialsSize = kPaddedPatternCount * kPaddedStateCount * kCategoryCount;
     kMatrixSize = kPaddedStateCount * kPaddedStateCount;
 
-    if (kFlags & BEAGLE_FLAG_EIGEN_COMPLEX)
+    // Spectral kernels (kernelsSpectralIfDef*.cu / kernelsSpectral*.cu) always read a
+    // real+imaginary pair per eigenstate (SPECTRAL_EIGENVALS_GPU unconditionally indexes
+    // eigenValues[PADDED_STATE_COUNT + state]) regardless of BEAGLE_FLAG_EIGEN_COMPLEX, so the
+    // device buffer must be sized/copied as complex whenever spectral representation is in use,
+    // otherwise that read runs past the data actually written by setEigenDecomposition.
+    if ((kFlags & BEAGLE_FLAG_EIGEN_COMPLEX) || (kFlags & BEAGLE_FLAG_SPECTRAL_REPRESENTATION))
         kEigenValuesSize = 2 * kPaddedStateCount;
     else
         kEigenValuesSize = kPaddedStateCount;
@@ -4743,6 +4756,7 @@ void BeagleGPUImpl<BEAGLE_GPU_GENERIC>::dispatchPruneSS(GPUPtr s1, GPUPtr s2, GP
                                                          GPUPtr scalingFactors, GPUPtr cumulativeScaling,
                                                          unsigned int startPattern, unsigned int endPattern,
                                                          int rescale, int streamIndex, int waitIndex) {
+    fprintf(stderr, "[DISPATCH] BASE BeagleGPUImpl::dispatchPruneSS called\n"); fflush(stderr);
     kernels->StatesStatesPruningDynamicScaling(s1, s2, p3,
         dMatrices[c1MatIdx], dMatrices[c2MatIdx],
         scalingFactors, cumulativeScaling,
diff --git a/libhmsbeagle/GPU/BeagleGPUSpectralImpl.h b/libhmsbeagle/GPU/BeagleGPUSpectralImpl.h
index 6ff2f311..56c8baa9 100644
--- a/libhmsbeagle/GPU/BeagleGPUSpectralImpl.h
+++ b/libhmsbeagle/GPU/BeagleGPUSpectralImpl.h
@@ -38,6 +38,12 @@ class BeagleGPUSpectralImpl : public BeagleGPUImpl<Real> {
     GPUPtr  dSpectralDistancesOrigin;
     GPUPtr* dSpectralDistances;
     int*    hEigenIndexForMatrix;
+    /* Per-matrix element stride within dSpectralDistancesOrigin, i.e.
+     * AlignMemOffset(kCategoryCount * sizeof(Real)) / sizeof(Real). Device
+     * alignment padding can make this larger than kCategoryCount, so any code
+     * that indexes into the pooled origin buffer by matrix index (e.g. the
+     * adjoint offset-queue) must multiply by this, not by kCategoryCount. */
+    unsigned int kSpectralDistanceStrideElements;
 
     /* Backward eigenvector buffers for the parent branch in pre-order.
      * dEvecT[ei]  = U   stored row-major (dEvecT[j*S+k] = U[j,k]).
diff --git a/libhmsbeagle/GPU/BeagleGPUSpectralImpl.hpp b/libhmsbeagle/GPU/BeagleGPUSpectralImpl.hpp
index b818ceca..10e7a418 100644
--- a/libhmsbeagle/GPU/BeagleGPUSpectralImpl.hpp
+++ b/libhmsbeagle/GPU/BeagleGPUSpectralImpl.hpp
@@ -88,7 +88,8 @@ int BeagleGPUSpectralImpl<BEAGLE_GPU_GENERIC>::createInstance(
     GPUInterface* gpuIf = this->gpu;
 
     dSpectralDistances = (GPUPtr*) malloc(sizeof(GPUPtr) * this->kMatrixCount);
-    size_t distStride = this->kCategoryCount * sizeof(Real);
+    size_t distStride = gpuIf->AlignMemOffset(this->kCategoryCount * sizeof(Real));
+    kSpectralDistanceStrideElements = (unsigned int)(distStride / sizeof(Real));
     dSpectralDistancesOrigin = gpuIf->AllocateMemory(this->kMatrixCount * distStride);
     for (int i = 0; i < this->kMatrixCount; i++) {
         dSpectralDistances[i] = gpuIf->CreateSubPointer(dSpectralDistancesOrigin, distStride * i, distStride);
@@ -96,7 +97,7 @@ int BeagleGPUSpectralImpl<BEAGLE_GPU_GENERIC>::createInstance(
 
     // Backward eigenvector arrays: one S*S block per eigen decomposition.
     int S = this->kPaddedStateCount;
-    size_t matStride = (size_t)S * S * sizeof(Real);
+    size_t matStride = gpuIf->AlignMemOffset((size_t)S * S * sizeof(Real));
     dEvecT  = (GPUPtr*) malloc(sizeof(GPUPtr) * eigenDecompositionCount);
     dIevcT  = (GPUPtr*) malloc(sizeof(GPUPtr) * eigenDecompositionCount);
     dEvecTOrigin  = gpuIf->AllocateMemory(eigenDecompositionCount * matStride);
@@ -191,6 +192,7 @@ void BeagleGPUSpectralImpl<BEAGLE_GPU_GENERIC>::dispatchPruneSS(
         GPUPtr scalingFactors, GPUPtr cumulativeScaling,
         unsigned int startPattern, unsigned int endPattern,
         int rescale, int streamIndex, int waitIndex) {
+    fprintf(stderr, "[DISPATCH] SPECTRAL BeagleGPUSpectralImpl::dispatchPruneSS called\n"); fflush(stderr);
     int ei1 = hEigenIndexForMatrix[c1MatIdx];
     int ei2 = hEigenIndexForMatrix[c2MatIdx];
     this->kernels->StatesStatesPruningSpectral(s1, s2, p3,
@@ -215,7 +217,12 @@ int BeagleGPUSpectralImpl<BEAGLE_GPU_GENERIC>::setEigenDecomposition(
     const int SC = this->kStateCount;
     const int SS = S * S;
     // kEigenValuesSize is private to BeagleGPUImpl; recompute the same way it does.
-    const int eigenValuesSize = (this->kFlags & BEAGLE_FLAG_EIGEN_COMPLEX) ? 2 * S : S;
+    // Spectral kernels always read a real+imaginary pair per eigenstate regardless of
+    // BEAGLE_FLAG_EIGEN_COMPLEX (see matching comment/fix in BeagleGPUImpl::createInstance),
+    // so this local copy must widen for BEAGLE_FLAG_SPECTRAL_REPRESENTATION too, or the
+    // imaginary half of the device buffer is left uninitialized.
+    const int eigenValuesSize = ((this->kFlags & BEAGLE_FLAG_EIGEN_COMPLEX) ||
+                                  (this->kFlags & BEAGLE_FLAG_SPECTRAL_REPRESENTATION)) ? 2 * S : S;
 
     // Forward:  dEvec[row*S+col]  = U[col,row];    dIevc[row*S+col]  = U^-1[col,row]
     // Backward: dEvecT[row*S+col] = U[row,col];    dIevcT[row*S+col] = U^-1[row,col]
@@ -280,6 +287,13 @@ int BeagleGPUSpectralImpl<BEAGLE_GPU_GENERIC>::setEigenDecomposition(
             Eval[S + i] = (Real) inEigenValues[SC + i];
     }
 
+    if (getenv("BEAGLE_DEBUG_EIGEN")) {
+        fprintf(stderr, "[GPU setEigenDecomposition] eigenIndex=%d S=%d SC=%d\n", eigenIndex, S, SC);
+        fprintf(stderr, "[GPU] Eval: "); for (int i=0;i<eigenValuesSize;i++) fprintf(stderr, "%.6f ", (double)Eval[i]); fprintf(stderr, "\n");
+        fprintf(stderr, "[GPU] Evec:\n"); for (int i=0;i<S;i++){ for(int j=0;j<S;j++) fprintf(stderr, "%.6f ", (double)Evec[i*S+j]); fprintf(stderr, "\n"); }
+        fprintf(stderr, "[GPU] Ievc:\n"); for (int i=0;i<S;i++){ for(int j=0;j<S;j++) fprintf(stderr, "%.6f ", (double)Ievc[i*S+j]); fprintf(stderr, "\n"); }
+        fflush(stderr);
+    }
     GPUInterface* gpuIf = this->gpu;
     gpuIf->MemcpyHostToDevice(this->dEvec[eigenIndex],        Evec,  sizeof(Real) * SS);
     gpuIf->MemcpyHostToDevice(this->dIevc[eigenIndex],        Ievc,  sizeof(Real) * SS);
@@ -599,7 +613,7 @@ int BeagleGPUSpectralImpl<BEAGLE_GPU_GENERIC>::calculateAdjointCrossProducts(
             rec[3] = (unsigned int)ei * SS;
             rec[4] = (unsigned int)ei * this->getEvecStrideElements();
             rec[5] = (unsigned int)ei * this->getEvalStrideElements();
-            rec[6] = (unsigned int)matIdx * this->kCategoryCount;
+            rec[6] = (unsigned int)matIdx * kSpectralDistanceStrideElements;
             rec[7] = (unsigned int)ei * SS;
             rec[8] = isAllReal ? 1u : 0u;
         }
@@ -625,13 +639,24 @@ int BeagleGPUSpectralImpl<BEAGLE_GPU_GENERIC>::calculateAdjointCrossProducts(
     }
 
     /* Download gradient to host.  The first eigen decomp's gradient is the
-     * primary output.  Convert Real→double. */
+     * primary output.  Convert Real→double.
+     *
+     * The device buffer is laid out S×S with S=kPaddedStateCount (padding
+     * rows/cols beyond kStateCount hold unused/garbage values from the
+     * padding eigenstates), but callers (matching BeagleCPUImpl's
+     * convention) allocate outSumDerivatives sized kStateCount*kStateCount.
+     * Copying the full padded SS run would overflow that buffer whenever
+     * kPaddedStateCount > kStateCount — only copy the top-left
+     * kStateCount×kStateCount submatrix, row by row with the correct
+     * strides on each side. */
     if (count > 0 && outSumDerivatives != NULL) {
         const int ei0 = hEigenIndexForMatrix[eigenIndices[0]];
+        const int SC = this->kStateCount;
         Real* hGrad = (Real*) gpuIf->CallocHost(sizeof(Real), SS);
         gpuIf->MemcpyDeviceToHost(hGrad, dGradient[ei0], SS * sizeof(Real));
-        for (int k = 0; k < SS; k++)
-            outSumDerivatives[k] = (double)hGrad[k];
+        for (int i = 0; i < SC; i++)
+            for (int j = 0; j < SC; j++)
+                outSumDerivatives[i * SC + j] = (double)hGrad[i * S + j];
         gpuIf->FreeHostMemory(hGrad);
     }
 
diff --git a/libhmsbeagle/GPU/CMake_CUDA/CMakeLists.txt b/libhmsbeagle/GPU/CMake_CUDA/CMakeLists.txt
index dee041db..6d5308c0 100644
--- a/libhmsbeagle/GPU/CMake_CUDA/CMakeLists.txt
+++ b/libhmsbeagle/GPU/CMake_CUDA/CMakeLists.txt
@@ -62,3 +62,5 @@ SET_TARGET_PROPERTIES(hmsbeagle-cuda
     SOVERSION "${BEAGLE_PLUGIN_VERSION_EXTENDED}"
 	SUFFIX "${BEAGLE_PLUGIN_SUFFIX}"
     )
+
+target_compile_options(hmsbeagle-cuda PRIVATE -fvisibility=hidden -fvisibility-inlines-hidden)
diff --git a/libhmsbeagle/GPU/CMake_CUDASpectral/CMakeLists.txt b/libhmsbeagle/GPU/CMake_CUDASpectral/CMakeLists.txt
index 4fb5d308..b43649d4 100644
--- a/libhmsbeagle/GPU/CMake_CUDASpectral/CMakeLists.txt
+++ b/libhmsbeagle/GPU/CMake_CUDASpectral/CMakeLists.txt
@@ -70,3 +70,5 @@ SET_TARGET_PROPERTIES(hmsbeagle-cuda-spectral
     SOVERSION "${BEAGLE_PLUGIN_VERSION_EXTENDED}"
     SUFFIX "${BEAGLE_PLUGIN_SUFFIX}"
     )
+
+target_compile_options(hmsbeagle-cuda-spectral PRIVATE -fvisibility=hidden -fvisibility-inlines-hidden)
diff --git a/libhmsbeagle/GPU/CMake_OpenCL/CMakeLists.txt b/libhmsbeagle/GPU/CMake_OpenCL/CMakeLists.txt
index 129a1492..c61007ff 100644
--- a/libhmsbeagle/GPU/CMake_OpenCL/CMakeLists.txt
+++ b/libhmsbeagle/GPU/CMake_OpenCL/CMakeLists.txt
@@ -53,3 +53,5 @@ SET_TARGET_PROPERTIES(hmsbeagle-opencl
     SOVERSION "${BEAGLE_PLUGIN_VERSION_EXTENDED}"
 	SUFFIX "${BEAGLE_PLUGIN_SUFFIX}"
     )
+
+target_compile_options(hmsbeagle-opencl PRIVATE -fvisibility=hidden -fvisibility-inlines-hidden)
diff --git a/libhmsbeagle/GPU/CMake_OpenCLSpectral/CMakeLists.txt b/libhmsbeagle/GPU/CMake_OpenCLSpectral/CMakeLists.txt
index 87f68ec0..010930fc 100644
--- a/libhmsbeagle/GPU/CMake_OpenCLSpectral/CMakeLists.txt
+++ b/libhmsbeagle/GPU/CMake_OpenCLSpectral/CMakeLists.txt
@@ -62,3 +62,5 @@ SET_TARGET_PROPERTIES(hmsbeagle-opencl-spectral
     SOVERSION "${BEAGLE_PLUGIN_VERSION_EXTENDED}"
     SUFFIX "${BEAGLE_PLUGIN_SUFFIX}"
     )
+
+target_compile_options(hmsbeagle-opencl-spectral PRIVATE -fvisibility=hidden -fvisibility-inlines-hidden)
diff --git a/libhmsbeagle/GPU/CUDAPlugin.cpp b/libhmsbeagle/GPU/CUDAPlugin.cpp
index 03ec13c8..3481613c 100644
--- a/libhmsbeagle/GPU/CUDAPlugin.cpp
+++ b/libhmsbeagle/GPU/CUDAPlugin.cpp
@@ -77,7 +77,7 @@ CUDAPlugin::~CUDAPlugin()
 
 
 extern "C" {
-void* plugin_init(void){
+__attribute__((visibility("default"))) void* plugin_init(void){
 	return new beagle::gpu::CUDAPlugin();
 }
 }
diff --git a/libhmsbeagle/GPU/CUDASpectralPlugin.cpp b/libhmsbeagle/GPU/CUDASpectralPlugin.cpp
index b0c36e10..92f2bd61 100644
--- a/libhmsbeagle/GPU/CUDASpectralPlugin.cpp
+++ b/libhmsbeagle/GPU/CUDASpectralPlugin.cpp
@@ -74,7 +74,7 @@ CUDASpectralPlugin::~CUDASpectralPlugin()
 
 
 extern "C" {
-void* plugin_init(void) {
+__attribute__((visibility("default"))) void* plugin_init(void) {
     return new beagle::gpu::CUDASpectralPlugin();
 }
 }
diff --git a/libhmsbeagle/GPU/GPUInterfaceOpenCL.cpp b/libhmsbeagle/GPU/GPUInterfaceOpenCL.cpp
index 62a4fd34..f26f11e0 100644
--- a/libhmsbeagle/GPU/GPUInterfaceOpenCL.cpp
+++ b/libhmsbeagle/GPU/GPUInterfaceOpenCL.cpp
@@ -544,15 +544,30 @@ void GPUInterface::LaunchKernel(GPUFunction deviceFunction,
     fprintf(stderr,"\t\t\tEntering GPUInterface::LaunchKernel\n");
 #endif
 
+    char kernelNameBuf[256] = {0};
+    clGetKernelInfo(deviceFunction, CL_KERNEL_FUNCTION_NAME, sizeof(kernelNameBuf), kernelNameBuf, NULL);
+    if (getenv("BEAGLE_DEBUG_KERNEL_ARGS")) {
+        fprintf(stderr, "[LaunchKernel] %s: parameterCountV=%d totalParameterCount=%d\n",
+                kernelNameBuf, parameterCountV, totalParameterCount);
+    }
+
     va_list parameters;
     va_start(parameters, totalParameterCount);
     for(int i = 0; i < parameterCountV; i++) {
         void* param = (void*)(size_t)va_arg(parameters, GPUPtr);
+        if (getenv("BEAGLE_DEBUG_KERNEL_ARGS")) {
+            fprintf(stderr, "[LaunchKernel] %s: ptr arg[%d] = %p\n", kernelNameBuf, i, param);
+            fflush(stderr);
+        }
 
         SAFE_CL(clSetKernelArg(deviceFunction, i, sizeof(param), &param));
     }
     for(int i = parameterCountV; i < totalParameterCount; i++) {
         unsigned int param = va_arg(parameters, unsigned int);
+        if (getenv("BEAGLE_DEBUG_KERNEL_ARGS")) {
+            fprintf(stderr, "[LaunchKernel] %s: int arg[%d] = %u (signed: %d)\n", kernelNameBuf, i, param, (int)param);
+            fflush(stderr);
+        }
 
         SAFE_CL(clSetKernelArg(deviceFunction, i, sizeof(unsigned int), &param));
     }
@@ -569,6 +584,14 @@ void GPUInterface::LaunchKernel(GPUFunction deviceFunction,
     globalWorkSize[1] = block.y * grid.y;
     globalWorkSize[2] = block.z * grid.z;
 
+    if (getenv("BEAGLE_DEBUG_KERNEL_ARGS")) {
+        fprintf(stderr, "[LaunchKernel] %s: block=(%d,%d,%d) grid=(%d,%d,%d) local=(%zu,%zu,%zu) global=(%zu,%zu,%zu)\n",
+                kernelNameBuf, block.x, block.y, block.z, grid.x, grid.y, grid.z,
+                localWorkSize[0], localWorkSize[1], localWorkSize[2],
+                globalWorkSize[0], globalWorkSize[1], globalWorkSize[2]);
+        fflush(stderr);
+    }
+
 #ifdef BEAGLE_DEBUG_VALUES
     for (int i=0; i<3; i++) {
         printf("localWorkSize[%d]  = %lu\n", i, localWorkSize[i]);
@@ -579,6 +602,23 @@ void GPUInterface::LaunchKernel(GPUFunction deviceFunction,
     printf("local = %lu\n\n", local);
 #endif
 
+    if (getenv("BEAGLE_DEBUG_KERNEL_ARGS")) {
+        size_t maxLocal = 0;
+        clGetKernelWorkGroupInfo(deviceFunction, openClDeviceId, CL_KERNEL_WORK_GROUP_SIZE,
+                                  sizeof(maxLocal), &maxLocal, NULL);
+        size_t prefMultiple = 0;
+        clGetKernelWorkGroupInfo(deviceFunction, openClDeviceId, CL_KERNEL_PREFERRED_WORK_GROUP_SIZE_MULTIPLE,
+                                  sizeof(prefMultiple), &prefMultiple, NULL);
+        fprintf(stderr, "[LaunchKernel] %s: CL_KERNEL_WORK_GROUP_SIZE=%zu preferredMultiple=%zu (requested local=%zu)\n",
+                kernelNameBuf, maxLocal, prefMultiple, localWorkSize[0]);
+        fflush(stderr);
+    }
+
+    if (getenv("BEAGLE_DEBUG_KERNEL_ARGS")) {
+        fprintf(stderr, "[LaunchKernel] %s: about to enqueue\n", kernelNameBuf);
+        fflush(stderr);
+    }
+
     if (globalWorkSize[1] == 1 && globalWorkSize[2] == 1) {
         SAFE_CL(clEnqueueNDRangeKernel(openClCommandQueues[0], deviceFunction, 1, NULL,
                                        globalWorkSize, localWorkSize, 0, NULL, NULL));
@@ -590,6 +630,14 @@ void GPUInterface::LaunchKernel(GPUFunction deviceFunction,
                                        globalWorkSize, localWorkSize, 0, NULL, NULL));
     }
 
+    if (getenv("BEAGLE_DEBUG_KERNEL_ARGS")) {
+        fprintf(stderr, "[LaunchKernel] %s: enqueued OK, calling clFinish\n", kernelNameBuf);
+        fflush(stderr);
+        clFinish(openClCommandQueues[0]);
+        fprintf(stderr, "[LaunchKernel] %s: clFinish returned OK\n", kernelNameBuf);
+        fflush(stderr);
+    }
+
 #ifdef BEAGLE_DEBUG_FLOW
     fprintf(stderr,"\t\t\tLeaving  GPUInterface::LaunchKernel\n");
 #endif
@@ -607,15 +655,30 @@ void GPUInterface::LaunchKernelConcurrent(GPUFunction deviceFunction,
     fprintf(stderr,"\t\t\tEntering GPUInterface::LaunchKernel\n");
 #endif
 
+    char kernelNameBuf[256] = {0};
+    clGetKernelInfo(deviceFunction, CL_KERNEL_FUNCTION_NAME, sizeof(kernelNameBuf), kernelNameBuf, NULL);
+    if (getenv("BEAGLE_DEBUG_KERNEL_ARGS")) {
+        fprintf(stderr, "[LaunchKernelConcurrent] %s: parameterCountV=%d totalParameterCount=%d streamIndex=%d waitIndex=%d\n",
+                kernelNameBuf, parameterCountV, totalParameterCount, streamIndex, waitIndex);
+    }
+
     va_list parameters;
     va_start(parameters, totalParameterCount);
     for(int i = 0; i < parameterCountV; i++) {
         void* param = (void*)(size_t)va_arg(parameters, GPUPtr);
+        if (getenv("BEAGLE_DEBUG_KERNEL_ARGS")) {
+            fprintf(stderr, "[LaunchKernelConcurrent] %s: ptr arg[%d] = %p\n", kernelNameBuf, i, param);
+            fflush(stderr);
+        }
 
         SAFE_CL(clSetKernelArg(deviceFunction, i, sizeof(param), &param));
     }
     for(int i = parameterCountV; i < totalParameterCount; i++) {
         unsigned int param = va_arg(parameters, unsigned int);
+        if (getenv("BEAGLE_DEBUG_KERNEL_ARGS")) {
+            fprintf(stderr, "[LaunchKernelConcurrent] %s: int arg[%d] = %u (signed: %d)\n", kernelNameBuf, i, param, (int)param);
+            fflush(stderr);
+        }
 
         SAFE_CL(clSetKernelArg(deviceFunction, i, sizeof(unsigned int), &param));
     }
@@ -632,6 +695,14 @@ void GPUInterface::LaunchKernelConcurrent(GPUFunction deviceFunction,
     globalWorkSize[1] = block.y * grid.y;
     globalWorkSize[2] = block.z * grid.z;
 
+    if (getenv("BEAGLE_DEBUG_KERNEL_ARGS")) {
+        fprintf(stderr, "[LaunchKernelConcurrent] %s: block=(%d,%d,%d) grid=(%d,%d,%d) local=(%zu,%zu,%zu) global=(%zu,%zu,%zu)\n",
+                kernelNameBuf, block.x, block.y, block.z, grid.x, grid.y, grid.z,
+                localWorkSize[0], localWorkSize[1], localWorkSize[2],
+                globalWorkSize[0], globalWorkSize[1], globalWorkSize[2]);
+        fflush(stderr);
+    }
+
 #ifdef BEAGLE_DEBUG_VALUES
     for (int i=0; i<3; i++) {
         printf("localWorkSize[%d]  = %lu\n", i, localWorkSize[i]);
@@ -649,6 +720,11 @@ void GPUInterface::LaunchKernelConcurrent(GPUFunction deviceFunction,
         dims = 2;
     }
 
+    if (getenv("BEAGLE_DEBUG_KERNEL_ARGS")) {
+        fprintf(stderr, "[LaunchKernelConcurrent] %s: about to enqueue (dims=%d)\n", kernelNameBuf, dims);
+        fflush(stderr);
+    }
+
     // if (streamIndex != -1) {
     //     streamIndex /= 2;
     //     // if (waitIndex != -1) {
@@ -688,6 +764,13 @@ void GPUInterface::LaunchKernelConcurrent(GPUFunction deviceFunction,
                                        globalWorkSize, localWorkSize,
                                        0, NULL, NULL));
     // }
+    if (getenv("BEAGLE_DEBUG_KERNEL_ARGS")) {
+        fprintf(stderr, "[LaunchKernelConcurrent] %s: enqueued OK, calling clFinish\n", kernelNameBuf);
+        fflush(stderr);
+        clFinish(openClCommandQueues[0]);
+        fprintf(stderr, "[LaunchKernelConcurrent] %s: clFinish returned OK\n", kernelNameBuf);
+        fflush(stderr);
+    }
 #ifdef BEAGLE_DEBUG_FLOW
     fprintf(stderr,"\t\t\tLeaving  GPUInterface::LaunchKernel\n");
 #endif
diff --git a/libhmsbeagle/GPU/OpenCLPlugin.cpp b/libhmsbeagle/GPU/OpenCLPlugin.cpp
index a36e49d9..3bfdce1d 100644
--- a/libhmsbeagle/GPU/OpenCLPlugin.cpp
+++ b/libhmsbeagle/GPU/OpenCLPlugin.cpp
@@ -88,7 +88,7 @@ OpenCLPlugin::~OpenCLPlugin()
 
 
 extern "C" {
-void* plugin_init(void){
+__attribute__((visibility("default"))) void* plugin_init(void){
 	return new beagle::gpu::OpenCLPlugin();
 }
 }
diff --git a/libhmsbeagle/GPU/OpenCLSpectralPlugin.cpp b/libhmsbeagle/GPU/OpenCLSpectralPlugin.cpp
index 65d52b5a..7d0544d1 100644
--- a/libhmsbeagle/GPU/OpenCLSpectralPlugin.cpp
+++ b/libhmsbeagle/GPU/OpenCLSpectralPlugin.cpp
@@ -84,7 +84,7 @@ OpenCLSpectralPlugin::~OpenCLSpectralPlugin()
 
 
 extern "C" {
-void* plugin_init(void) {
+__attribute__((visibility("default"))) void* plugin_init(void) {
     return new beagle::gpu::OpenCLSpectralPlugin();
 }
 }
diff --git a/libhmsbeagle/GPU/kernels/kernelsSpectralIfDef.cu b/libhmsbeagle/GPU/kernels/kernelsSpectralIfDef.cu
index 7b359fcf..96feff1c 100644
--- a/libhmsbeagle/GPU/kernels/kernelsSpectralIfDef.cu
+++ b/libhmsbeagle/GPU/kernels/kernelsSpectralIfDef.cu
@@ -859,25 +859,44 @@ KW_GLOBAL_KERNEL void kernelPartialsStatesGrowingTopRootSpectral(
 
 #define ADJOINT_BLOCK_SP_N  128
 
+/* CAS-retry float/double atomic add, in its own KW_DEVICE_FUNC-free OpenCL
+ * helper function rather than inlined as a macro at the call site.
+ *
+ * This split is load-bearing, not stylistic: with the CAS loop inlined
+ * directly into kernelAdjointMergedN (compiled at PADDED_STATE_COUNT=32,
+ * i.e. any state count that pads above 16 states, e.g. the 17-state case),
+ * clEnqueueNDRangeKernel succeeds but clFinish() never returns — reproduced
+ * down to a single block (no cross-block contention at all) and shown by
+ * bisection to require *both* things at once: the atomic_cmpxchg call AND a
+ * retry loop around it (a lone, non-looping atomic_cmpxchg — or the same
+ * loop shape swapped for a different atomic like atomic_xchg — completes
+ * fine; only "loop whose exit condition reads back an atomic's result"
+ * hangs, regardless of do/while vs for(;;) syntax). That is consistent with
+ * a known class of ROCm/LLVM AMDGPU backend bugs (this is ROCm 4.0.1, circa
+ * 2020) around convergent-instruction handling for loops inlined into large,
+ * register-heavy kernels; it is not a BEAGLE logic bug — the 16-state
+ * (unpadded, S=16) build of the identical macro works correctly. Moving the
+ * loop into its own `__attribute__((noinline))` function keeps it out of
+ * the miscompiled inlined context and reliably avoids the hang. */
 #ifndef ADJOINT_ATOMIC_ADD_GPU  /* may already be defined if both IfDef files merged */
 #ifdef CUDA
 #define ADJOINT_ATOMIC_ADD_GPU(ptr, val)  atomicAdd((ptr), (REAL)(val))
 #elif defined(DOUBLE_PRECISION)
+__attribute__((noinline)) void adjointAtomicAddGpuDPHelper(__global long* _anp, double val) {
+    long _ano, _ann;
+    do { _ano = *_anp; _ann = as_long(as_double(_ano) + val); }
+    while (atom_cmpxchg(_anp, _ano, _ann) != _ano);
+}
 #define ADJOINT_ATOMIC_ADD_GPU(ptr, val) \
-    do { \
-        __global long* _anp = (__global long*)(ptr); \
-        long _ano, _ann; \
-        do { _ano = *_anp; _ann = as_long(as_double(_ano) + (double)(val)); } \
-        while (atom_cmpxchg(_anp, _ano, _ann) != _ano); \
-    } while(0)
+    adjointAtomicAddGpuDPHelper((__global long*)(ptr), (double)(val))
 #else
+__attribute__((noinline)) void adjointAtomicAddGpuSPHelper(__global int* _anp, float val) {
+    int _ano, _ann;
+    do { _ano = *_anp; _ann = as_int(as_float(_ano) + val); }
+    while (atomic_cmpxchg(_anp, _ano, _ann) != _ano);
+}
 #define ADJOINT_ATOMIC_ADD_GPU(ptr, val) \
-    do { \
-        __global int* _anp = (__global int*)(ptr); \
-        int _ano, _ann; \
-        do { _ano = *_anp; _ann = as_int(as_float(_ano) + (float)(val)); } \
-        while (atomic_cmpxchg(_anp, _ano, _ann) != _ano); \
-    } while(0)
+    adjointAtomicAddGpuSPHelper((__global int*)(ptr), (float)(val))
 #endif
 #endif /* ADJOINT_ATOMIC_ADD_GPU */
 
@@ -910,37 +929,56 @@ KW_GLOBAL_KERNEL void kernelPartialsStatesGrowingTopRootSpectral(
  * ievc rotation itself happens exactly once per block afterward
  * (ADJOINTN_ROTATE_ROW) instead of once per (block, pattern) — valid by
  * linearity since `ievc` is constant across every pattern in a segment. */
-#define ADJOINTN_ACCUM_ROW_RAW(EVECT_COL, IS_STATES, REGOP, RAWROW) \
-    { \
-        for (int _rs_a = 0; _rs_a < PADDED_STATE_COUNT; _rs_a++) (REGOP)[_rs_a] = (REAL)0; \
-        for (int _k_a = tid; _k_a < totalPatterns; _k_a += ADJOINT_BLOCK_SP_N) { \
-            REAL _lhs_a = (REAL)0; \
-            for (int _j_a = 0; _j_a < PADDED_STATE_COUNT; _j_a++) \
-                SPECTRAL_FMA((EVECT_COL)[_j_a], prePartials[catOff + _k_a*PADDED_STATE_COUNT + _j_a], _lhs_a); \
-            const REAL _lv_a = _lhs_a * (patternWeights[_k_a] * sCatW / exp(perSiteLikelihoods[_k_a])); \
-            if (IS_STATES) { \
-                const int _st_a = tipStates[_k_a]; \
-                if (_st_a >= PADDED_STATE_COUNT) { \
-                    for (int _rs_a = 0; _rs_a < PADDED_STATE_COUNT; _rs_a++) (REGOP)[_rs_a] += _lv_a; \
-                } else { \
-                    (REGOP)[_st_a] += _lv_a; \
-                } \
-            } else { \
-                for (int _rs_a = 0; _rs_a < PADDED_STATE_COUNT; _rs_a++) \
-                    (REGOP)[_rs_a] += _lv_a * postPartials[catOff + _k_a*PADDED_STATE_COUNT + _rs_a]; \
-            } \
-        } \
-        for (int _rs_a = 0; _rs_a < PADDED_STATE_COUNT; _rs_a++) { \
-            sRedBuf[tid] = (REGOP)[_rs_a]; \
-            KW_LOCAL_FENCE; \
-            for (int _str_a = ADJOINT_BLOCK_SP_N >> 1; _str_a >= 1; _str_a >>= 1) { \
-                if (tid < _str_a) sRedBuf[tid] += sRedBuf[tid + _str_a]; \
-                KW_LOCAL_FENCE; \
-            } \
-            if (tid == 0) (RAWROW)[_rs_a] = sRedBuf[0]; \
-            KW_LOCAL_FENCE; \
-        } \
+/* Extracted into its own `__attribute__((noinline))` helper (§5.4 of
+ * GPU_ADJOINT_PLAN.md) — this is the single densest concentration of
+ * `tid`-predicated, barrier-synchronized code in the kernel (a 32x-repeated
+ * 128-wide tree reduction), shared by every code path (isAllReal, singleton,
+ * pairleader alike) and never previously isolated by the noinline treatment
+ * already proven for the SINGLETON/PAIRLEADER bodies above. Giving it its
+ * own function/register-allocation scope is a diagnostic for whether it is
+ * the specific miscompiled hot spot behind the PADDED_STATE_COUNT=32
+ * get_local_id(0) corruption documented in CLUSTER_AGENT_FINDINGS.md §3. */
+__attribute__((noinline)) void adjointAccumRowRaw(
+        int tid, __local REAL* evectCol, bool isStates,
+        REAL* regOp, __local REAL* rawRow,
+        __global REAL* prePartials, __global REAL* postPartials,
+        __global int* tipStates, __global REAL* patternWeights, REAL sCatW,
+        __global REAL* perSiteLikelihoods, __local REAL* sRedBuf,
+        int totalPatterns, int catOff) {
+    for (int _rs_a = 0; _rs_a < PADDED_STATE_COUNT; _rs_a++) regOp[_rs_a] = (REAL)0;
+    for (int _k_a = tid; _k_a < totalPatterns; _k_a += ADJOINT_BLOCK_SP_N) {
+        REAL _lhs_a = (REAL)0;
+        for (int _j_a = 0; _j_a < PADDED_STATE_COUNT; _j_a++)
+            SPECTRAL_FMA(evectCol[_j_a], prePartials[catOff + _k_a*PADDED_STATE_COUNT + _j_a], _lhs_a);
+        const REAL _lv_a = _lhs_a * (patternWeights[_k_a] * sCatW / exp(perSiteLikelihoods[_k_a]));
+        if (isStates) {
+            const int _st_a = tipStates[_k_a];
+            if (_st_a >= PADDED_STATE_COUNT) {
+                for (int _rs_a = 0; _rs_a < PADDED_STATE_COUNT; _rs_a++) regOp[_rs_a] += _lv_a;
+            } else {
+                regOp[_st_a] += _lv_a;
+            }
+        } else {
+            for (int _rs_a = 0; _rs_a < PADDED_STATE_COUNT; _rs_a++)
+                regOp[_rs_a] += _lv_a * postPartials[catOff + _k_a*PADDED_STATE_COUNT + _rs_a];
+        }
     }
+    for (int _rs_a = 0; _rs_a < PADDED_STATE_COUNT; _rs_a++) {
+        sRedBuf[tid] = regOp[_rs_a];
+        KW_LOCAL_FENCE;
+        for (int _str_a = ADJOINT_BLOCK_SP_N >> 1; _str_a >= 1; _str_a >>= 1) {
+            if (tid < _str_a) sRedBuf[tid] += sRedBuf[tid + _str_a];
+            KW_LOCAL_FENCE;
+        }
+        if (tid == 0) rawRow[_rs_a] = sRedBuf[0];
+        KW_LOCAL_FENCE;
+    }
+}
+
+#define ADJOINTN_ACCUM_ROW_RAW(EVECT_COL, IS_STATES, REGOP, RAWROW) \
+    adjointAccumRowRaw(tid, EVECT_COL, IS_STATES, REGOP, RAWROW, \
+        prePartials, postPartials, tipStates, patternWeights, sCatW, \
+        perSiteLikelihoods, sRedBuf, totalPatterns, catOff);
 
 /* Rotate a raw (un-projected) row RAWROW[PADDED_STATE_COUNT] through `ievc`
  * once: ROTATED[rs'] = sum_rs RAWROW[rs] * ievc[rs*S + rs']. Valid for every
@@ -950,20 +988,26 @@ KW_GLOBAL_KERNEL void kernelPartialsStatesGrowingTopRootSpectral(
  * block (not per pattern) — mathematically equivalent to the old
  * per-pattern `ievc` projection by linearity, since `ievc` is constant
  * across every pattern within one segment/branch. */
-#define ADJOINTN_ROTATE_ROW(RAWROW, ROTATED) \
-    { \
-        const int _rrpt = (PADDED_STATE_COUNT + ADJOINT_BLOCK_SP_N - 1) / ADJOINT_BLOCK_SP_N; \
-        for (int _m = 0; _m < _rrpt; _m++) { \
-            const int _rsp = tid + _m * ADJOINT_BLOCK_SP_N; \
-            if (_rsp < PADDED_STATE_COUNT) { \
-                REAL _acc = (REAL)0; \
-                for (int _rs = 0; _rs < PADDED_STATE_COUNT; _rs++) \
-                    SPECTRAL_FMA((RAWROW)[_rs], ievc[_rs*PADDED_STATE_COUNT + _rsp], _acc); \
-                (ROTATED)[_rsp] = _acc; \
-            } \
-        } \
-        KW_LOCAL_FENCE; \
+/* Extracted into its own noinline helper for the same reason and as part of
+ * the same §5.4 diagnostic as adjointAccumRowRaw() above — smaller than that
+ * one but shares the `tid`-indexed pattern with its own KW_LOCAL_FENCE. */
+__attribute__((noinline)) void adjointRotateRow(
+        int tid, __local REAL* rawRow, __global REAL* ievc, __local REAL* rotated) {
+    const int _rrpt = (PADDED_STATE_COUNT + ADJOINT_BLOCK_SP_N - 1) / ADJOINT_BLOCK_SP_N;
+    for (int _m = 0; _m < _rrpt; _m++) {
+        const int _rsp = tid + _m * ADJOINT_BLOCK_SP_N;
+        if (_rsp < PADDED_STATE_COUNT) {
+            REAL _acc = (REAL)0;
+            for (int _rs = 0; _rs < PADDED_STATE_COUNT; _rs++)
+                SPECTRAL_FMA(rawRow[_rs], ievc[_rs*PADDED_STATE_COUNT + _rsp], _acc);
+            rotated[_rsp] = _acc;
+        }
     }
+    KW_LOCAL_FENCE;
+}
+
+#define ADJOINTN_ROTATE_ROW(RAWROW, ROTATED) \
+    adjointRotateRow(tid, RAWROW, ievc, ROTATED);
 
 /* Apply the all-real integral transform to a single reduced row DEST[S] and
  * atomicAdd into dGradient row `ROW`. Verbatim math from
@@ -983,36 +1027,133 @@ KW_GLOBAL_KERNEL void kernelPartialsStatesGrowingTopRootSpectral(
 /* Apply the complex-eigenvalue integral transform. Verbatim math from
  * kernelAdjointPhase2ComplexN's two branches, tid==0 only, reading the
  * shared reduced rows sOpRow0[S] (row `ls`) and, for the pair-leader case,
- * sOpRow1[S] (row `ls+1`) instead of dOpBuf. */
+ * sOpRow1[S] (row `ls+1`) instead of dOpBuf.
+ *
+ * Each branch's per-iteration body is pulled into its own
+ * `__attribute__((noinline))` helper function, the same treatment already
+ * applied to the atomic CAS-retry loop above and for the same reason: these
+ * bodies carry a lot of REAL temporaries (the complex-pair branch alone has
+ * ~30), and inlining them directly into the `for` loop inside
+ * kernelAdjointMergedN (compiled at PADDED_STATE_COUNT=32) was observed to
+ * produce silently-wrong (NaN) results that turned into a hard GPU page
+ * fault the moment unrelated code (e.g. a debug printf) was added nearby —
+ * i.e. correctness that depended on incidental register allocation/code
+ * layout, the signature of the ROCm 4.0.1 AMDGPU-backend miscompilation
+ * class documented above, not a BEAGLE logic bug. Keeping only the cheap
+ * loop control (the `for` and the `_ri` branch dispatch) inlined and moving
+ * the heavy arithmetic out-of-line avoids it, matching the CAS-loop fix. */
+__attribute__((noinline)) void adjointSingletonRealStep(
+        __global REAL* dGradient, int _S_c, int ls, int _rs_c, REAL t,
+        REAL _ea_c, REAL _la_c,
+        __local REAL* sEvalR, __local REAL* sExpat, __local REAL* sOpRow0) {
+    const REAL _co_c = (t * fabs(_la_c - sEvalR[_rs_c]) < (REAL)1e-12)
+        ? t * _ea_c : (_ea_c - sExpat[_rs_c]) / (_la_c - sEvalR[_rs_c]);
+    ADJOINT_ATOMIC_ADD_GPU(&dGradient[ls*_S_c+_rs_c], sOpRow0[_rs_c] * _co_c);
+}
+
+__attribute__((noinline)) void adjointSingletonComplexStep(
+        __global REAL* dGradient, int _S_c, int ls, int _rs_c, REAL t,
+        REAL _ea_c, REAL _la_c, REAL _ri_c,
+        __local REAL* sEvalR, __local REAL* sExpat, __local REAL* sCosbt,
+        __local REAL* sSinbt, __local REAL* sOpRow0) {
+    const REAL _sr_c = sEvalR[_rs_c] - _la_c;
+    const REAL _dn_c = _sr_c*_sr_c + _ri_c*_ri_c;
+    REAL _i0_c, _i1_c;
+    if (_dn_c < (REAL)1e-12) { _i0_c = t; _i1_c = (REAL)0; }
+    else {
+        const REAL _ex_c = sExpat[_rs_c] / _ea_c;
+        _i0_c = (_ex_c*(_sr_c*sCosbt[_rs_c]+_ri_c*sSinbt[_rs_c])-_sr_c)/_dn_c;
+        _i1_c = (_ex_c*(_sr_c*sSinbt[_rs_c]-_ri_c*sCosbt[_rs_c])+_ri_c)/_dn_c;
+    }
+    const REAL _c0_c = _ea_c*_i0_c, _c1_c = _ea_c*_i1_c;
+    const REAL _n0_c = sOpRow0[_rs_c], _n1_c = sOpRow0[_rs_c+1];
+    ADJOINT_ATOMIC_ADD_GPU(&dGradient[ls*_S_c+_rs_c],     _c0_c*_n0_c+_c1_c*_n1_c);
+    ADJOINT_ATOMIC_ADD_GPU(&dGradient[ls*_S_c+_rs_c+1], -_c1_c*_n0_c+_c0_c*_n1_c);
+}
+
+__attribute__((noinline)) void adjointSingletonLoop(
+        int tid, int ls, REAL t, __global REAL* dGradient,
+        __local REAL* sEvalR, __local REAL* sEvalI, __local REAL* sExpat,
+        __local REAL* sCosbt, __local REAL* sSinbt, __local REAL* sOpRow0) {
+    if (tid != 0) return;
+    const int _S_c = PADDED_STATE_COUNT;
+    const REAL _ea_c = sExpat[ls], _la_c = sEvalR[ls];
+    for (int _rs_c = 0; _rs_c < _S_c; ) {
+        const REAL _ri_c = sEvalI[_rs_c];
+        if (_ri_c == (REAL)0) {
+            adjointSingletonRealStep(dGradient, _S_c, ls, _rs_c, t, _ea_c, _la_c,
+                                      sEvalR, sExpat, sOpRow0);
+            _rs_c++;
+        } else {
+            adjointSingletonComplexStep(dGradient, _S_c, ls, _rs_c, t, _ea_c, _la_c, _ri_c,
+                                         sEvalR, sExpat, sCosbt, sSinbt, sOpRow0);
+            _rs_c += 2;
+        }
+    }
+}
+
 #define ADJOINTN_APPLY_COMPLEX_SINGLETON() \
-    if (tid == 0) { \
-        const int _S_c = PADDED_STATE_COUNT; \
-        const REAL _ea_c = sExpat[ls], _la_c = sEvalR[ls]; \
-        for (int _rs_c = 0; _rs_c < _S_c; ) { \
-            const REAL _ri_c = sEvalI[_rs_c]; \
-            if (_ri_c == (REAL)0) { \
-                const REAL _co_c = (t * fabs(_la_c - sEvalR[_rs_c]) < (REAL)1e-12) \
-                    ? t * _ea_c : (_ea_c - sExpat[_rs_c]) / (_la_c - sEvalR[_rs_c]); \
-                ADJOINT_ATOMIC_ADD_GPU(&dGradient[ls*_S_c+_rs_c], sOpRow0[_rs_c] * _co_c); \
-                _rs_c++; \
-            } else { \
-                const REAL _sr_c = sEvalR[_rs_c] - _la_c; \
-                const REAL _dn_c = _sr_c*_sr_c + _ri_c*_ri_c; \
-                REAL _i0_c, _i1_c; \
-                if (_dn_c < (REAL)1e-12) { _i0_c = t; _i1_c = (REAL)0; } \
-                else { \
-                    const REAL _ex_c = sExpat[_rs_c] / _ea_c; \
-                    _i0_c = (_ex_c*(_sr_c*sCosbt[_rs_c]+_ri_c*sSinbt[_rs_c])-_sr_c)/_dn_c; \
-                    _i1_c = (_ex_c*(_sr_c*sSinbt[_rs_c]-_ri_c*sCosbt[_rs_c])+_ri_c)/_dn_c; \
-                } \
-                const REAL _c0_c = _ea_c*_i0_c, _c1_c = _ea_c*_i1_c; \
-                const REAL _n0_c = sOpRow0[_rs_c], _n1_c = sOpRow0[_rs_c+1]; \
-                ADJOINT_ATOMIC_ADD_GPU(&dGradient[ls*_S_c+_rs_c],     _c0_c*_n0_c+_c1_c*_n1_c); \
-                ADJOINT_ATOMIC_ADD_GPU(&dGradient[ls*_S_c+_rs_c+1], -_c1_c*_n0_c+_c0_c*_n1_c); \
-                _rs_c += 2; \
-            } \
-        } \
+    adjointSingletonLoop(tid, ls, t, dGradient, sEvalR, sEvalI, sExpat, sCosbt, sSinbt, sOpRow0);
+
+__attribute__((noinline)) void adjointPairleaderRealStep(
+        __global REAL* dGradient, int _S_q, int ls, int _rs_q, REAL t,
+        REAL _lr_q, REAL _li_q, REAL _ec_q, REAL _es_q, REAL _cI_q, REAL _sI_q, REAL _ea_q,
+        __local REAL* sEvalR, __local REAL* sExpat,
+        __local REAL* sOpRow0, __local REAL* sOpRow1) {
+    const REAL _sr_q = sEvalR[_rs_q] - _lr_q;
+    const REAL _dn_q = _sr_q*_sr_q + _li_q*_li_q;
+    REAL _i0_q, _i1_q;
+    if (_dn_q < (REAL)1e-12) { _i0_q = t; _i1_q = (REAL)0; }
+    else {
+        const REAL _ex_q = sExpat[_rs_q] / _ea_q;
+        _i0_q = (_ex_q*(_sr_q*_cI_q+_li_q*_sI_q)-_sr_q)/_dn_q;
+        _i1_q = (_ex_q*(_sr_q*_sI_q-_li_q*_cI_q)+_li_q)/_dn_q;
     }
+    const REAL _p0_q=_ec_q*_i0_q+_es_q*_i1_q, _p1_q=_ec_q*_i1_q-_es_q*_i0_q;
+    const REAL _p2_q=_es_q*_i0_q-_ec_q*_i1_q, _p3_q=_es_q*_i1_q+_ec_q*_i0_q;
+    const REAL _n0_q=sOpRow0[_rs_q], _n1_q=sOpRow1[_rs_q];
+    ADJOINT_ATOMIC_ADD_GPU(&dGradient[ls*_S_q+_rs_q],     _p0_q*_n0_q+_p1_q*_n1_q);
+    ADJOINT_ATOMIC_ADD_GPU(&dGradient[(ls+1)*_S_q+_rs_q], _p2_q*_n0_q+_p3_q*_n1_q);
+}
+
+__attribute__((noinline)) void adjointPairleaderComplexStep(
+        __global REAL* dGradient, int _S_q, int ls, int _rs_q, REAL t,
+        REAL _lr_q, REAL _li_q, REAL _ec_q, REAL _es_q, REAL _cI_q, REAL _sI_q, REAL _ea_q,
+        REAL _ri_q,
+        __local REAL* sEvalR, __local REAL* sExpat, __local REAL* sCosbt, __local REAL* sSinbt,
+        __local REAL* sOpRow0, __local REAL* sOpRow1) {
+    const REAL _rr_q=sEvalR[_rs_q], _ri2_q=_ri_q;
+    const REAL _sr_q=_rr_q-_lr_q, _si1_q=_li_q+_ri2_q, _si2_q=_ri2_q-_li_q;
+    const REAL _sr2_q=_sr_q*_sr_q;
+    const REAL _d1_q=_sr2_q+_si1_q*_si1_q, _d2_q=_sr2_q+_si2_q*_si2_q;
+    const REAL _ex_q=(_d1_q>=(REAL)1e-12||_d2_q>=(REAL)1e-12) ? sExpat[_rs_q]/_ea_q : (REAL)0;
+    const REAL _clcr_q=_cI_q*sCosbt[_rs_q], _slsr_q=_sI_q*sSinbt[_rs_q];
+    const REAL _clsr_q=_cI_q*sSinbt[_rs_q], _slcr_q=_sI_q*sCosbt[_rs_q];
+    REAL _i1r_q,_i1i_q;
+    if (_d1_q<(REAL)1e-12){_i1r_q=t;_i1i_q=(REAL)0;}
+    else{
+        const REAL _cs1_q=_clcr_q-_slsr_q, _sn1_q=_slcr_q+_clsr_q;
+        _i1r_q=(_sr_q*(_ex_q*_cs1_q-(REAL)1)+_si1_q*_ex_q*_sn1_q)/_d1_q;
+        _i1i_q=(_sr_q*_ex_q*_sn1_q-_si1_q*(_ex_q*_cs1_q-(REAL)1))/_d1_q;
+    }
+    REAL _i2r_q,_i2i_q;
+    if (_d2_q<(REAL)1e-12){_i2r_q=t;_i2i_q=(REAL)0;}
+    else{
+        const REAL _cs2_q=_clcr_q+_slsr_q, _sn2_q=_clsr_q-_slcr_q;
+        _i2r_q=(_sr_q*(_ex_q*_cs2_q-(REAL)1)+_si2_q*_ex_q*_sn2_q)/_d2_q;
+        _i2i_q=(_sr_q*_ex_q*_sn2_q-_si2_q*(_ex_q*_cs2_q-(REAL)1))/_d2_q;
+    }
+    const REAL _pr_q=_ec_q*_i1r_q+_es_q*_i1i_q, _pi_q=_ec_q*_i1i_q-_es_q*_i1r_q;
+    const REAL _mr_q=_ec_q*_i2r_q-_es_q*_i2i_q, _mi_q=_ec_q*_i2i_q+_es_q*_i2r_q;
+    const REAL _A_q=(REAL)0.5*(_mr_q+_pr_q), _B_q=(REAL)0.5*(_mi_q+_pi_q);
+    const REAL _C_q=(REAL)0.5*(_pi_q-_mi_q), _D_q=(REAL)0.5*(_mr_q-_pr_q);
+    const REAL _n00_q=sOpRow0[_rs_q],  _n01_q=sOpRow0[_rs_q+1];
+    const REAL _n10_q=sOpRow1[_rs_q], _n11_q=sOpRow1[_rs_q+1];
+    ADJOINT_ATOMIC_ADD_GPU(&dGradient[ls*_S_q+_rs_q],          _A_q*_n00_q+_B_q*_n01_q+_C_q*_n10_q+_D_q*_n11_q);
+    ADJOINT_ATOMIC_ADD_GPU(&dGradient[ls*_S_q+_rs_q+1],       -_B_q*_n00_q+_A_q*_n01_q-_D_q*_n10_q+_C_q*_n11_q);
+    ADJOINT_ATOMIC_ADD_GPU(&dGradient[(ls+1)*_S_q+_rs_q],     -_C_q*_n00_q-_D_q*_n01_q+_A_q*_n10_q+_B_q*_n11_q);
+    ADJOINT_ATOMIC_ADD_GPU(&dGradient[(ls+1)*_S_q+_rs_q+1],    _D_q*_n00_q-_C_q*_n01_q-_B_q*_n10_q+_A_q*_n11_q);
+}
 
 #define ADJOINTN_APPLY_COMPLEX_PAIRLEADER() \
     if (tid == 0) { \
@@ -1024,54 +1165,17 @@ KW_GLOBAL_KERNEL void kernelPartialsStatesGrowingTopRootSpectral(
         for (int _rs_q = 0; _rs_q < _S_q; ) { \
             const REAL _ri_q = sEvalI[_rs_q]; \
             if (_ri_q == (REAL)0) { \
-                const REAL _sr_q = sEvalR[_rs_q] - _lr_q; \
-                const REAL _dn_q = _sr_q*_sr_q + _li_q*_li_q; \
-                REAL _i0_q, _i1_q; \
-                if (_dn_q < (REAL)1e-12) { _i0_q = t; _i1_q = (REAL)0; } \
-                else { \
-                    const REAL _ex_q = sExpat[_rs_q] / _ea_q; \
-                    _i0_q = (_ex_q*(_sr_q*_cI_q+_li_q*_sI_q)-_sr_q)/_dn_q; \
-                    _i1_q = (_ex_q*(_sr_q*_sI_q-_li_q*_cI_q)+_li_q)/_dn_q; \
-                } \
-                const REAL _p0_q=_ec_q*_i0_q+_es_q*_i1_q, _p1_q=_ec_q*_i1_q-_es_q*_i0_q; \
-                const REAL _p2_q=_es_q*_i0_q-_ec_q*_i1_q, _p3_q=_es_q*_i1_q+_ec_q*_i0_q; \
-                const REAL _n0_q=sOpRow0[_rs_q], _n1_q=sOpRow1[_rs_q]; \
-                ADJOINT_ATOMIC_ADD_GPU(&dGradient[ls*_S_q+_rs_q],     _p0_q*_n0_q+_p1_q*_n1_q); \
-                ADJOINT_ATOMIC_ADD_GPU(&dGradient[(ls+1)*_S_q+_rs_q], _p2_q*_n0_q+_p3_q*_n1_q); \
+                adjointPairleaderRealStep(dGradient, _S_q, ls, _rs_q, t, \
+                                           _lr_q, _li_q, _ec_q, _es_q, _cI_q, _sI_q, _ea_q, \
+                                           sEvalR, sExpat, sOpRow0, sOpRow1); \
                 _rs_q++; \
-            } else { \
-                const REAL _rr_q=sEvalR[_rs_q], _ri2_q=_ri_q; \
-                const REAL _sr_q=_rr_q-_lr_q, _si1_q=_li_q+_ri2_q, _si2_q=_ri2_q-_li_q; \
-                const REAL _sr2_q=_sr_q*_sr_q; \
-                const REAL _d1_q=_sr2_q+_si1_q*_si1_q, _d2_q=_sr2_q+_si2_q*_si2_q; \
-                const REAL _ex_q=(_d1_q>=(REAL)1e-12||_d2_q>=(REAL)1e-12) ? sExpat[_rs_q]/_ea_q : (REAL)0; \
-                const REAL _clcr_q=_cI_q*sCosbt[_rs_q], _slsr_q=_sI_q*sSinbt[_rs_q]; \
-                const REAL _clsr_q=_cI_q*sSinbt[_rs_q], _slcr_q=_sI_q*sCosbt[_rs_q]; \
-                REAL _i1r_q,_i1i_q; \
-                if (_d1_q<(REAL)1e-12){_i1r_q=t;_i1i_q=(REAL)0;} \
-                else{ \
-                    const REAL _cs1_q=_clcr_q-_slsr_q, _sn1_q=_slcr_q+_clsr_q; \
-                    _i1r_q=(_sr_q*(_ex_q*_cs1_q-(REAL)1)+_si1_q*_ex_q*_sn1_q)/_d1_q; \
-                    _i1i_q=(_sr_q*_ex_q*_sn1_q-_si1_q*(_ex_q*_cs1_q-(REAL)1))/_d1_q; \
-                } \
-                REAL _i2r_q,_i2i_q; \
-                if (_d2_q<(REAL)1e-12){_i2r_q=t;_i2i_q=(REAL)0;} \
-                else{ \
-                    const REAL _cs2_q=_clcr_q+_slsr_q, _sn2_q=_clsr_q-_slcr_q; \
-                    _i2r_q=(_sr_q*(_ex_q*_cs2_q-(REAL)1)+_si2_q*_ex_q*_sn2_q)/_d2_q; \
-                    _i2i_q=(_sr_q*_ex_q*_sn2_q-_si2_q*(_ex_q*_cs2_q-(REAL)1))/_d2_q; \
-                } \
-                const REAL _pr_q=_ec_q*_i1r_q+_es_q*_i1i_q, _pi_q=_ec_q*_i1i_q-_es_q*_i1r_q; \
-                const REAL _mr_q=_ec_q*_i2r_q-_es_q*_i2i_q, _mi_q=_ec_q*_i2i_q+_es_q*_i2r_q; \
-                const REAL _A_q=(REAL)0.5*(_mr_q+_pr_q), _B_q=(REAL)0.5*(_mi_q+_pi_q); \
-                const REAL _C_q=(REAL)0.5*(_pi_q-_mi_q), _D_q=(REAL)0.5*(_mr_q-_pr_q); \
-                const REAL _n00_q=sOpRow0[_rs_q],  _n01_q=sOpRow0[_rs_q+1]; \
-                const REAL _n10_q=sOpRow1[_rs_q], _n11_q=sOpRow1[_rs_q+1]; \
-                ADJOINT_ATOMIC_ADD_GPU(&dGradient[ls*_S_q+_rs_q],          _A_q*_n00_q+_B_q*_n01_q+_C_q*_n10_q+_D_q*_n11_q); \
-                ADJOINT_ATOMIC_ADD_GPU(&dGradient[ls*_S_q+_rs_q+1],       -_B_q*_n00_q+_A_q*_n01_q-_D_q*_n10_q+_C_q*_n11_q); \
-                ADJOINT_ATOMIC_ADD_GPU(&dGradient[(ls+1)*_S_q+_rs_q],     -_C_q*_n00_q-_D_q*_n01_q+_A_q*_n10_q+_B_q*_n11_q); \
-                ADJOINT_ATOMIC_ADD_GPU(&dGradient[(ls+1)*_S_q+_rs_q+1],    _D_q*_n00_q-_C_q*_n01_q-_B_q*_n10_q+_A_q*_n11_q); \
+            } else if (_ri_q > (REAL)0) { \
+                adjointPairleaderComplexStep(dGradient, _S_q, ls, _rs_q, t, \
+                                              _lr_q, _li_q, _ec_q, _es_q, _cI_q, _sI_q, _ea_q, _ri_q, \
+                                              sEvalR, sExpat, sCosbt, sSinbt, sOpRow0, sOpRow1); \
                 _rs_q += 2; \
+            } else { \
+                _rs_q++; \
             } \
         } \
     }