1
// Copyright 2011 The Closure Library Authors. All Rights Reserved.
2
//
3
// Licensed under the Apache License, Version 2.0 (the "License");
4
// you may not use this file except in compliance with the License.
5
// You may obtain a copy of the License at
6
//
7
//      http://www.apache.org/licenses/LICENSE-2.0
8
//
9
// Unless required by applicable law or agreed to in writing, software
10
// distributed under the License is distributed on an "AS-IS" BASIS,
11
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
12
// See the License for the specific language governing permissions and
13
// limitations under the License.
14

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

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

25

26
/**
27
 * Solves a linear system where the matrix is square tri-diagonal. That is,
28
 * given a system of equations:
29
 *
30
 * A * result = vecRight,
31
 *
32
 * this class computes result = inv(A) * vecRight, where A has the special form
33
 * of a tri-diagonal matrix:
34
 *
35
 *    |dia(0) sup(0)   0    0     ...   0|
36
 *    |sub(0) dia(1) sup(1) 0     ...   0|
37
 * A =|                ...               |
38
 *    |0 ... 0 sub(n-2) dia(n-1) sup(n-1)|
39
 *    |0 ... 0    0     sub(n-1)   dia(n)|
40
 *
41
 * @param {!Array<number>} subDiag The sub diagonal of the matrix.
42
 * @param {!Array<number>} mainDiag The main diagonal of the matrix.
43
 * @param {!Array<number>} supDiag The super diagonal of the matrix.
44
 * @param {!Array<number>} vecRight The right vector of the system
45
 *     of equations.
46
 * @param {Array<number>=} opt_result The optional array to store the result.
47
 * @return {!Array<number>} The vector that is the solution to the system.
48
 */
49
goog.math.tdma.solve = function(
50
    subDiag, mainDiag, supDiag, vecRight, opt_result) {
51
  // Make a local copy of the main diagonal and the right vector.
52
  mainDiag = mainDiag.slice();
53
  vecRight = vecRight.slice();
54

55
  // The dimension of the matrix.
56
  var nDim = mainDiag.length;
57

58
  // Construct a modified linear system of equations with the same solution
59
  // as the input one.
60
  for (var i = 1; i < nDim; ++i) {
61
    var m = subDiag[i - 1] / mainDiag[i - 1];
62
    mainDiag[i] = mainDiag[i] - m * supDiag[i - 1];
63
    vecRight[i] = vecRight[i] - m * vecRight[i - 1];
64
  }
65

66
  // Solve the new system of equations by simple back-substitution.
67
  var result = opt_result || new Array(vecRight.length);
68
  result[nDim - 1] = vecRight[nDim - 1] / mainDiag[nDim - 1];
69
  for (i = nDim - 2; i >= 0; --i) {
70
    result[i] = (vecRight[i] - supDiag[i] * result[i + 1]) / mainDiag[i];
71
  }
72
  return result;
73
};
74

75