! F08XAF Example Program Text ! Mark 23 Release. NAG Copyright 2011. MODULE f08xafe_mod ! F08XAF 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 DGGES (F08XAF) ! 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 f08xafe_mod PROGRAM f08xafe ! F08XAF Example Main Program ! .. Use Statements .. USE nag_library, ONLY : dgemm, dgges, nag_wp, x04caf USE f08xafe_mod, ONLY : nb, nin, nout, selctg ! .. Implicit None Statement .. IMPLICIT NONE ! .. Local Scalars .. REAL (KIND=nag_wp) :: alph, bet INTEGER :: i, ifail, info, lda, ldb, ldc, & ldd, lde, ldvsl, ldvsr, 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) :: dummy(1) LOGICAL, ALLOCATABLE :: bwork(:) ! .. Intrinsic Functions .. INTRINSIC max, nint ! .. Executable Statements .. WRITE (nout,*) 'F08XAF 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 ! The NAG name equivalent of dgges is f08xaf CALL dgges('Vectors (left)','Vectors (right)','Sort',selctg,n,a,lda,b, & ldb,sdim,alphar,alphai,beta,vsl,ldvsl,vsr,ldvsr,dummy,lwork,bwork, & info) ! Make sure that there is enough workspace for blocksize nb. lwork = max(8*n+16+n*nb,nint(dummy(1))) ALLOCATE (work(lwork)) ! 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 generalized Schur form ! The NAG name equivalent of dgges is f08xaf CALL dgges('Vectors (left)','Vectors (right)','Sort',selctg,n,a,lda,b, & ldb,sdim,alphar,alphai,beta,vsl,ldvsl,vsr,ldvsr,work,lwork,bwork, & info) IF (info>0 .AND. info/=(n+2)) THEN WRITE (nout,99999) 'Failure in DGGES. 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 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) END IF 99999 FORMAT (1X,A,I4) 99998 FORMAT (1X,2A/1X,A) END PROGRAM f08xafe