! GAUSSIAN ELIMINATION DEMONSTRATION - BETTER VERSION ! PROGRAM P142 IMPLICIT NONE REAL :: M1(5,6),SOL(5),TOL INTEGER :: N INTERFACE SUBROUTINE MATSIN(M1,N) IMPLICIT NONE REAL ,INTENT(IN OUT) :: M1(:,:) INTEGER ,INTENT(IN OUT) :: N END SUBROUTINE MATSIN SUBROUTINE GAUSS(M1,SOL,N,TOL) IMPLICIT NONE REAL ,INTENT(IN OUT) :: M1(:,:),SOL(:),TOL INTEGER ,INTENT(IN OUT) :: N SUBROUTINE PRNMAT(M3,N) IMPLICIT NONE REAL ,INTENT(IN OUT) :: M3(:,:) INTEGER ,INTENT(IN OUT) :: N END SUBROUTINE PRNMAT END SUBROUTINE GAUSS END INTERFACE ! N=5 ! NUMBER OF EQUATIONS TOL=.001 ! PRINT *,' This is Program P142 - Gaussian elimination' PRINT *,'PROGRAM IS READING DATA INTO ARRAYS' CALL MATSIN(M1,N) PRINT *,'SOLVING SYSTEM OF EQUATIONS' CALL GAUSS(M1,SOL,N,TOL) IF(TOL<0) THEN PRINT *,'SYSTEM OF EQUATIONS HAS NO SINGLE SOLUTION' STOP ENDIF PRINT *,'SOLUTION:' PRINT 90,SOL 90 FORMAT(' | ',F8.3,' |') STOP END PROGRAM P142 SUBROUTINE MATSIN(M1,N) IMPLICIT NONE REAL ,INTENT(IN OUT) :: M1(:,:) INTEGER ,INTENT(IN OUT) :: N INTEGER :: I,J ! ! Tell program where data for READ is coming from OPEN(UNIT=5, FILE='P142.DAT') ! ! READ IN M1 ! ONE ROW PER CARD ! L1: DO I=1,N READ 27,(M1(I,J),J=1,N+1) END DO L1 27 FORMAT(10(F5.2)) RETURN END SUBROUTINE MATSIN ! SUBROUTINE PRNMAT(M3,N) IMPLICIT NONE REAL ,INTENT(IN OUT) :: M3(:,:) INTEGER ,INTENT(IN OUT) :: N INTEGER :: I,J ! ! L3: DO I=1,N PRINT 202, (M3(I,J),J=1,N+1) END DO L3 202 FORMAT(10(' ',F7.3)) ! RETURN END SUBROUTINE PRNMAT ! SUBROUTINE GAUSS(M1,SOL,N,TOL) IMPLICIT NONE REAL ,INTENT(IN OUT) :: M1(:,:),SOL(:),TOL INTEGER ,INTENT(IN OUT) :: N ! ! THIS ROUTINE PERFORMS GAUSSIAN ELIMINATION AND BACKSUBSTITUTION. ! IN THIS VERSION, WE CONSIDER THE CASE WHERE THE SYSTEM OF EQUATIONS ! HAS NO SINGLE SOLUTION ( INFINITELY MANY OR NONE ). TO MAKE IT EASY ! TO CHECK FOR THIS POSSIBILITY, THE COMPUTER CHOOSES WHICH ROW HAS THE ! GREATEST LEADING CO-EFFICIENT, AND USES THIS ROW IN THE ELIMINATION ! PROCESS. IF IT CANNOT FIND A NON-ZERO ROW (ZERO WITHIN A TOLERANCE ! SET BY THE ROUTINE WHICH CALLS THIS SUBROUTINE), IT RETURNS TO THE ! CALLING ROUTINE WITH A TOL=-90 FLAG, AND THE USER IS TOLD THAT THE ! SYSTEM HAS NO SOLUTION. ! REAL :: M2(:,:),TEMP,MAX INTEGER :: SWAP,I,J,K ! ! ! INSTEAD OF MODIFYING THE ORIGINAL ARRAY, WE WILL PRODUCE A WORKING COPY ! OF IT ! M2 = M1 ! L1: DO I=1,N MAX=-3.0 DO K=I,N IF(ABS(M2(K,I))>MAX) THEN MAX=ABS(M2(K,I)) SWAP=K ENDIF END DO IF(MAXI) THEN DO M=I,N+1 TEMP=M2(I,M) M2(I,M)=M2(SWAP,M) M2(SWAP,M)=TEMP ENDDO ENDIF ! L2: DO J=I+1,N TEMP=M2(J,I)/M2(I,I) L3: DO K=I,N+1 M2(J,K)= M2(J,K) - TEMP*M2(I,K) END DO L3 END DO L2 END DO L1 ! PRINT *,'TRIANGULARIZED MATRIX' CALL PRNMAT(M2,5) ! ! MATRIX IS NOW TRIANGULAR. USE BACKSUBSTITUTION TO SOLVE ! SOL(N)=M2(N,N+1)/M2(N,N) L4: DO I=N-1,1,-1 TEMP=0.0 L5: DO K=N,I+1,-1 TEMP=TEMP+M2(I,K)*SOL(K) END DO L5 SOL(I)=(M2(I,N+1) - TEMP)/M2(I,I) END DO L4 RETURN END SUBROUTINE GAUSS