Program f12jsfe

!     F12JSF 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, f12jsf, 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, ze
      Real (Kind=nag_wp)               :: eps, r
      Integer                          :: i, ifail, irevcm, iter, lda, ldb,    &
                                          ldx, ldy, ldz, lwork, m0, n, nconv
!     .. Local Arrays ..
      Complex (Kind=nag_wp), Allocatable :: a(:,:), az(:,:), b(:,:), d(:),     &
                                          work(:), x(:,:), y(:,:), z(:,:)
      Real (Kind=nag_wp), Allocatable  :: resid(:)
      Integer, Allocatable             :: ipiv(:)
!     .. Executable Statements ..
      Write (nout,*) 'F12JSF Example Program Results'
      Write (nout,*)
!     Skip heading in data file
      Read (nin,*)
      Read (nin,*) n

      lda = n
      ldb = 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),        &
        b(ldb,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)

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

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

!     Set the ellipse rotation angle and minor/major axis ratio
      Call f12jbf(handle,'Ellipse Contour Ratio = 10',ifail)
      Call f12jbf(handle,'Ellipse Rotation Angle = 45',ifail)

      emid = (0.0_nag_wp,0.0_nag_wp)
      r = 1.0_nag_wp
!     Generate the contour
      Call f12jff(handle,emid,r,ifail)

      irevcm = 0
revcomm: Do
        Call f12jsf(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 B - A
          Do i = 1, n
            az(1:i,i) = ze*b(1:i,i) - cone*a(1:i,i)
          End Do
!         Compute a Bunch-Kaufman factorization of ze B -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
!         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
          Call zsymm('Left','Upper',n,m0,cone,a,lda,z,ldz,czero,x,ldx)
          Cycle revcomm
        Case (4)
!         Compute x <- Bz
          Call zsymm('Left','Upper',n,m0,cone,b,ldb,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)

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

      Else
        Write (nout,99998) 'Failure in f12jsf 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 f12jsfe