Skip to content

Commit 120fb07

Browse files
termoshtttermoshtt
authored
Merge pull request #347 from rust-ndarray/lax-rcond-opnorm-work
Merge `Rcond_` and `OperatorNorm_` into `Lapack` trait
2 parents acd7858 + 1abe240 commit 120fb07

File tree

3 files changed

+172
-66
lines changed

3 files changed

+172
-66
lines changed

‎lax/src/lib.rs

+57-4
Original file line numberDiff line numberDiff line change
@@ -93,23 +93,22 @@ pub mod eig;
9393
pub mod eigh;
9494
pub mod eigh_generalized;
9595
pub mod least_squares;
96+
pub mod opnorm;
9697
pub mod qr;
98+
pub mod rcond;
9799
pub mod solve;
98100
pub mod solveh;
99101
pub mod svd;
100102
pub mod svddc;
101103

102104
mod alloc;
103-
mod opnorm;
104-
mod rcond;
105105
mod triangular;
106106
mod tridiagonal;
107107

108108
pub use self::cholesky::*;
109109
pub use self::flags::*;
110110
pub use self::least_squares::LeastSquaresOwned;
111111
pub use self::opnorm::*;
112-
pub use self::rcond::*;
113112
pub use self::svd::{SvdOwned, SvdRef};
114113
pub use self::triangular::*;
115114
pub use self::tridiagonal::*;
@@ -122,7 +121,7 @@ pub type Pivot = Vec<i32>;
122121

123122
#[cfg_attr(doc, katexit::katexit)]
124123
/// Trait for primitive types which implements LAPACK subroutines
125-
pub trait Lapack: OperatorNorm_ + Triangular_ + Tridiagonal_ + Rcond_ {
124+
pub trait Lapack: Triangular_ + Tridiagonal_ {
126125
/// Compute right eigenvalue and eigenvectors for a general matrix
127126
fn eig(
128127
calc_v: bool,
@@ -257,6 +256,48 @@ pub trait Lapack: OperatorNorm_ + Triangular_ + Tridiagonal_ + Rcond_ {
257256

258257
/// Solve linear equation $Ax = b$ using $U$ or $L$ calculated by [Lapack::cholesky]
259258
fn solve_cholesky(l: MatrixLayout, uplo: UPLO, a: &[Self], b: &mut [Self]) -> Result<()>;
259+
260+
/// Estimates the the reciprocal of the condition number of the matrix in 1-norm.
261+
///
262+
/// `anorm` should be the 1-norm of the matrix `a`.
263+
fn rcond(l: MatrixLayout, a: &[Self], anorm: Self::Real) -> Result<Self::Real>;
264+
265+
/// Compute norm of matrices
266+
///
267+
/// For a $n \times m$ matrix
268+
/// $$
269+
/// A = \begin{pmatrix}
270+
/// a_{11} & \cdots & a_{1m} \\\\
271+
/// \vdots & \ddots & \vdots \\\\
272+
/// a_{n1} & \cdots & a_{nm}
273+
/// \end{pmatrix}
274+
/// $$
275+
/// LAPACK can compute three types of norms:
276+
///
277+
/// - Operator norm based on 1-norm for its domain linear space:
278+
/// $$
279+
/// \Vert A \Vert_1 = \sup_{\Vert x \Vert_1 = 1} \Vert Ax \Vert_1
280+
/// = \max_{1 \le j \le m } \sum_{i=1}^n |a_{ij}|
281+
/// $$
282+
/// where
283+
/// $\Vert x\Vert_1 = \sum_{j=1}^m |x_j|$
284+
/// is 1-norm for a vector $x$.
285+
///
286+
/// - Operator norm based on $\infty$-norm for its domain linear space:
287+
/// $$
288+
/// \Vert A \Vert_\infty = \sup_{\Vert x \Vert_\infty = 1} \Vert Ax \Vert_\infty
289+
/// = \max_{1 \le i \le n } \sum_{j=1}^m |a_{ij}|
290+
/// $$
291+
/// where
292+
/// $\Vert x\Vert_\infty = \max_{j=1}^m |x_j|$
293+
/// is $\infty$-norm for a vector $x$.
294+
///
295+
/// - Frobenious norm
296+
/// $$
297+
/// \Vert A \Vert_F = \sqrt{\mathrm{Tr} \left(AA^\dagger\right)} = \sqrt{\sum_{i=1}^n \sum_{j=1}^m |a_{ij}|^2}
298+
/// $$
299+
///
300+
fn opnorm(t: NormType, l: MatrixLayout, a: &[Self]) -> Self::Real;
260301
}
261302

262303
macro_rules! impl_lapack {
@@ -418,6 +459,18 @@ macro_rules! impl_lapack {
418459
use cholesky::*;
419460
SolveCholeskyImpl::solve_cholesky(l, uplo, a, b)
420461
}
462+
463+
fn rcond(l: MatrixLayout, a: &[Self], anorm: Self::Real) -> Result<Self::Real> {
464+
use rcond::*;
465+
let mut work = RcondWork::<$s>::new(l);
466+
work.calc(a, anorm)
467+
}
468+
469+
fn opnorm(t: NormType, l: MatrixLayout, a: &[Self]) -> Self::Real {
470+
use opnorm::*;
471+
let mut work = OperatorNormWork::<$s>::new(t, l);
472+
work.calc(a)
473+
}
421474
}
422475
};
423476
}

‎lax/src/opnorm.rs

+37-23
Original file line numberDiff line numberDiff line change
@@ -1,27 +1,42 @@
1-
//! Operator norms of matrices
1+
//! Operator norm
22
33
use super::{AsPtr, NormType};
44
use crate::{layout::MatrixLayout, *};
55
use cauchy::*;
66

7-
pub trait OperatorNorm_: Scalar {
8-
fn opnorm(t: NormType, l: MatrixLayout, a: &[Self]) -> Self::Real;
7+
pub struct OperatorNormWork<T: Scalar> {
8+
pub ty: NormType,
9+
pub layout: MatrixLayout,
10+
pub work: Vec<MaybeUninit<T::Real>>,
911
}
1012

11-
macro_rules! impl_opnorm {
12-
($scalar:ty, $lange:path) => {
13-
impl OperatorNorm_ for $scalar {
14-
fn opnorm(t: NormType, l: MatrixLayout, a: &[Self]) -> Self::Real {
15-
let m = l.lda();
16-
let n = l.len();
17-
let t = match l {
18-
MatrixLayout::F { .. } => t,
19-
MatrixLayout::C { .. } => t.transpose(),
13+
pub trait OperatorNormWorkImpl {
14+
type Elem: Scalar;
15+
fn new(t: NormType, l: MatrixLayout) -> Self;
16+
fn calc(&mut self, a: &[Self::Elem]) -> <Self::Elem as Scalar>::Real;
17+
}
18+
19+
macro_rules! impl_operator_norm {
20+
($s:ty, $lange:path) => {
21+
impl OperatorNormWorkImpl for OperatorNormWork<$s> {
22+
type Elem = $s;
23+
24+
fn new(ty: NormType, layout: MatrixLayout) -> Self {
25+
let m = layout.lda();
26+
let work = match (ty, layout) {
27+
(NormType::Infinity, MatrixLayout::F { .. })
28+
| (NormType::One, MatrixLayout::C { .. }) => vec_uninit(m as usize),
29+
_ => Vec::new(),
2030
};
21-
let mut work: Vec<MaybeUninit<Self::Real>> = if matches!(t, NormType::Infinity) {
22-
vec_uninit(m as usize)
23-
} else {
24-
Vec::new()
31+
OperatorNormWork { ty, layout, work }
32+
}
33+
34+
fn calc(&mut self, a: &[Self::Elem]) -> <Self::Elem as Scalar>::Real {
35+
let m = self.layout.lda();
36+
let n = self.layout.len();
37+
let t = match self.layout {
38+
MatrixLayout::F { .. } => self.ty,
39+
MatrixLayout::C { .. } => self.ty.transpose(),
2540
};
2641
unsafe {
2742
$lange(
@@ -30,15 +45,14 @@ macro_rules! impl_opnorm {
3045
&n,
3146
AsPtr::as_ptr(a),
3247
&m,
33-
AsPtr::as_mut_ptr(&mut work),
48+
AsPtr::as_mut_ptr(&mut self.work),
3449
)
3550
}
3651
}
3752
}
3853
};
39-
} // impl_opnorm!
40-
41-
impl_opnorm!(f64, lapack_sys::dlange_);
42-
impl_opnorm!(f32, lapack_sys::slange_);
43-
impl_opnorm!(c64, lapack_sys::zlange_);
44-
impl_opnorm!(c32, lapack_sys::clange_);
54+
}
55+
impl_operator_norm!(c64, lapack_sys::zlange_);
56+
impl_operator_norm!(c32, lapack_sys::clange_);
57+
impl_operator_norm!(f64, lapack_sys::dlange_);
58+
impl_operator_norm!(f32, lapack_sys::slange_);

‎lax/src/rcond.rs

+78-39
Original file line numberDiff line numberDiff line change
@@ -1,85 +1,124 @@
1+
//! Reciprocal conditional number
2+
13
use crate::{error::*, layout::MatrixLayout, *};
24
use cauchy::*;
35
use num_traits::Zero;
46

5-
pub trait Rcond_: Scalar + Sized {
6-
/// Estimates the the reciprocal of the condition number of the matrix in 1-norm.
7-
///
8-
/// `anorm` should be the 1-norm of the matrix `a`.
9-
fn rcond(l: MatrixLayout, a: &[Self], anorm: Self::Real) -> Result<Self::Real>;
7+
pub 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>>>,
1012
}
1113

12-
macro_rules! impl_rcond_real {
13-
($scalar:ty, $gecon:path) => {
14-
impl Rcond_ for $scalar {
15-
fn rcond(l: MatrixLayout, a: &[Self], anorm: Self::Real) -> Result<Self::Real> {
16-
let (n, _) = l.size();
17-
let mut rcond = Self::Real::zero();
18-
let mut info = 0;
14+
pub 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+
24+
macro_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+
}
1940

20-
let mut work: Vec<MaybeUninit<Self>> = vec_uninit(4 * n as usize);
21-
let mut iwork: Vec<MaybeUninit<i32>> = vec_uninit(n as usize);
22-
let norm_type = match l {
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 {
2350
MatrixLayout::C { .. } => NormType::Infinity,
2451
MatrixLayout::F { .. } => NormType::One,
2552
};
2653
unsafe {
27-
$gecon(
54+
$con(
2855
norm_type.as_ptr(),
2956
&n,
3057
AsPtr::as_ptr(a),
31-
&l.lda(),
58+
&self.layout.lda(),
3259
&anorm,
3360
&mut rcond,
34-
AsPtr::as_mut_ptr(&mut work),
35-
AsPtr::as_mut_ptr(&mut iwork),
61+
AsPtr::as_mut_ptr(&mut self.work),
62+
AsPtr::as_mut_ptr(self.rwork.as_mut().unwrap()),
3663
&mut info,
3764
)
3865
};
3966
info.as_lapack_result()?;
40-
4167
Ok(rcond)
4268
}
4369
}
4470
};
4571
}
72+
impl_rcond_work_c!(c64, lapack_sys::zgecon_);
73+
impl_rcond_work_c!(c32, lapack_sys::cgecon_);
74+
75+
macro_rules! impl_rcond_work_r {
76+
($r:ty, $con:path) => {
77+
impl RcondWorkImpl for RcondWork<$r> {
78+
type Elem = $r;
4679

47-
impl_rcond_real!(f32, lapack_sys::sgecon_);
48-
impl_rcond_real!(f64, lapack_sys::dgecon_);
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+
}
4991

50-
macro_rules! impl_rcond_complex {
51-
($scalar:ty, $gecon:path) => {
52-
impl Rcond_ for $scalar {
53-
fn rcond(l: MatrixLayout, a: &[Self], anorm: Self::Real) -> Result<Self::Real> {
54-
let (n, _) = l.size();
55-
let mut rcond = Self::Real::zero();
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();
5699
let mut info = 0;
57-
let mut work: Vec<MaybeUninit<Self>> = vec_uninit(2 * n as usize);
58-
let mut rwork: Vec<MaybeUninit<Self::Real>> = vec_uninit(2 * n as usize);
59-
let norm_type = match l {
100+
let norm_type = match self.layout {
60101
MatrixLayout::C { .. } => NormType::Infinity,
61102
MatrixLayout::F { .. } => NormType::One,
62103
};
63104
unsafe {
64-
$gecon(
105+
$con(
65106
norm_type.as_ptr(),
66107
&n,
67108
AsPtr::as_ptr(a),
68-
&l.lda(),
109+
&self.layout.lda(),
69110
&anorm,
70111
&mut rcond,
71-
AsPtr::as_mut_ptr(&mut work),
72-
AsPtr::as_mut_ptr(&mut rwork),
112+
AsPtr::as_mut_ptr(&mut self.work),
113+
AsPtr::as_mut_ptr(self.iwork.as_mut().unwrap()),
73114
&mut info,
74115
)
75116
};
76117
info.as_lapack_result()?;
77-
78118
Ok(rcond)
79119
}
80120
}
81121
};
82122
}
83-
84-
impl_rcond_complex!(c32, lapack_sys::cgecon_);
85-
impl_rcond_complex!(c64, lapack_sys::zgecon_);
123+
impl_rcond_work_r!(f64, lapack_sys::dgecon_);
124+
impl_rcond_work_r!(f32, lapack_sys::sgecon_);

0 commit comments

Comments
 (0)