Program f11bd_a1w_fe

!     F11BD_A1W_F Example Program Text
!     Mark 27 Release. NAG Copyright 2019.

!     .. Use Statements ..
      Use iso_c_binding, Only: c_ptr
      Use nagad_library, Only: f11bd_a1w_f, f11be_a1w_f,                       &
                               nagad_a1w_get_derivative,                       &
                               nagad_a1w_inc_derivative,                       &
                               nagad_a1w_ir_interpret_adjoint_sparse,          &
                               nagad_a1w_ir_register_variable,                 &
                               nagad_a1w_ir_remove, nagad_a1w_ir_zero_adjoints &
                               , nagad_a1w_w_rtype, x10aa_a1w_f, x10ab_a1w_f,  &
                               x10za_a1w_f, Operator (-), Assignment (=),      &
                               Operator (*), Operator (+), Operator (/)
      Use nag_library, Only: nag_wp, x04caf
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter               :: nin = 5, nout = 6
!     .. Local Scalars ..
      Type (nagad_a1w_w_rtype)         :: a, alpha, anorm, b1, bb, c, sigmax,  &
                                          tol
      Type (c_ptr)                     :: ad_handle
      Real (Kind=nag_wp)               :: alphar
      Integer                          :: i, ifail, irevcm, iterm, k, lwork,   &
                                          lwreq, m, maxitn, monit, n
!     .. Local Arrays ..
      Type (nagad_a1w_w_rtype), Allocatable :: b(:), work(:), x(:)
      Type (nagad_a1w_w_rtype)         :: wgt(1)
      Real (Kind=nag_wp), Allocatable  :: dx(:,:)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: real
!     .. Executable Statements ..
      Write (nout,*) 'F11BD_A1W_F Example Program Results'
!     Skip heading in data file
      Read (nin,*)
      Read (nin,*) n, m
      Read (nin,*) alphar

      lwork = 2*m*n + 1000

      Allocate (b(n),x(n),work(lwork),dx(n,2))

!     Create AD tape
      Call x10za_a1w_f

!     Create AD configuration data object
      ifail = 0
      Call x10aa_a1w_f(ad_handle,ifail)


      alpha = alphar
      b1 = 12.0_nag_wp
      a = 1.0_nag_wp
      bb = b1 - 2.0_nag_wp
      c = 1.0_nag_wp

!     Register variables to differentiate w.r.t.
      Call nagad_a1w_ir_register_variable(alpha)
      Call nagad_a1w_ir_register_variable(bb)

      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_a1w_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_a1w_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)%value,b(i)%value,i=1,n)

      Write (nout,*)
      Write (nout,*) ' Derivatives calculated: First order adjoints'
      Write (nout,*) ' Computational mode    : algorithmic'

      Write (nout,*)
      Write (nout,*) ' Derivatives of solution w.r.t. alpha and diag:'
!     Setup evaluation of derivatives via adjoints
      k = 0
      Do i = 1, n
        k = k + 1
        Call nagad_a1w_ir_zero_adjoints
        Call nagad_a1w_inc_derivative(x(i),1.0_nag_wp)
        ifail = 0
        Call nagad_a1w_ir_interpret_adjoint_sparse(ifail)
!       Get derivatives
        dx(k,1) = nagad_a1w_get_derivative(alpha)
        dx(k,2) = nagad_a1w_get_derivative(bb)
      End Do
      Call x04caf('General',' ',k,2,dx,n,'      d/dalpha    d/ddiag',ifail)

!     Remove computational data object and tape
100   Call x10ab_a1w_f(ad_handle,ifail)
      Call nagad_a1w_ir_remove

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_a1w_fe