ブロック三重対称正定値係数行列を含む連立線形方程式を解く

目的

コレスキー因数分解された対称正定値ブロック三重対角係数行列を含む連立線形方程式を解く。

解決方法

係数対称正定値ブロック三重対角行列 (それぞれ同じサイズNB x NB の正方形ブロック) を LLT 因数分解する場合、次の 2 段階で解きます。

N x N ブロック (サイズ NB x NB) で、対角ブロックが下三角行列の下二重対角係数行列を含む連立線形方程式を解きます:
1. 右辺ベクトルの係数行列として使用されるサイズ NB x NB の下三角対角ブロックを含む N 連立方程式を解きます。
2. 右辺を更新します。
N x N ブロック (サイズ NB x NB) で、対角ブロックが上三角行列の上二重対角係数行列を含む連立線形方程式を解きます。
1. 連立方程式を解きます。
2. 右辺を更新します。

ソースコード: サンプル (https://www.intel.com/content/dam/develop/external/us/en/documents/mkl-cookbook-samples-120115.zip (英語)) の BlockTDS_SPD/source/dpbltrf.f ファイルを参照してください。GitHub の block_cholesky_decomposition (英語) も参照してください。

対称正定値ブロック三重対角行列のコレスキー因数分解

   … 
    … 
    CALL DTRSM('L', 'L', 'N', 'N', NB, NRHS, 1D0, D, LDD, F, LDF) 
    DO K = 2, N 
        CALL DGEMM('N', 'N', NB, NRHS, NB, -1D0, B(1,(K-2)*NB+1), LDB, F((K-2)*NB+1,1), LDF, 1D0, F((K-1)*NB+1,1), LDF) 
        CALL DTRSM('L','L', 'N', 'N', NB, NRHS, 1D0, D(1,(K-1)*NB+1), LDD, F((K-1)*NB+1,1), LDF) 
    END DO 

    CALL DTRSM('L', 'L', 'T', 'N', NB, NRHS, 1D0, D(1,(N-1)*NB+1), LDD, F((N-1)*NB+1,1), LDF) 
    DO K = N-1, 1, -1 
        CALL DGEMM('T', 'N', NB, NRHS, NB, -1D0, B(1,(K-1)*NB+1), LDB, F(K*NB+1,1), LDF, 1D0, F((K-1)*NB+1,1), LDF) 
        CALL DTRSM('L','L', 'T', 'N', NB, NRHS, 1D0, D(1,(K-1)*NB+1), LDD, F((K-1)*NB+1,1), LDF) 
    END DO 
    …

使用されるルーチン

タスク	ルーチン	説明
連立線形方程式を解く	DTRSM	三角行列を含む方程式を解きます。
右辺を更新する	DGEMM	一般的な行列を含む行列-行列の積を計算します。

説明

次の連立線形方程式について考えます。

行列 A は、すべてのブロックがサイズ NB x NB の対称正定値ブロック三重対角係数行列であると仮定します。A は、ブロック三重対称正定値行列の因数分解で説明しているように、次のように因数分解できます。

連立方程式を解くアルゴリズムは次のとおりです。

対角ブロックが下三角行列の下二重対角係数行列を含む連立線形方程式を解きます。
```
Y₁=L₁^-1 F₁ // trsm()
do i=2,N
     G_i=F_i - C_{i - 1} Y_{i - 1} // gemm()
     Y_i=L_i^-1 G_i // trsm()
end do
```
対角ブロックが上三角行列の上二重対角係数行列を含む連立線形方程式を解きます。
```
X_N=L_N^-T Y_N // trsm()
do i=N-1,1,-1
     Z_i=F_i-C_i^T X_{i + 1} // gemm()
     X_i=L_i^-T Z_i // trsm()
end do
```