    PROGRAM f12fgfe

!      F12FGF Example Program Text

!      Mark 23 Release. NAG Copyright 2011.

!      .. 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(:,:), mb(:,:), resid(:),     &
                                          v(:,:)
       INTEGER, ALLOCATABLE            :: icomm(:)
!      .. Intrinsic Functions ..
       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 super-diagonal.

       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 = 2*kl + ku + 1
       ldv = n
       ALLOCATE (mb(ldmb,n),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 20
       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)

20     CONTINUE
    END PROGRAM f12fgfe
