NAG Library Manual, Mark 30.1
Interfaces:  FL   CL   CPP   AD 

NAG FL Interface Introduction
Example description
    Program f12jkfe

!     F12JKF Example Program Text

!     Mark 30.1 Release. NAG Copyright 2024.

!      Intrinsic Procedures ..
!     .. Use Statements ..
      Use, Intrinsic                   :: iso_c_binding, Only: c_ptr
      Use nag_library, Only: f11dnf, f11dqf, f11xaf, f11znf, f12jaf, f12jbf,   &
                             f12jff, f12jkf, f12jzf, nag_wp, x02ajf, x04daf
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Complex (Kind=nag_wp), Parameter :: cone = (1.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, rnorm, tol
      Integer                          :: i, ifail, irevcm, iter, itn, la,     &
                                          ldx, ldy, ldz, liwork, lwork, m, m0, &
                                          n, nconv, nnz, nnzc, nnzch, nnzz,    &
                                          nnzzh, npivm
!     .. Local Arrays ..
      Complex (Kind=nag_wp), Allocatable :: az(:), azh(:), d(:), w(:),         &
                                          work(:), y(:,:)
      Real (Kind=nag_wp), Allocatable  :: a(:), resid(:), x(:,:), z(:,:)
      Integer, Allocatable             :: icol(:), icolz(:), icolzh(:),        &
                                          idiag(:), idiagh(:), ipiv(:),        &
                                          ipivp(:), ipivph(:), ipivq(:),       &
                                          ipivqh(:), irow(:), irowz(:),        &
                                          irowzh(:), istr(:), istrh(:),        &
                                          iwork(:)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: cmplx, conjg, max, min, sqrt
!     .. Executable Statements ..
      Write (nout,*) 'F12JKF Example Program Results'
      Write (nout,*)
!     Skip heading in data file
      Read (nin,*)
      Read (nin,*) n
      Read (nin,*) nnz

      m0 = n
      ldx = n
      ldy = n
      ldz = n
      m = min(n,50)
      lwork = max(64*n,4*n+m*(m+n+5)+121)
      la = 2*(nnz+n)
      liwork = 7*n + 2
      tol = sqrt(x02ajf())

      Allocate (x(ldx,m0),y(ldy,2*m0),z(ldz,m0),resid(2*m0),d(m0),a(nnz),w(n), &
        icol(nnz),irow(nnz),az(la),azh(la),icolz(la),irowz(la),icolzh(la),     &
        irowzh(la),ipiv(n),iwork(liwork),idiag(n),ipivp(n),ipivq(n),ipivph(n), &
        ipivqh(n),istrh(n+1),idiagh(n),istr(n+1),work(lwork))

!     Read the matrix A from data file
      Do i = 1, nnz
        Read (nin,*) a(i), irow(i), icol(i)
      End Do

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

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

      emid = (2.5_nag_wp,1.2_nag_wp)
      r = 1.4_nag_wp
!     Generate the contour
      Call f12jff(handle,emid,r,ifail)

      irevcm = 0
revcomm: Do
        Call f12jkf(handle,irevcm,ze,n,x,ldx,y,ldy,m0,nconv,d,z,ldz,eps,iter,  &
          resid,ifail)
        Select Case (irevcm)
        Case (1)
!         Form the sparse matrix ze I - A
          nnzz = nnz + n
          az(1:nnz) = cmplx(-a(1:nnz),kind=nag_wp)
          irowz(1:nnz) = irow(1:nnz)
          icolz(1:nnz) = icol(1:nnz)
          Do i = 1, n
            irowz(nnz+i) = i
            icolz(nnz+i) = i
            az(nnz+i) = ze
          End Do
!         Sort the elements into correct coordinate storage format
          Call f11znf(n,nnzz,az,irowz,icolz,'S','R',istr,iwork,ifail)
!         Form incomplete LU factorization of ze I - A
          Call f11dnf(n,nnzz,az,la,irowz,icolz,0,0.0_nag_wp,'P','N',ipivp,     &
            ipivq,istr,idiag,nnzc,npivm,iwork,liwork,ifail)
          If (ifail/=0) Then
            Exit revcomm
          End If
          Cycle revcomm
        Case (2)
!         Solve the linear system (ze I - A)w = y, with m0 righthand sides
          Do i = 1, m0
!           We will need to store the final result in y
            w(1:n) = y(1:n,i)
!           Initial guess
            y(1:n,i) = cone
            Call f11dqf('RGMRES',n,nnzz,az,la,irowz,icolz,ipivp,ipivq,istr,    &
              idiag,w,m,tol,500,y(1,i),rnorm,itn,work,lwork,ifail)
          End Do
          If (ifail/=0) Then
            Exit revcomm
          End If
          Cycle revcomm
        Case (3)
!         Form the sparse matrix (ze I - A)^H
          nnzzh = nnz + n
          azh(1:nnz) = cmplx(-a(1:nnz),kind=nag_wp)
          irowzh(1:nnz) = icol(1:nnz)
          icolzh(1:nnz) = irow(1:nnz)
          Do i = 1, n
            irowzh(nnz+i) = i
            icolzh(nnz+i) = i
            azh(nnz+i) = conjg(ze)
          End Do
!         Sort the elements into correct coordinate storage format
          Call f11znf(n,nnzzh,azh,irowzh,icolzh,'S','R',istrh,iwork,ifail)
!         Form incomplete LU factorization of (ze I - A)^H
          Call f11dnf(n,nnzzh,azh,la,irowzh,icolzh,0,0.0_nag_wp,'P','N',       &
            ipivph,ipivqh,istrh,idiagh,nnzch,npivm,iwork,liwork,ifail)
          If (ifail/=0) Then
            Exit revcomm
          End If
          Cycle revcomm
        Case (4)
!         Solve the linear system (ze I - A)^H w = y, with m0 righthand sides
          Do i = 1, m0
!           We will need to store the final result in y
            w(1:n) = y(1:n,i)
!           Initial guess
            y(1:n,i) = cone
            Call f11dqf('RGMRES',n,nnzzh,azh,la,irowzh,icolzh,ipivph,ipivqh,   &
              istrh,idiagh,w,m,tol,500,y(1,i),rnorm,itn,work,lwork,ifail)
          End Do
          If (ifail/=0) Then
            Exit revcomm
          End If
          Cycle revcomm
        Case (5)
!         Compute x <- Az
          Do i = 1, m0
            Call f11xaf('N',n,nnz,a,irow,icol,'N',z(1,i),x(1,i),ifail)
          End Do
          If (ifail/=0) Then
            Exit revcomm
          End If
          Cycle revcomm
        Case (6)
!         Compute x <- A^T z
          Do i = 1, m0
            Call f11xaf('T',n,nnz,a,irow,icol,'N',z(1,i),x(1,i),ifail)
          End Do
          If (ifail/=0) Then
            Exit revcomm
          End If
          Cycle revcomm
        Case (7)
!         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 (8)
!         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)

        Call x04daf('General',' ',n,nconv,y,ldz,'Right eigenvectors',ifail)
        Call x04daf('General',' ',n,nconv,y(1,m0+1),ldz,'Left eigenvectors',   &
          ifail)

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