diff --git a/.claude/CLAUDE.md b/.claude/CLAUDE.md index 5eefa2c..b62564f 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; `-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 @@ -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, and Accelerate on macOS through `.github/workflows/macos.yml`): ```sh - CXX=clang++ cmake -S . -B build-clang && 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 @@ -201,10 +211,15 @@ 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 - (`gen_spd`, `gen_tri`, `ref_potf2`, ...), and index ranges + `solve_errors`, the `ref_*` procedures), squareness where required + (`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/.claude/hooks/session-start.sh b/.claude/hooks/session-start.sh index afe62d5..688ed82 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 -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 @@ -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 -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 64088e0..f351021 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -12,41 +12,47 @@ concurrency: jobs: build-and-test: - name: ${{ matrix.cxx }} / MKL ${{ matrix.mkl }}${{ 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] - mkl: ["ON", "OFF"] + # 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] 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). - - { cc: gcc, cxx: g++, mkl: "ON", interface: ilp64 } + - { 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. The portable - # (MKL OFF) build needs nothing beyond the C++17 toolchain the image ships. - - 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 + 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_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 new file mode 100644 index 0000000..829cd07 --- /dev/null +++ b/.github/workflows/macos.yml @@ -0,0 +1,88 @@ +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 interface; its configure needs a Fortran +# compiler, which Homebrew's gcc supplies) into its own prefix, then +# 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 +# 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] + tags: ['v*'] + pull_request: + workflow_dispatch: + +concurrency: + group: ${{ github.workflow }}-${{ github.ref }} + cancel-in-progress: true + +jobs: + accelerate: + name: AppleClang / Accelerate + runs-on: macos-latest + env: + ACCELERATE_LAPACKE_REF: v2.0.0 + steps: + - uses: actions/checkout@v4 + + - name: Identify the SDK + id: sdk + run: | + sdk=$(xcrun --show-sdk-version) + echo "path=$(xcrun --show-sdk-path)" >> "$GITHUB_OUTPUT" + 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: ${{ 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 + # system does, even though only its C LAPACKE sources are compiled and + # 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: | + 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="${{ steps.sdk.outputs.path }}" \ + -DCMAKE_INSTALL_PREFIX="$RUNNER_TEMP/accelerate-lapacke-install" + cmake --build --preset accelerate-lapacke32-install + 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="$RUNNER_TEMP/accelerate-lapacke-install" + + - 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..651a2f8 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,6 +5,20 @@ 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` 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 + now needs such a stack; the library does not. A macOS workflow tests + against Accelerate through accelerate-lapacke. + ## [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..1d1a5ef 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -27,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. @@ -42,7 +48,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 +115,30 @@ 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): 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 + # 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..24089dd --- /dev/null +++ b/cmake/FindLAPACKE.cmake @@ -0,0 +1,124 @@ +# FindLAPACKE.cmake +# +# The LAPACKE + CBLAS stack the test suites validate against, as one target: +# +# 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 / ). +# +# 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: +# +# 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. +# +# 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 + +include(CMakePushCheckState) +include(CheckCXXSourceCompiles) + +set(_lapacke_includes "") +set(_lapacke_libs "") +set(_lapacke_defs "") +set(_lapacke_options "") + +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() + 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() + 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() + +# 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_LINKS + REASON_FAILURE_MESSAGE + "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_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_libs})") +endif() diff --git a/docs/PLANS.md b/docs/PLANS.md index ee2c9c2..bea94e0 100644 --- a/docs/PLANS.md +++ b/docs/PLANS.md @@ -93,11 +93,21 @@ listed in `examples.md`. ## Project-wide -- **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. +- **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`; the unpivoted + 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 + 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) + 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`; only `?geqr2` stays an elementwise `~eps` comparison. - **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..18331ba 100644 --- a/docs/building.md +++ b/docs/building.md @@ -3,9 +3,25 @@ ## 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'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 @@ -13,13 +29,23 @@ 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 (its configure needs a Fortran compiler, `brew install gcc`), +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`) `-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 +151,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. | +| `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/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 dc0ae8e..580b499 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_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 +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 6ed04a1..3fcf7ea 100644 --- a/docs/cbk_dsytrfnp_compact_design.md +++ b/docs/cbk_dsytrfnp_compact_design.md @@ -389,9 +389,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, @@ -439,7 +442,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 64af39e..8bd69d7 100644 --- a/docs/cbk_dtrsm_compact_design.md +++ b/docs/cbk_dtrsm_compact_design.md @@ -248,10 +248,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 @@ -260,9 +261,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..892add9 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_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. | | `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..bb8c90b 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); @@ -137,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) @@ -147,13 +151,15 @@ 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); - 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_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..07b15c6 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; @@ -67,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)); } @@ -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,131 @@ 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: 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; + gen_boosted(A, s == DENSE ? T(2) : T(0)); + if (s == DENSE) { + 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_orgqr(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, 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)); + } + } + + const double rtol_res = 20.0 * n * eps, rtol_orth = 100.0 * n * eps; + 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" + : 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; +} + +// 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 @@ -161,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; @@ -228,5 +359,10 @@ 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 + fails += run_invariant_set(); + fails += run_invariant_set(); + 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..029c27c --- /dev/null +++ b/tests/test_lapack_util.hpp @@ -0,0 +1,405 @@ +// 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_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 +// 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, 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 +// 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) +// 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 and ACCELERATE_NEW_LAPACK (the LAPACK 3.9.1 +// interface the shim is built for) are on the LAPACKE::LAPACKE target. +#include +#include +#else +#include +#include +#endif + +#include +#include +#include +#include + +namespace cbk::test { + +// ----------------------- LAPACKE / CBLAS, by scalar type ------------- +// 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; + +// 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 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_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_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_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); } \ +}; +// 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 (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 +// 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(M.ld(), 1); +} + +template int lapack_layout(const Mv &M) +{ + return is_rowmajor(M) ? LAPACK_ROW_MAJOR : LAPACK_COL_MAJOR; +} + +template CBLAS_ORDER 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; +} + +// 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), use_(M) + { + if (is_rowmajor(M) == rowmajor) return; + buf_.resize((std::size_t)M.rows * M.cols); + use_ = + mat_view(buf_.data(), M.rows, M.cols, rowmajor ? M.cols : M.rows, rowmajor); + copy_matrix(M, use_); + } + ~Staged() + { + if (!buf_.empty()) 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_, use_; + std::vector buf_; +}; + +// ----------------------- 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 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) +{ + 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_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)); + 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 +// reflectors in columns 0..k-1 of A (m x n, m >= n >= k). +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); +} + +// ?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_potrf(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 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; + 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(A0), trans, m, n, nrhs, A0.data, lapack_ld(A0), + Bs.data(), Bs.ld()); + + if (m >= n) + lapack::geqrf(lapack_layout(A), m, n, A.data, lapack_ld(A), tau); + else + 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 +// 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; + 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); + 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..db40bf6 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_orgqr, 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; @@ -91,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); @@ -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..7324067 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); @@ -94,19 +53,20 @@ 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); 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..d505fa7 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,61 +28,26 @@ #include #include -#include "test_compact_util.hpp" // compact, frand, gen_spd, MatrixBatch, pack/unpack +#include "test_lapack_util.hpp" // compact, ref_potrf, 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)); + ref_potrf(uplo, Aref.view(idx)); } // pack the full symmetric A, factor with the routine under test, unpack @@ -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 - 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); - 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); } @@ -245,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)); } } @@ -254,23 +222,26 @@ 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); 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; - e_spd = std::max(e_spd, (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)); + } 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 +289,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..d05bdd4 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,42 @@ 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. 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'); + 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); + 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( + "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 +403,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