Program f11bd_t1w_fe

!     F11BD_T1W_F Example Program Text
!     Mark 30.1 Release. NAG Copyright 2024.

!     .. Use Statements ..
      Use iso_c_binding, Only: c_ptr
      Use nagad_library, Only: f11bd_t1w_f, f11be_t1w_f,                       &
                               nagad_t1w_set_derivative, nagad_t1w_w_rtype,    &
                               x10aa_t1w_f, x10ab_t1w_f, Assignment (=),       &
                               Operator (-), 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_t1w_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_t1w_w_rtype), Allocatable :: b(:), work(:), x(:)
      Type (nagad_t1w_w_rtype)         :: wgt(1)
      Real (Kind=nag_wp), Allocatable  :: dx(:,:)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: real
!     .. Executable Statements ..
      Write (nout,*) 'F11BD_T1W_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 configuration data object
      ifail = 0
      Call x10aa_t1w_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

!     Call F11BDF to initialize the solver

      iterm = 2
      maxitn = 800
      sigmax = 0.0_nag_wp
      tol = 1.0E-10_nag_wp
      monit = 0

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

        If (k==1) Then
          Call nagad_t1w_set_derivative(alpha,1.0_nag_wp)
        Else
          Call nagad_t1w_set_derivative(bb,1.0_nag_wp)
        End If

        Do i = 1, n
          x(1:n) = 3.0_nag_wp
        End Do
        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

        irevcm = 0
        lwreq = lwork

        ifail = 1
loop:   Do
          Call f11be_t1w_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
        If (k==1) Then
          alpha%tangent = 0.0_nag_wp
        Else
          bb%tangent = 0.0_nag_wp
        End If
        dx(1:n,k) = x(1:n)%tangent
      End Do

!     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 tangents'
      Write (nout,*) ' Computational mode    : algorithmic'

      Write (nout,*)
      Write (nout,*) ' Derivatives of solution w.r.t. alpha and bb:'
      Call x04caf('General',' ',n,2,dx,n,'      d/dalpha    d/dbb',ifail)

!     Remove computational data object
100   Call x10ab_t1w_f(ad_handle,ifail)

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

    End Program f11bd_t1w_fe