From 2738454e1d50e22856acf7a056e36cd4675d4ab3 Mon Sep 17 00:00:00 2001 From: fil-monti Date: Thu, 9 Jul 2026 22:10:29 -0700 Subject: [PATCH 01/15] CPU plugins: hide default symbol visibility, export only plugin_init BEAGLE dlopen()s each plugin .so and resolves plugin_init via dlsym; every other symbol is an implementation detail. Without -fvisibility=hidden, identically-named symbols from a shared header-only template base class compiled into multiple co-loaded plugin .so's can interpose on each other at dynamic-link time, corrupting virtual dispatch across plugins loaded in the same process (observed as a SIGSEGV in clSetKernelArg when the OpenCL and OpenCL-Spectral plugins were both loaded). Compile hmsbeagle-cpu/-cpu-sse/-cpu-spectral with hidden default visibility and explicitly re-export plugin_init as the one symbol dlsym needs. --- libhmsbeagle/CPU/BeagleCPUPlugin.cpp | 2 +- libhmsbeagle/CPU/BeagleCPUSSEPlugin.cpp | 2 +- libhmsbeagle/CPU/BeagleCPUSpectralPlugin.cpp | 2 +- libhmsbeagle/CPU/CMakeLists.txt | 6 ++++++ 4 files changed, 9 insertions(+), 3 deletions(-) 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) From 2dd99a2911a3992bd326abb9869f1b045575d8c2 Mon Sep 17 00:00:00 2001 From: fil-monti Date: Thu, 9 Jul 2026 22:10:40 -0700 Subject: [PATCH 02/15] CUDA plugin: hide default symbol visibility, export only plugin_init Same cross-plugin symbol-interposition fix as the CPU plugins: hide all symbols in hmsbeagle-cuda by default and explicitly re-export plugin_init, the only symbol BEAGLE's dlopen()/dlsym() loader actually needs. --- libhmsbeagle/GPU/CMake_CUDA/CMakeLists.txt | 2 ++ libhmsbeagle/GPU/CUDAPlugin.cpp | 2 +- 2 files changed, 3 insertions(+), 1 deletion(-) 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/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(); } } From 26d7844f2e4650a7822df2c696d3ff44d28bb201 Mon Sep 17 00:00:00 2001 From: fil-monti Date: Thu, 9 Jul 2026 22:10:52 -0700 Subject: [PATCH 03/15] CUDA-Spectral plugin: hide default symbol visibility, export only plugin_init Same cross-plugin symbol-interposition fix as the CPU plugins: hide all symbols in hmsbeagle-cuda-spectral by default and explicitly re-export plugin_init, the only symbol BEAGLE's dlopen()/dlsym() loader actually needs. --- libhmsbeagle/GPU/CMake_CUDASpectral/CMakeLists.txt | 2 ++ libhmsbeagle/GPU/CUDASpectralPlugin.cpp | 2 +- 2 files changed, 3 insertions(+), 1 deletion(-) 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/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(); } } From fea8497c342b4427483f2ec977e56451876e2941 Mon Sep 17 00:00:00 2001 From: fil-monti Date: Thu, 9 Jul 2026 22:11:03 -0700 Subject: [PATCH 04/15] OpenCL plugin: hide default symbol visibility, export only plugin_init Same cross-plugin symbol-interposition fix as the CPU plugins: hide all symbols in hmsbeagle-opencl by default and explicitly re-export plugin_init, the only symbol BEAGLE's dlopen()/dlsym() loader actually needs. This plugin is one of the two directly implicated in the original SIGSEGV (co-loaded with OpenCL-Spectral, whose shared template base class was interposing across the two .so's before this fix). --- libhmsbeagle/GPU/CMake_OpenCL/CMakeLists.txt | 2 ++ libhmsbeagle/GPU/OpenCLPlugin.cpp | 2 +- 2 files changed, 3 insertions(+), 1 deletion(-) 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/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(); } } From 8f07e0f79ce6f78b6c365608169237837ddb05fa Mon Sep 17 00:00:00 2001 From: fil-monti Date: Thu, 9 Jul 2026 22:11:15 -0700 Subject: [PATCH 05/15] OpenCL-Spectral plugin: hide default symbol visibility, export only plugin_init Same cross-plugin symbol-interposition fix as the CPU plugins: hide all symbols in hmsbeagle-opencl-spectral by default and explicitly re-export plugin_init, the only symbol BEAGLE's dlopen()/dlsym() loader actually needs. This plugin is the other one directly implicated in the original SIGSEGV in clSetKernelArg (co-loaded with plain OpenCL; a shared header-only template base class was interposing across the two .so's before this fix, corrupting virtual dispatch). --- libhmsbeagle/GPU/CMake_OpenCLSpectral/CMakeLists.txt | 2 ++ libhmsbeagle/GPU/OpenCLSpectralPlugin.cpp | 2 +- 2 files changed, 3 insertions(+), 1 deletion(-) 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/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(); } } From 27c46226792813eab754a9fb5e048251ed5c7338 Mon Sep 17 00:00:00 2001 From: fil-monti Date: Thu, 9 Jul 2026 22:12:43 -0700 Subject: [PATCH 06/15] Spectral GPU impl: align pooled device buffer strides to fix CL_MISALIGNED_SUB_BUFFER_OFFSET dSpectralDistancesOrigin and dEvecTOrigin/dIevcTOrigin are single large device allocations sliced into per-matrix/per-eigendecomposition sub-buffers via CreateSubPointer. AMD's OpenCL implementation enforces sub-buffer offset alignment strictly (CL_MISALIGNED_SUB_BUFFER_OFFSET on GPU instance creation); NVIDIA/Apple's do not, which is why this was invisible until tested on AMD hardware. Round both strides up with AlignMemOffset() so every sub-buffer offset is a multiple of the device's required alignment. Since the aligned stride can now be wider than kCategoryCount elements, record it in the new kSpectralDistanceStrideElements member and use that (not kCategoryCount) in the adjoint offset-queue's per-matrix index computation, which walks the same pooled buffer. --- libhmsbeagle/GPU/BeagleGPUSpectralImpl.h | 6 ++++++ libhmsbeagle/GPU/BeagleGPUSpectralImpl.hpp | 7 ++++--- 2 files changed, 10 insertions(+), 3 deletions(-) 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 { 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..1d7bac35 100644 --- a/libhmsbeagle/GPU/BeagleGPUSpectralImpl.hpp +++ b/libhmsbeagle/GPU/BeagleGPUSpectralImpl.hpp @@ -88,7 +88,8 @@ int BeagleGPUSpectralImpl::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::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); @@ -599,7 +600,7 @@ int BeagleGPUSpectralImpl::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; } From 73e7ddc8b11ec846166c9b4e0c897786aa3e784f Mon Sep 17 00:00:00 2001 From: fil-monti Date: Thu, 9 Jul 2026 22:14:00 -0700 Subject: [PATCH 07/15] GPU spectral: widen eigenvalue buffer for SPECTRAL_REPRESENTATION, not just EIGEN_COMPLEX Spectral kernels (kernelsSpectralIfDef*.cu / kernelsSpectral*.cu) always read a real+imaginary pair per eigenstate, regardless of whether the model actually has complex eigenvalues: SPECTRAL_EIGENVALS_GPU unconditionally indexes eigenValues[PADDED_STATE_COUNT + state]. But BeagleGPUImpl::createInstance only sized/copied the eigenvalue device buffer as complex (2x width) when BEAGLE_FLAG_EIGEN_COMPLEX was set, so for real-eigenvalue models under the spectral representation, that kernel read ran past the data setEigenDecomposition actually wrote, pulling in uninitialized device memory and producing wrong postorder/preorder partials (~1e-2 to 1e-1 error vs. CPU). Fix: widen kEigenValuesSize whenever either BEAGLE_FLAG_EIGEN_COMPLEX or BEAGLE_FLAG_SPECTRAL_REPRESENTATION is set. This also required copying BEAGLE_FLAG_SPECTRAL_REPRESENTATION into kFlags in createInstance, which rebuilds kFlags from scratch and was dropping that bit entirely (everything downstream keys off kFlags, not the caller's raw preference/requirement flags). BeagleGPUSpectralImpl::setEigenDecomposition recomputes this same sizing locally (kEigenValuesSize itself is private to BeagleGPUImpl) and needed the identical widening to stay in sync. --- libhmsbeagle/GPU/BeagleGPUImpl.hpp | 15 ++++++++++++++- libhmsbeagle/GPU/BeagleGPUSpectralImpl.hpp | 7 ++++++- 2 files changed, 20 insertions(+), 2 deletions(-) diff --git a/libhmsbeagle/GPU/BeagleGPUImpl.hpp b/libhmsbeagle/GPU/BeagleGPUImpl.hpp index 446a8b43..f0165ac4 100644 --- a/libhmsbeagle/GPU/BeagleGPUImpl.hpp +++ b/libhmsbeagle/GPU/BeagleGPUImpl.hpp @@ -462,6 +462,14 @@ int BeagleGPUImpl::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::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; diff --git a/libhmsbeagle/GPU/BeagleGPUSpectralImpl.hpp b/libhmsbeagle/GPU/BeagleGPUSpectralImpl.hpp index 1d7bac35..f09ddb4f 100644 --- a/libhmsbeagle/GPU/BeagleGPUSpectralImpl.hpp +++ b/libhmsbeagle/GPU/BeagleGPUSpectralImpl.hpp @@ -216,7 +216,12 @@ int BeagleGPUSpectralImpl::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] From a96a75bfd0cd736abdf0c12db0d457f89cd600f5 Mon Sep 17 00:00:00 2001 From: fil-monti Date: Thu, 9 Jul 2026 22:14:36 -0700 Subject: [PATCH 08/15] Add BEAGLE_DEBUG_EIGEN env-gated dump of eigendecomposition buffers (CPU+GPU) Prints eigenvalues/eigenvectors/inverse-eigenvectors as they're set on both the CPU spectral path (EigenDecompositionSpectral::setEigenDecomposition) and the GPU spectral path (BeagleGPUSpectralImpl::setEigenDecomposition), gated behind the BEAGLE_DEBUG_EIGEN environment variable so it's silent by default. Used to compare CPU vs. GPU eigendecomposition inputs side by side when diagnosing spectral GPU numerical bugs; kept as a permanent debugging aid for the same purpose going forward. --- libhmsbeagle/CPU/EigenDecompositionSpectral.hpp | 7 +++++++ libhmsbeagle/GPU/BeagleGPUSpectralImpl.hpp | 7 +++++++ 2 files changed, 14 insertions(+) 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::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>( eigenValuesStorage[eigenIndex].data(), kStateCount, !isComplex); } diff --git a/libhmsbeagle/GPU/BeagleGPUSpectralImpl.hpp b/libhmsbeagle/GPU/BeagleGPUSpectralImpl.hpp index f09ddb4f..e7015177 100644 --- a/libhmsbeagle/GPU/BeagleGPUSpectralImpl.hpp +++ b/libhmsbeagle/GPU/BeagleGPUSpectralImpl.hpp @@ -286,6 +286,13 @@ int BeagleGPUSpectralImpl::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;igpu; gpuIf->MemcpyHostToDevice(this->dEvec[eigenIndex], Evec, sizeof(Real) * SS); gpuIf->MemcpyHostToDevice(this->dIevc[eigenIndex], Ievc, sizeof(Real) * SS); From 29f49b4c644ec68e3e2375ef43a64f5deb8b1c72 Mon Sep 17 00:00:00 2001 From: fil-monti Date: Thu, 9 Jul 2026 22:14:51 -0700 Subject: [PATCH 09/15] Add BEAGLE_DEBUG_KERNEL_ARGS env-gated kernel-launch tracing LaunchKernel/LaunchKernelConcurrent now print, when BEAGLE_DEBUG_KERNEL_ARGS is set: the kernel name (via clGetKernelInfo), every pointer/int argument, block/grid/work-group dimensions, and an explicit clFinish() with before/after prints bracketing the enqueue. This was the key tool for isolating which specific kernel a GPU hang was in: since clEnqueueNDRangeKernel can succeed while clFinish() never returns, the last kernel that doesn't print "clFinish returned OK" is the one that hung. Silent by default; kept as a permanent debugging aid. --- libhmsbeagle/GPU/GPUInterfaceOpenCL.cpp | 71 +++++++++++++++++++++++++ 1 file changed, 71 insertions(+) diff --git a/libhmsbeagle/GPU/GPUInterfaceOpenCL.cpp b/libhmsbeagle/GPU/GPUInterfaceOpenCL.cpp index 62a4fd34..3d774a72 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), ¶m)); } 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), ¶m)); } @@ -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,11 @@ void GPUInterface::LaunchKernel(GPUFunction deviceFunction, printf("local = %lu\n\n", local); #endif + 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 +618,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 +643,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), ¶m)); } 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), ¶m)); } @@ -632,6 +683,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 +708,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 +752,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 From 92511d388623e4ba9c680d715595c62e98924ab9 Mon Sep 17 00:00:00 2001 From: fil-monti Date: Thu, 9 Jul 2026 22:15:11 -0700 Subject: [PATCH 10/15] Add unconditional dispatchPruneSS trace prints (temporary, not env-gated) Prints "[DISPATCH] BASE" / "[DISPATCH] SPECTRAL" on every postorder pruning dispatch, unconditionally (unlike the other debug instrumentation in this branch, this one isn't gated behind an env var). Used to empirically confirm which BeagleGPUImpl vs. BeagleGPUSpectralImpl dispatch path a given plugin/instance was actually taking while diagnosing the cross-plugin symbol-interposition SIGSEGV. Noisy on every pruning call; should be removed or gated behind an env var (e.g. BEAGLE_DEBUG_FLOW, already used elsewhere in this file) before this branch is considered done. --- libhmsbeagle/GPU/BeagleGPUImpl.hpp | 1 + libhmsbeagle/GPU/BeagleGPUSpectralImpl.hpp | 1 + 2 files changed, 2 insertions(+) diff --git a/libhmsbeagle/GPU/BeagleGPUImpl.hpp b/libhmsbeagle/GPU/BeagleGPUImpl.hpp index f0165ac4..234a8fa6 100644 --- a/libhmsbeagle/GPU/BeagleGPUImpl.hpp +++ b/libhmsbeagle/GPU/BeagleGPUImpl.hpp @@ -4756,6 +4756,7 @@ void BeagleGPUImpl::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.hpp b/libhmsbeagle/GPU/BeagleGPUSpectralImpl.hpp index e7015177..06d8b14e 100644 --- a/libhmsbeagle/GPU/BeagleGPUSpectralImpl.hpp +++ b/libhmsbeagle/GPU/BeagleGPUSpectralImpl.hpp @@ -192,6 +192,7 @@ void BeagleGPUSpectralImpl::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, From 512e7cc87ac59d94c99998e3d7a7e477e1d1bf63 Mon Sep 17 00:00:00 2001 From: fil-monti Date: Thu, 9 Jul 2026 22:18:40 -0700 Subject: [PATCH 11/15] GPU adjoint kernel: fix brace mismatch, work around ROCm miscompilation of heavy inlined code MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Three related fixes to the GPU adjoint gradient kernel machinery (kernelAdjointMergedN and the SINGLETON/PAIRLEADER macros it dispatches to), all addressing the same underlying defect class: this ROCm 4.0.1/gfx906 setup's LLVM-AMDGPU backend miscompiles large, register-heavy inlined code when the kernel is built at PADDED_STATE_COUNT=32 (any state count that pads above 16, e.g. 17 states) — confirmed not to be a BEAGLE logic bug, since the identical macros at PADDED_STATE_COUNT=16 work correctly. 1. ADJOINT_ATOMIC_ADD_GPU's CAS-retry loop (float/double atomic add) was inlined directly at every call site. At PADDED_STATE_COUNT=32 this caused clEnqueueNDRangeKernel to succeed while clFinish() never returned. Bisected to require both the atomic_cmpxchg call AND a retry loop reading its own result. Fixed by moving the loop into its own __attribute__((noinline)) helper function per precision. 2. Restructuring ADJOINTN_APPLY_COMPLEX_PAIRLEADER's if/else into if/else-if (to add an explicit `_ri_q < 0` guard) had dropped the macro's final closing brace for its outer `if (tid == 0) { ... }` wrapper, breaking OpenCL kernel compilation entirely ("expected '}'"). Verified by a running brace-balance count against the sibling macro ADJOINTN_APPLY_COMPLEX_SINGLETON (untouched, correct) as ground truth: SINGLETON was 6-open/6-close balanced, PAIRLEADER was missing one close. Fixed by restoring the missing brace. 3. Once compiling and no longer hanging, 17-state gradients still showed NaN. Extracted the per-iteration bodies of both SINGLETON and PAIRLEADER (each carrying ~30 REAL temporaries) into noinline helper functions, the same treatment as (1), on the theory that inlining this much register pressure into the loop was the miscompilation trigger. This did NOT fix the residual NaN (root cause turned out to be deeper — see the GPUInterfaceOpenCL.cpp diagnostic addition and CLUSTER_AGENT_FINDINGS.md for the ongoing investigation), but is kept in place as a harmless, consistent mitigation for this defect class. --- .../GPU/kernels/kernelsSpectralIfDef.cu | 251 ++++++++++++------ 1 file changed, 165 insertions(+), 86 deletions(-) diff --git a/libhmsbeagle/GPU/kernels/kernelsSpectralIfDef.cu b/libhmsbeagle/GPU/kernels/kernelsSpectralIfDef.cu index 7b359fcf..9285d86c 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 */ @@ -983,36 +1002,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 +1140,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++; \ } \ } \ } From 181c09a891803e2e776b12f594944ab5729e7cbb Mon Sep 17 00:00:00 2001 From: fil-monti Date: Thu, 9 Jul 2026 22:19:16 -0700 Subject: [PATCH 12/15] GPU spectral: fix heap buffer overflow in calculateAdjointCrossProducts MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The device gradient buffer is laid out S×S with S=kPaddedStateCount (32 for a 17-state model), but every caller (matching BeagleCPUImpl's convention, confirmed against the real downstream consumer, SpectralBeagleCrossProductDelegate.java, which allocates stateCount*stateCount) sizes outSumDerivatives as kStateCount*kStateCount (289 for 17 states) — unpadded. The old code downloaded and copied all S*S = 1024 values into that 289-element buffer verbatim: a ~5.9 KB heap overflow whenever kPaddedStateCount > kStateCount, i.e. whenever the state count isn't already a "natural" width (4, 16, ...). This was never hit before because the padded (e.g. 17-state) adjoint path always hung or produced garbage earlier in the pipeline before reaching this code. Fixed by copying only the top-left kStateCount×kStateCount submatrix, using the correct stride on each side (S on the device-buffer source, kStateCount on the destination). --- libhmsbeagle/GPU/BeagleGPUSpectralImpl.hpp | 17 ++++++++++++++--- 1 file changed, 14 insertions(+), 3 deletions(-) diff --git a/libhmsbeagle/GPU/BeagleGPUSpectralImpl.hpp b/libhmsbeagle/GPU/BeagleGPUSpectralImpl.hpp index 06d8b14e..10e7a418 100644 --- a/libhmsbeagle/GPU/BeagleGPUSpectralImpl.hpp +++ b/libhmsbeagle/GPU/BeagleGPUSpectralImpl.hpp @@ -639,13 +639,24 @@ int BeagleGPUSpectralImpl::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); } From d1053d79195bcfe69219946278031e031cbf9828 Mon Sep 17 00:00:00 2001 From: fil-monti Date: Thu, 9 Jul 2026 22:19:53 -0700 Subject: [PATCH 13/15] adjointtest4: remove padded-buffer workaround now that the library fixes it MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit This test previously allocated a kPaddedStateCount-sized gradient buffer and re-indexed down to kStateCount×kStateCount, compensating from the test side for calculateAdjointCrossProducts overflowing an unpadded buffer. Now that the library itself only writes the correct kStateCount² submatrix (see the BeagleGPUSpectralImpl.hpp fix), this workaround is redundant and was actually incorrect against the real calling convention. Removed it so the test allocates and reads exactly like the real beast-mcmc caller (SpectralBeagleCrossProductDelegate.java) and BeagleCPUImpl both do. --- examples/hmctest/adjointtest4.cpp | 16 +++++----------- 1 file changed, 5 insertions(+), 11 deletions(-) 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 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 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 gbuf(Sp * Sp, 0.0); + std::vector 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"); } } From 65b35383ee8326a15cbc379ea85e857a9e07cdcf Mon Sep 17 00:00:00 2001 From: fil-monti Date: Thu, 9 Jul 2026 22:20:27 -0700 Subject: [PATCH 14/15] Add CL_KERNEL_WORK_GROUP_SIZE reporting to BEAGLE_DEBUG_KERNEL_ARGS tracing Reports the driver-queried CL_KERNEL_WORK_GROUP_SIZE and CL_KERNEL_PREFERRED_WORK_GROUP_SIZE_MULTIPLE for each kernel launch alongside the requested local work size, under the existing BEAGLE_DEBUG_KERNEL_ARGS gate. Added while investigating whether a GPU adjoint kernel's cooperative reduction was silently running with a smaller workgroup than requested; the driver's self-reported capacity ruled that out as the sole explanation but this remains a useful diagnostic for whoever picks up that investigation next (see CLUSTER_AGENT_FINDINGS.md, "Residual 17-state NaN"). --- libhmsbeagle/GPU/GPUInterfaceOpenCL.cpp | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/libhmsbeagle/GPU/GPUInterfaceOpenCL.cpp b/libhmsbeagle/GPU/GPUInterfaceOpenCL.cpp index 3d774a72..f26f11e0 100644 --- a/libhmsbeagle/GPU/GPUInterfaceOpenCL.cpp +++ b/libhmsbeagle/GPU/GPUInterfaceOpenCL.cpp @@ -602,6 +602,18 @@ 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); From 8a3e5b6a2a8d1487a255bdebd67e7ba65493a6b4 Mon Sep 17 00:00:00 2001 From: fil-monti Date: Thu, 9 Jul 2026 23:31:28 -0700 Subject: [PATCH 15/15] GPU adjoint kernel: fix get_local_id(0) corruption by extracting the shared reduction MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Fixes the residual NaN gradients for padded (PADDED_STATE_COUNT=32, e.g. 17-state) models that remained after the brace-mismatch and buffer-overflow fixes: get_local_id(0) was reading as a stuck constant for every work-item in this compiled kernel variant, a ROCm 4.0.1/LLVM-AMDGPU backend miscompilation (not a BEAGLE logic bug — confirmed via direct atomic-counter probes in a prior session, see CLUSTER_AGENT_FINDINGS.md §3). The already-proven noinline mitigation (used earlier for the ADJOINT_ATOMIC_ADD_GPU CAS-retry loop) had been tried on the SINGLETON/PAIRLEADER macro bodies without effect, which ruled those out but left the true hot spot unidentified. ADJOINTN_ACCUM_ROW_RAW and ADJOINTN_ROTATE_ROW were never tested: they're the shared prefix code executed on every single code path (isAllReal, singleton, and pairleader alike), and contain the kernel's densest concentration of tid-predicated, barrier-synchronized code (a 32x-repeated, 128-wide tree reduction in ADJOINTN_ACCUM_ROW_RAW alone). Extracting both into their own __attribute__((noinline)) helper functions (adjointAccumRowRaw, adjointRotateRow), giving them their own register-allocation scope, fixes the corruption entirely. Verified: adjointtest4 --nstates 17's gradient comparison goes from NaN/FAIL to PASS (max|diff|=5.58e-01 vs. threshold 22.1); the 16-state and 4-state cases are unaffected (16-state's pre-existing ~1e-3 precision-level diff confirmed unchanged by temporarily reverting just this fix and comparing); and, independently, the full rabies_smoke.xml through beast-mcmc now shows Pr(accept)=0.8 for VanillaHMC(host.logRates), up from 0.0 in the prior (buggy) run — real-world confirmation this was a genuine gradient-correctness bug feeding the HMC leapfrog integrator, not just a synthetic test artifact. Full investigation trail (including a negative-control experiment on regOp's dynamic indexing that turned out to be unable to observe this bug either way, since the corrupted tid breaks the reduction's write-out predicate before regOp's own value could ever surface) is in beagle-bugs/gpu-adjoint-workgroup-corruption/PLAN.md and CLUSTER_AGENT_FINDINGS2.md. --- .../GPU/kernels/kernelsSpectralIfDef.cu | 111 +++++++++++------- 1 file changed, 68 insertions(+), 43 deletions(-) diff --git a/libhmsbeagle/GPU/kernels/kernelsSpectralIfDef.cu b/libhmsbeagle/GPU/kernels/kernelsSpectralIfDef.cu index 9285d86c..96feff1c 100644 --- a/libhmsbeagle/GPU/kernels/kernelsSpectralIfDef.cu +++ b/libhmsbeagle/GPU/kernels/kernelsSpectralIfDef.cu @@ -929,37 +929,56 @@ __attribute__((noinline)) void adjointAtomicAddGpuSPHelper(__global int* _anp, f * 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 @@ -969,20 +988,26 @@ __attribute__((noinline)) void adjointAtomicAddGpuSPHelper(__global int* _anp, f * 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