dune_core/linalg/backend/
faer_dense.rs

1//! Dense LU backend over `faer` (partial pivoting). Corpus-sized
2//! circuits are tiny; sparsity (faer sparse or the xingo-sparse
3//! circuit profile) is Phase-3 work with its own benchmark (P-3).
4
5use faer::linalg::solvers::{PartialPivLu, Solve};
6use faer::{Mat, MatMut, c64};
7
8use super::super::{LinalgError, LinearSolver, Pattern};
9
10pub struct FaerDense {
11    n: usize,
12    entries: Vec<(usize, usize)>,
13    /// Diagonal positions present in the pattern (gmin landing spots,
14    /// mirroring `LoadGmin`: only existing diagonals get gmin).
15    diags: Vec<usize>,
16    mat: Mat<f64>,
17    lu: Option<PartialPivLu<f64>>,
18    cmat: Mat<c64>,
19    clu: Option<PartialPivLu<c64>>,
20}
21
22impl Default for FaerDense {
23    fn default() -> Self {
24        FaerDense {
25            n: 0,
26            entries: Vec::new(),
27            diags: Vec::new(),
28            mat: Mat::new(),
29            lu: None,
30            cmat: Mat::new(),
31            clu: None,
32        }
33    }
34}
35
36impl LinearSolver for FaerDense {
37    fn analyze(&mut self, pattern: &Pattern) {
38        self.n = pattern.n;
39        self.entries = pattern.entries.clone();
40        self.diags = pattern
41            .entries
42            .iter()
43            .enumerate()
44            .filter_map(|(k, (r, c))| (r == c).then_some(k))
45            .map(|k| self.entries[k].0)
46            .collect();
47        self.mat = Mat::zeros(self.n, self.n);
48        self.lu = None;
49        self.cmat = Mat::zeros(self.n, self.n);
50        self.clu = None;
51    }
52
53    fn factor(&mut self, values: &[f64], diag_gmin: f64) -> Result<(), LinalgError> {
54        debug_assert_eq!(values.len(), self.entries.len());
55        self.mat.fill(0.0);
56        for ((r, c), v) in self.entries.iter().zip(values) {
57            self.mat[(*r, *c)] += v;
58        }
59        if diag_gmin != 0.0 {
60            for &d in &self.diags {
61                self.mat[(d, d)] += diag_gmin;
62            }
63        }
64        // A structurally empty row or column is the classic "no DC
65        // path to ground" / floating-node singularity; report it
66        // before LU turns it into NaNs.
67        for i in 0..self.n {
68            let row_zero = (0..self.n).all(|j| self.mat[(i, j)] == 0.0);
69            let col_zero = (0..self.n).all(|j| self.mat[(j, i)] == 0.0);
70            if row_zero || col_zero {
71                return Err(LinalgError::Singular { row: i });
72            }
73        }
74        self.lu = Some(self.mat.partial_piv_lu());
75        Ok(())
76    }
77
78    fn solve_in_place(&mut self, rhs: &mut [f64]) -> Result<(), LinalgError> {
79        let lu = self.lu.as_ref().expect("factor before solve");
80        debug_assert_eq!(rhs.len(), self.n);
81        let m = MatMut::from_column_major_slice_mut(rhs, self.n, 1);
82        lu.solve_in_place(m);
83        if let Some(row) = rhs.iter().position(|v| !v.is_finite()) {
84            return Err(LinalgError::Singular { row });
85        }
86        Ok(())
87    }
88
89    fn factor_complex(&mut self, re: &[f64], im: &[f64]) -> Result<(), LinalgError> {
90        debug_assert_eq!(re.len(), self.entries.len());
91        self.cmat.fill(c64::new(0.0, 0.0));
92        for (k, (r, c)) in self.entries.iter().enumerate() {
93            self.cmat[(*r, *c)] += c64::new(re[k], im[k]);
94        }
95        for i in 0..self.n {
96            let row_zero = (0..self.n).all(|j| self.cmat[(i, j)] == c64::new(0.0, 0.0));
97            let col_zero = (0..self.n).all(|j| self.cmat[(j, i)] == c64::new(0.0, 0.0));
98            if row_zero || col_zero {
99                return Err(LinalgError::Singular { row: i });
100            }
101        }
102        self.clu = Some(self.cmat.partial_piv_lu());
103        Ok(())
104    }
105
106    fn solve_complex_in_place(
107        &mut self,
108        rhs_re: &mut [f64],
109        rhs_im: &mut [f64],
110    ) -> Result<(), LinalgError> {
111        let lu = self.clu.as_ref().expect("factor before solve");
112        debug_assert_eq!(rhs_re.len(), self.n);
113        let mut rhs: Vec<c64> = rhs_re
114            .iter()
115            .zip(rhs_im.iter())
116            .map(|(r, i)| c64::new(*r, *i))
117            .collect();
118        let m = MatMut::from_column_major_slice_mut(&mut rhs, self.n, 1);
119        lu.solve_in_place(m);
120        for (k, v) in rhs.iter().enumerate() {
121            if !v.re.is_finite() || !v.im.is_finite() {
122                return Err(LinalgError::Singular { row: k });
123            }
124            rhs_re[k] = v.re;
125            rhs_im[k] = v.im;
126        }
127        Ok(())
128    }
129
130    /// `SMPcaSolve`: solve `Aᵀ x = b` from the matrix stored by the
131    /// last [`LinearSolver::factor_complex`] (a fresh LU of the
132    /// transpose; noise matrices are corpus-sized).
133    fn solve_complex_transposed_in_place(
134        &mut self,
135        rhs_re: &mut [f64],
136        rhs_im: &mut [f64],
137    ) -> Result<(), LinalgError> {
138        debug_assert_eq!(rhs_re.len(), self.n);
139        let at: Mat<c64> = Mat::from_fn(self.n, self.n, |i, j| self.cmat[(j, i)]);
140        let lu = at.partial_piv_lu();
141        let mut rhs: Vec<c64> = rhs_re
142            .iter()
143            .zip(rhs_im.iter())
144            .map(|(r, i)| c64::new(*r, *i))
145            .collect();
146        let m = MatMut::from_column_major_slice_mut(&mut rhs, self.n, 1);
147        lu.solve_in_place(m);
148        for (k, v) in rhs.iter().enumerate() {
149            if !v.re.is_finite() || !v.im.is_finite() {
150                return Err(LinalgError::Singular { row: k });
151            }
152            rhs_re[k] = v.re;
153            rhs_im[k] = v.im;
154        }
155        Ok(())
156    }
157
158    /// Complex determinant: the backend-neutral dense routine of
159    /// [`super::dense_det`] (pole-zero matrices are corpus-sized).
160    fn det_complex(&mut self, re: &[f64], im: &[f64]) -> (num_complex::Complex64, i32) {
161        debug_assert_eq!(re.len(), self.entries.len());
162        super::dense_det::det_complex(self.n, &self.entries, re, im)
163    }
164}
dune_core/linalg/backend/faer_dense.rs

dune_core/linalg/backend/
faer_dense.rs
