Refactor: extract shared compact kernel machinery to common header

CMakeLists.txt

@@ -30,7 +30,7 @@ endif()
 # Core kernel library (portable; no BLAS dependency)
 # ---------------------------------------------------------------------------
 add_library(cqr_compact
-  src/cqr_compact_dispatch.cpp
+  src/cqr_ormqr_compact_dispatch.cpp
   src/cqr_geqrf_compact_dispatch.cpp
   src/cqr_potrf_compact_dispatch.cpp
   src/cqr_trsm_compact_dispatch.cpp)
@@ -47,7 +47,7 @@ if(CQR_WITH_MKL)
   find_package(MKLCompact REQUIRED)
 
   add_library(cqr_mkl_ext
-    src/cqr_mkl_ext.cpp
+    src/cqr_mkl_ormqr.cpp
     src/cqr_mkl_geqrf.cpp
     src/cqr_mkl_potrf.cpp
     src/cqr_mkl_trsm.cpp)
@@ -98,10 +98,10 @@ endif()
 if(CQR_BUILD_TESTS)
   enable_testing()
 
-  # Portable self-contained kernel test (no BLAS).
-  add_executable(test_cqr_compact src/test_cqr_compact.cpp)
-  target_link_libraries(test_cqr_compact PRIVATE cqr_compact)
-  add_test(NAME portable_kernel COMMAND test_cqr_compact)
+  # Portable self-contained ormqr (apply-Q) test (no BLAS).
+  add_executable(test_cqr_ormqr_compact src/test_cqr_ormqr_compact.cpp)
+  target_link_libraries(test_cqr_ormqr_compact PRIVATE cqr_compact)
+  add_test(NAME portable_ormqr COMMAND test_cqr_ormqr_compact)
 
   # Portable self-contained QR factorization test (no BLAS): the compact geqrf
   # kernel vs a scalar geqr2/larfg reference, plus C API validation.
@@ -124,10 +124,10 @@ if(CQR_BUILD_TESTS)
 
   # MKL-backed validation of cqr_mkl_dormqr_compact (design doc section 7).
   if(CQR_WITH_MKL)
-    add_executable(test_cqr_mkl_ext src/test_cqr_mkl_ext.cpp)
-    target_link_libraries(test_cqr_mkl_ext
+    add_executable(test_cqr_ormqr_mkl src/test_cqr_ormqr_mkl.cpp)
+    target_link_libraries(test_cqr_ormqr_mkl
       PRIVATE cqr_mkl_ext MKL::Compact)
-    add_test(NAME mkl_compact_suites COMMAND test_cqr_mkl_ext)
+    add_test(NAME mkl_ormqr_suites COMMAND test_cqr_ormqr_mkl)
 
     # MKL + dense LAPACK validation of cqr_mkl_dgeqrf_compact (design doc
     # section 7): factorization invariants, cross-check vs mkl_dgeqrf_compact,

PLANS.md

@@ -42,7 +42,7 @@ batches. Status vs. its design document:
   Column-major lower is the tuned contiguous path; row-major upper folds onto it
   by transpose duality, and the other two `(layout, uplo)` combinations route
   through the stride-generalized kernel. The shared `vsqrt<T,V>` helper now lives
-  in `cqr_compact.hpp` alongside `pack`/`BatchView`.
+  in `cqr_compact_common.hpp` alongside `pack`/`BatchView`.
 - **Validated (design section 7):** a BLAS-free test vs. a scalar `potf2`
   reference (both `uplo`, all layouts, padded final packs, and a non-SPD
   lane-isolation case gating that a poisoned lane never contaminates its
@@ -118,7 +118,7 @@ Gaps between the `ormqr` implementation and its design document
 
 - **Complex precisions (`cunmqr`/`zunmqr`).** Only real precisions exist; the
   family is real-only and `trans='C'` is folded to `'T'`
-  (in `src/cqr_mkl_ext.cpp`). Document the real-only scope, or
+  (in `src/cqr_mkl_ormqr.cpp`). Document the real-only scope, or
   add genuine complex specializations.
 - **Stress-test matrix (sections 7.3/7.4) absent.** Tests use only
   well-conditioned `frand` + diagonal boost. Missing: the `cond` scaling knob
@@ -141,7 +141,7 @@ What the implementation provides, mapped to the design document:
 - **Dispatcher (section 8.1):** unwraps the format to the interleave width `V`
   (SSE/AVX/AVX-512 -> 2/4/8 for FP64, 4/8/16 for FP32) and forwards to the
   templated kernel.
-- **Validation (section 7):** `test_cqr_mkl_ext` runs Suite 1 (isolated
+- **Validation (section 7):** `test_cqr_ormqr_mkl` runs Suite 1 (isolated
   `op(Q)*C` vs dense LAPACK, gate `20*s*eps`) over the full feature matrix --
   `side in {L,R} x layout in {col,row} x trans in {N,T}` -- and Suite 2
   (end-to-end `AX=B`: `mkl_dgeqrf_compact -> cqr_mkl_dormqr_compact ->

README.md

@@ -98,26 +98,27 @@ not part of the supported interface:
 
 | File | Role |
 |------|------|
+| `src/cqr_compact_common.hpp` | Shared machinery for the compact kernels: the `pack<T,V>` packed vector type plus the `BatchView` / `vsqrt` / `broadcast` helpers, included by every routine header below. |
 | `src/cqr_geqrf_compact.hpp` | Templated SIMD QR-factorization kernel (vectorized `geqr2`; scalar `T`, interleave width `V`). |
 | `src/cqr_potrf_compact.hpp` | Templated SIMD Cholesky-factorization kernel (vectorized `potf2`; scalar `T`, interleave width `V`). |
-| `src/cqr_compact.hpp` | Templated SIMD kernel `B := op(Q)*B` (scalar `T`, interleave width `V`); also hosts the shared `pack<T,V>` / `BatchView` / `vsqrt` / `broadcast` machinery. |
+| `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_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_compact_dispatch.cpp` | Portable ormqr 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_mkl_geqrf.cpp` | Unwraps `MKL_COMPACT_PACK` -> `V` and calls the geqrf kernel. |
+| `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_ext.cpp` | Unwraps `MKL_COMPACT_PACK` -> `V` and calls the ormqr kernel. |
-| `src/cqr_mkl_trsm.cpp` | Unwraps the MKL enums + `MKL_COMPACT_PACK` -> `V` and calls the trsm kernel (drop-in for `mkl_?trsm_compact`; no `work`/`info`). |
+| `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_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). |
 | `src/test_cqr_geqrf_mkl.cpp` | MKL + dense-LAPACK validation of `cqr_mkl_dgeqrf_compact` (residual, orthogonality, solve). |
 | `src/test_cqr_potrf_compact.cpp` | Self-contained potrf correctness test vs a scalar `potf2` reference (no BLAS). |
 | `src/test_cqr_potrf_mkl.cpp` | MKL + dense-LAPACK validation of `cqr_mkl_dpotrf_compact` (residual, untouched triangle, uniqueness, cross-check, solve). |
-| `src/test_cqr_compact.cpp` | Self-contained ormqr correctness/bench test (no BLAS). |
-| `src/test_cqr_mkl_ext.cpp` | MKL-backed validation through the real compact pipeline. |
+| `src/test_cqr_ormqr_compact.cpp` | Self-contained ormqr correctness/bench test (no BLAS). |
+| `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/cqr_compact_common.hpp

@@ -0,0 +1,158 @@
+/* cqr_compact_common.hpp
+ *
+ * Shared machinery for the compact (interleaved-batch) kernels: the V-wide
+ * packed vector type and the small lane-wise / addressing helpers built on it.
+ * Included by every routine header (geqrf, ormqr, potrf, trsm) so the pack type
+ * and views are defined once:
+ *
+ *   pack<T,V>                    -- the V-wide SIMD element (GNU vector type).
+ *   vsqrt / broadcast            -- lane-wise helpers (geqrf/potrf norms and
+ *                                   pivots, trsm's alpha scaling).
+ *   BatchView                    -- a strided 2-D view of one group of V
+ *                                   interleaved matrices.
+ *   make_view / make_const_view  -- reinterpret a packed T buffer as that view.
+ *
+ * V is the compact-format interleave width (the number of matrices whose
+ * element (i,j) is stored contiguously). It does NOT need to match the hardware
+ * vector width:
+ *   x86:   V*sizeof(T) = 16/32/64 bytes maps exactly to XMM/YMM/ZMM.
+ *   NEON:  128-bit registers; V=4 or V=8 doubles lower to short unrolled bursts
+ *          of 2/4 independent fmla v*.2d chains, which wide cores (Apple
+ *          M-series) execute very well.
+ *   SVE:   compile fixed-width with -msve-vector-bits=512 on A64FX to map V=8
+ *          doubles onto one SVE register.
+ *
+ * Assisted-by: Claude:claude-fable-5 Claude:claude-opus-4.8
+ */
+
+#ifndef CQR_COMPACT_COMMON_HPP
+#define CQR_COMPACT_COMMON_HPP
+
+#include <cstddef>
+#include <cmath>
+
+namespace cqr {
+namespace detail {
+
+/* ------------------------------------------------------------------ */
+/* pack<T,V>::type : the V-wide SIMD element                          */
+/* ------------------------------------------------------------------ */
+
+/* Define the V-wide element as a GNU vector type when the compiler provides
+ * the vector_size and may_alias attributes, detected directly via
+ * __has_attribute (itself guarded for preprocessors that predate it). A
+ * compiler that supplies these attributes -- GCC, Clang, Intel icpx/icpc --
+ * uses the vector type; any other stops at the #error below. */
+#if defined(__has_attribute)
+#if __has_attribute(vector_size) && __has_attribute(__may_alias__)
+#define CQR_HAS_GNU_VECTORS 1
+#endif
+#endif
+
+#if defined(CQR_HAS_GNU_VECTORS)
+
+template <typename T, int V> struct pack {
+    /* GNU vector_size requires a power-of-two byte width; the supported
+     * interleave widths are 2/4/8/16, matching the C API. Check it here -- the
+     * single chokepoint -- so a bad width fails with this message instead of a
+     * cryptic error inside the attribute instantiation. */
+    static_assert(V == 2 || V == 4 || V == 8 || V == 16,
+                  "interleave width V must be 2, 4, 8, or 16");
+    /* aligned(alignof(T)) relaxes the alignment requirement so the type
+     * is valid on any T-aligned buffer (unaligned vector loads are free
+     * on all modern hardware); may_alias exempts it from strict-aliasing
+     * violations when viewing a plain T array. */
+    using type
+        __attribute__((vector_size(V * sizeof(T)), aligned(alignof(T)), may_alias)) = T;
+};
+
+#else
+#error "cqr compact kernels require the GNU vector extensions " \
+       "(__attribute__((vector_size)) with may_alias); compile with a " \
+       "compiler that supports them (GCC, Clang, Intel icpx/icpc). These " \
+       "attributes are available under strict -std=c++17, not only GNU mode."
+#endif
+
+/* ------------------------------------------------------------------ */
+/* Lane-wise vector helpers shared by the compact kernels.             */
+/*                                                                     */
+/* These V-wide helpers take and return their vectors by reference.    */
+/* Passing a GNU vector by value would, without -march, commit the     */
+/* base-ISA vector argument/return ABI, which GCC and Clang (rightly)  */
+/* flag via -Wpsabi; a reference is just a pointer, so there is no such */
+/* boundary -- and once inlined the codegen is identical -- keeping the */
+/* build warning-clean with no compiler flag. Results are written      */
+/* through an out-parameter (named first).                             */
+/* ------------------------------------------------------------------ */
+
+/* r := sqrt(x), lane-wise. The short loop lowers to one vsqrt* on GCC/Clang; it
+ * runs once per column, negligible next to the O(n^2)/O(n^3) vector arithmetic.
+ * Used by geqrf's larfg (column norm) and potrf's pivot. */
+template <typename T, int V>
+inline void vsqrt(typename pack<T, V>::type &r,
+                  const typename pack<T, V>::type &x) noexcept
+{
+    for (int v = 0; v < V; ++v)
+        r[v] = std::sqrt(x[v]);
+}
+
+/* v := x broadcast to all V lanes. GNU vector types broadcast a scalar in
+ * arithmetic but not in assignment (`v = x;` is a compile error); `x - VT{}`
+ * subtracts an all-zero vector, leaving x in every lane (and, unlike `VT{} + x`,
+ * it preserves the sign of a zero x). Used by trsm to scale B by alpha. */
+template <typename T, int V>
+inline void broadcast(typename pack<T, V>::type &v, T x) noexcept
+{
+    v = x - typename pack<T, V>::type{};
+}
+
+/* ------------------------------------------------------------------ */
+/* BatchView: a strided 2-D view of one group of V interleaved matrices*/
+/*                                                                     */
+/* The element type is the V-wide pack (use a const pack for read-only */
+/* operands such as the reflector batch A). Indices are in elements;   */
+/* strides are in units of the V-wide pack, so one BatchView addresses */
+/* element (i,p) of every matrix in the group at once. The two axes    */
+/* are named for their role in the reflector sweep, not for row/col:   */
+/*   special -- the axis the Householder vector runs along             */
+/*              (rows of A and, for side='L', of C; columns for 'R'),  */
+/*   panel   -- the orthogonal axis, register-blocked 4 at a time.     */
+/* ------------------------------------------------------------------ */
+
+template <typename VT, typename Int = int> struct BatchView {
+    VT *const data = nullptr;
+    const Int special = 0; /* stride along the swept (reflector) axis */
+    const Int panel = 0;   /* stride along the orthogonal panel axis  */
+
+    /* The batch and matrix dimensions fit in Int (compact targets many small
+     * matrices), so the element offset is formed in Int with no widening of the
+     * induction variables -- which keeps the strided sweep vectorizable. The
+     * per-group base offset, which can exceed Int, is applied to the pointer by
+     * the caller before the view is built. */
+    inline VT &operator()(Int i, Int p) const noexcept
+    {
+        return data[i * special + p * panel];
+    }
+};
+
+/* Reinterpret a packed T buffer as a group view of V-wide pack elements. The
+ * strides are per-matrix leading dimensions (in pack units), so they take the
+ * kernel's Int like the dimensions do. */
+template <typename T, int V, typename Int = int>
+BatchView<const typename pack<T, V>::type, Int> make_const_view(const T *p, Int special,
+                                                                Int panel) noexcept
+{
+    using VT = typename pack<T, V>::type;
+    return {reinterpret_cast<const VT *>(p), special, panel};
+}
+template <typename T, int V, typename Int = int>
+BatchView<typename pack<T, V>::type, Int> make_view(T *p, Int special, Int panel) noexcept
+{
+    using VT = typename pack<T, V>::type;
+    return {reinterpret_cast<VT *>(p), special, panel};
+}
+
+} /* namespace detail */
+} /* namespace cqr */
+
+#endif /* CQR_COMPACT_COMMON_HPP */

src/cqr_geqrf_compact.hpp

@@ -10,7 +10,7 @@
  * geqrf_compact_general() -- column- or row-major.
  *
  * A portable, vectorized mkl_?geqrf_compact. It is the factorization companion
- * to cqr_compact.hpp's ormqr_compact (which applies the reflectors produced
+ * to cqr_ormqr_compact.hpp's ormqr_compact (which applies the reflectors produced
  * here), and reuses the same pack<T,V> / BatchView machinery.
  *
  * Algorithm: the unblocked LAPACK geqr2 (dlarfg to build each reflector, dlarf
@@ -36,7 +36,7 @@
 #ifndef CQR_GEQRF_COMPACT_HPP
 #define CQR_GEQRF_COMPACT_HPP
 
-#include "cqr_compact.hpp" /* pack<T,V>, BatchView, make_view, make_const_view, vsqrt */
+#include "cqr_compact_common.hpp" /* pack<T,V>, BatchView, make_view, make_const_view, vsqrt */
 
 #include <cstddef>
 #include <cstdint>
@@ -319,8 +319,8 @@ void geqrf_compact_general(bool rowmajor, Int m, Int n, T *ap, Int ldap, T *taup
     const std::size_t str_t = static_cast<std::size_t>(k) * V;
 
     /* element strides (in VT units) and per-matrix group stride (in T units) */
-    const std::size_t a_special = rowmajor ? (std::size_t)ldap : 1; /* down a col */
-    const std::size_t a_panel = rowmajor ? 1 : (std::size_t)ldap;   /* across cols */
+    const Int a_special = rowmajor ? ldap : 1; /* down a col  */
+    const Int a_panel = rowmajor ? 1 : ldap;   /* across cols */
     const std::size_t str_a =
         (rowmajor ? (std::size_t)ldap * m : (std::size_t)ldap * n) * V;
 

src/cqr_mkl_ext.cpp → src/cqr_mkl_ormqr.cpp

@@ -1,4 +1,4 @@
-/* cqr_mkl_ext.cpp
+/* cqr_mkl_ormqr.cpp
  *
  * Implementation of cqr_mkl_?ormqr_compact (design document section 8.1):
  * a thin C-linkage adapter that
@@ -30,7 +30,7 @@
  */
 
 #include "cqr_mkl_ext.h"
-#include "cqr_compact.hpp"
+#include "cqr_ormqr_compact.hpp"
 
 namespace {