    Program c05pdae

!     C05PDA Example Program Text

!     Mark 24 Release. NAG Copyright 2012.

!     .. Use Statements ..
      Use nag_library, Only: c05pda, dnrm2, nag_wp, x02ajf
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter               :: n = 9, nout = 6
      Integer, Parameter               :: ldfjac = n
!     .. Local Scalars ..
      Real (Kind=nag_wp)               :: factor, fnorm, xtol
      Integer                          :: icount, ifail, irevcm, j, k, lr, mode
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: diag(:), fjac(:,:), fvec(:), qtf(:), &
                                          r(:), w(:,:), x(:)
      Real (Kind=nag_wp)               :: rwsav(10)
      Integer                          :: iwsav(15)
      Logical                          :: lwsav(2)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: sqrt
!     .. Executable Statements ..
      Write (nout,*) 'C05PDA Example Program Results'

      Allocate (diag(n),fjac(ldfjac,n),fvec(n),qtf(n),r(n*(n+ &
        1)/2),w(n,4),x(n))

!     The following starting values provide a rough solution.

      x(1:n) = -1.0E0_nag_wp

      xtol = sqrt(x02ajf())
      diag(1:n) = 1.0E0_nag_wp
      mode = 2
      factor = 100.0E0_nag_wp
      icount = 0
      irevcm = 0
      ifail = -1

revcomm: Do

        Call c05pda(irevcm,n,x,fvec,fjac,ldfjac,xtol,diag,mode,factor,r,lr, &
          qtf,w,lwsav,iwsav,rwsav,ifail)

        Select Case (irevcm)
        Case (1)
          icount = icount + 1

!         Insert print statements here to monitor progess if desired

          Cycle revcomm
        Case (2)

!         Evaluate functions at current point

          fvec(1:n) = (3.0E0_nag_wp-2.0E0_nag_wp*x(1:n))*x(1:n) + 1.0E0_nag_wp
          fvec(2:n) = fvec(2:n) - x(1:(n-1))
          fvec(1:(n-1)) = fvec(1:(n-1)) - 2.0E0_nag_wp*x(2:n)
          Cycle revcomm
        Case (3)

!         Evaluate Jacobian at current point

          fjac(1:n,1:n) = 0.0E0_nag_wp

          Do k = 1, n
            fjac(k,k) = 3.0E0_nag_wp - 4.0E0_nag_wp*x(k)

            If (k/=1) Then
              fjac(k,k-1) = -1.0E0_nag_wp
            End If

            If (k/=n) Then
              fjac(k,k+1) = -2.0E0_nag_wp
            End If

          End Do

          Cycle revcomm
        Case Default
          Exit revcomm
        End Select

      End Do revcomm

      Select Case (ifail)
      Case (0)
!       The NAG name equivalent of dnrm2 is f06ejf
        fnorm = dnrm2(n,fvec,1)
        Write (nout,*)
        Write (nout,99999) 'Final 2 norm of the residuals after', icount, &
          ' iterations is ', fnorm
        Write (nout,*)
        Write (nout,*) 'Final approximate solution'
        Write (nout,99998)(x(j),j=1,n)
      Case (3:)
        Write (nout,*)
        Write (nout,*) 'Approximate solution'
        Write (nout,99998)(x(j),j=1,n)
      End Select

99999 Format (1X,A,I4,A,E12.4)
99998 Format (5X,3F12.4)
    End Program c05pdae
