diff --git a/.claude/CLAUDE.md b/.claude/CLAUDE.md index b62564f..46a2f71 100644 --- a/.claude/CLAUDE.md +++ b/.claude/CLAUDE.md @@ -251,6 +251,30 @@ 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 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`, + `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 API (`cbk.h`) validates LAPACK-style, returning `-j` for a bad j-th diff --git a/CHANGELOG.md b/CHANGELOG.md index 651a2f8..f7a3fec 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,8 +7,26 @@ release may change the API. ## [Unreleased] +## [0.2.0] - 2026-09-15 + ### 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`) @@ -56,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/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..75432a0 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 @@ -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) @@ -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 4e30c39..d935b11 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,90 @@ 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. +// +// 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) +{ + return err / (std::max(1, dim) * std::max(norm, norm_floor) * + (double)std::numeric_limits::epsilon()); +} + +template double residual_ratio(double anorm, Rv R, Xv X) +{ + 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) +{ + 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) { @@ -104,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; @@ -157,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); } \ @@ -343,12 +460,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 +561,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..f722ad8 100644 --- a/tests/test_gels_compact.cpp +++ b/tests/test_gels_compact.cpp @@ -8,17 +8,20 @@ // 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 // 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 @@ -47,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 @@ -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 */ @@ -77,20 +79,26 @@ 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(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: 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: 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); const auto R = mat_view(Rs.data(), rows_op, nrhs); const auto S = mat_view(Ss.data(), cols_op, nrhs); @@ -98,8 +106,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 @@ -107,12 +114,13 @@ 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) - e_prop = std::max(e_prop, std::abs((double)Ss[e]) / 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) - 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); @@ -120,55 +128,51 @@ 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(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 + // 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(Xmn, 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; } @@ -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..3eb9b81 100644 --- a/tests/test_geqrf_compact.cpp +++ b/tests/test_geqrf_compact.cpp @@ -6,18 +6,20 @@ // 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) -// plus LAPACK-style argument validation of the C API. +// 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 +// API. // // Assisted-by: Claude:claude-opus-4-8 Claude:claude-fable-5 @@ -38,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); @@ -49,37 +50,42 @@ 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 - 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(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)); - // 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(diff_norm1(Rec, A.view(idx)), 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); @@ -87,29 +93,27 @@ 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('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(), n, V); - e_solve = 0; + unpack_compact(Bo, bp.data(), lda, V); + 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(Bo.view(idx), 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; } @@ -149,7 +153,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); @@ -163,7 +166,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), @@ -176,20 +179,19 @@ 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||_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))); - worst_res = std::max(worst_res, resid / std::max(norm1(Am), norm_floor)); + 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 + // 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 @@ -200,23 +202,22 @@ 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)); - worst_el = std::max(worst_el, el / std::max(norm1(Am), norm_floor)); + 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))); } } - 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; } @@ -255,7 +256,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) { @@ -272,7 +272,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 @@ -293,14 +293,13 @@ 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(diff_norm1(Rec, A.view(idx)), 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; } @@ -355,6 +354,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..219799f 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); } \ @@ -115,12 +121,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,28 +215,50 @@ 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); } -// ?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)); - 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 @@ -235,6 +266,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 +289,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 +339,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); @@ -372,10 +413,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 / 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; +}; + +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()); + 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) { + 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; +} + +// 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 db40bf6..780ef8f 100644 --- a/tests/test_orgqr_compact.cpp +++ b/tests/test_orgqr_compact.cpp @@ -8,17 +8,20 @@ * 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), - * 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 */ @@ -41,12 +44,8 @@ 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 = 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 */ @@ -78,18 +77,20 @@ 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(diff_norm1(Qout.view(kk), Qref.view(kk)), 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||_1 / (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_norm1(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); @@ -99,8 +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); - 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(diff_norm1(Rec, A0.view(kk)), m, norm1(A0.view(kk)))); } } @@ -118,13 +119,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; } @@ -192,6 +192,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..4aaf31f 100644 --- a/tests/test_ormqr_compact.cpp +++ b/tests/test_ormqr_compact.cpp @@ -1,19 +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): - * 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 - * 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 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 */ @@ -27,68 +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 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 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)); } - 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); + 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); - /* 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); - 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) / - std::max(norm1(Bref.view(kk)), norm_floor)); - - /* check 2: solve recovers X */ - double e2 = 0; + double r = 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)); + r = std::max(r, forward_ratio(Bout.view(kk), X, rcond1(A.view(kk)))); } - /* 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); - 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)); - - /* 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"); - 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 ------------ */ @@ -101,21 +161,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]); } @@ -124,11 +181,13 @@ 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 */ - 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) { @@ -137,12 +196,12 @@ 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(x, 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; } @@ -164,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); @@ -184,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 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 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); diff --git a/tests/test_potrf_compact.cpp b/tests/test_potrf_compact.cpp index d505fa7..2327e3f 100644 --- a/tests/test_potrf_compact.cpp +++ b/tests/test_potrf_compact.cpp @@ -7,18 +7,25 @@ // (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 // ?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. +// 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 +// validation of the three C APIs. // // Assisted-by: Claude:claude-opus-4-8 Claude:claude-fable-5 @@ -37,7 +44,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'); @@ -51,34 +57,40 @@ 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; + 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) { 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))); } - e_fac = std::max(e_fac, el / std::max(norm1(Ref), norm_floor)); + 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 // ?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); @@ -86,23 +98,17 @@ 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(diff_norm1(Rec, Ain), 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); } @@ -114,9 +120,8 @@ 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 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,32 +132,85 @@ 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); - 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; +} + +// ------------- 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; } @@ -222,12 +280,13 @@ 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) { 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; @@ -237,20 +296,21 @@ 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 - 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); } @@ -306,6 +366,25 @@ 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) + 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..3d76c55 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 @@ -20,7 +21,11 @@ // ?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: +// 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: // - 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 @@ -75,12 +80,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); @@ -94,14 +98,11 @@ 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(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", + compact::name, uplo, n, r, verdict(r)); return !ok; } @@ -109,7 +110,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'); @@ -122,47 +122,51 @@ 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; + 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 (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 - for (int j = 0; j < n; ++j) + double e_fac = 0, e_rec = 0; + 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)); } - - // check 2: reconstruction of A from the named triangle's (D, L|U) - 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))); + 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)); + 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); } @@ -174,9 +178,8 @@ 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 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,34 +190,92 @@ 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); - 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, + 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; +} + +// ------------- 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; } @@ -290,23 +351,27 @@ 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); - for (int j = 0; j < n; ++j) + double e = 0; + 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 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; } @@ -349,12 +414,14 @@ 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); - for (int j = 0; j < n; ++j) + double e_sib = 0; + 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,20 +431,24 @@ 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)))); } - 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); } @@ -436,6 +507,25 @@ 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) + 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..37916ea 100644 --- a/tests/test_trsm_compact.cpp +++ b/tests/test_trsm_compact.cpp @@ -9,10 +9,15 @@ // 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 -// 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). +// interleave widths (column-major, the tuned path), as test ratios +// 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). // // Assisted-by: Claude:claude-opus-4-8 Claude:claude-fable-5 @@ -41,7 +46,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, @@ -51,54 +55,69 @@ 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)); } // 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); - - // 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); - const auto R = mat_view(Rs.data(), m, n), aB = mat_view(aBs.data(), m, n); + unpack_compact(Bout, bp.data(), ldb, V); + + // 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), 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) { - 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) + 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]; - worst_res = - std::max(worst_res, max_abs_diff(Rs.data(), aBs.data(), (size_t)m * n) / - std::max(norm1(aB), norm_floor)); + 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 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_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 %.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 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_ref, verdict(r_ref), r_res, verdict(r_res), info, + (ok && info == 0) ? "OK" : "FAIL"); return (info != 0) + !ok; } @@ -177,5 +196,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); }