opensrdk-linear-algebra 0.8.8

Standard linear algebra library using blas and lapack for OpenSRDK toolchain.
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65

use crate::matrix::*;
use crate::number::*;
use crate::DiagonalMatrix;
use crate::{bd::BidiagonalMatrix, matrix::st::SymmetricTridiagonalMatrix};
use lapack::dpttrf;
use lapack::zpttrf;
use rayon::prelude::*;
use serde::Deserialize;
use serde::Serialize;

#[derive(Clone, Debug, Serialize, Deserialize)]
pub struct PTTRF<T = f64>(pub BidiagonalMatrix<T>, pub DiagonalMatrix)
where
    T: Number;

impl SymmetricTridiagonalMatrix {
    /// # Cholesky decomposition
    /// for tridiagonal matrix
    /// `T = L * D * L^T`
    pub fn pttrf(self) -> Result<PTTRF, MatrixError> {
        let (mut d, mut e) = self.eject();
        let n = d.len() as i32;
        let mut info = 0;

        unsafe { dpttrf(n, &mut d, &mut e, &mut info) }

        if info != 0 {
            return Err(MatrixError::LapackRoutineError {
                routine: "dpttrf".to_owned(),
                info,
            });
        }

        let bd = BidiagonalMatrix::from(vec![1.0; n as usize], e)?;
        let d = DiagonalMatrix::new(d);

        Ok(PTTRF(bd, d))
    }
}

impl SymmetricTridiagonalMatrix<c64> {
    /// # Cholesky decomposition
    /// for tridiagonal matrix
    /// `T = L * D * L^T`
    pub fn pttrf(self) -> Result<PTTRF<c64>, MatrixError> {
        let (d, mut e) = self.eject();
        let mut d = d.into_par_iter().map(|di| di.re).collect::<Vec<_>>();
        let n = d.len() as i32;
        let mut info = 0;

        unsafe { zpttrf(n, &mut d, &mut e, &mut info) }

        if info != 0 {
            return Err(MatrixError::LapackRoutineError {
                routine: "zpttrf".to_owned(),
                info,
            });
        }

        let bd = BidiagonalMatrix::<c64>::from(vec![c64::one(); n as usize], e)?;
        let d = DiagonalMatrix::new(d);

        Ok(PTTRF::<c64>(bd, d))
    }
}