目的
一般的なブロック三重対角行列の LU 因数分解を行う。
解決方法
インテル® oneMKL LAPACK は、密行列、帯行列、三重対角行列を含む、一般的な行列の LU 因数分解用のさまざまなサブルーチンを提供しています。この手法は、すべてのブロックが正方形で同位という条件を前提として、一般的なブロック三重対角行列の機能の範囲を拡張します。
サイズ NB x NB の正方形ブロックでブロック三重対角行列の LU 因数分解を行うには:
部分 LU 因数分解を、行列の最初の 2 つのブロック行と最初の 3 つのブロック列 (M = 2NB、N = 3NB、K = NB) で構成されたサイズ M x N の長方形ブロックにシーケンシャルに適用し、最後の 1 つのブロック列が処理されるまで対角線に沿って下に移動します。
部分 LU 因数分解: 一般的なブロック三重対角行列の LU 因数分解では、長方形の M x N 行列の部分 LU 因数分解用に別の機能が用意されていると便利です。パラメーター K (ここで、K≤ min(M, N)) の部分 LU 因数分解アルゴリズムは、次のステップからなります。
M x K の部分行列の LU 因数分解を実行します。
三角係数行列を含む方程式を解きます。
右下の (M - K) x (N - K) ブロックを更新します。
生成される行列は A = P(LU + A1) です。ここで、L は下台形 M x K 行列、U は上台形行列、P は置換 (ピボット) 行列、A1 は最後の M - K 行と N - K 列の交差領域のみの非ゼロ要素行列です。
一般的な LU 因数分解を最後の (2NB) x (2NB) ブロックに適用します。
ソースコード: サンプル (https://www.intel.com/content/dam/develop/external/us/en/documents/mkl-cookbook-samples-120115.zip (英語)) の BlockTDS_GE/source/dgeblttrf.f ファイルを参照してください。GitHub の block_lu_decomposition (英語) も参照してください。
部分 LU 因数分解の実行
SUBROUTINE PTLDGETRF(M, N, K, A, LDA, IPIV, INFO)
…
CALL DGETRF( M, K, A, LDA, IPIV, INFO )
…
DO I=1,K
IF(IPIV(I).NE.I)THEN
CALL DSWAP(N-K, A(I,K+1), LDA, A(IPIV(I), K+1), LDA)
END IF
END DO
CALL DTRSM('L','L','N','U',K,N-K,1D0, A, LDA, A(1,K+1), LDA)
CALL DGEMM('N', 'N', M-K, N-K, K, -1D0, A(K+1,1), LDA,
& A(1,K+1), LDA, 1D0, A(K+1,K+1), LDA)
…ブロック三角行列の因数分解
…
DO K=1,N-2
C ブロック構造の2*NB x 3*NB の部分行列A を構成
C (D_K C_K 0 )
C (B_K D_K+1 C_K+1)
…
C 部分行列を部分的に因数分解
CALL PTLDGETRF(2*NB, 3*NB, NB, A, 2*NB, IPIV(1,K), INFO)
C 因数分解の結果を三重対角行列のブロックを格納する配列に戻す
…
END DO
C ループを抜けて最後の2*NB x 2*NB の部分行列を因数分解
CALL DGETRF(2*NB, 2*NB, A, 2*NB, IPIV(1,N-1), INFO)
C 最後の結果を三重対角行列のブロックを格納する配列に戻す
…使用されるルーチン
タスク |
ルーチン |
説明 |
|---|---|---|
M x K の部分行列の LU 因数分解 |
DGETRF |
一般的な m x n 行列の LU 因数分解を計算します。 |
行列の行の置換 |
DSWAP |
2 つのベクトルをスワップします。 |
三角係数行列を含む方程式を解く |
DTRSM |
三角行列を含む方程式を解きます。 |
右下の (M - K) x (N - K) ブロックを更新する |
DGEMM |
一般的な行列を含む行列-行列の積を計算します。 |
説明
部分 LU 因数分解では、A を長方形の m x n 行列にします。
注
読みやすいように、ここでは m、n、k、nb のように小文字のインデックスを使用しています。これらのインデックスは、Fortran ソリューションおよびコードサンプルで使用されている大文字のインデックスに相当します。
行列は、m x k (ここで、0 k≤n) の部分行列の LU 因数分解を使用して分解できます。このアプリケーションでは、m≤k≤n または n = k≤m の場合は行列の因数分解に ?getrf を直接使用できるため、k< min(m, n) です。
A はブロック行列として表現できます。
ここで、A11 は k x k の部分行列、A12 は k x (n - k) の部分行列、A21 は (m - k) x k の部分行列、A22 は (m - k) x (n - k) の部分行列です。
m x k のパネル A1 は次のように定義できます。
A1 は (?getrf を使用して) A1 = PLU として LU 因数分解できます。ここで、P は置換 (ピボット) 行列、L は対角要素を含む下台形行列、U は上三角行列です。
注
L の対角要素は格納する必要がないため、A1 を格納するために使用する配列を、L と U の要素を格納するために使用できます。
PT を A の 2 つ目のパネルに適用すると次のようになります。
次の方程式が得られます。
A''12 を次のように定義します。
PTA の方程式を変更します。
前の方程式と P を乗算すると次のようになります。
この方程式は最初の行列の部分 LU 因数分解と考えることができます。
注
積 L-111A'12 は ?trsm を呼び出して計算し、A12 に使用される配列の代わりに格納できます。更新 A'22 - L21(L-111A'12) は ?gemm を呼び出して計算し、A22 に使用される配列の代わりに格納できます。
部分行列にすべてのランクが含まれない場合、LU 因数分解に失敗するため、このメソッドは適用できません。
一般的な行列の LU 因数分解とは異なり、下記に示す一般的なブロック三重対角行列の因数分解 A = LU は、A = PLU (P は置換行列) 形式で記述できません。ピボットのため、左の係数 L は置換を含みます。また、右の係数 U も複雑になります (2 つの対角線ではなく 3 つの対角線を含む)。
ロック三重対角行列の LU 因数分解では、A を、すべてのブロックが正方形で同位のブロック三重対角行列 nb にします。
行列は A = LU として因数分解されます。
最初に、2 x 3 ブロックの部分行列を考えます。
この部分行列は次のように分解できます。
この分解は前述の部分 LU 因数分解を適用して取得できます。ここで、PT1 は nb 要素の順列の積であり、2nb x nb 行列として表現できます。
すべてのブロックがサイズ nb x nb の N x N ブロック行列を適用します。
分解は次のように表現できます。
次に、方程式の右辺の行列の、2 番目と 3 番目の行の 2 x 3 ブロック行列を因数分解します。
ここで、PT2 は次のように定義されます。
分解は次のようになります。
この表記をピボット行列に適用して方程式を単純化します。
ここで、PTj は 2nb x 2nb で j 番目と (j+1) 番目の行と列の交差領域にあります。分解は次のように簡潔に表記されます。
ステップ N - 2 のローカル因数分解は次のようになります。
このステップの後、ピボット行列で乗算します。
分解は次のようになります。
最後の (N - 1) 番目のステップでは、行列は正方形で因数分解は完了です。
最後のステップはそれまでのステップとピボットの構造が異なります。前の PTj (j = 1, 2, ..., N - 2) はすべて nb 置換の積 (nb 整数パラメーターに依存) でしたが、PTN-1 は次数 2nb の正方行列 (2nb パラメーターに依存) に適用されます。そのため、ピボット・インデックスをすべて格納するには、長さ (N - 2) nb + 2nb = Nnb の整数配列が必要です。
左から前の分解と ΠTN-1 を乗算すると、最終的な分解が得られます。
この分解と Π1Π2...ΠN-1 を乗算すると、LU 因数分解形式で記述できます。
この式を適用するときは、Πj (j = 1, 2, …, N-2) はインデックス j と j+1 のブロック行に適用される nb 転置の積ですが、ΠN-1 は最後の 2 つのブロック行 N-1 と N に適用される 2nb 転置の積であることに注意してください。