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

//! Compute eigenvalues and eigenvectors of symmetric matrices.
//!
//! Symmetric (or Hermitian, for complex) matrices are guaranteed to
//! have real eigenvalues.

use lapack::c::{ssyev, dsyev, cheev, zheev};
use super::types::{Solution, EigenError};
use impl_prelude::*;

/// Scalar trait for computing eigenvalues of a symmetric matrix.
///
/// In order to extract eigenvalues or eigenvectors from a matrix,
/// that matrix with must have  entries implementing the `Eigen` trait.
pub trait SymEigen: Sized {
    /// Scalar type of the resulting eigenvalues.
    type Eigenvalue;

    /// Solution structure from the eigenvalue computation.
    type Solution;

    /// Return the real eigenvalues of a symmetric matrix.
    ///
    /// If `with_vectors` is true, the right eigenvectors of 'V' are
    /// stored in the input matrix.
    fn compute_mut<D>(mat: &mut ArrayBase<D, Ix2>,
                      uplo: Symmetric,
                      with_vectors: bool)
                      -> Result<Array<Self::Eigenvalue, Ix1>, EigenError>
        where D: DataMut<Elem = Self>;

    /// Return the real eigenvalues of a symmetric matrix.
    fn compute_into<D>(mut mat: ArrayBase<D, Ix2>,
                       uplo: Symmetric)
                       -> Result<Array<Self::Eigenvalue, Ix1>, EigenError>
        where D: DataMut<Elem = Self> + DataOwned<Elem = Self> {
        Self::compute_mut(&mut mat, uplo, false)
    }

    /// Return the eigenvalues and, optionally, the eigenvectors of a
    /// symmetric matrix.
    ///
    /// # Remarks
    ///
    /// The input matrix is copied before the calculation takes place.
    fn compute<D>(mat: &ArrayBase<D, Ix2>,
                  uplo: Symmetric,
                  with_vectors: bool)
                  -> Result<Self::Solution, EigenError>
        where D: Data<Elem = Self>;
}

macro_rules! impl_sym_eigen {
    ($impl_type:ident, $eigen_type:ident, $func:ident) => (
        impl SymEigen for $impl_type {
            type Eigenvalue = $eigen_type;
            type Solution = Solution<Self, Self::Eigenvalue>;

            fn compute_mut<D>(mat: &mut ArrayBase<D, Ix2>, uplo: Symmetric, with_vectors: bool) ->
                Result<Array<Self::Eigenvalue, Ix1>, EigenError>
                where D: DataMut<Elem=Self>
            {
                let dim = mat.dim();
                if dim.0 != dim.1 {
                    return Err(EigenError::NotSquare);
                }

                let n = dim.0 as i32;

                let (mut data_slice, layout, ld) = match slice_and_layout_mut(mat) {
                    Some(x) => x,
                    None => return Err(EigenError::BadLayout)
                };

                let mut values = Array::default(n as Ix);
                let job = if with_vectors { b'V' } else { b'N' };

                let info = $func(layout, job, uplo as u8, n, data_slice,
                                 ld as i32, values.as_slice_mut().unwrap());

                if info  == 0 {
                    Ok(values)
                } else if info < 0 {
                    Err(EigenError::IllegalParameter(-info))
                } else {
                    Err(EigenError::Failed)
                }
            }

            fn compute<D>(mat: &ArrayBase<D, Ix2>,
                          uplo: Symmetric,
                          with_vectors: bool) -> Result<Self::Solution, EigenError>
                where D: Data<Elem=Self> {
                let vec: Vec<Self> = mat.iter().cloned().collect();
                let mut new_mat = Array::from_shape_vec(mat.dim(), vec).unwrap();
                let r = Self::compute_mut(&mut new_mat, uplo, with_vectors);
                r.map(|values| {
                    Solution {
                        values: values,
                        left_vectors: None,
                        right_vectors: if with_vectors { Some(new_mat) } else { None }
                    }
                })
            }
        }
    )
}

impl_sym_eigen!(f32, f32, ssyev);
impl_sym_eigen!(c32, f32, cheev);
impl_sym_eigen!(f64, f64, dsyev);
impl_sym_eigen!(c64, f64, zheev);

#[cfg(test)]
mod tests {}