Merge Rcond_ into Lapack trait

Author: termoshtt

lax/src/lib.rs

@@ -94,22 +94,21 @@ 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;
 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: 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<Self::Real>;
 }
 
 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<Self::Real> {
+                use rcond::*;
+                let mut work = RcondWork::<$s>::new(l);
+                work.calc(a, anorm)
+            }
         }
     };
 }

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<Self::Real>;
-}
-
-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;
-
-                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 {
-                    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<Self::Real> {
-                let (n, _) = l.size();
-                let mut rcond = Self::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 {
-                    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_);