Skip to main content

faer/sparse/linalg/
mod.rs

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