Program e04ufae

!     E04UFA Example Program Text

!     Mark 28.6 Release. NAG Copyright 2022.

!     .. Use Statements ..
      Use nag_library, Only: dgemv, e04ufa, e04wbf, nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Real (Kind=nag_wp), Parameter    :: one = 1.0_nag_wp
      Real (Kind=nag_wp), Parameter    :: zero = 0.0_nag_wp
      Integer, Parameter               :: inc1 = 1, lcwsav = 5, liwsav = 610,  &
                                          llwsav = 120, lrwsav = 475, nin = 5, &
                                          nout = 6
!     .. Local Scalars ..
      Real (Kind=nag_wp)               :: objf
      Integer                          :: i, ifail, irevcm, iter, j, lda,      &
                                          ldcj, ldr, liwork, lwork, n, nclin,  &
                                          ncnln, sda, sdcjac
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable  :: a(:,:), bl(:), bu(:), c(:),          &
                                          cjac(:,:), clamda(:), objgrd(:),     &
                                          r(:,:), work(:), x(:)
      Real (Kind=nag_wp)               :: rwsav(lrwsav)
      Integer, Allocatable             :: istate(:), iwork(:), needc(:)
      Integer                          :: iwsav(liwsav)
      Logical                          :: lwsav(llwsav)
      Character (80)                   :: cwsav(lcwsav)
!     .. Intrinsic Procedures ..
      Intrinsic                        :: max
!     .. Executable Statements ..
      Write (nout,*) 'E04UFA Example Program Results'

!     Skip heading in data file.
      Read (nin,*)

      Read (nin,*) n, nclin, ncnln
      liwork = 3*n + nclin + 2*ncnln
      lda = max(1,nclin)

      If (nclin>0) Then
        sda = n
      Else
        sda = 1
      End If

      ldcj = max(1,ncnln)

      If (ncnln>0) Then
        sdcjac = n
      Else
        sdcjac = 1
      End If

      ldr = n

      If (ncnln==0 .And. nclin>0) Then
        lwork = 2*n**2 + 21*n + 11*nclin + 2
      Else If (ncnln>0 .And. nclin>=0) Then
        lwork = 2*n**2 + n*nclin + 2*n*ncnln + 21*n + 11*nclin + 22*ncnln + 1
      Else
        lwork = 21*n + 2
      End If

      Allocate (istate(n+nclin+ncnln),iwork(liwork),a(lda,sda),                &
        bl(n+nclin+ncnln),bu(n+nclin+ncnln),c(max(1,                           &
        ncnln)),cjac(ldcj,sdcjac),clamda(n+nclin+ncnln),objgrd(n),r(ldr,n),    &
        x(n),work(lwork),needc(max(1,ncnln)))

      If (nclin>0) Then
        Read (nin,*)(a(i,1:n),i=1,nclin)
      End If

      Read (nin,*) bl(1:(n+nclin+ncnln))
      Read (nin,*) bu(1:(n+nclin+ncnln))
      Read (nin,*) x(1:n)

!     Set all constraint Jacobian elements to zero.
!     Note that this will only work when 'Derivative Level = 3'
!     (the default; see Section 11.2).

      cjac(1:ncnln,1:n) = zero

!     Initialise E04UFA

      ifail = 0
      Call e04wbf('E04UFA',cwsav,lcwsav,lwsav,llwsav,iwsav,liwsav,rwsav,       &
        lrwsav,ifail)

!     Solve the problem.

      irevcm = 0
      ifail = -1

revcomm: Do

        Call e04ufa(irevcm,n,nclin,ncnln,lda,ldcj,ldr,a,bl,bu,iter,istate,c,   &
          cjac,clamda,objf,objgrd,r,x,needc,iwork,liwork,work,lwork,cwsav,     &
          lwsav,iwsav,rwsav,ifail)

!       On intermediate exit IFAIL should not have been changed
!       and IREVCM should be > 0.

        If (irevcm==0) Then
          Exit revcomm
        End If

        If (irevcm==1 .Or. irevcm==3) Then

!         Evaluate the objective function.

          objf = x(1)*x(4)*(x(1)+x(2)+x(3)) + x(3)
        End If

        If (irevcm==2 .Or. irevcm==3) Then

!         Evaluate the objective gradient.

          objgrd(1) = x(4)*(x(1)+x(1)+x(2)+x(3))
          objgrd(2) = x(1)*x(4)
          objgrd(3) = x(1)*x(4) + one
          objgrd(4) = x(1)*(x(1)+x(2)+x(3))
        End If

        If (irevcm==4 .Or. irevcm==6) Then

!         Evaluate the nonlinear constraint functions.

          If (needc(1)>0) Then
            c(1) = x(1)**2 + x(2)**2 + x(3)**2 + x(4)**2
          End If

          If (needc(2)>0) Then
            c(2) = x(1)*x(2)*x(3)*x(4)
          End If

        End If

        If (irevcm==5 .Or. irevcm==6) Then

!         Evaluate the constraint Jacobian.

          If (needc(1)>0) Then
            cjac(1,1) = x(1) + x(1)
            cjac(1,2) = x(2) + x(2)
            cjac(1,3) = x(3) + x(3)
            cjac(1,4) = x(4) + x(4)
          End If

          If (needc(2)>0) Then
            cjac(2,1) = x(2)*x(3)*x(4)
            cjac(2,2) = x(1)*x(3)*x(4)
            cjac(2,3) = x(1)*x(2)*x(4)
            cjac(2,4) = x(1)*x(2)*x(3)
          End If

        End If

      End Do revcomm

      Select Case (ifail)
      Case (0:6,8)
        Write (nout,*)
        Write (nout,99999)
        Write (nout,*)

        Do i = 1, n
          Write (nout,99998) i, istate(i), x(i), clamda(i)
        End Do

        If (nclin>0) Then

!         A*x --> work.
!         The NAG name equivalent of dgemv is f06paf
          Call dgemv('N',nclin,n,one,a,lda,x,inc1,zero,work,inc1)

          Write (nout,*)
          Write (nout,*)
          Write (nout,99997)
          Write (nout,*)

          Do i = n + 1, n + nclin
            j = i - n
            Write (nout,99996) j, istate(i), work(j), clamda(i)
          End Do

        End If

        If (ncnln>0) Then
          Write (nout,*)
          Write (nout,*)
          Write (nout,99995)
          Write (nout,*)

          Do i = n + nclin + 1, n + nclin + ncnln
            j = i - n - nclin
            Write (nout,99994) j, istate(i), c(j), clamda(i)
          End Do

        End If

        Write (nout,*)
        Write (nout,*)
        Write (nout,99993) objf
      End Select

99999 Format (1X,'Varbl',2X,'Istate',3X,'Value',9X,'Lagr Mult')
99998 Format (1X,'V',2(1X,I3),4X,1P,G14.6,2X,1P,G12.4)
99997 Format (1X,'L Con',2X,'Istate',3X,'Value',9X,'Lagr Mult')
99996 Format (1X,'L',2(1X,I3),4X,1P,G14.6,2X,1P,G12.4)
99995 Format (1X,'N Con',2X,'Istate',3X,'Value',9X,'Lagr Mult')
99994 Format (1X,'N',2(1X,I3),4X,1P,G14.6,2X,1P,G12.4)
99993 Format (1X,'Final objective value = ',1P,G15.7)
    End Program e04ufae