From d6daf5890bd91d21b6d3afecfc23a116514ff232 Mon Sep 17 00:00:00 2001 From: Claude Date: Thu, 6 Aug 2026 13:37:14 +0000 Subject: [PATCH 1/4] refactor: rename ?ormqr files and split shared kernel helpers Rename the ormqr implementation files to the project's per-routine naming scheme, matching geqrf/potrf/trsm: cqr_compact.hpp -> cqr_ormqr_compact.hpp cqr_compact_dispatch.cpp -> cqr_ormqr_compact_dispatch.cpp cqr_mkl_ext.cpp -> cqr_mkl_ormqr.cpp test_cqr_compact.cpp -> test_cqr_ormqr_compact.cpp test_cqr_mkl_ext.cpp -> test_cqr_ormqr_mkl.cpp The shared C API headers (cqr_mkl_ext.h, cqr_compact.h) and the cqr_mkl_ext library target keep their umbrella names. Split the common helpers and the packed vector type out of the old cqr_compact.hpp into a new cqr_compact_common.hpp, shared by every routine header: pack, vsqrt, broadcast, BatchView, and the make_view/make_const_view factories. cqr_ormqr_compact.hpp now holds only the ormqr-specific Direction enum and kernels and includes the common header the same way geqrf/potrf/trsm already do. Update the CMake sources and ctest names (portable_kernel -> portable_ormqr, mkl_compact_suites -> mkl_ormqr_suites) and the README/PLANS references. Verified: full MKL build (11/11 ctest) and portable -DCQR_WITH_MKL=OFF build (4/4 ctest) pass; clang-format clean. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01CbYs39BXPvFUejjDq7kYsD --- CMakeLists.txt | 18 +-- PLANS.md | 6 +- README.md | 11 +- src/cqr_compact_common.hpp | 153 ++++++++++++++++++ src/cqr_geqrf_compact.hpp | 4 +- src/{cqr_mkl_ext.cpp => cqr_mkl_ormqr.cpp} | 4 +- ...{cqr_compact.hpp => cqr_ormqr_compact.hpp} | 136 ++-------------- ...tch.cpp => cqr_ormqr_compact_dispatch.cpp} | 4 +- src/cqr_potrf_compact.hpp | 4 +- src/cqr_trsm_compact.hpp | 2 +- ...compact.cpp => test_cqr_ormqr_compact.cpp} | 6 +- ...cqr_mkl_ext.cpp => test_cqr_ormqr_mkl.cpp} | 2 +- 12 files changed, 193 insertions(+), 157 deletions(-) create mode 100644 src/cqr_compact_common.hpp rename src/{cqr_mkl_ext.cpp => cqr_mkl_ormqr.cpp} (98%) rename src/{cqr_compact.hpp => cqr_ormqr_compact.hpp} (69%) rename src/{cqr_compact_dispatch.cpp => cqr_ormqr_compact_dispatch.cpp} (97%) rename src/{test_cqr_compact.cpp => test_cqr_ormqr_compact.cpp} (98%) rename src/{test_cqr_mkl_ext.cpp => test_cqr_ormqr_mkl.cpp} (99%) diff --git a/CMakeLists.txt b/CMakeLists.txt index 30b2fbe..fa6bbfe 100644 --- a/CMakeLists.txt +++ b/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, diff --git a/PLANS.md b/PLANS.md index b225d32..182dbe5 100644 --- a/PLANS.md +++ b/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` 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 -> diff --git a/README.md b/README.md index 4cba0a6..7330f2e 100644 --- a/README.md +++ b/README.md @@ -98,17 +98,18 @@ not part of the supported interface: | File | Role | |------|------| +| `src/cqr_compact_common.hpp` | Shared machinery for the compact kernels: the `pack` 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` / `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_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_ormqr.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_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. | @@ -116,8 +117,8 @@ not part of the supported interface: | `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. | diff --git a/src/cqr_compact_common.hpp b/src/cqr_compact_common.hpp new file mode 100644 index 0000000..0712334 --- /dev/null +++ b/src/cqr_compact_common.hpp @@ -0,0 +1,153 @@ +/* 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 -- 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 +#include + +namespace cqr { +namespace detail { + +/* ------------------------------------------------------------------ */ +/* pack::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 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 +inline void vsqrt(typename pack::type &r, + const typename pack::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 +inline void broadcast(typename pack::type &v, T x) noexcept +{ + v = x - typename pack::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 struct BatchView { + VT *const data = nullptr; + const std::size_t special = 0; /* stride along the swept (reflector) axis */ + const std::size_t panel = 0; /* stride along the orthogonal panel axis */ + + VT &operator()(Int i, Int p) const noexcept + { + return data[static_cast(i) * special + + static_cast(p) * panel]; + } +}; + +/* Reinterpret a packed T buffer as a group view of V-wide pack elements. */ +template +BatchView::type, Int> +make_const_view(const T *p, std::size_t special, std::size_t panel) noexcept +{ + using VT = typename pack::type; + return {reinterpret_cast(p), special, panel}; +} +template +BatchView::type, Int> make_view(T *p, std::size_t special, + std::size_t panel) noexcept +{ + using VT = typename pack::type; + return {reinterpret_cast(p), special, panel}; +} + +} /* namespace detail */ +} /* namespace cqr */ + +#endif /* CQR_COMPACT_COMMON_HPP */ diff --git a/src/cqr_geqrf_compact.hpp b/src/cqr_geqrf_compact.hpp index 2e0f44b..8281a5a 100644 --- a/src/cqr_geqrf_compact.hpp +++ b/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 / 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, BatchView, make_view, make_const_view, vsqrt */ +#include "cqr_compact_common.hpp" /* pack, BatchView, make_view, make_const_view, vsqrt */ #include #include diff --git a/src/cqr_mkl_ext.cpp b/src/cqr_mkl_ormqr.cpp similarity index 98% rename from src/cqr_mkl_ext.cpp rename to src/cqr_mkl_ormqr.cpp index f701c67..b2bffc7 100644 --- a/src/cqr_mkl_ext.cpp +++ b/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 { diff --git a/src/cqr_compact.hpp b/src/cqr_ormqr_compact.hpp similarity index 69% rename from src/cqr_compact.hpp rename to src/cqr_ormqr_compact.hpp index 230dd36..dedb269 100644 --- a/src/cqr_compact.hpp +++ b/src/cqr_ormqr_compact.hpp @@ -1,4 +1,4 @@ -/* cqr_compact.hpp +/* cqr_ormqr_compact.hpp * * Compact (interleaved-batch) application of Householder reflectors, * templated on scalar type T and interleave width V: @@ -12,15 +12,9 @@ * The missing mkl_?ormqr_compact, in portable form. * * Design: - * - V is the compact-format interleave width (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. + * - The pack element and the BatchView addressing machinery live in the + * shared cqr_compact_common.hpp (V is the compact-format interleave width; + * see there for the width/hardware mapping). * - The algorithm is unblocked dorm2r: pivot-free and branch-free, so * the scalar code lifts verbatim with double -> V-wide vector. * - RHS columns are register-blocked (JB=4) so each reflector load @@ -72,89 +66,18 @@ * Assisted-by: Claude:claude-fable-5 Claude:claude-opus-4.8 */ -#ifndef CQR_COMPACT_HPP -#define CQR_COMPACT_HPP +#ifndef CQR_ORMQR_COMPACT_HPP +#define CQR_ORMQR_COMPACT_HPP + +#include "cqr_compact_common.hpp" /* pack, BatchView, make_view, make_const_view */ #include #include -#include #include namespace cqr { namespace detail { -/* ------------------------------------------------------------------ */ -/* pack::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 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 "ormqr_compact requires 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 -inline void vsqrt(typename pack::type &r, - const typename pack::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 -inline void broadcast(typename pack::type &v, T x) noexcept -{ - v = x - typename pack::type{}; -} - /* ------------------------------------------------------------------ */ /* Reflector sweep direction (internal control flag). */ /* */ @@ -167,47 +90,6 @@ inline void broadcast(typename pack::type &v, T x) noexcept enum class Direction { Forward, Backward }; -/* ------------------------------------------------------------------ */ -/* 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 struct BatchView { - VT *const data = nullptr; - const std::size_t special = 0; /* stride along the swept (reflector) axis */ - const std::size_t panel = 0; /* stride along the orthogonal panel axis */ - - VT &operator()(Int i, Int p) const noexcept - { - return data[static_cast(i) * special + - static_cast(p) * panel]; - } -}; - -/* Reinterpret a packed T buffer as a group view of V-wide pack elements. */ -template -BatchView::type, Int> -make_const_view(const T *p, std::size_t special, std::size_t panel) noexcept -{ - using VT = typename pack::type; - return {reinterpret_cast(p), special, panel}; -} -template -BatchView::type, Int> make_view(T *p, std::size_t special, - std::size_t panel) noexcept -{ - using VT = typename pack::type; - return {reinterpret_cast(p), special, panel}; -} - /* ------------------------------------------------------------------ */ /* One group of V interleaved matrices */ /* ------------------------------------------------------------------ */ @@ -452,4 +334,4 @@ void ormqr_compact_general(bool left, bool rowmajor, char trans, Int m, Int n, I } /* namespace detail */ } /* namespace cqr */ -#endif /* CQR_COMPACT_HPP */ +#endif /* CQR_ORMQR_COMPACT_HPP */ diff --git a/src/cqr_compact_dispatch.cpp b/src/cqr_ormqr_compact_dispatch.cpp similarity index 97% rename from src/cqr_compact_dispatch.cpp rename to src/cqr_ormqr_compact_dispatch.cpp index 62570bc..9eb25b3 100644 --- a/src/cqr_compact_dispatch.cpp +++ b/src/cqr_ormqr_compact_dispatch.cpp @@ -1,4 +1,4 @@ -/* cqr_compact_dispatch.cpp +/* cqr_ormqr_compact_dispatch.cpp * * extern "C" wrappers around the templated implementation; dispatch on * the runtime interleave width V to a compile-time instantiation. @@ -7,7 +7,7 @@ */ #include "cqr_compact.h" -#include "cqr_compact.hpp" +#include "cqr_ormqr_compact.hpp" namespace { diff --git a/src/cqr_potrf_compact.hpp b/src/cqr_potrf_compact.hpp index 8237b4c..0735a04 100644 --- a/src/cqr_potrf_compact.hpp +++ b/src/cqr_potrf_compact.hpp @@ -11,7 +11,7 @@ * potrf_compact_general() -- column- or row-major, lower or upper. * * A portable, vectorized mkl_?potrf_compact, reusing the pack / vsqrt / - * BatchView machinery from cqr_compact.hpp (shared with geqrf/ormqr). Paired with + * BatchView machinery from cqr_compact_common.hpp (shared with geqrf/ormqr). Paired with * MKL's mkl_?trsm_compact it factors and solves batched SPD systems. * * Algorithm: the unblocked LAPACK potf2, right-looking, V matrices at a time. @@ -50,7 +50,7 @@ #ifndef CQR_POTRF_COMPACT_HPP #define CQR_POTRF_COMPACT_HPP -#include "cqr_compact.hpp" /* pack, BatchView, make_view, vsqrt */ +#include "cqr_compact_common.hpp" /* pack, BatchView, make_view, vsqrt */ #include #include diff --git a/src/cqr_trsm_compact.hpp b/src/cqr_trsm_compact.hpp index f189db3..4aa3f70 100644 --- a/src/cqr_trsm_compact.hpp +++ b/src/cqr_trsm_compact.hpp @@ -33,7 +33,7 @@ #ifndef CQR_TRSM_COMPACT_HPP #define CQR_TRSM_COMPACT_HPP -#include "cqr_compact.hpp" /* pack, BatchView, make_view, make_const_view */ +#include "cqr_compact_common.hpp" /* pack, BatchView, make_view, make_const_view */ #include #include diff --git a/src/test_cqr_compact.cpp b/src/test_cqr_ormqr_compact.cpp similarity index 98% rename from src/test_cqr_compact.cpp rename to src/test_cqr_ormqr_compact.cpp index 2dc13d7..671cfff 100644 --- a/src/test_cqr_compact.cpp +++ b/src/test_cqr_ormqr_compact.cpp @@ -1,4 +1,4 @@ -/* test_cqr_compact.cpp +/* test_cqr_ormqr_compact.cpp * * Self-contained validation of the templated compact ormqr. * Reference: unblocked Householder QR (dgeqr2-style, LAPACK reflector @@ -12,7 +12,7 @@ * 2. back substitution recovers X * 3. applying 'N' after 'T' recovers the original B (Q Q^T = I) * - * Build (native): g++ -O3 -march=native -std=c++17 cqr_compact_dispatch.cpp test_cqr_compact.cpp -o test_cqr + * Build (native): g++ -O3 -march=native -std=c++17 cqr_ormqr_compact_dispatch.cpp test_cqr_ormqr_compact.cpp -o test_cqr * Build (AArch64): aarch64-linux-gnu-g++ -O3 -march=armv8.2-a -std=c++17 -static ... * * Assisted-by: Claude:claude-fable-5 Claude:claude-opus-4.8 @@ -28,7 +28,7 @@ #include #include "cqr_compact.h" -#include "cqr_compact.hpp" +#include "cqr_ormqr_compact.hpp" #include "test_compact_util.hpp" // MatrixBatch, pack/unpack, frand, max_abs_diff using namespace cqr::test; diff --git a/src/test_cqr_mkl_ext.cpp b/src/test_cqr_ormqr_mkl.cpp similarity index 99% rename from src/test_cqr_mkl_ext.cpp rename to src/test_cqr_ormqr_mkl.cpp index 0f50268..e408dac 100644 --- a/src/test_cqr_mkl_ext.cpp +++ b/src/test_cqr_ormqr_mkl.cpp @@ -1,4 +1,4 @@ -/* test_cqr_mkl_ext.cpp +/* test_cqr_ormqr_mkl.cpp * * Validation of cqr_mkl_dormqr_compact against real Intel MKL, implementing * the two test suites of cqr_mkl_dormqr_compact_design.md section 7 through From e36072494ef2c8da79b1d8d70ca23ed9454739e4 Mon Sep 17 00:00:00 2001 From: Claude Date: Thu, 6 Aug 2026 16:01:54 +0000 Subject: [PATCH 2/4] docs: correct the cqr_mkl_ormqr.cpp dispatch note in the README table The MKL adapter switches on MKL_COMPACT_PACK directly (like potrf's row already states) rather than unwrapping it to a runtime interleave width V; align the file-table description with the code. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01CbYs39BXPvFUejjDq7kYsD --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 7330f2e..3b1da2b 100644 --- a/README.md +++ b/README.md @@ -109,7 +109,7 @@ not part of the supported interface: | `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_potrf.cpp` | Dispatches on `MKL_COMPACT_PACK` directly and maps `MKL_UPLO`/`MKL_LAYOUT`, then calls the potrf kernel. | -| `src/cqr_mkl_ormqr.cpp` | Unwraps `MKL_COMPACT_PACK` -> `V` and calls the ormqr 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` | 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_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. | From edcee77391662582e762e1463d8918524685d8f4 Mon Sep 17 00:00:00 2001 From: Claude Date: Thu, 6 Aug 2026 16:02:03 +0000 Subject: [PATCH 3/4] refactor: form BatchView element offsets in Int, not size_t BatchView::operator() cast both indices to size_t before scaling by the strides, forcing the strided sweep's address arithmetic to 64-bit. The batch and matrix dimensions all fit in the kernel's Int (the compact format targets many small matrices), so type the element strides as Int and compute i*special + p*panel in Int -- no widening of the loop induction variables, keeping the strided sweep vectorizable. make_view / make_const_view take Int strides to match, and the four routine drivers drop their (size_t) element-stride casts. The large per-group base offset stays size_t and is applied to the pointer before the view is built. Mark operator() inline. No behavior change; full MKL suite (11/11) and portable -DCQR_WITH_MKL=OFF suite (4/4) pass, clang-format clean. Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01CbYs39BXPvFUejjDq7kYsD --- src/cqr_compact_common.hpp | 27 ++++++++++++++++----------- src/cqr_geqrf_compact.hpp | 4 ++-- src/cqr_ormqr_compact.hpp | 18 +++++++++--------- src/cqr_potrf_compact.hpp | 3 +-- src/cqr_trsm_compact.hpp | 8 ++++---- 5 files changed, 32 insertions(+), 28 deletions(-) diff --git a/src/cqr_compact_common.hpp b/src/cqr_compact_common.hpp index 0712334..249aef7 100644 --- a/src/cqr_compact_common.hpp +++ b/src/cqr_compact_common.hpp @@ -121,27 +121,32 @@ inline void broadcast(typename pack::type &v, T x) noexcept template struct BatchView { VT *const data = nullptr; - const std::size_t special = 0; /* stride along the swept (reflector) axis */ - const std::size_t panel = 0; /* stride along the orthogonal panel axis */ - - VT &operator()(Int i, Int p) const noexcept + 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[static_cast(i) * special + - static_cast(p) * panel]; + return data[i * special + p * panel]; } }; -/* Reinterpret a packed T buffer as a group view of V-wide pack elements. */ +/* 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 -BatchView::type, Int> -make_const_view(const T *p, std::size_t special, std::size_t panel) noexcept +BatchView::type, Int> make_const_view(const T *p, Int special, + Int panel) noexcept { using VT = typename pack::type; return {reinterpret_cast(p), special, panel}; } template -BatchView::type, Int> make_view(T *p, std::size_t special, - std::size_t panel) noexcept +BatchView::type, Int> make_view(T *p, Int special, Int panel) noexcept { using VT = typename pack::type; return {reinterpret_cast(p), special, panel}; diff --git a/src/cqr_geqrf_compact.hpp b/src/cqr_geqrf_compact.hpp index 8281a5a..da1154d 100644 --- a/src/cqr_geqrf_compact.hpp +++ b/src/cqr_geqrf_compact.hpp @@ -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(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; diff --git a/src/cqr_ormqr_compact.hpp b/src/cqr_ormqr_compact.hpp index dedb269..b250e7b 100644 --- a/src/cqr_ormqr_compact.hpp +++ b/src/cqr_ormqr_compact.hpp @@ -201,9 +201,9 @@ void ormqr_compact_group_strided(Direction dir, Int spec_len, Int panel_cnt, Int Int p = 0; /* main loop: 4 panel slices at a time; A(i,kk) loaded once, used 4x. - * Both index*stride products are evaluated in 64-bit inside operator(), - * and the panel offset p*stride is loop-invariant across i, so the - * codegen matches the hand-strided version. */ + * The index*stride products are formed in Int inside operator() (the + * batch/matrix dims fit int), and the panel offset p*stride is loop- + * invariant across i, so the codegen matches the hand-strided version. */ for (; p + 4 <= panel_cnt; p += 4) { VT w0 = C(kk, p + 0), w1 = C(kk, p + 1), w2 = C(kk, p + 2), w3 = C(kk, p + 3); for (Int i = kk + 1; i < spec_len; ++i) { @@ -302,12 +302,12 @@ void ormqr_compact_general(bool left, bool rowmajor, char trans, Int m, Int n, I * axis) with kk along its columns; for C the special axis is rows when * side='L' and columns when side='R'. Column-major: a row step is 1 and a * column step is ld; row-major flips that. */ - const std::size_t c_row = rowmajor ? (std::size_t)ldcp : 1; - const std::size_t c_col = rowmajor ? 1 : (std::size_t)ldcp; - const std::size_t a_special = rowmajor ? (std::size_t)ldap : 1; - const std::size_t a_panel = rowmajor ? 1 : (std::size_t)ldap; - const std::size_t c_special = left ? c_row : c_col; - const std::size_t c_panel = left ? c_col : c_row; + const Int c_row = rowmajor ? ldcp : 1; + const Int c_col = rowmajor ? 1 : ldcp; + const Int a_special = rowmajor ? ldap : 1; + const Int a_panel = rowmajor ? 1 : ldap; + const Int c_special = left ? c_row : c_col; + const Int c_panel = left ? c_col : c_row; /* group strides (in scalar T units): elements packed per matrix is * ldap*(complementary extent) -- the column count k for col-major (A is diff --git a/src/cqr_potrf_compact.hpp b/src/cqr_potrf_compact.hpp index 0735a04..d72f198 100644 --- a/src/cqr_potrf_compact.hpp +++ b/src/cqr_potrf_compact.hpp @@ -231,8 +231,7 @@ void potrf_compact_general(bool rowmajor, bool upper, Int n, T *ap, Int ldap, In else /* strided: sweep the reflectorless potf2 with row stride ldap, * column stride 1 -- column-major upper and row-major lower. */ - potrf_compact_group_strided( - n, make_view(a, static_cast(ldap), 1)); + potrf_compact_group_strided(n, make_view(a, ldap, 1)); } } diff --git a/src/cqr_trsm_compact.hpp b/src/cqr_trsm_compact.hpp index 4aa3f70..12b3758 100644 --- a/src/cqr_trsm_compact.hpp +++ b/src/cqr_trsm_compact.hpp @@ -243,10 +243,10 @@ void trsm_compact_general(bool left, bool upper, bool rowmajor, bool tran, bool /* element strides (in VT units). Column-major: a row step is 1 and a column * step is ld; row-major flips that. A is s x s, B is m x n. */ - const std::size_t a_row = rowmajor ? (std::size_t)ldap : 1; - const std::size_t a_col = rowmajor ? 1 : (std::size_t)ldap; - const std::size_t b_row = rowmajor ? (std::size_t)ldbp : 1; - const std::size_t b_col = rowmajor ? 1 : (std::size_t)ldbp; + const Int a_row = rowmajor ? ldap : 1; + const Int a_col = rowmajor ? 1 : ldap; + const Int b_row = rowmajor ? ldbp : 1; + const Int b_col = rowmajor ? 1 : ldbp; /* group strides (in scalar T units): elements packed per matrix is * ld*(complementary extent) -- for A (s x s) that is ldap*s either way; From d380ba2238255b599361d8ee1c11732174e36b25 Mon Sep 17 00:00:00 2001 From: Claude Date: Thu, 6 Aug 2026 16:50:34 +0000 Subject: [PATCH 4/4] docs: align the geqrf/trsm dispatch notes with ormqr/potrf Both MKL adapters switch on MKL_COMPACT_PACK directly, like potrf and ormqr; drop the older "unwraps -> V" wording so all four cqr_mkl_*.cpp file-table rows describe the dispatch the same way. Each row still names the specific modes its adapter maps (geqrf: MKL_LAYOUT; trsm: the full MKL enum set). Co-Authored-By: Claude Opus 4.8 Claude-Session: https://claude.ai/code/session_01CbYs39BXPvFUejjDq7kYsD --- README.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/README.md b/README.md index 3b1da2b..7d82758 100644 --- a/README.md +++ b/README.md @@ -107,10 +107,10 @@ not part of the supported interface: | `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_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_ormqr.cpp` | Dispatches on `MKL_COMPACT_PACK` directly and maps `side`/`trans`/`MKL_LAYOUT`, then 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_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). |