From 61115facca2b56d249d0e7456e1776f664b2e230 Mon Sep 17 00:00:00 2001 From: wegamekinglc Date: Sun, 28 Jun 2026 17:28:51 +0800 Subject: [PATCH 1/2] Fuse CG/BCG Krylov sweeps and banded TriMultiply neighbour passes - P7: PrepareDirection_ and UpdateSolution_ in bcg.cpp collapse their separate scale/add (and LinearIncrement) sweeps into single AXPY passes over p/z/x/r (plus the bi-conjugate pp/zz/rr shadows). Operation order is preserved (p*multiply + z; dst + scale*src) so -ffp-contract=fast forms the same FMAs as the unfused two-pass code, keeping the result numerically faithful on the FP-sensitive joint-calibration path. - P8: TriMultiply in banded.cpp keeps the dominant diag*x Transform pass (vectorization-identical to the unfused code) and fuses the two strided neighbour-correction loops into one sweep. A full single-pass interleave was tested and rejected: it re-vectorizes the diag*x term, perturbing FP rounding enough to flip JointCalibrationTest.TestJointOisCurveAgreesWithStagedOis from 7e-6 (unfused) to 5.7e-4 (past the 1e-5 bar); the hybrid stays at 4.6e-6. - Add residual-bound regression tests for CG (symmetric) and BCG (asymmetric) Krylov solves, and an asymmetric TriMultiply test covering all three bands and both boundaries. Benchmarks (GCC -O3 -march=native -ffp-contract=fast): CGSolve 500x500 15.24us -> 13.95us (-8.5%) BCGSolve 500x500 19.40us -> 17.64us (-9.0%) TriDecomp MultiplyLeft 10K 5.53us -> 4.34us (-21.6%) TriDiagonal MultiplyLeft 10K 5.33us -> 4.27us (-19.7%) Co-Authored-By: Claude --- dal-cpp/dal/math/matrix/banded.cpp | 31 ++++++++--- dal-cpp/dal/math/matrix/bcg.cpp | 54 +++++++++++++++---- dal-cpp/tests/math/matrix/test_banded.cpp | 33 ++++++++++++ dal-cpp/tests/math/matrix/test_bcg.cpp | 66 +++++++++++++++++++++++ 4 files changed, 167 insertions(+), 17 deletions(-) diff --git a/dal-cpp/dal/math/matrix/banded.cpp b/dal-cpp/dal/math/matrix/banded.cpp index 8c99514eb..cd472fe32 100644 --- a/dal-cpp/dal/math/matrix/banded.cpp +++ b/dal-cpp/dal/math/matrix/banded.cpp @@ -17,13 +17,32 @@ namespace Dal { const Vector_<>& x, const Vector_<>& diag, const Vector_<>& above, const Vector_<>& below, Vector_<>* r) { REQUIRE(x.size() == diag.size(), ""); r->Resize(x.size()); + // Pass 1 (vectorizable elementwise, identical to the unfused Transform) seeds r = diag * x. Transform(x, diag, std::multiplies<>(), r); - auto pr = r->begin(); - for (auto px = x.begin() + 1, pa = above.begin(); pa != above.end(); ++px, ++pa, ++pr) - *pr += *px * *pa; - pr = r->begin() + 1; - for (auto px = x.begin(), pb = below.begin(); pb != below.end(); ++px, ++pb, ++pr) - *pr += *px * *pb; + // Pass 2 fuses the two strided neighbour corrections into one sweep: + // r[i] += above[i] * x[i+1] (superdiagonal, i = 0..n-2) + // r[i] += below[i-1] * x[i-1] (subdiagonal, i = 1..n-1) + // The accumulator order (above term first, then below) matches the unfused code so + // -ffp-contract=fast forms the same FMAs. + const int n = static_cast(x.size()); + if (n < 2) + return; + auto xNext = x.begin() + 1; // x[i+1], walks with i + auto xPrev = x.begin(); // x[i-1], lags by one; valid for i >= 1 + auto rIt = r->begin(); + auto aboveIt = above.begin(); + auto belowPrev = below.begin(); // below[i-1], lags by one; valid for i >= 1 + // i = 0: superdiagonal only (below has no i-1 term yet). + *rIt += *aboveIt * *xNext; + ++rIt; + ++aboveIt; + ++xNext; + for (int ii = 1; ii < n - 1; ++ii, ++rIt, ++aboveIt, ++belowPrev, ++xNext, ++xPrev) { + *rIt += *aboveIt * *xNext; + *rIt += *belowPrev * *xPrev; + } + // i = n-1: subdiagonal only (above has no term at the last row). + *rIt += *belowPrev * *xPrev; } Vector_<> TridagBetaInverse(const Vector_<>& diag, const Vector_<>& above, const Vector_<>& below) { diff --git a/dal-cpp/dal/math/matrix/bcg.cpp b/dal-cpp/dal/math/matrix/bcg.cpp index 5d07759ae..877229a3f 100644 --- a/dal-cpp/dal/math/matrix/bcg.cpp +++ b/dal-cpp/dal/math/matrix/bcg.cpp @@ -4,11 +4,11 @@ #include #include +#include #include #include #include #include -#include namespace Dal { namespace { @@ -71,12 +71,24 @@ namespace Dal { s.precondition.Right(s.rr, &s.zz); const double beta = InnerProduct(s.zzRef, s.r); const double multiply = ii > 0 ? beta / s.betaPrev : 0.0; - s.p *= multiply; - if (s.biConjugate) - s.pp *= multiply; - s.p += s.z; - if (s.biConjugate) - s.pp += s.zz; + // Fused AXPY sweep mirroring the unfused operation order (p *= multiply; p += z) + // so -ffp-contract=fast forms the same FMAs as the original two-pass code. + { + auto pIt = s.p.begin(); + auto zIt = s.z.begin(); + const auto pEnd = s.p.end(); + if (s.biConjugate) { + auto ppIt = s.pp.begin(); + auto zzIt = s.zz.begin(); + for (; pIt != pEnd; ++pIt, ++zIt, ++ppIt, ++zzIt) { + *pIt = *pIt * multiply + *zIt; + *ppIt = *ppIt * multiply + *zzIt; + } + } else { + for (; pIt != pEnd; ++pIt, ++zIt) + *pIt = *pIt * multiply + *zIt; + } + } s.betaPrev = beta; return beta; } @@ -86,10 +98,30 @@ namespace Dal { if (s.biConjugate) s.A.MultiplyRight(s.pp, &s.zz); const double alphaK = beta / InnerProduct(s.z, s.ppRef); - Transform(x, s.p, LinearIncrement(alphaK)); - Transform(&s.r, s.z, LinearIncrement(-alphaK)); - if (s.biConjugate) - Transform(&s.rr, s.zz, LinearIncrement(-alphaK)); + // Fused AXPY sweep mirroring LinearIncrement(scale) = dst + scale * src, so the + // FMA pattern matches the unfused Transform calls (r + (-alphaK)*z, not r - alphaK*z). + const double negAlphaK = -alphaK; + { + auto xIt = x->begin(); + auto pIt = s.p.begin(); + auto rIt = s.r.begin(); + auto zIt = s.z.begin(); + const auto xEnd = x->end(); + if (s.biConjugate) { + auto rrIt = s.rr.begin(); + auto zzIt = s.zz.begin(); + for (; xIt != xEnd; ++xIt, ++pIt, ++rIt, ++zIt, ++rrIt, ++zzIt) { + *xIt = *xIt + alphaK * *pIt; + *rIt = *rIt + negAlphaK * *zIt; + *rrIt = *rrIt + negAlphaK * *zzIt; + } + } else { + for (; xIt != xEnd; ++xIt, ++pIt, ++rIt, ++zIt) { + *xIt = *xIt + alphaK * *pIt; + *rIt = *rIt + negAlphaK * *zIt; + } + } + } } void ValidateKrylovParams_(int n, const Vector_<>& b, const Vector_<>* x, double tolRel, double tolAbs, int maxIterations) { diff --git a/dal-cpp/tests/math/matrix/test_banded.cpp b/dal-cpp/tests/math/matrix/test_banded.cpp index e08e0d1fa..fed974d5f 100644 --- a/dal-cpp/tests/math/matrix/test_banded.cpp +++ b/dal-cpp/tests/math/matrix/test_banded.cpp @@ -108,3 +108,36 @@ TEST(MatrixTest, TestTriDiagonalSolve) { for (int i = 0; i < n; ++i) ASSERT_NEAR(calculated[i], expected[i], 1e-6); } + +// Varying x with distinct, asymmetric above/below bands exercises every branch of +// the fused TriMultiply (interior points combine all three bands; both boundaries +// drop one term). Hand-computed expected values lock the fused contract. +TEST(MatrixTest, TestTriDiagonalMultiplyAsymmetricFused) { + const auto n = 5; + Sparse::TriDiagonal_ trig(n); + trig.Set(0, 0, 2.0); + trig.Set(0, 1, 3.0); // above[0] + trig.Set(1, 0, 4.0); // below[0] + trig.Set(1, 1, 5.0); + trig.Set(1, 2, 6.0); // above[1] + trig.Set(2, 1, 7.0); // below[1] + trig.Set(2, 2, 8.0); + trig.Set(2, 3, 9.0); // above[2] + trig.Set(3, 2, 10.0); // below[2] + trig.Set(3, 3, 11.0); + trig.Set(3, 4, 12.0); // above[3] + trig.Set(4, 3, 13.0); // below[3] + trig.Set(4, 4, 14.0); + + Vector_<> x = {1.0, 2.0, 3.0, 4.0, 5.0}; + Vector_<> calculated(n); + trig.MultiplyLeft(x, &calculated); + // r[0] = diag[0]*x[0] + above[0]*x[1] + // r[1] = diag[1]*x[1] + above[1]*x[2] + below[0]*x[0] + // r[2] = diag[2]*x[2] + above[2]*x[3] + below[1]*x[1] + // r[3] = diag[3]*x[3] + above[3]*x[4] + below[2]*x[2] + // r[4] = diag[4]*x[4] + below[3]*x[3] + Vector_<> expected = {2 * 1 + 3 * 2, 5 * 2 + 6 * 3 + 4 * 1, 8 * 3 + 9 * 4 + 7 * 2, 11 * 4 + 12 * 5 + 10 * 3, 14 * 5 + 13 * 4}; + for (int i = 0; i < n; ++i) + ASSERT_DOUBLE_EQ(calculated[i], expected[i]); +} diff --git a/dal-cpp/tests/math/matrix/test_bcg.cpp b/dal-cpp/tests/math/matrix/test_bcg.cpp index deeec8e7a..2cf0b2b16 100644 --- a/dal-cpp/tests/math/matrix/test_bcg.cpp +++ b/dal-cpp/tests/math/matrix/test_bcg.cpp @@ -3,6 +3,8 @@ // #include +#include +#include #include #include #include @@ -43,3 +45,67 @@ TEST(MatrixTest, TestBCGSolve) { ASSERT_NEAR(result[9], 0.07446809, 1e-8); } +// Tri-diagonal systems with distinct band values exercise every branch of the +// fused Krylov sweeps. CG requires a symmetric positive-definite matrix, so it +// gets a symmetric system; BCG handles the asymmetric case (and its shadow path). +// The contract is residual ||Ax - b||_inf near machine precision after a tight solve. +namespace { + void BuildSymmetricTridiag(Sparse::Square_* mat, int n) { + for (int i = 0; i < n; ++i) { + mat->Set(i, i, 4.0 + 0.1 * static_cast(i)); + if (i > 0) { + const double off = -0.5 - 0.01 * static_cast(i); + mat->Set(i, i - 1, off); + mat->Set(i - 1, i, off); + } + } + } + + void BuildAsymmetricTridiag(Sparse::Square_* mat, int n) { + for (int i = 0; i < n; ++i) { + mat->Set(i, i, 4.0 + 0.1 * static_cast(i)); + if (i > 0) + mat->Set(i, i - 1, -0.5 - 0.01 * static_cast(i)); + if (i + 1 < n) + mat->Set(i, i + 1, -0.3 + 0.02 * static_cast(i)); + } + } + + double ResidualInfNorm(const Sparse::Square_& A, const Vector_<>& x, const Vector_<>& b) { + Vector_<> ax(x.size()); + A.MultiplyLeft(x, &ax); + double worst = 0.0; + for (int i = 0; i < static_cast(ax.size()); ++i) + worst = std::max(worst, std::fabs(ax[i] - b[i])); + return worst; + } +} // namespace + +TEST(MatrixTest, TestCGSolveSymmetricLowResidual) { + const int n = 40; + Sparse::Square_* mat = Sparse::NewBandDiagonal(n, 1, 1); + BuildSymmetricTridiag(mat, n); + + Vector_<> b(n); + for (int i = 0; i < n; ++i) + b[i] = 1.0 + 0.1 * static_cast(i % 7); + Vector_<> x(n, 0.0); + + Sparse::CGSolve(*mat, b, 1e-12, 1e-14, 500, &x); + ASSERT_LT(ResidualInfNorm(*mat, x, b), 1e-8); +} + +TEST(MatrixTest, TestBCGSolveAsymmetricLowResidual) { + const int n = 40; + Sparse::Square_* mat = Sparse::NewBandDiagonal(n, 1, 1); + BuildAsymmetricTridiag(mat, n); + + Vector_<> b(n); + for (int i = 0; i < n; ++i) + b[i] = 1.0 + 0.1 * static_cast(i % 7); + Vector_<> x(n, 0.0); + + Sparse::BCGSolve(*mat, b, 1e-12, 1e-14, 500, &x); + ASSERT_LT(ResidualInfNorm(*mat, x, b), 1e-8); +} + From 923c0d024f86278c802bc4af2be937c75258d4b7 Mon Sep 17 00:00:00 2001 From: wegamekinglc Date: Sun, 28 Jun 2026 18:30:29 +0800 Subject: [PATCH 2/2] revert: drop P7 CG/BCG Krylov fusion (breaks clang auto-vectorization FP) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit P7 changes clang's auto-vectorization pattern (different from GCC even with -ffp-contract=fast on both), pushing JointCalibrationTest.TestJointOisCurveAgreesWithStagedOis drift from 7e-6 to 2.17e-5 (over the 1e-5 bar). P8 (banded TriMultiply hybrid fusion) is kept — it passes on all compilers. Co-Authored-By: Claude --- dal-cpp/dal/math/matrix/bcg.cpp | 54 +++++++-------------------------- 1 file changed, 11 insertions(+), 43 deletions(-) diff --git a/dal-cpp/dal/math/matrix/bcg.cpp b/dal-cpp/dal/math/matrix/bcg.cpp index 877229a3f..5d07759ae 100644 --- a/dal-cpp/dal/math/matrix/bcg.cpp +++ b/dal-cpp/dal/math/matrix/bcg.cpp @@ -4,11 +4,11 @@ #include #include -#include #include #include #include #include +#include namespace Dal { namespace { @@ -71,24 +71,12 @@ namespace Dal { s.precondition.Right(s.rr, &s.zz); const double beta = InnerProduct(s.zzRef, s.r); const double multiply = ii > 0 ? beta / s.betaPrev : 0.0; - // Fused AXPY sweep mirroring the unfused operation order (p *= multiply; p += z) - // so -ffp-contract=fast forms the same FMAs as the original two-pass code. - { - auto pIt = s.p.begin(); - auto zIt = s.z.begin(); - const auto pEnd = s.p.end(); - if (s.biConjugate) { - auto ppIt = s.pp.begin(); - auto zzIt = s.zz.begin(); - for (; pIt != pEnd; ++pIt, ++zIt, ++ppIt, ++zzIt) { - *pIt = *pIt * multiply + *zIt; - *ppIt = *ppIt * multiply + *zzIt; - } - } else { - for (; pIt != pEnd; ++pIt, ++zIt) - *pIt = *pIt * multiply + *zIt; - } - } + s.p *= multiply; + if (s.biConjugate) + s.pp *= multiply; + s.p += s.z; + if (s.biConjugate) + s.pp += s.zz; s.betaPrev = beta; return beta; } @@ -98,30 +86,10 @@ namespace Dal { if (s.biConjugate) s.A.MultiplyRight(s.pp, &s.zz); const double alphaK = beta / InnerProduct(s.z, s.ppRef); - // Fused AXPY sweep mirroring LinearIncrement(scale) = dst + scale * src, so the - // FMA pattern matches the unfused Transform calls (r + (-alphaK)*z, not r - alphaK*z). - const double negAlphaK = -alphaK; - { - auto xIt = x->begin(); - auto pIt = s.p.begin(); - auto rIt = s.r.begin(); - auto zIt = s.z.begin(); - const auto xEnd = x->end(); - if (s.biConjugate) { - auto rrIt = s.rr.begin(); - auto zzIt = s.zz.begin(); - for (; xIt != xEnd; ++xIt, ++pIt, ++rIt, ++zIt, ++rrIt, ++zzIt) { - *xIt = *xIt + alphaK * *pIt; - *rIt = *rIt + negAlphaK * *zIt; - *rrIt = *rrIt + negAlphaK * *zzIt; - } - } else { - for (; xIt != xEnd; ++xIt, ++pIt, ++rIt, ++zIt) { - *xIt = *xIt + alphaK * *pIt; - *rIt = *rIt + negAlphaK * *zIt; - } - } - } + Transform(x, s.p, LinearIncrement(alphaK)); + Transform(&s.r, s.z, LinearIncrement(-alphaK)); + if (s.biConjugate) + Transform(&s.rr, s.zz, LinearIncrement(-alphaK)); } void ValidateKrylovParams_(int n, const Vector_<>& b, const Vector_<>* x, double tolRel, double tolAbs, int maxIterations) {