!   Mark 24 Release. NAG Copyright 2012.

    Module e05ucfe_mod

!     E05UCF Example Program Module:
!            Parameters and User-defined Routines

!     .. Use Statements ..
      Use nag_library, Only: nag_wp
!     .. Implicit None Statement ..
      Implicit None
    Contains
      Subroutine schwefel_obj(mode,n,x,objf,objgrd,nstate,iuser,ruser)

!       .. Scalar Arguments ..
        Real (Kind=nag_wp), Intent (Out)     :: objf
        Integer, Intent (Inout)              :: mode
        Integer, Intent (In)                 :: n, nstate
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Inout)   :: objgrd(n), ruser(*)
        Real (Kind=nag_wp), Intent (In)      :: x(n)
        Integer, Intent (Inout)              :: iuser(*)
!       .. Intrinsic Procedures ..
        Intrinsic                            :: abs, cos, sin, sqrt, sum
!       .. Executable Statements ..
        If (mode==0 .Or. mode==2) Then
!         Evaluate the objective function.
          objf = sum(x(1:n)*sin(sqrt(abs(x(1:n)))))
        End If

        If (mode==1 .Or. mode==2) Then
!         Calculate the gradient of the objective function.
          objgrd(1:n) = sin(sqrt(abs(x(1:n)))) + 0.5_nag_wp*sqrt(abs(x(1:n)))* &
            cos(sqrt(abs(x(1:n))))
        End If

        Return
      End Subroutine schwefel_obj
      Subroutine schwefel_confun(mode,ncnln,n,ldcjsl,needc,x,c,cjsl,nstate, &
        iuser,ruser)

!       .. Scalar Arguments ..
        Integer, Intent (In)                 :: ldcjsl, n, ncnln, nstate
        Integer, Intent (Inout)              :: mode
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (Out)     :: c(ncnln)
        Real (Kind=nag_wp), Intent (Inout)   :: cjsl(ldcjsl,n), ruser(*)
        Real (Kind=nag_wp), Intent (In)      :: x(n)
        Integer, Intent (Inout)              :: iuser(*)
        Integer, Intent (In)                 :: needc(ncnln)
!       .. Local Scalars ..
        Real (Kind=nag_wp)                   :: t1, t2
        Integer                              :: k
!       .. Intrinsic Procedures ..
        Intrinsic                            :: cos, sin
!       .. Executable Statements ..
        If (mode==0 .Or. mode==2) Then
!         Constraint values are required.
!         Only those for which needc is non-zero need be set.
          Do k = 1, ncnln
            If (needc(k)>0) Then
              Select Case (k)
              Case (1)
                c(k) = x(1)**2 - x(2)**2 + 3.0_nag_wp*x(1)*x(2)
              Case (2)
                c(k) = cos((x(1)/200.0_nag_wp)**2+(x(2)/100.0_nag_wp))
              End Select
            End If
          End Do
        End If

        If (mode==1 .Or. mode==2) Then
!         Constraint derivatives are required.
          Do k = 1, ncnln
            Select Case (k)
            Case (1)
              cjsl(k,1) = 2.0_nag_wp*x(1) + 3.0_nag_wp*x(2)
              cjsl(k,2) = -2.0_nag_wp*x(2) + 3.0_nag_wp*x(1)
            Case (2)
              t1 = x(1)/200.0_nag_wp
              t2 = x(2)/100.0_nag_wp
              cjsl(k,1) = -sin(t1**2+t2)*2.0_nag_wp*t1/200.0_nag_wp
              cjsl(k,2) = -sin(t1**2+t2)/100.0_nag_wp
            End Select
          End Do
        End If

        Return
      End Subroutine schwefel_confun
      Subroutine mystart(npts,quas,n,repeat,bl,bu,iuser,ruser,mode)
!       Sets an initial point.
!       A typical user-defined start procedure (not used here)

!       .. Scalar Arguments ..
        Integer, Intent (Inout)              :: mode
        Integer, Intent (In)                 :: n, npts
        Logical, Intent (In)                 :: repeat
!       .. Array Arguments ..
        Real (Kind=nag_wp), Intent (In)      :: bl(n), bu(n)
        Real (Kind=nag_wp), Intent (Inout)   :: quas(n,npts), ruser(*)
        Integer, Intent (Inout)              :: iuser(*)
!       .. Executable Statements ..
!       quas(1:n,1:npts) is pre-assigned to zero.
        quas(1,1) = 420.9687_nag_wp
        quas(2,2) = 1.0_nag_wp !etc
        Return
      End Subroutine mystart
    End Module e05ucfe_mod
    Program e05ucfe

!     E05UCF Example Main Program

!     .. Use Statements ..
      Use nag_library, Only: dgemv, e05ucf, e05ucz, e05zkf, nag_wp
      Use e05ucfe_mod, Only: schwefel_confun, schwefel_obj
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter                   :: n = 2, nclin = 1, ncnln = 2,     &
                                              nin = 5, nout = 6
!     .. Local Scalars ..
      Integer                              :: i, ifail, j, k, l, lda, ldc,     &
                                              ldcjac, ldclda, ldobjd, ldr,     &
                                              ldx, liopts, listat, lopts, nb,  &
                                              npts, sda, sdcjac, sdr
      Logical                              :: repeat
!     .. Local Arrays ..
      Real (Kind=nag_wp), Allocatable      :: a(:,:), bl(:), bu(:), c(:,:),    &
                                              cjac(:,:,:), clamda(:,:),        &
                                              objf(:), objgrd(:,:), opts(:),   &
                                              r(:,:,:), work(:), x(:,:)
      Real (Kind=nag_wp)                   :: ruser(1)
      Integer, Allocatable                 :: info(:), iopts(:), istate(:,:),  &
                                              iter(:)
      Integer                              :: iuser(1)
!     .. Executable Statements ..
      Write (nout,*) 'E05UCF Example Program Results'
      Flush (nout)

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

      Read (nin,*) nb, npts
      Read (nin,*) repeat

      lda = nclin

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

      ldx = n
      ldobjd = n
      ldc = ncnln
      ldcjac = ncnln

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

      ldr = n
      sdr = n
      ldclda = n + nclin + ncnln
      listat = n + nclin + ncnln
      liopts = 740
      lopts = 485
      Allocate (a(lda,sda),bl(n+nclin+ncnln),bu(n+nclin+ncnln),x(ldx,nb), &
        objf(nb),objgrd(ldobjd,nb),iter(nb),c(ldc,nb),cjac(ldcjac,sdcjac,nb), &
        r(ldr,sdr,nb),clamda(ldclda,nb),istate(listat,nb),iopts(liopts), &
        opts(lopts),info(nb),work(nclin))

      bl(1:n+nclin+ncnln) = (/-500.0_nag_wp,-500.0_nag_wp,-10000.0_nag_wp, &
        -1.0_nag_wp,-0.9_nag_wp/)
      bu(1:n+nclin+ncnln) = (/500.0_nag_wp,500.0_nag_wp,10.0_nag_wp, &
        500000.0_nag_wp,0.9_nag_wp/)

      a(1,1) = 3.0_nag_wp
      a(1,2) = -2.0_nag_wp

!     Initialize the solver.

      ifail = 0
      Call e05zkf('Initialize = E05UCF',iopts,liopts,opts,lopts,ifail)

!     Solve the problem.
!     USE mystart from e05ucfe_mod and pass this instead of e05ucz to e05ucf
!     to set your own start points.

      ifail = -1
      Call e05ucf(n,nclin,ncnln,a,lda,bl,bu,schwefel_confun,schwefel_obj,npts, &
        x,ldx,e05ucz,repeat,nb,objf,objgrd,ldobjd,iter,c,ldc,cjac,ldcjac, &
        sdcjac,r,ldr,sdr,clamda,ldclda,istate,listat,iopts,opts,iuser,ruser, &
        info,ifail)

      Select Case (ifail)
      Case (0)
        l = nb
      Case (8)
        l = info(nb)
        Write (nout,'(1X,I16,A)') iter(nb), ' starting points converged'
      Case Default
        Go To 100
      End Select

loop: Do i = 1, l
        Write (nout,99999) i
        Write (nout,99998) info(i)
        Write (nout,99997) 'Varbl'
        Do j = 1, n
          Write (nout,99996) 'V', j, istate(j,i), x(j,i), clamda(j,i)
        End Do
        If (nclin>0) Then
          Write (nout,99997) 'L Con'

!         Below is a call to the level 2 BLAS routine DGEMV.
!         This performs the matrix vector multiplication A*X
!         (linear constraint values) and puts the result in
!         the first NCLIN locations of WORK.

          Call dgemv('N',nclin,n,1.0_nag_wp,a,lda,x(1,i),1,0.0_nag_wp,work,1)

          Do k = n + 1, n + nclin
            j = k - n
            Write (nout,99996) 'L', j, istate(k,i), work(j), clamda(k,i)
          End Do
        End If
        If (ncnln>0) Then
          Write (nout,99997) 'N Con'
          Do k = n + nclin + 1, n + nclin + ncnln
            j = k - n - nclin
            Write (nout,99996) 'N', j, istate(k,i), c(j,i), clamda(k,i)
          End Do
        End If
        Write (nout,99995) objf(i)
        Write (nout,99994)
        Write (nout,99993)(clamda(k,i),k=1,n+nclin+ncnln)

        If (l==1) Then
          Exit loop
        End If

        Write (nout,*)

        Write (nout,*) &
          ' ------------------------------------------------------ '

      End Do loop

100   Continue

99999 Format (/1X,'Solution number',I16)
99998 Format (/1X,'Local minimization exited with code',I5)
99997 Format (/1X,A,2X,'Istate',3X,'Value',9X,'Lagr Mult'/)
99996 Format (1X,A,2(1X,I3),4X,F12.4,2X,F12.4)
99995 Format (/1X,'Final objective value = ',1X,F12.4)
99994 Format (/1X,'QP multipliers')
99993 Format (1X,F12.4)
    End Program e05ucfe
