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..2cde6dd2a 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,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 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; +} + class GemmBatchStrideUsmTests : public ::testing::TestWithParam> {}; @@ -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::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 */