```    Program f08wcfe

!     F08WCF Example Program Text

!     Mark 26 Release. NAG Copyright 2016.

!     .. Use Statements ..
Use nag_library, Only: dggev3, m01def, m01edf, nag_wp, x02ajf, x04caf,   &
x04daf
!     .. Implicit None Statement ..
Implicit None
!     .. Parameters ..
Real (Kind=nag_wp), Parameter    :: zero = 0.0_nag_wp
Integer, Parameter               :: nin = 5, nout = 6
!     .. Local Scalars ..
Complex (Kind=nag_wp)            :: scal
Integer                          :: i, ifail, info, j, k, lda, ldb,      &
ldvr, lwork, n
!     .. Local Arrays ..
Complex (Kind=nag_wp), Allocatable :: eigval(:), eigvec(:,:)
Real (Kind=nag_wp), Allocatable  :: a(:,:), alphai(:), alphar(:),        &
b(:,:), beta(:), vr(:,:), work(:)
Real (Kind=nag_wp)               :: dummy(1,1)
Integer, Allocatable             :: irank(:)
!     .. Intrinsic Procedures ..
Intrinsic                        :: abs, all, cmplx, conjg, maxloc, nint
!     .. Executable Statements ..
Write (nout,*) 'F08WCF Example Program Results'
!     Skip heading in data file
lda = n
ldb = n
ldvr = n
Allocate (a(lda,n),alphai(n),alphar(n),b(ldb,n),beta(n),vr(ldvr,n),      &
eigvec(n,n),eigval(n),irank(n))

!     Use routine workspace query to get optimal workspace.
lwork = -1
!     The NAG name equivalent of dggev3 is f08wcf
Call dggev3('No left vectors','Vectors (right)',n,a,lda,b,ldb,alphar,    &
alphai,beta,dummy,1,vr,ldvr,dummy,lwork,info)

lwork = nint(dummy(1,1))
Allocate (work(lwork))

!     Read in the matrices A and B

!     Solve the generalized eigenvalue problem

!     The NAG name equivalent of dggev3 is f08wcf
Call dggev3('No left vectors','Vectors (right)',n,a,lda,b,ldb,alphar,    &
alphai,beta,dummy,1,vr,ldvr,work,lwork,info)
If (info>0) Then
Write (nout,*)
Write (nout,99999) 'Failure in DGGEV3. INFO =', info
Go To 100
End If

!     Re-normalize the eigenvectors, largest absolute element real
j = 0
Do i = 1, n
If (alphai(i)==zero) Then
eigvec(1:n,i) = cmplx(vr(1:n,i),zero,kind=nag_wp)
Else If (j==0) Then
eigvec(1:n,i) = cmplx(vr(1:n,i),vr(1:n,i+1),kind=nag_wp)
j = 1
Else
eigvec(1:n,i) = cmplx(vr(1:n,i-1),-vr(1:n,i),kind=nag_wp)
j = 0
End If
work(1:n) = abs(eigvec(1:n,i))
k = maxloc(work(1:n),1)
scal = conjg(eigvec(k,i))/abs(eigvec(k,i))
eigvec(1:n,i) = eigvec(1:n,i)*scal
End Do

!     If eigenvalues are finite, order by descending absolute values
If (all(abs(beta(1:n))>x02ajf())) Then
!       add small amount to alphai to distinguish conjugates
alphai(1:n) = alphai(1:n) + x02ajf()*10.0_nag_wp
eigval(1:n) = cmplx(alphar(1:n),alphai(1:n),kind=nag_wp)
eigval(1:n) = eigval(1:n)/beta(1:n)
work(1:n) = abs(eigval(1:n))
ifail = 0
Call m01def(work,n,1,n,1,1,'Descending',irank,ifail)
Call m01edf(eigval,1,n,irank,ifail)

!       Print ordered eigenvalues
ifail = 0
Call x04daf('Gen',' ',1,n,eigval,1,'Eigenvalues:',ifail)

!       Order the eigenvectors in the same way and print
Do j = 1, n
eigval(1:n) = eigvec(j,1:n)
Call m01edf(eigval,1,n,irank,ifail)
eigvec(j,1:n) = eigval(1:n)
End Do

Write (nout,*)
ifail = 0
Call x04daf('Gen',' ',n,n,eigvec,n,'Right Eigenvectors (columns):',    &
ifail)
Else
Write (nout,*) 'Some of the eigenvalues are infinite'
Write (nout,*)
ifail = 0
Call x04caf('Gen',' ',1,n,alphar,1,'Alpha (real):',ifail)
Call x04caf('Gen',' ',1,n,alphai,1,'Alpha (imag):',ifail)
Call x04caf('Gen',' ',1,n,beta,1,'Beta:',ifail)
End If
100   Continue

99999 Format (1X,A,I4)
End Program f08wcfe
```