usolve.js 4.4 KB

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