FORTRAN Ax=b Block Solve Program-1

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

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

1. 概要

　実密行列を係数とする連立一次方程式Ax=bをブロックガウス消去法でLU分解し、
LUx=bの数値解を前進後退代入で求める。
軸交換も特異性のチェックもしていない最も単純なプログラムである。

2. プログラム

C=================================================================C
      SUBROUTINE BLU1(A,B,N,ND,MB)
C=================================================================C
C  Real-Dense LU Decomposition by Block Gauss Elimination         C
C   and Solve Ax=b by Substitution                                C
C    A ---> LU Decomposition without Partial Pivoting             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-----------------------------------------------------------------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)
C----- Block Gauss Elimination Step ------------
C  Gauss Elimination
      do kk=1,N,MB
        MBK = min(kk+MB-1,N)
C   Block Process-1 
        do k=kk,MBK
C    A(*,k)=A(*,k)/A(k,k)
          do i=k+1,N
            A(i,k) = A(i,k)/A(k,k)
          end do
C    Lower-Block Elimination
          do j=k+1,MBK
            do i=k+1,N
              A(i,j) = A(i,j) - A(k,j)*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
            do i=j+1,MBK
              A(i,k) = A(i,k) - A(j,k)*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
            do i=MBK+1,N
              A(i,j) = A(i,j) - A(i,k)*A(k,j)
            end do 
          end do
        end do
      end do
C---- Solve LUx=b by Substitution -------------
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  Backword 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
      RETURN
      END