Address PR #242 review comments (clGetKernelInfo gating, imaginary-eigenvalue read gating)

.gitignore

@@ -209,6 +209,7 @@ examples/tinytest/.libs/
 examples/tinytest/tinytest.o
 libhmsbeagle/GPU/kernels/BeagleOpenCL_kernels.h
 libhmsbeagle/GPU/kernels/BeagleOpenCL_kernels_xcode.h
+libhmsbeagle/GPU/kernels/BeagleOpenCLSpectral_kernels.h
 libhmsbeagle/GPU/libhmsbeagle-opencl.la
 libhmsbeagle/GPU/libhmsbeagle_opencl_la-GPUImplHelper.lo
 libhmsbeagle/GPU/libhmsbeagle_opencl_la-GPUInterfaceOpenCL.lo

libhmsbeagle/GPU/AdjointBlockSize.h

@@ -0,0 +1,7 @@
+#ifndef BEAGLE_GPU_ADJOINT_BLOCK_SIZE_H
+#define BEAGLE_GPU_ADJOINT_BLOCK_SIZE_H
+
+/* Shared by the host launcher and the generated OpenCL generic-N source. */
+#define ADJOINT_BLOCK_SP_N 64
+
+#endif

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;
@@ -4756,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,

libhmsbeagle/GPU/BeagleGPUSpectralImpl.h

@@ -38,13 +38,42 @@ class BeagleGPUSpectralImpl : public BeagleGPUImpl<Real> {
     GPUPtr  dSpectralDistancesOrigin;
     GPUPtr* dSpectralDistances;
     int*    hEigenIndexForMatrix;
+
+    /* Lazy dense-matrix materialization (see PLAN.md in
+     * beagle-bugs/gpu-spectral-skip-redundant-dense-matrix/): probability-only
+     * updateTransitionMatrices calls no longer eagerly build the dense
+     * kStateCount x kStateCount transition matrix, since the spectral
+     * pruning/gradient paths never read it. Instead the eigenIndex/edgeLength
+     * needed to rebuild it are cached here, and the dense matrix is
+     * materialized on first actual getTransitionMatrix() call. All sized
+     * kMatrixCount. */
+    bool*   hDenseMatrixValid;
+    int*    hPendingDenseEigenIndex;
+    double* hPendingDenseEdgeLength;
+
     /* Per-matrix element stride within dSpectralDistancesOrigin, i.e.
      * AlignMemOffset(kCategoryCount * sizeof(Real)) / sizeof(Real). Device
      * alignment padding can make this larger than kCategoryCount, so any code
      * that indexes into the pooled origin buffer by matrix index (e.g. the
      * adjoint offset-queue) must multiply by this, not by kCategoryCount. */
     unsigned int kSpectralDistanceStrideElements;
 
+    /* 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).
@@ -62,6 +91,17 @@ class BeagleGPUSpectralImpl : public BeagleGPUImpl<Real> {
     GPUPtr  dGradientOrigin;
     GPUPtr* dGradient;
     bool*   hEigenDecompIsAllReal;
+    int*    hEigenDecompComplexPairLeaders;
+    int*    hEigenDecompSingletonRows;
+
+    /* Per-eigendecomposition element/byte stride shared by dEvecTOrigin,
+     * dIevcTOrigin, and dGradientOrigin, i.e. AlignMemOffset(S*S*sizeof(Real))
+     * (as elements / as bytes). Device alignment padding can make this
+     * larger than S*S, so any code that indexes into these pooled origin
+     * buffers by eigendecomposition index (e.g. the adjoint offset-queue's
+     * rec[3]/rec[7], or the gradient MemsetZero) must use this, not S*S. */
+    unsigned int kSpectralMatrixStrideElements;
+    size_t       kSpectralMatrixStrideBytes;
 
     /* Offset-queue for the batched adjoint kernel: one 9-unsigned-int record
      * per branch (see STATUS.md for the field layout), addressed via
@@ -109,6 +149,9 @@ class BeagleGPUSpectralImpl : public BeagleGPUImpl<Real> {
                                  const double* edgeLengths,
                                  int count) override;
 
+    int getTransitionMatrix(int matrixIndex,
+                            double* outMatrix) override;
+
     int updatePrePartials(const int* operations,
                           int operationCount,
                           int cumulativeScaleIndex,