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

目的

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

ソリューション

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

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

前方置換。ピボットで連立方程式 LY=F (L は下三角係数行列) を解きます。因数分解されたブロック三重対角行列では、最後のブロックを除く Y のブロックはすべて、次の方法によりループ内で見つかります。
1. 右辺にピボット置換を適用します。
2. 下三角係数行列の NB 線形方程式を解きます (NB はブロックの次数)。
3. 次のステップのために右辺を更新します。
最後のピボットの構造 (2 つのブロック置換を連続して適用する必要がある) と係数行列の構造により、最後の 2 つのブロック・コンポーネントはループの外で見つかります。
後方置換。式 UX=Y を解きます。ピボットを含まないため、このステップはより単純です。プロシージャーは、最初のステップに似ています。
1. 三角係数行列を含む方程式を解きます。
2. 右辺のブロックを更新します。
前のステップとの違いは、係数行列が下三角ではなく上三角で、ループの向きが逆になっていることです。

ソースコード: サンプル (http://software.intel.com/en-us/mkl_cookbook_samples) の BlockTDS_GE/source/dgeblttrs.f ファイルを参照してください。

前方置換

! 前方置換
! ループ内で配列 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 の帯行列として扱い、インテル® MKL LAPACK の帯行列を因数分解するサブルーチン (?gbtrf) と、帯行列を解くサブルーチン (?gbtrs) を呼び出して解くことができます。しかし、ブロック行列を帯行列として格納すると、帯の多くのゼロ要素が非ゼロとして扱われて計算中に処理されるため、この手法で説明したアプローチで必要な浮動小数点計算は少なくなります。より大きな NB では影響も大きくなります。帯行列をブロック三重対角行列として扱うこともできますが、ブロックに非ゼロとして扱われる多くのゼロが含まれるため、この格納手法は効率的ではありません。このため、帯格納手法とブロック三重対角格納手法、およびそれらのソルバーは、互いに補完的な手法として考えるべきです。

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