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

C     Last update of any of the component of this module: 
C
C     October 10, 2006.

      subroutine abb(n,x,l,u,m,lambda,rho,equatn,linear,linvar,ml,meq,
     +alin,acol,aval,annz,b,epsgap,epsgpsn,iprint,ncomp,macheps,bignum,
     +f,g,gpsupn,gap,ldh,H,nnodes,inform,ind,gbest,nun,ls,ltmp,us,utmp,
     +xbest)

      implicit none

C     SCALAR ARGUMENTS
      integer annz,ldh,m,ml,meq,n,ncomp,nnodes,inform,iprint
      double precision bignum,epsgpsn,f,epsgap,gap,gpsupn,macheps

C     ARRAY ARGUMENTS
      integer alin(annz),acol(annz),ind(ml)
      logical equatn(m),linear(m),linvar(n)
      double precision aval(annz),b(ml),g(n),gbest(n),H(ldh,n),l(n),
     +        lambda(m),ls(n),ltmp(n),nun(n),rho(m),u(n),us(n),utmp(n),
     +        x(n),xbest(n)

C     This subroutines finds the 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 listtype,utype
      integer qlim,qi,qf

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

C     LOCAL SCALARS
      logical overflow,underflow,debugging,infeas,reduce
      integer braind,i,mla,annza,genlinfail
      double precision fbest,fmin,gpsupnbest,nunpsupn,larrat,ratio,un

C     LOCAL ARRAYS
      double precision fint(2)

C     COMMON BLOCKS
      common /undata/ alpha,beta,gamma,y,nsub,under,utype
      common /list/ llist,ulist,qi,qf,qlim,listtype

C     EXTERNAL FUNCTION
      double precision iinf

C     PRESENTATION
      if ( iprint .ge. 1 ) then
          write(*, 1000)
          write(10,1000)
      end if

C     SET OPTIONS

C     Underestimatore type
      utype = 'PLAIN'

C     Debugging option to test interval analysis and valid additional
C     constraints
      debugging = .false.

C     Subdomains reduction
      reduce = .true.

C     INITIALIZATION

      genlinfail = 0

      gap = 0.0d0

C     SOLVE (LOCALLY) THE ORIGINAL NON-CONVEX PROBLEM
      fmin = - 1.0d+20

      under = .false.
      call easygenlin(n,x,l,u,m,lambda,rho,equatn,linear,ml,meq,alin,
     +acol,aval,annz,b,epsgpsn,fmin,iprint,ncomp,macheps,bignum,f,g,
     +gpsupn,ldh,H,inform)
      if ( inform .ne. 0 .and. inform .ne. 1 ) then
          genlinfail = genlinfail + 1
C         write(*,*) genlinfail,nnodes,'*'
      end if
      under = .true.

      fbest = f
      gpsupnbest = gpsupn
      do i = 1,n
          xbest(i) = x(i)
          gbest(i) = g(i)
      end do

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

C     INITIALIZE BRANCH AND BOUND STRUCTURES
      listtype = 'QUEUE'
      call inilist(n)

C     PUT THE FIRST NODE IN THE LIST OF OPEN NODES
      call putnode(n,l,u,overflow)
      if ( overflow ) then
          write(*,*) 'List of nodes is full!'
          stop
      end if

      nnodes = 1

C     MAIN LOOP
 10   continue

C     GET A NODE FORM THE LIST OF OPEN NODES
      call getnode(n,ls,us,underflow)
      if ( underflow ) then
          go to 500
      end if

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

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

          if ( debugging ) then

C             TEST ADITIONAL CONSTRAINTS
              call testaddc(n,ls,us,m,lambda,rho,equatn,linear,ml,meq,
     +        alin,acol,aval,annz,b,mla,annza,ind,fbest,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,ls,us,mla,meq,alin,acol,aval,annza,b,infeas)

          if ( infeas ) then
              go to 10
          end if

      end if

C     SET INTERVAL ANALYSIS PARAMETERS AND INITIAL VALUES
      call unini(n,ls,us,m,lambda,rho,equatn,linear,inform)
      if ( inform .lt. 0 ) then
          write(*,*) 'Inipi return inform = ',inform
          stop
      end if

C     TEST FOR DEBUGGING
      if ( debugging ) then

C         TEST INTERVALAR DEFINITION OF THE PROBLEM
          call testi(n,m,ls,us,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,ls,us,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,ls,us,m,lambda,rho,equatn,linear,inform)
          if ( inform .lt. 0 ) then
              write(*,*) 'testun return inform = ',inform
              stop
          end if

      end if

C     TEST INTERVAL OBJECTIVE FUNCTION
      call evalali(n,y,m,lambda,rho,equatn,linear,fint,inform)

C     IF THE SUBDOMAIN WILL BE DISCARDED, COMPUTE THE REAL GAP
      if ( iinf(fint) .ge. fbest - epsgap ) then

          if ( iinf(fint) .lt. fbest ) then
              gap = max( gap, fbest - iinf(fint) )
          end if

C     IF THE SUBDOMAIN WILL NOT BE DISCARDED YET
      else   

C         SOLVE THE CONVEX UNDERESTIMATING PROBLEM

C         SET THE INITIAL POINT
          do i = 1,n
              x(i) = ls(i) + 0.5d0 * ( us(i) - ls(i) )
          end do

C         CALL THE SOLVER
          fmin = fbest - epsgap

          call easygenlin(n,x,ls,us,m,lambda,rho,equatn,linear,ml,meq,
     +    alin,acol,aval,annz,b,epsgpsn,fmin,iprint,ncomp,macheps,
     +    bignum,un,nun,nunpsupn,ldh,H,inform)
          if ( inform .ne. 0 .and. inform .ne. 1 .and. 
     +         inform .ne. 4 ) then
              genlinfail = genlinfail + 1
C             write(*,*) genlinfail,nnodes
          end if

C         IF THE SUBDOMAIN WILL BE DISCARDED, COMPUTE THE REAL GAP
          if ( un .ge. fbest - epsgap ) then

              if ( un .lt. fbest ) then
                  gap = max( gap, fbest - un )
              end if

C         IF IT STILL WILL NOT BE DISCARDED
          else

C             SOLVE (LOCALLY) THE ORIGINAL NON-CONVEX PROBLEM
              fmin = - 1.0d+20

              under = .false.
              call easygenlin(n,x,ls,us,m,lambda,rho,equatn,linear,ml,
     +        meq,alin,acol,aval,annz,b,epsgpsn,fmin,iprint,ncomp,
     +        macheps,bignum,f,g,gpsupn,ldh,H,inform)
              if ( inform .ne. 0 .and. inform .ne. 1 ) then
                  genlinfail = genlinfail + 1
C                 write(*,*) genlinfail,nnodes,'*'
              end if
              under = .true.

C             SAVE BEST SOLUTION
              if ( f .lt. fbest ) then

                  if ( iprint .ge. 1 ) then
                      write(* ,2000) f,un
                      write(10,2000) f,un
                  end if

                  fbest = f
                  gpsupnbest = gpsupn
                  do i = 1,n
                      xbest(i) = x(i)
                      gbest(i) = g(i)
                  end do

              end if

C             SELECT A VARIABLE TO BRANCH (LEAST REDUCED AXIS RULE)
              braind = 0
              larrat = 0.0d0
              do i = 1,n
                  if ( .not. linvar(i) ) then
                      ratio = ( us(i) - ls(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
                  write(*,*) 'There is no variable to branch!'
                  stop
              end if

C             SPLIT INTO TWO SUBDOMAINS

              do i = 1,n
                  ltmp(i) = ls(i)
                  utmp(i) = us(i)
              end do

C             SET THE SUBDOMAIN BOUNDS
              ltmp(braind) = ls(braind)
              utmp(braind) = 0.5d0 * ( ls(braind) + us(braind) )

C             ADD TO THE LIST OF OPEN NODES
              call putnode(n,ltmp,utmp,overflow)
              if ( overflow ) then
                  write(*,*) 'List of nodes is full!'
                  stop
              end if

              nnodes = nnodes + 1

C             SET THE SUBDOMAIN BOUNDS
              ltmp(braind) = 0.5d0 * ( ls(braind) + us(braind) )
              utmp(braind) = us(braind)

C             ADD TO THE LIST OF OPEN NODES
              call putnode(n,ltmp,utmp,overflow)
              if ( overflow ) then
                  write(*,*) 'List of nodes is full!'
                  stop
              end if

              nnodes = nnodes + 1

          end if

      end if

C     ITERATE
      go to 10

C     END OF MAIN LOOP

 500  continue

C     SET SOLUTION AND RETURN
      f = fbest
      gpsupn = gpsupnbest
      do i = 1,n
          x(i) = xbest(i)
          g(i) = gbest(i)
      end do

C     PRINT FINAL INFORMATION
      inform = 0
      if ( iprint .ge. 1 ) then
          write(* ,3000) inform,epsgap,gap
          write(10,3000) inform,epsgap,gap
      end if

C     NONEXECUTABLE STATEMENTS
 1000 format( /,' Entry to ABB.',/)
 1500 format(   ' F = ',1P,D24.16)
 2000 format(   ' F = ',1P,D24.16,' (un = ',1P,D24.16,')')
 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 redbds(n,l,u,ml,meq,alin,acol,aval,annz,b,infeas)

C     SCALAR ARGUMENTS
      logical infeas
      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
      logical update
      integer icase,i,ivar,k,maxup,upcnt
      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

      maxup = 1

C     MAIN LOOP

      upcnt = 0

 100  continue

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

      k = 1

 200  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 200

      end if

      if ( .not. infeas ) then

C         Check wether the bounds were modified
          update = .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
                  update = .true.
              end if
          end do

          if ( update ) then

              upcnt = upcnt + 1

C             Update the bounds
              do i = 1,n
                  l(i) = lnew(i)
                  u(i) = unew(i)
              end do

C             Do it again
              if ( upcnt .lt. maxup ) then
                  go to 100
              end if

          end if

      end if

      end