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
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149

//! Solve singular value decomposition (SVD) of arbitrary matrices.

use lapack::c::{sgesvd, sgesdd, dgesvd, dgesdd, cgesvd, cgesdd, zgesvd, zgesdd};
use super::types::{SVDSolution, SVDError, SingularValue};
use impl_prelude::*;

const SVD_NORMAL_LIMIT: usize = 200;

/// Trait for scalars that can implement SVD.
pub trait SVD<SV: SingularValue>: Sized + Clone {
    /// Compute the singular value decomposition of a matrix.
    ///
    /// Use `Self::compute` when you don't wnat to consume the input
    /// matrix.
    ///
    /// On success, returns an `SVDSolution`, which always contains the
    /// singular values and optionally contains the left and right
    /// singular vectors. The left vectors (via the matrix `u`) are
    /// returned iff `compute_u` is true, and similarly for `vt` and
    /// `compute_vt`.
    fn compute_into<D>(mat: ArrayBase<D, Ix2>,
                       compute_u: bool,
                       compute_vt: bool)
                      -> Result<SVDSolution<Self, SV>, SVDError>
        where D: DataMut<Elem = Self> + DataOwned<Elem = Self>;

    /// Comptue the singular value decomposition of a matrix.
    ///
    /// Similar to [`SVD::compute_into`](#tymethod.compute_into), but
    /// the values are copied beforehand. leaving the original matrix
    /// un-modified.
    fn compute<D>(mat: &ArrayBase<D, Ix2>,
                  compute_u: bool,
                  compute_vt: bool)
                  -> Result<SVDSolution<Self, SV>, SVDError>
        where D: Data<Elem = Self>
    {
        let vec: Vec<Self> = mat.iter().cloned().collect();
        let m = Array::from_shape_vec(mat.dim(), vec).unwrap();
        Self::compute_into(m, compute_u, compute_vt)
    }
}


#[derive(Debug, PartialEq)]
enum SVDMethod {
    Normal,
    DivideAndConquer,
}


/// Choose a method based on the problem.
fn select_svd_method(d: &Ix2, compute_either: bool) -> SVDMethod {
    let mx = cmp::max(d.0, d.1);

    // When we're computing one of them singular vector sets, we have
    // to compute both with divide and conquer. So, we're bound by the
    // maximum size of the array.
    if compute_either {
        if mx > SVD_NORMAL_LIMIT {
            SVDMethod::Normal
        } else {
            SVDMethod::DivideAndConquer
        }
    } else {
        SVDMethod::DivideAndConquer
    }
}


macro_rules! impl_svd {
    ($impl_type:ident, $sv_type:ident, $svd_func:ident, $sdd_func:ident) => (
        impl SVD<$sv_type> for $impl_type {

            fn compute_into<D>(mut mat: ArrayBase<D, Ix2>,
                               mut compute_u: bool,
                               mut compute_vt: bool)
                               -> Result<SVDSolution<$impl_type, $sv_type>, SVDError>
                where D: DataMut<Elem=Self> + DataOwned<Elem = Self>{

                let dim = mat.dim();
                let (m, n) = dim;
                let mut s = Array::default(cmp::min(m, n));

                let (slice, layout, lda) = match slice_and_layout_mut(&mut mat) {
                    Some(x) => x,
                    None => return Err(SVDError::BadLayout)
                };

                let compute_either = compute_u || compute_vt;
                let method = select_svd_method(&dim, compute_either);
                if method == SVDMethod::DivideAndConquer {
                    compute_u = compute_either;
                    compute_vt = compute_either;
                }

                let mut u = matrix_with_layout(if compute_u { (m, m) } else { (0, 0) }, layout);
                let mut vt = matrix_with_layout(if compute_vt { (n, n) } else { (0, 0) }, layout);

                let job_u = if compute_u { b'A' } else { b'N' };
                let job_vt = if compute_vt { b'A' } else { b'N' };

                let info = match method {
                    SVDMethod::Normal => {
                        let mut superb = Array::default(cmp::min(m, n) - 2);

                        $svd_func(layout, job_u, job_vt, m as i32, n as i32, slice,
                                  lda as i32, s.as_slice_mut().expect("bad s implementation"),
                                  u.as_slice_mut().expect("bad u implementation"), m as i32,
                                  vt.as_slice_mut().expect("bad vt implementation"), n as i32,
                                  superb.as_slice_mut().expect("bad superb implementation"))
                    },
                    SVDMethod::DivideAndConquer => {
                        let job_z = if compute_u || compute_vt { b'A' } else { b'N' };
                        $sdd_func(layout, job_z, m as i32, n as i32, slice,
                                  lda as i32,
                                  s.as_slice_mut().expect("bad s implementation"),
                                  u.as_slice_mut().expect("bad u implementation"), m as i32,
                                  vt.as_slice_mut().expect("bad vt implementation"), n as i32)
                    }
                };

                match info {
                    0 => {
                        Ok(SVDSolution {
                            values: s,
                            left_vectors: if compute_u { Some(u) } else { None },
                            right_vectors: if compute_vt { Some(vt) } else { None }
                        })
                    },
                    x if x < 0 => {
                        Err(SVDError::IllegalParameter(-x - 1))
                    },
                    x if x > 0 => {
                        Err(SVDError::Unconverged)
                    },
                    _ => {
                        unreachable!();
                    }
                }
            }
        }
    )
}

impl_svd!(f32, f32, sgesvd, sgesdd);
impl_svd!(f64, f64, dgesvd, dgesdd);
impl_svd!(c32, f32, cgesvd, cgesdd);
impl_svd!(c64, f64, zgesvd, zgesdd);