SPIKEアルゴリズム

SPIKEアルゴリズムは、エリック・ポリッツィとアハメド・サメフによって開発されたバンド線形システム用のハイブリッド並列ソルバーである[1] ^ [2]

概要

SPIKEアルゴリズムは線形システムAX = Fを扱います。ここで、Aはよりはるかに小さい帯域幅の帯行列であり、Fは右辺を 含む行列です。このアルゴリズムは前処理段階と後処理段階に分かれています。n×n{\displaystyle n\times n}n{\displaystyle n}n×s{\displaystyle n\times s}s{\displaystyle s}

前処理段階

前処理段階では、線形システムAX = Fはブロック三角形式 に分割される。

[1B1C22B2Cp1p1Bp1Cpp][X1X2Xp1Xp][F1F2Fp1Fp]{\displaystyle {\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 jj = 1,..., p 、ただしp ≥ 2)が非特異であると仮定する。ブロック対角行列 を定義する。

D = diag( A 1 ,..., A p ),

するとDも非特異となる。この系の両辺に D −1を左乗すると

[V1W2V2Wp1Vp1Wp][X1X2Xp1Xp][G1G2Gp1Gp]{\displaystyle {\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 −1による左乗算は、以下の形式の系を 解くことと等価であるp{\displaystyle p}

A j [ V j W j G j ] = [ B j C j F j ]

(のW 1C 1、 のV pB pを省略)、これらは並列に実行できます。 j=1{\displaystyle j=1}j=p{\displaystyle j=p}

Aの帯状性により、各V jの左端の列と各W jの右端の列のうち、わずかしか非ゼロになりません。これらの列はスパイクと呼ばれます。

後処理段階

一般性を損なうことなく、各スパイクにはちょうど列(は よりはるかに小さい)が含まれると仮定する(必要であればスパイクをゼロの列で埋める)。すべてのV jW jのスパイクを次のように 分割する。m{\displaystyle m}m{\displaystyle m}n{\displaystyle n}

[Vj(t)VjVj(b)]{\displaystyle {\begin{bmatrix}{\boldsymbol {V}}_{j}^{(t)}\\{\boldsymbol {V}}_{j}'\\{\boldsymbol {V}}_{j}^{(b)}\end{bmatrix}}}そして[Wj(t)WjWj(b)]{\displaystyle {\begin{bmatrix}{\boldsymbol {W}}_{j}^{(t)}\\{\boldsymbol {W}}_{j}'\\{\boldsymbol {W}}_{j}^{(b)}\\\end{bmatrix}}}

ここで、V tj V bj W tj そしてW bj 次元である。同様にすべてのX jG jを 分割する。m×m{\displaystyle m\times m}

[Xj(t)XjXj(b)]{\displaystyle {\begin{bmatrix}{\boldsymbol {X}}_{j}^{(t)}\\{\boldsymbol {X}}_{j}'\\{\boldsymbol {X}}_{j}^{(b)}\end{bmatrix}}}そして[Gj(t)GjGj(b)].{\displaystyle {\begin{bmatrix}{\boldsymbol {G}}_{j}^{(t)}\\{\boldsymbol {G}}_{j}'\\{\boldsymbol {G}}_{j}^{(b)}\\\end{bmatrix}}.}

前処理段階で生成されたシステムは、はるかに小さいサイズのブロック五角形システムに縮小できることに注目してください(は よりもはるかに小さいことを思い出してください)。 m{\displaystyle m}n{\displaystyle n}

[Im0V1(t)0ImV1(b)00W2(t)Im0V2(t)W2(b)0ImV2(b)00Wp1(t)Im0Vp1(t)Wp1(b)0ImVp1(b)00Wp(t)Im0Wp(b)0Im][X1(t)X1(b)X2(t)X2(b)Xp1(t)Xp1(b)Xp(t)Xp(b)]=[G1(t)G1(b)G2(t)G2(b)Gp1(t)Gp1(b)Gp(t)Gp(b)],{\displaystyle {\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̃ = と表記します。

一度すべてのX tj およびX bj が見つかったら 、すべてのX′j

{X1=G1V1X2(t),Xj=GjVjXj+1(t)WjXj1(b),j=2,,p1,Xp=GpWpXp1(b).{\displaystyle {\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

論理的には 2 つの段階に分かれていますが、計算上は SPIKE アルゴリズムは次の 3 つの段階で構成されます。

  1. 対角ブロックを因数分解し、
  2. スパイクを計算し、
  3. 縮小されたシステムを解く。

これらの各段階は複数の方法で実行でき、多様なバリエーションが考えられます。注目すべきバリエーションとして、非対角優勢の場合の再帰SPIKEアルゴリズムと、対角優勢の場合の打ち切りSPIKEアルゴリズムがあります。これらのバリエーションに応じて、系は正確に解くことも近似的に解くこともできます。後者の場合、SPIKEはクリロフ部分空間法反復改良法などの反復法の前処理として用いられます。

再帰SPIKE

前処理段階

前処理段階の最初のステップは、対角ブロックA jを因数分解することです。数値安定性を確保するため、 LAPACKXGBTRFルーチンを用いて、部分ピボットを用いたLU分解を行うことができます。あるいは、部分ピボットを用いずに「対角ブースティング」戦略を用いて分解することも可能です。後者の方法は、特異な対角ブロックの問題に対処します。

具体的には、対角ブースティング戦略は以下の通りである。0 ε設定可能な「マシンゼロ」と表記する。LU分解の各ステップにおいて、ピボットが以下の条件を満たすことを要求する。

|ピボット| > 0 εA1

ピボットが条件を満たさない場合は、

pivot={pivot+ϵAj1if pivot0,pivotϵAj1if pivot<0{\displaystyle \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。対角ブロックが因数分解された後、スパイクが計算され、後処理段階に渡されます。

後処理段階

2つの分割の場合

2分割の場合、すなわちp = 2のとき、簡約されたシステムS̃X̃ = は次の形をとる。

[Im0V1(t)0ImV1(b)00W2(t)Im0W2(b)0Im][X1(t)X1(b)X2(t)X2(b)]=[G1(t)G1(b)G2(t)G2(b)].{\displaystyle {\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{.}}}

中心からさらに小さなシステムを抽出することができます。

[ImV1(b)W2(t)Im][X1(b)X2(t)]=[G1(b)G2(t)],{\displaystyle {\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分解を使って解くことができる。

[ImV1(b)W2(t)Im]=[ImW2(t)Im][ImV1(b)ImW2(t)V1(b)].{\displaystyle {\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{.}}}

一度X )1 およびX トン)2 見つかった場合、X トン)1 およびX b)2 計算は次のように行えます

X トン)1 = G トン)1 V トン)1 X トン)2 
X b)2 = G b)2 西 b)2 X )1 
複数パーティションの場合

pは2のべき乗、すなわちp = 2 dであると仮定する。ブロック対角行列を考える。

1 = 診断( [1] 1 、...、 [1] p /2 

どこ

D~k[1]=[Im0V2k1(t)0ImV2k1(b)00W2k(t)Im0W2k(b)0Im]{\displaystyle {\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}}}

k = 1,..., p /2に対して、 1 は本質的にから抽出された4 mの対角ブロックから構成されていることに注意してください。ここでS̃ を因数分解すると、

= 1 2

新しい行列S̃2は次の形をとる

[I3m0V1[2](t)0ImV1[2](b)00W2[2](t)Im0V2[2](t)W2[2](b)0I3mV2[2](b)00Wp/21[2](t)I3m0Vp/21[2](t)Wp/21[2](b)0ImVp/21[2](b)00Wp/2[2](t)Im0Wp/2[2](b)0I3m].{\displaystyle {\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̃ 2と非常によく似ており、スパイクの数と高さ(幅はmで一定)が異なるだけです。したがって、 2に対して 同様の因数分解を行うと、

2 = 2 3

そして

= 1 2 3

このような因数分解のステップは再帰的に実行することができる。d − 1ステップ実行後、次の因数分解が得られる

= 1 d −1 d

ここでS̃d2つのスパイクのみを持つ。縮約された系は次のように解かれる。

= −1  −1 d −1  −1 1 

2 分割の場合のブロック LU 分解手法は、 1、...、d −1、およびdを含む解決ステップを処理するために使用できます。これは、基本的に、一般化された 2 分割形式の複数の独立したシステムを解くためです。

pが 2 の累乗ではない 場合への一般化はほぼ簡単です。

切り捨てられたSPIKE

Aが対角優位の場合、簡約システムでは

[Im0V1(t)0ImV1(b)00W2(t)Im0V2(t)W2(b)0ImV2(b)00Wp1(t)Im0Vp1(t)Wp1(b)0ImVp1(b)00Wp(t)Im0Wp(b)0Im][X1(t)X1(b)X2(t)X2(b)Xp1(t)Xp1(b)Xp(t)Xp(b)]=[G1(t)G1(b)G2(t)G2(b)Gp1(t)Gp1(b)Gp(t)Gp(b)],{\displaystyle {\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 tj そしてW bj は無視できることが多い。これらを省略すると、簡約されたシステムはブロック対角になる。

[ImImV1(b)W2(t)ImImV2(b)Wp1(t)ImImVp1(b)Wp(t)ImIm][X1(t)X1(b)X2(t)X2(b)Xp1(t)Xp1(b)Xp(t)Xp(b)]=[G1(t)G1(b)G2(t)G2(b)Gp1(t)Gp1(b)Gp(t)Gp(b)]{\displaystyle {\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 GPU用に設計されました[5] 。特殊なブロック対角ピボット戦略に基づくGPU向けのSPIKEベースアルゴリズムは、 [6]で説明されています。

前処理としてのSPIKE

SPIKEアルゴリズムは、線形システムを解く反復法の前処理としても機能します。SPIKE前処理付き反復ソルバーを用いて線形システムAx = bを解くには、 Aから中心バンドを抽出してバンド前処理Mを作成し、各反復においてMを含む線形システムをSPIKEアルゴリズムで解きます。

前処理子を効果的に機能させるには、通常、行および/または列の置換によってAの「重い」要素を対角線に近づけ、前処理子の対象となるようにする必要があります。これは、 A重み付きスペクトル並べ替えを計算することで実現できます。

SPIKEアルゴリズムは、前処理を厳密に帯状化することを制限しないことで一般化できます。特に、各パーティションの対角ブロックは一般行列とすることができるため、帯状化ソルバーではなく、直接的な一般線形システムソルバーで処理できます。これにより前処理が強化され、収束の可能性が高まり、反復回数が削減されます。

実装

Intelは、SPIKEアルゴリズムの実装をIntel Adaptive Spike-Based Solver [7]という名称で提供しています。三角関数ソルバーは、NVIDIA GPU [8] [9]およびXeon Phiコプロセッサ向けにも開発されています。 [10]の手法は、 cuSPARSEライブラリの三角関数ソルバーの基礎となっています。[ 1 ]ギブンズ回転に基づくソルバーもGPUおよびIntel Xeon Phi向けに実装されています。[ 2 ]

参考文献

  1. ^ NVIDIA、2014 年 10 月 28 日にアクセス。CUDA ツールキット ドキュメント v. 6.5: cuSPARSE、 http://docs.nvidia.com/cuda/cusparse
  2. ^ベネティス、イオアニス;ソプチク、アレクサンドロス。クーリス、アレクサンドロス。ナコス、アレクサンドロス。ニコラオス、ニコラウサコス。ガロプロス、エフストラティオス (2015-09-03)。「コプロセッサー用の一般的な三重対角ソルバー: Intel Xeon Phi 向けの g-Spike の適応」 – ResearchGateより。
  1. ^ポリッツィ、E.;サメ、AH (2006)。 「並列ハイブリッド帯域システム ソルバー: SPIKE アルゴリズム」。並列コンピューティング 32 (2): 177–194土井: 10.1016/j.parco.2005.07.005
  2. ^ Polizzi, E.; Sameh, AH (2007). 「SPIKE: 帯状線形システムを解くための並列環境」. Computers & Fluids . 36 : 113– 141. doi : 10.1016/j.compfluid.2005.07.005 .
  3. ^ Mikkelsen, CCK; Manguoglu, M. (2008). 「Truncated SPIKEアルゴリズムの解析」. SIAM J. Matrix Anal. Appl. 30 (4): 1500– 1519. CiteSeerX 10.1.1.514.8748 . doi : 10.1137/080719571 . 
  4. ^ Manguoglu, M.; Sameh, AH; Schenk, O. (2009). 「PSPIKE: 並列ハイブリッドスパース線形システムソルバー」. Euro-Par 2009 Parallel Processing . Lecture Notes in Computer Science. Vol. 5704. pp.  797– 808. Bibcode : 2009LNCS.5704..797M . doi : 10.1007/978-3-642-03869-3_74 . ISBN 978-3-642-03868-6
  5. ^「Intel Adaptive Spike-Based Solver - Intel Software Network」 。 2009年3月23日閲覧
  6. ^ Sameh, AH; Kuck, DJ (1978).「安定な並列線形システムソルバーについて」 . Journal of the ACM . 25 : 81–91 . doi : 10.1145/322047.322054 . S2CID 17109524 . 
  7. ^ベネティス、アイオワ州;クーリス、A.ソプチク、A.ガロプロス、E.サメ、AH (2015)。 「GPU アーキテクチャ用のギブンス回転に基づく直接三重対角ソルバー」。並列コンピューティング 25 : 101–116土井: 10.1016/j.parco.2015.03.008
  8. ^ Chang, L.-W.; Stratton, J.; Kim, H.; Hwu, W.-M. (2012). 「GPUを用いたスケーラブルで数値的に安定した高性能三角関数ソルバー」 Proc. Int'l. Conf. High Performance Computing, Networking Storage and Analysis (SC'12) . Los Alamitos, CA, USA: IEEE Computer Soc. Press: 27:1–27:11. ISBN 978-1-4673-0804-5

さらに読む