Program f12agfe

!     F12AGF Example Program Text

!     Mark 26.2 Release. NAG Copyright 2017.

!     .. Use Statements ..
      Use nag_library, Only: dgbmv, dnrm2, f06bnf, f12adf, f12aff, f12agf,     &
                             nag_wp, x04abf, x04caf
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Real (Kind=nag_wp), Parameter    :: one = 1.0_nag_wp
      Real (Kind=nag_wp), Parameter    :: zero = 0.0_nag_wp
      Integer, Parameter               :: iset = 1, nin = 5, nout = 6
      Logical, Parameter               :: printr = .False.
!     .. Local Scalars ..
      Real (Kind=nag_wp)               :: h, rho, sigmai, sigmar
      Integer                          :: i, idiag, ifail, isub, isup, j, kl,  &
                                          ku, lcomm, ldab, ldmb, ldv, licomm,  &
                                          lo, n, ncol, nconv, ncv, nev, nx,    &
                                          outchn
      Logical                          :: first
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: ab(:,:), ax(:), comm(:), di(:),      &
                                          dr(:), d_print(:,:), mb(:,:), mx(:), &
                                          resid(:), v(:,:)
      Integer, Allocatable             :: icomm(:)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: abs, int, max, real
!     .. Executable Statements ..
      Write (nout,*) 'F12AGF Example Program Results'
      Write (nout,*)
      Flush (nout)

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

      Read (nin,*) nx, nev, ncv, sigmar, sigmai
      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 f12aff(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 f12aff(n,nev,ncv,icomm,licomm,comm,lcomm,ifail)

!     Set the mode.

      ifail = 0
      Call f12adf('SHIFTED IMAGINARY',icomm,comm,ifail)

!     Set problem type

      ifail = 0
      Call f12adf('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) = 0.0_nag_wp
      mb(1:ldmb,1:n) = 0.0_nag_wp

!     Main diagonal of A.

      idiag = kl + ku + 1
      ab(idiag,1:n) = 4.0_nag_wp
      mb(idiag,1:n) = 4.0_nag_wp

!     First subdiagonal and superdiagonal of A.

      isup = kl + ku
      isub = kl + ku + 2
      rho = 100.0_nag_wp
      h = one/real(nx+1,kind=nag_wp)

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

        Do j = lo + 1, lo + nx - 1
          ab(isub,j+1) = -one + 0.5_nag_wp*h*rho
          ab(isup,j) = -one - 0.5_nag_wp*h*rho
        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 (dr(nev+1),di(nev+1),v(ldv,ncv),resid(n))

      ifail = -1
      Call f12agf(kl,ku,ab,ldab,mb,ldmb,sigmar,sigmai,nconv,dr,di,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||.

      first = .True.
      Allocate (ax(n),mx(n),d_print(nconv,3))
      d_print(1:nconv,1) = dr(1:nconv)
      d_print(1:nconv,2) = di(1:nconv)

      Do j = 1, nconv

        If (di(j)==zero) Then

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

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

          ax(1:n) = -dr(j)*mx(1:n) + ax(1:n)
          d_print(j,3) = dnrm2(n,ax,1)/abs(dr(j))
        Else If (first) Then

          Call dgbmv('N',n,n,kl,ku,one,ab(kl+1,1),ldab,v(1,j),1,zero,ax,1)

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

          ax(1:n) = -dr(j)*mx(1:n) + ax(1:n)

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

          ax(1:n) = di(j)*mx(1:n) + ax(1:n)
          d_print(j,3) = dnrm2(n,ax,1)

          Call dgbmv('N',n,n,kl,ku,one,ab(kl+1,1),ldab,v(1,j+1),1,zero,ax,1)

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

          ax(1:n) = -dr(j)*mx(1:n) + ax(1:n)

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

          ax(1:n) = -di(j)*mx(1:n) + ax(1:n)
!         The NAG name equivalent of dnrm2 is f06ejf
          d_print(j,3) = f06bnf(d_print(j,3),dnrm2(n,ax,1))
          d_print(j,3) = d_print(j,3)/f06bnf(dr(j),di(j))
          d_print(j+1,3) = d_print(j,3)
          first = .False.
        Else
          first = .True.
        End If

      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 f12agfe