! Last change: DOS 24 Jul 2000 8:41 pm ! *** copyright 2000 *** ! *** filename cdnxpmt.f95 *** John F. Monahan ** ! ********************** program cdnxpmt ! cdnxpmt test program for simulation study of matrix condition ! estimator -- mimicking Stewart's but with infinity norm ! implicit none integer, parameter :: nmax = 50 ! max dimension integer, parameter :: nrep = 10 ! number of replications real, dimension(nmax,nmax) :: a,ai real, dimension(nmax) :: xu,xv real ran,gnroul real c,s,su,sv,zu,zv,rkap,r,anrm,ainrm,rax,rcond integer i,j,k,kn,n,lkap,meth,krep integer, dimension(nmax) :: pi integer, dimension(4) :: nlst ! ! 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 END INTERFACE ! !* 21 format(5f12.6) !* 22 format(10f8.4) 23 format(4i4,f12.6) ! ! four sizes of matrices data nlst/ 5, 10, 25, 50/ ! open( unit=6, file='cdnxpmt.out' ) ! ! initialize uniform random number generator r = ran(5151917) ! ! ** this is a three factor experiment ** ! - method of conditioning ! - magnitude of condition ! - size of matrix ! ! two methods ! method =1 is sharp break, method =2 is expo decay do meth = 1,2 ! four sizes of condition, in powers of 10 do lkap = 1,4 ! four sizes of matrices from list do kn = 1,4 n = nlst(kn) ! !w** write(6,23) n,lkap,meth,nrep rkap = 10.**lkap ! r=1 is sharp break, other is exponential decay r = 1. if( meth .eq. 2 ) r = exp( - log(rkap)/(n-1) ) ! ! start monte carlo loop do krep = 1,nrep ! ! create initial diagonal matrix a = 0. ai = 0. c = 1. s = 1. ! scale the diagonal do i = 1,n a(i,i) = c ai(i,i) = s c = c * r s = s / r end do ! loop on i ! this is real different for meth = 1 a(n,n) = 1./rkap ai(n,n) = rkap ! diagonal matrix now done, pre and post mult to get a and ! its inverse do i = 2,n ! generate Householder matrix with this size su = 0. sv = 0. do k = 1,i xu(k) = gnroul(k) su = su + xu(k)*xu(k) xv(k) = gnroul(-k) ! beware optimizing compilers sv = sv + xv(k)*xv(k) end do ! loop on k ! choose sign in usual way -- does not affect ! infinity norm computed below su = sign( sqrt(su), xu(1) ) sv = sign( sqrt(sv), xv(1) ) xu(1) = xu(1) + su su = su*xu(1) xv(1) = xv(1) + sv sv = sv*xv(1) ! ready to multiply Householder tfms ! first U on left of a, V on left of ai do j = 1,n zu = 0. zv = 0. do k = 1,i zu = zu + xu(k)*a(n-i+k,j) zv = zv + xv(k)*ai(n-i+k,j) end do ! loop on k zu = zu / su zv = zv / sv do k = 1,i a(n-i+k,j) = a(n-i+k,j) - zu*xu(k) ai(n-i+k,j) = ai(n-i+k,j) - zv*xv(k) end do ! loop on k end do ! loop on j -- columns of a and ai ! now U on right of ai, V on right of a do j = 1,n zu = 0. zv = 0. do k = 1,i zu = zu + xu(k)*ai(j,n-i+k) zv = zv + xv(k)*a(j,n-i+k) end do ! loop on k zu = zu / su zv = zv / sv do k = 1,i ai(j,n-i+k) = ai(j,n-i+k) - zu*xu(k) a(j,n-i+k) = a(j,n-i+k) - zv*xv(k) end do ! loop on k end do ! loop on j -- rows of a and ai ! end of iteration i = size of Householder tfm end do ! loop on i ! check out matrices -- is ai inverse of a? !* aia = matmul(a,ai) !* do i = 1,n !* write(6,22) (aia(i,j),j= 1,n) !* end do ! write out matrices !w do i = 1,n !w write(6,22) (a(i,j),j=1,n),(ai(i,j),j=1,n) !w end do ! loop on i ! ! matrix norms below unaffected by row/col sign changes ! ! compute matrix norms for a and ai anrm = pimnrm(a,n) ainrm = pimnrm(ai,n) ! compute condition number rax = anrm*ainrm !w** write(6,21) anrm,ainrm,rax ! ! get estimate of reciprocal condition from gaucpp call gaucpp(a,n,pi,rcond) !w** write(6,21) anrm,ainrm,rax,rcond ! compute ratio of true / estimate c = (rax*rcond) write(6,23) n,lkap,meth,krep,c end do ! loop on krep (replications) ! end do ! loop on kn -- chooses n = size of matrix end do ! loop on lkap -- size of condition end do ! loop on meth(od) stop contains real function pimnrm(a,n) ! computes the p=infinity norm of an n by n matrix a ! a real matrix, dimensioned 50 by 50 for convenience only ! n integer, true dimension of matrix real, dimension(nmax,nmax) :: a integer n real s integer i,j ! pimnrm = 0. do j = 1,n s = 0. do i = 1,n s = s + abs( a(j,i) ) end do ! loop on i pimnrm = max( s, pimnrm ) end do ! loop on i return end function pimnrm end program cdnxpmt REAL FUNCTION GNROUL(IX) ! RATIO OF UNIFORMS ALGORITHM FOR GENERATING NORMAL(0,1) RV'S ! USING LEVA'S QUADRATIC INNER AND OUTER BOUNDS ! ! A J KINDERMAN AND J F MONAHAN (1977) "COMPUTER GENERATION OF RANDOM ! VARIABLES USING THE RATIO OF UNIFORM DEVIATES," ACM TRANSACTIONS ON ! MATHEMATICAL SOFTWARE, VOLUME 3, PP.257-260 ! ! J L LEVA (1992) "A FAST NORMAL RANDOM NUMBER GENERATOR," ACM ! TRANSACTIONS ON MATHEMATICAL SOFTWARE, VOLUME 18, PP. 449-453. ! IMPLICIT NONE INTEGER, INTENT(IN) :: IX ! DUMMY ARGUMENT REAL, PARAMETER :: R = 1.7155277 ! FIRST CONSTANT R = SQRT(2/E) ! CONSTANTS S,T CENTER OF ELLIPSES REAL, PARAMETER :: S = .449871 REAL, PARAMETER :: T = -.386595 ! SHAPE PARAMETERS OF ELLIPSES REAL, PARAMETER :: A = .19600 REAL, PARAMETER :: B = .25472 ! RADII OF ELLIPSES REAL, PARAMETER :: R1 = .27597 REAL, PARAMETER :: R2 = .27846 REAL RAN REAL U,V,X,Y,Q ! START / RESTART HERE DO ! NOTE UNRESTRICTED DO ! GENERATE (U,V) UNIFORMLY IN BOX U = RAN(1) V = R*(RAN(2) - 0.5) X = U - S Y = ABS(V) - T ! COMPUTE QUADRATIC PIECE Q = X*X + Y*(A*Y - B*X) ! COMPUTE DEVIATE BEFORE TESTS GNROUL = V / U ! INNER BOUND -- QUICK ACCEPT CHECK IF( Q .LE. R1 ) RETURN ! OUTER BOUND -- QUICK REJECT CHECK IF( A .GT. R2 ) CYCLE ! FINAL RATIO OF UNIFORMS CHECK IF( GNROUL*GNROUL .LE. -4.*LOG(U) ) RETURN ! REJECT -- START OVER END DO ! END OF UNRESTRICTED DO END FUNCTION GNROUL REAL FUNCTION RAN(IDUM) ! PORTABLE IMPLEMENTATION OF UNIFORM PSEUDORANDOM NUMBER GENERATOR ! LEWIS, GOODMAN, & MILLER MULTIPLICATIVE GENERATOR ! X(N+1) = MOD( 16807*X(N), 2**31-1 ) ! ! P. BRANTLEY, B.L. FOX, L. SCHRAGE (1983) A GUIDE TO SIMULATION ! SPRINGER-VERLAG, NEW YORK. PP. 200-202. ! UPDATED VERSION OF ! LINUS SCHRAGE (1979)'A MORE PORTABLE FORTRAN RANDOM NUMBER GENERATOR' ! ACM TRANSACTIONS ON MATHEMATICAL SOFTWARE, VOLUME 5, PP. 132-138. ! ! ARGUMENT ! IDUM INTEGER FIRST CALL SETS SEED, IGNORED IN SUBSEQUENT CALLS ! IMPLICIT NONE INTEGER, INTENT(IN) :: IDUM REAL, PARAMETER :: RP = 4.656612875E-10 ! 1/P INTEGER, PARAMETER :: A = 16807 ! MULTIPLIER INTEGER, PARAMETER :: B = 127773 ! B = P / A INTEGER, PARAMETER :: C = 2836 ! C = P MOD A INTEGER, PARAMETER :: P = 2147483647 ! MODULUS 2**31 - 1 INTEGER, SAVE :: IX = 0 INTEGER K1 ! ! IF NOT FIRST CALL, THEN SKIP SETTING SEED IF( IX .EQ. 0) IX = IDUM ! WRITE NUMBER AS ALPHA*2**16 + BETA K1 = IX / B IX = A*( IX - K1*B) - K1*C ! ABOVE DOES A*IX MOD B -K1*C IF( IX .LT. 0 ) IX = IX + P ! RP BELOW IS RECIPROCAL OF P RAN = REAL(IX)*RP RETURN END FUNCTION RAN 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 cdnxpmt.f95 *********************