Reworked upper Hessenberg decomposition

benches/lib.rs

@@ -16,5 +16,6 @@ pub mod linalg {
 	mod norm;
 	mod triangular;
 	mod permutation;
+	mod hessenberg;
 	pub mod util;
 }

benches/linalg/hessenberg.rs

@@ -0,0 +1,73 @@
+use rulinalg::matrix::decomposition::{Decomposition, HessenbergDecomposition};
+
+use linalg::util::reproducible_random_matrix;
+
+use test::Bencher;
+
+#[bench]
+fn hessenberg_decomposition_decompose_100x100(b: &mut Bencher) {
+    let x = reproducible_random_matrix(100, 100);
+    b.iter(|| {
+        HessenbergDecomposition::decompose(x.clone())
+    })
+}
+
+#[bench]
+fn hessenberg_decomposition_decompose_unpack_100x100(b: &mut Bencher) {
+    let x = reproducible_random_matrix(100, 100);
+    b.iter(|| {
+        HessenbergDecomposition::decompose(x.clone()).unpack()
+    })
+}
+
+#[bench]
+fn hessenberg_decomposition_decompose_200x200(b: &mut Bencher) {
+    let x = reproducible_random_matrix(200, 200);
+    b.iter(|| {
+        HessenbergDecomposition::decompose(x.clone())
+    })
+}
+
+#[bench]
+fn hessenberg_decomposition_decompose_unpack_200x200(b: &mut Bencher) {
+    let x = reproducible_random_matrix(200, 200);
+    b.iter(|| {
+        HessenbergDecomposition::decompose(x.clone()).unpack()
+    })
+}
+
+#[bench]
+#[allow(deprecated)]
+fn upper_hessenberg_100x100(b: &mut Bencher) {
+    let x = reproducible_random_matrix(100, 100);
+    b.iter(|| {
+        x.clone().upper_hessenberg()
+    })
+}
+
+#[bench]
+#[allow(deprecated)]
+fn upper_hess_decomp_100x100(b: &mut Bencher) {
+    let x = reproducible_random_matrix(100, 100);
+    b.iter(|| {
+        x.clone().upper_hess_decomp()
+    })
+}
+
+#[bench]
+#[allow(deprecated)]
+fn upper_hessenberg_200x200(b: &mut Bencher) {
+    let x = reproducible_random_matrix(200, 200);
+    b.iter(|| {
+        x.clone().upper_hessenberg()
+    })
+}
+
+#[bench]
+#[allow(deprecated)]
+fn upper_hess_decomp_200x200(b: &mut Bencher) {
+    let x = reproducible_random_matrix(200, 200);
+    b.iter(|| {
+        x.clone().upper_hess_decomp()
+    })
+}

src/internal_utils.rs

@@ -1,6 +1,6 @@
 use matrix::{BaseMatrix, BaseMatrixMut};
 use libnum::{Zero, Num};
-use utils::in_place_vec_bin_op;
+use utils::{dot, in_place_vec_bin_op};
 
 pub fn nullify_lower_triangular_part<T, M>(matrix: &mut M)
     where T: Zero, M: BaseMatrixMut<T> {
@@ -20,6 +20,30 @@ pub fn nullify_upper_triangular_part<T, M>(matrix: &mut M)
     }
 }
 
+/// Given a vector `x` and a `m x n` matrix `A`, compute
+/// `y = A x`.
+///
+/// This is a stopgap solution until we have a more proper
+/// BLIS/BLAS-like API.
+pub fn gemv<T, M>(a: &M, x: &[T], y: &mut [T])
+    where M: BaseMatrix<T>, T: Num + Copy
+{
+    let m = a.rows();
+    let n = a.cols();
+
+    assert!(x.len() == n, "A and x must be dimensionally compatible.");
+    assert!(y.len() == m, "A and y must be dimensionally compatible.");
+
+    for element in y.iter_mut() {
+        *element = T::zero();
+    }
+
+    for i in 0 .. m {
+        let a_i = a.row(i).raw_slice();
+        y[i] = dot(a_i, x);
+    }
+}
+
 /// Given a vector `x` and a `m x n` matrix `A`, compute
 /// `y = A^T x`.
 ///
@@ -88,7 +112,7 @@ mod tests {
 
     use super::nullify_lower_triangular_part;
     use super::nullify_upper_triangular_part;
-    use super::transpose_gemv;
+    use super::{gemv, transpose_gemv};
     use super::ger;
 
     #[test]
@@ -117,6 +141,21 @@ mod tests {
         ]);
     }
 
+    #[test]
+    fn gemv_examples() {
+        {
+            let a = matrix![3.0, 4.0;
+                            5.0, 2.0;
+                            3.0, 1.0];
+            let x = vec![2.0, 3.0];
+            let mut y = vec![0.0; 3];
+            gemv(&a, &x, &mut y);
+
+            let y = Vector::new(y);
+            assert_vector_eq!(y, vector![18.0, 16.0, 9.0]);
+        }
+    }
+
     #[test]
     fn transpose_gemv_examples() {
         {

src/matrix/decomposition/eigen.rs

@@ -82,6 +82,7 @@ impl<T: Any + Float + Signed> Matrix<T> {
                       "Francis shift only works on matrices greater than 2x2.");
         debug_assert!(n == self.cols, "Matrix must be square for Francis shift.");
 
+        #[allow(deprecated)]
         let mut h = try!(self
             .upper_hessenberg()
             .map_err(|_| Error::new(ErrorKind::DecompFailure, "Could not compute eigenvalues.")));
@@ -223,6 +224,7 @@ impl<T: Any + Float + Signed> Matrix<T> {
                       "Francis shift only works on matrices greater than 2x2.");
         debug_assert!(n == self.cols, "Matrix must be square for Francis shift.");
 
+        #[allow(deprecated)]
         let (u, mut h) = try!(self.upper_hess_decomp().map_err(|_| {
             Error::new(ErrorKind::DecompFailure,
                        "Could not compute eigen decomposition.")