csSpsolve.js (3078B)
1 "use strict"; 2 3 Object.defineProperty(exports, "__esModule", { 4 value: true 5 }); 6 exports.createCsSpsolve = void 0; 7 8 var _csReach = require("./csReach.js"); 9 10 var _factory = require("../../../utils/factory.js"); 11 12 var name = 'csSpsolve'; 13 var dependencies = ['divideScalar', 'multiply', 'subtract']; 14 var createCsSpsolve = /* #__PURE__ */(0, _factory.factory)(name, dependencies, function (_ref) { 15 var divideScalar = _ref.divideScalar, 16 multiply = _ref.multiply, 17 subtract = _ref.subtract; 18 19 /** 20 * The function csSpsolve() computes the solution to G * x = bk, where bk is the 21 * kth column of B. When lo is true, the function assumes G = L is lower triangular with the 22 * diagonal entry as the first entry in each column. When lo is true, the function assumes G = U 23 * is upper triangular with the diagonal entry as the last entry in each column. 24 * 25 * @param {Matrix} g The G matrix 26 * @param {Matrix} b The B matrix 27 * @param {Number} k The kth column in B 28 * @param {Array} xi The nonzero pattern xi[top] .. xi[n - 1], an array of size = 2 * n 29 * The first n entries is the nonzero pattern, the last n entries is the stack 30 * @param {Array} x The soluton to the linear system G * x = b 31 * @param {Array} pinv The inverse row permutation vector, must be null for L * x = b 32 * @param {boolean} lo The lower (true) upper triangular (false) flag 33 * 34 * @return {Number} The index for the nonzero pattern 35 * 36 * Reference: http://faculty.cse.tamu.edu/davis/publications.html 37 */ 38 return function csSpsolve(g, b, k, xi, x, pinv, lo) { 39 // g arrays 40 var gvalues = g._values; 41 var gindex = g._index; 42 var gptr = g._ptr; 43 var gsize = g._size; // columns 44 45 var n = gsize[1]; // b arrays 46 47 var bvalues = b._values; 48 var bindex = b._index; 49 var bptr = b._ptr; // vars 50 51 var p, p0, p1, q; // xi[top..n-1] = csReach(B(:,k)) 52 53 var top = (0, _csReach.csReach)(g, b, k, xi, pinv); // clear x 54 55 for (p = top; p < n; p++) { 56 x[xi[p]] = 0; 57 } // scatter b 58 59 60 for (p0 = bptr[k], p1 = bptr[k + 1], p = p0; p < p1; p++) { 61 x[bindex[p]] = bvalues[p]; 62 } // loop columns 63 64 65 for (var px = top; px < n; px++) { 66 // x array index for px 67 var j = xi[px]; // apply permutation vector (U x = b), j maps to column J of G 68 69 var J = pinv ? pinv[j] : j; // check column J is empty 70 71 if (J < 0) { 72 continue; 73 } // column value indeces in G, p0 <= p < p1 74 75 76 p0 = gptr[J]; 77 p1 = gptr[J + 1]; // x(j) /= G(j,j) 78 79 x[j] = divideScalar(x[j], gvalues[lo ? p0 : p1 - 1]); // first entry L(j,j) 80 81 p = lo ? p0 + 1 : p0; 82 q = lo ? p1 : p1 - 1; // loop 83 84 for (; p < q; p++) { 85 // row 86 var i = gindex[p]; // x(i) -= G(i,j) * x(j) 87 88 x[i] = subtract(x[i], multiply(gvalues[p], x[j])); 89 } 90 } // return top of stack 91 92 93 return top; 94 }; 95 }); 96 exports.createCsSpsolve = createCsSpsolve;