From 7f67361ad738622b9a50beef105e4a4f7ab21736 Mon Sep 17 00:00:00 2001 From: Claude Date: Mon, 14 Sep 2026 23:37:12 +0000 Subject: [PATCH 1/6] tests: sweep the small-dimension cross product in every portable suite LAPACK's linear-equation suites take M and N from a list that starts 0 1 2 3 5 and run every path over the full M x N cross product, so the empty operand, the single row or column, and the orders around a block boundary are exercised in every routine. The portable suites only ran them where someone had thought of them (gels had m = 0 and 1 x n, geqrf and potrf started at 3, ormqr ran one fixed m). Every portable suite now sweeps its run_case over small_dims = {0, 1, 2, 3, 5} (test_compact_util.hpp): geqrf, trsm and gels over m x n (gels per layout and trans), ormqr over m x nrhs, orgqr over m >= n with LAPACK's K values {n, 0, 1, n/2}, potrf and sytrfnp over n through the factorization checks and n x nrhs through the solve pipeline -- each with one full group plus a padded one, and the existing checks unchanged, so a zero extent is checked as the routine's quick return (info = 0, every error 0). What that needed: - run_case passes ld = max(1, extent), as LAPACK's ld >= max(1, m) requires and the C API validates. - compact_buffer never returns an empty vector: a zero-extent operand is still an argument the routine may require present (posv's bp at nrhs = 0, gels's A at m = 0, asserted in cbk.cpp), so it keeps the one element a Fortran caller's dummy would be. - is_rowmajor accepts an empty view (mat_view of a 0 x n matrix has sj = 0), leaving ld >= max(1, extent) to lapack_ld. - The LAPACKE-backed references take LAPACK's own quick return on an empty operand before the call, so an empty batch's possibly null pointer never reaches a stack's NaN check; ref_gels's is B := 0, the ?gels contract, and the gels suite's residual-rows identity is skipped at min(m, n) = 0 for the same reason. Passes gcc Release, clang Debug (the view asserts on the empty shapes) and the MKL-extension build. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_014ZZ76AgYwQ3vtuyZ46tyMX --- tests/test_compact_util.hpp | 20 ++++++++++++++++++- tests/test_gels_compact.cpp | 26 ++++++++++++++++++++----- tests/test_geqrf_compact.cpp | 27 +++++++++++++++++++------- tests/test_lapack_util.hpp | 32 ++++++++++++++++++++++++++++--- tests/test_orgqr_compact.cpp | 20 +++++++++++++++++-- tests/test_ormqr_compact.cpp | 24 +++++++++++++++++------ tests/test_potrf_compact.cpp | 34 ++++++++++++++++++++++++--------- tests/test_sytrfnp_compact.cpp | 35 ++++++++++++++++++++++++---------- tests/test_trsm_compact.cpp | 24 +++++++++++++++++------ 9 files changed, 193 insertions(+), 49 deletions(-) diff --git a/tests/test_compact_util.hpp b/tests/test_compact_util.hpp index 4e30c39..452dc14 100644 --- a/tests/test_compact_util.hpp +++ b/tests/test_compact_util.hpp @@ -343,12 +343,18 @@ template void gen_tri(MatrixView A, bool upper) // dimension ldp -- the one place the ng * gstride sizing is written. For a // buffer the kernel fills (tau, say), use it directly; for one packed from a // batch, the pack_compact/pack_tau overloads below return it filled. +// +// Never empty: an operand with a zero extent is still an argument the routine +// may require present (Fortran semantics -- ?posv's bp at nrhs = 0, ?gels's A +// at m = 0 -- asserted in cbk.cpp), so the buffer keeps one element, the +// dummy a Fortran caller would pass. template std::vector compact_buffer(int nm, int rows, int cols, int ldp, int V, bool rowmajor = false) { const int ng = (nm + V - 1) / V; - return std::vector((std::size_t)ng * group_stride(rowmajor, ldp, rows, cols, V)); + return std::vector(std::max( + 1, (std::size_t)ng * group_stride(rowmajor, ldp, rows, cols, V))); } template @@ -438,6 +444,18 @@ template void unpack_tau(MatrixBatch &tau, const T *tp, int V) unpack_compact(tau, tp, tau.rows(), V); } +// ----------------------- the small-dimension sweep ------------------- +// LAPACK's linear-equation suites take M and N from a list that starts +// 0 1 2 3 5 (TESTING/dtest.in) and run every path over the full M x N cross +// product, so the empty operand, the single row or column, and the orders +// just below and above a register block are exercised in every routine, not +// only where someone thought of them. The suites' sweeps take their extents +// from this list and hand each shape to their own run_case: a zero extent is +// the routine's quick return, and the case checks exactly that (info = 0 and +// every error 0). The leading dimension of an empty operand is 1, as LAPACK's +// ld >= max(1, m) requires; the sweeps spell it max(1, extent). +constexpr int small_dims[] = {0, 1, 2, 3, 5}; + // ----------------------- shared checks and epilogues ----------------- // One row of a C API validation table: the call made, what it returned, and diff --git a/tests/test_gels_compact.cpp b/tests/test_gels_compact.cpp index bb8c90b..d0ded44 100644 --- a/tests/test_gels_compact.cpp +++ b/tests/test_gels_compact.cpp @@ -17,8 +17,10 @@ // residual sums of squares in rows n..m-1 of B equal ||B - op(A) X||^2 // minimum norm: op(A) X = B, and X equals the minimum-norm solution formed // the other way, X = op(A)^T Z with (op(A) op(A)^T) Z = B -// plus LAPACK-style argument validation of the C API and the min(m,n) = 0 case -// (B := 0). +// over the suite's shapes and the small-dimension cross product (small_dims: +// m x n from 0 to 5 per (layout, trans), so every over-/underdetermined case +// meets the empty operand and the single row or column), plus LAPACK-style +// argument validation of the C API and the min(m,n) = 0 case (B := 0). // // Assisted-by: Claude:claude-fable-5 @@ -63,7 +65,7 @@ int run_case(char layout, char trans, int nm, int m, int n, int nrhs) } // pack in the requested layout, solve, unpack - const int lda = row ? n : m, ldb = row ? nrhs : p; + const int lda = std::max(1, row ? n : m), ldb = std::max(1, row ? nrhs : p); std::vector ap = pack_compact(A, lda, V, row); std::vector bp = pack_compact(B, ldb, V, row); std::vector tp = compact_buffer(nm, q, 1, q, V); /* the kernel fills it */ @@ -111,8 +113,10 @@ int run_case(char layout, char trans, int nm, int m, int n, int nrhs) e_prop = std::max(e_prop, std::abs((double)Ss[e]) / scale); // rows cols_op..rows_op-1 of B hold the residual: squared column // norms == ||r_j||^2 (relative to ||b_j||^2: a square system has no - // residual rows and a zero residual) - for (int j = 0; j < nrhs; ++j) { + // residual rows and a zero residual). Not at min(m,n) = 0: that is + // ?gels's quick return, B := 0 (checked against ref_gels above), + // and the residual is B itself. + for (int j = 0; q > 0 && j < nrhs; ++j) { double rss = 0, rss_b = 0, bb = 0; for (int i = 0; i < rows_op; ++i) { rss += (double)R(i, j) * R(i, j); @@ -265,5 +269,17 @@ int main() fails += run_case('C', 'N', 8, 64, 20, 6); // wider RHS blocks (4+2 tail) fails += run_case('C', 'T', 8, 64, 20, 7); // (4+2+1 tail) + // the small-dimension cross product (small_dims): m x n per (layout, + // trans), so every over-/underdetermined case meets the empty operand, a + // single row or column, and the orders around the register block; a + // padded group + for (int mm : small_dims) + for (int nn : small_dims) { + for (char lay : {'C', 'R'}) + for (char tr : {'N', 'T'}) + fails += run_case(lay, tr, 5, mm, nn, 2); + fails += run_case('R', 'T', 9, mm, nn, 1); + } + return finish(fails); } diff --git a/tests/test_geqrf_compact.cpp b/tests/test_geqrf_compact.cpp index 07b15c6..429bac8 100644 --- a/tests/test_geqrf_compact.cpp +++ b/tests/test_geqrf_compact.cpp @@ -17,7 +17,10 @@ // ||Q^T Q - I||_1 <= 100 n eps, on dense, column-scaled, rank-deficient // and near-collinear inputs; (H, tau) elementwise vs the blocked // LAPACKE_?geqrf on the dense ones (relative, 100 n eps) -// plus LAPACK-style argument validation of the C API. +// over the suite's shapes and the small-dimension cross product (small_dims: +// m x n from 0 to 5, so the empty operand and the single row or column run +// through the same checks), plus LAPACK-style argument validation of the C +// API. // // Assisted-by: Claude:claude-opus-4-8 Claude:claude-fable-5 @@ -49,13 +52,14 @@ template static int run_case(int nm, int m, int n) } // pack A, factor with the routine under test, unpack (H, tau) - std::vector ap = pack_compact(A, m, V); + const int lda = std::max(1, m); + std::vector ap = pack_compact(A, lda, V); std::vector tp = compact_buffer(nm, k, 1, k, V); /* the kernel fills it */ - int info = compact::geqrf('C', m, n, ap.data(), m, tp.data(), V, nm); + int info = compact::geqrf('C', m, n, ap.data(), lda, tp.data(), V, nm); MatrixBatch Aout(nm, m, n), tau_out(nm, k, 1); - unpack_compact(Aout, ap.data(), m, V); + unpack_compact(Aout, ap.data(), lda, V); unpack_tau(tau_out, tp.data(), V); // check 1: (H, tau) match the scalar reference elementwise @@ -87,10 +91,11 @@ template static int run_case(int nm, int m, int n) MatrixBatch B(nm, n, nrhs); for (int idx = 0; idx < nm; ++idx) matmul(A.view(idx), X, B.view(idx)); /* B = A X */ - std::vector bp = pack_compact(B, n, V); - compact::ormqr('T', n, nrhs, k, ap.data(), n, tp.data(), bp.data(), n, V, nm); + std::vector bp = pack_compact(B, lda, V); + compact::ormqr('T', n, nrhs, k, ap.data(), lda, tp.data(), bp.data(), lda, V, + nm); MatrixBatch Bo(nm, n, nrhs); - unpack_compact(Bo, bp.data(), n, V); + unpack_compact(Bo, bp.data(), lda, V); e_solve = 0; for (int idx = 0; idx < nm; ++idx) { ref_trsm_upper(Aout.view(idx), Bo.view(idx)); @@ -355,6 +360,14 @@ int main() // float fails += run_case(16, 30, 30); fails += run_case(32, 43, 17); + // the small-dimension cross product (small_dims): the empty operand, a + // single row or column, and the orders around the register block, one + // full group plus a padded one + for (int m : small_dims) + for (int n : small_dims) { + fails += run_case(5, m, n); + fails += run_case(9, m, n); + } // numerical scope: an underflowing column reads as already triangular fails += test_underflow(); fails += test_underflow(); diff --git a/tests/test_lapack_util.hpp b/tests/test_lapack_util.hpp index 029c27c..2d9f233 100644 --- a/tests/test_lapack_util.hpp +++ b/tests/test_lapack_util.hpp @@ -115,12 +115,15 @@ CBK_TEST_LAPACK_DISPATCH(float, s) // without a copy when one stride is 1 and the other is at least the extent // along it -- column-major (si == 1, ld = sj) or row-major (sj == 1, ld = si); // the ambiguous degenerate shapes (a single row or column) are read as -// whichever description satisfies the library's ld >= extent rule. +// whichever description satisfies the library's ld >= extent rule. An empty +// extent constrains nothing (mat_view of a 0 x n matrix has sj = 0); the +// library's ld >= max(1, extent) is then lapack_ld's business. template bool is_rowmajor(const Mv &M) { - if (M.si == 1 && M.sj >= std::max(M.rows, 1)) return false; - assert(M.sj == 1 && M.si >= std::max(M.cols, 1) && "view is not a BLAS layout"); + if (M.si == 1 && (M.rows == 0 || M.sj >= std::max(M.rows, 1))) return false; + assert(M.sj == 1 && (M.cols == 0 || M.si >= std::max(M.cols, 1)) && + "view is not a BLAS layout"); return true; } @@ -206,16 +209,28 @@ template class Staged { // sequence: the suites gate such comparisons relative to the operand norms, // at a multiple of n * eps. Only ?geqr2 is unblocked by definition, and only // it is compared elementwise at ~eps. +// +// Each takes LAPACK's own quick return on an empty operand before the call: +// the routine would do nothing, but LAPACKE's NaN check runs first and an +// empty batch's operand is a possibly null pointer (MatrixBatch hands out +// std::vector storage), which not every stack's check guards. + +template bool is_empty(const Mv &M) +{ + return M.rows == 0 || M.cols == 0; +} // ?geqr2: unblocked Householder QR, (H, tau) in the LAPACK convention. template void ref_geqr2(MatrixView A, T *tau) { + if (is_empty(A)) return; lapack::geqr2(lapack_layout(A), A.rows, A.cols, A.data, lapack_ld(A), tau); } // ?geqrf: the blocked driver, for the invariants suites. template void ref_geqrf(MatrixView A, T *tau) { + if (is_empty(A)) return; lapack::geqrf(lapack_layout(A), A.rows, A.cols, A.data, lapack_ld(A), tau); } @@ -225,6 +240,7 @@ template void ref_ormqr(char trans, int k, Av A, const T *tau, MatrixView B) { assert(A.rows == B.rows && k <= std::min(A.rows, A.cols)); + if (is_empty(B) || k == 0) return; const Staged Bs(B, is_rowmajor(A)); lapack::ormqr(lapack_layout(A), 'L', trans, B.rows, B.cols, k, A.data, lapack_ld(A), tau, Bs.data(), Bs.ld()); @@ -235,6 +251,7 @@ void ref_ormqr(char trans, int k, Av A, const T *tau, MatrixView B) template void ref_orgqr(int k, MatrixView A, const T *tau) { assert(k <= A.cols && A.cols <= A.rows); + if (is_empty(A)) return; /* k = 0 is not empty: Q = I(:, 0:n-1) is written */ lapack::orgqr(lapack_layout(A), A.rows, A.cols, k, A.data, lapack_ld(A), tau); } @@ -257,6 +274,7 @@ template void ref_geqp3(MatrixView A, lapack_int *jpvt, T *tau) template lapack_int ref_potrf(char uplo, MatrixView A) { assert(A.rows == A.cols); + if (is_empty(A)) return 0; return lapack::potrf(lapack_layout(A), uplo, A.rows, A.data, lapack_ld(A)); } @@ -306,10 +324,18 @@ template void ref_trsm_upper(Rv R, MatrixView B) // ?gelqf storage (m < n) -- and tau its min(m,n) reflector scalars. LAPACK's // ?gels keeps tau in its workspace, so it runs on a copy of A for X, and // (A, tau) come from ?geqrf / ?gelqf on A itself, the same call ?gels makes. +// min(m, n) = 0 is ?gels's quick return, B := 0 over all max(m, n) rows (no +// factorization, no least-squares solution: the residual rows hold zeros too). template void ref_gels(char trans, MatrixView A, MatrixView B, T *tau) { const int m = A.rows, n = A.cols, nrhs = B.cols; assert(B.rows == std::max(m, n)); + if (is_empty(A)) { + for (int j = 0; j < nrhs; ++j) + for (int i = 0; i < B.rows; ++i) + B(i, j) = T(0); + return; + } const bool row = is_rowmajor(A); std::vector A0s((std::size_t)m * n); diff --git a/tests/test_orgqr_compact.cpp b/tests/test_orgqr_compact.cpp index db40bf6..6d09867 100644 --- a/tests/test_orgqr_compact.cpp +++ b/tests/test_orgqr_compact.cpp @@ -18,7 +18,9 @@ * 3. Q(:, 0:k-1) R = A with R = triu of the factorization (k > 0) * 4. padded lanes of a partial last group come out exactly identity * plus k < n (unit-seeded extra columns), k = 0 (Q = leading columns of I), - * row-major, and the C API's LAPACK-style argument validation. + * row-major, the small-dimension sweep (small_dims: every m >= n from 0 to 5 + * with k in {n, 0, 1, n/2}, LAPACK's K values, the empty operand included), + * and the C API's LAPACK-style argument validation. * * Assisted-by: Claude:claude-fable-5 */ @@ -46,7 +48,7 @@ static int run_case(int nm, int m, int n, int k, bool rowmajor = false) const double tol_orth = 100.0 * eps * m; /* independent invariants */ const double tol_rec = 100.0 * eps * m; - const int ld = rowmajor ? n : m; + const int ld = std::max(1, rowmajor ? n : m); /* factor a random m x k batch (empty at k = 0); stage its reflectors * into m x n */ @@ -192,6 +194,20 @@ int main() /* smallest size */ fails += run_case(3, 1, 1, 1); + /* the small-dimension sweep (small_dims): every m >= n, k over LAPACK's + * {n, 0, 1, n/2} (each value once), a padded group, one layout each */ + for (int mm : small_dims) + for (int nn : small_dims) { + if (nn > mm) continue; + const int ks[] = {nn, 0, 1, nn / 2}; + for (int t = 0; t < 4; ++t) { + const int kk = ks[t]; + if (kk > nn || std::find(ks, ks + t, kk) != ks + t) continue; + fails += run_case(5, mm, nn, kk); + fails += run_case(9, mm, nn, kk, true); + } + } + /* row-major */ fails += run_case(8, 43, 43, 43, true); fails += run_case(11, 43, 20, 20, true); /* padded partial group */ diff --git a/tests/test_ormqr_compact.cpp b/tests/test_ormqr_compact.cpp index 7324067..4fb5ea2 100644 --- a/tests/test_ormqr_compact.cpp +++ b/tests/test_ormqr_compact.cpp @@ -14,6 +14,9 @@ * the same reflectors, LAPACK's blocked above its crossover) * 2. back substitution recovers X * 3. applying 'N' after 'T' recovers the original B (Q Q^T = I) + * over the suite's shapes and the small-dimension cross product (small_dims: + * m x nrhs from 0 to 5, the empty operand and the 1/2/3-column right-hand + * side blocks included). * * Assisted-by: Claude:claude-fable-5 Claude:claude-opus-4-8 */ @@ -56,13 +59,14 @@ template static int run_case(int nm, int m, int nrhs) ref_ormqr('T', k, Afac.view(kk), tau[kk], Bref.view(kk)); } - std::vector ap = pack_compact(Afac, m, V); + const int ld = std::max(1, m); + std::vector ap = pack_compact(Afac, ld, V); std::vector tp = pack_tau(tau, V); - std::vector bp = pack_compact(B, m, V); + std::vector bp = pack_compact(B, ld, V); /* check 1: compact Q^T B vs LAPACKE_?ormqr, relative to ||Q^T B||_1 */ - compact::ormqr('T', m, nrhs, k, ap.data(), m, tp.data(), bp.data(), m, V, nm); - unpack_compact(Bout, bp.data(), m, V); + compact::ormqr('T', m, nrhs, k, ap.data(), ld, tp.data(), bp.data(), ld, V, nm); + unpack_compact(Bout, bp.data(), ld, V); double e1 = 0; for (int kk = 0; kk < nm; ++kk) e1 = std::max(e1, max_abs_diff(Bout[kk], Bref[kk], (size_t)m * nrhs) / @@ -76,8 +80,8 @@ template static int run_case(int nm, int m, int nrhs) } /* check 3: 'N' undoes 'T' */ - compact::ormqr('N', m, nrhs, k, ap.data(), m, tp.data(), bp.data(), m, V, nm); - unpack_compact(Bout, bp.data(), m, V); + compact::ormqr('N', m, nrhs, k, ap.data(), ld, tp.data(), bp.data(), ld, V, nm); + unpack_compact(Bout, bp.data(), ld, V); double e3 = 0; for (int kk = 0; kk < nm; ++kk) e3 = std::max(e3, max_abs_diff(Bout[kk], B[kk], (size_t)m * nrhs)); @@ -231,6 +235,14 @@ int main(int argc, char **) fails += run_case(16, m, nrhs); fails += run_case(32, m, nrhs); + /* the small-dimension cross product (small_dims): m x nrhs, the empty + * operand and the 1/2/3-column RHS blocks included, a padded group */ + for (int mm : small_dims) + for (int nr : small_dims) { + fails += run_case(5, mm, nr); + fails += run_case(9, mm, nr); + } + /* column-pivoted QR + back-permutation solve (salvaged cross-check) */ fails += run_case_pivoted(8, m, nrhs); fails += run_case_pivoted(11, m, nrhs); /* padded partial group */ diff --git a/tests/test_potrf_compact.cpp b/tests/test_potrf_compact.cpp index d505fa7..03dc448 100644 --- a/tests/test_potrf_compact.cpp +++ b/tests/test_potrf_compact.cpp @@ -18,7 +18,10 @@ // ?posv_compact reproduces that factor and X bit-for-bit // 5. nrhs = 0: ?posv_compact still factors (LAPACK ?posv), bit-identical to // ?potrf_compact, with a 1-element dummy bp (never referenced) -// plus LAPACK-style argument validation of the three C APIs. +// over the suite's orders and the small-dimension sweep (small_dims: every +// order from 0 to 5 through the factorization checks, n x nrhs through the +// solve pipeline, the empty operand included), plus LAPACK-style argument +// validation of the three C APIs. // // Assisted-by: Claude:claude-opus-4-8 Claude:claude-fable-5 @@ -51,10 +54,11 @@ static int run_case(int nm, int n, char uplo, char layout, double cond = 0.0) } // pack the full symmetric A, factor with the routine under test, unpack - std::vector ap = pack_compact(A, n, V, rowmajor); - int info = compact::potrf(layout, uplo, n, ap.data(), n, V, nm); + const int lda = std::max(1, n); + std::vector ap = pack_compact(A, lda, V, rowmajor); + int info = compact::potrf(layout, uplo, n, ap.data(), lda, V, nm); MatrixBatch Aout(nm, n, n); - unpack_compact(Aout, ap.data(), n, V, rowmajor); + unpack_compact(Aout, ap.data(), lda, V, rowmajor); double e_fac = 0, e_rec = 0, e_untouched = 0; std::vector Fs((size_t)n * n), Recs((size_t)n * n); @@ -116,7 +120,7 @@ static int run_solve(int nm, int n, int nrhs, char uplo, char layout) { const T eps = std::numeric_limits::epsilon(); const bool rowmajor = (layout == 'R' || layout == 'r'); - const int ldb = rowmajor ? nrhs : n; + const int lda = std::max(1, n), ldb = std::max(1, rowmajor ? nrhs : n); // known X, B = A X densely MatrixBatch A(nm, n, n), B(nm, n, nrhs); @@ -127,15 +131,15 @@ static int run_solve(int nm, int n, int nrhs, char uplo, char layout) matmul(A.view(idx), X, B.view(idx)); } - std::vector ap = pack_compact(A, n, V, rowmajor); + std::vector ap = pack_compact(A, lda, V, rowmajor); std::vector bp = pack_compact(B, ldb, V, rowmajor); std::vector ap2 = ap, bp2 = bp; // the fused call's copies - int info_f = compact::potrf(layout, uplo, n, ap.data(), n, V, nm); + int info_f = compact::potrf(layout, uplo, n, ap.data(), lda, V, nm); int info_s = - compact::potrs(layout, uplo, n, nrhs, ap.data(), n, bp.data(), ldb, V, nm); + compact::potrs(layout, uplo, n, nrhs, ap.data(), lda, bp.data(), ldb, V, nm); int info_v = - compact::posv(layout, uplo, n, nrhs, ap2.data(), n, bp2.data(), ldb, V, nm); + compact::posv(layout, uplo, n, nrhs, ap2.data(), lda, bp2.data(), ldb, V, nm); MatrixBatch Xhat(nm, n, nrhs); unpack_compact(Xhat, bp.data(), ldb, V, rowmajor); @@ -306,6 +310,18 @@ int main() fails += run_solve(16, 24, 3, u, l); } + // The small-dimension sweep (small_dims): every order over uplo x layout, + // then n x nrhs through the solve pipeline, a padded group each. + for (int n : small_dims) + for (char u : uplos) + for (char l : lays) + fails += run_case(5, n, u, l); + for (int n : small_dims) + for (int nrhs : small_dims) { + fails += run_solve(5, n, nrhs, 'L', 'C'); + fails += run_solve(9, n, nrhs, 'U', 'R'); + } + // nrhs = 0 must factor anyway (LAPACK ?posv), bit-identical to potrf, // bp a never-referenced dummy. fails += run_nrhs0(6, 20, 'L', 'C'); diff --git a/tests/test_sytrfnp_compact.cpp b/tests/test_sytrfnp_compact.cpp index d05bdd4..feb8874 100644 --- a/tests/test_sytrfnp_compact.cpp +++ b/tests/test_sytrfnp_compact.cpp @@ -20,7 +20,9 @@ // ?sysvnp_compact reproduces that factor and X bit-for-bit // 5. nrhs = 0: ?sysvnp_compact still factors (LAPACK ?sysv), bit-identical // to ?sytrfnp_compact, with a 1-element dummy bp (never referenced) -// plus: +// over the suite's orders and the small-dimension sweep (small_dims: every +// order from 0 to 5 through the factorization checks, n x nrhs through the +// solve pipeline, the empty operand included), plus: // - a zero on the *input* diagonal with nonsingular leading minors factors // fine (the pivots are the updated Schur-complement entries), // - a zero-*pivot* lane poisons itself with Inf/NaN without contaminating @@ -122,10 +124,11 @@ template static int run_case(int nm, int n, char uplo, char lay } // pack the full symmetric A, factor with the routine under test, unpack - std::vector ap = pack_compact(A, n, V, rowmajor); - int info = compact::sytrfnp(layout, uplo, n, ap.data(), n, V, nm); + const int lda = std::max(1, n); + std::vector ap = pack_compact(A, lda, V, rowmajor); + int info = compact::sytrfnp(layout, uplo, n, ap.data(), lda, V, nm); MatrixBatch Aout(nm, n, n); - unpack_compact(Aout, ap.data(), n, V, rowmajor); + unpack_compact(Aout, ap.data(), lda, V, rowmajor); double e_fac = 0, e_rec = 0, e_untouched = 0, a_norm = 1; for (int idx = 0; idx < nm; ++idx) { @@ -176,7 +179,7 @@ static int run_solve(int nm, int n, int nrhs, char uplo, char layout) { const T eps = std::numeric_limits::epsilon(); const bool rowmajor = (layout == 'R' || layout == 'r'); - const int ldb = rowmajor ? nrhs : n; + const int lda = std::max(1, n), ldb = std::max(1, rowmajor ? nrhs : n); // known X, B = A X densely MatrixBatch A(nm, n, n), B(nm, n, nrhs); @@ -187,15 +190,15 @@ static int run_solve(int nm, int n, int nrhs, char uplo, char layout) matmul(A.view(idx), X, B.view(idx)); } - std::vector ap = pack_compact(A, n, V, rowmajor); + std::vector ap = pack_compact(A, lda, V, rowmajor); std::vector bp = pack_compact(B, ldb, V, rowmajor); std::vector ap2 = ap, bp2 = bp; // the fused call's copies - int info_f = compact::sytrfnp(layout, uplo, n, ap.data(), n, V, nm); + int info_f = compact::sytrfnp(layout, uplo, n, ap.data(), lda, V, nm); int info_s = - compact::sytrsnp(layout, uplo, n, nrhs, ap.data(), n, bp.data(), ldb, V, nm); - int info_v = - compact::sysvnp(layout, uplo, n, nrhs, ap2.data(), n, bp2.data(), ldb, V, nm); + compact::sytrsnp(layout, uplo, n, nrhs, ap.data(), lda, bp.data(), ldb, V, nm); + int info_v = compact::sysvnp(layout, uplo, n, nrhs, ap2.data(), lda, bp2.data(), + ldb, V, nm); MatrixBatch Xhat(nm, n, nrhs); unpack_compact(Xhat, bp.data(), ldb, V, rowmajor); @@ -436,6 +439,18 @@ int main() fails += run_solve(16, 24, 3, u, l); } + // The small-dimension sweep (small_dims): every order over uplo x layout, + // then n x nrhs through the solve pipeline, a padded group each. + for (int n : small_dims) + for (char u : uplos) + for (char l : lays) + fails += run_case(5, n, u, l); + for (int n : small_dims) + for (int nrhs : small_dims) { + fails += run_solve(5, n, nrhs, 'L', 'C'); + fails += run_solve(9, n, nrhs, 'U', 'R'); + } + // nrhs = 0 must factor anyway (LAPACK ?sysv), bit-identical to sytrfnp, // bp a never-referenced dummy. fails += run_nrhs0(6, 20, 'L', 'C'); diff --git a/tests/test_trsm_compact.cpp b/tests/test_trsm_compact.cpp index 7057cde..e4c84aa 100644 --- a/tests/test_trsm_compact.cpp +++ b/tests/test_trsm_compact.cpp @@ -12,7 +12,9 @@ // interleave widths (column-major, the tuned path): forward error vs the // library's solve, plus the solve's own residual ||op(A) X - alpha B|| // formed with the library's triangular multiply (so a bug shared by the -// reference and the kernel cannot pass unseen). +// reference and the kernel cannot pass unseen), over the suite's shapes +// and the small-dimension cross product (small_dims: m x n from 0 to 5, +// the empty operand and the 1/2/3-column route included). // // Assisted-by: Claude:claude-opus-4-8 Claude:claude-fable-5 @@ -61,14 +63,15 @@ static int run_case(char side, char uplo, char transa, char diag, int nm, int m, } // pack, solve with the routine under test, unpack - std::vector ap = pack_compact(A, s, V); - std::vector bp = pack_compact(B, m, V); + const int lda = std::max(1, s), ldb = std::max(1, m); + std::vector ap = pack_compact(A, lda, V); + std::vector bp = pack_compact(B, ldb, V); - int info = compact::trsm('C', side, uplo, transa, diag, m, n, alpha, ap.data(), s, - bp.data(), m, V, nm); + int info = compact::trsm('C', side, uplo, transa, diag, m, n, alpha, ap.data(), + lda, bp.data(), ldb, V, nm); MatrixBatch Bout(nm, m, n); - unpack_compact(Bout, bp.data(), m, V); + unpack_compact(Bout, bp.data(), ldb, V); // Two gates: (1) forward error vs the library's solve, and (2) the // solve's own defining residual ||op(A) X - alpha B||, formed with the @@ -177,5 +180,14 @@ int main() fails += run_case('L', 'U', 'N', 'N', 16, 16, 1); fails += run_case('L', 'L', 'N', 'N', 32, 12, 2); + // the small-dimension cross product (small_dims): m x n for one flag + // combination per side, the empty operand and the 1/2/3-column route + // included, a padded group + for (int mm : small_dims) + for (int nn : small_dims) { + fails += run_case('L', 'U', 'N', 'N', 5, mm, nn); + fails += run_case('R', 'L', 'T', 'U', 9, mm, nn); + } + return finish(fails); } From 413cf8bba7b0f2789c248920fb80ebcec42e150c Mon Sep 17 00:00:00 2001 From: Claude Date: Mon, 14 Sep 2026 23:52:00 +0000 Subject: [PATCH 2/6] tests: gate every portable check on a LAPACK test ratio and one THRESH The portable suites gated each check on its own hand-picked tolerance -- 200 eps m, 20 n eps, 1e5 eps m, 1e4 eps p, 40 eps n ||A||, ... -- so no two checks were comparable, a near miss was invisible, and the cond(A)-dependent forward errors carried a guessed headroom. LAPACK's TESTING suites use one convention instead: every check is a dimensionless ratio, the error over what backward stability allows it to be, (order) * (operand norm) * eps, and a single threshold decides. test_compact_util.hpp now carries that convention: test_ratio(err, dim, norm), the residual form of dqrt01 / dpot01 / dget02; forward_ratio(err, normx, rcond), dget04's forward error discounted by rcond; THRESH = 30, the value of LAPACK's dtest.in; and passes(), written so NaN fails. test_lapack_util.hpp adds rcond1 -- 1 / (||A|| ||A^-1||) by LAPACKE ?getrf + ?gecon on the dense, unpacked input, since a condition number cannot be taken from the packed format -- and solve_ratios, the dget02 + dget04 pair the solve pipelines gate on. Every check in the seven portable suites is now such a ratio, printed with its case: the factorization residuals and reconstructions relative to ||A||, the elementwise comparisons against LAPACKE relative to the operand (1 for an orthogonal Q), the solve residuals to ||A|| ||X||, and the forward errors -- geqrf's and ormqr's back substitution, potrf's and sytrfnp's solve pipelines, trsm's solve vs the library's (rcond of the operator actually applied, the unit diagonal honored), gels's Gram-formed minimum-norm solution (rcond of the Gram matrix, the cond(A)^2 the old 1e4 headroom stood for) -- discounted by rcond. The structural contracts (bit-for-bit fused == two-step, untouched triangles, identity padding lanes) stay exact. Observed ratios: at most 2.2 across the suites, in both precisions, against THRESH = 30. The MKL suites keep their relative gates for now (solve_errors stays for them). CLAUDE.md records the convention. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_014ZZ76AgYwQ3vtuyZ46tyMX --- .claude/CLAUDE.md | 12 ++++ tests/test_compact_util.hpp | 43 ++++++++++++ tests/test_gels_compact.cpp | 84 ++++++++++++------------ tests/test_geqrf_compact.cpp | 114 ++++++++++++++++---------------- tests/test_lapack_util.hpp | 67 +++++++++++++++++-- tests/test_orgqr_compact.cpp | 50 +++++++------- tests/test_ormqr_compact.cpp | 78 +++++++++++----------- tests/test_potrf_compact.cpp | 78 +++++++++++----------- tests/test_sytrfnp_compact.cpp | 115 ++++++++++++++++----------------- tests/test_trsm_compact.cpp | 51 ++++++++------- 10 files changed, 402 insertions(+), 290 deletions(-) diff --git a/.claude/CLAUDE.md b/.claude/CLAUDE.md index b62564f..4a05b49 100644 --- a/.claude/CLAUDE.md +++ b/.claude/CLAUDE.md @@ -251,6 +251,18 @@ workspace contract, or the benchmarks' threading. own kernel. Register blocking is written as a `JB`-templated block helper with `for (c < JB)` loops the compiler unrolls, not as hand-expanded `w0..w3` copies. +- **Test ratios, one threshold.** The portable suites gate every numerical + check the way LAPACK's TESTING suites do: a dimensionless test ratio -- the + error over what backward stability allows, `(order) * (operand norm) * eps` + (`test_ratio` in `tests/test_compact_util.hpp`; `forward_ratio` is + dget04's, a forward error discounted by `rcond` from `rcond1`, LAPACKE's + `?getrf` + `?gecon` on the dense, unpacked input) -- compared against the + single `THRESH = 30` of LAPACK's `dtest.in`. A correct kernel's ratios are + O(1) (observed ≤ 2.2 across the suites) in either precision; each case prints + them, so a near miss is visible. Add a check as a ratio, never as a hand-picked + `k * eps * n` tolerance; only structural contracts (bit-for-bit reproduction, + untouched storage, identity padding lanes) are exact. The MKL suites still + carry their own relative gates. - **Argument checking.** The MKL-style API (`cbk_*`) skips validation like MKL's own compact routines (`info` is a scalar, `0` on success). The portable C API (`cbk.h`) validates LAPACK-style, returning `-j` for a bad j-th diff --git a/tests/test_compact_util.hpp b/tests/test_compact_util.hpp index 452dc14..f98798f 100644 --- a/tests/test_compact_util.hpp +++ b/tests/test_compact_util.hpp @@ -22,6 +22,7 @@ #include #include #include +#include #include #include #include @@ -75,6 +76,48 @@ template T frand() // a divide-by-zero guard (a normal double near DBL_MIN), not a tolerance. constexpr double norm_floor = 1e-300; +// ----------------------- test ratios and the threshold ---------------- +// Every numerical check is a dimensionless test ratio in the convention of +// LAPACK's TESTING suites: the error divided by what backward stability +// allows it to be, (order) * (operand norm) * eps, and one threshold decides +// -- THRESH = 30, the value of LAPACK's dtest.in. A correct kernel's ratio is +// O(1) in either precision; a case prints its ratios next to its verdict, so +// a near miss is visible, and a ratio >= THRESH fails (so does NaN: the +// comparison is written to fail it). One knob tunes every suite; there are +// no per-check tolerances. Only structural contracts -- a bit-for-bit +// reproduction, storage the routine must not touch, the identity in the +// padding lanes -- stay exact. +// +// test_ratio(err, dim, norm) err / (max(1, dim) * norm * eps) +// the residual form (dqrt01, dpot01, dget02): err the worst entry of +// R - Q^T A, L L^T - A, B - A X, ...; dim the order it accumulates over; +// norm the operand(s) it is relative to (1 for an orthogonal Q) +// forward_ratio(err, normx, rcond) err * rcond / (normx * eps) +// the forward-error form (dget04): err = ||x - x_true||, discounted by +// rcond = 1 / (||A|| ||A^-1||), the amplification a solve is allowed +constexpr double THRESH = 30.0; + +template double test_ratio(double err, int dim, double norm = 1.0) +{ + return err / (std::max(1, dim) * std::max(norm, norm_floor) * + (double)std::numeric_limits::epsilon()); +} + +template double forward_ratio(double err, double normx, double rcond) +{ + return test_ratio(err * rcond, 1, normx); +} + +inline bool passes(double ratio) +{ + return ratio < THRESH; /* false for NaN */ +} + +inline const char *verdict(double ratio) +{ + return passes(ratio) ? "OK" : "FAIL"; +} + // max |a - b| over n elements. template double max_abs_diff(const T *a, const T *b, size_t n) { diff --git a/tests/test_gels_compact.cpp b/tests/test_gels_compact.cpp index d0ded44..5cbe453 100644 --- a/tests/test_gels_compact.cpp +++ b/tests/test_gels_compact.cpp @@ -8,11 +8,12 @@ // time, blocked by LAPACK above its crossover -- and, independently of it, // the properties that define the two solutions. // -// Checks per (T, V, layout, trans, shape), the design's 7.1 gates: +// Checks per (T, V, layout, trans, shape), the design's 7.1 gates, each a +// test ratio against THRESH (test_compact_util.hpp): // 1. X (all max(m,n) rows of B, the residual rows included) == ref_gels, -// relative to ||X_lapack||_1 at 100 max(m,n) eps +// relative to ||X_lapack||_1 // 2. the factorization left in A and the tau left in taup == ref_gels, -// relative to ||A||_1 at 100 max(m,n) eps +// relative to ||A||_1 // 3. least squares: the normal equations op(A)^T (B - op(A) X) = 0, and the // residual sums of squares in rows n..m-1 of B equal ||B - op(A) X||^2 // minimum norm: op(A) X = B, and X equals the minimum-norm solution formed @@ -49,7 +50,6 @@ int run_case(char layout, char trans, int nm, int m, int n, int nrhs) const bool overdet = (tall != tran); const int rows_op = tran ? n : m, cols_op = tran ? m : n; // op(A) is rows_op x cols_op - const T eps = std::numeric_limits::epsilon(); // random A (diagonal-boosted, so op(A) is well conditioned) and random B; // the reference works on copies @@ -79,20 +79,24 @@ int run_case(char layout, char trans, int nm, int m, int n, int nrhs) unpack_tau(tau_out, tp.data(), V); // check 1 & 2: X (all p rows), the factorization, and tau vs the - // reference, each relative to the reference operand's L1 norm - double e_x = 0, e_h = 0, e_t = 0; + // reference, relative to the reference operand's L1 norm (tau, O(1) + // reflector scalars, to 1) + double r_x = 0, r_h = 0, r_t = 0; for (int idx = 0; idx < nm; ++idx) { - const double na = std::max(norm1(A.view(idx)), norm_floor); - e_x = std::max(e_x, max_abs_diff(Bout[idx], Bref[idx], B.stride()) / - std::max(norm1(Bref.view(idx)), norm_floor)); - e_h = std::max(e_h, max_abs_diff(Aout[idx], Aref[idx], A.stride()) / na); - e_t = std::max(e_t, max_abs_diff(tau_out[idx], tau_ref[idx], (size_t)q) / na); + const double na = norm1(A.view(idx)); + r_x = std::max(r_x, test_ratio(max_abs_diff(Bout[idx], Bref[idx], B.stride()), + p, norm1(Bref.view(idx)))); + r_h = std::max( + r_h, test_ratio(max_abs_diff(Aout[idx], Aref[idx], A.stride()), p, na)); + r_t = std::max( + r_t, test_ratio(max_abs_diff(tau_out[idx], tau_ref[idx], (size_t)q), p)); } - // check 3: the defining properties, formed without the reference: e_prop is - // the normal equations (least squares) or the residual (minimum norm), - // e_prop2 the residual-sum-of-squares rows or the Gram-formed solution - double e_prop = 0, e_prop2 = 0; + // check 3: the defining properties, formed without the reference: r_prop + // is the normal equations (least squares) or the residual (minimum norm), + // relative to ||op(A)|| ||B|| (dqrt17's scale), r_prop2 the + // residual-sum-of-squares rows or the Gram-formed solution + double r_prop = 0, r_prop2 = 0; std::vector Rs((size_t)rows_op * nrhs), Ss((size_t)cols_op * nrhs); const auto R = mat_view(Rs.data(), rows_op, nrhs); const auto S = mat_view(Ss.data(), cols_op, nrhs); @@ -100,8 +104,7 @@ int run_case(char layout, char trans, int nm, int m, int n, int nrhs) const auto Aop = tran ? A.view(idx).transposed() : A.view(idx); const auto Bin = leading(B.view(idx), rows_op, nrhs); // the right-hand sides const auto X = leading(Bout.view(idx), cols_op, nrhs); // the solution rows - const double scale = - std::max(norm1(Aop), norm_floor) * std::max(norm1(Bin), norm_floor); + const double scale = norm1(Aop) * norm1(Bin); matmul(Aop, X, R); // R = op(A) X if (overdet) { // r = B - op(A) X; normal equations op(A)^T r = 0 @@ -110,7 +113,8 @@ int run_case(char layout, char trans, int nm, int m, int n, int nrhs) R(i, j) = Bin(i, j) - R(i, j); matmul(Aop.transposed(), R, S); for (size_t e = 0; e < Ss.size(); ++e) - e_prop = std::max(e_prop, std::abs((double)Ss[e]) / scale); + r_prop = + std::max(r_prop, test_ratio(std::abs((double)Ss[e]), p, scale)); // rows cols_op..rows_op-1 of B hold the residual: squared column // norms == ||r_j||^2 (relative to ||b_j||^2: a square system has no // residual rows and a zero residual). Not at min(m,n) = 0: that is @@ -124,55 +128,55 @@ int run_case(char layout, char trans, int nm, int m, int n, int nrhs) } for (int i = cols_op; i < rows_op; ++i) rss_b += (double)Bout(idx, i, j) * Bout(idx, i, j); - e_prop2 = - std::max(e_prop2, std::abs(rss - rss_b) / std::max(bb, norm_floor)); + r_prop2 = std::max( + r_prop2, + test_ratio(std::abs(rss - rss_b) / std::max(bb, norm_floor), p)); } } else { // op(A) X = B for (int j = 0; j < nrhs; ++j) for (int i = 0; i < rows_op; ++i) - e_prop = - std::max(e_prop, std::abs((double)R(i, j) - Bin(i, j)) / scale); + r_prop = std::max( + r_prop, + test_ratio(std::abs((double)R(i, j) - Bin(i, j)), p, scale)); // minimum norm the other way: X = op(A)^T Z, G Z = B with the Gram // matrix G = op(A) op(A)^T (rows_op x rows_op, SPD), solved by the - // LAPACKE QR references + // LAPACKE QR references. The two solutions differ by what the + // Gram solve amplifies, cond(G) = cond(op(A))^2: dget04's forward + // error, discounted by rcond(G). std::vector Gs((size_t)rows_op * rows_op), Zs((size_t)rows_op * nrhs), tg(rows_op), Xmns((size_t)cols_op * nrhs); const auto G = mat_view(Gs.data(), rows_op, rows_op); const auto Z = mat_view(Zs.data(), rows_op, nrhs); const auto Xmn = mat_view(Xmns.data(), cols_op, nrhs); matmul(Aop, Aop.transposed(), G); + const double rcond_g = rcond1(G); copy_matrix(Bin, Z); ref_geqr2(G, tg.data()); ref_ormqr('T', rows_op, G, tg.data(), Z); ref_trsm_upper(G, Z); matmul(Aop.transposed(), Z, Xmn); - for (int j = 0; j < nrhs; ++j) - for (int i = 0; i < cols_op; ++i) - e_prop2 = std::max(e_prop2, std::abs((double)Xmn(i, j) - X(i, j)) / - std::max(norm1(X), norm_floor)); + r_prop2 = + std::max(r_prop2, forward_ratio(max_abs_diff(Xmns.data(), X.data, + (size_t)cols_op * nrhs), + norm1(X), rcond_g)); } } // One gate for X, the factorization and the defining properties: two // backward-stable implementations (LAPACK's blocked above its crossover) - // of a well-conditioned op(A)'s factorization and solve, relative to the - // operand norms (design 7.1). The Gram-formed minimum-norm solution - // carries cond(G) = cond(op(A))^2, so its gate gets that headroom. - const double tol = 100.0 * eps * p; - const double tol_prop2 = overdet ? tol : 1e4 * eps * p; - const bool ok_x = e_x <= tol, ok_h = e_h <= tol, ok_t = e_t <= tol; - const bool ok_p = e_prop <= tol, ok_p2 = e_prop2 <= tol_prop2; + // of a well-conditioned op(A)'s factorization and solve (design 7.1). + const bool ok_x = passes(r_x), ok_h = passes(r_h), ok_t = passes(r_t); + const bool ok_p = passes(r_prop), ok_p2 = passes(r_prop2); std::printf( - "T=%-6s V=%-2d %s trans=%c nm=%-2d m=%-3d n=%-3d nrhs=%d %-6s | X:%.1e %s " - "H:%.1e %s tau:%.1e %s | %s:%.1e %s %s:%.1e %s | info=%d\n", + "T=%-6s V=%-2d %s trans=%c nm=%-2d m=%-3d n=%-3d nrhs=%d %-6s | X:%-5.2g %s " + "H:%-5.2g %s tau:%-5.2g %s | %s:%-5.2g %s %s:%-5.2g %s | info=%d\n", compact::name, V, row ? "row" : "col", trans, nm, m, n, nrhs, - overdet ? "lstsq" : "minnrm", e_x, ok_x ? "OK" : "FAIL", e_h, - ok_h ? "OK" : "FAIL", e_t, ok_t ? "OK" : "FAIL", overdet ? "normal" : "resid", - e_prop, ok_p ? "OK" : "FAIL", overdet ? "rss" : "minnorm", e_prop2, - ok_p2 ? "OK" : "FAIL", info); + overdet ? "lstsq" : "minnrm", r_x, verdict(r_x), r_h, verdict(r_h), r_t, + verdict(r_t), overdet ? "normal" : "resid", r_prop, verdict(r_prop), + overdet ? "rss" : "minnorm", r_prop2, verdict(r_prop2), info); return (info != 0) + !ok_x + !ok_h + !ok_t + !ok_p + !ok_p2; } diff --git a/tests/test_geqrf_compact.cpp b/tests/test_geqrf_compact.cpp index 429bac8..50efced 100644 --- a/tests/test_geqrf_compact.cpp +++ b/tests/test_geqrf_compact.cpp @@ -6,17 +6,16 @@ // vectorized kernel executes V lanes at a time, so a correct kernel matches it // to working precision. // -// Checks per (T, V): -// 1. compact (H, tau) == LAPACKE_?geqr2 (H, tau) (elementwise, ~eps*scale) -// 2. reconstruction Q * triu(H) == A by ?ormqr (valid factorization) +// Checks per (T, V), each a test ratio against THRESH (test_compact_util.hpp): +// 1. compact (H, tau) == LAPACKE_?geqr2 (H, tau) (elementwise, ~eps) +// 2. reconstruction Q * triu(H) == A by ?ormqr (dqrt01's R - Q^T A) // 3. reflectors are usable: ormqr_compact('T') then ?trsm recovers a known -// X from B = A X (in-situ with ormqr) +// X from B = A X (dget04's forward error) // and, per (T, structure, cond) on the dense-LAPACK contract (design 7.1): -// 4. Q = ?orgqr(H, tau), R = triu(H): the factorization residual -// ||R - Q^T A||_1 / ||A||_1 <= 20 n eps and the orthogonality -// ||Q^T Q - I||_1 <= 100 n eps, on dense, column-scaled, rank-deficient -// and near-collinear inputs; (H, tau) elementwise vs the blocked -// LAPACKE_?geqrf on the dense ones (relative, 100 n eps) +// 4. Q = ?orgqr(H, tau), R = triu(H): dqrt01's factorization residual +// ||R - Q^T A|| / (m ||A|| eps) and orthogonality ||Q^T Q - I|| / (m eps), +// on dense, column-scaled, rank-deficient and near-collinear inputs; +// (H, tau) elementwise vs the blocked LAPACKE_?geqrf on the dense ones // over the suite's shapes and the small-dimension cross product (small_dims: // m x n from 0 to 5, so the empty operand and the single row or column run // through the same checks), plus LAPACK-style argument validation of the C @@ -41,7 +40,6 @@ using namespace cbk::test; template static int run_case(int nm, int m, int n) { const int k = std::min(m, n); - const T eps = std::numeric_limits::epsilon(); // random A batch, diagonal-boosted so the columns stay well conditioned MatrixBatch A(nm, m, n), Aref(nm, m, n), tau_ref(nm, k, 1); @@ -62,28 +60,33 @@ template static int run_case(int nm, int m, int n) unpack_compact(Aout, ap.data(), lda, V); unpack_tau(tau_out, tp.data(), V); - // check 1: (H, tau) match the scalar reference elementwise - double e_h = 0, e_t = 0; + // check 1: (H, tau) match the scalar reference elementwise -- the same + // operation sequence, so at ~eps; H relative to ||A||, tau (O(1) + // reflector scalars) to 1 + // check 2: reconstruction Q * triu(H) == A (dqrt01's R - Q^T A, by ?ormqr) + double r_h = 0, r_t = 0, r_rec = 0; for (int idx = 0; idx < nm; ++idx) { - e_h = std::max(e_h, max_abs_diff(Aout[idx], Aref[idx], (size_t)m * n)); - e_t = std::max(e_t, max_abs_diff(tau_out[idx], tau_ref[idx], (size_t)k)); - } + const double na = norm1(A.view(idx)); + r_h = std::max( + r_h, test_ratio(max_abs_diff(Aout[idx], Aref[idx], (size_t)m * n), m, na)); + r_t = std::max( + r_t, test_ratio(max_abs_diff(tau_out[idx], tau_ref[idx], (size_t)k), m)); - // check 2: reconstruction Q * triu(H) == A (valid factorization) - double e_rec = 0; - for (int idx = 0; idx < nm; ++idx) { std::vector Recs((size_t)m * n, T(0)); // start from R const auto Rec = mat_view(Recs.data(), m, n); for (int j = 0; j < n; ++j) for (int i = 0; i <= std::min(j, k - 1); ++i) Rec(i, j) = Aout(idx, i, j); ref_ormqr('N', k, Aout.view(idx), tau_out[idx], Rec); - e_rec = std::max(e_rec, max_abs_diff(Recs.data(), A[idx], (size_t)m * n)); + r_rec = std::max( + r_rec, + test_ratio(max_abs_diff(Recs.data(), A[idx], (size_t)m * n), m, na)); } // check 3: solve A X = B with the produced reflectors (square only), using - // the project's ormqr kernel + a triangular solve. X(:,j) = j+1. - double e_solve = -1; + // the project's ormqr kernel + a triangular solve. X(:,j) = j+1; dget04's + // forward error, discounted by rcond(A). + double r_solve = -1; if (m == n) { const int nrhs = 3; const std::vector Xs = known_solution(n, nrhs); @@ -96,25 +99,23 @@ template static int run_case(int nm, int m, int n) nm); MatrixBatch Bo(nm, n, nrhs); unpack_compact(Bo, bp.data(), lda, V); - e_solve = 0; + r_solve = 0; for (int idx = 0; idx < nm; ++idx) { ref_trsm_upper(Aout.view(idx), Bo.view(idx)); - e_solve = - std::max(e_solve, max_abs_diff(Bo[idx], Xs.data(), (size_t)n * nrhs)); + r_solve = std::max(r_solve, forward_ratio(max_abs_diff(Bo[idx], Xs.data(), + (size_t)n * nrhs), + norm1(X), rcond1(A.view(idx)))); } } - const double scale = std::max(1, m); - const double tol_fac = 200.0 * eps * scale; // same op sequence as ?geqr2, ~eps - const double tol_rec = 200.0 * eps * scale; // ?ormqr's rounding differs, same bound - const double tol_sol = 1e5 * eps * m; // cond(A)-dependent - bool ok_h = e_h <= tol_fac, ok_t = e_t <= tol_fac, ok_r = e_rec <= tol_rec; - bool ok_s = (e_solve < 0) || (e_solve <= tol_sol); - - std::printf("T=%-6s V=%-2d nm=%-2d m=%-3d n=%-3d | H:%.1e %s tau:%.1e %s rec:%.1e %s", - compact::name, V, nm, m, n, e_h, ok_h ? "OK" : "FAIL", e_t, - ok_t ? "OK" : "FAIL", e_rec, ok_r ? "OK" : "FAIL"); - if (e_solve >= 0) std::printf(" solve:%.1e %s", e_solve, ok_s ? "OK" : "FAIL"); + const bool ok_h = passes(r_h), ok_t = passes(r_t), ok_r = passes(r_rec); + const bool ok_s = (r_solve < 0) || passes(r_solve); + + std::printf( + "T=%-6s V=%-2d nm=%-2d m=%-3d n=%-3d | H:%-5.2g %s tau:%-5.2g %s rec:%-5.2g %s", + compact::name, V, nm, m, n, r_h, verdict(r_h), r_t, verdict(r_t), r_rec, + verdict(r_rec)); + if (r_solve >= 0) std::printf(" solve:%-5.2g %s", r_solve, verdict(r_solve)); std::printf(" | info=%d\n", info); return (info != 0) + !ok_h + !ok_t + !ok_r + !ok_s; } @@ -154,7 +155,6 @@ template static void gen_matrix(MatrixView A, double cond, Structur template static int run_invariants(int nm, int m, int n, double cond, Structure structure = DENSE) { - const double eps = std::numeric_limits::epsilon(); const int k = std::min(m, n); MatrixBatch A(nm, m, n); @@ -168,7 +168,7 @@ static int run_invariants(int nm, int m, int n, double cond, Structure structure unpack_compact(H, ap.data(), m, V); unpack_tau(tau, tp.data(), V); - double worst_res = 0, worst_orth = 0, worst_el = 0; + double r_res = 0, r_orth = 0, r_el = 0; std::vector Qs((size_t)m * k), QtAs((size_t)k * n), QtQs((size_t)k * k), Hrefs(A.stride()), tauref(k); const auto Q = mat_view(Qs.data(), m, k), QtA = mat_view(QtAs.data(), k, n), @@ -181,20 +181,21 @@ static int run_invariants(int nm, int m, int n, double cond, Structure structure copy_matrix(leading(Hm, m, k), Q); ref_orgqr(k, Q, tau[idx]); - // residual R - Q^T A, R = triu(H) (k x n) + // residual R - Q^T A, R = triu(H) (k x n): dqrt01's + // ||R - Q^T A|| / (m ||A|| eps) matmul(Q.transposed(), Am, QtA); double resid = 0; for (int j = 0; j < n; ++j) for (int i = 0; i < k; ++i) resid = std::max( resid, (double)std::abs(((i <= j) ? Hm(i, j) : T(0)) - QtA(i, j))); - worst_res = std::max(worst_res, resid / std::max(norm1(Am), norm_floor)); + r_res = std::max(r_res, test_ratio(resid, m, norm1(Am))); - // orthogonality Q^T Q - I + // orthogonality Q^T Q - I: dqrt01's ||I - Q^T Q|| / (m eps) matmul(Q.transposed(), Q, QtQ); for (int d = 0; d < k; ++d) QtQ(d, d) -= T(1); - worst_orth = std::max(worst_orth, norm1(QtQ)); + r_orth = std::max(r_orth, test_ratio(norm1(QtQ), m)); // elementwise vs the blocked LAPACKE_?geqrf, on the dense inputs only, // where the reflectors are essentially unique: it catches a @@ -207,21 +208,19 @@ static int run_invariants(int nm, int m, int n, double cond, Structure structure ref_geqrf(Href, tauref.data()); const double el = std::max(max_abs_diff(H[idx], Hrefs.data(), A.stride()), max_abs_diff(tau[idx], tauref.data(), (size_t)k)); - worst_el = std::max(worst_el, el / std::max(norm1(Am), norm_floor)); + r_el = std::max(r_el, test_ratio(el, n, norm1(Am))); } } - const double rtol_res = 20.0 * n * eps, rtol_orth = 100.0 * n * eps; - const double rtol_el = 100.0 * n * eps; - const bool ok = (info == 0) && (worst_res <= rtol_res) && (worst_orth <= rtol_orth) && - (worst_el <= rtol_el); + const bool ok = (info == 0) && passes(r_res) && passes(r_orth) && passes(r_el); const char *sname = structure == DENSE ? "dense" : structure == RANK_DEFICIENT ? "rankdef" : "collin"; - std::printf("T=%-6s V=%-2d nm=%-2d m=%-3d n=%-3d cond=%.0f %-8s| res %.2e (%.1e) " - "orth %.2e (%.1e) el %.1e | info=%d %s\n", - compact::name, V, nm, m, n, cond, sname, worst_res, rtol_res, - worst_orth, rtol_orth, worst_el, info, ok ? "OK" : "FAIL"); + std::printf( + "T=%-6s V=%-2d nm=%-2d m=%-3d n=%-3d cond=%.0f %-8s| res %-5.2g orth %-5.2g " + "el %-5.2g | info=%d %s\n", + compact::name, V, nm, m, n, cond, sname, r_res, r_orth, r_el, info, + ok ? "OK" : "FAIL"); return !ok; } @@ -260,7 +259,6 @@ template static int test_underflow() { const int nm = V, m = 6, n = 4, tiny_col = 0; /* no earlier reflector touches it */ const T tiny = std::is_same_v ? T(1e-170) : T(1e-25); - const T eps = std::numeric_limits::epsilon(); MatrixBatch A(nm, m, n); for (int idx = 0; idx < nm; ++idx) { @@ -277,7 +275,7 @@ template static int test_underflow() unpack_tau(tau, tp.data(), V); int fails = (info != 0); - double e_rec = 0; + double r_rec = 0; for (int idx = 0; idx < nm; ++idx) { /* column 0 is the first reflector's, so nothing reduces it first: its * tail underflows squared, tau = 0, the diagonal is kept, the body @@ -298,14 +296,14 @@ template static int test_underflow() for (int i = 0; i <= std::min(j, n - 1); ++i) Rec(i, j) = Aout(idx, i, j); ref_ormqr('N', n, Aout.view(idx), tau[idx], Rec); - e_rec = std::max(e_rec, max_abs_diff(Recs.data(), A[idx], (size_t)m * n)); + r_rec = std::max(r_rec, + test_ratio(max_abs_diff(Recs.data(), A[idx], (size_t)m * n), + m, norm1(A.view(idx)))); } - const double tol_rec = 200.0 * eps * m; - fails += !(e_rec <= tol_rec); + fails += !passes(r_rec); std::printf("T=%-6s V=%-2d underflow column: tau=0, diag kept, body 0, finite %s | " - "rec:%.1e %s\n", - compact::name, V, fails ? "FAIL" : "OK", e_rec, - e_rec <= tol_rec ? "OK" : "FAIL"); + "rec:%-5.2g %s\n", + compact::name, V, fails ? "FAIL" : "OK", r_rec, verdict(r_rec)); return fails; } diff --git a/tests/test_lapack_util.hpp b/tests/test_lapack_util.hpp index 2d9f233..ca7ec1a 100644 --- a/tests/test_lapack_util.hpp +++ b/tests/test_lapack_util.hpp @@ -53,8 +53,8 @@ namespace cbk::test { // ----------------------- LAPACKE / CBLAS, by scalar type ------------- -// lapack::geqr2 / geqrf / gelqf / orgqr / ormqr / geqp3 / potrf / gels -// forward to LAPACKE_?*, gemm / trsm / trmm to cblas_?*. One macro +// lapack::geqr2 / geqrf / gelqf / orgqr / ormqr / geqp3 / potrf / getrf / +// gecon / gels forward to LAPACKE_?*, gemm / trsm / trmm to cblas_?*. One macro // generates both specializations so the two cannot drift apart. template struct lapack; @@ -84,6 +84,12 @@ template <> struct lapack { { return LAPACKE_##p##geqp3(layout, m, n, a, lda, jpvt, tau); } \ static lapack_int potrf(int layout, char uplo, lapack_int n, T *a, lapack_int lda) \ { return LAPACKE_##p##potrf(layout, uplo, n, a, lda); } \ + static lapack_int getrf(int layout, lapack_int m, lapack_int n, T *a, lapack_int lda, \ + lapack_int *ipiv) \ + { return LAPACKE_##p##getrf(layout, m, n, a, lda, ipiv); } \ + static lapack_int gecon(int layout, char norm, lapack_int n, const T *a, \ + lapack_int lda, T anorm, T *rcond) \ + { return LAPACKE_##p##gecon(layout, norm, n, a, lda, anorm, rcond); } \ static lapack_int gels(int layout, char trans, lapack_int m, lapack_int n, \ lapack_int nrhs, T *a, lapack_int lda, T *b, lapack_int ldb) \ { return LAPACKE_##p##gels(layout, trans, m, n, nrhs, a, lda, b, ldb); } \ @@ -398,10 +404,59 @@ void tri_apply(char side, char uplo, char transa, char diag, Av A, Xv X, Rv R) // ----------------------- shared checks ------------------------------- -// The end-to-end solve gate every solving suite closes with: the worst -// relative forward error of Xhat against the known X, and the worst relative -// residual ||A Xhat - B|| formed by ?gemm -- so a solver bug cannot hide -// behind the factorization that produced Xhat. +// rcond = 1 / (||A||_1 ||A^-1||_1) of a square matrix, by ?getrf + ?gecon: +// the RCOND LAPACK's dget04 discounts a forward error by (forward_ratio). +// An empty matrix is perfectly conditioned; a singular one returns 0, which +// makes any forward error pass -- the ratio then says nothing, as in LAPACK. +template double rcond1(Av A) +{ + using T = elem_t; + assert(A.rows == A.cols); + const int n = A.rows; + if (n == 0) return 1.0; + std::vector LUs((std::size_t)n * n); + std::vector ipiv(n); + const auto LU = mat_view(LUs.data(), n, n); + copy_matrix(A, LU); + const double anorm = norm1(LU); + lapack::getrf(LAPACK_COL_MAJOR, n, n, LU.data, n, ipiv.data()); + T rcond = 0; + lapack::gecon(LAPACK_COL_MAJOR, '1', n, LU.data, n, (T)anorm, &rcond); + return (double)rcond; +} + +// The end-to-end solve gate every solving suite closes with, as test ratios: +// dget02's residual ||B - A Xhat|| / (n ||A|| ||Xhat|| eps), formed by ?gemm +// so a solver bug cannot hide behind the factorization that produced Xhat, +// and dget04's forward error ||Xhat - X|| rcond(A) / (||X|| eps) against the +// known X. The worst over the batch of each. +struct SolveRatios { + double res, fwd; +}; + +template +SolveRatios solve_ratios(const MatrixBatch &A, const MatrixBatch &B, + const MatrixBatch &Xhat, ConstMatrixView X) +{ + assert(B.rows() == Xhat.rows() && B.cols() == Xhat.cols() && X.rows == Xhat.rows() && + X.cols == Xhat.cols()); + const int n = A.rows(); + const size_t sB = Xhat.stride(); + std::vector AXs(sB); + const auto AX = mat_view(AXs.data(), Xhat.rows(), Xhat.cols()); + const double nX = norm1(X); + SolveRatios r{0, 0}; + for (int v = 0; v < Xhat.count(); ++v) { + matmul(A.view(v), Xhat.view(v), AX); + r.res = std::max(r.res, test_ratio(max_abs_diff(AXs.data(), B[v], sB), n, + norm1(A.view(v)) * norm1(Xhat.view(v)))); + r.fwd = std::max(r.fwd, forward_ratio(max_abs_diff(Xhat[v], X.data, sB), nX, + rcond1(A.view(v)))); + } + return r; +} + +// The same two, as relative errors (the MKL suites' gates). struct SolveErrors { double fwd, res; }; diff --git a/tests/test_orgqr_compact.cpp b/tests/test_orgqr_compact.cpp index 6d09867..998165f 100644 --- a/tests/test_orgqr_compact.cpp +++ b/tests/test_orgqr_compact.cpp @@ -8,13 +8,14 @@ * reflectors into the first k columns of an m x n buffer whose remaining * columns hold garbage (they must be overwritten), and generate Q. * - * Checks per (T, V, layout): - * 1. compact Q == LAPACKE_?orgqr Q (elementwise, 100 m eps: Q's entries - * are O(1) and, for full-rank A, Q is unique up to the shared sign - * convention, so two backward-stable routes agree to cond(A) * eps; + * Checks per (T, V, layout), each a test ratio against THRESH + * (test_compact_util.hpp): + * 1. compact Q == LAPACKE_?orgqr Q (elementwise, relative to 1: Q's + * entries are O(1) and, for full-rank A, Q is unique up to the shared + * sign convention, so two backward-stable routes agree to cond(A) * eps; * LAPACK blocks its accumulation above k = 128, so the rounding order * differs) - * 2. Q^T Q = I (orthonormal columns, formed independently) + * 2. Q^T Q = I (orthonormal columns, formed independently: dqrt02's) * 3. Q(:, 0:k-1) R = A with R = triu of the factorization (k > 0) * 4. padded lanes of a partial last group come out exactly identity * plus k < n (unit-seeded extra columns), k = 0 (Q = leading columns of I), @@ -43,10 +44,6 @@ template static int run_case(int nm, int m, int n, int k, bool rowmajor = false) { assert(m >= n && n >= k); - const T eps = std::numeric_limits::epsilon(); - const double tol_el = 100.0 * eps * m; /* vs LAPACKE_?orgqr, see the header */ - const double tol_orth = 100.0 * eps * m; /* independent invariants */ - const double tol_rec = 100.0 * eps * m; const int ld = std::max(1, rowmajor ? n : m); @@ -80,18 +77,21 @@ static int run_case(int nm, int m, int n, int k, bool rowmajor = false) } unpack_compact(Qout, ap.data(), ld, V, rowmajor); - /* check 1: elementwise vs LAPACKE_?orgqr */ - double e1 = 0; + /* check 1: elementwise vs LAPACKE_?orgqr (Q's entries are O(1): relative + * to 1, the norm of an orthogonal matrix) */ + double r1 = 0; for (int kk = 0; kk < nm; ++kk) - e1 = std::max(e1, max_abs_diff(Qout[kk], Qref[kk], Qout.stride())); + r1 = std::max(r1, + test_ratio(max_abs_diff(Qout[kk], Qref[kk], Qout.stride()), m)); - /* check 2: Q^T Q = I, formed densely */ - double e2 = 0; + /* check 2: Q^T Q = I, formed densely (dqrt02's ||I - Q^T Q|| / (m eps)) */ + double r2 = 0; for (int kk = 0; kk < nm; ++kk) - e2 = std::max(e2, orth_error(Qout.view(kk))); + r2 = std::max(r2, test_ratio(orth_error(Qout.view(kk)), m)); - /* check 3: Q(:, 0:k-1) R = A0, R = triu of the factorization */ - double e3 = 0; + /* check 3: Q(:, 0:k-1) R = A0, R = triu of the factorization, relative + * to ||A0|| */ + double r3 = 0; if (k > 0) { std::vector Rs((size_t)k * k), Recs((size_t)m * k); const auto R = mat_view(Rs.data(), k, k); @@ -101,8 +101,9 @@ static int run_case(int nm, int m, int n, int k, bool rowmajor = false) for (int i = 0; i < k; ++i) R(i, j) = (i <= j) ? Afac(kk, i, j) : T(0); matmul(leading(Qout.view(kk), m, k), R, Rec); - e3 = std::max(e3, max_abs_diff(Recs.data(), A0[kk], (size_t)m * k) / - std::max(norm1(A0.view(kk)), norm_floor)); + r3 = std::max(r3, + test_ratio(max_abs_diff(Recs.data(), A0[kk], (size_t)m * k), + m, norm1(A0.view(kk)))); } } @@ -120,13 +121,12 @@ static int run_case(int nm, int m, int n, int k, bool rowmajor = false) (double)std::abs(P(i, j)[v] - (i == j ? T(1) : T(0)))); } - bool ok1 = e1 <= tol_el, ok2 = e2 <= tol_orth, ok3 = e3 <= tol_rec, ok4 = e4 == 0; + const bool ok1 = passes(r1), ok2 = passes(r2), ok3 = passes(r3), ok4 = e4 == 0; fails += !ok1 + !ok2 + !ok3 + !ok4; - std::printf("T=%-6s V=%-2d nm=%-2d %s m=%-3d n=%-3d k=%-3d | ref: %.2e %s | QtQ=I: " - "%.2e %s | QR=A: %.2e %s | pad: %.1e %s\n", - compact::name, V, nm, rowmajor ? "row" : "col", m, n, k, e1, - ok1 ? "OK" : "FAIL", e2, ok2 ? "OK" : "FAIL", e3, ok3 ? "OK" : "FAIL", e4, - ok4 ? "OK" : "FAIL"); + std::printf("T=%-6s V=%-2d nm=%-2d %s m=%-3d n=%-3d k=%-3d | ref: %-5.2g %s | QtQ=I: " + "%-5.2g %s | QR=A: %-5.2g %s | pad: %.1e %s\n", + compact::name, V, nm, rowmajor ? "row" : "col", m, n, k, r1, + verdict(r1), r2, verdict(r2), r3, verdict(r3), e4, ok4 ? "OK" : "FAIL"); return fails; } diff --git a/tests/test_ormqr_compact.cpp b/tests/test_ormqr_compact.cpp index 4fb5ea2..ea5a415 100644 --- a/tests/test_ormqr_compact.cpp +++ b/tests/test_ormqr_compact.cpp @@ -8,11 +8,11 @@ * Test design: X(:,j) = j+1 (ones, twos, threes, ...), B = A*X, so the * b columns are scaled row sums of A; the solve must recover X. * - * Checks per (T, V): - * 1. compact Q^T B == LAPACKE_?ormqr Q^T B, relative to ||Q^T B||_1 at - * 20 m eps (the design's 7.1 gate: two backward-stable applications of - * the same reflectors, LAPACK's blocked above its crossover) - * 2. back substitution recovers X + * Checks per (T, V), each a test ratio against THRESH (test_compact_util.hpp): + * 1. compact Q^T B == LAPACKE_?ormqr Q^T B, relative to ||Q^T B||_1 (the + * design's 7.1 gate: two backward-stable applications of the same + * reflectors, LAPACK's blocked above its crossover) + * 2. back substitution recovers X (dget04's forward error, rcond(A)) * 3. applying 'N' after 'T' recovers the original B (Q Q^T = I) * over the suite's shapes and the small-dimension cross product (small_dims: * m x nrhs from 0 to 5, the empty operand and the 1/2/3-column right-hand @@ -39,10 +39,6 @@ using namespace cbk::test; template static int run_case(int nm, int m, int nrhs) { const int k = m; - const T eps = std::numeric_limits::epsilon(); - const double tol_el = 20.0 * eps * m; /* vs LAPACKE_?ormqr, relative */ - const double tol_exact = 100.0 * eps; /* the 'N' after 'T' round trip */ - const double tol_solve = 1e5 * eps * m; /* cond(A)-dependent */ const std::vector Xs = known_solution(m, nrhs); const auto X = mat_view(Xs.data(), m, nrhs); @@ -67,31 +63,35 @@ template static int run_case(int nm, int m, int nrhs) /* check 1: compact Q^T B vs LAPACKE_?ormqr, relative to ||Q^T B||_1 */ compact::ormqr('T', m, nrhs, k, ap.data(), ld, tp.data(), bp.data(), ld, V, nm); unpack_compact(Bout, bp.data(), ld, V); - double e1 = 0; + double r1 = 0; for (int kk = 0; kk < nm; ++kk) - e1 = std::max(e1, max_abs_diff(Bout[kk], Bref[kk], (size_t)m * nrhs) / - std::max(norm1(Bref.view(kk)), norm_floor)); + r1 = + std::max(r1, test_ratio(max_abs_diff(Bout[kk], Bref[kk], (size_t)m * nrhs), + m, norm1(Bref.view(kk)))); - /* check 2: solve recovers X */ - double e2 = 0; + /* check 2: solve recovers X (dget04, discounted by rcond(A)) */ + double r2 = 0; for (int kk = 0; kk < nm; ++kk) { ref_trsm_upper(Afac.view(kk), Bout.view(kk)); - e2 = std::max(e2, max_abs_diff(Bout[kk], Xs.data(), (size_t)m * nrhs)); + r2 = std::max( + r2, forward_ratio(max_abs_diff(Bout[kk], Xs.data(), (size_t)m * nrhs), + norm1(X), rcond1(A.view(kk)))); } - /* check 3: 'N' undoes 'T' */ + /* check 3: 'N' undoes 'T', relative to ||B||_1 */ compact::ormqr('N', m, nrhs, k, ap.data(), ld, tp.data(), bp.data(), ld, V, nm); unpack_compact(Bout, bp.data(), ld, V); - double e3 = 0; + double r3 = 0; for (int kk = 0; kk < nm; ++kk) - e3 = std::max(e3, max_abs_diff(Bout[kk], B[kk], (size_t)m * nrhs)); - - /* scale-aware: B entries are O(m), QQ^t roundtrip accumulates a bit */ - bool ok1 = e1 <= tol_el, ok2 = e2 <= tol_solve, ok3 = e3 <= tol_exact * m * 10; - std::printf("T=%-6s V=%-2d nm=%-2d m=%-3d nrhs=%d | QtB: %.2e %s | solve X: %.2e %s " - "| QQt=I: %.2e %s\n", - compact::name, V, nm, m, nrhs, e1, ok1 ? "OK" : "FAIL", e2, - ok2 ? "OK" : "FAIL", e3, ok3 ? "OK" : "FAIL"); + r3 = std::max(r3, test_ratio(max_abs_diff(Bout[kk], B[kk], (size_t)m * nrhs), + m, norm1(B.view(kk)))); + + const bool ok1 = passes(r1), ok2 = passes(r2), ok3 = passes(r3); + std::printf( + "T=%-6s V=%-2d nm=%-2d m=%-3d nrhs=%d | QtB: %-5.2g %s | solve X: %-5.2g %s " + "| QQt=I: %-5.2g %s\n", + compact::name, V, nm, m, nrhs, r1, verdict(r1), r2, verdict(r2), r3, + verdict(r3)); return !ok1 + !ok2 + !ok3; } @@ -105,21 +105,18 @@ template static int run_case(int nm, int m, int nrhs) template static int run_case_pivoted(int nm, int m, int nrhs) { const int k = m; /* square, full rank */ - const T eps = std::numeric_limits::epsilon(); - const double tol_solve = 1e5 * eps * m; /* cond(A)-dependent */ const std::vector Xs = known_solution(m, nrhs); const auto X = mat_view(Xs.data(), m, nrhs); - MatrixBatch Afac(nm, m, m), B(nm, m, nrhs), Bout(nm, m, nrhs), tau(nm, m, 1); + MatrixBatch A(nm, m, m), Afac(nm, m, m), B(nm, m, nrhs), Bout(nm, m, nrhs), + tau(nm, m, 1); MatrixBatch jpvt(nm, m, 1); for (int kk = 0; kk < nm; ++kk) { - std::vector As((size_t)m * m); - const auto A = mat_view(As.data(), m, m); - gen_boosted(A); /* diagonal boost tames cond */ - matmul(A, X, B.view(kk)); /* B = A X */ + gen_boosted(A.view(kk)); /* diagonal boost tames cond */ + matmul(A.view(kk), X, B.view(kk)); /* B = A X */ - std::copy(As.begin(), As.end(), Afac[kk]); + std::copy(A[kk], A[kk] + A.stride(), Afac[kk]); ref_geqp3(Afac.view(kk), jpvt[kk], tau[kk]); } @@ -131,8 +128,9 @@ template static int run_case_pivoted(int nm, int m, int nrhs) compact::ormqr('T', m, nrhs, k, ap.data(), m, tp.data(), bp.data(), m, V, nm); unpack_compact(Bout, bp.data(), m, V); - /* R y = c, then back-permute x(jpvt(j)) = y(j); compare against X */ - double e = 0; + /* R y = c, then back-permute x(jpvt(j)) = y(j); compare against X (dget04, + * discounted by rcond(A)) */ + double r = 0; std::vector xs((size_t)m * nrhs); const auto x = mat_view(xs.data(), m, nrhs); for (int kk = 0; kk < nm; ++kk) { @@ -141,12 +139,14 @@ template static int run_case_pivoted(int nm, int m, int nrhs) for (int j = 0; j < nrhs; ++j) for (int i = 0; i < m; ++i) x(jpvt[kk][i], j) = y(i, j); - e = std::max(e, max_abs_diff(xs.data(), Xs.data(), (size_t)m * nrhs)); + r = std::max( + r, forward_ratio(max_abs_diff(xs.data(), Xs.data(), (size_t)m * nrhs), + norm1(X), rcond1(A.view(kk)))); } - bool ok = e <= tol_solve; - std::printf("T=%-6s V=%-2d nm=%-2d m=%-3d nrhs=%d | pivoted solve X: %.2e %s\n", - compact::name, V, nm, m, nrhs, e, ok ? "OK" : "FAIL"); + const bool ok = passes(r); + std::printf("T=%-6s V=%-2d nm=%-2d m=%-3d nrhs=%d | pivoted solve X: %-5.2g %s\n", + compact::name, V, nm, m, nrhs, r, verdict(r)); return !ok; } diff --git a/tests/test_potrf_compact.cpp b/tests/test_potrf_compact.cpp index 03dc448..6411d25 100644 --- a/tests/test_potrf_compact.cpp +++ b/tests/test_potrf_compact.cpp @@ -7,11 +7,13 @@ // (recursive / blocked, so its rounding differs from the kernel's // right-looking sweep, but the SPD factor is unique). // -// Checks per (T, V, uplo, layout): +// Checks per (T, V, uplo, layout), each a test ratio against THRESH +// (test_compact_util.hpp): // 1. named-triangle factor == LAPACKE_?potrf factor, relative to the -// reference factor's L1 norm at 20 n eps (design 7.1: the unique factor -// agrees to a few n eps between backward-stable implementations) +// reference factor's L1 norm (design 7.1: the unique factor agrees to a +// few n eps between backward-stable implementations) // 2. reconstruction L L^T == A (lower) / U^T U == A (upper), by ?trmm +// (dpot01's residual) // 3. the strictly-opposite triangle of the compact buffer is bit-for-bit // unchanged from the input (the routine must not reference or write it) // 4. end-to-end solve: factor + ?potrs_compact recovers a known X, and @@ -40,7 +42,6 @@ using namespace cbk::test; template static int run_case(int nm, int n, char uplo, char layout, double cond = 0.0) { - const T eps = std::numeric_limits::epsilon(); const bool rowmajor = (layout == 'R' || layout == 'r'); const bool upper = (uplo == 'U' || uplo == 'u'); @@ -60,7 +61,7 @@ static int run_case(int nm, int n, char uplo, char layout, double cond = 0.0) MatrixBatch Aout(nm, n, n); unpack_compact(Aout, ap.data(), lda, V, rowmajor); - double e_fac = 0, e_rec = 0, e_untouched = 0; + double r_fac = 0, r_rec = 0, e_untouched = 0; std::vector Fs((size_t)n * n), Recs((size_t)n * n); const auto F = mat_view(Fs.data(), n, n), Rec = mat_view(Recs.data(), n, n); for (int idx = 0; idx < nm; ++idx) { @@ -78,11 +79,12 @@ static int run_case(int nm, int n, char uplo, char layout, double cond = 0.0) e_untouched = std::max(e_untouched, (double)std::abs(Fac(i, j) - Ain(i, j))); } - e_fac = std::max(e_fac, el / std::max(norm1(Ref), norm_floor)); + r_fac = std::max(r_fac, test_ratio(el, n, norm1(Ref))); // check 2: reconstruction of A from the named triangle's factor, by // ?trmm: F is the factor with the other triangle zeroed, so - // L L^T = F L^T (side 'R') and U^T U = U^T F (side 'L') + // L L^T = F L^T (side 'R') and U^T U = U^T F (side 'L') -- dpot01's + // ||L L^T - A|| / (n ||A|| eps) for (int j = 0; j < n; ++j) for (int i = 0; i < n; ++i) F(i, j) = (upper ? (i <= j) : (i >= j)) ? Fac(i, j) : T(0); @@ -90,23 +92,19 @@ static int run_case(int nm, int n, char uplo, char layout, double cond = 0.0) tri_apply('R', 'L', 'T', 'N', Fac, F, Rec); else tri_apply('L', 'U', 'T', 'N', Fac, F, Rec); - e_rec = std::max(e_rec, max_abs_diff(Recs.data(), A[idx], (size_t)n * n) / - std::max(norm1(Ain), norm_floor)); + r_rec = std::max(r_rec, + test_ratio(max_abs_diff(Recs.data(), A[idx], (size_t)n * n), + n, norm1(Ain))); } - // design 7.1: the unique factor vs LAPACK's, and ||L L^T - A|| / ||A||, - // both backward-stable quantities at n eps - const double tol = 20.0 * eps * std::max(1, n); - bool ok_f = e_fac <= tol; - bool ok_r = e_rec <= tol; - bool ok_u = e_untouched == 0.0; // must be bit-for-bit unchanged - bool ok_i = (info == 0); - - std::printf("T=%-6s V=%-2d uplo=%c lay=%c nm=%-2d n=%-3d cond=%.0f | fac:%.1e %s " - "rec:%.1e %s untouched:%s | info=%d %s\n", - compact::name, V, uplo, layout, nm, n, cond, e_fac, - ok_f ? "OK" : "FAIL", e_rec, ok_r ? "OK" : "FAIL", ok_u ? "OK" : "FAIL", - info, ok_i ? "OK" : "FAIL"); + const bool ok_f = passes(r_fac), ok_r = passes(r_rec); + const bool ok_u = e_untouched == 0.0; // must be bit-for-bit unchanged + const bool ok_i = (info == 0); + + std::printf("T=%-6s V=%-2d uplo=%c lay=%c nm=%-2d n=%-3d cond=%.0f | fac:%-5.2g %s " + "rec:%-5.2g %s untouched:%s | info=%d %s\n", + compact::name, V, uplo, layout, nm, n, cond, r_fac, verdict(r_fac), + r_rec, verdict(r_rec), ok_u ? "OK" : "FAIL", info, ok_i ? "OK" : "FAIL"); return (!ok_f) + (!ok_r) + (!ok_u) + (!ok_i); } @@ -118,7 +116,6 @@ static int run_case(int nm, int n, char uplo, char layout, double cond = 0.0) template static int run_solve(int nm, int n, int nrhs, char uplo, char layout) { - const T eps = std::numeric_limits::epsilon(); const bool rowmajor = (layout == 'R' || layout == 'r'); const int lda = std::max(1, n), ldb = std::max(1, rowmajor ? nrhs : n); @@ -144,19 +141,19 @@ static int run_solve(int nm, int n, int nrhs, char uplo, char layout) MatrixBatch Xhat(nm, n, nrhs); unpack_compact(Xhat, bp.data(), ldb, V, rowmajor); - const auto [e_fwd, e_res] = solve_errors(A, B, Xhat, X); + // dget02's residual and dget04's forward error + const auto [r_res, r_fwd] = solve_ratios(A, B, Xhat, X); // fused vs two-step, on the raw compact buffers (padded lanes included) const bool fused_same = (ap2 == ap) && (bp2 == bp); - // The tame SPD batches keep cond(A) modest, so residual and forward error - // share one backward-stability-sized gate. - const double rtol = 100.0 * std::max(1, n) * eps; - bool ok = (e_fwd <= rtol) && (e_res <= rtol) && fused_same && (info_f == 0) && - (info_s == 0) && (info_v == 0); - std::printf("T=%-6s V=%-2d uplo=%c lay=%c nm=%-2d n=%-3d nrhs=%d solve | fwd:%.1e " - "res:%.1e (%.1e) posv==trf+trs:%s info=%d/%d/%d %s\n", - compact::name, V, uplo, layout, nm, n, nrhs, e_fwd, e_res, rtol, - fused_same ? "yes" : "NO", info_f, info_s, info_v, ok ? "OK" : "FAIL"); + const bool ok = passes(r_res) && passes(r_fwd) && fused_same && (info_f == 0) && + (info_s == 0) && (info_v == 0); + std::printf( + "T=%-6s V=%-2d uplo=%c lay=%c nm=%-2d n=%-3d nrhs=%d solve | res:%-5.2g %s " + "fwd:%-5.2g %s posv==trf+trs:%s info=%d/%d/%d %s\n", + compact::name, V, uplo, layout, nm, n, nrhs, r_res, verdict(r_res), r_fwd, + verdict(r_fwd), fused_same ? "yes" : "NO", info_f, info_s, info_v, + ok ? "OK" : "FAIL"); return !ok; } @@ -226,7 +223,7 @@ template static int run_nonspd(int n, char uplo, char layout) MatrixBatch Aout(nm, n, n); unpack_compact(Aout, ap.data(), n, V, rowmajor); - double e_spd = 0; // worst relative error over the SPD sibling lanes + double r_spd = 0; // worst ratio over the SPD sibling lanes bool spd_finite = true; bool bad_poisoned = false; // the non-SPD lane must carry NaN/Inf for (int idx = 0; idx < nm; ++idx) { @@ -244,17 +241,16 @@ template static int run_nonspd(int n, char uplo, char layout) el = std::max(el, (double)std::abs(x - Aref(idx, i, j))); } if (idx != badlane) // vs LAPACKE_?potrf, relative to its factor's norm - e_spd = std::max(e_spd, el / std::max(norm1(Aref.view(idx)), norm_floor)); + r_spd = std::max(r_spd, test_ratio(el, n, norm1(Aref.view(idx)))); } - const double tol = 20.0 * std::numeric_limits::epsilon() * std::max(1, n); - bool ok_spd = spd_finite && (e_spd <= tol); // siblings uncontaminated & correct - bool ok_bad = bad_poisoned; // poison confined but present - bool ok_info = (info == 0); // GIGO: no early-exit, no error + const bool ok_spd = spd_finite && passes(r_spd); // siblings uncontaminated & correct + const bool ok_bad = bad_poisoned; // poison confined but present + const bool ok_info = (info == 0); // GIGO: no early-exit, no error - std::printf("T=%-6s V=%-2d uplo=%c lay=%c n=%-3d non-SPD lane | siblings:%.1e %s " + std::printf("T=%-6s V=%-2d uplo=%c lay=%c n=%-3d non-SPD lane | siblings:%-5.2g %s " "poison:%s info=%d %s\n", - compact::name, V, uplo, layout, n, e_spd, ok_spd ? "OK" : "FAIL", + compact::name, V, uplo, layout, n, r_spd, ok_spd ? "OK" : "FAIL", ok_bad ? "OK" : "FAIL", info, ok_info ? "OK" : "FAIL"); return (!ok_spd) + (!ok_bad) + (!ok_info); } diff --git a/tests/test_sytrfnp_compact.cpp b/tests/test_sytrfnp_compact.cpp index feb8874..def3c3d 100644 --- a/tests/test_sytrfnp_compact.cpp +++ b/tests/test_sytrfnp_compact.cpp @@ -11,7 +11,8 @@ // reference itself is validated first, against the library: its (L, D) must // reconstruct A through cblas_?trmm / ?gemm (test_reference). // -// Checks per (T, V, uplo, layout), on symmetric *indefinite* batches: +// Checks per (T, V, uplo, layout), on symmetric *indefinite* batches, each a +// test ratio against THRESH (test_compact_util.hpp): // 1. named-triangle factor == scalar reference factor (elementwise) // 2. reconstruction L D L^T == A (lower) / U^T D U == A (upper) // 3. the strictly-opposite triangle of the compact buffer is bit-for-bit @@ -77,12 +78,11 @@ template static void ref_sytf2np(char uplo, MatrixView As) // (L, D) must reconstruct A. W = L D (the columns of the unit-lower L scaled // by D), then W := W L^T by ?trmm (side 'R', unit diagonal); the upper // factor U^T D U is the same product read through the transposed view, as in -// ref_sytf2np itself, and A is symmetric. Gated at -// ||L D L^T - A||_1 / ||A||_1 <= 20 n eps, the design's 7.1 residual. +// ref_sytf2np itself, and A is symmetric. Gated on the test ratio +// ||L D L^T - A|| / (n ||A|| eps), the design's 7.1 residual. template static int test_reference(int n, char uplo) { - const T eps = std::numeric_limits::epsilon(); const bool upper = (uplo == 'U' || uplo == 'u'); std::vector As((size_t)n * n), Fs((size_t)n * n), Ws((size_t)n * n), Rs((size_t)n * n); @@ -96,14 +96,12 @@ template static int test_reference(int n, char uplo) for (int i = 0; i < n; ++i) W(i, j) = (i == j) ? L(j, j) : (i > j) ? L(i, j) * L(j, j) : T(0); // L D tri_apply('R', 'L', 'T', 'U', L, W, R); // (L D) L^T - const double worst = - max_abs_diff(Rs.data(), As.data(), As.size()) / std::max(norm1(A), norm_floor); - const double tol = 20.0 * eps * std::max(1, n); - const bool ok = worst <= tol; - std::printf( - "T=%-6s reference ref_sytf2np uplo=%c n=%-3d | ||LDL^T-A||/||A|| %.1e (%.1e) " - "%s\n", - compact::name, uplo, n, worst, tol, ok ? "OK" : "FAIL"); + const double r = + test_ratio(max_abs_diff(Rs.data(), As.data(), As.size()), n, norm1(A)); + const bool ok = passes(r); + std::printf("T=%-6s reference ref_sytf2np uplo=%c n=%-3d | ||LDL^T-A||/(n ||A|| eps) " + "%-5.2g %s\n", + compact::name, uplo, n, r, verdict(r)); return !ok; } @@ -111,7 +109,6 @@ template static int test_reference(int n, char uplo) template static int run_case(int nm, int n, char uplo, char layout) { - const T eps = std::numeric_limits::epsilon(); const bool rowmajor = (layout == 'R' || layout == 'r'); const bool upper = (uplo == 'U' || uplo == 'u'); @@ -130,12 +127,14 @@ template static int run_case(int nm, int n, char uplo, char lay MatrixBatch Aout(nm, n, n); unpack_compact(Aout, ap.data(), lda, V, rowmajor); - double e_fac = 0, e_rec = 0, e_untouched = 0, a_norm = 1; + double r_fac = 0, r_rec = 0, e_untouched = 0; for (int idx = 0; idx < nm; ++idx) { const auto Ain = A.view(idx), Fac = Aout.view(idx), Ref = Aref.view(idx); - a_norm = std::max(a_norm, norm1(Ain)); - // check 1: named-triangle factor vs scalar reference (elementwise) + const double na = norm1(Ain); + // check 1: named-triangle factor vs scalar reference (elementwise, the + // same operation sequence; relative to ||A||) // check 3: strictly-opposite triangle unchanged from the input A + double e_fac = 0, e_rec = 0; for (int j = 0; j < n; ++j) for (int i = 0; i < n; ++i) { const bool named = upper ? (i <= j) : (i >= j); @@ -145,27 +144,26 @@ template static int run_case(int nm, int n, char uplo, char lay e_untouched = std::max(e_untouched, (double)std::abs(Fac(i, j) - Ain(i, j))); } + r_fac = std::max(r_fac, test_ratio(e_fac, n, na)); - // check 2: reconstruction of A from the named triangle's (D, L|U) + // check 2: reconstruction of A from the named triangle's (D, L|U): + // ||L D L^T - A|| / (n ||A|| eps), dsyt01's form for (int i = 0; i < n; ++i) for (int j = 0; j < n; ++j) e_rec = std::max(e_rec, std::abs(ldlt_reconstruct(Fac, i, j, upper) - (double)Ain(i, j))); + r_rec = std::max(r_rec, test_ratio(e_rec, n, na)); } - const double scale = std::max(1, n); - const double tol_fac = 20.0 * eps * scale; // same op sequence - const double tol_rec = 40.0 * eps * scale * a_norm; // O(n) accumulation in recon - bool ok_f = e_fac <= tol_fac; - bool ok_r = e_rec <= tol_rec; - bool ok_u = e_untouched == 0.0; // must be bit-for-bit unchanged - bool ok_i = (info == 0); - - std::printf("T=%-6s V=%-2d uplo=%c lay=%c nm=%-2d n=%-3d | fac:%.1e %s rec:%.1e %s " - "untouched:%s | info=%d %s\n", - compact::name, V, uplo, layout, nm, n, e_fac, ok_f ? "OK" : "FAIL", - e_rec, ok_r ? "OK" : "FAIL", ok_u ? "OK" : "FAIL", info, - ok_i ? "OK" : "FAIL"); + const bool ok_f = passes(r_fac), ok_r = passes(r_rec); + const bool ok_u = e_untouched == 0.0; // must be bit-for-bit unchanged + const bool ok_i = (info == 0); + + std::printf( + "T=%-6s V=%-2d uplo=%c lay=%c nm=%-2d n=%-3d | fac:%-5.2g %s rec:%-5.2g %s " + "untouched:%s | info=%d %s\n", + compact::name, V, uplo, layout, nm, n, r_fac, verdict(r_fac), r_rec, + verdict(r_rec), ok_u ? "OK" : "FAIL", info, ok_i ? "OK" : "FAIL"); return (!ok_f) + (!ok_r) + (!ok_u) + (!ok_i); } @@ -177,7 +175,6 @@ template static int run_case(int nm, int n, char uplo, char lay template static int run_solve(int nm, int n, int nrhs, char uplo, char layout) { - const T eps = std::numeric_limits::epsilon(); const bool rowmajor = (layout == 'R' || layout == 'r'); const int lda = std::max(1, n), ldb = std::max(1, rowmajor ? nrhs : n); @@ -203,22 +200,20 @@ static int run_solve(int nm, int n, int nrhs, char uplo, char layout) MatrixBatch Xhat(nm, n, nrhs); unpack_compact(Xhat, bp.data(), ldb, V, rowmajor); - const auto [e_fwd, e_res] = solve_errors(A, B, Xhat, X); + // dget02's residual and dget04's forward error (the latter discounted by + // rcond(A), which the indefinite batch's conditioning enters through) + const auto [r_res, r_fwd] = solve_ratios(A, B, Xhat, X); // fused vs two-step, on the raw compact buffers (padded lanes included) const bool fused_same = (ap2 == ap) && (bp2 == bp); - // The residual gate is what the (backward-stable) sweeps control; the - // forward error additionally carries cond(A) of the indefinite batch, so - // its gate gets conditioning headroom. - const double rtol_res = 100.0 * std::max(1, n) * eps; - const double rtol_fwd = 500.0 * std::max(1, n) * eps; - bool ok = (e_fwd <= rtol_fwd) && (e_res <= rtol_res) && fused_same && (info_f == 0) && - (info_s == 0) && (info_v == 0); - std::printf("T=%-6s V=%-2d uplo=%c lay=%c nm=%-2d n=%-3d nrhs=%d solve | fwd:%.1e " - "(%.1e) res:%.1e (%.1e) sysv==trf+trs:%s info=%d/%d/%d %s\n", - compact::name, V, uplo, layout, nm, n, nrhs, e_fwd, rtol_fwd, e_res, - rtol_res, fused_same ? "yes" : "NO", info_f, info_s, info_v, - ok ? "OK" : "FAIL"); + const bool ok = passes(r_res) && passes(r_fwd) && fused_same && (info_f == 0) && + (info_s == 0) && (info_v == 0); + std::printf( + "T=%-6s V=%-2d uplo=%c lay=%c nm=%-2d n=%-3d nrhs=%d solve | res:%-5.2g %s " + "fwd:%-5.2g %s sysv==trf+trs:%s info=%d/%d/%d %s\n", + compact::name, V, uplo, layout, nm, n, nrhs, r_res, verdict(r_res), r_fwd, + verdict(r_fwd), fused_same ? "yes" : "NO", info_f, info_s, info_v, + ok ? "OK" : "FAIL"); return !ok; } @@ -293,10 +288,11 @@ template static int run_zerodiag(char uplo, char layout) MatrixBatch Aout(nm, n, n); unpack_compact(Aout, ap.data(), n, V, rowmajor); - double e = 0; + double r = 0; bool finite = true; for (int idx = 0; idx < nm; ++idx) { const auto Fac = Aout.view(idx), Ref = Aref.view(idx); + double e = 0; for (int j = 0; j < n; ++j) for (int i = 0; i < n; ++i) { const bool named = upper ? (i <= j) : (i >= j); @@ -304,12 +300,12 @@ template static int run_zerodiag(char uplo, char layout) if (!std::isfinite((double)Fac(i, j))) finite = false; e = std::max(e, (double)std::abs(Fac(i, j) - Ref(i, j))); } + r = std::max(r, test_ratio(e, n, norm1(A.view(idx)))); } - const double tol = 20.0 * std::numeric_limits::epsilon() * n; - bool ok = finite && (e <= tol) && (info == 0); - std::printf("T=%-6s V=%-2d uplo=%c lay=%c zero-diagonal A(1,1)=0 | err:%.1e " + const bool ok = finite && passes(r) && (info == 0); + std::printf("T=%-6s V=%-2d uplo=%c lay=%c zero-diagonal A(1,1)=0 | err:%-5.2g " "finite:%s %s\n", - compact::name, V, uplo, layout, e, finite ? "yes" : "NO", + compact::name, V, uplo, layout, r, finite ? "yes" : "NO", ok ? "OK" : "FAIL"); return !ok; } @@ -352,11 +348,12 @@ template static int run_zeropivot(int n, char uplo, char layout MatrixBatch Aout(nm, n, n); unpack_compact(Aout, ap.data(), n, V, rowmajor); - double e_sib = 0; // worst error over the factorable sibling lanes + double r_sib = 0; // worst ratio over the factorable sibling lanes bool sib_finite = true; bool bad_poisoned = false; // the zero-pivot lane must carry Inf/NaN for (int idx = 0; idx < nm; ++idx) { const auto Fac = Aout.view(idx), Ref = Aref.view(idx); + double e_sib = 0; for (int j = 0; j < n; ++j) for (int i = 0; i < n; ++i) { const bool named = upper ? (i <= j) : (i >= j); @@ -370,17 +367,19 @@ template static int run_zeropivot(int n, char uplo, char layout e_sib = std::max(e_sib, (double)std::abs(x - Ref(i, j))); } } + if (idx != badlane) + r_sib = std::max(r_sib, test_ratio(e_sib, n, norm1(A.view(idx)))); } - const double tol = 20.0 * std::numeric_limits::epsilon() * std::max(1, n); - bool ok_sib = sib_finite && (e_sib <= tol); // siblings uncontaminated & correct - bool ok_bad = bad_poisoned; // poison confined but present - bool ok_info = (info == 0); // GIGO: no early-exit, no error + const bool ok_sib = sib_finite && passes(r_sib); // siblings uncontaminated & correct + const bool ok_bad = bad_poisoned; // poison confined but present + const bool ok_info = (info == 0); // GIGO: no early-exit, no error - std::printf("T=%-6s V=%-2d uplo=%c lay=%c n=%-3d zero-pivot lane | siblings:%.1e %s " - "poison:%s info=%d %s\n", - compact::name, V, uplo, layout, n, e_sib, ok_sib ? "OK" : "FAIL", - ok_bad ? "OK" : "FAIL", info, ok_info ? "OK" : "FAIL"); + std::printf( + "T=%-6s V=%-2d uplo=%c lay=%c n=%-3d zero-pivot lane | siblings:%-5.2g %s " + "poison:%s info=%d %s\n", + compact::name, V, uplo, layout, n, r_sib, ok_sib ? "OK" : "FAIL", + ok_bad ? "OK" : "FAIL", info, ok_info ? "OK" : "FAIL"); return (!ok_sib) + (!ok_bad) + (!ok_info); } diff --git a/tests/test_trsm_compact.cpp b/tests/test_trsm_compact.cpp index e4c84aa..02288fe 100644 --- a/tests/test_trsm_compact.cpp +++ b/tests/test_trsm_compact.cpp @@ -9,7 +9,8 @@ // 1. C API argument validation -- exercises the LAPACK/BLAS-style info = -j // contract of the portable entry points. // 2. Numerical correctness over side / uplo / transa / diag, precisions and -// interleave widths (column-major, the tuned path): forward error vs the +// interleave widths (column-major, the tuned path), as test ratios +// against THRESH (test_compact_util.hpp): dget04's forward error vs the // library's solve, plus the solve's own residual ||op(A) X - alpha B|| // formed with the library's triangular multiply (so a bug shared by the // reference and the kernel cannot pass unseen), over the suite's shapes @@ -43,7 +44,6 @@ static int run_case(char side, char uplo, char transa, char diag, int nm, int m, { const bool left = (side == 'L'); const int s = left ? m : n; - const T eps = std::numeric_limits::epsilon(); const T alpha = T(0.5) + frand(); // a non-trivial, non-zero scalar // triangular A (s x s), column-major: random in the referenced triangle, @@ -73,35 +73,40 @@ static int run_case(char side, char uplo, char transa, char diag, int nm, int m, MatrixBatch Bout(nm, m, n); unpack_compact(Bout, bp.data(), ldb, V); - // Two gates: (1) forward error vs the library's solve, and (2) the - // solve's own defining residual ||op(A) X - alpha B||, formed with the - // library's triangular multiply -- so a bug shared by ref_trsm and the - // kernel cannot slip through (the ?trsm analogue of the reconstruction / - // round-trip identities the geqrf/potrf/ormqr self-tests check). - double worst_fwd = 0, worst_res = 0; - std::vector Rs((size_t)m * n), aBs((size_t)m * n); + // Two gates: (1) the forward error vs the library's solve -- dget04's, + // discounted by rcond of the operator actually applied, op(A) with the + // unit diagonal honored (E, formed by the library's triangular multiply + // of the identity) -- and (2) the solve's own defining residual + // ||op(A) X - alpha B|| / (s ||op(A)|| ||X|| eps), dtrt02's, formed with + // the same multiply -- so a bug shared by ref_trsm and the kernel cannot + // slip through (the ?trsm analogue of the reconstruction / round-trip + // identities the geqrf/potrf/ormqr self-tests check). + double r_fwd = 0, r_res = 0; + std::vector Rs((size_t)m * n), aBs((size_t)m * n), Is((size_t)s * s, T(0)), + Es((size_t)s * s); const auto R = mat_view(Rs.data(), m, n), aB = mat_view(aBs.data(), m, n); + const auto I = mat_view(Is.data(), s, s), E = mat_view(Es.data(), s, s); + for (int d = 0; d < s; ++d) + I(d, d) = T(1); for (int idx = 0; idx < nm; ++idx) { - worst_fwd = - std::max(worst_fwd, max_abs_diff(Bout[idx], Xref[idx], (size_t)m * n) / - std::max(norm1(Xref.view(idx)), norm_floor)); + tri_apply('L', uplo, transa, diag, A.view(idx), I, E); // E = op(A) + r_fwd = std::max( + r_fwd, forward_ratio(max_abs_diff(Bout[idx], Xref[idx], (size_t)m * n), + norm1(Xref.view(idx)), rcond1(E))); tri_apply(side, uplo, transa, diag, A.view(idx), Bout.view(idx), R); for (size_t e = 0; e < (size_t)m * n; ++e) aBs[e] = alpha * B[idx][e]; - worst_res = - std::max(worst_res, max_abs_diff(Rs.data(), aBs.data(), (size_t)m * n) / - std::max(norm1(aB), norm_floor)); + r_res = std::max(r_res, + test_ratio(max_abs_diff(Rs.data(), aBs.data(), (size_t)m * n), + s, norm1(E) * norm1(Bout.view(idx)))); } - const double worst = std::max(worst_fwd, worst_res); - // Backward-stable triangular solve on a diagonal-boosted A: both the forward - // error and the residual stay near working precision. - const double rtol = 1e3 * s * (double)eps; - const bool ok = (worst <= rtol); + const bool ok = passes(r_fwd) && passes(r_res); std::printf(" T=%-6s V=%-2d side=%c uplo=%c tr=%c diag=%c nm=%-2d m=%-3d n=%-3d | " - "fwd %.1e res %.1e (rtol %.1e) info=%d %s\n", - compact::name, V, side, uplo, transa, diag, nm, m, n, worst_fwd, - worst_res, rtol, info, (ok && info == 0) ? "OK" : "FAIL"); + "fwd %-5.2g %s res %-5.2g %s info=%d %s\n", + compact::name, V, side, uplo, transa, diag, nm, m, n, r_fwd, + verdict(r_fwd), r_res, verdict(r_res), info, + (ok && info == 0) ? "OK" : "FAIL"); return (info != 0) + !ok; } From a379fc857f777ffe37a2d7cb69722787851af98c Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 15 Sep 2026 00:06:03 +0000 Subject: [PATCH 3/6] tests: LAPACK's exact ratio forms, and the PO/SY/TR paths' invariants The ratio refactor kept LAPACK's denominators but not its numerators: max-abs entries where dqrt01 / dpot01 / dsyt01 take the 1-norm of the residual matrix, and a whole-matrix ||X||_1 with an order factor where dget02 / dpot02 / dtrt02 are per right-hand side without one -- up to n times more lenient than the checkers they cite. Every ratio now has the checker's exact form (test_compact_util.hpp): test_ratio over diff_norm1 / orth_norm1 (the 1-norm of R - Q^T A, L L^T - A, Q^T Q - I, ...); residual_ratio, dget02's max_j ||b_j - A x_j||_1 / (||A||_1 ||x_j||_1 eps); forward_ratio, dget04's max_j ||x_j - x_j^true||_inf / ||x_j^true||_inf * rcond / eps, on views rather than a scalar error. solve_ratios forms the residual explicitly and hands it to those. Observed ratios rose accordingly, to at most 4.9 against THRESH = 30. The PO, SY and TR paths were audited against dchkpo / dchksy / dchktr. Present already: dpot01 / dsyt01 reconstruction, dpot02 residual, dget04 forward error, dtrt02 residual, the fused drivers bit-for-bit against their two-step forms, plus the lane-isolation and untouched-triangle contracts LAPACK has no analogue for. Not applicable: dpot03 / dtrt01 (no inverse routine), dpocon / dtrcon (no condition estimator), the refinement bounds (no ?rfs). Missing, and added: - trsm had no exact-solution forward error, only agreement with cblas_?trsm: the right-hand side is now B = op(A) X for a random X (dlarhs), so dget04 runs against the exact alpha X as well as against the library's solution, and dtrt02's residual reads side 'R' as the transposed system its rows are the right-hand sides of. - No check tied potrs / sytrsnp to the trsm sweeps they are built from: run_compose in both suites solves the same system through the public ?trsm_compact (two non-unit sweeps for potrs; two unit sweeps around D^-1 for sytrsnp) and gates dget04 of the one solution against the other. potrs comes out bit-identical to its two sweeps in every configuration (reported, not gated); sytrsnp's ratio is ~0.03. CLAUDE.md's convention entry names the three forms and the cross-routine checks. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_014ZZ76AgYwQ3vtuyZ46tyMX --- .claude/CLAUDE.md | 32 ++++++---- tests/test_compact_util.hpp | 92 ++++++++++++++++++++++++--- tests/test_gels_compact.cpp | 32 +++++----- tests/test_geqrf_compact.cpp | 36 +++++------ tests/test_lapack_util.hpp | 28 ++++----- tests/test_orgqr_compact.cpp | 12 ++-- tests/test_ormqr_compact.cpp | 17 ++--- tests/test_potrf_compact.cpp | 85 ++++++++++++++++++++++--- tests/test_sytrfnp_compact.cpp | 110 ++++++++++++++++++++++++++++----- tests/test_trsm_compact.cpp | 67 +++++++++++--------- 10 files changed, 367 insertions(+), 144 deletions(-) diff --git a/.claude/CLAUDE.md b/.claude/CLAUDE.md index 4a05b49..4203f80 100644 --- a/.claude/CLAUDE.md +++ b/.claude/CLAUDE.md @@ -252,17 +252,27 @@ workspace contract, or the benchmarks' threading. written as a `JB`-templated block helper with `for (c < JB)` loops the compiler unrolls, not as hand-expanded `w0..w3` copies. - **Test ratios, one threshold.** The portable suites gate every numerical - check the way LAPACK's TESTING suites do: a dimensionless test ratio -- the - error over what backward stability allows, `(order) * (operand norm) * eps` - (`test_ratio` in `tests/test_compact_util.hpp`; `forward_ratio` is - dget04's, a forward error discounted by `rcond` from `rcond1`, LAPACKE's - `?getrf` + `?gecon` on the dense, unpacked input) -- compared against the - single `THRESH = 30` of LAPACK's `dtest.in`. A correct kernel's ratios are - O(1) (observed ≤ 2.2 across the suites) in either precision; each case prints - them, so a near miss is visible. Add a check as a ratio, never as a hand-picked - `k * eps * n` tolerance; only structural contracts (bit-for-bit reproduction, - untouched storage, identity padding lanes) are exact. The MKL suites still - carry their own relative gates. + check the way LAPACK's TESTING suites do: a dimensionless test ratio in the + exact form of LAPACK's own checker, compared against the single + `THRESH = 30` of LAPACK's `dtest.in`. The three forms are in + `tests/test_compact_util.hpp`: `test_ratio` for a factorization residual + (dqrt01 / dpot01 / dsyt01: the *1-norm* of `R - Q^T A`, `L L^T - A`, ... + over `(order) * ||A||_1 * eps` -- `diff_norm1`, `orth_norm1` form the + numerators), `residual_ratio` for a solve (dget02 / dpot02 / dtrt02: per + right-hand side, `||b_j - A x_j||_1 / (||A||_1 ||x_j||_1 eps)`, no order + factor), and `forward_ratio` for a forward error (dget04: per right-hand + side, `||x_j - x_j^true||_inf / ||x_j^true||_inf`, discounted by `rcond` from + `rcond1`, LAPACKE's `?getrf` + `?gecon` on the dense, unpacked input -- + condition numbers are never taken from the packed format). A correct + kernel's ratios are O(1) (observed ≤ 5 across the suites) in either + precision; each case prints them, so a near miss is visible. Add a check as + one of these ratios, never as a hand-picked `k * eps * n` tolerance or a + max-abs entry; only structural contracts (bit-for-bit reproduction, + untouched storage, identity padding lanes) are exact. Related routines are + also checked against each other, as LAPACK's paths do: `posv == potrf + + potrs` and `sysvnp == sytrfnp + sytrsnp` bit-for-bit, `potrs` against the two + `trsm` sweeps it is built from, `sytrsnp` against `trsm`, `D^-1`, `trsm`. + The MKL suites still carry their own relative gates. - **Argument checking.** The MKL-style API (`cbk_*`) skips validation like MKL's own compact routines (`info` is a scalar, `0` on success). The portable C API (`cbk.h`) validates LAPACK-style, returning `-j` for a bad j-th diff --git a/tests/test_compact_util.hpp b/tests/test_compact_util.hpp index f98798f..4012a14 100644 --- a/tests/test_compact_util.hpp +++ b/tests/test_compact_util.hpp @@ -88,13 +88,20 @@ constexpr double norm_floor = 1e-300; // reproduction, storage the routine must not touch, the identity in the // padding lanes -- stay exact. // -// test_ratio(err, dim, norm) err / (max(1, dim) * norm * eps) -// the residual form (dqrt01, dpot01, dget02): err the worst entry of -// R - Q^T A, L L^T - A, B - A X, ...; dim the order it accumulates over; -// norm the operand(s) it is relative to (1 for an orthogonal Q) -// forward_ratio(err, normx, rcond) err * rcond / (normx * eps) -// the forward-error form (dget04): err = ||x - x_true||, discounted by +// The three forms are LAPACK's own checkers', norms included: +// test_ratio(err, dim, norm) err / (max(1, dim) * norm * eps) +// the factorization form (dqrt01, dpot01, dsyt01, dqrt16): err the +// 1-norm of the residual matrix R - Q^T A, L L^T - A, ... (diff_norm1, +// orth_norm1); dim the order it accumulates over; norm the operand(s) it +// is relative to (1 for an orthogonal Q) +// residual_ratio(anorm, R, X) max_j ||r_j||_1 / (anorm ||x_j||_1 eps) +// the solve form (dget02, dpot02, dtrt02): per right-hand side, R the +// residual B - op(A) X, anorm = ||op(A)||_1 -- no order factor +// forward_ratio(Xhat, X, rcond) max_j (||xhat_j - x_j||_inf / ||x_j||_inf) rcond / eps +// the forward-error form (dget04): per right-hand side, discounted by // rcond = 1 / (||A|| ||A^-1||), the amplification a solve is allowed +// A zero solution column with a nonzero error is 1 / eps (a failure), as +// in LAPACK; an empty system is 0. constexpr double THRESH = 30.0; template double test_ratio(double err, int dim, double norm = 1.0) @@ -103,9 +110,44 @@ template double test_ratio(double err, int dim, double norm = 1.0) (double)std::numeric_limits::epsilon()); } -template double forward_ratio(double err, double normx, double rcond) +template double residual_ratio(double anorm, Rv R, Xv X) { - return test_ratio(err * rcond, 1, normx); + using T = elem_t; + assert(R.cols == X.cols); + const double eps = std::numeric_limits::epsilon(); + if (R.rows == 0 || R.cols == 0) return 0; + if (anorm <= 0) return 1.0 / eps; + double r = 0; + for (int j = 0; j < R.cols; ++j) { + double bnorm = 0, xnorm = 0; + for (int i = 0; i < R.rows; ++i) + bnorm += std::abs((double)R(i, j)); + for (int i = 0; i < X.rows; ++i) + xnorm += std::abs((double)X(i, j)); + r = std::max(r, xnorm <= 0 ? 1.0 : (bnorm / anorm) / xnorm); + } + return r / eps; +} + +template double forward_ratio(Xv Xhat, Xa X, double rcond) +{ + using T = elem_t; + assert(Xhat.rows == X.rows && Xhat.cols == X.cols); + const double eps = std::numeric_limits::epsilon(); + double r = 0; + for (int j = 0; j < X.cols; ++j) { + double diff = 0, xnorm = 0; + for (int i = 0; i < X.rows; ++i) { + diff = std::max(diff, std::abs((double)Xhat(i, j) - (double)X(i, j))); + xnorm = std::max(xnorm, std::abs((double)X(i, j))); + } + if (xnorm <= 0) { + if (diff > 0) r = std::max(r, 1.0); + } + else + r = std::max(r, (diff / xnorm) * rcond); + } + return r / eps; } inline bool passes(double ratio) @@ -147,9 +189,41 @@ template double norm1(Mv M) return mx; } +// ||A - B||_1 of two same-shaped views, in double: the residual norm LAPACK's +// checkers take of a reconstruction (dqrt01, dpot01, dsyt01) or a comparison. +template double diff_norm1(Av A, Bv B) +{ + assert(A.rows == B.rows && A.cols == B.cols); + double mx = 0; + for (int j = 0; j < A.cols; ++j) { + double s = 0; + for (int i = 0; i < A.rows; ++i) + s += std::abs((double)A(i, j) - (double)B(i, j)); + mx = std::max(mx, s); + } + return mx; +} + +// ||Q^T Q - I||_1, formed in double: dqrt02's orthogonality residual. +template double orth_norm1(Qv Q) +{ + double mx = 0; + for (int j = 0; j < Q.cols; ++j) { + double s = 0; + for (int i = 0; i < Q.cols; ++i) { + double d = 0; + for (int l = 0; l < Q.rows; ++l) + d += (double)Q(l, i) * Q(l, j); + s += std::abs(d - (i == j ? 1.0 : 0.0)); + } + mx = std::max(mx, s); + } + return mx; +} + // Deviation of Q's columns from orthonormality, max |(Q^T Q - I)(i,j)| over // the upper triangle (Q^T Q is symmetric), accumulated in double so the gate -// does not inherit float rounding. Shared by the orgqr suites. +// does not inherit float rounding. The MKL orgqr suite's gate. template double orth_error(Qv Q) { double e = 0; diff --git a/tests/test_gels_compact.cpp b/tests/test_gels_compact.cpp index 5cbe453..f722ad8 100644 --- a/tests/test_gels_compact.cpp +++ b/tests/test_gels_compact.cpp @@ -84,17 +84,19 @@ int run_case(char layout, char trans, int nm, int m, int n, int nrhs) double r_x = 0, r_h = 0, r_t = 0; for (int idx = 0; idx < nm; ++idx) { const double na = norm1(A.view(idx)); - r_x = std::max(r_x, test_ratio(max_abs_diff(Bout[idx], Bref[idx], B.stride()), - p, norm1(Bref.view(idx)))); - r_h = std::max( - r_h, test_ratio(max_abs_diff(Aout[idx], Aref[idx], A.stride()), p, na)); - r_t = std::max( - r_t, test_ratio(max_abs_diff(tau_out[idx], tau_ref[idx], (size_t)q), p)); + r_x = std::max(r_x, test_ratio(diff_norm1(Bout.view(idx), Bref.view(idx)), p, + norm1(Bref.view(idx)))); + r_h = std::max(r_h, + test_ratio(diff_norm1(Aout.view(idx), Aref.view(idx)), p, na)); + r_t = std::max(r_t, test_ratio(diff_norm1(mat_view(tau_out[idx], q, 1), + mat_view(tau_ref[idx], q, 1)), + p)); } // check 3: the defining properties, formed without the reference: r_prop - // is the normal equations (least squares) or the residual (minimum norm), - // relative to ||op(A)|| ||B|| (dqrt17's scale), r_prop2 the + // is the normal equations (least squares: dqrt17's ||op(A)^T r|| / + // (||op(A)|| ||B|| max(m,n,nrhs) eps)) or the residual (minimum norm: + // dqrt16's ||B - op(A) X|| / (max(m,n) ||op(A)|| ||X|| eps)), r_prop2 the // residual-sum-of-squares rows or the Gram-formed solution double r_prop = 0, r_prop2 = 0; std::vector Rs((size_t)rows_op * nrhs), Ss((size_t)cols_op * nrhs); @@ -112,9 +114,7 @@ int run_case(char layout, char trans, int nm, int m, int n, int nrhs) for (int i = 0; i < rows_op; ++i) R(i, j) = Bin(i, j) - R(i, j); matmul(Aop.transposed(), R, S); - for (size_t e = 0; e < Ss.size(); ++e) - r_prop = - std::max(r_prop, test_ratio(std::abs((double)Ss[e]), p, scale)); + r_prop = std::max(r_prop, test_ratio(norm1(S), std::max(p, nrhs), scale)); // rows cols_op..rows_op-1 of B hold the residual: squared column // norms == ||r_j||^2 (relative to ||b_j||^2: a square system has no // residual rows and a zero residual). Not at min(m,n) = 0: that is @@ -137,9 +137,8 @@ int run_case(char layout, char trans, int nm, int m, int n, int nrhs) // op(A) X = B for (int j = 0; j < nrhs; ++j) for (int i = 0; i < rows_op; ++i) - r_prop = std::max( - r_prop, - test_ratio(std::abs((double)R(i, j) - Bin(i, j)), p, scale)); + R(i, j) = Bin(i, j) - R(i, j); + r_prop = std::max(r_prop, test_ratio(norm1(R), p, norm1(Aop) * norm1(X))); // minimum norm the other way: X = op(A)^T Z, G Z = B with the Gram // matrix G = op(A) op(A)^T (rows_op x rows_op, SPD), solved by the // LAPACKE QR references. The two solutions differ by what the @@ -157,10 +156,7 @@ int run_case(char layout, char trans, int nm, int m, int n, int nrhs) ref_ormqr('T', rows_op, G, tg.data(), Z); ref_trsm_upper(G, Z); matmul(Aop.transposed(), Z, Xmn); - r_prop2 = - std::max(r_prop2, forward_ratio(max_abs_diff(Xmns.data(), X.data, - (size_t)cols_op * nrhs), - norm1(X), rcond_g)); + r_prop2 = std::max(r_prop2, forward_ratio(Xmn, X, rcond_g)); } } diff --git a/tests/test_geqrf_compact.cpp b/tests/test_geqrf_compact.cpp index 50efced..dba32d8 100644 --- a/tests/test_geqrf_compact.cpp +++ b/tests/test_geqrf_compact.cpp @@ -67,10 +67,11 @@ template static int run_case(int nm, int m, int n) double r_h = 0, r_t = 0, r_rec = 0; for (int idx = 0; idx < nm; ++idx) { const double na = norm1(A.view(idx)); - r_h = std::max( - r_h, test_ratio(max_abs_diff(Aout[idx], Aref[idx], (size_t)m * n), m, na)); - r_t = std::max( - r_t, test_ratio(max_abs_diff(tau_out[idx], tau_ref[idx], (size_t)k), m)); + r_h = std::max(r_h, + test_ratio(diff_norm1(Aout.view(idx), Aref.view(idx)), m, na)); + r_t = std::max(r_t, test_ratio(diff_norm1(mat_view(tau_out[idx], k, 1), + mat_view(tau_ref[idx], k, 1)), + m)); std::vector Recs((size_t)m * n, T(0)); // start from R const auto Rec = mat_view(Recs.data(), m, n); @@ -78,9 +79,7 @@ template static int run_case(int nm, int m, int n) for (int i = 0; i <= std::min(j, k - 1); ++i) Rec(i, j) = Aout(idx, i, j); ref_ormqr('N', k, Aout.view(idx), tau_out[idx], Rec); - r_rec = std::max( - r_rec, - test_ratio(max_abs_diff(Recs.data(), A[idx], (size_t)m * n), m, na)); + r_rec = std::max(r_rec, test_ratio(diff_norm1(Rec, A.view(idx)), m, na)); } // check 3: solve A X = B with the produced reflectors (square only), using @@ -102,9 +101,8 @@ template static int run_case(int nm, int m, int n) r_solve = 0; for (int idx = 0; idx < nm; ++idx) { ref_trsm_upper(Aout.view(idx), Bo.view(idx)); - r_solve = std::max(r_solve, forward_ratio(max_abs_diff(Bo[idx], Xs.data(), - (size_t)n * nrhs), - norm1(X), rcond1(A.view(idx)))); + r_solve = + std::max(r_solve, forward_ratio(Bo.view(idx), X, rcond1(A.view(idx)))); } } @@ -182,14 +180,12 @@ static int run_invariants(int nm, int m, int n, double cond, Structure structure ref_orgqr(k, Q, tau[idx]); // residual R - Q^T A, R = triu(H) (k x n): dqrt01's - // ||R - Q^T A|| / (m ||A|| eps) + // ||R - Q^T A||_1 / (m ||A||_1 eps) matmul(Q.transposed(), Am, QtA); - double resid = 0; for (int j = 0; j < n; ++j) for (int i = 0; i < k; ++i) - resid = std::max( - resid, (double)std::abs(((i <= j) ? Hm(i, j) : T(0)) - QtA(i, j))); - r_res = std::max(r_res, test_ratio(resid, m, norm1(Am))); + QtA(i, j) = ((i <= j) ? Hm(i, j) : T(0)) - QtA(i, j); + r_res = std::max(r_res, test_ratio(norm1(QtA), m, norm1(Am))); // orthogonality Q^T Q - I: dqrt01's ||I - Q^T Q|| / (m eps) matmul(Q.transposed(), Q, QtQ); @@ -206,8 +202,9 @@ static int run_invariants(int nm, int m, int n, double cond, Structure structure if (structure == DENSE) { copy_matrix(Am, Href); ref_geqrf(Href, tauref.data()); - const double el = std::max(max_abs_diff(H[idx], Hrefs.data(), A.stride()), - max_abs_diff(tau[idx], tauref.data(), (size_t)k)); + const double el = + std::max(diff_norm1(Hm, Href), diff_norm1(mat_view(tau[idx], k, 1), + mat_view(tauref.data(), k, 1))); r_el = std::max(r_el, test_ratio(el, n, norm1(Am))); } } @@ -296,9 +293,8 @@ template static int test_underflow() for (int i = 0; i <= std::min(j, n - 1); ++i) Rec(i, j) = Aout(idx, i, j); ref_ormqr('N', n, Aout.view(idx), tau[idx], Rec); - r_rec = std::max(r_rec, - test_ratio(max_abs_diff(Recs.data(), A[idx], (size_t)m * n), - m, norm1(A.view(idx)))); + r_rec = std::max( + r_rec, test_ratio(diff_norm1(Rec, A.view(idx)), m, norm1(A.view(idx)))); } fails += !passes(r_rec); std::printf("T=%-6s V=%-2d underflow column: tau=0, diag kept, body 0, finite %s | " diff --git a/tests/test_lapack_util.hpp b/tests/test_lapack_util.hpp index ca7ec1a..9cf1d66 100644 --- a/tests/test_lapack_util.hpp +++ b/tests/test_lapack_util.hpp @@ -426,10 +426,11 @@ template double rcond1(Av A) } // The end-to-end solve gate every solving suite closes with, as test ratios: -// dget02's residual ||B - A Xhat|| / (n ||A|| ||Xhat|| eps), formed by ?gemm -// so a solver bug cannot hide behind the factorization that produced Xhat, -// and dget04's forward error ||Xhat - X|| rcond(A) / (||X|| eps) against the -// known X. The worst over the batch of each. +// dget02's / dpot02's residual, per right-hand side, ||b_j - A xhat_j||_1 / +// (||A||_1 ||xhat_j||_1 eps), the residual formed by ?gemm so a solver bug +// cannot hide behind the factorization that produced Xhat; and dget04's +// forward error against the known X, discounted by rcond(A). The worst over +// the batch of each. struct SolveRatios { double res, fwd; }; @@ -440,18 +441,17 @@ SolveRatios solve_ratios(const MatrixBatch &A, const MatrixBatch &B, { assert(B.rows() == Xhat.rows() && B.cols() == Xhat.cols() && X.rows == Xhat.rows() && X.cols == Xhat.cols()); - const int n = A.rows(); - const size_t sB = Xhat.stride(); - std::vector AXs(sB); - const auto AX = mat_view(AXs.data(), Xhat.rows(), Xhat.cols()); - const double nX = norm1(X); + std::vector Rs(Xhat.stride()); + const auto R = mat_view(Rs.data(), Xhat.rows(), Xhat.cols()); SolveRatios r{0, 0}; for (int v = 0; v < Xhat.count(); ++v) { - matmul(A.view(v), Xhat.view(v), AX); - r.res = std::max(r.res, test_ratio(max_abs_diff(AXs.data(), B[v], sB), n, - norm1(A.view(v)) * norm1(Xhat.view(v)))); - r.fwd = std::max(r.fwd, forward_ratio(max_abs_diff(Xhat[v], X.data, sB), nX, - rcond1(A.view(v)))); + const auto Av = A.view(v), Bv = B.view(v), Xv = Xhat.view(v); + matmul(Av, Xv, R); // R := B - A Xhat + for (int j = 0; j < R.cols; ++j) + for (int i = 0; i < R.rows; ++i) + R(i, j) = Bv(i, j) - R(i, j); + r.res = std::max(r.res, residual_ratio(norm1(Av), R, Xv)); + r.fwd = std::max(r.fwd, forward_ratio(Xv, X, rcond1(Av))); } return r; } diff --git a/tests/test_orgqr_compact.cpp b/tests/test_orgqr_compact.cpp index 998165f..780ef8f 100644 --- a/tests/test_orgqr_compact.cpp +++ b/tests/test_orgqr_compact.cpp @@ -81,13 +81,12 @@ static int run_case(int nm, int m, int n, int k, bool rowmajor = false) * to 1, the norm of an orthogonal matrix) */ double r1 = 0; for (int kk = 0; kk < nm; ++kk) - r1 = std::max(r1, - test_ratio(max_abs_diff(Qout[kk], Qref[kk], Qout.stride()), m)); + r1 = std::max(r1, test_ratio(diff_norm1(Qout.view(kk), Qref.view(kk)), m)); - /* check 2: Q^T Q = I, formed densely (dqrt02's ||I - Q^T Q|| / (m eps)) */ + /* check 2: Q^T Q = I, formed densely (dqrt02's ||I - Q^T Q||_1 / (m eps)) */ double r2 = 0; for (int kk = 0; kk < nm; ++kk) - r2 = std::max(r2, test_ratio(orth_error(Qout.view(kk)), m)); + r2 = std::max(r2, test_ratio(orth_norm1(Qout.view(kk)), m)); /* check 3: Q(:, 0:k-1) R = A0, R = triu of the factorization, relative * to ||A0|| */ @@ -101,9 +100,8 @@ static int run_case(int nm, int m, int n, int k, bool rowmajor = false) for (int i = 0; i < k; ++i) R(i, j) = (i <= j) ? Afac(kk, i, j) : T(0); matmul(leading(Qout.view(kk), m, k), R, Rec); - r3 = std::max(r3, - test_ratio(max_abs_diff(Recs.data(), A0[kk], (size_t)m * k), - m, norm1(A0.view(kk)))); + r3 = std::max( + r3, test_ratio(diff_norm1(Rec, A0.view(kk)), m, norm1(A0.view(kk)))); } } diff --git a/tests/test_ormqr_compact.cpp b/tests/test_ormqr_compact.cpp index ea5a415..bf04f3a 100644 --- a/tests/test_ormqr_compact.cpp +++ b/tests/test_ormqr_compact.cpp @@ -65,17 +65,14 @@ template static int run_case(int nm, int m, int nrhs) unpack_compact(Bout, bp.data(), ld, V); double r1 = 0; for (int kk = 0; kk < nm; ++kk) - r1 = - std::max(r1, test_ratio(max_abs_diff(Bout[kk], Bref[kk], (size_t)m * nrhs), - m, norm1(Bref.view(kk)))); + r1 = std::max(r1, test_ratio(diff_norm1(Bout.view(kk), Bref.view(kk)), m, + norm1(Bref.view(kk)))); /* check 2: solve recovers X (dget04, discounted by rcond(A)) */ double r2 = 0; for (int kk = 0; kk < nm; ++kk) { ref_trsm_upper(Afac.view(kk), Bout.view(kk)); - r2 = std::max( - r2, forward_ratio(max_abs_diff(Bout[kk], Xs.data(), (size_t)m * nrhs), - norm1(X), rcond1(A.view(kk)))); + r2 = std::max(r2, forward_ratio(Bout.view(kk), X, rcond1(A.view(kk)))); } /* check 3: 'N' undoes 'T', relative to ||B||_1 */ @@ -83,8 +80,8 @@ template static int run_case(int nm, int m, int nrhs) unpack_compact(Bout, bp.data(), ld, V); double r3 = 0; for (int kk = 0; kk < nm; ++kk) - r3 = std::max(r3, test_ratio(max_abs_diff(Bout[kk], B[kk], (size_t)m * nrhs), - m, norm1(B.view(kk)))); + r3 = std::max(r3, test_ratio(diff_norm1(Bout.view(kk), B.view(kk)), m, + norm1(B.view(kk)))); const bool ok1 = passes(r1), ok2 = passes(r2), ok3 = passes(r3); std::printf( @@ -139,9 +136,7 @@ template static int run_case_pivoted(int nm, int m, int nrhs) for (int j = 0; j < nrhs; ++j) for (int i = 0; i < m; ++i) x(jpvt[kk][i], j) = y(i, j); - r = std::max( - r, forward_ratio(max_abs_diff(xs.data(), Xs.data(), (size_t)m * nrhs), - norm1(X), rcond1(A.view(kk)))); + r = std::max(r, forward_ratio(x, X, rcond1(A.view(kk)))); } const bool ok = passes(r); diff --git a/tests/test_potrf_compact.cpp b/tests/test_potrf_compact.cpp index 6411d25..2327e3f 100644 --- a/tests/test_potrf_compact.cpp +++ b/tests/test_potrf_compact.cpp @@ -20,6 +20,8 @@ // ?posv_compact reproduces that factor and X bit-for-bit // 5. nrhs = 0: ?posv_compact still factors (LAPACK ?posv), bit-identical to // ?potrf_compact, with a 1-element dummy bp (never referenced) +// 6. ?potrs_compact == the two ?trsm_compact sweeps it is built from +// (dget04 of the one solution against the other) // over the suite's orders and the small-dimension sweep (small_dims: every // order from 0 to 5 through the factorization checks, n x nrhs through the // solve pipeline, the empty operand included), plus LAPACK-style argument @@ -66,19 +68,23 @@ static int run_case(int nm, int n, char uplo, char layout, double cond = 0.0) const auto F = mat_view(Fs.data(), n, n), Rec = mat_view(Recs.data(), n, n); for (int idx = 0; idx < nm; ++idx) { const auto Fac = Aout.view(idx), Ref = Aref.view(idx), Ain = A.view(idx); - // check 1: named-triangle factor vs LAPACKE_?potrf, relative to the - // reference factor's L1 norm (the unique SPD factor: a sharp signal) + // check 1: named-triangle factor vs LAPACKE_?potrf, the 1-norm of the + // difference over the triangle relative to the reference factor's + // (the unique SPD factor: a sharp signal) // check 3: strictly-opposite triangle unchanged from the input A double el = 0; - for (int j = 0; j < n; ++j) + for (int j = 0; j < n; ++j) { + double cs = 0; for (int i = 0; i < n; ++i) { const bool named = upper ? (i <= j) : (i >= j); if (named) - el = std::max(el, (double)std::abs(Fac(i, j) - Ref(i, j))); + cs += std::abs((double)Fac(i, j) - Ref(i, j)); else e_untouched = std::max(e_untouched, (double)std::abs(Fac(i, j) - Ain(i, j))); } + el = std::max(el, cs); + } r_fac = std::max(r_fac, test_ratio(el, n, norm1(Ref))); // check 2: reconstruction of A from the named triangle's factor, by @@ -92,9 +98,7 @@ static int run_case(int nm, int n, char uplo, char layout, double cond = 0.0) tri_apply('R', 'L', 'T', 'N', Fac, F, Rec); else tri_apply('L', 'U', 'T', 'N', Fac, F, Rec); - r_rec = std::max(r_rec, - test_ratio(max_abs_diff(Recs.data(), A[idx], (size_t)n * n), - n, norm1(Ain))); + r_rec = std::max(r_rec, test_ratio(diff_norm1(Rec, Ain), n, norm1(Ain))); } const bool ok_f = passes(r_fac), ok_r = passes(r_rec); @@ -157,6 +161,59 @@ static int run_solve(int nm, int n, int nrhs, char uplo, char layout) return !ok; } +// ------------- potrs vs the two trsm sweeps it is built from ---------- +// ?potrs_compact is two ?trsm_compact sweeps on the factor (lower: L z = B, +// L^T X = z; upper: U^T z = B, U X = z), so the public trsm entry point must +// solve the same system from the same factor -- the consistency between +// related routines LAPACK's drivers check through their shared ratios. +// Gated as dget04's forward error of the one solution against the other +// (rcond(A)); whether the two are bit-identical is reported, not gated. + +template +static int run_compose(int nm, int n, int nrhs, char uplo, char layout) +{ + const bool rowmajor = (layout == 'R' || layout == 'r'); + const bool upper = (uplo == 'U' || uplo == 'u'); + const int lda = std::max(1, n), ldb = std::max(1, rowmajor ? nrhs : n); + + MatrixBatch A(nm, n, n), B(nm, n, nrhs); + const std::vector Xs = known_solution(n, nrhs); + const auto X = mat_view(Xs.data(), n, nrhs); + for (int idx = 0; idx < nm; ++idx) { + gen_spd(A.view(idx)); + matmul(A.view(idx), X, B.view(idx)); + } + std::vector ap = pack_compact(A, lda, V, rowmajor); + std::vector bp = pack_compact(B, ldb, V, rowmajor); + std::vector bp2 = bp; + + int bad = (compact::potrf(layout, uplo, n, ap.data(), lda, V, nm) != 0); + bad += (compact::potrs(layout, uplo, n, nrhs, ap.data(), lda, bp.data(), ldb, V, + nm) != 0); + // the same solve as two trsm sweeps on the factor + const char t1 = upper ? 'T' : 'N', t2 = upper ? 'N' : 'T'; + bad += (compact::trsm(layout, 'L', uplo, t1, 'N', n, nrhs, T(1), ap.data(), lda, + bp2.data(), ldb, V, nm) != 0); + bad += (compact::trsm(layout, 'L', uplo, t2, 'N', n, nrhs, T(1), ap.data(), lda, + bp2.data(), ldb, V, nm) != 0); + + MatrixBatch X1(nm, n, nrhs), X2(nm, n, nrhs); + unpack_compact(X1, bp.data(), ldb, V, rowmajor); + unpack_compact(X2, bp2.data(), ldb, V, rowmajor); + double r = 0; + for (int idx = 0; idx < nm; ++idx) + r = std::max(r, forward_ratio(X1.view(idx), X2.view(idx), rcond1(A.view(idx)))); + const bool same = (bp == bp2); + + const bool ok = passes(r) && (bad == 0); + std::printf( + "T=%-6s V=%-2d uplo=%c lay=%c nm=%-2d n=%-3d nrhs=%d potrs vs trsm+trsm | " + "%-5.2g %s bit-identical:%s %s\n", + compact::name, V, uplo, layout, nm, n, nrhs, r, verdict(r), + same ? "yes" : "no", ok ? "OK" : "FAIL"); + return !ok; +} + // ------------- nrhs = 0: the fused driver still factors -------------- // LAPACK ?posv calls ?potrf unconditionally -- the nrhs = 0 quick return is // ?potrs's -- so the fused driver must factor ap even with no right-hand @@ -228,7 +285,8 @@ template static int run_nonspd(int n, char uplo, char layout) bool bad_poisoned = false; // the non-SPD lane must carry NaN/Inf for (int idx = 0; idx < nm; ++idx) { double el = 0; - for (int j = 0; j < n; ++j) + for (int j = 0; j < n; ++j) { + double cs = 0; for (int i = 0; i < n; ++i) { const bool named = upper ? (i <= j) : (i >= j); if (!named) continue; @@ -238,8 +296,10 @@ template static int run_nonspd(int n, char uplo, char layout) continue; } spd_finite = spd_finite && std::isfinite((double)x); - el = std::max(el, (double)std::abs(x - Aref(idx, i, j))); + cs += std::abs((double)x - Aref(idx, i, j)); } + el = std::max(el, cs); + } if (idx != badlane) // vs LAPACKE_?potrf, relative to its factor's norm r_spd = std::max(r_spd, test_ratio(el, n, norm1(Aref.view(idx)))); } @@ -306,6 +366,13 @@ int main() fails += run_solve(16, 24, 3, u, l); } + // potrs against the two trsm sweeps it is built from, over uplo/layout. + for (char u : uplos) + for (char l : lays) { + fails += run_compose(8, 30, 5, u, l); + fails += run_compose(11, 24, 3, u, l); // padded partial group + } + // The small-dimension sweep (small_dims): every order over uplo x layout, // then n x nrhs through the solve pipeline, a padded group each. for (int n : small_dims) diff --git a/tests/test_sytrfnp_compact.cpp b/tests/test_sytrfnp_compact.cpp index def3c3d..3d76c55 100644 --- a/tests/test_sytrfnp_compact.cpp +++ b/tests/test_sytrfnp_compact.cpp @@ -21,6 +21,8 @@ // ?sysvnp_compact reproduces that factor and X bit-for-bit // 5. nrhs = 0: ?sysvnp_compact still factors (LAPACK ?sysv), bit-identical // to ?sytrfnp_compact, with a 1-element dummy bp (never referenced) +// 6. ?sytrsnp_compact == the unit-diagonal ?trsm_compact sweeps around +// D^-1 it is built from (dget04 of the one solution against the other) // over the suite's orders and the small-dimension sweep (small_dims: every // order from 0 to 5 through the factorization checks, n x nrhs through the // solve pipeline, the empty operand included), plus: @@ -96,8 +98,7 @@ template static int test_reference(int n, char uplo) for (int i = 0; i < n; ++i) W(i, j) = (i == j) ? L(j, j) : (i > j) ? L(i, j) * L(j, j) : T(0); // L D tri_apply('R', 'L', 'T', 'U', L, W, R); // (L D) L^T - const double r = - test_ratio(max_abs_diff(Rs.data(), As.data(), As.size()), n, norm1(A)); + const double r = test_ratio(diff_norm1(R, A), n, norm1(A)); const bool ok = passes(r); std::printf("T=%-6s reference ref_sytf2np uplo=%c n=%-3d | ||LDL^T-A||/(n ||A|| eps) " "%-5.2g %s\n", @@ -131,27 +132,29 @@ template static int run_case(int nm, int n, char uplo, char lay for (int idx = 0; idx < nm; ++idx) { const auto Ain = A.view(idx), Fac = Aout.view(idx), Ref = Aref.view(idx); const double na = norm1(Ain); - // check 1: named-triangle factor vs scalar reference (elementwise, the - // same operation sequence; relative to ||A||) + // check 1: named-triangle factor vs scalar reference (the same + // operation sequence; the 1-norm of the difference over the triangle, + // relative to ||A||) + // check 2: reconstruction of A from the named triangle's (D, L|U): + // dsyt01's ||L D L^T - A||_1 / (n ||A||_1 eps) // check 3: strictly-opposite triangle unchanged from the input A double e_fac = 0, e_rec = 0; - for (int j = 0; j < n; ++j) + for (int j = 0; j < n; ++j) { + double cs_fac = 0, cs_rec = 0; for (int i = 0; i < n; ++i) { const bool named = upper ? (i <= j) : (i >= j); if (named) - e_fac = std::max(e_fac, (double)std::abs(Fac(i, j) - Ref(i, j))); + cs_fac += std::abs((double)Fac(i, j) - Ref(i, j)); else e_untouched = std::max(e_untouched, (double)std::abs(Fac(i, j) - Ain(i, j))); + cs_rec += + std::abs(ldlt_reconstruct(Fac, i, j, upper) - (double)Ain(i, j)); } + e_fac = std::max(e_fac, cs_fac); + e_rec = std::max(e_rec, cs_rec); + } r_fac = std::max(r_fac, test_ratio(e_fac, n, na)); - - // check 2: reconstruction of A from the named triangle's (D, L|U): - // ||L D L^T - A|| / (n ||A|| eps), dsyt01's form - for (int i = 0; i < n; ++i) - for (int j = 0; j < n; ++j) - e_rec = std::max(e_rec, std::abs(ldlt_reconstruct(Fac, i, j, upper) - - (double)Ain(i, j))); r_rec = std::max(r_rec, test_ratio(e_rec, n, na)); } @@ -217,6 +220,66 @@ static int run_solve(int nm, int n, int nrhs, char uplo, char layout) return !ok; } +// ------------- sytrsnp vs the trsm sweeps it is built from ------------ +// ?sytrsnp_compact is two unit-diagonal ?trsm_compact sweeps around a +// division by D (lower: L z = B, z := D^-1 z, L^T X = z; upper: U^T z = B, +// z := D^-1 z, U X = z), so the public trsm entry point, with the diagonal +// scaling done densely in between, must solve the same system from the same +// factor. Gated as dget04's forward error of the one solution against the +// other (rcond(A)); the two are not expected bit-identical (the kernel's +// scaling is its own), so only the ratio is gated. + +template +static int run_compose(int nm, int n, int nrhs, char uplo, char layout) +{ + const bool rowmajor = (layout == 'R' || layout == 'r'); + const bool upper = (uplo == 'U' || uplo == 'u'); + const int lda = std::max(1, n), ldb = std::max(1, rowmajor ? nrhs : n); + + MatrixBatch A(nm, n, n), B(nm, n, nrhs); + const std::vector Xs = known_solution(n, nrhs); + const auto X = mat_view(Xs.data(), n, nrhs); + for (int idx = 0; idx < nm; ++idx) { + gen_sym_ldlt(A.view(idx)); + matmul(A.view(idx), X, B.view(idx)); + } + std::vector ap = pack_compact(A, lda, V, rowmajor); + std::vector bp = pack_compact(B, ldb, V, rowmajor); + std::vector bp2 = bp; + + int bad = (compact::sytrfnp(layout, uplo, n, ap.data(), lda, V, nm) != 0); + bad += (compact::sytrsnp(layout, uplo, n, nrhs, ap.data(), lda, bp.data(), ldb, V, + nm) != 0); + // the same solve as two unit-diagonal trsm sweeps around D^-1 + const char t1 = upper ? 'T' : 'N', t2 = upper ? 'N' : 'T'; + bad += (compact::trsm(layout, 'L', uplo, t1, 'U', n, nrhs, T(1), ap.data(), lda, + bp2.data(), ldb, V, nm) != 0); + MatrixBatch F(nm, n, n), Z(nm, n, nrhs); + unpack_compact(F, ap.data(), lda, V, rowmajor); + unpack_compact(Z, bp2.data(), ldb, V, rowmajor); + for (int idx = 0; idx < nm; ++idx) + for (int j = 0; j < nrhs; ++j) + for (int i = 0; i < n; ++i) + Z(idx, i, j) /= F(idx, i, i); + pack_compact(Z, bp2.data(), ldb, V, rowmajor); + bad += (compact::trsm(layout, 'L', uplo, t2, 'U', n, nrhs, T(1), ap.data(), lda, + bp2.data(), ldb, V, nm) != 0); + + MatrixBatch X1(nm, n, nrhs), X2(nm, n, nrhs); + unpack_compact(X1, bp.data(), ldb, V, rowmajor); + unpack_compact(X2, bp2.data(), ldb, V, rowmajor); + double r = 0; + for (int idx = 0; idx < nm; ++idx) + r = std::max(r, forward_ratio(X1.view(idx), X2.view(idx), rcond1(A.view(idx)))); + + const bool ok = passes(r) && (bad == 0); + std::printf("T=%-6s V=%-2d uplo=%c lay=%c nm=%-2d n=%-3d nrhs=%d sytrsnp vs " + "trsm+D+trsm | %-5.2g %s %s\n", + compact::name, V, uplo, layout, nm, n, nrhs, r, verdict(r), + ok ? "OK" : "FAIL"); + return !ok; +} + // ------------- nrhs = 0: the fused driver still factors -------------- // LAPACK ?sysv calls ?sytrf unconditionally -- the nrhs = 0 quick return is // ?sytrs's -- so the fused driver must factor ap even with no right-hand @@ -293,13 +356,16 @@ template static int run_zerodiag(char uplo, char layout) for (int idx = 0; idx < nm; ++idx) { const auto Fac = Aout.view(idx), Ref = Aref.view(idx); double e = 0; - for (int j = 0; j < n; ++j) + for (int j = 0; j < n; ++j) { + double cs = 0; for (int i = 0; i < n; ++i) { const bool named = upper ? (i <= j) : (i >= j); if (!named) continue; if (!std::isfinite((double)Fac(i, j))) finite = false; - e = std::max(e, (double)std::abs(Fac(i, j) - Ref(i, j))); + cs += std::abs((double)Fac(i, j) - Ref(i, j)); } + e = std::max(e, cs); + } r = std::max(r, test_ratio(e, n, norm1(A.view(idx)))); } const bool ok = finite && passes(r) && (info == 0); @@ -354,7 +420,8 @@ template static int run_zeropivot(int n, char uplo, char layout for (int idx = 0; idx < nm; ++idx) { const auto Fac = Aout.view(idx), Ref = Aref.view(idx); double e_sib = 0; - for (int j = 0; j < n; ++j) + for (int j = 0; j < n; ++j) { + double cs = 0; for (int i = 0; i < n; ++i) { const bool named = upper ? (i <= j) : (i >= j); if (!named) continue; @@ -364,9 +431,11 @@ template static int run_zeropivot(int n, char uplo, char layout } else { if (!std::isfinite((double)x)) sib_finite = false; - e_sib = std::max(e_sib, (double)std::abs(x - Ref(i, j))); + cs += std::abs((double)x - Ref(i, j)); } } + e_sib = std::max(e_sib, cs); + } if (idx != badlane) r_sib = std::max(r_sib, test_ratio(e_sib, n, norm1(A.view(idx)))); } @@ -438,6 +507,13 @@ int main() fails += run_solve(16, 24, 3, u, l); } + // sytrsnp against the trsm sweeps it is built from, over uplo/layout. + for (char u : uplos) + for (char l : lays) { + fails += run_compose(8, 30, 5, u, l); + fails += run_compose(11, 24, 3, u, l); // padded partial group + } + // The small-dimension sweep (small_dims): every order over uplo x layout, // then n x nrhs through the solve pipeline, a padded group each. for (int n : small_dims) diff --git a/tests/test_trsm_compact.cpp b/tests/test_trsm_compact.cpp index 02288fe..37916ea 100644 --- a/tests/test_trsm_compact.cpp +++ b/tests/test_trsm_compact.cpp @@ -10,10 +10,12 @@ // contract of the portable entry points. // 2. Numerical correctness over side / uplo / transa / diag, precisions and // interleave widths (column-major, the tuned path), as test ratios -// against THRESH (test_compact_util.hpp): dget04's forward error vs the -// library's solve, plus the solve's own residual ||op(A) X - alpha B|| -// formed with the library's triangular multiply (so a bug shared by the -// reference and the kernel cannot pass unseen), over the suite's shapes +// against THRESH (test_compact_util.hpp), LAPACK's TR path's: dget04's +// forward error vs the exact solution (B = op(A) X for a random X, as +// dlarhs) and vs the library's solve, plus dtrt02's residual +// ||alpha B - op(A) X|| formed with the library's triangular multiply +// (so a bug shared by the reference and the kernel cannot pass unseen), +// over the suite's shapes // and the small-dimension cross product (small_dims: m x n from 0 to 5, // the empty operand and the 1/2/3-column route included). // @@ -53,11 +55,15 @@ static int run_case(char side, char uplo, char transa, char diag, int nm, int m, for (int idx = 0; idx < nm; ++idx) gen_tri(A.view(idx), up); - // random B (m x n) and its reference solution - MatrixBatch B(nm, m, n), Xref(nm, m, n); + // a random solution X (dlarhs's), the right-hand side B = op(A) X (side + // 'L') or X op(A) ('R') by the library's triangular multiply -- so the + // solve of op(A) X = alpha B has the exact solution alpha X -- and the + // library's own solution of the same system + MatrixBatch X(nm, m, n), B(nm, m, n), Xref(nm, m, n); for (int idx = 0; idx < nm; ++idx) { for (size_t e = 0; e < (size_t)m * n; ++e) - B[idx][e] = frand(); + X[idx][e] = frand(); + tri_apply(side, uplo, transa, diag, A.view(idx), X.view(idx), B.view(idx)); std::copy(B[idx], B[idx] + (size_t)m * n, Xref[idx]); ref_trsm(side, uplo, transa, diag, alpha, A.view(idx), Xref.view(idx)); } @@ -73,39 +79,44 @@ static int run_case(char side, char uplo, char transa, char diag, int nm, int m, MatrixBatch Bout(nm, m, n); unpack_compact(Bout, bp.data(), ldb, V); - // Two gates: (1) the forward error vs the library's solve -- dget04's, - // discounted by rcond of the operator actually applied, op(A) with the - // unit diagonal honored (E, formed by the library's triangular multiply - // of the identity) -- and (2) the solve's own defining residual - // ||op(A) X - alpha B|| / (s ||op(A)|| ||X|| eps), dtrt02's, formed with - // the same multiply -- so a bug shared by ref_trsm and the kernel cannot - // slip through (the ?trsm analogue of the reconstruction / round-trip - // identities the geqrf/potrf/ormqr self-tests check). - double r_fwd = 0, r_res = 0; - std::vector Rs((size_t)m * n), aBs((size_t)m * n), Is((size_t)s * s, T(0)), + // Three gates, LAPACK's TR path's: dget04's forward error against the + // exact solution alpha X and against the library's solution, both + // discounted by rcond of the operator actually applied -- E = op(A) with + // the unit diagonal honored, formed by the library's triangular multiply + // of the identity -- and dtrt02's residual, per right-hand side, + // ||alpha b_j - op(A) x_j||_1 / (||op(A)||_1 ||x_j||_1 eps), formed with + // the same multiply, so a bug shared by ref_trsm and the kernel cannot + // slip through. Side 'R' is the transposed system op(A)^T X^T = alpha B^T: + // its right-hand sides are the rows, read through transposed views, and + // its operator op(A)^T. + double r_fwd = 0, r_ref = 0, r_res = 0; + std::vector Rs((size_t)m * n), aXs((size_t)m * n), Is((size_t)s * s, T(0)), Es((size_t)s * s); - const auto R = mat_view(Rs.data(), m, n), aB = mat_view(aBs.data(), m, n); + const auto R = mat_view(Rs.data(), m, n), aX = mat_view(aXs.data(), m, n); const auto I = mat_view(Is.data(), s, s), E = mat_view(Es.data(), s, s); for (int d = 0; d < s; ++d) I(d, d) = T(1); + const auto sys = [left](MatrixView M) { return left ? M : M.transposed(); }; for (int idx = 0; idx < nm; ++idx) { tri_apply('L', uplo, transa, diag, A.view(idx), I, E); // E = op(A) - r_fwd = std::max( - r_fwd, forward_ratio(max_abs_diff(Bout[idx], Xref[idx], (size_t)m * n), - norm1(Xref.view(idx)), rcond1(E))); + const auto Eop = sys(E); + const double rc = rcond1(Eop), anorm = norm1(Eop); + for (size_t e = 0; e < (size_t)m * n; ++e) + aXs[e] = alpha * X[idx][e]; + r_fwd = std::max(r_fwd, forward_ratio(sys(Bout.view(idx)), sys(aX), rc)); + r_ref = + std::max(r_ref, forward_ratio(sys(Bout.view(idx)), sys(Xref.view(idx)), rc)); tri_apply(side, uplo, transa, diag, A.view(idx), Bout.view(idx), R); for (size_t e = 0; e < (size_t)m * n; ++e) - aBs[e] = alpha * B[idx][e]; - r_res = std::max(r_res, - test_ratio(max_abs_diff(Rs.data(), aBs.data(), (size_t)m * n), - s, norm1(E) * norm1(Bout.view(idx)))); + Rs[e] = alpha * B[idx][e] - Rs[e]; // alpha B - op(A) Xhat + r_res = std::max(r_res, residual_ratio(anorm, sys(R), sys(Bout.view(idx)))); } - const bool ok = passes(r_fwd) && passes(r_res); + const bool ok = passes(r_fwd) && passes(r_ref) && passes(r_res); std::printf(" T=%-6s V=%-2d side=%c uplo=%c tr=%c diag=%c nm=%-2d m=%-3d n=%-3d | " - "fwd %-5.2g %s res %-5.2g %s info=%d %s\n", + "fwd %-5.2g %s ref %-5.2g %s res %-5.2g %s info=%d %s\n", compact::name, V, side, uplo, transa, diag, nm, m, n, r_fwd, - verdict(r_fwd), r_res, verdict(r_res), info, + verdict(r_fwd), r_ref, verdict(r_ref), r_res, verdict(r_res), info, (ok && info == 0) ? "OK" : "FAIL"); return (info != 0) + !ok; From 33ddb9085b7a9a37d67b060c69e198c734438386 Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 15 Sep 2026 00:18:42 +0000 Subject: [PATCH 4/6] api: give the portable ?ormqr_compact layout, side and trans The kernel (cbk_ormqr_compact.hpp) covers every layout, side and trans in one strided body -- side 'R' is the kernel over the transposed view -- and the MKL-style cbk_?ormqr_compact exposes all of it, but the portable ?ormqr_compact hardcoded side 'L', column-major, so the portable surface (the one every stack builds) could only apply Q^T from the left, and its suite could only test that. ?ormqr_compact now takes (layout, side, trans, m, n, k, ap, ldap, taup, cp, ldcp, V, nm), LAPACK ?ormqr's argument order behind the layout, like the other portable entry points and the MKL-style one. Validation in argument order: -1 layout, -2 side, -3 trans, -4 m, -5 n, -6 k (k <= nq, nq = m for side 'L' and n for 'R', the order of Q and the rows of A), -8 ldap (>= nq column-major, >= k row-major), -11 ldcp, -12 V, -13 nm. The compat and compact dispatchers, ref_ormqr (a side-aware overload, the side-'L' one forwarding to it) and the geqrf suite's call follow. tests/test_ormqr_compact.cpp runs the full layout x side x trans product on every stack: vs LAPACKE_?ormqr, vs ?gemm with the Q that ?orgqr forms explicitly -- LAPACK's dqrt03, ||CC - op(Q) C||_1 / (nq ||C||_1 eps) -- and the round trip op(Q)^T op(Q) C = C, at the suite's shapes (k < nq and a padded group among them), the small-dimension cross product with LAPACK's K values on the tuned route and k = nq on the transposed one, and the solve pipelines (?geqr2's and ?geqp3's reflectors) as before. Observed ratios: at most 1.6, the tiny shapes; 0.06-0.6 at the representative ones. PLANS.md and CLAUDE.md record the generalized surface and the dqrt03 check. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_014ZZ76AgYwQ3vtuyZ46tyMX --- .claude/CLAUDE.md | 4 +- docs/PLANS.md | 5 +- include/cbk.h | 40 +++--- src/cbk.cpp | 40 +++--- tests/test_compact_util.hpp | 6 +- tests/test_geqrf_compact.cpp | 4 +- tests/test_lapack_util.hpp | 23 ++- tests/test_ormqr_compact.cpp | 262 +++++++++++++++++++++++------------ 8 files changed, 250 insertions(+), 134 deletions(-) diff --git a/.claude/CLAUDE.md b/.claude/CLAUDE.md index 4203f80..46a2f71 100644 --- a/.claude/CLAUDE.md +++ b/.claude/CLAUDE.md @@ -271,7 +271,9 @@ workspace contract, or the benchmarks' threading. untouched storage, identity padding lanes) are exact. Related routines are also checked against each other, as LAPACK's paths do: `posv == potrf + potrs` and `sysvnp == sytrfnp + sytrsnp` bit-for-bit, `potrs` against the two - `trsm` sweeps it is built from, `sytrsnp` against `trsm`, `D^-1`, `trsm`. + `trsm` sweeps it is built from, `sytrsnp` against `trsm`, `D^-1`, `trsm`, + `ormqr` in every layout / side / trans against the explicit Q of `orgqr` + (dqrt03). The MKL suites still carry their own relative gates. - **Argument checking.** The MKL-style API (`cbk_*`) skips validation like MKL's own compact routines (`info` is a scalar, `0` on success). The portable C diff --git a/docs/PLANS.md b/docs/PLANS.md index bea94e0..63cec84 100644 --- a/docs/PLANS.md +++ b/docs/PLANS.md @@ -21,7 +21,10 @@ listed in `examples.md`. ## ormqr - **Done:** design 2-6 and 8.1; `side = 'R'` is the kernel over the - transposed view, `trans = 'C'` folds to `'T'`. + transposed view, `trans = 'C'` folds to `'T'`. The portable `?ormqr_compact` + takes layout, side and trans like the MKL-style entry point, and its suite + runs the layout x side x trans product against LAPACKE `?ormqr` and the + explicit Q of `?orgqr` (LAPACK's dqrt03) on every stack. - **Gap:** the stress structures of design 7.3 (`cond` knob, banded / row-scaled / clustered-scale inputs) are only partly covered; `geqrf`'s suite has the rank-deficient and near-collinear cases. The LAPACK-style diff --git a/include/cbk.h b/include/cbk.h index 65a2b42..0041d80 100644 --- a/include/cbk.h +++ b/include/cbk.h @@ -8,7 +8,7 @@ * interleave width V. * * dgeqrf_compact / sgeqrf_compact -- QR factorization A = Q R - * dormqr_compact / sormqr_compact -- apply Q or Q^T from the left, B := op(Q) B + * dormqr_compact / sormqr_compact -- apply Q or Q^T, C := op(Q) C or C op(Q) * dorgqr_compact / sorgqr_compact -- form the explicit Q of a QR, in place * dpotrf_compact / spotrf_compact -- Cholesky factorization A = L L^T or U^T U * dpotrs_compact / spotrs_compact -- solve A X = B from a Cholesky factor @@ -94,27 +94,33 @@ int dgeqrf_compact(char layout, int m, int n, double *ap, int ldap, double *taup int sgeqrf_compact(char layout, int m, int n, float *ap, int ldap, float *taup, int V, int nm); -/* Apply Q (or Q^T) of a compact QR to a compact RHS block from the left, - * B := op(Q) B -- the reflector-application (side='L') step between a compact - * QR factorization (?geqrf_compact) and a triangular solve. - * trans 'T' (Q^T B, the solve case) or 'N' (Q B) - * m, nrhs rows of B (and A); columns of B - * k number of reflectors (min(m,n) of the factorization) - * ap compact reflectors from ?geqrf_compact, (ldap, k) per matrix - * ldap compact leading dimension of A (>= m) +/* Apply Q (or Q^T) of a compact QR to a compact matrix, from either side: + * C := op(Q) C (side='L') or C := C op(Q) (side='R'), op(Q) = Q ('N') or Q^T + * ('T') -- LAPACK ?ormqr, one matrix per compact lane, in either layout. + * Q = H(0) H(1) ... H(k-1) is the product of the k reflectors ?geqrf_compact + * left in ap and taup; it is nq x nq with nq = m (side='L') or n (side='R'). + * layout 'C' column-major or 'R' row-major (A and C alike) + * side 'L' (C := op(Q) C) or 'R' (C := C op(Q)) + * trans 'T' (Q^T, the solve case) or 'N' (Q) + * m, n rows and columns of C + * k number of reflectors (0 <= k <= nq) + * ap compact reflectors from ?geqrf_compact, (nq x k) per matrix + * ldap compact leading dimension of A (>= nq column-major, >= k row-major) * taup compact tau (k per matrix, ld = k) - * bp compact B (m x nrhs), overwritten with op(Q) B - * ldbp compact leading dimension of B (>= m) + * cp compact C (m x n), overwritten with op(Q) C or C op(Q) + * ldcp compact leading dimension of C (>= m column-major, >= n row-major) * V, nm interleave width; total number of matrices (padded last group) * Returns 0, or -j for an illegal j-th argument: - * -1 trans -2 m (<0) -3 nrhs (<0) -4 k (<0 or >m) -6 ldap - * -9 ldbp -10 V (not 2/4/8/16) -11 nm (<0) + * -1 layout -2 side -3 trans -4 m (<0) -5 n (<0) -6 k (<0 or >nq) + * -8 ldap -11 ldcp -12 V (not 2/4/8/16) -13 nm (<0) */ -int dormqr_compact(char trans, int m, int nrhs, int k, const double *ap, int ldap, - const double *taup, double *bp, int ldbp, int V, int nm); +int dormqr_compact(char layout, char side, char trans, int m, int n, int k, + const double *ap, int ldap, const double *taup, double *cp, int ldcp, + int V, int nm); -int sormqr_compact(char trans, int m, int nrhs, int k, const float *ap, int ldap, - const float *taup, float *bp, int ldbp, int V, int nm); +int sormqr_compact(char layout, char side, char trans, int m, int n, int k, + const float *ap, int ldap, const float *taup, float *cp, int ldcp, + int V, int nm); /* Generate the explicit orthogonal factor of a compact QR: the first n columns * of Q = H(0) H(1) ... H(k-1) (m >= n >= k), formed in place over the diff --git a/src/cbk.cpp b/src/cbk.cpp index 1eb77b6..fca2cc5 100644 --- a/src/cbk.cpp +++ b/src/cbk.cpp @@ -78,23 +78,28 @@ int geqrf(char layout, int m, int n, T *ap, int ldap, T *taup, int V, int nm) } template -int ormqr(char trans, int m, int nrhs, int k, const T *ap, int ldap, const T *taup, T *bp, - int ldbp, int V, int nm) +int ormqr(char layout, char side, char trans, int m, int n, int k, const T *ap, int ldap, + const T *taup, T *cp, int ldcp, int V, int nm) { - if (!opt(trans, 'T') && !opt(trans, 'N')) return -1; - if (m < 0) return -2; - if (nrhs < 0) return -3; - if (k < 0 || k > m) return -4; - if (ldap < max1(m)) return -6; - if (ldbp < max1(m)) return -9; - if (!vlen_ok(V)) return -10; - if (nm < 0) return -11; - if (m == 0 || nrhs == 0 || k == 0 || nm == 0) return 0; - assert(ap != nullptr && taup != nullptr && bp != nullptr); + const bool col = opt(layout, 'C'), row = opt(layout, 'R'); + const bool left = opt(side, 'L'), right = opt(side, 'R'); + const int nq = left ? m : n; /* the order of Q, and the rows of A */ + if (!col && !row) return -1; + if (!left && !right) return -2; + if (!opt(trans, 'T') && !opt(trans, 'N')) return -3; + if (m < 0) return -4; + if (n < 0) return -5; + if (k < 0 || k > nq) return -6; + if (ldap < max1(row ? k : nq)) return -8; + if (ldcp < max1(row ? n : m)) return -11; + if (!vlen_ok(V)) return -12; + if (nm < 0) return -13; + if (m == 0 || n == 0 || k == 0 || nm == 0) return 0; + assert(ap != nullptr && taup != nullptr && cp != nullptr); for_vlen(V, [&](auto v) { - cbk::detail::ormqr_compact(true, false, trans, m, nrhs, k, - ap, ldap, taup, bp, ldbp, nm); + cbk::detail::ormqr_compact(left, row, trans, m, n, k, ap, + ldap, taup, cp, ldcp, nm); }); return 0; } @@ -326,10 +331,11 @@ int gels(char layout, char trans, int m, int n, int nrhs, T *ap, int ldap, T *bp { \ return geqrf(layout, m, n, ap, ldap, taup, V, nm); \ } \ - int p##ormqr_compact(char trans, int m, int nrhs, int k, const T *ap, int ldap, \ - const T *taup, T *bp, int ldbp, int V, int nm) \ + int p##ormqr_compact(char layout, char side, char trans, int m, int n, int k, \ + const T *ap, int ldap, const T *taup, T *cp, int ldcp, int V, \ + int nm) \ { \ - return ormqr(trans, m, nrhs, k, ap, ldap, taup, bp, ldbp, V, nm); \ + return ormqr(layout, side, trans, m, n, k, ap, ldap, taup, cp, ldcp, V, nm); \ } \ int p##orgqr_compact(char layout, int m, int n, int k, T *ap, int ldap, \ const T *taup, int V, int nm) \ diff --git a/tests/test_compact_util.hpp b/tests/test_compact_util.hpp index 4012a14..d935b11 100644 --- a/tests/test_compact_util.hpp +++ b/tests/test_compact_util.hpp @@ -274,9 +274,9 @@ template <> struct compact { static constexpr const char *name = label; \ static int geqrf(char lay, int m, int n, T *a, int ld, T *tau, int V, int nm) \ { return p##geqrf_compact(lay, m, n, a, ld, tau, V, nm); } \ - static int ormqr(char tr, int m, int nrhs, int k, const T *a, int lda, const T *tau, \ - T *b, int ldb, int V, int nm) \ - { return p##ormqr_compact(tr, m, nrhs, k, a, lda, tau, b, ldb, V, nm); } \ + static int ormqr(char lay, char si, char tr, int m, int n, int k, const T *a, \ + int lda, const T *tau, T *c, int ldc, int V, int nm) \ + { return p##ormqr_compact(lay, si, tr, m, n, k, a, lda, tau, c, ldc, V, nm); } \ static int orgqr(char lay, int m, int n, int k, T *a, int lda, const T *tau, \ int V, int nm) \ { return p##orgqr_compact(lay, m, n, k, a, lda, tau, V, nm); } \ diff --git a/tests/test_geqrf_compact.cpp b/tests/test_geqrf_compact.cpp index dba32d8..3eb9b81 100644 --- a/tests/test_geqrf_compact.cpp +++ b/tests/test_geqrf_compact.cpp @@ -94,8 +94,8 @@ template static int run_case(int nm, int m, int n) for (int idx = 0; idx < nm; ++idx) matmul(A.view(idx), X, B.view(idx)); /* B = A X */ std::vector bp = pack_compact(B, lda, V); - compact::ormqr('T', n, nrhs, k, ap.data(), lda, tp.data(), bp.data(), lda, V, - nm); + compact::ormqr('C', 'L', 'T', n, nrhs, k, ap.data(), lda, tp.data(), bp.data(), + lda, V, nm); MatrixBatch Bo(nm, n, nrhs); unpack_compact(Bo, bp.data(), lda, V); r_solve = 0; diff --git a/tests/test_lapack_util.hpp b/tests/test_lapack_util.hpp index 9cf1d66..219799f 100644 --- a/tests/test_lapack_util.hpp +++ b/tests/test_lapack_util.hpp @@ -240,16 +240,25 @@ template void ref_geqrf(MatrixView A, T *tau) lapack::geqrf(lapack_layout(A), A.rows, A.cols, A.data, lapack_ld(A), tau); } -// ?ormqr, side='L': B := Q^T B (trans 'T') or Q B ('N') from (H, tau), with -// k reflectors read from A (m x k or wider) and applied to B (m x nrhs). +// ?ormqr: C := op(Q) C (side 'L') or C op(Q) ('R'), op(Q) = Q ('N') or Q^T +// ('T'), from (H, tau) with k reflectors read from A (nq x k or wider, nq = m +// for 'L' and n for 'R') and applied to C (m x n). +template +void ref_ormqr(char side, char trans, int k, Av A, const T *tau, MatrixView C) +{ + const bool left = (side == 'L' || side == 'l'); + assert(A.rows == (left ? C.rows : C.cols) && k <= std::min(A.rows, A.cols)); + if (is_empty(C) || k == 0) return; + const Staged Cs(C, is_rowmajor(A)); + lapack::ormqr(lapack_layout(A), side, trans, C.rows, C.cols, k, A.data, + lapack_ld(A), tau, Cs.data(), Cs.ld()); +} + +// The common case, side 'L': B := Q^T B (trans 'T') or Q B ('N'). template void ref_ormqr(char trans, int k, Av A, const T *tau, MatrixView B) { - assert(A.rows == B.rows && k <= std::min(A.rows, A.cols)); - if (is_empty(B) || k == 0) return; - const Staged Bs(B, is_rowmajor(A)); - lapack::ormqr(lapack_layout(A), 'L', trans, B.rows, B.cols, k, A.data, - lapack_ld(A), tau, Bs.data(), Bs.ld()); + ref_ormqr('L', trans, k, A, tau, B); } // ?orgqr: generate the first n columns of Q = H(0)..H(k-1) in place over the diff --git a/tests/test_ormqr_compact.cpp b/tests/test_ormqr_compact.cpp index bf04f3a..4aaf31f 100644 --- a/tests/test_ormqr_compact.cpp +++ b/tests/test_ormqr_compact.cpp @@ -1,22 +1,28 @@ /* test_ormqr_compact.cpp * - * Validation of the templated compact ormqr against LAPACKE - * (test_lapack_util.hpp): the reflectors come from LAPACKE_?geqr2 (or - * LAPACKE_?geqp3 for the pivoted pipeline), the reference application from - * LAPACKE_?ormqr, the back substitution from cblas_?trsm. + * Validation of the templated compact ormqr (dormqr_compact / sormqr_compact) + * against LAPACKE (test_lapack_util.hpp): the reflectors come from + * LAPACKE_?geqr2 (or LAPACKE_?geqp3 for the pivoted pipeline), the reference + * application from LAPACKE_?ormqr, the explicit Q from LAPACKE_?orgqr, the + * back substitution from cblas_?trsm. * - * Test design: X(:,j) = j+1 (ones, twos, threes, ...), B = A*X, so the - * b columns are scaled row sums of A; the solve must recover X. - * - * Checks per (T, V), each a test ratio against THRESH (test_compact_util.hpp): - * 1. compact Q^T B == LAPACKE_?ormqr Q^T B, relative to ||Q^T B||_1 (the - * design's 7.1 gate: two backward-stable applications of the same - * reflectors, LAPACK's blocked above its crossover) - * 2. back substitution recovers X (dget04's forward error, rcond(A)) - * 3. applying 'N' after 'T' recovers the original B (Q Q^T = I) + * Checks per (T, V, layout, side, trans, m, n, k), each a test ratio against + * THRESH (test_compact_util.hpp), LAPACK's QR path's: + * 1. compact op(Q) C == LAPACKE_?ormqr op(Q) C, relative to its norm (two + * backward-stable applications of the same reflectors, LAPACK's blocked + * above its crossover) + * 2. dqrt03: op(Q) C (side 'L') or C op(Q) ('R') against ?gemm with the Q + * that ?orgqr forms explicitly, ||CC - op(Q) C||_1 / (nq ||C||_1 eps), + * nq the order of Q -- the four side/trans combinations in both layouts + * 3. applying op(Q)^T after op(Q) recovers the original C (Q Q^T = I) + * and, in the solve pipeline (side 'L', trans 'T', column-major, the tuned + * route), X(:,j) = j+1, B = A X: + * 4. Q^T B then back substitution recovers X (dget04's forward error, + * rcond(A)), from ?geqr2's reflectors and, through a jpvt + * back-permutation, from ?geqp3's * over the suite's shapes and the small-dimension cross product (small_dims: - * m x nrhs from 0 to 5, the empty operand and the 1/2/3-column right-hand - * side blocks included). + * m x n from 0 to 5 with LAPACK's K values, the empty operand included), plus + * LAPACK-style argument validation of the C API. * * Assisted-by: Claude:claude-fable-5 Claude:claude-opus-4-8 */ @@ -30,66 +36,119 @@ #include #include -#include "test_lapack_util.hpp" // compact, ref_geqr2/orm2r/geqp3/trsm_upper, pack/unpack +#include "test_lapack_util.hpp" // compact, ref_geqr2/ormqr/orgqr/geqp3/trsm_upper, pack/unpack using namespace cbk::test; -/* --------------------------- one test case -------------------------- */ +/* ------------------- one (layout, side, trans) case ------------------ */ +/* A is the nq x k reflector batch, nq = m (side 'L') or n ('R'); C is m x n. */ -template static int run_case(int nm, int m, int nrhs) +template +static int run_case(char layout, char side, char trans, int nm, int m, int n, int k) { - const int k = m; + const bool row = (layout == 'R'); + const bool left = (side == 'L'); + const bool tran = (trans == 'T'); + const int nq = left ? m : n; /* the order of Q, and the rows of A */ + + /* (H, tau) of a random nq x k matrix, a random C, and the two references: + * LAPACKE_?ormqr's op(Q) C, and the explicit Q of ?orgqr applied by ?gemm + * (dqrt03) */ + MatrixBatch H(nm, nq, k), tau(nm, k, 1), C(nm, m, n), Cref(nm, m, n), QC(nm, m, n); + std::vector Qs((size_t)nq * nq); + const auto Q = mat_view(Qs.data(), nq, nq); + for (int idx = 0; idx < nm; ++idx) { + gen_boosted(H.view(idx)); /* diagonal boost tames cond for float */ + ref_geqr2(H.view(idx), tau[idx]); + for (size_t e = 0; e < C.stride(); ++e) + C[idx][e] = frand(); + std::copy(C[idx], C[idx] + C.stride(), Cref[idx]); + ref_ormqr(side, trans, k, H.view(idx), tau[idx], Cref.view(idx)); + + std::fill(Qs.begin(), Qs.end(), T(0)); + copy_matrix(leading(H.view(idx), nq, k), leading(Q, nq, k)); + ref_orgqr(k, Q, tau[idx]); /* Q = H(0)..H(k-1), nq x nq */ + const auto opQ = tran ? Q.transposed() : Q; + if (left) + matmul(opQ, C.view(idx), QC.view(idx)); + else + matmul(C.view(idx), opQ, QC.view(idx)); + } + + /* pack in the requested layout, apply, unpack */ + const int lda = std::max(1, row ? k : nq), ldc = std::max(1, row ? n : m); + std::vector ap = pack_compact(H, lda, V, row); + std::vector tp = pack_tau(tau, V); + std::vector cp = pack_compact(C, ldc, V, row); + int bad = (compact::ormqr(layout, side, trans, m, n, k, ap.data(), lda, tp.data(), + cp.data(), ldc, V, nm) != 0); + MatrixBatch Cout(nm, m, n); + unpack_compact(Cout, cp.data(), ldc, V, row); + + /* checks 1 and 2: vs LAPACKE_?ormqr, and vs the explicit Q (dqrt03) */ + double r_ref = 0, r_q = 0; + for (int idx = 0; idx < nm; ++idx) { + r_ref = std::max(r_ref, test_ratio(diff_norm1(Cout.view(idx), Cref.view(idx)), + nq, norm1(Cref.view(idx)))); + r_q = std::max(r_q, test_ratio(diff_norm1(Cout.view(idx), QC.view(idx)), nq, + norm1(C.view(idx)))); + } + + /* check 3: the inverse application recovers C */ + bad += (compact::ormqr(layout, side, tran ? 'N' : 'T', m, n, k, ap.data(), lda, + tp.data(), cp.data(), ldc, V, nm) != 0); + unpack_compact(Cout, cp.data(), ldc, V, row); + double r_rt = 0; + for (int idx = 0; idx < nm; ++idx) + r_rt = std::max(r_rt, test_ratio(diff_norm1(Cout.view(idx), C.view(idx)), nq, + norm1(C.view(idx)))); + const bool ok = passes(r_ref) && passes(r_q) && passes(r_rt) && (bad == 0); + std::printf( + "T=%-6s V=%-2d nm=%-2d %s side=%c trans=%c m=%-3d n=%-3d k=%-3d | ref %-5.2g " + "%s | Q*C %-5.2g %s | QQt=I %-5.2g %s | %s\n", + compact::name, V, nm, row ? "row" : "col", side, trans, m, n, k, r_ref, + verdict(r_ref), r_q, verdict(r_q), r_rt, verdict(r_rt), ok ? "OK" : "FAIL"); + return !ok; +} + +/* ------------------- the solve pipeline: Q^T B, then R X = Q^T B --------- */ +/* From ?geqr2's reflectors of a square A, B = A X for a known X: the kernel + * applies Q^T, cblas_?trsm back-substitutes, and X must come back (dget04). */ + +template static int run_solve(int nm, int m, int nrhs) +{ + const int k = m; const std::vector Xs = known_solution(m, nrhs); const auto X = mat_view(Xs.data(), m, nrhs); - MatrixBatch A(nm, m, m), Afac(nm, m, m), B(nm, m, nrhs), Bref(nm, m, nrhs), - Bout(nm, m, nrhs), tau(nm, m, 1); + MatrixBatch A(nm, m, m), Afac(nm, m, m), B(nm, m, nrhs), Bout(nm, m, nrhs), + tau(nm, m, 1); for (int kk = 0; kk < nm; ++kk) { gen_boosted(A.view(kk)); /* diagonal boost tames cond for float */ matmul(A.view(kk), X, B.view(kk)); /* B = A X */ - - std::copy(A[kk], A[kk] + (size_t)m * m, Afac[kk]); + std::copy(A[kk], A[kk] + A.stride(), Afac[kk]); ref_geqr2(Afac.view(kk), tau[kk]); - std::copy(B[kk], B[kk] + (size_t)m * nrhs, Bref[kk]); - ref_ormqr('T', k, Afac.view(kk), tau[kk], Bref.view(kk)); } const int ld = std::max(1, m); std::vector ap = pack_compact(Afac, ld, V); std::vector tp = pack_tau(tau, V); std::vector bp = pack_compact(B, ld, V); - - /* check 1: compact Q^T B vs LAPACKE_?ormqr, relative to ||Q^T B||_1 */ - compact::ormqr('T', m, nrhs, k, ap.data(), ld, tp.data(), bp.data(), ld, V, nm); + const int info = compact::ormqr('C', 'L', 'T', m, nrhs, k, ap.data(), ld, + tp.data(), bp.data(), ld, V, nm); unpack_compact(Bout, bp.data(), ld, V); - double r1 = 0; - for (int kk = 0; kk < nm; ++kk) - r1 = std::max(r1, test_ratio(diff_norm1(Bout.view(kk), Bref.view(kk)), m, - norm1(Bref.view(kk)))); - /* check 2: solve recovers X (dget04, discounted by rcond(A)) */ - double r2 = 0; + double r = 0; for (int kk = 0; kk < nm; ++kk) { ref_trsm_upper(Afac.view(kk), Bout.view(kk)); - r2 = std::max(r2, forward_ratio(Bout.view(kk), X, rcond1(A.view(kk)))); + r = std::max(r, forward_ratio(Bout.view(kk), X, rcond1(A.view(kk)))); } - /* check 3: 'N' undoes 'T', relative to ||B||_1 */ - compact::ormqr('N', m, nrhs, k, ap.data(), ld, tp.data(), bp.data(), ld, V, nm); - unpack_compact(Bout, bp.data(), ld, V); - double r3 = 0; - for (int kk = 0; kk < nm; ++kk) - r3 = std::max(r3, test_ratio(diff_norm1(Bout.view(kk), B.view(kk)), m, - norm1(B.view(kk)))); - - const bool ok1 = passes(r1), ok2 = passes(r2), ok3 = passes(r3); - std::printf( - "T=%-6s V=%-2d nm=%-2d m=%-3d nrhs=%d | QtB: %-5.2g %s | solve X: %-5.2g %s " - "| QQt=I: %-5.2g %s\n", - compact::name, V, nm, m, nrhs, r1, verdict(r1), r2, verdict(r2), r3, - verdict(r3)); - return !ok1 + !ok2 + !ok3; + const bool ok = passes(r) && (info == 0); + std::printf("T=%-6s V=%-2d nm=%-2d m=%-3d nrhs=%d | solve X: %-5.2g %s | info=%d\n", + compact::name, V, nm, m, nrhs, r, verdict(r), info); + return !ok; } /* ------------------- pivoted-QR + back-permutation solve ------------ */ @@ -122,7 +181,8 @@ template static int run_case_pivoted(int nm, int m, int nrhs) std::vector bp = pack_compact(B, m, V); /* kernel: c := Q^T b */ - compact::ormqr('T', m, nrhs, k, ap.data(), m, tp.data(), bp.data(), m, V, nm); + compact::ormqr('C', 'L', 'T', m, nrhs, k, ap.data(), m, tp.data(), bp.data(), m, V, + nm); unpack_compact(Bout, bp.data(), m, V); /* R y = c, then back-permute x(jpvt(j)) = y(j); compare against X (dget04, @@ -163,8 +223,8 @@ template static void bench(int V, int nm, int m, int nrhs, int reps) timespec t0, t1; clock_gettime(CLOCK_MONOTONIC, &t0); for (int r = 0; r < reps; ++r) - compact::ormqr((r & 1) ? 'N' : 'T', m, nrhs, k, ap.data(), m, tp.data(), - bp.data(), m, V, nm); + compact::ormqr('C', 'L', (r & 1) ? 'N' : 'T', m, nrhs, k, ap.data(), m, + tp.data(), bp.data(), m, V, nm); clock_gettime(CLOCK_MONOTONIC, &t1); double sec = (t1.tv_sec - t0.tv_sec) + 1e-9 * (t1.tv_nsec - t0.tv_nsec); @@ -183,28 +243,37 @@ template static void bench(int V, int nm, int m, int nrhs, int reps) * (return -j for the j-th argument), not assert or miscompute. */ static int test_validation() { - const int m = 8, nrhs = 2, k = 8, V = 4, nm = 4, ld = 8; - std::vector ap((size_t)ld * k * V, 0), tau((size_t)k * V, 0), - bp((size_t)ld * nrhs * V, 0); - auto call = [&](char tr, int m_, int nrhs_, int k_, int ldap_, int ldbp_, int V_, - int nm_) { - return dormqr_compact(tr, m_, nrhs_, k_, ap.data(), ldap_, tau.data(), bp.data(), - ldbp_, V_, nm_); + const int m = 8, n = 6, k = 5, V = 4, nm = 4; /* A is nq x k, nq = m or n */ + std::vector ap((size_t)m * m * V, 0), tau((size_t)m * V, 0), + cp((size_t)m * m * V, 0); + auto call = [&](char lay, char si, char tr, int m_, int n_, int k_, int ldap_, + int ldcp_, int V_, int nm_) { + return dormqr_compact(lay, si, tr, m_, n_, k_, ap.data(), ldap_, tau.data(), + cp.data(), ldcp_, V_, nm_); }; // clang-format off const ApiCheck t[] = { - {"valid", call('T', m, nrhs, k, ld, ld, V, nm), 0}, - {"bad trans", call('X', m, nrhs, k, ld, ld, V, nm), -1}, - {"m<0", call('T', -1, nrhs, k, ld, ld, V, nm), -2}, - {"nrhs<0", call('T', m, -1, k, ld, ld, V, nm), -3}, - {"k>m", call('T', m, nrhs, m+1, ld, ld, V, nm), -4}, - {"ldapm (L)", call('C', 'L', 'T', m, n, m+1, m, m, V, nm), -6}, + {"k>n (R)", call('C', 'R', 'T', m, n, n+1, n, m, V, nm), -6}, + {"ldap(4, m, nrhs); - fails += run_case(8, m, nrhs); - fails += run_case(16, m, nrhs); - fails += run_case(11, m, nrhs); /* padded partial group */ - fails += - run_case(40, m, nrhs); /* 10 groups: OpenMP path when threads <= 10 */ - fails += run_case(8, m, nrhs); - fails += run_case(16, m, nrhs); - fails += run_case(32, m, nrhs); - - /* the small-dimension cross product (small_dims): m x nrhs, the empty - * operand and the 1/2/3-column RHS blocks included, a padded group */ + /* the full layout x side x trans product (dqrt03): A and Q are m x m + * (k <= m) on the left, n x n (k <= n) on the right; a padded partial + * group and k < nq among the shapes */ + for (char lay : {'C', 'R'}) + for (char side : {'L', 'R'}) + for (char tr : {'N', 'T'}) { + fails += + run_case(lay, side, tr, 8, 43, 5, side == 'L' ? 43 : 5); + fails += + run_case(lay, side, tr, 11, 32, 6, side == 'L' ? 20 : 4); + fails += + run_case(lay, side, tr, 16, 32, 8, side == 'L' ? 32 : 8); + } + /* the other widths, and 10 groups (the OpenMP path when threads <= 10), + * on the tuned route */ + fails += run_case('C', 'L', 'T', 4, m, nrhs, m); + fails += run_case('C', 'L', 'T', 40, m, nrhs, m); + fails += run_case('C', 'L', 'T', 32, m, nrhs, m); + + /* the small-dimension cross product (small_dims): m x n with LAPACK's K + * values {nq, 0, 1, nq/2} on the tuned route, k = nq on the transposed + * (side 'R', row-major) one; the empty operand and the 1/2/3-column + * blocks included, a padded group */ for (int mm : small_dims) - for (int nr : small_dims) { - fails += run_case(5, mm, nr); - fails += run_case(9, mm, nr); + for (int nn : small_dims) { + const int ks[] = {mm, 0, 1, mm / 2}; + for (int t = 0; t < 4; ++t) { + const int kk = ks[t]; + if (kk > mm || std::find(ks, ks + t, kk) != ks + t) continue; + fails += run_case('C', 'L', 'T', 5, mm, nn, kk); + } + fails += run_case('R', 'R', 'N', 9, mm, nn, nn); } + /* the solve pipeline: Q^T B then back substitution recovers X */ + fails += run_solve(8, m, nrhs); + fails += run_solve(11, m, nrhs); /* padded partial group */ + fails += run_solve(16, m, nrhs); + /* column-pivoted QR + back-permutation solve (salvaged cross-check) */ fails += run_case_pivoted(8, m, nrhs); fails += run_case_pivoted(11, m, nrhs); /* padded partial group */ From 4fc740b4903defd0cf1b6c4c84375e0cd2eb07ef Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 15 Sep 2026 00:19:02 +0000 Subject: [PATCH 5/6] changelog: the ?ormqr_compact signature and the test-suite conventions Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_014ZZ76AgYwQ3vtuyZ46tyMX --- CHANGELOG.md | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 651a2f8..cd71771 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -9,6 +9,22 @@ release may change the API. ### Changed +- **Breaking:** the portable `?ormqr_compact` takes `layout`, `side` and + `trans` -- `(layout, side, trans, m, n, k, ap, ldap, taup, cp, ldcp, V, nm)`, + LAPACK `?ormqr`'s argument order behind the layout, as the MKL-style + `cbk_?ormqr_compact` already did -- so Q can be applied from either side in + either layout; it was side `'L'`, column-major only. A caller of the old + signature `(trans, m, nrhs, k, ...)` now writes `('C', 'L', trans, m, nrhs, + k, ...)`. Its argument validation renumbers accordingly (`-1 layout`, + `-2 side`, `-3 trans`, ..., `-13 nm`). +- Tests: the portable suites gate every numerical check on a dimensionless + test ratio in the exact form of the corresponding LAPACK checker (dqrt01, + dqrt03, dpot01, dsyt01, dget02, dtrt02, dget04, ...) against one threshold, + `THRESH = 30`, instead of per-check tolerances; sweep the small-dimension + cross product (`m, n` from 0 to 5, the empty operand included) in every + routine; and check related routines against each other (`potrs` and + `sytrsnp` against the `trsm` sweeps they are built from, `ormqr` in every + layout, side and trans against the explicit Q of `orgqr`). - Tests: every suite validates against a real LAPACKE + CBLAS stack (OpenBLAS, Netlib, or MKL's own with the MKL extension; `cmake/FindLAPACKE.cmake` over `find_package(LAPACK)` and `BLA_VENDOR`) From 27d1e41fa7f5a3950e72fbe0434e617ec0065ff9 Mon Sep 17 00:00:00 2001 From: Claude Date: Tue, 15 Sep 2026 00:37:35 +0000 Subject: [PATCH 6/6] release: 0.2.0 The portable ?ormqr_compact changed signature (layout, side and trans), which under the changelog's pre-1.0 rule is a minor release. CBK_VERSION in cbk.h and project(VERSION) in CMakeLists.txt move to 0.2.0 together (cbk.cpp's static_assert keeps them in step; the shared library becomes libcbk.so.0.2.0, soname unchanged), and the changelog's Unreleased section -- the LAPACKE-backed suites, the test-ratio convention, and the ?ormqr_compact change -- becomes the 0.2.0 entry. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_014ZZ76AgYwQ3vtuyZ46tyMX --- CHANGELOG.md | 3 +++ CMakeLists.txt | 4 ++-- include/cbk.h | 4 ++-- 3 files changed, 7 insertions(+), 4 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index cd71771..f7a3fec 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,8 @@ release may change the API. ## [Unreleased] +## [0.2.0] - 2026-09-15 + ### Changed - **Breaking:** the portable `?ormqr_compact` takes `layout`, `side` and @@ -72,4 +74,5 @@ and the solves built on them, for many small matrices in the compact assumed well conditioned and within range. Deferred to a later release. - Tested on x86-64 Linux with GCC and Clang, against Intel MKL 2020.4 (lp64 and ilp64). +[0.2.0]: https://github.com/ivan-pi/cbk/releases/tag/v0.2.0 [0.1.0]: https://github.com/ivan-pi/cbk/releases/tag/v0.1.0 diff --git a/CMakeLists.txt b/CMakeLists.txt index 1d1a5ef..0bb53c5 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -3,7 +3,7 @@ cmake_minimum_required(VERSION 3.28) project(cbk - VERSION 0.1.0 + VERSION 0.2.0 DESCRIPTION "Compact Batch Kernels: batched QR, Cholesky and LDL^T for small matrices in the compact (interleaved) format" LANGUAGES CXX) @@ -73,7 +73,7 @@ target_compile_definitions(cbk PRIVATE if(CBK_WITH_OPENMP AND OpenMP_CXX_FOUND) target_link_libraries(cbk PRIVATE OpenMP::OpenMP_CXX) endif() -# libcbk.so.<0.1.0> with soname libcbk.so.<0>, from project(... VERSION) above. +# libcbk.so.<0.2.0> with soname libcbk.so.<0>, from project(... VERSION) above. set_target_properties(cbk PROPERTIES VERSION ${PROJECT_VERSION} SOVERSION ${PROJECT_VERSION_MAJOR}) add_library(cbk::cbk ALIAS cbk) diff --git a/include/cbk.h b/include/cbk.h index 0041d80..75432a0 100644 --- a/include/cbk.h +++ b/include/cbk.h @@ -58,12 +58,12 @@ /* Version. The macros describe the header; cbk_get_version() the library * linked, as CBK_VERSION_MAJOR * 10000 + CBK_VERSION_MINOR * 100 + - * CBK_VERSION_PATCH (0.1.0 -> 100), so a consumer can check the two agree. + * CBK_VERSION_PATCH (0.2.0 -> 200), so a consumer can check the two agree. * Macros, not an enum: this is a C header, and a consumer must be able to * test the version in the preprocessor (#if CBK_VERSION >= ...). */ /* NOLINTBEGIN(modernize-macro-to-enum) */ #define CBK_VERSION_MAJOR 0 -#define CBK_VERSION_MINOR 1 +#define CBK_VERSION_MINOR 2 #define CBK_VERSION_PATCH 0 #define CBK_VERSION \ (CBK_VERSION_MAJOR * 10000 + CBK_VERSION_MINOR * 100 + CBK_VERSION_PATCH)