M-Labs
/
nalgebra
3
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

generic implementation of QR decomposition 

But static matrices can't use it yet, they need to implement the
Row/Col slicing traits.
This commit is contained in:
Vincent Barrielle 2014-05-11 21:20:41 +02:00
parent 6ad11edf9b
commit dcf7b8ad01
1 changed files with 8 additions and 6 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

14
src/linalg/decompositions.rs
View File

@ -1,5 +1,4 @@
use std::num::{Zero, Float}; use std::num::{Zero, Float};
use na::DMat;
use traits::operations::Transpose; use traits::operations::Transpose;
use traits::structure::{ColSlice, Eye, Indexable}; use traits::structure::{ColSlice, Eye, Indexable};
use traits::geometry::Norm; use traits::geometry::Norm;
@ -34,11 +33,14 @@ fn householder_matrix<N: Float,
/// QR decomposition using Householder reflections /// QR decomposition using Householder reflections
/// # Arguments /// # Arguments
/// * `m` - matrix to decompose /// * `m` - matrix to decompose
pub fn decomp_qr<N: Clone + Float>(m: &DMat<N>) -> (DMat<N>, DMat<N>) { pub fn decomp_qr<N: Float,
let rows = m.nrows(); V: Indexable<uint, N> + Norm<N>,
let cols = m.ncols(); Mat: Clone + Eye + ColSlice<V> + Transpose
+ Indexable<(uint, uint), N> + Mul<Mat, Mat>>
(m: &Mat) -> (Mat, Mat) {
let (rows, cols) = m.shape();
assert!(rows >= cols); assert!(rows >= cols);
let mut q : DMat<N> = Eye::new_identity(rows); let mut q : Mat = Eye::new_identity(rows);
let mut r = m.clone(); let mut r = m.clone();
let iterations = min(rows - 1, cols); let iterations = min(rows - 1, cols);
@ -57,7 +59,7 @@ pub fn decomp_qr<N: Clone + Float>(m: &DMat<N>) -> (DMat<N>, DMat<N>) {
v.unsafe_set(0, x - alpha); v.unsafe_set(0, x - alpha);
} }
let _ = v.normalize(); let _ = v.normalize();
let qk: DMat<N> = householder_matrix(rows, 0, v); let qk: Mat = householder_matrix(rows, 0, v);
r = qk * r; r = qk * r;
q = q * Transpose::transpose_cpy(&qk); q = q * Transpose::transpose_cpy(&qk);
} }