From 4e83546dc3b80bd899d6cc29acf036e7696baa50 Mon Sep 17 00:00:00 2001 From: Richard Janis Goldschmidt Date: Mon, 28 Jan 2019 15:31:45 +0100 Subject: [PATCH 1/3] Update to 2018 and dependencies --- Cargo.toml | 15 +++++++-------- src/cgs.rs | 12 ------------ src/cgs2.rs | 11 ----------- src/lib.rs | 17 +---------------- src/mgs.rs | 13 +------------ src/test_macros.rs | 12 +++++++----- 6 files changed, 16 insertions(+), 64 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index c687e28..de3b5b5 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,7 @@ [package] name = "gramschmidt" version = "0.4.1" +edition = "2018" authors = ["Richard Janis Goldschmidt "] license = "MIT" description = "Classical, Modified, Reorthogonalized Gram Schmidt Orthogonalization and QR decompostion" @@ -10,16 +11,14 @@ repository = "https://github.com/SuperFluffy/gramschmidt-rs" readme = "README.md" [dependencies] -ndarray = "0.11" +ndarray = "0.12.1" [dependencies.cblas] -version = "0.1.5" +version = "0.2.0" default-features = false [dev-dependencies] -lazy_static = "1.0" -ndarray-rand = "0.7" -rand = "0.4" - -[dev-dependencies.openblas-src] -version = "0.5.6" +lazy_static = "1.2.0" +ndarray-rand = "0.9.0" +rand = "0.6.5" +openblas-src = "0.7.0" diff --git a/src/cgs.rs b/src/cgs.rs index 8abe136..a93c1d7 100644 --- a/src/cgs.rs +++ b/src/cgs.rs @@ -22,11 +22,6 @@ impl ClassicalGramSchmidt { /// # Example /// /// ``` - /// extern crate gramschmidt; - /// extern crate ndarray; - /// extern crate ndarray_rand; - /// extern crate rand; - /// /// use gramschmidt::ClassicalGramSchmidt; /// use ndarray::Array2; /// use ndarray_rand::RandomExt; @@ -67,8 +62,6 @@ impl ClassicalGramSchmidt { /// # Example /// /// ``` - /// extern crate gramschmidt; - /// /// use gramschmidt::ClassicalGramSchmidt; /// /// # fn main() { @@ -101,11 +94,6 @@ impl ClassicalGramSchmidt { /// previously configured. Panics otherwise. /// /// ``` - /// extern crate gramschmidt; - /// extern crate ndarray; - /// extern crate ndarray_rand; - /// extern crate rand; - /// /// use gramschmidt::ClassicalGramSchmidt; /// use ndarray::Array2; /// use ndarray_rand::RandomExt; diff --git a/src/cgs2.rs b/src/cgs2.rs index 387fcfb..e1a8777 100644 --- a/src/cgs2.rs +++ b/src/cgs2.rs @@ -26,11 +26,6 @@ impl ReorthogonalizedGramSchmidt { /// # Example /// /// ``` - /// extern crate gramschmidt; - /// extern crate ndarray; - /// extern crate ndarray_rand; - /// extern crate rand; - /// /// use gramschmidt::ReorthogonalizedGramSchmidt; /// use ndarray::Array2; /// use ndarray_rand::RandomExt; @@ -72,8 +67,6 @@ impl ReorthogonalizedGramSchmidt { /// # Example /// /// ``` - /// extern crate gramschmidt; - /// /// use gramschmidt::ReorthogonalizedGramSchmidt; /// /// # fn main() { @@ -109,11 +102,7 @@ impl ReorthogonalizedGramSchmidt { /// previously configured. Panics otherwise. /// /// ``` - /// extern crate gramschmidt; - /// extern crate ndarray; - /// extern crate ndarray_rand; /// extern crate openblas_src; - /// extern crate rand; /// /// use gramschmidt::ReorthogonalizedGramSchmidt; /// use ndarray::Array2; diff --git a/src/lib.rs b/src/lib.rs index c89d410..697ce06 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -1,21 +1,6 @@ -extern crate cblas; -extern crate ndarray; - -#[cfg(test)] -#[macro_use] -extern crate lazy_static; - -#[cfg(test)] -extern crate ndarray_rand; - -#[cfg(test)] -extern crate openblas_src; - -#[cfg(test)] -extern crate rand; - #[macro_use] mod test_macros; + mod cgs; mod cgs2; mod mgs; diff --git a/src/mgs.rs b/src/mgs.rs index 7c4697c..448095e 100644 --- a/src/mgs.rs +++ b/src/mgs.rs @@ -1,4 +1,5 @@ use cblas; + use ndarray::{ Data, IntoDimension, @@ -22,11 +23,6 @@ impl ModifiedGramSchmidt { /// # Example /// /// ``` - /// extern crate gramschmidt; - /// extern crate ndarray; - /// extern crate ndarray_rand; - /// extern crate rand; - /// /// use gramschmidt::ModifiedGramSchmidt; /// use ndarray::Array2; /// use ndarray_rand::RandomExt; @@ -64,8 +60,6 @@ impl ModifiedGramSchmidt { /// # Example /// /// ``` - /// extern crate gramschmidt; - /// /// use gramschmidt::ModifiedGramSchmidt; /// /// # fn main() { @@ -97,11 +91,6 @@ impl ModifiedGramSchmidt { /// previously configured. Panics otherwise. /// /// ``` - /// extern crate gramschmidt; - /// extern crate ndarray; - /// extern crate ndarray_rand; - /// extern crate rand; - /// /// use gramschmidt::ModifiedGramSchmidt; /// use ndarray::Array2; /// use ndarray_rand::RandomExt; diff --git a/src/test_macros.rs b/src/test_macros.rs index bb855cc..828be46 100644 --- a/src/test_macros.rs +++ b/src/test_macros.rs @@ -2,9 +2,11 @@ macro_rules! generate_tests { ($method:ident, $tolerance:expr) => { #[cfg(test)] mod tests { + extern crate openblas_src; + + use lazy_static::lazy_static; use ndarray::prelude::*; use super::*; - use utils; lazy_static!( static ref UNITY: Array2 = arr2( @@ -80,7 +82,7 @@ macro_rules! generate_tests { fn small_orthogonal() { let mut method = $method::from_matrix(&*SMALL); method.compute(&*SMALL); - assert!(utils::orthogonal(method.q(),$tolerance)); + assert!(crate::utils::orthogonal(method.q(),$tolerance)); } #[test] @@ -94,7 +96,7 @@ macro_rules! generate_tests { fn large_orthogonal() { let mut method = $method::from_matrix(&*LARGE); method.compute(&*LARGE); - assert!(utils::orthogonal(method.q(),$tolerance)); + assert!(crate::utils::orthogonal(method.q(),$tolerance)); } #[test] @@ -116,7 +118,7 @@ macro_rules! generate_tests { fn f_order_small_orthogonal() { let mut method = $method::from_matrix(&*F_SMALL); method.compute(&*F_SMALL); - assert!(utils::orthogonal(method.q(),$tolerance)); + assert!(crate::utils::orthogonal(method.q(),$tolerance)); } #[test] @@ -130,7 +132,7 @@ macro_rules! generate_tests { fn f_order_large_orthogonal() { let mut method = $method::from_matrix(&*F_LARGE); method.compute(&*F_LARGE); - assert!(utils::orthogonal(method.q(),$tolerance)); + assert!(crate::utils::orthogonal(method.q(),$tolerance)); } #[test] From bbba1daaab0e4bd8b415cef280c17db2f6b0c191 Mon Sep 17 00:00:00 2001 From: Richard Janis Goldschmidt Date: Tue, 12 Mar 2019 12:24:42 +0100 Subject: [PATCH 2/3] Refactor and update to edition 2018 --- Cargo.toml | 9 +- benches/matrix_sizes.rs | 45 ++--- src/cgs.rs | 340 ++++++++++++++------------------- src/cgs2.rs | 411 +++++++++++++++++----------------------- src/lib.rs | 161 +++++++++++++++- src/mgs.rs | 154 ++++++--------- src/test_macros.rs | 40 ++-- src/utils.rs | 29 +++ 8 files changed, 618 insertions(+), 571 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index de3b5b5..ee50663 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -11,14 +11,11 @@ repository = "https://github.com/SuperFluffy/gramschmidt-rs" readme = "README.md" [dependencies] +cblas = "0.2.0" ndarray = "0.12.1" -[dependencies.cblas] -version = "0.2.0" -default-features = false - [dev-dependencies] -lazy_static = "1.2.0" +lazy_static = "1.3.0" ndarray-rand = "0.9.0" -rand = "0.6.5" openblas-src = "0.7.0" +rand = "0.6.5" diff --git a/benches/matrix_sizes.rs b/benches/matrix_sizes.rs index 3a14784..98a1ca8 100644 --- a/benches/matrix_sizes.rs +++ b/benches/matrix_sizes.rs @@ -2,16 +2,17 @@ #![allow(non_snake_case)] -extern crate gram_schmidt; -extern crate ndarray; +// extern crate gram_schmidt; +// extern crate ndarray; extern crate openblas_src; extern crate test; // Built-in crate for benchmarking. -use gram_schmidt::{ - ClassicalGramSchmidt, - ModifiedGramSchmidt, - ReorthogonalizedGramSchmidt, +use gramschmidt::{ + GramSchmidt, + Classical, + Modified, + Reorthogonalized, }; use ndarray::prelude::*; @@ -51,32 +52,32 @@ use ndarray::prelude::*; // } macro_rules! create_bench { - (c $n:expr, $name:ident, $method_constructor:path) => { + (c $n:expr, $name:ident, $method:ty) => { #[bench] fn $name(bench: &mut test::Bencher) { let n = $n; let matrix = Array2::eye(n); - let mut method = $method_constructor(&matrix); + let mut method = <$method>::from_matrix(&matrix).unwrap(); let method = test::black_box(&mut method); bench.iter(|| { - method.compute(&matrix); + method.compute(&matrix).unwrap(); }); } }; - (f $n:expr, $name:ident, $method_constructor:path) => { + (f $n:expr, $name:ident, $method:ty) => { #[bench] fn $name(bench: &mut test::Bencher) { let n = $n; let matrix = Array2::eye(n).reversed_axes(); - let mut method = $method_constructor(&matrix); + let mut method = <$method>::from_matrix(&matrix).unwrap(); let method = test::black_box(&mut method); bench.iter(|| { - method.compute(&matrix); + method.compute(&matrix).unwrap(); }); } }; @@ -84,28 +85,28 @@ macro_rules! create_bench { macro_rules! bench_sizes { (c $n:expr, $name_cgs:ident, $name_mgs: ident, $name_cgs2: ident) => { - create_bench!(c $n, $name_cgs, ClassicalGramSchmidt::from_matrix); - create_bench!(c $n, $name_cgs2, ReorthogonalizedGramSchmidt::from_matrix); - create_bench!(c $n, $name_mgs, ModifiedGramSchmidt::from_matrix); + create_bench!(c $n, $name_cgs, Classical); + create_bench!(c $n, $name_cgs2, Reorthogonalized); + create_bench!(c $n, $name_mgs, Modified); }; (f $n:expr, $name_cgs:ident, $name_mgs: ident, $name_cgs2: ident) => { - create_bench!(f $n, $name_cgs, ClassicalGramSchmidt::from_matrix); - create_bench!(f $n, $name_cgs2, ReorthogonalizedGramSchmidt::from_matrix); - create_bench!(f $n, $name_mgs, ModifiedGramSchmidt::from_matrix); + create_bench!(f $n, $name_cgs, Classical); + create_bench!(f $n, $name_cgs2, Reorthogonalized); + create_bench!(f $n, $name_mgs, Modified); }; } bench_sizes!(c 256, c_cgs__256, c_mgs__256, c_cgs2__256); bench_sizes!(c 512, c_cgs__512, c_mgs__512, c_cgs2__512); -bench_sizes!(c 768, c_cgs__768, c_mgs__768, c_cgs2__768); -bench_sizes!(c 1024, c_cgs_1024, c_mgs_1024, c_cgs2_1024); +// bench_sizes!(c 768, c_cgs__768, c_mgs__768, c_cgs2__768); +// bench_sizes!(c 1024, c_cgs_1024, c_mgs_1024, c_cgs2_1024); // bench_sizes!(c 1536, c_cgs_1536, c_mgs_1536, c_cgs2_1536); // bench_sizes!(c 2048, c_cgs_2048, c_mgs_2048, c_cgs2_2048); bench_sizes!(f 256, f_cgs__256, f_mgs__256, f_cgs2__256); bench_sizes!(f 512, f_cgs__512, f_mgs__512, f_cgs2__512); -bench_sizes!(f 768, f_cgs__768, f_mgs__768, f_cgs2__768); -bench_sizes!(f 1024, f_cgs_1024, f_mgs_1024, f_cgs2_1024); +// bench_sizes!(f 768, f_cgs__768, f_mgs__768, f_cgs2__768); +// bench_sizes!(f 1024, f_cgs_1024, f_mgs_1024, f_cgs2_1024); // bench_sizes!(f 1536, f_cgs_1536, f_mgs_1536, f_cgs2_1536); // bench_sizes!(f 2048, f_cgs_2048, f_mgs_2048, f_cgs2_2048); diff --git a/src/cgs.rs b/src/cgs.rs index a93c1d7..d7ffc4a 100644 --- a/src/cgs.rs +++ b/src/cgs.rs @@ -1,248 +1,202 @@ use cblas; use ndarray::{ Data, - IntoDimension, + ShapeBuilder, }; use ndarray::prelude::*; use std::slice; +use crate::{ + Error, + GramSchmidt, + Result, + utils::{ + as_slice_with_layout, + get_layout, + }, +}; + +/// A classical Gram Schmidt factorization. See the [Gram Schmidt Wikipedia entry] for more information. +/// +/// Use this struct via the [`GramSchmidt` trait]. +/// +/// [Gram Schmidt Wikipedia entry]: https://en.wikipedia.org/wiki/Gram-Schmidt_process +/// [`GramSchmidt` trait]: GramSchmidt #[derive(Clone, Debug)] -pub struct ClassicalGramSchmidt { +pub struct Classical { q: Array2, r: Array2, - memory_order: cblas::Layout, + memory_layout: cblas::Layout, } -impl ClassicalGramSchmidt { - /// Reserves the memory for a QR decomposition via a classical Gram Schmidt orthogonalization - /// using the dimensions of a sample matrix. - /// - /// The resulting object can be used to orthogonalize matrices of the same dimensions. - /// - /// # Example - /// - /// ``` - /// use gramschmidt::ClassicalGramSchmidt; - /// use ndarray::Array2; - /// use ndarray_rand::RandomExt; - /// use rand::distributions::Normal; - /// - /// # fn main() { - /// - /// let matrix = Array2::random((10,10), Normal::new(0.0, 1.0)); - /// let mut cgs = ClassicalGramSchmidt::from_matrix(&matrix); - /// - /// # } - /// ``` - pub fn from_matrix(a: &ArrayBase) -> Self - where S: Data +impl GramSchmidt for Classical { + fn from_shape(shape: T) -> Result + where T: ShapeBuilder, { - let (is_fortran, memory_order) = if let Some(_) = a.as_slice() { - (false, cblas::Layout::RowMajor) - } else if let Some(_) = a.as_slice_memory_order() { - (true, cblas::Layout::ColumnMajor) - } else { - panic!("Array not contiguous!") + // Unfortunately we cannot check the shape itself to see if it's + // in ColumnMajor or RowMajor layout. So we need to first construct + // an array and then check that. + let shape = shape.into_shape(); + let q = Array2::zeros(shape); + let memory_layout = match get_layout(&q) { + Some(layout) => layout, + None => Err(Error::NonContiguous)?, }; - - let array_shape = a.raw_dim().set_f(is_fortran); - - ClassicalGramSchmidt { - q: Array2::zeros(array_shape), - r: Array2::zeros(array_shape), - memory_order, - } + let r = q.clone(); + Ok(Self { + q, + r, + memory_layout, + }) } - /// Reserves the memory for a QR decomposition via a classical Gram Schmidt orthogonalization - /// using a shape. - /// - /// The resulting object can be used to orthogonalize matrices of the same dimensions. - /// - /// # Example - /// - /// ``` - /// use gramschmidt::ClassicalGramSchmidt; - /// - /// # fn main() { - /// let fortran_order = false; - /// let mut cgs = ClassicalGramSchmidt::from_shape((10,10), fortran_order); - /// - /// # } - /// ``` - pub fn from_shape(shape: T, fortran_order: bool) -> Self - where T: IntoDimension, + fn compute(&mut self, a: &ArrayBase) -> Result<()> + where S: Data { - let dimension = shape.into_dimension(); - let array_shape = dimension.set_f(fortran_order); - let memory_order = if fortran_order { - cblas::Layout::ColumnMajor - } else { - cblas::Layout::RowMajor - }; - ClassicalGramSchmidt { - q: Array2::zeros(array_shape), - r: Array2::zeros(array_shape), - memory_order, - } - } + use cblas::Layout::*; + use Error::*; - /// Computes a QR decomposition using the classical Gram Schmidt orthonormalization of the - /// matrix `a`. - /// - /// The input matrix `a` has to have exactly the same dimension and memory layout as was - /// previously configured. Panics otherwise. - /// - /// ``` - /// use gramschmidt::ClassicalGramSchmidt; - /// use ndarray::Array2; - /// use ndarray_rand::RandomExt; - /// use rand::distributions::Normal; - /// - /// # fn main() { - /// - /// let matrix = Array2::random((10,10), Normal::new(0.0, 1.0)); - /// let mut cgs = ClassicalGramSchmidt::from_matrix(&matrix); - /// - /// # } - /// ``` - pub fn compute(&mut self, a: &ArrayBase) - where S: Data, - { assert_eq!(a.shape(), self.q.shape()); - // assert_eq!(a.strides(), self.q.strides()); let (n_rows, n_cols) = self.q.dim(); - let (a_column_inc, a_layout) = if let Some(_) = a.as_slice() { - (n_cols as i32, cblas::Layout::RowMajor) - } else if let Some(_) = a.as_slice_memory_order() { - (1, cblas::Layout::ColumnMajor) - } else { - panic!("Array not contiguous!") + let a_slice = match (self.memory_layout, as_slice_with_layout(&a)) { + (a, Some((_, b))) if a != b => Err(IncompatibleLayouts)?, + (_, Some((a_slice, _))) => a_slice, + (_, None) => Err(NonContiguous)?, + }; + + // leading_dim: the number of elements in the leading dimension + // next_elem: how many elements to jump to get to the next element in a column + // next_col: how many elements in the array to jump to get to the next column + let (leading_dim, next_elem, next_col) = match self.memory_layout { + ColumnMajor => (n_rows as i32, 1, n_rows), + RowMajor => (n_cols as i32, n_cols as i32, 1), }; for i in 0..n_cols { self.q.column_mut(i).assign(&a.column(i)); - let a_column = match a_layout { - cblas::Layout::RowMajor => &a.as_slice_memory_order().unwrap()[i..], - cblas::Layout::ColumnMajor => &a.as_slice_memory_order().unwrap()[n_rows * i..], + // The unsafe blocks below are because we need several overlapping slices into the + // q matrix. The mutable `q_column` is the i-th orthogonal vector which is currently + // being sought. The immutable `q` represents the already found vectors (vectors 0 + // to i-1). Appropriately choosing the offset, lda, and increment makes sure that + // blas does not access the same memory location. + // + // This is only necessary for row major formats. In column major formats there + // won't be overlaps. + + let q_len = self.q.len(); + let q_ptr = self.q.as_mut_ptr(); + let q_matrix = unsafe { + slice::from_raw_parts(q_ptr, q_len) + }; + + let q_column = match self.memory_layout { + ColumnMajor => { + let offset = n_rows * i; + unsafe { + slice::from_raw_parts_mut(q_ptr.offset(offset as isize), q_len - offset) + } + }, + + RowMajor => { + let offset = i as isize; + unsafe { + slice::from_raw_parts_mut(q_ptr.offset(offset), q_len - i) + } + }, + }; if i > 0 { - let (lda, r_column, r_inc) = match self.memory_order { - cblas::Layout::RowMajor => { - let r_column = &mut self.r.as_slice_memory_order_mut().unwrap()[i..]; - (n_cols as i32, r_column, n_cols as i32) - }, - - cblas::Layout::ColumnMajor => { - let r_column = &mut self.r.as_slice_memory_order_mut().unwrap()[n_rows * i..]; - (n_rows as i32, r_column, 1) - }, - }; - - // The unsafe blocks below are because we need several overlapping slices into the - // q matrix. The mutable `q_column` is the i-th orthogonal vector which is currently - // being sought. The immutable `q` represents the already found vectors (vectors 0 - // to i-1). Appropriately choosing the offset, lda, and increment makes sure that - // blas does not access the same memory location. - // - // This is only necessary for row major formats. In column major formats there - // won't be overlaps. - let len = self.q.len(); - let q_ptr = self.q.as_mut_ptr(); - let q = unsafe { - slice::from_raw_parts(q_ptr, len) - }; - let (q_column, q_inc) = match self.memory_order { - cblas::Layout::RowMajor => { - let offset = i as isize; - let q_column = unsafe { - slice::from_raw_parts_mut(q_ptr.offset(offset), len - i) - }; - (q_column, n_cols as i32) - }, - - cblas::Layout::ColumnMajor => { - let offset = n_rows * i; - let q_column = unsafe { - slice::from_raw_parts_mut(q_ptr.offset(offset as isize), len - offset) - }; - (q_column, 1) - }, - }; + let a_column = &a_slice[next_col * i..]; + // NOTE: This unwrap is save, because we have made sure at creation that r_slice is + // contiguous. + // + // NOTE: Unlike a_slice above which is defined outside the loop, we are mutating r at the + // end of the loop, which invalidates the mutable borrow. We thus have to pull the + // slice definition into the loop. + let r_slice = self.r.as_slice_memory_order_mut().unwrap(); + let r_column = &mut r_slice[next_col * i..]; + + // Calculate the product R_(i) = Q^T·A_(i), where A_(i) is the i-th column of the matrix A, + // and R_(i) is the i-th column of matrix R. unsafe { cblas::dgemv( - self.memory_order, // layout - cblas::Transpose::Ordinary, // transa - n_rows as i32, // m - i as i32, // n - 1.0, // alpha - q, // a - lda, // lda - a_column, // x - a_column_inc, // incx - 0.0, // beta - r_column, //y - r_inc // incy + self.memory_layout, + cblas::Transpose::Ordinary, + n_rows as i32, + i as i32, + 1.0, + q_matrix, + leading_dim, + a_column, + next_elem, + 0.0, + r_column, + next_elem, ); + // Calculate Q_(i) = A_(i) - Q · R_(i) = A_(i) - Q · (Q^T · A_(i)), where + // Q · (Q^T ·A_(i)) is the projection of the i-th column of A onto the already + // orthonormalized basis vectors Q_{0..i}. cblas::dgemv( - self.memory_order, // layout - cblas::Transpose::None, // transa - n_rows as i32, // m - i as i32, // n - -1.0, // alpha - q, // a - lda, // lda - r_column, // y - r_inc, // incx - 1.0, // beta - q_column, // y - q_inc, // incy + self.memory_layout, + cblas::Transpose::None, + n_rows as i32, + i as i32, + -1.0, + q_matrix, + leading_dim, + r_column, + next_elem, + 1.0, + q_column, + next_elem, ); } }; - let norm = { - let len = self.q.len(); - let q_ptr = self.q.as_mut_ptr(); - unsafe { - let (q_column, q_inc) = match self.memory_order { - cblas::Layout::RowMajor => { - let offset = i as isize; - let q_column = slice::from_raw_parts_mut(q_ptr.offset(offset), len - i); - (q_column, n_cols as i32) - }, - - cblas::Layout::ColumnMajor => { - let offset = n_rows * i; - let q_column = slice::from_raw_parts_mut(q_ptr.offset(offset as isize), len - offset); - (q_column, 1) - }, - }; - cblas::dnrm2(n_rows as i32, q_column, q_inc) - } + let norm = unsafe { + cblas::dnrm2(n_rows as i32, q_column, next_elem) }; let mut v = self.q.column_mut(i); v /= norm; self.r[(i,i)] = a.column(i).dot(&v); } + + Ok(()) } - /// Return a reference to the matrix q. - pub fn q(&self) -> &Array2 { + fn q(&self) -> &Array2 { &self.q } - /// Return a reference to the matrix q. - pub fn r(&self) -> &Array2 { + fn r(&self) -> &Array2 { &self.r } } -generate_tests!(ClassicalGramSchmidt, 1e-12); +/// Convenience function that calculates a [Classical Gram Schmidt] QR factorization, returning a +/// tuple `(Q,R)`. +/// +/// If you want to repeatedly calculate QR factorizations, then prefer constructing a [`Classical`] +/// struct and calling its [`GramSchmidt::compute`] method implemented through the [`GramSchmidt`] trait. +/// +/// [Classical Gram Schmidt]: https://en.wikipedia.org/wiki/Gram-Schmidt_process +/// [`Classical`]: Classical +/// [GramSchmidt]: GramSchmidt +/// [`GramSchmidt::compute`]: trait.GramSchmidt.html#tymethod.compute +pub fn cgs(a: &ArrayBase) -> Result<(Array, Array)> + where S: Data +{ + let mut cgs = Classical::from_matrix(a)?; + cgs.compute(a)?; + Ok((cgs.q().clone(), cgs.r().clone())) +} + +#[cfg(test)] +generate_tests!(Classical, 1e-12); diff --git a/src/cgs2.rs b/src/cgs2.rs index e1a8777..a05ddf9 100644 --- a/src/cgs2.rs +++ b/src/cgs2.rs @@ -1,301 +1,246 @@ use cblas; use ndarray::{ Data, - IntoDimension, + ShapeBuilder, }; use ndarray::prelude::*; use std::slice; +use crate::{ + Error, + GramSchmidt, + Result, + utils::{ + as_slice_with_layout, + get_layout, + } +}; + +/// A reorthogonalized Gram Schmidt factorization, also known as `CGS2` in the literature. See +/// [Giraud et al.] for a definition. It performs two successive classical Gram Schmidt procedures, +/// which has a higher performance than modified Gram Schmidt while providing a similar numerical +/// stability. +/// +/// Use this struct via the [`GramSchmidt` trait]. +/// +/// [Giraud et al.]: https://doi.org/10.1007/s00211-005-0615-4 +/// [`GramSchmidt` trait]: GramSchmidt #[derive(Clone, Debug)] -pub struct ReorthogonalizedGramSchmidt { +pub struct Reorthogonalized { q: Array2, r: Array2, work_vector: Array1, - memory_order: cblas::Layout, + memory_layout: cblas::Layout, } -impl ReorthogonalizedGramSchmidt { - /// Reserves the memory for a QR decomposition via a classical, reorthogonalized Gram Schmidt - /// orthogonalization using the dimensions of a sample matrix. This is also known as CGS2 in - /// the literature. See [Giraud et al.] for a definition. - /// - /// The resulting object can be used to orthogonalize matrices of the same dimensions. - /// - /// [Giraud et al.]: https://doi.org/10.1007/s00211-005-0615-4 - /// - /// # Example - /// - /// ``` - /// use gramschmidt::ReorthogonalizedGramSchmidt; - /// use ndarray::Array2; - /// use ndarray_rand::RandomExt; - /// use rand::distributions::Normal; - /// - /// # fn main() { - /// - /// let matrix = Array2::random((10,10), Normal::new(0.0, 1.0)); - /// let mut cgs2 = ReorthogonalizedGramSchmidt::from_matrix(&matrix); - /// - /// # } - /// ``` - pub fn from_matrix(a: &ArrayBase) -> Self - where S: Data +impl GramSchmidt for Reorthogonalized { + fn from_shape(shape: T) -> Result + where T: ShapeBuilder, { - let (is_fortran, memory_order) = if let Some(_) = a.as_slice() { - (false, cblas::Layout::RowMajor) - } else if let Some(_) = a.as_slice_memory_order() { - (true, cblas::Layout::ColumnMajor) - } else { - panic!("Array not contiguous!") + // Unfortunately we cannot check the shape itself to see if it's + // in ColumnMajor or RowMajor layout. So we need to first construct + // an array and then check that. + let shape = shape.into_shape(); + let q = Array2::zeros(shape); + let memory_layout = match get_layout(&q) { + Some(layout) => layout, + None => Err(Error::NonContiguous)?, }; - - let array_shape = a.raw_dim().set_f(is_fortran); - - ReorthogonalizedGramSchmidt { - q: Array2::zeros(array_shape), - r: Array2::zeros(array_shape), - work_vector: Array1::zeros(a.dim().0), - memory_order, - } + let r = q.clone(); + + // Similarly to the layout, we don't have direct access to the array dimensions via + // `Shape`, and thus need to go via `Dim::Pattern` of the already constructed arrays. + let work_vector = Array1::zeros(q.dim().0); + + Ok(Self { + q, + r, + work_vector, + memory_layout, + }) } - /// Reserves the memory for a QR decomposition via a classical, reorthogonalized Gram Schmidt - /// orthogonalization using a shape. - /// - /// The resulting object can be used to orthogonalize matrices of the same dimensions. - /// - /// # Example - /// - /// ``` - /// use gramschmidt::ReorthogonalizedGramSchmidt; - /// - /// # fn main() { - /// let fortran_order = false; - /// let mut cgs2 = ReorthogonalizedGramSchmidt::from_shape((10,10), fortran_order); - /// - /// # } - /// ``` - pub fn from_shape(shape: T, fortran_order: bool) -> Self - where T: IntoDimension, - { - let dimension = shape.into_dimension(); - let rows = dimension.into_pattern().0; - let array_shape = dimension.set_f(fortran_order); - let memory_order = if fortran_order { - cblas::Layout::ColumnMajor - } else { - cblas::Layout::RowMajor - }; - - ReorthogonalizedGramSchmidt { - q: Array2::zeros(array_shape), - r: Array2::zeros(array_shape), - work_vector: Array1::zeros(rows), - memory_order, - } - } - - /// Computes a QR decomposition using the classical, reorthogonalized Gram Schmidt - /// orthonormalization of the matrix `a`. - /// - /// The input matrix `a` has to have exactly the same dimension and memory layout as was - /// previously configured. Panics otherwise. - /// - /// ``` - /// extern crate openblas_src; - /// - /// use gramschmidt::ReorthogonalizedGramSchmidt; - /// use ndarray::Array2; - /// use ndarray_rand::RandomExt; - /// use rand::distributions::Normal; - /// - /// # fn main() { - /// - /// let matrix = Array2::random((10,10), Normal::new(0.0, 1.0)); - /// let mut cgs2 = ReorthogonalizedGramSchmidt::from_matrix(&matrix); - /// cgs2.compute(&matrix); - /// - /// # } - /// ``` - pub fn compute(&mut self, a: &ArrayBase) + fn compute(&mut self, a: &ArrayBase) -> Result<()> where S: Data, { + use cblas::Layout::*; + use Error::*; + assert_eq!(a.shape(), self.q.shape()); - // assert_eq!(a.strides(), self.q.strides()); let (n_rows, n_cols) = self.q.dim(); - let (a_column_inc, a_layout) = if let Some(_) = a.as_slice() { - (n_cols as i32, cblas::Layout::RowMajor) - } else if let Some(_) = a.as_slice_memory_order() { - (1, cblas::Layout::ColumnMajor) - } else { - panic!("Array not contiguous!") + let a_slice = match (self.memory_layout, as_slice_with_layout(&a)) { + (a, Some((_, b))) if a != b => Err(IncompatibleLayouts)?, + (_, Some((a_slice, _))) => a_slice, + (_, None) => Err(NonContiguous)?, + }; + + // leading_dim: the number of elements in the leading dimension + // next_elem: how many elements to jump to get to the next element in a column + // next_col: how many elements in the array to jump to get to the next column + let (leading_dim, next_elem, next_col) = match self.memory_layout { + ColumnMajor => (n_rows as i32, 1, n_rows), + RowMajor => (n_cols as i32, n_cols as i32, 1), }; + for i in 0..n_cols { self.q.column_mut(i).assign(&a.column(i)); - - let a_column = match a_layout { - cblas::Layout::RowMajor => &a.as_slice_memory_order().unwrap()[i..], - cblas::Layout::ColumnMajor => &a.as_slice_memory_order().unwrap()[n_rows * i..], + + let len = self.q.len(); + let q_ptr = self.q.as_mut_ptr(); + let q_matrix = unsafe { + slice::from_raw_parts(q_ptr, len) + }; + + let q_column = match self.memory_layout { + ColumnMajor => { + let offset = n_rows * i; + unsafe { + slice::from_raw_parts_mut(q_ptr.offset(offset as isize), len - offset) + } + }, + + RowMajor => { + let offset = i as isize; + unsafe { + slice::from_raw_parts_mut(q_ptr.offset(offset), len - i) + } + }, + }; if i > 0 { - let (lda, r_column, r_inc) = match self.memory_order { - cblas::Layout::RowMajor => { - let r_column = &mut self.r.as_slice_memory_order_mut().unwrap()[i..]; - (n_cols as i32, r_column, n_cols as i32) - }, - - cblas::Layout::ColumnMajor => { - let r_column = &mut self.r.as_slice_memory_order_mut().unwrap()[n_rows * i..]; - (n_rows as i32, r_column, 1) - }, - }; - - // The unsafe blocks below are because we need several overlapping slices into the - // q matrix. The mutable `q_column` is the i-th orthogonal vector which is currently - // being sought. The immutable `q` represents the already found vectors (vectors 0 - // to i-1). Appropriately choosing the offset, lda, and increment makes sure that - // blas does not access the same memory location. + let a_column = &a_slice[next_col * i..]; + + // NOTE: This unwrap is save, because we have made sure at creation that r_slice is + // contiguous. // - // This is only necessary for row major formats. In column major formats there - // won't be overlaps. - let len = self.q.len(); - let q_ptr = self.q.as_mut_ptr(); - let q = unsafe { - slice::from_raw_parts(q_ptr, len) - }; - let (q_column, q_inc) = match self.memory_order { - cblas::Layout::RowMajor => { - let offset = i as isize; - let q_column = unsafe { - slice::from_raw_parts_mut(q_ptr.offset(offset), len - i) - }; - (q_column, n_cols as i32) - }, - - cblas::Layout::ColumnMajor => { - let offset = n_rows * i; - let q_column = unsafe { - slice::from_raw_parts_mut(q_ptr.offset(offset as isize), len - offset) - }; - (q_column, 1) - }, - }; - - let work_vector = self.work_vector.as_slice_memory_order_mut().unwrap(); + // NOTE: Unlike a_slice above which is defined outside the loop, we are mutating r at the + // end of the loop, which invalidates the mutable borrow. We thus have to pull the + // slice definition into the loop. + let r_slice = self.r.as_slice_memory_order_mut().unwrap(); + let r_column = &mut r_slice[next_col * i..]; + + let work_slice = self.work_vector.as_slice_memory_order_mut().unwrap(); unsafe { + // First orthogonalization + // ======================= cblas::dgemv( - self.memory_order, // layout - cblas::Transpose::Ordinary, // transa - n_rows as i32, // m - i as i32, // n - 1.0, // alpha - q, // a - lda, // lda - a_column, // x - a_column_inc, // incx - 0.0, // beta - r_column, //y - r_inc // incy + self.memory_layout, + cblas::Transpose::Ordinary, + n_rows as i32, + i as i32, + 1.0, + q_matrix, + leading_dim, + a_column, + next_elem, + 0.0, + r_column, + next_elem ); cblas::dgemv( - self.memory_order, // layout - cblas::Transpose::None, // transa - n_rows as i32, // m - i as i32, // n - -1.0, // alpha - q, // a - lda, // lda - r_column, // y - r_inc, // incx - 1.0, // beta - q_column, // y - q_inc, // incy + self.memory_layout, + cblas::Transpose::None, + n_rows as i32, + i as i32, + -1.0, + q_matrix, + leading_dim, + r_column, + next_elem, + 1.0, + q_column, + next_elem, ); + // Second orthogonalization + // ======================= cblas::dgemv( - self.memory_order, // layout - cblas::Transpose::Ordinary, // transa - n_rows as i32, // m - i as i32, // n - 1.0, // alpha - q, // a - lda, // lda - q_column, // x - q_inc, // incx - 0.0, // beta - work_vector, //y - 1 // incy + self.memory_layout, + cblas::Transpose::Ordinary, + n_rows as i32, + i as i32, + 1.0, + q_matrix, + leading_dim, + q_column, + next_elem, + 0.0, + work_slice, + 1 // Always 1 from the definition of the work_slice/work_vector ); cblas::dgemv( - self.memory_order, // layout - cblas::Transpose::None, // transa - n_rows as i32, // m - i as i32, // n - -1.0, // alpha - q, // a - lda, // lda - work_vector, // y - 1, // incx - 1.0, // beta - q_column, // y - q_inc, // incy + self.memory_layout, + cblas::Transpose::None, + n_rows as i32, + i as i32, + -1.0, + q_matrix, + leading_dim, + work_slice, + 1, + 1.0, + q_column, + next_elem, ); cblas::daxpy( n_rows as i32, // n 1.0, // alpha - work_vector, // x - 1, // incx - r_column, // y - r_inc, // incy + work_slice, // x + 1, // Always 1 from the definition of the work_slice/work_vector + r_column, + next_elem, ); } }; - let norm = { - let len = self.q.len(); - let q_ptr = self.q.as_mut_ptr(); - unsafe { - let (q_column, q_inc) = match self.memory_order { - cblas::Layout::RowMajor => { - let offset = i as isize; - let q_column = slice::from_raw_parts_mut(q_ptr.offset(offset), len - i); - (q_column, n_cols as i32) - }, - - cblas::Layout::ColumnMajor => { - let offset = n_rows * i; - let q_column = slice::from_raw_parts_mut(q_ptr.offset(offset as isize), len - offset); - (q_column, 1) - }, - }; - cblas::dnrm2(n_rows as i32, q_column, q_inc) - } + let norm = unsafe { + cblas::dnrm2(n_rows as i32, q_column, next_elem) }; let mut v = self.q.column_mut(i); v /= norm; self.r[(i,i)] = a.column(i).dot(&v); } + + Ok(()) } - /// Return a reference to the matrix q. - pub fn q(&self) -> &Array2 { + fn q(&self) -> &Array2 { &self.q } - /// Return a reference to the matrix q. - pub fn r(&self) -> &Array2 { + fn r(&self) -> &Array2 { &self.r } } -generate_tests!(ReorthogonalizedGramSchmidt, 1e-13); +/// Convenience function that calculates a Reorthogonalized Gram Schmmidt QR factorization (see +/// [Giraud et al.] for details), returning a tuple `(Q,R)`. +/// +/// If you want to repeatedly calculate QR factorizations, then prefer constructing a +/// [`Reorthogonalized`] struct and calling its [`GramSchmidt::compute`] method implemented through +/// the [`GramSchmidt`] trait. +/// +/// [Giraud et al.]: https://doi.org/10.1007/s00211-005-0615-4 +/// [`Reorthogonalized`]: Reorthogonalized +/// [`GramSchmidt`]: GramSchmidt +/// [`GramSchmidt::compute`]: trait.GramSchmidt.html#method.compute +pub fn cgs2(a: &ArrayBase) -> Result<(Array, Array)> + where S: Data +{ + let mut cgs2 = Reorthogonalized::from_matrix(a)?; + cgs2.compute(a)?; + Ok((cgs2.q().clone(), cgs2.r().clone())) +} + +#[cfg(test)] +generate_tests!(Reorthogonalized, 1e-13); diff --git a/src/lib.rs b/src/lib.rs index 697ce06..fd78b01 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -1,3 +1,29 @@ +//! # Gram Schmidt procedures for Rust and `ndarray` +//! +//! This crate implements three different Gram Schmidt procedures in the form of QR decompositions: +//! +//! + The [classical Gram Schmidt] procedure, `[cgs]`; +//! + the [modified or stabilized Gram Schmidt] procedure, `[mgs]`; +//! + the [reorthogonalized Gram Schmidt procedure], `[cgs2]`. +//! +//! [ndarray]: https://github.com/rust-ndarray/ndarray +//! [classical Gram Schmidt]: https://en.wikipedia.org/wiki/Gram-Schmidt_process +//! [modified or stabilized Gram Schmidt]: https://en.wikipedia.org/wiki/Gram-Schmidt_process#Numerical_stabilty +//! [reorthogonalized Gram Schmidt procedure]: https://doi.org/10.1007/s00211-005-0615-4 +//! [cgs]: struct.Classical#method.cgs + +use ndarray::{ + ArrayBase, + Array2, + Data, + Ix2, + ShapeBuilder, +}; +use std::error; +use std::result; +use std::fmt; + +#[cfg(test)] #[macro_use] mod test_macros; @@ -5,9 +31,136 @@ mod cgs; mod cgs2; mod mgs; -#[cfg(test)] pub(crate) mod utils; -pub use cgs::ClassicalGramSchmidt; -pub use cgs2::ReorthogonalizedGramSchmidt; -pub use mgs::ModifiedGramSchmidt; +// Reexports +pub use cgs::{cgs, Classical}; +pub use cgs2::{cgs2, Reorthogonalized}; +pub use mgs::{mgs, Modified}; + +/// Errors that occur during a initialization of a Gram Schmidt factorization. +#[derive(Debug)] +pub enum Error { + /// The layout of the matrix to be factorized is incompatible with the layout the GramSchmidt + /// procedure was configured for. It means that the GramSchmidt procedure is configured to + /// work with either column major (Fortran layout) or row major (C layout) matrices. + IncompatibleLayouts, + + /// The array to be factorized is not contiguous. At the moment, all arrays to be factorized + /// have to be contiguous. + NonContiguous, +} + +pub type Result = result::Result; + +impl fmt::Display for Error { + fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { + use Error::*; + match self { + IncompatibleLayouts => write!(f, "The arrays representing the matrices don't have the same layouts."), + NonContiguous => write!(f, "Array shape is not contiguous"), + } + } +} + +impl error::Error for Error { + fn source(&self) -> Option<&(dyn error::Error + 'static)> { + None + } +} + +pub trait GramSchmidt: Sized { + /// Reserves the memory for a QR decomposition via a classical Gram Schmidt orthogonalization + /// using a shape. + /// + /// The resulting object can be used to orthogonalize matrices of the same dimensions. + /// + /// # Example + /// + /// ``` + /// use gramschmidt::{ + /// Classical, + /// GramSchmidt, + /// Result, + /// }; + /// + /// # fn main() -> Result<()> { + /// + /// let mut cgs = Classical::from_shape((10,10))?; + /// + /// # Ok(()) + /// # } + /// ``` + fn from_shape(shape: T) -> Result + where T: ShapeBuilder; + + /// Computes a QR decomposition using the classical Gram Schmidt orthonormalization of the + /// matrix `a`. + /// + /// The input matrix `a` has to have exactly the same dimension and memory layout as was + /// previously configured. Panics otherwise. + /// + /// ``` + /// extern crate openblas_src; + /// + /// use gramschmidt::{GramSchmidt, Classical}; + /// use ndarray::Array2; + /// use ndarray_rand::RandomExt; + /// use rand::distributions::Normal; + /// + /// # fn main() { + /// + /// let matrix = Array2::random((10,10), Normal::new(0.0, 1.0)); + /// let mut cgs = Classical::from_matrix(&matrix).unwrap(); + /// cgs.compute(&matrix); + /// + /// # } + /// ``` + fn compute(&mut self, a: &ArrayBase) -> Result<()> + where S: Data; + + /// Return a reference to the matrix q. + fn q(&self) -> &Array2; + + /// Return a reference to the matrix q. + fn r(&self) -> &Array2; + + // Blanket impls + + /// Uses a matrix to reserve memory for a QR decomposition via a classical Gram Schmidt. + /// + /// The resulting object can be used to orthogonalize matrices of the same dimensions. + /// + /// # Example + /// + /// ``` + /// use ndarray::Array; + /// use gramschmidt::{ + /// Classical, + /// GramSchmidt, + /// Result, + /// }; + /// + /// # fn main() -> Result<()> { + /// + /// let a = Array::zeros((10, 10)); + /// let mut cgs = Classical::from_matrix(&a)?; + /// + /// # Ok(()) + /// # } + /// ``` + fn from_matrix(a: &ArrayBase) -> Result + where S: Data + { + use cblas::Layout::*; + let dim = a.dim(); + let shape = match utils::get_layout(a) { + Some(ColumnMajor) => dim.f(), + Some(RowMajor) => dim.into_shape(), + None => Err(Error::NonContiguous)?, + }; + + Self::from_shape(shape) + } + +} diff --git a/src/mgs.rs b/src/mgs.rs index 448095e..e54b855 100644 --- a/src/mgs.rs +++ b/src/mgs.rs @@ -2,106 +2,54 @@ use cblas; use ndarray::{ Data, - IntoDimension, + ShapeBuilder, }; use ndarray::prelude::*; use std::slice; +use crate::{ + Error, + GramSchmidt, + Result, + utils::get_layout, +}; + +/// A modified Gram Schmidt factorization, which has a better numerical stability compared to +/// the classical Gram Schmidt procedure. See its [Wikipedia entry] for more information. +/// +/// Use this struct via the [`GramSchmidt` trait]. +/// +/// [Wikipedia entry]: https://en.wikipedia.org/wiki/Gram-Schmidt_process#Numerical_stabilty +/// [`GramSchmidt` trait]: GramSchmidt #[derive(Clone, Debug)] -pub struct ModifiedGramSchmidt { +pub struct Modified { q: Array2, r: Array2, - memory_order: cblas::Layout, + memory_layout: cblas::Layout, } -impl ModifiedGramSchmidt { - /// Reserves the memory for a QR decomposition via a modified Gram Schmidt orthogonalization - /// using the dimensions of a sample matrix. - /// - /// The resulting object can be used to orthogonalize matrices of the same dimensions. - /// - /// # Example - /// - /// ``` - /// use gramschmidt::ModifiedGramSchmidt; - /// use ndarray::Array2; - /// use ndarray_rand::RandomExt; - /// use rand::distributions::Normal; - /// - /// # fn main() { - /// let matrix = Array2::random((10,10), Normal::new(0.0, 1.0)); - /// let mut cgs = ModifiedGramSchmidt::from_matrix(&matrix); - /// # } - /// ``` - pub fn from_matrix(a: &ArrayBase) -> Self - where S: Data - { - let (is_fortran, memory_order) = if let Some(_) = a.as_slice() { - (false, cblas::Layout::RowMajor) - } else if let Some(_) = a.as_slice_memory_order() { - (true, cblas::Layout::ColumnMajor) - } else { - panic!("Array not contiguous!") - }; - - let array_shape = a.raw_dim().set_f(is_fortran); - ModifiedGramSchmidt { - q: Array2::zeros(array_shape), - r: Array2::zeros(array_shape), - memory_order, - } - } - - /// Reserves the memory for a QR decomposition via a modified Gram Schmidt orthogonalization - /// using a shape. - /// - /// The resulting object can be used to orthogonalize matrices of the same dimensions. - /// - /// # Example - /// - /// ``` - /// use gramschmidt::ModifiedGramSchmidt; - /// - /// # fn main() { - /// let fortran_order = false; - /// let mut cgs = ModifiedGramSchmidt::from_shape((10,10), fortran_order); - /// # } - /// ``` - pub fn from_shape(shape: T, fortran_order: bool) -> Self - where T: IntoDimension, +impl GramSchmidt for Modified { + fn from_shape(shape: T) -> Result + where T: ShapeBuilder, { - let dimension = shape.into_dimension(); - let array_shape = dimension.set_f(fortran_order); - let memory_order = if fortran_order { - cblas::Layout::ColumnMajor - } else { - cblas::Layout::RowMajor + // Unfortunately we cannot check the shape itself to see if it's + // in ColumnMajor or RowMajor layout. So we need to first construct + // an array and then check that. + let shape = shape.into_shape(); + let q = Array2::zeros(shape); + let memory_layout = match get_layout(&q) { + Some(layout) => layout, + None => Err(Error::NonContiguous)?, }; - ModifiedGramSchmidt { - q: Array2::zeros(array_shape), - r: Array2::zeros(array_shape), - memory_order, - } + let r = q.clone(); + Ok(Self { + q, + r, + memory_layout, + }) } - /// Computes a QR decomposition using the modified Gram Schmidt orthonormalization of the - /// matrix `a`. - /// - /// The input matrix `a` has to have exactly the same dimension and memory layout as was - /// previously configured. Panics otherwise. - /// - /// ``` - /// use gramschmidt::ModifiedGramSchmidt; - /// use ndarray::Array2; - /// use ndarray_rand::RandomExt; - /// use rand::distributions::Normal; - /// - /// # fn main() { - /// let matrix = Array2::random((10,10), Normal::new(0.0, 1.0)); - /// let mut cgs = ModifiedGramSchmidt::from_matrix(&matrix); - /// # } - /// ``` - pub fn compute(&mut self, a: &ArrayBase) + fn compute(&mut self, a: &ArrayBase) -> Result<()> where S: Data, { let (n_rows, n_cols) = a.dim(); @@ -128,7 +76,7 @@ impl ModifiedGramSchmidt { let len = self.q.len(); let q_ptr = self.q.as_mut_ptr(); unsafe { - let (q_column, q_inc) = match self.memory_order { + let (q_column, q_inc) = match self.memory_layout { cblas::Layout::RowMajor => { let offset = i as isize; let q_column = slice::from_raw_parts_mut(q_ptr.offset(offset), len - i); @@ -149,17 +97,37 @@ impl ModifiedGramSchmidt { let mut q_column = self.q.column_mut(i); q_column /= norm; } + + Ok(()) } - /// Return a reference to the matrix q. - pub fn q(&self) -> &Array2 { + fn q(&self) -> &Array2 { &self.q } - /// Return a reference to the matrix q. - pub fn r(&self) -> &Array2 { + fn r(&self) -> &Array2 { &self.r } } -generate_tests!(ModifiedGramSchmidt, 1e-13); +/// Convenience function that calculates a [Modified Gram Schmidt] QR factorization, returning a +/// tuple `(Q,R)`. +/// +/// If you want to repeatedly calculate QR factorizations, then prefer constructing a +/// [`Modified`] struct and calling its [`GramSchmidt::compute`] method implemented through +/// the [`GramSchmidt`] trait. +/// +/// [Modified Gram Schmidt]: https://en.wikipedia.org/wiki/Gram-Schmidt_process#Numerical_stabilty +/// [`Modified`]: Modified +/// [`GramSchmidt`]: GramSchmidt +/// [`GramSchmidt::compute`]: trait.GramSchmidt.html#tymethod.compute +pub fn mgs(a: &ArrayBase) -> Result<(Array, Array)> + where S: Data +{ + let mut mgs = Modified::from_matrix(a)?; + mgs.compute(a)?; + Ok((mgs.q().clone(), mgs.r().clone())) +} + +#[cfg(test)] +generate_tests!(Modified, 1e-13); diff --git a/src/test_macros.rs b/src/test_macros.rs index 828be46..69bbdbe 100644 --- a/src/test_macros.rs +++ b/src/test_macros.rs @@ -72,73 +72,73 @@ macro_rules! generate_tests { #[test] fn unity_stays_unity() { - let mut method = $method::from_matrix(&*UNITY); - method.compute(&*UNITY); + let mut method = $method::from_matrix(&*UNITY).unwrap(); + assert!(method.compute(&*UNITY).is_ok()); assert_eq!(&*UNITY, &method.q().dot(method.r())); } #[test] fn small_orthogonal() { - let mut method = $method::from_matrix(&*SMALL); - method.compute(&*SMALL); + let mut method = $method::from_matrix(&*SMALL).unwrap(); + assert!(method.compute(&*SMALL).is_ok()); assert!(crate::utils::orthogonal(method.q(),$tolerance)); } #[test] fn small_qr_returns_original() { - let mut method = $method::from_matrix(&*SMALL); - method.compute(&*SMALL); + let mut method = $method::from_matrix(&*SMALL).unwrap(); + assert!(method.compute(&*SMALL).is_ok()); assert!(SMALL.all_close(&method.q().dot(method.r()), $tolerance)); } #[test] fn large_orthogonal() { - let mut method = $method::from_matrix(&*LARGE); - method.compute(&*LARGE); + let mut method = $method::from_matrix(&*LARGE).unwrap(); + assert!(method.compute(&*LARGE).is_ok()); assert!(crate::utils::orthogonal(method.q(),$tolerance)); } #[test] fn large_qr_returns_original() { - let mut method = $method::from_matrix(&*LARGE); - method.compute(&*LARGE); + let mut method = $method::from_matrix(&*LARGE).unwrap(); + assert!(method.compute(&*LARGE).is_ok()); assert!(LARGE.all_close(&method.q().dot(method.r()), $tolerance)); } #[test] fn f_order_unity_stays_unity() { - let mut method = $method::from_matrix(&*F_UNITY); - method.compute(&*F_UNITY); + let mut method = $method::from_matrix(&*F_UNITY).unwrap(); + assert!(method.compute(&*F_UNITY).is_ok()); assert_eq!(&*F_UNITY, &method.q().dot(method.r())); } #[test] fn f_order_small_orthogonal() { - let mut method = $method::from_matrix(&*F_SMALL); - method.compute(&*F_SMALL); + let mut method = $method::from_matrix(&*F_SMALL).unwrap(); + assert!(method.compute(&*F_SMALL).is_ok()); assert!(crate::utils::orthogonal(method.q(),$tolerance)); } #[test] fn f_order_small_qr_returns_original() { - let mut method = $method::from_matrix(&*F_SMALL); - method.compute(&*F_SMALL); + let mut method = $method::from_matrix(&*F_SMALL).unwrap(); + assert!(method.compute(&*F_SMALL).is_ok()); assert!(F_SMALL.all_close(&method.q().dot(method.r()), $tolerance)); } #[test] fn f_order_large_orthogonal() { - let mut method = $method::from_matrix(&*F_LARGE); - method.compute(&*F_LARGE); + let mut method = $method::from_matrix(&*F_LARGE).unwrap(); + assert!(method.compute(&*F_LARGE).is_ok()); assert!(crate::utils::orthogonal(method.q(),$tolerance)); } #[test] fn f_order_large_qr_returns_original() { - let mut method = $method::from_matrix(&*F_LARGE); - method.compute(&*F_LARGE); + let mut method = $method::from_matrix(&*F_LARGE).unwrap(); + assert!(method.compute(&*F_LARGE).is_ok()); assert!(F_LARGE.all_close(&method.q().dot(method.r()), $tolerance)); } } diff --git a/src/utils.rs b/src/utils.rs index 915e8f3..4e7d5eb 100644 --- a/src/utils.rs +++ b/src/utils.rs @@ -1,9 +1,38 @@ use ndarray::Data; use ndarray::prelude::*; +#[cfg(test)] pub(crate) fn orthogonal(a: &ArrayBase, tol: f64) -> bool where S: Data { let b = a.dot(&a.t()); b.all_close(&Array2::eye(b.shape()[0]), tol) } + +/// Returns slice and layout underlying an array `a`. +pub(crate) fn get_layout(a: &ArrayBase) -> Option + where S: Data, + D: Dimension +{ + if let Some(_) = a.as_slice() { + Some(cblas::Layout::RowMajor) + } else if let Some(_) = a.as_slice_memory_order() { + Some(cblas::Layout::ColumnMajor) + } else { + None + } +} + +/// Returns slice and layout underlying an array `a`. +pub(crate) fn as_slice_with_layout(a: &ArrayBase) -> Option<(&[T], cblas::Layout)> + where S: Data, + D: Dimension +{ + if let Some(a_slice) = a.as_slice() { + Some((a_slice, cblas::Layout::RowMajor)) + } else if let Some(a_slice) = a.as_slice_memory_order() { + Some((a_slice, cblas::Layout::ColumnMajor)) + } else { + None + } +} From d449fd32d181e1dbec1f9e42edd106665ffadee1 Mon Sep 17 00:00:00 2001 From: Richard Janis Goldschmidt Date: Tue, 12 Mar 2019 12:25:12 +0100 Subject: [PATCH 3/3] Version 0.5.0 --- Cargo.toml | 2 +- README.md | 21 +++++++++++++++------ 2 files changed, 16 insertions(+), 7 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index ee50663..3d60897 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "gramschmidt" -version = "0.4.1" +version = "0.5.0" edition = "2018" authors = ["Richard Janis Goldschmidt "] license = "MIT" diff --git a/README.md b/README.md index d839180..7c12e8b 100644 --- a/README.md +++ b/README.md @@ -11,27 +11,36 @@ This crate provides the following methods: # Usage ```rust -extern crate gramschmidt; -extern crate ndarray; - // Import openblas_src or another blas source to have the linker find all symbols. extern crate openblas_src; -fn main() { +use gramschmidt::{ + GramSchmidt, + Reorthogonalized, + Result, +}; +use ndarray::arr2; + +fn main() -> Result<()> { let small_matrix = arr2( &[[2.0, 0.5, 0.0, 0.0], [0.0, 0.3, 0.0, 0.0], [0.0, 1.0, 0.7, 0.0], [0.0, 0.0, 0.0, 3.0]] ); - let mut cgs2 = ReorthogonalizedGramSchmidt::from_matrix(&small_matrix); - cgs2.compute(&small_matrix); + let mut cgs2 = Reorthogonalized::from_matrix(&small_matrix)?; + cgs2.compute(&small_matrix)?; assert!(small_matrix.all_close(&cgs2.q().dot(cgs2.r()), 1e-14)); + Ok(()) } ``` # Recent versions ++ `0.5.0`: Refactored the library and updated for edition 2018 + + the Gram Schmidt factorizations are now all implemented via the `GramSchmidt` trait; + + introduce some error handling; + + provide convenience functions `cgs`, `cgs2`, and `mgs`. + `0.4.1`: Fixed doc tests and expanded + simplified tests. + `0.4.0`: Major rework of the library structure: + The algorithms are now configured via structs, the traits are dropped.