1 | import { factory } from '../../../utils/factory'
|
2 | import { createSolveValidation } from './utils/solveValidation'
|
3 |
|
4 | const name = 'lsolve'
|
5 | const dependencies = [
|
6 | 'typed',
|
7 | 'matrix',
|
8 | 'divideScalar',
|
9 | 'multiplyScalar',
|
10 | 'subtract',
|
11 | 'equalScalar',
|
12 | 'DenseMatrix'
|
13 | ]
|
14 |
|
15 | export const createLsolve = /* #__PURE__ */ factory(name, dependencies, ({ typed, matrix, divideScalar, multiplyScalar, subtract, equalScalar, DenseMatrix }) => {
|
16 | const solveValidation = createSolveValidation({ DenseMatrix })
|
17 |
|
18 | /**
|
19 | * Solves the linear equation system by forwards substitution. Matrix must be a lower triangular matrix.
|
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 | * 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 |
|
44 | 'SparseMatrix, Array | Matrix': function (m, b) {
|
45 | // process matrix
|
46 | return _sparseForwardSubstitution(m, b)
|
47 | },
|
48 |
|
49 | 'DenseMatrix, Array | Matrix': function (m, b) {
|
50 | // process matrix
|
51 | return _denseForwardSubstitution(m, b)
|
52 | },
|
53 |
|
54 | 'Array, Array | Matrix': function (a, b) {
|
55 | // create dense matrix from array
|
56 | const m = matrix(a)
|
57 | // use matrix implementation
|
58 | const r = _denseForwardSubstitution(m, b)
|
59 | // result
|
60 | return r.valueOf()
|
61 | }
|
62 | })
|
63 |
|
64 | function _denseForwardSubstitution (m, b) {
|
65 | // validate matrix and vector, return copy of column vector b
|
66 | b = solveValidation(m, b, true)
|
67 | // column vector data
|
68 | const bdata = b._data
|
69 | // rows & columns
|
70 | const rows = m._size[0]
|
71 | const columns = m._size[1]
|
72 | // result
|
73 | const x = []
|
74 | // data
|
75 | const data = m._data
|
76 | // forward solve m * x = b, loop columns
|
77 | for (let j = 0; j < columns; j++) {
|
78 | // b[j]
|
79 | const bj = bdata[j][0] || 0
|
80 | // x[j]
|
81 | let xj
|
82 | // forward substitution (outer product) avoids inner looping when bj === 0
|
83 | if (!equalScalar(bj, 0)) {
|
84 | // value @ [j, j]
|
85 | const vjj = data[j][j]
|
86 | // check vjj
|
87 | if (equalScalar(vjj, 0)) {
|
88 | // system cannot be solved
|
89 | throw new Error('Linear system cannot be solved since matrix is singular')
|
90 | }
|
91 | // calculate xj
|
92 | xj = divideScalar(bj, vjj)
|
93 | // loop rows
|
94 | for (let i = j + 1; i < rows; i++) {
|
95 | // update copy of b
|
96 | bdata[i] = [subtract(bdata[i][0] || 0, multiplyScalar(xj, data[i][j]))]
|
97 | }
|
98 | } else {
|
99 | // zero @ j
|
100 | xj = 0
|
101 | }
|
102 | // update x
|
103 | x[j] = [xj]
|
104 | }
|
105 | // return vector
|
106 | return new DenseMatrix({
|
107 | data: x,
|
108 | size: [rows, 1]
|
109 | })
|
110 | }
|
111 |
|
112 | function _sparseForwardSubstitution (m, b) {
|
113 | // validate matrix and vector, return copy of column vector b
|
114 | b = solveValidation(m, b, true)
|
115 | // column vector data
|
116 | const bdata = b._data
|
117 | // rows & columns
|
118 | const rows = m._size[0]
|
119 | const columns = m._size[1]
|
120 | // matrix arrays
|
121 | const values = m._values
|
122 | const index = m._index
|
123 | const ptr = m._ptr
|
124 | // vars
|
125 | let i, k
|
126 | // result
|
127 | const x = []
|
128 | // forward solve m * x = b, loop columns
|
129 | for (let j = 0; j < columns; j++) {
|
130 | // b[j]
|
131 | const bj = bdata[j][0] || 0
|
132 | // forward substitution (outer product) avoids inner looping when bj === 0
|
133 | if (!equalScalar(bj, 0)) {
|
134 | // value @ [j, j]
|
135 | let vjj = 0
|
136 | // lower triangular matrix values & index (column j)
|
137 | const jvalues = []
|
138 | const jindex = []
|
139 | // last index in column
|
140 | let l = ptr[j + 1]
|
141 | // values in column, find value @ [j, j]
|
142 | for (k = ptr[j]; k < l; k++) {
|
143 | // row
|
144 | i = index[k]
|
145 | // check row (rows are not sorted!)
|
146 | if (i === j) {
|
147 | // update vjj
|
148 | vjj = values[k]
|
149 | } else if (i > j) {
|
150 | // store lower triangular
|
151 | jvalues.push(values[k])
|
152 | jindex.push(i)
|
153 | }
|
154 | }
|
155 | // at this point we must have a value @ [j, j]
|
156 | if (equalScalar(vjj, 0)) {
|
157 | // system cannot be solved, there is no value @ [j, j]
|
158 | throw new Error('Linear system cannot be solved since matrix is singular')
|
159 | }
|
160 | // calculate xj
|
161 | const xj = divideScalar(bj, vjj)
|
162 | // loop lower triangular
|
163 | for (k = 0, l = jindex.length; k < l; k++) {
|
164 | // row
|
165 | i = jindex[k]
|
166 | // update copy of b
|
167 | bdata[i] = [subtract(bdata[i][0] || 0, multiplyScalar(xj, jvalues[k]))]
|
168 | }
|
169 | // update x
|
170 | x[j] = [xj]
|
171 | } else {
|
172 | // update x
|
173 | x[j] = [0]
|
174 | }
|
175 | }
|
176 | // return vector
|
177 | return new DenseMatrix({
|
178 | data: x,
|
179 | size: [rows, 1]
|
180 | })
|
181 | }
|
182 | })
|