simple-squiggle

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

lsolve.js (4725B)


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