! Updated 10/24/2001. ! !!!!!!!!!!!!!!!!!!!!!!!!!!! Program 4.2 !!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! ! 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. ! ! ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! SUBROUTINE DTRM (A,N,D,INDX) ! ! Subroutine for evaluating the determinant of a matrix using ! the partial-pivoting Gaussian elimination scheme. ! Copyright (c) Tao Pang 2001. ! IMPLICIT NONE INTEGER, INTENT (IN) :: N INTEGER :: I,J,MSGN INTEGER, INTENT (OUT), DIMENSION (N) :: INDX REAL, INTENT (OUT) :: D REAL, INTENT (INOUT), DIMENSION (N,N) :: A ! CALL ELGS(A,N,INDX) ! D = 1.0 DO I = 1, N D = D*A(INDX(I),I) END DO ! MSGN = 1 DO I = 1, N DO WHILE (I.NE.INDX(I)) MSGN = -MSGN J = INDX(I) INDX(I) = INDX(J) INDX(J) = J END DO END DO D = MSGN*D END SUBROUTINE DTRM ! 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