1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103

//! Operator norm

use lax::Tridiagonal;
use ndarray::*;

use crate::error::*;
use crate::layout::*;
use crate::types::*;

pub use lax::NormType;

/// Operator norm using `*lange` LAPACK routines
///
/// [Wikipedia article on operator norm](https://en.wikipedia.org/wiki/Operator_norm)
pub trait OperationNorm {
    /// the value of norm
    type Output: Scalar;

    fn opnorm(&self, t: NormType) -> Result<Self::Output>;

    /// the one norm of a matrix (maximum column sum)
    fn opnorm_one(&self) -> Result<Self::Output> {
        self.opnorm(NormType::One)
    }

    /// the infinity norm of a matrix (maximum row sum)
    fn opnorm_inf(&self) -> Result<Self::Output> {
        self.opnorm(NormType::Infinity)
    }

    /// the Frobenius norm of a matrix (square root of sum of squares)
    fn opnorm_fro(&self) -> Result<Self::Output> {
        self.opnorm(NormType::Frobenius)
    }
}

impl<A, S> OperationNorm for ArrayBase<S, Ix2>
where
    A: Scalar + Lapack,
    S: Data<Elem = A>,
{
    type Output = A::Real;

    fn opnorm(&self, t: NormType) -> Result<Self::Output> {
        let l = self.layout()?;
        let a = self.as_allocated()?;
        Ok(A::opnorm(t, l, a))
    }
}

impl<A> OperationNorm for Tridiagonal<A>
where
    A: Scalar + Lapack,
{
    type Output = A::Real;

    fn opnorm(&self, t: NormType) -> Result<Self::Output> {
        // `self` is a tridiagonal matrix like,
        // [d0, u1,  0,   ...,       0,
        //  l1, d1, u2,            ...,
        //   0, l2, d2,
        //  ...           ...,  u{n-1},
        //   0,  ...,  l{n-1},  d{n-1},]
        let arr = match t {
            // opnorm_one() calculates muximum column sum.
            // Therefore, This part align the columns and make a (3 x n) matrix like,
            // [ 0, u1, u2, ..., u{n-1},
            //  d0, d1, d2, ..., d{n-1},
            //  l1, l2, l3, ...,      0,]
            NormType::One => {
                let zl: Array1<A> = Array::zeros(1);
                let zu: Array1<A> = Array::zeros(1);
                let dl = concatenate![Axis(0), &self.dl, zl]; // n
                let du = concatenate![Axis(0), zu, &self.du]; // n
                stack![Axis(0), du, &self.d, dl] // 3 x n
            }
            // opnorm_inf() calculates muximum row sum.
            // Therefore, This part align the rows and make a (n x 3) matrix like,
            // [     0,     d0,     u1,
            //      l1,     d1,     u2,
            //      l2,     d2,     u3,
            //     ...,    ...,    ...,
            //  l{n-1}, d{n-1},      0,]
            NormType::Infinity => {
                let zl: Array1<A> = Array::zeros(1);
                let zu: Array1<A> = Array::zeros(1);
                let dl = concatenate![Axis(0), zl, &self.dl]; // n
                let du = concatenate![Axis(0), &self.du, zu]; // n
                stack![Axis(1), dl, &self.d, du] // n x 3
            }
            // opnorm_fro() calculates square root of sum of squares.
            // Because it is independent of the shape of matrix,
            // this part make a (1 x (3n-2)) matrix like,
            // [l1, ..., l{n-1}, d0, ..., d{n-1}, u1, ..., u{n-1}]
            NormType::Frobenius => {
                concatenate![Axis(0), &self.dl, &self.d, &self.du].insert_axis(Axis(0))
            }
        };
        let l = arr.layout()?;
        let a = arr.as_allocated()?;
        Ok(A::opnorm(t, l, a))
    }
}