!   F04QAF Example Program Text
!   Mark 23 Release. NAG Copyright 2011.

    MODULE f04qafe_mod
!      F04QAF Example Program Module:
!             Parameters and User-defined Routines

!      .. Use Statements ..
       USE nag_library, ONLY : nag_wp
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Parameters ..
       INTEGER, PARAMETER                  :: iset = 1, liuser = 1,            &
                                              lruser = 1, nin = 5, nout = 6
!      .. Local Scalars ..
       INTEGER                             :: ncols, nrows
    CONTAINS
       SUBROUTINE atimes(n,x,y)
!         Called by routine aprod. Returns Y = Y + A*X,
!         where A is not stored explicitly.

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Scalar Arguments ..
          INTEGER, INTENT (IN)                :: n
!         .. Array Arguments ..
          REAL (KIND=nag_wp), INTENT (IN)     :: x(n)
          REAL (KIND=nag_wp), INTENT (INOUT)  :: y(n)
!         .. Local Scalars ..
          INTEGER                             :: i, i1, i2, i3, il, j
!         .. Executable Statements ..
          DO j = 1, nrows - 2
             y(j) = y(j) + x(j) - x(j+nrows-1)
          END DO
          DO j = 1, ncols - 2
             i = j*nrows - 1
             y(i) = y(i) + x(i) - x(i+1)
             i1 = i + 1
             il = i1 + nrows - 3
             DO i = i1, il
                i2 = i - nrows
                IF (j==1) i2 = i2 + 1
                i3 = i + nrows
                IF (j==ncols-2) i3 = i3 - 1
                y(i) = y(i) - x(i2) - x(i-1) + 4.0_nag_wp*x(i) - x(i+1) - &
                   x(i3)
             END DO
             i = il + 1
             y(i) = y(i) - x(i-1) + x(i)
          END DO
          DO j = n - nrows + 3, n
             y(j) = y(j) - x(j-nrows+1) + x(j)
          END DO
          RETURN
       END SUBROUTINE atimes
       SUBROUTINE aprod(mode,m,n,x,y,ruser,lruser,iuser,liuser)

!         APROD returns
!         Y = Y + A*X          when MODE = 1
!         X = X + ( A**T )*Y   when MODE = 2
!         for a given X and Y.

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Scalar Arguments ..
          INTEGER, INTENT (IN)                :: liuser, lruser, m, n
          INTEGER, INTENT (INOUT)             :: mode
!         .. Array Arguments ..
          REAL (KIND=nag_wp), INTENT (INOUT)  :: ruser(lruser), x(n), y(m)
          INTEGER, INTENT (INOUT)             :: iuser(liuser)
!         .. Local Scalars ..
          INTEGER                             :: j, j1, j2
!         .. Executable Statements ..
          IF (mode/=2) THEN
             CALL atimes(n,x,y)
             DO j = 1, nrows - 2
                y(m) = y(m) + x(j)
             END DO
             DO j = 1, ncols - 2
                y(m) = y(m) + x(j*nrows-1) + x(j*nrows+nrows-2)
             END DO
             DO j = m - nrows + 2, n
                y(m) = y(m) + x(j)
             END DO
          ELSE
             CALL atimes(n,y,x)
             DO j = 1, nrows - 2
                x(j) = x(j) + y(m)
             END DO
             DO j = 1, ncols - 2
                j1 = j*nrows - 1
                j2 = j1 + nrows - 1
                x(j1) = x(j1) + y(m)
                x(j2) = x(j2) + y(m)
             END DO
             DO j = m - nrows + 2, n
                x(j) = x(j) + y(m)
             END DO
          END IF
          RETURN
       END SUBROUTINE aprod
    END MODULE f04qafe_mod
    PROGRAM f04qafe

!      F04QAF Example Main Program

!      .. Use Statements ..
       USE nag_library, ONLY : f04qaf, nag_wp, x04abf
       USE f04qafe_mod, ONLY : aprod, iset, liuser, lruser, ncols, nin, nout,  &
                               nrows
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Local Scalars ..
       REAL (KIND=nag_wp)                  :: acond, anorm, arnorm, atol,      &
                                              btol, c, conlim, damp, h, rnorm, &
                                              xnorm
       INTEGER                             :: i1, ifail, inform, itn, itnlim,  &
                                              k, m, msglvl, n, outchn
!      .. Local Arrays ..
       REAL (KIND=nag_wp), ALLOCATABLE     :: b(:), se(:), work(:,:), x(:)
       REAL (KIND=nag_wp)                  :: ruser(lruser)
       INTEGER                             :: iuser(liuser)
!      .. Executable Statements ..
       WRITE (nout,*) 'F04QAF Example Program Results'
       WRITE (nout,*)
       FLUSH (nout)
!      Skip heading in data file
       READ (nin,*)
       READ (nin,*) nrows, ncols
       n = ncols*nrows - 4
       m = n + 1
       ALLOCATE (b(m),se(n),work(n,2),x(n))
       outchn = nout
       CALL x04abf(iset,outchn)

       h = 0.1_nag_wp
!      Initialize rhs and other quantities required by F04QAF.
!      Convergence will be sooner if we do not regard A as exact,
!      so atol is not set to zero.
       b(1:n) = 0.0_nag_wp
       c = -h**2
       i1 = nrows
       DO k = 3, ncols
          b(i1:(i1+nrows-3)) = c
          i1 = i1 + nrows
       END DO
       b(m) = 1.0_nag_wp/h
       damp = 0.0_nag_wp
       atol = 1.0E-5_nag_wp
       btol = 1.0E-4_nag_wp
       conlim = 1.0_nag_wp/atol
       itnlim = 100
!      * Set msglvl to 2 to get output at each iteration *
       msglvl = 1

!      ifail: behaviour on error exit
!             =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
       ifail = 0
       CALL f04qaf(m,n,b,x,se,aprod,damp,atol,btol,conlim,itnlim,msglvl,itn, &
          anorm,acond,rnorm,arnorm,xnorm,work,ruser,lruser,iuser,liuser, &
          inform,ifail)

       WRITE (nout,*)
       WRITE (nout,*) 'Solution returned by F04QAF'
       WRITE (nout,99999) x(1:n)
       WRITE (nout,*)
       WRITE (nout,99998) 'Norm of the residual = ', rnorm

99999  FORMAT (1X,5F9.3)
99998  FORMAT (1X,A,1P,E12.2)
    END PROGRAM f04qafe