1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
|
// Copyright 2011 The Closure Library Authors. All Rights Reserved.
//
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS-IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
/**
* @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) {
// Make a local copy of the main diagonal and the right vector.
mainDiag = mainDiag.slice();
vecRight = vecRight.slice();
// The dimension of the matrix.
var nDim = mainDiag.length;
// Construct a modified linear system of equations with the same solution
// as the input one.
for (var i = 1; i < nDim; ++i) {
var 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.
var 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;
};
|
