diff --git a/.github/actions/build-randlapack-stack/action.yml b/.github/actions/build-randlapack-stack/action.yml new file mode 100644 index 0000000..f720d3e --- /dev/null +++ b/.github/actions/build-randlapack-stack/action.yml @@ -0,0 +1,128 @@ +name: build-randlapack-stack +description: > + Build and install the RandLAPACK stack (BLAS++, Random123, LAPACK++, RandLAPACK) + as siblings of the checked-out bindings workspace, so the bindings can + find_package(RandLAPACK). Linux + macOS. Mirrors the dependency bring-up in + RandLAPACK's own core-linux.yaml / core-macos.yaml. + +inputs: + randlapack-ref: + description: RandLAPACK git ref to build the bindings against. + # Bump to main once FunNystromPP lands there; for now the branch that has it. + default: funnystrompp-v2-baseline + +runs: + using: composite + steps: + - name: base toolchain + BLAS + shell: bash + run: | + if [ "$RUNNER_OS" = "Linux" ]; then + export DEBIAN_FRONTEND=noninteractive + # Intel MKL (ILP64) is the BLAS on Linux, NOT OpenBLAS: the bindings run + # inside MATLAB, whose runtime BLAS is MKL ILP64 (64-bit ints). blaspp must + # match that ABI or MATLAB's MKL services the MEX's BLAS calls with the wrong + # integer width (garbage leading dimension -> "Parameter N incorrect on entry + # to DGEMM" -> segfault). This mirrors the working local dev setup. + sudo apt-get update -qq + sudo apt-get install -qq -y git-core gcc g++ gfortran cmake m4 \ + libgtest-dev wget gpg + wget -qO- https://apt.repos.intel.com/intel-gpg-keys/GPG-PUB-KEY-INTEL-SW-PRODUCTS.PUB \ + | sudo gpg --dearmor -o /usr/share/keyrings/oneapi-archive-keyring.gpg + echo "deb [signed-by=/usr/share/keyrings/oneapi-archive-keyring.gpg] https://apt.repos.intel.com/oneapi all main" \ + | sudo tee /etc/apt/sources.list.d/oneAPI.list > /dev/null + sudo apt-get update -qq + sudo apt-get install -qq -y intel-oneapi-mkl-devel + # Source Intel's own MKL env script (robust across versions and the + # 'latest' symlink) so MKLROOT / LIBRARY_PATH / CPATH / LD_LIBRARY_PATH are + # set correctly, then hand them to later steps AND the calling job (build/ + # run MEX) via $GITHUB_ENV. blaspp's link test needs the MKL libs on + # LIBRARY_PATH and its headers on CPATH; the MEX needs the ILP64 MKL on + # LD_LIBRARY_PATH at runtime. mkl vars.sh is MKL-only (no Intel iomp5): we + # use mkl_gnu_thread + system libgomp, so MATLAB's runtime stays unpolluted. + MKL_ENV="$(ls /opt/intel/oneapi/mkl/*/env/vars.sh | sort -V | tail -1)" + source "$MKL_ENV" + MKL_LIB="$(dirname "$(find "$MKLROOT" -name 'libmkl_core.so*' | head -1)")" + { + echo "MKLROOT=$MKLROOT" + echo "MKL_LIB=$MKL_LIB" + echo "CPATH=$CPATH" + echo "LIBRARY_PATH=$LIBRARY_PATH" + echo "LD_LIBRARY_PATH=$LD_LIBRARY_PATH" + } >> "$GITHUB_ENV" + elif [ "$RUNNER_OS" = "macOS" ]; then + # OpenBLAS, NOT Apple Accelerate: RandLAPACK refuses to run BQRRP when + # blaspp is linked against Accelerate (rl guard), and MKL has no Apple + # Silicon build. macOS two-level namespaces keep the MEX bound to its own + # OpenBLAS at runtime (no MATLAB-BLAS interposition), so LP64 is fine here. + brew install googletest openblas + OB="$(brew --prefix openblas)" + { + echo "OPENBLAS_LIB=$OB/lib" + echo "LIBRARY_PATH=$OB/lib" + echo "LD_LIBRARY_PATH=$OB/lib" + echo "DYLD_LIBRARY_PATH=$OB/lib" + echo "CPATH=$OB/include" + } >> "$GITHUB_ENV" + fi + + - name: clone RandLAPACK (+ RandBLAS submodule) + shell: bash + run: | + cd "${{ github.workspace }}/.." + git clone https://github.com/BallisticLA/RandLAPACK.git + cd RandLAPACK + git checkout "${{ inputs.randlapack-ref }}" + git submodule update --init --recursive + + - name: build + install BLAS++ + shell: bash + run: | + cd "${{ github.workspace }}/.." + git clone https://github.com/icl-utk-edu/blaspp.git + # Linux: force ILP64 so blaspp calls MKL with 64-bit ints, matching MATLAB's + # runtime MKL (MKLROOT is set, so blaspp auto-detects MKL). macOS: force + # OpenBLAS (LP64) so blaspp does NOT pick Accelerate, which RandLAPACK's BQRRP + # guard rejects; two-level namespaces keep the MEX on its own OpenBLAS. + BLASPP_OPTS="" + if [ "$RUNNER_OS" = "Linux" ]; then + BLASPP_OPTS="-Dblas_int=int64" + elif [ "$RUNNER_OS" = "macOS" ]; then + BLASPP_OPTS="-Dblas=OpenBLAS" + fi + cmake -S blaspp -B blaspp-build -DCMAKE_BUILD_TYPE=Release \ + -DCMAKE_INSTALL_PREFIX="$PWD/blaspp-install" -Dbuild_tests=OFF $BLASPP_OPTS + cmake --build blaspp-build -j2 --target install + + - name: install Random123 + shell: bash + run: | + cd "${{ github.workspace }}/.." + git clone https://github.com/DEShawResearch/Random123.git + if [ "$RUNNER_OS" = "Linux" ]; then + make -C Random123 prefix="$PWD/Random123-install" install-include + else + mkdir -p Random123-install/include + cp -rp Random123/include/Random123 Random123-install/include/ + fi + + - name: build + install LAPACK++ + shell: bash + run: | + cd "${{ github.workspace }}/.." + git clone https://github.com/icl-utk-edu/lapackpp.git + cmake -S lapackpp -B lapackpp-build -DCMAKE_BUILD_TYPE=Release \ + -Dblaspp_DIR="$PWD/blaspp-install/lib/cmake/blaspp" \ + -DCMAKE_INSTALL_PREFIX="$PWD/lapackpp-install" -Dbuild_tests=OFF + cmake --build lapackpp-build -j2 --target install + + - name: build + install RandLAPACK (headers + cmake config; no tests built) + shell: bash + run: | + cd "${{ github.workspace }}/.." + cmake -S RandLAPACK -B RandLAPACK-build -DCMAKE_BUILD_TYPE=Release \ + -Dblaspp_DIR="$PWD/blaspp-install/lib/cmake/blaspp" \ + -Dlapackpp_DIR="$PWD/lapackpp-install/lib/cmake/lapackpp" \ + -DRandom123_DIR="$PWD/Random123-install/include/" \ + -DCMAKE_INSTALL_PREFIX="$PWD/RandLAPACK-install" + cmake --build RandLAPACK-build -j2 --target install diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml new file mode 100644 index 0000000..7a556e4 --- /dev/null +++ b/.github/workflows/ci.yml @@ -0,0 +1,83 @@ +name: ci + +on: + push: + pull_request: + workflow_dispatch: + +# Tests the BINDING layer (the glue), not RandLAPACK's numerics — RandLAPACK has +# its own CI. We build the RandLAPACK stack (composite action, mirroring +# RandLAPACK's core-linux/core-macos), then build + exercise the bindings. +# +# MATLAB: matlab-actions/setup-matlab runs license-free on PUBLIC GitHub repos +# (this repo is public). Linux + macOS via the job matrix. + +jobs: + matlab: + strategy: + fail-fast: false + matrix: + os: [ubuntu-latest, macos-latest] + runs-on: ${{ matrix.os }} + steps: + - uses: actions/checkout@v4 + + - name: build RandLAPACK stack + uses: ./.github/actions/build-randlapack-stack + + - name: set up MATLAB + uses: matlab-actions/setup-matlab@v2 + + - name: build the MEX bindings + shell: bash + run: | + DEPS="$(cd "${{ github.workspace }}/.." && pwd)" + MATLAB_ROOT="$(dirname "$(dirname "$(command -v matlab)")")" + cmake -S . -B build \ + -DCMAKE_BUILD_TYPE=Release \ + -DRandLAPACK_DIR="$DEPS/RandLAPACK-install/lib/cmake/RandLAPACK" \ + -Dblaspp_DIR="$DEPS/blaspp-install/lib/cmake/blaspp" \ + -Dlapackpp_DIR="$DEPS/lapackpp-install/lib/cmake/lapackpp" \ + -DMatlab_ROOT_DIR="$MATLAB_ROOT" + cmake --build build -j2 + + - name: run MATLAB binding tests + uses: matlab-actions/run-command@v2 + with: + command: "addpath('matlab'); addpath('matlab/tests'); test_bqrrp; test_fun_nystrom_pp; disp('matlab binding tests passed');" + env: + # the MEX dynamically links blaspp/lapackpp; set BOTH search-path vars + # (Linux uses LD_LIBRARY_PATH, macOS uses DYLD_LIBRARY_PATH). The MEX also + # needs its vendor BLAS: Linux DT_NEEDEDs the ILP64 MKL (MKL_LIB), macOS + # loads OpenBLAS (OPENBLAS_LIB); both are exported by the stack action and + # prepended. The unused one is empty, collapsing to a leading ':' the + # loader ignores. + LD_LIBRARY_PATH: "${{ env.MKL_LIB }}:${{ github.workspace }}/../blaspp-install/lib:${{ github.workspace }}/../lapackpp-install/lib" + DYLD_LIBRARY_PATH: "${{ env.OPENBLAS_LIB }}:${{ github.workspace }}/../blaspp-install/lib:${{ github.workspace }}/../lapackpp-install/lib" + + python: + # Python bindings are a BQRRP-only skeleton for now; keep CI green while it + # matures (the job still reports build/test status). + continue-on-error: true + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v4 + + - name: build RandLAPACK stack + uses: ./.github/actions/build-randlapack-stack + + - uses: actions/setup-python@v5 + with: + python-version: '3.11' + + - name: build + test the Python bindings + shell: bash + run: | + DEPS="$(cd "${{ github.workspace }}/.." && pwd)" + export RandLAPACK_DIR="$DEPS/RandLAPACK-install/lib/cmake/RandLAPACK" + export blaspp_DIR="$DEPS/blaspp-install/lib/cmake/blaspp" + export lapackpp_DIR="$DEPS/lapackpp-install/lib/cmake/lapackpp" + export LD_LIBRARY_PATH="$DEPS/blaspp-install/lib:$DEPS/lapackpp-install/lib:${LD_LIBRARY_PATH:-}" + python -m pip install --upgrade pip pytest numpy pybind11 + python -m pip install ./python -v + pytest python/tests -q diff --git a/.gitignore b/.gitignore index 80624c4..b45b3d1 100644 --- a/.gitignore +++ b/.gitignore @@ -13,3 +13,6 @@ __pycache__/ *.egg-info/ .pytest_cache/ *.so + +# generated by CMake FetchContent (machine-specific paths) +matlab/benchmark/funnystrompp/external_paths.m diff --git a/CMakeLists.txt b/CMakeLists.txt index b41faba..187eed4 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -13,7 +13,37 @@ set(CMAKE_POSITION_INDEPENDENT_CODE ON) # Primary: find_package against an installed RandLAPACK. The user is expected # to set -DRandLAPACK_DIR=/path/to/install/lib/cmake/RandLAPACK when invoking # cmake. CPM fallback is a future addition. +# +# RandLAPACK_DIR must be the directory holding RandLAPACKConfig.cmake, but a +# common mistake is passing the RandNLA-project root or the install prefix +# instead. Accept those by searching the known install layouts beneath the +# given directory; anything else gets an actionable error instead of +# find_package's generic one. # --------------------------------------------------------------------------- +if(DEFINED RandLAPACK_DIR AND NOT EXISTS "${RandLAPACK_DIR}/RandLAPACKConfig.cmake") + file(GLOB _rl_configs + "${RandLAPACK_DIR}/lib/cmake/RandLAPACK/RandLAPACKConfig.cmake" + "${RandLAPACK_DIR}/lib64/cmake/RandLAPACK/RandLAPACKConfig.cmake" + "${RandLAPACK_DIR}/RandLAPACK-install/lib/cmake/RandLAPACK/RandLAPACKConfig.cmake" + "${RandLAPACK_DIR}/RandLAPACK-install/lib64/cmake/RandLAPACK/RandLAPACKConfig.cmake" + "${RandLAPACK_DIR}/install/RandLAPACK-install/lib/cmake/RandLAPACK/RandLAPACKConfig.cmake" + "${RandLAPACK_DIR}/install/RandLAPACK-install/lib64/cmake/RandLAPACK/RandLAPACKConfig.cmake") + if(_rl_configs) + list(GET _rl_configs 0 _rl_config) + get_filename_component(_rl_dir "${_rl_config}" DIRECTORY) + message(STATUS "RandLAPACK_DIR does not itself contain RandLAPACKConfig.cmake; " + "using ${_rl_dir} found beneath it") + set(RandLAPACK_DIR "${_rl_dir}" CACHE PATH "Directory containing RandLAPACKConfig.cmake" FORCE) + else() + message(FATAL_ERROR + "RandLAPACK_DIR=${RandLAPACK_DIR} does not contain RandLAPACKConfig.cmake, " + "and none was found in the usual install layouts beneath it. " + "If RandLAPACK's install.sh has not completed successfully yet, run it first -- " + "it prints the exact -DRandLAPACK_DIR to use. Otherwise locate the file with\n" + " find /install -name RandLAPACKConfig.cmake\n" + "and pass its directory.") + endif() +endif() find_package(RandLAPACK REQUIRED) message(STATUS "RandLAPACK found at ${RandLAPACK_DIR}") @@ -34,3 +64,15 @@ else() "${CMAKE_CURRENT_SOURCE_DIR}/matlab to the MATLAB path, " "but calls into bqrrp_mex will fail until the MEX is built.") endif() + +# --------------------------------------------------------------------------- +# Benchmark external baselines (off by default). When ON, FetchContent pulls +# the external reference implementations (Epperly, Persson) used by the +# funNystrom++ benchmark at configure time — NOT vendored in this repo. +# See matlab/benchmark/funnystrompp/CMakeLists.txt and THIRD_PARTY_LICENSES.md. +# --------------------------------------------------------------------------- +option(RANDLAPACK_BINDINGS_BENCHMARKS + "Fetch external benchmark baselines (Epperly, Persson) via FetchContent" OFF) +if(RANDLAPACK_BINDINGS_BENCHMARKS) + add_subdirectory(matlab/benchmark/funnystrompp) +endif() diff --git a/README.md b/README.md index 9f6f96f..c004924 100644 --- a/README.md +++ b/README.md @@ -2,9 +2,15 @@ MATLAB, and eventually Python, bindings for [RandLAPACK](https://github.com/BallisticLA/RandLAPACK). -## Scope (v0) +## Scope -Proof-of-concept binding for a single driver, `BQRRP` (randomized blocked QR with column pivoting), in single and double precision, real matrices, column-major. More drivers as their C++ APIs stabilize. Python next. +Bindings for selected RandLAPACK drivers, single and double precision, real matrices, column-major: + +* **`BQRRP`** (randomized blocked QR with column pivoting) — MATLAB + Python. +* **`FunNystromPP`** (matrix-function trace estimation, `tr f(A)`) — MATLAB. The + Python binding is deferred. + +More drivers as their C++ APIs stabilize. ## Requirements @@ -16,12 +22,34 @@ Proof-of-concept binding for a single driver, `BQRRP` (randomized blocked QR wit ```sh cmake -S . -B build \ + -DCMAKE_BUILD_TYPE=Release \ -DRandLAPACK_DIR=/path/to/RandLAPACK-install/lib/cmake/RandLAPACK \ -DMatlab_ROOT_DIR=/path/to/MATLAB cmake --build build -j ``` -If MATLAB is not found, MEX targets are skipped with a clear warning; the `.m` files still ship and can be used once the MEX is built. On success, `matlab/+randlapack/private/bqrrp_mex.` is in place. +Always pass `-DCMAKE_BUILD_TYPE=Release`: with no build type the MEX compiles +without optimization (and Ninja+MSVC even defaults to Debug, measured +2026-08-19), which is functionally correct but useless for timing. + +On Windows, run from an x64 developer shell and pass two extra things: +`-G Ninja` and `-DRANDLAPACK_RUNTIME_DLL_DIRS=` (the +RandLAPACK installer prints that directory). `-G Ninja` is REQUIRED, not a +preference: without it CMake selects the multi-config Visual Studio generator, +which appends a `Debug\`/`Release\` subdirectory to the MEX output path, so the +MEX lands in `matlab/+randlapack/private/Debug/` where MATLAB's package loader +never looks (verified 2026-08-19; Ninja ships with the VS "C++ CMake tools" +component). The DLL-dirs flag stages the BLAS runtime next to the MEX; without +it MATLAB fails to load the MEX with "The specified module could not be found". + +If MATLAB is not found, MEX targets are skipped with a clear warning; the `.m` files still ship and can be used once the MEX is built. On success, the per-driver MEX (`matlab/+randlapack/private/{bqrrp,fun_nystrom_pp}_mex.`) is in place. + +### Prebuilt release (no compiler, no CMake) + +To avoid building RandLAPACK and the MEX yourself, a release archive ships the +`.m` files plus a prebuilt `matlab/+randlapack/private/*.mexa64`. Unzip, +`addpath('.../matlab')`, and call the functions directly. Prebuilt binaries are +currently **Linux x86-64 only**; on macOS/Windows, build from source as above. ## Usage @@ -42,15 +70,21 @@ Two output modes, chosen by a trailing string (qr-style): Run `help randlapack.bqrrp` for the full signature, including the optional `b_sz`, `d_factor`, and RNG `state` arguments. -## Runtime notes - -On **Linux MATLAB R2026a + Ubuntu 24.04** (and similar new-glibc hosts), MATLAB ships an older `libstdc++.so.6` than your compiler's. MEX files built against the system `libstdc++` fail at load with `version GLIBCXX_3.4.32 not found`. Workaround: +For matrix-function trace estimation: -```sh -LD_PRELOAD=/usr/lib/x86_64-linux-gnu/libstdc++.so.6 matlab +```matlab +n = 1000; +G = randn(n, 200); A = G * G.' / 200; % symmetric PSD test matrix +est = randlapack.fun_nystrom_pp(A, 100, 50, 'Func', 'sqrt'); % k=100, s=50 ``` -Alias it in your shell profile to make it transparent. Not needed on macOS or older Ubuntu hosts. +Run `help randlapack.fun_nystrom_pp` for the full name-value signature (`Func`, `Q`, `LFAType`, `Depth`, `Sketch`, `Reorth`, …). + +## Runtime notes + +The MEX targets statically link `libstdc++`/`libgcc` (`-static-libstdc++ -static-libgcc`), so they load in any MATLAB regardless of the `GLIBCXX` version MATLAB ships — **no `LD_PRELOAD` is needed**. (Earlier versions required preloading the system `libstdc++.so.6`; that is no longer the case.) + +The MEX does dynamically link BLAS/LAPACK (via blaspp/lapackpp). When running outside the build environment, point `LD_LIBRARY_PATH` at the blaspp/lapackpp install `lib`/`lib64` dirs, and set `MKL_INTERFACE_LAYER=ILP64`, `MKL_THREADING_LAYER=GNU` if using MKL. ## Conventions diff --git a/THIRD_PARTY_LICENSES.md b/THIRD_PARTY_LICENSES.md new file mode 100644 index 0000000..8390613 --- /dev/null +++ b/THIRD_PARTY_LICENSES.md @@ -0,0 +1,46 @@ +# Third-Party Licenses / Acknowledgements + +The funNystrom++ benchmark (`matlab/benchmark/funnystrompp/`) compares +RandLAPACK's `FunNystromPP` against external reference implementations. **These +are not vendored in this repository.** They are fetched at build time via CMake +`FetchContent`, pinned to a commit each, only when the benchmark is enabled +(`-DRANDLAPACK_BINDINGS_BENCHMARKS=ON`). See +`matlab/benchmark/funnystrompp/CMakeLists.txt`. + +## Ethan Epperly — trace-estimation code (XTrace / XNysTrace / resphering) + +- Repository: https://github.com/eepperly/Ethan-Epperly-Thesis +- Pinned commit: `2d33ee40577caff18b16fed3454a147d08cbe16e` (branch `main`) +- License: **MIT**. + +``` +MIT License + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND ... +``` +(Full text in the repository's `LICENSE` file at the pinned commit.) + +## David Persson — funNystrom++ reference + +- Repository: https://github.com/davpersson/funNystrom +- Pinned commit: `0427bde1e3c238ce53cc184600482481490c2a48` (branch `master`) +- Files used: `Other/nystrom.m`, `Other/funnystrompp.m`, `Other/block_lanczos.m`. +- License: **none stated** at the time of writing (under default copyright, all + rights reserved). Used here **with attribution**, fetched (not redistributed) + for benchmark comparison only. **A license is being requested from the + author**; this acknowledgement is interim. +- Note: `block_lanczos.m` in that repository is Christopher Meyer's + HutchPlusPlus code, vendored by Persson with provenance noted in its header. + +We thank David Persson and Ethan Epperly for making their reference code +available. diff --git a/matlab/+randlapack/fun_nystrom_pp.m b/matlab/+randlapack/fun_nystrom_pp.m new file mode 100644 index 0000000..cd899d1 --- /dev/null +++ b/matlab/+randlapack/fun_nystrom_pp.m @@ -0,0 +1,522 @@ +function [est, t1, t2, times] = fun_nystrom_pp(A, k, Omega2, varargin) +%FUN_NYSTROM_PP Trace estimator tr(f(A)) via RandLAPACK FunNystromPP. +% Preferred four-method interface: randlapack.trace_estimate. +% LFAType scalar_auto/block_auto share target AutoEps and hard Budget, +% including the pilot. Both use sphere probes, independent pilot samples, +% the same rank/probe allocation, and Gauss output. Certification flags +% concern quadrature only, never the total stochastic trace error. +% +% est = randlapack.fun_nystrom_pp(A, k, s, ...) +% est = randlapack.fun_nystrom_pp(A, k, Omega2, ...) % explicit Phase-2 probes +% [est, t1, t2, times] = randlapack.fun_nystrom_pp(...) +% the 4th output is a wall-clock instrumentation struct (fields: +% marshal_in_ms (for sparse A including the wrapper's compressed-column +% conversion), phase1_ms, phase2_ms, fafun_ms, assembly_ms, +% specrec_ms, nystrom_us (1x11; only slots 1,2,3,7,11 populated, slot 7 +% = the whole spectral-recovery block), lfa_us (1x6; slot 6 = reorth +% time, real for scalar/block/block_qfa, 0 for the basis-free +% scalar_qfa/auto/adaptive), d_used (NaN for LFAType 'exact', which has +% no Lanczos depth), oracle_mv, auto_k, auto_s, +% probe_mv (also populated for 'adaptive', reusing these same generic +% fields — see LFAType 'adaptive' below), plus path-specific fields +% that are NaN when their path did not run: tr_U, tr_L, certified +% ('block_qfa'/'adaptive' Gauss/Radau traces + bracket flag) and +% probe_ms, probe_converged, phase2_certified ('auto'/'adaptive' +% depth-probe wall-clock + certification flags), and five numerical- +% degradation diagnostics that are always populated for a compiled +% call: ritz_clamped (negative Ritz values clamped to zero before f, +% summed over every Lanczos object the call used), rank_deficient_steps +% and min_diag_ratio (block QR steps whose smallest |diag R| fell below +% 16 eps times the largest, and the smallest such ratio seen; Inf when +% no block QR ran), nystrom_clamped (head eigenvalues driven negative +% by the shift removal and clamped), and bracket_evaluated (true only +% when a Gauss/Gauss-Radau pair was actually computed; when false, +% tr_U == tr_L carries no information). All are cumulative over the +% call and nothing branches on them. See fun_nystrom_pp_mex.cc for the +% full slot definitions. +% +% A n x n matrix, single or double, REQUIRED symmetric. The MEX +% reads only the upper triangle and mirrors it into the lower +% triangle unconditionally (the lower triangle is never read, not +% "used"). By default this wrapper checks symmetry with a +% tolerance-gated maxabs(A-A') <= sqrt(eps(class(A)))*maxabs(A) +% test and raises 'randlapack:fun_nystrom_pp:NotSymmetric' on +% failure; pass 'SkipSymCheck', true to skip the O(n^2) check for +% callers who already know A is symmetric (e.g. a large-matrix +% benchmark loop). +% +% VECTOR MODE: A may instead be an n-vector (either orientation), +% read as the diagonal entries of A = diag(A). The n x n matrix is +% never formed, in MATLAB or in the MEX, so this is the way to run +% at sizes where the dense matrix does not fit. The symmetry check +% is bypassed (there is no lower triangle), and LFAType 'adaptive' +% switches its Rademacher probes to the sphere, because a +-1 +% quadratic form is exact on a diagonal matrix and would make that +% tier's error degenerate rather than sampled. +% k Phase-1 rank (SCALAR). The Phase-1 sketch is a SASO +% generated INSIDE RandLAPACK::NystromEVD from (SketchSeed, VecNnz), +% matching the paper's Algorithm 1 line 1. Explicit Omega1 matrices +% are no longer accepted (the dense-sketch mode was removed from +% the kernel). +% arg3 Phase-2 Hutchinson probes. Pass a SCALAR s and the n x s sphere +% probes are sampled internally (seed = SketchSeed + 1000); OR pass +% an explicit n x s matrix (escape hatch for validation, e.g. +% feeding every estimator in a comparison the SAME probes). +% Skipped when k == n (then arg3 may be a 0/[] placeholder). Also +% never read for LFAType 'auto'/'adaptive' (those tiers derive +% their own probes from a certified depth probe): the MEX skips +% the marshal of an explicit Omega2 entirely in that case, so it +% is validated here but never copied. +% +% A knob explicitly passed but not consumed by the chosen 'LFAType' raises +% warning 'randlapack:fun_nystrom_pp:ignored_knob' (fired HERE, using +% inputParser's UsingDefaults so a knob left at its own default is silent +% while one explicitly repeated at its default value still warns; see +% each 'LFAType' entry below for what it reads). This replaces a prior +% MEX-side value-comparison warning that could not tell "the caller left +% this alone" from "the caller passed the default value on purpose", and +% that fired unconditionally whenever a caller forwarded a knob's value +% verbatim (e.g. every harness call). +% +% Optional Name-Value pairs: +% 'Func' char in {'sqrt', 'log', 'poly', 'effdim', 'square', 'identity', 'entropy'} +% entropy = -x*log(x), with f(0)=0; the Gauss/Gauss-Radau +% certificate applies from depth 2 (the scalar_auto and +% block_auto tiers accept it). +% (default 'sqrt'). 'poly' is f(x) = x(x + PolyLambda); +% 'effdim' is f(x) = x/(x + PolyLambda) (the "effective +% dimension" function; operator monotone). +% 'Q' subspace-iter count (default 2) +% 'PolyLambda' lambda used by Func 'poly' and 'effdim' (default 10) +% 'LFAType' {'exact', 'scalar', 'scalar_qfa', 'block', 'block_qfa', +% 'auto', 'adaptive'} oracle for f(A)*X (default 'block' — +% the matrix-free Krylov oracle). +% 'block_qfa' = block Lanczos-QFA: forms the s×s quadratic +% form Ω₂ᵀf(A)Ω₂ directly (no f(A)·Ω₂ mapback); cheapest. +% 'Adaptive' selects its depth rule: 0 = fixed Depth; +% 1 = Radau-certified stop (the block Gauss/Gauss-Radau +% bracket closes within 'AdaptiveTol' — the same certified +% meaning as scalar_qfa's adaptive mode; no delay window); +% 2 = legacy window rule (honors AdaptiveDelay/AdaptiveMin). +% With Adaptive = 1, 'RadauReturn' picks the returned value +% (0 = block Gauss, 1 = (Gauss+Radau)/2 midpoint) and the +% times struct reports the tr_U / tr_L bracket plus the +% certified flag. tr_U >= tr_L holds in the operator- +% monotone regime; the default MaxBoth scaling does not +% assume it. +% 'scalar_qfa' = scalar Lanczos-QFA: the per-probe quadratic +% forms directly (basis-free, O(n·s) memory). With +% 'Adaptive',1 each probe stops at its own depth via the +% Gauss-Radau certificate and 'AdaptiveTol' is a CERTIFIED +% per-probe relative error (no window/floor); Depth is the +% cap. Reorth/AdaptiveDelay/AdaptiveMin are ignored. The +% times struct reports d_used (max per-probe depth) and +% oracle_mv (actual Σ per-probe matvecs). +% 'exact' builds a full eigendecomposition of A (O(n^3)) and +% is intended for validation/reference use, not production. +% 'scalar' runs one Lanczos recurrence per probe (equivalent +% to Lanczos quadrature on each quadratic form). +% 'auto' = knob-free tier: pass 'Budget' (total A-matvec +% budget) and 'AutoEps'; the driver picks k, s, and the +% oracle depth itself (positional k/s become placeholders and +% Q/Depth/Reorth/Adaptive* are ignored: q is hardcoded to 1 +% internally; this wrapper warns if any of these were +% explicitly passed), running the certified +% scalar QFA for both the depth probe and Phase 2. The times +% struct reports the choices (auto_k, auto_s, d_used, +% probe_mv, oracle_mv), the depth-probe wall-clock +% (probe_ms; the probe runs before phase1_ms's clock, so +% phase1_ms + phase2_ms excludes it) and the certification +% flags (probe_converged, phase2_certified); spend closes +% as an upper bound +% probe_mv + q*auto_k + oracle_mv <= Budget (Phase 1 costs +% q*auto_k matvecs; q = 1 here). A non-positive Budget (< 1) +% raises error id 'randlapack:fun_nystrom_pp:Budget'; an +% otherwise-infeasible Budget (too small once feasibility is +% checked against eps/n) raises +% 'randlapack:fun_nystrom_pp:infeasibleBudget'. +% 'adaptive' = fully knob-free, EPS-TARGETED tier: no +% upfront Budget; pass 'AutoEps' (reused here as the target +% accuracy eps) and the driver picks k, s, and the oracle +% depth itself from a certified BLOCK Gauss-Radau depth +% probe (positional k/s become placeholders and +% Q/Depth/Reorth/Adaptive* are ignored, same as 'auto'). +% 'AutoDepthCap' caps the probe depth, shared with the +% 'auto' tier. 'AdaptiveMatvecCap' (default 0 = no cap) +% optionally bounds the total matvec spend. The times +% struct reports the choices in the SAME generic fields as +% 'auto' (auto_k, auto_s, d_used, probe_mv, oracle_mv, +% probe_ms, probe_converged, phase2_certified) plus tr_U / +% tr_L / certified (the block oracle's final Gauss/ +% Gauss-Radau trace bracket and certification flag, the +% same fields 'block_qfa' populates). An infeasible +% AdaptiveMatvecCap, or an eps/n regime the block-Krylov +% probe count cannot fund, raises the SAME error id as +% 'auto': 'randlapack:fun_nystrom_pp:infeasibleBudget'. +% 'Depth' Lanczos depth for 'scalar' / 'scalar_qfa' / 'block' / +% 'block_qfa' LFAType (default 200 for scalar/scalar_qfa, +% 20 for block AND block_qfa; a CAP in the adaptive QFA +% modes; ignored for 'exact') +% 'RadauReturn' block_qfa + Adaptive 1 only: value returned on a certified +% stop. 0 = block Gauss (default; matches the scalar +% oracle), 1 = (Gauss + Radau)/2 midpoint (for operator- +% monotone f the two quadratures err on opposite sides, so +% the midpoint halves the one-sided Gauss bias for free). +% 'AutoDepthCap' LFAType 'auto' only: fixed cap on the depth probe +% (default 0 = no fixed cap; the probe is then bounded only +% by n and by AutoProbeFrac). +% 'AutoProbeFrac' LFAType 'auto' only: fraction of the matvec Budget the +% depth probe may spend, in (0, 1) (default 0.125). +% 'AdaptiveKConst' / 'AdaptiveSConst' LFAType 'adaptive' only: the leading +% constants of the split rule k = c_k*sqrt(t)/eps, +% s = max(1, c_s/(sqrt(t)*eps)). Both default to 1, which is +% the paper's unit-constant allocation; the asymptotic rate is +% unchanged by either constant, only the allocation at a given +% (t, eps). +% 'CapRankFraction' Target-based tiers ('adaptive', 'scalar_auto', +% 'block_auto'): maximum rank share of post-pilot work, +% applied only when the target allocation exceeds the cap. +% In (0,1]; default 1 preserves historical rank-first clipping. +% 'QuadratureFraction' Same tiers: multiply AutoEps by this fraction for +% pilot and final Lanczos stopping, leaving the rank/probe +% target unchanged. In (0,1], default 1. This does not certify +% total stochastic error. Intended for controlled diagnostics. +% 'SpendCap' 'scalar_auto'/'block_auto' only, default false. True spends +% the whole Budget: rank and probe count are chosen to use the +% cap (at least 4 probes), and every Phase-2 probe runs to the +% probe's depth with no early stop, so times.phase2_checked +% is 0 and times.phase2_certified, times.certified and +% times.tr_L are NaN (nothing was checked). +% 'FirstRowQL' Optional scalar projected evaluator, default false. Single +% precision and nonconverged QL solves use stevd. Returned +% times record the request and any pilot/final fallback. +% 'AdaptiveMatvecCap' LFAType 'adaptive' only: optional total matvec cap +% (default 0 = no cap). Ignored by every other LFAType. +% 'Sketch' 'saso' (default): the kernel-internal sparse SASO (density VecNnz). +% 'gaussian': a dense Gaussian sketch applied by ordinary matrix +% products (comparison series; VecNnz is ignored). Other values are +% rejected before calling MEX. +% 'VecNnz' nonzeros per ROW of the SASO sketch (default 8; +% 0 = auto, resolved to ~log(k) inside the kernel) +% 'SketchSeed' RNG seed for the Phase-1 sketch (default 42) +% 'SkipSymCheck' skip the default symmetry validation of A (default +% false). Set true only when the caller already knows A is +% symmetric (e.g. a large-n benchmark loop where the O(n^2) +% check is unwanted overhead); does not change what the MEX +% does with A (upper triangle mirrored into lower either way). +% +% Returns: +% est trace estimate t1 + t2 +% t1 Phase 1 contribution sum f(lambda_hat_i) +% t2 Phase 2 Hutchinson correction (0 when k == n) +% +% Algorithm: two-phase funNystrom++ (rank-k shifted Nystrom approximation +% per arXiv:2508.21189 Alg. 2 + Hutchinson correction on the residual). +% MEX wrapper over RandLAPACK::FunNystromPP, T matching the class of A. +% NB the Phase-1 sketch is drawn by RandBLAS inside the kernel, so runs are +% reproducible at fixed SketchSeed but are NOT sketch-for-sketch identical +% to the Persson MATLAB reference (which draws its own Gaussian sketch). + + % --- Required-argument validation (friendly MATLAB-side errors) --- + % A is either a square n x n matrix or an n-vector read as diag(A). The + % 'square' attribute would reject the vector, so the shape is checked here + % instead; 'numel' is the only safe size read (size(A,1) is 1 for a row + % vector). + % Sparse A (double; MATLAB has no sparse single) must hold BOTH triangles: the MEX applies the whole matrix with + % RandBLAS sparse products and never infers a triangle. + if issparse(A) + validateattributes(A, {'double'}, {'2d', 'real', 'finite'}, mfilename, 'A', 1); + else + validateattributes(A, {'single', 'double'}, ... + {'2d', 'real', 'finite', 'nonsparse'}, ... + mfilename, 'A', 1); + end + vector_mode = ~issparse(A) && isvector(A) && ~isscalar(A); + if ~vector_mode && size(A, 1) ~= size(A, 2) + error('randlapack:fun_nystrom_pp:A_shape', ... + ['A must be a square n x n matrix or an n-vector of diagonal ' ... + 'entries; got %d x %d.'], size(A, 1), size(A, 2)); + end + cls = class(A); + if vector_mode + n = numel(A); + else + n = size(A, 1); + end + + if ~isscalar(k) + error('randlapack:fun_nystrom_pp:Omega1Explicit', ... + ['arg 2 must be a scalar rank k. Explicit Omega1 matrices are no ' ... + 'longer supported: the Phase-1 sketch is a SASO ' ... + 'generated inside RandLAPACK (control it via ''SketchSeed'' and ' ... + '''VecNnz'').']); + end + validateattributes(k, {'numeric'}, {'integer', 'positive', '<=', n}, ... + mfilename, 'k', 2); + k = double(k); + a1 = k; + + if k == n + a2 = double(0); % Phase 2 skipped; MEX leaves arg 3 unread + elseif isscalar(Omega2) + validateattributes(Omega2, {'numeric'}, {'integer', 'positive'}, ... + mfilename, 's', 3); + a2 = double(Omega2); + elseif isempty(Omega2) + error('randlapack:fun_nystrom_pp:Omega2Empty', ... + ['Omega2/s may be empty only when k == n (Phase 2 skipped). ' ... + 'Here k = %d and n = %d; pass a scalar s or an n x s sketch.'], k, n); + else + validateattributes(Omega2, {'single', 'double'}, ... + {'2d', 'real', 'finite', 'nonsparse'}, ... + mfilename, 'Omega2', 3); + if size(Omega2, 1) ~= n + error('randlapack:fun_nystrom_pp:Omega2Shape', ... + 'Omega2 must have n = %d rows to match A; got %d.', n, size(Omega2, 1)); + end + a2 = cast(Omega2, cls); + end + + % --- Name-Value options --- + p = inputParser; + addParameter(p, 'Func', 'sqrt', @(x) ischar(x) || isstring(x)); + addParameter(p, 'Q', 2, @(x) isnumeric(x) && isscalar(x) && x >= 1); + addParameter(p, 'PolyLambda', 10, @(x) isnumeric(x) && isscalar(x)); + addParameter(p, 'LFAType', 'block', @(x) ischar(x) || isstring(x)); + addParameter(p, 'Depth', [], @(x) isempty(x) || (isnumeric(x) && isscalar(x) && x >= 1)); + addParameter(p, 'Sketch', 'saso', @(x) ischar(x) || isstring(x)); + % VecNnz 0 = auto (~log k, resolved inside the kernel). + addParameter(p, 'VecNnz', 8, @(x) isnumeric(x) && isscalar(x) && x >= 0); + addParameter(p, 'SketchSeed', 42, @(x) isnumeric(x) && isscalar(x) && x >= 0); + % C++ side treats Reorth as a plain boolean (full reorth vs none); restrict + % to {0, 1} like the neighboring Adaptive/RadauReturn guards so e.g. + % Reorth = -3 (which would otherwise silently behave like Reorth = 1) is + % rejected instead. + addParameter(p, 'Reorth', 1, @(x) isnumeric(x) && isscalar(x) && any(x == [0 1])); + % QFA depth rule: 0 = fixed Depth; 1 = Radau-certified; 2 = legacy window + % (block_qfa only; scalar_qfa treats any nonzero as certified adaptive). + addParameter(p, 'Adaptive', 0, @(x) isnumeric(x) && isscalar(x) && any(x == [0 1 2])); + addParameter(p, 'AdaptiveTol',1e-2, @(x) isnumeric(x) && isscalar(x) && x > 0); + % block_qfa certified return value: 0 = block Gauss, 1 = midpoint. + addParameter(p, 'RadauReturn', 0, @(x) isnumeric(x) && isscalar(x) && any(x == [0 1])); + % Certificate window for block_qfa adaptive (0 = MEX/library default). The + % first convergence test is at depth AdaptiveMin + AdaptiveDelay. + addParameter(p, 'AdaptiveDelay', 0, @(x) isnumeric(x) && isscalar(x) && x >= 0); + addParameter(p, 'AdaptiveMin', 0, @(x) isnumeric(x) && isscalar(x) && x >= 0); + % Knob-free tier (LFAType 'auto'): total A-matvec budget + target accuracy; + % the driver picks k, s, and the oracle depth (positional k/s and + % Depth/Reorth/Adaptive are ignored in this mode). + addParameter(p, 'Budget', 0, @(x) isnumeric(x) && isscalar(x) && x >= 0); + addParameter(p, 'AutoEps', 1e-3, @(x) isnumeric(x) && isscalar(x) && x > 0 && x < 1); + % Auto-tier probe knobs: a fixed cap on the probe depth (0 = no fixed cap) + % and the fraction of the Budget the probe may spend, in (0, 1). + addParameter(p, 'AutoDepthCap', 0, @(x) isnumeric(x) && isscalar(x) && x >= 0); + addParameter(p, 'AutoProbeFrac', 0.125, @(x) isnumeric(x) && isscalar(x) && x > 0 && x < 1); + % Eps-targeted tier (LFAType 'adaptive') only: optional total matvec cap + % (0 = no cap). + addParameter(p, 'AdaptiveMatvecCap', 0, @(x) isnumeric(x) && isscalar(x) && x >= 0); + % Eps-targeted tier only: leading constants of the (k, s) split rule; 1 = the + % paper's unit constants. + addParameter(p, 'AdaptiveKConst', 1, @(x) isnumeric(x) && isscalar(x) && x > 0); + addParameter(p, 'AdaptiveSConst', 1, @(x) isnumeric(x) && isscalar(x) && x > 0); + % Allocation diagnostics for target-driven modes. Defaults reproduce the + % September 15 policy; high-level trace_estimate chooses the tested policy. + addParameter(p, 'CapRankFraction', 1, @(x) isnumeric(x) && isscalar(x) && x>0 && x<=1); + addParameter(p, 'QuadratureFraction', 1, @(x) isnumeric(x) && isscalar(x) && x>0 && x<=1); + % Equivalent scalar projected evaluator; float and failed QL use stevd. + addParameter(p, 'SpendCap', false, @(x) isscalar(x) && (islogical(x) || isnumeric(x)) && any(x==[0 1])); + addParameter(p, 'FirstRowQL', false, @(x) isscalar(x) && (islogical(x) || isnumeric(x)) && any(x==[0 1])); + % Default-on symmetry check (see the 'A' doc above); set true to skip it. + addParameter(p, 'SkipSymCheck', false, @(x) isscalar(x) && (islogical(x) || isnumeric(x))); + parse(p, varargin{:}); + + func = char(p.Results.Func); + q = double(p.Results.Q); + pl = double(p.Results.PolyLambda); + lfa_type = char(p.Results.LFAType); + sketch = lower(char(p.Results.Sketch)); + if ~any(strcmpi(sketch,{'saso','gaussian'})) + error('randlapack:fun_nystrom_pp:unsupportedSketch', ... + ['Sketch must be ''saso'' (sparse, default) or ''gaussian'' (dense, for comparison). ' ... + 'Use VecNnz to change its density; other Sketch values previously had no effect.']); + end + vec_nnz = double(p.Results.VecNnz); + sk_seed = double(p.Results.SketchSeed); + reorth = double(p.Results.Reorth); + adaptive = double(p.Results.Adaptive); % QFA types only; Depth is the cap + adapt_tol = double(p.Results.AdaptiveTol); + adapt_dl = double(p.Results.AdaptiveDelay); + adapt_mn = double(p.Results.AdaptiveMin); + radau_ret = double(p.Results.RadauReturn); + budget = double(p.Results.Budget); + auto_eps = double(p.Results.AutoEps); + auto_dcap = double(p.Results.AutoDepthCap); + auto_pfr = double(p.Results.AutoProbeFrac); + adapt_mvc = double(p.Results.AdaptiveMatvecCap); + adapt_kc = double(p.Results.AdaptiveKConst); + adapt_sc = double(p.Results.AdaptiveSConst); + cap_rank_fraction = double(p.Results.CapRankFraction); + quadrature_fraction = double(p.Results.QuadratureFraction); + skip_sym_check = logical(p.Results.SkipSymCheck); + + % --- Ignored-knob warnings (moved here from the MEX; see the doc note + % above 'Optional Name-Value pairs'). The consumption matrix below is + % built from this docstring's own 'LFAType' entries plus + % rl_fun_nystrom_pp.hh (q hardcoded to 1 and Adaptive*/Depth/Reorth + % unread by the 'auto'/'adaptive' overloads; scalar_qfa's own + % Reorth/AdaptiveDelay/AdaptiveMin no-op documented above; + % AdaptiveMatvecCap read only by the 'adaptive' overload's matvec_cap + % argument). p.UsingDefaults lists knobs NOT explicitly passed, so this + % fires iff a knob was explicitly passed (any value, including its own + % default) and the chosen LFAType does not consume it. --- + switch lfa_type + case 'exact' + ignored = {'Depth', 'Reorth', 'Adaptive', 'AdaptiveTol', ... + 'AdaptiveDelay', 'AdaptiveMin', 'RadauReturn', 'AdaptiveMatvecCap', ... + 'AdaptiveKConst', 'AdaptiveSConst'}; + case 'scalar' + ignored = {'Adaptive', 'AdaptiveTol', 'AdaptiveDelay', 'AdaptiveMin', ... + 'RadauReturn', 'AdaptiveMatvecCap', 'AdaptiveKConst', 'AdaptiveSConst'}; + case 'scalar_qfa' + ignored = {'Reorth', 'AdaptiveDelay', 'AdaptiveMin', 'RadauReturn', ... + 'AdaptiveMatvecCap', 'AdaptiveKConst', 'AdaptiveSConst'}; + case 'block' + ignored = {'Adaptive', 'AdaptiveTol', 'AdaptiveDelay', 'AdaptiveMin', ... + 'RadauReturn', 'AdaptiveMatvecCap', 'AdaptiveKConst', 'AdaptiveSConst'}; + case 'block_qfa' + ignored = {'AdaptiveMatvecCap', 'AdaptiveKConst', 'AdaptiveSConst'}; + case 'auto' + ignored = {'Q', 'Depth', 'Reorth', 'Adaptive', 'AdaptiveTol', ... + 'AdaptiveDelay', 'AdaptiveMin', 'RadauReturn', 'AdaptiveMatvecCap', ... + 'AdaptiveKConst', 'AdaptiveSConst'}; + case {'adaptive', 'scalar_auto', 'block_auto'} + ignored = {'Q', 'Depth', 'Reorth', 'Adaptive', 'AdaptiveTol', ... + 'AdaptiveDelay', 'AdaptiveMin', 'RadauReturn'}; + otherwise + % Unknown lfa_type: leave the diagnostic to the MEX's own + % 'randlapack:fun_nystrom_pp_mex:lfa_type' error. + ignored = {}; + end + if ~ismember(lfa_type, {'adaptive','scalar_auto','block_auto'}) + ignored = [ignored, {'CapRankFraction','QuadratureFraction'}]; + end + spend_cap = logical(p.Results.SpendCap) && ismember(lfa_type, {'scalar_auto','block_auto'}); + if ~ismember(lfa_type, {'scalar_auto','block_auto'}) + ignored = [ignored, {'SpendCap'}]; + end + explicitly_ignored = ignored(~ismember(ignored, p.UsingDefaults)); + if ~isempty(explicitly_ignored) + warning('randlapack:fun_nystrom_pp:ignored_knob', ... + ['LFAType ''%s'' ignores: %s (explicitly passed but not consumed ' ... + 'by this tier -- see the ''LFAType'' doc above for what each tier reads).'], ... + lfa_type, strjoin(explicitly_ignored, ', ')); + end + % Explicit nonempty Omega2 MATRIX with auto/adaptive: both tiers derive + % their own Phase-2 probes from a certified depth probe and never read a + % caller-supplied Omega2 (the MEX skips the marshal entirely for these + % two tiers, below). Scalar arg3 is exempted: it is the conventional + % placeholder spelling used by every call site (including this + % project's own harness), not a sign of a caller expecting their probes + % to be used. + if any(strcmp(lfa_type, {'auto', 'adaptive', 'scalar_auto', 'block_auto'})) && ~isscalar(Omega2) && ~isempty(Omega2) + warning('randlapack:fun_nystrom_pp:ignored_knob', ... + ['explicit Omega2 matrix ignored for LFAType ''%s'': this tier derives ' ... + 'its own Phase-2 probes from a certified depth probe and never reads a ' ... + 'caller-supplied Omega2.'], lfa_type); + end + + % Same error id the MEX's own auto_budget < 1 guard raises (D moderate + % finding: previously the two checks used different, undocumented ids; + % unified so direct fun_nystrom_pp_mex(...) callers see the identical + % identifier this wrapper would have raised). + if any(strcmp(lfa_type, {'auto', 'scalar_auto', 'block_auto'})) && budget < 1 + error('randlapack:fun_nystrom_pp:Budget', ... + 'LFAType ''auto'' requires a positive ''Budget'' (total A-matvec budget).'); + end + if isempty(p.Results.Depth) + if any(strcmp(lfa_type, {'block', 'block_qfa'})), d = 20; else, d = 200; end + else + d = double(p.Results.Depth); + end + + % --- Symmetry validation (default on; tolerance-gated). The MEX treats + % the upper triangle as authoritative; an asymmetric A has its lower triangle + % discarded with no other signal. Checked here, not in the MEX, so a + % direct fun_nystrom_pp_mex(...) call remains a raw low-level entry point + % (matches the header comment's "low-level" framing) while every call + % through this wrapper is protected by default. --- + % Vector mode has no lower triangle to discard, and A - A.' would broadcast + % an n x 1 vector into an n x n array (20 GB at n = 50,000), so the check is + % bypassed outright there rather than through SkipSymCheck. + if ~skip_sym_check && ~vector_mode + if issparse(A) % on the nonzeros: A(:) of a large sparse A would be an n^2-row vector + A_diff = A - A.'; + asym = max([0; abs(nonzeros(A_diff))]); + ascale = max([0; abs(nonzeros(A))]); + clear A_diff + else + % By column blocks of about 32 MB: A - A.' and abs(A(:)) would each be a full n x n temporary (two extra + % copies of A at peak). The maxima are exact, so the result equals the whole-matrix form bit for bit. + nb = max(1, floor(4e6 / max(n, 1))); + asym = zeros(1, 1, cls); ascale = zeros(1, 1, cls); + for j0 = 1:nb:n + J = j0:min(n, j0 + nb - 1); + Aj = A(:, J); + D = abs(Aj - A(J, :).'); % (:) not 'all': max(...,'all') needs R2018b + asym = max(asym, max(D(:))); + ascale = max(ascale, max(abs(Aj(:)))); + end + clear Aj D + end + if ascale == 0, ascale = 1; end + sym_tol = sqrt(eps(cls)); + if asym > sym_tol * ascale + error('randlapack:fun_nystrom_pp:NotSymmetric', ... + ['A must be symmetric (checked within tolerance): ' ... + 'max(abs(A-A'')) = %.3e exceeds tol*max(abs(A)) = %.3e ' ... + '(tol = sqrt(eps(''%s'')) = %.1e). The MEX reads only the ' ... + 'upper triangle as the symmetric matrix; a ' ... + 'non-symmetric A silently has its lower triangle discarded. ' ... + 'If A is known to be symmetric (e.g. by construction) and ' ... + 'this O(n^2) check is unwanted overhead at large n, pass ' ... + '''SkipSymCheck'', true.'], ... + asym, sym_tol * ascale, cls, sym_tol); + end + end + + % A is templated on by the MEX; keep it in the working precision. + A = cast(A, cls); + % Sparse A goes to the MEX as zero-based compressed-column arrays: find() returns the nonzeros column by column with + % rows ascending (CSC order) in O(nnz). The MEX cannot read a MATLAB sparse array's own buffers, and walking it element + % by element there costs O(n) per nonzero. + csc_ms = 0; + if issparse(A) + csc_timer = tic; + [ri, ci, vv] = find(A); + cnt = accumarray(ci(:), 1, [size(A, 2), 1]); + A = struct('sparse_n', int64(size(A, 1)), 'colptr', int64([0; cumsum(cnt)]), ... + 'rowidx', int64(ri(:) - 1), 'vals', double(vv(:))); + clear ri ci vv cnt + csc_ms = 1000 * toc(csc_timer); + end + + mex_args = {A, a1, a2, func, q, pl, lfa_type, d, ... + sketch, vec_nnz, sk_seed, reorth, adaptive, adapt_tol, ... + adapt_dl, adapt_mn, budget, auto_eps, radau_ret, ... + auto_dcap, auto_pfr, adapt_mvc, adapt_kc, adapt_sc, ... + cap_rank_fraction, quadrature_fraction}; + % Preserve the legacy positional call when this optional route is disabled. + if p.Results.FirstRowQL || spend_cap, mex_args{end+1} = double(logical(p.Results.FirstRowQL)); end + if spend_cap, mex_args{end+1} = 1; end + if nargout <= 1 + est = fun_nystrom_pp_mex(mex_args{:}); + elseif nargout == 2 + [est, t1] = fun_nystrom_pp_mex(mex_args{:}); + elseif nargout == 3 + [est, t1, t2] = fun_nystrom_pp_mex(mex_args{:}); + else + [est, t1, t2, times] = fun_nystrom_pp_mex(mex_args{:}); + % The compressed-column conversion is part of getting a sparse A into the MEX. + times.marshal_in_ms = times.marshal_in_ms + csc_ms; + end +end diff --git a/matlab/+randlapack/load_dense.m b/matlab/+randlapack/load_dense.m new file mode 100644 index 0000000..59a2371 --- /dev/null +++ b/matlab/+randlapack/load_dense.m @@ -0,0 +1,33 @@ +function [A, stats] = load_dense(file, n_rows, n_cols, varargin) +%LOAD_DENSE Read a float64 matrix from a raw file with many threads, so its pages spread over the NUMA nodes (Linux). +% [A, stats] = randlapack.load_dense(file, n_rows, n_cols, 'Offset', bytes, 'Divisor', d, 'Threads', t, ... +% 'Placement', 'interleave' | 'first_touch') +% +% Reads n_rows*n_cols float64 values stored contiguously at byte Offset (default 0) of file, in MATLAB's column-major +% order (as fread(fid, [n_rows n_cols], 'double') would), into a new array, and divides every value by Divisor +% (default 1: no division) with the same correctly rounded IEEE division as MATLAB's A/d. Each of Threads threads +% (default 0: OMP_NUM_THREADS, else all cores) reads its own contiguous span with pread into storage nobody has +% touched yet, and only ONE copy of the matrix exists at peak (fread into zeros(n) followed by A = A + 0 holds two). +% Placement 'interleave' (default) gives the empty array the MPOL_INTERLEAVE policy first, so its pages fault in +% round robin over the NUMA nodes (what numactl --interleave=all does for the whole process); 'first_touch' +% leaves each page on the node of the thread that reads it. A .npy file of a C-ordered symmetric matrix +% can be read this way too: its row-major data is the transpose, which equals the matrix when it is symmetric. +% +% stats: threads, ms, rss_before_mib, rss_after_alloc_mib (about rss_before: the pages are untouched until read), +% rss_after_fill_mib, interleave_set (false on a one-node machine, where there is nothing to interleave), +% interleave_errno, nodes (allowed node ids). +% +% See also randlapack.numa_spread, randlapack.numa_query. + p = inputParser; + addParameter(p, 'Offset', 0, @(x) isnumeric(x) && isscalar(x) && x >= 0 && x == floor(x)); + addParameter(p, 'Divisor', 1, @(x) isnumeric(x) && isscalar(x) && isfinite(x) && x ~= 0); + addParameter(p, 'Threads', 0, @(x) isnumeric(x) && isscalar(x) && x >= 0 && x == floor(x)); + addParameter(p, 'Placement', 'interleave', @(x) any(strcmp(x, {'interleave', 'first_touch'}))); + parse(p, varargin{:}); + o = p.Results; + if exist('numa_place_mex', 'file') ~= 3 + error('randlapack:load_dense:unavailable', 'numa_place_mex is not built for this platform.'); + end + [A, stats] = numa_place_mex('load', char(file), double(o.Offset), double(n_rows), double(n_cols), ... + double(o.Divisor), double(o.Threads), double(strcmp(o.Placement, 'interleave'))); +end diff --git a/matlab/+randlapack/numa_query.m b/matlab/+randlapack/numa_query.m new file mode 100644 index 0000000..a11df8d --- /dev/null +++ b/matlab/+randlapack/numa_query.m @@ -0,0 +1,13 @@ +function h = numa_query(A, max_pages) +%NUMA_QUERY Memory node of the pages of a dense array, sampled evenly (Linux; a query, nothing moves). +% h = randlapack.numa_query(A) up to 4,096 sampled pages +% h = randlapack.numa_query(A, max_pages) +% h.counts(k+1) = sampled pages on node k; h.absent = pages not faulted in; h.errors = other failures; +% h.sampled = pages sampled (one per stratum); h.page_offset, h.align64 = A's first byte modulo 4096 and 64. +% See also randlapack.numa_spread. + if nargin < 2, max_pages = 4096; end + if exist('numa_place_mex', 'file') ~= 3 + h = struct('counts', [], 'absent', NaN, 'sampled', 0, 'errors', NaN, 'page_offset', NaN, 'align64', NaN); return; + end + h = numa_place_mex('query', A, double(max_pages)); +end diff --git a/matlab/+randlapack/numa_spread.m b/matlab/+randlapack/numa_spread.m new file mode 100644 index 0000000..6ce5b11 --- /dev/null +++ b/matlab/+randlapack/numa_spread.m @@ -0,0 +1,32 @@ +function stats = numa_spread(A) +%NUMA_SPREAD Interleave the memory pages of a large dense array over the NUMA nodes, in place (Linux). +% stats = randlapack.numa_spread(A) +% +% MATLAB fills a matrix read from disk (fread, load) or built by assignment with one thread, so all of its pages sit on +% one memory node and every multithreaded product with it is limited by that node's bandwidth (2.4x slower products +% at n = 50,000 on a 4-node Xeon 6248R bigmem host, 2026-09-24). numa_spread moves A's pages in place so they interleave over +% the nodes this process may use, round robin by 2 MiB granule (move_pages with an explicit destination per page; +% mbind's MPOL_INTERLEAVE + MPOL_MF_MOVE would leave pages that already sit on an allowed node unmoved). A's values +% are not touched and no copy of A is made, unlike re-touching it with A = A + 0. A and every variable sharing its +% data keep the new placement. +% +% stats fields: supported, nodes (allowed node ids), bytes, threads (that ran), ms, pages (pages of A's span submitted +% for movement; 0 when nothing needs to move, i.e. one allowed node), retry_calls (move_pages retries), on_target (pages +% on their destination node at the end), not_moved (pages still elsewhere or refused, after one retry), unresolved +% (pages the kernel never reported on: a failed call or an early stop), failed_batches, absent (pages never +% touched), errno (the first error seen), and before/after page histograms (counts per node over up to 4,096 pages +% sampled one per stratum; absent = never touched, errors = other failures). Placement is complete when +% on_target + absent == pages. On a one-node machine nothing moves. Where the helper MEX is absent (not Linux, or +% not built) stats.supported is false and A is unchanged. +% +% See also randlapack.numa_query, randlapack.load_dense. + if ~(isnumeric(A) && ~issparse(A) && isreal(A) && (isa(A,'double') || isa(A,'single'))) + error('randlapack:numa_spread:class', 'A must be a full real single or double array.'); + end + if exist('numa_place_mex', 'file') ~= 3 + stats = struct('supported', false, 'nodes', [], 'bytes', 0, 'threads', 0, 'ms', 0, 'pages', 0, 'retry_calls', 0, 'on_target', 0, ... + 'failed_batches', 0, 'not_moved', 0, 'unresolved', 0, 'absent', 0, 'errno', 0, 'before', [], 'after', []); + return; + end + stats = numa_place_mex('spread', A); +end diff --git a/matlab/+randlapack/trace_estimate.m b/matlab/+randlapack/trace_estimate.m new file mode 100644 index 0000000..dbe62c6 --- /dev/null +++ b/matlab/+randlapack/trace_estimate.m @@ -0,0 +1,109 @@ +function [estimate, info] = trace_estimate(A, varargin) +%TRACE_ESTIMATE Four funNystrom++ methods with explicit work limits. +% [estimate, info] = randlapack.trace_estimate(A, 'Method', 'scalar_auto', ... +% 'TargetRelError', 1e-3, 'MaxMatvecs', 4000) +% +% Method: scalar, block, scalar_auto, block_auto. Fixed methods require Depth. +% All methods use sphere probes and count block products by their column count. +% Auto methods discover depth and choose rank/probe count from the target. +% If that allocation exceeds the cap, at most half of the post-pilot work +% goes to rank; the rest funds residual probes. The expert fun_nystrom_pp +% interface retains the old policy and experimental splits for reproduction. +% The target guides allocation; it is NOT a total-error guarantee. Returned +% quadrature_certified reports only Lanczos stopping, not sampling accuracy. +% oracle_mean_depth is actual Phase-2 work per final probe. depth_limit is +% the selected limit; legacy d_used reports this limit for automatic methods. +% The fixed methods split the budget equally between rank and quadrature, +% subject to rank <= n/2 and the block Krylov dimension <= n. +% +% A is a dense symmetric matrix or a vector representing a diagonal operator. +% Func supports sqrt, log (= log(1+x)), effdim. Lambda sets the denominator +% shift for effdim x/(x+Lambda), default 10 (same as the low-level interface). +% VecNnz controls SASO nonzeros per row (default 8, clipped to rank). +% VecNnz=0 uses the kernel's logarithmic density rule. Sketch type is SASO. +% Entropy is available with every method. See fun_nystrom_pp. +% SpendCap (auto methods only, default false) makes the method spend the whole +% cap: rank and probe count use it, every probe runs to the pilot's depth, and +% quadrature_certified then reports the pilot's certificate, since nothing +% later is checked. +% FirstRowQL optionally accelerates scalar projected solves. Single precision +% and nonconverged QL solves retain the established evaluator. Its use and any +% fallback, including the pilot, are recorded in info. + p = inputParser; + addParameter(p, 'Method', 'scalar_auto'); + addParameter(p, 'Func', 'sqrt'); + addParameter(p, 'Lambda', 10); + addParameter(p, 'VecNnz', 8); + addParameter(p, 'TargetRelError', 1e-3); + addParameter(p, 'MaxMatvecs', []); + addParameter(p, 'Depth', []); + addParameter(p, 'Seed', 42); + addParameter(p, 'SkipSymCheck', false); + addParameter(p, 'SpendCap', false, @(x) isscalar(x) && (islogical(x) || isnumeric(x)) && any(x==[0 1])); + addParameter(p, 'FirstRowQL', false, @(x) isscalar(x) && (islogical(x) || isnumeric(x)) && any(x==[0 1])); + parse(p, varargin{:}); + o = p.Results; + method = validatestring(o.Method, {'scalar','block','scalar_auto','block_auto'}); + func = validatestring(o.Func, {'sqrt','log','effdim','identity','entropy','poly','square'}); + validateattributes(o.MaxMatvecs, {'numeric'}, {'scalar','finite','integer','positive'}); + validateattributes(o.TargetRelError, {'numeric'}, {'scalar','finite','>',0,'<',1}); + validateattributes(o.Lambda, {'numeric'}, {'scalar','real','finite','positive'}); + validateattributes(o.VecNnz, {'numeric'}, {'scalar','real','finite','integer','nonnegative'}); + B = o.MaxMatvecs; + n = size(A, 1); + if isvector(A), n = numel(A); end + args = {'Func', func, 'VecNnz', o.VecNnz, 'SketchSeed', o.Seed, 'SkipSymCheck', o.SkipSymCheck}; + if o.FirstRowQL,args=[args,{'FirstRowQL',true}];end + if ismember(func,{'effdim','poly'}), args=[args,{'PolyLambda',o.Lambda}]; end + if endsWith(method, '_auto') + % The Gauss / Radau-at-0 bracket needs the derivatives of f from order two up to + % alternate in sign, not operator monotonicity. -x log x has the same pattern as + % log(1+x), so entropy qualifies; x^2 has a second derivative of +2 and does not. + validatestring(func, {'sqrt','log','effdim','identity','entropy'}); + [estimate, ~, ~, t] = randlapack.fun_nystrom_pp(A, 1, 1, args{:}, ... + 'LFAType', method, 'Budget', B, 'AutoEps', o.TargetRelError, ... + 'CapRankFraction', 0.5, 'SpendCap', logical(o.SpendCap)); + k = t.auto_k; s = t.auto_s; pilot = t.probe_mv; + if o.SpendCap, cert = logical(t.probe_converged); else, cert = logical(t.phase2_certified); end + pilot_ms = t.probe_ms; + else + if o.SpendCap + error('randlapack:trace_estimate:spendCap', 'SpendCap applies only to scalar_auto and block_auto.'); + end + validateattributes(o.Depth, {'numeric'}, {'scalar','finite','integer','positive'}); + d = min(n, o.Depth); + k = max(1, min(floor(n/2), floor(B/2))); + s = min(n-k, floor((B-k)/d)); + if strcmp(method, 'block'), s = min(s, floor(n/d)); end + if s < 1 + error('randlapack:trace_estimate:infeasibleBudget', ... + 'Budget %d cannot fund rank %d and one depth-%d probe.', B, k, d); + end + [estimate, ~, ~, t] = randlapack.fun_nystrom_pp(A, k, s, args{:}, ... + 'LFAType', [method '_qfa'], 'Q', 1, 'Depth', d, 'Adaptive', 0); + pilot = 0; cert = false; pilot_ms = 0; + end + info = t; + info.method = method; + info.function_name = func; + info.rank = k; + info.samples = s; + info.probe_ms = pilot_ms; + info.matvecs = pilot + k + t.oracle_mv; + info.algorithm_ms = pilot_ms + t.phase1_ms + t.phase2_ms; + info.oracle_mean_depth = NaN; + if s > 0, info.oracle_mean_depth = t.oracle_mv / s; end + if endsWith(method, '_auto'), info.depth_limit = t.d_used; else, info.depth_limit = d; end + info.quadrature_certified = cert; + info.target_rel_error = o.TargetRelError; + info.max_matvecs = B; + info.total_error_certified = false; + info.cap_rank_fraction = 0.5; + info.spend_cap = logical(o.SpendCap); + info.sketch = 'saso'; + info.vec_nnz_requested = o.VecNnz; + info.vec_nnz_effective = min(k,o.VecNnz); + if o.VecNnz==0, info.vec_nnz_effective=min(k,max(4,ceil(log(k)))); end + assert(info.matvecs <= B, 'randlapack:trace_estimate:budgetExceeded', ... + 'Spent %d matvecs with a hard cap of %d.', info.matvecs, B); +end diff --git a/matlab/CMakeLists.txt b/matlab/CMakeLists.txt index a85f967..6381a8e 100644 --- a/matlab/CMakeLists.txt +++ b/matlab/CMakeLists.txt @@ -7,6 +7,23 @@ set(MEX_OUTPUT_DIR ${CMAKE_CURRENT_SOURCE_DIR}/+randlapack/private) file(MAKE_DIRECTORY ${MEX_OUTPUT_DIR}) +# Statically link libstdc++ into each MEX. MATLAB R2026a ships its own +# (older) libstdc++.so.6 under sys/os/glnxa64/ which lacks newer GLIBCXX +# symbols required by code built with GCC 13+. Statically linking avoids +# the runtime dependency entirely — the MEX file ships its own libstdc++ +# symbols and is loadable in any MATLAB without LD_PRELOAD gymnastics. +# Applied via target_link_options on each matlab_add_mex target below. +# +# These are GCC/Linux-only flags. AppleClang on macOS rejects -static-libgcc +# ("unsupported option"): there is no libgcc to statically link (macOS uses +# libc++, not libstdc++), and the GLIBCXX-version problem above does not exist +# there. Leave the list empty off GNU so macOS/other compilers link normally. +if(CMAKE_CXX_COMPILER_ID STREQUAL "GNU") + set(MEX_STATIC_LIBSTDCXX_FLAGS -static-libstdc++ -static-libgcc) +else() + set(MEX_STATIC_LIBSTDCXX_FLAGS "") +endif() + matlab_add_mex( NAME bqrrp_mex SRC src/bqrrp_mex.cc @@ -26,3 +43,74 @@ set_target_properties(bqrrp_mex PROPERTIES RUNTIME_OUTPUT_DIRECTORY ${MEX_OUTPUT_DIR} ) target_compile_features(bqrrp_mex PRIVATE cxx_std_20) +target_link_options(bqrrp_mex PRIVATE ${MEX_STATIC_LIBSTDCXX_FLAGS}) +# Windows: stage the BLAS backend's runtime DLLs next to the MEX. MATLAB loads +# a MEX with the MEX file's own directory on the DLL search path, so app-local +# staging (the idiomatic Windows layout RandLAPACK's installer uses for its own +# executables) makes the MEX loadable without PATH edits. The helper is a +# no-op off Windows and only exists in RandLAPACK installs that export it. +if(COMMAND randlapack_stage_runtime_dlls) + randlapack_stage_runtime_dlls(bqrrp_mex) +endif() + +matlab_add_mex( + NAME fun_nystrom_pp_mex + SRC src/fun_nystrom_pp_mex.cc + LINK_TO RandLAPACK::RandLAPACK + R2018a +) +target_include_directories(fun_nystrom_pp_mex PRIVATE + $ +) +set_target_properties(fun_nystrom_pp_mex PROPERTIES + LIBRARY_OUTPUT_DIRECTORY ${MEX_OUTPUT_DIR} + RUNTIME_OUTPUT_DIRECTORY ${MEX_OUTPUT_DIR} +) +target_compile_features(fun_nystrom_pp_mex PRIVATE cxx_std_20) +option(RANDLAPACK_BINDINGS_SYMMETRIC_SKETCH + "Use triangle-aware sketch products from a supporting RandBLAS installation" OFF) +if(RANDLAPACK_BINDINGS_SYMMETRIC_SKETCH) + target_compile_definitions(fun_nystrom_pp_mex PRIVATE RANDLAPACK_SYMMETRIC_SKETCH=1) +endif() +target_link_options(fun_nystrom_pp_mex PRIVATE ${MEX_STATIC_LIBSTDCXX_FLAGS}) +if(COMMAND randlapack_stage_runtime_dlls) + randlapack_stage_runtime_dlls(fun_nystrom_pp_mex) +endif() + +# NUMA page placement for large dense inputs (randlapack.numa_spread / numa_query / load_dense). Linux only; it needs +# no RandLAPACK, only the MATLAB Data API. Built WITHOUT -ffast-math: load_dense's division must be the IEEE one. +if(UNIX AND NOT APPLE) + matlab_add_mex( + NAME numa_place_mex + SRC src/numa_place_mex.cc + R2018a + ) + set_target_properties(numa_place_mex PROPERTIES + LIBRARY_OUTPUT_DIRECTORY ${MEX_OUTPUT_DIR} + RUNTIME_OUTPUT_DIRECTORY ${MEX_OUTPUT_DIR} + ) + target_compile_features(numa_place_mex PRIVATE cxx_std_20) + if(CMAKE_CXX_COMPILER_ID MATCHES "GNU|Clang") + target_compile_options(numa_place_mex PRIVATE -fno-fast-math -fno-reciprocal-math) + endif() + # std::thread needs pthread_create; glibc < 2.34 (RHEL 8) keeps it in libpthread. Plain signature, as FindMatlab + # uses the plain target_link_libraries form on MEX targets. + find_package(Threads REQUIRED) + target_link_libraries(numa_place_mex Threads::Threads) + target_link_options(numa_place_mex PRIVATE ${MEX_STATIC_LIBSTDCXX_FLAGS}) +endif() + +# Windows: MATLAB's process already hosts Intel's OpenMP runtime (libiomp5md), +# and a MEX compiled with /openmp:llvm injects MSVC's libomp140; two OpenMP +# runtimes in one process abort at load (OMP Error #15). Intel's runtime and +# LLVM's implement the same kmpc ABI (LLVM's libomp descends from Intel's), so +# the supported MEX configuration is: satisfy the OpenMP references from +# MATLAB's own import library (shipped in bin/win64 for exactly this purpose) +# and suppress the compiler's default runtime. Plain target_link_libraries +# form, matching FindMatlab's internal usage (keyword/plain forms cannot mix). +if(WIN32) + foreach(mex_tgt IN ITEMS bqrrp_mex fun_nystrom_pp_mex) + target_link_libraries(${mex_tgt} "${Matlab_ROOT_DIR}/bin/win64/libiomp5md.lib") + target_link_options(${mex_tgt} PRIVATE "/nodefaultlib:libomp140.x86_64.lib") + endforeach() +endif() diff --git a/matlab/benchmark/funnystrompp/CMakeLists.txt b/matlab/benchmark/funnystrompp/CMakeLists.txt new file mode 100644 index 0000000..fb083c6 --- /dev/null +++ b/matlab/benchmark/funnystrompp/CMakeLists.txt @@ -0,0 +1,39 @@ +# funNystrom++ benchmark — external baselines via FetchContent. +# +# The benchmark compares RandLAPACK's FunNystromPP against external reference +# implementations. Those are NOT vendored in this repository; they are fetched +# at CMake-configure time, pinned to a commit each for reproducibility. The +# sources are MATLAB (.m), so there is no build step — FetchContent only +# downloads them, and we write their locations into external_paths.m, which the +# benchmark scripts addpath(). +# +# Acknowledgements / licenses (see ../../../THIRD_PARTY_LICENSES.md): +# Epperly https://github.com/eepperly/Ethan-Epperly-Thesis MIT +# XTrace / XNysTrace / resphering trace estimators. +# Persson https://github.com/davpersson/funNystrom license pending +# funNystrom++ reference (nystrom.m, funnystrompp.m, block_lanczos.m). +# No license file at time of writing: used here with attribution, +# fetched (not redistributed); a license is being requested from the +# author. block_lanczos.m therein is Meyer's HutchPlusPlus. + +include(FetchContent) + +FetchContent_Declare(epperly_trace + GIT_REPOSITORY https://github.com/eepperly/Ethan-Epperly-Thesis.git + GIT_TAG 2d33ee40577caff18b16fed3454a147d08cbe16e) # main, 2026-06 +FetchContent_Declare(persson_funnystrom + GIT_REPOSITORY https://github.com/davpersson/funNystrom.git + GIT_TAG 0427bde1e3c238ce53cc184600482481490c2a48) # master, 2026-06 + +# Download only (no add_subdirectory: these are MATLAB sources, no CMake build). +FetchContent_MakeAvailable(epperly_trace persson_funnystrom) + +set(EPPERLY_DIR "${epperly_trace_SOURCE_DIR}") +set(PERSSON_DIR "${persson_funnystrom_SOURCE_DIR}") +configure_file(external_paths.m.in + "${CMAKE_CURRENT_SOURCE_DIR}/external_paths.m" @ONLY) + +message(STATUS "funNystrom++ benchmark baselines fetched (not vendored):") +message(STATUS " Epperly (MIT): ${EPPERLY_DIR}") +message(STATUS " Persson (license pending): ${PERSSON_DIR}") +message(STATUS " -> external_paths.m written for the benchmark scripts") diff --git a/matlab/benchmark/funnystrompp/external_paths.m.in b/matlab/benchmark/funnystrompp/external_paths.m.in new file mode 100644 index 0000000..f1b8a12 --- /dev/null +++ b/matlab/benchmark/funnystrompp/external_paths.m.in @@ -0,0 +1,22 @@ +function p = external_paths() +%EXTERNAL_PATHS Locations of the fetched (NOT vendored) benchmark baselines, +% and put them on the MATLAB path. GENERATED by CMake (FetchContent) from +% external_paths.m.in — do not edit by hand; re-run cmake to regenerate. +% +% Acknowledgements (see ../../../THIRD_PARTY_LICENSES.md): +% Epperly github.com/eepperly/Ethan-Epperly-Thesis (MIT) — XTrace / XNysTrace +% / resphering trace estimators. +% Persson github.com/davpersson/funNystrom (license pending; used with +% attribution) — nystrom.m, funnystrompp.m, block_lanczos.m +% (block_lanczos.m is Meyer's HutchPlusPlus, vendored by Persson). +% +% Fetched at the pinned commits recorded in CMakeLists.txt. + + p.epperly = '@EPPERLY_DIR@'; + p.persson = '@PERSSON_DIR@'; + + % genpath: the upstream repos nest the .m files in subdirectories + % (Epperly under code/..., Persson under Other/), so add recursively. + addpath(genpath(p.epperly)); + addpath(genpath(p.persson)); +end diff --git a/matlab/src/fun_nystrom_pp_mex.cc b/matlab/src/fun_nystrom_pp_mex.cc new file mode 100644 index 0000000..8e1d117 --- /dev/null +++ b/matlab/src/fun_nystrom_pp_mex.cc @@ -0,0 +1,1325 @@ +// MEX wrapper for RandLAPACK::FunNystromPP using MATLAB's C++ Data API +// (R2018a+). +// +// MATLAB-side signature (low-level; users normally call randlapack.fun_nystrom_pp.m): +// +// [est, t1, t2] = fun_nystrom_pp_mex(A, arg1, arg2, func, q, poly_lambda, lfa_type, d, ... +// sketch_type, vec_nnz, sketch_seed, reorth, ... +// adaptive, adaptive_tol, adaptive_delay, ... +// adaptive_min, budget, auto_eps, radau_return, ... +// auto_depth_cap, auto_probe_frac, adaptive_matvec_cap) +// +// where: +// A n x n matrix, single or double, column-major. Upper triangle +// defines the matrix; legacy builds mirror it for general sketch +// products, while symmetric-sketch builds read it directly. +// arg1 Phase-1 sketch SIZE: a SCALAR k (the rank). The sketch itself is +// a SASO sketch generated INSIDE RandLAPACK::NystromEVD from +// (sketch_seed, vec_nnz) — matching the paper's Algorithm 1 line 1. +// Passing an explicit matrix is an error (the dense-sketch mode was +// removed from the kernel). +// arg2 Phase-2 Hutchinson probes: pass a SCALAR s to sample n x s Gaussian probes +// internally (seed = sketch_seed + 1000), OR an explicit n x s matrix. +// Phase 2 is skipped when k == n (arg2 then unread). +// func 'sqrt' | 'log' | 'poly' | 'effdim' | 'square' | 'identity' +// poly: f(x) = x(x + poly_lambda) +// effdim: f(x) = x/(x + poly_lambda) ("effective dimension"; operator monotone) +// q subspace-iter count (>= 1). q = 1 is single-pass Nystrom. +// poly_lambda lambda used by func 'poly' and 'effdim' +// lfa_type 'exact' | 'scalar' | 'scalar_qfa' | 'block' | 'block_qfa' | 'auto' | 'adaptive' +// exact: build f(A) once via syevd, every f(A)*X is a GEMM (validation oracle) +// scalar: per-column scalar Lanczos-FA at depth d (= Lanczos quadrature per probe) +// scalar_qfa: per-column scalar Lanczos-QFA — the per-probe +// quadratic forms ωⱼᵀf(A)ωⱼ directly (diagonal of +// Ω₂ᵀf(A)Ω₂; basis-free, O(n·s) memory). With +// adaptive = 1 each probe stops at its own depth via +// the Gauss-Radau certificate and adaptive_tol is a +// CERTIFIED per-probe relative error; reorth / +// adaptive_delay / adaptive_min are IGNORED (the +// certificate replaces the windowed rule; the +// recurrence is intrinsically no-reorth). +// block: block Lanczos-FA at depth d (BLAS-3 accelerated variant) +// block_qfa: block Lanczos-QFA — forms Ω₂ᵀf(A)Ω₂ (s×s) directly, +// skipping the f(A)·Ω₂ mapback (driver.use_qfa = true). +// Input 13 (adaptive) selects the depth rule: +// 0 = fixed depth d; +// 1 = Radau-certified stop (stop_rule = Radau: the +// block Gauss / Gauss-Radau bracket closes +// within adaptive_tol — the same certified +// meaning as scalar_qfa's adaptive mode; no +// delay window); +// 2 = legacy window rule (stop_rule = Window; +// honors adaptive_delay / adaptive_min). +// Input 19 (radau_return) selects the value returned +// on a certified stop: 0 = block Gauss (default), +// 1 = (Gauss + Radau)/2 midpoint. +// auto: knob-free tier. The driver picks k, s, and the oracle +// depth from (budget, auto_eps) [inputs 17, 18], running +// the certified scalar QFA for both the depth probe and +// Phase 2; the positional k/s/d/reorth/adaptive inputs +// are IGNORED. Inputs 20/21 (auto_depth_cap, +// auto_probe_frac) tune the depth probe. An infeasible +// budget raises MATLAB error id +// 'randlapack:fun_nystrom_pp:infeasibleBudget'. +// adaptive: fully knob-free, EPS-TARGETED tier (no upfront +// matvec budget). The driver picks k, s, and the +// oracle depth t itself from a certified BLOCK +// Gauss-Radau depth probe run at rtol = auto_eps +// [input 18, reused as eps]; the positional +// k/s/d/reorth/adaptive inputs are IGNORED. Input 20 +// (auto_depth_cap) caps the probe depth, shared with +// the 'auto' tier. Input 22 (adaptive_matvec_cap, 0 = +// no cap) optionally bounds the total matvec spend; +// when it clamps the rank/probe split infeasibly +// small, or when the eps/n regime cannot fund the +// block-Krylov-coupled probe count (s*t <= n), the +// driver throws and the MEX raises the SAME error id +// as the 'auto' tier's infeasible budget: +// 'randlapack:fun_nystrom_pp:infeasibleBudget'. +// d Lanczos depth (used for lfa_type in {scalar, scalar_qfa, block, +// block_qfa}; the adaptive QFA modes treat it as a depth CAP; +// ignored, like the positional k/s, for 'auto'/'adaptive') +// sketch_type IGNORED (accepted for call-site compatibility; the Phase-1 +// sketch is always the kernel-internal SASO). A MATLAB warning is +// issued if anything other than 'saso' is passed explicitly. +// vec_nnz nonzeros per ROW of the SASO sketch (default 8; 0 = auto, +// resolved to ~log(k) inside NystromEVD) [optional, input 10] +// sketch_seed RNG seed for the Phase-1 sketch (default 42) [optional, input 11] +// ... optional inputs 12-18 (reorth, adaptive, adaptive_tol, +// adaptive_delay, adaptive_min, budget, auto_eps) are documented +// inline where they are read in run_fun_nystrom below. +// radau_return block_qfa certified return value: 0 = Gauss (default), +// 1 = midpoint [optional, input 19] +// auto_depth_cap 'auto' tier: fixed cap on the probe depth (0 = no fixed +// cap, the default) [optional, input 20] +// auto_probe_frac 'auto' tier: fraction of the matvec budget the depth +// probe may spend, in (0, 1) (default 0.125) [optional, input 21] +// adaptive_matvec_cap 'adaptive' tier only: optional total matvec cap (0 = +// no cap, the default). Ignored by every other lfa_type. +// [optional, input 22] +// ... optional inputs 23-27 (adaptive_k_const, adaptive_s_const, +// cap_rank_fraction, quadrature_fraction, first_row_ql) are +// documented inline where they are read. +// spend_cap 'scalar_auto'/'block_auto' only: 1 = spend the whole Budget +// (driver.adaptive_spend_cap: rank and probe count chosen to use +// the cap, Phase-2 probes at the probe's depth with no early +// stop); 0 = shipped behaviour (default) [optional, input 28] +// +// Outputs: +// est trace estimate t1 + t2 (scalar double) +// t1 Phase 1 contribution (scalar double) +// t2 Phase 2 contribution (scalar double; 0 when k == n) +// times (optional 4th output) struct of wall-clock instrumentation: +// marshal_in_ms reading A/Omega2 out of MATLAB arrays: validation only by default (the +// driver reads MATLAB's buffers in place), plus the private copy when +// RANDLAPACK_FNPP_COPY_INPUT=1 +// input_copied true when that private copy was made (see copy_input_requested) +// phase1_ms NystromEVD total (driver field) +// phase2_ms Phase 2 total (driver field) +// fafun_ms time inside the f(A)*X oracle (subset of phase2) +// assembly_ms phase2_ms - fafun_ms (trace assembly) +// specrec_ms spectral-recovery block inside NystromEVD +// nystrom_us 1x11 NystromEVD breakdown (microseconds). +// Only slots 0, 1, 2, 6, 10 (C++ 0-based; MATLAB +// indices 1, 2, 3, 7, 11) are populated: +// slot 0 = alloc, 1 = syrf (QR stabilization), +// 2 = matvec, 6 = the WHOLE shifted spectral- +// recovery block (Alg. 2 lines 3-8, not just an +// svd), 10 = total. The remaining slots are 0. +// lfa_us 1x6 Lanczos-oracle breakdown (microseconds): +// [matvec run_lanczos apply rest total reorth] +// for lfa_type 'scalar'/'scalar_qfa'/'block'/ +// 'block_qfa'/'auto' (for the QFA types, apply = +// certificate checks + final quadrature evals; +// zeros for 'exact'). Slot 6 (reorth) is the +// recurrence's reorthogonalization time; it is +// real for 'scalar'/'block'/'block_qfa' (the +// block QFA reuses the FA recurrence and pays +// reorth whenever reorth = 1) and 0 for +// 'scalar_qfa'/'auto' (basis-free recurrence, +// no reorth by design). +// d_used Actual depth for fixed/QFA modes; a selected +// depth limit for automatic modes. For +// 'block_qfa' + adaptive this is the online-chosen +// depth (<= d cap); for 'scalar_qfa' + adaptive it +// is the MAX per-probe certified depth; for 'auto' +// it is the probe-discovered depth cap t; for +// 'adaptive' it is the depth probe's certified (or +// reached) depth t (driver.adaptive_t); for 'exact' +// it is NaN (no Lanczos recurrence runs; the fixed +// d input is not a meaningful cost figure there); +// otherwise (scalar/block/block_qfa fixed-depth) it +// equals the fixed d. +// oracle_mv total A-matvecs the Phase-2 oracle actually +// spent: Σ per-probe depths for 'scalar_qfa' and +// 'auto'; s*d_used for 'block_qfa'; the block +// oracle's matvecs member for 'adaptive' +// (driver.adaptive_oracle_matvecs); 0 for +// 'exact'/'scalar'/'block' (their count is the +// analytic s*d_used). The benchmark should prefer +// this over d_used-based re-costing. +// auto_k 'auto'/'adaptive' only: chosen Nystrom rank +// (driver.auto_k / driver.adaptive_k; else 0) +// auto_s 'auto'/'adaptive' only: chosen probe count +// (driver.auto_s / driver.adaptive_s; else 0) +// probe_mv 'auto'/'adaptive' only: matvecs the depth probe +// actually spent (else 0). With the certified +// oracle the budget closes as an upper bound: +// probe_mv + q*auto_k + oracle_mv <= budget +// (Phase 1 costs q*auto_k matvecs, q=1 here; the +// 'adaptive' tier has no upfront budget, so this +// is an accounting identity, not a constraint). +// +// New fields below use the NaN convention: NaN when the path +// that produces them did not run (instead of a fake 0, which +// would be indistinguishable from a real measurement). +// tr_U 'block_qfa'/'adaptive' only: final block Gauss +// trace of Ω₂ᵀf(A)Ω₂ (upper side of the Radau +// bracket); for 'adaptive' this is +// driver.adaptive_tr_U; NaN otherwise. +// tr_L 'block_qfa'/'adaptive' only: final block +// Gauss-Radau trace (lower side; equals tr_U when +// no certificate ran); for 'adaptive' this is +// driver.adaptive_tr_L, NaN under spend_cap (no +// bracket evaluated in Phase 2); NaN otherwise. +// certified 'block_qfa'/'adaptive'/'scalar_qfa' (with +// Adaptive=1) only: 1 if the Radau bracket +// closed within (adaptive_tol / auto_eps), 0 if +// not; for 'adaptive' this is +// driver.adaptive_phase2_certified (NaN under +// spend_cap); for +// 'scalar_qfa' this is scalar_qfa.all_certified +// (every probe column certified before the +// depth cap); NaN otherwise (including +// 'scalar_qfa' with Adaptive=0, where no +// certificate was checked). +// probe_ms 'auto'/'adaptive' only: wall-clock of the depth +// probe (runs BEFORE phase1_ms's clock starts, so +// phase1_ms + phase2_ms excludes it; 'adaptive' +// reads driver.t_adaptive_probe_ms); NaN +// otherwise. +// probe_converged 'auto'/'adaptive' only: 1 if the depth probe +// certified before the probe cap ('adaptive' +// reads driver.adaptive_probe_certified); NaN +// otherwise. +// phase2_certified 'auto'/'adaptive' only: 1 if every Phase-2 +// oracle column/block certified at its depth cap +// (distinct from probe_converged, which names +// ONLY the probe; 'adaptive' reads +// driver.adaptive_phase2_certified); NaN +// otherwise. +// phase2_checked 'adaptive' only: 1 if Phase 2 evaluated its +// bracket, 0 under spend_cap (phase2_certified, +// certified and tr_L are then NaN because +// nothing was checked); NaN otherwise. +// +// A may be single or double; the computation runs in that precision and the +// scalar outputs are returned as double. Omega1/Omega2 must match the class +// of A (the randlapack.fun_nystrom_pp.m wrapper enforces this). + +#include "mex.hpp" +#include "mexAdapter.hpp" + +#include +#include "rl_blaspp.hh" +#include "rl_lapackpp.hh" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace linops = RandLAPACK::linops; +namespace testing = RandLAPACK::testing; + +namespace { +// Carrier for MATLAB errors raised inside the MEX. raise() throws this and +// operator() converts it to feval('error', id, msg) as the LAST action. The +// former design called feval directly from raise(): the engine's +// MATLABException then re-entered operator()'s catch chain, where the typed +// catch (const matlab::Exception&) fails to match it (this MEX statically +// links libstdc++ while MATLAB's process carries its own copy, and typed +// catches of foreign classes are unreliable across the two RTTI domains), so +// the exception fell into catch (const std::exception&) and EVERY error id +// degraded to the generic StdError. A TU-local type in an anonymous +// namespace has its typeinfo wholly inside this MEX, so its catch always +// matches, and the terminal feval exception escapes operator() with no +// handler left to clobber the id. +struct MexError { + std::string id; + std::string msg; +}; +} // namespace + +using matlab::mex::ArgumentList; +using matlab::data::Array; +using matlab::data::ArrayFactory; +using matlab::data::ArrayType; +using matlab::data::CharArray; +using matlab::data::SparseArray; +template using TypedArray = matlab::data::TypedArray; + + +class MexFunction : public matlab::mex::Function { +private: + std::shared_ptr matlabPtr; + ArrayFactory factory; + + // Raise a MATLAB error with a given identifier and message. Throws the + // TU-local MexError; operator() catches it and emits the actual MATLAB + // error (see the MexError comment for why the feval must happen there). + [[noreturn]] void raise(const std::string& id, const std::string& msg) { + throw MexError{id, msg}; + } + + // Terminal error emission: feval('error', id, msg) raises a MATLAB-side + // error whose engine exception propagates out of the MEX unhandled, + // which is exactly what delivers the id to the MATLAB caller. Called + // only from operator()'s catch handlers. + [[noreturn]] void emit_matlab_error(const std::string& id, const std::string& msg) { + matlabPtr->feval(u"error", 0, std::vector{ + factory.createCharArray(id), + factory.createCharArray(msg) + }); + std::terminate(); // feval(error) does not return; this is a backstop. + } + + std::string read_string(const Array& a, const std::string& field) { + if (a.getType() != ArrayType::CHAR) { + raise("randlapack:fun_nystrom_pp_mex:string", + field + " must be a char array (single quotes), not a MATLAB string"); + } + CharArray ca = a; + return ca.toAscii(); + } + + // Any real numeric scalar class is accepted (double, single, and the + // signed/unsigned integer classes): a direct MEX call with int32(60) or + // single(1e-2) is legitimate, and the Data API's TypedArray + // conversion throws InvalidArrayTypeException on anything but double. + double read_double(const Array& a, const std::string& field) { + if (a.getNumberOfElements() != 1) { + raise("randlapack:fun_nystrom_pp_mex:scalar", field + " must be a scalar"); + } + switch (a.getType()) { + case ArrayType::DOUBLE: { TypedArray v = a; return static_cast(v[0]); } + case ArrayType::SINGLE: { TypedArray v = a; return static_cast(v[0]); } + case ArrayType::INT8: { TypedArray v = a; return static_cast(v[0]); } + case ArrayType::INT16: { TypedArray v = a; return static_cast(v[0]); } + case ArrayType::INT32: { TypedArray v = a; return static_cast(v[0]); } + case ArrayType::INT64: { TypedArray v = a; return static_cast(v[0]); } + case ArrayType::UINT8: { TypedArray v = a; return static_cast(v[0]); } + case ArrayType::UINT16: { TypedArray v = a; return static_cast(v[0]); } + case ArrayType::UINT32: { TypedArray v = a; return static_cast(v[0]); } + case ArrayType::UINT64: { TypedArray v = a; return static_cast(v[0]); } + default: + raise("randlapack:fun_nystrom_pp_mex:scalar", + field + " must be a real numeric scalar"); + } + } + + int64_t read_int(const Array& a, const std::string& field) { + return static_cast(read_double(a, field)); + } + + // Copy a TypedArray column-major into a fresh heap buffer, owned by a + // std::unique_ptr so that every raise() (which throws through this + // frame) and any exception out of driver.call frees it — the raw + // new[]/delete[] form leaked n^2 + n*s elements on every error path after + // the marshal. MATLAB arrays are + // column-major and contiguous, so we bulk-copy the whole buffer in one + // shot; &*begin() is the contiguous storage of a full array, so a single + // memcpy replaces a per-element loop (a measured marshal bottleneck at + // large n). `new T[]` without () default-initializes — no wasted zeroing + // pass before the memcpy overwrites every entry (std::vector::resize + // value-initializes, i.e. touches the n^2 buffer twice). + // + // `const TypedArray` is essential: a non-const TypedArray selects the + // mutable begin(), which forces MATLAB's copy-on-write to UNSHARE (deep- + // copy) the whole n^2 array before the memcpy copies it again. The const + // overload reads the shared buffer in place, so the input A is copied once, + // not twice — this is the dominant piece of the flat marshal cost at n=3000. + template + std::unique_ptr copy_into(const Array& in, size_t n_elems) { + std::unique_ptr out(new T[n_elems]); + const TypedArray typed = in; + if (n_elems == 0) return out; + // Opt-in performance switch, off by default: RANDLAPACK_PERF_PARCOPY=1 splits the copy over + // threads. A fresh buffer's pages land on the memory node of the thread that first writes them, so a + // one-thread memcpy puts the whole matrix on one node and every later product streams it from there. + const char* pc = std::getenv("RANDLAPACK_PERF_PARCOPY"); + if (pc != nullptr && pc[0] == '1' && n_elems >= ((size_t)1 << 20)) { + // Never more threads than cores or than 2^20-element chunks, and never let a failed + // thread creation escape: a joinable std::thread destroyed by an exception would + // terminate MATLAB. On failure the started threads are joined and the plain copy runs. + const char* ot = std::getenv("OMP_NUM_THREADS"); + const size_t hw = (size_t)std::max(1u, std::thread::hardware_concurrency()); + size_t nt = ot ? (size_t)std::max(1, std::atoi(ot)) : hw; + nt = std::min({nt, hw, n_elems >> 20}); + const T* src = &*typed.cbegin(); T* dst = out.get(); + auto chunk = [=](size_t t) { + size_t lo = n_elems * t / nt, hi = n_elems * (t + 1) / nt; + std::memcpy(dst + lo, src + lo, (hi - lo) * sizeof(T)); + }; + std::vector pool; + try { + pool.reserve(nt); + for (size_t t = 1; t < nt; ++t) pool.emplace_back(chunk, t); + } catch (const std::exception&) { + for (auto& th : pool) th.join(); + std::memcpy(dst, src, n_elems * sizeof(T)); + return out; + } + chunk(0); + for (auto& th : pool) th.join(); + return out; + } + std::memcpy(out.get(), &*typed.cbegin(), n_elems * sizeof(T)); + return out; + } + + // Zero-copy input (the default): the address of MATLAB's own column-major buffer, read in place. The const + // TypedArray selects the const iterator, so MATLAB's copy-on-write never unshares the array (see copy_into). + // The temporary TypedArray shares the buffer with `in`, which is an element of the MEX's ArgumentList and + // keeps the data alive for the whole call, so the pointer stays valid after `typed` goes out of scope. + // Every consumer reads A through const T*: ExplicitSymLinOp::A_buff, DiagSymLinOp::lambda, + // make_exact_fa_oracle (which copies before its destructive syevd) and FunNystromPP::call's Omega2; there is + // no const_cast or pointer cast anywhere on those paths (RandLAPACK fc4423e, RandBLAS 86ed633), so a write + // into MATLAB's array would not compile. The read-only proof is in the 2026-10-07 memopt report. + template + static const T* matlab_data(const Array& in) { + const TypedArray typed = in; + if (typed.getNumberOfElements() == 0) return nullptr; + return &*typed.cbegin(); + } + + // Input marshalling: RANDLAPACK_FNPP_COPY_INPUT=1 restores the private copy of A, Omega2 and the sparse CSC arrays + // (the code the paper ran; RANDLAPACK_PERF_PARCOPY=1 still threads that copy). Off by default: the MEX reads + // MATLAB's buffers in place, so a dense n x n problem holds one copy of A instead of two. Builds without + // RANDLAPACK_SYMMETRIC_SKETCH mirror the upper triangle into the buffer and therefore always copy. + static bool copy_input_requested() { +#ifndef RANDLAPACK_SYMMETRIC_SKETCH + return true; +#else + const char* v = std::getenv("RANDLAPACK_FNPP_COPY_INPUT"); + return v != nullptr && v[0] == '1'; +#endif + } + + // RANDLAPACK_PERF_PARCOPY=1 threads the private copy; without the copy it has nothing to do. Scripts from before + // 2026-10-07 set it to get a spread copy, so say once per MATLAB session (per MEX load) that it no longer applies: + // placement is now the caller's job (harness utils/numa_place.m, randlapack.numa_spread / load_dense). + void warn_parcopy_unused_once() { + static bool warned = false; + const char* pc = std::getenv("RANDLAPACK_PERF_PARCOPY"); + if (warned || pc == nullptr || pc[0] != '1') return; + warned = true; + matlabPtr->feval(u"warning", 0, std::vector{ + factory.createCharArray("randlapack:fun_nystrom_pp_mex:parcopyUnused"), + factory.createCharArray("RANDLAPACK_PERF_PARCOPY=1 has no effect: the MEX reads MATLAB's input arrays in place " + "(no copy to thread). For a dense A, spread its pages with randlapack.numa_spread " + "or load_dense; RANDLAPACK_FNPP_COPY_INPUT=1 restores the threaded private copy.")}); + } + + // Sparse A reaches the MEX from the wrapper as a 1x1 struct of zero-based compressed-column arrays: sparse_n (int64 + // scalar), colptr (int64, n+1), rowidx (int64, nnz), vals (double, nnz), extracted in MATLAB with find() in O(nnz). + // A raw MATLAB sparse array is refused: the Data API hides its buffers, and walking its nonzeros with getIndex costs + // O(n) per nonzero (the copy grew 4x per doubling of n: 14 s at n = 40,000, 55 s at 80,000; about an hour at 500k). + static bool is_sparse_struct(const Array& a) { + if (a.getType() != ArrayType::STRUCT || a.getNumberOfElements() != 1) return false; + const matlab::data::StructArray S = a; + for (const auto& f : S.getFieldNames()) + if (std::string(f) == "sparse_n") return true; + return false; + } + + Array sparse_field(const Array& in, const char* name, ArrayType want) { + const matlab::data::StructArray S = in; + bool has = false; + for (const auto& f : S.getFieldNames()) + if (std::string(f) == name) has = true; + if (!has) + raise("randlapack:fun_nystrom_pp_mex:A_sparse_field", + std::string("sparse A struct lacks field '") + name + "' (build it with randlapack.fun_nystrom_pp)"); + Array a = S[0][name]; + if (a.getType() != want) + raise("randlapack:fun_nystrom_pp_mex:A_sparse_field", + std::string("sparse A field '") + name + "' has the wrong class (int64 for sparse_n, colptr, rowidx; double for vals)"); + return a; + } + + int64_t sparse_n_of(const Array& in) { + const TypedArray nA = sparse_field(in, "sparse_n", ArrayType::INT64); + if (nA.getNumberOfElements() != 1 || nA[0] < 1) + raise("randlapack:fun_nystrom_pp_mex:A_sparse_n", "sparse A: sparse_n must be one positive int64"); + return nA[0]; + } + + // Validates the struct's arrays, then either copies them (copy == true) or points at MATLAB's own buffers: colptr + // starts at 0, never decreases and ends at nnz; every row index is in [0, n) and rows strictly increase within a + // column (the order element access and the sparse kernels rely on); every value is finite. The checks read the + // MATLAB arrays in place, so a refused matrix costs no copy. The zero-copy pointers are const, but the CSCMatrix + // built from them (a non-owning view, own_memory = false) stores plain T* / int64_t* members, so const SpMat& does + // NOT protect them: the sparse path is read-only by inspection, not by type. On this path RandBLAS's left_spmm / + // right_spmm, its CSC kernels and its sparse-times-sparse sketch only read the arrays; the MKL backend wraps them + // in mkl_sparse_?_create_csc handles for mkl_sparse_?_mm / _spmmd, which leave them unchanged; RandBLAS never calls + // mkl_sparse_order or reindex_inplace (RandBLAS 86ed633). Through randlapack.fun_nystrom_pp the arrays are fresh + // find() temporaries in any case, so a caller's sparse A could not be reached. + void marshal_csc(const Array& in, int64_t n, bool copy, std::unique_ptr& colptr_own, + std::unique_ptr& rowidx_own, std::unique_ptr& vals_own, + const int64_t*& colptr, const int64_t*& rowidx, const double*& vals, int64_t& nnz) { + const TypedArray cp = sparse_field(in, "colptr", ArrayType::INT64); + const TypedArray ri = sparse_field(in, "rowidx", ArrayType::INT64); + const TypedArray vv = sparse_field(in, "vals", ArrayType::DOUBLE); + nnz = static_cast(ri.getNumberOfElements()); + if (static_cast(cp.getNumberOfElements()) != n + 1) + raise("randlapack:fun_nystrom_pp_mex:A_sparse_colptr", "sparse A: colptr must have n+1 entries"); + if (static_cast(vv.getNumberOfElements()) != nnz) + raise("randlapack:fun_nystrom_pp_mex:A_sparse_nnz", "sparse A: rowidx and vals must have the same length"); + colptr = matlab_data(cp); + rowidx = matlab_data(ri); + vals = matlab_data(vv); + validate_csc(n, nnz, colptr, rowidx, vals); + if (copy) { + colptr_own.reset(new int64_t[n + 1]); + rowidx_own.reset(new int64_t[std::max(nnz, 1)]); + vals_own.reset(new double[std::max(nnz, 1)]); + std::memcpy(colptr_own.get(), colptr, static_cast(n + 1) * sizeof(int64_t)); + if (nnz > 0) { + std::memcpy(rowidx_own.get(), rowidx, static_cast(nnz) * sizeof(int64_t)); + std::memcpy(vals_own.get(), vals, static_cast(nnz) * sizeof(double)); + } + colptr = colptr_own.get(); rowidx = rowidx_own.get(); vals = vals_own.get(); + } else if (nnz == 0) { + // An empty MATLAB array has no buffer; the kernels never dereference these when nnz == 0, but keep them + // non-null, as the copy path's one-element allocations are. + rowidx_own.reset(new int64_t[1]); vals_own.reset(new double[1]); + rowidx = rowidx_own.get(); vals = vals_own.get(); + } + } + + void validate_csc(int64_t n, int64_t nnz, const int64_t* colptr, const int64_t* rowidx, const double* vals) { + if (colptr[0] != 0 || colptr[n] != nnz) + raise("randlapack:fun_nystrom_pp_mex:A_sparse_colptr", "sparse A: colptr must start at 0 and end at nnz"); + for (int64_t j = 0; j < n; ++j) // all of colptr first: a decrease would merge columns in the row scan below + if (colptr[j + 1] < colptr[j]) + raise("randlapack:fun_nystrom_pp_mex:A_sparse_colptr", + "sparse A: colptr decreases at column " + std::to_string(j + 1)); + for (int64_t j = 0; j < n; ++j) { + for (int64_t p = colptr[j]; p < colptr[j + 1]; ++p) { + if (rowidx[p] < 0 || rowidx[p] >= n) + raise("randlapack:fun_nystrom_pp_mex:A_sparse_row", + "sparse A: row index out of range in column " + std::to_string(j + 1)); + if (p > colptr[j] && rowidx[p] <= rowidx[p - 1]) + raise("randlapack:fun_nystrom_pp_mex:A_sparse_order", + "sparse A has unsorted or repeated row indices in column " + std::to_string(j + 1)); + if (!std::isfinite(vals[p])) + raise("randlapack:fun_nystrom_pp_mex:A_sparse_value", + "sparse A has a non-finite value in column " + std::to_string(j + 1)); + } + } + } + + // Templated worker: reads the T-typed matrices, builds the f(A)*X oracle, + // and drives FunNystromPP. T is double or float (dispatched in operator()). + template + void run_fun_nystrom(ArgumentList& outputs, ArgumentList& inputs) { + // --- input prep (timed). Seed-driven API: inputs 2 and 3 are EITHER a + // scalar (the sketch SIZE, sampled internally here from sketch_seed) OR + // an explicit dense sketch matrix (escape hatch for validation/repro). + // The driver always receives plain buffers, so it is untouched. --- + const auto t_marshal_start = std::chrono::steady_clock::now(); + using RNG = r123::Philox4x32; + + // A: n x n, or an n-vector holding diag(A) (RAII-owned; freed on every + // exit path, including raise() and exceptions out of driver.call). + // Vector mode copies n values, not n^2, and skips the triangle mirror + // below; the shape was validated in operator(). + const Array& A_in = inputs[0]; + const auto A_dims = A_in.getDimensions(); + // Sparse mode (sparse double A, both triangles stored, passed as the wrapper's CSC struct): a copy of the nonzeros, + // no n^2 buffer; A_buf stays null. + const bool sparse_mode = is_sparse_struct(A_in); + const bool vector_mode = !sparse_mode && (A_dims.size() == 2) && + ((A_dims[0] == 1 && A_dims[1] > 1) || (A_dims[1] == 1 && A_dims[0] > 1)); + const int64_t n = sparse_mode ? sparse_n_of(A_in) + : vector_mode ? static_cast(A_in.getNumberOfElements()) + : static_cast(A_dims[0]); + // Zero-copy by default (A_buf points into MATLAB's array); RANDLAPACK_FNPP_COPY_INPUT=1 or a legacy build + // copies into A_own as before. + const bool copy_input = copy_input_requested(); + if (!copy_input) warn_parcopy_unused_once(); + std::unique_ptr A_own; + std::unique_ptr sp_colptr_own, sp_rowidx_own; + std::unique_ptr sp_vals_own; + const int64_t* sp_colptr = nullptr; + const int64_t* sp_rowidx = nullptr; + const double* sp_vals = nullptr; + int64_t sp_nnz = 0; + const T* A_buf = nullptr; + if (sparse_mode) { + marshal_csc(A_in, n, copy_input, sp_colptr_own, sp_rowidx_own, sp_vals_own, + sp_colptr, sp_rowidx, sp_vals, sp_nnz); + } else if (copy_input) { + A_own = copy_into(A_in, vector_mode ? static_cast(n) : static_cast(n) * n); + A_buf = A_own.get(); + } else { + A_buf = matlab_data(A_in); + } + + // scalar params (read before sketch handling: internal sampling needs them) + const std::string func = read_string(inputs[3], "func"); + const int64_t q = read_int (inputs[4], "q"); + const T poly_lambda = static_cast(read_double(inputs[5], "poly_lambda")); + const std::string requested_lfa = read_string(inputs[6], "lfa_type"); + const bool target_auto = (requested_lfa == "scalar_auto" || requested_lfa == "block_auto"); + const std::string lfa_type = target_auto ? "adaptive" : requested_lfa; + // The Gauss / Radau-at-0 bracket needs the derivatives of f from order two up to + // alternate in sign, not f to be operator monotone. -x log x satisfies that with the + // same pattern as log(1+x): its nth derivative is (-1)^(n-1) (n-2)! / x^(n-1), so even + // orders are negative and odd orders positive. poly and square do not: x^2 has a second + // derivative of +2, the opposite sign, and the bracket does not hold for them. + if (target_auto && func != "sqrt" && func != "log" && func != "effdim" + && func != "identity" && func != "entropy") + raise("randlapack:fun_nystrom_pp:unsupportedAutoFunction", + "scalar_auto/block_auto support sqrt, log, effdim, identity and entropy; " + "poly and square have no Radau certificate"); + const int64_t d = read_int (inputs[7], "d"); + const std::string sketch_type = (inputs.size() >= 9) ? read_string(inputs[8], "sketch_type") : "saso"; + const int64_t vec_nnz_in = (inputs.size() >= 10) ? read_int(inputs[9], "vec_nnz") : 8; + // sketch_type 'gaussian' selects NystromEVD's dense Gaussian sketch (vec_nnz = -1). + const int64_t vec_nnz = (sketch_type == "gaussian") ? -1 : vec_nnz_in; + const int64_t sketch_seed = (inputs.size() >= 11) ? read_int(inputs[10], "sketch_seed") : 42; + // Optional input 12: Lanczos reorthogonalization flag (default 1 = full). + const int64_t reorth_flag = (inputs.size() >= 12) ? read_int(inputs[11], "reorth") : 1; + // Adaptive Lanczos-QFA depth (block_qfa only): choose the Lanczos depth + // online from the qfa certificate instead of the fixed d (which becomes + // the cap). adaptive_tol is the relative-change tolerance on tr(M_k). + const int64_t adaptive_fl = (inputs.size() >= 13) ? read_int(inputs[12], "adaptive") : 0; + const T adaptive_tl = (inputs.size() >= 14) + ? static_cast(read_double(inputs[13], "adaptive_tol")) : (T)1e-2; + // Certificate window (block_qfa adaptive only). The first convergence + // test is at depth adaptive_min + adaptive_delay, so these set the floor + // on d_used. 0 = keep the library default (see BlockLanczosQFA). + const int64_t adaptive_dl = (inputs.size() >= 15) ? read_int(inputs[14], "adaptive_delay") : 0; + const int64_t adaptive_mn = (inputs.size() >= 16) ? read_int(inputs[15], "adaptive_min") : 0; + // Knob-free tier (lfa_type == "auto"): total A-matvec budget and target + // accuracy; the driver picks k, s, and the oracle depth itself. The + // positional k/s/d/reorth/adaptive inputs are ignored in this mode. + const int64_t auto_budget = (inputs.size() >= 17) ? read_int(inputs[16], "budget") : 0; + const double auto_eps = (inputs.size() >= 18) ? read_double(inputs[17], "auto_eps") : 1e-3; + // Certified-return selector (block_qfa adaptive = 1 only): 0 = block + // Gauss, 1 = (Gauss + Radau)/2 midpoint. + const int64_t radau_ret = (inputs.size() >= 19) ? read_int(inputs[18], "radau_return") : 0; + // Auto-tier probe knobs: a fixed cap on the probe depth (0 = no fixed + // cap) and the fraction of the budget the probe may spend, in (0, 1). + const int64_t auto_dcap = (inputs.size() >= 20) ? read_int(inputs[19], "auto_depth_cap") : 0; + const double auto_pfrac = (inputs.size() >= 21) ? read_double(inputs[20], "auto_probe_frac") : 0.125; + // Eps-targeted tier (lfa_type == "adaptive") only: optional total + // matvec cap, 0 (default) => nullopt (no cap) at the driver call site. + const int64_t adaptive_mvcap = target_auto ? auto_budget + : ((inputs.size() >= 22) ? read_int(inputs[21], "adaptive_matvec_cap") : 0); + if (target_auto && auto_budget < 1) + raise("randlapack:fun_nystrom_pp:Budget", "scalar_auto/block_auto require a positive Budget"); + // Eps-targeted tier only: the leading constants of the paper's split + // rule k = c_k sqrt(t)/eps, s = max(1, c_s/(sqrt(t) eps)). Both default + // to 1, the driver's own default, so an omitted argument reproduces the + // paper's unit-constant allocation exactly. + const double adaptive_kc = (inputs.size() >= 23) ? read_double(inputs[22], "adaptive_k_const") : 1.0; + const double adaptive_sc = (inputs.size() >= 24) ? read_double(inputs[23], "adaptive_s_const") : 1.0; + const double cap_rank_fraction = (inputs.size() >= 25) ? read_double(inputs[24], "cap_rank_fraction") : 1.0; + const double quadrature_fraction = (inputs.size() >= 26) ? read_double(inputs[25], "quadrature_fraction") : 1.0; + const double first_row_ql = (inputs.size() >= 27) ? read_double(inputs[26], "first_row_ql") : 0.0; + if (first_row_ql != 0.0 && first_row_ql != 1.0) + raise("randlapack:fun_nystrom_pp_mex:first_row_ql", "first_row_ql must be 0 or 1"); + const double spend_cap = (inputs.size() >= 28) ? read_double(inputs[27], "spend_cap") : 0.0; + if (spend_cap != 0.0 && spend_cap != 1.0) + raise("randlapack:fun_nystrom_pp_mex:spend_cap", "spend_cap must be 0 or 1"); + if (spend_cap != 0.0 && !target_auto) + raise("randlapack:fun_nystrom_pp_mex:spend_cap", + "spend_cap applies only to lfa_type 'scalar_auto' or 'block_auto'"); + + if (!std::isfinite(adaptive_kc) || !std::isfinite(adaptive_sc) || + !(adaptive_kc > 0.0) || !(adaptive_sc > 0.0)) { + raise("randlapack:fun_nystrom_pp_mex:adaptive_split_const", + "adaptive_k_const (input 23) and adaptive_s_const (input 24) must be > 0"); + } + if (adaptive_fl < 0 || adaptive_fl > 2) { + raise("randlapack:fun_nystrom_pp_mex:adaptive", + "adaptive (input 13) must be 0 (fixed depth), 1 (Radau-certified), " + "or 2 (legacy window; block_qfa only)"); + } + if (radau_ret != 0 && radau_ret != 1) { + raise("randlapack:fun_nystrom_pp_mex:radau_return", + "radau_return (input 19) must be 0 (Gauss) or 1 (midpoint)"); + } + if (auto_dcap < 0) { + raise("randlapack:fun_nystrom_pp_mex:auto_depth_cap", + "auto_depth_cap (input 20) must be >= 0 (0 = no fixed cap)"); + } + if (!(auto_pfrac > 0.0 && auto_pfrac < 1.0)) { + raise("randlapack:fun_nystrom_pp_mex:auto_probe_frac", + "auto_probe_frac (input 21) must lie in (0, 1)"); + } + if (adaptive_mvcap < 0) { + raise("randlapack:fun_nystrom_pp_mex:adaptive_matvec_cap", + "adaptive_matvec_cap (input 22) must be >= 0 (0 = no cap)"); + } + if (vec_nnz_in < 0) { + raise("randlapack:fun_nystrom_pp_mex:vec_nnz", + "vec_nnz (input 10) must be >= 0 (0 = auto, ~log(k))"); + } + + if (inputs.size() >= 9 && sketch_type != "saso" && sketch_type != "gaussian") { + matlabPtr->feval(u"warning", 0, std::vector{ + factory.createCharArray("randlapack:fun_nystrom_pp_mex:sketch_type"), + factory.createCharArray( + "sketch_type '" + sketch_type + "' ignored: the Phase-1 sketch " + "is always a kernel-internal SASO now.") + }); + } + + // Ignored-knob warnings for 'auto'/'adaptive' (q hardcoded to 1 + // internally; Adaptive*/Depth/Reorth/RadauReturn all unread by + // these two tiers) used to live here as a value-comparison against + // the .m wrapper's own defaults -- which could not distinguish "the + // caller never touched this" from "the caller explicitly passed + // the default value", and fired on every call from a caller that + // always forwards a value verbatim (e.g. this repo's own MATLAB + // harness). Moved to randlapack.fun_nystrom_pp.m, which has access + // to inputParser's UsingDefaults and can tell the two apart; this + // low-level MEX entry point no longer warns on ignored knobs at all + // (matches its "low-level" framing in the header comment above). + + // --- Phase-1 sketch (input 2): scalar k. The sketch itself (SASO) is + // generated inside NystromEVD from (sketch_seed, vec_nnz); explicit + // matrices were rejected by the validation block in operator(). --- + const int64_t k = read_int(inputs[1], "k"); + + // --- Phase-2 Hutchinson probes (input 3): scalar s -> generated INSIDE + // the kernel (FunNystromPP draws a Gaussian n x s block from `state` and + // normalizes each column to ‖·‖₂ = √n; O2_buf stays nullptr). A matrix + // arg overrides with an explicit Omega2. Skipped if k == n. --- + const Array& O2_in = inputs[2]; + const bool phase2_skipped = (k == n); + int64_t s = 0; + std::unique_ptr O2_own; // RAII: freed on every exit path (copy path only) + const T* O2_buf = nullptr; + // 'auto'/'adaptive' never read Omega2/O2_buf: their driver.call + // overloads (below) don't even take an Omega2 argument, generating + // their own internal probes instead. Skip the marshal entirely for + // those two tiers, including the O(n*s) heap copy branch below, + // rather than paying a wasted copy of a caller-supplied explicit + // Omega2 (e.g. the "same probes for every estimator" comparison + // escape hatch) on every such call. + const bool skip_omega2_marshal = + (lfa_type == "auto" || lfa_type == "adaptive"); + if (!phase2_skipped && !skip_omega2_marshal) { + if (O2_in.getNumberOfElements() == 1) { + s = read_int(O2_in, "s"); + if (s < 1) { + raise("randlapack:fun_nystrom_pp_mex:s_range", + "s (arg 3 scalar) must be >= 1 when k < n"); + } + // O2_buf = nullptr -> the driver generates the normalized + // Gaussian probes internally (point 3 of the 2026-07-09 plan). + } else { + s = static_cast(O2_in.getDimensions()[1]); + if (s < 1) { + // An n x 0 array would otherwise flow into t2 = .../0 = NaN + // with no error raised. + raise("randlapack:fun_nystrom_pp_mex:Omega2_empty", + "explicit Omega2 (arg 3 matrix) must have s >= 1 " + "columns when k < n; got an n x 0 array"); + } + // The driver reads Omega2 through const T* (FunNystromPP::call's Omega2 and every oracle's B). + if (copy_input) { + O2_own = copy_into(O2_in, static_cast(n) * s); + O2_buf = O2_own.get(); + } else { + O2_buf = matlab_data(O2_in); + } + } + } + +#ifndef RANDLAPACK_SYMMETRIC_SKETCH + // The legacy general sketch product requires both triangles. Legacy builds always copy (copy_input_requested), + // so this writes the MEX's own buffer, never MATLAB's. + if (!vector_mode && !sparse_mode) { + T* A_w = A_own.get(); + for (int64_t j = 0; j < n; ++j) + for (int64_t i = j + 1; i < n; ++i) + A_w[i + j * n] = A_w[j + i * n]; + } +#endif + + const double marshal_in_ms = std::chrono::duration( + std::chrono::steady_clock::now() - t_marshal_start).count(); + + // --- scalar function f --- + std::function fscalar; + if (func == "sqrt") fscalar = [](T x) { return std::sqrt(std::max(x, (T)0)); }; + else if (func == "log") fscalar = [](T x) { return std::log(x + (T)1); }; // log(x+1) = tr(log(A+I)) + else if (func == "poly") fscalar = [poly_lambda](T x) { return x * (x + poly_lambda); }; + else if (func == "effdim") fscalar = [poly_lambda](T x) { return x / (x + poly_lambda); }; + else if (func == "square") fscalar = [](T x) { return x * x; }; + else if (func == "identity") fscalar = [](T x) { return x; }; + else if (func == "entropy") fscalar = [](T x) { return x > (T)0 ? -x * std::log(x) : (T)0; }; + else { + raise("randlapack:fun_nystrom_pp_mex:func", + "unknown func '" + func + "' (use sqrt|log|poly|effdim|square|identity|entropy)"); + } + + // --- f(A)*X oracle --- + using FAFun = std::function; + RandLAPACK::LanczosFA scalar_lfa; + RandLAPACK::LanczosQFA scalar_qfa; + RandLAPACK::BlockLanczosFA block_lfa; + RandLAPACK::BlockLanczosQFA block_qfa; + // Lanczos-QFA fills the quadratic form Ω₂ᵀf(A)Ω₂ directly (no f(A)·Ω₂ + // mapback); the driver takes its trace, reading only the diagonal. + // Signalled to the driver below. + const bool qfa_mode = (lfa_type == "block_qfa" || lfa_type == "scalar_qfa"); + + // --- drive funNystrom++ --- + // The driver and the four oracle objects are operator-independent and + // are read after the tier body (the timing outputs below), so they are + // declared here and the tier body is a generic lambda over the operator + // type: one implementation serves the dense and the diagonal operator. + RandLAPACK::FunNystromPP driver; + // Sub-timers always on: a handful of steady_clock reads per call, + // negligible next to any BLAS work; consumed by the 4th output. + driver.nystrom_ws.times_enabled = true; + scalar_lfa.timing = true; + scalar_qfa.timing = true; + block_lfa.timing = true; + block_qfa.timing = true; + driver.auto_sqfa.timing = true; + scalar_qfa.use_first_row_ql = (first_row_ql != 0.0); + driver.auto_sqfa.use_first_row_ql = (first_row_ql != 0.0); + driver.adaptive_bqfa.timing = true; + scalar_lfa.reorth = reorth_flag; + block_lfa.reorth = reorth_flag; + block_qfa.reorth = reorth_flag; + // scalar_qfa: intrinsically no-reorth (basis-free recurrence); reorth / + // adaptive_delay / adaptive_min do not apply: the Gauss-Radau + // certificate has no window. adaptive_tol is its CERTIFIED tolerance. + // Any nonzero adaptive selects it (the scalar class has no window rule). + scalar_qfa.adaptive = (adaptive_fl != 0); + scalar_qfa.adaptive_rtol = adaptive_tl; + // block_qfa adaptive semantics: 0 = fixed depth; 1 = Radau-certified + // (consistent with scalar_qfa's meaning); 2 = legacy window rule + // (honors adaptive_delay / adaptive_min). stop_rule and return_mode + // are assigned unconditionally: the oracle may be cached across calls + // (persistent-handle path), so every mode must be set by name, never + // inherited from a previous call. + block_qfa.adaptive = (adaptive_fl != 0); + block_qfa.adaptive_rtol = adaptive_tl; + block_qfa.stop_rule = (adaptive_fl == 2) + ? RandLAPACK::BlockQFAStop::Window : RandLAPACK::BlockQFAStop::Radau; + block_qfa.return_mode = (radau_ret == 1) + ? RandLAPACK::BlockQFAReturn::Midpoint : RandLAPACK::BlockQFAReturn::Gauss; + // Assign UNCONDITIONALLY, restoring the library default by name when the + // caller passed 0. The former `if (x > 0)` form depended on the oracle + // being freshly constructed every call to supply the default; once the + // oracle is cached (persistent-handle path) that assumption is false and + // the previous call's window leaks forward, changing d_used and hence + // both the estimate and the matvec count. The oracle is non-copyable and + // non-assignable, so "reset by reconstruction" is not available. + block_qfa.adaptive_delay = (adaptive_dl > 0) + ? adaptive_dl : RandLAPACK::BlockLanczosQFA::default_adaptive_delay; + block_qfa.adaptive_min = (adaptive_mn > 0) + ? adaptive_mn : RandLAPACK::BlockLanczosQFA::default_adaptive_min; + driver.vec_nnz = vec_nnz; + driver.use_qfa = qfa_mode; + driver.auto_depth_cap = auto_dcap; + driver.auto_probe_frac = static_cast(auto_pfrac); + // A +-1 quadratic form is exact on a diagonal matrix, which would make + // the eps-targeted tier's error degenerate rather than sampled, so + // vector mode draws that tier's probes from the sphere instead. + driver.adaptive_rademacher = !target_auto && !vector_mode; + driver.adaptive_use_scalar = (requested_lfa == "scalar_auto"); + driver.adaptive_reuse_pilot = !target_auto; + driver.adaptive_gauss_return = target_auto; + driver.adaptive_k_const = static_cast(adaptive_kc); + driver.adaptive_s_const = static_cast(adaptive_sc); + driver.adaptive_cap_rank_fraction = static_cast(cap_rank_fraction); + driver.adaptive_quadrature_fraction = static_cast(quadrature_fraction); + driver.adaptive_spend_cap = (spend_cap != 0.0); + + T t1 = (T)0, t2 = (T)0, est = (T)0; + const T *Omega2_ptr = phase2_skipped ? nullptr : O2_buf; + RandBLAS::RNGState state(static_cast(sketch_seed)); + // f(lambda) for the vector-mode exact oracle; unused otherwise. + std::unique_ptr flam_own; + + auto drive_tiers = [&](auto& A_op) { + FAFun fAfun; + if (lfa_type == "exact") { + if (vector_mode) { + // f(A)B is a row scaling by f(lambda). make_exact_fa_oracle's + // syevd needs the dense A and is not called here. + flam_own.reset(new T[n]); + for (int64_t i = 0; i < n; ++i) flam_own[i] = fscalar(A_buf[i]); + fAfun = [flam = flam_own.get()] + (int64_t m_, int64_t s_, const T *B, T *Y) { + for (int64_t jj = 0; jj < s_; ++jj) + for (int64_t ii = 0; ii < m_; ++ii) + Y[ii + jj * m_] = flam[ii] * B[ii + jj * m_]; + }; + } else if (sparse_mode) { + raise("randlapack:fun_nystrom_pp_mex:exact_sparse", + "lfa_type 'exact' needs the dense matrix (one-shot syevd); it is not available for sparse A"); + } else { + // V*diag(f(lambda))*V^T*B via a one-shot syevd. Shared implementation + // with the RandLAPACK test + benchmark (single point of correctness). + fAfun = testing::make_exact_fa_oracle(n, A_buf, fscalar); + } + } else if (lfa_type == "scalar") { + fAfun = [&scalar_lfa, &A_op, &fscalar, d] + (int64_t m_, int64_t s_, const T *B, T *Y) { + scalar_lfa.call(A_op, B, m_, s_, fscalar, d, Y); + }; + } else if (lfa_type == "block") { + fAfun = [&block_lfa, &A_op, &fscalar, d] + (int64_t m_, int64_t s_, const T *B, T *Y) { + block_lfa.call(A_op, B, m_, s_, fscalar, d, Y); + }; + } else if (lfa_type == "scalar_qfa") { + // Per-column quadratic forms onto the DIAGONAL of the s×s Y (the + // only part the driver's use_qfa trace read touches; off-diagonals + // are left unwritten — same convention as the auto tier). + fAfun = [&scalar_qfa, &A_op, &fscalar, d] + (int64_t m_, int64_t s_, const T *B, T *Y) { + // unique_ptr: an exception out of scalar_qfa.call must not + // leak the per-call buffer. + std::unique_ptr qv(new T[s_]); + scalar_qfa.call(A_op, B, m_, s_, fscalar, d, qv.get()); + for (int64_t jj = 0; jj < s_; ++jj) + Y[jj + jj * s_] = qv[jj]; + }; + } else if (lfa_type == "block_qfa") { + // Y is the s×s matrix M = Ω₂ᵀ f(A) Ω₂ (driver sizes fAOmega to s²). + fAfun = [&block_qfa, &A_op, &fscalar, d] + (int64_t m_, int64_t s_, const T *B, T *Y) { + block_qfa.call(A_op, B, m_, s_, fscalar, d, Y); + }; + } else if (lfa_type == "auto") { + // Knob-free tier: no oracle to build here; the driver's auto call + // owns its QFA oracle and picks k/s/depth from (budget, auto_eps). + // Same id as randlapack.fun_nystrom_pp.m's own pre-check (unified: + // previously this used a different, undocumented + // 'randlapack:fun_nystrom_pp_mex:budget' id, unreachable through + // the documented .m API since .m always validates first; a direct + // fun_nystrom_pp_mex(...) call now raises the same id the .m + // wrapper would have). + if (auto_budget < 1) + raise("randlapack:fun_nystrom_pp:Budget", + "lfa_type 'auto' requires a positive matvec budget (input 17)"); + } else if (lfa_type == "adaptive") { + // Eps-targeted knob-free tier: no oracle to build here either; the + // driver's eps-targeted call owns its BLOCK QFA oracle and picks + // k/s/depth from a certified depth probe at rtol = auto_eps. No + // upfront budget is required (unlike 'auto'); adaptive_matvec_cap + // (0 = none) is passed straight through below. + } else { + raise("randlapack:fun_nystrom_pp_mex:lfa_type", + "unknown lfa_type '" + lfa_type + + "' (use exact|scalar|scalar_qfa|block|block_qfa|auto|adaptive)"); + } + + if (lfa_type == "auto") { + // The auto overload throws std::invalid_argument for an infeasible + // matvec budget. Surface that as a dedicated MATLAB error id so + // the benchmark can catch it and write a skip row; other + // invalid_argument reasons (eps / probe_frac range) keep the + // generic StdError id. + try { + est = driver.call(A_op, fscalar, auto_budget, static_cast(auto_eps), + state, t1, t2); + } catch (const std::invalid_argument& e) { + const std::string msg = e.what(); + if (msg.find("infeasible") != std::string::npos) { + raise("randlapack:fun_nystrom_pp:infeasibleBudget", msg); + } + raise("randlapack:fun_nystrom_pp_mex:StdError", msg); + } + } else if (lfa_type == "adaptive") { + // Eps-targeted overload: auto_eps doubles as eps, auto_dcap (shared + // with the 'auto' tier's driver.auto_depth_cap member, already set + // above) bounds the probe depth. adaptive_mvcap == 0 => nullopt (no + // cap). The driver throws std::invalid_argument for two distinct + // infeasibility cases: the matvec_cap case and the block-Krylov + // s*t <= n case reachable at non-default probe block widths. Both + // driver messages contain "infeasible", so the OR-match above on + // "s*t <= n" is redundant-but-harmless; it is kept for robustness + // in case the driver message text changes. Both cases route to + // the SAME MATLAB error id as the 'auto' tier's infeasible-budget + // case, so callers can catch one id regardless of which + // knob-free tier they picked. + const std::optional mv_cap = (adaptive_mvcap > 0) + ? std::optional(adaptive_mvcap) : std::nullopt; + try { + est = driver.call(A_op, fscalar, static_cast(auto_eps), + state, t1, t2, mv_cap); + } catch (const std::invalid_argument& e) { + const std::string msg = e.what(); + if (msg.find("infeasible") != std::string::npos || + msg.find("s*t <= n") != std::string::npos) { + raise("randlapack:fun_nystrom_pp:infeasibleBudget", msg); + } + raise("randlapack:fun_nystrom_pp_mex:StdError", msg); + } + } else { + est = driver.call(A_op, fAfun, fscalar, + k, s, q, + state, Omega2_ptr, + t1, t2); + } + }; // end drive_tiers + + if (vector_mode) { + linops::DiagSymLinOp A_op(n, A_buf); + drive_tiers(A_op); + } else if (sparse_mode) { + if constexpr (std::is_same_v) { + using CSC = RandBLAS::sparse_data::CSCMatrix; + // CSCMatrix's non-owning constructor takes non-const pointers; the view is only read (see marshal_csc), + // so dropping const here never lets a write reach MATLAB's arrays. + CSC A_csc(n, n, sp_nnz, const_cast(sp_vals), const_cast(sp_rowidx), + const_cast(sp_colptr)); + linops::SparseSymLinOp A_op(A_csc); + drive_tiers(A_op); + } else { + raise("randlapack:fun_nystrom_pp_mex:A_sparse_single", "sparse A is double only"); + } + } else { + linops::ExplicitSymLinOp A_op(n, blas::Uplo::Upper, A_buf, n, Layout::ColMajor); + A_op.both_triangles = true; // A_buf is MATLAB's full symmetric matrix (or, with copy_input, a copy of it) + drive_tiers(A_op); + } + + outputs[0] = factory.createScalar(static_cast(est)); + if (outputs.size() >= 2) outputs[1] = factory.createScalar(static_cast(t1)); + if (outputs.size() >= 3) outputs[2] = factory.createScalar(static_cast(t2)); + + // --- optional 4th output: wall-clock instrumentation struct --- + if (outputs.size() >= 4) { + // NystromEVD 11-slot breakdown (microseconds, as doubles). + std::vector nys_us(driver.nystrom_ws.times.begin(), + driver.nystrom_ws.times.end()); + if (nys_us.size() != 11) nys_us.assign(11, 0.0); + // Lanczos-oracle 6-slot breakdown (whichever oracle ran); zeros for + // the exact oracle, which has no Lanczos phase to instrument. + std::vector lfa_us(6, 0.0); + if (lfa_type == "scalar" && scalar_lfa.times.size() >= 5) { + lfa_us.assign(scalar_lfa.times.begin(), scalar_lfa.times.end()); + } else if (lfa_type == "scalar_qfa" && scalar_qfa.times.size() >= 5) { + lfa_us.assign(scalar_qfa.times.begin(), scalar_qfa.times.end()); + } else if (lfa_type == "block" && block_lfa.times.size() >= 5) { + lfa_us.assign(block_lfa.times.begin(), block_lfa.times.end()); + } else if (lfa_type == "block_qfa" && block_qfa.times.size() >= 5) { + lfa_us.assign(block_qfa.times.begin(), block_qfa.times.end()); + } else if (lfa_type == "auto" && driver.auto_sqfa.times.size() >= 5) { + lfa_us.assign(driver.auto_sqfa.times.begin(), driver.auto_sqfa.times.end()); + } else if (requested_lfa == "scalar_auto" && driver.auto_sqfa.times.size() >= 5) { + lfa_us.assign(driver.auto_sqfa.times.begin(), driver.auto_sqfa.times.end()); + } else if (lfa_type == "adaptive" && driver.adaptive_bqfa.times.size() >= 5) { + lfa_us.assign(driver.adaptive_bqfa.times.begin(), driver.adaptive_bqfa.times.end()); + } + // NaN convention (see the header comment): used below for both the + // path-specific fields and 'exact''s d_used (declared here, before + // its first use, rather than where the path-specific block below + // used to declare it). + const double nan_v = std::numeric_limits::quiet_NaN(); + // Lanczos depth actually used by the f(A) oracle. For block_qfa with + // adaptive stopping this is the online-chosen depth (< the d cap); + // for 'auto' it is the probe-discovered depth t; for 'adaptive' it is + // the depth probe's certified (or reached) depth t + // (driver.adaptive_t); for every other lfa_type it is just the fixed + // d, EXCEPT 'exact', which never runs a Lanczos recurrence at all (it + // builds f(A) once via syevd) and so has no meaningful depth: NaN, + // not the leftover default-depth input `d`, which would otherwise + // read as a real (if constant and meaningless) cost figure in a + // benchmark plot. Exposed so the benchmark can count matvecs + // (matvecs of A are proportional to this depth, for oracles that + // have one). + double d_used; + if (lfa_type == "block_qfa") d_used = static_cast(block_qfa.d_used); + else if (lfa_type == "scalar_qfa") d_used = static_cast(scalar_qfa.d_used); + else if (lfa_type == "auto") d_used = static_cast(driver.auto_t); + else if (lfa_type == "adaptive") d_used = static_cast(driver.adaptive_t); + else if (lfa_type == "exact") d_used = nan_v; + else d_used = static_cast(d); + // Actual Phase-2 oracle matvecs: Σ per-probe certified depths for + // the scalar-QFA-backed types, s*d_used (the class's matvecs + // member) for block_qfa, the block oracle's matvecs member for + // 'adaptive'; 0 otherwise (analytic s*d_used). + double oracle_mv = 0.0; + if (lfa_type == "scalar_qfa") oracle_mv = static_cast(scalar_qfa.matvecs); + else if (lfa_type == "block_qfa") oracle_mv = static_cast(block_qfa.matvecs); + else if (lfa_type == "auto") oracle_mv = static_cast(driver.auto_oracle_matvecs); + else if (lfa_type == "adaptive") oracle_mv = static_cast(driver.adaptive_oracle_matvecs); + // Knob-free bookkeeping (zeros unless lfa_type == 'auto'/'adaptive'): + // the chosen rank / probe count and the matvecs the depth probe + // actually spent. With the certified oracle the budget closes as an + // upper bound: probe_mv + q*auto_k + oracle_mv <= budget (q = 1 in + // the auto tier; 'adaptive' has no upfront budget so this is an + // accounting identity, not a constraint enforced against an input). + const bool is_auto = (lfa_type == "auto"); + const bool is_adaptive = (lfa_type == "adaptive"); + const double auto_k_out = is_auto ? static_cast(driver.auto_k) + : is_adaptive ? static_cast(driver.adaptive_k) : 0.0; + const double auto_s_out = is_auto ? static_cast(driver.auto_s) + : is_adaptive ? static_cast(driver.adaptive_s) : 0.0; + const double probe_mv = is_auto ? static_cast(driver.auto_probe_matvecs) + : is_adaptive ? static_cast(driver.adaptive_probe_matvecs) : 0.0; + // Path-specific fields, NaN when the path did not run (see the + // header comment): block_qfa's/adaptive's Gauss/Radau traces + + // certification, and the auto/adaptive tiers' probe wall-clock + + // certification flags. nan_v declared above, next to d_used. + double tr_U_out = nan_v, tr_L_out = nan_v, cert_out = nan_v; + if (lfa_type == "block_qfa") { + tr_U_out = static_cast(block_qfa.tr_U); + tr_L_out = static_cast(block_qfa.tr_L); + cert_out = block_qfa.certified ? 1.0 : 0.0; + } else if (is_adaptive) { + tr_U_out = static_cast(driver.adaptive_tr_U); + tr_L_out = static_cast(driver.adaptive_tr_L); + cert_out = driver.adaptive_phase2_certified ? 1.0 : 0.0; + // scalar_auto: the driver's block traces are NaN by construction. The + // scalar oracle keeps per-column Gauss/Radau values for its last call, + // which is the Phase-2 call over driver.adaptive_s columns; report + // their sums so the bracket is available for every self-tuning tier. + if (driver.adaptive_use_scalar && driver.adaptive_s > 0 && + driver.auto_sqfa.gauss_val_sz >= driver.adaptive_s) { + double su = 0.0, sl = 0.0; + for (int64_t j = 0; j < driver.adaptive_s; ++j) { + su += static_cast(driver.auto_sqfa.gauss_val[j]); + sl += static_cast(driver.auto_sqfa.radau_val[j]); + } + tr_U_out = su; tr_L_out = sl; + } + // spend_cap: Phase 2 ran at fixed depth and evaluated no bracket, so the lower + // side and the certificate are not measurements (the oracles leave 0 or tr_U there). + if (!driver.adaptive_phase2_checked) { tr_L_out = nan_v; cert_out = nan_v; } + } else if (lfa_type == "scalar_qfa" && scalar_qfa.adaptive) { + // NaN when Adaptive=0 (fixed depth): no certificate was ever + // checked, so all_certified's default-false would misreport + // "uncertified" rather than "not applicable". + cert_out = scalar_qfa.all_certified ? 1.0 : 0.0; + // Sum of the per-column Gauss and Gauss-Radau values over the s + // Phase-2 columns: the scalar analogue of the block tier's tr_U/tr_L. + if (s > 0 && scalar_qfa.gauss_val_sz >= s) { + double su = 0.0, sl = 0.0; + for (int64_t j = 0; j < s; ++j) { + su += static_cast(scalar_qfa.gauss_val[j]); + sl += static_cast(scalar_qfa.radau_val[j]); + } + tr_U_out = su; tr_L_out = sl; + } + } + double probe_ms_out = nan_v, probe_conv_out = nan_v, ph2_cert_out = nan_v, ph2_checked_out = nan_v; + if (is_auto) { + probe_ms_out = driver.t_probe_ms; + probe_conv_out = driver.auto_probe_converged ? 1.0 : 0.0; + ph2_cert_out = driver.auto_phase2_certified ? 1.0 : 0.0; + } else if (is_adaptive) { + probe_ms_out = driver.t_adaptive_probe_ms; + probe_conv_out = driver.adaptive_probe_certified ? 1.0 : 0.0; + ph2_cert_out = driver.adaptive_phase2_checked ? (driver.adaptive_phase2_certified ? 1.0 : 0.0) : nan_v; + ph2_checked_out = driver.adaptive_phase2_checked ? 1.0 : 0.0; + } + matlab::data::StructArray ts = factory.createStructArray({1, 1}, + {"marshal_in_ms", "phase1_ms", "phase2_ms", "fafun_ms", + "assembly_ms", "specrec_ms", "nystrom_us", "lfa_us", "d_used", + "oracle_mv", "auto_k", "auto_s", "probe_mv", + "tr_U", "tr_L", "certified", + "probe_ms", "probe_converged", "phase2_certified", + "first_row_ql_requested", "first_row_ql_fallback", "symmetric_sketch", + "ritz_clamped", "rank_deficient_steps", "min_diag_ratio", + "nystrom_clamped", "bracket_evaluated", "phase2_checked", "input_copied"}); + ts[0]["marshal_in_ms"] = factory.createScalar(marshal_in_ms); + ts[0]["phase1_ms"] = factory.createScalar(driver.t_phase1_ms); + ts[0]["phase2_ms"] = factory.createScalar(driver.t_phase2_ms); + ts[0]["fafun_ms"] = factory.createScalar(driver.t_fafun_ms); + ts[0]["assembly_ms"] = factory.createScalar(driver.t_phase2_ms - driver.t_fafun_ms); + ts[0]["specrec_ms"] = factory.createScalar(driver.t_specrec_ms); + ts[0]["nystrom_us"] = factory.createArray({1, nys_us.size()}, + nys_us.data(), nys_us.data() + nys_us.size()); + ts[0]["lfa_us"] = factory.createArray({1, lfa_us.size()}, + lfa_us.data(), lfa_us.data() + lfa_us.size()); + ts[0]["d_used"] = factory.createScalar(d_used); + ts[0]["oracle_mv"] = factory.createScalar(oracle_mv); + ts[0]["auto_k"] = factory.createScalar(auto_k_out); + ts[0]["auto_s"] = factory.createScalar(auto_s_out); + ts[0]["probe_mv"] = factory.createScalar(probe_mv); + ts[0]["tr_U"] = factory.createScalar(tr_U_out); + ts[0]["tr_L"] = factory.createScalar(tr_L_out); + ts[0]["certified"] = factory.createScalar(cert_out); + ts[0]["probe_ms"] = factory.createScalar(probe_ms_out); + ts[0]["probe_converged"] = factory.createScalar(probe_conv_out); + ts[0]["phase2_certified"] = factory.createScalar(ph2_cert_out); + ts[0]["phase2_checked"] = factory.createScalar(ph2_checked_out); + // true when A (and an explicit Omega2) were copied into MEX-owned buffers (RANDLAPACK_FNPP_COPY_INPUT=1 or a + // legacy build); false when the driver read MATLAB's arrays in place. + ts[0]["input_copied"] = factory.createScalar(copy_input); + ts[0]["first_row_ql_requested"] = factory.createScalar(first_row_ql != 0.0); + ts[0]["first_row_ql_fallback"] = factory.createScalar( + scalar_qfa.first_row_ql_fallback_seen || driver.auto_sqfa.first_row_ql_fallback_seen); +#ifdef RANDLAPACK_SYMMETRIC_SKETCH + ts[0]["symmetric_sketch"] = factory.createScalar(true); +#else + ts[0]["symmetric_sketch"] = factory.createScalar(false); +#endif + // ---- Numerical-degradation diagnostics ---- + // Until now the kernel's guards were unobservable: a clamped Ritz value, a + // rank-deficient block QR or a clamped Nystrom eigenvalue left no trace in any + // exported field, so a degraded run and a healthy one produced identical telemetry. + // These five fields are the channel out. Exactly one code path runs per call, so + // summing across the objects reports that path's totals and leaves the rest at zero. + const double ritz_clamped_total = + (double)scalar_qfa.ritz_clamped + + (double)block_qfa.ritz_clamped + (double)block_qfa.fa.ritz_clamped + + (double)driver.auto_sqfa.ritz_clamped + + (double)driver.adaptive_bqfa.ritz_clamped + + (double)driver.adaptive_bqfa.fa.ritz_clamped; + const double rank_def_total = + (double)block_qfa.fa.rank_deficient_steps + + (double)driver.adaptive_bqfa.fa.rank_deficient_steps; + // Smallest diag ratio seen by whichever block QR ran; stays Inf when none did, so a + // reader can tell "no block QR" from "a block QR that was comfortably full rank". + double min_ratio = (double)block_qfa.fa.min_diag_ratio; + if ((double)driver.adaptive_bqfa.fa.min_diag_ratio < min_ratio) + min_ratio = (double)driver.adaptive_bqfa.fa.min_diag_ratio; + ts[0]["ritz_clamped"] = factory.createScalar(ritz_clamped_total); + ts[0]["rank_deficient_steps"] = factory.createScalar(rank_def_total); + ts[0]["min_diag_ratio"] = factory.createScalar(min_ratio); + ts[0]["nystrom_clamped"] = factory.createScalar((double)driver.nystrom_ws.clamped_eigenvalues); + ts[0]["bracket_evaluated"] = factory.createScalar( + block_qfa.bracket_evaluated || driver.adaptive_bqfa.bracket_evaluated); + outputs[3] = std::move(ts); + } + // A_own / O2_own release their buffers here (RAII). + } + +public: + MexFunction() : matlabPtr(getEngine()) {} + + void operator()(ArgumentList outputs, ArgumentList inputs) { + try { + if (inputs.size() < 8 || inputs.size() > 28) { + raise("randlapack:fun_nystrom_pp_mex:nargin", + "Expected 8 to 28 inputs: A, Omega1, Omega2, func, q, poly_lambda, " + "lfa_type, d [, sketch_type, vec_nnz, sketch_seed, reorth, " + "adaptive, adaptive_tol, adaptive_delay, adaptive_min, " + "budget, auto_eps, radau_return, auto_depth_cap, " + "auto_probe_frac, adaptive_matvec_cap, adaptive_k_const, " + "adaptive_s_const, cap_rank_fraction, quadrature_fraction, first_row_ql, spend_cap]"); + } + if (outputs.size() < 1 || outputs.size() > 4) { + raise("randlapack:fun_nystrom_pp_mex:nargout", + "Expected 1 to 4 outputs: [est, t1, t2, times]"); + } + + // --- A: square n x n, OR an n-vector read as diag(A), real, + // single/double. Vector mode never forms the n x n matrix. --- + const Array& A_in = inputs[0]; + auto A_dims = A_in.getDimensions(); + if (A_in.getType() == ArrayType::SPARSE_DOUBLE) { + raise("randlapack:fun_nystrom_pp_mex:A_sparse_raw", + "pass sparse A through randlapack.fun_nystrom_pp, which hands the MEX its compressed-column arrays " + "(reading a MATLAB sparse array element by element here costs O(n) per nonzero)"); + } + const bool sparse_in = is_sparse_struct(A_in); + const bool vector_mode = !sparse_in && (A_dims.size() == 2) && + ((A_dims[0] == 1 && A_dims[1] > 1) || (A_dims[1] == 1 && A_dims[0] > 1)); + if (!sparse_in && (A_dims.size() != 2 || (!vector_mode && A_dims[0] != A_dims[1]))) { + raise("randlapack:fun_nystrom_pp_mex:A_shape", + "A must be a square n x n matrix or an n-vector of diagonal entries"); + } + const ArrayType A_type = sparse_in ? ArrayType::DOUBLE : A_in.getType(); + if (A_type != ArrayType::DOUBLE && A_type != ArrayType::SINGLE) { + raise("randlapack:fun_nystrom_pp_mex:A_dtype", + "A must be single or double precision, or a sparse double matrix"); + } + // A_dims[0] is 1 for a row vector, so numel is the only safe read. + const int64_t n = sparse_in ? sparse_n_of(A_in) + : vector_mode ? static_cast(A_in.getNumberOfElements()) + : static_cast(A_dims[0]); + + // --- arg 2 (Phase-1): scalar k only. The Phase-1 sketch is a SASO + // generated inside RandLAPACK::NystromEVD from (sketch_seed, vec_nnz); + // the explicit-matrix escape hatch was removed with the kernel's + // dense-sketch mode. --- + const Array& O1_in = inputs[1]; + if (O1_in.getNumberOfElements() != 1) { + raise("randlapack:fun_nystrom_pp_mex:Omega1_explicit", + "arg 2 must be a scalar k. Explicit Omega1 matrices are no " + "longer supported: the Phase-1 sketch is a SASO " + "generated inside RandLAPACK (control it via sketch_seed and " + "vec_nnz)."); + } + const int64_t k = read_int(O1_in, "k"); + if (k < 1 || k > n) { + raise("randlapack:fun_nystrom_pp_mex:k_range", + "k (arg 2 scalar) must satisfy 1 <= k <= n"); + } + + // --- arg 3 (Phase-2): scalar s OR explicit n x s matrix; only when k < n --- + const Array& O2_in = inputs[2]; + if (k != n && O2_in.getNumberOfElements() != 1) { + auto O2_dims = O2_in.getDimensions(); + if (O2_dims.size() != 2 || static_cast(O2_dims[0]) != n) { + raise("randlapack:fun_nystrom_pp_mex:Omega2_shape", + "Omega2 (arg 3 matrix) must be n x s when k < n"); + } + if (O2_in.getType() != A_type) { + raise("randlapack:fun_nystrom_pp_mex:Omega2_dtype", + "Omega2 must be the same class (single/double) as A"); + } + } + + // --- dispatch on precision --- + if (A_type == ArrayType::DOUBLE) { + run_fun_nystrom(outputs, inputs); + } else { + run_fun_nystrom(outputs, inputs); + } + } + catch (const MexError& r) { + emit_matlab_error(r.id, r.msg); + } + catch (const RandLAPACK::Error& e) { + emit_matlab_error("randlapack:fun_nystrom_pp_mex:RandLAPACKError", e.what()); + } + catch (const RandBLAS::Error& e) { + emit_matlab_error("randlapack:fun_nystrom_pp_mex:RandBLASError", e.what()); + } + catch (const matlab::Exception&) { + // Engine exceptions (e.g. an interrupt) propagate unchanged. NB in + // the dual-libstdc++ setup this clause may fail to match, in which + // case the std::exception clause below re-labels the error as + // StdError with the message preserved. + throw; + } + catch (const std::exception& e) { + emit_matlab_error("randlapack:fun_nystrom_pp_mex:StdError", e.what()); + } + catch (...) { + emit_matlab_error("randlapack:fun_nystrom_pp_mex:Unknown", + "Unknown exception in fun_nystrom_pp_mex"); + } + } +}; diff --git a/matlab/src/numa_place_mex.cc b/matlab/src/numa_place_mex.cc new file mode 100644 index 0000000..c2e52a0 --- /dev/null +++ b/matlab/src/numa_place_mex.cc @@ -0,0 +1,538 @@ +// MEX helper for NUMA page placement of large dense MATLAB arrays (Linux only), using MATLAB's C++ Data API. +// +// Why: Linux places an anonymous page on the memory node of the thread that first writes it. MATLAB fills a matrix it +// reads from disk (fread, load) with one thread, so every page of A lands on one node, and every multithreaded product +// with A then streams it through that node's memory controller. On ISAAC's 4-node Xeon 6248R hosts that made the +// funNystrom++ products 2.4x slower at n = 50,000 (2026-09-24, on a 4-node bigmem host, clrm1217). The fixes below spread A over the nodes while holding +// ONE copy of it: no transient second matrix, unlike re-touching it with A = A + 0 or copying it in the MEX. +// +// stats = numa_place_mex('spread', A) +// Interleaves the pages of A over the memory nodes this process may use (get_mempolicy MPOL_F_MEMS_ALLOWED), in +// place: move_pages(MPOL_MF_MOVE) with an explicit destination for every page of the page-aligned span of A's +// buffer, round robin over the nodes by 2 MiB granule (the transparent-huge-page size, so a huge page gets one +// destination), split over threads. mbind(MPOL_INTERLEAVE, MPOL_MF_MOVE) is NOT used: it leaves every page that +// already sits on one of the allowed nodes where it is (mbind(2): "pages already residing on the specified nodes +// will not be moved such that they are interleaved"), which is exactly the one-node matrix this must fix. The +// values of A are untouched (the kernel moves pages, it does not rewrite them) and nothing is copied in MATLAB. +// The span may include parts of the first and last page that belong to neighbouring allocations; moving those +// is harmless. +// hist = numa_place_mex('query', A [, max_pages]) +// Memory node of up to max_pages (default 4096) pages sampled evenly across A, via move_pages with nodes = NULL +// (a query; nothing moves). hist.counts(k+1) = pages on node k; hist.absent = pages not yet faulted in. +// [A, stats] = numa_place_mex('load', file, offset_bytes, n_rows, n_cols, divisor, n_threads, interleave) +// Reads n_rows*n_cols float64 values stored contiguously (column-major as MATLAB reads them) at offset_bytes of +// file into a NEW uninitialized MATLAB array, with n_threads threads (0 = OMP_NUM_THREADS, else all cores), each +// pread-ing its own contiguous span. With interleave = 1 (the wrapper's default) the empty buffer first gets the +// MPOL_INTERLEAVE policy (mbind, nothing to move yet), so the reads fault its pages in round robin over the +// allowed nodes, the placement numactl --interleave=all gives; with interleave = 0 the pages land on the nodes +// of the threads that first touch them (best effort: the threads are not pinned, so the scheduler decides). (On ISAAC's 6248R nodes interleave measured faster than threaded first +// touch for products with A: gemv 100 vs 77 GB/s, stream 3, 2026-10-07.) Every +// value is then divided by divisor unless divisor == 1: a correctly rounded IEEE division, the same operation as +// MATLAB's blk/divisor (this file must not be built with -ffast-math or -freciprocal-math). One copy of the +// matrix exists at peak; stats reports the resident-set growth after allocation (~0: pages untouched) and fill. +// +// The MATLAB-side wrappers are randlapack.numa_spread, randlapack.numa_query and randlapack.load_dense. + +#include "mex.hpp" +#include "mexAdapter.hpp" + +#if defined(__FAST_MATH__) +#error "numa_place_mex must not be built with -ffast-math: load_dense's division has to be the IEEE one (bit-identical to MATLAB's A/d)" +#endif + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#ifdef __linux__ +#include +#include +#include +#include +#endif + +namespace { + +struct MexError { + std::string id; + std::string msg; +}; + +#ifdef __linux__ +// Constants of , spelled out so the MEX builds without libnuma's headers. +constexpr int kMpolInterleave = 3; +constexpr unsigned kMpolMfMove = 1u << 1; +constexpr unsigned long kMpolFMemsAllowed = 1ul << 2; +constexpr unsigned long kMaxNodes = 1024; // bits in the node masks below +constexpr size_t kMaskWords = kMaxNodes / (8 * sizeof(unsigned long)); + +long sys_mbind(void* addr, unsigned long len, int mode, const unsigned long* mask, unsigned long maxnode, + unsigned flags) { + return syscall(SYS_mbind, addr, len, mode, mask, maxnode, flags); +} +long sys_get_mempolicy(int* mode, unsigned long* mask, unsigned long maxnode, void* addr, unsigned long flags) { + return syscall(SYS_get_mempolicy, mode, mask, maxnode, addr, flags); +} +long sys_move_pages(int pid, unsigned long count, void** pages, const int* nodes, int* status, int flags) { + return syscall(SYS_move_pages, pid, count, pages, nodes, status, flags); +} + +size_t page_size() { return static_cast(sysconf(_SC_PAGESIZE)); } + +double rss_mib() { + long pages = 0, rss = 0; + FILE* f = std::fopen("/proc/self/statm", "r"); + if (f) { if (std::fscanf(f, "%ld %ld", &pages, &rss) != 2) rss = 0; std::fclose(f); } + return static_cast(rss) * static_cast(page_size()) / 1048576.0; +} +#endif + +// CPUs this process may run on (a SLURM job's cpuset), not the whole host. +size_t usable_cpus() { +#ifdef __linux__ + cpu_set_t set; + CPU_ZERO(&set); + if (sched_getaffinity(0, sizeof(set), &set) == 0) { + const int c = CPU_COUNT(&set); + if (c > 0) return static_cast(c); + } +#endif + return std::max(1, std::thread::hardware_concurrency()); +} + +size_t thread_count(size_t requested, size_t work_units) { + size_t hw = usable_cpus(); + size_t nt = requested; + if (nt == 0) { + const char* ot = std::getenv("OMP_NUM_THREADS"); + const int v = ot ? std::atoi(ot) : 0; // unset, empty, zero or negative: all usable CPUs + nt = v > 0 ? static_cast(v) : hw; + } + return std::max(1, std::min({nt, hw, std::max(1, work_units)})); +} + +// Runs body(t) for t = 0..nt-1, t = 0 on the calling thread. A failed thread creation joins the started threads and +// runs the remaining indices on the calling thread, so no joinable std::thread is ever destroyed (that would +// terminate MATLAB). An exception out of body (e.g. bad_alloc) is caught in the thread that raised it, since one +// escaping a std::thread would terminate MATLAB too; the return value is false if any body threw. *ran, if given, +// receives the number of threads that actually ran bodies (the caller plus those created). +template +bool run_threads(size_t nt, F&& body, size_t* ran = nullptr) { + std::atomic ok{true}; + auto safe = [&](size_t t) { try { body(t); } catch (...) { ok.store(false); } }; + std::vector pool; + size_t started = 1; + try { + pool.reserve(nt); + for (size_t t = 1; t < nt; ++t) { pool.emplace_back(safe, t); ++started; } + } catch (...) { + } + safe(0); + for (auto& th : pool) th.join(); + for (size_t t = started; t < nt; ++t) safe(t); + if (ran) *ran = started; + return ok.load(); +} + +} // namespace + +using matlab::mex::ArgumentList; +using matlab::data::Array; +using matlab::data::ArrayFactory; +using matlab::data::ArrayType; + +class MexFunction : public matlab::mex::Function { + std::shared_ptr engine = getEngine(); + ArrayFactory factory; + + [[noreturn]] static void raise(const std::string& id, const std::string& msg) { throw MexError{id, msg}; } + + // Address and byte length of a full numeric array's buffer, read through the const iterator (no unsharing). + static void span_of(const Array& a, const void*& ptr, size_t& bytes) { + ptr = nullptr; bytes = 0; + const size_t ne = a.getNumberOfElements(); + if (ne == 0) return; + switch (a.getType()) { + case ArrayType::DOUBLE: { const matlab::data::TypedArray t = a; ptr = &*t.cbegin(); bytes = ne * 8; break; } + case ArrayType::SINGLE: { const matlab::data::TypedArray t = a; ptr = &*t.cbegin(); bytes = ne * 4; break; } + default: raise("randlapack:numa_place:class", "A must be a full single or double array"); + } + } + + double scalar(const Array& a, const char* what) { + if (a.getNumberOfElements() != 1 || a.getType() != ArrayType::DOUBLE) + raise("randlapack:numa_place:arg", std::string(what) + " must be a double scalar"); + const matlab::data::TypedArray t = a; + return t[0]; + } + + std::string text(const Array& a, const char* what) { + if (a.getType() != ArrayType::CHAR) raise("randlapack:numa_place:arg", std::string(what) + " must be a char array"); + const matlab::data::CharArray c = a; + return c.toAscii(); + } + +#ifdef __linux__ + // Allowed memory nodes of this process (the cpuset of a SLURM job, all online nodes otherwise). + // Placement is only a performance helper: where the NUMA syscalls are unavailable (no CONFIG_NUMA, or a seccomp + // profile denying them) this returns no nodes and sets err, and the callers skip placement instead of failing. + std::vector allowed_nodes(std::vector& mask, int& err) { + mask.assign(kMaskWords, 0ul); + err = 0; + int mode = 0; + if (sys_get_mempolicy(&mode, mask.data(), kMaxNodes, nullptr, kMpolFMemsAllowed) != 0) { + err = errno; + return {}; + } + std::vector nodes; + for (unsigned long b = 0; b < kMaxNodes; ++b) + if (mask[b / (8 * sizeof(unsigned long))] & (1ul << (b % (8 * sizeof(unsigned long))))) nodes.push_back(static_cast(b)); + return nodes; + } + + // Node histogram of up to max_pages pages sampled evenly over [base, base + bytes). + // One page from each of m equal strata of the span, at a deterministic pseudo-random offset inside its stratum: + // evenly spaced samples could alias with spread's round robin (a stride of a whole number of 2 MiB granules times + // the node count would sample one node only). absent counts -ENOENT (never touched); errors, any other failure. + void histogram(const void* base, size_t bytes, size_t max_pages, std::vector& counts, double& absent, + double& sampled, double& errors) { + counts.clear(); absent = 0; sampled = 0; errors = 0; + if (bytes == 0) return; + const size_t ps = page_size(); + const uintptr_t lo = reinterpret_cast(base) & ~(uintptr_t)(ps - 1); + const uintptr_t hi = reinterpret_cast(base) + bytes; + const size_t npages = (hi - lo + ps - 1) / ps; + const size_t m = std::max(1, std::min(max_pages, npages)); + std::vector pages(m); + std::vector status(m, -1); + for (size_t i = 0; i < m; ++i) { + const size_t s0 = npages * i / m, s1 = npages * (i + 1) / m; + const size_t width = std::max(1, s1 - s0); + const size_t jitter = static_cast((static_cast(i) * 2654435761ull + 12345ull) % width); + pages[i] = reinterpret_cast(lo + std::min(npages - 1, s0 + jitter) * ps); + } + if (sys_move_pages(0, m, pages.data(), nullptr, status.data(), 0) != 0) { + errors = static_cast(m); sampled = static_cast(m); // query unavailable: all unknown + return; + } + for (size_t i = 0; i < m; ++i) { + if (status[i] >= 0) { + if (static_cast(status[i]) >= counts.size()) counts.resize(status[i] + 1, 0.0); + counts[status[i]] += 1.0; + } else if (status[i] == -ENOENT) { + absent += 1.0; // not faulted in yet + } else { + errors += 1.0; + } + } + sampled = static_cast(m); + } + + matlab::data::StructArray hist_struct(const std::vector& counts, double absent, double sampled, + double errors) { + auto s = factory.createStructArray({1, 1}, {"counts", "absent", "sampled", "errors"}); + s[0]["errors"] = factory.createScalar(errors); + s[0]["counts"] = factory.createArray({1, counts.size()}, counts.data(), counts.data() + counts.size()); + s[0]["absent"] = factory.createScalar(absent); + s[0]["sampled"] = factory.createScalar(sampled); + return s; + } + + void spread(ArgumentList& outputs, ArgumentList& inputs) { + if (inputs.size() != 2) raise("randlapack:numa_place:nargin", "spread takes one array"); + const void* ptr; size_t bytes; + span_of(inputs[1], ptr, bytes); + std::vector mask; + int nodes_errno = 0; + const std::vector nodes = allowed_nodes(mask, nodes_errno); + std::vector before, after; double ab = 0, aa = 0, sb = 0, sa = 0, eb = 0, ea = 0; + histogram(ptr, bytes, 4096, before, ab, sb, eb); + const auto t0 = std::chrono::steady_clock::now(); + std::atomic failed_batches{0}, not_moved{0}, unresolved{0}, absent{0}, pages_total{0}, moved_ok{0}, retry_calls{0}; + std::atomic first_errno{0}; + size_t nt = 1, nt_ran = 0; // 0 threads when the spread is skipped (one node) + // RANDLAPACK_NUMA_SPREAD_FORCE=1 runs the move_pages calls on a one-node machine too (every destination is node + // 0 there, so nothing moves): a test hook for the syscall path, which a laptop cannot otherwise reach. + const char* force = std::getenv("RANDLAPACK_NUMA_SPREAD_FORCE"); + const bool forced = force != nullptr && force[0] == '1'; + if (bytes > 0 && !nodes.empty() && (nodes.size() > 1 || forced)) { + const size_t ps = page_size(); + const uintptr_t lo = reinterpret_cast(ptr) & ~(uintptr_t)(ps - 1); + const uintptr_t hi = (reinterpret_cast(ptr) + bytes + ps - 1) & ~(uintptr_t)(ps - 1); + const uintptr_t granule = (uintptr_t)2 << 20; + const uintptr_t g0 = lo / granule, g1 = (hi + granule - 1) / granule; + const size_t ngran = static_cast(g1 - g0); + const size_t nn = nodes.size(); + // Threads own whole 2 MiB granules (a transparent huge page is never split between two concurrent + // move_pages calls), at least 32 granules (64 MiB) each; batches of 65,536 pages per call. + nt = thread_count(0, ngran / 32); + // pages = the whole span, counted up front: a worker that fails before or while submitting leaves its + // remainder in no category, so on_target + absent < pages and the placement reads as incomplete. + pages_total.store(static_cast((hi - lo) / ps)); + constexpr int kUnset = std::numeric_limits::min(); + const bool all_ran = run_threads(nt, [&](size_t t) { + const uintptr_t a0 = std::max(lo, (g0 + ngran * t / nt) * granule); + const uintptr_t a1 = std::min(hi, (g0 + ngran * (t + 1) / nt) * granule); + if (a1 <= a0) return; + const size_t tp = static_cast((a1 - a0) / ps); + constexpr size_t kBatch = 65536; + std::vector pages, rpages; + std::vector dst, status, rdst, rstatus; + try { pages.reserve(kBatch); dst.reserve(kBatch); status.resize(kBatch); } + catch (...) { failed_batches.fetch_add(1); unresolved.fetch_add(static_cast(tp)); return; } + auto note_errno = [&](int e) { int z = 0; first_errno.compare_exchange_strong(z, e); }; + for (size_t b0 = 0; b0 < tp; b0 += kBatch) { + const size_t cnt = std::min(kBatch, tp - b0); + pages.clear(); dst.clear(); + for (size_t i = 0; i < cnt; ++i) { + const uintptr_t a = a0 + (b0 + i) * ps; + pages.push_back(reinterpret_cast(a)); + dst.push_back(nodes[(a / granule) % nn]); + } + // The kernel writes a status only for the pages it got to: a positive return (pages not migrated, + // Linux >= 4.17) can come back with the rest of the list unprocessed. A sentinel marks those. + std::fill(status.begin(), status.begin() + cnt, kUnset); + long r = sys_move_pages(0, cnt, pages.data(), dst.data(), status.data(), kMpolMfMove); + if (r < 0) { + note_errno(errno); failed_batches.fetch_add(1); unresolved.fetch_add(static_cast(cnt)); + continue; + } + // One retry for every page not yet on its destination (unprocessed, or failed other than -ENOENT), + // one 2 MiB granule per call: the kernel stops a call at the first same-destination run (granule) it + // cannot migrate, so one stuck granule must not hold back the rest of the batch. + auto off_target = [&](size_t i) { + return status[i] == kUnset || (status[i] < 0 && status[i] != -ENOENT) || (status[i] >= 0 && status[i] != dst[i]); + }; + std::vector ridx; + for (size_t i = 0; i < cnt; ++i) if (off_target(i)) ridx.push_back(i); + // First ask where those pages are (nodes = NULL: a query, no LRU drain). With transparent huge pages + // a tail page typically reports -EBUSY while its head moves the whole folio, so most "off target" + // entries are already on target; only the rest are retried. + if (!ridx.empty()) { + rpages.clear(); + for (size_t i : ridx) rpages.push_back(pages[i]); + rstatus.assign(rpages.size(), kUnset); + if (sys_move_pages(0, rpages.size(), rpages.data(), nullptr, rstatus.data(), 0) == 0) { + std::vector still; + for (size_t m = 0; m < ridx.size(); ++m) { + if (rstatus[m] >= 0 || rstatus[m] == -ENOENT) status[ridx[m]] = rstatus[m]; + if (off_target(ridx[m])) still.push_back(ridx[m]); + } + ridx.swap(still); + } + } + for (size_t j = 0; j < ridx.size();) { + const uintptr_t gr = reinterpret_cast(pages[ridx[j]]) / granule; + size_t k = j; + rpages.clear(); rdst.clear(); + while (k < ridx.size() && reinterpret_cast(pages[ridx[k]]) / granule == gr) { + rpages.push_back(pages[ridx[k]]); rdst.push_back(dst[ridx[k]]); ++k; + } + rstatus.assign(rpages.size(), kUnset); + retry_calls.fetch_add(1); + r = sys_move_pages(0, rpages.size(), rpages.data(), rdst.data(), rstatus.data(), kMpolMfMove); + if (r < 0) { // the retried pages' outcome is unknown: unresolved, whatever the first call said + note_errno(errno); failed_batches.fetch_add(1); + for (size_t m = j; m < k; ++m) status[ridx[m]] = kUnset; + } else { + for (size_t m = j; m < k; ++m) { + status[ridx[m]] = rstatus[m - j]; + if (rstatus[m - j] < 0 && rstatus[m - j] != -ENOENT && rstatus[m - j] != kUnset) note_errno(-rstatus[m - j]); + } + } + j = k; + } + // The kernel's last word: where is every page still not on target? move_pages with nodes = NULL moves + // nothing and reports each page's node. This also makes the counts exact on 4.17-5.5 kernels that + // drop the status of a final run of pages already on their target. + ridx.clear(); + for (size_t i = 0; i < cnt; ++i) if (off_target(i)) ridx.push_back(i); + if (!ridx.empty()) { + rpages.clear(); + for (size_t i : ridx) rpages.push_back(pages[i]); + rstatus.assign(rpages.size(), kUnset); + if (sys_move_pages(0, rpages.size(), rpages.data(), nullptr, rstatus.data(), 0) == 0) + for (size_t m = 0; m < ridx.size(); ++m) + if (rstatus[m] >= 0 || rstatus[m] == -ENOENT || status[ridx[m]] == kUnset) status[ridx[m]] = rstatus[m]; + } + for (size_t i = 0; i < cnt; ++i) { + if (status[i] == kUnset) unresolved.fetch_add(1); + else if (status[i] == -ENOENT) absent.fetch_add(1); // never touched: no page to move + else if (status[i] < 0) { not_moved.fetch_add(1); note_errno(-status[i]); } + else if (status[i] != dst[i]) not_moved.fetch_add(1); // still on another node + else moved_ok.fetch_add(1); + } + } + }, &nt_ran); + if (!all_ran) failed_batches.fetch_add(1); // a worker threw (out of memory): its span may be unmoved + } + const double ms = std::chrono::duration(std::chrono::steady_clock::now() - t0).count(); + histogram(ptr, bytes, 4096, after, aa, sa, ea); + auto s = factory.createStructArray({1, 1}, {"supported", "nodes", "bytes", "threads", "ms", "pages", "on_target", + "failed_batches", "not_moved", "unresolved", "absent", "errno", + "retry_calls", "before", "after"}); + s[0]["retry_calls"] = factory.createScalar(static_cast(retry_calls.load())); + std::vector nd(nodes.begin(), nodes.end()); + s[0]["supported"] = factory.createScalar(true); + s[0]["nodes"] = factory.createArray({1, nd.size()}, nd.data(), nd.data() + nd.size()); + s[0]["bytes"] = factory.createScalar(static_cast(bytes)); + s[0]["threads"] = factory.createScalar(static_cast(nt_ran)); // threads that ran, not asked for + s[0]["ms"] = factory.createScalar(ms); + s[0]["pages"] = factory.createScalar(static_cast(pages_total.load())); + s[0]["on_target"] = factory.createScalar(static_cast(moved_ok.load())); + s[0]["failed_batches"] = factory.createScalar(static_cast(failed_batches.load())); + s[0]["unresolved"] = factory.createScalar(static_cast(unresolved.load())); + s[0]["not_moved"] = factory.createScalar(static_cast(not_moved.load())); + s[0]["absent"] = factory.createScalar(static_cast(absent.load())); + s[0]["errno"] = factory.createScalar(static_cast(first_errno.load() ? first_errno.load() : nodes_errno)); + s[0]["before"] = hist_struct(before, ab, sb, eb); + s[0]["after"] = hist_struct(after, aa, sa, ea); + outputs[0] = std::move(s); + } + + void query(ArgumentList& outputs, ArgumentList& inputs) { + if (inputs.size() < 2 || inputs.size() > 3) raise("randlapack:numa_place:nargin", "query takes an array and an optional page count"); + const void* ptr; size_t bytes; + span_of(inputs[1], ptr, bytes); + size_t mp = 4096; + if (inputs.size() == 3) { + const double v = scalar(inputs[2], "max_pages"); + if (!std::isfinite(v) || v < 1 || v > 1.0e9 || v != std::floor(v)) + raise("randlapack:numa_place:arg", "max_pages must be a positive integer <= 1e9"); + mp = static_cast(v); + } + std::vector c; double ab = 0, sm = 0, er = 0; + histogram(ptr, bytes, mp, c, ab, sm, er); + auto h = hist_struct(c, ab, sm, er); + // Offset of A's first byte within its page and modulo 64 bytes: MKL's rounding can depend on alignment, so a + // bit-level difference between two placements of the same matrix is read against this. + auto s = factory.createStructArray({1, 1}, {"counts", "absent", "sampled", "errors", "page_offset", "align64"}); + const uintptr_t a = reinterpret_cast(ptr); + s[0]["counts"] = h[0]["counts"]; s[0]["absent"] = h[0]["absent"]; s[0]["sampled"] = h[0]["sampled"]; + s[0]["errors"] = h[0]["errors"]; + s[0]["page_offset"] = factory.createScalar(static_cast(a % page_size())); + s[0]["align64"] = factory.createScalar(static_cast(a % 64)); + outputs[0] = std::move(s); + } + + void load(ArgumentList& outputs, ArgumentList& inputs) { + if (inputs.size() != 8) raise("randlapack:numa_place:nargin", "load takes file, offset_bytes, n_rows, n_cols, divisor, n_threads, interleave"); + const std::string file = text(inputs[1], "file"); + const double off_d = scalar(inputs[2], "offset_bytes"), nr_d = scalar(inputs[3], "n_rows"), + nc_d = scalar(inputs[4], "n_cols"), divisor = scalar(inputs[5], "divisor"), + nt_d = scalar(inputs[6], "n_threads"), il_d = scalar(inputs[7], "interleave"); + if (il_d != 0 && il_d != 1) raise("randlapack:numa_place:arg", "interleave must be 0 or 1"); + for (double v : {off_d, nr_d, nc_d, nt_d}) + if (!(v >= 0) || v != std::floor(v) || v > 9.0e15) raise("randlapack:numa_place:arg", "offset, sizes and thread count must be nonnegative integers"); + if (!std::isfinite(divisor) || divisor == 0) raise("randlapack:numa_place:arg", "divisor must be finite and nonzero"); + const size_t nr = static_cast(nr_d), nc = static_cast(nc_d); + const off_t off = static_cast(off_d); + const size_t ne = nr * nc; + if (nr != 0 && ne / nr != nc) raise("randlapack:numa_place:arg", "n_rows * n_cols overflows"); + const int fd = ::open(file.c_str(), O_RDONLY | O_CLOEXEC); + if (fd < 0) raise("randlapack:numa_place:open", "cannot open " + file + ": " + std::strerror(errno)); + struct FdGuard { int fd; ~FdGuard() { ::close(fd); } } guard{fd}; + const off_t end = ::lseek(fd, 0, SEEK_END); + if (end < 0 || static_cast(end) < off_d + 8.0 * static_cast(ne)) + raise("randlapack:numa_place:short", file + " is shorter than offset + 8*n_rows*n_cols bytes"); + const double rss0 = rss_mib(); + // createBuffer returns uninitialized storage: its pages are first touched by the pread calls below. + matlab::data::buffer_ptr_t buf = factory.createBuffer(ne); + const double rss_alloc = rss_mib(); + double* dst = buf.get(); + // Interleave policy on the untouched buffer (rounded out to whole pages; a boundary page shared with a + // neighbouring allocation is already present and simply stays where it is). + int il_errno = 0; bool il_set = false; + std::vector mask; + const std::vector nodes = allowed_nodes(mask, il_errno); // reported in stats for either placement + if (il_d == 1 && ne > 0) { + const char* force = std::getenv("RANDLAPACK_NUMA_SPREAD_FORCE"); + if (nodes.size() > 1 || (force != nullptr && force[0] == '1')) { + const size_t ps = page_size(); + const uintptr_t lo = reinterpret_cast(dst) & ~(uintptr_t)(ps - 1); + const uintptr_t hi = (reinterpret_cast(dst + ne) + ps - 1) & ~(uintptr_t)(ps - 1); + if (sys_mbind(reinterpret_cast(lo), hi - lo, kMpolInterleave, mask.data(), kMaxNodes, 0) != 0) + il_errno = errno; + else + il_set = true; + } + } + const auto t0 = std::chrono::steady_clock::now(); + // Spans of at least 32 MiB per thread. + const size_t nt = thread_count(static_cast(nt_d), (ne * 8) / (32u << 20)); + std::atomic err{0}; + std::atomic short_read{false}; + size_t nt_ran = 1; + const bool all_ran = run_threads(nt, [&](size_t t) { + const size_t lo = ne * t / nt, hi = ne * (t + 1) / nt; + size_t done = 0, want = (hi - lo) * 8; + char* p = reinterpret_cast(dst + lo); + while (done < want && err.load() == 0 && !short_read.load()) { + const size_t step = std::min(want - done, 64u << 20); + const ssize_t got = ::pread(fd, p + done, step, off + static_cast(lo * 8 + done)); + if (got < 0) { if (errno == EINTR) continue; int z = 0; err.compare_exchange_strong(z, errno); break; } + if (got == 0) { short_read = true; break; } + done += static_cast(got); + } + if (divisor != 1.0 && err.load() == 0 && !short_read.load()) + for (size_t i = lo; i < hi; ++i) dst[i] = dst[i] / divisor; + }, &nt_ran); + if (!all_ran) raise("randlapack:numa_place:read", "a reader thread failed while loading " + file); + if (err.load() != 0) raise("randlapack:numa_place:read", "read of " + file + " failed: " + std::strerror(err.load())); + if (short_read.load()) raise("randlapack:numa_place:short", "unexpected end of " + file); + const double ms = std::chrono::duration(std::chrono::steady_clock::now() - t0).count(); + const double rss_fill = rss_mib(); + outputs[0] = factory.createArrayFromBuffer({nr, nc}, std::move(buf)); + if (outputs.size() >= 2) { + auto s = factory.createStructArray({1, 1}, {"threads", "ms", "rss_before_mib", "rss_after_alloc_mib", "rss_after_fill_mib", + "interleave_set", "interleave_errno", "nodes"}); + std::vector nd(nodes.begin(), nodes.end()); + s[0]["interleave_set"] = factory.createScalar(il_set); + s[0]["interleave_errno"] = factory.createScalar(static_cast(il_errno)); + s[0]["nodes"] = factory.createArray({1, nd.size()}, nd.data(), nd.data() + nd.size()); + s[0]["threads"] = factory.createScalar(static_cast(nt_ran)); // threads that ran + s[0]["ms"] = factory.createScalar(ms); + s[0]["rss_before_mib"] = factory.createScalar(rss0); + s[0]["rss_after_alloc_mib"] = factory.createScalar(rss_alloc); + s[0]["rss_after_fill_mib"] = factory.createScalar(rss_fill); + outputs[1] = std::move(s); + } + } +#endif + +public: + void operator()(ArgumentList outputs, ArgumentList inputs) { + try { + if (inputs.size() < 1) raise("randlapack:numa_place:nargin", "usage: numa_place_mex('spread'|'query'|'load', ...)"); + const std::string cmd = text(inputs[0], "command"); +#ifdef __linux__ + if (cmd != "spread" && cmd != "query" && cmd != "load") + raise("randlapack:numa_place:command", "unknown command '" + cmd + "'"); + const size_t max_out = (cmd == "load") ? 2 : 1; + if (outputs.size() < 1 || outputs.size() > max_out) + raise("randlapack:numa_place:nargout", "'" + cmd + "' returns " + (max_out == 1 ? "one output" : "one or two outputs")); + if (cmd == "spread") spread(outputs, inputs); + else if (cmd == "query") query(outputs, inputs); + else if (cmd == "load") load(outputs, inputs); + else raise("randlapack:numa_place:command", "unknown command '" + cmd + "'"); +#else + raise("randlapack:numa_place:unsupported", "numa_place_mex is Linux only"); +#endif + } catch (const MexError& e) { + engine->feval(u"error", 0, std::vector{factory.createCharArray(e.id), factory.createCharArray(e.msg)}); + } catch (const std::exception& e) { + engine->feval(u"error", 0, std::vector{factory.createCharArray("randlapack:numa_place:StdError"), + factory.createCharArray(e.what())}); + } + } +}; diff --git a/matlab/test_trace_effdim_parameter.m b/matlab/test_trace_effdim_parameter.m new file mode 100644 index 0000000..7f43cae --- /dev/null +++ b/matlab/test_trace_effdim_parameter.m @@ -0,0 +1,22 @@ +function test_trace_effdim_parameter +% Public and expert interfaces must estimate the SAME shifted function. + lam=logspace(-2,2,80)'; + for shift=[.1 1 10] + for method={'scalar','block','scalar_auto','block_auto'} + [e,t]=randlapack.trace_estimate(lam,'Method',method{1},'Func','effdim', ... + 'Lambda',shift,'Depth',10,'MaxMatvecs',200,'Seed',91); + if endsWith(method{1},'_auto') + expected=randlapack.fun_nystrom_pp(lam,1,1,'LFAType',method{1}, ... + 'Func','effdim','PolyLambda',shift,'Budget',200,'AutoEps',1e-3, ... + 'CapRankFraction',.5,'SketchSeed',91); + else + expected=randlapack.fun_nystrom_pp(lam,t.rank,t.samples,'LFAType',[method{1} '_qfa'], ... + 'Func','effdim','PolyLambda',shift,'Depth',10,'Adaptive',0,'Q',1,'SketchSeed',91); + end + assert(abs(e-expected)<1e-12*max(1,abs(expected))); + end + end + f=@(x) x./(x+1);g=@(x) x./(x+10); + assert(abs(sum(f(lam))-sum(g(lam)))>1,'Test must distinguish the two regularizations.'); + fprintf('Effective-dimension parameter forwarding passed for all four methods.\n'); +end diff --git a/matlab/test_trace_estimate.m b/matlab/test_trace_estimate.m new file mode 100644 index 0000000..511bed5 --- /dev/null +++ b/matlab/test_trace_estimate.m @@ -0,0 +1,56 @@ +function test_trace_estimate +% Public contract, finite caps, and representation-independent probe law. + for precision={'double','single'} + lam=cast(logspace(0,-4,96)',precision{1}); + for method={'scalar','block','scalar_auto','block_auto'} + args={'Method',method{1},'TargetRelError',1e-3, ... + 'MaxMatvecs',400,'Depth',20,'Seed',77}; + [e,t]=randlapack.trace_estimate(lam,args{:}); + [ed,td]=randlapack.trace_estimate(diag(lam),args{:}); + assert(isfinite(e) && t.matvecs<=400); + assert(t.matvecs==td.matvecs); + assert(abs(e-ed)<100*eps(precision{1})*max(abs(e),1)); + assert(~t.total_error_certified); + assert(t.oracle_mean_depth <= t.depth_limit + 1e-9); + assert(abs(t.oracle_mean_depth*t.samples + t.rank+t.probe_mv-t.matvecs)<1e-8); + assert(t.oracle_mean_depth==td.oracle_mean_depth); + end + for method={'scalar_auto','block_auto'} + for fraction=[1 .5 .25 .125] + for qfraction=[1 .25] + [e,~,~,t]=randlapack.fun_nystrom_pp(lam,1,1, ... + 'LFAType',method{1},'Budget',180,'AutoEps',1e-3, ... + 'CapRankFraction',fraction,'QuadratureFraction',qfraction); + assert(isfinite(e) && t.probe_mv+t.auto_k+t.oracle_mv<=180); + end + end + for cap=[1 16 17 32 100] + try + [e,t]=randlapack.trace_estimate(lam,'Method',method{1}, ... + 'MaxMatvecs',cap,'TargetRelError',1e-3); + assert(isfinite(e) && t.matvecs<=cap && cap>=17); + catch ME + assert(cap<17 && contains(ME.identifier,'infeasibleBudget'),ME.message); + end + end + end + end + for method={'scalar_auto','block_auto'} + try + randlapack.fun_nystrom_pp(lam,1,1,'LFAType',method{1}, ... + 'Func','entropy','Budget',100); + error('test_trace_estimate:acceptedEntropy','Unsupported certificate accepted.'); + catch ME + assert(contains(ME.identifier,'unsupportedAutoFunction'),ME.message); + end + end + % Canonicalization must happen before the function-specific Lambda branch. + A=logspace(0,-3,48)'; + for method={'scalar','block_auto'} + args={'Method',method{1},'Lambda',.25,'Depth',8,'MaxMatvecs',96,'Seed',7}; + [a,ta]=randlapack.trace_estimate(A,args{:},'Func','effdim'); + [b,tb]=randlapack.trace_estimate(A,args{:},'Func','EFFDIM'); + assert(a==b && ta.matvecs==tb.matvecs && strcmp(tb.function_name,'effdim')); + end + fprintf('Four-method target/cap and dense/diagonal tests passed in both precisions.\n'); +end diff --git a/matlab/test_trace_sketch_density.m b/matlab/test_trace_sketch_density.m new file mode 100644 index 0000000..51a4ec8 --- /dev/null +++ b/matlab/test_trace_sketch_density.m @@ -0,0 +1,30 @@ +function test_trace_sketch_density +% Verify that the public knob changes the actual head, including auto modes. + A=logspace(0,-4,128)'; + for method={'scalar','block','scalar_auto','block_auto'} + estimates=[]; + for vnz=[8 16 32 0] + [e,t]=randlapack.trace_estimate(A,'Method',method{1},'Depth',8, ... + 'MaxMatvecs',256,'Seed',701,'VecNnz',vnz); + args={'Func','sqrt','VecNnz',vnz,'SketchSeed',701}; + if endsWith(method{1},'_auto') + ref=randlapack.fun_nystrom_pp(A,1,1,args{:},'LFAType',method{1}, ... + 'Budget',256,'AutoEps',1e-3,'CapRankFraction',.5); + else + ref=randlapack.fun_nystrom_pp(A,t.rank,t.samples,args{:}, ... + 'LFAType',[method{1} '_qfa'],'Q',1,'Depth',8,'Adaptive',0); + end + assert(abs(e-ref)<1e-12*max(1,abs(ref))); + assert(t.vec_nnz_effective>=1 && t.vec_nnz_effective<=t.rank); + estimates(end+1)=e; %#ok + end + assert(max(estimates)-min(estimates)>1e-8,'Density control did not affect this head.'); + end + try + randlapack.fun_nystrom_pp(A,8,8,'Sketch','gaussian'); + error('test_trace_sketch_density:acceptedUnsupported','Unsupported sketch was accepted.'); + catch ME + assert(strcmp(ME.identifier,'randlapack:fun_nystrom_pp:unsupportedSketch'),ME.message); + end + fprintf('Sketch-density forwarding/effect and unsupported-sketch rejection passed.\n'); +end diff --git a/matlab/tests/test_fun_nystrom_pp.m b/matlab/tests/test_fun_nystrom_pp.m new file mode 100644 index 0000000..97e67c4 --- /dev/null +++ b/matlab/tests/test_fun_nystrom_pp.m @@ -0,0 +1,597 @@ +function test_fun_nystrom_pp() +%test_fun_nystrom_pp Correctness checks for randlapack.fun_nystrom_pp. +% +% Seed-driven API; the Phase-1 sketch is always the kernel-internal SASO. +% For an SPD A with known spectrum, tr(f(A)) computed by funNystrom++ should +% match the analytic trace within Hutchinson sampling error. +% Checks, for both single and double: +% * all six LFATypes ('exact','scalar','scalar_qfa','block','block_qfa', +% 'auto') round-trip within tolerance of the analytic trace +% * the times struct carries the campaign's cost accounting: d_used (NaN +% for 'exact', which never runs a Lanczos recurrence); +% oracle_mv (= s*d_used for block_qfa, = sum of per-column depths for +% scalar_qfa); auto_k/auto_s/probe_mv with the budget closure +% probe_mv + q*auto_k + oracle_mv <= Budget (q = 1) for 'auto'; +% probe_converged/phase2_certified present; path-specific fields are NaN +% off their path +% * block_qfa + Adaptive=1 (Radau-certified) with f = log: tr_U >= tr_L +% bracket, certified flag set; RadauReturn = 1 (midpoint) runs +% * block_qfa + Adaptive=2 (legacy window) runs and never certifies +% * k == n => exact (rank-n Nystrom captures A; Phase 2 skipped) ~ machine eps +% * LFAType 'auto' with an infeasible Budget raises +% 'randlapack:fun_nystrom_pp:infeasibleBudget' +% * LFAType 'adaptive' (eps-targeted, no upfront Budget): success path with +% a loose AutoEps round-trips within tolerance, telemetry fields (auto_k, +% auto_s, d_used, probe_mv, oracle_mv, tr_U, tr_L, certified, +% probe_converged, phase2_certified) populated; a tiny AdaptiveMatvecCap +% raises 'randlapack:fun_nystrom_pp:infeasibleBudget' +% * VecNnz = 0 ("auto" ~log k) is accepted +% * explicit Omega2 works; an empty (n x 0) Omega2 with k < n raises +% * a non-symmetric A raises 'randlapack:fun_nystrom_pp:NotSymmetric' by +% default; 'SkipSymCheck', true bypasses the check +% * a warning ('randlapack:fun_nystrom_pp:ignored_knob') fires from the .m +% wrapper (inputParser.UsingDefaults) when 'Q'/'Adaptive*' are +% EXPLICITLY passed together with LFAType 'auto'/'adaptive' (which +% ignore them) -- including when the explicit value equals that knob's +% own default (UsingDefaults tracks "was it passed", not "does the +% value differ") -- and stays silent when the knob is left unset +% * LFAType 'auto'/'adaptive' still work when given an explicit Omega2 +% matrix (the marshal is skipped for these tiers, but the call must not +% break) +% * LFAType 'auto' with Budget < 1 raises 'randlapack:fun_nystrom_pp:Budget' +% (the id unified between the .m pre-check and the MEX's own guard) +% * ignored-knob warning coverage is UNIFORM across LFATypes, not just +% 'auto'/'adaptive': 'scalar_qfa' explicitly given 'Reorth' warns (it is +% intrinsically no-reorth); 'exact' explicitly given 'Depth' warns (no +% Lanczos recurrence runs) +% * 'Reorth' outside {0, 1} is rejected by inputParser validation +% ('MATLAB:InputParser:ArgumentFailedValidation') +% * a sparse double A is accepted and agrees with the dense call at the +% same seed to rounding; 'exact' on a sparse A is refused with the +% specific id 'randlapack:fun_nystrom_pp_mex:exact_sparse' +% +% Run: test_fun_nystrom_pp (errors via assert on failure; for CI). + + n = 200; + k = 60; + s = 40; + budget = 600; + for T = {'double', 'single'} + cls = T{1}; + rng(7); + [Q, ~] = qr(randn(n), 0); + lam = logspace(0, -3, n)'; + A = cast(Q * (lam .* Q'), cls); A = (A + A') / 2; + true_sqrt = sum(sqrt(lam)); + true_log = sum(log(lam + 1)); + + % --- 1. all six LFATypes round-trip against the sqrt fixture --- + lfa_types = {'exact', 'scalar', 'scalar_qfa', 'block', 'block_qfa', 'auto'}; + depth_for = struct('exact', 0, 'scalar', 60, 'scalar_qfa', 60, ... + 'block', 20, 'block_qfa', 20, 'auto', 0); + for it = 1:numel(lfa_types) + lfa = lfa_types{it}; + if strcmp(lfa, 'auto') + % 'Q' is ignored by 'auto' (D-I2); omit it here rather than + % pass a non-default value that would now warn. + extra = {'Budget', budget, 'AutoEps', 1e-2}; + [est, ~, ~, times] = randlapack.fun_nystrom_pp(A, k, s, ... + 'Func', 'sqrt', 'LFAType', lfa, 'SketchSeed', 11, extra{:}); + else + if depth_for.(lfa) > 0 + extra = {'Depth', depth_for.(lfa)}; + else + extra = {}; + end + [est, ~, ~, times] = randlapack.fun_nystrom_pp(A, k, s, ... + 'Func', 'sqrt', 'Q', 1, 'LFAType', lfa, 'SketchSeed', 11, extra{:}); + end + rel = abs(est - true_sqrt) / true_sqrt; + assert(isfinite(est) && rel < 5e-2, ... + '[%s/%s] rel_err = %.2e exceeds 5e-2', cls, lfa, rel); + + % d_used is the campaign's depth accounting; present for every + % type, but NaN for 'exact' (no Lanczos recurrence runs, so the + % leftover default-depth input is not a meaningful cost figure). + if strcmp(lfa, 'exact') + assert(isfield(times, 'd_used') && isnan(times.d_used), ... + '[%s/%s] d_used should be NaN for exact, got %g', ... + cls, lfa, times.d_used); + else + assert(isfield(times, 'd_used') && times.d_used >= 1, ... + '[%s/%s] d_used missing or < 1', cls, lfa); + end + + switch lfa + case 'scalar_qfa' + % Fixed depth: every column runs d steps, so the sum of + % per-column depths is exactly s * d_used. + assert(times.oracle_mv == s * times.d_used, ... + '[%s/%s] oracle_mv %g ~= s*d_used %g', ... + cls, lfa, times.oracle_mv, s * times.d_used); + case 'block_qfa' + % The class reports matvecs = s * d_used always. + assert(times.oracle_mv == s * times.d_used, ... + '[%s/%s] oracle_mv %g ~= s*d_used %g', ... + cls, lfa, times.oracle_mv, s * times.d_used); + assert(~isnan(times.tr_U) && ~isnan(times.tr_L) && ... + ~isnan(times.certified), ... + '[%s/%s] tr_U/tr_L/certified must be real for block_qfa', cls, lfa); + case 'auto' + assert(times.auto_k >= 1 && times.auto_s >= 1 && times.probe_mv >= 1, ... + '[%s/%s] auto_k/auto_s/probe_mv accounting missing', cls, lfa); + spend = times.probe_mv + 1 * times.auto_k + times.oracle_mv; % q = 1 + assert(spend <= budget, ... + '[%s/%s] budget closure violated: %g > %g', cls, lfa, spend, budget); + assert(~isnan(times.probe_converged) && ~isnan(times.phase2_certified), ... + '[%s/%s] probe_converged/phase2_certified must be real for auto', cls, lfa); + assert(~isnan(times.probe_ms) && times.probe_ms >= 0, ... + '[%s/%s] probe_ms missing for auto', cls, lfa); + otherwise + % exact/scalar/block report no oracle matvec count (0). + assert(times.oracle_mv == 0, ... + '[%s/%s] oracle_mv should be 0, got %g', cls, lfa, times.oracle_mv); + end + % Path-specific fields must be NaN off their path. + if ~strcmp(lfa, 'block_qfa') + assert(isnan(times.tr_U) && isnan(times.tr_L) && isnan(times.certified), ... + '[%s/%s] block_qfa-only fields must be NaN here', cls, lfa); + end + if ~strcmp(lfa, 'auto') + assert(isnan(times.probe_ms) && isnan(times.probe_converged) && ... + isnan(times.phase2_certified), ... + '[%s/%s] auto-only fields must be NaN here', cls, lfa); + end + end + + % --- 2. k == n: exact rank-n Nystrom, Phase 2 skipped -> near machine precision --- + est_full = randlapack.fun_nystrom_pp(A, n, 0, 'Func', 'sqrt', 'Q', 1, ... + 'LFAType', 'exact'); + rel_full = abs(est_full - true_sqrt) / true_sqrt; + tol_full = strcmp(cls, 'single') * 1e-3 + strcmp(cls, 'double') * 1e-8; + assert(rel_full < tol_full, ... + '[%s] k==n rel_err = %.2e exceeds %.0e (Phase-2-skip broken?)', ... + cls, rel_full, tol_full); + + % --- 3. block_qfa Radau-certified (Adaptive = 1) with f = log --- + % log is operator monotone, so Gauss (tr_U) and Radau (tr_L) bracket + % the true quadratic form from opposite sides: tr_U >= tr_L. + [est_bq, ~, ~, tq] = randlapack.fun_nystrom_pp(A, k, s, 'Func', 'log', ... + 'Q', 1, 'LFAType', 'block_qfa', 'Adaptive', 1, 'Depth', 60, ... + 'SketchSeed', 11); + rel_bq = abs(est_bq - true_log) / true_log; + assert(isfinite(est_bq) && rel_bq < 5e-2, ... + '[%s] block_qfa Radau rel_err = %.2e exceeds 5e-2', cls, rel_bq); + assert(tq.certified == 1, '[%s] Radau bracket did not certify', cls); + assert(tq.tr_U >= tq.tr_L - 1e-6 * abs(tq.tr_U), ... + '[%s] bracket inverted: tr_U = %.6e < tr_L = %.6e', cls, tq.tr_U, tq.tr_L); + assert(tq.d_used <= 60 && tq.oracle_mv == s * tq.d_used, ... + '[%s] certified block_qfa depth/matvec accounting wrong', cls); + + % RadauReturn = 1 (midpoint) on the same call: runs and stays accurate. + est_mid = randlapack.fun_nystrom_pp(A, k, s, 'Func', 'log', ... + 'Q', 1, 'LFAType', 'block_qfa', 'Adaptive', 1, 'Depth', 60, ... + 'RadauReturn', 1, 'SketchSeed', 11); + rel_mid = abs(est_mid - true_log) / true_log; + assert(isfinite(est_mid) && rel_mid < 5e-2, ... + '[%s] block_qfa midpoint rel_err = %.2e exceeds 5e-2', cls, rel_mid); + + % --- 4. block_qfa legacy window (Adaptive = 2): runs, never certifies --- + [est_w, ~, ~, tw] = randlapack.fun_nystrom_pp(A, k, s, 'Func', 'sqrt', ... + 'Q', 1, 'LFAType', 'block_qfa', 'Adaptive', 2, 'Depth', 60, ... + 'AdaptiveDelay', 3, 'AdaptiveMin', 4, 'SketchSeed', 11); + rel_w = abs(est_w - true_sqrt) / true_sqrt; + assert(isfinite(est_w) && rel_w < 5e-2, ... + '[%s] block_qfa window rel_err = %.2e exceeds 5e-2', cls, rel_w); + assert(tw.certified == 0, '[%s] window rule must not report certified', cls); + assert(tw.d_used <= 60, '[%s] window d_used exceeds cap', cls); + + % --- 5. scalar_qfa adaptive: oracle_mv = sum of per-column certified + % depths, bounded by [d_used, s * d_used] --- + [~, ~, ~, tsq] = randlapack.fun_nystrom_pp(A, k, s, 'Func', 'sqrt', ... + 'Q', 1, 'LFAType', 'scalar_qfa', 'Adaptive', 1, 'Depth', 80, ... + 'SketchSeed', 11); + assert(tsq.oracle_mv >= tsq.d_used && tsq.oracle_mv <= s * tsq.d_used, ... + '[%s] adaptive scalar_qfa oracle_mv %g outside [%g, %g]', ... + cls, tsq.oracle_mv, tsq.d_used, s * tsq.d_used); + % (I1 fix) certified must now be a real 0/1 for scalar_qfa+Adaptive=1 + % (scalar_qfa.all_certified marshalled through), and stay NaN for the + % fixed-depth (Adaptive=0) path since no certificate was ever checked. + assert(~isnan(tsq.certified) && any(tsq.certified == [0 1]), ... + '[%s] adaptive scalar_qfa certified must be logical (0/1), got %g', ... + cls, tsq.certified); + [~, ~, ~, tsq_fixed] = randlapack.fun_nystrom_pp(A, k, s, 'Func', 'sqrt', ... + 'Q', 1, 'LFAType', 'scalar_qfa', 'Adaptive', 0, 'Depth', 60, ... + 'SketchSeed', 11); + assert(isnan(tsq_fixed.certified), ... + '[%s] fixed-depth scalar_qfa certified should stay NaN, got %g', ... + cls, tsq_fixed.certified); + + % --- 6. infeasible auto Budget raises the dedicated error id --- + got_id = ''; + try + randlapack.fun_nystrom_pp(A, k, s, 'Func', 'sqrt', ... + 'LFAType', 'auto', 'Budget', 8); + catch err + got_id = err.identifier; + end + assert(strcmp(got_id, 'randlapack:fun_nystrom_pp:infeasibleBudget'), ... + '[%s] Budget=8 raised ''%s'', expected infeasibleBudget', cls, got_id); + + % --- 7. VecNnz = 0 ("auto" ~log k) is accepted --- + est_v0 = randlapack.fun_nystrom_pp(A, k, s, 'Func', 'sqrt', 'Q', 1, ... + 'LFAType', 'exact', 'VecNnz', 0, 'SketchSeed', 11); + rel_v0 = abs(est_v0 - true_sqrt) / true_sqrt; + assert(isfinite(est_v0) && rel_v0 < 5e-2, ... + '[%s] VecNnz=0 rel_err = %.2e exceeds 5e-2', cls, rel_v0); + + % --- 8. explicit Omega2 path works; empty Omega2 with k < n raises --- + rng(13); + Om2 = cast(randn(n, s), cls); + est_om = randlapack.fun_nystrom_pp(A, k, Om2, 'Func', 'sqrt', 'Q', 1, ... + 'LFAType', 'exact', 'SketchSeed', 11); + rel_om = abs(est_om - true_sqrt) / true_sqrt; + assert(isfinite(est_om) && rel_om < 5e-2, ... + '[%s] explicit-Omega2 rel_err = %.2e exceeds 5e-2', cls, rel_om); + got_id = ''; + try + randlapack.fun_nystrom_pp(A, k, zeros(n, 0, cls), 'Func', 'sqrt', ... + 'Q', 1, 'LFAType', 'exact'); + catch err + got_id = err.identifier; + end + assert(strcmp(got_id, 'randlapack:fun_nystrom_pp:Omega2Empty'), ... + '[%s] n x 0 Omega2 raised ''%s'', expected Omega2Empty', cls, got_id); + + % --- 9. LFAType 'adaptive': eps-targeted tier, no upfront Budget. + % Positional k/s are placeholders (ignored, like 'auto'); the + % driver picks them itself from a certified depth probe. --- + [est_ad, ~, ~, tad] = randlapack.fun_nystrom_pp(A, 1, 1, 'Func', 'sqrt', ... + 'LFAType', 'adaptive', 'AutoEps', 1e-1, 'SketchSeed', 11); + rel_ad = abs(est_ad - true_sqrt) / true_sqrt; + assert(isfinite(est_ad) && rel_ad < 5e-2, ... + '[%s] adaptive rel_err = %.2e exceeds 5e-2', cls, rel_ad); + assert(isfield(tad, 'd_used') && tad.d_used >= 1, ... + '[%s] adaptive d_used missing or < 1', cls); + assert(tad.oracle_mv > 0, ... + '[%s] adaptive oracle_mv should be > 0, got %g', cls, tad.oracle_mv); + assert(tad.auto_k > 0 && tad.auto_s > 0, ... + '[%s] adaptive auto_k/auto_s accounting missing (k=%g, s=%g)', ... + cls, tad.auto_k, tad.auto_s); + assert(~isnan(tad.probe_mv) && tad.probe_mv > 0, ... + '[%s] adaptive probe_mv missing or non-positive', cls); + assert(~isnan(tad.tr_U) && isfinite(tad.tr_U) && ... + ~isnan(tad.tr_L) && isfinite(tad.tr_L), ... + '[%s] adaptive tr_U/tr_L must be finite, got tr_U=%g tr_L=%g', ... + cls, tad.tr_U, tad.tr_L); + assert(~isnan(tad.certified) && any(tad.certified == [0 1]), ... + '[%s] adaptive certified must be logical (0/1), got %g', cls, tad.certified); + assert(~isnan(tad.probe_converged) && any(tad.probe_converged == [0 1]), ... + '[%s] adaptive probe_converged must be logical (0/1)', cls); + assert(~isnan(tad.phase2_certified) && any(tad.phase2_certified == [0 1]), ... + '[%s] adaptive phase2_certified must be logical (0/1)', cls); + + % --- 10. LFAType 'adaptive' with a tiny AdaptiveMatvecCap raises the + % SAME error id as the 'auto' tier's infeasible Budget --- + got_id = ''; + try + randlapack.fun_nystrom_pp(A, 1, 1, 'Func', 'sqrt', ... + 'LFAType', 'adaptive', 'AutoEps', 1e-1, 'AdaptiveMatvecCap', 5); + catch err + got_id = err.identifier; + end + assert(strcmp(got_id, 'randlapack:fun_nystrom_pp:infeasibleBudget'), ... + '[%s] AdaptiveMatvecCap=5 raised ''%s'', expected infeasibleBudget', ... + cls, got_id); + + % --- 11. Symmetry validation (D-I1): default-on, tolerance-gated, + % skippable. The MEX only reads the upper triangle and + % mirrors it into the lower; breaking A's upper/lower + % agreement by a modest absolute amount (well past the + % sqrt(eps) tolerance, but small next to A's own O(1) scale + % so the perturbed matrix is still well-conditioned enough + % for 'exact' to run cleanly under SkipSymCheck) must raise + % before ever reaching the MEX; 'SkipSymCheck', true bypasses + % it. SketchSeed pinned to 11 throughout, matching every + % other randomized call in this file (a couple of other + % seeds are known to hit an unrelated rare SASO + % rank-deficiency failure in NystromEVD at this k/vec_nnz). --- + A_asym = A; + A_asym(1, end) = A_asym(1, end) + 1e-2 * max(abs(A(:))); + got_id = ''; + try + randlapack.fun_nystrom_pp(A_asym, k, s, 'Func', 'sqrt', 'Q', 1, ... + 'LFAType', 'exact', 'SketchSeed', 11); + catch err + got_id = err.identifier; + end + assert(strcmp(got_id, 'randlapack:fun_nystrom_pp:NotSymmetric'), ... + '[%s] asymmetric A raised ''%s'', expected NotSymmetric', cls, got_id); + + est_skip = randlapack.fun_nystrom_pp(A_asym, k, s, 'Func', 'sqrt', ... + 'Q', 1, 'LFAType', 'exact', 'SketchSeed', 11, 'SkipSymCheck', true); + assert(isfinite(est_skip), ... + '[%s] SkipSymCheck=true run did not produce a finite estimate', cls); + + % The fixture A used throughout this test is explicitly symmetrized + % ((A+A')/2 above) and must NOT trip the default check. + est_sym_ok = randlapack.fun_nystrom_pp(A, k, s, 'Func', 'sqrt', ... + 'Q', 1, 'LFAType', 'exact', 'SketchSeed', 11); + assert(isfinite(est_sym_ok), ... + '[%s] symmetric fixture A incorrectly failed the symmetry check', cls); + + % --- 12. Ignored-knob warning (D-I2, fix batch 2-M): moved to the .m + % wrapper, driven by inputParser.UsingDefaults. 'Q'/ + % 'Adaptive*' EXPLICITLY passed together with LFAType + % 'auto'/'adaptive' warns; the knob left UNSET stays silent. --- + lastwarn(''); + [~] = randlapack.fun_nystrom_pp(A, k, s, 'Func', 'sqrt', 'SketchSeed', 11, ... + 'LFAType', 'auto', 'Budget', budget, 'AutoEps', 1e-2, 'Q', 3); + [~, warn_id] = lastwarn(); + assert(strcmp(warn_id, 'randlapack:fun_nystrom_pp:ignored_knob'), ... + '[%s] Q=3 with LFAType=auto did not warn ignored_knob (got ''%s'')', ... + cls, warn_id); + + lastwarn(''); + [~] = randlapack.fun_nystrom_pp(A, k, s, 'Func', 'sqrt', 'SketchSeed', 11, ... + 'LFAType', 'auto', 'Budget', budget, 'AutoEps', 1e-2); + [~, warn_id] = lastwarn(); + assert(isempty(warn_id), ... + '[%s] Q left unset on an auto call unexpectedly warned ''%s''', cls, warn_id); + + % UsingDefaults tracks "was it passed", not "does the value differ": + % explicitly repeating 'Q''s own default (2) must still warn here, + % unlike the old MEX-side value-comparison check it replaced (which + % could not tell this apart from "never touched"). + lastwarn(''); + [~] = randlapack.fun_nystrom_pp(A, k, s, 'Func', 'sqrt', 'SketchSeed', 11, ... + 'LFAType', 'auto', 'Budget', budget, 'AutoEps', 1e-2, 'Q', 2); + [~, warn_id] = lastwarn(); + assert(strcmp(warn_id, 'randlapack:fun_nystrom_pp:ignored_knob'), ... + '[%s] explicit Q=2 (its own default) on an auto call did not warn (got ''%s'')', ... + cls, warn_id); + + lastwarn(''); + [~] = randlapack.fun_nystrom_pp(A, 1, 1, 'Func', 'sqrt', 'SketchSeed', 11, ... + 'LFAType', 'adaptive', 'AutoEps', 1e-1, 'AdaptiveDelay', 5); + [~, warn_id] = lastwarn(); + assert(strcmp(warn_id, 'randlapack:fun_nystrom_pp:ignored_knob'), ... + '[%s] AdaptiveDelay=5 with LFAType=adaptive did not warn ignored_knob (got ''%s'')', ... + cls, warn_id); + + lastwarn(''); + [~] = randlapack.fun_nystrom_pp(A, 1, 1, 'Func', 'sqrt', 'SketchSeed', 11, ... + 'LFAType', 'adaptive', 'AutoEps', 1e-1); + [~, warn_id] = lastwarn(); + assert(isempty(warn_id), ... + '[%s] default-knob adaptive call unexpectedly warned ''%s''', cls, warn_id); + + % --- 12b. Uniform warning coverage (M5): the SAME ignored-knob + % mechanism now covers every LFAType, not just auto/adaptive. + % 'scalar_qfa' is intrinsically no-reorth; an explicit + % 'Reorth' warns. 'exact' never runs a Lanczos recurrence; + % an explicit 'Depth' warns. --- + lastwarn(''); + [~] = randlapack.fun_nystrom_pp(A, k, s, 'Func', 'sqrt', 'Q', 1, ... + 'SketchSeed', 11, 'LFAType', 'scalar_qfa', 'Depth', 60, 'Reorth', 0); + [~, warn_id] = lastwarn(); + assert(strcmp(warn_id, 'randlapack:fun_nystrom_pp:ignored_knob'), ... + '[%s] explicit Reorth with LFAType=scalar_qfa did not warn (got ''%s'')', ... + cls, warn_id); + + lastwarn(''); + [~] = randlapack.fun_nystrom_pp(A, k, s, 'Func', 'sqrt', 'Q', 1, ... + 'SketchSeed', 11, 'LFAType', 'exact', 'Depth', 60); + [~, warn_id] = lastwarn(); + assert(strcmp(warn_id, 'randlapack:fun_nystrom_pp:ignored_knob'), ... + '[%s] explicit Depth with LFAType=exact did not warn (got ''%s'')', ... + cls, warn_id); + + % --- 13. LFAType 'auto'/'adaptive' with an explicit Omega2 matrix + % (D-I3: the MEX skips marshalling it for these two tiers, + % since neither driver.call overload reads it: confirm the + % call still succeeds and returns a correct, fully-populated + % estimate rather than silently breaking). --- + rng(29); + Om2_probe = cast(randn(n, s), cls); + [est_auto_om, ~, ~, times_auto_om] = randlapack.fun_nystrom_pp(A, k, Om2_probe, ... + 'Func', 'sqrt', 'SketchSeed', 11, 'LFAType', 'auto', 'Budget', budget, 'AutoEps', 1e-2); + rel_auto_om = abs(est_auto_om - true_sqrt) / true_sqrt; + assert(isfinite(est_auto_om) && rel_auto_om < 5e-2, ... + '[%s] auto with explicit Omega2 (marshal-skip path) rel_err = %.2e exceeds 5e-2', ... + cls, rel_auto_om); + assert(times_auto_om.auto_k >= 1 && times_auto_om.auto_s >= 1, ... + '[%s] auto with explicit Omega2 telemetry missing (marshal-skip path)', cls); + + [est_ad_om, ~, ~, times_ad_om] = randlapack.fun_nystrom_pp(A, k, Om2_probe, ... + 'Func', 'sqrt', 'SketchSeed', 11, 'LFAType', 'adaptive', 'AutoEps', 1e-1); + rel_ad_om = abs(est_ad_om - true_sqrt) / true_sqrt; + assert(isfinite(est_ad_om) && rel_ad_om < 5e-2, ... + '[%s] adaptive with explicit Omega2 (marshal-skip path) rel_err = %.2e exceeds 5e-2', ... + cls, rel_ad_om); + assert(times_ad_om.auto_k > 0 && times_ad_om.auto_s > 0, ... + '[%s] adaptive with explicit Omega2 telemetry missing (marshal-skip path)', cls); + + % --- 14. Error-id unification (M1): LFAType 'auto' with Budget < 1 + % raises 'randlapack:fun_nystrom_pp:Budget' (previously the + % .m pre-check and the MEX's own guard used two different, + % undocumented ids; a direct fun_nystrom_pp_mex(...) call + % now raises the identical id this wrapper does). --- + got_id = ''; + try + randlapack.fun_nystrom_pp(A, k, s, 'Func', 'sqrt', 'Q', 1, ... + 'LFAType', 'auto', 'Budget', 0, 'AutoEps', 1e-2); + catch err + got_id = err.identifier; + end + assert(strcmp(got_id, 'randlapack:fun_nystrom_pp:Budget'), ... + '[%s] Budget=0 raised ''%s'', expected the unified Budget id', cls, got_id); + + % --- 15. Sparse input: a sparse double A (both triangles stored) + % is accepted and must agree with the dense call at the same + % seed to rounding; 'exact' needs a dense eigensolve and is + % refused with its specific id. MATLAB's sparse() only + % supports double storage, so this only runs for + % cls == 'double'. --- + if strcmp(cls, 'double') + got_id = ''; + try + randlapack.fun_nystrom_pp(sparse(A), k, s, 'Func', 'sqrt', ... + 'LFAType', 'exact'); + catch err + got_id = err.identifier; + end + assert(strcmp(got_id, 'randlapack:fun_nystrom_pp_mex:exact_sparse'), ... + '[%s] exact on sparse A raised ''%s'', expected randlapack:fun_nystrom_pp_mex:exact_sparse', ... + cls, got_id); + e_sparse = randlapack.fun_nystrom_pp(sparse(A), k, s, 'Func', 'sqrt', ... + 'LFAType', 'block_qfa', 'Depth', 20, 'Q', 1, 'SketchSeed', 7); + e_dense = randlapack.fun_nystrom_pp(A, k, s, 'Func', 'sqrt', ... + 'LFAType', 'block_qfa', 'Depth', 20, 'Q', 1, 'SketchSeed', 7); + assert(abs(e_sparse - e_dense) <= 1e-10 * abs(e_dense), ... + '[%s] sparse and dense A disagree: %.17g vs %.17g', cls, e_sparse, e_dense); + end + + % --- 16. Reorth-range rejection: 'Reorth' outside {0, 1} is rejected + % by inputParser's validation function before reaching the + % MEX (e.g. Reorth = -3, which would otherwise silently + % behave like Reorth = 1). --- + got_id = ''; + try + randlapack.fun_nystrom_pp(A, k, s, 'Func', 'sqrt', 'Q', 1, ... + 'LFAType', 'scalar', 'Reorth', -3); + catch err + got_id = err.identifier; + end + assert(strcmp(got_id, 'MATLAB:InputParser:ArgumentFailedValidation'), ... + '[%s] Reorth=-3 raised ''%s'', expected ArgumentFailedValidation', cls, got_id); + + % --- 17. Vector-as-diagonal mode: passing the n diagonal entries as a + % vector is the same computation as passing the dense + % diag(lambda), at the same seed, for every LFAType. The point + % is large n, where the dense matrix does not fit; correctness + % is checked here at the fixture size. --- + lam_d = cast(lam, cls); + A_diag = diag(lam_d); + vec_tol = strcmp(cls, 'single') * 1e-4 + strcmp(cls, 'double') * 1e-10; + for it = 1:numel(lfa_types) + lfa = lfa_types{it}; + if strcmp(lfa, 'auto') + extra = {'Budget', budget, 'AutoEps', 1e-2}; + elseif depth_for.(lfa) > 0 + extra = {'Q', 1, 'Depth', depth_for.(lfa)}; + else + extra = {'Q', 1}; + end + est_dense = randlapack.fun_nystrom_pp(A_diag, k, s, 'Func', 'sqrt', ... + 'LFAType', lfa, 'SketchSeed', 11, 'SkipSymCheck', true, extra{:}); + est_vec = randlapack.fun_nystrom_pp(lam_d, k, s, 'Func', 'sqrt', ... + 'LFAType', lfa, 'SketchSeed', 11, extra{:}); + assert(isfinite(est_vec), '[%s/%s] vector-mode estimate is not finite', cls, lfa); + assert(abs(est_vec - est_dense) <= vec_tol * abs(est_dense), ... + '[%s/%s] vector mode %.12e differs from dense diag %.12e', ... + cls, lfa, est_vec, est_dense); + end + + % Row-vector orientation is accepted and gives the same answer. + est_row = randlapack.fun_nystrom_pp(lam_d.', k, s, 'Func', 'sqrt', ... + 'Q', 1, 'LFAType', 'scalar_qfa', 'Depth', 60, 'SketchSeed', 11); + est_col = randlapack.fun_nystrom_pp(lam_d, k, s, 'Func', 'sqrt', ... + 'Q', 1, 'LFAType', 'scalar_qfa', 'Depth', 60, 'SketchSeed', 11); + assert(abs(est_row - est_col) <= vec_tol * abs(est_col), ... + '[%s] row-vector A gave %.12e, column vector gave %.12e', cls, est_row, est_col); + + % The eps-targeted tier runs in vector mode (its Rademacher probes are + % switched to sphere there, since a +-1 quadratic form is exact on a + % diagonal matrix) and produces a genuinely sampled, nonzero error. + est_ad = randlapack.fun_nystrom_pp(lam_d, k, s, 'Func', 'sqrt', ... + 'LFAType', 'adaptive', 'AutoEps', 1e-2, 'SketchSeed', 11); + rel_ad = abs(est_ad - true_sqrt) / true_sqrt; + assert(isfinite(est_ad) && rel_ad < 5e-2, ... + '[%s] adaptive vector-mode rel_err = %.2e exceeds 5e-2', cls, rel_ad); + assert(rel_ad > 0, '[%s] adaptive vector-mode error is exactly zero', cls); + + % Marshal cost is O(n), not O(n^2): a 50,000-entry diagonal marshals in + % milliseconds, where the dense matrix would be 20 GB. + lam_big = cast(logspace(0, -3, 50000)', cls); + [~, ~, ~, times_big] = randlapack.fun_nystrom_pp(lam_big, 8, 4, ... + 'Func', 'sqrt', 'Q', 1, 'LFAType', 'scalar_qfa', 'Depth', 4, 'SketchSeed', 11); + assert(times_big.marshal_in_ms < 100, ... + '[%s] n = 50000 vector marshal took %.1f ms, expected under 100 ms', ... + cls, times_big.marshal_in_ms); + + % Bad shapes are rejected with the documented id. + got_id = ''; + try + randlapack.fun_nystrom_pp(cast(ones(3, 4), cls), 2, 2, 'Func', 'sqrt', ... + 'Q', 1, 'LFAType', 'exact'); + catch err + got_id = err.identifier; + end + assert(strcmp(got_id, 'randlapack:fun_nystrom_pp:A_shape'), ... + '[%s] 3 x 4 A raised ''%s'', expected A_shape', cls, got_id); + + % SpendCap (2026-09-23): the self-tuning tiers spend the cap up to one probe, run every + % probe to the pilot's depth (phase2_checked = 0) and report the pilot's certificate; the + % default call is the shipped tier, with Phase 2 checked. + if strcmp(cls, 'double') + lam_s = logspace(0, -3, 3000)'; + for m = {'scalar_auto', 'block_auto'} + [~, i0] = randlapack.trace_estimate(lam_s, 'Method', m{1}, 'Func', 'log', ... + 'TargetRelError', 0.3, 'MaxMatvecs', 1024, 'Seed', 5); + [~, i1] = randlapack.trace_estimate(lam_s, 'Method', m{1}, 'Func', 'log', ... + 'TargetRelError', 0.3, 'MaxMatvecs', 1024, 'Seed', 5, 'SpendCap', true); + assert(i0.phase2_checked == 1 && i1.phase2_checked == 0, ... + '%s: phase2_checked %g / %g, expected 1 / 0', m{1}, i0.phase2_checked, i1.phase2_checked); + assert(i1.matvecs <= 1024 && 1024 - i1.matvecs < max(i1.d_used, 1), ... + '%s: SpendCap spent %d of 1024 at depth %d', m{1}, i1.matvecs, i1.d_used); + assert(i1.matvecs > i0.matvecs, '%s: SpendCap spent no more than shipped (%d vs %d)', ... + m{1}, i1.matvecs, i0.matvecs); + assert(i1.quadrature_certified == logical(i1.probe_converged), ... + '%s: SpendCap must report the pilot certificate', m{1}); + % the allocation is detail::spend_cap_split's (ported below), not merely some split + % that spends the cap: at n = 3000, cap 4096 and pilot depth 93 it is (1500, 23), where + % the rejected spend-maximizing rule gave (1492, 24) + for nb = 3000 + lam_b = logspace(0, 4, nb)'; + for cap = [1024, 4096] + [~, ib] = randlapack.trace_estimate(lam_b, 'Method', m{1}, 'Func', 'log', ... + 'TargetRelError', 1e-3, 'MaxMatvecs', cap, 'Seed', 5, 'SpendCap', true); + [ke, se] = spend_cap_split(cap - ib.probe_mv, ib.d_used, nb, 1, 4, strcmp(m{1}, 'block_auto')); + assert(ib.auto_k == ke && ib.auto_s == se, ... + '%s n=%d cap=%d t=%d: (k, s) = (%d, %d), spend_cap_split gives (%d, %d)', ... + m{1}, nb, cap, ib.d_used, ib.auto_k, ib.auto_s, ke, se); + end + end + % nothing was checked in Phase 2, so its lower side and certificate are NaN + assert(isnan(i1.tr_L) && isnan(i1.phase2_certified) && isfinite(i1.tr_U) && ... + isfinite(i0.tr_L) && ~isnan(i0.phase2_certified), ... + '%s: SpendCap tr_L %g, phase2_certified %g (expected NaN); shipped tr_L %g', ... + m{1}, i1.tr_L, i1.phase2_certified, i0.tr_L); + end + got_id = ''; + try + randlapack.trace_estimate(lam_s, 'Method', 'block', 'Func', 'log', ... + 'MaxMatvecs', 256, 'Depth', 8, 'SpendCap', true); + catch err + got_id = err.identifier; + end + assert(strcmp(got_id, 'randlapack:trace_estimate:spendCap'), ... + 'SpendCap on a fixed method raised ''%s''', got_id); + end + + fprintf('test_fun_nystrom_pp [%6s]: OK\n', cls); + end +end + +function [k, s] = spend_cap_split(avail, t, n, q, s_min, block_krylov) +% MATLAB copy of RandLAPACK detail::spend_cap_split (rl_fun_nystrom_pp.hh): half of the post-probe +% budget to the rank (at most n - s_min), at least s_min probes, the rank taking what dimension +% limits leave unspent. + t = max(1, t); q = max(1, q); kcap = max(1, fix(n / 2)); + k0 = max(1, min([kcap, n - s_min, fix(avail / (2 * q)), fix((avail - s_min * t) / q)])); + s = min(fix((avail - q * k0) / t), n - k0); + if block_krylov, s = min(s, fix(n / t)); end + s = max(0, s); + k = max(k0, min([kcap, n - s, fix((avail - s * t) / q)])); +end diff --git a/matlab/tests/test_parcopy_threads.m b/matlab/tests/test_parcopy_threads.m new file mode 100644 index 0000000..426ba6a --- /dev/null +++ b/matlab/tests/test_parcopy_threads.m @@ -0,0 +1,12 @@ +function test_parcopy_threads() +% The threaded input copy must never take MATLAB down: with an absurd OMP_NUM_THREADS it has to +% clamp or fall back and still return the same estimate as the plain copy. + rng(3); n = 1024; [Q,~] = qr(randn(n)); A = Q*diag(logspace(0,-4,n))*Q'; A = (A+A')/2; + args = {'Func','log','LFAType','block_qfa','Depth',10,'Q',1,'SketchSeed',5,'SkipSymCheck',true}; + setenv('RANDLAPACK_PERF_PARCOPY','0'); e0 = randlapack.fun_nystrom_pp(A, 64, 16, args{:}); + setenv('RANDLAPACK_PERF_PARCOPY','1'); setenv('OMP_NUM_THREADS','200000'); + e1 = randlapack.fun_nystrom_pp(A, 64, 16, args{:}); + setenv('OMP_NUM_THREADS','4'); setenv('RANDLAPACK_PERF_PARCOPY','0'); + assert(e1 == e0, 'estimate changed: %.17g vs %.17g', e1, e0); + fprintf('TEST_PARCOPY_THREADS PASSED (%.17g)\n', e1); +end