diff --git a/benches/lib.rs b/benches/lib.rs index f693a75..d2a6e11 100644 --- a/benches/lib.rs +++ b/benches/lib.rs @@ -16,5 +16,6 @@ pub mod linalg { mod norm; mod triangular; mod permutation; + mod hessenberg; pub mod util; } diff --git a/benches/linalg/hessenberg.rs b/benches/linalg/hessenberg.rs new file mode 100644 index 0000000..3086f91 --- /dev/null +++ b/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() + }) +} diff --git a/src/internal_utils.rs b/src/internal_utils.rs index dd2818f..5a4d8f5 100644 --- a/src/internal_utils.rs +++ b/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(matrix: &mut M) where T: Zero, M: BaseMatrixMut { @@ -20,6 +20,30 @@ pub fn nullify_upper_triangular_part(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(a: &M, x: &[T], y: &mut [T]) + where M: BaseMatrix, 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() { { diff --git a/src/matrix/decomposition/eigen.rs b/src/matrix/decomposition/eigen.rs index eb328f2..b6466e6 100644 --- a/src/matrix/decomposition/eigen.rs +++ b/src/matrix/decomposition/eigen.rs @@ -82,6 +82,7 @@ impl Matrix { "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 Matrix { "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.") diff --git a/src/matrix/decomposition/hessenberg.rs b/src/matrix/decomposition/hessenberg.rs index b18ff35..ef1366e 100644 --- a/src/matrix/decomposition/hessenberg.rs +++ b/src/matrix/decomposition/hessenberg.rs @@ -1,9 +1,171 @@ use matrix::{Matrix, BaseMatrix, BaseMatrixMut, MatrixSlice, MatrixSliceMut}; +use matrix::decomposition::{Decomposition, HouseholderReflection, + HouseholderComposition}; +use matrix::decomposition::householder; use error::{Error, ErrorKind}; +use vector::Vector; use std::any::Any; -use libnum::{Float}; +use libnum::{Zero, Float}; + +/// Result of unpacking a +/// [Hessenberg decomposition](struct.HessenbergDecomposition.html). +#[derive(Debug, Clone)] +pub struct QH { + /// The orthogonal factor `Q` in the decomposition. + pub q: Matrix, + /// The upper Hessenberg factor `H` in the decomposition. + pub h: Matrix +} + +/// Upper Hessenberg decomposition of real square matrices. +/// +/// Given any real square `m x m` matrix `A`, there exists an orthogonal +/// `m x m` matrix `Q` and an `m x m` *upper* Hessenberg matrix `H` such that +/// +/// ```text +/// A = Q H Qᵀ. (1) +/// ``` +/// +/// An Upper Hessenberg matrix `H` is characterized by the property that +/// `H[i, j] = 0` for any `i > j + 1`. For example, for a 4x4 matrix `H`, this +/// means that `H` has a the following pattern of zeros and (possibly) +/// non-zeros `x`: +/// +/// ```text +/// [ x x x x ] +/// [ x x x x ] +/// [ 0 x x x ] +/// [ 0 0 x x ] +/// ``` +/// +/// Hessenberg matrices are of particular interest because their eigenvalues +/// and eigenvectors can be computed much more efficiently than for a general +/// non-symmetric matrix, and because the Hessenberg matrix `H` corresponding +/// to any matrix `A` is unitarily similar to `A`, the eigenvalues of `H` +/// coincide with the eigenvalues of `A`. +/// +/// As in the case of [HouseholderQr](struct.HouseholderQr.html), +/// the `Q` factor is internally stored implicitly in terms of a sequence of +/// Householder transformations. The implicit `Q` matrix can be accessed by calling +/// the `.q()` method on the `HessenbergDecomposition` instance. +/// +/// # Examples +/// +/// ```rust +/// # #[macro_use] extern crate rulinalg; fn main() { +/// use rulinalg::matrix::decomposition::{Decomposition, HessenbergDecomposition, QH}; +/// +/// let x = matrix![5.0, 2.0, 0.0, -2.0; +/// 3.0, 3.0, 4.0, 3.0; +/// -3.0, 1.0, 2.0, 1.0; +/// 2.0, 3.0, 2.0, -3.0]; +/// +/// // We only clone `x` here so that we can compare for equality later +/// let QH { q, h } = HessenbergDecomposition::decompose(x.clone()).unpack(); +/// +/// assert_matrix_eq!(x, &q * h * q.transpose(), comp = float, ulp = 100); +/// # } +/// ``` +/// +/// # Internal storage format +/// +/// The internal storage format is very similar to that of +/// [HouseholderQr](struct.HouseholderQr.html). Each Householder reflector +/// is stored compactly in the elements of the matrix which correspond +/// to zeros in the Heissenberg factor `H`. In addition, a vector of length +/// `n - 1` is allocated to hold the multipliers for the Householder vectors. +#[derive(Debug, Clone)] +pub struct HessenbergDecomposition { + qh: Matrix, + tau: Vec +} + +impl HessenbergDecomposition where T: Float { + /// Computes the upper Hessenberg decomposition + /// for a given square matrix. + /// + /// # Panics + /// + /// - The matrix is not square. + pub fn decompose(matrix: Matrix) -> HessenbergDecomposition { + assert!(matrix.rows() == matrix.cols(), "Matrix must be square for Hessenberg decomposition."); + let n = matrix.rows(); + + // n - 1 is the number of elements along the sub-1 diagonal + // (which characterizes the Hessenberg matrix) + let n1 = n.saturating_sub(1); + + let mut qh = matrix; + let mut tau = vec![T::zero(); n1]; + + let mut buffer = vec![T::zero(); n1]; + let mut multiply_buffer = vec![T::zero(); n]; + + for j in 0 .. n1 { + buffer.truncate(n1 - j); + + // Compute the Householder matrix Q_j and left-apply it to the + // bottom right corner (we don't need to consider the columns + // before j, since they are implicitly zero) + let house = { + let mut bottom_right = qh.sub_slice_mut([j + 1, j], n1 - j, n - j); + bottom_right.col(0).clone_into_slice(&mut buffer); + + let house = HouseholderReflection::compute(Vector::new(buffer)); + house.buffered_left_multiply_into(&mut bottom_right, &mut multiply_buffer); + house.store_in_col(&mut bottom_right.col_mut(0)); + house + }; + + // Apply Q_j to the trailing columns + let mut trailing_cols = qh.sub_slice_mut([0, j + 1], n, n1 - j); + house.buffered_right_multiply_into(&mut trailing_cols, &mut multiply_buffer); + + tau[j] = house.tau(); + buffer = house.into_vector().into_vec(); + + } + + HessenbergDecomposition { + qh: qh, + tau: tau + } + } + + /// Returns the orthogonal factor Q as an instance + /// of a [HouseholderComposition](struct.HouseholderComposition.html) + /// operator. + pub fn q(&self) -> HouseholderComposition { + householder::create_composition(&self.qh, &self.tau, 1) + } +} + +impl Decomposition for HessenbergDecomposition where T: Float { + type Factors = QH; + + fn unpack(self) -> Self::Factors { + let q = self.q().first_k_columns(self.qh.rows()); + let mut h = self.qh; + nullify_lower_portion(&mut h); + + QH { + q: q, + h: h + } + } +} + +/// Makes the portion of the matrix which is not part of the +/// Upper Hessenberg part identically zero. +fn nullify_lower_portion(matrix: &mut Matrix) where T: Zero { + for (i, mut row) in matrix.row_iter_mut().enumerate().skip(1) { + for element in row.raw_slice_mut().iter_mut().take(i - 1) { + *element = T::zero(); + } + } +} impl Matrix { /// Returns H, where H is the upper hessenberg form. @@ -11,6 +173,11 @@ impl Matrix { /// If the transformation matrix is also required, you should /// use `upper_hess_decomp`. /// + /// Note: This function is deprecated and will be removed in a future + /// release. See + /// [HessenbergDecomposition](decomposition/struct.HessenbergDecomposition.html) + /// for its replacement. + /// /// # Examples /// /// ``` @@ -34,6 +201,7 @@ impl Matrix { /// # Failures /// /// - The matrix cannot be reduced to upper hessenberg form. + #[deprecated] pub fn upper_hessenberg(mut self) -> Result, Error> { let n = self.rows; assert!(n == self.cols, @@ -87,6 +255,11 @@ impl Matrix { /// /// Note: The current transform matrix seems broken... /// + /// Note: This function is deprecated and will be removed in a future + /// release. See + /// [HessenbergDecomposition](decomposition/struct.HessenbergDecomposition.html) + /// for its replacement. + /// /// # Examples /// /// ``` @@ -111,6 +284,7 @@ impl Matrix { /// # Failures /// /// - The matrix cannot be reduced to upper hessenberg form. + #[deprecated] pub fn upper_hess_decomp(self) -> Result<(Matrix, Matrix), Error> { let n = self.rows; assert!(n == self.cols, @@ -138,16 +312,23 @@ impl Matrix { } // Now we reduce to upper hessenberg + #[allow(deprecated)] Ok((transform, try!(self.upper_hessenberg()))) } } #[cfg(test)] mod tests { - use matrix::Matrix; + use super::QH; + use super::HessenbergDecomposition; + + use matrix::{BaseMatrix, Matrix}; + use matrix::decomposition::Decomposition; + use testsupport::is_upper_hessenberg; #[test] #[should_panic] + #[allow(deprecated)] fn test_non_square_upper_hessenberg() { let a: Matrix = Matrix::ones(2, 3); @@ -156,9 +337,70 @@ mod tests { #[test] #[should_panic] + #[allow(deprecated)] fn test_non_square_upper_hess_decomp() { let a: Matrix = Matrix::ones(2, 3); let _ = a.upper_hess_decomp(); } + + fn verify_hessenberg(x: Matrix) { + let QH {ref h, ref q } = HessenbergDecomposition::decompose(x.clone()).unpack(); + + assert!(is_upper_hessenberg(h)); + + let identity = Matrix::identity(h.rows()); + + // Orthogonality. For simplicity, we use a fixed absolute tolerance + // and expect elements to be in the vicinity of 1.0 + assert_matrix_eq!(q.transpose() * q, identity, comp = abs, tol = 1e-14); + assert_matrix_eq!(q * q.transpose(), identity, comp = abs, tol = 1e-14); + + // Reconstruction + assert_matrix_eq!(q * h * q.transpose(), x, comp = abs, tol = 1e-14); + } + + #[test] + fn hessenberg_decomposition() { + { + let x: Matrix = matrix![]; + verify_hessenberg(x); + } + + { + let x = matrix![3.0]; + verify_hessenberg(x); + } + + { + let x = matrix![3.0, 2.0; + 5.0, 1.0]; + verify_hessenberg(x); + } + + { + let x = matrix![3.0, 2.0, -5.0; + 5.0, 1.0, 2.0; + 3.0, -2.0, 0.0]; + verify_hessenberg(x); + } + + { + let x = matrix![3.0, 2.0, -5.0, 2.0; + 5.0, 1.0, 2.0, -2.0; + 3.0, -2.0, 0.0, 0.0; + 4.0, 3.0, -1.0, 3.0]; + verify_hessenberg(x); + } + + { + let x = matrix![1.0, 3.0, -5.0, 2.0, 2.0; + 0.0, 0.0, 5.0, -1.0, 4.0; + 3.0, -4.0, 0.0, 0.0, 5.0; + 0.0, 3.0, -1.0, 3.0, 1.0; + 2.0, -4.0, 3.0, 2.0, 3.0]; + verify_hessenberg(x); + } + } + } diff --git a/src/matrix/decomposition/householder.rs b/src/matrix/decomposition/householder.rs index 4b26b20..735454c 100644 --- a/src/matrix/decomposition/householder.rs +++ b/src/matrix/decomposition/householder.rs @@ -1,6 +1,7 @@ use matrix::{Matrix, BaseMatrix, BaseMatrixMut, Column, ColumnMut}; use vector::Vector; use utils; +use internal_utils::{gemv, transpose_gemv, ger}; use libnum::Float; @@ -95,13 +96,12 @@ impl HouseholderReflection { /// this function computes the product `HA` and stores the result /// back in `A`. /// - /// The user must provide a buffer of size `A.cols()` which is used + /// The user must provide a buffer of size at least `A.cols()` which is used /// to store intermediate results. pub fn buffered_left_multiply_into(&self, matrix: &mut M, buffer: &mut [T]) where M: BaseMatrixMut { - use internal_utils::{transpose_gemv, ger}; - assert!(buffer.len() == matrix.cols()); + assert!(buffer.len() >= matrix.cols()); // Recall that the Householder reflection is represented by // H = I - τ v vᵀ, @@ -120,7 +120,7 @@ impl HouseholderReflection { // Instead, we will use the provided buffer to hold the result of the // matrix-vector product. let ref v = self.v.data(); - let mut u = buffer; + let mut u = &mut buffer[0 .. matrix.cols()]; // u = A^T v transpose_gemv(matrix, v, u); @@ -129,6 +129,35 @@ impl HouseholderReflection { ger(matrix, - self.tau, v, u); } + /// Right-multiplies the given matrix by this Householder reflection. + /// + /// More precisely, let `H` denote this Householder reflection matrix, + /// and let `A` be a dimensionally compatible matrix. Then + /// this function computes the product `AH` and stores the result + /// back in `A`. + /// + /// The user must provide a buffer of size at least `A.rows()` which is used + /// to store intermediate results. + pub fn buffered_right_multiply_into(&self, matrix: &mut M, buffer: &mut [T]) + where M: BaseMatrixMut + { + assert!(buffer.len() >= matrix.rows()); + + // See `buffered_left_multiply_into`. The implementation here + // is almost exactly the same, except we instead have: + // AH = A - τ (A v) vᵀ = A - τ u vᵀ, + // where u = Av. + + let ref v = self.v.data(); + let mut u = &mut buffer[0..matrix.rows()]; + + // u = A v + gemv(matrix, v, u); + + // A <- A - τ v uᵀ + ger(matrix, - self.tau, u, v); + } + pub fn as_vector(&self) -> &Vector { &self.v } @@ -172,12 +201,18 @@ impl HouseholderReflection { /// Q = Q_1 * Q_2 * ... * Q_p /// ``` /// -/// as explained in the documentation for -/// [HouseholderQr](struct.HouseholderQr.html). +/// where `p` is the number of transformations. See the +/// documentation for [HouseholderQr](struct.HouseholderQr.html) +/// and [HessenbergDecomposition](struct.HessenbergDecomposition.html) +/// for more details. #[derive(Debug, Clone)] pub struct HouseholderComposition<'a, T> where T: 'a { storage: &'a Matrix, - tau: &'a [T] + tau: &'a [T], + + // The subdiagonal index determines the subdiagonal to pull Householder + // reflectors from. This is 0 for QR and 1 for Hessenberg. + subdiagonal: usize } /// Instantiates a HouseholderComposition with the given @@ -188,12 +223,13 @@ pub struct HouseholderComposition<'a, T> where T: 'a { /// a HouseholderComposition by themselves, which is desirable /// because we want to have the freedom to change details /// of the internal representation if necessary. -pub fn create_composition<'a, T>(storage: &'a Matrix, tau: &'a [T]) +pub fn create_composition<'a, T>(storage: &'a Matrix, tau: &'a [T], subdiagonal: usize) -> HouseholderComposition<'a, T> { HouseholderComposition { storage: storage, - tau: tau + tau: tau, + subdiagonal: subdiagonal } } @@ -209,15 +245,17 @@ impl<'a, T> HouseholderComposition<'a, T> where T: Float { let n = self.storage.cols(); let p = min(m, n); let q = matrix.cols(); + let s = self.subdiagonal; assert!(matrix.rows() == m, "Matrix does not have compatible dimensions."); let mut house_buffer = Vec::with_capacity(m); let mut multiply_buffer = vec![T::zero(); q]; - for j in (0 .. p).rev() { - house_buffer.resize(m - j, T::zero()); - let storage_block = self.storage.sub_slice([j, j], m - j, n - j); - let mut matrix_block = matrix.sub_slice_mut([j, 0], m - j, q); + for j in (0 .. p.saturating_sub(s)).rev() { + house_buffer.resize(m - j - s, T::zero()); + + let storage_block = self.storage.sub_slice([j + s, j], m - j - s, n - j); + let mut matrix_block = matrix.sub_slice_mut([j + s, 0], m - j - s, q); let house = load_house_from_col(&storage_block.col(0), self.tau[j], house_buffer); house.buffered_left_multiply_into(&mut matrix_block, @@ -237,6 +275,7 @@ impl<'a, T> HouseholderComposition<'a, T> where T: Float { let m = self.storage.rows(); let n = self.storage.cols(); let p = min(m, n); + let s = self.subdiagonal; assert!(k <= self.storage.rows(), "k cannot exceed m, the number of rows of Q"); @@ -252,12 +291,12 @@ impl<'a, T> HouseholderComposition<'a, T> where T: Float { // to reduce the number of operations // (note the size of the "q_k_block") let mut buffer = Vec::with_capacity(m); - let mut multiply_buffer = Vec::with_capacity(k); - for j in (0 .. min(p, k)).rev() { - buffer.resize(m - j, T::zero()); - multiply_buffer.resize(k - j, T::zero()); - let storage_block = self.storage.sub_slice([j, j], m - j, n - j); - let mut q_k_block = q_k.sub_slice_mut([j, j], m - j, k - j); + let mut multiply_buffer = vec![T::zero(); k]; + for j in (0 .. min(p.saturating_sub(s), k)).rev() { + buffer.resize(m - j - s, T::zero()); + + let storage_block = self.storage.sub_slice([j + s, j], m - j - s, n - j); + let mut q_k_block = q_k.sub_slice_mut([j + s, j], m - j - s, k - j); let house = load_house_from_col(&storage_block.col(0), self.tau[j], buffer); house.buffered_left_multiply_into(&mut q_k_block, @@ -377,6 +416,29 @@ mod tests { assert_matrix_eq!(x, expected, comp = abs, tol = 1e-3); } + #[test] + fn householder_reflection_right_multiply() { + let mut x = matrix![ 0.0, 1.0, 2.0, 3.0; + 4.0, 5.0, 6.0, 7.0; + 8.0, 9.0, 10.0, 11.0; + 12.0, 13.0, 14.0, 15.0 ]; + + let h = HouseholderReflection { + tau: 1.0 / 15.0, + v: vector![1.0, 2.0, 3.0, 4.0] + }; + + let mut buffer = vec![0.0; 4]; + + h.buffered_right_multiply_into(&mut x, &mut buffer); + + let expected = matrix![ -1.3333, -1.6667, - 2.0000, - 2.3333; + 0.0000, -3.0000, - 6.0000, - 9.0000; + 1.3333, -4.3333, -10.0000, -15.6667; + 2.6667, -5.6667, -14.0000, -22.3333 ]; + assert_matrix_eq!(x, expected, comp = abs, tol = 1e-3); + } + #[test] fn householder_composition_left_multiply() { let storage = matrix![ 5.0, 3.0, 2.0; @@ -391,7 +453,7 @@ mod tests { // let q = matrix![7.0/9.0, -28.0/45.0, 4.0/45.0; // -4.0/9.0, - 4.0/ 9.0, 7.0/ 9.0; // 4.0/9.0, 29.0/45.0, 28.0/45.0]; - let composition = create_composition(&storage, &tau); + let composition = create_composition(&storage, &tau, 0); { // Square @@ -437,7 +499,7 @@ mod tests { 2.0, 1.0, 3.0; -2.0, 3.0, -2.0]; let tau = vec![2.0/9.0, 1.0 / 5.0, 2.0]; - let composition = create_composition(&storage, &tau); + let composition = create_composition(&storage, &tau, 0); // This corresponds to the following `Q` matrix let q = matrix![7.0/9.0, -28.0/45.0, 4.0/45.0; diff --git a/src/matrix/decomposition/mod.rs b/src/matrix/decomposition/mod.rs index d2f9d93..e737c4c 100644 --- a/src/matrix/decomposition/mod.rs +++ b/src/matrix/decomposition/mod.rs @@ -110,6 +110,12 @@ //! //! //! +//! +//! HessenbergDecomposition +//! Square +//! +//! +//! //! //! @@ -147,6 +153,7 @@ pub use self::householder::HouseholderComposition; pub use self::lu::{PartialPivLu, LUP, FullPivLu, LUPQ}; pub use self::cholesky::Cholesky; pub use self::qr::{HouseholderQr, QR, ThinQR}; +pub use self::hessenberg::{HessenbergDecomposition, QH}; use libnum::{Float}; diff --git a/src/matrix/decomposition/qr.rs b/src/matrix/decomposition/qr.rs index 31dd9f6..412c1d6 100644 --- a/src/matrix/decomposition/qr.rs +++ b/src/matrix/decomposition/qr.rs @@ -206,7 +206,7 @@ impl HouseholderQr where T: Float { // gets shorter for each iteration, so we truncate the buffer // to the appropriate length. buffer.truncate(m - j); - multiply_buffer.truncate(bottom_right.cols()); + bottom_right.col(0).clone_into_slice(&mut buffer); let house = HouseholderReflection::compute(Vector::new(buffer)); @@ -226,7 +226,7 @@ impl HouseholderQr where T: Float { /// [HouseholderComposition](struct.HouseholderComposition.html) /// operator. pub fn q(&self) -> HouseholderComposition { - householder::create_composition(&self.qr, &self.tau) + householder::create_composition(&self.qr, &self.tau, 0) } /// Computes the *thin* (or reduced) QR decomposition. @@ -264,7 +264,7 @@ impl HouseholderQr where T: Float { r1: qr.r } } else { - let composition = householder::create_composition(&self.qr, &self.tau); + let composition = householder::create_composition(&self.qr, &self.tau, 0); let q1 = composition.first_k_columns(n); let r1 = extract_r1(&self.qr); ThinQR { @@ -291,7 +291,7 @@ impl Decomposition for HouseholderQr { fn assemble_q(qr: &Matrix, tau: &Vec) -> Matrix { let m = qr.rows(); - let q_operator = householder::create_composition(qr, tau); + let q_operator = householder::create_composition(qr, tau, 0); q_operator.first_k_columns(m) } diff --git a/src/testsupport/constraints.rs b/src/testsupport/constraints.rs index 4078645..3ca58dd 100644 --- a/src/testsupport/constraints.rs +++ b/src/testsupport/constraints.rs @@ -30,10 +30,28 @@ pub fn is_upper_triangular(m: &M) -> bool .all(|x| x == &T::zero())) } +/// Returns true if the matrix is upper Hessenberg, otherwise false. +/// Note that if the matrix is not square, it cannot be +/// upper Hessenberg, and so this function returns false. +pub fn is_upper_hessenberg(m: &M) -> bool + where T: Zero + PartialEq, M: BaseMatrix +{ + if m.rows() == m.cols() { + m.row_iter() + .enumerate() + .skip(1) + .all(|(i, row)| row.iter().take(i - 1) + .all(|x| x == &T::zero())) + } else { + false + } +} + #[cfg(test)] mod tests { use super::is_lower_triangular; use super::is_upper_triangular; + use super::is_upper_hessenberg; use matrix::Matrix; #[test] @@ -122,18 +140,73 @@ mod tests { } } + #[test] + fn is_upper_hessenberg_valid_input() { + { + let x: Matrix = matrix![]; + assert!(is_upper_hessenberg(&x)); + } + + { + let x = matrix![3]; + assert!(is_upper_hessenberg(&x)); + } + + { + let x = matrix![3, 2; + 4, 0]; + assert!(is_upper_hessenberg(&x)); + } + + { + let x = matrix![3, 2, 0; + 4, 0, 2; + 0, 1, 3]; + assert!(is_upper_hessenberg(&x)); + } + + { + let x = matrix![3, 2, 0, 2; + 4, 0, 2, 1; + 0, 1, 3, 4; + 0, 0, 2, 1]; + assert!(is_upper_hessenberg(&x)); + } + } + + #[test] + fn is_upper_hessenberg_invalid_input() { + { + let x = matrix![3, 2, 0; + 4, 0, 2; + 3, 1, 3]; + assert!(!is_upper_hessenberg(&x)); + } + + { + let x = matrix![3, 2, 0, 5; + 4, 0, 2, 2; + 3, 1, 3, 1; + 0, 3, 2, 0]; + assert!(!is_upper_hessenberg(&x)); + } + } + quickcheck! { fn property_zero_is_lower_triangular(m: usize, n: usize) -> bool { let x = Matrix::::zeros(m, n); is_lower_triangular(&x) } - } - quickcheck! { fn property_zero_is_upper_triangular(m: usize, n: usize) -> bool { let x = Matrix::::zeros(m, n); is_upper_triangular(&x) } + + fn property_zero_is_upper_hessenberg(n: usize) -> bool { + let x = Matrix::::zeros(n, n); + is_upper_hessenberg(&x) + } } } diff --git a/src/testsupport/mod.rs b/src/testsupport/mod.rs index 3277d1f..9356a0e 100644 --- a/src/testsupport/mod.rs +++ b/src/testsupport/mod.rs @@ -2,4 +2,5 @@ //! tests involving data structures provided by rulinalg. mod constraints; -pub use self::constraints::{is_lower_triangular, is_upper_triangular}; +pub use self::constraints::{is_lower_triangular, is_upper_triangular, + is_upper_hessenberg};