simple-squiggle

lsolveAll.js (6125B)

---

 "use strict";
 
 var _interopRequireDefault = require("@babel/runtime/helpers/interopRequireDefault");
 
 Object.defineProperty(exports, "__esModule", {
   value: true
 });
 exports.createLsolveAll = void 0;
 
 var _toConsumableArray2 = _interopRequireDefault(require("@babel/runtime/helpers/toConsumableArray"));
 
 var _factory = require("../../../utils/factory.js");
 
 var _solveValidation = require("./utils/solveValidation.js");
 
 var name = 'lsolveAll';
 var dependencies = ['typed', 'matrix', 'divideScalar', 'multiplyScalar', 'subtract', 'equalScalar', 'DenseMatrix'];
 var createLsolveAll = /* #__PURE__ */(0, _factory.factory)(name, dependencies, function (_ref) {
   var typed = _ref.typed,
       matrix = _ref.matrix,
       divideScalar = _ref.divideScalar,
       multiplyScalar = _ref.multiplyScalar,
       subtract = _ref.subtract,
       equalScalar = _ref.equalScalar,
       DenseMatrix = _ref.DenseMatrix;
   var solveValidation = (0, _solveValidation.createSolveValidation)({
     DenseMatrix: DenseMatrix
   });
   /**
    * Finds all solutions of a linear equation system by forwards substitution. Matrix must be a lower triangular matrix.
    *
    * `L * x = b`
    *
    * Syntax:
    *
    *    math.lsolveAll(L, b)
    *
    * Examples:
    *
    *    const a = [[-2, 3], [2, 1]]
    *    const b = [11, 9]
    *    const x = lsolveAll(a, b)  // [ [[-5.5], [20]] ]
    *
    * See also:
    *
    *    lsolve, lup, slu, usolve, lusolve
    *
    * @param {Matrix, Array} L       A N x N matrix or array (L)
    * @param {Matrix, Array} b       A column vector with the b values
    *
    * @return {DenseMatrix[] | Array[]}  An array of affine-independent column vectors (x) that solve the linear system
    */
 
   return typed(name, {
     'SparseMatrix, Array | Matrix': function SparseMatrixArrayMatrix(m, b) {
       return _sparseForwardSubstitution(m, b);
     },
     'DenseMatrix, Array | Matrix': function DenseMatrixArrayMatrix(m, b) {
       return _denseForwardSubstitution(m, b);
     },
     'Array, Array | Matrix': function ArrayArrayMatrix(a, b) {
       var m = matrix(a);
 
       var R = _denseForwardSubstitution(m, b);
 
       return R.map(function (r) {
         return r.valueOf();
       });
     }
   });
 
   function _denseForwardSubstitution(m, b_) {
     // the algorithm is derived from
     // https://www.overleaf.com/read/csvgqdxggyjv
     // array of right-hand sides
     var B = [solveValidation(m, b_, true)._data.map(function (e) {
       return e[0];
     })];
     var M = m._data;
     var rows = m._size[0];
     var columns = m._size[1]; // loop columns
 
     for (var i = 0; i < columns; i++) {
       var L = B.length; // loop right-hand sides
 
       for (var k = 0; k < L; k++) {
         var b = B[k];
 
         if (!equalScalar(M[i][i], 0)) {
           // non-singular row
           b[i] = divideScalar(b[i], M[i][i]);
 
           for (var j = i + 1; j < columns; j++) {
             // b[j] -= b[i] * M[j,i]
             b[j] = subtract(b[j], multiplyScalar(b[i], M[j][i]));
           }
         } else if (!equalScalar(b[i], 0)) {
           // singular row, nonzero RHS
           if (k === 0) {
             // There is no valid solution
             return [];
           } else {
             // This RHS is invalid but other solutions may still exist
             B.splice(k, 1);
             k -= 1;
             L -= 1;
           }
         } else if (k === 0) {
           // singular row, RHS is zero
           var bNew = (0, _toConsumableArray2.default)(b);
           bNew[i] = 1;
 
           for (var _j = i + 1; _j < columns; _j++) {
             bNew[_j] = subtract(bNew[_j], M[_j][i]);
           }
 
           B.push(bNew);
         }
       }
     }
 
     return B.map(function (x) {
       return new DenseMatrix({
         data: x.map(function (e) {
           return [e];
         }),
         size: [rows, 1]
       });
     });
   }
 
   function _sparseForwardSubstitution(m, b_) {
     // array of right-hand sides
     var B = [solveValidation(m, b_, true)._data.map(function (e) {
       return e[0];
     })];
     var rows = m._size[0];
     var columns = m._size[1];
     var values = m._values;
     var index = m._index;
     var ptr = m._ptr; // loop columns
 
     for (var i = 0; i < columns; i++) {
       var L = B.length; // loop right-hand sides
 
       for (var k = 0; k < L; k++) {
         var b = B[k]; // values & indices (column i)
 
         var iValues = [];
         var iIndices = []; // first & last indeces in column
 
         var firstIndex = ptr[i];
         var lastIndex = ptr[i + 1]; // find the value at [i, i]
 
         var Mii = 0;
 
         for (var j = firstIndex; j < lastIndex; j++) {
           var J = index[j]; // check row
 
           if (J === i) {
             Mii = values[j];
           } else if (J > i) {
             // store lower triangular
             iValues.push(values[j]);
             iIndices.push(J);
           }
         }
 
         if (!equalScalar(Mii, 0)) {
           // non-singular row
           b[i] = divideScalar(b[i], Mii);
 
           for (var _j2 = 0, _lastIndex = iIndices.length; _j2 < _lastIndex; _j2++) {
             var _J = iIndices[_j2];
             b[_J] = subtract(b[_J], multiplyScalar(b[i], iValues[_j2]));
           }
         } else if (!equalScalar(b[i], 0)) {
           // singular row, nonzero RHS
           if (k === 0) {
             // There is no valid solution
             return [];
           } else {
             // This RHS is invalid but other solutions may still exist
             B.splice(k, 1);
             k -= 1;
             L -= 1;
           }
         } else if (k === 0) {
           // singular row, RHS is zero
           var bNew = (0, _toConsumableArray2.default)(b);
           bNew[i] = 1;
 
           for (var _j3 = 0, _lastIndex2 = iIndices.length; _j3 < _lastIndex2; _j3++) {
             var _J2 = iIndices[_j3];
             bNew[_J2] = subtract(bNew[_J2], iValues[_j3]);
           }
 
           B.push(bNew);
         }
       }
     }
 
     return B.map(function (x) {
       return new DenseMatrix({
         data: x.map(function (e) {
           return [e];
         }),
         size: [rows, 1]
       });
     });
   }
 });
 exports.createLsolveAll = createLsolveAll;