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/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.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..7b479155a 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,327 @@ 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, 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; +} + +// 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> {}; @@ -283,6 +617,16 @@ 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, 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))); 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 */