simple-squiggle

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

usolve.js (4799B)


      1 "use strict";
      2 
      3 Object.defineProperty(exports, "__esModule", {
      4   value: true
      5 });
      6 exports.createUsolve = void 0;
      7 
      8 var _factory = require("../../../utils/factory.js");
      9 
     10 var _solveValidation = require("./utils/solveValidation.js");
     11 
     12 var name = 'usolve';
     13 var dependencies = ['typed', 'matrix', 'divideScalar', 'multiplyScalar', 'subtract', 'equalScalar', 'DenseMatrix'];
     14 var createUsolve = /* #__PURE__ */(0, _factory.factory)(name, dependencies, function (_ref) {
     15   var typed = _ref.typed,
     16       matrix = _ref.matrix,
     17       divideScalar = _ref.divideScalar,
     18       multiplyScalar = _ref.multiplyScalar,
     19       subtract = _ref.subtract,
     20       equalScalar = _ref.equalScalar,
     21       DenseMatrix = _ref.DenseMatrix;
     22   var solveValidation = (0, _solveValidation.createSolveValidation)({
     23     DenseMatrix: DenseMatrix
     24   });
     25   /**
     26    * Finds one solution of a linear equation system by backward substitution. Matrix must be an upper triangular matrix. Throws an error if there's no solution.
     27    *
     28    * `U * x = b`
     29    *
     30    * Syntax:
     31    *
     32    *    math.usolve(U, b)
     33    *
     34    * Examples:
     35    *
     36    *    const a = [[-2, 3], [2, 1]]
     37    *    const b = [11, 9]
     38    *    const x = usolve(a, b)  // [[8], [9]]
     39    *
     40    * See also:
     41    *
     42    *    usolveAll, lup, slu, usolve, lusolve
     43    *
     44    * @param {Matrix, Array} U       A N x N matrix or array (U)
     45    * @param {Matrix, Array} b       A column vector with the b values
     46    *
     47    * @return {DenseMatrix | Array}  A column vector with the linear system solution (x)
     48    */
     49 
     50   return typed(name, {
     51     'SparseMatrix, Array | Matrix': function SparseMatrixArrayMatrix(m, b) {
     52       return _sparseBackwardSubstitution(m, b);
     53     },
     54     'DenseMatrix, Array | Matrix': function DenseMatrixArrayMatrix(m, b) {
     55       return _denseBackwardSubstitution(m, b);
     56     },
     57     'Array, Array | Matrix': function ArrayArrayMatrix(a, b) {
     58       var m = matrix(a);
     59 
     60       var r = _denseBackwardSubstitution(m, b);
     61 
     62       return r.valueOf();
     63     }
     64   });
     65 
     66   function _denseBackwardSubstitution(m, b) {
     67     // make b into a column vector
     68     b = solveValidation(m, b, true);
     69     var bdata = b._data;
     70     var rows = m._size[0];
     71     var columns = m._size[1]; // result
     72 
     73     var x = [];
     74     var mdata = m._data; // loop columns backwards
     75 
     76     for (var j = columns - 1; j >= 0; j--) {
     77       // b[j]
     78       var bj = bdata[j][0] || 0; // x[j]
     79 
     80       var xj = void 0;
     81 
     82       if (!equalScalar(bj, 0)) {
     83         // value at [j, j]
     84         var vjj = mdata[j][j];
     85 
     86         if (equalScalar(vjj, 0)) {
     87           // system cannot be solved
     88           throw new Error('Linear system cannot be solved since matrix is singular');
     89         }
     90 
     91         xj = divideScalar(bj, vjj); // loop rows
     92 
     93         for (var i = j - 1; i >= 0; i--) {
     94           // update copy of b
     95           bdata[i] = [subtract(bdata[i][0] || 0, multiplyScalar(xj, mdata[i][j]))];
     96         }
     97       } else {
     98         // zero value at j
     99         xj = 0;
    100       } // update x
    101 
    102 
    103       x[j] = [xj];
    104     }
    105 
    106     return new DenseMatrix({
    107       data: x,
    108       size: [rows, 1]
    109     });
    110   }
    111 
    112   function _sparseBackwardSubstitution(m, b) {
    113     // make b into a column vector
    114     b = solveValidation(m, b, true);
    115     var bdata = b._data;
    116     var rows = m._size[0];
    117     var columns = m._size[1];
    118     var values = m._values;
    119     var index = m._index;
    120     var ptr = m._ptr; // result
    121 
    122     var x = []; // loop columns backwards
    123 
    124     for (var j = columns - 1; j >= 0; j--) {
    125       var bj = bdata[j][0] || 0;
    126 
    127       if (!equalScalar(bj, 0)) {
    128         // non-degenerate row, find solution
    129         var vjj = 0; // upper triangular matrix values & index (column j)
    130 
    131         var jValues = [];
    132         var jIndices = []; // first & last indeces in column
    133 
    134         var firstIndex = ptr[j];
    135         var lastIndex = ptr[j + 1]; // values in column, find value at [j, j], loop backwards
    136 
    137         for (var k = lastIndex - 1; k >= firstIndex; k--) {
    138           var i = index[k]; // check row (rows are not sorted!)
    139 
    140           if (i === j) {
    141             vjj = values[k];
    142           } else if (i < j) {
    143             // store upper triangular
    144             jValues.push(values[k]);
    145             jIndices.push(i);
    146           }
    147         } // at this point we must have a value in vjj
    148 
    149 
    150         if (equalScalar(vjj, 0)) {
    151           throw new Error('Linear system cannot be solved since matrix is singular');
    152         }
    153 
    154         var xj = divideScalar(bj, vjj);
    155 
    156         for (var _k = 0, _lastIndex = jIndices.length; _k < _lastIndex; _k++) {
    157           var _i = jIndices[_k];
    158           bdata[_i] = [subtract(bdata[_i][0], multiplyScalar(xj, jValues[_k]))];
    159         }
    160 
    161         x[j] = [xj];
    162       } else {
    163         // degenerate row, we can choose any value
    164         x[j] = [0];
    165       }
    166     }
    167 
    168     return new DenseMatrix({
    169       data: x,
    170       size: [rows, 1]
    171     });
    172   }
    173 });
    174 exports.createUsolve = createUsolve;