C     ******************************************************************
C     ******************************************************************

      subroutine newtonkkt(n,x,l,u,m,lambda,rho,equatn,linear,gtype,
     +epsact,epsfeas,epsopt,iprint,ncomp,macheps,bignum,xa,lambdaa,fa,
     +ca,snorma,nala,nalpsupna,iter,fcnt,gcnt,hcnt,dcnt,ns,accinfo,
     +inform)

      implicit none

C     PARAMETERS
      integer hnnzmax,jcnnzmax,lwi1,mmax,nmax

      parameter ( hnnzmax   =   15000000 )
      parameter ( jcnnzmax  =    5000000 )
      parameter ( lwi1      =   15000000 )
      parameter ( mmax      =     500000 )
      parameter ( nmax      =     500000 )

C     SCALAR ARGUMENTS
      integer accinfo,dcnt,fcnt,gcnt,hcnt,gtype,inform,iprint,iter,m,n,
     +        ncomp,ns
      double precision bignum,fa,epsact,epsfeas,epsopt,macheps,
     +        nalpsupna,snorma

C     ARRAY ARGUMENTS
      logical equatn(m),linear(m)
      double precision ca(m),l(n),lambda(m),lambdaa(m),nala(n),rho(m),
     +        u(n),x(n),xa(n)

C     LOCAL SCALARS
      logical scale
      integer col,divit,fun,hnnz,i,j,lin,maxfrt,maxit,minp,ninn,nsteps,
     +        var,ifail,na
      double precision nalpsupnprev,ops,snormprev,val,nalpsupn0,snorm0,
     +        snormref,nalpsuref,epsadd

C     LOCAL ARRAYS
      integer icntl(30),ikeep(3*nmax),info(20),inn(nmax),inp(mmax),
     +        ins(mmax),wi1(lwi1),wi2(2*nmax)
      double precision b(nmax+mmax),cntl(5),d(nmax+mmax),s(nmax+mmax),
     +        wd(nmax),wd2(4*(nmax+mmax))

C     COMMON ARRAYS
      integer hcol(hnnzmax),hlen(0:mmax),hlin(hnnzmax),hsta(0:hnnzmax)
      double precision hval(hnnzmax)

C     COMMON BLOCKS
      common /hessdata/ hval,hlin,hcol,hsta,hlen

#ifndef CUTEr

C     COMMON SCALARS
      logical constrc

C     COMMON ARRAYS
      integer jcvar(jcnnzmax),jcsta(mmax),jclen(mmax)
      double precision dpdc(mmax),g(nmax),jcval(jcnnzmax)

C     LOCAL SCALAR
      integer k

C     COMMON BLOCKS
      common /graddata/ g,dpdc,jcval,jcvar,jcsta,jclen,constrc

#else

C     COMMON SCALARS
      integer gnnz,jcnnz
      logical constrc

C     COMMON ARRAYS
      integer jcfun(jcnnzmax),jcvar(jcnnzmax)
      double precision dpdc(mmax),g(nmax),jcval(jcnnzmax)

C     LOCAL SCALAR
      integer flag,hlnnz

C     COMMON BLOCKS
      common /graddata/ g,dpdc,jcval,jcfun,jcvar,jcnnz,gnnz,constrc

#endif

C     ==================================================================
C     PRESENTATION
C     ==================================================================

      if ( iprint .ge. 1 ) then
          write(* ,8000)
          write(10,8000)
      end if

C     ==================================================================
C     INITIALIZE
C     ==================================================================

      maxit = 5

      scale = .true.

      epsadd = 1.0d-08

      iter  =  0
      fcnt  =  0
      gcnt  =  0
      hcnt  =  0
      dcnt  =  0
      divit =  0

C     ==================================================================
C     SET CURRENT POINT
C     ==================================================================

      na = n

      do i = 1,n
          xa(i) = max( l(i), min( x(i), u(i) ) )
      end do

      do j = 1,m
          lambdaa(j) = lambda(j)
      end do

C     ==================================================================
C     DETECT ACTIVE CONSTRAINTS AND ADD SQUARED SLACKS
C     ==================================================================

      call dact(na,xa,l,u,m,lambdaa,equatn,epsact,nmax,ns,ins,ninn,inn,
     +minp,inp,inform)

      if ( inform .lt. 0 ) then

          if ( iprint .ge. 1 ) then
              write(*, 9100) inform
              write(10,9100) inform
          end if

          return
      end if

C     ==================================================================
C     MAIN LOOP
C     ==================================================================

 100  continue

C     ==================================================================
C     COMPUTE FIRST DERIVATIVES
C     ==================================================================

C     OBJECTIVE FUNCTION AND CONSTRAINTS (JUST TO REPORT THEM)

      call evalobjcsq(na,xa,ns,m,ins,fa,ca,inform)
      fcnt = fcnt + 1

      if ( inform .lt. 0 ) then

          if ( iprint .ge. 1 ) then
              write(*, 9100) inform
              write(10,9100) inform
          end if

          return
      end if

#ifndef CUTEr

C     GRADIENT OF THE OBJECTIVE FUNCTION

      if ( gtype .eq. 0 ) then
          call evalgsq(na,xa,ns,m,ins,nala,inform)
      else
          call evalgdiff(na,xa,nala,macheps,inform)
      end if

      if ( inform .lt. 0 ) then

          if ( iprint .ge. 1 ) then
              write(*, 9100) inform
              write(10,9100) inform
          end if

          return
      end if

C     JACOBIAN OF THE CONSTRAINTS

      k = 0

      do j = 1,m

C         COMPUTE THE GRADIENT OF THE j-TH CONSTRAINT

          jcsta(j) = k + 1
          if ( gtype .eq. 0 ) then
              call evaljacsq(na,xa,ns,m,ins,j,jcvar(k+1),jcval(k+1),
     +        jclen(j),inform)
          else
              call evaljacdiff(na,xa,j,jcvar(k+1),jcval(k+1),jclen(j),
     +        macheps,inform)
          end if

          if ( inform .lt. 0 ) then

              if ( iprint .ge. 1 ) then
                  write(*, 9100) inform
                  write(10,9100) inform
              end if

              return
          end if

          k = k + jclen(j)

          do i = jcsta(j),jcsta(j) + jclen(j) - 1
              nala(jcvar(i)) = nala(jcvar(i)) + lambdaa(j) * jcval(i)
          end do

      end do

#else

C     COMPUTE THE GRADIENT OF THE OBJECTIVE FUNCTION AND THE GRADIENTS
C     OF THE CONSTRAINTS

      call evalgjc(na,xa,m,jcfun,jcvar,jcval,jcnnz,gnnz,flag)

      if ( flag .ne. 0 ) then
          inform = - 92

          if ( iprint .ge. 1 ) then
              write(*, 9100) inform
              write(10,9100) inform
          end if

          return
      end if

C     COMPUTE GRADIENT OF THE LAGRANGIAN

      do i = 1,na
          nala(i) = 0.0d0
      end do

      do i = 1,gnnz
          var       = jcvar(i)
          val       = jcval(i)
          nala(var) = nala(var) + val
      end do

      do i = gnnz + 1,jcnnz
          fun       = jcfun(i)
          var       = jcvar(i)
          val       = jcval(i)
          nala(var) = nala(var) + lambdaa(fun) * val
      end do

#endif

      gcnt = gcnt + 1

C     ==================================================================
C     ASSEMBLE KKT SYSTEM RHS AND COMPUTE ITS NORM
C     ==================================================================

      if ( iter .eq. 0 ) then
          nalpsupnprev = bignum
          snormprev    = bignum
      else
          nalpsupnprev = nalpsupna
          snormprev    = snorma
      end if

      nalpsupna = 0.0d0
      do i = 1,na
          if ( inn(i) .le. ninn ) then 
              d(inn(i)) = - nala(i)
          else
              d(inn(i)) = 0.0d0
          end if
          nalpsupna = max( nalpsupna, abs( d(inn(i)) ) ) 
      end do

C     This is true if nalpsupn is too big or NaN.
      if ( .not. nalpsupna .le. bignum ) then
          nalpsupna = bignum
      end if

      snorma = 0.0d0
      do j = 1,m
          if ( inp(j) .le. minp ) then
              d(na+inp(j)) = - ca(j)
          else
              d(na+inp(j)) = 0.0d0
          end if
          snorma = max( snorma, abs( d(na+inp(j)) ) )
      end do

C     This is true if snorm is too big or NaN.
      if ( .not. snorma .le. bignum ) then
          snorma = bignum
      end if

      if ( iter .eq. 0 ) then
          nalpsupn0 = nalpsupna
          snorm0    = snorma
      end if

C     ==================================================================
C     WRITE INFORMATION OF THE CURRENT POINT
C     ==================================================================

      if ( iprint .ge. 1 ) then
          write(* ,8010) iter,fa,snorma,nalpsupna
          write(10,8010) iter,fa,snorma,nalpsupna
      end if

C     write(*,*) 'Current point: ',(xa(i),i=1,na),(lambdaa(j),j=1,m)

C     ==================================================================
C     TEST STOPPING CRITERIA
C     ==================================================================

      if ( snorma .le. epsfeas .and. nalpsupna .le. epsopt ) then
        ! THE POINT IS FEASIBLE AND OPTIMAL

          accinfo = 0

          if ( iprint .ge. 1 ) then
              write(* ,9000)
              write(10,9000)
          end if

          go to 500
      end if

      snormref  = max(snorm0,   max(snormprev,   max(epsfeas, 1.0d+01)))
      nalpsuref = max(nalpsupn0,max(nalpsupnprev,max(epsopt,  1.0d+01)))

      if ( snorma .gt. snormref .or. nalpsupna .gt. nalpsuref ) then
        ! IT SEEMS TO BE DIVERGING

          divit = divit + 1

          if ( divit .ge. 1 ) then

              accinfo = 1

              if ( iprint .ge. 1 ) then
                  write(* ,9010)
                  write(10,9010)
              end if

              go to 500

          end if

      else
          divit = 0
      end if

      if ( iter .ge. maxit ) then
        ! MAXIMUM NUMBER OF ITERATIONS EXCEEDED

          accinfo = 2

          if ( iprint .ge. 1 ) then
              write(* ,9020)
              write(10,9020)
          end if

          go to 500
      end if

      if ( minp .gt. ninn ) then
        ! TOO MANY ``RELEVANT'' CONSTRAINTS: JACOBIAN WILL BE SINGULAR

          accinfo = 3

          if ( iprint .ge. 1 ) then
              write(* ,9030)
              write(10,9030)
          end if

          epsadd = max( 1.0d-08, epsadd * 10.0d0 )
          write(*,*) 'epsadd was increased to: ',epsadd
C         go to 500
      end if

C     ==================================================================
C     DO AN ITERATION
C     ==================================================================

      iter = iter + 1

C     ==================================================================
C     COMPUTE SECOND DERIVATIVES
C     ==================================================================

#ifndef CUTEr

C     COMPUTE HESSIAN OF THE OBJECTIVE FUNCTION

      hsta(0) = 1

      call evalhsq(na,xa,ns,m,ins,hlin,hcol,hval,hlen(0),inform)

      if ( inform .lt. 0 ) then

          if ( iprint .ge. 1 ) then
              write(*, 9100) inform
              write(10,9100) inform
          end if

          return
      end if

      k = hlen(0)

C     COMPUTE HESSIANS OF THE CONSTRAINTS

      do j = 1,m

          hsta(j) = k + 1

          call evalhcsq(na,xa,ns,m,ins,j,hlin(k+1),hcol(k+1),hval(k+1),
     +    hlen(j),inform)

          if ( inform .lt. 0 ) then

              if ( iprint .ge. 1 ) then
                  write(*, 9100) inform
                  write(10,9100) inform
              end if

              return
          end if

          k = k + hlen(j)

      end do

#else

      call evalhl(na,xa,m,lambdaa,hlin,hcol,hval,hlnnz,flag)

      if ( flag .ne. 0 ) then
          inform = - 96

          if ( iprint .ge. 1 ) then
              write(*, 9100) inform
              write(10,9100) inform
          end if

          return
      end if

#endif

      hcnt = hcnt + 1

C     ==================================================================
C     ASSEMBLE THE REDUCED JACOBIAN
C     ==================================================================

C     R'_x

C     We assume that just a triangle (lower or upper) of the symmetric 
C     Hessian of the Lagrangian was computed.

 400  continue

#ifndef CUTEr

      hnnz = 0

      do j = 0,m
          do i = hsta(j),hsta(j) + hlen(j) - 1
              if ( hlin(i) .ge. hcol(i) ) then

                  lin = inn(hlin(i))
                  col = inn(hcol(i))

                  if ( j .ne. 0 ) then
                      val = hval(i) * lambdaa(j)
                  else
                      val = hval(i)
                  end if

                  if ( lin .le. ninn .and. col .le. ninn ) then
                      if ( val .ne. 0.0d0 ) then
C                         A(lin,col) = A(lin,col) + val
                          hnnz = hnnz + 1
                          hlin(hnnz) = lin
                          hcol(hnnz) = col
                          hval(hnnz) = val
                      end if
                  end if

              end if
          end do
      end do

#else

      hnnz = 0

      do i = 1,hlnnz
          if ( hlin(i) .ge. hcol(i) ) then

              lin = inn(hlin(i))
              col = inn(hcol(i))
              val = hval(i)

              if ( lin .le. ninn .and. col .le. ninn ) then
                  if ( val .ne. 0.0d0 ) then
C                     A(lin,col) = A(lin,col) + val
                      hnnz = hnnz + 1
                      hlin(hnnz) = lin
                      hcol(hnnz) = col
                      hval(hnnz) = val
                  end if
              end if

          end if
      end do

#endif

C     R'_\lambda and R'_\mu

#ifndef CUTEr

      do j = 1,m
          do i = jcsta(j),jcsta(j) + jclen(j) - 1

              fun = inp(j)
              var = inn(jcvar(i))
              val = jcval(i)

              if ( var .le. ninn .and. fun .le. minp ) then
                  if ( val .ne. 0.0d0 ) then
C                     A(fun+ninn,var) = A(fun+ninn,var) + val
                      hnnz = hnnz + 1
                      hlin(hnnz) = fun + ninn
                      hcol(hnnz) = var
                      hval(hnnz) = val
                  end if
              end if

          end do
      end do

#else

      do i = gnnz + 1,jcnnz

          fun = inp(jcfun(i))
          var = inn(jcvar(i))
          val = jcval(i)

          if ( var .le. ninn .and. fun .le. minp ) then
              if ( val .ne. 0.0d0 ) then
C                 A(fun+ninn,var) = A(fun+ninn,var) + val
                  hnnz = hnnz + 1
                  hlin(hnnz) = fun + ninn
                  hcol(hnnz) = var
                  hval(hnnz) = val
              end if
          end if

      end do

#endif

C     Add a small number to the diagonal

      do i = 1,ninn + minp
          hlin(hnnz+i) = i
          hcol(hnnz+i) = i
          hval(hnnz+i) = epsadd
      end do

      hnnz = hnnz + ninn + minp

C     ==================================================================
C     SOLVE THE NEWTONIAN SYSTEM
C     ==================================================================

      if ( ninn + minp .ne. 0 ) then

C         SCALE MATRIX

          if ( scale ) then

              call mc30ad(ninn+minp,hnnz,hval,hlin,hcol,s,wd2,0,ifail)

              if ( ifail .ne. 0 ) then
                ! WARNING: SCALING FAILURE

                  scale = .false.

                  if ( iprint .ge. 1 ) then
                      write(* ,9090)
                      write(10,9090)
                  end if
              end if

C             write(*,*) 'Scaling factors= ',(exp(s(i)),i=1,ninn + minp)

          end if

          if ( scale ) then
              do i = 1,hnnz
                  hval(i) = hval(i) * exp( s(hlin(i)) + s(hcol(i)) )
              end do
          end if

C         INITIALIZE MA27

          call ma27id(icntl,cntl)

C         Suppress diagnosting, warning and error messages
          icntl(1) = 0
          icntl(2) = 0

C         ANALYSE SPARSITY PATTERN

          call ma27ad(ninn+minp,hnnz,hlin,hcol,wi1,lwi1,ikeep,wi2,
     +    nsteps,0,icntl,cntl,info,ops)

          if ( info(1) .ne. 0 ) then
            ! SINGULAR JACOBIAN

              accinfo = 4

              if ( iprint .ge. 1 ) then
                  write(* ,9040)
                  write(10,9040)
              end if

              epsadd = max( 1.0d-08, epsadd * 10.0d0 )
              write(*,*) 'epsadd was increased to: ',epsadd

              go to 400
C             go to 500
          end if

C         FACTORIZE MATRIX

          call ma27bd(ninn+minp,hnnz,hlin,hcol,hval,hnnzmax,wi1,lwi1,
     +    ikeep,nsteps,maxfrt,wi2,icntl,cntl,info)

          if ( info(1) .eq. -4 ) then
            ! INSUFFICIENT DOUBLE PRECISION WORKING SPACE

              accinfo = 6

              if ( iprint .ge. 1 ) then
                  write(* ,9060) hnnzmax,info(2)
                  write(10,9060) hnnzmax,info(2)
              end if

              go to 500

          else if ( info(1) .eq. -3 ) then
            ! INSUFFICIENT INTEGER WORKING SPACE

              accinfo = 7

              if ( iprint .ge. 1 ) then
                  write(* ,9070) lwi1,info(2)
                  write(10,9070) lwi1,info(2)
              end if

              go to 500

          else if ( info(1) .eq. 3 .or. info(1) .eq. -5 ) then
            ! SINGULAR JACOBIAN

              accinfo = 5

              if ( iprint .ge. 1 ) then
                  write(* ,9050)
                  write(10,9050)
              end if

              epsadd = max( 1.0d-08, epsadd * 10.0d0 )
              write(*,*) 'epsadd was increased to: ',epsadd

              go to 400
c             go to 500

C         else if ( info(15) .ne. minp ) then
C           ! WRONG INERTIA (SEE NOCEDAL AND WRIGHT)

C             Lemma 16.3 [pg. 447]: Suppose the the Jacobian of the 
C             constraints has full rank and that the reduced Hessian
C             Z^T H Z is positive definite. Then the Jacobian of the
C             KKT system has ninn positive eigenvalues, minp negative
C             eigenvalues, and no zero eigenvalues.

C             Note that at this point we know that the matrix has no 
C             zero eigenvalues. info(15) gives the number of negative
C             eigenvalues. 

C             accinfo = 8

C             if ( iprint .ge. 1 ) then
C                 write(* ,9080) ninn,minp,ninn+minp-info(15),info(15)
C                 write(10,9080) ninn,minp,ninn+minp-info(15),info(15)
C             end if

C             go to 500
          end if

C         SET RHS

          do i = 1,ninn
              b(i) = d(i)
          end do

          do j = 1,minp
              b(ninn+j) = d(na+j)
          end do

C         SOLVE THE EQUATIONS

          if ( scale ) then
              do i = 1,ninn + minp
                  b(i) = b(i) * exp( s(i) )
              end do
          end if

          call ma27cd(ninn+minp,hval,hnnzmax,wi1,lwi1,wd,maxfrt,b,wi2,
     +    nsteps,icntl,info)

          if ( scale ) then
              do i = 1,ninn + minp
                  b(i) = b(i) * exp( s(i) )
              end do
          end if

C         write(*,*),'Solution: ',(b(i),i=1,ninn+minp)

          dcnt = dcnt + 1

      end if

C     ==================================================================
C     ASSEMBLE d AND UPDATE x AND lambda
C     ==================================================================

      do i = 1,ninn
          d(i) = b(i)
      end do

      do j = 1,minp
          d(na+j) = b(ninn+j)
      end do

C     X VARIABLES 

      do i = 1,na
          xa(i) = xa(i) + d(inn(i))
      end do

C     LAMBDA VARIABLES

      do j = 1,m
          lambdaa(j) = lambdaa(j) + d(na+inp(j))
      end do

C     ==================================================================
C     ITERATE
C     ==================================================================

      go to 100

C     ==================================================================
C     END OF MAIN LOOP
C     ==================================================================

 500  continue

C     PRINT FINAL POINT INFORMATION

      if ( iprint .ge. 1 ) then
          write(* ,8015) iter,fa,snorma,nalpsupna
          write(10,8015) iter,fa,snorma,nalpsupna
      end if

C     COMPUTE FEASIBILITY AND COMPLEMENTARITY CONSIDERING ALL THE 
C     CONSTRAINTS

      snorma = 0.0d0
      do j = 1,m
          if ( equatn(j) .or. ca(j) .ge. - epsfeas ) then
              snorma = max( snorma, abs( ca(j) ) )
          else
              snorma = max( snorma, abs( lambdaa(j) ) )
          end if
      end do

      if ( iprint .ge. 1 ) then
          write(* ,8020) snorma
          write(10,8020) snorma
      end if

C     NON-EXECUTABLE STATEMENTS

 8000 format(/1X,'NEWTON-KKT scheme in action!')
 8010 format(/1X,'NEWTON-KKT Iteration',42X,' = ',5X,I6,
     +       /1X,'Objective function value',38X,' = ',1PD11.4,
     +       /1X,'Maximal violation of constraints and complementarity',
     +            10X,' = ',1PD11.4,
     +       /1X,'Sup-norm of the projected gradient of the Lagrangian',
     +            10X,' = ',1PD11.4)
 8015 format(/1X,'NEWTON-KKT Iteration (Final point)',28X,' = ',5X,I6,
     +       /1X,'Objective function value',38X,' = ',1PD11.4,
     +       /1X,'Maximal violation of constraints and complementarity',
     +            10X,' = ',1PD11.4,
     +       /1X,'Sup-norm of the projected gradient of the Lagrangian',
     +            10X,' = ',1PD11.4)
 8020 format(/1X,'Maximal violation of constraints and complementarity',
     +            10X,' = ',1PD11.4,
     +       /1X,'(considering all the constraints)')
 9000 format(/1X,'Flag of NEWTON-KKT = KKT system solved!')
 9010 format(/1X,'Flag of NEWTON-KKT = Newton seems to be diverging.')
 9020 format(/1X,'Flag of NEWTON-KKT = Maximum of iterations ',
     +           'reached.')
 9030 format(/1X,'Flag of NEWTON-KKT = Too many relevant ',
     +           'constraints, singular Jacobian.')
 9040 format(/1X,'Flag of NEWTON-KKT = Structurally singular ',
     +           'Jacobian.')
 9050 format(/1X,'Flag of NEWTON-KKT = Singular Jacobian.')
 9060 format(/1X,'Flag of NEWTON-KKT = Insufficient double ',
     +           'precision working space for the',/1X,'direct linear ',
     +           'solver ma27. Increase parameter hnnzmax from ',I8,
     +       /1X,'to at least ',I16,' or skip Newton-KKT.')
 9070 format(/1X,'Flag of NEWTON-KKT = Insufficient integer ',
     +           'working space for the',/1X,'direct linear solver ',
     +           'ma27. Increase parameter hnnzmax from ',I8,
     +       /1X,'to at least ',I16,' or skip Newton-KKT.')
C9080 format(/1X,'Flag of NEWTON-KKT = Wrong Jacobian inertia. ',
C    +       /1X,'Desired POS = ',I7,' NEG = ',I7,
C    +        1X,'Actual POS = ',I7,' NEG = ',I7)
 9090 format(/1X,'Warning: Scaling failure. It is being turned off.')
 9100 format(/1X,'Flag of NEWTON-KKT = ',I3,' Fatal Error')

      end


C     ******************************************************************
C     ******************************************************************

      subroutine dact(n,x,l,u,m,lambda,equatn,epsact,nmax,ns,ins,ninn,
     +inn,minp,inp,inform)

C     SCALAR ARGUMENTS
      integer inform,m,minp,n,ninn,nmax,ns
      double precision epsact

C     ARRAY ARGUMENTS
      logical equatn(m)
      integer inn(n),inp(m),ins(m)
      double precision l(n),lambda(m),u(n),x(n)

C     This subroutine detects active bounds, active constraints and
C     add squared slacks to the active inequalities.

C     LOCAL SCALARS
      integer i,j,flag,mnop,nnon
      double precision c

C     DETECT POTENTIALLY ACTIVE CONSTRAINTS, SET P = { INDICES OF 
C     ACTIVE CONSTRAINTS } AND ADD SQUARED SLACKS

      ns = 0
      minp = 0
      mnop = 0
      do j = 1,m
C         Compute constraint
          call evalc(n,x,j,c,flag)
          if ( flag .ne. 0 ) then
              inform = - 91
              return
          end if

C         Detect potentially active constraint

C         Active
          if ( equatn(j) .or. c .ge. - epsact ) then
              minp   = minp + 1
              inp(j) = minp

C             Add a squared slack for potentially active inequalities
              if ( equatn(j) ) then
                  ins(j) = 0

              else
                  ns = ns + 1
                  if ( n + ns .gt. nmax ) then
                      write(*,*) 'Fatal error: increase nmax.'
                      stop
                  end if

                  ins(j) = ns
                  x(n+ns) =   sqrt( 2.0d0 * max( 0.0d0, - c ) )
                  l(n+ns) = - 1.0d+20
                  u(n+ns) =   1.0d+20
              end if

C         Non-active
          else
              mnop   = mnop + 1
              inp(j) = m + 1 - mnop

              ins(j) = 0
              lambda(j) = 0.0d0
          end if
      end do

C     write(*,*) 'In P(',minp,mnop,'): ',(inp(j),j=1,m)

C     SET NEW NUMBER OF VARIABLES

      n  = n + ns

C     DETECT ACTIVE BOUNDS AND SET N = { INDICES OF FREE VARIABLES }

C     For x(i), inn(i) <= ninn means that i belongs to N and that it is
C     the inn(i)-th variable in the reduced space. On the other hand, 
C     inn(i) > ninn means that i does not belong to N and that the value
C     of x(i) will be saved in the inn(i)-th position (bigger than ninn).

      ninn = 0
      nnon = 0
      do i = 1,n
          if ( x(i) .gt. l(i) .and. x(i) .lt. u(i) ) then
              ninn   = ninn + 1
              inn(i) = ninn
          else
              nnon   = nnon + 1
              inn(i) = n + 1 - nnon
          end if
      end do

C     write(*,*) 'In N(',ninn,nnon,'): ',(inn(i),i=1,n)

      end

C     ******************************************************************
C     ******************************************************************

      subroutine evalobjcsq(n,x,ns,m,ins,f,c,inform)

C     SCALAR ARGUMENTS
      integer inform,m,n,ns
      double precision f

C     ARRAY ARGUMENTS
      integer ins(m)
      double precision c(m),x(n)

C     This subroutine adds the effect of the squared slacks to the
C     output of evalobjc subroutine.

C     LOCAL SCALARS
      integer j,no

      no = n - ns

      call evalobjc(no,x,f,m,c,inform)

      do j = 1,m
          if ( ins(j) .ne. 0 ) then
              c(j) = c(j) + 0.5d0 * x(no+ins(j)) ** 2
          end if
      end do

      end

C     ******************************************************************
C     ******************************************************************

      subroutine evalgsq(n,x,ns,m,ins,g,inform)

C     SCALAR ARGUMENTS
      integer inform,m,n,ns

C     ARRAY ARGUMENTS
      integer ins(m)
      double precision g(n),x(n)

C     This subroutine adds the effect of the squared slacks to the
C     output of evalg subroutine.

C     LOCAL SCALARS
      integer i,flag,no

      no = n - ns

      call evalg(no,x,g,flag)
      if ( flag .ne. 0 ) then
          inform = - 92
          return
      end if

      do i = 1,ns
          g(no+i) = 0.0d0
      end do

      end

C     ******************************************************************
C     ******************************************************************

      subroutine evalhsq(n,x,ns,m,ins,hlin,hcol,hval,hnnz,inform)

C     SCALAR ARGUMENTS
      integer inform,hnnz,m,n,ns

C     ARRAY ARGUMENTS
      integer hcol(*),hlin(*),ins(m)
      double precision hval(*),x(n)

C     This subroutine adds the effect of the squared slacks to the
C     output of evalh subroutine.

C     LOCAL SCALARS
      integer flag,no

      no = n - ns

      call evalh(no,x,hlin,hcol,hval,hnnz,flag)
      if ( flag .ne. 0 ) then
          inform = - 94
          return
      end if

      end

C     ******************************************************************
C     ******************************************************************

      subroutine evaljacsq(n,x,ns,m,ins,ind,jcvar,jcval,jcnnz,inform)

C     SCALAR ARGUMENTS
      integer ind,inform,jcnnz,m,n,ns

C     ARRAY ARGUMENTS
      integer ins(m),jcvar(*)
      double precision jcval(*),x(n)

C     This subroutine adds the effect of the squared slacks to the
C     output of evaljac subroutine.

C     LOCAL SCALARS
      integer flag,no

      no = n - ns

      call evaljac(no,x,ind,jcvar,jcval,jcnnz,flag)
      if ( flag .ne. 0 ) then
          inform = - 93
          return
      end if

      if ( ins(ind) .ne. 0 ) then
          jcnnz = jcnnz + 1
          jcvar(jcnnz) = no + ins(ind)
          jcval(jcnnz) = x(no+ins(ind))
      end if

      end
      
C     ******************************************************************
C     ******************************************************************

      subroutine evalhcsq(n,x,ns,m,ins,ind,hclin,hccol,hcval,hcnnz,
     +inform)

C     SCALAR ARGUMENTS
      integer hcnnz,ind,inform,m,n,ns

C     ARRAY ARGUMENTS
      integer hccol(*),hclin(*),ins(m)
      double precision hcval(*),x(n)

C     This subroutine adds the effect of the squared slacks to the
C     output of evalhc subroutine.

C     LOCAL SCALARS
      integer flag,no

      no = n - ns

      call evalhc(no,x,ind,hclin,hccol,hcval,hcnnz,flag)
      if ( flag .ne. 0 ) then
          inform = - 95
          return
      end if

      if ( ins(ind) .ne. 0 ) then
          hcnnz = hcnnz + 1
          hclin(hcnnz) = no + ins(ind)
          hccol(hcnnz) = no + ins(ind)
          hcval(hcnnz) = 1.0d0
      end if

      end
