csSpsolve.js 3.1 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788
  1. "use strict";
  2. Object.defineProperty(exports, "__esModule", {
  3. value: true
  4. });
  5. exports.createCsSpsolve = void 0;
  6. var _csReach = require("./csReach.js");
  7. var _factory = require("../../../utils/factory.js");
  8. var name = 'csSpsolve';
  9. var dependencies = ['divideScalar', 'multiply', 'subtract'];
  10. var createCsSpsolve = /* #__PURE__ */(0, _factory.factory)(name, dependencies, function (_ref) {
  11. var divideScalar = _ref.divideScalar,
  12. multiply = _ref.multiply,
  13. subtract = _ref.subtract;
  14. /**
  15. * The function csSpsolve() computes the solution to G * x = bk, where bk is the
  16. * kth column of B. When lo is true, the function assumes G = L is lower triangular with the
  17. * diagonal entry as the first entry in each column. When lo is true, the function assumes G = U
  18. * is upper triangular with the diagonal entry as the last entry in each column.
  19. *
  20. * @param {Matrix} g The G matrix
  21. * @param {Matrix} b The B matrix
  22. * @param {Number} k The kth column in B
  23. * @param {Array} xi The nonzero pattern xi[top] .. xi[n - 1], an array of size = 2 * n
  24. * The first n entries is the nonzero pattern, the last n entries is the stack
  25. * @param {Array} x The soluton to the linear system G * x = b
  26. * @param {Array} pinv The inverse row permutation vector, must be null for L * x = b
  27. * @param {boolean} lo The lower (true) upper triangular (false) flag
  28. *
  29. * @return {Number} The index for the nonzero pattern
  30. *
  31. * Reference: http://faculty.cse.tamu.edu/davis/publications.html
  32. */
  33. return function csSpsolve(g, b, k, xi, x, pinv, lo) {
  34. // g arrays
  35. var gvalues = g._values;
  36. var gindex = g._index;
  37. var gptr = g._ptr;
  38. var gsize = g._size;
  39. // columns
  40. var n = gsize[1];
  41. // b arrays
  42. var bvalues = b._values;
  43. var bindex = b._index;
  44. var bptr = b._ptr;
  45. // vars
  46. var p, p0, p1, q;
  47. // xi[top..n-1] = csReach(B(:,k))
  48. var top = (0, _csReach.csReach)(g, b, k, xi, pinv);
  49. // clear x
  50. for (p = top; p < n; p++) {
  51. x[xi[p]] = 0;
  52. }
  53. // scatter b
  54. for (p0 = bptr[k], p1 = bptr[k + 1], p = p0; p < p1; p++) {
  55. x[bindex[p]] = bvalues[p];
  56. }
  57. // loop columns
  58. for (var px = top; px < n; px++) {
  59. // x array index for px
  60. var j = xi[px];
  61. // apply permutation vector (U x = b), j maps to column J of G
  62. var J = pinv ? pinv[j] : j;
  63. // check column J is empty
  64. if (J < 0) {
  65. continue;
  66. }
  67. // column value indeces in G, p0 <= p < p1
  68. p0 = gptr[J];
  69. p1 = gptr[J + 1];
  70. // x(j) /= G(j,j)
  71. x[j] = divideScalar(x[j], gvalues[lo ? p0 : p1 - 1]);
  72. // first entry L(j,j)
  73. p = lo ? p0 + 1 : p0;
  74. q = lo ? p1 : p1 - 1;
  75. // loop
  76. for (; p < q; p++) {
  77. // row
  78. var i = gindex[p];
  79. // x(i) -= G(i,j) * x(j)
  80. x[i] = subtract(x[i], multiply(gvalues[p], x[j]));
  81. }
  82. }
  83. // return top of stack
  84. return top;
  85. };
  86. });
  87. exports.createCsSpsolve = createCsSpsolve;