1 | import { factory } from '../../../utils/factory'
|
2 | import { csEreach } from './csEreach'
|
3 | import { createCsSymperm } from './csSymperm'
|
4 |
|
5 | const name = 'csChol'
|
6 | const dependencies = [
|
7 | 'divideScalar',
|
8 | 'sqrt',
|
9 | 'subtract',
|
10 | 'multiply',
|
11 | 'im',
|
12 | 're',
|
13 | 'conj',
|
14 | 'equal',
|
15 | 'smallerEq',
|
16 | 'SparseMatrix'
|
17 | ]
|
18 |
|
19 | export const createCsChol = factory(name, dependencies, (
|
20 | {
|
21 | divideScalar,
|
22 | sqrt,
|
23 | subtract,
|
24 | multiply,
|
25 | im,
|
26 | re,
|
27 | conj,
|
28 | equal,
|
29 | smallerEq,
|
30 | SparseMatrix
|
31 | }
|
32 | ) => {
|
33 | const csSymperm = createCsSymperm({ conj, SparseMatrix })
|
34 |
|
35 | |
36 |
|
37 |
|
38 |
|
39 |
|
40 |
|
41 |
|
42 |
|
43 |
|
44 |
|
45 |
|
46 | return function csChol (m, s) {
|
47 |
|
48 | if (!m) { return null }
|
49 |
|
50 | const size = m._size
|
51 |
|
52 | const n = size[1]
|
53 |
|
54 | const parent = s.parent
|
55 | const cp = s.cp
|
56 | const pinv = s.pinv
|
57 |
|
58 | const lvalues = []
|
59 | const lindex = []
|
60 | const lptr = []
|
61 |
|
62 | const L = new SparseMatrix({
|
63 | values: lvalues,
|
64 | index: lindex,
|
65 | ptr: lptr,
|
66 | size: [n, n]
|
67 | })
|
68 |
|
69 | const c = []
|
70 | const x = []
|
71 |
|
72 | const cm = pinv ? csSymperm(m, pinv, 1) : m
|
73 |
|
74 | const cvalues = cm._values
|
75 | const cindex = cm._index
|
76 | const cptr = cm._ptr
|
77 |
|
78 | let k, p
|
79 |
|
80 | for (k = 0; k < n; k++) { lptr[k] = c[k] = cp[k] }
|
81 |
|
82 | for (k = 0; k < n; k++) {
|
83 |
|
84 | let top = csEreach(cm, k, parent, c)
|
85 |
|
86 | x[k] = 0
|
87 |
|
88 | for (p = cptr[k]; p < cptr[k + 1]; p++) {
|
89 | if (cindex[p] <= k) { x[cindex[p]] = cvalues[p] }
|
90 | }
|
91 |
|
92 | let d = x[k]
|
93 |
|
94 | x[k] = 0
|
95 |
|
96 | for (; top < n; top++) {
|
97 |
|
98 | const i = s[top]
|
99 |
|
100 | const lki = divideScalar(x[i], lvalues[lptr[i]])
|
101 |
|
102 | x[i] = 0
|
103 | for (p = lptr[i] + 1; p < c[i]; p++) {
|
104 |
|
105 | const r = lindex[p]
|
106 |
|
107 | x[r] = subtract(x[r], multiply(lvalues[p], lki))
|
108 | }
|
109 |
|
110 | d = subtract(d, multiply(lki, conj(lki)))
|
111 | p = c[i]++
|
112 |
|
113 | lindex[p] = k
|
114 | lvalues[p] = conj(lki)
|
115 | }
|
116 |
|
117 | if (smallerEq(re(d), 0) || !equal(im(d), 0)) {
|
118 |
|
119 | return null
|
120 | }
|
121 | p = c[k]++
|
122 |
|
123 | lindex[p] = k
|
124 | lvalues[p] = sqrt(d)
|
125 | }
|
126 |
|
127 | lptr[n] = cp[n]
|
128 |
|
129 | let P
|
130 |
|
131 | if (pinv) {
|
132 |
|
133 | const pvalues = []
|
134 | const pindex = []
|
135 | const pptr = []
|
136 |
|
137 | for (p = 0; p < n; p++) {
|
138 |
|
139 | pptr[p] = p
|
140 |
|
141 | pindex.push(pinv[p])
|
142 |
|
143 | pvalues.push(1)
|
144 | }
|
145 |
|
146 | pptr[n] = n
|
147 |
|
148 | P = new SparseMatrix({
|
149 | values: pvalues,
|
150 | index: pindex,
|
151 | ptr: pptr,
|
152 | size: [n, n]
|
153 | })
|
154 | }
|
155 |
|
156 | return {
|
157 | L: L,
|
158 | P: P
|
159 | }
|
160 | }
|
161 | })
|