simple-squiggle

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

csSymperm.js (2806B)

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