    PROGRAM c05pdfe

!      C05PDF Example Program Text

!      Mark 23 Release. NAG Copyright 2011.

!      .. Use Statements ..
       USE nag_library, ONLY : c05pdf, 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(:)
!      .. Intrinsic Functions ..
       INTRINSIC                          sqrt
!      .. Executable Statements ..
       WRITE (nout,*) 'C05PDF 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 c05pdf(irevcm,n,x,fvec,fjac,ldfjac,xtol,diag,mode,factor,r,lr, &
             qtf,w,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 c05pdfe
