! Last change: DOS 24 Jul 2000 8:45 pm ! *** copyright 2000 *** ! *** filename chex35.f95 *** John F. Monahan ** ! ********************** program chex35 ! compute condition numbers for Hilbert matrices and compare ! with results from gaucpp ! implicit none INTEGER, parameter :: length = 10 INTEGER, DIMENSION(length,length) :: ib INTEGER, DIMENSION(length) :: pi integer i,j,n,ijb INTEGER ibicof REAL, DIMENSION(length,length) :: b real hn,bn,s,rcond1,rcond2 ! 20 format(/'Condition of Hilbert matrix of size',i4/'inverse'/) 21 FORMAT(8i10) 22 format(/'Infinity norms of matrix and inverse',2e16.6) 23 FORMAT('reciprocal condition number and estimate',2e16.6) ! 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 END interface ! open ( UNIT=6, FILE='chex35.out' ) ! initialize harmonic sum hn = 0. ! do matrices in some small sizes do n = 1,7 ! infinity norm is harmonic sum - not perfect here hn = hn + 1./real(n) ! write(6,20) n ! create inverse of Hilbert matrix -- this one ! is integer and can be constructed exactly ! ** constructing Hilbert matrix entails some ! roundoff error immediately - unfair ** do i = 1,n do j = 1,n ijb = j*ibicof(i+j-2,i-1)*ibicof(i+n-1,i-1) & & *ibicof(j+n-1,n-i)*ibicof(n,j) if( mod(i+j,2) .eq. 1 ) ijb = - ijb ib(i,j) = ijb b(i,j) = real( ib(i,j) ) end do ! loop on j write(6,21) (ib(i,j),j=1,n) end do ! loop on i ! compute norm of b ! bn = 0. do j = 1,n s = 0. do i = 1,n s = s + abs( b(j,i) ) end do ! loop on i bn = max( s, bn ) end do ! loop on j ! write out norms of Hilbert matrix & inverse write(6,22) hn,bn ! reciprocal condition number rcond1 = 1./(hn*bn) ! compute condition estimator from gaucpp call gaucpp(b,n,pi,rcond2) write(6,23) rcond1,rcond2 end do ! loop on n stop END program chex35 INTEGER FUNCTION IBICOF(I,J) ! COMPUTES BINOMIAL COEFFICIENTS, USUALLY BY TABLE ! INTEGER, INTENT(IN) :: I,J INTEGER K,L,N INTEGER, DIMENSION(21) :: PUSH ! TABLE OF POINTERS INTEGER, DIMENSION(122) :: BICFTB ! ! TEST WHETHER IT IS IN TABLE - ! REMEMBER THE RESULT WILL GO HAYWIRE IN INTEGERS IF TOO BIG ! DATA PUSH/1,2,4,6,8,11,14,18,22,27,32,38,44,51,58,66,74,83,92, & & 102,112/ DATA BICFTB/1,1,1,1,2,1,3,1,4,6,1,5,10,1,6,15,20,1,7,21,35, & & 1,8,28,56,70,1,9,36,84,126,1,10,45,120,210,252, & & 1,11,55,165,330,462,1,12,66,220,495,792,924, & & 1,13,78,286,715,1287,1716,1,14,91,364,1001,2002,3003,3432,& & 1,15,105,455,1365,3003,5005,6435,1,16,120,560,1820,4368,8008, & & 11440,12870, 1,17,136,680,2380,6188,12376,19448,24310,1,18,153,& & 816,3060,8568,18564,31824,43758,48620, 1,19,171,969,3876,11628,& & 27132,50388,75582,92378, 1,20,190,1140,4845,15504,38760,77520, & & 125970, 167960, 184756 / ! ! TAKE CARE OF SPECIAL CASES IBICOF = 0 IF( ( J .LT. 0 ) .OR. ( J .GT. I ) ) RETURN IBICOF = 1 IF( ( J .EQ. 0 ) .OR. ( J .EQ. I ) ) RETURN ! GENERAL CASE L = MIN(J,I-J) IF( I .LE. 20 ) THEN ! TABLE LOOKUP K = PUSH(I+1) ! FIND START POINT FOR I IBICOF = BICFTB(K+L) ! PULL NUMBER FROM TABLE ELSE ! MULTIPLY ! FOR BIG VALUES COMPUTE BY MULTIPLICATION IBICOF = 1 DO N = 1,L IBICOF = (IBICOF*(I-N+1)) / N END DO ! LOOP ON N END IF ! ( I .LE. 20 ) RETURN END FUNCTION IBICOF 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 chex35.f95 *********************