Program f12fgfe

!     F12FGF Example Program Text

!     Mark 30.3 Release. nAG Copyright 2024.

!     .. Use Statements ..
      Use nag_library, Only: daxpy, dgbmv, dnrm2, f12fff, f12fgf, 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               :: inc1 = 1, iset = 1, nin = 5,         &
                                          nout = 6
!     .. Local Scalars ..
      Real (Kind=nag_wp)               :: h2, sigma
      Integer                          :: i, idiag, ifail, isub, isup, j, kl,  &
                                          ku, lcomm, ldab, ldmb, ldv, licomm,  &
                                          lo, n, nconv, ncv, nev, nx, outchn
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: ab(:,:), ax(:), comm(:), d(:),       &
                                          d_print(:,:), resid(:), v(:,:)
      Real (Kind=nag_wp)               :: mb(0)
      Integer, Allocatable             :: icomm(:)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: abs, int, max, real
!     .. Executable Statements ..
      Write (nout,*) 'F12FGF Example Program Results'
      Write (nout,*)

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

      Read (nin,*) nx, nev, ncv
      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 f12fff(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 f12fff(n,nev,ncv,icomm,licomm,comm,lcomm,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
      Allocate (ab(ldab,n))

!     Zero out AB.

      ab(1:ldab,1:n) = 0.0_nag_wp

!     Main diagonal of A.

      h2 = one/real((nx+1)*(nx+1),kind=nag_wp)
      idiag = kl + ku + 1
      ab(idiag,1:n) = 4.0_nag_wp/h2

!     First subdiagonal and superdiagonal of A.

      isup = kl + ku
      isub = kl + ku + 2

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

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

      End Do

!     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/h2
          ab(isub,j) = -one/h2
        End Do

      End Do

!     Find eigenvalues of largest magnitude and the corresponding
!     eigenvectors.

      ldmb = 0
      ldv = n
      Allocate (d(ncv),v(ldv,ncv+1),resid(n))

      ifail = -1
      Call f12fgf(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 (d_print(nconv,2),ax(n))
      d_print(1:nconv,1) = d(1:nconv)

      Do j = 1, nconv

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

!       The NAG name equivalent of daxpy is f06ecf
        Call daxpy(n,-d_print(j,1),v(1,j),inc1,ax,inc1)

!       The NAG name equivalent of dnrm2 is f06ejf
        d_print(j,2) = dnrm2(n,ax,1)
      End Do

      d_print(1:nconv,2) = d_print(1:nconv,2)/abs(d_print(1:nconv,1))

      Write (nout,*)
      Flush (nout)

      outchn = nout
      Call x04abf(iset,outchn)

      ifail = 0
      Call x04caf('G','N',nconv,2,d_print,nconv,' Ritz values and residuals',  &
        ifail)

100   Continue
    End Program f12fgfe