diff --git a/RandLAPACK.hh b/RandLAPACK.hh index 5de765e39..8b872f98d 100644 --- a/RandLAPACK.hh +++ b/RandLAPACK.hh @@ -6,8 +6,11 @@ #include "RandLAPACK/rl_lapackpp.hh" #include "RandBLAS.hh" +// util +#include "RandLAPACK/util/rl_util.hh" +#include "RandLAPACK/util/rl_util_linop.hh" + // misc -#include "RandLAPACK/misc/rl_util.hh" #include "RandLAPACK/misc/rl_linops.hh" #include "RandLAPACK/misc/rl_gen.hh" #include "RandLAPACK/misc/rl_pdkernels.hh" @@ -22,6 +25,7 @@ #include "RandLAPACK/comps/rl_syrf.hh" #include "RandLAPACK/comps/rl_orth.hh" #include "RandLAPACK/comps/rl_rpchol.hh" +#include "RandLAPACK/comps/rl_bk.hh" // Drivers #include "RandLAPACK/drivers/rl_rsvd.hh" diff --git a/RandLAPACK/CMakeLists.txt b/RandLAPACK/CMakeLists.txt index ec4776850..c1edfda54 100644 --- a/RandLAPACK/CMakeLists.txt +++ b/RandLAPACK/CMakeLists.txt @@ -16,10 +16,12 @@ set(RandLAPACK_cxx_sources rl_syps.hh rl_syrf.hh rl_rpchol.hh + rl_bk.hh rl_gen.hh rl_blaspp.hh rl_linops.hh rl_pdkernels.hh + rl_util_linop.hh rl_cusolver.hh rl_cuda_kernels.cuh @@ -61,12 +63,14 @@ target_include_directories( $ $ $ + $ $ $ $ $ $ $ + $ $ ) diff --git a/RandLAPACK/comps/rl_bk.hh b/RandLAPACK/comps/rl_bk.hh new file mode 100644 index 000000000..701fc80ba --- /dev/null +++ b/RandLAPACK/comps/rl_bk.hh @@ -0,0 +1,748 @@ +#pragma once + +#include "rl_util.hh" +#include "rl_blaspp.hh" +#include "rl_lapackpp.hh" +#include "rl_hqrrp.hh" +#include "rl_cqrrt.hh" +#include "rl_linops.hh" + +#include +#include +#include +#include +#include +#include +#include + +using namespace std::chrono; + +namespace RandLAPACK { + +/// BK (Block Krylov) is the computational routine underlying the ABRIK driver. +/// It builds left and right Krylov subspaces (X_ev, Y_od) and band matrices (R, S) +/// via block Krylov iterations with double reorthogonalization. +/// +/// The ABRIK driver calls BK to obtain these factored intermediates, then performs +/// SVD on R or S and reconstructs the final U, Sigma, V. +/// +/// This follows the same pattern as QB (comps) + RSVD (driver). + +// Struct outside of BK class to make symbols shorter +struct BKSubroutines { + enum QR_explicit {geqrf_ungqr, cqrrt}; +}; + +/// Reason BK terminated its main loop. +enum class BKTermination { + max_iters_reached, ///< Reached max_krylov_iters without convergence (resumable). + norm_converged, ///< norm_R exceeded threshold (A's spectral content exhausted). + rank_deficient ///< Near-zero diagonal entry in R or S (subspace can't grow). +}; + +template +class BK { + public: + using Subroutines = BKSubroutines; + Subroutines::QR_explicit qr_exp; + + bool verbose; + bool timing; + T tol; + int num_krylov_iters; + int max_krylov_iters; + std::vector times; + T norm_R_end; + BKTermination termination_reason; + + BK( + bool verb, + bool time_subroutines, + T ep + ) { + qr_exp = Subroutines::QR_explicit::geqrf_ungqr; + verbose = verb; + timing = time_subroutines; + tol = ep; + max_krylov_iters = INT_MAX; + } + + /// Builds the block Krylov subspaces and band matrices for a truncated SVD. + /// + /// @param[in] m + /// The number of rows in the matrix A. + /// + /// @param[in] n + /// The number of columns in the matrix A. + /// + /// @param[in] A + /// Pointer to the m-by-n matrix A, stored in a column-major format. + /// + /// @param[in] lda + /// Leading dimension of A. + /// + /// @param[in] k + /// Block size for Krylov iterations. + /// + /// @param[out] X_ev + /// Left Krylov basis (m x end_rows), allocated internally with calloc. + /// Caller must free(). + /// + /// @param[out] Y_od + /// Right Krylov basis (n x end_cols), allocated internally with calloc. + /// Caller must free(). + /// + /// @param[out] R + /// Upper band matrix (stored transposed), allocated internally with calloc. + /// Caller must free(). + /// + /// @param[out] S + /// Lower Hessenberg band matrix, allocated internally with calloc. + /// Caller must free(). + /// + /// @param[out] end_rows + /// Number of rows in the band matrix for SVD. + /// + /// @param[out] end_cols + /// Number of columns in the band matrix for SVD. + /// + /// @param[out] final_iter_is_odd + /// True if the last iteration was odd (use R for SVD), false if even (use S). + /// + /// @param[in] state + /// RNG state parameter, required for sketching operator generation. + /// + /// @return = 0: successful exit, -1: realloc failure + + // BK call that accepts a general dense matrix. + int call( + int64_t m, + int64_t n, + T* A, + int64_t lda, + int64_t k, + T* &X_ev, + T* &Y_od, + T* &R, + T* &S, + int64_t &end_rows, + int64_t &end_cols, + bool &final_iter_is_odd, + RandBLAS::RNGState &state + ) { + linops::GenLinOp A_linop(m, n, A, lda, Layout::ColMajor); + return this->call(A_linop, k, X_ev, Y_od, R, S, end_rows, end_cols, final_iter_is_odd, state); + } + + // BK call that accepts sparse matrix. + template + int call( + int64_t m, + int64_t n, + SpMat &A, + int64_t lda, + int64_t k, + T* &X_ev, + T* &Y_od, + T* &R, + T* &S, + int64_t &end_rows, + int64_t &end_cols, + bool &final_iter_is_odd, + RandBLAS::RNGState &state + ) { + linops::SpLinOp A_linop(m, n, A); + return this->call(A_linop, k, X_ev, Y_od, R, S, end_rows, end_cols, final_iter_is_odd, state); + } + + /// Resume a previous BK computation with more iterations. + /// X_ev, Y_od, R, S must be non-null from a prior call(). + /// Increase max_krylov_iters before calling. + template + int resume( + GLO& A, + int64_t k, + T* &X_ev, + T* &Y_od, + T* &R, + T* &S, + int64_t &end_rows, + int64_t &end_cols, + bool &final_iter_is_odd, + RandBLAS::RNGState &state + ) { + return this->call_impl(A, k, X_ev, Y_od, R, S, end_rows, end_cols, final_iter_is_odd, state, true); + } + + template + int call( + GLO& A, + int64_t k, + T* &X_ev, + T* &Y_od, + T* &R, + T* &S, + int64_t &end_rows, + int64_t &end_cols, + bool &final_iter_is_odd, + RandBLAS::RNGState &state + ) { + return this->call_impl(A, k, X_ev, Y_od, R, S, end_rows, end_cols, final_iter_is_odd, state, false); + } + + private: + template + int call_impl( + GLO& A, + int64_t k, + T* &X_ev, + T* &Y_od, + T* &R, + T* &S, + int64_t &end_rows, + int64_t &end_cols, + bool &final_iter_is_odd, + RandBLAS::RNGState &state, + bool resuming + ){ + steady_clock::time_point allocation_t_start; + steady_clock::time_point allocation_t_stop; + steady_clock::time_point ungqr_t_start; + steady_clock::time_point ungqr_t_stop; + steady_clock::time_point reorth_t_start; + steady_clock::time_point reorth_t_stop; + steady_clock::time_point qr_t_start; + steady_clock::time_point qr_t_stop; + steady_clock::time_point gemm_A_t_start; + steady_clock::time_point gemm_A_t_stop; + steady_clock::time_point main_loop_t_start; + steady_clock::time_point main_loop_t_stop; + steady_clock::time_point sketching_t_start; + steady_clock::time_point sketching_t_stop; + steady_clock::time_point r_cpy_t_start; + steady_clock::time_point r_cpy_t_stop; + steady_clock::time_point s_cpy_t_start; + steady_clock::time_point s_cpy_t_stop; + steady_clock::time_point norm_t_start; + steady_clock::time_point norm_t_stop; + steady_clock::time_point bk_total_t_start; + steady_clock::time_point bk_total_t_stop; + + long allocation_t_dur = 0; + long ungqr_t_dur = 0; + long reorth_t_dur = 0; + long qr_t_dur = 0; + long gemm_A_t_dur = 0; + long main_loop_t_dur = 0; + long sketching_t_dur = 0; + long r_cpy_t_dur = 0; + long s_cpy_t_dur = 0; + long norm_t_dur = 0; + long bk_total_t_dur = 0; + + if(this -> timing) + bk_total_t_start = steady_clock::now(); + + int64_t m = A.n_rows; + int64_t n = A.n_cols; + int max_iters = this->max_krylov_iters; + + // Loop state — initialized differently for fresh start vs resume. + int64_t iter, iter_od, iter_ev; + int64_t curr_X_cols, curr_Y_cols; + T norm_R; + T* Y_i; + T* X_i; + T* R_i; + T* R_ii; + T* S_i; + T* S_ii; + + // Pre-allocation: when max_krylov_iters is known, allocate all + // buffers upfront to avoid per-iteration realloc + memset. + bool prealloc = (max_iters != INT_MAX); + int64_t max_X_cols = 0, max_Y_cols = 0; + if (prealloc) { + // After max_iters iterations: + // odd iters (1,3,...) grow X_ev; even iters (2,4,...) grow Y_od + // Initial: k cols each. Each relevant iter adds k cols. + int64_t n_odd = (max_iters + 1) / 2; // ceil(max_iters/2) + int64_t n_even = max_iters / 2; + max_X_cols = k * (1 + n_odd); + max_Y_cols = k * (1 + n_even); + } + + if (!resuming) { + // --- Fresh start: allocate output buffers and initialize state --- + if(this -> timing) + allocation_t_start = steady_clock::now(); + + iter = 0; iter_od = 0; iter_ev = 0; + end_rows = 0; end_cols = 0; + norm_R = 0; + + if (prealloc) { + // Allocate to maximum size upfront — no realloc needed in loop. + Y_od = ( T * ) calloc( n * max_Y_cols, sizeof( T ) ); + X_ev = ( T * ) calloc( m * max_X_cols, sizeof( T ) ); + R = ( T * ) calloc( n * max_X_cols, sizeof( T ) ); + S = ( T * ) calloc( (n + k) * max_Y_cols, sizeof( T ) ); + } else { + // Tolerance-based: start small, realloc as needed. + Y_od = ( T * ) calloc( n * k, sizeof( T ) ); + X_ev = ( T * ) calloc( m * k, sizeof( T ) ); + R = ( T * ) calloc( n * k, sizeof( T ) ); + S = ( T * ) calloc( (n + k) * k, sizeof( T ) ); + } + curr_Y_cols = k; + curr_X_cols = k; + + // Initialize pointers. + Y_i = Y_od; + X_i = X_ev; + R_i = NULL; + R_ii = R; + S_i = S; + S_ii = &S[k]; + + if(this -> timing) { + allocation_t_stop = steady_clock::now(); + allocation_t_dur = duration_cast(allocation_t_stop - allocation_t_start).count(); + } + } else { + // --- Resume: reconstruct loop state from stored members --- + // Only valid after a prior call() that terminated with max_iters_reached. + iter = this->num_krylov_iters; + norm_R = this->norm_R_end; + iter_od = 1 + iter / 2; + iter_ev = (iter + 1) / 2; + curr_X_cols = (1 + iter_ev) * k; + curr_Y_cols = iter_od * k; + + // Grow buffers if the new max_krylov_iters requires more space + // than was allocated in the previous call. + if (prealloc) { + if (max_X_cols > curr_X_cols) { + X_ev = ( T * ) realloc(X_ev, m * max_X_cols * sizeof( T )); + R = ( T * ) realloc(R, n * max_X_cols * sizeof( T )); + if (!X_ev || !R) { + free(X_ev); free(Y_od); free(R); free(S); + X_ev = nullptr; Y_od = nullptr; R = nullptr; S = nullptr; + return -1; + } + std::fill(&X_ev[m * curr_X_cols], &X_ev[m * max_X_cols], T(0)); + std::fill(&R[n * curr_X_cols], &R[n * max_X_cols], T(0)); + } + if (max_Y_cols > curr_Y_cols) { + Y_od = ( T * ) realloc(Y_od, n * max_Y_cols * sizeof( T )); + S = ( T * ) realloc(S, (n + k) * max_Y_cols * sizeof( T )); + if (!Y_od || !S) { + free(X_ev); free(Y_od); free(R); free(S); + X_ev = nullptr; Y_od = nullptr; R = nullptr; S = nullptr; + return -1; + } + std::fill(&Y_od[n * curr_Y_cols], &Y_od[n * max_Y_cols], T(0)); + std::fill(&S[(n + k) * curr_Y_cols], &S[(n + k) * max_Y_cols], T(0)); + } + } + + // Reconstruct pointers into (potentially reallocated) buffers. + X_i = &X_ev[m * (curr_X_cols - k)]; + Y_i = &Y_od[n * (curr_Y_cols - k)]; + R_i = &R[iter_ev * k]; + R_ii = &R[(n * k * iter_ev) + k + (k * (iter_ev - 1))]; + S_i = &S[(n + k) * k * (iter_od - 1)]; + S_ii = &S[(n + k) * k * (iter_od - 1) + k + ((iter_od - 1) * k)]; + + // Advance past the completed iteration so the while-loop starts at the next one. + ++iter; + } + + // Internal temporaries — shared for both paths. + // These are pure scratch buffers (beta=0.0 GEMM outputs), no need to zero-initialize. + T* Y_orth_buf = ( T * ) malloc( k * n * sizeof( T ) ); + T* X_orth_buf = ( T * ) malloc( k * (n + k) * sizeof( T ) ); + // tau space for QR (geqrf fully overwrites it) + T* tau = ( T * ) malloc( k * sizeof( T ) ); + // Declared here (before cleanup lambda) so cleanup can free it. + // Conditionally allocated below only when CQRRT is used. + T* R_11_trans = nullptr; + + // Cleanup lambda for realloc failure — frees all buffers and nulls output pointers. + // free(nullptr) is a no-op, so no guards needed. + auto cleanup_and_fail = [&]() -> int { + free(Y_od); Y_od = nullptr; + free(X_ev); X_ev = nullptr; + free(R); R = nullptr; + free(S); S = nullptr; + free(tau); + free(Y_orth_buf); + free(X_orth_buf); + free(R_11_trans); + return -1; + }; + + // Pre-compute Fro norm of an input matrix. + T norm_A = A.fro_nrm(); + T sq_tol = std::pow(this->tol, 2); + T threshold = std::sqrt(1 - sq_tol) * norm_A; + + // Creating the CQRRT object in case it is to be used for explicit QR. + std::optional> CQRRT; + T d_factor = 1.25; + // Conditional initialization + if(this -> qr_exp == Subroutines::QR_explicit::cqrrt) { + CQRRT.emplace(false, tol); + CQRRT->nnz = 2; + R_11_trans = ( T * ) calloc( k * k, sizeof( T ) ); + } + + if (!resuming) { + // --- Fresh start: sketch generation, first GEMM, first QR --- + if(this -> timing) + sketching_t_start = steady_clock::now(); + + // Generate a dense Gaussian random matrix. + RandBLAS::DenseDist D(n, k); + state = RandBLAS::fill_dense(D, Y_i, state); + + if(this -> timing) { + sketching_t_stop = steady_clock::now(); + sketching_t_dur = duration_cast(sketching_t_stop - sketching_t_start).count(); + gemm_A_t_start = steady_clock::now(); + } + + // [X_ev, ~] = qr(A * Y_i, 0) + A(Side::Left, Layout::ColMajor, Op::NoTrans, Op::NoTrans, m, k, n, 1.0, Y_i, n, 0.0, X_i, m); + + if(this -> timing) { + gemm_A_t_stop = steady_clock::now(); + gemm_A_t_dur = duration_cast(gemm_A_t_stop - gemm_A_t_start).count(); + } + + if(this -> qr_exp == Subroutines::QR_explicit::cqrrt) { + if(this -> timing) + qr_t_start = steady_clock::now(); + + CQRRT -> call(m, k, X_i, m, R_11_trans, k, d_factor, state); + + if(this -> timing) { + qr_t_stop = steady_clock::now(); + qr_t_dur = duration_cast(qr_t_stop - qr_t_start).count(); + } + } else { + + if(this -> timing) + qr_t_start = steady_clock::now(); + + lapack::geqrf(m, k, X_i, m, tau); + + if(this -> timing) { + qr_t_stop = steady_clock::now(); + qr_t_dur = duration_cast(qr_t_stop - qr_t_start).count(); + ungqr_t_start = steady_clock::now(); + } + + // Convert X_i into an explicit form. It is now stored in X_ev as it should be. + lapack::ungqr(m, k, k, X_i, m, tau); + + if(this -> timing) { + ungqr_t_stop = steady_clock::now(); + ungqr_t_dur += duration_cast(ungqr_t_stop - ungqr_t_start).count(); + } + } + + // Advance odd iteration count. + ++iter_od; + // Advance iteration count. + ++iter; + } + + // Main loop — shared for both fresh start and resume. + while(1) { + if(this -> timing) + main_loop_t_start = steady_clock::now(); + + if (iter % 2 != 0) { + if(this -> timing) + gemm_A_t_start = steady_clock::now(); + // Y_i = A' * X_i + A(Side::Left, Layout::ColMajor, Op::Trans, Op::NoTrans, n, k, m, 1.0, X_i, m, 0.0, Y_i, n); + + if(this -> timing) { + gemm_A_t_stop = steady_clock::now(); + gemm_A_t_dur += duration_cast(gemm_A_t_stop - gemm_A_t_start).count(); + allocation_t_start = steady_clock::now(); + } + + // Grow X_ev buffer + curr_X_cols += k; + if (!prealloc) { + X_ev = ( T * ) realloc(X_ev, m * curr_X_cols * sizeof( T )); + } + // Move the X_i pointer + X_i = &X_ev[m * (curr_X_cols - k)]; + + + if(this -> timing) { + allocation_t_stop = steady_clock::now(); + allocation_t_dur += duration_cast(allocation_t_stop - allocation_t_start).count(); + reorth_t_start = steady_clock::now(); + } + + if (iter != 1) { + // R_i' = Y_i' * Y_od + blas::gemm(Layout::ColMajor, Op::Trans, Op::NoTrans, k, iter_ev * k, n, 1.0, Y_i, n, Y_od, n, 0.0, R_i, n); + + // Y_i = Y_i - Y_od * R_i + blas::gemm(Layout::ColMajor, Op::NoTrans, Op::Trans, n, k, iter_ev * k, -1.0, Y_od, n, R_i, n, 1.0, Y_i, n); + + // Reorthogonalization + blas::gemm(Layout::ColMajor, Op::Trans, Op::NoTrans, k, iter_ev * k, n, 1.0, Y_i, n, Y_od, n, 0.0, Y_orth_buf, k); + blas::gemm(Layout::ColMajor, Op::NoTrans, Op::Trans, n, k, iter_ev * k, -1.0, Y_od, n, Y_orth_buf, k, 1.0, Y_i, n); + } + + if(this -> timing) { + reorth_t_stop = steady_clock::now(); + reorth_t_dur += duration_cast(reorth_t_stop - reorth_t_start).count(); + } + + // Perform explicit QR via a method of choice + if(this -> qr_exp == Subroutines::QR_explicit::cqrrt) { + if(this -> timing) + qr_t_start = steady_clock::now(); + + CQRRT -> call(n, k, Y_i, n, R_11_trans, k, d_factor, state); + // Copy R_ii over to R's (in transposed format). + + util::transposition(0, k, R_11_trans, k, R_ii, n, 1); + if(this -> timing) { + qr_t_stop = steady_clock::now(); + qr_t_dur += duration_cast(qr_t_stop - qr_t_start).count(); + } + } else { + // [Y_i, R_ii] = qr(Y_i, 0) + if(this -> timing) + qr_t_start = steady_clock::now(); + lapack::geqrf(n, k, Y_i, n, tau); + + if(this -> timing) { + qr_t_stop = steady_clock::now(); + qr_t_dur += duration_cast(qr_t_stop - qr_t_start).count(); + r_cpy_t_start = steady_clock::now(); + } + + // Copy R_ii over to R's (in transposed format). + util::transposition(0, k, Y_i, n, R_ii, n, 1); + + if(this -> timing) { + r_cpy_t_stop = steady_clock::now(); + r_cpy_t_dur += duration_cast(r_cpy_t_stop - r_cpy_t_start).count(); + ungqr_t_start = steady_clock::now(); + } + + // Convert Y_i into an explicit form. It is now stored in Y_odd as it should be. + lapack::ungqr(n, k, k, Y_i, n, tau); + + if(this -> timing) { + ungqr_t_stop = steady_clock::now(); + ungqr_t_dur += duration_cast(ungqr_t_stop - ungqr_t_start).count(); + } + } + + // Early termination + // if (abs(R(end)) <= sqrt(eps('T'))) + if(std::abs(R_ii[(n + 1) * (k - 1)]) < std::sqrt(std::numeric_limits::epsilon())) { + this->termination_reason = BKTermination::rank_deficient; + break; + } + + // Grow R buffer + if (!prealloc) { + T* R_new = ( T * ) realloc(R, n * curr_X_cols * sizeof( T )); + if (!R_new) return cleanup_and_fail(); + R = R_new; + T* temp_r = &R[n * (curr_X_cols - k)]; + std::fill(temp_r, temp_r + n*k, 0.0); + } + + // Advance R pointers + R_i = &R[(iter_ev + 1) * k]; + R_ii = &R[(n * k * (iter_ev + 1)) + k + (k * (iter_ev))]; + + // Advance even iteration count; + ++iter_ev; + } + else { + if(this -> timing) + gemm_A_t_start = steady_clock::now(); + + // X_i = A * Y_i + A(Side::Left, Layout::ColMajor, Op::NoTrans, Op::NoTrans, m, k, n, 1.0, Y_i, n, 0.0, X_i, m); + + if(this -> timing) { + gemm_A_t_stop = steady_clock::now(); + gemm_A_t_dur += duration_cast(gemm_A_t_stop - gemm_A_t_start).count(); + allocation_t_start = steady_clock::now(); + } + + // Grow Y_od buffer + curr_Y_cols += k; + if (!prealloc) { + Y_od = ( T * ) realloc(Y_od, n * curr_Y_cols * sizeof( T )); + } + // Move the Y_i pointer + Y_i = &Y_od[n * (curr_Y_cols - k)]; + + if(this -> timing) { + allocation_t_stop = steady_clock::now(); + allocation_t_dur += duration_cast(allocation_t_stop - allocation_t_start).count(); + reorth_t_start = steady_clock::now(); + } + + // S_i = X_ev' * X_i + blas::gemm(Layout::ColMajor, Op::Trans, Op::NoTrans, iter_od * k, k, m, 1.0, X_ev, m, X_i, m, 0.0, S_i, n + k); + + //X_i = X_i - X_ev * S_i; + blas::gemm(Layout::ColMajor, Op::NoTrans, Op::NoTrans, m, k, iter_od * k, -1.0, X_ev, m, S_i, n + k, 1.0, X_i, m); + + // Reorthogonalization + blas::gemm(Layout::ColMajor, Op::Trans, Op::NoTrans, iter_od * k, k, m, 1.0, X_ev, m, X_i, m, 0.0, X_orth_buf, n + k); + blas::gemm(Layout::ColMajor, Op::NoTrans, Op::NoTrans, m, k, iter_od * k, -1.0, X_ev, m, X_orth_buf, n + k, 1.0, X_i, m); + + if(this -> timing) { + reorth_t_stop = steady_clock::now(); + reorth_t_dur += duration_cast(reorth_t_stop - reorth_t_start).count(); + } + + // Perform explicit QR via a method of choice + if(this -> qr_exp == Subroutines::QR_explicit::cqrrt) { + if(this -> timing) + qr_t_start = steady_clock::now(); + + CQRRT -> call(m, k, X_i, m, S_ii, n + k, d_factor, state); + + if(this -> timing) { + qr_t_stop = steady_clock::now(); + qr_t_dur += duration_cast(qr_t_stop - qr_t_start).count(); + } + + } else { + // [X_i, S_ii] = qr(X_i, 0); + if(this -> timing) + qr_t_start = steady_clock::now(); + + lapack::geqrf(m, k, X_i, m, tau); + + if(this -> timing) { + qr_t_stop = steady_clock::now(); + qr_t_dur += duration_cast(qr_t_stop - qr_t_start).count(); + s_cpy_t_start = steady_clock::now(); + } + + // Copy S_ii over to S's space under S_i (offset down by iter_od * k) + lapack::lacpy(MatrixType::Upper, k, k, X_i, m, S_ii, n + k); + + if(this -> timing) { + s_cpy_t_stop = steady_clock::now(); + s_cpy_t_dur += duration_cast(s_cpy_t_stop - s_cpy_t_start).count(); + ungqr_t_start = steady_clock::now(); + } + + // Convert X_i into an explicit form. It is now stored in X_ev as it should be + lapack::ungqr(m, k, k, X_i, m, tau); + + if(this -> timing) { + ungqr_t_stop = steady_clock::now(); + ungqr_t_dur += duration_cast(ungqr_t_stop - ungqr_t_start).count(); + } + } + + // Early termination + // if (abs(S(end)) <= sqrt(eps('T'))) + if(std::abs(S_ii[((n + k) + 1) * (k - 1)]) < std::sqrt(std::numeric_limits::epsilon())) { + this->termination_reason = BKTermination::rank_deficient; + break; + } + + if(this -> timing) { + allocation_t_start = steady_clock::now(); + } + + // Grow S buffer + if (!prealloc) { + T* S_new = ( T * ) realloc(S, (n + k) * curr_Y_cols * sizeof( T )); + if (!S_new) return cleanup_and_fail(); + S = S_new; + T* temp_s = &S[(n + k)* (curr_Y_cols - k)]; + std::fill(temp_s, temp_s + (n + k) * k, 0.0); + } + + // Advance S pointers + S_i = &S[(n + k) * k * iter_od]; + S_ii = &S[(n + k) * k * iter_od + k + (iter_od * k)]; + + // Advance odd iteration count; + ++iter_od; + + if(this -> timing) { + allocation_t_stop = steady_clock::now(); + allocation_t_dur += duration_cast(allocation_t_stop - allocation_t_start).count(); + } + } + + if(this -> timing) + norm_t_start = steady_clock::now(); + + // This is only changed on odd iters + if (iter % 2 != 0) + norm_R = lapack::lantr(Norm::Fro, Uplo::Upper, Diag::NonUnit, iter_ev * k, iter_ev * k, R, n); + + if(this -> timing) { + norm_t_stop = steady_clock::now(); + norm_t_dur += duration_cast(norm_t_stop - norm_t_start).count(); + main_loop_t_stop = steady_clock::now(); + main_loop_t_dur += duration_cast(main_loop_t_stop - main_loop_t_start).count(); + } + + if (iter >= max_iters) { + this->termination_reason = BKTermination::max_iters_reached; + break; + } + + ++iter; + //norm(R, 'fro') > sqrt(1 - sq_tol) * norm_A + if(norm_R > threshold) { + this->termination_reason = BKTermination::norm_converged; + break; + } + } + + // Set output state + this->norm_R_end = norm_R; + this->num_krylov_iters = iter; + end_cols = iter * k / 2; + iter % 2 == 0 ? end_rows = end_cols + k : end_rows = end_cols; + final_iter_is_odd = (iter % 2 != 0); + + // Free internal temporaries (NOT X_ev, Y_od, R, S — those are returned to caller) + free(tau); + free(Y_orth_buf); + free(X_orth_buf); + if(R_11_trans != nullptr) { + free(R_11_trans); + } + + if(this -> timing) { + bk_total_t_stop = steady_clock::now(); + bk_total_t_dur = duration_cast(bk_total_t_stop - bk_total_t_start).count(); + + this -> times.resize(10); + this -> times = {allocation_t_dur, ungqr_t_dur, reorth_t_dur, qr_t_dur, + gemm_A_t_dur, main_loop_t_dur, sketching_t_dur, + r_cpy_t_dur, s_cpy_t_dur, norm_t_dur}; + } + return 0; + } + }; +} diff --git a/RandLAPACK/comps/rl_qb.hh b/RandLAPACK/comps/rl_qb.hh index 91e43619f..a689679a5 100644 --- a/RandLAPACK/comps/rl_qb.hh +++ b/RandLAPACK/comps/rl_qb.hh @@ -159,16 +159,15 @@ int QB::call( BT = ( T * ) calloc(n * b_sz, sizeof( T ) ); // Allocate buffers T* QtQi = ( T * ) calloc( b_sz * b_sz, sizeof( T ) ); - T* A_cpy = ( T * ) calloc( m * n, sizeof( T ) ); // Declate pointers to the iteration buffers. T* Q_i; T* BT_i; - // pre-compute nrom + // pre-compute norm T norm_A = lapack::lange(Norm::Fro, m, n, A, m); - // Copy the initial data to avoid unwanted modification - lapack::lacpy(MatrixType::General, m, n, A, m, A_cpy, m); + // NOTE: A is modified in-place by the deflation step (A = A - Q_i * B_i). + // Callers who need to preserve A must make their own copy before calling QB. while(curr_sz < k) { // Dynamically changing block size. @@ -188,10 +187,9 @@ int QB::call( BT_i = &BT[n * curr_sz]; // Calling RangeFinder - if(this->rf.call(m, n, A_cpy, b_sz, Q_i, state)) { + if(this->rf.call(m, n, A, b_sz, Q_i, state)) { // RF failed k = curr_sz; - free(A_cpy); free(QtQi); return 6; } @@ -200,7 +198,6 @@ int QB::call( if (util::orthogonality_check(m, b_sz, Q_i, this->verbose)) { // Lost orthonormality of Q k = curr_sz; - free(A_cpy); free(QtQi); return 4; } @@ -215,7 +212,7 @@ int QB::call( } //B_i' = A' * Q_i' - blas::gemm(Layout::ColMajor, Op::Trans, Op::NoTrans, n, b_sz, m, 1.0, A_cpy, m, Q_i, m, 0.0, BT_i, n); + blas::gemm(Layout::ColMajor, Op::Trans, Op::NoTrans, n, b_sz, m, 1.0, A, m, Q_i, m, 0.0, BT_i, n); // Updating B norm estimation T norm_B_i = lapack::lange(Norm::Fro, n, b_sz, BT_i, n); @@ -228,7 +225,6 @@ int QB::call( if ((curr_sz > 0) && (approx_err > prev_err)) { // Early termination - error has grown. k = curr_sz; - free(A_cpy); free(QtQi); return 2; } @@ -237,7 +233,6 @@ int QB::call( if (util::orthogonality_check(m, next_sz, Q, this->verbose)) { // Lost orthonormality of Q k = curr_sz; - free(A_cpy); free(QtQi); return 5; } @@ -250,17 +245,15 @@ int QB::call( if (approx_err < tol) { // Reached the required error tol k = curr_sz; - free(A_cpy); free(QtQi); return 0; } // This step is only necessary for the next iteration // A = A - Q_i * B_i - blas::gemm(Layout::ColMajor, Op::NoTrans, Op::Trans, m, n, b_sz, -1.0, Q_i, m, BT_i, n, 1.0, A_cpy, m); + blas::gemm(Layout::ColMajor, Op::NoTrans, Op::Trans, m, n, b_sz, -1.0, Q_i, m, BT_i, n, 1.0, A, m); } - free(A_cpy); free(QtQi); // Reached expected rank without achieving the tolerance diff --git a/RandLAPACK/drivers/rl_abrik.hh b/RandLAPACK/drivers/rl_abrik.hh index 18c76d9fc..0ef823e4f 100644 --- a/RandLAPACK/drivers/rl_abrik.hh +++ b/RandLAPACK/drivers/rl_abrik.hh @@ -1,18 +1,17 @@ #pragma once -#include "rl_util.hh" +#include "rl_bk.hh" #include "rl_blaspp.hh" #include "rl_lapackpp.hh" -#include "rl_hqrrp.hh" +#include "rl_util.hh" #include "rl_linops.hh" +#include "rl_util_linop.hh" #include #include #include #include -#include #include -#include using namespace std::chrono; @@ -20,22 +19,24 @@ namespace RandLAPACK { /// ABRIK algorithm is a method for finding truncated SVD based on block Krylov iterations. /// This algorithm is a version of Algroithm A.1 from https://arxiv.org/pdf/2306.12418.pdf - /// - /// The main difference is in the fact that an economy SVD is performed only once at the very end + /// + /// The main difference is in the fact that an economy SVD is performed only once at the very end /// of the algorithm run and that the termination criteria is not based on singular vectir residual evaluation. /// Instead, the scheme terminates if: - /// 1. ||R||_F > sqrt(1 - eps^2) ||A||_F, which ensures that we've exhausted all vectors and doing more + /// 1. ||R||_F > sqrt(1 - eps^2) ||A||_F, which ensures that we've exhausted all vectors and doing more /// iterations would bring no benefit or that ||A - hat(A)||_F < eps * ||A||_F. /// 2. Stop if the bottom right entry of R or S is numerically close to zero (up to square root of machine eps). - /// + /// /// The main cost of this algorithm comes from large GEMMs with the input matrix A. /// /// The algorithm optionally times all of its subcomponents through a user-defined 'timing' parameter. + /// + /// ABRIK is a driver that delegates the block Krylov iteration to the BK computational routine, + /// then performs SVD on the resulting band matrix and reconstructs the final U, Sigma, V factors. + /// This follows the same pattern as RSVD (driver) + QB (comp). -// Struct outside of ABRIK class to make symbols shorter -struct ABRIKSubroutines { - enum QR_explicit {geqrf_ungqr, cqrrt}; -}; +// Backward compatibility alias +using ABRIKSubroutines = BKSubroutines; template class ABRIK { @@ -52,26 +53,27 @@ class ABRIK { std::vector times; T norm_R_end; - // Numbr of threads that will be used in - // functions where parallelism can tank performance. - int num_threads_min; - // Number of threads used in the rest of the code. - int num_threads_max; int64_t singular_triplets_found; + // Adaptive mode: check SVD residual after BK and resume if needed. + bool adaptive; // Enable adaptive residual checking (default: false). + int adaptive_increment; // Extra BK iterations per retry (0 = use max_krylov_iters). + int adaptive_max_retries; // Hard limit on resume attempts (default: 10). + ABRIK( bool verb, bool time_subroutines, T ep - ) { + ) : bk_obj(verb, time_subroutines, ep) { qr_exp = Subroutines::QR_explicit::geqrf_ungqr; verbose = verb; timing = time_subroutines; tol = ep; max_krylov_iters = INT_MAX; - num_threads_min = util::get_omp_threads(); - num_threads_max = util::get_omp_threads(); singular_triplets_found = 0; + adaptive = false; + adaptive_increment = 0; + adaptive_max_retries = 10; } /// Computes an SVD of the form: @@ -111,7 +113,7 @@ class ABRIK { /// Stores n by ((num_iters / 2) * k) orthonormal matrix of right singular vectors. /// /// @param[out] Sigma - /// Stores ((num_iters / 2) * k) singular values. + /// Stores ((num_iters / 2) * k) singular values. /// /// @return = 0: successful exit /// @@ -158,547 +160,202 @@ class ABRIK { T* &Sigma, RandBLAS::RNGState &state ){ - steady_clock::time_point allocation_t_start; - steady_clock::time_point allocation_t_stop; - steady_clock::time_point get_factors_t_start; - steady_clock::time_point get_factors_t_stop; - steady_clock::time_point ungqr_t_start; - steady_clock::time_point ungqr_t_stop; - steady_clock::time_point reorth_t_start; - steady_clock::time_point reorth_t_stop; - steady_clock::time_point qr_t_start; - steady_clock::time_point qr_t_stop; - steady_clock::time_point gemm_A_t_start; - steady_clock::time_point gemm_A_t_stop; - steady_clock::time_point main_loop_t_start; - steady_clock::time_point main_loop_t_stop; - steady_clock::time_point sketching_t_start; - steady_clock::time_point sketching_t_stop; - steady_clock::time_point r_cpy_t_start; - steady_clock::time_point r_cpy_t_stop; - steady_clock::time_point s_cpy_t_start; - steady_clock::time_point s_cpy_t_stop; - steady_clock::time_point norm_t_start; - steady_clock::time_point norm_t_stop; steady_clock::time_point total_t_start; steady_clock::time_point total_t_stop; - - long allocation_t_dur = 0; + steady_clock::time_point get_factors_t_start; + steady_clock::time_point get_factors_t_stop; + steady_clock::time_point allocation_t_start; + steady_clock::time_point allocation_t_stop; long get_factors_t_dur = 0; - long ungqr_t_dur = 0; - long reorth_t_dur = 0; - long qr_t_dur = 0; - long gemm_A_t_dur = 0; - long main_loop_t_dur = 0; - long sketching_t_dur = 0; - long r_cpy_t_dur = 0; - long s_cpy_t_dur = 0; - long norm_t_dur = 0; - long total_t_dur = 0; + long driver_alloc_t_dur = 0; + long total_t_dur = 0; - if(this -> timing) { + if(this -> timing) total_t_start = steady_clock::now(); - allocation_t_start = steady_clock::now(); - } - - int64_t m = A.n_rows; - int64_t n = A.n_cols; - int64_t iter = 0, iter_od = 0, iter_ev = 0, end_rows = 0, end_cols = 0; - T norm_R = 0; - int max_iters = this->max_krylov_iters;//std::min(this->max_krylov_iters, (int) (n / (T) k)); - - // We need a full copy of X and Y all the way through the algorithm - // due to an operation with X_odd and Y_odd happening at the end. - // Below pointers stay the same throughout the alg; the space will be alloacted iteratively - // Space for Y_i and Y_odd. - T* Y_od = ( T * ) calloc( n * k, sizeof( T ) ); - int64_t curr_Y_cols = k; - // Space for X_i and X_ev. - T* X_ev = ( T * ) calloc( m * k, sizeof( T ) ); - int64_t curr_X_cols = k; - - // While R and S matrices are structured (both band), we cannot make use of this structure through - // BLAS-level functions. - // Note also that we store a transposed version of R. - // - // At each iterations, matrices R and S grow by b_sz. - // At the end, size of R would by d x d and size of S would - // be (d + 1) x d, where d = numiters_complete * b_sz, d <= n. - // Note that the total amount of iterations will always be numiters <= n * 2 / block_size - T* R = ( T * ) calloc( n * k, sizeof( T ) ); - T* S = ( T * ) calloc( (n + k) * k, sizeof( T ) ); - - // These buffers are of constant size - T* Y_orth_buf = ( T * ) calloc( k * n, sizeof( T ) ); - T* X_orth_buf = ( T * ) calloc( k * (n + k), sizeof( T ) ); - - // Pointers allocation - // Below pointers will be offset by (n or m) * k at every even iteration. - T* Y_i = Y_od; - T* X_i = X_ev; - // S and S pointers are offset at every step. - T* R_i = NULL; - T* R_ii = R; - T* S_i = S; - T* S_ii = &S[k]; - // Pre-decloration of SVD-related buffers. - T* U_hat = NULL; - T* VT_hat = NULL; - // tau space for QR - T* tau = ( T * ) calloc( k, sizeof( T ) ); - - if(this -> timing) { - allocation_t_stop = steady_clock::now(); - allocation_t_dur = duration_cast(allocation_t_stop - allocation_t_start).count(); - } - // Pre-compute Fro norm of an input matrix. - //T norm_A = lapack::lange(Norm::Fro, m, n, A.A_buff, lda); - T norm_A = A.fro_nrm(); - T sq_tol = std::pow(this->tol, 2); - T threshold = std::sqrt(1 - sq_tol) * norm_A; - - // Creating the CQRRT object in case it is to be used for explicit QR. - std::optional> CQRRT; - T* R_11_trans = nullptr; - T d_factor = 1.25; - // Conditional initialization - if(this -> qr_exp == Subroutines::QR_explicit::cqrrt) { - CQRRT.emplace(false, tol); - CQRRT->nnz = 2; - R_11_trans = ( T * ) calloc( k * k, sizeof( T ) ); - } + // Forward config to BK + bk_obj.qr_exp = this->qr_exp; + bk_obj.tol = this->tol; + bk_obj.max_krylov_iters = this->max_krylov_iters; + bk_obj.verbose = this->verbose; + bk_obj.timing = this->timing; - if(this -> timing) - sketching_t_start = steady_clock::now(); - - // Generate a dense Gaussian random matrix. - // We are using the plain dense operator instead of DenseSkOp here since - // the space in which teh dense operator is stored will be reused later, and - // also needs to be used together with the input's abstract linear operator form. - // OMP_NUM_THREADS=4 seems to be the best option for dense sketch generation. - #ifdef RandBLAS_HAS_OpenMP - omp_set_num_threads(this->num_threads_min); - #endif - RandBLAS::DenseDist D(n, k); - state = RandBLAS::fill_dense(D, Y_i, state); - #ifdef RandBLAS_HAS_OpenMP - omp_set_num_threads(this->num_threads_max); - #endif + // Call BK to build Krylov subspaces and band matrices + T* X_ev = nullptr; + T* Y_od = nullptr; + T* R = nullptr; + T* S = nullptr; + int64_t end_rows = 0, end_cols = 0; + bool final_iter_is_odd = false; - if(this -> timing) { - sketching_t_stop = steady_clock::now(); - sketching_t_dur = duration_cast(sketching_t_stop - sketching_t_start).count(); - gemm_A_t_start = steady_clock::now(); - } + int status = bk_obj.call(A, k, X_ev, Y_od, R, S, + end_rows, end_cols, final_iter_is_odd, state); - // [X_ev, ~] = qr(A * Y_i, 0) - A(Side::Left, Layout::ColMajor, Op::NoTrans, Op::NoTrans, m, k, n, 1.0, Y_i, n, 0.0, X_i, m); + // Read back BK outputs + this->num_krylov_iters = bk_obj.num_krylov_iters; + this->norm_R_end = bk_obj.norm_R_end; - if(this -> timing) { - gemm_A_t_stop = steady_clock::now(); - gemm_A_t_dur = duration_cast(gemm_A_t_stop - gemm_A_t_start).count(); - } + if (status != 0) return status; - if(this -> qr_exp == Subroutines::QR_explicit::cqrrt) { - if(this -> timing) - qr_t_start = steady_clock::now(); - - CQRRT -> call(m, k, X_i, m, R_11_trans, k, d_factor, state); + int64_t m = A.n_rows; + int64_t n = A.n_cols; + int increment = (this->adaptive_increment > 0) + ? this->adaptive_increment : this->max_krylov_iters; - if(this -> timing) { - qr_t_stop = steady_clock::now(); - qr_t_dur = duration_cast(qr_t_stop - qr_t_start).count(); - } - } else { + T* U_hat = nullptr; + T* VT_hat = nullptr; + int retries = 0; + // SVD + reconstruction loop (runs once in non-adaptive mode). + while (true) { + // Phase: SVD on band matrix + factor reconstruction if(this -> timing) - qr_t_start = steady_clock::now(); + allocation_t_start = steady_clock::now(); - lapack::geqrf(m, k, X_i, m, tau); + // Internal SVD workspace — freed in this function. + U_hat = ( T * ) malloc( end_rows * end_cols * sizeof( T ) ); + VT_hat = ( T * ) malloc( end_cols * end_cols * sizeof( T ) ); - if(this -> timing) { - qr_t_stop = steady_clock::now(); - qr_t_dur = duration_cast(qr_t_stop - qr_t_start).count(); - ungqr_t_start = steady_clock::now(); - } - - // Convert X_i into an explicit form. It is now stored in X_ev as it should be. - lapack::ungqr(m, k, k, X_i, m, tau); + // Output arrays — ownership transfers to caller (use delete[]). + Sigma = new T[std::min(end_cols, end_rows)](); + U = new T[m * end_cols](); + V = new T[n * end_cols](); if(this -> timing) { - ungqr_t_stop = steady_clock::now(); - ungqr_t_dur += duration_cast(ungqr_t_stop - ungqr_t_start).count(); + allocation_t_stop = steady_clock::now(); + driver_alloc_t_dur += duration_cast(allocation_t_stop - allocation_t_start).count(); + get_factors_t_start = steady_clock::now(); } - } - - // Advance odd iteration count. - ++iter_od; - // Advance iteration count. - ++iter; - - // Iterate until in-loop termination criteria is met. - while(1) { - if(this -> timing) - main_loop_t_start = steady_clock::now(); - - if (iter % 2 != 0) { - if(this -> timing) - gemm_A_t_start = steady_clock::now(); - // Y_i = A' * X_i - A(Side::Left, Layout::ColMajor, Op::Trans, Op::NoTrans, n, k, m, 1.0, X_i, m, 0.0, Y_i, n); - - if(this -> timing) { - gemm_A_t_stop = steady_clock::now(); - gemm_A_t_dur += duration_cast(gemm_A_t_stop - gemm_A_t_start).count(); - allocation_t_start = steady_clock::now(); - } - - // Allocate more space for Y_od - curr_X_cols += k; - X_ev = ( T * ) realloc(X_ev, m * curr_X_cols * sizeof( T )); - // Move the X_i pointer; - X_i = &X_ev[m * (curr_X_cols - k)]; - - - if(this -> timing) { - allocation_t_stop = steady_clock::now(); - allocation_t_dur += duration_cast(allocation_t_stop - allocation_t_start).count(); - reorth_t_start = steady_clock::now(); - } - - if (iter != 1) { - // R_i' = Y_i' * Y_od - blas::gemm(Layout::ColMajor, Op::Trans, Op::NoTrans, k, iter_ev * k, n, 1.0, Y_i, n, Y_od, n, 0.0, R_i, n); - - // Y_i = Y_i - Y_od * R_i - blas::gemm(Layout::ColMajor, Op::NoTrans, Op::Trans, n, k, iter_ev * k, -1.0, Y_od, n, R_i, n, 1.0, Y_i, n); - - // Reorthogonalization - blas::gemm(Layout::ColMajor, Op::Trans, Op::NoTrans, k, iter_ev * k, n, 1.0, Y_i, n, Y_od, n, 0.0, Y_orth_buf, k); - blas::gemm(Layout::ColMajor, Op::NoTrans, Op::Trans, n, k, iter_ev * k, -1.0, Y_od, n, Y_orth_buf, k, 1.0, Y_i, n); - } - if(this -> timing) { - reorth_t_stop = steady_clock::now(); - reorth_t_dur += duration_cast(reorth_t_stop - reorth_t_start).count(); - } - - // Perform explicit QR via a method of choice - if(this -> qr_exp == Subroutines::QR_explicit::cqrrt) { - if(this -> timing) - qr_t_start = steady_clock::now(); - - CQRRT -> call(n, k, Y_i, n, R_11_trans, k, d_factor, state); - // Copy R_ii over to R's (in transposed format). - - util::transposition(0, k, R_11_trans, k, R_ii, n, 1); - if(this -> timing) { - qr_t_stop = steady_clock::now(); - qr_t_dur += duration_cast(qr_t_stop - qr_t_start).count(); - } + if (this->adaptive) { + // Adaptive: run gesdd on a copy to preserve R/S for potential resume. + T* svd_input = ( T * ) malloc( end_rows * end_cols * sizeof( T ) ); + if (final_iter_is_odd) { + lapack::lacpy(MatrixType::General, end_rows, end_cols, R, n, svd_input, end_rows); } else { - // [Y_i, R_ii] = qr(Y_i, 0) - std::fill(&tau[0], &tau[k], 0.0); - - if(this -> timing) - qr_t_start = steady_clock::now(); - lapack::geqrf(n, k, Y_i, n, tau); - - if(this -> timing) { - qr_t_stop = steady_clock::now(); - qr_t_dur += duration_cast(qr_t_stop - qr_t_start).count(); - r_cpy_t_start = steady_clock::now(); - } - - // Copy R_ii over to R's (in transposed format). - #ifdef RandBLAS_HAS_OpenMP - omp_set_num_threads(this->num_threads_min); - #endif - util::transposition(0, k, Y_i, n, R_ii, n, 1); - #ifdef RandBLAS_HAS_OpenMP - omp_set_num_threads(this->num_threads_max); - #endif - - if(this -> timing) { - r_cpy_t_stop = steady_clock::now(); - r_cpy_t_dur += duration_cast(r_cpy_t_stop - r_cpy_t_start).count(); - ungqr_t_start = steady_clock::now(); - } - - // Convert Y_i into an explicit form. It is now stored in Y_odd as it should be. - lapack::ungqr(n, k, k, Y_i, n, tau); - - if(this -> timing) { - ungqr_t_stop = steady_clock::now(); - ungqr_t_dur += duration_cast(ungqr_t_stop - ungqr_t_start).count(); - } - } - - // Early termination - // if (abs(R(end)) <= sqrt(eps('T'))) - if(std::abs(R_ii[(n + 1) * (k - 1)]) < std::sqrt(std::numeric_limits::epsilon())) { - //printf("TERMINATION 1 at iteration %ld\n", iter); - break; - } - - // Allocate more space for R - T* R_new = ( T * ) realloc(R, n * curr_X_cols * sizeof( T )); - if (!R_new) { - // Handle realloc failure. - free(Y_od); - free(X_ev); - free(tau); - free(R); - free(S); - free(U_hat); - free(VT_hat); - free(Y_orth_buf); - free(X_orth_buf); - if(R_11_trans != nullptr) { - free(R_11_trans); - } - return -1; + lapack::lacpy(MatrixType::General, end_rows, end_cols, S, n + k, svd_input, end_rows); } - // Need to make sure the newly-allocated space is empty - R = R_new; - T* temp_r = &R[n * (curr_X_cols - k)]; - std::fill(temp_r, temp_r + n*k, 0.0); - - // Advance R pointers - R_i = &R[(iter_ev + 1) * k]; - R_ii = &R[(n * k * (iter_ev + 1)) + k + (k * (iter_ev))]; - - // Advance even iteration count; - ++iter_ev; - } - else { - if(this -> timing) - gemm_A_t_start = steady_clock::now(); - - // X_i = A * Y_i - A(Side::Left, Layout::ColMajor, Op::NoTrans, Op::NoTrans, m, k, n, 1.0, Y_i, n, 0.0, X_i, m); - - if(this -> timing) { - gemm_A_t_stop = steady_clock::now(); - gemm_A_t_dur += duration_cast(gemm_A_t_stop - gemm_A_t_start).count(); - allocation_t_start = steady_clock::now(); - } - - // Allocate more spece for Y_od - curr_Y_cols += k; - Y_od = ( T * ) realloc(Y_od, n * curr_Y_cols * sizeof( T )); - // Move the X_i pointer; - Y_i = &Y_od[n * (curr_Y_cols - k)]; - - if(this -> timing) { - allocation_t_stop = steady_clock::now(); - allocation_t_dur += duration_cast(allocation_t_stop - allocation_t_start).count(); - reorth_t_start = steady_clock::now(); - } - - // S_i = X_ev' * X_i - blas::gemm(Layout::ColMajor, Op::Trans, Op::NoTrans, iter_od * k, k, m, 1.0, X_ev, m, X_i, m, 0.0, S_i, n + k); - - //X_i = X_i - X_ev * S_i; - blas::gemm(Layout::ColMajor, Op::NoTrans, Op::NoTrans, m, k, iter_od * k, -1.0, X_ev, m, S_i, n + k, 1.0, X_i, m); - - // Reorthogonalization - blas::gemm(Layout::ColMajor, Op::Trans, Op::NoTrans, iter_od * k, k, m, 1.0, X_ev, m, X_i, m, 0.0, X_orth_buf, n + k); - blas::gemm(Layout::ColMajor, Op::NoTrans, Op::NoTrans, m, k, iter_od * k, -1.0, X_ev, m, X_orth_buf, n + k, 1.0, X_i, m); - - if(this -> timing) { - reorth_t_stop = steady_clock::now(); - reorth_t_dur += duration_cast(reorth_t_stop - reorth_t_start).count(); - } - - // Perform explicit QR via a method of choice - if(this -> qr_exp == Subroutines::QR_explicit::cqrrt) { - if(this -> timing) - qr_t_start = steady_clock::now(); - - CQRRT -> call(m, k, X_i, m, S_ii, n + k, d_factor, state); - - if(this -> timing) { - qr_t_stop = steady_clock::now(); - qr_t_dur += duration_cast(qr_t_stop - qr_t_start).count(); - } - + lapack::gesdd(Job::SomeVec, end_rows, end_cols, svd_input, end_rows, + Sigma, U_hat, end_rows, VT_hat, end_cols); + free(svd_input); + } else { + // Non-adaptive: gesdd overwrites R or S directly (they're freed below). + if (final_iter_is_odd) { + lapack::gesdd(Job::SomeVec, end_rows, end_cols, R, n, + Sigma, U_hat, end_rows, VT_hat, end_cols); } else { - // [X_i, S_ii] = qr(X_i, 0); - std::fill(&tau[0], &tau[k], 0.0); - - if(this -> timing) - qr_t_start = steady_clock::now(); - - lapack::geqrf(m, k, X_i, m, tau); - - if(this -> timing) { - qr_t_stop = steady_clock::now(); - qr_t_dur += duration_cast(qr_t_stop - qr_t_start).count(); - s_cpy_t_start = steady_clock::now(); - } - - // Copy S_ii over to S's space under S_i (offset down by iter_od * k) - lapack::lacpy(MatrixType::Upper, k, k, X_i, m, S_ii, n + k); - - if(this -> timing) { - s_cpy_t_stop = steady_clock::now(); - s_cpy_t_dur += duration_cast(s_cpy_t_stop - s_cpy_t_start).count(); - ungqr_t_start = steady_clock::now(); - } - - // Convert X_i into an explicit form. It is now stored in X_ev as it should be - lapack::ungqr(m, k, k, X_i, m, tau); - - if(this -> timing) { - ungqr_t_stop = steady_clock::now(); - ungqr_t_dur += duration_cast(ungqr_t_stop - ungqr_t_start).count(); - } - } - - // Early termination - // if (abs(S(end)) <= sqrt(eps('T'))) - if(std::abs(S_ii[((n + k) + 1) * (k - 1)]) < std::sqrt(std::numeric_limits::epsilon())) { - //printf("TERMINATION 2 at iteration %ld\n", iter); - break; - } - - if(this -> timing) { - allocation_t_start = steady_clock::now(); + lapack::gesdd(Job::SomeVec, end_rows, end_cols, S, n + k, + Sigma, U_hat, end_rows, VT_hat, end_cols); } + } - // Allocate more space for S - T* S_new = ( T * ) realloc(S, (n + k) * curr_Y_cols * sizeof( T )); - if (!S_new) { - // Handle realloc failure. - free(Y_od); - free(X_ev); - free(tau); - free(R); - free(S); - free(U_hat); - free(VT_hat); - free(Y_orth_buf); - free(X_orth_buf); - if(R_11_trans != nullptr) { - free(R_11_trans); - } - return -1; - } - // Need to make sure the newly-allocated space is empty - S = S_new; - T* temp_s = &S[(n + k)* (curr_Y_cols - k)]; - std::fill(temp_s, temp_s + (n + k) * k, 0.0); - - // Advance S pointers - S_i = &S[(n + k) * k * iter_od]; - S_ii = &S[(n + k) * k * iter_od + k + (iter_od * k)]; + // U = X_ev * U_hat + blas::gemm(Layout::ColMajor, Op::NoTrans, Op::NoTrans, m, end_cols, end_rows, + 1.0, X_ev, m, U_hat, end_rows, 0.0, U, m); + // V = Y_od * V_hat + blas::gemm(Layout::ColMajor, Op::NoTrans, Op::Trans, n, end_cols, end_cols, + 1.0, Y_od, n, VT_hat, end_cols, 0.0, V, n); - // Advance odd iteration count; - ++iter_od; + this->singular_triplets_found = end_cols; - if(this -> timing) { - allocation_t_stop = steady_clock::now(); - allocation_t_dur += duration_cast(allocation_t_stop - allocation_t_start).count(); - } + if(this -> timing) { + get_factors_t_stop = steady_clock::now(); + get_factors_t_dur += duration_cast(get_factors_t_stop - get_factors_t_start).count(); } - if(this -> timing) - norm_t_start = steady_clock::now(); + if (!this->adaptive) break; - // This is only changed on odd iters - if (iter % 2 != 0) - norm_R = lapack::lantr(Norm::Fro, Uplo::Upper, Diag::NonUnit, iter_ev * k, iter_ev * k, R, n); + // --- Adaptive residual check --- + T residual = linops::svd_residual(A, U, V, Sigma, end_cols); - if(this -> timing) { - norm_t_stop = steady_clock::now(); - norm_t_dur += duration_cast(norm_t_stop - norm_t_start).count(); - main_loop_t_stop = steady_clock::now(); - main_loop_t_dur += duration_cast(main_loop_t_stop - main_loop_t_start).count(); + if (residual <= this->tol) { + if (this->verbose) + printf("ABRIK adaptive: converged, residual %e <= tol %e after %d retries.\n", + residual, this->tol, retries); + break; } - if (iter >= max_iters) { + if (bk_obj.termination_reason == BKTermination::norm_converged) { + std::cerr << "ABRIK adaptive: BK terminated via norm convergence. " + << "Cannot improve further. Residual = " << residual + << ", tol = " << this->tol << std::endl; break; } - - ++iter; - //norm(R, 'fro') > sqrt(1 - sq_tol) * norm_A - if(norm_R > threshold) { - // Threshold termination. + if (bk_obj.termination_reason == BKTermination::rank_deficient) { + std::cerr << "ABRIK adaptive: BK terminated due to rank deficiency. " + << "Cannot improve further. Residual = " << residual + << ", tol = " << this->tol << std::endl; + break; + } + if (retries >= this->adaptive_max_retries) { + std::cerr << "ABRIK adaptive: reached max retries (" << this->adaptive_max_retries + << "). Residual = " << residual << ", tol = " << this->tol << std::endl; break; } - } - this->norm_R_end = norm_R; - this->num_krylov_iters = iter; - end_cols = num_krylov_iters * k / 2; - iter % 2 == 0 ? end_rows = end_cols + k : end_rows = end_cols; - + // Not satisfied, BK stopped at max_iters: discard current factors, resume BK. + delete[] U; U = nullptr; + delete[] V; V = nullptr; + delete[] Sigma; Sigma = nullptr; + free(U_hat); U_hat = nullptr; + free(VT_hat); VT_hat = nullptr; - if(this -> timing) { - allocation_t_start = steady_clock::now(); - } + bk_obj.max_krylov_iters += increment; + status = bk_obj.resume(A, k, X_ev, Y_od, R, S, + end_rows, end_cols, final_iter_is_odd, state); - U_hat = ( T * ) calloc( end_rows * end_cols, sizeof( T ) ); - VT_hat = ( T * ) calloc( end_cols * end_cols, sizeof( T ) ); + this->num_krylov_iters = bk_obj.num_krylov_iters; + this->norm_R_end = bk_obj.norm_R_end; - if(this -> timing) { - allocation_t_stop = steady_clock::now(); - allocation_t_dur += duration_cast(allocation_t_stop - allocation_t_start).count(); - get_factors_t_start = steady_clock::now(); - } - - - Sigma = new T[std::min(end_cols, end_rows)](); - U = new T[m * end_cols](); - V = new T[n * end_cols](); + if (status != 0) { + // BK resume failed (realloc failure); BK already cleaned up its buffers. + return status; + } - if (iter % 2 != 0) { - // [U_hat, Sigma, V_hat] = svd(R') - lapack::gesdd(Job::SomeVec, end_rows, end_cols, R, n, Sigma, U_hat, end_rows, VT_hat, end_cols); - } else { - // [U_hat, Sigma, V_hat] = svd(S) - lapack::gesdd(Job::SomeVec, end_rows, end_cols, S, n + k, Sigma, U_hat, end_rows, VT_hat, end_cols); + ++retries; } - // U = X_ev * U_hat - blas::gemm(Layout::ColMajor, Op::NoTrans, Op::NoTrans, m, end_cols, end_rows, 1.0, X_ev, m, U_hat, end_rows, 0.0, U, m); - // V = Y_od * V_hat - blas::gemm(Layout::ColMajor, Op::NoTrans, Op::Trans, n, end_cols, end_cols, 1.0, Y_od, n, VT_hat, end_cols, 0.0, V, n); - - this->singular_triplets_found = end_cols; - - if(this -> timing) { - get_factors_t_stop = steady_clock::now(); - get_factors_t_dur = duration_cast(get_factors_t_stop - get_factors_t_start).count(); - allocation_t_start = steady_clock::now(); - } + if(this -> timing) + allocation_t_start = steady_clock::now(); + // Free BK-allocated buffers and SVD workspace free(Y_od); free(X_ev); - free(tau); free(R); free(S); free(U_hat); free(VT_hat); - free(Y_orth_buf); - free(X_orth_buf); - if(R_11_trans != nullptr) { - free(R_11_trans); - } if(this -> timing) { - allocation_t_stop = steady_clock::now(); - allocation_t_dur += duration_cast(allocation_t_stop - allocation_t_start).count(); + allocation_t_stop = steady_clock::now(); + driver_alloc_t_dur += duration_cast(allocation_t_stop - allocation_t_start).count(); } + // Assemble the 13-entry timing vector (same layout as before) 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 - (allocation_t_dur + get_factors_t_dur + ungqr_t_dur + reorth_t_dur + qr_t_dur + gemm_A_t_dur + sketching_t_dur + r_cpy_t_dur + s_cpy_t_dur + norm_t_dur); - this -> times.resize(13); - this -> times = {allocation_t_dur, get_factors_t_dur, ungqr_t_dur, reorth_t_dur, qr_t_dur, gemm_A_t_dur, main_loop_t_dur, sketching_t_dur, r_cpy_t_dur, s_cpy_t_dur, norm_t_dur, t_rest, total_t_dur}; + + // BK times: [0]=alloc, [1]=ungqr, [2]=reorth, [3]=qr, [4]=gemm_A, + // [5]=main_loop, [6]=sketching, [7]=r_cpy, [8]=s_cpy, [9]=norm + auto& bt = bk_obj.times; + long allocation_t_dur = bt[0] + driver_alloc_t_dur; + long ungqr_t_dur = bt[1]; + long reorth_t_dur = bt[2]; + long qr_t_dur = bt[3]; + long gemm_A_t_dur = bt[4]; + long main_loop_t_dur = bt[5]; + long sketching_t_dur = bt[6]; + long r_cpy_t_dur = bt[7]; + long s_cpy_t_dur = bt[8]; + long norm_t_dur = bt[9]; + + long t_rest = total_t_dur - (allocation_t_dur + get_factors_t_dur + ungqr_t_dur + reorth_t_dur + + qr_t_dur + gemm_A_t_dur + sketching_t_dur + r_cpy_t_dur + s_cpy_t_dur + norm_t_dur); + + this -> times = {allocation_t_dur, get_factors_t_dur, ungqr_t_dur, reorth_t_dur, qr_t_dur, + gemm_A_t_dur, main_loop_t_dur, sketching_t_dur, r_cpy_t_dur, s_cpy_t_dur, + norm_t_dur, t_rest, total_t_dur}; if (this -> verbose) { printf("\n\n/------------ABRIK TIMING RESULTS BEGIN------------/\n"); @@ -726,12 +383,15 @@ class ABRIK { printf("S_ii cpy takes %22.2f%% of runtime.\n", 100 * ((T) s_cpy_t_dur / (T) total_t_dur)); printf("Norm R takes %22.2f%% of runtime.\n", 100 * ((T) norm_t_dur / (T) total_t_dur)); printf("Rest takes %22.2f%% of runtime.\n", 100 * ((T) t_rest / (T) total_t_dur)); - + printf("\nMain loop takes %22.2f%% of runtime.\n", 100 * ((T) main_loop_t_dur / (T) total_t_dur)); printf("/-------------ABRIK TIMING RESULTS END-------------/\n\n"); } } return 0; } + + private: + BK bk_obj; }; -} \ No newline at end of file +} diff --git a/RandLAPACK/misc/rl_gen.hh b/RandLAPACK/misc/rl_gen.hh index 714e6433b..5af699a6e 100644 --- a/RandLAPACK/misc/rl_gen.hh +++ b/RandLAPACK/misc/rl_gen.hh @@ -439,18 +439,57 @@ void process_input_mat( // Exit querying mod. workspace_query_mod = 0; } else { - double value; - int i, j; - // Read input file - std::ifstream inputMat(filename); + // Fast matrix read: slurp file into memory, locate row boundaries, + // then parse rows in parallel with strtod + OpenMP. + // + // Why not use ifstream >> ? + // C++ formatted stream extraction (operator>>) is extremely slow for + // large matrices — parsing 100M doubles takes ~30s. This approach: + // 1. fread the entire file into a char buffer (~2s for 800MB, I/O-bound) + // 2. Scan for newlines to find row boundaries (~0.1s, memory-bandwidth) + // 3. Parse each row independently via strtod (~1s with OpenMP) + // Total: ~3s vs ~30s — roughly 10x faster. + // + // Note: input file is row-major text (whitespace-separated doubles, + // one row per line). We store in column-major: A[m * col + row]. + + // Step 1: Read entire file into memory. + FILE* fp = std::fopen(filename, "r"); + if (!fp) + throw std::runtime_error(std::string("Cannot open file: ") + filename); + std::fseek(fp, 0, SEEK_END); + long file_size = std::ftell(fp); + std::fseek(fp, 0, SEEK_SET); + + std::vector buf(file_size + 1); + std::fread(buf.data(), 1, file_size, fp); + buf[file_size] = '\0'; + std::fclose(fp); + + // Step 2: Find the starting byte offset of each row. + // Each row ends with '\n'; row 0 starts at byte 0. + // Text numbers have variable width, so we must scan for newlines + // (there is no closed-form formula for row offsets in a text file). + std::vector row_starts(m); + row_starts[0] = buf.data(); + int64_t row_idx = 1; + for (long pos = 0; pos < file_size && row_idx < m; ++pos) { + if (buf[pos] == '\n') { + row_starts[row_idx++] = buf.data() + pos + 1; + } + } - // Place the contents of a file into the matrix space. - // Matrix is input in a row-major order, we process data in column-major. - // Reads here are, unfortunately, sequential; - for(j = 0; j < m; ++j) { - for(i = 0; i < n; ++i) { - inputMat >> value; - A[m * i + j] = value; + // Step 3: Parse each row in parallel. + // strtod is thread-safe (no shared state). Each thread handles a + // contiguous block of rows and writes to non-overlapping columns of A + // (column-major storage: row j of input → A[m*i + j] for each col i). + #pragma omp parallel for schedule(static) + for (int64_t j = 0; j < m; ++j) { + char* ptr = row_starts[j]; + char* end; + for (int64_t i = 0; i < n; ++i) { + A[m * i + j] = (T) std::strtod(ptr, &end); + ptr = end; } } } diff --git a/RandLAPACK/misc/rl_util.hh b/RandLAPACK/util/rl_util.hh similarity index 100% rename from RandLAPACK/misc/rl_util.hh rename to RandLAPACK/util/rl_util.hh diff --git a/RandLAPACK/util/rl_util_linop.hh b/RandLAPACK/util/rl_util_linop.hh new file mode 100644 index 000000000..07fd6f450 --- /dev/null +++ b/RandLAPACK/util/rl_util_linop.hh @@ -0,0 +1,43 @@ +#pragma once + +#include "rl_linops.hh" + +namespace RandLAPACK::linops { + +/// Computes the SVD residual: +/// sqrt(||AV - U * diag(Sigma)||^2_F + ||A'U - V * diag(Sigma)||^2_F) / sigma_k. +/// U is m x k (col-major, ld m), V is n x k (col-major, ld n), Sigma is length k. +template +T svd_residual(GLO& A, T* U, T* V, T* Sigma, int64_t k) { + int64_t m = A.n_rows; + int64_t n = A.n_cols; + + T* U_cpy = new T[m * k](); + T* V_cpy = new T[n * k](); + + // U_cpy = U * diag(Sigma) + lapack::lacpy(MatrixType::General, m, k, U, m, U_cpy, m); + for (int64_t i = 0; i < k; ++i) + blas::scal(m, Sigma[i], &U_cpy[m * i], 1); + + // U_cpy = AV - U*diag(Sigma) + A(Layout::ColMajor, Op::NoTrans, Op::NoTrans, m, k, n, (T)1.0, V, n, (T)-1.0, U_cpy, m); + + // V_cpy = V * diag(Sigma) + lapack::lacpy(MatrixType::General, n, k, V, n, V_cpy, n); + for (int64_t i = 0; i < k; ++i) + blas::scal(n, Sigma[i], &V_cpy[n * i], 1); + + // V_cpy = A'U - V*diag(Sigma) + A(Layout::ColMajor, Op::Trans, Op::NoTrans, n, k, m, (T)1.0, U, m, (T)-1.0, V_cpy, n); + + T nrm1 = lapack::lange(Norm::Fro, m, k, U_cpy, m); + T nrm2 = lapack::lange(Norm::Fro, n, k, V_cpy, n); + + delete[] U_cpy; + delete[] V_cpy; + + return std::hypot(nrm1, nrm2) / Sigma[k - 1]; +} + +} // end namespace RandLAPACK::linops diff --git a/benchmark/CMakeLists.txt b/benchmark/CMakeLists.txt index b058c06af..eda4b2c73 100644 --- a/benchmark/CMakeLists.txt +++ b/benchmark/CMakeLists.txt @@ -75,7 +75,7 @@ include(FetchContent) FetchContent_Declare( eigen_lib GIT_REPOSITORY https://gitlab.com/libeigen/eigen - GIT_TAG master + GIT_TAG 3.4.1 GIT_SHALLOW TRUE ) FetchContent_MakeAvailable(eigen_lib) diff --git a/benchmark/bench_ABRIK/ABRIK_runtime_breakdown.cc b/benchmark/bench_ABRIK/ABRIK_runtime_breakdown.cc index 4a266d736..16dba05d1 100644 --- a/benchmark/bench_ABRIK/ABRIK_runtime_breakdown.cc +++ b/benchmark/bench_ABRIK/ABRIK_runtime_breakdown.cc @@ -1,18 +1,17 @@ /* ABRIK runtime breakdown benchmark - assesses the time taken by each subcomponent of ABRIK. -Records all, data, not just the best. -There are 10 things that we time: - 1.Allocate and free time. - 2.Time to acquire the SVD factors. - 3.UNGQR time. - 4.Reorthogonalization time. - 5.QR time. - 6.GEMM A time. - 7.Sketching time. - 8.R_ii cpy time. - 9.S_ii cpy time. - 10.Norm R time. +Precision (float or double) is specified as the first CLI argument. +Records all data, not just the best. + +Output: CSV file with '#'-prefixed metadata header, column names, then data rows. +ABRIK allocates U/V/Sigma with new[] internally -> cleanup with delete[]. + +Subroutine timings (microseconds): + 1. Allocation/free 2. SVD factors 3. UNGQR 4. Reorthogonalization + 5. QR 6. GEMM A 7. Main loop 8. Sketching + 9. R_ii copy 10. S_ii copy 11. Norm R 12. Rest 13. Total */ + #include "RandLAPACK.hh" #include "rl_blaspp.hh" #include "rl_lapackpp.hh" @@ -20,6 +19,7 @@ There are 10 things that we time: #include #include +#include using Subroutines = RandLAPACK::ABRIKSubroutines; @@ -44,168 +44,164 @@ struct ABRIK_benchmark_data { Sigma = nullptr; } - ~ABRIK_benchmark_data(){ + ~ABRIK_benchmark_data() { delete[] A; - delete[] U; - delete[] V; - delete[] Sigma; } }; -// Re-generate and clear data -template -static void data_regen(RandLAPACK::gen::mat_gen_info m_info, - ABRIK_benchmark_data &all_data, - RandBLAS::RNGState &state, int overwrite_A) { - auto m = all_data.row; - auto n = all_data. col; - - if (overwrite_A) { - RandLAPACK::gen::mat_gen(m_info, all_data.A, state); - } - delete[] all_data.U; - delete[] all_data.V; - delete[] all_data.Sigma; - all_data.U = nullptr; - all_data.V = nullptr; - all_data.Sigma = nullptr; -} - template static void call_all_algs( - RandLAPACK::gen::mat_gen_info m_info, int64_t num_runs, int64_t k, int64_t num_matmuls, ABRIK_benchmark_data &all_data, RandBLAS::RNGState &state, - std::string output_filename) { + std::ofstream &outfile) { auto m = all_data.row; - auto n = all_data. col; + auto n = all_data.col; auto tol = all_data.tolerance; bool time_subroutines = true; // Additional params setup. - RandLAPACK::ABRIK ABRIK(false, time_subroutines, tol); + RandLAPACK::ABRIK ABRIK(false, time_subroutines, tol); ABRIK.max_krylov_iters = num_matmuls; - ABRIK.num_threads_min = 4; - //ABRIK.qr_exp = Subroutines::QR_explicit::cqrrt; - ABRIK.num_threads_max = RandLAPACK::util::get_omp_threads(); - // Making sure the states are unchanged - auto state_gen = state; auto state_alg = state; - - // Timing vars std::vector inner_timing; for (int i = 0; i < num_runs; ++i) { printf("\nBlock size %ld, num matmuls %ld. Iteration %d start.\n", k, num_matmuls, i); ABRIK.call(m, n, all_data.A, m, k, all_data.U, all_data.V, all_data.Sigma, state_alg); - - // Update timing vector - inner_timing = ABRIK.times; - // Add info about the run - inner_timing.insert (inner_timing.begin(), num_matmuls); - inner_timing.insert (inner_timing.begin(), k); - - std::ofstream file(output_filename, std::ios::app); - std::copy(inner_timing.begin(), inner_timing.end(), std::ostream_iterator(file, ", ")); - file << "\n"; - - // Clear and re-generate data - data_regen(m_info, all_data, state_gen, 0); - state_gen = state; + + // Write CSV data row: b_sz, num_matmuls, then timing columns + outfile << k << ", " << num_matmuls; + for (const auto &t : ABRIK.times) + outfile << ", " << t; + outfile << "\n"; + outfile.flush(); + + // Cleanup ABRIK outputs (new[]) + delete[] all_data.U; all_data.U = nullptr; + delete[] all_data.V; all_data.V = nullptr; + delete[] all_data.Sigma; all_data.Sigma = nullptr; + state_alg = state; } } -int main(int argc, char *argv[]) { +template +static void run_benchmark(int argc, char *argv[]) { - if (argc < 11) { - // Expected input into this benchmark. - std::cerr << "Usage: " << argv[0] << " " << std::endl; - return 1; + if (argc < 12) { + std::cerr << "Usage: " << argv[0] << " " << std::endl; + return; } - int num_runs = std::stol(argv[3]); - int64_t m_expected = std::stol(argv[4]); - int64_t n_expected = std::stol(argv[5]); - int64_t custom_rank = std::stol(argv[6]); + int num_runs = std::stol(argv[4]); + int64_t m_expected = std::stol(argv[5]); + int64_t n_expected = std::stol(argv[6]); + int64_t custom_rank = std::stol(argv[7]); std::vector b_sz; - for (int i = 0; i < std::stol(argv[7]); ++i) - b_sz.push_back(std::stoi(argv[i + 9])); - // Save elements in string for logging purposes + for (int i = 0; i < std::stol(argv[8]); ++i) + b_sz.push_back(std::stoi(argv[i + 10])); std::ostringstream oss1; for (const auto &val : b_sz) oss1 << val << ", "; std::string b_sz_string = oss1.str(); std::vector matmuls; - for (int i = 0; i < std::stol(argv[8]); ++i) - matmuls.push_back(std::stoi(argv[i + 9 + std::stol(argv[7])])); - // Save elements in string for logging purposes + for (int i = 0; i < std::stol(argv[9]); ++i) + matmuls.push_back(std::stoi(argv[i + 10 + std::stol(argv[8])])); std::ostringstream oss2; for (const auto &val : matmuls) oss2 << val << ", "; std::string matmuls_string = oss2.str(); - double tol = std::pow(std::numeric_limits::epsilon(), 0.85); + T tol = std::pow(std::numeric_limits::epsilon(), (T)0.85); auto state = RandBLAS::RNGState(); auto state_constant = state; int64_t m = 0, n = 0; // Generate the input matrix. - RandLAPACK::gen::mat_gen_info m_info(m, n, RandLAPACK::gen::custom_input); - m_info.filename = argv[2]; + RandLAPACK::gen::mat_gen_info m_info(m, n, RandLAPACK::gen::custom_input); + m_info.filename = argv[3]; m_info.workspace_query_mod = 1; - // Workspace query; - RandLAPACK::gen::mat_gen(m_info, NULL, state); + RandLAPACK::gen::mat_gen(m_info, NULL, state); // Update basic params. m = m_info.rows; n = m_info.cols; if (m_expected != m || n_expected != n) { - printf("Expected m: %ld, actual m: %ld\n", m_expected, m); - printf("Expected n: %ld, actual n: %ld\n", n_expected, n); - std::cerr << "Expected input size did not matrch actual input size. Aborting." << std::endl; - return 1; + std::cerr << "Expected input size (" << m_expected << ", " << n_expected << ") did not match actual input size (" << m << ", " << n << "). Aborting." << std::endl; + return; } // Allocate basic workspace. - ABRIK_benchmark_data all_data(m, n, tol); - - // Fill the data matrix; + ABRIK_benchmark_data all_data(m, n, tol); RandLAPACK::gen::mat_gen(m_info, all_data.A, state); printf("Finished data preparation\n"); - // Declare a data file - std::string output_filename = "_ABRIK_runtime_breakdown_num_info_lines_" + std::to_string(6) + ".txt"; + + // Generate date/time prefix + std::time_t now = std::time(nullptr); + char date_prefix[20]; + std::strftime(date_prefix, sizeof(date_prefix), "%Y%m%d_%H%M%S_", std::localtime(&now)); + + // Build output file path + std::string output_filename = std::string(date_prefix) + "ABRIK_runtime_breakdown.csv"; std::string path; - if (std::string(argv[1]) != ".") { - path = argv[1] + output_filename; + if (std::string(argv[2]) != ".") { + path = std::string(argv[2]) + "/" + output_filename; } else { path = output_filename; } std::ofstream file(path, std::ios::out | std::ios::app); - // Writing important data into file - file << "Description: Results from the ABRIK runtime breakdown benchmark, recording the time it takes to perform every subroutine in ABRIK." - "\nFile format: 13 data columns, each corresponding to a given ABRIK subroutine: allocation_t_dur, get_factors_t_dur, ungqr_t_dur, reorth_t_dur, qr_t_dur, gemm_A_t_dur, main_loop_t_dur, sketching_t_dur, r_cpy_t_dur, s_cpy_t_dur, norm_t_dur, t_rest, total_t_dur" - " rows correspond to ABRIK runs with block sizes varying as specified, with numruns repititions of each block size" - "\nInput type:" + std::string(argv[2]) + - "\nInput size:" + std::to_string(m) + " by " + std::to_string(n) + - "\nAdditional parameters: Krylov block sizes " + b_sz_string + - " matmuls: " + matmuls_string + - " num runs per size " + std::to_string(num_runs) + - " num singular values and vectors approximated " + std::to_string(custom_rank) + - "\n"; + // Write metadata header (prefixed with # for easy parsing) + file << "# ABRIK Runtime Breakdown Benchmark\n" + << "# Precision: " << argv[1] << "\n" + << "# Input matrix: " << argv[3] << "\n" + << "# Input size: " << m << " x " << n << "\n" + << "# Target rank: " << custom_rank << "\n" + << "# Krylov block sizes: " << b_sz_string << "\n" + << "# Matmul counts: " << matmuls_string << "\n" + << "# Runs per configuration: " << num_runs << "\n" + << "# Tolerance: " << tol << "\n" + << "# Timings in microseconds\n"; + // Write CSV column header + file << "b_sz, num_matmuls, " + << "allocation_t, get_factors_t, ungqr_t, reorth_t, qr_t, gemm_A_t, " + << "main_loop_t, sketching_t, r_cpy_t, s_cpy_t, norm_t, t_rest, total_t\n"; file.flush(); + auto start_total = steady_clock::now(); + size_t i = 0, j = 0; - for (;i < b_sz.size(); ++i) { - for (;j < matmuls.size(); ++j) { - call_all_algs(m_info, num_runs, b_sz[i], matmuls[j], all_data, state_constant, path); + for (; i < b_sz.size(); ++i) { + for (; j < matmuls.size(); ++j) { + call_all_algs(num_runs, b_sz[i], matmuls[j], all_data, state_constant, file); } j = 0; } + + auto stop_total = steady_clock::now(); + long total_us = duration_cast(stop_total - start_total).count(); + printf("\nTOTAL BENCHMARK TIME: %.2f seconds\n", total_us / 1e6); +} + +int main(int argc, char *argv[]) { + if (argc < 2) { + std::cerr << "Usage: " << argv[0] << " " << std::endl; + return 1; + } + + std::string precision = argv[1]; + if (precision == "double") { + run_benchmark(argc, argv); + } else if (precision == "float") { + run_benchmark(argc, argv); + } else { + std::cerr << "Error: precision must be 'double' or 'float', got '" << precision << "'" << std::endl; + return 1; + } + return 0; } diff --git a/benchmark/bench_ABRIK/ABRIK_runtime_breakdown_sparse.cc b/benchmark/bench_ABRIK/ABRIK_runtime_breakdown_sparse.cc index 627383865..64ee423c4 100644 --- a/benchmark/bench_ABRIK/ABRIK_runtime_breakdown_sparse.cc +++ b/benchmark/bench_ABRIK/ABRIK_runtime_breakdown_sparse.cc @@ -1,18 +1,17 @@ /* -ABRIK runtime breakdown benchmark - assesses the time taken by each subcomponent of ABRIK. -Records all, data, not just the best. -There are 10 things that we time: - 1.Allocate and free time. - 2.Time to acquire the SVD factors. - 3.UNGQR time. - 4.Reorthogonalization time. - 5.QR time. - 6.GEMM A time. - 7.Sketching time. - 8.R_ii cpy time. - 9.S_ii cpy time. - 10.Norm R time. +ABRIK runtime breakdown benchmark (sparse) - assesses the time taken by each subcomponent of ABRIK +on sparse input matrices read from Matrix Market format. +Precision (float or double) is specified as the first CLI argument. + +Output: CSV file with '#'-prefixed metadata header, column names, then data rows. +ABRIK allocates U/V/Sigma with new[] internally -> cleanup with delete[]. + +Subroutine timings (microseconds): + 1. Allocation/free 2. SVD factors 3. UNGQR 4. Reorthogonalization + 5. QR 6. GEMM A 7. Main loop 8. Sketching + 9. R_ii copy 10. S_ii copy 11. Norm R 12. Rest 13. Total */ + #include "RandLAPACK.hh" #include "rl_blaspp.hh" #include "rl_lapackpp.hh" @@ -20,6 +19,7 @@ There are 10 things that we time: #include #include +#include #include @@ -45,16 +45,11 @@ struct ABRIK_benchmark_data { Sigma = nullptr; } - ~ABRIK_benchmark_data(){ - delete[] U; - delete[] V; - delete[] Sigma; - } + ~ABRIK_benchmark_data() {} }; template RandBLAS::sparse_data::coo::COOMatrix from_matrix_market(std::string fn) { - int64_t n_rows, n_cols = 0; std::vector rows{}; std::vector cols{}; @@ -62,11 +57,11 @@ RandBLAS::sparse_data::coo::COOMatrix from_matrix_market(std::string fn) { std::ifstream file_stream(fn); fast_matrix_market::read_matrix_market_triplet( - file_stream, n_rows, n_cols, rows, cols, vals + file_stream, n_rows, n_cols, rows, cols, vals ); RandBLAS::sparse_data::coo::COOMatrix out(n_rows, n_cols); - reserve_coo(vals.size(),out); + reserve_coo(vals.size(), out); for (int i = 0; i < out.nnz; ++i) { out.rows[i] = rows[i]; out.cols[i] = cols[i]; @@ -79,18 +74,16 @@ RandBLAS::sparse_data::coo::COOMatrix from_matrix_market(std::string fn) { template RandBLAS::CSCMatrix format_conversion(int64_t m, int64_t n, RandBLAS::COOMatrix& input_mat_coo) { - // Grab the leading principal submatrix of the size of half the input - RandBLAS::COOMatrix input_mat_transformed(m, n); + // Grab the leading principal submatrix + RandBLAS::COOMatrix input_mat_transformed(m, n); - // check how many nonzeros are in the left principal submatrix + // Count nonzeros in the submatrix int64_t nnz_sub = 0; for (int64_t i = 0; i < input_mat_coo.nnz; ++i) { - if (input_mat_coo.rows[i] < m && input_mat_coo.cols[i] < n) { + if (input_mat_coo.rows[i] < m && input_mat_coo.cols[i] < n) ++nnz_sub; - } } - // Allocate reserve_coo(nnz_sub, input_mat_transformed); int64_t ell = 0; @@ -103,25 +96,12 @@ RandBLAS::CSCMatrix format_conversion(int64_t m, int64_t n, RandBLAS::COOMatr } } - // Convert the sparse matrix format for performance - RandBLAS::CSCMatrix input_mat_csc(m, n); - //RandBLAS::conversions::coo_to_csc(input_mat_coo, input_mat_csc); + RandBLAS::CSCMatrix input_mat_csc(m, n); RandBLAS::conversions::coo_to_csc(input_mat_transformed, input_mat_csc); return input_mat_csc; } -// Re-generate and clear data -template -static void data_regen(ABRIK_benchmark_data &all_data) { - delete[] all_data.U; - delete[] all_data.V; - delete[] all_data.Sigma; - all_data.U = nullptr; - all_data.V = nullptr; - all_data.Sigma = nullptr; -} - template static void call_all_algs( int64_t num_runs, @@ -129,126 +109,149 @@ static void call_all_algs( int64_t num_matmuls, ABRIK_benchmark_data &all_data, RandBLAS::RNGState &state, - std::string output_filename, + std::ofstream &outfile, int use_cqrrt) { auto m = all_data.row; - auto n = all_data. col; + auto n = all_data.col; auto tol = all_data.tolerance; bool time_subroutines = true; // Additional params setup. - RandLAPACK::ABRIK ABRIK(false, time_subroutines, tol); + RandLAPACK::ABRIK ABRIK(false, time_subroutines, tol); ABRIK.max_krylov_iters = num_matmuls; - ABRIK.num_threads_min = 4; if (use_cqrrt) ABRIK.qr_exp = Subroutines::QR_explicit::cqrrt; - ABRIK.num_threads_max = RandLAPACK::util::get_omp_threads(); - // Making sure the states are unchanged auto state_alg = state; - // Timing vars - std::vector inner_timing; - for (int i = 0; i < num_runs; ++i) { printf("\nBlock size %ld, num matmuls %ld. Iteration %d start.\n", k, num_matmuls, i); ABRIK.call(m, n, *all_data.A_input, m, k, all_data.U, all_data.V, all_data.Sigma, state_alg); - - // Update timing vector - inner_timing = ABRIK.times; - // Add info about the run - inner_timing.insert (inner_timing.begin(), num_matmuls); - inner_timing.insert (inner_timing.begin(), k); - - std::ofstream file(output_filename, std::ios::app); - std::copy(inner_timing.begin(), inner_timing.end(), std::ostream_iterator(file, ", ")); - file << "\n"; - - // Clear and re-generate data - data_regen(all_data); + + // Write CSV data row: b_sz, num_matmuls, then timing columns + outfile << k << ", " << num_matmuls; + for (const auto &t : ABRIK.times) + outfile << ", " << t; + outfile << "\n"; + outfile.flush(); + + // Cleanup ABRIK outputs (new[]) + delete[] all_data.U; all_data.U = nullptr; + delete[] all_data.V; all_data.V = nullptr; + delete[] all_data.Sigma; all_data.Sigma = nullptr; + state_alg = state; } } -int main(int argc, char *argv[]) { -/* +template +static void run_benchmark(int argc, char *argv[]) { + if (argc < 9) { - // Expected input into this benchmark. - std::cerr << "Usage: " << argv[0] << " " << std::endl; - return 1; + std::cerr << "Usage: " << argv[0] << " " << std::endl; + return; } -*/ - double submatrix_dim_ratio = 0.5; + int num_runs = std::stol(argv[4]); + int64_t custom_rank = std::stol(argv[5]); + T submatrix_dim_ratio = (T)std::stod(argv[6]); - int num_runs = std::stol(argv[3]); - int64_t custom_rank = std::stol(argv[4]); std::vector b_sz; - for (int i = 0; i < std::stol(argv[5]); ++i) - b_sz.push_back(std::stoi(argv[i + 7])); - // Save elements in string for logging purposes + for (int i = 0; i < std::stol(argv[7]); ++i) + b_sz.push_back(std::stoi(argv[i + 9])); std::ostringstream oss1; for (const auto &val : b_sz) oss1 << val << ", "; std::string b_sz_string = oss1.str(); std::vector matmuls; - for (int i = 0; i < std::stol(argv[6]); ++i) - matmuls.push_back(std::stoi(argv[i + 7 + std::stol(argv[5])])); - // Save elements in string for logging purposes + for (int i = 0; i < std::stol(argv[8]); ++i) + matmuls.push_back(std::stoi(argv[i + 9 + std::stol(argv[7])])); std::ostringstream oss2; for (const auto &val : matmuls) oss2 << val << ", "; std::string matmuls_string = oss2.str(); - double tol = std::pow(std::numeric_limits::epsilon(), 0.85); + T tol = std::pow(std::numeric_limits::epsilon(), (T)0.85); auto state = RandBLAS::RNGState(); auto state_constant = state; - // Read the input fast matrix market data - // The idea is that input_mat_coo will be automatically freed at the end of function execution - auto input_mat_coo = from_matrix_market(std::string(argv[2])); - auto m = input_mat_coo.n_rows * submatrix_dim_ratio; - auto n = input_mat_coo.n_cols * submatrix_dim_ratio; + // Read the input matrix market data + auto input_mat_coo = from_matrix_market(std::string(argv[3])); + auto m = (int64_t)(input_mat_coo.n_rows * submatrix_dim_ratio); + auto n = (int64_t)(input_mat_coo.n_cols * submatrix_dim_ratio); - // Convert coo into csc matrix - this will grab the leading principal submatrix - // depending on what m and n were set to. - auto input_mat_csc = format_conversion(m, n, input_mat_coo); + // Convert COO to CSC (grabs leading principal submatrix) + auto input_mat_csc = format_conversion(m, n, input_mat_coo); // Allocate basic workspace. - ABRIK_benchmark_data> all_data(m, n, tol); + ABRIK_benchmark_data> all_data(m, n, tol); all_data.A_input = &input_mat_csc; printf("Finished data preparation\n"); - // Declare a data file - std::string output_filename = "_ABRIK_runtime_breakdown_sparse_num_info_lines_" + std::to_string(6) + ".txt"; + + // Generate date/time prefix + std::time_t now = std::time(nullptr); + char date_prefix[20]; + std::strftime(date_prefix, sizeof(date_prefix), "%Y%m%d_%H%M%S_", std::localtime(&now)); + + // Build output file path + std::string output_filename = std::string(date_prefix) + "ABRIK_runtime_breakdown_sparse.csv"; std::string path; - if (std::string(argv[1]) != ".") { - path = argv[1] + output_filename; + if (std::string(argv[2]) != ".") { + path = std::string(argv[2]) + "/" + output_filename; } else { path = output_filename; } std::ofstream file(path, std::ios::out | std::ios::app); - // Writing important data into file - file << "Description: Results from the saprse ABRIK runtime breakdown benchmark, recording the time it takes to perform every subroutine in ABRIK." - "\nFile format: 13 data columns, each corresponding to a given ABRIK subroutine: allocation_t_dur, get_factors_t_dur, ungqr_t_dur, reorth_t_dur, qr_t_dur, gemm_A_t_dur, main_loop_t_dur, sketching_t_dur, r_cpy_t_dur, s_cpy_t_dur, norm_t_dur, t_rest, total_t_dur" - " rows correspond to ABRIK runs with block sizes varying as specified, with numruns repititions of each block size" - "\nInput type:" + std::string(argv[2]) + - "\nInput size:" + std::to_string(m) + " by " + std::to_string(n) + - "\nAdditional parameters: Krylov block sizes " + b_sz_string + - " matmuls: " + matmuls_string + - " num runs per size " + std::to_string(num_runs) + - " num singular values and vectors approximated " + std::to_string(custom_rank) + - "\n"; + // Write metadata header (prefixed with # for easy parsing) + file << "# ABRIK Runtime Breakdown Benchmark (Sparse)\n" + << "# Precision: " << argv[1] << "\n" + << "# Input matrix: " << argv[3] << "\n" + << "# Input size: " << m << " x " << n << " (submatrix ratio: " << submatrix_dim_ratio << ")\n" + << "# Target rank: " << custom_rank << "\n" + << "# Krylov block sizes: " << b_sz_string << "\n" + << "# Matmul counts: " << matmuls_string << "\n" + << "# Runs per configuration: " << num_runs << "\n" + << "# Tolerance: " << tol << "\n" + << "# Timings in microseconds\n"; + // Write CSV column header + file << "b_sz, num_matmuls, " + << "allocation_t, get_factors_t, ungqr_t, reorth_t, qr_t, gemm_A_t, " + << "main_loop_t, sketching_t, r_cpy_t, s_cpy_t, norm_t, t_rest, total_t\n"; file.flush(); + auto start_total = steady_clock::now(); + int use_cqrrt = 1; size_t i = 0, j = 0; - for (;i < b_sz.size(); ++i) { - for (;j < matmuls.size(); ++j) { - call_all_algs(num_runs, b_sz[i], matmuls[j], all_data, state_constant, path, use_cqrrt); + for (; i < b_sz.size(); ++i) { + for (; j < matmuls.size(); ++j) { + call_all_algs(num_runs, b_sz[i], matmuls[j], all_data, state_constant, file, use_cqrrt); } j = 0; } + + auto stop_total = steady_clock::now(); + long total_us = duration_cast(stop_total - start_total).count(); + printf("\nTOTAL BENCHMARK TIME: %.2f seconds\n", total_us / 1e6); +} + +int main(int argc, char *argv[]) { + if (argc < 2) { + std::cerr << "Usage: " << argv[0] << " " << std::endl; + return 1; + } + + std::string precision = argv[1]; + if (precision == "double") { + run_benchmark(argc, argv); + } else if (precision == "float") { + run_benchmark(argc, argv); + } else { + std::cerr << "Error: precision must be 'double' or 'float', got '" << precision << "'" << std::endl; + return 1; + } + return 0; } diff --git a/benchmark/bench_ABRIK/ABRIK_speed_comparisons.cc b/benchmark/bench_ABRIK/ABRIK_speed_comparisons.cc index 01ad32fca..3d32adb67 100644 --- a/benchmark/bench_ABRIK/ABRIK_speed_comparisons.cc +++ b/benchmark/bench_ABRIK/ABRIK_speed_comparisons.cc @@ -1,9 +1,16 @@ /* -Additional ABRIK speed comparison benchmark - runs ABRIK, RSVD and SVDS from Spectra library. -The user is required to provide a matrix file to be read, set min and max numbers of large gemms (Krylov iterations) that the algorithm is allowed to perform min and max block sizes that ABRIK is to use; -furthermore, the user is to provide a 'custom rank' parameter (number of singular vectors to approximate by ABRIK). -The benchmark outputs the basic data of a given run, as well as the ABRIK runtime and singular vector residual error, -which is computed as "sqrt(||AV - SU||^2_F + ||A'U - VS||^2_F / sqrt(target_rank)" (for "custom rank" singular vectors and values). +ABRIK speed comparison benchmark - runs ABRIK, RSVD, SVDS (Spectra), and optionally full SVD (GESDD). +Precision (float or double) is specified as the first CLI argument. +The user provides a matrix file, numbers of Krylov iterations, block sizes, and a target rank. + +Output: CSV file with '#'-prefixed metadata header, column names, then data rows. +Each algorithm section manages its own output allocation/deallocation: + - ABRIK allocates with new[] internally -> cleanup with delete[] + - RSVD allocates with calloc internally -> cleanup with free() + - SVDS/SVD are allocated by the benchmark -> cleanup with delete[] + +Residual metric: sqrt(||S^{-1}AV - U||^2_F + ||(A'U)S^{-1} - V||^2_F). +Timings in microseconds. */ #include "RandLAPACK.hh" @@ -14,12 +21,24 @@ which is computed as "sqrt(||AV - SU||^2_F + ||A'U - VS||^2_F / sqrt(target_rank #include #include #include +#include +#include // External libs includes #include #include -using Matrix = Eigen::MatrixXf; -using Vector = Eigen::VectorXf; +#include "BudgetedSVDSolver.hh" + +// Traits struct mapping scalar type T to Eigen matrix/vector types. +template struct EigenTypes; +template <> struct EigenTypes { + using Matrix = Eigen::MatrixXd; + using Vector = Eigen::VectorXd; +}; +template <> struct EigenTypes { + using Matrix = Eigen::MatrixXf; + using Vector = Eigen::VectorXf; +}; template struct ABRIK_benchmark_data { @@ -28,45 +47,23 @@ struct ABRIK_benchmark_data { T tolerance; T* A; T* U; - T* VT; - T* V; + T* VT; + T* V; T* Sigma; - T* A_lowrank_svd; - T* A_lowrank_svd_const; - T* Buffer; - T* Sigma_cpy; - T* U_cpy; - T* V_cpy; - Matrix A_spectra; - - ABRIK_benchmark_data(int64_t m, int64_t n, T tol) : - A_spectra(m, n) + ABRIK_benchmark_data(int64_t m, int64_t n, T tol) { - A = new T[m * n](); - U = nullptr; - VT = nullptr; - V = nullptr; - Sigma = nullptr; - U_cpy = nullptr; - V_cpy = nullptr; - - A_lowrank_svd = nullptr; - A_lowrank_svd_const = nullptr; - row = m; - col = n; - tolerance = tol; + A = new T[m * n](); + U = nullptr; + VT = nullptr; + V = nullptr; + Sigma = nullptr; + row = m; + col = n; + tolerance = tol; } ~ABRIK_benchmark_data() { delete[] A; - delete[] U; - delete[] VT; - delete[] V; - delete[] Sigma; - delete[] U_cpy; - delete[] V_cpy; - delete[] A_lowrank_svd; - delete[] A_lowrank_svd_const; } }; @@ -82,12 +79,12 @@ struct ABRIK_algorithm_objects { RandLAPACK::ABRIK ABRIK; ABRIK_algorithm_objects( - bool verbosity, - bool cond_check, - bool orth_check, - bool time_subroutines, - int64_t p, - int64_t passes_per_iteration, + bool verbosity, + bool cond_check, + bool orth_check, + bool time_subroutines, + int64_t p, + int64_t passes_per_iteration, int64_t block_sz, T tol ) : @@ -102,95 +99,42 @@ struct ABRIK_algorithm_objects { {} }; -// Re-generate and clear data +// Re-generate input matrix A (re-reads from file). Only needed after GESDD destroys A. template -static void data_regen(RandLAPACK::gen::mat_gen_info m_info, - ABRIK_benchmark_data &all_data, - RandBLAS::RNGState &state, int overwrite_A) { - - auto m = all_data.row; - auto n = all_data.col; - - if (overwrite_A) { - RandLAPACK::gen::mat_gen(m_info, all_data.A, state); - Eigen::Map(all_data.A_spectra.data(), all_data.A_spectra.rows(), all_data.A_spectra.cols()) = Eigen::Map(all_data.A, m, n); - if (all_data.A_lowrank_svd != nullptr) - lapack::lacpy(MatrixType::General, m, n, all_data.A_lowrank_svd_const, m, all_data.A_lowrank_svd, m); - } - - delete[] all_data.U; - delete[] all_data.VT; - delete[] all_data.V; - delete[] all_data.Sigma; - delete[] all_data.U_cpy; - delete[] all_data.V_cpy; - - all_data.U = nullptr; - all_data.VT = nullptr; - all_data.V = nullptr; - all_data.Sigma = nullptr; - all_data.U_cpy = nullptr; - all_data.V_cpy = nullptr; +static void regen_input(RandLAPACK::gen::mat_gen_info m_info, + ABRIK_benchmark_data &all_data, + RandBLAS::RNGState state) { + RandLAPACK::gen::mat_gen(m_info, all_data.A, state); } -// This routine computes the residual norm error, consisting of two parts (one of which) vanishes -// in exact precision. Target_rank defines size of U, V as returned by ABRIK; target_rank <= target_rank. -template +// Computes the residual norm error: sqrt(||S^{-1}AV - U||^2_F + ||(A'U)S^{-1} - V||^2_F). +// Scratch buffers are allocated and freed locally. +template static T -residual_error_comp(TestData &all_data, int64_t target_rank) { - auto m = all_data.row; - auto n = all_data.col; +residual_error_comp(T* A, int64_t m, int64_t n, + T* U, T* V, T* Sigma, int64_t target_rank, + T* U_scratch, T* V_scratch) { - all_data.U_cpy = new T[m * target_rank](); - all_data.V_cpy = new T[n * target_rank](); - - lapack::lacpy(MatrixType::General, m, target_rank, all_data.U, m, all_data.U_cpy, m); - lapack::lacpy(MatrixType::General, n, target_rank, all_data.V, n, all_data.V_cpy, n); - - // AV - US - // Scale columns of U by S + // S^{-1}AV - U + blas::gemm(Layout::ColMajor, Op::NoTrans, Op::NoTrans, m, target_rank, n, + (T)1, A, m, V, n, (T)0, U_scratch, m); for (int i = 0; i < target_rank; ++i) - blas::scal(m, all_data.Sigma[i], &all_data.U_cpy[m * i], 1); + blas::scal(m, (T)1 / Sigma[i], &U_scratch[m * i], 1); + blas::axpy(m * target_rank, (T)-1, U, 1, U_scratch, 1); - // Compute AV(:, 1:target_rank) - SU(1:target_rank) - blas::gemm(Layout::ColMajor, Op::NoTrans, Op::NoTrans, m, target_rank, n, 1.0, all_data.A, m, all_data.V, n, -1.0, all_data.U_cpy, m); - - // A'U - VS - // Scale columns of V by S + // (A'U)S^{-1} - V + blas::gemm(Layout::ColMajor, Op::Trans, Op::NoTrans, n, target_rank, m, + (T)1, A, m, U, m, (T)0, V_scratch, n); for (int i = 0; i < target_rank; ++i) - blas::scal(n, all_data.Sigma[i], &all_data.V_cpy[i * n], 1); - // Compute A'U(:, 1:target_rank) - VS(1:target_rank). - blas::gemm(Layout::ColMajor, Op::Trans, Op::NoTrans, n, target_rank, m, 1.0, all_data.A, m, all_data.U, m, -1.0, all_data.V_cpy, n); + blas::scal(n, (T)1 / Sigma[i], &V_scratch[i * n], 1); + blas::axpy(n * target_rank, (T)-1, V, 1, V_scratch, 1); - T nrm1 = lapack::lange(Norm::Fro, m, target_rank, all_data.U_cpy, m); - T nrm2 = lapack::lange(Norm::Fro, n, target_rank, all_data.V_cpy, n); + T nrm1 = lapack::lange(Norm::Fro, m, target_rank, U_scratch, m); + T nrm2 = lapack::lange(Norm::Fro, n, target_rank, V_scratch, n); return std::hypot(nrm1, nrm2); } -template -static T -approx_error_comp(ABRIK_benchmark_data &all_data, int64_t target_rank, T norm_A_lowrank) { - - auto m = all_data.row; - auto n = all_data.col; - - all_data.U_cpy = new T[m * target_rank](); - lapack::lacpy(MatrixType::General, m, target_rank, all_data.U, m, all_data.U_cpy, m); - - // U * S; scale the columns of U by S - for (int i = 0; i < target_rank; ++i) - blas::scal(m, all_data.Sigma[i], &all_data.U_cpy[i * m], 1); - - // U * S * V' - A_cpy ~= 0? - blas::gemm(Layout::ColMajor, Op::NoTrans, Op::Trans, m, n, target_rank, 1.0, all_data.U_cpy, m, all_data.V, n, -1.0, all_data.A_lowrank_svd, m); - - T nrm = lapack::lange(Norm::Fro, m, n, all_data.A_lowrank_svd, m); - printf("||A_hat_cursom_rank - A_svd_target_rank||_F / ||A_svd_target_rank||_F: %e\n", nrm / norm_A_lowrank); - - return nrm / norm_A_lowrank; -} - template static void call_all_algs( RandLAPACK::gen::mat_gen_info m_info, @@ -198,151 +142,146 @@ static void call_all_algs( int64_t b_sz, int64_t num_matmuls, int64_t target_rank, + bool run_gesdd, ABRIK_algorithm_objects &all_algs, ABRIK_benchmark_data &all_data, RandBLAS::RNGState &state, - std::string output_filename, - T norm_A_lowrank) { + std::ofstream &outfile) { + + using EMatrix = typename EigenTypes::Matrix; + using EVector = typename EigenTypes::Vector; - int i; auto m = all_data.row; auto n = all_data.col; auto tol = all_data.tolerance; // Additional params setup. all_algs.RSVD.block_sz = b_sz; - // Matrices R or S that give us the singular value spectrum returned by ABRIK will be of size b_sz * num_krylov_iters / 2. - // These matrices will be full-rank. - // Hence, target_rank = b_sz * num_krylov_iters / 2 - // ABRIK.max_krylov_iters = (int) ((target_rank * 2) / b_sz); - // - // Instead of the above approach, we now pre-specify the maximum number of Krylov iters that we allow for in num_matmuls. + // Instead of computing max_krylov_iters from target_rank, we pre-specify + // the maximum number of Krylov iters via num_matmuls. all_algs.ABRIK.max_krylov_iters = (int) num_matmuls; - all_algs.ABRIK.num_threads_min = 4; - all_algs.ABRIK.num_threads_max = RandLAPACK::util::get_omp_threads(); - + // timing vars long dur_ABRIK = 0; long dur_rsvd = 0; long dur_svds = 0; long dur_svd = 0; - // Making sure the states are unchanged - auto state_gen = state; + // Pre-allocate scratch buffers for residual_error_comp (avoid alloc/free per call). + T* U_scratch = new T[m * target_rank]; + T* V_scratch = new T[n * target_rank]; + auto state_alg = state; - T residual_err_custom_SVD = 0; + T residual_err_custom_SVD = 0; T residual_err_custom_ABRIK = 0; - T residual_err_custom_RSVD = 0; - T residual_err_custom_SVDS = 0; - - T lowrank_err_SVD = 0; - T lowrank_err_ABRIK = 0; - T lowrank_err_RSVD = 0; - T lowrank_err_SVDS = 0; + T residual_err_custom_RSVD = 0; + T residual_err_custom_SVDS = 0; int64_t singular_triplets_target_ABRIK = 0; int64_t singular_triplets_found_RSVD = 0; int64_t singular_triplets_target_RSVD = 0; - int64_t singular_triplets_found_SVDS = 0; int64_t singular_triplets_target_SVDS = 0; - for (i = 0; i < num_runs; ++i) { + for (int i = 0; i < num_runs; ++i) { printf("\nBlock size %ld, num matmuls %ld. Iteration %d start.\n", b_sz, num_matmuls, i); - - // Running ABRIK + + // ---- ABRIK ---- + // ABRIK allocates U, V, Sigma with new[] internally. auto start_ABRIK = steady_clock::now(); all_algs.ABRIK.call(m, n, all_data.A, m, b_sz, all_data.U, all_data.V, all_data.Sigma, state_alg); auto stop_ABRIK = steady_clock::now(); dur_ABRIK = duration_cast(stop_ABRIK - start_ABRIK).count(); printf("TOTAL TIME FOR ABRIK %ld\n", dur_ABRIK); - // This is in case the number of singular triplets is smaller than the target rank singular_triplets_target_ABRIK = std::min(target_rank, all_algs.ABRIK.singular_triplets_found); + residual_err_custom_ABRIK = residual_error_comp(all_data.A, m, n, all_data.U, all_data.V, all_data.Sigma, singular_triplets_target_ABRIK, U_scratch, V_scratch); + printf("ABRIK residual error: %.16e\n", residual_err_custom_ABRIK); - residual_err_custom_ABRIK = residual_error_comp(all_data, singular_triplets_target_ABRIK); - printf("ABRIK sqrt(||AV - SU||^2_F + ||A'U - VS||^2_F) / sqrt(target_rank): %.16e\n", residual_err_custom_ABRIK); + // Cleanup ABRIK outputs (new[]) + delete[] all_data.U; all_data.U = nullptr; + delete[] all_data.V; all_data.V = nullptr; + delete[] all_data.Sigma; all_data.Sigma = nullptr; - if (all_data.A_lowrank_svd != nullptr) - lowrank_err_ABRIK = approx_error_comp(all_data, singular_triplets_target_ABRIK, norm_A_lowrank); - state_alg = state; - state_gen = state; - data_regen(m_info, all_data, state_gen, 1); - - // Running RSVD - auto start_rsvd = steady_clock::now(); - // Below should technically be the same as - // all_algs.ABRIK.singular_triplets_found, unless ABRIK terminated early. - singular_triplets_found_RSVD = (int64_t ) (b_sz * num_matmuls / 2); + // Note: ABRIK does not modify A, so no regen_input needed. + + // ---- RSVD ---- + // RSVD (via QB) now modifies A in-place (deflation). We must copy A + // beforehand and restore it after, since ABRIK/SVDS need the original. + // RSVD allocates U, V, Sigma with calloc internally (via QB realloc chain). + singular_triplets_found_RSVD = (int64_t) (b_sz * num_matmuls / 2); - all_data.U = new T[m * singular_triplets_found_RSVD](); - all_data.V = new T[n * singular_triplets_found_RSVD](); - all_data.Sigma = new T[singular_triplets_found_RSVD](); + T* A_rsvd_copy = new T[m * n]; + lapack::lacpy(MatrixType::General, m, n, all_data.A, m, A_rsvd_copy, m); - all_algs.RSVD.call(m, n, all_data.A, singular_triplets_found_RSVD, tol, all_data.U, all_data.Sigma, all_data.V, state_alg); + auto start_rsvd = steady_clock::now(); + all_algs.RSVD.call(m, n, A_rsvd_copy, singular_triplets_found_RSVD, tol, all_data.U, all_data.Sigma, all_data.V, state_alg); auto stop_rsvd = steady_clock::now(); dur_rsvd = duration_cast(stop_rsvd - start_rsvd).count(); printf("TOTAL TIME FOR RSVD %ld\n", dur_rsvd); - // This is in case the number of singular triplets is smaller than the target rank singular_triplets_target_RSVD = std::min(target_rank, singular_triplets_found_RSVD); + residual_err_custom_RSVD = residual_error_comp(all_data.A, m, n, all_data.U, all_data.V, all_data.Sigma, singular_triplets_target_RSVD, U_scratch, V_scratch); + printf("RSVD residual error: %.16e\n", residual_err_custom_RSVD); - residual_err_custom_RSVD = residual_error_comp(all_data, singular_triplets_target_RSVD); - printf("RSVD sqrt(||AV - SU||^2_F + ||A'U - VS||^2_F) / sqrt(target_rank): %.16e\n", residual_err_custom_RSVD); + delete[] A_rsvd_copy; + + // Cleanup RSVD outputs (calloc) + free(all_data.U); all_data.U = nullptr; + free(all_data.V); all_data.V = nullptr; + free(all_data.Sigma); all_data.Sigma = nullptr; - if (all_data.A_lowrank_svd != nullptr) - lowrank_err_RSVD = approx_error_comp(all_data, singular_triplets_target_RSVD, norm_A_lowrank); - - state_alg = state; - state_gen = state; - data_regen(m_info, all_data, state_gen, 1); - - // Running SVDS - auto start_svds = steady_clock::now(); - - // Despite my earlier expectations, estimating a larger number of - // singular triplets via SVDS does improve the quality of the first singular triplets. - // As such, aiming for just the "target rank" would be unfair. - - // Below line also accounts for the case when number of singular triplets is smaller than the target rank. - singular_triplets_found_SVDS = std::min((int64_t ) (b_sz * num_matmuls / 2), n-2); - - printf("nev: %ld, nvc: %ld\n", singular_triplets_found_SVDS, std::min(2 * singular_triplets_found_SVDS, n-1)); - Spectra::PartialSVDSolver svds(all_data.A_spectra, singular_triplets_found_SVDS, std::min(2 * singular_triplets_found_SVDS, n-1)); - svds.compute(); - auto stop_svds = steady_clock::now(); - dur_svds = duration_cast(stop_svds - start_svds).count(); - printf("TOTAL TIME FOR SVDS %ld\n", dur_svds); - - // Copy data from Spectra (Eigen) format to the nomal C++. - Matrix U_spectra = svds.matrix_U(singular_triplets_found_SVDS); - Matrix V_spectra = svds.matrix_V(singular_triplets_found_SVDS); - Vector S_spectra = svds.singular_values(); - - all_data.U = new T[m * singular_triplets_found_SVDS](); - all_data.V = new T[n * singular_triplets_found_SVDS](); - all_data.Sigma = new T[singular_triplets_found_SVDS](); - - Eigen::Map(all_data.U, m, singular_triplets_found_SVDS) = U_spectra; - Eigen::Map(all_data.V, n, singular_triplets_found_SVDS) = V_spectra; - Eigen::Map(all_data.Sigma, singular_triplets_found_SVDS) = S_spectra; - - singular_triplets_target_SVDS = std::min(target_rank, singular_triplets_found_SVDS); - - residual_err_custom_SVDS = residual_error_comp(all_data, singular_triplets_target_SVDS); - printf("SVDS sqrt(||AV - SU||^2_F + ||A'U - VS||^2_F) / sqrt(target_rank): %.16e\n", residual_err_custom_SVDS); - - if (all_data.A_lowrank_svd != nullptr) - lowrank_err_SVDS = approx_error_comp(all_data, singular_triplets_target_SVDS, norm_A_lowrank); - state_alg = state; - state_gen = state; - data_regen(m_info, all_data, state_gen, 1); - + + // ---- SVDS (Spectra, budgeted) ---- + // Uses BudgetedPartialSVDSolver with fixed nev=target_rank and a + // matvec budget derived from (b_sz, num_matmuls) for fair comparison. + { + int64_t nev_svds = target_rank; + int64_t ncv_default = std::min(2 * nev_svds + 1, n - 1); + int64_t svds_budget = b_sz * num_matmuls; + int64_t ncv_svds = BenchmarkUtil::effective_ncv(svds_budget, nev_svds, ncv_default); + int64_t max_restarts = BenchmarkUtil::budget_to_restarts(svds_budget, nev_svds, ncv_svds); + + auto start_svds = steady_clock::now(); + Eigen::Map A_map(all_data.A, m, n); + BenchmarkUtil::BudgetedPartialSVDSolver svds(A_map, nev_svds, ncv_svds); + svds.compute(max_restarts); + auto stop_svds = steady_clock::now(); + dur_svds = duration_cast(stop_svds - start_svds).count(); + printf("SVDS: budget=%ld matvecs, max_restarts=%ld, A'A_ops=%ld, time=%ld us\n", + svds_budget, max_restarts, svds.num_operations(), dur_svds); + + // Extract ALL Ritz approximations (even unconverged) + EMatrix U_spectra = svds.matrix_U(nev_svds); + EMatrix V_spectra = svds.matrix_V(nev_svds); + EVector S_spectra = svds.singular_values(); + + all_data.U = new T[m * nev_svds](); + all_data.V = new T[n * nev_svds](); + all_data.Sigma = new T[nev_svds](); + + Eigen::Map(all_data.U, m, nev_svds) = U_spectra; + Eigen::Map(all_data.V, n, nev_svds) = V_spectra; + Eigen::Map(all_data.Sigma, nev_svds) = S_spectra; + + singular_triplets_target_SVDS = std::min(target_rank, nev_svds); + residual_err_custom_SVDS = residual_error_comp(all_data.A, m, n, all_data.U, all_data.V, all_data.Sigma, singular_triplets_target_SVDS, U_scratch, V_scratch); + printf("SVDS residual error: %.16e\n", residual_err_custom_SVDS); + + // Cleanup SVDS outputs (new[]) + delete[] all_data.U; all_data.U = nullptr; + delete[] all_data.V; all_data.V = nullptr; + delete[] all_data.Sigma; all_data.Sigma = nullptr; + + state_alg = state; + // Note: SVDS operates on A_spectra (const ref), does not modify A. + } + + // ---- SVD (GESDD) ---- // There is no reason to run SVD many times, as it always outputs the same result. - if ((b_sz == 16) && (num_matmuls == 4) && ((i == 0) || (i == 1))) { - // Running SVD + if (run_gesdd && (i == 0)) { auto start_svd = steady_clock::now(); all_data.U = new T[m * n](); all_data.Sigma = new T[n](); @@ -353,253 +292,161 @@ static void call_all_algs( dur_svd = duration_cast(stop_svd - start_svd).count(); printf("TOTAL TIME FOR SVD %ld\n", dur_svd); - // Standard SVD destorys matrix A, need to re-read it before running accuracy tests. - state_gen = state; - RandLAPACK::gen::mat_gen(m_info, all_data.A, state_gen); + // GESDD destroys A, re-read before residual computation. + regen_input(m_info, all_data, state); RandLAPACK::util::transposition(n, n, all_data.VT, n, all_data.V, n, 0); - residual_err_custom_SVD = residual_error_comp(all_data, target_rank); - printf("SVD sqrt(||AV - US||^2_F + ||A'U - VS||^2_F) / sqrt(target_rank): %.16e\n", residual_err_custom_SVD); + residual_err_custom_SVD = residual_error_comp(all_data.A, m, n, all_data.U, all_data.V, all_data.Sigma, target_rank, U_scratch, V_scratch); + printf("SVD residual error: %.16e\n", residual_err_custom_SVD); - if (all_data.A_lowrank_svd != nullptr) - lowrank_err_SVD = approx_error_comp(all_data, target_rank, norm_A_lowrank); + // Cleanup SVD outputs (new[]) + delete[] all_data.U; all_data.U = nullptr; + delete[] all_data.VT; all_data.VT = nullptr; + delete[] all_data.V; all_data.V = nullptr; + delete[] all_data.Sigma; all_data.Sigma = nullptr; state_alg = state; - state_gen = state; - data_regen(m_info, all_data, state_gen, 1); + regen_input(m_info, all_data, state); } - std::ofstream file(output_filename, std::ios::app); - file << b_sz << ", " << all_algs.ABRIK.max_krylov_iters << ", " << target_rank << ", " - << residual_err_custom_ABRIK << ", " << lowrank_err_ABRIK << ", " << dur_ABRIK << ", " - << residual_err_custom_RSVD << ", " << lowrank_err_RSVD << ", " << dur_rsvd << ", " - << residual_err_custom_SVDS << ", " << lowrank_err_SVDS << ", " << dur_svds << ", " - << residual_err_custom_SVD << ", " << lowrank_err_SVD << ", " << dur_svd << ",\n"; + // Write CSV data row + outfile << b_sz << ", " << all_algs.ABRIK.max_krylov_iters << ", " << target_rank << ", " + << residual_err_custom_ABRIK << ", " << dur_ABRIK << ", " + << residual_err_custom_RSVD << ", " << dur_rsvd << ", " + << residual_err_custom_SVDS << ", " << dur_svds << ", " + << residual_err_custom_SVD << ", " << dur_svd << "\n"; + outfile.flush(); } + + delete[] U_scratch; + delete[] V_scratch; } -/* -int main(int argc, char *argv[]) { +template +static void run_benchmark(int argc, char *argv[]) { + using EMatrix = typename EigenTypes::Matrix; if (argc < 12) { - // Expected input into this benchmark. - std::cerr << "Usage: " << argv[0] << " " << std::endl; - return 1; + std::cerr << "Usage: " << argv[0] << " " << std::endl; + return; } int num_runs = std::stol(argv[4]); int64_t m_expected = std::stol(argv[5]); int64_t n_expected = std::stol(argv[6]); int64_t target_rank = std::stol(argv[7]); + bool run_gesdd = (std::stoi(argv[8]) != 0); std::vector b_sz; - for (int i = 0; i < std::stol(argv[8]); ++i) - b_sz.push_back(std::stoi(argv[i + 10])); - // Save elements in string for logging purposes + for (int i = 0; i < std::stol(argv[9]); ++i) + b_sz.push_back(std::stoi(argv[i + 11])); std::ostringstream oss1; for (const auto &val : b_sz) oss1 << val << ", "; std::string b_sz_string = oss1.str(); std::vector matmuls; - for (int i = 0; i < std::stol(argv[9]); ++i) - matmuls.push_back(std::stoi(argv[i + 10 + std::stol(argv[8])])); - // Save elements in string for logging purposes + for (int i = 0; i < std::stol(argv[10]); ++i) + matmuls.push_back(std::stoi(argv[i + 11 + std::stol(argv[9])])); std::ostringstream oss2; for (const auto &val : matmuls) oss2 << val << ", "; std::string matmuls_string = oss2.str(); - double tol = std::pow(std::numeric_limits::epsilon(), 0.85); + T tol = std::pow(std::numeric_limits::epsilon(), (T)0.85); auto state = RandBLAS::RNGState(); auto state_constant = state; - double norm_A_lowrank = 0; int64_t m = 0, n = 0; // Generate the input matrix. - RandLAPACK::gen::mat_gen_info m_info(m, n, RandLAPACK::gen::custom_input); - m_info.filename = argv[2]; + RandLAPACK::gen::mat_gen_info m_info(m, n, RandLAPACK::gen::custom_input); + m_info.filename = argv[3]; m_info.workspace_query_mod = 1; - // Workspace query; - RandLAPACK::gen::mat_gen(m_info, NULL, state); + RandLAPACK::gen::mat_gen(m_info, NULL, state); // Update basic params. m = m_info.rows; n = m_info.cols; if (m_expected != m || n_expected != n) { - std::cerr << "Expected input size (" << m_expected << ", " << n_expected << ") did not matrch actual input size (" << m << ", " << n << "). Aborting." << std::endl; - return 1; + std::cerr << "Expected input size (" << m_expected << ", " << n_expected << ") did not match actual input size (" << m << ", " << n << "). Aborting." << std::endl; + return; } // Allocate basic workspace. - ABRIK_benchmark_data all_data(m, n, tol); - // Fill the data matrix; + ABRIK_benchmark_data all_data(m, n, tol); RandLAPACK::gen::mat_gen(m_info, all_data.A, state); // Declare objects for RSVD and ABRIK - int64_t p = 5; + int64_t p = 2; int64_t passes_per_iteration = 1; - // Block size will need to be altered. int64_t block_sz = 0; - ABRIK_algorithm_objects all_algs(false, false, false, false, p, passes_per_iteration, block_sz, tol); + ABRIK_algorithm_objects all_algs(false, false, false, false, p, passes_per_iteration, block_sz, tol); // Copying input data into a Spectra (Eigen) matrix object - Eigen::Map(all_data.A_spectra.data(), all_data.A_spectra.rows(), all_data.A_spectra.cols()) = Eigen::Map(all_data.A, m, n); - - // Optional pass of lowrank SVD matrix into the benchmark - if (std::string(argv[3]) != ".") { - printf("Lowrank A input.\n"); - RandLAPACK::gen::mat_gen_info m_info_A_svd(m, n, RandLAPACK::gen::custom_input); - m_info_A_svd.filename = argv[3]; - m_info_A_svd.workspace_query_mod = 0; - all_data.A_lowrank_svd = new double[m * n](); - all_data.A_lowrank_svd_const = new double[m * n](); - RandLAPACK::gen::mat_gen(m_info_A_svd, all_data.A_lowrank_svd_const, state); - lapack::lacpy(MatrixType::General, m, n, all_data.A_lowrank_svd_const, m, all_data.A_lowrank_svd, m); - - // Pre-compute norm(A lowrank) for future benchmarking - norm_A_lowrank = lapack::lange(Norm::Fro, m, n, all_data.A_lowrank_svd, m); - } - printf("Finished data preparation\n"); - // Declare a data file - std::string output_filename = "_ABRIK_speed_comparisons_num_info_lines_" + std::to_string(6) + ".txt"; + + // Generate date/time prefix + std::time_t now = std::time(nullptr); + char date_prefix[20]; + std::strftime(date_prefix, sizeof(date_prefix), "%Y%m%d_%H%M%S_", std::localtime(&now)); + + // Build output file path + std::string output_filename = std::string(date_prefix) + "ABRIK_speed_comparisons.csv"; std::string path; - if (std::string(argv[1]) != ".") { - path = argv[1] + output_filename; + if (std::string(argv[2]) != ".") { + path = std::string(argv[2]) + "/" + output_filename; } else { path = output_filename; } std::ofstream file(path, std::ios::out | std::ios::app); - // Writing important data into file - file << "Description: Results from the ABRIK speed comparison benchmark, recording the time it takes to perform ABRIK and alternative methods for low-rank SVD." - "\nFile format: 15 columns, showing krylov block size, nummber of matmuls permitted, and num svals and svecs to approximate, followed by the residual error, standard lowrank error and execution time for all algorithms (ABRIK, RSVD, SVDS, SVD)" - "\n Rows correspond to algorithm runs with Krylov block sizes varying as specified, and numbers of matmuls varying as specified per eah block size, with num_runs repititions of each number of matmuls." - "\nInput type:" + std::string(argv[2]) + - "\nInput size:" + std::to_string(m) + " by " + std::to_string(n) + - "\nAdditional parameters: Krylov block sizes " + b_sz_string + - " matmuls: " + matmuls_string + - " num runs per size " + std::to_string(num_runs) + - " num singular values and vectors approximated " + std::to_string(target_rank) + - "\n"; + // Write metadata header (prefixed with # for easy parsing) + file << "# ABRIK Speed Comparison Benchmark\n" + << "# Precision: " << argv[1] << "\n" + << "# Input matrix: " << argv[3] << "\n" + << "# Input size: " << m << " x " << n << "\n" + << "# Target rank: " << target_rank << "\n" + << "# Krylov block sizes: " << b_sz_string << "\n" + << "# Matmul counts: " << matmuls_string << "\n" + << "# Runs per configuration: " << num_runs << "\n" + << "# Tolerance: " << tol << "\n" + << "# Run GESDD: " << (run_gesdd ? "yes" : "no") << "\n" + << "# Residual metric: sqrt(||S^{-1}AV - U||^2_F + ||(A'U)S^{-1} - V||^2_F)\n" + << "# Timings in microseconds\n"; + // Write CSV column header + file << "b_sz, num_matmuls, target_rank, " + << "err_ABRIK, dur_ABRIK, " + << "err_RSVD, dur_RSVD, " + << "err_SVDS, dur_SVDS, " + << "err_SVD, dur_SVD\n"; file.flush(); + auto start_total = steady_clock::now(); + size_t i = 0, j = 0; - for (;i < b_sz.size(); ++i) { - for (;j < matmuls.size(); ++j) { - call_all_algs(m_info, num_runs, b_sz[i], matmuls[j], target_rank, all_algs, all_data, state_constant, path, norm_A_lowrank); + for (; i < b_sz.size(); ++i) { + for (; j < matmuls.size(); ++j) { + call_all_algs(m_info, num_runs, b_sz[i], matmuls[j], target_rank, run_gesdd, all_algs, all_data, state_constant, file); } j = 0; } + + auto stop_total = steady_clock::now(); + long total_us = duration_cast(stop_total - start_total).count(); + printf("\nTOTAL BENCHMARK TIME: %.2f seconds\n", total_us / 1e6); } -*/ int main(int argc, char *argv[]) { - - if (argc < 12) { - // Expected input into this benchmark. - std::cerr << "Usage: " << argv[0] << " " << std::endl; + if (argc < 2) { + std::cerr << "Usage: " << argv[0] << " " << std::endl; return 1; } - int num_runs = std::stol(argv[4]); - int64_t m_expected = std::stol(argv[5]); - int64_t n_expected = std::stol(argv[6]); - int64_t target_rank = std::stol(argv[7]); - std::vector b_sz; - for (int i = 0; i < std::stol(argv[8]); ++i) - b_sz.push_back(std::stoi(argv[i + 10])); - // Save elements in string for logging purposes - std::ostringstream oss1; - for (const auto &val : b_sz) - oss1 << val << ", "; - std::string b_sz_string = oss1.str(); - std::vector matmuls; - for (int i = 0; i < std::stol(argv[9]); ++i) - matmuls.push_back(std::stoi(argv[i + 10 + std::stol(argv[8])])); - // Save elements in string for logging purposes - std::ostringstream oss2; - for (const auto &val : matmuls) - oss2 << val << ", "; - std::string matmuls_string = oss2.str(); - float tol = std::pow(std::numeric_limits::epsilon(), 0.85); - auto state = RandBLAS::RNGState(); - auto state_constant = state; - float norm_A_lowrank = 0; - int64_t m = 0, n = 0; - - // Generate the input matrix. - RandLAPACK::gen::mat_gen_info m_info(m, n, RandLAPACK::gen::custom_input); - m_info.filename = argv[2]; - m_info.workspace_query_mod = 1; - // Workspace query; - RandLAPACK::gen::mat_gen(m_info, NULL, state); - - // Update basic params. - m = m_info.rows; - n = m_info.cols; - if (m_expected != m || n_expected != n) { - std::cerr << "Expected input size (" << m_expected << ", " << n_expected << ") did not matrch actual input size (" << m << ", " << n << "). Aborting." << std::endl; - return 1; - } - - // Allocate basic workspace. - ABRIK_benchmark_data all_data(m, n, tol); - // Fill the data matrix; - RandLAPACK::gen::mat_gen(m_info, all_data.A, state); - - // Declare objects for RSVD and ABRIK - int64_t p = 2; - int64_t passes_per_iteration = 1; - // Block size will need to be altered. - int64_t block_sz = 0; - ABRIK_algorithm_objects all_algs(false, false, false, false, p, passes_per_iteration, block_sz, tol); - - // Copying input data into a Spectra (Eigen) matrix object - Eigen::Map(all_data.A_spectra.data(), all_data.A_spectra.rows(), all_data.A_spectra.cols()) = Eigen::Map(all_data.A, m, n); - - // Optional pass of lowrank SVD matrix into the benchmark - if (std::string(argv[3]) != ".") { - printf("Lowrank A input.\n"); - RandLAPACK::gen::mat_gen_info m_info_A_svd(m, n, RandLAPACK::gen::custom_input); - m_info_A_svd.filename = argv[3]; - m_info_A_svd.workspace_query_mod = 0; - all_data.A_lowrank_svd = new float[m * n](); - all_data.A_lowrank_svd_const = new float[m * n](); - RandLAPACK::gen::mat_gen(m_info_A_svd, all_data.A_lowrank_svd_const, state); - lapack::lacpy(MatrixType::General, m, n, all_data.A_lowrank_svd_const, m, all_data.A_lowrank_svd, m); - - // Pre-compute norm(A lowrank) for future benchmarking - norm_A_lowrank = lapack::lange(Norm::Fro, m, n, all_data.A_lowrank_svd, m); - } - - printf("Finished data preparation\n"); - // Declare a data file - std::string output_filename = "_ABRIK_speed_comparisons_num_info_lines_" + std::to_string(6) + ".txt"; - std::string path; - if (std::string(argv[1]) != ".") { - path = argv[1] + output_filename; + std::string precision = argv[1]; + if (precision == "double") { + run_benchmark(argc, argv); + } else if (precision == "float") { + run_benchmark(argc, argv); } else { - path = output_filename; - } - std::ofstream file(path, std::ios::out | std::ios::app); - - // Writing important data into file - file << "Description: Results from the ABRIK speed comparison benchmark, recording the time it takes to perform ABRIK and alternative methods for low-rank SVD." - "\nFile format: 15 columns, showing krylov block size, nummber of matmuls permitted, and num svals and svecs to approximate, followed by the residual error, standard lowrank error and execution time for all algorithms (ABRIK, RSVD, SVDS, SVD)" - "\n Rows correspond to algorithm runs with Krylov block sizes varying as specified, and numbers of matmuls varying as specified per eah block size, with num_runs repititions of each number of matmuls." - "\nInput type:" + std::string(argv[2]) + - "\nInput size:" + std::to_string(m) + " by " + std::to_string(n) + - "\nAdditional parameters: Krylov block sizes " + b_sz_string + - " matmuls: " + matmuls_string + - " num runs per size " + std::to_string(num_runs) + - " num singular values and vectors approximated " + std::to_string(target_rank) + - "\n"; - file.flush(); - - size_t i = 0, j = 0; - for (;i < b_sz.size(); ++i) { - for (;j < matmuls.size(); ++j) { - call_all_algs(m_info, num_runs, b_sz[i], matmuls[j], target_rank, all_algs, all_data, state_constant, path, norm_A_lowrank); - } - j = 0; + std::cerr << "Error: precision must be 'double' or 'float', got '" << precision << "'" << std::endl; + return 1; } -} \ No newline at end of file + return 0; +} diff --git a/benchmark/bench_ABRIK/ABRIK_speed_comparisons_sparse.cc b/benchmark/bench_ABRIK/ABRIK_speed_comparisons_sparse.cc index e49063d67..5eb7b4799 100644 --- a/benchmark/bench_ABRIK/ABRIK_speed_comparisons_sparse.cc +++ b/benchmark/bench_ABRIK/ABRIK_speed_comparisons_sparse.cc @@ -1,9 +1,13 @@ /* -Additional ABRIK speed comparison benchmark - runs ABRIK, RSVD and SVDS from Spectra library. -The user is required to provide a matrix file to be read, set min and max numbers of large gemms (Krylov iterations) that the algorithm is allowed to perform min and max block sizes that ABRIK is to use; -furthermore, the user is to provide a 'custom rank' parameter (number of singular vectors to approximate by ABRIK). -The benchmark outputs the basic data of a given run, as well as the ABRIK runtime and singular vector residual error, -which is computed as "sqrt(||AV - SU||^2_F + ||A'U - VS||^2_F / sqrt(target_rank)" (for "custom rank" singular vectors and values). +ABRIK speed comparison benchmark (sparse) - runs ABRIK and optionally SVDS (Spectra) / GESDD +on sparse input matrices read from Matrix Market format. +Precision (float or double) is specified as the first CLI argument. + +Output: CSV file with '#'-prefixed metadata header, column names, then data rows. +ABRIK allocates U/V/Sigma with new[] internally -> cleanup with delete[]. + +Residual metric: sqrt(||S^{-1}AV - U||^2_F + ||(A'U)S^{-1} - V||^2_F). +Timings in microseconds. */ #include "RandLAPACK.hh" @@ -15,6 +19,11 @@ which is computed as "sqrt(||AV - SU||^2_F + ||A'U - VS||^2_F / sqrt(target_rank #include #include #include +#include +#include +#include +#include +#include // External libs includes #include @@ -22,13 +31,43 @@ which is computed as "sqrt(||AV - SU||^2_F + ||A'U - VS||^2_F / sqrt(target_rank #include #include #include - -using SpMatrix = Eigen::SparseMatrix; -using Matrix = Eigen::MatrixXd; -using Vector = Eigen::VectorXd; +#include "BudgetedSVDSolver.hh" + +// Traits struct mapping scalar type T to Eigen matrix/vector/sparse types. +template struct EigenTypes; +template <> struct EigenTypes { + using Matrix = Eigen::MatrixXd; + using Vector = Eigen::VectorXd; + using SpMatrix = Eigen::SparseMatrix; +}; +template <> struct EigenTypes { + using Matrix = Eigen::MatrixXf; + using Vector = Eigen::VectorXf; + using SpMatrix = Eigen::SparseMatrix; +}; using Subroutines = RandLAPACK::ABRIKSubroutines; +// Helper function to ensure directory exists (creates parent directories if needed) +void ensure_directory_exists(const std::string& path) { + struct stat info; + if (stat(path.c_str(), &info) == 0) { + if (info.st_mode & S_IFDIR) return; + std::cerr << "Error: " << path << " exists but is not a directory" << std::endl; + return; + } + size_t pos = path.find_last_of('/'); + if (pos != std::string::npos && pos > 0) { + ensure_directory_exists(path.substr(0, pos)); + } + if (mkdir(path.c_str(), 0755) != 0) { + std::cerr << "Warning: Could not create directory " << path + << " (error: " << strerror(errno) << ")" << std::endl; + } else { + std::cout << "Created output directory: " << path << std::endl; + } +} + template struct ABRIK_benchmark_data { int64_t row; @@ -37,34 +76,21 @@ struct ABRIK_benchmark_data { SpMat* A_input; RandLAPACK::linops::SpLinOp A_linop; T* U; - T* V; + T* V; T* Sigma; - T* U_cpy; - T* V_cpy; - SpMatrix A_spectra; ABRIK_benchmark_data(SpMat& input_mat, int64_t m, int64_t n, T tol) : - A_spectra(m, n), A_linop(m, n, input_mat) { U = nullptr; V = nullptr; Sigma = nullptr; - U_cpy = nullptr; - V_cpy = nullptr; - row = m; col = n; tolerance = tol; } - ~ABRIK_benchmark_data() { - delete[] U; - delete[] V; - delete[] Sigma; - delete[] U_cpy; - delete[] V_cpy; - } + ~ABRIK_benchmark_data() {} }; template @@ -72,8 +98,8 @@ struct ABRIK_algorithm_objects { RandLAPACK::ABRIK ABRIK; ABRIK_algorithm_objects( - bool verbosity, - bool time_subroutines, + bool verbosity, + bool time_subroutines, T tol ) : ABRIK(verbosity, time_subroutines, tol) @@ -82,7 +108,6 @@ struct ABRIK_algorithm_objects { template RandBLAS::sparse_data::coo::COOMatrix from_matrix_market(std::string fn) { - int64_t n_rows, n_cols = 0; std::vector rows{}; std::vector cols{}; @@ -90,11 +115,11 @@ RandBLAS::sparse_data::coo::COOMatrix from_matrix_market(std::string fn) { std::ifstream file_stream(fn); fast_matrix_market::read_matrix_market_triplet( - file_stream, n_rows, n_cols, rows, cols, vals + file_stream, n_rows, n_cols, rows, cols, vals ); RandBLAS::sparse_data::coo::COOMatrix out(n_rows, n_cols); - reserve_coo(vals.size(),out); + reserve_coo(vals.size(), out); for (int i = 0; i < out.nnz; ++i) { out.rows[i] = rows[i]; out.cols[i] = cols[i]; @@ -105,8 +130,7 @@ RandBLAS::sparse_data::coo::COOMatrix from_matrix_market(std::string fn) { } template -void from_matrix_market(std::string fn, SpMatrix& A) { - +void from_matrix_market(std::string fn, typename EigenTypes::SpMatrix& A) { int64_t n_rows, n_cols = 0; std::vector rows{}; std::vector cols{}; @@ -114,18 +138,18 @@ void from_matrix_market(std::string fn, SpMatrix& A) { std::ifstream file_stream(fn); fast_matrix_market::read_matrix_market_triplet( - file_stream, n_rows, n_cols, rows, cols, vals + file_stream, n_rows, n_cols, rows, cols, vals ); std::vector> tripletList; - for (int i = 0; i < vals.size(); ++i) + for (size_t i = 0; i < vals.size(); ++i) tripletList.push_back(Eigen::Triplet(rows[i], cols[i], vals[i])); A.setFromTriplets(tripletList.begin(), tripletList.end()); } template -void from_matrix_market_leading_submatrix(std::string fn, SpMatrix& A, T submatrix_dim_ratio) { +void from_matrix_market_leading_submatrix(std::string fn, typename EigenTypes::SpMatrix& A, T submatrix_dim_ratio) { int64_t n_rows, n_cols = 0; std::vector rows{}; std::vector cols{}; @@ -136,16 +160,13 @@ void from_matrix_market_leading_submatrix(std::string fn, SpMatrix& A, T submatr file_stream, n_rows, n_cols, rows, cols, vals ); - // Use half-size for the leading principal submatrix int64_t m = n_rows * submatrix_dim_ratio; int64_t n = n_cols * submatrix_dim_ratio; - // Create triplets only for entries within the submatrix std::vector> tripletList; for (size_t i = 0; i < vals.size(); ++i) { - if (rows[i] < m && cols[i] < n) { + if (rows[i] < m && cols[i] < n) tripletList.emplace_back(rows[i], cols[i], vals[i]); - } } A.setFromTriplets(tripletList.begin(), tripletList.end()); @@ -154,18 +175,14 @@ void from_matrix_market_leading_submatrix(std::string fn, SpMatrix& A, T submatr template RandBLAS::CSCMatrix format_conversion(int64_t m, int64_t n, RandBLAS::COOMatrix& input_mat_coo) { - // Grab the leading principal submatrix of the size of half the input - RandBLAS::COOMatrix input_mat_transformed(m, n); + RandBLAS::COOMatrix input_mat_transformed(m, n); - // check how many nonzeros are in the left principal submatrix int64_t nnz_sub = 0; for (int64_t i = 0; i < input_mat_coo.nnz; ++i) { - if (input_mat_coo.rows[i] < m && input_mat_coo.cols[i] < n) { + if (input_mat_coo.rows[i] < m && input_mat_coo.cols[i] < n) ++nnz_sub; - } } - // Allocate reserve_coo(nnz_sub, input_mat_transformed); int64_t ell = 0; @@ -178,65 +195,123 @@ RandBLAS::CSCMatrix format_conversion(int64_t m, int64_t n, RandBLAS::COOMatr } } - // Convert the sparse matrix format for performance - RandBLAS::CSCMatrix input_mat_csc(m, n); - //RandBLAS::conversions::coo_to_csc(input_mat_coo, input_mat_csc); + RandBLAS::CSCMatrix input_mat_csc(m, n); RandBLAS::conversions::coo_to_csc(input_mat_transformed, input_mat_csc); return input_mat_csc; } +// Computes the residual norm error: sqrt(||S^{-1}AV - U||^2_F + ||(A'U)S^{-1} - V||^2_F). +// Uses LinOp for sparse matvec. Pre-allocated scratch buffers U_scratch, V_scratch are used. +template +static T +residual_error_vectors_comp(LinOp& A_linop, int64_t m, int64_t n, + T* U, T* V, T* Sigma, int64_t target_rank, + T* U_scratch, T* V_scratch) { -// Re-generate and clear data -template -static void data_regen(ABRIK_benchmark_data &all_data, - RandBLAS::RNGState &state) { - - delete[] all_data.U; - delete[] all_data.V; - delete[] all_data.Sigma; - delete[] all_data.U_cpy; - delete[] all_data.V_cpy; - all_data.U = nullptr; - all_data.V = nullptr; - all_data.Sigma = nullptr; - all_data.U_cpy = nullptr; - all_data.V_cpy = nullptr; + // S^{-1}AV - U + A_linop(Side::Left, Layout::ColMajor, Op::NoTrans, Op::NoTrans, m, target_rank, n, + (T)1, V, n, (T)0, U_scratch, m); + for (int i = 0; i < target_rank; ++i) + blas::scal(m, (T)1 / Sigma[i], &U_scratch[m * i], 1); + blas::axpy(m * target_rank, (T)-1, U, 1, U_scratch, 1); + + // (A'U)S^{-1} - V + A_linop(Side::Left, Layout::ColMajor, Op::Trans, Op::NoTrans, n, target_rank, m, + (T)1, U, m, (T)0, V_scratch, n); + for (int i = 0; i < target_rank; ++i) + blas::scal(n, (T)1 / Sigma[i], &V_scratch[i * n], 1); + blas::axpy(n * target_rank, (T)-1, V, 1, V_scratch, 1); + + T nrm1 = lapack::lange(Norm::Fro, m, target_rank, U_scratch, m); + T nrm2 = lapack::lange(Norm::Fro, n, target_rank, V_scratch, n); + + return std::hypot(nrm1, nrm2); } -// This routine computes the residual norm error, consisting of two parts (one of which) vanishes -// in exact precision. Target_rank defines size of U, V as returned by ABRIK; target_rank <= target_rank. -template -static T -residual_error_comp(TestData &all_data, int64_t target_rank) { - auto m = all_data.row; - auto n = all_data.col; +// Helper function to write matrices to Matrix Market format (MATLAB-compatible) +template +void write_matrix_to_file(const std::string& filename, const T* matrix, int64_t rows, int64_t cols, bool is_vector = false) { + std::ofstream file(filename); + if (!file.is_open()) { + std::cerr << "Error: Could not open file " << filename << " for writing." << std::endl; + return; + } + file << "%%MatrixMarket matrix array real general\n"; + file << rows << " " << cols << "\n"; + file << std::scientific << std::setprecision(16); + for (int64_t j = 0; j < cols; ++j) { + for (int64_t i = 0; i < rows; ++i) { + file << matrix[i + j * rows] << "\n"; + } + } + file.close(); + std::cout << "Successfully wrote " << (is_vector ? "vector" : "matrix") + << " (" << rows << " x " << cols << ") to " << filename << std::endl; +} - all_data.U_cpy = new T[m * target_rank](); - all_data.V_cpy = new T[n * target_rank](); +// Perform direct SVD using LAPACK's GESDD on dense matrix +template +void compute_direct_gesdd(const typename EigenTypes::SpMatrix& A_sparse, int64_t m, int64_t n, int64_t target_rank, + std::string output_dir, bool write_output_matrices) { + printf("\n========== Running Direct GESDD ==========\n"); - lapack::lacpy(MatrixType::General, m, target_rank, all_data.U, m, all_data.U_cpy, m); - lapack::lacpy(MatrixType::General, n, target_rank, all_data.V, n, all_data.V_cpy, n); + // Convert sparse to dense + printf("Converting sparse matrix to dense...\n"); + T* A_dense = new T[m * n](); - // AV - US - // Scale columns of U by S - for (int i = 0; i < target_rank; ++i) - blas::scal(m, all_data.Sigma[i], &all_data.U_cpy[m * i], 1); + for (int k = 0; k < A_sparse.outerSize(); ++k) { + for (typename EigenTypes::SpMatrix::InnerIterator it(A_sparse, k); it; ++it) { + A_dense[it.row() + it.col() * m] = it.value(); + } + } - // Compute AV(:, 1:target_rank) - SU(1:target_rank) - all_data.A_linop(Side::Left, Layout::ColMajor, Op::NoTrans, Op::NoTrans, m, target_rank, n, 1.0, all_data.V, n, -1.0, all_data.U_cpy, m); + int64_t min_mn = std::min(m, n); - // A'U - VS - // Scale columns of V by S - for (int i = 0; i < target_rank; ++i) - blas::scal(n, all_data.Sigma[i], &all_data.V_cpy[i * n], 1); - // Compute A'U(:, 1:target_rank) - VS(1:target_rank). - all_data.A_linop(Side::Left, Layout::ColMajor, Op::Trans, Op::NoTrans, n, target_rank, m, 1.0, all_data.U, m, -1.0, all_data.V_cpy, n); + printf("Computing full SVD with GESDD...\n"); + auto start_gesdd = steady_clock::now(); - T nrm1 = lapack::lange(Norm::Fro, m, target_rank, all_data.U_cpy, m); - T nrm2 = lapack::lange(Norm::Fro, n, target_rank, all_data.V_cpy, n); + T* S_full = new T[min_mn](); + T* U_full = new T[m * min_mn](); + T* VT_full = new T[min_mn * n](); - return std::hypot(nrm1, nrm2); + T* A_copy = new T[m * n](); + lapack::lacpy(MatrixType::General, m, n, A_dense, m, A_copy, m); + + lapack::gesdd(lapack::Job::SomeVec, m, n, A_copy, m, S_full, U_full, m, VT_full, min_mn); + + auto stop_gesdd = steady_clock::now(); + long dur_gesdd = duration_cast(stop_gesdd - start_gesdd).count(); + printf("TOTAL TIME FOR GESDD: %ld microseconds\n", dur_gesdd); + + // Transpose VT to get V + T* V_full = new T[n * min_mn](); + for (int64_t i = 0; i < min_mn; ++i) { + for (int64_t j = 0; j < n; ++j) { + V_full[j + i * n] = VT_full[i + j * min_mn]; + } + } + + printf("GESDD completed. First 5 singular values:\n"); + for (int64_t i = 0; i < std::min((int64_t)5, min_mn); ++i) { + printf(" Sigma[%ld] = %.16e\n", i, S_full[i]); + } + + if (write_output_matrices) { + std::string prefix = output_dir + "/GESDD"; + write_matrix_to_file(prefix + "_U.mtx", U_full, m, min_mn, false); + write_matrix_to_file(prefix + "_V.mtx", V_full, n, min_mn, false); + write_matrix_to_file(prefix + "_Sigma.mtx", S_full, min_mn, 1, true); + } + + delete[] A_dense; + delete[] A_copy; + delete[] VT_full; + delete[] V_full; + delete[] S_full; + delete[] U_full; + + printf("==========================================\n\n"); } template @@ -247,191 +322,260 @@ static void call_all_algs( int64_t target_rank, ABRIK_algorithm_objects &all_algs, ABRIK_benchmark_data &all_data, + const typename EigenTypes::SpMatrix &A_spectra_ref, RandBLAS::RNGState &state, - std::string output_filename, - std::string input_path) { + std::ofstream &outfile, + std::string input_path, + std::string output_dir, + bool write_output_matrices) { + + using ESpMatrix = typename EigenTypes::SpMatrix; + using EMatrix = typename EigenTypes::Matrix; + using EVector = typename EigenTypes::Vector; - int i; auto m = all_data.row; auto n = all_data.col; - auto tol = all_data.tolerance; // Additional params setup. - // Matrices R or S that give us the singular value spectrum returned by ABRIK will be of size b_sz * num_krylov_iters / 2. - // These matrices will be full-rank. - // Hence, target_rank = b_sz * num_krylov_iters / 2 - // ABRIK.max_krylov_iters = (int) ((target_rank * 2) / b_sz); - // - // Instead of the above approach, we now pre-specify the maximum number of Krylov iters that we allow for in num_matmuls. all_algs.ABRIK.max_krylov_iters = (int) num_matmuls; - all_algs.ABRIK.num_threads_min = 4; - all_algs.ABRIK.num_threads_max = RandLAPACK::util::get_omp_threads(); - // Useful for all sparse matrices except 0. - all_algs.ABRIK.qr_exp = Subroutines::QR_explicit::cqrrt; - + // timing vars long dur_ABRIK = 0; long dur_svds = 0; - // Making sure the states are unchanged - auto state_gen = state; + // Pre-allocate scratch buffers for residual_error_vectors_comp (avoid alloc/free per call). + T* U_scratch = new T[m * target_rank]; + T* V_scratch = new T[n * target_rank]; + auto state_alg = state; - T residual_err_custom_SVDS = 0; - T residual_err_custom_ABRIK = 0; + T residual_err_vec_SVDS = 0; + T residual_err_vec_ABRIK = 0; int64_t singular_triplets_target_ABRIK = 0; int64_t singular_triplets_found_SVDS = 0; int64_t singular_triplets_target_SVDS = 0; - for (i = 0; i < num_runs; ++i) { + for (int i = 0; i < num_runs; ++i) { printf("\nBlock size %ld, num matmuls %ld. Iteration %d start.\n", b_sz, num_matmuls, i); - // Running ABRIK + // ---- ABRIK ---- + // ABRIK allocates U, V, Sigma with new[] internally. auto start_ABRIK = steady_clock::now(); all_algs.ABRIK.call(m, n, *all_data.A_input, m, b_sz, all_data.U, all_data.V, all_data.Sigma, state_alg); auto stop_ABRIK = steady_clock::now(); dur_ABRIK = duration_cast(stop_ABRIK - start_ABRIK).count(); printf("TOTAL TIME FOR ABRIK %ld\n", dur_ABRIK); - // This is in case the number of singular triplets is smaller than the target rank singular_triplets_target_ABRIK = std::min(target_rank, all_algs.ABRIK.singular_triplets_found); printf("Singular triplets: %ld\n", singular_triplets_target_ABRIK); - residual_err_custom_ABRIK = residual_error_comp(all_data, singular_triplets_target_ABRIK); - printf("ABRIK sqrt(||AV - SU||^2_F + ||A'U - VS||^2_F) / sqrt(target_rank): %.16e\n", residual_err_custom_ABRIK); - - state_alg = state; - state_gen = state; - data_regen(all_data, state_gen); - - - // Running SVDS - auto start_svds = steady_clock::now(); - - // Despite my earlier expectations, estimating a larger number of - // singular triplets via SVDS does improve the quality of the first singular triplets. - // As such, aiming for just the "target rank" would be unfair. - - // Below line also accounts for the case when number of singular triplets is smaller than the target rank. - singular_triplets_found_SVDS = std::min((int64_t ) (b_sz * num_matmuls / 2), n-2); - - Spectra::PartialSVDSolver svds(all_data.A_spectra, singular_triplets_found_SVDS, std::min(2 * singular_triplets_found_SVDS, n-1)); - svds.compute(); - auto stop_svds = steady_clock::now(); - dur_svds = duration_cast(stop_svds - start_svds).count(); - printf("TOTAL TIME FOR SVDS %ld\n", dur_svds); - - // Copy data from Spectra (Eigen) format to the nomal C++. - Matrix U_spectra = svds.matrix_U(singular_triplets_found_SVDS); - Matrix V_spectra = svds.matrix_V(singular_triplets_found_SVDS); - Vector S_spectra = svds.singular_values(); + residual_err_vec_ABRIK = residual_error_vectors_comp(all_data.A_linop, m, n, all_data.U, all_data.V, all_data.Sigma, singular_triplets_target_ABRIK, U_scratch, V_scratch); + printf("ABRIK residual error: %.16e\n", residual_err_vec_ABRIK); - all_data.U = new T[m * singular_triplets_found_SVDS](); - all_data.V = new T[n * singular_triplets_found_SVDS](); - all_data.Sigma = new T[m * singular_triplets_found_SVDS](); - - Eigen::Map(all_data.U, m, singular_triplets_found_SVDS) = U_spectra; - Eigen::Map(all_data.V, n, singular_triplets_found_SVDS) = V_spectra; - Eigen::Map(all_data.Sigma, singular_triplets_found_SVDS) = S_spectra; + // Write ABRIK output matrices to files if requested (only on first run) + if (write_output_matrices && i == 0) { + int64_t full_abrik_output_size = all_algs.ABRIK.singular_triplets_found; + std::string prefix = output_dir + "/ABRIK_bsz" + std::to_string(b_sz) + "_mm" + std::to_string(num_matmuls); + write_matrix_to_file(prefix + "_U.mtx", all_data.U, m, full_abrik_output_size, false); + write_matrix_to_file(prefix + "_V.mtx", all_data.V, n, full_abrik_output_size, false); + write_matrix_to_file(prefix + "_Sigma.mtx", all_data.Sigma, full_abrik_output_size, 1, true); + } - singular_triplets_target_SVDS = std::min(target_rank, singular_triplets_found_SVDS); + // Cleanup ABRIK outputs (new[]) + delete[] all_data.U; all_data.U = nullptr; + delete[] all_data.V; all_data.V = nullptr; + delete[] all_data.Sigma; all_data.Sigma = nullptr; - residual_err_custom_SVDS = residual_error_comp(all_data, singular_triplets_target_SVDS); - printf("SVDS sqrt(||AV - SU||^2_F + ||A'U - VS||^2_F) / sqrt(target_rank): %.16e\n", residual_err_custom_SVDS); - state_alg = state; - state_gen = state; - data_regen(all_data, state_gen); - std::ofstream file(output_filename, std::ios::app); - file << b_sz << ", " << all_algs.ABRIK.max_krylov_iters << ", " << target_rank << ", " - << residual_err_custom_ABRIK << ", " << dur_ABRIK << ", " - << residual_err_custom_SVDS << ", " << dur_svds << ",\n"; + // ---- SVDS (Spectra, sparse, budgeted) ---- + // Uses BudgetedPartialSVDSolver with fixed nev=target_rank and a + // matvec budget derived from (b_sz, num_matmuls) for fair comparison. + { + int64_t nev_svds = target_rank; + int64_t ncv_default = std::min(2 * nev_svds + 1, n - 1); + int64_t svds_budget = b_sz * num_matmuls; + int64_t ncv_svds = BenchmarkUtil::effective_ncv(svds_budget, nev_svds, ncv_default); + int64_t max_restarts = BenchmarkUtil::budget_to_restarts(svds_budget, nev_svds, ncv_svds); + + auto start_svds = steady_clock::now(); + BenchmarkUtil::BudgetedPartialSVDSolver svds(A_spectra_ref, nev_svds, ncv_svds); + svds.compute(max_restarts); + auto stop_svds = steady_clock::now(); + dur_svds = duration_cast(stop_svds - start_svds).count(); + printf("SVDS: budget=%ld matvecs, max_restarts=%ld, A'A_ops=%ld, time=%ld us\n", + svds_budget, max_restarts, svds.num_operations(), dur_svds); + + singular_triplets_found_SVDS = nev_svds; + + EMatrix U_spectra = svds.matrix_U(singular_triplets_found_SVDS); + EMatrix V_spectra = svds.matrix_V(singular_triplets_found_SVDS); + EVector S_spectra = svds.singular_values(); + + all_data.U = new T[m * singular_triplets_found_SVDS](); + all_data.V = new T[n * singular_triplets_found_SVDS](); + all_data.Sigma = new T[singular_triplets_found_SVDS](); + + Eigen::Map(all_data.U, m, singular_triplets_found_SVDS) = U_spectra; + Eigen::Map(all_data.V, n, singular_triplets_found_SVDS) = V_spectra; + Eigen::Map(all_data.Sigma, singular_triplets_found_SVDS) = S_spectra; + + singular_triplets_target_SVDS = std::min(target_rank, singular_triplets_found_SVDS); + + residual_err_vec_SVDS = residual_error_vectors_comp(all_data.A_linop, m, n, all_data.U, all_data.V, all_data.Sigma, singular_triplets_target_SVDS, U_scratch, V_scratch); + printf("SVDS residual error: %.16e\n", residual_err_vec_SVDS); + + // Cleanup SVDS outputs (new[]) + delete[] all_data.U; all_data.U = nullptr; + delete[] all_data.V; all_data.V = nullptr; + delete[] all_data.Sigma; all_data.Sigma = nullptr; + + state_alg = state; + } + + // Write CSV data row + outfile << b_sz << ", " << all_algs.ABRIK.max_krylov_iters << ", " << target_rank << ", " + << residual_err_vec_ABRIK << ", " << dur_ABRIK << ", " + << residual_err_vec_SVDS << ", " << dur_svds << "\n"; + outfile.flush(); } -} -int main(int argc, char *argv[]) { + delete[] U_scratch; + delete[] V_scratch; +} - if (argc < 9) { - // Expected input into this benchmark. - std::cerr << "Usage: " << argv[0] << " " << std::endl; - return 1; +template +static void run_benchmark(int argc, char *argv[]) { + + if (argc < 12) { + std::cerr << "Usage: " << argv[0] << " " << std::endl; + std::cerr << " run_gesdd: 1 to run direct GESDD, 0 to skip" << std::endl; + std::cerr << " write_matrices: 1 to write U,V,Sigma to files, 0 to skip" << std::endl; + std::cerr << " submatrix_dim_ratio: ratio of input matrix to use (e.g., 0.5 for half, 1.0 for full matrix)" << std::endl; + return; } - double submatrix_dim_ratio = 0.5; + T submatrix_dim_ratio = (T)std::stod(argv[8]); - int num_runs = std::stol(argv[3]); - int64_t target_rank = std::stol(argv[4]); + int num_runs = std::stol(argv[4]); + int64_t target_rank = std::stol(argv[5]); + bool run_gesdd = (std::stoi(argv[6]) != 0); + bool write_matrices = (std::stoi(argv[7]) != 0); std::vector b_sz; - for (int i = 0; i < std::stol(argv[5]); ++i) - b_sz.push_back(std::stoi(argv[i + 7])); - // Save elements in string for logging purposes + for (int i = 0; i < std::stol(argv[9]); ++i) + b_sz.push_back(std::stoi(argv[i + 11])); std::ostringstream oss1; for (const auto &val : b_sz) oss1 << val << ", "; std::string b_sz_string = oss1.str(); std::vector matmuls; - for (int i = 0; i < std::stol(argv[6]); ++i) - matmuls.push_back(std::stoi(argv[i + 7 + std::stol(argv[5])])); - // Save elements in string for logging purposes + for (int i = 0; i < std::stol(argv[10]); ++i) + matmuls.push_back(std::stoi(argv[i + 11 + std::stol(argv[9])])); std::ostringstream oss2; for (const auto &val : matmuls) oss2 << val << ", "; std::string matmuls_string = oss2.str(); - double tol = std::pow(std::numeric_limits::epsilon(), 0.85); + T tol = std::pow(std::numeric_limits::epsilon(), (T)0.85); auto state = RandBLAS::RNGState(); auto state_constant = state; - // Read the input fast matrix market data - // The idea is that input_mat_coo will be automatically freed at the end of function execution - auto input_mat_coo = from_matrix_market(std::string(argv[2])); - auto m = input_mat_coo.n_rows * submatrix_dim_ratio; - auto n = input_mat_coo.n_cols * submatrix_dim_ratio; + // Ensure output directory exists if we're writing matrices + if (write_matrices) { + ensure_directory_exists(std::string(argv[2])); + } + + // Read the input matrix market data + auto input_mat_coo = from_matrix_market(std::string(argv[3])); + auto m = (int64_t)(input_mat_coo.n_rows * submatrix_dim_ratio); + auto n = (int64_t)(input_mat_coo.n_cols * submatrix_dim_ratio); - // Convert coo into csc matrix - this will grab the leading principal submatrix - // depending on what m and n were set to. - auto input_mat_csc = format_conversion(m, n, input_mat_coo); + // Convert COO to CSC (grabs leading principal submatrix) + auto input_mat_csc = format_conversion(m, n, input_mat_coo); // Allocate basic workspace. - ABRIK_benchmark_data> all_data(input_mat_csc, m, n, tol); + ABRIK_benchmark_data> all_data(input_mat_csc, m, n, tol); all_data.A_input = &input_mat_csc; - // Read matrix into spectra format - //from_matrix_market(std::string(argv[2]), all_data.A_spectra); - from_matrix_market_leading_submatrix(std::string(argv[2]), all_data.A_spectra, submatrix_dim_ratio); + + // Read matrix into Spectra (Eigen sparse) format + using ESpMatrix = typename EigenTypes::SpMatrix; + ESpMatrix A_spectra(m, n); + from_matrix_market_leading_submatrix(std::string(argv[3]), A_spectra, submatrix_dim_ratio); // Declare ABRIK object - ABRIK_algorithm_objects all_algs(false, false, tol); + ABRIK_algorithm_objects all_algs(false, false, tol); printf("Finished data preparation\n"); - // Declare a data file - std::string output_filename = "_ABRIK_speed_comparisons_sparse_num_info_lines_" + std::to_string(6) + ".txt"; + + // Generate date/time prefix + std::time_t now = std::time(nullptr); + char date_prefix[20]; + std::strftime(date_prefix, sizeof(date_prefix), "%Y%m%d_%H%M%S_", std::localtime(&now)); + + // Build output file path + std::string output_filename = std::string(date_prefix) + "ABRIK_speed_comparisons_sparse.csv"; std::string path; - if (std::string(argv[1]) != ".") { - path = argv[1] + output_filename; + if (std::string(argv[2]) != ".") { + path = std::string(argv[2]) + "/" + output_filename; } else { path = output_filename; } std::ofstream file(path, std::ios::out | std::ios::app); - // Writing important data into file - file << "Description: Results from the ABRIK speed comparison benchmark, recording the time it takes to perform ABRIK and alternative methods for low-rank SVD, specifically on sparse matrices." - "\nFile format: 15 columns, showing krylov block size, nummber of matmuls permitted, and num svals and svecs to approximate, followed by the residual error, standard lowrank error and execution time for all algorithms (ABRIK, SVDS)" - "\n Rows correspond to algorithm runs with Krylov block sizes varying as specified, and numbers of matmuls varying as specified per each block size, with num_runs repititions of each number of matmuls." - "\nInput type:" + std::string(argv[2]) + - "\nInput size:" + std::to_string(m) + " by " + std::to_string(n) + - "\nAdditional parameters: Krylov block sizes " + b_sz_string + - " matmuls: " + matmuls_string + - " num runs per size " + std::to_string(num_runs) + - " num singular values and vectors approximated " + std::to_string(target_rank) + - "\n"; + // Write metadata header (prefixed with # for easy parsing) + file << "# ABRIK Speed Comparison Benchmark (Sparse)\n" + << "# Precision: " << argv[1] << "\n" + << "# Input matrix: " << argv[3] << "\n" + << "# Input size: " << m << " x " << n << " (submatrix ratio: " << submatrix_dim_ratio << ")\n" + << "# Target rank: " << target_rank << "\n" + << "# Krylov block sizes: " << b_sz_string << "\n" + << "# Matmul counts: " << matmuls_string << "\n" + << "# Runs per configuration: " << num_runs << "\n" + << "# Tolerance: " << tol << "\n" + << "# Run GESDD: " << (run_gesdd ? "yes" : "no") << "\n" + << "# Write matrices: " << (write_matrices ? "yes" : "no") << "\n" + << "# Residual metric: sqrt(||S^{-1}AV - U||^2_F + ||(A'U)S^{-1} - V||^2_F)\n" + << "# Timings in microseconds\n"; + // Write CSV column header + file << "b_sz, num_matmuls, target_rank, " + << "err_ABRIK, dur_ABRIK, " + << "err_SVDS, dur_SVDS\n"; file.flush(); + // Run direct GESDD if requested + if (run_gesdd) { + compute_direct_gesdd(A_spectra, m, n, target_rank, std::string(argv[2]), write_matrices); + } + + auto start_total = steady_clock::now(); + size_t i = 0, j = 0; - for (;i < b_sz.size(); ++i) { - for (;j < matmuls.size(); ++j) { - call_all_algs(num_runs, b_sz[i], matmuls[j], target_rank, all_algs, all_data, state_constant, path, std::string(argv[2])); + for (; i < b_sz.size(); ++i) { + for (; j < matmuls.size(); ++j) { + call_all_algs(num_runs, b_sz[i], matmuls[j], target_rank, all_algs, all_data, A_spectra, state_constant, file, std::string(argv[3]), std::string(argv[2]), write_matrices); } j = 0; } + + auto stop_total = steady_clock::now(); + long total_us = duration_cast(stop_total - start_total).count(); + printf("\nTOTAL BENCHMARK TIME: %.2f seconds\n", total_us / 1e6); +} + +int main(int argc, char *argv[]) { + if (argc < 2) { + std::cerr << "Usage: " << argv[0] << " " << std::endl; + return 1; + } + + std::string precision = argv[1]; + if (precision == "double") { + run_benchmark(argc, argv); + } else if (precision == "float") { + run_benchmark(argc, argv); + } else { + std::cerr << "Error: precision must be 'double' or 'float', got '" << precision << "'" << std::endl; + return 1; + } + return 0; } diff --git a/benchmark/bench_ABRIK/BudgetedSVDSolver.hh b/benchmark/bench_ABRIK/BudgetedSVDSolver.hh new file mode 100644 index 000000000..78428eab7 --- /dev/null +++ b/benchmark/bench_ABRIK/BudgetedSVDSolver.hh @@ -0,0 +1,216 @@ +/* +BudgetedSVDSolver — a wrapper around Spectra's PartialSVDSolver that exposes +a "max_restarts" parameter and returns ALL Ritz approximations (not just +converged ones). + +This allows benchmarking SVDS with a controlled computational budget, +making it directly comparable to ABRIK's (b_sz, num_matmuls) parameterization. + +Algorithm: Implicitly Restarted Lanczos (single-vector) for eigenvalues of A'A. + - nev eigenvalues requested (= target_rank) + - ncv Krylov subspace dimension (= 2*nev + 1) + - max_restarts controls the number of restart iterations + - Total A'A operations ≈ ncv + max_restarts * (ncv - nev) + - Total matvecs with A ≈ 2 * (total A'A operations) +*/ + +#ifndef BUDGETED_SVD_SOLVER_HH +#define BUDGETED_SVD_SOLVER_HH + +#include +#include +#include +#include +#include +#include + +namespace BenchmarkUtil { + +// Subclass of SymEigsSolver that can extract ALL Ritz pairs (even unconverged) +// by recomputing eigenvectors from the tridiagonal matrix H. +template +class AllRitzSymEigsSolver : public Spectra::SymEigsSolver +{ +public: + using Base = Spectra::SymEigsSolver; + using Scalar = typename OpType::Scalar; + using Index = Eigen::Index; + using RealScalar = typename Eigen::NumTraits::Real; + using Matrix = Eigen::Matrix; + using RealMatrix = Eigen::Matrix; + using RealVector = Eigen::Matrix; + + AllRitzSymEigsSolver(OpType& op, Index nev, Index ncv) : + Base(op, nev, ncv) {} + + // Return ALL nev Ritz values (eigenvalue approximations), sorted largest first. + // These are available after compute() regardless of convergence status. + RealVector all_eigenvalues() const + { + // m_ritz_val is protected, contains ncv Ritz values sorted by selection rule + // First nev entries are the "wanted" ones (largest for SVD). + return this->m_ritz_val.head(this->m_nev); + } + + // Return ALL nev Ritz vectors in the original space, sorted to match all_eigenvalues(). + // Recomputes eigenvectors of H from scratch (O(ncv^3), negligible vs matvecs). + Matrix all_eigenvectors() const + { + // Get the tridiagonal matrix H from the Lanczos factorization + Spectra::TridiagEigen decomp(this->m_fac.matrix_H().real()); + const RealVector& evals = decomp.eigenvalues(); + const RealMatrix& evecs = decomp.eigenvectors(); + + // Sort eigenvalues largest first (matching all_eigenvalues() order) + std::vector ind(evals.size()); + std::iota(ind.begin(), ind.end(), 0); + std::sort(ind.begin(), ind.end(), + [&evals](Index a, Index b) { return evals[a] > evals[b]; }); + + // Extract first nev eigenvectors in sorted order + Index nev = this->m_nev; + Index ncv = this->m_ncv; + RealMatrix ritz_vec(ncv, nev); + for (Index i = 0; i < nev; i++) + ritz_vec.col(i) = evecs.col(ind[i]); + + // Transform to original space: V * ritz_vec + // V is the Krylov basis (n x ncv), ritz_vec is (ncv x nev) + return this->m_fac.matrix_V() * ritz_vec; + } +}; + +// Budgeted Partial SVD Solver. +// Like Spectra::PartialSVDSolver but with explicit max_restarts control +// and access to ALL Ritz approximations. +template > +class BudgetedPartialSVDSolver +{ +public: + using Scalar = typename MatrixType::Scalar; + using Index = Eigen::Index; + using Matrix = Eigen::Matrix; + using Vector = Eigen::Matrix; + using ConstGenericMatrix = const Eigen::Ref; + +private: + ConstGenericMatrix m_mat; + const Index m_m; + const Index m_n; + Spectra::SVDMatOp* m_op; + AllRitzSymEigsSolver>* m_eigs; + Index m_nconv; + Index m_nev; + +public: + BudgetedPartialSVDSolver(ConstGenericMatrix& mat, Index ncomp, Index ncv) : + m_mat(mat), m_m(mat.rows()), m_n(mat.cols()), m_nev(ncomp) + { + if (m_m > m_n) + m_op = new Spectra::SVDTallMatOp(mat); + else + m_op = new Spectra::SVDWideMatOp(mat); + + m_eigs = new AllRitzSymEigsSolver>(*m_op, ncomp, ncv); + } + + ~BudgetedPartialSVDSolver() + { + delete m_eigs; + delete m_op; + } + + // Run with explicit max_restarts. + // Uses tiny tolerance to force all restart iterations. + // Returns number of (formally) converged eigenvalues. + Index compute(Index max_restarts) + { + m_eigs->init(); + // Use extremely small tol to prevent early stopping + m_nconv = m_eigs->compute(Spectra::SortRule::LargestAlge, max_restarts, 1e-100); + return m_nconv; + } + + // Number of matrix operations performed (A'A applications) + Index num_operations() const { return m_eigs->num_operations(); } + + // ALL singular value approximations (nev values, even unconverged) + Vector singular_values() const + { + Vector evals = m_eigs->all_eigenvalues(); + // Clamp negative eigenvalues to small positive (numerical noise) + for (Index i = 0; i < evals.size(); i++) + evals[i] = std::sqrt(std::max(evals[i], Scalar(0))); + return evals; + } + + // ALL left singular vector approximations + Matrix matrix_U(Index nu) + { + nu = (std::min)(nu, m_nev); + Matrix evecs = m_eigs->all_eigenvectors().leftCols(nu); + Vector evals = m_eigs->all_eigenvalues().head(nu); + + if (m_m <= m_n) + return evecs; + + // U = A * V * diag(1/sigma) + Matrix result = m_mat * evecs; + for (Index i = 0; i < nu; i++) { + Scalar sigma = std::sqrt(std::max(evals[i], Scalar(0))); + if (sigma > 0) + result.col(i) /= sigma; + } + return result; + } + + // ALL right singular vector approximations + Matrix matrix_V(Index nv) + { + nv = (std::min)(nv, m_nev); + Matrix evecs = m_eigs->all_eigenvectors().leftCols(nv); + Vector evals = m_eigs->all_eigenvalues().head(nv); + + if (m_m > m_n) + return evecs; + + // V = A' * U * diag(1/sigma) + Matrix result = m_mat.transpose() * evecs; + for (Index i = 0; i < nv; i++) { + Scalar sigma = std::sqrt(std::max(evals[i], Scalar(0))); + if (sigma > 0) + result.col(i) /= sigma; + } + return result; + } +}; + +// Compute effective ncv for a given matvec budget. +// When budget is small, reduce ncv so Spectra builds a smaller Krylov subspace +// (less work, lower accuracy). ncv must satisfy nev < ncv <= budget/2. +inline int64_t effective_ncv(int64_t budget, int64_t nev, int64_t ncv_default) +{ + int64_t ata_ops = budget / 2; + // If budget allows full ncv, use it + if (ata_ops >= ncv_default) + return ncv_default; + // Otherwise shrink ncv to fit budget, but keep ncv > nev + return std::max(nev + 1, ata_ops); +} + +// Convert a matvec budget to max_restarts for the Lanczos solver. +// budget = total matvecs with A (comparable to ABRIK's b_sz * num_matmuls) +// Each A'A application = 2 matvecs with A. +// Initial Lanczos: ncv A'A ops. Each restart: ~(ncv - nev) A'A ops. +// ncv should be the effective ncv (from effective_ncv()). +inline int64_t budget_to_restarts(int64_t budget, int64_t nev, int64_t ncv) +{ + int64_t ata_ops = budget / 2; // A'A operations from matvec budget + if (ata_ops <= ncv) + return 0; // not enough budget for even one restart — just initial factorization + return (ata_ops - ncv) / (ncv - nev); +} + +} // namespace BenchmarkUtil + +#endif // BUDGETED_SVD_SOLVER_HH diff --git a/test/comps/test_qb.cc b/test/comps/test_qb.cc index 264ee42d5..a77c7e4c6 100644 --- a/test/comps/test_qb.cc +++ b/test/comps/test_qb.cc @@ -110,7 +110,6 @@ class TestQB : public ::testing::Test auto n = all_data.col; auto k = all_data.rank; - T* A_dat = all_data.A.data(); T* A_hat_dat = all_data.A_hat.data(); T* A_k_dat = all_data.A_k.data(); @@ -119,10 +118,14 @@ class TestQB : public ::testing::Test T* S_dat = all_data.S.data(); T* VT_dat = all_data.VT.data(); + // Save a copy of A before QB call (QB now modifies A in-place via deflation). + T* A_dat = new T[m * n]; + lapack::lacpy(MatrixType::General, m, n, all_data.A.data(), m, A_dat, m); + T* Q = nullptr; T* BT = nullptr; - // Regular QB2 call + // Regular QB2 call — NOTE: this modifies all_data.A in-place. all_algs.QB.call(m, n, all_data.A.data(), k, block_sz, tol, Q, BT, state); // Reassing pointers because Q, B have been resized @@ -168,6 +171,7 @@ class TestQB : public ::testing::Test T norm_test_4 = lapack::lange(Norm::Fro, m, n, A_hat_dat, m); printf("FRO NORM OF A_k - QB: %e\n", norm_test_4); ASSERT_NEAR(norm_test_4, 0, test_tol); + delete[] A_dat; free(Q); free(BT); } @@ -189,15 +193,18 @@ class TestQB : public ::testing::Test int64_t k_est = std::min(m, n); - T* A_dat = all_data.A.data(); T* Q_dat = all_data.Q.data(); T* BT_dat = all_data.BT.data(); T* A_hat_dat = all_data.A_hat.data(); + // Save a copy of A before QB call (QB now modifies A in-place via deflation). + T* A_dat = new T[m * n]; + lapack::lacpy(MatrixType::General, m, n, all_data.A.data(), m, A_dat, m); + T* Q = nullptr; T* BT = nullptr; - // Regular QB2 call + // Regular QB2 call — NOTE: this modifies all_data.A in-place. all_algs.QB.call(m, n, all_data.A.data(), k_est, block_sz, tol, Q, BT, state); // Reassing pointers because Q, B have been resized @@ -224,6 +231,7 @@ class TestQB : public ::testing::Test printf("FRO NORM OF A: %e\n", norm_A); EXPECT_TRUE(norm_test_1 <= (tol * norm_A)); } + delete[] A_dat; free(Q); free(BT); } diff --git a/test/drivers/test_abrik.cc b/test/drivers/test_abrik.cc index 36b1181b0..70885179e 100644 --- a/test/drivers/test_abrik.cc +++ b/test/drivers/test_abrik.cc @@ -158,8 +158,8 @@ TEST_F(TestABRIK, ABRIK_basic1) { ABRIKTestData all_data(m, n); RandLAPACK::ABRIK ABRIK(false, false, tol); - ABRIK.num_threads_max = RandLAPACK::util::get_omp_threads(); - ABRIK.num_threads_min = 1; + + RandLAPACK::gen::mat_gen_info m_info(m, n, RandLAPACK::gen::gaussian); RandLAPACK::gen::mat_gen(m_info, all_data.A, state); @@ -179,8 +179,8 @@ TEST_F(TestABRIK, ABRIK_basic) { ABRIKTestData all_data(m, n); RandLAPACK::ABRIK ABRIK(false, false, tol); - ABRIK.num_threads_max = RandLAPACK::util::get_omp_threads(); - ABRIK.num_threads_min = 1; + + RandLAPACK::gen::mat_gen_info m_info(m, n, RandLAPACK::gen::gaussian); RandLAPACK::gen::mat_gen(m_info, all_data.A, state); @@ -200,8 +200,8 @@ TEST_F(TestABRIK, ABRIK_sparse_csc) { ABRIKTestDataSparse> all_data(m, n); RandLAPACK::ABRIK ABRIK(false, false, tol); - ABRIK.num_threads_min = 1; - ABRIK.num_threads_max = RandLAPACK::util::get_omp_threads(); + + RandLAPACK::gen::mat_gen_info m_info(m, n, RandLAPACK::gen::gaussian); test::test_datastructures::test_spmats::iid_sparsify_random_dense(m, n, Layout::ColMajor, all_data.A_buff, 0.9, 0); @@ -221,8 +221,8 @@ TEST_F(TestABRIK, ABRIK_sparse_csr) { ABRIKTestDataSparse> all_data(m, n); RandLAPACK::ABRIK ABRIK(false, false, tol); - ABRIK.num_threads_min = 1; - ABRIK.num_threads_max = RandLAPACK::util::get_omp_threads(); + + RandLAPACK::gen::mat_gen_info m_info(m, n, RandLAPACK::gen::gaussian); test::test_datastructures::test_spmats::iid_sparsify_random_dense(m, n, Layout::ColMajor, all_data.A_buff, 0.9, 0); @@ -242,8 +242,8 @@ TEST_F(TestABRIK, ABRIK_sparse_coo) { ABRIKTestDataSparse> all_data(m, n); RandLAPACK::ABRIK ABRIK(false, false, tol); - ABRIK.num_threads_min = 1; - ABRIK.num_threads_max = RandLAPACK::util::get_omp_threads(); + + RandLAPACK::gen::mat_gen_info m_info(m, n, RandLAPACK::gen::gaussian); test::test_datastructures::test_spmats::iid_sparsify_random_dense(m, n, Layout::ColMajor, all_data.A_buff, 0.9, 0); @@ -263,13 +263,177 @@ TEST_F(TestABRIK, ABRIK_sparse_coo_cqrrt) { ABRIKTestDataSparse> all_data(m, n); RandLAPACK::ABRIK ABRIK(false, false, tol); - ABRIK.num_threads_min = 1; + ABRIK.qr_exp = Subroutines::QR_explicit::cqrrt; - ABRIK.num_threads_max = RandLAPACK::util::get_omp_threads(); + RandLAPACK::gen::mat_gen_info m_info(m, n, RandLAPACK::gen::gaussian); test::test_datastructures::test_spmats::iid_sparsify_random_dense(m, n, Layout::ColMajor, all_data.A_buff, 0.9, 0); RandBLAS::sparse_data::coo::dense_to_coo(Layout::ColMajor, all_data.A_buff, 0.0, all_data.A); test_ABRIK_general(b_sz, target_rank, custom_rank, all_data, ABRIK, state); +} + +// ========== Adaptive mode tests ========== + +// Adaptive mode converges from a small initial max_krylov_iters. +TEST_F(TestABRIK, ABRIK_adaptive_converges) { + int64_t m = 200; + int64_t n = 100; + int64_t b_sz = 10; + double tol = std::pow(std::numeric_limits::epsilon(), 0.85); + auto state = RandBLAS::RNGState(); + + ABRIKTestData all_data(m, n); + RandLAPACK::ABRIK ABRIK(false, false, tol); + ABRIK.adaptive = true; + ABRIK.max_krylov_iters = 4; // Start with few iterations + + RandLAPACK::gen::mat_gen_info m_info(m, n, RandLAPACK::gen::gaussian); + RandLAPACK::gen::mat_gen(m_info, all_data.A, state); + lapack::lacpy(MatrixType::General, m, n, all_data.A, m, all_data.A_buff, m); + + ABRIK.call(m, n, all_data.A, m, b_sz, all_data.U, all_data.V, all_data.Sigma, state); + + auto k = ABRIK.singular_triplets_found; + double residual = residual_error_comp(all_data, k); + printf("adaptive_converges: residual %e, k=%ld, iters=%d\n", residual, k, ABRIK.num_krylov_iters); + ASSERT_LE(residual, 10 * std::pow(std::numeric_limits::epsilon(), 0.825)); + ASSERT_GT(ABRIK.num_krylov_iters, 4); // Should have extended beyond initial +} + +// Adaptive mode with unreasonable tolerance — BK norm converges, ABRIK stops gracefully. +TEST_F(TestABRIK, ABRIK_adaptive_norm_converged) { + int64_t m = 200; + int64_t n = 100; + int64_t b_sz = 10; + double tol = 1e-20; // Unreachable in double precision + auto state = RandBLAS::RNGState(); + + ABRIKTestData all_data(m, n); + RandLAPACK::ABRIK ABRIK(false, false, tol); + ABRIK.adaptive = true; + ABRIK.max_krylov_iters = 4; + + RandLAPACK::gen::mat_gen_info m_info(m, n, RandLAPACK::gen::gaussian); + RandLAPACK::gen::mat_gen(m_info, all_data.A, state); + lapack::lacpy(MatrixType::General, m, n, all_data.A, m, all_data.A_buff, m); + + ABRIK.call(m, n, all_data.A, m, b_sz, all_data.U, all_data.V, all_data.Sigma, state); + + // Should terminate gracefully despite unreasonable tolerance. + auto k = ABRIK.singular_triplets_found; + printf("adaptive_norm_converged: iters=%d, k=%ld\n", ABRIK.num_krylov_iters, k); + ASSERT_GT(k, (int64_t)0); + // Result should still be reasonable even though tol wasn't met. + double residual = residual_error_comp(all_data, std::min(k, (int64_t)50)); + printf("adaptive_norm_converged: residual %e\n", residual); + ASSERT_LE(residual, 10 * std::pow(std::numeric_limits::epsilon(), 0.825)); +} + +// Adaptive mode with a rank-deficient matrix — BK detects rank deficiency, ABRIK stops. +TEST_F(TestABRIK, ABRIK_adaptive_rank_deficient) { + int64_t m = 100; + int64_t n = 50; + int64_t b_sz = 10; + int64_t true_rank = 5; + double tol = 1e-20; // Unreachable + auto state = RandBLAS::RNGState(); + + ABRIKTestData all_data(m, n); + RandLAPACK::ABRIK ABRIK(false, false, tol); + ABRIK.adaptive = true; + ABRIK.max_krylov_iters = 4; + + // Create a rank-5 matrix: A = L * R + double* L = new double[m * true_rank](); + double* R_mat = new double[true_rank * n](); + RandBLAS::DenseDist DL(m, true_rank); + state = RandBLAS::fill_dense(DL, L, state); + RandBLAS::DenseDist DR(true_rank, n); + state = RandBLAS::fill_dense(DR, R_mat, state); + blas::gemm(Layout::ColMajor, Op::NoTrans, Op::NoTrans, m, n, true_rank, + 1.0, L, m, R_mat, true_rank, 0.0, all_data.A, m); + lapack::lacpy(MatrixType::General, m, n, all_data.A, m, all_data.A_buff, m); + delete[] L; + delete[] R_mat; + + ABRIK.call(m, n, all_data.A, m, b_sz, all_data.U, all_data.V, all_data.Sigma, state); + + auto k = ABRIK.singular_triplets_found; + printf("adaptive_rank_deficient: iters=%d, k=%ld\n", ABRIK.num_krylov_iters, k); + ASSERT_GT(k, (int64_t)0); +} + +// Adaptive mode with max_retries=1 — verifies the retry limit is respected. +TEST_F(TestABRIK, ABRIK_adaptive_max_retries) { + int64_t m = 200; + int64_t n = 100; + int64_t b_sz = 10; + double tol = 1e-20; // Unreachable + auto state = RandBLAS::RNGState(); + + ABRIKTestData all_data(m, n); + RandLAPACK::ABRIK ABRIK(false, false, tol); + ABRIK.adaptive = true; + ABRIK.max_krylov_iters = 4; + ABRIK.adaptive_max_retries = 1; + + RandLAPACK::gen::mat_gen_info m_info(m, n, RandLAPACK::gen::gaussian); + RandLAPACK::gen::mat_gen(m_info, all_data.A, state); + lapack::lacpy(MatrixType::General, m, n, all_data.A, m, all_data.A_buff, m); + + ABRIK.call(m, n, all_data.A, m, b_sz, all_data.U, all_data.V, all_data.Sigma, state); + + printf("adaptive_max_retries: iters=%d, k=%ld\n", ABRIK.num_krylov_iters, ABRIK.singular_triplets_found); + // Initial call: 4 iters. After 1 retry: 4 more iters = 8 total. + ASSERT_GT(ABRIK.num_krylov_iters, 4); + ASSERT_LE(ABRIK.num_krylov_iters, 8); + ASSERT_GT(ABRIK.singular_triplets_found, (int64_t)0); +} + +// Adaptive mode produces comparable quality to non-adaptive with enough iterations. +TEST_F(TestABRIK, ABRIK_adaptive_matches_nonadaptive) { + int64_t m = 200; + int64_t n = 100; + int64_t b_sz = 10; + double tol = std::pow(std::numeric_limits::epsilon(), 0.85); + + // Generate the matrix once. + ABRIKTestData data1(m, n); + auto state = RandBLAS::RNGState(); + RandLAPACK::gen::mat_gen_info m_info(m, n, RandLAPACK::gen::gaussian); + RandLAPACK::gen::mat_gen(m_info, data1.A, state); + lapack::lacpy(MatrixType::General, m, n, data1.A, m, data1.A_buff, m); + + // Copy for second run. + ABRIKTestData data2(m, n); + lapack::lacpy(MatrixType::General, m, n, data1.A_buff, m, data2.A, m); + lapack::lacpy(MatrixType::General, m, n, data1.A_buff, m, data2.A_buff, m); + + // Run 1: non-adaptive with generous iterations. + auto state1 = RandBLAS::RNGState(); + RandLAPACK::ABRIK ABRIK1(false, false, tol); + ABRIK1.max_krylov_iters = 20; + ABRIK1.call(m, n, data1.A, m, b_sz, data1.U, data1.V, data1.Sigma, state1); + + auto k1 = ABRIK1.singular_triplets_found; + double residual1 = residual_error_comp(data1, std::min(k1, (int64_t)50)); + + // Run 2: adaptive with small initial iterations. + auto state2 = RandBLAS::RNGState(); + RandLAPACK::ABRIK ABRIK2(false, false, tol); + ABRIK2.adaptive = true; + ABRIK2.max_krylov_iters = 4; + ABRIK2.call(m, n, data2.A, m, b_sz, data2.U, data2.V, data2.Sigma, state2); + + auto k2 = ABRIK2.singular_triplets_found; + double residual2 = residual_error_comp(data2, std::min(k2, (int64_t)50)); + + printf("non-adaptive: residual %e, k=%ld, iters=%d\n", residual1, k1, ABRIK1.num_krylov_iters); + printf("adaptive: residual %e, k=%ld, iters=%d\n", residual2, k2, ABRIK2.num_krylov_iters); + + // Both should achieve good quality. + ASSERT_LE(residual1, 10 * std::pow(std::numeric_limits::epsilon(), 0.825)); + ASSERT_LE(residual2, 10 * std::pow(std::numeric_limits::epsilon(), 0.825)); } \ No newline at end of file