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
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
use crate::amd::*;
use crate::amd_2::amd_2;
use crate::internal::*;
use crate::valid::valid;
use num_traits::PrimInt;
use std::fmt::Display;
pub fn amd_1<I: PrimInt + Display>(
n: I,
a_p: &[I],
a_i: &[I],
mut len: &mut [usize],
iwlen: usize,
control: &Control,
mut info: &mut Info,
) -> (Vec<I>, Vec<I>) {
let un = n.to_usize().unwrap();
// Construct the matrix for amd_2.
debug_assert!(un > 0);
let mut p_e: Vec<isize> = vec![0; un];
let mut s_p: Vec<usize> = vec![0; un];
let mut t_p: Vec<usize> = vec![0; un];
let mut i_w: Vec<isize> = vec![0; iwlen];
debug_assert!(valid(n, n, a_p, a_i) == Status::OK);
// Construct the pointers for A+A'.
let mut pfree: usize = 0;
for j in 0..un {
p_e[j] = pfree as isize;
s_p[j] = pfree;
pfree += len[j];
}
// Note that this restriction on iwlen is slightly more restrictive than
// what is strictly required in amd_2. amd_2 can operate with no elbow
// room at all, but it will be very slow. For better performance, at
// least size-n elbow room is enforced.
debug_assert!(iwlen >= pfree + un);
#[cfg(feature = "debug1")]
for p in 0..iwlen {
i_w[p] = EMPTY;
}
for k in 0..un {
debug1_print!("Construct row/column k = {} of A+A'\n", k);
let p1 = a_p[k].to_usize().unwrap();
let p2 = a_p[k + 1].to_usize().unwrap();
// Construct A+A'.
let mut p = p1;
while p < p2 {
// Scan the upper triangular part of A.
let j = a_i[p].to_usize().unwrap();
debug_assert!(/*j >= 0 &&*/ j < un);
if j < k {
// Entry A(j,k) in the strictly upper triangular part.
#[cfg(feature = "debug1")]
if j == (un - 1) {
debug_assert!(s_p[j] < pfree);
} else {
debug_assert!((s_p[j] as isize) < p_e[j + 1]);
}
#[cfg(feature = "debug1")]
if k == (un - 1) {
debug_assert!(s_p[k] < pfree);
} else {
debug_assert!((s_p[k] as isize) < p_e[k + 1]);
}
i_w[s_p[j]] = k as isize;
s_p[j] += 1;
i_w[s_p[k]] = j as isize;
s_p[k] += 1;
p += 1;
} else if j == k {
// Skip the diagonal.
p += 1;
break;
} else {
// j > k
// First entry below the diagonal.
break;
}
// Scan lower triangular part of A, in column j until reaching
// row k. Start where last scan left off.
debug_assert!(
a_p[j].to_usize().unwrap() <= t_p[j] && t_p[j] <= a_p[j + 1].to_usize().unwrap()
);
let mut pj = t_p[j];
let pj2 = a_p[j + 1].to_usize().unwrap();
while pj < pj2 {
let i = a_i[pj].to_usize().unwrap();
debug_assert!(/*i >= 0 &&*/ i < un);
if i < k {
// A (i,j) is only in the lower part, not in upper.
#[cfg(feature = "debug1")]
if i == un - 1 {
debug_assert!(s_p[i] < pfree);
} else {
debug_assert!((s_p[i] as isize) < p_e[i + 1]);
}
#[cfg(feature = "debug1")]
if j == un - 1 {
debug_assert!(s_p[j] < pfree);
} else {
debug_assert!((s_p[j] as isize) < p_e[j + 1]);
}
i_w[s_p[i]] = j as isize;
s_p[i] += 1;
i_w[s_p[j]] = i as isize;
s_p[j] += 1;
pj += 1;
} else if i == k {
// Entry A(k,j) in lower part and A(j,k) in upper.
pj += 1;
break;
} else {
// i > k
// Consider this entry later, when k advances to i.
break;
}
}
t_p[j] = pj;
}
t_p[k] = p;
}
// Clean up, for remaining mismatched entries.
for j in 0..un {
let p2 = a_p[j + 1].to_usize().unwrap();
for pj in t_p[j]..p2 {
let i = a_i[pj].to_usize().unwrap();
debug_assert!(/*i >= 0 &&*/ i < un);
// A(i,j) is only in the lower part, not in upper.
#[cfg(feature = "debug1")]
if i == un - 1 {
debug_assert!(s_p[i] < pfree);
} else {
debug_assert!((s_p[i] as isize) < p_e[i + 1]);
}
#[cfg(feature = "debug1")]
if j == un - 1 {
debug_assert!(s_p[j] < pfree);
} else {
debug_assert!((s_p[j] as isize) < p_e[j + 1]);
}
i_w[s_p[i]] = j as isize;
s_p[i] += 1;
i_w[s_p[j]] = i as isize;
s_p[j] += 1;
}
}
#[cfg(feature = "debug1")]
for j in 0..un - 1 {
debug_assert!((s_p[j] as isize) == p_e[j + 1])
}
debug_assert!(s_p[un - 1] == pfree);
// Tp and Sp no longer needed.
// Order the matrix.
let (_nv, p_inv, p, _e_len) = amd_2::<I>(
n, &mut p_e, &mut i_w, &mut len, iwlen, pfree, control, &mut info,
);
(p, p_inv)
}