A14b: quick-wins batch 2 (#168 J zero-fill, #124 materialize guard)

RandLAPACK/drivers/rl_bqrrp.hh

@@ -319,7 +319,8 @@ int BQRRP<T, RNG>::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);

RandLAPACK/drivers/rl_cqrrpt.hh

@@ -12,6 +12,7 @@
 #include <vector>
 #include <chrono>
 #include <numeric>
+#include <algorithm>
 
 using namespace std::chrono;
 
@@ -207,9 +208,13 @@ int CQRRPT<T, RNG>::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<T, RNG> S(DS, state);

RandLAPACK/drivers/rl_cqrrpt_gpu.hh

@@ -15,6 +15,7 @@
 #include <vector>
 #include <chrono>
 #include <numeric>
+#include <algorithm>
 
 using namespace std::chrono;
 
@@ -195,6 +196,10 @@ int CQRRPT_GPU<T, RNG>::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();
     /***********************************************************************************/

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 <typename LinOp>
 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<SpMat>& 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 <typename LinOp1, typename LinOp2>
+void materialize(linops::CompositeOperator<LinOp1, LinOp2>& 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

RandLAPACK/testing/rl_test_utils.hh

@@ -173,20 +173,14 @@ void initialize_test_buffers(::std::vector<T>& C_test, ::std::vector<T>& 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 <typename T, typename LinOp>
 void materialize_linop(LinOp& A_linop, T* A_dense) {
     int64_t m = A_linop.n_rows;

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<double>::epsilon(), 0.85);
+
+    // Generate the rank-deficient input once.
+    auto gen_state = RandBLAS::RNGState();
+    std::vector<double> A_orig(m * n, 0.0);
+    RandLAPACK::gen::mat_gen_info<double> 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<double> 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<double, r123::Philox4x32> 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) {

test/drivers/test_cqrrpt_gpu.cu

@@ -158,3 +158,50 @@ TEST_F(TestCQRRPT, CQRRPT_GPU_full_rank_no_hqrrp) {
     norm_and_copy_computational_helper<double, r123::Philox4x32>(norm_A, all_data);
     test_CQRRPT_general<double, r123::Philox4x32, RandLAPACK::CQRRPT_GPU<double, r123::Philox4x32>>(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<double>::epsilon(), 0.85);
+
+    // Generate the rank-deficient input once.
+    auto gen_state = RandBLAS::RNGState();
+    std::vector<double> A_orig(m * n, 0.0);
+    RandLAPACK::gen::mat_gen_info<double> 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<double, r123::Philox4x32>(m_info, A_orig.data(), gen_state);
+
+    int64_t rank_ref = -1;
+    for (int trial = 0; trial < 3; ++trial) {
+        CQRRPTTestData<double> 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<double, r123::Philox4x32> 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<double, r123::Philox4x32>(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<double, r123::Philox4x32, RandLAPACK::CQRRPT_GPU<double, r123::Philox4x32>>(d_factor, norm_A, all_data, CQRRPT_GPU, state);
+
+        if (trial == 0)
+            rank_ref = all_data.rank;
+        ASSERT_EQ(all_data.rank, rank_ref);
+    }
+}

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<float>(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<double> 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<double> 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<double> L_op(m, k, L_buf.data(), m, Layout::ColMajor);
+    RandLAPACK::linops::DenseLinOp<double> 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<double> buf(ldb * n, -42.0);
+    RandLAPACK::materialize(comp, m, n, buf.data(), ldb);
+
+    vector<double> 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<double> 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<double> L_op(m, k, L_buf.data(), m, Layout::ColMajor);
+    RandLAPACK::linops::DenseLinOp<double> R_op(k, n, R_buf.data(), k, Layout::ColMajor);
+    RandLAPACK::linops::CompositeOperator comp(m, n, L_op, R_op);
+
+    vector<double> 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);
+    }
+}