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