! Last change: DOS 24 Jul 2000 8:47 pm ! *** copyright 2000 *** ! *** filename chex36.f95 *** John F. Monahan ** ! ********************** program chex36 ! chex36 demonstration of condition number of a matrix ! example is interpolation problem -- just amplification ! integer, parameter :: msize = 10 ! declared size matrix, vector real, dimension(msize,msize) :: a real, dimension(msize) :: b,bh,bd real xi,eps,rcond,xnorm,dnorm real, parameter :: pie = 3.14159265 integer, dimension(msize) :: pi integer i,j,k,n,method ! 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 ! 20 format(' n = ',i2,' method: equal spacing',4x,' rcond =',f9.6) 21 format(' n = ',i2,' method: Cheby spacing',4x,' rcond =',f9.6) 22 format(' || x || =',f12.6,' eps =',f12.6,4x,'||x-x*|| =',f12.6) ! OPEN( UNIT=6, FILE='chex36.out' ) ! ! set n, initial eps, vary method n = 4 do method = 1,2 eps = 1.e-5 ! ! set up vandermonde matrix for interpolation do i = 1,n if(method .eq. 1) then ! two choices of x's -- first is equally spaced on (-1,1) xi = real(2*i-1-n)/real(2*n) ! second choice is chebyshev abscissas else xi = cos( real(2*i-1)*pie/real(2*n) ) endif ! function interpolated is simple-looking b(i) = 1./(1.+xi*xi) ! make a copy bh(i) = b(i) a(i,1) = 1. do j = 2,n a(i,j) = xi*a(i,j-1) enddo ! loop on j enddo ! loop on i ! have matrix -- factor and give condition estimate call gaucpp(a,n,pi,rcond) ! solve equations call gausse(a,n,pi,b) if( method .eq. 1 ) then write(6,20) n,rcond else write(6,21) n,rcond end if ! ( method .eq. 1 ) ! compute norm of solution vector - norm of rhs near 1 xnorm = 0. do i = 1,n xnorm = max( xnorm, abs( b(i) ) ) enddo ! add disturbances on the order of eps do k = 1,5 do i = 1,n eps = - eps bd(i) = bh(i) + eps enddo ! loop on i ! solve equations with perturbed rhs call gausse(a,n,pi,bd) ! compute norm of difference on solution dnorm = 0. do i = 1,n dnorm = max( dnorm, abs( b(i) - bd(i) ) ) enddo ! loop on i write(6,22) xnorm,eps,dnorm ! make eps bigger for next run through eps = 4.*eps enddo ! loop on k -- 5 reps with increasing eps enddo ! loop on method ( equal / Chebyshev spacing ) stop end program chex36 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 ! 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 ! 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 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 ! 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 ! 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 chex36.f95 *********************