1 | import { csCumsum } from './csCumsum'
|
2 | import { factory } from '../../../utils/factory'
|
3 |
|
4 | const name = 'csSymperm'
|
5 | const dependencies = ['conj', 'SparseMatrix']
|
6 |
|
7 | export const createCsSymperm = /* #__PURE__ */ factory(name, dependencies, ({ conj, SparseMatrix }) => {
|
8 | /**
|
9 | * Computes the symmetric permutation of matrix A accessing only
|
10 | * the upper triangular part of A.
|
11 | *
|
12 | * C = P * A * P'
|
13 | *
|
14 | * @param {Matrix} a The A matrix
|
15 | * @param {Array} pinv The inverse of permutation vector
|
16 | * @param {boolean} values Process matrix values (true)
|
17 | *
|
18 | * @return {Matrix} The C matrix, C = P * A * P'
|
19 | *
|
20 | * Reference: http://faculty.cse.tamu.edu/davis/publications.html
|
21 | */
|
22 | return function csSymperm (a, pinv, values) {
|
23 | // A matrix arrays
|
24 | const avalues = a._values
|
25 | const aindex = a._index
|
26 | const aptr = a._ptr
|
27 | const asize = a._size
|
28 | // columns
|
29 | const n = asize[1]
|
30 | // C matrix arrays
|
31 | const cvalues = values && avalues ? [] : null
|
32 | const cindex = [] // (nz)
|
33 | const cptr = [] // (n + 1)
|
34 | // variables
|
35 | let i, i2, j, j2, p, p0, p1
|
36 | // create workspace vector
|
37 | const w = [] // (n)
|
38 | // count entries in each column of C
|
39 | for (j = 0; j < n; j++) {
|
40 | // column j of A is column j2 of C
|
41 | j2 = pinv ? pinv[j] : j
|
42 | // loop values in column j
|
43 | for (p0 = aptr[j], p1 = aptr[j + 1], p = p0; p < p1; p++) {
|
44 | // row
|
45 | i = aindex[p]
|
46 | // skip lower triangular part of A
|
47 | if (i > j) { continue }
|
48 | // row i of A is row i2 of C
|
49 | i2 = pinv ? pinv[i] : i
|
50 | // column count of C
|
51 | w[Math.max(i2, j2)]++
|
52 | }
|
53 | }
|
54 | // compute column pointers of C
|
55 | csCumsum(cptr, w, n)
|
56 | // loop columns
|
57 | for (j = 0; j < n; j++) {
|
58 | // column j of A is column j2 of C
|
59 | j2 = pinv ? pinv[j] : j
|
60 | // loop values in column j
|
61 | for (p0 = aptr[j], p1 = aptr[j + 1], p = p0; p < p1; p++) {
|
62 | // row
|
63 | i = aindex[p]
|
64 | // skip lower triangular part of A
|
65 | if (i > j) { continue }
|
66 | // row i of A is row i2 of C
|
67 | i2 = pinv ? pinv[i] : i
|
68 | // C index for column j2
|
69 | const q = w[Math.max(i2, j2)]++
|
70 | // update C index for entry q
|
71 | cindex[q] = Math.min(i2, j2)
|
72 | // check we need to process values
|
73 | if (cvalues) { cvalues[q] = (i2 <= j2) ? avalues[p] : conj(avalues[p]) }
|
74 | }
|
75 | }
|
76 | // return C matrix
|
77 | return new SparseMatrix({
|
78 | values: cvalues,
|
79 | index: cindex,
|
80 | ptr: cptr,
|
81 | size: [n, n]
|
82 | })
|
83 | }
|
84 | })
|