ndarray-linalg 0.1.3

Linear algebra package for rust-ndarray using LAPACK
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

//! Implement QR decomposition

use std::cmp::min;
use lapack::fortran::*;
use num_traits::Zero;

use error::LapackError;

pub trait ImplQR: Sized {
    fn qr(n: usize, m: usize, mut a: Vec<Self>) -> Result<(Vec<Self>, Vec<Self>), LapackError>;
    fn lq(n: usize, m: usize, mut a: Vec<Self>) -> Result<(Vec<Self>, Vec<Self>), LapackError>;
}

macro_rules! impl_qr {
    ($geqrf:path, $orgqr:path, $gelqf:path, $orglq:path) => {
// XXX These codes are most same, but the argument of $orgqr and $orglq are different!
fn qr(n: usize, m: usize, mut a: Vec<Self>) -> Result<(Vec<Self>, Vec<Self>), LapackError> {
    let n = n as i32;
    let m = m as i32;
    let mut info = 0;
    let k = min(m, n);
    let lda = m;
    let lw_default = 1000;
    let mut tau = vec![Self::zero(); k as usize];
    let mut work = vec![Self::zero(); lw_default];
// estimate lwork
    $geqrf(m, n, &mut a, lda, &mut tau, &mut work, -1, &mut info);
    let lwork_r = work[0] as i32;
    if lwork_r > lw_default as i32 {
        work = vec![Self::zero(); lwork_r as usize];
    }
// calc R
    $geqrf(m, n, &mut a, lda, &mut tau, &mut work, lwork_r, &mut info);
    if info != 0 {
        return Err(From::from(info));
    }
    let r = a.clone();
// re-estimate lwork
    $orgqr(m, k, k, &mut a, lda, &mut tau, &mut work, -1, &mut info);
    let lwork_q = work[0] as i32;
    if lwork_q > lwork_r {
        work = vec![Self::zero(); lwork_q as usize];
    }
// calc Q
    $orgqr(m,
           k,
           k,
           &mut a,
           lda,
           &mut tau,
           &mut work,
           lwork_q,
           &mut info);
    if info == 0 {
        Ok((a, r))
    } else {
        Err(From::from(info))
    }
}
fn lq(n: usize, m: usize, mut a: Vec<Self>) -> Result<(Vec<Self>, Vec<Self>), LapackError> {
    let n = n as i32;
    let m = m as i32;
    let mut info = 0;
    let k = min(m, n);
    let lda = m;
    let lw_default = 1000;
    let mut tau = vec![Self::zero(); k as usize];
    let mut work = vec![Self::zero(); lw_default];
// estimate lwork
    $gelqf(m, n, &mut a, lda, &mut tau, &mut work, -1, &mut info);
    let lwork_r = work[0] as i32;
    if lwork_r > lw_default as i32 {
        work = vec![Self::zero(); lwork_r as usize];
    }
// calc R
    $gelqf(m, n, &mut a, lda, &mut tau, &mut work, lwork_r, &mut info);
    if info != 0 {
        return Err(From::from(info));
    }
    let r = a.clone();
// re-estimate lwork
    $orglq(k, n, k, &mut a, lda, &mut tau, &mut work, -1, &mut info);
    let lwork_q = work[0] as i32;
    if lwork_q > lwork_r {
        work = vec![Self::zero(); lwork_q as usize];
    }
// calc Q
    $orglq(k,
           n,
           k,
           &mut a,
           lda,
           &mut tau,
           &mut work,
           lwork_q,
           &mut info);
    if info == 0 {
        Ok((a, r))
    } else {
        Err(From::from(info))
    }
}
}} // endmacro

impl ImplQR for f64 {
    impl_qr!(dgeqrf, dorgqr, dgelqf, dorglq);
}

impl ImplQR for f32 {
    impl_qr!(sgeqrf, sorgqr, sgelqf, sorglq);
}