gridoxide_sparse/
lu.rs

1//! The sparse LU backend, chosen at run time.
2//!
3//! Every Newton loop in the crate — power flow, state estimation, continuation,
4//! dynamics — factors a real system whose sparsity pattern is fixed and whose
5//! values change every iteration, which is the contract
6//! [`LinearSolver`] states. Each backend implements that trait; [`SparseLu`]
7//! picks one of them from a [`LuBackend`] value, so a caller holds one field
8//! and writes no `match` over the backends, and the `klu`/`pardiso` feature
9//! gates live here and nowhere else.
10//!
11//! The cost is one `match` per factorization, against the factorization
12//! itself.
13
14use crate::klu_native::KluNativeSystem;
15use crate::sparse::{LinearSolver, RealSparseSystem};
16#[cfg(feature = "klu")]
17use crate::sparse_klu::KluRealSystem;
18#[cfg(feature = "pardiso")]
19use crate::sparse_pardiso::PardisoRealSystem;
20
21/// Which sparse LU factors the system.
22#[derive(Clone, Copy, Debug, Default, PartialEq, Eq)]
23pub enum LuBackend {
24    /// `faer`'s sparse LU ([`RealSparseSystem`]). The default.
25    #[default]
26    Faer,
27    /// The pure-Rust port of KLU ([`KluNativeSystem`]). Always compiled.
28    KluNative,
29    /// The vendored SuiteSparse KLU over FFI. Needs the `klu` feature.
30    #[cfg(feature = "klu")]
31    Klu,
32    /// Intel oneMKL's PARDISO over FFI. Needs the `pardiso` feature.
33    #[cfg(feature = "pardiso")]
34    Pardiso,
35}
36
37enum System {
38    Faer(RealSparseSystem),
39    KluNative(KluNativeSystem),
40    #[cfg(feature = "klu")]
41    Klu(KluRealSystem),
42    #[cfg(feature = "pardiso")]
43    Pardiso(PardisoRealSystem),
44}
45
46/// A sparse LU that analyzes a pattern once and refactors it on every solve,
47/// with its backend picked by a [`LuBackend`].
48///
49/// Starts empty. [`analyze`](Self::analyze) builds the symbolic factorization;
50/// the solves reuse it until [`reset`](Self::reset). The same
51/// precondition as [`LinearSolver`] holds: every solve must present the
52/// pattern that was analyzed, in the same order.
53pub struct SparseLu {
54    backend: LuBackend,
55    system: Option<System>,
56    /// What [`adopt_storage_order`](Self::adopt_storage_order) returned, for
57    /// a caller that rebuilds its values' source while keeping this
58    /// factorization.
59    order: Option<Vec<usize>>,
60}
61
62impl SparseLu {
63    pub fn new(backend: LuBackend) -> Self {
64        Self {
65            backend,
66            system: None,
67            order: None,
68        }
69    }
70
71    pub fn backend(&self) -> LuBackend {
72        self.backend
73    }
74
75    /// Whether a symbolic factorization is cached.
76    pub fn is_analyzed(&self) -> bool {
77        self.system.is_some()
78    }
79
80    /// Analyzes `entries`' pattern and factors it, replacing anything cached.
81    /// `false` if the pattern is structurally invalid or the matrix singular,
82    /// in which case nothing is cached.
83    pub fn analyze(&mut self, n: usize, entries: &[(usize, usize, f64)]) -> bool {
84        self.system = match self.backend {
85            LuBackend::Faer => RealSparseSystem::new(n, entries).map(System::Faer),
86            LuBackend::KluNative => KluNativeSystem::new(n, entries).map(System::KluNative),
87            #[cfg(feature = "klu")]
88            LuBackend::Klu => KluRealSystem::new(n, entries).map(System::Klu),
89            #[cfg(feature = "pardiso")]
90            LuBackend::Pardiso => PardisoRealSystem::new(n, entries).map(System::Pardiso),
91        };
92        self.order = None;
93        self.system.is_some()
94    }
95
96    /// Asks the analyzed backend to take values in its own storage order
97    /// from now on, so a caller that can write them that way spares it a
98    /// copy per solve. Returns `order`, where `order[k]` is the index, in the
99    /// entries [`analyze`](Self::analyze) was given, of the value position
100    /// `k` takes; the caller reorders its values to match, once, and every
101    /// later [`factor_and_solve_values`](Self::factor_and_solve_values) and
102    /// [`factor_and_solve`](Self::factor_and_solve) presents them so.
103    ///
104    /// `None`, changing nothing, for a backend without one fixed layout
105    /// (faer, PARDISO), for duplicate positions, or before `analyze`. The
106    /// two KLUs say yes: copying the values into their column-compressed
107    /// layout cost a tenth of a warm Newton solve.
108    pub fn adopt_storage_order(&mut self) -> Option<Vec<usize>> {
109        let order = match self.system.as_mut()? {
110            System::KluNative(s) => s.adopt_storage_order(),
111            #[cfg(feature = "klu")]
112            System::Klu(s) => s.adopt_storage_order(),
113            _ => None,
114        }?;
115        self.order = Some(order.clone());
116        Some(order)
117    }
118
119    /// The order adopted by [`adopt_storage_order`](Self::adopt_storage_order),
120    /// if one is in force.
121    pub fn storage_order(&self) -> Option<&[usize]> {
122        self.order.as_deref()
123    }
124
125    /// Discards the cached factorization, so the next solve needs a fresh
126    /// [`analyze`](Self::analyze).
127    pub fn reset(&mut self) {
128        self.system = None;
129        self.order = None;
130    }
131
132    /// See [`LinearSolver::factor_and_solve`]. `None` if nothing is analyzed.
133    pub fn factor_and_solve(
134        &mut self,
135        entries: &[(usize, usize, f64)],
136        rhs: &[f64],
137    ) -> Option<Vec<f64>> {
138        match self.system.as_mut()? {
139            System::Faer(s) => s.factor_and_solve(entries, rhs),
140            System::KluNative(s) => s.factor_and_solve(entries, rhs),
141            #[cfg(feature = "klu")]
142            System::Klu(s) => s.factor_and_solve(entries, rhs),
143            #[cfg(feature = "pardiso")]
144            System::Pardiso(s) => s.factor_and_solve(entries, rhs),
145        }
146    }
147
148    /// See [`LinearSolver::factor_and_solve_values`]. `None` if nothing is
149    /// analyzed.
150    pub fn factor_and_solve_values(&mut self, values: &[f64], rhs: &[f64]) -> Option<Vec<f64>> {
151        match self.system.as_mut()? {
152            System::Faer(s) => s.factor_and_solve_values(values, rhs),
153            System::KluNative(s) => s.factor_and_solve_values(values, rhs),
154            #[cfg(feature = "klu")]
155            System::Klu(s) => s.factor_and_solve_values(values, rhs),
156            #[cfg(feature = "pardiso")]
157            System::Pardiso(s) => s.factor_and_solve_values(values, rhs),
158        }
159    }
160}
161
162#[cfg(test)]
163mod tests {
164    use super::*;
165
166    /// A 3×3 system given in an order that is not column-major, so the
167    /// adopted order is a real permutation.
168    const ENTRIES: [(usize, usize, f64); 7] = [
169        (0, 0, 4.0),
170        (0, 1, 1.0),
171        (1, 1, 5.0),
172        (2, 1, 2.0),
173        (1, 0, 1.0),
174        (2, 2, 3.0),
175        (1, 2, 1.0),
176    ];
177    const RHS: [f64; 3] = [1.0, 2.0, 3.0];
178
179    fn adopted_order_solves_the_same_system(backend: LuBackend) {
180        let values: Vec<f64> = ENTRIES.iter().map(|e| e.2).collect();
181        let mut lu = SparseLu::new(backend);
182        assert!(lu.analyze(3, &ENTRIES));
183        let before = lu.factor_and_solve_values(&values, &RHS).unwrap();
184
185        let order = lu.adopt_storage_order().expect("KLU adopts its order");
186        assert_ne!(
187            order,
188            (0..ENTRIES.len()).collect::<Vec<_>>(),
189            "not a real permutation"
190        );
191        assert_eq!(lu.storage_order(), Some(&order[..]));
192        let reordered: Vec<f64> = order.iter().map(|&i| values[i]).collect();
193        let after = lu.factor_and_solve_values(&reordered, &RHS).unwrap();
194        assert_eq!(before, after);
195
196        lu.reset();
197        assert_eq!(lu.storage_order(), None);
198    }
199
200    #[test]
201    fn klu_native_adopts_its_storage_order() {
202        adopted_order_solves_the_same_system(LuBackend::KluNative);
203    }
204
205    #[cfg(feature = "klu")]
206    #[test]
207    fn klu_adopts_its_storage_order() {
208        adopted_order_solves_the_same_system(LuBackend::Klu);
209    }
210
211    #[test]
212    fn faer_keeps_the_callers_order() {
213        let mut lu = SparseLu::new(LuBackend::Faer);
214        assert!(lu.analyze(3, &ENTRIES));
215        assert_eq!(lu.adopt_storage_order(), None);
216    }
217
218    #[test]
219    fn duplicate_positions_keep_the_callers_order() {
220        let mut entries = ENTRIES.to_vec();
221        entries.push((0, 0, 1.0));
222        let mut lu = SparseLu::new(LuBackend::KluNative);
223        assert!(lu.analyze(3, &entries));
224        assert_eq!(lu.adopt_storage_order(), None);
225    }
226}
gridoxide_sparse/lu.rs

gridoxide_sparse/
lu.rs
