!   F04QAF Example Program Text
!   Mark 28.4 Release. NAG Copyright 2022.

    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 = 0, lruser = 0,    &
                                          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) Then
              i2 = i2 + 1
            End If
            i3 = i + nrows
            If (j==ncols-2) Then
              i3 = i3 - 1
            End If
            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 f04qafe_mod, Only: aprod, iset, liuser, lruser, ncols, nin, nout,    &
                             nrows
      Use nag_library, Only: f04qaf, nag_wp, x04abf
!     .. 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