!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
! File ctmex.F90
! Mex function for snctrl.
!
! 11 Feb 2010: 32/64-bit Matlab gateway routine for SNCTRL.
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

#include "fintrf.h"

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

  subroutine mexFunction ( nlhs, plhs, nrhs, prhs )
    use precision,   only : ip, rp
    use control,     only : sncInit, snctrlD, sncSpec
    use ctmexModule
    implicit none

    integer(4) :: nlhs, nrhs
    mwPointer  :: prhs(*), plhs(*)
    !===========================================================================
    ! Mex function for snctrl.
    !
    ! 13 Jan 2009: First version of mexFunction for snctrl.
    ! 09 Feb 2010: 64-bit version.
    !===========================================================================
    ! Matlab
    mwPointer   :: mxGetM, mxGetN, mxGetPr, mxCreateNumericMatrix
    integer(ip) :: mxClassIDFromClassName, clsInt, clsDbl

    ! Temporary Real Variables
    real(rp)          :: rOpt, rStart, r1, r2, rnY, rnU, rnP, rnC, rnPhs
    real(rp), pointer :: robjL(:), rytype(:,:), rctype(:,:), rnpInt(:)


    ! SNOPT/SNCTRL Local Variables
    integer(ip)   :: Start, INFO, mincw, miniw, minrw, nS, nInf, &
                     iOpt, iPrt, iSum, iSpc, nY, nU, nP, nC, nPhs, s, strlen
    real(rp)      :: sInf
    character(8)  :: probname
    character(120):: filename
    integer(ip),  pointer :: ctype(:,:), hs(:)
    real(rp),     pointer :: x(:)
    external      :: mxodecon, mxalgcon, mxvarbds

    !---------------------------------------------------------------------------

    if ( nrhs < 1 ) &
         call mexErrMsgTxt ( 'Need at least 1 input argument.' )

    clsDbl = mxClassIDFromClassName('double')
    clsInt = mxClassIDFromClassName('int32')

    ! Get option.
    call mxCopyPtrToReal8 ( mxGetPr(prhs(1)), rOpt, 1 )
    iOpt = int(rOpt,ip)


    if ( iOpt == 1 ) then
       iSum = 6

       strlen = mxGetM(prhs(2)) * mxGetN(prhs(2))
       if ( strlen > 0 ) then
          iPrt = 9
          call mxGetString ( prhs(2), filename, strlen )
          open ( iPrt, file=filename, status='UNKNOWN' )
       else
          iPrt = 0

       end if

       call sncInit ( iPrt, iSum, prob, cw, lencw, iw, leniw, rw, lenrw )

    else if ( iOpt == 2 ) then
       iSpc = 4

       strlen = mxGetM(prhs(2)) * mxGetN(prhs(2))
       if ( strlen > 0 ) then
          call mxGetString ( prhs(2), filename, strlen )
          open ( iSpc, file=filename, status='old' )
          call sncSpec ( iSpc, INFO, cw, lencw, iw, leniw, rw, lenrw )

          if ( nlhs == 1 ) then
             plhs(1) = mxCreateNumericMatrix( 1, 1, clsInt, 0 )
             call mxCopyInteger4ToPtr ( int(INFO,4), mxGetPr(plhs(1)), 1 )
          end if

          close ( iSpc )

       else
          call mexErrMsgTxt ( 'Specs filename is invalid.' )

       end if

    else if ( iOpt == 3 ) then
       call mxCopyPtrToReal8 ( mxGetPr(prhs(2)),   rStart, 1 )

       strlen = mxGetM(prhs(3)) * mxGetN(prhs(3))
       call mxGetString      ( prhs(3), probname, strlen )

       call mxCopyPtrToReal8 ( mxGetPr(prhs(4)),   rnY, 1 )
       call mxCopyPtrToReal8 ( mxGetPr(prhs(5)),   rnU, 1 )
       call mxCopyPtrToReal8 ( mxGetPr(prhs(6)),   rnP, 1 )
       call mxCopyPtrToReal8 ( mxGetPr(prhs(7)),   rnC, 1 )
       call mxCopyPtrToReal8 ( mxGetPr(prhs(8)), rnPhs, 1 )

       nY   = int(rnY,ip)
       nU   = int(rnU,ip)
       nP   = int(rnP,ip)
       nC   = int(rnC,ip)
       nPhs = int(rnPhs,ip)

       prob%probName = probname
       prob%nY = nY
       prob%nU = nU
       prob%nP = nP
       prob%nC = nC
       prob%nPhs = nPhs

       allocate ( robjL(nPhs), rnpInt(nPhs), prob%phsPt(nPhs+1), rytype(nY,nPhs) )

       call mxCopyPtrToReal8 ( mxGetPr(prhs(9)),       robjL, nPhs )
       call mxCopyPtrToReal8 ( mxGetPr(prhs(10)),     rnpInt, nPhs )
       call mxCopyPtrToReal8 ( mxGetPr(prhs(11)), prob%phsPt, nPhs+1 )
       call mxCopyPtrToReal8 ( mxGetPr(prhs(12)),     rytype, nPhs*nY )

       allocate ( prob%objL(nPhs), prob%npInt(nPhs), prob%ytype(nY,nPhs) )
       prob%objL  = int(robjL,ip)
       prob%npInt = int(rnpInt,ip)
       prob%ytype = int(rytype,ip)

       deallocate ( robjL, rnpInt, rytype )

       if ( nC > 0 ) then
          allocate ( rctype(nC,nPhs) )
          call mxCopyPtrToReal8 ( mxGetPr(prhs(13)), rctype, nPhs*nC )

          allocate ( prob%ctype(nC,nPhs) )
          prob%ctype = int(rctype,ip)
          deallocate ( rctype )

       else
          allocate ( prob%ctype(1,nPhs) )
          prob%ctype = 0

       end if

       strlen = mxGetM(prhs(14)) * mxGetN(prhs(14))
       call mxGetString ( prhs(14), varbds, strlen )
       strlen = mxGetM(prhs(15)) * mxGetN(prhs(15))
       call mxGetString ( prhs(15), odecon, strlen )

       if ( nC > 0 ) then
          strlen = mxGetM(prhs(16)) * mxGetN(prhs(16))
          call mxGetString ( prhs(16), algcon, strlen )
       end if

       call snctrlD ( Start, prob, x, hs, mxodecon, mxalgcon, mxvarbds, &
                      INFO, mincw, miniw, minrw, nS, nInf, sInf, &
                      cu, lencu, iu, leniu, ru, lenru, &
                      cw, lencw, iw, leniw, rw, lenrw )

       ! Return output:
       if ( nlhs == 2 ) then
          s = size(x)
          plhs(1) = mxCreateNumericMatrix(s, 1, clsDbl, 0)
          plhs(2) = mxCreateNumericMatrix(s, 1, clsInt, 0)
          call mxCopyReal8ToPtr ( x, mxGetPr(plhs(1)), s )
          call mxCopyInteger4ToPtr ( int(hs,4), mxGetPr(plhs(2)), s )

       end if

       deallocate ( prob%objL, prob%ytype, prob%npInt, prob%phsPt, prob%ctype )
       deallocate ( x, hs )

    else if ( 4 <= iOpt .and. iOpt <= 6 ) then
       call mxsncSets ( iOpt, nlhs, plhs, nrhs, prhs )

    else if ( 7 <= iOpt .and. iOpt <= 9 ) then
       call mxsncGets ( iOpt, nlhs, plhs, nrhs, prhs )

    end if

  end subroutine mexFunction

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

  subroutine mxsncSets ( iOpt, nlhs, plhs, nrhs, prhs )
    use precision, only : ip, rp
    use control,   only : sncSet, sncSeti, sncSetr
    use ctmexModule
    implicit none

    integer(4)  :: nlhs, nrhs
    integer(ip) :: iOpt
    mwPointer   :: plhs(*), prhs(*)
    !===========================================================================
    mwPointer     :: mxGetM, mxGetN, mxGetPr
    integer(ip)   :: strlen, ivalue, Errors
    real(rp)      :: rvalue
    character(72) :: buffer

    strlen = mxGetM(prhs(2)) * mxGetN(prhs(2))
    call mxGetString ( prhs(2), buffer, strlen )

    if ( iOpt == 4 ) then
       call sncSet ( buffer, 0_ip, 0_ip, Errors, cw, lencw, iw, leniw, rw, lenrw )

    else if ( iOpt == 5 ) then
       call mxCopyPtrToReal8 ( mxGetPr(prhs(3)), rvalue, 1 )
       ivalue = int(rvalue,ip)

       call sncSeti ( buffer, ivalue, 0_ip, 0_ip, Errors, &
                      cw, lencw, iw, leniw, rw, lenrw )

    else if (iOpt == 6 ) then
       call mxCopyPtrToReal8 ( mxGetPr(prhs(3)), rvalue, 1 )
       call sncSetr ( buffer, rvalue, 0_ip, 0_ip, Errors, &
                      cw, lencw, iw, leniw, rw, lenrw )

    end if

  end subroutine mxsncSets

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

  subroutine mxsncGets ( iOpt, nlhs, plhs, nrhs, prhs )
    use precision, only : ip, rp
    use control,   only : sncGet, sncGetc, sncGeti, sncGetr
    use ctmexModule
    implicit none

    integer(4)  :: nlhs, nrhs
    integer(ip) :: iOpt
    mwPointer   :: plhs(*), prhs(*)
    !===========================================================================
    mwPointer     :: mxGetM, mxGetN, mxGetPr, mxCreateNumericMatrix, &
                     mxCreateCharArray
    integer(ip)   :: mxClassIDFromClassName, clsInt, clsDbl
    integer(ip)   :: strlen, ivalue, Errors
    real(rp)      :: rvalue
    character(8)  :: cvalue
    character(72) :: buffer


    clsDbl = mxClassIDFromClassName('double')
    clsInt = mxClassIDFromClassName('int32')

    strlen = mxGetM(prhs(2)) * mxGetN(prhs(2))
    call mxGetString ( prhs(2), buffer, strlen )


    if ( iOpt == 7 ) then
       ivalue = sncGet ( buffer, Errors, cw, lencw, iw, leniw, rw, lenrw )

       if ( nlhs == 1 ) then
          plhs(1) = mxCreateNumericMatrix ( 1, 1, clsInt, 0 )
          call mxCopyInteger4ToPtr ( int(ivalue,4), mxGetPr(plhs(1)), 1 )
       end if

    else if ( iOpt == 8 ) then
       call sncGetc ( buffer, cvalue, Errors, cw, lencw, iw, leniw, rw, lenrw )

       if ( nlhs == 1 ) then
          plhs(1) = mxCreateCharArray ( 1, 8 )
          call mxCopyCharacterToPtr ( cvalue, mxGetPr(plhs(1)), 8 )
       end if

    else if ( iOpt == 9 ) then
       call sncGeti ( buffer, ivalue, Errors, cw, lencw, iw, leniw, rw, lenrw )

       if ( nlhs == 1 ) then
          plhs(1) = mxCreateNumericMatrix ( 1, 1, clsInt, 0 )
          call mxCopyInteger4ToPtr ( int(ivalue,4), mxGetPr(plhs(1)), 1 )
       end if

    else if ( iOpt == 10 ) then
       call sncGetr ( buffer, rvalue, Errors, cw, lencw, iw, leniw, rw, lenrw )

       if ( nlhs == 1 ) then
          plhs(1) = mxCreateNumericMatrix ( 1, 1, clsInt, 0 )
          call mxCopyReal8ToPtr ( rvalue, mxGetPr(plhs(1)), 1 )
       end if

    end if

  end subroutine mxsncGets

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

  subroutine mxvarbds ( curPhs, nPhs, nY, nU, nP, nC, nNodes, &
                        lbds, ubds, x, plbds, pubds, p, clbds, cubds )
    use precision,   only : ip, rp
    use ctmexModule, only : varbds
    implicit none
    integer(ip) :: curPhs, nPhs, nY, nU, nP, nC, nNodes
    real(rp)    :: lbds(nY+nU,nNodes), ubds(nY+nU,nNodes), x(nY+nU,nNodes), &
                   plbds(nP), pubds(nP), p(nP), clbds(nC), cubds(nC)
    !===========================================================================
    ! Wrapper for varbds for the mex interface.
    ! Does the Matlab callback to varbds.m to get the variable bounds.
    !
    ! 13 Jan 2009: First version of mxvarbds.
    ! 09 Feb 2010: v5.
    !===========================================================================
    integer(4)  :: nlhs, nrhs
    mwPointer   :: mxGetPr, mxCreateNumericMatrix
    integer(ip) :: mxClassIDFromClassName, clsInt
    mwPointer   :: plhs(8), prhs(7)

    nlhs = 3
    nrhs = 7
    if ( nP > 0 ) &
         nlhs = nlhs + 3
    if ( nC > 0 ) &
         nlhs = nlhs + 2

    clsInt = mxClassIDFromClassName ( 'int32' )

    prhs(1) = mxCreateNumericMatrix ( 1, 1, clsInt, 0 )
    prhs(2) = mxCreateNumericMatrix ( 1, 1, clsInt, 0 )
    prhs(3) = mxCreateNumericMatrix ( 1, 1, clsInt, 0 )
    prhs(4) = mxCreateNumericMatrix ( 1, 1, clsInt, 0 )
    prhs(5) = mxCreateNumericMatrix ( 1, 1, clsInt, 0 )
    prhs(6) = mxCreateNumericMatrix ( 1, 1, clsInt, 0 )
    prhs(7) = mxCreateNumericMatrix ( 1, 1, clsInt, 0 )

    call mxCopyInteger4ToPtr ( int(curPhs,4), mxGetPr(prhs(1)), 1 )
    call mxCopyInteger4ToPtr (   int(nPhs,4), mxGetPr(prhs(2)), 1 )
    call mxCopyInteger4ToPtr (     int(nY,4), mxGetPr(prhs(3)), 1 )
    call mxCopyInteger4ToPtr (     int(nU,4), mxGetPr(prhs(4)), 1 )
    call mxCopyInteger4ToPtr (     int(nP,4), mxGetPr(prhs(5)), 1 )
    call mxCopyInteger4ToPtr (     int(nC,4), mxGetPr(prhs(6)), 1 )
    call mxCopyInteger4ToPtr ( int(nNodes,4), mxGetPr(prhs(7)), 1 )

    call mexCallMATLAB ( nlhs, plhs, nrhs, prhs, varbds )

    call mxCopyPtrToReal8 ( mxGetPr(plhs(1)), lbds, (nY+nU)*nNodes )
    call mxCopyPtrToReal8 ( mxGetPr(plhs(2)), ubds, (nY+nU)*nNodes )
    call mxCopyPtrToReal8 ( mxGetPr(plhs(3)),    x, (nY+nU)*nNodes )

    if ( nP > 0 ) then
       call mxCopyPtrToReal8 ( mxGetPr(plhs(4)), plbds, nP )
       call mxCopyPtrToReal8 ( mxGetPr(plhs(5)), pubds, nP )
       call mxCopyPtrToReal8 ( mxGetPr(plhs(6)),     p, nP )
       if ( nC > 0 ) then
          call mxCopyPtrToReal8 ( mxGetPr(plhs(7)), clbds, nC )
          call mxCopyPtrToReal8 ( mxGetPr(plhs(8)), cubds, nC )
       end if

    else
       if ( nC > 0 ) then
          call mxCopyPtrToReal8 ( mxGetPr(plhs(4)), clbds, nC )
          call mxCopyPtrToReal8 ( mxGetPr(plhs(5)), cubds, nC )
       end if

    end if

  end subroutine mxvarbds

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

  subroutine mxodecon ( snStat, curPhs, nPhs, nY, nU, nP, nNodes, F, J, &
                        dvar, pvar, needF, needJ, cu, lencu, iu, leniu, ru, lenru )
    use precision,   only : ip, rp
    use ctmexModule, only : odecon
    implicit none
    integer(ip)  :: snStat, curPhs, nPhs, nY, nU, nP, nNodes, &
                    needF, needJ, lencu, leniu, lenru, iu(leniu)
    real(rp)     :: F(nY,nNodes), J(nY,nY+nU+nP,nNodes), &
                    dvar(nY+nU,nNodes), pvar(nP), ru(lenru)
    character(8) :: cu(lencu)
    !===================================================================
    ! Wrapper for odecon for the mex interface.
    ! Does the Matlab callback to odecon.m to get state equations and
    ! derivatives.
    !
    ! 13 Jan 2009: First version of mx_odecon.
    ! 09 Feb 2010: v5.
    !===================================================================
    integer(4), parameter :: nlhs = 2, nrhs = 11
    mwPointer   :: mxGetPr, mxCreateNumericMatrix
    integer(ip) :: mxClassIDFromClassName, clsInt, clsDbl
    mwPointer   :: plhs(2), prhs(11)


    clsInt = mxClassIDFromClassName('int32')
    clsDbl = mxClassIDFromClassName('double')

    prhs(1) = mxCreateNumericMatrix ( 1, 1, clsInt, 0 )
    prhs(2) = mxCreateNumericMatrix ( 1, 1, clsInt, 0 )
    prhs(3) = mxCreateNumericMatrix ( 1, 1, clsInt, 0 )
    prhs(4) = mxCreateNumericMatrix ( 1, 1, clsInt, 0 )
    prhs(5) = mxCreateNumericMatrix ( 1, 1, clsInt, 0 )
    prhs(6) = mxCreateNumericMatrix ( 1, 1, clsInt, 0 )
    prhs(7) = mxCreateNumericMatrix ( 1, 1, clsInt, 0 )
    prhs(8) = mxCreateNumericMatrix ( nY+nU,nNodes, clsDbl, 0 )
    prhs(9) = mxCreateNumericMatrix ( nP, 1, clsDbl, 0 )
    prhs(10) = mxCreateNumericMatrix ( 1, 1, clsInt, 0 )
    prhs(11) = mxCreateNumericMatrix ( 1, 1, clsInt, 0 )

    call mxCopyInteger4ToPtr ( int(snStat,4), mxGetPr(prhs(1)), 1 )
    call mxCopyInteger4ToPtr ( int(curPhs,4), mxGetPr(prhs(2)), 1 )
    call mxCopyInteger4ToPtr (   int(nPhs,4), mxGetPr(prhs(3)), 1 )
    call mxCopyInteger4ToPtr (     int(nY,4), mxGetPr(prhs(4)), 1 )
    call mxCopyInteger4ToPtr (     int(nU,4), mxGetPr(prhs(5)), 1 )
    call mxCopyInteger4ToPtr (     int(nP,4), mxGetPr(prhs(6)), 1 )
    call mxCopyInteger4ToPtr ( int(nNodes,4), mxGetPr(prhs(7)), 1 )
    call mxCopyReal8ToPtr    (          dvar, mxGetPr(prhs(8)), (nY+nU)*nNodes )
    call mxCopyReal8ToPtr    (          pvar, mxGetPr(prhs(9)), nP )
    call mxCopyInteger4ToPtr (  int(needF,4), mxGetPr(prhs(10)), 1 )
    call mxCopyInteger4ToPtr (  int(needJ,4), mxGetPr(prhs(11)), 1 )

    call mexCallMATLAB ( nlhs, plhs, nrhs, prhs, odecon )

    if ( needF > 0 ) then
       call mxCopyPtrToReal8 ( mxGetPr(plhs(1)), F, nY*nNodes )
    end if

    if ( needJ > 0 ) then
       call mxCopyPtrToReal8 ( mxGetPr(plhs(2)), J, nY*(nY+nU+nP)*nNodes )
    end if

  end subroutine mxodecon

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

  subroutine mxalgcon ( snStat, curPhs, nPhs, nC, nY, nU, nP, nNodes, C, G, &
                        dvar, pvar, needC, needG, cu, lencu, iu, leniu, ru, lenru )
    use precision,   only : ip, rp
    use ctmexModule, only : algcon
    implicit none
    integer(ip)  :: snStat, curPhs, nPhs, nC, nY, nU, nP, nNodes, &
                    needC, needG, lencu, leniu, lenru, iu(leniu)
    real(rp)     :: C(nC,nNodes), G(nC,nY+nU+nP,nNodes), &
                    dvar(nY+nU,nNodes), pvar(nP), ru(lenru)
    character(8) :: cu(lencu)
    !===================================================================
    ! Wrapper for algbds for the mex interface.
    ! Does the Matlab callback to algbds.m to get algebraic constraints
    ! and derivatives.
    !
    ! 13 Jan 2009: First version of mx_algbds.
    ! 09 Feb 2010: v5.
    !===================================================================
    integer(4), parameter :: nlhs = 2, nrhs = 12
    mwPointer   :: mxGetPr, mxCreateNumericMatrix
    integer(ip) :: mxClassIDFromClassName, clsInt, clsDbl
    mwPointer   :: plhs(2), prhs(12)


    clsInt = mxClassIDFromClassName ( 'int32' )
    clsDbl = mxClassIDFromClassName ( 'double' )

    prhs(1) = mxCreateNumericMatrix ( 1, 1, clsInt, 0 )
    prhs(2) = mxCreateNumericMatrix ( 1, 1, clsInt, 0 )
    prhs(3) = mxCreateNumericMatrix ( 1, 1, clsInt, 0 )
    prhs(4) = mxCreateNumericMatrix ( 1, 1, clsInt, 0 )
    prhs(5) = mxCreateNumericMatrix ( 1, 1, clsInt, 0 )
    prhs(6) = mxCreateNumericMatrix ( 1, 1, clsInt, 0 )
    prhs(7) = mxCreateNumericMatrix ( 1, 1, clsInt, 0 )
    prhs(8) = mxCreateNumericMatrix ( 1, 1, clsInt, 0 )
    prhs(9) = mxCreateNumericMatrix ( nY+nU,nNodes, clsDbl, 0 )
    prhs(10) = mxCreateNumericMatrix ( nP, 1, clsDbl, 0 )
    prhs(11) = mxCreateNumericMatrix ( 1, 1, clsInt, 0 )
    prhs(12) = mxCreateNumericMatrix ( 1, 1, clsInt, 0 )

    call mxCopyInteger4ToPtr ( int(snStat,4), mxGetPr(prhs(1)), 1 )
    call mxCopyInteger4ToPtr ( int(curPhs,4), mxGetPr(prhs(2)), 1 )
    call mxCopyInteger4ToPtr (   int(nPhs,4), mxGetPr(prhs(3)), 1 )
    call mxCopyInteger4ToPtr (     int(nC,4), mxGetPr(prhs(4)), 1 )
    call mxCopyInteger4ToPtr (     int(nY,4), mxGetPr(prhs(5)), 1 )
    call mxCopyInteger4ToPtr (     int(nU,4), mxGetPr(prhs(6)), 1 )
    call mxCopyInteger4ToPtr (     int(nP,4), mxGetPr(prhs(7)), 1 )
    call mxCopyInteger4ToPtr ( int(nNodes,4), mxGetPr(prhs(8)), 1 )
    call mxCopyReal8ToPtr    (          dvar, mxGetPr(prhs(9)), (nY+nU)*nNodes )
    call mxCopyReal8ToPtr    (          pvar, mxGetPr(prhs(10)), nP )
    call mxCopyInteger4ToPtr (  int(needC,4), mxGetPr(prhs(11)), 1 )
    call mxCopyInteger4ToPtr (  int(needG,4), mxGetPr(prhs(12)), 1 )

    call mexCallMATLAB ( nlhs, plhs, nrhs, prhs, algcon )

    if ( needC > 0 ) then
       call mxCopyPtrToReal8 ( mxGetPr(plhs(1)), C, nC*nNodes )
    end if

    if ( needG > 0 ) then
       call mxCopyPtrToReal8 ( mxGetPr(plhs(2)), G, nC*(nY+nU+nP)*nNodes )
    end if

  end subroutine mxalgcon

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