c---------------------------------------------------------------Line 1 c Date assambled: July 31, 1998 c Date modified: February 19, 2001 c c c c --------- ======= SPENBAR Package ======------- c c c c User's information: c To use SPENBAR package you must download in separate directory c the SPENBAR.FOR program, as well as the subroutines IP and PROB c contained in the files from PROB directory of this package. c Compile these subroutines. c Link the object files in a task and then run the created task. c c For exaple, if you want to solve the alkilation problem which c is described in the Fortran file ALKIL.FOR you must take the c following steps: c 1) Download SPENBAR.FOR and ALKIL.FOR in a separate directory. c 2) Compile these subroutines to get: SPENBAR.OBJ and ALKIL.OBJ. c 3) Link these object files into a task, let say SPENBAR.EXE. c 4) Run SPENBAR.EXE. You get the solution of the problem in c SOLUTION.DAT file. Some other information concerning the c optimization process you may find in TNSIMPLE.DAT, TN.OUT c and TN.ST8 files. c c The program is taylored to solve problems with: 1000 variables, c 500 inequalities, 500 equalities, 2000 upper and lower bounds, c 800 nonzeros on inequality constraints and 10000 nonzeros in c equality constraints. c For larger problems you should modify these elements acordingly. c c Good luck ! c Remember: Everything for free comes with no guarantee ! c--------------------------------------------------------------- * Date last modified: December 29, 1997 (by Neculai Andrei) * SPENBAR is a penalty-barrier algorithm. *----------------------------------------------------------- * double precision tau,beta,mu,epsmch logical epsup,lamone,eithor integer*2 gh,gm,gs,gc,ght,gmt,gst,gct, wctime common /nbloc/ nb * * 'tn.out': output file for PENBAR * open(1,file='tn.out') open(3,file='tn.st8') * open(unit=4,file='tnsimple.dat',status='unknown') open(unit=2,file='solution.dat',status='unknown') * tau=1.D-8 beta=.9d0 mu=1.d-1 lamone=.false. eithor=.false. epsup=.false. levout=3 * epsmch=1.d0 1 epsmch=epsmch/2.d0 if(1.d0+epsmch.ne.1.d0) go to 1 epsmch=epsmch*2.d0 * * write(4,27) 27 format(/,5x,'Penalty-Barrier Method for Nonlinear Optimization') write(4,28) 28 format(5x,'Truncated Newton method (simple) TN') write(4,30) 30 format(5x,'Program SPENBAR (Sparse Penalty-Barrier Method)'/) write(4,29) 29 format(5x,'February 19, 2001') * * write(3,30) write(3,20)beta write(3,21)mu write(3,22) write(3,23) write(3,24) write(3,25)epsmch 20 format(5x,'beta = ',e14.7) 21 format(5x,'mu = ',e14.7) 22 format(5x,'lamone = .false.') 23 format(5x,'eithor = .false.') 24 format(5x,'epsup = .false.') 25 format(5x,'epsmch = ',e35.20) write(3,11)tau 11 format(/,5x,'Accuracy = ',e14.7,/) write(3,12) 12 format(5x,'n m me mb o-it in-it f-eval ', 1 'quadr inf Execution time') * do i=1,1 nb=i write(*,15)nb, tau 15 format(10x,'nblocks =',i4,6x,'tau = ',e14.7) * *-------------------------------------------------- call gettim(gh,gm,gs,gc) * call nlp(levout,eithor,lamone,epsup,tau,beta,mu,epsmch, 1 n,m,me,mb,itout,init,nfun,nquad,inf,nb) * ======================================================= * call gettim(ght,gmt,gst,gct) call exetim(gh,gm,gs,gc,ght,gmt,gst,gct) write(3,14)n,m,me,mb,itout,init,nfun,nquad,inf,ght,gmt,gst,gct write(4,18) 18 format(/,5x,'Synthesis of Optimization with SPENBAR:',/) write(4,12) write(4,14)n,m,me,mb,itout,init,nfun,nquad,inf,ght,gmt,gst,gct 14 format(2x,i4,2x,i3,2x,i3,1x,i4,3x, 1 i2,5x,i5,4x,i6,3x,i5,4x,i3,6x,3(i2,':'),i2) write(*,12) write(*,14)n,m,me,mb,itout,init,nfun,nquad,inf,ght,gmt,gst,gct wctime=gct+100*(gst+60*(gmt+60*ght)) write(1,16)wctime write(4,16)wctime write(*,16)wctime 16 format(/,5x,'Execution Time = ',i8) write(4,17) 17 format(4x,75('=')) * end do * close(1) close(3) close(4) * end c############################################################## c SUBROUTINES c =============== c############################################################ * Date last modified February 7, 1996 *------------------------------------- * Solving nonlinear constrained minimization problems by * using a penalty-barrier method. * * Set the dimensions according to your problem: * ndim = upper bound on the number of variables * mdim = upper bound on the number of general inequality constraints * (without simple bounds) * medim = upper bound on the number of equality constraints * mbdim = upper bound on the number of simple bound constraints * (need not be bigger than 2*ndim) * nszcdim = upper bound on the number of structural non-zero elements * in Jacobian of the Inequality constraints. * nszedim = upper bound on the number of structural non-zero elements * in Jacobian of the Equality constraints. *---------------------------------------------------------- * subroutine nlp(lo,eithor,lamone,epsup,tau,bbeta,bmu,epsmch, 1 n,m,me,mb,itout,init,nfun,nquad,inf,nb) * parameter (ndim=1000, mdim=500, medim=500, mbdim=2000, 1 nszcdim=800, nszedim=10000, 1 mgdim=mdim+mbdim, 1 nfpar=1+mdim+mbdim+medim, 1 nmult=mdim+mbdim+medim, 1 nwdim=14*ndim, 1 nworkdim=2*nmult+ndim) * common /three/ beta * *----------------------- Summary on the arrays and their dimensions: * * Dimension Arrays: Total: * ---------------------------------------------------------------- * n x,g,gobj,w1,w2 5*n * m c,lambda,c1 3*m * me e,e1 2*me * mb cb,lamsb,cb1,sb 3*mb + mb * n+1 ipgc,ipge 2*(n+1) * nszc gc,gvv,icgc 2*(nszc) * nsze ge,gev,icge 2*(nsze) * m+mb (=mg) dell,c0,qa,qb,qc,quadr 5*(m+mb) + (m+mb) * m+me+mb+1 (=nfpar) fpar m+me+mb+1 * 14*n (=nwdim) w 14*n * 5*(m+mb+me)+2 gg 5*(m+mb+me)+2 * 2*(m+mb+me)+n work 2*(m+mb+me)+n * ---------------------------------------------------------------- * TOTAL: * real*8: 20*n+16*(m+mb)+10*me+2*nsze+2*nszc+3 * integer: 2*n+mb+nszc+nsze+2 * logical: m+mb *------------------------------------------------------------------- * double precision stepmx,redmu,sneps,gnorm, 1 objf,cmin,bmu,bbeta,accrcy,epsmch,f,eps,eta, 1 aw,bw,an,bn,w0,n0,omega,aleta, 1 xnorm,eleps,beta,objfp,tau,tmp, 1 en,lamn,swtol,mup, 1 x(ndim),g(ndim),gobj(ndim),w1(ndim),w2(ndim), 1 c(mdim),c1(mdim),lambda(mdim), 1 e(medim),e1(medim), 1 cb(mbdim),cb1(mbdim),lamsb(mbdim), 1 w(nwdim), 1 fpar(nfpar), 1 gg(nmult*5+2), 1 work(nworkdim), 1 ge(nszedim),gev(nszedim), 1 gc(nszcdim),gvv(nszcdim), 1 qa(mgdim),qb(mgdim),qc(mgdim), 1 dell(mgdim),c0(mgdim) * logical fcrit,strong,quadr(mgdim),lamup,epsup,savepre,eithor, 1 lamone,muup,muupp * integer ipgc(ndim+1),ipge(ndim+1) integer icgc(nszcdim),icge(nszedim) integer sb(mbdim) * * *----------------------- Initialization * * itout: number of outer iterations * init: number of inner iterations * nfun: number of total function calls * nquad: number of function calls that used * the quadratic extrapolation * nterm=1: stop algorithm * beta: quadratic extrapolation parameter * itout=0 init=0 nfun=0 nquad=0 nterm=0 beta=bbeta * * eleps: is the epsilon for cancellation * in floating point operations * strong=.true.: strong Wolfe conditions are used in line search * (recommended for truncated-Newton, set eta=0.25) * eta: constant in the Wolfe conditions of the line search * stepmx: maximal allowable steplength in line search * redmu: factor by which the barrier-penalty parameter * mu is reduced in each outer iteration * savepre: save preconditioner in truncated-Newton (recommended) * swtol: if mu < swtol, then switch to either-or update strategy * (by Andy Conn, Nick Gould, and Philippe Toint) * accrcy: is the accuracy of the computed function * values (needed for TN) * MXINIT and MSGLVL see subroutine lmqn (truncated-Newton) * eleps=1.d2*epsmch strong=.true. eta=.25d0 stepmx=1.D2 redmu=10.d0 savepre=.true. * if (eithor) then swtol=1.d1*bmu else swtol=1.d-20 endif * accrcy=1.D2*epsmch MXINIT = N/2 MSGLVL = -3 if (lo.gt.2) MSGLVL = 1 * *----------------------- Initialize the problem for PENBAR * * n = number of variables * m = number of general inequalities constraints * (without simple bounds) * mb = number of simple bounds * me = number of equalities constraints * if(lo.gt.0) then write(1,900) write(4,900) 900 format(/5x,'Penalty-Barrier Method for Nonlinear Optimization') write(1,1900) write(4,1900) 1900 format( 5x,'Truncated Newton method (simple) TN') write(1,901) write(4,901) 901 format( 5x,49('-'),/) end if * call ini(n,m,mb,me,sb,x,icgc,ipgc,icge,ipge,nszc,nsze,nb) * if (lo.gt.0) then write(1,902) 902 format(/,5x,'Synopsis of the Problem:',/, 1 5x,'------------------------',/) write(1,903)n 903 format(5x,'Number of variables: n =',i4) write(1,904)m 904 format(5x,'Number of inequality constraints: m =',i4) write(1,905)me 905 format(5x,'Number of equality constraints: me =',i4) write(1,906)mb 906 format(5x,'Number of bounded variables: mb =',i4) write(1,907) tau 907 format(5x,'Accuracy of the solution: tau =',e14.7) write(1,908) beta 908 format(5x,'Quadratic extrapolation parameter: beta =',e14.7) write(1,909)bmu 909 format(5x,'Initial penalty-barrier parameter: mu =',e14.7) write(1,910) epsup 910 format(5x,'Updating eps in unconstrained min.:epsup =',L3) write(1,911) epsmch 911 format(5x,'Machine accuracy: epsmch =',e35.25,/) endif write(4,902) ! write in 'tnsimple' for statistics. write(4,903)n write(4,904)m write(4,905)me write(4,906)mb write(4,907)tau write(4,908)beta write(4,909)bmu write(4,910)epsup write(4,911)epsmch * if ((n.gt.ndim).or.(m.gt.mdim).or.(me.gt.medim) * .or.(mb.gt.mbdim)) then write(1,*) 'ndim, mdim, medim or mbdim too small!!' write(1,*) 'n=',n,', ndim=',ndim write(1,*) 'm=',m,', mdim=',mdim write(1,*) 'me=',me,', medim=',medim write(1,*) 'mb=',mb,', mbdim=',mbdim return endif * if (nszc.gt.nszcdim) write(1,*) 'nszcdim too small !' if (nsze.gt.nszedim) write(1,*) 'nszedim too small !' * call prob(n,m,mb,me,sb,x,objf,gobj,c,gc,cb,e,ge,nszc,nsze,nb) * * *----------------------- Print the initial problem function values * if (lo.gt.1) then write(1,920) 920 format(5x,'Print the initial problem function values:',/, 1 5x,'------------------------------------------',/) write(1,921)objf 921 format(5x,'Initial Objective Function Value =',e14.7) write(1,922) 922 format(/,5x,'Initial x: (Initial values of variables)') write(1,923)(i,x(i),i=1,n) 923 format(3(2x,i4,4x,e14.7)) if(m.ne.0) then write(1,924) 924 format(/,5x,'Initial c: (Initial values of inequalities)') write(1,923)(i,c(i),i=1,m) end if if(me.ne.0) then write(1,925) 925 format(/,5x,'Initial e: (Initial values of equalities)') write(1,923)(i,e(i),i=1,me) end if if(mb.ne.0) then write(1,926) 926 format(/,5x,'Initial cb: (Initial values of var. bounds)') write(1,923)(i,cb(i),i=1,mb) end if endif * *----------------------- Initialization has to be such that the * initial point is strictly feasible with * respect to the simple bounds * cmin=1.d6 do 5 i=1,mb 5 if (cb(i).lt.cmin) cmin=cb(i) if (cmin.le.0.d0) then write(1,*) * 'Initial point NOT strictly feasible w.r.t. simple bounds' write(1,*) 'nonfeasibility =',cmin stop endif * *----------------------- Fix the dimension parameters: * * nfp = dimension of array fpar: fpar = [mu, lambda], * where lambda is the vector of the Lagrangian * multiplier estimates of the inequalities, the simple * bounds, and the equalities -- in this order * nlam = dimension of multiplier vector lambda * the first mgsb simple bounds are handled as general constraints * nw = the dimension of the working vector w * nwork= the dimension of the working vector work * (needed for TN and to pass vectors through subroutines) * nfp=1+m+mb+me nlam=nfp-1 mg=m mgsb=mg-m nw=14*n nwork=2*(m+mb+me)+n * * maxit = maximum number of outer iterations * (when using either-or-update strategy, the algorithm might * fail due to never updating lambda, or the upper bound on * the number of outer iterations might be too small) * mxfun = max. number of function evals * mxline = max. number of function evals in a single line search * mxinit = max. number of inner iterations * maxit=2*n if (maxit.gt.50) maxit=50 if (maxit.le.20) maxit=20 if (eithor) maxit=maxit*2 mxfun=5000*n mxline=50 mxinit=n/2 if (mxinit.gt.50) mxinit=50 if (mxinit.le.0) mxinit=1 * * inf = number of outer iterations with infeasible iterate * cmin = minimum value of all constraints * c0: shifts within the modified log-terms (set such that * ill-conditioning caused by the quadratic extrapolation * is avoided) * maximal c0 = 100; other possible choices for c0: c0=1 (Polyak) * inf=0 cmin=1.D6 do 10 i=1,m if (c(i).lt.cmin) cmin=c(i) c c0(i)=1.D0 c0(i)=dmax1(1.d0,-c(i)) if (c0(i).gt.1.d2) c0(i)=1.d2 10 continue do 20 i=1,mb if (cb(i).lt.cmin) cmin=cb(i) 20 continue * if (lo.gt.1.and.m.ne.0)then write(1,927) 927 format(/,5x,'Initial Shifts: (to avoid ill-conditioning)') write(1,923)(i,c0(i),i=1,mg) end if * *----------------------- Initialize the vector: fpar = [mu, lambda] * If lamone=.false. then calculate lambda * to minimize the norm of the 1st order conditions, * else set the initial lambda to all one. * * In fparini arrays lambda, work and d are * distributed in array work as: * work[ lambda(nlam) | work(nlam) | d(n) ] * call fparini(n,m,mb,me,bmu,gobj,ge,gc,sb,work(1),fpar, 1 c0,mg,qa,qb,qc,lamone,nfp,nlam,work(nlam+1), 1 gg,work(2*nlam+1),eleps,lo, 1 ipgc,ipge,icgc,icge,nszc,nsze) * * if (lo.gt.1) then write(1,928) 928 format(/,5x,'Initial Lagrange multipliers:') write(1,923)(i-1,fpar(i),i=2,nfp) end if * *----------------------- Initialize the parameters used by * Conn, Gould, and Toint * (see FUNDP report 92/16 p.12) * aw=1.d0 bw=1.d0 an=.1d0 bn=.9d0 w0=1.d0 n0=.1258925d0 omega=w0*(fpar(1)**aw) aleta=n0*(fpar(1)**an) * *----------------------- End of initialization. * Start the outer iterations * 100 itout=itout+1 write(*,101) itout 101 format(5x,'itout =',i4) nfunp=nfun objfp=objf if (itout.gt.maxit) then write(1,*) 'STOP - more than',maxit,' outer iterations' return endif * * * eps: convergence tolerance for unconstrained minimization * epsup = .true. : update eps to allow for relatively * loose convergence at the beginning * epsup = .false.: eps=tau, i.e. same criteria as for KKT-conditions * sneps: criteria for norm of search vector (if small, * quit inner iteration) * * if (epsup) then tmp=1.+itout eps=10.**(-tmp) if (eps.lt.tau) eps=tau else eps=tau endif sneps=dmax1(1.d-10,eps**3) * *----------------------- Inner iteration loop: * Call truncated-Newton by Stephen Nash * if(lo.gt.1) write(1,946) itout 946 format(/,1x,'Print out the results of outer iteration:',i3,/, 1 1x,'--------------------------------------------') if(msglvl.ge.1) write(1,945) 945 format(1x,'Truncated Newton minimization iterations:'1x, 1 ' Penalty-Barrier minimization') * CALL LMQN (IERROR, N, X, F, G, W, nw, * MSGLVL, MXINIT, MXFUN, ETA, STEPMX, ACCRCY,strong, * eps,epsmch,iter,gnorm,xnorm,eleps,dell, * objf,gobj,fpar,c,gc,c0,sb,cb,work,nwork, 1 m,mb,me,mg,mgsb,nfp,nfun,nquad,quadr, * savepre,itout,e,ge,sneps,qa,qb,qc, 1 w1,w2,c1,cb1,e1,gvv,gev, 1 ipgc,ipge,icgc,icge,nszc,nsze) * init=init+iter * *----------------------- Printout the result of this outer iteration * if (lo.gt.1) then write(1,930) itout 930 format(/,5x,'Outer Iteration Number: ',i4) write(1,931) iter 931 format(5x,'Number of inner iterations: ',i5) write(1,932) nfun-nfunp 932 format(5x,'Number of function calls in TN-Newton:',i5) write(1,933) fpar(1) 933 format(5x,'Parameter mu = ',e14.7) * if (ierror.ne.0) then write(1,934) ierror 934 format(5x,'Truncated Newton LMQN-ierror =',i4) end if * write(1,935) objf 935 format(5x,'Original Objective function value: ',e14.7) write(1,936) 936 format(/,5x,'Variables values. Array x:') write(1,923) (i,x(i),i=1,n) if(m.ne.0) then write(1,937) 937 format(/,5x,'Value of inequalities constraints. Array c:') write(1,923) (i,c(i),i=1,m) end if if(me.ne.0) then write(1,938) 938 format(/,5x,'Value of equalities constraints. Array e:') write(1,923) (i,e(i),i=1,me) end if * if(mb.ne.0) then !!!!! 2001, February 15. * write(1,939) *939 format(/,5x,'Value of bounds of variables. Array cb:') * write(1,923) (i,cb(i),i=1,mb) * end if endif * *----------------------- Decision on the way. * * Determine whether to update the penalty-barrier parameter * and/or the multipliers estimates: * Update the barrier/penalty parameter mu and the multipliers every * iteration as long as mu is not too small. * If mu is smaller than a threshold (swtol), we update either the * multipliers or mu, depending on the progress we made (following * a procedure given by Conn, Gould, and Toint). * With muupp we save whether mu was reduced from the previous to this * iteration (needed for the fractional-change-in-the-objective- * function stopping-criteria). * * Remark: mu is updated in the following lines, whereas the update of * lambda is done in the subroutine 'stop' (if necessary) * muupp=muup muup=.true. mup=fpar(1) if (fpar(1).ge.swtol) then lamup=.true. fpar(1)=fpar(1)/redmu do 110 i=1,m c0(i)=dmax1(1.d0,-c(i)) if (c0(i).gt.1.d2) c0(i)=1.d2 110 continue * call quad(fpar(1),c0,mg,qa,qb,qc) * omega=w0*(fpar(1)**bw) aleta=n0*(fpar(1)**bn) else if (lo.gt.1) write(1,*) * 'either-or update-strategy' en=0.d0 do 510 i=1,me 510 en=en+e(i)*e(i) en=dsqrt(en) lamn=0.d0 do 530 i=1,m lambda(i)=fpar(i+1)*dell(i) 530 lamn=lamn+(fpar(i+1)-lambda(i))*(fpar(i+1)-lambda(i)) do 540 i=1,mb lamsb(i)=fpar(m+1+i)*fpar(1)/cb(i) 540 lamn=lamn+(fpar(m+i+1)-lamsb(i))*(fpar(m+i+1)-lamsb(i)) lamn=dsqrt(lamn)*fpar(1) lamn=lamn+en if (lamn.le.aleta) then lamup=.true. muup=.false. omega=omega*(fpar(1)**bw) aleta=aleta*(fpar(1)**bn) else lamup=.false. fpar(1)=fpar(1)/redmu do 210 i=1,m c0(i)=dmax1(1.d0,-c(i)) if (c0(i).gt.1.d2) c0(i)=1.d2 210 continue * call quad(fpar(1),c0,mg,qa,qb,qc) * omega=w0*(fpar(1)**aw) aleta=n0*(fpar(1)**an) endif endif * c----------------------- Check stopping criteria whether * to keep on iterating c call stop(n,m,mb,me,mg,mgsb,nfp,itout, * c,cb,sb,objf,gobj,gc,cmin,fpar,nterm,tau,xnorm,gnorm, * objfp,dell,eleps,fcrit,e,ge,lamup,mup,muupp,lo, * ipgc,ipge,icgc,icge,nszc,nsze) * if (cmin.lt.0.d0) inf=inf+1 * if (nterm.eq.0) goto 100 c c----------------------- Final printout c if (lo.gt.0) then write(1,950) 950 format(//,5x,'Solution of the Problem:',/, 1 5x,'========================') write(1,951) itout 951 format(5x,'Total number of Outer Iterations:',i2) write(1,952) init 952 format(5x,'Total number of Inner Iterations:',i5) write(1,953) nfun 953 format(5x,'Total number of function calls :',i5) write(1,954) nquad 954 format(5x,'Total number of quadratic extrapolations :',i5) write(1,955) inf 955 format(5x,'Total number of infeasible outer iters: ',i5) write(1,957) objf 957 format(5x,'Optimal value of objective function:',e14.7) write(1,958) fpar(1) 958 format(5x,'Final value of parameter mu:',e14.7) * write(4,950) write(4,951)itout write(4,952)init write(4,953)nfun write(4,954)nquad write(4,955)inf write(4,957)objf write(4,958)fpar(1) * write(1,956) 956 format(/,5x,'Optimal value of variables:') write(1,923)(i,x(i),i=1,n) * write(4,956) write(4,923)(i,x(i),i=1,n) do i=1,n write(2,1101) x(i) 1101 format(e20.13,',') end do * if(m.ne.0) then write(1,959) write(4,959) 959 format(/,5x,'Final value of inequality constraints:') write(1,923)(i,c(i),i=1,m) write(4,923)(i,c(i),i=1,m) end if * if(me.ne.0) then write(1,960) write(4,960) 960 format(/,5x,'Final value of equality constraints:') write(1,923)(i,e(i),i=1,me) write(4,923)(i,e(i),i=1,me) end if * if(mb.ne.0) then write(1,961) 961 format(/,5x,'Final value of boundary constraints:') write(1,923)(i,cb(i),i=1,mb) end if * write(1,962) write(4,962) 962 format(/,5x,'Final value of Lagrange multipliers:') write(1,923)(i-1,fpar(i),i=2,nfp) write(4,923)(i-1,fpar(i),i=2,nfp) * write(1,963) 963 format(//,5x,'Synthesis of Optimization',/, 1 5x,'-------------------------',/) write(1,*) ' o-it in-it f-eval ', 1 ' quadr #inf objf KKT mu0' * write(1,964) itout,init,nfun,nquad,inf,objf, 1 (.not.fcrit),bmu 964 format(i5,i7,i7,i7,i6,F25.9,L3,F10.4) endif * return end c############################################################## *------------------------------------------------------------- * Date last modified: February 7, 1996 * ----------------------------------- * Calculate the penalty-barrier function and its gradient * at the iterate x (safeguarded w.r.t. elimination) * * f is the penalty-barrier function. * g is the gradient of penalty-barrier function. *------------------------------------------------------------- * subroutine barfct(n,m,mb,me,mg,mgsb,nfp,nfun,nq, 1 x,f,g,objf,gobj,c,gc,cb,sb,c0,fpar, 1 eleps,dell,delq,l,lsb,quadr,e,ge,qa,qb,qc, 1 ipgc,ipge,icgc,icge,nszc,nsze) * common /three/ beta common /nbloc/ nb * double precision x(n),f,g(n),c(m),gc(nszc),objf,gobj(n), 1 fpar(nfp),l(mg),sl,delq(mg),dell(mg),cb(mb),lsb(mb-mgsb), 1 c0(mg),u1,u2,u3,u4,u5,u6,u7,eleps,beta, 1 e(me),ge(nsze),qa(mg),qb(mg),qc(mg),epsmch * integer ipgc(n+1),ipge(n+1) integer icgc(nszc),icge(nsze) integer sb(mb) * logical quadr(mg) * double precision dlog,dabs * *----------------------- Get problem info at point x * call prob(n,m,mb,me,sb,x,objf,gobj,c,gc,cb,e,ge, 1 nszc,nsze,nb) * * * Remark: * The first mgsb components of the simple bounds cb are NOT handled * separately (mgsb=mg-m) - the others (mgsb+1,..,mb) are handeled * separately and are implemented with c0=0 and NO quadratic * extrapolation (fct. calculation is simplified: gc(i) for i=1,mb is * not stored since we know sb...) * *---------------------------- Calculation of the penalty-barrier * function value. It is stored in f. * epsmch=2.220D-15 do 10 i=mgsb+1,mb if(cb(i).lt.epsmch) cb(i)=epsmch lsb(i-mgsb)=dlog(cb(i)/fpar(1)) 10 continue * nquadr=0 do 20 i=1,m if (c(i).ge.-beta*c0(i)*fpar(1)) then quadr(i)=.false. l(i)=dlog(c0(i)+c(i)/fpar(1)) dell(i)=1.D0/(c0(i)+c(i)/fpar(1)) else nquadr=nquadr+1 quadr(i)=.true. l(i)=.5D0*qa(i)*c(i)*c(i)+qb(i)*c(i)+qc(i) if (dabs(l(i)).lt.eleps*dmax1(dabs(qb(i)*c(i)),dabs(qc(i)))) 1 l(i)=0.D0 delq(i)=qa(i)*c(i)+qb(i) if (dabs(delq(i)).lt.eleps*dabs(qb(i))) delq(i)=0.D0 dell(i)=fpar(1)*delq(i) endif 20 continue * do 30 i=m+1,m+mgsb if (cb(i-m).ge.-beta*c0(i)*fpar(1)) then quadr(i)=.false. l(i)=dlog(c0(i)+cb(i-m)/fpar(1)) dell(i)=1.D0/(c0(i)+cb(i-m)/fpar(1)) else nquadr=nquadr+1 quadr(i)=.true. l(i)=.5D0*qa(i)*cb(i-m)*cb(i-m)+qb(i)*cb(i-m)+qc(i) if (dabs(l(i)).lt.eleps*dmax1(dabs(qb(i)*cb(i-m)), * dabs(qc(i)))) l(i)=0.D0 delq(i)=qa(i)*cb(i-m)+qb(i) if (dabs(delq(i)).lt.eleps*dabs(qb(i))) delq(i)=0.D0 dell(i)=fpar(1)*delq(i) endif 30 continue * *---------------------------- Contribution of the inequality * constraints to f. * if (nquadr.gt.0) nq=nq+1 f=objf u2=0.D0 u3=0.D0 do 40 i=1,mg u1=fpar(i+1)*l(i) if (u1.lt.0.D0) then u2=u2+fpar(i+1)*l(i) else u3=u3+fpar(i+1)*l(i) endif 40 continue if (mg.gt.0) then if (dabs(u2+u3).lt.eleps*dabs(u2)) then u2=0.D0 u3=0.D0 endif sl=u2+u3 f=f-fpar(1)*sl if (dabs(objf-fpar(1)*sl).lt.eleps*dabs(objf)) f=0.D0 endif * *---------------------------- Contribution of the simple bounds * constraints to f. * u2=0.D0 u3=0.D0 do 50 i=1,mb-mgsb u1=fpar(mg+i+1)*lsb(i) if (u1.lt.0.D0) then u2=u2+fpar(mg+i+1)*lsb(i) else u3=u3+fpar(mg+i+1)*lsb(i) endif 50 continue if ((mb-mgsb).gt.0) then if (dabs(u2+u3).lt.eleps*dabs(u2)) then u2=0.D0 u3=0.D0 endif sl=u2+u3 f=f-fpar(1)*sl if (dabs(f).lt.eleps*dabs(fpar(1)*sl)) f=0.D0 endif * *---------------------------- Contribution of the equality * constraints to f. * u2=0.D0 u3=0.D0 u4=0.D0 do 55 i=1,me u1=fpar(1+m+mb+i)*e(i) if (u1.lt.0.D0) then u2=u2+fpar(1+m+mb+i)*e(i) else u3=u3+fpar(1+m+mb+i)*e(i) endif u4=u4+e(i)*e(i) 55 continue if (me.gt.0) then if (dabs(u2+u3).lt.eleps*dabs(u2)) then u2=0.D0 u3=0.D0 endif f=f-u2-u3 if (dabs(f).lt.eleps*dabs(u2+u3)) f=0.D0 f=f+u4/fpar(1)/2.d0 if (dabs(f).lt.eleps*dabs(u4/fpar(1)/2.d0)) f=0.D0 endif * *---------------------------- Calculation of the gradient of the * penalty-barrier function. * It is stored in array g. do 90 i=1,n g(i)=gobj(i) u2=0.D0 u3=0.D0 u4=0.D0 u5=0.D0 u6=0.D0 u7=0.D0 * *---------------------------- Cycle 60: * Multiply column i of gc with Lagrange * multipliers of inequality constraints * stored in fpar(1+1),...,fpar(m+1) kp=ipgc(i) if(kp.eq.0) go to 61 ii=i 59 ku=ipgc(ii+1) if(ku.eq.0) then ii=ii+1 go to 59 end if ku=ku-1 do 60 k=kp,ku j=icgc(k) if (quadr(j)) then u1=-fpar(1)*fpar(j+1)*delq(j)*gc(k) else u1=-fpar(1+j)*dell(j)*gc(k) endif if (u1.lt.0.D0) then u2=u2+u1 else u3=u3+u1 endif 60 continue 61 do 70 j=m+1,m+mgsb if (sb(j-m).eq.i) then if (quadr(j)) then u1=-fpar(1)*fpar(j+1)*delq(j) else u1=-fpar(j+1)*dell(j) endif if (u1.lt.0.D0) then u2=u2+u1 else u3=u3+u1 endif else if (sb(j-m).eq.-i) then if (quadr(j)) then u1=fpar(1)*fpar(j+1)*delq(j) else u1=fpar(j+1)*dell(j) endif if (u1.lt.0.D0) then u2=u2+u1 else u3=u3+u1 endif endif 70 continue if (mg.gt.0) then if (dabs(u2+u3).lt.eleps*dabs(u2)) then u2=0.D0 u3=0.D0 endif g(i)=g(i)+u2+u3 if (dabs(gobj(i)+u2+u3).lt.eleps*dabs(gobj(i))) g(i)=0.D0 endif u2=0.D0 u3=0.D0 u4=0.d0 do 75 j=mgsb+1,mb if (sb(j).eq.i) then u1=-fpar(1)*fpar(m+j+1)/cb(j) else if (sb(j).eq.-i) then u1=fpar(1)*fpar(m+j+1)/cb(j) else u1=0.d0 endif if (u1.lt.0.D0) then u2=u2+u1 else u3=u3+u1 endif 75 continue if ((mb-mgsb).gt.0) then g(i)=g(i)+u2+u3 if (dabs(g(i)).lt.eleps*dabs(u2+u3)) g(i)=0.D0 endif * *---------------------------- Cycle 80: * Multiply column i of ge with Lagrange * multipliers of equality constraints * stored in fpar(1+m+mb+1),..., * fpar(1+m+mb+me). * u2=0.d0 u3=0.d0 kp=ipge(i) if(kp.eq.0) go to 81 ii=i 76 ku=ipge(ii+1) if(ku.eq.0) then ii=ii+1 go to 76 end if ku=ku-1 do 80 k=kp,ku j=icge(k) u1=(-fpar(1+m+mb+j)+e(j)/fpar(1))*ge(k) if (u1.lt.0.d0) then u2=u2+u1 else u3=u3+u1 endif 80 continue 81 if (me.gt.0) then if (dabs(u2+u3).lt.eleps*dabs(u2)) then u2=0.D0 u3=0.D0 endif g(i)=g(i)+u2+u3 if (dabs(g(i)).lt.eleps*dabs(u2+u3)) g(i)=0.D0 endif 90 continue nfun=nfun+1 return end * c############################################################ *---------------------------------------------------------- * Date last modified: February 7, 1996 *---------------------------------------------------------- * nvar Number of variables. * m Number of inequalities constraints. * me Number of equalities constraints. * mb Number of bounded variables. * lambda Current input vector. * ge Gradient of the equality constraints. * gc Gradient of the inequality constraints. * sb Index vector for simpe bounds on vareiables. * gobj Gradient of the objective function. * (of the original problem) * w Working area of dimension n. * * f Value of the function. (output) * g Gradient of the function. (output) * * Please see the Technical Report RRR 12-94, April 1994 of * Marc G. Breitfeld and David F. Shanno, page 5. *---------------------------------------------------------- * subroutine calfgpen(nvar,m,me,mb,lambda,f,g, 1 ge,gc,sb,gobj,w,ipgc,ipge,icgc,icge,nszc,nsze) * real*8 lambda(me+m+mb),g(me+m+mb),w(nvar),f real*8 ge(nsze),gc(nszc),gobj(nvar) * integer ipgc(nvar+1),ipge(nvar+1) integer icgc(nszc),icge(nsze) integer sb(mb) * *---------------------------- Calculation of the function value. * do i=1,nvar w(i)=gobj(i) end do * do 20 j=1,me * *---------------------------- Cycle 10: * Multiply row j of ge with Lagrange * multipliers of equality constraints. * do 10 i=1,nvar kp=ipge(i) if(kp.eq.0) go to 10 ii=i 9 ku=ipge(ii+1) if(ku.eq.0) then ii=ii+1 go to 9 end if ku=ku-1 do k=kp,ku if(icge(k).eq.j) then w(i)=w(i)-lambda(j)*ge(k) go to 10 end if end do 10 continue 20 continue * * do 40 j=1,m * *---------------------------- Cycle 30: * Multiply row j of gc with Lagrange * multipliers of inequality constraints. * do 30 i=1,nvar kp=ipgc(i) if(kp.eq.0) go to 30 ii=i 29 ku=ipgc(ii+1) if(ku.eq.0) then ii=ii+1 go to 29 end if ku=ku-1 do k=kp,ku if(icgc(k).eq.j) then w(i)=w(i)-lambda(j+me)*gc(k) go to 30 end if end do 30 continue 40 continue * do j=1,mb if(sb(j).gt.0) then w(sb(j))= w(sb(j)) - lambda(me+m+j) else w(-sb(j))=w(-sb(j)) + lambda(me+m+j) end if end do * f=0.d0 do i=1,nvar f=f+w(i)*w(i) end do * *---------------------------- Calculation of the gradient. * do j=1,me+m+mb g(j)=0.d0 end do * do 60 j=1,me do 50 i=1,nvar kp=ipge(i) if(kp.eq.0) go to 50 ii=i 49 ku=ipge(ii+1) if(ku.eq.0) then ii=ii+1 go to 49 end if ku=ku-1 do k=kp,ku if(icge(k).eq.j) then g(j)=g(j)-2.0*w(i)*ge(k) go to 50 end if end do 50 continue 60 continue * do 80 j=1,m do 70 i=1,nvar kp=ipgc(i) if(kp.eq.0) go to 70 ii=i 69 ku=ipgc(ii+1) if(ku.eq.0) then ii=ii+1 go to 69 end if ku=ku-1 do k=kp,ku if(icgc(k).eq.j) then g(j+me)=g(j+me)-2.0*w(i)*gc(k) go to 70 end if end do 70 continue 80 continue * do j=1,mb if(sb(j).gt.0) then g(me+m+j)=-2.0*w(sb(j)) else g(me+m+j)= 2.0*w(-sb(j)) end if end do * return end * c############################################################ * Date created: March 25, 1995 * Date last modified: February 7, 1996 *------------------------------------------------------------- * Conjugate Gradient - Shanno - Phua Algorithm. * This subroutine minimizes an unconstrained nonlinear function * of n variables using a variant of Beale restart conjugate * gradient algorithm. (modified by N. Andrei) * * Parameters: n,x,f,g,nfun,iter,eps,nexit,mxfun,work,mdim,acc * are described in subroutine cgsp.for, in: * Technical Report: N. Andrei, Computational Experience with * Conjugate Gradient Algorithms for Large-Scale Unconstrained * Optimization. * ICI, July 21, 1995. * * Parameters: nvar,me,m,mb,gobj,ge,gc,sb,d are used for * computation the value of the minimizing function, in this * case the first order optimality conditions for Lagrange * multipliers computation. *------------------------------------------------------------ * subroutine cgsppen(n,x,f,g,nfun,iter,eps,nexit,mxfun, 1 work,mdim,acc, 1 nvar,me,m,mb,gobj,ge,gc,sb,d, 1 ipgc,ipge,icgc,icge,nszc,nsze) * *------------------------------------------------------------ * double precision x(n),g(n),work(mdim) double precision f,fp,fmin,alpha,at,ap,gsq,dg,dg1 double precision dp,step,acc,dal,u1,u2,u3,u4,eps double precision xsq,rtst,dsqrt,dmin1,dmax1,dabs * real*8 gobj(nvar),ge(nsze),gc(nszc),d(nvar) * integer ipgc(nvar+1),ipge(nvar+1) integer icgc(nszc),icge(nsze) integer sb(mb) * logical rsw * * Initialization of iter,nfun,nexit,ioutk which counts output * iterations. * iter=0 nfun=0 nexit=0 alpha=1.0 dg=1.0 * nx=n ng=nx+n * nry=ng+n nrd=nry+n ncons=5*n ncons1=ncons+1 ncons2=ncons+2 * *------------------------------------------------------------ 20 call calfgpen(nvar,m,me,mb,x,f,g,ge,gc,sb,gobj,d, 1 ipgc,ipge,icgc,icge,nszc,nsze) * nfun=nfun+1 nrst=n rsw=.true. * dg1=0.0 xsq=0.0 do 30 i=1,n work(i)=-g(i) xsq=xsq+x(i)*x(i) 30 dg1=dg1-g(i)*g(i) gsq=-dg1 * if(gsq.le.eps*eps*dmax1(1.0d0,xsq)) go to 600 * * Begin the major iteration loop. * =============================== * 40 fmin=f ncalls=nfun * alpha=alpha*dg/dg1 * if(nrst.eq.1) alpha=1.0 * if(rsw) alpha=1.0/dsqrt(gsq) * ap=0.0 fp=fmin dp=dg1 * dg=dg1 * iter=iter+1 * * Calculate the current steplength and store the current x and g. * step=0.0 do i=1,n step=step+work(i)*work(i) nxpi=nx+i ngpi=ng+i work(nxpi)=x(i) work(ngpi)=g(i) end do step=dsqrt(step) * 80 if(alpha*step.gt.acc) go to 90 * if(.not.rsw) go to 20 nexit=2 go to 600 * * Calculate the trial point. * 90 do 100 i=1,n nxpi=nx+i 100 x(i)=work(nxpi)+alpha*work(i) * call calfgpen(nvar,m,me,mb,x,f,g,ge,gc,sb,gobj,d, 1 ipgc,ipge,icgc,icge,nszc,nsze) * nfun=nfun+1 if(nfun.le.mxfun) go to 110 nexit=1 go to 600 * 110 dal=0.0 do i=1,n dal=dal+g(i)*work(i) end do * if(f.gt.fmin.and.dal.lt.0.0) go to 160 * if(f.gt.(fmin+0.0001*alpha*dg).or.dabs(dal/dg) 1 .gt.(0.9)) go to 130 * if((nfun-ncalls).le.1.and.dabs(dal/dg).gt.eps) 1 go to 130 go to 170 * * 130 u1=dp+dal-3.0*(fp-f)/(ap-alpha) u2=u1*u1-dp*dal if(u2.lt.0.0) u2=0.0 u2=dsqrt(u2) at=alpha-(alpha-ap)*(dal+u2-u1)/(dal-dp+2.0*u2) * * Test whether the line minimum has been bracketed. * if((dal/dp).gt.0.0) go to 140 * if(at.lt.(1.01*dmin1(alpha,ap)).or.at.gt.(0.99*dmax1 1 (alpha,ap))) at=(alpha+ap)/2.0 * go to 150 * 140 if(dal.gt.0.0.and.0.0.lt.at.and.at.lt. 1 (0.99*dmin1(ap,alpha))) go to 150 * if(dal.le.0.0.and.at.gt.(1.01*dmax1(ap,alpha))) go to 150 * if(dal.le.0.0) at=2.0*dmax1(ap,alpha) * if(dal.gt.0.0) at=dmin1(ap,alpha)/2.0 * * Set ap=alpha, alpha=at, and continue search. * 150 ap=alpha fp=f dp=dal alpha=at go to 80 * 160 alpha=alpha/3.0 ap=0.0 fp=fmin dp=dg go to 80 * * The line search has converged. Stop linear search iterations. * 170 continue * *------------------------------------------------------------ * Test for convergence of the algorithm. * gsq=0.0 xsq=0.0 do 180 i=1,n gsq=gsq+g(i)*g(i) 180 xsq=xsq+x(i)*x(i) * if(gsq.le.eps*eps*dmax1(1.0d0,xsq)) go to 600 *------------------------------------------------------------ * * Search continues. Set work(i)=alpha*work(i), the full step vector. * do 190 i=1,n 190 work(i)=alpha*work(i) * rtst=0.0 do 200 i=1,n ngpi=ng+i 200 rtst=rtst+g(i)*work(ngpi) if(dabs(rtst/gsq).gt.0.2) nrst=n * if(nrst.ne.n) go to 220 work(ncons+1)=0.0 work(ncons+2)=0.0 do 210 i=1,n nrdpi=nrd+i nrypi=nry+i ngpi=ng+i work(nrypi)=g(i)-work(ngpi) work(nrdpi)=work(i) work(ncons1)=work(ncons1)+work(nrypi)*work(nrypi) 210 work(ncons2)=work(ncons2)+work(i)*work(nrypi) * * Calculate the restart Hessian times the current gradient. * 220 u1=0.0 u2=0.0 do 230 i=1,n nrdpi=nrd+i nrypi=nry+i u1=u1-work(nrdpi)*g(i)/work(ncons1) 230 u2=u2+work(nrdpi)*g(i)*2.0/work(ncons2)- 1 work(nrypi)*g(i)/work(ncons1) u3=work(ncons2)/work(ncons1) do 240 i=1,n nxpi=nx+i nrdpi=nrd+i nrypi=nry+i 240 work(nxpi)=-u3*g(i)-u1*work(nrypi)-u2*work(nrdpi) * if(nrst.eq.n) go to 300 * 250 u1=0.0 u2=0.0 u3=0.0 u4=0.0 do 260 i=1,n ngpi=ng+i nrdpi=nrd+i nrypi=nry+i u1=u1-(g(i)-work(ngpi))*work(nrdpi)/work(ncons1) u2=u2-(g(i)-work(ngpi))*work(nrypi)/work(ncons1) 1 +2.0*work(nrdpi)*(g(i)-work(ngpi))/work(ncons2) 260 u3=u3+work(i)*(g(i)-work(ngpi)) step=0.0 do 270 i=1,n ngpi=ng+i nrdpi=nrd+i nrypi=nry+i step=(work(ncons2)/work(ncons1))*(g(i)-work(ngpi))+ 1 u1*work(nrypi)+u2*work(nrdpi) u4=u4+step*(g(i)-work(ngpi)) 270 work(ngpi)=step * * Calculate the doubly updated Hessian times the current * gradient to obtain the search vector. * u1=0.0 u2=0.0 do 280 i=1,n u1=u1-work(i)*g(i)/u3 ngpi=ng+i 280 u2=u2+(1.0+u4/u3)*work(i)*g(i)/u3-work(ngpi)*g(i)/u3 do 290 i=1,n ngpi=ng+i nxpi=nx+i 290 work(nxpi)=work(nxpi)-u1*work(ngpi)-u2*work(i) * * Calculate the derivative along the new search vector. * 300 dg1=0.0 do 310 i=1,n nxpi=nx+i work(i)=work(nxpi) 310 dg1=dg1+work(i)*g(i) * * If the new direction is not a descent direction, STOP. * if(dg1.gt.0.0) go to 320 * if(nrst.eq.n) nrst=0 nrst=nrst+1 rsw=.false. go to 40 * * Roundoff has produced a bad direction. * 320 nexit=3 * 600 continue * return end c############################################################ c-------------------------------------------------------------- SUBROUTINE CHKUCP(LWTEST,MAXFUN,NWHY,N,ALPHA,EPSMCH, * ETA,PEPS,RTEPS,RTOL,RTOLSQ,STEPMX,TEST, * XTOL,XNORM,X,LW,SMALL,TINY,ACCRCY) IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER LW,LWTEST,MAXFUN,NWHY,N DOUBLE PRECISION ACCRCY,ALPHA,EPSMCH,ETA,PEPS,RTEPS,RTOL, * RTOLSQ,STEPMX,TEST,XTOL,XNORM,SMALL,TINY DOUBLE PRECISION X(N) C C CHECKS PARAMETERS AND SETS CONSTANTS WHICH ARE COMMON TO BOTH C DERIVATIVE AND NON-DERIVATIVE ALGORITHMS C DOUBLE PRECISION DABS,DSQRT,dnrm2 SMALL = EPSMCH*EPSMCH TINY = SMALL NWHY = -1 RTEPS = DSQRT(EPSMCH) RTOL = XTOL IF (DABS(RTOL) .LT. ACCRCY) RTOL = 1.D1*RTEPS C C CHECK FOR ERRORS IN THE INPUT PARAMETERS C IF (LW .LT. LWTEST * .OR. N .LT. 1 .OR. RTOL .LT. 0.D0 .OR. ETA .GE. 1.D0 .OR. * ETA .LT. 0.D0 .OR. STEPMX .LT. RTOL .OR. * MAXFUN .LT. 1) RETURN NWHY = 0 C C SET CONSTANTS FOR LATER C RTOLSQ = RTOL*RTOL PEPS = ACCRCY**0.6666D0 XNORM = DNRM2(N,X,1) ALPHA = 0.D0 TEST = 0.D0 RETURN END C c############################################################ C%% TRUNCATED-NEWTON METHOD: BLAS C NOTE: ALL ROUTINES HERE ARE FROM LINPACK WITH THE EXCEPTION C OF DXPY (A VERSION OF DAXPY WITH A=1.0) C WRITTEN BY: STEPHEN G. NASH C OPERATIONS RESEARCH AND APPLIED STATISTICS DEPT. C GEORGE MASON UNIVERSITY C FAIRFAX, VA 22030 C**************************************************************** SUBROUTINE DAXPY(N,DA,DX,INCX,DY,INCY) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C C CONSTANT TIMES A VECTOR PLUS A VECTOR. C USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE. C JACK DONGARRA, LINPACK, 3/11/78. C DOUBLE PRECISION DX(1),DY(1),DA INTEGER I,INCX,INCY,IX,IY,M,MP1,N C IF(N.LE.0)RETURN IF (DA .EQ. 0.0D0) RETURN IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20 C C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS C NOT EQUAL TO 1 C IX = 1 IY = 1 IF(INCX.LT.0)IX = (-N+1)*INCX + 1 IF(INCY.LT.0)IY = (-N+1)*INCY + 1 DO 10 I = 1,N DY(IY) = DY(IY) + DA*DX(IX) IX = IX + INCX IY = IY + INCY 10 CONTINUE RETURN C C CODE FOR BOTH INCREMENTS EQUAL TO 1 C C C CLEAN-UP LOOP C 20 M = MOD(N,4) IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M DY(I) = DY(I) + DA*DX(I) 30 CONTINUE IF( N .LT. 4 ) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,4 DY(I) = DY(I) + DA*DX(I) DY(I + 1) = DY(I + 1) + DA*DX(I + 1) DY(I + 2) = DY(I + 2) + DA*DX(I + 2) DY(I + 3) = DY(I + 3) + DA*DX(I + 3) 50 CONTINUE RETURN END c c############################################################ SUBROUTINE DCOPY(N,DX,INCX,DY,INCY) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C C COPIES A VECTOR, X, TO A VECTOR, Y. C USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE. C JACK DONGARRA, LINPACK, 3/11/78. C DOUBLE PRECISION DX(1),DY(1) INTEGER I,INCX,INCY,IX,IY,M,MP1,N C IF(N.LE.0)RETURN IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20 C C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS C NOT EQUAL TO 1 C IX = 1 IY = 1 IF(INCX.LT.0)IX = (-N+1)*INCX + 1 IF(INCY.LT.0)IY = (-N+1)*INCY + 1 DO 10 I = 1,N DY(IY) = DX(IX) IX = IX + INCX IY = IY + INCY 10 CONTINUE RETURN C C CODE FOR BOTH INCREMENTS EQUAL TO 1 C C C CLEAN-UP LOOP C 20 M = MOD(N,7) IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M DY(I) = DX(I) 30 CONTINUE IF( N .LT. 7 ) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,7 DY(I) = DX(I) DY(I + 1) = DX(I + 1) DY(I + 2) = DX(I + 2) DY(I + 3) = DX(I + 3) DY(I + 4) = DX(I + 4) DY(I + 5) = DX(I + 5) DY(I + 6) = DX(I + 6) 50 CONTINUE RETURN END c c############################################################ DOUBLE PRECISION FUNCTION ddots(N,DX,INCX,DY,INCY,eleps) IMPLICIT DOUBLE PRECISION (A-H,O-Z) * * FORMS THE DOT PRODUCT OF TWO VECTORS. * Sum up positive and negative part separately and check * for cancellation. * (tolerance = eleps) * DOUBLE PRECISION DX(1),DY(1),u1,u2,eleps INTEGER I,INCX,INCY,n * if (incx.gt.1 .or. incy.gt.1) then write(1,*) 'ddot not set up for incx or incy > 1' stop endif * ddots = 0.D0 * IF(N.LE.0)RETURN * u1 = 0.0D0 u2 = 0.0D0 DO 10 I = 1,n if ((dx(i).ge.0.d0 .and. dy(i).ge.0.d0) .or. 1 (dx(i).le.0.d0 .and. dy(i).le.0.d0)) then u1=u1+DX(I)*DY(I) else u2=u2+DX(I)*DY(I) endif 10 CONTINUE * ddots=u1+u2 * if (dabs(ddots).lt.eleps*dabs(u1)) ddots=0.d0 * RETURN END * c############################################################ DOUBLE PRECISION FUNCTION DNRM2 ( N, DX, INCX) IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER NEXT DOUBLE PRECISION DX(1),CUTLO,CUTHI,HITEST,SUM,XMAX,ZERO,ONE DATA ZERO, ONE /0.0D0, 1.0D0/ C C EUCLIDEAN NORM OF THE N-VECTOR STORED IN DX() WITH STORAGE C INCREMENT INCX . C IF N .LE. 0 RETURN WITH RESULT = 0. C IF N .GE. 1 THEN INCX MUST BE .GE. 1 C C C.L.LAWSON, 1978 JAN 08 C C FOUR PHASE METHOD USING TWO BUILT-IN CONSTANTS THAT ARE C HOPEFULLY APPLICABLE TO ALL MACHINES. C CUTLO = MAXIMUM OF DSQRT(U/EPS) OVER ALL KNOWN MACHINES. C CUTHI = MINIMUM OF DSQRT(V) OVER ALL KNOWN MACHINES. C WHERE C EPS = SMALLEST NO. SUCH THAT EPS + 1. .GT. 1. C U = SMALLEST POSITIVE NO. (UNDERFLOW LIMIT) C V = LARGEST NO. (OVERFLOW LIMIT) C C BRIEF OUTLINE OF ALGORITHM.. C C PHASE 1 SCANS ZERO COMPONENTS. C MOVE TO PHASE 2 WHEN A COMPONENT IS NONZERO AND .LE. CUTLO C MOVE TO PHASE 3 WHEN A COMPONENT IS .GT. CUTLO C MOVE TO PHASE 4 WHEN A COMPONENT IS .GE. CUTHI/M C WHERE M = N FOR X() REAL AND M = 2*N FOR COMPLEX. C C VALUES FOR CUTLO AND CUTHI.. C FROM THE ENVIRONMENTAL PARAMETERS LISTED IN THE IMSL CONVERTER C DOCUMENT THE LIMITING VALUES ARE AS FOLLOWS.. C CUTLO, S.P. U/EPS = 2**(-102) FOR HONEYWELL. CLOSE SECONDS ARE C UNIVAC AND DEC AT 2**(-103) C THUS CUTLO = 2**(-51) = 4.44089E-16 C CUTHI, S.P. V = 2**127 FOR UNIVAC, HONEYWELL, AND DEC. C THUS CUTHI = 2**(63.5) = 1.30438E19 C CUTLO, D.P. U/EPS = 2**(-67) FOR HONEYWELL AND DEC. C THUS CUTLO = 2**(-33.5) = 8.23181D-11 C CUTHI, D.P. SAME AS S.P. CUTHI = 1.30438D19 C DATA CUTLO, CUTHI / 8.232D-11, 1.304D19 / C DATA CUTLO, CUTHI / 4.441E-16, 1.304E19 / DATA CUTLO, CUTHI / 8.232D-11, 1.304D19 / C IF(N .GT. 0) GO TO 10 DNRM2 = ZERO GO TO 300 C 10 ASSIGN 30 TO NEXT SUM = ZERO NN = N * INCX C BEGIN MAIN LOOP I = 1 20 GO TO NEXT,(30, 50, 70, 110) 30 IF( DABS(DX(I)) .GT. CUTLO) GO TO 85 ASSIGN 50 TO NEXT XMAX = ZERO C C PHASE 1. SUM IS ZERO C 50 IF( DX(I) .EQ. ZERO) GO TO 200 IF( DABS(DX(I)) .GT. CUTLO) GO TO 85 C C PREPARE FOR PHASE 2. ASSIGN 70 TO NEXT GO TO 105 C C PREPARE FOR PHASE 4. C 100 I = J ASSIGN 110 TO NEXT SUM = (SUM / DX(I)) / DX(I) 105 XMAX = DABS(DX(I)) GO TO 115 C C PHASE 2. SUM IS SMALL. C SCALE TO AVOID DESTRUCTIVE UNDERFLOW. C 70 IF( DABS(DX(I)) .GT. CUTLO ) GO TO 75 C C COMMON CODE FOR PHASES 2 AND 4. C IN PHASE 4 SUM IS LARGE. SCALE TO AVOID OVERFLOW. C 110 IF( DABS(DX(I)) .LE. XMAX ) GO TO 115 SUM = ONE + SUM * (XMAX / DX(I))**2 XMAX = DABS(DX(I)) GO TO 200 C 115 SUM = SUM + (DX(I)/XMAX)**2 GO TO 200 C C C PREPARE FOR PHASE 3. C 75 SUM = (SUM * XMAX) * XMAX C C C FOR REAL OR D.P. SET HITEST = CUTHI/N C FOR COMPLEX SET HITEST = CUTHI/(2*N) C 85 HITEST = CUTHI/FLOAT( N ) C C PHASE 3. SUM IS MID-RANGE. NO SCALING. C DO 95 J =I,NN,INCX IF(DABS(DX(J)) .GE. HITEST) GO TO 100 95 SUM = SUM + DX(J)**2 DNRM2 = DSQRT( SUM ) GO TO 300 C 200 CONTINUE I = I + INCX IF ( I .LE. NN ) GO TO 20 C C END OF MAIN LOOP. C C COMPUTE SQUARE ROOT AND ADJUST FOR SCALING. C DNRM2 = XMAX * DSQRT(SUM) 300 CONTINUE RETURN END C c############################################################ C****************************************************************** C SPECIAL BLAS FOR Y = X+Y C****************************************************************** SUBROUTINE DXPY(N,DX,INCX,DY,INCY) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C C VECTOR PLUS A VECTOR. C USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE. C STEPHEN G. NASH 5/30/89. C DOUBLE PRECISION DX(1),DY(1) INTEGER I,INCX,INCY,IX,IY,M,MP1,N C IF(N.LE.0)RETURN IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20 C C CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS C NOT EQUAL TO 1 C IX = 1 IY = 1 IF(INCX.LT.0)IX = (-N+1)*INCX + 1 IF(INCY.LT.0)IY = (-N+1)*INCY + 1 DO 10 I = 1,N DY(IY) = DY(IY) + DX(IX) IX = IX + INCX IY = IY + INCY 10 CONTINUE RETURN C C CODE FOR BOTH INCREMENTS EQUAL TO 1 C C C CLEAN-UP LOOP C 20 M = MOD(N,4) IF( M .EQ. 0 ) GO TO 40 DO 30 I = 1,M DY(I) = DY(I) + DX(I) 30 CONTINUE IF( N .LT. 4 ) RETURN 40 MP1 = M + 1 DO 50 I = MP1,N,4 DY(I) = DY(I) + DX(I) DY(I + 1) = DY(I + 1) + DX(I + 1) DY(I + 2) = DY(I + 2) + DX(I + 2) DY(I + 3) = DY(I + 3) + DX(I + 3) 50 CONTINUE RETURN END c c############################################################ *---------------------------------------------------------------- * Date created : May 30, 1995 * Date last modified : May 30, 1995 * * Subroutine for execution time computation. * *---------------------------------------------------------------- * subroutine exetim(tih,tim,tis,tic,tfh,tfm,tfs,tfc) * integer*2 tih,tim,tis,tic integer*2 tfh,tfm,tfs,tfc * integer*4 ti,tf integer*4 ch,cm,cs data ch,cm,cs/360000,6000,100/ * ti=tih*ch+tim*cm+tis*cs+tic tf=tfh*ch+tfm*cm+tfs*cs+tfc tf=tf-ti tfh=tf/ch tf=tf-tfh*ch tfm=tf/cm tf=tf-tfm*cm tfs=tf/cs tfc=tf-tfs*cs * return end *---------------------------------------------------------------- c############################################################ *------------------------------------------------------------- * Date last modified: February 7, 1996 * ----------------------------------- * Subroutine for initialization of vector fpar. * fpar(1)=mu * fpar(i), for i=1,me+m+mb+1, are the Lagrangian multipliers * corresponding to the m general constraints, * mb simple bounds, and me equalities in this order. * *------------------------------------------------------------- subroutine fparini(n,m,mb,me,bmu,gobj,ge,gc,sb, * lambda,fpar, * c0,mg,qa,qb,qc,lamone,nfp,ngg, * work,gg,d,eleps,lo, 1 ipgc,ipge,icgc,icge,nszc,nsze) * double precision fpar(nfp),bmu,gobj(n),valf, * gc(nszc),ge(nsze),lambda(ngg),c0(mg), * eleps,qa(mg),qb(mg),qc(mg),d(n), * work(ngg),gg(5*ngg+2) * real*8 acci,epsi * integer ipgc(n+1),ipge(n+1) integer icgc(nszc),icge(nsze) integer sb(mb) logical lamone * fpar(1)=bmu * call quad(fpar(1),c0,mg,qa,qb,qc) * *----------------------- If lamone=.true. then set the * initial lambda to all one. * if (lamone) then do 10 i=2,nfp 10 fpar(i)=1.d0 return endif * *----------------------- If lamone=.false. then calculate * the initial lambda by minimizing the * first order optimality conditions. * do i=1,ngg lambda(i)=0.d0 end do epsi=0.00000001 mxfun=1000 mdim=5*ngg+2 acci=eleps*100.0 * call cgsppen(ngg,lambda,valf,work,nfun,iteri,epsi,nexit, 1 mxfun,gg,mdim,acci,n,me,m,mb,gobj,ge,gc,sb,d, 1 ipgc,ipge,icgc,icge,nszc,nsze) * if(lo.gt.1) then write(1,299) 299 format(/,5x,'Initial lagrange multipliers computation:') write(1,300)valf 300 format(5x,'First order cond. value =',e25.15) write(1,310)nexit,iteri,nfun 310 format(5x,'nexit =',i2,' iter =',i6,' nfun =',i6) end if write(*,310)nexit,iteri,nfun * *----------------------- Update the Lagrange multipliers * do 190 i=1,me if (lambda(i).lt.0.d0 .and. lambda(i).gt.-1.d0) lambda(i)=-1.d0 if (lambda(i).ge.0.d0 .and. lambda(i).lt.1.d0) lambda(i)=1.d0 if (lambda(i).lt.-1.D2) lambda(i)=-1.D2 190 if (lambda(i).gt.1.D2) lambda(i)=1.D2 * do 195 i=me+1,me+m+mb if (lambda(i).lt.1.D0) lambda(i)=1.D0 195 if (lambda(i).gt.1.D2) lambda(i)=1.D2 * do 200 i=1,me 200 fpar(1+m+mb+i)=lambda(i) * do 210 i=1,m+mb 210 fpar(1+i)=lambda(me+i) * return end * c############################################################ SUBROUTINE GETPTC(BIG,SMALL,RTSMLL,RELTOL,ABSTOL,TNYTOL, * FPRESN,ETA,RMU,XBND,U,FU,GU,XMIN,FMIN,GMIN, * XW,FW,GW,A,B,OLDF,B1,SCXBND,E,STEP,FACTOR, * BRAKTD,GTEST1,GTEST2,TOL,IENTRY,ITEST,strong) IMPLICIT DOUBLE PRECISION (A-H,O-Z) LOGICAL BRAKTD,strong INTEGER IENTRY,ITEST DOUBLE PRECISION BIG,SMALL,RTSMLL,RELTOL,ABSTOL,TNYTOL, * FPRESN,ETA,RMU,XBND,U,FU,GU,XMIN,FMIN,GMIN, * XW,FW,GW,A,B,OLDF,B1,SCXBND,E,STEP,FACTOR, * GTEST1,GTEST2,TOL,DENOM C C ************************************************************ C GETPTC, AN ALGORITHM FOR FINDING A STEPLENGTH, CALLED REPEATEDLY BY C ROUTINES WHICH REQUIRE A STEP LENGTH TO BE COMPUTED USING CUBIC C INTERPOLATION. THE PARAMETERS CONTAIN INFORMATION ABOUT THE INTERVAL C IN WHICH A LOWER POINT IS TO BE FOUND AND FROM THIS GETPTC COMPUTES A C POINT AT WHICH THE FUNCTION CAN BE EVALUATED BY THE CALLING PROGRAM. C THE VALUE OF THE INTEGER PARAMETERS IENTRY DETERMINES THE PATH TAKEN C THROUGH THE CODE. C ************************************************************ C LOGICAL CONVRG DOUBLE PRECISION ABGMIN,ABGW,ABSR,A1,CHORDM,CHORDU, * D1,D2,P,Q,R,S,SCALE,SUMSQ,TWOTOL,XMIDPT DOUBLE PRECISION ZERO, POINT1,HALF,ONE,THREE,FIVE,ELEVEN C C THE FOLLOWING STANDARD FUNCTIONS AND SYSTEM FUNCTIONS ARE CALLED C WITHIN GETPTC C DOUBLE PRECISION DABS, DSQRT C ZERO = 0.D0 POINT1 = 1.D-1 HALF = 5.D-1 ONE = 1.D0 THREE = 3.D0 FIVE = 5.D0 ELEVEN = 11.D0 C C BRANCH TO APPROPRIATE SECTION OF CODE DEPENDING ON THE C VALUE OF IENTRY. C GOTO (10,20), IENTRY C C IENTRY=1 C CHECK INPUT PARAMETERS C 10 ITEST = 2 IF (U .LE. ZERO .OR. XBND .LE. TNYTOL .OR. GU .GT. ZERO) * RETURN ITEST = 1 IF (XBND .LT. ABSTOL) ABSTOL = XBND TOL = ABSTOL TWOTOL = TOL + TOL C C A AND B DEFINE THE INTERVAL OF UNCERTAINTY, X AND XW ARE POINTS C WITH LOWEST AND SECOND LOWEST FUNCTION VALUES SO FAR OBTAINED. C INITIALIZE A,SMIN,XW AT ORIGIN AND CORRESPONDING VALUES OF C FUNCTION AND PROJECTION OF THE GRADIENT ALONG DIRECTION OF SEARCH C AT VALUES FOR LATEST ESTIMATE AT MINIMUM. C A = ZERO XW = ZERO XMIN = ZERO OLDF = FU FMIN = FU FW = FU GW = GU GMIN = GU STEP = U FACTOR = FIVE C C THE MINIMUM HAS NOT YET BEEN BRACKETED. C BRAKTD = .FALSE. C C SET UP XBND AS A BOUND ON THE STEP TO BE TAKEN. (XBND IS NOT COMPUTED C EXPLICITLY BUT SCXBND IS ITS SCALED VALUE.) SET THE UPPER BOUND C ON THE INTERVAL OF UNCERTAINTY INITIALLY TO XBND + TOL(XBND). C SCXBND = XBND B = SCXBND + RELTOL*DABS(SCXBND) + ABSTOL E = B + B B1 = B C C COMPUTE THE CONSTANTS REQUIRED FOR THE TWO CONVERGENCE CRITERIA. C GTEST1 = -RMU*GU gtest2 = eta*gu if (strong) gtest2 = -eta*gu c GTEST2 = -ETA*GU C C SET IENTRY TO INDICATE THAT THIS IS THE FIRST ITERATION C IENTRY = 2 GO TO 210 C C IENTRY = 2 C C UPDATE A,B,XW, AND XMIN C 20 IF (FU .GT. FMIN) GO TO 60 C C IF FUNCTION VALUE NOT INCREASED, NEW POINT BECOMES NEXT C ORIGIN AND OTHER POINTS ARE SCALED ACCORDINGLY. C CHORDU = OLDF - (XMIN + U)*GTEST1 IF (FU .LE. CHORDU) GO TO 30 C C THE NEW FUNCTION VALUE DOES NOT SATISFY THE SUFFICIENT DECREASE C CRITERION. PREPARE TO MOVE THE UPPER BOUND TO THIS POINT AND C FORCE THE INTERPOLATION SCHEME TO EITHER BISECT THE INTERVAL OF C UNCERTAINTY OR TAKE THE LINEAR INTERPOLATION STEP WHICH ESTIMATES C THE ROOT OF F(ALPHA)=CHORD(ALPHA). C CHORDM = OLDF - XMIN*GTEST1 GU = -GMIN DENOM = CHORDM-FMIN IF (DABS(DENOM) .GE. 1.D-15) GO TO 25 DENOM = 1.D-15 IF (CHORDM-FMIN .LT. 0.D0) DENOM = -DENOM 25 CONTINUE IF (XMIN .NE. ZERO) GU = GMIN*(CHORDU-FU)/DENOM FU = HALF*U*(GMIN+GU) + FMIN IF (FU .LT. FMIN) FU = FMIN GO TO 60 30 FW = FMIN FMIN = FU GW = GMIN GMIN = GU XMIN = XMIN + U A = A-U B = B-U XW = -U SCXBND = SCXBND - U IF (GU .LE. ZERO) GO TO 40 B = ZERO BRAKTD = .TRUE. GO TO 50 40 A = ZERO 50 TOL = DABS(XMIN)*RELTOL + ABSTOL GO TO 90 C C IF FUNCTION VALUE INCREASED, ORIGIN REMAINS UNCHANGED C BUT NEW POINT MAY NOW QUALIFY AS W. C 60 IF (U .LT. ZERO) GO TO 70 B = U BRAKTD = .TRUE. GO TO 80 70 A = U 80 XW = U FW = FU GW = GU 90 TWOTOL = TOL + TOL XMIDPT = HALF*(A + B) c c check termination criteria of line search c (the 2nd corresponds to the STRONG Wolfe (strong=.true.), c the first one to the regular ones (no absolute value)) c if (.not.strong) convrg = dabs(xmidpt) .le. twotol - half*(b-a) * .or. gmin .ge. gtest2 .and. fmin .lt. oldf .and. * (dabs(xmin - xbnd) .gt. tol .or. .not. braktd) if (strong) convrg = dabs(xmidpt) .le. twotol - half*(b-a) .or. * dabs(gmin) .le. gtest2 .and. fmin .lt. oldf .and. * (dabs(xmin - xbnd) .gt. tol .or. .not. braktd) if (.not. convrg) go to 100 itest = 0 if (xmin .ne. zero) return C C IF THE FUNCTION HAS NOT BEEN REDUCED, CHECK TO SEE THAT THE RELATIVE C CHANGE IN F(X) IS CONSISTENT WITH THE ESTIMATE OF THE DELTA- C UNIMODALITY CONSTANT, TOL. IF THE CHANGE IN F(X) IS LARGER THAN C EXPECTED, REDUCE THE VALUE OF TOL. C ITEST = 3 IF (DABS(OLDF-FW) .LE. FPRESN*(ONE + DABS(OLDF))) RETURN TOL = POINT1*TOL IF (TOL .LT. TNYTOL) RETURN RELTOL = POINT1*RELTOL ABSTOL = POINT1*ABSTOL TWOTOL = POINT1*TWOTOL C C CONTINUE WITH THE COMPUTATION OF A TRIAL STEP LENGTH C 100 R = ZERO Q = ZERO S = ZERO IF (DABS(E) .LE. TOL) GO TO 150 C C FIT CUBIC THROUGH XMIN AND XW C R = THREE*(FMIN-FW)/XW + GMIN + GW ABSR = DABS(R) Q = ABSR IF (GW .EQ. ZERO .OR. GMIN .EQ. ZERO) GO TO 140 C C COMPUTE THE SQUARE ROOT OF (R*R - GMIN*GW) IN A WAY C WHICH AVOIDS UNDERFLOW AND OVERFLOW. C ABGW = DABS(GW) ABGMIN = DABS(GMIN) S = DSQRT(ABGMIN)*DSQRT(ABGW) IF ((GW/ABGW)*GMIN .GT. ZERO) GO TO 130 C C COMPUTE THE SQUARE ROOT OF R*R + S*S. C SUMSQ = ONE P = ZERO IF (ABSR .GE. S) GO TO 110 C C THERE IS A POSSIBILITY OF OVERFLOW. C IF (S .GT. RTSMLL) P = S*RTSMLL IF (ABSR .GE. P) SUMSQ = ONE +(ABSR/S)**2 SCALE = S GO TO 120 C C THERE IS A POSSIBILITY OF UNDERFLOW. C 110 IF (ABSR .GT. RTSMLL) P = ABSR*RTSMLL IF (S .GE. P) SUMSQ = ONE + (S/ABSR)**2 SCALE = ABSR 120 SUMSQ = DSQRT(SUMSQ) Q = BIG IF (SCALE .LT. BIG/SUMSQ) Q = SCALE*SUMSQ GO TO 140 C C COMPUTE THE SQUARE ROOT OF R*R - S*S C 130 Q = DSQRT(DABS(R+S))*DSQRT(DABS(R-S)) IF (R .GE. S .OR. R .LE. (-S)) GO TO 140 R = ZERO Q = ZERO GO TO 150 C C COMPUTE THE MINIMUM OF FITTED CUBIC C 140 IF (XW .LT. ZERO) Q = -Q S = XW*(GMIN - R - Q) Q = GW - GMIN + Q + Q IF (Q .GT. ZERO) S = -S IF (Q .LE. ZERO) Q = -Q R = E IF (B1 .NE. STEP .OR. BRAKTD) E = STEP C C CONSTRUCT AN ARTIFICIAL BOUND ON THE ESTIMATED STEPLENGTH C 150 A1 = A B1 = B STEP = XMIDPT IF (BRAKTD) GO TO 160 STEP = -FACTOR*XW IF (STEP .GT. SCXBND) STEP = SCXBND IF (STEP .NE. SCXBND) FACTOR = FIVE*FACTOR GO TO 170 C C IF THE MINIMUM IS BRACKETED BY 0 AND XW THE STEP MUST LIE C WITHIN (A,B). C 160 IF ((A .NE. ZERO .OR. XW .GE. ZERO) .AND. (B .NE. ZERO .OR. * XW .LE. ZERO)) GO TO 180 C C IF THE MINIMUM IS NOT BRACKETED BY 0 AND XW THE STEP MUST LIE C WITHIN (A1,B1). C D1 = XW D2 = A IF (A .EQ. ZERO) D2 = B C THIS LINE MIGHT BE C IF (A .EQ. ZERO) D2 = E U = - D1/D2 STEP = FIVE*D2*(POINT1 + ONE/U)/ELEVEN IF (U .LT. ONE) STEP = HALF*D2*DSQRT(U) 170 IF (STEP .LE. ZERO) A1 = STEP IF (STEP .GT. ZERO) B1 = STEP C C REJECT THE STEP OBTAINED BY INTERPOLATION IF IT LIES OUTSIDE THE C REQUIRED INTERVAL OR IT IS GREATER THAN HALF THE STEP OBTAINED C DURING THE LAST-BUT-ONE ITERATION. C 180 IF (DABS(S) .LE. DABS(HALF*Q*R) .OR. * S .LE. Q*A1 .OR. S .GE. Q*B1) GO TO 200 C C A CUBIC INTERPOLATION STEP C STEP = S/Q C C THE FUNCTION MUST NOT BE EVALUTATED TOO CLOSE TO A OR B. C IF (STEP - A .GE. TWOTOL .AND. B - STEP .GE. TWOTOL) GO TO 210 IF (XMIDPT .GT. ZERO) GO TO 190 STEP = -TOL GO TO 210 190 STEP = TOL GO TO 210 200 E = B-A C C IF THE STEP IS TOO LARGE, REPLACE BY THE SCALED BOUND (SO AS TO C COMPUTE THE NEW POINT ON THE BOUNDARY). C 210 IF (STEP .LT. SCXBND) GO TO 220 STEP = SCXBND C C MOVE SXBD TO THE LEFT SO THAT SBND + TOL(XBND) = XBND. C SCXBND = SCXBND - (RELTOL*DABS(XBND)+ABSTOL)/(ONE + RELTOL) 220 U = STEP IF (DABS(STEP) .LT. TOL .AND. STEP .LT. ZERO) U = -TOL IF (DABS(STEP) .LT. TOL .AND. STEP .GE. ZERO) U = TOL ITEST = 1 RETURN END c c############################################################ *------------------------------------------------------------- * Date last modified: February 9, 1996 * ----------------------------------- * subroutine gtims1 (v,bv,n,x,m,mb,me,mg,mgsb,nfp,nfun, 1 c,cb,sb,fpar,c0,gc,gv,gobj,accrcy,quadr,nquad,eleps, 1 e,ge,gev,qa,qb,w1,w2,c1,cb1,e1, 1 ipgc,ipge,icgc,icge,nszc,nsze) * common /three/ beta * double precision v(n),bv(n),x(n),fpar(nfp),c0(mg), 1 c(m),gc(nszc),gv(nszc),cb(mb),beta, 1 cv,term,gobj(n),accrcy,u1,u2,u3,eleps, 1 e(me),ge(nsze),gev(nsze),qa(mg),qb(mg), 1 w1(n),w2(n),c1(m),cb1(mb),e1(me) * integer ipgc(n+1),ipge(n+1) integer icgc(nszc),icge(nsze) integer sb(mb) * logical quadr(mg) * *----------------------- Compute Hessian-vector product for * barrier function. * First compute Hessian-vector product for * nonlinear functions. * call gtims2(x,v,bv,gv,n,m,mb,me,gc,gobj, 1 w1,w2,c1,cb1,mgsb,cb,sb,accrcy, 1 e1,ge,gev,nszc,nsze) * * We needed a fct. call in gtims2, i.e. update nfun * update nquad - only statistics, not necessary for running of code * c(w1) in gtims2 is stored in w(2*n+1) and cb(w1) in w(2*n+m+1). * nfun=nfun+1 nq=0 do 5 i=1,m 5 if (c1(i).lt.-fpar(1)*beta*c0(i)) nq=nq+1 do 7 i=1,mgsb 7 if (cb1(i).lt.-fpar(1)*beta*c0(i)) nq=nq+1 if (nq.gt.0) nquad=nquad+1 * *----------------------- Compute B*v (term involving gradients only) * safeguarding necessary!! * do 30 i = 1,m u2=0.d0 u3=0.d0 * *---------------------------- Cycle 10: * Multiply row i of gc with vector v * do 10 j=1,n kp=ipgc(j) if(kp.eq.0) go to 10 ii=j 9 ku=ipgc(ii+1) if(ku.eq.0) then ii=ii+1 go to 9 end if ku=ku-1 do k=kp,ku if(icgc(k).eq.i) then u1=v(j)*gc(k) if (u1.gt.0.d0) then u2=u2+u1 else u3=u3+u1 endif go to 10 end if end do 10 continue cv=u2+u3 if (dabs(cv).lt.eleps*dabs(u2)) cv=0.D0 if (quadr(i)) then term = - fpar(1)*qa(i)* cv * fpar(i+1) else term=cv*fpar(i+1)/(c(i)/fpar(1)+c0(i))/ * (c(i)/fpar(1)+c0(i))/fpar(1) end if * *---------------------------- Cycle 20: * Multiply row i of gc with a * constant term. * do 20 j = 1,n kp=ipgc(j) if(kp.eq.0) go to 20 ii=j 11 ku=ipgc(ii+1) if(ku.eq.0) then ii=ii+1 go to 11 end if ku=ku-1 do k=kp,ku if(icgc(k).eq.i) then bv(j) = bv(j) + term*gc(k) if (dabs(bv(j)).lt.eleps*dabs(term*gc(k))) 1 bv(j)=0.d0 go to 20 end if end do 20 continue 30 continue * do 40 i = 1,mgsb if (quadr(m+i)) then term=-fpar(1)*qa(i)*v(abs(sb(i)))*fpar(m+i+1) else term=v(abs(sb(i)))*fpar(m+i+1)/(cb(i)/fpar(1)+c0(m+i))/ * (cb(i)/fpar(1)+c0(m+i))/fpar(1) endif bv(abs(sb(i))) = bv(abs(sb(i))) + term if (dabs(bv(abs(sb(i)))).lt.eleps*dabs(term)) * bv(abs(sb(i)))=0.d0 40 continue * *----------------------- Classic terms. * do 50 i=mgsb+1,mb term=fpar(1)*v(abs(sb(i)))*fpar(m+i+1)/cb(i)/cb(i) bv(abs(sb(i)))=bv(abs(sb(i)))+term if (dabs(bv(abs(sb(i)))).lt.eleps*dabs(term)) * bv(abs(sb(i)))=0.d0 50 continue * *----------------------- Compute B*v (term involving Hessian, * Hessian of classic terms = 0). * do 80 i = 1,m if (quadr(i)) then term=(qa(i)*c(i)+qb(i))*fpar(1)*fpar(1+i) else term=fpar(1+i)/(c(i)/fpar(1)+c0(i)) end if * *---------------------------- Cycle 70: * Multiply row i of gv with a * constant term. do 70 j = 1,n kp=ipgc(j) if(kp.eq.0) go to 70 ii=j 69 ku=ipgc(ii+1) if(ku.eq.0) then ii=ii+1 go to 69 end if ku=ku-1 do k=kp,ku if(icgc(k).eq.i) then bv(j) = bv(j) - term * gv(k) if(dabs(bv(j)).lt.eleps*dabs(term*gv(k))) * bv(j)=0.d0 go to 70 end if end do 70 continue 80 continue * *----------------------- Augmented Lagrangian terms, saveguard * calculation of bv. * do 120 i=1,me u2=0.d0 u3=0.d0 * *---------------------------- Cycle 90: * Multiply row i of ge with vector v. * do 90 j=1,n kp=ipge(j) if(kp.eq.0) go to 90 ii=j 89 ku=ipge(ii+1) if(ku.eq.0) then ii=ii+1 go to 89 end if ku=ku-1 do k=kp,ku if(icge(k).eq.i) then u1=v(j)*ge(k) if(u1.gt.0.d0) then u2=u2+u1 else u3=u3+u1 endif go to 90 end if end do 90 continue cv=u2+u3 if (dabs(cv).lt.eleps*dabs(u2)) cv=0.D0 cv=cv/fpar(1) term=fpar(1+m+mb+i)-e(i)/fpar(1) * *---------------------------- Cycle 100: * Multiply row i of ge, as well as * row i of gev with a constant term. * do 100 j=1,n kp=ipge(j) if(kp.eq.0) go to 100 ii=j 99 ku=ipge(ii+1) if(ku.eq.0) then ii=ii+1 go to 99 end if ku=ku-1 do k=kp,ku if(icge(k).eq.i) then bv(j)=bv(j) - term*gev(k) + cv*ge(k) if(dabs(bv(j)).lt.eleps*dabs(-term*gev(k)+ * cv*ge(k))) bv(j)=0.d0 go to 100 end if end do 100 continue 120 continue return end * c############################################################ *-------------------------------------------------------------- * Date last modified: February 7, 1996 * ----------------------------------- * This routine computes the product of the matrix g times the * vector v and stores the result in the vector gv. * (finite-difference version). * First make sure delta won't make us leave the feasibility * region defined by the simple bounds. *-------------------------------------------------------------- * subroutine gtims2(x,v,bv,gv,n,m,mb,me,gc,gobj,w1,w2,c1,cb1, 1 mgsb,cb,sb,accrcy,e1,ge,gev,nszc,nsze) * double precision v(n),bv(n),x(n),gobj(n),w1(n),w2(n), 1 gc(nszc),gv(nszc),c1(m),cb1(mb),delta,dinv,cb(mb), 1 accrcy,objf,e1(me),ge(nsze),gev(nsze) * common /nbloc/ nb integer sb(mb) * delta = 1.d-8 * call maxal(n,mb,mgsb,delta,cb,v,sb,accrcy) * dinv = 1.d0/delta do 30 i = 1,n w1(i) = x(i) + delta*v(i) 30 continue * call prob(n,m,mb,me,sb,w1,objf,w2,c1,gv,cb1,e1,gev, 1 nszc,nsze,nb) * do 40 i = 1,n bv(i) = (w2(i) - gobj(i))*dinv 40 continue * * do 60 j = 1,nszc gv(j) = (gv(j) - gc(j))*dinv 60 continue * do 70 j = 1,nsze gev(j) = (gev(j) - ge(j))*dinv 70 continue * return end * c############################################################ SUBROUTINE INITP3(DIAGB,EMAT,N,LRESET,YKSK,YRSR,BSK, * SK,YK,SR,YR,MODET,UPD1,eleps) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION DIAGB(N),EMAT(N),YKSK,YRSR,BSK(N),SK(N), * YK(N),COND,SR(N),YR(N),SDS,SRDS,YRSK,TD,D1,DN, * ddots,eleps LOGICAL LRESET,UPD1 IF (UPD1) GO TO 90 IF (LRESET) GO TO 60 DO 10 I = 1,N BSK(I) = DIAGB(I)*SR(I) 10 CONTINUE SDS = DDOTs(N,SR,1,BSK,1,eleps) SRDS = DDOTs(N,SK,1,BSK,1,eleps) YRSK = DDOTs(N,YR,1,SK,1,eleps) DO 20 I = 1,N TD = DIAGB(I) BSK(I) = TD*SK(I) - BSK(I)*SRDS/SDS+YR(I)*YRSK/YRSR EMAT(I) = TD-TD*TD*SR(I)*SR(I)/SDS+YR(I)*YR(I)/YRSR 20 CONTINUE SDS = DDOTs(N,SK,1,BSK,1,eleps) DO 30 I = 1,N EMAT(I) = EMAT(I) - BSK(I)*BSK(I)/SDS+YK(I)*YK(I)/YKSK 30 CONTINUE GO TO 110 60 CONTINUE DO 70 I = 1,N BSK(I) = DIAGB(I)*SK(I) 70 CONTINUE SDS = DDOTs(N,SK,1,BSK,1,eleps) DO 80 I = 1,N TD = DIAGB(I) EMAT(I) = TD - TD*TD*SK(I)*SK(I)/SDS + YK(I)*YK(I)/YKSK 80 CONTINUE GO TO 110 90 CONTINUE CALL DCOPY(N,DIAGB,1,EMAT,1) 110 CONTINUE IF (MODET .LT. 1) RETURN D1 = EMAT(1) DN = EMAT(1) DO 120 I = 1,N IF (EMAT(I) .LT. D1) D1 = EMAT(I) IF (EMAT(I) .GT. DN) DN = EMAT(I) 120 CONTINUE COND = DN/D1 WRITE(1,800) D1,DN,COND 800 FORMAT(' ',//8X,'DMIN =',1PD12.4,' DMAX =',1PD12.4, * ' COND =',1PD12.4,/) RETURN END C c############################################################ SUBROUTINE INITPC(DIAGB,EMAT,N,W,LW,MODET, * UPD1,YKSK,GSK,YRSR,LRESET,eleps) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION DIAGB(N),EMAT(N),W(LW),eleps DOUBLE PRECISION YKSK,GSK,YRSR LOGICAL LRESET,UPD1 COMMON/SUBSCR/ LGV,LZ1,LZK,LV,LSK,LYK,LDIAGB,LSR,LYR, * LHYR,LHG,LHYK,LPK,LEMAT,LWTEST CALL INITP3(DIAGB,EMAT,N,LRESET,YKSK,YRSR,W(LHYK), * W(LSK),W(LYK),W(LSR),W(LYR),MODET,UPD1,eleps) RETURN END C c############################################################ *------------------------------------------------------------ * Date last modified: February 7, 1996 *------------------------------------------------------------ * SUBROUTINE LINDER(N,SMALL,EPSMCH,RELTOL,ABSTOL, * TNYTOL,ETA,SFTBND,XBND,P,GTP,X,F,ALPHA,G,NFTOTL, * IFLAG,W,LW,strong, * m,mb,me,mg,mgsb,nfp,nfun,nquad, * objf,gobj,c,gc,cb,sb,c0,fpar,eleps,dell,work,quadr, * e,ge,qa,qb,qc,ipgc,ipge,icgc,icge,nszc,nsze) * INTEGER N,NFTOTL,IFLAG,LW,sb(mb) * DOUBLE PRECISION SMALL,EPSMCH,RELTOL,ABSTOL,TNYTOL,ETA, 1 SFTBND,XBND,GTP,F,ALPHA,P(N),X(N),G(N),W(LW), 1 c(m),cb(mb),c0(mg),fpar(nfp),objf,gobj(n),gc(nszc), 1 dell(mg),work(2*mg+mb-mgsb),eleps,e(me),ge(nsze), 1 qa(mg),qb(mg),qc(mg) * integer ipgc(n+1),ipge(n+1) integer icgc(nszc),icge(nsze) * INTEGER I,IENTRY,ITEST,L,LG,LX,NUMF,ITCNT DOUBLE PRECISION A,B,B1,BIG,EE,FACTOR,FMIN,FPRESN,FU, 1 FW,GMIN,GTEST1,GTEST2,GU,GW,OLDF,SCXBND,STEP, 1 TOL,U,XMIN,XW,RMU,RTSMLL,UALPHA * LOGICAL BRAKTD,quadr(mg),strong * * THE FOLLOWING STANDARD FUNCTIONS AND SYSTEM FUNCTIONS ARE * CALLED WITHIN LINDER * DOUBLE PRECISION DSQRT,ddots * * ALLOCATE THE ADDRESSES FOR LOCAL WORKSPACE * LX = 1 LG = LX + N LSPRNT = 0 NPRNT = 10000 RTSMLL = DSQRT(SMALL) BIG = 1.D0/SMALL ITCNT = 0 * * SET THE ESTIMATED RELATIVE PRECISION IN F(X). * FPRESN = 10.D0*EPSMCH NUMF = 0 U = ALPHA FU = F FMIN = F GU = GTP RMU = 1.0D-4 * * FIRST ENTRY SETS UP THE INITIAL INTERVAL OF UNCERTAINTY. * IENTRY = 1 10 CONTINUE * * TEST FOR TOO MANY ITERATIONS * ITCNT = ITCNT + 1 IFLAG = 1 IF (ITCNT .GT. 20) GO TO 50 IFLAG = 0 * CALL GETPTC(BIG,SMALL,RTSMLL,RELTOL,ABSTOL,TNYTOL, 1 FPRESN,ETA,RMU,XBND,U,FU,GU,XMIN,FMIN,GMIN, 1 XW,FW,GW,A,B,OLDF,B1,SCXBND,EE,STEP,FACTOR, 1 BRAKTD,GTEST1,GTEST2,TOL,IENTRY,ITEST,strong) * IF (LSPRNT .GE. NPRNT) CALL LSOUT(IENTRY,ITEST,XMIN,FMIN,GMIN, 1 XW,FW,GW,U,A,B,TOL,RELTOL,SCXBND,XBND) * * IF ITEST=1, THE ALGORITHM REQUIRES THE FUNCTION VALUE TO BE * CALCULATED. * IF (ITEST .NE. 1) GO TO 30 UALPHA = XMIN + U L = LX DO 20 I = 1,N W(L) = X(I) + UALPHA*P(I) L = L + 1 20 CONTINUE * * In barfct the array work is structured as follows: * work[ delg(mg) | l(mg) | lsb(mb) ]. * call barfct(n,m,mb,me,mg,mgsb,nfp,nfun,nquad, 1 w(lx),fu,w(lg),objf,gobj,c,gc,cb,sb,c0,fpar, 1 eleps,dell,work(1),work(mg+1),work(2*mg+1),quadr, 1 e,ge,qa,qb,qc,ipgc,ipge,icgc,icge,nszc,nsze) * NUMF = NUMF + 1 GU = DDOTs(N,W(LG),1,P,1,eleps) * * THE GRADIENT VECTOR CORRESPONDING TO THE BEST POINT IS * OVERWRITTEN IF FU IS LESS THAN FMIN AND FU IS SUFFICIENTLY * LOWER THAN F AT THE ORIGIN. * IF (FU .LE. FMIN .AND. FU .LE. OLDF-UALPHA*GTEST1) 1 CALL DCOPY(N,W(LG),1,G,1) GOTO 10 * * IF ITEST=2 OR 3 A LOWER POINT COULD NOT BE FOUND * 30 CONTINUE NFTOTL = NUMF IFLAG = 1 IF (ITEST .NE. 0) GO TO 50 * * IF ITEST=0 A SUCCESSFUL SEARCH HAS BEEN MADE * IFLAG = 0 F = FMIN ALPHA = XMIN DO 40 I = 1,N X(I) = X(I) + ALPHA*P(I) 40 CONTINUE * * fct call is not necessary (since without it, f and g is known), * but easier (to avoid we'd have to save together with fmin the * corresponding c's, gc's, ge's,...). * nq=nquad * call barfct(n,m,mb,me,mg,mgsb,nfp,nfun,nquad, * x,f,g,objf,gobj,c,gc,cb,sb,c0,fpar, * eleps,dell,work(1),work(mg+1),work(2*mg+1),quadr, * e,ge,qa,qb,qc,ipgc,ipge,icgc,icge,nszc,nsze) * nfun=nfun-1 nquad=nq 50 RETURN END * c############################################################ *--------------------------------------------------------------- * Date last modified: February 7, 1996. (N. Andrei) *--------------------------------------------------------------- * TRUNCATED-NEWTON METHOD: SUBROUTINES * WRITTEN BY: STEPHEN G. NASH * OPERATIONS RESEARCH AND APPLIED STATISTICS DEPT. * GEORGE MASON UNIVERSITY * FAIRFAX, VA 22030 * Remarks: * * 1. Changes made by Marc G. Breitfeld to link to PENBAR * original code in upper case, changes in lower case * most important changes: * 1) special matrix-vector product for penalty-barier function * gtims1 and gtims2 (as suggested by A. Sofer and S. Nash) * 2) safeguarded scalar product (ddots: sums up positive and * negative part separately and checks for cancellation) * 3) routine calls barfct, which calculates the penalty-barrier function * (instead of SFUN as in the original TN) * * 2. Changes made by Neculai Andrei to consider the sparsity of * the Jacobians of the inequality and equality constraints. * Only the structural nonzero elements of these Jacobians are * stored in column order. A zero column of Jacobian has zero * element in array ipgc (or ipge). * ******************************************************************* * * THIS ROUTINE SOLVES THE OPTIMIZATION PROBLEM * * MINIMIZE F(X) * X * * WHERE X IS A VECTOR OF N REAL VARIABLES. THE METHOD USED IS * A TRUNCATED-NEWTON ALGORITHM (SEE "NEWTON-TYPE MINIMIZATION VIA * THE LANCZOS METHOD" BY S.G. NASH (SIAM J. NUMER. ANAL. 21 (1984), * PP. 770-778). THIS ALGORITHM FINDS A LOCAL MINIMUM OF F(X). IT DOES * NOT ASSUME THAT THE FUNCTION F IS CONVEX (AND SO CANNOT GUARANTEE A * GLOBAL SOLUTION), BUT DOES ASSUME THAT THE FUNCTION IS BOUNDED BELOW. * IT CAN SOLVE PROBLEMS HAVING ANY NUMBER OF VARIABLES, BUT IT IS * ESPECIALLY USEFUL WHEN THE NUMBER OF VARIABLES (N) IS LARGE. * * SUBROUTINE PARAMETERS: * * IERROR - (INTEGER) ERROR CODE * ( 0 => NORMAL RETURN) * ( 2 => MORE THAN MAXFUN EVALUATIONS) * ( 3 => LINE SEARCH FAILED TO FIND * ( LOWER POINT (MAY NOT BE SERIOUS) * (-1 => ERROR IN INPUT PARAMETERS) * N - (INTEGER) NUMBER OF VARIABLES * X - (REAL*8) VECTOR OF LENGTH AT LEAST N; ON INPUT, AN INITIAL * ESTIMATE OF THE SOLUTION; ON OUTPUT, THE COMPUTED SOLUTION. * G - (REAL*8) VECTOR OF LENGTH AT LEAST N; ON OUTPUT, THE FINAL * VALUE OF THE GRADIENT * F - (REAL*8) ON INPUT, A ROUGH ESTIMATE OF THE VALUE OF THE * OBJECTIVE FUNCTION AT THE SOLUTION; ON OUTPUT, THE VALUE * OF THE OBJECTIVE FUNCTION AT THE SOLUTION * W - (REAL*8) WORK VECTOR OF LENGTH AT LEAST 14*N * LW - (INTEGER) THE DECLARED DIMENSION OF W * ETA - SEVERITY OF THE LINESEARCH * MAXFUN - MAXIMUM ALLOWABLE NUMBER OF FUNCTION EVALUATIONS * XTOL - DESIRED ACCURACY FOR THE SOLUTION X* * STEPMX - MAXIMUM ALLOWABLE STEP IN THE LINESEARCH * ACCRCY - ACCURACY OF COMPUTED FUNCTION VALUES * MSGLVL - DETERMINES QUANTITY OF PRINTED OUTPUT * 0 = NONE, 1 = ONE LINE PER MAJOR ITERATION. * MAXIT - MAXIMUM NUMBER OF INNER ITERATIONS PER STEP * *---------------------------------------------------------------- * SUBROUTINE LMQN (IFAIL, N, X, F, G, W, LW, 1 MSGLVL, MAXIT, MAXFUN, ETA, STEPMX, ACCRCY,strong, 1 XTOL,epsmch,niter,gnorm,xnorm,eleps,dell, 1 objf,gobj,fpar,c,gc,c0,sb,cb,work,nwork, 1 m,mb,me,mg,mgsb,nfp,nfun,nquad,quadr, 1 savepre,itout,e,ge,sneps,qa,qb,qc, 1 w1,w2,c1,cb1,e1,gvv,gev, 1 ipgc,ipge,icgc,icge,nszc,nsze) * INTEGER MSGLVL,N,MAXFUN,IFAIL,LW,sb(mb) * DOUBLE PRECISION X(N),G(N),W(LW),ETA,XTOL,STEPMX,F,ACCRCY, 1 eleps,dell(mg),eps,e(me),ge(nsze),sneps, 1 objf,gobj(n),fpar(nfp),c(m),gc(nszc),c0(mg),cb(mb),work(nwork), 1 qa(mg),qb(mg),qc(mg), 1 w1(n),w2(n),c1(m),cb1(mb),e1(me),gvv(nszc),gev(nsze) * integer ipgc(n+1),ipge(n+1) integer icgc(nszc),icge(nsze) * logical quadr(mg),strong,savepre * * THIS ROUTINE IS A TRUNCATED-NEWTON METHOD. * THE TRUNCATED-NEWTON METHOD IS PRECONDITIONED BY A LIMITED-MEMORY * QUASI-NEWTON METHOD (THIS PRECONDITIONING STRATEGY IS DEVELOPED * IN THIS ROUTINE) WITH A FURTHER DIAGONAL SCALING (SEE ROUTINE NDIA3). * FOR FURTHER DETAILS ON THE PARAMETERS, SEE ROUTINE TN. * INTEGER I, ICYCLE, IOLDG, IPK, IYK, LOLDG, LPK, LSR, 1 LWTEST, LYK, LYR, NFTOTL, NITER, NM1, NUMF, NWHY DOUBLE PRECISION ABSTOL, ALPHA, DIFNEW, DIFOLD, EPSMCH, 1 EPSRED, FKEEP, FM, FNEW, FOLD, FSTOP, FTEST, GNORM, GSK, 1 GTG, GTPNEW, OLDF, OLDGTP, ONE, PE, PEPS, PNORM, RELTOL, 1 RTEPS, RTLEPS, RTOL, RTOLSQ, SMALL, SPE, TINY, 1 TNYTOL, TOLEPS, XNORM, YKSK, YRSR, ZERO LOGICAL LRESET, UPD1 * * THE FOLLOWING IMSL AND STANDARD FUNCTIONS ARE USED * DOUBLE PRECISION DABS, DSQRT, STEP1, DNRM2,ddots COMMON /SUBSCR/ LGV,LZ1,LZK,LV,LSK,LYK,LDIAGB,LSR,LYR, 1 LOLDG,LHG,LHYK,LPK,LEMAT,LWTEST * * INITIALIZE PARAMETERS AND CONSTANTS * eps=xtol IF (MSGLVL .GE. -2) WRITE(1,800) CALL SETPAR(N) UPD1 = .TRUE. IRESET = 0 NFEVAL = 0 NMODIF = 0 NLINCG = 0 FSTOP = F ZERO = 0.D0 ONE = 1.D0 NM1 = N - 1 * * WITHIN THIS ROUTINE THE ARRAY W(LOLDG) IS SHARED BY W(LHYR) * LHYR = LOLDG * * CHECK PARAMETERS AND SET CONSTANTS * CALL CHKUCP(LWTEST,MAXFUN,NWHY,N,ALPHA,EPSMCH, 1 ETA,PEPS,RTEPS,RTOL,RTOLSQ,STEPMX,FTEST, 1 XTOL,XNORM,X,LW,SMALL,TINY,ACCRCY) * IF (NWHY .LT. 0) GO TO 120 * CALL SETUCR(SMALL,NFTOTL,NITER,N,F,FNEW, 1 FM,GTG,OLDF,G,X, 1 m,mb,me,mg,mgsb,nfp,nfun,nquad, 1 objf,gobj,c,gc,cb,sb,c0,fpar,eleps,dell,work(1),quadr, 1 e,ge,qa,qb,qc,ipgc,ipge,icgc,icge,nszc,nsze) * FOLD = FNEW IF (MSGLVL .GE. 1) WRITE(1,810) NITER,NFTOTL,NLINCG,FNEW,GTG * * CHECK FOR SMALL GRADIENT AT THE STARTING POINT. * test as in PENBAR with modified Newton is added * FTEST = ONE + DABS(FNEW) c IF (GTG .LT. 1.D-4*EPSMCH*FTEST*FTEST) GO TO 90 IF ((GTG .LT. 1.D-4*EPSMCH*FTEST*FTEST).or. 1 (gtg.le.eps*eps*dmax1(xnorm*xnorm,1.d0))) GO TO 90 * * SET INITIAL VALUES TO OTHER PARAMETERS * ICYCLE = NM1 TOLEPS = RTOL + RTEPS RTLEPS = RTOLSQ + EPSMCH GNORM = DSQRT(GTG) DIFNEW = ZERO EPSRED = 5.0D-2 FKEEP = FNEW * * SET THE DIAGONAL OF THE APPROXIMATE HESSIAN TO UNITY. * if (itout.eq.1 .or. .not.savepre) then IDIAGB = LDIAGB DO 10 I = 1,N W(IDIAGB) = ONE IDIAGB = IDIAGB + 1 10 CONTINUE end if * * ..................START OF MAIN ITERATIVE LOOP.......... * * COMPUTE THE NEW SEARCH DIRECTION * MODET = MSGLVL - 3 * CALL MODLNP(MODET,W(LPK),W(LGV),W(LZ1),W(LV), * W(LDIAGB),W(LEMAT),X,G,W(LZK), * N,W,LW,NITER,MAXIT,NFEVAL,NMODIF, * NLINCG,UPD1,YKSK,GSK,YRSR,LRESET,.FALSE.,IPIVOT, * ACCRCY,GTPNEW,GNORM,XNORM,gobj,nfun,nquad,quadr,eleps, * m,mb,me,mg,mgsb,nfp,c,cb,sb,gc,fpar,c0, * e,ge,qa,qb,w1,w2,c1,cb1,e1,gvv,gev, * ipgc,ipge,icgc,icge,nszc,nsze) * 20 CONTINUE * CALL DCOPY(N,G,1,W(LOLDG),1) PNORM = DNRM2(N,W(LPK),1) OLDF = FNEW OLDGTP = GTPNEW * * PREPARE TO COMPUTE THE STEP LENGTH * PE = PNORM + EPSMCH * * COMPUTE THE ABSOLUTE AND RELATIVE TOLERANCES FOR THE LINEAR SEARCH * RELTOL = RTEPS*(XNORM + ONE)/PE ABSTOL = - EPSMCH*FTEST/(OLDGTP - EPSMCH) * * COMPUTE THE SMALLEST ALLOWABLE SPACING BETWEEN POINTS IN * THE LINEAR SEARCH * spe is the maximal steplength * TNYTOL = EPSMCH*(XNORM + ONE)/PE spe = STEPMX/PE call maxal(n,mb,mgsb,spe,cb,w(lpk),sb,accrcy) * * SET THE INITIAL STEP LENGTH. * ALPHA = STEP1(FNEW,FM,OLDGTP,SPE,epsmch) if (alpha.gt.spe) alpha=spe * * PERFORM THE LINEAR SEARCH * CALL LINDER(N,SMALL,EPSMCH,RELTOL,ABSTOL,TNYTOL, 1 ETA,ZERO,SPE,W(LPK),OLDGTP,X,FNEW,ALPHA,G,NUMF, 1 NWHY,W,LW,strong, 1 m,mb,me,mg,mgsb,nfp,nfun,nquad, 1 objf,gobj,c,gc,cb,sb,c0,fpar,eleps,dell,work(1),quadr, 1 e,ge,qa,qb,qc,ipgc,ipge,icgc,icge,nszc,nsze) * FOLD = FNEW NITER = NITER + 1 NFTOTL = NFTOTL + NUMF GTG = DDOTs(N,G,1,G,1,eleps) IF (MSGLVL .GE. 1) WRITE(1,810) NITER,NFTOTL,NLINCG,FNEW,GTG IF (NWHY .LT. 0) GO TO 120 IF (NWHY .EQ. 0 .OR. NWHY .EQ. 2) GO TO 30 * * THE LINEAR SEARCH HAS FAILED TO FIND A LOWER POINT * NWHY = 3 GO TO 100 30 IF (NWHY .LE. 1) GO TO 40 * call barfct(n,m,mb,me,mg,mgsb,nfp,nfun,nquad, 1 x,fnew,g,objf,gobj,c,gc,cb,sb,c0,fpar, 1 eleps,dell,work(1),work(mg+1),work(2*mg+1),quadr, 1 e,ge,qa,qb,qc,ipgc,ipge,icgc,icge,nszc,nsze) * NFTOTL = NFTOTL + 1 * * TERMINATE IF MORE THAN MAXFUN EVALUTATIONS HAVE BEEN MADE * 40 NWHY = 2 IF (NFTOTL .GT. MAXFUN) GO TO 110 NWHY = 0 * * SET UP PARAMETERS USED IN CONVERGENCE AND RESETTING TESTS * DIFOLD = DIFNEW DIFNEW = OLDF - FNEW * * IF THIS IS THE FIRST ITERATION OF A NEW CYCLE, COMPUTE THE * PERCENTAGE REDUCTION FACTOR FOR THE RESETTING TEST. * IF (ICYCLE .NE. 1) GO TO 50 IF (DIFNEW .GT. 2.0D0 *DIFOLD) EPSRED = EPSRED + EPSRED IF (DIFNEW .LT. 5.0D-1*DIFOLD) EPSRED = 5.0D-1*EPSRED 50 CONTINUE GNORM = DSQRT(GTG) FTEST = ONE + DABS(FNEW) XNORM = DNRM2(N,X,1) * * TEST FOR CONVERGENCE - included: norm of search vector * and the criteria used in PENBAR with the modified Newton * IF ((ALPHA*PNORM .LT. TOLEPS*(ONE + XNORM) 1 .AND. DABS(DIFNEW) .LT. RTLEPS*FTEST 1 .AND. GTG .LT. PEPS*FTEST*FTEST) 1 .OR. GTG .LT. 1.D-4*ACCRCY*FTEST*FTEST 1 .or. gtg.le.eps*eps*dmax1(xnorm*xnorm,1.d0) 1 .or. (pnorm/xnorm).le.sneps) goto 90 * * COMPUTE THE CHANGE IN THE ITERATES AND THE CORRESPONDING CHANGE * IN THE GRADIENTS * ISK = LSK IPK = LPK IYK = LYK IOLDG = LOLDG DO 60 I = 1,N W(IYK) = G(I) - W(IOLDG) W(ISK) = ALPHA*W(IPK) IPK = IPK + 1 ISK = ISK + 1 IYK = IYK + 1 IOLDG = IOLDG + 1 60 CONTINUE * * SET UP PARAMETERS USED IN UPDATING THE DIRECTION OF SEARCH. * YKSK = DDOTs(N,W(LYK),1,W(LSK),1,eleps) LRESET = .FALSE. IF (ICYCLE .EQ. NM1 .OR. DIFNEW .LT. 1 EPSRED*(FKEEP-FNEW)) LRESET = .TRUE. IF (LRESET) GO TO 70 YRSR = DDOTs(N,W(LYR),1,W(LSR),1,eleps) IF (YRSR .LE. ZERO) LRESET = .TRUE. 70 CONTINUE UPD1 = .FALSE. * * COMPUTE THE NEW SEARCH DIRECTION * MODET = MSGLVL - 3 CALL MODLNP(MODET,W(LPK),W(LGV),W(LZ1),W(LV), * W(LDIAGB),W(LEMAT),X,G,W(LZK), * N,W,LW,NITER,MAXIT,NFEVAL,NMODIF, * NLINCG,UPD1,YKSK,GSK,YRSR,LRESET,.FALSE.,IPIVOT, * ACCRCY,GTPNEW,GNORM,XNORM,gobj,nfun,nquad,quadr,eleps, * m,mb,me,mg,mgsb,nfp,c,cb,sb,gc,fpar,c0, * e,ge,qa,qb,w1,w2,c1,cb1,e1,gvv,gev, * ipgc,ipge,icgc,icge,nszc,nsze) * IF (LRESET) GO TO 80 * * STORE THE ACCUMULATED CHANGE IN THE POINT AND GRADIENT AS AN * "AVERAGE" DIRECTION FOR PRECONDITIONING. * CALL DXPY(N,W(LSK),1,W(LSR),1) CALL DXPY(N,W(LYK),1,W(LYR),1) ICYCLE = ICYCLE + 1 GOTO 20 * * RESET * 80 IRESET = IRESET + 1 * * INITIALIZE THE SUM OF ALL THE CHANGES IN X. * CALL DCOPY(N,W(LSK),1,W(LSR),1) CALL DCOPY(N,W(LYK),1,W(LYR),1) FKEEP = FNEW ICYCLE = 1 GO TO 20 * * ...............END OF MAIN ITERATION....................... * 90 IFAIL = 0 F = FNEW RETURN 100 OLDF = FNEW * * LOCAL SEARCH HERE COULD BE INSTALLED HERE * 110 F = OLDF * * SET IFAIL * 120 IFAIL = NWHY RETURN 800 FORMAT(/' NIT NF CG', 10X, 'F', 21X, 'GTG',/) 810 FORMAT(' ',I4,1X,I5,1X,I5,1X,1PD22.15,2X,1PD15.8) END * c############################################################ SUBROUTINE LSOUT(ILOC,ITEST,XMIN,FMIN,GMIN,XW,FW,GW,U,A, * B,TOL,EPS,SCXBD,XLAMDA) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION XMIN,FMIN,GMIN,XW,FW,GW,U,A,B, * TOL,EPS,SCXBD,XLAMDA C C ERROR PRINTOUTS FOR GETPTC C DOUBLE PRECISION YA,YB,YBND,YW,YU YU = XMIN + U YA = A + XMIN YB = B + XMIN YW = XW + XMIN YBND = SCXBD + XMIN WRITE(1,800) WRITE(1,810) TOL,EPS WRITE(1,820) YA,YB WRITE(1,830) YBND WRITE(1,840) YW,FW,GW WRITE(1,850) XMIN,FMIN,GMIN WRITE(1,860) YU WRITE(1,870) ILOC,ITEST RETURN 800 FORMAT(///' OUTPUT FROM LINEAR SEARCH') 810 FORMAT(' TOL AND EPS'/2D25.14) 820 FORMAT(' CURRENT UPPER AND LOWER BOUNDS'/2D25.14) 830 FORMAT(' STRICT UPPER BOUND'/D25.14) 840 FORMAT(' XW, FW, GW'/3D25.14) 850 FORMAT(' XMIN, FMIN, GMIN'/3D25.14) 860 FORMAT(' NEW ESTIMATE'/2D25.14) 870 FORMAT(' ILOC AND ITEST'/2I3) END C c############################################################ *------------------------------------------------------------- * Date last modified: January 4, 1996 * ----------------------------------- * Calculates the maximal alpha such that the simple bounds are * NOT violated. *------------------------------------------------------------- * subroutine maxal(n,mb,mgsb,almax,cb,s,sb,accrcy) * double precision almax,cb(mb),s(n),temp double precision accrcy,eps,delta,alini integer sb(mb) * * delta is how close we go to the singularity * delta=.9995d0 eps=1.D4*accrcy alini=almax do 10 i=mgsb+1,mb if (sb(i).gt.0 .and. s(sb(i)).lt.0.D0) then temp=-cb(i)/s(sb(i)) if (cb(i).lt.eps) then s(sb(i))=0.D0 else if (almax.gt.temp) then almax=temp endif else if (sb(i).lt.0 .and. s(-sb(i)).gt.0.D0) then temp=cb(i)/s(-sb(i)) if (cb(i).lt.eps) then s(-sb(i))=0.D0 else if (almax.gt.temp) then almax=temp endif endif 10 continue if (almax.le.0.D0) write(1,*) almax,'= almax <= 0' * if (almax.lt.alini) almax=delta*almax * return end * c############################################################ *------------------------------------------------------------ * Date last modified: February 7, 1996 *------------------------------------------------------------ * SUBROUTINE MODLNP(MODET,ZSOL,GV,R,V,DIAGB,EMAT, * X,G,ZK,N,W,LW,NITER,MAXIT,NFEVAL,NMODIF,NLINCG, * UPD1,YKSK,GSK,YRSR,LRESET,BOUNDS,IPIVOT,ACCRCY, * GTP,GNORM,XNORM,gobj,nfun,nquad,quadr,eleps, * m,mb,me,mg,mgsb,nfp,c,cb,sb,gc,fpar,c0, * e,ge,qa,qb,w1,w2,c1,cb1,e1,gvv,gev, 1 ipgc,ipge,icgc,icge,nszc,nsze) c INTEGER MODET,N,NITER,IPIVOT(1),sb(mb) * DOUBLE PRECISION ZSOL(N),G(N),GV(N),R(N),V(N),DIAGB(N),W(LW), * EMAT(N),ZK(N),X(N),ACCRCY,ALPHA,BETA1,GSK,GTP,PR, * QOLD,QNEW,QTEST,RHSNRM,RZ,RZOLD,TOL,VGV,YKSK,YRSR, * GNORM,XNORM,c(m),cb(mb),gc(nszc),fpar(nfp),c0(mg), * gobj(n),eleps,e(me),ge(nsze),qa(mg),qb(mg), * w1(n),w2(n),c1(m),cb1(mb),e1(me), * gvv(nszc),gev(nsze) * integer ipgc(n+1),ipge(n+1) integer icgc(nszc),icge(nsze) * DOUBLE PRECISION ddots LOGICAL FIRST,UPD1,LRESET,BOUNDS,quadr(mg) C C THIS ROUTINE PERFORMS A PRECONDITIONED CONJUGATE-GRADIENT C ITERATION IN ORDER TO SOLVE THE NEWTON EQUATIONS FOR A SEARCH C DIRECTION FOR A TRUNCATED-NEWTON ALGORITHM. WHEN THE VALUE OF THE C QUADRATIC MODEL IS SUFFICIENTLY REDUCED, C THE ITERATION IS TERMINATED. C C PARAMETERS C C MODET - INTEGER WHICH CONTROLS AMOUNT OF OUTPUT C ZSOL - COMPUTED SEARCH DIRECTION C G - CURRENT GRADIENT C GV,GZ1,V - SCRATCH VECTORS C R - RESIDUAL C DIAGB,EMAT - DIAGONAL PRECONDITONING MATRIX C NITER - NONLINEAR ITERATION # C FEVAL - VALUE OF QUADRATIC FUNCTION C C ************************************************************* C INITIALIZATION C ************************************************************* C C GENERAL INITIALIZATION C IF (MODET .GT. 0) WRITE(1,800) IF (MAXIT .EQ. 0) RETURN FIRST = .TRUE. RHSNRM = GNORM TOL = 1.D-12 QOLD = 0.D0 C C INITIALIZATION FOR PRECONDITIONED CONJUGATE-GRADIENT ALGORITHM C CALL INITPC(DIAGB,EMAT,N,W,LW,MODET, * UPD1,YKSK,GSK,YRSR,LRESET,eleps) DO 10 I = 1,N R(I) = -G(I) V(I) = 0.D0 ZSOL(I) = 0.D0 10 CONTINUE C C ************************************************************ C MAIN ITERATION C ************************************************************ C DO 30 K = 1,MAXIT NLINCG = NLINCG + 1 IF (MODET .GT. 1) WRITE(1,810) K C C CG ITERATION TO SOLVE SYSTEM OF EQUATIONS C IF (BOUNDS) then do i=1,n if(ipivot(i).ne.0) r(i)=0.d0 end do end if * CALL MSOLVE(R,ZK,N,W,LW,UPD1,YKSK,GSK, * YRSR,LRESET,FIRST,eleps) IF (BOUNDS) then do i=1,n if(ipivot(i).ne.0) zk(i)=0.d0 end do end if * RZ = DDOTs(N,R,1,ZK,1,eleps) IF (RZ/RHSNRM .LT. TOL) GO TO 80 IF (K .EQ. 1) BETA1 = 0.D0 IF (K .GT. 1) BETA1 = RZ/RZOLD DO 20 I = 1,N V(I) = ZK(I) + BETA1*V(I) 20 CONTINUE * IF (BOUNDS) then do i=1,n if(ipivot(i).ne.0) v(i)=0.d0 end do end if * CALL GTIMS1(v,gv,n,x,m,mb,me,mg,mgsb,nfp,nfun, * c,cb,sb,fpar,c0,gc,gvv,gobj, * accrcy,quadr,nquad,eleps, * e,ge,gev,qa,qb,w1,w2,c1,cb1,e1, 1 ipgc,ipge,icgc,icge,nszc,nsze) * IF (BOUNDS) then do i=1,n if(ipivot(i).ne.0) gv(i)=0.d0 end do end if * NFEVAL = NFEVAL + 1 VGV = DDOTs(N,V,1,GV,1,eleps) IF (VGV/RHSNRM .LT. TOL) GO TO 50 CALL NDIA3(N,EMAT,V,GV,R,VGV,MODET,eleps) C C COMPUTE LINEAR STEP LENGTH C ALPHA = RZ / VGV IF (MODET .GE. 1) WRITE(1,820) ALPHA C C COMPUTE CURRENT SOLUTION AND RELATED VECTORS C CALL DAXPY(N,ALPHA,V,1,ZSOL,1) CALL DAXPY(N,-ALPHA,GV,1,R,1) C C TEST FOR CONVERGENCE C GTP = DDOTs(N,ZSOL,1,G,1,eleps) PR = DDOTs(N,R,1,ZSOL,1,eleps) QNEW = 5.D-1 * (GTP + PR) QTEST = K * (1.D0 - QOLD/QNEW) IF (QTEST .LT. 0.D0) GO TO 70 QOLD = QNEW IF (QTEST .LE. 5.D-1) GO TO 70 C C PERFORM CAUTIONARY TEST C IF (GTP .GT. 0) GO TO 40 RZOLD = RZ 30 CONTINUE C C TERMINATE ALGORITHM C K = K-1 GO TO 70 C C TRUNCATE ALGORITHM IN CASE OF AN EMERGENCY C 40 IF (MODET .GE. -1) WRITE(1,830) K CALL DAXPY(N,-ALPHA,V,1,ZSOL,1) GTP = DDOTs(N,ZSOL,1,G,1,eleps) GO TO 90 50 CONTINUE IF (MODET .GT. -2) WRITE(1,840) 60 IF (K .GT. 1) GO TO 70 CALL MSOLVE(G,ZSOL,N,W,LW,UPD1,YKSK,GSK,YRSR,LRESET,FIRST,eleps) * do i=1,n zsol(i)=-zsol(i) end do * IF (BOUNDS) then do i=1,n if(ipivot(i).ne.0) zsol(i)=0.d0 end do end if * GTP = DDOTs(N,ZSOL,1,G,1,eleps) 70 CONTINUE IF (MODET .GE. -1) WRITE(1,850) K GO TO 90 80 CONTINUE IF (MODET .GE. -1) WRITE(1,860) IF (K .GT. 1) GO TO 70 CALL DCOPY(N,G,1,ZSOL,1) * do i=1,n zsol(i)=-zsol(i) end do * IF (BOUNDS) then do i=1,n if(ipivot(i).ne.0) zsol(i)=0.d0 end do end if * GTP = DDOTs(N,ZSOL,1,G,1,eleps) GO TO 70 C C STORE (OR RESTORE) DIAGONAL PRECONDITIONING C 90 CONTINUE CALL DCOPY(N,EMAT,1,DIAGB,1) RETURN 800 FORMAT(' ',//,' ENTERING MODLNP') 810 FORMAT(' ',//,' ### ITERATION ',I2,' ###') 820 FORMAT(' ALPHA',1PD16.8) 830 FORMAT(' G(T)Z POSITIVE AT ITERATION ',I2, * ' - TRUNCATING METHOD',/) 840 FORMAT(' ',10X,'HESSIAN NOT POSITIVE-DEFINITE') 850 FORMAT(' ',/,8X,'MODLAN TRUNCATED AFTER ',I3,' ITERATIONS') 860 FORMAT(' PRECONDITIONING NOT POSITIVE-DEFINITE') END C c############################################################ SUBROUTINE MSLV(G,Y,N,SK,YK,DIAGB,SR,YR,HYR,HG,HYK, * UPD1,YKSK,GSK,YRSR,LRESET,FIRST,eleps) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION G(N),Y(N) C C THIS ROUTINE ACTS AS A PRECONDITIONING STEP FOR THE C LINEAR CONJUGATE-GRADIENT ROUTINE. IT IS ALSO THE C METHOD OF COMPUTING THE SEARCH DIRECTION FROM THE C GRADIENT FOR THE NON-LINEAR CONJUGATE-GRADIENT CODE. C IT REPRESENTS A TWO-STEP SELF-SCALED BFGS FORMULA. C DOUBLE PRECISION YKSK,GSK,YRSR,RDIAGB,YKHYK,GHYK, * YKSR,YKHYR,YRHYR,GSR,GHYR,ddots,eleps DOUBLE PRECISION SK(N),YK(N),DIAGB(N),SR(N),YR(N),HYR(N),HG(N), * HYK(N),ONE LOGICAL LRESET,UPD1,FIRST IF (UPD1) GO TO 100 ONE = 1.D0 GSK = DDOTs(N,G,1,SK,1,eleps) IF (LRESET) GO TO 60 C C COMPUTE HG AND HY WHERE H IS THE INVERSE OF THE DIAGONALS C DO 57 I = 1,N RDIAGB = 1.0D0/DIAGB(I) HG(I) = G(I)*RDIAGB IF (FIRST) HYK(I) = YK(I)*RDIAGB IF (FIRST) HYR(I) = YR(I)*RDIAGB 57 CONTINUE IF (FIRST) YKSR = DDOTs(N,YK,1,SR,1,eleps) IF (FIRST) YKHYR = DDOTs(N,YK,1,HYR,1,eleps) GSR = DDOTs(N,G,1,SR,1,eleps) GHYR = DDOTs(N,G,1,HYR,1,eleps) IF (FIRST) YRHYR = DDOTs(N,YR,1,HYR,1,eleps) CALL SSBFGS(N,ONE,SR,YR,HG,HYR,YRSR, * YRHYR,GSR,GHYR,HG) IF (FIRST) CALL SSBFGS(N,ONE,SR,YR,HYK,HYR,YRSR, * YRHYR,YKSR,YKHYR,HYK) YKHYK = DDOTs(N,HYK,1,YK,1,eleps) GHYK = DDOTs(N,HYK,1,G,1,eleps) CALL SSBFGS(N,ONE,SK,YK,HG,HYK,YKSK, * YKHYK,GSK,GHYK,Y) RETURN 60 CONTINUE C C COMPUTE GH AND HY WHERE H IS THE INVERSE OF THE DIAGONALS C DO 65 I = 1,N RDIAGB = 1.D0/DIAGB(I) HG(I) = G(I)*RDIAGB IF (FIRST) HYK(I) = YK(I)*RDIAGB 65 CONTINUE IF (FIRST) YKHYK = DDOTs(N,YK,1,HYK,1,eleps) GHYK = DDOTs(N,G,1,HYK,1,eleps) CALL SSBFGS(N,ONE,SK,YK,HG,HYK,YKSK, * YKHYK,GSK,GHYK,Y) RETURN 100 CONTINUE DO 110 I = 1,N 110 Y(I) = G(I) / DIAGB(I) RETURN END C c############################################################ SUBROUTINE MSOLVE(G,Y,N,W,LW,UPD1,YKSK,GSK, * YRSR,LRESET,FIRST,eleps) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION G(N),Y(N),W(LW),YKSK,GSK,YRSR,eleps LOGICAL UPD1,LRESET,FIRST C C THIS ROUTINE SETS UPT THE ARRAYS FOR MSLV C COMMON/SUBSCR/ LGV,LZ1,LZK,LV,LSK,LYK,LDIAGB,LSR,LYR, * LHYR,LHG,LHYK,LPK,LEMAT,LWTEST CALL MSLV(G,Y,N,W(LSK),W(LYK),W(LDIAGB),W(LSR),W(LYR),W(LHYR), * W(LHG),W(LHYK),UPD1,YKSK,GSK,YRSR,LRESET,FIRST,eleps) RETURN END c c############################################################ SUBROUTINE NDIA3(N,E,V,GV,R,VGV,MODET,eleps) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION E(N),V(N),GV(N),R(N),VGV,VR,ddots, * eleps C C UPDATE THE PRECONDITIOING MATRIX BASED ON A DIAGONAL VERSION C OF THE BFGS QUASI-NEWTON UPDATE. C VR = DDOTs(N,V,1,R,1,eleps) DO 10 I = 1,N E(I) = E(I) - R(I)*R(I)/VR + GV(I)*GV(I)/VGV IF (E(I) .GT. 1.D-6) GO TO 10 IF (MODET .GT. -2) WRITE(1,800) E(I) E(I) = 1.D0 10 CONTINUE RETURN 800 FORMAT(' *** EMAT NEGATIVE: ',1PD16.8) END c############################################################ *------------------------------------------------------------- * Date last modified: January 4, 1996 *------------------------------------ * Calculation of the coefficients of the quadratic approximation * to the logarithmic terms in the modified barrier function, i.e. * q(t) = 1/2*qa*t^2 + qb*t + qc. * Remark: the vectors qa, qb, and qc need to be updated whenever * mu has changed. *-------------------------------------------------------------- * subroutine quad(mu,c0,mg,qa,qb,qc) * common /three/ beta * double precision beta,qa(mg),qb(mg),qc(mg),mu,z,c0(mg) * z=(1.D0-beta)*(1.D0-beta) * do 10 i=1,mg qa(i)=-1.D0/mu/mu/c0(i)/c0(i)/z qb(i)=(1.D0-2.D0*beta)/mu/c0(i)/z qc(i)=beta*(2.D0-3.D0*beta)/2.D0/z+dlog(c0(i)*(1.D0-beta)) 10 continue return end * c############################################################ SUBROUTINE SETPAR(N) IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER LSUB(14) COMMON/SUBSCR/ LSUB,LWTEST C C SET UP PARAMETERS FOR THE OPTIMIZATION ROUTINE C DO 10 I = 1,14 LSUB(I) = (I-1)*N + 1 10 CONTINUE LWTEST = LSUB(14) + N - 1 RETURN END C c############################################################ *----------------------------------------------------------- * Date last modified: February 7, 1996 *----------------------------------------------------------- * SUBROUTINE SETUCR(SMALL,NFTOTL,NITER,N,F,FNEW, * FM,GTG,OLDF,G,X, * m,mb,me,mg,mgsb,nfp,nfun,nquad, * objf,gobj,c,gc,cb,sb,c0,fpar,eleps,dell,work,quadr, * e,ge,qa,qb,qc,ipgc,ipge,icgc,icge,nszc,nsze) c IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER sb(mb) DOUBLE PRECISION F,FNEW,FM,GTG,OLDF,SMALL,G(N),X(N), * objf,gobj(n),c(m),gc(nszc),cb(mb),c0(mg),fpar(nfp), * eleps,dell(mg),work(2*mg+mb-mgsb),ddots,e(me),ge(nsze), * qa(mg),qb(mg),qc(mg) * integer ipgc(n+1),ipge(n+1) integer icgc(nszc),icge(nsze) * logical quadr(mg) C C CHECK INPUT PARAMETERS, COMPUTE THE INITIAL FUNCTION VALUE, SET C CONSTANTS FOR THE SUBSEQUENT MINIMIZATION C c changed: f is not set yet!! c FM = F C C COMPUTE THE INITIAL FUNCTION VALUE C call barfct(n,m,mb,me,mg,mgsb,nfp,nfun,nquad, * x,fnew,g,objf,gobj,c,gc,cb,sb,c0,fpar, * eleps,dell,work(1),work(mg+1),work(2*mg+1),quadr, * e,ge,qa,qb,qc,ipgc,ipge,icgc,icge,nszc,nsze) * c changed: now, fnew is set FM = fnew NFTOTL = 1 C C SET CONSTANTS FOR LATER C NITER = 0 OLDF = FNEW GTG = DDOTs(N,G,1,G,1,eleps) RETURN END C c############################################################ subroutine sing(n,m,mb,mg,mgsb,nfp,c,c0,cb,sb,fpar,gc, * alpha,s,dg1,eleps) c c Estimation of a first trial step such that the log stays defined c (by Murray/Wright '92) c double precision c(m),gc(m*n),c0(mg),cb(mb), 1 fpar(nfp),s(n),eleps,u1,u2,u3, 1 dgmin,dg,sigmin,sig,dgsigm,dg1,alpha,atmp logical gencon,sbgen integer sb(mb) c c Calculate directional derivatives of c along search vector s. If they c are all nonnegative, we take initial step alpha and return. If not, we c find the max. alpha via a lin. approx. of c. c We memorize with gencon whether the minimum occured at a general c constraint or a simple bound. sbgen is true if the minimum occured c at a simple bound that was NOT handled separately. c c write(1,*) 'sing' gencon=.true. sbgen=.false. dgmin=0.D0 sigmin=1.D8 do 20 i=1,m u2=0.D0 u3=0.D0 do 10 j=1,n u1=gc(i+(j-1)*m)*s(j) if (u1.gt.0.D0) then u2=u2+gc(i+(j-1)*m)*s(j) else u3=u3+gc(i+(j-1)*m)*s(j) endif 10 continue if (dabs(u2+u3).lt.eleps*dabs(u2)) then u2=0.D0 u3=0.D0 endif dg=u2+u3 if (dg.ge.0.) goto 20 if (dg.lt.dgmin) dgmin=dg sig=(-fpar(1)*c0(i)-c(i))/dg if (sig.lt.sigmin) then sigmin=sig dgsigm=dg imin=i else if (sig.eq.sigmin) then if (dg.lt.dgsigm) then dgsigm=dg imin=i endif endif 20 continue do 30 i=1,mgsb if (sb(i).gt.0) then dg=s(sb(i)) else dg=-s(-sb(i)) endif if (dg.ge.0.) goto 30 if (dg.lt.dgmin) dgmin=dg sig=(-fpar(1)*c0(i)-cb(i))/dg if (sig.lt.sigmin) then sigmin=sig dgsigm=dg imin=i sbgen=.true. else if (sig.eq.sigmin) then if (dg.lt.dgsigm) then dgsigm=dg imin=i sbgen=.true. endif endif 30 continue do 40 i=mgsb+1,mb if (sb(i).gt.0) then dg=s(sb(i)) else dg=-s(-sb(i)) endif if (dg.ge.0.) goto 40 if (dg.lt.dgmin) dgmin=dg sig=-cb(i)/dg if (sig.lt.sigmin) then sigmin=sig dgsigm=dg imin=i gencon=.false. else if (sig.eq.sigmin) then if (dg.lt.dgsigm) then dgsigm=dg imin=i gencon=.false. endif endif 40 continue c c If all directional derivatives nonnegative, return c if (dgmin.ge.0.D0) return c c Estimate the line-minimum of the barrier function by only considering c constraint imin and using the linear approximation of f and c. c dg1 is the directional derivative of the objective. c 50 continue if (.not.gencon) then if (dgsigm.ne.0.D0) u1=fpar(m+imin+1)/dg1 if ((u1.lt.0.D0).and.(dgsigm.ne.0.D0)) then atmp=sigmin+fpar(1)*u1 else atmp=.9D0*sigmin endif else if (sbgen) imin=imin+m if (dgsigm.ne.0.D0) u1=fpar(imin+1)/dg1-c0(imin)/dgsigm if ((u1.lt.0.D0).and.(dgsigm.ne.0.D0)) then atmp=sigmin+fpar(1)*u1 else atmp=.9D0*sigmin endif endif if (atmp.lt.1.D-14) atmp=alpha c c If the suggested trial step is less than the initial step from 'ucmin' c then take it. c if (atmp.lt.alpha) then alpha=atmp endif return end c c############################################################ SUBROUTINE SSBFGS(N,GAMMA,SJ,YJ,HJV,HJYJ,YJSJ,YJHYJ, * VSJ,VHYJ,HJP1V) IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER N DOUBLE PRECISION GAMMA,YJSJ,YJHYJ,VSJ,VHYJ DOUBLE PRECISION SJ(N),YJ(N),HJV(N),HJYJ(N),HJP1V(N) C C SELF-SCALED BFGS C INTEGER I DOUBLE PRECISION BETA1,DELTA DELTA = (1.D0 + GAMMA*YJHYJ/YJSJ)*VSJ/YJSJ * - GAMMA*VHYJ/YJSJ BETA1 = -GAMMA*VSJ/YJSJ DO 10 I = 1,N HJP1V(I) = GAMMA*HJV(I) + DELTA*SJ(I) + BETA1*HJYJ(I) 10 CONTINUE RETURN END C c############################################################ DOUBLE PRECISION FUNCTION STEP1(FNEW,FM,GTP,SMAX,epsmch) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION FNEW,FM,GTP,SMAX C C ******************************************************** C STEP1 RETURNS THE LENGTH OF THE INITIAL STEP TO BE TAKEN C ALONG THE VECTOR P IN THE NEXT LINEAR SEARCH. C ******************************************************** C DOUBLE PRECISION ALPHA,D,epsmch DOUBLE PRECISION DABS D = DABS(FNEW-FM) ALPHA = 1.D0 IF (2.D0*D .LE. (-GTP) .AND. D .GE. EPSMCH) 1 ALPHA = -2.D0*D/GTP IF (ALPHA .GE. SMAX) ALPHA = SMAX STEP1 = ALPHA RETURN END * c############################################################ * Date last modified: February 9, 1996 * ----------------------------------- * Stop, i.e.set nterm=1 as soon as v is small enough. * v is the max of 1) 1st order of the original problem * 2) complementarity * 3) infeasibility. * Or we stop if the fractional change in the objective is small * enough. * * The subroutine considres the sparsity of the Jacobians. * ------------------------------------------------------------- * subroutine stop(n,m,mb,me,mg,mgsb,nfp,itout, * c,cb,sb,objf,gobj,gc,cmin,fpar,nterm,tau,xn,gn,objfp, * dell,eleps,fcrit,e,ge,lamup,mup,muupp,lo, 1 ipgc,ipge,icgc,icge,nszc,nsze) * double precision c(m),gc(nszc),gobj(n),mup, 1 dell(mg),cb(mb),xn,u1,u2,u3,gn, 1 objf,fpar(nfp),cmin,compl,tau,v,eleps, 1 ord1,ord1n,objfp,fraco,tau2,e(me),ge(nsze),emax * integer sb(mb) integer ipgc(n+1),ipge(n+1) integer icgc(nszc),icge(nsze) * logical fcrit,lamup,muupp * double precision dabs * * Use the updated multipliers for the complementarity and first order. * If we do not update the multipliers, the non-updated ones are used. * fcrit=.true. * *----------------------- xn = Norm of vector x xn=xn+1.d0 * *----------------------- cmin = Minimum value of inequalities * constraints including bounds of * variables. cmin=1.D6 do 5 i=1,m 5 if (cmin.gt.c(i)) cmin=c(i) do 10 i=1,mb 10 if (cmin.gt.cb(i)) cmin=cb(i) if (m+mb.eq.0) cmin=0.d0 * *----------------------- emax = Maximum absolute value of * equalities. * emax=0.d0 do 15 i=1,me 15 if (emax.lt.dabs(e(i))) emax=dabs(e(i)) * *----------------------- compl = Value of complementary * conditions. * compl=0.D0 if (lamup) then * *----------------------- Update de Lagrange multipliers. * do 20 i=1,mg fpar(i+1)=fpar(i+1)*dell(i) compl=compl+fpar(i+1)*dabs(c(i)) 20 continue * do 25 i=mgsb+1 , mb fpar(m+1+i)=fpar(m+1+i)*mup/cb(i) compl=compl+fpar(m+1+i)*dabs(cb(i)) 25 continue * do 30 i=1,me fpar(1+m+mb+i)=fpar(1+m+mb+i)-e(i)/mup 30 compl=compl+dabs(fpar(1+m+mb+i)*e(i)) * else do 33 i=1,m 33 compl=compl+fpar(i+1)*dabs(c(i)) do 35 i=1,mb 35 compl=compl+fpar(m+1+i)*dabs(cb(i)) do 38 i=1,me 38 compl=compl+dabs(fpar(1+m+mb+i)*e(i)) endif * *----------------------- ord1n = Value of the first order * optimality conditions. * ord1n=0.D0 do 50 i=1,n u2=0.D0 u3=0.D0 kp=ipgc(i) if(kp.eq.0) go to 41 ii=i 39 ku=ipgc(ii+1) if(ku.eq.0) then ii=ii+1 go to 39 end if ku=ku-1 do 40 k=kp,ku j=icgc(k) u1=-fpar(j+1)*gc(k) if (u1.gt.0.d0) then u2=u2+u1 else u3=u3+u1 endif 40 continue * 41 do 45 j=1,mb if (sb(j).eq.i) then u2=u2-fpar(1+m+j) else if (sb(j).eq.-i) then u3=u3+fpar(1+m+j) endif 45 continue * kp=ipge(i) if(kp.eq.0) go to 48 ii=i 46 ku=ipge(ii+1) if(ku.eq.0) then ii=ii+1 go to 46 end if ku=ku-1 do 47 k=kp,ku j=icge(k) u1=-fpar(1+m+mb+j)*ge(k) if (u1.gt.0.d0) then u2=u2+u1 else u3=u3+u1 endif 47 continue 48 continue if (dabs(u2+u3).lt.eleps*dabs(u2)) then u2=0.d0 u3=0.d0 endif ord1=gobj(i)+u2+u3 if (dabs(gobj(i)+u2+u3).lt.eleps*dabs(u2+u3)) ord1=0.D0 50 if (dabs(ord1).gt.ord1n) ord1n=dabs(ord1) * fraco=dabs(objf-objfp)/(1.d0+dabs(objfp)) if (itout.eq.1) fraco=1.d0 tau2=1.d-2*tau if (tau2.lt.1.d-8) tau2=1.d-8 * * The fractional-change-in-the-objective criteria may only be * used if the penalty-barrier parameter was reduced from the * previous to this iteration. * if (muupp.and.(emax.le.tau).and.(-cmin.le.tau).and. * (fraco.le.tau2)) nterm=1 * v=dmax1(emax,-cmin,compl/xn,ord1n/xn) * if (v.le.tau) then nterm=1 fcrit=.false. endif * if (lo.gt.1) then write(1,100) 100 format(/,5x,'Stopping criteria (scaled by norm of x):') write(1,110) ord1n/xn 110 format(5x,'First order conditions: ',e20.13) write(1,120) compl/xn 120 format(5x,'Complementarity conditions: ',e20.13) write(1,130) cmin 130 format(5x,'Minimum value of inequalities: ',e20.13) write(1,140) emax 140 format(5x,'Maximum abs. value of equalities: ',e20.13) write(1,150) gn/xn 150 format(5x,'Scaled norm of penbar function: ',e20.13) endif return end * c############################################################ * Date assembled: July 31, 1998. * * Dr. Neculai Andrei * Research Institute for Informatics * Center for Advanced Modeling and Optimization * 8-10, Averescu Avenue, Bucharest 1, Romania * E-mail: nandrei@ici.ro c############################################################ * * * Last Line of SPENBAR (line 4158)