simple-squiggle

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

csChol.js (3839B)

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