1 use na::Matrix3;
2
3 #[test]
4 #[rustfmt::skip]
full_piv_lu_simple()5 fn full_piv_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.full_piv_lu();
12 assert_eq!(lu.determinant(), 4.0);
13
14 let (p, l, u, q) = lu.unpack();
15
16 let mut lu = l * u;
17 p.inv_permute_rows(&mut lu);
18 q.inv_permute_columns(&mut lu);
19
20 assert!(relative_eq!(m, lu, epsilon = 1.0e-7));
21 }
22
23 #[test]
24 #[rustfmt::skip]
full_piv_lu_simple_with_pivot()25 fn full_piv_lu_simple_with_pivot() {
26 let m = Matrix3::new(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.full_piv_lu();
31 assert_eq!(lu.determinant(), -4.0);
32
33 let (p, l, u, q) = lu.unpack();
34
35 let mut lu = l * u;
36 p.inv_permute_rows(&mut lu);
37 q.inv_permute_columns(&mut lu);
38
39 assert!(relative_eq!(m, lu, epsilon = 1.0e-7));
40 }
41
42 #[cfg(feature = "arbitrary")]
43 mod proptest_tests {
44 macro_rules! gen_tests(
45 ($module: ident, $scalar: expr, $scalar_type: ty) => {
46 mod $module {
47 use std::cmp;
48 use num::One;
49 use na::{DMatrix, Matrix4x3, DVector, Vector4};
50 #[allow(unused_imports)]
51 use crate::core::helper::{RandScalar, RandComplex};
52
53 use crate::proptest::*;
54 use proptest::{prop_assert, proptest};
55
56 proptest! {
57 #[test]
58 fn full_piv_lu(m in dmatrix_($scalar)) {
59 let lu = m.clone().full_piv_lu();
60 let (p, l, u, q) = lu.unpack();
61 let mut lu = l * u;
62 p.inv_permute_rows(&mut lu);
63 q.inv_permute_columns(&mut lu);
64
65 prop_assert!(relative_eq!(m, lu, epsilon = 1.0e-7))
66 }
67
68 #[test]
69 fn full_piv_lu_static_3_5(m in matrix3x5_($scalar)) {
70 let lu = m.full_piv_lu();
71 let (p, l, u, q) = lu.unpack();
72 let mut lu = l * u;
73 p.inv_permute_rows(&mut lu);
74 q.inv_permute_columns(&mut lu);
75
76 prop_assert!(relative_eq!(m, lu, epsilon = 1.0e-7))
77 }
78
79 #[test]
80 fn full_piv_lu_static_5_3(m in matrix5x3_($scalar)) {
81 let lu = m.full_piv_lu();
82 let (p, l, u, q) = lu.unpack();
83 let mut lu = l * u;
84 p.inv_permute_rows(&mut lu);
85 q.inv_permute_columns(&mut lu);
86
87 prop_assert!(relative_eq!(m, lu, epsilon = 1.0e-7))
88 }
89
90 #[test]
91 fn full_piv_lu_static_square(m in matrix4_($scalar)) {
92 let lu = m.full_piv_lu();
93 let (p, l, u, q) = lu.unpack();
94 let mut lu = l * u;
95 p.inv_permute_rows(&mut lu);
96 q.inv_permute_columns(&mut lu);
97
98 prop_assert!(relative_eq!(m, lu, epsilon = 1.0e-7))
99 }
100
101 #[test]
102 fn full_piv_lu_solve(n in PROPTEST_MATRIX_DIM, nb in PROPTEST_MATRIX_DIM) {
103 let m = DMatrix::<$scalar_type>::new_random(n, n).map(|e| e.0);
104
105 let lu = m.clone().full_piv_lu();
106 let b1 = DVector::<$scalar_type>::new_random(n).map(|e| e.0);
107 let b2 = DMatrix::<$scalar_type>::new_random(n, nb).map(|e| e.0);
108
109 let sol1 = lu.solve(&b1);
110 let sol2 = lu.solve(&b2);
111
112 prop_assert!(sol1.is_none() || relative_eq!(&m * sol1.unwrap(), b1, epsilon = 1.0e-6));
113 prop_assert!(sol2.is_none() || relative_eq!(&m * sol2.unwrap(), b2, epsilon = 1.0e-6));
114 }
115
116 #[test]
117 fn full_piv_lu_solve_static(m in matrix4_($scalar)) {
118 let lu = m.full_piv_lu();
119 let b1 = Vector4::<$scalar_type>::new_random().map(|e| e.0);
120 let b2 = Matrix4x3::<$scalar_type>::new_random().map(|e| e.0);
121
122 let sol1 = lu.solve(&b1);
123 let sol2 = lu.solve(&b2);
124
125 prop_assert!(sol1.is_none() || relative_eq!(&m * sol1.unwrap(), b1, epsilon = 1.0e-6));
126 prop_assert!(sol2.is_none() || relative_eq!(&m * sol2.unwrap(), b2, epsilon = 1.0e-6));
127 }
128
129 #[test]
130 fn full_piv_lu_inverse(n in PROPTEST_MATRIX_DIM) {
131 let n = cmp::max(1, cmp::min(n, 15)); // To avoid slowing down the test too much.
132 let m = DMatrix::<$scalar_type>::new_random(n, n).map(|e| e.0);
133
134 let mut l = m.lower_triangle();
135 let mut u = m.upper_triangle();
136
137 // Ensure the matrix is well conditioned for inversion.
138 l.fill_diagonal(One::one());
139 u.fill_diagonal(One::one());
140 let m = l * u;
141
142 let m1 = m.clone().full_piv_lu().try_inverse().unwrap();
143 let id1 = &m * &m1;
144 let id2 = &m1 * &m;
145
146 prop_assert!(id1.is_identity(1.0e-5));
147 prop_assert!(id2.is_identity(1.0e-5));
148 }
149
150 #[test]
151 fn full_piv_lu_inverse_static(m in matrix4_($scalar)) {
152 let lu = m.full_piv_lu();
153
154 if let Some(m1) = lu.try_inverse() {
155 let id1 = &m * &m1;
156 let id2 = &m1 * &m;
157
158 prop_assert!(id1.is_identity(1.0e-5));
159 prop_assert!(id2.is_identity(1.0e-5));
160 }
161 }
162 }
163 }
164 }
165 );
166
167 gen_tests!(complex, complex_f64(), RandComplex<f64>);
168 gen_tests!(f64, PROPTEST_F64, RandScalar<f64>);
169 }
170
171 /*
172 #[test]
173 fn swap_rows() {
174 let mut m = Matrix5x3::new(
175 11.0, 12.0, 13.0,
176 21.0, 22.0, 23.0,
177 31.0, 32.0, 33.0,
178 41.0, 42.0, 43.0,
179 51.0, 52.0, 53.0);
180
181 let expected = Matrix5x3::new(
182 11.0, 12.0, 13.0,
183 41.0, 42.0, 43.0,
184 31.0, 32.0, 33.0,
185 21.0, 22.0, 23.0,
186 51.0, 52.0, 53.0);
187
188 m.swap_rows(1, 3);
189
190 assert_eq!(m, expected);
191 }
192
193 #[test]
194 fn swap_columns() {
195 let mut m = Matrix3x5::new(
196 11.0, 12.0, 13.0, 14.0, 15.0,
197 21.0, 22.0, 23.0, 24.0, 25.0,
198 31.0, 32.0, 33.0, 34.0, 35.0);
199
200 let expected = Matrix3x5::new(
201 11.0, 14.0, 13.0, 12.0, 15.0,
202 21.0, 24.0, 23.0, 22.0, 25.0,
203 31.0, 34.0, 33.0, 32.0, 35.0);
204
205 m.swap_columns(1, 3);
206
207 assert_eq!(m, expected);
208 }
209
210 #[test]
211 fn remove_columns() {
212 let m = Matrix3x5::new(
213 11, 12, 13, 14, 15,
214 21, 22, 23, 24, 25,
215 31, 32, 33, 34, 35);
216
217 let expected1 = Matrix3x4::new(
218 12, 13, 14, 15,
219 22, 23, 24, 25,
220 32, 33, 34, 35);
221
222 let expected2 = Matrix3x4::new(
223 11, 12, 13, 14,
224 21, 22, 23, 24,
225 31, 32, 33, 34);
226
227 let expected3 = Matrix3x4::new(
228 11, 12, 14, 15,
229 21, 22, 24, 25,
230 31, 32, 34, 35);
231
232 assert_eq!(m.remove_column(0), expected1);
233 assert_eq!(m.remove_column(4), expected2);
234 assert_eq!(m.remove_column(2), expected3);
235
236 let expected1 = Matrix3::new(
237 13, 14, 15,
238 23, 24, 25,
239 33, 34, 35);
240
241 let expected2 = Matrix3::new(
242 11, 12, 13,
243 21, 22, 23,
244 31, 32, 33);
245
246 let expected3 = Matrix3::new(
247 11, 12, 15,
248 21, 22, 25,
249 31, 32, 35);
250
251 assert_eq!(m.remove_fixed_columns::<U2>(0), expected1);
252 assert_eq!(m.remove_fixed_columns::<U2>(3), expected2);
253 assert_eq!(m.remove_fixed_columns::<U2>(2), expected3);
254
255 // The following is just to verify that the return type dimensions is correctly inferred.
256 let computed: Matrix<_, U3, Dynamic, _> = m.remove_columns(3, 2);
257 assert!(computed.eq(&expected2));
258 }
259
260
261 #[test]
262 fn remove_rows() {
263 let m = Matrix5x3::new(
264 11, 12, 13,
265 21, 22, 23,
266 31, 32, 33,
267 41, 42, 43,
268 51, 52, 53);
269
270 let expected1 = Matrix4x3::new(
271 21, 22, 23,
272 31, 32, 33,
273 41, 42, 43,
274 51, 52, 53);
275
276 let expected2 = Matrix4x3::new(
277 11, 12, 13,
278 21, 22, 23,
279 31, 32, 33,
280 41, 42, 43);
281
282 let expected3 = Matrix4x3::new(
283 11, 12, 13,
284 21, 22, 23,
285 41, 42, 43,
286 51, 52, 53);
287
288 assert_eq!(m.remove_row(0), expected1);
289 assert_eq!(m.remove_row(4), expected2);
290 assert_eq!(m.remove_row(2), expected3);
291
292 let expected1 = Matrix3::new(
293 31, 32, 33,
294 41, 42, 43,
295 51, 52, 53);
296
297 let expected2 = Matrix3::new(
298 11, 12, 13,
299 21, 22, 23,
300 31, 32, 33);
301
302 let expected3 = Matrix3::new(
303 11, 12, 13,
304 21, 22, 23,
305 51, 52, 53);
306
307 assert_eq!(m.remove_fixed_rows::<U2>(0), expected1);
308 assert_eq!(m.remove_fixed_rows::<U2>(3), expected2);
309 assert_eq!(m.remove_fixed_rows::<U2>(2), expected3);
310
311 // The following is just to verify that the return type dimensions is correctly inferred.
312 let computed: Matrix<_, Dynamic, U3, _> = m.remove_rows(3, 2);
313 assert!(computed.eq(&expected2));
314 }
315
316
317 #[test]
318 fn insert_columns() {
319 let m = Matrix5x3::new(
320 11, 12, 13,
321 21, 22, 23,
322 31, 32, 33,
323 41, 42, 43,
324 51, 52, 53);
325
326 let expected1 = Matrix5x4::new(
327 0, 11, 12, 13,
328 0, 21, 22, 23,
329 0, 31, 32, 33,
330 0, 41, 42, 43,
331 0, 51, 52, 53);
332
333 let expected2 = Matrix5x4::new(
334 11, 12, 13, 0,
335 21, 22, 23, 0,
336 31, 32, 33, 0,
337 41, 42, 43, 0,
338 51, 52, 53, 0);
339
340 let expected3 = Matrix5x4::new(
341 11, 12, 0, 13,
342 21, 22, 0, 23,
343 31, 32, 0, 33,
344 41, 42, 0, 43,
345 51, 52, 0, 53);
346
347 assert_eq!(m.insert_column(0, 0), expected1);
348 assert_eq!(m.insert_column(3, 0), expected2);
349 assert_eq!(m.insert_column(2, 0), expected3);
350
351 let expected1 = Matrix5::new(
352 0, 0, 11, 12, 13,
353 0, 0, 21, 22, 23,
354 0, 0, 31, 32, 33,
355 0, 0, 41, 42, 43,
356 0, 0, 51, 52, 53);
357
358 let expected2 = Matrix5::new(
359 11, 12, 13, 0, 0,
360 21, 22, 23, 0, 0,
361 31, 32, 33, 0, 0,
362 41, 42, 43, 0, 0,
363 51, 52, 53, 0, 0);
364
365 let expected3 = Matrix5::new(
366 11, 12, 0, 0, 13,
367 21, 22, 0, 0, 23,
368 31, 32, 0, 0, 33,
369 41, 42, 0, 0, 43,
370 51, 52, 0, 0, 53);
371
372 assert_eq!(m.insert_fixed_columns::<U2>(0, 0), expected1);
373 assert_eq!(m.insert_fixed_columns::<U2>(3, 0), expected2);
374 assert_eq!(m.insert_fixed_columns::<U2>(2, 0), expected3);
375
376 // The following is just to verify that the return type dimensions is correctly inferred.
377 let computed: Matrix<_, U5, Dynamic, _> = m.insert_columns(3, 2, 0);
378 assert!(computed.eq(&expected2));
379 }
380
381
382 #[test]
383 fn insert_rows() {
384 let m = Matrix3x5::new(
385 11, 12, 13, 14, 15,
386 21, 22, 23, 24, 25,
387 31, 32, 33, 34, 35);
388
389 let expected1 = Matrix4x5::new(
390 0, 0, 0, 0, 0,
391 11, 12, 13, 14, 15,
392 21, 22, 23, 24, 25,
393 31, 32, 33, 34, 35);
394
395 let expected2 = Matrix4x5::new(
396 11, 12, 13, 14, 15,
397 21, 22, 23, 24, 25,
398 31, 32, 33, 34, 35,
399 0, 0, 0, 0, 0);
400
401 let expected3 = Matrix4x5::new(
402 11, 12, 13, 14, 15,
403 21, 22, 23, 24, 25,
404 0, 0, 0, 0, 0,
405 31, 32, 33, 34, 35);
406
407 assert_eq!(m.insert_row(0, 0), expected1);
408 assert_eq!(m.insert_row(3, 0), expected2);
409 assert_eq!(m.insert_row(2, 0), expected3);
410
411 let expected1 = Matrix5::new(
412 0, 0, 0, 0, 0,
413 0, 0, 0, 0, 0,
414 11, 12, 13, 14, 15,
415 21, 22, 23, 24, 25,
416 31, 32, 33, 34, 35);
417
418 let expected2 = Matrix5::new(
419 11, 12, 13, 14, 15,
420 21, 22, 23, 24, 25,
421 31, 32, 33, 34, 35,
422 0, 0, 0, 0, 0,
423 0, 0, 0, 0, 0);
424
425 let expected3 = Matrix5::new(
426 11, 12, 13, 14, 15,
427 21, 22, 23, 24, 25,
428 0, 0, 0, 0, 0,
429 0, 0, 0, 0, 0,
430 31, 32, 33, 34, 35);
431
432 assert_eq!(m.insert_fixed_rows::<2>(0, 0), expected1);
433 assert_eq!(m.insert_fixed_rows::<2>(3, 0), expected2);
434 assert_eq!(m.insert_fixed_rows::<2>(2, 0), expected3);
435
436 // The following is just to verify that the return type dimensions is correctly inferred.
437 let computed: Matrix<_, Dynamic, U5, _> = m.insert_rows(3, 2, 0);
438 assert!(computed.eq(&expected2));
439 }
440
441 #[test]
442 fn resize() {
443 let m = Matrix3x5::new(
444 11, 12, 13, 14, 15,
445 21, 22, 23, 24, 25,
446 31, 32, 33, 34, 35);
447
448 let add_add = DMatrix::from_row_slice(5, 6, &[
449 11, 12, 13, 14, 15, 42,
450 21, 22, 23, 24, 25, 42,
451 31, 32, 33, 34, 35, 42,
452 42, 42, 42, 42, 42, 42,
453 42, 42, 42, 42, 42, 42]);
454
455 let del_del = DMatrix::from_row_slice(1, 2, &[11, 12]);
456
457 let add_del = DMatrix::from_row_slice(5, 2, &[
458 11, 12,
459 21, 22,
460 31, 32,
461 42, 42,
462 42, 42]);
463
464 let del_add = DMatrix::from_row_slice(1, 8, &[
465 11, 12, 13, 14, 15, 42, 42, 42]);
466
467 assert_eq!(del_del, m.resize(1, 2, 42));
468 assert_eq!(add_add, m.resize(5, 6, 42));
469 assert_eq!(add_del, m.resize(5, 2, 42));
470 assert_eq!(del_add, m.resize(1, 8, 42));
471 }
472 */
473