Test: Add verifyProductUpdate and bound product_sme's update checks absolutely libeigen/eigen!3384 Co-authored-by: Rasmus Munk Larsen <rmlarsen@gmail.com>
diff --git a/test/product_small.cpp b/test/product_small.cpp index 3fe40a3..c152e60 100644 --- a/test/product_small.cpp +++ b/test/product_small.cpp
@@ -446,6 +446,27 @@ VERIFY(!verifyProduct(nan, zero, abs_a, abs_b, 4)); } +// verifyProductUpdate() bounds C0 + alpha*A*B by 5*sqrt(k+2)*eps*2*|| |alpha|*|A|*|B| + |C0| ||_F, here +// 20*eps*(2e-6*|alpha| + 1): an error of a few eps*|C0| passes although the bound on A*B alone rejects it. +template <typename Scalar> +void product_verify_update(const Scalar &large_alpha) { + using Mat = Matrix<Scalar, Dynamic, Dynamic>; + using Real = typename NumTraits<Scalar>::Real; + const Real eps = NumTraits<Scalar>::epsilon(); + const Mat A = Mat::Constant(1, 2, Scalar(1e-3f)), B = Mat::Constant(2, 1, Scalar(1e-3f)), C0 = Mat::Ones(1, 1); + for (const Scalar alpha : {Scalar(-2), large_alpha}) { + const Mat expected = C0 + alpha * A * B; + const Mat off4 = expected.array() + Scalar(4 * eps), off40 = expected.array() + Scalar(40 * eps); + VERIFY(verifyProductUpdate(off4, expected, alpha, A, B, C0)); + VERIFY(!verifyProduct(off4, expected, A, B)); + VERIFY(verifyProductUpdate(off40, expected, alpha, A, B, C0) == (numext::abs(alpha) > Real(2))); + VERIFY(!verifyProductUpdate(Mat::Constant(1, 1, Scalar(std::numeric_limits<Real>::infinity())), expected, alpha, A, + B, C0)); + VERIFY(!verifyProductUpdate(Mat::Constant(1, 1, Scalar(std::numeric_limits<Real>::quiet_NaN())), expected, alpha, A, + B, C0)); + } +} + EIGEN_DECLARE_TEST(product_small) { for (int i = 0; i < g_repeat; i++) { CALL_SUBTEST_1(product(Matrix<float, 3, 2>())); @@ -492,6 +513,8 @@ CALL_SUBTEST_6(product_small_regressions<0>()); CALL_SUBTEST_6(product_transposed_triple_bound<0>()); CALL_SUBTEST_6(product_verify_rejects_nonfinite<0>()); + CALL_SUBTEST_6(product_verify_update(1e6f)); + CALL_SUBTEST_6(product_verify_update(std::complex<float>(0.0f, 1e6f))); // Deterministic sweep at transition boundaries (outside g_repeat). CALL_SUBTEST_54(product_transition_sizes<float>());
diff --git a/test/product_sme.cpp b/test/product_sme.cpp index 7a43822..d2c91da 100644 --- a/test/product_sme.cpp +++ b/test/product_sme.cpp
@@ -279,19 +279,19 @@ const SmeColMajorMat<Scalar> reference = deepA.lazyProduct(deepB); C = initial; C.noalias() += deepA * deepB; - VERIFY_IS_APPROX(SmeColMajorMat<Scalar>(C), SmeColMajorMat<Scalar>(initial + reference)); + VERIFY(verifyProductUpdate(C, initial + reference, Scalar(1), deepA, deepB, initial)); C = initial; C.noalias() += alpha * (deepA * deepB); - VERIFY_IS_APPROX(SmeColMajorMat<Scalar>(C), SmeColMajorMat<Scalar>(initial + alpha * reference)); + VERIFY(verifyProductUpdate(C, initial + alpha * reference, alpha, deepA, deepB, initial)); // A real alpha takes the complex store's fused path: -1 from -=, and a // non-unit real scale. const Scalar real_alpha(typename NumTraits<Scalar>::Real(1.375)); C = initial; C.noalias() -= deepA * deepB; - VERIFY_IS_APPROX(SmeColMajorMat<Scalar>(C), SmeColMajorMat<Scalar>(initial - reference)); + VERIFY(verifyProductUpdate(C, initial - reference, Scalar(-1), deepA, deepB, initial)); C = initial; C.noalias() += real_alpha * (deepA * deepB); - VERIFY_IS_APPROX(SmeColMajorMat<Scalar>(C), SmeColMajorMat<Scalar>(initial + real_alpha * reference)); + VERIFY(verifyProductUpdate(C, initial + real_alpha * reference, real_alpha, deepA, deepB, initial)); #ifdef EIGEN_VECTORIZE_SME_F64F64 // RowMajor GEMM normally transposes the computation. Call the packed kernel @@ -309,7 +309,7 @@ internal::sme_gebp_dispatch<Scalar, false, false>::run( C.data(), Index(ResMat::IsRowMajor ? n : 1), Index(ResMat::IsRowMajor ? 1 : n), packedA.data(), packedB.data(), Index(n), Index(65), Index(n), scale, Index(65), Index(65), Index(0), Index(0)); - VERIFY_IS_APPROX(SmeColMajorMat<Scalar>(C), SmeColMajorMat<Scalar>(initial + scale * reference)); + VERIFY(verifyProductUpdate(C, initial + scale * reference, scale, deepA, deepB, initial)); } } #endif @@ -510,16 +510,16 @@ MC C = MC::Random(m, n); const MC C0 = C; C.noalias() = A * B; - VERIFY_IS_APPROX(C, ref); + VERIFY(verifyProduct(C, ref, A, B)); C = C0; C.noalias() += A * B; - VERIFY_IS_APPROX(C, C0 + ref); + VERIFY(verifyProductUpdate(C, C0 + ref, Scalar(1), A, B, C0)); C = C0; C.noalias() -= Scalar(2) * A * B; - VERIFY_IS_APPROX(C, C0 - Scalar(2) * ref); + VERIFY(verifyProductUpdate(C, C0 - Scalar(2) * ref, Scalar(-2), A, B, C0)); const MA At = A.transpose(); C.noalias() = At.transpose() * B; - VERIFY_IS_APPROX(C, ref); + VERIFY(verifyProduct(C, ref, A, B)); } } @@ -776,7 +776,7 @@ const SmeColMajorMat<Scalar> Ac = ConjLhs ? A.conjugate().eval() : A; const SmeColMajorMat<Scalar> Bc = ConjRhs ? B.conjugate().eval() : B; const SmeColMajorMat<Scalar> ref = c_before + alpha * Ac.lazyProduct(Bc); - VERIFY_IS_APPROX(SmeColMajorMat<Scalar>(C), ref); + VERIFY(verifyProductUpdate(C, ref, alpha, Ac, Bc, c_before)); // Cells the strided map skips stay untouched. if (strided_c) {
diff --git a/test/product_test_helpers.h b/test/product_test_helpers.h index 746de37..55f63f2 100644 --- a/test/product_test_helpers.h +++ b/test/product_test_helpers.h
@@ -33,40 +33,63 @@ // lambda: probability parameter; P(lambda) = 1 - 2*exp(-lambda^2/2). // lambda=5 gives P > 0.9999 per inner product. -// Returns lambda * sqrt(k) * epsilon * num_products * || |A|*|B| ||_F, with epsilon that of the -// operand scalar but the value kept in double: narrowed to half, a triple-product bound already -// overflows at k = 100 for all-one operands. +// Returns || |alpha|*|A|*|B| + |C0| ||_F, with |C0| given coefficient-wise by c0(i, j). // -// || |A|*|B| ||_F is accumulated in double with an explicit loop rather than formed as +// The norm is accumulated in double with an explicit loop rather than formed as // (A.cwiseAbs() * B.cwiseAbs()).norm(): a bound built with the product implementation would inherit // that implementation's defects, and one that overflowed to infinity would admit any result at all. -template <typename DerivedA, typename DerivedB> -double product_error_bound(const MatrixBase<DerivedA>& A, const MatrixBase<DerivedB>& B, int num_products = 1, - double lambda = 5) { - using Scalar = typename DerivedA::Scalar; - const Index k = A.cols(); +template <typename DerivedA, typename DerivedB, typename C0Coeff> +double abs_product_update_norm(double abs_alpha, const MatrixBase<DerivedA>& A, const MatrixBase<DerivedB>& B, + C0Coeff c0) { double squared_norm = 0.0; for (Index i = 0; i < A.rows(); ++i) for (Index j = 0; j < B.cols(); ++j) { double sum = 0.0; - for (Index l = 0; l < k; ++l) + for (Index l = 0; l < A.cols(); ++l) sum += static_cast<double>(numext::abs(A.coeff(i, l))) * static_cast<double>(numext::abs(B.coeff(l, j))); + sum = abs_alpha * sum + c0(i, j); squared_norm += sum * sum; } - return lambda * numext::sqrt(double(k)) * static_cast<double>(NumTraits<Scalar>::epsilon()) * double(num_products) * - numext::sqrt(squared_norm); + return numext::sqrt(squared_norm); } -// Verify that two computations of A*B agree within the Higham-Mary bound. -// Returns true if ||actual - expected||_F <= product_error_bound(A, B, ...). +// Returns lambda * sqrt(k) * epsilon * num_products * || |A|*|B| ||_F, with epsilon that of the +// operand scalar but the value kept in double: narrowed to half, a triple-product bound already +// overflows at k = 100 for all-one operands. +template <typename DerivedA, typename DerivedB> +double product_error_bound(const MatrixBase<DerivedA>& A, const MatrixBase<DerivedB>& B, int num_products = 1, + double lambda = 5) { + using Scalar = typename DerivedA::Scalar; + const double norm = abs_product_update_norm(1.0, A, B, [](Index, Index) { return 0.0; }); + return lambda * numext::sqrt(double(A.cols())) * static_cast<double>(NumTraits<Scalar>::epsilon()) * + double(num_products) * norm; +} + +// Returns, for C = C0 + alpha*A*B (for beta*C0 + alpha*A*B, pass beta*C0 as C0), +// lambda * sqrt(k+2) * epsilon * num_products * || |alpha|*|A|*|B| + |C0| ||_F. +// k+2: a k-term inner product, then the alpha scale and the sum with c0 (real arithmetic; complex +// scalars use the same depth by extension, as in product_error_bound). The |C0| term is needed: a +// bound on A*B alone fails once |alpha|*|A|*|B| << |C0|. +template <typename DerivedA, typename DerivedB, typename DerivedC0> +double product_update_error_bound(const typename DerivedA::Scalar& alpha, const MatrixBase<DerivedA>& A, + const MatrixBase<DerivedB>& B, const MatrixBase<DerivedC0>& C0, int num_products = 1, + double lambda = 5) { + using Scalar = typename DerivedA::Scalar; + eigen_assert(C0.rows() == A.rows() && C0.cols() == B.cols()); + const double norm = abs_product_update_norm(static_cast<double>(numext::abs(alpha)), A, B, [&C0](Index i, Index j) { + return static_cast<double>(numext::abs(C0.coeff(i, j))); + }); + return lambda * numext::sqrt(double(A.cols() + 2)) * static_cast<double>(NumTraits<Scalar>::epsilon()) * + double(num_products) * norm; +} + +// Returns true if ||actual - expected||_F <= bound. // // The error norm is accumulated in double like the bound, so that it overflows only when the // difference itself does. A non-finite error or bound fails: inf <= inf would otherwise accept // an overflowed result. -template <typename D1, typename D2, typename DA, typename DB> -inline bool verifyProduct(const MatrixBase<D1>& actual, const MatrixBase<D2>& expected, const MatrixBase<DA>& A, - const MatrixBase<DB>& B, int num_products = 2, double lambda = 5) { - const double bound = product_error_bound(A, B, num_products, lambda); +template <typename D1, typename D2> +bool verify_within_product_bound(const MatrixBase<D1>& actual, const MatrixBase<D2>& expected, double bound) { const typename D1::PlainObject diff = actual - expected; double squared_error = 0.0; for (Index j = 0; j < diff.cols(); ++j) @@ -82,6 +105,24 @@ return true; } +// Verify that two computations of A*B agree within the Higham-Mary bound. +// Returns true if ||actual - expected||_F <= product_error_bound(A, B, ...). +template <typename D1, typename D2, typename DA, typename DB> +inline bool verifyProduct(const MatrixBase<D1>& actual, const MatrixBase<D2>& expected, const MatrixBase<DA>& A, + const MatrixBase<DB>& B, int num_products = 2, double lambda = 5) { + return verify_within_product_bound(actual, expected, product_error_bound(A, B, num_products, lambda)); +} + +// Verify that two computations of C0 + alpha*A*B agree within the Higham-Mary bound. +// Returns true if ||actual - expected||_F <= product_update_error_bound(alpha, A, B, C0, ...). +template <typename D1, typename D2, typename DA, typename DB, typename DC0> +bool verifyProductUpdate(const MatrixBase<D1>& actual, const MatrixBase<D2>& expected, const typename DA::Scalar& alpha, + const MatrixBase<DA>& A, const MatrixBase<DB>& B, const MatrixBase<DC0>& C0, + int num_products = 2, double lambda = 5) { + return verify_within_product_bound(actual, expected, + product_update_error_bound(alpha, A, B, C0, num_products, lambda)); +} + } // namespace Eigen #endif // EIGEN_TEST_PRODUCT_TEST_HELPERS_H