1 #![cfg_attr(rustfmt, rustfmt_skip)]
2
3 use na::Matrix3;
4
5 #[test]
lu_simple()6 fn lu_simple() {
7 let m = Matrix3::new(
8 2.0, -1.0, 0.0,
9 -1.0, 2.0, -1.0,
10 0.0, -1.0, 2.0);
11
12 let lu = m.lu();
13 assert_eq!(lu.determinant(), 4.0);
14
15 let (p, l, u) = lu.unpack();
16
17 let mut lu = l * u;
18 p.inv_permute_rows(&mut lu);
19
20 assert!(relative_eq!(m, lu, epsilon = 1.0e-7));
21 }
22
23 #[test]
lu_simple_with_pivot()24 fn lu_simple_with_pivot() {
25 let m = Matrix3::new(
26 0.0, -1.0, 2.0,
27 -1.0, 2.0, -1.0,
28 2.0, -1.0, 0.0);
29
30 let lu = m.lu();
31 assert_eq!(lu.determinant(), -4.0);
32
33 let (p, l, u) = lu.unpack();
34
35 let mut lu = l * u;
36 p.inv_permute_rows(&mut lu);
37
38 assert!(relative_eq!(m, lu, epsilon = 1.0e-7));
39 }
40
41 #[cfg(feature = "arbitrary")]
42 mod quickcheck_tests {
43 #[allow(unused_imports)]
44 use crate::core::helper::{RandScalar, RandComplex};
45
46 macro_rules! gen_tests(
47 ($module: ident, $scalar: ty) => {
48 mod $module {
49 use std::cmp;
50 use na::{DMatrix, Matrix4, Matrix4x3, Matrix5x3, Matrix3x5, DVector, Vector4};
51 #[allow(unused_imports)]
52 use crate::core::helper::{RandScalar, RandComplex};
53
54 quickcheck! {
55 fn lu(m: DMatrix<$scalar>) -> bool {
56 let mut m = m;
57 if m.len() == 0 {
58 m = DMatrix::<$scalar>::new_random(1, 1);
59 }
60
61 let m = m.map(|e| e.0);
62
63 let lu = m.clone().lu();
64 let (p, l, u) = lu.unpack();
65 let mut lu = l * u;
66 p.inv_permute_rows(&mut lu);
67
68 relative_eq!(m, lu, epsilon = 1.0e-7)
69 }
70
71 fn lu_static_3_5(m: Matrix3x5<$scalar>) -> bool {
72 let m = m.map(|e| e.0);
73 let lu = m.lu();
74 let (p, l, u) = lu.unpack();
75 let mut lu = l * u;
76 p.inv_permute_rows(&mut lu);
77
78 relative_eq!(m, lu, epsilon = 1.0e-7)
79 }
80
81 fn lu_static_5_3(m: Matrix5x3<$scalar>) -> bool {
82 let m = m.map(|e| e.0);
83 let lu = m.lu();
84 let (p, l, u) = lu.unpack();
85 let mut lu = l * u;
86 p.inv_permute_rows(&mut lu);
87
88 relative_eq!(m, lu, epsilon = 1.0e-7)
89 }
90
91 fn lu_static_square(m: Matrix4<$scalar>) -> bool {
92 let m = m.map(|e| e.0);
93 let lu = m.lu();
94 let (p, l, u) = lu.unpack();
95 let mut lu = l * u;
96 p.inv_permute_rows(&mut lu);
97
98 relative_eq!(m, lu, epsilon = 1.0e-7)
99 }
100
101 fn lu_solve(n: usize, nb: usize) -> bool {
102 if n != 0 && nb != 0 {
103 let n = cmp::min(n, 50); // To avoid slowing down the test too much.
104 let nb = cmp::min(nb, 50); // To avoid slowing down the test too much.
105 let m = DMatrix::<$scalar>::new_random(n, n).map(|e| e.0);
106
107 let lu = m.clone().lu();
108 let b1 = DVector::<$scalar>::new_random(n).map(|e| e.0);
109 let b2 = DMatrix::<$scalar>::new_random(n, nb).map(|e| e.0);
110
111 let sol1 = lu.solve(&b1);
112 let sol2 = lu.solve(&b2);
113
114 return (sol1.is_none() || relative_eq!(&m * sol1.unwrap(), b1, epsilon = 1.0e-6)) &&
115 (sol2.is_none() || relative_eq!(&m * sol2.unwrap(), b2, epsilon = 1.0e-6))
116 }
117
118 return true;
119 }
120
121 fn lu_solve_static(m: Matrix4<$scalar>) -> bool {
122 let m = m.map(|e| e.0);
123 let lu = m.lu();
124 let b1 = Vector4::<$scalar>::new_random().map(|e| e.0);
125 let b2 = Matrix4x3::<$scalar>::new_random().map(|e| e.0);
126
127 let sol1 = lu.solve(&b1);
128 let sol2 = lu.solve(&b2);
129
130 return (sol1.is_none() || relative_eq!(&m * sol1.unwrap(), b1, epsilon = 1.0e-6)) &&
131 (sol2.is_none() || relative_eq!(&m * sol2.unwrap(), b2, epsilon = 1.0e-6))
132 }
133
134 fn lu_inverse(n: usize) -> bool {
135 let n = cmp::max(1, cmp::min(n, 15)); // To avoid slowing down the test too much.
136 let m = DMatrix::<$scalar>::new_random(n, n).map(|e| e.0);
137
138 let mut l = m.lower_triangle();
139 let mut u = m.upper_triangle();
140
141 // Ensure the matrix is well conditioned for inversion.
142 l.fill_diagonal(na::one());
143 u.fill_diagonal(na::one());
144 let m = l * u;
145
146 let m1 = m.clone().lu().try_inverse().unwrap();
147 let id1 = &m * &m1;
148 let id2 = &m1 * &m;
149
150 return id1.is_identity(1.0e-5) && id2.is_identity(1.0e-5);
151 }
152
153 fn lu_inverse_static(m: Matrix4<$scalar>) -> bool {
154 let m = m.map(|e| e.0);
155 let lu = m.lu();
156
157 if let Some(m1) = lu.try_inverse() {
158 let id1 = &m * &m1;
159 let id2 = &m1 * &m;
160
161 id1.is_identity(1.0e-5) && id2.is_identity(1.0e-5)
162 }
163 else {
164 true
165 }
166 }
167 }
168 }
169 }
170 );
171
172 gen_tests!(complex, RandComplex<f64>);
173 gen_tests!(f64, RandScalar<f64>);
174 }
175