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

NAG FL Interface Introduction
Example description
    Program f12jrfe

!     F12JRF Example Program Text

!     Mark 28.5 Release. NAG Copyright 2022.

!      Intrinsic Procedures ..
!     .. Use Statements ..
      Use, Intrinsic                   :: iso_c_binding, Only: c_ptr
      Use nag_library, Only: f11dnf, f11dqf, f11xsf, f11znf, f12jaf, f12jbf,   &
                             f12jef, f12jrf, 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)               :: emax, emin, eps, rnorm, tol
      Integer                          :: i, ifail, irevcm, iter, itn, la,     &
                                          ldx, ldy, ldz, liwork, lwork, m, m0, &
                                          n, nconv, nnza, nnzb, nnzc, nnzch,   &
                                          nnzz, nnzzh, npivm
!     .. Local Arrays ..
      Complex (Kind=nag_wp), Allocatable :: a(:), az(:), azh(:), b(:), w(:),   &
                                          work(:), x(:,:), y(:,:), z(:,:)
      Real (Kind=nag_wp), Allocatable  :: d(:), resid(:)
      Integer, Allocatable             :: icola(:), icolb(:), icolz(:),        &
                                          icolzh(:), idiag(:), idiagh(:),      &
                                          ipiv(:), ipivp(:), ipivph(:),        &
                                          ipivq(:), ipivqh(:), irowa(:),       &
                                          irowb(:), irowz(:), irowzh(:),       &
                                          istr(:), istrh(:), iwork(:)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: conjg, max, min, sqrt
!     .. Executable Statements ..
      Write (nout,*) 'F12JRF Example Program Results'
      Write (nout,*)
!     Skip heading in data file
      Read (nin,*)
      Read (nin,*) n
      Read (nin,*) nnza, nnzb

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

      Allocate (x(ldx,m0),y(ldy,m0),z(ldz,m0),resid(m0),d(m0),a(nnza),b(nnzb), &
        w(n),icola(nnza),irowa(nnza),icolb(nnzb),irowb(nnzb),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, nnza
        Read (nin,*) a(i), irowa(i), icola(i)
      End Do

!     Read the matrix B from data file
      Do i = 1, nnzb
        Read (nin,*) b(i), irowb(i), icolb(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 = Zolotarev',            &
        ifail=ifail)

      emin = -1.0_nag_wp
      emax = 0.0_nag_wp
!     Generate the contour
      Call f12jef(handle=handle,emin=emin,emax=emax,ifail=ifail)

      irevcm = 0
revcomm: Do
        Call f12jrf(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 B - A
          nnzz = 2*(nnza+nnzb)
!         We store it in the arrays az, irowz and icolz
          az(1:nnza) = -a(1:nnza)
          irowz(1:nnza) = irowa(1:nnza)
          icolz(1:nnza) = icola(1:nnza)
          az(nnza+1:2*nnza) = -a(1:nnza)
          irowz(nnza+1:2*nnza) = icola(1:nnza)
          icolz(nnza+1:2*nnza) = irowa(1:nnza)
!         Add the elements of ze B
          irowz(2*nnza+1:2*nnza+nnzb) = irowb(1:nnzb)
          icolz(2*nnza+1:2*nnza+nnzb) = icolb(1:nnzb)
          az(2*nnza+1:2*nnza+nnzb) = ze*b(1:nnzb)
          irowz(2*nnza+nnzb+1:2*nnza+2*nnzb) = icolb(1:nnzb)
          icolz(2*nnza+nnzb+1:2*nnza+2*nnzb) = irowb(1:nnzb)
          az(2*nnza+nnzb+1:2*nnza+2*nnzb) = conjg(ze)*b(1:nnzb)
!         Sort the elements into correct coordinate storage format
          Call f11znf(n,nnzz,az,irowz,icolz,'S','R',istr,iwork,ifail)
!         We have double counted the diagonal elements, so halve them
          Do i = 1, nnzz
            If (irowz(i)==icolz(i)) Then
              az(i) = 0.5_nag_wp*az(i)
            End If
          End Do

!         Form incomplete LU factorization of ze B - 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 B - A)^H
          nnzzh = 2*(nnza+nnzb)
!         We store it in the arrays azh, irowzh and icolzh
          azh(1:nnza) = -a(1:nnza)
          irowzh(1:nnza) = irowa(1:nnza)
          icolzh(1:nnza) = icola(1:nnza)
          azh(nnza+1:2*nnza) = -a(1:nnza)
          irowzh(nnza+1:2*nnza) = icola(1:nnza)
          icolzh(nnza+1:2*nnza) = irowa(1:nnza)
!         Add the elements of (ze B)^H
          irowzh(2*nnza+1:2*nnza+nnzb) = irowb(1:nnzb)
          icolzh(2*nnza+1:2*nnza+nnzb) = icolb(1:nnzb)
          azh(2*nnza+1:2*nnza+nnzb) = conjg(ze)*b(1:nnzb)
          irowzh(2*nnza+nnzb+1:2*nnza+2*nnzb) = icolb(1:nnzb)
          icolzh(2*nnza+nnzb+1:2*nnza+2*nnzb) = irowb(1:nnzb)
          azh(2*nnza+nnzb+1:2*nnza+2*nnzb) = ze*b(1:nnzb)
!         Sort the elements into correct coordinate storage format
          Call f11znf(n,nnzzh,azh,irowzh,icolzh,'S','R',istr,iwork,ifail)
!         We have double counted the diagonal elements, so halve them
          Do i = 1, nnzzh
            If (irowzh(i)==icolzh(i)) Then
              azh(i) = 0.5_nag_wp*azh(i)
            End If
          End Do

!         Form incomplete LU factorization of (ze B - 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 f11xsf(n,nnza,a,irowa,icola,'N',z(1,i),x(1,i),ifail)
          End Do
          If (ifail/=0) Then
            Exit revcomm
          End If
          Cycle revcomm
        Case (6)
!         Compute x <- Bz
          Do i = 1, m0
            Call f11xsf(n,nnzb,b,irowb,icolb,'N',z(1,i),x(1,i),ifail)
          End Do
          If (ifail/=0) Then
            Exit revcomm
          End If
          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)

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

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

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