From f27c8aff35bc8eb324e9b024ed42684b6a676fda Mon Sep 17 00:00:00 2001 From: fil-monti Date: Fri, 10 Jul 2026 10:08:38 -0700 Subject: [PATCH 01/13] Only call clGetKernelInfo when BEAGLE_DEBUG_KERNEL_ARGS is set In both LaunchKernel and LaunchKernelConcurrent, the clGetKernelInfo() call that fetches the kernel name for debug tracing was unconditional, paying for an OpenCL driver round-trip on every single kernel launch even when BEAGLE_DEBUG_KERNEL_ARGS isn't set and the name is never used. Move it inside the existing getenv() gate, matching how every other use of kernelNameBuf in both functions is already gated. Addresses PR #242 review comments (msuchard) on GPUInterfaceOpenCL.cpp:548,659. --- libhmsbeagle/GPU/GPUInterfaceOpenCL.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/libhmsbeagle/GPU/GPUInterfaceOpenCL.cpp b/libhmsbeagle/GPU/GPUInterfaceOpenCL.cpp index f26f11e0..d2a4aada 100644 --- a/libhmsbeagle/GPU/GPUInterfaceOpenCL.cpp +++ b/libhmsbeagle/GPU/GPUInterfaceOpenCL.cpp @@ -545,8 +545,8 @@ void GPUInterface::LaunchKernel(GPUFunction deviceFunction, #endif char kernelNameBuf[256] = {0}; - clGetKernelInfo(deviceFunction, CL_KERNEL_FUNCTION_NAME, sizeof(kernelNameBuf), kernelNameBuf, NULL); if (getenv("BEAGLE_DEBUG_KERNEL_ARGS")) { + clGetKernelInfo(deviceFunction, CL_KERNEL_FUNCTION_NAME, sizeof(kernelNameBuf), kernelNameBuf, NULL); fprintf(stderr, "[LaunchKernel] %s: parameterCountV=%d totalParameterCount=%d\n", kernelNameBuf, parameterCountV, totalParameterCount); } @@ -656,8 +656,8 @@ void GPUInterface::LaunchKernelConcurrent(GPUFunction deviceFunction, #endif char kernelNameBuf[256] = {0}; - clGetKernelInfo(deviceFunction, CL_KERNEL_FUNCTION_NAME, sizeof(kernelNameBuf), kernelNameBuf, NULL); if (getenv("BEAGLE_DEBUG_KERNEL_ARGS")) { + clGetKernelInfo(deviceFunction, CL_KERNEL_FUNCTION_NAME, sizeof(kernelNameBuf), kernelNameBuf, NULL); fprintf(stderr, "[LaunchKernelConcurrent] %s: parameterCountV=%d totalParameterCount=%d streamIndex=%d waitIndex=%d\n", kernelNameBuf, parameterCountV, totalParameterCount, streamIndex, waitIndex); } From 741f8a8b1802f54bcbd945e64e1d1fd9dee68a27 Mon Sep 17 00:00:00 2001 From: fil-monti Date: Fri, 10 Jul 2026 13:05:22 -0700 Subject: [PATCH 02/13] GPU spectral kernels: gate imaginary-eigenvalue reads on isAllReal1/isAllReal2 SPECTRAL_EIGENVALS_GPU() and SPECTRAL_EIGENVALS_SIB_ONLY_GPU() (in both kernelsSpectralIfDef.cu, the general N-state file, and kernelsSpectralIfDef4.cu, the 4-state specialization -- these carry byte-identical copies of both macros) unconditionally read eigenValues[PADDED_STATE_COUNT + state], regardless of whether that eigendecomposition is actually complex, unlike the adjoint kernel (kernelAdjointMergedN), which already gates this read correctly behind a per-branch isAllReal flag. Add isAllReal1/isAllReal2 parameters (one per child) to both macros and every one of the twelve kernel functions that call them, and skip the imaginary read + SPECTRAL_SINCOS call entirely when real, writing sDs=e, sCs=0 directly (the same real-eigenvalue convention already documented elsewhere in this file) rather than relying on cos(0)=1/sin(0)=0 falling out of an angle that would otherwise never need to be read. This is a prerequisite for narrowing the host-side eigenvalue buffer sizing back to EIGEN_COMPLEX-only (next commit) -- without this gating, that narrowing would reopen the original out-of-bounds read. See beagle-bugs/gpu-spectral-imaginary-eigenvalue-gating/PLAN.md for the full analysis (PR #242 review comment from msuchard on BeagleGPUImpl.hpp:513). Mirrors the adjoint kernel's already-proven isAllReal pattern rather than a compile-time kernel-variant split, since real/complex status is a per-eigendecomposition runtime property (recomputed on every setEigenDecomposition call, can change between MCMC iterations) rather than a fixed per-instance one. Note: the CUDA-side kernels (kernelsSpectral.cu/kernelsSpectral4.cu) have the identical unconditional-read issue but are not touched here -- BUILD_CUDA is off on the hardware available for this work, so this couldn't be built or verified there. --- .../GPU/kernels/kernelsSpectralIfDef.cu | 98 +++++++++++++------ .../GPU/kernels/kernelsSpectralIfDef4.cu | 98 +++++++++++++------ 2 files changed, 134 insertions(+), 62 deletions(-) diff --git a/libhmsbeagle/GPU/kernels/kernelsSpectralIfDef.cu b/libhmsbeagle/GPU/kernels/kernelsSpectralIfDef.cu index 96feff1c..dccb322b 100644 --- a/libhmsbeagle/GPU/kernels/kernelsSpectralIfDef.cu +++ b/libhmsbeagle/GPU/kernels/kernelsSpectralIfDef.cu @@ -150,23 +150,41 @@ /* patIdx-0 threads: one thread per eigenstate (state = k). * eigenValues layout: [realEV_0..realEV_{S-1} | imagEV_0..imagEV_{S-1}] * distances[matrix] = branchLength * categoryRate[matrix]. + * ISALLREAL1/ISALLREAL2: per-eigendecomposition runtime flags (from + * hEigenDecompIsAllReal, threaded through as ordinary kernel arguments — see + * kernelAdjointMergedN's isAllReal for the precedent this mirrors). When a + * decomposition is all-real, the imaginary half of eigenValues[] is never + * read: it may not even be allocated (kEigenValuesSize is only widened for + * BEAGLE_FLAG_EIGEN_COMPLEX, not merely for spectral-representation use), so + * reading it unconditionally is an out-of-bounds access, not just wasted + * work. sCs == 0 (the real-eigenvalue convention already documented on + * SPECTRAL_COMMON_SMEM_GPU above) is produced explicitly rather than relying + * on cos(0)=1/sin(0)=0 falling out of an unread angle. * Ends with KW_LOCAL_FENCE so sP*, sScale, sDs*, sCs* are all visible. */ -#define SPECTRAL_EIGENVALS_GPU() \ +#define SPECTRAL_EIGENVALS_GPU(ISALLREAL1, ISALLREAL2) \ if (patIdx == 0) { \ REAL t1 = distances1[matrix]; \ REAL e1 = exp(eigenValues1[state] * t1); \ - REAL bt1 = eigenValues1[PADDED_STATE_COUNT + state] * t1; \ - REAL cv1, sv1; \ - SPECTRAL_SINCOS(bt1, sv1, cv1); \ - sDs1[state] = e1 * cv1; \ - sCs1[state] = e1 * sv1; \ + sDs1[state] = e1; \ + sCs1[state] = (REAL)0; \ + if (!(ISALLREAL1)) { \ + REAL bt1 = eigenValues1[PADDED_STATE_COUNT + state] * t1; \ + REAL cv1, sv1; \ + SPECTRAL_SINCOS(bt1, sv1, cv1); \ + sDs1[state] = e1 * cv1; \ + sCs1[state] = e1 * sv1; \ + } \ REAL t2 = distances2[matrix]; \ REAL e2 = exp(eigenValues2[state] * t2); \ - REAL bt2 = eigenValues2[PADDED_STATE_COUNT + state] * t2; \ - REAL cv2, sv2; \ - SPECTRAL_SINCOS(bt2, sv2, cv2); \ - sDs2[state] = e2 * cv2; \ - sCs2[state] = e2 * sv2; \ + sDs2[state] = e2; \ + sCs2[state] = (REAL)0; \ + if (!(ISALLREAL2)) { \ + REAL bt2 = eigenValues2[PADDED_STATE_COUNT + state] * t2; \ + REAL cv2, sv2; \ + SPECTRAL_SINCOS(bt2, sv2, cv2); \ + sDs2[state] = e2 * cv2; \ + sCs2[state] = e2 * sv2; \ + } \ } \ KW_LOCAL_FENCE; @@ -370,6 +388,7 @@ KW_DEVICE_FUNC void kernelSpectralBody( KW_GLOBAL_VAR REAL* KW_RESTRICT eigenValues2, KW_GLOBAL_VAR REAL* KW_RESTRICT distances2, KW_GLOBAL_VAR REAL* KW_RESTRICT scalingFactors, /* NULL if !USE_SCALING */ + int isAllReal1, int isAllReal2, int totalPatterns) { SPECTRAL_INDICES_GPU() @@ -387,7 +406,7 @@ KW_DEVICE_FUNC void kernelSpectralBody( SPECTRAL_LOAD_SCALE_GPU() #endif - SPECTRAL_EIGENVALS_GPU() /* ends with KW_LOCAL_FENCE */ + SPECTRAL_EIGENVALS_GPU(isAllReal1, isAllReal2) /* ends with KW_LOCAL_FENCE */ /* Phase 1: project to eigenspace. * PP case fuses both children into one peel loop (half the barriers). @@ -443,12 +462,13 @@ KW_GLOBAL_KERNEL void kernelPartialsPartialsNoScaleSpectral( KW_GLOBAL_VAR REAL* KW_RESTRICT evec2, KW_GLOBAL_VAR REAL* KW_RESTRICT eigenValues2, KW_GLOBAL_VAR REAL* KW_RESTRICT distances2, + int isAllReal1, int isAllReal2, int totalPatterns) { SPECTRAL_INDICES_GPU() SPECTRAL_COMMON_SMEM_GPU() SPECTRAL_LOAD_PARTIALS1_GPU() SPECTRAL_LOAD_PARTIALS2_GPU() - SPECTRAL_EIGENVALS_GPU() + SPECTRAL_EIGENVALS_GPU(isAllReal1, isAllReal2) SPECTRAL_PHASE1_PARTIALS_DUAL_GPU(sBuf1, sP1, sQ1, ievc1, sBuf2, sP2, sQ2, ievc2) SPECTRAL_PHASE2_GPU() SPECTRAL_PHASE3_GPU() @@ -468,13 +488,14 @@ KW_GLOBAL_KERNEL void kernelPartialsPartialsFixedScaleSpectral( KW_GLOBAL_VAR REAL* KW_RESTRICT eigenValues2, KW_GLOBAL_VAR REAL* KW_RESTRICT distances2, KW_GLOBAL_VAR REAL* KW_RESTRICT scalingFactors, + int isAllReal1, int isAllReal2, int totalPatterns) { SPECTRAL_INDICES_GPU() SPECTRAL_COMMON_SMEM_GPU() SPECTRAL_LOAD_PARTIALS1_GPU() SPECTRAL_LOAD_PARTIALS2_GPU() SPECTRAL_LOAD_SCALE_GPU() - SPECTRAL_EIGENVALS_GPU() + SPECTRAL_EIGENVALS_GPU(isAllReal1, isAllReal2) SPECTRAL_PHASE1_PARTIALS_DUAL_GPU(sBuf1, sP1, sQ1, ievc1, sBuf2, sP2, sQ2, ievc2) SPECTRAL_PHASE2_GPU() SPECTRAL_PHASE3_GPU() @@ -497,11 +518,12 @@ KW_GLOBAL_KERNEL void kernelStatesPartialsNoScaleSpectral( KW_GLOBAL_VAR REAL* KW_RESTRICT evec2, KW_GLOBAL_VAR REAL* KW_RESTRICT eigenValues2, KW_GLOBAL_VAR REAL* KW_RESTRICT distances2, + int isAllReal1, int isAllReal2, int totalPatterns) { SPECTRAL_INDICES_GPU() SPECTRAL_COMMON_SMEM_GPU() SPECTRAL_LOAD_PARTIALS2_GPU() /* no sP1: child 1 is States */ - SPECTRAL_EIGENVALS_GPU() + SPECTRAL_EIGENVALS_GPU(isAllReal1, isAllReal2) SPECTRAL_PHASE1_STATES_GPU(sQ1, ievc1, states1) SPECTRAL_PHASE1_PARTIALS_GPU(sBuf2, sP2, sQ2, ievc2) SPECTRAL_PHASE2_GPU() @@ -522,12 +544,13 @@ KW_GLOBAL_KERNEL void kernelStatesPartialsFixedScaleSpectral( KW_GLOBAL_VAR REAL* KW_RESTRICT eigenValues2, KW_GLOBAL_VAR REAL* KW_RESTRICT distances2, KW_GLOBAL_VAR REAL* KW_RESTRICT scalingFactors, + int isAllReal1, int isAllReal2, int totalPatterns) { SPECTRAL_INDICES_GPU() SPECTRAL_COMMON_SMEM_GPU() SPECTRAL_LOAD_PARTIALS2_GPU() SPECTRAL_LOAD_SCALE_GPU() - SPECTRAL_EIGENVALS_GPU() + SPECTRAL_EIGENVALS_GPU(isAllReal1, isAllReal2) SPECTRAL_PHASE1_STATES_GPU(sQ1, ievc1, states1) SPECTRAL_PHASE1_PARTIALS_GPU(sBuf2, sP2, sQ2, ievc2) SPECTRAL_PHASE2_GPU() @@ -550,10 +573,11 @@ KW_GLOBAL_KERNEL void kernelStatesStatesNoScaleSpectral( KW_GLOBAL_VAR REAL* KW_RESTRICT evec2, KW_GLOBAL_VAR REAL* KW_RESTRICT eigenValues2, KW_GLOBAL_VAR REAL* KW_RESTRICT distances2, + int isAllReal1, int isAllReal2, int totalPatterns) { SPECTRAL_INDICES_GPU() SPECTRAL_COMMON_SMEM_GPU() /* no LOAD_PARTIALS: both children are States */ - SPECTRAL_EIGENVALS_GPU() + SPECTRAL_EIGENVALS_GPU(isAllReal1, isAllReal2) SPECTRAL_PHASE1_STATES_GPU(sQ1, ievc1, states1) SPECTRAL_PHASE1_STATES_GPU(sQ2, ievc2, states2) SPECTRAL_PHASE2_GPU() @@ -574,11 +598,12 @@ KW_GLOBAL_KERNEL void kernelStatesStatesFixedScaleSpectral( KW_GLOBAL_VAR REAL* KW_RESTRICT eigenValues2, KW_GLOBAL_VAR REAL* KW_RESTRICT distances2, KW_GLOBAL_VAR REAL* KW_RESTRICT scalingFactors, + int isAllReal1, int isAllReal2, int totalPatterns) { SPECTRAL_INDICES_GPU() SPECTRAL_COMMON_SMEM_GPU() SPECTRAL_LOAD_SCALE_GPU() - SPECTRAL_EIGENVALS_GPU() + SPECTRAL_EIGENVALS_GPU(isAllReal1, isAllReal2) SPECTRAL_PHASE1_STATES_GPU(sQ1, ievc1, states1) SPECTRAL_PHASE1_STATES_GPU(sQ2, ievc2, states2) SPECTRAL_PHASE2_GPU() @@ -650,16 +675,21 @@ KW_GLOBAL_KERNEL void kernelStatesStatesFixedScaleSpectral( partials3[u] = sum; \ } -/* Load only sibling (child2) eigenvalue exponentials — for Top Root kernels. */ -#define SPECTRAL_EIGENVALS_SIB_ONLY_GPU() \ +/* Load only sibling (child2) eigenvalue exponentials — for Top Root kernels. + * ISALLREAL2: see SPECTRAL_EIGENVALS_GPU's comment above — same gating. */ +#define SPECTRAL_EIGENVALS_SIB_ONLY_GPU(ISALLREAL2) \ if (patIdx == 0) { \ REAL t2 = distances2[matrix]; \ REAL e2 = exp(eigenValues2[state] * t2); \ - REAL bt2 = eigenValues2[PADDED_STATE_COUNT + state] * t2; \ - REAL cv2, sv2; \ - SPECTRAL_SINCOS(bt2, sv2, cv2); \ - sDs2[state] = e2 * cv2; \ - sCs2[state] = e2 * sv2; \ + sDs2[state] = e2; \ + sCs2[state] = (REAL)0; \ + if (!(ISALLREAL2)) { \ + REAL bt2 = eigenValues2[PADDED_STATE_COUNT + state] * t2; \ + REAL cv2, sv2; \ + SPECTRAL_SINCOS(bt2, sv2, cv2); \ + sDs2[state] = e2 * cv2; \ + sCs2[state] = e2 * sv2; \ + } \ } \ KW_LOCAL_FENCE; @@ -694,12 +724,13 @@ KW_GLOBAL_KERNEL void kernelPartialsPartialsAutoScaleSpectral( KW_GLOBAL_VAR REAL* KW_RESTRICT eigenValues2, KW_GLOBAL_VAR REAL* KW_RESTRICT distances2, KW_GLOBAL_VAR signed char* KW_RESTRICT scalingFactors, + int isAllReal1, int isAllReal2, int totalPatterns) { SPECTRAL_INDICES_GPU() SPECTRAL_COMMON_SMEM_GPU() SPECTRAL_LOAD_PARTIALS1_GPU() SPECTRAL_LOAD_PARTIALS2_GPU() - SPECTRAL_EIGENVALS_GPU() + SPECTRAL_EIGENVALS_GPU(isAllReal1, isAllReal2) SPECTRAL_PHASE1_PARTIALS_DUAL_GPU(sBuf1, sP1, sQ1, ievc1, sBuf2, sP2, sQ2, ievc2) SPECTRAL_PHASE2_GPU() SPECTRAL_PHASE3_GPU() @@ -724,12 +755,13 @@ KW_GLOBAL_KERNEL void kernelPartialsPartialsGrowingSpectral( KW_GLOBAL_VAR REAL* KW_RESTRICT evec2, KW_GLOBAL_VAR REAL* KW_RESTRICT eigenValues2, KW_GLOBAL_VAR REAL* KW_RESTRICT distances2, + int isAllReal1, int isAllReal2, int totalPatterns) { SPECTRAL_INDICES_GPU() SPECTRAL_COMMON_SMEM_GPU() SPECTRAL_LOAD_PARTIALS1_GPU() /* sP1 = parent pre-order */ SPECTRAL_LOAD_PARTIALS2_GPU() /* sP2 = sibling post-order */ - SPECTRAL_EIGENVALS_GPU() + SPECTRAL_EIGENVALS_GPU(isAllReal1, isAllReal2) SPECTRAL_PHASE1_PARTIALS_GPU(sBuf2, sP2, sQ2, ievc2) SPECTRAL_PHASE2_SIB_GPU() SPECTRAL_PHASE3_SIB_HADAMARD_GPU() /* sQ2 = P_sib·p_sib ⊙ p_par */ @@ -750,11 +782,12 @@ KW_GLOBAL_KERNEL void kernelPartialsStatesGrowingSpectral( KW_GLOBAL_VAR REAL* KW_RESTRICT evec2, KW_GLOBAL_VAR REAL* KW_RESTRICT eigenValues2, KW_GLOBAL_VAR REAL* KW_RESTRICT distances2, + int isAllReal1, int isAllReal2, int totalPatterns) { SPECTRAL_INDICES_GPU() SPECTRAL_COMMON_SMEM_GPU() SPECTRAL_LOAD_PARTIALS1_GPU() /* sP1 = parent pre-order */ - SPECTRAL_EIGENVALS_GPU() + SPECTRAL_EIGENVALS_GPU(isAllReal1, isAllReal2) SPECTRAL_PHASE1_STATES_GPU(sQ2, ievc2, states2) SPECTRAL_PHASE2_SIB_GPU() SPECTRAL_PHASE3_SIB_HADAMARD_GPU() /* sQ2 = P_sib·e_s ⊙ p_par */ @@ -785,11 +818,12 @@ KW_GLOBAL_KERNEL void kernelPartialsStatesGrowingTopSpectral( KW_GLOBAL_VAR REAL* KW_RESTRICT evec2, KW_GLOBAL_VAR REAL* KW_RESTRICT eigenValues2, KW_GLOBAL_VAR REAL* KW_RESTRICT distances2, + int isAllReal1, int isAllReal2, int totalPatterns) { SPECTRAL_INDICES_GPU() SPECTRAL_COMMON_SMEM_GPU() SPECTRAL_LOAD_PARTIALS1_GPU() - SPECTRAL_EIGENVALS_GPU() + SPECTRAL_EIGENVALS_GPU(isAllReal1, isAllReal2) SPECTRAL_PHASE1_PARTIALS_GPU(sBuf1, sP1, sQ1, ievc1) /* parent backward Ph1 */ SPECTRAL_PHASE1_STATES_GPU(sQ2, ievc2, states2) /* sibling forward Ph1 */ SPECTRAL_PHASE2_GPU() @@ -805,12 +839,13 @@ KW_GLOBAL_KERNEL void kernelPartialsPartialsGrowingTopRootSpectral( KW_GLOBAL_VAR REAL* KW_RESTRICT evec2, KW_GLOBAL_VAR REAL* KW_RESTRICT eigenValues2, KW_GLOBAL_VAR REAL* KW_RESTRICT distances2, + int isAllReal2, int totalPatterns) { SPECTRAL_INDICES_GPU() SPECTRAL_COMMON_SMEM_GPU() SPECTRAL_LOAD_PARTIALS1_GPU() /* sP1 = root pre-order (Hadamard factor) */ SPECTRAL_LOAD_PARTIALS2_GPU() - SPECTRAL_EIGENVALS_SIB_ONLY_GPU() + SPECTRAL_EIGENVALS_SIB_ONLY_GPU(isAllReal2) SPECTRAL_PHASE1_PARTIALS_GPU(sBuf2, sP2, sQ2, ievc2) SPECTRAL_PHASE2_SIB_GPU() SPECTRAL_PHASE3_SIB_HADAMARD_WRITE_GPU() @@ -824,11 +859,12 @@ KW_GLOBAL_KERNEL void kernelPartialsStatesGrowingTopRootSpectral( KW_GLOBAL_VAR REAL* KW_RESTRICT evec2, KW_GLOBAL_VAR REAL* KW_RESTRICT eigenValues2, KW_GLOBAL_VAR REAL* KW_RESTRICT distances2, + int isAllReal2, int totalPatterns) { SPECTRAL_INDICES_GPU() SPECTRAL_COMMON_SMEM_GPU() SPECTRAL_LOAD_PARTIALS1_GPU() /* sP1 = root pre-order (Hadamard factor) */ - SPECTRAL_EIGENVALS_SIB_ONLY_GPU() + SPECTRAL_EIGENVALS_SIB_ONLY_GPU(isAllReal2) SPECTRAL_PHASE1_STATES_GPU(sQ2, ievc2, states2) SPECTRAL_PHASE2_SIB_GPU() SPECTRAL_PHASE3_SIB_HADAMARD_WRITE_GPU() diff --git a/libhmsbeagle/GPU/kernels/kernelsSpectralIfDef4.cu b/libhmsbeagle/GPU/kernels/kernelsSpectralIfDef4.cu index 1ca91d79..af212c3f 100644 --- a/libhmsbeagle/GPU/kernels/kernelsSpectralIfDef4.cu +++ b/libhmsbeagle/GPU/kernels/kernelsSpectralIfDef4.cu @@ -152,23 +152,41 @@ /* patIdx-0 threads: one thread per eigenstate (state = k). * eigenValues layout: [realEV_0..realEV_{S-1} | imagEV_0..imagEV_{S-1}] * distances[matrix] = branchLength * categoryRate[matrix]. + * ISALLREAL1/ISALLREAL2: per-eigendecomposition runtime flags (from + * hEigenDecompIsAllReal, threaded through as ordinary kernel arguments — see + * kernelAdjointMerged4's isAllReal for the precedent this mirrors). When a + * decomposition is all-real, the imaginary half of eigenValues[] is never + * read: it may not even be allocated (kEigenValuesSize is only widened for + * BEAGLE_FLAG_EIGEN_COMPLEX, not merely for spectral-representation use), so + * reading it unconditionally is an out-of-bounds access, not just wasted + * work. sCs == 0 (the real-eigenvalue convention already documented on + * SPECTRAL_COMMON_SMEM_GPU above) is produced explicitly rather than relying + * on cos(0)=1/sin(0)=0 falling out of an unread angle. * Ends with KW_LOCAL_FENCE so sP*, sScale, sDs*, sCs* are all visible. */ -#define SPECTRAL_EIGENVALS_GPU() \ +#define SPECTRAL_EIGENVALS_GPU(ISALLREAL1, ISALLREAL2) \ if (patIdx == 0) { \ REAL t1 = distances1[matrix]; \ REAL e1 = exp(eigenValues1[state] * t1); \ - REAL bt1 = eigenValues1[PADDED_STATE_COUNT + state] * t1; \ - REAL cv1, sv1; \ - SPECTRAL_SINCOS(bt1, sv1, cv1); \ - sDs1[state] = e1 * cv1; \ - sCs1[state] = e1 * sv1; \ + sDs1[state] = e1; \ + sCs1[state] = (REAL)0; \ + if (!(ISALLREAL1)) { \ + REAL bt1 = eigenValues1[PADDED_STATE_COUNT + state] * t1; \ + REAL cv1, sv1; \ + SPECTRAL_SINCOS(bt1, sv1, cv1); \ + sDs1[state] = e1 * cv1; \ + sCs1[state] = e1 * sv1; \ + } \ REAL t2 = distances2[matrix]; \ REAL e2 = exp(eigenValues2[state] * t2); \ - REAL bt2 = eigenValues2[PADDED_STATE_COUNT + state] * t2; \ - REAL cv2, sv2; \ - SPECTRAL_SINCOS(bt2, sv2, cv2); \ - sDs2[state] = e2 * cv2; \ - sCs2[state] = e2 * sv2; \ + sDs2[state] = e2; \ + sCs2[state] = (REAL)0; \ + if (!(ISALLREAL2)) { \ + REAL bt2 = eigenValues2[PADDED_STATE_COUNT + state] * t2; \ + REAL cv2, sv2; \ + SPECTRAL_SINCOS(bt2, sv2, cv2); \ + sDs2[state] = e2 * cv2; \ + sCs2[state] = e2 * sv2; \ + } \ } \ KW_LOCAL_FENCE; @@ -375,6 +393,7 @@ KW_DEVICE_FUNC void kernelSpectralBody( KW_GLOBAL_VAR REAL* KW_RESTRICT eigenValues2, KW_GLOBAL_VAR REAL* KW_RESTRICT distances2, KW_GLOBAL_VAR REAL* KW_RESTRICT scalingFactors, /* NULL if !USE_SCALING */ + int isAllReal1, int isAllReal2, int totalPatterns) { SPECTRAL_INDICES_GPU() @@ -392,7 +411,7 @@ KW_DEVICE_FUNC void kernelSpectralBody( SPECTRAL_LOAD_SCALE_GPU() #endif - SPECTRAL_EIGENVALS_GPU() /* ends with KW_LOCAL_FENCE */ + SPECTRAL_EIGENVALS_GPU(isAllReal1, isAllReal2) /* ends with KW_LOCAL_FENCE */ /* Phase 1: project to eigenspace. * #ifdef replaces if constexpr (is_same) from the C++ version. */ @@ -444,12 +463,13 @@ KW_GLOBAL_KERNEL void kernelPartialsPartialsNoScaleSpectral( KW_GLOBAL_VAR REAL* KW_RESTRICT evec2, KW_GLOBAL_VAR REAL* KW_RESTRICT eigenValues2, KW_GLOBAL_VAR REAL* KW_RESTRICT distances2, + int isAllReal1, int isAllReal2, int totalPatterns) { SPECTRAL_INDICES_GPU() SPECTRAL_COMMON_SMEM_GPU() SPECTRAL_LOAD_PARTIALS1_GPU() SPECTRAL_LOAD_PARTIALS2_GPU() - SPECTRAL_EIGENVALS_GPU() + SPECTRAL_EIGENVALS_GPU(isAllReal1, isAllReal2) SPECTRAL_PHASE1_PARTIALS_DUAL_GPU(sBuf1, sP1, sQ1, ievc1, sBuf2, sP2, sQ2, ievc2) SPECTRAL_PHASE2_GPU() SPECTRAL_PHASE3_GPU() @@ -469,13 +489,14 @@ KW_GLOBAL_KERNEL void kernelPartialsPartialsFixedScaleSpectral( KW_GLOBAL_VAR REAL* KW_RESTRICT eigenValues2, KW_GLOBAL_VAR REAL* KW_RESTRICT distances2, KW_GLOBAL_VAR REAL* KW_RESTRICT scalingFactors, + int isAllReal1, int isAllReal2, int totalPatterns) { SPECTRAL_INDICES_GPU() SPECTRAL_COMMON_SMEM_GPU() SPECTRAL_LOAD_PARTIALS1_GPU() SPECTRAL_LOAD_PARTIALS2_GPU() SPECTRAL_LOAD_SCALE_GPU() - SPECTRAL_EIGENVALS_GPU() + SPECTRAL_EIGENVALS_GPU(isAllReal1, isAllReal2) SPECTRAL_PHASE1_PARTIALS_DUAL_GPU(sBuf1, sP1, sQ1, ievc1, sBuf2, sP2, sQ2, ievc2) SPECTRAL_PHASE2_GPU() SPECTRAL_PHASE3_GPU() @@ -498,11 +519,12 @@ KW_GLOBAL_KERNEL void kernelStatesPartialsNoScaleSpectral( KW_GLOBAL_VAR REAL* KW_RESTRICT evec2, KW_GLOBAL_VAR REAL* KW_RESTRICT eigenValues2, KW_GLOBAL_VAR REAL* KW_RESTRICT distances2, + int isAllReal1, int isAllReal2, int totalPatterns) { SPECTRAL_INDICES_GPU() SPECTRAL_COMMON_SMEM_GPU() SPECTRAL_LOAD_PARTIALS2_GPU() /* no sP1: child 1 is States */ - SPECTRAL_EIGENVALS_GPU() + SPECTRAL_EIGENVALS_GPU(isAllReal1, isAllReal2) SPECTRAL_PHASE1_STATES_GPU(sQ1, ievc1, states1) SPECTRAL_PHASE1_PARTIALS_GPU(sBuf2, sP2, sQ2, ievc2) SPECTRAL_PHASE2_GPU() @@ -523,12 +545,13 @@ KW_GLOBAL_KERNEL void kernelStatesPartialsFixedScaleSpectral( KW_GLOBAL_VAR REAL* KW_RESTRICT eigenValues2, KW_GLOBAL_VAR REAL* KW_RESTRICT distances2, KW_GLOBAL_VAR REAL* KW_RESTRICT scalingFactors, + int isAllReal1, int isAllReal2, int totalPatterns) { SPECTRAL_INDICES_GPU() SPECTRAL_COMMON_SMEM_GPU() SPECTRAL_LOAD_PARTIALS2_GPU() SPECTRAL_LOAD_SCALE_GPU() - SPECTRAL_EIGENVALS_GPU() + SPECTRAL_EIGENVALS_GPU(isAllReal1, isAllReal2) SPECTRAL_PHASE1_STATES_GPU(sQ1, ievc1, states1) SPECTRAL_PHASE1_PARTIALS_GPU(sBuf2, sP2, sQ2, ievc2) SPECTRAL_PHASE2_GPU() @@ -551,10 +574,11 @@ KW_GLOBAL_KERNEL void kernelStatesStatesNoScaleSpectral( KW_GLOBAL_VAR REAL* KW_RESTRICT evec2, KW_GLOBAL_VAR REAL* KW_RESTRICT eigenValues2, KW_GLOBAL_VAR REAL* KW_RESTRICT distances2, + int isAllReal1, int isAllReal2, int totalPatterns) { SPECTRAL_INDICES_GPU() SPECTRAL_COMMON_SMEM_GPU() /* no LOAD_PARTIALS: both children are States */ - SPECTRAL_EIGENVALS_GPU() + SPECTRAL_EIGENVALS_GPU(isAllReal1, isAllReal2) SPECTRAL_PHASE1_STATES_GPU(sQ1, ievc1, states1) SPECTRAL_PHASE1_STATES_GPU(sQ2, ievc2, states2) SPECTRAL_PHASE2_GPU() @@ -575,11 +599,12 @@ KW_GLOBAL_KERNEL void kernelStatesStatesFixedScaleSpectral( KW_GLOBAL_VAR REAL* KW_RESTRICT eigenValues2, KW_GLOBAL_VAR REAL* KW_RESTRICT distances2, KW_GLOBAL_VAR REAL* KW_RESTRICT scalingFactors, + int isAllReal1, int isAllReal2, int totalPatterns) { SPECTRAL_INDICES_GPU() SPECTRAL_COMMON_SMEM_GPU() SPECTRAL_LOAD_SCALE_GPU() - SPECTRAL_EIGENVALS_GPU() + SPECTRAL_EIGENVALS_GPU(isAllReal1, isAllReal2) SPECTRAL_PHASE1_STATES_GPU(sQ1, ievc1, states1) SPECTRAL_PHASE1_STATES_GPU(sQ2, ievc2, states2) SPECTRAL_PHASE2_GPU() @@ -664,16 +689,21 @@ KW_GLOBAL_KERNEL void kernelStatesStatesFixedScaleSpectral( } /* Load only sibling (child2) eigenvalue exponentials — for Top Root kernels - * where child1 is the root and carries no branch transform. */ -#define SPECTRAL_EIGENVALS_SIB_ONLY_GPU() \ + * where child1 is the root and carries no branch transform. + * ISALLREAL2: see SPECTRAL_EIGENVALS_GPU's comment above — same gating. */ +#define SPECTRAL_EIGENVALS_SIB_ONLY_GPU(ISALLREAL2) \ if (patIdx == 0) { \ REAL t2 = distances2[matrix]; \ REAL e2 = exp(eigenValues2[state] * t2); \ - REAL bt2 = eigenValues2[PADDED_STATE_COUNT + state] * t2; \ - REAL cv2, sv2; \ - SPECTRAL_SINCOS(bt2, sv2, cv2); \ - sDs2[state] = e2 * cv2; \ - sCs2[state] = e2 * sv2; \ + sDs2[state] = e2; \ + sCs2[state] = (REAL)0; \ + if (!(ISALLREAL2)) { \ + REAL bt2 = eigenValues2[PADDED_STATE_COUNT + state] * t2; \ + REAL cv2, sv2; \ + SPECTRAL_SINCOS(bt2, sv2, cv2); \ + sDs2[state] = e2 * cv2; \ + sCs2[state] = e2 * sv2; \ + } \ } \ KW_LOCAL_FENCE; @@ -708,12 +738,13 @@ KW_GLOBAL_KERNEL void kernelPartialsPartialsAutoScaleSpectral( KW_GLOBAL_VAR REAL* KW_RESTRICT eigenValues2, KW_GLOBAL_VAR REAL* KW_RESTRICT distances2, KW_GLOBAL_VAR signed char* KW_RESTRICT scalingFactors, + int isAllReal1, int isAllReal2, int totalPatterns) { SPECTRAL_INDICES_GPU() SPECTRAL_COMMON_SMEM_GPU() SPECTRAL_LOAD_PARTIALS1_GPU() SPECTRAL_LOAD_PARTIALS2_GPU() - SPECTRAL_EIGENVALS_GPU() + SPECTRAL_EIGENVALS_GPU(isAllReal1, isAllReal2) SPECTRAL_PHASE1_PARTIALS_DUAL_GPU(sBuf1, sP1, sQ1, ievc1, sBuf2, sP2, sQ2, ievc2) SPECTRAL_PHASE2_GPU() SPECTRAL_PHASE3_GPU() @@ -744,12 +775,13 @@ KW_GLOBAL_KERNEL void kernelPartialsPartialsGrowingSpectral( KW_GLOBAL_VAR REAL* KW_RESTRICT evec2, KW_GLOBAL_VAR REAL* KW_RESTRICT eigenValues2, KW_GLOBAL_VAR REAL* KW_RESTRICT distances2, + int isAllReal1, int isAllReal2, int totalPatterns) { SPECTRAL_INDICES_GPU() SPECTRAL_COMMON_SMEM_GPU() SPECTRAL_LOAD_PARTIALS1_GPU() /* sP1 = parent pre-order */ SPECTRAL_LOAD_PARTIALS2_GPU() /* sP2 = sibling post-order */ - SPECTRAL_EIGENVALS_GPU() + SPECTRAL_EIGENVALS_GPU(isAllReal1, isAllReal2) /* Sibling: forward (P_sib · p_sib) */ SPECTRAL_PHASE1_PARTIALS_GPU(sBuf2, sP2, sQ2, ievc2) SPECTRAL_PHASE2_SIB_GPU() @@ -772,11 +804,12 @@ KW_GLOBAL_KERNEL void kernelPartialsStatesGrowingSpectral( KW_GLOBAL_VAR REAL* KW_RESTRICT evec2, KW_GLOBAL_VAR REAL* KW_RESTRICT eigenValues2, KW_GLOBAL_VAR REAL* KW_RESTRICT distances2, + int isAllReal1, int isAllReal2, int totalPatterns) { SPECTRAL_INDICES_GPU() SPECTRAL_COMMON_SMEM_GPU() SPECTRAL_LOAD_PARTIALS1_GPU() /* sP1 = parent pre-order */ - SPECTRAL_EIGENVALS_GPU() + SPECTRAL_EIGENVALS_GPU(isAllReal1, isAllReal2) /* Sibling: forward (P_sib · e_s) using States Phase 1 */ SPECTRAL_PHASE1_STATES_GPU(sQ2, ievc2, states2) SPECTRAL_PHASE2_SIB_GPU() @@ -815,11 +848,12 @@ KW_GLOBAL_KERNEL void kernelPartialsStatesGrowingTopSpectral( KW_GLOBAL_VAR REAL* KW_RESTRICT evec2, KW_GLOBAL_VAR REAL* KW_RESTRICT eigenValues2, KW_GLOBAL_VAR REAL* KW_RESTRICT distances2, + int isAllReal1, int isAllReal2, int totalPatterns) { SPECTRAL_INDICES_GPU() SPECTRAL_COMMON_SMEM_GPU() SPECTRAL_LOAD_PARTIALS1_GPU() /* sP1 = parent pre-order */ - SPECTRAL_EIGENVALS_GPU() + SPECTRAL_EIGENVALS_GPU(isAllReal1, isAllReal2) SPECTRAL_PHASE1_PARTIALS_GPU(sBuf1, sP1, sQ1, ievc1) /* parent backward Ph1 */ SPECTRAL_PHASE1_STATES_GPU(sQ2, ievc2, states2) /* sibling forward Ph1 */ SPECTRAL_PHASE2_GPU() /* scale sQ1 and sQ2 independently */ @@ -837,12 +871,13 @@ KW_GLOBAL_KERNEL void kernelPartialsPartialsGrowingTopRootSpectral( KW_GLOBAL_VAR REAL* KW_RESTRICT evec2, KW_GLOBAL_VAR REAL* KW_RESTRICT eigenValues2, KW_GLOBAL_VAR REAL* KW_RESTRICT distances2, + int isAllReal2, int totalPatterns) { SPECTRAL_INDICES_GPU() SPECTRAL_COMMON_SMEM_GPU() SPECTRAL_LOAD_PARTIALS1_GPU() /* sP1 = root pre-order (Hadamard factor) */ SPECTRAL_LOAD_PARTIALS2_GPU() /* sP2 = sibling post-order */ - SPECTRAL_EIGENVALS_SIB_ONLY_GPU() + SPECTRAL_EIGENVALS_SIB_ONLY_GPU(isAllReal2) SPECTRAL_PHASE1_PARTIALS_GPU(sBuf2, sP2, sQ2, ievc2) SPECTRAL_PHASE2_SIB_GPU() SPECTRAL_PHASE3_SIB_HADAMARD_WRITE_GPU() /* evec2·sQ2 ⊙ sP1 → partials3 */ @@ -857,11 +892,12 @@ KW_GLOBAL_KERNEL void kernelPartialsStatesGrowingTopRootSpectral( KW_GLOBAL_VAR REAL* KW_RESTRICT evec2, KW_GLOBAL_VAR REAL* KW_RESTRICT eigenValues2, KW_GLOBAL_VAR REAL* KW_RESTRICT distances2, + int isAllReal2, int totalPatterns) { SPECTRAL_INDICES_GPU() SPECTRAL_COMMON_SMEM_GPU() SPECTRAL_LOAD_PARTIALS1_GPU() /* sP1 = root pre-order (Hadamard factor) */ - SPECTRAL_EIGENVALS_SIB_ONLY_GPU() + SPECTRAL_EIGENVALS_SIB_ONLY_GPU(isAllReal2) SPECTRAL_PHASE1_STATES_GPU(sQ2, ievc2, states2) SPECTRAL_PHASE2_SIB_GPU() SPECTRAL_PHASE3_SIB_HADAMARD_WRITE_GPU() From 91c1fe816920e23fdb0c258cb0bac9fa758e5847 Mon Sep 17 00:00:00 2001 From: fil-monti Date: Fri, 10 Jul 2026 13:12:23 -0700 Subject: [PATCH 03/13] GPU spectral: thread isAllReal1/isAllReal2 through host dispatch Adds isAllReal1/isAllReal2 int parameters to the eight KernelLauncher wrapper functions covering the pruning and growing/pre-order spectral kernels (PartialsPartialsPruningSpectral, StatesPartialsPruningSpectral, StatesStatesPruningSpectral, PartialsPartialsGrowingSpectral(Top), PartialsStatesGrowingSpectral(Top), and the two GrowingSpectralTopRoot variants which need only isAllReal2, one child), threading them into the new kernel arguments added in the previous commit and bumping parameterCountV/totalParameterCount at each LaunchKernel/LaunchKernelConcurrent call site accordingly. BeagleGPUSpectralImpl.hpp's dispatchPrunePP/SP/SS and dispatchGrowingSpectral compute the actual flags from hEigenDecompIsAllReal[ei1]/[ei2] -- already populated per-eigendecomposition from the real eigenvalue data in setEigenDecomposition -- and pass them through at every call site. A plain scalar kernel argument (rather than a device offset-queue record, as the adjoint kernel's isAllReal uses) is the right mechanism here: the adjoint kernel is a single merged launch covering many branches at once, while these pruning/growing kernels launch one call per edge/pair, so a plain argument looked up host-side immediately before each launch is simpler and sufficient. --- libhmsbeagle/GPU/BeagleGPUSpectralImpl.hpp | 21 ++++++-- libhmsbeagle/GPU/KernelLauncher.cpp | 60 ++++++++++++++++------ libhmsbeagle/GPU/KernelLauncher.h | 24 +++++++-- 3 files changed, 78 insertions(+), 27 deletions(-) diff --git a/libhmsbeagle/GPU/BeagleGPUSpectralImpl.hpp b/libhmsbeagle/GPU/BeagleGPUSpectralImpl.hpp index 10e7a418..f83444c2 100644 --- a/libhmsbeagle/GPU/BeagleGPUSpectralImpl.hpp +++ b/libhmsbeagle/GPU/BeagleGPUSpectralImpl.hpp @@ -165,6 +165,7 @@ void BeagleGPUSpectralImpl::dispatchPrunePP( this->dIevc[ei2], this->dEvec[ei2], this->dEigenValues[ei2], dSpectralDistances[c2MatIdx], scalingFactors, cumulativeScaling, this->kPaddedPatternCount, this->kCategoryCount, + hEigenDecompIsAllReal[ei1] ? 1 : 0, hEigenDecompIsAllReal[ei2] ? 1 : 0, rescale, streamIndex, waitIndex); } @@ -182,6 +183,7 @@ void BeagleGPUSpectralImpl::dispatchPruneSP( this->dIevc[ei2], this->dEvec[ei2], this->dEigenValues[ei2], dSpectralDistances[c2MatIdx], scalingFactors, cumulativeScaling, this->kPaddedPatternCount, this->kCategoryCount, + hEigenDecompIsAllReal[ei1] ? 1 : 0, hEigenDecompIsAllReal[ei2] ? 1 : 0, rescale, streamIndex, waitIndex); } @@ -200,6 +202,7 @@ void BeagleGPUSpectralImpl::dispatchPruneSS( this->dIevc[ei2], this->dEvec[ei2], this->dEigenValues[ei2], dSpectralDistances[c2MatIdx], scalingFactors, cumulativeScaling, this->kPaddedPatternCount, this->kCategoryCount, + hEigenDecompIsAllReal[ei1] ? 1 : 0, hEigenDecompIsAllReal[ei2] ? 1 : 0, rescale, streamIndex, waitIndex); } @@ -320,18 +323,21 @@ void BeagleGPUSpectralImpl::dispatchGrowingSpectral( GPUPtr evec2 = this->dEvec[ei2]; GPUPtr eval2 = this->dEigenValues[ei2]; GPUPtr dist2 = dSpectralDistances[c2MatIdx]; + int isAllReal2 = hEigenDecompIsAllReal[ei2] ? 1 : 0; if (isRoot) { if (sibIsStates) { this->kernels->PartialsStatesGrowingSpectralTopRoot( partials1, c2, partials3, ievc2, evec2, eval2, dist2, - this->kPaddedPatternCount, this->kCategoryCount); + this->kPaddedPatternCount, this->kCategoryCount, + isAllReal2); } else { this->kernels->PartialsPartialsGrowingSpectralTopRoot( partials1, c2, partials3, ievc2, evec2, eval2, dist2, - this->kPaddedPatternCount, this->kCategoryCount); + this->kPaddedPatternCount, this->kCategoryCount, + isAllReal2); } } else { int ei1 = hEigenIndexForMatrix[c1MatIdx]; @@ -341,6 +347,7 @@ void BeagleGPUSpectralImpl::dispatchGrowingSpectral( GPUPtr bEvec1 = dIevcT[ei1]; GPUPtr eval1 = this->dEigenValues[ei1]; GPUPtr dist1 = dSpectralDistances[c1MatIdx]; + int isAllReal1 = hEigenDecompIsAllReal[ei1] ? 1 : 0; if (isTop) { if (sibIsStates) { @@ -348,7 +355,8 @@ void BeagleGPUSpectralImpl::dispatchGrowingSpectral( partials1, c2, partials3, bIevc1, bEvec1, eval1, dist1, ievc2, evec2, eval2, dist2, - this->kPaddedPatternCount, this->kCategoryCount); + this->kPaddedPatternCount, this->kCategoryCount, + isAllReal1, isAllReal2); } else { // Top NotRoot PP: reuse no-scale PP pruning kernel with backward matrices this->kernels->PartialsPartialsPruningSpectral( @@ -357,6 +365,7 @@ void BeagleGPUSpectralImpl::dispatchGrowingSpectral( ievc2, evec2, eval2, dist2, nullptr, nullptr, this->kPaddedPatternCount, this->kCategoryCount, + isAllReal1, isAllReal2, -1, -1, -1); } } else { @@ -365,13 +374,15 @@ void BeagleGPUSpectralImpl::dispatchGrowingSpectral( partials1, c2, partials3, bIevc1, bEvec1, eval1, dist1, ievc2, evec2, eval2, dist2, - this->kPaddedPatternCount, this->kCategoryCount); + this->kPaddedPatternCount, this->kCategoryCount, + isAllReal1, isAllReal2); } else { this->kernels->PartialsPartialsGrowingSpectral( partials1, c2, partials3, bIevc1, bEvec1, eval1, dist1, ievc2, evec2, eval2, dist2, - this->kPaddedPatternCount, this->kCategoryCount); + this->kPaddedPatternCount, this->kCategoryCount, + isAllReal1, isAllReal2); } } } diff --git a/libhmsbeagle/GPU/KernelLauncher.cpp b/libhmsbeagle/GPU/KernelLauncher.cpp index f9de48c0..e62bcc49 100644 --- a/libhmsbeagle/GPU/KernelLauncher.cpp +++ b/libhmsbeagle/GPU/KernelLauncher.cpp @@ -1041,13 +1041,16 @@ void KernelLauncher::PartialsPartialsGrowingSpectral(GPUPtr partials1, GPUPtr ievc2, GPUPtr evec2, GPUPtr eigenValues2, GPUPtr distances2, unsigned int patternCount, - unsigned int categoryCount) { + unsigned int categoryCount, + int isAllReal1, + int isAllReal2) { gpu->LaunchKernel(fPartialsPartialsGrowingSpectral, bgSpectralPeelingBlock, bgSpectralPeelingGrid, - 11, 12, + 11, 14, partials1, partials2, partials3, ievc1, evec1, eigenValues1, distances1, ievc2, evec2, eigenValues2, distances2, + isAllReal1, isAllReal2, patternCount); gpu->SynchronizeDevice(); } @@ -1060,13 +1063,16 @@ void KernelLauncher::PartialsStatesGrowingSpectral(GPUPtr partials1, GPUPtr ievc2, GPUPtr evec2, GPUPtr eigenValues2, GPUPtr distances2, unsigned int patternCount, - unsigned int categoryCount) { + unsigned int categoryCount, + int isAllReal1, + int isAllReal2) { gpu->LaunchKernel(fPartialsStatesGrowingSpectral, bgSpectralPeelingBlock, bgSpectralPeelingGrid, - 11, 12, + 11, 14, partials1, states2, partials3, ievc1, evec1, eigenValues1, distances1, ievc2, evec2, eigenValues2, distances2, + isAllReal1, isAllReal2, patternCount); gpu->SynchronizeDevice(); } @@ -1079,13 +1085,16 @@ void KernelLauncher::PartialsStatesGrowingSpectralTop(GPUPtr partials1, GPUPtr ievc2, GPUPtr evec2, GPUPtr eigenValues2, GPUPtr distances2, unsigned int patternCount, - unsigned int categoryCount) { + unsigned int categoryCount, + int isAllReal1, + int isAllReal2) { gpu->LaunchKernel(fPartialsStatesGrowingTopSpectral, bgSpectralPeelingBlock, bgSpectralPeelingGrid, - 11, 12, + 11, 14, partials1, states2, partials3, ievc1, evec1, eigenValues1, distances1, ievc2, evec2, eigenValues2, distances2, + isAllReal1, isAllReal2, patternCount); gpu->SynchronizeDevice(); } @@ -1096,12 +1105,14 @@ void KernelLauncher::PartialsPartialsGrowingSpectralTopRoot(GPUPtr partials1, GPUPtr ievc2, GPUPtr evec2, GPUPtr eigenValues2, GPUPtr distances2, unsigned int patternCount, - unsigned int categoryCount) { + unsigned int categoryCount, + int isAllReal2) { gpu->LaunchKernel(fPartialsPartialsGrowingTopRootSpectral, bgSpectralPeelingBlock, bgSpectralPeelingGrid, - 7, 8, + 7, 9, partials1, partials2, partials3, ievc2, evec2, eigenValues2, distances2, + isAllReal2, patternCount); gpu->SynchronizeDevice(); } @@ -1112,12 +1123,14 @@ void KernelLauncher::PartialsStatesGrowingSpectralTopRoot(GPUPtr partials1, GPUPtr ievc2, GPUPtr evec2, GPUPtr eigenValues2, GPUPtr distances2, unsigned int patternCount, - unsigned int categoryCount) { + unsigned int categoryCount, + int isAllReal2) { gpu->LaunchKernel(fPartialsStatesGrowingTopRootSpectral, bgSpectralPeelingBlock, bgSpectralPeelingGrid, - 7, 8, + 7, 9, partials1, states2, partials3, ievc2, evec2, eigenValues2, distances2, + isAllReal2, patternCount); gpu->SynchronizeDevice(); } @@ -1779,26 +1792,30 @@ void KernelLauncher::PartialsPartialsPruningSpectral(GPUPtr partials1, GPUPtr cumulativeScaling, unsigned int patternCount, unsigned int categoryCount, + int isAllReal1, + int isAllReal2, int doRescaling, int streamIndex, int waitIndex) { if (doRescaling == 2) { // auto-rescaling gpu->LaunchKernel(fPartialsPartialsByPatternBlockAutoScaling, bgSpectralPeelingBlock, bgSpectralPeelingGrid, - 12, 13, + 12, 15, partials1, partials2, partials3, ievc1, evec1, eigenValues1, distances1, ievc2, evec2, eigenValues2, distances2, scalingFactors, + isAllReal1, isAllReal2, patternCount); } else if (doRescaling != 0) { gpu->LaunchKernelConcurrent(fPartialsPartialsByPatternBlockCoherent, bgSpectralPeelingBlock, bgSpectralPeelingGrid, streamIndex, waitIndex, - 11, 12, + 11, 14, partials1, partials2, partials3, ievc1, evec1, eigenValues1, distances1, ievc2, evec2, eigenValues2, distances2, + isAllReal1, isAllReal2, patternCount); if (doRescaling > 0) { KernelLauncher::RescalePartials(partials3, scalingFactors, cumulativeScaling, @@ -1808,11 +1825,12 @@ void KernelLauncher::PartialsPartialsPruningSpectral(GPUPtr partials1, gpu->LaunchKernelConcurrent(fPartialsPartialsByPatternBlockFixedScaling, bgSpectralPeelingBlock, bgSpectralPeelingGrid, streamIndex, waitIndex, - 12, 13, + 12, 15, partials1, partials2, partials3, ievc1, evec1, eigenValues1, distances1, ievc2, evec2, eigenValues2, distances2, scalingFactors, + isAllReal1, isAllReal2, patternCount); } } @@ -1828,6 +1846,8 @@ void KernelLauncher::StatesPartialsPruningSpectral(GPUPtr states1, GPUPtr cumulativeScaling, unsigned int patternCount, unsigned int categoryCount, + int isAllReal1, + int isAllReal2, int doRescaling, int streamIndex, int waitIndex) { @@ -1835,10 +1855,11 @@ void KernelLauncher::StatesPartialsPruningSpectral(GPUPtr states1, gpu->LaunchKernelConcurrent(fStatesPartialsByPatternBlockCoherent, bgSpectralPeelingBlock, bgSpectralPeelingGrid, streamIndex, waitIndex, - 11, 12, + 11, 14, states1, partials2, partials3, ievc1, evec1, eigenValues1, distances1, ievc2, evec2, eigenValues2, distances2, + isAllReal1, isAllReal2, patternCount); if (doRescaling > 0) { KernelLauncher::RescalePartials(partials3, scalingFactors, cumulativeScaling, @@ -1848,11 +1869,12 @@ void KernelLauncher::StatesPartialsPruningSpectral(GPUPtr states1, gpu->LaunchKernelConcurrent(fStatesPartialsByPatternBlockFixedScaling, bgSpectralPeelingBlock, bgSpectralPeelingGrid, streamIndex, waitIndex, - 12, 13, + 12, 15, states1, partials2, partials3, ievc1, evec1, eigenValues1, distances1, ievc2, evec2, eigenValues2, distances2, scalingFactors, + isAllReal1, isAllReal2, patternCount); } } @@ -1868,6 +1890,8 @@ void KernelLauncher::StatesStatesPruningSpectral(GPUPtr states1, GPUPtr cumulativeScaling, unsigned int patternCount, unsigned int categoryCount, + int isAllReal1, + int isAllReal2, int doRescaling, int streamIndex, int waitIndex) { @@ -1875,10 +1899,11 @@ void KernelLauncher::StatesStatesPruningSpectral(GPUPtr states1, gpu->LaunchKernelConcurrent(fStatesStatesByPatternBlockCoherent, bgSpectralPeelingBlock, bgSpectralPeelingGrid, streamIndex, waitIndex, - 11, 12, + 11, 14, states1, states2, partials3, ievc1, evec1, eigenValues1, distances1, ievc2, evec2, eigenValues2, distances2, + isAllReal1, isAllReal2, patternCount); if (doRescaling > 0) { KernelLauncher::RescalePartials(partials3, scalingFactors, cumulativeScaling, @@ -1888,11 +1913,12 @@ void KernelLauncher::StatesStatesPruningSpectral(GPUPtr states1, gpu->LaunchKernelConcurrent(fStatesStatesByPatternBlockFixedScaling, bgSpectralPeelingBlock, bgSpectralPeelingGrid, streamIndex, waitIndex, - 12, 13, + 12, 15, states1, states2, partials3, ievc1, evec1, eigenValues1, distances1, ievc2, evec2, eigenValues2, distances2, scalingFactors, + isAllReal1, isAllReal2, patternCount); } } diff --git a/libhmsbeagle/GPU/KernelLauncher.h b/libhmsbeagle/GPU/KernelLauncher.h index 2c76ccc1..d8af125b 100644 --- a/libhmsbeagle/GPU/KernelLauncher.h +++ b/libhmsbeagle/GPU/KernelLauncher.h @@ -320,6 +320,8 @@ class KernelLauncher { GPUPtr cumulativeScaling, unsigned int patternCount, unsigned int categoryCount, + int isAllReal1, + int isAllReal2, int doRescaling, int streamIndex, int waitIndex); @@ -335,6 +337,8 @@ class KernelLauncher { GPUPtr cumulativeScaling, unsigned int patternCount, unsigned int categoryCount, + int isAllReal1, + int isAllReal2, int doRescaling, int streamIndex, int waitIndex); @@ -350,6 +354,8 @@ class KernelLauncher { GPUPtr cumulativeScaling, unsigned int patternCount, unsigned int categoryCount, + int isAllReal1, + int isAllReal2, int doRescaling, int streamIndex, int waitIndex); @@ -386,7 +392,9 @@ class KernelLauncher { GPUPtr ievc2, GPUPtr evec2, GPUPtr eigenValues2, GPUPtr distances2, unsigned int patternCount, - unsigned int categoryCount); + unsigned int categoryCount, + int isAllReal1, + int isAllReal2); void PartialsStatesGrowingSpectral(GPUPtr partials1, GPUPtr states2, @@ -396,7 +404,9 @@ class KernelLauncher { GPUPtr ievc2, GPUPtr evec2, GPUPtr eigenValues2, GPUPtr distances2, unsigned int patternCount, - unsigned int categoryCount); + unsigned int categoryCount, + int isAllReal1, + int isAllReal2); void PartialsStatesGrowingSpectralTop(GPUPtr partials1, GPUPtr states2, @@ -406,7 +416,9 @@ class KernelLauncher { GPUPtr ievc2, GPUPtr evec2, GPUPtr eigenValues2, GPUPtr distances2, unsigned int patternCount, - unsigned int categoryCount); + unsigned int categoryCount, + int isAllReal1, + int isAllReal2); void PartialsPartialsGrowingSpectralTopRoot(GPUPtr partials1, GPUPtr partials2, @@ -414,7 +426,8 @@ class KernelLauncher { GPUPtr ievc2, GPUPtr evec2, GPUPtr eigenValues2, GPUPtr distances2, unsigned int patternCount, - unsigned int categoryCount); + unsigned int categoryCount, + int isAllReal2); void PartialsStatesGrowingSpectralTopRoot(GPUPtr partials1, GPUPtr states2, @@ -422,7 +435,8 @@ class KernelLauncher { GPUPtr ievc2, GPUPtr evec2, GPUPtr eigenValues2, GPUPtr distances2, unsigned int patternCount, - unsigned int categoryCount); + unsigned int categoryCount, + int isAllReal2); /* Adjoint cross-product launchers — merged single launch covers every * branch in a call at once (grid spans branchCount), isStates/isAllReal From 5b2c92db89962136c3577ab39e367a3b6da8a661 Mon Sep 17 00:00:00 2001 From: fil-monti Date: Fri, 10 Jul 2026 13:12:31 -0700 Subject: [PATCH 04/13] GPU spectral: narrow eigenvalue buffer back to EIGEN_COMPLEX-only, matching CPU Now that the pruning kernels gate their imaginary-eigenvalue reads on a per-decomposition isAllReal flag (previous two commits), the device eigenvalue buffer no longer needs to be widened for every BEAGLE_FLAG_SPECTRAL_REPRESENTATION instance regardless of whether it's actually complex -- only EIGEN_COMPLEX requires the extra width, exactly matching EigenDecompositionSpectral.hpp's CPU implementation (`kEigenValuesSize(isComplex ? kStateCount * 2 : kStateCount)`). Addresses PR #242 review comment (msuchard) on BeagleGPUImpl.hpp:513: "this is not right. if isAllReal then the imaginary parts should never be written or read." Widening on SPECTRAL_REPRESENTATION alone had no CPU analog and wasted memory/copy work for spectral-mode instances that are structurally guaranteed real. This change alone would reopen the original out-of-bounds read bug without the kernel-side gating from the previous commits landing first -- see beagle-bugs/gpu-spectral-imaginary-eigenvalue-gating/PLAN.md for the full analysis of why the two changes are coupled. Verified: adjointtest4 --gpu 1 (4/16/17-state) unchanged from baseline; new --realonly test confirms (via BEAGLE_DEBUG_EIGEN=1) the device buffer is now genuinely narrow, not just coincidentally safe; full rabies_smoke.xml through beast-mcmc unaffected (that model already requests EIGEN_COMPLEX, so it was never relying on the removed over-widening). --- libhmsbeagle/GPU/BeagleGPUImpl.hpp | 15 +++++++++------ libhmsbeagle/GPU/BeagleGPUSpectralImpl.hpp | 12 ++++++------ 2 files changed, 15 insertions(+), 12 deletions(-) diff --git a/libhmsbeagle/GPU/BeagleGPUImpl.hpp b/libhmsbeagle/GPU/BeagleGPUImpl.hpp index 234a8fa6..ab16b838 100644 --- a/libhmsbeagle/GPU/BeagleGPUImpl.hpp +++ b/libhmsbeagle/GPU/BeagleGPUImpl.hpp @@ -505,12 +505,15 @@ int BeagleGPUImpl::createInstance(int tipCount, kPartialsSize = kPaddedPatternCount * kPaddedStateCount * kCategoryCount; kMatrixSize = kPaddedStateCount * kPaddedStateCount; - // 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)) + // Spectral kernels (kernelsSpectralIfDef*.cu / kernelsSpectral*.cu) read a real+imaginary + // pair per eigenstate only when that eigendecomposition is not all-real + // (SPECTRAL_EIGENVALS_GPU is now gated on the isAllReal1/isAllReal2 kernel arguments, threaded + // from hEigenDecompIsAllReal — see BeagleGPUSpectralImpl.hpp::setEigenDecomposition), matching + // the CPU implementation's EigenDecompositionSpectral.hpp, which only widens on + // BEAGLE_FLAG_EIGEN_COMPLEX. Widening for BEAGLE_FLAG_SPECTRAL_REPRESENTATION alone (regardless + // of EIGEN_COMPLEX) has no CPU analog and would waste memory/copy work for spectral-mode + // instances that are structurally guaranteed real, so this only widens for EIGEN_COMPLEX. + if (kFlags & BEAGLE_FLAG_EIGEN_COMPLEX) kEigenValuesSize = 2 * kPaddedStateCount; else kEigenValuesSize = kPaddedStateCount; diff --git a/libhmsbeagle/GPU/BeagleGPUSpectralImpl.hpp b/libhmsbeagle/GPU/BeagleGPUSpectralImpl.hpp index f83444c2..cddc45db 100644 --- a/libhmsbeagle/GPU/BeagleGPUSpectralImpl.hpp +++ b/libhmsbeagle/GPU/BeagleGPUSpectralImpl.hpp @@ -220,12 +220,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. - // 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; + // Spectral kernels read a real+imaginary pair per eigenstate only when that decomposition is + // not all-real (see matching comment/fix in BeagleGPUImpl::createInstance) -- gated at the + // kernel level on isAllReal1/isAllReal2, computed below from inEigenValues and cached in + // hEigenDecompIsAllReal. This local copy only widens for BEAGLE_FLAG_EIGEN_COMPLEX, matching + // CPU's EigenDecompositionSpectral.hpp. + const int eigenValuesSize = (this->kFlags & BEAGLE_FLAG_EIGEN_COMPLEX) ? 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 59ea066e2a8cfcececdbf72db41b558ac3b12000 Mon Sep 17 00:00:00 2001 From: fil-monti Date: Fri, 10 Jul 2026 13:13:16 -0700 Subject: [PATCH 05/13] adjointtest4: add real-only and dynamic real/complex transition tests Two new regression tests for the isAllReal gating/buffer-narrowing change: - runTestRealOnly (--realonly): a symmetric/reversible 17-state circulant (r_fwd == r_bkd, so every eigenvalue is exactly real, not just numerically close to zero) that never requests BEAGLE_FLAG_EIGEN_COMPLEX, at a padded state count (17 -> kPaddedStateCount=32). The device eigenvalue buffer is allocated at exactly kPaddedStateCount width, no imaginary slack at all -- if the kernel's imaginary-eigenvalue read weren't properly gated, this configuration reads past the end of that buffer. This is the one existing-test gap the PR #242 review comment identified: every previous test instance requested EIGEN_COMPLEX, so none of them exercised the narrow-buffer code path this change enables. - runTestDynamicTransition (--dynamic): a single BEAGLE_FLAG_EIGEN_COMPLEX 16-state instance with three sequential setEigenDecomposition calls on the same eigenIndex -- real, then complex, then real again -- exercising hEigenDecompIsAllReal's per-call (not per-instance) nature: it's recomputed from the actual eigenvalues on every call, since a model's eigendecomposition can genuinely flip between real and complex across MCMC iterations. Includes a step0-vs-step2 consistency check (identical real input data, with a complex decomposition processed on the same instance in between) to catch stale kernel/device state from the transition; the small residual diff this surfaced was tracked down to pre-existing GPU atomic-add reduction non-associativity noise between separate kernel launches, unrelated to this change, and documented inline. buildCirculantN() gained optional r_fwd/r_bkd parameters (defaulting to the existing asymmetric 1.0/0.5 values, so all prior callers are unaffected) to support constructing the symmetric/real-only case above without duplicating the eigendecomposition-construction logic. Both new modes are included in --all alongside the existing 4/16/17-state tests. --- examples/hmctest/adjointtest4.cpp | 345 +++++++++++++++++++++++++++++- 1 file changed, 336 insertions(+), 9 deletions(-) diff --git a/examples/hmctest/adjointtest4.cpp b/examples/hmctest/adjointtest4.cpp index c28ee171..55532a6f 100644 --- a/examples/hmctest/adjointtest4.cpp +++ b/examples/hmctest/adjointtest4.cpp @@ -91,13 +91,18 @@ static double edgeLengths[4] = { 0.1, 0.1, 0.2, 0.5 }; * Evec[j, 2m-1] = (1/N)*cos(2πjm/N) Ievc[2m-1, j] = 2*cos(2πjm/N) * Evec[j, 2m ] = (1/N)*sin(2πjm/N) Ievc[2m, j] = 2*sin(2πjm/N) */ +/* r_fwd/r_bkd default to the asymmetric (complex-eigenvalue) case used by the + * existing 16-/17-state tests. Passing r_fwd == r_bkd instead makes Q + * symmetric/reversible: b = (r_fwd-r_bkd)*sin(theta_m) is then identically + * zero for every mode, so every eigenvalue comes out genuinely real (not + * just numerically close to zero) -- used by the real-only test below. */ static void buildCirculantN( int N, std::vector& evec, /* N×N */ std::vector& ivec, /* N×N */ - std::vector& eval) /* 2N */ + std::vector& eval, /* 2N */ + double r_fwd = 1.0, double r_bkd = 0.5) { - const double r_fwd = 1.0, r_bkd = 0.5; const double two_pi_over_N = 2.0 * M_PI / N; const bool evenN = (N % 2 == 0); const int nPairs = evenN ? N/2 - 1 : (N-1)/2; @@ -934,6 +939,318 @@ static int runTest17(int gpuDevice) { return gradOk ? 0 : 1; } +/* ═══════════════════════════════════════════════════════════════════════ + * Real-only 17-state test (BEAGLE_FLAG_EIGEN_COMPLEX NOT requested) + * + * Regression test for the msuchard PR #242 review comment ("if isAllReal + * then the imaginary parts should never be written or read"): the GPU + * spectral pruning kernels now gate their imaginary-eigenvalue reads on a + * per-decomposition isAllReal flag, and the device eigenvalue buffer is + * sized EIGEN_COMPLEX-only (matching CPU), not widened just because + * BEAGLE_FLAG_SPECTRAL_REPRESENTATION is in use. This test uses a + * *symmetric* circulant (r_fwd == r_bkd, so eigenvalues are exactly real, + * not just numerically close to zero) at a padded state count + * (17 -> kPaddedStateCount=32) and never requests + * BEAGLE_FLAG_EIGEN_COMPLEX, so the GPU device eigenvalue buffer is + * allocated at exactly kPaddedStateCount width -- no imaginary slack at + * all. If the kernel's imaginary-eigenvalue read were not properly gated, + * this configuration reads past the end of that buffer: the exact bug + * class PR #242 fixed and this follow-up refines. Run with + * BEAGLE_DEBUG_EIGEN=1 to see the device Eval buffer printed at its + * (narrow) width directly. + * ═══════════════════════════════════════════════════════════════════════*/ + +static int runTestRealOnly(int gpuDevice) { + const int S = 17, nPat = 4, nCats = 2; + const double tol = 1e-3; /* absolute tolerance for entries with |value|<=1 */ + const double relTol = 1e-2; /* relative tolerance for entries with |value|>1 (single vs double precision), same policy as runTest17 */ + + fprintf(stdout, "\n=== Real-only 17-state test (no BEAGLE_FLAG_EIGEN_COMPLEX) ===\n"); + fprintf(stdout, " r_fwd=r_bkd=1.0 -> symmetric/reversible circulant, all eigenvalues exactly real\n"); + fprintf(stdout, " Device eigenvalue buffer allocated at exactly kPaddedStateCount width (no imaginary slack)\n"); + + std::vector evec, ivec, eval; + buildCirculantN(S, evec, ivec, eval, /*r_fwd=*/1.0, /*r_bkd=*/1.0); + + fprintf(stdout, " Eigenvalue imag parts (must all be exactly 0): "); + bool anyImag = false; + for (int k = 0; k < S; k++) { + fprintf(stdout, "%.3f ", eval[S + k]); + if (eval[S + k] != 0.0) anyImag = true; + } + fprintf(stdout, "\n"); + if (anyImag) { + fprintf(stderr, " ERROR: expected all-zero imaginary eigenvalues for a symmetric circulant\n"); + return 1; + } + + std::vector freqs(S, 1.0 / S); + double rates_[2] = { 0.5, 1.5 }; + double catWts_[2] = { 0.5, 0.5 }; + double edgeLens_[4] = { 0.1, 0.1, 0.2, 0.5 }; + + int hSt[4] = { 3, 0, 5, 7 }; + int cSt[4] = { 3, 0, 5, 5 }; + int gSt[4] = { 0, 0, 0, 12 }; + + auto runOne = [&](bool useGpu, int dev, bool singlePrec) -> std::pair> { + int inst = createInstance(useGpu, dev, singlePrec, S, nPat, nCats, /*complexEigen=*/false); + if (inst < 0) return { 0.0, {} }; + + beagleSetTipStates(inst, 0, hSt); + beagleSetTipStates(inst, 1, cSt); + beagleSetTipStates(inst, 2, gSt); + + beagleSetCategoryRates(inst, rates_); + std::vector pw(nPat, 1.0); + beagleSetPatternWeights(inst, pw.data()); + beagleSetStateFrequencies(inst, 0, freqs.data()); + beagleSetCategoryWeights(inst, 0, catWts_); + beagleSetEigenDecomposition(inst, 0, evec.data(), ivec.data(), eval.data()); + + int nodeIdx[4] = { 0, 1, 2, 3 }; + beagleUpdateTransitionMatrices(inst, 0, nodeIdx, NULL, NULL, edgeLens_, 4); + + BeagleOperation postOps[2] = { + { 3, BEAGLE_OP_NONE, BEAGLE_OP_NONE, 0, 0, 1, 1 }, + { 4, BEAGLE_OP_NONE, BEAGLE_OP_NONE, 2, 2, 3, 3 } + }; + beagleUpdatePartials(inst, postOps, 2, BEAGLE_OP_NONE); + + double logL = 0.0; + int rootIdx = 4, cwIdx = 0, sfIdx = 0, csi = BEAGLE_OP_NONE; + beagleCalculateRootLogLikelihoods(inst, &rootIdx, &cwIdx, &sfIdx, &csi, 1, &logL); + + std::vector rootPre((size_t)nCats * nPat * S); + for (int c = 0; c < nCats; c++) + for (int p = 0; p < nPat; p++) + for (int s = 0; s < S; s++) + rootPre[(size_t)c * nPat * S + p * S + s] = freqs[s]; + beagleSetPartials(inst, 5, rootPre.data()); + + BeagleOperation preOps[4] = { + { 6, BEAGLE_OP_NONE, BEAGLE_OP_NONE, 5, BEAGLE_OP_NONE, 2, 2 }, + { 7, BEAGLE_OP_NONE, BEAGLE_OP_NONE, 5, BEAGLE_OP_NONE, 3, 3 }, + { 8, BEAGLE_OP_NONE, BEAGLE_OP_NONE, 6, BEAGLE_OP_NONE, 0, 0 }, + { 9, BEAGLE_OP_NONE, BEAGLE_OP_NONE, 6, BEAGLE_OP_NONE, 1, 1 } + }; + beagleUpdatePrePartials_v5(inst, preOps, 4, BEAGLE_OP_NONE, BEAGLE_PARTIALS_TOP); + + BeagleBranchOperation branchOps[4] = { + { 1, 8, 1, 0 }, + { 0, 9, 0, 0 }, + { 2, 7, 2, 0 }, + { 3, 6, 3, 0 } + }; + std::vector grad(S * S, 0.0); + beagleCalculateAdjointDerivative(inst, branchOps, 0, 0, 4, 0, 4, grad.data(), NULL); + + beagleFinalizeInstance(inst); + return { logL, grad }; + }; + + fprintf(stdout, "\nCreating CPU (double) instance:\n"); + auto [cpuLogL, cpuGrad] = runOne(false, -1, false); + if (cpuGrad.empty()) return 1; + fprintf(stdout, "CPU log-likelihood: %.6f\n", cpuLogL); + + if (gpuDevice < 0) return 0; + + fprintf(stdout, "\nCreating GPU (single) instance:\n"); + auto [gpuLogL, gpuGrad] = runOne(true, gpuDevice, true); + if (gpuGrad.empty()) return 1; + fprintf(stdout, "GPU log-likelihood: %.6f\n", gpuLogL); + + /* Gradient magnitudes here range up to ~190 (same character as runTest17's + * complex-eigenvalue case), so single-precision GPU vs double-precision + * CPU error scales with the *largest* magnitude in the computation -- + * threshold against the matrix-wide max magnitude instead of a flat + * absolute tolerance (same policy as runTest17). */ + double maxMagnitude = 0.0; + for (int k = 0; k < S*S; k++) { + if (!std::isnan(cpuGrad[k])) maxMagnitude = std::max(maxMagnitude, fabs(cpuGrad[k])); + if (!std::isnan(gpuGrad[k])) maxMagnitude = std::max(maxMagnitude, fabs(gpuGrad[k])); + } + const double threshold = std::max(tol, relTol * maxMagnitude); + + bool gradOk = true; + double maxDiff = 0.0; + for (int k = 0; k < S*S; k++) { + double d = fabs(cpuGrad[k] - gpuGrad[k]); + if (d > maxDiff || std::isnan(d)) maxDiff = d; + if (d > threshold || std::isnan(d)) gradOk = false; + } + fprintf(stdout, "Real-only 17-state gradient max|diff|=%.3e threshold=%.3e (rel to max|value|=%.3e) %s\n", + maxDiff, threshold, maxMagnitude, gradOk ? "PASS" : "FAIL"); + return gradOk ? 0 : 1; +} + +/* ═══════════════════════════════════════════════════════════════════════ + * Dynamic real <-> complex eigendecomposition transition test (16-state) + * + * A single BEAGLE_FLAG_EIGEN_COMPLEX-flagged instance (device eigenvalue + * buffer wide enough for either case) has beagleSetEigenDecomposition + * called three times in sequence for the same eigenIndex, alternating + * between a genuinely real-only decomposition and a genuinely complex one + * -- mirroring how one MCMC run can produce either kind of decomposition + * from one iteration to the next for a general asymmetric-rate model (see + * hEigenDecompIsAllReal's recomputation on every setEigenDecomposition + * call). Confirms the isAllReal gating doesn't leave stale sDs/sCs kernel + * state, or stale device Eval data, bleeding from one call into the next. + * ═══════════════════════════════════════════════════════════════════════*/ + +static int runTestDynamicTransition(int gpuDevice) { + const int S = 16, nPat = 4, nCats = 2; + const double tol = 1e-3; + + fprintf(stdout, "\n=== Dynamic real<->complex eigendecomposition transition test (16-state) ===\n"); + + std::vector evecReal, ivecReal, evalReal; + buildCirculantN(S, evecReal, ivecReal, evalReal, /*r_fwd=*/1.0, /*r_bkd=*/1.0); + std::vector evecCplx, ivecCplx, evalCplx; + buildCirculantN(S, evecCplx, ivecCplx, evalCplx, /*r_fwd=*/1.0, /*r_bkd=*/0.5); + + std::vector freqs(S, 1.0 / S); + double rates_[2] = { 0.5, 1.5 }; + double catWts_[2] = { 0.5, 0.5 }; + double edgeLens_[4] = { 0.1, 0.1, 0.2, 0.5 }; + int hSt[4] = { 3, 0, 5, 7 }; + int cSt[4] = { 3, 0, 5, 5 }; + int gSt[4] = { 0, 0, 0, 7 }; + + auto setupAndRun = [&](int inst, const std::vector& evec, + const std::vector& ivec, const std::vector& eval) + -> std::pair> { + beagleSetEigenDecomposition(inst, 0, evec.data(), ivec.data(), eval.data()); + int nodeIdx[4] = { 0, 1, 2, 3 }; + beagleUpdateTransitionMatrices(inst, 0, nodeIdx, NULL, NULL, edgeLens_, 4); + BeagleOperation postOps[2] = { + { 3, BEAGLE_OP_NONE, BEAGLE_OP_NONE, 0, 0, 1, 1 }, + { 4, BEAGLE_OP_NONE, BEAGLE_OP_NONE, 2, 2, 3, 3 } + }; + beagleUpdatePartials(inst, postOps, 2, BEAGLE_OP_NONE); + double logL = 0.0; + int rootIdx = 4, cwIdx = 0, sfIdx = 0, csi = BEAGLE_OP_NONE; + beagleCalculateRootLogLikelihoods(inst, &rootIdx, &cwIdx, &sfIdx, &csi, 1, &logL); + std::vector rootPre((size_t)nCats * nPat * S); + for (int c = 0; c < nCats; c++) + for (int p = 0; p < nPat; p++) + for (int s = 0; s < S; s++) + rootPre[(size_t)c * nPat * S + p * S + s] = freqs[s]; + beagleSetPartials(inst, 5, rootPre.data()); + BeagleOperation preOps[4] = { + { 6, BEAGLE_OP_NONE, BEAGLE_OP_NONE, 5, BEAGLE_OP_NONE, 2, 2 }, + { 7, BEAGLE_OP_NONE, BEAGLE_OP_NONE, 5, BEAGLE_OP_NONE, 3, 3 }, + { 8, BEAGLE_OP_NONE, BEAGLE_OP_NONE, 6, BEAGLE_OP_NONE, 0, 0 }, + { 9, BEAGLE_OP_NONE, BEAGLE_OP_NONE, 6, BEAGLE_OP_NONE, 1, 1 } + }; + beagleUpdatePrePartials_v5(inst, preOps, 4, BEAGLE_OP_NONE, BEAGLE_PARTIALS_TOP); + BeagleBranchOperation branchOps[4] = { + { 1, 8, 1, 0 }, + { 0, 9, 0, 0 }, + { 2, 7, 2, 0 }, + { 3, 6, 3, 0 } + }; + std::vector grad(S * S, 0.0); + beagleCalculateAdjointDerivative(inst, branchOps, 0, 0, 4, 0, 4, grad.data(), NULL); + return { logL, grad }; + }; + + auto makeInst = [&](bool useGpu, int dev, bool singlePrec) -> int { + int inst = createInstance(useGpu, dev, singlePrec, S, nPat, nCats, /*complexEigen=*/true); + if (inst < 0) return inst; + beagleSetTipStates(inst, 0, hSt); + beagleSetTipStates(inst, 1, cSt); + beagleSetTipStates(inst, 2, gSt); + beagleSetCategoryRates(inst, rates_); + std::vector pw(nPat, 1.0); + beagleSetPatternWeights(inst, pw.data()); + beagleSetStateFrequencies(inst, 0, freqs.data()); + beagleSetCategoryWeights(inst, 0, catWts_); + return inst; + }; + + fprintf(stdout, "\nCreating CPU (double) instance (BEAGLE_FLAG_EIGEN_COMPLEX requested):\n"); + int cpuInst = makeInst(false, -1, false); + if (cpuInst < 0) { fprintf(stderr, "Failed CPU instance\n"); return 1; } + + const char* stepNames[3] = { "real", "complex", "real (again)" }; + const std::vector* stepEvec[3] = { &evecReal, &evecCplx, &evecReal }; + const std::vector* stepIvec[3] = { &ivecReal, &ivecCplx, &ivecReal }; + const std::vector* stepEval[3] = { &evalReal, &evalCplx, &evalReal }; + + double cpuLogL[3]; + std::vector cpuGrad[3]; + for (int i = 0; i < 3; i++) { + auto r = setupAndRun(cpuInst, *stepEvec[i], *stepIvec[i], *stepEval[i]); + cpuLogL[i] = r.first; + cpuGrad[i] = r.second; + fprintf(stdout, " [CPU step %d: %-13s] logL=%.6f\n", i, stepNames[i], cpuLogL[i]); + } + beagleFinalizeInstance(cpuInst); + + if (gpuDevice < 0) return 0; + + fprintf(stdout, "\nCreating GPU (single) instance (BEAGLE_FLAG_EIGEN_COMPLEX requested):\n"); + int gpuInst = makeInst(true, gpuDevice, true); + if (gpuInst < 0) { fprintf(stderr, "Failed GPU instance\n"); return 1; } + + /* The "complex" step reuses runTest16's exact asymmetric circulant + * (r_fwd=1.0, r_bkd=0.5), which has a documented pre-existing marginal + * single-vs-double-precision diff of ~1.02e-03 (CLUSTER_AGENT_FINDINGS2.md + * section 8.2) unrelated to eigen-real/complex gating. Give that one step + * a matching looser tolerance so this test isn't tripped by that same + * known, already-accepted noise floor; the "real" steps (the actual + * subject of this test) are held to the tighter tol. */ + const double tolComplexStep = 1.1e-3; + std::vector gpuGradSteps[3]; + bool allOk = true; + for (int i = 0; i < 3; i++) { + auto r = setupAndRun(gpuInst, *stepEvec[i], *stepIvec[i], *stepEval[i]); + double gpuLogL = r.first; + gpuGradSteps[i] = r.second; + fprintf(stdout, " [GPU step %d: %-13s] logL=%.6f\n", i, stepNames[i], gpuLogL); + + double maxDiff = 0.0; + for (int k = 0; k < S*S; k++) { + double d = fabs(cpuGrad[i][k] - gpuGradSteps[i][k]); + if (d > maxDiff || std::isnan(d)) maxDiff = d; + } + double stepTol = (i == 1) ? tolComplexStep : tol; + bool stepOk = (maxDiff <= stepTol) && !std::isnan(maxDiff); + fprintf(stdout, " step %d (%s) gradient max|diff|=%.3e (tol=%.1e) %s\n", + i, stepNames[i], maxDiff, stepTol, stepOk ? "PASS" : "FAIL"); + allOk &= stepOk; + } + + /* Consistency check: step 0 and step 2 use the IDENTICAL real-only + * decomposition, with a genuinely complex decomposition (step 1) + * processed on the same instance/eigenIndex in between. If isAllReal + * gating left any stale sDs/sCs kernel state or stale device Eval data + * from the complex step bleeding into the next real step, step 2 would + * diverge from step 0 -- this is the most direct test for that failure + * mode, independent of any CPU-vs-GPU precision noise. + * Tolerance is not zero: verified empirically (by rerunning with all + * three steps set to the same real decomposition, i.e. no complex step + * at all) that back-to-back identical GPU launches already differ by + * this same ~1e-5 magnitude on this hardware -- ADJOINT_ATOMIC_ADD_GPU's + * float addition is not associative, so summation order across + * concurrent atomics is not guaranteed bit-identical between launches. + * That noise floor is unrelated to eigen real/complex gating, so the + * threshold here is set an order of magnitude above it. */ + double stepConsistencyDiff = 0.0; + for (int k = 0; k < S*S; k++) + stepConsistencyDiff = std::max(stepConsistencyDiff, fabs(gpuGradSteps[0][k] - gpuGradSteps[2][k])); + bool consistencyOk = stepConsistencyDiff < 1e-4; + fprintf(stdout, " GPU step0 vs step2 (both real, complex processed in between) max|diff|=%.3e %s\n", + stepConsistencyDiff, consistencyOk ? "PASS (no stale state)" : "FAIL (possible stale state)"); + allOk &= consistencyOk; + beagleFinalizeInstance(gpuInst); + + fprintf(stdout, "Dynamic real<->complex transition test: %s\n", allOk ? "PASS" : "FAIL"); + return allOk ? 0 : 1; +} + /* ── main ────────────────────────────────────────────────────────────── */ int main(int argc, const char *argv[]) { @@ -941,9 +1258,11 @@ int main(int argc, const char *argv[]) { bool useGpu = false; int whichDevice = -1; - bool run4 = true; - bool run16 = false; - bool run17 = false; + bool run4 = true; + bool run16 = false; + bool run17 = false; + bool runRealOnly = false; + bool runDynamic = false; for (int i = 1; i < argc; i++) { if (strcmp(argv[i], "--gpu") == 0 && i+1 < argc) { @@ -954,17 +1273,25 @@ int main(int argc, const char *argv[]) { run4 = (ns == 4); run16 = (ns == 16); run17 = (ns == 17); + } else if (strcmp(argv[i], "--realonly") == 0) { + run4 = false; + runRealOnly = true; + } else if (strcmp(argv[i], "--dynamic") == 0) { + run4 = false; + runDynamic = true; } else if (strcmp(argv[i], "--all") == 0) { - run4 = run16 = run17 = true; + run4 = run16 = run17 = runRealOnly = runDynamic = true; } } int gpuDev = useGpu ? whichDevice : -1; int result = 0; - if (run4) result |= runTest4(gpuDev); - if (run16) result |= runTest16(gpuDev); - if (run17) result |= runTest17(gpuDev); + if (run4) result |= runTest4(gpuDev); + if (run16) result |= runTest16(gpuDev); + if (run17) result |= runTest17(gpuDev); + if (runRealOnly) result |= runTestRealOnly(gpuDev); + if (runDynamic) result |= runTestDynamicTransition(gpuDev); return result; } From cc73bdf25b30a89a07716b10ede76a05592ee986 Mon Sep 17 00:00:00 2001 From: fil-monti Date: Fri, 10 Jul 2026 22:58:09 -0700 Subject: [PATCH 06/13] Add BEAGLE_PROFILE_SUMMARY diagnostic instrumentation, hoist redundant getenv() Investigating why GPU BEAGLE was ~16-17x slower than CPU on a real analysis (see beagle-bugs/gpu-spectral-per-branch-distance-batching/PLAN.md), no ROCm profiler was available on the target machine, so this adds targeted accumulating-counter instrumentation instead: call counts and elapsed host-side wall time for LaunchKernel, LaunchKernelConcurrent, SynchronizeHost (clFinish), and MemcpyHostToDevice, gated behind a new BEAGLE_PROFILE_SUMMARY env var and printed once via atexit() (not per-call, to avoid repeating an earlier session's mistake of an unconditional per-call print flooding a benchmark log with millions of lines). A warmup window discards the first 300 SynchronizeHost calls before any counter starts accumulating, so the summary reflects steady-state behavior rather than one-time startup costs. Also fixes a real, if minor, inefficiency found while reading this code: LaunchKernel/LaunchKernelConcurrent were calling getenv("BEAGLE_DEBUG_KERNEL_ARGS") repeatedly (11-15 times per single kernel launch) instead of once; hoisted to a single check per call. This instrumentation is what identified the root cause fixed in the next commit (6.5 million individually-blocking MemcpyHostToDevice calls dominating GPU wall-clock time) and remains available as a standing, default-off profiling knob. --- libhmsbeagle/GPU/GPUInterfaceOpenCL.cpp | 157 ++++++++++++++++++++++-- 1 file changed, 144 insertions(+), 13 deletions(-) diff --git a/libhmsbeagle/GPU/GPUInterfaceOpenCL.cpp b/libhmsbeagle/GPU/GPUInterfaceOpenCL.cpp index d2a4aada..0543a1b4 100644 --- a/libhmsbeagle/GPU/GPUInterfaceOpenCL.cpp +++ b/libhmsbeagle/GPU/GPUInterfaceOpenCL.cpp @@ -22,6 +22,7 @@ #include #include #include +#include #include "libhmsbeagle/beagle.h" #include "libhmsbeagle/GPU/GPUImplDefs.h" @@ -29,6 +30,78 @@ #include "libhmsbeagle/GPU/GPUInterface.h" #include "libhmsbeagle/GPU/KernelResource.h" +// --------------------------------------------------------------------------- +// Diagnostic-only profiling instrumentation, gated behind BEAGLE_PROFILE_SUMMARY. +// Accumulates call counts + elapsed host-side time into static counters (no +// per-call printing) and dumps one summary at process exit. A warmup window +// (first BEAGLE_PROFILE_WARMUP_SYNC SynchronizeHost calls, i.e. host<->device +// round trips) is discarded before any counter starts accumulating, so the +// summary reflects steady-state behavior rather than one-time clBuildProgram/ +// initial-allocation/GPU-clock-ramp costs. Not intended to be a permanent +// feature; safe to leave disabled (default off, near-zero overhead: one +// getenv + one bool check per call) or revert entirely. +// --------------------------------------------------------------------------- +namespace { + +const long BEAGLE_PROFILE_WARMUP_SYNC = 300; + +struct BeagleProfileStats { + bool registered = false; + bool enabled = false; + bool warm = false; + long warmupRemaining = BEAGLE_PROFILE_WARMUP_SYNC; + + long launchKernelCalls = 0; + double launchKernelSeconds = 0.0; + long launchKernelConcurrentCalls = 0; + double launchKernelConcurrentSeconds = 0.0; + long syncHostCalls = 0; + double syncHostSeconds = 0.0; + long memcpyH2DCalls = 0; + double memcpyH2DSeconds = 0.0; + size_t memcpyH2DBytes = 0; +}; + +BeagleProfileStats gBeagleProfile; + +void beagleProfilePrintSummary() { + if (!gBeagleProfile.enabled) return; + BeagleProfileStats& p = gBeagleProfile; + double totalAccounted = p.launchKernelSeconds + p.launchKernelConcurrentSeconds + + p.syncHostSeconds + p.memcpyH2DSeconds; + fprintf(stderr, "\n===== BEAGLE_PROFILE_SUMMARY (GPUInterfaceOpenCL) =====\n"); + fprintf(stderr, "warmup discarded: first %ld SynchronizeHost calls\n", BEAGLE_PROFILE_WARMUP_SYNC); + fprintf(stderr, "LaunchKernel : calls=%-8ld total=%.6fs avg=%.3fus\n", + p.launchKernelCalls, p.launchKernelSeconds, + p.launchKernelCalls ? p.launchKernelSeconds * 1e6 / p.launchKernelCalls : 0.0); + fprintf(stderr, "LaunchKernelConcurrent : calls=%-8ld total=%.6fs avg=%.3fus\n", + p.launchKernelConcurrentCalls, p.launchKernelConcurrentSeconds, + p.launchKernelConcurrentCalls ? p.launchKernelConcurrentSeconds * 1e6 / p.launchKernelConcurrentCalls : 0.0); + fprintf(stderr, "SynchronizeHost(clFinish): calls=%-8ld total=%.6fs avg=%.3fus\n", + p.syncHostCalls, p.syncHostSeconds, + p.syncHostCalls ? p.syncHostSeconds * 1e6 / p.syncHostCalls : 0.0); + fprintf(stderr, "MemcpyHostToDevice : calls=%-8ld total=%.6fs avg=%.3fus bytes=%zu\n", + p.memcpyH2DCalls, p.memcpyH2DSeconds, + p.memcpyH2DCalls ? p.memcpyH2DSeconds * 1e6 / p.memcpyH2DCalls : 0.0, + p.memcpyH2DBytes); + fprintf(stderr, "Total host-side time accounted above (post-warmup): %.6fs\n", totalAccounted); + fprintf(stderr, "===== END BEAGLE_PROFILE_SUMMARY =====\n\n"); + fflush(stderr); +} + +inline bool beagleProfilingEnabled() { + if (!gBeagleProfile.registered) { + gBeagleProfile.enabled = (getenv("BEAGLE_PROFILE_SUMMARY") != NULL); + gBeagleProfile.registered = true; + if (gBeagleProfile.enabled) { + atexit(beagleProfilePrintSummary); + } + } + return gBeagleProfile.enabled; +} + +} // anonymous namespace + #define SAFE_CL(call) { \ int error = call; \ if(error != CL_SUCCESS) { \ @@ -477,8 +550,24 @@ void GPUInterface::SynchronizeHost() { // SAFE_CL(clFinish(openClCommandQueues[i])); // } + bool prof = beagleProfilingEnabled(); + std::chrono::high_resolution_clock::time_point t0; + if (prof) t0 = std::chrono::high_resolution_clock::now(); + SAFE_CL(clFinish(openClCommandQueues[0])); + if (prof) { + double dt = std::chrono::duration( + std::chrono::high_resolution_clock::now() - t0).count(); + if (gBeagleProfile.warmupRemaining > 0) { + gBeagleProfile.warmupRemaining--; + if (gBeagleProfile.warmupRemaining == 0) gBeagleProfile.warm = true; + } else { + gBeagleProfile.syncHostCalls++; + gBeagleProfile.syncHostSeconds += dt; + } + } + #ifdef BEAGLE_DEBUG_FLOW fprintf(stderr,"\t\t\tLeaving GPUInterface::SynchronizeHost\n"); #endif @@ -545,17 +634,24 @@ void GPUInterface::LaunchKernel(GPUFunction deviceFunction, #endif char kernelNameBuf[256] = {0}; - if (getenv("BEAGLE_DEBUG_KERNEL_ARGS")) { + // Hoisted once instead of calling getenv() again on every parameter in the loops + // below (previously re-checked per-arg, tens of getenv() calls per kernel launch). + const bool debugArgs = (getenv("BEAGLE_DEBUG_KERNEL_ARGS") != NULL); + if (debugArgs) { clGetKernelInfo(deviceFunction, CL_KERNEL_FUNCTION_NAME, sizeof(kernelNameBuf), kernelNameBuf, NULL); fprintf(stderr, "[LaunchKernel] %s: parameterCountV=%d totalParameterCount=%d\n", kernelNameBuf, parameterCountV, totalParameterCount); } + bool prof = beagleProfilingEnabled(); + std::chrono::high_resolution_clock::time_point t0; + if (prof) t0 = std::chrono::high_resolution_clock::now(); + 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")) { + if (debugArgs) { fprintf(stderr, "[LaunchKernel] %s: ptr arg[%d] = %p\n", kernelNameBuf, i, param); fflush(stderr); } @@ -564,7 +660,7 @@ void GPUInterface::LaunchKernel(GPUFunction deviceFunction, } for(int i = parameterCountV; i < totalParameterCount; i++) { unsigned int param = va_arg(parameters, unsigned int); - if (getenv("BEAGLE_DEBUG_KERNEL_ARGS")) { + if (debugArgs) { fprintf(stderr, "[LaunchKernel] %s: int arg[%d] = %u (signed: %d)\n", kernelNameBuf, i, param, (int)param); fflush(stderr); } @@ -584,7 +680,7 @@ void GPUInterface::LaunchKernel(GPUFunction deviceFunction, globalWorkSize[1] = block.y * grid.y; globalWorkSize[2] = block.z * grid.z; - if (getenv("BEAGLE_DEBUG_KERNEL_ARGS")) { + if (debugArgs) { 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], @@ -602,7 +698,7 @@ void GPUInterface::LaunchKernel(GPUFunction deviceFunction, printf("local = %lu\n\n", local); #endif - if (getenv("BEAGLE_DEBUG_KERNEL_ARGS")) { + if (debugArgs) { size_t maxLocal = 0; clGetKernelWorkGroupInfo(deviceFunction, openClDeviceId, CL_KERNEL_WORK_GROUP_SIZE, sizeof(maxLocal), &maxLocal, NULL); @@ -614,7 +710,7 @@ void GPUInterface::LaunchKernel(GPUFunction deviceFunction, fflush(stderr); } - if (getenv("BEAGLE_DEBUG_KERNEL_ARGS")) { + if (debugArgs) { fprintf(stderr, "[LaunchKernel] %s: about to enqueue\n", kernelNameBuf); fflush(stderr); } @@ -630,7 +726,14 @@ void GPUInterface::LaunchKernel(GPUFunction deviceFunction, globalWorkSize, localWorkSize, 0, NULL, NULL)); } - if (getenv("BEAGLE_DEBUG_KERNEL_ARGS")) { + if (prof && gBeagleProfile.warm) { + double dt = std::chrono::duration( + std::chrono::high_resolution_clock::now() - t0).count(); + gBeagleProfile.launchKernelCalls++; + gBeagleProfile.launchKernelSeconds += dt; + } + + if (debugArgs) { fprintf(stderr, "[LaunchKernel] %s: enqueued OK, calling clFinish\n", kernelNameBuf); fflush(stderr); clFinish(openClCommandQueues[0]); @@ -656,17 +759,24 @@ void GPUInterface::LaunchKernelConcurrent(GPUFunction deviceFunction, #endif char kernelNameBuf[256] = {0}; - if (getenv("BEAGLE_DEBUG_KERNEL_ARGS")) { + // Hoisted once instead of calling getenv() again on every parameter in the loops + // below (previously re-checked per-arg, tens of getenv() calls per kernel launch). + const bool debugArgs = (getenv("BEAGLE_DEBUG_KERNEL_ARGS") != NULL); + if (debugArgs) { clGetKernelInfo(deviceFunction, CL_KERNEL_FUNCTION_NAME, sizeof(kernelNameBuf), kernelNameBuf, NULL); fprintf(stderr, "[LaunchKernelConcurrent] %s: parameterCountV=%d totalParameterCount=%d streamIndex=%d waitIndex=%d\n", kernelNameBuf, parameterCountV, totalParameterCount, streamIndex, waitIndex); } + bool prof = beagleProfilingEnabled(); + std::chrono::high_resolution_clock::time_point t0; + if (prof) t0 = std::chrono::high_resolution_clock::now(); + 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")) { + if (debugArgs) { fprintf(stderr, "[LaunchKernelConcurrent] %s: ptr arg[%d] = %p\n", kernelNameBuf, i, param); fflush(stderr); } @@ -675,7 +785,7 @@ void GPUInterface::LaunchKernelConcurrent(GPUFunction deviceFunction, } for(int i = parameterCountV; i < totalParameterCount; i++) { unsigned int param = va_arg(parameters, unsigned int); - if (getenv("BEAGLE_DEBUG_KERNEL_ARGS")) { + if (debugArgs) { fprintf(stderr, "[LaunchKernelConcurrent] %s: int arg[%d] = %u (signed: %d)\n", kernelNameBuf, i, param, (int)param); fflush(stderr); } @@ -695,7 +805,7 @@ void GPUInterface::LaunchKernelConcurrent(GPUFunction deviceFunction, globalWorkSize[1] = block.y * grid.y; globalWorkSize[2] = block.z * grid.z; - if (getenv("BEAGLE_DEBUG_KERNEL_ARGS")) { + if (debugArgs) { 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], @@ -720,7 +830,7 @@ void GPUInterface::LaunchKernelConcurrent(GPUFunction deviceFunction, dims = 2; } - if (getenv("BEAGLE_DEBUG_KERNEL_ARGS")) { + if (debugArgs) { fprintf(stderr, "[LaunchKernelConcurrent] %s: about to enqueue (dims=%d)\n", kernelNameBuf, dims); fflush(stderr); } @@ -764,7 +874,14 @@ void GPUInterface::LaunchKernelConcurrent(GPUFunction deviceFunction, globalWorkSize, localWorkSize, 0, NULL, NULL)); // } - if (getenv("BEAGLE_DEBUG_KERNEL_ARGS")) { + if (prof && gBeagleProfile.warm) { + double dt = std::chrono::duration( + std::chrono::high_resolution_clock::now() - t0).count(); + gBeagleProfile.launchKernelConcurrentCalls++; + gBeagleProfile.launchKernelConcurrentSeconds += dt; + } + + if (debugArgs) { fprintf(stderr, "[LaunchKernelConcurrent] %s: enqueued OK, calling clFinish\n", kernelNameBuf); fflush(stderr); clFinish(openClCommandQueues[0]); @@ -1034,9 +1151,23 @@ void GPUInterface::MemcpyHostToDevice(GPUPtr dest, fprintf(stderr, "\t\t\tEntering GPUInterface::MemcpyHostToDevice\n"); #endif + bool prof = beagleProfilingEnabled(); + std::chrono::high_resolution_clock::time_point t0; + if (prof) t0 = std::chrono::high_resolution_clock::now(); + + // CL_TRUE => blocking write: this call does not return until the transfer + // completes, regardless of payload size. SAFE_CL(clEnqueueWriteBuffer(openClCommandQueues[0], dest, CL_TRUE, 0, memSize, src, 0, NULL, NULL)); + if (prof && gBeagleProfile.warm) { + double dt = std::chrono::duration( + std::chrono::high_resolution_clock::now() - t0).count(); + gBeagleProfile.memcpyH2DCalls++; + gBeagleProfile.memcpyH2DSeconds += dt; + gBeagleProfile.memcpyH2DBytes += memSize; + } + #ifdef BEAGLE_DEBUG_FLOW fprintf(stderr, "\t\t\tLeaving GPUInterface::MemcpyHostToDevice\n"); #endif From 9f6c2bcddb827169e2e42e0aa4a6417bc33a46ae Mon Sep 17 00:00:00 2001 From: fil-monti Date: Fri, 10 Jul 2026 22:59:43 -0700 Subject: [PATCH 07/13] Remove unconditional dispatchPruneSS trace prints These were added as temporary diagnostics while tracking down a cross-plugin symbol-interposition bug (see the "Add unconditional dispatchPruneSS trace prints (temporary, not env-gated)" commit) and flagged there as needing cleanup before this branch was done -- they print on every single postorder pruning dispatch, unconditionally (no env-var gate, unlike every other debug aid in this codebase), and were confirmed to be actively drowning out real signal in a benchmark log during later performance work. No longer needed; whatever they were diagnosing has long since been fixed and verified. --- libhmsbeagle/GPU/BeagleGPUImpl.hpp | 1 - libhmsbeagle/GPU/BeagleGPUSpectralImpl.hpp | 1 - 2 files changed, 2 deletions(-) diff --git a/libhmsbeagle/GPU/BeagleGPUImpl.hpp b/libhmsbeagle/GPU/BeagleGPUImpl.hpp index ab16b838..f8eac334 100644 --- a/libhmsbeagle/GPU/BeagleGPUImpl.hpp +++ b/libhmsbeagle/GPU/BeagleGPUImpl.hpp @@ -4759,7 +4759,6 @@ 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 cddc45db..ed154972 100644 --- a/libhmsbeagle/GPU/BeagleGPUSpectralImpl.hpp +++ b/libhmsbeagle/GPU/BeagleGPUSpectralImpl.hpp @@ -194,7 +194,6 @@ 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 224bd6cd00c488b42e3b93d30922794c7ea92611 Mon Sep 17 00:00:00 2001 From: fil-monti Date: Fri, 10 Jul 2026 23:01:08 -0700 Subject: [PATCH 08/13] GPU spectral: batch per-branch distance uploads instead of one memcpy per branch Profiling found GPU BEAGLE spending 82.6% of measured GPU-interface time (100.4s of 121.4s accounted, on a 3000-iteration real analysis) in 6.5 million individually-blocking ~100-byte MemcpyHostToDevice calls -- one per tree branch per HMC leapfrog step, updating that branch's distance (effective branch length x category rate) in BeagleGPUSpectralImpl:: updateTransitionMatrices. Actual GPU compute (clFinish wait) was only 3.5% of that time; the bottleneck was round-trip latency paid millions of times over for tiny transfers, not bandwidth. Fixes it by mirroring a pattern already proven in this same codebase: BeagleGPUImpl::updateTransitionMatrices (the base class, called directly above the code this changes) already solves the identical class of problem for its dense-matrix update via gather-batch-scatter -- build one flat host-side queue of (destination-offset, value) pairs, upload it with exactly two MemcpyHostToDevice calls regardless of branch count, then dispatch one GPU kernel that scatters each value into its device-memory destination. Applies that same pattern here: - New dedicated queue buffers (dSpectralPtrQueue/dSpectralDistanceQueue + host mirrors), sized kMatrixCount * kCategoryCount. Not reusing the base class's own queue buffers: they're declared `private:` in BeagleGPUImpl, unreachable from this subclass, so dedicated buffers were required, not just safer. - New trivial kernelScatterSpectralDistances kernel (one thread per queue entry, no PADDED_STATE_COUNT/STATE_COUNT dependence), added identically to both kernelsSpectralIfDef.cu and kernelsSpectralIfDef4.cu (matching this codebase's existing precedent of duplicating spectral-only kernel logic into both files rather than the shared dense-kernel file). - Destination offsets computed via the existing kSpectralDistanceStrideElements member (probabilityIndices[i] * kSpectralDistanceStrideElements + j) -- no new addressing arithmetic needed. - New KernelLauncher::ScatterSpectralDistances wrapper, deliberately without a synchronizing call after the launch (matching GetTransitionProbabilitiesSquare's own precedent) -- adding one would reintroduce exactly the kind of round-trip stall this fix removes. Verified correctness (adjointtest4 --gpu 1 --all unchanged from baseline, including a controlled A/B confirming a pre-existing marginal 16-state precision-tolerance FAIL is unrelated to this change) and performance (same 3000-iteration benchmark, BEAGLE_PROFILE_SUMMARY-instrumented): MemcpyHostToDevice calls dropped 6,546,621 -> 105,772 (61.9x fewer), MemcpyHostToDevice time 100.4s -> 2.1s (47.7x less), total BEAST wall time 267.9s -> 173.3s (1.55x faster). CPU-vs-GPU ratio on this analysis improved from ~15.9x GPU-slower to ~10.3x GPU-slower -- GPU still slower overall for a problem this small, but the specific round-trip-latency bottleneck this session targeted is resolved; most of the remaining gap is JVM/Java-side work outside BEAGLE's GPU interface entirely (only ~15% of BEAST's wall time is now spent in the instrumented GPU-interface layer at all). See beagle-bugs/gpu-spectral-per-branch-distance-batching/PLAN.md for the full design rationale. --- libhmsbeagle/GPU/BeagleGPUSpectralImpl.h | 16 ++++++++++ libhmsbeagle/GPU/BeagleGPUSpectralImpl.hpp | 29 ++++++++++++++++--- libhmsbeagle/GPU/KernelLauncher.cpp | 23 +++++++++++++++ libhmsbeagle/GPU/KernelLauncher.h | 15 ++++++++++ .../GPU/kernels/kernelsSpectralIfDef.cu | 18 ++++++++++++ .../GPU/kernels/kernelsSpectralIfDef4.cu | 19 ++++++++++++ 6 files changed, 116 insertions(+), 4 deletions(-) diff --git a/libhmsbeagle/GPU/BeagleGPUSpectralImpl.h b/libhmsbeagle/GPU/BeagleGPUSpectralImpl.h index 56c8baa9..607cd72b 100644 --- a/libhmsbeagle/GPU/BeagleGPUSpectralImpl.h +++ b/libhmsbeagle/GPU/BeagleGPUSpectralImpl.h @@ -45,6 +45,22 @@ class BeagleGPUSpectralImpl : public BeagleGPUImpl { * adjoint offset-queue) must multiply by this, not by kCategoryCount. */ unsigned int kSpectralDistanceStrideElements; + /* Gather-batch-scatter queue for updateTransitionMatrices' per-branch + * distance uploads: one flat host-side (destination-offset, value) pair + * per (branch, category), uploaded with two MemcpyHostToDevice calls + * regardless of branch count, then scattered into dSpectralDistancesOrigin + * by one GPU kernel -- mirrors BeagleGPUImpl's own dPtrQueue/dDistanceQueue + * pattern for its dense-matrix update. Dedicated buffers rather than + * reusing that base-class pair: those are declared in BeagleGPUImpl's + * `private:` section, not `protected:`, so this derived class cannot + * reach them without changing the base class's access specifiers. + * Sized kMatrixCount * kCategoryCount, the largest a single + * updateTransitionMatrices call can need. */ + GPUPtr dSpectralPtrQueue; + GPUPtr dSpectralDistanceQueue; + unsigned int* hSpectralPtrQueue; + Real* hSpectralDistanceQueue; + /* Backward eigenvector buffers for the parent branch in pre-order. * dEvecT[ei] = U stored row-major (dEvecT[j*S+k] = U[j,k]). * Used as ievc1 in Growing kernels (backward Phase 1: U^T·p). diff --git a/libhmsbeagle/GPU/BeagleGPUSpectralImpl.hpp b/libhmsbeagle/GPU/BeagleGPUSpectralImpl.hpp index ed154972..34bf34d0 100644 --- a/libhmsbeagle/GPU/BeagleGPUSpectralImpl.hpp +++ b/libhmsbeagle/GPU/BeagleGPUSpectralImpl.hpp @@ -34,6 +34,8 @@ namespace gpu { BEAGLE_GPU_TEMPLATE BeagleGPUSpectralImpl::BeagleGPUSpectralImpl() : dSpectralDistancesOrigin(0), dSpectralDistances(NULL), hEigenIndexForMatrix(NULL), + dSpectralPtrQueue(0), dSpectralDistanceQueue(0), + hSpectralPtrQueue(NULL), hSpectralDistanceQueue(NULL), dEvecTOrigin(0), dEvecT(NULL), dIevcTOrigin(0), dIevcT(NULL), dGradientOrigin(0), dGradient(NULL), hEigenDecompIsAllReal(NULL), dAdjointQueue(0), hAdjointQueue(NULL), kAdjointQueueCapacity(0), @@ -45,6 +47,10 @@ BeagleGPUSpectralImpl::~BeagleGPUSpectralImpl() { GPUInterface* gpuIf = this->gpu; if (gpuIf) { if (dSpectralDistancesOrigin) gpuIf->FreeMemory(dSpectralDistancesOrigin); + if (dSpectralPtrQueue) gpuIf->FreeMemory(dSpectralPtrQueue); + if (dSpectralDistanceQueue) gpuIf->FreeMemory(dSpectralDistanceQueue); + if (hSpectralPtrQueue) gpuIf->FreeHostMemory(hSpectralPtrQueue); + if (hSpectralDistanceQueue) gpuIf->FreeHostMemory(hSpectralDistanceQueue); if (dEvecTOrigin) gpuIf->FreeMemory(dEvecTOrigin); if (dIevcTOrigin) gpuIf->FreeMemory(dIevcTOrigin); if (dGradientOrigin) gpuIf->FreeMemory(dGradientOrigin); @@ -95,6 +101,17 @@ int BeagleGPUSpectralImpl::createInstance( dSpectralDistances[i] = gpuIf->CreateSubPointer(dSpectralDistancesOrigin, distStride * i, distStride); } + // Gather-batch-scatter queue for the per-branch distance upload in + // updateTransitionMatrices (see BeagleGPUSpectralImpl.h for rationale). + // Sized for the largest possible single call: every matrix, every category. + int spectralQueueLength = this->kMatrixCount * this->kCategoryCount; + dSpectralPtrQueue = gpuIf->AllocateMemory(sizeof(unsigned int) * spectralQueueLength); + hSpectralPtrQueue = (unsigned int*) gpuIf->MallocHost(sizeof(unsigned int) * spectralQueueLength); + checkHostMemory(hSpectralPtrQueue); + dSpectralDistanceQueue = gpuIf->AllocateMemory(sizeof(Real) * spectralQueueLength); + hSpectralDistanceQueue = (Real*) gpuIf->MallocHost(sizeof(Real) * spectralQueueLength); + checkHostMemory(hSpectralDistanceQueue); + // Backward eigenvector arrays: one S*S block per eigen decomposition. int S = this->kPaddedStateCount; size_t matStride = gpuIf->AlignMemOffset((size_t)S * S * sizeof(Real)); @@ -136,17 +153,21 @@ int BeagleGPUSpectralImpl::updateTransitionMatrices( if (count > 0) { const double* categoryRates = this->hCategoryRates[0]; int nCat = this->kCategoryCount; - Real* hDist = (Real*) malloc(sizeof(Real) * nCat); GPUInterface* gpuIf = this->gpu; + int totalCount = 0; for (int i = 0; i < count; i++) { int matIdx = probabilityIndices[i]; hEigenIndexForMatrix[matIdx] = eigenIndex; for (int j = 0; j < nCat; j++) { - hDist[j] = (Real)(edgeLengths[i] * categoryRates[j]); + hSpectralPtrQueue[totalCount] = matIdx * kSpectralDistanceStrideElements + j; + hSpectralDistanceQueue[totalCount] = (Real)(edgeLengths[i] * categoryRates[j]); + totalCount++; } - gpuIf->MemcpyHostToDevice(dSpectralDistances[matIdx], hDist, sizeof(Real) * nCat); } - free(hDist); + gpuIf->MemcpyHostToDevice(dSpectralPtrQueue, hSpectralPtrQueue, sizeof(unsigned int) * totalCount); + gpuIf->MemcpyHostToDevice(dSpectralDistanceQueue, hSpectralDistanceQueue, sizeof(Real) * totalCount); + this->kernels->ScatterSpectralDistances(dSpectralDistancesOrigin, dSpectralPtrQueue, + dSpectralDistanceQueue, totalCount); } return BEAGLE_SUCCESS; } diff --git a/libhmsbeagle/GPU/KernelLauncher.cpp b/libhmsbeagle/GPU/KernelLauncher.cpp index e62bcc49..6adba33c 100644 --- a/libhmsbeagle/GPU/KernelLauncher.cpp +++ b/libhmsbeagle/GPU/KernelLauncher.cpp @@ -318,6 +318,8 @@ void KernelLauncher::LoadKernels() { } else { fAdjointMergedN = gpu->GetFunction("kernelAdjointMergedN"); } + + fScatterSpectralDistances = gpu->GetFunction("kernelScatterSpectralDistances"); } fPartialsPartialsEdgeFirstDerivatives = gpu->GetFunction( @@ -1184,6 +1186,27 @@ void KernelLauncher::AdjointCrossProductMergedSpectralN(GPUPtr partialsOrigin, gpu->SynchronizeDevice(); } +void KernelLauncher::ScatterSpectralDistances(GPUPtr distOrigin, + GPUPtr destOffsets, + GPUPtr values, + unsigned int totalCount) { +#ifdef BEAGLE_DEBUG_FLOW + fprintf(stderr, "\t\tEntering KernelLauncher::ScatterSpectralDistances\n"); +#endif + + const unsigned int blockSize = 128; + Dim3Int block(blockSize); + Dim3Int grid((totalCount + blockSize - 1) / blockSize); + + gpu->LaunchKernel(fScatterSpectralDistances, block, grid, 3, 4, + distOrigin, destOffsets, values, + totalCount); + +#ifdef BEAGLE_DEBUG_FLOW + fprintf(stderr, "\t\tLeaving KernelLauncher::ScatterSpectralDistances\n"); +#endif +} + void KernelLauncher::PartialsPartialsPruningDynamicCheckScaling(GPUPtr partials1, GPUPtr partials2, diff --git a/libhmsbeagle/GPU/KernelLauncher.h b/libhmsbeagle/GPU/KernelLauncher.h index d8af125b..0a3666fa 100644 --- a/libhmsbeagle/GPU/KernelLauncher.h +++ b/libhmsbeagle/GPU/KernelLauncher.h @@ -71,6 +71,12 @@ class KernelLauncher { GPUFunction fAdjointMergedN; /* Adjoint cross-product kernel — 4-state, merged single launch */ GPUFunction fAdjointMerged4; + /* Scatter kernel for the batched spectral-distance queue: one thread per + * (branch, category) queue entry, writes value to a device-computed + * offset. State-count-independent, so the same function name resolves + * regardless of which spectral kernel program (4-state or generic-N) + * was compiled. */ + GPUFunction fScatterSpectralDistances; GPUFunction fStatesPartialsByPatternBlockCoherentMulti; GPUFunction fStatesPartialsByPatternBlockCoherentPartition; GPUFunction fStatesPartialsByPatternBlockCoherent; @@ -473,6 +479,15 @@ class KernelLauncher { unsigned int categoryCount, int branchCount); + /* Batched scatter of per-branch spectral distances: destOffsets[k]/ + * values[k] pairs (one per (branch, category) queue entry) are written + * to distOrigin[destOffsets[k]] in a single launch, replacing what was + * previously one MemcpyHostToDevice call per branch. */ + void ScatterSpectralDistances(GPUPtr distOrigin, + GPUPtr destOffsets, + GPUPtr values, + unsigned int totalCount); + void PartialsStatesEdgeFirstDerivatives(GPUPtr out, GPUPtr states0, GPUPtr partials0, diff --git a/libhmsbeagle/GPU/kernels/kernelsSpectralIfDef.cu b/libhmsbeagle/GPU/kernels/kernelsSpectralIfDef.cu index dccb322b..038d1291 100644 --- a/libhmsbeagle/GPU/kernels/kernelsSpectralIfDef.cu +++ b/libhmsbeagle/GPU/kernels/kernelsSpectralIfDef.cu @@ -1326,6 +1326,24 @@ KW_GLOBAL_KERNEL void kernelAdjointMergedN( } } +/* ═══════════════════════════════════════════════════════════════════════════ + * Batched scatter for updateTransitionMatrices' per-branch distance queue. + * One thread per (branch, category) queue entry: replaces what used to be + * one MemcpyHostToDevice call per branch with two flat uploads (destOffsets, + * values) plus this single on-device scatter. No PADDED_STATE_COUNT/ + * STATE_COUNT dependence, so this is byte-identical in kernelsSpectralIfDef4.cu. + * ═══════════════════════════════════════════════════════════════════════════*/ +KW_GLOBAL_KERNEL void kernelScatterSpectralDistances( + KW_GLOBAL_VAR REAL* KW_RESTRICT distOrigin, + KW_GLOBAL_VAR unsigned int* KW_RESTRICT destOffsets, + KW_GLOBAL_VAR REAL* KW_RESTRICT values, + int totalCount) { + int idx = KW_GROUP_ID_0 * KW_LOCAL_SIZE_0 + KW_LOCAL_ID_0; + if (idx < totalCount) { + distOrigin[destOffsets[idx]] = values[idx]; + } +} + #ifdef CUDA } /* extern "C" */ #endif diff --git a/libhmsbeagle/GPU/kernels/kernelsSpectralIfDef4.cu b/libhmsbeagle/GPU/kernels/kernelsSpectralIfDef4.cu index af212c3f..2cf39f7d 100644 --- a/libhmsbeagle/GPU/kernels/kernelsSpectralIfDef4.cu +++ b/libhmsbeagle/GPU/kernels/kernelsSpectralIfDef4.cu @@ -1265,6 +1265,25 @@ KW_GLOBAL_KERNEL void kernelAdjointMerged4( } } +/* ═══════════════════════════════════════════════════════════════════════════ + * Batched scatter for updateTransitionMatrices' per-branch distance queue. + * One thread per (branch, category) queue entry: replaces what used to be + * one MemcpyHostToDevice call per branch with two flat uploads (destOffsets, + * values) plus this single on-device scatter. No PADDED_STATE_COUNT/ + * STATE_COUNT dependence, so this is byte-identical to kernelsSpectralIfDef.cu's + * copy. + * ═══════════════════════════════════════════════════════════════════════════*/ +KW_GLOBAL_KERNEL void kernelScatterSpectralDistances( + KW_GLOBAL_VAR REAL* KW_RESTRICT distOrigin, + KW_GLOBAL_VAR unsigned int* KW_RESTRICT destOffsets, + KW_GLOBAL_VAR REAL* KW_RESTRICT values, + int totalCount) { + int idx = KW_GROUP_ID_0 * KW_LOCAL_SIZE_0 + KW_LOCAL_ID_0; + if (idx < totalCount) { + distOrigin[destOffsets[idx]] = values[idx]; + } +} + #ifdef CUDA } /* extern "C" */ #endif From 9cf962616fe505e4945cb834d72d9e881be3e765 Mon Sep 17 00:00:00 2001 From: fil-monti Date: Sat, 11 Jul 2026 00:27:01 -0700 Subject: [PATCH 09/13] Extend BEAGLE_PROFILE_SUMMARY to cover MemcpyDeviceToHost and MemsetZero Follow-up to the distance-batching fix: after that fix, only ~25.3s of a 171.4s real-analysis BEAST run was accounted for by the four originally profiled metrics (LaunchKernel, LaunchKernelConcurrent, SynchronizeHost, MemcpyHostToDevice) -- an unexplained ~148s (85%) gap that couldn't be generic JVM overhead, since the CPU-backed run of the identical analysis takes only ~16.8s total. Adding MemcpyDeviceToHost instrumentation (same accumulating-counter, warmup-gated, atexit-printed style as the existing metrics) found it accounts for 91.9s (53.6% of wall time) on its own -- bigger than the entire MemcpyHostToDevice bottleneck was before its fix. Adding MemsetZero instrumentation found a further 35-37s, previously assumed "likely small" and never actually measured. Together these bring total accounted time to ~90% of wall time, closing the mystery gap from ~53s down to ~16s. Both new metrics were then diagnosed (via a temporary, since-reverted SynchronizeHost() probe inserted at each call site -- not part of this commit) to be genuine, unavoidable GPU-kernel-completion wait time, not fixable memcpy/write-level overhead: OpenCL kernel launches are asynchronous, so the real compute time doesn't show up at the kernel enqueue call, it surfaces at whatever blocking call happens next. No software fix exists at this level; the remaining GPU-vs-CPU gap for this problem size is now understood to be substantially explained by real compute time, not a bug. Full investigation and conclusion in beagle-bugs/gpu-spectral-device-to-host-readback/PLAN.md and PROJECT_STATUS.md. MemsetZero's existing blocking clEnqueueWriteBuffer mechanism (deliberately not clEnqueueFillBuffer, per the pre-existing comment documenting a real, previously-observed intermittent data race with that alternative) is left completely unchanged -- only a timing wrapper was added around it. --- libhmsbeagle/GPU/GPUInterfaceOpenCL.cpp | 43 ++++++++++++++++++++++++- 1 file changed, 42 insertions(+), 1 deletion(-) diff --git a/libhmsbeagle/GPU/GPUInterfaceOpenCL.cpp b/libhmsbeagle/GPU/GPUInterfaceOpenCL.cpp index 0543a1b4..eae6124d 100644 --- a/libhmsbeagle/GPU/GPUInterfaceOpenCL.cpp +++ b/libhmsbeagle/GPU/GPUInterfaceOpenCL.cpp @@ -60,6 +60,12 @@ struct BeagleProfileStats { long memcpyH2DCalls = 0; double memcpyH2DSeconds = 0.0; size_t memcpyH2DBytes = 0; + long memcpyD2HCalls = 0; + double memcpyD2HSeconds = 0.0; + size_t memcpyD2HBytes = 0; + long memsetZeroCalls = 0; + double memsetZeroSeconds = 0.0; + size_t memsetZeroBytes = 0; }; BeagleProfileStats gBeagleProfile; @@ -68,7 +74,8 @@ void beagleProfilePrintSummary() { if (!gBeagleProfile.enabled) return; BeagleProfileStats& p = gBeagleProfile; double totalAccounted = p.launchKernelSeconds + p.launchKernelConcurrentSeconds + - p.syncHostSeconds + p.memcpyH2DSeconds; + p.syncHostSeconds + p.memcpyH2DSeconds + p.memcpyD2HSeconds + + p.memsetZeroSeconds; fprintf(stderr, "\n===== BEAGLE_PROFILE_SUMMARY (GPUInterfaceOpenCL) =====\n"); fprintf(stderr, "warmup discarded: first %ld SynchronizeHost calls\n", BEAGLE_PROFILE_WARMUP_SYNC); fprintf(stderr, "LaunchKernel : calls=%-8ld total=%.6fs avg=%.3fus\n", @@ -84,6 +91,14 @@ void beagleProfilePrintSummary() { p.memcpyH2DCalls, p.memcpyH2DSeconds, p.memcpyH2DCalls ? p.memcpyH2DSeconds * 1e6 / p.memcpyH2DCalls : 0.0, p.memcpyH2DBytes); + fprintf(stderr, "MemcpyDeviceToHost : calls=%-8ld total=%.6fs avg=%.3fus bytes=%zu\n", + p.memcpyD2HCalls, p.memcpyD2HSeconds, + p.memcpyD2HCalls ? p.memcpyD2HSeconds * 1e6 / p.memcpyD2HCalls : 0.0, + p.memcpyD2HBytes); + fprintf(stderr, "MemsetZero : calls=%-8ld total=%.6fs avg=%.3fus bytes=%zu\n", + p.memsetZeroCalls, p.memsetZeroSeconds, + p.memsetZeroCalls ? p.memsetZeroSeconds * 1e6 / p.memsetZeroCalls : 0.0, + p.memsetZeroBytes); fprintf(stderr, "Total host-side time accounted above (post-warmup): %.6fs\n", totalAccounted); fprintf(stderr, "===== END BEAGLE_PROFILE_SUMMARY =====\n\n"); fflush(stderr); @@ -1123,6 +1138,10 @@ void GPUInterface::MemsetZero(GPUPtr dest, fprintf(stderr, "\t\t\tEntering GPUInterface::MemsetZero\n"); #endif + bool prof = beagleProfilingEnabled(); + std::chrono::high_resolution_clock::time_point t0; + if (prof) t0 = std::chrono::high_resolution_clock::now(); + /* Deliberately not clEnqueueFillBuffer: on this OpenCL implementation * (observed on Apple M1 Max), a kernel launched afterward on the same * in-order queue intermittently read stale/garbage data from a @@ -1138,6 +1157,14 @@ void GPUInterface::MemsetZero(GPUPtr dest, NULL, NULL)); free(hZero); + if (prof && gBeagleProfile.warm) { + double dt = std::chrono::duration( + std::chrono::high_resolution_clock::now() - t0).count(); + gBeagleProfile.memsetZeroCalls++; + gBeagleProfile.memsetZeroSeconds += dt; + gBeagleProfile.memsetZeroBytes += memSize; + } + #ifdef BEAGLE_DEBUG_FLOW fprintf(stderr, "\t\t\tLeaving GPUInterface::MemsetZero\n"); #endif @@ -1180,9 +1207,23 @@ void GPUInterface::MemcpyDeviceToHost(void* dest, fprintf(stderr, "\t\t\tEntering GPUInterface::MemcpyDeviceToHost\n"); #endif + bool prof = beagleProfilingEnabled(); + std::chrono::high_resolution_clock::time_point t0; + if (prof) t0 = std::chrono::high_resolution_clock::now(); + + // CL_TRUE => blocking read: this call does not return until the transfer + // completes, regardless of payload size. SAFE_CL(clEnqueueReadBuffer(openClCommandQueues[0], src, CL_TRUE, 0, memSize, dest, 0, NULL, NULL)); + if (prof && gBeagleProfile.warm) { + double dt = std::chrono::duration( + std::chrono::high_resolution_clock::now() - t0).count(); + gBeagleProfile.memcpyD2HCalls++; + gBeagleProfile.memcpyD2HSeconds += dt; + gBeagleProfile.memcpyD2HBytes += memSize; + } + #ifdef BEAGLE_DEBUG_FLOW fprintf(stderr, "\t\t\tLeaving GPUInterface::MemcpyDeviceToHost\n"); #endif From 56e3239a774ab508c9d7f550d57764d83b2d66f0 Mon Sep 17 00:00:00 2001 From: fil-monti Date: Sun, 12 Jul 2026 08:32:44 -0700 Subject: [PATCH 10/13] Add BeagleOpenCLSpectral_kernels.h to .gitignore MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Auto-generated build artifact (same category as the already-ignored BeagleOpenCL_kernels.h) — the .gitignore had a gap for the newer spectral-kernel filename. --- .gitignore | 1 + 1 file changed, 1 insertion(+) diff --git a/.gitignore b/.gitignore index c00d9ba5..2dc47505 100644 --- a/.gitignore +++ b/.gitignore @@ -209,6 +209,7 @@ examples/tinytest/.libs/ examples/tinytest/tinytest.o libhmsbeagle/GPU/kernels/BeagleOpenCL_kernels.h libhmsbeagle/GPU/kernels/BeagleOpenCL_kernels_xcode.h +libhmsbeagle/GPU/kernels/BeagleOpenCLSpectral_kernels.h libhmsbeagle/GPU/libhmsbeagle-opencl.la libhmsbeagle/GPU/libhmsbeagle_opencl_la-GPUImplHelper.lo libhmsbeagle/GPU/libhmsbeagle_opencl_la-GPUInterfaceOpenCL.lo From b8090e5a860a9423a9559138cc2dcda43cf9a328 Mon Sep 17 00:00:00 2001 From: fil-monti Date: Sun, 12 Jul 2026 08:34:00 -0700 Subject: [PATCH 11/13] GPU spectral: skip redundant dense transition-matrix construction BeagleGPUSpectralImpl::updateTransitionMatrices unconditionally built a full dense P(t) matrix via the base class's construction on every probability-only branch-length update, even though spectral pruning/gradient never reads it (dispatchPrunePP/SP/SS only touch dIevc/dEvec/dEigenValues/dSpectralDistances). Not simply removable, though -- getTransitionMatrix() is real, live API used elsewhere (e.g. DiscreteTraitNodeHeightDelegate.java), so there's no way to know in advance whether some caller will need it. Fixed via lazy materialization: cache eigenIndex/edgeLength per matrix index instead of building eagerly, and add a getTransitionMatrix override that rebuilds on first actual access by calling the unmodified base-class construction for just that one index. No kernel code touched. Verified via a bit-identical trajectory check (config A never calls getTransitionMatrix, so this must be a pure no-op for its results) plus a standalone test built to exercise the new lazy-rebuild path directly, comparing output byte-for-byte against the unmodified eager construction across several edge cases -- identical in every case. ~3% wall-clock improvement on the main benchmark, larger (~20-27%) on workloads that call updateTransitionMatrices more frequently. Full writeup: beagle-bugs/gpu-spectral-skip-redundant-dense-matrix/. --- libhmsbeagle/GPU/BeagleGPUSpectralImpl.h | 16 +++++++ libhmsbeagle/GPU/BeagleGPUSpectralImpl.hpp | 49 +++++++++++++++++++--- 2 files changed, 60 insertions(+), 5 deletions(-) diff --git a/libhmsbeagle/GPU/BeagleGPUSpectralImpl.h b/libhmsbeagle/GPU/BeagleGPUSpectralImpl.h index 607cd72b..b33e25a5 100644 --- a/libhmsbeagle/GPU/BeagleGPUSpectralImpl.h +++ b/libhmsbeagle/GPU/BeagleGPUSpectralImpl.h @@ -38,6 +38,19 @@ class BeagleGPUSpectralImpl : public BeagleGPUImpl { GPUPtr dSpectralDistancesOrigin; GPUPtr* dSpectralDistances; int* hEigenIndexForMatrix; + + /* Lazy dense-matrix materialization (see PLAN.md in + * beagle-bugs/gpu-spectral-skip-redundant-dense-matrix/): probability-only + * updateTransitionMatrices calls no longer eagerly build the dense + * kStateCount x kStateCount transition matrix, since the spectral + * pruning/gradient paths never read it. Instead the eigenIndex/edgeLength + * needed to rebuild it are cached here, and the dense matrix is + * materialized on first actual getTransitionMatrix() call. All sized + * kMatrixCount. */ + bool* hDenseMatrixValid; + int* hPendingDenseEigenIndex; + double* hPendingDenseEdgeLength; + /* 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 @@ -125,6 +138,9 @@ class BeagleGPUSpectralImpl : public BeagleGPUImpl { const double* edgeLengths, int count) override; + int getTransitionMatrix(int matrixIndex, + double* outMatrix) override; + int updatePrePartials(const int* operations, int operationCount, int cumulativeScaleIndex, diff --git a/libhmsbeagle/GPU/BeagleGPUSpectralImpl.hpp b/libhmsbeagle/GPU/BeagleGPUSpectralImpl.hpp index 34bf34d0..dfc2c94c 100644 --- a/libhmsbeagle/GPU/BeagleGPUSpectralImpl.hpp +++ b/libhmsbeagle/GPU/BeagleGPUSpectralImpl.hpp @@ -34,6 +34,7 @@ namespace gpu { BEAGLE_GPU_TEMPLATE BeagleGPUSpectralImpl::BeagleGPUSpectralImpl() : dSpectralDistancesOrigin(0), dSpectralDistances(NULL), hEigenIndexForMatrix(NULL), + hDenseMatrixValid(NULL), hPendingDenseEigenIndex(NULL), hPendingDenseEdgeLength(NULL), dSpectralPtrQueue(0), dSpectralDistanceQueue(0), hSpectralPtrQueue(NULL), hSpectralDistanceQueue(NULL), dEvecTOrigin(0), dEvecT(NULL), dIevcTOrigin(0), dIevcT(NULL), @@ -62,6 +63,9 @@ BeagleGPUSpectralImpl::~BeagleGPUSpectralImpl() { free(dIevcT); free(dGradient); free(hEigenIndexForMatrix); + free(hDenseMatrixValid); + free(hPendingDenseEigenIndex); + free(hPendingDenseEdgeLength); free(hEigenDecompIsAllReal); } @@ -91,6 +95,10 @@ int BeagleGPUSpectralImpl::createInstance( hEigenIndexForMatrix = (int*) calloc(this->kMatrixCount, sizeof(int)); + hDenseMatrixValid = (bool*) calloc(this->kMatrixCount, sizeof(bool)); // false = stale/never built + hPendingDenseEigenIndex = (int*) calloc(this->kMatrixCount, sizeof(int)); + hPendingDenseEdgeLength = (double*) calloc(this->kMatrixCount, sizeof(double)); + GPUInterface* gpuIf = this->gpu; dSpectralDistances = (GPUPtr*) malloc(sizeof(GPUPtr) * this->kMatrixCount); @@ -144,11 +152,26 @@ int BeagleGPUSpectralImpl::updateTransitionMatrices( const int* secondDerivativeIndices, const double* edgeLengths, int count) { - int rc = BeagleGPUImpl::updateTransitionMatrices(eigenIndex, probabilityIndices, - firstDerivativeIndices, - secondDerivativeIndices, - edgeLengths, count); - if (rc != BEAGLE_SUCCESS || firstDerivativeIndices != NULL) return rc; + // Only defer the dense (probability-only) case. Derivative requests are rare in this + // project's usage and the existing eager-build behavior for them is left completely + // untouched, to keep this change minimal and scoped to the one case that's actually + // wasteful (see PLAN.md in beagle-bugs/gpu-spectral-skip-redundant-dense-matrix/). + if (firstDerivativeIndices != NULL) { + return BeagleGPUImpl::updateTransitionMatrices(eigenIndex, probabilityIndices, + firstDerivativeIndices, + secondDerivativeIndices, + edgeLengths, count); + } + + // Probability-only update: DO NOT build the dense matrix now. Cache what's needed to + // build it lazily (see getTransitionMatrix override below), in case some caller asks + // for it later. Spectral pruning/gradient paths never read dMatrices. + for (int i = 0; i < count; i++) { + int matIdx = probabilityIndices[i]; + hPendingDenseEigenIndex[matIdx] = eigenIndex; + hPendingDenseEdgeLength[matIdx] = edgeLengths[i]; + hDenseMatrixValid[matIdx] = false; + } if (count > 0) { const double* categoryRates = this->hCategoryRates[0]; @@ -172,6 +195,22 @@ int BeagleGPUSpectralImpl::updateTransitionMatrices( return BEAGLE_SUCCESS; } +BEAGLE_GPU_TEMPLATE +int BeagleGPUSpectralImpl::getTransitionMatrix(int matrixIndex, + double* outMatrix) { + if (!hDenseMatrixValid[matrixIndex]) { + // Materialize on demand, reusing the existing, already-correct base-class + // construction unchanged -- just deferred to first actual use instead of built + // eagerly and thrown away. + int rc = BeagleGPUImpl::updateTransitionMatrices( + hPendingDenseEigenIndex[matrixIndex], &matrixIndex, + NULL, NULL, &hPendingDenseEdgeLength[matrixIndex], 1); + if (rc != BEAGLE_SUCCESS) return rc; + hDenseMatrixValid[matrixIndex] = true; + } + return BeagleGPUImpl::getTransitionMatrix(matrixIndex, outMatrix); +} + BEAGLE_GPU_TEMPLATE void BeagleGPUSpectralImpl::dispatchPrunePP( GPUPtr p1, GPUPtr p2, GPUPtr p3, From 858a5825059dd5d545c7ea30de0eee55370535e1 Mon Sep 17 00:00:00 2001 From: fil-monti Date: Sun, 12 Jul 2026 08:35:28 -0700 Subject: [PATCH 12/13] GPU spectral: fix misaligned adjoint-queue/MemsetZero offsets for multi-eigendecomposition models dEvecTOrigin/dIevcTOrigin/dGradientOrigin are allocated with a device-aligned per-eigendecomposition stride (AlignMemOffset does real, non-trivial rounding on non-NVIDIA/non-Apple platforms, confirmed matStride=128 bytes vs logical S*S*sizeof(Real)=64 bytes for S=4 on this AMD GPU), but the adjoint queue-building code in calculateAdjointCrossProducts indexed into them using the raw unaligned S*S for two of nine fields (rec[3], rec[7]) and for the MemsetZero call size, while three sibling fields in the same block correctly use a dedicated aligned-stride accessor/member. Invisible in every existing benchmark, which all use exactly one eigendecomposition. Fixed by adding kSpectralMatrixStrideElements/Bytes and using it for rec[3]/rec[7]/the MemsetZero size instead of the raw S*S. Verified with a new standalone test (adjointtest4 --multieigen): creates 2 distinct eigendecompositions, binds different branches to each, compares GPU gradient output against an independent CPU reference. Pre-fix, eigenIndex=1's GPU gradient came back as uniform garbage against a correct, structured CPU reference. Post-fix, both eigendecompositions match the CPU reference closely. Zero regression on the existing single-eigendecomposition benchmark and no performance change, as expected for a pure address-computation fix. Full writeup: beagle-bugs/gpu-spectral-misaligned-adjoint-queue-stride/. --- examples/hmctest/adjointtest4.cpp | 185 ++++++++++++++++++++- libhmsbeagle/GPU/BeagleGPUSpectralImpl.h | 9 + libhmsbeagle/GPU/BeagleGPUSpectralImpl.hpp | 13 +- 3 files changed, 200 insertions(+), 7 deletions(-) diff --git a/examples/hmctest/adjointtest4.cpp b/examples/hmctest/adjointtest4.cpp index 55532a6f..b3d2260f 100644 --- a/examples/hmctest/adjointtest4.cpp +++ b/examples/hmctest/adjointtest4.cpp @@ -207,7 +207,8 @@ static void printResources() { static int createInstance(bool useGpu, int whichDevice, bool singlePrec, int stateCount, int nPatterns, int nCats, - bool complexEigen = false) { + bool complexEigen = false, + int nEigenDecomp = 1, int nTransitionMatrices = 4) { long prefFlags = BEAGLE_FLAG_SCALERS_RAW; prefFlags |= useGpu ? BEAGLE_FLAG_PROCESSOR_GPU : BEAGLE_FLAG_PROCESSOR_CPU; prefFlags |= singlePrec ? BEAGLE_FLAG_PRECISION_SINGLE : BEAGLE_FLAG_PRECISION_DOUBLE; @@ -223,8 +224,8 @@ static int createInstance(bool useGpu, int whichDevice, bool singlePrec, 3, /* nCompact */ stateCount, nPatterns, - 1, /* nEigenDecomp */ - 4, /* nTransitionMatrices */ + nEigenDecomp, + nTransitionMatrices, nCats, 0, /* nScaleBuffers */ whichDevice >= 0 ? &whichDevice : NULL, @@ -1251,6 +1252,177 @@ static int runTestDynamicTransition(int gpuDevice) { return allOk ? 0 : 1; } +/* ═══════════════════════════════════════════════════════════════════════ + * Multi-eigendecomposition test (S=4, nEigenDecomp=2) + * + * Regression test for the misaligned adjoint-queue/MemsetZero stride bug + * (beagle-bugs/gpu-spectral-misaligned-adjoint-queue-stride/PLAN.md). + * BeagleGPUSpectralImpl allocates dEvecTOrigin/dIevcTOrigin/dGradientOrigin + * with a per-eigendecomposition stride that is *aligned* + * (AlignMemOffset(S*S*sizeof(Real))), which device alignment rules can make + * strictly larger than the logical S*S*sizeof(Real). The queue-building code + * in calculateAdjointCrossProducts (rec[3]/rec[7]) and the MemsetZero call + * that clears the gradient accumulator must use that same aligned stride to + * index into these pooled buffers by eigendecomposition index -- using the + * raw unaligned S*S instead is invisible whenever eigenIndex==0 (multiplied + * away by ei*anything==0), which is true of every other test/benchmark in + * this project. This test creates a real nEigenDecomp=2 instance and runs + * the SAME 4-edge tree/topology twice on the SAME 4 transition-matrix slots: + * once with all edges bound to eigenIndex=0 (JC69 -- expected to pass + * regardless of the bug, since ei=0 can never expose it) and once with all + * edges rebound to eigenIndex=1 (a distinct, non-degenerate real + * eigendecomposition -- exercises the exact rec[3]/rec[7]/MemsetZero code + * path this bug lives in, since ei=1 is where "ei*SS" and "ei*matStride" + * first diverge whenever matStride>SS*sizeof(Real) on this hardware). + * GPU output for each step is compared against an independently-computed + * CPU (double-precision) reference for the same step -- the CPU + * implementation has no analogous aligned/unaligned pooled-buffer stride + * concept, so it is unaffected by this bug and serves as ground truth. + * ═══════════════════════════════════════════════════════════════════════*/ + +struct MultiEigenResult { + double logLA, logLB; + std::vector gradA, gradB; /* S*S each */ +}; + +static MultiEigenResult runMultiEigenOne(bool useGpu, int dev, bool singlePrec, + const std::vector& evecB, + const std::vector& ivecB, + const std::vector& evalB) { + const int S = 4, nPat = 4, nCats = 2; + MultiEigenResult R; + + int inst = createInstance(useGpu, dev, singlePrec, S, nPat, nCats, + /*complexEigen=*/false, /*nEigenDecomp=*/2); + if (inst < 0) { R.logLA = R.logLB = NAN; return R; } + + int *hSt = dnaToStates(human); + int *cSt = dnaToStates(chimp); + int *gSt = dnaToStates(gorilla); + beagleSetTipStates(inst, 0, hSt); + beagleSetTipStates(inst, 1, cSt); + beagleSetTipStates(inst, 2, gSt); + free(hSt); free(cSt); free(gSt); + + beagleSetCategoryRates(inst, rates2); + std::vector pw(nPat, 1.0); + beagleSetPatternWeights(inst, pw.data()); + beagleSetStateFrequencies(inst, 0, freqs4); + beagleSetCategoryWeights(inst, 0, catWeights2); + + int nodeIdx[4] = { 0, 1, 2, 3 }; + BeagleOperation postOps[2] = { + { 3, BEAGLE_OP_NONE, BEAGLE_OP_NONE, 0, 0, 1, 1 }, + { 4, BEAGLE_OP_NONE, BEAGLE_OP_NONE, 2, 2, 3, 3 } + }; + std::vector rootPre((size_t)nCats * nPat * S); + for (int c = 0; c < nCats; c++) + for (int p = 0; p < nPat; p++) + for (int s = 0; s < S; s++) + rootPre[(size_t)c * nPat * S + p * S + s] = freqs4[s]; + BeagleOperation preOps[4] = { + { 6, BEAGLE_OP_NONE, BEAGLE_OP_NONE, 5, BEAGLE_OP_NONE, 2, 2 }, + { 7, BEAGLE_OP_NONE, BEAGLE_OP_NONE, 5, BEAGLE_OP_NONE, 3, 3 }, + { 8, BEAGLE_OP_NONE, BEAGLE_OP_NONE, 6, BEAGLE_OP_NONE, 0, 0 }, + { 9, BEAGLE_OP_NONE, BEAGLE_OP_NONE, 6, BEAGLE_OP_NONE, 1, 1 } + }; + BeagleBranchOperation branchOps[4] = { + { 1, 8, 1, 0 }, + { 0, 9, 0, 0 }, + { 2, 7, 2, 0 }, + { 3, 6, 3, 0 } + }; + + auto runStep = [&](int eigenIndex, const double* evec, const double* ivec, + const double* eval) -> std::pair> { + beagleSetEigenDecomposition(inst, eigenIndex, evec, ivec, eval); + beagleUpdateTransitionMatrices(inst, eigenIndex, nodeIdx, NULL, NULL, edgeLengths, 4); + beagleUpdatePartials(inst, postOps, 2, BEAGLE_OP_NONE); + double logL = 0.0; + int rootIdx = 4, cwIdx = 0, sfIdx = 0, csi = BEAGLE_OP_NONE; + beagleCalculateRootLogLikelihoods(inst, &rootIdx, &cwIdx, &sfIdx, &csi, 1, &logL); + beagleSetPartials(inst, 5, rootPre.data()); + beagleUpdatePrePartials_v5(inst, preOps, 4, BEAGLE_OP_NONE, BEAGLE_PARTIALS_TOP); + std::vector grad(S * S, 0.0); + beagleCalculateAdjointDerivative(inst, branchOps, 0, 0, 4, 0, 4, grad.data(), NULL); + return { logL, grad }; + }; + + /* Step A: eigenIndex=0 (JC69) -- ei==0, bug is provably invisible here. */ + auto a = runStep(0, jcEvec, jcIvec, jcEval); + R.logLA = a.first; R.gradA = a.second; + + /* Step B: eigenIndex=1 (distinct real model) -- exercises the bug. */ + auto b = runStep(1, evecB.data(), ivecB.data(), evalB.data()); + R.logLB = b.first; R.gradB = b.second; + + beagleFinalizeInstance(inst); + return R; +} + +static int runTestMultiEigen(int gpuDevice) { + const int S = 4; + const double tol = 1e-4; + + fprintf(stdout, "\n=== Multi-eigendecomposition test (S=4, nEigenDecomp=2) ===\n"); + setenv("BEAGLE_DEBUG_STRIDE", "1", 1); + + /* Model B: symmetric (r_fwd=r_bkd=1.0) 4-state circulant -- exactly + * real eigenvalues, numerically distinct from JC69. */ + std::vector evecB, ivecB, evalB; + buildCirculantN(4, evecB, ivecB, evalB, /*r_fwd=*/1.0, /*r_bkd=*/1.0); + + fprintf(stdout, "\nCreating CPU (double) instance (nEigenDecomp=2):\n"); + MultiEigenResult cpuR = runMultiEigenOne(false, -1, false, evecB, ivecB, evalB); + if (std::isnan(cpuR.logLA)) { fprintf(stderr, "Failed CPU instance\n"); return 1; } + fprintf(stdout, " CPU step A (eigenIndex=0, JC69): logL=%.6f\n", cpuR.logLA); + fprintf(stdout, " CPU step B (eigenIndex=1, model B): logL=%.6f\n", cpuR.logLB); + + if (gpuDevice < 0) { + fprintf(stdout, "(no --gpu device given; CPU-only run, nothing to compare)\n"); + return 0; + } + + fprintf(stdout, "\nCreating GPU (single) instance (nEigenDecomp=2):\n"); + MultiEigenResult gpuR = runMultiEigenOne(true, gpuDevice, true, evecB, ivecB, evalB); + if (std::isnan(gpuR.logLA)) { fprintf(stderr, "Failed GPU instance\n"); return 1; } + fprintf(stdout, " GPU step A (eigenIndex=0, JC69): logL=%.6f\n", gpuR.logLA); + fprintf(stdout, " GPU step B (eigenIndex=1, model B): logL=%.6f\n", gpuR.logLB); + + auto compareGrad = [&](const char* label, const std::vector& cpu, + const std::vector& gpu) -> bool { + double maxDiff = 0.0; + for (int k = 0; k < S * S; k++) { + double d = fabs(cpu[k] - gpu[k]); + if (d > maxDiff || std::isnan(d)) maxDiff = d; + } + bool ok = (maxDiff <= tol) && !std::isnan(maxDiff); + fprintf(stdout, " %-42s max|CPU-GPU diff|=%.3e (tol=%.0e) %s\n", + label, maxDiff, tol, ok ? "PASS" : "FAIL"); + if (!ok) { + fprintf(stdout, " CPU grad:\n"); + for (int i = 0; i < S; i++) { + fprintf(stdout, " "); + for (int j = 0; j < S; j++) fprintf(stdout, " %10.6f", cpu[i*S+j]); + fprintf(stdout, "\n"); + } + fprintf(stdout, " GPU grad:\n"); + for (int i = 0; i < S; i++) { + fprintf(stdout, " "); + for (int j = 0; j < S; j++) fprintf(stdout, " %10.6f", gpu[i*S+j]); + fprintf(stdout, "\n"); + } + } + return ok; + }; + + bool okA = compareGrad("Step A gradient (eigenIndex=0, JC69)", cpuR.gradA, gpuR.gradA); + bool okB = compareGrad("Step B gradient (eigenIndex=1, model B)", cpuR.gradB, gpuR.gradB); + + fprintf(stdout, "Multi-eigendecomposition test: %s\n", (okA && okB) ? "PASS" : "FAIL"); + return (okA && okB) ? 0 : 1; +} + /* ── main ────────────────────────────────────────────────────────────── */ int main(int argc, const char *argv[]) { @@ -1263,6 +1435,7 @@ int main(int argc, const char *argv[]) { bool run17 = false; bool runRealOnly = false; bool runDynamic = false; + bool runMultiEigen = false; for (int i = 1; i < argc; i++) { if (strcmp(argv[i], "--gpu") == 0 && i+1 < argc) { @@ -1279,8 +1452,11 @@ int main(int argc, const char *argv[]) { } else if (strcmp(argv[i], "--dynamic") == 0) { run4 = false; runDynamic = true; + } else if (strcmp(argv[i], "--multieigen") == 0) { + run4 = false; + runMultiEigen = true; } else if (strcmp(argv[i], "--all") == 0) { - run4 = run16 = run17 = runRealOnly = runDynamic = true; + run4 = run16 = run17 = runRealOnly = runDynamic = runMultiEigen = true; } } @@ -1292,6 +1468,7 @@ int main(int argc, const char *argv[]) { if (run17) result |= runTest17(gpuDev); if (runRealOnly) result |= runTestRealOnly(gpuDev); if (runDynamic) result |= runTestDynamicTransition(gpuDev); + if (runMultiEigen) result |= runTestMultiEigen(gpuDev); return result; } diff --git a/libhmsbeagle/GPU/BeagleGPUSpectralImpl.h b/libhmsbeagle/GPU/BeagleGPUSpectralImpl.h index b33e25a5..a52a1504 100644 --- a/libhmsbeagle/GPU/BeagleGPUSpectralImpl.h +++ b/libhmsbeagle/GPU/BeagleGPUSpectralImpl.h @@ -92,6 +92,15 @@ class BeagleGPUSpectralImpl : public BeagleGPUImpl { GPUPtr* dGradient; bool* hEigenDecompIsAllReal; + /* Per-eigendecomposition element/byte stride shared by dEvecTOrigin, + * dIevcTOrigin, and dGradientOrigin, i.e. AlignMemOffset(S*S*sizeof(Real)) + * (as elements / as bytes). Device alignment padding can make this + * larger than S*S, so any code that indexes into these pooled origin + * buffers by eigendecomposition index (e.g. the adjoint offset-queue's + * rec[3]/rec[7], or the gradient MemsetZero) must use this, not S*S. */ + unsigned int kSpectralMatrixStrideElements; + size_t kSpectralMatrixStrideBytes; + /* Offset-queue for the batched adjoint kernel: one 9-unsigned-int record * per branch (see STATUS.md for the field layout), addressed via * integer offsets into the pooled partials/states/evec/ievc/eigenvalue/ diff --git a/libhmsbeagle/GPU/BeagleGPUSpectralImpl.hpp b/libhmsbeagle/GPU/BeagleGPUSpectralImpl.hpp index dfc2c94c..f529b985 100644 --- a/libhmsbeagle/GPU/BeagleGPUSpectralImpl.hpp +++ b/libhmsbeagle/GPU/BeagleGPUSpectralImpl.hpp @@ -123,6 +123,13 @@ int BeagleGPUSpectralImpl::createInstance( // Backward eigenvector arrays: one S*S block per eigen decomposition. int S = this->kPaddedStateCount; size_t matStride = gpuIf->AlignMemOffset((size_t)S * S * sizeof(Real)); + kSpectralMatrixStrideBytes = matStride; + kSpectralMatrixStrideElements = (unsigned int)(matStride / sizeof(Real)); + if (getenv("BEAGLE_DEBUG_STRIDE")) { + fprintf(stderr, "[BEAGLE_DEBUG_STRIDE] S=%d matStride=%zu bytes, S*S*sizeof(Real)=%zu bytes, " + "aligned>logical=%s\n", S, matStride, (size_t)S * S * sizeof(Real), + (matStride > (size_t)S * S * sizeof(Real)) ? "YES" : "no"); + } dEvecT = (GPUPtr*) malloc(sizeof(GPUPtr) * eigenDecompositionCount); dIevcT = (GPUPtr*) malloc(sizeof(GPUPtr) * eigenDecompositionCount); dEvecTOrigin = gpuIf->AllocateMemory(eigenDecompositionCount * matStride); @@ -642,7 +649,7 @@ int BeagleGPUSpectralImpl::calculateAdjointCrossProducts( * through the parent view (intermittent large garbage gradients, * OpenCL, Apple M1 Max) even on the in-order command queue — filling * through the same view the kernel uses avoids that hazard. */ - gpuIf->MemsetZero(dGradientOrigin, (size_t)kSpectralEigenDecompCount * SS * sizeof(Real)); + gpuIf->MemsetZero(dGradientOrigin, (size_t)kSpectralEigenDecompCount * kSpectralMatrixStrideBytes); GPUPtr catW = this->dWeights[cwIdx]; @@ -680,11 +687,11 @@ int BeagleGPUSpectralImpl::calculateAdjointCrossProducts( rec[1] = isStates ? this->getStatesOffsetElements(postIdx) : this->getPartialsOffsetElements(postIdx); rec[2] = isStates ? 1u : 0u; - rec[3] = (unsigned int)ei * SS; + rec[3] = (unsigned int)ei * kSpectralMatrixStrideElements; rec[4] = (unsigned int)ei * this->getEvecStrideElements(); rec[5] = (unsigned int)ei * this->getEvalStrideElements(); rec[6] = (unsigned int)matIdx * kSpectralDistanceStrideElements; - rec[7] = (unsigned int)ei * SS; + rec[7] = (unsigned int)ei * kSpectralMatrixStrideElements; rec[8] = isAllReal ? 1u : 0u; } From b33c7ac76822047de4a08da55a00b8d2cdd8c6b8 Mon Sep 17 00:00:00 2001 From: fil-monti Date: Sun, 16 Aug 2026 16:10:00 +0200 Subject: [PATCH 13/13] GPU spectral: tune adjoint kernel block size to match AMD wavefront Retunes kernelAdjointMergedN launch block width from 128 to ADJOINT_BLOCK_SP_N=64, matching this AMD hardware native wavefront and avoiding inactive-thread overhead in the per-block pattern reduction at moderate pattern counts. Cuts the kernel own device time by about 35 percent and overall wall time by about 18 percent on the isolated spectral exact-gradient benchmark, the largest real win found in this line of investigation. AdjointBlockSize.h is the single shared source of truth for the block-size constant, consumed by both the host launcher and the generated OpenCL kernel source. Also adds BEAGLE_DEBUG_KERNEL_ARGS reporting of requested and max work-group size, preferred multiple, and local and private memory usage, and per-branch composition counters (states vs partials, real vs complex eigendecomposition) used to interpret the block-size sweep. --- libhmsbeagle/GPU/AdjointBlockSize.h | 7 + libhmsbeagle/GPU/BeagleGPUSpectralImpl.h | 2 + libhmsbeagle/GPU/BeagleGPUSpectralImpl.hpp | 102 ++++++++++++ .../GPU/CMake_OpenCLSpectral/CMakeLists.txt | 1 + libhmsbeagle/GPU/GPUInterfaceOpenCL.cpp | 150 ++++++++++++++++-- libhmsbeagle/GPU/KernelLauncher.cpp | 4 +- .../GPU/kernels/kernelsSpectralIfDef.cu | 4 +- .../kernels/make_opencl_spectral_kernels.sh | 2 + 8 files changed, 260 insertions(+), 12 deletions(-) create mode 100644 libhmsbeagle/GPU/AdjointBlockSize.h diff --git a/libhmsbeagle/GPU/AdjointBlockSize.h b/libhmsbeagle/GPU/AdjointBlockSize.h new file mode 100644 index 00000000..82b9637b --- /dev/null +++ b/libhmsbeagle/GPU/AdjointBlockSize.h @@ -0,0 +1,7 @@ +#ifndef BEAGLE_GPU_ADJOINT_BLOCK_SIZE_H +#define BEAGLE_GPU_ADJOINT_BLOCK_SIZE_H + +/* Shared by the host launcher and the generated OpenCL generic-N source. */ +#define ADJOINT_BLOCK_SP_N 64 + +#endif diff --git a/libhmsbeagle/GPU/BeagleGPUSpectralImpl.h b/libhmsbeagle/GPU/BeagleGPUSpectralImpl.h index a52a1504..57e4f5e5 100644 --- a/libhmsbeagle/GPU/BeagleGPUSpectralImpl.h +++ b/libhmsbeagle/GPU/BeagleGPUSpectralImpl.h @@ -91,6 +91,8 @@ class BeagleGPUSpectralImpl : public BeagleGPUImpl { GPUPtr dGradientOrigin; GPUPtr* dGradient; bool* hEigenDecompIsAllReal; + int* hEigenDecompComplexPairLeaders; + int* hEigenDecompSingletonRows; /* Per-eigendecomposition element/byte stride shared by dEvecTOrigin, * dIevcTOrigin, and dGradientOrigin, i.e. AlignMemOffset(S*S*sizeof(Real)) diff --git a/libhmsbeagle/GPU/BeagleGPUSpectralImpl.hpp b/libhmsbeagle/GPU/BeagleGPUSpectralImpl.hpp index f529b985..e95b8533 100644 --- a/libhmsbeagle/GPU/BeagleGPUSpectralImpl.hpp +++ b/libhmsbeagle/GPU/BeagleGPUSpectralImpl.hpp @@ -18,6 +18,59 @@ #define __BeagleGPUSpectralImpl_hpp__ #include "libhmsbeagle/beagle.h" +#include +#include +#include + +namespace { +struct BeagleAdjointDetailStats { + long calls = 0; + double total = 0.0, integrate = 0.0, memset = 0.0, queueBuild = 0.0; + double queueH2D = 0.0, mergedLaunch = 0.0, gradientD2H = 0.0, hostCopy = 0.0; + ~BeagleAdjointDetailStats() { + if (!calls) return; + fprintf(stderr, "\n===== BEAGLE_PROFILE_ADJOINT_DETAIL =====\n"); +#define BEAGLE_ADJOINT_LINE(label, value) fprintf(stderr, "%-28s: calls=%-8ld total=%.6fs avg=%.3fus\n", label, calls, value, value * 1e6 / calls) + BEAGLE_ADJOINT_LINE("calculateAdjoint total", total); + BEAGLE_ADJOINT_LINE("IntegrateLikelihoods enqueue", integrate); + BEAGLE_ADJOINT_LINE("MemsetZero host wall", memset); + BEAGLE_ADJOINT_LINE("hAdjointQueue construction", queueBuild); + BEAGLE_ADJOINT_LINE("adjoint queue H2D", queueH2D); + BEAGLE_ADJOINT_LINE("merged adjoint enqueue", mergedLaunch); + BEAGLE_ADJOINT_LINE("gradient D2H", gradientD2H); + BEAGLE_ADJOINT_LINE("padded gradient host copy", hostCopy); +#undef BEAGLE_ADJOINT_LINE + fprintf(stderr, "===== END BEAGLE_PROFILE_ADJOINT_DETAIL =====\n\n"); + } +}; +struct BeagleAdjointCompositionStats { + long calls = 0, branches = 0, statesBranches = 0, allRealBranches = 0; + long complexPairLeaders = 0, singletonRows = 0; + int paddedPatterns = 0, paddedStates = 0, categories = 0; + ~BeagleAdjointCompositionStats() { + if (!calls) return; + fprintf(stderr, "\n===== BEAGLE_PROFILE_ADJOINT_COMPOSITION =====\n"); + fprintf(stderr, "calls=%ld branches=%ld isStates=%ld isAllReal=%ld " + "complexPairLeaders=%ld singletonRows=%ld " + "kPaddedPatternCount=%d kPaddedStateCount=%d kCategoryCount=%d\n", + calls, branches, statesBranches, allRealBranches, + complexPairLeaders, singletonRows, + paddedPatterns, paddedStates, categories); + fprintf(stderr, "===== END BEAGLE_PROFILE_ADJOINT_COMPOSITION =====\n\n"); + } +}; +inline BeagleAdjointCompositionStats& beagleAdjointCompositionStats() { + static BeagleAdjointCompositionStats stats; + return stats; +} +inline BeagleAdjointDetailStats& beagleAdjointDetailStats() { + static BeagleAdjointDetailStats stats; + return stats; +} +inline double beagleAdjointElapsed(std::chrono::high_resolution_clock::time_point t) { + return std::chrono::duration(std::chrono::high_resolution_clock::now() - t).count(); +} +} namespace beagle { namespace gpu { @@ -39,6 +92,7 @@ BeagleGPUSpectralImpl::BeagleGPUSpectralImpl() hSpectralPtrQueue(NULL), hSpectralDistanceQueue(NULL), dEvecTOrigin(0), dEvecT(NULL), dIevcTOrigin(0), dIevcT(NULL), dGradientOrigin(0), dGradient(NULL), hEigenDecompIsAllReal(NULL), + hEigenDecompComplexPairLeaders(NULL), hEigenDecompSingletonRows(NULL), dAdjointQueue(0), hAdjointQueue(NULL), kAdjointQueueCapacity(0), kSpectralEigenDecompCount(0) { } @@ -67,6 +121,8 @@ BeagleGPUSpectralImpl::~BeagleGPUSpectralImpl() { free(hPendingDenseEigenIndex); free(hPendingDenseEdgeLength); free(hEigenDecompIsAllReal); + free(hEigenDecompComplexPairLeaders); + free(hEigenDecompSingletonRows); } BEAGLE_GPU_TEMPLATE @@ -94,6 +150,8 @@ int BeagleGPUSpectralImpl::createInstance( kSpectralEigenDecompCount = eigenDecompositionCount; hEigenIndexForMatrix = (int*) calloc(this->kMatrixCount, sizeof(int)); + hEigenDecompComplexPairLeaders = (int*) calloc(eigenDecompositionCount, sizeof(int)); + hEigenDecompSingletonRows = (int*) calloc(eigenDecompositionCount, sizeof(int)); hDenseMatrixValid = (bool*) calloc(this->kMatrixCount, sizeof(bool)); // false = stale/never built hPendingDenseEigenIndex = (int*) calloc(this->kMatrixCount, sizeof(int)); @@ -350,6 +408,11 @@ int BeagleGPUSpectralImpl::setEigenDecomposition( for (int i = 0; i < SC && allReal; i++) if (inEigenValues[SC + i] != 0.0) allReal = false; hEigenDecompIsAllReal[eigenIndex] = allReal; + int pairLeaders = 0; + for (int i = 0; i < SC; i++) + if (inEigenValues[SC + i] > 0.0) pairLeaders++; + hEigenDecompComplexPairLeaders[eigenIndex] = pairLeaders; + hEigenDecompSingletonRows[eigenIndex] = S - 2 * pairLeaders; if (this->kFlags & BEAGLE_FLAG_EIGEN_COMPLEX) { for (int i = 0; i < SC; i++) @@ -623,12 +686,18 @@ int BeagleGPUSpectralImpl::calculateAdjointCrossProducts( double* outSumDerivatives, double* outSumSquaredDerivatives) { + const bool detail = (getenv("BEAGLE_PROFILE_ADJOINT_DETAIL") != NULL); + const bool composition = (getenv("BEAGLE_PROFILE_ADJOINT_COMPOSITION") != NULL); + BeagleAdjointDetailStats& detailStats = beagleAdjointDetailStats(); + std::chrono::high_resolution_clock::time_point stageStart, callStart; + if (detail) callStart = std::chrono::high_resolution_clock::now(); GPUInterface* gpuIf = this->gpu; const int S = this->kPaddedStateCount; const int SS = S * S; const int cwIdx = categoryWeightsIndices[0]; /* Compute per-site marginal likelihoods into dIntegrationTmp. */ + if (detail) stageStart = std::chrono::high_resolution_clock::now(); this->kernels->IntegrateLikelihoods( this->dIntegrationTmp, this->dPartials[rootPostOrderIndex], @@ -636,6 +705,7 @@ int BeagleGPUSpectralImpl::calculateAdjointCrossProducts( this->dFrequencies[stateFrequenciesIndex], this->kPaddedPatternCount, this->kCategoryCount); + if (detail) detailStats.integrate += beagleAdjointElapsed(stageStart); /* Zero the gradient accumulator on-device, in one fill over the whole * pooled origin buffer (covers every eigen decomposition, not just @@ -649,11 +719,14 @@ int BeagleGPUSpectralImpl::calculateAdjointCrossProducts( * through the parent view (intermittent large garbage gradients, * OpenCL, Apple M1 Max) even on the in-order command queue — filling * through the same view the kernel uses avoids that hazard. */ + if (detail) stageStart = std::chrono::high_resolution_clock::now(); gpuIf->MemsetZero(dGradientOrigin, (size_t)kSpectralEigenDecompCount * kSpectralMatrixStrideBytes); + if (detail) detailStats.memset += beagleAdjointElapsed(stageStart); GPUPtr catW = this->dWeights[cwIdx]; if (count > 0) { + if (detail) stageStart = std::chrono::high_resolution_clock::now(); /* Merged single-launch path (both S=4 and generic-N): build one * offset-queue record per branch (single buffer, integer offsets * into pooled origins — never per-branch device pointers), then one @@ -681,6 +754,14 @@ int BeagleGPUSpectralImpl::calculateAdjointCrossProducts( const int ei = hEigenIndexForMatrix[matIdx]; const bool isStates = (this->dStates[postIdx] != 0); const bool isAllReal = hEigenDecompIsAllReal[ei]; + if (composition) { + BeagleAdjointCompositionStats& s = beagleAdjointCompositionStats(); + s.branches++; + if (isStates) s.statesBranches++; + if (isAllReal) s.allRealBranches++; + s.complexPairLeaders += hEigenDecompComplexPairLeaders[ei]; + s.singletonRows += hEigenDecompSingletonRows[ei]; + } unsigned int* rec = hAdjointQueue + (size_t)i * kAdjointQueueFieldsPerBranch; rec[0] = this->getPartialsOffsetElements(preIdx); @@ -694,10 +775,14 @@ int BeagleGPUSpectralImpl::calculateAdjointCrossProducts( rec[7] = (unsigned int)ei * kSpectralMatrixStrideElements; rec[8] = isAllReal ? 1u : 0u; } + if (detail) detailStats.queueBuild += beagleAdjointElapsed(stageStart); + if (detail) stageStart = std::chrono::high_resolution_clock::now(); gpuIf->MemcpyHostToDevice(dAdjointQueue, hAdjointQueue, sizeof(unsigned int) * count * kAdjointQueueFieldsPerBranch); + if (detail) detailStats.queueH2D += beagleAdjointElapsed(stageStart); + if (detail) stageStart = std::chrono::high_resolution_clock::now(); if (S == 4) { this->kernels->AdjointCrossProductMergedSpectral4( this->getPartialsOrigin(), this->getStatesOrigin(), dEvecTOrigin, @@ -713,6 +798,7 @@ int BeagleGPUSpectralImpl::calculateAdjointCrossProducts( this->dIntegrationTmp, dGradientOrigin, dAdjointQueue, this->kPaddedPatternCount, this->kCategoryCount, count); } + if (detail) detailStats.mergedLaunch += beagleAdjointElapsed(stageStart); } /* Download gradient to host. The first eigen decomp's gradient is the @@ -730,13 +816,29 @@ int BeagleGPUSpectralImpl::calculateAdjointCrossProducts( const int ei0 = hEigenIndexForMatrix[eigenIndices[0]]; const int SC = this->kStateCount; Real* hGrad = (Real*) gpuIf->CallocHost(sizeof(Real), SS); + if (detail) stageStart = std::chrono::high_resolution_clock::now(); gpuIf->MemcpyDeviceToHost(hGrad, dGradient[ei0], SS * sizeof(Real)); + if (detail) detailStats.gradientD2H += beagleAdjointElapsed(stageStart); + if (detail) stageStart = std::chrono::high_resolution_clock::now(); for (int i = 0; i < SC; i++) for (int j = 0; j < SC; j++) outSumDerivatives[i * SC + j] = (double)hGrad[i * S + j]; + if (detail) detailStats.hostCopy += beagleAdjointElapsed(stageStart); gpuIf->FreeHostMemory(hGrad); } + if (detail) { + detailStats.calls++; + detailStats.total += beagleAdjointElapsed(callStart); + } + if (composition) { + BeagleAdjointCompositionStats& s = beagleAdjointCompositionStats(); + s.calls++; + s.paddedPatterns = this->kPaddedPatternCount; + s.paddedStates = this->kPaddedStateCount; + s.categories = this->kCategoryCount; + } + return BEAGLE_SUCCESS; } diff --git a/libhmsbeagle/GPU/CMake_OpenCLSpectral/CMakeLists.txt b/libhmsbeagle/GPU/CMake_OpenCLSpectral/CMakeLists.txt index 010930fc..57f8b5ca 100644 --- a/libhmsbeagle/GPU/CMake_OpenCLSpectral/CMakeLists.txt +++ b/libhmsbeagle/GPU/CMake_OpenCLSpectral/CMakeLists.txt @@ -22,6 +22,7 @@ add_custom_command( "../kernels/kernelsXDerivatives.cu" "../kernels/kernelsAll.cu" "../kernels/kernelsSpectralIfDef.cu" + "../AdjointBlockSize.h" "../kernels/kernelsSpectralIfDef4.cu" COMMENT "Building OpenCL spectral kernels" VERBATIM) diff --git a/libhmsbeagle/GPU/GPUInterfaceOpenCL.cpp b/libhmsbeagle/GPU/GPUInterfaceOpenCL.cpp index eae6124d..b5be7667 100644 --- a/libhmsbeagle/GPU/GPUInterfaceOpenCL.cpp +++ b/libhmsbeagle/GPU/GPUInterfaceOpenCL.cpp @@ -22,6 +22,8 @@ #include #include #include +#include +#include #include #include "libhmsbeagle/beagle.h" @@ -68,7 +70,110 @@ struct BeagleProfileStats { size_t memsetZeroBytes = 0; }; +struct BeagleOpenCLEventStat { + long calls = 0; + double seconds = 0.0; +}; + +struct BeaglePendingEvent { + cl_event event; + std::string name; +}; + BeagleProfileStats gBeagleProfile; +bool gBeagleOpenCLEventsInitialized = false; +bool gBeagleOpenCLEventsEnabled = false; +std::vector gBeaglePendingEvents; +std::map gBeagleOpenCLEventStats; + +void beagleOpenCLEventPrintSummary() { + if (!gBeagleOpenCLEventsEnabled) return; + fprintf(stderr, "\n===== BEAGLE_PROFILE_OPENCL_EVENTS =====\n"); + for (std::map::const_iterator it = + gBeagleOpenCLEventStats.begin(); it != gBeagleOpenCLEventStats.end(); ++it) { + fprintf(stderr, "%-28s: calls=%-8ld device=%.6fs avg=%.3fus\n", + it->first.c_str(), it->second.calls, it->second.seconds, + it->second.calls ? it->second.seconds * 1e6 / it->second.calls : 0.0); + } + fprintf(stderr, "pending events not yet drained: %zu\n", gBeaglePendingEvents.size()); + fprintf(stderr, "===== END BEAGLE_PROFILE_OPENCL_EVENTS =====\n\n"); +} + +inline bool beagleOpenCLEventsEnabled() { + if (!gBeagleOpenCLEventsInitialized) { + gBeagleOpenCLEventsEnabled = (getenv("BEAGLE_PROFILE_OPENCL_EVENTS") != NULL); + gBeagleOpenCLEventsInitialized = true; + if (gBeagleOpenCLEventsEnabled) atexit(beagleOpenCLEventPrintSummary); + } + return gBeagleOpenCLEventsEnabled; +} + +void beagleRecordCompletedEvent(const char* name, cl_event event) { + cl_ulong start = 0, end = 0; + int error = clGetEventProfilingInfo(event, CL_PROFILING_COMMAND_START, sizeof(start), &start, NULL); + if (error != CL_SUCCESS) { fprintf(stderr, "OpenCL event profiling start query failed: %d\n", error); exit(-1); } + error = clGetEventProfilingInfo(event, CL_PROFILING_COMMAND_END, sizeof(end), &end, NULL); + if (error != CL_SUCCESS) { fprintf(stderr, "OpenCL event profiling end query failed: %d\n", error); exit(-1); } + BeagleOpenCLEventStat& s = gBeagleOpenCLEventStats[name]; + s.calls++; + s.seconds += (double)(end - start) * 1e-9; + error = clReleaseEvent(event); + if (error != CL_SUCCESS) { fprintf(stderr, "OpenCL event release failed: %d\n", error); exit(-1); } +} + +void beagleDrainPendingEvents() { + for (size_t i = 0; i < gBeaglePendingEvents.size(); ++i) + beagleRecordCompletedEvent(gBeaglePendingEvents[i].name.c_str(), + gBeaglePendingEvents[i].event); + gBeaglePendingEvents.clear(); +} + +bool beagleTargetKernel(GPUFunction kernel, char* name, size_t size) { + if (!beagleOpenCLEventsEnabled()) return false; + int error = clGetKernelInfo(kernel, CL_KERNEL_FUNCTION_NAME, size, name, NULL); + if (error != CL_SUCCESS) { fprintf(stderr, "OpenCL kernel-name query failed: %d\n", error); exit(-1); } + return strcmp(name, "kernelAdjointMergedN") == 0 || + strcmp(name, "kernelAdjointMerged4") == 0 || + strcmp(name, "kernelIntegrateLikelihoods") == 0; +} + +void beagleReportAdjointKernelInfoOnce(GPUFunction kernel, cl_device_id device, + size_t requestedLocalSize) { + static bool reported = false; + if (reported || getenv("BEAGLE_PROFILE_ADJOINT_KERNEL_INFO") == NULL) return; + char name[256] = {0}; + int error = clGetKernelInfo(kernel, CL_KERNEL_FUNCTION_NAME, sizeof(name), name, NULL); + if (error != CL_SUCCESS) { fprintf(stderr, "OpenCL adjoint kernel-name query failed: %d\n", error); exit(-1); } + if (strcmp(name, "kernelAdjointMergedN") != 0) return; + size_t maxWorkGroup = 0, preferredMultiple = 0; + cl_ulong localMemory = 0, privateMemory = 0; + error = clGetKernelWorkGroupInfo(kernel, device, CL_KERNEL_WORK_GROUP_SIZE, + sizeof(maxWorkGroup), &maxWorkGroup, NULL); + if (error != CL_SUCCESS) { fprintf(stderr, "OpenCL work-group-size query failed: %d\n", error); exit(-1); } + error = clGetKernelWorkGroupInfo(kernel, device, + CL_KERNEL_PREFERRED_WORK_GROUP_SIZE_MULTIPLE, + sizeof(preferredMultiple), &preferredMultiple, NULL); + if (error != CL_SUCCESS) { fprintf(stderr, "OpenCL preferred-multiple query failed: %d\n", error); exit(-1); } + error = clGetKernelWorkGroupInfo(kernel, device, CL_KERNEL_LOCAL_MEM_SIZE, + sizeof(localMemory), &localMemory, NULL); + if (error != CL_SUCCESS) { fprintf(stderr, "OpenCL local-memory query failed: %d\n", error); exit(-1); } + error = clGetKernelWorkGroupInfo(kernel, device, CL_KERNEL_PRIVATE_MEM_SIZE, + sizeof(privateMemory), &privateMemory, NULL); + if (error != CL_SUCCESS) { fprintf(stderr, "OpenCL private-memory query failed: %d\n", error); exit(-1); } + fprintf(stderr, "\n===== BEAGLE_PROFILE_ADJOINT_KERNEL_INFO =====\n"); + fprintf(stderr, "kernel=%s requestedLocalSize=%zu workGroupSize=%zu " + "preferredMultiple=%zu localMemory=%llu privateMemory=%llu " + "valid=%s\n", + name, requestedLocalSize, maxWorkGroup, preferredMultiple, + (unsigned long long)localMemory, (unsigned long long)privateMemory, + requestedLocalSize <= maxWorkGroup ? "YES" : "NO"); + fprintf(stderr, "===== END BEAGLE_PROFILE_ADJOINT_KERNEL_INFO =====\n\n"); + if (requestedLocalSize > maxWorkGroup) { + fprintf(stderr, "kernelAdjointMergedN local size exceeds CL_KERNEL_WORK_GROUP_SIZE\n"); + exit(-1); + } + reported = true; +} void beagleProfilePrintSummary() { if (!gBeagleProfile.enabled) return; @@ -77,7 +182,8 @@ void beagleProfilePrintSummary() { p.syncHostSeconds + p.memcpyH2DSeconds + p.memcpyD2HSeconds + p.memsetZeroSeconds; fprintf(stderr, "\n===== BEAGLE_PROFILE_SUMMARY (GPUInterfaceOpenCL) =====\n"); - fprintf(stderr, "warmup discarded: first %ld SynchronizeHost calls\n", BEAGLE_PROFILE_WARMUP_SYNC); + fprintf(stderr, "warmup discarded: %s\n", getenv("BEAGLE_PROFILE_NO_WARMUP") + ? "disabled by BEAGLE_PROFILE_NO_WARMUP" : "first 300 SynchronizeHost calls"); fprintf(stderr, "LaunchKernel : calls=%-8ld total=%.6fs avg=%.3fus\n", p.launchKernelCalls, p.launchKernelSeconds, p.launchKernelCalls ? p.launchKernelSeconds * 1e6 / p.launchKernelCalls : 0.0); @@ -107,6 +213,10 @@ void beagleProfilePrintSummary() { inline bool beagleProfilingEnabled() { if (!gBeagleProfile.registered) { gBeagleProfile.enabled = (getenv("BEAGLE_PROFILE_SUMMARY") != NULL); + if (gBeagleProfile.enabled && getenv("BEAGLE_PROFILE_NO_WARMUP") != NULL) { + gBeagleProfile.warm = true; + gBeagleProfile.warmupRemaining = 0; + } gBeagleProfile.registered = true; if (gBeagleProfile.enabled) { atexit(beagleProfilePrintSummary); @@ -413,7 +523,8 @@ void GPUInterface::SetDevice(int deviceNumber, openClCommandQueues = (cl_command_queue*) malloc(sizeof(cl_command_queue) * BEAGLE_STREAM_COUNT); openClEvents = (cl_event*) malloc(sizeof(cl_event) * BEAGLE_STREAM_COUNT); - cl_command_queue_properties queueProperties = 0;//CL_QUEUE_OUT_OF_ORDER_EXEC_MODE_ENABLE; + cl_command_queue_properties queueProperties = beagleOpenCLEventsEnabled() + ? CL_QUEUE_PROFILING_ENABLE : 0;//CL_QUEUE_OUT_OF_ORDER_EXEC_MODE_ENABLE; for (int i=0; i < BEAGLE_STREAM_COUNT; i++) { openClCommandQueues[i] = clCreateCommandQueue(openClContext, openClDeviceId, queueProperties, &err); @@ -730,16 +841,23 @@ void GPUInterface::LaunchKernel(GPUFunction deviceFunction, fflush(stderr); } + beagleReportAdjointKernelInfoOnce(deviceFunction, openClDeviceId, localWorkSize[0]); + char eventKernelName[256] = {0}; + const bool targetEvent = beagleTargetKernel(deviceFunction, eventKernelName, + sizeof(eventKernelName)); + cl_event kernelEvent = NULL; + cl_event* kernelEventPtr = targetEvent ? &kernelEvent : NULL; if (globalWorkSize[1] == 1 && globalWorkSize[2] == 1) { SAFE_CL(clEnqueueNDRangeKernel(openClCommandQueues[0], deviceFunction, 1, NULL, - globalWorkSize, localWorkSize, 0, NULL, NULL)); + globalWorkSize, localWorkSize, 0, NULL, kernelEventPtr)); } else if (globalWorkSize[2] == 1) { SAFE_CL(clEnqueueNDRangeKernel(openClCommandQueues[0], deviceFunction, 2, NULL, - globalWorkSize, localWorkSize, 0, NULL, NULL)); + globalWorkSize, localWorkSize, 0, NULL, kernelEventPtr)); } else { SAFE_CL(clEnqueueNDRangeKernel(openClCommandQueues[0], deviceFunction, 3, NULL, - globalWorkSize, localWorkSize, 0, NULL, NULL)); + globalWorkSize, localWorkSize, 0, NULL, kernelEventPtr)); } + if (targetEvent) gBeaglePendingEvents.push_back(BeaglePendingEvent{kernelEvent, eventKernelName}); if (prof && gBeagleProfile.warm) { double dt = std::chrono::duration( @@ -1153,8 +1271,13 @@ void GPUInterface::MemsetZero(GPUPtr dest, * clEnqueueWriteBuffer path MemcpyHostToDevice already uses elsewhere * in this file, which has no such issue. */ void* hZero = calloc(1, memSize); + cl_event event = NULL; SAFE_CL(clEnqueueWriteBuffer(openClCommandQueues[0], dest, CL_TRUE, 0, memSize, hZero, 0, - NULL, NULL)); + NULL, beagleOpenCLEventsEnabled() ? &event : NULL)); + if (event) { + beagleDrainPendingEvents(); + beagleRecordCompletedEvent("MemsetZero", event); + } free(hZero); if (prof && gBeagleProfile.warm) { @@ -1184,8 +1307,13 @@ void GPUInterface::MemcpyHostToDevice(GPUPtr dest, // CL_TRUE => blocking write: this call does not return until the transfer // completes, regardless of payload size. + cl_event event = NULL; SAFE_CL(clEnqueueWriteBuffer(openClCommandQueues[0], dest, CL_TRUE, 0, memSize, src, 0, - NULL, NULL)); + NULL, beagleOpenCLEventsEnabled() ? &event : NULL)); + if (event) { + beagleDrainPendingEvents(); + beagleRecordCompletedEvent("MemcpyHostToDevice", event); + } if (prof && gBeagleProfile.warm) { double dt = std::chrono::duration( @@ -1213,8 +1341,13 @@ void GPUInterface::MemcpyDeviceToHost(void* dest, // CL_TRUE => blocking read: this call does not return until the transfer // completes, regardless of payload size. + cl_event event = NULL; SAFE_CL(clEnqueueReadBuffer(openClCommandQueues[0], src, CL_TRUE, 0, memSize, dest, 0, - NULL, NULL)); + NULL, beagleOpenCLEventsEnabled() ? &event : NULL)); + if (event) { + beagleDrainPendingEvents(); + beagleRecordCompletedEvent("MemcpyDeviceToHost", event); + } if (prof && gBeagleProfile.warm) { double dt = std::chrono::duration( @@ -1612,4 +1745,3 @@ void GPUInterface::ReleaseDevice(int deviceNumber) { }; // namespace - diff --git a/libhmsbeagle/GPU/KernelLauncher.cpp b/libhmsbeagle/GPU/KernelLauncher.cpp index 6adba33c..73463387 100644 --- a/libhmsbeagle/GPU/KernelLauncher.cpp +++ b/libhmsbeagle/GPU/KernelLauncher.cpp @@ -24,6 +24,7 @@ #include "libhmsbeagle/beagle.h" #include "libhmsbeagle/GPU/GPUImplDefs.h" #include "libhmsbeagle/GPU/KernelLauncher.h" +#include "libhmsbeagle/GPU/AdjointBlockSize.h" /**************CODE***********/ @@ -1176,7 +1177,7 @@ void KernelLauncher::AdjointCrossProductMergedSpectralN(GPUPtr partialsOrigin, unsigned int patternCount, unsigned int categoryCount, int branchCount) { - Dim3Int block(128); + Dim3Int block(ADJOINT_BLOCK_SP_N); Dim3Int grid(kPaddedStateCount, categoryCount, branchCount); gpu->LaunchKernel(fAdjointMergedN, block, grid, 11, 12, partialsOrigin, statesOrigin, evecTOrigin, ievcOrigin, evalOrigin, @@ -2760,4 +2761,3 @@ void KernelLauncher::SumSites3(GPUPtr dArray1, }; // namespace - diff --git a/libhmsbeagle/GPU/kernels/kernelsSpectralIfDef.cu b/libhmsbeagle/GPU/kernels/kernelsSpectralIfDef.cu index 038d1291..82885eec 100644 --- a/libhmsbeagle/GPU/kernels/kernelsSpectralIfDef.cu +++ b/libhmsbeagle/GPU/kernels/kernelsSpectralIfDef.cu @@ -893,7 +893,9 @@ KW_GLOBAL_KERNEL void kernelPartialsStatesGrowingTopRootSpectral( #pragma OPENCL EXTENSION cl_khr_int64_base_atomics : enable #endif -#define ADJOINT_BLOCK_SP_N 128 +#ifndef ADJOINT_BLOCK_SP_N +#error "ADJOINT_BLOCK_SP_N must be supplied by AdjointBlockSize.h" +#endif /* 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. diff --git a/libhmsbeagle/GPU/kernels/make_opencl_spectral_kernels.sh b/libhmsbeagle/GPU/kernels/make_opencl_spectral_kernels.sh index b5bc5ed1..19ef629a 100755 --- a/libhmsbeagle/GPU/kernels/make_opencl_spectral_kernels.sh +++ b/libhmsbeagle/GPU/kernels/make_opencl_spectral_kernels.sh @@ -37,6 +37,7 @@ for s in $STATE_COUNT_LIST; do cat $srcdir/kernelsAll.cu | sed 's/\\/\\\\/g' | sed 's/\"/\\"/g' | sed 's/$/\\n\\/' >> BeagleOpenCLSpectral_kernels.h cat $srcdir/kernelsX.cu | sed 's/\\/\\\\/g' | sed 's/\"/\\"/g' | sed 's/$/\\n\\/' >> BeagleOpenCLSpectral_kernels.h cat $srcdir/kernelsXDerivatives.cu | sed 's/\\/\\\\/g' | sed 's/\"/\\"/g' | sed 's/$/\\n\\/' >> BeagleOpenCLSpectral_kernels.h + cat $srcdir/../AdjointBlockSize.h | sed 's/\\/\\\\/g' | sed 's/\"/\\"/g' | sed 's/$/\\n\\/' >> BeagleOpenCLSpectral_kernels.h cat $srcdir/kernelsSpectralIfDef.cu | sed 's/\\/\\\\/g' | sed 's/\"/\\"/g' | sed 's/$/\\n\\/' >> BeagleOpenCLSpectral_kernels.h echo "\"" >> BeagleOpenCLSpectral_kernels.h done @@ -66,6 +67,7 @@ for s in $STATE_COUNT_LIST; do cat $srcdir/kernelsAll.cu | sed 's/\\/\\\\/g' | sed 's/\"/\\"/g' | sed 's/$/\\n\\/' >> BeagleOpenCLSpectral_kernels.h cat $srcdir/kernelsX.cu | sed 's/\\/\\\\/g' | sed 's/\"/\\"/g' | sed 's/$/\\n\\/' >> BeagleOpenCLSpectral_kernels.h cat $srcdir/kernelsXDerivatives.cu | sed 's/\\/\\\\/g' | sed 's/\"/\\"/g' | sed 's/$/\\n\\/' >> BeagleOpenCLSpectral_kernels.h + cat $srcdir/../AdjointBlockSize.h | sed 's/\\/\\\\/g' | sed 's/\"/\\"/g' | sed 's/$/\\n\\/' >> BeagleOpenCLSpectral_kernels.h cat $srcdir/kernelsSpectralIfDef.cu | sed 's/\\/\\\\/g' | sed 's/\"/\\"/g' | sed 's/$/\\n\\/' >> BeagleOpenCLSpectral_kernels.h echo "\"" >> BeagleOpenCLSpectral_kernels.h done