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> OperationNorm for ArrayRef<A, Ix2>
38where
39    A: Scalar + Lapack,
40{
41    type Output = A::Real;
42
43    fn opnorm(&self, t: NormType) -> Result<Self::Output> {
44        let l = self.layout()?;
45        let a = self.as_allocated()?;
46        Ok(A::opnorm(t, l, a))
47    }
48}
49
50impl<A> OperationNorm for Tridiagonal<A>
51where
52    A: Scalar + Lapack,
53{
54    type Output = A::Real;
55
56    fn opnorm(&self, t: NormType) -> Result<Self::Output> {
57        // `self` is a tridiagonal matrix like,
58        // [d0, u1,  0,   ...,       0,
59        //  l1, d1, u2,            ...,
60        //   0, l2, d2,
61        //  ...           ...,  u{n-1},
62        //   0,  ...,  l{n-1},  d{n-1},]
63        let arr = match t {
64            // opnorm_one() calculates muximum column sum.
65            // Therefore, This part align the columns and make a (3 x n) matrix like,
66            // [ 0, u1, u2, ..., u{n-1},
67            //  d0, d1, d2, ..., d{n-1},
68            //  l1, l2, l3, ...,      0,]
69            NormType::One => {
70                let zl: Array1<A> = Array::zeros(1);
71                let zu: Array1<A> = Array::zeros(1);
72                let dl = concatenate![Axis(0), &self.dl, zl]; // n
73                let du = concatenate![Axis(0), zu, &self.du]; // n
74                stack![Axis(0), du, &self.d, dl] // 3 x n
75            }
76            // opnorm_inf() calculates muximum row sum.
77            // Therefore, This part align the rows and make a (n x 3) matrix like,
78            // [     0,     d0,     u1,
79            //      l1,     d1,     u2,
80            //      l2,     d2,     u3,
81            //     ...,    ...,    ...,
82            //  l{n-1}, d{n-1},      0,]
83            NormType::Infinity => {
84                let zl: Array1<A> = Array::zeros(1);
85                let zu: Array1<A> = Array::zeros(1);
86                let dl = concatenate![Axis(0), zl, &self.dl]; // n
87                let du = concatenate![Axis(0), &self.du, zu]; // n
88                stack![Axis(1), dl, &self.d, du] // n x 3
89            }
90            // opnorm_fro() calculates square root of sum of squares.
91            // Because it is independent of the shape of matrix,
92            // this part make a (1 x (3n-2)) matrix like,
93            // [l1, ..., l{n-1}, d0, ..., d{n-1}, u1, ..., u{n-1}]
94            NormType::Frobenius => {
95                concatenate![Axis(0), &self.dl, &self.d, &self.du].insert_axis(Axis(0))
96            }
97        };
98        let l = arr.layout()?;
99        let a = arr.as_allocated()?;
100        Ok(A::opnorm(t, l, a))
101    }
102}
ndarray_linalg/opnorm.rs

ndarray_linalg/
opnorm.rs
