lup.js (10203B)
1 import { clone } from '../../../utils/object.js'; 2 import { factory } from '../../../utils/factory.js'; 3 var name = 'lup'; 4 var dependencies = ['typed', 'matrix', 'abs', 'addScalar', 'divideScalar', 'multiplyScalar', 'subtract', 'larger', 'equalScalar', 'unaryMinus', 'DenseMatrix', 'SparseMatrix', 'Spa']; 5 export var createLup = /* #__PURE__ */factory(name, dependencies, _ref => { 6 var { 7 typed, 8 matrix, 9 abs, 10 addScalar, 11 divideScalar, 12 multiplyScalar, 13 subtract, 14 larger, 15 equalScalar, 16 unaryMinus, 17 DenseMatrix, 18 SparseMatrix, 19 Spa 20 } = _ref; 21 22 /** 23 * Calculate the Matrix LU decomposition with partial pivoting. Matrix `A` is decomposed in two matrices (`L`, `U`) and a 24 * row permutation vector `p` where `A[p,:] = L * U` 25 * 26 * Syntax: 27 * 28 * math.lup(A) 29 * 30 * Example: 31 * 32 * const m = [[2, 1], [1, 4]] 33 * const r = math.lup(m) 34 * // r = { 35 * // L: [[1, 0], [0.5, 1]], 36 * // U: [[2, 1], [0, 3.5]], 37 * // P: [0, 1] 38 * // } 39 * 40 * See also: 41 * 42 * slu, lsolve, lusolve, usolve 43 * 44 * @param {Matrix | Array} A A two dimensional matrix or array for which to get the LUP decomposition. 45 * 46 * @return {{L: Array | Matrix, U: Array | Matrix, P: Array.<number>}} The lower triangular matrix, the upper triangular matrix and the permutation matrix. 47 */ 48 return typed(name, { 49 DenseMatrix: function DenseMatrix(m) { 50 return _denseLUP(m); 51 }, 52 SparseMatrix: function SparseMatrix(m) { 53 return _sparseLUP(m); 54 }, 55 Array: function Array(a) { 56 // create dense matrix from array 57 var m = matrix(a); // lup, use matrix implementation 58 59 var r = _denseLUP(m); // result 60 61 62 return { 63 L: r.L.valueOf(), 64 U: r.U.valueOf(), 65 p: r.p 66 }; 67 } 68 }); 69 70 function _denseLUP(m) { 71 // rows & columns 72 var rows = m._size[0]; 73 var columns = m._size[1]; // minimum rows and columns 74 75 var n = Math.min(rows, columns); // matrix array, clone original data 76 77 var data = clone(m._data); // l matrix arrays 78 79 var ldata = []; 80 var lsize = [rows, n]; // u matrix arrays 81 82 var udata = []; 83 var usize = [n, columns]; // vars 84 85 var i, j, k; // permutation vector 86 87 var p = []; 88 89 for (i = 0; i < rows; i++) { 90 p[i] = i; 91 } // loop columns 92 93 94 for (j = 0; j < columns; j++) { 95 // skip first column in upper triangular matrix 96 if (j > 0) { 97 // loop rows 98 for (i = 0; i < rows; i++) { 99 // min i,j 100 var min = Math.min(i, j); // v[i, j] 101 102 var s = 0; // loop up to min 103 104 for (k = 0; k < min; k++) { 105 // s = l[i, k] - data[k, j] 106 s = addScalar(s, multiplyScalar(data[i][k], data[k][j])); 107 } 108 109 data[i][j] = subtract(data[i][j], s); 110 } 111 } // row with larger value in cvector, row >= j 112 113 114 var pi = j; 115 var pabsv = 0; 116 var vjj = 0; // loop rows 117 118 for (i = j; i < rows; i++) { 119 // data @ i, j 120 var v = data[i][j]; // absolute value 121 122 var absv = abs(v); // value is greater than pivote value 123 124 if (larger(absv, pabsv)) { 125 // store row 126 pi = i; // update max value 127 128 pabsv = absv; // value @ [j, j] 129 130 vjj = v; 131 } 132 } // swap rows (j <-> pi) 133 134 135 if (j !== pi) { 136 // swap values j <-> pi in p 137 p[j] = [p[pi], p[pi] = p[j]][0]; // swap j <-> pi in data 138 139 DenseMatrix._swapRows(j, pi, data); 140 } // check column is in lower triangular matrix 141 142 143 if (j < rows) { 144 // loop rows (lower triangular matrix) 145 for (i = j + 1; i < rows; i++) { 146 // value @ i, j 147 var vij = data[i][j]; 148 149 if (!equalScalar(vij, 0)) { 150 // update data 151 data[i][j] = divideScalar(data[i][j], vjj); 152 } 153 } 154 } 155 } // loop columns 156 157 158 for (j = 0; j < columns; j++) { 159 // loop rows 160 for (i = 0; i < rows; i++) { 161 // initialize row in arrays 162 if (j === 0) { 163 // check row exists in upper triangular matrix 164 if (i < columns) { 165 // U 166 udata[i] = []; 167 } // L 168 169 170 ldata[i] = []; 171 } // check we are in the upper triangular matrix 172 173 174 if (i < j) { 175 // check row exists in upper triangular matrix 176 if (i < columns) { 177 // U 178 udata[i][j] = data[i][j]; 179 } // check column exists in lower triangular matrix 180 181 182 if (j < rows) { 183 // L 184 ldata[i][j] = 0; 185 } 186 187 continue; 188 } // diagonal value 189 190 191 if (i === j) { 192 // check row exists in upper triangular matrix 193 if (i < columns) { 194 // U 195 udata[i][j] = data[i][j]; 196 } // check column exists in lower triangular matrix 197 198 199 if (j < rows) { 200 // L 201 ldata[i][j] = 1; 202 } 203 204 continue; 205 } // check row exists in upper triangular matrix 206 207 208 if (i < columns) { 209 // U 210 udata[i][j] = 0; 211 } // check column exists in lower triangular matrix 212 213 214 if (j < rows) { 215 // L 216 ldata[i][j] = data[i][j]; 217 } 218 } 219 } // l matrix 220 221 222 var l = new DenseMatrix({ 223 data: ldata, 224 size: lsize 225 }); // u matrix 226 227 var u = new DenseMatrix({ 228 data: udata, 229 size: usize 230 }); // p vector 231 232 var pv = []; 233 234 for (i = 0, n = p.length; i < n; i++) { 235 pv[p[i]] = i; 236 } // return matrices 237 238 239 return { 240 L: l, 241 U: u, 242 p: pv, 243 toString: function toString() { 244 return 'L: ' + this.L.toString() + '\nU: ' + this.U.toString() + '\nP: ' + this.p; 245 } 246 }; 247 } 248 249 function _sparseLUP(m) { 250 // rows & columns 251 var rows = m._size[0]; 252 var columns = m._size[1]; // minimum rows and columns 253 254 var n = Math.min(rows, columns); // matrix arrays (will not be modified, thanks to permutation vector) 255 256 var values = m._values; 257 var index = m._index; 258 var ptr = m._ptr; // l matrix arrays 259 260 var lvalues = []; 261 var lindex = []; 262 var lptr = []; 263 var lsize = [rows, n]; // u matrix arrays 264 265 var uvalues = []; 266 var uindex = []; 267 var uptr = []; 268 var usize = [n, columns]; // vars 269 270 var i, j, k; // permutation vectors, (current index -> original index) and (original index -> current index) 271 272 var pvCo = []; 273 var pvOc = []; 274 275 for (i = 0; i < rows; i++) { 276 pvCo[i] = i; 277 pvOc[i] = i; 278 } // swap indices in permutation vectors (condition x < y)! 279 280 281 var swapIndeces = function swapIndeces(x, y) { 282 // find pv indeces getting data from x and y 283 var kx = pvOc[x]; 284 var ky = pvOc[y]; // update permutation vector current -> original 285 286 pvCo[kx] = y; 287 pvCo[ky] = x; // update permutation vector original -> current 288 289 pvOc[x] = ky; 290 pvOc[y] = kx; 291 }; // loop columns 292 293 294 var _loop = function _loop() { 295 // sparse accumulator 296 var spa = new Spa(); // check lower triangular matrix has a value @ column j 297 298 if (j < rows) { 299 // update ptr 300 lptr.push(lvalues.length); // first value in j column for lower triangular matrix 301 302 lvalues.push(1); 303 lindex.push(j); 304 } // update ptr 305 306 307 uptr.push(uvalues.length); // k0 <= k < k1 where k0 = _ptr[j] && k1 = _ptr[j+1] 308 309 var k0 = ptr[j]; 310 var k1 = ptr[j + 1]; // copy column j into sparse accumulator 311 312 for (k = k0; k < k1; k++) { 313 // row 314 i = index[k]; // copy column values into sparse accumulator (use permutation vector) 315 316 spa.set(pvCo[i], values[k]); 317 } // skip first column in upper triangular matrix 318 319 320 if (j > 0) { 321 // loop rows in column j (above diagonal) 322 spa.forEach(0, j - 1, function (k, vkj) { 323 // loop rows in column k (L) 324 SparseMatrix._forEachRow(k, lvalues, lindex, lptr, function (i, vik) { 325 // check row is below k 326 if (i > k) { 327 // update spa value 328 spa.accumulate(i, unaryMinus(multiplyScalar(vik, vkj))); 329 } 330 }); 331 }); 332 } // row with larger value in spa, row >= j 333 334 335 var pi = j; 336 var vjj = spa.get(j); 337 var pabsv = abs(vjj); // loop values in spa (order by row, below diagonal) 338 339 spa.forEach(j + 1, rows - 1, function (x, v) { 340 // absolute value 341 var absv = abs(v); // value is greater than pivote value 342 343 if (larger(absv, pabsv)) { 344 // store row 345 pi = x; // update max value 346 347 pabsv = absv; // value @ [j, j] 348 349 vjj = v; 350 } 351 }); // swap rows (j <-> pi) 352 353 if (j !== pi) { 354 // swap values j <-> pi in L 355 SparseMatrix._swapRows(j, pi, lsize[1], lvalues, lindex, lptr); // swap values j <-> pi in U 356 357 358 SparseMatrix._swapRows(j, pi, usize[1], uvalues, uindex, uptr); // swap values in spa 359 360 361 spa.swap(j, pi); // update permutation vector (swap values @ j, pi) 362 363 swapIndeces(j, pi); 364 } // loop values in spa (order by row) 365 366 367 spa.forEach(0, rows - 1, function (x, v) { 368 // check we are above diagonal 369 if (x <= j) { 370 // update upper triangular matrix 371 uvalues.push(v); 372 uindex.push(x); 373 } else { 374 // update value 375 v = divideScalar(v, vjj); // check value is non zero 376 377 if (!equalScalar(v, 0)) { 378 // update lower triangular matrix 379 lvalues.push(v); 380 lindex.push(x); 381 } 382 } 383 }); 384 }; 385 386 for (j = 0; j < columns; j++) { 387 _loop(); 388 } // update ptrs 389 390 391 uptr.push(uvalues.length); 392 lptr.push(lvalues.length); // return matrices 393 394 return { 395 L: new SparseMatrix({ 396 values: lvalues, 397 index: lindex, 398 ptr: lptr, 399 size: lsize 400 }), 401 U: new SparseMatrix({ 402 values: uvalues, 403 index: uindex, 404 ptr: uptr, 405 size: usize 406 }), 407 p: pvCo, 408 toString: function toString() { 409 return 'L: ' + this.L.toString() + '\nU: ' + this.U.toString() + '\nP: ' + this.p; 410 } 411 }; 412 } 413 });