! Last change: DOS 24 Jul 2000 9:03 pm ! *** copyright 2000 *** ! *** filename gaucpp.f95 *** John F. Monahan ** ! ********************** PROGRAM PGAUCPP ! gaucpp -- gauss elimination with partial pivoting with condition ! estimate ! SOLVE SYSTEM OF LINEAR EQUATIONS A * X = B BY GAUSSIAN ELIMINATION ! GAUSPP COMPUTES THE L U FACTORIZATION OF THE MATRIX A ! GAUSSE SOLVES FOR A PARTICULAR RHS B ! ! DECLARATIONS IMPLICIT NONE INTEGER, PARAMETER :: LENGTH = 10 REAL, DIMENSION(LENGTH,LENGTH) :: A REAL, DIMENSION(LENGTH) :: B REAL rcond INTEGER, DIMENSION(LENGTH) :: PI INTEGER I,J,N ! INTERFACE BLOCK INTERFACE subroutine gaucpp(a,n,pi,rcond) IMPLICIT NONE REAL, DIMENSION(:,:), INTENT(IN OUT) :: a REAL, INTENT(OUT) :: rcond INTEGER, INTENT (IN) :: n INTEGER, DIMENSION(:), INTENT(OUT) :: pi end subroutine gaucpp SUBROUTINE GAUSSE(A,N,PI,B) IMPLICIT NONE REAL, DIMENSION(:,:), INTENT(IN) :: A REAL, DIMENSION(:), INTENT(IN OUT) :: B INTEGER, INTENT (IN) :: N INTEGER, DIMENSION(:), INTENT(IN) :: PI END SUBROUTINE GAUSSE END INTERFACE ! FORMAT STATEMENTS 20 FORMAT(2X,4F6.2) 21 FORMAT(1X,10F13.8) 22 FORMAT(2X,'PERMUTATION VECTOR ',10I4) 25 format(' reciprocal of condition number',f13.8) 26 FORMAT(' INPUT MATRIX ') 27 FORMAT(' FACTORED MATRIX -- DIAGONAL IS PART OF U ') 28 FORMAT(' RIGHT HAND SIDE ') 29 FORMAT(' SOLUTION VECTOR ') ! OPEN( UNIT=5, FILE='gauspp.dat' ) OPEN( UNIT=6, FILE='gaucpp.out' ) ! READ THE MATRIX A N = 4 DO I = 1,N READ(5,20) (A(I,J),J=1,N) ENDDO ! WRITE IT OUT WRITE(6,26) DO I = 1,N WRITE(6,21) (A(I,J),J=1,N) ENDDO ! DO THE GAUSSIAN ELIMINATION FACTORIZATION ! with condition estimate CALL GAUCPP(A,N,PI,rcond) WRITE(6,22) (PI(J),J=1,N) write(6,25) rcond ! WRITE(6,27) DO I = 1,N WRITE(6,21) (A(I,J),J=1,N) ENDDO ! ! READ IN THE RHS B READ(5,20) (B(J),J=1,N) WRITE(6,28) WRITE(6,21) (B(J),J=1,N) ! CALL GAUSSE(A,N,PI,B) ! WRITE(6,29) WRITE(6,21) (B(J),J=1,N) ! ! STOP END PROGRAM PGAUCPP subroutine gaucpp(a,n,pi,rcond) ! computes the LU factorization of the matrix a using gaussian ! elimination with partial pivoting and computes condtion estimate IMPLICIT NONE REAL, DIMENSION(:,:), INTENT(IN OUT) :: a REAL, INTENT(OUT) :: rcond INTEGER, INTENT (IN) :: n INTEGER, DIMENSION(:), INTENT(OUT) :: pi ! REAL, DIMENSION(n) :: p,y real anorm,rnorm,znorm,ykp,ykm,sp,sm integer i,k,kk,kp1 ! ! a real input: matrix of real dimension n by n ! output: LU factorization, PA = LU ! n integer dimension of matrix ! pi integer vector of permutations -- forms P ! rcond real reciprocal of condition number of a using ! infinity norm ! ! first step is to compute the infinity norm of the matrix A ! then the LU factorization is computed, overwriting A, using GAUSPP ! then compute the vector y which leads to large norm estimate ! using Golub and van Loan Algorithm 4.5-1 ! then solve Az=y to finish using GAUSSE ! ! See Golub and van Loan, Matrix Computations, pp.76-78 ! Cline, Moler, Stewart, and Wilkinson (1979) An Estimate of the ! Condition Number of a Matrix, SIAM J Num Anal, 16, p368ff ! Dongarra, Bunch, Moler and Stewart (1979) Linpack User's Guide ! ! J F Monahan (June, 1992) Dept of Stat, NCSU, Raleigh, NC 27695-8203 ! revised June 1999, April 2000 for Fortran 95 ! ! INTERFACE BLOCK INTERFACE SUBROUTINE GAUSPP(A,N,PI) IMPLICIT NONE REAL, DIMENSION(:,:), INTENT(IN OUT) :: A INTEGER, INTENT (IN) :: N INTEGER, DIMENSION(:), INTENT(OUT) :: PI END SUBROUTINE GAUSPP SUBROUTINE GAUSSE(A,N,PI,B) IMPLICIT NONE REAL, DIMENSION(:,:), INTENT(IN) :: A REAL, DIMENSION(:), INTENT(IN OUT) :: B INTEGER, INTENT (IN) :: N INTEGER, DIMENSION(:), INTENT(IN) :: PI END SUBROUTINE GAUSSE END INTERFACE rcond = 0. ! compute the norm of A anorm = 0. do k = 1,n sm = 0. do i = 1,n sm = sm + abs( a(k,i) ) enddo ! loop on i anorm = max( anorm, sm ) enddo ! loop on k ! call gauspp to compute the factorization P A = L U call gauspp(a,n,pi) ! if singular, then return as is with rcond = 0. if( pi(n) .ne. n ) return ! have L and U factors stored in A, permutation in PI ! apply algorithm 4.5-1 to transpose(U) do i = 1,n p(i) = 0. enddo ! loop on i ! now the main loop rnorm = 1. do k = 1,n ykp = (1./rnorm - p(k) ) / a(k,k) ykm = - (1./rnorm + p(k) ) / a(k,k) sp = abs( ykp ) sm = abs( ykm ) kp1 = k + 1 if( kp1 .le. n ) then do i = kp1,n sp = sp + abs( (p(i) + a(i,k)*ykp)/a(i,i) ) sm = sm + abs( (p(i) + a(i,k)*ykm)/a(i,i) ) enddo ! loop on i endif ! have the best y(k) y(k) = ykp if( sp .ge. sm ) y(k) = ykm if( kp1 .le. n ) then do i = kp1,n p(i) = p(i) + a(i,k)*y(k) enddo ! loop on i endif ! normalize every time -- these grow fast! sm = abs( y(k) ) if( sm .gt. rnorm ) then rnorm = sm do i = 1,n y(i) = y(i) / rnorm p(i) = p(i) / rnorm enddo ! loop on i endif ! end of main loop enddo ! loop on k ! now solve transpose(L) r = y but overwrite do kk = 1,n k = n + 1 - kk if( k .ne. n ) then kp1 = k + 1 do i = kp1,n y(k) = y(k) - a(i,k)*y(i) enddo ! loop on i endif enddo ! loop on kk, k ! compute norm of r rnorm = 0. do i = 1,n rnorm = max( rnorm, abs( y(i) ) ) enddo ! loop on i ! solve a system in A call gausse(a,n,pi,y) ! compute norm of this solution (z) znorm = 0. do i = 1,n znorm = max( znorm, abs( y(i) ) ) enddo ! compute the reciprocal of condition number rcond = rnorm / (anorm*znorm) return end subroutine gaucpp SUBROUTINE GAUSPP(A,N,PI) ! SUBROUTINE TO DO THE LU FACTORIZATION OF GAUSSIAN ELIMINATION ! FACTORS P * A = L U WHERE P IS THE PERMUTATION MATRIX STORED IN PI ! A IS DESTROYED AS L AND U ARE OVERWRITTEN ON IT ! L IS UNIT LOWER AND U UPPER TRIANGULAR ! ! USE GAUSSE TO SOLVE A SYSTEM OF LINEAR EQUATIONS ! ! A IN MATRIX OF COEFFICIENTS N BY N ! OUT FACTORS L AND U OVERLAID ! N IN TRUE SIZE OF MATRIX A (NOT DECLARED SIZE) ! OUT UNCHANGED ! PI IN ANYTHING ! OUT PERMUTATION LIST OF N-1 PIVOT LOCATIONS ! PI(N) GIVES COMPUTED RANK, PI(N) = N IF NONSINGULAR ! ! J F MONAHAN, SEPT 1985, N C S U, RALEIGH, NC 27695-8203 USA ! REVISED JUNE 1999, APRIL 2000 FOR FORTRAN 95 ! IMPLICIT NONE REAL, DIMENSION(:,:), INTENT(IN OUT) :: A INTEGER, INTENT (IN) :: N INTEGER, DIMENSION(:), INTENT(OUT) :: PI REAL S,BIG INTEGER I,J,K,KP1,NM1,IBIG ! formats and writes (lower case) for demo only 21 format(8x,10f13.8) 22 format(2x,'within gauspp, switch rows',i3,' and',i3 & & /2x,'within gauspp, current matrix is') ! PI(N) = 0 IF( N .GT. 1 ) THEN ! DO NM1 STEPS TO FORM A TRIANGULAR SYSTEM NM1 = N - 1 ! DO K = 1,NM1 ! SEARCH FOR BEST PIVOT IN COLUMN K BIG = 0. IBIG = 0 DO I = K,N S = ABS( A(I,K) ) IF( S .GT. BIG ) THEN ! KEEP TRACK OF BEST I AND ITS SIZE BIG = S IBIG = I ENDIF ENDDO ! LOOP ON I ! TEST FOR SINGULARITY ! CHANGE 0. TO AN EPSILON IF DESIRED IF( BIG .EQ. 0. ) RETURN ! KEEP TRACK OF THE RANK IN PI(N) PI(N) = K ! STORE PIVOT LOCATION IN PERMUTATION LIST PI(K) = IBIG write(6,22) k,ibig ! SWITCH ONLY IF NECESSARY IF( IBIG .NE. I ) THEN DO J = 1,N S = A(K,J) A(K,J) = A(IBIG,J) A(IBIG,J) = S ENDDO ! LOOP ON J ! ENDIF ! ELIMINATE REMAINDER OF COLUMN K AND CREATE ONE KP1 = K + 1 DO I = KP1,N A(I,K) = A(I,K) / A(K,K) DO J = KP1,N A(I,J) = A(I,J) - A(I,K) * A(K,J) ENDDO ! LOOP ON J ENDDO ! LOOP ON I do i = 1,n write(6,21) (a(i,j),j=1,n) enddo ! loop on i ENDDO ! END OF MAIN LOOP ON K ENDIF ! TEST LAST (AND MAYBE ONLY) ONE FOR SINGULARITY IF ( ABS( A(N,N) ) .GT. 0. ) THEN ! ONLY (COMPUTATIONALLY) NONSINGULAR MATRICES ! REACH HERE PI(N) = N ENDIF RETURN END SUBROUTINE GAUSPP SUBROUTINE GAUSSE(A,N,PI,B) ! ! SUBROUTINE TO SOLVE A SYSTEM OF LINEAR EQUATIONS A * X = B ! ! A IN HOLDS FACTORS L AND U OVERLAID FROM GAUSPP ! OUT UNCHANGED ! N IN TRUE SIZE OF MATRIX A (NOT DECLARED SIZE) ! OUT UNCHANGED ! PI IN LIST OF PIVOTS -- HOLDS PERMUTATION MATRIX P ! IN FACTORED FORM ! OUT UNCHANGED ! ! *** GAUSPP MUST HAVE BEEN CALLED PREVIOUSLY *** ! ! GAUSPP FACTORS P * A = L U ! TO SOLVE A X = B, FIRST PERMUTE B TO GET P A X = L U X = P B ! THEN SOLVE TWO TRIANGULAR SYSTEMS OF EQUATIONS ! ! J F MONAHAN, SEPT 1985, N C S U, RALEIGH, NC 27695 - 8203 USA ! REVISED JUNE, 1999 FOR FORTRAN 95 ! ! DECLARATIONS FIRST IMPLICIT NONE REAL, DIMENSION(:,:), INTENT(IN OUT) :: A REAL, DIMENSION(:), INTENT(IN OUT) :: B INTEGER, INTENT (IN) :: N INTEGER, DIMENSION(:), INTENT(IN) :: PI REAL S INTEGER I,J,K,II,IP1,IM1,NM1 ! IF( N .GT. 1 ) THEN NM1 = N - 1 ! PERMUTE B FIRST DO K = 1,NM1 S = B(K) I = PI(K) B(K) = B(I) B(I) = S ENDDO ! LOOP ON K ! formats and writes (lower case) for demo only 21 format(8x,'within gausse, intermediate results '/10f13.8) write(6,21) (b(j),j=1,k) ! L IS UNIT TRIANGULAR MATRIX ! SOLVE L * Y = P * B FROM THE TOP DOWN DO I = 2,N IM1 = I - 1 DO J = 1,IM1 B(I) = B(I) - A(I,J) * B(J) ENDDO ! LOOP ON J ENDDO ! LOOP ON I write(6,21) (b(j),j=1,k) ! U IS UPPER TRIANGULAR ! SOLVE U * X = Y ENDIF B(N) = B(N) / A(N,N) ! LAST ONE DONE FIRST IF( N .EQ. 1 ) RETURN ! NOW FROM THE BOTTOM UP, I = N-1, ... , 2, 1 DO II = 2,N I = N + 1 - II IP1 = I + 1 DO J = IP1,N B(I) = B(I) - A(I,J) * B(J) ENDDO ! LOOP ON J B(I) = B(I) / A(I,I) ENDDO ! LOOP ON II, I RETURN END SUBROUTINE GAUSSE ! *** end of filename gaucpp.f95 *********************