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