■LU Decompositionのプログラミング: 正方行列の場合■
有限要素法のソフトで生成される剛性マトリクスは、メモリーを節約するためにバンドマトリクス化されています。 それに伴い三角行列分解のアルゴリズムにも多少の手を入れる必要があります。 まず、ここでは、正方マトリクスを三角行列分解を行う方法を紹介します。後ほどバンドマトリクスの三角行列分解 を紹介します。

正方マトリクスの場合、プログラムにすると三角行列分解は以下のサブプログラムで計算できます。 Ｗｅｂサイトに掲載されているCholesky分解法の解説書をベースにプログラムしてみました。

      SUBROUTINE DECOMP ( MXN, NNODE, A )
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION A(MXN,MXN)
C------- COMPUTATION OF L(I,J)
C-------I = 1
      DO J = 2, NNODE
         A(1,J) = A(1,J)/A(1,1)
      END DO
C------ I >= 2
      DO I = 2 , NNODE
         SUM = 0.D0
            DO M = 1, I-1
              SUM = SUM + A(M,I)**2*A(M,M)
            END DO
         A(I,I) = A(I,I) - SUM
      DO J = I+1, NNODE
         SUM = 0.D0
            DO M = 1, I-1
              SUM = SUM + A(M,I)*A(M,J)*A(M,M)
            END DO
            A(I,J) = (A(I,J)-SUM)/A(I,I)
      END DO
      END DO
C ----------- PUT ZERO INTO LOWER [A]
      DO I = 1 , NNODE
      DO J = 1, I-1
      A(I,J) = 0.D0
      END DO
      END DO
      RETURN
      END

上のサブプログラムを用いて連立方程式を解く方法を紹介します。まず、[A]{x}={B}の[A]と{B}を読み込みます。 次に、[L]と[D]を次の行で計算します。これが、連立方程式の計算を行う第1ステップになります。

      CALL DECOMP ( MXN, NNODE, A )

すると、[A]の上半分には[L]^Tが、対角には[D]が入ります。したがって、 計算しなければならない連立方程式は[L][D][L]^T{x}={B}でしたね。 そして、[D][L]^T{x}={x'}と置き、[L][D][L]^T{x}={B} は、[L]{x'}={B}で書けました。 {x'}は代入法で得られるということでしたよね。 下がそのプログラムです。そして、{x'}の計算結果は、{B}に入ります。

      DO I = 2 , NNODE
      SUM = 0.D0
      DO J = 1, I-1
      SUM = SUM + A(I,J)*B(J)
      END DO
      B(I) = B(I) - SUM
      END DO

      DO I = 1 , NNODE
      B(I) = B(I) / A(I,I)
      END DO

      DO I = NNODE-1, 1, -1
      SUM = 0.D0
      DO J = I+1 , NNODE
      SUM = SUM + A(I,J)*B(J)
      END DO
      B(I) = B(I) - SUM
      END DO