! F08XAF Example Program Text ! Mark 24 Release. NAG Copyright 2012. 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 ! .. 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, dlange => f06raf, nag_wp, x02ajf, & x04caf Use f08xafe_mod, Only: nb, nin, nout, selctg ! .. Implicit None Statement .. Implicit None ! .. Local Scalars .. Real (Kind=nag_wp) :: alph, bet, normd, norme 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 Procedures .. 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) ! 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 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 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 warn if either is too large ! f06raf is the NAG name equivalent of the LAPACK auxiliary dlange 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 Write (nout,99999) & 'Number of eigenvalues for which SELCTG is true = ', sdim 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 End If Flush (nout) End If 99999 Format (1X,A,I4) 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) End Program f08xafe