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