[blas][cublas] Support int8 inputs with float output in gemm_batch

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<float>, std::complex<float>, std::complex<float>,
@@ -208,6 +209,8 @@ GEMM_STRIDED_BATCH_LAUNCHER(std::complex<double>, std::complex<double>, 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<TYPE_A, 1>& a, int64_t lda,             \
@@ -220,7 +223,6 @@ GEMM_STRIDED_BATCH_LAUNCHER(std::complex<double>, std::complex<double>, std::com
                                 dtype_string<TYPE_C>() + "," + dtype_string<TYPE_S>() + ">");     \
     }
 
-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<float>, std::complex<float>, std::complex<float>,
@@ -689,7 +692,6 @@ GEMM_STRIDED_BATCH_LAUNCHER_USM(std::complex<double>, std::complex<double>, std:
                                 dtype_string<TYPE_C>() + "," + dtype_string<TYPE_S>() + ">");  \
     }
 
-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<float>, std::complex<float>, std::complex<float>,
@@ -782,7 +785,6 @@ GEMM_BATCH_LAUNCHER_USM(std::complex<double>, std::complex<double>, std::complex
                                 dtype_string<TYPE_C>() + "," + dtype_string<TYPE_S>() + ">");      \
     }
 
-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

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<Tc, int32_t>)
         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<Ta, std::int8_t> &&
+                                   std::is_same_v<Tb, std::int8_t> && std::is_same_v<Tc, float> &&
+                                   std::is_same_v<Ts, float>;
+    double abs_error_bound = 0.0;
+    if constexpr (int8_to_float)
+        abs_error_bound = std::numeric_limits<float>::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;
 }

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<Tc, int32_t>)
         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<Ta, std::int8_t> &&
+                                   std::is_same_v<Tb, std::int8_t> && std::is_same_v<Tc, float> &&
+                                   std::is_same_v<Ts, float>;
+    double abs_error_bound = 0.0;
+    if constexpr (int8_to_float)
+        abs_error_bound = std::numeric_limits<float>::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,232 @@ 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<event> 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<std::int8_t, usm::alloc::shared, 64>(cxt, *dev);
+    auto uc = usm_allocator<float, usm::alloc::shared, 64>(cxt, *dev);
+    vector<std::int8_t, decltype(ua)> A(stride_a * batch_size, ua), B(stride_b * batch_size, ua);
+    vector<float, decltype(uc)> C(stride_c * batch_size, uc);
+    std::vector<float> 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<float>::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, 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;
+
+                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<std::tuple<sycl::device*, oneapi::math::layout>> {};
 
@@ -283,6 +522,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::int8_t, std::int8_t, std::int32_t, float>(
         std::get<0>(GetParam()), std::get<1>(GetParam()), 5)));