Program f12jufe

!     F12JUF Example Program Text

!     Mark 30.3 Release. nAG Copyright 2024.

!     .. Use Statements ..
      Use, Intrinsic                   :: iso_c_binding, Only: c_ptr
      Use nag_library, Only: f12jaf, f12jbf, f12jff, f12juf, f12jzf, nag_wp,   &
                             x04daf, zsymm, zsytrf, zsytrs
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Complex (Kind=nag_wp), Parameter :: cone = (1.0_nag_wp,0.0_nag_wp)
      Complex (Kind=nag_wp), Parameter :: czero = (0.0_nag_wp,0.0_nag_wp)
      Integer, Parameter               :: nin = 5, nout = 6
!     .. Local Scalars ..
      Type (c_ptr)                     :: handle
      Complex (Kind=nag_wp)            :: emid, tmp, ze
      Real (Kind=nag_wp)               :: eps, r
      Integer                          :: deg, i, ifail, irevcm, iter, j, k,   &
                                          lda, ldx, ldy, ldz, lwork, m0, n,    &
                                          nconv
!     .. Local Arrays ..
      Complex (Kind=nag_wp), Allocatable :: a(:,:,:), az(:,:), d(:), work(:),  &
                                          x(:,:), y(:,:), z(:,:)
      Real (Kind=nag_wp), Allocatable  :: resid(:)
      Integer, Allocatable             :: ipiv(:)
!     .. Executable Statements ..
      Write (nout,*) 'F12JUF Example Program Results'
      Write (nout,*)
!     Skip heading in data file
      Read (nin,*)
      Read (nin,*) n, deg

      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,deg+1),  &
        az(n,n),ipiv(n),work(lwork))

!     Read the upper triangular part of the matrices from data file
      Do j = 1, deg + 1
        Read (nin,*)(a(i,i:n,j),i=1,n)
      End Do

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

!     Set the ellipse minor/major axis ratio
      Call f12jbf(handle,'Ellipse Contour Ratio = 80',ifail)

      emid = (-2.0_nag_wp,-2.0_nag_wp)
      r = 1.5_nag_wp
!     Generate the contour
      Call f12jff(handle,emid,r,ifail)

      irevcm = 0
revcomm: Do
        Call f12juf(handle,irevcm,deg,ze,n,x,ldx,y,ldy,k,m0,nconv,d,z,ldz,eps, &
          iter,resid,ifail)
        Select Case (irevcm)
        Case (1)
!         Form the matrix \sum ze^i A_i
          Do j = 1, n
            az(1:j,j) = a(1:j,j,1)
          End Do
          Do i = 1, deg
            Do j = 1, n
              az(1:j,j) = az(1:j,j) + (ze**i)*a(1:j,j,i+1)
            End Do
          End Do
!         Compute a Bunch-Kaufman factorization of \sum ze^i A_i
!         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 (\sum ze^i A_i)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 <- A_k z
          Call zsymm('Left','Upper',n,m0,cone,a(1,1,k+1),lda,z,ldz,czero,x,    &
            ldx)
          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: unit first element
        Do i = 1, nconv
          tmp = z(1,i)
          z(1:n,i) = z(1:n,i)/tmp
        End Do

        Call x04daf('General',' ',n,nconv,z,ldz,'Eigenvectors',ifail)
      Else
        Write (nout,99998) 'Failure in f12juf IFAIL =', ifail
      End If

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

99999 Format ((3X,4(' (',F7.4,',',F7.4,')',:)))
99998 Format (1X,A,I4)
    End Program f12jufe