CQRRT diagnostics - #129
Open
mmelnich wants to merge 71 commits into
Open
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.
The regularization was mu = mu_factor * u(precond) -- with no ||A|| or size
scaling it is ~mn times too small past kappa(A) ~ 1/sqrt(u), so A^T A + mu^2 I
stays numerically singular at high kappa and even the double/double cell fails
(observed on the 2026-06-15 FEM2 run). Replace with the textbook shifted-
CholeskyQR shift applied through the augmented operator:
mu = mu_factor * ||A||_2 * sqrt((mn + n^2) * u(precond))
so mu^2 is the Fukaya shift (~ ||A||^2 * u * mn), keeping the Gram PD for
kappa(A) up to ~1/u. mu_factor ~ sqrt(11) ~ 3.3 is the principled value; lower it
to reproduce the under-regularized "mu = 10u" regime on purpose. ||A||_2 (already
estimated for the Higham metric) is hoisted above the mu computation and reused.
Also: usage strings now list the irlsq_reg mode.
The A2norm/sqrt(mn)-scaled "Fukaya" shift (e2eb9fb) was the wrong objective: it sizes mu for Q-factor ORTHOGONALITY, but IR-LSQ needs R to be a good right PRECONDITIONER, which requires mu <~ sigma_min. The Fukaya mu (~8e-4*||A|| at FEM2 scale) is orders above sigma_min, so the preconditioned CG matrix M = R^-T A^T A R^-1 had kappa(M) ~ (mu/sigma_min)^2 ~ 1e7+ and CG stalled at the iteration cap (fwd_err ~ 0.9, all methods bit-identical). Restore the collaborator's spec exactly: mu = mu_factor * u(precond), mu_factor=10 => mu = 10u, no ||A|| or size scaling. The augmented operator A_hat = [A; mu*I], its Q-less CholeskyQR R, and R-as-right-preconditioner are unchanged and match the spec. ||A||_2 is still estimated for the Higham metric only.
…olesky Pass the same unbounded eps*trace(G) shift-retry to cholqr_primitive that CholQR/CholQR2 use (unshifted first attempt, max_retries=-1, growth 10), so an ill-conditioned or single-precision Gram is rescued instead of failing at potrf. Double precision is unchanged (shift never activates).
…nt reporting New rl_lsqr.hh (matrix-free Paige-Saunders LSQR, optional right preconditioner) and rl_blendenpik.hh (sparse SASO sketch -> Householder QR -> LSQR); wired into the irlsq_reg benchmark as method bit 32. cholqr_primitive/cholqr_iterate now report the adaptive-shift retry count, surfaced per method (n_chol_retries) and in a new chol_retries CSV column.
Instrument IterRefineLSQ so a capped inner solve is distinguishable from a converged one (it previously returned success either way): per-step InnerCGStatus, achieved/best relative residual and best iteration, surfaced as new CSV columns plus cond(J R^-1). Expose max_inner and inner_tol on the CLI; the 400-iteration ceiling was max_inner=200 times two outer steps. Blendenpik: add the sketch-and-solve initialization required for forward stability (arXiv:2406.03468 sec. 3.1), applied as an equivalent shift outside rl_lsqr so the five methods sharing it are unaffected; equalize its iteration budget with the other methods; report convergence and solver residual instead of a sentinel. Correctness: bound the previously unbounded Cholesky retry loop and guard against non-finite shifts; free the owned Q buffer before overwriting it (five sites); replace hard-coded times[] total indices with total_us(); stop reading C when beta == 0 in ToeplitzLinOp. Convert seven rank-1 trsm calls to trsv. Export rl_lsqr.hh and rl_blendenpik.hh from RandLAPACK.hh and add their first unit tests. Suite: 336/336.
The matrix built to break CholeskyQR was in fact the most benign input possible: gen_bad_cholqr_singvals returned all ones, so the condition number was 1 for every value of cond requested. Three faults compounded. 'int offset = k' made the decay loop 'for (i = offset; i < k; ++i)' empty. The rate log(1e8/cond)/(1-(n-offset)) was written for a block of n-offset entries while the loop wrote into a length-k vector. And the signature could not express the intent: the dispatcher passes info.rank as k, and rank defaults to n, so both arguments are equal by default and neither could supply the count of leading ones. There were no callers outside gen_bad_cholqr_mat, no test and no benchmark, which is why it survived. The unusable second dimension argument is dropped. The leading half of the spectrum is one; the trailing half drops to 1e-8 and decays geometrically to 1/cond, giving condition number exactly cond. Guards throw for k < 2 and for cond < 1e8, below which the trailing block would rise rather than decay, leaving a non-monotone spectrum whose condition number is 1e8 rather than the request. That threshold is also where the modeled failure begins, since unshifted CholeskyQR loses orthogonality past eps^(-1/2). Verified: condition number exact at 1e8, 1e10 and 1e12, monotone in each case. End to end through mat_type::bad_cholqr, CholQR now loses orthogonality entirely (9.86e-01 at cond 1e8) while CholQR2 and sCholQR3 degrade gracefully and sCholQR3 needs no shift retry. Note the half-and-half split between the two blocks is a choice; nothing in the original pinned it down, since the parameter that would have carried it was unusable.
CholQR_dense, CholQR2_dense and sCholQR3_dense take a raw column-major buffer
instead of a LinearOperator, for callers working through the regular BLAS API.
They add no numerics of their own: each wraps its input in a linops::DenseLinOp
and delegates to the shared cholqr_iterate engine, so the pass count, the shift
policy and the adaptive-shift retry are inherited from the LinOp drivers rather
than reimplemented. sCholQR3_dense carries the same eps shift on iteration 1 and
zero on iterations 2 and 3.
A is const and is not modified. R is the output factor. Q = A R^{-1} is
materialized only when a buffer is supplied, outside the timing region, matching
the LinOp drivers' test mode.
Validated against the expected stability ordering on a rotated spectrum at three
condition numbers. Note that a column-scaled generator is unsuitable here:
CholeskyQR factors column scaling out, so the family does not separate and every
method reports machine precision regardless of the nominal condition number.
rl_lsqr already computed times[0..2] but the benchmark recorded only the total, leaving solve_time_us uninterpretable: the 07-29 campaign showed per-iteration solve cost falling as the circulant FFT length grew, which no cost model allows. Emitting the parts separates operator cost from preconditioner cost from LSQR's own vector work and setup. Also corrects the applications-benchmark usage comment: the accepted modes are irlsq_reg|rspec, and an unknown mode is silently undiagnosed.
The inner CG could reach its residual floor early and then grind to the iteration cap with no progress, reporting success either way. On the FEM2 operator at kappa^colnorm=1e10 the CholQR preconditioner is unusable (cond(J R^-1)=7.8e4 against ~1.000 for the other methods) and its CG hit bottom at iteration 17 of each 200-iteration step. A paired diagnostic run settled the mechanism and ruled out the obvious fix: cap 200/step: best_relres 3.483414e-09 @ iter 17, solution error 48.7, 13 s cap 2000/step: best_relres 3.483414e-09 @ iter 17, solution error 547.0, 166 s Ten times the budget left the best residual bit-identical and made the outer solution 11x worse, so the cap was never the binding constraint and raising it is a regression. The last iterate is also worse than the best one, which is the mechanism behind that degradation. Adds InnerCGStatus::Stagnated, an inner_stag_window / inner_stag_rel_improve pair (default: no 0.1% drop for 20 consecutive iterations), and a best-iterate snapshot returned on both Stagnated and HitCap exits. Convergence is tested before stagnation, so a converging solve is unaffected; window 0 disables. Also fixes record_inner_cg_diagnosis, which selected the worst outer step by numeric enum value. Stagnated=3 would have outranked and masked a genuine Breakdown=2; it now ranks by explicit severity.
…er_restarts, default 1) After the inner CG terminates (Converged/Stagnated/HitCap), rerun it once from the iterate it returned, recomputing the TRUE residual c - M z (one extra M apply). CG's recursive residual drifts from the true residual in finite precision, and both the convergence test and the stagnation window read the recursive one; the restart discards the drifted recurrence and gives fresh conjugacy. A genuinely converged solve pays only the entry check (0 iterations). New apply_M helper shared by the iteration body and the warm entry so the two can never compute M differently. Per-step diagnostics aggregate attempts: iters summed, best tracked across attempts, status/relres from the final one. Tests: the two single-attempt-mechanics tests now pin inner_restarts = 0; new test pins the summed count (cap 2 -> exactly 4), the never-lose-ground best residual, and the near-free converged restart. Suite 338/338.
…e order-independent RSS is process-cumulative (glibc keeps freed arenas mapped), so the FIRST tracked method in a benchmark absorbed the whole process ramp-up while later methods reused already-faulted pages and reported ~0 (the 2026-07-29 peak_rss_kb=4 effect, and the inflated CQRRT_linop storage bars in the Toeplitz figures: measured 471 MB vs 184 MB analytical, only because it ran first). Trimming freed heap before the baseline makes each method measure from live memory only.
… num_runs repeats Policy (2026-08-05): the sketch-and-solve x0 warm start is Blendenpik-only. Method mask bit 32 now emits TWO rows, Blendenpik (warm) and Blendenpik_cold (x0 = 0); every Q-less QR method is unconditionally cold. The qless/IR warm CLI knobs are removed; CQRRT_linop_applications REJECTS the old 15/16-arg invocations so stale job scripts fail loudly. Timing hygiene (single-run solves at 4-7 iterations sit at the wall-clock noise floor): untimed CPU warmups before any measurement (toeplitz had none; irlsq_reg/sparse modes warmed only the QR build, never the solve chain), and the Toeplitz benchmark gains [num_runs] with a run CSV column and per-row flush; the MATLAB plotters aggregate (best run by default).
Closes the 2026-07-14 deferred item: the Toeplitz reference benchmark's second solver (restarted PCG on the right-preconditioned normal equations) now exists as RandLAPACK::restarted_pcg_ne and is the Toeplitz benchmark's DEFAULT solver (deliberate divergence from the reference's lsqr default; pass "lsqr" to reproduce the 08-05-era campaigns). Both least-squares drivers now run the SAME inner solver: new rl_pcg_inner.hh holds the instrumented CG kernel (stagnation window + best-iterate return from 07-29, true-residual warm-start entry from 08-05), extracted verbatim from IterRefineLSQ and shared with restarted_pcg_ne. Stagnation in a pcg_ne round is not terminal (the next true-residual round decides); breakdown still is. Outer contracts aligned: IterRefineLSQ gains outer_tol (stop refining once ||b - Jx||/||b|| meets it, capped at n_refine_steps; 0 = historical fixed-step behaviour), and both benchmarks default to 3 TOTAL outer rounds. Accuracy rationale (2026-08-06 controlled probe): tighter inner tolerances and inner restarts do NOT improve final accuracy (bit-identical across four decades of tolerance); outer refinement steps do. From x0 = 0, 3 outer steps reach machine-precision backward error at kappa ~ 5e10 where the previous hardcoded 2 steps plateaued at ~5e-13. Epperly et al. (arXiv:2406.03468) Alg. 1 proves 2 steps suffice only WITH sketch-and-solve initialization, which we deliberately do not use (collaborator request, 06-09). Benchmark CLI: Toeplitz gains [solver] [pcg_restart_maxit] [pcg_max_restarts] + solver/pcg_restarts CSV columns (appended; name-based MATLAB readers are unaffected). CQRRT_linop_applications gains [ir_inner_restarts] (default 1) [ir_n_steps] (default 3) [ir_outer_tol] (default 10*eps of solve precision); Blendenpik's equalized budget scales with n_steps; the stale warm-start script guard moves past the new positions. Tests: 7 new (restarted_pcg_ne correctness/restart/cap/stagnation semantics, IterRefineLSQ outer_tol early exit); suite 345/345.
The reference MATLAB recomputes each round's normal-equation residual as
g - H z, a difference of two large kappa-contaminated quantities. Replaced with
the mathematically identical stable form R^{-T}(A^T(b - A x)), mapping the SMALL
true LS residual (Epperly et al. arXiv:2406.03468, Alg. 1 line 5). Second
deliberate deviation from the reference, and a measured one: A/B on the m=800
prolate benchmark case (FFT operator, lambda_rel 1e-20, everything else
identical) has g - H z stalling every method's LS relres at 1.75e-6 with
recovery error 1.75e4, while this form reaches the 1e-10 noise floor with
recovery 1.9e-3, matching warm Blendenpik. CQRRT does it in 4 inner iterations.
Two new tests: dense column-scaled kappa ~ 1e11 and dense augmented prolate
(noisy, sqrt(lambda) = 1e-10 rows, shifted-CholQR factor) both converge -- note
they do NOT discriminate the residual forms (dense replicas lack the FFT
apply's rounding, which is a necessary ingredient); the benchmark case is the
discriminating harness, documented in the driver. Suite 347/347.
Seven paths now: [2] solve R_sk*X=I via Side::Left (the shipping PCholQRPrecondMethod::TRSM_IDENTITY default, with laset parity) and [3] the old Side::Right reversed ordering, kept to document the solve-ordering effect. Part B step-divergence updated to match the current default. CSV columns renumbered; path names embedded as header comments. Measured on photogrammetry2 (d=4n, nnz=4, medians of 5): default 1.2e-9, reversed 1.1e-4, TRTRI 1.2e-4.
restarted_pcg_ne is now the single engine behind both benchmarks: restart_drop default 1e-2 -> 1e-4 (Oleg's pacing); inner_abs_tol guard so rounds stop once below the absolute target; per-round PCGRoundHistory diagnostics with inner/outer time split; outer stagnation exit (2 rounds without true-residual improvement) so methods at a noise floor stop instead of grinding the round cap. IterRefineLSQ is a thin adapter over the engine: inner_restarts removed (every round is a true-residual restart), round_drop added (0 = legacy fixed-tol rounds), always cold-started; a new test pins bitwise-identical solutions between adapter and engine. Round caps 20 in both benchmarks (FEM2 ir_n_steps, Toeplitz pcg_max_restarts; the old 3-4 rounds truncated the unpreconditioned baseline five orders above tol). FEM2 CLI slot 13 is ir_round_drop now; values >= 1 rejected so stale scripts fail loudly. Suite 347 -> 348. Local m=800 prolate smoke: CQRRT 4 inner iters, CholQR 39, unprec 278, all at the 1e-10 data noise floor.
MKL's threaded dtrsv collapses on the small dense factors these solvers apply: measured on an n=2000 factor, 0.448 ms at 1 thread, 0.162 ms at 4, and 31.3 ms at 16. It is not a triangularity problem, so blocking does not fix it (a plain dgemv of the same size degrades the same way, and a hand-blocked solve still measured 15.1 ms at 16 threads), and it does not improve with size: at 16 threads dtrsv costs 46/84/162 ms for n=2000/4000/8000 against 0.16/1.15/4.16 ms at 4 threads. The preconditioned solvers apply the preconditioner twice per inner iteration while the unpreconditioned baseline applies it zero times, so the overhead taxed exactly the methods that converge in few iterations and inverted the wall-clock ranking against unpreconditioned. New Blas2ThreadGuard (RAII, MKL-only, no-op elsewhere) wraps all nine trsv sites in the narrowest scope, leaving the surrounding operator applies at full width. Cap defaults to 4 and is overridable at runtime via RANDLAPACK_BLAS2_THREADS. Includes mkl_service.h only; the umbrella mkl.h redeclares LAPACK at MKL's integer width and collides with LAPACK++. Measured on 16000x2000 prolate-Toeplitz at full threads, guard off then on: CQRRT solve 1531 -> 22.8 ms, CholQR 8293 -> 115 ms, unpreconditioned 106 -> 90 ms. Suite 348/348.
Measured dtrsv on a Xeon Gold 6430 (64 cores, exclusive) rather than on
a desktop: threading helps up to 8-16 threads and degrades past that,
and the optimum moves with n.
threads 1 4 8 16 32 64
n = 2000 0.562 0.408 0.392 0.774 0.909 0.887
n = 8256 18.147 6.232 4.915 4.735 5.828 7.130
n = 20000 108.209 45.121 27.110 19.276 19.564 24.354
So the cap is a peak-seeker, now two-tier (8 for n <= 4000, 16 above)
and passed the factor dimension, rather than the flat 4 taken from a
desktop. It also corrects the expected magnitude: unguarded costs
1.5-2.3x here, not the 100x seen on an oversubscribed WSL2 box.
The same job confirms a tall gemv scales cleanly to 64 threads
(225 -> 23 -> 21 ms), so rl_determiter/KRILL is deliberately left
unguarded: one barrier per call over a large operand is a different
regime from one barrier per column.
Suite 348/348.
MKL's threaded FFT intermittently stalls at full width. A stage-resolved
probe caught single DftiComputeForward calls taking 16-32 ms instead of
~0.1 ms, with ~99% of the stall inside the FFT call itself while fill,
pointwise multiply and the backward transform stayed at a steady
80-100 us.
A solver converging in a handful of applies cannot average those stalls
out, so the artifact scaled INVERSELY with preconditioner quality and
inverted the wall-clock ranking. Measured on 16000x2000 prolate-Toeplitz,
3 repeats each:
threads CQRRT solve unprec solve
8 19.5 / 18.8 / 18.5 ms 200 / 200 / 199 ms
16 26.4 / 26.4 / 25.1 ms 177 / 172 / 174 ms
64 344 / 342 / 359 ms 303 / 313 / 413 ms
At 64 the two look equal; at 8 preconditioning wins 10.6x, which is what
the iteration counts (5 vs 267) always implied.
Capping is the vendor-recommended remedy for single small transforms
(batching via DFTI_NUMBER_OF_TRANSFORMS is the alternative and is
unavailable: CG produces one right-hand side at a time). The operator
sweep shows nothing is lost -- 64 threads buys no speed at these
transform sizes, only stall risk.
Guard generalized with an explicit-cap constructor; default 16,
overridable via RANDLAPACK_FFT_THREADS. Applied at all three DFTI
compute sites (the only FFT user in the repo). Suite 348/348.
As published, Blendenpik is sketch + QR + LSQR with no refinement, so the
suite compared its preconditioner through a different solver than every
Q-less method uses -- conflating preconditioner quality with solver
structure. The published rows are untouched; two rows are added that hand
Blendenpik's own R and its own answer to the shared restarted-PCG engine,
started from its warm and cold solutions respectively.
restarted_pcg_ne gains an optional x0: z is seeded with R*x0 so that
x = R^{-1} z reproduces it, and the first round's normal-equation residual
comes from the true residual b - A x0 in the stable form (g - H z would
reintroduce the cancellation removed on 2026-08-06). IterRefineLSQ exposes
it as warm_x0; every Q-less path still cold-starts.
New mask bits: 128 (Toeplitz), 64 (FEM2). Generators updated to 255 / 127.
Measured on 16000x2000 prolate-Toeplitz -- refinement fully rescues the
cold case, and the preconditioner was never at fault (cond 3.21 in all
four rows):
row solver relres recovery error
Blendenpik (warm) 1.005e-10 7.440e-04
Blendenpik_cold 3.834e-07 2.741e+03
Blendenpik_refine 1.005e-10 7.440e-04
Blendenpik_cold_refine 1.005e-10 7.441e-04
Suite 348 -> 349 (new warm-start test pins: converged x0 costs zero inner
iterations, a perturbed x0 is refined back to the reference, warm and cold
reach the same fixed point).
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)
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.
Adds
CQRRT_orth_gap.cc— a standalone benchmark that compares the orthogonality quality of CQRRT_linop and CQRRT_expl on the same input matrix. Five input modes are supported:generate— synthetic sparse matrix with controlled condition number and densityfile— any Matrix Market (.mtx) file from diskcomposite— composite operator from aK.mtx+V.mtxpair (buildsCholSolver(K) ∘ Sparse(V) = L⁻¹V), for generalized LS / FEM-type problemsdiag— step-by-step diagnostic on a Matrix Market file, reimplements CQRRT manually to pinpoint where numerical divergence between the two paths arisesdiag_gen— same step-by-step diagnostic on a synthetic matrixUsage examples