nalgebra/linalg/
svd2.rs

1use crate::{Matrix2, RealField, Vector2, SVD, U2};
2
3// Implementation of the 2D SVD from https://ieeexplore.ieee.org/document/486688
4// See also https://scicomp.stackexchange.com/questions/8899/robust-algorithm-for-2-times-2-svd
5pub fn svd_ordered2<T: RealField>(
6    m: &Matrix2<T>,
7    compute_u: bool,
8    compute_v: bool,
9) -> SVD<T, U2, U2> {
10    let half: T = crate::convert(0.5);
11    let one: T = crate::convert(1.0);
12
13    let e = (m.m11.clone() + m.m22.clone()) * half.clone();
14    let f = (m.m11.clone() - m.m22.clone()) * half.clone();
15    let g = (m.m21.clone() + m.m12.clone()) * half.clone();
16    let h = (m.m21.clone() - m.m12.clone()) * half.clone();
17    let q = (e.clone() * e.clone() + h.clone() * h.clone()).sqrt();
18    let r = (f.clone() * f.clone() + g.clone() * g.clone()).sqrt();
19
20    // Note that the singular values are always sorted because sx >= sy
21    // because q >= 0 and r >= 0.
22    let sx = q.clone() + r.clone();
23    let sy = q - r;
24    let sy_sign = if sy < T::zero() { -one } else { one };
25    let singular_values = Vector2::new(sx, sy * sy_sign.clone());
26
27    if compute_u || compute_v {
28        let a1 = g.atan2(f);
29        let a2 = h.atan2(e);
30        let theta = (a2.clone() - a1.clone()) * half.clone();
31        let phi = (a2 + a1) * half;
32        let (st, ct) = theta.sin_cos();
33        let (sp, cp) = phi.sin_cos();
34
35        let u = Matrix2::new(cp.clone(), -sp.clone(), sp, cp);
36        let v_t = Matrix2::new(ct.clone(), -st.clone(), st * sy_sign.clone(), ct * sy_sign);
37
38        SVD {
39            u: if compute_u { Some(u) } else { None },
40            singular_values,
41            v_t: if compute_v { Some(v_t) } else { None },
42        }
43    } else {
44        SVD {
45            u: None,
46            singular_values,
47            v_t: None,
48        }
49    }
50}
nalgebra/linalg/svd2.rs

nalgebra/linalg/
svd2.rs
