Program f07cd_t1w_fe

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

!     .. Use Statements ..
      Use iso_c_binding, Only: c_ptr
      Use nagad_library, Only: f07cd_t1w_f, f07ce_t1w_f,                       &
                               nagad_t1w_inc_derivative, nagad_t1w_w_rtype,    &
                               x10aa_t1w_f, x10ab_t1w_f, Assignment (=)
      Use nag_library, Only: nag_wp, x04caf
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter               :: nin = 5, nout = 6
!     .. Local Scalars ..
      Type (c_ptr)                     :: ad_handle
      Integer                          :: i, ifail, k, ldb, n, n1, nrhs
!     .. Local Arrays ..
      Type (nagad_t1w_w_rtype), Allocatable :: b(:,:), d(:), df(:), dl(:),     &
                                          dlf(:), du(:), du2(:), duf(:),       &
                                          x(:,:)
      Real (Kind=nag_wp), Allocatable  :: dxdb(:,:), dxdd(:,:), dxdl(:,:),     &
                                          dxdu(:,:), xv(:,:)
      Integer, Allocatable             :: ipiv(:)
!     .. Executable Statements ..
      Write (nout,*) 'F07CD_T1W_F Example Program Results'
      Write (nout,*)
      Flush (nout)
!     Skip heading in data file
      Read (nin,*)
      Read (nin,*) n, nrhs

      n1 = n - 1
      ldb = n
      Allocate (b(ldb,nrhs),d(n),dl(n1),du(n1),du2(n-2),ipiv(n))
      Allocate (df(n),dlf(n1),duf(n1))
      Allocate (x(ldb,nrhs))
      Allocate (dxdb(n,n),dxdd(n,n),dxdl(n,n1),dxdu(n,n1))
      Allocate (xv(ldb,nrhs))

!     Read the tridiagonal matrix A from data file

      Read (nin,*) dxdu(1:n1,1)
      Read (nin,*) dxdd(1:n,1)
      Read (nin,*) dxdl(1:n1,1)
      du(1:n-1) = dxdu(1:n1,1)
      d(1:n) = dxdd(1:n,1)
      dl(1:n1) = dxdl(1:n1,1)

!     Read the right hand matrix B
      Read (nin,*)(xv(i,1:nrhs),i=1,n)
      b(1:n,1:nrhs) = xv(1:n,1:nrhs)

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

      Do i = 1, 4*n - 2
        If (i<=n1) Then
          k = i
          Call nagad_t1w_inc_derivative(du(k),1.0_nag_wp)
        Else If (i<=n+n1) Then
          k = i - n1
          Call nagad_t1w_inc_derivative(d(k),1.0_nag_wp)
        Else If (i<=n+n1+n1) Then
          k = i - n1 - n
          Call nagad_t1w_inc_derivative(dl(k),1.0_nag_wp)
        Else
          k = i - n1 - n - n1
          Call nagad_t1w_inc_derivative(b(k,1),1.0_nag_wp)
        End If
        df = d
        duf = du
        dlf = dl

!       Factorize the tridiagonal matrix A
        ifail = 0
        Call f07cd_t1w_f(ad_handle,n,dlf,df,duf,du2,ipiv,ifail)

!       Solve the equations AX = B
        x = b
        ifail = 0
        Call f07ce_t1w_f(ad_handle,'N',n,nrhs,dlf,df,duf,du2,ipiv,x,ldb,ifail)

        If (i<=n1) Then
          du(k)%tangent = 0.0_nag_wp
          dxdu(1:n,k) = x(1:n,1)%tangent
        Else If (i<=n+n1) Then
          d(k)%tangent = 0.0_nag_wp
          dxdd(1:n,k) = x(1:n,1)%tangent
        Else If (i<=n+n1+n1) Then
          dl(k)%tangent = 0.0_nag_wp
          dxdl(1:n,k) = x(1:n,1)%tangent
        Else
          b(k,1)%tangent = 0.0_nag_wp
          dxdb(1:n,k) = x(1:n,1)%tangent
        End If
      End Do

!     Print the solution
      xv(1:n,1:nrhs) = x(1:n,1:nrhs)%value
      ifail = 0
      Call x04caf('General',' ',n,nrhs,xv,ldb,'Solution(s)',ifail)

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

      Write (nout,*)
      Write (nout,*) ' Derivatives of first solution column w.r.t. inputs:'

      ifail = 0
      Call x04caf('General',' ',n,n1,dxdu,n,'dx/d(du)',ifail)
      Call x04caf('General',' ',n,n,dxdd,n,'dx/d(d)',ifail)
      Call x04caf('General',' ',n,n1,dxdl,n,'dx/d(dl)',ifail)
      Call x04caf('General',' ',n,n,dxdb,n,'dx/d(b)',ifail)

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

    End Program f07cd_t1w_fe