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

use impl_prelude::*;
use lapack::c::{slaswp, dlaswp, zlaswp, claswp};

/// A compressed representation of matrix permutations.
///
/// The internal representation of a matrix permutation mimics the
/// output of various LAPACK functions.
#[derive(Debug, Clone)]
pub struct MatrixPermutation {
    ipiv: Vec<i32>
}

impl MatrixPermutation {
    /// Return permutation from ipivot result from a LAPACK[E] call.
    pub fn from_ipiv(ipiv: Vec<i32>) -> MatrixPermutation {
        MatrixPermutation { ipiv: ipiv  }
    }

    /// Permute an input matrix `mat` by this permutation.
    pub fn permute_into<T, D1>(&self, mat: ArrayBase<D1, Ix2>) -> Option<ArrayBase<D1, Ix2>>
        where T: Permutes,
              D1: DataOwned<Elem=T> + DataMut<Elem=T> {
        Permutes::permute_into(mat, &self.ipiv)
    }

    /// Permute a matrix `mat` by this permutation.
    pub fn permute<T, D1>(&self, mat: &ArrayBase<D1, Ix2>) -> Option<Array<T, Ix2>>
        where T: Permutes, D1: Data<Elem=T> {
        Permutes::permute(mat, &self.ipiv)
    }

    /// Return the native LAPACKE representation of the permutation.
    ///
    /// This will only be useful for interfacing with other LAPACKE
    /// functions.
    pub fn ipiv(&self) -> &[i32] {
        &self.ipiv
    }
}

/// A `Permutes` is a scalar that one can apply a
/// `MatrixPermutation` to.
///
/// The methods of permuatble are not meant to be used by end-users.
pub trait Permutes: Sized + Clone {
    /// Permute a matrix by the pivot representation of a permutation.
    fn permute_into<D1>(mat: ArrayBase<D1, Ix2>, ipiv: &[i32]) -> Option<ArrayBase<D1, Ix2>>
        where D1: DataOwned<Elem=Self> + DataMut<Elem=Self>;

    /// Permute a matrix by the pivot representation of a permutation.
    fn permute<D1>(mat: &ArrayBase<D1, Ix2>, ipiv: &[i32]) -> Option<Array<Self, Ix2>>
        where D1: Data<Elem=Self> {
        Self::permute_into(mat.to_owned(), ipiv)
    }
}

macro_rules! impl_perm {
    ($perm_type:ty, $perm_func:ident) => (
        impl Permutes for $perm_type {
            fn permute_into<D1>(mut mat: ArrayBase<D1, Ix2>, ipiv: &[i32]) -> Option<ArrayBase<D1, Ix2>>
                where D1: DataOwned<Elem=Self> + DataMut<Elem=Self> {

                let (m, n) = mat.dim();

                if ipiv.len() > m {
                    return None
                }


                let info = {
                    let (mut slice, layout, lda) = match slice_and_layout_mut(&mut mat) {
                        None => return None,
                        Some(fwd) => fwd
                    };

                    $perm_func(layout, n as i32, slice, lda as i32, 1, ipiv.len() as i32,
                               ipiv, 1)
                };

                if info == 0 {
                    Some(mat)
                } else {
                    None
                }
            }
        }
    )
}

impl_perm!(f32, slaswp);
impl_perm!(f64, dlaswp);
impl_perm!(c32, claswp);
impl_perm!(c64, zlaswp);