! ! PROGRAM GAUSS SOLVES A SYSTEM OF EQUATIONS BY ! USING GAUSS ELIMINATION AND U-L DECOMPOSITION. ! PROGRAM GAUSS IMPLICIT NONE REAL, ALLOCATABLE, DIMENSION(:,:) :: A, UL INTEGER, ALLOCATABLE, DIMENSION(:) :: IPS REAL, ALLOCATABLE, DIMENSION(:) :: B, X INTEGER :: I, J, N OPEN(6,FILE='GAUSS.OUT') ! ! SET THE NUMBER OF SIMULTANEOUS EQUATIONS ! N=5 ! ! ALLOCATE ARRAY DIMENSIONS ! ALLOCATE (A(N,N), UL(N,N), IPS(N), B(N), X(N)) ! ! WRITE OUT HEADER DESCRIBING PROBLEM ! WRITE(6,200)N 200 FORMAT(' RON KRIZ TEMP-GRAD F90 HOMEWORK, MSE2094',& &/,TR1,' SOLVES UP TO 50X50 SYSTEM OF EQUATIONS BY LU ',& &/,TR1,' DECOMPOSITION, WHERE IN THIS CASE N=',I2) ! ! ZERO [A] {B} ARRAYS ! DO I=1,N B(I)=0.0 DO J=1,N A(I,J)=0.0 END DO END DO ! ! SET NONZERO ARRAY COMPONENTS ! A(1,1)=1.0 A(2,2)=-33.0/16.0 A(3,3)=-33.0/16.0 A(4,4)=-33.0/16.0 A(5,5)=1.0 A(2,1)=1.0 A(2,3)=1.0 A(3,2)=1.0 A(3,4)=1.0 A(4,3)=1.0 A(4,5)=1.0 ! B(1)=5.0 ! ! SOLVE SYSTEM OF EQUATIONS FOR UNKNOWNS ! CALL DECOMP(N,A,UL,IPS) CALL SOLVE(N,B,UL,IPS,X) ! ! WRITE OUT RESULTS ! DO J=1,N WRITE(6,300)J,X(J) END DO 300 FORMAT(' X(',I2,')=',E15.8) STOP ! ! USING THE CONTAINS STATEMENT YOU CAN AVOID BUILDING A ! FORTRAN90 INTERFACE AT THE BEGINING OF THE MAIN PROGRAM. ! THE CONTAINS STATEMENT INCLUDES THE SUBROUTINE WITH THE ! MAIN PROGRAM, HENCE THE "END PROGRAM GAUSS" IS AT THE ! END OF THIS FILE. ! CONTAINS ! !---------------------------------------------------------------------------------------------------------------- ! ! SUBROUTINE DECOMP, DECOMPOSES THE ! [A] ARRAY INTO [U] AND [L] SUBMATRICES ! SUBROUTINE DECOMP(N,A,UL,IPS) IMPLICIT NONE INTEGER, INTENT(IN) :: N REAL, DIMENSION(:,:), INTENT(IN) :: A REAL, DIMENSION(:,:), INTENT(OUT) :: UL INTEGER, DIMENSION(:), INTENT(OUT) :: IPS REAL, ALLOCATABLE, DIMENSION(:) :: SCALES REAL :: ROWNRM,SIZE,BIG,PIVOT,EM INTEGER :: I,J,K,NM1,IP,IDXP1V,KP,KP1 ALLOCATE (SCALES(N)) DO I=1,N IPS(I)=I ROWNRM=0.0 DO J=1,N UL(I,J)=A(I,J) IF ( ( ROWNRM - ABS( UL(I,J) ) ) .LT. 0.0 ) ROWNRM=ABS(UL(I,J)) END DO IF (ROWNRM .EQ. 0.0) CALL SING(1) IF (ROWNRM .EQ. 0.0) SCALES(I)=0.0 SCALES(I)=1.0/ROWNRM END DO NM1=N-1 DO K=1,NM1 BIG = 0.0 DO I=K,N IP = IPS(I) SIZE = ABS(UL(IP,K))*SCALES(IP) IF ((SIZE -BIG) .LE. 0.0) GO TO 11 BIG = SIZE IDXP1V=I 11 END DO IF (BIG .EQ. 0.0) CALL SING(2) IF (BIG .EQ. 0.0) GO TO 17 IF ((IDXP1V - K) .EQ. 0) GO TO 15 J = IPS(K) IPS(K) = IPS(IDXP1V) IPS(IDXP1V) = J 15 KP=IPS(K) PIVOT = UL(KP,K) KP1 = K+1 DO I=KP1,N IP=IPS(I) EM = -UL(IP,K)/PIVOT UL(IP,K) = -EM DO J=KP1,N UL(IP,J) = UL(IP,J) + EM*UL(KP,J) END DO END DO 17 END DO KP = IPS(N) IF (UL(KP,N) .EQ. 0) CALL SING(2) RETURN END SUBROUTINE DECOMP ! !--------------------------------------------------------------------------------------------------------------- ! ! SUBROUTINE SOLVE SOLVES FOR THE UNKNOWNS BY ! BACK-SUBSTITUTION USING THE DECOMPOSED [U] AND [L] MATRICES. ! SUBROUTINE SOLVE(N,B,UL,IPS,X) INTEGER, INTENT(IN) :: N REAL, DIMENSION(:), INTENT(IN) :: B REAL, DIMENSION(:,:), INTENT(IN) :: UL INTEGER, DIMENSION(:), INTENT(IN) :: IPS REAL, DIMENSION(:), INTENT(OUT) :: X REAL :: SUM INTEGER :: I,J,NP1,IP,IM1,IP1,IBACK NP1=N+1 IP = IPS(1) X(1) = B(IP) DO I=2,N IP = IPS(I) IM1 = I-1 SUM = 0.0 DO J=1,IM1 SUM = SUM + UL(IP,J)*X(J) END DO X(I) = B(IP) - SUM END DO IP = IPS(N) X(N) = X(N)/UL(IP,N) DO IBACK = 2,N I = NP1 - IBACK IP = IPS(I) IP1 = I+1 SUM = 0.0 DO J=IP1,N SUM = SUM + UL(IP,J)*X(J) END DO X(I) = (X(I)-SUM)/UL(IP,I) END DO RETURN END SUBROUTINE SOLVE ! !----------------------------------------------------------------------------------------------------------- ! ! SUBROUTINE SING FLAGS SINGULARITIES IN THE CALCULATIONS ! WHEN SOLVING FOR THE UNKNOWNS AND WRITES MESSAGES ! TO THE OUTPUT FILE INDICATING THE TYPE OF SINGULARITY. ! SUBROUTINE SING(IWHY) INTEGER, INTENT(IN) :: IWHY IF (IWHY .EQ. 1) WRITE(6,11) IF (IWHY .EQ. 2) WRITE(6,12) IF (IWHY .EQ. 3) WRITE(6,13) 11 FORMAT(' MATRIX WITH ZERO ROW IN DECOMPOSE') 12 FORMAT(' SINGULAR MATRIX IN DECOMPOSE. ZERO DIVIDE IN SOLVE') 13 FORMAT(' NO CONVERGENCE IN IMPRUV. MATRIX IS NEARLY SINGULAR') RETURN END SUBROUTINE SING END PROGRAM GAUSS