lax/
rcond.rs

1//! Reciprocal conditional number
2
3use crate::{error::*, layout::MatrixLayout, *};
4use cauchy::*;
5use num_traits::Zero;
6
7pub struct RcondWork<T: Scalar> {
8    pub layout: MatrixLayout,
9    pub work: Vec<MaybeUninit<T>>,
10    pub rwork: Option<Vec<MaybeUninit<T::Real>>>,
11    pub iwork: Option<Vec<MaybeUninit<i32>>>,
12}
13
14pub trait RcondWorkImpl {
15    type Elem: Scalar;
16    fn new(l: MatrixLayout) -> Self;
17    fn calc(
18        &mut self,
19        a: &[Self::Elem],
20        anorm: <Self::Elem as Scalar>::Real,
21    ) -> Result<<Self::Elem as Scalar>::Real>;
22}
23
24macro_rules! impl_rcond_work_c {
25    ($c:ty, $con:path) => {
26        impl RcondWorkImpl for RcondWork<$c> {
27            type Elem = $c;
28
29            fn new(layout: MatrixLayout) -> Self {
30                let (n, _) = layout.size();
31                let work = vec_uninit(2 * n as usize);
32                let rwork = vec_uninit(2 * n as usize);
33                RcondWork {
34                    layout,
35                    work,
36                    rwork: Some(rwork),
37                    iwork: None,
38                }
39            }
40
41            fn calc(
42                &mut self,
43                a: &[Self::Elem],
44                anorm: <Self::Elem as Scalar>::Real,
45            ) -> Result<<Self::Elem as Scalar>::Real> {
46                let (n, _) = self.layout.size();
47                let mut rcond = <Self::Elem as Scalar>::Real::zero();
48                let mut info = 0;
49                let norm_type = match self.layout {
50                    MatrixLayout::C { .. } => NormType::Infinity,
51                    MatrixLayout::F { .. } => NormType::One,
52                };
53                unsafe {
54                    $con(
55                        norm_type.as_ptr(),
56                        &n,
57                        AsPtr::as_ptr(a),
58                        &self.layout.lda(),
59                        &anorm,
60                        &mut rcond,
61                        AsPtr::as_mut_ptr(&mut self.work),
62                        AsPtr::as_mut_ptr(self.rwork.as_mut().unwrap()),
63                        &mut info,
64                    )
65                };
66                info.as_lapack_result()?;
67                Ok(rcond)
68            }
69        }
70    };
71}
72impl_rcond_work_c!(c64, lapack_sys::zgecon_);
73impl_rcond_work_c!(c32, lapack_sys::cgecon_);
74
75macro_rules! impl_rcond_work_r {
76    ($r:ty, $con:path) => {
77        impl RcondWorkImpl for RcondWork<$r> {
78            type Elem = $r;
79
80            fn new(layout: MatrixLayout) -> Self {
81                let (n, _) = layout.size();
82                let work = vec_uninit(4 * n as usize);
83                let iwork = vec_uninit(n as usize);
84                RcondWork {
85                    layout,
86                    work,
87                    rwork: None,
88                    iwork: Some(iwork),
89                }
90            }
91
92            fn calc(
93                &mut self,
94                a: &[Self::Elem],
95                anorm: <Self::Elem as Scalar>::Real,
96            ) -> Result<<Self::Elem as Scalar>::Real> {
97                let (n, _) = self.layout.size();
98                let mut rcond = <Self::Elem as Scalar>::Real::zero();
99                let mut info = 0;
100                let norm_type = match self.layout {
101                    MatrixLayout::C { .. } => NormType::Infinity,
102                    MatrixLayout::F { .. } => NormType::One,
103                };
104                unsafe {
105                    $con(
106                        norm_type.as_ptr(),
107                        &n,
108                        AsPtr::as_ptr(a),
109                        &self.layout.lda(),
110                        &anorm,
111                        &mut rcond,
112                        AsPtr::as_mut_ptr(&mut self.work),
113                        AsPtr::as_mut_ptr(self.iwork.as_mut().unwrap()),
114                        &mut info,
115                    )
116                };
117                info.as_lapack_result()?;
118                Ok(rcond)
119            }
120        }
121    };
122}
123impl_rcond_work_r!(f64, lapack_sys::dgecon_);
124impl_rcond_work_r!(f32, lapack_sys::sgecon_);
lax/rcond.rs

lax/
rcond.rs
