From b7e5435d251d5be3bf935f2c89697964a1cef4c5 Mon Sep 17 00:00:00 2001 From: Zheming Jin Date: Fri, 14 Aug 2026 14:44:31 -0700 Subject: [PATCH 1/3] [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/3] [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; From a081fe0d7442b98c68e2b6429cc1dfac5c4c4969 Mon Sep 17 00:00:00 2001 From: Zheming Jin Date: Fri, 14 Aug 2026 19:18:20 -0700 Subject: [PATCH 3/3] [blas][rocblas] Support int8 inputs with float output in gemm_batch rocBLAS reaches int8 inputs only with an int32 output and compute type, which takes int32 alpha and beta, whereas oneMath specifies a float output and float scalars for this combination. Accumulate the products exactly in an int32 workspace and apply the float scalars in a scaling kernel afterwards, for the buffer strided, USM strided and USM grouped entry points in both layouts. Bound k so the int32 accumulator cannot wrap, reject sizes whose workspace or kernel range would overflow, and reject an ldc below the row count, which the scaling kernel would otherwise fold onto the next column. Grouped scaling stays a single kernel by locating the group that owns each entry with a binary search over per-group metadata. The scaling kernels launch over a flat range, since HIP cannot map every large multi-dimensional range onto its grid, and the added regression holds a large prime in n to cover that. Co-authored-by: Cursor --- src/blas/backends/rocblas/rocblas_batch.cpp | 418 +++++++++++++++++- .../blas/batch/gemm_batch_stride_usm.cpp | 100 +++++ 2 files changed, 512 insertions(+), 6 deletions(-) diff --git a/src/blas/backends/rocblas/rocblas_batch.cpp b/src/blas/backends/rocblas/rocblas_batch.cpp index b6e550724..43ad62c15 100644 --- a/src/blas/backends/rocblas/rocblas_batch.cpp +++ b/src/blas/backends/rocblas/rocblas_batch.cpp @@ -25,6 +25,9 @@ #include "oneapi/math/exceptions.hpp" #include "oneapi/math/blas/detail/rocblas/onemath_blas_rocblas.hpp" +#include +#include + // Helper Functions template @@ -69,6 +72,120 @@ namespace blas { namespace rocblas { namespace column_major { +inline void check_int8_float_nonnegative(const char* name, int64_t value) { + if (value < 0) { + throw invalid_argument("blas", "gemm_batch", std::string(name) + " must be nonnegative"); + } +} + +// rocBLAS reaches int8 inputs only with an int32 output and compute type, so the int8-to-float +// combination accumulates exactly in int32 and applies oneMath's float alpha and beta afterwards. +// Every int8 magnitude is at most 128, so 128 * 128 bounds a single product and a larger k than +// this would wrap the int32 accumulator rocBLAS requires. +inline void check_int8_float_accumulation_size(int64_t k) { + constexpr int64_t max_product = 128 * 128; + constexpr int64_t max_safe_k = std::numeric_limits::max() / max_product; + if (k > max_safe_k) { + throw unimplemented("blas", "gemm_batch", + "for int8 inputs with a float output and k above " + + std::to_string(max_safe_k) + + ", which would overflow the int32 accumulator"); + } +} + +// The scaling kernels below reach an entry of C at column * ldc + row, so an ldc below the row +// count would fold one column onto the next and read past the workspace rocBLAS filled. rocBLAS +// rejects such a call itself, but only from its host task, which does not hold back the kernel. +inline void check_int8_float_output_leading_dimension(int64_t rows, int64_t ld) { + if (ld < rows) { + throw invalid_argument("blas", "gemm_batch", "ldc is smaller than the number of rows of C"); + } +} + +inline int64_t checked_int8_float_product(int64_t lhs, int64_t rhs, const char* description) { + if (lhs != 0 && rhs > std::numeric_limits::max() / lhs) { + throw invalid_argument("blas", "gemm_batch", + std::string(description) + " exceeds the supported size"); + } + return lhs * rhs; +} + +inline int64_t checked_int8_float_sum(int64_t lhs, int64_t rhs, const char* description) { + if (rhs > std::numeric_limits::max() - lhs) { + throw invalid_argument("blas", "gemm_batch", + std::string(description) + " exceeds the supported size"); + } + return lhs + rhs; +} + +inline int64_t checked_int8_float_matrix_elements(int64_t rows, int64_t columns) { + return checked_int8_float_product(rows, columns, "matrix element count"); +} + +inline std::size_t checked_int8_float_size_t(int64_t value, const char* description) { + if (static_cast(value) > + static_cast(std::numeric_limits::max())) { + throw invalid_argument("blas", "gemm_batch", + std::string(description) + " exceeds the supported size"); + } + return static_cast(value); +} + +inline int64_t add_int8_float_workspace_elements(int64_t total, int64_t batch_count, int64_t ld, + int64_t columns) { + const int64_t per_matrix = checked_int8_float_product(ld, columns, "workspace matrix size"); + const int64_t group_elements = + checked_int8_float_product(batch_count, per_matrix, "workspace group size"); + return checked_int8_float_sum(total, group_elements, "workspace size"); +} + +inline std::size_t checked_int8_float_entries(int64_t rows, int64_t columns, int64_t batch_count) { + const int64_t per_matrix = checked_int8_float_matrix_elements(rows, columns); + const int64_t entries = + checked_int8_float_product(per_matrix, batch_count, "scaling kernel range"); + return checked_int8_float_size_t(entries, "scaling kernel range"); +} + +// A call with no output entries still owes the caller an event that tracks its dependencies. +inline sycl::event int8_float_empty_event(sycl::queue& queue, + const std::vector& dependencies) { + return queue.submit([&](sycl::handler& cgh) { + cgh.depends_on(dependencies); + cgh.host_task([]() {}); + }); +} + +// The scaling kernels below use a single flat range. Some HIP runtimes cannot map particular large +// multi-dimensional global ranges to their per-dimension grid limits. Each work item maps its +// linear index back to the padded layout of C. +struct int8_float_entry { + int64_t batch; + int64_t element; +}; + +inline int8_float_entry int8_float_entry_of(std::size_t linear_index, int64_t rows, int64_t columns, + int64_t ld) { + const int64_t linear = static_cast(linear_index); + const int64_t per_matrix = rows * columns; + const int64_t batch = linear / per_matrix; + const int64_t within = linear - batch * per_matrix; + const int64_t column = within / rows; + const int64_t row = within - column * rows; + return { batch, column * ld + row }; +} + +struct int8_float_group_metadata { + int64_t entry_begin; + int64_t entry_end; + int64_t matrix_begin; + int64_t workspace_begin; + int64_t rows; + int64_t columns; + int64_t ld; + float alpha; + float beta; +}; + // Buffer APIs template @@ -265,6 +382,62 @@ inline void gemm_batch_impl(sycl::queue& queue, transpose transa, transpose tran }); } +inline void gemm_batch_int8_float_impl(sycl::queue& queue, transpose transa, transpose transb, + int64_t m, int64_t n, int64_t k, float alpha, + sycl::buffer& a, int64_t lda, + int64_t stridea, sycl::buffer& b, + int64_t ldb, int64_t strideb, float beta, + sycl::buffer& c, int64_t ldc, int64_t stridec, + int64_t batch_size) { + check_int8_float_nonnegative("m", m); + check_int8_float_nonnegative("n", n); + check_int8_float_nonnegative("k", k); + check_int8_float_nonnegative("lda", lda); + check_int8_float_nonnegative("ldb", ldb); + check_int8_float_nonnegative("ldc", ldc); + check_int8_float_nonnegative("stridea", stridea); + check_int8_float_nonnegative("strideb", strideb); + check_int8_float_nonnegative("stridec", stridec); + check_int8_float_nonnegative("batch_size", batch_size); + overflow_check(m, n, k, lda, ldb, ldc, stridea, strideb, stridec, batch_size); + check_int8_float_accumulation_size(k); + if (m == 0 || n == 0 || batch_size == 0) { + return; + } + check_int8_float_output_leading_dimension(m, ldc); + + // The int32 workspace has to outlive this call: a plain local buffer would wait for the scaling + // kernel in its destructor and make this one type combination synchronous, so a host task holds + // the last reference to it instead. + auto accum = std::make_shared>(c.get_range()); + const std::int32_t accumulate = alpha == 0.0f ? 0 : 1; + constexpr std::int32_t discard_c = 0; + gemm_batch_impl(queue, transa, transb, m, n, k, accumulate, a, lda, stridea, b, ldb, strideb, + discard_c, *accum, ldc, stridec, batch_size); + + const auto entries = checked_int8_float_entries(m, n, batch_size); + auto done = queue.submit([&](sycl::handler& cgh) { + auto accum_acc = accum->get_access(cgh); + auto c_acc = c.get_access(cgh); + cgh.parallel_for(sycl::range<1>{ entries }, [=](sycl::id<1> index) { + const auto entry = int8_float_entry_of(index[0], m, n, ldc); + const auto offset = entry.batch * stridec + entry.element; + float result = 0.0f; + if (alpha != 0.0f) { + result = alpha * static_cast(accum_acc[offset]); + } + if (beta != 0.0f) { + result += beta * c_acc[offset]; + } + c_acc[offset] = result; + }); + }); + queue.submit([&](sycl::handler& cgh) { + cgh.depends_on(done); + cgh.host_task([accum]() {}); + }); +} + #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, \ @@ -286,6 +459,15 @@ GEMM_STRIDED_BATCH_LAUNCHER(sycl::half, sycl::half, float, float) #undef GEMM_STRIDED_BATCH_LAUNCHER +void gemm_batch(sycl::queue& queue, transpose transa, transpose transb, int64_t m, int64_t n, + int64_t k, float alpha, sycl::buffer& a, int64_t lda, + int64_t stridea, sycl::buffer& b, int64_t ldb, int64_t strideb, + float beta, sycl::buffer& c, int64_t ldc, int64_t stridec, + int64_t batch_size) { + gemm_batch_int8_float_impl(queue, transa, transb, m, n, k, alpha, a, lda, stridea, b, ldb, + strideb, beta, c, ldc, stridec, batch_size); +} + #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, \ @@ -298,7 +480,8 @@ GEMM_STRIDED_BATCH_LAUNCHER(sycl::half, sycl::half, float, float) dtype_string() + "," + dtype_string() + ">"); \ } -GEMM_STRIDED_BATCH_LAUNCHER(std::int8_t, std::int8_t, float, float) +// An int32 output reaches rocBLAS only with an int32 compute type, which takes int32 alpha and beta, +// whereas oneMath specifies float scalars for this combination. GEMM_STRIDED_BATCH_LAUNCHER(std::int8_t, std::int8_t, std::int32_t, float) #undef GEMM_STRIDED_BATCH_LAUNCHER @@ -885,6 +1068,67 @@ inline sycl::event gemm_batch_strided_usm_impl(sycl::queue& queue, transpose tra return done; } +inline sycl::event gemm_batch_strided_usm_int8_float_impl( + sycl::queue& queue, transpose transa, transpose transb, int64_t m, int64_t n, int64_t k, + float alpha, const std::int8_t* a, int64_t lda, int64_t stridea, const std::int8_t* b, + int64_t ldb, int64_t strideb, float beta, float* c, int64_t ldc, int64_t stridec, + int64_t batch_size, const std::vector& dependencies) { + check_int8_float_nonnegative("m", m); + check_int8_float_nonnegative("n", n); + check_int8_float_nonnegative("k", k); + check_int8_float_nonnegative("lda", lda); + check_int8_float_nonnegative("ldb", ldb); + check_int8_float_nonnegative("ldc", ldc); + check_int8_float_nonnegative("stridea", stridea); + check_int8_float_nonnegative("strideb", strideb); + check_int8_float_nonnegative("stridec", stridec); + check_int8_float_nonnegative("batch_size", batch_size); + overflow_check(m, n, k, lda, ldb, ldc, stridea, strideb, stridec, batch_size); + check_int8_float_accumulation_size(k); + if (m == 0 || n == 0 || batch_size == 0) { + return int8_float_empty_event(queue, dependencies); + } + check_int8_float_output_leading_dimension(m, ldc); + + // The workspace holds one int32 matrix per batch at the stride the caller uses for C. + const auto previous_batches = + checked_int8_float_product(stridec, batch_size - 1, "strided workspace size"); + const auto accum_size = + add_int8_float_workspace_elements(previous_batches, /*batch_count*/ 1, ldc, n); + auto* accum = sycl::malloc_device( + checked_int8_float_size_t(accum_size, "workspace size"), queue); + if (accum == nullptr) { + throw device_bad_alloc("blas", "gemm_batch", queue.get_device()); + } + + const std::int32_t accumulate = alpha == 0.0f ? 0 : 1; + constexpr std::int32_t discard_c = 0; + auto gemm_done = gemm_batch_strided_usm_impl(queue, transa, transb, m, n, k, accumulate, a, lda, + stridea, b, ldb, strideb, discard_c, accum, ldc, + stridec, batch_size, dependencies); + const auto entries = checked_int8_float_entries(m, n, batch_size); + auto done = queue.submit([&](sycl::handler& cgh) { + cgh.depends_on(gemm_done); + cgh.parallel_for(sycl::range<1>{ entries }, [=](sycl::id<1> index) { + const auto entry = int8_float_entry_of(index[0], m, n, ldc); + const auto offset = entry.batch * stridec + entry.element; + float result = 0.0f; + if (alpha != 0.0f) { + result = alpha * static_cast(accum[offset]); + } + if (beta != 0.0f) { + result += beta * c[offset]; + } + c[offset] = result; + }); + }); + queue.submit([&](sycl::handler& cgh) { + cgh.depends_on(done); + cgh.host_task([=]() { sycl::free(accum, queue); }); + }); + return done; +} + #define GEMM_STRIDED_BATCH_LAUNCHER_USM(TYPE_A, TYPE_B, TYPE_C, TYPE_S) \ sycl::event gemm_batch(sycl::queue& queue, transpose transa, transpose transb, int64_t m, \ int64_t n, int64_t k, TYPE_S alpha, const TYPE_A* a, int64_t lda, \ @@ -907,6 +1151,16 @@ GEMM_STRIDED_BATCH_LAUNCHER_USM(sycl::half, sycl::half, float, float) #undef GEMM_STRIDED_BATCH_LAUNCHER_USM +sycl::event gemm_batch(sycl::queue& queue, transpose transa, transpose transb, int64_t m, int64_t n, + int64_t k, float alpha, const std::int8_t* a, int64_t lda, int64_t stridea, + const std::int8_t* b, int64_t ldb, int64_t strideb, float beta, float* c, + int64_t ldc, int64_t stridec, int64_t batch_size, + const std::vector& dependencies) { + return gemm_batch_strided_usm_int8_float_impl(queue, transa, transb, m, n, k, alpha, a, lda, + stridea, b, ldb, strideb, beta, c, ldc, stridec, + batch_size, dependencies); +} + #define GEMM_STRIDED_BATCH_LAUNCHER_USM(TYPE_A, TYPE_B, TYPE_C, TYPE_S) \ sycl::event gemm_batch(sycl::queue& queue, transpose transa, transpose transb, int64_t m, \ int64_t n, int64_t k, TYPE_S alpha, const TYPE_A* a, int64_t lda, \ @@ -919,7 +1173,6 @@ GEMM_STRIDED_BATCH_LAUNCHER_USM(sycl::half, sycl::half, float, float) 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 @@ -968,6 +1221,151 @@ inline sycl::event gemm_batch_usm_impl(sycl::queue& queue, transpose* transa, tr return done; } +inline sycl::event gemm_batch_usm_int8_float_impl(sycl::queue& queue, transpose* transa, + transpose* transb, int64_t* m, int64_t* n, + int64_t* k, float* alpha, const std::int8_t** a, + int64_t* lda, const std::int8_t** b, int64_t* ldb, + float* beta, float** c, int64_t* ldc, + int64_t group_count, int64_t* group_size, + const std::vector& dependencies) { + check_int8_float_nonnegative("group_count", group_count); + overflow_check(group_count); + int64_t batch_count = 0; + int64_t accum_size = 0; + int64_t total_entries = 0; + for (int64_t group = 0; group < group_count; ++group) { + check_int8_float_nonnegative("m", m[group]); + check_int8_float_nonnegative("n", n[group]); + check_int8_float_nonnegative("k", k[group]); + check_int8_float_nonnegative("lda", lda[group]); + check_int8_float_nonnegative("ldb", ldb[group]); + check_int8_float_nonnegative("ldc", ldc[group]); + check_int8_float_nonnegative("group_size", group_size[group]); + overflow_check(m[group], n[group], k[group], lda[group], ldb[group], ldc[group], + group_size[group]); + check_int8_float_accumulation_size(k[group]); + if (m[group] > 0 && n[group] > 0 && group_size[group] > 0) { + check_int8_float_output_leading_dimension(m[group], ldc[group]); + } + batch_count = checked_int8_float_sum(batch_count, group_size[group], "total batch count"); + accum_size = + add_int8_float_workspace_elements(accum_size, group_size[group], ldc[group], n[group]); + const int64_t group_entries = + checked_int8_float_product(checked_int8_float_matrix_elements(m[group], n[group]), + group_size[group], "scaling kernel group range"); + total_entries = + checked_int8_float_sum(total_entries, group_entries, "scaling kernel range"); + } + if (total_entries == 0) { + return int8_float_empty_event(queue, dependencies); + } + + // rocBLAS takes an array of pointers for the output of this entry point, so the workspace is one + // allocation split into a matrix per batch. The int32 scalars stand in for oneMath's float alpha + // and beta, which the single scaling kernel below applies instead. + const auto group_count_size = checked_int8_float_size_t(group_count, "group count"); + const auto batch_count_size = checked_int8_float_size_t(batch_count, "total batch count"); + const auto accum_size_size = checked_int8_float_size_t(accum_size, "workspace size"); + const auto total_entries_size = + checked_int8_float_size_t(total_entries, "scaling kernel range"); + auto* alpha_int = sycl::malloc_shared(group_count_size, queue); + auto* beta_int = sycl::malloc_shared(group_count_size, queue); + auto** accum = sycl::malloc_shared(batch_count_size, queue); + auto* metadata = sycl::malloc_shared(group_count_size, queue); + auto* accum_data = + sycl::malloc_device(std::max(accum_size_size, 1), queue); + if (alpha_int == nullptr || beta_int == nullptr || accum == nullptr || metadata == nullptr || + accum_data == nullptr) { + sycl::free(alpha_int, queue); + sycl::free(beta_int, queue); + sycl::free(accum, queue); + sycl::free(metadata, queue); + sycl::free(accum_data, queue); + throw device_bad_alloc("blas", "gemm_batch", queue.get_device()); + } + + int64_t matrix_offset = 0; + int64_t workspace_offset = 0; + int64_t entry_offset = 0; + for (int64_t group = 0; group < group_count; ++group) { + alpha_int[group] = alpha[group] == 0.0f ? 0 : 1; + beta_int[group] = 0; + + const int64_t matrix_size = + checked_int8_float_product(ldc[group], n[group], "workspace matrix size"); + const int64_t group_workspace = + checked_int8_float_product(group_size[group], matrix_size, "workspace group size"); + const int64_t group_entries = + checked_int8_float_product(checked_int8_float_matrix_elements(m[group], n[group]), + group_size[group], "scaling kernel group range"); + const int64_t next_entry_offset = + checked_int8_float_sum(entry_offset, group_entries, "scaling kernel range"); + metadata[group] = { entry_offset, next_entry_offset, matrix_offset, + workspace_offset, m[group], n[group], + ldc[group], alpha[group], beta[group] }; + for (int64_t batch = 0; batch < group_size[group]; ++batch) { + accum[matrix_offset + batch] = accum_data + workspace_offset + batch * matrix_size; + } + matrix_offset = + checked_int8_float_sum(matrix_offset, group_size[group], "total batch count"); + workspace_offset = + checked_int8_float_sum(workspace_offset, group_workspace, "workspace size"); + entry_offset = next_entry_offset; + } + + auto done = gemm_batch_usm_impl(queue, transa, transb, m, n, k, alpha_int, a, lda, b, ldb, + beta_int, accum, ldc, group_count, group_size, dependencies); + + // Locate the group owning each flat entry with a binary search over the immutable metadata. + // This keeps grouped scaling to one kernel submission regardless of group or batch count. + const auto gemm_done = done; + done = queue.submit([&](sycl::handler& cgh) { + cgh.depends_on(gemm_done); + cgh.parallel_for(sycl::range<1>{ total_entries_size }, [=](sycl::id<1> index) { + const int64_t linear = static_cast(index[0]); + int64_t first = 0; + int64_t last = group_count; + while (first < last) { + const int64_t middle = first + (last - first) / 2; + if (linear < metadata[middle].entry_end) { + last = middle; + } + else { + first = middle + 1; + } + } + + const auto group = metadata[first]; + const auto entry = + int8_float_entry_of(static_cast(linear - group.entry_begin), + group.rows, group.columns, group.ld); + const int64_t matrix_size = group.ld * group.columns; + float* output = c[group.matrix_begin + entry.batch]; + const std::int32_t* input = + accum_data + group.workspace_begin + entry.batch * matrix_size; + float result = 0.0f; + if (group.alpha != 0.0f) { + result = group.alpha * static_cast(input[entry.element]); + } + if (group.beta != 0.0f) { + result += group.beta * output[entry.element]; + } + output[entry.element] = result; + }); + }); + queue.submit([&](sycl::handler& cgh) { + cgh.depends_on(done); + cgh.host_task([=]() { + sycl::free(alpha_int, queue); + sycl::free(beta_int, queue); + sycl::free(accum, queue); + sycl::free(metadata, queue); + sycl::free(accum_data, queue); + }); + }); + return done; +} + #define GEMM_BATCH_LAUNCHER_USM(TYPE_A, TYPE_B, TYPE_C, TYPE_S) \ sycl::event gemm_batch(sycl::queue& queue, transpose* transa, transpose* transb, int64_t* m, \ int64_t* n, int64_t* k, TYPE_S* alpha, const TYPE_A** a, int64_t* lda, \ @@ -989,6 +1387,15 @@ GEMM_BATCH_LAUNCHER_USM(sycl::half, sycl::half, float, float) #undef GEMM_BATCH_LAUNCHER_USM +sycl::event gemm_batch(sycl::queue& queue, transpose* transa, transpose* transb, int64_t* m, + int64_t* n, int64_t* k, float* alpha, const std::int8_t** a, int64_t* lda, + const std::int8_t** b, int64_t* ldb, float* beta, float** c, int64_t* ldc, + int64_t group_count, int64_t* group_size, + const std::vector& dependencies) { + return gemm_batch_usm_int8_float_impl(queue, transa, transb, m, n, k, alpha, a, lda, b, ldb, + beta, c, ldc, group_count, group_size, dependencies); +} + #define GEMM_BATCH_LAUNCHER_USM(TYPE_A, TYPE_B, TYPE_C, TYPE_S) \ sycl::event gemm_batch(sycl::queue& queue, transpose* transa, transpose* transb, int64_t* m, \ int64_t* n, int64_t* k, TYPE_S* alpha, const TYPE_A** a, int64_t* lda, \ @@ -1001,7 +1408,6 @@ GEMM_BATCH_LAUNCHER_USM(sycl::half, sycl::half, float, float) 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 @@ -1563,6 +1969,7 @@ GEMM_STRIDED_BATCH_LAUNCHER(std::complex, std::complex, std::com std::complex) 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) #undef GEMM_STRIDED_BATCH_LAUNCHER @@ -1578,7 +1985,6 @@ GEMM_STRIDED_BATCH_LAUNCHER(sycl::half, sycl::half, float, float) 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 @@ -2083,6 +2489,7 @@ GEMM_STRIDED_BATCH_LAUNCHER_USM(std::complex, std::complex, std: std::complex) 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) #undef GEMM_STRIDED_BATCH_LAUNCHER_USM @@ -2098,7 +2505,6 @@ GEMM_STRIDED_BATCH_LAUNCHER_USM(sycl::half, sycl::half, float, float) 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 @@ -2135,6 +2541,7 @@ GEMM_BATCH_LAUNCHER_USM(std::complex, std::complex, std::complex std::complex) 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) #undef GEMM_BATCH_LAUNCHER_USM @@ -2150,7 +2557,6 @@ GEMM_BATCH_LAUNCHER_USM(sycl::half, sycl::half, float, float) 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_usm.cpp b/tests/unit_tests/blas/batch/gemm_batch_stride_usm.cpp index 2cde6dd2a..7b479155a 100644 --- a/tests/unit_tests/blas/batch/gemm_batch_stride_usm.cpp +++ b/tests/unit_tests/blas/batch/gemm_batch_stride_usm.cpp @@ -504,6 +504,101 @@ int int8_accumulation_error_model(device* dev, oneapi::math::layout layout) { return good; } +// The rocBLAS int8-to-float fallback scales its int32 accumulator in a kernel of its own, and +// launches it over a flat range because HIP cannot map every large multi-dimensional range onto +// its grid. A large prime n keeps that flat launch covered. +int int8_flat_range_large_prime(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(); + std::vector dependencies; + const int64_t m = 3, n = 65537, k = 1, batch_size = 1; + const bool col = layout == oneapi::math::layout::col_major; + const int64_t lda = col ? m : k; + const int64_t ldb = col ? k : n; + const int64_t 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 ua = usm_allocator(cxt, *dev); + auto uc = usm_allocator(cxt, *dev); + vector A(stride_a, ua), B(stride_b, ua); + vector C(stride_c, uc); + for (int64_t i = 0; i < m; ++i) + A[col ? i : i * lda] = std::int8_t(i + 1); + for (int64_t j = 0; j < n; ++j) + B[col ? j * ldb : j] = std::int8_t(j % 7 - 3); + std::fill(C.begin(), C.end(), -1.0f); + + try { + event done; +#ifdef CALL_RT_API + if (col) { + done = oneapi::math::blas::column_major::gemm_batch( + main_queue, oneapi::math::transpose::nontrans, oneapi::math::transpose::nontrans, m, + n, k, 1.0f, &A[0], lda, stride_a, &B[0], ldb, stride_b, 0.0f, &C[0], ldc, stride_c, + batch_size, dependencies); + } + else { + done = oneapi::math::blas::row_major::gemm_batch( + main_queue, oneapi::math::transpose::nontrans, oneapi::math::transpose::nontrans, m, + n, k, 1.0f, &A[0], lda, stride_a, &B[0], ldb, stride_b, 0.0f, &C[0], ldc, stride_c, + batch_size, dependencies); + } + done.wait_and_throw(); +#else + if (col) { + TEST_RUN_BLAS_CT_SELECT(main_queue, oneapi::math::blas::column_major::gemm_batch, + oneapi::math::transpose::nontrans, + oneapi::math::transpose::nontrans, m, n, k, 1.0f, &A[0], lda, + stride_a, &B[0], ldb, stride_b, 0.0f, &C[0], ldc, stride_c, + batch_size, dependencies); + } + else { + TEST_RUN_BLAS_CT_SELECT(main_queue, oneapi::math::blas::row_major::gemm_batch, + oneapi::math::transpose::nontrans, + oneapi::math::transpose::nontrans, m, n, k, 1.0f, &A[0], lda, + stride_a, &B[0], ldb, stride_b, 0.0f, &C[0], ldc, stride_c, + batch_size, dependencies); + } + main_queue.wait_and_throw(); +#endif + } + catch (const oneapi::math::unimplemented&) { + return test_skipped; + } + catch (const std::exception& error) { + std::cout << "Error raised during large-prime GEMM_BATCH_STRIDE:\n" + << error.what() << std::endl; + return false; + } + + for (int64_t j = 0; j < n; ++j) + for (int64_t i = 0; i < m; ++i) { + const int64_t index = col ? i + j * ldc : i * ldc + j; + const float expected = float((i + 1) * (j % 7 - 3)); + if (C[index] != expected) { + std::cout << "Difference in entry (" << i << ',' << j << "): DPC++ " << C[index] + << " vs. Reference " << expected << std::endl; + return false; + } + } + return true; +} + class GemmBatchStrideUsmTests : public ::testing::TestWithParam> {}; @@ -527,6 +622,11 @@ TEST_P(GemmBatchStrideUsmTests, Int8Int8SinglePrecisionErrorModel) { (int8_accumulation_error_model(std::get<0>(GetParam()), std::get<1>(GetParam())))); } +TEST_P(GemmBatchStrideUsmTests, Int8Int8SinglePrecisionLargePrimeRange) { + EXPECT_TRUEORSKIP( + (int8_flat_range_large_prime(std::get<0>(GetParam()), std::get<1>(GetParam())))); +} + TEST_P(GemmBatchStrideUsmTests, Int8Int8Int32Precision) { EXPECT_TRUEORSKIP((test( std::get<0>(GetParam()), std::get<1>(GetParam()), 5)));