usolveAll.js 5.3 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190
  1. import { factory } from '../../../utils/factory.js';
  2. import { createSolveValidation } from './utils/solveValidation.js';
  3. var name = 'usolveAll';
  4. var dependencies = ['typed', 'matrix', 'divideScalar', 'multiplyScalar', 'subtract', 'equalScalar', 'DenseMatrix'];
  5. export var createUsolveAll = /* #__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 all solutions of a linear equation system by backward substitution. Matrix must be an upper triangular matrix.
  20. *
  21. * `U * x = b`
  22. *
  23. * Syntax:
  24. *
  25. * math.usolveAll(U, b)
  26. *
  27. * Examples:
  28. *
  29. * const a = [[-2, 3], [2, 1]]
  30. * const b = [11, 9]
  31. * const x = usolveAll(a, b) // [ [[8], [9]] ]
  32. *
  33. * See also:
  34. *
  35. * usolve, 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[]} An array of affine-independent column vectors (x) that solve the linear system
  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.map(r => r.valueOf());
  53. }
  54. });
  55. function _denseBackwardSubstitution(m, b_) {
  56. // the algorithm is derived from
  57. // https://www.overleaf.com/read/csvgqdxggyjv
  58. // array of right-hand sides
  59. var B = [solveValidation(m, b_, true)._data.map(e => e[0])];
  60. var M = m._data;
  61. var rows = m._size[0];
  62. var columns = m._size[1];
  63. // loop columns backwards
  64. for (var i = columns - 1; i >= 0; i--) {
  65. var L = B.length;
  66. // loop right-hand sides
  67. for (var k = 0; k < L; k++) {
  68. var b = B[k];
  69. if (!equalScalar(M[i][i], 0)) {
  70. // non-singular row
  71. b[i] = divideScalar(b[i], M[i][i]);
  72. for (var j = i - 1; j >= 0; j--) {
  73. // b[j] -= b[i] * M[j,i]
  74. b[j] = subtract(b[j], multiplyScalar(b[i], M[j][i]));
  75. }
  76. } else if (!equalScalar(b[i], 0)) {
  77. // singular row, nonzero RHS
  78. if (k === 0) {
  79. // There is no valid solution
  80. return [];
  81. } else {
  82. // This RHS is invalid but other solutions may still exist
  83. B.splice(k, 1);
  84. k -= 1;
  85. L -= 1;
  86. }
  87. } else if (k === 0) {
  88. // singular row, RHS is zero
  89. var bNew = [...b];
  90. bNew[i] = 1;
  91. for (var _j = i - 1; _j >= 0; _j--) {
  92. bNew[_j] = subtract(bNew[_j], M[_j][i]);
  93. }
  94. B.push(bNew);
  95. }
  96. }
  97. }
  98. return B.map(x => new DenseMatrix({
  99. data: x.map(e => [e]),
  100. size: [rows, 1]
  101. }));
  102. }
  103. function _sparseBackwardSubstitution(m, b_) {
  104. // array of right-hand sides
  105. var B = [solveValidation(m, b_, true)._data.map(e => e[0])];
  106. var rows = m._size[0];
  107. var columns = m._size[1];
  108. var values = m._values;
  109. var index = m._index;
  110. var ptr = m._ptr;
  111. // loop columns backwards
  112. for (var i = columns - 1; i >= 0; i--) {
  113. var L = B.length;
  114. // loop right-hand sides
  115. for (var k = 0; k < L; k++) {
  116. var b = B[k];
  117. // values & indices (column i)
  118. var iValues = [];
  119. var iIndices = [];
  120. // first & last indeces in column
  121. var firstIndex = ptr[i];
  122. var lastIndex = ptr[i + 1];
  123. // find the value at [i, i]
  124. var Mii = 0;
  125. for (var j = lastIndex - 1; j >= firstIndex; j--) {
  126. var J = index[j];
  127. // check row
  128. if (J === i) {
  129. Mii = values[j];
  130. } else if (J < i) {
  131. // store upper triangular
  132. iValues.push(values[j]);
  133. iIndices.push(J);
  134. }
  135. }
  136. if (!equalScalar(Mii, 0)) {
  137. // non-singular row
  138. b[i] = divideScalar(b[i], Mii);
  139. // loop upper triangular
  140. for (var _j2 = 0, _lastIndex = iIndices.length; _j2 < _lastIndex; _j2++) {
  141. var _J = iIndices[_j2];
  142. b[_J] = subtract(b[_J], multiplyScalar(b[i], iValues[_j2]));
  143. }
  144. } else if (!equalScalar(b[i], 0)) {
  145. // singular row, nonzero RHS
  146. if (k === 0) {
  147. // There is no valid solution
  148. return [];
  149. } else {
  150. // This RHS is invalid but other solutions may still exist
  151. B.splice(k, 1);
  152. k -= 1;
  153. L -= 1;
  154. }
  155. } else if (k === 0) {
  156. // singular row, RHS is zero
  157. var bNew = [...b];
  158. bNew[i] = 1;
  159. // loop upper triangular
  160. for (var _j3 = 0, _lastIndex2 = iIndices.length; _j3 < _lastIndex2; _j3++) {
  161. var _J2 = iIndices[_j3];
  162. bNew[_J2] = subtract(bNew[_J2], iValues[_j3]);
  163. }
  164. B.push(bNew);
  165. }
  166. }
  167. }
  168. return B.map(x => new DenseMatrix({
  169. data: x.map(e => [e]),
  170. size: [rows, 1]
  171. }));
  172. }
  173. });