chromium/third_party/google-closure-library/closure/goog/math/tdma.js

/**
 * @license
 * Copyright The Closure Library Authors.
 * SPDX-License-Identifier: Apache-2.0
 */

/**
 * @fileoverview The Tridiagonal matrix algorithm solver solves a special
 * version of a sparse linear system Ax = b where A is tridiagonal.
 *
 * See http://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm
 */

goog.provide('goog.math.tdma');


/**
 * Solves a linear system where the matrix is square tri-diagonal. That is,
 * given a system of equations:
 *
 * A * result = vecRight,
 *
 * this class computes result = inv(A) * vecRight, where A has the special form
 * of a tri-diagonal matrix:
 *
 *    |dia(0) sup(0)   0    0     ...   0|
 *    |sub(0) dia(1) sup(1) 0     ...   0|
 * A =|                ...               |
 *    |0 ... 0 sub(n-2) dia(n-1) sup(n-1)|
 *    |0 ... 0    0     sub(n-1)   dia(n)|
 *
 * @param {!Array<number>} subDiag The sub diagonal of the matrix.
 * @param {!Array<number>} mainDiag The main diagonal of the matrix.
 * @param {!Array<number>} supDiag The super diagonal of the matrix.
 * @param {!Array<number>} vecRight The right vector of the system
 *     of equations.
 * @param {Array<number>=} opt_result The optional array to store the result.
 * @return {!Array<number>} The vector that is the solution to the system.
 */
goog.math.tdma.solve = function(
    subDiag, mainDiag, supDiag, vecRight, opt_result) {
  'use strict';
  // Make a local copy of the main diagonal and the right vector.
  mainDiag = mainDiag.slice();
  vecRight = vecRight.slice();

  // The dimension of the matrix.
  const nDim = mainDiag.length;

  // Construct a modified linear system of equations with the same solution
  // as the input one.
  let i;
  for (i = 1; i < nDim; ++i) {
    const m = subDiag[i - 1] / mainDiag[i - 1];
    mainDiag[i] = mainDiag[i] - m * supDiag[i - 1];
    vecRight[i] = vecRight[i] - m * vecRight[i - 1];
  }

  // Solve the new system of equations by simple back-substitution.
  const result = opt_result || new Array(vecRight.length);
  result[nDim - 1] = vecRight[nDim - 1] / mainDiag[nDim - 1];
  for (i = nDim - 2; i >= 0; --i) {
    result[i] = (vecRight[i] - supDiag[i] * result[i + 1]) / mainDiag[i];
  }
  return result;
};