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
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
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145

// Copyright (C) 2016-2018 ERGO-Code
// Copyright (C) 2022-2023 Richard Lincoln
//
// Depth first search in a graph.

use crate::LUInt;

// Compute `reach(i)` in a graph by depth first search.
//
// `begin`, `end`, `index` define the graph. When `end` is not None, then node `j`
// has neighbours `index[begin[j]..end[j]-1]`. When `end` is None, then the
// neighbour list is terminated by a negative index.
//
// On return `xi[newtop..top-1]` hold `reach(i)` in topological order; `newtop` is
// the function return value. Nodes that were already marked are excluded from
// the reach.
//
// `pstack` is size `m` workspace (`m` the number of nodes in the graph); the
// contents of `pstack` is undefined on entry/return.
//
// `marked` is size m array. Node `j` is marked iff `marked[j] == m`.
// On return nodes `xi[newtop..top-1]` are marked.
//
// If node `i` is marked on entry, the function does nothing.
pub(crate) fn dfs(
    i: usize,
    begin: &[LUInt],
    end: Option<&[LUInt]>,
    index: &[LUInt],
    top: usize,
    xi: &mut [LUInt],
    // pstack: &mut [lu_int],
    pstack: &mut [f64],
    marked: &mut [LUInt],
    m: LUInt,
) -> usize {
    if marked[i as usize] == m as LUInt {
        return top;
    }

    if let Some(end) = end {
        dfs_end(i, begin, end, index, top, xi, pstack, marked, m)
    } else {
        dfs_begin(i, begin, index, top, xi, pstack, marked, m)
    }
}

// adapted from T. Davis, CSPARSE
fn dfs_end(
    mut i: usize,
    begin: &[LUInt],
    end: &[LUInt],
    index: &[LUInt],
    mut top: usize,
    xi: &mut [LUInt],
    // pstack: &mut [lu_int],
    pstack: &mut [f64],
    marked: &mut [LUInt],
    m: LUInt,
) -> usize {
    let mut head: LUInt = 0;
    assert_ne!(marked[i as usize], m as LUInt);

    xi[0] = i as LUInt;
    while head >= 0 {
        i = xi[head as usize] as usize;
        if marked[i] != m {
            // node i has not been visited
            marked[i] = m;
            pstack[head as usize] = begin[i] as f64;
        }
        let mut done = 1;
        // continue dfs at node i
        for p in (pstack[head as usize] as LUInt)..end[i] {
            let inext = index[p as usize];
            if marked[inext as usize] == m {
                continue; // skip visited node
            }
            pstack[head as usize] = (p + 1) as f64;
            // xi[++head] = inext; /* start dfs at node inext */
            head += 1;
            xi[head as usize] = inext; // start dfs at node inext
            done = 0;
            break;
        }
        if done != 0 {
            // node i has no unvisited neighbours
            head -= 1;
            // xi[--top] = i;
            top -= 1;
            xi[top] = i as LUInt;
        }
    }
    top
}

// adapted from T. Davis, CSPARSE
fn dfs_begin(
    mut i: usize,
    begin: &[LUInt],
    index: &[LUInt],
    mut top: usize,
    xi: &mut [LUInt],
    // pstack: &mut [lu_int],
    pstack: &mut [f64],
    marked: &mut [LUInt],
    m: LUInt,
) -> usize {
    let mut head: LUInt = 0;
    assert_ne!(marked[i as usize], m as LUInt);

    xi[0] = i as LUInt;
    while head >= 0 {
        i = xi[head as usize] as usize;
        if marked[i] != m {
            // node i has not been visited
            marked[i] = m;
            pstack[head as usize] = begin[i] as f64;
        }
        let mut done = 1;
        // continue dfs at node i
        // for (p = pstack[head]; (inext = index[p]) >= 0; p++)
        let mut p = pstack[head as usize] as LUInt;
        while index[p as usize] >= 0 {
            let inext = index[p as usize];
            if marked[inext as usize] == m {
                p += 1; // TODO: check
                continue; // skip visited node
            }
            pstack[head as usize] = (p + 1) as f64;
            // xi[++head] = inext; /* start dfs at node inext */
            head += 1;
            xi[head as usize] = inext; // start dfs at node inext
            done = 0;
            break;
        }
        if done != 0 {
            // node i has no unvisited neighbours
            head -= 1;
            top -= 1;
            xi[top] = i as LUInt;
        }
    }
    top
}