Program f08wnfe

!     F08WNF Example Program Text

!     Mark 26.2 Release. NAG Copyright 2017.

!     .. Use Statements ..
      Use nag_library, Only: m01daf, m01edf, nag_wp, x02ajf, x04daf, zggev
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Real (Kind=nag_wp), Parameter    :: one = 1.0_nag_wp
      Real (Kind=nag_wp), Parameter    :: zero = 0.0_nag_wp
      Integer, Parameter               :: nb = 64, nin = 5, nout = 6
      Complex (Kind=nag_wp), Parameter :: cone = (one,zero)
!     .. 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 :: a(:,:), alpha(:), b(:,:), beta(:), &
                                          vr(:,:), work(:)
      Complex (Kind=nag_wp)            :: dummy(1,1)
      Real (Kind=nag_wp), Allocatable  :: rwork(:)
      Integer, Allocatable             :: irank(:)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: abs, all, max, maxloc, nint, real
!     .. Executable Statements ..
      Write (nout,*) 'F08WNF 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),alpha(n),b(ldb,n),beta(n),vr(ldvr,n),rwork(8*n))

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

!     Make sure that there is enough workspace for block size nb.
      lwork = max((nb+1)*n,nint(real(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 zggev is f08wnf
      Call zggev('No left vectors','Vectors (right)',n,a,lda,b,ldb,alpha,beta, &
        dummy,1,vr,ldvr,work,lwork,rwork,info)

      If (info>0) Then
        Write (nout,*)
        Write (nout,99999) 'Failure in ZGGEV. INFO =', info
      Else
!       Re-normalize the eigenvectors, largest absolute element real (=1)
        Do i = 1, n
          rwork(1:n) = abs(vr(1:n,i))
          k = maxloc(rwork(1:n),1)
          scal = cone/vr(k,i)
          vr(1:n,i) = vr(1:n,i)*scal
          vr(k,i) = cone
        End Do

        Write (nout,*)
        Flush (nout)
        If (all(abs(beta(1:n))>x02ajf())) Then
!         Reorder eigenvalues by descending absolute value and print
          alpha(1:n) = alpha(1:n)/beta(1:n)
          rwork(1:n) = abs(alpha(1:n))
          Allocate (irank(n))
          ifail = 0
          Call m01daf(rwork,1,n,'Descending',irank,ifail)
          Call m01edf(alpha,1,n,irank,ifail)
          ifail = 0
          Call x04daf('Gen',' ',1,n,alpha,1,'Eigenvalues:',ifail)

!         Reorder eigenvectors accordingly
          Do j = 1, n
            beta(1:n) = vr(j,1:n)
            Call m01edf(beta,1,n,irank,ifail)
            vr(j,1:n) = beta(1:n)
          End Do
        Else
          Write (nout,*)                                                       &
            'Some of the eigenvalues are infinite or undetermined'
          Write (nout,*)
          Flush (nout)
          ifail = 0
          Call x04daf('Gen',' ',1,n,alpha,1,'Alpha:',ifail)
          Call x04daf('Gen',' ',1,n,beta,1,'Beta:',ifail)
        End If
        Write (nout,*)
        Flush (nout)
        ifail = 0
        Call x04daf('Gen',' ',n,n,vr,ldvr,'Eigenvectors (columns):',ifail)
      End If

99999 Format (1X,A,I4)
    End Program f08wnfe