simple-squiggle

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

csEtree.js (1618B)


      1 /**
      2  * Computes the elimination tree of Matrix A (using triu(A)) or the
      3  * elimination tree of A'A without forming A'A.
      4  *
      5  * @param {Matrix}  a               The A Matrix
      6  * @param {boolean} ata             A value of true the function computes the etree of A'A
      7  *
      8  * Reference: http://faculty.cse.tamu.edu/davis/publications.html
      9  */
     10 export function csEtree(a, ata) {
     11   // check inputs
     12   if (!a) {
     13     return null;
     14   } // a arrays
     15 
     16 
     17   var aindex = a._index;
     18   var aptr = a._ptr;
     19   var asize = a._size; // rows & columns
     20 
     21   var m = asize[0];
     22   var n = asize[1]; // allocate result
     23 
     24   var parent = []; // (n)
     25   // allocate workspace
     26 
     27   var w = []; // (n + (ata ? m : 0))
     28 
     29   var ancestor = 0; // first n entries in w
     30 
     31   var prev = n; // last m entries (ata = true)
     32 
     33   var i, inext; // check we are calculating A'A
     34 
     35   if (ata) {
     36     // initialize workspace
     37     for (i = 0; i < m; i++) {
     38       w[prev + i] = -1;
     39     }
     40   } // loop columns
     41 
     42 
     43   for (var k = 0; k < n; k++) {
     44     // node k has no parent yet
     45     parent[k] = -1; // nor does k have an ancestor
     46 
     47     w[ancestor + k] = -1; // values in column k
     48 
     49     for (var p0 = aptr[k], p1 = aptr[k + 1], p = p0; p < p1; p++) {
     50       // row
     51       var r = aindex[p]; // node
     52 
     53       i = ata ? w[prev + r] : r; // traverse from i to k
     54 
     55       for (; i !== -1 && i < k; i = inext) {
     56         // inext = ancestor of i
     57         inext = w[ancestor + i]; // path compression
     58 
     59         w[ancestor + i] = k; // check no anc., parent is k
     60 
     61         if (inext === -1) {
     62           parent[i] = k;
     63         }
     64       }
     65 
     66       if (ata) {
     67         w[prev + r] = k;
     68       }
     69     }
     70   }
     71 
     72   return parent;
     73 }