From b7e5435d251d5be3bf935f2c89697964a1cef4c5 Mon Sep 17 00:00:00 2001 From: Zheming Jin Date: Fri, 14 Aug 2026 14:44:31 -0700 Subject: [PATCH 1/2] [blas][cublas] Support int8 inputs with float output in gemm_batch cuBLAS reaches this combination through cublasGemmStridedBatchedEx and cublasGemmBatchedEx, which already accept the datatypes the existing launchers forward, so the column-major buffer, USM strided and USM group entry points only needed routing to the implementation instead of throwing unimplemented. The int32 output combination stays unimplemented because cuBLAS produces it only under CUBLAS_COMPUTE_32I, which takes int32 alpha and beta, whereas oneMath specifies float scalars. A float output accumulated from int8 inputs is rounded at the magnitude of the terms summed rather than at the magnitude of the output entry, so an entry whose sum cancels cannot meet any relative bound. The shared checker takes an optional absolute tolerance, defaulted to zero so that existing callers are unaffected, and the int8-to-float gemm_batch tests pass eps times k * 128 * 128, an upper bound on the accumulated magnitude sum|a*b|. Int8Int8SinglePrecisionErrorModel covers that path with fixed data whose leading rows and columns cancel exactly. Co-authored-by: Cursor --- src/blas/backends/cublas/cublas_batch.cpp | 8 +- .../blas/batch/gemm_batch_stride.cpp | 14 +- .../blas/batch/gemm_batch_stride_usm.cpp | 245 +++++++++++++++++- .../unit_tests/blas/batch/gemm_batch_usm.cpp | 18 +- tests/unit_tests/blas/include/test_common.hpp | 36 ++- 5 files changed, 301 insertions(+), 20 deletions(-) diff --git a/src/blas/backends/cublas/cublas_batch.cpp b/src/blas/backends/cublas/cublas_batch.cpp index 4481195af..b515292fc 100644 --- a/src/blas/backends/cublas/cublas_batch.cpp +++ b/src/blas/backends/cublas/cublas_batch.cpp @@ -199,6 +199,7 @@ inline void gemm_batch_impl(sycl::queue& queue, transpose transa, transpose tran GEMM_STRIDED_BATCH_LAUNCHER(sycl::half, sycl::half, sycl::half, sycl::half) GEMM_STRIDED_BATCH_LAUNCHER(sycl::half, sycl::half, float, float) +GEMM_STRIDED_BATCH_LAUNCHER(std::int8_t, std::int8_t, float, float) GEMM_STRIDED_BATCH_LAUNCHER(float, float, float, float) GEMM_STRIDED_BATCH_LAUNCHER(double, double, double, double) GEMM_STRIDED_BATCH_LAUNCHER(std::complex, std::complex, std::complex, @@ -208,6 +209,8 @@ GEMM_STRIDED_BATCH_LAUNCHER(std::complex, std::complex, std::com #undef GEMM_STRIDED_BATCH_LAUNCHER +// cuBLAS computes an int32 output only with CUBLAS_COMPUTE_32I, which requires int32 alpha and +// beta, whereas oneMath specifies float scalars for this combination. #define GEMM_STRIDED_BATCH_LAUNCHER(TYPE_A, TYPE_B, TYPE_C, TYPE_S) \ void gemm_batch(sycl::queue& queue, transpose transa, transpose transb, int64_t m, int64_t n, \ int64_t k, TYPE_S alpha, sycl::buffer& a, int64_t lda, \ @@ -220,7 +223,6 @@ GEMM_STRIDED_BATCH_LAUNCHER(std::complex, std::complex, std::com dtype_string() + "," + dtype_string() + ">"); \ } -GEMM_STRIDED_BATCH_LAUNCHER(std::int8_t, std::int8_t, float, float) GEMM_STRIDED_BATCH_LAUNCHER(std::int8_t, std::int8_t, std::int32_t, float) #undef GEMM_STRIDED_BATCH_LAUNCHER @@ -668,6 +670,7 @@ inline sycl::event gemm_batch_strided_usm_impl(sycl::queue& queue, transpose tra GEMM_STRIDED_BATCH_LAUNCHER_USM(sycl::half, sycl::half, sycl::half, sycl::half) GEMM_STRIDED_BATCH_LAUNCHER_USM(sycl::half, sycl::half, float, float) +GEMM_STRIDED_BATCH_LAUNCHER_USM(std::int8_t, std::int8_t, float, float) GEMM_STRIDED_BATCH_LAUNCHER_USM(float, float, float, float) GEMM_STRIDED_BATCH_LAUNCHER_USM(double, double, double, double) GEMM_STRIDED_BATCH_LAUNCHER_USM(std::complex, std::complex, std::complex, @@ -689,7 +692,6 @@ GEMM_STRIDED_BATCH_LAUNCHER_USM(std::complex, std::complex, std: dtype_string() + "," + dtype_string() + ">"); \ } -GEMM_STRIDED_BATCH_LAUNCHER_USM(std::int8_t, std::int8_t, float, float) GEMM_STRIDED_BATCH_LAUNCHER_USM(std::int8_t, std::int8_t, std::int32_t, float) #undef GEMM_STRIDED_BATCH_LAUNCHER_USM @@ -761,6 +763,7 @@ inline sycl::event gemm_batch_usm_impl(sycl::queue& queue, transpose* transa, tr GEMM_BATCH_LAUNCHER_USM(sycl::half, sycl::half, sycl::half, sycl::half) GEMM_BATCH_LAUNCHER_USM(sycl::half, sycl::half, float, float) +GEMM_BATCH_LAUNCHER_USM(std::int8_t, std::int8_t, float, float) GEMM_BATCH_LAUNCHER_USM(float, float, float, float) GEMM_BATCH_LAUNCHER_USM(double, double, double, double) GEMM_BATCH_LAUNCHER_USM(std::complex, std::complex, std::complex, @@ -782,7 +785,6 @@ GEMM_BATCH_LAUNCHER_USM(std::complex, std::complex, std::complex dtype_string() + "," + dtype_string() + ">"); \ } -GEMM_BATCH_LAUNCHER_USM(std::int8_t, std::int8_t, float, float) GEMM_BATCH_LAUNCHER_USM(std::int8_t, std::int8_t, std::int32_t, float) #undef GEMM_BATCH_LAUNCHER_USM diff --git a/tests/unit_tests/blas/batch/gemm_batch_stride.cpp b/tests/unit_tests/blas/batch/gemm_batch_stride.cpp index 50e90ccbb..d64f4e540 100644 --- a/tests/unit_tests/blas/batch/gemm_batch_stride.cpp +++ b/tests/unit_tests/blas/batch/gemm_batch_stride.cpp @@ -219,13 +219,25 @@ int test(device* dev, oneapi::math::layout layout, int64_t batch_size) { if (std::is_same_v) error_mag = 1; + // A float output accumulated from int8 inputs is rounded at the magnitude of the terms summed, + // |alpha| * sum|a*b|, which k * 128 * 128 bounds from above. An entry whose sum cancels is far + // smaller than that and so cannot meet any relative bound, so allow an absolute error of eps + // times the accumulated magnitude instead. + constexpr bool int8_to_float = std::is_same_v && + std::is_same_v && std::is_same_v && + std::is_same_v; + double abs_error_bound = 0.0; + if constexpr (int8_to_float) + abs_error_bound = std::numeric_limits::epsilon() * std::abs(double(alpha)) * + double(k) * 128.0 * 128.0; + for (size_t i = 0; i < C_ref.size(); ++i) { C_cast_ref[i] = C_ref[i]; } auto C_accessor = C_buffer.get_host_access(read_only); bool good = check_almost_equal_matrix(C_accessor, C_cast_ref, oneapi::math::layout::col_major, stride_c * batch_size, 1, stride_c * batch_size, - error_mag, std::cout); + error_mag, std::cout, abs_error_bound); return (int)good; } diff --git a/tests/unit_tests/blas/batch/gemm_batch_stride_usm.cpp b/tests/unit_tests/blas/batch/gemm_batch_stride_usm.cpp index 1f46e1d68..e438741a3 100644 --- a/tests/unit_tests/blas/batch/gemm_batch_stride_usm.cpp +++ b/tests/unit_tests/blas/batch/gemm_batch_stride_usm.cpp @@ -250,12 +250,25 @@ int test(device* dev, oneapi::math::layout layout, int64_t batch_size) { if (std::is_same_v) error_mag = 1; + // A float output accumulated from int8 inputs is rounded at the magnitude of the terms summed, + // |alpha| * sum|a*b|, which k * 128 * 128 bounds from above. An entry whose sum cancels is far + // smaller than that and so cannot meet any relative bound, so allow an absolute error of eps + // times the accumulated magnitude instead. Int8Int8SinglePrecisionErrorModel checks that the + // error really does stay inside eps * |alpha| * sum|a*b| on fixed data. + constexpr bool int8_to_float = std::is_same_v && + std::is_same_v && std::is_same_v && + std::is_same_v; + double abs_error_bound = 0.0; + if constexpr (int8_to_float) + abs_error_bound = std::numeric_limits::epsilon() * std::abs(double(alpha)) * + double(k) * 128.0 * 128.0; + for (size_t i = 0; i < C_ref.size(); ++i) { C_cast_ref[i] = C_ref[i]; } bool good = check_almost_equal_matrix(C, C_cast_ref, oneapi::math::layout::col_major, stride_c * batch_size, 1, stride_c * batch_size, - error_mag, std::cout); + error_mag, std::cout, abs_error_bound); oneapi::math::free_shared(a_array, cxt); oneapi::math::free_shared(b_array, cxt); @@ -265,6 +278,231 @@ int test(device* dev, oneapi::math::layout layout, int64_t batch_size) { return (int)good; } +// Regression test for the int8-to-float tolerance above. The sizes and data are fixed rather than +// drawn from std::rand(), so this does not depend on the order the tests run in, and the expected +// result is accumulated exactly in integers, so it does not depend on the reference BLAS either. +// The leading rows of A and columns of B are built to cancel exactly, which puts those entries out +// of reach of any relative bound and leaves the absolute bound as the only one that can accept +// them. Every entry is checked against that pair of bounds, and against the accumulation error +// model the absolute bound is calibrated from: eps times the magnitude of the terms summed. +int int8_accumulation_error_model(device* dev, oneapi::math::layout layout) { + auto exception_handler = [](exception_list exceptions) { + for (std::exception_ptr const& e : exceptions) { + try { + std::rethrow_exception(e); + } + catch (exception const& e) { + std::cout << "Caught asynchronous SYCL exception during GEMM_BATCH_STRIDE:\n" + << e.what() << std::endl; + print_error_code(e); + } + } + }; + + queue main_queue(*dev, exception_handler); + context cxt = main_queue.get_context(); + event done; + std::vector dependencies; + + const auto transa = oneapi::math::transpose::nontrans; + const auto transb = oneapi::math::transpose::nontrans; + // Shapes are not interchangeable here: a backend may accumulate some of them exactly, in which + // case no tolerance is needed and nothing exercises this one. This shape leaves a rounding + // error large enough for the cancelling entries below to fall back on the absolute bound. + const int64_t m = 466, n = 15, batch_size = 2; + const int64_t k = 141; // a multiple of three, for the cancelling triples built below + const int64_t cancelling = 8; // leading rows of A and columns of B that cancel exactly + // alpha is not a power of two, so that scaling the terms of the sum rounds. + const float alpha = 0.3f, beta = 0.25f; + + const bool col = layout == oneapi::math::layout::col_major; + const int64_t lda = col ? m : k, ldb = col ? k : n, ldc = col ? m : n; + const int64_t stride_a = col ? lda * k : lda * m; + const int64_t stride_b = col ? ldb * n : ldb * k; + const int64_t stride_c = col ? ldc * n : ldc * m; + auto a_at = [=](int64_t i, int64_t l) { + return col ? i + l * lda : i * lda + l; + }; + auto b_at = [=](int64_t l, int64_t j) { + return col ? l + j * ldb : l * ldb + j; + }; + auto c_at = [=](int64_t i, int64_t j) { + return col ? i + j * ldc : i * ldc + j; + }; + + auto ua = usm_allocator(cxt, *dev); + auto uc = usm_allocator(cxt, *dev); + vector A(stride_a * batch_size, ua), B(stride_b * batch_size, ua); + vector C(stride_c * batch_size, uc); + std::vector C_in(stride_c * batch_size); + + std::uint32_t seed = 20250814u; + auto next = [&seed]() { + seed = seed * 1664525u + 1013904223u; + return seed >> 16; + }; + for (int64_t b = 0; b < batch_size; b++) { + for (int64_t i = 0; i < m; i++) + for (int64_t l = 0; l < k; l++) + A[b * stride_a + a_at(i, l)] = std::int8_t(int(next() % 254) - 127); + for (int64_t l = 0; l < k; l++) + for (int64_t j = 0; j < n; j++) + B[b * stride_b + b_at(l, j)] = std::int8_t(int(next() % 254) - 127); + // The leading rows of A and columns of B are built from triples whose products are 5x, -3x + // and -2x, so the exact dot product of any such row with any such column is zero while the + // terms summed stay large. None of the three is a power of two times another, so rounding + // the scaled terms does not cancel along with the terms themselves. + for (int64_t i = 0; i < cancelling; i++) + for (int64_t l = 0; l < k; l += 3) { + const int g = int(next() % 25) + 1; + A[b * stride_a + a_at(i, l)] = std::int8_t(5 * g); + A[b * stride_a + a_at(i, l + 1)] = std::int8_t(3 * g); + A[b * stride_a + a_at(i, l + 2)] = std::int8_t(2 * g); + } + for (int64_t j = 0; j < cancelling; j++) + for (int64_t l = 0; l < k; l += 3) { + const int h = (int(next() % 127) + 1) * (next() % 2 ? 1 : -1); + B[b * stride_b + b_at(l, j)] = std::int8_t(h); + B[b * stride_b + b_at(l + 1, j)] = std::int8_t(-h); + B[b * stride_b + b_at(l + 2, j)] = std::int8_t(-h); + } + for (int64_t j = 0; j < n; j++) + for (int64_t i = 0; i < m; i++) { + const auto idx = b * stride_c + c_at(i, j); + C[idx] = float(next() % 1024) / 512.0f - 1.0f; + C_in[idx] = C[idx]; + } + } + + try { +#ifdef CALL_RT_API + switch (layout) { + case oneapi::math::layout::col_major: + done = oneapi::math::blas::column_major::gemm_batch( + main_queue, transa, transb, m, n, k, alpha, &A[0], lda, stride_a, &B[0], ldb, + stride_b, beta, &C[0], ldc, stride_c, batch_size, dependencies); + break; + case oneapi::math::layout::row_major: + done = oneapi::math::blas::row_major::gemm_batch( + main_queue, transa, transb, m, n, k, alpha, &A[0], lda, stride_a, &B[0], ldb, + stride_b, beta, &C[0], ldc, stride_c, batch_size, dependencies); + break; + default: break; + } + done.wait_and_throw(); +#else + switch (layout) { + case oneapi::math::layout::col_major: + TEST_RUN_BLAS_CT_SELECT(main_queue, oneapi::math::blas::column_major::gemm_batch, + transa, transb, m, n, k, alpha, &A[0], lda, stride_a, &B[0], + ldb, stride_b, beta, &C[0], ldc, stride_c, batch_size, + dependencies); + break; + case oneapi::math::layout::row_major: + TEST_RUN_BLAS_CT_SELECT(main_queue, oneapi::math::blas::row_major::gemm_batch, + transa, transb, m, n, k, alpha, &A[0], lda, stride_a, &B[0], + ldb, stride_b, beta, &C[0], ldc, stride_c, batch_size, + dependencies); + break; + default: break; + } + main_queue.wait_and_throw(); +#endif + } + catch (exception const& e) { + std::cout << "Caught synchronous SYCL exception during GEMM_BATCH_STRIDE:\n" + << e.what() << std::endl; + print_error_code(e); + } + + catch (const oneapi::math::unimplemented& e) { + return test_skipped; + } + + catch (const std::runtime_error& error) { + std::cout << "Error raised during execution of GEMM_BATCH_STRIDE:\n" + << error.what() << std::endl; + } + + const double eps = std::numeric_limits::epsilon(); + // The same pair of bounds the int8-to-float tests above apply, evaluated here against an exact + // integer reference: a relative bound of 10 * k * eps, or an absolute one of eps times the + // bound k * 128 * 128 on the accumulated magnitude. + const double relative_bound = double(10 * k) * eps; + const double absolute_bound = eps * std::abs(double(alpha)) * double(k) * 128.0 * 128.0; + double worst_model_usage = 0.0, worst_absolute_usage = 0.0; + double worst_cancelling_relative_allowance = 0.0; + int64_t entries_missing_relative_bound = 0, cancelling_missing_relative_bound = 0; + bool good = true; + for (int64_t b = 0; b < batch_size; b++) { + const std::int8_t* Ab = &A[b * stride_a]; + const std::int8_t* Bb = &B[b * stride_b]; + for (int64_t j = 0; j < n; j++) + for (int64_t i = 0; i < m; i++) { + std::int64_t dot = 0, abs_sum = 0; + for (int64_t l = 0; l < k; l++) { + const std::int64_t a = Ab[a_at(i, l)], bb = Bb[b_at(l, j)]; + dot += a * bb; + abs_sum += std::abs(a * bb); + } + const auto idx = b * stride_c + c_at(i, j); + const double expected = + double(alpha) * double(dot) + double(beta) * double(C_in[idx]); + const double error = std::abs(double(C[idx]) - expected); + const bool cancels = i < cancelling && j < cancelling; + if (cancels && dot != 0) { + std::cout << "test bug: entry (" << i << "," << j + << ") was built to cancel but its dot product is " << dot + << std::endl; + return false; + } + + // The error the accumulation is allowed: eps times the magnitude of the terms + // summed, plus the scaling of C and the final addition. Exceeding this means the + // calibration the absolute tolerance rests on no longer describes the backend. + const double model_bound = + eps * (std::abs(double(alpha)) * double(abs_sum) + + std::abs(double(beta) * double(C_in[idx])) + std::abs(expected)); + worst_model_usage = std::max(worst_model_usage, error / model_bound); + worst_absolute_usage = std::max(worst_absolute_usage, error / absolute_bound); + if (error > model_bound) + good = false; + + if (cancels) + worst_cancelling_relative_allowance = std::max( + worst_cancelling_relative_allowance, relative_bound * std::abs(expected)); + if (error > relative_bound * std::abs(expected)) { + entries_missing_relative_bound++; + if (cancels) + cancelling_missing_relative_bound++; + if (error > absolute_bound) + good = false; + } + } + } + + // The cancelling entries are the ones the absolute bound exists for: whatever error the + // backend makes on them, the relative bound can only accept a fraction of what the model + // permits, so they rest on the absolute bound alone. + if (worst_cancelling_relative_allowance >= absolute_bound) { + std::cout << "test bug: the relative bound already covers the cancelling entries, so they " + "do not exercise the absolute tolerance" + << std::endl; + return false; + } + + std::cout << "int8 accumulation error reached " << worst_model_usage + << " of the accumulated magnitude the model allows and " << worst_absolute_usage + << " of the absolute tolerance; " << entries_missing_relative_bound + << " entries missed the relative bound, " << cancelling_missing_relative_bound + << " of them cancelling" << std::endl; + if (!good) + std::cout << "int8 accumulation error exceeded the tolerance the int8-to-float gemm_batch " + "tests rely on" + << std::endl; + return good; +} + class GemmBatchStrideUsmTests : public ::testing::TestWithParam> {}; @@ -283,6 +521,11 @@ TEST_P(GemmBatchStrideUsmTests, Int8Int8SinglePrecision) { std::get<1>(GetParam()), 5))); } +TEST_P(GemmBatchStrideUsmTests, Int8Int8SinglePrecisionErrorModel) { + EXPECT_TRUEORSKIP( + (int8_accumulation_error_model(std::get<0>(GetParam()), std::get<1>(GetParam())))); +} + TEST_P(GemmBatchStrideUsmTests, Int8Int8Int32Precision) { EXPECT_TRUEORSKIP((test( std::get<0>(GetParam()), std::get<1>(GetParam()), 5))); diff --git a/tests/unit_tests/blas/batch/gemm_batch_usm.cpp b/tests/unit_tests/blas/batch/gemm_batch_usm.cpp index 8c4fd6a37..012132c33 100644 --- a/tests/unit_tests/blas/batch/gemm_batch_usm.cpp +++ b/tests/unit_tests/blas/batch/gemm_batch_usm.cpp @@ -324,6 +324,14 @@ int test(device* dev, oneapi::math::layout layout, int64_t group_count) { // Compare the results of reference implementation and DPC++ implementation. int tol_scalar = 10; + // A float output accumulated from int8 inputs is rounded at the magnitude of the terms summed, + // |alpha| * sum|a*b|, which k * 128 * 128 bounds from above. An entry whose sum cancels is far + // smaller than that and so cannot meet any relative bound, so allow an absolute error of eps + // times the accumulated magnitude instead. + constexpr bool int8_to_float = std::is_same_v && + std::is_same_v && std::is_same_v && + std::is_same_v; + idx = 0; for (i = 0; i < group_count; i++) { for (j = 0; j < group_size[i]; j++) { @@ -331,10 +339,16 @@ int test(device* dev, oneapi::math::layout layout, int64_t group_count) { if (std::is_same_v) error_mag = 1; + double abs_error_bound = 0.0; + if constexpr (int8_to_float) + abs_error_bound = std::numeric_limits::epsilon() * + std::abs(double(alpha[i])) * double(k[i]) * 128.0 * 128.0; + copy_matrix(c_ref_array[idx], layout, oneapi::math::transpose::nontrans, m[i], n[i], ldc[i], c_cast_ref_array[idx]); - good = good && check_almost_equal_matrix(c_array[idx], c_cast_ref_array[idx], layout, - m[i], n[i], ldc[i], error_mag, std::cout); + good = good && + check_almost_equal_matrix(c_array[idx], c_cast_ref_array[idx], layout, m[i], + n[i], ldc[i], error_mag, std::cout, abs_error_bound); idx++; } } diff --git a/tests/unit_tests/blas/include/test_common.hpp b/tests/unit_tests/blas/include/test_common.hpp index 64df0bd76..d4182d9d0 100644 --- a/tests/unit_tests/blas/include/test_common.hpp +++ b/tests/unit_tests/blas/include/test_common.hpp @@ -438,26 +438,36 @@ void rand_tbsv_matrix(vec& M, oneapi::math::layout layout, oneapi::math::uplo up } // Correctness checking. +// A mixed-precision operation can accumulate at magnitudes far above the size of its output +// entries. The rounding error of an entry whose sum cancels is then set by the accumulation scale +// rather than by the entry itself, and no relative bound can cover it. Such callers pass that +// scale as abs_bound, which only ever widens the absolute part of the check. template -typename std::enable_if::value, bool>::type check_equal(fp x, fp x_ref, - int error_mag) { +typename std::enable_if::value, bool>::type check_equal( + fp x, fp x_ref, int error_mag, double abs_bound = 0.0) { using fp_real = typename complex_info::real_type; fp_real bound = (error_mag * num_components() * std::numeric_limits::epsilon()); + fp_real abs_limit = std::max(bound, fp_real(abs_bound)); bool ok; fp_real aerr = std::abs(x - x_ref); fp_real rerr = aerr / std::abs(x_ref); - ok = (rerr <= bound) || (aerr <= bound); - if (!ok) + ok = (rerr <= bound) || (aerr <= abs_limit); + if (!ok) { std::cout << "relative error = " << rerr << " absolute error = " << aerr - << " limit = " << bound << std::endl; + << " limit = " << bound; + if (abs_limit > bound) + std::cout << " absolute limit = " << abs_limit; + std::cout << std::endl; + } return ok; } +// An integer result must match exactly, so both tolerances are ignored here. template -typename std::enable_if::value, bool>::type check_equal(fp x, fp x_ref, - int error_mag) { +typename std::enable_if::value, bool>::type check_equal( + fp x, fp x_ref, int error_mag, double abs_bound = 0.0) { return (x == x_ref); } @@ -566,13 +576,13 @@ bool check_equal_trsv_vector(vec1& v, vec2& v_ref, int n, int inc, int error_mag template bool check_equal_matrix(acc1& M, acc2& M_ref, oneapi::math::layout layout, int m, int n, int ld, - int error_mag, std::ostream& out) { + int error_mag, std::ostream& out, double abs_bound = 0.0) { bool good = true; int idx, count = 0; for (int j = 0; j < n; j++) { for (int i = 0; i < m; i++) { idx = (layout == oneapi::math::layout::col_major) ? i + j * ld : j + i * ld; - if (!check_equal(M[idx], M_ref[idx], error_mag)) { + if (!check_equal(M[idx], M_ref[idx], error_mag, abs_bound)) { out << "Difference in entry (" << i << ',' << j << "): DPC++ " << M[idx] << " vs. Reference " << M_ref[idx] << std::endl; good = false; @@ -588,13 +598,13 @@ bool check_equal_matrix(acc1& M, acc2& M_ref, oneapi::math::layout layout, int m template bool check_equal_matrix(const fp* M, const fp* M_ref, oneapi::math::layout layout, int m, int n, - int ld, int error_mag, std::ostream& out) { + int ld, int error_mag, std::ostream& out, double abs_bound = 0.0) { bool good = true; int idx, count = 0; for (int j = 0; j < n; j++) { for (int i = 0; i < m; i++) { idx = (layout == oneapi::math::layout::col_major) ? i + j * ld : j + i * ld; - if (!check_equal(M[idx], M_ref[idx], error_mag)) { + if (!check_equal(M[idx], M_ref[idx], error_mag, abs_bound)) { out << "Difference in entry (" << i << ',' << j << "): DPC++ " << M[idx] << " vs. Reference " << M_ref[idx] << std::endl; good = false; @@ -702,11 +712,11 @@ bool check_almost_equal_matrix_int(Ta& M, Tb& M_ref, oneapi::math::layout layout template bool check_almost_equal_matrix(Ta& M, Tb& M_ref, oneapi::math::layout layout, int m, int n, int ld, - int error_mag, std::ostream& out) { + int error_mag, std::ostream& out, double abs_bound = 0.0) { // Only call if returned dtype is integral if constexpr (is_matrix_type_integral() && is_matrix_type_integral()) return check_almost_equal_matrix_int(M, M_ref, layout, m, n, ld, error_mag, out); - return check_equal_matrix(M, M_ref, layout, m, n, ld, error_mag, out); + return check_equal_matrix(M, M_ref, layout, m, n, ld, error_mag, out, abs_bound); } #endif /* header guard */ From f95fdbbc56cf790c4e76238f5a4936c1a567cc83 Mon Sep 17 00:00:00 2001 From: Zheming Jin Date: Fri, 14 Aug 2026 19:15:54 -0700 Subject: [PATCH 2/2] [blas][tests] Guard the int8 error-model diagnostic against a zero model bound An entry whose terms and stored C value are all zero gives a zero model bound, which the reported usage ratio would divide by. Co-authored-by: Cursor --- tests/unit_tests/blas/batch/gemm_batch_stride_usm.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/tests/unit_tests/blas/batch/gemm_batch_stride_usm.cpp b/tests/unit_tests/blas/batch/gemm_batch_stride_usm.cpp index e438741a3..2cde6dd2a 100644 --- a/tests/unit_tests/blas/batch/gemm_batch_stride_usm.cpp +++ b/tests/unit_tests/blas/batch/gemm_batch_stride_usm.cpp @@ -463,7 +463,8 @@ int int8_accumulation_error_model(device* dev, oneapi::math::layout layout) { const double model_bound = eps * (std::abs(double(alpha)) * double(abs_sum) + std::abs(double(beta) * double(C_in[idx])) + std::abs(expected)); - worst_model_usage = std::max(worst_model_usage, error / model_bound); + worst_model_usage = + std::max(worst_model_usage, model_bound > 0.0 ? error / model_bound : 0.0); worst_absolute_usage = std::max(worst_absolute_usage, error / absolute_bound); if (error > model_bound) good = false;