xingo_sparse/
csr.rs

1//! # Chapter 5.8.2 — Compressed sparse row storage
2//!
3//! ## Purpose
4//! [`CsrMatrix`] stores `A` by rows. The simplex method computes the pivot
5//! row `α_r = ρ_r A_N` row-wise when `ρ_r` is sparse (Hall & McKinnon 2005),
6//! and presolve and domain propagation walk constraints row by row. The CSR
7//! copy is derived from the canonical [`CscMatrix`] on demand by
8//! [`CscMatrix::to_csr`]; whoever needs it repeatedly (the solver) owns and
9//! caches it, so the solver's `Model` stays a plain value.
10//!
11//! ## Notation
12//! Row `i` holds the pairs `(col_idx[p], val[p])` for
13//! `p ∈ [row_ptr[i], row_ptr[i+1])`. The arrays `(row_ptr, col_idx, val)` of
14//! `A` are exactly the CSC arrays of `Aᵀ`.
15//!
16//! ## Invariants
17//! (C1)–(C4) of Chapter 5.8 with `ptr = row_ptr`, `idx = col_idx`,
18//! `n_major = m`, `n_minor = n`; checked by [`CsrMatrix::check_invariants`].
19//!
20//! ## References
21//! Hall & McKinnon, "Hyper-sparsity in the revised simplex method and how to
22//! exploit it", Comput. Optim. Appl. 32, 2005; Davis (2006), §2.5.
23
24use crate::compressed::{check_compressed, transpose};
25use crate::index::{Col, Row};
26use crate::real::Real;
27use crate::{CscMatrix, SparseError};
28
29/// Sparse matrix in compressed sparse row (CSR) format.
30#[derive(Clone, Debug, PartialEq)]
31pub struct CsrMatrix<T: Real = f64> {
32    /// Number of columns `n`.
33    pub n_cols: usize,
34    /// Row pointers, length `m + 1`; row `i` occupies
35    /// `row_ptr[i]..row_ptr[i + 1]`.
36    pub row_ptr: Vec<usize>,
37    /// Column index of each stored entry, strictly increasing within a row.
38    pub col_idx: Vec<u32>,
39    /// Value of each stored entry; finite and nonzero.
40    pub val: Vec<T>,
41}
42
43impl<T: Real> CsrMatrix<T> {
44    /// Builds a matrix from raw compressed arrays, validating every invariant.
45    pub fn from_raw_parts(
46        n_cols: usize,
47        row_ptr: Vec<usize>,
48        col_idx: Vec<u32>,
49        val: Vec<T>,
50    ) -> Result<Self, SparseError> {
51        let r = Self {
52            n_cols,
53            row_ptr,
54            col_idx,
55            val,
56        };
57        r.check_invariants()?;
58        Ok(r)
59    }
60
61    /// Number of rows `m`.
62    #[inline]
63    pub fn n_rows(&self) -> usize {
64        self.row_ptr.len() - 1
65    }
66
67    /// Number of columns `n`.
68    #[inline]
69    pub fn n_cols(&self) -> usize {
70        self.n_cols
71    }
72
73    /// Number of stored entries.
74    #[inline]
75    pub fn nnz(&self) -> usize {
76        self.val.len()
77    }
78
79    /// Column indices and values of row `i`.
80    #[inline]
81    pub fn row(&self, i: Row) -> (&[u32], &[T]) {
82        let r = self.row_ptr[i.index()]..self.row_ptr[i.index() + 1];
83        (&self.col_idx[r.clone()], &self.val[r])
84    }
85
86    /// Entry `a_ij`, or `0` if not stored; `O(log nnz(row i))`.
87    pub fn get(&self, i: Row, j: Col) -> T {
88        let (cols, vals) = self.row(i);
89        cols.binary_search(&j.0).map_or(T::ZERO, |p| vals[p])
90    }
91
92    /// The same matrix in column-wise (CSC) storage, by Algorithm T.
93    pub fn to_csc(&self) -> CscMatrix<T> {
94        let (col_ptr, row_idx, val) =
95            transpose(self.n_cols, &self.row_ptr, &self.col_idx, &self.val);
96        let a = CscMatrix {
97            n_rows: self.n_rows(),
98            col_ptr,
99            row_idx,
100            val,
101        };
102        debug_assert!(a.check_invariants().is_ok(), "{:?}", a.check_invariants());
103        a
104    }
105
106    /// Checks invariants (C1)–(C4).
107    pub fn check_invariants(&self) -> Result<(), SparseError> {
108        check_compressed(
109            "row",
110            "column",
111            self.n_cols,
112            &self.row_ptr,
113            &self.col_idx,
114            &self.val,
115        )
116    }
117}
118
119#[cfg(test)]
120mod tests {
121    use super::*;
122    use crate::real::{approx_eq, is_exact_zero};
123
124    #[test]
125    fn row_access() {
126        let r = CsrMatrix::from_raw_parts(3, vec![0, 2, 2, 3], vec![0, 2, 1], vec![1.0, 2.0, 3.0])
127            .unwrap();
128        assert_eq!(r.n_rows(), 3);
129        assert_eq!(r.row(Row(0)), (&[0_u32, 2][..], &[1.0, 2.0][..]));
130        assert!(approx_eq(r.get(Row(2), Col(1)), 3.0, 0.0));
131        assert!(is_exact_zero(r.get(Row(1), Col(1))));
132        assert!(CsrMatrix::from_raw_parts(3, vec![0, 1], vec![3], vec![1.0_f64]).is_err());
133    }
134}
xingo_sparse/csr.rs

xingo_sparse/
csr.rs
