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

use crate::algebra::abstr::AbsDiffEq;
use crate::{
    algebra::{
        abstr::{Field, Scalar},
        linear::matrix::{General, QRDec, QRDecomposition, Transpose, UpperTriangular},
    },
    elementary::Power,
};

impl<T> QRDecomposition<T> for General<T>
where
    T: Field + Scalar + Power + AbsDiffEq,
{
    /// QR Decomposition with Givens rotations
    ///
    /// A = QR \
    /// Q is an orthogonal matrix \
    /// R is an upper triangular matrix \
    ///
    /// # Panics
    ///
    /// if A is not a square matrix
    ///
    /// # Example
    ///
    /// ```
    /// use mathru::algebra::linear::matrix::{General, UpperTriangular, QRDecomposition};
    ///
    /// let a: General<f64> = General::new(2, 2, vec![1.0, -2.0, 3.0, -7.0]);
    ///
    /// let (q, r): (General<f64>, UpperTriangular<f64>) = a.dec_qr().unwrap().qr();
    /// ```
    fn dec_qr(&self) -> Result<QRDec<T>, ()> {
        let (m, n) = self.dim();
        debug_assert!(m >= n);

        let mut q: General<T> = General::one(self.m);
        let mut r: General<T> = self.clone();

        for j in 0..self.n {
            for i in (j + 1..self.m).rev() {
                let a_jj: T = r[[j, j]];
                let a_ij: T = r[[i, j]];

                let p: T = (a_jj * a_jj + a_ij * a_ij).sqrt();

                if p.abs_diff_ne(&T::zero(), T::default_epsilon())
                    && a_jj.abs_diff_ne(&T::zero(), T::default_epsilon())
                    && a_ij.abs_diff_ne(&T::zero(), T::default_epsilon())
                {
                    let c: T = a_jj / p;
                    let s: T = -a_ij / p;
                    let g_ij: General<T> = General::givens(r.m, i, j, c, s);

                    r = &g_ij * &r;
                    q = &g_ij * &q;
                }
            }
        }
        q = q.transpose();
        Ok(QRDec::new(q, UpperTriangular::new(r)))
    }
}