C     =================================================================
C     Module: Global bound-constraints AlphaBB solver (using GENLIN)
C     =================================================================

C     Last update of any of the component of this module: 
C
C     May 18, 2012

      subroutine abb(n,l,u,m,lambda,rho,equatn,linear,linvar,ml,meq,
     +alin,acol,aval,annz,b,fknown,epsgap,epsfeas,epsopt,sparse,laub,
     +xub,gub,gpsupnub,gap,nnodes,inform)

      implicit none

C     SCALAR ARGUMENTS
      logical sparse
      integer annz,m,ml,meq,n,nnodes,inform
      double precision epsfeas,epsopt,epsgap,fknown,gap,gpsupnub,laub

C     ARRAY ARGUMENTS
      integer alin(annz),acol(annz)
      logical equatn(m),linear(m),linvar(n)
      double precision aval(annz),b(ml),gub(n),l(n),lambda(m),rho(m),
     +        u(n),xub(n)

C     This subroutines finds a global minimizer of a bound-constrained
C     problem using a Branch and Bound strategy based on the alphaBB
C     algorithm.

C     PARAMETERS
      include "dim.f"

C     COMMON SCALARS
      logical under
      character * 5 utype

C     COMMON ARRAYS
      integer nsub(nmax)
      double precision alpha(nmax,nsubmax),beta(nmax,nsubmax),
     +        gamma(nmax,nsubmax),y(2*nmax)

C     LOCAL SCALARS
      logical overflow,underflow
      integer braind,genlinfail,i
      double precision cnorm,gammakell,larrat,lalbsb,luhalf,ratio,wkell

C     LOCAL ARRAYS
      double precision lsb(nmax),ltmp(nmax),usb(nmax),utmp(nmax)

C     COMMON BLOCKS
      common /undata/ alpha,beta,gamma,y,nsub,under,utype

C     ==================================================================
C     PRESENTATION AND OPTIONS
C     ==================================================================

      write(*, 1000)
      write(10,1000)

C     Underestimator type
      utype = 'PLAIN'

C     ==================================================================
C     INITIALIZATION
C     ==================================================================

C     Number of failures of the local solver

      genlinfail = 0

C     Upper bound initial value

      laub = 1.0d+20

C     Lower and upper bounds for the original box

      do i = 1,n
          ltmp(i) = l(i)
          utmp(i) = u(i)
      end do

      call procsubbox(n,ltmp,utmp,m,lambda,rho,equatn,linear,ml,meq,
     +alin,acol,aval,annz,b,fknown,epsopt,sparse,lalbsb,laub,xub,gub,
     +gpsupnub,wkell,gammakell,cnorm,genlinfail,inform)

C     Add original box to S

      nnodes = 0

      call iniset()

      if ( lalbsb .le. laub ) then
          call putsubbox(n,ltmp,utmp,lalbsb,overflow)
          if ( overflow ) then
              write(*,*) 'Subboxes set S is full!'
              stop
          end if

          nnodes = nnodes + 1
      else
          gap = 1.0d+20
          write(* ,*) 'SC: Original box discarded (Infeas subproblem)!'
          write(10,*) 'SC: Original box discarded (Infeas subproblem)!'
          go to 500
      end if

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

 10   continue

C     ==================================================================
C     GET A NODE FROM THE LIST OF OPEN NODES
C     ==================================================================

      call getsubbox(n,lsb,usb,lalbsb,underflow)
      if ( underflow ) then
          gap = 0.0d0
          write(* ,*) 'SC satisfied: Empty suboxes set S (gap=0)!'
          write(10,*) 'SC satisfied: Empty suboxes set S (gap=0)!'
          go to 500
      end if

      gap = laub - lalbsb

      do i = 1,n
          if ( lsb(i) .gt. usb(i) ) then
              write(*,*) 'Wrong subbox get from S'
              write(*,*) i,lsb(i),usb(i)
              stop
          end if
      end do

C     ==================================================================
C     STOPPING CRITERION
C     ==================================================================

c     if ( gap .le. epsgap ) then
      if ( gap .le. wkell .or. ( cnorm .le. epsfeas .and. 
     +     gammakell + gap .le. epsopt ) ) then

          write(* ,*) 'SC satisfied: gap = ',gap
          write(10,*) 'SC satisfied: gap = ',gap

          go to 500
      end if

C     ==================================================================
C     SELECT A VARIABLE TO BRANCH (LEAST REDUCED AXIS RULE) AND SPLIT 
C     THE SUBBOX INTO TWO SUBBOXES 
C     ==================================================================

C     SELECT A VARIABLE TO BRANCH

      braind = 0
      larrat = 0.0d0
      do i = 1,n
          if ( .not. linvar(i) ) then
              ratio = ( usb(i) - lsb(i) ) / ( u(i) - l(i) )
              if ( ratio .gt. larrat ) then
                  braind = i
                  larrat = ratio
              end if
          end if
      end do

      if ( braind .eq. 0 ) then
          braind = 0
          larrat = 0.0d0
          do i = 1,n
              if ( .not. linvar(i) ) then
                  luhalf = 0.5d0 * ( lsb(i) + usb(i) )
                  if ( luhalf .gt. lsb(i) .and. 
     +                 luhalf .lt. usb(i) ) then
                      ratio = usb(i) - lsb(i)
                      if ( ratio .gt. larrat ) then
                          braind = i
                          larrat = ratio
                      end if
                  end if
              end if
          end do

      end if

      if ( braind .eq. 0 ) then
          write(*,*) 'There is no variable to branch!'
          write(*,*) gap,laub,cnorm
          write(*,*) 'Subbox: '
          write(*,*) 'lsb: ',(lsb(i),i=1,n)
          write(*,*) 'usb: ',(usb(i),i=1,n)
          write(*,*) 'lalbsb: ',lalbsb
          stop
      end if

C     SPLIT SUBBOX INTO TWO SUBBOXES
 
C     SUBBOX 1

      do i = 1,n
          ltmp(i) = lsb(i)
          utmp(i) = usb(i)
      end do

      utmp(braind) = 0.5d0 * ( lsb(braind) + usb(braind) )

      call procsubbox(n,ltmp,utmp,m,lambda,rho,equatn,linear,ml,meq,
     +alin,acol,aval,annz,b,fknown,epsopt,sparse,lalbsb,laub,xub,gub,
     +gpsupnub,wkell,gammakell,cnorm,genlinfail,inform)

      if ( lalbsb .lt. laub ) then
          call putsubbox(n,ltmp,utmp,lalbsb,overflow)
          if ( overflow ) then
              write(*,*) 'Subboxes set S is full!'
              stop
          end if
      
          do i = 1,n
              if ( ltmp(i) .gt. utmp(i) ) then
                  write(*,*) 'Wrong subbox put in S'
                  write(*,*) i,ltmp(i),utmp(i)
                  stop
              end if
          end do

          nnodes = nnodes + 1
      end if

C     SUBBOX 2

      do i = 1,n
          ltmp(i) = lsb(i)
          utmp(i) = usb(i)
      end do

      ltmp(braind) = 0.5d0 * ( lsb(braind) + usb(braind) )

      call procsubbox(n,ltmp,utmp,m,lambda,rho,equatn,linear,ml,meq,
     +alin,acol,aval,annz,b,fknown,epsopt,sparse,lalbsb,laub,xub,gub,
     +gpsupnub,wkell,gammakell,cnorm,genlinfail,inform)

      if ( lalbsb .lt. laub ) then
          call putsubbox(n,ltmp,utmp,lalbsb,overflow)
          if ( overflow ) then
              write(*,*) 'Subboxes set S is full!'
              stop
          end if
      
          do i = 1,n
              if ( ltmp(i) .gt. utmp(i) ) then
                  write(*,*) 'Wrong subbox put in S'
                  write(*,*) i,ltmp(i),utmp(i)
                  stop
              end if
          end do

          nnodes = nnodes + 1
      end if

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

      go to 10

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

 500  continue

C     PRINT FINAL INFORMATION

      inform = 0

      write(* ,3000) inform,epsgap,gap
      write(10,3000) inform,epsgap,gap

C     ==================================================================
C     NONEXECUTABLE STATEMENTS
C     ==================================================================

 1000 format( /,' Entry to ABB.',/)
 3000 format( /,' Flag of ABB = ',I2,
     +          ' (Convergence with an absolute gap between',
     +          ' the lower and',/,' upper bounds smaller than ',
     +            1PD11.4,'. In fact, the real gap is ',1PD11.4,'.)')

      end

C     ******************************************************************
C     ******************************************************************
      subroutine procsubbox(n,l,u,m,lambda,rho,equatn,linear,ml,meq,
     +alin,acol,aval,annz,b,fknown,epsopt,sparse,lalbsb,laub,xub,gub,
     +gpsupnub,wkell,gammakell,cnorm,genlinfail,inform)

      implicit none

C     SCALAR ARGUMENTS
      logical sparse
      integer annz,genlinfail,inform,m,meq,ml,n
      double precision cnorm,epsopt,fknown,gammakell,gpsupnub,lalbsb,
     +        laub,wkell

C     ARRAY ARGUMENTS
      logical equatn(m),linear(m)
      integer alin(annz),acol(annz)
      double precision aval(annz),b(ml),gub(n),l(n),lambda(m),rho(m),
     +        u(n),xub(n)

C     PARAMETERS
      include "dim.f"

C     COMMON SCALARS
      logical under
      character * 5 utype

C     COMMON ARRAYS
      integer nsub(nmax)
      double precision alpha(nmax,nsubmax),beta(nmax,nsubmax),
     +        gamma(nmax,nsubmax),y(2*nmax)

C     LOCAL SCALARS
      logical debugging,infeas,updtd
      integer annza,flag,i,icase,iprint,j,maxupdtd,mla,ncomp,ninf,
     +        updtdcnt
      double precision bignum,cbesti,epsgpsn,fmin,gpsupn,lalbsbtmp,
     +        laubtmp,macheps,sinf,wkella,wkellb,wkellc

C     LOCAL ARRAYS
      integer ind(mlmax)
      double precision cint(2),fint(2),g(nmax),H(nmax,nmax),x(nmax)

C     COMMON BLOCKS
      common /undata/ alpha,beta,gamma,y,nsub,under,utype

C     EXTERNAL FUNCTION
      logical irlei
      double precision iinf

C     Debugging option to test interval analysis and valid additional
C     constraints

      debugging = .false.

C     ==================================================================
C     REDUCE BOUNDS: MINIMIZING AND MAXIMIZING EACH X(i) SUBJECT TO 
C     THE LINEAR CONSTRAINTS, TRY TO REDUCE THE CURRENT SUBBOX OR 
C     EVEN DETECT INFEASIBILITY
C     ==================================================================

      maxupdtd =  1
      updtdcnt =  0

 100  continue 

C     ADD VALID LINEAR CONSTRAINTS (THAT DEPEND ON L AND U)
      call addc(n,l,u,laub,ml,meq,alin,acol,aval,annz,b,fknown,mla,
     +annza,ind)

      if ( debugging ) then

C         TEST ADITIONAL CONSTRAINTS
          call testaddc(n,l,u,m,lambda,rho,equatn,linear,ml,meq,alin,
     +    acol,aval,annz,b,mla,annza,ind,laub,inform)
          if ( inform .ne. 0 ) then
              write(*,*) 'Testaddc return inform = ',inform
              stop
          end if

      end if

C     SOLVE THE 2N LP'S OR DETECT INFEASIBILITY
      call redbds(n,l,u,mla,meq,alin,acol,aval,annza,b,infeas,updtd)

      if ( infeas ) then
c         write(*,*) 'salio por null reduced subbox en (1)'
          lalbsb = 1.0d+20
          return
      end if

      if ( updtd ) then
          updtdcnt = updtdcnt + 1
          if ( updtdcnt .lt. maxupdtd ) then
              go to 100
          end if
      end if

      do i = 1,n
          if ( l(i) .gt. u(i) ) then
              write(*,*) 'Wrongly updated bounds in redbds'
              write(*,*) i,l(i),u(i)
              stop
          end if
      end do

C     ==================================================================
C     CHECK FEASIBILITY OF ORIGINAL CONSTRAINTS WITHIN THE REDUCED SUBBOX
C     ==================================================================

C     SET INTERVAL ANALYSIS PARAMETERS AND INITIAL VALUES

C     (At this point, y receives the interval representation of 
C      the subbox [l,u] that is being used below)

      call unini(n,l,u,m,lambda,rho,equatn,linear,inform)
      if ( inform .lt. 0 ) then
          write(*,*) 'unini return inform = ',inform
          stop
      end if

C     TEST FOR DEBUGGING
      if ( debugging ) then

C         TEST INTERVALAR DEFINITION OF THE PROBLEM
          call testi(n,m,l,u,inform)
          if ( inform .lt. 0 ) then
              write(*,*) 'Testi return inform = ',inform
              stop
          end if

C         TEST INTERVAL AUGMENTED LAGRANGIAN AND ITS HESSIAN
          call testali(n,l,u,m,lambda,rho,equatn,linear,inform)
          if ( inform .lt. 0 ) then 
              write(*,*) 'testhali return inform = ',inform
              stop
          end if

C         TEST UNDERESTIMATOR CONVEXITY 
          call testun(n,l,u,m,lambda,rho,equatn,linear,inform)
          if ( inform .lt. 0 ) then
              write(*,*) 'testun return inform = ',inform
              stop
          end if

      end if

C     CHECK ORIIGINAL CONSTRAINTS

      infeas = .false.

      j = 1
 50   if ( j .le. m .and. .not. infeas ) then
          call evalci(n,y,j,cint,flag)
          if ( flag .ne. 0 ) then
              inform = - 91
              return
          end if

          if ( ( equatn(j) .and. .not. irlei(0.0d0,cint) ) .or.
     +         ( .not. equatn(j) .and. iinf(cint) .gt. 0.0d0 ) ) then
              infeas = .true.
          end if

          j = j + 1
          go to 50
      end if

      if ( infeas ) then
c         write(*,*) 'salio por infeas en (2)'
          lalbsb = 1.0d+20
          return
      end if

C     CHECK ORIIGINAL OBJECTIVE FUNCTION

      call evalfi(n,y,fint,flag)
      if ( flag .ne. 0 ) then
          inform = - 90
          return
      end if

C     fknown represents the smallest value of f at a known feasible
C     points. It is assumed fknown is a parameter (given by pre-
C     processing local minimizations of the original problem) and no 
C     effort is done to update fknown in case other feasbible points 
C     are discovered on the fly. If no feasible points is known, it
C     is assumed that fknown = 1.0d+20.

      if ( fknown .ne. 1.0d+20 .and. fknown .lt. iinf(fint) ) then
c         write(*,*) 'salio por interval analysis de la f en (2b)'
          lalbsb = 1.0d+20
          return
      end if

C     ==================================================================
C     COMPUTE INTERVAL AUGMENTED LAGRANGIAN WITHIN THE REDUCED SUBBOX
C     ==================================================================

      call evalali(n,y,m,lambda,rho,equatn,linear,fint,inform)
      lalbsb = iinf(fint)

      if ( laub .le. lalbsb ) then
c         write(*,*) 'salio por interval analysis del LA en (3)'
          return
      end if

C     ==================================================================
C     SOLVE THE CONVEX UNDERESTIMATING PROBLEM WITHIN THE SUBREGION 
C     DETERMINED BY THE REDUCED SUBBOX PLUS THE ADDITIONAL CONSTRAINTS
C     ==================================================================

C     SET THE INITIAL POINT

      do i = 1,n
          x(i) = 0.5d0 * ( l(i) + u(i) )
      end do

C     CALL THE SOLVER

      if ( .not. sparse ) then

          fmin    = - 1.0d+20
          epsgpsn =   epsopt
          iprint  =   0
          ncomp   =   6
          macheps =   1.0d-16
          bignum  =   1.0d+99

          call easygenlin(n,x,l,u,m,lambda,rho,equatn,linear,ml,meq,
     +    alin,acol,aval,annz,b,epsgpsn,fmin,iprint,ncomp,macheps,
     +    bignum,lalbsbtmp,g,gpsupn,nmax,H,flag)

          if ( flag .ne. 0 .and. flag .ne. 1 .and. flag .ne. 4 ) then
              genlinfail = genlinfail + 1
          end if

      else

          fmin    = - 1.0d+20
          epsgpsn =   epsopt

          call easyminos2(n,x,l,u,m,lambda,rho,equatn,linear,ml,meq,
     +    alin,acol,aval,annz,b,epsgpsn,fmin,lalbsbtmp,g,gpsupn,icase,
     +    ninf,sinf,flag)

      end if

      lalbsb = max( lalbsb, lalbsbtmp )

      if ( laub .le. lalbsb ) then
c         write(*,*) 'salio por alpha-convex underestimator en (4)'
          return
      end if

C     ==================================================================
C     PERFORM LOCAL MINIMIZATION OF THE NON-CONVEX AUGMENTED LAGRANGIAN
C     SUBPROBLEM WITHIN THE SUBREGION DETERMINED BY THE REDUCED SUBBOX  
C     PLUS THE ADDITIONAL CONSTRAINTS
C     ==================================================================

C     LOCAL MINIMIZATION

      under = .false.

      if ( .not. sparse ) then

          fmin    = - 1.0d+20
          epsgpsn =   epsopt
          iprint  =   0
          ncomp   =   6
          macheps =   1.0d-16
          bignum  =   1.0d+99

          call easygenlin(n,x,l,u,m,lambda,rho,equatn,linear,ml,meq,
     +    alin,acol,aval,annz,b,epsgpsn,fmin,iprint,ncomp,macheps,
     +    bignum,laubtmp,g,gpsupn,nmax,H,flag)

          if ( flag .ne. 0 .and. flag .ne. 1 ) then
              genlinfail = genlinfail + 1
          end if

      else

          fmin    = - 1.0d+20
          epsgpsn =   epsopt

          call easyminos2(n,x,l,u,m,lambda,rho,equatn,linear,ml,meq,
     +    alin,acol,aval,annz,b,epsgpsn,fmin,laubtmp,g,gpsupn,icase,
     +    ninf,sinf,flag)

      end if

c     write(*,*) 'hizo local minimization y encontro f = ',laubtmp

      under = .true.

C     CHECK IMPROVEMENT IN THE INCUMBENT SOLUTION

      if ( laubtmp .lt. laub ) then
c         write(*,*) 'Upper bound improved from ',laub,' to ',laubtmp

          laub = laubtmp
          gpsupnub = gpsupn
          do i = 1,n
              xub(i) = x(i)
              gub(i) = g(i)
          end do

          wkella = 0.0d0
          wkellb = 0.0d0
          wkellc = 0.0d0
          gammakell = 0.0d0
          do i = 1,m
              call evalc(n,xub,i,cbesti,flag)

              if ( equatn(i) ) then
                  wkella = wkella + cbesti ** 2
                  gammakell = gammakell + ( lambda(i) / rho(1) ) ** 2
     +                      - ( lambda(i) / rho(1) + cbesti ) ** 2
              else
                  wkellb = wkellb + max( 0.0d0, cbesti ) ** 2
                  wkellc = wkellc + abs( min( - cbesti, lambda(i) / 
     +                     rho(1) ) )
                  gammakell = gammakell + ( lambda(i) / rho(1) ) ** 2
     +                      - max( 0.0d0, ( lambda(i) / rho(1) + 
     +                        cbesti ) ) ** 2
              end if
          end do
          cnorm = sqrt( wkella ) + sqrt( wkellb )
          wkell = sqrt( wkella ) + sqrt( wkellb ) + wkellc
          gammakell = 0.5d0 * rho(1) * gammakell

          if ( iprint .ge. 1 ) then
              write(* ,2000) laub,wkell,gammakell,cnorm
              write(10,2000) laub,wkell,gammakell,cnorm
          end if
      end if

C     ==================================================================
C     RETURN
C     ==================================================================

      return

C     ==================================================================
C     NON-EXECUTABLE STATEMENTS
C     ==================================================================

 2000 format(   ' LAub = ',1P,D24.16,' w = ',1P,D7.1,
     +          ' gamma = ',1P,D8.1,' feas = ',1P,D7.1)

      end

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

      subroutine redbds(n,l,u,ml,meq,alin,acol,aval,annz,b,infeas,updtd)

C     SCALAR ARGUMENTS
      logical infeas,updtd
      integer annz,meq,ml,n

C     ARRAY ARGUMENTS
      integer alin(annz),acol(annz)
      double precision aval(annz),b(ml),l(n),u(n)

C     PARAMETER
      include "dim.f"

C     COMMON SCALARS
      logical sparse

C     LOCAL SCALARS
      integer icase,i,ivar,k
      double precision ctx,tmp
      character * 5 bdstype

C     LOCAL ARRAYS
      double precision lnew(nmax),unew(nmax)

C     COMMON BLOCKS
      common /innsol/ sparse

      if ( sparse ) then
          call redbds2(n,l,u,ml,meq,alin,acol,aval,annz,b,infeas,updtd)
          return
      end if

      infeas = .false.

      if ( ml .eq. 0 ) then
          return
      end if

C     WHILE INFEASIBILITY IS NOT DETECTED REDUCE THE LOWER AND UPPER 
C     BOUND OF EACH VARIABLE

      k = 1

 100  if ( k .le. 2 * n .and. .not. infeas ) then

          if ( k .le. n ) then
              bdstype = 'LOWER'
              ivar = k
          else
              bdstype = 'UPPER'
              ivar = k - n
          end if

C         SOLVE AN LP USING SIMPLEX

          call easysimplex(n,l,u,ml,meq,alin,acol,aval,annz,b,bdstype,
     +    ivar,ctx,icase)

C         ANALIZE SIMPLEX OUTPUT

C         Solution found
          if ( icase .eq. 0 ) then
              if ( bdstype .eq. 'LOWER' ) then
                  lnew(ivar) = min( u(ivar), max( l(ivar), - ctx ) )
              else ! if ( bdstype .eq. 'UPPER' ) then
                  unew(ivar) = max( l(ivar), min( u(ivar),   ctx ) )
                  if ( unew(ivar) .lt. lnew(ivar) ) then
                      tmp        = lnew(ivar)
                      lnew(ivar) = unew(ivar)
                      unew(ivar) = tmp
                  end if
              end if

C         Infeasible
          else if ( icase .eq. - 1 ) then
              infeas = .true.

C         Unbounded
          else
              write(*,*) 'Unbounded LP !!!'
              write(*,*) 'All variables must have bounds !!!'
              stop
          end if

          k = k + 1
          go to 100

      end if

C     Check wether the bounds were modified
      if ( .not. infeas ) then
          updtd = .false.
          do i = 1,n
              if ( abs( l(i) - lnew(i) ) .gt. 1.0d-08 .or.
     +             abs( u(i) - unew(i) ) .gt. 1.0d-08 ) then
                  updtd = .true.
              end if
              l(i) = lnew(i)
              u(i) = unew(i)
          end do
      end if

      end

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

      subroutine redbds2(n,l,u,ml,meq,alin,acol,aval,annz,b,infeas,
     +updtd)

C     SCALAR ARGUMENTS
      logical infeas,updtd
      integer annz,meq,ml,n

C     ARRAY ARGUMENTS
      integer alin(annz),acol(annz)
      double precision aval(annz),b(ml),l(n),u(n)

C     PARAMETER
      include "dim.f"

C     LOCAL SCALARS
      integer icase,i,ivar,k
      double precision ctx,tmp
      character * 5 bdstype

C     LOCAL ARRAYS
      double precision lnew(nmax),unew(nmax)

      infeas = .false.

      if ( ml .eq. 0 ) then
          return
      end if

C     WHILE INFEASIBILITY IS NOT DETECTED REDUCE THE LOWER AND UPPER 
C     BOUND OF EACH VARIABLE

      k = 1

 100  if ( k .le. 2 * n .and. .not. infeas ) then

          if ( k .le. n ) then
              bdstype = 'LOWER'
              ivar = k
          else
              bdstype = 'UPPER'
              ivar = k - n
          end if

C         SOLVE AN LP USING SIMPLEX

c         call easysimplex(n,l,u,ml,meq,alin,acol,aval,annz,b,bdstype,
c    +    ivar,ctx,icase)

          call easyminos(n,l,u,ml,meq,alin,acol,aval,annz,b,bdstype,
     +    ivar,ctx,icase)

C         ANALIZE SIMPLEX OUTPUT

C         Solution found
          if ( icase .eq. 0 ) then
              if ( bdstype .eq. 'LOWER' ) then
c                 lnew(ivar) = min( u(ivar), max( l(ivar), - ctx ) )
                  lnew(ivar) = min( u(ivar), max( l(ivar),   ctx ) )
              else ! if ( bdstype .eq. 'UPPER' ) then
                  unew(ivar) = max( l(ivar), min( u(ivar),   ctx ) )
                  if ( unew(ivar) .lt. lnew(ivar) ) then
                      tmp        = lnew(ivar)
                      lnew(ivar) = unew(ivar)
                      unew(ivar) = tmp
                  end if
              end if

C         Infeasible
c         else if ( icase .eq. - 1 ) then
          else if ( icase .eq.   1 ) then
              infeas = .true.

C         Unbounded
c         else
c             write(*,*) 'Unbounded LP !!!'
c             write(*,*) 'All variables must have bounds !!!'
c             stop
c         end if

C         Other cases: leave the box as it is.
          else
              return
          end if

          k = k + 1
          go to 100

      end if

C     Check wether the bounds were modified
      if ( .not. infeas ) then
          updtd = .false.
          do i = 1,n
              if ( abs( l(i) - lnew(i) ) .gt. 1.0d-08 .or.
     +             abs( u(i) - unew(i) ) .gt. 1.0d-08 ) then
                  updtd = .true.
              end if
              l(i) = lnew(i)
              u(i) = unew(i)
          end do
      end if

      end

