From d7cf5037100c477474e97fc132b9f60fba2c3281 Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Wed, 15 Jul 2020 17:49:33 +0900 Subject: [PATCH 01/11] Split impl for real based on LAPACK --- lax/src/least_squares.rs | 167 ++++++++++++++++++++++++++++++++++++++- 1 file changed, 164 insertions(+), 3 deletions(-) diff --git a/lax/src/least_squares.rs b/lax/src/least_squares.rs index 69553a44..269fa952 100644 --- a/lax/src/least_squares.rs +++ b/lax/src/least_squares.rs @@ -2,7 +2,7 @@ use crate::{error::*, layout::MatrixLayout}; use cauchy::*; -use num_traits::Zero; +use num_traits::{ToPrimitive, Zero}; /// Result of LeastSquares pub struct LeastSquaresOutput { @@ -28,6 +28,169 @@ pub trait LeastSquaresSvdDivideConquer_: Scalar { ) -> Result>; } +/// Eval iwork size +/// +/// - NLVL = INT( LOG_2( MIN( M, N ) / ( SMLSIZ + 1 ) ) ) + 1 +/// - LIWORK = 3 * MIN( M, N ) * NLVL + 11 * MIN( M, N ) +/// +/// where SMLSIZ is returned by ILAENV and is equal to the maximum size of the subproblems +/// at the bottom of the computation tree (usually about 25). +/// +/// We put SMLSIZ=0 to estimate NLVL as large as possible +/// because its size will usually very small. +fn iwork_size(m: i32, n: i32) -> usize { + let mn = m.min(n) as usize; + let nlvl = (mn.to_f32().unwrap().log2() + 1.0) + .max(0.0) + .to_usize() + .unwrap(); + (3 * mn * nlvl + 11 * mn).max(1) +} + +macro_rules! impl_least_squares_real { + ($scalar:ty, $gelsd:path) => { + impl LeastSquaresSvdDivideConquer_ for $scalar { + unsafe fn least_squares( + l: MatrixLayout, + a: &mut [Self], + b: &mut [Self], + ) -> Result> { + let m = l.lda(); + let n = l.len(); + let k = m.min(n); + if (m as usize) > b.len() || (n as usize) > b.len() { + return Err(Error::InvalidShape); + } + let rcond: Self::Real = -1.; + let mut singular_values: Vec = vec![Self::Real::zero(); k as usize]; + let mut rank: i32 = 0; + + let mut iwork = vec![0; iwork_size(m, n)]; + + // eval work size + let mut info = 0; + let mut work_size = [Self::zero()]; + $gelsd( + m, + n, + 1, // nrhs + a, + m, + b, + b.len() as i32, + &mut singular_values, + rcond, + &mut rank, + &mut work_size, + -1, + &mut iwork, + &mut info, + ); + info.as_lapack_result()?; + + // calc + let lwork = work_size[0].to_usize().unwrap(); + let mut work = vec![Self::zero(); lwork]; + $gelsd( + m, + n, + 1, // nrhs + a, + m, + b, + b.len() as i32, + &mut singular_values, + rcond, + &mut rank, + &mut work, + lwork as i32, + &mut iwork, + &mut info, + ); + info.as_lapack_result()?; + + Ok(LeastSquaresOutput { + singular_values, + rank, + }) + } + + unsafe fn least_squares_nrhs( + a_layout: MatrixLayout, + a: &mut [Self], + b_layout: MatrixLayout, + b: &mut [Self], + ) -> Result> { + let m = a_layout.lda(); + let n = a_layout.len(); + let k = m.min(n); + if (m as usize) > b.len() + || (n as usize) > b.len() + || a_layout.lapacke_layout() != b_layout.lapacke_layout() + { + return Err(Error::InvalidShape); + } + let (b_lda, nrhs) = b_layout.size(); + let rcond: Self::Real = -1.; + let mut singular_values: Vec = vec![Self::Real::zero(); k as usize]; + let mut rank: i32 = 0; + + let mut iwork = vec![0; iwork_size(m, n)]; + + // eval work size + let mut info = 0; + let mut work_size = [Self::zero()]; + $gelsd( + m, + n, + nrhs, + a, + m, + b, + b_lda, + &mut singular_values, + rcond, + &mut rank, + &mut work_size, + -1, + &mut iwork, + &mut info, + ); + info.as_lapack_result()?; + + // calc + let lwork = work_size[0].to_usize().unwrap(); + let mut work = vec![Self::zero(); lwork]; + $gelsd( + m, + n, + nrhs, + a, + m, + b, + b_lda, + &mut singular_values, + rcond, + &mut rank, + &mut work, + lwork as i32, + &mut iwork, + &mut info, + ); + info.as_lapack_result()?; + + Ok(LeastSquaresOutput { + singular_values, + rank, + }) + } + } + }; +} + +impl_least_squares_real!(f64, lapack::dgelsd); +impl_least_squares_real!(f32, lapack::sgelsd); + macro_rules! impl_least_squares { ($scalar:ty, $gelsd:path) => { impl LeastSquaresSvdDivideConquer_ for $scalar { @@ -113,7 +276,5 @@ macro_rules! impl_least_squares { }; } -impl_least_squares!(f64, lapacke::dgelsd); -impl_least_squares!(f32, lapacke::sgelsd); impl_least_squares!(c64, lapacke::zgelsd); impl_least_squares!(c32, lapacke::cgelsd); From 45c717031d18115eb504c22eb2bf78d4dcadb068 Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Wed, 15 Jul 2020 19:54:25 +0900 Subject: [PATCH 02/11] Query liwork --- lax/src/least_squares.rs | 33 ++++++++------------------------- 1 file changed, 8 insertions(+), 25 deletions(-) diff --git a/lax/src/least_squares.rs b/lax/src/least_squares.rs index 269fa952..102ac261 100644 --- a/lax/src/least_squares.rs +++ b/lax/src/least_squares.rs @@ -28,25 +28,6 @@ pub trait LeastSquaresSvdDivideConquer_: Scalar { ) -> Result>; } -/// Eval iwork size -/// -/// - NLVL = INT( LOG_2( MIN( M, N ) / ( SMLSIZ + 1 ) ) ) + 1 -/// - LIWORK = 3 * MIN( M, N ) * NLVL + 11 * MIN( M, N ) -/// -/// where SMLSIZ is returned by ILAENV and is equal to the maximum size of the subproblems -/// at the bottom of the computation tree (usually about 25). -/// -/// We put SMLSIZ=0 to estimate NLVL as large as possible -/// because its size will usually very small. -fn iwork_size(m: i32, n: i32) -> usize { - let mn = m.min(n) as usize; - let nlvl = (mn.to_f32().unwrap().log2() + 1.0) - .max(0.0) - .to_usize() - .unwrap(); - (3 * mn * nlvl + 11 * mn).max(1) -} - macro_rules! impl_least_squares_real { ($scalar:ty, $gelsd:path) => { impl LeastSquaresSvdDivideConquer_ for $scalar { @@ -65,11 +46,10 @@ macro_rules! impl_least_squares_real { let mut singular_values: Vec = vec![Self::Real::zero(); k as usize]; let mut rank: i32 = 0; - let mut iwork = vec![0; iwork_size(m, n)]; - // eval work size let mut info = 0; let mut work_size = [Self::zero()]; + let mut iwork_size = [0]; $gelsd( m, n, @@ -83,7 +63,7 @@ macro_rules! impl_least_squares_real { &mut rank, &mut work_size, -1, - &mut iwork, + &mut iwork_size, &mut info, ); info.as_lapack_result()?; @@ -91,6 +71,8 @@ macro_rules! impl_least_squares_real { // calc let lwork = work_size[0].to_usize().unwrap(); let mut work = vec![Self::zero(); lwork]; + let liwork = iwork_size[0].to_usize().unwrap(); + let mut iwork = vec![0; liwork]; $gelsd( m, n, @@ -135,11 +117,10 @@ macro_rules! impl_least_squares_real { let mut singular_values: Vec = vec![Self::Real::zero(); k as usize]; let mut rank: i32 = 0; - let mut iwork = vec![0; iwork_size(m, n)]; - // eval work size let mut info = 0; let mut work_size = [Self::zero()]; + let mut iwork_size = [0]; $gelsd( m, n, @@ -153,7 +134,7 @@ macro_rules! impl_least_squares_real { &mut rank, &mut work_size, -1, - &mut iwork, + &mut iwork_size, &mut info, ); info.as_lapack_result()?; @@ -161,6 +142,8 @@ macro_rules! impl_least_squares_real { // calc let lwork = work_size[0].to_usize().unwrap(); let mut work = vec![Self::zero(); lwork]; + let liwork = iwork_size[0].to_usize().unwrap(); + let mut iwork = vec![0; liwork]; $gelsd( m, n, From 6df2779080f70755ff51938e720ba6adb0c3040a Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Wed, 15 Jul 2020 21:11:35 +0900 Subject: [PATCH 03/11] Use least_squares_nrhs for 1-dim b case --- lax/src/least_squares.rs | 61 ++-------------------------------------- 1 file changed, 2 insertions(+), 59 deletions(-) diff --git a/lax/src/least_squares.rs b/lax/src/least_squares.rs index 102ac261..468a759b 100644 --- a/lax/src/least_squares.rs +++ b/lax/src/least_squares.rs @@ -36,65 +36,8 @@ macro_rules! impl_least_squares_real { a: &mut [Self], b: &mut [Self], ) -> Result> { - let m = l.lda(); - let n = l.len(); - let k = m.min(n); - if (m as usize) > b.len() || (n as usize) > b.len() { - return Err(Error::InvalidShape); - } - let rcond: Self::Real = -1.; - let mut singular_values: Vec = vec![Self::Real::zero(); k as usize]; - let mut rank: i32 = 0; - - // eval work size - let mut info = 0; - let mut work_size = [Self::zero()]; - let mut iwork_size = [0]; - $gelsd( - m, - n, - 1, // nrhs - a, - m, - b, - b.len() as i32, - &mut singular_values, - rcond, - &mut rank, - &mut work_size, - -1, - &mut iwork_size, - &mut info, - ); - info.as_lapack_result()?; - - // calc - let lwork = work_size[0].to_usize().unwrap(); - let mut work = vec![Self::zero(); lwork]; - let liwork = iwork_size[0].to_usize().unwrap(); - let mut iwork = vec![0; liwork]; - $gelsd( - m, - n, - 1, // nrhs - a, - m, - b, - b.len() as i32, - &mut singular_values, - rcond, - &mut rank, - &mut work, - lwork as i32, - &mut iwork, - &mut info, - ); - info.as_lapack_result()?; - - Ok(LeastSquaresOutput { - singular_values, - rank, - }) + let b_layout = l.resized(b.len() as i32, 1); + Self::least_squares_nrhs(l, a, b_layout, b) } unsafe fn least_squares_nrhs( From f417f6f15e19e8c09abecc9a6ac88ba77ea222ca Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Fri, 24 Jul 2020 21:48:12 +0900 Subject: [PATCH 04/11] Fix test scalar types --- ndarray-linalg/tests/least_squares.rs | 12 +++--- ndarray-linalg/tests/least_squares_nrhs.rs | 48 +++++++++++----------- 2 files changed, 30 insertions(+), 30 deletions(-) diff --git a/ndarray-linalg/tests/least_squares.rs b/ndarray-linalg/tests/least_squares.rs index c388c9d7..e2df3370 100644 --- a/ndarray-linalg/tests/least_squares.rs +++ b/ndarray-linalg/tests/least_squares.rs @@ -27,13 +27,13 @@ macro_rules! impl_exact { paste::item! { #[test] fn []() { - let a: Array2 = random((3, 3)); + let a: Array2<$scalar> = random((3, 3)); test_exact(a) } #[test] fn []() { - let a: Array2 = random((3, 3).f()); + let a: Array2<$scalar> = random((3, 3).f()); test_exact(a) } } @@ -73,13 +73,13 @@ macro_rules! impl_overdetermined { paste::item! { #[test] fn []() { - let a: Array2 = random((4, 3)); + let a: Array2<$scalar> = random((4, 3)); test_overdetermined(a) } #[test] fn []() { - let a: Array2 = random((4, 3).f()); + let a: Array2<$scalar> = random((4, 3).f()); test_overdetermined(a) } } @@ -110,13 +110,13 @@ macro_rules! impl_underdetermined { paste::item! { #[test] fn []() { - let a: Array2 = random((3, 4)); + let a: Array2<$scalar> = random((3, 4)); test_underdetermined(a) } #[test] fn []() { - let a: Array2 = random((3, 4).f()); + let a: Array2<$scalar> = random((3, 4).f()); test_underdetermined(a) } } diff --git a/ndarray-linalg/tests/least_squares_nrhs.rs b/ndarray-linalg/tests/least_squares_nrhs.rs index 4c964697..8072c409 100644 --- a/ndarray-linalg/tests/least_squares_nrhs.rs +++ b/ndarray-linalg/tests/least_squares_nrhs.rs @@ -31,8 +31,8 @@ macro_rules! impl_exact { paste::item! { #[test] fn []() { - let a: Array2 = random((3, 3)); - let b: Array2 = random((3, 2)); + let a: Array2<$scalar> = random((3, 3)); + let b: Array2<$scalar> = random((3, 2)); test_exact(a, b) } @@ -40,15 +40,15 @@ macro_rules! impl_exact { #[test] fn []() { - let a: Array2 = random((3, 3)); - let b: Array2 = random((3, 2).f()); + let a: Array2<$scalar> = random((3, 3)); + let b: Array2<$scalar> = random((3, 2).f()); test_exact(a, b) } #[test] fn []() { - let a: Array2 = random((3, 3).f()); - let b: Array2 = random((3, 2)); + let a: Array2<$scalar> = random((3, 3).f()); + let b: Array2<$scalar> = random((3, 2)); test_exact(a, b) } @@ -56,8 +56,8 @@ macro_rules! impl_exact { #[test] fn []() { - let a: Array2 = random((3, 3).f()); - let b: Array2 = random((3, 2).f()); + let a: Array2<$scalar> = random((3, 3).f()); + let b: Array2<$scalar> = random((3, 2).f()); test_exact(a, b) } } @@ -103,8 +103,8 @@ macro_rules! impl_overdetermined { paste::item! { #[test] fn []() { - let a: Array2 = random((4, 3)); - let b: Array2 = random((4, 2)); + let a: Array2<$scalar> = random((4, 3)); + let b: Array2<$scalar> = random((4, 2)); test_overdetermined(a, b) } @@ -112,15 +112,15 @@ macro_rules! impl_overdetermined { #[test] fn []() { - let a: Array2 = random((4, 3).f()); - let b: Array2 = random((4, 2)); + let a: Array2<$scalar> = random((4, 3).f()); + let b: Array2<$scalar> = random((4, 2)); test_overdetermined(a, b) } #[test] fn []() { - let a: Array2 = random((4, 3)); - let b: Array2 = random((4, 2).f()); + let a: Array2<$scalar> = random((4, 3)); + let b: Array2<$scalar> = random((4, 2).f()); test_overdetermined(a, b) } @@ -128,8 +128,8 @@ macro_rules! impl_overdetermined { #[test] fn []() { - let a: Array2 = random((4, 3).f()); - let b: Array2 = random((4, 2).f()); + let a: Array2<$scalar> = random((4, 3).f()); + let b: Array2<$scalar> = random((4, 2).f()); test_overdetermined(a, b) } } @@ -162,8 +162,8 @@ macro_rules! impl_underdetermined { paste::item! { #[test] fn []() { - let a: Array2 = random((3, 4)); - let b: Array2 = random((3, 2)); + let a: Array2<$scalar> = random((3, 4)); + let b: Array2<$scalar> = random((3, 2)); test_underdetermined(a, b) } @@ -171,15 +171,15 @@ macro_rules! impl_underdetermined { #[test] fn []() { - let a: Array2 = random((3, 4).f()); - let b: Array2 = random((3, 2)); + let a: Array2<$scalar> = random((3, 4).f()); + let b: Array2<$scalar> = random((3, 2)); test_underdetermined(a, b) } #[test] fn []() { - let a: Array2 = random((3, 4)); - let b: Array2 = random((3, 2).f()); + let a: Array2<$scalar> = random((3, 4)); + let b: Array2<$scalar> = random((3, 2).f()); test_underdetermined(a, b) } @@ -187,8 +187,8 @@ macro_rules! impl_underdetermined { #[test] fn []() { - let a: Array2 = random((3, 4).f()); - let b: Array2 = random((3, 2).f()); + let a: Array2<$scalar> = random((3, 4).f()); + let b: Array2<$scalar> = random((3, 2).f()); test_underdetermined(a, b) } } From 8ff555d2ca17fcccae3af15e6b9444547996f0dc Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Sat, 25 Jul 2020 02:29:21 +0900 Subject: [PATCH 05/11] Out-place transpose --- lax/src/layout.rs | 87 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 87 insertions(+) diff --git a/lax/src/layout.rs b/lax/src/layout.rs index 9dad70e6..43b7ee87 100644 --- a/lax/src/layout.rs +++ b/lax/src/layout.rs @@ -97,6 +97,37 @@ impl MatrixLayout { MatrixLayout::F { col, lda } => MatrixLayout::C { row: lda, lda: col }, } } + + /// Transpose without changing memory representation + /// + /// C-contigious row=2, lda=3 + /// + /// ```text + /// [[1, 2, 3] + /// [4, 5, 6]] + /// ``` + /// + /// and F-contigious col=2, lda=3 + /// + /// ```text + /// [[1, 4] + /// [2, 5] + /// [3, 6]] + /// ``` + /// + /// have same memory representation `[1, 2, 3, 4, 5, 6]`, and this toggles them. + /// + /// ``` + /// # use lax::layout::*; + /// let layout = MatrixLayout::C { row: 2, lda: 3 }; + /// assert_eq!(layout.t(), MatrixLayout::F { col: 2, lda: 3 }); + /// ``` + pub fn t(&self) -> Self { + match *self { + MatrixLayout::C { row, lda } => MatrixLayout::F { col: row, lda }, + MatrixLayout::F { col, lda } => MatrixLayout::C { row: col, lda }, + } + } } /// In-place transpose of a square matrix by keeping F/C layout @@ -139,3 +170,59 @@ pub fn square_transpose(layout: MatrixLayout, a: &mut [T]) { } } } + +/// Out-place transpose for general matrix +/// +/// Inplace transpose of non-square matrices is hard. +/// See also: https://en.wikipedia.org/wiki/In-place_matrix_transposition +/// +/// ```rust +/// # use lax::layout::*; +/// let layout = MatrixLayout::C { row: 2, lda: 3 }; +/// let a = vec![1., 2., 3., 4., 5., 6.]; +/// let mut b = vec![0.0; a.len()]; +/// let l = transpose(layout, &a, &mut b); +/// assert_eq!(l, MatrixLayout::F { col: 3, lda: 2 }); +/// assert_eq!(b, &[1., 4., 2., 5., 3., 6.]); +/// ``` +/// +/// ```rust +/// # use lax::layout::*; +/// let layout = MatrixLayout::F { col: 2, lda: 3 }; +/// let a = vec![1., 2., 3., 4., 5., 6.]; +/// let mut b = vec![0.0; a.len()]; +/// let l = transpose(layout, &a, &mut b); +/// assert_eq!(l, MatrixLayout::C { row: 3, lda: 2 }); +/// assert_eq!(b, &[1., 4., 2., 5., 3., 6.]); +/// ``` +/// +/// Panics +/// ------ +/// - If size of `a` and `layout` size mismatch +/// +pub fn transpose(layout: MatrixLayout, from: &[T], to: &mut [T]) -> MatrixLayout { + let (m, n) = layout.size(); + let transposed = layout.resized(n, m).t(); + let m = m as usize; + let n = n as usize; + assert_eq!(from.len(), m * n); + assert_eq!(to.len(), m * n); + + match layout { + MatrixLayout::C { .. } => { + for i in 0..m { + for j in 0..n { + to[j * m + i] = from[i * n + j]; + } + } + } + MatrixLayout::F { .. } => { + for i in 0..m { + for j in 0..n { + to[i * n + j] = from[j * m + i]; + } + } + } + } + transposed +} From 200b9d435a3a0d1adb562e7956841cbf7a0b9689 Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Sat, 25 Jul 2020 17:07:06 +0900 Subject: [PATCH 06/11] Take transpose --- lax/src/least_squares.rs | 67 +++++++++++++++------- ndarray-linalg/src/least_squares.rs | 1 + ndarray-linalg/tests/least_squares_nrhs.rs | 1 + 3 files changed, 49 insertions(+), 20 deletions(-) diff --git a/lax/src/least_squares.rs b/lax/src/least_squares.rs index 468a759b..64349051 100644 --- a/lax/src/least_squares.rs +++ b/lax/src/least_squares.rs @@ -1,6 +1,6 @@ //! Least squares -use crate::{error::*, layout::MatrixLayout}; +use crate::{error::*, layout::*}; use cauchy::*; use num_traits::{ToPrimitive, Zero}; @@ -46,16 +46,37 @@ macro_rules! impl_least_squares_real { b_layout: MatrixLayout, b: &mut [Self], ) -> Result> { - let m = a_layout.lda(); - let n = a_layout.len(); + // Minimize |b - Ax|_2 + // + // where + // A : (m, n) + // b : (m, p) + // x : (n, p) + let (m, n) = a_layout.size(); + let (m_, p) = b_layout.size(); let k = m.min(n); - if (m as usize) > b.len() - || (n as usize) > b.len() - || a_layout.lapacke_layout() != b_layout.lapacke_layout() - { - return Err(Error::InvalidShape); - } - let (b_lda, nrhs) = b_layout.size(); + assert_eq!(m, m_); + + // Transpose if a is C-continuous + let mut a_t = None; + let a_layout = match a_layout { + MatrixLayout::C { .. } => { + a_t = Some(vec![Self::zero(); a.len()]); + transpose(a_layout, a, a_t.as_mut().unwrap()) + } + MatrixLayout::F { .. } => a_layout, + }; + + // Transpose if b is C-continuous + let mut b_t = None; + let b_layout = match b_layout { + MatrixLayout::C { .. } => { + b_t = Some(vec![Self::zero(); b.len()]); + transpose(b_layout, b, b_t.as_mut().unwrap()) + } + MatrixLayout::F { .. } => b_layout, + }; + let rcond: Self::Real = -1.; let mut singular_values: Vec = vec![Self::Real::zero(); k as usize]; let mut rank: i32 = 0; @@ -67,11 +88,11 @@ macro_rules! impl_least_squares_real { $gelsd( m, n, - nrhs, - a, - m, - b, - b_lda, + p, + a_t.as_mut().map(|v| v.as_mut_slice()).unwrap_or(a), + a_layout.lda(), + b_t.as_mut().map(|v| v.as_mut_slice()).unwrap_or(b), + b_layout.lda(), &mut singular_values, rcond, &mut rank, @@ -90,11 +111,11 @@ macro_rules! impl_least_squares_real { $gelsd( m, n, - nrhs, - a, - m, - b, - b_lda, + p, + a_t.as_mut().map(|v| v.as_mut_slice()).unwrap_or(a), + a_layout.lda(), + b_t.as_mut().map(|v| v.as_mut_slice()).unwrap_or(b), + b_layout.lda(), &mut singular_values, rcond, &mut rank, @@ -105,6 +126,12 @@ macro_rules! impl_least_squares_real { ); info.as_lapack_result()?; + // Skip a_t -> a transpose because A has been destroyed + // Re-transpose b + if let Some(b_t) = b_t { + transpose(b_layout, &b_t, b); + } + Ok(LeastSquaresOutput { singular_values, rank, diff --git a/ndarray-linalg/src/least_squares.rs b/ndarray-linalg/src/least_squares.rs index 18d2033f..1a6e3a47 100644 --- a/ndarray-linalg/src/least_squares.rs +++ b/ndarray-linalg/src/least_squares.rs @@ -76,6 +76,7 @@ use crate::types::*; /// is a `m x 1` column vector. If `I` is `Ix2`, the RHS is a `n x k` matrix /// (which can be seen as solving `Ax = b` k times for different b) and /// the solution is a `m x k` matrix. +#[derive(Debug, Clone)] pub struct LeastSquaresResult { /// The singular values of the matrix A in `Ax = b` pub singular_values: Array1, diff --git a/ndarray-linalg/tests/least_squares_nrhs.rs b/ndarray-linalg/tests/least_squares_nrhs.rs index 8072c409..95737873 100644 --- a/ndarray-linalg/tests/least_squares_nrhs.rs +++ b/ndarray-linalg/tests/least_squares_nrhs.rs @@ -9,6 +9,7 @@ fn test_exact(a: Array2, b: Array2) { assert_eq!(b.layout().unwrap().size(), (3, 2)); let result = a.least_squares(&b).unwrap(); + dbg!(&result); // unpack result let x: Array2 = result.solution; let residual_l2_square: Array1 = result.residual_sum_of_squares.unwrap(); From 87f0cea7ea4e7f0bdbbd8cd20bb25e2ce73f12f4 Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Sat, 25 Jul 2020 17:15:49 +0900 Subject: [PATCH 07/11] Support complex case --- lax/src/least_squares.rs | 114 ++++++++------------------------------- 1 file changed, 22 insertions(+), 92 deletions(-) diff --git a/lax/src/least_squares.rs b/lax/src/least_squares.rs index 64349051..2c963158 100644 --- a/lax/src/least_squares.rs +++ b/lax/src/least_squares.rs @@ -28,8 +28,15 @@ pub trait LeastSquaresSvdDivideConquer_: Scalar { ) -> Result>; } -macro_rules! impl_least_squares_real { - ($scalar:ty, $gelsd:path) => { +macro_rules! impl_least_squares { + (@real, $scalar:ty, $gelsd:path) => { + impl_least_squares!(@body, $scalar, $gelsd, ); + }; + (@complex, $scalar:ty, $gelsd:path) => { + impl_least_squares!(@body, $scalar, $gelsd, rwork); + }; + + (@body, $scalar:ty, $gelsd:path, $($rwork:ident),*) => { impl LeastSquaresSvdDivideConquer_ for $scalar { unsafe fn least_squares( l: MatrixLayout, @@ -85,6 +92,9 @@ macro_rules! impl_least_squares_real { let mut info = 0; let mut work_size = [Self::zero()]; let mut iwork_size = [0]; + $( + let mut $rwork = [Self::Real::zero()]; + )* $gelsd( m, n, @@ -98,6 +108,7 @@ macro_rules! impl_least_squares_real { &mut rank, &mut work_size, -1, + $(&mut $rwork,)* &mut iwork_size, &mut info, ); @@ -108,6 +119,10 @@ macro_rules! impl_least_squares_real { let mut work = vec![Self::zero(); lwork]; let liwork = iwork_size[0].to_usize().unwrap(); let mut iwork = vec![0; liwork]; + $( + let lrwork = $rwork[0].to_usize().unwrap(); + let mut $rwork = vec![Self::Real::zero(); lrwork]; + )* $gelsd( m, n, @@ -121,6 +136,7 @@ macro_rules! impl_least_squares_real { &mut rank, &mut work, lwork as i32, + $(&mut $rwork,)* &mut iwork, &mut info, ); @@ -141,93 +157,7 @@ macro_rules! impl_least_squares_real { }; } -impl_least_squares_real!(f64, lapack::dgelsd); -impl_least_squares_real!(f32, lapack::sgelsd); - -macro_rules! impl_least_squares { - ($scalar:ty, $gelsd:path) => { - impl LeastSquaresSvdDivideConquer_ for $scalar { - unsafe fn least_squares( - a_layout: MatrixLayout, - a: &mut [Self], - b: &mut [Self], - ) -> Result> { - let (m, n) = a_layout.size(); - if (m as usize) > b.len() || (n as usize) > b.len() { - return Err(Error::InvalidShape); - } - let k = ::std::cmp::min(m, n); - let nrhs = 1; - let ldb = match a_layout { - MatrixLayout::F { .. } => m.max(n), - MatrixLayout::C { .. } => 1, - }; - let rcond: Self::Real = -1.; - let mut singular_values: Vec = vec![Self::Real::zero(); k as usize]; - let mut rank: i32 = 0; - - $gelsd( - a_layout.lapacke_layout(), - m, - n, - nrhs, - a, - a_layout.lda(), - b, - ldb, - &mut singular_values, - rcond, - &mut rank, - ) - .as_lapack_result()?; - - Ok(LeastSquaresOutput { - singular_values, - rank, - }) - } - - unsafe fn least_squares_nrhs( - a_layout: MatrixLayout, - a: &mut [Self], - b_layout: MatrixLayout, - b: &mut [Self], - ) -> Result> { - let (m, n) = a_layout.size(); - if (m as usize) > b.len() - || (n as usize) > b.len() - || a_layout.lapacke_layout() != b_layout.lapacke_layout() - { - return Err(Error::InvalidShape); - } - let k = ::std::cmp::min(m, n); - let nrhs = b_layout.size().1; - let rcond: Self::Real = -1.; - let mut singular_values: Vec = vec![Self::Real::zero(); k as usize]; - let mut rank: i32 = 0; - - $gelsd( - a_layout.lapacke_layout(), - m, - n, - nrhs, - a, - a_layout.lda(), - b, - b_layout.lda(), - &mut singular_values, - rcond, - &mut rank, - ) - .as_lapack_result()?; - Ok(LeastSquaresOutput { - singular_values, - rank, - }) - } - } - }; -} - -impl_least_squares!(c64, lapacke::zgelsd); -impl_least_squares!(c32, lapacke::cgelsd); +impl_least_squares!(@real, f64, lapack::dgelsd); +impl_least_squares!(@real, f32, lapack::sgelsd); +impl_least_squares!(@complex, c64, lapack::zgelsd); +impl_least_squares!(@complex, c32, lapack::cgelsd); From b742f6fe65907cc6af79734457ec3d2b85ce8149 Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Sat, 25 Jul 2020 17:45:35 +0900 Subject: [PATCH 08/11] Fix assert to support under-determined case --- lax/src/least_squares.rs | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/lax/src/least_squares.rs b/lax/src/least_squares.rs index 2c963158..02f682fb 100644 --- a/lax/src/least_squares.rs +++ b/lax/src/least_squares.rs @@ -57,12 +57,12 @@ macro_rules! impl_least_squares { // // where // A : (m, n) - // b : (m, p) - // x : (n, p) + // b : (max(m, n), nrhs) // `b` has to store `x` on exit + // x : (n, nrhs) let (m, n) = a_layout.size(); - let (m_, p) = b_layout.size(); + let (m_, nrhs) = b_layout.size(); let k = m.min(n); - assert_eq!(m, m_); + assert!(m_ >= m); // Transpose if a is C-continuous let mut a_t = None; @@ -98,7 +98,7 @@ macro_rules! impl_least_squares { $gelsd( m, n, - p, + nrhs, a_t.as_mut().map(|v| v.as_mut_slice()).unwrap_or(a), a_layout.lda(), b_t.as_mut().map(|v| v.as_mut_slice()).unwrap_or(b), @@ -126,7 +126,7 @@ macro_rules! impl_least_squares { $gelsd( m, n, - p, + nrhs, a_t.as_mut().map(|v| v.as_mut_slice()).unwrap_or(a), a_layout.lda(), b_t.as_mut().map(|v| v.as_mut_slice()).unwrap_or(b), From 5ad434371cdf7644b9008e51f4a10fac7113f612 Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Sat, 25 Jul 2020 18:26:12 +0900 Subject: [PATCH 09/11] Enable tests for C/F mixed cases #234 --- ndarray-linalg/tests/least_squares_nrhs.rs | 12 ------------ 1 file changed, 12 deletions(-) diff --git a/ndarray-linalg/tests/least_squares_nrhs.rs b/ndarray-linalg/tests/least_squares_nrhs.rs index 95737873..dd7d283c 100644 --- a/ndarray-linalg/tests/least_squares_nrhs.rs +++ b/ndarray-linalg/tests/least_squares_nrhs.rs @@ -37,8 +37,6 @@ macro_rules! impl_exact { test_exact(a, b) } - /* Unsupported currently. See https://github.com/rust-ndarray/ndarray-linalg/issues/234 - #[test] fn []() { let a: Array2<$scalar> = random((3, 3)); @@ -53,8 +51,6 @@ macro_rules! impl_exact { test_exact(a, b) } - */ - #[test] fn []() { let a: Array2<$scalar> = random((3, 3).f()); @@ -109,8 +105,6 @@ macro_rules! impl_overdetermined { test_overdetermined(a, b) } - /* Unsupported currently. See https://github.com/rust-ndarray/ndarray-linalg/issues/234 - #[test] fn []() { let a: Array2<$scalar> = random((4, 3).f()); @@ -125,8 +119,6 @@ macro_rules! impl_overdetermined { test_overdetermined(a, b) } - */ - #[test] fn []() { let a: Array2<$scalar> = random((4, 3).f()); @@ -168,8 +160,6 @@ macro_rules! impl_underdetermined { test_underdetermined(a, b) } - /* Unsupported currently. See https://github.com/rust-ndarray/ndarray-linalg/issues/234 - #[test] fn []() { let a: Array2<$scalar> = random((3, 4).f()); @@ -184,8 +174,6 @@ macro_rules! impl_underdetermined { test_underdetermined(a, b) } - */ - #[test] fn []() { let a: Array2<$scalar> = random((3, 4).f()); From 00b044b7cca26968b4671b978288cb26fcd1ff97 Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Sat, 25 Jul 2020 18:26:42 +0900 Subject: [PATCH 10/11] Fix error handling --- ndarray-linalg/src/least_squares.rs | 31 ++++++++++------------------- 1 file changed, 11 insertions(+), 20 deletions(-) diff --git a/ndarray-linalg/src/least_squares.rs b/ndarray-linalg/src/least_squares.rs index 1a6e3a47..0431a822 100644 --- a/ndarray-linalg/src/least_squares.rs +++ b/ndarray-linalg/src/least_squares.rs @@ -267,6 +267,9 @@ where &mut self, rhs: &mut ArrayBase, ) -> Result> { + if self.shape()[0] != rhs.shape()[0] { + return Err(ShapeError::from_kind(ErrorKind::IncompatibleShape).into()); + } let (m, n) = (self.shape()[0], self.shape()[1]); if n > m { // we need a new rhs b/c it will be overwritten with the solution @@ -285,7 +288,7 @@ fn compute_least_squares_srhs( rhs: &mut ArrayBase, ) -> Result> where - E: Scalar + Lapack + LeastSquaresSvdDivideConquer_, + E: Scalar + Lapack, D1: DataMut, D2: DataMut, { @@ -293,7 +296,7 @@ where singular_values, rank, } = unsafe { - ::least_squares( + E::least_squares( a.layout()?, a.as_allocated_mut()?, rhs.as_slice_memory_order_mut() @@ -348,6 +351,9 @@ where &mut self, rhs: &mut ArrayBase, ) -> Result> { + if self.shape()[0] != rhs.shape()[0] { + return Err(ShapeError::from_kind(ErrorKind::IncompatibleShape).into()); + } let (m, n) = (self.shape()[0], self.shape()[1]); if n > m { // we need a new rhs b/c it will be overwritten with the solution @@ -550,28 +556,13 @@ mod tests { // // Testing error cases // - #[test] fn incompatible_shape_error_on_mismatching_num_rows() { let a: Array2 = array![[1., 2.], [4., 5.], [3., 4.]]; let b: Array1 = array![1., 2.]; - let res = a.least_squares(&b); - match res { - Err(LinalgError::Lapack(err)) if matches!(err, lax::error::Error::InvalidShape) => {} - _ => panic!("Expected Err()"), - } - } - - #[test] - fn incompatible_shape_error_on_mismatching_layout() { - let a: Array2 = array![[1., 2.], [4., 5.], [3., 4.]]; - let b = array![[1.], [2.]].t().to_owned(); - assert_eq!(b.layout().unwrap(), MatrixLayout::F { col: 2, lda: 1 }); - - let res = a.least_squares(&b); - match res { - Err(LinalgError::Lapack(err)) if matches!(err, lax::error::Error::InvalidShape) => {} - _ => panic!("Expected Err()"), + match a.least_squares(&b) { + Err(LinalgError::Shape(e)) if e.kind() == ErrorKind::IncompatibleShape => {} + _ => panic!("Should be raise IncompatibleShape"), } } } From 962497eaecdeff4a5c0a2c5d614534217b037d79 Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Sat, 25 Jul 2020 18:32:13 +0900 Subject: [PATCH 11/11] Drop unsafe --- lax/src/least_squares.rs | 80 +++++++++++++++-------------- ndarray-linalg/src/least_squares.rs | 28 +++++----- 2 files changed, 54 insertions(+), 54 deletions(-) diff --git a/lax/src/least_squares.rs b/lax/src/least_squares.rs index 02f682fb..d684c9b8 100644 --- a/lax/src/least_squares.rs +++ b/lax/src/least_squares.rs @@ -14,13 +14,13 @@ pub struct LeastSquaresOutput { /// Wraps `*gelsd` pub trait LeastSquaresSvdDivideConquer_: Scalar { - unsafe fn least_squares( + fn least_squares( a_layout: MatrixLayout, a: &mut [Self], b: &mut [Self], ) -> Result>; - unsafe fn least_squares_nrhs( + fn least_squares_nrhs( a_layout: MatrixLayout, a: &mut [Self], b_layout: MatrixLayout, @@ -38,7 +38,7 @@ macro_rules! impl_least_squares { (@body, $scalar:ty, $gelsd:path, $($rwork:ident),*) => { impl LeastSquaresSvdDivideConquer_ for $scalar { - unsafe fn least_squares( + fn least_squares( l: MatrixLayout, a: &mut [Self], b: &mut [Self], @@ -47,7 +47,7 @@ macro_rules! impl_least_squares { Self::least_squares_nrhs(l, a, b_layout, b) } - unsafe fn least_squares_nrhs( + fn least_squares_nrhs( a_layout: MatrixLayout, a: &mut [Self], b_layout: MatrixLayout, @@ -95,23 +95,25 @@ macro_rules! impl_least_squares { $( let mut $rwork = [Self::Real::zero()]; )* - $gelsd( - m, - n, - nrhs, - a_t.as_mut().map(|v| v.as_mut_slice()).unwrap_or(a), - a_layout.lda(), - b_t.as_mut().map(|v| v.as_mut_slice()).unwrap_or(b), - b_layout.lda(), - &mut singular_values, - rcond, - &mut rank, - &mut work_size, - -1, - $(&mut $rwork,)* - &mut iwork_size, - &mut info, - ); + unsafe { + $gelsd( + m, + n, + nrhs, + a_t.as_mut().map(|v| v.as_mut_slice()).unwrap_or(a), + a_layout.lda(), + b_t.as_mut().map(|v| v.as_mut_slice()).unwrap_or(b), + b_layout.lda(), + &mut singular_values, + rcond, + &mut rank, + &mut work_size, + -1, + $(&mut $rwork,)* + &mut iwork_size, + &mut info, + ) + }; info.as_lapack_result()?; // calc @@ -123,23 +125,25 @@ macro_rules! impl_least_squares { let lrwork = $rwork[0].to_usize().unwrap(); let mut $rwork = vec![Self::Real::zero(); lrwork]; )* - $gelsd( - m, - n, - nrhs, - a_t.as_mut().map(|v| v.as_mut_slice()).unwrap_or(a), - a_layout.lda(), - b_t.as_mut().map(|v| v.as_mut_slice()).unwrap_or(b), - b_layout.lda(), - &mut singular_values, - rcond, - &mut rank, - &mut work, - lwork as i32, - $(&mut $rwork,)* - &mut iwork, - &mut info, - ); + unsafe { + $gelsd( + m, + n, + nrhs, + a_t.as_mut().map(|v| v.as_mut_slice()).unwrap_or(a), + a_layout.lda(), + b_t.as_mut().map(|v| v.as_mut_slice()).unwrap_or(b), + b_layout.lda(), + &mut singular_values, + rcond, + &mut rank, + &mut work, + lwork as i32, + $(&mut $rwork,)* + &mut iwork, + &mut info, + ); + } info.as_lapack_result()?; // Skip a_t -> a transpose because A has been destroyed diff --git a/ndarray-linalg/src/least_squares.rs b/ndarray-linalg/src/least_squares.rs index 0431a822..0ff518ad 100644 --- a/ndarray-linalg/src/least_squares.rs +++ b/ndarray-linalg/src/least_squares.rs @@ -295,14 +295,12 @@ where let LeastSquaresOutput:: { singular_values, rank, - } = unsafe { - E::least_squares( - a.layout()?, - a.as_allocated_mut()?, - rhs.as_slice_memory_order_mut() - .ok_or_else(|| LinalgError::MemoryNotCont)?, - )? - }; + } = E::least_squares( + a.layout()?, + a.as_allocated_mut()?, + rhs.as_slice_memory_order_mut() + .ok_or_else(|| LinalgError::MemoryNotCont)?, + )?; let (m, n) = (a.shape()[0], a.shape()[1]); let solution = rhs.slice(s![0..n]).to_owned(); @@ -385,14 +383,12 @@ where let LeastSquaresOutput:: { singular_values, rank, - } = unsafe { - E::least_squares_nrhs( - a_layout, - a.as_allocated_mut()?, - rhs_layout, - rhs.as_allocated_mut()?, - )? - }; + } = E::least_squares_nrhs( + a_layout, + a.as_allocated_mut()?, + rhs_layout, + rhs.as_allocated_mut()?, + )?; let solution: Array2 = rhs.slice(s![..a.shape()[1], ..]).to_owned(); let singular_values = Array::from_shape_vec((singular_values.len(),), singular_values)?;