simple-squiggle

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

csChol.js (4198B)

      1 "use strict";
      2 
      3 Object.defineProperty(exports, "__esModule", {
      4   value: true
      5 });
      6 exports.createCsChol = void 0;
      7 
      8 var _factory = require("../../../utils/factory.js");
      9 
     10 var _csEreach = require("./csEreach.js");
     11 
     12 var _csSymperm = require("./csSymperm.js");
     13 
     14 var name = 'csChol';
     15 var dependencies = ['divideScalar', 'sqrt', 'subtract', 'multiply', 'im', 're', 'conj', 'equal', 'smallerEq', 'SparseMatrix'];
     16 var createCsChol = /* #__PURE__ */(0, _factory.factory)(name, dependencies, function (_ref) {
     17   var divideScalar = _ref.divideScalar,
     18       sqrt = _ref.sqrt,
     19       subtract = _ref.subtract,
     20       multiply = _ref.multiply,
     21       im = _ref.im,
     22       re = _ref.re,
     23       conj = _ref.conj,
     24       equal = _ref.equal,
     25       smallerEq = _ref.smallerEq,
     26       SparseMatrix = _ref.SparseMatrix;
     27   var csSymperm = (0, _csSymperm.createCsSymperm)({
     28     conj: conj,
     29     SparseMatrix: SparseMatrix
     30   });
     31   /**
     32    * Computes the Cholesky factorization of matrix A. It computes L and P so
     33    * L * L' = P * A * P'
     34    *
     35    * @param {Matrix}  m               The A Matrix to factorize, only upper triangular part used
     36    * @param {Object}  s               The symbolic analysis from cs_schol()
     37    *
     38    * @return {Number}                 The numeric Cholesky factorization of A or null
     39    *
     40    * Reference: http://faculty.cse.tamu.edu/davis/publications.html
     41    */
     42 
     43   return function csChol(m, s) {
     44     // validate input
     45     if (!m) {
     46       return null;
     47     } // m arrays
     48 
     49 
     50     var size = m._size; // columns
     51 
     52     var n = size[1]; // symbolic analysis result
     53 
     54     var parent = s.parent;
     55     var cp = s.cp;
     56     var pinv = s.pinv; // L arrays
     57 
     58     var lvalues = [];
     59     var lindex = [];
     60     var lptr = []; // L
     61 
     62     var L = new SparseMatrix({
     63       values: lvalues,
     64       index: lindex,
     65       ptr: lptr,
     66       size: [n, n]
     67     }); // vars
     68 
     69     var c = []; // (2 * n)
     70 
     71     var x = []; // (n)
     72     // compute C = P * A * P'
     73 
     74     var cm = pinv ? csSymperm(m, pinv, 1) : m; // C matrix arrays
     75 
     76     var cvalues = cm._values;
     77     var cindex = cm._index;
     78     var cptr = cm._ptr; // vars
     79 
     80     var k, p; // initialize variables
     81 
     82     for (k = 0; k < n; k++) {
     83       lptr[k] = c[k] = cp[k];
     84     } // compute L(k,:) for L*L' = C
     85 
     86 
     87     for (k = 0; k < n; k++) {
     88       // nonzero pattern of L(k,:)
     89       var top = (0, _csEreach.csEreach)(cm, k, parent, c); // x (0:k) is now zero
     90 
     91       x[k] = 0; // x = full(triu(C(:,k)))
     92 
     93       for (p = cptr[k]; p < cptr[k + 1]; p++) {
     94         if (cindex[p] <= k) {
     95           x[cindex[p]] = cvalues[p];
     96         }
     97       } // d = C(k,k)
     98 
     99 
    100       var d = x[k]; // clear x for k+1st iteration
    101 
    102       x[k] = 0; // solve L(0:k-1,0:k-1) * x = C(:,k)
    103 
    104       for (; top < n; top++) {
    105         // s[top..n-1] is pattern of L(k,:)
    106         var i = s[top]; // L(k,i) = x (i) / L(i,i)
    107 
    108         var lki = divideScalar(x[i], lvalues[lptr[i]]); // clear x for k+1st iteration
    109 
    110         x[i] = 0;
    111 
    112         for (p = lptr[i] + 1; p < c[i]; p++) {
    113           // row
    114           var r = lindex[p]; // update x[r]
    115 
    116           x[r] = subtract(x[r], multiply(lvalues[p], lki));
    117         } // d = d - L(k,i)*L(k,i)
    118 
    119 
    120         d = subtract(d, multiply(lki, conj(lki)));
    121         p = c[i]++; // store L(k,i) in column i
    122 
    123         lindex[p] = k;
    124         lvalues[p] = conj(lki);
    125       } // compute L(k,k)
    126 
    127 
    128       if (smallerEq(re(d), 0) || !equal(im(d), 0)) {
    129         // not pos def
    130         return null;
    131       }
    132 
    133       p = c[k]++; //  store L(k,k) = sqrt(d) in column k
    134 
    135       lindex[p] = k;
    136       lvalues[p] = sqrt(d);
    137     } // finalize L
    138 
    139 
    140     lptr[n] = cp[n]; // P matrix
    141 
    142     var P; // check we need to calculate P
    143 
    144     if (pinv) {
    145       // P arrays
    146       var pvalues = [];
    147       var pindex = [];
    148       var pptr = []; // create P matrix
    149 
    150       for (p = 0; p < n; p++) {
    151         // initialize ptr (one value per column)
    152         pptr[p] = p; // index (apply permutation vector)
    153 
    154         pindex.push(pinv[p]); // value 1
    155 
    156         pvalues.push(1);
    157       } // update ptr
    158 
    159 
    160       pptr[n] = n; // P
    161 
    162       P = new SparseMatrix({
    163         values: pvalues,
    164         index: pindex,
    165         ptr: pptr,
    166         size: [n, n]
    167       });
    168     } // return L & P
    169 
    170 
    171     return {
    172       L: L,
    173       P: P
    174     };
    175   };
    176 });
    177 exports.createCsChol = createCsChol;