blu 0.2.1

LU factorization with dynamic Markowitz search and columnwise threshold pivoting
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

// Copyright (C) 2016-2018 ERGO-Code
// Copyright (C) 2022-2023 Richard Lincoln

use crate::lu::lu::*;

// Computes the 1-norm and infinity-norm of the matrix that was freshly
// factorized. Unit cols inserted by the factorization are handled implicitly.
pub(crate) fn matrix_norm(
    lu: &mut LU,
    b_begin: &[usize],
    b_end: &[usize],
    b_i: &[usize],
    b_x: &[f64],
) {
    let m = lu.m;
    let rank = lu.rank;
    let pivotcol = &pivotcol!(lu);
    let pivotrow = &pivotrow!(lu);
    let rowsum = &mut lu.work1;

    assert_eq!(lu.nupdate.unwrap(), 0);

    for i in 0..m {
        rowsum[i as usize] = 0.0;
    }
    let mut onenorm = 0.0;
    let mut infnorm = 0.0;
    for k in 0..rank {
        let jpivot = pivotcol[k as usize] as usize;
        let mut colsum = 0.0;
        for pos in b_begin[jpivot]..b_end[jpivot] {
            colsum += b_x[pos as usize].abs();
            rowsum[b_i[pos as usize] as usize] += b_x[pos as usize].abs();
        }
        onenorm = f64::max(onenorm, colsum);
    }
    for k in rank..m {
        let ipivot = pivotrow[k as usize] as usize;
        rowsum[ipivot] += 1.0;
        onenorm = f64::max(onenorm, 1.0);
    }
    for i in 0..m {
        infnorm = f64::max(infnorm, rowsum[i as usize]);
    }

    lu.onenorm = onenorm;
    lu.infnorm = infnorm;
}