diff --git a/CMakeLists.txt b/CMakeLists.txt index 3b6c4ca..eba64db 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -33,7 +33,8 @@ add_library(cqr_compact src/cqr_ormqr_compact_dispatch.cpp src/cqr_geqrf_compact_dispatch.cpp src/cqr_potrf_compact_dispatch.cpp - src/cqr_trsm_compact_dispatch.cpp) + src/cqr_trsm_compact_dispatch.cpp + src/cqr_syrk_compact_dispatch.cpp) target_include_directories(cqr_compact PUBLIC "${CMAKE_CURRENT_SOURCE_DIR}/src") target_compile_features(cqr_compact PUBLIC cxx_std_17) add_library(cqr::compact ALIAS cqr_compact) @@ -50,7 +51,8 @@ if(CQR_WITH_MKL) src/cqr_mkl_ormqr.cpp src/cqr_mkl_geqrf.cpp src/cqr_mkl_potrf.cpp - src/cqr_mkl_trsm.cpp) + src/cqr_mkl_trsm.cpp + src/cqr_mkl_syrk.cpp) target_include_directories(cqr_mkl_ext PUBLIC "${CMAKE_CURRENT_SOURCE_DIR}/src") target_compile_features(cqr_mkl_ext PUBLIC cxx_std_17) @@ -132,6 +134,13 @@ if(CQR_BUILD_TESTS) target_link_libraries(test_cqr_trsm_compact PRIVATE cqr_compact) add_test(NAME portable_trsm COMMAND test_cqr_trsm_compact) + # Portable self-contained symmetric rank-k update test (no BLAS): C API + # argument validation + the compact syrk kernel vs a scalar ?syrk reference, + # over the full uplo x trans x layout matrix. + add_executable(test_cqr_syrk_compact src/test_cqr_syrk_compact.cpp) + target_link_libraries(test_cqr_syrk_compact PRIVATE cqr_compact) + add_test(NAME portable_syrk COMMAND test_cqr_syrk_compact) + # MKL-backed validation of cqr_mkl_dormqr_compact (design doc section 7). if(CQR_WITH_MKL) add_executable(test_cqr_ormqr_mkl src/test_cqr_ormqr_mkl.cpp) @@ -165,6 +174,15 @@ if(CQR_BUILD_TESTS) PRIVATE cqr_mkl_ext cqr_compact MKL::Compact) add_test(NAME mkl_trsm_suites COMMAND test_cqr_trsm_mkl) + # MKL cross-check of cqr_mkl_?syrk_compact: vs per-matrix cblas_?syrk and vs + # mkl_?gemm_compact over the full feature matrix, plus an end-to-end + # Cholesky-QR (syrk Gram matrix -> potrf -> trsm) driven entirely by cqr + # kernels (no MKL compute). + add_executable(test_cqr_syrk_mkl src/test_cqr_syrk_mkl.cpp) + target_link_libraries(test_cqr_syrk_mkl + PRIVATE cqr_mkl_ext cqr_compact MKL::Compact) + add_test(NAME mkl_syrk_suites COMMAND test_cqr_syrk_mkl) + # The compact QR-solve example self-validates, so register it too. add_test(NAME example_solve_qr_compact COMMAND solve_qr_compact) diff --git a/PLANS.md b/PLANS.md index f391c95..384a78c 100644 --- a/PLANS.md +++ b/PLANS.md @@ -115,6 +115,42 @@ columns). Remaining performance-only opportunities: the small-`n` (`~10`) per-group overhead (~`0.6-0.9x`), reciprocal-multiplying the diagonal in the blocked paths, and SIMD-tuning the strided kernel. +## syrk (`cqr_mkl_dsyrk_compact`) + +The compact batched symmetric rank-k update (`cqr_mkl_dsyrk_compact_design.md`): a +portable, vectorized `mkl_?syrk_compact` (which MKL omits), completing the compact +BLAS-3 set and forming the Gram matrix of a Cholesky QR. Status vs. its design +document: + +- **Implemented (design sections 2-6, 8.1):** both API surfaces -- the MKL-style + `cqr_mkl_?syrk_compact` (drop-in style, no `work`/`info`, dispatching on the + `MKL_COMPACT_PACK` format enum) and the portable `dsyrk_compact`/`ssyrk_compact` + (LAPACK/BLAS-style `info = -j` validation) -- over the vectorized dot-product + rank-k update. Column-major `trans='T'` (`C = A^T A`, the Cholesky-QR Gram + matrix in LAPACK's native layout) is the tuned path: a `JB = 4` register-blocked + contiguous column-dot; the other three trans/layout combinations route through a + stride-generalized kernel. Full `uplo x trans` in FP64/FP32, `beta = 0` handled + as the BLAS overwrite (C not read), `k = 0` / `alpha = 0` reducing to + `C := beta C`. +- **Validated (design section 7):** a BLAS-free portable test vs a scalar `?syrk` + reference over the full `uplo x trans x layout` matrix (the reference accumulates + in a different order, so a bug shared by reference and kernel cannot pass), plus + an MKL test cross-checking vs per-matrix `cblas_?syrk` (whole matrix: active + triangle correct + opposite untouched) and vs `mkl_?gemm_compact` (active + triangle), and closing an end-to-end Cholesky QR (`cqr_mkl_dsyrk_compact -> + cqr_mkl_dpotrf_compact -> cqr_mkl_dtrsm_compact`) that recovers `A = Q R` with + `Q^T Q = I`. All are CTest-registered. +- **Known gaps / scoped out (design section 6.6):** the strided (`trans='N'` and + both row-major) inner sweep is correctness-first, not separately SIMD-tuned -- + `trans='N'` row-major is itself a contiguous case and a natural future second + tuned path, mirroring `potrf`'s row-major-upper dual. No overflow/underflow-safe + scaling. The complex (`c`/`z`) symmetric update is `?herk`, out of scope. The + Cholesky-QR orthogonality caveat (`cond(A)^2 * eps`) is a property of that + algorithm, not of `?syrk`. +- **Deferred:** a throughput benchmark (`cqr_mkl_?syrk_compact` vs the general + `mkl_?gemm_compact`), and a worked Cholesky-QR solve example/benchmark mirroring + `solve_qr_compact`/`bench_qr_compact`, are left for a future change. + ## Known gaps Gaps between the `ormqr` implementation and its design document diff --git a/README.md b/README.md index a4e73e5..3ec4b52 100644 --- a/README.md +++ b/README.md @@ -21,6 +21,12 @@ an Intel MKL-style API: `mkl_?trsm_compact` (the batched triangular solve), so the whole `AX = B` pipeline runs with no MKL compute kernel. See its [design document](cqr_mkl_dtrsm_compact_design.md). +* **`cqr_mkl_?syrk_compact`** - the batched **symmetric rank-k update** + (`C := alpha A A^T + beta C`), the Compact-format `?syrk` MKL omits. It exploits + the symmetry MKL's compact `gemm` cannot (half the flops, a triangle of writes) + and forms the Gram matrix `A^T A` of a **Cholesky QR**, feeding + `cqr_mkl_?potrf_compact` then `cqr_mkl_?trsm_compact` to factor `A = Q R` with no + MKL compute kernel. See its [design document](cqr_mkl_dsyrk_compact_design.md). All routines come in single and double precision. Together they factor and solve batched systems entirely in the compact format, with no MKL compute kernel. @@ -99,8 +105,8 @@ include: | File | Role | |------|------| -| `src/cqr_mkl_ext.h` | The MKL-style public API: `cqr_mkl_?geqrf_compact` (QR factorization, drop-in for `mkl_?geqrf_compact`), `cqr_mkl_?ormqr_compact` (apply `Q`/`Q^T`, the missing `mkl_?ormqr_compact`), `cqr_mkl_?potrf_compact` (Cholesky, drop-in for `mkl_?potrf_compact`), and `cqr_mkl_?trsm_compact` (triangular solve, drop-in for `mkl_?trsm_compact`). Takes `MKL_COMPACT_PACK` formats. | -| `src/cqr_compact.h` | The portable C API, all eight exported functions: `dgeqrf_compact` / `sgeqrf_compact` (QR factorization), `dormqr_compact` / `sormqr_compact` (apply `Q` / `Q^T`), `dpotrf_compact` / `spotrf_compact` (Cholesky), and `dtrsm_compact` / `strsm_compact` (triangular solve), with an explicit interleave width `V` and no MKL dependency. | +| `src/cqr_mkl_ext.h` | The MKL-style public API: `cqr_mkl_?geqrf_compact` (QR factorization, drop-in for `mkl_?geqrf_compact`), `cqr_mkl_?ormqr_compact` (apply `Q`/`Q^T`, the missing `mkl_?ormqr_compact`), `cqr_mkl_?potrf_compact` (Cholesky, drop-in for `mkl_?potrf_compact`), `cqr_mkl_?trsm_compact` (triangular solve, drop-in for `mkl_?trsm_compact`), and `cqr_mkl_?syrk_compact` (symmetric rank-k update, the missing `mkl_?syrk_compact`). Takes `MKL_COMPACT_PACK` formats. | +| `src/cqr_compact.h` | The portable C API, all ten exported functions: `dgeqrf_compact` / `sgeqrf_compact` (QR factorization), `dormqr_compact` / `sormqr_compact` (apply `Q` / `Q^T`), `dpotrf_compact` / `spotrf_compact` (Cholesky), `dtrsm_compact` / `strsm_compact` (triangular solve), and `dsyrk_compact` / `ssyrk_compact` (symmetric rank-k update), with an explicit interleave width `V` and no MKL dependency. | Everything else under `src/` is internal - implementation details and tests, not part of the supported interface: @@ -112,14 +118,17 @@ not part of the supported interface: | `src/cqr_potrf_compact.hpp` | Templated SIMD Cholesky-factorization kernel (vectorized `potf2`; scalar `T`, interleave width `V`). | | `src/cqr_ormqr_compact.hpp` | Templated SIMD kernel `B := op(Q)*B` (scalar `T`, interleave width `V`). | | `src/cqr_trsm_compact.hpp` | Templated compact triangular-solve kernels (tuned column-major/left + general strided) and group driver (scalar `T`, interleave width `V`). | +| `src/cqr_syrk_compact.hpp` | Templated compact symmetric rank-k update kernels (tuned column-major `A^T A` + general strided) and group driver (scalar `T`, interleave width `V`). | | `src/cqr_geqrf_compact_dispatch.cpp` | Portable geqrf C entry points (runtime `V` -> compile-time dispatch). | | `src/cqr_potrf_compact_dispatch.cpp` | Portable potrf C entry points (runtime `V` -> compile-time dispatch). | | `src/cqr_ormqr_compact_dispatch.cpp` | Portable ormqr C entry points (runtime `V` -> compile-time dispatch). | | `src/cqr_trsm_compact_dispatch.cpp` | Portable trsm C entry points with LAPACK/BLAS-style `info = -j` validation (runtime `V` -> compile-time dispatch). | +| `src/cqr_syrk_compact_dispatch.cpp` | Portable syrk C entry points with LAPACK/BLAS-style `info = -j` validation (runtime `V` -> compile-time dispatch). | | `src/cqr_mkl_geqrf.cpp` | Dispatches on `MKL_COMPACT_PACK` directly and maps `MKL_LAYOUT`, then calls the geqrf kernel. | | `src/cqr_mkl_potrf.cpp` | Dispatches on `MKL_COMPACT_PACK` directly and maps `MKL_UPLO`/`MKL_LAYOUT`, then calls the potrf kernel. | | `src/cqr_mkl_ormqr.cpp` | Dispatches on `MKL_COMPACT_PACK` directly and maps `side`/`trans`/`MKL_LAYOUT`, then calls the ormqr kernel. | | `src/cqr_mkl_trsm.cpp` | Dispatches on `MKL_COMPACT_PACK` directly and maps the MKL enums (`MKL_SIDE`/`MKL_UPLO`/`MKL_TRANSPOSE`/`MKL_DIAG`/`MKL_LAYOUT`), then calls the trsm kernel (drop-in for `mkl_?trsm_compact`; no `work`/`info`). | +| `src/cqr_mkl_syrk.cpp` | Dispatches on `MKL_COMPACT_PACK` directly and maps `MKL_UPLO`/`MKL_TRANSPOSE`/`MKL_LAYOUT`, then calls the syrk kernel (the missing `mkl_?syrk_compact`; no `work`/`info`). | | `src/cqr_mkl_alloc.h` | Optional RAII buffer helpers (`mkl_alloc_bytes`, `mkl_buffer`) wrapping `mkl_malloc`/`mkl_free`. | | `src/test_compact_util.hpp` | Shared test helpers (seeded RNG, error metrics, SPD generation, Compact pack/unpack); header-only, no MKL. | | `src/test_cqr_geqrf_compact.cpp` | Self-contained geqrf correctness test vs a scalar `geqr2` reference (no BLAS). | @@ -130,6 +139,8 @@ not part of the supported interface: | `src/test_cqr_ormqr_mkl.cpp` | MKL-backed validation through the real compact pipeline. | | `src/test_cqr_trsm_compact.cpp` | Self-contained trsm test (no BLAS): C API validation + numerical vs a scalar `?trsm` reference. | | `src/test_cqr_trsm_mkl.cpp` | MKL-backed cross-check of `cqr_mkl_?trsm_compact` vs `mkl_?trsm_compact` + an end-to-end MKL-compute-free solve. | +| `src/test_cqr_syrk_compact.cpp` | Self-contained syrk test (no BLAS): C API validation + numerical vs a scalar `?syrk` reference. | +| `src/test_cqr_syrk_mkl.cpp` | MKL-backed validation of `cqr_mkl_?syrk_compact` vs `cblas_?syrk` and `mkl_?gemm_compact` + an end-to-end MKL-compute-free Cholesky QR. | ### Examples diff --git a/cqr_mkl_dsyrk_compact_design.md b/cqr_mkl_dsyrk_compact_design.md new file mode 100644 index 0000000..ac0e31e --- /dev/null +++ b/cqr_mkl_dsyrk_compact_design.md @@ -0,0 +1,354 @@ +# API Design Document: `cqr_mkl_dsyrk_compact` + +> Assisted-by: Claude:claude-opus-4.8 + +## 1. Overview + +This extension provides a batched **symmetric rank-k update** for a set of +matrices stored in Intel MKL's Compact (interleaved-batch) format. It is the +Compact-format counterpart of BLAS `?syrk`, completing the compact BLAS-3 set +alongside MKL's own `mkl_?gemm_compact` and this project's +`cqr_mkl_?trsm_compact`. Unlike forming the product through a general +`mkl_?gemm_compact`, it exploits the symmetry of the result -- only one triangle +of the symmetric `n x n` output is referenced -- for roughly half the flops and +half the writes. It mirrors an `mkl_?syrk_compact` in signature and semantics, +but MKL ships no such routine, so like `cqr_mkl_?ormqr_compact` this fills a gap +in the compact ecosystem, as a fully portable, open implementation built on GNU +vector types. + +The symmetric rank-k update is the defining kernel of **Cholesky QR** -- the +factorization this project is named for. For a tall matrix `A` (`m x n`, +`m >= n`), Cholesky QR forms the Gram matrix, factors it, and recovers `Q`: + +``` +cqr_mkl_dsyrk_compact('U','T', A -> G); // G = A^T A (rank-k update) +cqr_mkl_dpotrf_compact('U', G -> R); // G = R^T R (Cholesky) +cqr_mkl_dtrsm_compact ('R','U','N', R, A->Q); // Q = A R^{-1} (triangular solve) +``` + +so `A = Q R` with `Q^T Q = I`. With `cqr_mkl_?syrk_compact` in place, this entire +pipeline runs in the compact format with no MKL compute kernel, exactly as the +QR-based `geqrf -> ormqr -> trsm` solver already does. MKL's compact pack/unpack +helpers (`mkl_?gepack_compact`, `mkl_get_format_compact`, ...) are still used to +move data in and out of the interleaved layout; only the *computation* is cqr's. + +The primary target is many small-to-medium matrices (both dimensions in `3..500`, +emphasis below 170), matching the rest of the toolkit; the tuned path is +column-major, `trans = 'T'` (`C = A^T A`) -- the Gram matrix above, in the +column-major layout LAPACK and MKL produce their factors in. + +## 2. Syntax + +```c +void cqr_mkl_dsyrk_compact ( + MKL_LAYOUT layout, MKL_UPLO uplo, MKL_TRANSPOSE trans, + MKL_INT n, MKL_INT k, double alpha, + const double * ap, MKL_INT ldap, + double beta, double * cp, MKL_INT ldcp, + MKL_COMPACT_PACK format, MKL_INT nm +); +``` + +The signature is exactly BLAS +[`?syrk`](https://www.intel.com/content/www/us/en/docs/onemkl/developer-reference-c/2025-2/syrk.html) +plus the three arguments MKL's compact routines add (`layout`, `format`, `nm`) -- +and, like the other compact BLAS-3 routines (`mkl_?gemm_compact`, +`mkl_?trsm_compact`), no `work`/`lwork`/`info` (there is no workspace and, per the +compact convention, no argument checking). Both real precisions are provided: +`cqr_mkl_dsyrk_compact` (double) and `cqr_mkl_ssyrk_compact` (single). + +## 3. Description + +For each matrix in the batch the routine forms, in place, one of the symmetric +rank-k updates + +``` +C := alpha * A * A^T + beta * C (trans = MKL_NOTRANS, A is n x k) +C := alpha * A^T * A + beta * C (trans = MKL_TRANS, A is k x n) +``` + +where `alpha` and `beta` are scalars, `A` is a general rectangular matrix, and +`C` is the symmetric `n x n` result. Only the triangle of `C` named by `uplo` is +referenced and updated; the opposite triangle is neither read nor written, so on +exit it holds exactly its entry value (the caller may leave it uninitialised). + +`op(A) = A^H` (`MKL_CONJTRANS`) is accepted and, for the real types handled here, +is identical to `A^T`; it is folded to the transpose path. (The genuinely +conjugated update `C := alpha A A^H + beta C` is `?herk`, a distinct routine, and +is out of scope here.) When `beta = 0`, `C` is not read on entry -- its prior +contents, even `NaN`/`Inf`, are overwritten -- exactly as reference BLAS `?syrk` +defines it. `k = 0` and `alpha = 0` are *not* special-cased: both reduce the +update to `C := beta * C` (an empty contraction contributes the zero rank-k +term), which the kernel computes directly. + +**Constraint note.** As with all Compact routines, every matrix in the call +shares the same dimensions (`n`, `k`), leading dimensions (`ldap`, `ldcp`), +storage `layout`, `format`, and the same `uplo`/`trans`. The batch is processed +one *pack* (group of `V` interleaved matrices) at a time; `V` is derived from +`format`. + +**Relationship to MKL and ArmPL.** Neither production interleave-batch library +provides a symmetric rank-k update; both offer only the general batched matrix +multiply, so a batched `syrk` today must be formed through that general multiply, +forgoing the symmetry. `cqr_mkl_?syrk_compact` supplies the missing +specialization, expressed through MKL's compact interface: + +* **Intel MKL** exposes the rank-k update only as the general + `mkl_?gemm_compact` (there is no `mkl_?syrk_compact`); `cqr_mkl_?syrk_compact` + adds the symmetry-exploiting specialization over the same Compact + (interleaved) format, selected by an opaque `MKL_COMPACT_PACK`, with the BLAS + `?syrk` parameter set (`uplo`, `trans`, `n`, `k`, `alpha`, `beta`). +* **Arm Performance Libraries** likewise offer the general batched multiply + (`armpl_?gemm_interleave_batch`) but no interleave-batch `syrk`, so the + symmetric update is unavailable there too and would have to be formed with the + general `gemm`. Where ArmPL's interleave-batch API exposes the layout + explicitly through `ninter`/`nbatch` and arbitrary strides, this routine keeps + MKL's abstraction instead, hiding the interleave width and strides behind + `MKL_COMPACT_PACK` + `MKL_LAYOUT` + the compact leading dimensions, for symmetry + with MKL's native compact ecosystem. ArmPL's real `gemm` interface has no + conjugate transpose; neither do the real types here. + +## 4. Input Parameters + +* **`layout`** (`MKL_LAYOUT`): the in-memory storage order of every matrix in the + batch, `MKL_COL_MAJOR` (tuned path) or `MKL_ROW_MAJOR`. +* **`uplo`** (`MKL_UPLO`): `MKL_UPPER` -- reference and update the upper triangle + of `C`; `MKL_LOWER` -- the lower triangle. The opposite triangle is untouched. +* **`trans`** (`MKL_TRANSPOSE`): `MKL_NOTRANS` -> `C := alpha A A^T + beta C` + (`A` is `n x k`); `MKL_TRANS` -> `C := alpha A^T A + beta C` (`A` is `k x n`); + `MKL_CONJTRANS` -> `A^H` (== `A^T` for the real types, folded to the transpose + path). +* **`n`** (`MKL_INT`): the order of the symmetric matrix `C`, and the + non-contracted dimension of `A` (`n >= 0`). +* **`k`** (`MKL_INT`): the contracted dimension -- the number of columns of `A` + for `trans = MKL_NOTRANS`, or rows for `MKL_TRANS` (`k >= 0`). +* **`alpha`** (`double`): scalar multiplying the rank-k product. +* **`ap`** (`const double *`): the compact buffer of `nm` matrices `A` + (`n x k` for `MKL_NOTRANS`, `k x n` for `MKL_TRANS`), packed with + `mkl_?gepack_compact`. +* **`ldap`** (`MKL_INT`): leading dimension of each `A` within the compact buffer, + `>= max(1, rows(A))` for column-major or `>= max(1, cols(A))` for row-major. +* **`beta`** (`double`): scalar multiplying `C`. When `beta = 0`, `C` is not read. +* **`cp`** (`double *`): the compact buffer of `nm` symmetric matrices `C` + (`n x n`); the `uplo` triangle is overwritten in place. +* **`ldcp`** (`MKL_INT`): leading dimension of each `C` within the compact buffer, + `>= max(1, n)`. +* **`format`** (`MKL_COMPACT_PACK`): the pack format from + `mkl_get_format_compact()`; selects the interleave width `V` + (SSE/AVX/AVX-512 -> 2/4/8 for FP64, 4/8/16 for FP32). +* **`nm`** (`MKL_INT`): total number of matrices in the batch. + +**Buffer alignment.** Any base alignment of `ap`/`cp` is correct. For full speed, +align each buffer's base to the pack width (64 B covers every format) so each +SIMD access stays on one cache line -- worth up to ~40% on small, cache-resident +sizes. `mkl_malloc(bytes, 64)` (the default of this project's `mkl_alloc_bytes`) +already does this. + +## 5. Output Parameters + +* **`cp`**: the `uplo` triangle of each `C` is overwritten with + `alpha op(A) op(A)^T + beta C`, in Compact format; the opposite triangle is + left untouched. + +An unrecognized `format` selects no kernel, so the call is a silent no-op (as for +`cqr_mkl_?trsm_compact`, `?syrk` has no `info` to report a dispatch-level status). + +## 6. Design Considerations & Compatibility + +### 6.1 The algorithm: vectorized dot-product rank-k update + +Each update is formed by the dot-product (inner-product) variant of `?syrk`, +executed `V` matrices at a time. Because Compact format interleaves the `V` +matrices so that element `(i,j)` of all `V` is contiguous, the scalar algorithm +lifts verbatim with `double -> V`-wide vector: every `*`, `+`, and the final +`alpha * . + beta * .` becomes a lane-wise SIMD operation over `V` independent +matrices. There is no data-dependent branch to mask, so the lift is direct. + +Each output `C(i,j)` is a length-`k` inner product of two "vectors" of `A` along +the contraction axis: + +``` +C(i,j) = alpha * sum_{p=0..k-1} A_i(p) * A_j(p) [ + beta * C(i,j) ] +``` + +where `A_i` is row `i` of `A` for `trans = 'N'` (contracting over columns) or +column `i` of `A` for `trans = 'T'` (contracting over rows). The dot-product form +is preferred over the rank-1 (outer-product) form for three reasons: `C` is +written exactly once (minimal `C` traffic, and the `beta` fold happens on that +single store); the arithmetic is branch-free; and it **blocks cleanly against the +triangle** -- for a fixed row `i`, every column `j` in that row's triangle is a +*full* length-`k` dot, with no diagonal-corner peeling. + +### 6.2 Register blocking (`JB = 4`) + +The output columns are register-blocked four at a time (`JB = 4`): for a fixed +row `i`, one load of the `A`-vector `A_i(p)` is reused across the four columns +`A_{j..j+3}(p)`, halving the dominant `A` traffic -- the same reuse trick the +`ormqr` kernel applies across its RHS columns and `trsm` across its RHS block. The +one-, two-, and three-column tails of the triangle are handled by a scalar +remainder loop (the triangle length varies per row, so unlike `trsm` there is no +single fixed tail to specialize; the remainder is at most three columns per row). + +### 6.3 Layouts and transpose: one tuned contiguous kernel, one strided + +There are four `(trans, layout)` combinations. In the dot-product form the +contraction axis of `A` should be unit-stride for the inner loop to stream +contiguously, and exactly one combination achieves that in the natural layout: + +* **Contiguous (tuned).** Column-major + `trans = 'T'` (`C = A^T A`): the + contracted vectors are the *columns* of `A` (`A` is `k x n`), which are + contiguous in the column-major compact buffer, so the length-`k` dot streams + both operands down contiguous packs. This is the Cholesky-QR Gram matrix in the + layout LAPACK/MKL produce factors in, so it is the case the toolkit most wants + fast -- mirroring how `trsm` tunes the `side='L'` back-substitution and `potrf` + the column-major lower factor. +* **Strided.** Column-major + `trans = 'N'` (the contracted vectors are *rows* of + `A`, strided by `ldap`) and both row-major cases are supported for MKL + compatibility through a stride-generalized kernel over the same dot-product math + (correctness-first; the non-contiguous inner sweep is not separately + SIMD-tuned), reusing the existing `BatchView` addressing. `trans = 'N'` + row-major is itself a contiguous case (rows of `A` are contiguous row-major) and + is a natural candidate for a future second tuned path, exactly as `potrf`'s + row-major-upper dual is; it is left on the strided path here to keep the tuned + surface small, as `trsm` does. + +The `uplo` triangle is defined on the mathematical indices `(i,j)` of the +symmetric `C` regardless of `trans` or `layout`; only the strides of the two +`BatchView`s differ between the strided cases. + +### 6.4 Padding and SIMD semantics + +When `nm` is not a multiple of `V`, `mkl_?gepack_compact` fills the unused slots +of the last pack with identity matrices. The rank-k update needs no special +handling for them: the kernel runs the whole final pack unmasked at full width, +computing some (garbage) rank-k update in the padded lanes, and those lanes are +simply never unpacked or read back. Because every arithmetic step is +unconditional (no pivot, no `larfg`-style mask), the padding cannot produce a +`NaN` that contaminates a real lane -- each lane is independent. This is the +`?syrk` analogue of the `potrf` observation that the identity flows through the +unconditional pivot path unmasked. + +### 6.5 No argument checking (Compact convention) + +Like MKL's own compact routines -- which "skip error checking for performance +reasons" and make "the user responsible for passing correct parameters" +(Intel MKL, *Numerical Limitations for Compact BLAS and Compact LAPACK +Routines*) -- `cqr_mkl_?syrk_compact` validates nothing. Because `?syrk` has no +`info`, there is not even a dispatch-level status: an unrecognized `format` is a +silent no-op. It is the caller's responsibility to pass valid parameters, and to +validate them beforehand if defensive checking is wanted -- for which the portable +`dsyrk_compact` / `ssyrk_compact` (`cqr_compact.h`) provide LAPACK/BLAS-style +`info = -j` argument validation. + +### 6.6 Numerical scope + +The rank-k update is a sequence of fused multiply-adds, backward stable to working +precision, matching reference BLAS `?syrk` element for element up to the +accumulation order (the dot-product summation rounds differently from, e.g., a +rank-1 accumulation, but both are backward stable). As with the rest of the +toolkit, no overflow/underflow-safe scaling is attempted, and complex precisions +(`c`/`z`, whose symmetric update is `?herk`) are out of scope. When +`cqr_mkl_?syrk_compact` feeds a Cholesky QR (section 1), the well-known +conditioning caveat of Cholesky QR applies: the orthogonality of the recovered +`Q` degrades like `cond(A)^2 * eps`, so for ill-conditioned `A` a reorthogonaliza- +tion pass (CholeskyQR2) or the QR-based pipeline is preferable. The rank-k update +itself is unaffected -- that is a property of the algorithm it feeds, not of +`?syrk`. + +## 7. Testing and Validation Methodology + +Matching standard BLAS `?syrk` to working precision is the minimum bar. SIMD, +blocking, and layout handling are internal strategies only: the returned `C` must +satisfy the same numerical invariants as an unbatched `?syrk`. Tolerances are +purely relative to the working precision, scaled by the accumulation length `k`. +All suites below are CTest-registered. + +### 7.1 Suite A -- vs dense per-matrix `cblas_?syrk` + +`test_cqr_syrk_mkl.cpp` packs a random `A` and a random (non-symmetric) `C` with +the genuine MKL Compact API, runs `cqr_mkl_?syrk_compact`, and compares the +unpacked result against a per-matrix `cblas_?syrk` from the same seed, element for +element over the *whole* `n x n` `C`. Because `?syrk` writes only the `uplo` +triangle, this simultaneously gates the active-triangle result and that the +opposite triangle is left untouched. Run over the full feature matrix +(`precision x layout x uplo x trans`), a spread of `alpha`/`beta` (including +`beta = 0` overwrite and `alpha = 0`), and shapes covering wide/tall factors, +`n = 1`, `k = 1`, `nm = 1`, and padded partial last groups. + +### 7.2 Suite B -- vs `mkl_?gemm_compact` (triangle only) + +The same packed inputs are run through `cqr_mkl_?syrk_compact` and, independently, +`mkl_?gemm_compact` forming the full product (`(NoTrans, Trans)` for `A A^T`, +`(Trans, NoTrans)` for `A^T A`). Since gemm writes the entire matrix while syrk +writes one triangle, the comparison is restricted to the active `uplo` triangle. +This cross-checks against a second, independent MKL kernel over the full feature +matrix -- the direct analogue of the `trsm`/`geqrf` suites' cross-checks against +their native MKL counterparts. + +### 7.3 Suite C -- end-to-end Cholesky QR, no MKL compute kernel + +The capstone: for a tall, well-conditioned random `A`, the fully open pipeline + +``` +cqr_mkl_dsyrk_compact('U','T') // G = A^T A +cqr_mkl_dpotrf_compact('U') // G = R^T R +cqr_mkl_dtrsm_compact('R','U','N') // Q = A R^{-1} +``` + +must recover a valid QR factorization. Gate the reconstruction +`|| Q R - A ||_1 / ||A||_1` at `50 * n * eps` (backward stable, independent of +conditioning) and the orthogonality `|| Q^T Q - I ||` at `1e3 * n * eps` (small +for the tall, well-conditioned inputs used, per the `cond(A)^2` caveat of +section 6.6). This exercises `cqr_mkl_?syrk_compact` closing a complete batched +Cholesky QR with no MKL compute kernel -- the syrk analogue of the `trsm` suite's +end-to-end QR solve. + +### 7.4 Portable self-test (no BLAS) + +`test_cqr_syrk_compact.cpp` validates the templated kernel directly against a +scalar `?syrk` reference over the full `uplo x trans x layout` matrix, across +precisions, interleave widths, and padded final packs, plus the LAPACK/BLAS-style +argument validation of the portable C API. The reference accumulates the rank-k +product in a **different summation order** (rank-1 outer products, contraction +index outermost) than the kernel's inner-product form, so a bug common to both +cannot pass unseen while a correct kernel still agrees to working precision. All +four `(trans, layout)` combinations are covered so the tuned (`trans='T'`, +column-major) and strided kernels are both exercised. It needs no external +libraries at all; only Suites A-C require an MKL installation (for the Compact +API). + +## 8. Implementation Strategy + +Modern C++ (C++17) templated on scalar type `T` and interleave width `V`, exposed +through `extern "C"` for the FFI-stable surfaces, reusing the existing +`cqr::detail::pack` / `BatchView` GNU-vector machinery. + +### 8.1 API boundary + +Two C-linkage interfaces wrap the same templated kernel: + +* **MKL-style API** (`cqr_mkl_ext.h`): `cqr_mkl_dsyrk_compact` / + `cqr_mkl_ssyrk_compact` -- taking the MKL enums and `MKL_COMPACT_PACK`, with no + argument checking (Compact convention). It **dispatches on the + `MKL_COMPACT_PACK` format enum** (`MKL_COMPACT_SSE`/`AVX`/`AVX512` -> + `V = 16/32/64 bytes / sizeof(T)`), a compile-time constant per case, exactly as + `cqr_mkl_?trsm_compact` / `cqr_mkl_?potrf_compact` do. +* **Portable C API** (`cqr_compact.h`): `dsyrk_compact` / `ssyrk_compact` -- an + MKL-independent surface taking an explicit interleave width `V` and `char` + selectors, with LAPACK/BLAS-style `info = -j` argument validation. + +Both instantiate the kernel on `MKL_INT` (or `int`) so ILP64 dimensions are not +narrowed; the algorithm and its internal routines are described in section 6. + +### 8.2 Source layout + +| File | Role | +|------|------| +| `src/cqr_syrk_compact.hpp` | Templated SIMD rank-k kernel: `syrk_compact_group` (tuned `trans='T'` column-major), `syrk_compact_group_strided` (general via `BatchView`), and the batch driver `syrk_compact_general` (scalar `T`, width `V`). | +| `src/cqr_syrk_compact_dispatch.cpp` | Portable `?syrk_compact` C entry points (runtime `V` -> compile-time dispatch, `info = -j`). | +| `src/cqr_mkl_syrk.cpp` | Unwraps the MKL enums + dispatches on `MKL_COMPACT_PACK` -> `V`, calls the kernel (drop-in style; no `work`/`info`). | +| `src/test_cqr_syrk_compact.cpp` | Self-contained correctness test vs a scalar `?syrk` reference (no BLAS). | +| `src/test_cqr_syrk_mkl.cpp` | MKL validation: vs `cblas_?syrk`, vs `mkl_?gemm_compact`, and the end-to-end Cholesky QR. | + +The `cqr_mkl_?syrk_compact` prototypes are added to `cqr_mkl_ext.h` and the +portable `?syrk_compact` prototypes to `cqr_compact.h`. diff --git a/src/cqr_compact.h b/src/cqr_compact.h index ae31b80..c0f454b 100644 --- a/src/cqr_compact.h +++ b/src/cqr_compact.h @@ -10,6 +10,7 @@ * dormqr_compact / sormqr_compact -- apply Q or Q^T from the left, B := op(Q) B * dpotrf_compact / spotrf_compact -- Cholesky factorization A = L L^T or A = U^T U * dtrsm_compact / strsm_compact -- triangular solve op(A) X = alpha B, etc. + * dsyrk_compact / ssyrk_compact -- symmetric rank-k update C := alpha A op(A) + beta C * * Compact layout; group g = idx/V, slot v = idx%V: * A_v(i,j) = ap [ g*ldap*ncol*V + (j*ldap + i)*V + v ] (column-major) @@ -131,6 +132,39 @@ int strsm_compact(char layout, char side, char uplo, char transa, char diag, int float alpha, const float *ap, int ldap, float *bp, int ldbp, int V, int nm); +/* Symmetric rank-k update -- the portable form of mkl_?syrk_compact, forming for + * every matrix in the batch, in place, + * C := alpha op(A) op(A)^T + beta C, op(A) = A ('N', A is n x k) + * or A^T ('T'/'C', A is k x n), + * with C the symmetric n x n result whose uplo triangle is the only part + * referenced and updated. Its headline use is the Gram matrix A^T A of a Cholesky + * QR (trans='T'), which then feeds ?potrf_compact and ?trsm_compact. + * layout 'C'/'c' column-major (tuned when trans='T') or 'R'/'r' row-major + * uplo 'U' update the upper triangle of C or 'L' the lower + * trans 'N' (A A^T) or 'T'/'C' (A^T A; 'C' == 'T' for the real types) + * n order of C (and the non-contracted dimension of A) + * k the contracted dimension (columns of A for 'N', rows for 'T') + * alpha scalar multiplying the rank-k product + * ap compact A (n x k for 'N', k x n for 'T') + * ldap compact leading dimension of A (>= max(1, rows(A)) col-major, + * >= max(1, cols(A)) row-major) + * beta scalar multiplying C; beta = 0 overwrites C (its prior value, even + * NaN, is not read) + * cp compact symmetric C (n x n); its uplo triangle is overwritten + * ldcp compact leading dimension of C (>= max(1, n)) + * V, nm interleave width; total number of matrices (padded last group) + * Returns 0, or -j for an illegal j-th argument: + * -1 layout -2 uplo -3 trans -4 n (<0) -5 k (<0) -8 ldap + * -11 ldcp -12 V (not 2/4/8/16) -13 nm (<0) + * (alpha, beta, ap and cp are never inspected, matching LAPACK/BLAS.) */ +int dsyrk_compact(char layout, char uplo, char trans, int n, int k, double alpha, + const double *ap, int ldap, double beta, double *cp, int ldcp, int V, + int nm); + +int ssyrk_compact(char layout, char uplo, char trans, int n, int k, float alpha, + const float *ap, int ldap, float beta, float *cp, int ldcp, int V, + int nm); + #ifdef __cplusplus } #endif diff --git a/src/cqr_mkl_ext.h b/src/cqr_mkl_ext.h index fb1b6ee..47917e4 100644 --- a/src/cqr_mkl_ext.h +++ b/src/cqr_mkl_ext.h @@ -8,6 +8,7 @@ * cqr_mkl_?ormqr_compact -- apply Q (or Q^T) of a Compact-format QR * cqr_mkl_?potrf_compact -- Cholesky factorization of an SPD Compact-format batch * cqr_mkl_?trsm_compact -- triangular solve op(A) X = alpha B (and variants) + * cqr_mkl_?syrk_compact -- symmetric rank-k update C := alpha A op(A) + beta C * * All use MKL's MKL_LAYOUT + MKL_COMPACT_PACK interface, so they drop into the * MKL compact ecosystem, but are backed by this project's own portable SIMD @@ -23,9 +24,12 @@ * they mix freely with MKL's native compact routines. cqr_mkl_?trsm_compact is a * portable, open alternative to mkl_?trsm_compact (identical signature), the step * that closes a batched solve, so it runs end to end with no MKL compute kernel. - * The API mirrors MKL's native compact ecosystem (MKL_LAYOUT + MKL_COMPACT_PACK); - * see the full parameter reference in the per-routine design docs - * (cqr_mkl_d{geqrf,ormqr,potrf,trsm}_compact_design.md). + * cqr_mkl_?syrk_compact is the missing mkl_?syrk_compact: the symmetric rank-k + * update, which forms the Gram matrix A^T A of a Cholesky QR (feeding + * cqr_mkl_?potrf_compact then cqr_mkl_?trsm_compact) at half the flops of a + * general gemm. The API mirrors MKL's native compact ecosystem + * (MKL_LAYOUT + MKL_COMPACT_PACK); see the full parameter reference in the + * per-routine design docs (cqr_mkl_d{geqrf,ormqr,potrf,trsm,syrk}_compact_design.md). * * Typical use -- the batched AX = B solver (now MKL-compute-free): * cqr_mkl_dgeqrf_compact (..., A -> H, tau); // A = Q R @@ -144,6 +148,34 @@ void cqr_mkl_strsm_compact(MKL_LAYOUT layout, MKL_SIDE side, MKL_UPLO uplo, float alpha, const float *ap, MKL_INT ldap, float *bp, MKL_INT ldbp, MKL_COMPACT_PACK format, MKL_INT nm); +/* Symmetric rank-k update. For every matrix in the batch, forms in place + * + * C := alpha * A * A^T + beta * C (trans = MKL_NOTRANS, A is n x k) or + * C := alpha * A^T * A + beta * C (trans = MKL_TRANS, A is k x n), + * + * where C is the symmetric n x n result and only its uplo triangle (MKL_UPPER or + * MKL_LOWER) is referenced and updated. op(A) = A^H (MKL_CONJTRANS) folds to A^T + * for the real types (the conjugated update is ?herk, out of scope here). As in + * BLAS ?syrk, beta = 0 overwrites C (its prior contents, even NaN, are not read). + * + * The Compact-format counterpart of BLAS ?syrk, completing the compact BLAS-3 set + * alongside mkl_?gemm_compact and mkl_?trsm_compact; unlike forming A*op(A) with + * mkl_?gemm_compact it exploits the symmetry (half the flops, a triangle of + * writes). Its headline use is the Gram matrix of a Cholesky QR: A^T A -> potrf + * -> trsm. Like the BLAS ?syrk it batches, it takes no workspace and reports no + * info, and does no argument checking (Compact convention) -- use dsyrk_compact / + * ssyrk_compact (cqr_compact.h) for LAPACK/BLAS-style validation. See + * cqr_mkl_dsyrk_compact_design.md. */ +void cqr_mkl_dsyrk_compact(MKL_LAYOUT layout, MKL_UPLO uplo, MKL_TRANSPOSE trans, + MKL_INT n, MKL_INT k, double alpha, const double *ap, + MKL_INT ldap, double beta, double *cp, MKL_INT ldcp, + MKL_COMPACT_PACK format, MKL_INT nm); + +void cqr_mkl_ssyrk_compact(MKL_LAYOUT layout, MKL_UPLO uplo, MKL_TRANSPOSE trans, + MKL_INT n, MKL_INT k, float alpha, const float *ap, + MKL_INT ldap, float beta, float *cp, MKL_INT ldcp, + MKL_COMPACT_PACK format, MKL_INT nm); + #ifdef __cplusplus } #endif diff --git a/src/cqr_mkl_syrk.cpp b/src/cqr_mkl_syrk.cpp new file mode 100644 index 0000000..e774617 --- /dev/null +++ b/src/cqr_mkl_syrk.cpp @@ -0,0 +1,80 @@ +/* cqr_mkl_syrk.cpp + * + * cqr_mkl_?syrk_compact (design doc section 8.1): a thin C-linkage adapter that + * unwraps the MKL enums to bool flags and MKL_COMPACT_PACK to the interleave + * width V, then forwards to cqr::detail::syrk_compact_general (on + * MKL_INT so ILP64 dimensions are not narrowed). + * + * The Compact-format counterpart of BLAS ?syrk, completing the compact BLAS-3 + * set alongside mkl_?gemm_compact and mkl_?trsm_compact. Like the BLAS ?syrk it + * batches -- and like the compact BLAS-3 routines it mirrors -- it takes no + * workspace and reports no info, and (Compact convention) does not validate its + * arguments: use the portable dsyrk_compact / ssyrk_compact (cqr_compact.h) for + * LAPACK/BLAS-style info = -j checking. Semantics and compact storage layout are + * documented in cqr_syrk_compact.hpp. + * + * Assisted-by: Claude:claude-opus-4.8 + */ + +#include "cqr_mkl_ext.h" +#include "cqr_syrk_compact.hpp" + +namespace { + +template +void run(MKL_LAYOUT layout, MKL_UPLO uplo, MKL_TRANSPOSE trans, MKL_INT n, MKL_INT k, + T alpha, const T *ap, MKL_INT ldap, T beta, T *cp, MKL_INT ldcp, + MKL_COMPACT_PACK format, MKL_INT nm) +{ + /* Empty problem: no output to produce. n == 0 leaves nothing to update; + * nm == 0 means no matrices (also keeps the kernel's nm >= 1 invariant). + * k == 0 and alpha == 0 are NOT short-circuited: both still scale C by beta, + * which the kernel performs (an empty contraction yields the zero update). */ + if (n == 0 || nm == 0) return; + + const bool rowmajor = (layout == MKL_ROW_MAJOR); + const bool upper = (uplo == MKL_UPPER); + /* real ?syrk: A^T only -- MKL_CONJTRANS folds to MKL_TRANS (conjugation is + * ?herk). NOTRANS gives A*A^T, anything else A^T*A. */ + const bool tran = (trans != MKL_NOTRANS); + + /* Dispatch on the pack format; the interleave width V = (SIMD bytes)/sizeof(T) + * is a compile-time constant in each case (SSE = 16 B, AVX = 32 B, + * AVX512 = 64 B). An unrecognized format selects no kernel; with no info to + * report that is a silent no-op (BLAS syrk convention). */ + switch (format) { + case MKL_COMPACT_SSE: + cqr::detail::syrk_compact_general( + upper, tran, rowmajor, n, k, alpha, ap, ldap, beta, cp, ldcp, nm); + break; + case MKL_COMPACT_AVX: + cqr::detail::syrk_compact_general( + upper, tran, rowmajor, n, k, alpha, ap, ldap, beta, cp, ldcp, nm); + break; + case MKL_COMPACT_AVX512: + cqr::detail::syrk_compact_general( + upper, tran, rowmajor, n, k, alpha, ap, ldap, beta, cp, ldcp, nm); + break; + default: break; /* unrecognized pack format: no kernel, no-op */ + } +} + +} /* anonymous namespace */ + +extern "C" void cqr_mkl_dsyrk_compact(MKL_LAYOUT layout, MKL_UPLO uplo, + MKL_TRANSPOSE trans, MKL_INT n, MKL_INT k, + double alpha, const double *ap, MKL_INT ldap, + double beta, double *cp, MKL_INT ldcp, + MKL_COMPACT_PACK format, MKL_INT nm) +{ + run(layout, uplo, trans, n, k, alpha, ap, ldap, beta, cp, ldcp, format, nm); +} + +extern "C" void cqr_mkl_ssyrk_compact(MKL_LAYOUT layout, MKL_UPLO uplo, + MKL_TRANSPOSE trans, MKL_INT n, MKL_INT k, + float alpha, const float *ap, MKL_INT ldap, + float beta, float *cp, MKL_INT ldcp, + MKL_COMPACT_PACK format, MKL_INT nm) +{ + run(layout, uplo, trans, n, k, alpha, ap, ldap, beta, cp, ldcp, format, nm); +} diff --git a/src/cqr_syrk_compact.hpp b/src/cqr_syrk_compact.hpp new file mode 100644 index 0000000..f79d6da --- /dev/null +++ b/src/cqr_syrk_compact.hpp @@ -0,0 +1,264 @@ +/* cqr_syrk_compact.hpp + * + * Compact (interleaved-batch) symmetric rank-k update, templated on scalar type + * T and interleave width V: + * + * C := alpha * A * A^T + beta * C (trans='N', A is n x k) + * C := alpha * A^T * A + beta * C (trans='T', A is k x n) + * + * C is the symmetric n x n result; only the uplo triangle (lower or upper) is + * read and written. This is the missing mkl_?syrk_compact, in portable form -- + * the companion to mkl_?gemm_compact / mkl_?trsm_compact, exploiting the + * symmetry MKL's gemm_compact cannot (half the flops, a triangle of writes). Its + * headline use is the Gram matrix of a Cholesky QR: A^T A -> potrf -> trsm. + * + * syrk_compact_group() -- tuned trans='T' (A^T A), column-major path. + * syrk_compact_group_strided() -- any uplo / trans / layout, via BatchView. + * syrk_compact_general() -- all groups; the sole batch entry point. + * + * Algorithm: the dot-product form. Each output C(i,j) is a length-k inner + * product of two "vectors" of A along the contraction axis, accumulated in + * registers and combined with beta*C on store. Compact format stores element + * (i,j) of all V matrices contiguously, so the scalar code lifts verbatim with + * T -> V-wide vector, one lane per matrix (no data-dependent branch). Preferring + * the dot form over the rank-1 (outer-product) form writes C exactly once and + * blocks cleanly against the triangle: for a fixed row i, every column j in that + * row's triangle is a *full* length-k dot (no diagonal-corner peeling). Output + * columns are register-blocked (JB=4) so each A vector load is reused across four + * columns -- the same reuse the ormqr/trsm kernels apply across their columns. + * beta is folded in on store; when beta==0 the element is overwritten rather than + * read, so an uninitialised / NaN C is allowed, matching reference BLAS ?syrk. + * + * Compact storage convention (matches MKL Compact / mkl_?gepack_compact); group + * g = idx/V, slot v = idx%V, C the symmetric n x n batch. For the tuned trans='T' + * column-major path A is the (ldap, n) batch, so + * A_v(p,i) = ap[ g*ldap*n*V + (i*ldap + p)*V + v ] (column i is contiguous) + * C_v(i,j) = cp[ g*ldcp*n*V + (j*ldcp + i)*V + v ] + * The strided path expresses every other uplo/trans/layout solely by which + * physical stride is the n-index axis, the contraction axis, and the C + * row/column -- exactly as side='L'/'R' is a stride choice in the ormqr kernel. + * + * Scope: real types (s/d). Complex would be ?herk (conjugated), a separate + * routine, not a trans='C' mode here -- mirroring the real-only ormqr/trsm scope. + * + * Assisted-by: Claude:claude-opus-4.8 + */ + +#ifndef CQR_SYRK_COMPACT_HPP +#define CQR_SYRK_COMPACT_HPP + +#include "cqr_compact_common.hpp" /* pack, BatchView, make_view, make_const_view */ + +#include +#include +#include + +namespace cqr { +namespace detail { + +/* ------------------------------------------------------------------ */ +/* One group of V interleaved matrices: tuned trans='T', column-major. */ +/* */ +/* C := alpha A^T A + beta C, A the (ldap, n) batch, C the (ldcp, n) */ +/* batch. Column i of A starts at a_ + i*ldap and its k entries are */ +/* contiguous, so the length-k dot that forms C(i,j) streams both */ +/* operands down contiguous packs -- the Cholesky-QR Gram matrix in the */ +/* layout LAPACK produces its factors in. JB=4 reuses each A(:,i) load */ +/* across four output columns j. */ +/* ------------------------------------------------------------------ */ + +template +void syrk_compact_group(bool upper, Int n, Int k, T alpha, T beta, const T *a_, Int ldap, + T *c_, Int ldcp) +{ + using VT = typename pack::type; + static_assert(std::is_floating_point::value, + "syrk_compact is defined for real float/double"); + assert(ldap >= k && ldcp >= n && k >= 0); + + const VT *A = reinterpret_cast(a_); + VT *C = reinterpret_cast(c_); + /* beta==0 overwrites C (must not read it: it may be uninitialised/NaN). */ + const bool overwrite = (beta == T(0)); + + for (Int i = 0; i < n; ++i) { + const VT *Ai = A + static_cast(i) * ldap; /* column i of A */ + const Int jlo = upper ? i : 0; /* row i's triangle */ + const Int jhi = upper ? n : i + 1; + + Int j = jlo; + + /* main loop: 4 output columns at a time; A(:,i) loaded once, used 4x */ + for (; j + 4 <= jhi; j += 4) { + const VT *Aj0 = A + static_cast(j + 0) * ldap; + const VT *Aj1 = A + static_cast(j + 1) * ldap; + const VT *Aj2 = A + static_cast(j + 2) * ldap; + const VT *Aj3 = A + static_cast(j + 3) * ldap; + VT w0{}, w1{}, w2{}, w3{}; + for (Int p = 0; p < k; ++p) { + const VT aip = Ai[p]; + w0 += aip * Aj0[p]; + w1 += aip * Aj1[p]; + w2 += aip * Aj2[p]; + w3 += aip * Aj3[p]; + } + VT &c0 = C[i + static_cast(j + 0) * ldcp]; + VT &c1 = C[i + static_cast(j + 1) * ldcp]; + VT &c2 = C[i + static_cast(j + 2) * ldcp]; + VT &c3 = C[i + static_cast(j + 3) * ldcp]; + if (overwrite) { + c0 = alpha * w0; + c1 = alpha * w1; + c2 = alpha * w2; + c3 = alpha * w3; + } + else { + c0 = alpha * w0 + beta * c0; + c1 = alpha * w1 + beta * c1; + c2 = alpha * w2 + beta * c2; + c3 = alpha * w3 + beta * c3; + } + } + + /* remainder columns */ + for (; j < jhi; ++j) { + const VT *Aj = A + static_cast(j) * ldap; + VT w{}; + for (Int p = 0; p < k; ++p) + w += Ai[p] * Aj[p]; + VT &cij = C[i + static_cast(j) * ldcp]; + cij = overwrite ? (alpha * w) : (alpha * w + beta * cij); + } + } +} + +/* ------------------------------------------------------------------ */ +/* One group, fully general: any uplo / trans / layout via BatchView. */ +/* */ +/* A is viewed as A(idx, p): idx is the C-index (n) axis, p the */ +/* contraction (k) axis -- so trans='N' vs 'T' is just which physical */ +/* stride plays each role, exactly as side='L'/'R' is in the ormqr */ +/* kernel. C is viewed as C(i, j) with its own row/column strides. The */ +/* dot-product math is identical to the tuned kernel above; only the */ +/* addressing, now carried by the two BatchViews, changes. */ +/* ------------------------------------------------------------------ */ + +template +void syrk_compact_group_strided(bool upper, Int n, Int k, T alpha, T beta, + BatchView::type, Int> A, + BatchView::type, Int> C) +{ + using VT = typename pack::type; + static_assert(std::is_floating_point::value, + "syrk_compact is defined for real float/double"); + /* strides must be non-degenerate so distinct indices map to distinct slots */ + assert(A.special && A.panel && C.special && C.panel); + + const bool overwrite = (beta == T(0)); + + for (Int i = 0; i < n; ++i) { + const Int jlo = upper ? i : 0; + const Int jhi = upper ? n : i + 1; + + Int j = jlo; + + /* main loop: 4 output columns at a time; A(i,:) loaded once, used 4x */ + for (; j + 4 <= jhi; j += 4) { + VT w0{}, w1{}, w2{}, w3{}; + for (Int p = 0; p < k; ++p) { + const VT aip = A(i, p); + w0 += aip * A(j + 0, p); + w1 += aip * A(j + 1, p); + w2 += aip * A(j + 2, p); + w3 += aip * A(j + 3, p); + } + VT &c0 = C(i, j + 0); + VT &c1 = C(i, j + 1); + VT &c2 = C(i, j + 2); + VT &c3 = C(i, j + 3); + if (overwrite) { + c0 = alpha * w0; + c1 = alpha * w1; + c2 = alpha * w2; + c3 = alpha * w3; + } + else { + c0 = alpha * w0 + beta * c0; + c1 = alpha * w1 + beta * c1; + c2 = alpha * w2 + beta * c2; + c3 = alpha * w3 + beta * c3; + } + } + + /* remainder columns */ + for (; j < jhi; ++j) { + VT w{}; + for (Int p = 0; p < k; ++p) + w += A(i, p) * A(j, p); + VT &cij = C(i, j); + cij = overwrite ? (alpha * w) : (alpha * w + beta * cij); + } + } +} + +/* ------------------------------------------------------------------ */ +/* All groups, fully general: uplo, trans in {N,T}, col- or row-major. */ +/* This is the sole batch entry point -- both C adapters call it for */ +/* every case, and it dispatches the tuned path itself (below). */ +/* */ +/* trans='N': A is n x k; trans='T': A is k x n. The trans='T', */ +/* column-major case routes to the tuned contiguous kernel; the other */ +/* three combinations use the strided kernel (same dot-product math, */ +/* correctness-first addressing). A padded partial last group is */ +/* processed too -- its padding slots carry whatever mkl_?gepack_compact */ +/* wrote, and their (garbage) C triangle is simply never read back. */ +/* ------------------------------------------------------------------ */ + +template +void syrk_compact_general(bool upper, bool trans, bool rowmajor, Int n, Int k, T alpha, + const T *ap, Int ldap, T beta, T *cp, Int ldcp, Int nm) +{ + assert(n >= 0 && k >= 0 && nm >= 1); + + /* A element strides (in VT units): one along the C-index (n) axis, one along + * the contraction (k) axis. trans='T' runs the n index down A's columns + * (A is k x n), trans='N' down its rows (A is n x k); column-major makes the + * first physical axis unit-stride, row-major the second. */ + const Int a_nidx = trans ? (rowmajor ? 1 : ldap) : (rowmajor ? ldap : 1); + const Int a_kidx = trans ? (rowmajor ? ldap : 1) : (rowmajor ? 1 : ldap); + + /* C is symmetric n x n: row index i, column index j. The uplo triangle is + * defined on the math indices (i,j) regardless of layout; only the strides + * differ. */ + const Int c_row = rowmajor ? ldcp : 1; + const Int c_col = rowmajor ? 1 : ldcp; + + /* group strides (scalar T units). A's packed per-matrix extent is ldap times + * the count of its non-leading axis; C's is ldcp*n. These span the whole + * batch, so widen to size_t before the product to avoid overflow. */ + const Int a_lines = rowmajor ? (trans ? k : n) : (trans ? n : k); + const std::size_t str_a = static_cast(ldap) * a_lines * V; + const std::size_t str_c = static_cast(ldcp) * n * V; + + /* trans='T', column-major (A^T A) routes to the tuned contiguous kernel; the + * other three trans/layout combinations use the strided kernel. The choice is + * loop-invariant across groups. */ + const bool tuned = (trans && !rowmajor); + const Int ngroups = (nm + V - 1) / V; + for (Int g = 0; g < ngroups; ++g) { + if (tuned) + syrk_compact_group(upper, n, k, alpha, beta, + ap + (std::size_t)g * str_a, ldap, + cp + (std::size_t)g * str_c, ldcp); + else + syrk_compact_group_strided( + upper, n, k, alpha, beta, + make_const_view(ap + (std::size_t)g * str_a, a_nidx, a_kidx), + make_view(cp + (std::size_t)g * str_c, c_row, c_col)); + } +} + +} /* namespace detail */ +} /* namespace cqr */ + +#endif /* CQR_SYRK_COMPACT_HPP */ diff --git a/src/cqr_syrk_compact_dispatch.cpp b/src/cqr_syrk_compact_dispatch.cpp new file mode 100644 index 0000000..64c4a81 --- /dev/null +++ b/src/cqr_syrk_compact_dispatch.cpp @@ -0,0 +1,107 @@ +/* cqr_syrk_compact_dispatch.cpp + * + * extern "C" wrappers around the templated compact symmetric rank-k update; + * validate the arguments LAPACK/BLAS-style and dispatch on the runtime interleave + * width V to a compile-time instantiation (cqr::detail::syrk_compact_general). + * + * Assisted-by: Claude:claude-opus-4.8 + */ + +#include "cqr_compact.h" +#include "cqr_syrk_compact.hpp" + +#include + +namespace { + +/* Validate the arguments LAPACK/BLAS-style and dispatch on the interleave width + * V. Returns 0 on success, or -j if the j-th argument (1-based, in signature + * order) had an illegal value. Scalar (alpha, beta) and pointer (ap, cp) + * arguments are not inspected, matching LAPACK/BLAS; never aborts the host + * process. */ +template +int dispatch(char layout, char uplo, char trans, int n, int k, T alpha, const T *ap, + int ldap, T beta, T *cp, int ldcp, int V, int nm) +{ + const bool col = (layout == 'C' || layout == 'c'); + const bool row = (layout == 'R' || layout == 'r'); + const bool up = (uplo == 'U' || uplo == 'u'); + const bool lo = (uplo == 'L' || uplo == 'l'); + const bool trans_ok = (trans == 'N' || trans == 'n' || trans == 'T' || trans == 't' || + trans == 'C' || trans == 'c'); + const bool tran = (trans == 'T' || trans == 't' || trans == 'C' || trans == 'c'); + + /* A is n x k (trans='N') or k x n (trans='T'); the leading dim bounds its + * stored-axis extent -- rows for column-major, columns for row-major. */ + const int arows = tran ? k : n; + const int acols = tran ? n : k; + const int ldamin = row ? (acols < 1 ? 1 : acols) : (arows < 1 ? 1 : arows); + const int ldcmin = (n < 1 ? 1 : n); /* C is n x n, layout-independent */ + + if (!col && !row) return -1; + if (!up && !lo) return -2; + if (!trans_ok) return -3; + if (n < 0) return -4; + if (k < 0) return -5; + /* -6 alpha, -9 beta: scalars; every value (including 0) is legal. */ + if (ldap < ldamin) return -8; + if (ldcp < ldcmin) return -11; + if (V != 2 && V != 4 && V != 8 && V != 16) return -12; + if (nm < 0) return -13; + + /* Nothing to produce for an empty problem (also keeps the kernel's nm >= 1 + * invariant satisfied below). k == 0 is NOT empty: it means C := beta*C (an + * empty contraction yields the zero update), so it flows to the kernel like + * any other value, as does alpha == 0. */ + if (n == 0 || nm == 0) return 0; + + /* Non-empty problem: the buffers are about to be dereferenced. LAPACK does + * not inspect pointers, and neither does a release build, but a debug assert + * catches an accidental null before it becomes a wild write. */ + assert(ap != nullptr && cp != nullptr); + + const bool rowmajor = row; + const bool upper = up; + + /* V is the compact interleave width, not necessarily one hardware register: + * pack is a GNU vector the compiler maps to registers or short unrolled + * bursts, so every width is valid for both types. MKL's format -> V mapping + * only ever selects 2/4/8 for double and 4/8/16 for float. */ + switch (V) { + case 2: + cqr::detail::syrk_compact_general(upper, tran, rowmajor, n, k, alpha, ap, + ldap, beta, cp, ldcp, nm); + break; + case 4: + cqr::detail::syrk_compact_general(upper, tran, rowmajor, n, k, alpha, ap, + ldap, beta, cp, ldcp, nm); + break; + case 8: + cqr::detail::syrk_compact_general(upper, tran, rowmajor, n, k, alpha, ap, + ldap, beta, cp, ldcp, nm); + break; + case 16: + cqr::detail::syrk_compact_general(upper, tran, rowmajor, n, k, alpha, ap, + ldap, beta, cp, ldcp, nm); + break; + } + return 0; +} + +} /* anonymous namespace */ + +extern "C" int dsyrk_compact(char layout, char uplo, char trans, int n, int k, + double alpha, const double *ap, int ldap, double beta, + double *cp, int ldcp, int V, int nm) +{ + return dispatch(layout, uplo, trans, n, k, alpha, ap, ldap, beta, cp, ldcp, V, + nm); +} + +extern "C" int ssyrk_compact(char layout, char uplo, char trans, int n, int k, + float alpha, const float *ap, int ldap, float beta, + float *cp, int ldcp, int V, int nm) +{ + return dispatch(layout, uplo, trans, n, k, alpha, ap, ldap, beta, cp, ldcp, V, + nm); +} diff --git a/src/test_compact_util.hpp b/src/test_compact_util.hpp index a50e731..c22c09a 100644 --- a/src/test_compact_util.hpp +++ b/src/test_compact_util.hpp @@ -42,6 +42,15 @@ template double max_abs_diff(const T *a, const T *b, size_t n) return d; } +// max |a| over n elements (BLAS i?amax magnitude, without the index). +template double maxabs(const T *a, size_t n) +{ + double d = 0; + for (size_t i = 0; i < n; ++i) + d = std::max(d, (double)std::abs(a[i])); + return d; +} + // L1 (max column sum) norm of a column-major m x n matrix. Templated so the // portable FP32/FP64 suites can reuse it; existing double callers deduce // T = double and are unaffected. diff --git a/src/test_cqr_syrk_compact.cpp b/src/test_cqr_syrk_compact.cpp new file mode 100644 index 0000000..2485da7 --- /dev/null +++ b/src/test_cqr_syrk_compact.cpp @@ -0,0 +1,230 @@ +// test_cqr_syrk_compact.cpp +// +// Self-contained validation of the templated compact symmetric rank-k update +// (dsyrk_compact / ssyrk_compact), with no BLAS dependency. The reference is a +// scalar ?syrk implemented directly -- the same terms the vectorized kernel +// accumulates V lanes at a time, but summed in a deliberately different order +// (rank-1 outer products, p outermost) than the kernel's per-element dot (p +// innermost), so a shared arithmetic bug cannot pass unseen while a correct +// kernel still agrees to working precision. +// +// Two parts: +// 1. C API argument validation -- exercises the LAPACK/BLAS-style info = -j +// contract of the portable entry points. +// 2. Numerical correctness over the full uplo x trans x layout matrix, +// precisions, and interleave widths (including padded partial groups): +// the whole n x n C is compared, so the check simultaneously gates the +// active-triangle result and that the opposite triangle is left untouched. +// All four (trans, layout) combinations are covered so the tuned +// (trans='T', column-major) kernel and the strided kernel are both exercised. +// +// Assisted-by: Claude:claude-opus-4.8 + +#include +#include +#include +#include +#include + +#include "cqr_compact.h" // dsyrk_compact / ssyrk_compact (C entry points) +#include "test_compact_util.hpp" // rng/frand, max_abs_diff, norm1, MatrixBatch, pack/unpack + +using namespace cqr::test; + +// ----------------------- reference kernel (scalar) ------------------ +// Dense column-major BLAS ?syrk: C := alpha op(A) op(A)^T + beta C, writing only +// the uplo triangle of the symmetric n x n C. op(A) = A (trans='N', A is n x k) +// or A^T (trans='T'/'C', A is k x n). The rank-k product is accumulated as a sum +// of rank-1 outer products (contraction index p outermost) -- a different +// summation order from the kernel's inner-product form, so the two round +// differently and a bug common to both cannot hide. beta = 0 overwrites (the +// prior C is not read), matching the kernel and reference BLAS ?syrk. + +template +static void ref_syrk(char uplo, char trans, int n, int k, T alpha, const T *A, int lda, + T beta, T *C, int ldc) +{ + const bool upper = (uplo == 'U' || uplo == 'u'); + const bool tran = (trans == 'T' || trans == 't' || trans == 'C' || trans == 'c'); + // op(A)(i,p): A(i,p) for 'N', A(p,i) for 'T' (both from column-major A). + auto Aop = [&](int i, int p) -> T { + return tran ? A[p + (size_t)i * lda] : A[i + (size_t)p * lda]; + }; + auto Ce = [&](int i, int j) -> T & { return C[i + (size_t)j * ldc]; }; + + // C := beta C (or 0) on the active triangle first, then accumulate. + for (int j = 0; j < n; ++j) { + const int ilo = upper ? 0 : j, ihi = upper ? j + 1 : n; + for (int i = ilo; i < ihi; ++i) + Ce(i, j) = (beta == T(0)) ? T(0) : beta * Ce(i, j); + } + for (int p = 0; p < k; ++p) + for (int j = 0; j < n; ++j) { + const int ilo = upper ? 0 : j, ihi = upper ? j + 1 : n; + const T ajp = alpha * Aop(j, p); /* loop-invariant across the i sweep */ + for (int i = ilo; i < ihi; ++i) + Ce(i, j) += ajp * Aop(i, p); + } +} + +// precision-overloaded shim: pick d/s by the pointer type +static int syrk_c(char lay, char up, char tr, int n, int k, double al, const double *a, + int lda, double be, double *c, int ldc, int V, int nm) +{ + return dsyrk_compact(lay, up, tr, n, k, al, a, lda, be, c, ldc, V, nm); +} +static int syrk_c(char lay, char up, char tr, int n, int k, float al, const float *a, + int lda, float be, float *c, int ldc, int V, int nm) +{ + return ssyrk_compact(lay, up, tr, n, k, al, a, lda, be, c, ldc, V, nm); +} + +// --------------------------- one numerical case --------------------- +// A is n x k (trans='N') or k x n (trans='T'); C is the symmetric n x n result. +// The MatrixBatch backing is column-major; pack_compact serializes it in the +// requested compact layout, so a single reference (column-major) covers both. + +template +static int run_case(bool rowmajor, char uplo, char trans, int nm, int n, int k) +{ + const bool tran = (trans == 'T'); + const T eps = std::numeric_limits::epsilon(); + const T alpha = T(0.5) + frand(); // non-trivial, non-zero coefficients + const T beta = T(0.25) + frand(); + + const int Arows = tran ? k : n, Acols = tran ? n : k; + const int ldAd = Arows; // dense (col-major) leading dim + const int ldAp = rowmajor ? Acols : Arows; // compact leading dim + const int ldC = n; + + // random A and a random (non-symmetric) C; syrk touches only the uplo + // triangle, so checking the whole matrix confirms the opposite triangle is + // left intact. + MatrixBatch A(nm, Arows, Acols), C(nm, n, n), Cref(nm, n, n); + for (int idx = 0; idx < nm; ++idx) { + for (size_t e = 0; e < (size_t)Arows * Acols; ++e) + A[idx][e] = frand(); + for (size_t e = 0; e < (size_t)n * n; ++e) + C[idx][e] = frand(); + std::copy(C[idx], C[idx] + (size_t)n * n, Cref[idx]); + ref_syrk(uplo, trans, n, k, alpha, A[idx], ldAd, beta, Cref[idx], ldC); + } + + // pack, run the routine under test, unpack + const int ng = (nm + V - 1) / V; + std::vector ap((size_t)ng * Arows * Acols * V), cp((size_t)ng * n * n * V); + pack_compact(A, ap.data(), ldAp, V, rowmajor); + pack_compact(C, cp.data(), ldC, V, rowmajor); + + const int info = syrk_c(rowmajor ? 'R' : 'C', uplo, trans, n, k, alpha, ap.data(), + ldAp, beta, cp.data(), ldC, V, nm); + + MatrixBatch Cout(nm, n, n); + unpack_compact(Cout, cp.data(), ldC, V, rowmajor); + + // whole-matrix relative error: active triangle correct + opposite untouched. + double worst = 0; + for (int idx = 0; idx < nm; ++idx) + worst = std::max(worst, max_abs_diff(Cout[idx], Cref[idx], (size_t)n * n) / + std::max(norm1(Cref[idx], n, n), 1e-300)); + + // The reference sums k terms in a different order than the kernel, so the two + // agree only to working precision (scaled by the accumulation length k). + const double rtol = 32.0 * (k + 1) * (double)eps; + const bool ok = (worst <= rtol) && (info == 0); + std::printf(" T=%-6s V=%-2d %s uplo=%c trans=%c nm=%-2d n=%-3d k=%-3d | " + "rel %.2e (rtol %.2e) info=%d %s\n", + sizeof(T) == 8 ? "double" : "float", V, rowmajor ? "row" : "col", uplo, + trans, nm, n, k, worst, rtol, info, ok ? "OK" : "FAIL"); + return ok ? 0 : 1; +} + +// --------------------- C API argument validation -------------------- + +static int test_validation() +{ + const int n = 6, k = 4, V = 4, nm = 4; + // buffers sized for the largest valid case (A up to 6x6, C is 6x6), zeroed so + // the valid calls that actually run the kernel produce a well-defined result. + std::vector A((size_t)n * n * V, 0), C((size_t)n * n * V, 0); + auto call = [&A, &C](char lay, char up, char tr, int n_, int k_, int ldap_, int ldcp_, + int V_, int nm_) { + return dsyrk_compact(lay, up, tr, n_, k_, 1.0, A.data(), ldap_, 0.0, C.data(), + ldcp_, V_, nm_); + }; + // clang-format off + struct { const char *what; int got, want; } t[] = { + {"valid N col", call('C','U','N', n, k, n, n, V, nm), 0}, + {"valid T col", call('C','L','T', n, k, k, n, V, nm), 0}, + {"valid N row", call('R','U','N', n, k, k, n, V, nm), 0}, + {"valid T row", call('R','L','T', n, k, n, n, V, nm), 0}, + {"trans=C", call('C','U','C', n, k, k, n, V, nm), 0}, + {"k=0 (C:=bC)", call('C','U','N', n, 0, n, n, V, nm), 0}, + {"bad layout", call('X','U','N', n, k, n, n, V, nm), -1}, + {"bad uplo", call('C','X','N', n, k, n, n, V, nm), -2}, + {"bad trans", call('C','U','X', n, k, n, n, V, nm), -3}, + {"n<0", call('C','U','N', -1, k, n, n, V, nm), -4}, + {"k<0", call('C','U','N', n, -1, n, n, V, nm), -5}, + {"ldap(rowmajor, uplo, trans, 8, 9, 5); + + // precisions, widths, tall/wide factors, and padded partial final groups + fails += run_case(false, 'L', 'T', 6, 12, 4); + fails += run_case(false, 'U', 'T', 16, 10, 20); // wide (k > n) + fails += run_case(false, 'L', 'N', 11, 12, 3); // padded last group + fails += run_case(true, 'U', 'N', 11, 7, 9); // row-major, padded + fails += run_case(false, 'U', 'T', 16, 16, 4); + fails += run_case(false, 'L', 'T', 32, 10, 7); + fails += run_case(true, 'L', 'N', 32, 8, 6); // row-major strided + + // corner cases: single matrix, n = 1, k = 1 (rank-1), n = k = JB width + for (char trans : {'N', 'T'}) + for (char uplo : {'U', 'L'}) { + fails += run_case(false, uplo, trans, 1, 5, 3); // nm = 1 + fails += run_case(false, uplo, trans, 5, 1, 4); // n = 1 + fails += run_case(false, uplo, trans, 4, 6, 1); // k = 1 + fails += run_case(false, uplo, trans, 8, 4, 4); // n = k = 4 + } + + if (fails) { + std::printf("\n%d CHECK(S) FAILED\n", fails); + return 1; + } + std::printf("\nall checks passed\n"); + return 0; +} diff --git a/src/test_cqr_syrk_mkl.cpp b/src/test_cqr_syrk_mkl.cpp new file mode 100644 index 0000000..0eacd7e --- /dev/null +++ b/src/test_cqr_syrk_mkl.cpp @@ -0,0 +1,426 @@ +/* test_cqr_syrk_mkl.cpp + * + * Validation of cqr_mkl_?syrk_compact against real Intel MKL, through the genuine + * MKL Compact pipeline (mkl_?gepack_compact / mkl_?geunpack_compact), in three + * independent ways over the full feature matrix (precision x layout x uplo x + * trans) and a range of alpha/beta and batch shapes (including padded partial + * last groups). See cqr_mkl_dsyrk_compact_design.md. + * + * Suite A (design doc 7.1) -- vs dense per-matrix cblas_?syrk: + * The compact result is checked, element for element over the *whole* n x n C, + * against a per-matrix reference from the same seed. Because ?syrk writes only + * the uplo triangle, this simultaneously gates the active-triangle result and + * that the opposite triangle is left untouched. + * Suite B (design doc 7.2) -- vs mkl_?gemm_compact (same compact pipeline): + * The full product alpha*A*op(A) + beta*C is formed with mkl_?gemm_compact into + * a second packed C; since gemm writes the whole matrix and syrk only one + * triangle, the comparison is restricted to the active uplo triangle. + * Suite C (design doc 7.3) -- end-to-end Cholesky QR, no MKL compute kernel: + * G = A^T A (cqr_mkl_?syrk_compact) -> R = chol(G) (cqr_mkl_?potrf_compact) -> + * Q = A R^{-1} (cqr_mkl_?trsm_compact) must satisfy Q R = A and Q^T Q = I. + * + * Build: needs Intel MKL; wired up by CMakeLists.txt. + * + * Assisted-by: Claude:claude-opus-4.8 + */ + +#include +#include + +#include "cqr_mkl_ext.h" +#include "cqr_mkl_alloc.h" /* mkl_alloc_bytes (calls mkl_malloc; links MKL) */ +#include "test_compact_util.hpp" /* rng/frand, max_abs_diff, norm1, batch_ptrs */ + +#include +#include +#include +#include +#include +#include + +using namespace cqr::test; + +namespace { + +/* max_abs_diff, norm1, maxabs, batch_ptrs and frand come from + * test_compact_util.hpp (shared across the compact test suites). */ + +/* ------------------------------------------------------------------ */ +/* Precision-dispatched MKL / CBLAS wrappers: overload on the scalar */ +/* pointer type, or specialize the size query (which carries no pointer).*/ +/* ------------------------------------------------------------------ */ + +template +MKL_INT compact_size(MKL_INT m, MKL_INT n, MKL_COMPACT_PACK f, MKL_INT nm); +template <> +MKL_INT compact_size(MKL_INT m, MKL_INT n, MKL_COMPACT_PACK f, MKL_INT nm) +{ + return mkl_dget_size_compact(m, n, f, nm); +} +template <> +MKL_INT compact_size(MKL_INT m, MKL_INT n, MKL_COMPACT_PACK f, MKL_INT nm) +{ + return mkl_sget_size_compact(m, n, f, nm); +} + +void pack(MKL_LAYOUT l, MKL_INT m, MKL_INT n, const double *const *a, MKL_INT lda, + double *ap, MKL_INT ldap, MKL_COMPACT_PACK f, MKL_INT nm) +{ + mkl_dgepack_compact(l, m, n, a, lda, ap, ldap, f, nm); +} +void pack(MKL_LAYOUT l, MKL_INT m, MKL_INT n, const float *const *a, MKL_INT lda, + float *ap, MKL_INT ldap, MKL_COMPACT_PACK f, MKL_INT nm) +{ + mkl_sgepack_compact(l, m, n, a, lda, ap, ldap, f, nm); +} + +void unpack(MKL_LAYOUT l, MKL_INT m, MKL_INT n, double **a, MKL_INT lda, const double *ap, + MKL_INT ldap, MKL_COMPACT_PACK f, MKL_INT nm) +{ + mkl_dgeunpack_compact(l, m, n, a, lda, ap, ldap, f, nm); +} +void unpack(MKL_LAYOUT l, MKL_INT m, MKL_INT n, float **a, MKL_INT lda, const float *ap, + MKL_INT ldap, MKL_COMPACT_PACK f, MKL_INT nm) +{ + mkl_sgeunpack_compact(l, m, n, a, lda, ap, ldap, f, nm); +} + +void gemm_compact(MKL_LAYOUT l, MKL_TRANSPOSE ta, MKL_TRANSPOSE tb, MKL_INT m, MKL_INT n, + MKL_INT k, double alpha, const double *ap, MKL_INT ldap, + const double *bp, MKL_INT ldbp, double beta, double *cp, MKL_INT ldcp, + MKL_COMPACT_PACK f, MKL_INT nm) +{ + mkl_dgemm_compact(l, ta, tb, m, n, k, alpha, ap, ldap, bp, ldbp, beta, cp, ldcp, f, + nm); +} +void gemm_compact(MKL_LAYOUT l, MKL_TRANSPOSE ta, MKL_TRANSPOSE tb, MKL_INT m, MKL_INT n, + MKL_INT k, float alpha, const float *ap, MKL_INT ldap, const float *bp, + MKL_INT ldbp, float beta, float *cp, MKL_INT ldcp, MKL_COMPACT_PACK f, + MKL_INT nm) +{ + mkl_sgemm_compact(l, ta, tb, m, n, k, alpha, ap, ldap, bp, ldbp, beta, cp, ldcp, f, + nm); +} + +/* dense per-matrix reference */ +void syrk(CBLAS_LAYOUT l, CBLAS_UPLO u, CBLAS_TRANSPOSE t, MKL_INT n, MKL_INT k, + double alpha, const double *a, MKL_INT lda, double beta, double *c, MKL_INT ldc) +{ + cblas_dsyrk(l, u, t, n, k, alpha, a, lda, beta, c, ldc); +} +void syrk(CBLAS_LAYOUT l, CBLAS_UPLO u, CBLAS_TRANSPOSE t, MKL_INT n, MKL_INT k, + float alpha, const float *a, MKL_INT lda, float beta, float *c, MKL_INT ldc) +{ + cblas_ssyrk(l, u, t, n, k, alpha, a, lda, beta, c, ldc); +} + +/* routine under test */ +void syrk_compact(MKL_LAYOUT l, MKL_UPLO u, MKL_TRANSPOSE t, MKL_INT n, MKL_INT k, + double alpha, const double *ap, MKL_INT ldap, double beta, double *cp, + MKL_INT ldcp, MKL_COMPACT_PACK f, MKL_INT nm) +{ + cqr_mkl_dsyrk_compact(l, u, t, n, k, alpha, ap, ldap, beta, cp, ldcp, f, nm); +} +void syrk_compact(MKL_LAYOUT l, MKL_UPLO u, MKL_TRANSPOSE t, MKL_INT n, MKL_INT k, + float alpha, const float *ap, MKL_INT ldap, float beta, float *cp, + MKL_INT ldcp, MKL_COMPACT_PACK f, MKL_INT nm) +{ + cqr_mkl_ssyrk_compact(l, u, t, n, k, alpha, ap, ldap, beta, cp, ldcp, f, nm); +} + +const char *pname(bool dbl) +{ + return dbl ? "d" : "s"; +} + +/* element (i,j) of a dense matrix stored in the given layout (ld = column stride + * col-major / row stride row-major) */ +template T elem(const T *p, int ld, int i, int j, bool rowmajor) +{ + return rowmajor ? p[(size_t)i * ld + j] : p[(size_t)j * ld + i]; +} + +/* max |X - Y| over the active uplo triangle of an n x n matrix (both stored in + * the same layout) */ +template +double tri_maxdiff(const T *X, const T *Y, int n, int ld, bool lower, bool rowmajor) +{ + double d = 0; + for (int i = 0; i < n; ++i) { + const int jlo = lower ? 0 : i, jhi = lower ? i + 1 : n; + for (int j = jlo; j < jhi; ++j) + d = std::max(d, (double)std::abs(elem(X, ld, i, j, rowmajor) - + elem(Y, ld, i, j, rowmajor))); + } + return d; +} + +/* ------------------------------------------------------------------ */ +/* Suite A: cqr_mkl_?syrk_compact vs per-matrix dense cblas_?syrk */ +/* ------------------------------------------------------------------ */ + +template +int suiteA(bool rowmajor, bool lower, bool trans, int nm, int n, int k, T alpha, T beta) +{ + const MKL_COMPACT_PACK fmt = mkl_get_format_compact(); + const int V = cqr::detail::vlen_for_format(fmt); + const T eps = std::numeric_limits::epsilon(); + + const MKL_LAYOUT ml = rowmajor ? MKL_ROW_MAJOR : MKL_COL_MAJOR; + const CBLAS_LAYOUT cl = rowmajor ? CblasRowMajor : CblasColMajor; + const MKL_UPLO mu = lower ? MKL_LOWER : MKL_UPPER; + const CBLAS_UPLO cu = lower ? CblasLower : CblasUpper; + const MKL_TRANSPOSE mtr = trans ? MKL_TRANS : MKL_NOTRANS; + const CBLAS_TRANSPOSE ctr = trans ? CblasTrans : CblasNoTrans; + + /* A is n x k (notrans) or k x n (trans); the dense == compact leading dim is + * the stored-axis extent (rows col-major, cols row-major). */ + const int Arows = trans ? k : n, Acols = trans ? n : k; + const int ldA = rowmajor ? Acols : Arows; + const int ldC = n; /* C is n x n in both layouts */ + const size_t sA = (size_t)Arows * Acols, sC = (size_t)n * n; + + std::vector A(nm * sA), C(nm * sC), Cref(nm * sC), Cout(nm * sC); + std::generate(A.begin(), A.end(), frand); + std::generate(C.begin(), C.end(), frand); + Cref = C; + for (int v = 0; v < nm; ++v) + syrk(cl, cu, ctr, n, k, alpha, A.data() + v * sA, ldA, beta, Cref.data() + v * sC, + ldC); + + auto Ap = batch_ptrs(A.data(), nm, sA); + auto Cp = batch_ptrs(C.data(), nm, sC); + auto ap = cqr::detail::mkl_alloc_bytes(compact_size(Arows, Acols, fmt, nm)); + auto cp = cqr::detail::mkl_alloc_bytes(compact_size(n, n, fmt, nm)); + pack(ml, Arows, Acols, Ap.data(), ldA, ap.get(), ldA, fmt, nm); + pack(ml, n, n, Cp.data(), ldC, cp.get(), ldC, fmt, nm); + + syrk_compact(ml, mu, mtr, n, k, alpha, ap.get(), ldA, beta, cp.get(), ldC, fmt, nm); + + auto Op = batch_ptrs(Cout.data(), nm, sC); + unpack(ml, n, n, Op.data(), ldC, cp.get(), ldC, fmt, nm); + + double worst = 0; + for (int v = 0; v < nm; ++v) { + /* whole matrix: active triangle correct + opposite triangle intact */ + double rel = max_abs_diff(Cout.data() + v * sC, Cref.data() + v * sC, sC) / + std::max(maxabs(Cref.data() + v * sC, sC), 1e-300); + worst = std::max(worst, rel); + } + const double rtol = 32.0 * (k + 1) * (double)eps; + const bool ok = (worst <= rtol); + std::printf(" [A:cblas] %s%s uplo=%c trans=%c V=%-2d nm=%-2d n=%-3d k=%-3d " + "a=%+.1f b=%+.1f | rel %.2e (rtol %.1e) %s\n", + pname(std::is_same::value), rowmajor ? "/row" : "/col", + lower ? 'L' : 'U', trans ? 'T' : 'N', V, nm, n, k, (double)alpha, + (double)beta, worst, rtol, ok ? "OK" : "FAIL"); + return ok ? 0 : 1; +} + +/* ------------------------------------------------------------------ */ +/* Suite B: cqr_mkl_?syrk_compact vs mkl_?gemm_compact (triangle only) */ +/* ------------------------------------------------------------------ */ + +template +int suiteB(bool rowmajor, bool lower, bool trans, int nm, int n, int k, T alpha, T beta) +{ + const MKL_COMPACT_PACK fmt = mkl_get_format_compact(); + const int V = cqr::detail::vlen_for_format(fmt); + const T eps = std::numeric_limits::epsilon(); + + const MKL_LAYOUT ml = rowmajor ? MKL_ROW_MAJOR : MKL_COL_MAJOR; + const MKL_UPLO mu = lower ? MKL_LOWER : MKL_UPPER; + const MKL_TRANSPOSE mtr = trans ? MKL_TRANS : MKL_NOTRANS; + /* gemm forms the same product: A*A^T via (NoTrans, Trans); A^T*A via + * (Trans, NoTrans). */ + const MKL_TRANSPOSE ga = trans ? MKL_TRANS : MKL_NOTRANS; + const MKL_TRANSPOSE gb = trans ? MKL_NOTRANS : MKL_TRANS; + + const int Arows = trans ? k : n, Acols = trans ? n : k; + const int ldA = rowmajor ? Acols : Arows; + const int ldC = n; + const size_t sA = (size_t)Arows * Acols, sC = (size_t)n * n; + + std::vector A(nm * sA), C(nm * sC); + std::generate(A.begin(), A.end(), frand); + std::generate(C.begin(), C.end(), frand); + + auto Ap = batch_ptrs(A.data(), nm, sA); + auto Cp = batch_ptrs(C.data(), nm, sC); + auto ap = cqr::detail::mkl_alloc_bytes(compact_size(Arows, Acols, fmt, nm)); + auto cs = cqr::detail::mkl_alloc_bytes(compact_size(n, n, fmt, nm)); /* syrk */ + auto cg = cqr::detail::mkl_alloc_bytes(compact_size(n, n, fmt, nm)); /* gemm */ + pack(ml, Arows, Acols, Ap.data(), ldA, ap.get(), ldA, fmt, nm); + pack(ml, n, n, Cp.data(), ldC, cs.get(), ldC, fmt, nm); + pack(ml, n, n, Cp.data(), ldC, cg.get(), ldC, fmt, nm); + + syrk_compact(ml, mu, mtr, n, k, alpha, ap.get(), ldA, beta, cs.get(), ldC, fmt, nm); + gemm_compact(ml, ga, gb, n, n, k, alpha, ap.get(), ldA, ap.get(), ldA, beta, cg.get(), + ldC, fmt, nm); + + std::vector Sout(nm * sC), Gout(nm * sC); + auto Sp = batch_ptrs(Sout.data(), nm, sC); + auto Gp = batch_ptrs(Gout.data(), nm, sC); + unpack(ml, n, n, Sp.data(), ldC, cs.get(), ldC, fmt, nm); + unpack(ml, n, n, Gp.data(), ldC, cg.get(), ldC, fmt, nm); + + double worst = 0; + for (int v = 0; v < nm; ++v) { + const T *sp = Sout.data() + v * sC, *gp = Gout.data() + v * sC; + /* compare only the triangle syrk wrote (gemm filled the whole matrix) */ + double diff = tri_maxdiff(sp, gp, n, ldC, lower, rowmajor); + worst = std::max(worst, diff / std::max(maxabs(gp, sC), 1e-300)); + } + const double rtol = 32.0 * (k + 1) * (double)eps; + const bool ok = (worst <= rtol); + std::printf(" [B:gemm ] %s%s uplo=%c trans=%c V=%-2d nm=%-2d n=%-3d k=%-3d " + "a=%+.1f b=%+.1f | rel %.2e (rtol %.1e) %s\n", + pname(std::is_same::value), rowmajor ? "/row" : "/col", + lower ? 'L' : 'U', trans ? 'T' : 'N', V, nm, n, k, (double)alpha, + (double)beta, worst, rtol, ok ? "OK" : "FAIL"); + return ok ? 0 : 1; +} + +/* ------------------------------------------------------------------ */ +/* Suite C: end-to-end Cholesky QR, no MKL compute kernel */ +/* G = A^T A (syrk, upper) -> R = chol(G) (potrf, upper) -> */ +/* Q = A R^{-1} (trsm, right/upper) ; check Q R = A and Q^T Q = I. */ +/* ------------------------------------------------------------------ */ + +int suiteC(int nm, int m, int n) +{ + const MKL_COMPACT_PACK fmt = mkl_get_format_compact(); + const int V = cqr::detail::vlen_for_format(fmt); + const double eps = std::numeric_limits::epsilon(); + + const size_t sA = (size_t)m * n, sG = (size_t)n * n; + std::vector A(nm * sA); + /* tall random A (m >= n) is well-conditioned */ + std::generate(A.begin(), A.end(), frand); + + auto Ap = batch_ptrs(A.data(), nm, sA); + auto a_buf = + cqr::detail::mkl_alloc_bytes(mkl_dget_size_compact(m, n, fmt, nm)); + auto g_buf = + cqr::detail::mkl_alloc_bytes(mkl_dget_size_compact(n, n, fmt, nm)); + double *ap = a_buf.get(), *gp = g_buf.get(); + mkl_dgepack_compact(MKL_COL_MAJOR, m, n, Ap.data(), m, ap, m, fmt, nm); + + MKL_INT info = 99; + + /* 1. Gram matrix G = A^T A (upper triangle) via our compact syrk */ + cqr_mkl_dsyrk_compact(MKL_COL_MAJOR, MKL_UPPER, MKL_TRANS, n, m, 1.0, ap, m, 0.0, gp, + n, fmt, nm); + + /* 2. R = chol(G): G = R^T R, R upper (potrf reads the upper triangle syrk wrote) */ + cqr_mkl_dpotrf_compact(MKL_COL_MAJOR, MKL_UPPER, n, gp, n, &info, fmt, nm); + + /* 3. Q = A R^{-1}: solve Q R = A (right side, upper, no-trans), A -> Q in place */ + cqr_mkl_dtrsm_compact(MKL_COL_MAJOR, MKL_RIGHT, MKL_UPPER, MKL_NOTRANS, MKL_NONUNIT, + m, n, 1.0, gp, n, ap, m, fmt, nm); + + std::vector Q(nm * sA), R(nm * sG); + auto Qp = batch_ptrs(Q.data(), nm, sA); + auto Rp = batch_ptrs(R.data(), nm, sG); + mkl_dgeunpack_compact(MKL_COL_MAJOR, m, n, Qp.data(), m, ap, m, fmt, nm); + mkl_dgeunpack_compact(MKL_COL_MAJOR, n, n, Rp.data(), n, gp, n, fmt, nm); + + double worst_recon = 0, worst_orth = 0; + std::vector QR(sA), QtQ(sG); + for (int v = 0; v < nm; ++v) { + const double *Av = A.data() + v * sA, *Qv = Q.data() + v * sA; + const double *Rv = R.data() + v * sG; + /* reconstruction Q R (R upper triangular) vs the original A */ + for (int j = 0; j < n; ++j) + for (int i = 0; i < m; ++i) { + double s = 0; + for (int l = 0; l <= j; ++l) /* R upper: R(l,j) nonzero for l <= j */ + s += Qv[i + (size_t)l * m] * Rv[l + (size_t)j * n]; + QR[i + (size_t)j * m] = s; + } + worst_recon = std::max(worst_recon, max_abs_diff(QR.data(), Av, sA) / + std::max(norm1(Av, m, n), 1e-300)); + /* orthogonality Q^T Q vs I */ + for (int j = 0; j < n; ++j) + for (int i = 0; i < n; ++i) { + double s = 0; + for (int l = 0; l < m; ++l) + s += Qv[l + (size_t)i * m] * Qv[l + (size_t)j * m]; + QtQ[i + (size_t)j * n] = s - (i == j ? 1.0 : 0.0); + } + worst_orth = std::max(worst_orth, maxabs(QtQ.data(), sG)); + } + /* reconstruction is backward stable (independent of conditioning); the + * Cholesky-QR orthogonality error grows like cond(A)^2 * eps, which stays + * small for these tall, well-conditioned random inputs. */ + const double rtol_recon = 50.0 * n * eps; + const double rtol_orth = 1e3 * n * eps; + const bool ok = + (info == 0) && (worst_recon <= rtol_recon) && (worst_orth <= rtol_orth); + std::printf(" [C:cqr ] V=%-2d nm=%-2d m=%-3d n=%-3d | recon %.2e (rtol %.1e) " + "orth %.2e (rtol %.1e) info=%d %s\n", + V, nm, m, n, worst_recon, rtol_recon, worst_orth, rtol_orth, (int)info, + ok ? "OK" : "FAIL"); + return ok ? 0 : 1; +} + +/* Run Suites A and B over the full feature matrix for one precision T. */ +template int run_precision() +{ + int fails = 0; + /* (nm, n, k), spanning ordinary shapes and corner cases */ + const int shapes[][3] = { + {8, 9, 4}, /* baseline; one full group at V=8 (double) */ + {16, 5, 7}, /* wide factor (k > n), nm a multiple of V */ + {11, 12, 3}, /* tall factor (k < n), padded partial last group */ + {1, 5, 3}, /* nm = 1: smallest batch (a single partial group) */ + {5, 1, 4}, /* n = 1: degenerate single-element triangle */ + {4, 6, 1}, /* k = 1: rank-1 update */ + {17, 4, 4}, /* nm one past a full group; n = k = the JB=4 width */ + }; + /* (alpha, beta): identity, scaled-accumulate, beta=0 overwrite, alpha=0. */ + const T coeffs[][2] = { + {T(1), T(0)}, + {T(2), T(0.5)}, + {T(-1.5), T(1)}, + {T(0), T(0.5)}, + }; + + for (bool rowmajor : {false, true}) + for (bool lower : {false, true}) + for (bool trans : {false, true}) + for (auto &[nm, n, k] : shapes) + for (auto &[alpha, beta] : coeffs) { + fails += suiteA(rowmajor, lower, trans, nm, n, k, alpha, beta); + fails += suiteB(rowmajor, lower, trans, nm, n, k, alpha, beta); + } + return fails; +} + +} /* anonymous namespace */ + +int main() +{ + std::printf("MKL compact format = %d, V(double) = %d, V(float) = %d\n", + (int)mkl_get_format_compact(), + cqr::detail::vlen_for_format(mkl_get_format_compact()), + cqr::detail::vlen_for_format(mkl_get_format_compact())); + + int fails = 0; + std::printf("-- double (Suites A vs cblas, B vs gemm_compact) --\n"); + fails += run_precision(); + std::printf("-- single (Suites A vs cblas, B vs gemm_compact) --\n"); + fails += run_precision(); + + std::printf("-- Cholesky QR end-to-end (Suite C: syrk -> potrf -> trsm) --\n"); + fails += suiteC(8, 32, 5); + fails += suiteC(8, 64, 8); + fails += suiteC(7, 40, 6); /* padded partial last group */ + + if (fails) { + std::printf("\n%d CHECK(S) FAILED\n", fails); + return 1; + } + std::printf("\nall checks passed\n"); + return 0; +}