! 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
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)
! 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