ldpc_toolbox/
linalg.rs

1use ndarray::{Array2, LinalgScalar, s};
2
3#[derive(Debug, Copy, Clone, Eq, PartialEq, Hash)]
4pub enum Error {
5    NotInvertible,
6}
7
8pub fn gauss_reduction<A: LinalgScalar + PartialEq>(array: &mut Array2<A>) -> Result<(), Error> {
9    let (n, m) = array.dim();
10    assert!(n <= m);
11
12    // Reduce to upper triangular with ones on diagonal
13    for j in 0..n {
14        // Find non-zero element in current column
15        let Some(k) = array
16            .slice(s![j.., j])
17            .iter()
18            .enumerate()
19            .find_map(|(t, x)| if x.is_zero() { None } else { Some(j + t) })
20        else {
21            return Err(Error::NotInvertible);
22        };
23
24        if k != j {
25            // Swap rows j and k
26            for t in j..m {
27                array.swap([j, t], [k, t]);
28            }
29        }
30
31        // Make a 1 by dividing
32        let x = array[[j, j]];
33        if !x.is_one() {
34            for t in j..m {
35                array[[j, t]] = array[[j, t]] / x;
36            }
37        }
38
39        // Subtract to rows below to make zeros below diagonal
40        for t in (j + 1)..n {
41            let x = array[[t, j]];
42            if !x.is_zero() {
43                // avoid calculations if we're subtracting zero
44                for u in j..m {
45                    array[[t, u]] = array[[t, u]] - x * array[[j, u]];
46                }
47            }
48        }
49    }
50
51    // Reduce to identity
52    for j in (0..n).rev() {
53        // Subtract to rows above to make zeros above diagonal
54        for t in 0..j {
55            let x = array[[t, j]];
56            if !x.is_zero() {
57                // avoid calculations if we're subtracting zero
58                for u in j..m {
59                    array[[t, u]] = array[[t, u]] - x * array[[j, u]];
60                }
61            }
62        }
63    }
64
65    Ok(())
66}
67
68pub fn row_echelon_form<A: LinalgScalar + PartialEq>(array: &mut Array2<A>) {
69    let (n, m) = array.dim();
70
71    let mut j = 0;
72    let mut k = 0;
73    while j < m && k < n {
74        // Find non-zero element in current column, at or below row k
75        let Some(s) = array
76            .slice(s![k.., j])
77            .iter()
78            .enumerate()
79            .find_map(|(t, x)| if x.is_zero() { None } else { Some(k + t) })
80        else {
81            // All the elements at or below row k are zero. Done with this
82            // column.
83            j += 1;
84            continue;
85        };
86
87        if s != k {
88            // Swap rows s and k
89            for t in j..m {
90                array.swap([s, t], [k, t]);
91            }
92        }
93
94        let x = array[[k, j]];
95
96        // Subtract to rows below to make zeros below row k
97        for t in (k + 1)..n {
98            let y = array[[t, j]];
99            if !y.is_zero() {
100                // avoid calculations if we're subtracting zero
101                for u in j..m {
102                    array[[t, u]] = array[[t, u]] - y * array[[k, u]] / x;
103                }
104            }
105        }
106
107        j += 1;
108        k += 1;
109    }
110}
111
112#[cfg(test)]
113mod test {
114    use super::*;
115    use crate::gf2::GF2;
116    use ndarray::arr2;
117    use num_traits::{One, Zero};
118
119    #[test]
120    fn gauss() {
121        let i = GF2::one();
122        let o = GF2::zero();
123        let mut a = arr2(&[
124            [i, o, i, i, i, o, i, o, i],
125            [i, i, o, o, i, i, o, i, o],
126            [i, i, i, o, o, i, i, o, i],
127        ]);
128        gauss_reduction(&mut a).unwrap();
129        let expected = arr2(&[
130            [i, o, o, i, o, o, o, i, o],
131            [o, i, o, i, i, i, o, o, o],
132            [o, o, i, o, i, o, i, i, i],
133        ]);
134        assert_eq!(&a, &expected);
135    }
136
137    #[test]
138    fn row_echelon() {
139        let i = GF2::one();
140        let o = GF2::zero();
141        let mut a = arr2(&[
142            [i, i, o, o, i, o, i, o, i],
143            [i, o, o, i, i, i, o, i, o],
144            [i, i, o, o, o, i, i, o, i],
145        ]);
146        row_echelon_form(&mut a);
147        let expected = arr2(&[
148            [i, i, o, o, i, o, i, o, i],
149            [o, i, o, i, o, i, i, i, i],
150            [o, o, o, o, i, i, o, o, o],
151        ]);
152        assert_eq!(&a, &expected);
153    }
154}
ldpc_toolbox/linalg.rs

ldpc_toolbox/
linalg.rs
