FORTRAN Ax=b Solve Program-1

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

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

1. 概要

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

2. プログラム

C=================================================================C
      SUBROUTINE GLU1(A,B,N,ND)
C=================================================================C
C  Real-Dense LU Decomposition by Gauss Elimination               C
C   and Solve Ax=b by Substitution                                C
C     A ---> LU Decomposition Only                                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-----------------------------------------------------------------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)
C----- Gauss Elimination Step ------------      
      do k=1,N-1
        do i=k+1,N
          A(i,k) = A(i,k)/A(k,k)
        end do
C   Main Elimination
        do j=k+1,N
          do i=k+1,N
            A(i,j) = A(i,j) - A(k,j)*A(i,k)
          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) - A(i,j)*B(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
      RETURN
      END