! F08XPF Example Program Text ! Mark 23 Release. NAG Copyright 2011. 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 ! .. Implicit None Statement .. IMPLICIT NONE ! .. Function Return Value .. LOGICAL :: selctg ! .. Scalar Arguments .. COMPLEX (KIND=nag_wp), INTENT (IN) :: a, b ! .. Local Scalars .. LOGICAL :: d ! .. Intrinsic Functions .. INTRINSIC abs ! .. Executable Statements .. IF (abs(a)<6.0E0_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, f06uaf, nag_wp, x02ajf, x04dbf, zggesx 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, ifail, info, lda, ldb, ldvsl, & ldvsr, liwork, lwork, n, sdim ! .. Local Arrays .. COMPLEX (KIND=nag_wp), ALLOCATABLE :: a(:,:), alpha(:), b(:,:), & beta(:), 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(:) CHARACTER (1) :: clabs(1), rlabs(1) ! .. Intrinsic Functions .. 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 ldvsl = n ldvsr = n ALLOCATE (a(lda,n),alpha(n),b(ldb,n),beta(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 anorm = f06uaf('Frobenius',n,n,a,lda,rwork) bnorm = f06uaf('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 WRITE (nout,99999) & 'Number of eigenvalues for which SELCTG is true = ', sdim, & '(dimension of deflating subspaces)' 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 out the factors of the generalized Schur factorization ! ifail: behaviour on error exit ! =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft ifail = 0 CALL x04dbf('General',' ',n,n,a,lda,'Bracketed','F7.2', & 'Generalized Schur matrix S','Integer',rlabs,'Integer',clabs,80, & 0,ifail) WRITE (nout,*) FLUSH (nout) CALL x04dbf('General',' ',n,n,b,ldb,'Bracketed','F7.2', & 'Generalized Schur matrix T','Integer',rlabs,'Integer',clabs,80, & 0,ifail) WRITE (nout,*) FLUSH (nout) CALL x04dbf('General',' ',n,n,vsl,ldvsl,'Bracketed','F7.4', & 'Matrix of left generalized Schur vectors','Integer',rlabs, & 'Integer',clabs,80,0,ifail) WRITE (nout,*) FLUSH (nout) CALL x04dbf('General',' ',n,n,vsr,ldvsr,'Bracketed','F7.4', & 'Matrix of right generalized Schur vectors','Integer',rlabs, & 'Integer',clabs,80,0,ifail) ! Print out the reciprocal condition numbers WRITE (nout,*) WRITE (nout,99997) & 'Reciprocals of left and right projection norms onto', & 'the deflating subspaces for the selected eigenvalues', & 'RCONDE(1) = ', rconde(1), ', RCONDE(2) = ', rconde(2) WRITE (nout,*) WRITE (nout,99997) & 'Reciprocal condition numbers for the left and right', & 'deflating subspaces', 'RCONDV(1) = ', rcondv(1), & ', RCONDV(2) = ', rcondv(2) FLUSH (nout) ! Compute the machine precision and sqrt(anorm**2+bnorm**2) eps = x02ajf() abnorm = f06bnf(anorm,bnorm) tol = eps*abnorm ! Print out the approximate asymptotic error bound on the ! average absolute error of the selected eigenvalues given by ! eps*norm((A, B))/PL, where PL = RCONDE(1) WRITE (nout,*) WRITE (nout,99996) & 'Approximate asymptotic error bound for selected ', & 'eigenvalues = ', tol/rconde(1) ! Print out an approximate asymptotic bound on the maximum ! angular error in the computed deflating subspaces given by ! eps*norm((A, B))/DIF(2), where DIF(2) = RCONDV(2) WRITE (nout,99996) & 'Approximate asymptotic error bound for the deflating ', & 'subspaces = ', 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,2(A,1P,E8.1)) 99996 FORMAT (1X,2A,1P,E8.1) END PROGRAM f08xpfe