csAmd.js (17630B)
1 import { factory } from '../../../utils/factory.js'; 2 import { csFkeep } from './csFkeep.js'; 3 import { csFlip } from './csFlip.js'; 4 import { csTdfs } from './csTdfs.js'; 5 var name = 'csAmd'; 6 var dependencies = ['add', 'multiply', 'transpose']; 7 export var createCsAmd = /* #__PURE__ */factory(name, dependencies, _ref => { 8 var { 9 add, 10 multiply, 11 transpose 12 } = _ref; 13 14 /** 15 * Approximate minimum degree ordering. The minimum degree algorithm is a widely used 16 * heuristic for finding a permutation P so that P*A*P' has fewer nonzeros in its factorization 17 * than A. It is a gready method that selects the sparsest pivot row and column during the course 18 * of a right looking sparse Cholesky factorization. 19 * 20 * Reference: http://faculty.cse.tamu.edu/davis/publications.html 21 * 22 * @param {Number} order 0: Natural, 1: Cholesky, 2: LU, 3: QR 23 * @param {Matrix} m Sparse Matrix 24 */ 25 return function csAmd(order, a) { 26 // check input parameters 27 if (!a || order <= 0 || order > 3) { 28 return null; 29 } // a matrix arrays 30 31 32 var asize = a._size; // rows and columns 33 34 var m = asize[0]; 35 var n = asize[1]; // initialize vars 36 37 var lemax = 0; // dense threshold 38 39 var dense = Math.max(16, 10 * Math.sqrt(n)); 40 dense = Math.min(n - 2, dense); // create target matrix C 41 42 var cm = _createTargetMatrix(order, a, m, n, dense); // drop diagonal entries 43 44 45 csFkeep(cm, _diag, null); // C matrix arrays 46 47 var cindex = cm._index; 48 var cptr = cm._ptr; // number of nonzero elements in C 49 50 var cnz = cptr[n]; // allocate result (n+1) 51 52 var P = []; // create workspace (8 * (n + 1)) 53 54 var W = []; 55 var len = 0; // first n + 1 entries 56 57 var nv = n + 1; // next n + 1 entries 58 59 var next = 2 * (n + 1); // next n + 1 entries 60 61 var head = 3 * (n + 1); // next n + 1 entries 62 63 var elen = 4 * (n + 1); // next n + 1 entries 64 65 var degree = 5 * (n + 1); // next n + 1 entries 66 67 var w = 6 * (n + 1); // next n + 1 entries 68 69 var hhead = 7 * (n + 1); // last n + 1 entries 70 // use P as workspace for last 71 72 var last = P; // initialize quotient graph 73 74 var mark = _initializeQuotientGraph(n, cptr, W, len, head, last, next, hhead, nv, w, elen, degree); // initialize degree lists 75 76 77 var nel = _initializeDegreeLists(n, cptr, W, degree, elen, w, dense, nv, head, last, next); // minimum degree node 78 79 80 var mindeg = 0; // vars 81 82 var i, j, k, k1, k2, e, pj, ln, nvi, pk, eln, p1, p2, pn, h, d; // while (selecting pivots) do 83 84 while (nel < n) { 85 // select node of minimum approximate degree. amd() is now ready to start eliminating the graph. It first 86 // finds a node k of minimum degree and removes it from its degree list. The variable nel keeps track of thow 87 // many nodes have been eliminated. 88 for (k = -1; mindeg < n && (k = W[head + mindeg]) === -1; mindeg++) { 89 ; 90 } 91 92 if (W[next + k] !== -1) { 93 last[W[next + k]] = -1; 94 } // remove k from degree list 95 96 97 W[head + mindeg] = W[next + k]; // elenk = |Ek| 98 99 var elenk = W[elen + k]; // # of nodes k represents 100 101 var nvk = W[nv + k]; // W[nv + k] nodes of A eliminated 102 103 nel += nvk; // Construct a new element. The new element Lk is constructed in place if |Ek| = 0. nv[i] is 104 // negated for all nodes i in Lk to flag them as members of this set. Each node i is removed from the 105 // degree lists. All elements e in Ek are absorved into element k. 106 107 var dk = 0; // flag k as in Lk 108 109 W[nv + k] = -nvk; 110 var p = cptr[k]; // do in place if W[elen + k] === 0 111 112 var pk1 = elenk === 0 ? p : cnz; 113 var pk2 = pk1; 114 115 for (k1 = 1; k1 <= elenk + 1; k1++) { 116 if (k1 > elenk) { 117 // search the nodes in k 118 e = k; // list of nodes starts at cindex[pj] 119 120 pj = p; // length of list of nodes in k 121 122 ln = W[len + k] - elenk; 123 } else { 124 // search the nodes in e 125 e = cindex[p++]; 126 pj = cptr[e]; // length of list of nodes in e 127 128 ln = W[len + e]; 129 } 130 131 for (k2 = 1; k2 <= ln; k2++) { 132 i = cindex[pj++]; // check node i dead, or seen 133 134 if ((nvi = W[nv + i]) <= 0) { 135 continue; 136 } // W[degree + Lk] += size of node i 137 138 139 dk += nvi; // negate W[nv + i] to denote i in Lk 140 141 W[nv + i] = -nvi; // place i in Lk 142 143 cindex[pk2++] = i; 144 145 if (W[next + i] !== -1) { 146 last[W[next + i]] = last[i]; 147 } // check we need to remove i from degree list 148 149 150 if (last[i] !== -1) { 151 W[next + last[i]] = W[next + i]; 152 } else { 153 W[head + W[degree + i]] = W[next + i]; 154 } 155 } 156 157 if (e !== k) { 158 // absorb e into k 159 cptr[e] = csFlip(k); // e is now a dead element 160 161 W[w + e] = 0; 162 } 163 } // cindex[cnz...nzmax] is free 164 165 166 if (elenk !== 0) { 167 cnz = pk2; 168 } // external degree of k - |Lk\i| 169 170 171 W[degree + k] = dk; // element k is in cindex[pk1..pk2-1] 172 173 cptr[k] = pk1; 174 W[len + k] = pk2 - pk1; // k is now an element 175 176 W[elen + k] = -2; // Find set differences. The scan1 function now computes the set differences |Le \ Lk| for all elements e. At the start of the 177 // scan, no entry in the w array is greater than or equal to mark. 178 // clear w if necessary 179 180 mark = _wclear(mark, lemax, W, w, n); // scan 1: find |Le\Lk| 181 182 for (pk = pk1; pk < pk2; pk++) { 183 i = cindex[pk]; // check if W[elen + i] empty, skip it 184 185 if ((eln = W[elen + i]) <= 0) { 186 continue; 187 } // W[nv + i] was negated 188 189 190 nvi = -W[nv + i]; 191 var wnvi = mark - nvi; // scan Ei 192 193 for (p = cptr[i], p1 = cptr[i] + eln - 1; p <= p1; p++) { 194 e = cindex[p]; 195 196 if (W[w + e] >= mark) { 197 // decrement |Le\Lk| 198 W[w + e] -= nvi; 199 } else if (W[w + e] !== 0) { 200 // ensure e is a live element, 1st time e seen in scan 1 201 W[w + e] = W[degree + e] + wnvi; 202 } 203 } 204 } // degree update 205 // The second pass computes the approximate degree di, prunes the sets Ei and Ai, and computes a hash 206 // function h(i) for all nodes in Lk. 207 // scan2: degree update 208 209 210 for (pk = pk1; pk < pk2; pk++) { 211 // consider node i in Lk 212 i = cindex[pk]; 213 p1 = cptr[i]; 214 p2 = p1 + W[elen + i] - 1; 215 pn = p1; // scan Ei 216 217 for (h = 0, d = 0, p = p1; p <= p2; p++) { 218 e = cindex[p]; // check e is an unabsorbed element 219 220 if (W[w + e] !== 0) { 221 // dext = |Le\Lk| 222 var dext = W[w + e] - mark; 223 224 if (dext > 0) { 225 // sum up the set differences 226 d += dext; // keep e in Ei 227 228 cindex[pn++] = e; // compute the hash of node i 229 230 h += e; 231 } else { 232 // aggressive absorb. e->k 233 cptr[e] = csFlip(k); // e is a dead element 234 235 W[w + e] = 0; 236 } 237 } 238 } // W[elen + i] = |Ei| 239 240 241 W[elen + i] = pn - p1 + 1; 242 var p3 = pn; 243 var p4 = p1 + W[len + i]; // prune edges in Ai 244 245 for (p = p2 + 1; p < p4; p++) { 246 j = cindex[p]; // check node j dead or in Lk 247 248 var nvj = W[nv + j]; 249 250 if (nvj <= 0) { 251 continue; 252 } // degree(i) += |j| 253 254 255 d += nvj; // place j in node list of i 256 257 cindex[pn++] = j; // compute hash for node i 258 259 h += j; 260 } // check for mass elimination 261 262 263 if (d === 0) { 264 // absorb i into k 265 cptr[i] = csFlip(k); 266 nvi = -W[nv + i]; // |Lk| -= |i| 267 268 dk -= nvi; // |k| += W[nv + i] 269 270 nvk += nvi; 271 nel += nvi; 272 W[nv + i] = 0; // node i is dead 273 274 W[elen + i] = -1; 275 } else { 276 // update degree(i) 277 W[degree + i] = Math.min(W[degree + i], d); // move first node to end 278 279 cindex[pn] = cindex[p3]; // move 1st el. to end of Ei 280 281 cindex[p3] = cindex[p1]; // add k as 1st element in of Ei 282 283 cindex[p1] = k; // new len of adj. list of node i 284 285 W[len + i] = pn - p1 + 1; // finalize hash of i 286 287 h = (h < 0 ? -h : h) % n; // place i in hash bucket 288 289 W[next + i] = W[hhead + h]; 290 W[hhead + h] = i; // save hash of i in last[i] 291 292 last[i] = h; 293 } 294 } // finalize |Lk| 295 296 297 W[degree + k] = dk; 298 lemax = Math.max(lemax, dk); // clear w 299 300 mark = _wclear(mark + lemax, lemax, W, w, n); // Supernode detection. Supernode detection relies on the hash function h(i) computed for each node i. 301 // If two nodes have identical adjacency lists, their hash functions wil be identical. 302 303 for (pk = pk1; pk < pk2; pk++) { 304 i = cindex[pk]; // check i is dead, skip it 305 306 if (W[nv + i] >= 0) { 307 continue; 308 } // scan hash bucket of node i 309 310 311 h = last[i]; 312 i = W[hhead + h]; // hash bucket will be empty 313 314 W[hhead + h] = -1; 315 316 for (; i !== -1 && W[next + i] !== -1; i = W[next + i], mark++) { 317 ln = W[len + i]; 318 eln = W[elen + i]; 319 320 for (p = cptr[i] + 1; p <= cptr[i] + ln - 1; p++) { 321 W[w + cindex[p]] = mark; 322 } 323 324 var jlast = i; // compare i with all j 325 326 for (j = W[next + i]; j !== -1;) { 327 var ok = W[len + j] === ln && W[elen + j] === eln; 328 329 for (p = cptr[j] + 1; ok && p <= cptr[j] + ln - 1; p++) { 330 // compare i and j 331 if (W[w + cindex[p]] !== mark) { 332 ok = 0; 333 } 334 } // check i and j are identical 335 336 337 if (ok) { 338 // absorb j into i 339 cptr[j] = csFlip(i); 340 W[nv + i] += W[nv + j]; 341 W[nv + j] = 0; // node j is dead 342 343 W[elen + j] = -1; // delete j from hash bucket 344 345 j = W[next + j]; 346 W[next + jlast] = j; 347 } else { 348 // j and i are different 349 jlast = j; 350 j = W[next + j]; 351 } 352 } 353 } 354 } // Finalize new element. The elimination of node k is nearly complete. All nodes i in Lk are scanned one last time. 355 // Node i is removed from Lk if it is dead. The flagged status of nv[i] is cleared. 356 357 358 for (p = pk1, pk = pk1; pk < pk2; pk++) { 359 i = cindex[pk]; // check i is dead, skip it 360 361 if ((nvi = -W[nv + i]) <= 0) { 362 continue; 363 } // restore W[nv + i] 364 365 366 W[nv + i] = nvi; // compute external degree(i) 367 368 d = W[degree + i] + dk - nvi; 369 d = Math.min(d, n - nel - nvi); 370 371 if (W[head + d] !== -1) { 372 last[W[head + d]] = i; 373 } // put i back in degree list 374 375 376 W[next + i] = W[head + d]; 377 last[i] = -1; 378 W[head + d] = i; // find new minimum degree 379 380 mindeg = Math.min(mindeg, d); 381 W[degree + i] = d; // place i in Lk 382 383 cindex[p++] = i; 384 } // # nodes absorbed into k 385 386 387 W[nv + k] = nvk; // length of adj list of element k 388 389 if ((W[len + k] = p - pk1) === 0) { 390 // k is a root of the tree 391 cptr[k] = -1; // k is now a dead element 392 393 W[w + k] = 0; 394 } 395 396 if (elenk !== 0) { 397 // free unused space in Lk 398 cnz = p; 399 } 400 } // Postordering. The elimination is complete, but no permutation has been computed. All that is left 401 // of the graph is the assembly tree (ptr) and a set of dead nodes and elements (i is a dead node if 402 // nv[i] is zero and a dead element if nv[i] > 0). It is from this information only that the final permutation 403 // is computed. The tree is restored by unflipping all of ptr. 404 // fix assembly tree 405 406 407 for (i = 0; i < n; i++) { 408 cptr[i] = csFlip(cptr[i]); 409 } 410 411 for (j = 0; j <= n; j++) { 412 W[head + j] = -1; 413 } // place unordered nodes in lists 414 415 416 for (j = n; j >= 0; j--) { 417 // skip if j is an element 418 if (W[nv + j] > 0) { 419 continue; 420 } // place j in list of its parent 421 422 423 W[next + j] = W[head + cptr[j]]; 424 W[head + cptr[j]] = j; 425 } // place elements in lists 426 427 428 for (e = n; e >= 0; e--) { 429 // skip unless e is an element 430 if (W[nv + e] <= 0) { 431 continue; 432 } 433 434 if (cptr[e] !== -1) { 435 // place e in list of its parent 436 W[next + e] = W[head + cptr[e]]; 437 W[head + cptr[e]] = e; 438 } 439 } // postorder the assembly tree 440 441 442 for (k = 0, i = 0; i <= n; i++) { 443 if (cptr[i] === -1) { 444 k = csTdfs(i, k, W, head, next, P, w); 445 } 446 } // remove last item in array 447 448 449 P.splice(P.length - 1, 1); // return P 450 451 return P; 452 }; 453 /** 454 * Creates the matrix that will be used by the approximate minimum degree ordering algorithm. The function accepts the matrix M as input and returns a permutation 455 * vector P. The amd algorithm operates on a symmetrix matrix, so one of three symmetric matrices is formed. 456 * 457 * Order: 0 458 * A natural ordering P=null matrix is returned. 459 * 460 * Order: 1 461 * Matrix must be square. This is appropriate for a Cholesky or LU factorization. 462 * P = M + M' 463 * 464 * Order: 2 465 * Dense columns from M' are dropped, M recreated from M'. This is appropriatefor LU factorization of unsymmetric matrices. 466 * P = M' * M 467 * 468 * Order: 3 469 * This is best used for QR factorization or LU factorization is matrix M has no dense rows. A dense row is a row with more than 10*sqr(columns) entries. 470 * P = M' * M 471 */ 472 473 function _createTargetMatrix(order, a, m, n, dense) { 474 // compute A' 475 var at = transpose(a); // check order = 1, matrix must be square 476 477 if (order === 1 && n === m) { 478 // C = A + A' 479 return add(a, at); 480 } // check order = 2, drop dense columns from M' 481 482 483 if (order === 2) { 484 // transpose arrays 485 var tindex = at._index; 486 var tptr = at._ptr; // new column index 487 488 var p2 = 0; // loop A' columns (rows) 489 490 for (var j = 0; j < m; j++) { 491 // column j of AT starts here 492 var p = tptr[j]; // new column j starts here 493 494 tptr[j] = p2; // skip dense col j 495 496 if (tptr[j + 1] - p > dense) { 497 continue; 498 } // map rows in column j of A 499 500 501 for (var p1 = tptr[j + 1]; p < p1; p++) { 502 tindex[p2++] = tindex[p]; 503 } 504 } // finalize AT 505 506 507 tptr[m] = p2; // recreate A from new transpose matrix 508 509 a = transpose(at); // use A' * A 510 511 return multiply(at, a); 512 } // use A' * A, square or rectangular matrix 513 514 515 return multiply(at, a); 516 } 517 /** 518 * Initialize quotient graph. There are four kind of nodes and elements that must be represented: 519 * 520 * - A live node is a node i (or a supernode) that has not been selected as a pivot nad has not been merged into another supernode. 521 * - A dead node i is one that has been removed from the graph, having been absorved into r = flip(ptr[i]). 522 * - A live element e is one that is in the graph, having been formed when node e was selected as the pivot. 523 * - A dead element e is one that has benn absorved into a subsequent element s = flip(ptr[e]). 524 */ 525 526 527 function _initializeQuotientGraph(n, cptr, W, len, head, last, next, hhead, nv, w, elen, degree) { 528 // Initialize quotient graph 529 for (var k = 0; k < n; k++) { 530 W[len + k] = cptr[k + 1] - cptr[k]; 531 } 532 533 W[len + n] = 0; // initialize workspace 534 535 for (var i = 0; i <= n; i++) { 536 // degree list i is empty 537 W[head + i] = -1; 538 last[i] = -1; 539 W[next + i] = -1; // hash list i is empty 540 541 W[hhead + i] = -1; // node i is just one node 542 543 W[nv + i] = 1; // node i is alive 544 545 W[w + i] = 1; // Ek of node i is empty 546 547 W[elen + i] = 0; // degree of node i 548 549 W[degree + i] = W[len + i]; 550 } // clear w 551 552 553 var mark = _wclear(0, 0, W, w, n); // n is a dead element 554 555 556 W[elen + n] = -2; // n is a root of assembly tree 557 558 cptr[n] = -1; // n is a dead element 559 560 W[w + n] = 0; // return mark 561 562 return mark; 563 } 564 /** 565 * Initialize degree lists. Each node is placed in its degree lists. Nodes of zero degree are eliminated immediately. Nodes with 566 * degree >= dense are alsol eliminated and merged into a placeholder node n, a dead element. Thes nodes will appera last in the 567 * output permutation p. 568 */ 569 570 571 function _initializeDegreeLists(n, cptr, W, degree, elen, w, dense, nv, head, last, next) { 572 // result 573 var nel = 0; // loop columns 574 575 for (var i = 0; i < n; i++) { 576 // degree @ i 577 var d = W[degree + i]; // check node i is empty 578 579 if (d === 0) { 580 // element i is dead 581 W[elen + i] = -2; 582 nel++; // i is a root of assembly tree 583 584 cptr[i] = -1; 585 W[w + i] = 0; 586 } else if (d > dense) { 587 // absorb i into element n 588 W[nv + i] = 0; // node i is dead 589 590 W[elen + i] = -1; 591 nel++; 592 cptr[i] = csFlip(n); 593 W[nv + n]++; 594 } else { 595 var h = W[head + d]; 596 597 if (h !== -1) { 598 last[h] = i; 599 } // put node i in degree list d 600 601 602 W[next + i] = W[head + d]; 603 W[head + d] = i; 604 } 605 } 606 607 return nel; 608 } 609 610 function _wclear(mark, lemax, W, w, n) { 611 if (mark < 2 || mark + lemax < 0) { 612 for (var k = 0; k < n; k++) { 613 if (W[w + k] !== 0) { 614 W[w + k] = 1; 615 } 616 } 617 618 mark = 2; 619 } // at this point, W [0..n-1] < mark holds 620 621 622 return mark; 623 } 624 625 function _diag(i, j) { 626 return i !== j; 627 } 628 });