diff --git a/lax/src/lib.rs b/lax/src/lib.rs index e673d261..4989fe15 100644 --- a/lax/src/lib.rs +++ b/lax/src/lib.rs @@ -93,15 +93,15 @@ 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; pub mod solveh; pub mod svd; 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: Triangular_ + Tridiagonal_ { /// Compute right eigenvalue and eigenvectors for a general matrix fn eig( calc_v: bool, @@ -257,6 +256,48 @@ 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; + + /// 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; } macro_rules! impl_lapack { @@ -418,6 +459,18 @@ 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) + } + + 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 60933489..1789f385 100644 --- a/lax/src/opnorm.rs +++ b/lax/src/opnorm.rs @@ -1,27 +1,42 @@ -//! Operator norms of matrices +//! Operator norm use super::{AsPtr, NormType}; use crate::{layout::MatrixLayout, *}; use cauchy::*; -pub trait OperatorNorm_: Scalar { - fn opnorm(t: NormType, l: MatrixLayout, a: &[Self]) -> Self::Real; +pub struct OperatorNormWork { + pub ty: NormType, + pub layout: MatrixLayout, + pub work: Vec>, } -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(), +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(), }; - let mut work: Vec> = if matches!(t, NormType::Infinity) { - vec_uninit(m as usize) - } else { - 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( @@ -30,15 +45,14 @@ macro_rules! impl_opnorm { &n, AsPtr::as_ptr(a), &m, - AsPtr::as_mut_ptr(&mut work), + AsPtr::as_mut_ptr(&mut self.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_); +} +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_); diff --git a/lax/src/rcond.rs b/lax/src/rcond.rs index 6cc72749..4d4a4c92 100644 --- a/lax/src/rcond.rs +++ b/lax/src/rcond.rs @@ -1,85 +1,124 @@ +//! Reciprocal conditional number + use crate::{error::*, layout::MatrixLayout, *}; use cauchy::*; use num_traits::Zero; -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; +pub struct RcondWork { + pub layout: MatrixLayout, + pub work: Vec>, + pub rwork: Option>>, + pub iwork: Option>>, } -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; +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, + } + } - let mut work: Vec> = vec_uninit(4 * n as usize); - let mut iwork: Vec> = vec_uninit(n as usize); - let norm_type = match l { + 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 { - $gecon( + $con( norm_type.as_ptr(), &n, AsPtr::as_ptr(a), - &l.lda(), + &self.layout.lda(), &anorm, &mut rcond, - AsPtr::as_mut_ptr(&mut work), - AsPtr::as_mut_ptr(&mut iwork), + 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; -impl_rcond_real!(f32, lapack_sys::sgecon_); -impl_rcond_real!(f64, lapack_sys::dgecon_); + 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), + } + } -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(); + 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 mut work: Vec> = vec_uninit(2 * n as usize); - let mut rwork: Vec> = vec_uninit(2 * n as usize); - let norm_type = match l { + let norm_type = match self.layout { MatrixLayout::C { .. } => NormType::Infinity, MatrixLayout::F { .. } => NormType::One, }; unsafe { - $gecon( + $con( norm_type.as_ptr(), &n, AsPtr::as_ptr(a), - &l.lda(), + &self.layout.lda(), &anorm, &mut rcond, - AsPtr::as_mut_ptr(&mut work), - AsPtr::as_mut_ptr(&mut rwork), + 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_complex!(c32, lapack_sys::cgecon_); -impl_rcond_complex!(c64, lapack_sys::zgecon_); +impl_rcond_work_r!(f64, lapack_sys::dgecon_); +impl_rcond_work_r!(f32, lapack_sys::sgecon_);