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

NAG FL Interface Introduction
Example description
    Program f12jvfe

!     F12JVF 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, f11xnf, f11znf, f12jaf, f12jgf,   &
                             f12jvf, 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, deg, i, ifail, ind, irevcm,     &
                                          iter, itn, j, k, la, ldx, ldy, ldz,  &
                                          liwork, lwork, m, m0, n, nconv,      &
                                          nnzc, nnzch, nnzmax, nnzsum, 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(:), nnza(:),         &
                                          tedge(:)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: conjg, max, maxval, min, sqrt, sum
!     .. Executable Statements ..
      Write (nout,*) 'F12JVF Example Program Results'
      Write (nout,*)
!     Skip heading in data file
      Read (nin,*)
      Read (nin,*) n, deg
      Allocate (nnza(deg+1))

      Read (nin,*) nnza(1:deg+1)
      nnzmax = maxval(nnza)
      nnzsum = sum(nnza)
      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*(nnzsum)
      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(nnzmax,deg+1),w(n),icol(nnzmax,deg+1),irow(nnzmax,deg+1),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 matrices from data file
      Do j = 1, deg + 1
        Do i = 1, nnza(j)
          Read (nin,*) a(i,j), irow(i,j), icol(i,j)
        End Do
      End Do

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

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

      zedge(1) = (-1.5_nag_wp,-0.5_nag_wp)
      zedge(2) = (0.5_nag_wp,1.5_nag_wp)

      tedge(1) = 100
      tedge(2) = 0

      nedge(1) = 20
      nedge(2) = 10

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

      irevcm = 0
revcomm: Do
        Call f12jvf(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 sparse matrix \sum ze^i A_i
          nnzz = nnzsum
          ind = 1
          Do i = 0, deg
            az(ind:ind+nnza(i+1)-1) = (ze**i)*a(1:nnza(i+1),i+1)
            irowz(ind:ind+nnza(i+1)-1) = irow(1:nnza(i+1),i+1)
            icolz(ind:ind+nnza(i+1)-1) = icol(1:nnza(i+1),i+1)
            ind = ind + nnza(i+1)
          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 \sum ze^i A_i
          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 (\sum ze^i A_i) 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 (\sum ze^i A_i)^H
          nnzzh = nnzsum
          ind = 1
          Do i = 0, deg
            azh(ind:ind+nnza(i+1)-1) = conjg((ze**i)*a(1:nnza(i+1),i+1))
            irowzh(ind:ind+nnza(i+1)-1) = icol(1:nnza(i+1),i+1)
            icolzh(ind:ind+nnza(i+1)-1) = irow(1:nnza(i+1),i+1)
            ind = ind + nnza(i+1)
          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 (\sum ze^i A_i)^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 (\sum ze^i A_i)^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 <- A_k z
          Do i = 1, m0
            Call f11xnf('N',n,nnza(k+1),a(1,k+1),irow(1,k+1),icol(1,k+1),'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_k^H z
          Do i = 1, m0
            Call f11xnf('T',n,nnza(k+1),a(1,k+1),irow(1,k+1),icol(1,k+1),'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)
        Call x04daf('General',' ',n,nconv,z(1,m0+1),ldz,'Left eigenvectors',   &
          ifail)

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

!     Destroy data handle
      Call f12jzf(handle,ifail)

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