*     ------------------------------------------------------------------
*     This is a simple MATLAB gateway for SQOPT. It should be sufficient
*     for most problems. If your FORTRAN compiler does not understand
*     the %val construct or if you are afraid of compiler warnings
*     related to %val parameters, please use a version without the %val.
*     
*     SQOPT solves the quadratic programming problem
*
*                                              (  x  )
*     min (1/2 x^T H x + c^T x) such that l <= (     ) <= u
*      x                                       ( A x )
*
*     SYNTAX: Provided the name of your MEX-file is
*             sqopt.<OS secific ending, e.g. mexsol on Sun Solaris>,
*             you can call sqopt the following ways:
*
*     x = sqopt(H)
*     x = sqopt(H,c)
*     x = sqopt(H,c,A,l,u)
*     x = sqopt(H,c,A,l,u,x0)
*     x = sqopt(H,c,A,l,u,x0,print)
*     [x,fval] = sqopt(...)
*     [x,fval,exitflag] = sqopt(...)
*     [x,fval,exitflag,lambda] = sqopt(...)
*     [x,fval,exitflag,lambda, nS] = sqopt(...)
*
*     INPUT: at least H or c must be provided
*     H     : sparse matrix,
*             dimension may be less than the number n of variables
*     c     : explicit objective vector (dense), dimension n
*     A     : (m,n)-matrix A (sparse) of constraints
*     l     : lower bounds (dense), dimension n+m
*     u     : upper bounds (dense), dimension n+m
*     x0    : vector of initial values (dense), dimension n
*     print : print Output to file 'sqopt.out' (unit 9) if print > 0
*
*     OUTPUT: x must be provided, others are optional
*     x        : final variables and slacks, dimension n+m
*     fval     : value of the objective function
*     exitflag : result of the call of to sqopt
*     lambda   : vector of reduced costs, dimension n+m,
*                the last m entries are the Lagrange multipliers for
*                the general constraints
*     nS       : number of degrees of freedom (number of variables minus
*                the number of active constraints)
*
*     This gateway was written by
*        Jens Bastian    bastian@eas.iis.fhg.de
*     and tested and revised by
*        Marion Bastian  bastian@math.tu-dresden.de
*
*     This version dated 2000-11-20.
*     ------------------------------------------------------------------

      subroutine mexFunction(nlhs, plhs, nrhs, prhs)

      implicit none

      integer plhs(*), prhs(*)
      integer nlhs, nrhs

      integer mxCalloc, mxCreateFull, mxGetIr, mxGetJc, mxGetPr
      integer mxGetM, mxGetN, mxGetNzmax
      integer mxIsSparse

      integer ind_H, ind_c, ind_A, ind_l, ind_u, ind_x0, ind_print
      integer ind_x, ind_fval, ind_exitflag, ind_lambda, ind_nS
      parameter (ind_H  = 1,
     &           ind_c  = 2,
     &           ind_A  = 3,
     &           ind_l  = 4,
     &           ind_u  = 5,
     &           ind_x0 = 6,
     &           ind_print = 7,
     &           ind_x  = 1,
     &           ind_fval = 2,
     &           ind_exitflag = 3,
     &           ind_lambda = 4,
     &           ind_nS = 5)

      integer c_ptr, a_ptr, ha_ptr, ka_ptr, bl_ptr, bu_ptr
      integer hh_ptr, kh_ptr, ira_ptr, jca_ptr
      integer xs_ptr, fval_ptr, exitflag_ptr, nS_ptr, rc_ptr
      integer helast_ptr, hs_ptr, pi_ptr, ptr
      integer rw_ptr, iw_ptr, ru_ptr, iu_ptr

      integer n, m, lena, lenc, ncolH, iObj
      integer leniw, lenrw, leniu, lenru
      integer iPrint, iSumm, inform, maxS
      integer mincw, miniw, minrw, nS, nInf
      real*8  Obj, ObjAdd, sInf
      character*8 Prob

      integer lencw, lencu, nName
      parameter (lencw = 500,
     &           lencu = 1,
     &           nName = 1)
      character*8 cw(lencw), cu(lencu), Names(nName)

      integer lenh, na, lx, ll, lu
      real*8 temp(1)

      external userHx

*     Check input and output.

      if (nrhs .lt. 1)
     &     call mexErrMsgTxt('At least one input is required.')
      if (nlhs .lt. 1)
     &     call mexErrMsgTxt('At least one output is required.')

*     Store H in user workspace, indexing remains zero based.

      ncolH = mxGetN(prhs(ind_H))
      if (ncolH .ne. 0) then
         if (mxIsSparse(prhs(ind_H)) .eq. 1) then
            ru_ptr = mxGetPr(prhs(ind_H))
            lenh   = mxGetNzmax(prhs(ind_H))
            hh_ptr = mxGetIr(prhs(ind_H))
            kh_ptr = mxGetJc(prhs(ind_H))
            iu_ptr = mxCalloc(lenh+ncolH+2, 4)
            call mxCopyPtrToInteger4(kh_ptr+4*ncolH, %val(iu_ptr), 1)
            call mxCopyPtrToInteger4(hh_ptr, %val(iu_ptr+4), lenh)
            call mxCopyPtrToInteger4(kh_ptr, %val(iu_ptr+4*(lenh+1)),
     &                               ncolH+1)
            lenru = lenh
            leniu = lenh + ncolH + 2
         else
            call mexErrMsgTxt('H must be sparse.')
         end if
      else
         ru_ptr = mxCalloc(1, 8)
         iu_ptr = mxCalloc(1, 4)
         lenru  = 1
         leniu  = 1
      end if

*     Get pointer and dimension of c.

      lenc = 0
      if (nrhs .ge. ind_c) then
         if (mxIsSparse(prhs(ind_c)) .eq. 1)
     &        call mexErrMsgTxt('c can not be sparse.')
         lenc = mxGetM(prhs(ind_c))
         n    = max(ncolH, lenc)
         if (n .eq. 0)
     &        call mexErrMsgTxt('At least H or c must be provided.')
      end if
      if (lenc .ne. 0) then
         if (lenc .lt. n) call mexErrMsgTxt
     &        ('Dimension of c is less than dimension of H.')
         c_ptr = mxGetPr(prhs(ind_c))
      else
         c_ptr = mxCalloc(1, 8)
      end if

*     Get pointers and dimensions of A.
*     A is a dummy free row if there are no constraints or just bounds
*     on the variables.
*     NOTE: Index arrays of A are zero-based in Matlab,
*     but must be one-based for SQOPT.
      
      if (nrhs .ge. ind_A) then
         if (nrhs .lt. ind_u) call mexErrMsgTxt
     &        ('Lower and upper bounds are required together with A.')
         na = mxGetN(prhs(ind_A))
      else
         na = 0
      end if
      if (na .ne. 0) then
         if (mxIsSparse(prhs(ind_A)) .eq. 0)
     &        call mexErrMsgTxt('A must be sparse.')
         if (lenc .eq. 0) then
            if (na .lt. n)
     &           call mexErrMsgTxt('A has too few columns.')
         else
            if (na .ne. n) call mexErrMsgTxt
     &           ('Dimensions of A and c are incompatible.')
         end if
         n       = na
         m       = mxGetM(prhs(ind_A))
         lena    = mxGetNzmax(prhs(ind_A))
         a_ptr   = mxGetPr(prhs(ind_A))
         ira_ptr = mxGetIr(prhs(ind_A))
         jca_ptr = mxGetJc(prhs(ind_A))
         ha_ptr  = mxCalloc(lena, 4)
         ka_ptr  = mxCalloc(n+1, 4)
         call mxCopyPtrToInteger4(ira_ptr, %val(ha_ptr), lena)
         call mxCopyPtrToInteger4(jca_ptr, %val(ka_ptr), n+1)
         call InitA1(%val(ha_ptr), %val(ka_ptr), lena, n)
      else
         m      = 1
         lena   = 1
         a_ptr  = mxCalloc(lena, 8)
         ha_ptr = mxCalloc(lena, 4)
         ka_ptr = mxCalloc(n+1, 4)
         call InitA0(%val(a_ptr), %val(ha_ptr), %val(ka_ptr), n)
      end if

*     Get pointers and dimensions of l, u.
*     Handle the dummy free row if there is one.

      if (nrhs .eq. ind_l)
     &     call mexErrMsgTxt('Upper bound must be specified, too.')
      ll = mxGetM(prhs(ind_l))
      if (ll .ne. 0) then         
         if (ll .ne. m+n)
     &        call mexErrMsgTxt('Dimension of l is incompatible.')
         if (na .ne. 0) then
            bl_ptr = mxGetPr(prhs(ind_l))
         else
            ptr = mxCalloc(n+m, 8)
            call mxCopyPtrToReal8(bl_ptr, %val(ptr), n)
            call SetInf(%val(ptr+n), 1, -1)
            bl_ptr = ptr
         end if
      else
         bl_ptr = mxCalloc(n+m, 8)
         call SetInf(%val(bl_ptr), n+m, -1)
      end if
      lu = mxGetM(prhs(ind_u))
      if (lu .ne. 0) then
         if (lu .ne. n+m)
     &        call mexErrMsgTxt('Dimension of u is incompatible.')
         if (na .ne. 0) then
            bu_ptr = mxGetPr(prhs(ind_u))
         else
            ptr = mxCalloc(n+m, 8)
            call mxCopyPtrToReal8(bu_ptr, %val(ptr), n)
            call SetInf(%val(ptr+n), 1, 1)
            bu_ptr = ptr
         end if
      else
         bu_ptr = mxCalloc(n+m, 8)
         call SetInf(%val(bu_ptr), n+m, 1)
      end if

*     Create arrays for return values and
*     get pointers to these arrays.

      plhs(ind_x) = mxCreateFull(n+m, 1, 0)
      xs_ptr      = mxGetPr(plhs(ind_x))
      if (nlhs .ge. ind_fval) then
         plhs(ind_fval) = mxCreateFull(1, 1, 0)
         fval_ptr       = mxGetPr(plhs(ind_fval))
      end if
      if (nlhs .ge. ind_exitflag) then
         plhs(ind_exitflag) = mxCreateFull(1, 1, 0)
         exitflag_ptr       = mxGetPr(plhs(ind_exitflag))
      end if
      if (nlhs .ge. ind_lambda) then
         plhs(ind_lambda) = mxCreateFull(n+m, 1, 0)
         rc_ptr           = mxGetPr(plhs(ind_lambda))
      else
         rc_ptr = mxCalloc(n+m, 8)
      end if
      if (nlhs .ge. ind_nS) then
         plhs(ind_nS) = mxCreateFull(1, 1, 0)
         nS_ptr       = mxGetPr(plhs(ind_nS))
      end if

*     Get pointers and dimensions of x0.

      if (nrhs .ge. ind_x0) then
         lx = mxGetM(prhs(ind_x0))
         if (lx .gt. n+m)
     &        call mexErrMsgTxt('Dimension of x0 is too big.')
         call mxCopyPtrToReal8(mxGetPr(prhs(ind_x0)), %val(xs_ptr), lx)
      end if

*     Check if the user wants a PRINT file.

      if (nrhs .ge. ind_print) then
         call mxCopyPtrToReal8(mxGetPr(prhs(ind_print)), temp, 1)
         iPrint = temp(1)
         if (iPrint .gt. 0) iPrint = 9
      else
         iPrint = 0
      end if

*     Allocate needed arrays.

      helast_ptr = mxCalloc(n+m, 4)
      hs_ptr     = mxCalloc(n+m, 4)
      pi_ptr     = mxCalloc(m, 8)
      lenrw      = max(1000, 200*(m+n))
      rw_ptr     = mxCalloc(lenrw, 8)
      leniw      = max(1000, 100*(m+n))
      iw_ptr     = mxCalloc(leniw, 4)

      call InitHsHelast(%val(hs_ptr), %val(helast_ptr), n, m)

      iObj   = 0
      ObjAdd = 0.0d+0

*     Name the Problem.

      Prob = 'sqoptmex'

*     Specify some of the SQOPT files.

      iSumm  = 0
      if (iprint .gt. 0)
     &     open(unit=9, file='sqopt.out', status='unknown')

*     ------------------------------------------------------------------
*     First,  sqInit MUST be called to initialize optional parameters 
*     to their default values.
*     ------------------------------------------------------------------
 100  call sqInit(iPrint, iSumm, cw, lencw, %val(iw_ptr), leniw,
     &     %val(rw_ptr), lenrw )

*     ------------------------------------------------------------------
*     Specify options that were not set in the Specs file.
*     ------------------------------------------------------------------
      maxS   = ncolH + 1
      call sqseti( 'Superbasics Limit ', maxS, iPrint, iSumm, inform,
     &             cw, lencw, %val(iw_ptr), leniw, %val(rw_ptr), lenrw )

*     ------------------------------------------------------------------
*     Go for it, using a Cold start.
*     ------------------------------------------------------------------
      call sqopt('Cold', userHx, m, n, lena, nName, lenc, ncolH, iObj,
     &           ObjAdd, Prob, %val(a_ptr), %val(ha_ptr), %val(ka_ptr),
     &           %val(bl_ptr), %val(bu_ptr), %val(c_ptr), Names, 
     &           %val(helast_ptr), %val(hs_ptr), %val(xs_ptr),
     &           %val(pi_ptr), %val(rc_ptr), inform,
     &           mincw, miniw, minrw, nS, ninf, sinf, Obj,
     &           cu, lencu, %val(iu_ptr), leniu, %val(ru_ptr), lenru,
     &           cw, lencw, %val(iw_ptr), leniw, %val(rw_ptr), lenrw)

      if (inform .eq. 43) then
         call mxFree(iw_ptr)
         leniw  = miniw
         iw_ptr = mxCalloc(leniw, 4)
         go to 100
      end if
      if (inform .eq. 44) then
         call mxFree(rw_ptr)
         lenrw  = minrw
         rw_ptr = mxCalloc(lenrw, 8)
         go to 100
      end if
      if (inform .eq. 20) then
         call mxFree(iw_ptr)
         leniw = leniw * 5
         iw_ptr = mxCalloc(leniw, 4)
         call mxFree(rw_ptr)
         lenrw = lenrw * 5
         rw_ptr = mxCalloc(lenrw, 8)
         go to 100
      end if

      if (nlhs .ge. ind_fval) then
         temp(1) = Obj
         call mxCopyReal8ToPtr(temp, fval_ptr, 1)
      end if
      if (nlhs .ge. ind_exitflag) then
         temp(1) = inform
         call mxCopyReal8ToPtr(temp, exitflag_ptr, 1)
      end if
      if (nlhs .ge. ind_nS) then
         temp(1) = nS
         call mxCopyReal8ToPtr(temp, nS_ptr, 1)
      end if

      if (iprint .gt. 0) close(unit=9)

      return
*     end of gateway routine
      end

*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

      subroutine userHx( ncolH, x, Hx, nState, 
     &     cu, lencu, iu, leniu, ru, lenru )

      implicit none

      integer     ncolH, nState, lencu, leniu, lenru
      real*8      x(ncolH), Hx(ncolH)
      character*8 cu(lencu)
      integer     iu(leniu)
      real*8      ru(lenru)

      integer i, j, ir, jc, jc1, iir, ijc

*     ------------------------------------------------------------------
*     calculates H*x for sparse matrix H and full vector x
*     cu: unused
*     iu: used for sparse matrix information
*         iu(1)                           : nzmax
*         iu(2)...iu(nzmax+1)             : ir(*)
*         iu(nzmax+2)...iu(nzmax+ncolH+2) : jc(*)
*     ru: used for sparse matrix information
*         ru(1)...ru(nzmax) : H(*)
*     ------------------------------------------------------------------

      do 10 i=1, ncolH
         Hx(i) = 0.0
 10   continue

*     pointer into work arrays
      iir = 2
      ijc = iu(1) + 2

*     NOTE: sparse indexing is zero based in MATLAB
      do 20 j=1, ncolH
         jc  = iu(ijc+j-1) + 1
         jc1 = iu(ijc+j)
         do 30 i=jc, jc1
            ir = iu(iir+i-1) + 1
            Hx(ir) = Hx(ir) + ru(i) * x(j)
 30      continue
 20   continue

      return
*     end of userHx
      end

*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

      subroutine InitA1(ha, ka, lena, n)

      implicit none

      integer ha(*), ka(*)
      integer lena, n
      integer i

      do 10 i=1, lena
         ha(i) = ha(i) + 1
 10   continue
      do 20 i=1, n+1
         ka(i) = ka(i) + 1
 20   continue

      return
      end

*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

      subroutine InitA0(a, ha, ka, n)

      implicit none

      real*8 a(*)
      integer ha(*), ka(*)
      integer n
      integer i

      a(1)  = 0.0
      ha(1) = 1
      ka(1) = 1
      do 10 i=2, n+1
         ka(i) = 2
 10   continue

      return
      end

*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

      subroutine SetInf(b, n, sign)

      implicit none

      real*8 b(*)
      integer n, sign
      real*8 mexGetInf
      real*8 inf
      integer i

      inf = mexGetInf()
      if (sign .lt. 0) inf = -inf
      do 10 i=1, n
         b(i) = inf
 10   continue

      return
      end

*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

      subroutine InitHsHelast(hs, helast, n, m)

      implicit none

      integer hs(*), helast(*)
      integer m, n
      integer j

*     ------------------------------------------------------------------
*     Set the initial value and status of each variable.
*     For want of something better to do, make the variables xs(1:n)
*     temporarily fixed at their current values. 
*     The crash can set the rest.
*     ------------------------------------------------------------------
      do 10 j = 1, n
         hs(j) = 0
 10   continue

*     ------------------------------------------------------------------
*     Fix the column variables to be non-elastic and the row  variables 
*     to be elastic.
*     ------------------------------------------------------------------
      do 20 j = 1, n
         helast(j) = 0
 20   continue
 
      do 30 j = n+1, n+m
         helast(j) = 3
 30   continue

      return
      end
