ndarray_linalg/
qr.rs

1//! QR decomposition
2//!
3//! [Wikipedia article on QR decomposition](https://en.wikipedia.org/wiki/QR_decomposition)
4
5use ndarray::*;
6use num_traits::Zero;
7
8use crate::convert::*;
9use crate::error::*;
10use crate::layout::*;
11use crate::triangular::*;
12use crate::types::*;
13
14pub use lax::UPLO;
15
16/// QR decomposition for matrix reference
17///
18/// This creates copy due for reshaping array.
19/// To avoid copy and the matrix is square, please use `QRSquare*` traits.
20pub trait QR {
21    type Q;
22    type R;
23    fn qr(&self) -> Result<(Self::Q, Self::R)>;
24}
25
26/// QR decomposition
27///
28/// This creates copy due for reshaping array.
29/// To avoid copy and the matrix is square, please use `QRSquare*` traits.
30pub trait QRInto: Sized {
31    type Q;
32    type R;
33    fn qr_into(self) -> Result<(Self::Q, Self::R)>;
34}
35
36/// QR decomposition for square matrix reference
37pub trait QRSquare: Sized {
38    type Q;
39    type R;
40    fn qr_square(&self) -> Result<(Self::Q, Self::R)>;
41}
42
43/// QR decomposition for square matrix
44pub trait QRSquareInto: Sized {
45    type R;
46    fn qr_square_into(self) -> Result<(Self, Self::R)>;
47}
48
49/// QR decomposition for mutable reference of square matrix
50pub trait QRSquareInplace: Sized {
51    type R;
52    fn qr_square_inplace(&mut self) -> Result<(&mut Self, Self::R)>;
53}
54
55impl<A, S> QRSquareInplace for ArrayBase<S, Ix2>
56where
57    A: Scalar + Lapack,
58    S: DataMut<Elem = A>,
59{
60    type R = Array2<A>;
61
62    fn qr_square_inplace(&mut self) -> Result<(&mut Self, Self::R)> {
63        let l = self.square_layout()?;
64        let r = A::qr(l, self.as_allocated_mut()?)?;
65        let r: Array2<_> = into_matrix(l, r)?;
66        Ok((self, r.into_triangular(UPLO::Upper)))
67    }
68}
69
70impl<A, S> QRSquareInto for ArrayBase<S, Ix2>
71where
72    A: Scalar + Lapack,
73    S: DataMut<Elem = A>,
74{
75    type R = Array2<A>;
76
77    fn qr_square_into(mut self) -> Result<(Self, Self::R)> {
78        let (_, r) = self.qr_square_inplace()?;
79        Ok((self, r))
80    }
81}
82
83impl<A, S> QRSquare for ArrayBase<S, Ix2>
84where
85    A: Scalar + Lapack,
86    S: Data<Elem = A>,
87{
88    type Q = Array2<A>;
89    type R = Array2<A>;
90
91    fn qr_square(&self) -> Result<(Self::Q, Self::R)> {
92        let a = self.to_owned();
93        a.qr_square_into()
94    }
95}
96
97impl<A, S> QRInto for ArrayBase<S, Ix2>
98where
99    A: Scalar + Lapack,
100    S: DataMut<Elem = A>,
101{
102    type Q = Array2<A>;
103    type R = Array2<A>;
104
105    fn qr_into(mut self) -> Result<(Self::Q, Self::R)> {
106        let n = self.nrows();
107        let m = self.ncols();
108        let k = ::std::cmp::min(n, m);
109        let l = self.layout()?;
110        let r = A::qr(l, self.as_allocated_mut()?)?;
111        let r: Array2<_> = into_matrix(l, r)?;
112        let q = self;
113        Ok((take_slice(&q, n, k), take_slice_upper(&r, k, m)))
114    }
115}
116
117impl<A, S> QR for ArrayBase<S, Ix2>
118where
119    A: Scalar + Lapack,
120    S: Data<Elem = A>,
121{
122    type Q = Array2<A>;
123    type R = Array2<A>;
124
125    fn qr(&self) -> Result<(Self::Q, Self::R)> {
126        let a = self.to_owned();
127        a.qr_into()
128    }
129}
130
131fn take_slice<A, S1, S2>(a: &ArrayBase<S1, Ix2>, n: usize, m: usize) -> ArrayBase<S2, Ix2>
132where
133    A: Copy,
134    S1: Data<Elem = A>,
135    S2: DataMut<Elem = A> + DataOwned,
136{
137    let av = a.slice(s![..n as isize, ..m as isize]);
138    replicate(&av)
139}
140
141fn take_slice_upper<A, S1, S2>(a: &ArrayBase<S1, Ix2>, n: usize, m: usize) -> ArrayBase<S2, Ix2>
142where
143    A: Copy + Zero,
144    S1: Data<Elem = A>,
145    S2: DataMut<Elem = A> + DataOwned,
146{
147    let av = a.slice(s![..n, ..m]);
148    let mut a = replicate(&av);
149    Zip::indexed(&mut a).for_each(|(i, j), elt| {
150        if i > j {
151            *elt = A::zero()
152        }
153    });
154    a
155}
ndarray_linalg/qr.rs

ndarray_linalg/
qr.rs
