| 1 | /**
|
| 2 | * Computes the elimination tree of Matrix A (using triu(A)) or the
|
| 3 | * elimination tree of A'A without forming A'A.
|
| 4 | *
|
| 5 | * @param {Matrix} a The A Matrix
|
| 6 | * @param {boolean} ata A value of true the function computes the etree of A'A
|
| 7 | *
|
| 8 | * Reference: http://faculty.cse.tamu.edu/davis/publications.html
|
| 9 | */
|
| 10 | export function csEtree (a, ata) {
|
| 11 | // check inputs
|
| 12 | if (!a) { return null }
|
| 13 | // a arrays
|
| 14 | const aindex = a._index
|
| 15 | const aptr = a._ptr
|
| 16 | const asize = a._size
|
| 17 | // rows & columns
|
| 18 | const m = asize[0]
|
| 19 | const n = asize[1]
|
| 20 |
|
| 21 | // allocate result
|
| 22 | const parent = [] // (n)
|
| 23 |
|
| 24 | // allocate workspace
|
| 25 | const w = [] // (n + (ata ? m : 0))
|
| 26 | const ancestor = 0 // first n entries in w
|
| 27 | const prev = n // last m entries (ata = true)
|
| 28 |
|
| 29 | let i, inext
|
| 30 |
|
| 31 | // check we are calculating A'A
|
| 32 | if (ata) {
|
| 33 | // initialize workspace
|
| 34 | for (i = 0; i < m; i++) { w[prev + i] = -1 }
|
| 35 | }
|
| 36 | // loop columns
|
| 37 | for (let k = 0; k < n; k++) {
|
| 38 | // node k has no parent yet
|
| 39 | parent[k] = -1
|
| 40 | // nor does k have an ancestor
|
| 41 | w[ancestor + k] = -1
|
| 42 | // values in column k
|
| 43 | for (let p0 = aptr[k], p1 = aptr[k + 1], p = p0; p < p1; p++) {
|
| 44 | // row
|
| 45 | const r = aindex[p]
|
| 46 | // node
|
| 47 | i = ata ? (w[prev + r]) : r
|
| 48 | // traverse from i to k
|
| 49 | for (; i !== -1 && i < k; i = inext) {
|
| 50 | // inext = ancestor of i
|
| 51 | inext = w[ancestor + i]
|
| 52 | // path compression
|
| 53 | w[ancestor + i] = k
|
| 54 | // check no anc., parent is k
|
| 55 | if (inext === -1) { parent[i] = k }
|
| 56 | }
|
| 57 | if (ata) { w[prev + r] = k }
|
| 58 | }
|
| 59 | }
|
| 60 | return parent
|
| 61 | }
|