FunNystrom++ V2 - #138
Open
mmelnich wants to merge 31 commits into
Open
Conversation
…entize
Reference-aligned Persson-Kressner port with bit-exact MATLAB↔C++
cross-validation at LAPACK noise floor.
- Phase 5: vendored RandLAPACK::LanczosFA + BlockLanczosFA from
funnystrompp; util::upsize raw-pointer overload, util::symmetrize.
- Phase 6: SASO sketch (RandBLAS::SparseSkOp) + SkOp-aware first matvec
via right_spmm through ExplicitSymLinOp overload.
- Phase 7a: Cholesky-fast / eig(YᵀY) HMT §5.1 path in NystromEVD_v2.
5-rep median speedup vs SVD-pinv fall-back: 2.71×-3.80× across
k∈{500,1000,1500,2000} at n=2000 (3.80× at k=n=2000, exceeds the
~2× claim in funnystrompp PR #132).
- Phase 7a-perf: force_fallback knob + tight t_specrec_ms inner timer.
- Gap port (#1,#2,#3,#5): optional f_zero correction, Phase 2 skip at
k==n, gesdd(G)→syevd(G) in fall-back, SkOp overload on
FunNystromPP_v2::call.
- Componentize: NystromEVD_v2 + workspace moved to drivers/rl_nystrom_evd_v2.hh.
12 cross-validation fixtures pass at tol=1e-11 (10/12 at strict
tol=1e-12; 2 misses are alg_poly.t2 ~1.4e-12 long-standing cancellation
noise). v2 unit tests 4/4. Full RandLAPACK suite 309/309.
Closed
…pace, 11-slot timing Bundle #2 of the v2-polish cycle: * Rename rl_nystrom_evd_v2.hh → rl_nystrom_evd.hh; NystromEVD_v2 free function → NystromEVD; NystromEVD_v2_workspace → NystromEVD_workspace. * Convert workspace from std::vector<T> members to raw T* + int64_t pairs grown via util::upsize and freed in a destructor, matching PR #132 / rl_revd2.hh patterns. Output params now use (T*&, int64_t&) references-to-pointer. * Add 11-slot timing vector (alloc / syrf / matvec / gram / potrf / trsm / svd / post_svd / error_est=0 / rest / total) with per- phase chrono instrumentation; gated on ws.times_enabled. Layout mirrors PR #132's NystromEVD for cross-driver tooling parity. Slot 8 (error_est) is unused in this fixed-k variant; kept at index 8 to preserve the V1 layout. * FunNystromPP_v2: U, lambda, Y_2, fAOmega, Y0 promoted from stack std::vectors to persistent raw-pointer class members; nystrom_ws promoted to a class member so allocations are reused across repeated call() invocations. Destructor frees all five buffers. FunNystromPP_v2's public call signature is unchanged. Verification: * make -j8 + make install -j8 + benchmark rebuild: clean * v2 unit tests: 4/4 PASS (FullRankCapture err_t1 = 2.426e-16) * 12-fixture cross-validation at tol=1e-12: 10/12 pass, same 2 alg_poly.t2 cancellation-regime cells as pre-refactor (rel_Δ = 4.241e-13, 1.399e-12). est at noise floor across all 12 fixtures (rel_Δ ∈ [0.0, 3.835e-15]).
**Silence all RandLAPACK compiler warnings under `-Wall -Wextra`** Resolves all 22 warnings emitted when building the workshop docker image (`ghcr.io/rileyjmurray/randlapack-python-workshop:latest`). **Categories fixed** (-Wunused-but-set-variable, -Wunused-parameter, -Wsign-compare, -Wvolatile, -Wdangling-else): - Removed dead code: `prevnormNR` in `rl_determiter.hh`; unused `[rows_B, cols_B]` structured bindings in `rl_dense_linop.hh` / `rl_composite_linop.hh`. - Removed genuinely vestigial parameters: `lda` from sparse `ABRIK::call` (LDA is meaningless for sparse — `test_ABRIK_general` uses `if constexpr` to dispatch dense vs. sparse), `k` from `test_csr_col_split`, `norm_A` from `test_CQRRPT_orthogonalization`. All callers updated. - For parameters that are unused only in one specialization but required by template-overload symmetry (`density`/`layout` in `make_operator`), used the existing `UNUSED(x)` macro from `RandBLAS/exceptions.hh` with explanatory comments. - Misc: braced `ASSERT_FALSE` in `test_rpchol.cc` for `-Wdangling-else`; rewrote `x += i` on `volatile int` as `x = x + i` to avoid C++20 deprecation in `test_memory_tracker.cc`. Added issue #137 to deal with clang warnings/errors (will involve updating RandBLAS submodule).
… for UW Madison workshop
) Two RandLAPACK sources triggered clang `-Wunused-but-set-variable` and `-Wmisleading-indentation` under `-Wall -Wextra`. GCC 13.3 heuristics do not flag either, which is why they passed through the prior warnings-cleanup PR #136. `RandLAPACK/testing/rl_gen.hh:174`: removed dead `int idx = offset;` and `++idx;` inside `gen_exp_singvals`; the variable was written but never read. `RandLAPACK/drivers/rl_abrik.hh:281, 286, 415, 419`: de-indented the bodies of four `#ifdef RandBLAS_HAS_OpenMP` blocks so they align with the surrounding function-body indent. The previous extra indent visually chained off the preceding single-statement `if(this->timing)`, which is what triggers `-Wmisleading-indentation`. No control-flow change. Verified inside ghcr.io/rileyjmurray/randlapack-python-workshop with clang 20.1.8: 0 warnings under `-Wall -Wextra` (down from 46), 708/708 ctest pass. GCC 13.3 build unchanged.
Two related bugs were present in both the CPU (rl_bqrrp.hh) and GPU (rl_bqrrp_gpu.hh) implementations of BQRRP. They are independent and the fix is in two parts: 1. Termination condition (CPU and GPU) The loop termination used `curr_sz >= n`. For wide A (m < n), curr_sz tops out at min(m, n) = m < n, so this condition never fires on full-rank wide input. The main loop exhausts maxiter iterations and falls through `return 0` without setting `this->rank` or running the termination handler. Result: `rank` is left at its default (0). Fix: use `curr_sz >= std::min(m, n)`. 2. reconcile_and_cleanup wide-tail copy (GPU only) The GPU's pointer-swap dance (Algorithm 7 of the BQRRP paper, used instead of the CPU's in-place column permutation) needs an explicit reconciliation step at termination to bring data home from whichever buffer holds the latest writes. The existing reconciliation: * Walks block-aligned column ranges [0, processed) where processed = (iter + 1) * b_sz_const (or smaller if block_rank != b_sz_const). * Conditionally copies the trailing region past `processed` only when `iter != maxiter - 1` (intended to handle early termination). For wide A at full completion: processed = m < n, so a valid R12 tail exists in [m, n) but the existing condition `iter != maxiter - 1` skips the copy (full completion has iter == maxiter - 1). Additionally, the J pivot vector has no trailing copy at all on the odd-iter path, leaving stale iter-(k-1) values in J[m:n] for wide A. Since iter k's col_swap pulls values from those positions into earlier slots, the final J ends up with duplicates. Fix: replace the `iter != maxiter - 1` gate with `processed < n`, which handles both early termination and wide A at full completion. Add a symmetric J wide-tail copy on the odd-iter path. The same fix applies in both the block_zero early-exit and the normal termination paths. The CPU implementation does not need the second fix because it uses in-place column permutation -- the factorization data ends up directly in caller's A across iterations, so loop fallthrough only loses the `rank` field. This explains why the wide CPU performance results in the BQRRP paper (Fig 3.12) were unaffected by the bug -- timing didn't depend on the rank output. Tests added: * test_bqrrp.cc::BQRRP_wide_aspect (CPU, m=2000, n=5000) * test_bqrrp_gpu.cu::BQRRP_GPU_wide_aspect (GPU, m=1000, n=2000) The upstream test suite previously had no wide-A coverage (all tests were tall m=5000, n=2800 or square 1000 x 1000), which is why this bug went undetected. Verified via an equivalent fix in a standalone CUDA port of BQRRP: residual went from ~5e-1 (with both bugs) to ~1.5e-15 (machine epsilon) on a m=256, n=512 Gaussian input. The reverse-engineered fix in that port matches the patches here. Co-authored-by: Maksim Melnichenko <mmelnichenko@icsi.berkeley.edu>
Route every PR to @mmelnich as code owner. Pairs with forthcoming branch protection on main that requires Code Owner review.
Converted 150+ printf statements to std::cout across library, tests, and benchmarks (not using streams where applicable yet). Issue #75.
…ion (#140) Four coordinated changes preparing RandLAPACK for the bindings repo: install.sh dependency discovery via BLASPP_INSTALL_DIR / LAPACKPP_INSTALL_DIR / RANDOM123_INSTALL_DIR env vars; skip clones and reuse existing installs when set. Also fixes a CUDA flag conflict when external blaspp depends on CUDAToolkit. New RandLAPACK/rl_exceptions.hh defines RandLAPACK::Error plus a single `randlapack_require(cond) << "msg " << x;` macro backed by a RAII StreamThrower whose destructor throws at end of full expression. Registered in the umbrella header and CMakeLists.txt. The macro uses `for(; !(cond); )` rather than `if(cond);else` to short circuit without triggering -Wdangling-else inside unbraced outer ifs. 211 call sites across 19 RandLAPACK source files migrate from randblas_require to randlapack_require with hand readable messages that include the actual parameter values, e.g. `randlapack_require(b_sz >= 1) << "BQRRP::call: b_sz=" << b_sz << " must be positive"`. Polarity flips from "error if BAD" to "require GOOD"; (long long) casts drop since operator<< handles int64_t directly. Driver entry point validation in every driver's call() and constructor, closing the segfault paths that bindings would expose: BQRRP (8 checks), RSVD (5), CQRRT (7), CQRRPT (8), REVD2 (6 across 2 overloads), ABRIK (9 across 3 overloads), krill_full_rpchol (6 new + 1 retained mu_size).
Coordinates RandLAPACK with RandBLAS PRs #158, #169, and #172. Each change stands alone; the commit history reflects that grouping. ### Submodule bump `6abdfcf` → `8417f4b` (current `origin/main`), covering #158 (test infrastructure refactor), #160, #164, #169 (CMake improvements), #170, #171, and #172. ### Required follow-on changes from #169 **`Random123::Random123` rename.** RandBLAS #169 retired the bare `Random123` INTERFACE target in favor of the namespaced `Random123::Random123 IMPORTED GLOBAL`. The `target_link_libraries` entry in `RandLAPACK/CMakeLists.txt` is updated to match. **`RandBLAS_DIR` layout.** RandBLAS #169 moved the installed config from `lib/cmake/RandBLASConfig.cmake` to `lib/cmake/RandBLAS/RandBLASConfig.cmake`. The fallback in `RandLAPACKConfig.cmake.in` now points at `${CMAKE_CURRENT_LIST_DIR}/../RandBLAS` instead of `${CMAKE_CURRENT_LIST_DIR}/..`. ### Test migration (forced by #158) RandBLAS #158 promoted the test-only comparison and sparse-data utilities from `RandBLAS/test/` headers into the library proper under `RandBLAS/RandBLAS/testing/`. They install with the library; consumers reach them via the `RandLAPACK::RandLAPACK` transitive include path. Eleven files updated (10 RandLAPACK tests + 1 extras test): | Before | After | |---|---| | `"../../RandBLAS/test/comparison.hh"` | `<RandBLAS/testing/comparison.hh>` | | `"../../RandBLAS/test/test_datastructures/test_spmats/common.hh"` | `<RandBLAS/testing/sparse_data.hh>` | | `test::comparison::` | `RandBLAS::testing::` | | `test::test_datastructures::test_spmats::` | `RandBLAS::testing::` | `extras/test/CMakeLists.txt` drops the now-obsolete probe that added the RandBLAS source tree to the include path to resolve a test-only header that wasn't installed. Post-#158 those headers are part of the library install and reach consumers through the package's normal include path; re-pointing the probe at the new source location would re-introduce the same anti-pattern, so it's deleted. ### Consumer-facing CMake target improvements (ported from #169) - Build-time-only flags wrapped in `$<BUILD_INTERFACE:...>` so `-Wall -Wextra` (and `-fsanitize=address` when `SANITIZE_ADDRESS=ON`) apply in-tree but disappear from the installed export. - `target_compile_features(RandLAPACK INTERFACE cxx_std_20)` advertises the language requirement on the exported target. - `RandLAPACKConfig.cmake` calls `find_dependency(OpenMP COMPONENTS CXX)` directly when RandLAPACK was built with OpenMP, rather than relying on RandBLAS' Config to incidentally pull it into scope. - `EXPORT_LINK_INTERFACE_LIBRARIES` (deprecated since CMake 2.8.12) removed from `install(EXPORT)`. ### `RandLAPACK_BUILD_TESTS` option New `option(RandLAPACK_BUILD_TESTS "Build RandLAPACK's test suite" ON)` gates `add_subdirectory(test)` at the top level. The name is prefixed so it can coexist with RandBLAS's own `BUILD_TESTS` when RandBLAS is embedded as a subdirectory; the two options control their respective test subdirs independently. `install.sh` is updated to pass `-DBUILD_TESTS=OFF -DRandLAPACK_BUILD_TESTS=ON` to the RandLAPACK cmake invocation: the end-user installation path doesn't run RandBLAS's tests anyway (they're covered by RandBLAS's own CI) and skipping them shaves measurable time off the install.
rileyjmurray
requested changes
May 30, 2026
rileyjmurray
left a comment
Contributor
There was a problem hiding this comment.
I recommend you have Claude pull all changes from the benchmark/ files on this branch and move them to a new PR. Then you can reduce the scope of this PR to actually be about FunNystrom++.
…e_bin helpers to testing utils, add LanczosFA/NystromEVD timing + reorth flag
…Alg 1 / Alg 2 / Lanczos-FA listings); code simplifications
…o-copy q=1 sketch, ping-pong subspace iters, lower-only block tridiag), pseudocode-tagged comments
…-sketch mode and the FunNystromPP SkOp overload; document fscalar/fAfun contract and symmetrize rationale
…-sketch mode and the FunNystromPP SkOp overload; RandBLAS-only RNG + raw buffers in test/benchmark; document fscalar/fAfun contract and symmetrize rationale
… at q=1, no densification), ws.Q rename kills the sketch alias, SkOp matvec required at compile time; drop wasted zero-inits
… ws.Q rename, SkOp matvec required; BlockLanczosFA: assemble T_k in place during the recurrence (drop A_blk/B_blk staging); drop wasted zero-inits
…ernalize normalized-Gaussian Omega2; FunNystromPP use_qfa switch + QFA/FA consistency test
…&C gesdd bug per BALLISTIC thread), loosen CQRRT orthogonality atol eps^0.7->eps^0.65
…owed relative-change stop, d_used); run_lanczos gains an optional per-step stop callback + steps_run; adaptive test
…inal syevd on adaptive stop and the (d*s)^2 copy in the fixed path (preserve_T flag); drop prose em-dashes in QFA/driver comments
…ccuracy, ~3 fewer steps; floor adaptive_min+delay = 4)
…tive-QFA depth probe + budget-closing allocation, rank cap n/2, s <= n-k, batched Phase-2 QFA); BlockLanczosFA: throw on block wider than n; BlockLanczosQFA: adaptive_delay default 5 -> 2
…stop plus fixed-depth mode; audit fixes (scalar-QFA include, non-adaptive pivot guard, auto-tier doc, k->n potrf message)
…e-2 skip; name the block-QFA adaptive defaults; time reorthogonalization explicitly
Four fixes, three of them found while auditing for solver reuse.
1. NystromEVD was always called with vec_nnz (default 8) regardless of k. A
SASO with more nonzeros per column than the column has entries is rejected
by RandBLAS ("(vec_nnz <= dim_major) was required"), so ANY call with k < 8
threw instead of degrading. This is not an exotic corner: the benchmark's
smallest budgets give k = B/2 = 5 at B = 10, and the knob-free auto tier
picks its own k and lands below 8 at small budgets (47 of 221 rungs in a
full-campaign rehearsal died this way). Clamp to min(vec_nnz, k); vec_nnz = k
is a dense column, the correct degenerate limit of a SASO.
2. t_fafun_ms was assigned only inside the k < n branch, so on a k == n call
(Phase 2 skipped) it kept the previous call's value and consumers computing
assembly = t_phase2_ms - t_fafun_ms got a negative number. Latent while every
call constructed a fresh driver; a wrong-answer bug under reuse.
3. BlockLanczosQFA::adaptive_delay / adaptive_min are now exposed as named
default_adaptive_* constants so a caller reusing an oracle can restore them
explicitly. Relying on construction to supply the defaults leaks the previous
call's convergence window forward, changing d_used, the estimate and the
matvec count.
4. LanczosFA / BlockLanczosFA record _t_reorth_us and widen `times` to six slots
(matvec, run_lanczos, apply_f, rest, total, reorth). Full reorthogonalization
is O(d^2 n s) against the recurrence's O(d n s) and is the entire difference
between this oracle and the basis-free LanczosQFA, so it is now a reported
number rather than something inferred by differencing two tiers. The QFA
classes emit 0 in that slot by design (no reorth).
Tests: SmallRankBelowVecNnzDoesNotThrow, FafunTimerResetAtKEqualsN and
ReuseAcrossCallsIsBitIdentical (a grow/shrink/grow sequence of (k, s) against
one reused driver, compared bit-for-bit against fresh instances) are new; the
second fails against the pre-fix driver. 18/18 TestFunNystromPPv2 pass.
# Conflicts: # RandLAPACK/CMakeLists.txt # RandLAPACK/drivers/rl_abrik.hh # install.sh # test/comps/test_determiter.cc # test/comps/test_qb.cc # test/comps/test_rpchol.cc # test/drivers/test_cqrrt.cc # test/drivers/test_krill.cc
…oes not exist on MSVC) BinaryIoRoundTrip hardcoded /tmp and POSIX mkstemps/close, the only compile failure in the whole branch under MSVC. save_dense_bin creates the file itself, so no create-then-close dance is needed. Verified: Linux suite 331/331; Windows (MSVC 19.44 + oneMKL ILP64) 767/767 via install.ps1.
…eans budget/n-bounded) The fixed default silently floored the oracle bias whenever the certified depth exceeded 200: in the 2026-07-29 campaign (kappa >= 1e6, certified depth 700-2200) the auto tier stagnated 10-38x above the 0.1 percent target across the entire budget ladder. The probe is now bounded only by n and half the budget; a positive auto_depth_cap remains an opt-in fixed cap. Regression test pins the old failure (probe must exceed 200 on a geometric kappa=1e6 spectrum at eps=1e-6; budget closure unchanged). Suite 332/332.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
WIP