! F08XBF Example Program Text ! Mark 23 Release. NAG Copyright 2011. 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 ! .. Implicit None Statement .. IMPLICIT NONE ! .. 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, f06bnf, f06raf, 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, 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 Functions .. 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) ! 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 anorm = f06raf('Frobenius',n,n,a,lda,work) bnorm = f06raf('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 Q * S * Z**T from Schur factorization of (A,B) alph = 1._nag_wp bet = 0._nag_wp CALL dgemm('N','N',n,n,n,alph,vsl,ldvsl,a,lda,bet,c,ldc) CALL dgemm('N','T',n,n,n,alph,c,ldc,vsr,ldvsr,bet,d,ldd) ! Compute Q * T * Z**T from Schur factorization of (A,B) CALL dgemm('N','N',n,n,n,alph,vsl,ldvsl,b,ldb,bet,c,ldc) CALL dgemm('N','T',n,n,n,alph,c,ldc,vsr,ldvsr,bet,e,lde) ! Print matrix A and B, as computed from Schur factorization ifail = 0 CALL x04caf('General',' ',n,n,d,ldd, & 'Matrix A recreated from Schur factorization Q*S*Z^T',ifail) WRITE (nout,*) CALL x04caf('General',' ',n,n,e,lde, & 'Matrix B recreated from Schur factorization Q*T*Z^T',ifail) WRITE (nout,*) 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 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 f08xbfe