lsolveAll.js 5.2 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186
  1. import { factory } from '../../../utils/factory.js';
  2. import { createSolveValidation } from './utils/solveValidation.js';
  3. var name = 'lsolveAll';
  4. var dependencies = ['typed', 'matrix', 'divideScalar', 'multiplyScalar', 'subtract', 'equalScalar', 'DenseMatrix'];
  5. export var createLsolveAll = /* #__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 forwards substitution. Matrix must be a lower triangular matrix.
  20. *
  21. * `L * x = b`
  22. *
  23. * Syntax:
  24. *
  25. * math.lsolveAll(L, b)
  26. *
  27. * Examples:
  28. *
  29. * const a = [[-2, 3], [2, 1]]
  30. * const b = [11, 9]
  31. * const x = lsolveAll(a, b) // [ [[-5.5], [20]] ]
  32. *
  33. * See also:
  34. *
  35. * lsolve, lup, slu, usolve, lusolve
  36. *
  37. * @param {Matrix, Array} L A N x N matrix or array (L)
  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 _sparseForwardSubstitution(m, b);
  45. },
  46. 'DenseMatrix, Array | Matrix': function DenseMatrixArrayMatrix(m, b) {
  47. return _denseForwardSubstitution(m, b);
  48. },
  49. 'Array, Array | Matrix': function ArrayArrayMatrix(a, b) {
  50. var m = matrix(a);
  51. var R = _denseForwardSubstitution(m, b);
  52. return R.map(r => r.valueOf());
  53. }
  54. });
  55. function _denseForwardSubstitution(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
  64. for (var i = 0; i < columns; 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 < columns; 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 < columns; _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 _sparseForwardSubstitution(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
  112. for (var i = 0; i < columns; 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 = firstIndex; j < lastIndex; j++) {
  126. var J = index[j];
  127. // check row
  128. if (J === i) {
  129. Mii = values[j];
  130. } else if (J > i) {
  131. // store lower 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. for (var _j2 = 0, _lastIndex = iIndices.length; _j2 < _lastIndex; _j2++) {
  140. var _J = iIndices[_j2];
  141. b[_J] = subtract(b[_J], multiplyScalar(b[i], iValues[_j2]));
  142. }
  143. } else if (!equalScalar(b[i], 0)) {
  144. // singular row, nonzero RHS
  145. if (k === 0) {
  146. // There is no valid solution
  147. return [];
  148. } else {
  149. // This RHS is invalid but other solutions may still exist
  150. B.splice(k, 1);
  151. k -= 1;
  152. L -= 1;
  153. }
  154. } else if (k === 0) {
  155. // singular row, RHS is zero
  156. var bNew = [...b];
  157. bNew[i] = 1;
  158. for (var _j3 = 0, _lastIndex2 = iIndices.length; _j3 < _lastIndex2; _j3++) {
  159. var _J2 = iIndices[_j3];
  160. bNew[_J2] = subtract(bNew[_J2], iValues[_j3]);
  161. }
  162. B.push(bNew);
  163. }
  164. }
  165. }
  166. return B.map(x => new DenseMatrix({
  167. data: x.map(e => [e]),
  168. size: [rows, 1]
  169. }));
  170. }
  171. });