三重対角行列アルゴリズム

Improved reduction for specific matrices

数値線形代数において、三重対角行列アルゴリズム（トーマスアルゴリズムとも呼ばれる）は、ガウス消去法の簡略化された形式で、三重対角連立方程式を解くのに使用できます。n個の未知数を持つ三重対角連立方程式は次のように表すことができます。

${\displaystyle a_{i}x_{i-1}+b_{i}x_{i}+c_{i}x_{i+1}=d_{i},}$

ここで、および ${\displaystyle a_{1}=0}$ ${\displaystyle c_{n}=0}$

${\displaystyle {\begin{bmatrix}b_{1}&c_{1}&&&0\\a_{2}&b_{2}&c_{2}&&\\&a_{3}&b_{3}&\ddots &\\&&\ddots &\ddots &c_{n-1}\\0&&&a_{n}&b_{n}\end{bmatrix}}{\begin{bmatrix}x_{1}\\x_{2}\\x_{3}\\\vdots \\x_{n}\end{bmatrix}}={\begin{bmatrix}d_{1}\\d_{2}\\d_{3}\\\vdots \\d_{n}\end{bmatrix}}.}$

このようなシステムでは、ガウス消去法ではなく、演算で解を得ることができます。最初のスイープで を消去し、次に（省略された）後退代入によって解を生成します。このような行列の例は、 1次元ポアソン方程式の離散化や自然3次スプライン補間からよく生じます。 ${\displaystyle O(n)}$ ${\displaystyle O(n^{3})}$ ${\displaystyle a_{i}}$

トーマスのアルゴリズムは一般に安定ではありませんが、行列が対角優勢（行または列のいずれか）である場合や対称正定値である場合など、いくつかの特殊なケースでは安定します。^{[1] [2]}トーマスのアルゴリズムの安定性のより正確な特徴付けについては、ハイアムの定理9.12を参照してください。^[3]一般的なケースで安定性が必要な場合は、代わりに部分ピボット付きガウス消去法（GEPP）が推奨されます。^[2]

方法

前進スイープは、次のように新しい係数を計算することから成り、新しい係数はプライムで表します。

${\displaystyle c'_{i}={\begin{cases}{\cfrac {c_{i}}{b_{i}}},&i=1,\\{\cfrac {c_{i}}{b_{i}-a_{i}c'_{i-1}}},&i=2,3,\dots ,n-1\end{cases}}}$

と

${\displaystyle d'_{i}={\begin{cases}{\cfrac {d_{i}}{b_{i}}},&i=1,\\{\cfrac {d_{i}-a_{i}d'_{i-1}}{b_{i}-a_{i}c'_{i-1}}},&i=2,3,\dots ,n.\end{cases}}}$

解は後退代入によって得られます。

${\displaystyle x_{n}=d'_{n},}$

${\displaystyle x_{i}=d'_{i}-c'_{i}x_{i+1},\quad i=n-1,n-2,\ldots ,1.}$

上記の方法では、元の係数ベクトルは変更されませんが、新しい係数も記録する必要があります。係数ベクトルが変更される可能性がある場合、記録の少ないアルゴリズムは次のようになります

doの場合 ${\displaystyle i=2,3,\dots ,n,}$

${\displaystyle w:={\cfrac {a_{i}}{b_{i-1}}},}$

${\displaystyle b_{i}:=b_{i}-wc_{i-1},}$

${\displaystyle d_{i}:=d_{i}-wd_{i-1},}$

後退代入が続きます

${\displaystyle x_{n}={\cfrac {d_{n}}{b_{n}}},}$

${\displaystyle x_{i}={\cfrac {d_{i}-c_{i}x_{i+1}}{b_{i}}}\quad {\text{for }}i=n-1,n-2,\dots ,1.}$

C関数としての実装では、スクラッチ空間を使用してacの入力を変更しないようにし、再利用できるようにします。

void thomas ( const int X , double x [ restrict X ] , const double a [ restrict X ], const double b [ restrict X ] , const double c [ restrict X ], double scratch [ restrict X ] ) { / * Ax = dを解きます。ここで、 Aはベクトルa、b、cからなる三角行列です 。X = 方程式の数 x[] = 最初に入力dを含み、xを返しますインデックスは [0, ..., X - 1] です。 a[] = 下対角線、[1, ..., X - 1] からインデックスされます。 b[] = 主対角線、[0, ..., X - 1] からインデックスされます 。 c[] = 上対角線、[0, ..., X - 2] からインデックスされます。 scratch[] = 呼び出し元によって提供される長さ X のスクラッチ スペース。a、b、c を const にすることができます。 この例では実行されません。手動による高価な共通部分式の削除 */ scratch [ 0 ] = c [ 0 ] / b [ 0 ]; x [ 0 ] = x [ 0 ] / b [ 0 ];                                  /* 1からX - 1まで（両端を含む）のループ */ for ( int ix = 1 ; ix < X ; ix ++ ) { if ( ix < X -1 ){ scratch [ ix ] = c [ ix ] / ( b [ ix ] - a [ ix ] * scratch [ ix - 1 ]); } x [ ix ] = ( x [ ix ] - a [ ix ] * x [ ix - 1 ]) / ( b [ ix ] - a [ ix ] * scratch [ ix - 1 ]); }                                             /* X - 2から0まで（両端を含む）のループ */ for ( int ix = X - 2 ; ix >= 0 ; ix -- ) x [ ix ] -= scratch [ ix ] * x [ ix + 1 ]; }                  

導出

三重対角行列アルゴリズムの導出は、ガウス消去法の特殊なケースです。

未知数が であり、解くべき方程式が であると仮定します。 ${\displaystyle x_{1},\ldots ,x_{n}}$

${\displaystyle {\begin{alignedat}{4}&&&b_{1}x_{1}&&+c_{1}x_{2}&&=d_{1}\\&a_{i}x_{i-1}&&+b_{i}x_{i}&&+c_{i}x_{i+1}&&=d_{i}\,,\quad i=2,\ldots ,n-1\\&a_{n}x_{n-1}&&+b_{n}x_{n}&&&&=d_{n}\,.\end{alignedat}}}$

2番目の方程式（）を1番目の方程式で次のように修正することを考えます。 ${\displaystyle i=2}$

${\displaystyle ({\mbox{equation 2}})\cdot b_{1}-({\mbox{equation 1}})\cdot a_{2}}$

これは次のようになります。

${\displaystyle (b_{2}b_{1}-c_{1}a_{2})x_{2}+c_{2}b_{1}x_{3}=d_{2}b_{1}-d_{1}a_{2}.}$

2番目の方程式から が消去されていることに注意してください。修正された2番目の方程式に同様の手法を適用すると、3番目の方程式は次のようになります 。 ${\displaystyle x_{1}}$

${\displaystyle (b_{3}(b_{2}b_{1}-c_{1}a_{2})-c_{2}b_{1}a_{3})x_{3}+c_{3}(b_{2}b_{1}-c_{1}a_{2})x_{4}=d_{3}(b_{2}b_{1}-c_{1}a_{2})-(d_{2}b_{1}-d_{1}a_{2})a_{3}.\,}$

今回はが消去されました。この手順を まで繰り返すと、（修正された）方程式には未知数が1つだけ含まれます。これを について解き、次に を使って方程式を解き、すべての未知数が解かれるまでこれを繰り返します ${\displaystyle x_{2}}$ ${\displaystyle n^{th}}$ ${\displaystyle n^{th}}$ ${\displaystyle x_{n}}$ ${\displaystyle (n-1)^{th}}$

明らかに、修正された方程式の係数は、明示的に記述するとますます複雑になります。手順を調べると、修正された係数（チルダで表記）は再帰的に定義できます。

${\displaystyle {\tilde {a}}_{i}=0\,}$

${\displaystyle {\tilde {b}}_{1}=b_{1}\,}$

${\displaystyle {\tilde {b}}_{i}=b_{i}{\tilde {b}}_{i-1}-{\tilde {c}}_{i-1}a_{i}\,}$

${\displaystyle {\tilde {c}}_{1}=c_{1}\,}$

${\displaystyle {\tilde {c}}_{i}=c_{i}{\tilde {b}}_{i-1}\,}$

${\displaystyle {\tilde {d}}_{1}=d_{1}\,}$

${\displaystyle {\tilde {d}}_{i}=d_{i}{\tilde {b}}_{i-1}-{\tilde {d}}_{i-1}a_{i}.\,}$

解法をさらに迅速化するために、ゼロ除算のリスクがない場合は分割できます。新しい修正係数（それぞれプライム記号で表記）は次のようになります。 ${\displaystyle {\tilde {b}}_{i}}$

${\displaystyle a'_{i}=0\,}$

${\displaystyle b'_{i}=1\,}$

${\displaystyle c'_{1}={\frac {c_{1}}{b_{1}}}\,}$

${\displaystyle c'_{i}={\frac {c_{i}}{b_{i}-c'_{i-1}a_{i}}}\,}$

${\displaystyle d'_{1}={\frac {d_{1}}{b_{1}}}\,}$

${\displaystyle d'_{i}={\frac {d_{i}-d'_{i-1}a_{i}}{b_{i}-c'_{i-1}a_{i}}}.\,}$

これにより、上記の元の方程式で定義された同じ未知数と係数を持つ次の方程式が得られます。

${\displaystyle {\begin{array}{lcl}x_{i}+c'_{i}x_{i+1}=d'_{i}\qquad &;&\ i=1,\ldots ,n-1\\x_{n}=d'_{n}\qquad &;&\ i=n.\\\end{array}}\,}$

最後の方程式には未知数が1つだけ含まれています。これを解くと、次の最後の方程式が未知数が1つに減少するため、この後方代入を使用してすべての未知数を求めることができます。

${\displaystyle x_{n}=d'_{n}\,}$

${\displaystyle x_{i}=d'_{i}-c'_{i}x_{i+1}\qquad ;\ i=n-1,n-2,\ldots ,1.}$

変形

状況によっては、特に周期境界条件が関係する場合、三重対角行列系のわずかに摂動された形を解く必要がある場合があります。

$${\displaystyle {\begin{alignedat}{4}&a_{1}x_{n}&&+b_{1}x_{1}&&+c_{1}x_{2}&&=d_{1}\\&a_{i}x_{i-1}&&+b_{i}x_{i}&&+c_{i}x_{i+1}&&=d_{i}\,,\quad i=2,\ldots ,n-1\\&a_{n}x_{n-1}&&+b_{n}x_{n}&&+c_{n}x_{1}&&=d_{n}\,.\end{alignedat}}}$$

この場合、シャーマン・モリソンの公式を使用してガウス消去法の追加操作を回避し、トーマスアルゴリズムを使用することができます。この方法では、入力とスパース補正ベクトルの両方について、修正された非巡回バージョンのシステムを解き、解を組み合わせる必要があります。純粋な三重対角行列アルゴリズムの順方向部分を共有できるため、両方の解を一度に計算すれば効率的に実行できます

次のように表すと、 $${\displaystyle A={\begin{bmatrix}b_{1}&c_{1}&&&a_{1}\\a_{2}&b_{2}&c_{2}&&\\&a_{3}&b_{3}&\ddots &\\&&\ddots &\ddots &c_{n-1}\\c_{n}&&&a_{n}&b_{n}\end{bmatrix}},x={\begin{bmatrix}x_{1}\\x_{2}\\x_{3}\\\vdots \\x_{n}\end{bmatrix}},d={\begin{bmatrix}d_{1}\\d_{2}\\d_{3}\\\vdots \\d_{n}\end{bmatrix}}}$$

解くべき方程式は次のようになります。 $${\displaystyle Ax=d}$$

この場合、係数とは一般的にゼロではないため、それらの存在によりトーマスアルゴリズムを直接適用することはできません。したがって、と を次のように考えることができます。ここで、 は選択するパラメータです。行列Aは として再構成できます。解は次のように得られます。^[4]まず、トーマスアルゴリズムを適用して2つの三角対角方程式系を解きます。 ${\displaystyle a_{1}}$ ${\displaystyle c_{n}}$ ${\displaystyle B\in \mathbb {R} ^{n\times n}}$ ${\displaystyle u,v\in \mathbb {R} ^{n}}$ $${\displaystyle B={\begin{bmatrix}b_{1}-\gamma &c_{1}&&&0\\a_{2}&b_{2}&c_{2}&&\\&a_{3}&b_{3}&\ddots &\\&&\ddots &\ddots &c_{n-1}\\0&&&a_{n}&b_{n}-{\frac {c_{n}a_{1}}{\gamma }}\end{bmatrix}},u={\begin{bmatrix}\gamma \\0\\0\\\vdots \\c_{n}\end{bmatrix}},v={\begin{bmatrix}1\\0\\0\\\vdots \\a_{1}/\gamma \end{bmatrix}}.}$$ ${\displaystyle \gamma \in \mathbb {R} }$ ${\displaystyle A=B+uv^{\mathsf {T}}}$ $${\displaystyle By=d\qquad \qquad Bq=u}$$

次に、シャーマン・モリソンの公式を使用して解xを再構成します。 $${\displaystyle {\begin{aligned}x&=A^{-1}d=(B+uv^{T})^{-1}d=B^{-1}d-{\frac {B^{-1}uv^{T}B^{-1}}{1+v^{T}B^{-1}u}}d=y-{\frac {qv^{T}y}{1+v^{T}q}}\end{aligned}}}$$

C関数 としての実装。これは、acの入力を変更することを避けるためにスクラッチスペースを使用し、再利用できるようにします

void cyclic_thomas ( const int X , double x [ restrict X ], const double a [ restrict X ], const double b [ restrict X ], const double c [ restrict X ], double cmod [ restrict X ], double u [ restrict X ]) {                          /* Ax = v を解きます。ここで、A はベクトル a、b、c からなる巡回三角行列です。 X = 方程式の数 x[] = 最初に入力 v を含み、x を返します。インデックスは [0, ..., X - 1] です。 a[] = 下対角線、[1, ..., X - 1] から規則的にインデックスされ、a[0] は左下隅です。 b[] = 主対角線、[0, ..., X - 1] から規則的にインデックスされ、c[X - 1] は右上隅です。 cmod[], u[] = それぞれ長さ X のスクラッチベクトル */ /* 巡回三角対角系の左下隅と右上隅です。*/ const double alpha = a [ 0 ]; const double alpha = a[0]; const double beta = c[X - 1]; /* arbitrary, but chosen such that division by zero is avoided */ const double gamma = -b[0]; cmod[0] = c[0] / (b[0] - gamma); u[0] = gamma / (b[0] - gamma); x[0] /= (b[0] - gamma); /* loop from 1 to X - 2 inclusive */ for (int ix = 1; ix + 1 < X; ix++) { const double m = 1.0 / (b[ix] - a[ix] * cmod[ix - 1]); cmod[ix] = c[ix] * m; u[ix] = (0.0f - a[ix] * u[ix - 1]) * m; x[ix] = (x[ix] - a[ix] * x[ix - 1]) * m; } /* handle X - 1 */ const double m = 1.0 / ( b [ X - 1 ] - alpha * beta / gamma - a [ X - 1 ] * cmod [ X - 2 ]);                       u [ X - 1 ] = ( alpha - a [ X - 1 ] * u [ X - 2 ]) * m ;               x [ X - 1 ] = ( x [ X - 1 ] - a [ X - 1 ] * x [ X - 2 ]) * m ;                 /* X - 2 から 0 まで（両端を含む）ループ */ for ( int ix = X - 2 ; ix >= 0 ; ix -- ) {            u [ ix ] -= cmod [ ix ] * u [ ix + 1 ];       x [ ix ] -= cmod [ ix ] * x [ ix + 1 ];       } const double fact = ( x [ 0 ] + x [ X - 1 ] * alpha / gamma ) / ( 1.0 + u [ 0 ] + u [ X - 1 ] * alpha / gamma );                         /* 0からX - 1まで（両端を含む）ループ */ for ( int ix = 0 ; ix < X ; ix ++ )         x [ ix ] -= fact * u [ ix ];    }

上で考察した三重対角線系のわずかに摂動した形を解く別の方法もある。^[5]次元の2つの補助線形系を考えてみよう。 ${\displaystyle (n-1)\times (n-1)}$ $${\displaystyle {\begin{aligned}\qquad \ \ \ \ \ b_{2}u_{2}+c_{2}u_{3}&=d_{2}\\a_{3}u_{2}+b_{3}u_{3}+c_{3}u_{4}&=d_{3}\\a_{i}u_{i-1}+b_{i}u_{i}+c_{i}u_{i+1}&=d_{i}\\\dots \\a_{n}u_{n-1}+b_{n}u_{n}\qquad &=d_{n}\,.\end{aligned}}\quad i=4,\ldots ,n-1\qquad \qquad {\begin{aligned}\qquad \ \ \ \ \ b_{2}v_{2}+c_{2}v_{3}&=-a_{2}\\a_{3}v_{2}+b_{3}v_{3}+c_{3}v_{4}&=0\\a_{i}u_{i-1}+b_{i}u_{i}+c_{i}u_{i+1}&=0\\\dots \\a_{n}v_{n-1}+b_{n}v_{n}\qquad &=-c_{n}\,.\end{aligned}}\quad i=4,\ldots ,n-1}$$

便宜上、さらにとを定義します。これで、2つの補助三角方程式系にトーマスアルゴリズムを適用することで、解を求めることができます。 ${\displaystyle u_{1}=0}$ ${\displaystyle v_{1}=1}$ ${\displaystyle \{u_{2},u_{3}\dots ,u_{n}\}}$ ${\displaystyle \{v_{2},v_{3}\dots ,v_{n}\}}$

解は次のように表すことができます。 ${\displaystyle \{x_{1},x_{2}\dots ,x_{n}\}}$ $${\displaystyle x_{i}=u_{i}+x_{1}v_{i}\qquad i=1,2,\dots ,n}$$

実際、2番目の補助方程式系の各方程式に を掛け、1番目の補助方程式系の対応する方程式と加算し、表現 を使用すると、方程式番号 ${\displaystyle x_{1}}$ ${\displaystyle x_{i}=u_{i}+x_{1}v_{i}}$2からnが満たされていることがすぐにわかります。残っているのは方程式番号1を満たすことだけです。そのためには、 と の式を考え、とを元の方程式系の最初の方程式に代入します。これにより、 の1つのスカラー方程式が得られます。 ${\displaystyle i=2}$ ${\displaystyle i=n}$ ${\displaystyle x_{2}=u_{2}+x_{1}v_{2}}$ ${\displaystyle x_{n}=u_{n}+x_{1}v_{n}}$ ${\displaystyle x_{1}}$ $${\displaystyle b_{1}x_{1}+c_{1}(u_{2}+x_{1}v_{2})+a_{1}(u_{n}+x_{1}v_{n})=d_{1}}$$

したがって、次式が得られます。 $${\displaystyle x_{1}={\frac {d_{1}-a_{1}u_{n}-c_{1}u_{2}}{b_{1}+a_{1}v_{n}+c_{1}v_{2}}}}$$

C関数としての実装では、スクラッチスペースを使用してacへの入力を変更しないようにし、再利用できるようにします

void cyclic_thomas ( const int X , double x [ restrict X ], const double a [ restrict X ], const double b [ restrict X ], const double c [ restrict X ], double cmod [ restrict X ], double v [ restrict X ]) {                          /* まず、長さ X - 1 の連立方程式を、右辺2つについて解きます。ix == 0 は無視します */ cmod [ 1 ] = c [ 1 ] / b [ 1 ];     v [ 1 ] = - a [ 1 ] / b [ 1 ];     x [ 1 ] = x [ 1 ] / b [ 1 ];     /* 2 から X - 1 まで（両端を含む）ループします */ for ( int ix = 2 ; ix < X - 1 ; ix ++ ) {            定数倍精度浮動小数点数型(const double m) = 1.0 / ( b [ ix ] - a [ ix ] * cmod [ ix - 1 ] );             cmod [ ix ] = c [ ix ] * m ;     v [ ix ] = ( 0.0f - a [ ix ] * v [ ix - 1 ]) * m ;           x [ ix ] = ( x [ ix ] - a [ ix ] * x [ ix - 1 ]) * m ;           } /* X - 1 を扱う */ const double m = 1.0 / ( b [ X - 1 ] - a [ X - 1 ] * cmod [ X - 2 ] );                 cmod [ X - 1 ] = c [ X - 1 ] * m ;         v [ X - 1 ] = ( - c [ 0 ] - a [ X - 1 ] * v [ X - 2 ]) * m ;               x [ X - 1 ] = ( x [ X - 1 ] - a [ X - 1 ] * x [ X - 2 ]) * m ;                 /* X - 2 から 1 まで（両端を含む）ループ */ for ( int ix = X - 2 ; ix >= 1 ; ix -- ) {            v [ ix ] -= cmod [ ix ] * v [ ix + 1 ];       x [ ix ] -= cmod [ ix ] * x [ ix + 1 ];       } x [ 0 ] = ( x [ 0 ] - a [ 0 ] * x [ X - 1 ] - c [ 0 ] * x [ 1 ]) / ( b [ 0 ] + a [ 0 ] * v [ X - 1 ] + c [ 0 ] * v [ 1 ]);                         /* 1 から X - 1 まで（両端を含む）ループ */ for ( int ix = 1 ; ix < X ; ix ++ )         x [ ix ] += x [ 0 ] * v [ ix ];    }

どちらの場合も、解くべき補助系は真に三重対角系であるため、系を解く全体的な計算量は、系の次元n、つまり算術演算に対して線形のままです ${\displaystyle Ax=d}$ ${\displaystyle O(n)}$

他の状況では、連立方程式はブロック三角行列（ブロック行列を参照）であり、より小さな部分行列が上記の行列システムの個々の要素として配置される場合があります（例：2次元ポアソン問題）。このような状況のために、簡略化されたガウス消去法が開発されています。^[6]

Alfio Quarteroni 、Sacco、Saleriによる教科書『Numerical Mathematics』には、一部の除算を回避（代わりに乗算を使用）するアルゴリズムの修正版が掲載されており、これは一部のコンピュータアーキテクチャで有益です。

GPUを含む多くのベクトルおよび並列アーキテクチャ向けに、並列三角ソルバーが公開されています^{[7] [8]}

並列三角ソルバーとブロック三角ソルバーの詳細な説明については、 ^{[9]を参照してください。}

参考文献

1. Pradip Niyogi (2006).数値流体力学入門. ピアソン・エデュケーション・インディア. p. 76. ISBN 978-81-7758-764-7。
2. Biswa Nath Datta (2010).数値線形代数とその応用 第2版. SIAM. p. 162. ISBN 978-0-89871-765-5。
3. Nicholas J. Higham (2002).数値アルゴリズムの精度と安定性：第2版. SIAM. p. 175. ISBN 978-0-89871-802-7。
4. Batista, Milan; Ibrahim Karawia, Abdel Rahman A. (2009). 「巡回ブロック三対角線および巡回ブロック五対角線線形連立方程式の解法におけるシャーマン・モリソン・ウッドベリー公式の使用」 .応用数学と計算. 210 (2): 558– 563. doi :10.1016/j.amc.2009.01.003. ISSN  0096-3003
5. Ryaben'kii, Victor S.; Tsynkov, Semyon V. (2006-11-02)、「序論」、数値解析の理論的入門、Chapman and Hall/CRC、pp.  1-19、doi :10.1201/9781420011166-1、ISBN 978-0-429-14339-7、2022年5月25日閲覧
6. Quarteroni, Alfio ; Sacco, Riccardo; Saleri, Fausto (2007). 「Section 3.8」、数値数学、Springer、ニューヨーク、ISBN 978-3-540-34658-6。
7. Chang, L.-W.; Hwu, W,-M. (2014). 「GPU上で三角行列ソルバーを実装するためのガイド」V. Kidratenko (編). GPUによる数値計算. Springer. ISBN 978-3-319-06548-9。{{cite conference}}: CS1 maint: multiple names: authors list (link)
8. Venetis, IE; Kouris, A.; Sobczyk, A.; Gallopoulos, E.; Sameh, A. (2015). 「GPUアーキテクチャ向けのギブンズ回転に基づく直接三角関数ソルバー」Parallel Computing . 49 : 101–116 . doi :10.1016/j.parco.2015.03.008
9. Gallopoulos, E.; Philippe, B.; Sameh, AH (2016). 「第5章行列計算における並列処理」 . Springer. ISBN 978-94-017-7188-7。

• Conte, SD; de Boor, C. (1972).初等数値解析. McGraw-Hill, New York. ISBN 0070124469。
• この記事には、CFD-WikiのGFDLライセンスに基づく記事「Tridiagonal_matrix_algorithm_-_TDMA_(Thomas_algorithm)」のテキストが組み込まれています。

• Press, WH; Teukolsky, SA; Vetterling, WT; Flannery, BP (2007). 「第2.4節」. Numerical Recipes: The Art of Scientific Computing (第3版). ニューヨーク: Cambridge University Press. ISBN 978-0-521-88068-8。

Retrieved from "https://en.wikipedia.org/w/index.php?title=Tridiagonal_matrix_algorithm&oldid=1316098762"