diff --git a/examples/hmctest/adjointtest4.cpp b/examples/hmctest/adjointtest4.cpp index c28ee171..55532a6f 100644 --- a/examples/hmctest/adjointtest4.cpp +++ b/examples/hmctest/adjointtest4.cpp @@ -91,13 +91,18 @@ static double edgeLengths[4] = { 0.1, 0.1, 0.2, 0.5 }; * Evec[j, 2m-1] = (1/N)*cos(2πjm/N) Ievc[2m-1, j] = 2*cos(2πjm/N) * Evec[j, 2m ] = (1/N)*sin(2πjm/N) Ievc[2m, j] = 2*sin(2πjm/N) */ +/* r_fwd/r_bkd default to the asymmetric (complex-eigenvalue) case used by the + * existing 16-/17-state tests. Passing r_fwd == r_bkd instead makes Q + * symmetric/reversible: b = (r_fwd-r_bkd)*sin(theta_m) is then identically + * zero for every mode, so every eigenvalue comes out genuinely real (not + * just numerically close to zero) -- used by the real-only test below. */ static void buildCirculantN( int N, std::vector& evec, /* N×N */ std::vector& ivec, /* N×N */ - std::vector& eval) /* 2N */ + std::vector& eval, /* 2N */ + double r_fwd = 1.0, double r_bkd = 0.5) { - const double r_fwd = 1.0, r_bkd = 0.5; const double two_pi_over_N = 2.0 * M_PI / N; const bool evenN = (N % 2 == 0); const int nPairs = evenN ? N/2 - 1 : (N-1)/2; @@ -934,6 +939,318 @@ static int runTest17(int gpuDevice) { return gradOk ? 0 : 1; } +/* ═══════════════════════════════════════════════════════════════════════ + * Real-only 17-state test (BEAGLE_FLAG_EIGEN_COMPLEX NOT requested) + * + * Regression test for the msuchard PR #242 review comment ("if isAllReal + * then the imaginary parts should never be written or read"): the GPU + * spectral pruning kernels now gate their imaginary-eigenvalue reads on a + * per-decomposition isAllReal flag, and the device eigenvalue buffer is + * sized EIGEN_COMPLEX-only (matching CPU), not widened just because + * BEAGLE_FLAG_SPECTRAL_REPRESENTATION is in use. This test uses a + * *symmetric* circulant (r_fwd == r_bkd, so eigenvalues are exactly real, + * not just numerically close to zero) at a padded state count + * (17 -> kPaddedStateCount=32) and never requests + * BEAGLE_FLAG_EIGEN_COMPLEX, so the GPU device eigenvalue buffer is + * allocated at exactly kPaddedStateCount width -- no imaginary slack at + * all. If the kernel's imaginary-eigenvalue read were not properly gated, + * this configuration reads past the end of that buffer: the exact bug + * class PR #242 fixed and this follow-up refines. Run with + * BEAGLE_DEBUG_EIGEN=1 to see the device Eval buffer printed at its + * (narrow) width directly. + * ═══════════════════════════════════════════════════════════════════════*/ + +static int runTestRealOnly(int gpuDevice) { + const int S = 17, nPat = 4, nCats = 2; + const double tol = 1e-3; /* absolute tolerance for entries with |value|<=1 */ + const double relTol = 1e-2; /* relative tolerance for entries with |value|>1 (single vs double precision), same policy as runTest17 */ + + fprintf(stdout, "\n=== Real-only 17-state test (no BEAGLE_FLAG_EIGEN_COMPLEX) ===\n"); + fprintf(stdout, " r_fwd=r_bkd=1.0 -> symmetric/reversible circulant, all eigenvalues exactly real\n"); + fprintf(stdout, " Device eigenvalue buffer allocated at exactly kPaddedStateCount width (no imaginary slack)\n"); + + std::vector evec, ivec, eval; + buildCirculantN(S, evec, ivec, eval, /*r_fwd=*/1.0, /*r_bkd=*/1.0); + + fprintf(stdout, " Eigenvalue imag parts (must all be exactly 0): "); + bool anyImag = false; + for (int k = 0; k < S; k++) { + fprintf(stdout, "%.3f ", eval[S + k]); + if (eval[S + k] != 0.0) anyImag = true; + } + fprintf(stdout, "\n"); + if (anyImag) { + fprintf(stderr, " ERROR: expected all-zero imaginary eigenvalues for a symmetric circulant\n"); + return 1; + } + + std::vector freqs(S, 1.0 / S); + double rates_[2] = { 0.5, 1.5 }; + double catWts_[2] = { 0.5, 0.5 }; + double edgeLens_[4] = { 0.1, 0.1, 0.2, 0.5 }; + + int hSt[4] = { 3, 0, 5, 7 }; + int cSt[4] = { 3, 0, 5, 5 }; + int gSt[4] = { 0, 0, 0, 12 }; + + auto runOne = [&](bool useGpu, int dev, bool singlePrec) -> std::pair> { + int inst = createInstance(useGpu, dev, singlePrec, S, nPat, nCats, /*complexEigen=*/false); + if (inst < 0) return { 0.0, {} }; + + beagleSetTipStates(inst, 0, hSt); + beagleSetTipStates(inst, 1, cSt); + beagleSetTipStates(inst, 2, gSt); + + beagleSetCategoryRates(inst, rates_); + std::vector pw(nPat, 1.0); + beagleSetPatternWeights(inst, pw.data()); + beagleSetStateFrequencies(inst, 0, freqs.data()); + beagleSetCategoryWeights(inst, 0, catWts_); + beagleSetEigenDecomposition(inst, 0, evec.data(), ivec.data(), eval.data()); + + int nodeIdx[4] = { 0, 1, 2, 3 }; + beagleUpdateTransitionMatrices(inst, 0, nodeIdx, NULL, NULL, edgeLens_, 4); + + BeagleOperation postOps[2] = { + { 3, BEAGLE_OP_NONE, BEAGLE_OP_NONE, 0, 0, 1, 1 }, + { 4, BEAGLE_OP_NONE, BEAGLE_OP_NONE, 2, 2, 3, 3 } + }; + beagleUpdatePartials(inst, postOps, 2, BEAGLE_OP_NONE); + + double logL = 0.0; + int rootIdx = 4, cwIdx = 0, sfIdx = 0, csi = BEAGLE_OP_NONE; + beagleCalculateRootLogLikelihoods(inst, &rootIdx, &cwIdx, &sfIdx, &csi, 1, &logL); + + std::vector rootPre((size_t)nCats * nPat * S); + for (int c = 0; c < nCats; c++) + for (int p = 0; p < nPat; p++) + for (int s = 0; s < S; s++) + rootPre[(size_t)c * nPat * S + p * S + s] = freqs[s]; + beagleSetPartials(inst, 5, rootPre.data()); + + BeagleOperation preOps[4] = { + { 6, BEAGLE_OP_NONE, BEAGLE_OP_NONE, 5, BEAGLE_OP_NONE, 2, 2 }, + { 7, BEAGLE_OP_NONE, BEAGLE_OP_NONE, 5, BEAGLE_OP_NONE, 3, 3 }, + { 8, BEAGLE_OP_NONE, BEAGLE_OP_NONE, 6, BEAGLE_OP_NONE, 0, 0 }, + { 9, BEAGLE_OP_NONE, BEAGLE_OP_NONE, 6, BEAGLE_OP_NONE, 1, 1 } + }; + beagleUpdatePrePartials_v5(inst, preOps, 4, BEAGLE_OP_NONE, BEAGLE_PARTIALS_TOP); + + BeagleBranchOperation branchOps[4] = { + { 1, 8, 1, 0 }, + { 0, 9, 0, 0 }, + { 2, 7, 2, 0 }, + { 3, 6, 3, 0 } + }; + std::vector grad(S * S, 0.0); + beagleCalculateAdjointDerivative(inst, branchOps, 0, 0, 4, 0, 4, grad.data(), NULL); + + beagleFinalizeInstance(inst); + return { logL, grad }; + }; + + fprintf(stdout, "\nCreating CPU (double) instance:\n"); + auto [cpuLogL, cpuGrad] = runOne(false, -1, false); + if (cpuGrad.empty()) return 1; + fprintf(stdout, "CPU log-likelihood: %.6f\n", cpuLogL); + + if (gpuDevice < 0) return 0; + + fprintf(stdout, "\nCreating GPU (single) instance:\n"); + auto [gpuLogL, gpuGrad] = runOne(true, gpuDevice, true); + if (gpuGrad.empty()) return 1; + fprintf(stdout, "GPU log-likelihood: %.6f\n", gpuLogL); + + /* Gradient magnitudes here range up to ~190 (same character as runTest17's + * complex-eigenvalue case), so single-precision GPU vs double-precision + * CPU error scales with the *largest* magnitude in the computation -- + * threshold against the matrix-wide max magnitude instead of a flat + * absolute tolerance (same policy as runTest17). */ + double maxMagnitude = 0.0; + for (int k = 0; k < S*S; k++) { + if (!std::isnan(cpuGrad[k])) maxMagnitude = std::max(maxMagnitude, fabs(cpuGrad[k])); + if (!std::isnan(gpuGrad[k])) maxMagnitude = std::max(maxMagnitude, fabs(gpuGrad[k])); + } + const double threshold = std::max(tol, relTol * maxMagnitude); + + bool gradOk = true; + double maxDiff = 0.0; + for (int k = 0; k < S*S; k++) { + double d = fabs(cpuGrad[k] - gpuGrad[k]); + if (d > maxDiff || std::isnan(d)) maxDiff = d; + if (d > threshold || std::isnan(d)) gradOk = false; + } + fprintf(stdout, "Real-only 17-state gradient max|diff|=%.3e threshold=%.3e (rel to max|value|=%.3e) %s\n", + maxDiff, threshold, maxMagnitude, gradOk ? "PASS" : "FAIL"); + return gradOk ? 0 : 1; +} + +/* ═══════════════════════════════════════════════════════════════════════ + * Dynamic real <-> complex eigendecomposition transition test (16-state) + * + * A single BEAGLE_FLAG_EIGEN_COMPLEX-flagged instance (device eigenvalue + * buffer wide enough for either case) has beagleSetEigenDecomposition + * called three times in sequence for the same eigenIndex, alternating + * between a genuinely real-only decomposition and a genuinely complex one + * -- mirroring how one MCMC run can produce either kind of decomposition + * from one iteration to the next for a general asymmetric-rate model (see + * hEigenDecompIsAllReal's recomputation on every setEigenDecomposition + * call). Confirms the isAllReal gating doesn't leave stale sDs/sCs kernel + * state, or stale device Eval data, bleeding from one call into the next. + * ═══════════════════════════════════════════════════════════════════════*/ + +static int runTestDynamicTransition(int gpuDevice) { + const int S = 16, nPat = 4, nCats = 2; + const double tol = 1e-3; + + fprintf(stdout, "\n=== Dynamic real<->complex eigendecomposition transition test (16-state) ===\n"); + + std::vector evecReal, ivecReal, evalReal; + buildCirculantN(S, evecReal, ivecReal, evalReal, /*r_fwd=*/1.0, /*r_bkd=*/1.0); + std::vector evecCplx, ivecCplx, evalCplx; + buildCirculantN(S, evecCplx, ivecCplx, evalCplx, /*r_fwd=*/1.0, /*r_bkd=*/0.5); + + std::vector freqs(S, 1.0 / S); + double rates_[2] = { 0.5, 1.5 }; + double catWts_[2] = { 0.5, 0.5 }; + double edgeLens_[4] = { 0.1, 0.1, 0.2, 0.5 }; + int hSt[4] = { 3, 0, 5, 7 }; + int cSt[4] = { 3, 0, 5, 5 }; + int gSt[4] = { 0, 0, 0, 7 }; + + auto setupAndRun = [&](int inst, const std::vector& evec, + const std::vector& ivec, const std::vector& eval) + -> std::pair> { + beagleSetEigenDecomposition(inst, 0, evec.data(), ivec.data(), eval.data()); + int nodeIdx[4] = { 0, 1, 2, 3 }; + beagleUpdateTransitionMatrices(inst, 0, nodeIdx, NULL, NULL, edgeLens_, 4); + BeagleOperation postOps[2] = { + { 3, BEAGLE_OP_NONE, BEAGLE_OP_NONE, 0, 0, 1, 1 }, + { 4, BEAGLE_OP_NONE, BEAGLE_OP_NONE, 2, 2, 3, 3 } + }; + beagleUpdatePartials(inst, postOps, 2, BEAGLE_OP_NONE); + double logL = 0.0; + int rootIdx = 4, cwIdx = 0, sfIdx = 0, csi = BEAGLE_OP_NONE; + beagleCalculateRootLogLikelihoods(inst, &rootIdx, &cwIdx, &sfIdx, &csi, 1, &logL); + std::vector rootPre((size_t)nCats * nPat * S); + for (int c = 0; c < nCats; c++) + for (int p = 0; p < nPat; p++) + for (int s = 0; s < S; s++) + rootPre[(size_t)c * nPat * S + p * S + s] = freqs[s]; + beagleSetPartials(inst, 5, rootPre.data()); + BeagleOperation preOps[4] = { + { 6, BEAGLE_OP_NONE, BEAGLE_OP_NONE, 5, BEAGLE_OP_NONE, 2, 2 }, + { 7, BEAGLE_OP_NONE, BEAGLE_OP_NONE, 5, BEAGLE_OP_NONE, 3, 3 }, + { 8, BEAGLE_OP_NONE, BEAGLE_OP_NONE, 6, BEAGLE_OP_NONE, 0, 0 }, + { 9, BEAGLE_OP_NONE, BEAGLE_OP_NONE, 6, BEAGLE_OP_NONE, 1, 1 } + }; + beagleUpdatePrePartials_v5(inst, preOps, 4, BEAGLE_OP_NONE, BEAGLE_PARTIALS_TOP); + BeagleBranchOperation branchOps[4] = { + { 1, 8, 1, 0 }, + { 0, 9, 0, 0 }, + { 2, 7, 2, 0 }, + { 3, 6, 3, 0 } + }; + std::vector grad(S * S, 0.0); + beagleCalculateAdjointDerivative(inst, branchOps, 0, 0, 4, 0, 4, grad.data(), NULL); + return { logL, grad }; + }; + + auto makeInst = [&](bool useGpu, int dev, bool singlePrec) -> int { + int inst = createInstance(useGpu, dev, singlePrec, S, nPat, nCats, /*complexEigen=*/true); + if (inst < 0) return inst; + beagleSetTipStates(inst, 0, hSt); + beagleSetTipStates(inst, 1, cSt); + beagleSetTipStates(inst, 2, gSt); + beagleSetCategoryRates(inst, rates_); + std::vector pw(nPat, 1.0); + beagleSetPatternWeights(inst, pw.data()); + beagleSetStateFrequencies(inst, 0, freqs.data()); + beagleSetCategoryWeights(inst, 0, catWts_); + return inst; + }; + + fprintf(stdout, "\nCreating CPU (double) instance (BEAGLE_FLAG_EIGEN_COMPLEX requested):\n"); + int cpuInst = makeInst(false, -1, false); + if (cpuInst < 0) { fprintf(stderr, "Failed CPU instance\n"); return 1; } + + const char* stepNames[3] = { "real", "complex", "real (again)" }; + const std::vector* stepEvec[3] = { &evecReal, &evecCplx, &evecReal }; + const std::vector* stepIvec[3] = { &ivecReal, &ivecCplx, &ivecReal }; + const std::vector* stepEval[3] = { &evalReal, &evalCplx, &evalReal }; + + double cpuLogL[3]; + std::vector cpuGrad[3]; + for (int i = 0; i < 3; i++) { + auto r = setupAndRun(cpuInst, *stepEvec[i], *stepIvec[i], *stepEval[i]); + cpuLogL[i] = r.first; + cpuGrad[i] = r.second; + fprintf(stdout, " [CPU step %d: %-13s] logL=%.6f\n", i, stepNames[i], cpuLogL[i]); + } + beagleFinalizeInstance(cpuInst); + + if (gpuDevice < 0) return 0; + + fprintf(stdout, "\nCreating GPU (single) instance (BEAGLE_FLAG_EIGEN_COMPLEX requested):\n"); + int gpuInst = makeInst(true, gpuDevice, true); + if (gpuInst < 0) { fprintf(stderr, "Failed GPU instance\n"); return 1; } + + /* The "complex" step reuses runTest16's exact asymmetric circulant + * (r_fwd=1.0, r_bkd=0.5), which has a documented pre-existing marginal + * single-vs-double-precision diff of ~1.02e-03 (CLUSTER_AGENT_FINDINGS2.md + * section 8.2) unrelated to eigen-real/complex gating. Give that one step + * a matching looser tolerance so this test isn't tripped by that same + * known, already-accepted noise floor; the "real" steps (the actual + * subject of this test) are held to the tighter tol. */ + const double tolComplexStep = 1.1e-3; + std::vector gpuGradSteps[3]; + bool allOk = true; + for (int i = 0; i < 3; i++) { + auto r = setupAndRun(gpuInst, *stepEvec[i], *stepIvec[i], *stepEval[i]); + double gpuLogL = r.first; + gpuGradSteps[i] = r.second; + fprintf(stdout, " [GPU step %d: %-13s] logL=%.6f\n", i, stepNames[i], gpuLogL); + + double maxDiff = 0.0; + for (int k = 0; k < S*S; k++) { + double d = fabs(cpuGrad[i][k] - gpuGradSteps[i][k]); + if (d > maxDiff || std::isnan(d)) maxDiff = d; + } + double stepTol = (i == 1) ? tolComplexStep : tol; + bool stepOk = (maxDiff <= stepTol) && !std::isnan(maxDiff); + fprintf(stdout, " step %d (%s) gradient max|diff|=%.3e (tol=%.1e) %s\n", + i, stepNames[i], maxDiff, stepTol, stepOk ? "PASS" : "FAIL"); + allOk &= stepOk; + } + + /* Consistency check: step 0 and step 2 use the IDENTICAL real-only + * decomposition, with a genuinely complex decomposition (step 1) + * processed on the same instance/eigenIndex in between. If isAllReal + * gating left any stale sDs/sCs kernel state or stale device Eval data + * from the complex step bleeding into the next real step, step 2 would + * diverge from step 0 -- this is the most direct test for that failure + * mode, independent of any CPU-vs-GPU precision noise. + * Tolerance is not zero: verified empirically (by rerunning with all + * three steps set to the same real decomposition, i.e. no complex step + * at all) that back-to-back identical GPU launches already differ by + * this same ~1e-5 magnitude on this hardware -- ADJOINT_ATOMIC_ADD_GPU's + * float addition is not associative, so summation order across + * concurrent atomics is not guaranteed bit-identical between launches. + * That noise floor is unrelated to eigen real/complex gating, so the + * threshold here is set an order of magnitude above it. */ + double stepConsistencyDiff = 0.0; + for (int k = 0; k < S*S; k++) + stepConsistencyDiff = std::max(stepConsistencyDiff, fabs(gpuGradSteps[0][k] - gpuGradSteps[2][k])); + bool consistencyOk = stepConsistencyDiff < 1e-4; + fprintf(stdout, " GPU step0 vs step2 (both real, complex processed in between) max|diff|=%.3e %s\n", + stepConsistencyDiff, consistencyOk ? "PASS (no stale state)" : "FAIL (possible stale state)"); + allOk &= consistencyOk; + beagleFinalizeInstance(gpuInst); + + fprintf(stdout, "Dynamic real<->complex transition test: %s\n", allOk ? "PASS" : "FAIL"); + return allOk ? 0 : 1; +} + /* ── main ────────────────────────────────────────────────────────────── */ int main(int argc, const char *argv[]) { @@ -941,9 +1258,11 @@ int main(int argc, const char *argv[]) { bool useGpu = false; int whichDevice = -1; - bool run4 = true; - bool run16 = false; - bool run17 = false; + bool run4 = true; + bool run16 = false; + bool run17 = false; + bool runRealOnly = false; + bool runDynamic = false; for (int i = 1; i < argc; i++) { if (strcmp(argv[i], "--gpu") == 0 && i+1 < argc) { @@ -954,17 +1273,25 @@ int main(int argc, const char *argv[]) { run4 = (ns == 4); run16 = (ns == 16); run17 = (ns == 17); + } else if (strcmp(argv[i], "--realonly") == 0) { + run4 = false; + runRealOnly = true; + } else if (strcmp(argv[i], "--dynamic") == 0) { + run4 = false; + runDynamic = true; } else if (strcmp(argv[i], "--all") == 0) { - run4 = run16 = run17 = true; + run4 = run16 = run17 = runRealOnly = runDynamic = true; } } int gpuDev = useGpu ? whichDevice : -1; int result = 0; - if (run4) result |= runTest4(gpuDev); - if (run16) result |= runTest16(gpuDev); - if (run17) result |= runTest17(gpuDev); + if (run4) result |= runTest4(gpuDev); + if (run16) result |= runTest16(gpuDev); + if (run17) result |= runTest17(gpuDev); + if (runRealOnly) result |= runTestRealOnly(gpuDev); + if (runDynamic) result |= runTestDynamicTransition(gpuDev); return result; } diff --git a/libhmsbeagle/GPU/BeagleGPUImpl.hpp b/libhmsbeagle/GPU/BeagleGPUImpl.hpp index 234a8fa6..f8eac334 100644 --- a/libhmsbeagle/GPU/BeagleGPUImpl.hpp +++ b/libhmsbeagle/GPU/BeagleGPUImpl.hpp @@ -505,12 +505,15 @@ int BeagleGPUImpl::createInstance(int tipCount, kPartialsSize = kPaddedPatternCount * kPaddedStateCount * kCategoryCount; kMatrixSize = kPaddedStateCount * kPaddedStateCount; - // Spectral kernels (kernelsSpectralIfDef*.cu / kernelsSpectral*.cu) always read a - // real+imaginary pair per eigenstate (SPECTRAL_EIGENVALS_GPU unconditionally indexes - // eigenValues[PADDED_STATE_COUNT + state]) regardless of BEAGLE_FLAG_EIGEN_COMPLEX, so the - // device buffer must be sized/copied as complex whenever spectral representation is in use, - // otherwise that read runs past the data actually written by setEigenDecomposition. - if ((kFlags & BEAGLE_FLAG_EIGEN_COMPLEX) || (kFlags & BEAGLE_FLAG_SPECTRAL_REPRESENTATION)) + // Spectral kernels (kernelsSpectralIfDef*.cu / kernelsSpectral*.cu) read a real+imaginary + // pair per eigenstate only when that eigendecomposition is not all-real + // (SPECTRAL_EIGENVALS_GPU is now gated on the isAllReal1/isAllReal2 kernel arguments, threaded + // from hEigenDecompIsAllReal — see BeagleGPUSpectralImpl.hpp::setEigenDecomposition), matching + // the CPU implementation's EigenDecompositionSpectral.hpp, which only widens on + // BEAGLE_FLAG_EIGEN_COMPLEX. Widening for BEAGLE_FLAG_SPECTRAL_REPRESENTATION alone (regardless + // of EIGEN_COMPLEX) has no CPU analog and would waste memory/copy work for spectral-mode + // instances that are structurally guaranteed real, so this only widens for EIGEN_COMPLEX. + if (kFlags & BEAGLE_FLAG_EIGEN_COMPLEX) kEigenValuesSize = 2 * kPaddedStateCount; else kEigenValuesSize = kPaddedStateCount; @@ -4756,7 +4759,6 @@ void BeagleGPUImpl::dispatchPruneSS(GPUPtr s1, GPUPtr s2, GP GPUPtr scalingFactors, GPUPtr cumulativeScaling, unsigned int startPattern, unsigned int endPattern, int rescale, int streamIndex, int waitIndex) { - fprintf(stderr, "[DISPATCH] BASE BeagleGPUImpl::dispatchPruneSS called\n"); fflush(stderr); kernels->StatesStatesPruningDynamicScaling(s1, s2, p3, dMatrices[c1MatIdx], dMatrices[c2MatIdx], scalingFactors, cumulativeScaling, diff --git a/libhmsbeagle/GPU/BeagleGPUSpectralImpl.h b/libhmsbeagle/GPU/BeagleGPUSpectralImpl.h index 56c8baa9..607cd72b 100644 --- a/libhmsbeagle/GPU/BeagleGPUSpectralImpl.h +++ b/libhmsbeagle/GPU/BeagleGPUSpectralImpl.h @@ -45,6 +45,22 @@ class BeagleGPUSpectralImpl : public BeagleGPUImpl { * adjoint offset-queue) must multiply by this, not by kCategoryCount. */ unsigned int kSpectralDistanceStrideElements; + /* Gather-batch-scatter queue for updateTransitionMatrices' per-branch + * distance uploads: one flat host-side (destination-offset, value) pair + * per (branch, category), uploaded with two MemcpyHostToDevice calls + * regardless of branch count, then scattered into dSpectralDistancesOrigin + * by one GPU kernel -- mirrors BeagleGPUImpl's own dPtrQueue/dDistanceQueue + * pattern for its dense-matrix update. Dedicated buffers rather than + * reusing that base-class pair: those are declared in BeagleGPUImpl's + * `private:` section, not `protected:`, so this derived class cannot + * reach them without changing the base class's access specifiers. + * Sized kMatrixCount * kCategoryCount, the largest a single + * updateTransitionMatrices call can need. */ + GPUPtr dSpectralPtrQueue; + GPUPtr dSpectralDistanceQueue; + unsigned int* hSpectralPtrQueue; + Real* hSpectralDistanceQueue; + /* Backward eigenvector buffers for the parent branch in pre-order. * dEvecT[ei] = U stored row-major (dEvecT[j*S+k] = U[j,k]). * Used as ievc1 in Growing kernels (backward Phase 1: U^T·p). diff --git a/libhmsbeagle/GPU/BeagleGPUSpectralImpl.hpp b/libhmsbeagle/GPU/BeagleGPUSpectralImpl.hpp index 10e7a418..34bf34d0 100644 --- a/libhmsbeagle/GPU/BeagleGPUSpectralImpl.hpp +++ b/libhmsbeagle/GPU/BeagleGPUSpectralImpl.hpp @@ -34,6 +34,8 @@ namespace gpu { BEAGLE_GPU_TEMPLATE BeagleGPUSpectralImpl::BeagleGPUSpectralImpl() : dSpectralDistancesOrigin(0), dSpectralDistances(NULL), hEigenIndexForMatrix(NULL), + dSpectralPtrQueue(0), dSpectralDistanceQueue(0), + hSpectralPtrQueue(NULL), hSpectralDistanceQueue(NULL), dEvecTOrigin(0), dEvecT(NULL), dIevcTOrigin(0), dIevcT(NULL), dGradientOrigin(0), dGradient(NULL), hEigenDecompIsAllReal(NULL), dAdjointQueue(0), hAdjointQueue(NULL), kAdjointQueueCapacity(0), @@ -45,6 +47,10 @@ BeagleGPUSpectralImpl::~BeagleGPUSpectralImpl() { GPUInterface* gpuIf = this->gpu; if (gpuIf) { if (dSpectralDistancesOrigin) gpuIf->FreeMemory(dSpectralDistancesOrigin); + if (dSpectralPtrQueue) gpuIf->FreeMemory(dSpectralPtrQueue); + if (dSpectralDistanceQueue) gpuIf->FreeMemory(dSpectralDistanceQueue); + if (hSpectralPtrQueue) gpuIf->FreeHostMemory(hSpectralPtrQueue); + if (hSpectralDistanceQueue) gpuIf->FreeHostMemory(hSpectralDistanceQueue); if (dEvecTOrigin) gpuIf->FreeMemory(dEvecTOrigin); if (dIevcTOrigin) gpuIf->FreeMemory(dIevcTOrigin); if (dGradientOrigin) gpuIf->FreeMemory(dGradientOrigin); @@ -95,6 +101,17 @@ int BeagleGPUSpectralImpl::createInstance( dSpectralDistances[i] = gpuIf->CreateSubPointer(dSpectralDistancesOrigin, distStride * i, distStride); } + // Gather-batch-scatter queue for the per-branch distance upload in + // updateTransitionMatrices (see BeagleGPUSpectralImpl.h for rationale). + // Sized for the largest possible single call: every matrix, every category. + int spectralQueueLength = this->kMatrixCount * this->kCategoryCount; + dSpectralPtrQueue = gpuIf->AllocateMemory(sizeof(unsigned int) * spectralQueueLength); + hSpectralPtrQueue = (unsigned int*) gpuIf->MallocHost(sizeof(unsigned int) * spectralQueueLength); + checkHostMemory(hSpectralPtrQueue); + dSpectralDistanceQueue = gpuIf->AllocateMemory(sizeof(Real) * spectralQueueLength); + hSpectralDistanceQueue = (Real*) gpuIf->MallocHost(sizeof(Real) * spectralQueueLength); + checkHostMemory(hSpectralDistanceQueue); + // Backward eigenvector arrays: one S*S block per eigen decomposition. int S = this->kPaddedStateCount; size_t matStride = gpuIf->AlignMemOffset((size_t)S * S * sizeof(Real)); @@ -136,17 +153,21 @@ int BeagleGPUSpectralImpl::updateTransitionMatrices( if (count > 0) { const double* categoryRates = this->hCategoryRates[0]; int nCat = this->kCategoryCount; - Real* hDist = (Real*) malloc(sizeof(Real) * nCat); GPUInterface* gpuIf = this->gpu; + int totalCount = 0; for (int i = 0; i < count; i++) { int matIdx = probabilityIndices[i]; hEigenIndexForMatrix[matIdx] = eigenIndex; for (int j = 0; j < nCat; j++) { - hDist[j] = (Real)(edgeLengths[i] * categoryRates[j]); + hSpectralPtrQueue[totalCount] = matIdx * kSpectralDistanceStrideElements + j; + hSpectralDistanceQueue[totalCount] = (Real)(edgeLengths[i] * categoryRates[j]); + totalCount++; } - gpuIf->MemcpyHostToDevice(dSpectralDistances[matIdx], hDist, sizeof(Real) * nCat); } - free(hDist); + gpuIf->MemcpyHostToDevice(dSpectralPtrQueue, hSpectralPtrQueue, sizeof(unsigned int) * totalCount); + gpuIf->MemcpyHostToDevice(dSpectralDistanceQueue, hSpectralDistanceQueue, sizeof(Real) * totalCount); + this->kernels->ScatterSpectralDistances(dSpectralDistancesOrigin, dSpectralPtrQueue, + dSpectralDistanceQueue, totalCount); } return BEAGLE_SUCCESS; } @@ -165,6 +186,7 @@ void BeagleGPUSpectralImpl::dispatchPrunePP( this->dIevc[ei2], this->dEvec[ei2], this->dEigenValues[ei2], dSpectralDistances[c2MatIdx], scalingFactors, cumulativeScaling, this->kPaddedPatternCount, this->kCategoryCount, + hEigenDecompIsAllReal[ei1] ? 1 : 0, hEigenDecompIsAllReal[ei2] ? 1 : 0, rescale, streamIndex, waitIndex); } @@ -182,6 +204,7 @@ void BeagleGPUSpectralImpl::dispatchPruneSP( this->dIevc[ei2], this->dEvec[ei2], this->dEigenValues[ei2], dSpectralDistances[c2MatIdx], scalingFactors, cumulativeScaling, this->kPaddedPatternCount, this->kCategoryCount, + hEigenDecompIsAllReal[ei1] ? 1 : 0, hEigenDecompIsAllReal[ei2] ? 1 : 0, rescale, streamIndex, waitIndex); } @@ -192,7 +215,6 @@ void BeagleGPUSpectralImpl::dispatchPruneSS( GPUPtr scalingFactors, GPUPtr cumulativeScaling, unsigned int startPattern, unsigned int endPattern, int rescale, int streamIndex, int waitIndex) { - fprintf(stderr, "[DISPATCH] SPECTRAL BeagleGPUSpectralImpl::dispatchPruneSS called\n"); fflush(stderr); int ei1 = hEigenIndexForMatrix[c1MatIdx]; int ei2 = hEigenIndexForMatrix[c2MatIdx]; this->kernels->StatesStatesPruningSpectral(s1, s2, p3, @@ -200,6 +222,7 @@ void BeagleGPUSpectralImpl::dispatchPruneSS( this->dIevc[ei2], this->dEvec[ei2], this->dEigenValues[ei2], dSpectralDistances[c2MatIdx], scalingFactors, cumulativeScaling, this->kPaddedPatternCount, this->kCategoryCount, + hEigenDecompIsAllReal[ei1] ? 1 : 0, hEigenDecompIsAllReal[ei2] ? 1 : 0, rescale, streamIndex, waitIndex); } @@ -217,12 +240,12 @@ int BeagleGPUSpectralImpl::setEigenDecomposition( const int SC = this->kStateCount; const int SS = S * S; // kEigenValuesSize is private to BeagleGPUImpl; recompute the same way it does. - // Spectral kernels always read a real+imaginary pair per eigenstate regardless of - // BEAGLE_FLAG_EIGEN_COMPLEX (see matching comment/fix in BeagleGPUImpl::createInstance), - // so this local copy must widen for BEAGLE_FLAG_SPECTRAL_REPRESENTATION too, or the - // imaginary half of the device buffer is left uninitialized. - const int eigenValuesSize = ((this->kFlags & BEAGLE_FLAG_EIGEN_COMPLEX) || - (this->kFlags & BEAGLE_FLAG_SPECTRAL_REPRESENTATION)) ? 2 * S : S; + // Spectral kernels read a real+imaginary pair per eigenstate only when that decomposition is + // not all-real (see matching comment/fix in BeagleGPUImpl::createInstance) -- gated at the + // kernel level on isAllReal1/isAllReal2, computed below from inEigenValues and cached in + // hEigenDecompIsAllReal. This local copy only widens for BEAGLE_FLAG_EIGEN_COMPLEX, matching + // CPU's EigenDecompositionSpectral.hpp. + const int eigenValuesSize = (this->kFlags & BEAGLE_FLAG_EIGEN_COMPLEX) ? 2 * S : S; // Forward: dEvec[row*S+col] = U[col,row]; dIevc[row*S+col] = U^-1[col,row] // Backward: dEvecT[row*S+col] = U[row,col]; dIevcT[row*S+col] = U^-1[row,col] @@ -320,18 +343,21 @@ void BeagleGPUSpectralImpl::dispatchGrowingSpectral( GPUPtr evec2 = this->dEvec[ei2]; GPUPtr eval2 = this->dEigenValues[ei2]; GPUPtr dist2 = dSpectralDistances[c2MatIdx]; + int isAllReal2 = hEigenDecompIsAllReal[ei2] ? 1 : 0; if (isRoot) { if (sibIsStates) { this->kernels->PartialsStatesGrowingSpectralTopRoot( partials1, c2, partials3, ievc2, evec2, eval2, dist2, - this->kPaddedPatternCount, this->kCategoryCount); + this->kPaddedPatternCount, this->kCategoryCount, + isAllReal2); } else { this->kernels->PartialsPartialsGrowingSpectralTopRoot( partials1, c2, partials3, ievc2, evec2, eval2, dist2, - this->kPaddedPatternCount, this->kCategoryCount); + this->kPaddedPatternCount, this->kCategoryCount, + isAllReal2); } } else { int ei1 = hEigenIndexForMatrix[c1MatIdx]; @@ -341,6 +367,7 @@ void BeagleGPUSpectralImpl::dispatchGrowingSpectral( GPUPtr bEvec1 = dIevcT[ei1]; GPUPtr eval1 = this->dEigenValues[ei1]; GPUPtr dist1 = dSpectralDistances[c1MatIdx]; + int isAllReal1 = hEigenDecompIsAllReal[ei1] ? 1 : 0; if (isTop) { if (sibIsStates) { @@ -348,7 +375,8 @@ void BeagleGPUSpectralImpl::dispatchGrowingSpectral( partials1, c2, partials3, bIevc1, bEvec1, eval1, dist1, ievc2, evec2, eval2, dist2, - this->kPaddedPatternCount, this->kCategoryCount); + this->kPaddedPatternCount, this->kCategoryCount, + isAllReal1, isAllReal2); } else { // Top NotRoot PP: reuse no-scale PP pruning kernel with backward matrices this->kernels->PartialsPartialsPruningSpectral( @@ -357,6 +385,7 @@ void BeagleGPUSpectralImpl::dispatchGrowingSpectral( ievc2, evec2, eval2, dist2, nullptr, nullptr, this->kPaddedPatternCount, this->kCategoryCount, + isAllReal1, isAllReal2, -1, -1, -1); } } else { @@ -365,13 +394,15 @@ void BeagleGPUSpectralImpl::dispatchGrowingSpectral( partials1, c2, partials3, bIevc1, bEvec1, eval1, dist1, ievc2, evec2, eval2, dist2, - this->kPaddedPatternCount, this->kCategoryCount); + this->kPaddedPatternCount, this->kCategoryCount, + isAllReal1, isAllReal2); } else { this->kernels->PartialsPartialsGrowingSpectral( partials1, c2, partials3, bIevc1, bEvec1, eval1, dist1, ievc2, evec2, eval2, dist2, - this->kPaddedPatternCount, this->kCategoryCount); + this->kPaddedPatternCount, this->kCategoryCount, + isAllReal1, isAllReal2); } } } diff --git a/libhmsbeagle/GPU/GPUInterfaceOpenCL.cpp b/libhmsbeagle/GPU/GPUInterfaceOpenCL.cpp index f26f11e0..eae6124d 100644 --- a/libhmsbeagle/GPU/GPUInterfaceOpenCL.cpp +++ b/libhmsbeagle/GPU/GPUInterfaceOpenCL.cpp @@ -22,6 +22,7 @@ #include #include #include +#include #include "libhmsbeagle/beagle.h" #include "libhmsbeagle/GPU/GPUImplDefs.h" @@ -29,6 +30,93 @@ #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; + long memcpyD2HCalls = 0; + double memcpyD2HSeconds = 0.0; + size_t memcpyD2HBytes = 0; + long memsetZeroCalls = 0; + double memsetZeroSeconds = 0.0; + size_t memsetZeroBytes = 0; +}; + +BeagleProfileStats gBeagleProfile; + +void beagleProfilePrintSummary() { + if (!gBeagleProfile.enabled) return; + BeagleProfileStats& p = gBeagleProfile; + double totalAccounted = p.launchKernelSeconds + p.launchKernelConcurrentSeconds + + 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", + 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, "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); +} + +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 +565,24 @@ void GPUInterface::SynchronizeHost() { // SAFE_CL(clFinish(openClCommandQueues[i])); // } + bool prof = beagleProfilingEnabled(); + std::chrono::high_resolution_clock::time_point t0; + if (prof) t0 = std::chrono::high_resolution_clock::now(); + SAFE_CL(clFinish(openClCommandQueues[0])); + if (prof) { + double dt = std::chrono::duration( + std::chrono::high_resolution_clock::now() - t0).count(); + if (gBeagleProfile.warmupRemaining > 0) { + gBeagleProfile.warmupRemaining--; + if (gBeagleProfile.warmupRemaining == 0) gBeagleProfile.warm = true; + } else { + gBeagleProfile.syncHostCalls++; + gBeagleProfile.syncHostSeconds += dt; + } + } + #ifdef BEAGLE_DEBUG_FLOW fprintf(stderr,"\t\t\tLeaving GPUInterface::SynchronizeHost\n"); #endif @@ -545,17 +649,24 @@ 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")) { + // 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 +675,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 +695,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 +713,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 +725,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 +741,14 @@ void GPUInterface::LaunchKernel(GPUFunction deviceFunction, globalWorkSize, localWorkSize, 0, NULL, NULL)); } - if (getenv("BEAGLE_DEBUG_KERNEL_ARGS")) { + if (prof && gBeagleProfile.warm) { + double dt = std::chrono::duration( + std::chrono::high_resolution_clock::now() - t0).count(); + gBeagleProfile.launchKernelCalls++; + gBeagleProfile.launchKernelSeconds += dt; + } + + if (debugArgs) { fprintf(stderr, "[LaunchKernel] %s: enqueued OK, calling clFinish\n", kernelNameBuf); fflush(stderr); clFinish(openClCommandQueues[0]); @@ -656,17 +774,24 @@ 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")) { + // 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 +800,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 +820,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 +845,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 +889,14 @@ void GPUInterface::LaunchKernelConcurrent(GPUFunction deviceFunction, globalWorkSize, localWorkSize, 0, NULL, NULL)); // } - if (getenv("BEAGLE_DEBUG_KERNEL_ARGS")) { + if (prof && gBeagleProfile.warm) { + double dt = std::chrono::duration( + std::chrono::high_resolution_clock::now() - t0).count(); + gBeagleProfile.launchKernelConcurrentCalls++; + gBeagleProfile.launchKernelConcurrentSeconds += dt; + } + + if (debugArgs) { fprintf(stderr, "[LaunchKernelConcurrent] %s: enqueued OK, calling clFinish\n", kernelNameBuf); fflush(stderr); clFinish(openClCommandQueues[0]); @@ -1006,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 @@ -1021,6 +1157,14 @@ void GPUInterface::MemsetZero(GPUPtr dest, NULL, NULL)); free(hZero); + if (prof && gBeagleProfile.warm) { + double dt = std::chrono::duration( + std::chrono::high_resolution_clock::now() - t0).count(); + gBeagleProfile.memsetZeroCalls++; + gBeagleProfile.memsetZeroSeconds += dt; + gBeagleProfile.memsetZeroBytes += memSize; + } + #ifdef BEAGLE_DEBUG_FLOW fprintf(stderr, "\t\t\tLeaving GPUInterface::MemsetZero\n"); #endif @@ -1034,9 +1178,23 @@ void GPUInterface::MemcpyHostToDevice(GPUPtr dest, fprintf(stderr, "\t\t\tEntering GPUInterface::MemcpyHostToDevice\n"); #endif + bool prof = beagleProfilingEnabled(); + std::chrono::high_resolution_clock::time_point t0; + if (prof) t0 = std::chrono::high_resolution_clock::now(); + + // CL_TRUE => blocking write: this call does not return until the transfer + // completes, regardless of payload size. SAFE_CL(clEnqueueWriteBuffer(openClCommandQueues[0], dest, CL_TRUE, 0, memSize, src, 0, NULL, NULL)); + if (prof && gBeagleProfile.warm) { + double dt = std::chrono::duration( + std::chrono::high_resolution_clock::now() - t0).count(); + gBeagleProfile.memcpyH2DCalls++; + gBeagleProfile.memcpyH2DSeconds += dt; + gBeagleProfile.memcpyH2DBytes += memSize; + } + #ifdef BEAGLE_DEBUG_FLOW fprintf(stderr, "\t\t\tLeaving GPUInterface::MemcpyHostToDevice\n"); #endif @@ -1049,9 +1207,23 @@ void GPUInterface::MemcpyDeviceToHost(void* dest, fprintf(stderr, "\t\t\tEntering GPUInterface::MemcpyDeviceToHost\n"); #endif + bool prof = beagleProfilingEnabled(); + std::chrono::high_resolution_clock::time_point t0; + if (prof) t0 = std::chrono::high_resolution_clock::now(); + + // CL_TRUE => blocking read: this call does not return until the transfer + // completes, regardless of payload size. SAFE_CL(clEnqueueReadBuffer(openClCommandQueues[0], src, CL_TRUE, 0, memSize, dest, 0, NULL, NULL)); + if (prof && gBeagleProfile.warm) { + double dt = std::chrono::duration( + std::chrono::high_resolution_clock::now() - t0).count(); + gBeagleProfile.memcpyD2HCalls++; + gBeagleProfile.memcpyD2HSeconds += dt; + gBeagleProfile.memcpyD2HBytes += memSize; + } + #ifdef BEAGLE_DEBUG_FLOW fprintf(stderr, "\t\t\tLeaving GPUInterface::MemcpyDeviceToHost\n"); #endif diff --git a/libhmsbeagle/GPU/KernelLauncher.cpp b/libhmsbeagle/GPU/KernelLauncher.cpp index f9de48c0..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( @@ -1041,13 +1043,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 +1065,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 +1087,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 +1107,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 +1125,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(); } @@ -1171,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, @@ -1779,26 +1815,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 +1848,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 +1869,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 +1878,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 +1892,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 +1913,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 +1922,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 +1936,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..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; @@ -320,6 +326,8 @@ class KernelLauncher { GPUPtr cumulativeScaling, unsigned int patternCount, unsigned int categoryCount, + int isAllReal1, + int isAllReal2, int doRescaling, int streamIndex, int waitIndex); @@ -335,6 +343,8 @@ class KernelLauncher { GPUPtr cumulativeScaling, unsigned int patternCount, unsigned int categoryCount, + int isAllReal1, + int isAllReal2, int doRescaling, int streamIndex, int waitIndex); @@ -350,6 +360,8 @@ class KernelLauncher { GPUPtr cumulativeScaling, unsigned int patternCount, unsigned int categoryCount, + int isAllReal1, + int isAllReal2, int doRescaling, int streamIndex, int waitIndex); @@ -386,7 +398,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 +410,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 +422,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 +432,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 +441,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 @@ -459,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 96feff1c..038d1291 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() @@ -1290,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 1ca91d79..2cf39f7d 100644 --- a/libhmsbeagle/GPU/kernels/kernelsSpectralIfDef4.cu +++ b/libhmsbeagle/GPU/kernels/kernelsSpectralIfDef4.cu @@ -152,23 +152,41 @@ /* patIdx-0 threads: one thread per eigenstate (state = k). * eigenValues layout: [realEV_0..realEV_{S-1} | imagEV_0..imagEV_{S-1}] * distances[matrix] = branchLength * categoryRate[matrix]. + * ISALLREAL1/ISALLREAL2: per-eigendecomposition runtime flags (from + * hEigenDecompIsAllReal, threaded through as ordinary kernel arguments — see + * kernelAdjointMerged4's isAllReal for the precedent this mirrors). When a + * decomposition is all-real, the imaginary half of eigenValues[] is never + * read: it may not even be allocated (kEigenValuesSize is only widened for + * BEAGLE_FLAG_EIGEN_COMPLEX, not merely for spectral-representation use), so + * reading it unconditionally is an out-of-bounds access, not just wasted + * work. sCs == 0 (the real-eigenvalue convention already documented on + * SPECTRAL_COMMON_SMEM_GPU above) is produced explicitly rather than relying + * on cos(0)=1/sin(0)=0 falling out of an unread angle. * Ends with KW_LOCAL_FENCE so sP*, sScale, sDs*, sCs* are all visible. */ -#define SPECTRAL_EIGENVALS_GPU() \ +#define SPECTRAL_EIGENVALS_GPU(ISALLREAL1, ISALLREAL2) \ if (patIdx == 0) { \ REAL t1 = distances1[matrix]; \ REAL e1 = exp(eigenValues1[state] * t1); \ - REAL bt1 = eigenValues1[PADDED_STATE_COUNT + state] * t1; \ - REAL cv1, sv1; \ - SPECTRAL_SINCOS(bt1, sv1, cv1); \ - sDs1[state] = e1 * cv1; \ - sCs1[state] = e1 * sv1; \ + sDs1[state] = e1; \ + sCs1[state] = (REAL)0; \ + if (!(ISALLREAL1)) { \ + REAL bt1 = eigenValues1[PADDED_STATE_COUNT + state] * t1; \ + REAL cv1, sv1; \ + SPECTRAL_SINCOS(bt1, sv1, cv1); \ + sDs1[state] = e1 * cv1; \ + sCs1[state] = e1 * sv1; \ + } \ REAL t2 = distances2[matrix]; \ REAL e2 = exp(eigenValues2[state] * t2); \ - REAL bt2 = eigenValues2[PADDED_STATE_COUNT + state] * t2; \ - REAL cv2, sv2; \ - SPECTRAL_SINCOS(bt2, sv2, cv2); \ - sDs2[state] = e2 * cv2; \ - sCs2[state] = e2 * sv2; \ + sDs2[state] = e2; \ + sCs2[state] = (REAL)0; \ + if (!(ISALLREAL2)) { \ + REAL bt2 = eigenValues2[PADDED_STATE_COUNT + state] * t2; \ + REAL cv2, sv2; \ + SPECTRAL_SINCOS(bt2, sv2, cv2); \ + sDs2[state] = e2 * cv2; \ + sCs2[state] = e2 * sv2; \ + } \ } \ KW_LOCAL_FENCE; @@ -375,6 +393,7 @@ KW_DEVICE_FUNC void kernelSpectralBody( KW_GLOBAL_VAR REAL* KW_RESTRICT eigenValues2, KW_GLOBAL_VAR REAL* KW_RESTRICT distances2, KW_GLOBAL_VAR REAL* KW_RESTRICT scalingFactors, /* NULL if !USE_SCALING */ + int isAllReal1, int isAllReal2, int totalPatterns) { SPECTRAL_INDICES_GPU() @@ -392,7 +411,7 @@ KW_DEVICE_FUNC void kernelSpectralBody( SPECTRAL_LOAD_SCALE_GPU() #endif - SPECTRAL_EIGENVALS_GPU() /* ends with KW_LOCAL_FENCE */ + SPECTRAL_EIGENVALS_GPU(isAllReal1, isAllReal2) /* ends with KW_LOCAL_FENCE */ /* Phase 1: project to eigenspace. * #ifdef replaces if constexpr (is_same) from the C++ version. */ @@ -444,12 +463,13 @@ KW_GLOBAL_KERNEL void kernelPartialsPartialsNoScaleSpectral( KW_GLOBAL_VAR REAL* KW_RESTRICT evec2, KW_GLOBAL_VAR REAL* KW_RESTRICT eigenValues2, KW_GLOBAL_VAR REAL* KW_RESTRICT distances2, + int isAllReal1, int isAllReal2, int totalPatterns) { SPECTRAL_INDICES_GPU() SPECTRAL_COMMON_SMEM_GPU() SPECTRAL_LOAD_PARTIALS1_GPU() SPECTRAL_LOAD_PARTIALS2_GPU() - SPECTRAL_EIGENVALS_GPU() + SPECTRAL_EIGENVALS_GPU(isAllReal1, isAllReal2) SPECTRAL_PHASE1_PARTIALS_DUAL_GPU(sBuf1, sP1, sQ1, ievc1, sBuf2, sP2, sQ2, ievc2) SPECTRAL_PHASE2_GPU() SPECTRAL_PHASE3_GPU() @@ -469,13 +489,14 @@ KW_GLOBAL_KERNEL void kernelPartialsPartialsFixedScaleSpectral( KW_GLOBAL_VAR REAL* KW_RESTRICT eigenValues2, KW_GLOBAL_VAR REAL* KW_RESTRICT distances2, KW_GLOBAL_VAR REAL* KW_RESTRICT scalingFactors, + int isAllReal1, int isAllReal2, int totalPatterns) { SPECTRAL_INDICES_GPU() SPECTRAL_COMMON_SMEM_GPU() SPECTRAL_LOAD_PARTIALS1_GPU() SPECTRAL_LOAD_PARTIALS2_GPU() SPECTRAL_LOAD_SCALE_GPU() - SPECTRAL_EIGENVALS_GPU() + SPECTRAL_EIGENVALS_GPU(isAllReal1, isAllReal2) SPECTRAL_PHASE1_PARTIALS_DUAL_GPU(sBuf1, sP1, sQ1, ievc1, sBuf2, sP2, sQ2, ievc2) SPECTRAL_PHASE2_GPU() SPECTRAL_PHASE3_GPU() @@ -498,11 +519,12 @@ KW_GLOBAL_KERNEL void kernelStatesPartialsNoScaleSpectral( KW_GLOBAL_VAR REAL* KW_RESTRICT evec2, KW_GLOBAL_VAR REAL* KW_RESTRICT eigenValues2, KW_GLOBAL_VAR REAL* KW_RESTRICT distances2, + int isAllReal1, int isAllReal2, int totalPatterns) { SPECTRAL_INDICES_GPU() SPECTRAL_COMMON_SMEM_GPU() SPECTRAL_LOAD_PARTIALS2_GPU() /* no sP1: child 1 is States */ - SPECTRAL_EIGENVALS_GPU() + SPECTRAL_EIGENVALS_GPU(isAllReal1, isAllReal2) SPECTRAL_PHASE1_STATES_GPU(sQ1, ievc1, states1) SPECTRAL_PHASE1_PARTIALS_GPU(sBuf2, sP2, sQ2, ievc2) SPECTRAL_PHASE2_GPU() @@ -523,12 +545,13 @@ KW_GLOBAL_KERNEL void kernelStatesPartialsFixedScaleSpectral( KW_GLOBAL_VAR REAL* KW_RESTRICT eigenValues2, KW_GLOBAL_VAR REAL* KW_RESTRICT distances2, KW_GLOBAL_VAR REAL* KW_RESTRICT scalingFactors, + int isAllReal1, int isAllReal2, int totalPatterns) { SPECTRAL_INDICES_GPU() SPECTRAL_COMMON_SMEM_GPU() SPECTRAL_LOAD_PARTIALS2_GPU() SPECTRAL_LOAD_SCALE_GPU() - SPECTRAL_EIGENVALS_GPU() + SPECTRAL_EIGENVALS_GPU(isAllReal1, isAllReal2) SPECTRAL_PHASE1_STATES_GPU(sQ1, ievc1, states1) SPECTRAL_PHASE1_PARTIALS_GPU(sBuf2, sP2, sQ2, ievc2) SPECTRAL_PHASE2_GPU() @@ -551,10 +574,11 @@ KW_GLOBAL_KERNEL void kernelStatesStatesNoScaleSpectral( KW_GLOBAL_VAR REAL* KW_RESTRICT evec2, KW_GLOBAL_VAR REAL* KW_RESTRICT eigenValues2, KW_GLOBAL_VAR REAL* KW_RESTRICT distances2, + int isAllReal1, int isAllReal2, int totalPatterns) { SPECTRAL_INDICES_GPU() SPECTRAL_COMMON_SMEM_GPU() /* no LOAD_PARTIALS: both children are States */ - SPECTRAL_EIGENVALS_GPU() + SPECTRAL_EIGENVALS_GPU(isAllReal1, isAllReal2) SPECTRAL_PHASE1_STATES_GPU(sQ1, ievc1, states1) SPECTRAL_PHASE1_STATES_GPU(sQ2, ievc2, states2) SPECTRAL_PHASE2_GPU() @@ -575,11 +599,12 @@ KW_GLOBAL_KERNEL void kernelStatesStatesFixedScaleSpectral( KW_GLOBAL_VAR REAL* KW_RESTRICT eigenValues2, KW_GLOBAL_VAR REAL* KW_RESTRICT distances2, KW_GLOBAL_VAR REAL* KW_RESTRICT scalingFactors, + int isAllReal1, int isAllReal2, int totalPatterns) { SPECTRAL_INDICES_GPU() SPECTRAL_COMMON_SMEM_GPU() SPECTRAL_LOAD_SCALE_GPU() - SPECTRAL_EIGENVALS_GPU() + SPECTRAL_EIGENVALS_GPU(isAllReal1, isAllReal2) SPECTRAL_PHASE1_STATES_GPU(sQ1, ievc1, states1) SPECTRAL_PHASE1_STATES_GPU(sQ2, ievc2, states2) SPECTRAL_PHASE2_GPU() @@ -664,16 +689,21 @@ KW_GLOBAL_KERNEL void kernelStatesStatesFixedScaleSpectral( } /* Load only sibling (child2) eigenvalue exponentials — for Top Root kernels - * where child1 is the root and carries no branch transform. */ -#define SPECTRAL_EIGENVALS_SIB_ONLY_GPU() \ + * where child1 is the root and carries no branch transform. + * ISALLREAL2: see SPECTRAL_EIGENVALS_GPU's comment above — same gating. */ +#define SPECTRAL_EIGENVALS_SIB_ONLY_GPU(ISALLREAL2) \ if (patIdx == 0) { \ REAL t2 = distances2[matrix]; \ REAL e2 = exp(eigenValues2[state] * t2); \ - REAL bt2 = eigenValues2[PADDED_STATE_COUNT + state] * t2; \ - REAL cv2, sv2; \ - SPECTRAL_SINCOS(bt2, sv2, cv2); \ - sDs2[state] = e2 * cv2; \ - sCs2[state] = e2 * sv2; \ + sDs2[state] = e2; \ + sCs2[state] = (REAL)0; \ + if (!(ISALLREAL2)) { \ + REAL bt2 = eigenValues2[PADDED_STATE_COUNT + state] * t2; \ + REAL cv2, sv2; \ + SPECTRAL_SINCOS(bt2, sv2, cv2); \ + sDs2[state] = e2 * cv2; \ + sCs2[state] = e2 * sv2; \ + } \ } \ KW_LOCAL_FENCE; @@ -708,12 +738,13 @@ KW_GLOBAL_KERNEL void kernelPartialsPartialsAutoScaleSpectral( KW_GLOBAL_VAR REAL* KW_RESTRICT eigenValues2, KW_GLOBAL_VAR REAL* KW_RESTRICT distances2, KW_GLOBAL_VAR signed char* KW_RESTRICT scalingFactors, + int isAllReal1, int isAllReal2, int totalPatterns) { SPECTRAL_INDICES_GPU() SPECTRAL_COMMON_SMEM_GPU() SPECTRAL_LOAD_PARTIALS1_GPU() SPECTRAL_LOAD_PARTIALS2_GPU() - SPECTRAL_EIGENVALS_GPU() + SPECTRAL_EIGENVALS_GPU(isAllReal1, isAllReal2) SPECTRAL_PHASE1_PARTIALS_DUAL_GPU(sBuf1, sP1, sQ1, ievc1, sBuf2, sP2, sQ2, ievc2) SPECTRAL_PHASE2_GPU() SPECTRAL_PHASE3_GPU() @@ -744,12 +775,13 @@ KW_GLOBAL_KERNEL void kernelPartialsPartialsGrowingSpectral( KW_GLOBAL_VAR REAL* KW_RESTRICT evec2, KW_GLOBAL_VAR REAL* KW_RESTRICT eigenValues2, KW_GLOBAL_VAR REAL* KW_RESTRICT distances2, + int isAllReal1, int isAllReal2, int totalPatterns) { SPECTRAL_INDICES_GPU() SPECTRAL_COMMON_SMEM_GPU() SPECTRAL_LOAD_PARTIALS1_GPU() /* sP1 = parent pre-order */ SPECTRAL_LOAD_PARTIALS2_GPU() /* sP2 = sibling post-order */ - SPECTRAL_EIGENVALS_GPU() + SPECTRAL_EIGENVALS_GPU(isAllReal1, isAllReal2) /* Sibling: forward (P_sib · p_sib) */ SPECTRAL_PHASE1_PARTIALS_GPU(sBuf2, sP2, sQ2, ievc2) SPECTRAL_PHASE2_SIB_GPU() @@ -772,11 +804,12 @@ KW_GLOBAL_KERNEL void kernelPartialsStatesGrowingSpectral( KW_GLOBAL_VAR REAL* KW_RESTRICT evec2, KW_GLOBAL_VAR REAL* KW_RESTRICT eigenValues2, KW_GLOBAL_VAR REAL* KW_RESTRICT distances2, + int isAllReal1, int isAllReal2, int totalPatterns) { SPECTRAL_INDICES_GPU() SPECTRAL_COMMON_SMEM_GPU() SPECTRAL_LOAD_PARTIALS1_GPU() /* sP1 = parent pre-order */ - SPECTRAL_EIGENVALS_GPU() + SPECTRAL_EIGENVALS_GPU(isAllReal1, isAllReal2) /* Sibling: forward (P_sib · e_s) using States Phase 1 */ SPECTRAL_PHASE1_STATES_GPU(sQ2, ievc2, states2) SPECTRAL_PHASE2_SIB_GPU() @@ -815,11 +848,12 @@ KW_GLOBAL_KERNEL void kernelPartialsStatesGrowingTopSpectral( KW_GLOBAL_VAR REAL* KW_RESTRICT evec2, KW_GLOBAL_VAR REAL* KW_RESTRICT eigenValues2, KW_GLOBAL_VAR REAL* KW_RESTRICT distances2, + int isAllReal1, int isAllReal2, int totalPatterns) { SPECTRAL_INDICES_GPU() SPECTRAL_COMMON_SMEM_GPU() SPECTRAL_LOAD_PARTIALS1_GPU() /* sP1 = parent pre-order */ - SPECTRAL_EIGENVALS_GPU() + SPECTRAL_EIGENVALS_GPU(isAllReal1, isAllReal2) SPECTRAL_PHASE1_PARTIALS_GPU(sBuf1, sP1, sQ1, ievc1) /* parent backward Ph1 */ SPECTRAL_PHASE1_STATES_GPU(sQ2, ievc2, states2) /* sibling forward Ph1 */ SPECTRAL_PHASE2_GPU() /* scale sQ1 and sQ2 independently */ @@ -837,12 +871,13 @@ KW_GLOBAL_KERNEL void kernelPartialsPartialsGrowingTopRootSpectral( KW_GLOBAL_VAR REAL* KW_RESTRICT evec2, KW_GLOBAL_VAR REAL* KW_RESTRICT eigenValues2, KW_GLOBAL_VAR REAL* KW_RESTRICT distances2, + int isAllReal2, int totalPatterns) { SPECTRAL_INDICES_GPU() SPECTRAL_COMMON_SMEM_GPU() SPECTRAL_LOAD_PARTIALS1_GPU() /* sP1 = root pre-order (Hadamard factor) */ SPECTRAL_LOAD_PARTIALS2_GPU() /* sP2 = sibling post-order */ - SPECTRAL_EIGENVALS_SIB_ONLY_GPU() + SPECTRAL_EIGENVALS_SIB_ONLY_GPU(isAllReal2) SPECTRAL_PHASE1_PARTIALS_GPU(sBuf2, sP2, sQ2, ievc2) SPECTRAL_PHASE2_SIB_GPU() SPECTRAL_PHASE3_SIB_HADAMARD_WRITE_GPU() /* evec2·sQ2 ⊙ sP1 → partials3 */ @@ -857,11 +892,12 @@ KW_GLOBAL_KERNEL void kernelPartialsStatesGrowingTopRootSpectral( KW_GLOBAL_VAR REAL* KW_RESTRICT evec2, KW_GLOBAL_VAR REAL* KW_RESTRICT eigenValues2, KW_GLOBAL_VAR REAL* KW_RESTRICT distances2, + int isAllReal2, int totalPatterns) { SPECTRAL_INDICES_GPU() SPECTRAL_COMMON_SMEM_GPU() SPECTRAL_LOAD_PARTIALS1_GPU() /* sP1 = root pre-order (Hadamard factor) */ - SPECTRAL_EIGENVALS_SIB_ONLY_GPU() + SPECTRAL_EIGENVALS_SIB_ONLY_GPU(isAllReal2) SPECTRAL_PHASE1_STATES_GPU(sQ2, ievc2, states2) SPECTRAL_PHASE2_SIB_GPU() SPECTRAL_PHASE3_SIB_HADAMARD_WRITE_GPU() @@ -1229,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