Add .rcond() and .rcond_into()

This commit also changes the character representations of `NormType`
to capital letters to be compatible with `*gecon`.

Author: jturner314

src/lapack_traits/opnorm.rs

@@ -8,9 +8,9 @@ use types::*;
 
 #[repr(u8)]
 pub enum NormType {
-    One = b'o',
-    Infinity = b'i',
-    Frobenius = b'f',
+    One = b'O',
+    Infinity = b'I',
+    Frobenius = b'F',
 }
 
 impl NormType {

src/lapack_traits/solve.rs

@@ -4,12 +4,14 @@ use lapack::c;
 
 use error::*;
 use layout::MatrixLayout;
+use num_traits::Zero;
 use types::*;
 
 use super::{Pivot, Transpose, into_result};
+use super::opnorm::NormType;
 
 /// Wraps `*getrf`, `*getri`, and `*getrs`
-pub trait Solve_: Sized {
+pub trait Solve_: AssociatedReal + Sized {
     /// Computes the LU factorization of a general `m x n` matrix `a` using
     /// partial pivoting with row interchanges.
     ///
@@ -20,11 +22,12 @@ pub trait Solve_: Sized {
     /// if it is used to solve a system of equations.
     unsafe fn lu(MatrixLayout, a: &mut [Self]) -> Result<Pivot>;
     unsafe fn inv(MatrixLayout, a: &mut [Self], &Pivot) -> Result<()>;
+    unsafe fn rcond(MatrixLayout, a: &[Self], anorm: Self::Real) -> Result<Self::Real>;
     unsafe fn solve(MatrixLayout, Transpose, a: &[Self], &Pivot, b: &mut [Self]) -> Result<()>;
 }
 
 macro_rules! impl_solve {
-    ($scalar:ty, $getrf:path, $getri:path, $getrs:path) => {
+    ($scalar:ty, $getrf:path, $getri:path, $gecon:path, $getrs:path) => {
 
 impl Solve_ for $scalar {
     unsafe fn lu(l: MatrixLayout, a: &mut [Self]) -> Result<Pivot> {
@@ -41,6 +44,13 @@ impl Solve_ for $scalar {
         into_result(info, ())
     }
 
+    unsafe fn rcond(l: MatrixLayout, a: &[Self], anorm: Self::Real) -> Result<Self::Real> {
+        let (n, _) = l.size();
+        let mut rcond = Self::Real::zero();
+        let info = $gecon(l.lapacke_layout(), NormType::One as u8, n, a, l.lda(), anorm, &mut rcond);
+        into_result(info, rcond)
+    }
+
     unsafe fn solve(l: MatrixLayout, t: Transpose, a: &[Self], ipiv: &Pivot, b: &mut [Self]) -> Result<()> {
         let (n, _) = l.size();
         let nrhs = 1;
@@ -52,7 +62,7 @@ impl Solve_ for $scalar {
 
 }} // impl_solve!
 
-impl_solve!(f64, c::dgetrf, c::dgetri, c::dgetrs);
-impl_solve!(f32, c::sgetrf, c::sgetri, c::sgetrs);
-impl_solve!(c64, c::zgetrf, c::zgetri, c::zgetrs);
-impl_solve!(c32, c::cgetrf, c::cgetri, c::cgetrs);
+impl_solve!(f64, c::dgetrf, c::dgetri, c::dgecon, c::dgetrs);
+impl_solve!(f32, c::sgetrf, c::sgetri, c::sgecon, c::sgetrs);
+impl_solve!(c64, c::zgetrf, c::zgetri, c::zgecon, c::zgetrs);
+impl_solve!(c32, c::cgetrf, c::cgetri, c::cgecon, c::cgetrs);

src/solve.rs

@@ -49,7 +49,9 @@ use ndarray::*;
 
 use super::convert::*;
 use super::error::*;
+use super::lapack_traits::opnorm::NormType;
 use super::layout::*;
+use super::opnorm::OperationNorm;
 use super::types::*;
 
 pub use lapack_traits::{Pivot, Transpose};
@@ -419,3 +421,77 @@ where
         }
     }
 }
+
+/// An interface for *estimating* the reciprocal condition number of matrix refs.
+pub trait ReciprocalConditionNum<A: Scalar> {
+    /// *Estimates* the reciprocal of the condition number of the matrix in
+    /// 1-norm.
+    ///
+    /// This method uses the LAPACK `*gecon` routines, which *estimate*
+    /// `self.inv().opnorm(One)` and then compute `rcond = 1. /
+    /// (self.opnorm(One) * self.inv().opnorm(One))`.
+    ///
+    /// * If `rcond` is near `0.`, the matrix is badly conditioned.
+    /// * If `rcond` is near `1.`, the matrix is well conditioned.
+    fn rcond(&self) -> Result<A::Real>;
+}
+
+/// An interface for *estimating* the reciprocal condition number of matrices.
+pub trait ReciprocalConditionNumInto<A: Scalar> {
+    /// *Estimates* the reciprocal of the condition number of the matrix in
+    /// 1-norm.
+    ///
+    /// This method uses the LAPACK `*gecon` routines, which *estimate*
+    /// `self.inv().opnorm(One)` and then compute `rcond = 1. /
+    /// (self.opnorm(One) * self.inv().opnorm(One))`.
+    ///
+    /// * If `rcond` is near `0.`, the matrix is badly conditioned.
+    /// * If `rcond` is near `1.`, the matrix is well conditioned.
+    fn rcond_into(self) -> Result<A::Real>;
+}
+
+impl<A, S> ReciprocalConditionNum<A> for LUFactorized<S>
+where
+    A: Scalar,
+    S: Data<Elem = A>,
+{
+    fn rcond(&self) -> Result<A::Real> {
+        unsafe {
+            A::rcond(
+                self.a.layout()?,
+                self.a.as_allocated()?,
+                self.a.opnorm(NormType::One)?,
+            )
+        }
+    }
+}
+
+impl<A, S> ReciprocalConditionNumInto<A> for LUFactorized<S>
+where
+    A: Scalar,
+    S: Data<Elem = A>,
+{
+    fn rcond_into(self) -> Result<A::Real> {
+        self.rcond()
+    }
+}
+
+impl<A, S> ReciprocalConditionNum<A> for ArrayBase<S, Ix2>
+where
+    A: Scalar,
+    S: Data<Elem = A>,
+{
+    fn rcond(&self) -> Result<A::Real> {
+        self.factorize()?.rcond_into()
+    }
+}
+
+impl<A, S> ReciprocalConditionNumInto<A> for ArrayBase<S, Ix2>
+where
+    A: Scalar,
+    S: DataMut<Elem = A>,
+{
+    fn rcond_into(self) -> Result<A::Real> {
+        self.factorize_into()?.rcond_into()
+    }
+}

tests/solve.rs

@@ -5,6 +5,7 @@ extern crate num_traits;
 
 use ndarray::*;
 use ndarray_linalg::*;
+use ndarray_linalg::lapack_traits::opnorm::NormType;
 
 #[test]
 fn solve_random() {
@@ -23,3 +24,51 @@ fn solve_random_t() {
     let y = a.solve_into(b).unwrap();
     assert_close_l2!(&x, &y, 1e-7);
 }
+
+#[test]
+fn rcond() {
+    macro_rules! rcond {
+        ($elem:ty, $rows:expr, $atol:expr) => {
+            let a: Array2<$elem> = random(($rows, $rows));
+            let rcond = 1. / (a.opnorm(NormType::One).unwrap() * a.inv().unwrap().opnorm(NormType::One).unwrap());
+            assert_aclose!(a.rcond().unwrap(), rcond, $atol);
+            assert_aclose!(a.rcond_into().unwrap(), rcond, $atol);
+        }
+    }
+    for rows in 1..6 {
+        rcond!(f64, rows, 0.2);
+        rcond!(f32, rows, 0.5);
+        rcond!(c64, rows, 0.2);
+        rcond!(c32, rows, 0.5);
+    }
+}
+
+#[test]
+fn rcond_hilbert() {
+    macro_rules! rcond_hilbert {
+        ($elem:ty, $rows:expr, $atol:expr) => {
+            let a = Array2::<$elem>::from_shape_fn(($rows, $rows), |(i, j)| 1. / (i as $elem + j as $elem - 1.));
+            assert_aclose!(a.rcond().unwrap(), 0., $atol);
+            assert_aclose!(a.rcond_into().unwrap(), 0., $atol);
+        }
+    }
+    rcond_hilbert!(f64, 10, 1e-9);
+    rcond_hilbert!(f32, 10, 1e-3);
+}
+
+#[test]
+fn rcond_identity() {
+    macro_rules! rcond_identity {
+        ($elem:ty, $rows:expr, $atol:expr) => {
+            let a = Array2::<$elem>::eye($rows);
+            assert_aclose!(a.rcond().unwrap(), 1., $atol);
+            assert_aclose!(a.rcond_into().unwrap(), 1., $atol);
+        }
+    }
+    for rows in 1..6 {
+        rcond_identity!(f64, rows, 1e-9);
+        rcond_identity!(f32, rows, 1e-3);
+        rcond_identity!(c64, rows, 1e-9);
+        rcond_identity!(c32, rows, 1e-3);
+    }
+}