diff --git a/.gitignore b/.gitignore index c00d9ba5..2dc47505 100644 --- a/.gitignore +++ b/.gitignore @@ -209,6 +209,7 @@ examples/tinytest/.libs/ examples/tinytest/tinytest.o libhmsbeagle/GPU/kernels/BeagleOpenCL_kernels.h libhmsbeagle/GPU/kernels/BeagleOpenCL_kernels_xcode.h +libhmsbeagle/GPU/kernels/BeagleOpenCLSpectral_kernels.h libhmsbeagle/GPU/libhmsbeagle-opencl.la libhmsbeagle/GPU/libhmsbeagle_opencl_la-GPUImplHelper.lo libhmsbeagle/GPU/libhmsbeagle_opencl_la-GPUInterfaceOpenCL.lo diff --git a/examples/hmctest/adjointtest4.cpp b/examples/hmctest/adjointtest4.cpp index c28ee171..b3d2260f 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; @@ -202,7 +207,8 @@ static void printResources() { static int createInstance(bool useGpu, int whichDevice, bool singlePrec, int stateCount, int nPatterns, int nCats, - bool complexEigen = false) { + bool complexEigen = false, + int nEigenDecomp = 1, int nTransitionMatrices = 4) { long prefFlags = BEAGLE_FLAG_SCALERS_RAW; prefFlags |= useGpu ? BEAGLE_FLAG_PROCESSOR_GPU : BEAGLE_FLAG_PROCESSOR_CPU; prefFlags |= singlePrec ? BEAGLE_FLAG_PRECISION_SINGLE : BEAGLE_FLAG_PRECISION_DOUBLE; @@ -218,8 +224,8 @@ static int createInstance(bool useGpu, int whichDevice, bool singlePrec, 3, /* nCompact */ stateCount, nPatterns, - 1, /* nEigenDecomp */ - 4, /* nTransitionMatrices */ + nEigenDecomp, + nTransitionMatrices, nCats, 0, /* nScaleBuffers */ whichDevice >= 0 ? &whichDevice : NULL, @@ -934,6 +940,489 @@ 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; +} + +/* ═══════════════════════════════════════════════════════════════════════ + * Multi-eigendecomposition test (S=4, nEigenDecomp=2) + * + * Regression test for the misaligned adjoint-queue/MemsetZero stride bug + * (beagle-bugs/gpu-spectral-misaligned-adjoint-queue-stride/PLAN.md). + * BeagleGPUSpectralImpl allocates dEvecTOrigin/dIevcTOrigin/dGradientOrigin + * with a per-eigendecomposition stride that is *aligned* + * (AlignMemOffset(S*S*sizeof(Real))), which device alignment rules can make + * strictly larger than the logical S*S*sizeof(Real). The queue-building code + * in calculateAdjointCrossProducts (rec[3]/rec[7]) and the MemsetZero call + * that clears the gradient accumulator must use that same aligned stride to + * index into these pooled buffers by eigendecomposition index -- using the + * raw unaligned S*S instead is invisible whenever eigenIndex==0 (multiplied + * away by ei*anything==0), which is true of every other test/benchmark in + * this project. This test creates a real nEigenDecomp=2 instance and runs + * the SAME 4-edge tree/topology twice on the SAME 4 transition-matrix slots: + * once with all edges bound to eigenIndex=0 (JC69 -- expected to pass + * regardless of the bug, since ei=0 can never expose it) and once with all + * edges rebound to eigenIndex=1 (a distinct, non-degenerate real + * eigendecomposition -- exercises the exact rec[3]/rec[7]/MemsetZero code + * path this bug lives in, since ei=1 is where "ei*SS" and "ei*matStride" + * first diverge whenever matStride>SS*sizeof(Real) on this hardware). + * GPU output for each step is compared against an independently-computed + * CPU (double-precision) reference for the same step -- the CPU + * implementation has no analogous aligned/unaligned pooled-buffer stride + * concept, so it is unaffected by this bug and serves as ground truth. + * ═══════════════════════════════════════════════════════════════════════*/ + +struct MultiEigenResult { + double logLA, logLB; + std::vector gradA, gradB; /* S*S each */ +}; + +static MultiEigenResult runMultiEigenOne(bool useGpu, int dev, bool singlePrec, + const std::vector& evecB, + const std::vector& ivecB, + const std::vector& evalB) { + const int S = 4, nPat = 4, nCats = 2; + MultiEigenResult R; + + int inst = createInstance(useGpu, dev, singlePrec, S, nPat, nCats, + /*complexEigen=*/false, /*nEigenDecomp=*/2); + if (inst < 0) { R.logLA = R.logLB = NAN; return R; } + + int *hSt = dnaToStates(human); + int *cSt = dnaToStates(chimp); + int *gSt = dnaToStates(gorilla); + beagleSetTipStates(inst, 0, hSt); + beagleSetTipStates(inst, 1, cSt); + beagleSetTipStates(inst, 2, gSt); + free(hSt); free(cSt); free(gSt); + + beagleSetCategoryRates(inst, rates2); + std::vector pw(nPat, 1.0); + beagleSetPatternWeights(inst, pw.data()); + beagleSetStateFrequencies(inst, 0, freqs4); + beagleSetCategoryWeights(inst, 0, catWeights2); + + int nodeIdx[4] = { 0, 1, 2, 3 }; + BeagleOperation postOps[2] = { + { 3, BEAGLE_OP_NONE, BEAGLE_OP_NONE, 0, 0, 1, 1 }, + { 4, BEAGLE_OP_NONE, BEAGLE_OP_NONE, 2, 2, 3, 3 } + }; + std::vector rootPre((size_t)nCats * nPat * S); + for (int c = 0; c < nCats; c++) + for (int p = 0; p < nPat; p++) + for (int s = 0; s < S; s++) + rootPre[(size_t)c * nPat * S + p * S + s] = freqs4[s]; + BeagleOperation preOps[4] = { + { 6, BEAGLE_OP_NONE, BEAGLE_OP_NONE, 5, BEAGLE_OP_NONE, 2, 2 }, + { 7, BEAGLE_OP_NONE, BEAGLE_OP_NONE, 5, BEAGLE_OP_NONE, 3, 3 }, + { 8, BEAGLE_OP_NONE, BEAGLE_OP_NONE, 6, BEAGLE_OP_NONE, 0, 0 }, + { 9, BEAGLE_OP_NONE, BEAGLE_OP_NONE, 6, BEAGLE_OP_NONE, 1, 1 } + }; + BeagleBranchOperation branchOps[4] = { + { 1, 8, 1, 0 }, + { 0, 9, 0, 0 }, + { 2, 7, 2, 0 }, + { 3, 6, 3, 0 } + }; + + auto runStep = [&](int eigenIndex, const double* evec, const double* ivec, + const double* eval) -> std::pair> { + beagleSetEigenDecomposition(inst, eigenIndex, evec, ivec, eval); + beagleUpdateTransitionMatrices(inst, eigenIndex, nodeIdx, NULL, NULL, edgeLengths, 4); + beagleUpdatePartials(inst, postOps, 2, BEAGLE_OP_NONE); + double logL = 0.0; + int rootIdx = 4, cwIdx = 0, sfIdx = 0, csi = BEAGLE_OP_NONE; + beagleCalculateRootLogLikelihoods(inst, &rootIdx, &cwIdx, &sfIdx, &csi, 1, &logL); + beagleSetPartials(inst, 5, rootPre.data()); + beagleUpdatePrePartials_v5(inst, preOps, 4, BEAGLE_OP_NONE, BEAGLE_PARTIALS_TOP); + std::vector grad(S * S, 0.0); + beagleCalculateAdjointDerivative(inst, branchOps, 0, 0, 4, 0, 4, grad.data(), NULL); + return { logL, grad }; + }; + + /* Step A: eigenIndex=0 (JC69) -- ei==0, bug is provably invisible here. */ + auto a = runStep(0, jcEvec, jcIvec, jcEval); + R.logLA = a.first; R.gradA = a.second; + + /* Step B: eigenIndex=1 (distinct real model) -- exercises the bug. */ + auto b = runStep(1, evecB.data(), ivecB.data(), evalB.data()); + R.logLB = b.first; R.gradB = b.second; + + beagleFinalizeInstance(inst); + return R; +} + +static int runTestMultiEigen(int gpuDevice) { + const int S = 4; + const double tol = 1e-4; + + fprintf(stdout, "\n=== Multi-eigendecomposition test (S=4, nEigenDecomp=2) ===\n"); + setenv("BEAGLE_DEBUG_STRIDE", "1", 1); + + /* Model B: symmetric (r_fwd=r_bkd=1.0) 4-state circulant -- exactly + * real eigenvalues, numerically distinct from JC69. */ + std::vector evecB, ivecB, evalB; + buildCirculantN(4, evecB, ivecB, evalB, /*r_fwd=*/1.0, /*r_bkd=*/1.0); + + fprintf(stdout, "\nCreating CPU (double) instance (nEigenDecomp=2):\n"); + MultiEigenResult cpuR = runMultiEigenOne(false, -1, false, evecB, ivecB, evalB); + if (std::isnan(cpuR.logLA)) { fprintf(stderr, "Failed CPU instance\n"); return 1; } + fprintf(stdout, " CPU step A (eigenIndex=0, JC69): logL=%.6f\n", cpuR.logLA); + fprintf(stdout, " CPU step B (eigenIndex=1, model B): logL=%.6f\n", cpuR.logLB); + + if (gpuDevice < 0) { + fprintf(stdout, "(no --gpu device given; CPU-only run, nothing to compare)\n"); + return 0; + } + + fprintf(stdout, "\nCreating GPU (single) instance (nEigenDecomp=2):\n"); + MultiEigenResult gpuR = runMultiEigenOne(true, gpuDevice, true, evecB, ivecB, evalB); + if (std::isnan(gpuR.logLA)) { fprintf(stderr, "Failed GPU instance\n"); return 1; } + fprintf(stdout, " GPU step A (eigenIndex=0, JC69): logL=%.6f\n", gpuR.logLA); + fprintf(stdout, " GPU step B (eigenIndex=1, model B): logL=%.6f\n", gpuR.logLB); + + auto compareGrad = [&](const char* label, const std::vector& cpu, + const std::vector& gpu) -> bool { + double maxDiff = 0.0; + for (int k = 0; k < S * S; k++) { + double d = fabs(cpu[k] - gpu[k]); + if (d > maxDiff || std::isnan(d)) maxDiff = d; + } + bool ok = (maxDiff <= tol) && !std::isnan(maxDiff); + fprintf(stdout, " %-42s max|CPU-GPU diff|=%.3e (tol=%.0e) %s\n", + label, maxDiff, tol, ok ? "PASS" : "FAIL"); + if (!ok) { + fprintf(stdout, " CPU grad:\n"); + for (int i = 0; i < S; i++) { + fprintf(stdout, " "); + for (int j = 0; j < S; j++) fprintf(stdout, " %10.6f", cpu[i*S+j]); + fprintf(stdout, "\n"); + } + fprintf(stdout, " GPU grad:\n"); + for (int i = 0; i < S; i++) { + fprintf(stdout, " "); + for (int j = 0; j < S; j++) fprintf(stdout, " %10.6f", gpu[i*S+j]); + fprintf(stdout, "\n"); + } + } + return ok; + }; + + bool okA = compareGrad("Step A gradient (eigenIndex=0, JC69)", cpuR.gradA, gpuR.gradA); + bool okB = compareGrad("Step B gradient (eigenIndex=1, model B)", cpuR.gradB, gpuR.gradB); + + fprintf(stdout, "Multi-eigendecomposition test: %s\n", (okA && okB) ? "PASS" : "FAIL"); + return (okA && okB) ? 0 : 1; +} + /* ── main ────────────────────────────────────────────────────────────── */ int main(int argc, const char *argv[]) { @@ -941,9 +1430,12 @@ 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; + bool runMultiEigen = false; for (int i = 1; i < argc; i++) { if (strcmp(argv[i], "--gpu") == 0 && i+1 < argc) { @@ -954,17 +1446,29 @@ 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], "--multieigen") == 0) { + run4 = false; + runMultiEigen = true; } else if (strcmp(argv[i], "--all") == 0) { - run4 = run16 = run17 = true; + run4 = run16 = run17 = runRealOnly = runDynamic = runMultiEigen = 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); + if (runMultiEigen) result |= runTestMultiEigen(gpuDev); return result; } diff --git a/libhmsbeagle/GPU/AdjointBlockSize.h b/libhmsbeagle/GPU/AdjointBlockSize.h new file mode 100644 index 00000000..82b9637b --- /dev/null +++ b/libhmsbeagle/GPU/AdjointBlockSize.h @@ -0,0 +1,7 @@ +#ifndef BEAGLE_GPU_ADJOINT_BLOCK_SIZE_H +#define BEAGLE_GPU_ADJOINT_BLOCK_SIZE_H + +/* Shared by the host launcher and the generated OpenCL generic-N source. */ +#define ADJOINT_BLOCK_SP_N 64 + +#endif diff --git a/libhmsbeagle/GPU/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..57e4f5e5 100644 --- a/libhmsbeagle/GPU/BeagleGPUSpectralImpl.h +++ b/libhmsbeagle/GPU/BeagleGPUSpectralImpl.h @@ -38,6 +38,19 @@ class BeagleGPUSpectralImpl : public BeagleGPUImpl { GPUPtr dSpectralDistancesOrigin; GPUPtr* dSpectralDistances; int* hEigenIndexForMatrix; + + /* Lazy dense-matrix materialization (see PLAN.md in + * beagle-bugs/gpu-spectral-skip-redundant-dense-matrix/): probability-only + * updateTransitionMatrices calls no longer eagerly build the dense + * kStateCount x kStateCount transition matrix, since the spectral + * pruning/gradient paths never read it. Instead the eigenIndex/edgeLength + * needed to rebuild it are cached here, and the dense matrix is + * materialized on first actual getTransitionMatrix() call. All sized + * kMatrixCount. */ + bool* hDenseMatrixValid; + int* hPendingDenseEigenIndex; + double* hPendingDenseEdgeLength; + /* Per-matrix element stride within dSpectralDistancesOrigin, i.e. * AlignMemOffset(kCategoryCount * sizeof(Real)) / sizeof(Real). Device * alignment padding can make this larger than kCategoryCount, so any code @@ -45,6 +58,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). @@ -62,6 +91,17 @@ class BeagleGPUSpectralImpl : public BeagleGPUImpl { GPUPtr dGradientOrigin; GPUPtr* dGradient; bool* hEigenDecompIsAllReal; + int* hEigenDecompComplexPairLeaders; + int* hEigenDecompSingletonRows; + + /* Per-eigendecomposition element/byte stride shared by dEvecTOrigin, + * dIevcTOrigin, and dGradientOrigin, i.e. AlignMemOffset(S*S*sizeof(Real)) + * (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 { const double* edgeLengths, int count) override; + int getTransitionMatrix(int matrixIndex, + double* outMatrix) override; + int updatePrePartials(const int* operations, int operationCount, int cumulativeScaleIndex, diff --git a/libhmsbeagle/GPU/BeagleGPUSpectralImpl.hpp b/libhmsbeagle/GPU/BeagleGPUSpectralImpl.hpp index 10e7a418..e95b8533 100644 --- a/libhmsbeagle/GPU/BeagleGPUSpectralImpl.hpp +++ b/libhmsbeagle/GPU/BeagleGPUSpectralImpl.hpp @@ -18,6 +18,59 @@ #define __BeagleGPUSpectralImpl_hpp__ #include "libhmsbeagle/beagle.h" +#include +#include +#include + +namespace { +struct BeagleAdjointDetailStats { + long calls = 0; + double total = 0.0, integrate = 0.0, memset = 0.0, queueBuild = 0.0; + double queueH2D = 0.0, mergedLaunch = 0.0, gradientD2H = 0.0, hostCopy = 0.0; + ~BeagleAdjointDetailStats() { + if (!calls) return; + fprintf(stderr, "\n===== BEAGLE_PROFILE_ADJOINT_DETAIL =====\n"); +#define BEAGLE_ADJOINT_LINE(label, value) fprintf(stderr, "%-28s: calls=%-8ld total=%.6fs avg=%.3fus\n", label, calls, value, value * 1e6 / calls) + BEAGLE_ADJOINT_LINE("calculateAdjoint total", total); + BEAGLE_ADJOINT_LINE("IntegrateLikelihoods enqueue", integrate); + BEAGLE_ADJOINT_LINE("MemsetZero host wall", memset); + BEAGLE_ADJOINT_LINE("hAdjointQueue construction", queueBuild); + BEAGLE_ADJOINT_LINE("adjoint queue H2D", queueH2D); + BEAGLE_ADJOINT_LINE("merged adjoint enqueue", mergedLaunch); + BEAGLE_ADJOINT_LINE("gradient D2H", gradientD2H); + BEAGLE_ADJOINT_LINE("padded gradient host copy", hostCopy); +#undef BEAGLE_ADJOINT_LINE + fprintf(stderr, "===== END BEAGLE_PROFILE_ADJOINT_DETAIL =====\n\n"); + } +}; +struct BeagleAdjointCompositionStats { + long calls = 0, branches = 0, statesBranches = 0, allRealBranches = 0; + long complexPairLeaders = 0, singletonRows = 0; + int paddedPatterns = 0, paddedStates = 0, categories = 0; + ~BeagleAdjointCompositionStats() { + if (!calls) return; + fprintf(stderr, "\n===== BEAGLE_PROFILE_ADJOINT_COMPOSITION =====\n"); + fprintf(stderr, "calls=%ld branches=%ld isStates=%ld isAllReal=%ld " + "complexPairLeaders=%ld singletonRows=%ld " + "kPaddedPatternCount=%d kPaddedStateCount=%d kCategoryCount=%d\n", + calls, branches, statesBranches, allRealBranches, + complexPairLeaders, singletonRows, + paddedPatterns, paddedStates, categories); + fprintf(stderr, "===== END BEAGLE_PROFILE_ADJOINT_COMPOSITION =====\n\n"); + } +}; +inline BeagleAdjointCompositionStats& beagleAdjointCompositionStats() { + static BeagleAdjointCompositionStats stats; + return stats; +} +inline BeagleAdjointDetailStats& beagleAdjointDetailStats() { + static BeagleAdjointDetailStats stats; + return stats; +} +inline double beagleAdjointElapsed(std::chrono::high_resolution_clock::time_point t) { + return std::chrono::duration(std::chrono::high_resolution_clock::now() - t).count(); +} +} namespace beagle { namespace gpu { @@ -34,8 +87,12 @@ namespace gpu { BEAGLE_GPU_TEMPLATE BeagleGPUSpectralImpl::BeagleGPUSpectralImpl() : dSpectralDistancesOrigin(0), dSpectralDistances(NULL), hEigenIndexForMatrix(NULL), + hDenseMatrixValid(NULL), hPendingDenseEigenIndex(NULL), hPendingDenseEdgeLength(NULL), + dSpectralPtrQueue(0), dSpectralDistanceQueue(0), + hSpectralPtrQueue(NULL), hSpectralDistanceQueue(NULL), dEvecTOrigin(0), dEvecT(NULL), dIevcTOrigin(0), dIevcT(NULL), dGradientOrigin(0), dGradient(NULL), hEigenDecompIsAllReal(NULL), + hEigenDecompComplexPairLeaders(NULL), hEigenDecompSingletonRows(NULL), dAdjointQueue(0), hAdjointQueue(NULL), kAdjointQueueCapacity(0), kSpectralEigenDecompCount(0) { } @@ -45,6 +102,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); @@ -56,7 +117,12 @@ BeagleGPUSpectralImpl::~BeagleGPUSpectralImpl() { free(dIevcT); free(dGradient); free(hEigenIndexForMatrix); + free(hDenseMatrixValid); + free(hPendingDenseEigenIndex); + free(hPendingDenseEdgeLength); free(hEigenDecompIsAllReal); + free(hEigenDecompComplexPairLeaders); + free(hEigenDecompSingletonRows); } BEAGLE_GPU_TEMPLATE @@ -84,6 +150,12 @@ int BeagleGPUSpectralImpl::createInstance( kSpectralEigenDecompCount = eigenDecompositionCount; hEigenIndexForMatrix = (int*) calloc(this->kMatrixCount, sizeof(int)); + hEigenDecompComplexPairLeaders = (int*) calloc(eigenDecompositionCount, sizeof(int)); + hEigenDecompSingletonRows = (int*) calloc(eigenDecompositionCount, sizeof(int)); + + hDenseMatrixValid = (bool*) calloc(this->kMatrixCount, sizeof(bool)); // false = stale/never built + hPendingDenseEigenIndex = (int*) calloc(this->kMatrixCount, sizeof(int)); + hPendingDenseEdgeLength = (double*) calloc(this->kMatrixCount, sizeof(double)); GPUInterface* gpuIf = this->gpu; @@ -95,9 +167,27 @@ 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)); + kSpectralMatrixStrideBytes = matStride; + kSpectralMatrixStrideElements = (unsigned int)(matStride / sizeof(Real)); + if (getenv("BEAGLE_DEBUG_STRIDE")) { + fprintf(stderr, "[BEAGLE_DEBUG_STRIDE] S=%d matStride=%zu bytes, S*S*sizeof(Real)=%zu bytes, " + "aligned>logical=%s\n", S, matStride, (size_t)S * S * sizeof(Real), + (matStride > (size_t)S * S * sizeof(Real)) ? "YES" : "no"); + } dEvecT = (GPUPtr*) malloc(sizeof(GPUPtr) * eigenDecompositionCount); dIevcT = (GPUPtr*) malloc(sizeof(GPUPtr) * eigenDecompositionCount); dEvecTOrigin = gpuIf->AllocateMemory(eigenDecompositionCount * matStride); @@ -127,30 +217,65 @@ int BeagleGPUSpectralImpl::updateTransitionMatrices( const int* secondDerivativeIndices, const double* edgeLengths, int count) { - int rc = BeagleGPUImpl::updateTransitionMatrices(eigenIndex, probabilityIndices, - firstDerivativeIndices, - secondDerivativeIndices, - edgeLengths, count); - if (rc != BEAGLE_SUCCESS || firstDerivativeIndices != NULL) return rc; + // Only defer the dense (probability-only) case. Derivative requests are rare in this + // project's usage and the existing eager-build behavior for them is left completely + // untouched, to keep this change minimal and scoped to the one case that's actually + // wasteful (see PLAN.md in beagle-bugs/gpu-spectral-skip-redundant-dense-matrix/). + if (firstDerivativeIndices != NULL) { + return BeagleGPUImpl::updateTransitionMatrices(eigenIndex, probabilityIndices, + firstDerivativeIndices, + secondDerivativeIndices, + edgeLengths, count); + } + + // Probability-only update: DO NOT build the dense matrix now. Cache what's needed to + // build it lazily (see getTransitionMatrix override below), in case some caller asks + // for it later. Spectral pruning/gradient paths never read dMatrices. + for (int i = 0; i < count; i++) { + int matIdx = probabilityIndices[i]; + hPendingDenseEigenIndex[matIdx] = eigenIndex; + hPendingDenseEdgeLength[matIdx] = edgeLengths[i]; + hDenseMatrixValid[matIdx] = false; + } if (count > 0) { const double* categoryRates = this->hCategoryRates[0]; 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; } +BEAGLE_GPU_TEMPLATE +int BeagleGPUSpectralImpl::getTransitionMatrix(int matrixIndex, + double* outMatrix) { + if (!hDenseMatrixValid[matrixIndex]) { + // Materialize on demand, reusing the existing, already-correct base-class + // construction unchanged -- just deferred to first actual use instead of built + // eagerly and thrown away. + int rc = BeagleGPUImpl::updateTransitionMatrices( + hPendingDenseEigenIndex[matrixIndex], &matrixIndex, + NULL, NULL, &hPendingDenseEdgeLength[matrixIndex], 1); + if (rc != BEAGLE_SUCCESS) return rc; + hDenseMatrixValid[matrixIndex] = true; + } + return BeagleGPUImpl::getTransitionMatrix(matrixIndex, outMatrix); +} + BEAGLE_GPU_TEMPLATE void BeagleGPUSpectralImpl::dispatchPrunePP( GPUPtr p1, GPUPtr p2, GPUPtr p3, @@ -165,6 +290,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 +308,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 +319,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 +326,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 +344,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] @@ -281,6 +408,11 @@ int BeagleGPUSpectralImpl::setEigenDecomposition( for (int i = 0; i < SC && allReal; i++) if (inEigenValues[SC + i] != 0.0) allReal = false; hEigenDecompIsAllReal[eigenIndex] = allReal; + int pairLeaders = 0; + for (int i = 0; i < SC; i++) + if (inEigenValues[SC + i] > 0.0) pairLeaders++; + hEigenDecompComplexPairLeaders[eigenIndex] = pairLeaders; + hEigenDecompSingletonRows[eigenIndex] = S - 2 * pairLeaders; if (this->kFlags & BEAGLE_FLAG_EIGEN_COMPLEX) { for (int i = 0; i < SC; i++) @@ -320,18 +452,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 +476,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 +484,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 +494,7 @@ void BeagleGPUSpectralImpl::dispatchGrowingSpectral( ievc2, evec2, eval2, dist2, nullptr, nullptr, this->kPaddedPatternCount, this->kCategoryCount, + isAllReal1, isAllReal2, -1, -1, -1); } } else { @@ -365,13 +503,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); } } } @@ -546,12 +686,18 @@ int BeagleGPUSpectralImpl::calculateAdjointCrossProducts( double* outSumDerivatives, double* outSumSquaredDerivatives) { + const bool detail = (getenv("BEAGLE_PROFILE_ADJOINT_DETAIL") != NULL); + const bool composition = (getenv("BEAGLE_PROFILE_ADJOINT_COMPOSITION") != NULL); + BeagleAdjointDetailStats& detailStats = beagleAdjointDetailStats(); + std::chrono::high_resolution_clock::time_point stageStart, callStart; + if (detail) callStart = std::chrono::high_resolution_clock::now(); GPUInterface* gpuIf = this->gpu; const int S = this->kPaddedStateCount; const int SS = S * S; const int cwIdx = categoryWeightsIndices[0]; /* Compute per-site marginal likelihoods into dIntegrationTmp. */ + if (detail) stageStart = std::chrono::high_resolution_clock::now(); this->kernels->IntegrateLikelihoods( this->dIntegrationTmp, this->dPartials[rootPostOrderIndex], @@ -559,6 +705,7 @@ int BeagleGPUSpectralImpl::calculateAdjointCrossProducts( this->dFrequencies[stateFrequenciesIndex], this->kPaddedPatternCount, this->kCategoryCount); + if (detail) detailStats.integrate += beagleAdjointElapsed(stageStart); /* Zero the gradient accumulator on-device, in one fill over the whole * pooled origin buffer (covers every eigen decomposition, not just @@ -572,11 +719,14 @@ int BeagleGPUSpectralImpl::calculateAdjointCrossProducts( * through the parent view (intermittent large garbage gradients, * OpenCL, Apple M1 Max) even on the in-order command queue — filling * through the same view the kernel uses avoids that hazard. */ - gpuIf->MemsetZero(dGradientOrigin, (size_t)kSpectralEigenDecompCount * SS * sizeof(Real)); + if (detail) stageStart = std::chrono::high_resolution_clock::now(); + gpuIf->MemsetZero(dGradientOrigin, (size_t)kSpectralEigenDecompCount * kSpectralMatrixStrideBytes); + if (detail) detailStats.memset += beagleAdjointElapsed(stageStart); GPUPtr catW = this->dWeights[cwIdx]; if (count > 0) { + if (detail) stageStart = std::chrono::high_resolution_clock::now(); /* Merged single-launch path (both S=4 and generic-N): build one * offset-queue record per branch (single buffer, integer offsets * into pooled origins — never per-branch device pointers), then one @@ -604,23 +754,35 @@ int BeagleGPUSpectralImpl::calculateAdjointCrossProducts( const int ei = hEigenIndexForMatrix[matIdx]; const bool isStates = (this->dStates[postIdx] != 0); const bool isAllReal = hEigenDecompIsAllReal[ei]; + if (composition) { + BeagleAdjointCompositionStats& s = beagleAdjointCompositionStats(); + s.branches++; + if (isStates) s.statesBranches++; + if (isAllReal) s.allRealBranches++; + s.complexPairLeaders += hEigenDecompComplexPairLeaders[ei]; + s.singletonRows += hEigenDecompSingletonRows[ei]; + } unsigned int* rec = hAdjointQueue + (size_t)i * kAdjointQueueFieldsPerBranch; rec[0] = this->getPartialsOffsetElements(preIdx); rec[1] = isStates ? this->getStatesOffsetElements(postIdx) : this->getPartialsOffsetElements(postIdx); rec[2] = isStates ? 1u : 0u; - rec[3] = (unsigned int)ei * SS; + rec[3] = (unsigned int)ei * kSpectralMatrixStrideElements; rec[4] = (unsigned int)ei * this->getEvecStrideElements(); rec[5] = (unsigned int)ei * this->getEvalStrideElements(); rec[6] = (unsigned int)matIdx * kSpectralDistanceStrideElements; - rec[7] = (unsigned int)ei * SS; + rec[7] = (unsigned int)ei * kSpectralMatrixStrideElements; rec[8] = isAllReal ? 1u : 0u; } + if (detail) detailStats.queueBuild += beagleAdjointElapsed(stageStart); + if (detail) stageStart = std::chrono::high_resolution_clock::now(); gpuIf->MemcpyHostToDevice(dAdjointQueue, hAdjointQueue, sizeof(unsigned int) * count * kAdjointQueueFieldsPerBranch); + if (detail) detailStats.queueH2D += beagleAdjointElapsed(stageStart); + if (detail) stageStart = std::chrono::high_resolution_clock::now(); if (S == 4) { this->kernels->AdjointCrossProductMergedSpectral4( this->getPartialsOrigin(), this->getStatesOrigin(), dEvecTOrigin, @@ -636,6 +798,7 @@ int BeagleGPUSpectralImpl::calculateAdjointCrossProducts( this->dIntegrationTmp, dGradientOrigin, dAdjointQueue, this->kPaddedPatternCount, this->kCategoryCount, count); } + if (detail) detailStats.mergedLaunch += beagleAdjointElapsed(stageStart); } /* Download gradient to host. The first eigen decomp's gradient is the @@ -653,13 +816,29 @@ int BeagleGPUSpectralImpl::calculateAdjointCrossProducts( const int ei0 = hEigenIndexForMatrix[eigenIndices[0]]; const int SC = this->kStateCount; Real* hGrad = (Real*) gpuIf->CallocHost(sizeof(Real), SS); + if (detail) stageStart = std::chrono::high_resolution_clock::now(); gpuIf->MemcpyDeviceToHost(hGrad, dGradient[ei0], SS * sizeof(Real)); + if (detail) detailStats.gradientD2H += beagleAdjointElapsed(stageStart); + if (detail) stageStart = std::chrono::high_resolution_clock::now(); for (int i = 0; i < SC; i++) for (int j = 0; j < SC; j++) outSumDerivatives[i * SC + j] = (double)hGrad[i * S + j]; + if (detail) detailStats.hostCopy += beagleAdjointElapsed(stageStart); gpuIf->FreeHostMemory(hGrad); } + if (detail) { + detailStats.calls++; + detailStats.total += beagleAdjointElapsed(callStart); + } + if (composition) { + BeagleAdjointCompositionStats& s = beagleAdjointCompositionStats(); + s.calls++; + s.paddedPatterns = this->kPaddedPatternCount; + s.paddedStates = this->kPaddedStateCount; + s.categories = this->kCategoryCount; + } + return BEAGLE_SUCCESS; } diff --git a/libhmsbeagle/GPU/CMake_OpenCLSpectral/CMakeLists.txt b/libhmsbeagle/GPU/CMake_OpenCLSpectral/CMakeLists.txt index 010930fc..57f8b5ca 100644 --- a/libhmsbeagle/GPU/CMake_OpenCLSpectral/CMakeLists.txt +++ b/libhmsbeagle/GPU/CMake_OpenCLSpectral/CMakeLists.txt @@ -22,6 +22,7 @@ add_custom_command( "../kernels/kernelsXDerivatives.cu" "../kernels/kernelsAll.cu" "../kernels/kernelsSpectralIfDef.cu" + "../AdjointBlockSize.h" "../kernels/kernelsSpectralIfDef4.cu" COMMENT "Building OpenCL spectral kernels" VERBATIM) diff --git a/libhmsbeagle/GPU/GPUInterfaceOpenCL.cpp b/libhmsbeagle/GPU/GPUInterfaceOpenCL.cpp index f26f11e0..b5be7667 100644 --- a/libhmsbeagle/GPU/GPUInterfaceOpenCL.cpp +++ b/libhmsbeagle/GPU/GPUInterfaceOpenCL.cpp @@ -22,6 +22,9 @@ #include #include #include +#include +#include +#include #include "libhmsbeagle/beagle.h" #include "libhmsbeagle/GPU/GPUImplDefs.h" @@ -29,6 +32,201 @@ #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; +}; + +struct BeagleOpenCLEventStat { + long calls = 0; + double seconds = 0.0; +}; + +struct BeaglePendingEvent { + cl_event event; + std::string name; +}; + +BeagleProfileStats gBeagleProfile; +bool gBeagleOpenCLEventsInitialized = false; +bool gBeagleOpenCLEventsEnabled = false; +std::vector gBeaglePendingEvents; +std::map gBeagleOpenCLEventStats; + +void beagleOpenCLEventPrintSummary() { + if (!gBeagleOpenCLEventsEnabled) return; + fprintf(stderr, "\n===== BEAGLE_PROFILE_OPENCL_EVENTS =====\n"); + for (std::map::const_iterator it = + gBeagleOpenCLEventStats.begin(); it != gBeagleOpenCLEventStats.end(); ++it) { + fprintf(stderr, "%-28s: calls=%-8ld device=%.6fs avg=%.3fus\n", + it->first.c_str(), it->second.calls, it->second.seconds, + it->second.calls ? it->second.seconds * 1e6 / it->second.calls : 0.0); + } + fprintf(stderr, "pending events not yet drained: %zu\n", gBeaglePendingEvents.size()); + fprintf(stderr, "===== END BEAGLE_PROFILE_OPENCL_EVENTS =====\n\n"); +} + +inline bool beagleOpenCLEventsEnabled() { + if (!gBeagleOpenCLEventsInitialized) { + gBeagleOpenCLEventsEnabled = (getenv("BEAGLE_PROFILE_OPENCL_EVENTS") != NULL); + gBeagleOpenCLEventsInitialized = true; + if (gBeagleOpenCLEventsEnabled) atexit(beagleOpenCLEventPrintSummary); + } + return gBeagleOpenCLEventsEnabled; +} + +void beagleRecordCompletedEvent(const char* name, cl_event event) { + cl_ulong start = 0, end = 0; + int error = clGetEventProfilingInfo(event, CL_PROFILING_COMMAND_START, sizeof(start), &start, NULL); + if (error != CL_SUCCESS) { fprintf(stderr, "OpenCL event profiling start query failed: %d\n", error); exit(-1); } + error = clGetEventProfilingInfo(event, CL_PROFILING_COMMAND_END, sizeof(end), &end, NULL); + if (error != CL_SUCCESS) { fprintf(stderr, "OpenCL event profiling end query failed: %d\n", error); exit(-1); } + BeagleOpenCLEventStat& s = gBeagleOpenCLEventStats[name]; + s.calls++; + s.seconds += (double)(end - start) * 1e-9; + error = clReleaseEvent(event); + if (error != CL_SUCCESS) { fprintf(stderr, "OpenCL event release failed: %d\n", error); exit(-1); } +} + +void beagleDrainPendingEvents() { + for (size_t i = 0; i < gBeaglePendingEvents.size(); ++i) + beagleRecordCompletedEvent(gBeaglePendingEvents[i].name.c_str(), + gBeaglePendingEvents[i].event); + gBeaglePendingEvents.clear(); +} + +bool beagleTargetKernel(GPUFunction kernel, char* name, size_t size) { + if (!beagleOpenCLEventsEnabled()) return false; + int error = clGetKernelInfo(kernel, CL_KERNEL_FUNCTION_NAME, size, name, NULL); + if (error != CL_SUCCESS) { fprintf(stderr, "OpenCL kernel-name query failed: %d\n", error); exit(-1); } + return strcmp(name, "kernelAdjointMergedN") == 0 || + strcmp(name, "kernelAdjointMerged4") == 0 || + strcmp(name, "kernelIntegrateLikelihoods") == 0; +} + +void beagleReportAdjointKernelInfoOnce(GPUFunction kernel, cl_device_id device, + size_t requestedLocalSize) { + static bool reported = false; + if (reported || getenv("BEAGLE_PROFILE_ADJOINT_KERNEL_INFO") == NULL) return; + char name[256] = {0}; + int error = clGetKernelInfo(kernel, CL_KERNEL_FUNCTION_NAME, sizeof(name), name, NULL); + if (error != CL_SUCCESS) { fprintf(stderr, "OpenCL adjoint kernel-name query failed: %d\n", error); exit(-1); } + if (strcmp(name, "kernelAdjointMergedN") != 0) return; + size_t maxWorkGroup = 0, preferredMultiple = 0; + cl_ulong localMemory = 0, privateMemory = 0; + error = clGetKernelWorkGroupInfo(kernel, device, CL_KERNEL_WORK_GROUP_SIZE, + sizeof(maxWorkGroup), &maxWorkGroup, NULL); + if (error != CL_SUCCESS) { fprintf(stderr, "OpenCL work-group-size query failed: %d\n", error); exit(-1); } + error = clGetKernelWorkGroupInfo(kernel, device, + CL_KERNEL_PREFERRED_WORK_GROUP_SIZE_MULTIPLE, + sizeof(preferredMultiple), &preferredMultiple, NULL); + if (error != CL_SUCCESS) { fprintf(stderr, "OpenCL preferred-multiple query failed: %d\n", error); exit(-1); } + error = clGetKernelWorkGroupInfo(kernel, device, CL_KERNEL_LOCAL_MEM_SIZE, + sizeof(localMemory), &localMemory, NULL); + if (error != CL_SUCCESS) { fprintf(stderr, "OpenCL local-memory query failed: %d\n", error); exit(-1); } + error = clGetKernelWorkGroupInfo(kernel, device, CL_KERNEL_PRIVATE_MEM_SIZE, + sizeof(privateMemory), &privateMemory, NULL); + if (error != CL_SUCCESS) { fprintf(stderr, "OpenCL private-memory query failed: %d\n", error); exit(-1); } + fprintf(stderr, "\n===== BEAGLE_PROFILE_ADJOINT_KERNEL_INFO =====\n"); + fprintf(stderr, "kernel=%s requestedLocalSize=%zu workGroupSize=%zu " + "preferredMultiple=%zu localMemory=%llu privateMemory=%llu " + "valid=%s\n", + name, requestedLocalSize, maxWorkGroup, preferredMultiple, + (unsigned long long)localMemory, (unsigned long long)privateMemory, + requestedLocalSize <= maxWorkGroup ? "YES" : "NO"); + fprintf(stderr, "===== END BEAGLE_PROFILE_ADJOINT_KERNEL_INFO =====\n\n"); + if (requestedLocalSize > maxWorkGroup) { + fprintf(stderr, "kernelAdjointMergedN local size exceeds CL_KERNEL_WORK_GROUP_SIZE\n"); + exit(-1); + } + reported = true; +} + +void beagleProfilePrintSummary() { + if (!gBeagleProfile.enabled) return; + 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: %s\n", getenv("BEAGLE_PROFILE_NO_WARMUP") + ? "disabled by BEAGLE_PROFILE_NO_WARMUP" : "first 300 SynchronizeHost calls"); + fprintf(stderr, "LaunchKernel : calls=%-8ld total=%.6fs avg=%.3fus\n", + p.launchKernelCalls, p.launchKernelSeconds, + p.launchKernelCalls ? p.launchKernelSeconds * 1e6 / p.launchKernelCalls : 0.0); + 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); + if (gBeagleProfile.enabled && getenv("BEAGLE_PROFILE_NO_WARMUP") != NULL) { + gBeagleProfile.warm = true; + gBeagleProfile.warmupRemaining = 0; + } + gBeagleProfile.registered = true; + if (gBeagleProfile.enabled) { + atexit(beagleProfilePrintSummary); + } + } + return gBeagleProfile.enabled; +} + +} // anonymous namespace + #define SAFE_CL(call) { \ int error = call; \ if(error != CL_SUCCESS) { \ @@ -325,7 +523,8 @@ void GPUInterface::SetDevice(int deviceNumber, openClCommandQueues = (cl_command_queue*) malloc(sizeof(cl_command_queue) * BEAGLE_STREAM_COUNT); openClEvents = (cl_event*) malloc(sizeof(cl_event) * BEAGLE_STREAM_COUNT); - cl_command_queue_properties queueProperties = 0;//CL_QUEUE_OUT_OF_ORDER_EXEC_MODE_ENABLE; + cl_command_queue_properties queueProperties = beagleOpenCLEventsEnabled() + ? CL_QUEUE_PROFILING_ENABLE : 0;//CL_QUEUE_OUT_OF_ORDER_EXEC_MODE_ENABLE; for (int i=0; i < BEAGLE_STREAM_COUNT; i++) { openClCommandQueues[i] = clCreateCommandQueue(openClContext, openClDeviceId, queueProperties, &err); @@ -477,8 +676,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 +760,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 +786,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 +806,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 +824,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,23 +836,37 @@ 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); } + beagleReportAdjointKernelInfoOnce(deviceFunction, openClDeviceId, localWorkSize[0]); + char eventKernelName[256] = {0}; + const bool targetEvent = beagleTargetKernel(deviceFunction, eventKernelName, + sizeof(eventKernelName)); + cl_event kernelEvent = NULL; + cl_event* kernelEventPtr = targetEvent ? &kernelEvent : NULL; if (globalWorkSize[1] == 1 && globalWorkSize[2] == 1) { SAFE_CL(clEnqueueNDRangeKernel(openClCommandQueues[0], deviceFunction, 1, NULL, - globalWorkSize, localWorkSize, 0, NULL, NULL)); + globalWorkSize, localWorkSize, 0, NULL, kernelEventPtr)); } else if (globalWorkSize[2] == 1) { SAFE_CL(clEnqueueNDRangeKernel(openClCommandQueues[0], deviceFunction, 2, NULL, - globalWorkSize, localWorkSize, 0, NULL, NULL)); + globalWorkSize, localWorkSize, 0, NULL, kernelEventPtr)); } else { SAFE_CL(clEnqueueNDRangeKernel(openClCommandQueues[0], deviceFunction, 3, NULL, - globalWorkSize, localWorkSize, 0, NULL, NULL)); + globalWorkSize, localWorkSize, 0, NULL, kernelEventPtr)); + } + if (targetEvent) gBeaglePendingEvents.push_back(BeaglePendingEvent{kernelEvent, eventKernelName}); + + if (prof && gBeagleProfile.warm) { + double dt = std::chrono::duration( + std::chrono::high_resolution_clock::now() - t0).count(); + gBeagleProfile.launchKernelCalls++; + gBeagleProfile.launchKernelSeconds += dt; } - if (getenv("BEAGLE_DEBUG_KERNEL_ARGS")) { + if (debugArgs) { fprintf(stderr, "[LaunchKernel] %s: enqueued OK, calling clFinish\n", kernelNameBuf); fflush(stderr); clFinish(openClCommandQueues[0]); @@ -656,17 +892,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 +918,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 +938,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 +963,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 +1007,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 +1256,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 @@ -1017,10 +1271,23 @@ void GPUInterface::MemsetZero(GPUPtr dest, * clEnqueueWriteBuffer path MemcpyHostToDevice already uses elsewhere * in this file, which has no such issue. */ void* hZero = calloc(1, memSize); + cl_event event = NULL; SAFE_CL(clEnqueueWriteBuffer(openClCommandQueues[0], dest, CL_TRUE, 0, memSize, hZero, 0, - NULL, NULL)); + NULL, beagleOpenCLEventsEnabled() ? &event : NULL)); + if (event) { + beagleDrainPendingEvents(); + beagleRecordCompletedEvent("MemsetZero", event); + } free(hZero); + if (prof && gBeagleProfile.warm) { + 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,8 +1301,27 @@ 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. + cl_event event = NULL; SAFE_CL(clEnqueueWriteBuffer(openClCommandQueues[0], dest, CL_TRUE, 0, memSize, src, 0, - NULL, NULL)); + NULL, beagleOpenCLEventsEnabled() ? &event : NULL)); + if (event) { + beagleDrainPendingEvents(); + beagleRecordCompletedEvent("MemcpyHostToDevice", event); + } + + if (prof && gBeagleProfile.warm) { + double dt = std::chrono::duration( + 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"); @@ -1049,8 +1335,27 @@ 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. + cl_event event = NULL; SAFE_CL(clEnqueueReadBuffer(openClCommandQueues[0], src, CL_TRUE, 0, memSize, dest, 0, - NULL, NULL)); + NULL, beagleOpenCLEventsEnabled() ? &event : NULL)); + if (event) { + beagleDrainPendingEvents(); + beagleRecordCompletedEvent("MemcpyDeviceToHost", event); + } + + if (prof && gBeagleProfile.warm) { + double dt = std::chrono::duration( + 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"); @@ -1440,4 +1745,3 @@ void GPUInterface::ReleaseDevice(int deviceNumber) { }; // namespace - diff --git a/libhmsbeagle/GPU/KernelLauncher.cpp b/libhmsbeagle/GPU/KernelLauncher.cpp index f9de48c0..73463387 100644 --- a/libhmsbeagle/GPU/KernelLauncher.cpp +++ b/libhmsbeagle/GPU/KernelLauncher.cpp @@ -24,6 +24,7 @@ #include "libhmsbeagle/beagle.h" #include "libhmsbeagle/GPU/GPUImplDefs.h" #include "libhmsbeagle/GPU/KernelLauncher.h" +#include "libhmsbeagle/GPU/AdjointBlockSize.h" /**************CODE***********/ @@ -318,6 +319,8 @@ void KernelLauncher::LoadKernels() { } else { fAdjointMergedN = gpu->GetFunction("kernelAdjointMergedN"); } + + fScatterSpectralDistances = gpu->GetFunction("kernelScatterSpectralDistances"); } fPartialsPartialsEdgeFirstDerivatives = gpu->GetFunction( @@ -1041,13 +1044,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 +1066,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 +1088,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 +1108,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 +1126,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(); } @@ -1161,7 +1177,7 @@ void KernelLauncher::AdjointCrossProductMergedSpectralN(GPUPtr partialsOrigin, unsigned int patternCount, unsigned int categoryCount, int branchCount) { - Dim3Int block(128); + Dim3Int block(ADJOINT_BLOCK_SP_N); Dim3Int grid(kPaddedStateCount, categoryCount, branchCount); gpu->LaunchKernel(fAdjointMergedN, block, grid, 11, 12, partialsOrigin, statesOrigin, evecTOrigin, ievcOrigin, evalOrigin, @@ -1171,6 +1187,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 +1816,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 +1849,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 +1870,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 +1879,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 +1893,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 +1914,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 +1923,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 +1937,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); } } @@ -2711,4 +2761,3 @@ void KernelLauncher::SumSites3(GPUPtr dArray1, }; // namespace - 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..82885eec 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() @@ -857,7 +893,9 @@ KW_GLOBAL_KERNEL void kernelPartialsStatesGrowingTopRootSpectral( #pragma OPENCL EXTENSION cl_khr_int64_base_atomics : enable #endif -#define ADJOINT_BLOCK_SP_N 128 +#ifndef ADJOINT_BLOCK_SP_N +#error "ADJOINT_BLOCK_SP_N must be supplied by AdjointBlockSize.h" +#endif /* CAS-retry float/double atomic add, in its own KW_DEVICE_FUNC-free OpenCL * helper function rather than inlined as a macro at the call site. @@ -1290,6 +1328,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 diff --git a/libhmsbeagle/GPU/kernels/make_opencl_spectral_kernels.sh b/libhmsbeagle/GPU/kernels/make_opencl_spectral_kernels.sh index b5bc5ed1..19ef629a 100755 --- a/libhmsbeagle/GPU/kernels/make_opencl_spectral_kernels.sh +++ b/libhmsbeagle/GPU/kernels/make_opencl_spectral_kernels.sh @@ -37,6 +37,7 @@ for s in $STATE_COUNT_LIST; do cat $srcdir/kernelsAll.cu | sed 's/\\/\\\\/g' | sed 's/\"/\\"/g' | sed 's/$/\\n\\/' >> BeagleOpenCLSpectral_kernels.h cat $srcdir/kernelsX.cu | sed 's/\\/\\\\/g' | sed 's/\"/\\"/g' | sed 's/$/\\n\\/' >> BeagleOpenCLSpectral_kernels.h cat $srcdir/kernelsXDerivatives.cu | sed 's/\\/\\\\/g' | sed 's/\"/\\"/g' | sed 's/$/\\n\\/' >> BeagleOpenCLSpectral_kernels.h + cat $srcdir/../AdjointBlockSize.h | sed 's/\\/\\\\/g' | sed 's/\"/\\"/g' | sed 's/$/\\n\\/' >> BeagleOpenCLSpectral_kernels.h cat $srcdir/kernelsSpectralIfDef.cu | sed 's/\\/\\\\/g' | sed 's/\"/\\"/g' | sed 's/$/\\n\\/' >> BeagleOpenCLSpectral_kernels.h echo "\"" >> BeagleOpenCLSpectral_kernels.h done @@ -66,6 +67,7 @@ for s in $STATE_COUNT_LIST; do cat $srcdir/kernelsAll.cu | sed 's/\\/\\\\/g' | sed 's/\"/\\"/g' | sed 's/$/\\n\\/' >> BeagleOpenCLSpectral_kernels.h cat $srcdir/kernelsX.cu | sed 's/\\/\\\\/g' | sed 's/\"/\\"/g' | sed 's/$/\\n\\/' >> BeagleOpenCLSpectral_kernels.h cat $srcdir/kernelsXDerivatives.cu | sed 's/\\/\\\\/g' | sed 's/\"/\\"/g' | sed 's/$/\\n\\/' >> BeagleOpenCLSpectral_kernels.h + cat $srcdir/../AdjointBlockSize.h | sed 's/\\/\\\\/g' | sed 's/\"/\\"/g' | sed 's/$/\\n\\/' >> BeagleOpenCLSpectral_kernels.h cat $srcdir/kernelsSpectralIfDef.cu | sed 's/\\/\\\\/g' | sed 's/\"/\\"/g' | sed 's/$/\\n\\/' >> BeagleOpenCLSpectral_kernels.h echo "\"" >> BeagleOpenCLSpectral_kernels.h done