simple-squiggle

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

lusolve.js (4049B)

      1 "use strict";
      2 
      3 Object.defineProperty(exports, "__esModule", {
      4   value: true
      5 });
      6 exports.createLusolve = void 0;
      7 
      8 var _is = require("../../../utils/is.js");
      9 
     10 var _factory = require("../../../utils/factory.js");
     11 
     12 var _solveValidation = require("./utils/solveValidation.js");
     13 
     14 var _csIpvec = require("../sparse/csIpvec.js");
     15 
     16 var name = 'lusolve';
     17 var dependencies = ['typed', 'matrix', 'lup', 'slu', 'usolve', 'lsolve', 'DenseMatrix'];
     18 var createLusolve = /* #__PURE__ */(0, _factory.factory)(name, dependencies, function (_ref) {
     19   var typed = _ref.typed,
     20       matrix = _ref.matrix,
     21       lup = _ref.lup,
     22       slu = _ref.slu,
     23       usolve = _ref.usolve,
     24       lsolve = _ref.lsolve,
     25       DenseMatrix = _ref.DenseMatrix;
     26   var solveValidation = (0, _solveValidation.createSolveValidation)({
     27     DenseMatrix: DenseMatrix
     28   });
     29   /**
     30    * Solves the linear system `A * x = b` where `A` is an [n x n] matrix and `b` is a [n] column vector.
     31    *
     32    * Syntax:
     33    *
     34    *    math.lusolve(A, b)     // returns column vector with the solution to the linear system A * x = b
     35    *    math.lusolve(lup, b)   // returns column vector with the solution to the linear system A * x = b, lup = math.lup(A)
     36    *
     37    * Examples:
     38    *
     39    *    const m = [[1, 0, 0, 0], [0, 2, 0, 0], [0, 0, 3, 0], [0, 0, 0, 4]]
     40    *
     41    *    const x = math.lusolve(m, [-1, -1, -1, -1])        // x = [[-1], [-0.5], [-1/3], [-0.25]]
     42    *
     43    *    const f = math.lup(m)
     44    *    const x1 = math.lusolve(f, [-1, -1, -1, -1])       // x1 = [[-1], [-0.5], [-1/3], [-0.25]]
     45    *    const x2 = math.lusolve(f, [1, 2, 1, -1])          // x2 = [[1], [1], [1/3], [-0.25]]
     46    *
     47    *    const a = [[-2, 3], [2, 1]]
     48    *    const b = [11, 9]
     49    *    const x = math.lusolve(a, b)  // [[2], [5]]
     50    *
     51    * See also:
     52    *
     53    *    lup, slu, lsolve, usolve
     54    *
     55    * @param {Matrix | Array | Object} A      Invertible Matrix or the Matrix LU decomposition
     56    * @param {Matrix | Array} b               Column Vector
     57    * @param {number} [order]                 The Symbolic Ordering and Analysis order, see slu for details. Matrix must be a SparseMatrix
     58    * @param {Number} [threshold]             Partial pivoting threshold (1 for partial pivoting), see slu for details. Matrix must be a SparseMatrix.
     59    *
     60    * @return {DenseMatrix | Array}           Column vector with the solution to the linear system A * x = b
     61    */
     62 
     63   return typed(name, {
     64     'Array, Array | Matrix': function ArrayArrayMatrix(a, b) {
     65       a = matrix(a);
     66       var d = lup(a);
     67 
     68       var x = _lusolve(d.L, d.U, d.p, null, b);
     69 
     70       return x.valueOf();
     71     },
     72     'DenseMatrix, Array | Matrix': function DenseMatrixArrayMatrix(a, b) {
     73       var d = lup(a);
     74       return _lusolve(d.L, d.U, d.p, null, b);
     75     },
     76     'SparseMatrix, Array | Matrix': function SparseMatrixArrayMatrix(a, b) {
     77       var d = lup(a);
     78       return _lusolve(d.L, d.U, d.p, null, b);
     79     },
     80     'SparseMatrix, Array | Matrix, number, number': function SparseMatrixArrayMatrixNumberNumber(a, b, order, threshold) {
     81       var d = slu(a, order, threshold);
     82       return _lusolve(d.L, d.U, d.p, d.q, b);
     83     },
     84     'Object, Array | Matrix': function ObjectArrayMatrix(d, b) {
     85       return _lusolve(d.L, d.U, d.p, d.q, b);
     86     }
     87   });
     88 
     89   function _toMatrix(a) {
     90     if ((0, _is.isMatrix)(a)) {
     91       return a;
     92     }
     93 
     94     if ((0, _is.isArray)(a)) {
     95       return matrix(a);
     96     }
     97 
     98     throw new TypeError('Invalid Matrix LU decomposition');
     99   }
    100 
    101   function _lusolve(l, u, p, q, b) {
    102     // verify decomposition
    103     l = _toMatrix(l);
    104     u = _toMatrix(u); // apply row permutations if needed (b is a DenseMatrix)
    105 
    106     if (p) {
    107       b = solveValidation(l, b, true);
    108       b._data = (0, _csIpvec.csIpvec)(p, b._data);
    109     } // use forward substitution to resolve L * y = b
    110 
    111 
    112     var y = lsolve(l, b); // use backward substitution to resolve U * x = y
    113 
    114     var x = usolve(u, y); // apply column permutations if needed (x is a DenseMatrix)
    115 
    116     if (q) {
    117       x._data = (0, _csIpvec.csIpvec)(q, x._data);
    118     }
    119 
    120     return x;
    121   }
    122 });
    123 exports.createLusolve = createLusolve;