LU分解FORTRANプログラムNO.3

    2003/08/18 日立製作所 & 早稲田大学 後 保範 ( Ushiro Yasunori )
--------------------------------------------------------------

1. 概要

 実密行列を係数とする連立一次方程式Ax=bをガウス消去法でLU分解し、
LUx=bの数値解を前進後退代入で求める。
部分軸交換付きで軸交換対象行は全て交換するプログラム。

2. プログラム

C=================================================================C
      SUBROUTINE GLU3(A,B,N,ND,EPS,IP,IER)
C=================================================================C
C  Real-Dense LU Decomposition by Gauss Elimination               C
C   and Solve Ax=b by Substitution                                C
C     A ---> LU Decomposition with Partial Pivoting               C
C            Type-1 Partial Pivoting (Changing all rows)          C
C-----------------------------------------------------------------C
C    A(ND,N)  R*8, I/O, A Coefficient Matrix                      C
C    B(N)     R*8, I/O, A right-hand vector(b) and solution(x)    C
C    N        I*4, In,  Matrix Size of A                          C
C    ND       I*4, In,  Array Size of A ( ND >= N )               C
C    EPS      R*8, In,  Value for Singularity  Check              C
C    IP(N)    I*4, Out, Pivot Number                              C
C    IER      I*4, Out, 0 : Normal Execution                      C
C                       1 : Singular Stop                         C
C                       2 ; Parameter Error                       C
C-----------------------------------------------------------------C
C    Written by Yasunori Ushiro,   2003/08/17                     C
C        ( Hitachi Ltd. and Waseda University )                   C
C=================================================================C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION A(ND,N), B(N), IP(N)
C----- Gauss Elimination Step ------------      
C  Check Parameter
      IER = 0
      if(ND.lt.N) then
        IER = 2
        go to 100
      end if
C  Gauss Elimination
      do k=1,N
C   Search Maximum Value in k's column
        KPIV = k
        PIV  = abs(A(k,k))          
        do i=k+1,N
          if(abs(A(i,k)).gt.PIV) then
             KPIV = i
             PIV  = abs(A(i,k))
          end if
        end do
C   Check Singularity
        if(PIV.lt.EPS) then
          IER = 1
          go to 100
        end if
        IP(k) = KPIV
C   Change A(k,*) <--> A(KPIV,*), All rows
        if(KPIV.ne.k) then 
          do j=1,N
            AKJ       = A(k,j)
            A(k,j)    = A(KPIV,j)
            A(KPIV,j) = AKJ
          end do 
        end if
C   Pivot Value
        DPIV   = 1.0/A(k,k) 
        A(k,k) = DPIV
C   A(*,k)=A(*,k)/A(k,k)
        do i=k+1,N
          A(i,k) = A(i,k)*DPIV
        end do
C   Main Elimination
        do j=k+1,N
          AKJ = A(k,j)
          do i=k+1,N
            A(i,j) = A(i,j) - AKJ*A(i,k)
          end do
        end do
      end do
C---- Solve LUx=b by Substitution -------------
C  Interchange Entries of B
      do k=1,N-1
        i    = IP(k)
        BW   = B(k)
        B(k) = B(i)
        B(i) = BW
      end do
C  Forward Substitution
      do j=1,N-1
        do i=j+1,N
          B(i) = B(i) - B(j)*A(i,j)
        end do
      end do
C  Back Substitution
      do j=N,1,-1
        B(j) = B(j)*A(j,j)
        do i=1,j-1
          B(i) = B(i) - A(i,j)*B(j)
        end do
      end do
C
  100 continue
C
      RETURN
      END