! F08PBF Example Program Text ! Mark 23 Release. NAG Copyright 2011. MODULE f08pbfe_mod ! F08PBF 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 select(wr,wi) ! Logical function select for use with DGEESX (F08PBF) ! Returns the value .TRUE. if the eigenvalue is real ! .. Implicit None Statement .. IMPLICIT NONE ! .. Function Return Value .. LOGICAL :: select ! .. Scalar Arguments .. REAL (KIND=nag_wp), INTENT (IN) :: wi, wr ! .. Local Scalars .. LOGICAL :: d ! .. Executable Statements .. IF (wi==0.0E0_nag_wp) THEN d = .TRUE. ELSE d = .FALSE. END IF select = d RETURN END FUNCTION select END MODULE f08pbfe_mod PROGRAM f08pbfe ! F08PBF Example Main Program ! .. Use Statements .. USE nag_library, ONLY : dgeesx, dgemm, f06raf, nag_wp, x02ajf, x04caf USE f08pbfe_mod, ONLY : nb, nin, nout, select ! .. Implicit None Statement .. IMPLICIT NONE ! .. Local Scalars .. REAL (KIND=nag_wp) :: alpha, anorm, beta, eps, rconde, & rcondv, tol INTEGER :: i, ifail, info, lda, ldc, ldd, & ldvs, liwork, lwork, n, sdim ! .. Local Arrays .. REAL (KIND=nag_wp), ALLOCATABLE :: a(:,:), c(:,:), d(:,:), vs(:,:), & wi(:), work(:), wr(:) REAL (KIND=nag_wp) :: dummy(1) INTEGER :: idum(1) INTEGER, ALLOCATABLE :: iwork(:) LOGICAL, ALLOCATABLE :: bwork(:) ! .. Intrinsic Functions .. INTRINSIC max, nint ! .. Executable Statements .. WRITE (nout,*) 'F08PBF Example Program Results' WRITE (nout,*) FLUSH (nout) ! Skip heading in data file READ (nin,*) READ (nin,*) n lda = n ldc = n ldd = n ldvs = n ALLOCATE (a(lda,n),c(ldc,n),d(ldd,n),vs(ldvs,n),wi(n),wr(n),bwork(n)) ! Use routine workspace query to get optimal workspace. lwork = -1 liwork = -1 ! The NAG name equivalent of dgeesx is f08pbf CALL dgeesx('Vectors (Schur)','Sort',select, & 'Both reciprocal condition numbers',n,a,lda,sdim,wr,wi,vs,ldvs, & rconde,rcondv,dummy,lwork,idum,liwork,bwork,info) ! Make sure that there is enough workspace for blocksize nb. liwork = max((n*n)/4,idum(1)) lwork = max(n*(nb+2+n/2),nint(dummy(1))) ALLOCATE (work(lwork),iwork(liwork)) ! Read in the matrix A READ (nin,*) (a(i,1:n),i=1,n) ! Print Matrix A ! 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) ! Find the Frobenius norms of A anorm = f06raf('Frobenius',n,n,a,lda,work) ! Find the Schur factorization of A ! The NAG name equivalent of dgeesx is f08pbf CALL dgeesx('Vectors (Schur)','Sort',select, & 'Both reciprocal condition numbers',n,a,lda,sdim,wr,wi,vs,ldvs, & rconde,rcondv,work,lwork,iwork,liwork,bwork,info) IF (info==0 .OR. info==(n+2)) THEN ! Compute Z*T*Z**T from Schur factorization of A alpha = 1._nag_wp beta = 0._nag_wp CALL dgemm('N','N',n,n,n,alpha,vs,ldvs,a,lda,beta,c,ldc) CALL dgemm('N','T',n,n,n,alpha,c,ldc,vs,ldvs,beta,d,ldd) ! Print matrix A, as computed from Schur factorization ifail = 0 CALL x04caf('General',' ',n,n,d,ldd, & 'Matrix A recreated from Schur factorization Z*T*Z^T',ifail) WRITE (nout,*) FLUSH (nout) ! Print solution WRITE (nout,99999) & 'Number of eigenvalues for which SELECT is true = ', sdim, & '(dimension of invariant subspace)' WRITE (nout,*) IF (info==(n+2)) THEN WRITE (nout,99998) '***Note that rounding errors mean ', & 'that leading eigenvalues in the Schur form', & 'no longer satisfy SELECT = .TRUE.' WRITE (nout,*) END IF FLUSH (nout) ! Print out the reciprocal condition numbers WRITE (nout,99997) & 'Reciprocal of projection norm onto the invariant', & 'subspace for the selected eigenvalues', 'RCONDE = ', rconde WRITE (nout,*) WRITE (nout,99996) & 'Reciprocal condition number for the invariant subspace', & 'RCONDV = ', rcondv ! Compute the machine precision eps = x02ajf() tol = eps*anorm ! Print out the approximate asymptotic error bound on the ! average absolute error of the selected eigenvalues given by ! eps*norm(A)/RCONDE WRITE (nout,*) WRITE (nout,99995) & 'Approximate asymptotic error bound for selected ', & 'eigenvalues = ', tol/rconde ! Print out an approximate asymptotic bound on the maximum ! angular error in the computed invariant subspace given by ! eps*norm(A)/RCONDV WRITE (nout,99995) & 'Approximate asymptotic error bound for the invariant ', & 'subspace = ', tol/rcondv ELSE WRITE (nout,99994) 'Failure in DGEESX. INFO =', info END IF 99999 FORMAT (1X,A,I4/1X,A) 99998 FORMAT (1X,2A/1X,A) 99997 FORMAT (1X,A/1X,A/1X,A,1P,E8.1) 99996 FORMAT (1X,A/1X,A,1P,E8.1) 99995 FORMAT (1X,2A,1P,E8.1) 99994 FORMAT (1X,A,I4) END PROGRAM f08pbfe