1 use na::Matrix3;
2
3 #[test]
4 #[rustfmt::skip]
lu_simple()5 fn lu_simple() {
6 let m = Matrix3::new(
7 2.0, -1.0, 0.0,
8 -1.0, 2.0, -1.0,
9 0.0, -1.0, 2.0);
10
11 let lu = m.lu();
12 assert_eq!(lu.determinant(), 4.0);
13
14 let (p, l, u) = lu.unpack();
15
16 let mut lu = l * u;
17 p.inv_permute_rows(&mut lu);
18
19 assert!(relative_eq!(m, lu, epsilon = 1.0e-7));
20 }
21
22 #[test]
23 #[rustfmt::skip]
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 = "proptest-support")]
42 mod proptest_tests {
43 macro_rules! gen_tests(
44 ($module: ident, $scalar: expr, $scalar_type: ty) => {
45 mod $module {
46 use na::{DMatrix, Matrix4x3, DVector, Vector4};
47 #[allow(unused_imports)]
48 use crate::core::helper::{RandScalar, RandComplex};
49 use crate::proptest::*;
50 use proptest::{prop_assert, proptest};
51
52 proptest! {
53 #[test]
54 fn lu(m in dmatrix_($scalar)) {
55 let lu = m.clone().lu();
56 let (p, l, u) = lu.unpack();
57 let mut lu = l * u;
58 p.inv_permute_rows(&mut lu);
59
60 prop_assert!(relative_eq!(m, lu, epsilon = 1.0e-7))
61 }
62
63 #[test]
64 fn lu_static_3_5(m in matrix3x5_($scalar)) {
65 let lu = m.lu();
66 let (p, l, u) = lu.unpack();
67 let mut lu = l * u;
68 p.inv_permute_rows(&mut lu);
69
70 prop_assert!(relative_eq!(m, lu, epsilon = 1.0e-7))
71 }
72
73 fn lu_static_5_3(m in matrix5x3_($scalar)) {
74 let lu = m.lu();
75 let (p, l, u) = lu.unpack();
76 let mut lu = l * u;
77 p.inv_permute_rows(&mut lu);
78
79 prop_assert!(relative_eq!(m, lu, epsilon = 1.0e-7));
80 }
81
82 #[test]
83 fn lu_static_square(m in matrix4_($scalar)) {
84 let lu = m.lu();
85 let (p, l, u) = lu.unpack();
86 let mut lu = l * u;
87 p.inv_permute_rows(&mut lu);
88
89 prop_assert!(relative_eq!(m, lu, epsilon = 1.0e-7));
90 }
91
92 #[test]
93 fn lu_solve(n in PROPTEST_MATRIX_DIM, nb in PROPTEST_MATRIX_DIM) {
94 let m = DMatrix::<$scalar_type>::new_random(n, n).map(|e| e.0);
95
96 let lu = m.clone().lu();
97 let b1 = DVector::<$scalar_type>::new_random(n).map(|e| e.0);
98 let b2 = DMatrix::<$scalar_type>::new_random(n, nb).map(|e| e.0);
99
100 let sol1 = lu.solve(&b1);
101 let sol2 = lu.solve(&b2);
102
103 prop_assert!(sol1.is_none() || relative_eq!(&m * sol1.unwrap(), b1, epsilon = 1.0e-6));
104 prop_assert!(sol2.is_none() || relative_eq!(&m * sol2.unwrap(), b2, epsilon = 1.0e-6));
105 }
106
107 #[test]
108 fn lu_solve_static(m in matrix4_($scalar)) {
109 let lu = m.lu();
110 let b1 = Vector4::<$scalar_type>::new_random().map(|e| e.0);
111 let b2 = Matrix4x3::<$scalar_type>::new_random().map(|e| e.0);
112
113 let sol1 = lu.solve(&b1);
114 let sol2 = lu.solve(&b2);
115
116 prop_assert!(sol1.is_none() || relative_eq!(&m * sol1.unwrap(), b1, epsilon = 1.0e-6));
117 prop_assert!(sol2.is_none() || relative_eq!(&m * sol2.unwrap(), b2, epsilon = 1.0e-6));
118 }
119
120 #[test]
121 fn lu_inverse(n in PROPTEST_MATRIX_DIM) {
122 let m = DMatrix::<$scalar_type>::new_random(n, n).map(|e| e.0);
123 let mut l = m.lower_triangle();
124 let mut u = m.upper_triangle();
125
126 // Ensure the matrix is well conditioned for inversion.
127 l.fill_diagonal(na::one());
128 u.fill_diagonal(na::one());
129 let m = l * u;
130
131 let m1 = m.clone().lu().try_inverse().unwrap();
132 let id1 = &m * &m1;
133 let id2 = &m1 * &m;
134
135 prop_assert!(id1.is_identity(1.0e-5));
136 prop_assert!(id2.is_identity(1.0e-5));
137 }
138
139 #[test]
140 fn lu_inverse_static(m in matrix4_($scalar)) {
141 let lu = m.lu();
142
143 if let Some(m1) = lu.try_inverse() {
144 let id1 = &m * &m1;
145 let id2 = &m1 * &m;
146
147 prop_assert!(id1.is_identity(1.0e-5));
148 prop_assert!(id2.is_identity(1.0e-5));
149 }
150 }
151 }
152 }
153 }
154 );
155
156 gen_tests!(complex, complex_f64(), RandComplex<f64>);
157 gen_tests!(f64, PROPTEST_F64, RandScalar<f64>);
158 }
159