simple-squiggle

A restricted subset of Squiggle
Log | Files | Refs | README

csReach.js (1658B)


      1 import { csMarked } from './csMarked.js';
      2 import { csMark } from './csMark.js';
      3 import { csDfs } from './csDfs.js';
      4 /**
      5  * The csReach function computes X = Reach(B), where B is the nonzero pattern of the n-by-1
      6  * sparse column of vector b. The function returns the set of nodes reachable from any node in B. The
      7  * nonzero pattern xi of the solution x to the sparse linear system Lx=b is given by X=Reach(B).
      8  *
      9  * @param {Matrix}  g               The G matrix
     10  * @param {Matrix}  b               The B matrix
     11  * @param {Number}  k               The kth column in B
     12  * @param {Array}   xi              The nonzero pattern xi[top] .. xi[n - 1], an array of size = 2 * n
     13  *                                  The first n entries is the nonzero pattern, the last n entries is the stack
     14  * @param {Array}   pinv            The inverse row permutation vector
     15  *
     16  * @return {Number}                 The index for the nonzero pattern
     17  *
     18  * Reference: http://faculty.cse.tamu.edu/davis/publications.html
     19  */
     20 
     21 export function csReach(g, b, k, xi, pinv) {
     22   // g arrays
     23   var gptr = g._ptr;
     24   var gsize = g._size; // b arrays
     25 
     26   var bindex = b._index;
     27   var bptr = b._ptr; // columns
     28 
     29   var n = gsize[1]; // vars
     30 
     31   var p, p0, p1; // initialize top
     32 
     33   var top = n; // loop column indeces in B
     34 
     35   for (p0 = bptr[k], p1 = bptr[k + 1], p = p0; p < p1; p++) {
     36     // node i
     37     var i = bindex[p]; // check node i is marked
     38 
     39     if (!csMarked(gptr, i)) {
     40       // start a dfs at unmarked node i
     41       top = csDfs(i, g, top, xi, pinv);
     42     }
     43   } // loop columns from top -> n - 1
     44 
     45 
     46   for (p = top; p < n; p++) {
     47     // restore G
     48     csMark(gptr, xi[p]);
     49   }
     50 
     51   return top;
     52 }