PROGRAM f12agfe

!      F12AGF Example Program Text

!      Mark 23 Release. NAG Copyright 2011.

!      .. 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 Functions ..
       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 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
             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 20
       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)

20     CONTINUE
    END PROGRAM f12agfe