Program f02jqfe

!     F02JQF Example Program Text

!     Mark 25 Release. NAG Copyright 2014.

!     .. Use Statements ..
      Use nag_library, Only: f02jqf, nag_wp, x02alf
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Real (Kind=nag_wp), Parameter    :: one = 1.0E+0_nag_wp
      Real (Kind=nag_wp), Parameter    :: tol = 0.0E0_nag_wp
      Real (Kind=nag_wp), Parameter    :: zero = 0.0E+0_nag_wp
      Integer, Parameter               :: nin = 5, nout = 6
!     .. Local Scalars ..
      Real (Kind=nag_wp)               :: inf, rbetaj, t0, t1, tmp
      Integer                          :: i, ifail, iwarn, j, lda, ldb, ldc,   &
                                          ldvl, ldvr, n, scal, sense, tdvl, tdvr
      Character (1)                    :: jobvl, jobvr
!     .. Local Arrays ..
      Complex (Kind=nag_wp), Allocatable :: a(:,:), alpha(:), b(:,:), beta(:), &
                                            c(:,:), cvr(:), e(:), vl(:,:),     &
                                            vr(:,:)
      Real (Kind=nag_wp), Allocatable  :: bevl(:), bevr(:), s(:)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: abs, aimag, cmplx, maxval, real
!     .. Executable Statements ..
      Write (nout,*) 'F02JQF Example Program Results'
!     Skip heading in data file and read in n, scal, sense, jobVL and jobVR
      Read (nin,*)
      Read (nin,*) n, scal, sense
      Read (nin,*) jobvl, jobvr

      lda = n
      ldb = n
      ldc = n
      ldvl = n
      ldvr = n
      tdvl = 2*n
      tdvr = 2*n
      Allocate (a(lda,n),b(ldb,n),c(ldc,n),alpha(2*n),beta(2*n),e(2*n), &
        vl(ldvl,tdvl),vr(ldvr,tdvr),s(2*n),bevr(2*n),bevl(2*n),cvr(n))

!     Read in the matrices A, B and C

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

!     Solve the quadratic eigenvalue problem
      ifail = -1
      Call f02jqf(scal,jobvl,jobvr,sense,tol,n,a,lda,b,ldb,c,ldc,alpha,beta, &
        vl,ldvl,vr,ldvr,s,bevl,bevr,iwarn,ifail)
      If (iwarn/=0) Then
        Write (nout,*)
        Write (nout,99999) 'Warning from f02jqf. IWARN =', iwarn
      End If

      Write (nout,*)
      If (ifail/=0) Then
        Write (nout,99999) 'Failure in f02jqf. IFAIL =', ifail
      Else
!       Infinity
        inf = x02alf()
        Do j = 1, 2*n
          rbetaj = real(beta(j))
          If (rbetaj>=one) Then
            e(j) = alpha(j)/rbetaj
          Else
            tmp = inf*rbetaj
            If ((abs(real(alpha(j)))<tmp) .And. (abs(aimag(alpha(j)))<tmp)) &
              Then
              e(j) = alpha(j)/rbetaj
            Else
              e(j) = cmplx(inf,zero,kind=nag_wp)
            End If
          End If
          If (abs(e(j))<inf) Then
            Write (nout,99998) 'Eigenvalue(', j, ') = (', real(e(j)), ', ', &
              aimag(e(j)), ')'
          Else
            Write (nout,99997) 'Eigenvalue(', j, ') is infinite'
          End If
        End Do

        If (jobvr=='V' .Or. jobvr=='v') Then
          Write (nout,*)
          Write (nout,99996) 'Right Eigenvectors'
          Do j = 1, 2*n, 3
            If (j<=2*n-2) Then
              Write (nout,99999) 'Eigenvector ', j, '          Eigenvector ', &
                j + 1, '          Eigenvector ', j + 2
              Write (nout,99995)(vr(i,j),vr(i,j+1),vr(i,j+2),i=1,n)
            Else If (j==2*n-1) Then
              Write (nout,99999) 'Eigenvector ', j, '          Eigenvector ', &
                j + 1
              Write (nout,99994)(vr(i,j),vr(i,j+1),i=1,n)
            Else If (j==2*n) Then
              Write (nout,99999) 'Eigenvector ', j
              Write (nout,99993)(vr(i,j),i=1,n)
            End If
          End Do
        End If

        If (jobvl=='V' .Or. jobvl=='v') Then
          Write (nout,*)
          Write (nout,*) 'Left Eigenvectors'
          Do j = 1, 2*n, 3
            If (j<=2*n-2) Then
              Write (nout,99999) 'Eigenvector ', j, '          Eigenvector ', &
                j + 1, '          Eigenvector ', j + 2
              Write (nout,99995)(vl(i,j),vl(i,j+1),vl(i,j+2),i=1,n)
            Else If (j==2*n-1) Then
              Write (nout,99999) 'Eigenvector ', j, '          Eigenvector ', &
                j + 1
              Write (nout,99994)(vl(i,j),vl(i,j+1),i=1,n)
            Else If (j==2*n) Then
              Write (nout,99999) 'Eigenvector ', j
              Write (nout,99993)(vl(i,j),i=1,n)
            End If
          End Do
        End If

        If (sense==1 .Or. sense==5 .Or. sense==6 .Or. sense==7) Then
          Write (nout,*)
          Write (nout,*) 'Eigenvalue Condition numbers'
          Do j = 1, 2*n
            Write (nout,99992) s(j)
          End Do
        End If

        If ((sense==3) .Or. (sense==4) .Or. (sense==6) .Or. (sense==7)) Then
          t0 = maxval(bevr)
          Write (nout,*)
          Write (nout,99996) &
            'Max backward error for eigenvalues and right eigenvectors', t0
        End If

        If ((sense==2) .Or. (sense==4) .Or. (sense==5) .Or. (sense==7)) Then
          t1 = maxval(bevl)
          Write (nout,*)
          Write (nout,99996) &
            'Max backward error for eigenvalues and left eigenvectors ', t1
        End If
      End If

99999 Format (1X,3(A,I4))
99998 Format (1X,A,I3,A,1P,E11.4,A,1P,E11.4,A)
99997 Format (1X,A,I3,A)
99996 Format (1X,A,1P,E11.4)
99995 Format (3(1X,'(',1P,E11.4,',',1P,E11.4,')'))
99994 Format (2(1X,'(',1P,E11.4,',',1P,E11.4,')'))
99993 Format (1X,'(',1P,E11.4,',',1P,E11.4,')')
99992 Format (1X,1P,E11.4)
    End Program f02jqfe