csSqr.js 5.0 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183
  1. "use strict";
  2. Object.defineProperty(exports, "__esModule", {
  3. value: true
  4. });
  5. exports.createCsSqr = void 0;
  6. var _csPermute = require("./csPermute.js");
  7. var _csPost = require("./csPost.js");
  8. var _csEtree = require("./csEtree.js");
  9. var _csAmd = require("./csAmd.js");
  10. var _csCounts = require("./csCounts.js");
  11. var _factory = require("../../../utils/factory.js");
  12. var name = 'csSqr';
  13. var dependencies = ['add', 'multiply', 'transpose'];
  14. var createCsSqr = /* #__PURE__ */(0, _factory.factory)(name, dependencies, function (_ref) {
  15. var add = _ref.add,
  16. multiply = _ref.multiply,
  17. transpose = _ref.transpose;
  18. var csAmd = (0, _csAmd.createCsAmd)({
  19. add: add,
  20. multiply: multiply,
  21. transpose: transpose
  22. });
  23. var csCounts = (0, _csCounts.createCsCounts)({
  24. transpose: transpose
  25. });
  26. /**
  27. * Symbolic ordering and analysis for QR and LU decompositions.
  28. *
  29. * @param {Number} order The ordering strategy (see csAmd for more details)
  30. * @param {Matrix} a The A matrix
  31. * @param {boolean} qr Symbolic ordering and analysis for QR decomposition (true) or
  32. * symbolic ordering and analysis for LU decomposition (false)
  33. *
  34. * @return {Object} The Symbolic ordering and analysis for matrix A
  35. *
  36. * Reference: http://faculty.cse.tamu.edu/davis/publications.html
  37. */
  38. return function csSqr(order, a, qr) {
  39. // a arrays
  40. var aptr = a._ptr;
  41. var asize = a._size;
  42. // columns
  43. var n = asize[1];
  44. // vars
  45. var k;
  46. // symbolic analysis result
  47. var s = {};
  48. // fill-reducing ordering
  49. s.q = csAmd(order, a);
  50. // validate results
  51. if (order && !s.q) {
  52. return null;
  53. }
  54. // QR symbolic analysis
  55. if (qr) {
  56. // apply permutations if needed
  57. var c = order ? (0, _csPermute.csPermute)(a, null, s.q, 0) : a;
  58. // etree of C'*C, where C=A(:,q)
  59. s.parent = (0, _csEtree.csEtree)(c, 1);
  60. // post order elimination tree
  61. var post = (0, _csPost.csPost)(s.parent, n);
  62. // col counts chol(C'*C)
  63. s.cp = csCounts(c, s.parent, post, 1);
  64. // check we have everything needed to calculate number of nonzero elements
  65. if (c && s.parent && s.cp && _vcount(c, s)) {
  66. // calculate number of nonzero elements
  67. for (s.unz = 0, k = 0; k < n; k++) {
  68. s.unz += s.cp[k];
  69. }
  70. }
  71. } else {
  72. // for LU factorization only, guess nnz(L) and nnz(U)
  73. s.unz = 4 * aptr[n] + n;
  74. s.lnz = s.unz;
  75. }
  76. // return result S
  77. return s;
  78. };
  79. /**
  80. * Compute nnz(V) = s.lnz, s.pinv, s.leftmost, s.m2 from A and s.parent
  81. */
  82. function _vcount(a, s) {
  83. // a arrays
  84. var aptr = a._ptr;
  85. var aindex = a._index;
  86. var asize = a._size;
  87. // rows & columns
  88. var m = asize[0];
  89. var n = asize[1];
  90. // initialize s arrays
  91. s.pinv = []; // (m + n)
  92. s.leftmost = []; // (m)
  93. // vars
  94. var parent = s.parent;
  95. var pinv = s.pinv;
  96. var leftmost = s.leftmost;
  97. // workspace, next: first m entries, head: next n entries, tail: next n entries, nque: next n entries
  98. var w = []; // (m + 3 * n)
  99. var next = 0;
  100. var head = m;
  101. var tail = m + n;
  102. var nque = m + 2 * n;
  103. // vars
  104. var i, k, p, p0, p1;
  105. // initialize w
  106. for (k = 0; k < n; k++) {
  107. // queue k is empty
  108. w[head + k] = -1;
  109. w[tail + k] = -1;
  110. w[nque + k] = 0;
  111. }
  112. // initialize row arrays
  113. for (i = 0; i < m; i++) {
  114. leftmost[i] = -1;
  115. }
  116. // loop columns backwards
  117. for (k = n - 1; k >= 0; k--) {
  118. // values & index for column k
  119. for (p0 = aptr[k], p1 = aptr[k + 1], p = p0; p < p1; p++) {
  120. // leftmost[i] = min(find(A(i,:)))
  121. leftmost[aindex[p]] = k;
  122. }
  123. }
  124. // scan rows in reverse order
  125. for (i = m - 1; i >= 0; i--) {
  126. // row i is not yet ordered
  127. pinv[i] = -1;
  128. k = leftmost[i];
  129. // check row i is empty
  130. if (k === -1) {
  131. continue;
  132. }
  133. // first row in queue k
  134. if (w[nque + k]++ === 0) {
  135. w[tail + k] = i;
  136. }
  137. // put i at head of queue k
  138. w[next + i] = w[head + k];
  139. w[head + k] = i;
  140. }
  141. s.lnz = 0;
  142. s.m2 = m;
  143. // find row permutation and nnz(V)
  144. for (k = 0; k < n; k++) {
  145. // remove row i from queue k
  146. i = w[head + k];
  147. // count V(k,k) as nonzero
  148. s.lnz++;
  149. // add a fictitious row
  150. if (i < 0) {
  151. i = s.m2++;
  152. }
  153. // associate row i with V(:,k)
  154. pinv[i] = k;
  155. // skip if V(k+1:m,k) is empty
  156. if (--nque[k] <= 0) {
  157. continue;
  158. }
  159. // nque[k] is nnz (V(k+1:m,k))
  160. s.lnz += w[nque + k];
  161. // move all rows to parent of k
  162. var pa = parent[k];
  163. if (pa !== -1) {
  164. if (w[nque + pa] === 0) {
  165. w[tail + pa] = w[tail + k];
  166. }
  167. w[next + w[tail + k]] = w[head + pa];
  168. w[head + pa] = w[next + i];
  169. w[nque + pa] += w[nque + k];
  170. }
  171. }
  172. for (i = 0; i < m; i++) {
  173. if (pinv[i] < 0) {
  174. pinv[i] = k++;
  175. }
  176. }
  177. return true;
  178. }
  179. });
  180. exports.createCsSqr = createCsSqr;