1 | import { isArray, isMatrix } from '../../../utils/is'
|
2 | import { factory } from '../../../utils/factory'
|
3 | import { createSolveValidation } from './utils/solveValidation'
|
4 | import { csIpvec } from '../sparse/csIpvec'
|
5 |
|
6 | const name = 'lusolve'
|
7 | const dependencies = [
|
8 | 'typed',
|
9 | 'matrix',
|
10 | 'lup',
|
11 | 'slu',
|
12 | 'usolve',
|
13 | 'lsolve',
|
14 | 'DenseMatrix'
|
15 | ]
|
16 |
|
17 | export const createLusolve = /* #__PURE__ */ factory(name, dependencies, ({ typed, matrix, lup, slu, usolve, lsolve, DenseMatrix }) => {
|
18 | const solveValidation = createSolveValidation({ 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 | return typed(name, {
|
54 |
|
55 | 'Array, Array | Matrix': function (a, b) {
|
56 | // convert a to matrix
|
57 | a = matrix(a)
|
58 | // matrix lup decomposition
|
59 | const d = lup(a)
|
60 | // solve
|
61 | const x = _lusolve(d.L, d.U, d.p, null, b)
|
62 | // convert result to array
|
63 | return x.valueOf()
|
64 | },
|
65 |
|
66 | 'DenseMatrix, Array | Matrix': function (a, b) {
|
67 | // matrix lup decomposition
|
68 | const d = lup(a)
|
69 | // solve
|
70 | return _lusolve(d.L, d.U, d.p, null, b)
|
71 | },
|
72 |
|
73 | 'SparseMatrix, Array | Matrix': function (a, b) {
|
74 | // matrix lup decomposition
|
75 | const d = lup(a)
|
76 | // solve
|
77 | return _lusolve(d.L, d.U, d.p, null, b)
|
78 | },
|
79 |
|
80 | 'SparseMatrix, Array | Matrix, number, number': function (a, b, order, threshold) {
|
81 | // matrix lu decomposition
|
82 | const d = slu(a, order, threshold)
|
83 | // solve
|
84 | return _lusolve(d.L, d.U, d.p, d.q, b)
|
85 | },
|
86 |
|
87 | 'Object, Array | Matrix': function (d, b) {
|
88 | // solve
|
89 | return _lusolve(d.L, d.U, d.p, d.q, b)
|
90 | }
|
91 | })
|
92 |
|
93 | function _toMatrix (a) {
|
94 | // check it is a matrix
|
95 | if (isMatrix(a)) { return a }
|
96 | // check array
|
97 | if (isArray(a)) { return matrix(a) }
|
98 | // throw
|
99 | throw new TypeError('Invalid Matrix LU decomposition')
|
100 | }
|
101 |
|
102 | function _lusolve (l, u, p, q, b) {
|
103 | // verify L, U, P
|
104 | l = _toMatrix(l)
|
105 | u = _toMatrix(u)
|
106 | // validate matrix and vector
|
107 | b = solveValidation(l, b, false)
|
108 | // apply row permutations if needed (b is a DenseMatrix)
|
109 | if (p) { b._data = csIpvec(p, b._data) }
|
110 | // use forward substitution to resolve L * y = b
|
111 | const y = lsolve(l, b)
|
112 | // use backward substitution to resolve U * x = y
|
113 | const x = usolve(u, y)
|
114 | // apply column permutations if needed (x is a DenseMatrix)
|
115 | if (q) { x._data = csIpvec(q, x._data) }
|
116 | // return solution
|
117 | return x
|
118 | }
|
119 | })
|