usolve.js 4.7 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166
  1. "use strict";
  2. Object.defineProperty(exports, "__esModule", {
  3. value: true
  4. });
  5. exports.createUsolve = void 0;
  6. var _factory = require("../../../utils/factory.js");
  7. var _solveValidation = require("./utils/solveValidation.js");
  8. var name = 'usolve';
  9. var dependencies = ['typed', 'matrix', 'divideScalar', 'multiplyScalar', 'subtract', 'equalScalar', 'DenseMatrix'];
  10. var createUsolve = /* #__PURE__ */(0, _factory.factory)(name, dependencies, function (_ref) {
  11. var typed = _ref.typed,
  12. matrix = _ref.matrix,
  13. divideScalar = _ref.divideScalar,
  14. multiplyScalar = _ref.multiplyScalar,
  15. subtract = _ref.subtract,
  16. equalScalar = _ref.equalScalar,
  17. DenseMatrix = _ref.DenseMatrix;
  18. var solveValidation = (0, _solveValidation.createSolveValidation)({
  19. DenseMatrix: DenseMatrix
  20. });
  21. /**
  22. * Finds one solution of a linear equation system by backward substitution. Matrix must be an upper triangular matrix. Throws an error if there's no solution.
  23. *
  24. * `U * x = b`
  25. *
  26. * Syntax:
  27. *
  28. * math.usolve(U, b)
  29. *
  30. * Examples:
  31. *
  32. * const a = [[-2, 3], [2, 1]]
  33. * const b = [11, 9]
  34. * const x = usolve(a, b) // [[8], [9]]
  35. *
  36. * See also:
  37. *
  38. * usolveAll, lup, slu, usolve, lusolve
  39. *
  40. * @param {Matrix, Array} U A N x N matrix or array (U)
  41. * @param {Matrix, Array} b A column vector with the b values
  42. *
  43. * @return {DenseMatrix | Array} A column vector with the linear system solution (x)
  44. */
  45. return typed(name, {
  46. 'SparseMatrix, Array | Matrix': function SparseMatrixArrayMatrix(m, b) {
  47. return _sparseBackwardSubstitution(m, b);
  48. },
  49. 'DenseMatrix, Array | Matrix': function DenseMatrixArrayMatrix(m, b) {
  50. return _denseBackwardSubstitution(m, b);
  51. },
  52. 'Array, Array | Matrix': function ArrayArrayMatrix(a, b) {
  53. var m = matrix(a);
  54. var r = _denseBackwardSubstitution(m, b);
  55. return r.valueOf();
  56. }
  57. });
  58. function _denseBackwardSubstitution(m, b) {
  59. // make b into a column vector
  60. b = solveValidation(m, b, true);
  61. var bdata = b._data;
  62. var rows = m._size[0];
  63. var columns = m._size[1];
  64. // result
  65. var x = [];
  66. var mdata = m._data;
  67. // loop columns backwards
  68. for (var j = columns - 1; j >= 0; j--) {
  69. // b[j]
  70. var bj = bdata[j][0] || 0;
  71. // x[j]
  72. var xj = void 0;
  73. if (!equalScalar(bj, 0)) {
  74. // value at [j, j]
  75. var vjj = mdata[j][j];
  76. if (equalScalar(vjj, 0)) {
  77. // system cannot be solved
  78. throw new Error('Linear system cannot be solved since matrix is singular');
  79. }
  80. xj = divideScalar(bj, vjj);
  81. // loop rows
  82. for (var i = j - 1; i >= 0; i--) {
  83. // update copy of b
  84. bdata[i] = [subtract(bdata[i][0] || 0, multiplyScalar(xj, mdata[i][j]))];
  85. }
  86. } else {
  87. // zero value at j
  88. xj = 0;
  89. }
  90. // update x
  91. x[j] = [xj];
  92. }
  93. return new DenseMatrix({
  94. data: x,
  95. size: [rows, 1]
  96. });
  97. }
  98. function _sparseBackwardSubstitution(m, b) {
  99. // make b into a column vector
  100. b = solveValidation(m, b, true);
  101. var bdata = b._data;
  102. var rows = m._size[0];
  103. var columns = m._size[1];
  104. var values = m._values;
  105. var index = m._index;
  106. var ptr = m._ptr;
  107. // result
  108. var x = [];
  109. // loop columns backwards
  110. for (var j = columns - 1; j >= 0; j--) {
  111. var bj = bdata[j][0] || 0;
  112. if (!equalScalar(bj, 0)) {
  113. // non-degenerate row, find solution
  114. var vjj = 0;
  115. // upper triangular matrix values & index (column j)
  116. var jValues = [];
  117. var jIndices = [];
  118. // first & last indeces in column
  119. var firstIndex = ptr[j];
  120. var lastIndex = ptr[j + 1];
  121. // values in column, find value at [j, j], loop backwards
  122. for (var k = lastIndex - 1; k >= firstIndex; k--) {
  123. var i = index[k];
  124. // check row (rows are not sorted!)
  125. if (i === j) {
  126. vjj = values[k];
  127. } else if (i < j) {
  128. // store upper triangular
  129. jValues.push(values[k]);
  130. jIndices.push(i);
  131. }
  132. }
  133. // at this point we must have a value in vjj
  134. if (equalScalar(vjj, 0)) {
  135. throw new Error('Linear system cannot be solved since matrix is singular');
  136. }
  137. var xj = divideScalar(bj, vjj);
  138. for (var _k = 0, _lastIndex = jIndices.length; _k < _lastIndex; _k++) {
  139. var _i = jIndices[_k];
  140. bdata[_i] = [subtract(bdata[_i][0], multiplyScalar(xj, jValues[_k]))];
  141. }
  142. x[j] = [xj];
  143. } else {
  144. // degenerate row, we can choose any value
  145. x[j] = [0];
  146. }
  147. }
  148. return new DenseMatrix({
  149. data: x,
  150. size: [rows, 1]
  151. });
  152. }
  153. });
  154. exports.createUsolve = createUsolve;