Program f08ngfe ! F08NGF Example Program Text ! Mark 24 Release. NAG Copyright 2012. ! .. Use Statements .. Use nag_library, Only: dgehrd, dhsein, dhseqr, dormhr, nag_wp, x04caf ! .. Implicit None Statement .. Implicit None ! .. Parameters .. Integer, Parameter :: nin = 5, nout = 6 ! .. Local Scalars .. Complex (Kind=nag_wp) :: eig, eig1 Real (Kind=nag_wp) :: thresh Integer :: i, ifail, info, j, k, lda, ldc, ldh, & ldvl, ldz, lwork, m, n ! .. Local Arrays .. Real (Kind=nag_wp), Allocatable :: a(:,:), c(:,:), h(:,:), tau(:), & vl(:,:), wi(:), work(:), wr(:), z(:,:) Integer, Allocatable :: ifaill(:), ifailr(:) Logical, Allocatable :: select(:) ! .. Intrinsic Procedures .. Intrinsic :: aimag, cmplx, real ! .. Executable Statements .. Write (nout,*) 'F08NGF Example Program Results' ! Skip heading in data file Read (nin,*) Read (nin,*) n ldz = 1 lda = n ldc = n ldh = n ldvl = n lwork = 64*n Allocate (a(lda,n),c(ldc,n),h(ldh,n),tau(n),vl(ldvl,n),wi(n), & work(lwork),wr(n),z(ldz,1),ifaill(n),ifailr(n),select(n)) ! Read A from data file Read (nin,*)(a(i,1:n),i=1,n) Read (nin,*) thresh ! Reduce A to upper Hessenberg form H = (Q**T)*A*Q ! The NAG name equivalent of dgehrd is f08nef Call dgehrd(n,1,n,a,lda,tau,work,lwork,info) ! Copy A to H h(1:n,1:n) = a(1:n,1:n) ! Calculate the eigenvalues of H (same as A) ! The NAG name equivalent of dhseqr is f08pef Call dhseqr('Eigenvalues','No vectors',n,1,n,h,ldh,wr,wi,z,ldz,work, & lwork,info) Write (nout,*) If (info>0) Then Write (nout,*) 'Failure to converge.' Else Write (nout,*) 'Eigenvalues' Write (nout,99999)(' (',wr(i),',',wi(i),')',i=1,n) Do i = 1, n select(i) = wr(i) < thresh End Do ! Calculate the eigenvectors of H (as specified by SELECT), ! storing the result in C ! The NAG name equivalent of dhsein is f08pkf Call dhsein('Right','QR','No initial vectors',select,n,a,lda,wr,wi,vl, & ldvl,c,ldc,n,m,work,ifaill,ifailr,info) ! Calculate the eigenvectors of A = Q * (eigenvectors of H) ! The NAG name equivalent of dormhr is f08ngf Call dormhr('Left','No transpose',n,m,1,n,a,lda,tau,c,ldc,work,lwork, & info) ! Print eigenvectors Write (nout,*) Flush (nout) ! Normalize selected eigenvectors j = 0 k = 1 Do While (k<=n) If (select(k)) Then j = j + 1 If (wi(k)==0.0_nag_wp) Then Do i = 2, n c(i,j) = c(i,j)/c(1,j) End Do c(1,j) = 1.0_nag_wp Else eig1 = cmplx(c(1,j),c(1,j+1),kind=nag_wp) c(1,j) = 1.0_nag_wp c(1,j+1) = 0.0_nag_wp Do i = 2, n eig = cmplx(c(i,j),c(i,j+1),kind=nag_wp) eig = eig/eig1 c(i,j) = real(eig) c(i,j+1) = aimag(eig) End Do j = j + 1 k = k + 1 End If End If k = k + 1 End Do ! ifail: behaviour on error exit ! =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft ifail = 0 Call x04caf('General',' ',n,m,c,ldc,'Contents of array C',ifail) End If 99999 Format (1X,A,F8.4,A,F8.4,A) End Program f08ngfe