conspire/math/sparse/matrix/amd/
mod.rs

1#[cfg(test)]
2mod test;
3
4use super::CscMatrix;
5
6const NONE: isize = -1;
7
8fn flip(i: isize) -> isize {
9    -i - 2
10}
11
12fn wclear(mark: isize, lemax: isize, w: &mut [isize], n: usize) -> isize {
13    if mark < 2 || mark + lemax < 0 {
14        w[..n].iter_mut().for_each(|w_i| {
15            if *w_i != 0 {
16                *w_i = 1
17            }
18        });
19        2
20    } else {
21        mark
22    }
23}
24
25fn tdfs(
26    root: isize,
27    mut count: usize,
28    head: &mut [isize],
29    next: &[isize],
30    post: &mut [isize],
31    stack: &mut [isize],
32) -> usize {
33    let mut top = 0_isize;
34    stack[0] = root;
35    while top >= 0 {
36        let p = stack[top as usize];
37        let i = head[p as usize];
38        if i == NONE {
39            top -= 1;
40            post[count] = p;
41            count += 1;
42        } else {
43            head[p as usize] = next[i as usize];
44            top += 1;
45            stack[top as usize] = i;
46        }
47    }
48    count
49}
50
51impl CscMatrix {
52    /// Fill-reducing ordering from approximate minimum degree on the pattern of A + Aᵀ.
53    pub fn amd(&self) -> Vec<usize> {
54        let n = self.height();
55        assert_eq!(n, self.width());
56        if n == 0 {
57            return Vec::new();
58        }
59        let transpose = self.transpose();
60        let mut cp = vec![0_isize; n + 1];
61        let mut ci = Vec::<isize>::new();
62        let mut len = vec![0_isize; n + 1];
63        (0..n).for_each(|j| {
64            cp[j] = ci.len() as isize;
65            let mut a = self.column(j).map(|(i, _)| i).peekable();
66            let mut b = transpose.column(j).map(|(i, _)| i).peekable();
67            loop {
68                let i = match (a.peek(), b.peek()) {
69                    (Some(&x), Some(&y)) => {
70                        if x <= y {
71                            if x == y {
72                                b.next();
73                            }
74                            a.next();
75                            x
76                        } else {
77                            b.next();
78                            y
79                        }
80                    }
81                    (Some(&x), None) => {
82                        a.next();
83                        x
84                    }
85                    (None, Some(&y)) => {
86                        b.next();
87                        y
88                    }
89                    (None, None) => break,
90                };
91                if i != j {
92                    ci.push(i as isize);
93                }
94            }
95            len[j] = ci.len() as isize - cp[j];
96        });
97        let mut cnz = ci.len() as isize;
98        let nzmax = cnz + cnz / 5 + 2 * n as isize;
99        ci.resize(nzmax as usize, 0);
100        let dense = 16
101            .max((10.0 * (n as f64).sqrt()) as isize)
102            .min(n as isize - 2);
103        let mut nv = vec![1_isize; n + 1];
104        let mut w = vec![1_isize; n + 1];
105        let mut elen = vec![0_isize; n + 1];
106        let mut degree = len.clone();
107        let mut head = vec![NONE; n + 1];
108        let mut next = vec![NONE; n + 1];
109        let mut last = vec![NONE; n + 1];
110        let mut hhead = vec![NONE; n + 1];
111        cp[n] = NONE;
112        elen[n] = -2;
113        w[n] = 0;
114        let mut nel = 0_isize;
115        let mut mindeg = 0_usize;
116        let mut lemax = 0_isize;
117        let mut mark = wclear(0, 0, &mut w, n);
118        (0..n).for_each(|i| {
119            let d = degree[i];
120            if d == 0 {
121                elen[i] = -2;
122                nel += 1;
123                cp[i] = NONE;
124                w[i] = 0;
125            } else if d > dense {
126                nv[i] = 0;
127                elen[i] = NONE;
128                nel += 1;
129                cp[i] = flip(n as isize);
130                nv[n] += 1;
131            } else {
132                if head[d as usize] != NONE {
133                    last[head[d as usize] as usize] = i as isize;
134                }
135                next[i] = head[d as usize];
136                head[d as usize] = i as isize;
137            }
138        });
139        while nel < n as isize {
140            let mut pivot = NONE;
141            while mindeg < n {
142                pivot = head[mindeg];
143                if pivot != NONE {
144                    break;
145                }
146                mindeg += 1;
147            }
148            let k = pivot as usize;
149            if next[k] != NONE {
150                last[next[k] as usize] = NONE;
151            }
152            head[mindeg] = next[k];
153            let elenk = elen[k];
154            let mut nvk = nv[k];
155            nel += nvk;
156            if elenk > 0 && cnz + mindeg as isize >= nzmax {
157                (0..n).for_each(|j| {
158                    let p = cp[j];
159                    if p >= 0 {
160                        cp[j] = ci[p as usize];
161                        ci[p as usize] = flip(j as isize);
162                    }
163                });
164                let mut q = 0_usize;
165                let mut p = 0_usize;
166                while (p as isize) < cnz {
167                    let j = flip(ci[p]);
168                    p += 1;
169                    if j >= 0 {
170                        let j = j as usize;
171                        ci[q] = cp[j];
172                        cp[j] = q as isize;
173                        q += 1;
174                        (0..len[j] - 1).for_each(|_| {
175                            ci[q] = ci[p];
176                            q += 1;
177                            p += 1;
178                        });
179                    }
180                }
181                cnz = q as isize;
182            }
183            let mut dk = 0_isize;
184            nv[k] = -nvk;
185            let mut p = cp[k];
186            let pk1 = if elenk == 0 { p } else { cnz };
187            let mut pk2 = pk1;
188            for k1 in 1..=(elenk + 1) {
189                let (e, mut pj, ln) = if k1 > elenk {
190                    (k, p, len[k] - elenk)
191                } else {
192                    let e = ci[p as usize] as usize;
193                    p += 1;
194                    (e, cp[e], len[e])
195                };
196                (0..ln).for_each(|_| {
197                    let i = ci[pj as usize] as usize;
198                    pj += 1;
199                    let nvi = nv[i];
200                    if nvi > 0 {
201                        dk += nvi;
202                        nv[i] = -nvi;
203                        ci[pk2 as usize] = i as isize;
204                        pk2 += 1;
205                        if next[i] != NONE {
206                            last[next[i] as usize] = last[i];
207                        }
208                        if last[i] != NONE {
209                            next[last[i] as usize] = next[i];
210                        } else {
211                            head[degree[i] as usize] = next[i];
212                        }
213                    }
214                });
215                if e != k {
216                    cp[e] = flip(k as isize);
217                    w[e] = 0;
218                }
219            }
220            if elenk != 0 {
221                cnz = pk2;
222            }
223            degree[k] = dk;
224            cp[k] = pk1;
225            len[k] = pk2 - pk1;
226            elen[k] = -2;
227            mark = wclear(mark, lemax, &mut w, n);
228            (pk1..pk2).for_each(|pk| {
229                let i = ci[pk as usize] as usize;
230                let eln = elen[i];
231                if eln > 0 {
232                    let nvi = -nv[i];
233                    let wnvi = mark - nvi;
234                    (cp[i]..cp[i] + eln).for_each(|p| {
235                        let e = ci[p as usize] as usize;
236                        if w[e] >= mark {
237                            w[e] -= nvi;
238                        } else if w[e] != 0 {
239                            w[e] = degree[e] + wnvi;
240                        }
241                    });
242                }
243            });
244            (pk1..pk2).for_each(|pk| {
245                let i = ci[pk as usize] as usize;
246                let p1 = cp[i];
247                let p2 = p1 + elen[i];
248                let mut pn = p1;
249                let mut hash = 0_usize;
250                let mut d = 0_isize;
251                (p1..p2).for_each(|p| {
252                    let e = ci[p as usize] as usize;
253                    if w[e] != 0 {
254                        let dext = w[e] - mark;
255                        if dext > 0 {
256                            d += dext;
257                            ci[pn as usize] = e as isize;
258                            pn += 1;
259                            hash += e;
260                        } else {
261                            cp[e] = flip(k as isize);
262                            w[e] = 0;
263                        }
264                    }
265                });
266                elen[i] = pn - p1 + 1;
267                let p3 = pn;
268                let p4 = p1 + len[i];
269                (p2..p4).for_each(|p| {
270                    let j = ci[p as usize] as usize;
271                    let nvj = nv[j];
272                    if nvj > 0 {
273                        d += nvj;
274                        ci[pn as usize] = j as isize;
275                        pn += 1;
276                        hash += j;
277                    }
278                });
279                if d == 0 {
280                    cp[i] = flip(k as isize);
281                    let nvi = -nv[i];
282                    dk -= nvi;
283                    nvk += nvi;
284                    nel += nvi;
285                    nv[i] = 0;
286                    elen[i] = NONE;
287                } else {
288                    degree[i] = degree[i].min(d);
289                    ci[pn as usize] = ci[p3 as usize];
290                    ci[p3 as usize] = ci[p1 as usize];
291                    ci[p1 as usize] = k as isize;
292                    len[i] = pn - p1 + 1;
293                    let hash = hash % n;
294                    next[i] = hhead[hash];
295                    hhead[hash] = i as isize;
296                    last[i] = hash as isize;
297                }
298            });
299            degree[k] = dk;
300            lemax = lemax.max(dk);
301            mark = wclear(mark + lemax, lemax, &mut w, n);
302            (pk1..pk2).for_each(|pk| {
303                let i = ci[pk as usize] as usize;
304                if nv[i] < 0 {
305                    let hash = last[i] as usize;
306                    let mut i = hhead[hash];
307                    hhead[hash] = NONE;
308                    while i != NONE && next[i as usize] != NONE {
309                        let iu = i as usize;
310                        let ln = len[iu];
311                        let eln = elen[iu];
312                        (cp[iu] + 1..cp[iu] + ln).for_each(|p| w[ci[p as usize] as usize] = mark);
313                        let mut jlast = iu;
314                        let mut j = next[iu];
315                        while j != NONE {
316                            let ju = j as usize;
317                            let mut ok = len[ju] == ln && elen[ju] == eln;
318                            let mut p = cp[ju] + 1;
319                            while ok && p < cp[ju] + ln {
320                                if w[ci[p as usize] as usize] != mark {
321                                    ok = false;
322                                }
323                                p += 1;
324                            }
325                            if ok {
326                                cp[ju] = flip(i);
327                                nv[iu] += nv[ju];
328                                nv[ju] = 0;
329                                elen[ju] = NONE;
330                                j = next[ju];
331                                next[jlast] = j;
332                            } else {
333                                jlast = ju;
334                                j = next[ju];
335                            }
336                        }
337                        i = next[iu];
338                        mark += 1;
339                    }
340                }
341            });
342            let mut p = pk1;
343            (pk1..pk2).for_each(|pk| {
344                let i = ci[pk as usize] as usize;
345                let nvi = -nv[i];
346                if nvi > 0 {
347                    nv[i] = nvi;
348                    let d = (degree[i] + dk - nvi).min(n as isize - nel - nvi);
349                    if head[d as usize] != NONE {
350                        last[head[d as usize] as usize] = i as isize;
351                    }
352                    next[i] = head[d as usize];
353                    last[i] = NONE;
354                    head[d as usize] = i as isize;
355                    mindeg = mindeg.min(d as usize);
356                    degree[i] = d;
357                    ci[p as usize] = i as isize;
358                    p += 1;
359                }
360            });
361            nv[k] = nvk;
362            len[k] = p - pk1;
363            if len[k] == 0 {
364                cp[k] = NONE;
365                w[k] = 0;
366            }
367            if elenk != 0 {
368                cnz = p;
369            }
370        }
371        (0..=n).for_each(|i| cp[i] = flip(cp[i]));
372        head.iter_mut().for_each(|head_i| *head_i = NONE);
373        (0..n).rev().for_each(|j| {
374            if nv[j] <= 0 {
375                next[j] = head[cp[j] as usize];
376                head[cp[j] as usize] = j as isize;
377            }
378        });
379        (0..=n).rev().for_each(|e| {
380            if nv[e] > 0 && cp[e] != NONE {
381                next[e] = head[cp[e] as usize];
382                head[cp[e] as usize] = e as isize;
383            }
384        });
385        let mut post = vec![0_isize; n + 1];
386        let mut stack = vec![0_isize; n + 1];
387        let mut count = 0_usize;
388        (0..=n).for_each(|i| {
389            if cp[i] == NONE {
390                count = tdfs(i as isize, count, &mut head, &next, &mut post, &mut stack);
391            }
392        });
393        post.into_iter()
394            .filter(|&i| i != n as isize)
395            .map(|i| i as usize)
396            .collect()
397    }
398}
conspire/math/sparse/matrix/amd/mod.rs

conspire/math/sparse/matrix/amd/
mod.rs
