Алгоритм SPIKE

Алгоритм SPIKE — это гибридный параллельный решатель для полосовых линейных систем, разработанный Эриком Полицци и Ахмедом Самехом. ^[1] ^ ^[2]

Обзор

Алгоритм SPIKE имеет дело с линейной системой AX = F , где A — полосчатая система. $n\times n$ матрица полосы пропускания намного меньше, чем $n$, а F представляет собой $n\times s$ матрица, содержащая $s$ правые стороны. Он разделен на этап предварительной обработки и этап постобработки.

Этап предварительной обработки

На этапе предварительной обработки линейная система AX = F разбивается на блок трехдиагональной формы

${\begin{bmatrix}{\boldsymbol {A}}_{1}&{\boldsymbol {B}}_{1}\\{\boldsymbol {C}}_{2}&{\boldsymbol {A}}_{2}&{\boldsymbol {B}}_{2}\\&\ddots &\ddots &\ddots \\&&{\boldsymbol {C}}_{p-1}&{\boldsymbol {A}}_{p-1}&{\boldsymbol {B}}_{p-1}\\&&&{\boldsymbol {C}}_{p}&{\boldsymbol {A}}_{p}\end{bmatrix}}{\begin{bmatrix}{\boldsymbol {X}}_{1}\\{\boldsymbol {X}}_{2}\\\vdots \\{\boldsymbol {X}}_{p-1}\\{\boldsymbol {X}}_{p}\end{bmatrix}}={\begin{bmatrix}{\boldsymbol {F}}_{1}\\{\boldsymbol {F}}_{2}\\\vdots \\{\boldsymbol {F}}_{p-1}\\{\boldsymbol {F}}_{p}\end{bmatrix}}.$

Предположим пока, что диагональные блоки A _j ( j = 1,..., p с p ≥ 2 ) неособы . Определите блочную диагональную матрицу

D = диаг( А ₁ ,..., А _п ) ,

тогда D также неособа. Умножение слева D ⁻¹ обеим сторонам системы дает

${\begin{bmatrix}{\boldsymbol {I}}&{\boldsymbol {V}}_{1}\\{\boldsymbol {W}}_{2}&{\boldsymbol {I}}&{\boldsymbol {V}}_{2}\\&\ddots &\ddots &\ddots \\&&{\boldsymbol {W}}_{p-1}&{\boldsymbol {I}}&{\boldsymbol {V}}_{p-1}\\&&&{\boldsymbol {W}}_{p}&{\boldsymbol {I}}\end{bmatrix}}{\begin{bmatrix}{\boldsymbol {X}}_{1}\\{\boldsymbol {X}}_{2}\\\vdots \\{\boldsymbol {X}}_{p-1}\\{\boldsymbol {X}}_{p}\end{bmatrix}}={\begin{bmatrix}{\boldsymbol {G}}_{1}\\{\boldsymbol {G}}_{2}\\\vdots \\{\boldsymbol {G}}_{p-1}\\{\boldsymbol {G}}_{p}\end{bmatrix}},$

которая должна быть решена на этапе постобработки. Умножение слева на D ⁻¹ эквивалентно решению $p$ системы вида

A _j [ V _j W _j G _j ] знак равно [ B _j C _j F _j ]

(опуская W ₁ и C ₁ для $j=1$, а также V _p и B _p для $j=p$), что может осуществляться параллельно.

Из-за полосовой природы A только несколько крайних левых столбцов каждого V _j и несколько крайних правых столбцов каждого W _j могут быть ненулевыми. Эти столбцы называются шипами .

Этап постобработки

Без ограничения общности будем считать, что каждый спайк содержит ровно $m$ столбцы ( $m$ намного меньше, чем $n$) (при необходимости дополните пик столбцами нулей). Разделите пики во всех V _j и W _j на

${\begin{bmatrix}{\boldsymbol {V}}_{j}^{(t)}\\{\boldsymbol {V}}_{j}'\\{\boldsymbol {V}}_{j}^{(b)}\end{bmatrix}}$ и ${\begin{bmatrix}{\boldsymbol {W}}_{j}^{(t)}\\{\boldsymbol {W}}_{j}'\\{\boldsymbol {W}}_{j}^{(b)}\\\end{bmatrix}}$

где V   ( т )
j  , V  ( b )
j   , Вт   ( т )
j   и W   ( б )
j   имеют размеры $m\times m$. Аналогично разбейте все X _j и G _j на

${\begin{bmatrix}{\boldsymbol {X}}_{j}^{(t)}\\{\boldsymbol {X}}_{j}'\\{\boldsymbol {X}}_{j}^{(b)}\end{bmatrix}}$ и ${\begin{bmatrix}{\boldsymbol {G}}_{j}^{(t)}\\{\boldsymbol {G}}_{j}'\\{\boldsymbol {G}}_{j}^{(b)}\\\end{bmatrix}}.$

Заметим, что систему, полученную на этапе предварительной обработки, можно свести к блочной пятидиагональной системе гораздо меньшего размера (напомним, что $m$ намного меньше, чем $n$)

${\begin{bmatrix}{\boldsymbol {I}}_{m}&{\boldsymbol {0}}&{\boldsymbol {V}}_{1}^{(t)}\\{\boldsymbol {0}}&{\boldsymbol {I}}_{m}&{\boldsymbol {V}}_{1}^{(b)}&{\boldsymbol {0}}\\{\boldsymbol {0}}&{\boldsymbol {W}}_{2}^{(t)}&{\boldsymbol {I}}_{m}&{\boldsymbol {0}}&{\boldsymbol {V}}_{2}^{(t)}\\&{\boldsymbol {W}}_{2}^{(b)}&{\boldsymbol {0}}&{\boldsymbol {I}}_{m}&{\boldsymbol {V}}_{2}^{(b)}&{\boldsymbol {0}}\\&&\ddots &\ddots &\ddots &\ddots &\ddots \\&&&{\boldsymbol {0}}&{\boldsymbol {W}}_{p-1}^{(t)}&{\boldsymbol {I}}_{m}&{\boldsymbol {0}}&{\boldsymbol {V}}_{p-1}^{(t)}\\&&&&{\boldsymbol {W}}_{p-1}^{(b)}&{\boldsymbol {0}}&{\boldsymbol {I}}_{m}&{\boldsymbol {V}}_{p-1}^{(b)}&{\boldsymbol {0}}\\&&&&&{\boldsymbol {0}}&{\boldsymbol {W}}_{p}^{(t)}&{\boldsymbol {I}}_{m}&{\boldsymbol {0}}\\&&&&&&{\boldsymbol {W}}_{p}^{(b)}&{\boldsymbol {0}}&{\boldsymbol {I}}_{m}\end{bmatrix}}{\begin{bmatrix}{\boldsymbol {X}}_{1}^{(t)}\\{\boldsymbol {X}}_{1}^{(b)}\\{\boldsymbol {X}}_{2}^{(t)}\\{\boldsymbol {X}}_{2}^{(b)}\\\vdots \\{\boldsymbol {X}}_{p-1}^{(t)}\\{\boldsymbol {X}}_{p-1}^{(b)}\\{\boldsymbol {X}}_{p}^{(t)}\\{\boldsymbol {X}}_{p}^{(b)}\end{bmatrix}}={\begin{bmatrix}{\boldsymbol {G}}_{1}^{(t)}\\{\boldsymbol {G}}_{1}^{(b)}\\{\boldsymbol {G}}_{2}^{(t)}\\{\boldsymbol {G}}_{2}^{(b)}\\\vdots \\{\boldsymbol {G}}_{p-1}^{(t)}\\{\boldsymbol {G}}_{p-1}^{(b)}\\{\boldsymbol {G}}_{p}^{(t)}\\{\boldsymbol {G}}_{p}^{(b)}\end{bmatrix}}{\text{,}}$

которую мы называем приведенной системой и обозначаем S̃X̃ = G̃ .

Как только все X   ( t )
j   и X   ( б )
j   найдены, все X ′ _j можно восстановить с идеальным параллелизмом с помощью

${\begin{cases}{\boldsymbol {X}}_{1}'={\boldsymbol {G}}_{1}'-{\boldsymbol {V}}_{1}'{\boldsymbol {X}}_{2}^{(t)}{\text{,}}\\{\boldsymbol {X}}_{j}'={\boldsymbol {G}}_{j}'-{\boldsymbol {V}}_{j}'{\boldsymbol {X}}_{j+1}^{(t)}-{\boldsymbol {W}}_{j}'{\boldsymbol {X}}_{j-1}^{(b)}{\text{,}}&j=2,\ldots ,p-1{\text{,}}\\{\boldsymbol {X}}_{p}'={\boldsymbol {G}}_{p}'-{\boldsymbol {W}}_{p}{\boldsymbol {X}}_{p-1}^{(b)}{\text{.}}\end{cases}}$

SPIKE как полиалгоритмический решатель полосовых линейных систем

Несмотря на то, что алгоритм SPIKE логически разделен на два этапа, в вычислительном отношении он состоит из трех этапов:

1. факторизация диагональных блоков,
2. вычисление спайков,
3. решение приведенной системы.

Каждый из этих этапов может быть реализован несколькими способами, что допускает множество вариантов. Двумя примечательными вариантами являются рекурсивный алгоритм SPIKE для случаев с недиагональным доминированием и усеченный алгоритм SPIKE для случаев с диагональным доминированием. В зависимости от варианта система может быть решена точно или приближенно. В последнем случае SPIKE используется как предобуславливатель для итерационных схем, таких как методы подпространств Крылова и итеративное уточнение .

Рекурсивный ШИП

Этап предварительной обработки

Первым шагом этапа предварительной обработки является факторизация диагональных блоков A _j . можно использовать LAPACK Для численной стабильности XGBTRF подпрограммы для LU-факторизации их с частичным поворотом. В качестве альтернативы их можно также факторизовать без частичного поворота, но с помощью стратегии «диагонального повышения». Последний метод решает проблему сингулярных диагональных блоков.

Конкретно стратегия диагонального повышения выглядит следующим образом. Пусть 0 _ε обозначает настраиваемый «ноль машины». На каждом этапе LU-факторизации мы требуем, чтобы опорная точка удовлетворяла условию

|поворот| > 0 _ε ‖ А ‖ ₁ .

Если опорная точка не удовлетворяет условию, она повышается за счет

$\mathrm {pivot} ={\begin{cases}\mathrm {pivot} +\epsilon \lVert {\boldsymbol {A}}_{j}\rVert _{1}&{\text{if }}\mathrm {pivot} \geq 0{\text{,}}\\\mathrm {pivot} -\epsilon \lVert {\boldsymbol {A}}_{j}\rVert _{1}&{\text{if }}\mathrm {pivot} <0\end{cases}}$

где ε машины — положительный параметр, зависящий от единицы округления , и факторизация продолжается с усиленным центром вращения. Этого можно достичь с помощью модифицированных версий ScaLAPACK . XDBTRF процедуры. После факторизации диагональных блоков пики вычисляются и передаются на этап постобработки.

Этап постобработки

Двухраздельный случай

В двухраздельном случае, т. е. когда p = 2 , приведенная система S̃X̃ = G̃ имеет вид

${\begin{bmatrix}{\boldsymbol {I}}_{m}&{\boldsymbol {0}}&{\boldsymbol {V}}_{1}^{(t)}\\{\boldsymbol {0}}&{\boldsymbol {I}}_{m}&{\boldsymbol {V}}_{1}^{(b)}&{\boldsymbol {0}}\\{\boldsymbol {0}}&{\boldsymbol {W}}_{2}^{(t)}&{\boldsymbol {I}}_{m}&{\boldsymbol {0}}\\&{\boldsymbol {W}}_{2}^{(b)}&{\boldsymbol {0}}&{\boldsymbol {I}}_{m}\end{bmatrix}}{\begin{bmatrix}{\boldsymbol {X}}_{1}^{(t)}\\{\boldsymbol {X}}_{1}^{(b)}\\{\boldsymbol {X}}_{2}^{(t)}\\{\boldsymbol {X}}_{2}^{(b)}\end{bmatrix}}={\begin{bmatrix}{\boldsymbol {G}}_{1}^{(t)}\\{\boldsymbol {G}}_{1}^{(b)}\\{\boldsymbol {G}}_{2}^{(t)}\\{\boldsymbol {G}}_{2}^{(b)}\end{bmatrix}}{\text{.}}$

Из центра можно извлечь еще меньшую систему:

${\begin{bmatrix}{\boldsymbol {I}}_{m}&{\boldsymbol {V}}_{1}^{(b)}\\{\boldsymbol {W}}_{2}^{(t)}&{\boldsymbol {I}}_{m}\end{bmatrix}}{\begin{bmatrix}{\boldsymbol {X}}_{1}^{(b)}\\{\boldsymbol {X}}_{2}^{(t)}\end{bmatrix}}={\begin{bmatrix}{\boldsymbol {G}}_{1}^{(b)}\\{\boldsymbol {G}}_{2}^{(t)}\end{bmatrix}}{\text{,}}$

которое можно решить с помощью блочной LU-факторизации

${\begin{bmatrix}{\boldsymbol {I}}_{m}&{\boldsymbol {V}}_{1}^{(b)}\\{\boldsymbol {W}}_{2}^{(t)}&{\boldsymbol {I}}_{m}\end{bmatrix}}={\begin{bmatrix}{\boldsymbol {I}}_{m}\\{\boldsymbol {W}}_{2}^{(t)}&{\boldsymbol {I}}_{m}\end{bmatrix}}{\begin{bmatrix}{\boldsymbol {I}}_{m}&{\boldsymbol {V}}_{1}^{(b)}\\&{\boldsymbol {I}}_{m}-{\boldsymbol {W}}_{2}^{(t)}{\boldsymbol {V}}_{1}^{(b)}\end{bmatrix}}{\text{.}}$

Однажды Х   ( б )
1   и Х   ( т )
2  найдено , X   ( t )
1   и Х   ( б )
2   можно вычислить через

Икс   ( т )
1   знак равно г   ( т )
1   − V  ( t )
1   х   ( т )
2  ,

Икс   ( б )
2   знак равно г   ( б )
2   - Вт   ( б )
2   х   ( б )
1  .

Случай с несколькими разделами

Предположим, что p — степень двойки, т. е. p = 2. ^д . Рассмотрим блочную диагональную матрицу

D̃ ₁ = диаг( D̃   [1]
1   ,..., Д̃   [1]
п /2   )

где

${\boldsymbol {\tilde {D}}}_{k}^{[1]}={\begin{bmatrix}{\boldsymbol {I}}_{m}&{\boldsymbol {0}}&{\boldsymbol {V}}_{2k-1}^{(t)}\\{\boldsymbol {0}}&{\boldsymbol {I}}_{m}&{\boldsymbol {V}}_{2k-1}^{(b)}&{\boldsymbol {0}}\\{\boldsymbol {0}}&{\boldsymbol {W}}_{2k}^{(t)}&{\boldsymbol {I}}_{m}&{\boldsymbol {0}}\\&{\boldsymbol {W}}_{2k}^{(b)}&{\boldsymbol {0}}&{\boldsymbol {I}}_{m}\end{bmatrix}}$

для к = 1,..., р /2 . Обратите внимание, что D̃ ₁ по существу состоит из диагональных блоков порядка 4 m, извлеченных из S̃ . Теперь мы факторизуем S как

S̃ знак равно D̃ ₁ S̃ ₂ .

Новая матрица S̃ ₂ имеет вид

${\begin{bmatrix}{\boldsymbol {I}}_{3m}&{\boldsymbol {0}}&{\boldsymbol {V}}_{1}^{[2](t)}\\{\boldsymbol {0}}&{\boldsymbol {I}}_{m}&{\boldsymbol {V}}_{1}^{[2](b)}&{\boldsymbol {0}}\\{\boldsymbol {0}}&{\boldsymbol {W}}_{2}^{[2](t)}&{\boldsymbol {I}}_{m}&{\boldsymbol {0}}&{\boldsymbol {V}}_{2}^{[2](t)}\\&{\boldsymbol {W}}_{2}^{[2](b)}&{\boldsymbol {0}}&{\boldsymbol {I}}_{3m}&{\boldsymbol {V}}_{2}^{[2](b)}&{\boldsymbol {0}}\\&&\ddots &\ddots &\ddots &\ddots &\ddots \\&&&{\boldsymbol {0}}&{\boldsymbol {W}}_{p/2-1}^{[2](t)}&{\boldsymbol {I}}_{3m}&{\boldsymbol {0}}&{\boldsymbol {V}}_{p/2-1}^{[2](t)}\\&&&&{\boldsymbol {W}}_{p/2-1}^{[2](b)}&{\boldsymbol {0}}&{\boldsymbol {I}}_{m}&{\boldsymbol {V}}_{p/2-1}^{[2](b)}&{\boldsymbol {0}}\\&&&&&{\boldsymbol {0}}&{\boldsymbol {W}}_{p/2}^{[2](t)}&{\boldsymbol {I}}_{m}&{\boldsymbol {0}}\\&&&&&&{\boldsymbol {W}}_{p/2}^{[2](b)}&{\boldsymbol {0}}&{\boldsymbol {I}}_{3m}\end{bmatrix}}{\text{.}}$

Его структура очень похожа на структуру S̃ ₂ , отличается только количеством шипов и их высотой (их ширина остается неизменной при m ). Таким образом, аналогичный шаг факторизации можно выполнить для S̃ _2, чтобы получить

С̃ ₂ = Д̃ ₂ С̃ ₃

и

S̃ знак равно D̃ ₁ D̃ ₂ S̃ ₃ .

Такие этапы факторизации могут выполняться рекурсивно. После d − 1 шагов мы получаем факторизацию

S̃ знак равно D̃ ₁ ⋯ D̃ _{d −1} S̃ _d ,

где у S̃ _d всего два шипа. Приведенная система будет тогда решена через

X̃ = S̃   −1
d   D̃   −1
d −1   ⋯ D̃   −1
1   г. ​

Техника блочной LU-факторизации в случае с двумя разделами может использоваться для обработки шагов решения, включающих D̃ ₁ , ..., D̃ _{d −1} и S̃ _d, поскольку они по существу решают несколько независимых систем обобщенных двухраздельных форм.

Обобщение на случаи, когда p не является степенью двойки, почти тривиально.

Усеченный ШИП

Когда A является диагонально-доминантным, в приведенной системе

${\begin{bmatrix}{\boldsymbol {I}}_{m}&{\boldsymbol {0}}&{\boldsymbol {V}}_{1}^{(t)}\\{\boldsymbol {0}}&{\boldsymbol {I}}_{m}&{\boldsymbol {V}}_{1}^{(b)}&{\boldsymbol {0}}\\{\boldsymbol {0}}&{\boldsymbol {W}}_{2}^{(t)}&{\boldsymbol {I}}_{m}&{\boldsymbol {0}}&{\boldsymbol {V}}_{2}^{(t)}\\&{\boldsymbol {W}}_{2}^{(b)}&{\boldsymbol {0}}&{\boldsymbol {I}}_{m}&{\boldsymbol {V}}_{2}^{(b)}&{\boldsymbol {0}}\\&&\ddots &\ddots &\ddots &\ddots &\ddots \\&&&{\boldsymbol {0}}&{\boldsymbol {W}}_{p-1}^{(t)}&{\boldsymbol {I}}_{m}&{\boldsymbol {0}}&{\boldsymbol {V}}_{p-1}^{(t)}\\&&&&{\boldsymbol {W}}_{p-1}^{(b)}&{\boldsymbol {0}}&{\boldsymbol {I}}_{m}&{\boldsymbol {V}}_{p-1}^{(b)}&{\boldsymbol {0}}\\&&&&&{\boldsymbol {0}}&{\boldsymbol {W}}_{p}^{(t)}&{\boldsymbol {I}}_{m}&{\boldsymbol {0}}\\&&&&&&{\boldsymbol {W}}_{p}^{(b)}&{\boldsymbol {0}}&{\boldsymbol {I}}_{m}\end{bmatrix}}{\begin{bmatrix}{\boldsymbol {X}}_{1}^{(t)}\\{\boldsymbol {X}}_{1}^{(b)}\\{\boldsymbol {X}}_{2}^{(t)}\\{\boldsymbol {X}}_{2}^{(b)}\\\vdots \\{\boldsymbol {X}}_{p-1}^{(t)}\\{\boldsymbol {X}}_{p-1}^{(b)}\\{\boldsymbol {X}}_{p}^{(t)}\\{\boldsymbol {X}}_{p}^{(b)}\end{bmatrix}}={\begin{bmatrix}{\boldsymbol {G}}_{1}^{(t)}\\{\boldsymbol {G}}_{1}^{(b)}\\{\boldsymbol {G}}_{2}^{(t)}\\{\boldsymbol {G}}_{2}^{(b)}\\\vdots \\{\boldsymbol {G}}_{p-1}^{(t)}\\{\boldsymbol {G}}_{p-1}^{(b)}\\{\boldsymbol {G}}_{p}^{(t)}\\{\boldsymbol {G}}_{p}^{(b)}\end{bmatrix}}{\text{,}}$

блоки V   ( t )
j   и W   ( б )
j   часто пренебрежимо малы. Без них приведенная система становится блочно-диагональной.

${\begin{bmatrix}{\boldsymbol {I}}_{m}\\&{\boldsymbol {I}}_{m}&{\boldsymbol {V}}_{1}^{(b)}\\&{\boldsymbol {W}}_{2}^{(t)}&{\boldsymbol {I}}_{m}\\&&&{\boldsymbol {I}}_{m}&{\boldsymbol {V}}_{2}^{(b)}\\&&&\ddots &\ddots &\ddots \\&&&&{\boldsymbol {W}}_{p-1}^{(t)}&{\boldsymbol {I}}_{m}\\&&&&&&{\boldsymbol {I}}_{m}&{\boldsymbol {V}}_{p-1}^{(b)}\\&&&&&&{\boldsymbol {W}}_{p}^{(t)}&{\boldsymbol {I}}_{m}\\&&&&&&&&{\boldsymbol {I}}_{m}\end{bmatrix}}{\begin{bmatrix}{\boldsymbol {X}}_{1}^{(t)}\\{\boldsymbol {X}}_{1}^{(b)}\\{\boldsymbol {X}}_{2}^{(t)}\\{\boldsymbol {X}}_{2}^{(b)}\\\vdots \\{\boldsymbol {X}}_{p-1}^{(t)}\\{\boldsymbol {X}}_{p-1}^{(b)}\\{\boldsymbol {X}}_{p}^{(t)}\\{\boldsymbol {X}}_{p}^{(b)}\end{bmatrix}}={\begin{bmatrix}{\boldsymbol {G}}_{1}^{(t)}\\{\boldsymbol {G}}_{1}^{(b)}\\{\boldsymbol {G}}_{2}^{(t)}\\{\boldsymbol {G}}_{2}^{(b)}\\\vdots \\{\boldsymbol {G}}_{p-1}^{(t)}\\{\boldsymbol {G}}_{p-1}^{(b)}\\{\boldsymbol {G}}_{p}^{(t)}\\{\boldsymbol {G}}_{p}^{(b)}\end{bmatrix}}$

и может быть легко решено параллельно ^[3] .

Усеченный алгоритм SPIKE можно обернуть внутрь какой-нибудь внешней итеративной схемы (например, BiCGSTAB или итеративного уточнения ) для повышения точности решения.

SPIKE для трехдиагональных систем

Первое разделение и алгоритм SPIKE были представлены в ^[4] и был разработан как средство улучшения свойств устойчивости параллельного решателя на основе вращений Гивенса для трехдиагональных систем. Версия алгоритма под названием g-Spike, основанная на последовательном вращении Гивенса, применяемом независимо к каждому блоку, была разработана для графического процессора NVIDIA. ^[5] . Алгоритм на основе SPIKE для графического процессора, основанный на специальной стратегии диагонального поворота блоков, описан в ^[6] .

SPIKE как предобуславливатель

Алгоритм SPIKE также может выступать в качестве предварительного условия для итерационных методов решения линейных систем. Чтобы решить линейную систему Ax = b с использованием итеративного решателя с предварительной обработкой SPIKE, из A извлекают центральные полосы , чтобы сформировать полосовой предобуславливатель M , и решают линейные системы, включающие M на каждой итерации, с помощью алгоритма SPIKE.

Чтобы предобуславливатель был эффективным, обычно необходима перестановка строк и/или столбцов, чтобы переместить «тяжелые» элементы A близко к диагонали, чтобы они были охвачены предобусловливателем. Это может быть достигнуто путем вычисления переупорядочения A взвешенного спектрального .

Алгоритм SPIKE можно обобщить, не ограничивая предобусловливатель строго полосами. В частности, диагональный блок в каждом разделе может быть общей матрицей и, таким образом, обрабатываться прямым решателем общей линейной системы, а не ленточным решателем. Это расширяет возможности предварительного обуславливателя и, следовательно, повышает вероятность сходимости и уменьшает количество итераций.

Реализации

Intel предлагает реализацию алгоритма SPIKE под названием Intel Adaptive Spike-Based Solver. ^[7] . Трехдиагональные решатели также были разработаны для графического процессора NVIDIA. ^{[8] [9]} и сопроцессоры Xeon Phi. Метод в ^[10] является основой трехдиагонального решателя в библиотеке cuSPARSE. Решатель на основе вращений Гивенса был также реализован для Графический процессор и Intel Xeon Phi.

Ссылки

1. ^ Полицци, Э.; Самех, АХ (2006). «Решатель параллельной гибридной полосовой системы: алгоритм SPIKE». Параллельные вычисления . 32 (2): 177–194. дои : 10.1016/j.parco.2005.07.005 .
2. ^ Полицци, Э.; Самех, АХ (2007). «SPIKE: параллельная среда для решения полосовых линейных систем». Компьютеры и жидкости . 36 : 113–141. doi : 10.1016/j.compfluid.2005.07.005 .
3. ^ Миккельсен, CCK; Мангуоглу, М. (2008). «Анализ алгоритма усеченного SPIKE». СИАМ Дж. Матричный анал. Прил. 30 (4): 1500–1519. CiteSeerX   10.1.1.514.8748 . дои : 10.1137/080719571 .
4. ^ Мангуоглу, М.; Самех, АХ; Шенк, О. (2009). «PSPIKE: параллельный гибридный решатель разреженных линейных систем». Параллельная обработка Euro-Par 2009 . Конспекты лекций по информатике. Том. 5704. стр. 797–808. Бибкод : 2009LNCS.5704..797M . дои : 10.1007/978-3-642-03869-3_74 . ISBN  978-3-642-03868-6 .
5. ^ «Адаптивный решатель Intel на основе шипов — сеть программного обеспечения Intel» . Проверено 23 марта 2009 г.
6. ^ Самех, АХ; Кук, диджей (1978). «Об устойчивых решателях параллельных линейных систем» . Журнал АКМ . 25 : 81–91. дои : 10.1145/322047.322054 . S2CID   17109524 .
7. ^ Венетис, IE; Курис, А.; Собчик А.; Галлопулос, Э.; Самех, А.Х. (2015). «Прямой трехдиагональный решатель, основанный на вращении Гивенса для архитектур графических процессоров». Параллельные вычисления . 25 : 101–116. дои : 10.1016/j.parco.2015.03.008 .
8. ^ Чанг, Л.-В.; Страттон, Дж.; Ким, Х.; Хву, В.-М. (2012). «Масштабируемый, численно стабильный, высокопроизводительный трехдиагональный решатель с использованием графических процессоров». Учеб. Международный Конф. Высокопроизводительные вычисления, сетевое хранение и анализ (SC'12) . Лос-Аламитос, Калифорния, США: IEEE Computer Soc. Пресса: 27:1–27:11. ISBN  978-1-4673-0804-5 .

Дальнейшее чтение

• Галлопулос, Э.; Филипп, Б.; Самех, А.Х. (2015). Параллелизм в матричных вычислениях . Спрингер. ISBN  978-94-017-7188-7 .