Алгоритм трехдиагональной матрицы

Улучшенное сокращение для определенных матриц

В числовой линейной алгебре алгоритм трехдиагональной матрицы , также известный как алгоритм Томаса (названный в честь Ллевеллина Томаса ), является упрощенной формой метода Гаусса , который может быть использован для решения трехдиагональных систем уравнений . Трехдиагональная система для 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}}.}$

Для таких систем решение может быть получено в операциях вместо требуемого исключения Гаусса . Первая развертка устраняет 's, а затем (сокращенная) обратная подстановка производит решение. Примеры таких матриц обычно возникают из дискретизации уравнения Пуассона 1D и естественной кубической сплайновой интерполяции . ${\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.}$

Метод выше не изменяет исходные векторы коэффициентов, но также должен отслеживать новые коэффициенты. Если векторы коэффициентов могут быть изменены, то алгоритм с меньшим учетом выглядит так:

Для делать ${\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[] = изначально содержит входные данные v и возвращает x. индексируется из [0, ..., X - 1]  a[] = поддиагональ, индексируется из [1, ..., X - 1]  b[] = главная диагональ, индексируется из [0, ..., X - 1]  c[] = наддиагональ, индексируется из [0, ..., X - 2]  scratch[] = пространство для царапин длиной X, предоставленное вызывающей стороной, позволяющее сделать a, b, c константами  в этом примере не выполняется: ручное дорогостоящее исключение общей подвыражения  */ scratch [ 0 ] = c [ 0 ] / b [ 0 ]; х [ 0 ] = х [ 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}}}$

Рассмотрим изменение второго ( ) уравнения с помощью первого уравнения следующим образом: ${\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}.}$

Обратите внимание, что было исключено из второго уравнения. Использование аналогичной тактики с измененным вторым уравнением в третьем уравнении дает: ${\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}.\,}$

Это время было устранено. Если эту процедуру повторять до строки; (измененное) уравнение будет включать только одно неизвестное, . Это может быть решено для и затем использовано для решения уравнения, и так далее, пока все неизвестные не будут решены для. ${\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}}\,}$

Последнее уравнение содержит только одно неизвестное. Решая его, мы сводим следующее последнее уравнение к одному неизвестному, так что эта обратная подстановка может быть использована для нахождения всех неизвестных:

${\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}$$

В этом случае коэффициенты и , вообще говоря, не равны нулю, поэтому их наличие не позволяет напрямую применить алгоритм Томаса. Поэтому мы можем рассматривать и следующим образом: Где — параметр, который необходимо выбрать. Матрица может быть восстановлена ​​как . Тогда решение получается следующим образом: ^[4] сначала решаем две трехдиагональные системы уравнений, применяя алгоритм Томаса: ${\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}$ ${\displaystyle A=B+uv^{\mathsf {T}}}$ $${\displaystyle By=d\qquad \qquad Bq=u}$$

Затем реконструируем решение, используя формулу Шермана-Моррисона : ${\displaystyle 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[] = наддиагональ, регулярно индексируется из [0, ..., X - 2], c[X - 1] - верхний правый  cmod[], u[] = векторные царапины, каждый из которых имеет длину X  */                           /* нижний левый и верхний правый углы циклической трехдиагональной системы соответственно */ const double alpha = a [ 0 ]; const double beta = c [ X - 1 ];             /* произвольно, но выбрано так, чтобы избежать деления на ноль */ const double gamma = - b [ 0 ];      cmod [ 0 ] = альфа / ( b [ 0 ] - гамма ); u [ 0 ] = гамма / ( b [ 0 ] - гамма ); x [ 0 ] /= ( b [ 0 ] - гамма );                   /* цикл от 1 до X - 2 включительно */ 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 ; }                                                      /* обработка X - 1 */ const double m = 1.0 / ( b [ X - 1 ] - альфа * бета / гамма - бета * cmod [ X - 2 ]); u [ X - 1 ] = ( альфа - 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 ] * бета / гамма ) / ( 1,0 + u [ 0 ] + u [ X - 1 ] * бета / гамма );                         /* цикл от 0 до X - 1 включительно */ for ( int ix = 0 ; ix < X ; ix ++ ) x [ ix ] -= fact * u [ ix ]; }              

Существует также другой способ решения слегка возмущенной формы трехдиагональной системы, рассмотренной выше. ^[5] Рассмотрим две вспомогательные линейные системы размерности : ${\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}$$

Для удобства дополнительно определим и . Теперь мы можем найти решения и , применив алгоритм Томаса к двум вспомогательным трехдиагональным системам. ${\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}$$

Действительно, умножая каждое уравнение второй вспомогательной системы на , складывая с соответствующим уравнением первой вспомогательной системы и используя представление , мы сразу видим, что уравнения с номерами по исходной системы удовлетворяются; осталось только удовлетворить уравнению с номером . Для этого рассмотрим формулы для и и подставим и в первое уравнение исходной системы. Это дает одно скалярное уравнение для : ${\displaystyle x_{1}}$ ${\displaystyle x_{i}=u_{i}+x_{1}v_{i}}$ ${\displaystyle 2}$ ${\displaystyle n}$ ${\displaystyle 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 для двух правых частей, игнорируя 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 ]; }              

В обоих случаях вспомогательные системы, которые необходимо решить, являются действительно трехдиагональными, поэтому общая вычислительная сложность решения системы остается линейной по отношению к размерности системы , то есть арифметическим операциям. ${\displaystyle Ax=d}$ ${\displaystyle n}$ ${\displaystyle O(n)}$

В других ситуациях система уравнений может быть блочно-трехдиагональной (см. блочную матрицу ), с меньшими подматрицами, расположенными как отдельные элементы в вышеприведенной матричной системе (например, двумерная задача Пуассона ). Для этих ситуаций были разработаны упрощенные формы гауссовского исключения. ^[6]

В учебнике «Численная математика» Альфио Квартерони , Сакко и Салери приводится модифицированная версия алгоритма, которая позволяет избежать некоторых делений (используя вместо них умножения), что полезно на некоторых архитектурах компьютеров.

Параллельные трехдиагональные решатели были опубликованы для многих векторных и параллельных архитектур, включая графические процессоры ^{[7] [8]}

Для подробного рассмотрения параллельных трехдиагональных и блочных трехдиагональных решателей см. ^[9]

Ссылки

1. Прадип Нийоги (2006). Введение в вычислительную гидродинамику . Pearson Education India. стр. 76. ISBN 978-81-7758-764-7.
2. Biswa Nath Datta (2010). Численная линейная алгебра и приложения, второе издание . SIAM. стр. 162. ISBN 978-0-89871-765-5.
3. Николас Дж. Хайэм (2002). Точность и устойчивость численных алгоритмов: Второе издание . SIAM. стр. 175. ISBN 978-0-89871-802-7.
4. Батиста, Милан; Ибрагим Каравия, Абдель Рахман А. (2009). «Использование формулы Шермана–Моррисона–Вудбери для решения циклических блочно-трехдиагональных и циклических блочно-пятидиагональных линейных систем уравнений». Прикладная математика и вычисления . 210 (2): 558–563. doi :10.1016/j.amc.2009.01.003. ISSN  0096-3003.
5. Рябенький, Виктор С.; Цынков, Семен В. (2006-11-02), "Введение", Теоретическое введение в численный анализ , Chapman and Hall/CRC, стр. 1–19, doi :10.1201/9781420011166-1, ISBN 978-0-429-14339-7, получено 2022-05-25
6. Квартерони, Альфио ; Сакко, Риккардо; Салери, Фаусто (2007). «Раздел 3.8». Численная математика . Спрингер, Нью-Йорк. ISBN 978-3-540-34658-6.
7. Chang, L.-W.; Hwu, W,-M. (2014). "Руководство по реализации трехдиагональных решателей на графических процессорах". В V. Kidratenko (ред.). Численные вычисления на графических процессорах . Springer. ISBN 978-3-319-06548-9.{{cite conference}}: CS1 maint: multiple names: authors list (link)
8. Венетис, IE; Курис, А.; Собчик А.; Галлопулос, Э.; Самех, А. (2015). «Прямой трехдиагональный решатель, основанный на вращении Гивенса для архитектур графических процессоров». Параллельные вычисления . 49 : 101–116. doi :10.1016/j.parco.2015.03.008.
9. Галлопулос, Э.; Филипп, Б.; Самех, А.Х. (2016). "Глава 5". Параллелизм в матричных вычислениях . Springer. ISBN 978-94-017-7188-7.

• Конте, SD; де Бур, C. (1972). Элементарный численный анализ . McGraw-Hill, Нью-Йорк. ISBN 0070124469.
• В данной статье использован текст из статьи Tridiagonal_matrix_algorithm_-_TDMA_(Thomas_algorithm) на CFD-Wiki, которая находится под лицензией GFDL .

• 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.