Program f08wcfe

!     F08WCF Example Program Text

!     Mark 26.1 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'
      Flush (nout)
!     Skip heading in data file
      Read (nin,*)
      Read (nin,*) n
      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

      Read (nin,*)(a(i,1:n),i=1,n)
      Read (nin,*)(b(i,1:n),i=1,n)

!     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,*)
        Flush (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,*)
        Flush (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