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