faer/linalg/cholesky/llt/
factor.rs

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

faer/linalg/cholesky/llt/
factor.rs
