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


use std::cmp::min;
use std::fmt::Debug;
use ndarray::LinalgScalar;
use num_traits::float::Float;

use error::LapackError;
use binding;

pub trait TruncatableFloat: Float {
    fn to_int(self) -> i32;
}

impl TruncatableFloat for f64 {
    fn to_int(self) -> i32 {
        self as i32
    }
}
impl TruncatableFloat for f32 {
    fn to_int(self) -> i32 {
        self as i32
    }
}

pub trait LapackScalar
    : Debug + TruncatableFloat + LinalgScalar + binding::LapackBinding {
    fn eigh(n: usize, mut a: Vec<Self>) -> Result<(Vec<Self>, Vec<Self>), LapackError> {
        let mut w = vec![Self::zero(); n];
        let mut work = vec![Self::zero(); 4 * n];
        let mut info = 0;
        Self::_syev(b'V',
                    b'U',
                    n as i32,
                    &mut a,
                    n as i32,
                    &mut w,
                    &mut work,
                    4 * n as i32,
                    &mut info);
        if info == 0 {
            Ok((w, a))
        } else {
            Err(From::from(info))
        }
    }
    fn inv(size: usize, mut a: Vec<Self>) -> Result<Vec<Self>, LapackError> {
        let n = size as i32;
        let lda = n;
        let mut ipiv = vec![0; size];
        let mut info = 0;
        Self::_getrf(n, n, &mut a, lda, &mut ipiv, &mut info);
        if info != 0 {
            return Err(From::from(info));
        }
        let lwork = n;
        let mut work = vec![Self::zero(); size];
        Self::_getri(n, &mut a, lda, &mut ipiv, &mut work, lwork, &mut info);
        if info == 0 {
            Ok(a)
        } else {
            Err(From::from(info))
        }
    }
    fn norm_1(m: usize, n: usize, mut a: Vec<Self>) -> Self {
        let mut work = Vec::<Self>::new();
        Self::_lange(b'o', m as i32, n as i32, &mut a, m as i32, &mut work)
    }
    fn norm_i(m: usize, n: usize, mut a: Vec<Self>) -> Self {
        let mut work = vec![Self::zero(); m];
        Self::_lange(b'i', m as i32, n as i32, &mut a, m as i32, &mut work)
    }
    fn norm_f(m: usize, n: usize, mut a: Vec<Self>) -> Self {
        let mut work = Vec::<Self>::new();
        Self::_lange(b'f', m as i32, n as i32, &mut a, m as i32, &mut work)
    }
    fn svd(n: usize,
           m: usize,
           mut a: Vec<Self>)
           -> Result<(Vec<Self>, Vec<Self>, Vec<Self>), LapackError> {
        let mut info = 0;
        let n = n as i32;
        let m = m as i32;
        let lda = m;
        let ldu = m;
        let ldvt = n;
        let lwmax = 1000; // XXX
        let lwork = -1;
        let mut u = vec![Self::zero(); (ldu * m) as usize];
        let mut vt = vec![Self::zero(); (ldvt * n) as usize];
        let mut s = vec![Self::zero(); n as usize];
        let mut work = vec![Self::zero(); lwmax];
        Self::_gesvd('A' as u8,
                     'A' as u8,
                     m,
                     n,
                     &mut a,
                     lda,
                     &mut s,
                     &mut u,
                     ldu,
                     &mut vt,
                     ldvt,
                     &mut work,
                     lwork,
                     &mut info); // calc optimal work
        let lwork = min(lwmax as i32, work[0].to_int());
        Self::_gesvd('A' as u8,
                     'A' as u8,
                     m,
                     n,
                     &mut a,
                     lda,
                     &mut s,
                     &mut u,
                     ldu,
                     &mut vt,
                     ldvt,
                     &mut work,
                     lwork,
                     &mut info);
        if info == 0 {
            Ok((u, s, vt))
        } else {
            Err(From::from(info))
        }
    }
}

impl LapackScalar for f64 {}
impl LapackScalar for f32 {}