Program f12aufe

!     F12AUF Example Program Text

!     Mark 26.2 Release. NAG Copyright 2017.

!     .. Use Statements ..
      Use nag_library, Only: dznrm2, f12arf, f12atf, f12auf, nag_wp, x04abf,   &
                             x04caf, zaxpy, zgbmv
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Complex (Kind=nag_wp), Parameter :: one = (1.0_nag_wp,0.0_nag_wp)
      Complex (Kind=nag_wp), Parameter :: zero = (0.0_nag_wp,0.0_nag_wp)
      Integer, Parameter               :: iset = 1, nin = 5, nout = 6
      Logical, Parameter               :: printr = .False.
!     .. Local Scalars ..
      Complex (Kind=nag_wp)            :: ch, sigma
      Real (Kind=nag_wp)               :: h, rho
      Integer                          :: i, idiag, ifail, isub, isup, j, kl,  &
                                          ku, lcomm, ldab, ldmb, ldv, licomm,  &
                                          lo, n, ncol, nconv, ncv, nev, nx,    &
                                          outchn
!     .. Local Arrays ..
      Complex (Kind=nag_wp), Allocatable :: ab(:,:), ax(:), comm(:), d(:),     &
                                          mb(:,:), mx(:), resid(:), v(:,:)
      Real (Kind=nag_wp), Allocatable  :: d_print(:,:)
      Integer, Allocatable             :: icomm(:)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: abs, aimag, cmplx, int, max, real
!     .. Executable Statements ..
      Write (nout,*) 'F12AUF Example Program Results'
      Write (nout,*)
      Flush (nout)

!     Skip heading in data file
      Read (nin,*)

      Read (nin,*) nx, nev, ncv
      Read (nin,*) sigma
      n = nx*nx

!     Initialize communication arrays.
!     Query the required sizes of the communication arrays.

      licomm = -1
      lcomm = -1
      Allocate (icomm(max(1,licomm)),comm(max(1,lcomm)))

      ifail = 0
      Call f12atf(n,nev,ncv,icomm,licomm,comm,lcomm,ifail)

      licomm = icomm(1)
      lcomm = int(comm(1))
      Deallocate (icomm,comm)
      Allocate (icomm(max(1,licomm)),comm(max(1,lcomm)))

      ifail = 0
      Call f12atf(n,nev,ncv,icomm,licomm,comm,lcomm,ifail)

!     Set the mode.

      ifail = 0
      Call f12arf('SHIFTED INVERSE',icomm,comm,ifail)

!     Set problem type

      ifail = 0
      Call f12arf('GENERALIZED',icomm,comm,ifail)

!     Construct the matrix A in banded form and store in AB.
!     KU, KL are number of superdiagonals and subdiagonals within
!     the band of matrices A and M.

      kl = nx
      ku = nx
      ldab = 2*kl + ku + 1
      ldmb = 2*kl + ku + 1
      Allocate (ab(ldab,n),mb(ldmb,n))

!     Zero out AB and MB.

      ab(1:ldab,1:n) = zero
      mb(1:ldmb,1:n) = zero

!     Main diagonal of A.

      idiag = kl + ku + 1
      ab(idiag,1:n) = cmplx(4.0_nag_wp,0.0_nag_wp,kind=nag_wp)
      mb(idiag,1:n) = ab(idiag,1)

!     First subdiagonal and superdiagonal of A.

      isup = kl + ku
      isub = kl + ku + 2
      rho = 100.0_nag_wp
      h = 1._nag_wp/real(nx+1,kind=nag_wp)
      ch = cmplx(0.5_nag_wp*h*rho,0.0_nag_wp,kind=nag_wp)

      Do i = 1, nx
        lo = (i-1)*nx

        Do j = lo + 1, lo + nx - 1
          ab(isub,j+1) = -one + ch
          ab(isup,j) = -one - ch
        End Do

      End Do

      mb(isub,2:n) = one
      mb(isup,1:n-1) = one

!     KL-th subdiagonal and KU-th superdiagonal.

      isup = kl + 1
      isub = 2*kl + ku + 1

      Do i = 1, nx - 1
        lo = (i-1)*nx

        Do j = lo + 1, lo + nx
          ab(isup,nx+j) = -one
          ab(isub,j) = -one
        End Do

      End Do

!     Find eigenvalues closest in value to SIGMA and corresponding
!     eigenvectors.

      ldv = n
      Allocate (d(nev),v(ldv,ncv),resid(n))

      ifail = -1
      Call f12auf(kl,ku,ab,ldab,mb,ldmb,sigma,nconv,d,v,ldv,resid,v,ldv,comm,  &
        icomm,ifail)

      If (ifail/=0) Then
        Go To 100
      End If

!     Compute the residual norm  ||A*x - lambda*x||.

      Allocate (ax(n),mx(n),d_print(nconv,3))
      d_print(1:nconv,1) = real(d(1:nconv))
      d_print(1:nconv,2) = aimag(d(1:nconv))

      Do j = 1, nconv

!       The NAG name equivalent of zgbmv is f06sbf
        Call zgbmv('N',n,n,kl,ku,one,ab(kl+1,1),ldab,v(1,j),1,zero,ax,1)

        Call zgbmv('N',n,n,kl,ku,one,mb(kl+1,1),ldmb,v(1,j),1,zero,mx,1)

        Call zaxpy(n,-d(j),mx,1,ax,1)
        d_print(j,3) = dznrm2(n,ax,1)/abs(d(j))

      End Do

      Write (nout,*)
      Flush (nout)

      outchn = nout
      Call x04abf(iset,outchn)

      If (printr) Then
!       Print residual associated with each Ritz value.
        ncol = 3
      Else
        ncol = 2
      End If
      ifail = 0
      Call x04caf('G','N',nconv,ncol,d_print,nconv,                            &
        ' Ritz values closest to sigma',ifail)

100   Continue
    End Program f12aufe