```!   F08XPF Example Program Text
!   Mark 25 Release. NAG Copyright 2014.

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
!     .. Accessibility Statements ..
Private
Public                               :: selctg
!     .. Parameters ..
Integer, Parameter, Public           :: 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
!       .. Intrinsic Procedures ..
Intrinsic                            :: abs
!       .. Executable Statements ..
selctg = (abs(a)<6.0_nag_wp*abs(b))
Return
End Function selctg
End Module f08xpfe_mod
Program f08xpfe

!     F08XPF Example Main Program

!     .. Use Statements ..
Use nag_library, Only: f06bnf, nag_wp, x02ajf, x04dbf, zgemm, zggesx,    &
zlange => f06uaf
Use f08xpfe_mod, Only: nb, nin, nout, selctg
!     .. Implicit None Statement ..
Implicit None
!     .. Local Scalars ..
Complex (Kind=nag_wp)                :: alph, bet
Real (Kind=nag_wp)                   :: abnorm, anorm, bnorm, eps,       &
normd, norme, tol
Integer                              :: i, ifail, info, 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(:)
Character (1)                        :: clabs(1), rlabs(1)
!     .. Intrinsic Procedures ..
Intrinsic                            :: cmplx, max, nint, real
!     .. Executable Statements ..
Write (nout,*) 'F08XPF Example Program Results'
Write (nout,*)
Flush (nout)
!     Skip heading in data file
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

!     Copy A and B into D and E respectively
d(1:n,1:n) = a(1:n,1:n)
e(1:n,1:n) = b(1:n,1:n)

!     Print matrices A and B
!     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','F8.4','Matrix A', &
'Integer',rlabs,'Integer',clabs,80,0,ifail)
Write (nout,*)
Flush (nout)

ifail = 0
Call x04dbf('General',' ',n,n,b,ldb,'Bracketed','F8.4','Matrix B', &
'Integer',rlabs,'Integer',clabs,80,0,ifail)
Write (nout,*)
Flush (nout)

!     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 .Or. info==(n+2)) Then

!       Compute A - Q*S*Z^H from the factorization of (A,B) and store in
!       matrix D
!       The NAG name equivelent of zgemm is f06zaf
alph = cmplx(1,kind=nag_wp)
bet = cmplx(0,kind=nag_wp)
Call zgemm('N','N',n,n,n,alph,vsl,ldvsl,a,lda,bet,c,ldc)
alph = cmplx(-1,kind=nag_wp)
bet = cmplx(1,kind=nag_wp)
Call zgemm('N','C',n,n,n,alph,c,ldc,vsr,ldvsr,bet,d,ldd)

!       Compute B - Q*T*Z^H from the factorization of (A,B) and store in
!       matrix E
alph = cmplx(1,kind=nag_wp)
bet = cmplx(0,kind=nag_wp)
Call zgemm('N','N',n,n,n,alph,vsl,ldvsl,b,ldb,bet,c,ldc)
alph = cmplx(-1,kind=nag_wp)
bet = cmplx(1,kind=nag_wp)
Call zgemm('N','C',n,n,n,alph,c,ldc,vsr,ldvsr,bet,e,lde)

!       Find norms of matrices D and E and warn if either is too large
normd = zlange('O',ldd,n,d,ldd,rwork)
norme = zlange('O',lde,n,e,lde,rwork)
If (normd>x02ajf()**0.75_nag_wp .Or. norme>x02ajf()**0.75_nag_wp) Then
Write (nout,*) 'Norm of A-(Q*S*Z^T) or norm of B-(Q*T*Z^T) &
&is much greater than 0.'
Write (nout,*) 'Schur factorization has failed.'
Else

!         Print solution
Write (nout,99999) &
'Number of eigenvalues for which SELCTG is true = ', sdim, &
'(dimension of deflating subspaces)'

Write (nout,*)
!         Print generalized eigenvalues
Write (nout,*) 'Selected generalized eigenvalues'

Do i = 1, sdim
If (beta(i)/=0.0_nag_wp) Then
Write (nout,99998) i, alpha(i)/beta(i)
Else
Write (nout,99997) i
End If
End Do

If (info==(n+2)) Then
Write (nout,99996) '***Note that rounding errors mean ', &
'that leading eigenvalues in the generalized', &
'Schur form no longer satisfy SELCTG = .TRUE.'
Write (nout,*)
End If
Flush (nout)

!         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,99995) &
'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,99995) &
'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

Else
Write (nout,99999) 'Failure in ZGGESX. INFO =', info
End If

99999 Format (1X,A,I4/1X,A)
99998 Format (1X,I2,1X,'(',1P,E11.4,',',E11.4,')')
99997 Format (1X,I4,'Eigenvalue is infinite')
99996 Format (1X,2A/1X,A)
99995 Format (1X,A/1X,A/1X,3(A,1P,E8.1))
End Program f08xpfe
```