Add cqr_mkl_?syrk_compact: batched symmetric rank-k update

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)
 

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

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