mathru 0.16.2

Fundamental algorithms for scientific computing in Rust
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
114
115
116
117
118
119
120
121
122
123

use crate::algebra::abstr::Zero;
use crate::algebra::linear::matrix::{LUDecomposition, UnitLowerTriangular, UpperTriangular};
use crate::algebra::{
    abstr::{Field, Scalar},
    linear::matrix::{General, LUDec},
};

impl<T> LUDecomposition<T> for General<T>
where
    T: Field + Scalar,
{
    /// Decomposes the matrix into a upper and a lower matrix
    ///
    /// PA = LU
    ///
    /// # Arguments
    ///
    ///
    /// # Example
    ///
    /// ```
    /// use mathru::algebra::linear::matrix::{General, UnitLowerTriangular, UpperTriangular, LUDecomposition};
    /// let a: General<f64> = General::new(2, 2, vec![1.0, -2.0, 3.0, -7.0]);
    ///
    /// let (l, u, p): (UnitLowerTriangular<f64>, UpperTriangular<f64>, General<f64>) = a.dec_lu().unwrap().lup();
    /// ```
    fn dec_lu(&self) -> Result<LUDec<T>, ()> {
        let (m, n): (usize, usize) = self.dim();
        debug_assert_eq!(m, n);

        let (m, n): (usize, usize) = self.dim();
        let m_i32: i32 = m as i32;
        let n_i32: i32 = n as i32;

        let dim_min: i32 = m_i32.min(n_i32);
        let mut ipiv: Vec<i32> = vec![Zero::zero(); dim_min as usize];

        let mut info: i32 = 0;

        let mut self_data = self.clone().data;

        T::xgetrf(
            m_i32,
            n_i32,
            self_data.as_mut_slice(),
            m_i32,
            ipiv.as_mut_slice(),
            &mut info,
        );

        if info != 0 {
            return Err(());
        }

        let mat: General<T> = General::new(m, n, self_data);
        let l: UnitLowerTriangular<T> = General::l(mat.clone());
        let u: UpperTriangular<T> = General::u(mat.clone());
        let p: General<T> = General::p(ipiv);

        Ok(LUDec::new(l, u, p))
    }
}

impl<T> General<T>
where
    T: Field + Scalar,
{
    fn l(mut mat: General<T>) -> UnitLowerTriangular<T> {
        let (m, n): (usize, usize) = mat.dim();

        //fill upper triangle with zero
        for i in 0..m {
            for k in i..n {
                mat[[i, k]] = T::zero();
            }
        }

        //set diagonal to 1
        for i in 0..m {
            mat[[i, i]] = T::one();
        }

        UnitLowerTriangular { matrix: mat }
    }

    fn u(mut mat: General<T>) -> UpperTriangular<T> {
        let (m, _n): (usize, usize) = mat.dim();

        //fill lower triangle with zero
        for i in 0..m {
            for k in 0..i {
                mat[[i, k]] = T::zero();
            }
        }
        UpperTriangular { matrix: mat }
    }

    /// transforms a sequence of permutations to a permutation matrix
    fn p(per: Vec<i32>) -> Self {
        let length = per.len();

        let mut perm: Vec<usize> = vec![0; length];

        for i in 0..length {
            perm[i] = i;
        }

        for i in 0..length - 1 {
            let temp = perm[(per[i] - 1) as usize];
            perm[(per[i] - 1) as usize] = perm[i];
            perm[i] = temp;
        }

        let mut p: General<T> = General::zero(length, length);

        for i in 0..length {
            let k = perm[i];
            p[[i, k]] = T::one();
        }

        p
    }
}