C ****************************************************************** C ****************************************************************** subroutine algencan(epsfeas,epsopt,efacc,eoacc,iprint,ncomp,n,x,l, +u,m,lambda,equatn,linear,coded,checkder,fu,cnormu,snorm,nlpsupn, +inform,ifile) implicit none C SCALAR ARGUMENTS logical checkder integer inform,iprint,m,n,ncomp double precision cnormu,efacc,eoacc,epsfeas,epsopt,fu,nlpsupn, + snorm character*(*) ifile C ARRAY ARGUMENTS logical coded(11),equatn(m),linear(m) double precision l(n),lambda(m),u(n),x(n) include "dim.par" include "machconst.inc" include "algconst.par" include "counters.inc" include "outtyp.inc" include "algparam.inc" include "scaling.inc" include "slacks.inc" include "fixvar.inc" C PARAMETERS integer n0,n1,n2 parameter ( n0 = 500 ) parameter ( n1 = 10000 ) parameter ( n2 = 20000 ) C LOCAL SCALARS logical innfail,lss,scl integer alinfo,fcsubt,geninfo,i,iter,j,maxit,nwcalls,nwtotit, + oprint,outiter,solinfo,totiter,msqcalls,msqtotit,toteccnt, + totejccnt,totehccnt,aveeccnt,aveejccnt,aveehccnt double precision cnorm,cnormb,cnormub,f,fb,fub,ncsupn,nlpsupnb, + rsupn real time C LOCAL ARRAYS character * 4 lsssub,sclsub double precision nl(nmax),c(mmax),cu(mmax),rho(mmax) real dum(2) C DATA STATEMENTS data dum/0.0,0.0/ C EXTERNAL FUNCTIONS AND SUBROUTINES external fparam,checkd,auglag,gencan,lss,scl external evalf,evalg,evalh,evalc,evaljac,evalhc,evalfc,evalgjac, + evalgjacp,evalhl,evalhlp C ================================================================== C Start timing C ================================================================== time = dtime(dum) C ================================================================== C Set safety mode C ================================================================== safemode = .false. C ================================================================== C Initialization C ================================================================== C Set machine-dependent constants bignum = 1.0d+99 macheps = 1.0d-16 macheps12 = sqrt( macheps ) macheps13 = macheps ** ( 1.0d0 / 3.0d0 ) macheps23 = macheps ** ( 2.0d0 / 3.0d0 ) C Set global counters fcnt = 0 efcnt = 0 efccnt = 0 egcnt = 0 egjccnt = 0 egjcpcnt = 0 ehcnt = 0 ehlcnt = 0 ehlpcnt = 0 do j = 1,m eccnt(j) = 0 ejccnt(j) = 0 ehccnt(j) = 0 end do C ================================================================== C Check user provided subroutines C ================================================================== C Set user-provided subroutines indicators fcoded = coded(1) gcoded = coded(2) hcoded = coded(3) ccoded = coded(4) jaccoded = coded(5) hccoded = coded(6) fccoded = coded(7) gjaccoded = coded(8) gjacpcoded = coded(9) hlcoded = coded(10) hlpcoded = coded(11) C Check whether mandatory subroutines are being properly provided. C For unconstrained and bound-constrained problems, EVALF must be C coded by the user. For constrained problems, EVALF and EVALC, or, C alternatively, EVALFC must be coded. (Note that EVALF/EVALC and C EVALFC should not be provided concurrently.) For feasibility C problems, a constant null objective function must be coded and the C problem solved with the IGNORE-OBJECTIVE-FUNCTION keyword. Coded C subroutines must be indicated by setting the entrances of array C named CODED within subroutine INIP. C Moreover, to avoid odd combinations, only the following choices C will be considered valid: C C If the objective function and the constraints are given by evalf C and evalc, respectively, first derivatives must be given by evalg C and evaljac, while second derivatives must be given by evalh and C evalhc. C C If the objective function and the constraints are given by evalfc C then first derivatives may be given by evalgjac or evalgjacp. If C first derivatives are given by evalgjac, second derivatives must C be given by evalhl. On the other hand, if first derivatives are C given by evalgjacp, second derivatives must be given by evalhlp. C C Any other odd combination will be ignored. C C Variables being set below have the following meaning: C C firstde: whether, following the rules dictated above, the user- C provided subroutines will allow the method to compute first C derivatives. C C seconde: whether, following the rules dictated above, the user- C provided subroutines will allow the method to compute the Hessian C matrix of the augmented Lagrangian (to minimize the augmented C Lagrangian subproblems using a second-order method) and/or the C Hessian of the Lagrangian plus the Jacobian of the constraints C (in order to solve the KKT system by Newton's method). C C truehpr: whether, following the rules dictated above, the user- C provided subroutines will allow the method to compute the C product of the true Hessian of the augmented Lagrangian times C a given vector. It would allow the method to solve the augmented C Lagrangian subproblems using a truncated-Newton method using C Conjugate Gradients to solve the Newtonian linear systems. fcsubt = 0 if ( fcoded .and. ( ccoded .or. m .eq. 0 ) ) then fcsubt = 1 firstde = .false. seconde = .false. truehpr = .false. if ( gcoded .and. ( jaccoded .or. m .eq. 0 ) ) then firstde = .true. if ( hcoded .and. ( hccoded .or. m .eq. 0 ) ) then seconde = .true. truehpr = .true. end if end if else if ( fccoded ) then fcsubt = 2 firstde = .false. seconde = .false. truehpr = .false. if ( gjaccoded ) then firstde = .true. if ( hlcoded ) then seconde = .true. truehpr = .true. end if else if ( gjacpcoded ) then firstde = .true. if ( hlpcoded ) then truehpr = .true. end if end if end if C ================================================================== C Set default values for algoritmic parameters C ================================================================== C Hessian-vector product strategy: HAPPRO, INCQUO or TRUEHP. TRUEHP C is the default option. If the proper subroutines were not coded by C the user, then HAPPRO is used instead. HAPPRO reduces to INCQUO in C the unconstrained and bound-constrained cases and switches between C INCQUO and the product by an approximation of the Hessian in the C constrained case. if ( truehpr ) then hptype = 'TRUEHP' else hptype = 'HAPPRO' end if C Inner solver method: TR (trust regions, i.e. BETRA), NW (Newtonian C system with the use of a direct solver, i.e. MA57AD), TN + TRUEHP C (Truncated Newton and true Hessian-vector product) or TN + HAPPRO C (Truncated Newton and product with Hessian approximation). If the C proper subroutines were not provided by the user then TN + HAPPRO C is used. if ( seconde .and. lss(lsssub) .and. n .le. n0 ) then innslvr = 'TR' else if ( seconde .and. lss(lsssub) .and. n .le. n1 ) then innslvr = 'NW' else innslvr = 'TN' if ( n .gt. n2 ) then hptype = 'HAPPRO' end if end if C Scaling of linear systems if ( lsssub .eq. 'MA57' .or. scl(sclsub) ) then sclsys = .true. else sclsys = .false. end if C Ignore objective function (to only find a feasible point by C minimizing 1/2 of the squared infeasibility) ignoref = .false. C Acceleration step if ( seconde .and. lss(lsssub) ) then skipacc = .false. else skipacc = .true. end if C Slacks for inequality constraints slacks = .false. C Remove fixed variables (with identical lower and upper bounds) rmfixv = .true. C Scale objective function and constraints if ( m .gt. 0 .and. .not. ignoref ) then scale = .true. else scale = .false. end if C Flags used to make inner calls to Gencan innercall = .false. useustp = .false. C ================================================================== C Main output control (silent-mode?) C ================================================================== iprintctl(1) = .true. ! Banner iprintctl(2) = .true. ! Parameters and problem processing iprintctl(3) = .true. ! Warnings and errors messages iprintctl(4) = .true. ! Screen-mirror file algencan.out iprintctl(5) = .true. ! Solution file solution.txt iprintctl(6) = .false. ! Statistics file with table line iprintctl(7) = .true. ! User-provided subs calls counters and timing oprint = iprint if (oprint .eq. 0 ) then do i = 1,7 iprintctl(i) = .false. end do end if if ( iprintctl(4) ) then open(unit=10,file=ifile,status='replace') else open(unit=10, status='scratch') end if C ================================================================== C Set solver arguments using the specification file C ================================================================== call fparam(epsfeas,epsopt,efacc,eoacc,oprint,ncomp) C Outer and inner iterations output detail iprintout = oprint / 10 iprintinn = mod( oprint, 10 ) C Error tracker inform = 0 C Stop if mandatory subroutines were not properly provided if ( fcsubt .eq. 0 ) then if ( iprintctl(3) ) then C write(* ,1000) write(10,1000) end if if ( safemode ) then inform = - 88 call reperr(inform) return end if end if C ================================================================== C Initialize problem data structures C ================================================================== call sinip(n,x,l,u,m,lambda,equatn,linear,coded,checkder,inform) if ( inform .lt. 0 ) return nprint = min( n, ncomp ) mprint = min( m, ncomp ) C ================================================================== C Call the solver C ================================================================== C ALGENCAN for PNL problems if ( .not. ignoref .and. m .gt. 0 ) then call auglag(n,x,l,u,m,lambda,equatn,linear,epsfeas,epsopt, + efacc,eoacc,f,c,cnorm,snorm,nl,nlpsupn,fu,cu,cnormu,fub, + cnormub,fb,cnormb,nlpsupnb,ncsupn,rsupn,outiter,totiter, + nwcalls,nwtotit,msqcalls,msqtotit,innfail,alinfo,inform) solinfo = alinfo C GENCAN for box-constrained problems and feasibility problems else maxit = 999999999 C Used in feasibility problems (ignoref=true). With lambda=0 and C rho=1, to minimize 1/2 of the squared infeasibility coincides C with minimizing the augmented Lagrangian. do j = 1,m lambda(j) = 0.0d0 rho(j) = 1.0d0 end do call gencan(n,x,l,u,m,lambda,equatn,linear,rho,epsfeas,epsopt, + maxit,iter,f,nl,nlpsupn,cnorm,cnormu,geninfo,inform) solinfo = geninfo ncsupn = 0.0d0 rsupn = 0.0d0 outiter = 0 totiter = iter nwcalls = 0 nwtotit = 0 msqcalls = 0 msqtotit = 0 innfail = .false. if ( ignoref ) then f = 0.0d0 end if fu = f fb = f fub = fu cnormb = cnorm cnormub = cnormu nlpsupnb = nlpsupn end if if ( inform .lt. 0 ) return C ================================================================== C End problem data structures C ================================================================== call sendp(n,x,l,u,m,lambda,equatn,linear,inform) if ( inform .lt. 0 ) return C ================================================================== C Stop timing C ================================================================== time = dtime(dum) time = dum(1) C ================================================================== C Write statistics C ================================================================== if ( iprintctl(7) ) then C write(* ,9000) time write(10,9000) time toteccnt = 0 totejccnt = 0 totehccnt = 0 do j = 1,m toteccnt = toteccnt + eccnt(j) totejccnt = totejccnt + ejccnt(j) totehccnt = totehccnt + ehccnt(j) end do if ( m .gt. 0 ) then aveeccnt = toteccnt / m aveejccnt = totejccnt / m aveehccnt = totehccnt / m else aveeccnt = 0 aveejccnt = 0 aveehccnt = 0 end if C write(* ,9010) fcoded,efcnt,gcoded,egcnt,hcoded,ehcnt,ccoded, C + toteccnt,aveeccnt,jaccoded,totejccnt,aveejccnt, C + hccoded,totehccnt,aveehccnt,fccoded,efccnt, C + gjaccoded,egjccnt,gjacpcoded,egjcpcnt,hlcoded, C + ehlcnt,hlpcoded,ehlpcnt write(10,9010) fcoded,efcnt,gcoded,egcnt,hcoded,ehcnt,ccoded, + toteccnt,aveeccnt,jaccoded,totejccnt,aveejccnt, + hccoded,totehccnt,aveehccnt,fccoded,efccnt, + gjaccoded,egjccnt,gjacpcoded,egjcpcnt,hlcoded, + ehlcnt,hlpcoded,ehlpcnt end if C ================================================================== C Close output file C ================================================================== close(10) C ================================================================== C Write statistics file with table line C ================================================================== if ( iprintctl(6) ) then open(20,file='algencan-tabline.out') write(20,9040) fu,cnormu,f,cnorm,nlpsupn,fub,cnormub,fb, + cnormb,nlpsupnb,ncsupn,rsupn,inform,solinfo, + innfail,n,m,outiter,totiter,fcnt,nwcalls, + nwtotit,msqcalls,msqtotit,time close(20) end if C ================================================================== C NON-EXECUTABLE STATEMENTS C ================================================================== 1000 format(/,1X,'*** Mandatory subroutines are not being ', + 'provided properly ***',/,1X,'For unconstrained and ', + 'bound-constrained problems, EVALF must be coded by ', + 'the',/,1X,'user. For constrained problems, EVALF ', + 'and EVALC, or, alternatively, EVALFC',/,1X,'must ', + 'be coded. (Note that EVALF/EVALC and EVALFC should ', + 'not be provided',/,1X,'concurrently.) For ', + 'feasibility problems, a constant ', + 'null objective function',/,1X,'must be coded and ', + 'the problem solved with the ', + 'IGNORE-OBJECTIVE-FUNCTION',/,1X,'keyword. Coded ', + 'subroutines must be indicated by setting ', + 'the entrances of array',/,1X,'named CODED ', + 'within subroutine INIP.') 9000 format(/,1X,'Total CPU time in seconds: ',F8.2) 9010 format(/,1X,'User-provided subroutines calls counters: ',/, + /,1X,'Subroutine evalf (coded=',L1,'): ',I6, + /,1X,'Subroutine evalg (coded=',L1,'): ',I6, + /,1X,'Subroutine evalh (coded=',L1,'): ',I6, + /,1X,'Subroutine evalc (coded=',L1,'): ',I6, + 1X,'(',I6,' calls per constraint in average)', + /,1X,'Subroutine evaljac (coded=',L1,'): ',I6, + 1X,'(',I6,' calls per constraint in average)', + /,1X,'Subroutine evalhc (coded=',L1,'): ',I6, + 1X,'(',I6,' calls per constraint in average)', + /,1X,'Subroutine evalfc (coded=',L1,'): ',I6, + /,1X,'Subroutine evalgjac (coded=',L1,'): ',I6, + /,1X,'Subroutine evalgjacp (coded=',L1,'): ',I6, + /,1X,'Subroutine evalhl (coded=',L1,'): ',I6, + /,1X,'Subroutine evalhlp (coded=',L1,'): ',I6,/) 9040 format(1X,1P,D24.16,1X,1P,D7.1,1X,1P,D24.16,1X,1P,D7.1,1X,1P,D7.1, + 1X,1P,D24.16,1X,1P,D7.1,1X,1P,D24.16,1X,1P,D7.1,1X,1P,D7.1, + 1X,1P,D7.1,1X,1P,D7.1,1X,I3,1X,I1,1X,L1,1X,I6,1X,I6,1X,I3, + 1X,I7,1X,I7,1X,I2,1X,I7,1X,I7,1X,I7,0P,F8.2) end