lsolve.js 4.4 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157
  1. import { factory } from '../../../utils/factory.js';
  2. import { createSolveValidation } from './utils/solveValidation.js';
  3. var name = 'lsolve';
  4. var dependencies = ['typed', 'matrix', 'divideScalar', 'multiplyScalar', 'subtract', 'equalScalar', 'DenseMatrix'];
  5. export var createLsolve = /* #__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 forwards substitution. Matrix must be a lower triangular matrix. Throws an error if there's no solution.
  20. *
  21. * `L * x = b`
  22. *
  23. * Syntax:
  24. *
  25. * math.lsolve(L, b)
  26. *
  27. * Examples:
  28. *
  29. * const a = [[-2, 3], [2, 1]]
  30. * const b = [11, 9]
  31. * const x = lsolve(a, b) // [[-5.5], [20]]
  32. *
  33. * See also:
  34. *
  35. * lsolveAll, 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} A column vector with the linear system solution (x)
  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.valueOf();
  53. }
  54. });
  55. function _denseForwardSubstitution(m, b) {
  56. // validate matrix and vector, return copy of column vector b
  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
  65. for (var j = 0; j < columns; j++) {
  66. var bj = bdata[j][0] || 0;
  67. var xj = void 0;
  68. if (!equalScalar(bj, 0)) {
  69. // non-degenerate row, find solution
  70. var vjj = mdata[j][j];
  71. if (equalScalar(vjj, 0)) {
  72. throw new Error('Linear system cannot be solved since matrix is singular');
  73. }
  74. xj = divideScalar(bj, vjj);
  75. // loop rows
  76. for (var i = j + 1; i < rows; i++) {
  77. bdata[i] = [subtract(bdata[i][0] || 0, multiplyScalar(xj, mdata[i][j]))];
  78. }
  79. } else {
  80. // degenerate row, we can choose any value
  81. xj = 0;
  82. }
  83. x[j] = [xj];
  84. }
  85. return new DenseMatrix({
  86. data: x,
  87. size: [rows, 1]
  88. });
  89. }
  90. function _sparseForwardSubstitution(m, b) {
  91. // validate matrix and vector, return copy of column vector b
  92. b = solveValidation(m, b, true);
  93. var bdata = b._data;
  94. var rows = m._size[0];
  95. var columns = m._size[1];
  96. var values = m._values;
  97. var index = m._index;
  98. var ptr = m._ptr;
  99. // result
  100. var x = [];
  101. // loop columns
  102. for (var j = 0; j < columns; j++) {
  103. var bj = bdata[j][0] || 0;
  104. if (!equalScalar(bj, 0)) {
  105. // non-degenerate row, find solution
  106. var vjj = 0;
  107. // matrix values & indices (column j)
  108. var jValues = [];
  109. var jIndices = [];
  110. // first and last index in the column
  111. var firstIndex = ptr[j];
  112. var lastIndex = ptr[j + 1];
  113. // values in column, find value at [j, j]
  114. for (var k = firstIndex; k < lastIndex; k++) {
  115. var i = index[k];
  116. // check row (rows are not sorted!)
  117. if (i === j) {
  118. vjj = values[k];
  119. } else if (i > j) {
  120. // store lower triangular
  121. jValues.push(values[k]);
  122. jIndices.push(i);
  123. }
  124. }
  125. // at this point we must have a value in vjj
  126. if (equalScalar(vjj, 0)) {
  127. throw new Error('Linear system cannot be solved since matrix is singular');
  128. }
  129. var xj = divideScalar(bj, vjj);
  130. for (var _k = 0, l = jIndices.length; _k < l; _k++) {
  131. var _i = jIndices[_k];
  132. bdata[_i] = [subtract(bdata[_i][0] || 0, multiplyScalar(xj, jValues[_k]))];
  133. }
  134. x[j] = [xj];
  135. } else {
  136. // degenerate row, we can choose any value
  137. x[j] = [0];
  138. }
  139. }
  140. return new DenseMatrix({
  141. data: x,
  142. size: [rows, 1]
  143. });
  144. }
  145. });