Conversation
Contributor
Author
|
Closing this & adding a reference to the paper. |
Contributor
Author
|
Temporarily re-opened to add functionality. |
When solving for R_sk^{-1} explicitly to precondition a linop, solve
R_sk * X = I (Side::Left) rather than X * R_sk = I (Side::Right).
The two are mathematically equivalent but the latter exhibits poor
backward error in the subsequent product A * R_sk^{-1} when R_sk is
ill-conditioned. On photogrammetry2 (κ(R_sk) ≈ 2.2e8) this changes
orth_err in CQRRT_linop from 1.6e-4 to 1.5e-9 -- matching the
GEQP3/BQRRP stabilized variants and obviating their need as the
default path.
Same fix applied to the initial M = R_1^{-1} formation in
sCholQR3_linops.
Slim CQRRT_linop_applications to run only the patched CQRRT_linop;
swap GETRI for BQRRP in CQRRT_diagnostic to show the stabilized
counterpart alongside trsm/trtri.
…b bugs
- Replace NMR/Kronecker benchmark with CQRRT_linop_irlsq.cc.
Loads a tall sparse .mtx, wraps as SparseLinOp, generates synthetic
x_true + b = J*x_true + noise. For each Q-less QR variant: draws a
fresh sparse sketch S2 (independent of CQRRT's S1), forms x_0 =
R^{-1} R^{-T} (S2 A)^T (S2 b) (paper Algorithm 1, line 3, Q-less
form), then runs 2-step IR with inner CG.
- Strip Tikhonov from IterRefineLSQ. Drop KroneckerOperator,
RegularizedLinOp, and their tests/includes/CMake entries. Drop the
GEQP3-stabilized variant from the benchmark and rename
CQRRT_linop_stb_bqrrp -> CQRRT_linop_bqrrp.
- Fix IR-LSQ populate_times double-counting bug. Inner CG's
TRSM/fwd/adj contributions were tracked into both t_inner_total
(wallclock of the inner_cg call) and the shared
t_{trsm,fwd,adj}_total counters, so 'other = outer_total - inner -
trsm - fwd - adj' went negative. Split into outer-only and
inner-only counters and report inner_ex = t_inner_total -
t_inner_{trsm,fwd,adj}, total_trsm = t_outer_trsm + t_inner_trsm
(similarly fwd, adj). 'other' is now a non-negative residue
(axpy/copy/nrm2 bookkeeping).
- Fix three analytical_kb formulas in rl_memory_tracker.hh.
cqrrt_linops_bqrrp_analytical_kb captured only the BQRRP-precond
moment (d*n + 4n^2) and missed the later Gram-loop moment (d*n + n
+ n^2 + m*b_eff). For tall inputs the latter dominates. Now
returns the max of the two; signature extended to (m, n, d_factor,
block_size). scholqr3_linops_analytical_kb forgot G3_factor at the
iter-3 peak (5n^2 -> 6n^2). scholqr3_linops_basic_analytical_kb
forgot all three G_i_factor members (3n^2 -> 6n^2).
Per collaborator's correction (Oleg, 2026-05-25), the FEM2 operator is the
Petrov-Galerkin form J = B^{-1} * A * U_h with B = chol(M), A = K, U_h = V
(mass-matrix Cholesky factor, stiffness, prolongation respectively).
CQRRT_linop_applications now takes three .mtx files in FEM mode and builds
a doubly-nested CompositeOperator:
J = CompositeOperator(L_inv_op,
CompositeOperator(K_op, V_op))
with L_inv_op = CholSolverLinOp<T>(M_file, half_solve=true). Because L
factors M (not K), the composite does NOT algebraically collapse; the LS
normal equations land on J' J = V' K M^{-1} K V (Petrov-Galerkin coarse-
grid mass-weighted stiffness-squared).
CLI changes:
FEM mode (NEW): <K_file> <M_file> <V_file> <d_factor> [...] (3 files)
Sparse mode: sparse <A_file> <d_factor> [...] (unchanged)
The previous 2-stage composite (L^{-1} * V with L = chol(K)) is gone --
that formulation was based on a misread of the generator's outputs. The
sparse mode is unchanged.
No driver, linop, test, or CMake changes.
- Single binary now selects between SVD post-processing, IR-LSQ refinement,
or both via a positional <mode> arg.
- 5-method dispatch (CQRRT_linop, CholQR, sCholQR3, sCholQR3_basic,
CQRRT_linop_bqrrp) restored to the FEM/sparse composite paths.
- Add power-iteration estimate of ||A||_2 and replace ls_residual_norm
with the Higham normwise backward error ||Ax-b||/(||A||*||x||+||b||),
drivable to machine epsilon for a backward-stable LS solver.
- Add memlite blocked-compute orth_err (O(n^2 + m*b)), runs for every
selected method in every mode; new orth_error column in irlsq CSV.
- FEM + irlsq: b = L^{-1} * Gaussian random vector (no x_true).
- CQRRT_linop_irlsq.cc removed (CMake target dropped).
…s spec)
Introduces comps/rl_cholqr.hh with three free-function templates:
- blocked_preconditioned_gram(A, R_pre, G, ...) : Layer 0; computes
G = R_pre^T A^T A R_pre via blocked linop calls. nullptr R_pre handles
the M=I case via a small per-block identity scratch.
- cholqr_primitive(A, R, shift_factor, ...) : Algorithm 1 (with optional
shift for sCholQR3 iter 1).
- pcholqr_primitive(A, P, R, method, ...) : Algorithm 2; the unifying
building block. Dispatches the P^{-1} step on PCholQRPrecondMethod
(TRSM_IDENTITY / TRTRI / GEQP3 / BQRRP). Adds the new TRTRI method.
Gram step in pcholqr_primitive dispatches on method:
- TRSM_IDENTITY/TRTRI: per-block writes A^T A R_pre into G, then a
single TRSM(P^T, G) applies the left factor. O(n^3/2) vs O(n^3),
stable since P is preserved (the optimization CQRRT_linops used to
have inline; now shared with sCholQR3 iters 2-3).
- GEQP3/BQRRP: per-block GEMM with explicit R_pre^T, preserving the
QRCP stability advantage.
Driver refactor:
- CholQR_linops, sCholQR3_linops (both variants), CQRRT_linops are now
thin wrappers over the primitives. ~33% LOC reduction across the
family (1946 -> 1310 lines).
- rl_cqrrt.hh now holds both dense CQRRT and CQRRT_linops in one file;
rl_cqrrt_linops.hh deleted (CMake had no separate target; 3 benchmark
.cc files updated to drop the include).
- Backwards-compat alias `using CQRRTLinopPrecond = PCholQRPrecondMethod;`
keeps existing benchmark code unchanged.
Memory tracker formulas updated to reflect the now-freed buffers:
scholqr3_linops_analytical_kb: 6n^2 + (m+n)*b -> 2n^2 + (m+n)*b
scholqr3_linops_basic_analytical_kb: m*n + 6n^2 -> m*n + 2n^2
Tests:
- 3 new tests in test_orth_linop.cc cover TRTRI / GEQP3 / BQRRP precond
paths (existing tests only exercised TRSM_IDENTITY).
- All 26 CholQR/sCholQR3/CQRRT tests pass.
PowerOp<InnerOp> (RandLAPACK/linops/rl_power_linop.hh) — generic wrapper representing A^j for a square base linear operator A. Chains j calls to the base op with two ping-pong scratch buffers; A^j is never materialized. j == 1 takes a no-scratch fast path; j == 2 allocates one scratch; j >= 3 allocates two. Op::Trans dispatches base(Op::Trans, ...) j times. Side::Left only (the only consumer pattern we have so far). sparse_axpby_shared_pattern (extras/misc/ext_sparse_axpy.hh) — computes C := alpha*A + beta*B for two CSRMatrix inputs whose sparsity patterns are bit-identical (rowptr + colidxs equal). O(nnz) value-only path. The target consumer is X = K - omega*M in the reduced-spectral application, where K and M from a single FEM mesh always share sparsity exactly. A general-purpose sparse_axpby for different patterns belongs upstream (RandBLAS issue; MKL has mkl_sparse_d_add, cuSPARSE has cusparseDcsrgeam2). Tests: 6 new PowerOp tests (j=1, j=3, multi-RHS, Op::Trans, alpha/beta, PowerOp wrapping CompositeOperator — the rspec usage pattern). 3 new sparse_axpby tests (tridiagonal, shifted-inverse pattern, float type). All 9 pass; full test suite still passes.
Sibling to CholSolverLinOp. Wraps Eigen::SparseLU with COLAMD ordering,
factor-once / solve-many pattern. Handles both SPD and indefinite sparse
matrices — used by the reduced-spectral application when omega is an
interior shift and X = K - omega*M becomes indefinite (Cholesky fails).
Scope:
- In-memory Eigen::SparseMatrix constructor (rspec mode computes X at
runtime via sparse_axpby; the file-based constructor mirroring
CholSolverLinOp's pattern can be added when a consumer needs it).
- operator(): Side::Left, ColMajor, Op::NoTrans on B, Op::NoTrans /
Op::Trans on A. RowMajor and Side::Right deferred.
Tests cover SPD tridiagonal (1D Laplacian), indefinite tridiagonal,
non-symmetric matrix (exercises trans dispatch, A^{-T} != A^{-1}), and
multi-RHS with alpha/beta accumulation. All 4 pass.
TransposedOp<InnerOp> (RandLAPACK/linops/rl_transposed_linop.hh) — implicit
transpose view of any LinearOperator. One-line dispatch wrapper: forwards
to base() with the trans flag flipped. Generic over the LinearOperator
concept so it composes with DenseLinOp, SparseLinOp, CompositeOperator,
PowerOp, and even itself (double-transpose tests pass).
Use case from the rspec application: build C = L^T * X^{-1} * L as
CompositeOperator(TransposedOp(L_op), CompositeOperator(X_inv_op, L_op))
without materializing a separate L^T sparse matrix. But TransposedOp is
intentionally generic — any future 3+ operand chain in our codebase that
needs a transpose-on-one-operand will reuse it.
Also: added the concept-required 12-arg operator() overload (no Side,
delegates to Side::Left) to both PowerOp and TransposedOp. Without this
overload, nesting these wrappers (e.g., PowerOp around PowerOp, or
TransposedOp around PowerOp) fails the LinearOperator concept check.
Other linops (DenseLinOp, SparseLinOp, CompositeOperator, CholSolverLinOp)
already had this overload; the new wrappers now match the convention.
Tests: 5 new TransposedOp tests (dense, double-transpose==identity,
CompositeOperator inner, transposed-of-transposed, around-PowerOp).
All 11 PowerOp + TransposedOp tests pass.
…mark
Implements Algorithm 4 from the collaborator's pseudocode: reduced-basis
Rayleigh-Ritz approximation of eigenvalues of the symmetric operator
C = L^T * (K - omega*M)^{-1} * L, where L L^T = M
on the subspace range(V_app), with V_app = C^j * V_FEM.
New mode "rspec" alongside the existing svd/irlsq/both. Two new positional
CLI args: <omega> (double, default 0.0) and <power_j> (int, default 1,
constrained to {1, 2, 3} per collaborator spec). FEM input only — sparse
mode rejected for rspec.
Operator chain assembly (rl_cqrrt_applications.cc, run_rspec_benchmark):
X = sparse_axpby_shared_pattern(K, -omega*M) // O(nnz)
X_eigen = convert(X) // triplets
X_inv_op = SparseLUSolverLinOp(X_eigen).factorize() // try/catch
L_op = SparseLinOp(L_inv_op.make_L_csc()) // L from chol(M)
C_op = CompositeOperator(TransposedOp(L_op),
CompositeOperator(X_inv_op, L_op)) // L^T X^{-1} L
Cj_op = PowerOp(C_op, power_j)
V_app_op = CompositeOperator(Cj_op, V_op) // implicit, no mat
Per-(algorithm, run): 5-method dispatch (CQRRT_linop, CholQR, sCholQR3,
sCholQR3_basic, CQRRT_linop_bqrrp) reuses the same QR drivers as the
svd/irlsq paths. After QR returns R, Rayleigh-Ritz forms the small n*n
matrix T = R^{-T} * V_app^T * C * V_app * R^{-1} block-by-block (no mxn
intermediate materialization), symmetrizes against rounding drift, then
syevd gives eigenvalues and eigenvectors. Ritz residuals computed for
top-k pairs as ||K v - lambda M v|| / (||K v|| + |lambda| ||M v||) where
v = V_FEM * R^{-1} * u.
Output: <ts>_rspec_results.csv with columns
algorithm, run, m, n, omega, power_j, qr_status, qr_time_us, peak_rss_kb,
analytical_kb, factor_time_us, rspec_total_us, eig_0..eig_{k-1},
resid_0..resid_{k-1}.
Singular-X handling: when omega is too close to an eigenvalue of (K, M),
SparseLU.factorize raises RandLAPACK::Error. Caught at the top of
run_rspec_benchmark; a single failure row with qr_status=-99 is written
and the run returns 0 (not an error — collaborator flagged this as an
expected case).
Fix to PowerOp + TransposedOp: const-correctness on the dense B input
(T* const B -> const T* B) so they compose under CompositeOperator,
which always passes B as const T* through its body.
ext_cholsolver_linop and ext_sparselu_linop minor cleanup.
Build: clean. 36 existing CholQR/sCholQR3/CQRRT/PowerOp/TransposedOp tests
pass. End-to-end smoke test on the small FEM2 problem (75824 x 8304,
omega=0, j=1, method_mask=1) progresses through matrix load + L factor
(310 ms) + X=K-omega*M assembly + SparseLU factor (660 ms) + composite
operator construction. PCholQR warmup is in progress — execution-time
proof that the data flow is correct; the actual benchmark needs the
ISAAC walltime budget (rspec at this size is dominated by the
SparseLU(75824).solve(8304 RHS) at each C-application).
A unit test that compares Ritz eigenvalues against lapack::sygv on a
small synthetic problem is left as a TODO; the FEM smoke test
exercises every new code path.
- CQRRT_linop_applications: remove GSVD post-processing (gesdd on R for
generalized singular values/vectors), upcast-orth diagnostic, and the
"svd"/"both" modes. Available modes are now {irlsq, rspec}. Drops
~410 lines: result-struct SVD fields, A_materialized + AtA_precomputed
setup, do_svd/do_irlsq branching, GSVD CSV writers, skip_svd/upcast_orth
CLI args, and the K_file/V_file params from the inner runner (unused
after the GSVD writer deletion).
- extras/test: delete test_ext_sparse_axpy.cc and drop it from
CMakeLists.txt; the helper it covered is not part of the API surface
we still care about.
6d026d5 added them only to the sparse-mode selector; run_irlsq_reg has its OWN method selector and Blendenpik dispatch, so the FEM2 campaign ran seven rows instead of nine despite method_mask=127 setting bit 64. Same multi-path trap that earlier made an env knob silently miss the code under test. Verified on the real 75466x8256 kc1e10 operator: cold Blendenpik's backward error 2.057e-08 (forward 5.164) becomes 2.255e-16 (7.799e-07) with refinement, matching the warm variant, at the same preconditioner. IterRefineLSQ gains inner_iters_total() so the refined rows report Blendenpik's LSQR iterations plus refinement, as the Toeplitz rows do.
…r-block explicit left factor in preconditioned Gram) and RANDLAPACK_SCHOLQR3_SHIFT=theory (11*eps*n*trace(G) first-pass shift)
…oor semantics Time the per-round NE recomputation and warm-seed ops; add status codes for round-budget and floor exits; advance the outer stagnation reference only on significant improvement and give it its own window knob; best-iterate return on CG breakdown; LSQR stop-test reporting; solve-scoped FFT width matching; batched multi-column Toeplitz FFT applies; corrected analytical storage formulas; input validation and dead-knob removal in the QR drivers.
…unting One dispatch for the refine rows in both LS benchmarks (init_only sketch-and- solve x0, all iterative work in the shared engine): fixes the warm row that ran cold and the phase time missing from every column. New columns: setup/x0 time, lsqr vs CG iteration counts, named stop reasons, inner-kernel vs restart overhead split, whole-row wall clock, x0 quality; per-round sidecar CSVs; full argv/env/commit provenance headers; unified mask decode; extended warmup.
mmelnich
added a commit
that referenced
this pull request
Aug 29, 2026
mmelnich
force-pushed
the
spring-2026-wip
branch
from
August 29, 2026 14:51
87594db to
db21de0
Compare
…ws to the fixed algorithm Cholesky breakdown now reports as a FAIL row (qr_status != 0) instead of silently switching the row to the shift-rescued variant. Applied at all 20 measured driver constructions in the two campaign binaries; warm-ups exempt; value echoed in the CSV env header. Unset keeps the library default (unbounded rescue); library behavior unchanged.
… optional Default unchanged (symmetrize, as the paper prescribes). Setting =0 factorizes the upper triangle as computed (pre-B5 behavior), echoed in the benchmark CSV env line, for A/B campaigns isolating the symmetrization's ULP-level effect on borderline pivots.
…red engine The FEM driver has always handed restarted_pcg_ne an absolute inner target (eps^0.85 times the first cycle's normal-equation right-hand side), so a cycle whose normal-equation residual is already at rounding level returns after one iteration. The Toeplitz driver passed 0.0, so the two stagnation-confirmation cycles that end every row at the data-noise floor each ran a full inner solve on rounding noise: 25 to 45 iterations for Blendenpik and CholQR, against one iteration per cycle on FEM. Reported totals were inflated by up to 5x for the weak preconditioners (CholQR on the middle cell: 92 reported, 1 productive) and were not comparable with the FEM benchmark. All three engine calls (unpreconditioned, Q-less rows, refined Blendenpik) now pass the same guard as FEM, the value is echoed in the CSV header as inner_abs_tol, and the file-header note on knobs differing from FEM is corrected. Convergence points are unchanged; only the cost of confirming stagnation drops.
…N24 eq. 4.2) Both least-squares benchmarks now report the backward error of every computed solution in the sense of Epperly, Meier and Nakatsukasa (2024): the smallest ||[dA, theta*db]||_F making x the exact least-squares solution of the perturbed problem, estimated by the Karlson-Walden formula (their Fact 4.1, within a factor sqrt(2)) with A^T A replaced by (SA)^T(SA) for a sparse sign sketch with d = 2n rows (their eq. 4.2). Two weights are written, both relative to ||A||_F: theta = ||A||_F/||b|| (the value their Algorithm 4 tests at run time) and theta = infinity (perturb A only). The residual-orthogonality ratio ||A^T r|| / (||A||_F ||r||) is recorded next to them. The reference (sketch of the operator plus one d x n SVD) is built once per problem after the last timed row, so no row's timing or peak-RSS window contains it; each row then costs one forward and one adjoint application. ||A||_F is accumulated exactly from the operator's columns while the sketch is formed. Results go to a sidecar CSV (*_backward_error.csv) keyed by (algorithm, run), leaving the streaming results CSV untouched. The FEM driver evaluates the base problem M x ~ b; the Toeplitz driver the augmented problem every row actually solves.
…benchmark Method-mask bit 128 (irlsq_reg only) runs the shared refinement engine on the raw operator with R = nullptr. No factor is built, so the build phase, the Cholesky records, the storage model and the orthogonality entry keep their no-value sentinels and the solve is the whole row. This is the row the Toeplitz benchmark has carried since the original port. The sparse irlsq path rejects the bit and rspec warns on it.
…LI (arg 18, default eps^0.85 as before)
…oor (paced mode); one helper for the three call sites
…tic target None of this had a consumer. Verified against both downstream trees before deleting: no SLURM script invokes CQRRT_diagnostic, no MATLAB plotter reads its output, and the one paper table it fed now lives in the manuscript's archive directory, un-inputted. The dense CholQR drivers (CholQR_dense, CholQR2_dense, sCholQR3_dense) had no caller outside their own test, which said as much. They were an interface veneer over the same Q-less engine, not separate functionality: the non-Q-less implementations (CholQRQ, HQRQ, PLUL in rl_orth.hh, plus CQRRPT/BQRRP/HQRRP) are untouched. Dense CQRRT is a different driver with real callers and stays, so the test file is renamed to match what it actually covers. The audit build directories do not belong in a shared .gitignore; they are local and now live in .git/info/exclude. - delete benchmark/bench_CQRRT_linops/CQRRT_diagnostic.cc and its add_benchmark entry - delete RandLAPACK/drivers/rl_cholqr_dense.hh and its RandLAPACK.hh include - drop the 5 dense-CholQR tests; rename test_cholqr_dense.cc -> test_cqrrt_dense.cc - revert .gitignore to match main
IterRefineLSQ::warm_x0 was never assigned by any caller, so the docstring calling it "the Blendenpik handoff" was false: that handoff goes through restarted_pcg_ne's own x0 parameter. The engine now receives an explicit nullptr. In rl_blendenpik.hh the laset over R's strict lower triangle writes zeros over zeros: R is value-initialized and lacpy(Upper) never touches that triangle.
…egrity guard Both live benchmarks took long positional argument lists (20 slots for the applications driver, 18 for Toeplitz). That interface cost a campaign: a value landed one slot early, the binary parsed and echoed it, and seven cluster jobs silently encoded a different experiment. Named flags make that impossible by construction: a value is bound to a name, an unknown or malformed flag is a hard error, and an omitted one takes a documented default. The positional form still works so existing job scripts keep running, and it now warns. Measurement drove the defaults: across every job script in the repo, only the output directory, m and n ever varied for Toeplitz, so the other 15 arguments became defaults. A full run is now --out=DIR --m=75466 --n=8256. kappa_target is no longer settable. It was 1 in every run that produced data, and the FEM2 generators warn that a larger value re-applies a column scaling already baked into V. The CSV column is retained. Two benchmarks were also dispatching the same five QR methods through copy-pasted blocks that had already diverged: one recorded chol_retries unconditionally, the other only on success. They now share one harvest, which records it unconditionally, since a retry count is meaningful whether or not the factorization succeeded. - add RandLAPACK::bench::BenchArgs, QRRun and run_cholqr_family - add check_csv_arity, run over all 7 CSV outputs: headers are string literals and rows are chained <<, so nothing otherwise ties their field counts together - drop the dead zero-fill from the generic materialize (it is overwritten by a beta=0 apply; the SparseLinOp overload genuinely needs its fill and keeps it)
mkl_set_num_threads_local sets the thread-local count outright, with no minimum against what the caller already has. The guard therefore RAISED the width to 8 (trsv) or 16 (FFT) whenever fewer threads were available, including under OMP_NUM_THREADS=1: it oversubscribed restricted allocations and made single-threaded runs not single-threaded. That was enough to make results irreproducible. Two runs of the same binary with identical arguments disagreed on the long refine rows, and in one comparison on the iteration count (44 vs 43), because the widened regions reduced in a different order. With the minimum in place those comparisons are bit-identical. Campaign allocations request 64, 32 or 16 CPUs against caps of 8 and 16, so the guard was already capping there and no published result changes. Also trims the header commentary to the constraints themselves: the hardware timing tables and the calibration narrative belong in the dev log, not in a library header.
…tput product
Two n-by-n calls were doing three times the necessary work. P^{-1} is formed by solving
P X = I, where X is upper triangular because P is, so column j is zero below row j and a
panel of columns needs only the leading rows of P. The output product R = R_chol * P has
the same structure. Both were running at full width: n^3 flops where n^3/3 suffices.
Panel-blocking them measures 1.653 s -> 0.840 s at n=8304, b=256, 8 threads.
NOT bit-identical, despite the arithmetic on the nonzero part being the same. BLAS selects
different internal blocking for a different operand shape, which reorders accumulation:
measured 1.4e-17 at n=5 rising to 1.4e-11 at n=512 against the full-width result. A panel
width at or above n degenerates to the original single call and IS exactly identical, which
the new tests pin.
Results from the affected benchmarks move at ULP level and need a rerun before publication.
- add invert_upper_into and trmm_upper_upper_left in comps/rl_cholqr.hh
- add test/comps/test_triangular_panels.cc: agreement to 1e-12 relative across 6 sizes and
4 panel widths, exact structural zeros, a P*P^-1 = I check, and the bit-identity boundary
- collapse six byte-identical total_us docstrings to one statement of the constraint
…e path Every right-hand side is real, the circulant embedding is real, and the result was taken as the real part, yet the transform was complex-to-complex: it computed a conjugate-symmetric spectrum and discarded half of it. A real-to-complex pair computes L/2+1 bins instead of L. Measured 2.54x on one forward/multiply/backward pair at L = 2^19 on 8 threads. Applied to the single-RHS path only, which is the CG and LSQR inner loop and issues by far the most applies. The batched multi-column path is untouched. The spectra stay full length and are still built by the complex descriptor: for a real embedding F[k] = conj(F[L-k]), so their leading L/2+1 entries ARE the half spectrum, which keeps one source of truth for the embedding. Correctness checked three ways: the existing dense-reference tests; a new test comparing the real-to-complex single-column path against the untouched complex-to-complex batched path at m=2000,n=400 and m=8000,n=1000 in both transpose directions (agreement better than 1e-14); and a standalone harness agreeing to 1.8e-12 on an output of magnitude ~700. Results are unchanged at the default inner floor: the unpreconditioned row, the most transform-sensitive one, gives the same 71 iterations and a last-digit residual difference. At an inner floor of 1e-16 the two transforms diverge (145 iterations to 2.6e-11 versus 111 to 8.9e-10), which is a property of that setting rather than of either transform: with the floor effectively disabled, how far the solve descends depends on rounding, and both residuals sit near the 1e-11 data noise level.
Adds BackwardErrorOracle<T> and be_tol. After each round, where the loop already holds b - A x and A^T r, the oracle is evaluated; a value <= be_tol ends the run with status 5. The LS tolerance keeps precedence on every path. Per-round oracle values, the x0 value and the oracle wall time land in PCGRoundHistory; the wall time is excluded from times[3]. IterRefineLSQ forwards the oracle and republishes the new fields; its return contract is unchanged. Inactive by default.
…gine's stop test Both least-squares drivers gain --be-tol-mult (positional slot 16 for FEM, arg 19 for Toeplitz; 0 = off, today's behaviour). When on, the Karlson-Walden reference is built once before the rows and wrapped as the oracle with tolerance mult*sqrt(n)*u; the post-pass sidecar reuses it. New CSV columns: be_kw (rounds), t_be_us and be_x0 (results); stop reason "be"; header echo of the resolved tolerance.
… exit, IterRefineLSQ forwarding
…ks the oracle exit Evaluating the oracle at the ambient MKL width inside the width-pinned solve forced two OpenMP team re-formations per round, whose cost landed in the next capped region's time. A CG breakdown in the triggering round now reports status 2. Tests: x identity and count consistency at a status-5 exit with the oracle's own time kept out of times[3]; an active but unmet oracle changes nothing; the LS tolerance outranks the oracle in the same round.
…nit roundoff; be_final column The reference sketch never applied the sparse-sign isometry scale, so every estimate was understated by up to sqrt(zeta) in the converged regime and the stop test was looser than stated. The oracle sketch is pinned at 8 nonzeros per column. be-tol-mult now scales sqrt(n)*u with u = eps/2. Both results CSVs gain be_final, the estimate of the returned iterate. Degenerate estimator values map to +inf; a float solve precision with the oracle on is refused.
…imator, float refusal at parse time be_final is +inf when the oracle was on but no round ran and no warm start was evaluated, so it can never read as below the target; the estimator maps a NaN to +inf in the sidecar as well. --be-tol-mult > 0 with a float solve precision is refused before the matrices are loaded. Both headers echo the oracle's sketch density. Test: a CG breakdown in the triggering round outranks the oracle.
mmelnich
changed the title
Conflicts resolved: - install.sh: take main's root wrapper; the installer now lives in installers/install.sh (the branch's edits were comment punctuation only). - RandLAPACK/drivers/rl_cqrrt_linops.hh: deleted; the operator driver was folded into rl_cqrrt.hh as CQRRTO_linops. Main's (T) scalar casts are applied there instead. - RandLAPACK/drivers/rl_cqrrt.hh: keep the branch layout, add main's (T) casts on the trsm/syrk/trmm scalars. - test/CMakeLists.txt: union of both sides (branch test sources and the DFTI/ILP64 block; main's test_config.cc, /bigobj and runtime DLL staging). Also retire the branch's last two __APPLE__ guards (cholqr_primitive's GEQP3/BQRRP refusal and the test that skipped precond_method_BQRRP), per main's "macOS runs the full library" policy; that refusal is what failed core-macos on every push since 09-01. RandLAPACK/linops/rl_materialize.hh: take main's generic fallback as is (zero-fill kept, identity allocation capped) and its guard test unchanged. The branch had dropped the fill for benchmark-scale callers, but the only such caller is CQRRTO_linop_basic's orthogonality check and the cap now bounds that path anyway; the era benchmark never routes through it. Suite: 425/425 locally.
…e augmentation
Engine: pcg_inner gains an optional poll hook (PCGInnerControls::poll_every,
::poll) and the exit status OracleMet; restarted_pcg_ne forwards a new trailing
be_poll_every and, when the backward-error oracle is active, evaluates it every
be_poll_every inner iterations on the trial iterate R^{-1}(z + dz), formed with
the round-end fold's own arithmetic in dedicated buffers, and ends the run with
status 5 mid-round. The poll's triangular solve and two operator applies count
as solve time inside the kernel slice; the oracle's own time stays excluded.
History records polls per round and the poll wall time. IterRefineLSQ and
run_refined_blendenpik forward the knob. Default 0 keeps every existing path
bit-identical.
Toeplitz driver: --be-poll-every (positional slot 20), rejected without a
positive --be-tol-mult; header echo; rounds CSV gains be_polls and the
OracleMet legend; results CSV gains t_poll_us.
FEM driver: --mu-factor=0 runs the Q-less factorizations on J itself through
one generic lambda over both operand types; the mu column stays (0) and the
header says so. Rounds CSV plumbing for the per-round poll counts.
Tests: two kernel, three engine, one adapter and one refined-Blendenpik test,
each written before its implementation; full suite 425 -> 432.
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
1 participant
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.
Summary
This PR turns the Q-less randomized QR work into a coherent line in RandLAPACK: the operator-based Q-less driver is now called
CQRRTO_linops(the paper's CQRRTO), the CholQR family is rebuilt on one shared primitive, a restarted preconditioned normal-equations solver with named exit conditions is added, and the two least-squares benchmarks that produce the paper's figures (a FEM regularized problem and a matrix-free Toeplitz autoregression problem) live inbenchmark/. Every experiment knob is on the command line or in an environment variable and is echoed into the CSV headers, so a result file identifies its own configuration.Behaviour on the existing
mainAPI: the denseCQRRTdriver and ABRIK'sQR_explicit::cqrrtkeep their names;CQRRT_linopsremains as a deprecated alias ofCQRRTO_linops. Relative tomain,rl_cqrrt_linops.hhis folded intorl_cqrrt.hhand the generalized-LS composite benchmark is removed.Changes
Library
rl_cqrrt.hh: denseCQRRTand operator-basedCQRRTO_linops(formerlyCQRRT_linops), both with the adaptive Cholesky shift retry and its shift record;dclamped to at leastn.comps/rl_cholqr.hh: onecholqr_primitivewith optional preconditioner and adaptive shift;CholQR,CholQR2,sCholQR3andsCholQR3_basicare thin wrappers over a sharedcholqr_iterateengine. Knobs:RANDLAPACK_CHOL_MAX_RETRIES(retry cap;0pins a row to the fixed algorithm),RANDLAPACK_CHOL_SYMMETRIZE(optional Gram symmetrization),RANDLAPACK_SCHOLQR3_SHIFT.drivers/rl_restarted_pcg_ne.hhandrl_pcg_inner.hh: restarted PCG on the right-preconditioned normal equations with a stable round residual (R^{-T} A^T (b - A x)), per-round restart pacing, inner stagnation exit returning the best iterate, an optional absolute inner floor, an outer LS-floor exit, and an optional backward-error oracle exit; every exit has its own status. Per-round history and a non-overlapping timing breakdown are returned.drivers/rl_iter_refine_lsq.hh:IterRefineLSQis now a thin adapter over that engine.drivers/rl_lsqr.hhandrl_blendenpik.hh: matrix-free LSQR and Blendenpik (LSQR with a sketch-QR preconditioner), warm or zero start, Cholesky retry-count reporting;benchmark/refined_blendenpik.hhruns Blendenpik's preconditioner and start through the shared engine.VStackOp([A; mu I], so the augmented operator can be handed to a sketch-based driver directly),PowerOp,TransposedOp,ScaledIdentityOp,materialize, plusextrasoperators for a Cholesky-solver composite (L^{-1} K V) and an FFT Toeplitz operator with a real-to-complex transform.misc/rl_blas2_threads.hh: a thread-width cap for level-2 BLAS and FFT regions, calibrated on the benchmark hardware, applied as a cap rather than a set.testing/rl_memory_tracker.hhandrl_gen.hh: peak-RSS tracker with a per-window baseline, analytical memory models per algorithm, and a fixedgen_bad_cholqr_singvals.Benchmarks (
benchmark/bench_CQRRTO_linops/,benchmark/bench_toeplitz_ls/)CQRRTO_linop_applications: FEM regularized least squares; rows for CQRRTO, the CholQR family, published and refined Blendenpik, and an unpreconditioned reference; multi-run repeats; per-round records, a timing breakdown, and a sketched Karlson-Walden backward-error sidecar.toeplitz_ls_benchmark: the same roster on the matrix-free augmented Toeplitz problem.cqrrto_bench_common.hh, echo host, argv, commit and every environment knob into the CSV header, and verify header and row field counts after writing.--be-tol-multstops a row once the sketched Karlson-Walden estimate of the iterate is belowmult * sqrt(n) * u * ||A||_F, checked once per round with the oracle's time kept out of every reported solve time; the estimate's sparse-sign sketch is scaled by its isometry factor and pinned at eight nonzeros per column.Tests
CQRRTdense andCQRRTO_linopsoperator paths,IterRefineLSQand the restarted engine (including the oracle exits and their precedence), LSQR and Blendenpik, the refined Blendenpik dispatch, every new operator, the thread cap, the generator, and shared-pattern sparseaxpby.Removed (relative to
main)rl_cqrrt_linops.hh: its operator-based driver now lives inrl_cqrrt.hhasCQRRTO_linops.CQRRT_linop_composite_applications: the generalized SVD / generalized LS benchmark; itsL^{-1} K Voperator is kept underextrasand covered by the operator tests.