program algencanma C ALGENCAN Test driver. This master file uses ALGENCAN (an Augmented C Lagrangian method for inequality constraints and bounds that uses C GENCAN to solve the subproblems) to solve problems from the CUTE C collection. 67 different Augmented Lagrangian methods are implemented. C C E. G. Birgin, R. Castillo and J. M. Martínez, "Numerical comparison C of Augmented Lagrangian algorithms for nonconvex problems", C Computational Optimization and Applications 31, pp. 31-56, 2005. C PARAMETERS logical nearlyq integer nmax,mmax,mininterp,maxextrap,ncomp,cgscre,gtype,trtype, + cgmaxit,maxit,maxfc,maxitnfp,maxitngp,maxitnqmp,iprint double precision eta,delmin,lammax,lammin,theta,gamma,beta,sigma1, + sigma2,nint,next,sterel,steabs,epsrel,epsabs,infrel, + infabs,fmin,delta0,cggpnf,cgepsi,cgepsf,epsnqmp,epsgpen, + epsgpsn,epsfc,epsnfp parameter (nmax = 100000) parameter (mmax = 100000) parameter (eta = 0.9d0) parameter (delmin = 0.1d0) parameter (lammax = 1.0d+02) parameter (lammin = 1.0d-02) parameter (theta = 1.0d-06) parameter (gamma = 1.0d-04) parameter (beta = 0.5d0) parameter (sigma1 = 0.1d0) parameter (sigma2 = 0.9d0) parameter (nint = 2.0d0) parameter (next = 2.0d0) parameter (sterel = 1.0d-07) parameter (steabs = 1.0d-10) parameter (epsrel = 1.0d-10) parameter (epsabs = 1.0d-20) parameter (infrel = 1.0d+20) parameter (infabs = 1.0d+99) parameter (mininterp = 4) parameter (maxextrap = 3) parameter (ncomp = 12) parameter (cgscre = 2) parameter (gtype = 0) parameter (trtype = 0) parameter (nearlyq = .false.) parameter (fmin = -infrel) parameter (delta0 = -1.0d0) parameter (cggpnf = 1.0d-4) parameter (cgepsi = 1.0d-1) parameter (cgepsf = 1.0d-5) parameter (epsnqmp = 1.0d-4) parameter (cgmaxit = -1) parameter (maxit = 5000) parameter (maxfc = 100000) parameter (maxitnfp = maxit) parameter (maxitngp = maxit) parameter (maxitnqmp = 5) parameter (iprint = 0) parameter (epsgpen = 1.0d-4) parameter (epsgpsn = 1.0d-4) parameter (epsfc = 1.0d-4) parameter (epsnfp = 0.0d0) C LOCAL SCALARS logical efirst,lfirst,nvfirst integer i,j,n,outiter,totiter,totfcnt,totgcnt,totcgit,inform, + totspgit,tottnit,ntimes,htvtype,inform2,icomb,icomb1, + icomb2 double precision f,al,gpeucn2,gpsupn,csupn,mumenor real elatime,dum,cputim C LOCAL ARRAYS logical equatn(mmax),linear(mmax) integer w2(nmax) double precision x(nmax),l(nmax),u(nmax),c(mmax),mu(mmax), + cprev(mmax),nal(nmax),w1(10*nmax+maxitngp),xini(nmax), + lini(nmax),uini(nmax) character * 10 xnames(nmax),gnames(mmax),pname character * 1 fact,opti,crit,erro C COMMON SCALAR integer nc,m,morig,ptype,ttype double precision tconst,umuini,urhoini,urhomul,urhofrac,umumin, + umumax C COMMON ARRAYS integer rcl(mmax),rcu(mmax) double precision xs(nmax,2),cs(mmax,2),mubar(mmax), + zeta(mmax),rho(mmax),cl(mmax),cu(mmax) C COMMON BLOCKS common /constr/xs,cs,nc common /auglag/mubar,zeta,rho,m common /cbds/cl,cu,rcl,rcu,morig common /p/ ptype,ttype,tconst common /user/umuini,urhoini,urhomul,urhofrac,umumin,umumax C EXTERNAL SUBROUTINES external csetup,cnames,cfn,cputim integer i1,i2,i3,i4,i5 double precision muini(162),rhoini(162),rhomul(162),rhofrac(162), + mumin(162),mumax(162),muiniop(2),rhoiniop(3),rhomultop(3), + rhofracop(3),muminop(3),mumaxop(3) muiniop(1) = 1.0d-6 muiniop(2) = 1.0d0 rhoiniop(1) = 1.0d-3 rhoiniop(2) = 1.0d0 rhoiniop(3) = 10.0d0 rhomultop(1) = 2.0d0 rhomultop(2) = 10.0d0 rhomultop(3) = 100.0d0 rhofracop(1) = 1.0d-2 rhofracop(2) = 1.0d-1 rhofracop(3) = 0.5d0 muminop(1) = 1.0d-6 mumaxop(1) = 1.0d+20 muminop(2) = 1.0d-6 mumaxop(2) = 1.0d+6 muminop(3) = 1.0d-3 mumaxop(3) = 1.0d+3 i = 1 do i1 = 1,2 do i2 = 1,3 do i3 = 1,3 do i4 = 1,3 do i5 = 1,3 muini (i) = muiniop (i1) rhoini (i) = rhoiniop (i2) rhomul (i) = rhomultop(i3) rhofrac(i) = rhofracop(i4) mumin (i) = muminop (i5) mumax (i) = mumaxop (i5) i = i + 1 end do end do end do end do end do C SET UP SIF DATA open(10,file='OUTSDIF.d',form='formatted',status='old') rewind 10 efirst = .true. lfirst = .false. nvfirst = .false. call csetup(10,20,n,morig,xini,lini,uini,nmax,equatn,linear,mu, + cl,cu,mmax,efirst,lfirst,nvfirst) do i = 1,n if ( lini(i) .gt. uini(i) ) then write(*,*) 'Wrong bounds. Some lb is greater than the ' write(*,*) 'corresponding ub. Infeasible problem. ' stop end if end do do i = 1,morig if ( equatn(i) ) then write(*,*) 'This master file is not ready to handle ', + 'equalities' stop end if end do do i = 1,morig if ( mu(i) .ne. 0.0d0 ) then write(*,*) 'Opa! Initial estimation for the Lagrangian ', + 'multipliers given by CUTEr is not zero (as ', + 'we assumed). Zero will be used instead of ', + mu(i) end if end do m = morig do i = 1,morig if ( cl(i) .le. - infrel .and. cu(i) .ge. infrel ) then write(*,*) 'This master file is not ready to handle ', + 'unreal inequaltities of type ', + '- infrel <= c(x) <= infrel' stop else if ( cl(i) .gt. -infrel .and. cu(i) .ge. infrel ) then rcl(i) = i rcu(i) = 0 else if ( cl(i) .le. -infrel .and. cu(i) .lt. infrel ) then rcl(i) = 0 rcu(i) = i else m = m + 1 rcl(i) = i rcu(i) = m end if end do if ( m .gt. mmax ) then write(*,*) 'Increase mmax!',' m = ',m,' mmax = ',mmax stop end if call cnames(n,morig,pname,xnames,gnames) close(10) C SOLVE THE PROBLEM WITH THE DIFFERENT METHODS open(30,file='phrline.txt') open(40,file='phrline2.txt') open(50,file='param.txt') 10 continue C SET UP METHOD read(50,*) ptype,ttype,tconst,icomb1,icomb2 if ( ptype .eq. 99 ) then go to 30 end if if ( ptype .eq. 4 .and. ttype .eq. 1 ) then htvtype = 0 else htvtype = 1 end if C SOLVE IT do icomb = icomb1,icomb2 umuini = muini(icomb) urhoini = rhoini(icomb) urhomul = rhomul(icomb) urhofrac = rhofrac(icomb) umumin = mumin(icomb) umumax = mumax(icomb) nc = 0 do i = 1,n x(i) = xini(i) l(i) = lini(i) u(i) = uini(i) end do elatime = cputim(dum) call algencan(pname,n,x,l,u,m,c,mu,mubar,rho,zeta,epsfc, + maxit,maxfc,epsnqmp,maxitnqmp,al,nal,gpeucn2,gpsupn,csupn, + outiter,totiter,totfcnt,totgcnt,totcgit,totspgit, + tottnit,inform,nearlyq,epsgpen,epsgpsn,epsnfp,fmin,delta0, + cgmaxit,cgscre,cggpnf,cgepsi,cgepsf,maxitnfp,maxitngp,gtype, + htvtype,trtype,iprint,w1,w2,cprev,eta,delmin,lammax,lammin, + theta,gamma,beta,sigma1,sigma2,nint,next,mininterp,maxextrap, + ncomp,sterel,steabs,epsrel,epsabs,infrel,infabs) elatime = cputim(dum) - elatime C CHECK IF MORE THAN ONE RUN IS NECESSARY if ( inform .ge. 0 .and. csupn .le. epsfc ) then if ( elatime .lt. 0.1d0 ) then ntimes = 1000 else if ( elatime .lt. 1.0d0 ) then ntimes = 100 else if ( elatime .lt. 10.0d0 ) then ntimes = 10 else ntimes = 1 end if else ntimes = 0 end if C IF IT IS, RUN AS MANY TIMES AS REQUIRED 20 continue if ( ntimes .gt. 1 ) then elatime = cputim(dum) do j = 1,ntimes nc= 0 do i = 1,n x(i) = xini(i) end do call algencan(pname,n,x,l,u,m,c,mu,mubar,rho,zeta, + epsfc,maxit,maxfc,epsnqmp,maxitnqmp,al,nal,gpeucn2, + gpsupn,csupn,outiter,totiter,totfcnt,totgcnt, + totcgit,totspgit,tottnit,inform,nearlyq,epsgpen, + epsgpsn,epsnfp,fmin,delta0,cgmaxit,cgscre,cggpnf, + cgepsi,cgepsf,maxitnfp,maxitngp,gtype,htvtype,trtype, + iprint,w1,w2,cprev,eta,delmin,lammax,lammin,theta, + gamma,beta,sigma1,sigma2,nint,next,mininterp, + maxextrap,ncomp,sterel,steabs,epsrel,epsabs,infrel, + infabs) end do elatime = cputim(dum) - elatime if ( elatime .lt. 0.3d0 .and. ntimes .ne. 10000 ) then ntimes = 10000 go to 20 end if elatime = elatime / ntimes end if if ( inform .ge. 0 ) then call evalfc(n,x,f,m,c,inform2) end if if ( inform .lt. 0 ) then erro = 'E' fact = 'N' opti = 'N' crit = 'N' csupn = 1.0d+99 gpsupn = 1.0d+99 f = 1.0d+99 else erro = 'N' if ( csupn .le. epsfc ) then fact = 'F' else fact = 'N' end if if ( gpsupn .le. epsgpsn ) then opti = 'O' else opti = 'N' end if if ( fact .eq. 'F' .and. opti .eq. 'O' ) then crit = 'C' else crit = 'N' end if end if mumenor = 1.0d+99 do i = 1,m mumenor = min( mumenor, abs( mu(i) ) ) end do write(30,997) ptype,ttype,pname,n,m,outiter,totiter,totspgit, + tottnit,totfcnt,totgcnt,totcgit,elatime,f,csupn, + gpsupn,mumenor,inform write(*, 998) pname,n,m,ptype,ttype,elatime,f,csupn, + gpsupn,fact,opti,crit,erro write(40,998) pname,n,m,ptype,ttype,elatime,f,csupn, + gpsupn,fact,opti,crit,erro end do C ITERATE go to 10 C CLOSE FILES AND STOP 30 continue close(30) close(40) close(50) stop 997 format(i2,'-',i2,'&',a8,' (',i6,' , ',i6,') &',i2,' &',i7,' &',i7, + ' & ',i7,' &',i7,' &',i7,' &',i7,' &',f8.2,' &',1PD12.4, + ' &',1PD11.3,' &',1PD11.3,' &',1PD11.3,' &',i3) 998 format(a8,' ',I4,' ',I5,' ',I1,' ',I2,' ',F15.10,' ', + 1PD14.7,' ',1PD7.1,' ',1PD7.1,' ',a1,' ',a1,' ',a1,' ',a1) 999 format(A8,1X,I2,1X,I2,1X,18(1X,I1), + /,15X,18(1X,I1), + /,15X,18(1X,I1)) end C ****************************************************** C ****************************************************** subroutine evalfc(n,x,f,m,c,flag) C This subroutine computes the objective function and C the constraints. C C On Entry: C C n integer, C number of variables, C C x double precision x(n), C current point, C C m integer, C number of constraints, C C On Return C C f double precision, C objective function value at x, C C c double precision c(m), C constraints at x, C C flag integer C 0 means "no errors", C 1 means "some error occured in the objective function evaluation", C 2 means "some error occured in the constraints evaluation". C PARAMETERS integer mmax parameter (mmax=100000) C COMMON SCALARS integer morig C COMMON ARRAYS integer rcl(mmax),rcu(mmax) double precision cl(mmax),cu(mmax) C SCALAR ARGUMENTS double precision f integer n,m,flag C ARRAY ARGUMENTS double precision x(n),c(m) C LOCAL SCALARS integer i double precision tmp C COMMON BLOCKS common /cbds/cl,cu,rcl,rcu,morig C EXTERNAL SUBROUTINES external cfn flag = 0 call cfn(n,morig,x,f,mmax,c) do i = 1,morig tmp = c(i) if ( rcl(i) .ne. 0 ) then c(rcl(i)) = - tmp + cl(i) end if if ( rcu(i) .ne. 0 ) then c(rcu(i)) = tmp - cu(i) end if end do end C ****************************************************** C ****************************************************** subroutine evalgjc(n,x,m,jcvar,jcfun,jcval,jcnze,flag) C This subroutine computes the gradient of the objective C function and the Jacobian of the constraints. C C On Entry: C C n integer, C number of variables, C C x double precision x(n), C current point, C C m integer, C number of constraints, C C On Return C C g double precision g(n), C gradient of the objective function at x, C C jcnze integer C number of non-null Jacobian elements, C C jcvar,jcfun integer jcvar(jcnze),jcfun(jcnze) C see below, C C jcval double precision jcval(jcnze) C the value of the non-zero Jacobian position C [jcfun(i),jcvar(i)] is jcval(i). C C NOTE: jcfun(i) equal zero means objective function. C C flag integer C 0 means "no errors", C 1 means "some error occured in the gradient evaluation", C 2 means "some error occured in the Jacobian evaluation". C PARAMETERS integer mmax,jcnzemax parameter (mmax=100000,jcnzemax=100000) C COMMON SCALARS integer morig C COMMON ARRAYS integer rcl(mmax),rcu(mmax) double precision cl(mmax),cu(mmax) C SCALAR ARGUMENTS integer n,m,flag,jcnze C ARRAY ARGUMENTS integer jcvar(jcnzemax),jcfun(jcnzemax) double precision x(n),jcval(jcnzemax) C LOCAL SCALARS integer i,dummy1,jcnzeorig double precision dummy2(1) C EXTERNAL SUBROUTINES external csgr C COMMON BLOCKS common /cbds/cl,cu,rcl,rcu,morig flag = 0 call csgr(n,morig,.false.,dummy1,dummy2,x,jcnze,jcnzemax, +jcval,jcvar,jcfun) if ( 2 * jcnze .gt. jcnzemax ) then write(*,*) 'jcnzemax must be increased up to ',2 * jcnze stop end if jcnzeorig = jcnze do i = 1,jcnzeorig if ( jcfun(i) .ne. 0 ) then if ( rcl(jcfun(i)) .ne. 0 ) then if ( rcu(jcfun(i)) .ne. 0 ) then jcnze = jcnze + 1 jcfun(jcnze) = rcu(jcfun(i)) jcvar(jcnze) = jcvar(i) jcval(jcnze) = jcval(i) end if jcval(i) = - jcval(i) end if end if end do end