lusolve.js (3727B)
1 import { isArray, isMatrix } from '../../../utils/is.js'; 2 import { factory } from '../../../utils/factory.js'; 3 import { createSolveValidation } from './utils/solveValidation.js'; 4 import { csIpvec } from '../sparse/csIpvec.js'; 5 var name = 'lusolve'; 6 var dependencies = ['typed', 'matrix', 'lup', 'slu', 'usolve', 'lsolve', 'DenseMatrix']; 7 export var createLusolve = /* #__PURE__ */factory(name, dependencies, _ref => { 8 var { 9 typed, 10 matrix, 11 lup, 12 slu, 13 usolve, 14 lsolve, 15 DenseMatrix 16 } = _ref; 17 var solveValidation = createSolveValidation({ 18 DenseMatrix 19 }); 20 /** 21 * Solves the linear system `A * x = b` where `A` is an [n x n] matrix and `b` is a [n] column vector. 22 * 23 * Syntax: 24 * 25 * math.lusolve(A, b) // returns column vector with the solution to the linear system A * x = b 26 * math.lusolve(lup, b) // returns column vector with the solution to the linear system A * x = b, lup = math.lup(A) 27 * 28 * Examples: 29 * 30 * const m = [[1, 0, 0, 0], [0, 2, 0, 0], [0, 0, 3, 0], [0, 0, 0, 4]] 31 * 32 * const x = math.lusolve(m, [-1, -1, -1, -1]) // x = [[-1], [-0.5], [-1/3], [-0.25]] 33 * 34 * const f = math.lup(m) 35 * const x1 = math.lusolve(f, [-1, -1, -1, -1]) // x1 = [[-1], [-0.5], [-1/3], [-0.25]] 36 * const x2 = math.lusolve(f, [1, 2, 1, -1]) // x2 = [[1], [1], [1/3], [-0.25]] 37 * 38 * const a = [[-2, 3], [2, 1]] 39 * const b = [11, 9] 40 * const x = math.lusolve(a, b) // [[2], [5]] 41 * 42 * See also: 43 * 44 * lup, slu, lsolve, usolve 45 * 46 * @param {Matrix | Array | Object} A Invertible Matrix or the Matrix LU decomposition 47 * @param {Matrix | Array} b Column Vector 48 * @param {number} [order] The Symbolic Ordering and Analysis order, see slu for details. Matrix must be a SparseMatrix 49 * @param {Number} [threshold] Partial pivoting threshold (1 for partial pivoting), see slu for details. Matrix must be a SparseMatrix. 50 * 51 * @return {DenseMatrix | Array} Column vector with the solution to the linear system A * x = b 52 */ 53 54 return typed(name, { 55 'Array, Array | Matrix': function ArrayArrayMatrix(a, b) { 56 a = matrix(a); 57 var d = lup(a); 58 59 var x = _lusolve(d.L, d.U, d.p, null, b); 60 61 return x.valueOf(); 62 }, 63 'DenseMatrix, Array | Matrix': function DenseMatrixArrayMatrix(a, b) { 64 var d = lup(a); 65 return _lusolve(d.L, d.U, d.p, null, b); 66 }, 67 'SparseMatrix, Array | Matrix': function SparseMatrixArrayMatrix(a, b) { 68 var d = lup(a); 69 return _lusolve(d.L, d.U, d.p, null, b); 70 }, 71 'SparseMatrix, Array | Matrix, number, number': function SparseMatrixArrayMatrixNumberNumber(a, b, order, threshold) { 72 var d = slu(a, order, threshold); 73 return _lusolve(d.L, d.U, d.p, d.q, b); 74 }, 75 'Object, Array | Matrix': function ObjectArrayMatrix(d, b) { 76 return _lusolve(d.L, d.U, d.p, d.q, b); 77 } 78 }); 79 80 function _toMatrix(a) { 81 if (isMatrix(a)) { 82 return a; 83 } 84 85 if (isArray(a)) { 86 return matrix(a); 87 } 88 89 throw new TypeError('Invalid Matrix LU decomposition'); 90 } 91 92 function _lusolve(l, u, p, q, b) { 93 // verify decomposition 94 l = _toMatrix(l); 95 u = _toMatrix(u); // apply row permutations if needed (b is a DenseMatrix) 96 97 if (p) { 98 b = solveValidation(l, b, true); 99 b._data = csIpvec(p, b._data); 100 } // use forward substitution to resolve L * y = b 101 102 103 var y = lsolve(l, b); // use backward substitution to resolve U * x = y 104 105 var x = usolve(u, y); // apply column permutations if needed (x is a DenseMatrix) 106 107 if (q) { 108 x._data = csIpvec(q, x._data); 109 } 110 111 return x; 112 } 113 });