Program f12jjfe

!     F12JJF Example Program Text

!     Mark 29.2 Release. NAG Copyright 2023.

!     .. Use Statements ..
      Use, Intrinsic                   :: iso_c_binding, Only: c_ptr
      Use nag_library, Only: dsymm, f12jaf, f12jbf, f12jef, f12jjf, f12jzf,    &
                             nag_wp, x04caf, zsytrf, zsytrs
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter               :: nin = 5, nout = 6
!     .. Local Scalars ..
      Type (c_ptr)                     :: handle
      Complex (Kind=nag_wp)            :: ze
      Real (Kind=nag_wp)               :: emax, emin, eps
      Integer                          :: i, ifail, irevcm, iter, lda, ldx,    &
                                          ldy, ldz, lwork, m0, n, nconv
!     .. Local Arrays ..
      Complex (Kind=nag_wp), Allocatable :: az(:,:), work(:), y(:,:)
      Real (Kind=nag_wp), Allocatable  :: a(:,:), d(:), resid(:), x(:,:),      &
                                          z(:,:)
      Integer, Allocatable             :: ipiv(:)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: cmplx
!     .. Executable Statements ..
      Write (nout,*) 'F12JJF Example Program Results'
      Write (nout,*)
!     Skip heading in data file
      Read (nin,*)
      Read (nin,*) n

      lda = n
      m0 = n
      ldx = n
      ldy = n
      ldz = n
      lwork = 64*n

      Allocate (x(ldx,m0),y(ldy,m0),z(ldz,m0),resid(m0),d(m0),a(lda,n),        &
        az(n,n),ipiv(n),work(lwork))

!     Read the upper triangular part of the matrix A from data file
      Read (nin,*)(a(i,i:n),i=1,n)

!     Initialize the data handle
      ifail = 0
      Call f12jaf(handle,ifail)

!     Set the number of contour points to use
      Call f12jbf(handle,'Contour Points Hermitian = 12',ifail)

      emin = -3.0_nag_wp
      emax = 3.0_nag_wp

!     Generate the contour
      Call f12jef(handle,emin,emax,ifail)

      irevcm = 0
revcomm: Do
        Call f12jjf(handle,irevcm,ze,n,x,ldx,y,ldy,m0,nconv,d,z,ldz,eps,iter,  &
          resid,ifail)
        Select Case (irevcm)
        Case (1)
!         Form the matrix ze I - A
          Do i = 1, n
            az(i,i:n) = cmplx(-a(i,i:n),kind=nag_wp)
            az(i,i) = ze + az(i,i)
          End Do
!         Compute a Bunch-Kaufman factorization of ze I -A
!         The NAG name equivalent of zsytrf is f07nrf
          Call zsytrf('Upper',n,az,n,ipiv,work,lwork,ifail)
          If (ifail/=0) Then
            Exit revcomm
          End If
          Cycle revcomm
        Case (2)
!         Solve the linear system (ze I -A)w = y, overwriting y with w
!         The NAG name equivalent of zsytrs is f07nsf
          Call zsytrs('Upper',n,m0,az,n,ipiv,y,ldy,ifail)
          If (ifail/=0) Then
            Exit revcomm
          End If
          Cycle revcomm
        Case (3)
!         Compute x <- Az
!         The NAG name equivalent of dsymm is f06ycf
          Call dsymm('Left','Upper',n,m0,1.0_nag_wp,a,lda,z,ldz,0.0_nag_wp,x,  &
            ldx)
          Cycle revcomm
        Case (4)
!         Since we are not solving a generalized eigenvalue problem set x = z
          x(1:n,1:m0) = z(1:n,1:m0)
          Cycle revcomm
        Case Default
          Exit revcomm
        End Select

      End Do revcomm

      If (ifail==0) Then

!       Print solution

        Write (nout,*) 'Eigenvalues'
        Write (nout,99999) d(1:nconv)
        Flush (nout)
!       Normalize the eigenvectors: first element positive
        Do i = 1, nconv
          If ((z(1,i)<0.0_nag_wp)) Then
            z(1:n,i) = -z(1:n,i)
          End If
        End Do

        Call x04caf('General',' ',n,nconv,z,ldz,'Eigenvectors',ifail)

      Else
        Write (nout,99998) 'Failure in f12jjf IFAIL =', ifail
      End If

!     Destroy the data handle
      Call f12jzf(handle,ifail)

99999 Format (3X,(8F8.4))
99998 Format (1X,A,I4)
    End Program f12jjfe