ブロックLU分解FORTRANプログラムNO.2

1. 概要

2. プログラム

C=================================================================C
      SUBROUTINE BLU2(A,B,N,ND,MB,EPS,IP,IER)
C=================================================================C
C  Real-Dense LU Decomposition by Block 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, Input and Output Matrix                   C
C    B(N)     R*8, I/O, Right vector and Solution Vector          C
C    N        I*4, In,  Matrix Order of A                         C
C    ND       I*4, In,  Array Size of A ( ND >= N )               C
C    MB       I*4, In.  Block Size                                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/09/04                     C
C        ( Hitachi Ltd. and Waseda University )                   C
C     Ver.1   No tuning Version                                   C
C=================================================================C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION A(ND,N), B(N), IP(N)
C----- Block Gauss Elimination Step ------------ 
C  Check Parameter
      IER = 0
      if(ND.lt.N .or. MB.lt.1) then
        IER = 2
        go to 100
      end if
C  Gauss Elimination
      do kk=1,N,MB
        MBK = MIN(kk+MB-1,N)
C   Block Process-1 
        do k=kk,MBK
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,*)
          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    Lower-Block Elimination
          do j=k+1,MBK
            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   Block Process-2 ( Upper-Block Elimination )
        do k=MBK+1,N
          do j=kk,MBK-1
            AJK = A(j,k)
            do i=j+1,MBK
              A(i,k) = A(i,k) - AJK*A(i,j)
            end do
          end do
        end do
C   Block Process-3 ( Main Elimination by C=C-A*B) )
        do j=MBK+1,N
          do k=kk,MBK
            AKJ = A(k,j)
            do i=MBK+1,N
              A(i,j) = A(i,j) - AKJ*A(i,k)
            end do 
          end do
        end do
      end do
C---- Solve LUx=b by Substitution -------------
C  Interchange Entries of B
      do j=1,N-1
        k    = IP(j)
        BW   = B(k)
        B(k) = B(j)
        B(j) = BW
      end do               
C  Forward Substitution
      do j=1,N-1
        BJ = B(j)
        do i=j+1,N
          B(i) = B(i) - BJ*A(i,j)
        end do
      end do
C  Backword Substitution
      do j=N,1,-1
        BJ   = B(j)*A(j,j)
        B(j) = BJ
        do i=1,j-1
          B(i) = B(i) - BJ*A(i,j)
        end do 
      end do
C
  100 continue
C
      RETURN
      END