!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
! File:  ct25opt.f90
! Subroutines related to SNOPT (solver, wrapper, etc).
!
! 24 Dec 2008: Version for v4 of snctrl.
! 09 Feb 2010: v5.
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

module ct25opt
  use precision,    only : ip, rp
  use ctWorkModule, only : ctWork, ctGrid
  use ctGlobal,     only : nY, nU, nP, nC, nPhs, iDisc
  use ct20eval,     only : sncEvalTR, sncEvalHS
  implicit none

  private
  public  :: c2solv
  private :: snoptBKer, s0fgctrl, ctrlObj, ctrlCon

  ! Grid pointers (only for use in ct25opt)
  real(rp),    private, pointer :: gstep(:)
  integer(ip), private, pointer :: gndPtr(:)

contains

  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  subroutine c2solv ( Start, prob, x, hs, work, grid, odecon, algcon, &
                      INFO, mincw, miniw, minrw, nS, nInf, sInf, &
                      cu, lencu, iu, leniu, ru, lenru, &
                      cw, lencw, iw, leniw, rw, lenrw )
    integer(ip), intent(in)    :: Start, lencu, leniu, lenru, lencw, leniw, lenrw
    character(8),intent(in)    :: prob
    integer(ip), intent(inout) :: nS, hs(:), iu(leniu), iw(leniw)
    real(rp),    intent(inout) :: x(:), ru(lenru), rw(lenrw)
    character(8),intent(inout) :: cw(lencw), cu(lencu)
    integer(ip), intent(out)   :: INFO, mincw, miniw, minrw, nInf
    real(rp),    intent(out)   :: sInf
    external     :: odecon, algcon
    type(ctWork) :: work
    type(ctGrid) :: grid

    !===========================================================================
    ! Calls the SNOPT(B) kernel.
    !
    ! 24 Dec 2008: Renamed sncSolve.
    ! 09 Feb 2010: v5.
    !===========================================================================
    integer(ip)  :: m, n, neJ, nnObj, nnJac, nnCon, iObj, nName
    parameter     ( nName = 1 )
    character(4) :: CStart
    character(8) :: Names(nName)
    real(rp)     :: ObjAdd, Obj
    external     :: snLog, snLog2, sqLog, snSTOP
    integer(ip), pointer :: indJ(:), locJ(:)
    real(rp),    pointer :: Jcol(:), bl(:), bu(:), rc(:), pi(:)

    ! Set global grid stuff (passed to ctrlCon below).
    gstep   => grid%step
    gndPtr  => grid%ndPtr

    ObjAdd = 0.0

    CStart = 'Cold'
    if ( Start == 2 ) then
       CStart = 'Warm'
    end if

    m     = work%m
    n     = work%n
    neJ   = work%lenJ
    nnObj = 0
    nnJac = work%n
    nnCon = work%nnCon
    iObj  = work%m

    Jcol => work%Jcol
    indJ => work%indJ
    locJ => work%locJ
    bl   => work%bl
    bu   => work%bu
    pi   => work%pi
    rc   => work%rc

    call snoptBKer ( CStart, m, n, neJ, nName, nnCon, nnObj, nnJac, &
                     iObj, ObjAdd, prob, odecon, algcon, &
                     snLog, snLog2, sqLog, snSTOP, &
                     Jcol, indJ, locJ, bl, bu, Names, hs, x, pi, rc, &
                     INFO, mincw, miniw, minrw, nS, nInf, sInf, Obj, &
                     cu, lencu, iu, leniu, ru, lenru, &
                     cw, lencw, iw, leniw, rw, lenrw )

    nullify ( gstep, gndPtr )

  end subroutine c2solv

  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  subroutine ctrlCon ( mode, nnCon, nnJac, neJac, x, fCon, gCon, Status, &
                       odecon, algcon, cu, lencu, iu, leniu, ru, lenru )
    integer(ip), intent(in)    :: mode, nnCon, nnJac, neJac, Status, &
                                  lencu, leniu, lenru
    real(rp),    intent(in)    :: x(nnJac)
    integer(ip), intent(inout) :: iu(leniu)
    real(rp),    intent(inout) :: ru(lenru)
    character(8),intent(inout) :: cu(lencu)
    real(rp),    intent(out)   :: fCon(nnCon), gCon(neJac)
    external     :: odecon, algcon

    !===========================================================================
    ! Returns the constraint/derivative values.
    !===========================================================================
    integer(ip) :: s

    s  = size(gstep)

    ! SNCTRLA
    if ( iDisc == 0 ) then
       call sncEvalTR ( Status, mode, nY, nU, nP, nC, nPhs, gndPtr, &
                        gstep, s, x, nnJac, fCon, nnCon, gCon, neJac, &
                        odecon, algcon, cu, lencu, iu, leniu, ru, lenru )
    else
       call sncEvalHS ( Status, mode, nY, nU, nP, nC, nPhs, gndPtr, &
                        gstep, s, x, nnJac, fCon, nnCon, gCon, neJac, &
                        odecon, algcon, cu, lencu, iu, leniu, ru, lenru )
    end if

  end subroutine ctrlCon

  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  subroutine ctrlObj ( mode, nnObj, x, fObj, gObj, Status, odecon, &
                       cu, lencu, iu, leniu, ru, lenru )
    integer(ip), intent(in)    :: mode, nnObj, Status, lencu, leniu, lenru
    real(rp),    intent(in)    :: x(nnObj)
    integer(ip), intent(inout) :: iu(leniu)
    real(rp),    intent(inout) :: ru(lenru)
    character(8),intent(inout) :: cu(lencu)
    real(rp),    intent(out)   :: fObj, gObj(nnObj)
    external     :: odecon

    !===========================================================================
    ! This is funObj for SNOPTB.
    ! Nothing to do!  Linear objective!
    !===========================================================================

  end subroutine ctrlObj

  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  subroutine snoptBKer ( Start, m, n, ne, nName, nnCon, nnObj, nnJac, &
                         iObj, ObjAdd, Prob, odecon, algcon, &
                         snLog, snLog2, sqLog, snSTOP, &
                         Jcol, indJ, locJ, bl, bu, Names, hs, x, pi, rc, &
                         INFO, mincw, miniw, minrw, nS, nInf, sInf, Obj, &
                         cu, lencu, iu, leniu, ru, lenru, &
                         cw, lencw, iw, leniw, rw, lenrw )
    external    :: odecon, algcon, snLog, snLog2, sqLog, snSTOP
    integer(ip) :: INFO, iObj, lencu, lencw, leniu, leniw, lenru, lenrw, m, &
                   mincw, miniw, minrw, n, ne, nInf, nName, nnCon, &
                   nnJac, nnObj, nS, hs(n+m), indJ(ne), iu(leniu), iw(leniw), &
                   locJ(n+1)
    real(rp)      :: Obj, ObjAdd, sInf, Jcol(ne), bl(n+m), bu(n+m), pi(m), &
                     rc(n+m), ru(lenru), rw(lenrw), x(n+m)
    character*(*) :: Start
    character(8)  :: Prob, cu(lencu), cw(lencw), Names(nName)
    ! ==========================================================================
    ! snKerB does the work for snOptB. (Kernel for snoptB)
    !
    ! Developers can call this version with customized versions of
    ! snLog, snLog2  and  snSTOP.
    !
    ! 06 Sep 2008: SNOPTB adapted for control interface.
    ! 24 Dec 2008: Renamed snoptBKer.
    ! ==========================================================================
    character    :: Solver*6, str*80, str2*80
    integer(ip)  :: Useriw(130)
    real(rp)     :: Userrw(130)
    character(8) :: Usercw(130)
    logical      :: FPonly, gotR, PrtMem
    integer(ip) :: &
         Errors, Hcalls, Htype, inform, lenR, lenx0, lgObj, lhEtyp, &
         lkx, llocG, liwEst, lrwEst, lvlHes, lvlSrt, lx0, maxcw, &
         maxiw, maxR, maxrw, maxS, minmax, mProb, mQNmod, nb, negCon, &
         ngObj, ngObj0, ngQP, nextcw, nextiw, nextrw, nkx, nlocG, &
         nlocJ, nMajor, nnH, nnH0, nrhs, nrhs0, nx0, stkyOp
    real(rp) :: &
         fObj, ObjTru, rhs(1), x0(1)
    external :: &
         sqHx, s8qpHx
    ! ------------------------------------------------------------------
    integer(ip) ::            COLD,       BASIS,      WARM,       HOT
    parameter         (COLD   = 0, BASIS  = 1, WARM   = 2, HOT    = 3)
    integer(ip) ::            StdIn
    parameter         (StdIn  = 2)
    integer(ip) ::            HUnset,     HNorml
    parameter         (HUnset =-1, HNorml = 0)

    parameter         (mProb  =  51) ! Problem name
    parameter         (lvlSrt =  69) ! cold:warm:basis:hot start
    parameter         (Hcalls = 188) ! number of Hessian-vector prods
    parameter         (Htype  = 202) ! Current Hessian type
    ! ------------------------------------------------------------------
    Solver = 'SNCTRL'
    INFO   = 0

    ! ------------------------------------------------------------------
    ! Check memory limits and fetch the workspace starting positions.
    ! ------------------------------------------------------------------
    call s2Mem0 &
         ( INFO, Solver, lencw, leniw, lenrw, iw, &
         mincw, miniw, minrw, maxcw, maxiw, maxrw, &
         nextcw, nextiw, nextrw )
    if (INFO > 0) go to 999 ! Exit without printing

    ! Save the user's option choices  (weird choices get overwritten).
    ! Initialize timers and the standard input file.

    call chcopy( 130, cw(51), 1, Usercw, 1 )
    call icopy ( 130, iw(51), 1, Useriw, 1 )
    call dcopy ( 130, rw(51), 1, Userrw, 1 )

    call s1time( 0, 0, iw, leniw, rw, lenrw  )
    call s1file( StdIn, iw, leniw )

    ! Check the input arguments.

    call s3argB &
         ( inform, Start, m, n, ne, nName, nS, &
         nnCon, nnObj, nnJac, iObj, &
         indJ, locJ, bl, bu, Names, hs, pi, iw(lvlSrt), Errors, &
         iw, leniw, rw, lenrw )
    if (inform > 0) then
       INFO = inform
       go to 800
    end if

    ! ------------------------------------------------------------------
    ! The obligatory call to snInit has already set the defaults.
    ! Check that the optional parameters have sensible values.
    ! Print the options.
    ! ------------------------------------------------------------------
    cw(mProb) = Prob

    call s8dflt &
         ( m, n, nnCon, nnJac, nnObj, &
         cw, lencw, iw, leniw, rw, lenrw )
    call s3prtB &
         ( m, n, nnCon, nnJac, nnObj, iw, leniw, rw, lenrw )

    ! ------------------------------------------------------------------
    ! Compute the storage requirements for SNOPT  from the following
    ! variables:
    ! m,      n,     ne
    ! lenR  , maxS
    ! nnCon , nnJac, nnObj,
    ! negCon
    ! All have to be known before calling s2Mem.
    ! The only one in doubt is negCon, the number of Jacobian elements.
    ! Count them in s8Gsiz.
    ! ------------------------------------------------------------------
    nb      = n + m
    nlocJ   = n + 1
    nkx     = nb

    call s8Gsiz &
         ( m, nnCon, nnJac, ne, nlocJ, locJ, indJ, negCon )

    maxR    = iw( 52) ! max columns of R.
    maxS    = iw( 53) ! max # of superbasics
    mQNmod  = iw( 54) ! (ge 0) max # of BFGS updates
    lvlHes  = iw( 72) ! 0, 1,  2  => LM, FM, Exact Hessian
    minmax  = iw( 87) ! 1, 0, -1  => MIN, FP, MAX

    lenR    = maxR*(maxR + 1)/2  +  (maxS - maxR)

    nnH     = max( nnJac, nnObj )
    ngObj   = nnObj   ! Local nnObj is altered for FP

    ! Load iw with various problem dimensions.

    iw( 15) = n       ! copy of the number of columns
    iw( 16) = m       ! copy of the number of rows
    iw( 17) = ne      ! copy of the number of nonzeros in Jcol
    iw( 20) = negCon  ! # of nonzeros in gCon
    iw( 21) = nnJac   ! # of Jacobian  variables
    iw( 22) = nnObj   ! # of objective variables
    iw( 23) = nnCon   ! # of nonlinear constraints
    iw( 24) = nnH     !   max( nnObj, nnJac )
    iw( 28) = lenR    ! R(lenR) is the reduced Hessian factor
    iw(204) = iObj    ! position of the objective row in J

    ! ------------------------------------------------------------------
    ! If only a feasible point is requested, save the base point for the
    ! objective function:  1/2 || x - x0 ||^2
    ! ------------------------------------------------------------------
    FPonly  = minmax .eq. 0
    if ( FPonly ) then
       ngObj  = nnH

       lx0    = nextiw
       nextrw = lx0    + nnH
       minrw  = nextrw - 1
       if (minrw <= lenrw) then
          iw(298) = lx0
          call dcopy ( nnH, x, 1, rw(lx0), 1 )
       end if
    end if

    ! ------------------------------------------------------------------
    ! Allocate the local arrays for snOpt.
    ! s8Map  maps snOpt integer and double arrays.
    ! s2BMap maps the arrays for the LU routines.
    ! s2Mem  checks what space is available and prints any messages.
    ! ------------------------------------------------------------------
    call s8Map &
         ( m, n, negCon, nkx, nnCon, nnJac, ngObj, &
         lenR, maxR, maxS, mQNmod, lvlHes, &
         nextcw, nextiw, nextrw, iw, leniw )
    call s2Bmap &
         ( m, n, ne, maxS, &
         nextiw, nextrw, maxiw, maxrw, liwEst, lrwEst, iw, leniw )
    PrtMem = .true.           ! Print all messages in s2Mem
    call s2Mem &
         ( inform, PrtMem, liwEst, lrwEst, &
         nextcw, nextiw, nextrw, &
         maxcw, maxiw, maxrw, lencw, leniw, lenrw, &
         mincw, miniw, minrw, iw )
    if (inform /= 0) then
       INFO = inform
       go to 800
    end if

    ! Define the row and column ordering for J.
    ! SNOPT  uses natural order throughout, so kx = kxN.

    iw(247) = nkx     ! dimension of kx and its inverse, kxN
    lkx     = iw(251) ! j  = kx (jN) => col j of Jcol is variable jN
    iw(252) = lkx     ! jN = kxN(j ) => col j of Jcol is variable jN

    call s1perm( n, iw(lkx) )
    call s1perm( m, iw(lkx+n) )

    ! ------------------------------------------------------------------
    ! Construct column pointers for the nonlinear part of the  Jacobian.
    ! ------------------------------------------------------------------
    if (nnCon > 0) then
       llocG = iw(260) ! locG(nlocG) = column pointers for indG
       nlocG = nnJac + 1

       call s8Gloc &
            ( nnCon, nnJac, &
            ne, nlocJ, locJ, indJ, negCon, nlocG, iw(llocG) )
    end if

    ! ------------------------------------------------------------------
    ! Solve the problem.
    ! ------------------------------------------------------------------
    if (nnH == 0) then

       ! The problem is a linear program.

       nrhs   = 0             ! No constraint rhs vector.
       nx0    = 0             ! No constant shift for x.
       nrhs0  = max( nrhs , 1   )
       lenx0  = max( nx0  , 1   )
       ngObj0 = max( ngObj, 1   )
       nnH0   = max( nnH  , 1   )
       ngQP   = max( ngObj, nnH )

       lhEtyp = iw(283)      ! hEtype(nb) definition of elastic vars
       lgObj  = iw(297)      ! gObj(ngObj) = Objective gradient

       call iload ( nb, 3, iw(lhEtyp), 1 )

       call s5solv &
            ( INFO, Solver, iw(lvlSrt), &
            sqHx, s8qpHx, sqLog, iw(Hcalls), gotR, &
            m, n, nb, nnH0, nnH, nName, ngQP, ngObj0, ngObj, &
            iObj, ObjAdd, fObj, ObjTru, nInf, sInf, &
            ne, nlocJ, locJ, indJ, Jcol, &
            bl, bu, rw(lgObj), Names, &
            nrhs0, nrhs, rhs, lenx0, nx0, x0, &
            iw(lhEtyp), hs, x, pi, rc, nS, &
            cu, lencu, iu, leniu, ru, lenru, &
            cw, lencw, iw, leniw, rw, lenrw )
    else

       ! The problem is nonlinear.
       ! Define the type of initial Hessian.

       if      (iw(lvlSrt) == COLD ) then
          iw(Htype)  = HUnset
       else if (iw(lvlSrt) == BASIS) then
          iw(Htype)  = HUnset
       else if (iw(lvlSrt) == WARM ) then
          iw(Htype)  = HUnset
       else if (iw(lvlSrt) == HOT  ) then
          iw(Htype)  = HNorml
       end if

       call s8solv &
            ( INFO, Solver, iw(lvlSrt), &
            s0fgctrl, odecon, algcon, &
            snLog, snLog2, snSTOP, gotR, &
            m, n, nb, nnCon, nnJac, ngObj, &
            nName, iObj, ObjAdd, fObj, ObjTru, nInf, sInf, &
            ne, nlocJ, locJ, indJ, Jcol, &
            bl, bu, Names, &
            hs, x, pi, rc, nMajor, nS, &
            cu, lencu, iu, leniu, ru, lenru, &
            cw, lencw, iw, leniw, rw, lenrw )
    end if

    Obj = fObj

    ! If "sticky parameters no",  restore the user-defined options

    stkyOp  = iw(116)

    if (stkyOp <= 0) then
       call chcopy &
            ( 130, Usercw, 1, cw(51), 1 )
       call icopy &
            ( 130, Useriw, 1, iw(51), 1 )
       call dcopy &
            ( 130, Userrw, 1, rw(51), 1 )
    end if

    ! Print times for all clocks (if lvlTim > 0).

    call s1time( 0, 2, iw, leniw, rw, lenrw )

    return

    ! Local exit messages.

800 call snWRAP( INFO, Solver, str, str2, iw, leniw )

999 return

  end subroutine snoptBKer

  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

  subroutine s0fgctrl ( iExit, modefg, getCon, getObj, &
                        n, negCon, nnCon0, nnCon, nnJac, nnL, ngObj0, ngObj, &
                        odecon, algcon, x, ne, nlocJ, locJ, indJ, &
                        fCon, fObj, gCon, gObj, &
                        cu, lencu, iu, leniu, ru, lenru, &
                        cw, lencw, iw, leniw, rw, lenrw )
    external :: &
         odecon, algcon
    logical :: &
         getCon, getObj
    integer(ip) :: &
         iExit, lencu, leniu, lenru, lencw, leniw, lenrw, modefg, n, &
         ne, negCon, nlocJ, nnCon0, nnCon, nnJac, nnL, &
         ngObj0, ngObj, indJ(ne), locJ(nlocJ), iu(leniu), &
         iw(leniw)
    real(rp) :: &
         fObj, fCon(nnCon0), gObj(ngObj0), gCon(negCon), x(n), &
         ru(lenru), rw(lenrw)
    character(8) :: &
         cu(lencu), cw(lencw)
    ! ==========================================================================
    ! s0fgctrl  is an instance of fgwrap that calls the routines
    !   ctrlCon   and   ctrlObj   to evaluate the problem functions
    ! and possibly their gradients.
    ! odecon and algcon are user-written routines which provide
    ! the ODEs and algebraic constraints and the Jacobians.
    !
    ! Arguments  ctrlCon  and  ctrlObj  are called using modefg to
    ! control the gradients as follows:
    !
    ! modefg        Task
    ! ------        ----
    ! 2     Assign fCon, fObj and all known elements of gCon and gObj.
    ! 1     Assign all known elements of gCon and gObj.
    ! (fObj and fCon are ignored).
    ! 0     Assign fObj, fCon.  (gCon and gObj are ignored).
    !
    ! If s0fgctrl is called with minmax = 0 (feasible point only) then
    ! ngObj = max(nnJac,nnObj)  and the user objective is not used.
    !
    ! 06 Sep 2008: s0fgB adapted for control interface.
    ! 27 Aug 2011: Call-status removed from the argument list.
    ! ==========================================================================
    character &
         str*80
    external :: &
         ddot
    logical :: &
         FPonly, scaled
    integer(ip) :: &
         gotFD, gotG, gotGl, l, lG, lgConu, lgObju, lvlScl, lvlTim, &
         lvlDer, lAscal, lx0, lxScal, minmax, modeC, modeF, ngrad, &
         nfCon1, nfCon2, nfObj1, nfObj2, nnGlin, nnObj, Status
    real(rp) :: &
         ddot, Gdummy
    ! ------------------------------------------------------------------
    parameter         (lvlDer =  70) ! = 0,1,2 or 3, deriv. level
    parameter         (gotFD  = 183) ! > 0 => some differences needed
    parameter         (gotG   = 184) ! > 0 => some exact derivs
    parameter         (gotGl  = 185) ! > 0 => constant Jacob elements
    parameter         (nfCon1 = 189) ! calls to fCon: mode = 0
    parameter         (nfCon2 = 190) ! calls to fCon  mode > 0
    parameter         (nfObj1 = 194) ! calls to fObj: mode = 0
    parameter         (nfObj2 = 195) ! calls to fObj: mode > 0

    real(rp) ::   half,            one
    parameter    (half   = 0.5d+0, one   = 1.0d+0)
    ! ------------------------------------------------------------------
    nnObj     = iw( 22) ! # of objective variables
    lvlScl    = iw( 75) ! scale option
    minmax    = iw( 87) ! 1, 0, -1  => MIN, FP, MAX
    lvlTim    = iw(182) ! Timing level

    lAscal    = iw(296) ! Ascale(nb) = row and column scales
    lx0       = iw(298) ! x0(nnL)    = proximal point starting point
    lxScal    = iw(302) ! xScal(n)   = copy of scaled  x
    lgConu    = iw(319) ! record of unknown derivatives and constants
    lgObju    = iw(323) ! record of unknown derivatives

    Gdummy    = rw( 69) ! definition of an 'unset' value

    iExit     = 0

    FPonly    = minmax .eq. 0
    scaled    = lvlScl .eq. 2

    modeC     = modefg
    modeF     = modefg

    ! Determine the user-function call-status.

    call s8Stat( Status, iw, leniw )

    if (Status == 1) then
       !--------------------------------------------------------------
       ! First evaluation of the problem functions.
       ! On entry, lvlScl = 0.
       !--------------------------------------------------------------
       iw(gotFD) =  0
       iw(gotG)  =  0
       iw(gotGl) =  0
       call snPRNT( 13, ' ', iw, leniw )
       call dload ( negCon, Gdummy, gCon, 1 )
       call dload ( ngObj , Gdummy, gObj, 1 )
    end if

    !-----------------------------------------------------------------
    ! Unscale x (never required for Status = 1)
    !-----------------------------------------------------------------
    if ( scaled ) then
       call dcopy &
            ( nnL, x         , 1, rw(lxScal), 1 )
       call ddscl &
            ( nnL, rw(lAscal), 1, x         , 1 )

       ! If the Jacobian has some constant elements, they are wrecked
       ! by the scaling.  Restore them from gConu.

       if ( getCon ) then
          if (modefg > 0  .AND.  iw(gotGl) > 0) then
             call dcopy ( negCon, rw(lgConu), 1, gCon, 1 )
          end if
       end if
    end if

    !=================================================================
    ! Compute the constraint functions.
    !=================================================================
    if ( getCon ) then
       if (lvlTim >= 2) call s1time( 4, 0, iw, leniw, rw, lenrw )
       call ctrlCon ( modeC, nnCon, nnJac, negCon, &
                      x, fCon, gCon, Status,       &
                      odecon, algcon,              &
                      cu, lencu, iu, leniu, ru, lenru )
       if (lvlTim >= 2) call s1time(-4, 0, iw, leniw, rw, lenrw )
       iw(nfCon1) = iw(nfCon1) + 1
       if (modefg > 0) &
            iw(nfCon2) = iw(nfCon2) + 1
    end if

    !=================================================================
    ! Compute the objective function.
    !=================================================================
    if (getObj  .AND.  modeC >= 0) then
       if ( FPonly ) then
          call dcopy &
               ( ngObj, x, 1, gObj, 1 )
          call daxpy &
               ( ngObj, (-one), rw(lx0), 1, gObj, 1 )
          fObj = half*ddot ( ngObj, gObj, 1, gObj, 1 )
       else ! ngObj = nnObj
          if (lvlTim >= 2) call s1time( 5, 0, iw, leniw, rw, lenrw )
          call ctrlObj &
               ( modeF, ngObj, x, fObj, gObj, Status, odecon, &
               cu, lencu, iu, leniu, ru, lenru )
          if (lvlTim >= 2) call s1time(-5, 0, iw, leniw, rw, lenrw )
          iw(nfObj1) = iw(nfObj1) + 1
          if (modefg > 0) &
               iw(nfObj2) = iw(nfObj2) + 1
       end if
    end if

    !-----------------------------------------------------------------
    ! Scale  x and the derivatives.
    !-----------------------------------------------------------------
    if ( scaled ) then
       call dcopy ( nnL, rw(lxScal), 1, x, 1 )

       if ( getCon ) then
          call dddiv ( nnCon, rw(lAscal+n), 1, fCon, 1 )
          if (modefg > 0  .AND.  iw(gotG) > 0) then
             call s8sclJ &
                  ( nnCon, nnJac, negCon, n, rw(lAscal), &
                  ne, nlocJ, locJ, indJ, gCon, rw, lenrw )
          end if
       end if

       if (getObj  .AND.  modeC >= 0) then
          if (modefg > 0  .AND.  iw(gotG) > 0) then
             call s8sclg &
                  ( ngObj, rw(lAscal), gObj, rw, lenrw )
          end if
       end if
    end if

    if (modeC < 0  .OR.  modeF < 0) then
       !--------------------------------------------------------------
       ! The user may be saying the function is undefined (mode = -1)
       ! or may just want to stop                         (mode < -1).
       !--------------------------------------------------------------
       if (modeC == -1  .OR.  modeF == -1) then
          iExit = -1
       else
          if (modeC < 0) then
             iExit = 72
          else
             iExit = 73
          end if
       end if
    end if

    !=================================================================
    ! Do some housekeeping after the first call.
    !=================================================================
    if (Status == 1  .AND.  iExit == 0) then
       if ( getCon ) then
          !-----------------------------------------------------------
          ! Count how many Jacobian elements are provided.
          !-----------------------------------------------------------
          nnGlin = 0
          ngrad  = 0
          do l = 1, negCon
             if (gCon(l) /= Gdummy) ngrad  = ngrad + 1
          end do

          write(str, 1100) ngrad, negCon
          call snPRNT( 3, str, iw, leniw )

          if (ngrad < negCon) then

             ! Some Jacobian elements are missing.

             if (iw(lvlDer) >= 2) then
                !-----------------------------------------------------
                ! All the Jacobian is known.  Any undefined elements
                ! are assumed constant, and are restored from gConu.
                !-----------------------------------------------------
                call snPRNT( 3, &
                     ' ==>  Some constraint derivatives are missing, ' &
                     //' assumed constant.', iw, leniw )
                call snPRNT( 3, ' ', iw, leniw )

                lG  = lgConu
                do l  = 1, negCon
                   if (gCon(l) == Gdummy) then
                      gCon(l) = rw(lG)
                      nnGlin  = nnGlin + 1
                   end if
                   lG = lG + 1
                end do
             else
                !-----------------------------------------------------
                ! Save a permanent copy of gCon in gConu so that we know
                ! which derivatives must be estimated.
                !-----------------------------------------------------
                call dcopy ( negCon, gCon, 1, rw(lgConu), 1 )
             end if
          end if ! ngrad < negCon
          if (ngrad + nnGlin < negCon) iw(gotFD) = 1
          if (ngrad          >      0) iw(gotG ) = 1
          if (nnGlin         >      0) iw(gotGl) = 1
       end if

       if ( getObj ) then
          !-----------------------------------------------------------
          ! Count how many working gradient elements are known.
          ! (These may be the gradients of the FP objective.)
          !-----------------------------------------------------------
          ngrad  = 0
          do l = 1, ngObj
             if (gObj(l) /= Gdummy) ngrad = ngrad + 1
          end do

          if ( FPonly ) then
             write(str, 2010) ngObj
             call snPRNT( 3, str, iw, leniw )
          else
             write(str, 2000) ngrad, nnObj
             call snPRNT( 3, str, iw, leniw )
          end if

          if (ngrad < ngObj) then

             ! Some objective gradients are missing.

             if (iw(lvlDer) == 1  .OR.  iw(lvlDer) == 3) then
                !-----------------------------------------------------
                ! The objective gradient was meant to be known.
                !-----------------------------------------------------
                iw(lvlDer) = iw(lvlDer) - 1
                write(str, 2100) iw(lvlDer)
                call snPRNT( 3, str, iw, leniw )
             end if
             !--------------------------------------------------------
             ! Copy gObj into gObju.
             !--------------------------------------------------------
             call dcopy ( ngObj, gObj, 1, rw(lgObju), 1 )
          end if
          if (ngrad < ngObj) iw(gotFD) = 1
          if (ngrad >     0) iw(gotG ) = 1
       end if
    end if

    return

1100 format(' The user has defined', i8, '   out of', i8, &
         '   constraint gradients.')
2000 format(' The user has defined', i8, '   out of', i8, &
         '   objective  gradients.')
2010 format(' SnOptB  will define ', i8, '   gradients for the ', &
         ' FP objective.')
2100 format(' XXX  Some objective  derivatives are missing ---', &
         ' derivative level reduced to', i3)

  end subroutine s0fgctrl

  !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

end module ct25opt
