xingo_sparse/
compressed.rs

1//! # Chapter 5.8 — Sparse matrices and vectors
2//!
3//! *Chapter 1.2 of the first edition, moved here from `xingo-core` with
4//! task L1.2 and made generic in the scalar `T: `[`Real`]
5//! (Chapter 5.6);
6//! `xingo-core` re-exports the `f64` instances at their old paths.*
7//!
8//! ## Purpose
9//! Constraint matrices of real LPs have a handful of nonzeros per column, so
10//! the solver stores only the nonzeros. This chapter defines the two
11//! compressed matrix formats and the work vector used by the simplex kernels:
12//!
13//! * [`CscMatrix`](crate::CscMatrix) — compressed sparse column, the canonical storage of `A`
14//!   (columns are what the simplex method prices and enters; Chapter
15//!   5.8.1);
16//! * [`CsrMatrix`](crate::CsrMatrix) — compressed sparse row, the same matrix stored by rows,
17//!   built on demand for row-wise operations (pivot row `α_r = ρ_r A_N`,
18//!   presolve, propagation; Chapter 5.8.2);
19//! * [`SparseVec`](crate::SparseVec) — a dense work array with an explicit nonzero pattern
20//!   (Chapter 5.8.3);
21//! * [`Triplets`](crate::Triplets) — the coordinate builder: entries pushed in any order,
22//!   duplicates summed, compressed by Algorithm C′ (Chapter 5.8.4; task
23//!   L1.7a);
24//! * [`BlockPartition`](crate::BlockPartition) — the row/column block boundaries of a
25//!   partitioned matrix such as `[A | ±I]` or `[−D Aᵀ; A δI]` (Chapter
26//!   5.8.5; task L1.7a).
27//!
28//! ## Notation
29//! `A ∈ R^{m×n}` with `nnz` stored entries. In CSC, the entries of column `j`
30//! are `row_idx[p], val[p]` for `p ∈ [col_ptr[j], col_ptr[j+1])`. CSR is the
31//! same with the roles of rows and columns exchanged; physically, the CSR
32//! arrays of `A` are exactly the CSC arrays of `Aᵀ`. We call the indexed
33//! dimension (columns for CSC) the *major* one and the other the *minor* one.
34//!
35//! ## Algorithm T (Transposition of a compressed matrix)
36//! Given major pointers `ptr` (length `n_major + 1`), minor indices `idx` and
37//! values `val`, produce the compressed arrays of the transpose.
38//!
39//! T1. [Count.] For each stored entry with minor index `i`, increment
40//!     `cnt[i]`.
41//! T2. [Prefix sums.] Set `tptr[i] = Σ_{k<i} cnt[k]`; then `tptr[n_minor] =
42//!     nnz`. Copy `next ← tptr[0..n_minor]`.
43//! T3. [Scatter.] For each major index `j` in increasing order and each
44//!     entry `(i, v)` of `j` in storage order, write `(j, v)` at position
45//!     `next[i]` and increment `next[i]`.
46//!
47//! Because T3 visits `j` in increasing order, the major indices inside every
48//! output list come out sorted, and entries that tie keep their input order
49//! (the sort is stable). Cost: `O(nnz + n_major + n_minor)` time, no
50//! comparisons.
51//!
52//! ## Invariants of a compressed matrix
53//! (C1) `ptr.len() = n_major + 1`, `ptr[0] = 0`, `ptr` non-decreasing and
54//!      `ptr[n_major] = idx.len() = val.len() = nnz`.
55//! (C2) Inside each major list, minor indices are strictly increasing (sorted,
56//!      no duplicates) and smaller than `n_minor`.
57//! (C3) No explicit zeros: every stored value is finite and nonzero. (The
58//!      contract's (S1) admits zeros; the first edition's stricter rule is
59//!      kept, since the solver relies on it — e.g. the factorization's
60//!      `nnz(B)` statistics and the presolve's counts.)
61//! (C4) `n_major` and `n_minor` fit in `u32`.
62//!
63//! Algorithm T maps a structure satisfying (C1)–(C4) to one satisfying them.
64//!
65//! ## References
66//! Gustavson, "Two fast algorithms for sparse matrices: multiplication and
67//! permuted transposition", ACM TOMS 4(3), 1978; Duff, Erisman & Reid,
68//! *Direct Methods for Sparse Matrices* (1986), §2.13; Davis, *Direct Methods
69//! for Sparse Linear Systems* (SIAM, 2006), §2.4–2.5 (`cs_transpose`,
70//! `cs_compress`, `cs_dupl`).
71
72use crate::SparseError;
73use crate::index::dim_to_u32;
74use crate::real::{Real, is_exact_zero};
75
76/// Compressed arrays `(ptr, idx, val)` of a matrix.
77pub(crate) type Compressed<T> = (Vec<usize>, Vec<u32>, Vec<T>);
78
79/// Algorithm T: transposes compressed arrays with `n_major = ptr.len() − 1`
80/// major lists and `n_minor` minor indices.
81///
82/// Preconditions: (C1), indices `< n_minor`, `n_major ≤ u32::MAX`.
83pub(crate) fn transpose<T: Real>(
84    n_minor: usize,
85    ptr: &[usize],
86    idx: &[u32],
87    val: &[T],
88) -> Compressed<T> {
89    let n_major = ptr.len() - 1;
90    let nnz = ptr[n_major];
91
92    // T1. Count entries per minor index.
93    let mut tptr = vec![0_usize; n_minor + 1];
94    for &i in &idx[..nnz] {
95        tptr[i as usize] += 1;
96    }
97
98    // T2. Exclusive prefix sums; tptr[n_minor] ends at nnz.
99    let mut acc = 0;
100    for c in &mut tptr {
101        let cnt = *c;
102        *c = acc;
103        acc += cnt;
104    }
105    debug_assert_eq!(tptr[n_minor], nnz);
106    let mut next = tptr[..n_minor].to_vec();
107
108    // T3. Scatter by increasing major index (stable).
109    let mut tidx = vec![0_u32; nnz];
110    let mut tval = vec![T::ZERO; nnz];
111    for (j, w) in ptr.windows(2).enumerate() {
112        let jj = u32::try_from(j).expect("major dimension exceeds u32::MAX");
113        for p in w[0]..w[1] {
114            let i = idx[p] as usize;
115            let q = next[i];
116            next[i] += 1;
117            tidx[q] = jj;
118            tval[q] = val[p];
119        }
120    }
121    (tptr, tidx, tval)
122}
123
124/// Checks invariants (C1)–(C4) of compressed arrays. `major`/`minor` name the
125/// dimensions in error messages (e.g. `"column"`, `"row"`).
126pub(crate) fn check_compressed<T: Real>(
127    major: &'static str,
128    minor: &'static str,
129    n_minor: usize,
130    ptr: &[usize],
131    idx: &[u32],
132    val: &[T],
133) -> Result<(), SparseError> {
134    // (C4)
135    dim_to_u32(minor, n_minor)?;
136    let Some((&first, _)) = ptr.split_first() else {
137        return Err(SparseError::Malformed(format!(
138            "{major} pointer array is empty"
139        )));
140    };
141    dim_to_u32(major, ptr.len() - 1)?;
142    // (C1)
143    if first != 0 {
144        return Err(SparseError::Malformed(format!(
145            "{major} pointer starts at {first}, not 0"
146        )));
147    }
148    let nnz = ptr[ptr.len() - 1];
149    if idx.len() != nnz {
150        return Err(SparseError::Dimension {
151            what: "index array",
152            expected: nnz,
153            found: idx.len(),
154        });
155    }
156    if val.len() != nnz {
157        return Err(SparseError::Dimension {
158            what: "value array",
159            expected: nnz,
160            found: val.len(),
161        });
162    }
163    for (j, w) in ptr.windows(2).enumerate() {
164        if w[0] > w[1] {
165            return Err(SparseError::Malformed(format!(
166                "{major} pointer decreases at {major} {j}: {} > {}",
167                w[0], w[1]
168            )));
169        }
170        // (C2)
171        let list = &idx[w[0]..w[1]];
172        if let Some(&last) = list.last()
173            && last as usize >= n_minor
174        {
175            return Err(SparseError::IndexOutOfRange {
176                what: minor,
177                index: last as usize,
178                bound: n_minor,
179            });
180        }
181        if let Some(k) = list.windows(2).position(|p| p[0] >= p[1]) {
182            return Err(SparseError::Malformed(format!(
183                "{minor} indices of {major} {j} not strictly increasing at offset {k}"
184            )));
185        }
186    }
187    // (C3)
188    for (p, &v) in val.iter().enumerate() {
189        if !v.is_finite() || is_exact_zero(v) {
190            return Err(SparseError::InvalidValue {
191                what: "stored value",
192                index: p,
193                value: v.to_f64(),
194            });
195        }
196    }
197    Ok(())
198}
xingo_sparse/compressed.rs

xingo_sparse/
compressed.rs
