simple-squiggle

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

csSymperm.js (2599B)

      1 import { csCumsum } from './csCumsum.js';
      2 import { factory } from '../../../utils/factory.js';
      3 var name = 'csSymperm';
      4 var dependencies = ['conj', 'SparseMatrix'];
      5 export var createCsSymperm = /* #__PURE__ */factory(name, dependencies, _ref => {
      6   var {
      7     conj,
      8     SparseMatrix
      9   } = _ref;
     10 
     11   /**
     12    * Computes the symmetric permutation of matrix A accessing only
     13    * the upper triangular part of A.
     14    *
     15    * C = P * A * P'
     16    *
     17    * @param {Matrix}  a               The A matrix
     18    * @param {Array}   pinv            The inverse of permutation vector
     19    * @param {boolean} values          Process matrix values (true)
     20    *
     21    * @return {Matrix}                 The C matrix, C = P * A * P'
     22    *
     23    * Reference: http://faculty.cse.tamu.edu/davis/publications.html
     24    */
     25   return function csSymperm(a, pinv, values) {
     26     // A matrix arrays
     27     var avalues = a._values;
     28     var aindex = a._index;
     29     var aptr = a._ptr;
     30     var asize = a._size; // columns
     31 
     32     var n = asize[1]; // C matrix arrays
     33 
     34     var cvalues = values && avalues ? [] : null;
     35     var cindex = []; // (nz)
     36 
     37     var cptr = []; // (n + 1)
     38     // variables
     39 
     40     var i, i2, j, j2, p, p0, p1; // create workspace vector
     41 
     42     var w = []; // (n)
     43     // count entries in each column of C
     44 
     45     for (j = 0; j < n; j++) {
     46       // column j of A is column j2 of C
     47       j2 = pinv ? pinv[j] : j; // loop values in column j
     48 
     49       for (p0 = aptr[j], p1 = aptr[j + 1], p = p0; p < p1; p++) {
     50         // row
     51         i = aindex[p]; // skip lower triangular part of A
     52 
     53         if (i > j) {
     54           continue;
     55         } // row i of A is row i2 of C
     56 
     57 
     58         i2 = pinv ? pinv[i] : i; // column count of C
     59 
     60         w[Math.max(i2, j2)]++;
     61       }
     62     } // compute column pointers of C
     63 
     64 
     65     csCumsum(cptr, w, n); // loop columns
     66 
     67     for (j = 0; j < n; j++) {
     68       // column j of A is column j2 of C
     69       j2 = pinv ? pinv[j] : j; // loop values in column j
     70 
     71       for (p0 = aptr[j], p1 = aptr[j + 1], p = p0; p < p1; p++) {
     72         // row
     73         i = aindex[p]; // skip lower triangular part of A
     74 
     75         if (i > j) {
     76           continue;
     77         } // row i of A is row i2 of C
     78 
     79 
     80         i2 = pinv ? pinv[i] : i; // C index for column j2
     81 
     82         var q = w[Math.max(i2, j2)]++; // update C index for entry q
     83 
     84         cindex[q] = Math.min(i2, j2); // check we need to process values
     85 
     86         if (cvalues) {
     87           cvalues[q] = i2 <= j2 ? avalues[p] : conj(avalues[p]);
     88         }
     89       }
     90     } // return C matrix
     91 
     92 
     93     return new SparseMatrix({
     94       values: cvalues,
     95       index: cindex,
     96       ptr: cptr,
     97       size: [n, n]
     98     });
     99   };
    100 });