Module compressed 

Chapter 5.8 — Sparse matrices and vectors

Chapter 1.2 of the first edition, moved here from xingo-core with task L1.2 and made generic in the scalar T: Real (Chapter 5.6); xingo-core re-exports the f64 instances at their old paths.

Purpose

Constraint matrices of real LPs have a handful of nonzeros per column, so the solver stores only the nonzeros. This chapter defines the two compressed matrix formats and the work vector used by the simplex kernels:

• CscMatrix — compressed sparse column, the canonical storage of A (columns are what the simplex method prices and enters; Chapter 5.8.1);
• CsrMatrix — compressed sparse row, the same matrix stored by rows, built on demand for row-wise operations (pivot row α_r = ρ_r A_N, presolve, propagation; Chapter 5.8.2);
• SparseVec — a dense work array with an explicit nonzero pattern (Chapter 5.8.3);
• Triplets — the coordinate builder: entries pushed in any order, duplicates summed, compressed by Algorithm C′ (Chapter 5.8.4; task L1.7a);
• BlockPartition — the row/column block boundaries of a partitioned matrix such as [A | ±I] or [−D Aᵀ; A δI] (Chapter 5.8.5; task L1.7a).

Notation

A ∈ R^{m×n} with nnz stored entries. In CSC, the entries of column j are row_idx[p], val[p] for p ∈ [col_ptr[j], col_ptr[j+1]). CSR is the same with the roles of rows and columns exchanged; physically, the CSR arrays of A are exactly the CSC arrays of Aᵀ. We call the indexed dimension (columns for CSC) the major one and the other the minor one.

Algorithm T (Transposition of a compressed matrix)

Given major pointers ptr (length n_major + 1), minor indices idx and values val, produce the compressed arrays of the transpose.

T1. [Count.] For each stored entry with minor index i, increment cnt[i]. T2. [Prefix sums.] Set tptr[i] = Σ_{k<i} cnt[k]; then tptr[n_minor] = nnz. Copy next ← tptr[0..n_minor]. T3. [Scatter.] For each major index j in increasing order and each entry (i, v) of j in storage order, write (j, v) at position next[i] and increment next[i].

Because T3 visits j in increasing order, the major indices inside every output list come out sorted, and entries that tie keep their input order (the sort is stable). Cost: O(nnz + n_major + n_minor) time, no comparisons.

Invariants of a compressed matrix

(C1) ptr.len() = n_major + 1, ptr[0] = 0, ptr non-decreasing and ptr[n_major] = idx.len() = val.len() = nnz. (C2) Inside each major list, minor indices are strictly increasing (sorted, no duplicates) and smaller than n_minor. (C3) No explicit zeros: every stored value is finite and nonzero. (The contract’s (S1) admits zeros; the first edition’s stricter rule is kept, since the solver relies on it — e.g. the factorization’s nnz(B) statistics and the presolve’s counts.) (C4) n_major and n_minor fit in u32.

Algorithm T maps a structure satisfying (C1)–(C4) to one satisfying them.

References

Gustavson, “Two fast algorithms for sparse matrices: multiplication and permuted transposition”, ACM TOMS 4(3), 1978; Duff, Erisman & Reid, Direct Methods for Sparse Matrices (1986), §2.13; Davis, Direct Methods for Sparse Linear Systems (SIAM, 2006), §2.4–2.5 (cs_transpose, cs_compress, cs_dupl).