usolveAll.js (5389B)
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 43 return typed(name, { 44 'SparseMatrix, Array | Matrix': function SparseMatrixArrayMatrix(m, b) { 45 return _sparseBackwardSubstitution(m, b); 46 }, 47 'DenseMatrix, Array | Matrix': function DenseMatrixArrayMatrix(m, b) { 48 return _denseBackwardSubstitution(m, b); 49 }, 50 'Array, Array | Matrix': function ArrayArrayMatrix(a, b) { 51 var m = matrix(a); 52 53 var R = _denseBackwardSubstitution(m, b); 54 55 return R.map(r => r.valueOf()); 56 } 57 }); 58 59 function _denseBackwardSubstitution(m, b_) { 60 // the algorithm is derived from 61 // https://www.overleaf.com/read/csvgqdxggyjv 62 // array of right-hand sides 63 var B = [solveValidation(m, b_, true)._data.map(e => e[0])]; 64 var M = m._data; 65 var rows = m._size[0]; 66 var columns = m._size[1]; // loop columns backwards 67 68 for (var i = columns - 1; i >= 0; i--) { 69 var L = B.length; // loop right-hand sides 70 71 for (var k = 0; k < L; k++) { 72 var b = B[k]; 73 74 if (!equalScalar(M[i][i], 0)) { 75 // non-singular row 76 b[i] = divideScalar(b[i], M[i][i]); 77 78 for (var j = i - 1; j >= 0; j--) { 79 // b[j] -= b[i] * M[j,i] 80 b[j] = subtract(b[j], multiplyScalar(b[i], M[j][i])); 81 } 82 } else if (!equalScalar(b[i], 0)) { 83 // singular row, nonzero RHS 84 if (k === 0) { 85 // There is no valid solution 86 return []; 87 } else { 88 // This RHS is invalid but other solutions may still exist 89 B.splice(k, 1); 90 k -= 1; 91 L -= 1; 92 } 93 } else if (k === 0) { 94 // singular row, RHS is zero 95 var bNew = [...b]; 96 bNew[i] = 1; 97 98 for (var _j = i - 1; _j >= 0; _j--) { 99 bNew[_j] = subtract(bNew[_j], M[_j][i]); 100 } 101 102 B.push(bNew); 103 } 104 } 105 } 106 107 return B.map(x => new DenseMatrix({ 108 data: x.map(e => [e]), 109 size: [rows, 1] 110 })); 111 } 112 113 function _sparseBackwardSubstitution(m, b_) { 114 // array of right-hand sides 115 var B = [solveValidation(m, b_, true)._data.map(e => e[0])]; 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; // loop columns backwards 121 122 for (var i = columns - 1; i >= 0; i--) { 123 var L = B.length; // loop right-hand sides 124 125 for (var k = 0; k < L; k++) { 126 var b = B[k]; // values & indices (column i) 127 128 var iValues = []; 129 var iIndices = []; // first & last indeces in column 130 131 var firstIndex = ptr[i]; 132 var lastIndex = ptr[i + 1]; // find the value at [i, i] 133 134 var Mii = 0; 135 136 for (var j = lastIndex - 1; j >= firstIndex; j--) { 137 var J = index[j]; // check row 138 139 if (J === i) { 140 Mii = values[j]; 141 } else if (J < i) { 142 // store upper triangular 143 iValues.push(values[j]); 144 iIndices.push(J); 145 } 146 } 147 148 if (!equalScalar(Mii, 0)) { 149 // non-singular row 150 b[i] = divideScalar(b[i], Mii); // loop upper triangular 151 152 for (var _j2 = 0, _lastIndex = iIndices.length; _j2 < _lastIndex; _j2++) { 153 var _J = iIndices[_j2]; 154 b[_J] = subtract(b[_J], multiplyScalar(b[i], iValues[_j2])); 155 } 156 } else if (!equalScalar(b[i], 0)) { 157 // singular row, nonzero RHS 158 if (k === 0) { 159 // There is no valid solution 160 return []; 161 } else { 162 // This RHS is invalid but other solutions may still exist 163 B.splice(k, 1); 164 k -= 1; 165 L -= 1; 166 } 167 } else if (k === 0) { 168 // singular row, RHS is zero 169 var bNew = [...b]; 170 bNew[i] = 1; // loop upper triangular 171 172 for (var _j3 = 0, _lastIndex2 = iIndices.length; _j3 < _lastIndex2; _j3++) { 173 var _J2 = iIndices[_j3]; 174 bNew[_J2] = subtract(bNew[_J2], iValues[_j3]); 175 } 176 177 B.push(bNew); 178 } 179 } 180 } 181 182 return B.map(x => new DenseMatrix({ 183 data: x.map(e => [e]), 184 size: [rows, 1] 185 })); 186 } 187 });