russell_lab 1.13.0

Scientific laboratory for linear algebra and numerical mathematics
Documentation 
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
104
105
106
107
108
109
110
111

use super::ComplexMatrix;
use crate::{to_i32, Complex64, Norm};

extern "C" {
    // Computes the matrix norm
    // <https://www.netlib.org/lapack/explore-html/d5/d8f/zlange_8f.html>
    fn c_zlange(
        norm_code: i32,
        m: *const i32,
        n: *const i32,
        a: *const Complex64,
        lda: *const i32,
        work: *mut f64,
    ) -> f64;
}

/// (zlange) Computes the matrix norm
///
/// Computes one of:
///
/// ```text
/// ‖a‖_1 = max_j ( Σ_i |aij| )
///
/// ‖a‖_∞ = max_i ( Σ_j |aij| )
///
/// ‖a‖_F = sqrt(Σ_i Σ_j |aij|²) == ‖a‖_2
///
/// ‖a‖_max = max_ij ( |aij| )
/// ```
///
/// See also: <https://www.netlib.org/lapack/explore-html/d5/d8f/zlange_8f.html>
///
/// # Examples
///
/// ```
/// use russell_lab::*;
///
/// fn main() {
///     let a = ComplexMatrix::from(&[
///         [-2.0,  2.0],
///         [ 1.0, -4.0],
///     ]);
///     approx_eq(complex_mat_norm(&a, Norm::One), 6.0, 1e-15);
///     approx_eq(complex_mat_norm(&a, Norm::Inf), 5.0, 1e-15);
///     approx_eq(complex_mat_norm(&a, Norm::Fro), 5.0, 1e-15);
///     approx_eq(complex_mat_norm(&a, Norm::Max), 4.0, 1e-15);
/// }
/// ```
pub fn complex_mat_norm(a: &ComplexMatrix, kind: Norm) -> f64 {
    let (m, n) = a.dims();
    if m == 0 || n == 0 {
        return 0.0;
    }

    let m_i32 = to_i32(m);
    let n_i32 = to_i32(n);
    let lda = m_i32;
    let mut work = if kind == Norm::Inf { vec![0.0; m] } else { Vec::new() };
    let norm_code = kind as i32;
    unsafe { c_zlange(norm_code, &m_i32, &n_i32, a.as_data().as_ptr(), &lda, work.as_mut_ptr()) }
}

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

#[cfg(test)]
mod tests {
    use super::{complex_mat_norm, ComplexMatrix};
    use crate::{approx_eq, cpx, Complex64, ComplexFloat, Norm};

    #[test]
    fn complex_mat_norm_works() {
        let a_0x0 = ComplexMatrix::new(0, 0);
        let a_0x1 = ComplexMatrix::new(0, 1);
        let a_1x0 = ComplexMatrix::new(1, 0);
        assert_eq!(complex_mat_norm(&a_0x0, Norm::One), 0.0);
        assert_eq!(complex_mat_norm(&a_0x0, Norm::Inf), 0.0);
        assert_eq!(complex_mat_norm(&a_0x0, Norm::Fro), 0.0);
        assert_eq!(complex_mat_norm(&a_0x0, Norm::Max), 0.0);
        assert_eq!(complex_mat_norm(&a_0x1, Norm::One), 0.0);
        assert_eq!(complex_mat_norm(&a_0x1, Norm::Inf), 0.0);
        assert_eq!(complex_mat_norm(&a_0x1, Norm::Fro), 0.0);
        assert_eq!(complex_mat_norm(&a_0x1, Norm::Max), 0.0);
        assert_eq!(complex_mat_norm(&a_1x0, Norm::One), 0.0);
        assert_eq!(complex_mat_norm(&a_1x0, Norm::Inf), 0.0);
        assert_eq!(complex_mat_norm(&a_1x0, Norm::Fro), 0.0);
        assert_eq!(complex_mat_norm(&a_1x0, Norm::Max), 0.0);
        #[rustfmt::skip]
        let a = ComplexMatrix::from(&[
            [cpx!( 5.0, 1.0), cpx!(-4.0, 2.0), cpx!(2.0, 3.0)],
            [cpx!(-1.0, 1.0), cpx!( 2.0, 2.0), cpx!(3.0, 3.0)],
            [cpx!(-2.0, 1.0), cpx!( 1.0, 2.0), cpx!(0.0, 3.0)],
        ]);
        approx_eq(
            complex_mat_norm(&a, Norm::One),
            a.get(0, 2).abs() + a.get(1, 2).abs() + a.get(2, 2).abs(),
            1e-14,
        );
        approx_eq(
            complex_mat_norm(&a, Norm::Inf),
            a.get(0, 0).abs() + a.get(0, 1).abs() + a.get(0, 2).abs(),
            1e-14,
        );
        let mut fro = 0.0;
        for v in a.as_data() {
            fro += v.abs() * v.abs();
        }
        fro = f64::sqrt(fro);
        approx_eq(complex_mat_norm(&a, Norm::Fro), fro, 1e-14);
        approx_eq(complex_mat_norm(&a, Norm::Max), cpx!(5.0, 1.0).abs(), 1e-14);
    }
}