From c6ccae484d4bad2aca998024b7e411d6475964cf Mon Sep 17 00:00:00 2001 From: Claude Date: Mon, 14 Sep 2026 22:07:27 +0000 Subject: [PATCH 1/8] Validate the test suites against a real LAPACKE + CBLAS stack The hand-rolled scalar references of the portable suites (ref_geqr2, ref_orm2r, ref_org2r, ref_geqp3, ref_potf2, ref_trsm, ref_gels, matmul, tri_apply) now forward to LAPACKE_?geqr2 / ?ormqr / ?orgqr / ?geqp3 / ?potrf / ?gels and cblas_?trsm / ?trmm / ?gemm through the new tests/test_lapack_util.hpp, which hands a MatrixView's data and leading dimension to the library and stages the few LAPACKE operands that have no transpose flag. The unpivoted LDL^T keeps its scalar reference (LAPACK has none) and validates it against cblas_?trmm first. Where the reference is LAPACK's blocked driver the gates are relative to the operand norms at a multiple of n eps (the design documents' 7.1 gates); only ?geqr2 stays an elementwise ~eps comparison. cmake/FindLAPACKE.cmake finds the stack (MKL's own with the MKL extension, otherwise OpenBLAS, Netlib, or any lapacke + find_package(LAPACK) pair; -DCBK_TEST_LAPACK pins one, which matters under Debian's alternatives symlinks) and exposes it as LAPACKE::LAPACKE. The dense-LAPACK invariants of the geqrf and potrf MKL suites moved into the portable suites, so they run on every stack; the MKL suites keep the compat API coverage and the cross-checks against MKL's compact kernels. CI runs OpenBLAS and Netlib under gcc and clang next to the MKL jobs; a tag-only macOS workflow targets Accelerate through accelerate-lapacke (unverified). Closes #27. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_01GdNaquWMN3b39gvTsXZyjn --- .claude/CLAUDE.md | 35 ++- .claude/hooks/session-start.sh | 9 + .github/workflows/ci.yml | 35 ++- .github/workflows/macos.yml | 40 +++ CHANGELOG.md | 12 + CMakeLists.txt | 45 ++- README.md | 4 +- cmake/FindLAPACKE.cmake | 252 ++++++++++++++++ docs/PLANS.md | 71 +++-- docs/building.md | 34 ++- docs/cbk_dgels_compact_design.md | 9 +- docs/cbk_dgeqrf_compact_design.md | 21 +- docs/cbk_dorgqr_compact_design.md | 15 +- docs/cbk_dpotrf_compact_design.md | 24 +- docs/cbk_dsytrfnp_compact_design.md | 11 +- docs/cbk_dtrsm_compact_design.md | 13 +- docs/layout.md | 17 +- tests/test_compact_util.hpp | 278 +----------------- tests/test_gels_compact.cpp | 47 +-- tests/test_gels_mkl.cpp | 2 +- tests/test_geqrf_compact.cpp | 161 ++++++++++- tests/test_geqrf_mkl.cpp | 177 +----------- tests/test_lapack_util.hpp | 427 ++++++++++++++++++++++++++++ tests/test_mkl_util.hpp | 34 +-- tests/test_orgqr_compact.cpp | 53 +--- tests/test_orgqr_mkl.cpp | 2 +- tests/test_ormqr_compact.cpp | 75 ++--- tests/test_ormqr_mkl.cpp | 2 +- tests/test_potrf_compact.cpp | 137 ++++----- tests/test_potrf_mkl.cpp | 124 +------- tests/test_sytrfnp_compact.cpp | 78 ++++- tests/test_sytrfnp_mkl.cpp | 2 +- tests/test_trsm_compact.cpp | 29 +- tests/test_trsm_mkl.cpp | 2 +- 34 files changed, 1346 insertions(+), 931 deletions(-) create mode 100644 .github/workflows/macos.yml create mode 100644 cmake/FindLAPACKE.cmake create mode 100644 tests/test_lapack_util.hpp diff --git a/.claude/CLAUDE.md b/.claude/CLAUDE.md index 1a8b377..e5d0a4f 100644 --- a/.claude/CLAUDE.md +++ b/.claude/CLAUDE.md @@ -8,8 +8,13 @@ optional Intel MKL-style API that drops into MKL's compact ecosystem. ## Build and test -The default configure builds the portable library and its BLAS-free tests and -needs nothing but CMake and a C++17 compiler: +The default configure builds the portable library and its tests. The library +needs nothing but CMake and a C++17 compiler; the tests validate against a +real LAPACKE + CBLAS (`cmake/FindLAPACKE.cmake`: OpenBLAS or Netlib without +MKL, `libopenblas-dev liblapacke-dev` / `liblapacke-dev liblapack-dev` on +Debian, MKL's own with the extension; `-DCBK_TEST_LAPACK=openblas|netlib|mkl` +pins one -- do that, since Debian's alternatives symlinks let either stand in +for the other): ```sh cmake -S . -B build -DCMAKE_BUILD_TYPE=Release @@ -18,8 +23,8 @@ ctest --test-dir build --output-on-failure ``` `-DCBK_WITH_MKL=ON` adds the MKL-style API, the MKL-backed suites (validated -against MKL's compact kernels and LAPACK/LAPACKE), the example and the -benchmarks; it needs Intel MKL (`sudo apt-get install libmkl-dev` on +against MKL's compact kernels; every suite's LAPACKE is then MKL's), the +example and the benchmarks; it needs Intel MKL (`sudo apt-get install libmkl-dev` on Debian/Ubuntu). Run that build too before pushing changes to the kernels or the wrappers: it is the one that cross-checks against MKL. @@ -61,9 +66,12 @@ src/ the templated kernels (cbk_*_compact.hpp, one per routine, on the and the owning MatrixBatch the tests, benchmarks and examples share (internal: src/ is on their include path, but the public API stays include/) -tests/ portable (no BLAS) and MKL-backed suites, templated on the scalar - type; test_compact_util.hpp / test_mkl_util.hpp hold the helpers and - the compact / compat / mkl / lapack dispatch structs +tests/ the portable C API's suites (vs LAPACKE/CBLAS on any stack) and + the MKL-backed ones, templated on the scalar type; + test_compact_util.hpp (compact, the generators, pack/unpack), + test_lapack_util.hpp (lapack and the ref_* procedures over + views, matmul, tri_apply, solve_errors) and test_mkl_util.hpp + (compat, mkl) hold the helpers and dispatch structs (tests/install/: the install check's consumer project and script) examples/ the worked solve and the benchmarks (BENCHMARKS.md), on bench_util.hpp docs/ one design document per routine, the guides README.md indexes @@ -175,10 +183,12 @@ workspace contract, or the benchmarks' threading. - **Build and test with both gcc and clang before pushing.** CI runs both, and the packs' alignment is exactly the kind of contract only one of them enforces: both alignment faults so far (issue #34 and the one above) were - invisible to gcc. The default (portable) tree needs nothing but the compiler: + invisible to gcc. The default (portable) tree needs only the compiler and a + LAPACKE stack for the tests (CI runs OpenBLAS and Netlib under both + compilers): ```sh - CXX=clang++ cmake -S . -B build-clang && cmake --build build-clang && ctest --test-dir build-clang + CXX=clang++ cmake -S . -B build-clang -DCBK_TEST_LAPACK=openblas && cmake --build build-clang && ctest --test-dir build-clang ``` - **Threading over groups.** Every all-groups driver is a call to @@ -201,9 +211,14 @@ workspace contract, or the benchmarks' threading. library and its tests are one internal codebase and share these views on purpose; what keeps the suites honest is that they compute the *answers* independently -- scalar LAPACK references, dense LAPACK/MKL cross-checks. + The dense references are the library's: `test_lapack_util.hpp` hands a + view's `data` and leading dimension to LAPACKE/CBLAS (`is_rowmajor`, + `lapack_ld`, `lapack_layout`), reads a BLAS-3 operand stored the other way + round as its transpose, and stages the few LAPACKE operands that have no + transpose flag (`Staged`). A routine that takes views (or batches) `assert`s what its contract assumes: dimension compatibility across operands (`matmul`, `tri_apply`, - `solve_errors`, the scalar references), squareness where required + `solve_errors`, the `ref_*` procedures), squareness where required (`gen_spd`, `gen_tri`, `ref_potf2`, ...), and index ranges (`MatrixView::operator()`, `MatrixBatch::operator[]`). Costs nothing in Release, and the Debug suite run is what exercises it. `BatchView` is the diff --git a/.claude/hooks/session-start.sh b/.claude/hooks/session-start.sh index afe62d5..f14b17d 100755 --- a/.claude/hooks/session-start.sh +++ b/.claude/hooks/session-start.sh @@ -8,6 +8,11 @@ # libmkl-dev Intel MKL (see the Prerequisites section of README.md): # the distro package, since Intel's own apt repo # (apt.repos.intel.com) is blocked by the web network policy. +# libopenblas-dev, liblapacke-dev, liblapack-dev +# the two non-MKL LAPACKE stacks the test suites validate +# against (cmake/FindLAPACKE.cmake), so a session can run +# the tests with -DCBK_TEST_LAPACK=openblas and =netlib as +# CI does. # libomp-dev clang's OpenMP runtime and its omp.h, which clang-tidy # needs (.pre-commit-config.yaml explains). # pre-commit the formatting and linting driver, with its hook @@ -38,6 +43,9 @@ fi # without them ever shows up. missing=() [ -f /usr/include/mkl/mkl_compact.h ] || missing+=(libmkl-dev) +dpkg -s libopenblas-dev >/dev/null 2>&1 || missing+=(libopenblas-dev) +dpkg -s liblapacke-dev >/dev/null 2>&1 || missing+=(liblapacke-dev) +dpkg -s liblapack-dev >/dev/null 2>&1 || missing+=(liblapack-dev) dpkg -s libomp-dev >/dev/null 2>&1 || missing+=(libomp-dev) command -v clang-format >/dev/null 2>&1 || missing+=(clang-format) command -v clang-tidy >/dev/null 2>&1 || missing+=(clang-tidy) @@ -68,5 +76,6 @@ else fi echo "session-start: build with 'cmake -S . -B build -DCBK_WITH_MKL=ON && cmake --build build && ctest --test-dir build' (MKL is off by default)" +echo "session-start: the tests validate against LAPACKE; without MKL pick the stack with -DCBK_TEST_LAPACK=openblas|netlib (default: the first found)" echo "session-start: 'pre-commit run --all-files' checks the style; see .claude/CLAUDE.md" echo "session-start: for clang-tidy, 'CXX=clang++ cmake -S . -B build-tidy' then 'pre-commit run --hook-stage manual clang-tidy --all-files'" diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 64088e0..39c883d 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -12,21 +12,30 @@ concurrency: jobs: build-and-test: - name: ${{ matrix.cxx }} / MKL ${{ matrix.mkl }}${{ matrix.interface == 'ilp64' && ' (ilp64)' || '' }} + name: ${{ matrix.cxx }} / ${{ matrix.mkl == 'ON' && 'MKL' || matrix.lapack }}${{ matrix.interface == 'ilp64' && ' (ilp64)' || '' }} runs-on: ubuntu-latest strategy: fail-fast: false matrix: cc: [gcc, clang] + # The tests validate against a real LAPACKE + CBLAS: MKL's own with + # the MKL extension, OpenBLAS or Netlib's reference LAPACK without it + # (cmake/FindLAPACKE.cmake; -DCBK_TEST_LAPACK pins the stack, so a + # job cannot silently pick up the other one through Debian's + # update-alternatives symlinks). mkl: ["ON", "OFF"] + lapack: [openblas, netlib] interface: [lp64] + exclude: + - { mkl: "ON", lapack: netlib } include: - { cc: gcc, cxx: g++ } - { cc: clang, cxx: clang++ } # MKL's 64-bit integer interface: MKL_INT is part of the MKL-style # API's ABI, so the wrappers and their suites must build and pass - # against it too (libmkl-dev ships both interface layers). - - { cc: gcc, cxx: g++, mkl: "ON", interface: ilp64 } + # against it too (libmkl-dev ships both interface layers). The + # LAPACKE the suites call is MKL's, at the same interface. + - { cc: gcc, cxx: g++, mkl: "ON", lapack: openblas, interface: ilp64 } env: CC: ${{ matrix.cc }} CXX: ${{ matrix.cxx }} @@ -34,19 +43,35 @@ jobs: - uses: actions/checkout@v4 # libmkl-dev supplies the Compact extension and MKL's LAPACK/LAPACKE, so - # the MKL-backed suites validate against real MKL + LAPACK. The portable - # (MKL OFF) build needs nothing beyond the C++17 toolchain the image ships. + # the MKL-backed suites validate against real MKL + LAPACK, and the + # portable suites against MKL's LAPACKE. - name: Install Intel MKL if: matrix.mkl == 'ON' run: | sudo apt-get update sudo apt-get install -y --no-install-recommends libmkl-dev + # Without MKL the suites need a LAPACKE stack: OpenBLAS (Debian ships + # its LAPACKE as Netlib's liblapacke, hence both packages), or Netlib's + # reference LAPACK/BLAS with its LAPACKE. + - name: Install OpenBLAS + if: matrix.mkl == 'OFF' && matrix.lapack == 'openblas' + run: | + sudo apt-get update + sudo apt-get install -y --no-install-recommends libopenblas-dev liblapacke-dev + + - name: Install Netlib LAPACK + if: matrix.mkl == 'OFF' && matrix.lapack == 'netlib' + run: | + sudo apt-get update + sudo apt-get install -y --no-install-recommends liblapacke-dev liblapack-dev libblas-dev + - name: Configure run: > cmake -S . -B build -G Ninja -DCMAKE_BUILD_TYPE=Release -DCBK_WITH_MKL=${{ matrix.mkl }} + -DCBK_TEST_LAPACK=${{ matrix.mkl == 'ON' && 'mkl' || matrix.lapack }} -DMKLCompact_INTERFACE=${{ matrix.interface }} - name: Build diff --git a/.github/workflows/macos.yml b/.github/workflows/macos.yml new file mode 100644 index 0000000..a42c1d5 --- /dev/null +++ b/.github/workflows/macos.yml @@ -0,0 +1,40 @@ +name: macOS + +# The test suites on macOS, against Apple's Accelerate through +# accelerate-lapacke (https://github.com/lepus2589/accelerate-lapacke), which +# builds the LAPACKE interface Accelerate lacks. macOS runners are expensive, +# so this runs only for release tags and on manual dispatch, never per push +# or pull request. Unverified as of its introduction (no macOS machine was at +# hand): expect to adjust the accelerate-lapacke preset or the CBLAS header +# mapping (CBK_LAPACKE_ACCELERATE in tests/test_lapack_util.hpp) on the first +# run. +on: + push: + tags: ['v*'] + workflow_dispatch: + +jobs: + accelerate: + name: AppleClang / Accelerate + runs-on: macos-latest + steps: + - uses: actions/checkout@v4 + + - name: Build accelerate-lapacke + run: | + git clone --depth 1 https://github.com/lepus2589/accelerate-lapacke.git "$RUNNER_TEMP/accelerate-lapacke" + cd "$RUNNER_TEMP/accelerate-lapacke" + cmake --workflow --preset accelerate-lapacke32-install + + - name: Configure + run: > + cmake -S . -B build + -DCMAKE_BUILD_TYPE=Release + -DCBK_TEST_LAPACK=accelerate + -DCMAKE_PREFIX_PATH="$HOME/.local" + + - name: Build + run: cmake --build build -j + + - name: Test + run: ctest --test-dir build --output-on-failure diff --git a/CHANGELOG.md b/CHANGELOG.md index a9e3ee4..bde4ee3 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,6 +5,18 @@ All notable changes to cbk are recorded here. The format follows numbers [Semantic Versioning](https://semver.org/): until 1.0.0, a minor release may change the API. +## [Unreleased] + +### Changed + +- Tests: every suite validates against a real LAPACKE + CBLAS stack + (OpenBLAS, Netlib, or MKL's own with the MKL extension; + `cmake/FindLAPACKE.cmake`, `-DCBK_TEST_LAPACK`) instead of hand-rolled + scalar references (issue #27). The dense-LAPACK cross-checks of the MKL + suites run in the portable suites, on any stack; the unpivoted LDL^T keeps + its scalar reference, itself validated against CBLAS. Building the tests + now needs such a stack; the library does not. + ## [0.1.0] - 2026-09-13 The first release: batched QR, Cholesky and unpivoted LDL^T factorizations, diff --git a/CMakeLists.txt b/CMakeLists.txt index c25c2a7..2ecff85 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -16,6 +16,13 @@ set(CMAKE_EXPORT_COMPILE_COMMANDS ON) option(CBK_WITH_MKL "Build the Intel MKL compact extension, examples and MKL tests" OFF) option(CBK_WITH_OPENMP "Thread each routine's loop over compact groups with OpenMP" ON) option(CBK_BUILD_TESTS "Build tests and register them with CTest" ON) +# The tests validate against a real LAPACKE + CBLAS (the library itself needs +# none): which stack, as cmake/FindLAPACKE.cmake's LAPACKE_VENDOR. `auto` +# takes MKL's with the MKL extension and otherwise the first of OpenBLAS, +# Netlib, any `lapacke` + find_package(LAPACK) pair; the other values pin one. +set(CBK_TEST_LAPACK "auto" CACHE STRING + "LAPACKE/CBLAS stack the tests validate against: auto | mkl | openblas | netlib | accelerate") +set_property(CACHE CBK_TEST_LAPACK PROPERTY STRINGS auto mkl openblas netlib accelerate) # The install rules belong to cbk's own build, not to a project that # add_subdirectory()s it. option(CBK_INSTALL "Generate the install target" ${PROJECT_IS_TOP_LEVEL}) @@ -42,7 +49,7 @@ endif() # --------------------------------------------------------------------------- # The library: one target, cbk (cbk::cbk). By default the portable C API of # include/cbk.h (no BLAS, no MKL), which is what a default configure -# builds and tests. Static by default; -DBUILD_SHARED_LIBS=ON builds it shared. +# builds and tests (the tests, not the library, need a LAPACKE stack). Static by default; -DBUILD_SHARED_LIBS=ON builds it shared. # --------------------------------------------------------------------------- add_library(cbk src/cbk.cpp) target_include_directories(cbk @@ -109,23 +116,47 @@ endif() if(CBK_BUILD_TESTS) enable_testing() - # Portable self-contained suites (no BLAS): each kernel vs a scalar reference - # of the same unblocked algorithm, plus the C API's LAPACK-style validation. + # The dense references of every suite come from a real LAPACKE + CBLAS + # (tests/test_lapack_util.hpp): with the MKL extension it must be MKL's own + # (one LAPACK stack per binary), otherwise any of OpenBLAS, Netlib, ... + if(CBK_WITH_MKL) + if(NOT CBK_TEST_LAPACK MATCHES "^(auto|mkl)$") + message(FATAL_ERROR "CBK_TEST_LAPACK=${CBK_TEST_LAPACK} cannot be combined with CBK_WITH_MKL=ON: the MKL-backed suites link MKL, so the dense references must be MKL's LAPACKE too") + endif() + set(LAPACKE_VENDOR MKL) + elseif(CBK_TEST_LAPACK STREQUAL "mkl") + set(LAPACKE_VENDOR MKL) + elseif(CBK_TEST_LAPACK STREQUAL "openblas") + set(LAPACKE_VENDOR OpenBLAS) + elseif(CBK_TEST_LAPACK STREQUAL "netlib") + set(LAPACKE_VENDOR Netlib) + elseif(CBK_TEST_LAPACK STREQUAL "accelerate") + set(LAPACKE_VENDOR Accelerate) + elseif(CBK_TEST_LAPACK STREQUAL "auto") + set(LAPACKE_VENDOR Any) + else() + message(FATAL_ERROR "CBK_TEST_LAPACK must be auto, mkl, openblas, netlib or accelerate (got '${CBK_TEST_LAPACK}')") + endif() + find_package(LAPACKE REQUIRED) + + # The portable C API's suites: each kernel vs LAPACKE/CBLAS (its unblocked + # LAPACK counterpart elementwise where that is the same algorithm, residual + # and invariant gates elsewhere), plus the C API's LAPACK-style validation. foreach(routine geqrf ormqr orgqr potrf sytrfnp trsm gels) add_executable(test_${routine}_compact tests/test_${routine}_compact.cpp) - target_link_libraries(test_${routine}_compact PRIVATE cbk) + target_link_libraries(test_${routine}_compact PRIVATE cbk LAPACKE::LAPACKE) target_include_directories(test_${routine}_compact PRIVATE "${CMAKE_CURRENT_SOURCE_DIR}/src") add_test(NAME portable_${routine} COMMAND test_${routine}_compact) endforeach() if(CBK_WITH_MKL) - # MKL + dense-LAPACK validation of each cbk_?*_compact routine (the - # design docs' section 7): invariants vs LAPACK, cross-checks vs MKL's own + # The MKL-style API's suites (the design docs' section 7): the compat + # entry points through MKL's own pack/unpack, cross-checked against MKL's # compact kernels, and end-to-end solves. foreach(routine geqrf ormqr orgqr potrf sytrfnp trsm gels) add_executable(test_${routine}_mkl tests/test_${routine}_mkl.cpp) - target_link_libraries(test_${routine}_mkl PRIVATE cbk MKL::Compact) + target_link_libraries(test_${routine}_mkl PRIVATE cbk MKL::Compact LAPACKE::LAPACKE) target_include_directories(test_${routine}_mkl PRIVATE "${CMAKE_CURRENT_SOURCE_DIR}/src") add_test(NAME mkl_${routine}_suites COMMAND test_${routine}_mkl) diff --git a/README.md b/README.md index 3b0dc2d..4fe84c4 100644 --- a/README.md +++ b/README.md @@ -34,7 +34,9 @@ groups of `V` matrices with OpenMP (`-DCBK_WITH_OPENMP=OFF` disables it). ## Getting started -Requires CMake >= 3.28 and a C++17 compiler (GCC/Clang); nothing else. +Requires CMake >= 3.28 and a C++17 compiler (GCC/Clang); the library needs +nothing else, and the tests a LAPACKE + CBLAS stack (OpenBLAS, Netlib or MKL) +to validate against. ```sh cmake -S . -B build -DCMAKE_BUILD_TYPE=Release -DCMAKE_CXX_FLAGS="-O3 -march=native" diff --git a/cmake/FindLAPACKE.cmake b/cmake/FindLAPACKE.cmake new file mode 100644 index 0000000..3befa75 --- /dev/null +++ b/cmake/FindLAPACKE.cmake @@ -0,0 +1,252 @@ +# FindLAPACKE.cmake +# +# Locate a LAPACKE + CBLAS stack for the test suites and expose it as +# +# LAPACKE::LAPACKE -- the include dirs, the link line, and one compile +# definition naming the stack's headers: +# CBK_LAPACKE_MKL , +# CBK_LAPACKE_ACCELERATE , +# (none) , +# +# LAPACKE_VENDOR Any (default) | MKL | OpenBLAS | Netlib | Accelerate +# which stack to look for; Any tries OpenBLAS, then Netlib, +# then whatever `lapacke` library and find_package(LAPACK) +# turn up. +# LAPACKE_VENDOR_FOUND the stack that was found, as one of the names above. +# +# Every candidate is link-tested with a program that calls LAPACKE_dgeqrf, +# LAPACKE_dgeqr2, cblas_dgemm and cblas_dtrsm, so a header without its +# library (or the other way round) is rejected rather than found. +# +# MKL reuses FindMKLCompact.cmake's MKL::Compact (the whole MKL link +# line, which carries LAPACKE and CBLAS). With -DCBK_WITH_MKL=ON +# this is the only acceptable stack: a test binary must never +# link two LAPACK implementations. +# OpenBLAS pkg-config `openblas` (the distro's include subdirectory -- +# Debian's openblas-pthread/, Fedora's openblas/ -- comes from +# there), else cblas.h + libopenblas by hand. Debian and Ubuntu +# strip LAPACKE out of libopenblas and ship Netlib's +# liblapacke separately, so when the link test fails on +# libopenblas alone it is retried with liblapacke added. +# Netlib liblapacke + liblapack + libblas. On Debian the plain +# liblapack.so / libblas.so are update-alternatives symlinks +# that usually point at OpenBLAS, so the reference libraries +# are looked for under the lapack/ and blas/ subdirectories +# first (PATH_SUFFIXES), which is where the distro keeps them. +# Accelerate macOS: Apple's Accelerate framework provides LAPACK and +# CBLAS but no LAPACKE; https://github.com/lepus2589/accelerate-lapacke +# builds one on top of it and installs a `lapacke` CMake +# package (find_package(LAPACKE CONFIG)). Unverified in CI: +# the macOS workflow runs on tags and by hand only. +# +# ILP64 is supported through MKL only (MKL_ILP64 comes with MKL::Compact); +# the other stacks are used at their default 32-bit integer interface. +# +# Assisted-by: Claude:claude-fable-5 + +set(LAPACKE_VENDOR "Any" CACHE STRING + "LAPACKE/CBLAS stack for the tests: Any | MKL | OpenBLAS | Netlib | Accelerate") +set_property(CACHE LAPACKE_VENDOR PROPERTY STRINGS Any MKL OpenBLAS Netlib Accelerate) +if(NOT LAPACKE_VENDOR MATCHES "^(Any|MKL|OpenBLAS|Netlib|Accelerate)$") + message(FATAL_ERROR "LAPACKE_VENDOR must be Any, MKL, OpenBLAS, Netlib or Accelerate (got '${LAPACKE_VENDOR}')") +endif() + +include(CMakePushCheckState) +include(CheckCXXSourceCompiles) + +# The link test. `defs` selects the header set exactly as the test util does. +function(_lapacke_link_test vendor includes libs defs outvar) + cmake_push_check_state(RESET) + set(CMAKE_REQUIRED_INCLUDES ${includes}) + set(CMAKE_REQUIRED_LIBRARIES ${libs}) + set(CMAKE_REQUIRED_DEFINITIONS ${defs}) + set(CMAKE_REQUIRED_QUIET TRUE) + check_cxx_source_compiles(" + #if defined(CBK_LAPACKE_MKL) + #include + #include + #elif defined(CBK_LAPACKE_ACCELERATE) + #include + #include + #else + #include + #include + #endif + int main() { + double a[4] = {2, 0, 0, 2}, tau[2], c[4] = {1, 0, 0, 1}; + lapack_int info = LAPACKE_dgeqrf(LAPACK_COL_MAJOR, 2, 2, a, 2, tau); + info += LAPACKE_dgeqr2(LAPACK_COL_MAJOR, 2, 2, a, 2, tau); + cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, 2, 2, 2, 1.0, a, 2, c, 2, 0.0, c, 2); + cblas_dtrsm(CblasColMajor, CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, 2, 2, 1.0, a, 2, c, 2); + return (int)info; + } + " LAPACKE_LINKS_${vendor}) + cmake_pop_check_state() + set(${outvar} ${LAPACKE_LINKS_${vendor}} PARENT_SCOPE) +endfunction() + +set(_lapacke_found_vendor "") +set(_lapacke_includes "") +set(_lapacke_libs "") +set(_lapacke_defs "") + +# ---- MKL ------------------------------------------------------------------- +if(NOT _lapacke_found_vendor AND LAPACKE_VENDOR STREQUAL "MKL") + find_package(MKLCompact REQUIRED) + _lapacke_link_test(MKL "${MKLCompact_INCLUDE_DIR}" "MKL::Compact" "-DCBK_LAPACKE_MKL" _ok) + if(_ok) + set(_lapacke_found_vendor MKL) + set(_lapacke_includes "") # MKL::Compact carries them + set(_lapacke_libs MKL::Compact) + set(_lapacke_defs CBK_LAPACKE_MKL) + endif() +endif() + +# ---- OpenBLAS --------------------------------------------------------------- +if(NOT _lapacke_found_vendor AND LAPACKE_VENDOR MATCHES "^(Any|OpenBLAS)$") + find_package(PkgConfig QUIET) + set(_ob_includes "") + set(_ob_libs "") + if(PKG_CONFIG_FOUND) + pkg_check_modules(_ob QUIET openblas) + if(_ob_FOUND) + set(_ob_includes ${_ob_INCLUDE_DIRS}) + set(_ob_libs ${_ob_LINK_LIBRARIES}) + endif() + endif() + if(NOT _ob_libs) + find_library(LAPACKE_OpenBLAS_LIBRARY NAMES openblas + PATH_SUFFIXES openblas openblas-pthread openblas-openmp openblas-serial) + find_path(LAPACKE_OpenBLAS_INCLUDE_DIR NAMES openblas_config.h + PATH_SUFFIXES openblas openblas-pthread openblas-openmp openblas-serial) + mark_as_advanced(LAPACKE_OpenBLAS_LIBRARY LAPACKE_OpenBLAS_INCLUDE_DIR) + if(LAPACKE_OpenBLAS_LIBRARY AND LAPACKE_OpenBLAS_INCLUDE_DIR) + set(_ob_libs "${LAPACKE_OpenBLAS_LIBRARY}") + set(_ob_includes "${LAPACKE_OpenBLAS_INCLUDE_DIR}") + endif() + endif() + if(_ob_libs) + # lapacke.h sits next to cblas.h in a full OpenBLAS install, or in the + # default include dir when the distro ships Netlib's LAPACKE separately. + find_path(LAPACKE_OpenBLAS_LAPACKE_INCLUDE_DIR NAMES lapacke.h + HINTS ${_ob_includes} + PATH_SUFFIXES openblas openblas-pthread openblas-openmp openblas-serial) + mark_as_advanced(LAPACKE_OpenBLAS_LAPACKE_INCLUDE_DIR) + if(LAPACKE_OpenBLAS_LAPACKE_INCLUDE_DIR) + list(APPEND _ob_includes "${LAPACKE_OpenBLAS_LAPACKE_INCLUDE_DIR}") + list(REMOVE_DUPLICATES _ob_includes) + _lapacke_link_test(OpenBLAS "${_ob_includes}" "${_ob_libs}" "" _ok) + if(NOT _ok) + # Debian/Ubuntu: LAPACKE lives in Netlib's liblapacke, which calls + # back into whichever LAPACK is loaded -- libopenblas here, since it + # is linked first and exports the LAPACK symbols itself. + find_library(LAPACKE_OpenBLAS_LAPACKE_LIBRARY NAMES lapacke) + mark_as_advanced(LAPACKE_OpenBLAS_LAPACKE_LIBRARY) + if(LAPACKE_OpenBLAS_LAPACKE_LIBRARY) + list(PREPEND _ob_libs "${LAPACKE_OpenBLAS_LAPACKE_LIBRARY}") + _lapacke_link_test(OpenBLAS_lapacke "${_ob_includes}" "${_ob_libs}" "" _ok) + endif() + endif() + if(_ok) + set(_lapacke_found_vendor OpenBLAS) + set(_lapacke_includes ${_ob_includes}) + set(_lapacke_libs ${_ob_libs}) + endif() + endif() + endif() +endif() + +# ---- Netlib ----------------------------------------------------------------- +if(NOT _lapacke_found_vendor AND LAPACKE_VENDOR MATCHES "^(Any|Netlib)$") + find_path(LAPACKE_Netlib_INCLUDE_DIR NAMES lapacke.h PATH_SUFFIXES lapacke) + find_path(LAPACKE_Netlib_CBLAS_INCLUDE_DIR NAMES cblas.h PATH_SUFFIXES cblas) + find_library(LAPACKE_Netlib_LAPACKE_LIBRARY NAMES lapacke) + # The reference libraries in their own subdirectories first (Debian), then + # the plain names. + find_library(LAPACKE_Netlib_LAPACK_LIBRARY NAMES lapack PATH_SUFFIXES lapack) + find_library(LAPACKE_Netlib_BLAS_LIBRARY NAMES blas PATH_SUFFIXES blas) + mark_as_advanced(LAPACKE_Netlib_INCLUDE_DIR LAPACKE_Netlib_CBLAS_INCLUDE_DIR + LAPACKE_Netlib_LAPACKE_LIBRARY LAPACKE_Netlib_LAPACK_LIBRARY LAPACKE_Netlib_BLAS_LIBRARY) + if(LAPACKE_Netlib_INCLUDE_DIR AND LAPACKE_Netlib_CBLAS_INCLUDE_DIR + AND LAPACKE_Netlib_LAPACKE_LIBRARY AND LAPACKE_Netlib_LAPACK_LIBRARY + AND LAPACKE_Netlib_BLAS_LIBRARY) + set(_nl_includes "${LAPACKE_Netlib_INCLUDE_DIR}" "${LAPACKE_Netlib_CBLAS_INCLUDE_DIR}") + list(REMOVE_DUPLICATES _nl_includes) + set(_nl_libs "${LAPACKE_Netlib_LAPACKE_LIBRARY}" "${LAPACKE_Netlib_LAPACK_LIBRARY}" + "${LAPACKE_Netlib_BLAS_LIBRARY}") + # A test binary calls LAPACKE and CBLAS, never LAPACK itself, so GNU ld's + # --as-needed (Debian's default) would drop liblapack from its NEEDED + # list and liblapacke would then load `liblapack.so.3` through the + # loader cache -- the alternatives symlink, i.e. OpenBLAS again. Keeping + # it a direct dependency makes the loader take the reference library + # from the RUNPATH CMake sets for the build tree, and liblapacke's own + # request for the same soname is then satisfied by the copy already + # loaded. + if(CMAKE_SYSTEM_NAME STREQUAL "Linux" AND NOT MSVC) + list(PREPEND _nl_libs "-Wl,--no-as-needed") + endif() + _lapacke_link_test(Netlib "${_nl_includes}" "${_nl_libs}" "" _ok) + if(_ok) + set(_lapacke_found_vendor Netlib) + set(_lapacke_includes ${_nl_includes}) + set(_lapacke_libs ${_nl_libs}) + endif() + endif() +endif() + +# ---- Accelerate (macOS) ----------------------------------------------------- +if(NOT _lapacke_found_vendor AND LAPACKE_VENDOR STREQUAL "Accelerate") + find_package(LAPACKE CONFIG QUIET) # accelerate-lapacke's package: target `lapacke` + find_library(LAPACKE_Accelerate_FRAMEWORK Accelerate) + mark_as_advanced(LAPACKE_Accelerate_FRAMEWORK) + if(TARGET lapacke AND LAPACKE_Accelerate_FRAMEWORK) + set(_ac_libs lapacke "${LAPACKE_Accelerate_FRAMEWORK}") + _lapacke_link_test(Accelerate "" "${_ac_libs}" "-DCBK_LAPACKE_ACCELERATE" _ok) + if(_ok) + set(_lapacke_found_vendor Accelerate) + set(_lapacke_libs ${_ac_libs}) + set(_lapacke_defs CBK_LAPACKE_ACCELERATE) + endif() + endif() +endif() + +# ---- Any other lapacke + LAPACK/BLAS pair ----------------------------------- +if(NOT _lapacke_found_vendor AND LAPACKE_VENDOR STREQUAL "Any") + find_path(LAPACKE_Generic_INCLUDE_DIR NAMES lapacke.h) + find_path(LAPACKE_Generic_CBLAS_INCLUDE_DIR NAMES cblas.h) + find_library(LAPACKE_Generic_LAPACKE_LIBRARY NAMES lapacke) + mark_as_advanced(LAPACKE_Generic_INCLUDE_DIR LAPACKE_Generic_CBLAS_INCLUDE_DIR + LAPACKE_Generic_LAPACKE_LIBRARY) + find_package(LAPACK QUIET) + if(LAPACKE_Generic_INCLUDE_DIR AND LAPACKE_Generic_CBLAS_INCLUDE_DIR + AND LAPACKE_Generic_LAPACKE_LIBRARY AND LAPACK_FOUND) + set(_ge_includes "${LAPACKE_Generic_INCLUDE_DIR}" "${LAPACKE_Generic_CBLAS_INCLUDE_DIR}") + list(REMOVE_DUPLICATES _ge_includes) + set(_ge_libs "${LAPACKE_Generic_LAPACKE_LIBRARY}" ${LAPACK_LIBRARIES}) + _lapacke_link_test(Generic "${_ge_includes}" "${_ge_libs}" "" _ok) + if(_ok) + set(_lapacke_found_vendor Generic) + set(_lapacke_includes ${_ge_includes}) + set(_lapacke_libs ${_ge_libs}) + endif() + endif() +endif() + +set(LAPACKE_VENDOR_FOUND "${_lapacke_found_vendor}") +set(LAPACKE_INCLUDE_DIRS "${_lapacke_includes}") +set(LAPACKE_LIBRARIES "${_lapacke_libs}") + +include(FindPackageHandleStandardArgs) +find_package_handle_standard_args(LAPACKE + REQUIRED_VARS LAPACKE_VENDOR_FOUND + REASON_FAILURE_MESSAGE + "No LAPACKE + CBLAS stack for LAPACKE_VENDOR=${LAPACKE_VENDOR}. Install one (Debian/Ubuntu: libopenblas-dev, or liblapacke-dev with liblapack-dev, or libmkl-dev with -DCBK_WITH_MKL=ON) or point CMAKE_PREFIX_PATH at it. The tests need one, the library itself does not (-DCBK_BUILD_TESTS=OFF).") + +if(LAPACKE_FOUND AND NOT TARGET LAPACKE::LAPACKE) + add_library(LAPACKE::LAPACKE INTERFACE IMPORTED) + set_target_properties(LAPACKE::LAPACKE PROPERTIES + INTERFACE_INCLUDE_DIRECTORIES "${LAPACKE_INCLUDE_DIRS}" + INTERFACE_LINK_LIBRARIES "${LAPACKE_LIBRARIES}" + INTERFACE_COMPILE_DEFINITIONS "${_lapacke_defs}") + message(STATUS "LAPACKE for the tests: ${LAPACKE_VENDOR_FOUND} (${LAPACKE_LIBRARIES})") +endif() diff --git a/docs/PLANS.md b/docs/PLANS.md index a4ca64e..c4b2eee 100644 --- a/docs/PLANS.md +++ b/docs/PLANS.md @@ -12,10 +12,11 @@ pivoting, and overflow/underflow-safe scaling are out of scope throughout. - **Implemented (design 6-8):** vectorized unblocked `geqr2` with a branch-free masked `larfg`; the trailing update is `ormqr`'s `larf`. One strided kernel for both layouts (column-major contiguous). -- **Validated (design 7):** BLAS-free test vs a scalar `geqr2`; MKL/LAPACK - test gating factorization residual (`20 n eps`) and orthogonality - (`100 n eps`), cross-checking `mkl_dgeqrf_compact` in both layouts, closing - `AX = B`, and covering rank-deficient / near-collinear inputs. +- **Validated (design 7):** portable test vs `LAPACKE_?geqr2` (elementwise) + gating factorization residual (`20 n eps`) and orthogonality (`100 n eps`) + of the `?orgqr`-formed `Q`, covering rank-deficient / near-collinear + inputs; MKL test cross-checking `mkl_dgeqrf_compact` in both layouts and + closing `AX = B`. - **Benchmarked:** `bench_geqrf_compact`. - **Scoped out (design 6.6):** `dlarfg` rescaling (the portable test pins the underflow behavior: `tau = 0`, diagonal kept, body zeroed), column pivoting, blocked @@ -27,8 +28,8 @@ pivoting, and overflow/underflow-safe scaling are out of scope throughout. sides, both layouts, `trans in {N, T}` (`C` folds to `T`). The shared `larf` is register-blocked four slices at a time; `side = 'R'` is the same kernel over the transposed view. Padded slots carry `tau = 0`, so they are no-ops. -- **Validated (design 7):** BLAS-free test vs a scalar `dorm2r` (including a - column-pivoted QR + back-permutation solve); MKL test of `op(Q) C` vs dense +- **Validated (design 7):** portable test vs `LAPACKE_?ormqr` (including a + `?geqp3` column-pivoted QR + back-permutation solve); MKL test of `op(Q) C` vs dense `LAPACKE_dormqr` over `layout x side x trans` (`20 s eps`) plus the `mkl_dgeqrf_compact -> cbk_dormqr_compact -> mkl_dtrsm_compact` solve (`100 n eps`). @@ -50,7 +51,7 @@ pivoting, and overflow/underflow-safe scaling are out of scope throughout. (`1 - tau` diagonal, `-tau`-scaled reflector body, zeros above). In place on the `?geqrf_compact` output (`m >= n >= k`); both layouts through the one strided kernel; no scratch (the MKL-style `lwork = -1` query answers 1). -- **Validated (design 7):** BLAS-free test vs a scalar `dorg2r` plus the +- **Validated (design 7):** portable test vs `LAPACKE_?orgqr` plus the independent invariants (`Q^T Q = I`, `Q R = A`, padded lanes exactly identity), covering `k < n`, `k = 0` and both layouts; MKL test vs dense `LAPACKE_dorgqr` (cross-check tolerance) gating orthogonality (`100 n eps`) @@ -75,14 +76,14 @@ compact `potrs` or `posv`. The solve as two non-unit `trsm` group sweeps (`?sytrsnp` minus the diagonal step); `posv` factoring and solving each group while its factor is cache-resident, bit-identical to the two calls. -- **Validated (design 7):** BLAS-free test vs a scalar `potf2` (both `uplo`, - both layouts, padding, non-SPD lane isolation), the end-to-end SPD solve +- **Validated (design 7):** portable test vs `LAPACKE_?potrf` (both `uplo`, + both layouts, padding, the `cond` knob, non-SPD lane isolation) gating the + factor (`20 n eps`, relative), the reconstruction residual (`20 n eps`) + and the untouched triangle (bit-for-bit), the end-to-end SPD solve (two-step and fused) and C-API validation of all three entry points; - MKL/LAPACK test gating the - reconstruction residual (`20 n eps`), the untouched triangle (bit-for-bit), - the factor vs `LAPACKE_dpotrf` (`20 n eps`), the cross-check vs - `mkl_dpotrf_compact` (bit-exact), the SPD solve through `mkl_?trsm_compact` - and through `potrs` (`100 n eps`), and the fused driver's bit-identity. + MKL test with the cross-check vs `mkl_dpotrf_compact` (bit-exact), the SPD + solve through `mkl_?trsm_compact` and through `potrs` (`100 n eps`), and + the fused driver's bit-identity. - **Benchmarked:** `bench_potrf_compact` (the factorization) and `bench_posv_compact` (the end-to-end solve: fused vs its own two-step calls vs MKL's compact `potrf + trsm x2` pipeline vs per-matrix `LAPACKE_dposv`). @@ -124,8 +125,10 @@ driver (`docs/cbk_dsytrfnp_compact_design.md`); MKL has no compact `sysvnp` factoring and solving each group while its factor is cache-resident, bit-identical to the two calls. The upper convention is `A = U^T D U` (design 6.3), not `?sytrf`'s `U D U^T`. -- **Validated (design 7):** BLAS-free test over `(T, V, uplo, layout)` with - padding, the end-to-end indefinite solve (two-step and fused), the zero-pivot +- **Validated (design 7):** portable test over `(T, V, uplo, layout)` with + padding (vs the one hand-rolled scalar reference left, itself validated + against `cblas_?trmm`), the end-to-end indefinite solve (two-step and + fused), the zero-pivot semantics of design 6.2, and C-API validation of all three entry points; MKL test gating reconstruction (`20 n eps`), the untouched triangle, an elementwise `(L, D)` cross-check vs `mkl_?getrfnp_compact` (`~3e-15` FP64), @@ -152,8 +155,8 @@ driver (`docs/cbk_dsytrfnp_compact_design.md`); MKL has no compact tuned (4/2/1 register-blocked row-dot, contiguous axpy for the leftover column); the group kernel routes on views, so a transposed column-major view reaches the tuned path too (used by gels). Other combinations are strided. -- **Validated (design 7):** BLAS-free test vs a scalar `?trsm` (forward error - and `||op(A) X - alpha B||`); MKL cross-check vs `mkl_?trsm_compact` over the +- **Validated (design 7):** portable test vs `cblas_?trsm` (forward error + and `||op(A) X - alpha B||` by `cblas_?trmm`); MKL cross-check vs `mkl_?trsm_compact` over the full feature matrix plus the end-to-end solve. - **Performance vs `mkl_?trsm_compact`** (single thread, AVX-512, orders 10-148, measured before the view-based routing): `nrhs = 1` ~`1.0x`, `nrhs` @@ -176,8 +179,8 @@ driver (`docs/cbk_dsytrfnp_compact_design.md`); MKL has no compact sides. A fused factor-and-apply kernel for the least-squares case was retired (design 6.2): no arithmetic saved, and it evicted `B` from L1 at mid sizes. -- **Validated (design 7):** BLAS-free test vs a scalar `ref_gels` (`X`, - factorization and `tau` elementwise) plus the defining properties formed +- **Validated (design 7):** portable test vs `LAPACKE_?gels` and + `?geqrf`/`?gelqf` (`X`, factorization and `tau`, relative) plus the defining properties formed independently, C-API validation, and the `min(m,n) = 0` quick return; MKL test vs per-matrix `LAPACKE_?gels` over every `(layout, trans)` and square/tall/wide shape (forward error `100 max(m,n) eps`, residual @@ -202,14 +205,26 @@ driver (`docs/cbk_dsytrfnp_compact_design.md`); MKL has no compact ## Project-wide - **Precision coverage.** Every suite is templated on the scalar type and runs - in FP64 and FP32 through the `compat` / `mkl` / `lapack` dispatch - of `tests/test_mkl_util.hpp`; FP32 cross-checks agree with `mkl_s*_compact` - to ~1e-6, gated at 1e-4. -- **Planned: test against a real BLAS/LAPACK (issue #27).** The portable - suites hand-roll their scalar reference routines. Instead, assume a library - is present for testing -- `find_package(LAPACK REQUIRED)` -- and validate - against it, rather than maintaining our own reference versions. This also - opens the dense cross-checks (today MKL-only) to any BLAS/LAPACK stack. + in FP64 and FP32 through the `compact` / `lapack` dispatch of + `tests/test_compact_util.hpp` / `test_lapack_util.hpp` and the `compat` + / `mkl` one of `tests/test_mkl_util.hpp`; FP32 cross-checks agree with + `mkl_s*_compact` to ~1e-6, gated at 1e-4. +- **Done: test against a real BLAS/LAPACK (issue #27).** The suites validate + against LAPACKE + CBLAS (`tests/test_lapack_util.hpp`): the `ref_*` + procedures forward to `LAPACKE_?geqr2` / `?ormqr` / `?orgqr` / `?geqp3` / + `?potrf` / `?gels` and `cblas_?trsm` / `?trmm` / `?gemm`, with the + unpivoted LDL^T (`ref_sytf2np`, which LAPACK lacks) the one hand-rolled + reference left, validated against `cblas_?trmm` in its suite. + `cmake/FindLAPACKE.cmake` takes MKL's LAPACKE with the MKL extension and + OpenBLAS, Netlib or any `lapacke` + `find_package(LAPACK)` pair without + (`-DCBK_TEST_LAPACK`); CI runs OpenBLAS and Netlib under gcc and clang. The + dense cross-checks of the MKL suites (`geqrf`'s and `potrf`'s section 7.1) + moved into the portable suites, on every stack. Where a reference is + LAPACK's blocked driver, the gate is relative to the operand norms at a + multiple of `n eps` (design 7.1), not elementwise at `~eps`; only `?geqr2` + stays an elementwise `~eps` comparison. macOS through Accelerate + + accelerate-lapacke has a tag-only workflow (`.github/workflows/macos.yml`), + not yet exercised. - **Planned: LAPACK-style test coverage.** Adopt the testing approaches of the reference LAPACK repository (its `TESTING/LIN` drivers): `?latms`-style generators with prescribed condition number and spectral distribution, and diff --git a/docs/building.md b/docs/building.md index 10b827b..a68f854 100644 --- a/docs/building.md +++ b/docs/building.md @@ -3,9 +3,23 @@ ## The portable library (default) Requires **CMake >= 3.28**, a **C++17 compiler** (GCC/Clang) and a build tool -(Make/Ninja). Nothing else: the default configure builds the `cbk::cbk` -library (the C API of `include/cbk.h`) and its BLAS-free test suites, -which check every kernel against a scalar reference of the same algorithm. +(Make/Ninja). The library needs nothing else. The test suites validate every +kernel against a real **LAPACKE + CBLAS** (issue #27), so the default +configure, which builds them, also needs one of + +* OpenBLAS -- Debian/Ubuntu: `sudo apt-get install libopenblas-dev liblapacke-dev` + (Debian ships OpenBLAS's LAPACKE as Netlib's `liblapacke`), or +* Netlib's reference LAPACK -- `sudo apt-get install liblapacke-dev liblapack-dev libblas-dev`, + or +* Intel MKL's own, with the MKL extension below (then it is the only stack the + tests may link). + +`cmake/FindLAPACKE.cmake` picks the first it finds (OpenBLAS, then Netlib, +then any `lapacke` library next to what `find_package(LAPACK)` returns) and +reports it as `LAPACKE for the tests: ...`; `-DCBK_TEST_LAPACK=openblas|netlib|mkl` +pins one, which matters on Debian, where `liblapack.so` is an +`update-alternatives` symlink and either stack can otherwise stand in for the +other. `-DCBK_BUILD_TESTS=OFF` drops the requirement with the tests. ```sh cmake -S . -B build -DCMAKE_BUILD_TYPE=Release @@ -13,13 +27,20 @@ cmake --build build -j ctest --test-dir build --output-on-failure ``` +On macOS, Accelerate provides LAPACK and CBLAS but no LAPACKE; +[accelerate-lapacke](https://github.com/lepus2589/accelerate-lapacke) builds +one on top of it, and `-DCBK_TEST_LAPACK=accelerate` with +`-DCMAKE_PREFIX_PATH` at its install looks for that (the tag-only +`.github/workflows/macos.yml`; not yet exercised on a macOS machine). + ## The Intel MKL extension (`-DCBK_WITH_MKL=ON`) `-DCBK_WITH_MKL=ON` adds the MKL-style API (the `cbk_?*_compact` entry points of `include/cbk_compat.h`, compiled into the same `cbk::cbk` library), the MKL-backed -test suites (cross-checked against MKL's own compact kernels and dense -LAPACK/LAPACKE), the worked example and the benchmarks. It needs **Intel MKL**, +test suites (cross-checked against MKL's own compact kernels), the worked +example and the benchmarks; every suite's dense references then come from +MKL's LAPACKE. It needs **Intel MKL**, which provides the Compact-format API (`mkl_compact.h`, `mkl_?gepack_compact`, `mkl_?geqrf_compact`, ...). Any MKL works: @@ -125,7 +146,8 @@ CI builds against a fresh install for every MKL / static / shared combination |--------|---------|--------| | `CBK_WITH_MKL` | `OFF` | Build the MKL-style API, the MKL-backed tests, the example and the benchmarks; requires Intel MKL. | | `CBK_WITH_OPENMP` | `ON` | Thread the loop over groups with OpenMP; `OFF` gives single-threaded routines. See [threading.md](threading.md). | -| `CBK_BUILD_TESTS` | `ON` | Build the tests and register them with CTest. | +| `CBK_BUILD_TESTS` | `ON` | Build the tests and register them with CTest; needs a LAPACKE + CBLAS stack. | +| `CBK_TEST_LAPACK` | `auto` | The stack the tests validate against: `auto` (MKL's with the MKL extension, else the first of OpenBLAS, Netlib, any), `mkl`, `openblas`, `netlib`, `accelerate`. | | `BUILD_SHARED_LIBS` | `OFF` | CMake's own switch: build `libcbk` shared instead of static. | | `CBK_INSTALL` | `ON` when top-level | Generate the install rules; off by default under a parent project's `add_subdirectory`. | | `MKLCompact_ROOT` | -- | A oneAPI MKL prefix, when `MKLROOT` is not set (MKL build only). | diff --git a/docs/cbk_dgels_compact_design.md b/docs/cbk_dgels_compact_design.md index fd0724b..6148f9c 100644 --- a/docs/cbk_dgels_compact_design.md +++ b/docs/cbk_dgels_compact_design.md @@ -212,10 +212,11 @@ in FP64 and FP32. `mkl_?geqrf_compact -> cbk_?ormqr_compact -> mkl_?trsm_compact` must agree in `ap`, `tau` (elementwise, cross-check tolerance) and `X` (`100 * n * eps`), and both must match the known `X`. -* **Portable self-test** (no BLAS): the templated kernel across - `(T, V, layout, trans, shape)`, padded final packs, `nrhs = 1` and the - smallest sizes, against a scalar `ref_gels` of the same unblocked steps - (`X`, factorization and `tau` gated at `~eps`); the defining properties +* **Portable self-test** (`tests/test_gels_compact.cpp`, on any LAPACKE + stack): the templated kernel across `(T, V, layout, trans, shape)`, padded + final packs, `nrhs = 1` and the smallest sizes, against `LAPACKE_?gels` and + `?geqrf` / `?gelqf` through `ref_gels` (`X`, factorization and `tau`, + relative, at `100 max(m,n) eps`); the defining properties formed without the reference (normal equations and residual rows for least squares; `op(A) X = B` and agreement with `X = op(A)^T Z`, `(op(A) op(A)^T) Z = B`, for minimum norm); the portable C API's argument diff --git a/docs/cbk_dgeqrf_compact_design.md b/docs/cbk_dgeqrf_compact_design.md index 755b076..0746ca7 100644 --- a/docs/cbk_dgeqrf_compact_design.md +++ b/docs/cbk_dgeqrf_compact_design.md @@ -258,16 +258,17 @@ must recover `X`. Gate the forward error `Xhat - X` and the residual ### 7.4 Portable self-test -A self-contained test validates the templated kernel directly against a scalar -reference (`ref_geqr2` / `ref_larfg`) across `(T, V)` combinations and partial -(padded) final packs, the underflow scope of 6.6 (a column whose tail -underflows squared comes back with `tau = 0`, its diagonal kept and its body -zeroed, every lane finite, and the factorization still exact to working -precision), plus the -LAPACK-style argument validation of the portable C API. It needs no external libraries at all; only the suites above require an -MKL installation (for the Compact API). BLAS and LAPACK themselves are assumed -available, as they are on most platforms -- it is the MKL Compact extension that -must be installed separately. +`tests/test_geqrf_compact.cpp` validates the templated kernel directly against +LAPACK's unblocked `?geqr2` (`LAPACKE_?geqr2` through `ref_geqr2`, +elementwise at `~eps`) across `(T, V)` combinations and partial (padded) final +packs, the underflow scope of 6.6 (a column whose tail underflows squared +comes back with `tau = 0`, its diagonal kept and its body zeroed, every lane +finite, and the factorization still exact to working precision), the +dense-LAPACK invariants of 7.1 (`Q` from `?orgqr`, residual and orthogonality +gates, over the structured inputs) on any LAPACKE stack, plus the LAPACK-style +argument validation of the portable C API. It needs a LAPACKE + CBLAS +(`cmake/FindLAPACKE.cmake`: OpenBLAS, Netlib, or MKL's own); only the suites +above require an MKL installation (for the Compact API). ## 8. Implementation Strategy diff --git a/docs/cbk_dorgqr_compact_design.md b/docs/cbk_dorgqr_compact_design.md index 01719c5..176db96 100644 --- a/docs/cbk_dorgqr_compact_design.md +++ b/docs/cbk_dorgqr_compact_design.md @@ -234,13 +234,14 @@ CTest-registered and run in FP64 and FP32. same reflectors, same application order, different structure exploitation. This pins the accumulation against an already-validated routine without any dense reference in the loop. -* **Portable self-test** (no BLAS): the templated kernel across - `(T, V, layout, shape)` including `k < n`, `k = 0` (`Q` = leading - columns of `I`), `min(m, n) = 0`, and padded final packs (padded lanes - come out identity, real lanes unperturbed), against a scalar - `ref_org2r` of the same steps (elementwise at `~eps`); orthogonality - and `Q R = A` formed without the reference; and the portable C API's - argument validation, `-j` per the header. +* **Portable self-test** (`tests/test_orgqr_compact.cpp`, on any LAPACKE + stack): the templated kernel across `(T, V, layout, shape)` including + `k < n`, `k = 0` (`Q` = leading columns of `I`), `min(m, n) = 0`, and + padded final packs (padded lanes come out identity, real lanes + unperturbed), against `LAPACKE_?orgqr` (elementwise at `100 m eps`: + `Q` is unique for full-rank `A`, LAPACK blocks its accumulation above + `k = 128`); orthogonality and `Q R = A` formed without the reference; and + the portable C API's argument validation, `-j` per the header. ## 8. Implementation diff --git a/docs/cbk_dpotrf_compact_design.md b/docs/cbk_dpotrf_compact_design.md index a942e9f..c895032 100644 --- a/docs/cbk_dpotrf_compact_design.md +++ b/docs/cbk_dpotrf_compact_design.md @@ -440,15 +440,17 @@ with the rest of the compact toolkit. ### 7.4 Portable self-test -A self-contained test validates the templated kernels directly against a scalar -reference (`ref_potf2`) across `(T, V)` combinations, both `uplo`, and partial -(padded) final packs, plus the end-to-end portable solve (`?potrf_compact` + -`?potrs_compact` recovering a known `X`), the fused `?posv_compact` -bit-identical to it, and the LAPACK-style argument validation of the three -portable C APIs. It needs no external libraries at all; only the suites above -require an MKL installation (for the Compact API). BLAS and LAPACK themselves -are assumed available, as they are on most platforms -- it is the MKL Compact -extension that must be installed separately. +`tests/test_potrf_compact.cpp` validates the templated kernels directly +against `LAPACKE_?potrf` (through `ref_potf2`; the unique SPD factor, +relative to its norm at `20 n eps`) across `(T, V)` combinations, both +`uplo`, both layouts, the `cond` knob, and partial (padded) final packs, with +the 7.1 invariants (reconstruction by `?trmm`, the untouched triangle) and +the non-SPD lane isolation of 6.2, plus the end-to-end portable solve +(`?potrf_compact` + `?potrs_compact` recovering a known `X`), the fused +`?posv_compact` bit-identical to it, and the LAPACK-style argument validation +of the three portable C APIs. It needs a LAPACKE + CBLAS +(`cmake/FindLAPACKE.cmake`: OpenBLAS, Netlib, or MKL's own); only the suites +above require an MKL installation (for the Compact API). ### 7.5 Suite 4 -- The packaged solve: `?potrs` and the fused `?posv` @@ -495,5 +497,5 @@ through `extern "C"` for the FFI-stable surfaces, built on the project's | `src/cbk_posv_compact.hpp` | The fused factor-and-solve driver over both group kernels. | | `src/cbk.cpp` | Portable `?potrf_compact` / `?potrs_compact` / `?posv_compact` C entry points (runtime `V` -> compile-time dispatch, `info = -j`). | | `src/cbk_compat.cpp` | Unwraps `MKL_COMPACT_PACK` -> `V` and `MKL_UPLO`/`MKL_LAYOUT`, calls the kernels. | -| `tests/test_potrf_compact.cpp` | Self-contained correctness test vs a scalar `potf2` reference (no BLAS), incl. the potrs/posv solve. | -| `tests/test_potrf_mkl.cpp` | MKL + dense-LAPACK validation (residual, uniqueness, cross-check, solves, fused bit-identity). | +| `tests/test_potrf_compact.cpp` | Correctness test vs `LAPACKE_?potrf` on any LAPACKE stack (residual, uniqueness, untouched triangle), incl. the potrs/posv solve. | +| `tests/test_potrf_mkl.cpp` | MKL validation (cross-check vs `mkl_?potrf_compact`, solves through MKL's `trsm` and `potrs`/`posv`, fused bit-identity). | diff --git a/docs/cbk_dsytrfnp_compact_design.md b/docs/cbk_dsytrfnp_compact_design.md index e84c23c..42230eb 100644 --- a/docs/cbk_dsytrfnp_compact_design.md +++ b/docs/cbk_dsytrfnp_compact_design.md @@ -391,9 +391,12 @@ suites also gate the `nrhs = 0` contract of section 3: `?sysvnp` must factor ### 7.4 Portable self-test -A self-contained test (no BLAS) validates the templated kernels directly -against a scalar reference across `(T, V)` combinations, both `uplo`, both -layouts, and padded final packs, plus: +`tests/test_sytrfnp_compact.cpp` validates the templated kernels directly +against a scalar reference (`ref_sytf2np`: LAPACK has no unpivoted LDL^T, so +this is the one reference the suites still hand-roll, and it is itself +checked first -- its `(L, D)` must reconstruct `A` through `cblas_?trmm` at +`20 n eps`) across `(T, V)` combinations, both `uplo`, both layouts, and +padded final packs, plus: * the **end-to-end portable solve** (`?sytrfnp_compact` + `?sytrsnp_compact`) and the **fused** `?sysvnp_compact`, bit-identical to it, @@ -441,7 +444,7 @@ kernels. | `src/cbk_sysvnp_compact.hpp` | The fused factor-and-solve driver over both group kernels. | | `src/cbk.cpp` | Portable `?sytrfnp_compact` / `?sytrsnp_compact` / `?sysvnp_compact` C entry points (runtime `V` dispatch, `info = -j`). | | `src/cbk_compat.cpp` | The MKL-style adapters (`MKL_COMPACT_PACK` -> `V`, `MKL_UPLO`/`MKL_LAYOUT` -> flags). | -| `tests/test_sytrfnp_compact.cpp` | Self-contained correctness test vs a scalar reference (no BLAS), FP64 and FP32. | +| `tests/test_sytrfnp_compact.cpp` | Correctness test vs the scalar reference (itself validated against CBLAS), FP64 and FP32. | | `tests/test_sytrfnp_mkl.cpp` | MKL validation (invariants, `getrfnp` cross-check, indefinite solve, fused bit-identity), FP64 and FP32. | The prototypes are added to `cbk_compat.h` (MKL-style) and `cbk.h` diff --git a/docs/cbk_dtrsm_compact_design.md b/docs/cbk_dtrsm_compact_design.md index 0ede3f3..ac86306 100644 --- a/docs/cbk_dtrsm_compact_design.md +++ b/docs/cbk_dtrsm_compact_design.md @@ -249,10 +249,11 @@ must satisfy the same numerical invariants as an unbatched `?trsm`. Tolerances are purely relative to the working precision. All suites below are CTest-registered. -### 7.1 Suite 1 -- Portable, vs a scalar `?trsm` reference (no BLAS) +### 7.1 Suite 1 -- Portable, vs `cblas_?trsm` -`tests/test_trsm_compact.cpp` validates the templated kernel directly against a -scalar `?trsm` (the same algorithm, one matrix at a time) over the full +`tests/test_trsm_compact.cpp` validates the templated kernel directly against +the library's `?trsm` (`cblas_?trsm` through `ref_trsm`, on any LAPACKE + +CBLAS stack: the same algorithm, one matrix at a time) over the full `side x uplo x transa x diag` matrix, across precisions (FP32/FP64), interleave widths, and padded final packs, column-major. The row-dot kernels run the identical operation sequence one lane per matrix and so match the reference to @@ -261,9 +262,9 @@ accumulates in a different (but equally backward-stable) order and agrees to working precision (~`1e-16`). Two gates, both at a generous multiple of `eps` on the diagonal-boosted (well conditioned) factors: the relative forward error of `X` against the reference -solve, and the solve's own residual `||op(A) X - alpha B||` formed with an -independent triangular multiply -- so a bug shared by the scalar reference and -the kernel cannot pass unseen (the `?trsm` analogue of the reconstruction check +solve, and the solve's own residual `||op(A) X - alpha B||` formed with the +library's triangular multiply (`cblas_?trmm`) -- so a bug shared by the +reference and the kernel cannot pass unseen (the `?trsm` analogue of the reconstruction check the `geqrf`/`potrf` self-tests apply). ### 7.2 Suite 2 -- MKL cross-check, vs `mkl_?trsm_compact` diff --git a/docs/layout.md b/docs/layout.md index dd44450..92275b5 100644 --- a/docs/layout.md +++ b/docs/layout.md @@ -3,11 +3,13 @@ ``` include/ public headers src/ kernels and the two adapter sources -tests/ portable (no BLAS) and MKL-backed test suites; tests/install/ checks - an installed cbk from a downstream project +tests/ the portable C API's suites (vs LAPACKE/CBLAS, any stack) and the + MKL-backed ones; tests/install/ checks an installed cbk from a + downstream project examples/ worked solve + benchmarks (see examples/BENCHMARKS.md) docs/ design documents and guides (see docs/README.md) -cmake/ FindMKLCompact.cmake, cbkConfig.cmake.in (the installed package) +cmake/ FindMKLCompact.cmake, FindLAPACKE.cmake (the tests' stack), + cbkConfig.cmake.in (the installed package) ``` ## Public interface @@ -45,10 +47,11 @@ scalar, `0` on success). | `src/cbk_gels_compact.hpp` | The one-call least-squares / minimum-norm solve: one driver over the geqrf, ormqr and trsm group kernels, on the tall view of `A` (transposed when `m < n`, which is the LQ case). | | `src/cbk.cpp` | The portable C API: argument validation and `V` dispatch for all twenty-two entry points. | | `src/cbk_compat.cpp` | The MKL-style API: MKL enum / `MKL_COMPACT_PACK` unwrapping for all twenty-two entry points. | -| `tests/test_compact_util.hpp` | Shared test helpers: RNG, error metrics, input generation, scalar reference kernels, Compact pack/unpack, and `compact` (the portable C API dispatched on the scalar type); header-only, no MKL. | -| `tests/test_mkl_util.hpp` | Scalar-type dispatch for the MKL-backed suites: `compat` (routines under test), `mkl` (MKL's compact API and kernels), `lapack` (LAPACKE/CBLAS references). | -| `tests/test_*_compact.cpp` | Portable self-contained suites (no BLAS): each kernel vs its scalar reference, plus C API validation, in FP64 and FP32. | -| `tests/test_*_mkl.cpp` | MKL + dense-LAPACK validation, templated on the scalar type and run in FP64 and FP32: invariants vs LAPACK, cross-checks vs MKL's compact kernels, end-to-end solves. | +| `tests/test_compact_util.hpp` | Shared test helpers: RNG, error metrics, input generation, Compact pack/unpack, and `compact` (the portable C API dispatched on the scalar type); header-only, library-free. | +| `tests/test_lapack_util.hpp` | The dense references on a real LAPACKE + CBLAS: `lapack` (the LAPACKE/CBLAS dispatch), the `ref_*` procedures over views (`ref_geqr2`, `ref_orm2r`, `ref_potf2`, `ref_trsm`, `ref_gels`, ...), `matmul`, `tri_apply`, `solve_errors`. The stack comes from `cmake/FindLAPACKE.cmake`. | +| `tests/test_mkl_util.hpp` | Scalar-type dispatch for the MKL-backed suites: `compat` (routines under test), `mkl` (MKL's compact API and kernels). | +| `tests/test_*_compact.cpp` | The portable C API's suites, on any LAPACKE stack: each kernel vs its LAPACK counterpart (elementwise where that is the same unblocked algorithm, residual and invariant gates elsewhere), plus C API validation, in FP64 and FP32. | +| `tests/test_*_mkl.cpp` | MKL validation, templated on the scalar type and run in FP64 and FP32: the MKL-style API through MKL's pack/unpack, cross-checks vs MKL's compact kernels, end-to-end solves. | | `examples/bench_util.hpp` | The benchmarks' shared harness (timing, aligned storage, command line). | The contributor conventions behind this layout -- one kernel per routine over diff --git a/tests/test_compact_util.hpp b/tests/test_compact_util.hpp index 46598ad..4e30c39 100644 --- a/tests/test_compact_util.hpp +++ b/tests/test_compact_util.hpp @@ -1,9 +1,10 @@ // test_compact_util.hpp // // Shared helpers for the compact-format test suites: a seeded RNG, error -// metrics, input generation, scalar LAPACK-style reference kernels, precision -// overloads of the portable C API, and Compact pack/unpack. Header-only and -// MKL-free, so the BLAS-free portable tests use it too. +// metrics, input generation, precision overloads of the portable C API, and +// Compact pack/unpack. Header-only and library-free; the dense references +// (LAPACKE / CBLAS through lapack) are in test_lapack_util.hpp, which +// includes this, and the MKL-side dispatch in test_mkl_util.hpp. // // Assisted-by: Claude:claude-opus-4-8 Claude:claude-fable-5 @@ -105,7 +106,7 @@ template double norm1(Mv M) // Deviation of Q's columns from orthonormality, max |(Q^T Q - I)(i,j)| over // the upper triangle (Q^T Q is symmetric), accumulated in double so the gate -// does not inherit float rounding. BLAS-free; shared by the orgqr suites. +// does not inherit float rounding. Shared by the orgqr suites. template double orth_error(Qv Q) { double e = 0; @@ -119,21 +120,6 @@ template double orth_error(Qv Q) return e; } -// C (m x n) := A (m x k) * B (k x n) -- the plain triple loop, for forming -// right-hand sides and residuals. The shapes come from the views, which also -// carry the layout, so the three operands need not share one. -template void matmul(Av A, Bv B, Cv C) -{ - assert(A.rows == C.rows && B.cols == C.cols && A.cols == B.rows); - for (int j = 0; j < C.cols; ++j) - for (int i = 0; i < C.rows; ++i) { - elem_t s = 0; - for (int l = 0; l < A.cols; ++l) - s += A(i, l) * B(l, j); - C(i, j) = s; - } -} - // Elementwise copy dst <- src, same shape, each side in its own layout -- the // staging step between a column-major batch and a row-major API surface. template void copy_matrix(Sv src, MatrixView dst) @@ -156,191 +142,6 @@ template std::vector known_solution(int n, int nrhs) return Xs; } -// ----------------------- reference kernels (scalar) ----------------- -// Column-major, in place: the unblocked LAPACK algorithms the vectorized -// kernels execute V lanes at a time, so a correct kernel matches them to -// working precision. Matrix operands are views; the vector operands of ?larfg -// stay pointers, since a reflector is a contiguous column segment and not a -// matrix (which is also why these references are column-major only -- the -// views they are handed must have unit row stride, as their col() asserts). - -// dlarfg: reflector from (alpha, x[0..m-2]); alpha := beta on exit. -template void ref_larfg(int m, T *alpha, T *x, T *tau) -{ - T xnorm = 0; - for (int i = 0; i < m - 1; ++i) - xnorm = std::hypot(xnorm, x[i]); - if (xnorm == T(0)) { - *tau = 0; - return; - } - T beta = -std::copysign(std::hypot(*alpha, xnorm), *alpha); - *tau = (beta - *alpha) / beta; - T scal = T(1) / (*alpha - beta); - for (int i = 0; i < m - 1; ++i) - x[i] *= scal; - *alpha = beta; -} - -// dgeqr2: unblocked Householder QR, (H, tau) in the LAPACK convention. -template void ref_geqr2(MatrixView A, T *tau) -{ - const int m = A.rows, n = A.cols, k = std::min(m, n); - for (int kk = 0; kk < k; ++kk) { - // The reflector runs down column kk from the diagonal; A.col() gives - // the contiguous column, so the tail below the diagonal is addressable - // even when it is empty (kk + 1 == m on the last row). - ref_larfg(m - kk, &A(kk, kk), A.col(kk) + kk + 1, &tau[kk]); - for (int j = kk + 1; j < n; ++j) { - T w = A(kk, j); - for (int i = kk + 1; i < m; ++i) - w += A(i, kk) * A(i, j); - A(kk, j) -= tau[kk] * w; - for (int i = kk + 1; i < m; ++i) - A(i, j) -= tau[kk] * A(i, kk) * w; - } - } -} - -// dorm2r, side='L': B := Q^T B (trans 'T') or Q B ('N') from (H, tau), with -// k reflectors read from A (m x k or wider) and applied to B (m x nrhs). -template -void ref_orm2r(char trans, int k, Av A, const T *tau, MatrixView B) -{ - const int m = B.rows, nrhs = B.cols; - assert(A.rows == m && k <= std::min(A.rows, A.cols)); - bool fwd = (trans == 'T'); - for (int s = 0; s < k; ++s) { - int kk = fwd ? s : k - 1 - s; - for (int j = 0; j < nrhs; ++j) { - T w = B(kk, j); - for (int i = kk + 1; i < m; ++i) - w += A(i, kk) * B(i, j); - B(kk, j) -= tau[kk] * w; - for (int i = kk + 1; i < m; ++i) - B(i, j) -= tau[kk] * A(i, kk) * w; - } - } -} - -// Back substitution R X = B with R the upper triangle of an n x n array. -template void ref_trsm_upper(Rv R, MatrixView B) -{ - const int n = B.rows, nrhs = B.cols; - assert(R.rows == n && R.cols == n); - for (int j = 0; j < nrhs; ++j) - for (int i = n - 1; i >= 0; --i) { - T s = B(i, j); - for (int l = i + 1; l < n; ++l) - s -= R(i, l) * B(l, j); - B(i, j) = s / R(i, i); - } -} - -// Dense BLAS ?trsm: solves op(A) X = alpha B (side='L') or X op(A) = alpha B -// (side='R') in place, A the order-s triangular factor. B is pre-scaled by -// alpha (so alpha == 0 gives B := 0), then a unit-alpha substitution runs. Only -// the referenced triangle of A is touched; the diagonal is skipped entirely -// when diag='U'. The ?trsm suite's reference, and the triangular steps of -// ref_gels. -template -void ref_trsm(char side, char uplo, char transa, char diag, T alpha, Av A, - MatrixView B) -{ - const bool left = (side == 'L' || side == 'l'); - const bool upper = (uplo == 'U' || uplo == 'u'); - const bool tran = (transa == 'T' || transa == 't' || transa == 'C' || transa == 'c'); - const bool unit = (diag == 'U' || diag == 'u'); - const int m = B.rows, n = B.cols; - assert(A.rows == A.cols && A.rows == (left ? m : n)); - - for (int j = 0; j < n; ++j) - for (int i = 0; i < m; ++i) - B(i, j) *= alpha; // B := alpha B (alpha == 0 -> B := 0) - - if (left) { - // solve op(A) X = B column by column; A is m x m - const bool back = (upper != tran); - for (int j = 0; j < n; ++j) - for (int t = 0; t < m; ++t) { - int i = back ? m - 1 - t : t; - T s = B(i, j); - if (back) - for (int l = i + 1; l < m; ++l) - s -= (tran ? A(l, i) : A(i, l)) * B(l, j); - else - for (int l = 0; l < i; ++l) - s -= (tran ? A(l, i) : A(i, l)) * B(l, j); - B(i, j) = unit ? s : s / A(i, i); - } - } - else { - // solve X op(A) = B column of X at a time; A is n x n - const bool fwd = (upper != tran); - for (int t = 0; t < n; ++t) { - int j = fwd ? t : n - 1 - t; - if (fwd) - for (int l = 0; l < j; ++l) { - T a = tran ? A(j, l) : A(l, j); - for (int i = 0; i < m; ++i) - B(i, j) -= a * B(i, l); - } - else - for (int l = j + 1; l < n; ++l) { - T a = tran ? A(j, l) : A(l, j); - for (int i = 0; i < m; ++i) - B(i, j) -= a * B(i, l); - } - if (!unit) { - T d = A(j, j); - for (int i = 0; i < m; ++i) - B(i, j) /= d; - } - } - } -} - -// dgels, unblocked: the least-squares (op(A) with more rows than columns) or -// minimum-norm (more columns than rows) solution of op(A) X = B, op(A) = A -// ('N') or A^T ('T'), A m x n and B max(m,n) x nrhs, each in the layout its -// view carries. Factors the tall orientation F (A, or A^T when m < n: that is -// the LQ of A in ?gelqf storage) with ref_geqr2 -- on exit A holds the -// factorization and tau its min(m,n) reflector scalars -- then B := Q^T B, -// R X = B (least squares: rows n..m-1 of B keep the residual) or R^T Y = B, -// X = Q [Y; 0] (minimum norm). The same steps the compact kernel runs V lanes -// at a time. -template void ref_gels(char trans, MatrixView A, MatrixView B, T *tau) -{ - const bool tran = (trans == 'T' || trans == 't' || trans == 'C' || trans == 'c'); - const int m = A.rows, n = A.cols, nrhs = B.cols; - const bool tall = (m >= n); - const int p = tall ? m : n, q = tall ? n : m; - const bool overdet = (tall != tran); - assert(B.rows == p); - - // F: the tall orientation, p x q, column-major (the references need it so) - std::vector Fs((size_t)p * q); - const auto F = mat_view(Fs.data(), p, q); - copy_matrix(tall ? A : A.transposed(), F); - - ref_geqr2(F, tau); - const auto R = leading(F, q, q), Bq = leading(B, q, nrhs); - if (overdet) { - ref_orm2r('T', q, F, tau, B); - ref_trsm('L', 'U', 'N', 'N', T(1), R, Bq); - } - else { - ref_trsm('L', 'U', 'T', 'N', T(1), R, Bq); - for (int j = 0; j < nrhs; ++j) - for (int i = q; i < p; ++i) - B(i, j) = T(0); - ref_orm2r('N', q, F, tau, B); - } - - // hand the factorization back in A's orientation - copy_matrix(F, tall ? A : A.transposed()); -} - // ----------------------- portable C API, by scalar type ------------ // compact::geqrf / ormqr / orgqr / potrf / potrs / posv / sytrfnp / sytrsnp // / sysvnp / trsm / gels forward to the d/s entry points of cbk.h, so @@ -484,7 +285,9 @@ template void gen_sym_ldlt(MatrixView A) // A(i,j) = sum_{l <= min(i,j)} F(i,l) D(l) F(j,l) with F the unit factor read // from the strict lower triangle (uplo 'L') or, for the transpose dual // A = U^T D U, from the strict upper one (F(i,l) = U(l,i)). `at(i,j)` reads the -// stored factor. Shared by the reconstruction checks of both ?sytrfnp suites. +// stored factor. Shared by the reconstruction checks of both ?sytrfnp suites +// (LAPACK has no unpivoted LDL^T to hand this to; the portable suite checks +// it against a ?trmm-formed reconstruction once). template double ldlt_reconstruct(const At &at, int i, int j, bool upper) { double s = 0; @@ -513,43 +316,6 @@ template void gen_tri(MatrixView A, bool upper) A(d, d) = (A(d, d) >= 0 ? T(1) : T(-1)) * (T(2) + std::abs(frand())); } -// Apply a triangular operator to a general matrix -- the "forward" direction of -// a ?trsm, for checking a solve's defining residual ||op(A) X - alpha B||. -// R (m x n) := op(A) X (side 'L') or X op(A) (side 'R'), with A the order-s -// (s = m for 'L', n for 'R') triangular factor: uplo 'U'/'L', op(A) = A ('N') -// or A^T ('T'/'C'), unit ('U') or non-unit ('N') diagonal. The shapes come -// from R and A; each operand carries its own layout in its view. -template -void tri_apply(char side, char uplo, char transa, char diag, Av A, Xv X, Rv R) -{ - using T = elem_t; - const bool left = (side == 'L' || side == 'l'); - const bool upper = (uplo == 'U' || uplo == 'u'); - const bool tran = (transa == 'T' || transa == 't' || transa == 'C' || transa == 'c'); - const bool unit = (diag == 'U' || diag == 'u'); - const int m = R.rows, n = R.cols, s = left ? m : n; - assert(A.rows == s && A.cols == s && X.rows == m && X.cols == n); - // op(A)(i,k): unit or A(i,i) on the diagonal; off-diagonal is the referenced - // entry of A (op = A) or of its transpose (op = A^T), else zero. - auto Mop = [&](int i, int k) -> T { - if (i == k) return unit ? T(1) : A(i, i); - const bool ref = tran ? (upper ? i > k : i < k) : (upper ? k > i : k < i); - if (!ref) return T(0); - return tran ? A(k, i) : A(i, k); - }; - for (int j = 0; j < n; ++j) - for (int i = 0; i < m; ++i) { - T acc = 0; - if (left) - for (int k = 0; k < s; ++k) - acc += Mop(i, k) * X(k, j); - else - for (int k = 0; k < s; ++k) - acc += X(i, k) * Mop(k, j); - R(i, j) = acc; - } -} - // ----------------------- dense batches and compact packing ----------- // The dense batches are MatrixBatch (src/cbk_matrix_batch.hpp), shared with // the benchmarks; the suites allocate with the default allocator. @@ -674,34 +440,6 @@ template void unpack_tau(MatrixBatch &tau, const T *tp, int V) // ----------------------- shared checks and epilogues ----------------- -// The end-to-end solve gate every solving suite closes with: the worst -// relative forward error of Xhat against the known X, and the worst relative -// residual ||A Xhat - B|| formed densely -- so a solver bug cannot hide behind -// the factorization that produced Xhat. -struct SolveErrors { - double fwd, res; -}; - -template -SolveErrors solve_errors(const MatrixBatch &A, const MatrixBatch &B, - const MatrixBatch &Xhat, ConstMatrixView X) -{ - assert(B.rows() == Xhat.rows() && B.cols() == Xhat.cols() && X.rows == Xhat.rows() && - X.cols == Xhat.cols()); - const size_t sB = Xhat.stride(); - std::vector AXs(sB); - const auto AX = mat_view(AXs.data(), Xhat.rows(), Xhat.cols()); - const double nX = std::max(norm1(X), norm_floor); - SolveErrors e{0, 0}; - for (int v = 0; v < Xhat.count(); ++v) { - e.fwd = std::max(e.fwd, max_abs_diff(Xhat[v], X.data, sB) / nX); - matmul(A.view(v), Xhat.view(v), AX); - e.res = std::max(e.res, max_abs_diff(AXs.data(), B[v], sB) / - std::max(norm1(B.view(v)), norm_floor)); - } - return e; -} - // One row of a C API validation table: the call made, what it returned, and // what LAPACK-style argument checking must return. struct ApiCheck { diff --git a/tests/test_gels_compact.cpp b/tests/test_gels_compact.cpp index 7dc3d0d..bbe9933 100644 --- a/tests/test_gels_compact.cpp +++ b/tests/test_gels_compact.cpp @@ -1,16 +1,18 @@ // test_gels_compact.cpp // -// Self-contained validation of the templated compact least-squares / -// minimum-norm solve (dgels_compact / sgels_compact), with no BLAS dependency. -// The reference is ref_gels: the unblocked algorithm (geqr2 of the tall -// orientation, orm2r, back substitution) in scalar form -- the same steps the -// vectorized kernel executes V lanes at a time, so a correct kernel matches it -// to working precision -- and, independently of it, the properties that define -// the two solutions. +// Validation of the templated compact least-squares / minimum-norm solve +// (dgels_compact / sgels_compact) against LAPACKE (test_lapack_util.hpp). +// The reference is ref_gels: LAPACKE_?gels for X and LAPACKE_?geqrf / +// ?gelqf for the factorization and tau (which is also what ?gels runs +// inside) -- the same steps the vectorized kernel executes V lanes at a +// time, blocked by LAPACK above its crossover -- and, independently of it, +// the properties that define the two solutions. // -// Checks per (T, V, layout, trans, shape): -// 1. X (all max(m,n) rows of B, the residual rows included) == ref_gels -// 2. the factorization left in A and the tau left in taup == ref_gels +// Checks per (T, V, layout, trans, shape), the design's 7.1 gates: +// 1. X (all max(m,n) rows of B, the residual rows included) == ref_gels, +// relative to ||X_lapack||_1 at 100 max(m,n) eps +// 2. the factorization left in A and the tau left in taup == ref_gels, +// relative to ||A||_1 at 100 max(m,n) eps // 3. least squares: the normal equations op(A)^T (B - op(A) X) = 0, and the // residual sums of squares in rows n..m-1 of B equal ||B - op(A) X||^2 // minimum norm: op(A) X = B, and X equals the minimum-norm solution formed @@ -27,7 +29,7 @@ #include #include -#include "test_compact_util.hpp" // compact, scalar references, MatrixBatch, pack/unpack +#include "test_lapack_util.hpp" // compact, ref_gels/geqr2/orm2r/trsm_upper, pack/unpack using namespace cbk::test; @@ -74,13 +76,15 @@ int run_case(char layout, char trans, int nm, int m, int n, int nrhs) unpack_compact(Bout, bp.data(), ldb, V, row); unpack_tau(tau_out, tp.data(), V); - // check 1 & 2: X (all p rows), the factorization, and tau vs the reference - double e_x = 0, e_h = 0, e_t = 0, nrm_b = 0; + // check 1 & 2: X (all p rows), the factorization, and tau vs the + // reference, each relative to the reference operand's L1 norm + double e_x = 0, e_h = 0, e_t = 0; for (int idx = 0; idx < nm; ++idx) { - e_x = std::max(e_x, max_abs_diff(Bout[idx], Bref[idx], B.stride())); - e_h = std::max(e_h, max_abs_diff(Aout[idx], Aref[idx], A.stride())); - e_t = std::max(e_t, max_abs_diff(tau_out[idx], tau_ref[idx], (size_t)q)); - nrm_b = std::max(nrm_b, norm1(Bref.view(idx))); + const double na = std::max(norm1(A.view(idx)), norm_floor); + e_x = std::max(e_x, max_abs_diff(Bout[idx], Bref[idx], B.stride()) / + std::max(norm1(Bref.view(idx)), norm_floor)); + e_h = std::max(e_h, max_abs_diff(Aout[idx], Aref[idx], A.stride()) / na); + e_t = std::max(e_t, max_abs_diff(tau_out[idx], tau_ref[idx], (size_t)q) / na); } // check 3: the defining properties, formed without the reference: e_prop is @@ -128,7 +132,7 @@ int run_case(char layout, char trans, int nm, int m, int n, int nrhs) std::max(e_prop, std::abs((double)R(i, j) - Bin(i, j)) / scale); // minimum norm the other way: X = op(A)^T Z, G Z = B with the Gram // matrix G = op(A) op(A)^T (rows_op x rows_op, SPD), solved by the - // scalar QR references + // LAPACKE QR references std::vector Gs((size_t)rows_op * rows_op), Zs((size_t)rows_op * nrhs), tg(rows_op), Xmns((size_t)cols_op * nrhs); const auto G = mat_view(Gs.data(), rows_op, rows_op); @@ -147,8 +151,11 @@ int run_case(char layout, char trans, int nm, int m, int n, int nrhs) } } - const double tol_ref = 200.0 * eps * p; // same op sequence, ~eps - const double tol_x = 200.0 * eps * p * std::max(nrm_b, 1.0); + // vs LAPACK: two backward-stable implementations (LAPACK's blocked above + // its crossover) of a well-conditioned op(A)'s factorization and solve, + // relative to the operand norms (design 7.1) + const double tol_ref = 100.0 * eps * p; + const double tol_x = 100.0 * eps * p; const double tol_prop = 100.0 * eps * p; // backward-stable quantities // the Gram-formed minimum-norm solution carries cond(G) = cond(op(A))^2 const double tol_prop2 = overdet ? tol_prop : 1e4 * eps * p; diff --git a/tests/test_gels_mkl.cpp b/tests/test_gels_mkl.cpp index b754d4b..2b341cc 100644 --- a/tests/test_gels_mkl.cpp +++ b/tests/test_gels_mkl.cpp @@ -32,7 +32,7 @@ * Assisted-by: Claude:claude-fable-5 */ -#include "test_mkl_util.hpp" /* compat, mkl, lapack + the MKL-free helpers */ +#include "test_mkl_util.hpp" /* compat, mkl, lapack + the shared helpers */ #include #include diff --git a/tests/test_geqrf_compact.cpp b/tests/test_geqrf_compact.cpp index e8402f0..3135e37 100644 --- a/tests/test_geqrf_compact.cpp +++ b/tests/test_geqrf_compact.cpp @@ -1,16 +1,22 @@ // test_geqrf_compact.cpp // -// Self-contained validation of the templated compact QR factorization -// (dgeqrf_compact / sgeqrf_compact), with no BLAS dependency. The reference is -// the unblocked Householder QR (LAPACK dgeqr2 / dlarfg) implemented in scalar -// form -- the same algorithm the vectorized kernel executes V lanes at a time, -// so a correct kernel matches it to working precision. +// Validation of the templated compact QR factorization (dgeqrf_compact / +// sgeqrf_compact) against LAPACKE (test_lapack_util.hpp). The reference is +// LAPACK's unblocked Householder QR, ?geqr2 / ?larfg -- the same algorithm the +// vectorized kernel executes V lanes at a time, so a correct kernel matches it +// to working precision. // // Checks per (T, V): -// 1. compact (H, tau) == scalar geqr2 (H, tau) (elementwise, ~eps*scale) -// 2. reconstruction Q * triu(H) == A (valid factorization) -// 3. reflectors are usable: ormqr_compact('T') then triangular solve -// recovers a known X from B = A X (in-situ with ormqr) +// 1. compact (H, tau) == LAPACKE_?geqr2 (H, tau) (elementwise, ~eps*scale) +// 2. reconstruction Q * triu(H) == A by ?ormqr (valid factorization) +// 3. reflectors are usable: ormqr_compact('T') then ?trsm recovers a known +// X from B = A X (in-situ with ormqr) +// and, per (T, structure, cond) on the dense-LAPACK contract (design 7.1): +// 4. Q = ?orgqr(H, tau), R = triu(H): the factorization residual +// ||R - Q^T A||_1 / ||A||_1 <= 20 n eps and the orthogonality +// ||Q^T Q - I||_1 <= 100 n eps, on dense, column-scaled, rank-deficient +// and near-collinear inputs; (H, tau) elementwise vs the blocked +// LAPACKE_?geqrf on the dense ones (relative, 100 n eps) // plus LAPACK-style argument validation of the C API. // // Assisted-by: Claude:claude-opus-4-8 Claude:claude-fable-5 @@ -23,7 +29,7 @@ #include #include -#include "test_compact_util.hpp" // compact, scalar references, MatrixBatch, pack/unpack +#include "test_lapack_util.hpp" // compact, lapack, ref_*, MatrixBatch, pack/unpack using namespace cbk::test; @@ -94,9 +100,9 @@ template static int run_case(int nm, int m, int n) } const double scale = std::max(1, m); - const double tol_fac = 200.0 * eps * scale; // same op sequence, ~eps - const double tol_rec = 200.0 * eps * scale; - const double tol_sol = 1e5 * eps * m; // cond(A)-dependent + const double tol_fac = 200.0 * eps * scale; // same op sequence as ?geqr2, ~eps + const double tol_rec = 200.0 * eps * scale; // ?ormqr's rounding differs, same bound + const double tol_sol = 1e5 * eps * m; // cond(A)-dependent bool ok_h = e_h <= tol_fac, ok_t = e_t <= tol_fac, ok_r = e_rec <= tol_rec; bool ok_s = (e_solve < 0) || (e_solve <= tol_sol); @@ -108,6 +114,115 @@ template static int run_case(int nm, int m, int n) return (info != 0) + !ok_h + !ok_t + !ok_r + !ok_s; } +// ------------------ invariants vs dense LAPACK (design 7.1) ---------- + +// Input structures (a subset of the GPU competition's stress set). The QR +// residual and orthogonality are backward-stable quantities, so they must hold +// to working precision for every structure -- rank deficiency and near-collinear +// columns included -- exactly as they do for dense LAPACK. +enum Structure { DENSE, RANK_DEFICIENT, NEAR_COLLINEAR }; + +// A random m x n matrix. For DENSE, a diagonal boost tames the conditioning +// and cond applies the competition's column scaling (columns *= +// logspace(0,-cond,n)). The structured variants make the last column a +// (near-)copy of the first, so the trailing reflector sees a (near-)zero +// sub-diagonal norm -- the branch the masked larfg must get right. +template static void gen_matrix(MatrixView A, double cond, Structure s) +{ + const int m = A.rows, n = A.cols; + for (int j = 0; j < n; ++j) + for (int i = 0; i < m; ++i) + A(i, j) = frand(); + if (s == DENSE) { + for (int d = 0; d < std::min(m, n); ++d) + A(d, d) += T(2); + if (cond > 0.0) + for (int j = 0; j < n; ++j) { + double sc = std::pow(10.0, -cond * (n > 1 ? (double)j / (n - 1) : 0.0)); + for (int i = 0; i < m; ++i) + A(i, j) *= sc; + } + } + else if (n >= 2) { + const double noise = (s == NEAR_COLLINEAR) ? 1e-9 : 0.0; /* exact dup if 0 */ + for (int i = 0; i < m; ++i) + A(i, n - 1) = A(i, 0) * (1.0 + noise * frand()); + } +} + +template +static int run_invariants(int nm, int m, int n, double cond, Structure structure = DENSE) +{ + const double eps = std::numeric_limits::epsilon(); + const int k = std::min(m, n); + + MatrixBatch A(nm, m, n); + for (int idx = 0; idx < nm; ++idx) + gen_matrix(A.view(idx), cond, structure); + + std::vector ap = pack_compact(A, m, V); + std::vector tp = compact_buffer(nm, k, 1, k, V); + const int info = compact::geqrf('C', m, n, ap.data(), m, tp.data(), V, nm); + MatrixBatch H(nm, m, n), tau(nm, k, 1); + unpack_compact(H, ap.data(), m, V); + unpack_tau(tau, tp.data(), V); + + double worst_res = 0, worst_orth = 0, worst_el = 0; + std::vector Qs((size_t)m * k), QtAs((size_t)k * n), QtQs((size_t)k * k), + Hrefs(A.stride()), tauref(k); + const auto Q = mat_view(Qs.data(), m, k), QtA = mat_view(QtAs.data(), k, n), + QtQ = mat_view(QtQs.data(), k, k), Href = mat_view(Hrefs.data(), m, n); + for (int idx = 0; idx < nm; ++idx) { + const auto Am = A.view(idx), Hm = H.view(idx); + + // Q = the first k columns of H(0)..H(k-1), from the reflectors in the + // first k columns of H + copy_matrix(leading(Hm, m, k), Q); + ref_org2r(k, Q, tau[idx]); + + // residual R - Q^T A, R = triu(H) (k x n) + matmul(Q.transposed(), Am, QtA); + double resid = 0; + for (int j = 0; j < n; ++j) + for (int i = 0; i < k; ++i) + resid = std::max( + resid, (double)std::abs(((i <= j) ? Hm(i, j) : T(0)) - QtA(i, j))); + worst_res = std::max(worst_res, resid / std::max(norm1(Am), norm_floor)); + + // orthogonality Q^T Q - I + matmul(Q.transposed(), Q, QtQ); + for (int d = 0; d < k; ++d) + QtQ(d, d) -= T(1); + worst_orth = std::max(worst_orth, norm1(QtQ)); + + // elementwise vs the blocked LAPACKE_?geqrf + copy_matrix(Am, Href); + ref_geqrf(Href, tauref.data()); + const double el = std::max(max_abs_diff(H[idx], Hrefs.data(), A.stride()), + max_abs_diff(tau[idx], tauref.data(), (size_t)k)); + worst_el = std::max(worst_el, el / std::max(norm1(Am), norm_floor)); + } + + // The elementwise difference vs LAPACKE_?geqrf is gated on the dense + // inputs only, where the reflectors are essentially unique: it catches a + // sign-convention regression the residual gates cannot see (observed + // ~n*eps). Rank-deficient and near-collinear inputs stay ungated (their + // reflectors are not unique; el ~ 1e-2 is expected) and it is printed as a + // diagnostic there (design doc 7.1). + const double rtol_res = 20.0 * n * eps, rtol_orth = 100.0 * n * eps; + const double rtol_el = (structure == DENSE) ? 100.0 * n * eps : HUGE_VAL; + const bool ok = (info == 0) && (worst_res <= rtol_res) && (worst_orth <= rtol_orth) && + (worst_el <= rtol_el); + const char *sname = structure == DENSE ? "dense" + : structure == RANK_DEFICIENT ? "rankdef" + : "collin"; + std::printf("T=%-6s V=%-2d nm=%-2d m=%-3d n=%-3d cond=%.0f %-8s| res %.2e (%.1e) " + "orth %.2e (%.1e) el %.1e | info=%d %s\n", + compact::name, V, nm, m, n, cond, sname, worst_res, rtol_res, + worst_orth, rtol_orth, worst_el, info, ok ? "OK" : "FAIL"); + return !ok; +} + // ------------------------ underflow scope (design 6.6) --------------- // A column whose sub-diagonal entries are so small that their squares @@ -228,5 +343,25 @@ int main() fails += test_underflow(); fails += test_underflow(); + // the dense-LAPACK contract (design 7.1): residual and orthogonality of + // the materialized Q, over conditioning and structure, in both precisions + for (int prec = 0; prec < 2; ++prec) { + auto inv = [&](int nm, int m, int n, double cond, Structure st = DENSE) { + return prec ? run_invariants(nm, m, n, cond, st) + : run_invariants(nm, m, n, cond, st); + }; + fails += inv(8, 30, 30, 0.0); + fails += inv(16, 60, 60, 0.0); + fails += inv(11, 43, 43, 0.0); // padded partial group + fails += inv(8, 64, 20, 0.0); // tall + fails += inv(8, 20, 64, 0.0); // wide + fails += inv(8, 30, 30, 4.0); // column-scaled (dynamic range) + fails += inv(8, 128, 128, 0.0); + // conditioning-robustness stress: backward-stable gates must still hold + fails += inv(8, 40, 40, 0.0, RANK_DEFICIENT); + fails += inv(8, 40, 40, 0.0, NEAR_COLLINEAR); + fails += inv(11, 60, 24, 0.0, RANK_DEFICIENT); // wide-ish, padded group + } + return finish(fails); } diff --git a/tests/test_geqrf_mkl.cpp b/tests/test_geqrf_mkl.cpp index b8f812e..88271db 100644 --- a/tests/test_geqrf_mkl.cpp +++ b/tests/test_geqrf_mkl.cpp @@ -1,19 +1,11 @@ /* test_geqrf_mkl.cpp * - * Validation of cbk_?geqrf_compact against real Intel MKL and dense LAPACK, - * through the genuine MKL Compact pipeline (mkl_dgepack_compact / - * mkl_dgeunpack_compact). This is the design document's "hard correctness gate - * against standard dense LAPACK equivalents" (cbk_dgeqrf_compact_design.md - * section 7). - * - * Suite 1 (section 7.1) -- Factorization invariants vs dense LAPACK: - * A random batch is factored by cbk_dgeqrf_compact, unpacked, and per - * matrix checked against the LAPACK-style QR contract used by the GPU - * competition checker: Q is materialized from (H, tau) via LAPACKE_dorgqr, - * R = triu(H), and the routine gates - * - factorization residual || R - Q^T A ||_1 / ||A||_1 <= 20 * n * eps, - * - orthogonality || Q^T Q - I ||_1 <= 100 * n * eps. - * (H, tau) are also compared elementwise to LAPACKE_dgeqrf as a diagnostic. + * Validation of cbk_?geqrf_compact against real Intel MKL, through the + * genuine MKL Compact pipeline (mkl_dgepack_compact / mkl_dgeunpack_compact): + * the MKL-side half of cbk_dgeqrf_compact_design.md section 7. The + * dense-LAPACK invariants of section 7.1 (Q from ?orgqr, residual and + * orthogonality gates over the structured inputs) run in the portable suite, + * test_geqrf_compact.cpp, on any LAPACKE stack. * * Suite 2 (section 7.2) -- Cross-check vs mkl_dgeqrf_compact: * The same packed batch is factored by both cbk_dgeqrf_compact and the @@ -29,7 +21,7 @@ * Assisted-by: Claude:claude-opus-4-8 Claude:claude-fable-5 */ -#include "test_mkl_util.hpp" /* compat, mkl, lapack + the MKL-free helpers */ +#include "test_mkl_util.hpp" /* compat, mkl, lapack + the shared helpers */ #include #include @@ -42,144 +34,6 @@ using namespace cbk::test; namespace { -/* Input structures (a subset of the GPU competition's stress set). The QR - * residual and orthogonality are backward-stable quantities, so they must hold - * to working precision for every structure -- rank deficiency and near-collinear - * columns included -- exactly as they do for dense LAPACK. */ -enum Structure { DENSE, RANK_DEFICIENT, NEAR_COLLINEAR }; - -/* build a random m x n matrix. For DENSE, a diagonal boost tames the - * conditioning and cond applies the competition's column scaling - * (columns *= logspace(0,-cond,n)). The structured variants make the last column - * a (near-)copy of the first, so the trailing reflector sees a (near-)zero - * sub-diagonal norm -- the branch the masked larfg must get right. */ -template void gen_matrix(MatrixView A, double cond, Structure s = DENSE) -{ - const int m = A.rows, n = A.cols; - for (int j = 0; j < n; ++j) - for (int i = 0; i < m; ++i) - A(i, j) = frand(); - if (s == DENSE) { - for (int d = 0; d < std::min(m, n); ++d) - A(d, d) += T(2); - if (cond > 0.0) - for (int j = 0; j < n; ++j) { - double sc = std::pow(10.0, -cond * (n > 1 ? (double)j / (n - 1) : 0.0)); - for (int i = 0; i < m; ++i) - A(i, j) *= sc; - } - } - else if (n >= 2) { - const double noise = (s == NEAR_COLLINEAR) ? 1e-9 : 0.0; /* exact dup if 0 */ - for (int i = 0; i < m; ++i) - A(i, n - 1) = A(i, 0) * (1.0 + noise * frand()); - } -} - -/* ---------------- Suite 1: invariants vs dense LAPACK (col-major) ------ */ - -template -int suite1(int nm, int m, int n, double cond, Structure structure = DENSE) -{ - const double eps = std::numeric_limits::epsilon(); - const MKL_COMPACT_PACK fmt = mkl_get_format_compact(); - const int V = mkl::vlen(fmt); - const int k = std::min(m, n); - - MatrixBatch A(nm, m, n); - for (int v = 0; v < nm; ++v) - gen_matrix(A.view(v), cond, structure); - - /* pack A, factor with the routine under test, unpack (H, tau) */ - auto Ap = A.base_ptrs(); - MKL_INT sz_a = mkl::get_size(m, n, fmt, nm); - MKL_INT sz_t = mkl::get_size(k, 1, fmt, nm); - auto ap_buf = cbk::detail::mkl_alloc_bytes(sz_a); - auto taup_buf = cbk::detail::mkl_alloc_bytes(sz_t); - T *ap = ap_buf.get(), *taup = taup_buf.get(); - mkl::gepack(MKL_COL_MAJOR, m, n, Ap.data(), m, ap, m, fmt, nm); - - MKL_INT info = 99; - T wq; - compat::geqrf(MKL_COL_MAJOR, m, n, ap, m, taup, &wq, -1, &info, fmt, nm); - compat::geqrf(MKL_COL_MAJOR, m, n, ap, m, taup, &wq, (MKL_INT)wq, &info, fmt, nm); - - MatrixBatch H(nm, m, n), tau(nm, k, 1); - auto Hp = H.base_ptrs(); - auto Tp = tau.base_ptrs(); - mkl::geunpack(MKL_COL_MAJOR, m, n, Hp.data(), m, ap, m, fmt, nm); - mkl::geunpack(MKL_COL_MAJOR, k, 1, Tp.data(), k, taup, k, fmt, nm); - - int fails = 0; - if (info != 0) { - ++fails; - std::printf(" info = %ld (expected 0)\n", (long)info); - } - - double worst_res = 0, worst_orth = 0, worst_el = 0; - std::vector Q((size_t)m * k), QtA((size_t)k * n), QtQ((size_t)k * k); - std::vector Href(A.stride()), tauref(tau.stride()); - const auto QtAm = mat_view(QtA.data(), k, n); - const auto QtQm = mat_view(QtQ.data(), k, k); - for (int v = 0; v < nm; ++v) { - const T *Av = A[v]; - const T *Hv = H[v], *tv = tau[v]; - const auto Am = A.view(v); /* the input, column-major */ - const auto Hm = H.view(v); /* the factor (H, tau) it produced */ - - /* Q = householder_product(H, tau): first k columns of Q (m x k). The - * reflectors occupy the first k columns of the m x n H, i.e. the first - * m*k contiguous (column-major) elements. */ - std::copy(Hv, Hv + (size_t)m * k, Q.begin()); - lapack::orgqr(LAPACK_COL_MAJOR, m, k, k, Q.data(), m, tv); - - /* residual R - Q^T A, R = triu(H) (k x n) */ - lapack::gemm(CblasColMajor, CblasTrans, CblasNoTrans, k, n, m, T(1), Q.data(), - m, Av, m, T(0), QtA.data(), k); - double resid = 0; - for (int j = 0; j < n; ++j) - for (int i = 0; i < k; ++i) { - double R = (i <= j) ? Hm(i, j) : 0.0; /* triu */ - resid = std::max(resid, std::abs(R - QtAm(i, j))); - } - worst_res = std::max(worst_res, resid / std::max(norm1(Am), norm_floor)); - - /* orthogonality Q^T Q - I (k x k) */ - lapack::gemm(CblasColMajor, CblasTrans, CblasNoTrans, k, k, m, T(1), Q.data(), - m, Q.data(), m, T(0), QtQ.data(), k); - for (int d = 0; d < k; ++d) - QtQm(d, d) -= T(1); - worst_orth = std::max(worst_orth, norm1(QtQm)); - - /* diagnostic: elementwise vs LAPACKE_dgeqrf */ - std::copy(Av, Av + A.stride(), Href.begin()); - lapack::geqrf(LAPACK_COL_MAJOR, m, n, Href.data(), m, tauref.data()); - double el = std::max(max_abs_diff(Hv, Href.data(), A.stride()), - max_abs_diff(tv, tauref.data(), tau.stride())); - worst_el = std::max(worst_el, el / std::max(norm1(Am), norm_floor)); - } - - /* The elementwise difference vs LAPACKE_?geqrf is gated on the dense - * inputs only, where the reflectors are essentially unique: it catches a - * sign-convention regression the residual gates cannot see (observed - * ~n*eps). Rank-deficient and near-collinear inputs stay ungated (their - * reflectors are not unique; el ~ 1e-2 is expected) and it is printed as a - * diagnostic there (design doc 7.1). */ - const double rtol_res = 20.0 * n * eps, rtol_orth = 100.0 * n * eps; - const double rtol_el = (structure == DENSE) ? 100.0 * n * eps : HUGE_VAL; - bool ok = - (worst_res <= rtol_res) && (worst_orth <= rtol_orth) && (worst_el <= rtol_el); - fails += !ok; - const char *sname = structure == DENSE ? "dense" - : structure == RANK_DEFICIENT ? "rankdef" - : "collin"; - std::printf(" [suite1] V=%-2d nm=%-2d m=%-3d n=%-3d cond=%.0f %-8s| res %.2e (%.1e) " - "orth %.2e (%.1e) el %.1e %s\n", - V, nm, m, n, cond, sname, worst_res, rtol_res, worst_orth, rtol_orth, - worst_el, ok ? "OK" : "FAIL"); - return fails; -} - /* ---------------- Suite 2: cross-check vs mkl_dgeqrf_compact ----------- */ template int suite2(MKL_LAYOUT layout, int nm, int m, int n) @@ -192,7 +46,7 @@ template int suite2(MKL_LAYOUT layout, int nm, int m, int n) MatrixBatch A(nm, m, n); for (int v = 0; v < nm; ++v) - gen_matrix(A.view(v), 0.0); + gen_boosted(A.view(v)); /* A was generated column-major; for a row-major run reinterpret the same * numbers as a row-major m x n (a genuinely different matrix, still fine). */ auto Ap = A.base_ptrs(); @@ -252,7 +106,7 @@ template int suite3(int nm, int n, int nrhs) MatrixBatch A(nm, n, n), B(nm, n, nrhs); for (int v = 0; v < nm; ++v) { - gen_matrix(A.view(v), 0.0); + gen_boosted(A.view(v)); matmul(A.view(v), X, B.view(v)); /* B = A X */ } auto Ap = A.base_ptrs(); @@ -320,19 +174,6 @@ template int run_suites() int fails = 0; - /* Suite 1: invariants vs dense LAPACK, square + rectangular, a range of cond */ - fails += suite1(8, 30, 30, 0.0); - fails += suite1(16, 60, 60, 0.0); - fails += suite1(11, 43, 43, 0.0); /* padded partial group */ - fails += suite1(8, 64, 20, 0.0); /* tall */ - fails += suite1(8, 20, 64, 0.0); /* wide */ - fails += suite1(8, 30, 30, 4.0); /* column-scaled (dynamic range) */ - fails += suite1(8, 128, 128, 0.0); - /* conditioning-robustness stress: backward-stable gates must still hold */ - fails += suite1(8, 40, 40, 0.0, RANK_DEFICIENT); - fails += suite1(8, 40, 40, 0.0, NEAR_COLLINEAR); - fails += suite1(11, 60, 24, 0.0, RANK_DEFICIENT); /* wide-ish, padded group */ - /* Suite 2: cross-check vs mkl_dgeqrf_compact, both layouts */ fails += suite2(MKL_COL_MAJOR, 8, 30, 30); fails += suite2(MKL_COL_MAJOR, 8, 64, 20); diff --git a/tests/test_lapack_util.hpp b/tests/test_lapack_util.hpp new file mode 100644 index 0000000..24f619f --- /dev/null +++ b/tests/test_lapack_util.hpp @@ -0,0 +1,427 @@ +// test_lapack_util.hpp +// +// The dense references of every test suite, on a real LAPACKE + CBLAS stack +// (issue #27): lapack, a struct of static forwarders to the LAPACKE_?* and +// cblas_?* entry points specialized for double (d) and float (s), and over it +// the reference procedures the suites call by name -- ref_geqr2, ref_orm2r, +// ref_potf2, ref_trsm, ref_gels, ..., matmul, tri_apply, solve_errors -- each +// taking the project's MatrixView operands and forwarding their `.data` and +// leading dimensions to the library. The only reference the suites still +// hand-roll is the unpivoted LDL^T of test_sytrfnp_compact.cpp (LAPACK has no +// such routine), and that one is validated against this library too. +// +// Which headers to include is the stack's business, decided by +// cmake/FindLAPACKE.cmake through one definition on the LAPACKE::LAPACKE +// target: MKL's own / when the tests are built +// with the MKL extension (a test binary must never link two LAPACK stacks), +// the plain / of OpenBLAS or Netlib otherwise. +// Integer arguments are lapack_int (MKL_INT under MKL, so an ilp64 MKL build +// gets its 64-bit integers); the CBLAS calls take the same lapack_int and +// convert to the stack's own BLAS integer, which matches on every stack the +// tests support (the other stacks are used at their 32-bit interface). +// +// Assisted-by: Claude:claude-fable-5 + +#ifndef TEST_LAPACK_UTIL_HPP +#define TEST_LAPACK_UTIL_HPP + +#include "test_compact_util.hpp" + +#if defined(CBK_LAPACKE_MKL) +#include +#include +#elif defined(CBK_LAPACKE_ACCELERATE) +#include +#include +#else +#include +#include +#endif + +#include +#include +#include +#include +#include + +namespace cbk::test { + +// ----------------------- LAPACKE / CBLAS, by scalar type ------------- +// lapack::geqr2 / geqrf / gelqf / orgqr / ormqr / geqp3 / potrf / gels / +// larfg forward to LAPACKE_?*, gemm / trsm / trmm to cblas_?*. One macro +// generates both specializations so the two cannot drift apart. + +template struct lapack; + +// clang-format off +// NOLINTBEGIN(bugprone-macro-parentheses): T is a type name, p a token to paste +#define CBK_TEST_LAPACK_DISPATCH(T, p) \ +template <> struct lapack { \ + static lapack_int larfg(lapack_int n, T *alpha, T *x, lapack_int incx, T *tau) \ + { return LAPACKE_##p##larfg(n, alpha, x, incx, tau); } \ + static lapack_int geqr2(int layout, lapack_int m, lapack_int n, T *a, lapack_int lda, \ + T *tau) \ + { return LAPACKE_##p##geqr2(layout, m, n, a, lda, tau); } \ + static lapack_int geqrf(int layout, lapack_int m, lapack_int n, T *a, lapack_int lda, \ + T *tau) \ + { return LAPACKE_##p##geqrf(layout, m, n, a, lda, tau); } \ + static lapack_int gelqf(int layout, lapack_int m, lapack_int n, T *a, lapack_int lda, \ + T *tau) \ + { return LAPACKE_##p##gelqf(layout, m, n, a, lda, tau); } \ + static lapack_int orgqr(int layout, lapack_int m, lapack_int n, lapack_int k, T *a, \ + lapack_int lda, const T *tau) \ + { return LAPACKE_##p##orgqr(layout, m, n, k, a, lda, tau); } \ + static lapack_int ormqr(int layout, char side, char trans, lapack_int m, \ + lapack_int n, lapack_int k, const T *a, lapack_int lda, \ + const T *tau, T *c, lapack_int ldc) \ + { return LAPACKE_##p##ormqr(layout, side, trans, m, n, k, a, lda, tau, c, ldc); } \ + static lapack_int geqp3(int layout, lapack_int m, lapack_int n, T *a, lapack_int lda, \ + lapack_int *jpvt, T *tau) \ + { return LAPACKE_##p##geqp3(layout, m, n, a, lda, jpvt, tau); } \ + static lapack_int potrf(int layout, char uplo, lapack_int n, T *a, lapack_int lda) \ + { return LAPACKE_##p##potrf(layout, uplo, n, a, lda); } \ + static lapack_int gels(int layout, char trans, lapack_int m, lapack_int n, \ + lapack_int nrhs, T *a, lapack_int lda, T *b, lapack_int ldb) \ + { return LAPACKE_##p##gels(layout, trans, m, n, nrhs, a, lda, b, ldb); } \ + static void gemm(CBLAS_LAYOUT layout, CBLAS_TRANSPOSE ta, CBLAS_TRANSPOSE tb, \ + lapack_int m, lapack_int n, lapack_int k, T alpha, const T *a, \ + lapack_int lda, const T *b, lapack_int ldb, T beta, T *c, \ + lapack_int ldc) \ + { cblas_##p##gemm(layout, ta, tb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc); } \ + static void trsm(CBLAS_LAYOUT layout, CBLAS_SIDE side, CBLAS_UPLO uplo, \ + CBLAS_TRANSPOSE ta, CBLAS_DIAG diag, lapack_int m, lapack_int n, \ + T alpha, const T *a, lapack_int lda, T *b, lapack_int ldb) \ + { cblas_##p##trsm(layout, side, uplo, ta, diag, m, n, alpha, a, lda, b, ldb); } \ + static void trmm(CBLAS_LAYOUT layout, CBLAS_SIDE side, CBLAS_UPLO uplo, \ + CBLAS_TRANSPOSE ta, CBLAS_DIAG diag, lapack_int m, lapack_int n, \ + T alpha, const T *a, lapack_int lda, T *b, lapack_int ldb) \ + { cblas_##p##trmm(layout, side, uplo, ta, diag, m, n, alpha, a, lda, b, ldb); } \ +}; +// NOLINTEND(bugprone-macro-parentheses) +// clang-format on + +CBK_TEST_LAPACK_DISPATCH(double, d) +CBK_TEST_LAPACK_DISPATCH(float, s) +#undef CBK_TEST_LAPACK_DISPATCH + +// ----------------------- views as LAPACKE / CBLAS arguments ---------- +// A MatrixView carries its layout as strides (si, sj); LAPACKE and CBLAS take +// one layout flag per call plus a leading dimension per operand. A view maps +// without a copy when one stride is 1 and the other is at least the extent +// along it -- column-major (si == 1, ld = sj) or row-major (sj == 1, ld = si); +// the ambiguous degenerate shapes (a single row or column) are read as +// whichever description satisfies the library's ld >= extent rule. + +template bool is_rowmajor(const Mv &M) +{ + if (M.si == 1 && M.sj >= std::max(M.rows, 1)) return false; + assert(M.sj == 1 && M.si >= std::max(M.cols, 1) && "view is not a BLAS layout"); + return true; +} + +template lapack_int lapack_ld(const Mv &M) +{ + return std::max(is_rowmajor(M) ? M.si : M.sj, 1); +} + +template int lapack_layout(const Mv &M) +{ + return is_rowmajor(M) ? LAPACK_ROW_MAJOR : LAPACK_COL_MAJOR; +} + +template CBLAS_LAYOUT cblas_layout(const Mv &M) +{ + return is_rowmajor(M) ? CblasRowMajor : CblasColMajor; +} + +// The BLAS character arguments as CBLAS enumerators. +inline CBLAS_TRANSPOSE cblas_trans(char t) +{ + return (t == 'T' || t == 't' || t == 'C' || t == 'c') ? CblasTrans : CblasNoTrans; +} +inline CBLAS_TRANSPOSE cblas_flip(CBLAS_TRANSPOSE t) +{ + return t == CblasNoTrans ? CblasTrans : CblasNoTrans; +} +inline CBLAS_UPLO cblas_uplo(char u) +{ + return (u == 'U' || u == 'u') ? CblasUpper : CblasLower; +} +inline CBLAS_UPLO cblas_flip(CBLAS_UPLO u) +{ + return u == CblasUpper ? CblasLower : CblasUpper; +} +inline CBLAS_SIDE cblas_side(char s) +{ + return (s == 'L' || s == 'l') ? CblasLeft : CblasRight; +} +inline CBLAS_DIAG cblas_diag(char d) +{ + return (d == 'U' || d == 'u') ? CblasUnit : CblasNonUnit; +} + +// An operand of a call whose layout is fixed by another operand: the view +// itself when its layout already matches, else a contiguous staging copy in +// the call's layout -- copied back on destruction when the operand is +// writable (T non-const). Only the LAPACKE routines without a transpose flag +// for the operand need it (ormqr's reflectors, gels's right-hand side); the +// BLAS-3 calls read a mismatched operand as its transpose instead. +template class Staged { + public: + Staged(MatrixView M, bool rowmajor) : orig_(M), copied_(is_rowmajor(M) != rowmajor) + { + if (!copied_) { + use_ = M; + return; + } + buf_.resize((std::size_t)M.rows * M.cols); + const auto w = mat_view(buf_.data(), M.rows, M.cols, rowmajor ? M.cols : M.rows, + rowmajor); + copy_matrix(M, w); + use_ = MatrixView{w.data, w.si, w.sj, w.rows, w.cols}; + } + ~Staged() + { + if constexpr (!std::is_const_v) { + if (copied_) copy_matrix(use_, orig_); + } + } + Staged(const Staged &) = delete; + Staged &operator=(const Staged &) = delete; + + T *data() const { return use_.data; } + lapack_int ld() const { return lapack_ld(use_); } + + private: + MatrixView orig_; + bool copied_; + std::vector> buf_; + MatrixView use_{}; +}; + +// A read-only view of the same storage, whatever the constness of the view +// handed in (distinct MatrixView instantiations do not convert). +template ConstMatrixView> as_const_view(const Mv &M) +{ + return {M.data, M.si, M.sj, M.rows, M.cols}; +} + +// ----------------------- reference procedures ------------------------ +// The names the suites call, each forwarded to the library. LAPACK's blocked +// drivers (?geqrf, ?ormqr, ?orgqr, ?potrf, ?gels) run their unblocked +// counterparts below the crossover order their ilaenv settings pick (which +// covers most of the suites' sizes) and a blocked variant above it, so a +// comparison against them is a comparison of two backward-stable +// implementations of the same factorization, not of the same operation +// sequence: the suites gate such comparisons relative to the operand norms, +// at a multiple of n * eps. Only ?geqr2 and ?larfg are unblocked by +// definition, and only those two are compared elementwise at ~eps. + +// ?larfg: reflector from (alpha, x[0..m-2]); alpha := beta on exit. +template void ref_larfg(int m, T *alpha, T *x, T *tau) +{ + lapack::larfg(m, alpha, x, 1, tau); +} + +// ?geqr2: unblocked Householder QR, (H, tau) in the LAPACK convention. +template void ref_geqr2(MatrixView A, T *tau) +{ + lapack::geqr2(lapack_layout(A), A.rows, A.cols, A.data, lapack_ld(A), tau); +} + +// ?geqrf: the blocked driver, for the invariants suites. +template void ref_geqrf(MatrixView A, T *tau) +{ + lapack::geqrf(lapack_layout(A), A.rows, A.cols, A.data, lapack_ld(A), tau); +} + +// ?ormqr, side='L': B := Q^T B (trans 'T') or Q B ('N') from (H, tau), with +// k reflectors read from A (m x k or wider) and applied to B (m x nrhs). +template +void ref_orm2r(char trans, int k, Av A, const T *tau, MatrixView B) +{ + assert(A.rows == B.rows && k <= std::min(A.rows, A.cols)); + const bool row = is_rowmajor(B); + const Staged> As(as_const_view(A), row); + lapack::ormqr(lapack_layout(B), 'L', trans, B.rows, B.cols, k, As.data(), As.ld(), + tau, B.data, lapack_ld(B)); +} + +// ?orgqr: generate the first n columns of Q = H(0)..H(k-1) in place over the +// reflectors in columns 0..k-1 of A (m x n, m >= n >= k). +template void ref_org2r(int k, MatrixView A, const T *tau) +{ + assert(k <= A.cols && A.cols <= A.rows); + lapack::orgqr(lapack_layout(A), A.rows, A.cols, k, A.data, lapack_ld(A), tau); +} + +// ?geqp3: column-pivoted QR. On exit A holds the reflectors below the diagonal +// and R on/above it for the *permuted* matrix A(:,jpvt); jpvt[j] is the +// ORIGINAL (0-based) column index placed at position j, so A(:,jpvt) = Q R. +// (LAPACK's jpvt is 1-based and takes "free column" markers on input; both +// are handled here.) +template void ref_geqp3(MatrixView A, lapack_int *jpvt, T *tau) +{ + std::fill(jpvt, jpvt + A.cols, 0); + lapack::geqp3(lapack_layout(A), A.rows, A.cols, A.data, lapack_ld(A), jpvt, tau); + for (int j = 0; j < A.cols; ++j) + jpvt[j] -= 1; +} + +// ?potrf on a dense n x n matrix, in place. Lower: A = L L^T, factor in the +// lower triangle. Upper: A = U^T U, factor in the upper triangle. Only the +// named triangle is read or written. Returns LAPACK's info (> 0: not SPD). +template lapack_int ref_potf2(char uplo, MatrixView A) +{ + assert(A.rows == A.cols); + return lapack::potrf(lapack_layout(A), uplo, A.rows, A.data, lapack_ld(A)); +} + +// The triangular operand of a BLAS-3 call in the layout of the call: an A +// stored the other way round is read as its transpose, with the transpose +// flag and the triangle flipped to compensate. +template +void tri_operand(bool rowmajor, Av &A, CBLAS_TRANSPOSE &ta, CBLAS_UPLO &uplo) +{ + if (is_rowmajor(A) != rowmajor) { + A = A.transposed(); + ta = cblas_flip(ta); + uplo = cblas_flip(uplo); + } +} + +// BLAS ?trsm: solves op(A) X = alpha B (side='L') or X op(A) = alpha B +// (side='R') in place, A the order-s triangular factor; alpha == 0 gives +// B := 0. Only the referenced triangle of A is touched; the diagonal is not +// referenced when diag='U'. The ?trsm suite's reference, and the triangular +// steps of ref_gels. +template +void ref_trsm(char side, char uplo, char transa, char diag, T alpha, Av A, + MatrixView B) +{ + const bool left = (side == 'L' || side == 'l'); + assert(A.rows == A.cols && A.rows == (left ? B.rows : B.cols)); + const bool row = is_rowmajor(B); + CBLAS_TRANSPOSE ta = cblas_trans(transa); + CBLAS_UPLO ul = cblas_uplo(uplo); + tri_operand(row, A, ta, ul); + lapack::trsm(cblas_layout(B), cblas_side(side), ul, ta, cblas_diag(diag), B.rows, + B.cols, alpha, A.data, lapack_ld(A), B.data, lapack_ld(B)); +} + +// Back substitution R X = B with R the upper triangle of an n x n array. +template void ref_trsm_upper(Rv R, MatrixView B) +{ + ref_trsm('L', 'U', 'N', 'N', T(1), R, B); +} + +// ?gels: the least-squares (op(A) with more rows than columns) or minimum-norm +// (more columns than rows) solution of op(A) X = B, op(A) = A ('N') or A^T +// ('T'), A m x n and B max(m,n) x nrhs, each in the layout its view carries. +// On exit B holds X (least squares: rows n..m-1 keep the residual), A the +// factorization of the tall orientation -- the QR of A (m >= n) or its LQ in +// ?gelqf storage (m < n) -- and tau its min(m,n) reflector scalars. LAPACK's +// ?gels leaves that factorization in A but keeps tau in its workspace, so +// (A, tau) are taken from ?geqrf / ?gelqf on a copy of the input, which is +// the same call ?gels makes. +template void ref_gels(char trans, MatrixView A, MatrixView B, T *tau) +{ + const int m = A.rows, n = A.cols, nrhs = B.cols; + assert(B.rows == std::max(m, n)); + const bool row = is_rowmajor(A); + + std::vector A0s((std::size_t)m * n); + const auto A0 = mat_view(A0s.data(), m, n, row ? n : m, row); + copy_matrix(A, A0); + + { + const Staged Bs(B, row); + lapack::gels(lapack_layout(A), trans, m, n, nrhs, A.data, lapack_ld(A), + Bs.data(), Bs.ld()); + } + if (m >= n) + lapack::geqrf(lapack_layout(A0), m, n, A0.data, lapack_ld(A0), tau); + else + lapack::gelqf(lapack_layout(A0), m, n, A0.data, lapack_ld(A0), tau); + copy_matrix(A0, A); +} + +// C (m x n) := A (m x k) * B (k x n) by ?gemm, for forming right-hand sides +// and residuals. The shapes come from the views, which also carry the +// layout: the call runs in C's, and an operand stored the other way round is +// passed as its transpose. +template void matmul(Av A, Bv B, Cv C) +{ + using T = elem_t; + assert(A.rows == C.rows && B.cols == C.cols && A.cols == B.rows); + const int k = A.cols; // the inner extent, before any operand is transposed + const bool row = is_rowmajor(C); + CBLAS_TRANSPOSE ta = CblasNoTrans, tb = CblasNoTrans; + if (is_rowmajor(A) != row) { + A = A.transposed(); + ta = CblasTrans; + } + if (is_rowmajor(B) != row) { + B = B.transposed(); + tb = CblasTrans; + } + lapack::gemm(cblas_layout(C), ta, tb, C.rows, C.cols, k, T(1), A.data, + lapack_ld(A), B.data, lapack_ld(B), T(0), C.data, lapack_ld(C)); +} + +// Apply a triangular operator to a general matrix by ?trmm -- the "forward" +// direction of a ?trsm, for checking a solve's defining residual +// ||op(A) X - alpha B||. R (m x n) := op(A) X (side 'L') or X op(A) (side +// 'R'), with A the order-s (s = m for 'L', n for 'R') triangular factor: +// uplo 'U'/'L', op(A) = A ('N') or A^T ('T'/'C'), unit ('U') or non-unit +// ('N') diagonal. Each operand carries its own layout in its view. +template +void tri_apply(char side, char uplo, char transa, char diag, Av A, Xv X, Rv R) +{ + using T = elem_t; + const bool left = (side == 'L' || side == 'l'); + const int s = left ? R.rows : R.cols; + assert(A.rows == s && A.cols == s && X.rows == R.rows && X.cols == R.cols); + (void)s; + copy_matrix(X, R); + const bool row = is_rowmajor(R); + CBLAS_TRANSPOSE ta = cblas_trans(transa); + CBLAS_UPLO ul = cblas_uplo(uplo); + tri_operand(row, A, ta, ul); + lapack::trmm(cblas_layout(R), cblas_side(side), ul, ta, cblas_diag(diag), R.rows, + R.cols, T(1), A.data, lapack_ld(A), R.data, lapack_ld(R)); +} + +// ----------------------- shared checks ------------------------------- + +// The end-to-end solve gate every solving suite closes with: the worst +// relative forward error of Xhat against the known X, and the worst relative +// residual ||A Xhat - B|| formed by ?gemm -- so a solver bug cannot hide +// behind the factorization that produced Xhat. +struct SolveErrors { + double fwd, res; +}; + +template +SolveErrors solve_errors(const MatrixBatch &A, const MatrixBatch &B, + const MatrixBatch &Xhat, ConstMatrixView X) +{ + assert(B.rows() == Xhat.rows() && B.cols() == Xhat.cols() && X.rows == Xhat.rows() && + X.cols == Xhat.cols()); + const size_t sB = Xhat.stride(); + std::vector AXs(sB); + const auto AX = mat_view(AXs.data(), Xhat.rows(), Xhat.cols()); + const double nX = std::max(norm1(X), norm_floor); + SolveErrors e{0, 0}; + for (int v = 0; v < Xhat.count(); ++v) { + e.fwd = std::max(e.fwd, max_abs_diff(Xhat[v], X.data, sB) / nX); + matmul(A.view(v), Xhat.view(v), AX); + e.res = std::max(e.res, max_abs_diff(AXs.data(), B[v], sB) / + std::max(norm1(B.view(v)), norm_floor)); + } + return e; +} + +} // namespace cbk::test + +#endif // TEST_LAPACK_UTIL_HPP diff --git a/tests/test_mkl_util.hpp b/tests/test_mkl_util.hpp index 2f0c6c9..b3ce0ec 100644 --- a/tests/test_mkl_util.hpp +++ b/tests/test_mkl_util.hpp @@ -8,12 +8,13 @@ // sysvnp,trsm,gels}_compact // mkl MKL's Compact API: sizes, pack/unpack, and MKL's own compact // kernels as references -// lapack dense LAPACKE / CBLAS references // // Each is a struct of static forwarders specialized for double (d/D) and float // (s/S); the specializations are generated by one macro per struct so the two -// cannot drift apart. Needs MKL headers; the MKL-free helpers (including -// compact, the portable C API's dispatch) are in test_compact_util.hpp. +// cannot drift apart. Needs MKL headers. The dense LAPACKE / CBLAS references +// (lapack, from MKL's own LAPACKE in this build) are in +// test_lapack_util.hpp, and the library-free helpers (including compact, +// the portable C API's dispatch) in test_compact_util.hpp; both come in here. // // Assisted-by: Claude:claude-opus-4-8 Claude:claude-fable-5 @@ -25,13 +26,12 @@ #include "cbk_compat.h" #include "cbk_mkl_alloc.h" -#include "test_compact_util.hpp" +#include "test_lapack_util.hpp" namespace cbk::test { template struct compat; template struct mkl; -template struct lapack; // clang-format off // NOLINTBEGIN(bugprone-macro-parentheses): T is a type name, p a token to paste @@ -115,30 +115,6 @@ template <> struct mkl { MKL_INT ldap, T *bp, MKL_INT ldbp, MKL_COMPACT_PACK fmt, MKL_INT nm) \ { mkl_##p##trsm_compact(layout, side, uplo, tr, diag, m, n, alpha, ap, ldap, bp, ldbp, \ fmt, nm); } \ -}; \ -template <> struct lapack { \ - static lapack_int geqrf(int layout, lapack_int m, lapack_int n, T *a, lapack_int lda, \ - T *tau) \ - { return LAPACKE_##p##geqrf(layout, m, n, a, lda, tau); } \ - static lapack_int orgqr(int layout, lapack_int m, lapack_int n, lapack_int k, T *a, \ - lapack_int lda, const T *tau) \ - { return LAPACKE_##p##orgqr(layout, m, n, k, a, lda, tau); } \ - static lapack_int ormqr(int layout, char side, char trans, lapack_int m, \ - lapack_int n, lapack_int k, const T *a, lapack_int lda, \ - const T *tau, T *c, lapack_int ldc) \ - { return LAPACKE_##p##ormqr(layout, side, trans, m, n, k, a, lda, tau, c, ldc); } \ - static lapack_int potrf(int layout, char uplo, lapack_int n, T *a, lapack_int lda) \ - { return LAPACKE_##p##potrf(layout, uplo, n, a, lda); } \ - static lapack_int gelqf(int layout, lapack_int m, lapack_int n, T *a, lapack_int lda, \ - T *tau) \ - { return LAPACKE_##p##gelqf(layout, m, n, a, lda, tau); } \ - static lapack_int gels(int layout, char trans, lapack_int m, lapack_int n, \ - lapack_int nrhs, T *a, lapack_int lda, T *b, lapack_int ldb) \ - { return LAPACKE_##p##gels(layout, trans, m, n, nrhs, a, lda, b, ldb); } \ - static void gemm(CBLAS_LAYOUT layout, CBLAS_TRANSPOSE ta, CBLAS_TRANSPOSE tb, \ - MKL_INT m, MKL_INT n, MKL_INT k, T alpha, const T *a, MKL_INT lda, \ - const T *b, MKL_INT ldb, T beta, T *c, MKL_INT ldc) \ - { cblas_##p##gemm(layout, ta, tb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc); } \ }; // NOLINTEND(bugprone-macro-parentheses) // clang-format on diff --git a/tests/test_orgqr_compact.cpp b/tests/test_orgqr_compact.cpp index 726b572..e0eb9ae 100644 --- a/tests/test_orgqr_compact.cpp +++ b/tests/test_orgqr_compact.cpp @@ -1,15 +1,19 @@ /* test_orgqr_compact.cpp * - * Self-contained validation of the templated compact orgqr. - * Reference: unblocked Householder QR (ref_geqr2, LAPACK reflector convention) - * + a scalar dorg2r, both templated on T. + * Validation of the templated compact orgqr against LAPACKE + * (test_lapack_util.hpp): the reflectors come from LAPACKE_?geqr2 and the + * reference Q from LAPACKE_?orgqr. * * Test design: factor a random m x k batch (diagonal-boosted), stage the * reflectors into the first k columns of an m x n buffer whose remaining * columns hold garbage (they must be overwritten), and generate Q. * * Checks per (T, V, layout): - * 1. compact Q == scalar dorg2r Q (elementwise, ~eps: same op sequence) + * 1. compact Q == LAPACKE_?orgqr Q (elementwise, 100 m eps: Q's entries + * are O(1) and, for full-rank A, Q is unique up to the shared sign + * convention, so two backward-stable routes agree to cond(A) * eps; + * LAPACK blocks its accumulation above k = 128, so the rounding order + * differs) * 2. Q^T Q = I (orthonormal columns, formed independently) * 3. Q(:, 0:k-1) R = A with R = triu of the factorization (k > 0) * 4. padded lanes of a partial last group come out exactly identity @@ -27,42 +31,10 @@ #include #include -#include "test_compact_util.hpp" // compact, scalar references, MatrixBatch, pack/unpack +#include "test_lapack_util.hpp" // compact, ref_geqr2, ref_org2r, MatrixBatch, pack/unpack using namespace cbk::test; -/* ----------------------- reference kernel (scalar) ------------------ */ - -/* dorg2r: generate the first n columns of Q = H(0)..H(k-1) in place over the - * reflectors of ref_geqr2 (columns 0..k-1 of A, m x n, m >= n >= k). Backward - * accumulation: unit-seed columns k..n-1, then apply each H(kk) to the - * trailing columns and form column kk = H(kk) e_kk. */ -template static void ref_org2r(int k, MatrixView A, const T *tau) -{ - const int m = A.rows, n = A.cols; - assert(k <= n && n <= m); - for (int j = k; j < n; ++j) { - for (int i = 0; i < m; ++i) - A(i, j) = T(0); - A(j, j) = T(1); - } - for (int kk = k - 1; kk >= 0; --kk) { - for (int j = kk + 1; j < n; ++j) { - T w = A(kk, j); - for (int i = kk + 1; i < m; ++i) - w += A(i, kk) * A(i, j); - A(kk, j) -= tau[kk] * w; - for (int i = kk + 1; i < m; ++i) - A(i, j) -= tau[kk] * A(i, kk) * w; - } - for (int i = kk + 1; i < m; ++i) - A(i, kk) *= -tau[kk]; - A(kk, kk) = T(1) - tau[kk]; - for (int i = 0; i < kk; ++i) - A(i, kk) = T(0); - } -} - /* --------------------------- one test case -------------------------- */ template @@ -70,7 +42,7 @@ static int run_case(int nm, int m, int n, int k, bool rowmajor = false) { assert(m >= n && n >= k); const T eps = std::numeric_limits::epsilon(); - const double tol_exact = 100.0 * eps; /* same op sequence as the reference */ + const double tol_el = 100.0 * eps * m; /* vs LAPACKE_?orgqr, see the header */ const double tol_orth = 100.0 * eps * m; /* independent invariants */ const double tol_rec = 100.0 * eps * m; @@ -106,7 +78,7 @@ static int run_case(int nm, int m, int n, int k, bool rowmajor = false) } unpack_compact(Qout, ap.data(), ld, V, rowmajor); - /* check 1: elementwise vs scalar dorg2r */ + /* check 1: elementwise vs LAPACKE_?orgqr */ double e1 = 0; for (int kk = 0; kk < nm; ++kk) e1 = std::max(e1, max_abs_diff(Qout[kk], Qref[kk], Qout.stride())); @@ -146,7 +118,7 @@ static int run_case(int nm, int m, int n, int k, bool rowmajor = false) (double)std::abs(P(i, j)[v] - (i == j ? T(1) : T(0)))); } - bool ok1 = e1 <= tol_exact, ok2 = e2 <= tol_orth, ok3 = e3 <= tol_rec, ok4 = e4 == 0; + bool ok1 = e1 <= tol_el, ok2 = e2 <= tol_orth, ok3 = e3 <= tol_rec, ok4 = e4 == 0; fails += !ok1 + !ok2 + !ok3 + !ok4; std::printf("T=%-6s V=%-2d nm=%-2d %s m=%-3d n=%-3d k=%-3d | ref: %.2e %s | QtQ=I: " "%.2e %s | QR=A: %.2e %s | pad: %.1e %s\n", @@ -203,6 +175,7 @@ int main() fails += run_case(8, 43, 43, 43); fails += run_case(16, 43, 43, 43); fails += run_case(32, 43, 43, 43); + fails += run_case(8, 128, 128, 128); /* LAPACK's blocked orgqr */ /* tall thin Q (n = k < m) */ fails += run_case(8, 43, 20, 20); diff --git a/tests/test_orgqr_mkl.cpp b/tests/test_orgqr_mkl.cpp index 4ac789a..9c3e435 100644 --- a/tests/test_orgqr_mkl.cpp +++ b/tests/test_orgqr_mkl.cpp @@ -29,7 +29,7 @@ * Assisted-by: Claude:claude-fable-5 */ -#include "test_mkl_util.hpp" /* compat, mkl, lapack + the MKL-free helpers */ +#include "test_mkl_util.hpp" /* compat, mkl, lapack + the shared helpers */ #include #include diff --git a/tests/test_ormqr_compact.cpp b/tests/test_ormqr_compact.cpp index d62eba1..38d2d7b 100644 --- a/tests/test_ormqr_compact.cpp +++ b/tests/test_ormqr_compact.cpp @@ -1,14 +1,17 @@ /* test_ormqr_compact.cpp * - * Self-contained validation of the templated compact ormqr. - * Reference: unblocked Householder QR (dgeqr2-style, LAPACK reflector - * convention) + scalar dorm2r + back substitution, all templated on T. + * Validation of the templated compact ormqr against LAPACKE + * (test_lapack_util.hpp): the reflectors come from LAPACKE_?geqr2 (or + * LAPACKE_?geqp3 for the pivoted pipeline), the reference application from + * LAPACKE_?ormqr, the back substitution from cblas_?trsm. * * Test design: X(:,j) = j+1 (ones, twos, threes, ...), B = A*X, so the * b columns are scaled row sums of A; the solve must recover X. * * Checks per (T, V): - * 1. compact Q^T B == scalar dorm2r Q^T B (elementwise, ~eps) + * 1. compact Q^T B == LAPACKE_?ormqr Q^T B, relative to ||Q^T B||_1 at + * 20 m eps (the design's 7.1 gate: two backward-stable applications of + * the same reflectors, LAPACK's blocked above its crossover) * 2. back substitution recovers X * 3. applying 'N' after 'T' recovers the original B (Q Q^T = I) * @@ -24,62 +27,18 @@ #include #include -#include "test_compact_util.hpp" // compact, scalar references, MatrixBatch, pack/unpack +#include "test_lapack_util.hpp" // compact, ref_geqr2/orm2r/geqp3/trsm_upper, pack/unpack using namespace cbk::test; -/* ----------------------- reference kernels (scalar) ----------------- */ -/* ref_larfg / ref_geqr2 / ref_orm2r / ref_trsm_upper come from test_compact_util.hpp. */ - -/* Column-pivoted Householder QR (dgeqp3-style, greedy max trailing-column - * norm). On exit A holds the reflectors below the diagonal and R on/above it - * for the *permuted* matrix A(:,jpvt); jpvt[j] is the ORIGINAL (0-based) - * column index placed at position j, so A(:,jpvt) = Q R. Used to drive the - * kernel from a rank-revealing factorization + back-permutation solve -- - * the production (RBF-FD) use case the kernel must support. */ -template static void ref_geqp3(MatrixView A, int *jpvt, T *tau) -{ - const int m = A.rows, n = A.cols, k = std::min(m, n); - for (int j = 0; j < n; ++j) - jpvt[j] = j; - - for (int kk = 0; kk < k; ++kk) { - int piv = kk; - T best = -1; - for (int j = kk; j < n; ++j) { - T s = 0; - for (int i = kk; i < m; ++i) - s += A(i, j) * A(i, j); - if (s > best) { - best = s; - piv = j; - } - } - if (piv != kk) { - for (int i = 0; i < m; ++i) - std::swap(A(i, kk), A(i, piv)); - std::swap(jpvt[kk], jpvt[piv]); - } - /* as in ref_geqr2: the reflector is the contiguous column tail */ - ref_larfg(m - kk, &A(kk, kk), A.col(kk) + kk + 1, &tau[kk]); - for (int j = kk + 1; j < n; ++j) { - T w = A(kk, j); - for (int i = kk + 1; i < m; ++i) - w += A(i, kk) * A(i, j); - A(kk, j) -= tau[kk] * w; - for (int i = kk + 1; i < m; ++i) - A(i, j) -= tau[kk] * A(i, kk) * w; - } - } -} - /* --------------------------- one test case -------------------------- */ template static int run_case(int nm, int m, int nrhs) { const int k = m; const T eps = std::numeric_limits::epsilon(); - const double tol_exact = 100.0 * eps; /* same op sequence */ + const double tol_el = 20.0 * eps * m; /* vs LAPACKE_?ormqr, relative */ + const double tol_exact = 100.0 * eps; /* the 'N' after 'T' round trip */ const double tol_solve = 1e5 * eps * m; /* cond(A)-dependent */ const std::vector Xs = known_solution(m, nrhs); @@ -101,12 +60,13 @@ template static int run_case(int nm, int m, int nrhs) std::vector tp = pack_tau(tau, V); std::vector bp = pack_compact(B, m, V); - /* check 1: compact Q^T B vs scalar */ + /* check 1: compact Q^T B vs LAPACKE_?ormqr, relative to ||Q^T B||_1 */ compact::ormqr('T', m, nrhs, k, ap.data(), m, tp.data(), bp.data(), m, V, nm); unpack_compact(Bout, bp.data(), m, V); double e1 = 0; for (int kk = 0; kk < nm; ++kk) - e1 = std::max(e1, max_abs_diff(Bout[kk], Bref[kk], (size_t)m * nrhs)); + e1 = std::max(e1, max_abs_diff(Bout[kk], Bref[kk], (size_t)m * nrhs) / + std::max(norm1(Bref.view(kk)), norm_floor)); /* check 2: solve recovers X */ double e2 = 0; @@ -123,7 +83,7 @@ template static int run_case(int nm, int m, int nrhs) e3 = std::max(e3, max_abs_diff(Bout[kk], B[kk], (size_t)m * nrhs)); /* scale-aware: B entries are O(m), QQ^t roundtrip accumulates a bit */ - bool ok1 = e1 <= tol_exact * m, ok2 = e2 <= tol_solve, ok3 = e3 <= tol_exact * m * 10; + bool ok1 = e1 <= tol_el, ok2 = e2 <= tol_solve, ok3 = e3 <= tol_exact * m * 10; std::printf("T=%-6s V=%-2d nm=%-2d m=%-3d nrhs=%d | QtB: %.2e %s | solve X: %.2e %s " "| QQt=I: %.2e %s\n", compact::name, V, nm, m, nrhs, e1, ok1 ? "OK" : "FAIL", e2, @@ -133,8 +93,9 @@ template static int run_case(int nm, int m, int nrhs) /* ------------------- pivoted-QR + back-permutation solve ------------ */ /* Salvaged from the former ArmPL cross-check: feed the kernel reflectors - * from a column-pivoted (rank-revealing) QR and recover X through a - * jpvt back-permutation, i.e. solve A x = b with A(:,jpvt) = Q R: + * from a column-pivoted (rank-revealing) QR, LAPACKE_?geqp3 through + * ref_geqp3 (0-based jpvt), and recover X through a jpvt back-permutation, + * i.e. solve A x = b with A(:,jpvt) = Q R: * R y = Q^T b (kernel applies Q^T), x(jpvt(j)) = y(j). * Unpivoted tests never exercise this end-to-end permuted pipeline. */ template static int run_case_pivoted(int nm, int m, int nrhs) @@ -147,7 +108,7 @@ template static int run_case_pivoted(int nm, int m, int nrhs) const auto X = mat_view(Xs.data(), m, nrhs); MatrixBatch Afac(nm, m, m), B(nm, m, nrhs), Bout(nm, m, nrhs), tau(nm, m, 1); - MatrixBatch jpvt(nm, m, 1); + MatrixBatch jpvt(nm, m, 1); for (int kk = 0; kk < nm; ++kk) { std::vector As((size_t)m * m); const auto A = mat_view(As.data(), m, m); diff --git a/tests/test_ormqr_mkl.cpp b/tests/test_ormqr_mkl.cpp index b9e2674..662cca6 100644 --- a/tests/test_ormqr_mkl.cpp +++ b/tests/test_ormqr_mkl.cpp @@ -24,7 +24,7 @@ * Assisted-by: Claude:claude-opus-4-8 Claude:claude-fable-5 */ -#include "test_mkl_util.hpp" /* compat, mkl, lapack + the MKL-free helpers */ +#include "test_mkl_util.hpp" /* compat, mkl, lapack + the shared helpers */ #include #include diff --git a/tests/test_potrf_compact.cpp b/tests/test_potrf_compact.cpp index c063a19..c6b009d 100644 --- a/tests/test_potrf_compact.cpp +++ b/tests/test_potrf_compact.cpp @@ -1,16 +1,17 @@ // test_potrf_compact.cpp // -// Self-contained validation of the templated compact Cholesky factorization -// (dpotrf_compact / spotrf_compact), its solve companion (dpotrs_compact / -// spotrs_compact) and the fused factor-and-solve (dposv_compact / -// sposv_compact), with no BLAS dependency. The reference is the unblocked -// Cholesky (LAPACK ?potf2) implemented in scalar form -- the same algorithm -// the vectorized kernel executes V lanes at a time, so a correct kernel -// matches it to working precision. +// Validation of the templated compact Cholesky factorization (dpotrf_compact +// / spotrf_compact), its solve companion (dpotrs_compact / spotrs_compact) +// and the fused factor-and-solve (dposv_compact / sposv_compact) against +// LAPACKE (test_lapack_util.hpp): the reference factor is LAPACKE_?potrf +// (recursive / blocked, so its rounding differs from the kernel's +// right-looking sweep, but the SPD factor is unique). // // Checks per (T, V, uplo, layout): -// 1. named-triangle factor == scalar potf2 factor (elementwise, ~eps*scale) -// 2. reconstruction L L^T == A (lower) / U^T U == A (upper) +// 1. named-triangle factor == LAPACKE_?potrf factor, relative to the +// reference factor's L1 norm at 20 n eps (design 7.1: the unique factor +// agrees to a few n eps between backward-stable implementations) +// 2. reconstruction L L^T == A (lower) / U^T U == A (upper), by ?trmm // 3. the strictly-opposite triangle of the compact buffer is bit-for-bit // unchanged from the input (the routine must not reference or write it) // 4. end-to-end solve: factor + ?potrs_compact recovers a known X, and @@ -27,59 +28,24 @@ #include #include -#include "test_compact_util.hpp" // compact, frand, gen_spd, MatrixBatch, pack/unpack +#include "test_lapack_util.hpp" // compact, ref_potf2, tri_apply, gen_spd, pack/unpack using namespace cbk::test; -// ----------------------- reference kernel (scalar) ------------------ -// Unblocked right-looking potf2 on a dense column-major n x n matrix, in place. -// Lower: A = L L^T, factor in the lower triangle. Upper: A = U^T U, factor in -// the upper triangle. Only the named triangle is read or written. No SPD check -// (mirrors the routine under test): a bad pivot yields NaN/Inf. - -template static void ref_potf2(char uplo, MatrixView A) -{ - assert(A.rows == A.cols); - const bool upper = (uplo == 'U' || uplo == 'u'); - const int n = A.rows; - if (!upper) { - for (int j = 0; j < n; ++j) { - T d = std::sqrt(A(j, j)); - A(j, j) = d; - T invd = T(1) / d; - for (int i = j + 1; i < n; ++i) - A(i, j) *= invd; // scale pivot column - for (int jj = j + 1; jj < n; ++jj) // rank-1 trailing update, lower - for (int i = jj; i < n; ++i) - A(i, jj) -= A(i, j) * A(jj, j); - } - } - else { - for (int j = 0; j < n; ++j) { - T d = std::sqrt(A(j, j)); - A(j, j) = d; - T invd = T(1) / d; - for (int c = j + 1; c < n; ++c) - A(j, c) *= invd; // scale pivot row - for (int c = j + 1; c < n; ++c) // rank-1 trailing update, upper - for (int r = j + 1; r <= c; ++r) - A(r, c) -= A(j, r) * A(j, c); - } - } -} - // --------------------------- one test case -------------------------- -template static int run_case(int nm, int n, char uplo, char layout) +template +static int run_case(int nm, int n, char uplo, char layout, double cond = 0.0) { const T eps = std::numeric_limits::epsilon(); const bool rowmajor = (layout == 'R' || layout == 'r'); const bool upper = (uplo == 'U' || uplo == 'u'); - // random SPD batch + scalar reference factor for this uplo + // random SPD batch (cond > 0 squeezes its spectrum) + LAPACKE's factor + // for this uplo MatrixBatch A(nm, n, n), Aref(nm, n, n); for (int idx = 0; idx < nm; ++idx) { - gen_spd(A.view(idx)); + gen_spd(A.view(idx), cond); std::copy(A[idx], A[idx] + (size_t)n * n, Aref[idx]); ref_potf2(uplo, Aref.view(idx)); } @@ -91,50 +57,52 @@ template static int run_case(int nm, int n, char uplo, char lay unpack_compact(Aout, ap.data(), n, V, rowmajor); double e_fac = 0, e_rec = 0, e_untouched = 0; - std::vector Recs((size_t)n * n); - const auto Rec = mat_view(Recs.data(), n, n); + std::vector Fs((size_t)n * n), Recs((size_t)n * n); + const auto F = mat_view(Fs.data(), n, n), Rec = mat_view(Recs.data(), n, n); for (int idx = 0; idx < nm; ++idx) { - // check 1: named-triangle factor vs scalar reference (elementwise) + const auto Fac = Aout.view(idx), Ref = Aref.view(idx), Ain = A.view(idx); + // check 1: named-triangle factor vs LAPACKE_?potrf, relative to the + // reference factor's L1 norm (the unique SPD factor: a sharp signal) // check 3: strictly-opposite triangle unchanged from the input A + double el = 0; for (int j = 0; j < n; ++j) for (int i = 0; i < n; ++i) { const bool named = upper ? (i <= j) : (i >= j); if (named) - e_fac = std::max(e_fac, - (double)std::abs(Aout(idx, i, j) - Aref(idx, i, j))); + el = std::max(el, (double)std::abs(Fac(i, j) - Ref(i, j))); else - e_untouched = std::max( - e_untouched, (double)std::abs(Aout(idx, i, j) - A(idx, i, j))); + e_untouched = + std::max(e_untouched, (double)std::abs(Fac(i, j) - Ain(i, j))); } + e_fac = std::max(e_fac, el / std::max(norm1(Ref), norm_floor)); - // check 2: reconstruction of A from the named triangle's factor - for (int i = 0; i < n; ++i) - for (int j = 0; j < n; ++j) { - double s = 0; - if (!upper) // A = L L^T : sum_l L(i,l) L(j,l), l <= min(i,j) - for (int l = 0; l <= std::min(i, j); ++l) - s += (double)Aout(idx, i, l) * (double)Aout(idx, j, l); - else // A = U^T U : sum_l U(l,i) U(l,j), l <= min(i,j) - for (int l = 0; l <= std::min(i, j); ++l) - s += (double)Aout(idx, l, i) * (double)Aout(idx, l, j); - Rec(i, j) = (T)s; - } - e_rec = std::max(e_rec, max_abs_diff(Recs.data(), A[idx], (size_t)n * n)); + // check 2: reconstruction of A from the named triangle's factor, by + // ?trmm: F is the factor with the other triangle zeroed, so + // L L^T = F L^T (side 'R') and U^T U = U^T F (side 'L') + for (int j = 0; j < n; ++j) + for (int i = 0; i < n; ++i) + F(i, j) = (upper ? (i <= j) : (i >= j)) ? Fac(i, j) : T(0); + if (!upper) + tri_apply('R', 'L', 'T', 'N', Fac, F, Rec); + else + tri_apply('L', 'U', 'T', 'N', Fac, F, Rec); + e_rec = std::max(e_rec, max_abs_diff(Recs.data(), A[idx], (size_t)n * n) / + std::max(norm1(Ain), norm_floor)); } const double scale = std::max(1, n); - const double tol_fac = 20.0 * eps * scale; // same op sequence, unique factor - const double tol_rec = 40.0 * eps * scale * scale; // O(n) accumulation in recon + const double tol_fac = 20.0 * eps * scale; // design 7.1: unique factor, n eps + const double tol_rec = 20.0 * eps * scale; // ||L L^T - A|| / ||A||: backward stable bool ok_f = e_fac <= tol_fac; bool ok_r = e_rec <= tol_rec; bool ok_u = e_untouched == 0.0; // must be bit-for-bit unchanged bool ok_i = (info == 0); - std::printf("T=%-6s V=%-2d uplo=%c lay=%c nm=%-2d n=%-3d | fac:%.1e %s rec:%.1e %s " - "untouched:%s | info=%d %s\n", - compact::name, V, uplo, layout, nm, n, e_fac, ok_f ? "OK" : "FAIL", - e_rec, ok_r ? "OK" : "FAIL", ok_u ? "OK" : "FAIL", info, - ok_i ? "OK" : "FAIL"); + std::printf("T=%-6s V=%-2d uplo=%c lay=%c nm=%-2d n=%-3d cond=%.0f | fac:%.1e %s " + "rec:%.1e %s untouched:%s | info=%d %s\n", + compact::name, V, uplo, layout, nm, n, cond, e_fac, + ok_f ? "OK" : "FAIL", e_rec, ok_r ? "OK" : "FAIL", ok_u ? "OK" : "FAIL", + info, ok_i ? "OK" : "FAIL"); return (!ok_f) + (!ok_r) + (!ok_u) + (!ok_i); } @@ -254,10 +222,11 @@ template static int run_nonspd(int n, char uplo, char layout) MatrixBatch Aout(nm, n, n); unpack_compact(Aout, ap.data(), n, V, rowmajor); - double e_spd = 0; // worst error over the SPD sibling lanes + double e_spd = 0; // worst relative error over the SPD sibling lanes bool spd_finite = true; bool bad_poisoned = false; // the non-SPD lane must carry NaN/Inf - for (int idx = 0; idx < nm; ++idx) + for (int idx = 0; idx < nm; ++idx) { + double el = 0; for (int j = 0; j < n; ++j) for (int i = 0; i < n; ++i) { const bool named = upper ? (i <= j) : (i >= j); @@ -268,9 +237,12 @@ template static int run_nonspd(int n, char uplo, char layout) } else { if (!std::isfinite((double)x)) spd_finite = false; - e_spd = std::max(e_spd, (double)std::abs(x - Aref(idx, i, j))); + el = std::max(el, (double)std::abs(x - Aref(idx, i, j))); } } + if (idx != badlane) // vs LAPACKE_?potrf, relative to its factor's norm + e_spd = std::max(e_spd, el / std::max(norm1(Aref.view(idx)), norm_floor)); + } const double tol = 20.0 * std::numeric_limits::epsilon() * std::max(1, n); bool ok_spd = spd_finite && (e_spd <= tol); // siblings uncontaminated & correct @@ -318,8 +290,11 @@ int main() fails += run_case(11, 43, u, l); // padded partial group fails += run_case(3, 3, u, l); // smallest, padded fails += run_case(16, 30, u, l); - fails += run_case(17, 24, u, l); // padded partial group + fails += run_case(17, 24, u, l); // padded partial group + fails += run_case(8, 40, u, l, 2.0); // dynamic range (cond knob) } + fails += run_case(8, 128, 'L', 'C'); // LAPACK's blocked potrf + fails += run_case(8, 128, 'U', 'R'); // End-to-end solve A X = B on SPD batches, closing the pipeline (factor + // potrs, and the fused posv), over uplo/layout/precision, padded groups, diff --git a/tests/test_potrf_mkl.cpp b/tests/test_potrf_mkl.cpp index f937873..7e6c163 100644 --- a/tests/test_potrf_mkl.cpp +++ b/tests/test_potrf_mkl.cpp @@ -1,18 +1,11 @@ /* test_potrf_mkl.cpp * - * Validation of cbk_?potrf_compact against real Intel MKL and dense LAPACK, - * through the genuine MKL Compact pipeline (mkl_dgepack_compact / - * mkl_dgeunpack_compact). This is the design document's correctness gate against - * standard dense LAPACK (cbk_dpotrf_compact_design.md section 7). - * - * Suite 1 (section 7.1) -- Factorization invariants vs dense LAPACK, for each - * (uplo, layout): a random SPD batch is factored by cbk_dpotrf_compact, - * unpacked, and per matrix checked against the LAPACK Cholesky contract: - * - reconstruction residual || L L^T - A ||_1 / ||A||_1 <= 20 n eps - * (U^T U for upper), - * - the strictly-opposite triangle is bit-for-bit unchanged from the input, - * - since the SPD factor is unique, elementwise vs LAPACKE_dpotrf, - * || L_cbk - L_lapack ||_1 / ||L_lapack||_1 <= 20 n eps. + * Validation of cbk_?potrf_compact against real Intel MKL, through the + * genuine MKL Compact pipeline (mkl_dgepack_compact / mkl_dgeunpack_compact): + * the MKL-side half of cbk_dpotrf_compact_design.md section 7. The + * dense-LAPACK invariants of section 7.1 (reconstruction, the untouched + * triangle, the factor vs LAPACKE_?potrf, over uplo / layout / cond) run in + * the portable suite, test_potrf_compact.cpp, on any LAPACKE stack. * * Suite 2 (section 7.2) -- Cross-check vs mkl_dpotrf_compact: the same packed * batch factored by both, compact buffers compared elementwise at a small @@ -37,7 +30,7 @@ * Assisted-by: Claude:claude-opus-4-8 Claude:claude-fable-5 */ -#include "test_mkl_util.hpp" /* compat, mkl, lapack + the MKL-free helpers */ +#include "test_mkl_util.hpp" /* compat, mkl + the shared helpers */ #include #include @@ -49,99 +42,6 @@ using namespace cbk::test; namespace { -/* ---------------- Suite 1: invariants vs dense LAPACK ------------------ */ - -template -int suite1(MKL_LAYOUT layout, MKL_UPLO uplo, int nm, int n, double cond) -{ - const double eps = std::numeric_limits::epsilon(); - const MKL_COMPACT_PACK fmt = mkl_get_format_compact(); - const int V = mkl::vlen(fmt); - const bool row = (layout == MKL_ROW_MAJOR); - const bool up = (uplo == MKL_UPPER); - const char ul = up ? 'U' : 'L'; - - MatrixBatch A(nm, n, n); - for (int v = 0; v < nm; ++v) - gen_spd(A.view(v), cond); - - /* pack the full symmetric A, factor with the routine under test, unpack */ - auto Ap = A.base_ptrs(); - MKL_INT sz_a = mkl::get_size(n, n, fmt, nm); - auto ap_buf = cbk::detail::mkl_alloc_bytes(sz_a); - T *ap = ap_buf.get(); - mkl::gepack(layout, n, n, Ap.data(), n, ap, n, fmt, nm); - - MKL_INT info = 99; - compat::potrf(layout, uplo, n, ap, n, &info, fmt, nm); - - MatrixBatch H(nm, n, n); - auto Hp = H.base_ptrs(); - mkl::geunpack(layout, n, n, Hp.data(), n, ap, n, fmt, nm); - - int fails = 0; - if (info != 0) { - ++fails; - std::printf(" info = %ld (expected 0)\n", (long)info); - } - - double worst_res = 0, worst_el = 0, worst_untouched = 0; - std::vector Lref(A.stride()); - for (int v = 0; v < nm; ++v) { - /* the factor is stored in `layout`; the input and the residual are the - * column-major dense side */ - const auto Hm = H.view(v, row); - const auto Am = A.view(v); - - /* reconstruction residual and untouched-triangle check */ - std::vector R(A.stride(), 0.0); - const auto Res = mat_view(R.data(), n, n); - for (int i = 0; i < n; ++i) - for (int j = 0; j < n; ++j) { - double s = 0; - int lmax = std::min(i, j); - if (!up) /* A = L L^T */ - for (int l = 0; l <= lmax; ++l) - s += Hm(i, l) * Hm(j, l); - else /* A = U^T U */ - for (int l = 0; l <= lmax; ++l) - s += Hm(l, i) * Hm(l, j); - Res(i, j) = (T)s - Am(i, j); - } - worst_res = std::max(worst_res, norm1(Res) / std::max(norm1(Am), norm_floor)); - - for (int i = 0; i < n; ++i) - for (int j = 0; j < n; ++j) { - const bool named = up ? (i <= j) : (i >= j); - if (!named) - worst_untouched = - std::max(worst_untouched, std::abs(Hm(i, j) - Am(i, j))); - } - - /* elementwise vs LAPACKE_dpotrf (unique SPD factor -> a sharp signal) */ - std::copy(A[v], A[v] + A.stride(), Lref.begin()); - lapack::potrf(LAPACK_COL_MAJOR, ul, n, Lref.data(), n); - const auto Lr = mat_view(Lref.data(), n, n); - double el = 0, lref_norm = 0; - for (int i = 0; i < n; ++i) - for (int j = 0; j < n; ++j) { - const bool named = up ? (i <= j) : (i >= j); - if (named) el = std::max(el, std::abs(Hm(i, j) - Lr(i, j))); - } - lref_norm = norm1(Lr); /* triangular factor L1 norm */ - worst_el = std::max(worst_el, el / std::max(lref_norm, norm_floor)); - } - - const double rtol = 20.0 * n * eps; - bool ok = (worst_res <= rtol) && (worst_el <= rtol) && (worst_untouched == 0.0); - fails += !ok; - std::printf(" [suite1] %s uplo=%c V=%-2d nm=%-2d n=%-3d cond=%.0f | res %.2e " - "el %.2e (%.1e) untouched %.0e %s\n", - row ? "row" : "col", ul, V, nm, n, cond, worst_res, worst_el, rtol, - worst_untouched, ok ? "OK" : "FAIL"); - return fails; -} - /* ---------------- Suite 2: cross-check vs mkl_dpotrf_compact ----------- */ template int suite2(MKL_LAYOUT layout, MKL_UPLO uplo, int nm, int n) @@ -358,18 +258,8 @@ template int run_suites() int fails = 0; - /* Suite 1: invariants vs dense LAPACK, over (uplo, layout), sizes and cond */ const MKL_LAYOUT lays[] = {MKL_COL_MAJOR, MKL_ROW_MAJOR}; const MKL_UPLO ups[] = {MKL_LOWER, MKL_UPPER}; - for (MKL_LAYOUT L : lays) - for (MKL_UPLO U : ups) { - fails += suite1(L, U, 8, 30, 0.0); - fails += suite1(L, U, 16, 60, 0.0); - fails += suite1(L, U, 11, 43, 0.0); /* padded partial group */ - fails += suite1(L, U, 8, 40, 2.0); /* dynamic range (cond knob) */ - } - fails += suite1(MKL_COL_MAJOR, MKL_LOWER, 8, 128, 0.0); - fails += suite1(MKL_COL_MAJOR, MKL_LOWER, 4, 3, 0.0); /* smallest, padded */ /* Suite 2: cross-check vs mkl_dpotrf_compact, both layouts and uplo */ for (MKL_LAYOUT L : lays) diff --git a/tests/test_sytrfnp_compact.cpp b/tests/test_sytrfnp_compact.cpp index 237b673..538a586 100644 --- a/tests/test_sytrfnp_compact.cpp +++ b/tests/test_sytrfnp_compact.cpp @@ -1,12 +1,15 @@ // test_sytrfnp_compact.cpp // -// Self-contained validation of the templated compact unpivoted LDL^T -// factorization (dsytrfnp_compact / ssytrfnp_compact), its solve companion -// (dsytrsnp_compact / ssytrsnp_compact) and the fused factor-and-solve -// (dsysvnp_compact / ssysvnp_compact), with no BLAS dependency. The reference is -// the unblocked unpivoted factorization implemented in scalar form -- the same -// algorithm the vectorized kernel executes V lanes at a time, so a correct -// kernel matches it to working precision. +// Validation of the templated compact unpivoted LDL^T factorization +// (dsytrfnp_compact / ssytrfnp_compact), its solve companion (dsytrsnp_compact +// / ssytrsnp_compact) and the fused factor-and-solve (dsysvnp_compact / +// ssysvnp_compact). LAPACK has no unpivoted LDL^T (?sytrf pivots, so its +// factors differ elementwise), so the reference is the one scalar routine the +// suites still hand-roll, ref_sytf2np below: the unblocked unpivoted +// factorization -- the same algorithm the vectorized kernel executes V lanes +// at a time, so a correct kernel matches it to working precision. The +// reference itself is validated first, against the library: its (L, D) must +// reconstruct A through cblas_?trmm / ?gemm (test_reference). // // Checks per (T, V, uplo, layout), on symmetric *indefinite* batches: // 1. named-triangle factor == scalar reference factor (elementwise) @@ -32,7 +35,7 @@ #include #include -#include "test_compact_util.hpp" // compact, gen_sym_ldlt, ldlt_reconstruct, pack/unpack +#include "test_lapack_util.hpp" // compact, matmul, tri_apply, gen_sym_ldlt, pack/unpack using namespace cbk::test; @@ -66,6 +69,54 @@ template static void ref_sytf2np(char uplo, MatrixView As) } } +// ------------- the reference itself, vs the library ------------------ +// ref_sytf2np is the only hand-rolled reference left in the suites, so it is +// checked against LAPACKE/CBLAS before anything is checked against it: its +// (L, D) must reconstruct A. Lower: W = L D (columns of the unit-lower L +// scaled by D), then W := W L^T by ?trmm (side 'R', unit diagonal). Upper: +// W = D U (rows scaled), then W := U^T W (side 'L'). Gated at +// ||L D L^T - A||_1 / ||A||_1 <= 20 n eps, the design's 7.1 residual. + +template static int test_reference(int n, char uplo) +{ + const T eps = std::numeric_limits::epsilon(); + const bool upper = (uplo == 'U' || uplo == 'u'); + const int trials = 4; + double worst = 0; + std::vector As((size_t)n * n), Fs((size_t)n * n), Ws((size_t)n * n), + Rs((size_t)n * n); + const auto A = mat_view(As.data(), n, n), F = mat_view(Fs.data(), n, n), + W = mat_view(Ws.data(), n, n), R = mat_view(Rs.data(), n, n); + for (int t = 0; t < trials; ++t) { + gen_sym_ldlt(A); + copy_matrix(A, F); + ref_sytf2np(uplo, F); + // W: the unit factor with D folded in, the other triangle zeroed + for (int j = 0; j < n; ++j) + for (int i = 0; i < n; ++i) { + if (i == j) + W(i, j) = F(i, i); + else if (upper) + W(i, j) = (i < j) ? F(i, i) * F(i, j) : T(0); // D U + else + W(i, j) = (i > j) ? F(i, j) * F(j, j) : T(0); // L D + } + if (upper) + tri_apply('L', 'U', 'T', 'U', F, W, R); // R = U^T (D U) + else + tri_apply('R', 'L', 'T', 'U', F, W, R); // R = (L D) L^T + worst = std::max(worst, max_abs_diff(Rs.data(), As.data(), As.size()) / + std::max(norm1(A), norm_floor)); + } + const double tol = 20.0 * eps * std::max(1, n); + const bool ok = worst <= tol; + std::printf( + "T=%-6s reference ref_sytf2np uplo=%c n=%-3d | ||LDL^T-A||/||A|| %.1e (%.1e) " + "%s\n", + compact::name, uplo, n, worst, tol, ok ? "OK" : "FAIL"); + return !ok; +} + // --------------------------- one test case -------------------------- template static int run_case(int nm, int n, char uplo, char layout) @@ -364,10 +415,17 @@ int main() int fails = 0; fails += test_validation(); - // Full feature matrix: both uplo x both layouts, several (V, n, nm), - // including padded partial final groups (nm not a multiple of V). + // The hand-rolled reference against the library first. const char uplos[] = {'L', 'U'}; const char lays[] = {'C', 'R'}; + for (char u : uplos) { + fails += test_reference(30, u); + fails += test_reference(43, u); + fails += test_reference(24, u); + } + + // Full feature matrix: both uplo x both layouts, several (V, n, nm), + // including padded partial final groups (nm not a multiple of V). for (char u : uplos) for (char l : lays) { fails += run_case(4, 16, u, l); diff --git a/tests/test_sytrfnp_mkl.cpp b/tests/test_sytrfnp_mkl.cpp index aba819c..5d42cd0 100644 --- a/tests/test_sytrfnp_mkl.cpp +++ b/tests/test_sytrfnp_mkl.cpp @@ -39,7 +39,7 @@ * Assisted-by: Claude:claude-fable-5 */ -#include "test_mkl_util.hpp" /* compat, mkl + the MKL-free helpers */ +#include "test_mkl_util.hpp" /* compat, mkl + the shared helpers */ #include #include diff --git a/tests/test_trsm_compact.cpp b/tests/test_trsm_compact.cpp index 4bdcbc3..7057cde 100644 --- a/tests/test_trsm_compact.cpp +++ b/tests/test_trsm_compact.cpp @@ -1,18 +1,17 @@ // test_trsm_compact.cpp // -// Self-contained validation of the templated compact triangular solve -// (dtrsm_compact / strsm_compact), with no BLAS dependency. The reference is a -// scalar BLAS-?trsm implemented directly -- the same math the vectorized kernel -// executes V lanes at a time, so a correct kernel matches it to working -// precision. +// Validation of the templated compact triangular solve (dtrsm_compact / +// strsm_compact) against CBLAS (test_lapack_util.hpp): the reference is +// cblas_?trsm -- the same math the vectorized kernel executes V lanes at a +// time -- and the residual is formed by cblas_?trmm. // // Two parts: // 1. C API argument validation -- exercises the LAPACK/BLAS-style info = -j // contract of the portable entry points. // 2. Numerical correctness over side / uplo / transa / diag, precisions and -// interleave widths (column-major, the tuned path): forward error vs a -// scalar reference solve, plus the solve's own residual ||op(A) X - alpha B|| -// formed with an independent triangular multiply (so a bug shared by the +// interleave widths (column-major, the tuned path): forward error vs the +// library's solve, plus the solve's own residual ||op(A) X - alpha B|| +// formed with the library's triangular multiply (so a bug shared by the // reference and the kernel cannot pass unseen). // // Assisted-by: Claude:claude-opus-4-8 Claude:claude-fable-5 @@ -25,12 +24,12 @@ #include #include -#include "test_compact_util.hpp" // compact, frand, gen_tri, tri_apply, MatrixBatch, pack/unpack +#include "test_lapack_util.hpp" // compact, ref_trsm, tri_apply, gen_tri, pack/unpack using namespace cbk::test; -// The scalar reference is ref_trsm of test_compact_util.hpp (dense BLAS ?trsm -// over views, in place); tri_apply forms the residual independently. +// The reference is ref_trsm of test_lapack_util.hpp (cblas_?trsm over views, +// in place); tri_apply (cblas_?trmm) forms the residual independently. // Padded pack slots carry the identity -- for a triangular A that is a unit // diagonal, so the kernel's divisions never hit a zero pivot in the padding. @@ -71,9 +70,9 @@ static int run_case(char side, char uplo, char transa, char diag, int nm, int m, MatrixBatch Bout(nm, m, n); unpack_compact(Bout, bp.data(), m, V); - // Two gates: (1) forward error vs the scalar reference solve, and (2) the - // solve's own defining residual ||op(A) X - alpha B||, formed with an - // independent triangular multiply -- so a bug shared by ref_trsm and the + // Two gates: (1) forward error vs the library's solve, and (2) the + // solve's own defining residual ||op(A) X - alpha B||, formed with the + // library's triangular multiply -- so a bug shared by ref_trsm and the // kernel cannot slip through (the ?trsm analogue of the reconstruction / // round-trip identities the geqrf/potrf/ormqr self-tests check). double worst_fwd = 0, worst_res = 0; @@ -152,7 +151,7 @@ int main() int fails = 0; fails += test_validation(); - std::printf("numerical cases (column-major, vs scalar reference):\n"); + std::printf("numerical cases (column-major, vs cblas_?trsm):\n"); // full side x uplo x transa x diag matrix at a representative shape/width for (char side : {'L', 'R'}) for (char uplo : {'U', 'L'}) diff --git a/tests/test_trsm_mkl.cpp b/tests/test_trsm_mkl.cpp index b6355c9..f8618fe 100644 --- a/tests/test_trsm_mkl.cpp +++ b/tests/test_trsm_mkl.cpp @@ -17,7 +17,7 @@ * Assisted-by: Claude:claude-opus-4-8 Claude:claude-fable-5 */ -#include "test_mkl_util.hpp" /* compat, mkl, lapack + the MKL-free helpers */ +#include "test_mkl_util.hpp" /* compat, mkl, lapack + the shared helpers */ #include #include From b9d2ac3068e8d7d2d67c8d99da955e14738ec713 Mon Sep 17 00:00:00 2001 From: Claude Date: Mon, 14 Sep 2026 22:17:43 +0000 Subject: [PATCH 2/8] Test against Accelerate on the macOS runners .github/workflows/macos.yml runs the suites on macos-latest on the same triggers as the Linux CI: it builds accelerate-lapacke (Netlib's LAPACKE over Accelerate's LAPACK; a Fortran compiler from Homebrew's gcc for its configure) into ~/.local and configures cbk with -DCBK_TEST_LAPACK=accelerate. FindLAPACKE.cmake's Accelerate branch looks up the installed `lapacke` package under its own lower-case config name (so the module is not re-entered), takes CBLAS from Accelerate's vecLib sub-framework through the SDK's framework directory (xcrun --show-sdk-path), and never includes the umbrella, whose LAPACK prototypes are the same C symbols lapacke.h declares with different parameter types. The test util spells the CBLAS layout enum CBLAS_ORDER, the one name every CBLAS header defines. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_01GdNaquWMN3b39gvTsXZyjn --- .claude/CLAUDE.md | 2 +- .github/workflows/macos.yml | 33 +++++++++++++++++++++++++-------- CHANGELOG.md | 3 ++- cmake/FindLAPACKE.cmake | 36 ++++++++++++++++++++++++++---------- docs/PLANS.md | 7 ++++--- docs/building.md | 8 +++++--- tests/test_lapack_util.hpp | 28 +++++++++++++++++++--------- 7 files changed, 82 insertions(+), 35 deletions(-) diff --git a/.claude/CLAUDE.md b/.claude/CLAUDE.md index e5d0a4f..2aa40c6 100644 --- a/.claude/CLAUDE.md +++ b/.claude/CLAUDE.md @@ -185,7 +185,7 @@ workspace contract, or the benchmarks' threading. enforces: both alignment faults so far (issue #34 and the one above) were invisible to gcc. The default (portable) tree needs only the compiler and a LAPACKE stack for the tests (CI runs OpenBLAS and Netlib under both - compilers): + compilers, and Accelerate on macOS through `.github/workflows/macos.yml`): ```sh CXX=clang++ cmake -S . -B build-clang -DCBK_TEST_LAPACK=openblas && cmake --build build-clang && ctest --test-dir build-clang diff --git a/.github/workflows/macos.yml b/.github/workflows/macos.yml index a42c1d5..059220e 100644 --- a/.github/workflows/macos.yml +++ b/.github/workflows/macos.yml @@ -1,18 +1,24 @@ name: macOS -# The test suites on macOS, against Apple's Accelerate through -# accelerate-lapacke (https://github.com/lepus2589/accelerate-lapacke), which -# builds the LAPACKE interface Accelerate lacks. macOS runners are expensive, -# so this runs only for release tags and on manual dispatch, never per push -# or pull request. Unverified as of its introduction (no macOS machine was at -# hand): expect to adjust the accelerate-lapacke preset or the CBLAS header -# mapping (CBK_LAPACKE_ACCELERATE in tests/test_lapack_util.hpp) on the first -# run. +# The test suites on Apple silicon against Accelerate: Apple's framework +# provides LAPACK and CBLAS but no LAPACKE, so the job first builds +# accelerate-lapacke (https://github.com/lepus2589/accelerate-lapacke, Netlib's +# LAPACKE over Accelerate's LAPACK 3.9.1 interface; its configure needs a +# Fortran compiler, which Homebrew's gcc supplies) into ~/.local, then +# configures cbk with -DCBK_TEST_LAPACK=accelerate (cmake/FindLAPACKE.cmake: +# the `lapacke` package plus from the SDK). macOS runners +# cost several times a Linux minute, so this is one job, not a matrix. on: push: + branches: [main] tags: ['v*'] + pull_request: workflow_dispatch: +concurrency: + group: ${{ github.workflow }}-${{ github.ref }} + cancel-in-progress: true + jobs: accelerate: name: AppleClang / Accelerate @@ -20,11 +26,22 @@ jobs: steps: - uses: actions/checkout@v4 + # accelerate-lapacke's project enables Fortran (Netlib's LAPACK build + # system does, even though only its C LAPACKE sources are compiled and + # the result links nothing but Accelerate). + - name: Install a Fortran compiler + run: | + brew install gcc + fc=$(ls "$(brew --prefix)"/bin/gfortran* | head -n 1) + echo "FC=$fc" >> "$GITHUB_ENV" + "$fc" --version | head -n 1 + - name: Build accelerate-lapacke run: | git clone --depth 1 https://github.com/lepus2589/accelerate-lapacke.git "$RUNNER_TEMP/accelerate-lapacke" cd "$RUNNER_TEMP/accelerate-lapacke" cmake --workflow --preset accelerate-lapacke32-install + ls "$HOME/.local/include/lapacke.h" "$HOME/.local/lib/" - name: Configure run: > diff --git a/CHANGELOG.md b/CHANGELOG.md index bde4ee3..4b3b1f7 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -15,7 +15,8 @@ release may change the API. scalar references (issue #27). The dense-LAPACK cross-checks of the MKL suites run in the portable suites, on any stack; the unpivoted LDL^T keeps its scalar reference, itself validated against CBLAS. Building the tests - now needs such a stack; the library does not. + now needs such a stack; the library does not. A macOS workflow tests + against Accelerate through accelerate-lapacke. ## [0.1.0] - 2026-09-13 diff --git a/cmake/FindLAPACKE.cmake b/cmake/FindLAPACKE.cmake index 3befa75..4f27e74 100644 --- a/cmake/FindLAPACKE.cmake +++ b/cmake/FindLAPACKE.cmake @@ -5,7 +5,7 @@ # LAPACKE::LAPACKE -- the include dirs, the link line, and one compile # definition naming the stack's headers: # CBK_LAPACKE_MKL , -# CBK_LAPACKE_ACCELERATE , +# CBK_LAPACKE_ACCELERATE , # (none) , # # LAPACKE_VENDOR Any (default) | MKL | OpenBLAS | Netlib | Accelerate @@ -35,9 +35,10 @@ # first (PATH_SUFFIXES), which is where the distro keeps them. # Accelerate macOS: Apple's Accelerate framework provides LAPACK and # CBLAS but no LAPACKE; https://github.com/lepus2589/accelerate-lapacke -# builds one on top of it and installs a `lapacke` CMake -# package (find_package(LAPACKE CONFIG)). Unverified in CI: -# the macOS workflow runs on tags and by hand only. +# builds one on top of it (Netlib's LAPACKE, so it needs a +# Fortran compiler to configure) and installs a `lapacke` +# CMake package. .github/workflows/macos.yml builds it on the +# macOS runners. # # ILP64 is supported through MKL only (MKL_ILP64 comes with MKL::Compact); # the other stacks are used at their default 32-bit integer interface. @@ -67,7 +68,7 @@ function(_lapacke_link_test vendor includes libs defs outvar) #include #elif defined(CBK_LAPACKE_ACCELERATE) #include - #include + #include #else #include #include @@ -89,6 +90,7 @@ set(_lapacke_found_vendor "") set(_lapacke_includes "") set(_lapacke_libs "") set(_lapacke_defs "") +set(_lapacke_options "") # ---- MKL ------------------------------------------------------------------- if(NOT _lapacke_found_vendor AND LAPACKE_VENDOR STREQUAL "MKL") @@ -195,17 +197,30 @@ if(NOT _lapacke_found_vendor AND LAPACKE_VENDOR MATCHES "^(Any|Netlib)$") endif() # ---- Accelerate (macOS) ----------------------------------------------------- +# accelerate-lapacke installs Netlib's lapacke-config.cmake (looked up under +# its own lower-case file name, so this module is not re-entered) with the +# `lapacke` target; CBLAS is Accelerate's, reached as through +# the framework directory of the SDK's Accelerate.framework, which +# `xcrun --show-sdk-path` locates. Only the vecLib sub-framework's header is +# included, never the umbrella: its LAPACK +# prototypes and the ones lapacke.h declares are the same C symbols with +# different parameter types. if(NOT _lapacke_found_vendor AND LAPACKE_VENDOR STREQUAL "Accelerate") - find_package(LAPACKE CONFIG QUIET) # accelerate-lapacke's package: target `lapacke` + find_package(lapacke CONFIG QUIET) find_library(LAPACKE_Accelerate_FRAMEWORK Accelerate) mark_as_advanced(LAPACKE_Accelerate_FRAMEWORK) - if(TARGET lapacke AND LAPACKE_Accelerate_FRAMEWORK) + execute_process(COMMAND xcrun --show-sdk-path + OUTPUT_VARIABLE _ac_sdk OUTPUT_STRIP_TRAILING_WHITESPACE ERROR_QUIET) + set(_ac_frameworks "${_ac_sdk}/System/Library/Frameworks/Accelerate.framework/Frameworks") + if(TARGET lapacke AND LAPACKE_Accelerate_FRAMEWORK AND EXISTS "${_ac_frameworks}/vecLib.framework") set(_ac_libs lapacke "${LAPACKE_Accelerate_FRAMEWORK}") - _lapacke_link_test(Accelerate "" "${_ac_libs}" "-DCBK_LAPACKE_ACCELERATE" _ok) + set(_ac_defs -DCBK_LAPACKE_ACCELERATE -DACCELERATE_NEW_LAPACK "-F${_ac_frameworks}") + _lapacke_link_test(Accelerate "" "${_ac_libs}" "${_ac_defs}" _ok) if(_ok) set(_lapacke_found_vendor Accelerate) set(_lapacke_libs ${_ac_libs}) - set(_lapacke_defs CBK_LAPACKE_ACCELERATE) + set(_lapacke_defs CBK_LAPACKE_ACCELERATE ACCELERATE_NEW_LAPACK) + set(_lapacke_options "-F${_ac_frameworks}") endif() endif() endif() @@ -247,6 +262,7 @@ if(LAPACKE_FOUND AND NOT TARGET LAPACKE::LAPACKE) set_target_properties(LAPACKE::LAPACKE PROPERTIES INTERFACE_INCLUDE_DIRECTORIES "${LAPACKE_INCLUDE_DIRS}" INTERFACE_LINK_LIBRARIES "${LAPACKE_LIBRARIES}" - INTERFACE_COMPILE_DEFINITIONS "${_lapacke_defs}") + INTERFACE_COMPILE_DEFINITIONS "${_lapacke_defs}" + INTERFACE_COMPILE_OPTIONS "${_lapacke_options}") message(STATUS "LAPACKE for the tests: ${LAPACKE_VENDOR_FOUND} (${LAPACKE_LIBRARIES})") endif() diff --git a/docs/PLANS.md b/docs/PLANS.md index c4b2eee..6382afc 100644 --- a/docs/PLANS.md +++ b/docs/PLANS.md @@ -222,9 +222,10 @@ driver (`docs/cbk_dsytrfnp_compact_design.md`); MKL has no compact moved into the portable suites, on every stack. Where a reference is LAPACK's blocked driver, the gate is relative to the operand norms at a multiple of `n eps` (design 7.1), not elementwise at `~eps`; only `?geqr2` - stays an elementwise `~eps` comparison. macOS through Accelerate + - accelerate-lapacke has a tag-only workflow (`.github/workflows/macos.yml`), - not yet exercised. + stays an elementwise `~eps` comparison. macOS runs against Accelerate + + accelerate-lapacke on the macOS runners (`.github/workflows/macos.yml`), on + the same triggers as the Linux CI; written without a macOS machine at hand, + so its first run is the check. - **Planned: LAPACK-style test coverage.** Adopt the testing approaches of the reference LAPACK repository (its `TESTING/LIN` drivers): `?latms`-style generators with prescribed condition number and spectral distribution, and diff --git a/docs/building.md b/docs/building.md index a68f854..019dde3 100644 --- a/docs/building.md +++ b/docs/building.md @@ -29,9 +29,11 @@ ctest --test-dir build --output-on-failure On macOS, Accelerate provides LAPACK and CBLAS but no LAPACKE; [accelerate-lapacke](https://github.com/lepus2589/accelerate-lapacke) builds -one on top of it, and `-DCBK_TEST_LAPACK=accelerate` with -`-DCMAKE_PREFIX_PATH` at its install looks for that (the tag-only -`.github/workflows/macos.yml`; not yet exercised on a macOS machine). +one on top of it (its configure needs a Fortran compiler, `brew install gcc`), +and `-DCBK_TEST_LAPACK=accelerate` with `-DCMAKE_PREFIX_PATH` at its install +(`~/.local` by default) looks for that: the `lapacke` package plus +Accelerate's own CBLAS through ``. `.github/workflows/macos.yml` +runs exactly this on the macOS runners, one job next to the Linux matrix. ## The Intel MKL extension (`-DCBK_WITH_MKL=ON`) diff --git a/tests/test_lapack_util.hpp b/tests/test_lapack_util.hpp index 24f619f..5f9fa2c 100644 --- a/tests/test_lapack_util.hpp +++ b/tests/test_lapack_util.hpp @@ -14,7 +14,8 @@ // cmake/FindLAPACKE.cmake through one definition on the LAPACKE::LAPACKE // target: MKL's own / when the tests are built // with the MKL extension (a test binary must never link two LAPACK stacks), -// the plain / of OpenBLAS or Netlib otherwise. +// the plain / of OpenBLAS or Netlib otherwise, and on +// macOS accelerate-lapacke's with Accelerate's CBLAS. // Integer arguments are lapack_int (MKL_INT under MKL, so an ilp64 MKL build // gets its 64-bit integers); the CBLAS calls take the same lapack_int and // convert to the stack's own BLAS integer, which matches on every stack the @@ -31,8 +32,16 @@ #include #include #elif defined(CBK_LAPACKE_ACCELERATE) -#include +// macOS: LAPACKE from accelerate-lapacke (Netlib's LAPACKE over Accelerate's +// LAPACK 3.9.1 interface), CBLAS from Accelerate's vecLib sub-framework +// directly -- not through the umbrella, whose LAPACK +// prototypes would collide with the ones lapacke.h pulls in. The framework +// search path for is on the LAPACKE::LAPACKE target. +#ifndef ACCELERATE_NEW_LAPACK +#define ACCELERATE_NEW_LAPACK +#endif #include +#include #else #include #include @@ -83,16 +92,16 @@ template <> struct lapack { static lapack_int gels(int layout, char trans, lapack_int m, lapack_int n, \ lapack_int nrhs, T *a, lapack_int lda, T *b, lapack_int ldb) \ { return LAPACKE_##p##gels(layout, trans, m, n, nrhs, a, lda, b, ldb); } \ - static void gemm(CBLAS_LAYOUT layout, CBLAS_TRANSPOSE ta, CBLAS_TRANSPOSE tb, \ + static void gemm(CBLAS_ORDER layout, CBLAS_TRANSPOSE ta, CBLAS_TRANSPOSE tb, \ lapack_int m, lapack_int n, lapack_int k, T alpha, const T *a, \ lapack_int lda, const T *b, lapack_int ldb, T beta, T *c, \ lapack_int ldc) \ { cblas_##p##gemm(layout, ta, tb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc); } \ - static void trsm(CBLAS_LAYOUT layout, CBLAS_SIDE side, CBLAS_UPLO uplo, \ + static void trsm(CBLAS_ORDER layout, CBLAS_SIDE side, CBLAS_UPLO uplo, \ CBLAS_TRANSPOSE ta, CBLAS_DIAG diag, lapack_int m, lapack_int n, \ T alpha, const T *a, lapack_int lda, T *b, lapack_int ldb) \ { cblas_##p##trsm(layout, side, uplo, ta, diag, m, n, alpha, a, lda, b, ldb); } \ - static void trmm(CBLAS_LAYOUT layout, CBLAS_SIDE side, CBLAS_UPLO uplo, \ + static void trmm(CBLAS_ORDER layout, CBLAS_SIDE side, CBLAS_UPLO uplo, \ CBLAS_TRANSPOSE ta, CBLAS_DIAG diag, lapack_int m, lapack_int n, \ T alpha, const T *a, lapack_int lda, T *b, lapack_int ldb) \ { cblas_##p##trmm(layout, side, uplo, ta, diag, m, n, alpha, a, lda, b, ldb); } \ @@ -106,7 +115,8 @@ CBK_TEST_LAPACK_DISPATCH(float, s) // ----------------------- views as LAPACKE / CBLAS arguments ---------- // A MatrixView carries its layout as strides (si, sj); LAPACKE and CBLAS take -// one layout flag per call plus a leading dimension per operand. A view maps +// one layout flag per call plus a leading dimension per operand (the CBLAS +// one is spelled CBLAS_ORDER, the name every CBLAS header agrees on). A view maps // without a copy when one stride is 1 and the other is at least the extent // along it -- column-major (si == 1, ld = sj) or row-major (sj == 1, ld = si); // the ambiguous degenerate shapes (a single row or column) are read as @@ -129,7 +139,7 @@ template int lapack_layout(const Mv &M) return is_rowmajor(M) ? LAPACK_ROW_MAJOR : LAPACK_COL_MAJOR; } -template CBLAS_LAYOUT cblas_layout(const Mv &M) +template CBLAS_ORDER cblas_layout(const Mv &M) { return is_rowmajor(M) ? CblasRowMajor : CblasColMajor; } @@ -175,8 +185,8 @@ template class Staged { return; } buf_.resize((std::size_t)M.rows * M.cols); - const auto w = mat_view(buf_.data(), M.rows, M.cols, rowmajor ? M.cols : M.rows, - rowmajor); + const auto w = + mat_view(buf_.data(), M.rows, M.cols, rowmajor ? M.cols : M.rows, rowmajor); copy_matrix(M, w); use_ = MatrixView{w.data, w.si, w.sj, w.rows, w.cols}; } From 050dad73f14521e2af8b43c9e976f2fe19418d58 Mon Sep 17 00:00:00 2001 From: Claude Date: Mon, 14 Sep 2026 22:21:09 +0000 Subject: [PATCH 3/8] macOS workflow: hand accelerate-lapacke the SDK path accelerate-lapack refuses to configure without CMAKE_OSX_SYSROOT (it matches the SDK's Accelerate stubs to the system), and a workflow preset takes no cache variables, so the build is a configure and a build call. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_01GdNaquWMN3b39gvTsXZyjn --- .github/workflows/macos.yml | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/.github/workflows/macos.yml b/.github/workflows/macos.yml index 059220e..51c2767 100644 --- a/.github/workflows/macos.yml +++ b/.github/workflows/macos.yml @@ -36,11 +36,16 @@ jobs: echo "FC=$fc" >> "$GITHUB_ENV" "$fc" --version | head -n 1 + # accelerate-lapack (its dependency) insists on an explicit SDK + # (CMAKE_OSX_SYSROOT) to match Accelerate's text-based stubs to the + # system; the workflow preset takes no cache variables, so configure + # and build are two calls. - name: Build accelerate-lapacke run: | git clone --depth 1 https://github.com/lepus2589/accelerate-lapacke.git "$RUNNER_TEMP/accelerate-lapacke" cd "$RUNNER_TEMP/accelerate-lapacke" - cmake --workflow --preset accelerate-lapacke32-install + cmake --preset accelerate-lapacke32 -DCMAKE_OSX_SYSROOT="$(xcrun --show-sdk-path)" + cmake --build --preset accelerate-lapacke32-install ls "$HOME/.local/include/lapacke.h" "$HOME/.local/lib/" - name: Configure From de20b9fc241930f76b2e78084738e1184c38d914 Mon Sep 17 00:00:00 2001 From: Claude Date: Mon, 14 Sep 2026 22:30:20 +0000 Subject: [PATCH 4/8] Tests ask for C++17 themselves; cache accelerate-lapacke on macOS AppleClang's default standard is older than gcc's and clang's on Linux, and the library's cxx_std_17 is a PRIVATE build requirement of the kernels that nothing propagates, so every test, example and benchmark target now requests the feature itself (reproduced on Linux with -DCMAKE_CXX_STANDARD=14; green with the change). The macOS workflow installs accelerate-lapacke, pinned to v2.0.0, into its own prefix and caches that under a key naming the release and the macOS SDK it was linked against; the Fortran compiler is installed only on a cache miss. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_01GdNaquWMN3b39gvTsXZyjn --- .github/workflows/macos.yml | 60 +++++++++++++++++++++++++------------ CMakeLists.txt | 8 ++++- docs/building.md | 2 +- 3 files changed, 49 insertions(+), 21 deletions(-) diff --git a/.github/workflows/macos.yml b/.github/workflows/macos.yml index 51c2767..bd05ddd 100644 --- a/.github/workflows/macos.yml +++ b/.github/workflows/macos.yml @@ -3,11 +3,17 @@ name: macOS # The test suites on Apple silicon against Accelerate: Apple's framework # provides LAPACK and CBLAS but no LAPACKE, so the job first builds # accelerate-lapacke (https://github.com/lepus2589/accelerate-lapacke, Netlib's -# LAPACKE over Accelerate's LAPACK 3.9.1 interface; its configure needs a -# Fortran compiler, which Homebrew's gcc supplies) into ~/.local, then +# LAPACKE over Accelerate's LAPACK interface; its configure needs a Fortran +# compiler, which Homebrew's gcc supplies) into its own prefix, then # configures cbk with -DCBK_TEST_LAPACK=accelerate (cmake/FindLAPACKE.cmake: # the `lapacke` package plus from the SDK). macOS runners -# cost several times a Linux minute, so this is one job, not a matrix. +# cost several times a Linux minute, so this is one job, not a matrix, and +# the LAPACKE build (several minutes, most of it configuring Netlib LAPACK's +# Fortran tree) is cached: the install prefix is restored by actions/cache +# under a key naming the pinned accelerate-lapacke release and the macOS +# SDK the dylib was linked against (its Accelerate stubs must match the +# system), so a runner image upgrade or a version bump rebuilds it and +# nothing else does. on: push: branches: [main] @@ -19,6 +25,10 @@ concurrency: group: ${{ github.workflow }}-${{ github.ref }} cancel-in-progress: true +env: + ACCELERATE_LAPACKE_REF: v2.0.0 + ACCELERATE_LAPACKE_PREFIX: ${{ github.workspace }}/../accelerate-lapacke + jobs: accelerate: name: AppleClang / Accelerate @@ -26,34 +36,46 @@ jobs: steps: - uses: actions/checkout@v4 + - name: Identify the SDK + id: sdk + run: | + echo "path=$(xcrun --show-sdk-path)" >> "$GITHUB_OUTPUT" + echo "version=$(xcrun --show-sdk-version)" >> "$GITHUB_OUTPUT" + echo "macOS $(sw_vers -productVersion), SDK $(xcrun --show-sdk-version)" + + - name: Cache accelerate-lapacke + id: lapacke-cache + uses: actions/cache@v4 + with: + path: ${{ env.ACCELERATE_LAPACKE_PREFIX }} + key: accelerate-lapacke-${{ env.ACCELERATE_LAPACKE_REF }}-${{ runner.os }}-${{ runner.arch }}-sdk${{ steps.sdk.outputs.version }} + # accelerate-lapacke's project enables Fortran (Netlib's LAPACK build # system does, even though only its C LAPACKE sources are compiled and - # the result links nothing but Accelerate). - - name: Install a Fortran compiler - run: | - brew install gcc - fc=$(ls "$(brew --prefix)"/bin/gfortran* | head -n 1) - echo "FC=$fc" >> "$GITHUB_ENV" - "$fc" --version | head -n 1 - - # accelerate-lapack (its dependency) insists on an explicit SDK - # (CMAKE_OSX_SYSROOT) to match Accelerate's text-based stubs to the - # system; the workflow preset takes no cache variables, so configure - # and build are two calls. + # the result links nothing but Accelerate). Its accelerate-lapack + # dependency insists on an explicit SDK (CMAKE_OSX_SYSROOT); the + # workflow preset takes no cache variables, so configure and build are + # two calls, and the install prefix overrides the preset's ~/.local. - name: Build accelerate-lapacke + if: steps.lapacke-cache.outputs.cache-hit != 'true' run: | - git clone --depth 1 https://github.com/lepus2589/accelerate-lapacke.git "$RUNNER_TEMP/accelerate-lapacke" + brew install gcc + export FC=$(ls "$(brew --prefix)"/bin/gfortran* | head -n 1) + "$FC" --version | head -n 1 + git clone --depth 1 --branch "$ACCELERATE_LAPACKE_REF" https://github.com/lepus2589/accelerate-lapacke.git "$RUNNER_TEMP/accelerate-lapacke" cd "$RUNNER_TEMP/accelerate-lapacke" - cmake --preset accelerate-lapacke32 -DCMAKE_OSX_SYSROOT="$(xcrun --show-sdk-path)" + cmake --preset accelerate-lapacke32 \ + -DCMAKE_OSX_SYSROOT="${{ steps.sdk.outputs.path }}" \ + -DCMAKE_INSTALL_PREFIX="$ACCELERATE_LAPACKE_PREFIX" cmake --build --preset accelerate-lapacke32-install - ls "$HOME/.local/include/lapacke.h" "$HOME/.local/lib/" + ls "$ACCELERATE_LAPACKE_PREFIX/include/lapacke.h" "$ACCELERATE_LAPACKE_PREFIX/lib/" - name: Configure run: > cmake -S . -B build -DCMAKE_BUILD_TYPE=Release -DCBK_TEST_LAPACK=accelerate - -DCMAKE_PREFIX_PATH="$HOME/.local" + -DCMAKE_PREFIX_PATH="$ACCELERATE_LAPACKE_PREFIX" - name: Build run: cmake --build build -j diff --git a/CMakeLists.txt b/CMakeLists.txt index 2ecff85..9d49217 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -101,6 +101,7 @@ if(CBK_WITH_MKL) bench_posv_compact bench_sysvnp_compact) add_executable(${prog} examples/${prog}.cpp) target_link_libraries(${prog} PRIVATE cbk MKL::Compact) + target_compile_features(${prog} PRIVATE cxx_std_17) # src/ for the internal header-only helpers the examples share with the # tests (cbk_matrix_view.hpp); the public API stays include/. target_include_directories(${prog} PRIVATE "${CMAKE_CURRENT_SOURCE_DIR}/src") @@ -111,7 +112,10 @@ if(CBK_WITH_MKL) endif() # --------------------------------------------------------------------------- -# Tests +# Tests. Like the examples above, each executable asks for C++17 itself: the +# library's cxx_std_17 is PRIVATE (a build requirement of the kernels, not of +# its C headers), so nothing propagates it, and AppleClang's default standard +# is older than gcc's and clang's on Linux. # --------------------------------------------------------------------------- if(CBK_BUILD_TESTS) enable_testing() @@ -145,6 +149,7 @@ if(CBK_BUILD_TESTS) foreach(routine geqrf ormqr orgqr potrf sytrfnp trsm gels) add_executable(test_${routine}_compact tests/test_${routine}_compact.cpp) target_link_libraries(test_${routine}_compact PRIVATE cbk LAPACKE::LAPACKE) + target_compile_features(test_${routine}_compact PRIVATE cxx_std_17) target_include_directories(test_${routine}_compact PRIVATE "${CMAKE_CURRENT_SOURCE_DIR}/src") add_test(NAME portable_${routine} COMMAND test_${routine}_compact) @@ -157,6 +162,7 @@ if(CBK_BUILD_TESTS) foreach(routine geqrf ormqr orgqr potrf sytrfnp trsm gels) add_executable(test_${routine}_mkl tests/test_${routine}_mkl.cpp) target_link_libraries(test_${routine}_mkl PRIVATE cbk MKL::Compact LAPACKE::LAPACKE) + target_compile_features(test_${routine}_mkl PRIVATE cxx_std_17) target_include_directories(test_${routine}_mkl PRIVATE "${CMAKE_CURRENT_SOURCE_DIR}/src") add_test(NAME mkl_${routine}_suites COMMAND test_${routine}_mkl) diff --git a/docs/building.md b/docs/building.md index 019dde3..03a36cf 100644 --- a/docs/building.md +++ b/docs/building.md @@ -31,7 +31,7 @@ On macOS, Accelerate provides LAPACK and CBLAS but no LAPACKE; [accelerate-lapacke](https://github.com/lepus2589/accelerate-lapacke) builds one on top of it (its configure needs a Fortran compiler, `brew install gcc`), and `-DCBK_TEST_LAPACK=accelerate` with `-DCMAKE_PREFIX_PATH` at its install -(`~/.local` by default) looks for that: the `lapacke` package plus +(its own prefix; the workflow caches it) looks for that: the `lapacke` package plus Accelerate's own CBLAS through ``. `.github/workflows/macos.yml` runs exactly this on the macOS runners, one job next to the Linux matrix. From c53a6ea8bf55d1a3574fc0526cdb067373845dec Mon Sep 17 00:00:00 2001 From: Claude Date: Mon, 14 Sep 2026 22:59:01 +0000 Subject: [PATCH 5/8] Simplify: FindLAPACKE on CMake's FindLAPACK, one C++ standard, review cleanups cmake/FindLAPACKE.cmake shrinks to what FindLAPACK does not know: with the MKL extension the stack is MKL::Compact's; otherwise find_package(LAPACK) picks the implementation under BLA_VENDOR (OpenBLAS, Generic, Apple) and the module adds lapacke.h, cblas.h (a distro's OpenBLAS subdirectory included) and liblapacke where LAPACKE is a library of its own, then runs one link test. An unset BLA_VENDOR means Generic, not FindLAPACK's "All": that took an installed MKL first and, under Netlib's liblapacke and OpenBLAS's cblas.h, put three LAPACK stacks in one binary (the suites failed that way). The CBK_TEST_LAPACK option goes; CI passes BLA_VENDOR from a single stack axis whose rows carry their apt packages, and the macOS job passes Apple. CMAKE_CXX_STANDARD 17 (no extensions) is set once for every target instead of a cxx_std_17 line per executable. From the review of the diff: ref_larfg and its dispatch entry were dead; lapack_ld is MatrixView::ld(); Staged stages only writable operands (ormqr's and gels's right-hand sides), so the const machinery goes; ref_gels factors A in place and runs ?gels on the copy; ACCELERATE_NEW_LAPACK lives on the target only; gen_matrix builds on gen_boosted; the geqrf invariants skip the ungated ?geqrf comparison on structured inputs and run per precision through one template; the gels and potrf suites keep one tolerance where three names held one value; sytrfnp's reference check runs one triangle through the transposed view. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_01GdNaquWMN3b39gvTsXZyjn --- .claude/CLAUDE.md | 8 +- .claude/hooks/session-start.sh | 6 +- .github/workflows/ci.yml | 53 ++---- .github/workflows/macos.yml | 12 +- CHANGELOG.md | 3 +- CMakeLists.txt | 48 ++--- cmake/FindLAPACKE.cmake | 326 +++++++++------------------------ docs/PLANS.md | 4 +- docs/building.md | 23 ++- tests/test_gels_compact.cpp | 19 +- tests/test_geqrf_compact.cpp | 81 ++++---- tests/test_lapack_util.hpp | 98 ++++------ tests/test_potrf_compact.cpp | 19 +- tests/test_sytrfnp_compact.cpp | 40 ++-- 14 files changed, 257 insertions(+), 483 deletions(-) diff --git a/.claude/CLAUDE.md b/.claude/CLAUDE.md index 709ba70..ff4a326 100644 --- a/.claude/CLAUDE.md +++ b/.claude/CLAUDE.md @@ -12,9 +12,9 @@ The default configure builds the portable library and its tests. The library needs nothing but CMake and a C++17 compiler; the tests validate against a real LAPACKE + CBLAS (`cmake/FindLAPACKE.cmake`: OpenBLAS or Netlib without MKL, `libopenblas-dev liblapacke-dev` / `liblapacke-dev liblapack-dev` on -Debian, MKL's own with the extension; `-DCBK_TEST_LAPACK=openblas|netlib|mkl` -pins one -- do that, since Debian's alternatives symlinks let either stand in -for the other): +Debian, MKL's own with the extension; `-DBLA_VENDOR=OpenBLAS|Generic` picks +one through CMake's FindLAPACK -- do that, since Debian's alternatives +symlinks let either stand in for the other): ```sh cmake -S . -B build -DCMAKE_BUILD_TYPE=Release @@ -188,7 +188,7 @@ workspace contract, or the benchmarks' threading. compilers, and Accelerate on macOS through `.github/workflows/macos.yml`): ```sh - CXX=clang++ cmake -S . -B build-clang -DCBK_TEST_LAPACK=openblas && cmake --build build-clang && ctest --test-dir build-clang + CXX=clang++ cmake -S . -B build-clang -DBLA_VENDOR=OpenBLAS && cmake --build build-clang && ctest --test-dir build-clang ``` - **Threading over groups.** Every all-groups driver is a call to diff --git a/.claude/hooks/session-start.sh b/.claude/hooks/session-start.sh index f14b17d..688ed82 100755 --- a/.claude/hooks/session-start.sh +++ b/.claude/hooks/session-start.sh @@ -11,8 +11,8 @@ # libopenblas-dev, liblapacke-dev, liblapack-dev # the two non-MKL LAPACKE stacks the test suites validate # against (cmake/FindLAPACKE.cmake), so a session can run -# the tests with -DCBK_TEST_LAPACK=openblas and =netlib as -# CI does. +# the tests with -DBLA_VENDOR=OpenBLAS and =Generic as CI +# does. # libomp-dev clang's OpenMP runtime and its omp.h, which clang-tidy # needs (.pre-commit-config.yaml explains). # pre-commit the formatting and linting driver, with its hook @@ -76,6 +76,6 @@ else fi echo "session-start: build with 'cmake -S . -B build -DCBK_WITH_MKL=ON && cmake --build build && ctest --test-dir build' (MKL is off by default)" -echo "session-start: the tests validate against LAPACKE; without MKL pick the stack with -DCBK_TEST_LAPACK=openblas|netlib (default: the first found)" +echo "session-start: the tests validate against LAPACKE; without MKL pick the stack with -DBLA_VENDOR=OpenBLAS|Generic (default: the first found)" echo "session-start: 'pre-commit run --all-files' checks the style; see .claude/CLAUDE.md" echo "session-start: for clang-tidy, 'CXX=clang++ cmake -S . -B build-tidy' then 'pre-commit run --hook-stage manual clang-tidy --all-files'" diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 39c883d..f351021 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -12,66 +12,47 @@ concurrency: jobs: build-and-test: - name: ${{ matrix.cxx }} / ${{ matrix.mkl == 'ON' && 'MKL' || matrix.lapack }}${{ matrix.interface == 'ilp64' && ' (ilp64)' || '' }} + name: ${{ matrix.cxx }} / ${{ matrix.stack }}${{ matrix.interface == 'ilp64' && ' (ilp64)' || '' }} runs-on: ubuntu-latest strategy: fail-fast: false matrix: cc: [gcc, clang] - # The tests validate against a real LAPACKE + CBLAS: MKL's own with - # the MKL extension, OpenBLAS or Netlib's reference LAPACK without it - # (cmake/FindLAPACKE.cmake; -DCBK_TEST_LAPACK pins the stack, so a - # job cannot silently pick up the other one through Debian's - # update-alternatives symlinks). - mkl: ["ON", "OFF"] - lapack: [openblas, netlib] + # The LAPACKE + CBLAS stack the tests validate against, and the + # packages that provide it: MKL's own with the MKL extension (which + # libmkl-dev also supplies the Compact API for), else OpenBLAS or the + # reference LAPACK through CMake's FindLAPACK, picked by BLA_VENDOR + # (cmake/FindLAPACKE.cmake). Debian ships OpenBLAS's LAPACKE as + # Netlib's liblapacke, hence that package in both rows. + stack: [MKL, OpenBLAS, Generic] interface: [lp64] - exclude: - - { mkl: "ON", lapack: netlib } include: - { cc: gcc, cxx: g++ } - { cc: clang, cxx: clang++ } + - { stack: MKL, packages: libmkl-dev } + - { stack: OpenBLAS, packages: libopenblas-dev liblapacke-dev } + - { stack: Generic, packages: liblapacke-dev liblapack-dev libblas-dev } # MKL's 64-bit integer interface: MKL_INT is part of the MKL-style # API's ABI, so the wrappers and their suites must build and pass - # against it too (libmkl-dev ships both interface layers). The - # LAPACKE the suites call is MKL's, at the same interface. - - { cc: gcc, cxx: g++, mkl: "ON", lapack: openblas, interface: ilp64 } + # against it too (libmkl-dev ships both interface layers). + - { cc: gcc, cxx: g++, stack: MKL, packages: libmkl-dev, interface: ilp64 } env: CC: ${{ matrix.cc }} CXX: ${{ matrix.cxx }} steps: - uses: actions/checkout@v4 - # libmkl-dev supplies the Compact extension and MKL's LAPACK/LAPACKE, so - # the MKL-backed suites validate against real MKL + LAPACK, and the - # portable suites against MKL's LAPACKE. - - name: Install Intel MKL - if: matrix.mkl == 'ON' + - name: Install ${{ matrix.stack }} run: | sudo apt-get update - sudo apt-get install -y --no-install-recommends libmkl-dev - - # Without MKL the suites need a LAPACKE stack: OpenBLAS (Debian ships - # its LAPACKE as Netlib's liblapacke, hence both packages), or Netlib's - # reference LAPACK/BLAS with its LAPACKE. - - name: Install OpenBLAS - if: matrix.mkl == 'OFF' && matrix.lapack == 'openblas' - run: | - sudo apt-get update - sudo apt-get install -y --no-install-recommends libopenblas-dev liblapacke-dev - - - name: Install Netlib LAPACK - if: matrix.mkl == 'OFF' && matrix.lapack == 'netlib' - run: | - sudo apt-get update - sudo apt-get install -y --no-install-recommends liblapacke-dev liblapack-dev libblas-dev + sudo apt-get install -y --no-install-recommends ${{ matrix.packages }} - name: Configure run: > cmake -S . -B build -G Ninja -DCMAKE_BUILD_TYPE=Release - -DCBK_WITH_MKL=${{ matrix.mkl }} - -DCBK_TEST_LAPACK=${{ matrix.mkl == 'ON' && 'mkl' || matrix.lapack }} + -DCBK_WITH_MKL=${{ matrix.stack == 'MKL' && 'ON' || 'OFF' }} + -DBLA_VENDOR=${{ matrix.stack }} -DMKLCompact_INTERFACE=${{ matrix.interface }} - name: Build diff --git a/.github/workflows/macos.yml b/.github/workflows/macos.yml index bd05ddd..bebae62 100644 --- a/.github/workflows/macos.yml +++ b/.github/workflows/macos.yml @@ -5,8 +5,9 @@ name: macOS # accelerate-lapacke (https://github.com/lepus2589/accelerate-lapacke, Netlib's # LAPACKE over Accelerate's LAPACK interface; its configure needs a Fortran # compiler, which Homebrew's gcc supplies) into its own prefix, then -# configures cbk with -DCBK_TEST_LAPACK=accelerate (cmake/FindLAPACKE.cmake: -# the `lapacke` package plus from the SDK). macOS runners +# configures cbk with -DBLA_VENDOR=Apple (cmake/FindLAPACKE.cmake: Accelerate +# from CMake's FindLAPACK, lapacke.h and liblapacke from that prefix, +# from the SDK). macOS runners # cost several times a Linux minute, so this is one job, not a matrix, and # the LAPACKE build (several minutes, most of it configuring Netlib LAPACK's # Fortran tree) is cached: the install prefix is restored by actions/cache @@ -39,9 +40,10 @@ jobs: - name: Identify the SDK id: sdk run: | + sdk=$(xcrun --show-sdk-version) echo "path=$(xcrun --show-sdk-path)" >> "$GITHUB_OUTPUT" - echo "version=$(xcrun --show-sdk-version)" >> "$GITHUB_OUTPUT" - echo "macOS $(sw_vers -productVersion), SDK $(xcrun --show-sdk-version)" + echo "version=$sdk" >> "$GITHUB_OUTPUT" + echo "macOS $(sw_vers -productVersion), SDK $sdk" - name: Cache accelerate-lapacke id: lapacke-cache @@ -74,7 +76,7 @@ jobs: run: > cmake -S . -B build -DCMAKE_BUILD_TYPE=Release - -DCBK_TEST_LAPACK=accelerate + -DBLA_VENDOR=Apple -DCMAKE_PREFIX_PATH="$ACCELERATE_LAPACKE_PREFIX" - name: Build diff --git a/CHANGELOG.md b/CHANGELOG.md index 4b3b1f7..651a2f8 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,7 +11,8 @@ release may change the API. - Tests: every suite validates against a real LAPACKE + CBLAS stack (OpenBLAS, Netlib, or MKL's own with the MKL extension; - `cmake/FindLAPACKE.cmake`, `-DCBK_TEST_LAPACK`) instead of hand-rolled + `cmake/FindLAPACKE.cmake` over `find_package(LAPACK)` and `BLA_VENDOR`) + instead of hand-rolled scalar references (issue #27). The dense-LAPACK cross-checks of the MKL suites run in the portable suites, on any stack; the unpivoted LDL^T keeps its scalar reference, itself validated against CBLAS. Building the tests diff --git a/CMakeLists.txt b/CMakeLists.txt index 9d49217..1d1a5ef 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -16,13 +16,6 @@ set(CMAKE_EXPORT_COMPILE_COMMANDS ON) option(CBK_WITH_MKL "Build the Intel MKL compact extension, examples and MKL tests" OFF) option(CBK_WITH_OPENMP "Thread each routine's loop over compact groups with OpenMP" ON) option(CBK_BUILD_TESTS "Build tests and register them with CTest" ON) -# The tests validate against a real LAPACKE + CBLAS (the library itself needs -# none): which stack, as cmake/FindLAPACKE.cmake's LAPACKE_VENDOR. `auto` -# takes MKL's with the MKL extension and otherwise the first of OpenBLAS, -# Netlib, any `lapacke` + find_package(LAPACK) pair; the other values pin one. -set(CBK_TEST_LAPACK "auto" CACHE STRING - "LAPACKE/CBLAS stack the tests validate against: auto | mkl | openblas | netlib | accelerate") -set_property(CACHE CBK_TEST_LAPACK PROPERTY STRINGS auto mkl openblas netlib accelerate) # The install rules belong to cbk's own build, not to a project that # add_subdirectory()s it. option(CBK_INSTALL "Generate the install target" ${PROJECT_IS_TOP_LEVEL}) @@ -34,10 +27,16 @@ if(NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES) endif() # C++17, no GNU extensions (the kernels' vector_size attribute compiles under -# strict -std=c++17). Optimization flags are the caller's to choose; for a -# performance build pass them through CMAKE_CXX_FLAGS, e.g. +# strict -std=c++17), for every target here: the library's own requirement +# below is PRIVATE (its headers are C), so nothing would pass it to the +# tests, examples and benchmarks, and AppleClang's default standard is older +# than gcc's and clang's on Linux. Optimization flags are the caller's to +# choose; for a performance build pass them through CMAKE_CXX_FLAGS, e.g. # cmake -S . -B build -DCMAKE_CXX_FLAGS="-O3 -march=native" # so the SIMD kernels target the host's widest vectors. +set(CMAKE_CXX_STANDARD 17) +set(CMAKE_CXX_STANDARD_REQUIRED ON) +set(CMAKE_CXX_EXTENSIONS OFF) # OpenMP threads the group loop inside every routine (see the threading note # in src/cbk_common.hpp); the benchmarks use it for their outer loops. @@ -101,7 +100,6 @@ if(CBK_WITH_MKL) bench_posv_compact bench_sysvnp_compact) add_executable(${prog} examples/${prog}.cpp) target_link_libraries(${prog} PRIVATE cbk MKL::Compact) - target_compile_features(${prog} PRIVATE cxx_std_17) # src/ for the internal header-only helpers the examples share with the # tests (cbk_matrix_view.hpp); the public API stays include/. target_include_directories(${prog} PRIVATE "${CMAKE_CURRENT_SOURCE_DIR}/src") @@ -112,35 +110,15 @@ if(CBK_WITH_MKL) endif() # --------------------------------------------------------------------------- -# Tests. Like the examples above, each executable asks for C++17 itself: the -# library's cxx_std_17 is PRIVATE (a build requirement of the kernels, not of -# its C headers), so nothing propagates it, and AppleClang's default standard -# is older than gcc's and clang's on Linux. +# Tests # --------------------------------------------------------------------------- if(CBK_BUILD_TESTS) enable_testing() # The dense references of every suite come from a real LAPACKE + CBLAS - # (tests/test_lapack_util.hpp): with the MKL extension it must be MKL's own - # (one LAPACK stack per binary), otherwise any of OpenBLAS, Netlib, ... - if(CBK_WITH_MKL) - if(NOT CBK_TEST_LAPACK MATCHES "^(auto|mkl)$") - message(FATAL_ERROR "CBK_TEST_LAPACK=${CBK_TEST_LAPACK} cannot be combined with CBK_WITH_MKL=ON: the MKL-backed suites link MKL, so the dense references must be MKL's LAPACKE too") - endif() - set(LAPACKE_VENDOR MKL) - elseif(CBK_TEST_LAPACK STREQUAL "mkl") - set(LAPACKE_VENDOR MKL) - elseif(CBK_TEST_LAPACK STREQUAL "openblas") - set(LAPACKE_VENDOR OpenBLAS) - elseif(CBK_TEST_LAPACK STREQUAL "netlib") - set(LAPACKE_VENDOR Netlib) - elseif(CBK_TEST_LAPACK STREQUAL "accelerate") - set(LAPACKE_VENDOR Accelerate) - elseif(CBK_TEST_LAPACK STREQUAL "auto") - set(LAPACKE_VENDOR Any) - else() - message(FATAL_ERROR "CBK_TEST_LAPACK must be auto, mkl, openblas, netlib or accelerate (got '${CBK_TEST_LAPACK}')") - endif() + # (tests/test_lapack_util.hpp): MKL's own with the MKL extension, otherwise + # whatever CMake's FindLAPACK picks under BLA_VENDOR (OpenBLAS, Generic for + # Netlib, Apple) plus its LAPACKE (cmake/FindLAPACKE.cmake). find_package(LAPACKE REQUIRED) # The portable C API's suites: each kernel vs LAPACKE/CBLAS (its unblocked @@ -149,7 +127,6 @@ if(CBK_BUILD_TESTS) foreach(routine geqrf ormqr orgqr potrf sytrfnp trsm gels) add_executable(test_${routine}_compact tests/test_${routine}_compact.cpp) target_link_libraries(test_${routine}_compact PRIVATE cbk LAPACKE::LAPACKE) - target_compile_features(test_${routine}_compact PRIVATE cxx_std_17) target_include_directories(test_${routine}_compact PRIVATE "${CMAKE_CURRENT_SOURCE_DIR}/src") add_test(NAME portable_${routine} COMMAND test_${routine}_compact) @@ -162,7 +139,6 @@ if(CBK_BUILD_TESTS) foreach(routine geqrf ormqr orgqr potrf sytrfnp trsm gels) add_executable(test_${routine}_mkl tests/test_${routine}_mkl.cpp) target_link_libraries(test_${routine}_mkl PRIVATE cbk MKL::Compact LAPACKE::LAPACKE) - target_compile_features(test_${routine}_mkl PRIVATE cxx_std_17) target_include_directories(test_${routine}_mkl PRIVATE "${CMAKE_CURRENT_SOURCE_DIR}/src") add_test(NAME mkl_${routine}_suites COMMAND test_${routine}_mkl) diff --git a/cmake/FindLAPACKE.cmake b/cmake/FindLAPACKE.cmake index 4f27e74..24089dd 100644 --- a/cmake/FindLAPACKE.cmake +++ b/cmake/FindLAPACKE.cmake @@ -1,268 +1,124 @@ # FindLAPACKE.cmake # -# Locate a LAPACKE + CBLAS stack for the test suites and expose it as +# The LAPACKE + CBLAS stack the test suites validate against, as one target: # -# LAPACKE::LAPACKE -- the include dirs, the link line, and one compile -# definition naming the stack's headers: -# CBK_LAPACKE_MKL , -# CBK_LAPACKE_ACCELERATE , -# (none) , +# LAPACKE::LAPACKE include dirs, link line, and the definition that tells +# tests/test_lapack_util.hpp which headers to include +# (CBK_LAPACKE_MKL, CBK_LAPACKE_ACCELERATE, or none for +# the plain / ). # -# LAPACKE_VENDOR Any (default) | MKL | OpenBLAS | Netlib | Accelerate -# which stack to look for; Any tries OpenBLAS, then Netlib, -# then whatever `lapacke` library and find_package(LAPACK) -# turn up. -# LAPACKE_VENDOR_FOUND the stack that was found, as one of the names above. +# With the MKL extension (MKL::Compact defined by FindMKLCompact.cmake) the +# stack is MKL's own -- a test binary must never link two LAPACK +# implementations. Otherwise CMake's FindLAPACK picks the implementation, +# steered by BLA_VENDOR as usual (OpenBLAS, Generic for the distro's plain +# liblapack, Apple for Accelerate). An unset BLA_VENDOR means Generic here, +# not FindLAPACK's "All": that would take an installed MKL first, and MKL's +# LAPACK under Netlib's liblapacke and OpenBLAS's cblas.h is three stacks in +# one binary (measured: the suites fail that way). This module then adds +# what FindLAPACK does not know about: # -# Every candidate is link-tested with a program that calls LAPACKE_dgeqrf, -# LAPACKE_dgeqr2, cblas_dgemm and cblas_dtrsm, so a header without its -# library (or the other way round) is rejected rather than found. +# lapacke.h and cblas.h in the default include dirs or a distro's +# OpenBLAS subdirectory (openblas/, openblas-pthread/, +# ...); with Accelerate, which has CBLAS but no +# LAPACKE, lapacke.h comes from an installed +# accelerate-lapacke (CMAKE_PREFIX_PATH at its +# prefix; https://github.com/lepus2589/accelerate-lapacke) +# and CBLAS from the SDK's vecLib sub-framework, whose +# framework directory xcrun locates. +# liblapacke where LAPACKE is a library of its own: Netlib's, +# which Debian also uses for its OpenBLAS, and +# accelerate-lapacke's. A full OpenBLAS carries the +# symbols itself and has no such library. # -# MKL reuses FindMKLCompact.cmake's MKL::Compact (the whole MKL link -# line, which carries LAPACKE and CBLAS). With -DCBK_WITH_MKL=ON -# this is the only acceptable stack: a test binary must never -# link two LAPACK implementations. -# OpenBLAS pkg-config `openblas` (the distro's include subdirectory -- -# Debian's openblas-pthread/, Fedora's openblas/ -- comes from -# there), else cblas.h + libopenblas by hand. Debian and Ubuntu -# strip LAPACKE out of libopenblas and ship Netlib's -# liblapacke separately, so when the link test fails on -# libopenblas alone it is retried with liblapacke added. -# Netlib liblapacke + liblapack + libblas. On Debian the plain -# liblapack.so / libblas.so are update-alternatives symlinks -# that usually point at OpenBLAS, so the reference libraries -# are looked for under the lapack/ and blas/ subdirectories -# first (PATH_SUFFIXES), which is where the distro keeps them. -# Accelerate macOS: Apple's Accelerate framework provides LAPACK and -# CBLAS but no LAPACKE; https://github.com/lepus2589/accelerate-lapacke -# builds one on top of it (Netlib's LAPACKE, so it needs a -# Fortran compiler to configure) and installs a `lapacke` -# CMake package. .github/workflows/macos.yml builds it on the -# macOS runners. -# -# ILP64 is supported through MKL only (MKL_ILP64 comes with MKL::Compact); -# the other stacks are used at their default 32-bit integer interface. +# One link test of LAPACKE_dgeqrf, LAPACKE_dgeqr2, cblas_dgemm and cblas_dtrsm +# accepts the result. ILP64 comes through MKL only. # # Assisted-by: Claude:claude-fable-5 -set(LAPACKE_VENDOR "Any" CACHE STRING - "LAPACKE/CBLAS stack for the tests: Any | MKL | OpenBLAS | Netlib | Accelerate") -set_property(CACHE LAPACKE_VENDOR PROPERTY STRINGS Any MKL OpenBLAS Netlib Accelerate) -if(NOT LAPACKE_VENDOR MATCHES "^(Any|MKL|OpenBLAS|Netlib|Accelerate)$") - message(FATAL_ERROR "LAPACKE_VENDOR must be Any, MKL, OpenBLAS, Netlib or Accelerate (got '${LAPACKE_VENDOR}')") -endif() - include(CMakePushCheckState) include(CheckCXXSourceCompiles) -# The link test. `defs` selects the header set exactly as the test util does. -function(_lapacke_link_test vendor includes libs defs outvar) - cmake_push_check_state(RESET) - set(CMAKE_REQUIRED_INCLUDES ${includes}) - set(CMAKE_REQUIRED_LIBRARIES ${libs}) - set(CMAKE_REQUIRED_DEFINITIONS ${defs}) - set(CMAKE_REQUIRED_QUIET TRUE) - check_cxx_source_compiles(" - #if defined(CBK_LAPACKE_MKL) - #include - #include - #elif defined(CBK_LAPACKE_ACCELERATE) - #include - #include - #else - #include - #include - #endif - int main() { - double a[4] = {2, 0, 0, 2}, tau[2], c[4] = {1, 0, 0, 1}; - lapack_int info = LAPACKE_dgeqrf(LAPACK_COL_MAJOR, 2, 2, a, 2, tau); - info += LAPACKE_dgeqr2(LAPACK_COL_MAJOR, 2, 2, a, 2, tau); - cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, 2, 2, 2, 1.0, a, 2, c, 2, 0.0, c, 2); - cblas_dtrsm(CblasColMajor, CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, 2, 2, 1.0, a, 2, c, 2); - return (int)info; - } - " LAPACKE_LINKS_${vendor}) - cmake_pop_check_state() - set(${outvar} ${LAPACKE_LINKS_${vendor}} PARENT_SCOPE) -endfunction() - -set(_lapacke_found_vendor "") set(_lapacke_includes "") set(_lapacke_libs "") set(_lapacke_defs "") set(_lapacke_options "") -# ---- MKL ------------------------------------------------------------------- -if(NOT _lapacke_found_vendor AND LAPACKE_VENDOR STREQUAL "MKL") - find_package(MKLCompact REQUIRED) - _lapacke_link_test(MKL "${MKLCompact_INCLUDE_DIR}" "MKL::Compact" "-DCBK_LAPACKE_MKL" _ok) - if(_ok) - set(_lapacke_found_vendor MKL) - set(_lapacke_includes "") # MKL::Compact carries them - set(_lapacke_libs MKL::Compact) - set(_lapacke_defs CBK_LAPACKE_MKL) - endif() -endif() - -# ---- OpenBLAS --------------------------------------------------------------- -if(NOT _lapacke_found_vendor AND LAPACKE_VENDOR MATCHES "^(Any|OpenBLAS)$") - find_package(PkgConfig QUIET) - set(_ob_includes "") - set(_ob_libs "") - if(PKG_CONFIG_FOUND) - pkg_check_modules(_ob QUIET openblas) - if(_ob_FOUND) - set(_ob_includes ${_ob_INCLUDE_DIRS}) - set(_ob_libs ${_ob_LINK_LIBRARIES}) - endif() - endif() - if(NOT _ob_libs) - find_library(LAPACKE_OpenBLAS_LIBRARY NAMES openblas - PATH_SUFFIXES openblas openblas-pthread openblas-openmp openblas-serial) - find_path(LAPACKE_OpenBLAS_INCLUDE_DIR NAMES openblas_config.h - PATH_SUFFIXES openblas openblas-pthread openblas-openmp openblas-serial) - mark_as_advanced(LAPACKE_OpenBLAS_LIBRARY LAPACKE_OpenBLAS_INCLUDE_DIR) - if(LAPACKE_OpenBLAS_LIBRARY AND LAPACKE_OpenBLAS_INCLUDE_DIR) - set(_ob_libs "${LAPACKE_OpenBLAS_LIBRARY}") - set(_ob_includes "${LAPACKE_OpenBLAS_INCLUDE_DIR}") - endif() - endif() - if(_ob_libs) - # lapacke.h sits next to cblas.h in a full OpenBLAS install, or in the - # default include dir when the distro ships Netlib's LAPACKE separately. - find_path(LAPACKE_OpenBLAS_LAPACKE_INCLUDE_DIR NAMES lapacke.h - HINTS ${_ob_includes} - PATH_SUFFIXES openblas openblas-pthread openblas-openmp openblas-serial) - mark_as_advanced(LAPACKE_OpenBLAS_LAPACKE_INCLUDE_DIR) - if(LAPACKE_OpenBLAS_LAPACKE_INCLUDE_DIR) - list(APPEND _ob_includes "${LAPACKE_OpenBLAS_LAPACKE_INCLUDE_DIR}") - list(REMOVE_DUPLICATES _ob_includes) - _lapacke_link_test(OpenBLAS "${_ob_includes}" "${_ob_libs}" "" _ok) - if(NOT _ok) - # Debian/Ubuntu: LAPACKE lives in Netlib's liblapacke, which calls - # back into whichever LAPACK is loaded -- libopenblas here, since it - # is linked first and exports the LAPACK symbols itself. - find_library(LAPACKE_OpenBLAS_LAPACKE_LIBRARY NAMES lapacke) - mark_as_advanced(LAPACKE_OpenBLAS_LAPACKE_LIBRARY) - if(LAPACKE_OpenBLAS_LAPACKE_LIBRARY) - list(PREPEND _ob_libs "${LAPACKE_OpenBLAS_LAPACKE_LIBRARY}") - _lapacke_link_test(OpenBLAS_lapacke "${_ob_includes}" "${_ob_libs}" "" _ok) - endif() - endif() - if(_ok) - set(_lapacke_found_vendor OpenBLAS) - set(_lapacke_includes ${_ob_includes}) - set(_lapacke_libs ${_ob_libs}) - endif() - endif() - endif() -endif() - -# ---- Netlib ----------------------------------------------------------------- -if(NOT _lapacke_found_vendor AND LAPACKE_VENDOR MATCHES "^(Any|Netlib)$") - find_path(LAPACKE_Netlib_INCLUDE_DIR NAMES lapacke.h PATH_SUFFIXES lapacke) - find_path(LAPACKE_Netlib_CBLAS_INCLUDE_DIR NAMES cblas.h PATH_SUFFIXES cblas) - find_library(LAPACKE_Netlib_LAPACKE_LIBRARY NAMES lapacke) - # The reference libraries in their own subdirectories first (Debian), then - # the plain names. - find_library(LAPACKE_Netlib_LAPACK_LIBRARY NAMES lapack PATH_SUFFIXES lapack) - find_library(LAPACKE_Netlib_BLAS_LIBRARY NAMES blas PATH_SUFFIXES blas) - mark_as_advanced(LAPACKE_Netlib_INCLUDE_DIR LAPACKE_Netlib_CBLAS_INCLUDE_DIR - LAPACKE_Netlib_LAPACKE_LIBRARY LAPACKE_Netlib_LAPACK_LIBRARY LAPACKE_Netlib_BLAS_LIBRARY) - if(LAPACKE_Netlib_INCLUDE_DIR AND LAPACKE_Netlib_CBLAS_INCLUDE_DIR - AND LAPACKE_Netlib_LAPACKE_LIBRARY AND LAPACKE_Netlib_LAPACK_LIBRARY - AND LAPACKE_Netlib_BLAS_LIBRARY) - set(_nl_includes "${LAPACKE_Netlib_INCLUDE_DIR}" "${LAPACKE_Netlib_CBLAS_INCLUDE_DIR}") - list(REMOVE_DUPLICATES _nl_includes) - set(_nl_libs "${LAPACKE_Netlib_LAPACKE_LIBRARY}" "${LAPACKE_Netlib_LAPACK_LIBRARY}" - "${LAPACKE_Netlib_BLAS_LIBRARY}") - # A test binary calls LAPACKE and CBLAS, never LAPACK itself, so GNU ld's - # --as-needed (Debian's default) would drop liblapack from its NEEDED - # list and liblapacke would then load `liblapack.so.3` through the - # loader cache -- the alternatives symlink, i.e. OpenBLAS again. Keeping - # it a direct dependency makes the loader take the reference library - # from the RUNPATH CMake sets for the build tree, and liblapacke's own - # request for the same soname is then satisfied by the copy already - # loaded. - if(CMAKE_SYSTEM_NAME STREQUAL "Linux" AND NOT MSVC) - list(PREPEND _nl_libs "-Wl,--no-as-needed") - endif() - _lapacke_link_test(Netlib "${_nl_includes}" "${_nl_libs}" "" _ok) - if(_ok) - set(_lapacke_found_vendor Netlib) - set(_lapacke_includes ${_nl_includes}) - set(_lapacke_libs ${_nl_libs}) - endif() +if(TARGET MKL::Compact) + set(_lapacke_libs MKL::Compact) + set(_lapacke_defs CBK_LAPACKE_MKL) + set(LAPACKE_VENDOR_FOUND "MKL") +else() + if(NOT BLA_VENDOR) + set(BLA_VENDOR Generic) endif() -endif() - -# ---- Accelerate (macOS) ----------------------------------------------------- -# accelerate-lapacke installs Netlib's lapacke-config.cmake (looked up under -# its own lower-case file name, so this module is not re-entered) with the -# `lapacke` target; CBLAS is Accelerate's, reached as through -# the framework directory of the SDK's Accelerate.framework, which -# `xcrun --show-sdk-path` locates. Only the vecLib sub-framework's header is -# included, never the umbrella: its LAPACK -# prototypes and the ones lapacke.h declares are the same C symbols with -# different parameter types. -if(NOT _lapacke_found_vendor AND LAPACKE_VENDOR STREQUAL "Accelerate") - find_package(lapacke CONFIG QUIET) - find_library(LAPACKE_Accelerate_FRAMEWORK Accelerate) - mark_as_advanced(LAPACKE_Accelerate_FRAMEWORK) - execute_process(COMMAND xcrun --show-sdk-path - OUTPUT_VARIABLE _ac_sdk OUTPUT_STRIP_TRAILING_WHITESPACE ERROR_QUIET) - set(_ac_frameworks "${_ac_sdk}/System/Library/Frameworks/Accelerate.framework/Frameworks") - if(TARGET lapacke AND LAPACKE_Accelerate_FRAMEWORK AND EXISTS "${_ac_frameworks}/vecLib.framework") - set(_ac_libs lapacke "${LAPACKE_Accelerate_FRAMEWORK}") - set(_ac_defs -DCBK_LAPACKE_ACCELERATE -DACCELERATE_NEW_LAPACK "-F${_ac_frameworks}") - _lapacke_link_test(Accelerate "" "${_ac_libs}" "${_ac_defs}" _ok) - if(_ok) - set(_lapacke_found_vendor Accelerate) - set(_lapacke_libs ${_ac_libs}) - set(_lapacke_defs CBK_LAPACKE_ACCELERATE ACCELERATE_NEW_LAPACK) - set(_lapacke_options "-F${_ac_frameworks}") - endif() + find_package(LAPACK REQUIRED) + set(_openblas_dirs openblas openblas-pthread openblas-openmp openblas-serial) + find_path(LAPACKE_INCLUDE_DIR NAMES lapacke.h PATH_SUFFIXES ${_openblas_dirs}) + find_library(LAPACKE_LIBRARY NAMES lapacke) + mark_as_advanced(LAPACKE_INCLUDE_DIR LAPACKE_LIBRARY) + set(_lapacke_includes ${LAPACKE_INCLUDE_DIR}) + set(_lapacke_libs LAPACK::LAPACK) + if(LAPACKE_LIBRARY) + list(PREPEND _lapacke_libs "${LAPACKE_LIBRARY}") endif() -endif() - -# ---- Any other lapacke + LAPACK/BLAS pair ----------------------------------- -if(NOT _lapacke_found_vendor AND LAPACKE_VENDOR STREQUAL "Any") - find_path(LAPACKE_Generic_INCLUDE_DIR NAMES lapacke.h) - find_path(LAPACKE_Generic_CBLAS_INCLUDE_DIR NAMES cblas.h) - find_library(LAPACKE_Generic_LAPACKE_LIBRARY NAMES lapacke) - mark_as_advanced(LAPACKE_Generic_INCLUDE_DIR LAPACKE_Generic_CBLAS_INCLUDE_DIR - LAPACKE_Generic_LAPACKE_LIBRARY) - find_package(LAPACK QUIET) - if(LAPACKE_Generic_INCLUDE_DIR AND LAPACKE_Generic_CBLAS_INCLUDE_DIR - AND LAPACKE_Generic_LAPACKE_LIBRARY AND LAPACK_FOUND) - set(_ge_includes "${LAPACKE_Generic_INCLUDE_DIR}" "${LAPACKE_Generic_CBLAS_INCLUDE_DIR}") - list(REMOVE_DUPLICATES _ge_includes) - set(_ge_libs "${LAPACKE_Generic_LAPACKE_LIBRARY}" ${LAPACK_LIBRARIES}) - _lapacke_link_test(Generic "${_ge_includes}" "${_ge_libs}" "" _ok) - if(_ok) - set(_lapacke_found_vendor Generic) - set(_lapacke_includes ${_ge_includes}) - set(_lapacke_libs ${_ge_libs}) - endif() + if(APPLE AND LAPACK_LIBRARIES MATCHES "Accelerate") + execute_process(COMMAND xcrun --show-sdk-path + OUTPUT_VARIABLE _sdk OUTPUT_STRIP_TRAILING_WHITESPACE ERROR_QUIET) + set(_lapacke_options "-F${_sdk}/System/Library/Frameworks/Accelerate.framework/Frameworks") + set(_lapacke_defs CBK_LAPACKE_ACCELERATE ACCELERATE_NEW_LAPACK) + else() + find_path(CBLAS_INCLUDE_DIR NAMES cblas.h PATH_SUFFIXES ${_openblas_dirs}) + mark_as_advanced(CBLAS_INCLUDE_DIR) + list(APPEND _lapacke_includes ${CBLAS_INCLUDE_DIR}) endif() + set(LAPACKE_VENDOR_FOUND "${BLA_VENDOR}") endif() -set(LAPACKE_VENDOR_FOUND "${_lapacke_found_vendor}") -set(LAPACKE_INCLUDE_DIRS "${_lapacke_includes}") -set(LAPACKE_LIBRARIES "${_lapacke_libs}") +# The link test, against exactly what the target will carry (cached per stack). +string(MAKE_C_IDENTIFIER "LAPACKE_LINK_TEST_${LAPACKE_VENDOR_FOUND}" _lapacke_check) +cmake_push_check_state(RESET) +set(CMAKE_REQUIRED_INCLUDES ${_lapacke_includes}) +set(CMAKE_REQUIRED_LIBRARIES ${_lapacke_libs}) +set(CMAKE_REQUIRED_FLAGS "${_lapacke_options}") +foreach(_d IN LISTS _lapacke_defs) + list(APPEND CMAKE_REQUIRED_DEFINITIONS "-D${_d}") +endforeach() +set(CMAKE_REQUIRED_QUIET TRUE) +check_cxx_source_compiles(" +#if defined(CBK_LAPACKE_MKL) +#include +#include +#elif defined(CBK_LAPACKE_ACCELERATE) +#include +#include +#else +#include +#include +#endif +int main() { + double a[4] = {2, 0, 0, 2}, tau[2], c[4] = {1, 0, 0, 1}; + lapack_int info = LAPACKE_dgeqrf(LAPACK_COL_MAJOR, 2, 2, a, 2, tau); + info += LAPACKE_dgeqr2(LAPACK_COL_MAJOR, 2, 2, a, 2, tau); + cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, 2, 2, 2, 1.0, a, 2, c, 2, 0.0, c, 2); + cblas_dtrsm(CblasColMajor, CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, 2, 2, 1.0, a, 2, c, 2); + return (int)info; +} +" ${_lapacke_check}) +cmake_pop_check_state() +set(LAPACKE_LINKS ${${_lapacke_check}}) include(FindPackageHandleStandardArgs) find_package_handle_standard_args(LAPACKE - REQUIRED_VARS LAPACKE_VENDOR_FOUND + REQUIRED_VARS LAPACKE_LINKS REASON_FAILURE_MESSAGE - "No LAPACKE + CBLAS stack for LAPACKE_VENDOR=${LAPACKE_VENDOR}. Install one (Debian/Ubuntu: libopenblas-dev, or liblapacke-dev with liblapack-dev, or libmkl-dev with -DCBK_WITH_MKL=ON) or point CMAKE_PREFIX_PATH at it. The tests need one, the library itself does not (-DCBK_BUILD_TESTS=OFF).") + "no LAPACKE + CBLAS next to the LAPACK found. Debian/Ubuntu: libopenblas-dev liblapacke-dev, or liblapacke-dev liblapack-dev libblas-dev with -DBLA_VENDOR=Generic. macOS: an installed accelerate-lapacke on CMAKE_PREFIX_PATH with -DBLA_VENDOR=Apple. Or -DCBK_WITH_MKL=ON for MKL's. The tests need one, the library does not (-DCBK_BUILD_TESTS=OFF).") if(LAPACKE_FOUND AND NOT TARGET LAPACKE::LAPACKE) add_library(LAPACKE::LAPACKE INTERFACE IMPORTED) set_target_properties(LAPACKE::LAPACKE PROPERTIES - INTERFACE_INCLUDE_DIRECTORIES "${LAPACKE_INCLUDE_DIRS}" - INTERFACE_LINK_LIBRARIES "${LAPACKE_LIBRARIES}" + INTERFACE_INCLUDE_DIRECTORIES "${_lapacke_includes}" + INTERFACE_LINK_LIBRARIES "${_lapacke_libs}" INTERFACE_COMPILE_DEFINITIONS "${_lapacke_defs}" INTERFACE_COMPILE_OPTIONS "${_lapacke_options}") - message(STATUS "LAPACKE for the tests: ${LAPACKE_VENDOR_FOUND} (${LAPACKE_LIBRARIES})") + message(STATUS "LAPACKE for the tests: ${LAPACKE_VENDOR_FOUND} (${_lapacke_libs})") endif() diff --git a/docs/PLANS.md b/docs/PLANS.md index 17ddbf6..bea94e0 100644 --- a/docs/PLANS.md +++ b/docs/PLANS.md @@ -100,8 +100,8 @@ listed in `examples.md`. LDL^T (`ref_sytf2np`, which LAPACK lacks) is the one hand-rolled reference left, validated against `cblas_?trmm` in its suite. `cmake/FindLAPACKE.cmake` takes MKL's LAPACKE with the MKL extension and - OpenBLAS, Netlib or any `lapacke` + `find_package(LAPACK)` pair without - (`-DCBK_TEST_LAPACK`); CI runs OpenBLAS and Netlib under gcc and clang, and + otherwise `find_package(LAPACK)`'s pick under `BLA_VENDOR` plus its + LAPACKE; CI runs OpenBLAS and Netlib under gcc and clang, and Accelerate + accelerate-lapacke on the macOS runners (`.github/workflows/macos.yml`, written without a macOS machine at hand). The dense cross-checks of the `geqrf` and `potrf` MKL suites (design 7.1) diff --git a/docs/building.md b/docs/building.md index 03a36cf..18331ba 100644 --- a/docs/building.md +++ b/docs/building.md @@ -14,12 +14,14 @@ configure, which builds them, also needs one of * Intel MKL's own, with the MKL extension below (then it is the only stack the tests may link). -`cmake/FindLAPACKE.cmake` picks the first it finds (OpenBLAS, then Netlib, -then any `lapacke` library next to what `find_package(LAPACK)` returns) and -reports it as `LAPACKE for the tests: ...`; `-DCBK_TEST_LAPACK=openblas|netlib|mkl` -pins one, which matters on Debian, where `liblapack.so` is an -`update-alternatives` symlink and either stack can otherwise stand in for the -other. `-DCBK_BUILD_TESTS=OFF` drops the requirement with the tests. +CMake's own `find_package(LAPACK)` picks the implementation, steered by +`BLA_VENDOR` as usual (`-DBLA_VENDOR=OpenBLAS`, `Generic` for Netlib's +reference libraries, unset for the first found), and `cmake/FindLAPACKE.cmake` +adds `lapacke.h`, `cblas.h` and `liblapacke` where the distro ships LAPACKE as +a library of its own, reporting the result as `LAPACKE for the tests: ...`. +On Debian `liblapack.so` is an `update-alternatives` symlink, so `Generic` +means whichever implementation the alternative points at. +`-DCBK_BUILD_TESTS=OFF` drops the requirement with the tests. ```sh cmake -S . -B build -DCMAKE_BUILD_TYPE=Release @@ -30,9 +32,10 @@ ctest --test-dir build --output-on-failure On macOS, Accelerate provides LAPACK and CBLAS but no LAPACKE; [accelerate-lapacke](https://github.com/lepus2589/accelerate-lapacke) builds one on top of it (its configure needs a Fortran compiler, `brew install gcc`), -and `-DCBK_TEST_LAPACK=accelerate` with `-DCMAKE_PREFIX_PATH` at its install -(its own prefix; the workflow caches it) looks for that: the `lapacke` package plus -Accelerate's own CBLAS through ``. `.github/workflows/macos.yml` +and `-DBLA_VENDOR=Apple` with `-DCMAKE_PREFIX_PATH` at its install (its own +prefix; the workflow caches it) takes Accelerate from `find_package(LAPACK)`, +`lapacke.h` and `liblapacke` from that prefix, and Accelerate's own CBLAS +through ``. `.github/workflows/macos.yml` runs exactly this on the macOS runners, one job next to the Linux matrix. ## The Intel MKL extension (`-DCBK_WITH_MKL=ON`) @@ -149,7 +152,7 @@ CI builds against a fresh install for every MKL / static / shared combination | `CBK_WITH_MKL` | `OFF` | Build the MKL-style API, the MKL-backed tests, the example and the benchmarks; requires Intel MKL. | | `CBK_WITH_OPENMP` | `ON` | Thread the loop over groups with OpenMP; `OFF` gives single-threaded routines. See [threading.md](threading.md). | | `CBK_BUILD_TESTS` | `ON` | Build the tests and register them with CTest; needs a LAPACKE + CBLAS stack. | -| `CBK_TEST_LAPACK` | `auto` | The stack the tests validate against: `auto` (MKL's with the MKL extension, else the first of OpenBLAS, Netlib, any), `mkl`, `openblas`, `netlib`, `accelerate`. | +| `BLA_VENDOR` | -- | CMake's own: the LAPACK the tests validate against (`OpenBLAS`, `Generic` for Netlib, `Apple`); ignored with `CBK_WITH_MKL=ON`, where it is MKL's. | | `BUILD_SHARED_LIBS` | `OFF` | CMake's own switch: build `libcbk` shared instead of static. | | `CBK_INSTALL` | `ON` when top-level | Generate the install rules; off by default under a parent project's `add_subdirectory`. | | `MKLCompact_ROOT` | -- | A oneAPI MKL prefix, when `MKLROOT` is not set (MKL build only). | diff --git a/tests/test_gels_compact.cpp b/tests/test_gels_compact.cpp index bbe9933..ca23733 100644 --- a/tests/test_gels_compact.cpp +++ b/tests/test_gels_compact.cpp @@ -151,16 +151,15 @@ int run_case(char layout, char trans, int nm, int m, int n, int nrhs) } } - // vs LAPACK: two backward-stable implementations (LAPACK's blocked above - // its crossover) of a well-conditioned op(A)'s factorization and solve, - // relative to the operand norms (design 7.1) - const double tol_ref = 100.0 * eps * p; - const double tol_x = 100.0 * eps * p; - const double tol_prop = 100.0 * eps * p; // backward-stable quantities - // the Gram-formed minimum-norm solution carries cond(G) = cond(op(A))^2 - const double tol_prop2 = overdet ? tol_prop : 1e4 * eps * p; - const bool ok_x = e_x <= tol_x, ok_h = e_h <= tol_ref, ok_t = e_t <= tol_ref; - const bool ok_p = e_prop <= tol_prop, ok_p2 = e_prop2 <= tol_prop2; + // One gate for X, the factorization and the defining properties: two + // backward-stable implementations (LAPACK's blocked above its crossover) + // of a well-conditioned op(A)'s factorization and solve, relative to the + // operand norms (design 7.1). The Gram-formed minimum-norm solution + // carries cond(G) = cond(op(A))^2, so its gate gets that headroom. + const double tol = 100.0 * eps * p; + const double tol_prop2 = overdet ? tol : 1e4 * eps * p; + const bool ok_x = e_x <= tol, ok_h = e_h <= tol, ok_t = e_t <= tol; + const bool ok_p = e_prop <= tol, ok_p2 = e_prop2 <= tol_prop2; std::printf( "T=%-6s V=%-2d %s trans=%c nm=%-2d m=%-3d n=%-3d nrhs=%d %-6s | X:%.1e %s " diff --git a/tests/test_geqrf_compact.cpp b/tests/test_geqrf_compact.cpp index 3135e37..b4cde8e 100644 --- a/tests/test_geqrf_compact.cpp +++ b/tests/test_geqrf_compact.cpp @@ -122,20 +122,16 @@ template static int run_case(int nm, int m, int n) // columns included -- exactly as they do for dense LAPACK. enum Structure { DENSE, RANK_DEFICIENT, NEAR_COLLINEAR }; -// A random m x n matrix. For DENSE, a diagonal boost tames the conditioning -// and cond applies the competition's column scaling (columns *= -// logspace(0,-cond,n)). The structured variants make the last column a -// (near-)copy of the first, so the trailing reflector sees a (near-)zero -// sub-diagonal norm -- the branch the masked larfg must get right. +// A random m x n matrix: gen_boosted's, with the diagonal boost taming the +// conditioning for DENSE only, where cond applies the competition's column +// scaling (columns *= logspace(0,-cond,n)). The structured variants make the +// last column a (near-)copy of the first, so the trailing reflector sees a +// (near-)zero sub-diagonal norm -- the branch the masked larfg must get right. template static void gen_matrix(MatrixView A, double cond, Structure s) { const int m = A.rows, n = A.cols; - for (int j = 0; j < n; ++j) - for (int i = 0; i < m; ++i) - A(i, j) = frand(); + gen_boosted(A, s == DENSE ? T(2) : T(0)); if (s == DENSE) { - for (int d = 0; d < std::min(m, n); ++d) - A(d, d) += T(2); if (cond > 0.0) for (int j = 0; j < n; ++j) { double sc = std::pow(10.0, -cond * (n > 1 ? (double)j / (n - 1) : 0.0)); @@ -195,22 +191,23 @@ static int run_invariants(int nm, int m, int n, double cond, Structure structure QtQ(d, d) -= T(1); worst_orth = std::max(worst_orth, norm1(QtQ)); - // elementwise vs the blocked LAPACKE_?geqrf - copy_matrix(Am, Href); - ref_geqrf(Href, tauref.data()); - const double el = std::max(max_abs_diff(H[idx], Hrefs.data(), A.stride()), - max_abs_diff(tau[idx], tauref.data(), (size_t)k)); - worst_el = std::max(worst_el, el / std::max(norm1(Am), norm_floor)); + // elementwise vs the blocked LAPACKE_?geqrf, on the dense inputs only, + // where the reflectors are essentially unique: it catches a + // sign-convention regression the residual gates cannot see (observed + // ~n*eps). Rank-deficient and near-collinear inputs have no unique + // reflectors (el ~ 1e-2 there), so the comparison is skipped (design + // doc 7.1). + if (structure == DENSE) { + copy_matrix(Am, Href); + ref_geqrf(Href, tauref.data()); + const double el = std::max(max_abs_diff(H[idx], Hrefs.data(), A.stride()), + max_abs_diff(tau[idx], tauref.data(), (size_t)k)); + worst_el = std::max(worst_el, el / std::max(norm1(Am), norm_floor)); + } } - // The elementwise difference vs LAPACKE_?geqrf is gated on the dense - // inputs only, where the reflectors are essentially unique: it catches a - // sign-convention regression the residual gates cannot see (observed - // ~n*eps). Rank-deficient and near-collinear inputs stay ungated (their - // reflectors are not unique; el ~ 1e-2 is expected) and it is printed as a - // diagnostic there (design doc 7.1). const double rtol_res = 20.0 * n * eps, rtol_orth = 100.0 * n * eps; - const double rtol_el = (structure == DENSE) ? 100.0 * n * eps : HUGE_VAL; + const double rtol_el = 100.0 * n * eps; const bool ok = (info == 0) && (worst_res <= rtol_res) && (worst_orth <= rtol_orth) && (worst_el <= rtol_el); const char *sname = structure == DENSE ? "dense" @@ -223,6 +220,25 @@ static int run_invariants(int nm, int m, int n, double cond, Structure structure return !ok; } +// The invariants cases, in one precision and width. +template static int run_invariant_set() +{ + int fails = 0; + fails += run_invariants(8, 30, 30, 0.0); + fails += run_invariants(16, 60, 60, 0.0); + fails += run_invariants(11, 43, 43, 0.0); // padded partial group + fails += run_invariants(8, 64, 20, 0.0); // tall + fails += run_invariants(8, 20, 64, 0.0); // wide + fails += run_invariants(8, 30, 30, 4.0); // column-scaled (dynamic range) + fails += run_invariants(8, 128, 128, 0.0); + // conditioning-robustness stress: backward-stable gates must still hold + fails += run_invariants(8, 40, 40, 0.0, RANK_DEFICIENT); + fails += run_invariants(8, 40, 40, 0.0, NEAR_COLLINEAR); + fails += + run_invariants(11, 60, 24, 0.0, RANK_DEFICIENT); // wide-ish, padded group + return fails; +} + // ------------------------ underflow scope (design 6.6) --------------- // A column whose sub-diagonal entries are so small that their squares @@ -345,23 +361,8 @@ int main() // the dense-LAPACK contract (design 7.1): residual and orthogonality of // the materialized Q, over conditioning and structure, in both precisions - for (int prec = 0; prec < 2; ++prec) { - auto inv = [&](int nm, int m, int n, double cond, Structure st = DENSE) { - return prec ? run_invariants(nm, m, n, cond, st) - : run_invariants(nm, m, n, cond, st); - }; - fails += inv(8, 30, 30, 0.0); - fails += inv(16, 60, 60, 0.0); - fails += inv(11, 43, 43, 0.0); // padded partial group - fails += inv(8, 64, 20, 0.0); // tall - fails += inv(8, 20, 64, 0.0); // wide - fails += inv(8, 30, 30, 4.0); // column-scaled (dynamic range) - fails += inv(8, 128, 128, 0.0); - // conditioning-robustness stress: backward-stable gates must still hold - fails += inv(8, 40, 40, 0.0, RANK_DEFICIENT); - fails += inv(8, 40, 40, 0.0, NEAR_COLLINEAR); - fails += inv(11, 60, 24, 0.0, RANK_DEFICIENT); // wide-ish, padded group - } + fails += run_invariant_set(); + fails += run_invariant_set(); return finish(fails); } diff --git a/tests/test_lapack_util.hpp b/tests/test_lapack_util.hpp index 5f9fa2c..4d7f5b5 100644 --- a/tests/test_lapack_util.hpp +++ b/tests/test_lapack_util.hpp @@ -36,10 +36,8 @@ // LAPACK 3.9.1 interface), CBLAS from Accelerate's vecLib sub-framework // directly -- not through the umbrella, whose LAPACK // prototypes would collide with the ones lapacke.h pulls in. The framework -// search path for is on the LAPACKE::LAPACKE target. -#ifndef ACCELERATE_NEW_LAPACK -#define ACCELERATE_NEW_LAPACK -#endif +// search path for and ACCELERATE_NEW_LAPACK (the LAPACK 3.9.1 +// interface the shim is built for) are on the LAPACKE::LAPACKE target. #include #include #else @@ -50,14 +48,13 @@ #include #include #include -#include #include namespace cbk::test { // ----------------------- LAPACKE / CBLAS, by scalar type ------------- -// lapack::geqr2 / geqrf / gelqf / orgqr / ormqr / geqp3 / potrf / gels / -// larfg forward to LAPACKE_?*, gemm / trsm / trmm to cblas_?*. One macro +// lapack::geqr2 / geqrf / gelqf / orgqr / ormqr / geqp3 / potrf / gels +// forward to LAPACKE_?*, gemm / trsm / trmm to cblas_?*. One macro // generates both specializations so the two cannot drift apart. template struct lapack; @@ -66,8 +63,6 @@ template struct lapack; // NOLINTBEGIN(bugprone-macro-parentheses): T is a type name, p a token to paste #define CBK_TEST_LAPACK_DISPATCH(T, p) \ template <> struct lapack { \ - static lapack_int larfg(lapack_int n, T *alpha, T *x, lapack_int incx, T *tau) \ - { return LAPACKE_##p##larfg(n, alpha, x, incx, tau); } \ static lapack_int geqr2(int layout, lapack_int m, lapack_int n, T *a, lapack_int lda, \ T *tau) \ { return LAPACKE_##p##geqr2(layout, m, n, a, lda, tau); } \ @@ -131,7 +126,7 @@ template bool is_rowmajor(const Mv &M) template lapack_int lapack_ld(const Mv &M) { - return std::max(is_rowmajor(M) ? M.si : M.sj, 1); + return std::max(M.ld(), 1); } template int lapack_layout(const Mv &M) @@ -170,31 +165,25 @@ inline CBLAS_DIAG cblas_diag(char d) return (d == 'U' || d == 'u') ? CblasUnit : CblasNonUnit; } -// An operand of a call whose layout is fixed by another operand: the view -// itself when its layout already matches, else a contiguous staging copy in -// the call's layout -- copied back on destruction when the operand is -// writable (T non-const). Only the LAPACKE routines without a transpose flag -// for the operand need it (ormqr's reflectors, gels's right-hand side); the -// BLAS-3 calls read a mismatched operand as its transpose instead. +// The writable operand of a call whose layout another operand fixes: the +// view itself when its layout already matches, else a contiguous staging +// copy in the call's layout, copied back on destruction. Only the LAPACKE +// routines without a transpose flag need it (ormqr's and gels's right-hand +// sides, staged into the layout of the factor); the BLAS-3 calls read a +// mismatched operand as its transpose instead. template class Staged { public: - Staged(MatrixView M, bool rowmajor) : orig_(M), copied_(is_rowmajor(M) != rowmajor) + Staged(MatrixView M, bool rowmajor) : orig_(M), use_(M) { - if (!copied_) { - use_ = M; - return; - } + if (is_rowmajor(M) == rowmajor) return; buf_.resize((std::size_t)M.rows * M.cols); - const auto w = + use_ = mat_view(buf_.data(), M.rows, M.cols, rowmajor ? M.cols : M.rows, rowmajor); - copy_matrix(M, w); - use_ = MatrixView{w.data, w.si, w.sj, w.rows, w.cols}; + copy_matrix(M, use_); } ~Staged() { - if constexpr (!std::is_const_v) { - if (copied_) copy_matrix(use_, orig_); - } + if (!buf_.empty()) copy_matrix(use_, orig_); } Staged(const Staged &) = delete; Staged &operator=(const Staged &) = delete; @@ -203,19 +192,10 @@ template class Staged { lapack_int ld() const { return lapack_ld(use_); } private: - MatrixView orig_; - bool copied_; - std::vector> buf_; - MatrixView use_{}; + MatrixView orig_, use_; + std::vector buf_; }; -// A read-only view of the same storage, whatever the constness of the view -// handed in (distinct MatrixView instantiations do not convert). -template ConstMatrixView> as_const_view(const Mv &M) -{ - return {M.data, M.si, M.sj, M.rows, M.cols}; -} - // ----------------------- reference procedures ------------------------ // The names the suites call, each forwarded to the library. LAPACK's blocked // drivers (?geqrf, ?ormqr, ?orgqr, ?potrf, ?gels) run their unblocked @@ -224,14 +204,8 @@ template ConstMatrixView> as_const_view(const Mv &M) // comparison against them is a comparison of two backward-stable // implementations of the same factorization, not of the same operation // sequence: the suites gate such comparisons relative to the operand norms, -// at a multiple of n * eps. Only ?geqr2 and ?larfg are unblocked by -// definition, and only those two are compared elementwise at ~eps. - -// ?larfg: reflector from (alpha, x[0..m-2]); alpha := beta on exit. -template void ref_larfg(int m, T *alpha, T *x, T *tau) -{ - lapack::larfg(m, alpha, x, 1, tau); -} +// at a multiple of n * eps. Only ?geqr2 is unblocked by definition, and only +// it is compared elementwise at ~eps. // ?geqr2: unblocked Householder QR, (H, tau) in the LAPACK convention. template void ref_geqr2(MatrixView A, T *tau) @@ -251,10 +225,9 @@ template void ref_orm2r(char trans, int k, Av A, const T *tau, MatrixView B) { assert(A.rows == B.rows && k <= std::min(A.rows, A.cols)); - const bool row = is_rowmajor(B); - const Staged> As(as_const_view(A), row); - lapack::ormqr(lapack_layout(B), 'L', trans, B.rows, B.cols, k, As.data(), As.ld(), - tau, B.data, lapack_ld(B)); + const Staged Bs(B, is_rowmajor(A)); + lapack::ormqr(lapack_layout(A), 'L', trans, B.rows, B.cols, k, A.data, + lapack_ld(A), tau, Bs.data(), Bs.ld()); } // ?orgqr: generate the first n columns of Q = H(0)..H(k-1) in place over the @@ -331,9 +304,8 @@ template void ref_trsm_upper(Rv R, MatrixView B) // On exit B holds X (least squares: rows n..m-1 keep the residual), A the // factorization of the tall orientation -- the QR of A (m >= n) or its LQ in // ?gelqf storage (m < n) -- and tau its min(m,n) reflector scalars. LAPACK's -// ?gels leaves that factorization in A but keeps tau in its workspace, so -// (A, tau) are taken from ?geqrf / ?gelqf on a copy of the input, which is -// the same call ?gels makes. +// ?gels keeps tau in its workspace, so it runs on a copy of A for X, and +// (A, tau) come from ?geqrf / ?gelqf on A itself, the same call ?gels makes. template void ref_gels(char trans, MatrixView A, MatrixView B, T *tau) { const int m = A.rows, n = A.cols, nrhs = B.cols; @@ -343,17 +315,14 @@ template void ref_gels(char trans, MatrixView A, MatrixView B, T std::vector A0s((std::size_t)m * n); const auto A0 = mat_view(A0s.data(), m, n, row ? n : m, row); copy_matrix(A, A0); + const Staged Bs(B, row); + lapack::gels(lapack_layout(A0), trans, m, n, nrhs, A0.data, lapack_ld(A0), + Bs.data(), Bs.ld()); - { - const Staged Bs(B, row); - lapack::gels(lapack_layout(A), trans, m, n, nrhs, A.data, lapack_ld(A), - Bs.data(), Bs.ld()); - } if (m >= n) - lapack::geqrf(lapack_layout(A0), m, n, A0.data, lapack_ld(A0), tau); + lapack::geqrf(lapack_layout(A), m, n, A.data, lapack_ld(A), tau); else - lapack::gelqf(lapack_layout(A0), m, n, A0.data, lapack_ld(A0), tau); - copy_matrix(A0, A); + lapack::gelqf(lapack_layout(A), m, n, A.data, lapack_ld(A), tau); } // C (m x n) := A (m x k) * B (k x n) by ?gemm, for forming right-hand sides @@ -389,10 +358,9 @@ template void tri_apply(char side, char uplo, char transa, char diag, Av A, Xv X, Rv R) { using T = elem_t; - const bool left = (side == 'L' || side == 'l'); - const int s = left ? R.rows : R.cols; - assert(A.rows == s && A.cols == s && X.rows == R.rows && X.cols == R.cols); - (void)s; + assert(A.rows == A.cols && + A.rows == ((side == 'L' || side == 'l') ? R.rows : R.cols) && + X.rows == R.rows && X.cols == R.cols); copy_matrix(X, R); const bool row = is_rowmajor(R); CBLAS_TRANSPOSE ta = cblas_trans(transa); diff --git a/tests/test_potrf_compact.cpp b/tests/test_potrf_compact.cpp index c6b009d..fff1823 100644 --- a/tests/test_potrf_compact.cpp +++ b/tests/test_potrf_compact.cpp @@ -90,11 +90,11 @@ static int run_case(int nm, int n, char uplo, char layout, double cond = 0.0) std::max(norm1(Ain), norm_floor)); } - const double scale = std::max(1, n); - const double tol_fac = 20.0 * eps * scale; // design 7.1: unique factor, n eps - const double tol_rec = 20.0 * eps * scale; // ||L L^T - A|| / ||A||: backward stable - bool ok_f = e_fac <= tol_fac; - bool ok_r = e_rec <= tol_rec; + // design 7.1: the unique factor vs LAPACK's, and ||L L^T - A|| / ||A||, + // both backward-stable quantities at n eps + const double tol = 20.0 * eps * std::max(1, n); + bool ok_f = e_fac <= tol; + bool ok_r = e_rec <= tol; bool ok_u = e_untouched == 0.0; // must be bit-for-bit unchanged bool ok_i = (info == 0); @@ -233,12 +233,11 @@ template static int run_nonspd(int n, char uplo, char layout) if (!named) continue; const T x = Aout(idx, i, j); if (idx == badlane) { - if (!std::isfinite((double)x)) bad_poisoned = true; - } - else { - if (!std::isfinite((double)x)) spd_finite = false; - el = std::max(el, (double)std::abs(x - Aref(idx, i, j))); + bad_poisoned = bad_poisoned || !std::isfinite((double)x); + continue; } + spd_finite = spd_finite && std::isfinite((double)x); + el = std::max(el, (double)std::abs(x - Aref(idx, i, j))); } if (idx != badlane) // vs LAPACKE_?potrf, relative to its factor's norm e_spd = std::max(e_spd, el / std::max(norm1(Aref.view(idx)), norm_floor)); diff --git a/tests/test_sytrfnp_compact.cpp b/tests/test_sytrfnp_compact.cpp index 538a586..d05bdd4 100644 --- a/tests/test_sytrfnp_compact.cpp +++ b/tests/test_sytrfnp_compact.cpp @@ -72,42 +72,30 @@ template static void ref_sytf2np(char uplo, MatrixView As) // ------------- the reference itself, vs the library ------------------ // ref_sytf2np is the only hand-rolled reference left in the suites, so it is // checked against LAPACKE/CBLAS before anything is checked against it: its -// (L, D) must reconstruct A. Lower: W = L D (columns of the unit-lower L -// scaled by D), then W := W L^T by ?trmm (side 'R', unit diagonal). Upper: -// W = D U (rows scaled), then W := U^T W (side 'L'). Gated at +// (L, D) must reconstruct A. W = L D (the columns of the unit-lower L scaled +// by D), then W := W L^T by ?trmm (side 'R', unit diagonal); the upper +// factor U^T D U is the same product read through the transposed view, as in +// ref_sytf2np itself, and A is symmetric. Gated at // ||L D L^T - A||_1 / ||A||_1 <= 20 n eps, the design's 7.1 residual. template static int test_reference(int n, char uplo) { const T eps = std::numeric_limits::epsilon(); const bool upper = (uplo == 'U' || uplo == 'u'); - const int trials = 4; - double worst = 0; std::vector As((size_t)n * n), Fs((size_t)n * n), Ws((size_t)n * n), Rs((size_t)n * n); const auto A = mat_view(As.data(), n, n), F = mat_view(Fs.data(), n, n), W = mat_view(Ws.data(), n, n), R = mat_view(Rs.data(), n, n); - for (int t = 0; t < trials; ++t) { - gen_sym_ldlt(A); - copy_matrix(A, F); - ref_sytf2np(uplo, F); - // W: the unit factor with D folded in, the other triangle zeroed - for (int j = 0; j < n; ++j) - for (int i = 0; i < n; ++i) { - if (i == j) - W(i, j) = F(i, i); - else if (upper) - W(i, j) = (i < j) ? F(i, i) * F(i, j) : T(0); // D U - else - W(i, j) = (i > j) ? F(i, j) * F(j, j) : T(0); // L D - } - if (upper) - tri_apply('L', 'U', 'T', 'U', F, W, R); // R = U^T (D U) - else - tri_apply('R', 'L', 'T', 'U', F, W, R); // R = (L D) L^T - worst = std::max(worst, max_abs_diff(Rs.data(), As.data(), As.size()) / - std::max(norm1(A), norm_floor)); - } + gen_sym_ldlt(A); + copy_matrix(A, F); + ref_sytf2np(uplo, F); + const auto L = upper ? F.transposed() : F; // the unit-lower factor, D on its diagonal + for (int j = 0; j < n; ++j) + for (int i = 0; i < n; ++i) + W(i, j) = (i == j) ? L(j, j) : (i > j) ? L(i, j) * L(j, j) : T(0); // L D + tri_apply('R', 'L', 'T', 'U', L, W, R); // (L D) L^T + const double worst = + max_abs_diff(Rs.data(), As.data(), As.size()) / std::max(norm1(A), norm_floor); const double tol = 20.0 * eps * std::max(1, n); const bool ok = worst <= tol; std::printf( From 8044adbd7979728a92660695886cb39308d7203d Mon Sep 17 00:00:00 2001 From: Claude Date: Mon, 14 Sep 2026 23:08:11 +0000 Subject: [PATCH 6/8] macOS workflow: cache prefix without relative path segments actions/cache rejects a path containing '..' ("Relative pathing '.' and '..' is not allowed") and saved nothing, so every run rebuilt accelerate-lapacke. The install prefix is now under the runner's temp directory, an absolute path, set at the job level where that context is available. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_01GdNaquWMN3b39gvTsXZyjn --- .github/workflows/macos.yml | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/.github/workflows/macos.yml b/.github/workflows/macos.yml index bebae62..2edb536 100644 --- a/.github/workflows/macos.yml +++ b/.github/workflows/macos.yml @@ -26,14 +26,16 @@ concurrency: group: ${{ github.workflow }}-${{ github.ref }} cancel-in-progress: true -env: - ACCELERATE_LAPACKE_REF: v2.0.0 - ACCELERATE_LAPACKE_PREFIX: ${{ github.workspace }}/../accelerate-lapacke - jobs: accelerate: name: AppleClang / Accelerate runs-on: macos-latest + # The install prefix must be an absolute path without `.` or `..` + # segments: actions/cache rejects those and then saves nothing (the + # runner's temp directory is only known at the job level, hence env here). + env: + ACCELERATE_LAPACKE_REF: v2.0.0 + ACCELERATE_LAPACKE_PREFIX: ${{ runner.temp }}/accelerate-lapacke-install steps: - uses: actions/checkout@v4 From a8f0380ccfb4473c87ef14433533216da758043f Mon Sep 17 00:00:00 2001 From: Claude Date: Mon, 14 Sep 2026 23:13:37 +0000 Subject: [PATCH 7/8] macOS workflow: the runner context is step-level only A job-level env cannot read runner.temp, and the workflow file failed to parse (a run with no jobs). The cache path names the runner's temp directory in the step itself and the shell steps use $RUNNER_TEMP. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_01GdNaquWMN3b39gvTsXZyjn --- .github/workflows/macos.yml | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/.github/workflows/macos.yml b/.github/workflows/macos.yml index 2edb536..829cd07 100644 --- a/.github/workflows/macos.yml +++ b/.github/workflows/macos.yml @@ -30,12 +30,8 @@ jobs: accelerate: name: AppleClang / Accelerate runs-on: macos-latest - # The install prefix must be an absolute path without `.` or `..` - # segments: actions/cache rejects those and then saves nothing (the - # runner's temp directory is only known at the job level, hence env here). env: ACCELERATE_LAPACKE_REF: v2.0.0 - ACCELERATE_LAPACKE_PREFIX: ${{ runner.temp }}/accelerate-lapacke-install steps: - uses: actions/checkout@v4 @@ -47,11 +43,15 @@ jobs: echo "version=$sdk" >> "$GITHUB_OUTPUT" echo "macOS $(sw_vers -productVersion), SDK $sdk" + # The install prefix, $RUNNER_TEMP/accelerate-lapacke-install below, is + # an absolute path without `.` or `..` segments: actions/cache rejects + # those and then saves nothing. (The runner context that names it is + # only available inside steps, so it is spelled out rather than an env.) - name: Cache accelerate-lapacke id: lapacke-cache uses: actions/cache@v4 with: - path: ${{ env.ACCELERATE_LAPACKE_PREFIX }} + path: ${{ runner.temp }}/accelerate-lapacke-install key: accelerate-lapacke-${{ env.ACCELERATE_LAPACKE_REF }}-${{ runner.os }}-${{ runner.arch }}-sdk${{ steps.sdk.outputs.version }} # accelerate-lapacke's project enables Fortran (Netlib's LAPACK build @@ -70,16 +70,16 @@ jobs: cd "$RUNNER_TEMP/accelerate-lapacke" cmake --preset accelerate-lapacke32 \ -DCMAKE_OSX_SYSROOT="${{ steps.sdk.outputs.path }}" \ - -DCMAKE_INSTALL_PREFIX="$ACCELERATE_LAPACKE_PREFIX" + -DCMAKE_INSTALL_PREFIX="$RUNNER_TEMP/accelerate-lapacke-install" cmake --build --preset accelerate-lapacke32-install - ls "$ACCELERATE_LAPACKE_PREFIX/include/lapacke.h" "$ACCELERATE_LAPACKE_PREFIX/lib/" + ls "$RUNNER_TEMP/accelerate-lapacke-install/include/lapacke.h" "$RUNNER_TEMP/accelerate-lapacke-install/lib/" - name: Configure run: > cmake -S . -B build -DCMAKE_BUILD_TYPE=Release -DBLA_VENDOR=Apple - -DCMAKE_PREFIX_PATH="$ACCELERATE_LAPACKE_PREFIX" + -DCMAKE_PREFIX_PATH="$RUNNER_TEMP/accelerate-lapacke-install" - name: Build run: cmake --build build -j From fa05c1c1bdf22ee87538bf26f647c8e10813ca53 Mon Sep 17 00:00:00 2001 From: Claude Date: Mon, 14 Sep 2026 23:13:37 +0000 Subject: [PATCH 8/8] Tests: name the references after the routines they call ref_ormqr, ref_orgqr and ref_potrf forward to LAPACKE_?ormqr, ?orgqr and ?potrf, the blocked drivers, so the unblocked names orm2r, org2r and potf2 no longer described them. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_01GdNaquWMN3b39gvTsXZyjn --- .claude/CLAUDE.md | 2 +- docs/cbk_dpotrf_compact_design.md | 2 +- docs/layout.md | 2 +- tests/test_gels_compact.cpp | 2 +- tests/test_geqrf_compact.cpp | 6 +++--- tests/test_lapack_util.hpp | 10 +++++----- tests/test_orgqr_compact.cpp | 4 ++-- tests/test_ormqr_compact.cpp | 2 +- tests/test_potrf_compact.cpp | 6 +++--- 9 files changed, 18 insertions(+), 18 deletions(-) diff --git a/.claude/CLAUDE.md b/.claude/CLAUDE.md index ff4a326..b62564f 100644 --- a/.claude/CLAUDE.md +++ b/.claude/CLAUDE.md @@ -219,7 +219,7 @@ workspace contract, or the benchmarks' threading. A routine that takes views (or batches) `assert`s what its contract assumes: dimension compatibility across operands (`matmul`, `tri_apply`, `solve_errors`, the `ref_*` procedures), squareness where required - (`gen_spd`, `gen_tri`, `ref_potf2`, ...), and index ranges + (`gen_spd`, `gen_tri`, `ref_potrf`, ...), and index ranges (`MatrixView::operator()`, `MatrixBatch::operator[]`). Costs nothing in Release, and the Debug suite run is what exercises it. `BatchView` is the exception by construction -- it carries strides, not extents, so a kernel diff --git a/docs/cbk_dpotrf_compact_design.md b/docs/cbk_dpotrf_compact_design.md index 84ad29b..580b499 100644 --- a/docs/cbk_dpotrf_compact_design.md +++ b/docs/cbk_dpotrf_compact_design.md @@ -441,7 +441,7 @@ with the rest of the compact toolkit. ### 7.4 Portable self-test `tests/test_potrf_compact.cpp` validates the templated kernels directly -against `LAPACKE_?potrf` (through `ref_potf2`; the unique SPD factor, +against `LAPACKE_?potrf` (through `ref_potrf`; the unique SPD factor, relative to its norm at `20 n eps`) across `(T, V)` combinations, both `uplo`, both layouts, the `cond` knob, and partial (padded) final packs, with the 7.1 invariants (reconstruction by `?trmm`, the untouched triangle) and diff --git a/docs/layout.md b/docs/layout.md index 92275b5..892add9 100644 --- a/docs/layout.md +++ b/docs/layout.md @@ -48,7 +48,7 @@ scalar, `0` on success). | `src/cbk.cpp` | The portable C API: argument validation and `V` dispatch for all twenty-two entry points. | | `src/cbk_compat.cpp` | The MKL-style API: MKL enum / `MKL_COMPACT_PACK` unwrapping for all twenty-two entry points. | | `tests/test_compact_util.hpp` | Shared test helpers: RNG, error metrics, input generation, Compact pack/unpack, and `compact` (the portable C API dispatched on the scalar type); header-only, library-free. | -| `tests/test_lapack_util.hpp` | The dense references on a real LAPACKE + CBLAS: `lapack` (the LAPACKE/CBLAS dispatch), the `ref_*` procedures over views (`ref_geqr2`, `ref_orm2r`, `ref_potf2`, `ref_trsm`, `ref_gels`, ...), `matmul`, `tri_apply`, `solve_errors`. The stack comes from `cmake/FindLAPACKE.cmake`. | +| `tests/test_lapack_util.hpp` | The dense references on a real LAPACKE + CBLAS: `lapack` (the LAPACKE/CBLAS dispatch), the `ref_*` procedures over views (`ref_geqr2`, `ref_ormqr`, `ref_potrf`, `ref_trsm`, `ref_gels`, ...), `matmul`, `tri_apply`, `solve_errors`. The stack comes from `cmake/FindLAPACKE.cmake`. | | `tests/test_mkl_util.hpp` | Scalar-type dispatch for the MKL-backed suites: `compat` (routines under test), `mkl` (MKL's compact API and kernels). | | `tests/test_*_compact.cpp` | The portable C API's suites, on any LAPACKE stack: each kernel vs its LAPACK counterpart (elementwise where that is the same unblocked algorithm, residual and invariant gates elsewhere), plus C API validation, in FP64 and FP32. | | `tests/test_*_mkl.cpp` | MKL validation, templated on the scalar type and run in FP64 and FP32: the MKL-style API through MKL's pack/unpack, cross-checks vs MKL's compact kernels, end-to-end solves. | diff --git a/tests/test_gels_compact.cpp b/tests/test_gels_compact.cpp index ca23733..bb8c90b 100644 --- a/tests/test_gels_compact.cpp +++ b/tests/test_gels_compact.cpp @@ -141,7 +141,7 @@ int run_case(char layout, char trans, int nm, int m, int n, int nrhs) matmul(Aop, Aop.transposed(), G); copy_matrix(Bin, Z); ref_geqr2(G, tg.data()); - ref_orm2r('T', rows_op, G, tg.data(), Z); + ref_ormqr('T', rows_op, G, tg.data(), Z); ref_trsm_upper(G, Z); matmul(Aop.transposed(), Z, Xmn); for (int j = 0; j < nrhs; ++j) diff --git a/tests/test_geqrf_compact.cpp b/tests/test_geqrf_compact.cpp index b4cde8e..07b15c6 100644 --- a/tests/test_geqrf_compact.cpp +++ b/tests/test_geqrf_compact.cpp @@ -73,7 +73,7 @@ template static int run_case(int nm, int m, int n) for (int j = 0; j < n; ++j) for (int i = 0; i <= std::min(j, k - 1); ++i) Rec(i, j) = Aout(idx, i, j); - ref_orm2r('N', k, Aout.view(idx), tau_out[idx], Rec); + ref_ormqr('N', k, Aout.view(idx), tau_out[idx], Rec); e_rec = std::max(e_rec, max_abs_diff(Recs.data(), A[idx], (size_t)m * n)); } @@ -174,7 +174,7 @@ static int run_invariants(int nm, int m, int n, double cond, Structure structure // Q = the first k columns of H(0)..H(k-1), from the reflectors in the // first k columns of H copy_matrix(leading(Hm, m, k), Q); - ref_org2r(k, Q, tau[idx]); + ref_orgqr(k, Q, tau[idx]); // residual R - Q^T A, R = triu(H) (k x n) matmul(Q.transposed(), Am, QtA); @@ -292,7 +292,7 @@ template static int test_underflow() for (int j = 0; j < n; ++j) for (int i = 0; i <= std::min(j, n - 1); ++i) Rec(i, j) = Aout(idx, i, j); - ref_orm2r('N', n, Aout.view(idx), tau[idx], Rec); + ref_ormqr('N', n, Aout.view(idx), tau[idx], Rec); e_rec = std::max(e_rec, max_abs_diff(Recs.data(), A[idx], (size_t)m * n)); } const double tol_rec = 200.0 * eps * m; diff --git a/tests/test_lapack_util.hpp b/tests/test_lapack_util.hpp index 4d7f5b5..029c27c 100644 --- a/tests/test_lapack_util.hpp +++ b/tests/test_lapack_util.hpp @@ -3,8 +3,8 @@ // The dense references of every test suite, on a real LAPACKE + CBLAS stack // (issue #27): lapack, a struct of static forwarders to the LAPACKE_?* and // cblas_?* entry points specialized for double (d) and float (s), and over it -// the reference procedures the suites call by name -- ref_geqr2, ref_orm2r, -// ref_potf2, ref_trsm, ref_gels, ..., matmul, tri_apply, solve_errors -- each +// the reference procedures the suites call by name -- ref_geqr2, ref_ormqr, +// ref_potrf, ref_trsm, ref_gels, ..., matmul, tri_apply, solve_errors -- each // taking the project's MatrixView operands and forwarding their `.data` and // leading dimensions to the library. The only reference the suites still // hand-roll is the unpivoted LDL^T of test_sytrfnp_compact.cpp (LAPACK has no @@ -222,7 +222,7 @@ template void ref_geqrf(MatrixView A, T *tau) // ?ormqr, side='L': B := Q^T B (trans 'T') or Q B ('N') from (H, tau), with // k reflectors read from A (m x k or wider) and applied to B (m x nrhs). template -void ref_orm2r(char trans, int k, Av A, const T *tau, MatrixView B) +void ref_ormqr(char trans, int k, Av A, const T *tau, MatrixView B) { assert(A.rows == B.rows && k <= std::min(A.rows, A.cols)); const Staged Bs(B, is_rowmajor(A)); @@ -232,7 +232,7 @@ void ref_orm2r(char trans, int k, Av A, const T *tau, MatrixView B) // ?orgqr: generate the first n columns of Q = H(0)..H(k-1) in place over the // reflectors in columns 0..k-1 of A (m x n, m >= n >= k). -template void ref_org2r(int k, MatrixView A, const T *tau) +template void ref_orgqr(int k, MatrixView A, const T *tau) { assert(k <= A.cols && A.cols <= A.rows); lapack::orgqr(lapack_layout(A), A.rows, A.cols, k, A.data, lapack_ld(A), tau); @@ -254,7 +254,7 @@ template void ref_geqp3(MatrixView A, lapack_int *jpvt, T *tau) // ?potrf on a dense n x n matrix, in place. Lower: A = L L^T, factor in the // lower triangle. Upper: A = U^T U, factor in the upper triangle. Only the // named triangle is read or written. Returns LAPACK's info (> 0: not SPD). -template lapack_int ref_potf2(char uplo, MatrixView A) +template lapack_int ref_potrf(char uplo, MatrixView A) { assert(A.rows == A.cols); return lapack::potrf(lapack_layout(A), uplo, A.rows, A.data, lapack_ld(A)); diff --git a/tests/test_orgqr_compact.cpp b/tests/test_orgqr_compact.cpp index e0eb9ae..db40bf6 100644 --- a/tests/test_orgqr_compact.cpp +++ b/tests/test_orgqr_compact.cpp @@ -31,7 +31,7 @@ #include #include -#include "test_lapack_util.hpp" // compact, ref_geqr2, ref_org2r, MatrixBatch, pack/unpack +#include "test_lapack_util.hpp" // compact, ref_geqr2, ref_orgqr, MatrixBatch, pack/unpack using namespace cbk::test; @@ -63,7 +63,7 @@ static int run_case(int nm, int m, int n, int k, bool rowmajor = false) H(i, j) = (j < k) ? Afac(kk, i, j) : frand(); /* garbage cols */ std::copy(Hn[kk], Hn[kk] + Hn.stride(), Qref[kk]); - ref_org2r(k, Qref.view(kk), tau[kk]); + ref_orgqr(k, Qref.view(kk), tau[kk]); } std::vector ap = pack_compact(Hn, ld, V, rowmajor); diff --git a/tests/test_ormqr_compact.cpp b/tests/test_ormqr_compact.cpp index 38d2d7b..7324067 100644 --- a/tests/test_ormqr_compact.cpp +++ b/tests/test_ormqr_compact.cpp @@ -53,7 +53,7 @@ template static int run_case(int nm, int m, int nrhs) std::copy(A[kk], A[kk] + (size_t)m * m, Afac[kk]); ref_geqr2(Afac.view(kk), tau[kk]); std::copy(B[kk], B[kk] + (size_t)m * nrhs, Bref[kk]); - ref_orm2r('T', k, Afac.view(kk), tau[kk], Bref.view(kk)); + ref_ormqr('T', k, Afac.view(kk), tau[kk], Bref.view(kk)); } std::vector ap = pack_compact(Afac, m, V); diff --git a/tests/test_potrf_compact.cpp b/tests/test_potrf_compact.cpp index fff1823..d505fa7 100644 --- a/tests/test_potrf_compact.cpp +++ b/tests/test_potrf_compact.cpp @@ -28,7 +28,7 @@ #include #include -#include "test_lapack_util.hpp" // compact, ref_potf2, tri_apply, gen_spd, pack/unpack +#include "test_lapack_util.hpp" // compact, ref_potrf, tri_apply, gen_spd, pack/unpack using namespace cbk::test; @@ -47,7 +47,7 @@ static int run_case(int nm, int n, char uplo, char layout, double cond = 0.0) for (int idx = 0; idx < nm; ++idx) { gen_spd(A.view(idx), cond); std::copy(A[idx], A[idx] + (size_t)n * n, Aref[idx]); - ref_potf2(uplo, Aref.view(idx)); + ref_potrf(uplo, Aref.view(idx)); } // pack the full symmetric A, factor with the routine under test, unpack @@ -213,7 +213,7 @@ template static int run_nonspd(int n, char uplo, char layout) else { gen_spd(A.view(idx)); std::copy(A[idx], A[idx] + (size_t)n * n, Aref[idx]); - ref_potf2(uplo, Aref.view(idx)); + ref_potrf(uplo, Aref.view(idx)); } }