LAPACK-style test suites (ratios, dimension sweeps, cross-routine checks) and a side/layout-general portable ?ormqr_compact (0.2.0)

.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<T>` 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

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

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)

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

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

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 <typename T>
-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<T, decltype(v)::value>(true, false, trans, m, nrhs, k,
-                                                          ap, ldap, taup, bp, ldbp, nm);
+        cbk::detail::ormqr_compact<T, decltype(v)::value>(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)                                   \