rlst/dense/linalg/lapack/
ldlt.rs

1//! Lapack Symmetric/Hermitian Factorization.
2
3use crate::base_types::{RlstResult, UpLo};
4use crate::dense::array::{Array, DynArray};
5use crate::dense::linalg::lapack::interface::sytrf::SytrfUplo;
6use crate::traits::linalg::decompositions::Ldl;
7use crate::traits::linalg::lapack::Lapack;
8use crate::{Shape, UnsafeRandom1DAccessByValue};
9
10impl<Item, ArrayImpl> Ldl for Array<ArrayImpl, 2>
11where
12    Item: Lapack,
13    ArrayImpl: UnsafeRandom1DAccessByValue<Item = Item> + Shape<2>,
14{
15    type Item = Item;
16    fn ldl(&self, uplo: UpLo) -> RlstResult<LdlDecomposition<Item>> {
17        let mut ldl_mat = DynArray::new_from(self);
18        let [m, n] = ldl_mat.shape();
19        assert_eq!(m, n, "Matrix must be square for LDL decomposition.");
20        let n = m;
21        let lda = n;
22        let mut ipiv = vec![0_i32; n];
23
24        let lapack_uplo = match uplo {
25            UpLo::Upper => SytrfUplo::Upper,
26            UpLo::Lower => SytrfUplo::Lower,
27        };
28
29        Item::sytrf(lapack_uplo, n, ldl_mat.data_mut().unwrap(), lda, &mut ipiv)?;
30
31        Ok(LdlDecomposition {
32            ldl: ldl_mat,
33            ipiv,
34            uplo,
35        })
36    }
37}
38
39/// Packed pivoted LDL factors of a square symmetric/Hermitian matrix.
40///
41/// The factors and LAPACK pivot vector are intentionally not exposed; use
42/// [`Self::solve`] to solve systems with the original matrix.
43pub struct LdlDecomposition<Item> {
44    ldl: DynArray<Item, 2>,
45    ipiv: Vec<i32>,
46    uplo: UpLo,
47}
48
49impl<Item> LdlDecomposition<Item>
50where
51    Item: Lapack,
52{
53    /// Solve `A x = b` from the stored factors without modifying `b`.
54    ///
55    /// A vector right-hand side must have length `n`; a matrix right-hand side must have `n`
56    /// rows. The underlying LAPACK wrapper checks the resulting flat storage length, but this
57    /// does not independently reject every incompatible shape when the number of right-hand
58    /// sides is zero. Singular pivots are reported when constructing the factorization.
59    pub fn solve<ArrayImpl, const NDIM: usize>(
60        &self,
61        b: &Array<ArrayImpl, NDIM>,
62    ) -> RlstResult<DynArray<Item, NDIM>>
63    where
64        ArrayImpl: UnsafeRandom1DAccessByValue<Item = Item> + Shape<NDIM>,
65    {
66        let mut sol = DynArray::new_from(b);
67        let n = self.ldl.shape()[0];
68        let nrhs = if NDIM == 1 { 1 } else { sol.shape()[1] };
69
70        let lapack_uplo = match self.uplo {
71            UpLo::Upper => SytrfUplo::Upper,
72            UpLo::Lower => SytrfUplo::Lower,
73        };
74
75        Item::sytrs(
76            lapack_uplo,
77            n,
78            nrhs,
79            self.ldl.data().unwrap(),
80            n,
81            &self.ipiv,
82            sol.data_mut().unwrap(),
83            n,
84        )?;
85
86        Ok(sol)
87    }
88}
89
90#[cfg(test)]
91mod test {
92    use super::*;
93    use crate::MultIntoResize;
94    use crate::dense::array::DynArray;
95    use crate::empty_array;
96    use num::complex::ComplexFloat;
97
98    #[test]
99    fn test_sytrf_f64_lower() {
100        let n = 5;
101        let mut arr = DynArray::<f64, 2>::from_shape([n, n]);
102        for i in 0..n {
103            for j in 0..n {
104                arr[[i, j]] = (i + j) as f64 + 1.0;
105            }
106            arr[[i, i]] += 10.0;
107        }
108
109        let ldl = arr.ldl(UpLo::Lower).unwrap();
110        let mut x_expected = DynArray::<f64, 1>::from_shape([n]);
111        for i in 0..n {
112            x_expected[[i]] = (i + 1) as f64;
113        }
114
115        let rhs = empty_array().simple_mult_into_resize(arr.r(), x_expected.r());
116        let x_actual = ldl.solve(&rhs).unwrap();
117
118        crate::assert_array_relative_eq!(x_actual, x_expected, 1e-12);
119    }
120
121    #[test]
122    fn test_sytrf_f64_upper() {
123        let n = 5;
124        let mut arr = DynArray::<f64, 2>::from_shape([n, n]);
125        for i in 0..n {
126            for j in 0..n {
127                arr[[i, j]] = (i + j) as f64 + 1.0;
128            }
129            arr[[i, i]] += 10.0;
130        }
131
132        let ldl = arr.ldl(UpLo::Upper).unwrap();
133        let mut x_expected = DynArray::<f64, 1>::from_shape([n]);
134        for i in 0..n {
135            x_expected[[i]] = (i + 1) as f64;
136        }
137
138        let rhs = empty_array().simple_mult_into_resize(arr.r(), x_expected.r());
139        let x_actual = ldl.solve(&rhs).unwrap();
140
141        crate::assert_array_relative_eq!(x_actual, x_expected, 1e-12);
142    }
143
144    #[test]
145    fn test_sytrf_c64_lower() {
146        let n = 5;
147        let mut arr = DynArray::<crate::base_types::c64, 2>::from_shape([n, n]);
148        for i in 0..n {
149            for j in 0..n {
150                let val = (i + j) as f64 + 1.0;
151                arr[[i, j]] = crate::base_types::c64::new(val, i as f64 - j as f64);
152            }
153            arr[[i, i]] = crate::base_types::c64::new(arr[[i, i]].re(), 0.0); // Must be Hermitian
154            arr[[i, i]] += 10.0;
155        }
156        // Ensure Hermitian
157        for i in 0..n {
158            for j in i + 1..n {
159                arr[[i, j]] = arr[[j, i]].conj();
160            }
161        }
162
163        let ldl = arr.ldl(UpLo::Lower).unwrap();
164        let mut x_expected = DynArray::<crate::base_types::c64, 1>::from_shape([n]);
165        for i in 0..n {
166            x_expected[[i]] = crate::base_types::c64::new((i + 1) as f64, (i + 1) as f64);
167        }
168
169        let rhs = empty_array().simple_mult_into_resize(arr.r(), x_expected.r());
170        let x_actual = ldl.solve(&rhs).unwrap();
171
172        crate::assert_array_relative_eq!(x_actual, x_expected, 1e-12);
173    }
174
175    #[test]
176    fn test_sytrf_c64_upper() {
177        let n = 5;
178        let mut arr = DynArray::<crate::base_types::c64, 2>::from_shape([n, n]);
179        for i in 0..n {
180            for j in 0..n {
181                let val = (i + j) as f64 + 1.0;
182                arr[[i, j]] = crate::base_types::c64::new(val, i as f64 - j as f64);
183            }
184            arr[[i, i]] = crate::base_types::c64::new(arr[[i, i]].re(), 0.0); // Must be Hermitian
185            arr[[i, i]] += 10.0;
186        }
187        // Ensure Hermitian
188        for i in 0..n {
189            for j in i + 1..n {
190                arr[[i, j]] = arr[[j, i]].conj();
191            }
192        }
193
194        let ldl = arr.ldl(UpLo::Upper).unwrap();
195        let mut x_expected = DynArray::<crate::base_types::c64, 1>::from_shape([n]);
196        for i in 0..n {
197            x_expected[[i]] = crate::base_types::c64::new((i + 1) as f64, (i + 1) as f64);
198        }
199
200        let rhs = empty_array().simple_mult_into_resize(arr.r(), x_expected.r());
201        let x_actual = ldl.solve(&rhs).unwrap();
202
203        crate::assert_array_relative_eq!(x_actual, x_expected, 1e-12);
204    }
205
206    #[test]
207    fn test_sytrf_f32_lower() {
208        let n = 5;
209        let mut arr = DynArray::<f32, 2>::from_shape([n, n]);
210        for i in 0..n {
211            for j in 0..n {
212                arr[[i, j]] = (i + j) as f32 + 1.0;
213            }
214            arr[[i, i]] += 10.0;
215        }
216
217        let ldl = arr.ldl(UpLo::Lower).unwrap();
218        let mut x_expected = DynArray::<f32, 1>::from_shape([n]);
219        for i in 0..n {
220            x_expected[[i]] = (i + 1) as f32;
221        }
222
223        let rhs = empty_array().simple_mult_into_resize(arr.r(), x_expected.r());
224        let x_actual = ldl.solve(&rhs).unwrap();
225
226        crate::assert_array_relative_eq!(x_actual, x_expected, 1e-5);
227    }
228
229    #[test]
230    fn test_sytrf_f32_upper() {
231        let n = 5;
232        let mut arr = DynArray::<f32, 2>::from_shape([n, n]);
233        for i in 0..n {
234            for j in 0..n {
235                arr[[i, j]] = (i + j) as f32 + 1.0;
236            }
237            arr[[i, i]] += 10.0;
238        }
239
240        let ldl = arr.ldl(UpLo::Upper).unwrap();
241        let mut x_expected = DynArray::<f32, 1>::from_shape([n]);
242        for i in 0..n {
243            x_expected[[i]] = (i + 1) as f32;
244        }
245
246        let rhs = empty_array().simple_mult_into_resize(arr.r(), x_expected.r());
247        let x_actual = ldl.solve(&rhs).unwrap();
248
249        crate::assert_array_relative_eq!(x_actual, x_expected, 1e-5);
250    }
251
252    #[test]
253    fn test_sytrf_c32_lower() {
254        let n = 5;
255        let mut arr = DynArray::<crate::base_types::c32, 2>::from_shape([n, n]);
256        for i in 0..n {
257            for j in 0..n {
258                let val = (i + j) as f32 + 1.0;
259                arr[[i, j]] = crate::base_types::c32::new(val, i as f32 - j as f32);
260            }
261            arr[[i, i]] = crate::base_types::c32::new(arr[[i, i]].re(), 0.0); // Must be Hermitian
262            arr[[i, i]] += 10.0;
263        }
264        // Ensure Hermitian
265        for i in 0..n {
266            for j in i + 1..n {
267                arr[[i, j]] = arr[[j, i]].conj();
268            }
269        }
270
271        let ldl = arr.ldl(UpLo::Lower).unwrap();
272        let mut x_expected = DynArray::<crate::base_types::c32, 1>::from_shape([n]);
273        for i in 0..n {
274            x_expected[[i]] = crate::base_types::c32::new((i + 1) as f32, (i + 1) as f32);
275        }
276
277        let rhs = empty_array().simple_mult_into_resize(arr.r(), x_expected.r());
278        let x_actual = ldl.solve(&rhs).unwrap();
279
280        crate::assert_array_relative_eq!(x_actual, x_expected, 1e-5);
281    }
282
283    #[test]
284    fn test_sytrf_c32_upper() {
285        let n = 5;
286        let mut arr = DynArray::<crate::base_types::c32, 2>::from_shape([n, n]);
287        for i in 0..n {
288            for j in 0..n {
289                let val = (i + j) as f32 + 1.0;
290                arr[[i, j]] = crate::base_types::c32::new(val, i as f32 - j as f32);
291            }
292            arr[[i, i]] = crate::base_types::c32::new(arr[[i, i]].re(), 0.0); // Must be Hermitian
293            arr[[i, i]] += 10.0;
294        }
295        // Ensure Hermitian
296        for i in 0..n {
297            for j in i + 1..n {
298                arr[[i, j]] = arr[[j, i]].conj();
299            }
300        }
301
302        let ldl = arr.ldl(UpLo::Upper).unwrap();
303        let mut x_expected = DynArray::<crate::base_types::c32, 1>::from_shape([n]);
304        for i in 0..n {
305            x_expected[[i]] = crate::base_types::c32::new((i + 1) as f32, (i + 1) as f32);
306        }
307
308        let rhs = empty_array().simple_mult_into_resize(arr.r(), x_expected.r());
309        let x_actual = ldl.solve(&rhs).unwrap();
310
311        crate::assert_array_relative_eq!(x_actual, x_expected, 1e-5);
312    }
313}
rlst/dense/linalg/lapack/ldlt.rs

rlst/dense/linalg/lapack/
ldlt.rs
