From 09f92262a66dc0e1dcf2b201d47be3f3fe84c651 Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Sat, 1 Oct 2022 17:21:39 +0900 Subject: [PATCH 1/5] Add RcondWork --- lax/src/rcond.rs | 119 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 119 insertions(+) diff --git a/lax/src/rcond.rs b/lax/src/rcond.rs index 6cc72749..09ce75a8 100644 --- a/lax/src/rcond.rs +++ b/lax/src/rcond.rs @@ -2,6 +2,125 @@ use crate::{error::*, layout::MatrixLayout, *}; use cauchy::*; use num_traits::Zero; +pub struct RcondWork { + pub layout: MatrixLayout, + pub work: Vec>, + pub rwork: Option>>, + pub iwork: Option>>, +} + +pub trait RcondWorkImpl { + type Elem: Scalar; + fn new(l: MatrixLayout) -> Self; + fn calc( + &mut self, + a: &[Self::Elem], + anorm: ::Real, + ) -> Result<::Real>; +} + +macro_rules! impl_rcond_work_c { + ($c:ty, $con:path) => { + impl RcondWorkImpl for RcondWork<$c> { + type Elem = $c; + + fn new(layout: MatrixLayout) -> Self { + let (n, _) = layout.size(); + let work = vec_uninit(2 * n as usize); + let rwork = vec_uninit(2 * n as usize); + RcondWork { + layout, + work, + rwork: Some(rwork), + iwork: None, + } + } + + fn calc( + &mut self, + a: &[Self::Elem], + anorm: ::Real, + ) -> Result<::Real> { + let (n, _) = self.layout.size(); + let mut rcond = ::Real::zero(); + let mut info = 0; + let norm_type = match self.layout { + MatrixLayout::C { .. } => NormType::Infinity, + MatrixLayout::F { .. } => NormType::One, + }; + unsafe { + $con( + norm_type.as_ptr(), + &n, + AsPtr::as_ptr(a), + &self.layout.lda(), + &anorm, + &mut rcond, + AsPtr::as_mut_ptr(&mut self.work), + AsPtr::as_mut_ptr(self.rwork.as_mut().unwrap()), + &mut info, + ) + }; + info.as_lapack_result()?; + Ok(rcond) + } + } + }; +} +impl_rcond_work_c!(c64, lapack_sys::zgecon_); +impl_rcond_work_c!(c32, lapack_sys::cgecon_); + +macro_rules! impl_rcond_work_r { + ($r:ty, $con:path) => { + impl RcondWorkImpl for RcondWork<$r> { + type Elem = $r; + + fn new(layout: MatrixLayout) -> Self { + let (n, _) = layout.size(); + let work = vec_uninit(4 * n as usize); + let iwork = vec_uninit(n as usize); + RcondWork { + layout, + work, + rwork: None, + iwork: Some(iwork), + } + } + + fn calc( + &mut self, + a: &[Self::Elem], + anorm: ::Real, + ) -> Result<::Real> { + let (n, _) = self.layout.size(); + let mut rcond = ::Real::zero(); + let mut info = 0; + let norm_type = match self.layout { + MatrixLayout::C { .. } => NormType::Infinity, + MatrixLayout::F { .. } => NormType::One, + }; + unsafe { + $con( + norm_type.as_ptr(), + &n, + AsPtr::as_ptr(a), + &self.layout.lda(), + &anorm, + &mut rcond, + AsPtr::as_mut_ptr(&mut self.work), + AsPtr::as_mut_ptr(self.iwork.as_mut().unwrap()), + &mut info, + ) + }; + info.as_lapack_result()?; + Ok(rcond) + } + } + }; +} +impl_rcond_work_r!(f64, lapack_sys::dgecon_); +impl_rcond_work_r!(f32, lapack_sys::sgecon_); + pub trait Rcond_: Scalar + Sized { /// Estimates the the reciprocal of the condition number of the matrix in 1-norm. /// From a118bf0b05308438e47887db0b3a482561b6d340 Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Sat, 1 Oct 2022 17:24:27 +0900 Subject: [PATCH 2/5] Merge Rcond_ into Lapack trait --- lax/src/lib.rs | 16 ++++++++-- lax/src/rcond.rs | 82 ------------------------------------------------ 2 files changed, 13 insertions(+), 85 deletions(-) diff --git a/lax/src/lib.rs b/lax/src/lib.rs index e673d261..dbe68578 100644 --- a/lax/src/lib.rs +++ b/lax/src/lib.rs @@ -94,6 +94,7 @@ pub mod eigh; pub mod eigh_generalized; pub mod least_squares; pub mod qr; +pub mod rcond; pub mod solve; pub mod solveh; pub mod svd; @@ -101,7 +102,6 @@ pub mod svddc; mod alloc; mod opnorm; -mod rcond; mod triangular; mod tridiagonal; @@ -109,7 +109,6 @@ pub use self::cholesky::*; pub use self::flags::*; pub use self::least_squares::LeastSquaresOwned; pub use self::opnorm::*; -pub use self::rcond::*; pub use self::svd::{SvdOwned, SvdRef}; pub use self::triangular::*; pub use self::tridiagonal::*; @@ -122,7 +121,7 @@ pub type Pivot = Vec; #[cfg_attr(doc, katexit::katexit)] /// Trait for primitive types which implements LAPACK subroutines -pub trait Lapack: OperatorNorm_ + Triangular_ + Tridiagonal_ + Rcond_ { +pub trait Lapack: OperatorNorm_ + Triangular_ + Tridiagonal_ { /// Compute right eigenvalue and eigenvectors for a general matrix fn eig( calc_v: bool, @@ -257,6 +256,11 @@ pub trait Lapack: OperatorNorm_ + Triangular_ + Tridiagonal_ + Rcond_ { /// Solve linear equation $Ax = b$ using $U$ or $L$ calculated by [Lapack::cholesky] fn solve_cholesky(l: MatrixLayout, uplo: UPLO, a: &[Self], b: &mut [Self]) -> Result<()>; + + /// Estimates the the reciprocal of the condition number of the matrix in 1-norm. + /// + /// `anorm` should be the 1-norm of the matrix `a`. + fn rcond(l: MatrixLayout, a: &[Self], anorm: Self::Real) -> Result; } macro_rules! impl_lapack { @@ -418,6 +422,12 @@ macro_rules! impl_lapack { use cholesky::*; SolveCholeskyImpl::solve_cholesky(l, uplo, a, b) } + + fn rcond(l: MatrixLayout, a: &[Self], anorm: Self::Real) -> Result { + use rcond::*; + let mut work = RcondWork::<$s>::new(l); + work.calc(a, anorm) + } } }; } diff --git a/lax/src/rcond.rs b/lax/src/rcond.rs index 09ce75a8..1cf86286 100644 --- a/lax/src/rcond.rs +++ b/lax/src/rcond.rs @@ -120,85 +120,3 @@ macro_rules! impl_rcond_work_r { } impl_rcond_work_r!(f64, lapack_sys::dgecon_); impl_rcond_work_r!(f32, lapack_sys::sgecon_); - -pub trait Rcond_: Scalar + Sized { - /// Estimates the the reciprocal of the condition number of the matrix in 1-norm. - /// - /// `anorm` should be the 1-norm of the matrix `a`. - fn rcond(l: MatrixLayout, a: &[Self], anorm: Self::Real) -> Result; -} - -macro_rules! impl_rcond_real { - ($scalar:ty, $gecon:path) => { - impl Rcond_ for $scalar { - fn rcond(l: MatrixLayout, a: &[Self], anorm: Self::Real) -> Result { - let (n, _) = l.size(); - let mut rcond = Self::Real::zero(); - let mut info = 0; - - let mut work: Vec> = vec_uninit(4 * n as usize); - let mut iwork: Vec> = vec_uninit(n as usize); - let norm_type = match l { - MatrixLayout::C { .. } => NormType::Infinity, - MatrixLayout::F { .. } => NormType::One, - }; - unsafe { - $gecon( - norm_type.as_ptr(), - &n, - AsPtr::as_ptr(a), - &l.lda(), - &anorm, - &mut rcond, - AsPtr::as_mut_ptr(&mut work), - AsPtr::as_mut_ptr(&mut iwork), - &mut info, - ) - }; - info.as_lapack_result()?; - - Ok(rcond) - } - } - }; -} - -impl_rcond_real!(f32, lapack_sys::sgecon_); -impl_rcond_real!(f64, lapack_sys::dgecon_); - -macro_rules! impl_rcond_complex { - ($scalar:ty, $gecon:path) => { - impl Rcond_ for $scalar { - fn rcond(l: MatrixLayout, a: &[Self], anorm: Self::Real) -> Result { - let (n, _) = l.size(); - let mut rcond = Self::Real::zero(); - let mut info = 0; - let mut work: Vec> = vec_uninit(2 * n as usize); - let mut rwork: Vec> = vec_uninit(2 * n as usize); - let norm_type = match l { - MatrixLayout::C { .. } => NormType::Infinity, - MatrixLayout::F { .. } => NormType::One, - }; - unsafe { - $gecon( - norm_type.as_ptr(), - &n, - AsPtr::as_ptr(a), - &l.lda(), - &anorm, - &mut rcond, - AsPtr::as_mut_ptr(&mut work), - AsPtr::as_mut_ptr(&mut rwork), - &mut info, - ) - }; - info.as_lapack_result()?; - - Ok(rcond) - } - } - }; -} - -impl_rcond_complex!(c32, lapack_sys::cgecon_); -impl_rcond_complex!(c64, lapack_sys::zgecon_); From dfc720af3eb286dedb2c7165f1f02fcba6c724c2 Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Sat, 1 Oct 2022 17:53:46 +0900 Subject: [PATCH 3/5] Add OperatorNormWork --- lax/src/opnorm.rs | 53 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 53 insertions(+) diff --git a/lax/src/opnorm.rs b/lax/src/opnorm.rs index 60933489..8a08c79f 100644 --- a/lax/src/opnorm.rs +++ b/lax/src/opnorm.rs @@ -4,6 +4,59 @@ use super::{AsPtr, NormType}; use crate::{layout::MatrixLayout, *}; use cauchy::*; +pub struct OperatorNormWork { + pub ty: NormType, + pub layout: MatrixLayout, + pub work: Vec>, +} + +pub trait OperatorNormWorkImpl { + type Elem: Scalar; + fn new(t: NormType, l: MatrixLayout) -> Self; + fn calc(&mut self, a: &[Self::Elem]) -> ::Real; +} + +macro_rules! impl_operator_norm { + ($s:ty, $lange:path) => { + impl OperatorNormWorkImpl for OperatorNormWork<$s> { + type Elem = $s; + + fn new(ty: NormType, layout: MatrixLayout) -> Self { + let m = layout.lda(); + let work = match (ty, layout) { + (NormType::Infinity, MatrixLayout::F { .. }) + | (NormType::One, MatrixLayout::C { .. }) => vec_uninit(m as usize), + _ => Vec::new(), + }; + OperatorNormWork { ty, layout, work } + } + + fn calc(&mut self, a: &[Self::Elem]) -> ::Real { + let m = self.layout.lda(); + let n = self.layout.len(); + let t = match self.layout { + MatrixLayout::F { .. } => self.ty, + MatrixLayout::C { .. } => self.ty.transpose(), + }; + unsafe { + $lange( + t.as_ptr(), + &m, + &n, + AsPtr::as_ptr(a), + &m, + AsPtr::as_mut_ptr(&mut self.work), + ) + } + } + } + }; +} +impl_operator_norm!(c64, lapack_sys::zlange_); +impl_operator_norm!(c32, lapack_sys::clange_); +impl_operator_norm!(f64, lapack_sys::dlange_); +impl_operator_norm!(f32, lapack_sys::slange_); + pub trait OperatorNorm_: Scalar { fn opnorm(t: NormType, l: MatrixLayout, a: &[Self]) -> Self::Real; } From 52a504fe1d23ac16136a34d0d7aa2d9ea51ac7df Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Sat, 1 Oct 2022 18:35:48 +0900 Subject: [PATCH 4/5] Merge OperatorNorm_ into Lapack trait --- lax/src/lib.rs | 13 +++++++++++-- lax/src/opnorm.rs | 39 --------------------------------------- 2 files changed, 11 insertions(+), 41 deletions(-) diff --git a/lax/src/lib.rs b/lax/src/lib.rs index dbe68578..2a7c28d2 100644 --- a/lax/src/lib.rs +++ b/lax/src/lib.rs @@ -93,6 +93,7 @@ pub mod eig; pub mod eigh; pub mod eigh_generalized; pub mod least_squares; +pub mod opnorm; pub mod qr; pub mod rcond; pub mod solve; @@ -101,7 +102,6 @@ pub mod svd; pub mod svddc; mod alloc; -mod opnorm; mod triangular; mod tridiagonal; @@ -121,7 +121,7 @@ pub type Pivot = Vec; #[cfg_attr(doc, katexit::katexit)] /// Trait for primitive types which implements LAPACK subroutines -pub trait Lapack: OperatorNorm_ + Triangular_ + Tridiagonal_ { +pub trait Lapack: Triangular_ + Tridiagonal_ { /// Compute right eigenvalue and eigenvectors for a general matrix fn eig( calc_v: bool, @@ -261,6 +261,9 @@ pub trait Lapack: OperatorNorm_ + Triangular_ + Tridiagonal_ { /// /// `anorm` should be the 1-norm of the matrix `a`. fn rcond(l: MatrixLayout, a: &[Self], anorm: Self::Real) -> Result; + + /// Compute operator norm of a matrix + fn opnorm(t: NormType, l: MatrixLayout, a: &[Self]) -> Self::Real; } macro_rules! impl_lapack { @@ -428,6 +431,12 @@ macro_rules! impl_lapack { let mut work = RcondWork::<$s>::new(l); work.calc(a, anorm) } + + fn opnorm(t: NormType, l: MatrixLayout, a: &[Self]) -> Self::Real { + use opnorm::*; + let mut work = OperatorNormWork::<$s>::new(t, l); + work.calc(a) + } } }; } diff --git a/lax/src/opnorm.rs b/lax/src/opnorm.rs index 8a08c79f..0bfc94f2 100644 --- a/lax/src/opnorm.rs +++ b/lax/src/opnorm.rs @@ -56,42 +56,3 @@ impl_operator_norm!(c64, lapack_sys::zlange_); impl_operator_norm!(c32, lapack_sys::clange_); impl_operator_norm!(f64, lapack_sys::dlange_); impl_operator_norm!(f32, lapack_sys::slange_); - -pub trait OperatorNorm_: Scalar { - fn opnorm(t: NormType, l: MatrixLayout, a: &[Self]) -> Self::Real; -} - -macro_rules! impl_opnorm { - ($scalar:ty, $lange:path) => { - impl OperatorNorm_ for $scalar { - fn opnorm(t: NormType, l: MatrixLayout, a: &[Self]) -> Self::Real { - let m = l.lda(); - let n = l.len(); - let t = match l { - MatrixLayout::F { .. } => t, - MatrixLayout::C { .. } => t.transpose(), - }; - let mut work: Vec> = if matches!(t, NormType::Infinity) { - vec_uninit(m as usize) - } else { - Vec::new() - }; - unsafe { - $lange( - t.as_ptr(), - &m, - &n, - AsPtr::as_ptr(a), - &m, - AsPtr::as_mut_ptr(&mut work), - ) - } - } - } - }; -} // impl_opnorm! - -impl_opnorm!(f64, lapack_sys::dlange_); -impl_opnorm!(f32, lapack_sys::slange_); -impl_opnorm!(c64, lapack_sys::zlange_); -impl_opnorm!(c32, lapack_sys::clange_); From 1abe240a6da8d83566a825aaa6c8693ad0ffe725 Mon Sep 17 00:00:00 2001 From: Toshiki Teramura Date: Sun, 2 Oct 2022 21:29:31 +0900 Subject: [PATCH 5/5] Add document for opnorm --- lax/src/lib.rs | 36 +++++++++++++++++++++++++++++++++++- lax/src/opnorm.rs | 2 +- lax/src/rcond.rs | 2 ++ 3 files changed, 38 insertions(+), 2 deletions(-) diff --git a/lax/src/lib.rs b/lax/src/lib.rs index 2a7c28d2..4989fe15 100644 --- a/lax/src/lib.rs +++ b/lax/src/lib.rs @@ -262,7 +262,41 @@ pub trait Lapack: Triangular_ + Tridiagonal_ { /// `anorm` should be the 1-norm of the matrix `a`. fn rcond(l: MatrixLayout, a: &[Self], anorm: Self::Real) -> Result; - /// Compute operator norm of a matrix + /// Compute norm of matrices + /// + /// For a $n \times m$ matrix + /// $$ + /// A = \begin{pmatrix} + /// a_{11} & \cdots & a_{1m} \\\\ + /// \vdots & \ddots & \vdots \\\\ + /// a_{n1} & \cdots & a_{nm} + /// \end{pmatrix} + /// $$ + /// LAPACK can compute three types of norms: + /// + /// - Operator norm based on 1-norm for its domain linear space: + /// $$ + /// \Vert A \Vert_1 = \sup_{\Vert x \Vert_1 = 1} \Vert Ax \Vert_1 + /// = \max_{1 \le j \le m } \sum_{i=1}^n |a_{ij}| + /// $$ + /// where + /// $\Vert x\Vert_1 = \sum_{j=1}^m |x_j|$ + /// is 1-norm for a vector $x$. + /// + /// - Operator norm based on $\infty$-norm for its domain linear space: + /// $$ + /// \Vert A \Vert_\infty = \sup_{\Vert x \Vert_\infty = 1} \Vert Ax \Vert_\infty + /// = \max_{1 \le i \le n } \sum_{j=1}^m |a_{ij}| + /// $$ + /// where + /// $\Vert x\Vert_\infty = \max_{j=1}^m |x_j|$ + /// is $\infty$-norm for a vector $x$. + /// + /// - Frobenious norm + /// $$ + /// \Vert A \Vert_F = \sqrt{\mathrm{Tr} \left(AA^\dagger\right)} = \sqrt{\sum_{i=1}^n \sum_{j=1}^m |a_{ij}|^2} + /// $$ + /// fn opnorm(t: NormType, l: MatrixLayout, a: &[Self]) -> Self::Real; } diff --git a/lax/src/opnorm.rs b/lax/src/opnorm.rs index 0bfc94f2..1789f385 100644 --- a/lax/src/opnorm.rs +++ b/lax/src/opnorm.rs @@ -1,4 +1,4 @@ -//! Operator norms of matrices +//! Operator norm use super::{AsPtr, NormType}; use crate::{layout::MatrixLayout, *}; diff --git a/lax/src/rcond.rs b/lax/src/rcond.rs index 1cf86286..4d4a4c92 100644 --- a/lax/src/rcond.rs +++ b/lax/src/rcond.rs @@ -1,3 +1,5 @@ +//! Reciprocal conditional number + use crate::{error::*, layout::MatrixLayout, *}; use cauchy::*; use num_traits::Zero;