algebra_linear 1.1.1

Pure-Rust BLAS, LAPACK, and sparse LU
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115

//! Permutation helpers.
//!
//! A permutation is a `[usize]` in which every index appears exactly once. The
//! routines here allocate their results and return [`SparseError`] rather than
//! panicking on malformed input.

use super::SparseError;

/// Returns `true` if every index in `0..perm.len()` appears exactly once.
///
/// The check allocates a `Vec<bool>`; it is not on a hot path.
pub fn is_permutation(perm: &[usize]) -> bool {
    let n = perm.len();
    let mut seen = vec![false; n];
    for &p in perm {
        if p >= n || seen[p] {
            return false;
        }
        seen[p] = true;
    }
    true
}

/// Returns the inverse of `perm`.
///
/// If `perm` maps a new position to an old index (`perm[new] = old`), the
/// inverse maps an old index to its new position (`inverse[old] = new`).
///
/// # Errors
///
/// [`SparseError::BadShape`] if `perm` is not a permutation.
pub fn invert(perm: &[usize]) -> Result<Vec<usize>, SparseError> {
    if !is_permutation(perm) {
        return Err(SparseError::BadShape);
    }
    let mut inverse = vec![0; perm.len()];
    for (new, &old) in perm.iter().enumerate() {
        inverse[old] = new;
    }
    Ok(inverse)
}

/// Composes two permutations: `compose(a, b)[i] = a[b[i]]`.
///
/// # Errors
///
/// [`SparseError::BadShape`] if the lengths differ or either argument is not a
/// permutation.
pub fn compose(a: &[usize], b: &[usize]) -> Result<Vec<usize>, SparseError> {
    if a.len() != b.len() || !is_permutation(a) || !is_permutation(b) {
        return Err(SparseError::BadShape);
    }
    Ok(b.iter().map(|&bi| a[bi]).collect())
}

/// Permutes `x` in place: `x[i] = old[perm[i]]`.
///
/// This is the direction the sparse solve uses for the row permutation, so it
/// matches the element order of a permuted right-hand side.
///
/// # Errors
///
/// [`SparseError::BadShape`] if `x.len() != perm.len()` or `perm` is not a
/// permutation.
pub fn permute<T: Copy>(x: &mut [T], perm: &[usize]) -> Result<(), SparseError> {
    if x.len() != perm.len() || !is_permutation(perm) {
        return Err(SparseError::BadShape);
    }
    let copy = x.to_vec();
    for (i, &p) in perm.iter().enumerate() {
        x[i] = copy[p];
    }
    Ok(())
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn detects_permutations() {
        assert!(is_permutation(&[]));
        assert!(is_permutation(&[0]));
        assert!(is_permutation(&[1, 0]));
        assert!(is_permutation(&[2, 0, 1]));
        assert!(!is_permutation(&[1, 1]));
        assert!(!is_permutation(&[0, 3]));
    }

    #[test]
    fn inverts_a_permutation() {
        let p = [2, 0, 1];
        let inv = invert(&p).unwrap();
        assert_eq!(inv, [1, 2, 0]);
        assert_eq!(invert(&[]).unwrap(), Vec::<usize>::new());
        assert_eq!(invert(&[0, 1, 1]).unwrap_err(), SparseError::BadShape);
    }

    #[test]
    fn composes_two_permutations() {
        let a = [2, 0, 1];
        let b = [1, 2, 0];
        assert_eq!(compose(&a, &b).unwrap(), [a[1], a[2], a[0]]);
        assert_eq!(compose(&a, &[0, 0, 1]).unwrap_err(), SparseError::BadShape);
        assert_eq!(compose(&a, &[0, 1]).unwrap_err(), SparseError::BadShape);
    }

    #[test]
    fn permutes_a_vector() {
        let mut x = [10, 20, 30];
        permute(&mut x, &[2, 0, 1]).unwrap();
        assert_eq!(x, [30, 10, 20]);
        assert_eq!(permute(&mut x, &[0, 1]).unwrap_err(), SparseError::BadShape);
    }
}