Skip to main content

faer/linalg/cholesky/llt/
factor.rs

1use crate::internal_prelude::*;
2use crate::linalg::cholesky::ldlt::factor::LdltParams;
3/// dynamic $LL^\top$ regularization.
4///
5/// values below `epsilon` in absolute value, or with the wrong sign are set to
6/// `delta` with their corrected sign
7#[derive(Copy, Clone, Debug)]
8pub struct LltRegularization<T> {
9	/// regularized value
10	pub dynamic_regularization_delta: T,
11	/// regularization threshold
12	pub dynamic_regularization_epsilon: T,
13}
14/// info about the result of the $LL^\top$ factorization
15#[derive(Copy, Clone, Debug)]
16pub struct LltInfo {
17	/// number of pivots whose value or sign had to be corrected
18	pub dynamic_regularization_count: usize,
19}
20/// error in the $LL^\top$ factorization
21#[derive(Copy, Clone, Debug)]
22pub enum LltError {
23	NonPositivePivot { index: usize },
24}
25impl core::fmt::Display for LltError {
26	fn fmt(&self, f: &mut core::fmt::Formatter<'_>) -> core::fmt::Result {
27		core::fmt::Debug::fmt(self, f)
28	}
29}
30impl core::error::Error for LltError {}
31impl<T: RealField> Default for LltRegularization<T> {
32	fn default() -> Self {
33		Self {
34			dynamic_regularization_delta: zero(),
35			dynamic_regularization_epsilon: zero(),
36		}
37	}
38}
39#[derive(Copy, Clone, Debug)]
40pub struct LltParams {
41	pub recursion_threshold: usize,
42	pub block_size: usize,
43	#[doc(hidden)]
44	pub non_exhaustive: NonExhaustive,
45}
46impl<T: ComplexField> Auto<T> for LltParams {
47	#[inline]
48	fn auto() -> Self {
49		let ldlt = <LdltParams as Auto<T>>::auto();
50		Self {
51			recursion_threshold: ldlt.recursion_threshold,
52			block_size: ldlt.block_size,
53			non_exhaustive: NonExhaustive(()),
54		}
55	}
56}
57#[inline]
58pub fn cholesky_in_place_scratch<T: ComplexField>(
59	dim: usize,
60	par: ParShape,
61	params: Spec<LltParams, T>,
62) -> StackReq {
63	_ = par;
64	_ = params;
65	temp_mat_scratch::<T>(dim, 1)
66}
67pub fn cholesky_in_place<T: ComplexField>(
68	A: MatMut<'_, T>,
69	regularization: LltRegularization<T::Real>,
70	par: Par,
71	stack: &mut MemStack,
72	params: Spec<LltParams, T>,
73) -> Result<LltInfo, LltError> {
74	let params = params.config;
75	let N = A.nrows();
76	let mut D = unsafe { temp_mat_uninit(N, 1, stack).0 };
77	let D = D.as_mat_mut();
78	match linalg::cholesky::ldlt::factor::cholesky_block_left_looking(
79		A,
80		D.col_mut(0).transpose_mut(),
81		params.block_size,
82		params.recursion_threshold,
83		params.block_size,
84		true,
85		regularization.dynamic_regularization_delta > zero()
86			&& regularization.dynamic_regularization_epsilon > zero(),
87		&regularization.dynamic_regularization_epsilon,
88		&regularization.dynamic_regularization_delta,
89		None,
90		par,
91	) {
92		Ok(count) => Ok(LltInfo {
93			dynamic_regularization_count: count,
94		}),
95		Err(index) => Err(LltError::NonPositivePivot { index }),
96	}
97}