! F08XBF Example Program Text
! Mark 29.2 Release. NAG Copyright 2023.
Module f08xbfe_mod
! F08XBF 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(ar,ai,b)
! Logical function selctg for use with DGGESX (F08XBF)
! Returns the value .TRUE. if the eigenvalue is real and positive
! .. Function Return Value ..
Logical :: selctg
! .. Scalar Arguments ..
Real (Kind=nag_wp), Intent (In) :: ai, ar, b
! .. Executable Statements ..
selctg = (ar>0._nag_wp .And. ai==0._nag_wp .And. b/=0._nag_wp)
Return
End Function selctg
End Module f08xbfe_mod
Program f08xbfe
! F08XBF Example Main Program
! .. Use Statements ..
Use f08xbfe_mod, Only: nb, nin, nout, selctg
Use nag_library, Only: dgemm, dggesx, dlange => f06raf, f06bnf, nag_wp, &
x02ajf, x04caf
! .. Implicit None Statement ..
Implicit None
! .. Local Scalars ..
Real (Kind=nag_wp) :: abnorm, alph, anorm, bet, bnorm, &
eps, normd, norme, tol
Integer :: i, ifail, info, lda, ldb, ldc, ldd, &
lde, ldvsl, ldvsr, liwork, lwork, n, &
sdim
! .. Local Arrays ..
Real (Kind=nag_wp), Allocatable :: a(:,:), alphai(:), alphar(:), &
b(:,:), beta(:), c(:,:), d(:,:), &
e(:,:), vsl(:,:), vsr(:,:), work(:)
Real (Kind=nag_wp) :: rconde(2), rcondv(2), rdum(1)
Integer :: idum(1)
Integer, Allocatable :: iwork(:)
Logical, Allocatable :: bwork(:)
! .. Intrinsic Procedures ..
Intrinsic :: max, nint
! .. Executable Statements ..
Write (nout,*) 'F08XBF 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),alphai(n),alphar(n),b(ldb,n),beta(n),vsl(ldvsl,n), &
vsr(ldvsr,n),bwork(n),c(ldc,n),d(ldd,n),e(lde,n))
! Use routine workspace query to get optimal workspace.
lwork = -1
liwork = -1
! The NAG name equivalent of dggesx is f08xbf
Call dggesx('Vectors (left)','Vectors (right)','Sort',selctg, &
'Both reciprocal condition numbers',n,a,lda,b,ldb,sdim,alphar,alphai, &
beta,vsl,ldvsl,vsr,ldvsr,rconde,rcondv,rdum,lwork,idum,liwork,bwork, &
info)
! Make sure that there is enough workspace for block size nb.
lwork = max(8*(n+1)+16+n*nb+n*n/2,nint(rdum(1)))
liwork = max(n+6,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 x04caf('General',' ',n,n,a,lda,'Matrix A',ifail)
Write (nout,*)
Flush (nout)
ifail = 0
Call x04caf('General',' ',n,n,b,ldb,'Matrix B',ifail)
Write (nout,*)
Flush (nout)
! Find the Frobenius norms of A and B
! The NAG name equivalent of the LAPACK auxiliary dlange is f06raf
anorm = dlange('Frobenius',n,n,a,lda,work)
bnorm = dlange('Frobenius',n,n,b,ldb,work)
! Find the generalized Schur form
! The NAG name equivalent of dggesx is f08xbf
Call dggesx('Vectors (left)','Vectors (right)','Sort',selctg, &
'Both reciprocal condition numbers',n,a,lda,b,ldb,sdim,alphar,alphai, &
beta,vsl,ldvsl,vsr,ldvsr,rconde,rcondv,work,lwork,iwork,liwork,bwork, &
info)
If (info==0 .Or. info==(n+2)) Then
! Compute A - Q*S*Z^T from the factorization of (A,B) and store in
! matrix D
! The NAG name equivalent of dgemm is f06yaf
alph = 1.0_nag_wp
bet = 0.0_nag_wp
Call dgemm('N','N',n,n,n,alph,vsl,ldvsl,a,lda,bet,c,ldc)
alph = -1.0_nag_wp
bet = 1.0_nag_wp
Call dgemm('N','T',n,n,n,alph,c,ldc,vsr,ldvsr,bet,d,ldd)
! Compute B - Q*T*Z^T from the factorization of (A,B) and store in
! matrix E
alph = 1.0_nag_wp
bet = 0.0_nag_wp
Call dgemm('N','N',n,n,n,alph,vsl,ldvsl,b,ldb,bet,c,ldc)
alph = -1.0_nag_wp
bet = 1.0_nag_wp
Call dgemm('N','T',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 = dlange('O',ldd,n,d,ldd,work)
norme = dlange('O',lde,n,e,lde,work)
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, '(', alphar(i)/beta(i), ',', &
alphai(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)
! Print out the reciprocal condition numbers
Write (nout,*)
Write (nout,99995) &
'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,99995) &
'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,99994) &
'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,99994) &
'Approximate asymptotic error bound for the deflating ', &
'subspaces = ', tol/rcondv(2)
End If
Else
Write (nout,99999) 'Failure in DGGESX. INFO =', info
End If
99999 Format (1X,A,I4,/,1X,A)
99998 Format (1X,I4,5X,A,F7.3,A,F7.3,A)
99997 Format (1X,I4,'Eigenvalue is infinite')
99996 Format (1X,2A,/,1X,A)
99995 Format (1X,A,/,1X,A,/,1X,2(A,1P,E8.1))
99994 Format (1X,2A,1P,E8.1)
End Program f08xbfe