< 目次

LU 因数分解されたブロック三重対角係数行列を含む連立線形方程式を解く

目的

oneMKL LAPACK のルーチンを使用してブロック三重対角係数行列を含む連立方程式の解を求める (LAPACK にはブロック三重対角係数行列を含む式を直接解くルーチンがないため)。

解決方法

oneMKL LAPACK は、LU 因数分解された係数行列を含む連立方程式を解くためのさまざまなサブルーチンを提供しています。(密行列、帯行列、三重対角行列など)。この手法は、すべてのブロックが正方形で同位という条件を前提として、このセットをブロック三重対角行列に拡張します。ブロック三重対角行列 A の形式は次のとおりです。




LU 因数分解された行列 A=LU と複数の右辺 (RHS) を含む式 AX=F は、2 段階で解きます (LU 因数分解の詳細は、一般的なブロック三重対角行列の因数分解を参照)。

  1. 前方置換。ピボットで連立方程式 LY=F (L は下三角係数行列) を解きます。因数分解されたブロック三重対角行列では、最後のブロックを除く Y のブロックはすべて、次の方法によりループ内で見つかります。

    1. 右辺にピボット置換を適用します。

    2. 下三角係数行列の NB 線形方程式を解きます (NB はブロックの次数)。

    3. 次のステップのために右辺を更新します。

    最後のピボットの構造 (2 つのブロック置換を連続して適用する必要がある) と係数行列の構造により、最後の 2 つのブロック・コンポーネントはループの外で見つかります。

  2. 後方置換。式 UX=Y を解きます。ピボットを含まないため、このステップはより単純です。プロシージャーは、最初のステップに似ています。

    1. 三角係数行列を含む方程式を解きます。

    2. 右辺のブロックを更新します。

    前のステップとの違いは、係数行列が下三角ではなく上三角で、ループの向きが逆になっていることです。

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

前方置換

! 前方置換
! ループ内で配列F に格納されているコンポーネントY_K を計算 
    DO K = 1, N-2 
        DO I = 1, NB IF (IPIV(I,K) .NE. I) THEN 
            CALL DSWAP(NRHS, F((K-1)*NB+I,1), LDF, F((K-1)*NB+IPIV(I,K),1), LDF) 
        END IF 
    END DO 
    CALL DTRSM('L', 'L', 'N', 'U', NB, NRHS, 1D0, D(1,(K-1)*NB+1), NB, F((K-1)*NB+1,1), LDF) 
    CALL DGEMM('N', 'N', NB, NRHS, NB, -1D0, DL(1,(K-1)*NB+1), NB, F((K-1)*NB+1,1), LDF, 1D0, + F(K*NB+1,1), LDF) 
    END DO 

! 最後の2 つのピボットを適用 
    DO I = 1, NB 
        IF (IPIV(I,N-1) .NE. I) THEN 
            CALL DSWAP(NRHS, F((N-2)*NB+I,1), LDF, F((N-2)*NB+IPIV(I,N-1),1), LDF) 
        END IF 
    END DO 

    DO I = 1, NB 
        IF(IPIV(I,N)-NB.NE.I) THEN 
            CALL DSWAP(NRHS, F((N-1)*NB+I,1), LDF, F((N-2)*NB+IPIV(I,N),1), LDF) 
        END IF 
    END DO 
! ループの外でY_N-1 とY_N を計算して配列F に格納 
    CALL DTRSM('L', 'L', 'N', 'U', NB, NRHS, 1D0, D(1,(N-2)*NB+1), NB, F((N-2)*NB+1,1), LDF) 
    CALL DGEMM('N', 'N', NB, NRHS, NB, -1D0, DL(1,(N-2)*NB +1), NB, F((N-2)*NB+1,1), LDF, 1D0, + F((N-1)*NB+1,1), LDF)

後方置換

… 
! ループの外でX_N を計算して配列F に格納 
    CALL DTRSM('L', 'U', 'N', 'N', NB, NRHS, 1D0, D(1,(N-1)*NB+1), NB, F((N-1)*NB+1,1), LDF) 
! ループの外でX_N-1 を計算して配列F に格納 
    CALL DGEMM('N', 'N', NB, NRHS, NB, -1D0, DU1(1,(N-2)*NB +1), NB, F((N-1)*NB+1,1), LDF, 1D0, + F((N-2)*NB+1,1), LDF) 
    CALL DTRSM('L', 'U', 'N', 'N', NB, NRHS, 1D0, D(1,(N-2)*NB+1), NB, F((N-2)*NB+1,1), LDF) 
! ループ内で配列F に格納されているコンポーネントX_K を計算 
    DO K = N-2, 1, -1 
        CALL DGEMM('N','N',NB, NRHS, NB, -1D0, DU1(1,(K-1)*NB +1), NB, F(K*NB+1,1), LDF, 1D0, + F((K-1)*NB+1,1), LDF) 
        CALL DGEMM('N','N',NB, NRHS, NB, -1D0, DU2(1,(K-1)*NB +1), NB, F((K+1)*NB+1,1), LDF, 1D0, + F((K-1)*NB+1,1), LDF) 
        CALL DTRSM('L', 'U', 'N', 'N', NB, NRHS, 1D0, D(1,(K-1)*NB+1), NB, F((K-1)*NB+1,1), LDF) 
    END DO 
…

使用されるルーチン

タスク

ルーチン

説明

ピボット置換を適用する

dswap

ベクトルを別のベクトルとスワップします。

下三角係数行列と上三角係数行列を用いて連立線形方程式を解く

dtrsm

三角行列を含む方程式を解きます。

右辺のブロックを更新する

dgemm

一般的な行列を含む行列-行列の積を計算します。

説明

注

サイズ NB x NB のブロックの一般的なブロック三重対角行列は、帯域幅 4*NB-1 の帯行列として扱い、oneMKL LAPACK の帯行列を因数分解するサブルーチン (?gbtrf) と、帯行列を解くサブルーチン (?gbtrs) を呼び出して解くことができます。しかし、ブロック行列を帯行列として格納すると、帯の多くのゼロ要素が非ゼロとして扱われて計算中に処理されるため、この手法で説明したアプローチで必要な浮動小数点計算は少なくなります。より大きな NB では影響も大きくなります。帯行列をブロック三重対角行列として扱うこともできますが、ブロックに非ゼロとして扱われる多くのゼロが含まれるため、この格納手法は効率的ではありません。このため、帯格納手法とブロック三重対角格納手法、およびそれらのソルバーは、互いに補完的な手法として考えるべきです。

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




ブロック三重対角係数行列 A は次のように因数分解されると仮定します。




使用している用語の定義は、一般的なブロック三重対角行列の因数分解を参照してください。

式は 2 つの連立線形方程式に分解されます。




2 つ目の方程式を拡張します。




Y1 を見つけるには、最初に置換 Π1T を適用する必要があります。この置換は、右辺の最初の 2 つのブロックのみ変更します。




置換をローカルに適用します。



Y1 が見つかりました。



Y1 を見つけた後、同様の計算を繰り返してほかの値 (Y2, Y3, ..., YN - 2) を見つけます。

注

ΠN - 1 の異なる構造 (一般的なブロック三重対角行列の因数分解を参照) は、同じ方程式を YN - 1 と YN の計算に使用できず、ループの外で計算する必要があることを意味します。

Y を見つける方程式を用いるアルゴリズムは次のとおりです。

do for k = 1 to N - 2 
                                   
     // ?trsm を使用                         // ?trsm 
     // ?gemm を使用して右側を更新    // ?gemm を使用して右側を更新
end do                                                // ?trsm

UX = Y 方程式は次のように表現できます。



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

 // ?trsm 
 // ?gemm に続く ?trsm 
do for k = N - 1 to 1 in steps of -1 
     // ?gemm 
     // ?gemm 
     // ?trsm end do
上位トピック: インテル® oneAPI マス・カーネル・ライブラリー・クックブック