simple-squiggle

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

slu.js (3554B)

      1 import { isInteger } from '../../../utils/number.js';
      2 import { factory } from '../../../utils/factory.js';
      3 import { createCsSqr } from '../sparse/csSqr.js';
      4 import { createCsLu } from '../sparse/csLu.js';
      5 var name = 'slu';
      6 var dependencies = ['typed', 'abs', 'add', 'multiply', 'transpose', 'divideScalar', 'subtract', 'larger', 'largerEq', 'SparseMatrix'];
      7 export var createSlu = /* #__PURE__ */factory(name, dependencies, _ref => {
      8   var {
      9     typed,
     10     abs,
     11     add,
     12     multiply,
     13     transpose,
     14     divideScalar,
     15     subtract,
     16     larger,
     17     largerEq,
     18     SparseMatrix
     19   } = _ref;
     20   var csSqr = createCsSqr({
     21     add,
     22     multiply,
     23     transpose
     24   });
     25   var csLu = createCsLu({
     26     abs,
     27     divideScalar,
     28     multiply,
     29     subtract,
     30     larger,
     31     largerEq,
     32     SparseMatrix
     33   });
     34   /**
     35    * Calculate the Sparse Matrix LU decomposition with full pivoting. Sparse Matrix `A` is decomposed in two matrices (`L`, `U`) and two permutation vectors (`pinv`, `q`) where
     36    *
     37    * `P * A * Q = L * U`
     38    *
     39    * Syntax:
     40    *
     41    *    math.slu(A, order, threshold)
     42    *
     43    * Examples:
     44    *
     45    *    const A = math.sparse([[4,3], [6, 3]])
     46    *    math.slu(A, 1, 0.001)
     47    *    // returns:
     48    *    // {
     49    *    //   L: [[1, 0], [1.5, 1]]
     50    *    //   U: [[4, 3], [0, -1.5]]
     51    *    //   p: [0, 1]
     52    *    //   q: [0, 1]
     53    *    // }
     54    *
     55    * See also:
     56    *
     57    *    lup, lsolve, usolve, lusolve
     58    *
     59    * @param {SparseMatrix} A              A two dimensional sparse matrix for which to get the LU decomposition.
     60    * @param {Number}       order          The Symbolic Ordering and Analysis order:
     61    *                                       0 - Natural ordering, no permutation vector q is returned
     62    *                                       1 - Matrix must be square, symbolic ordering and analisis is performed on M = A + A'
     63    *                                       2 - Symbolic ordering and analisis is performed on M = A' * A. Dense columns from A' are dropped, A recreated from A'.
     64    *                                           This is appropriatefor LU factorization of unsymmetric matrices.
     65    *                                       3 - Symbolic ordering and analisis is performed on M = A' * A. This is best used for LU factorization is matrix M has no dense rows.
     66    *                                           A dense row is a row with more than 10*sqr(columns) entries.
     67    * @param {Number}       threshold       Partial pivoting threshold (1 for partial pivoting)
     68    *
     69    * @return {Object} The lower triangular matrix, the upper triangular matrix and the permutation vectors.
     70    */
     71 
     72   return typed(name, {
     73     'SparseMatrix, number, number': function SparseMatrixNumberNumber(a, order, threshold) {
     74       // verify order
     75       if (!isInteger(order) || order < 0 || order > 3) {
     76         throw new Error('Symbolic Ordering and Analysis order must be an integer number in the interval [0, 3]');
     77       } // verify threshold
     78 
     79 
     80       if (threshold < 0 || threshold > 1) {
     81         throw new Error('Partial pivoting threshold must be a number from 0 to 1');
     82       } // perform symbolic ordering and analysis
     83 
     84 
     85       var s = csSqr(order, a, false); // perform lu decomposition
     86 
     87       var f = csLu(a, s, threshold); // return decomposition
     88 
     89       return {
     90         L: f.L,
     91         U: f.U,
     92         p: f.pinv,
     93         q: s.q,
     94         toString: function toString() {
     95           return 'L: ' + this.L.toString() + '\nU: ' + this.U.toString() + '\np: ' + this.p.toString() + (this.q ? '\nq: ' + this.q.toString() : '') + '\n';
     96         }
     97       };
     98     }
     99   });
    100 });