opensrdk-linear-algebra 0.8.5

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
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113

use super::trf::PPTRF;
use crate::matrix::ge::Matrix;
use crate::matrix::MatrixError;
use crate::number::c64;
use lapack::{dpptrs, zpptrs};

impl PPTRF {
    /// # Solve equation
    ///
    /// with matrix decomposed by potrf
    ///
    /// $$
    /// \mathbf{A} \mathbf{x} = \mathbf{b}
    /// $$
    ///
    /// $$
    /// \mathbf{x} = \mathbf{A}^{-1} \mathbf{b}
    /// $$
    pub fn pptrs(&self, b: Matrix) -> Result<Matrix, MatrixError> {
        let PPTRF(mat) = self;
        let n = mat.dim();

        let mut info = 0;

        let n = n as i32;
        let mut b = b;

        unsafe {
            dpptrs(
                'L' as u8,
                n,
                b.cols() as i32,
                &mat.elems,
                b.elems_mut(),
                n,
                &mut info,
            );
        }

        match info {
            0 => Ok(b),
            _ => Err(MatrixError::LapackRoutineError {
                routine: "dpptrs".to_owned(),
                info,
            }),
        }
    }
}

impl PPTRF<c64> {
    /// # Solve equation
    ///
    /// with matrix decomposed by potrf
    ///
    /// $$
    /// \mathbf{A} \mathbf{x} = \mathbf{b}
    /// $$
    ///
    /// $$
    /// \mathbf{x} = \mathbf{A}^{-1} \mathbf{b}
    /// $$
    pub fn pptrs(&self, b: Matrix<c64>) -> Result<Matrix<c64>, MatrixError> {
        let PPTRF::<c64>(mat) = self;
        let n = mat.dim();

        let mut info = 0;

        let n = n as i32;
        let mut b = b;

        unsafe {
            zpptrs(
                'L' as u8,
                n,
                b.cols() as i32,
                &mat.elems,
                b.elems_mut(),
                n,
                &mut info,
            );
        }

        match info {
            0 => Ok(b),
            _ => Err(MatrixError::LapackRoutineError {
                routine: "zpptrs".to_owned(),
                info,
            }),
        }
    }
}

#[cfg(test)]
mod tests {
    use crate::*;
    #[test]
    fn it_works() {
        let a = vec![2.0, 1.0, 2.0];
        let c = SymmetricPackedMatrix::from(2, a).unwrap();
        let b = mat![
            1.0, 3.0;
            2.0, 4.0
        ];
        let l = c.pptrf().unwrap();
        let x_t = l.pptrs(b).unwrap();

        println!("{:#?}", x_t);
        // assert_eq!(x_t[0][0], 0.0);
        // assert_eq!(x_t[0][1], 1.0);
        // assert_eq!(x_t[1][0], 5.0 / 3.0 - 1.0);
        // assert_eq!(x_t[1][1], 5.0 / 3.0);
    }
}