PROGRAM GAUSSIAN
PARAMETER (IN=20)
REAL:: A(IN,IN), B(IN)
PRINT *,'INPUT NUMBER OF EQUATIONS N='
READ (5,*) N
PRINT *,'INPUT MATRIX COEFFICIENTS A(I,J)='
READ (5,*) ((A(I,J),J=1,N),I=1,N)
PRINT *,'INPUT RIGHT-HAND SIDE VECTOR B(I)='
READ (5,*) (B(I),I=1,N)
WRITE (6,*) ('****GAUSSIAN ELIMINATION ****')
WRITE (6,*)
WRITE (6,*) ('COEFFICIENT MATRIX you inputed:')
CALL PRINTA(A,IN,N,N,6)
WRITE(6,*)
WRITE(6,*) ('right hand side vector you inputed:')
CALL PRINTV(B,N,6)
WRITE(6,*)
! CONVERT TO UPPER TRIANGULAR FORM
DO K = 1,N-1
IF (ABS(A(K,K)).GT.1.E-6) THEN
DO I = K+1, N
X = A(I,K)/A(K,K)
DO J = K+1, N
A(I,J) = A(I,J) -A(K,J)*X
ENDDO
B(I) = B(I) - B(K)*X
ENDDO
ELSE
WRITE (6,*) 'ZERO PIVOT FOUND IN LINE:'
WRITE (6,*) K
STOP
END IF
ENDDO
WRITE(6,*) 'MODIFIED MATRIX'
CALL PRINTA(A,IN,N,N,6)
WRITE(6,*)
WRITE(6,*) 'MODIFIED RIGHT HAND SIDE VECTOR'
CALL PRINTV (B,N,6)
WRITE(6,*)
! BACK SUBSTITUTION
DO I = N,1,-1
SUM = B(I)
IF (I.LT.N) THEN
DO J= I+1,N
SUM = SUM - A(I,J)*B(J)
ENDDO
END IF
B(I) = SUM/A(I,I)
ENDDO
! PRINT THE RESULTS
write(6,*) ('SOLUTION VECTOR')
CALL PRINTV(B,N,6)
!
END PROGRAM GAUSSIAN
!------------------------------------------
SUBROUTINE PRINTA(A,IA,M,N,ICH)
!
! WRITE A 2D ARRAY TO OUTPUT CHANNEL 'ICH'
!
REAL A(IA,*)
DO I =1,M
WRITE(ICH,2) (A(I,J),J=1,N)
ENDDO
2 FORMAT(1X,6E12.4)
!
END SUBROUTINE PRINTA
!-----------------------------------------
SUBROUTINE PRINTV(VEC,N,ICH)
!
! WRITE A COLUMN VECTOR TO CHANNEL 'ICH'
!
REAL VEC(*)
WRITE(ICH,1) (VEC(I),I=1,N)
1 FORMAT(1X,6E12.4)
!
END SUBROUTINE PRINTV
!-----------------------------------------
PARAMETER (IN=20)
REAL:: A(IN,IN), B(IN)
PRINT *,'INPUT NUMBER OF EQUATIONS N='
READ (5,*) N
PRINT *,'INPUT MATRIX COEFFICIENTS A(I,J)='
READ (5,*) ((A(I,J),J=1,N),I=1,N)
PRINT *,'INPUT RIGHT-HAND SIDE VECTOR B(I)='
READ (5,*) (B(I),I=1,N)
WRITE (6,*) ('****GAUSSIAN ELIMINATION ****')
WRITE (6,*)
WRITE (6,*) ('COEFFICIENT MATRIX you inputed:')
CALL PRINTA(A,IN,N,N,6)
WRITE(6,*)
WRITE(6,*) ('right hand side vector you inputed:')
CALL PRINTV(B,N,6)
WRITE(6,*)
! CONVERT TO UPPER TRIANGULAR FORM
DO K = 1,N-1
IF (ABS(A(K,K)).GT.1.E-6) THEN
DO I = K+1, N
X = A(I,K)/A(K,K)
DO J = K+1, N
A(I,J) = A(I,J) -A(K,J)*X
ENDDO
B(I) = B(I) - B(K)*X
ENDDO
ELSE
WRITE (6,*) 'ZERO PIVOT FOUND IN LINE:'
WRITE (6,*) K
STOP
END IF
ENDDO
WRITE(6,*) 'MODIFIED MATRIX'
CALL PRINTA(A,IN,N,N,6)
WRITE(6,*)
WRITE(6,*) 'MODIFIED RIGHT HAND SIDE VECTOR'
CALL PRINTV (B,N,6)
WRITE(6,*)
! BACK SUBSTITUTION
DO I = N,1,-1
SUM = B(I)
IF (I.LT.N) THEN
DO J= I+1,N
SUM = SUM - A(I,J)*B(J)
ENDDO
END IF
B(I) = SUM/A(I,I)
ENDDO
! PRINT THE RESULTS
write(6,*) ('SOLUTION VECTOR')
CALL PRINTV(B,N,6)
!
END PROGRAM GAUSSIAN
!------------------------------------------
SUBROUTINE PRINTA(A,IA,M,N,ICH)
!
! WRITE A 2D ARRAY TO OUTPUT CHANNEL 'ICH'
!
REAL A(IA,*)
DO I =1,M
WRITE(ICH,2) (A(I,J),J=1,N)
ENDDO
2 FORMAT(1X,6E12.4)
!
END SUBROUTINE PRINTA
!-----------------------------------------
SUBROUTINE PRINTV(VEC,N,ICH)
!
! WRITE A COLUMN VECTOR TO CHANNEL 'ICH'
!
REAL VEC(*)
WRITE(ICH,1) (VEC(I),I=1,N)
1 FORMAT(1X,6E12.4)
!
END SUBROUTINE PRINTV
!-----------------------------------------
Ada komentar, kritik, saran, atau request?