! Updated 10/24/2001.
!
!!!!!!!!!!!!!!!!!!!!!!!!!!! Program 4.3 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! !
! Please Note: !
! !
! (1) This computer program is written by Tao Pang in conjunction with !
! his book, "An Introduction to Computational Physics," published !
! by Cambridge University Press in 1997. !
! !
! (2) No warranties, express or implied, are made for this program. !
! !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
PROGRAM EX43
!
!
! An example of solving linear equation set A(N,N)*X(N) = B(N)
! with the partial-pivoting Gaussian elimination scheme. The
! numerical values are for the Wheatstone bridge example discussed
! in Section 4.1 in the book with all resistances being 100 ohms
! and the voltage 200 volts.
!
IMPLICIT NONE
INTEGER, PARAMETER :: N=3
INTEGER :: I,J
INTEGER, DIMENSION (N) :: INDX
REAL, DIMENSION (N) :: X,B
REAL, DIMENSION (N,N) :: A
DATA B /200.0,0.0,0.0/, &
((A(I,J), J=1,N),I=1,N) /100.0,100.0,100.0,-100.0, &
300.0,-100.0,-100.0,-100.0, 300.0/
!
CALL LEGS (A,N,B,X,INDX)
!
WRITE (6, "(F16.8)") (X(I), I=1,N)
END PROGRAM EX43
SUBROUTINE LEGS (A,N,B,X,INDX)
!
! Subroutine to solve the equation A(N,N)*X(N) = B(N) with the
! partial-pivoting Gaussian elimination scheme.
! Copyright (c) Tao Pang 2001.
!
IMPLICIT NONE
INTEGER, INTENT (IN) :: N
INTEGER :: I,J
INTEGER, INTENT (OUT), DIMENSION (N) :: INDX
REAL, INTENT (INOUT), DIMENSION (N,N) :: A
REAL, INTENT (INOUT), DIMENSION (N) :: B
REAL, INTENT (OUT), DIMENSION (N) :: X
!
CALL ELGS (A,N,INDX)
!
DO I = 1, N-1
DO J = I+1, N
B(INDX(J)) = B(INDX(J))-A(INDX(J),I)*B(INDX(I))
END DO
END DO
!
X(N) = B(INDX(N))/A(INDX(N),N)
DO I = N-1, 1, -1
X(I) = B(INDX(I))
DO J = I+1, N
X(I) = X(I)-A(INDX(I),J)*X(J)
END DO
X(I) = X(I)/A(INDX(I),I)
END DO
!
END SUBROUTINE LEGS
!
SUBROUTINE ELGS (A,N,INDX)
!
! Subroutine to perform the partial-pivoting Gaussian elimination.
! A(N,N) is the original matrix in the input and transformed matrix
! plus the pivoting element ratios below the diagonal in the output.
! INDX(N) records the pivoting order. Copyright (c) Tao Pang 2001.
!
IMPLICIT NONE
INTEGER, INTENT (IN) :: N
INTEGER :: I,J,K,ITMP
INTEGER, INTENT (OUT), DIMENSION (N) :: INDX
REAL :: C1,PI,PI1,PJ
REAL, INTENT (INOUT), DIMENSION (N,N) :: A
REAL, DIMENSION (N) :: C
!
! Initialize the index
!
DO I = 1, N
INDX(I) = I
END DO
!
! Find the rescaling factors, one from each row
!
DO I = 1, N
C1= 0.0
DO J = 1, N
C1 = AMAX1(C1,ABS(A(I,J)))
END DO
C(I) = C1
END DO
!
! Search the pivoting (largest) element from each column
!
DO J = 1, N-1
PI1 = 0.0
DO I = J, N
PI = ABS(A(INDX(I),J))/C(INDX(I))
IF (PI.GT.PI1) THEN
PI1 = PI
K = I
ENDIF
END DO
!
! Interchange the rows via INDX(N) to record pivoting order
!
ITMP = INDX(J)
INDX(J) = INDX(K)
INDX(K) = ITMP
DO I = J+1, N
PJ = A(INDX(I),J)/A(INDX(J),J)
!
! Record pivoting ratios below the diagonal
!
A(INDX(I),J) = PJ
!
! Modify other elements accordingly
!
DO K = J+1, N
A(INDX(I),K) = A(INDX(I),K)-PJ*A(INDX(J),K)
END DO
END DO
END DO
!
END SUBROUTINE ELGS