MODULE Precisions INTEGER,PARAMETER:: Long=SELECTED_REAL_KIND(15,307) END MODULE Precisions PROGRAM c2406 USE Precisions IMPLICIT NONE INTEGER :: I,N REAL (Long), ALLOCATABLE:: A(:,:),B(:),X(:) LOGICAL:: Singular INTERFACE SUBROUTINE Gaussian_Elimination(A,N,B,X,Singular) USE Precisions IMPLICIT NONE INTEGER, INTENT(IN)::N REAL (Long), INTENT (INOUT) :: A(:,:),B(:) REAL (Long), INTENT(OUT)::X(:) LOGICAL, INTENT(OUT) :: Singular END SUBROUTINE Gaussian_Elimination END INTERFACE PRINT *,'Number of equations?' READ *,N ALLOCATE(A(1:N,1:N),B(1:N),X(1:N)) DO I=1,N PRINT *,'Input elements of row ',I,' of A' READ*,A(I,1:N) PRINT*,'Input element ',I,' of B' READ *,B(I) END DO CALL Gaussian_Elimination(A,N,B,X,Singular) IF(Singular) THEN PRINT*, 'Matrix is singular' ELSE PRINT*, 'Solution X:' PRINT*,X(1:N) ENDIF END PROGRAM c2406 SUBROUTINE Gaussian_Elimination(A,N,B,X,Singular) ! Routine to solve a system Ax=b ! using Gaussian Elimination ! with partial pivoting ! The code is based on the Linpack routines ! SGEFA and SGESL ! and operates on columns rather than rows! USE Precisions IMPLICIT NONE ! Matrix A and vector B are over-written ! Arguments INTEGER, INTENT(IN):: N REAL (Long),INTENT(INOUT):: A(:,:),B(:) REAL (Long),INTENT(OUT)::X(:) LOGICAL,INTENT(OUT)::Singular ! Local variables INTEGER::I,J,K,Pivot_row REAL (Long):: Pivot,Multiplier,Sum,Element REAL (Long),PARAMETER::Eps=1.E-13_Long ! ! Work through the matrix column by column ! DO K=1,N-1 ! ! Find largest element in column K for pivot ! Pivot_row = MAXVAL( MAXLOC( ABS( A(K:N,K) ) ) ) & + K - 1 ! ! Test to see if A is singular ! if so return to main program ! IF(ABS(A(Pivot_row,K)) <= Eps) THEN Singular=.TRUE. RETURN ELSE Singular = .FALSE. ENDIF ! ! Exchange elements in column K if largest is ! not on the diagonal ! IF(Pivot_row /= K) THEN Element=A(Pivot_row,K) A(Pivot_Row,K)=A(K,K) A(K,K)=Element Element=B(Pivot_row) B(Pivot_row)=B(K) B(K)=Element ENDIF ! ! Compute multipliers ! elements of column K below diagonal ! are set to these multipliers for use ! in elimination later on ! A(K+1:N,K) = A(K+1:N,K)/A(K,K) ! ! Row elimination performed by columns for efficiency ! DO J=K+1,N Pivot = A(Pivot_row,J) IF(Pivot_row /= K) THEN ! Swap if pivot row is not K A(Pivot_row,J)=A(K,J) A(K,J)=Pivot ENDIF A(K+1:N,J)=A(K+1:N,J)-Pivot* A(K+1:N,K) END DO ! ! Apply same operations to B ! B(K+1:N)=B(K+1:N)-A(K+1:N,K)*B(K) END DO ! ! Backward substitution ! DO I=N,1,-1 Sum = 0.0 DO J= I+1,N Sum=Sum+A(I,J)*X(J) END DO X(I)=(B(I)-Sum)/A(I,I) END DO END SUBROUTINE Gaussian_Elimination