! F08XPF Example Program Text ! Mark 24 Release. NAG Copyright 2012. Module f08xpfe_mod ! F08XPF Example Program Module: ! Parameters and User-defined Routines ! .. Use Statements .. Use nag_library, Only: nag_wp ! .. Implicit None Statement .. Implicit None ! .. Parameters .. Integer, Parameter :: nb = 64, nin = 5, nout = 6 Contains Function selctg(a,b) ! Logical function selctg for use with ZGGESX (F08XPF) ! Returns the value .TRUE. if the absolute value of the eigenvalue ! a/b < 6.0 ! .. Function Return Value .. Logical :: selctg ! .. Scalar Arguments .. Complex (Kind=nag_wp), Intent (In) :: a, b ! .. Local Scalars .. Logical :: d ! .. Intrinsic Procedures .. Intrinsic :: abs ! .. Executable Statements .. If (abs(a)<6.0_nag_wp*abs(b)) Then d = .True. Else d = .False. End If selctg = d Return End Function selctg End Module f08xpfe_mod Program f08xpfe ! F08XPF Example Main Program ! .. Use Statements .. Use nag_library, Only: f06bnf, nag_wp, x02ajf, zggesx, zlange => f06uaf Use f08xpfe_mod, Only: nb, nin, nout, selctg ! .. Implicit None Statement .. Implicit None ! .. Local Scalars .. Real (Kind=nag_wp) :: abnorm, anorm, bnorm, eps, tol Integer :: i, info, j, lda, ldb, ldc, ldd, & lde, ldvsl, ldvsr, liwork, & lwork, n, sdim ! .. Local Arrays .. Complex (Kind=nag_wp), Allocatable :: a(:,:), alpha(:), b(:,:), & beta(:), c(:,:), d(:,:), e(:,:), & vsl(:,:), vsr(:,:), work(:) Complex (Kind=nag_wp) :: dummy(1) Real (Kind=nag_wp) :: rconde(2), rcondv(2) Real (Kind=nag_wp), Allocatable :: rwork(:) Integer :: idum(1) Integer, Allocatable :: iwork(:) Logical, Allocatable :: bwork(:) ! .. Intrinsic Procedures .. Intrinsic :: max, nint, real ! .. Executable Statements .. Write (nout,*) 'F08XPF Example Program Results' Write (nout,*) Flush (nout) ! Skip heading in data file Read (nin,*) Read (nin,*) n lda = n ldb = n ldc = n ldd = n lde = n ldvsl = n ldvsr = n Allocate (a(lda,n),alpha(n),b(ldb,n),beta(n),c(ldc,n),d(ldd,n),e(lde,n), & vsl(ldvsl,n),vsr(ldvsr,n),rwork(8*n),bwork(n)) ! Use routine workspace query to get optimal workspace. lwork = -1 liwork = -1 ! The NAG name equivalent of zggesx is f08xpf Call zggesx('Vectors (left)','Vectors (right)','Sort',selctg, & 'Both reciprocal condition numbers',n,a,lda,b,ldb,sdim,alpha,beta,vsl, & ldvsl,vsr,ldvsr,rconde,rcondv,dummy,lwork,rwork,idum,liwork,bwork, & info) ! Make sure that there is enough workspace for blocksize nb. lwork = max(n*nb+n*n/2,nint(real(dummy(1)))) liwork = max(n+2,idum(1)) Allocate (work(lwork),iwork(liwork)) ! Read in the matrices A and B Read (nin,*)(a(i,1:n),i=1,n) Read (nin,*)(b(i,1:n),i=1,n) ! Find the Frobenius norms of A and B ! The NAG name equivalent of the LAPACK auxiliary zlange is f06uaf anorm = zlange('Frobenius',n,n,a,lda,rwork) bnorm = zlange('Frobenius',n,n,b,ldb,rwork) ! Find the generalized Schur form ! The NAG name equivalent of zggesx is f08xpf Call zggesx('Vectors (left)','Vectors (right)','Sort',selctg, & 'Both reciprocal condition numbers',n,a,lda,b,ldb,sdim,alpha,beta,vsl, & ldvsl,vsr,ldvsr,rconde,rcondv,work,lwork,rwork,iwork,liwork,bwork, & info) If (info>0 .And. info/=(n+2)) Then Write (nout,99999) 'Failure in ZGGESX. INFO =', info Else ! Print Results Write (nout,99999) 'Number of selected eigenvalues = ', sdim Write (nout,*) If (info==(n+2)) Then Write (nout,99998) '***Note that rounding errors mean ', & 'that leading eigenvalues in the generalized', & 'Schur form no longer satisfy SELCTG = .TRUE.' Write (nout,*) End If Flush (nout) ! Print information on generalized eigenvalues. Write (nout,*) 'Selected Generalized Eigenvalues' Write (nout,*) Write (nout,99996)(j,alpha(j)/beta(j),j=1,sdim) ! Compute the machine precision and sqrt(anorm**2+bnorm**2) eps = x02ajf() abnorm = f06bnf(anorm,bnorm) tol = eps*abnorm ! Print out the reciprocal condition numbers and error bound for ! selected eigenvalues Write (nout,*) Write (nout,99997) & 'Reciprocal condition numbers for the average of the', & 'selected eigenvalues and their asymptotic error bound', & 'rcond-left = ', rconde(1), ', rcond-right = ', rconde(2), & ', error = ', tol/rconde(1) Write (nout,*) Write (nout,99997) & 'Reciprocal condition numbers for the deflating subspaces', & 'and their approximate asymptotic error bound', 'rcond-left = ', & rcondv(1), ', rcond-right = ', rcondv(2), ', error = ', & tol/rcondv(2) End If 99999 Format (1X,A,I4/1X,A) 99998 Format (1X,2A/1X,A) 99997 Format (1X,A/1X,A/1X,3(A,1P,E8.1)) 99996 Format (1X,I2,1X,'(',1P,E11.4,',',E11.4,')') End Program f08xpfe