From f27c8aff35bc8eb324e9b024ed42684b6a676fda Mon Sep 17 00:00:00 2001
From: fil-monti <fmonti@ucla.edu>
Date: Fri, 10 Jul 2026 10:08:38 -0700
Subject: [PATCH 1/9] 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 <fmonti@ucla.edu>
Date: Fri, 10 Jul 2026 13:05:22 -0700
Subject: [PATCH 2/9] 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<Child1, States>) 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 <fmonti@ucla.edu>
Date: Fri, 10 Jul 2026 13:12:23 -0700
Subject: [PATCH 3/9] 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<BEAGLE_GPU_GENERIC>::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<BEAGLE_GPU_GENERIC>::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<BEAGLE_GPU_GENERIC>::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<BEAGLE_GPU_GENERIC>::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<BEAGLE_GPU_GENERIC>::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<BEAGLE_GPU_GENERIC>::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<BEAGLE_GPU_GENERIC>::dispatchGrowingSpectral(
                     ievc2, evec2, eval2, dist2,
                     nullptr, nullptr,
                     this->kPaddedPatternCount, this->kCategoryCount,
+                    isAllReal1, isAllReal2,
                     -1, -1, -1);
             }
         } else {
@@ -365,13 +374,15 @@ void BeagleGPUSpectralImpl<BEAGLE_GPU_GENERIC>::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 <fmonti@ucla.edu>
Date: Fri, 10 Jul 2026 13:12:31 -0700
Subject: [PATCH 4/9] 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<BEAGLE_GPU_GENERIC>::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<BEAGLE_GPU_GENERIC>::setEigenDecomposition(
     const int SC = this->kStateCount;
     const int SS = S * S;
     // kEigenValuesSize is private to BeagleGPUImpl; recompute the same way it does.
-    // 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 <fmonti@ucla.edu>
Date: Fri, 10 Jul 2026 13:13:16 -0700
Subject: [PATCH 5/9] 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<double>& evec,   /* N×N */
     std::vector<double>& ivec,   /* N×N */
-    std::vector<double>& eval)   /* 2N  */
+    std::vector<double>& 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<double> 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<double> 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<double, std::vector<double>> {
+        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<double> 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<double> 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<double> 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<double> evecReal, ivecReal, evalReal;
+    buildCirculantN(S, evecReal, ivecReal, evalReal, /*r_fwd=*/1.0, /*r_bkd=*/1.0);
+    std::vector<double> evecCplx, ivecCplx, evalCplx;
+    buildCirculantN(S, evecCplx, ivecCplx, evalCplx, /*r_fwd=*/1.0, /*r_bkd=*/0.5);
+
+    std::vector<double> 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<double>& evec,
+                           const std::vector<double>& ivec, const std::vector<double>& eval)
+                           -> std::pair<double, std::vector<double>> {
+        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<double> 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<double> 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<double> 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<double>* stepEvec[3] = { &evecReal, &evecCplx, &evecReal };
+    const std::vector<double>* stepIvec[3] = { &ivecReal, &ivecCplx, &ivecReal };
+    const std::vector<double>* stepEval[3] = { &evalReal, &evalCplx, &evalReal };
+
+    double cpuLogL[3];
+    std::vector<double> 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<double> 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 <fmonti@ucla.edu>
Date: Fri, 10 Jul 2026 22:58:09 -0700
Subject: [PATCH 6/9] 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 <cstdarg>
 #include <cmath>
 #include <map>
+#include <chrono>
 
 #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<double>(
+            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<double>(
+            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<double>(
+            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<double>(
+            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 <fmonti@ucla.edu>
Date: Fri, 10 Jul 2026 22:59:43 -0700
Subject: [PATCH 7/9] 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<BEAGLE_GPU_GENERIC>::dispatchPruneSS(GPUPtr s1, GPUPtr s2, GP
                                                          GPUPtr scalingFactors, GPUPtr cumulativeScaling,
                                                          unsigned int startPattern, unsigned int endPattern,
                                                          int rescale, int streamIndex, int waitIndex) {
-    fprintf(stderr, "[DISPATCH] BASE BeagleGPUImpl::dispatchPruneSS called\n"); fflush(stderr);
     kernels->StatesStatesPruningDynamicScaling(s1, s2, p3,
         dMatrices[c1MatIdx], dMatrices[c2MatIdx],
         scalingFactors, cumulativeScaling,
diff --git a/libhmsbeagle/GPU/BeagleGPUSpectralImpl.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<BEAGLE_GPU_GENERIC>::dispatchPruneSS(
         GPUPtr scalingFactors, GPUPtr cumulativeScaling,
         unsigned int startPattern, unsigned int endPattern,
         int rescale, int streamIndex, int waitIndex) {
-    fprintf(stderr, "[DISPATCH] SPECTRAL BeagleGPUSpectralImpl::dispatchPruneSS called\n"); fflush(stderr);
     int ei1 = hEigenIndexForMatrix[c1MatIdx];
     int ei2 = hEigenIndexForMatrix[c2MatIdx];
     this->kernels->StatesStatesPruningSpectral(s1, s2, p3,

From 224bd6cd00c488b42e3b93d30922794c7ea92611 Mon Sep 17 00:00:00 2001
From: fil-monti <fmonti@ucla.edu>
Date: Fri, 10 Jul 2026 23:01:08 -0700
Subject: [PATCH 8/9] 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<Real> {
      * 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<BEAGLE_GPU_GENERIC>::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<BEAGLE_GPU_GENERIC>::~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<BEAGLE_GPU_GENERIC>::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<BEAGLE_GPU_GENERIC>::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 <fmonti@ucla.edu>
Date: Sat, 11 Jul 2026 00:27:01 -0700
Subject: [PATCH 9/9] 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<double>(
+            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<double>(
+            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