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
//! 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> OperationNorm for ArrayRef<A, Ix2>
where
A: Scalar + Lapack,
{
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))
}
}