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

      subroutine newtd(nind,x,l,u,g,m,rho,equatn,d,adsupn,maxelem,
     +memfail,inform)

      implicit none

C     SCALAR ARGUMENTS
      logical memfail
      integer inform,nind,m
      double precision adsupn,maxelem

C     ARRAY ARGUMENTS
      logical equatn(m)
      double precision d(nind),g(*),l(*),rho(m),u(*),x(*)

C     This subroutine solves the Newtonian system
C
C             ( H + rho A^T A ) x = b
C
C     by solving the Martinez-Santos system
C
C             H x + A^t y = b
C             A x - y/rho = 0

#include "dim.par"
#include "machconst.com"
#include "algparam.com"
#include "outtyp.com"
#include "itetyp.com"

C     PARAMETERS
      integer dimmax
      parameter ( dimmax = nmax + mmax )

C     COMMON ARRAYS
      double precision mindiag(nmax)

C     LOCAL SCALARS
      logical adddiff
      integer dim,hnnz,i,iter,lssinfo,nneigv,pind
      double precision diff,pval

C     LOCAL ARRAYS
      integer hdiag(dimmax),hlin(hnnzmax),hcol(hnnzmax)
      double precision adddiag(dimmax),adddiagprev(dimmax),
     +        hval(hnnzmax),sol(dimmax)

C     COMMON BLOCKS
      common /diadat/ mindiag
      save   /diadat/

C     ------------------------------------------------------------------
C     Presentation
C     ------------------------------------------------------------------

      if ( iprintinn .ge. 5 ) then
          write(* ,1000)
          write(10,1000)
      end if

C     ------------------------------------------------------------------
C     Initialization
C     ------------------------------------------------------------------

      iter    = 0
      memfail = .false.

      call lssini(sclsys,.true.,.false.)

C     ------------------------------------------------------------------
C     Compute ML matrix
C     ------------------------------------------------------------------

      call mlsyst(nind,x,g,m,rho,equatn,hlin,hcol,hval,hnnz,hdiag,sol,
     +dim,inform)

      if ( inform .lt. 0 ) return

      maxelem = 0.0d0
      do i = 1,hnnz
          maxelem = max( maxelem, abs( hval(i) ) )
      end do

C     ------------------------------------------------------------------
C     Analyse sparsity pattern
C     ------------------------------------------------------------------

      call lssana(dim,hnnz,hlin,hcol,hval,hdiag,lssinfo)

      if ( lssinfo .eq. 6 ) then
        ! INSUFFICIENT SPACE TO STORE THE LINEAR SYSTEM

          memfail = .true.
          return

      end if

C     ------------------------------------------------------------------
C     Main loop
C     ------------------------------------------------------------------

 100  continue

      iter = iter + 1

C     ------------------------------------------------------------------
C     Compute regularization
C     ------------------------------------------------------------------

      if ( iter .eq. 1 ) then

          if ( sameface .and. ittype .eq. 3 ) then

              do i = 1,nind
                  mindiag(i) = 0.1d0 * mindiag(i)
              end do

          else

              do i = 1,nind
                  if ( g(i) .eq. 0.0d0 ) then
                      mindiag(i) = 0.0d0
                  else
                      if ( g(i) .gt. 0.0d0 ) then
                          diff = x(i) - l(i)
                      else if ( g(i) .lt. 0.0d0 ) then
                          diff = u(i) - x(i)
                      end if
                      mindiag(i) = abs( g(i) / diff )
                  end if
              end do

          end if

          do i = 1,nind
              adddiag(i) = max( macheps23, mindiag(i) - hval(hdiag(i)) )
          end do

      else

          do i = 1,nind
              adddiagprev(i) = adddiag(i)
          end do

 110      continue

          do i = 1,nind
              if ( mindiag(i) .eq. 0.0d0 ) then
                  mindiag(i) = macheps23
              else
                  mindiag(i) = 10.0d0 * mindiag(i)
              end if
          end do

          do i = 1,nind
              adddiag(i) = max( macheps23, mindiag(i) - hval(hdiag(i)) )
          end do

          adddiff = .false.
          do i = 1,nind
              if ( adddiag(i) .gt. adddiagprev(i) ) then
                  adddiff = .true.
              end if
          end do

          if ( .not. adddiff ) then
              go to 110
          end if

      end if

      do i = nind + 1,dim
          adddiag(i) = 0.0d0
      end do

      adsupn = 0.0d0
      do i = 1,dim
          adsupn = max( adsupn, adddiag(i) )
      end do

      if ( iprintinn .ge. 5 ) then
          write(* ,1010) adsupn
          write(10,1010) adsupn
      end if

C     ------------------------------------------------------------------
C     Factorize matrix
C     ------------------------------------------------------------------

      call lssfac(dim,hnnz,hlin,hcol,hval,hdiag,adddiag,pind,pval,
     +nneigv,lssinfo)

      if ( lssinfo .eq. 0 .or. lssinfo .eq. 1 ) then

          if ( nneigv .ne. dim - nind ) then
C           ! WRONG INERTIA (SEE NOCEDAL AND WRIGHT)

C             Lemma 16.3 [pg. 447]: Assume that 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 n positive eigenvalues, m negative
C             eigenvalues, and no zero eigenvalues.

C             Note that at this point we know that the matrix has no
C             zero eigenvalues. nneigv gives the number of negative
C             eigenvalues.

              if ( iprintinn .ge. 5 ) then
                  write(* ,1020) nneigv,dim - nind
                  write(10,1020) nneigv,dim - nind
              end if

              go to 100

          else

              if ( iprintinn .ge. 5 ) then
                  write(* ,1030)
                  write(10,1030)
              end if

          end if

      else if ( lssinfo .eq. 2 ) then
        ! SINGULAR JACOBIAN

          if ( iprintinn .ge. 5 ) then
              write(* ,1040)
              write(10,1040)
          end if

          go to 100

      else if ( lssinfo .eq. 6 ) then
        ! INSUFFICIENT SPACE TO STORE THE LINEAR SYSTEM

          memfail = .true.
          return

      else if ( lssinfo .eq. 7 ) then
        ! INSUFFICIENT DOUBLE PRECISION WORKING SPACE

          memfail = .true.
          return

      else ! if ( lssinfo .eq. 8 ) then
        ! INSUFFICIENT INTEGER WORKING SPACE

          memfail = .true.
          return

      end if

C     ------------------------------------------------------------------
C     Solve
C     ------------------------------------------------------------------

      call lsssol(dim,sol)

      do i = 1,nind
          d(i) = sol(i)
      end do

      if ( iprintinn .ge. 5 .and. nprint .ne. 0 ) then
          write(*, 1050) min0(nind,nprint),(d(i),i=1,min0(nind,nprint))
          write(10,1050) min0(nind,nprint),(d(i),i=1,min0(nind,nprint))
      end if

C     NON-EXECUTABLE STATEMENTS

 1000 format(/,5X,'Sparse factorization of the ML system.')
 1010 format(  5X,'Maximum value added to the diagonal: ',1P,D24.16)
 1020 format(  5X,'ML-matrix with wrong inertia.',
     +       /,5X,'Actual number of negative eigenvalues  = ',I16,'.',
     +       /,5X,'Desired number of negative eigenvalues = ',I16,'.')
 1030 format(  5X,'Direct solver finished successfully.')
 1040 format(  5X,'ML-matrix numerically singular.')
 1050 format(/,5X,'Newton direction (first ',I7,' components): ',
     +       /,1(5X,6(1X,1P,D11.4)))

      end

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

      subroutine mlsyst(nind,x,nal,m,rho,equatn,ulin,ucol,uval,unnz,
     +udiag,b,dim,inform)

      implicit none

C     SCALAR ARGUMENTS
      integer dim,inform,m,nind,unnz

C     ARRAY ARGUMENTS
      logical equatn(m)
      integer ucol(*),udiag(*),ulin(*)
      double precision b(*),nal(nind),rho(m),uval(*),x(*)

#include "dim.par"
#include "graddat.com"
#include "rspace.com"

C     LOCAL SCALARS
      integer col,i,j,k,lin,var

C     LOCAL ARRAYS
      integer wi(nmax)

C     This subrotuine is called from the reduced space.

C     MATRIX

C     Compute Hessian of the Lagrangian

      do i = 1,nt - nind
          x(nind+i) = xcomplement(i)
      end do

      call expand(nind,x)

      call sevalhl(nt,x,m,dpdc,ulin,ucol,uval,unnz,inform)
      if ( inform .lt. 0 ) return

      call shrink(nind,x)

C     Preparation for shrink (wi indicates, for each free variable x_i,
C     its rank within the set of free variables. wi(i)=0 if x_i is not
C     a free variable)

      do i = 1,nt
         wi(i) = 0
      end do

      do i = 1,nind
         wi(ind(i)) = i
      end do

C     Shrink Hessian of the Lagrangian and set diagonal-elements indices

      do i = 1,nind
          udiag(i) = 0
      end do

      k = 0

      do i = 1,unnz
          lin = wi(ulin(i))
          col = wi(ucol(i))

          if ( lin .ne. 0 .and. col .ne. 0 ) then
              k = k + 1
              ulin(k) = lin
              ucol(k) = col
              uval(k) = uval(i)

              if ( lin .eq. col ) then
                  udiag(lin) = k
              end if
          end if
      end do

      do i = 1,nind
          if ( udiag(i) .eq. 0 ) then
              k = k + 1
              ulin(k) = i
              ucol(k) = i
              uval(k) = 0.0d0

              udiag(i) = k
          end if
      end do

C     Shrink Jacobian and add diagonal matrix - 1.0 / rho

      dim = nind

      do j = 1,m
          if ( equatn(j) .or. dpdc(j) .gt. 0.0d0 ) then

              dim = dim + 1

              do i = jcsta(j),jcsta(j) + jclen(j) - 1
                  var = wi(jcvar(i))

                  if ( var .ne. 0 ) then
                      k = k + 1
                      ulin(k) = dim
                      ucol(k) = var
                      uval(k) = jcval(i)
                  end if
              end do

              k = k + 1
              ulin(k) = dim
              ucol(k) = dim
              uval(k) = - 1.0d0 / rho(j)

              udiag(dim) = k
          end if
      end do

      unnz = k

C     RHS

      do i = 1,nind
          b(i) = - nal(i)
      end do

      do i = nind + 1,dim
          b(i) = 0.0d0
      end do

      end
