ndarray_linalg/
opnorm.rs

1//! Operator norm
2
3use lax::Tridiagonal;
4use ndarray::*;
5
6use crate::error::*;
7use crate::layout::*;
8use crate::types::*;
9
10pub use lax::NormType;
11
12/// Operator norm using `*lange` LAPACK routines
13///
14/// [Wikipedia article on operator norm](https://en.wikipedia.org/wiki/Operator_norm)
15pub trait OperationNorm {
16    /// the value of norm
17    type Output: Scalar;
18
19    fn opnorm(&self, t: NormType) -> Result<Self::Output>;
20
21    /// the one norm of a matrix (maximum column sum)
22    fn opnorm_one(&self) -> Result<Self::Output> {
23        self.opnorm(NormType::One)
24    }
25
26    /// the infinity norm of a matrix (maximum row sum)
27    fn opnorm_inf(&self) -> Result<Self::Output> {
28        self.opnorm(NormType::Infinity)
29    }
30
31    /// the Frobenius norm of a matrix (square root of sum of squares)
32    fn opnorm_fro(&self) -> Result<Self::Output> {
33        self.opnorm(NormType::Frobenius)
34    }
35}
36
37impl<A, S> OperationNorm for ArrayBase<S, Ix2>
38where
39    A: Scalar + Lapack,
40    S: Data<Elem = A>,
41{
42    type Output = A::Real;
43
44    fn opnorm(&self, t: NormType) -> Result<Self::Output> {
45        let l = self.layout()?;
46        let a = self.as_allocated()?;
47        Ok(A::opnorm(t, l, a))
48    }
49}
50
51impl<A> OperationNorm for Tridiagonal<A>
52where
53    A: Scalar + Lapack,
54{
55    type Output = A::Real;
56
57    fn opnorm(&self, t: NormType) -> Result<Self::Output> {
58        // `self` is a tridiagonal matrix like,
59        // [d0, u1,  0,   ...,       0,
60        //  l1, d1, u2,            ...,
61        //   0, l2, d2,
62        //  ...           ...,  u{n-1},
63        //   0,  ...,  l{n-1},  d{n-1},]
64        let arr = match t {
65            // opnorm_one() calculates muximum column sum.
66            // Therefore, This part align the columns and make a (3 x n) matrix like,
67            // [ 0, u1, u2, ..., u{n-1},
68            //  d0, d1, d2, ..., d{n-1},
69            //  l1, l2, l3, ...,      0,]
70            NormType::One => {
71                let zl: Array1<A> = Array::zeros(1);
72                let zu: Array1<A> = Array::zeros(1);
73                let dl = concatenate![Axis(0), &self.dl, zl]; // n
74                let du = concatenate![Axis(0), zu, &self.du]; // n
75                stack![Axis(0), du, &self.d, dl] // 3 x n
76            }
77            // opnorm_inf() calculates muximum row sum.
78            // Therefore, This part align the rows and make a (n x 3) matrix like,
79            // [     0,     d0,     u1,
80            //      l1,     d1,     u2,
81            //      l2,     d2,     u3,
82            //     ...,    ...,    ...,
83            //  l{n-1}, d{n-1},      0,]
84            NormType::Infinity => {
85                let zl: Array1<A> = Array::zeros(1);
86                let zu: Array1<A> = Array::zeros(1);
87                let dl = concatenate![Axis(0), zl, &self.dl]; // n
88                let du = concatenate![Axis(0), &self.du, zu]; // n
89                stack![Axis(1), dl, &self.d, du] // n x 3
90            }
91            // opnorm_fro() calculates square root of sum of squares.
92            // Because it is independent of the shape of matrix,
93            // this part make a (1 x (3n-2)) matrix like,
94            // [l1, ..., l{n-1}, d0, ..., d{n-1}, u1, ..., u{n-1}]
95            NormType::Frobenius => {
96                concatenate![Axis(0), &self.dl, &self.d, &self.du].insert_axis(Axis(0))
97            }
98        };
99        let l = arr.layout()?;
100        let a = arr.as_allocated()?;
101        Ok(A::opnorm(t, l, a))
102    }
103}
ndarray_linalg/opnorm.rs

ndarray_linalg/
opnorm.rs
