simple-squiggle

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

csEreach.js (1789B)


      1 import { csMark } from './csMark.js';
      2 import { csMarked } from './csMarked.js';
      3 /**
      4  * Find nonzero pattern of Cholesky L(k,1:k-1) using etree and triu(A(:,k))
      5  *
      6  * @param {Matrix}  a               The A matrix
      7  * @param {Number}  k               The kth column in A
      8  * @param {Array}   parent          The parent vector from the symbolic analysis result
      9  * @param {Array}   w               The nonzero pattern xi[top] .. xi[n - 1], an array of size = 2 * n
     10  *                                  The first n entries is the nonzero pattern, the last n entries is the stack
     11  *
     12  * @return {Number}                 The index for the nonzero pattern
     13  *
     14  * Reference: http://faculty.cse.tamu.edu/davis/publications.html
     15  */
     16 
     17 export function csEreach(a, k, parent, w) {
     18   // a arrays
     19   var aindex = a._index;
     20   var aptr = a._ptr;
     21   var asize = a._size; // columns
     22 
     23   var n = asize[1]; // initialize top
     24 
     25   var top = n; // vars
     26 
     27   var p, p0, p1, len; // mark node k as visited
     28 
     29   csMark(w, k); // loop values & index for column k
     30 
     31   for (p0 = aptr[k], p1 = aptr[k + 1], p = p0; p < p1; p++) {
     32     // A(i,k) is nonzero
     33     var i = aindex[p]; // only use upper triangular part of A
     34 
     35     if (i > k) {
     36       continue;
     37     } // traverse up etree
     38 
     39 
     40     for (len = 0; !csMarked(w, i); i = parent[i]) {
     41       // L(k,i) is nonzero, last n entries in w
     42       w[n + len++] = i; // mark i as visited
     43 
     44       csMark(w, i);
     45     }
     46 
     47     while (len > 0) {
     48       // decrement top & len
     49       --top;
     50       --len; // push path onto stack, last n entries in w
     51 
     52       w[n + top] = w[n + len];
     53     }
     54   } // unmark all nodes
     55 
     56 
     57   for (p = top; p < n; p++) {
     58     // use stack value, last n entries in w
     59     csMark(w, w[n + p]);
     60   } // unmark node k
     61 
     62 
     63   csMark(w, k); // s[top..n-1] contains pattern of L(k,:)
     64 
     65   return top;
     66 }