Program f11bd_p0w_fe

!     F11BD_P0W_F Example Program Text
!     Mark 29.3 Release. NAG Copyright 2023.

!     .. Use Statements ..
      Use iso_c_binding, Only: c_ptr
      Use nagad_library, Only: f11bd_p0w_f, f11be_p0w_f
      Use nag_library, Only: nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter               :: nin = 5, nout = 6
!     .. Local Scalars ..
      Type (c_ptr)                     :: ad_handle
      Real (Kind=nag_wp)               :: a, alpha, anorm, b1, bb, c, sigmax,  &
                                          tol
      Integer                          :: i, ifail, irevcm, iterm, lwork,      &
                                          lwreq, m, maxitn, monit, n
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: b(:), work(:), x(:)
      Real (Kind=nag_wp)               :: wgt(1)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: real
!     .. Executable Statements ..
      Write (nout,*) 'F11BD_P0W_F Example Program Results'
!     Skip heading in data file
      Read (nin,*)
      Read (nin,*) n, m
      Read (nin,*) alpha

      lwork = 2*m*n + 1000
      Allocate (b(n),x(n),work(lwork))

      ifail = 0
      b1 = 12.0_nag_wp
      a = 1.0_nag_wp
      bb = b1 - 2.0_nag_wp
      c = 1.0_nag_wp

      Do i = 1, n
        b(i) = b1*real(i,kind=nag_wp)
      End Do
      b(n) = b(n) - real(n+1,kind=nag_wp)

      b(1) = b(1) + (b1-1.0_nag_wp)*alpha
      b(2:n-1) = b(2:n-1) + b1*alpha
      b(n) = b(n) + (b1-1.0_nag_wp)*alpha

      Do i = 1, n
        x(1:n) = 3.0_nag_wp
      End Do

!     Call F11BDF to initialize the solver
      iterm = 2
      maxitn = 800
      sigmax = 0.0_nag_wp
      tol = 1.0E-10_nag_wp
      monit = 0

      ifail = 0
      Call f11bd_p0w_f(ad_handle,'RGMRES','P','2','N',iterm,n,m,tol,maxitn,    &
        anorm,sigmax,monit,lwreq,work,lwork,ifail)

      irevcm = 0
      lwreq = lwork

      ifail = 1
loop: Do
        Call f11be_p0w_f(ad_handle,irevcm,x,b,wgt,work,lwreq,ifail)
        If (irevcm/=4) Then
          ifail = -1
          Select Case (irevcm)
          Case (-1)
!           b = A^Tx
            b(1) = bb*x(1) + a*x(2)
            Do i = 2, n - 1
              b(i) = c*x(i-1) + bb*x(i) + a*x(i+1)
            End Do
            b(n) = c*x(n-1) + bb*x(n)
          Case (1)
!           b = Ax
            b(1) = bb*x(1) + c*x(2)
            Do i = 2, n - 1
              b(i) = a*x(i-1) + bb*x(i) + c*x(i+1)
            End Do
            b(n) = a*x(n-1) + bb*x(n)
          Case (2)
            b(1:n) = x(1:n)/bb
          End Select
        Else If (ifail/=0) Then
          Write (nout,99997) ifail
          Go To 100
        Else
          Exit loop
        End If
      End Do loop

!     Output x
      Write (nout,99999)
      Write (nout,99998)(x(i),b(i),i=1,n)

100   Continue
99999 Format (/,2X,'  Solution vector',2X,' Residual vector')
99998 Format (1X,1P,E16.4,1X,E16.4)
99997 Format (1X,/,1X,' ** F11BE_a1W_F returned with IFAIL = ',I5)

    End Program f11bd_p0w_fe