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<i32>;
#[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<Self::Real>;
+
+ /// 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<Self::Real> {
+ 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<T: Scalar> {
+ pub ty: NormType,
+ pub layout: MatrixLayout,
+ pub work: Vec<MaybeUninit<T::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(),
+pub trait OperatorNormWorkImpl {
+ type Elem: Scalar;
+ fn new(t: NormType, l: MatrixLayout) -> Self;
+ fn calc(&mut self, a: &[Self::Elem]) -> <Self::Elem as Scalar>::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<MaybeUninit<Self::Real>> = if matches!(t, NormType::Infinity) {
- vec_uninit(m as usize)
- } else {
- Vec::new()
+ OperatorNormWork { ty, layout, work }
+ }
+
+ fn calc(&mut self, a: &[Self::Elem]) -> <Self::Elem as Scalar>::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<Self::Real>;
+pub struct RcondWork<T: Scalar> {
+ pub layout: MatrixLayout,
+ pub work: Vec<MaybeUninit<T>>,
+ pub rwork: Option<Vec<MaybeUninit<T::Real>>>,
+ pub iwork: Option<Vec<MaybeUninit<i32>>>,
}
-macro_rules! impl_rcond_real {
- ($scalar:ty, $gecon:path) => {
- impl Rcond_ for $scalar {
- fn rcond(l: MatrixLayout, a: &[Self], anorm: Self::Real) -> Result<Self::Real> {
- 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: <Self::Elem as Scalar>::Real,
+ ) -> Result<<Self::Elem as Scalar>::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<MaybeUninit<Self>> = vec_uninit(4 * n as usize);
- let mut iwork: Vec<MaybeUninit<i32>> = vec_uninit(n as usize);
- let norm_type = match l {
+ fn calc(
+ &mut self,
+ a: &[Self::Elem],
+ anorm: <Self::Elem as Scalar>::Real,
+ ) -> Result<<Self::Elem as Scalar>::Real> {
+ let (n, _) = self.layout.size();
+ let mut rcond = <Self::Elem as Scalar>::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<Self::Real> {
- let (n, _) = l.size();
- let mut rcond = Self::Real::zero();
+ fn calc(
+ &mut self,
+ a: &[Self::Elem],
+ anorm: <Self::Elem as Scalar>::Real,
+ ) -> Result<<Self::Elem as Scalar>::Real> {
+ let (n, _) = self.layout.size();
+ let mut rcond = <Self::Elem as Scalar>::Real::zero();
let mut info = 0;
- let mut work: Vec<MaybeUninit<Self>> = vec_uninit(2 * n as usize);
- let mut rwork: Vec<MaybeUninit<Self::Real>> = 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_);