opensrdk_linear_algebra/matrix/gt/
trf.rs

1use super::TridiagonalMatrix;
2use crate::matrix::MatrixError;
3use crate::number::*;
4use lapack::dgttrf;
5use lapack::zgttrf;
6use serde::Deserialize;
7use serde::Serialize;
8
9#[derive(Clone, Debug, Serialize, Deserialize)]
10pub struct GTTRF<T = f64>(pub Vec<T>, pub [Vec<T>; 3], pub Vec<i32>)
11where
12    T: Number;
13
14impl TridiagonalMatrix {
15    /// # LU decomposition
16    /// for tridiagonal matrix
17    pub fn gttrf(self) -> Result<GTTRF, MatrixError> {
18        let (mut dl, mut d, mut du) = self.eject();
19        let n = d.len();
20        let mut du2 = vec![0.0; n.max(2) - 2];
21        let mut ipiv = vec![0; n];
22        let mut info = 0;
23
24        let n = n as i32;
25
26        unsafe { dgttrf(n, &mut dl, &mut d, &mut du, &mut du2, &mut ipiv, &mut info) }
27
28        if info != 0 {
29            return Err(MatrixError::LapackRoutineError {
30                routine: "dgttrf".to_owned(),
31                info,
32            });
33        }
34
35        let u = [d, du, du2];
36
37        Ok(GTTRF(dl, u, ipiv))
38    }
39}
40
41impl TridiagonalMatrix<c64> {
42    /// # LU decomposition
43    /// for tridiagonal matrix
44    pub fn gttrf(self) -> Result<GTTRF<c64>, MatrixError> {
45        let (mut dl, mut d, mut du) = self.eject();
46        let n = d.len();
47        let mut du2 = vec![c64::default(); n.max(2) - 2];
48        let mut ipiv = vec![0; n];
49        let mut info = 0;
50
51        let n = n as i32;
52
53        unsafe { zgttrf(n, &mut dl, &mut d, &mut du, &mut du2, &mut ipiv, &mut info) }
54
55        if info != 0 {
56            return Err(MatrixError::LapackRoutineError {
57                routine: "dgttrf".to_owned(),
58                info,
59            });
60        }
61
62        let u = [d, du, du2];
63
64        Ok(GTTRF::<c64>(dl, u, ipiv))
65    }
66}
opensrdk_linear_algebra/matrix/gt/trf.rs

opensrdk_linear_algebra/matrix/gt/
trf.rs
