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

NAG FL Interface Introduction
Example description
    Program f12jtfe

!     F12JTF Example Program Text

!     Mark 28.4 Release. NAG Copyright 2022.

!      Intrinsic Procedures ..
!     .. Use Statements ..
      Use, Intrinsic                   :: iso_c_binding, Only: c_ptr
      Use nag_library, Only: f11dnf, f11dqf, f11xnf, f11znf, f12jaf, f12jbf,   &
                             f12jgf, f12jtf, 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)            :: ze
      Real (Kind=nag_wp)               :: eps, rnorm, tol
      Integer                          :: ccn, 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 :: a(:), az(:), azh(:), d(:), w(:),   &
                                          work(:), x(:,:), y(:,:), z(:,:),     &
                                          zedge(:)
      Real (Kind=nag_wp), Allocatable  :: resid(:)
      Integer, Allocatable             :: icol(:), icolz(:), icolzh(:),        &
                                          idiag(:), idiagh(:), ipiv(:),        &
                                          ipivp(:), ipivph(:), ipivq(:),       &
                                          ipivqh(:), irow(:), irowz(:),        &
                                          irowzh(:), istr(:), istrh(:),        &
                                          iwork(:), nedge(:), tedge(:)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: conjg, max, min, sqrt
!     .. Executable Statements ..
      Write (nout,*) 'F12JTF 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,2*m0),y(ldy,m0),z(ldz,2*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=handle,ifail=ifail)

!     Set the integration type
      Call f12jbf(handle=handle,str='Integration Type = Trapezoidal',          &
        ifail=ifail)

!     Define arrays for the custom contour
      ccn = 4
      Allocate (nedge(ccn),tedge(ccn),zedge(ccn))

      zedge(1) = (0.0_nag_wp,1.0_nag_wp)
      zedge(2) = (0.0_nag_wp,-1.0_nag_wp)
      zedge(3) = (-1.0_nag_wp,-1.0_nag_wp)
      zedge(4) = (-1.0_nag_wp,1.0_nag_wp)

      tedge(1) = 50
      tedge(2) = 0
      tedge(3) = 0
      tedge(4) = 0

      nedge(1) = 20
      nedge(2) = 1
      nedge(3) = 1
      nedge(4) = 1

!     Generate the contour
      Call f12jgf(handle=handle,ccn=ccn,nedge=nedge,tedge=tedge,zedge=zedge,   &
        ifail=ifail)

      irevcm = 0
revcomm: Do
        Call f12jtf(handle=handle,irevcm=irevcm,ze=ze,n=n,x=x,ldx=ldx,y=y,     &
          ldy=ldy,m0=m0,nconv=nconv,d=d,z=z,ldz=ldz,eps=eps,iter=iter,         &
          resid=resid,ifail=ifail)
        Select Case (irevcm)
        Case (1)
!         Form the sparse matrix ze I - A
          nnzz = nnz + n
          az(1:nnz) = -a(1:nnz)
          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)y = w, 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) = -conjg(a(1:nnz))
          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 y = w, 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 f11xnf('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^H z
          Do i = 1, m0
            Call f11xnf('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 .Or. ifail==14) Then

!       Print solution

        Write (nout,*) 'Eigenvalues'
        Write (nout,99999) d(1:nconv)
        Flush (nout)

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

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

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

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