relp/algorithm/two_phase/tableau/inverse_maintenance/carry/lower_upper/decomposition/
mod.rs

1/// # LU Decomposition
2///
3
4use std::cmp::Ordering;
5use std::mem;
6use std::ops::{Mul, Neg, Sub};
7
8use num_traits::Zero;
9
10use crate::algorithm::two_phase::tableau::inverse_maintenance::carry::lower_upper::decomposition::pivoting::{Markowitz, PivotRule};
11use crate::algorithm::two_phase::tableau::inverse_maintenance::carry::lower_upper::LUDecomposition;
12use crate::algorithm::two_phase::tableau::inverse_maintenance::carry::lower_upper::permutation::{FullPermutation, Permutation, SwapPermutation};
13use crate::algorithm::two_phase::tableau::inverse_maintenance::ops;
14
15mod pivoting;
16
17impl<F> LUDecomposition<F>
18where
19    F: ops::Field + ops::FieldHR,
20{
21    /// Compute the factorization `PBQ = LU`.
22    ///
23    /// # Arguments
24    ///
25    /// * `rows`: A row major representation of the basis columns.
26    #[must_use]
27    pub fn rows(mut rows: Vec<Vec<(usize, F)>>) -> Self {
28        debug_assert!(rows.iter().all(|row| row.iter().is_sorted_by_key(|&(i, _)| i)));
29
30        // The number of columns gets smaller over time.
31        let m = rows.len();
32        debug_assert!(m > 1);
33
34        // These values will be returned.
35        // They stay the same size, although nothing will be changed in the last indices once they
36        // have been processed.
37        let mut row_permutation = (0..m).collect::<Vec<_>>();
38        let mut column_permutation = (0..m).collect::<Vec<_>>();
39        let mut lower_triangular_row_major = vec![Vec::new(); m - 1];
40
41        // These values get smaller over time.
42        let (mut nnz_row, mut nnz_column) = count_non_zeros(&rows);
43        let pivot_rule = Markowitz::new();
44        for k in 0..m {
45            let (pivot_row, pivot_column) = pivot_rule.choose_pivot(&nnz_row, &nnz_column, &rows, k);
46            // Administration for swapping (pivot_row, pivot_column) to (k, k)
47            swap(
48                pivot_row, pivot_column,
49                k,
50                &mut row_permutation, &mut column_permutation,
51                &mut nnz_row, &mut nnz_column,
52                &mut rows,
53                &mut lower_triangular_row_major,
54            );
55
56            // Update the nonzero counters for this row removal
57            for &(j, _) in &rows[k] {
58                nnz_row[k] -= 1;
59                nnz_column[j] -= 1;
60            }
61
62            // Rest of the loop: row reduction operations to zero the values below the pivot
63            let (current_row, remaining_rows) = {
64                let (current_row, remaining_rows) = rows[k..].split_first_mut().unwrap();
65                (&*current_row, remaining_rows)
66            };
67            let pivot_value = &current_row[0].1;
68
69            // Compute the factor with which row k should be multiplied and subtracted from the
70            // other rows
71            let mut ratios_to_subtract = Vec::with_capacity(nnz_column[k]);
72            for (i, row) in remaining_rows.iter_mut().enumerate() {
73                debug_assert!(!row.is_empty(), "The first item exists (invertibility).");
74                if row[0].0 == k {
75                    ratios_to_subtract.push((i, row.remove(0).1 / pivot_value));
76                    nnz_row[k + 1 + i] -= 1;
77                    nnz_column[k] -= 1;
78                }
79            }
80
81            // Eliminate the zeros below the pivot
82            for (i, ratio) in ratios_to_subtract {
83                let (nnz_row_net_difference, nnz_column_removed, nnz_column_added) = subtract_multiple_of_row_from_other_row(
84                    &mut remaining_rows[i],
85                    &ratio,
86                    &current_row[1..],
87                );
88
89                nnz_row[k + 1 + i] -= nnz_row_net_difference.0;
90                nnz_row[k + 1 + i] += nnz_row_net_difference.1;
91                for column in nnz_column_removed {
92                    nnz_column[column] -= 1;
93                }
94                for column in nnz_column_added {
95                    nnz_column[column] += 1;
96                }
97
98                lower_triangular_row_major[k + 1 + i - 1].push((k, ratio));
99            }
100
101            debug_assert_eq!(nnz_row[k], 0);
102            debug_assert_eq!(nnz_column[k], 0);
103            debug_assert!(nnz_row[(k + 1)..].iter()
104                .enumerate()
105                .all(|(i, count)| rows[k + 1 + i].len() == *count));
106        }
107
108        // We collected row-major but need column major, so transpose
109        let mut upper_triangular = vec![Vec::new(); m - 1];
110        let mut upper_diagonal = Vec::with_capacity(m);
111        for (i, row) in rows.into_iter().enumerate() {
112            let mut iter = row.into_iter();
113            let (index, diagonal) = iter.next().unwrap();
114            debug_assert_eq!(index, i);
115            upper_diagonal.push(diagonal);
116            for (j, value) in iter {
117                upper_triangular[j - 1].push((i, value));
118            }
119        }
120        let mut lower_triangular = vec![Vec::new(); m - 1];
121        for (i, row) in lower_triangular_row_major.into_iter().enumerate()
122            // We collected starting at row 1, but the lowest row index was 1.
123            .map(|(i, v)| (i + 1, v)) {
124            for (j, v) in row {
125                lower_triangular[j].push((i, v));
126            }
127        }
128
129        // Compute the inverse permutations and invert
130        let mut row_permutation = FullPermutation::new(row_permutation);
131        row_permutation.invert();
132        let mut column_permutation = FullPermutation::new(column_permutation);
133        column_permutation.invert();
134
135        Self {
136            row_permutation,
137            column_permutation,
138            lower_triangular,
139            upper_triangular,
140            upper_diagonal,
141            updates: Vec::new(),
142        }
143    }
144}
145
146fn subtract_multiple_of_row_from_other_row<T>(
147    to_edit: &mut Vec<(usize, T)>,
148    ratio: &T,
149    being_removed: &[(usize, T)],
150) -> ((usize, usize), Vec<usize>, Vec<usize>)
151where
152    for<'r> T: Mul<&'r T, Output=T> + Sub<T, Output=T> + Zero + Eq,
153    for<'r> &'r T: Neg<Output=T> + Mul<&'r T, Output=T>,
154{
155    debug_assert!(!ratio.is_zero());
156
157    let mut column_nnz_added = Vec::new();
158    if to_edit.is_empty() {
159        // TODO(PERFORMANCE): Is this case even possible?
160        let row_nnz_added = being_removed.len();
161        for (j, v) in being_removed {
162            to_edit.push((*j, -ratio * v));
163            column_nnz_added.push(*j);
164        }
165        ((0, row_nnz_added), Vec::with_capacity(0), column_nnz_added)
166    } else {
167        let old_row = mem::replace(to_edit, Vec::with_capacity(0));
168        let old_row_nnz = old_row.len();
169
170        let mut column_nnz_removed = Vec::new();
171        let mut new = Vec::new();
172
173        let mut index = 0;
174        for (j, old_value) in old_row {
175            while index < being_removed.len() && being_removed[index].0 < j {
176                new.push((being_removed[index].0, -ratio * &being_removed[index].1));
177                column_nnz_added.push(being_removed[index].0);
178                index += 1;
179            }
180
181            if index < being_removed.len() && being_removed[index].0 == j {
182                let product = ratio * &being_removed[index].1;
183                if product != old_value {
184                    new.push((j, old_value - product));
185                } else {
186                    column_nnz_removed.push(j);
187                }
188                index += 1;
189            } else {
190                new.push((j, old_value));
191            }
192        }
193        while index < being_removed.len() {
194            new.push((being_removed[index].0, -ratio * &being_removed[index].1));
195            column_nnz_added.push(being_removed[index].0);
196            index += 1;
197        }
198
199        let new_row_nnz = new.len();
200        *to_edit = new;
201
202        let (row_nnz_net_removed, row_nnz_net_added) = match new_row_nnz.cmp(&old_row_nnz) {
203            Ordering::Less => (old_row_nnz - new_row_nnz, 0),
204            Ordering::Equal => (0, 0),
205            Ordering::Greater => (0, new_row_nnz - old_row_nnz),
206        };
207
208        ((row_nnz_net_removed, row_nnz_net_added), column_nnz_removed, column_nnz_added)
209    }
210}
211
212/// Swap `pivot_row` and `pivot_column` to k.
213///
214/// # Arguments
215///
216/// * `pivot_row`: Index of row swapped with index k.
217/// * `pivot_column`: Index of column swapped with index k.
218/// * `k`: Step in the iteration and row and column index that will be swapped with.
219/// * `row_permutation`: Maps from current row index to the index the row had in the original
220/// matrix.
221/// * `row_permutation`: Maps from current column index to the index the column had in the original
222/// matrix.
223/// * `nnz_row`: Number of nonzero elements per row.
224fn swap<T>(
225    pivot_row: usize, pivot_column: usize,
226    k: usize,
227    row_permutation: &mut [usize], column_permutation: &mut [usize],
228    nnz_row: &mut [usize], nnz_column: &mut [usize],
229    rows: &mut [Vec<(usize, T)>],
230    lower_triangular_row_major: &mut [Vec<(usize, T)>],
231) {
232    let n = rows.len();
233    debug_assert!(pivot_row < n);
234    debug_assert!(pivot_column < n);
235    debug_assert!(k < n);
236    debug_assert_eq!(row_permutation.len(), n);
237    debug_assert_eq!(column_permutation.len(), n);
238    debug_assert_eq!(nnz_row.len(), n);
239    debug_assert_eq!(nnz_column.len(), n);
240    debug_assert_eq!(rows.len(), n);
241    debug_assert_eq!(lower_triangular_row_major.len(), n - 1);
242
243    // We work from low indices to high, and don't change the lower (< k) indices.
244    debug_assert!(pivot_row >= k && pivot_column >= k);
245
246    if pivot_row != k {
247        row_permutation.swap(pivot_row, k);
248        nnz_row.swap(pivot_row, k);
249        rows.swap(pivot_row, k);
250
251        if pivot_row == 0 && k > 0 {
252            debug_assert!(lower_triangular_row_major[k - 1].is_empty());
253        }
254        if pivot_row > 0 && k == 0 {
255            debug_assert!(lower_triangular_row_major[pivot_row - 1].is_empty());
256        }
257        if pivot_row > 0 && k > 0 {
258            // If the pivot_row > 0 && k == 0, there is nothing to swap anyway
259            lower_triangular_row_major.swap(pivot_row - 1, k - 1);
260        }
261    }
262
263    if pivot_column != k {
264        column_permutation.swap(pivot_column, k);
265        nnz_column.swap(pivot_column, k);
266
267        // TODO(ENHANCEMENT): Can some swapping by avoided by doing it all at the end?
268        let swap = SwapPermutation::new((pivot_column, k), n);
269        for row in rows {
270            swap.forward_sorted(row);
271        }
272    }
273}
274
275/// Count the number of nonzero values in each row and column.
276///
277/// Requires iterating over all values once.
278fn count_non_zeros<T>(rows: &[Vec<(usize, T)>]) -> (Vec<usize>, Vec<usize>) {
279    let n = rows.len();
280    debug_assert!(n > 0);
281    debug_assert!(rows.iter().all(|row| row.iter().all(|&(i, _)| i < n)));
282
283    let nnz_row = rows.iter()
284        .map(Vec::len)
285        .collect::<Vec<_>>();
286
287    let nnz_column = {
288        let mut counts = vec![0; n];
289        for column in rows {
290            for &(i, _) in column {
291                counts[i] += 1;
292            }
293        }
294        counts
295    };
296
297    debug_assert_eq!(nnz_column.len(), n);
298    debug_assert_eq!(nnz_row.len(), n);
299    // We should be working with an invertible matrix
300    debug_assert!(nnz_column.iter().all(|&count| count > 0));
301    debug_assert!(nnz_row.iter().all(|&count| count > 0));
302
303    (nnz_row, nnz_column)
304}
305
306#[cfg(test)]
307mod test {
308    use relp_num::{RB, R8, RationalBig, NonZero, Rational8};
309
310    use crate::algorithm::two_phase::matrix_provider::column::{Column, SparseSliceIterator, DenseSliceIterator};
311    use crate::algorithm::two_phase::matrix_provider::column::identity::IdentityColumn;
312    use crate::algorithm::two_phase::tableau::inverse_maintenance::carry::BasisInverse;
313    use crate::algorithm::two_phase::tableau::inverse_maintenance::carry::lower_upper::LUDecomposition;
314    use crate::algorithm::two_phase::tableau::inverse_maintenance::carry::lower_upper::permutation::FullPermutation;
315    use crate::algorithm::two_phase::tableau::inverse_maintenance::ColumnComputationInfo;
316    use crate::data::linear_algebra::vector::{SparseVector, Vector};
317    use std::collections::VecDeque;
318
319    #[test]
320    fn identity_2() {
321        let rows = vec![vec![(0, RB!(1))], vec![(1, RB!(1))]];
322        let result = LUDecomposition::rows(rows);
323        let expected = LUDecomposition {
324            row_permutation: FullPermutation::identity(2),
325            column_permutation: FullPermutation::identity(2),
326            lower_triangular: vec![vec![]],
327            upper_triangular: vec![vec![]],
328            upper_diagonal: vec![RB!(1), RB!(1)],
329            updates: vec![],
330        };
331
332        assert_eq!(result, expected);
333    }
334
335    #[test]
336    fn identity_3() {
337        let rows = vec![vec![(0, RB!(1))], vec![(1, RB!(1))], vec![(2, RB!(1))]];
338        let result = LUDecomposition::rows(rows);
339        let expected = LUDecomposition {
340            row_permutation: FullPermutation::identity(3),
341            column_permutation: FullPermutation::identity(3),
342            lower_triangular: vec![vec![], vec![]],
343            upper_triangular: vec![vec![], vec![]],
344            upper_diagonal: vec![RB!(1), RB!(1), RB!(1)],
345            updates: vec![],
346        };
347
348        assert_eq!(result, expected);
349    }
350
351    #[test]
352    fn offdiagonal_2_upper() {
353        let rows = vec![vec![(0, RB!(1)), (1, RB!(1))], vec![(1, RB!(1))]];
354        let result = LUDecomposition::rows(rows);
355        let expected = LUDecomposition {
356            row_permutation: FullPermutation::identity(2),
357            column_permutation: FullPermutation::identity(2),
358            lower_triangular: vec![vec![]],
359            upper_triangular: vec![vec![(0, RB!(1))]],
360            upper_diagonal: vec![RB!(1), RB!(1)],
361            updates: vec![],
362        };
363
364        assert_eq!(result, expected);
365    }
366
367    #[test]
368    fn offdiagonal_2_lower() {
369        let rows = vec![vec![(0, RB!(1))], vec![(0, RB!(1)), (1, RB!(1))]];
370        let result = LUDecomposition::rows(rows);
371        let expected = LUDecomposition {
372            row_permutation: FullPermutation::identity(2),
373            column_permutation: FullPermutation::identity(2),
374            lower_triangular: vec![vec![(1, RB!(1))]],
375            upper_triangular: vec![vec![]],
376            upper_diagonal: vec![RB!(1), RB!(1)],
377            updates: vec![],
378        };
379
380        assert_eq!(result, expected);
381    }
382
383    #[test]
384    fn offdiagonal_2_both() {
385        let rows = vec![vec![(0, RB!(1)), (1, RB!(1))], vec![(0, RB!(1))]];
386        let result = LUDecomposition::rows(rows);
387        let expected = LUDecomposition {
388            row_permutation: FullPermutation::new(vec![1, 0]),
389            column_permutation: FullPermutation::identity(2),
390            lower_triangular: vec![vec![(1, RB!(1))]],
391            upper_triangular: vec![vec![]],
392            upper_diagonal: vec![RB!(1), RB!(1)],
393            updates: vec![],
394        };
395
396        assert_eq!(result, expected);
397    }
398
399    #[test]
400    fn wikipedia_example() {
401        let columns = vec![vec![(0, RB!(4)), (1, RB!(3))], vec![(0, RB!(6)), (1, RB!(3))]];
402        let result = LUDecomposition::rows(columns);
403        let expected = LUDecomposition {
404            row_permutation: FullPermutation::identity(2),
405            column_permutation: FullPermutation::identity(2),
406            lower_triangular: vec![vec![(1, RB!(3, 2))]],
407            upper_triangular: vec![vec![(0, RB!(3))]],
408            upper_diagonal: vec![RB!(4), RB!(-3, 2)],
409            updates: vec![],
410        };
411
412        assert_eq!(result, expected);
413    }
414
415    #[test]
416    fn wikipedia_example2() {
417        let rows = vec![vec![(0, RB!(-1)), (1, RB!(3, 2))], vec![(0, RB!(1)), (1, RB!(-1))]];
418        let result = LUDecomposition::rows(rows);
419        let expected = LUDecomposition {
420            row_permutation: FullPermutation::identity(2),
421            column_permutation: FullPermutation::identity(2),
422            lower_triangular: vec![vec![(1, RB!(-1))]],
423            upper_triangular: vec![vec![(0, RB!(3, 2))]],
424            upper_diagonal: vec![RB!(-1), RB!(1, 2)],
425            updates: vec![],
426        };
427
428        assert_eq!(result, expected);
429
430        assert_eq!(
431            expected.left_multiply_by_basis_inverse(IdentityColumn::new(0).iter()).into_column(),
432            SparseVector::new(vec![(0, RB!(2)), (1, RB!(2))], 2),
433        );
434        assert_eq!(
435            expected.left_multiply_by_basis_inverse(IdentityColumn::new(1).iter()).into_column(),
436            SparseVector::new(vec![(0, RB!(3)), (1, RB!(2))], 2),
437        );
438    }
439
440    pub fn to_columns<const M: usize>(rows: &[[i32; M]; M]) -> VecDeque<Vec<(usize, Rational8)>> {
441        let mut columns = vec![vec![]; M].into_iter().collect::<VecDeque<_>>();
442
443        for (i, row) in rows.into_iter().enumerate() {
444            for (j, v) in row.into_iter().enumerate() {
445                if v.is_not_zero() {
446                    columns[j].push((i, R8!(*v)));
447                }
448            }
449        }
450
451        columns
452    }
453
454    fn test_matrix<const M: usize>(rows: [[i32; M]; M]) {
455        let columns = to_columns(&rows);
456        let result = LUDecomposition::<RationalBig>::rows(
457            rows.iter().map(|row| {
458                row.iter().enumerate()
459                    .filter(|(_, v)| v.is_not_zero())
460                    .map(|(i, v)| (i, v.into()))
461                    .collect()
462            }).collect()
463        );
464        for (j, column) in columns.iter().enumerate() {
465            assert_eq!(
466                result.left_multiply_by_basis_inverse(SparseSliceIterator::new(column)).into_column(),
467                SparseVector::standard_basis_vector(j, M),
468                "{}", j,
469            );
470        }
471        for (j, row) in rows.into_iter().enumerate() {
472            assert_eq!(
473                result.right_multiply_by_basis_inverse(DenseSliceIterator::new(&row)),
474                SparseVector::standard_basis_vector(j, M),
475                "{}", j,
476            );
477        }
478    }
479
480    #[test]
481    fn test_3x3() {
482        test_matrix([
483            [ 2,  3,  0],
484            [ 5,  0, 11],
485            [23, 29,  0],
486        ]);
487    }
488
489    #[test]
490    fn test_4x4_1() {
491        test_matrix([
492            [ 2,  3,  0,  5],
493            [ 5,  0, 11, 13],
494            [23, 29,  0, 57],
495            [31, 37, 41,  0],
496        ]);
497    }
498
499    #[test]
500    fn test_4x4_2() {
501        test_matrix([
502            [-101,    0,    0,   -5],
503            [-110,  -81,    0,    0],
504            [   0,    0,    1, -111],
505            [   0,   93,   69,    0],
506        ]);
507    }
508
509    #[test]
510    fn test_4x4_3() {
511        test_matrix([
512            [  0,   0, -84, 122],
513            [  0,   9,   0,   0],
514            [-39, 115,   0,  57],
515            [  0, -12, 121,   0],
516        ]);
517    }
518
519    #[test]
520    fn test_5x5_banded() {
521        test_matrix([
522            [2,  3,  0,  0,  0],
523            [5,  7, 11,  0,  0],
524            [0, 29, 13, 57,  0],
525            [0,  0, 41, 17,  0],
526            [0,  0,  0, 53, 51],
527        ]);
528    }
529
530    #[test]
531    fn test_5x5_1() {
532        test_matrix([
533            [29, 23,  0, 19, 0],
534            [ 0,  0, 17, 13, 0],
535            [ 0,  0,  7,  0, 0],
536            [ 5,  0,  0,  3, 0],
537            [ 0,  0,  0,  0, 2],
538        ]);
539    }
540
541    #[test]
542    fn test_5x5_2() {
543        test_matrix([
544            [29, 23,  0, 19, 0],
545            [ 0,  0, 17, 13, 0],
546            [ 0, 11,  7,  0, 0],
547            [ 5,  0,  0,  3, 0],
548            [ 0,  0,  0,  0, 2],
549        ]);
550    }
551
552    #[test]
553    fn test_5x5_3() {
554        test_matrix([
555            [ 2,  3,  0,  5,  7],
556            [ 5,  0, 11, 13, 17],
557            [23, 29,  0, 57, 59],
558            [31, 37, 41,  0,  0],
559            [43,  0, 47, 53, 51],
560        ]);
561    }
562
563    #[test]
564    fn test_5x5_4() {
565        test_matrix([
566            [ 2,  3,  0,  5,  7],
567            [ 5,  0, 11, 13, 17],
568            [23, 29,  0, 57, 59],
569            [31, 37, 41,  0,  0],
570            [43,  0, 47, 53, 51],
571        ]);
572    }
573
574    #[test]
575    fn test_5x5_5() {
576        test_matrix([
577            [   0,   54,   43,    0,   84],
578            [   4,    0,    0,    0,    0],
579            [   0, -111,  -27,    0,  -86],
580            [  -6,    0,    0,   17,  -62],
581            [-109,    0,    0,    0, -104],
582        ]);
583    }
584
585    #[test]
586    fn test_5x5_6() {
587        test_matrix([
588            [ -71, -124,    0,    0, -108],
589            [   0,   66, -121,  -74,  -53],
590            [   0,  104,    0,    0,    0],
591            [   0,   55,    0,    1,   -3],
592            [  93,    0,    0,    0,  104],
593        ]);
594    }
595
596    #[test]
597    fn test_6x6_1() {
598        test_matrix([
599            [   0,    0,    0,  -25,    0,    0],
600            [ -15,   79,    0,    0,    0,    0],
601            [   0,    0,    0,    0,    0,   14],
602            [   0,    0,    0, -114,  -61,    0],
603            [   0,    0,  109,    0,    0, -126],
604            [  46,    0,    0,   50,   21,    0],
605        ]);
606    }
607
608    #[test]
609    fn test_6x6_2() {
610        test_matrix([
611            [   0,    0,  -26,  -68,   84,    0],
612            [-125,   43,    0,    0,    0,  -63],
613            [   0,    0,    1,   90,    0,    0],
614            [   0,  -81,    0,    0,    0,    0],
615            [ -15,    0,    0,  -81,    0,    0],
616            [   0,  -12,    0,    0,    0,    1],
617        ]);
618    }
619
620    #[test]
621    fn test_10x10_1() {
622        test_matrix([
623            [   0,    0,    0,    0,    0,    0,  -60,    0,  -10,    0],
624            [   0,    0,    0,    0,    0,    0,    0,  -84,    0,    0],
625            [   0, -105,    0,    0,    0,    0,    0,    0,    0,    0],
626            [   0,    0,    0,  -25,    0,    0,    0,    0,  116,    0],
627            [   0,    0,    0,    0,  -18,    0,    0,    0,    0,    0],
628            [   0,    0,    0,  -72,    0,    0,    0,    0,    0,    0],
629            [   0,    0,   16,   48,    0,    0,    0,    0,    0,    0],
630            [ -57,    0,    0,  -88,  107,    0,    0,    0,    0,    0],
631            [-122, -108,    0,    0,    0,   91,    0,    0, -127,    0],
632            [   0,   85,    0,    0,  106,    0,    0,    0,    0, -121],
633        ]);
634    }
635
636    #[test]
637    fn test_11x11_1() {
638        test_matrix([
639            [   0,    0,    0,    0,    0,    0,    0,    0,    0,  -13,    0],
640            [   0,    0,    0,    0,    0,    0,    0,    0,    0,    0, -122],
641            [   0,    0,    0,    0,    0,  102,   82,    0,    0,    0,   13],
642            [   0,    0,    0,    0,    0, -107,  -39,    0,    0,    0,    0],
643            [   0,    0,    0,    0,    0,    0,  -39,   48, -113,    0,    0],
644            [  24,    0,    0,    0,    0,    0,    0,    0,    0,  -93, -120],
645            [-111,    0,    0,  -81,    0,    0,    0,    0,    0,    0,    0],
646            [   0,    0,    0,    0,    0,   82,    0,    0,   76,    0,    0],
647            [   0,    0,  -51,    0,    0,    0,  126,    0,    0,    0, -105],
648            [   0,  118,    0,    0,    0,    0,    0,    0,    0,    0,   27],
649            [   0,    0,  120,    0,  -31,    0,    0,    0,    0,    0,    0],
650        ]);
651    }
652
653    mod subtract_multiple_of_column_from_other_column {
654        use relp_num::R32;
655
656        use crate::algorithm::two_phase::tableau::inverse_maintenance::carry::lower_upper::decomposition::subtract_multiple_of_row_from_other_row;
657
658        #[test]
659        fn empty() {
660            let mut column1 = vec![];
661            let column2 = vec![];
662            subtract_multiple_of_row_from_other_row(&mut column1, &R32!(1), &column2);
663            assert_eq!(column1, vec![]);
664        }
665
666        #[test]
667        fn edited_empty() {
668            let mut column1 = vec![];
669            let column2 = vec![(1, R32!(1))];
670            subtract_multiple_of_row_from_other_row(&mut column1, &R32!(1), &column2);
671            assert_eq!(column1, vec![(1, R32!(-1))]);
672        }
673
674        #[test]
675        fn other_empty() {
676            let mut column1 = vec![(1, R32!(1))];
677            let column2 = vec![];
678            subtract_multiple_of_row_from_other_row(&mut column1, &R32!(1), &column2);
679            assert_eq!(column1, vec![(1, R32!(1))]);
680        }
681
682        #[test]
683        fn single_before() {
684            let mut column1 = vec![(1, R32!(1))];
685            let column2 = vec![(2, R32!(3))];
686            subtract_multiple_of_row_from_other_row(&mut column1, &R32!(1), &column2);
687            assert_eq!(column1, vec![(1, R32!(1)), (2, R32!(-3))]);
688        }
689
690        #[test]
691        fn single_at() {
692            let mut column1 = vec![(1, R32!(1))];
693            let column2 = vec![(1, R32!(3))];
694            subtract_multiple_of_row_from_other_row(&mut column1, &R32!(1), &column2);
695            assert_eq!(column1, vec![(1, R32!(-2))]);
696        }
697
698        #[test]
699        fn single_at_make_zero() {
700            let mut column1 = vec![(1, R32!(1))];
701            let column2 = vec![(1, R32!(3))];
702            subtract_multiple_of_row_from_other_row(&mut column1, &R32!(1, 3), &column2);
703            assert_eq!(column1, vec![]);
704        }
705
706        #[test]
707        fn single_after() {
708            let mut column1 = vec![(1, R32!(1))];
709            let column2 = vec![(0, R32!(3))];
710            subtract_multiple_of_row_from_other_row(&mut column1, &R32!(1), &column2);
711            assert_eq!(column1, vec![(0, R32!(-3)), (1, R32!(1))]);
712        }
713    }
714}
relp/algorithm/two_phase/tableau/inverse_maintenance/carry/lower_upper/decomposition/mod.rs

relp/algorithm/two_phase/tableau/inverse_maintenance/carry/lower_upper/decomposition/
mod.rs
