faer/sparse/linalg/
mod.rs

1use crate::internal_prelude_sp::*;
2
3const CHOLESKY_SUPERNODAL_RATIO_FACTOR: f64 = 40.0;
4const QR_SUPERNODAL_RATIO_FACTOR: f64 = 40.0;
5const LU_SUPERNODAL_RATIO_FACTOR: f64 = 40.0;
6
7/// tuning parameters for the supernodal factorizations
8#[derive(Copy, Clone, Debug)]
9pub struct SymbolicSupernodalParams<'a> {
10	/// supernode relaxation thresholds
11	///
12	/// let `n` be the total number of columns in two adjacent supernodes.
13	/// let `z` be the fraction of zero entries in the two supernodes if they
14	/// are merged (z includes zero entries from prior amalgamations). the
15	/// two supernodes are merged if:
16	///
17	/// `(n <= relax[0].0 && z < relax[0].1) || (n <= relax[1].0 && z < relax[1].1) || ...`
18	pub relax: Option<&'a [(usize, f64)]>,
19}
20
21const DEFAULT_RELAX: &'static [(usize, f64)] = &[(4, 1.0), (16, 0.8), (48, 0.1), (usize::MAX, 0.05)];
22
23impl Default for SymbolicSupernodalParams<'_> {
24	#[inline]
25	fn default() -> Self {
26		Self { relax: Some(DEFAULT_RELAX) }
27	}
28}
29
30/// nonnegative threshold controlling when the supernodal factorization is used
31///
32/// increasing it makes it more likely for the simplicial factorization to be used,
33/// while decreasing it makes it more likely for the supernodal factorization to be used
34///
35/// `1.0` is the default value
36#[derive(Copy, Clone, Debug, PartialEq)]
37pub struct SupernodalThreshold(pub f64);
38
39impl Default for SupernodalThreshold {
40	#[inline]
41	fn default() -> Self {
42		Self(1.0)
43	}
44}
45
46impl SupernodalThreshold {
47	/// determine automatically which variant to select
48	pub const AUTO: Self = Self(1.0);
49	/// simplicial factorization is always selected
50	pub const FORCE_SIMPLICIAL: Self = Self(f64::INFINITY);
51	/// supernodal factorization is always selected
52	pub const FORCE_SUPERNODAL: Self = Self(0.0);
53}
54
55/// sparse $ll^\top$ error
56#[derive(Copy, Clone, Debug)]
57pub enum LltError {
58	/// numerical error
59	Numeric(linalg::cholesky::llt::factor::LltError),
60	/// non algorithmic error
61	Generic(FaerError),
62}
63impl core::fmt::Display for LltError {
64	#[inline]
65	fn fmt(&self, f: &mut core::fmt::Formatter<'_>) -> core::fmt::Result {
66		core::fmt::Debug::fmt(self, f)
67	}
68}
69
70impl core::error::Error for LltError {}
71
72impl From<linalg::cholesky::llt::factor::LltError> for LltError {
73	fn from(value: linalg::cholesky::llt::factor::LltError) -> Self {
74		Self::Numeric(value)
75	}
76}
77/// sparse $lu$ error.
78#[derive(Copy, Clone, Debug)]
79pub enum LuError {
80	/// rank deficient symbolic structure
81	SymbolicSingular {
82		/// iteration at which a pivot could not be found
83		index: usize,
84	},
85	/// non algorithmic error
86	Generic(FaerError),
87}
88
89impl core::fmt::Display for LuError {
90	#[inline]
91	fn fmt(&self, f: &mut core::fmt::Formatter<'_>) -> core::fmt::Result {
92		core::fmt::Debug::fmt(self, f)
93	}
94}
95
96impl core::error::Error for LuError {}
97
98impl<T: Into<FaerError>> From<T> for LltError {
99	fn from(value: T) -> Self {
100		Self::Generic(value.into())
101	}
102}
103impl<T: Into<FaerError>> From<T> for LuError {
104	fn from(value: T) -> Self {
105		Self::Generic(value.into())
106	}
107}
108
109/// sparse matrix multiplication
110pub mod matmul;
111/// sparse matrix triangular solve
112pub mod triangular_solve;
113
114/// high-level sparse matrix solvers
115#[path = "../solvers.rs"]
116pub mod solvers;
117
118pub mod amd;
119pub mod colamd;
120
121pub mod cholesky;
122pub mod lu;
123pub mod qr;
124
125mod ghost {
126	use crate::Index;
127	pub use crate::utils::bound::*;
128
129	pub const NONE_BYTE: u8 = u8::MAX;
130
131	#[inline]
132	pub fn fill_zero<'n, 'a, I: Index>(slice: &'a mut [I], size: Dim<'n>) -> &'a mut [Idx<'n, I>] {
133		let len = slice.len();
134		if len > 0 {
135			assert!(*size > 0);
136		}
137		unsafe {
138			core::ptr::write_bytes(slice.as_mut_ptr(), 0u8, len);
139			&mut *(slice as *mut _ as *mut _)
140		}
141	}
142
143	#[inline]
144	pub fn fill_none<'n, 'a, I: Index>(slice: &'a mut [I::Signed], size: Dim<'n>) -> &'a mut [MaybeIdx<'n, I>] {
145		let _ = size;
146		let len = slice.len();
147		unsafe { core::ptr::write_bytes(slice.as_mut_ptr(), NONE_BYTE, len) };
148		unsafe { &mut *(slice as *mut _ as *mut _) }
149	}
150
151	#[inline]
152	pub fn copy_slice<'n, 'a, I: Index>(dst: &'a mut [I], src: &[Idx<'n, I>]) -> &'a mut [Idx<'n, I>] {
153		let dst: &mut [Idx<'_, I>] = unsafe { &mut *(dst as *mut _ as *mut _) };
154		dst.copy_from_slice(src);
155		dst
156	}
157}