diff --git a/RandLAPACK/drivers/rl_bqrrp.hh b/RandLAPACK/drivers/rl_bqrrp.hh index 743dbf1e..3f65c607 100644 --- a/RandLAPACK/drivers/rl_bqrrp.hh +++ b/RandLAPACK/drivers/rl_bqrrp.hh @@ -319,7 +319,8 @@ int BQRRP::call( internal_nb = std::min(internal_nb, b_sz); block_rank = b_sz; - // Zero-out data - may not be necessary + // Zero-out data. The J_buffer fill is required: geqp3 reads jpvt on + // entry and treats any nonzero entry as a fixed column. std::fill(&J_buffer[0], &J_buffer[n], 0); std::fill(&J_buffer_lu[0], &J_buffer_lu[std::min(d, n)], 0); std::fill(&Work2[0], &Work2[n], (T) 0.0); diff --git a/RandLAPACK/drivers/rl_cqrrpt.hh b/RandLAPACK/drivers/rl_cqrrpt.hh index d365732b..53de315e 100644 --- a/RandLAPACK/drivers/rl_cqrrpt.hh +++ b/RandLAPACK/drivers/rl_cqrrpt.hh @@ -12,6 +12,7 @@ #include #include #include +#include using namespace std::chrono; @@ -207,9 +208,13 @@ int CQRRPT::call( // Buffer for column pivoting. int64_t* J_buf = new int64_t[n](); + // geqp3 reads jpvt on entry (nonzero marks a fixed column), so zero the + // caller's J to keep its prior contents from steering the pivoting. + std::fill(J, J + n, (int64_t) 0); + if(this -> timing) saso_t_start = steady_clock::now(); - + /// Generating a SASO RandBLAS::SparseDist DS(d, m, this->nnz); RandBLAS::SparseSkOp S(DS, state); diff --git a/RandLAPACK/drivers/rl_cqrrpt_gpu.hh b/RandLAPACK/drivers/rl_cqrrpt_gpu.hh index 33c2703d..3b742b74 100644 --- a/RandLAPACK/drivers/rl_cqrrpt_gpu.hh +++ b/RandLAPACK/drivers/rl_cqrrpt_gpu.hh @@ -15,6 +15,7 @@ #include #include #include +#include using namespace std::chrono; @@ -195,6 +196,10 @@ int CQRRPT_GPU::call( // Buffer for column pivoting. int64_t* J_buf = new int64_t[n](); + // geqp3 reads jpvt on entry (nonzero marks a fixed column), so zero the + // caller's J to keep its prior contents from steering the pivoting. + std::fill(J, J + n, (int64_t) 0); + if(this -> timing) saso_t_start = steady_clock::now(); /***********************************************************************************/ diff --git a/RandLAPACK/linops/rl_materialize.hh b/RandLAPACK/linops/rl_materialize.hh index e640f9d2..720d7e6e 100644 --- a/RandLAPACK/linops/rl_materialize.hh +++ b/RandLAPACK/linops/rl_materialize.hh @@ -5,12 +5,14 @@ // // Generic fallback: multiply by the identity matrix. // Overloaded specializations avoid that cost when the underlying storage is -// directly accessible (DenseLinOp, SparseLinOp). +// directly accessible (DenseLinOp, SparseLinOp) or reachable through the +// operands (CompositeOperator). #include "rl_exceptions.hh" #include "rl_concepts.hh" #include "rl_dense_linop.hh" #include "rl_sparse_linop.hh" +#include "rl_composite_linop.hh" #include "rl_blaspp.hh" #include "rl_lapackpp.hh" #include "rl_util.hh" @@ -20,6 +22,12 @@ namespace RandLAPACK { +/// Cap on the dimension the generic materialize fallback will accept. The +/// fallback allocates an n by n identity (2 GiB in double at the cap) behind +/// the caller's buffer. The type-specific overloads below form no identity +/// and carry no cap. +inline constexpr int64_t MATERIALIZE_IDENTITY_MAX_DIM = 16384; + /// Materialize a linear operator into a dense column-major buffer. /// /// Generic fallback: forms buf = A * I by applying the operator to the @@ -33,6 +41,11 @@ namespace RandLAPACK { template void materialize(LinOp& A, int64_t m, int64_t n, typename LinOp::scalar_t* buf, int64_t ldb) { using T = typename LinOp::scalar_t; + randlapack_require(n <= MATERIALIZE_IDENTITY_MAX_DIM) << "n=" << n + << " exceeds MATERIALIZE_IDENTITY_MAX_DIM=" << MATERIALIZE_IDENTITY_MAX_DIM + << ": the generic materialize fallback would allocate an n*n identity (" + << (n * n * (int64_t) sizeof(T)) / (1024 * 1024) << " MiB)." + << " Use a type-specific materialize overload instead"; randlapack_require(ldb >= m) << "ldb=" << ldb << " < m=" << m << " (ldb must be >= m)"; // Zero the output buffer. for (int64_t j = 0; j < n; ++j) @@ -88,4 +101,25 @@ void materialize(linops::SparseLinOp& A, int64_t m, int64_t n, } } +/// Specialization for CompositeOperator: materialize each operand, then form +/// buf = L * R with one gemm. Builds no identity, carries no cap, and never +/// calls the composite's operator(); nested composites recurse. +template +void materialize(linops::CompositeOperator& A, int64_t m, int64_t n, + typename LinOp1::scalar_t* buf, int64_t ldb) { + using T = typename LinOp1::scalar_t; + randlapack_require(m == A.n_rows) << "m=" << m << " must equal A.n_rows=" << A.n_rows << " for materialize specialization"; + randlapack_require(n == A.n_cols) << "n=" << n << " must equal A.n_cols=" << A.n_cols << " for materialize specialization"; + randlapack_require(ldb >= m) << "ldb=" << ldb << " < m=" << m << " (ldb must be >= m)"; + int64_t k = A.left_op.n_cols; + T* L = new T[m * k](); + T* R = new T[k * n](); + materialize(A.left_op, m, k, L, m); + materialize(A.right_op, k, n, R, k); + blas::gemm(blas::Layout::ColMajor, blas::Op::NoTrans, blas::Op::NoTrans, + m, n, k, (T)1.0, L, m, R, k, (T)0.0, buf, ldb); + delete[] L; + delete[] R; +} + } // end namespace RandLAPACK diff --git a/RandLAPACK/testing/rl_test_utils.hh b/RandLAPACK/testing/rl_test_utils.hh index b20ecad8..41e6e5cf 100644 --- a/RandLAPACK/testing/rl_test_utils.hh +++ b/RandLAPACK/testing/rl_test_utils.hh @@ -173,20 +173,14 @@ void initialize_test_buffers(::std::vector& C_test, ::std::vector& C_refer // Linear Operator Materialization and Analysis // ============================================================================ -/// Materialize a linear operator into a dense column-major matrix. -/// Computes A_dense = A_linop * I by applying the operator to the identity. +/// Materialize a linear operator into a dense column-major matrix via +/// RandLAPACK::materialize, which dispatches on the operator type. DenseLinOp, +/// SparseLinOp and CompositeOperator materialize directly from their storage; +/// anything else goes through the generic fallback, A_dense = A_linop * I. /// -/// WARNING: This function uses operator() to materialize, which makes it -/// unsuitable for testing operator() itself (circular reasoning). For -/// correctness tests of individual linop types, prefer type-specific -/// materialization: -/// - DenseLinOp: copy A_buff directly -/// - SparseLinOp: use RandLAPACK::util::sparse_to_dense on A_sp -/// - CompositeOperator: materialize each operand independently, then -/// compute the product with blas::gemm -/// This function is appropriate for tests that assume operator() is correct -/// and need a dense representation for some other purpose (e.g., computing -/// singular values, comparing block views against the full operator). +/// WARNING: the generic fallback materializes through operator(), so for such +/// types this function cannot be used to test operator() itself (circular +/// reasoning); use it only where operator() is assumed correct. template void materialize_linop(LinOp& A_linop, T* A_dense) { int64_t m = A_linop.n_rows; diff --git a/test/drivers/test_cqrrpt.cc b/test/drivers/test_cqrrpt.cc index 7aebcd9a..0149799e 100644 --- a/test/drivers/test_cqrrpt.cc +++ b/test/drivers/test_cqrrpt.cc @@ -251,6 +251,55 @@ TEST_F(TestCQRRPT, CQRRPT_low_rank_with_bqrrp) { norm_and_copy_computational_helper(norm_A, all_data); test_CQRRPT_general(d_factor, norm_A, all_data, CQRRPT, state); } + +// geqp3 reads jpvt on entry (nonzero marks a fixed column), so the prior +// contents of the caller's J buffer must not influence pivoting. Run the same +// rank-deficient factorization with a clean J and with two dirtied J buffers, +// at the same RNG state each time; rank and factorization quality must match. +TEST_F(TestCQRRPT, CQRRPT_dirty_J_rank_deficient) { + int64_t m = 2000; + int64_t n = 50; + int64_t k = 40; + double d_factor = 2; + double tol = std::pow(std::numeric_limits::epsilon(), 0.85); + + // Generate the rank-deficient input once. + auto gen_state = RandBLAS::RNGState(); + std::vector A_orig(m * n, 0.0); + RandLAPACK::gen::mat_gen_info m_info(m, n, RandLAPACK::gen::polynomial); + m_info.cond_num = 2; + m_info.rank = k; + m_info.exponent = 2.0; + RandLAPACK::gen::mat_gen(m_info, A_orig.data(), gen_state); + + int64_t rank_ref = -1; + for (int trial = 0; trial < 3; ++trial) { + CQRRPTTestData all_data(m, n, k); + lapack::lacpy(MatrixType::General, m, n, A_orig.data(), m, all_data.A.data(), m); + + // Trial 0 keeps the zero-initialized J; trials 1 and 2 dirty it with + // different nonzero garbage. + if (trial > 0) { + for (int64_t i = 0; i < n; ++i) + all_data.J[i] = 1 + ((7919 * trial + 31 * i) % n); + } + + RandLAPACK::CQRRPT CQRRPT(false, tol); + CQRRPT.nnz = 2; + CQRRPT.qrcp = Subroutines::QRCP::geqp3; + + double norm_A = 0; + norm_and_copy_computational_helper(norm_A, all_data); + // Same RNG state in every trial, so any difference is due to J alone. + auto state = RandBLAS::RNGState(); + test_CQRRPT_general(d_factor, norm_A, all_data, CQRRPT, state); + + if (trial == 0) + rank_ref = all_data.rank; + ASSERT_EQ(all_data.rank, rank_ref); + } +} + // Using L2 norm rank estimation here is similar to using raive estimation. // Fro norm underestimates rank even worse. TEST_F(TestCQRRPT, CQRRPT_bad_orth) { diff --git a/test/drivers/test_cqrrpt_gpu.cu b/test/drivers/test_cqrrpt_gpu.cu index 1c145b24..3820e25b 100644 --- a/test/drivers/test_cqrrpt_gpu.cu +++ b/test/drivers/test_cqrrpt_gpu.cu @@ -158,3 +158,50 @@ TEST_F(TestCQRRPT, CQRRPT_GPU_full_rank_no_hqrrp) { norm_and_copy_computational_helper(norm_A, all_data); test_CQRRPT_general>(d_factor, norm_A, all_data, CQRRPT_GPU, state); } + +// GPU counterpart of CQRRPT_dirty_J_rank_deficient in test_cqrrpt.cc (the QRCP +// itself runs on the host): dirty J buffers must not influence pivoting. +TEST_F(TestCQRRPT, CQRRPT_GPU_dirty_J_rank_deficient) { + int64_t m = 2000; + int64_t n = 50; + int64_t k = 40; + double d_factor = 2.0; + double tol = std::pow(std::numeric_limits::epsilon(), 0.85); + + // Generate the rank-deficient input once. + auto gen_state = RandBLAS::RNGState(); + std::vector A_orig(m * n, 0.0); + RandLAPACK::gen::mat_gen_info m_info(m, n, RandLAPACK::gen::polynomial); + m_info.cond_num = 2; + m_info.rank = k; + m_info.exponent = 2.0; + RandLAPACK::gen::mat_gen(m_info, A_orig.data(), gen_state); + + int64_t rank_ref = -1; + for (int trial = 0; trial < 3; ++trial) { + CQRRPTTestData all_data(m, n, k); + lapack::lacpy(MatrixType::General, m, n, A_orig.data(), m, all_data.A.data(), m); + + // Trial 0 keeps the zero-initialized J; trials 1 and 2 dirty it with + // different nonzero garbage. + if (trial > 0) { + for (int64_t i = 0; i < n; ++i) + all_data.J[i] = 1 + ((7919 * trial + 31 * i) % n); + } + + RandLAPACK::CQRRPT_GPU CQRRPT_GPU(false, false, tol); + CQRRPT_GPU.nnz = 2; + CQRRPT_GPU.num_threads = 4; + CQRRPT_GPU.no_hqrrp = 1; + + double norm_A = 0; + norm_and_copy_computational_helper(norm_A, all_data); + // Same RNG state in every trial, so any difference is due to J alone. + auto state = RandBLAS::RNGState(); + test_CQRRPT_general>(d_factor, norm_A, all_data, CQRRPT_GPU, state); + + if (trial == 0) + rank_ref = all_data.rank; + ASSERT_EQ(all_data.rank, rank_ref); + } +} diff --git a/test/linops/test_linops.cc b/test/linops/test_linops.cc index 44b2c3e8..1fbd3b50 100644 --- a/test/linops/test_linops.cc +++ b/test/linops/test_linops.cc @@ -110,3 +110,93 @@ TEST_F(TestSpectralPrecondLinearOperator, test_diag_n5_k3) { TEST_F(TestSpectralPrecondLinearOperator, test_diag_n5_k4) { run_diag(5, 4, 0.1); } + +// The generic materialize fallback allocates an n by n identity behind the +// caller's back; that hidden allocation is capped. +namespace { +struct IdentityFallbackProbeOp { + using scalar_t = double; + void operator()(blas::Side, Layout, Op, Op, + int64_t, int64_t, int64_t, double, + const double*, int64_t, double, + double*, int64_t) {} +}; +} + +TEST(TestMaterializeGuard, generic_fallback_refuses_huge_identity) { + IdentityFallbackProbeOp A; + int64_t n_over = RandLAPACK::MATERIALIZE_IDENTITY_MAX_DIM + 1; + double stub = 0.0; + // The guard fires before the output buffer is touched, so a stub is safe. + EXPECT_THROW(RandLAPACK::materialize(A, n_over, n_over, &stub, n_over), RandLAPACK::Error); + + // Below the cap the generic path still works. + int64_t n_small = 4; + std::vector buf(n_small * n_small, 1.0); + RandLAPACK::materialize(A, n_small, n_small, buf.data(), n_small); + // The probe writes nothing, so only the zeroing of the buffer happened. + for (auto v : buf) + ASSERT_EQ(v, 0.0); +} + +// The CompositeOperator overload multiplies materialized operands with one +// gemm: no identity, no cap, no call through the composite's operator(). +TEST(TestMaterializeComposite, matches_reference_product) { + int64_t m = 7; + int64_t k = 4; + int64_t n = 5; + int64_t ldb = m + 2; + + // Deterministic operand fills. + vector L_buf(m * k), R_buf(k * n); + for (int64_t i = 0; i < m * k; ++i) + L_buf[i] = 0.5 + (double)((3 * i) % 11); + for (int64_t i = 0; i < k * n; ++i) + R_buf[i] = -1.0 + (double)((5 * i) % 7); + + RandLAPACK::linops::DenseLinOp L_op(m, k, L_buf.data(), m, Layout::ColMajor); + RandLAPACK::linops::DenseLinOp R_op(k, n, R_buf.data(), k, Layout::ColMajor); + RandLAPACK::linops::CompositeOperator comp(m, n, L_op, R_op); + + // Sentinel-fill past-the-column entries to check ldb handling. + vector buf(ldb * n, -42.0); + RandLAPACK::materialize(comp, m, n, buf.data(), ldb); + + vector ref(m * n, 0.0); + blas::gemm(Layout::ColMajor, Op::NoTrans, Op::NoTrans, m, n, k, + 1.0, L_buf.data(), m, R_buf.data(), k, 0.0, ref.data(), m); + + for (int64_t j = 0; j < n; ++j) { + for (int64_t i = 0; i < m; ++i) + ASSERT_DOUBLE_EQ(buf[i + j * ldb], ref[i + j * m]); + for (int64_t i = m; i < ldb; ++i) + ASSERT_EQ(buf[i + j * ldb], -42.0) << "materialize wrote past column height"; + } +} + +TEST(TestMaterializeComposite, wide_composite_is_exempt_from_identity_cap) { + // Wider than the generic fallback allows; the operand-product path handles + // it in m*k + k*n scratch. + int64_t m = 3; + int64_t k = 2; + int64_t n = RandLAPACK::MATERIALIZE_IDENTITY_MAX_DIM + 1; + + vector L_buf(m * k, 1.0), R_buf(k * n); + for (int64_t i = 0; i < k * n; ++i) + R_buf[i] = (double)(i % 5); + + RandLAPACK::linops::DenseLinOp L_op(m, k, L_buf.data(), m, Layout::ColMajor); + RandLAPACK::linops::DenseLinOp R_op(k, n, R_buf.data(), k, Layout::ColMajor); + RandLAPACK::linops::CompositeOperator comp(m, n, L_op, R_op); + + vector buf(m * n, 0.0); + RandLAPACK::materialize(comp, m, n, buf.data(), m); + + // With every row of L equal to ones, each column j of the product holds + // the column sum of R in every entry. + for (int64_t j = 0; j < n; ++j) { + double col_sum = R_buf[j * k] + R_buf[j * k + 1]; + for (int64_t i = 0; i < m; ++i) + ASSERT_DOUBLE_EQ(buf[i + j * m], col_sum); + } +}