!   F04QAF Example Program Text
!   Mark 25 Release. NAG Copyright 2014.

    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
!     .. Accessibility Statements ..
      Private
      Public                               :: aprod
!     .. Parameters ..
      Integer, Parameter, Public           :: iset = 1, liuser = 1,            &
                                              lruser = 1, nin = 5, nout = 6
!     .. Local Scalars ..
      Integer, Public, Save                :: ncols, nrows
    Contains
      Subroutine atimes(n,x,y)
!       Called by routine aprod. Returns Y = Y + A*X,
!       where A is not stored explicitly.

!       .. 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.

!       .. 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