gridoxide_sparse/klu_native/
types.rs

1//! Sentinel-index encoding and fixed solver configuration shared across
2//! every `klu_native` submodule — mirrors `vendor/suitesparse/KLU/Include/
3//! klu_internal.h`'s `EMPTY`/`FLIP`/`UNFLIP` (identical, by inspection, to
4//! `AMD/Include/amd_internal.h`'s and `BTF/Include/btf.h`'s own
5//! `BTF_FLIP`/`BTF_UNFLIP`/`BTF_ISFLIPPED`) and `KLU/Source/klu_defaults.c`'s
6//! `KLU_common` defaults.
7//!
8//! Ported deliberately literally, not "idiomatized" into an `Option`/`enum`
9//! at this layer — `lpivot`'s diagonal-preference bookkeeping (`kernel.rs`,
10//! landing in a later phase) reads `Pinv[diagrow] < 0` directly and recovers
11//! `FLIP(Pinv[pivrow])` to reroute it, and BTF's `Q` uses the same encoding
12//! for unmatched columns. This is exactly the kind of index bookkeeping
13//! where `block_sparse.rs`'s `adj.row`/`adj.col` bug happened — silently
14//! wrong only on inputs a casual test wouldn't cover. Conversion to
15//! `Option<usize>`/checked `Vec<usize>` happens only at each phase's own
16//! public API boundary, once the sentinel-form algorithm is done with it.
17
18/// Sentinel for "not yet assigned" — matches every one of the three ported
19/// packages' own `#define EMPTY (-1)` (confirmed identical across
20/// `KLU/Include/klu_internal.h:77`, `AMD/Include/amd_internal.h:71,93`, and
21/// implicitly in `BTF/Include/btf.h`'s `BTF_ISFLIPPED`/`BTF_UNFLIP`, which
22/// treat `-1` the same way).
23pub const EMPTY: i64 = -1;
24
25/// "Negation about -1": marks an otherwise-non-negative index as "flipped"
26/// (BTF uses this for unmatched columns; KLU's `lpivot` uses it to encode a
27/// row's earlier diagonal preference once that row has been claimed by a
28/// different column). `flip(flip(j)) == j` for all `j`, and
29/// `flip(EMPTY) == EMPTY`. Matches `klu_internal.h:78`'s `FLIP` /
30/// `amd_internal.h:72`'s `FLIP` / `btf.h`'s `BTF_FLIP` exactly (all three
31/// upstream definitions are identical: `-(j)-2`).
32pub const fn flip(i: i64) -> i64 {
33    -i - 2
34}
35
36/// `true` iff `i` is currently in "flipped" form (`i < EMPTY`). Matches
37/// `klu_internal.h`'s implicit `UNFLIP` condition / `btf.h`'s
38/// `BTF_ISFLIPPED` exactly.
39pub const fn is_flipped(i: i64) -> bool {
40    i < EMPTY
41}
42
43/// Recovers the original non-negative index whether or not `i` is currently
44/// flipped — a no-op if it wasn't. Matches `klu_internal.h:79`'s `UNFLIP` /
45/// `amd_internal.h:73`'s `UNFLIP` / `btf.h`'s `BTF_UNFLIP` exactly.
46pub const fn unflip(i: i64) -> i64 {
47    if is_flipped(i) { flip(i) } else { i }
48}
49
50/// Fixed solver configuration, mirroring the specific `KLU_common` fields
51/// `sparse_klu.rs::new_common()` (via `klu_defaults`) leaves at their
52/// defaults — confirmed directly from `KLU/Source/klu_defaults.c` — and that
53/// this port's algorithm actually consumes. Kept as a struct (not bare
54/// constants) so differential tests (`#[cfg(all(test, feature = "klu"))]`,
55/// comparing against the vendored C oracle) can vary `tol`/`scale` without
56/// touching any production code path, which always uses `Options::default()`.
57///
58/// Two `KLU_common` fields are deliberately **absent** here, not just
59/// defaulted, because gridoxide's call pattern makes the alternative
60/// branches they'd select genuinely unreachable (confirmed dead, not merely
61/// unused):
62/// - `ordering`: real KLU supports AMD (0), COLAMD (1), user-supplied P/Q
63///   (2), or a user callback (3). `klu_defaults.c` defaults to `0` (AMD) and
64///   `sparse_klu.rs` never overrides it, so this port only ever implements
65///   the AMD path — no `ordering` field exists to accidentally select
66///   another one.
67/// - `maxwork`: bounds the work `btf_maxtrans` may spend before giving up
68///   with a possibly-incomplete matching. Defaults to `0` ("no limit") and
69///   is never overridden, so this port's BTF implementation always runs
70///   `btf_maxtrans` to completion — the work-limit bookkeeping threaded
71///   through the upstream C is dropped outright, not merely defaulted off.
72#[derive(Clone, Copy, Debug, PartialEq)]
73pub struct Options {
74    /// Partial-pivoting threshold: a candidate pivot is accepted only if
75    /// its magnitude is at least `tol` times the largest candidate
76    /// available in that column (otherwise the largest is forced instead).
77    /// `klu_defaults.c`: `Common->tol = 0.001`.
78    pub tol: f64,
79    /// Row-scaling strategy: `1` = sum of absolute values per row, `2` =
80    /// max absolute value per row. `klu_defaults.c`: `Common->scale = 2`.
81    pub scale: i32,
82    /// Whether to BTF-preorder before per-block AMD ordering.
83    /// `klu_defaults.c`: `Common->btf = TRUE`. gridoxide never disables
84    /// this, but it's kept as a field (rather than assumed) so `analyze`
85    /// phase tests can exercise the single-block (`btf: false`) path too.
86    pub btf: bool,
87    /// Abort immediately (and report singular) on a pivot that's still zero
88    /// after the tolerance-based search, rather than substituting a tiny
89    /// nonzero value and continuing. `klu_defaults.c`:
90    /// `Common->halt_if_singular = TRUE`.
91    pub halt_if_singular: bool,
92}
93
94impl Default for Options {
95    fn default() -> Self {
96        // Matches KLU_defaults exactly for every field this port consumes.
97        Self {
98            tol: 0.001,
99            scale: 2,
100            btf: true,
101            halt_if_singular: true,
102        }
103    }
104}
105
106#[cfg(test)]
107mod tests {
108    use super::*;
109
110    #[test]
111    fn flip_is_its_own_inverse() {
112        for j in [-1i64, 0, 1, 5, -5, 1000, -1000] {
113            assert_eq!(flip(flip(j)), j, "flip(flip({j})) should equal {j}");
114        }
115    }
116
117    #[test]
118    fn flip_of_empty_is_empty() {
119        assert_eq!(flip(EMPTY), EMPTY);
120    }
121
122    #[test]
123    fn is_flipped_matches_upstream_condition() {
124        assert!(!is_flipped(EMPTY)); // -1 is not "flipped"
125        assert!(!is_flipped(0));
126        assert!(!is_flipped(5));
127        assert!(is_flipped(flip(0)));
128        assert!(is_flipped(flip(5)));
129    }
130
131    #[test]
132    fn unflip_recovers_original_regardless_of_flip_state() {
133        for j in [0i64, 1, 5, 1000] {
134            assert_eq!(unflip(j), j, "unflip of an already-plain index is a no-op");
135            assert_eq!(
136                unflip(flip(j)),
137                j,
138                "unflip should recover the original index"
139            );
140        }
141    }
142
143    #[test]
144    fn options_default_matches_klu_defaults_c() {
145        let opts = Options::default();
146        assert_eq!(opts.tol, 0.001);
147        assert_eq!(opts.scale, 2);
148        assert!(opts.btf);
149        assert!(opts.halt_if_singular);
150    }
151}
gridoxide_sparse/klu_native/types.rs

gridoxide_sparse/klu_native/
types.rs
