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