Gauss-Jordan Elimination

Gauss-Jordan Elimination

MODULE GJE

CONTAINS

SUBROUTINE simul(a, b, ndim, n, error)

IMPLICIT NONE

INTEGER, INTENT(IN) :: ndim
REAL, INTENT(INOUT), DIMENSION(ndim, ndim) :: a
REAL, INTENT(INOUT), DIMENSION(ndim) :: b
INTEGER, INTENT(IN) :: n
INTEGER, INTENT(OUT) :: error

REAL, PARAMETER :: EPSILON = 1.0E-6

REAL :: factor

INTEGER :: irow
INTEGER :: ipeak
INTEGER :: jrow
INTEGER :: kcol

REAL :: TEMP

INTEGER :: i, j

DO irow = 1, n

ipeak = irow

! 하삼각 행렬 중 해당 열에서 최대값 찾기
DO jrow = irow + 1,n
IF (ABS(a(jrow, irow)) > ABS(a(ipeak, irow))) THEN
ipeak = jrow
END IF
END DO

IF (ABS(a(ipeak, irow)) < EPSILON) THEN
error = 1
RETURN
END IF

! 최대값 행과 현재 행 바꾸기
IF (ipeak /= irow) THEN
DO kcol = 1, n
temp = a(ipeak, kcol)
a(ipeak, kcol) = a(irow, kcol)
a(irow, kcol) = temp
END DO
temp = b(ipeak)
b(ipeak) = b(irow)
b(irow) = temp
END IF

WRITE (*,*) "Coefficients after swapping rows :"
DO i = 1, n
WRITE (*,"(1X,7F11.4)") (a(i,j), j = 1, n), b(i)
END DO

! 현재 행을 제외하고 열을 0로 만들기
DO jrow = 1, n
IF (jrow /= irow) THEN
factor = -a(jrow, irow)/a(irow,irow)
DO kcol = 1, n
a(jrow, kcol) = a(irow, kcol) * factor + a(jrow, kcol)
END DO
b(jrow) = b(irow) * factor + b(jrow)
END IF
END DO

WRITE (*,*) "Coefficients after sweepping :"
DO i = 1, n
WRITE (*,"(1X,7F11.4)") (a(i,j), j = 1, n), b(i)
END DO

END DO

! 대각원소를 1로 만들기
DO irow = 1, n
b(irow) = b(irow) / a(irow, irow)
a(irow, irow) = 1.
END DO

error = 0

END SUBROUTINE simul

END MODULE GJE

위 소스를 gje.f95로 저장

PROGRAM test_simul

USE gje

IMPLICIT NONE

!constants
INTEGER, PARAMETER :: MAX_SIZE = 10 ! max number of equations

!local variables
REAL, DIMENSION(MAX_SIZE, MAX_SIZE) :: a ! left hand side
REAL, DIMENSION(MAX_SIZE) :: b ! right hand side
INTEGER :: error
CHARACTER(LEN = 20) :: file_name
INTEGER :: i, j ! loop index
INTEGER :: n ! number of simul equations
INTEGER :: istat ! I/O status

! get the file name
WRITE (*,*) "Enter the file name containing the equations:"
READ (*,'(A20)') file_name

! open file
OPEN (UNIT = 11, FILE = file_name, STATUS = 'OLD', ACTION = 'READ', IOSTAT = istat)

! was the OPEN successful?
IF (istat == 0) THEN
READ (11, *) n

IF (n <= MAX_SIZE) THEN
DO i = 1, n
READ (11, *) (a(i, j), j=1, n), b(i)
END DO

WRITE (*,*) "Coefficients before call:"
DO i = 1, n
WRITE (*, '(1X, 7F11.4)') (a(i, j), j = 1, n), b(i)
END DO

CALL simul(a, b, MAX_SIZE, n, error)

IF (error /= 0) THEN
WRITE (*, 1010)
1010 FORMAT (/1X, 'Zero pivot encountered!', &
//1X, 'There is no unique solution to this system.')
ELSE
WRITE (*,*) "Coefficients after call :"
DO i = 1, n
WRITE (*,"(1X,7F11.4)") (a(i,j), j = 1, n), b(i)
END DO

WRITE(*,*) "The Solutions are :"
DO i = 1, n
WRITE (*,"(3X,'X(', I2,') = ', F16.6)") i, b(i)
END DO

END IF
ELSE
WRITE (*,*) "The number of equatons is less than or equal to ", MAX_SIZE
END IF
ELSE
WRITE (*, 1020) istat
1020 FORMAT (1X, 'File open failed -- status = ', I6)
END IF

END PROGRAM test_simul

위 소스를 test_simul.f95로 저장

3
1. 1. 1. 1.
2. 1. 1. 2.
1. 3. 2. 4.

위 데이터를 inputs1.txt로 저장

컴파일 및 테스트