diff --git a/RandLAPACK.hh b/RandLAPACK.hh index c9741e3ee..8a0ddb33b 100644 --- a/RandLAPACK.hh +++ b/RandLAPACK.hh @@ -18,6 +18,7 @@ #include "RandLAPACK/testing/rl_test_utils.hh" // Computational routines +#include "RandLAPACK/comps/rl_matfun.hh" #include "RandLAPACK/comps/rl_determiter.hh" #include "RandLAPACK/comps/rl_preconditioners.hh" #include "RandLAPACK/comps/rl_qb.hh" @@ -27,6 +28,9 @@ #include "RandLAPACK/comps/rl_syrf.hh" #include "RandLAPACK/comps/rl_orth.hh" #include "RandLAPACK/comps/rl_rpchol.hh" +#include "RandLAPACK/comps/rl_lanczos_fa.hh" +#include "RandLAPACK/comps/rl_lanczos_fa_block.hh" +#include "RandLAPACK/comps/rl_hutchinson.hh" // Drivers #include "RandLAPACK/drivers/rl_rsvd.hh" @@ -36,9 +40,10 @@ #include "RandLAPACK/drivers/rl_scholqr3_linops.hh" #include "RandLAPACK/drivers/rl_cqrrpt.hh" #include "RandLAPACK/drivers/rl_bqrrp.hh" -#include "RandLAPACK/drivers/rl_revd2.hh" +#include "RandLAPACK/drivers/rl_nystrom_evd.hh" #include "RandLAPACK/drivers/rl_abrik.hh" #include "RandLAPACK/drivers/rl_krill.hh" +#include "RandLAPACK/drivers/rl_fun_nystrom_pp.hh" // Cuda functions - issues with linking/visibility when present if the below is uncommented. // A temporary fix is to add the below directly in the test/benchmark files. diff --git a/RandLAPACK/CMakeLists.txt b/RandLAPACK/CMakeLists.txt index 7d70a1d17..434f8edd4 100644 --- a/RandLAPACK/CMakeLists.txt +++ b/RandLAPACK/CMakeLists.txt @@ -6,7 +6,7 @@ set(RandLAPACK_cxx_sources rl_cqrrpt.hh rl_bqrrp.hh rl_rsvd.hh - rl_revd2.hh + rl_nystrom_evd.hh rl_qb.hh rl_orth.hh rl_util.hh diff --git a/RandLAPACK/comps/rl_hutchinson.hh b/RandLAPACK/comps/rl_hutchinson.hh new file mode 100644 index 000000000..b0791c5d9 --- /dev/null +++ b/RandLAPACK/comps/rl_hutchinson.hh @@ -0,0 +1,81 @@ +#pragma once + +#include "rl_blaspp.hh" +#include "rl_linops.hh" +#include "rl_util.hh" + +#include +#include +#include + +namespace RandLAPACK { + + +/// Hutchinson stochastic trace estimator. +/// +/// Estimates tr(M) using the identity E[ω^T M ω] = tr(M) for zero-mean +/// unit-variance ω. Draws s independent Rademacher vectors (iid Unif{±1}), +/// applies M to all at once, and returns (1/s) * <Ω, Z>_F where Z = M*Ω. +/// +/// M must satisfy linops::SymmetricLinearOperator. +/// +/// @tparam T Floating-point scalar type. +/// @tparam RNG Random number generator type. +template +class Hutchinson { +public: + // Internal buffers — grown with new/delete[], never shrunk. + T* Omega = nullptr; int64_t Omega_sz = 0; + T* Z = nullptr; int64_t Z_sz = 0; + + Hutchinson() = default; + Hutchinson(const Hutchinson&) = delete; + Hutchinson& operator=(const Hutchinson&) = delete; + + ~Hutchinson() { delete[] Omega; delete[] Z; } + + // ------------------------------------------------------------------ + /// Low-level estimator: given precomputed Ω (n×s) and Z = M*Ω (n×s), + /// returns <Ω, Z>_F / s = (1/s) * Σ_j ω_j^T (M ω_j). + /// Uses blas::dot on flattened n*s arrays — no allocation. + /// + /// @param[in] Omega_buf n×s sketch matrix (column-major). + /// @param[in] Z n×s result of M applied to Omega_buf (column-major). + /// @param[in] n Ambient dimension. + /// @param[in] s Number of samples. + /// @returns Frobenius inner product <Ω, Z>_F / s. + T estimate(const T* Omega_buf, const T* Z, int64_t n, int64_t s) const { + // Frobenius inner product of two n×s matrices, treated as flat n*s vectors + return blas::dot(n * s, Omega_buf, 1, Z, 1) / static_cast(s); + } + + // ------------------------------------------------------------------ + /// High-level estimator: draws Ω internally, applies M, returns trace estimate. + /// n is taken from M.dim. + /// + /// @param[in] M Operator satisfying SymmetricLinearOperator. + /// @param[in] s Number of Hutchinson samples. + /// @param[in] state RandBLAS RNG state; advanced on return. + template + T call(SLO& M, int64_t s, RandBLAS::RNGState& state) { + int64_t n = M.dim; + + util::resize(Omega, Omega_sz, n * s); + + // Draw Ω with iid Rademacher entries (Unif{±1}). + // RandBLAS has no ScalarDist::Rademacher, but ScalarDist::Uniform fills + // with Unif[-1,1] via r123ext::uneg11; sign-transforming gives exact ±1. + RandBLAS::DenseDist D(n, s, RandBLAS::ScalarDist::Uniform); + state = RandBLAS::fill_dense(D, Omega, state); + for (int64_t i = 0; i < n * s; ++i) + Omega[i] = (Omega[i] >= 0) ? (T)1 : (T)-1; + + util::resize(Z, Z_sz, n * s); + M(Layout::ColMajor, s, (T)1.0, Omega, n, (T)0.0, Z, n); + + return estimate(Omega, Z, n, s); + } +}; + + +} // end namespace RandLAPACK diff --git a/RandLAPACK/comps/rl_lanczos_fa.hh b/RandLAPACK/comps/rl_lanczos_fa.hh new file mode 100644 index 000000000..9a8128e16 --- /dev/null +++ b/RandLAPACK/comps/rl_lanczos_fa.hh @@ -0,0 +1,278 @@ +#pragma once + +#include "rl_blaspp.hh" +#include "rl_lapackpp.hh" +#include "rl_linops.hh" +#include "rl_util.hh" + +#include +#include +#include +#include +#include +#include + +#ifdef _OPENMP +#include +#endif + +namespace RandLAPACK { + + +/// d-step block Lanczos for matrix function application f(A)B. +/// Approximates f(A)B column-wise via independent Krylov subspaces of dimension d. +/// See: T. Chen, "A Lanczos-FA algorithm for matrix function approximation" (2022). +/// +/// @tparam T Floating-point scalar type. +template +class LanczosFA { +public: + /// Reorthogonalization control. + /// 1 = full (project out all previous Krylov vectors after each step). + /// 0 = none (vanilla Lanczos, per Persson's reference implementation). + /// Lanczos-FA tolerates loss of orthogonality better than eigenvalue + /// Lanczos (Paige-Greenbaum theory), so vanilla often works in practice. + /// Full reorthogonalization is the safe default for a numerical library. + int64_t reorth = 1; + + // Internal buffers — grown with new/delete[], never shrunk between calls. + // Dimension key: n = operator dimension, s = number of RHS vectors (columns of B), + // d = number of Lanczos steps. + // + // K: (d+1) × n × s — Krylov basis blocks. + // Layout: K[step * n*s + col * n + row] = row-th entry of step-th basis vector for column col. + // Storing steps as contiguous n×s slices keeps each batch matvec contiguous, + // while the per-column stride (n*s) lets apply() use strided gemv for reconstruction. + // alpha: s × d — tridiagonal diagonals, alpha[j*d + i] = α_{i,j}. + // beta: s × (d-1) — tridiagonal subdiagonals, beta[j*(d-1) + i] = β_{i+1,j}. + // lapack::stevd expects the diagonal and subdiagonal as separate arrays, + // so alpha and beta are stored separately rather than interleaved. + // normb: s — column norms of B before normalization. + T* K = nullptr; int64_t K_sz = 0; + T* alpha = nullptr; int64_t alpha_sz = 0; + T* beta = nullptr; int64_t beta_sz = 0; + T* normb = nullptr; int64_t normb_sz = 0; + T* workspace = nullptr; int64_t workspace_sz = 0; + + bool timing = false; + std::vector times; // populated after call() when timing==true + // Slots: matvec, run_lanczos, apply_f, rest, total + long _t_matvec_us = 0; // accumulated in run_lanczos() when timing==true + + LanczosFA() = default; + LanczosFA(const LanczosFA&) = delete; + LanczosFA& operator=(const LanczosFA&) = delete; + + ~LanczosFA() { delete[] K; delete[] alpha; delete[] beta; delete[] normb; delete[] workspace; } + + // ------------------------------------------------------------------ + /// Run the d-step block Lanczos recurrence on B. + /// Fills K, alpha, beta, normb from B (n×s column-major). + /// Calls A exactly d times, each application to an n×s matrix. + /// + /// @param[in] A SymmetricLinearOperator — matvec oracle. + /// @param[in] B n×s input matrix (column-major); not modified. + /// @param[in] n Dimension of A. + /// @param[in] s Number of right-hand sides (Hutchinson samples). + /// @param[in] d Number of Lanczos steps. + template + void run_lanczos(SLO& A, const T* B, int64_t n, int64_t s, int64_t d) { + using namespace std::chrono; + steady_clock::time_point _mv_t0, _mv_t1; + _t_matvec_us = 0; + + // Grow buffers if needed + util::resize(K, K_sz, (d + 1) * n * s); + util::resize(alpha, alpha_sz, d * s); + if (d > 1) util::resize(beta, beta_sz, (d - 1) * s); + util::resize(normb, normb_sz, s); + + // Step 0: q_1 = column-normalize B; store in K[:,:,0] + T* K0 = K; + lapack::lacpy(lapack::MatrixType::General, n, s, B, n, K0, n); +#pragma omp parallel for schedule(static) + for (int64_t j = 0; j < s; ++j) { + T nrm = blas::nrm2(n, K0 + j * n, 1); + normb[j] = nrm; + // Zero input column: skip normalization; normb[j]=0 so apply_f outputs zero. + if (nrm > (T)0) + blas::scal(n, (T)1.0 / nrm, K0 + j * n, 1); + } + + // Step 0 matvec: K[:,:,1] = A * K[:,:,0] + T* K1 = K + n * s; + if (this->timing) _mv_t0 = steady_clock::now(); + A(Layout::ColMajor, s, (T)1.0, K0, n, (T)0.0, K1, n); + if (this->timing) { _mv_t1 = steady_clock::now(); _t_matvec_us += duration_cast(_mv_t1 - _mv_t0).count(); } + + // α[0, j] = q_1[:,j] · (A q_1)[:,j] — s independent inner products, + // one tridiagonal diagonal entry per column. +#pragma omp parallel for schedule(static) + for (int64_t j = 0; j < s; ++j) + alpha[j * d + 0] = blas::dot(n, K1 + j * n, 1, K0 + j * n, 1); + + // Main Lanczos loop: steps 1..d-1 + // At the start of iteration i: + // K[:,:,i] = A*q_i - β_i*q_{i-1} (partial three-term, β part done last iter) + // K[:,:,i+1] is free workspace + // This iteration: + // (1) subtract α_i*q_i to complete three-term → K_{i+1} becomes unnormalized q_{i+1} + // (2) optional reorthogonalization + // (3) β_{i+1} = ||K_{i+1}||, (4) normalize → q_{i+1} + // (5) K_{i+2} = A*q_{i+1} - β_{i+1}*q_i (start of next three-term) + // (6) α_{i+1} = q_{i+1}[:,j] · K_{i+2}[:,j] + // All per-column loops below are independent across j and parallelized. + for (int64_t i = 0; i < d - 1; ++i) { + T* K_prev = K + i * n * s; // K[:,:,i] = q_i (normalized) + T* K_curr = K + (i + 1) * n * s; // K[:,:,i+1] = partial (A*q_i - β_i*q_{i-1}) + T* K_new = K + (i + 2) * n * s; // K[:,:,i+2] = workspace for A*q_{i+1} + + // (1) Complete three-term: K_curr -= α_i * K_prev + // Each column has a different scalar α_{i,j}, so GEMM would cost O(n·s²); axpy is optimal. +#pragma omp parallel for schedule(static) + for (int64_t j = 0; j < s; ++j) + blas::axpy(n, -alpha[j * d + i], K_prev + j * n, 1, K_curr + j * n, 1); + + // (2) Optional full reorthogonalization: project K_curr[:,j] out of all q_0..q_i. + // Outer loop over j is parallel (columns are independent); inner prev-loop is + // sequential per column (each projection modifies K_curr[:,j] in place). + int64_t reorth_steps = reorth ? (i + 1) : 0; +#pragma omp parallel for schedule(static) + for (int64_t j = 0; j < s; ++j) { + for (int64_t prev = 0; prev < reorth_steps; ++prev) { + T* K_p = K + prev * n * s; + T coeff = blas::dot(n, K_curr + j * n, 1, K_p + j * n, 1); + blas::axpy(n, -coeff, K_p + j * n, 1, K_curr + j * n, 1); + } + } + + // (3) β_{i+1} = column norms, (4) normalize → q_{i+1} + // Zero norm means the Krylov basis has collapsed for that column. + // Store β=0 (the tridiagonal subdiagonal entry) and skip normalization; + // stevd handles a zero subdiagonal correctly (independent 1×1 blocks). +#pragma omp parallel for schedule(static) + for (int64_t j = 0; j < s; ++j) { + T nrm = blas::nrm2(n, K_curr + j * n, 1); + beta[j * (d - 1) + i] = nrm; + if (nrm > (T)0) + blas::scal(n, (T)1.0 / nrm, K_curr + j * n, 1); + } + + // (5) K_new = A*q_{i+1} - β_{i+1}*q_i + // Different β per column — same reasoning as (1), axpy is optimal. + if (this->timing) _mv_t0 = steady_clock::now(); + A(Layout::ColMajor, s, (T)1.0, K_curr, n, (T)0.0, K_new, n); + if (this->timing) { _mv_t1 = steady_clock::now(); _t_matvec_us += duration_cast(_mv_t1 - _mv_t0).count(); } +#pragma omp parallel for schedule(static) + for (int64_t j = 0; j < s; ++j) + blas::axpy(n, -beta[j * (d - 1) + i], K_prev + j * n, 1, K_new + j * n, 1); + + // (6) α_{i+1} = q_{i+1}[:,j] · K_new[:,j] +#pragma omp parallel for schedule(static) + for (int64_t j = 0; j < s; ++j) + alpha[j * d + (i + 1)] = blas::dot(n, K_new + j * n, 1, K_curr + j * n, 1); + } + } + + // ------------------------------------------------------------------ + /// Evaluate f(A)B from precomputed Krylov data (K, alpha, beta, normb). + /// Per column j: eigendecompose T_j = S_j diag(θ_j) S_j^T via lapack::stev, then: + /// out[:,j] = normb[j] * Q_j * S_j * diag(f(θ_j)) * S_j[0,:]^T + /// where Q_j is the n×d Lanczos basis stored in K and S_j[0,:] is the first row + /// of the eigenvector matrix (Chen 2022, eq. 2.3). + /// Per-column stev calls are independent — parallelized with OpenMP. + /// + /// @tparam F Callable as T(T) — lambda, function pointer, or functor. + /// std::invocable (C++20) enforces this at the call site. + /// @param[in] f Instance of F applied elementwise to tridiagonal eigenvalues θ. + /// @param[in] n Dimension of A. + /// @param[in] s Number of right-hand sides. + /// @param[in] d Number of Lanczos steps (tridiagonal size). + /// @param[out] out n×s output matrix (column-major); overwritten. + template F> + void apply_f(F f, int64_t n, int64_t s, int64_t d, T* out) { + // Per-thread workspace: alpha_j(d) + beta_j(d-1) + Z_j(d*d) + c_j(d) + v_j(d) = d^2 + 4d - 1 + int64_t workspace_per_thread = d * d + 4 * d - 1; + int nthreads = 1; +#ifdef _OPENMP + nthreads = omp_get_max_threads(); +#endif + util::resize(workspace, workspace_sz, (int64_t)nthreads * workspace_per_thread); + +#pragma omp parallel for schedule(static) + for (int64_t j = 0; j < s; ++j) { + int tid = 0; +#ifdef _OPENMP + tid = omp_get_thread_num(); +#endif + T* base = workspace + tid * workspace_per_thread; + T* alpha_j = base; + T* beta_j = alpha_j + d; + T* Z_j = beta_j + std::max(d - 1, (int64_t)0); + T* c_j = Z_j + d * d; + T* v_j = c_j + d; + + // Copy per-column tridiagonal entries into per-thread workspace. + // stevd overwrites its alpha/beta arrays in-place (they become eigenvalues + // and workspace), so we must work on copies — otherwise the member alpha/beta + // would be destroyed and a second apply_f call (with a different f) would + // produce garbage without re-running run_lanczos. + // alpha[j*d .. j*d+d-1] and beta[j*(d-1) .. j*(d-1)+(d-2)] are flat 1D + // vectors of length d and d-1, not 2D matrices — blas::copy is correct here; + // lacpy is for 2D matrices with distinct leading dimensions. + blas::copy(d, alpha + j * d, 1, alpha_j, 1); + if (d > 1) + blas::copy(d - 1, beta + j * (d - 1), 1, beta_j, 1); + + // d×d tridiagonal eigendecomposition: T_j = Z_j * diag(θ) * Z_j^T + // alpha_j → eigenvalues θ (ascending); Z_j → eigenvectors (column-major) + lapack::stevd(lapack::Job::Vec, d, alpha_j, beta_j, Z_j, d); + + // c_j[i] = f(θ_i) * S_j[0, i] + // In column-major Z_j (d×d): entry (row=0, col=i) = Z_j[i*d + 0] + for (int64_t i = 0; i < d; ++i) + c_j[i] = f(alpha_j[i]) * Z_j[i * d + 0]; + + // v_j = Z_j * c_j (d×d matrix times d-vector) + blas::gemv(Layout::ColMajor, Op::NoTrans, d, d, + (T)1.0, Z_j, d, c_j, 1, (T)0.0, v_j, 1); + + // out[:,j] = normb[j] * Q_j * v_j + // Q_j is n×d with column stride n*s (strided view into K buffer) + blas::gemv(Layout::ColMajor, Op::NoTrans, n, d, + normb[j], K + j * n, n * s, v_j, 1, (T)0.0, out + j * n, 1); + } + } + + // ------------------------------------------------------------------ + /// Combined run + apply: compute f(A)B in one call. + /// + /// @param[in] A SymmetricLinearOperator. + /// @param[in] B n×s input matrix (column-major). + /// @param[in] n Dimension of A. + /// @param[in] s Number of right-hand sides. + /// @param[in] f Scalar function T→T. + /// @param[in] d Lanczos steps. + /// @param[out] out n×s output, overwritten with f(A)B approximation. + template F> + void call(SLO& A, const T* B, int64_t n, int64_t s, F f, int64_t d, T* out) { + using namespace std::chrono; + _t_matvec_us = 0; + steady_clock::time_point t_total_start, t_lanczos_end, t_end; + if (this->timing) t_total_start = steady_clock::now(); + run_lanczos(A, B, n, s, d); + if (this->timing) t_lanczos_end = steady_clock::now(); + apply_f(f, n, s, d, out); + if (this->timing) { + t_end = steady_clock::now(); + long total_us = duration_cast(t_end - t_total_start).count(); + long lanczos_us = duration_cast(t_lanczos_end - t_total_start).count(); + long apply_f_us = duration_cast(t_end - t_lanczos_end).count(); + long rest_us = total_us - lanczos_us - apply_f_us; + this->times = {_t_matvec_us, lanczos_us, apply_f_us, rest_us, total_us}; + } + } +}; + + +} // end namespace RandLAPACK diff --git a/RandLAPACK/comps/rl_lanczos_fa_block.hh b/RandLAPACK/comps/rl_lanczos_fa_block.hh new file mode 100644 index 000000000..ac6ab3602 --- /dev/null +++ b/RandLAPACK/comps/rl_lanczos_fa_block.hh @@ -0,0 +1,275 @@ +#pragma once + +#include "rl_blaspp.hh" +#include "rl_lapackpp.hh" +#include "rl_linops.hh" +#include "rl_util.hh" + +#include +#include +#include +#include +#include +#include +#include + +namespace RandLAPACK { + + +/// d-step block Lanczos for matrix function application f(A)B. +/// Builds a single joint block Krylov subspace using BLAS-3 throughout, +/// replacing the s independent scalar Lanczos sequences in LanczosFA. +/// +/// Algorithm reference: T. Chen, "A handbook for matrix-function-based Krylov methods", +/// arXiv:2410.11090 (2024). Algorithm 9.2 (block Lanczos recurrence), +/// Definition 9.6 (block Lanczos-FA). +/// +/// At each recurrence step all updates (alpha, beta, Z) are BLAS-3 (GEMM + GEQRF). +/// apply_f calls a single syevd on the (d*s)×(d*s) block tridiagonal instead of +/// s separate stevd calls of size d. +/// +/// Block Lanczos-FA formula: out ≈ Q_basis * f(T_k) * E₁ * R₀ +/// Q_basis = [Q₀|...|Q_{d-1}] (n×d*s), T_k = d*s×d*s block tridiagonal, +/// E₁ = first d*s×s columns of identity, B = Q₀*R₀ (initial QR). +/// +/// Known limitation (v1): deflation is not implemented. When the block Krylov +/// space fills before d steps (B_step develops near-zero singular values), accuracy +/// degrades but the algorithm does not crash. For problems where d*s approaches the +/// effective rank of A, use a smaller d or the scalar LanczosFA. +/// +/// @tparam T Floating-point scalar type. +template +class BlockLanczosFA { +public: + /// Reorthogonalization control. + /// 1 = full (project each new block Z out of all previous Krylov blocks). + /// 0 = none. + int64_t reorth = 1; + + // Internal buffers — grown with new/delete[], never shrunk between calls. + // Dimension key: n = operator dimension, s = block size, d = Lanczos steps. + // + // K_big: (d+1)*n*s — block Krylov basis + matvec scratch. + // Layout: K_big[step*n*s .. (step+1)*n*s-1] = Q_step (n×s col-major, ld=n). + // First d blocks form Q_basis (n×d*s col-major, ld=n) used by apply_f. + // Block d is scratch for the current-step matvec output. + // R0_buf: s*s — upper triangular factor from initial QR of B. + // tau_buf: n — Householder scalars; reused at each geqrf/orgqr call. + // A_blk: d*s*s — block alphas (s×s symmetric), A_blk[step*s*s] = A_step. + // B_blk: d*s*s — block betas (s×s upper triangular), B_blk[step*s*s] = B_step. + // Only d-1 entries are populated (no beta after the last step); + // size d is allocated to keep indexing uniform. + // workspace: apply_f scratch — T_dense (d*s×d*s) + eig_vals (d*s) + G (d*s×s) + C1 (d*s×s). + // proj_buf: s*s — reorthogonalization scratch (Q_p^T * Y projection); reused across steps. + T* K_big = nullptr; int64_t K_big_sz = 0; + T* R0_buf = nullptr; int64_t R0_sz = 0; + T* tau_buf = nullptr; int64_t tau_buf_sz = 0; + T* A_blk = nullptr; int64_t A_blk_sz = 0; + T* B_blk = nullptr; int64_t B_blk_sz = 0; + T* workspace = nullptr; int64_t workspace_sz = 0; + T* proj_buf = nullptr; int64_t proj_buf_sz = 0; + + bool timing = false; + std::vector times; + long _t_matvec_us = 0; + + BlockLanczosFA() = default; + BlockLanczosFA(const BlockLanczosFA&) = delete; + BlockLanczosFA& operator=(const BlockLanczosFA&) = delete; + + ~BlockLanczosFA() { + delete[] K_big; delete[] R0_buf; delete[] tau_buf; + delete[] A_blk; delete[] B_blk; delete[] workspace; delete[] proj_buf; + } + + // ------------------------------------------------------------------ + /// Run the d-step block Lanczos recurrence on B (n×s col-major). + /// Fills K_big, R0_buf, A_blk, B_blk. + /// Calls A exactly d times, each applied to an n×s block. + template + void run_lanczos(SLO& A, const T* B, int64_t n, int64_t s, int64_t d) { + using namespace std::chrono; + steady_clock::time_point _mv_t0, _mv_t1; + _t_matvec_us = 0; + + util::resize(K_big, K_big_sz, (d + 1) * n * s); + util::resize(R0_buf, R0_sz, s * s); + util::resize(tau_buf, tau_buf_sz, n); + util::resize(A_blk, A_blk_sz, d * s * s); + util::resize(B_blk, B_blk_sz, d * s * s); + if (reorth) util::resize(proj_buf, proj_buf_sz, s * s); + + // Initial QR: B = Q0 * R0. Q0 overwrites K_big[0..n*s-1]. + T* Q0 = K_big; + lapack::lacpy(lapack::MatrixType::General, n, s, B, n, Q0, n); + lapack::geqrf(n, s, Q0, n, tau_buf); + lapack::laset(lapack::MatrixType::General, s, s, (T)0, (T)0, R0_buf, s); + lapack::lacpy(lapack::MatrixType::Upper, s, s, Q0, n, R0_buf, s); + lapack::orgqr(n, s, s, Q0, n, tau_buf); + + for (int64_t step = 0; step < d; ++step) { + T* Q_step = K_big + step * n * s; + T* Q_prev = (step > 0) ? K_big + (step - 1) * n * s : nullptr; + T* Y = K_big + (step + 1) * n * s; // matvec output, then Z in-place + T* A_step = A_blk + step * s * s; + T* B_prev = (step > 0) ? B_blk + (step - 1) * s * s : nullptr; + + // Y = A * Q_step + if (this->timing) _mv_t0 = steady_clock::now(); + A(Layout::ColMajor, s, (T)1.0, Q_step, n, (T)0.0, Y, n); + if (this->timing) { _mv_t1 = steady_clock::now(); _t_matvec_us += duration_cast(_mv_t1 - _mv_t0).count(); } + + // Y -= Q_{step-1} * B_{step-1}^T + if (step > 0) + blas::gemm(Layout::ColMajor, Op::NoTrans, Op::Trans, + n, s, s, (T)-1.0, Q_prev, n, B_prev, s, (T)1.0, Y, n); + + // A_step = Q_step^T * Y (block alpha, s×s). Symmetric in + // exact arithmetic since Q_stepᵀ·A·Q_step is — symmetrize + // away the small finite-arithmetic asymmetry. + blas::gemm(Layout::ColMajor, Op::Trans, Op::NoTrans, + s, s, n, (T)1.0, Q_step, n, Y, n, (T)0.0, A_step, s); + util::symmetrize(s, A_step, s); + + // Y -= Q_step * A_step (Z = Y - Q_step*A_step, in-place) + blas::gemm(Layout::ColMajor, Op::NoTrans, Op::NoTrans, + n, s, s, (T)-1.0, Q_step, n, A_step, s, (T)1.0, Y, n); + + // Optional full reorthogonalization: Z -= Q_p * (Q_p^T * Z) for each prev block + if (reorth) { + for (int64_t prev = 0; prev <= step; ++prev) { + T* Q_p = K_big + prev * n * s; + blas::gemm(Layout::ColMajor, Op::Trans, Op::NoTrans, + s, s, n, (T)1.0, Q_p, n, Y, n, (T)0.0, proj_buf, s); + blas::gemm(Layout::ColMajor, Op::NoTrans, Op::NoTrans, + n, s, s, (T)-1.0, Q_p, n, proj_buf, s, (T)1.0, Y, n); + } + } + + // QR of Z → Q_{step+1} (in Y) and B_step (upper R factor). + // Skipped at the last step: Q_d is never needed by apply_f. + if (step < d - 1) { + T* B_step = B_blk + step * s * s; + lapack::geqrf(n, s, Y, n, tau_buf); + lapack::laset(lapack::MatrixType::General, s, s, (T)0, (T)0, B_step, s); + lapack::lacpy(lapack::MatrixType::Upper, s, s, Y, n, B_step, s); + lapack::orgqr(n, s, s, Y, n, tau_buf); + } + } + } + + // ------------------------------------------------------------------ + /// Evaluate f(A)B from precomputed Krylov data (K_big, R0_buf, A_blk, B_blk). + /// + /// Computation: + /// 1. Assemble T_dense (d*s × d*s block tridiagonal) from A_blk and B_blk. + /// 2. syevd: T_dense → eigenvectors V (in-place), eigenvalues λ. + /// 3. W (s×m): W[i,j] = f(λⱼ)*V[i,j] for i=0..s-1 (first s rows of V, col-scaled). + /// 4. C1 (d*s × s) = V * W^T — this equals f(T_k) * E₁. + /// 5. C1 *= R₀ (TRMM: right-multiply by upper-triangular R₀). + /// 6. out (n × s) = Q_basis * C1 (GEMM). + template F> + void apply_f(F f, int64_t n, int64_t s, int64_t d, T* out) { + int64_t m = d * s; + util::resize(workspace, workspace_sz, m * m + m + 2 * m * s); + + T* T_dense = workspace; + T* eig_vals = T_dense + m * m; + T* G = eig_vals + m; + T* C1 = G + m * s; + + // 1. Assemble T_dense (m×m, m = d·s) as block tridiagonal: + // + // ┌ A₀ B₀ᵀ ┐ + // │ B₀ A₁ B₁ᵀ │ + // T_d = │ B₁ A₂ … │ + // │ ⋱ ⋱ B_{d-2}ᵀ│ + // └ B_{d-2} A_{d-1}┘ + // + // where Aᵢ is the s×s diagonal block in A_blk[i] and Bᵢ is the s×s + // upper-triangular factor produced by the QR of Z at step i+1. + std::memset(T_dense, 0, m * m * sizeof(T)); + for (int64_t step = 0; step < d; ++step) { + int64_t b0 = step * s; + // Diagonal block + lapack::lacpy(lapack::MatrixType::General, s, s, + A_blk + step * s * s, s, T_dense + b0 * m + b0, m); + // Off-diagonal blocks (B_step in lower, B_step^T in upper) + if (step < d - 1) { + T* B_step = B_blk + step * s * s; + int64_t b1 = (step + 1) * s; + // Lower off-diagonal: T(b1:b1+s, b0:b0+s) = B_step + lapack::lacpy(lapack::MatrixType::General, s, s, + B_step, s, T_dense + b0 * m + b1, m); + // Upper off-diagonal: T(b0:b0+s, b1:b1+s) = B_step^T + for (int64_t j = 0; j < s; ++j) + for (int64_t i = 0; i < s; ++i) + T_dense[(b1 + j) * m + (b0 + i)] = B_step[i * s + j]; + } + } + // Symmetrize to eliminate any floating-point asymmetry between the + // explicit B_step lower block and its hand-transposed upper twin. + util::symmetrize(m, T_dense, m); + + // 2. Eigendecomposition: T_dense → V (eigenvectors overwrite T_dense), eig_vals → λ. + lapack::syevd(lapack::Job::Vec, blas::Uplo::Lower, m, T_dense, m, eig_vals); + + // 3. W (s×m col-major, ld=s; stored in the G buffer): W[i,j] = f(λⱼ)*V[i,j] + // for i=0..s-1, j=0..m-1. Each column j of W is the first s elements of + // column j of V, scaled by f(λⱼ). + // + // Derivation: f(T_k)*E₁ = V*diag(f(λ))*V^T*E₁. + // V^T*E₁ = first s columns of V^T = {row i of V} for i=0..s-1 stacked as cols. + // diag(f(λ))*(V^T*E₁): scale row j → f(λⱼ)*(row j of V^T*E₁) = f(λⱼ)*V[:,j][0:s]. + // That intermediate is W^T (m×s), so f(T_k)*E₁ = V * W^T. + T* W = G; // reuse G buffer; W is s×m col-major (ld=s) + for (int64_t j = 0; j < m; ++j) { + T fev = f(std::max(eig_vals[j], (T)0)); + const T* V_col = T_dense + j * m; // col j of V (contiguous, length m) + T* W_col = W + j * s; // col j of W (contiguous, length s) + for (int64_t i = 0; i < s; ++i) + W_col[i] = fev * V_col[i]; // first s rows of V[:,j], scaled + } + + // 4. C1 (m × s) = V * W^T — equals f(T_k) * E₁. + // GEMM(NoTrans, Trans, m, s, m): C = V(m×m) * W^T where W is s×m (ld=s). + blas::gemm(Layout::ColMajor, Op::NoTrans, Op::Trans, + m, s, m, (T)1.0, T_dense, m, W, s, (T)0.0, C1, m); + + // 5. C1 *= R₀ (TRMM: C1 = C1 * R₀, right upper triangular). + blas::trmm(Layout::ColMajor, Side::Right, Uplo::Upper, Op::NoTrans, Diag::NonUnit, + m, s, (T)1.0, R0_buf, s, C1, m); + + // 6. out (n × s) = Q_basis (n × m) * C1 (m × s). + // Q_basis = K_big[0..m*n-1] (n×m col-major, ld=n). + blas::gemm(Layout::ColMajor, Op::NoTrans, Op::NoTrans, + n, s, m, (T)1.0, K_big, n, C1, m, (T)0.0, out, n); + } + + // ------------------------------------------------------------------ + /// Combined run + apply. + /// + /// Drop-in replacement for LanczosFA::call — same signature, slots into + /// FunNystromPP and ResidualOp as the LanczosFA_t template parameter. + template F> + void call(SLO& A, const T* B, int64_t n, int64_t s, F f, int64_t d, T* out) { + using namespace std::chrono; + _t_matvec_us = 0; + steady_clock::time_point t_total_start, t_lanczos_end, t_end; + if (this->timing) t_total_start = steady_clock::now(); + run_lanczos(A, B, n, s, d); + if (this->timing) t_lanczos_end = steady_clock::now(); + apply_f(f, n, s, d, out); + if (this->timing) { + t_end = steady_clock::now(); + long total_us = duration_cast(t_end - t_total_start).count(); + long lanczos_us = duration_cast(t_lanczos_end - t_total_start).count(); + long apply_f_us = duration_cast(t_end - t_lanczos_end).count(); + long rest_us = total_us - lanczos_us - apply_f_us; + this->times = {_t_matvec_us, lanczos_us, apply_f_us, rest_us, total_us}; + } + } +}; + + +} // end namespace RandLAPACK diff --git a/RandLAPACK/comps/rl_matfun.hh b/RandLAPACK/comps/rl_matfun.hh new file mode 100644 index 000000000..e6b142c74 --- /dev/null +++ b/RandLAPACK/comps/rl_matfun.hh @@ -0,0 +1,36 @@ +#pragma once + +#include +#include + +namespace RandLAPACK { + +/// SqrtFun: f(x) = sqrt(max(x, 0)). +/// Clamps to zero to tolerate floating-point roundoff near zero eigenvalues. +template +struct SqrtFun { + T operator()(T x) const { return std::sqrt(std::max(x, (T)0)); } +}; + +/// LogFun: f(x) = log(x). +/// Undefined for x <= 0. The caller must ensure the matrix is strictly PD +/// (all eigenvalues > 0) before using this function. +template +struct LogFun { + T operator()(T x) const { return std::log(x); } +}; + +/// PolyFun: f(x) = x*(x + lambda) = x^2 + lambda*x. +/// +/// Models tr(K(K + λI)) = tr(K² + λK) for a symmetric PSD kernel matrix K. +/// LanczosFA is exact for this function in d=2 Lanczos steps, since Gauss +/// quadrature integrates polynomials of degree ≤ 2d-1 exactly (d=2 covers +/// degree 3, which contains all degree-2 polynomials). +template +struct PolyFun { + T lam; + explicit PolyFun(T lambda) : lam(lambda) {} + T operator()(T x) const { return x * (x + lam); } +}; + +} // end namespace RandLAPACK diff --git a/RandLAPACK/comps/rl_preconditioners.hh b/RandLAPACK/comps/rl_preconditioners.hh index afbdcd090..58c08badd 100644 --- a/RandLAPACK/comps/rl_preconditioners.hh +++ b/RandLAPACK/comps/rl_preconditioners.hh @@ -9,7 +9,7 @@ #include "rl_syps.hh" #include "rl_syrf.hh" #include "rl_rpchol.hh" -#include "rl_revd2.hh" +#include "rl_nystrom_evd.hh" #include #include @@ -305,7 +305,7 @@ RandBLAS::RNGState nystrom_pc_data( // ^ Define the symmetric rangefinder algorithm. // (*) Use power sketching followed by Householder orthogonalization. // (*) Do not check condition numbers or log to std::out. - RandLAPACK::REVD2 NystromAlg(syrf, num_steps_power_iter_error_est, false); + RandLAPACK::NystromEVD NystromAlg(syrf, num_steps_power_iter_error_est, false); // ^ Define the algorithm for low-rank approximation via Nystrom. // (*) Handle accuracy requests by estimating ||A - V diag(eigvals) V'|| // with "num_steps_power_iter_error_est" steps of power iteration. @@ -313,7 +313,15 @@ RandBLAS::RNGState nystrom_pc_data( T tol = mu_min / 5; // ^ Set tolerance to something materially smaller than the smallest // regularization parameter the user claims to need. - return NystromAlg.call(A, k, tol, V, eigvals, state); + T* V_buf = nullptr; int64_t V_sz = 0; + T* ev_buf = nullptr; int64_t ev_sz = 0; + NystromAlg.call(A, k, tol, V_buf, V_sz, ev_buf, ev_sz, state); + int64_t m = A.dim; + V.assign(V_buf, V_buf + m * k); + eigvals.assign(ev_buf, ev_buf + k); + delete[] V_buf; + delete[] ev_buf; + return state; } /** diff --git a/RandLAPACK/comps/rl_syps.hh b/RandLAPACK/comps/rl_syps.hh index 4bda9a60c..29e244136 100644 --- a/RandLAPACK/comps/rl_syps.hh +++ b/RandLAPACK/comps/rl_syps.hh @@ -27,6 +27,16 @@ class SYPS { bool verbose; bool cond_check; std::vector cond_nums; + // 0 = Gaussian (DenseDist, default); 1 = SASO (SparseDist, vec_nnz per column) + int sketch_type = 0; + int vec_nnz = 4; + // Internal stabilization (applied between power iterations). + // nullptr → use the historical hardcoded geqrf+ungqr path (HQRQ-equivalent). + // Non-null → invoke the supplied Stabilization (e.g. CholQRQ) and fall back to + // the hardcoded geqrf+ungqr if it returns nonzero (e.g. CholQR Cholesky failure). + // Power iteration produces increasingly column-correlated matrices; CholQR can + // hit potrf failure on ill-conditioned inputs, so the fallback is intentional. + Stabilization* internal_stab = nullptr; SYPS( int64_t p, // number of passes @@ -109,8 +119,6 @@ class SYPS { bool callers_skop_buff = skop_buff != nullptr; if (!callers_skop_buff) skop_buff = new T[m * k]; - RandBLAS::DenseDist D(m, k); - state = RandBLAS::fill_dense(D, skop_buff, state); bool callers_work_buff = work_buff != nullptr; if (!callers_work_buff) @@ -120,18 +128,64 @@ class SYPS { T *symm_out = work_buff; T *symm_in = skop_buff; T *tau = new T[k]{}; - while (p - p_done > 0) { - A(Layout::ColMajor, k, 1.0, symm_in, m, 0.0, symm_out, m); - ++p_done; - if (p_done % q == 0) { - if(lapack::geqrf(m, k, symm_out, m, tau)) { - delete [] tau; + + // Inline lambda for the per-iter stabilization (deduplicates the path used + // for both the sparse-first-iter shortcut and the regular dense loop). + auto stabilize = [&](T* buf) { + bool fallback = (this->internal_stab == nullptr); + if (!fallback) { + fallback = (this->internal_stab->call(m, k, buf) != 0); + } + if (fallback) { + if (lapack::geqrf(m, k, buf, m, tau)) { + delete[] tau; throw std::runtime_error("GEQRF failed."); } - lapack::ungqr(m, k, k, symm_out, m, tau); + lapack::ungqr(m, k, k, buf, m, tau); + } + }; + + if (sketch_type == 1) { + // Sparse sketch path — apply A to the sparse Ω directly via the + // SkOp overload of A (dense·sparse spmm), skipping the n×k + // densification of Ω. T and sint_t (int64_t) must be explicit + // since SparseDist::sample can't deduce T from the RNGState. + auto S = RandBLAS::SparseDist(m, k, this->vec_nnz).sample(state); + state = S.next_state; + RandBLAS::fill_sparse(S); + + // SLOs that implement the SkOp-taking operator() get the fast + // sparse-aware first matvec. Other SLOs (RegExplicitSymLinOp, + // SparseSymLinOp without the SkOp overload, etc.) fall back + // to the legacy densification. + if constexpr (requires (SLO& a, decltype(S)& s) { + a(Layout::ColMajor, k, (T)1.0, s, (T)0.0, symm_out, m); + }) { + A(Layout::ColMajor, k, (T)1.0, S, (T)0.0, symm_out, m); + } else { + auto Scoo = RandBLAS::coo_view_of_skop(S); + std::fill(skop_buff, skop_buff + m * k, (T)0.0); + RandLAPACK::util::sparse_to_dense(Scoo, Layout::ColMajor, skop_buff); + A(Layout::ColMajor, k, (T)1.0, skop_buff, m, (T)0.0, symm_out, m); } + ++p_done; + if (p_done % q == 0) stabilize(symm_out); + symm_out = (p_done % 2 == 1) ? skop_buff : work_buff; + symm_in = (p_done % 2 == 1) ? work_buff : skop_buff; + } else { + // Dense Gaussian sketch — fill skop_buff in place; first matvec + // happens uniformly in the loop below. + RandBLAS::DenseDist D(m, k); + state = RandBLAS::fill_dense(D, skop_buff, state); + } + + // Remaining iterations are pure dense matvecs through the SLO interface. + while (p - p_done > 0) { + A(Layout::ColMajor, k, 1.0, symm_in, m, 0.0, symm_out, m); + ++p_done; + if (p_done % q == 0) stabilize(symm_out); symm_out = (p_done % 2 == 1) ? skop_buff : work_buff; - symm_in = (p_done % 2 == 1) ? work_buff : skop_buff; + symm_in = (p_done % 2 == 1) ? work_buff : skop_buff; } delete[] tau; if (p % 2 == 1) diff --git a/RandLAPACK/comps/rl_syrf.hh b/RandLAPACK/comps/rl_syrf.hh index 66d4da01c..7b225a743 100644 --- a/RandLAPACK/comps/rl_syrf.hh +++ b/RandLAPACK/comps/rl_syrf.hh @@ -74,7 +74,7 @@ class SYRF { int64_t m, const T* A, int64_t k, - std::vector &Q, + T* Q, RandBLAS::RNGState &state, T* work_buff ) { @@ -86,7 +86,7 @@ class SYRF { int call( SLO &A, int64_t k, - std::vector &Q, + T* Q, RandBLAS::RNGState &state, T* work_buff ) { @@ -97,19 +97,18 @@ class SYRF { RandBLAS::util::safe_scal(m * k, (T) 0.0, work_buff); - T* Q_dat = util::upsize(m * k, Q); - syps.call(A, k, state, work_buff, Q_dat); + syps.call(A, k, state, work_buff, Q); // Q = orth(A * Omega) - A(Layout::ColMajor, k, (T) 1.0, work_buff, m, (T) 0.0, Q_dat, m); + A(Layout::ColMajor, k, (T) 1.0, work_buff, m, (T) 0.0, Q, m); if(this->cond_check) { - util::upsize(m * k, this->cond_work_mat); - util::upsize(k, this->cond_work_vec); + util::resize(m * k, this->cond_work_mat); + util::resize(k, this->cond_work_vec); this->cond_nums.push_back( - util::cond_num_check(m, k, Q.data(), this->verbose) + util::cond_num_check(m, k, Q, this->verbose) ); } - if(this->orth.call(m, k, Q.data())) + if(this->orth.call(m, k, Q)) throw std::runtime_error("Orthogonalization failed."); if (!callers_work_buff) diff --git a/RandLAPACK/drivers/rl_fun_nystrom_pp.hh b/RandLAPACK/drivers/rl_fun_nystrom_pp.hh new file mode 100644 index 000000000..a13434957 --- /dev/null +++ b/RandLAPACK/drivers/rl_fun_nystrom_pp.hh @@ -0,0 +1,223 @@ +#pragma once + +#include "rl_blaspp.hh" +#include "rl_linops.hh" +#include "rl_util.hh" + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace RandLAPACK { + +/// Per-call timing breakdown for FunNystromPP::call. +/// Pass a pointer to call() to collect Phase 1 and Phase 2 wall-clock times. +struct FunNystromPP_timing { + long phase1_us = 0; // NystromEVD Nyström approximation (Phase 1) + long phase2_us = 0; // Hutchinson trace correction (Phase 2) +}; + +/// funNyström++ trace estimator: estimates tr(f(A)) for symmetric PSD A. +/// +/// Algorithm (2 phases): +/// Phase 1 — Nyström approximation via NystromEVD (O(k) matvecs): +/// (V, λ) = NystromEVD(A, k) →  = V diag(λ) V^T +/// tr(f(Â)) = Σ f(λ_i) + (n-k) * f(0) [free from the eigenvalues] +/// +/// Phase 2 — Hutchinson correction on the residual (d*s matvecs): +/// Ω = Rademacher(n×s) +/// Z1 = f(A)*Ω via LanczosFA [d batch matvecs] +/// Z2 = f(Â)*Ω = V * (f(λ) ⊙ (V^T Ω)) [two GEMMs, free] +/// correction = <Ω, Z1 - Z2>_F / s +/// +/// return tr(f(Â)) + correction +/// +/// Complexity: O(√κ/ε) matvecs total, versus O(√κ/ε²) for naive Hutchinson. +/// The variance reduction comes from subtracting the control variate f(Â)Ω: +/// because  ≈ A in the dominant eigenspace, f(A)-f(Â) has small Frobenius +/// norm and the Hutchinson estimator converges in far fewer samples. +/// +/// Requirements on A: +/// - Symmetric positive semidefinite (λ_min ≥ 0) for most functions +/// (e.g., sqrt, x^α with α > 0 where f(0) is finite). +/// - For f = log, A must be strictly positive definite (λ_min > 0), +/// since log(0) = -∞. Pass any finite f_zero placeholder and ensure +/// λ_min > 0 at the call site. +/// +/// Parameter selection (f = sqrt, κ = λ_max/λ_min): +/// d = O(√κ log(n/ε)) Lanczos steps — caller's responsibility +/// s = O(1/(√κ ε)) Hutchinson samples — caller's responsibility +/// k = O(√κ/ε) Nyström rank — caller's responsibility +/// Use error_est_power_iters=0 in NystromEVD and tol=0 to get fixed-rank behavior +/// (k never changes; the adaptive loop exits after one iteration). +/// +/// @tparam NystromEVD_t Type of the NystromEVD Nyström component. +/// @tparam LanczosFA_t Type of the LanczosFA component. +/// @tparam Hutchinson_t Type of the Hutchinson trace estimator component. +template +class FunNystromPP { +public: + using T = typename NystromEVD_t::T; + using RNG = typename NystromEVD_t::RNG; + + NystromEVD_t& nystrom_evd; // Phase 1: Nyström approximation + LanczosFA_t& lanczos_fa; // Phase 2: f(A)*Ω oracle + Hutchinson_t& hutchinson; // Phase 2: trace correction estimator + + // Persistent output/working buffers — grown with new/delete[], never shrunk. + // V, eigvals: NystromEVD outputs reused across calls with same (n, k). + // F_vec: f applied elementwise to eigvals. + // tmp, Z1, Z2: workspace owned here and lent to linops::ResidualOp each call. + T* V_buf = nullptr; int64_t V_buf_sz = 0; + T* eigvals_buf = nullptr; int64_t eigvals_buf_sz = 0; + T* F_vec = nullptr; int64_t F_vec_sz = 0; + T* tmp = nullptr; int64_t tmp_sz = 0; + T* Z1 = nullptr; int64_t Z1_sz = 0; + T* Z2 = nullptr; int64_t Z2_sz = 0; + + FunNystromPP(const FunNystromPP&) = delete; + FunNystromPP& operator=(const FunNystromPP&) = delete; + + ~FunNystromPP() { + delete[] V_buf; delete[] eigvals_buf; + delete[] F_vec; delete[] tmp; delete[] Z1; delete[] Z2; + } + + // ------------------------------------------------------------------ + /// Construct a FunNystromPP estimator from pre-built components. + /// + /// @param r NystromEVD instance for Phase 1 (Nyström approximation). + /// @param l LanczosFA instance for computing f(A)*Ω in Phase 2. + /// @param h Hutchinson instance for Phase 2 trace correction. + /// + /// All three components are stored by reference and must outlive this object. + FunNystromPP(NystromEVD_t& r, LanczosFA_t& l, Hutchinson_t& h) + : nystrom_evd(r), lanczos_fa(l), hutchinson(h) {} + + // ------------------------------------------------------------------ + /// Estimate tr(f(A)). + /// + /// @param[in] A SymmetricLinearOperator — matvec oracle for A. + /// @param[in] f Scalar function T→T (e.g., sqrt, log, pow(x,α)). + /// @param[in] f_zero Optional f(0). When std::nullopt (the default), and + /// k < n, the value is auto-derived as f((T)0); for + /// functions like log where f(0) is non-finite, the + /// caller must pass an explicit value or guarantee + /// k == n. Must be finite if provided. When k == n + /// (Phase 1 captures the full spectrum), f_zero is + /// unused and may be omitted. + /// @param[in] k Nyström rank. Caller chooses based on κ, ε. + /// @param[in] s Hutchinson samples. Caller chooses based on κ, ε. + /// @param[in] d Lanczos steps per sample. Caller chooses based on κ. + /// @param[in] state RandBLAS RNG state; advanced on return. + /// @param[out] timing Optional. If non-null, receives Phase 1 and Phase 2 + /// wall-clock times in microseconds. + /// @returns Estimate of tr(f(A)). + template F> + T call( + SLO& A, + F f, + std::optional f_zero, + int64_t k, + int64_t s, + int64_t d, + RandBLAS::RNGState& state, + FunNystromPP_timing* timing = nullptr + ) { + if (f_zero.has_value() && !std::isfinite(*f_zero)) + throw std::invalid_argument("f_zero must be finite when provided"); + + int64_t n = A.dim; + + // ------------------------------------------------------------------ + // Phase 1: Nyström approximation + // NystromEVD with error_est_power_iters=0 and tol=0 runs a single pass + // with fixed rank k. V_buf (n×k_in) and eigvals_buf (k_in) are the outputs. + // ------------------------------------------------------------------ + // NystromEVD::call takes rank by reference and may increase it adaptively. + // We pass k_in (a copy) so the caller's k stays unchanged. + // All post-call uses reference k_in (not k) so they stay consistent + // with the actual rank used, even in adaptive mode. + int64_t k_in = k; + auto t_phase1_start = std::chrono::steady_clock::now(); + nystrom_evd.call(A, k_in, (T)0.0, V_buf, V_buf_sz, eigvals_buf, eigvals_buf_sz, state); + auto t_phase1_end = std::chrono::steady_clock::now(); + + // f(λ): apply scalar f to k_in eigenvalues + util::resize(F_vec, F_vec_sz, k_in); + std::transform(eigvals_buf, eigvals_buf + k_in, F_vec, f); + + // Resolve f_zero. The (n-k_in)*f(0) term is only used when k_in < n; + // when k_in == n the full spectrum is captured by Phase 1 and f_zero is + // unused (so we don't require the caller to provide one). + T fz; + if (f_zero.has_value()) { + fz = *f_zero; + } else if (k_in < n) { + fz = f((T)0); + if (!std::isfinite(fz)) + throw std::invalid_argument( + "f(0) is non-finite; pass explicit f_zero or ensure k captures full rank"); + } else { + fz = (T)0; // unused: (n - k_in) == 0 + } + + // tr(f(Â)) = Σ f(λ_i) + (n-k_in)*f(0) + T tr_Ahat = (T)0.0; + for (int64_t i = 0; i < k_in; ++i) + tr_Ahat += F_vec[i]; + tr_Ahat += static_cast(n - k_in) * fz; + + // ------------------------------------------------------------------ + // Phase 2: Hutchinson correction on residual f(A) - f(Â) + // linops::ResidualOp wraps the correction as a SymmetricLinearOperator + // so Hutchinson::call can draw Ω and compute <Ω, (f(A)-f(Â))Ω>_F / s. + // + // Skip Phase 2 when k_in == n: NystromEVD has captured the full spectrum + // exactly ( = A), so f(A) - f(Â) is analytically zero. Running Phase 2 + // anyway would compute <Ω, (Z1 - Z2)>/s where Z1 (LFA approximation) and + // Z2 (V·diag(f(λ))·Vᵀ·Ω) follow different floating-point paths, leaving + // a misleading O(ε_mach × s × LFA_residual) "noise floor" on the trace. + // The structural bound k + s ≤ n implies that when k = n, s must be 0 + // for a sensible matvec budget — this skip is its natural consequence. + // ------------------------------------------------------------------ + T correction = (T)0; + auto t_phase2_start = std::chrono::steady_clock::now(); + auto t_phase2_end = t_phase2_start; + + if (k_in < n) { + util::resize(tmp, tmp_sz, k_in * s); + util::resize(Z1, Z1_sz, n * s); + util::resize(Z2, Z2_sz, n * s); + + linops::ResidualOp res_op( + n, A, lanczos_fa, f, d, k_in, V_buf, F_vec, tmp, Z1, Z2 + ); + t_phase2_start = std::chrono::steady_clock::now(); + correction = hutchinson.call(res_op, s, state); + t_phase2_end = std::chrono::steady_clock::now(); + } else { + // Phase 2 skipped: lanczos_fa.call() never executed, so its + // .times vector stays empty. Populate with zeros (5 slots: matvec, + // lanczos, apply_f, rest, total) so callers indexing into it for + // timing breakdowns get well-defined zeros instead of UB. + lanczos_fa.times.assign(5, 0L); + } + + if (timing) { + using us = std::chrono::microseconds; + timing->phase1_us = std::chrono::duration_cast(t_phase1_end - t_phase1_start).count(); + timing->phase2_us = std::chrono::duration_cast(t_phase2_end - t_phase2_start).count(); + } + + return tr_Ahat + correction; + } +}; + + +} // end namespace RandLAPACK diff --git a/RandLAPACK/drivers/rl_nystrom_evd.hh b/RandLAPACK/drivers/rl_nystrom_evd.hh new file mode 100644 index 000000000..5389fa11e --- /dev/null +++ b/RandLAPACK/drivers/rl_nystrom_evd.hh @@ -0,0 +1,363 @@ +#pragma once + +#include "rl_syps.hh" +#include "rl_syrf.hh" +#include "rl_blaspp.hh" +#include "rl_lapackpp.hh" +#include "rl_util.hh" +#include "rl_linops.hh" + +#include +#include +#include +#include + +namespace RandLAPACK { + + +// ----------------------------------------------------------------------------- +/// Power scheme for error estimation, based on Algorithm E.1 from https://arxiv.org/pdf/2110.02820.pdf. +/// p - number of algorithm iterations +/// vector_buf - buffer for vector operations (must hold at least 4*m elements) +/// Mat_buf - buffer for matrix operations (must hold at least m*k elements) +/// All other parameters come from NystromEVD +template +T power_error_est( + SLO &A, + int64_t k, + int p, + T* vector_buf, + T* V, + T* Mat_buf, + T* eigvals +) { + int64_t m = A.dim; + T err = 0; + for(int iter = 0; iter < p; ++iter) { + T g_norm = blas::nrm2(m, vector_buf, 1); + blas::scal(m, (T)1 / g_norm, vector_buf, 1); + + // V' * g/||g|| — result in column 1 of vector_buf + gemv(Layout::ColMajor, Op::Trans, m, k, (T)1, V, m, vector_buf, 1, (T)0, &vector_buf[m], 1); + + // V * diag(eigvals) into Mat_buf — column-major: column col occupies rows [0, m) + for (int64_t col = 0; col < k; ++col) + for (int64_t row = 0; row < m; ++row) + Mat_buf[col * m + row] = V[col * m + row] * eigvals[col]; + + // V * diag(eigvals) * V' * g/||g|| — result in column 2 of vector_buf + gemv(Layout::ColMajor, Op::NoTrans, m, k, (T)1, Mat_buf, m, &vector_buf[m], 1, (T)0, &vector_buf[2 * m], 1); + // A * g/||g|| — result in column 3 of vector_buf + A(Layout::ColMajor, 1, (T)1, vector_buf, m, (T)0, &vector_buf[3*m], m); + + // w = A*g/||g|| - V*diag(eigvals)*V'*g/||g|| + blas::axpy(m, (T)-1, &vector_buf[2 * m], 1, &vector_buf[3 * m], 1); + // error estimate: (g/||g||)' * w + err = blas::dot(m, vector_buf, 1, &vector_buf[3 * m], 1); + std::copy(&vector_buf[3 * m], &vector_buf[4 * m], vector_buf); + } + return err; +} + + +template +class NystromEVD { + public: + using T = typename SYRF_t::T; + using RNG = typename SYRF_t::RNG; + SYRF_t &syrf; + int error_est_p; + bool verbose; + + // Internal working buffers — owned by this object, grown with new[]/delete[] via util::resize + T* Y = nullptr; int64_t Y_sz = 0; + T* Omega = nullptr; int64_t Omega_sz = 0; + T* R = nullptr; int64_t R_sz = 0; + T* S = nullptr; int64_t S_sz = 0; + T* symrf_work = nullptr; int64_t symrf_work_sz = 0; + + bool timing = false; + std::vector times; // populated after call() when timing==true + // Slots: alloc, syrf, matvec, gram, potrf, trsm, svd, post_svd, error_est, rest, total + + NystromEVD( + SYRF_t &syrf_obj, + int error_est_power_iters, + bool verb = false + ) : syrf(syrf_obj) { + error_est_p = error_est_power_iters; + verbose = verb; + } + + NystromEVD(const NystromEVD&) = delete; + NystromEVD& operator=(const NystromEVD&) = delete; + + ~NystromEVD() { + delete[] Y; + delete[] Omega; + delete[] R; + delete[] S; + delete[] symrf_work; + } + + /// Computes a rank-k approximation to an EVD of a symmetric positive semidefinite matrix: + /// A_hat = V diag(eigvals) V^*, + /// where V is a matrix of approximate eigenvectors and eigvals holds the eigenvalues. + /// + /// Adaptive: if tolerance is not met, doubles k until convergence or k == m. + /// Set error_est_power_iters=0 and tol=0 for fixed-rank single-pass behavior. + /// + /// This code is identical to Algorithm E2 from https://arxiv.org/pdf/2110.02820.pdf. + /// + /// @param[in] uplo Triangle of A that holds the data (Upper or Lower). + /// @param[in] m Number of rows/cols of A. + /// @param[in] A m×m matrix, column-major. + /// @param[in,out] k Sketch rank on entry; final rank on exit (may grow adaptively). + /// @param[in] tol Convergence tolerance. Use 0.0 for fixed-rank. + /// @param[in,out] V Caller-owned buffer; on exit holds m×k eigenvectors. + /// @param[in,out] V_sz Capacity of V in elements; grown by util::resize as needed. + /// @param[in,out] eigvals Caller-owned buffer; on exit holds k eigenvalues. + /// @param[in,out] eigvals_sz Capacity of eigvals; grown by util::resize as needed. + /// + int call( + Uplo uplo, + int64_t m, + const T* A, + int64_t &k, + T tol, + T*& V, int64_t& V_sz, + T*& eigvals, int64_t& eigvals_sz, + RandBLAS::RNGState &state + ) { + linops::ExplicitSymLinOp A_linop(m, uplo, A, m, Layout::ColMajor); + return this->call(A_linop, k, tol, V, V_sz, eigvals, eigvals_sz, state); + } + + template + int call( + SLO &A, + int64_t &k, + T tol, + T*& V, int64_t& V_sz, + T*& eigvals, int64_t& eigvals_sz, + RandBLAS::RNGState &state + ) { + using namespace std::chrono; + int64_t m = A.dim; + T err = 0; + RandBLAS::RNGState error_est_state(state.counter, state.key); + error_est_state.key.incr(1); + + steady_clock::time_point t0, t1, total_t_start, total_t_stop; + long alloc_t_dur = 0; + long syrf_t_dur = 0; + long matvec_t_dur = 0; + long gram_t_dur = 0; + long potrf_t_dur = 0; + long trsm_t_dur = 0; + long svd_t_dur = 0; + long post_svd_t_dur = 0; + long error_est_t_dur = 0; + long total_t_dur = 0; + + if (this->timing) total_t_start = steady_clock::now(); + + while(true) { + if (this->timing) t0 = steady_clock::now(); + util::resize(eigvals, eigvals_sz, k); + util::resize(V, V_sz, m * k); + T* V_dat = V; + util::resize(Y, Y_sz, m * k); + util::resize(Omega, Omega_sz, std::max(m * k, m * (int64_t)4)); + util::resize(R, R_sz, k * k); + util::resize(S, S_sz, k * k); + util::resize(symrf_work, symrf_work_sz, m * k); + if (this->timing) { t1 = steady_clock::now(); alloc_t_dur += duration_cast(t1 - t0).count(); } + + if (this->timing) t0 = steady_clock::now(); + this->syrf.call(A, k, this->Omega, state, symrf_work); + if (this->timing) { t1 = steady_clock::now(); syrf_t_dur += duration_cast(t1 - t0).count(); } + + if (this->timing) t0 = steady_clock::now(); + A(Layout::ColMajor, k, (T)1, Omega, m, (T)0, Y, m); + if (this->timing) { t1 = steady_clock::now(); matvec_t_dur += duration_cast(t1 - t0).count(); } + + // ----- Phase 1 spectral recovery ----- + // + // Two paths sharing the same Gram-matrix construction: + // + // FAST PATH (well-conditioned Gram — the common case): + // G = Ωᵀ · Y (k×k sym PSD) + // G = Lᵀ · L (Cholesky) + // Y := Y · L⁻¹ (TRSM) + // G_y = Yᵀ · Y (small Gram) + // G_y = V · Σ² · Vᵀ (syevd) + // U = Y · V · diag(1/σ) (gemm + scal) + // eigvals = Σ² + // + // The eig-of-Yᵀ·Y trick (Halko-Martinsson-Tropp §5.1) avoids + // the more expensive gesdd of n×k while keeping eigenvector + // accuracy fine because Y is well-conditioned post-TRSM. + // + // FALL-BACK (Cholesky failure on near-singular Gram): + // syevd(G) → V_G, D (pseudoinverse) + // B = Y · V_G · diag(1/√D) · V_Gᵀ (3 GEMMs) + // gesdd(B) → U, Σ + // eigvals = Σ² + // + // This is the safe path. It handles rank-deficient Gram via + // thresholded pseudoinverse. Mirrors the Persson-Kressner + // reference nystrom.m. Slightly more expensive than the + // fast path but always succeeds. + // + // No regularization shift in either path — both paths preserve + // the ε_mach accuracy of the SVD-pseudoinverse rewrite. + + T eps_mach = std::numeric_limits::epsilon(); + + // Form Gram G := Ωᵀ · Y into R (k×k). G = Ωᵀ·A·Ω in exact + // arithmetic and is symmetric, but the two non-trivial GEMMs + // (sketch + power-iter Y, then Ωᵀ Y) introduce small + // asymmetry; symmetrize so potrf and the later syevd see a + // numerically symmetric matrix. + if (this->timing) t0 = steady_clock::now(); + blas::gemm(Layout::ColMajor, Op::Trans, Op::NoTrans, k, k, m, + (T)1, Omega, m, Y, m, (T)0, R, k); + RandLAPACK::util::symmetrize(k, R, k); + // Backup G into symrf_work so we can restore it if Cholesky + // fails (potrf is destructive). symrf_work is m×k = enough for + // the k×k Gram. + lapack::lacpy(lapack::MatrixType::General, k, k, R, k, symrf_work, k); + if (this->timing) { t1 = steady_clock::now(); gram_t_dur += duration_cast(t1 - t0).count(); } + + int64_t r = 0; + + // Try Cholesky. + if (this->timing) t0 = steady_clock::now(); + int chol_status = lapack::potrf(Uplo::Upper, k, R, k); + if (this->timing) { t1 = steady_clock::now(); potrf_t_dur += duration_cast(t1 - t0).count(); } + + if (chol_status == 0) { + // -------- FAST PATH -------- + RandLAPACK::util::get_U(k, k, R, k); // zero strict lower + + // Y := Y · R⁻¹ (R is upper-tri Cholesky factor) + if (this->timing) t0 = steady_clock::now(); + blas::trsm(Layout::ColMajor, Side::Right, Uplo::Upper, + Op::NoTrans, Diag::NonUnit, m, k, + (T)1, R, k, Y, m); + if (this->timing) { t1 = steady_clock::now(); trsm_t_dur += duration_cast(t1 - t0).count(); } + + // G_y = Yᵀ · Y (k×k symmetric); reuse R buffer. syrk + // writes only the upper triangle; the strict lower stays + // unspecified, which is fine — the following syevd call + // is invoked with Uplo::Upper and reads only the upper + // triangle (then overwrites all of R with eigenvectors). + if (this->timing) t0 = steady_clock::now(); + blas::syrk(Layout::ColMajor, Uplo::Upper, Op::Trans, + k, m, (T)1, Y, m, (T)0, R, k); + // syevd(G_y) → R (eigvecs V), S (Σ², ASCENDING order). + lapack::syevd(lapack::Job::Vec, Uplo::Upper, k, R, k, S); + // Reverse to descending order to match the rest of the + // pipeline (post-SVD threshold logic expects S[0]=largest). + for (int64_t ii = 0; ii < k / 2; ++ii) std::swap(S[ii], S[k - 1 - ii]); + for (int64_t jj = 0; jj < k / 2; ++jj) { + blas::swap(k, R + jj * k, 1, R + (k - 1 - jj) * k, 1); + } + // U = Y · V (left singular vectors, unscaled). + blas::gemm(Layout::ColMajor, Op::NoTrans, Op::NoTrans, m, k, k, + (T)1, Y, m, R, k, (T)0, V_dat, m); + // Scale columns by 1/σ (= 1/sqrt(σ²)). Threshold tiny σ² as 0. + T sig_sq_max = (k > 0) ? std::max(S[0], (T)0) : (T)0; + T sig_sq_thresh = (T)4 * eps_mach * sig_sq_max; + for (int64_t jj = 0; jj < k; ++jj) { + if (S[jj] > sig_sq_thresh) { + T inv_sig = (T)1 / std::sqrt(S[jj]); + blas::scal(m, inv_sig, V_dat + jj * m, 1); + ++r; + } else { + std::fill(V_dat + jj * m, V_dat + (jj + 1) * m, (T)0); + } + } + // eigvals = σ² (already descending). + for (int64_t i = 0; i < k; ++i) + eigvals[i] = std::max(S[i], (T)0); + if (this->timing) { t1 = steady_clock::now(); svd_t_dur += duration_cast(t1 - t0).count(); } + } else { + // -------- FALL-BACK: syevd-pseudoinverse -------- + // Restore Gram into R from backup. + lapack::lacpy(lapack::MatrixType::General, k, k, symrf_work, k, R, k); + + if (this->timing) t0 = steady_clock::now(); + lapack::syevd(lapack::Job::Vec, Uplo::Lower, k, R, k, S); + if (this->timing) { t1 = steady_clock::now(); potrf_t_dur += duration_cast(t1 - t0).count(); } + + if (this->timing) t0 = steady_clock::now(); + T D_max = (k > 0) ? std::max(S[k - 1], (T)0) : (T)0; + T D_thresh = (T)2 * eps_mach * D_max; + // symrf_work was the Gram backup; safe to overwrite now. + blas::gemm(Layout::ColMajor, Op::NoTrans, Op::NoTrans, m, k, k, + (T)1, Y, m, R, k, (T)0, symrf_work, m); + for (int64_t jj = 0; jj < k; ++jj) { + T scale = (S[jj] > D_thresh) ? (T)1 / std::sqrt(S[jj]) : (T)0; + blas::scal(m, scale, symrf_work + jj * m, 1); + } + blas::gemm(Layout::ColMajor, Op::NoTrans, Op::Trans, m, k, k, + (T)1, symrf_work, m, R, k, (T)0, Y, m); + if (this->timing) { t1 = steady_clock::now(); trsm_t_dur += duration_cast(t1 - t0).count(); } + + if (this->timing) t0 = steady_clock::now(); + lapack::gesdd(Job::SomeVec, m, k, Y, m, S, V_dat, m, R, k); + if (this->timing) { t1 = steady_clock::now(); svd_t_dur += duration_cast(t1 - t0).count(); } + + T S_max = (k > 0) ? S[0] : (T)0; + T S_thresh = (T)2 * eps_mach * S_max; + for (int64_t i = 0; i < k; ++i) { + eigvals[i] = std::pow(S[i], 2); + if (S[i] > S_thresh) ++r; + } + std::fill(&V_dat[m * r], &V_dat[m * k], (T)0); + } + + if (this->timing) t0 = steady_clock::now(); + if (this->timing) { t1 = steady_clock::now(); post_svd_t_dur += duration_cast(t1 - t0).count(); } + + // Fixed-rank mode: skip error estimation entirely. + if (this->error_est_p == 0 && tol == (T)0) break; + + if (this->timing) t0 = steady_clock::now(); + RandBLAS::DenseDist g(m, 1); + error_est_state = RandBLAS::fill_dense(g, Omega, error_est_state); + err = power_error_est(A, k, this->error_est_p, Omega, V_dat, Y, eigvals); + if (this->timing) { t1 = steady_clock::now(); error_est_t_dur += duration_cast(t1 - t0).count(); } + + // Convergence: stop when the residual estimate hits the + // numerical noise floor (~ ε_mach · λ_max ≈ ε_mach · eigvals[0] + // since gesdd returned descending Σ → eigvals[0] is largest) + // or when the rank has saturated. + T noise_floor = eps_mach * (k > 0 ? eigvals[0] : (T)0); + if (err <= 5 * std::max(tol, noise_floor) || k == m) { + break; + } else if (2 * k > m) { + k = m; + } else { + k = 2 * k; + } + } + + if (this->timing) { + total_t_stop = steady_clock::now(); + total_t_dur = duration_cast(total_t_stop - total_t_start).count(); + long t_rest = total_t_dur - (alloc_t_dur + syrf_t_dur + matvec_t_dur + gram_t_dur + + potrf_t_dur + trsm_t_dur + svd_t_dur + post_svd_t_dur + error_est_t_dur); + this->times = {alloc_t_dur, syrf_t_dur, matvec_t_dur, gram_t_dur, + potrf_t_dur, trsm_t_dur, svd_t_dur, post_svd_t_dur, + error_est_t_dur, t_rest, total_t_dur}; + } + return 0; + } + +}; + + +} // end namespace RandLAPACK diff --git a/RandLAPACK/drivers/rl_revd2.hh b/RandLAPACK/drivers/rl_revd2.hh deleted file mode 100644 index 3667691a3..000000000 --- a/RandLAPACK/drivers/rl_revd2.hh +++ /dev/null @@ -1,241 +0,0 @@ -#pragma once - -#include "rl_syps.hh" -#include "rl_syrf.hh" -#include "rl_blaspp.hh" -#include "rl_lapackpp.hh" -#include "rl_util.hh" -#include "rl_linops.hh" - -#include -#include -#include - -namespace RandLAPACK { - - -// ----------------------------------------------------------------------------- -/// Power scheme for error estimation, based on Algorithm E.1 from https://arxiv.org/pdf/2110.02820.pdf. -/// This routine is too specialized to be included into RandLAPACK::utils -/// p - number of algorithm iterations -/// vector_buf - buffer for vector operations -/// Mat_buf - buffer for matrix operations -/// All other parameters come from REVD2 -template -T power_error_est( - SLO &A, - int64_t k, - int p, - T* vector_buf, - T* V, - T* Mat_buf, - T* eigvals -) { - int64_t m = A.dim; - T err = 0; - for(int i = 0; i < p; ++i) { - T g_norm = blas::nrm2(m, vector_buf, 1); - // Compute g = g / ||g|| - we need this because dot product does not take in an alpha - blas::scal(m, 1 / g_norm, vector_buf, 1); - - // Compute V' * g / ||g|| - // Using the second column of vector_buff as a buffer for matrix-vector product - gemv(Layout::ColMajor, Op::Trans, m, k, 1.0, V, m, vector_buf, 1, 0.0, &vector_buf[m], 1); - - // Compute V*E, eigvals diag - // Using Mat_buf as a buffer for V * diag(eigvals). - for (int i = 0, j = 0; i < m * k; ++i) { - Mat_buf[i] = V[i] * eigvals[j]; - if((i + 1) % m == 0 && i != 0) - ++j; - } - - // Compute V * diag(eigvals) * V' * g / ||g|| - // Using the third column of vector_buf as a buffer for matrix-vector product - gemv(Layout::ColMajor, Op::NoTrans, m, k, 1.0, Mat_buf, m, &vector_buf[m], 1, 0.0, &vector_buf[2 * m], 1); - // Compute A * g / ||g|| - // Using the forth column of vector_buff as a buffer for matrix-vector product - A(Layout::ColMajor, 1, 1.0, vector_buf, m, 0.0, &vector_buf[3*m], m); - // symv(Layout::ColMajor, uplo, m, 1.0, A, m, vector_buf, 1, 0.0, &vector_buf[3 * m], 1); - - // Compute w = (A * g / ||g|| - V * diag(eigvals) * V' * g / ||g||) - // Result is stored in the 4th column of vector_buf - blas::axpy(m, -1.0, &vector_buf[2 * m], 1, &vector_buf[3 * m], 1); - // Compute (g / ||g||)' * w - this is our measure for the error - err = blas::dot(m, vector_buf, 1, &vector_buf[3 * m], 1); - // v_0 <- v - std::copy(&vector_buf[3 * m], &vector_buf[4 * m], vector_buf); - } - return err; -} - - -template -class REVD2 { - public: - using T = typename SYRF_t::T; - using RNG = typename SYRF_t::RNG; - SYRF_t &syrf; - int error_est_p; - bool verbose; - - std::vector Y; - std::vector Omega; - std::vector R; - std::vector S; - std::vector symrf_work; - - // Constructor - REVD2( - SYRF_t &syrf_obj, - int error_est_power_iters, - bool verb = false - ) : syrf(syrf_obj) { - error_est_p = error_est_power_iters; - verbose = verb; - } - - /// Computes a rank-k approximation to an EVD of a symmetric positive semidefinite matrix: - /// A_hat = V diag(eigvals) V^*, - /// where V is a matrix of eigenvectors and eigvals is a vector of eigenvalues. - /// - /// This function is adaptive. If the tolerance is not met, increases the rank - /// estimation parameter by 2. Tolerance is the maximum of user-specified tol times - /// 5 and computed 'nu' times 5. This is motivated by the fact that the approximation - /// error will never be smaller than nu. - /// - /// This code is identical to Algorithm E2 from https://arxiv.org/pdf/2110.02820.pdf. - /// It uses a SymmetricRangeFinder for constructing a sketching operator. - /// It has a lot of potential in terms of storage space optimization, - /// which, however, will affect readability. - /// - /// @param[in] m - /// The number of rows in the matrix A. - /// - /// @param[in] A - /// The m-by-m matrix A, stored in a column-major format. - /// Must be SPD. - /// - /// @param[in] k - /// Column dimension of a sketch, k <= n. - /// - /// @param[in, out] V - /// On entry, is empty and may not have any space allocated for it. - /// On exit, stores m-by-k matrix matrix of (approximate) eigenvectors. - /// - /// @param[in, out] eigvals - /// On entry, is empty and may not have any space allocated for it. - /// On exit, stores k eigenvalues. - /// - int call( - Uplo uplo, - int64_t m, - const T* A, - int64_t &k, - T tol, - std::vector &V, - std::vector &eigvals, - RandBLAS::RNGState &state - ) { - linops::ExplicitSymLinOp A_linop(m, uplo, A, m, Layout::ColMajor); - return this->call(A_linop, k, tol, V, eigvals, state); - } - - template - int call( - SLO &A, - int64_t &k, - T tol, - std::vector &V, - std::vector &eigvals, - RandBLAS::RNGState &state - ) { - int64_t m = A.dim; - T err = 0; - RandBLAS::RNGState error_est_state(state.counter, state.key); - error_est_state.key.incr(1); - while(true) { - util::upsize(k, eigvals); - T* V_dat = util::upsize(m * k, V); - T* Y_dat = util::upsize(m * k, this->Y); - T* Omega_dat = util::upsize(m * k, this->Omega); - T* R_dat = util::upsize(k * k, this->R); - T* S_dat = util::upsize(k * k, this->S); - T* symrf_work_dat = util::upsize(m * k, this->symrf_work); - - // Construnct a sketching operator - // If CholeskyQR is used for stab/orth here, RF can fail - this->syrf.call(A, k, this->Omega, state, symrf_work_dat); - - // Y = A * Omega - A(Layout::ColMajor, k, 1.0, Omega_dat, m, 0.0, Y_dat, m); - - T nu = std::numeric_limits::epsilon() * lapack::lange(Norm::Fro, m, k, Y_dat, m); - - // We need Y = Y + v Omega - // We further need R = chol(Omega' Y) - // Solve this as R = chol(Omega' Y + v Omega'Omega) - // Compute v Omega' Omega; syrk only computes the lower triangular part. Need full. - blas::syrk(Layout::ColMajor, Uplo::Lower, Op::Trans, k, m, nu, Omega_dat, m, 0.0, R_dat, k); - for(int i = 1; i < k; ++i) - blas::copy(k - i, &R_dat[i + ((i-1) * k)], 1, &R_dat[(i - 1) + (i * k)], k); - // Compute Omega' Y + v Omega' Omega - blas::gemm(Layout::ColMajor, Op::Trans, Op::NoTrans, k, k, m, 1.0, Omega_dat, m, Y_dat, m, 1.0, R_dat, k); - - // Compute R = chol(Omega' Y + v Omega' Omega) - // Looks like if POTRF gets passed a non-triangular matrix, it will also output a non-triangular one - if(lapack::potrf(Uplo::Upper, k, R_dat, k)) - throw std::runtime_error("Cholesky decomposition failed."); - RandLAPACK::util::get_U(k, k, R_dat, k); - - // B = Y(R')^-1 - need to transpose R - blas::trsm(Layout::ColMajor, Side::Right, Uplo::Upper, Op::NoTrans, Diag::NonUnit, m, k, 1.0, R_dat, k, Y_dat, m); - - //[V, S, ~] = SVD(B) - // Although we don't need the right singular vectors, we need to give space for those. - // Use R as a buffer for that. - lapack::gesdd(Job::SomeVec, m, k, Y_dat, m, S_dat, V_dat, m, R_dat, k); - - // eigvals = diag(S^2) - T buf; - int64_t r = 0; - int i; - for(i = 0; i < k; ++i) { - buf = std::pow(S[i], 2); - eigvals[i] = buf; - // r = number of entries in eigvals that are greater than v - if(buf > nu) - ++r; - } - - // Undo regularlization - // Need to make sure no eigenvalue is negative - for(i = 0; i < r; ++i) - (eigvals[i] - nu < 0) ? 0 : eigvals[i] -=nu; - - std::fill(&V_dat[m * r], &V_dat[m * k], 0.0); - - // Error estimation - // Using the first column of Omega as a buffer for a random vector - // To perform the following safely, need to make sure Omega has at least 4 columns - Omega_dat = util::upsize(m * 4, this->Omega); - RandBLAS::DenseDist g(m, 1); - error_est_state = RandBLAS::fill_dense(g, Omega_dat, error_est_state); - - err = power_error_est(A, k, this->error_est_p, Omega_dat, V_dat, Y_dat, eigvals.data()); - - if(err <= 5 * std::max(tol, nu) || k == m) { - break; - } else if (2 * k > m) { - k = m; - } else { - k = 2 * k; - } - } - return 0; - } - -}; - - -} // end namespace RandLAPACK diff --git a/RandLAPACK/linops/rl_linops.hh b/RandLAPACK/linops/rl_linops.hh index b981f2c35..5eaf89100 100644 --- a/RandLAPACK/linops/rl_linops.hh +++ b/RandLAPACK/linops/rl_linops.hh @@ -16,3 +16,4 @@ #include "rl_composite_linop.hh" #include "rl_sym_linops.hh" #include "rl_materialize.hh" +#include "rl_matfun_linops.hh" diff --git a/RandLAPACK/linops/rl_matfun_linops.hh b/RandLAPACK/linops/rl_matfun_linops.hh new file mode 100644 index 000000000..705ac651d --- /dev/null +++ b/RandLAPACK/linops/rl_matfun_linops.hh @@ -0,0 +1,110 @@ +#pragma once + +#include "rl_blaspp.hh" +#include + +#ifdef _OPENMP +#include +#endif + +#include + +namespace RandLAPACK::linops { + + +/// Residual operator (f(A) - f(Â)) for use as a SymmetricLinearOperator. +/// +/// Computes C = α*(f(A)*B - f(Â)*B) + β*C via: +/// Z1 = f(A)*B using a matrix-function oracle (e.g. LanczosFA) +/// Z2 = f(Â)*B = V*(diag(F_vec)*(V^T B)) using two GEMMs +/// +/// Primary use: Hutchinson correction in funNyström++ (FunNystromPP driver). +///  = V diag(λ) V^T is the Nyström approximation; V (n×k) and F_vec (k, +/// holding f(λ)) are provided by the caller. +/// +/// Precision note: Z1 and Z2 follow different arithmetic paths even when f(A)=f(Â) +/// exactly (i.e. k equals the effective rank of A). Z1 uses an iterative Krylov +/// method; Z2 uses two dense GEMMs with the eigenvectors. The per-column disagreement +/// is ~1e-14, producing a systematic Hutchinson bias of roughly 1e-5 relative when +/// Phase 2's true value is zero. This is an inherent floating-point limitation of +/// combining an iterative oracle with a direct GEMM in the same residual; it is not a +/// bug. To reach machine precision at k = effective_rank, skip Phase 2 entirely. +/// +/// @tparam T Scalar type. +/// @tparam SLO_t Type of the operator A. +/// @tparam LFA_t Type of the matrix-function oracle (e.g. LanczosFA). +/// @tparam F_t Callable type T→T. +template +struct ResidualOp { + using scalar_t = T; + const int64_t dim; + + SLO_t& A; + LFA_t& matfun_oracle; // e.g. LanczosFA: computes f(A)*B + F_t f; + int64_t d; // Lanczos steps (or oracle depth) + int64_t k; // Nyström rank + T* V; // n×k Nyström eigenvectors + T* F_vec; // k-vector: f applied to Nyström eigenvalues + T* tmp; // k×s workspace (owned by caller) + T* Z1; // n×s workspace: f(A)*Ω (owned by caller) + T* Z2; // n×s workspace: f(Â)*Ω (owned by caller) + + ResidualOp(int64_t n_, SLO_t& A_, LFA_t& oracle_, F_t f_, + int64_t d_, int64_t k_, T* V_, T* Fv_, T* tmp_, T* Z1_, T* Z2_) + : dim(n_), A(A_), matfun_oracle(oracle_), f(f_), + d(d_), k(k_), V(V_), F_vec(Fv_), tmp(tmp_), Z1(Z1_), Z2(Z2_) {} + + // ------------------------------------------------------------------ + /// Apply: C = α*(f(A)*B - f(Â)*B) + β*C. + /// + /// @param[in] layout Ignored; the operator always uses ColMajor internally. + /// @param[in] n_vecs Number of right-hand sides. + /// @param[in] alpha Scalar multiplier for the residual. + /// @param[in] B n×n_vecs input matrix (col-major, ldb). + /// @param[in] ldb Leading dimension of B. + /// @param[in] beta Scalar multiplier for existing C; β=0 is handled without + /// reading C, so NaN in C does not propagate. + /// @param[in,out] C n×n_vecs output matrix (col-major, ldc); overwritten. + /// @param[in] ldc Leading dimension of C; may exceed dim. + void operator()([[maybe_unused]] Layout layout, int64_t n_vecs, T alpha, + T* const B, int64_t ldb, T beta, T* C, int64_t ldc) { + if (ldb != dim) + throw std::invalid_argument("ResidualOp: B must be contiguous (ldb == dim)"); + // Z1 = f(A)*B via matrix-function oracle + matfun_oracle.call(A, B, dim, n_vecs, f, d, Z1); + + // Z2 = f(Â)*B = V * (diag(F_vec) * (V^T B)) + blas::gemm(Layout::ColMajor, Op::Trans, Op::NoTrans, + k, n_vecs, dim, (T)1.0, V, dim, B, ldb, (T)0.0, tmp, k); +#ifdef _OPENMP + #pragma omp parallel for schedule(static) +#endif + for (int64_t j = 0; j < n_vecs; ++j) + for (int64_t i = 0; i < k; ++i) + tmp[j * k + i] *= F_vec[i]; + blas::gemm(Layout::ColMajor, Op::NoTrans, Op::NoTrans, + dim, n_vecs, k, (T)1.0, V, dim, tmp, k, (T)0.0, Z2, dim); + + // C = alpha*(Z1 - Z2) + beta*C column-by-column to respect ldc. + // Z1 and Z2 are always stored flat (stride dim); C may have ldc >= dim. + // beta==0 uses a fused read of Z1,Z2 to avoid scal(0,C) propagating NaN. + if (beta == (T)0) { +#ifdef _OPENMP + #pragma omp parallel for schedule(static) +#endif + for (int64_t j = 0; j < n_vecs; ++j) + for (int64_t i = 0; i < dim; ++i) + C[j * ldc + i] = alpha * (Z1[j * dim + i] - Z2[j * dim + i]); + } else { + for (int64_t j = 0; j < n_vecs; ++j) { + blas::scal(dim, beta, C + j * ldc, 1); + blas::axpy(dim, alpha, Z1 + j * dim, 1, C + j * ldc, 1); + blas::axpy(dim, -alpha, Z2 + j * dim, 1, C + j * ldc, 1); + } + } + } +}; + + +} // end namespace RandLAPACK::linops diff --git a/RandLAPACK/linops/rl_sparse_linop.hh b/RandLAPACK/linops/rl_sparse_linop.hh index 108c69827..5d1b707e9 100644 --- a/RandLAPACK/linops/rl_sparse_linop.hh +++ b/RandLAPACK/linops/rl_sparse_linop.hh @@ -464,4 +464,47 @@ public: } }; + +/*********************************************************/ +/* */ +/* SparseSymLinOp */ +/* */ +/*********************************************************/ + +/// @brief Symmetric sparse linear operator satisfying the SymmetricLinearOperator concept. +/// +/// Wraps a sparse matrix assumed to be symmetric with both triangles explicitly stored. +/// Delegates to `left_spmm` with NoTrans for SYMM-like semantics: C = alpha * A * B + beta * C. +/// +/// Satisfies SymmetricLinearOperator (provides `dim` and SYMM-like operator()). +/// Also provides `n_rows` and `n_cols` (both equal to `dim`) for LinearOperator compatibility. +/// +/// @tparam SpMat Sparse matrix type (CSR, CSC, or COO format). +/// +/// @note The sparse matrix must be square (n×n) and symmetric with both triangles stored. +/// The caller's SpMat must outlive this operator. +template +struct SparseSymLinOp { + using T = typename SpMat::scalar_t; + using scalar_t = T; + const int64_t dim; + const int64_t n_rows; + const int64_t n_cols; + SpMat& A_sp; + + /// @brief Construct from an existing symmetric sparse matrix. + /// @param n Matrix dimension (n × n); A_sp must be n × n. + /// @param src Sparse matrix; must be square, symmetric, both triangles stored. + SparseSymLinOp(int64_t n, SpMat& src) + : dim(n), n_rows(n), n_cols(n), A_sp(src) {} + + /// @brief Apply: C = alpha * A_sp * B + beta * C. + void operator()(Layout layout, int64_t n_vecs, T alpha, + T* const B, int64_t ldb, + T beta, T* C, int64_t ldc) { + RandBLAS::sparse_data::left_spmm(layout, Op::NoTrans, Op::NoTrans, + dim, n_vecs, dim, alpha, A_sp, 0, 0, B, ldb, beta, C, ldc); + } +}; + } // end namespace RandLAPACK::linops diff --git a/RandLAPACK/linops/rl_sym_linops.hh b/RandLAPACK/linops/rl_sym_linops.hh index 524774ab2..db31e3ebd 100644 --- a/RandLAPACK/linops/rl_sym_linops.hh +++ b/RandLAPACK/linops/rl_sym_linops.hh @@ -104,6 +104,51 @@ struct ExplicitSymLinOp { return A_buff[i + j*lda]; } } + + /// SkOp overload: apply A to a SketchingOperator S (sparse or dense). + /// Computes C = alpha * A * S + beta * C. + /// + /// Mirrors the pattern in `SparseLinOp::operator()(SkOp&)`: + /// - Dense SkOp → extract S.buff and dispatch to the existing SYMM path. + /// - Sparse SkOp → use RandBLAS::sparse_data::right_spmm directly, + /// bypassing the densification of S into a dense buffer. + /// + /// PRECONDITION (sparse path only): A_buff must have BOTH triangles populated. + /// `right_spmm` reads A as a generic dense matrix and does not exploit + /// symmetry — so it will read whatever is in the lower triangle. For + /// matrices generated via `gen_sym_psd_lowrank` (which only fills upper), + /// the caller must symmetrize before constructing this operator. + /// `SparseSymLinOp` imposes the same constraint by convention. + template + void operator()( + Layout layout, + int64_t n_vecs, + T alpha, + SkOp& S, + T beta, + T* C, + int64_t ldc + ) { + if constexpr (requires { S.buff; S.layout; S.dist; }) { + // Dense sketch — extract buffer, dispatch to existing dense matvec. + // The existing path uses blas::symm and works fine with one-triangle storage. + if (S.buff == nullptr) RandBLAS::fill_dense(S); + int64_t ldS = S.dist.dim_major; + randblas_require(S.layout == layout); // simplification — caller ensures match + (*this)(layout, n_vecs, alpha, S.buff, ldS, beta, C, ldc); + } else { + // Sparse sketch — dense·sparse via right_spmm. Requires both triangles of A. + if (S.nnz < 0) RandBLAS::fill_sparse(S); + auto S_coo = RandBLAS::coo_view_of_skop(S); + RandBLAS::sparse_data::right_spmm( + layout, Op::NoTrans, Op::NoTrans, + dim, n_vecs, dim, + alpha, this->A_buff, this->lda, + S_coo, 0, 0, + beta, C, ldc + ); + } + } }; /*********************************************************/ diff --git a/RandLAPACK/misc/rl_util.hh b/RandLAPACK/misc/rl_util.hh index 0c8012d6e..6e57676f2 100644 --- a/RandLAPACK/misc/rl_util.hh +++ b/RandLAPACK/misc/rl_util.hh @@ -124,6 +124,22 @@ void get_U( } } +/// In-place symmetrize: replace A[i,j] and A[j,i] (i ≠ j) with their +/// average. Used to clean up small finite-arithmetic asymmetry in matrices +/// that should be symmetric in exact arithmetic (e.g. Ωᵀ·A·Ω built via two +/// non-trivial GEMMs, or block-Lanczos diagonal blocks Qᵢᵀ·A·Qᵢ). +/// Column-major; touches the strict upper and strict lower triangles. +template +void symmetrize(int64_t n, T* A, int64_t lda) { + for (int64_t j = 0; j < n; ++j) { + for (int64_t i = j + 1; i < n; ++i) { + T avg = (T)0.5 * (A[i + j * lda] + A[j + i * lda]); + A[i + j * lda] = avg; + A[j + i * lda] = avg; + } + } +} + /// Returns true if all diagonal entries of an upper-triangular matrix are /// nonzero (safe to use in TRSM). Returns false if any diagonal is exactly /// zero, which would cause division by zero during triangular solve. @@ -186,10 +202,10 @@ void col_swap( } } -/// Checks if the given size is larger than available. -/// If so, resizes the vector. +/// Grow a std::vector to at least target_sz elements (zero-fills new entries). +/// Returns a raw pointer to the underlying data. template -T* upsize( +T* resize( int64_t target_sz, std::vector &A ) { @@ -199,6 +215,19 @@ T* upsize( return A.data(); } +/// Grow a raw buffer to at least `needed` elements. +/// Replaces the allocation; existing contents are not preserved. +/// New storage is uninitialized — callers must write before reading. +template +void resize(T*& buf, int64_t& buf_sz, int64_t needed) { + if (needed > buf_sz) { + delete[] buf; + buf = new T[needed]; + buf_sz = needed; + } +} + + /// Uses recursion to find the rank of the matrix pointed to by A_dat. /// Does so by attempting to find the smallest k such that /// ||A[k:, k:]||_F <= tau_trunk * ||A||. @@ -246,7 +275,7 @@ void normc( const std::vector &A, std::vector &A_norm ) { - util::upsize(m * n, A_norm); + util::resize(m * n, A_norm); T col_nrm = 0.0; for(int i = 0; i < n; ++i) { diff --git a/RandLAPACK/testing/rl_gen.hh b/RandLAPACK/testing/rl_gen.hh index f70c94482..a3b94c5bc 100644 --- a/RandLAPACK/testing/rl_gen.hh +++ b/RandLAPACK/testing/rl_gen.hh @@ -3,6 +3,10 @@ #include "rl_blaspp.hh" #include "rl_lapackpp.hh" +#ifdef _OPENMP +#include +#endif + #include #include #include @@ -885,4 +889,176 @@ std::vector gen_poly_mat_psd(int64_t m, T cond_num, T exponent, uint32_t seed return G; } +/// Write k algebraically-decaying spectral values into s[0..k-1]: +/// s[i] = scale * (i+1)^{-exponent}, i = 0, 1, …, k-1. +/// Matches Persson & Kressner (SIMAX 2023) Figure 5(a): scale=100, exponent=2. +/// +/// @param[in] k Number of values to generate. +/// @param[in] scale Overall scale factor. +/// @param[in] exponent Decay exponent (positive). +/// @param[out] s Output buffer of length k; overwritten. +template +void gen_alg_decay_singvals(int64_t k, T scale, T exponent, T* s) { + #ifdef _OPENMP + #pragma omp parallel for schedule(static) + #endif + for (int64_t i = 0; i < k; ++i) + s[i] = scale * std::pow((T)(i + 1), -exponent); +} + +/// Write k exponentially-decaying spectral values into s[0..k-1]: +/// s[i] = scale * exp(-(i+1)/decay_rate), i = 0, 1, …, k-1. +/// Matches Persson & Kressner (SIMAX 2023) Figure 5(b): scale=1, decay_rate=100. +/// +/// @param[in] k Number of values to generate. +/// @param[in] scale Overall scale factor. +/// @param[in] decay_rate Rate parameter (larger = slower decay). +/// @param[out] s Output buffer of length k; overwritten. +template +void gen_exp_decay_singvals(int64_t k, T scale, T decay_rate, T* s) { + #ifdef _OPENMP + #pragma omp parallel for schedule(static) + #endif + for (int64_t i = 0; i < k; ++i) + s[i] = scale * std::exp(-(T)(i + 1) / decay_rate); +} + +/// Form the upper triangle of A = Q diag(eigvals) Q' with a Haar-random n×k matrix Q. +/// +/// Q is obtained via QR factorization of an n×k Gaussian matrix. Column j of Q is +/// scaled by sqrt(eigvals[j]) in-place — no k×k spectrum matrix is formed. The upper +/// triangle of A is then filled via syrk. The n - k trailing eigenvalues are implicitly zero. +/// +/// @param[in] n Matrix dimension (n × n output) +/// @param[in] k Rank; k ≤ n +/// @param[out] A Output buffer, upper triangle on exit (col-major, leading dim lda) +/// @param[in] lda Leading dimension of A +/// @param[in] eigvals Pointer to k non-negative eigenvalues +/// @param[in,out] state RNG state +template +void gen_sym_psd_lowrank( + int64_t n, int64_t k, + T* A, int64_t lda, + const T* eigvals, + RandBLAS::RNGState& state +) { + T* Q = new T[n * k](); + T* tau = new T[k](); + + auto dist = RandBLAS::DenseDist(n, k); + state = RandBLAS::fill_dense(dist, Q, state); + + lapack::geqrf(n, k, Q, n, tau); + lapack::orgqr(n, k, k, Q, n, tau); + + #ifdef _OPENMP + #pragma omp parallel for schedule(static) + #endif + for (int64_t j = 0; j < k; ++j) + blas::scal(n, std::sqrt(eigvals[j]), Q + j * n, 1); + + blas::syrk(Layout::ColMajor, Uplo::Upper, Op::NoTrans, + n, k, (T)1.0, Q, n, (T)0.0, A, lda); + + delete[] Q; + delete[] tau; +} + +/// Form an n×n symmetric kernel matrix K from n random data points in R^d. +/// +/// Data points x_1,...,x_n are drawn i.i.d. from N(0, I_d). The upper +/// triangle of K is filled with K[i,j] = k(x_i, x_j); the lower triangle +/// is left unset (use Uplo::Upper with LAPACK symmetric routines). +/// +/// Pairwise dot products are computed via a single SYRK call (BLAS-3, +/// multi-threaded via MKL/OpenBLAS), so assembly is O(n²d) dominated. +/// +/// Two kernels are supported via kernel_type: +/// 0 = RBF: k(x,y) = exp(-||x-y||² / (2*bandwidth²)) +/// 1 = Polynomial: k(x,y) = (x·y + poly_c)^poly_p +/// +/// For RBF with Gaussian data, the expected squared distance is 2d, so +/// bandwidth = sqrt(d) (median heuristic) gives average kernel value ~exp(-1). +/// +/// @param[in] n Matrix dimension; also the number of data points +/// @param[in] d Data point dimension (ambient space) +/// @param[out] K Output buffer; upper triangle on exit (col-major, ldk) +/// @param[in] ldk Leading dimension of K +/// @param[in] kernel_type 0=RBF, 1=polynomial +/// @param[in] bandwidth RBF length scale σ; ignored for polynomial +/// @param[in] poly_c Polynomial additive constant c; ignored for RBF +/// @param[in] poly_p Polynomial degree p; ignored for RBF +/// @param[in,out] state RNG state; advanced on return +template +void gen_kernel_matrix( + int64_t n, int64_t d, + T* K, int64_t ldk, + int kernel_type, + T bandwidth, + T poly_c, + int poly_p, + RandBLAS::RNGState& state +) { + // X: n×d col-major; row i is data point x_i ~ N(0, I_d) + T* X = new T[n * d]; + state = RandBLAS::fill_dense(RandBLAS::DenseDist(n, d), X, state); + + // G = X*X^T (upper triangle, col-major): G[i + j*n] = x_i · x_j (i <= j) + T* G = new T[n * n](); + blas::syrk(Layout::ColMajor, Uplo::Upper, Op::NoTrans, + n, d, (T)1.0, X, n, (T)0.0, G, n); + delete[] X; + + if (kernel_type == 0) { + T inv_denom = (T)1.0 / ((T)2.0 * bandwidth * bandwidth); + #ifdef _OPENMP + #pragma omp parallel for schedule(static) + #endif + for (int64_t j = 0; j < n; ++j) { + for (int64_t i = 0; i <= j; ++i) { + T sq_dist = G[i + i*n] - (T)2.0 * G[i + j*n] + G[j + j*n]; + K[i + j*ldk] = std::exp(-sq_dist * inv_denom); + } + } + } else { + #ifdef _OPENMP + #pragma omp parallel for schedule(static) + #endif + for (int64_t j = 0; j < n; ++j) { + for (int64_t i = 0; i <= j; ++i) + K[i + j*ldk] = std::pow(G[i + j*n] + poly_c, (T)poly_p); + } + } + delete[] G; +} + +/// Form A = Q diag(eigvals) Q' + noise * I, upper triangle. +/// +/// Calls gen_sym_psd_lowrank for the low-rank component, then adds `noise` +/// to each diagonal entry. When noise > 0, every eigenvalue is shifted by +/// noise, making A strictly PD — required for f = log (λ_min > 0). +/// +/// Useful for log det(I + K) benchmarks where A = I + K: set noise = 1 and +/// use the kernel eigenvalues as `eigvals`. +/// +/// @param[in] n Matrix dimension (n × n output) +/// @param[in] k Rank of the low-rank component; k ≤ n +/// @param[out] A Output buffer, upper triangle on exit (col-major, leading dim lda) +/// @param[in] lda Leading dimension of A +/// @param[in] eigvals Pointer to k non-negative eigenvalues for the low-rank part +/// @param[in] noise Diagonal shift (≥ 0). Use 0 for no shift. +/// @param[in,out] state RNG state +template +void gen_shifted_lowrank_psd( + int64_t n, int64_t k, + T* A, int64_t lda, + const T* eigvals, + T noise, + RandBLAS::RNGState& state +) { + gen_sym_psd_lowrank(n, k, A, lda, eigvals, state); + for (int64_t i = 0; i < n; ++i) + A[i * lda + i] += noise; +} + } diff --git a/benchmark/CMakeLists.txt b/benchmark/CMakeLists.txt index 6aa7749d1..5d7969e27 100644 --- a/benchmark/CMakeLists.txt +++ b/benchmark/CMakeLists.txt @@ -165,6 +165,9 @@ add_benchmark(NAME ABRIK_runtime_breakdown_sparse CXX_SOURCES bench_ABRIK/ABRIK_ add_benchmark(NAME ABRIK_speed_comparisons CXX_SOURCES bench_ABRIK/ABRIK_speed_comparisons.cc LINK_LIBS ${Benchmark_libs_external}) add_benchmark(NAME ABRIK_speed_comparisons_sparse CXX_SOURCES bench_ABRIK/ABRIK_speed_comparisons_sparse.cc LINK_LIBS ${Benchmark_libs_external}) +# FunNystromPP benchmarks +add_benchmark(NAME FunNystromPP_benchmark CXX_SOURCES bench_FunNystromPP/FunNystromPP_benchmark.cc LINK_LIBS ${Benchmark_libs}) + # GPU benchmarks (if CUDA is available) if (HAVE_CUDA) message(STATUS "Building GPU benchmarks...") diff --git a/benchmark/bench_BQRRP/BQRRP_error_analysis.cc b/benchmark/bench_BQRRP/BQRRP_error_analysis.cc index 3060103b0..14299bd05 100644 --- a/benchmark/bench_BQRRP/BQRRP_error_analysis.cc +++ b/benchmark/bench_BQRRP/BQRRP_error_analysis.cc @@ -70,7 +70,7 @@ error_check(QR_benchmark_data &all_data, auto n = all_data.col; auto k = all_data.rank; - RandLAPACK::util::upsize(k * k, all_data.I_ref); + RandLAPACK::util::resize(k * k, all_data.I_ref); RandLAPACK::util::eye(k, k, all_data.I_ref); T* A_dat = all_data.A_cpy1.data(); diff --git a/benchmark/bench_CQRRPT/CQRRPT_error_analysis.cc b/benchmark/bench_CQRRPT/CQRRPT_error_analysis.cc index 2f7c816c6..a12a3877e 100644 --- a/benchmark/bench_CQRRPT/CQRRPT_error_analysis.cc +++ b/benchmark/bench_CQRRPT/CQRRPT_error_analysis.cc @@ -72,7 +72,7 @@ error_check(QR_benchmark_data &all_data, auto n = col_sz; auto k = all_data.rank; - RandLAPACK::util::upsize(k * k, all_data.I_ref); + RandLAPACK::util::resize(k * k, all_data.I_ref); RandLAPACK::util::eye(k, k, all_data.I_ref); T* A_dat = all_data.A_cpy1.data(); diff --git a/benchmark/bench_FunNystromPP/FunNystromPP_benchmark.cc b/benchmark/bench_FunNystromPP/FunNystromPP_benchmark.cc new file mode 100644 index 000000000..ae9f3f239 --- /dev/null +++ b/benchmark/bench_FunNystromPP/FunNystromPP_benchmark.cc @@ -0,0 +1,826 @@ +#if defined(__APPLE__) +int main() {return 0;} +#else +/* +FunNystromPP benchmark -- compares funNystrom++ against plain Hutchinson+LanczosFA +for estimating tr(f(A)). + +Two algorithms are timed and their matvec counts recorded for each matrix size: + 1. FunNystromPP(k, s, d): Nystrom rank-k approximation plus Hutchinson correction + with s samples, each using d Lanczos steps. + 2. Hutchinson+LanczosFA(s, d): plain Hutchinson with s Rademacher samples, + each using LanczosFA with d steps to apply f(A). No Nystrom component. + Uses the same s and d as FunNystromPP (equal correction-phase budget). + +Both algorithms use the Lanczos variant selected by lfa_type (see below). + +Matrix input (mat_or_file): + 1 = psd_alg: A = Q diag(λ) Q', λᵢ = 100 * i^{-2} (algebraic decay, κ ≈ n²/100) + 2 = psd_exp: A = Q diag(λ) Q', λᵢ = exp(-i/100) (exponential decay) + 3 = rbf_kernel: K[i,j] = exp(-||x_i-x_j||²/(2σ²)), x_i~N(0,I_d), σ=sqrt(d) + 4 = poly_kernel: K[i,j] = (x_i·x_j + 1)², x_i~N(0,I_d) + /path.txt = load dense square symmetric matrix from space-separated text file + (row-major, full matrix stored; upper triangle used by algorithms) + /path.mtx = load sparse symmetric matrix in Matrix Market format + (coordinate real symmetric; both triangles expanded internally) + + For types 1/2: k_mat eigenvalues are nonzero; exact tr(f(A)) is always available. + For types 3/4: k_mat_ratio reinterpreted as data dimension d = n/k_mat_ratio; + reference requires compute_ref=1 (syevd, only feasible for small n). + +Function types (func_type): + 0 = sqrt(x) [f_zero = 0; safe for SPSD matrices] + 1 = log(x) [requires all eigenvalues > 0; for generated types, noise=1 is applied + automatically (A = lowrank + I), so all eigenvalues >= 1] + 2 = x*(x + poly_lambda) [degree-2 polynomial; exact in d=2 Lanczos steps; use for tr(K(K+λI))] + +Lanczos variant (lfa_type): + 0 = scalar LanczosFA (s independent Krylov sequences; BLAS-1/2 dominant) + 1 = BlockLanczosFA (single joint block Krylov subspace; BLAS-3 dominant) + +SYRF orthogonalization (orth_type): + 0 = HouseholderQR (HQRQ) — backward-stable (geqrf+ungqr); ~2·n·k² flops + 1 = CholeskyQR (CholQRQ) — ~2× faster (syrk+potrf+trsm); fails on + ill-conditioned input (the SYRF output after + power iteration is generally well-conditioned). + +Sketch types (sketch_type): + 0 = Gaussian (DenseDist, default) + 1 = SASO (SparseDist; nonzeros per column set by vec_nnz argument) + +compute_ref: + 0 = skip reference (report -1 in output) + 1 = run syevd to compute exact tr(f(A)) + For types 1/2 this flag is ignored; the reference is always available from eigenvalues. + For types 3/4 and file input, syevd is run only when compute_ref=1. + +precision: + double = run in double precision (default) + float = run in single precision (sparse path is always double regardless) + +Both algorithms use a fixed seed so comparisons are reproducible. +Matvec counts reflect actual A applications inside each algorithm (not estimates). +The matrix is constructed/loaded once per n; that cost is not timed. + +Output format: 32 comma-separated columns per row (trailing comma): + n, k, k_mat, s, d, + matvec_pp, matvec_hutch, + time_pp_total (us), time_pp_phase1 (us), time_pp_phase2 (us), + time_hutch (us), + nystrom_alloc, nystrom_syrf, nystrom_matvec, nystrom_gram, + nystrom_potrf, nystrom_trsm, nystrom_svd, nystrom_post_svd, + nystrom_error_est, nystrom_rest, nystrom_total (all us), + lfa_matvec, lfa_run_lanczos, lfa_apply_f, lfa_rest, lfa_total (all us), + est_fun_pp, est_hutch, + true_tr (-1 if unavailable), err_fun_pp, err_hutch + + lfa_type is encoded in the output filename: _scalar_ or _block_. + +Usage: + FunNystromPP_benchmark + + [n_1 n_2 ...] + dir_path : output directory (use "." for current directory) + num_runs : number of timed repetitions per matrix size + k : Nystrom rank (constant; must satisfy k <= n) + k_mat : matrix rank for generated types (constant; for kernel types equals data dim d; + ignored for file input) + s : Hutchinson sample count (constant) + d_steps : Lanczos depth per Hutchinson sample + mat_or_file : 1=psd_alg, 2=psd_exp, 3=rbf_kernel, 4=poly_kernel, or path to .txt/.mtx file + func_type : 0=sqrt, 1=log, 2=poly (x*(x+poly_lambda)) + compute_ref : 0=skip eigendecomp reference, 1=run syevd (file input only) + sketch_type : 0=Gaussian, 1=SASO + vec_nnz : nonzeros per column for SASO sketch (sketch_type=1; ignored for Gaussian) + poly_lambda : lambda parameter for func_type=2 (ignored for func_type=0,1) + precision : double or float (sparse .mtx path is always double) + lfa_type : 0=scalar LanczosFA, 1=BlockLanczosFA + n_1 ... : matrix dimensions (generated types only; ignored for file input) +*/ + +#include "RandLAPACK.hh" +#include "rl_blaspp.hh" +#include "rl_lapackpp.hh" +#include "rl_gen.hh" +#include "load_mtx.hh" + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace linops = RandLAPACK::linops; +using RNG = r123::Philox4x32; +using std::chrono::steady_clock; +using std::chrono::duration_cast; +using std::chrono::microseconds; + +// ============================================================================ +// Benchmark data struct +// ============================================================================ + +template +struct FunNystromPP_benchmark_data { + int64_t n; + int64_t k; + int64_t k_mat; + int64_t s; + int64_t d; + bool from_file; + T noise; + int64_t d_dim; + std::vector A; + std::vector eigvals; + + FunNystromPP_benchmark_data(int64_t n_max, int64_t k_, int64_t k_mat_, + int64_t s_, int64_t d_, bool from_file_) : + n(n_max), k(k_), k_mat(k_mat_), s(s_), d(d_), from_file(from_file_), + noise((T)0.0), d_dim(0), A(n_max * n_max, 0.0) + {} +}; + +// ============================================================================ +// CountingLinOp: wraps any SLO and counts individual matrix-vector products. +// ============================================================================ + +template +struct CountingLinOp { + using scalar_t = typename SLO_t::scalar_t; + const int64_t dim; + SLO_t& base; + int64_t count = 0; + + CountingLinOp(int64_t n, SLO_t& slo) : dim(n), base(slo) {} + void reset() { count = 0; } + + void operator()(Layout layout, int64_t n_vecs, scalar_t alpha, + scalar_t* const B, int64_t ldb, + scalar_t beta, scalar_t* C, int64_t ldc) { + count += n_vecs; + base(layout, n_vecs, alpha, B, ldb, beta, C, ldc); + } +}; + +// ============================================================================ +// data_regen: build or load A for the current n. +// +// Optional disk cache: generated matrices (the synthetic-decay and kernel +// paths) are saved to /tmp/funny_A_mat_n_kmat_func_.bin +// after first generation. Subsequent invocations with the same parameters +// load from the cache instead of regenerating, saving the Haar-random QR + +// scale + syrk cost on every benchmark invocation. +// +// The cache file is just a flat binary dump of A (n*n entries, col-major). +// Eigenvalues, noise, and d_dim are recomputed deterministically each time — +// they're cheap to redo and avoid format/version concerns in the cache file. +// ============================================================================ + +template +static void data_regen( + FunNystromPP_benchmark_data& data, + RandBLAS::RNGState& state, + int mat_type, + int func_type, + const std::string& file_path +) { + int64_t n = data.n; + std::fill(data.A.begin(), data.A.begin() + n * n, (T)0.0); + + if (data.from_file) { + int64_t m = n, ncols = n; + int wq = 0; + RandLAPACK::gen::process_input_mat( + m, ncols, data.A.data(), + const_cast(file_path.c_str()), wq); + data.eigvals.clear(); + data.noise = (T)0.0; + data.d_dim = 0; + return; + } + + // ------ Set up small bookkeeping (eigvals / noise / d_dim) — deterministic + // and cheap, recomputed every call regardless of cache hit. ------ + if (mat_type >= 3) { + data.d_dim = data.k_mat; + data.noise = (func_type == 1) ? (T)1.0 : (T)0.0; + data.eigvals.clear(); + } else { + data.d_dim = 0; + data.eigvals.resize(data.k_mat); + if (mat_type == 1) + RandLAPACK::gen::gen_alg_decay_singvals(data.k_mat, (T)100.0, (T)2.0, + data.eigvals.data()); + else + RandLAPACK::gen::gen_exp_decay_singvals(data.k_mat, (T)1.0, (T)100.0, + data.eigvals.data()); + data.noise = (func_type == 1) ? (T)1.0 : (T)0.0; + } + + // ------ Try the matrix file cache. ------ + char cache_path[256]; + std::snprintf(cache_path, sizeof(cache_path), + "/tmp/funny_A_mat%d_n%lld_kmat%lld_func%d_%c.bin", + mat_type, (long long)n, (long long)data.k_mat, func_type, + (sizeof(T) == 8) ? 'd' : 'f'); + { + std::ifstream cache_in(cache_path, std::ios::binary); + if (cache_in.good()) { + cache_in.read(reinterpret_cast(data.A.data()), + static_cast(n) * n * sizeof(T)); + if (cache_in.gcount() == + static_cast(n) * n * sizeof(T)) { + return; // cache hit + } + } + } + + // ------ Cache miss: generate. ------ + if (mat_type >= 3) { + T bandwidth = std::sqrt((T)data.d_dim); + if (mat_type == 3) + RandLAPACK::gen::gen_kernel_matrix( + n, data.d_dim, data.A.data(), n, 0, bandwidth, (T)0.0, 0, state); + else + RandLAPACK::gen::gen_kernel_matrix( + n, data.d_dim, data.A.data(), n, 1, (T)0.0, (T)1.0, 2, state); + if (data.noise > (T)0) { + for (int64_t i = 0; i < n; ++i) + data.A[i + i * n] += data.noise; + } + } else { + if (func_type == 1) { + RandLAPACK::gen::gen_shifted_lowrank_psd(n, data.k_mat, data.A.data(), n, + data.eigvals.data(), data.noise, state); + } else { + RandLAPACK::gen::gen_sym_psd_lowrank(n, data.k_mat, data.A.data(), n, + data.eigvals.data(), state); + } + } + + // Mirror upper into lower for the SkOp sparse-spmm path + // (right_spmm reads A as a generic dense matrix; symmetry isn't exploited). + if (mat_type < 3) { + for (int64_t j = 0; j < n; ++j) + for (int64_t i = j + 1; i < n; ++i) + data.A[i + j * n] = data.A[j + i * n]; + } + + // ------ Save the freshly generated matrix to cache for next time. ------ + std::ofstream cache_out(cache_path, std::ios::binary | std::ios::trunc); + if (cache_out.good()) { + cache_out.write(reinterpret_cast(data.A.data()), + static_cast(n) * n * sizeof(T)); + } +} + +// ============================================================================ +// LFAOp: wraps LanczosFA as a SymmetricLinearOperator for plain Hutchinson. +// ============================================================================ + +template +struct LFAOp { + using scalar_t = typename SLO_t::scalar_t; + const int64_t dim; + SLO_t& A; + LFA_t& lfa; + F_t f; + int64_t d; + std::vector Z_buf; + + LFAOp(int64_t n, SLO_t& A_, LFA_t& lfa_, F_t f_, int64_t d_) + : dim(n), A(A_), lfa(lfa_), f(f_), d(d_), Z_buf() {} + + void operator()([[maybe_unused]] Layout layout, int64_t n_vecs, scalar_t alpha, + scalar_t* const B, [[maybe_unused]] int64_t ldb, + scalar_t beta, scalar_t* C, [[maybe_unused]] int64_t ldc) + { + Z_buf.assign(dim * n_vecs, (scalar_t)0); + lfa.call(A, B, dim, n_vecs, f, d, Z_buf.data()); + if (beta == (scalar_t)0) { + for (int64_t i = 0; i < dim * n_vecs; ++i) + C[i] = alpha * Z_buf[i]; + } else { + blas::scal(dim * n_vecs, beta, C, 1); + blas::axpy(dim * n_vecs, alpha, Z_buf.data(), 1, C, 1); + } + } +}; + +// ============================================================================ +// call_all_algs: timed + matvec-counted comparison for one matrix size. +// ============================================================================ + +template +static void call_all_algs( + int64_t numruns, + FunNystromPP_benchmark_data& data, + const RandBLAS::RNGState& state, + int mat_type, + int func_type, + double poly_lambda, + bool compute_ref, + int sketch_type, + int vec_nnz, + bool skip_hutchinson, + const std::string& output_filename +) { + int64_t n = data.n, k = data.k, s = data.s, d = data.d; + + std::function f; + const char* fname; + T f_zero = (T)0; + if (func_type == 0) { + f = [](T x){ return (T)std::sqrt((double)std::max(x, (T)0)); }; + fname = "sqrt"; + } else if (func_type == 1) { + f = [](T x){ return (T)std::log((double)x); }; + fname = "log"; + } else { + T lam = (T)poly_lambda; + f = [lam](T x){ return x * (x + lam); }; + fname = "poly"; + } + + bool ref_available = false; + T true_tr = (T)0; + if (!data.from_file && mat_type < 3) { + for (int64_t j = 0; j < data.k_mat; ++j) + true_tr += f(data.eigvals[j] + data.noise); + true_tr += (T)(n - data.k_mat) * f(data.noise); + ref_available = true; + } else if (compute_ref) { + std::vector A_copy(data.A.begin(), data.A.begin() + n * n); + std::vector ev(n); + lapack::syevd(lapack::Job::NoVec, lapack::Uplo::Upper, n, + A_copy.data(), n, ev.data()); + for (int64_t j = 0; j < n; ++j) + true_tr += f(ev[j]); + ref_available = true; + } + + RandLAPACK::SYPS syps(2, 1, false, false); + syps.sketch_type = sketch_type; + syps.vec_nnz = vec_nnz; + ORTH_t orth(false, false); + // Wire ORTH_t into SYPS's internal power-iter stabilization too — same + // user choice applies. SYPS auto-falls-back to HQRQ on CholQR failure. + ORTH_t orth_inner(false, false); + syps.internal_stab = &orth_inner; + RandLAPACK::SYRF, ORTH_t> syrf(syps, orth); + RandLAPACK::NystromEVD nystrom_evd(syrf, 0); + nystrom_evd.timing = true; + LFA_t lfa_pp; + lfa_pp.timing = true; + RandLAPACK::Hutchinson hutch_pp; + RandLAPACK::FunNystromPP + driver(nystrom_evd, lfa_pp, hutch_pp); + + LFA_t lfa_base; + RandLAPACK::Hutchinson hutch_base; + + using ESLO = linops::ExplicitSymLinOp; + using CSLO = CountingLinOp; + using FuncT = std::function; + + const char* mat_str = data.from_file ? "file" : + mat_type == 1 ? "psd_alg" : + mat_type == 2 ? "psd_exp" : + mat_type == 3 ? "rbf_kernel" : "poly_kernel"; + + for (int i = 0; i < numruns; ++i) { + printf("ITERATION %d, N=%lld, K=%lld, K_MAT=%lld, S=%lld, D=%lld [%s, %s]\n", + i, (long long)n, (long long)k, (long long)data.k_mat, + (long long)s, (long long)d, fname, mat_str); + + auto state_alg = state; + + // FunNystromPP + ESLO A_op(n, blas::Uplo::Upper, data.A.data(), n, Layout::ColMajor); + CSLO A_pp(n, A_op); + + RandLAPACK::FunNystromPP_timing pp_timing; + auto start_pp = steady_clock::now(); + T est_pp = driver.call(A_pp, f, f_zero, k, s, d, state_alg, &pp_timing); + auto stop_pp = steady_clock::now(); + long dur_pp = duration_cast(stop_pp - start_pp).count(); + int64_t mv_pp = A_pp.count; + + auto& nt = nystrom_evd.times; + auto& lt = lfa_pp.times; + + printf(" FunNystromPP: est=%.6e time=%lld us (ph1=%lld ph2=%lld) matvecs=%lld\n", + (double)est_pp, (long long)dur_pp, + (long long)pp_timing.phase1_us, (long long)pp_timing.phase2_us, + (long long)mv_pp); + printf(" NystromEVD: alloc=%lld syrf=%lld matvec=%lld gram=%lld potrf=%lld trsm=%lld svd=%lld post_svd=%lld err_est=%lld rest=%lld total=%lld\n", + (long long)nt[0], (long long)nt[1], (long long)nt[2], (long long)nt[3], + (long long)nt[4], (long long)nt[5], (long long)nt[6], (long long)nt[7], + (long long)nt[8], (long long)nt[9], (long long)nt[10]); + printf(" LanczosFA: matvec=%lld lanczos=%lld apply_f=%lld total=%lld\n", + (long long)lt[0], (long long)lt[1], (long long)lt[2], (long long)lt[4]); + + // Hutchinson + LanczosFA (skipped when --skip-hutchinson is set; + // useful for s-sweeps that only need the Hutchinson reference at one + // (lfa, d) variant — e.g. scalar d=200 — and not at every sweep + // point. CSV output preserves the column shape with zeros. + T est_h = (T)0; + long dur_h = 0; + int64_t mv_h = 0; + if (!skip_hutchinson) { + ESLO A_op2(n, blas::Uplo::Upper, data.A.data(), n, Layout::ColMajor); + CSLO A_hutch(n, A_op2); + LFAOp lfa_op(n, A_hutch, lfa_base, f, d); + + auto start_h = steady_clock::now(); + est_h = hutch_base.call(lfa_op, s, state_alg); + auto stop_h = steady_clock::now(); + dur_h = duration_cast(stop_h - start_h).count(); + mv_h = A_hutch.count; + printf(" Hutchinson+LFA: est=%.6e time=%lld us matvecs=%lld\n", + (double)est_h, (long long)dur_h, (long long)mv_h); + } else { + printf(" Hutchinson+LFA: SKIPPED\n"); + } + + double err_pp = ref_available ? std::abs((double)est_pp - (double)true_tr) / std::abs((double)true_tr) : -1.0; + double err_h = (!skip_hutchinson && ref_available) + ? std::abs((double)est_h - (double)true_tr) / std::abs((double)true_tr) : -1.0; + if (ref_available) + printf(" Reference: true=%.6e err_pp=%e err_h=%e\n", + (double)true_tr, err_pp, err_h); + else + printf(" Reference: N/A\n"); + + double out_tr = ref_available ? (double)true_tr : -1.0; + std::ofstream file(output_filename, std::ios::out | std::ios::app); + file << n << ", " << k << ", " << data.k_mat << ", " + << s << ", " << d << ", " + << mv_pp << ", " << mv_h << ", " + << dur_pp << ", " << pp_timing.phase1_us << ", " << pp_timing.phase2_us << ", " + << dur_h << ", " + << nt[0] << ", " << nt[1] << ", " << nt[2] << ", " << nt[3] << ", " + << nt[4] << ", " << nt[5] << ", " << nt[6] << ", " << nt[7] << ", " + << nt[8] << ", " << nt[9] << ", " << nt[10] << ", " + << lt[0] << ", " << lt[1] << ", " << lt[2] << ", " << lt[3] << ", " << lt[4] << ", " + << (double)est_pp << ", " << (double)est_h << ", " + << out_tr << ", " << err_pp << ", " << err_h << ",\n"; + } +} + +// ============================================================================ +// call_all_algs_sparse: sparse path (always double). +// ============================================================================ + +template +static void call_all_algs_sparse( + int64_t numruns, + int64_t n, + int64_t k, + int64_t s, + int64_t d, + RandBLAS::sparse_data::CSRMatrix& csr, + const RandBLAS::RNGState& state, + int func_type, + double poly_lambda, + int sketch_type, + int vec_nnz, + bool skip_hutchinson, + const std::string& output_filename +) { + std::function f; + const char* fname; + T f_zero = (T)0; + if (func_type == 0) { + f = [](T x){ return (T)std::sqrt((double)std::max(x, (T)0)); }; + fname = "sqrt"; + } else if (func_type == 1) { + f = [](T x){ return (T)std::log((double)x); }; + fname = "log"; + } else { + T lam = (T)poly_lambda; + f = [lam](T x){ return x * (x + lam); }; + fname = "poly"; + } + + RandLAPACK::SYPS syps(2, 1, false, false); + syps.sketch_type = sketch_type; + syps.vec_nnz = vec_nnz; + ORTH_t orth(false, false); + ORTH_t orth_inner(false, false); + syps.internal_stab = &orth_inner; + RandLAPACK::SYRF, ORTH_t> syrf(syps, orth); + RandLAPACK::NystromEVD nystrom_evd(syrf, 0); + nystrom_evd.timing = true; + LFA_t lfa_pp; + lfa_pp.timing = true; + RandLAPACK::Hutchinson hutch_pp; + RandLAPACK::FunNystromPP + driver(nystrom_evd, lfa_pp, hutch_pp); + + LFA_t lfa_base; + RandLAPACK::Hutchinson hutch_base; + + using SSLO = linops::SparseSymLinOp>; + using CSLO = CountingLinOp; + using FuncT = std::function; + + for (int i = 0; i < numruns; ++i) { + printf("ITERATION %d, N=%lld, K=%lld, K_MAT=0, S=%lld, D=%lld [%s, sparse_mtx]\n", + i, (long long)n, (long long)k, (long long)s, (long long)d, fname); + + auto state_alg = state; + + SSLO A_op_pp(n, csr); + CSLO A_pp(n, A_op_pp); + + RandLAPACK::FunNystromPP_timing pp_timing; + auto start_pp = steady_clock::now(); + T est_pp = driver.call(A_pp, f, f_zero, k, s, d, state_alg, &pp_timing); + auto stop_pp = steady_clock::now(); + long dur_pp = duration_cast(stop_pp - start_pp).count(); + int64_t mv_pp = A_pp.count; + + auto& nt = nystrom_evd.times; + auto& lt = lfa_pp.times; + + printf(" FunNystromPP: est=%.6e time=%lld us (ph1=%lld ph2=%lld) matvecs=%lld\n", + (double)est_pp, (long long)dur_pp, + (long long)pp_timing.phase1_us, (long long)pp_timing.phase2_us, + (long long)mv_pp); + printf(" NystromEVD: alloc=%lld syrf=%lld matvec=%lld gram=%lld potrf=%lld trsm=%lld svd=%lld post_svd=%lld err_est=%lld rest=%lld total=%lld\n", + (long long)nt[0], (long long)nt[1], (long long)nt[2], (long long)nt[3], + (long long)nt[4], (long long)nt[5], (long long)nt[6], (long long)nt[7], + (long long)nt[8], (long long)nt[9], (long long)nt[10]); + printf(" LanczosFA: matvec=%lld lanczos=%lld apply_f=%lld total=%lld\n", + (long long)lt[0], (long long)lt[1], (long long)lt[2], (long long)lt[4]); + + T est_h = (T)0; + long dur_h = 0; + int64_t mv_h = 0; + if (!skip_hutchinson) { + state_alg = state; + SSLO A_op_h(n, csr); + CSLO A_hutch(n, A_op_h); + LFAOp lfa_op(n, A_hutch, lfa_base, f, d); + + auto start_h = steady_clock::now(); + est_h = hutch_base.call(lfa_op, s, state_alg); + auto stop_h = steady_clock::now(); + dur_h = duration_cast(stop_h - start_h).count(); + mv_h = A_hutch.count; + printf(" Hutchinson+LFA: est=%.6e time=%lld us matvecs=%lld\n", + (double)est_h, (long long)dur_h, (long long)mv_h); + } else { + printf(" Hutchinson+LFA: SKIPPED\n"); + } + printf(" Reference: N/A (sparse matrix)\n"); + + std::ofstream file(output_filename, std::ios::out | std::ios::app); + file << n << ", " << k << ", " << 0 << ", " + << s << ", " << d << ", " + << mv_pp << ", " << mv_h << ", " + << dur_pp << ", " << pp_timing.phase1_us << ", " << pp_timing.phase2_us << ", " + << dur_h << ", " + << nt[0] << ", " << nt[1] << ", " << nt[2] << ", " << nt[3] << ", " + << nt[4] << ", " << nt[5] << ", " << nt[6] << ", " << nt[7] << ", " + << nt[8] << ", " << nt[9] << ", " << nt[10] << ", " + << lt[0] << ", " << lt[1] << ", " << lt[2] << ", " << lt[3] << ", " << lt[4] << ", " + << (double)est_pp << ", " << (double)est_h << ", " + << -1.0 << ", " << -1.0 << ", " << -1.0 << ",\n"; + } +} + +// ============================================================================ +// main +// ============================================================================ + +int main(int argc, char* argv[]) { + if (argc < 17) { + std::cerr << "Usage: " << argv[0] + << " " + " " + " " + " [n_1 n_2 ...]\n" + << " dir_path : output directory (use '.' for current dir)\n" + << " num_runs : timed repetitions per matrix size\n" + << " k : Nystrom rank (constant; must satisfy k <= n)\n" + << " k_mat : matrix rank for generated types (for kernel types = data dim)\n" + << " s : Hutchinson sample count (constant)\n" + << " d_steps : Lanczos depth per sample\n" + << " mat_or_file : 1=psd_alg, 2=psd_exp, 3=rbf_kernel, 4=poly_kernel, or path\n" + << " func_type : 0=sqrt, 1=log, 2=poly (x*(x+poly_lambda))\n" + << " compute_ref : 0=skip eigendecomp, 1=run syevd (file input only)\n" + << " sketch_type : 0=Gaussian, 1=SASO\n" + << " vec_nnz : nonzeros per column for SASO (sketch_type=1; ignored otherwise)\n" + << " poly_lambda : lambda for func_type=2 (ignored otherwise)\n" + << " precision : double or float (sparse .mtx path is always double)\n" + << " lfa_type : 0=scalar LanczosFA, 1=BlockLanczosFA\n" + << " orth_type : 0=HouseholderQR (HQRQ), 1=Cholesky QR (CholQRQ)\n" + << " skip_hutchinson : 0=run plain Hutchinson+LFA baseline,\n" + << " 1=skip it (CSV columns padded with zeros).\n" + << " Useful for s/k sweeps that only need the\n" + << " Hutchinson reference at one (lfa, d) variant.\n" + << " n_1 ... : matrix sizes (generated types only)\n"; + return 1; + } + + int64_t numruns = std::stol(argv[2]); + int64_t k_const = std::stol(argv[3]); + int64_t k_mat_const = std::stol(argv[4]); + int64_t s_const = std::stol(argv[5]); + int64_t d_steps = std::stol(argv[6]); + int func_type = std::stoi(argv[8]); + bool compute_ref = std::stoi(argv[9]) != 0; + int sketch_type = std::stoi(argv[10]); + int vec_nnz = std::stoi(argv[11]); + double poly_lambda = std::stod(argv[12]); + std::string precision_str = argv[13]; + bool use_float = (precision_str == "float"); + int lfa_type = std::stoi(argv[14]); + const char* lfa_str = lfa_type == 0 ? "scalar" : "block"; + int orth_type = std::stoi(argv[15]); + const char* orth_str = orth_type == 0 ? "hqrq" : "cholqr"; + bool skip_hutchinson = std::stoi(argv[16]) != 0; + + bool from_file = false; + bool is_mtx = false; + std::string mat_file_path; + int mat_type = 0; + try { + mat_type = std::stoi(argv[7]); + } catch (...) { + from_file = true; + mat_file_path = std::string(argv[7]); + auto& p = mat_file_path; + is_mtx = p.size() >= 4 && p.substr(p.size() - 4) == ".mtx"; + } + + const char* func_str = func_type == 0 ? "sqrt" : (func_type == 1 ? "log" : "poly"); + const char* mat_str = from_file ? mat_file_path.c_str() : + mat_type == 1 ? "psd_alg" : + mat_type == 2 ? "psd_exp" : + mat_type == 3 ? "rbf_kernel" : "poly_kernel"; + const char* sketch_str = sketch_type == 1 ? "SASO" : "Gaussian"; + + std::vector n_sizes; + if (is_mtx) { + int64_t mtx_n = 0; + { + std::ifstream peek(mat_file_path); + std::string ln; + while (std::getline(peek, ln)) + if (!ln.empty() && ln[0] != '%') break; + std::istringstream ss(ln); + int64_t nr, nc, nz; ss >> nr >> nc >> nz; + mtx_n = nr; + } + n_sizes = {mtx_n}; + } else if (from_file) { + int64_t fm = 0, fn = 0; + int wq = 1; + RandLAPACK::gen::process_input_mat( + fm, fn, nullptr, + const_cast(mat_file_path.c_str()), wq); + if (fm != fn) + throw std::runtime_error("File matrix must be square (got " + + std::to_string(fm) + "x" + std::to_string(fn) + ")"); + n_sizes = {fm}; + } else { + if (argc < 18) { + std::cerr << "Error: at least one matrix size (n_1) required for generated types.\n"; + return 1; + } + for (int i = 17; i < argc; ++i) + n_sizes.push_back(std::stol(argv[i])); + } + + std::ostringstream oss; + for (auto v : n_sizes) oss << v << ", "; + std::string sizes_str = oss.str(); + + auto state = RandBLAS::RNGState(); + auto state_constant = state; + + std::string output_filename = + "_FunNystromPP_benchmark_" + precision_str + "_" + std::string(lfa_str) + + "_" + std::string(orth_str) + "_num_info_lines_8.txt"; + std::string path; + if (std::string(argv[1]) != ".") + path = std::string(argv[1]) + output_filename; + else + path = output_filename; + + const std::string mat_construction_str = + is_mtx ? "loaded from .mtx (sparse; both triangles expanded)" : + from_file ? "loaded from .txt (dense)" : + mat_type >= 3 ? "gen_kernel_matrix (random data points x_i~N(0,I_d); cost not timed)" : + "gen_sym_psd_lowrank (Haar-random eigenvectors; cost not timed)"; + + std::ofstream file(path, std::ios::out | std::ios::app); + file << "Description: FunNystromPP vs Hutchinson+LanczosFA benchmark for tr(f(A))." + "\nFile format: 32 columns: n, k, k_mat, s, d, matvec_pp, matvec_hutch," + " time_pp_total (us), time_pp_phase1 (us), time_pp_phase2 (us), time_hutch (us)," + " nystrom_alloc, nystrom_syrf, nystrom_matvec, nystrom_gram," + " nystrom_potrf, nystrom_trsm, nystrom_svd, nystrom_post_svd," + " nystrom_error_est, nystrom_rest, nystrom_total," + " lfa_matvec, lfa_run_lanczos, lfa_apply_f, lfa_rest, lfa_total (all us)," + " est_fun_pp, est_hutch, true_tr (-1 if unavailable), err_fun_pp, err_hutch;" + " rows = numruns repetitions per matrix size." + "\nNum OMP threads: " + std::to_string(RandLAPACK::util::get_omp_threads()) + + "\nInput sizes: " + sizes_str + + "\nParameters: k=" + std::string(argv[3]) + + " k_mat=" + std::string(argv[4]) + + " s=" + std::string(argv[5]) + + " d=" + std::to_string(d_steps) + + " mat=" + std::string(mat_str) + + " func=" + func_str + + (func_type == 2 ? (" poly_lambda=" + std::to_string(poly_lambda)) : std::string("")) + + " sketch=" + sketch_str + + (sketch_type == 1 ? (" vec_nnz=" + std::to_string(vec_nnz)) : std::string("")) + + " precision=" + precision_str + + " lfa_type=" + std::string(lfa_str) + + " orth_type=" + std::string(orth_str) + + " skip_hutch=" + std::to_string((int)skip_hutchinson) + + "\nNum runs per size: " + std::to_string(numruns) + + "\nMatrix construction: " + mat_construction_str + + "\n"; + file.flush(); + + auto start_all = steady_clock::now(); + + using ScalarLFA_d = RandLAPACK::LanczosFA; + using BlockLFA_d = RandLAPACK::BlockLanczosFA; + using ScalarLFA_f = RandLAPACK::LanczosFA; + using BlockLFA_f = RandLAPACK::BlockLanczosFA; + using HQRQ_d = RandLAPACK::HQRQ; + using CholQR_d = RandLAPACK::CholQRQ; + using HQRQ_f = RandLAPACK::HQRQ; + using CholQR_f = RandLAPACK::CholQRQ; + + if (is_mtx) { + // Sparse path: always double. + if (use_float) + printf("Note: sparse .mtx path runs in double regardless of precision flag.\n"); + printf("Loading sparse matrix from %s ...\n", mat_file_path.c_str()); + int64_t n = n_sizes[0]; + auto csr = FunNystromPP_bench::load_mtx(mat_file_path, n); + printf(" Loaded: n=%lld nnz=%lld\n", (long long)n, (long long)csr.nnz); + auto run_sparse = [&]() { + call_all_algs_sparse( + numruns, n, k_const, s_const, d_steps, + csr, state_constant, + func_type, poly_lambda, sketch_type, vec_nnz, + skip_hutchinson, path); + }; + if (lfa_type == 0 && orth_type == 0) run_sparse.template operator()(); + else if (lfa_type == 0 && orth_type == 1) run_sparse.template operator()(); + else if (lfa_type == 1 && orth_type == 0) run_sparse.template operator()(); + else run_sparse.template operator()(); + } else { + // Dense path: dispatch on precision, lfa_type, orth_type (8-way). + auto run_dense = [&]() { + int64_t n_max = *std::max_element(n_sizes.begin(), n_sizes.end()); + int64_t k_mat_max = from_file ? 0 : k_mat_const; + FunNystromPP_benchmark_data all_data( + n_max, k_const, k_mat_max, s_const, d_steps, from_file); + + for (int64_t n : n_sizes) { + all_data.n = n; + all_data.k = k_const; + all_data.k_mat = from_file ? 0 : std::min(k_mat_const, n); + all_data.s = s_const; + all_data.d = d_steps; + + auto state_gen = state_constant; + data_regen(all_data, state_gen, mat_type, func_type, mat_file_path); + + call_all_algs( + numruns, all_data, state_constant, mat_type, + func_type, poly_lambda, compute_ref, sketch_type, vec_nnz, + skip_hutchinson, path); + } + }; + + if (!use_float && lfa_type == 0 && orth_type == 0) run_dense.template operator()(); + else if (!use_float && lfa_type == 0 && orth_type == 1) run_dense.template operator()(); + else if (!use_float && lfa_type == 1 && orth_type == 0) run_dense.template operator()(); + else if (!use_float && lfa_type == 1 && orth_type == 1) run_dense.template operator()(); + else if ( use_float && lfa_type == 0 && orth_type == 0) run_dense.template operator()(); + else if ( use_float && lfa_type == 0 && orth_type == 1) run_dense.template operator()(); + else if ( use_float && lfa_type == 1 && orth_type == 0) run_dense.template operator()(); + else run_dense.template operator()(); + } + + auto stop_all = steady_clock::now(); + long dur_all = duration_cast(stop_all - start_all).count(); + file << "Total benchmark execution time: " + std::to_string(dur_all) + "\n"; + file.flush(); + + return 0; +} +#endif diff --git a/benchmark/bench_FunNystromPP/load_mtx.hh b/benchmark/bench_FunNystromPP/load_mtx.hh new file mode 100644 index 000000000..0493c3c91 --- /dev/null +++ b/benchmark/bench_FunNystromPP/load_mtx.hh @@ -0,0 +1,113 @@ +#pragma once +// Minimal Matrix Market (.mtx) reader. +// Supports: coordinate real/integer, symmetric or general. +// Returns a CSRMatrix with both triangles stored (required by SparseSymLinOp). + +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace FunNystromPP_bench { + +inline RandBLAS::sparse_data::CSRMatrix +load_mtx(const std::string& filename, int64_t& out_n) { + std::ifstream f(filename); + if (!f.is_open()) + throw std::runtime_error("load_mtx: cannot open " + filename); + + // Parse banner line. + std::string banner; + std::getline(f, banner); + std::string blo = banner; + for (char& c : blo) c = (char)std::tolower((unsigned char)c); + if (blo.find("%%matrixmarket") == std::string::npos) + throw std::runtime_error("load_mtx: missing %%MatrixMarket banner"); + if (blo.find("coordinate") == std::string::npos) + throw std::runtime_error("load_mtx: only coordinate format supported"); + if (blo.find("real") == std::string::npos && + blo.find("integer") == std::string::npos && + blo.find("double") == std::string::npos) + throw std::runtime_error("load_mtx: only real/integer/double field type supported"); + bool is_symm = blo.find("symmetric") != std::string::npos; + bool is_skew = blo.find("skew") != std::string::npos; + + // Skip comment lines; find the size line. + std::string line; + while (std::getline(f, line)) { + std::string lo = line; + for (char& c : lo) c = (char)std::tolower((unsigned char)c); + if (lo.find("%%matrixmarket") != std::string::npos) continue; + if (!line.empty() && line[0] != '%') break; + } + + int64_t nrows, ncols, nnz_file; + { + std::istringstream ss(line); + ss >> nrows >> ncols >> nnz_file; + } + if (nrows != ncols) + throw std::runtime_error("load_mtx: matrix must be square (got " + + std::to_string(nrows) + "x" + std::to_string(ncols) + ")"); + out_n = nrows; + + // Read COO entries; expand symmetric/skew-symmetric off-diagonal entries. + std::vector ri, ci; + std::vector vi; + ri.reserve(is_symm ? 2*nnz_file : nnz_file); + ci.reserve(is_symm ? 2*nnz_file : nnz_file); + vi.reserve(is_symm ? 2*nnz_file : nnz_file); + + while (std::getline(f, line)) { + if (line.empty() || line[0] == '%') continue; + std::istringstream ls(line); + int64_t r, c; double v; + if (!(ls >> r >> c >> v)) continue; + --r; --c; // .mtx is 1-indexed; convert to 0-indexed + ri.push_back(r); ci.push_back(c); vi.push_back(v); + if ((is_symm || is_skew) && r != c) { + ri.push_back(c); ci.push_back(r); + vi.push_back(is_skew ? -v : v); + } + } + int64_t total_nnz = (int64_t)vi.size(); + + // Sort COO by (row, col) to build CSR. + std::vector ord(total_nnz); + std::iota(ord.begin(), ord.end(), (int64_t)0); + std::sort(ord.begin(), ord.end(), [&](int64_t a, int64_t b) { + return ri[a] != ri[b] ? ri[a] < ri[b] : ci[a] < ci[b]; + }); + + // Build owning CSRMatrix. + using CSR = RandBLAS::sparse_data::CSRMatrix; + CSR csr(nrows, ncols); + csr.reserve(total_nnz); // allocates vals[total_nnz], colidxs[total_nnz], rowptr[nrows+1]; sets nnz + + // First pass: count entries per row into rowptr[1..nrows]. + std::fill(csr.rowptr, csr.rowptr + nrows + 1, (int64_t)0); + for (int64_t k = 0; k < total_nnz; ++k) + ++csr.rowptr[ri[ord[k]] + 1]; + // Prefix sum → rowptr[0..nrows] = start of each row. + for (int64_t i = 0; i < nrows; ++i) + csr.rowptr[i+1] += csr.rowptr[i]; + + // Second pass: fill colidxs and vals; temporarily advance rowptr[row] to track position. + std::vector pos(csr.rowptr, csr.rowptr + nrows); // copy of row starts + for (int64_t k = 0; k < total_nnz; ++k) { + int64_t row = ri[ord[k]]; + int64_t p = pos[row]++; + csr.colidxs[p] = ci[ord[k]]; + csr.vals[p] = vi[ord[k]]; + } + + return csr; +} + +} // namespace FunNystromPP_bench diff --git a/install.sh b/install.sh index 0c0525173..3ad9a2a26 100755 --- a/install.sh +++ b/install.sh @@ -247,6 +247,7 @@ fi cmake -S $RANDNLA_PROJECT_DIR/lib/blaspp/ -B $RANDNLA_PROJECT_DIR/build/blaspp-build/ \ -Dgpu_backend=$RANDNLA_PROJECT_GPU_AVAIL \ -DCMAKE_BUILD_TYPE=Release \ + -DBUILD_TESTING=OFF \ -Dblas_int=$BLAS_INT \ -DCMAKE_INSTALL_PREFIX=$RANDNLA_PROJECT_DIR/install/blaspp-install/ \ $MACOS_BLAS_FLAGS $MACOS_OPENMP_FLAGS @@ -262,7 +263,7 @@ fi # Configure, build, and install LAPACK++ # Add "-DBLAS_LIBRARIES='-lflame -lblis'" if using AMD AOCL -cmake -S $RANDNLA_PROJECT_DIR/lib/lapackpp/ -B $RANDNLA_PROJECT_DIR/build/lapackpp-build/ -Dgpu_backend=$RANDNLA_PROJECT_GPU_AVAIL -DCMAKE_BUILD_TYPE=Release -Dblaspp_DIR=$RANDNLA_PROJECT_DIR/install/blaspp-install/$LIB_VAR/cmake/blaspp/ -DCMAKE_INSTALL_PREFIX=$RANDNLA_PROJECT_DIR/install/lapackpp-install -DCMAKE_INSTALL_RPATH_USE_LINK_PATH=ON $MACOS_LAPACK_FLAGS $MACOS_OPENMP_FLAGS +cmake -S $RANDNLA_PROJECT_DIR/lib/lapackpp/ -B $RANDNLA_PROJECT_DIR/build/lapackpp-build/ -Dgpu_backend=$RANDNLA_PROJECT_GPU_AVAIL -DCMAKE_BUILD_TYPE=Release -DBUILD_TESTING=OFF -Dblaspp_DIR=$RANDNLA_PROJECT_DIR/install/blaspp-install/$LIB_VAR/cmake/blaspp/ -DCMAKE_INSTALL_PREFIX=$RANDNLA_PROJECT_DIR/install/lapackpp-install -DCMAKE_INSTALL_RPATH_USE_LINK_PATH=ON $MACOS_LAPACK_FLAGS $MACOS_OPENMP_FLAGS make -C $RANDNLA_PROJECT_DIR/build/lapackpp-build/ -j20 install # Configure, build, and install RandLAPACK echo "==========================================" diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 815ae3586..4851b68cb 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -18,8 +18,10 @@ if (GTest_FOUND) drivers/test_cqrrpt.cc drivers/test_cqrrt.cc drivers/test_bqrrp.cc - drivers/test_revd2.cc + drivers/test_nystrom_evd.cc drivers/test_hqrrp.cc + comps/test_lanczos_fa.cc + drivers/test_fun_nystrom_pp.cc drivers/test_abrik.cc drivers/test_orth_linop.cc misc/test_util.cc diff --git a/test/comps/test_lanczos_fa.cc b/test/comps/test_lanczos_fa.cc new file mode 100644 index 000000000..90ea176dd --- /dev/null +++ b/test/comps/test_lanczos_fa.cc @@ -0,0 +1,542 @@ +#include "RandLAPACK.hh" +#include "rl_blaspp.hh" +#include "rl_lapackpp.hh" +#include "rl_gen.hh" + +#include +#include +#include +#include +#include + +namespace linops = RandLAPACK::linops; + +class TestLanczosFA : public ::testing::Test { +protected: + virtual void SetUp() {} + virtual void TearDown() {} + + static std::vector make_diag_matrix(const std::vector& diag) { + int64_t n = diag.size(); + std::vector A(n * n, 0.0); + for (int64_t i = 0; i < n; ++i) + A[i + i * n] = diag[i]; + return A; + } + + // Exact reference for f(A)B when A = diag(d): result[:,j] = f(d) .* B[:,j]. + static std::vector diag_fa_ref( + const std::vector& d, const std::vector& B, + int64_t n, int64_t s, std::function f + ) { + std::vector out(n * s, 0.0); + for (int64_t j = 0; j < s; ++j) + for (int64_t i = 0; i < n; ++i) + out[j * n + i] = f(d[i]) * B[j * n + i]; + return out; + } + + // Exact reference for f(A)B for general symmetric A via dense eigendecomposition: + // A = Q diag(λ) Q' (syevd, reads upper triangle), then f(A)B = Q diag(f(λ)) Q'B. + // A_copy passed by value because syevd overwrites it with eigenvectors. + static std::vector dense_fa_ref( + std::vector A_copy, + const std::vector& B, + int64_t n, int64_t s, + std::function f + ) { + std::vector eigvals(n); + lapack::syevd(lapack::Job::Vec, lapack::Uplo::Upper, n, + A_copy.data(), n, eigvals.data()); + std::vector tmp(n * s); + blas::gemm(Layout::ColMajor, Op::Trans, Op::NoTrans, + n, s, n, 1.0, A_copy.data(), n, B.data(), n, 0.0, tmp.data(), n); + for (int64_t j = 0; j < s; ++j) + for (int64_t i = 0; i < n; ++i) + tmp[j * n + i] *= f(eigvals[i]); + std::vector out(n * s); + blas::gemm(Layout::ColMajor, Op::NoTrans, Op::NoTrans, + n, s, n, 1.0, A_copy.data(), n, tmp.data(), n, 0.0, out.data(), n); + return out; + } + + // A = diag(1,...,n); reference is diag_fa_ref (free, exact). + // With d < n steps the output is approximate; tolerance encodes expected accuracy. + template + static void run_diagonal_fa_test(F f, int64_t n, int64_t s, int64_t d, + double tol, uint64_t seed) { + auto state = RandBLAS::RNGState(seed); + std::vector diag_vec(n); + std::iota(diag_vec.begin(), diag_vec.end(), 1.0); + std::vector A_mat = make_diag_matrix(diag_vec); + std::vector B(n * s); + RandBLAS::DenseDist DB(n, s); + state = RandBLAS::fill_dense(DB, B.data(), state); + std::vector out(n * s, 0.0); + linops::ExplicitSymLinOp A_op(n, blas::Uplo::Upper, A_mat.data(), n, Layout::ColMajor); + RandLAPACK::LanczosFA lfa; + lfa.reorth = 1; + lfa.call(A_op, B.data(), n, s, f, d, out.data()); + auto ref = diag_fa_ref(diag_vec, B, n, s, f); + double err = 0.0, ref_norm = 0.0; + for (int64_t i = 0; i < n * s; ++i) { + double r = out[i] - ref[i]; + err += r * r; + ref_norm += ref[i] * ref[i]; + } + double rel_err = std::sqrt(err / ref_norm); + printf("||f(A)B - ref||/||ref|| = %e\n", rel_err); + ASSERT_LT(rel_err, tol); + } + + // A = B'B + n*I (κ ≈ 5, λ_min ≥ n); reference is dense_fa_ref (syevd). + // Upper triangle only — ExplicitSymLinOp and syevd both read Uplo::Upper. + template + static void run_dense_fa_test(F f, int64_t n, int64_t s, int64_t d, + double tol, uint64_t seed) { + auto state = RandBLAS::RNGState(seed); + std::vector B_raw(n * n); + RandBLAS::DenseDist D1(n, n); + state = RandBLAS::fill_dense(D1, B_raw.data(), state); + std::vector A(n * n, 0.0); + blas::syrk(Layout::ColMajor, Uplo::Upper, Op::Trans, n, n, 1.0, B_raw.data(), n, 0.0, A.data(), n); + for (int64_t i = 0; i < n; ++i) + A[i + i * n] += (double)n; + std::vector B(n * s); + RandBLAS::DenseDist D2(n, s); + state = RandBLAS::fill_dense(D2, B.data(), state); + std::vector out(n * s, 0.0); + linops::ExplicitSymLinOp A_op(n, blas::Uplo::Upper, A.data(), n, Layout::ColMajor); + RandLAPACK::LanczosFA lfa; + lfa.reorth = 1; + lfa.call(A_op, B.data(), n, s, f, d, out.data()); + auto ref = dense_fa_ref(A, B, n, s, f); + double err = 0.0, ref_norm = 0.0; + for (int64_t i = 0; i < n * s; ++i) { + double r = out[i] - ref[i]; + err += r * r; + ref_norm += ref[i] * ref[i]; + } + double rel_err = std::sqrt(err / ref_norm); + printf("||f(A)B - ref||/||ref|| = %e\n", rel_err); + ASSERT_LT(rel_err, tol); + } +}; + + +// f=sqrt on diag(1,...,50), d=20 < n=50. Good but not exact (Krylov subspace smaller than spectrum). +TEST_F(TestLanczosFA, DiagonalSqrt) { + run_diagonal_fa_test([](double x){ return std::sqrt(x); }, 50, 5, 20, 1e-4, 42); +} + +// d=1 on A = 4*I: all eigenvalues equal, so the 1-step Krylov result is exact. +// Exercises the paths where beta is never allocated and the main recurrence loop +// (0..d-2) does not execute. stevd is called on a 1×1 tridiagonal. +TEST_F(TestLanczosFA, ScalarMatrixSqrt_d1) { + int64_t n = 8, s = 3, d = 1; + auto state = RandBLAS::RNGState(55); + + std::vector A_mat(n * n, 0.0); + for (int64_t i = 0; i < n; ++i) A_mat[i + i * n] = 4.0; + + std::vector B(n * s); + RandBLAS::DenseDist DB(n, s); + state = RandBLAS::fill_dense(DB, B.data(), state); + + std::vector out(n * s, 0.0); + linops::ExplicitSymLinOp A_op(n, blas::Uplo::Upper, A_mat.data(), n, Layout::ColMajor); + RandLAPACK::LanczosFA lfa; + lfa.reorth = 1; + lfa.call(A_op, B.data(), n, s, [](double x){ return std::sqrt(x); }, d, out.data()); + + // f(A)B = sqrt(4)*B = 2*B exactly for A = 4*I + double err = 0.0, ref_norm = 0.0; + for (int64_t i = 0; i < n * s; ++i) { + double r = out[i] - 2.0 * B[i]; + err += r * r; + ref_norm += 4.0 * B[i] * B[i]; + } + double rel_err = std::sqrt(err / ref_norm); + printf("d=1 on 4*I: ||f(A)B - ref||/||ref|| = %e\n", rel_err); + ASSERT_LT(rel_err, 1e-10); +} + +// f=log on diag(1,...,50). Entries ≥ 1 so log is well-defined. Tolerance relaxed to 1e-3 +// (log has slower polynomial convergence on [1,50] than sqrt). +TEST_F(TestLanczosFA, DiagonalLog) { + run_diagonal_fa_test([](double x){ return std::log(x); }, 50, 4, 20, 1e-3, 7); +} + +// f=poly x*(x+2) on diag(1,...,50). Degree-2 polynomial; d=20 >> 2 steps needed for exactness. +TEST_F(TestLanczosFA, DiagonalPoly) { + double lam = 2.0; + run_diagonal_fa_test([lam](double x){ return x * (x + lam); }, 50, 5, 20, 1e-4, 17); +} + +// f=sqrt on a random dense PSD matrix. d=n/2=25 steps → near-exact convergence on κ≈5. +TEST_F(TestLanczosFA, DensePSDSqrt) { + run_dense_fa_test([](double x){ return std::sqrt(std::max(x, 0.0)); }, 50, 4, 25, 1e-6, 31); +} + +// f=exp on the same dense PSD setup. Eigenvalues in [50,~250] stay within double range. +// Exercises a function not covered by the diagonal tests (per the reference implementation). +TEST_F(TestLanczosFA, DensePSDExp) { + run_dense_fa_test([](double x){ return std::exp(x); }, 50, 4, 25, 1e-6, 77); +} + +// f=sqrt on an RBF kernel matrix generated by gen_kernel_matrix (n=20, d_dim=5). +// All eigenvalues in (0,1] (Mercer; diagonal=1), so sqrt is well-defined. +// Reference via dense_fa_ref (syevd). Tolerance relaxed to 1e-2: RBF spectrum +// on 5-dim data does not decay fast enough for d=10 to reach tighter accuracy. +TEST_F(TestLanczosFA, KernelRBFSqrt) { + using RNG = r123::Philox4x32; + int64_t n = 20, d_dim = 5, s = 3, d = 10; + auto state = RandBLAS::RNGState(101); + + std::vector K(n * n, 0.0); + double bandwidth = std::sqrt((double)d_dim); + RandLAPACK::gen::gen_kernel_matrix( + n, d_dim, K.data(), n, 0, bandwidth, 0.0, 0, state); + + std::vector B(n * s); + RandBLAS::DenseDist DB(n, s); + state = RandBLAS::fill_dense(DB, B.data(), state); + + std::vector out(n * s, 0.0); + auto f = [](double x){ return std::sqrt(std::max(x, 0.0)); }; + linops::ExplicitSymLinOp A_op(n, blas::Uplo::Upper, K.data(), n, Layout::ColMajor); + RandLAPACK::LanczosFA lfa; + lfa.reorth = 1; + lfa.call(A_op, B.data(), n, s, f, d, out.data()); + + auto ref = dense_fa_ref(K, B, n, s, f); + double err = 0.0, ref_norm = 0.0; + for (int64_t i = 0; i < n * s; ++i) { + double r = out[i] - ref[i]; + err += r * r; + ref_norm += ref[i] * ref[i]; + } + double rel_err = std::sqrt(err / ref_norm); + printf("KernelRBFSqrt: ||f(A)B - ref||/||ref|| = %e\n", rel_err); + ASSERT_LT(rel_err, 1e-2); +} + +// Both reorth modes agree on a well-conditioned diagonal matrix (κ=40). +// Tests that reorthogonalization doesn't introduce a bug, not that it helps +// (for that one would need κ >> 1 where vanilla loses orthogonality). +TEST_F(TestLanczosFA, ReorthVsVanilla) { + int64_t n = 40, s = 3, d = 15; + auto state = RandBLAS::RNGState(99); + + std::vector diag(n); + std::iota(diag.begin(), diag.end(), 1.0); + std::vector A_mat = make_diag_matrix(diag); + + std::vector B(n * s); + RandBLAS::DenseDist DB(n, s); + state = RandBLAS::fill_dense(DB, B.data(), state); + + linops::ExplicitSymLinOp A_op(n, blas::Uplo::Upper, A_mat.data(), n, Layout::ColMajor); + auto f_sqrt = [](double x){ return std::sqrt(x); }; + + using LFA = RandLAPACK::LanczosFA; + std::vector out_full(n * s, 0.0), out_van(n * s, 0.0); + LFA lfa_full, lfa_van; + lfa_full.reorth = 1; + lfa_van.reorth = 0; + + lfa_full.call(A_op, B.data(), n, s, f_sqrt, d, out_full.data()); + lfa_van.call( A_op, B.data(), n, s, f_sqrt, d, out_van.data()); + + double diff = 0.0, norm = 0.0; + for (int64_t i = 0; i < n * s; ++i) { + double r = out_full[i] - out_van[i]; + diff += r * r; + norm += out_full[i] * out_full[i]; + } + double rel_diff = std::sqrt(diff / norm); + printf("Full reorth vs vanilla relative diff: %e\n", rel_diff); + ASSERT_LT(rel_diff, 1e-6); +} + + +// --------------------------------------------------------------------------- +// BlockLanczosFA tests: verify that the block algorithm matches the dense +// reference at sufficient depth and agrees with scalar LanczosFA. + +class TestBlockLanczosFA : public TestLanczosFA { +protected: + template + static void run_block_diagonal_fa_test(F f, int64_t n, int64_t s, int64_t d, + double tol, uint64_t seed) { + auto state = RandBLAS::RNGState(seed); + std::vector diag_vec(n); + std::iota(diag_vec.begin(), diag_vec.end(), 1.0); + std::vector A_mat = make_diag_matrix(diag_vec); + std::vector B(n * s); + RandBLAS::DenseDist DB(n, s); + state = RandBLAS::fill_dense(DB, B.data(), state); + std::vector out(n * s, 0.0); + linops::ExplicitSymLinOp A_op(n, blas::Uplo::Upper, A_mat.data(), n, Layout::ColMajor); + RandLAPACK::BlockLanczosFA blfa; + blfa.reorth = 1; + blfa.call(A_op, B.data(), n, s, f, d, out.data()); + auto ref = diag_fa_ref(diag_vec, B, n, s, f); + double err = 0.0, ref_norm = 0.0; + for (int64_t i = 0; i < n * s; ++i) { + double r = out[i] - ref[i]; + err += r * r; + ref_norm += ref[i] * ref[i]; + } + double rel_err = std::sqrt(err / ref_norm); + printf("Block ||f(A)B - ref||/||ref|| = %e\n", rel_err); + ASSERT_LT(rel_err, tol); + } + + template + static void run_block_dense_fa_test(F f, int64_t n, int64_t s, int64_t d, + double tol, uint64_t seed) { + auto state = RandBLAS::RNGState(seed); + std::vector B_raw(n * n); + RandBLAS::DenseDist D1(n, n); + state = RandBLAS::fill_dense(D1, B_raw.data(), state); + std::vector A(n * n, 0.0); + blas::syrk(Layout::ColMajor, Uplo::Upper, Op::Trans, n, n, 1.0, B_raw.data(), n, 0.0, A.data(), n); + for (int64_t i = 0; i < n; ++i) + A[i + i * n] += (double)n; + std::vector B(n * s); + RandBLAS::DenseDist D2(n, s); + state = RandBLAS::fill_dense(D2, B.data(), state); + std::vector out(n * s, 0.0); + linops::ExplicitSymLinOp A_op(n, blas::Uplo::Upper, A.data(), n, Layout::ColMajor); + RandLAPACK::BlockLanczosFA blfa; + blfa.reorth = 1; + blfa.call(A_op, B.data(), n, s, f, d, out.data()); + auto ref = dense_fa_ref(A, B, n, s, f); + double err = 0.0, ref_norm = 0.0; + for (int64_t i = 0; i < n * s; ++i) { + double r = out[i] - ref[i]; + err += r * r; + ref_norm += ref[i] * ref[i]; + } + double rel_err = std::sqrt(err / ref_norm); + printf("Block dense ||f(A)B - ref||/||ref|| = %e\n", rel_err); + ASSERT_LT(rel_err, tol); + } +}; + +// f=sqrt on diag(1,...,50), d=20 < n. +TEST_F(TestBlockLanczosFA, DiagonalSqrt) { + run_block_diagonal_fa_test([](double x){ return std::sqrt(x); }, 50, 5, 20, 1e-4, 42); +} + +// f=sqrt on a random dense PSD matrix, d=n/2 steps. +TEST_F(TestBlockLanczosFA, DensePSDSqrt) { + run_block_dense_fa_test([](double x){ return std::sqrt(std::max(x, 0.0)); }, 50, 4, 25, 1e-6, 31); +} + +// d=1 edge case: only the initial QR of B is performed; the recurrence loop runs one +// step (computes A_step) but no further blocks are appended. The block Krylov subspace +// is span(Q_0) = span(B), so out = Q_0 * f(Q_0^T A Q_0) * R_0 — a Galerkin approximation +// of f(A)B in the original input subspace. Quality is poor for general A but the path +// must run without crashing and produce a finite, meaningful result. +TEST_F(TestBlockLanczosFA, DStepEqualsOne) { + int64_t n = 30, s = 4, d = 1; + auto state = RandBLAS::RNGState(7); + + // A = diag(1..n); B random. + std::vector diag_vec(n); + std::iota(diag_vec.begin(), diag_vec.end(), 1.0); + std::vector A_mat = make_diag_matrix(diag_vec); + std::vector B(n * s); + RandBLAS::DenseDist DB(n, s); + state = RandBLAS::fill_dense(DB, B.data(), state); + + linops::ExplicitSymLinOp A_op(n, blas::Uplo::Upper, A_mat.data(), n, Layout::ColMajor); + auto f_sqrt = [](double x){ return std::sqrt(std::max(x, 0.0)); }; + + std::vector out(n * s, 0.0); + RandLAPACK::BlockLanczosFA blfa; + blfa.reorth = 1; + blfa.call(A_op, B.data(), n, s, f_sqrt, d, out.data()); + + // Sanity: output is finite, nonzero, and not catastrophically wrong. + double out_norm = 0.0; + bool all_finite = true; + for (int64_t i = 0; i < n * s; ++i) { + out_norm += out[i] * out[i]; + all_finite = all_finite && std::isfinite(out[i]); + } + ASSERT_TRUE(all_finite); + ASSERT_GT(std::sqrt(out_norm), 1e-3); + printf("d=1 block ||out|| = %e (sanity check, accuracy not asserted)\n", std::sqrt(out_norm)); +} + +// reorth=0 (vanilla block recurrence) vs reorth=1 (full reorthogonalization) on the +// same problem. With moderate d on a well-conditioned matrix, both should converge +// to a similar answer; this exercises the no-reorth code path and guards against +// regressions in it. +TEST_F(TestBlockLanczosFA, NoReorthAgreesWithFull) { + int64_t n = 50, s = 4, d = 25; + auto state = RandBLAS::RNGState(31); + + // PSD A = G^T G + n*I (well-conditioned). + std::vector B_raw(n * n); + RandBLAS::DenseDist D1(n, n); + state = RandBLAS::fill_dense(D1, B_raw.data(), state); + std::vector A(n * n, 0.0); + blas::syrk(Layout::ColMajor, Uplo::Upper, Op::Trans, n, n, 1.0, B_raw.data(), n, 0.0, A.data(), n); + for (int64_t i = 0; i < n; ++i) A[i + i * n] += (double)n; + + std::vector B(n * s); + RandBLAS::DenseDist D2(n, s); + state = RandBLAS::fill_dense(D2, B.data(), state); + + linops::ExplicitSymLinOp A_op(n, blas::Uplo::Upper, A.data(), n, Layout::ColMajor); + auto f_sqrt = [](double x){ return std::sqrt(std::max(x, 0.0)); }; + + std::vector out_full(n * s, 0.0), out_van(n * s, 0.0); + RandLAPACK::BlockLanczosFA blfa_full, blfa_van; + blfa_full.reorth = 1; + blfa_van.reorth = 0; + + blfa_full.call(A_op, B.data(), n, s, f_sqrt, d, out_full.data()); + blfa_van.call (A_op, B.data(), n, s, f_sqrt, d, out_van.data()); + + double diff = 0.0, norm = 0.0; + for (int64_t i = 0; i < n * s; ++i) { + double r = out_full[i] - out_van[i]; + diff += r * r; + norm += out_full[i] * out_full[i]; + } + double rel_diff = std::sqrt(diff / norm); + printf("Block reorth=1 vs reorth=0 relative diff: %e\n", rel_diff); + ASSERT_LT(rel_diff, 1e-6); +} + +// Zero-column input: if B has a column of zeros, the corresponding output column +// must be zero. The initial QR produces R_0 with a zero column at that index, which +// propagates to a zero output column via the formula out = Q_basis * f(T_k) * E_1 * R_0. +TEST_F(TestBlockLanczosFA, ZeroColumn) { + int64_t n = 40, s = 3, d = 10; + auto state = RandBLAS::RNGState(11); + + std::vector diag_vec(n); + std::iota(diag_vec.begin(), diag_vec.end(), 1.0); + std::vector A_mat = make_diag_matrix(diag_vec); + std::vector B(n * s); + RandBLAS::DenseDist DB(n, s); + state = RandBLAS::fill_dense(DB, B.data(), state); + // Zero out the middle column. + for (int64_t i = 0; i < n; ++i) B[i + 1 * n] = 0.0; + + linops::ExplicitSymLinOp A_op(n, blas::Uplo::Upper, A_mat.data(), n, Layout::ColMajor); + auto f_sqrt = [](double x){ return std::sqrt(std::max(x, 0.0)); }; + + std::vector out(n * s, 0.0); + RandLAPACK::BlockLanczosFA blfa; + blfa.reorth = 1; + blfa.call(A_op, B.data(), n, s, f_sqrt, d, out.data()); + + double zero_col_norm = 0.0, other_cols_norm = 0.0; + for (int64_t i = 0; i < n; ++i) zero_col_norm += out[i + 1 * n] * out[i + 1 * n]; + for (int64_t j = 0; j < s; ++j) { + if (j == 1) continue; + for (int64_t i = 0; i < n; ++i) + other_cols_norm += out[i + j * n] * out[i + j * n]; + } + zero_col_norm = std::sqrt(zero_col_norm); + other_cols_norm = std::sqrt(other_cols_norm); + printf("Zero-column out norm: %e (others: %e)\n", zero_col_norm, other_cols_norm); + // Output column for zero input column should be at roundoff level. + ASSERT_LT(zero_col_norm, 1e-10); + // Other columns should be substantially nonzero (sanity). + ASSERT_GT(other_cols_norm, 1e-3); +} + +// BlockLanczosFA vs LanczosFA on the same problem — outputs should agree at sufficient depth. +TEST_F(TestBlockLanczosFA, AgreeWithScalar) { + using RNG = r123::Philox4x32; + int64_t n = 50, s = 4, d = 25; + auto state = RandBLAS::RNGState(31); + + std::vector B_raw(n * n); + RandBLAS::DenseDist D1(n, n); + state = RandBLAS::fill_dense(D1, B_raw.data(), state); + std::vector A(n * n, 0.0); + blas::syrk(Layout::ColMajor, Uplo::Upper, Op::Trans, n, n, 1.0, B_raw.data(), n, 0.0, A.data(), n); + for (int64_t i = 0; i < n; ++i) + A[i + i * n] += (double)n; + + std::vector B(n * s); + RandBLAS::DenseDist D2(n, s); + state = RandBLAS::fill_dense(D2, B.data(), state); + + linops::ExplicitSymLinOp A_op(n, blas::Uplo::Upper, A.data(), n, Layout::ColMajor); + auto f_sqrt = [](double x){ return std::sqrt(std::max(x, 0.0)); }; + + std::vector out_scalar(n * s, 0.0), out_block(n * s, 0.0); + RandLAPACK::LanczosFA lfa; + lfa.reorth = 1; + lfa.call(A_op, B.data(), n, s, f_sqrt, d, out_scalar.data()); + + RandLAPACK::BlockLanczosFA blfa; + blfa.reorth = 1; + blfa.call(A_op, B.data(), n, s, f_sqrt, d, out_block.data()); + + double diff = 0.0, norm = 0.0; + for (int64_t i = 0; i < n * s; ++i) { + double r = out_block[i] - out_scalar[i]; + diff += r * r; + norm += out_scalar[i] * out_scalar[i]; + } + double rel_diff = std::sqrt(diff / norm); + printf("Block vs scalar relative diff: %e\n", rel_diff); + ASSERT_LT(rel_diff, 1e-4); +} + +// --------------------------------------------------------------------------- +// Hutchinson tests live here because Hutchinson is a building block of +// FunNystromPP alongside LanczosFA, not a standalone algorithm with its own file. + +class TestHutchinson : public ::testing::Test { +protected: + virtual void SetUp() {} + virtual void TearDown() {} +}; + +// M = diag(1,...,100): tr(M) = n*(n+1)/2 = 5050 exactly. +// s=500 Rademacher samples gives std/mean ≈ 0.5%; 5% tolerance is conservative. +// Uses a custom DiagonalOp to test the SymmetricLinearOperator interface independently +// of ExplicitSymLinOp. +TEST_F(TestHutchinson, DiagonalTrace) { + using RNG = r123::Philox4x32; + int64_t n = 100, s = 500; + auto state = RandBLAS::RNGState(13); + + std::vector d(n); + std::iota(d.begin(), d.end(), 1.0); + double true_trace = n * (n + 1.0) / 2.0; + + struct DiagonalOp { + using scalar_t = double; + const int64_t dim; + const std::vector& diag; + void operator()([[maybe_unused]] Layout layout, int64_t n_vecs, double alpha, + double* const B, int64_t ldb, double beta, double* C, int64_t ldc) { + for (int64_t j = 0; j < n_vecs; ++j) + for (int64_t i = 0; i < dim; ++i) { + double val = alpha * diag[i] * B[j * ldb + i]; + C[j * ldc + i] = (beta == 0.0) ? val : val + beta * C[j * ldc + i]; + } + } + }; + + DiagonalOp M{n, d}; + RandLAPACK::Hutchinson hutch; + double est = hutch.call(M, s, state); + + double rel_err = std::abs(est - true_trace) / true_trace; + printf("Hutchinson trace estimate: est=%e, true=%e, rel_err=%e\n", est, true_trace, rel_err); + ASSERT_LT(rel_err, 0.05); +} diff --git a/test/comps/test_syrf.cc b/test/comps/test_syrf.cc index b8afc19b7..638e640ab 100644 --- a/test/comps/test_syrf.cc +++ b/test/comps/test_syrf.cc @@ -99,7 +99,7 @@ class TestSYRF : public ::testing::Test auto m = all_data.row; auto k = all_data.rank; - all_algs.syrf.call(Uplo::Upper, m, all_data.A.data(), k, all_data.Q, state, NULL); + all_algs.syrf.call(Uplo::Upper, m, all_data.A.data(), k, all_data.Q.data(), state, NULL); // Reassing pointers because Q, B have been resized T* Q_dat = all_data.Q.data(); diff --git a/test/drivers/test_bqrrp.cc b/test/drivers/test_bqrrp.cc index f58b7ad43..db9e1d117 100644 --- a/test/drivers/test_bqrrp.cc +++ b/test/drivers/test_bqrrp.cc @@ -67,7 +67,7 @@ class TestBQRRP : public ::testing::Test auto n = all_data.col; auto k = all_data.rank; - RandLAPACK::util::upsize(k * k, all_data.I_ref); + RandLAPACK::util::resize(k * k, all_data.I_ref); RandLAPACK::util::eye(k, k, all_data.I_ref); T* A_dat = all_data.A_cpy1.data(); @@ -131,7 +131,7 @@ class TestBQRRP : public ::testing::Test all_data.rank = BQRRP.rank; printf("RANK AS RETURNED BY BQRRP %4ld\n", all_data.rank); - RandLAPACK::util::upsize(all_data.rank * n, all_data.R); + RandLAPACK::util::resize(all_data.rank * n, all_data.R); lapack::lacpy(MatrixType::Upper, all_data.rank, n, all_data.A.data(), m, all_data.R.data(), all_data.rank); diff --git a/test/drivers/test_bqrrp_gpu.cu b/test/drivers/test_bqrrp_gpu.cu index 62f187d9b..29ed87e41 100644 --- a/test/drivers/test_bqrrp_gpu.cu +++ b/test/drivers/test_bqrrp_gpu.cu @@ -118,7 +118,7 @@ class TestBQRRP : public ::testing::TestWithParam auto n = all_data.col; auto k = all_data.rank; - RandLAPACK::util::upsize(k * k, all_data.I_ref); + RandLAPACK::util::resize(k * k, all_data.I_ref); RandLAPACK::util::eye(k, k, all_data.I_ref); T* A_dat = all_data.A_cpy1.data(); @@ -188,7 +188,7 @@ class TestBQRRP : public ::testing::TestWithParam cudaMemcpy(all_data.J.data(), all_data.J_device, n * sizeof(int64_t), cudaMemcpyDeviceToHost); lapack::ungqr(m, n, n, all_data.Q.data(), m, all_data.tau.data()); - RandLAPACK::util::upsize(all_data.rank * n, all_data.R); + RandLAPACK::util::resize(all_data.rank * n, all_data.R); lapack::lacpy(MatrixType::Upper, all_data.rank, n, all_data.R_full.data(), m, all_data.R.data(), all_data.rank); RandLAPACK::util::col_swap(m, n, n, all_data.A_cpy1.data(), m, all_data.J); diff --git a/test/drivers/test_cqrrpt.cc b/test/drivers/test_cqrrpt.cc index 6a07ee3c5..1bc32e984 100644 --- a/test/drivers/test_cqrrpt.cc +++ b/test/drivers/test_cqrrpt.cc @@ -63,7 +63,7 @@ class TestCQRRPT : public ::testing::Test auto n = all_data.col; auto k = all_data.rank; - RandLAPACK::util::upsize(k * k, all_data.I_ref); + RandLAPACK::util::resize(k * k, all_data.I_ref); RandLAPACK::util::eye(k, k, all_data.I_ref); T* A_dat = all_data.A_cpy1.data(); diff --git a/test/drivers/test_cqrrpt_gpu.cu b/test/drivers/test_cqrrpt_gpu.cu index fd51e6b2e..4dd541586 100644 --- a/test/drivers/test_cqrrpt_gpu.cu +++ b/test/drivers/test_cqrrpt_gpu.cu @@ -69,7 +69,7 @@ class TestCQRRPT : public ::testing::Test auto n = all_data.col; auto k = all_data.rank; - RandLAPACK::util::upsize(k * k, all_data.I_ref); + RandLAPACK::util::resize(k * k, all_data.I_ref); RandLAPACK::util::eye(k, k, all_data.I_ref); T* A_dat = all_data.A_cpy1.data(); diff --git a/test/drivers/test_cqrrt.cc b/test/drivers/test_cqrrt.cc index fe8e65876..d0e682773 100644 --- a/test/drivers/test_cqrrt.cc +++ b/test/drivers/test_cqrrt.cc @@ -59,7 +59,7 @@ class TestCQRRT : public ::testing::Test auto n = all_data.col; auto k = n; - RandLAPACK::util::upsize(k * k, all_data.I_ref); + RandLAPACK::util::resize(k * k, all_data.I_ref); RandLAPACK::util::eye(k, k, all_data.I_ref); T* A_dat = all_data.A_cpy1.data(); diff --git a/test/drivers/test_fun_nystrom_pp.cc b/test/drivers/test_fun_nystrom_pp.cc new file mode 100644 index 000000000..15f3df6ce --- /dev/null +++ b/test/drivers/test_fun_nystrom_pp.cc @@ -0,0 +1,247 @@ +#include "RandLAPACK.hh" +#include "rl_blaspp.hh" +#include "rl_lapackpp.hh" +#include "rl_gen.hh" + +#include +#include +#include +#include +#include + +namespace linops = RandLAPACK::linops; + +class TestFunNystromPP : public ::testing::Test { +protected: + virtual void SetUp() {} + virtual void TearDown() {} + + // Exact reference for tr(f(A)) via dense syevd: sum of f(λ_i). + // A_copy passed by value because syevd overwrites it with eigenvectors. + static double true_trace_fa( + int64_t n, std::vector A_copy, + std::function f + ) { + std::vector eigvals(n); + lapack::syevd(lapack::Job::NoVec, lapack::Uplo::Upper, n, + A_copy.data(), n, eigvals.data()); + double tr = 0.0; + for (int64_t i = 0; i < n; ++i) + tr += f(eigvals[i]); + return tr; + } + + // Full algorithm stack. error_est_power_iters=0 → fixed-rank single pass in NystromEVD. + template + struct Algs { + using SYPS_t = RandLAPACK::SYPS; + using Orth_t = RandLAPACK::HQRQ; + using SYRF_t = RandLAPACK::SYRF; + using NystromEVD_t = RandLAPACK::NystromEVD; + using LFA_t = RandLAPACK::LanczosFA; + using Hutch_t = RandLAPACK::Hutchinson; + using Driver_t = RandLAPACK::FunNystromPP; + + SYPS_t syps; + Orth_t orth; + SYRF_t syrf; + NystromEVD_t nystrom_evd; + LFA_t lfa; + Hutch_t hutch; + Driver_t driver; + + Algs() : + syps(3, 1, false, false), + orth(false, false), + syrf(syps, orth), + nystrom_evd(syrf, 0), + driver(nystrom_evd, lfa, hutch) {} + }; + + // Common test body: run the full pipeline, print and assert relative error. + template + static void test_FunNystromPP_general( + Algs& algs, + linops::ExplicitSymLinOp& A_op, + F f, double f_zero, double true_tr, + int64_t k, int64_t s, int64_t d, double tol, + RandBLAS::RNGState& state + ) { + double est = algs.driver.call(A_op, f, f_zero, k, s, d, state); + double rel_err = std::abs(est - true_tr) / true_tr; + printf("FunNystromPP: est=%e, true=%e, rel_err=%e\n", est, true_tr, rel_err); + ASSERT_LT(rel_err, tol); + } +}; + + +// Random dense PSD matrix A = B'B + n*I, f=sqrt. Reference via syevd. +// k=10 << n=30: Nyström is coarse, Hutchinson correction carries real weight. Tol = 10%. +TEST_F(TestFunNystromPP, DensePSDSqrt) { + using RNG = r123::Philox4x32; + int64_t n = 30, k = 10, s = 200, d = 15; + auto state = RandBLAS::RNGState(0); + + std::vector B_raw(n * n); + RandBLAS::DenseDist DB(n, n); + state = RandBLAS::fill_dense(DB, B_raw.data(), state); + + std::vector A(n * n, 0.0); + blas::syrk(Layout::ColMajor, Uplo::Upper, Op::Trans, n, n, 1.0, B_raw.data(), n, 0.0, A.data(), n); + for (int64_t i = 0; i < n; ++i) + A[i + i * n] += (double)n; + // Symmetrize so both syevd (Uplo::Upper) and true_trace_fa see a consistent matrix. + for (int64_t j = 0; j < n; ++j) + for (int64_t i = j + 1; i < n; ++i) + A[i + j * n] = A[j + i * n]; + + auto f_sqrt = [](double x){ return std::sqrt(std::max(x, 0.0)); }; + double true_tr = true_trace_fa(n, A, f_sqrt); + + linops::ExplicitSymLinOp A_op(n, blas::Uplo::Upper, A.data(), n, Layout::ColMajor); + Algs algs; + test_FunNystromPP_general(algs, A_op, f_sqrt, 0.0, true_tr, k, s, d, 0.1, state); +} + + +// Diagonal A = diag(1,...,n), f=sqrt. True tr(sqrt(A)) = Σ sqrt(i) requires no eigensolver. +// k=15 of n=50 captures the dominant eigenvalues; tail is (n-k)*f(0)=0. Tol = 5%. +TEST_F(TestFunNystromPP, DiagonalSqrt) { + using RNG = r123::Philox4x32; + int64_t n = 50, k = 15, s = 300, d = 20; + auto state = RandBLAS::RNGState(1); + + std::vector A(n * n, 0.0); + double true_tr = 0.0; + for (int64_t i = 0; i < n; ++i) { + A[i + i * n] = (double)(i + 1); + true_tr += std::sqrt((double)(i + 1)); + } + + auto f_sqrt = [](double x){ return std::sqrt(x); }; + linops::ExplicitSymLinOp A_op(n, blas::Uplo::Upper, A.data(), n, Layout::ColMajor); + Algs algs; + test_FunNystromPP_general(algs, A_op, f_sqrt, 0.0, true_tr, k, s, d, 0.05, state); +} + + +// A = gen_shifted_lowrank_psd(noise=1): eigenvalues {lambda_j+1} ∪ {1 repeated n-k_mat times}. +// f=log; all eigenvalues ≥ 1 so log is well-defined. f(noise) = log(1) = 0 → f_zero=0. +// true_tr = Σ log(lambda_j + 1); (n-k_mat) tail terms vanish. Tol = 20%. +TEST_F(TestFunNystromPP, ShiftedLowRankLog) { + using RNG = r123::Philox4x32; + int64_t n = 40, k_mat = 20, k = 15, s = 500, d = 20; + double noise = 1.0; + auto state = RandBLAS::RNGState(2); + + std::vector eigvals(k_mat); + RandLAPACK::gen::gen_alg_decay_singvals(k_mat, 100.0, 2.0, eigvals.data()); + + std::vector A(n * n, 0.0); + RandLAPACK::gen::gen_shifted_lowrank_psd(n, k_mat, A.data(), n, eigvals.data(), noise, state); + + auto f_log = [](double x){ return std::log(x); }; + double true_tr = 0.0; + for (int64_t j = 0; j < k_mat; ++j) + true_tr += f_log(eigvals[j] + noise); + // tail: (n - k_mat) * log(1) = 0 + + linops::ExplicitSymLinOp A_op(n, blas::Uplo::Upper, A.data(), n, Layout::ColMajor); + Algs algs; + test_FunNystromPP_general(algs, A_op, f_log, 0.0, true_tr, k, s, d, 0.20, state); +} + + +// A = gen_sym_psd_lowrank; f=x*(x+1), degree-2 polynomial. d=2 Lanczos steps is exact for +// degree-2 trace estimation via Gauss quadrature (exact for degree ≤ 2d-1 = 3 quadratic forms). +// true_tr = Σ lambda_j*(lambda_j+1); (n-k_mat) tail terms: 0*(0+1)=0. Tol = 10%. +TEST_F(TestFunNystromPP, LowRankPoly) { + using RNG = r123::Philox4x32; + int64_t n = 40, k_mat = 20, k = 15, s = 200, d = 2; + double lam = 1.0; + auto state = RandBLAS::RNGState(3); + + std::vector eigvals(k_mat); + RandLAPACK::gen::gen_alg_decay_singvals(k_mat, 100.0, 2.0, eigvals.data()); + + std::vector A(n * n, 0.0); + RandLAPACK::gen::gen_sym_psd_lowrank(n, k_mat, A.data(), n, eigvals.data(), state); + + auto f_poly = [lam](double x){ return x * (x + lam); }; + double true_tr = 0.0; + for (int64_t j = 0; j < k_mat; ++j) + true_tr += f_poly(eigvals[j]); + // tail: (n - k_mat) * f(0) = 0 + + linops::ExplicitSymLinOp A_op(n, blas::Uplo::Upper, A.data(), n, Layout::ColMajor); + Algs algs; + test_FunNystromPP_general(algs, A_op, f_poly, 0.0, true_tr, k, s, d, 0.10, state); +} + + +// A = RBF kernel matrix via gen_kernel_matrix (n=30, d_dim=5). Eigenvalues in (0,1]. +// f=sqrt; reference via syevd. k=10, s=300, d=15. Tol = 15%. +TEST_F(TestFunNystromPP, KernelRBFSqrt) { + using RNG = r123::Philox4x32; + int64_t n = 30, d_dim = 5, k = 10, s = 300, d = 15; + auto state = RandBLAS::RNGState(4); + + std::vector K(n * n, 0.0); + double bandwidth = std::sqrt((double)d_dim); + RandLAPACK::gen::gen_kernel_matrix( + n, d_dim, K.data(), n, 0, bandwidth, 0.0, 0, state); + + auto f_sqrt = [](double x){ return std::sqrt(std::max(x, 0.0)); }; + double true_tr = true_trace_fa(n, K, f_sqrt); + + linops::ExplicitSymLinOp A_op(n, blas::Uplo::Upper, K.data(), n, Layout::ColMajor); + Algs algs; + test_FunNystromPP_general(algs, A_op, f_sqrt, 0.0, true_tr, k, s, d, 0.15, state); +} + + +// At k=k_mat on a rank-k_mat matrix, NystromEVD captures the entire spectrum, so Phase 1 +// computes tr(f(A)) to machine precision and Phase 2's true correction is zero. +// However, Phase 2 (Hutchinson+LanczosFA) produces a small nonzero t2 because LanczosFA +// and the GEMM-based f(Â)*Omega in ResidualOp follow different floating-point arithmetic +// paths, disagreeing at ~1e-14 per column. Hutchinson accumulates this into a systematic +// bias of roughly 1e-5 relative. This test documents that expected behavior: +// - Phase 1 (t1) is machine-precision accurate: err < 1e-10 +// - Total error is bounded by the Phase 2 floor: err < 1e-3 +// - t2 is nonzero but small (the two-path arithmetic residual) +// To hit machine precision on the total, skip Phase 2 when k >= effective rank. +TEST_F(TestFunNystromPP, Phase1MachinePrecisionAtKmat) { + using RNG = r123::Philox4x32; + int64_t n = 200, k_mat = 10, k = 10, s = 200, d = 20; + auto state = RandBLAS::RNGState(7); + + std::vector eigvals(k_mat); + RandLAPACK::gen::gen_alg_decay_singvals(k_mat, 100.0, 2.0, eigvals.data()); + + std::vector A(n * n, 0.0); + RandLAPACK::gen::gen_sym_psd_lowrank(n, k_mat, A.data(), n, eigvals.data(), state); + + auto f = [](double x){ return std::sqrt(std::max(x, 0.0)); }; + double f_zero = 0.0; + double true_tr = 0.0; + for (int64_t j = 0; j < k_mat; ++j) + true_tr += f(eigvals[j]); + + linops::ExplicitSymLinOp A_op(n, blas::Uplo::Upper, A.data(), n, Layout::ColMajor); + Algs algs; + double est = algs.driver.call(A_op, f, f_zero, k, s, d, state); + + double t1 = 0.0; + for (int64_t i = 0; i < k; ++i) + t1 += algs.driver.F_vec[i]; + double t2 = est - t1; + + double err_t1 = std::abs(t1 - true_tr) / true_tr; + double err_tot = std::abs(est - true_tr) / true_tr; + + printf("Phase1MachinePrecision: t1=%e t2=%e err_t1=%e err_tot=%e\n", + t1, t2, err_t1, err_tot); + + ASSERT_LT(err_t1, 1e-10); // Phase 1 machine precision at k=k_mat + ASSERT_LT(err_tot, 1e-3); // total bounded by Phase 2 LFA floor +} diff --git a/test/drivers/test_hqrrp.cc b/test/drivers/test_hqrrp.cc index b2a546bc6..419ffe288 100644 --- a/test/drivers/test_hqrrp.cc +++ b/test/drivers/test_hqrrp.cc @@ -66,7 +66,7 @@ class TestHQRRP : public ::testing::Test auto n = all_data.col; auto k = all_data.rank; - RandLAPACK::util::upsize(k * k, all_data.I_ref); + RandLAPACK::util::resize(k * k, all_data.I_ref); RandLAPACK::util::eye(k, k, all_data.I_ref); T* A_dat = all_data.A_cpy1.data(); @@ -124,7 +124,7 @@ class TestHQRRP : public ::testing::Test RandLAPACK::hqrrp(m, n, all_data.A.data(), m, all_data.J.data(), all_data.tau.data(), b_sz, (int64_t) (d_factor * b_sz), panel_pivoting, use_cholqr, state, (T**) nullptr); - RandLAPACK::util::upsize(all_data.rank * n, all_data.R); + RandLAPACK::util::resize(all_data.rank * n, all_data.R); lapack::lacpy(MatrixType::Upper, all_data.rank, n, all_data.A.data(), m, all_data.R.data(), all_data.rank); lapack::ungqr(m, n, n, all_data.A.data(), m, all_data.tau.data()); diff --git a/test/drivers/test_revd2.cc b/test/drivers/test_nystrom_evd.cc similarity index 73% rename from test/drivers/test_revd2.cc rename to test/drivers/test_nystrom_evd.cc index 6cfad98bb..269c87459 100644 --- a/test/drivers/test_revd2.cc +++ b/test/drivers/test_nystrom_evd.cc @@ -8,7 +8,7 @@ #include -class TestREVD2 : public ::testing::Test +class TestNystromEVD : public ::testing::Test { protected: @@ -17,61 +17,62 @@ class TestREVD2 : public ::testing::Test virtual void TearDown() {}; template - struct REVD2TestData { + struct NystromEVDTestData { int64_t dim; int64_t rank; std::vector A; std::vector A_cpy; - std::vector V; - std::vector eigvals; + T* V = nullptr; int64_t V_sz = 0; + T* eigvals = nullptr; int64_t eigvals_sz = 0; std::vector E; std::vector Buf; - REVD2TestData( + NystromEVDTestData( int64_t m, int64_t k - ) : - A(m * m, 0.0), - A_cpy(m * m, 0.0), - V(m * k, 0.0), - eigvals(k, 0.0), - E(k * k, 0.0), + ) : + A(m * m, 0.0), + A_cpy(m * m, 0.0), + E(k * k, 0.0), Buf(m * k, 0.0) { dim = m; rank = k; } + + ~NystromEVDTestData() { delete[] V; delete[] eigvals; } }; - template - struct REVD2UploTestData { + template + struct NystromEVDUploTestData { int64_t dim; int64_t rank; std::vector work; std::vector A_u; - std::vector V_u; - std::vector eigvals_u; + T* V_u = nullptr; int64_t V_u_sz = 0; + T* eigvals_u = nullptr; int64_t eigvals_u_sz = 0; std::vector A_l; - std::vector V_l; - std::vector eigvals_l; + T* V_l = nullptr; int64_t V_l_sz = 0; + T* eigvals_l = nullptr; int64_t eigvals_l_sz = 0; std::vector E_u; std::vector E_l; - REVD2UploTestData( + NystromEVDUploTestData( int64_t m, int64_t k - ) : + ) : work(m * m, 0.0), A_u(m * m, 0.0), - V_u(m * k, 0.0), - eigvals_u(k, 0.0), A_l(m * m, 0.0), - V_l(m * k, 0.0), - eigvals_l(k, 0.0), - E_u(k * k, 0.0), + E_u(k * k, 0.0), E_l(k * k, 0.0) { dim = m; rank = k; } + + ~NystromEVDUploTestData() { + delete[] V_u; delete[] eigvals_u; + delete[] V_l; delete[] eigvals_l; + } }; template @@ -82,7 +83,7 @@ class TestREVD2 : public ::testing::Test SYPS_t syps; Orth_t orth; SYRF_t syrf; - RandLAPACK::REVD2 revd2; + RandLAPACK::NystromEVD nystrom_evd; algorithm_objects( @@ -95,12 +96,12 @@ class TestREVD2 : public ::testing::Test syps(num_syps_passes, passes_per_syps_stabilization, verbose, cond_check), orth(cond_check, verbose), syrf(syps, orth, verbose, cond_check), - revd2(syrf, num_steps_power_iter_error_est, verbose) + nystrom_evd(syrf, num_steps_power_iter_error_est, verbose) {} }; template - static void symm_mat_and_copy_computational_helper(T &norm_A, REVD2TestData &all_data) { + static void symm_mat_and_copy_computational_helper(T &norm_A, NystromEVDTestData &all_data) { auto m = all_data.dim; // We're using Nystrom, the original must be positive semidefinite blas::syrk( @@ -116,7 +117,7 @@ class TestREVD2 : public ::testing::Test } template - static void uplo_computational_helper(REVD2UploTestData &all_data) { + static void uplo_computational_helper(NystromEVDUploTestData &all_data) { auto m = all_data.dim; T* A_u_dat = all_data.A_u.data(); T* A_l_dat = all_data.A_l.data(); @@ -138,13 +139,13 @@ class TestREVD2 : public ::testing::Test /// General test for REVD: /// Computes the decomposition factors, then checks A-U\Sigma\transpose{V}. template - static void test_REVD2_general( + static void test_NystromEVD_general( int64_t k_start, T tol, int rank_expectation, T err_expectation, T &norm_A, - REVD2TestData &all_data, + NystromEVDTestData &all_data, algorithm_objects &all_algs, RandBLAS::RNGState state ) { @@ -152,18 +153,19 @@ class TestREVD2 : public ::testing::Test auto m = all_data.dim; int64_t k = k_start; - all_algs.revd2.call(blas::Uplo::Upper, m, all_data.A.data(), k, tol, all_data.V, all_data.eigvals, state); + all_algs.nystrom_evd.call(blas::Uplo::Upper, m, all_data.A.data(), k, tol, + all_data.V, all_data.V_sz, all_data.eigvals, all_data.eigvals_sz, state); - T* E_dat = RandLAPACK::util::upsize(k * k, all_data.E); - T* Buf_dat = RandLAPACK::util::upsize(m * k, all_data.Buf); + T* E_dat = RandLAPACK::util::resize(k * k, all_data.E); + T* Buf_dat = RandLAPACK::util::resize(m * k, all_data.Buf); T* A_cpy_dat = all_data.A_cpy.data(); - T* V_dat = all_data.V.data(); + T* V_dat = all_data.V; - // Construnct A_hat = U1 * S1 * VT1 + // Construct A_hat = V * diag(eigvals) * V' - // Turn vector into diagonal matrix - RandLAPACK::util::diag(k, k, all_data.eigvals.data(), k, all_data.E.data()); + // Turn array into diagonal matrix + RandLAPACK::util::diag(k, k, all_data.eigvals, k, all_data.E.data()); // V * E = Buf blas::gemm(Layout::ColMajor, Op::NoTrans, Op::NoTrans, m, k, k, 1.0, V_dat, m, E_dat, k, 0.0, Buf_dat, m); // A - Buf * V' - should be close to 0 @@ -171,18 +173,18 @@ class TestREVD2 : public ::testing::Test T norm_0 = lapack::lange(Norm::Fro, m, m, A_cpy_dat, m); printf("||A - VEV'||_F / ||A||_F: %e\n", norm_0 / norm_A); - ASSERT_NEAR(norm_0 / norm_A, err_expectation, 10 * err_expectation); - ASSERT_NEAR(k, rank_expectation, std::numeric_limits::epsilon()); + ASSERT_LT(norm_0 / norm_A, 10 * err_expectation); + ASSERT_EQ(k, rank_expectation); } /// General test for REVD: /// Computes the decomposition factors, then checks A-U\Sigma\transpose{V}. template - static void test_REVD2_uplo( + static void test_NystromEVD_uplo( int64_t k_start, T tol, T err_expectation, - REVD2UploTestData &all_data, + NystromEVDUploTestData &all_data, algorithm_objects &all_algs, RandBLAS::RNGState state ) { @@ -190,19 +192,21 @@ class TestREVD2 : public ::testing::Test auto m = all_data.dim; int64_t k = k_start; - all_algs.revd2.call(blas::Uplo::Upper, m, all_data.A_u.data(), k, tol, all_data.V_u, all_data.eigvals_u, state); - all_algs.revd2.call(blas::Uplo::Lower, m, all_data.A_l.data(), k, tol, all_data.V_l, all_data.eigvals_l, state); - - T* E_u_dat = RandLAPACK::util::upsize(k * k, all_data.E_u); - T* E_l_dat = RandLAPACK::util::upsize(k * k, all_data.E_l); - T* V_u_dat = all_data.V_u.data(); - T* V_l_dat = all_data.V_l.data(); + all_algs.nystrom_evd.call(blas::Uplo::Upper, m, all_data.A_u.data(), k, tol, + all_data.V_u, all_data.V_u_sz, all_data.eigvals_u, all_data.eigvals_u_sz, state); + all_algs.nystrom_evd.call(blas::Uplo::Lower, m, all_data.A_l.data(), k, tol, + all_data.V_l, all_data.V_l_sz, all_data.eigvals_l, all_data.eigvals_l_sz, state); + + T* E_u_dat = RandLAPACK::util::resize(k * k, all_data.E_u); + T* E_l_dat = RandLAPACK::util::resize(k * k, all_data.E_l); + T* V_u_dat = all_data.V_u; + T* V_l_dat = all_data.V_l; T* work_u_dat = all_data.A_u.data(); T* work_l_dat = all_data.A_l.data(); T* A_approx_dat = all_data.work.data(); - RandLAPACK::util::diag(k, k, all_data.eigvals_u.data(), k, all_data.E_u.data()); - RandLAPACK::util::diag(k, k, all_data.eigvals_l.data(), k, all_data.E_l.data()); + RandLAPACK::util::diag(k, k, all_data.eigvals_u, k, all_data.E_u.data()); + RandLAPACK::util::diag(k, k, all_data.eigvals_l, k, all_data.E_l.data()); // Reconstruct factorizations, compare the result // V_u * E_u = work_u @@ -220,7 +224,7 @@ class TestREVD2 : public ::testing::Test } }; -TEST_F(TestREVD2, Underestimation1) { +TEST_F(TestNystromEVD, Underestimation1) { using RNG = r123::Philox4x32; int64_t m = 1000; @@ -239,7 +243,7 @@ TEST_F(TestREVD2, Underestimation1) { bool verbose = false; bool cond_check = false; - REVD2TestData all_data(m, k); + NystromEVDTestData all_data(m, k); algorithm_objects all_algs( verbose, cond_check, num_syps_passes, @@ -254,12 +258,12 @@ TEST_F(TestREVD2, Underestimation1) { RandLAPACK::gen::mat_gen(m_info, all_data.A_cpy.data(), state); symm_mat_and_copy_computational_helper(norm_A, all_data); - test_REVD2_general( + test_NystromEVD_general( k_start, tol, rank_expectation, err_expectation, norm_A, all_data, all_algs, state ); } -TEST_F(TestREVD2, Underestimation2) { +TEST_F(TestNystromEVD, Underestimation2) { using RNG = r123::Philox4x32; int64_t m = 1000; @@ -277,7 +281,7 @@ TEST_F(TestREVD2, Underestimation2) { bool verbose = false; bool cond_check = false; - REVD2TestData all_data(m, k); + NystromEVDTestData all_data(m, k); algorithm_objects all_algs( verbose, cond_check, num_syps_passes, @@ -292,12 +296,12 @@ TEST_F(TestREVD2, Underestimation2) { RandLAPACK::gen::mat_gen(m_info, all_data.A_cpy.data(), state); symm_mat_and_copy_computational_helper(norm_A, all_data); - test_REVD2_general( + test_NystromEVD_general( k_start, tol, rank_expectation, err_expectation, norm_A, all_data, all_algs, state ); } -TEST_F(TestREVD2, Overestimation1) { +TEST_F(TestNystromEVD, Overestimation1) { using RNG = r123::Philox4x32; int64_t m = 1000; @@ -315,7 +319,7 @@ TEST_F(TestREVD2, Overestimation1) { bool verbose = false; bool cond_check = false; - REVD2TestData all_data(m, k); + NystromEVDTestData all_data(m, k); algorithm_objects all_algs( verbose, cond_check, num_syps_passes, @@ -330,12 +334,12 @@ TEST_F(TestREVD2, Overestimation1) { RandLAPACK::gen::mat_gen(m_info, all_data.A_cpy.data(), state); symm_mat_and_copy_computational_helper(norm_A, all_data); - test_REVD2_general( + test_NystromEVD_general( k_start, tol, rank_expectation, err_expectation, norm_A, all_data, all_algs, state ); } -TEST_F(TestREVD2, Oversetimation2) { +TEST_F(TestNystromEVD, Overestimation2) { using RNG = r123::Philox4x32; int64_t m = 1000; @@ -353,7 +357,7 @@ TEST_F(TestREVD2, Oversetimation2) { bool verbose = false; bool cond_check = false; - REVD2TestData all_data(m, k); + NystromEVDTestData all_data(m, k); algorithm_objects all_algs( verbose, cond_check, num_syps_passes, @@ -368,12 +372,12 @@ TEST_F(TestREVD2, Oversetimation2) { RandLAPACK::gen::mat_gen(m_info, all_data.A_cpy.data(), state); symm_mat_and_copy_computational_helper(norm_A, all_data); - test_REVD2_general( + test_NystromEVD_general( k_start, tol, rank_expectation, err_expectation, norm_A, all_data, all_algs, state ); } -TEST_F(TestREVD2, Exactness) { +TEST_F(TestNystromEVD, Exactness) { using RNG = r123::Philox4x32; int64_t m = 100; @@ -391,7 +395,7 @@ TEST_F(TestREVD2, Exactness) { bool verbose = false; bool cond_check = false; - REVD2TestData all_data(m, k); + NystromEVDTestData all_data(m, k); algorithm_objects all_algs( verbose, cond_check, num_syps_passes, @@ -406,12 +410,35 @@ TEST_F(TestREVD2, Exactness) { RandLAPACK::gen::mat_gen(m_info, all_data.A_cpy.data(), state); symm_mat_and_copy_computational_helper(norm_A, all_data); - test_REVD2_general( + test_NystromEVD_general( k_start, tol, rank_expectation, err_expectation, norm_A, all_data, all_algs, state ); } -TEST_F(TestREVD2, Uplo) { +// Verify that error_est_power_iters=0 and tol=0 produce fixed-rank behavior: +// k must not increase, since FunNystromPP's (n-k)*f(0) tail correction depends on k staying fixed. +TEST_F(TestNystromEVD, FixedRank) { + using RNG = r123::Philox4x32; + int64_t m = 40, k = 10; + auto state = RandBLAS::RNGState(5); + + std::vector A(m * m, 0.0); + for (int64_t i = 0; i < m; ++i) + A[i + i * m] = (double)(i + 1); + + algorithm_objects all_algs(false, false, 3, 1, /*error_est_p=*/0); + double* V_out = nullptr; int64_t V_out_sz = 0; + double* eigvals_out = nullptr; int64_t eigvals_out_sz = 0; + int64_t k_before = k; + all_algs.nystrom_evd.call(blas::Uplo::Upper, m, A.data(), k, 0.0, V_out, V_out_sz, eigvals_out, eigvals_out_sz, state); + delete[] V_out; + delete[] eigvals_out; + + printf("NystromEVD fixed-k: k_before=%lld, k_after=%lld\n", (long long)k_before, (long long)k); + ASSERT_EQ(k, k_before); +} + +TEST_F(TestNystromEVD, Uplo) { using RNG = r123::Philox4x32; int64_t m = 100; @@ -427,7 +454,7 @@ TEST_F(TestREVD2, Uplo) { bool verbose = false; bool cond_check = false; - REVD2UploTestData all_data(m, k); + NystromEVDUploTestData all_data(m, k); algorithm_objects all_algs( verbose, cond_check, num_syps_passes, @@ -443,5 +470,5 @@ TEST_F(TestREVD2, Uplo) { uplo_computational_helper(all_data); - test_REVD2_uplo(k_start, tol, err_expectation, all_data, all_algs, state); + test_NystromEVD_uplo(k_start, tol, err_expectation, all_data, all_algs, state); }