csSymperm.js 2.8 KB

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