usolveAll.js 6.1 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211
  1. "use strict";
  2. var _interopRequireDefault = require("@babel/runtime/helpers/interopRequireDefault");
  3. Object.defineProperty(exports, "__esModule", {
  4. value: true
  5. });
  6. exports.createUsolveAll = void 0;
  7. var _toConsumableArray2 = _interopRequireDefault(require("@babel/runtime/helpers/toConsumableArray"));
  8. var _factory = require("../../../utils/factory.js");
  9. var _solveValidation = require("./utils/solveValidation.js");
  10. var name = 'usolveAll';
  11. var dependencies = ['typed', 'matrix', 'divideScalar', 'multiplyScalar', 'subtract', 'equalScalar', 'DenseMatrix'];
  12. var createUsolveAll = /* #__PURE__ */(0, _factory.factory)(name, dependencies, function (_ref) {
  13. var typed = _ref.typed,
  14. matrix = _ref.matrix,
  15. divideScalar = _ref.divideScalar,
  16. multiplyScalar = _ref.multiplyScalar,
  17. subtract = _ref.subtract,
  18. equalScalar = _ref.equalScalar,
  19. DenseMatrix = _ref.DenseMatrix;
  20. var solveValidation = (0, _solveValidation.createSolveValidation)({
  21. DenseMatrix: DenseMatrix
  22. });
  23. /**
  24. * Finds all solutions of a linear equation system by backward substitution. Matrix must be an upper triangular matrix.
  25. *
  26. * `U * x = b`
  27. *
  28. * Syntax:
  29. *
  30. * math.usolveAll(U, b)
  31. *
  32. * Examples:
  33. *
  34. * const a = [[-2, 3], [2, 1]]
  35. * const b = [11, 9]
  36. * const x = usolveAll(a, b) // [ [[8], [9]] ]
  37. *
  38. * See also:
  39. *
  40. * usolve, lup, slu, usolve, lusolve
  41. *
  42. * @param {Matrix, Array} U A N x N matrix or array (U)
  43. * @param {Matrix, Array} b A column vector with the b values
  44. *
  45. * @return {DenseMatrix[] | Array[]} An array of affine-independent column vectors (x) that solve the linear system
  46. */
  47. return typed(name, {
  48. 'SparseMatrix, Array | Matrix': function SparseMatrixArrayMatrix(m, b) {
  49. return _sparseBackwardSubstitution(m, b);
  50. },
  51. 'DenseMatrix, Array | Matrix': function DenseMatrixArrayMatrix(m, b) {
  52. return _denseBackwardSubstitution(m, b);
  53. },
  54. 'Array, Array | Matrix': function ArrayArrayMatrix(a, b) {
  55. var m = matrix(a);
  56. var R = _denseBackwardSubstitution(m, b);
  57. return R.map(function (r) {
  58. return r.valueOf();
  59. });
  60. }
  61. });
  62. function _denseBackwardSubstitution(m, b_) {
  63. // the algorithm is derived from
  64. // https://www.overleaf.com/read/csvgqdxggyjv
  65. // array of right-hand sides
  66. var B = [solveValidation(m, b_, true)._data.map(function (e) {
  67. return e[0];
  68. })];
  69. var M = m._data;
  70. var rows = m._size[0];
  71. var columns = m._size[1];
  72. // loop columns backwards
  73. for (var i = columns - 1; i >= 0; i--) {
  74. var L = B.length;
  75. // loop right-hand sides
  76. for (var k = 0; k < L; k++) {
  77. var b = B[k];
  78. if (!equalScalar(M[i][i], 0)) {
  79. // non-singular row
  80. b[i] = divideScalar(b[i], M[i][i]);
  81. for (var j = i - 1; j >= 0; j--) {
  82. // b[j] -= b[i] * M[j,i]
  83. b[j] = subtract(b[j], multiplyScalar(b[i], M[j][i]));
  84. }
  85. } else if (!equalScalar(b[i], 0)) {
  86. // singular row, nonzero RHS
  87. if (k === 0) {
  88. // There is no valid solution
  89. return [];
  90. } else {
  91. // This RHS is invalid but other solutions may still exist
  92. B.splice(k, 1);
  93. k -= 1;
  94. L -= 1;
  95. }
  96. } else if (k === 0) {
  97. // singular row, RHS is zero
  98. var bNew = (0, _toConsumableArray2["default"])(b);
  99. bNew[i] = 1;
  100. for (var _j = i - 1; _j >= 0; _j--) {
  101. bNew[_j] = subtract(bNew[_j], M[_j][i]);
  102. }
  103. B.push(bNew);
  104. }
  105. }
  106. }
  107. return B.map(function (x) {
  108. return new DenseMatrix({
  109. data: x.map(function (e) {
  110. return [e];
  111. }),
  112. size: [rows, 1]
  113. });
  114. });
  115. }
  116. function _sparseBackwardSubstitution(m, b_) {
  117. // array of right-hand sides
  118. var B = [solveValidation(m, b_, true)._data.map(function (e) {
  119. return e[0];
  120. })];
  121. var rows = m._size[0];
  122. var columns = m._size[1];
  123. var values = m._values;
  124. var index = m._index;
  125. var ptr = m._ptr;
  126. // loop columns backwards
  127. for (var i = columns - 1; i >= 0; i--) {
  128. var L = B.length;
  129. // loop right-hand sides
  130. for (var k = 0; k < L; k++) {
  131. var b = B[k];
  132. // values & indices (column i)
  133. var iValues = [];
  134. var iIndices = [];
  135. // first & last indeces in column
  136. var firstIndex = ptr[i];
  137. var lastIndex = ptr[i + 1];
  138. // find the value at [i, i]
  139. var Mii = 0;
  140. for (var j = lastIndex - 1; j >= firstIndex; j--) {
  141. var J = index[j];
  142. // check row
  143. if (J === i) {
  144. Mii = values[j];
  145. } else if (J < i) {
  146. // store upper triangular
  147. iValues.push(values[j]);
  148. iIndices.push(J);
  149. }
  150. }
  151. if (!equalScalar(Mii, 0)) {
  152. // non-singular row
  153. b[i] = divideScalar(b[i], Mii);
  154. // loop upper triangular
  155. for (var _j2 = 0, _lastIndex = iIndices.length; _j2 < _lastIndex; _j2++) {
  156. var _J = iIndices[_j2];
  157. b[_J] = subtract(b[_J], multiplyScalar(b[i], iValues[_j2]));
  158. }
  159. } else if (!equalScalar(b[i], 0)) {
  160. // singular row, nonzero RHS
  161. if (k === 0) {
  162. // There is no valid solution
  163. return [];
  164. } else {
  165. // This RHS is invalid but other solutions may still exist
  166. B.splice(k, 1);
  167. k -= 1;
  168. L -= 1;
  169. }
  170. } else if (k === 0) {
  171. // singular row, RHS is zero
  172. var bNew = (0, _toConsumableArray2["default"])(b);
  173. bNew[i] = 1;
  174. // loop upper triangular
  175. for (var _j3 = 0, _lastIndex2 = iIndices.length; _j3 < _lastIndex2; _j3++) {
  176. var _J2 = iIndices[_j3];
  177. bNew[_J2] = subtract(bNew[_J2], iValues[_j3]);
  178. }
  179. B.push(bNew);
  180. }
  181. }
  182. }
  183. return B.map(function (x) {
  184. return new DenseMatrix({
  185. data: x.map(function (e) {
  186. return [e];
  187. }),
  188. size: [rows, 1]
  189. });
  190. });
  191. }
  192. });
  193. exports.createUsolveAll = createUsolveAll;