Merge Rcond_ and OperatorNorm_ into Lapack trait

lax/src/lib.rs

@@ -93,23 +93,22 @@ 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;
 
 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)
+            }
         }
     };
 }

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_);

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_);