! F08XBF Example Program Text ! Mark 24 Release. NAG Copyright 2012. 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 ! .. Parameters .. Integer, Parameter :: 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 ! .. Function Return Value .. Logical :: selctg ! .. Scalar Arguments .. Real (Kind=nag_wp), Intent (In) :: ai, ar, b ! .. Local Scalars .. Logical :: d ! .. Executable Statements .. If (ai==0.0E0_nag_wp) Then d = .True. Else d = .False. End If selctg = d Return End Function selctg End Module f08xbfe_mod Program f08xbfe ! F08XBF Example Main Program ! .. Use Statements .. Use nag_library, Only: dgemm, dggesx, dlange => f06raf, f06bnf, nag_wp, & x02ajf, x04caf Use f08xbfe_mod, Only: nb, nin, nout, selctg ! .. 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 blocksize 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 Matrix A and Matrix 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,*) Call x04caf('General',' ',n,n,b,ldb,'Matrix B',ifail) Write (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 .And. info/=(n+2)) Then Write (nout,99999) 'Failure in DGGESX. INFO =', info Else ! Compute A - Q * S * Z^T from Schur factorization of (A,B) and store in ! matrix D ! The NAG name equivelent 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 Schur 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 print warning 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.8_nag_wp .Or. norme>x02ajf()**0.8_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 Results Write (nout,99999) & 'Number of eigenvalues for which SELCTG is true = ', sdim, & '(dimension of deflating subspaces)' Write (nout,*) ! Print the 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 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