c version albaali&fletcher step size selection c one sided powell-wolfe test c************************************************************ c initialization program for bfgs-minimization c************************************************************ c copyright = spellucci implicit none integer nx,nsmall,maxit parameter (nx=300,nsmall=34,maxit=4000) double precision x,d,gradx,fx,xnorm,dnorm,sig,x0,d0, 1 gradx0,fx0,x0norm,d0norm,sig0,difx,ftest,dirder common/xdat/x(nx),d(nx),gradx(nx),fx,xnorm,dnorm,sig, 1 x0(nx),d0(nx),gradx0(nx),fx0,x0norm,d0norm,sig0, 2 difx(nx),ftest,dirder double precision stepterm,sigsm,sigla,alpha,beta,delta,theta, 1 dscal,kappa,tl common/stepparam/stepterm,sigsm,sigla,alpha,beta,delta,theta, 1 dscal,kappa,tl double precision epsx,epsg,epsf logical test1,test2,test3,test4 common/param/epsx,epsg,epsf,test1,test2,test3,test4 double precision epsmac,tolmac,deldif common/mach/epsmac,tolmac,deldif integer n common/dim/n integer icf,icgf,iup,inoup,ires,imod common/count/icf,icgf,iup,inoup,ires,imod double precision a,diaga,a0,diaga0,accinf,z,y,tracea,traceb common/bfgsdat/a(nx,nx),diaga(nx),a0(nsmall,nsmall), 1 diaga0(nsmall),accinf(maxit,14),z(nx),y(nx) 2 ,tracea,traceb logical analyt common/gradinf/analyt integer niter,itermax double precision eta,thet2,bfgsterm common/bfgsparam/niter,itermax,eta,thet2,bfgsterm c for testing purposes only this kind of fdat double precision xst integer icase common/fdat/xst(nx),icase character*80 ident character*30 xifile,xofile common/ident/ident,xifile,xofile integer anz,wdh,inx,ih1,ih2,ih3,ih4 real cputim,cpu logical std,single character*12 list character*1 char real statf,statgf,statcpu double precision term,tol1 integer statfail,statit,statup,statno,statre integer i,j,k,idum double precision vecnor external vecnor save statfail=0 statit=0 statup=0 statno=0 statre=0 statf=0. statgf=0. statcpu=0. epsmac = 2.d0**(-20) 100 continue epsmac=epsmac/2.d0 term=1.d0+epsmac if ( term .ne. 1.d0 ) goto 100 epsmac=epsmac+epsmac tolmac=epsmac 200 continue tol1=tolmac tolmac=tolmac/16.d0 if ( tolmac .ne. 0.d0 ) goto 200 tolmac=tol1 do i=1,maxit do j=1,14 accinf(i,j)=0.d0 enddo enddo write(*,*) 'minimizer bfgs2 version spellucci 3/96' write(*,*) 'single test ? otherwise case from file ' read(*,*) single if ( .not. single ) then write(*,*) 'give name of inputfile (e.g. bfgssteu.dat)' read(*,*) list endif if ( single ) then write(*,*) 'input case ' read(*,*) icase write(*,*) ' parameter std? if t then from block data or setup2' read(*,*) std call setup2 c*********** parameter from block-data if not defined here if ( std ) then epsx=1.d-7 epsg=epsx epsf=1.d-10 eta=epsmac**.5d0 sigsm=1.d-10 sigla=2048.d0 beta=2.d0 alpha=.1d0 delta=.01d0 theta=.99d0 deldif=dexp(dlog(epsmac)/7.d0)*.25d0 itermax=maxit kappa=0.9d0 cold kappa=0.1d0 tl=0.1d0 else write(*,*) 'input epsx' read(*,*) epsx write(*,*) 'input epsg' read(*,*) epsg write(*,*) 'input epsf' read(*,*) epsf write(*,*) 'input sigsm' read(*,*) sigsm write(*,*) 'input sigla' read(*,*) sigla beta=2.d0 alpha=.1d0 write(*,*) 'input delta' read(*,*) delta theta=.99d0 deldif=dexp(dlog(epsmac)/7.d0)*.25d0 itermax=maxit cold kappa=0.1d0 cc kappa=0.9d0 write(*,*) 'input kappa' read(*,*) kappa tl=0.1d0 write(*,*) 'input test1, test2, test3, test4 :t or f' read(*,*) test1,test2,test3,test4 endif anz=1 write(*,*) 'input x:0 keyboard,1 blockdata or setup, 2 file' read(5,*) inx if ( inx .eq. 2 ) then write(*,*) 'input filename for x-input in std-format!' read(5,*) xifile open(20,file=xifile,status='old') read(20,*) c******* skip ident-text do i=1,n read(20,*) xst(i) enddo close(20) else if ( inx .eq. 0 ) then do i=1,n write(*,*) 'x',i,' =?' read(5,*) xst(i) enddo endif endif xofile='bfgs2'//ident(1:6)//'.prot' open(60,file=xofile,status='unknown') if ( test4 ) then open(50,file='bfgs2kurzres',status='unknown') write(50,fmt='('' it '',5x,''x1'',10x,''x2'',8x, f ''n(grad)'',5x,''trace_a'',5x,''trace_b'')') endif else open(30,file=list,status='old') open(40,file='bfgs2res.lis',status='new') write(40, f fmt='('' pure bfgs with albaali-fletcher-stepsize '')') anz=1000 endif do wdh = 1, anz if ( .not. single ) then read(30,*,end=4000) icase epsx=1.d-7 epsg=epsx epsf=1.d-10 eta=epsmac**.6666d0 sigsm=1.d-10 sigla=2048.d0 beta=2.d0 alpha=.1d0 delta=.01d0 kappa=0.9d0 cold kappa=0.1d0 tl=0.1d0 theta=.99d0 deldif=dexp(dlog(epsmac)/7.d0)*.25d0 itermax=maxit c********** setup2 may change parameters and xst call setup2 xofile='bfgs2'//ident(1:6)//'.prot' open(60,file=xofile,status='unknown') if ( wdh .eq. 1 ) then write(40,fmt='('' control parameters :'','// 1 '/1x,'' epsmac='',d12.5,/1x,'' tolmac='',d12.5,'// 2 '/1x,'' deldif='',d12.5,'// 3 '/1x,'' epsx='',d12.5,/1x,'' epsg='',d12.5,'// 4 '/1x,'' epsf='',d12.5,'// 5 '/1x,'' eta='',d12.5,/1x,'' sigmin='',d12.5,'// 6 '/1x,'' sigla='',d12.5,'// 7 '/1x,'' facx='',d12.5,/1x,'' kappa='',d12.5,'// 8 '/1x,'' delta='',d12.5,'// 9 '/1x,'' tl='',d12.5,/1x,'' facf='',d12.5)') a epsmac,tolmac,deldif,epsx,epsg,epsf,eta, b sigsm,sigla,beta,kappa,delta,tl,2.d0 c write(40,fmt='(1x,'' n ident-text optimalvalue '','// 1 ''' n(delta_x) n(grad) '','// 2 ''' cond(a) spb iter '','// 3 '''upda noup rest mod #f #grad time'')') endif endif write(60,500) epsmac,tolmac,deldif,epsx,epsg,epsf,eta, 1 sigsm,sigla,beta,alpha,delta,theta c epsx = required precision in x (relative) c itermax = 30 * n is sufficient as a rule c test1 : if true, a short iteration protocol is printed at c teh end of the run. the information is stored in accinf: c accinf, see below c test2 = .true. more detailed information during the run c c test3 = .true. print quasi-newton matrix each step c c the variable icase allows to run several examples in sequence c if f(x,fx) depends on icase icf=0 icgf=0 iup=0 inoup=0 ires=0 imod=0 cpu=cputim(idum) call bfgs2(0.d0) c start may change control parameters c param 0 for bfgs means reenter = false cpu=(cputim(idum)-cpu) if ( .not. analyt ) icf=icf-n*icgf*6 if ( .not. single ) then if ( traceb .gt. 0.d0 ) then char='+' else char='-' endif write(40,fmt='(1x,i3,1x,a12,1x,d21.14,3(1x,d11.4),1x,a1, 2 5i5,2i6,f8.2)') 3 n,ident,accinf(niter,2),accinf(niter,3),accinf(niter,4), 4 accinf(niter,8),char, 5 niter,iup,inoup,ires,imod,icf,icgf,cpu if ( bfgsterm .lt. 0 f .and. vecnor(1,n,gradx) .gt. 1.d-3) then statfail=statfail+1 else statf=real(icf)+statf statgf=statgf+real(icgf) statcpu=statcpu+cpu statit=statit+niter statup=statup+iup statno=statno+inoup statre=statre+ires endif endif write(60,fmt='(/,''example '',a40, 1 /1x,''termination '',d8.1,/1x,''fxopt = '',d21.14)') 2 ident,bfgsterm,fx do i=1,n write(60,fmt='(1x,''x('',i3,'')='',d21.14, 1 '' grad('',i3,'')='',d21.14)') i,x(i),i,gradx(i) enddo c icf = number of function calls (not counted those used c for the gradient c icgf = number of gradient computations (requiring 6*n function c values each) write(60,3000) icf,icgf write(60, fmt='(1x,''steps='',i5,'' updates='',i5, f '' noupdates='',i5, f '' restarts='',i5,'' modsteps='',i5)') f niter,iup,inoup,ires,imod c if ( bfsgterm .lt. 0.d0 ) test1=.true. if ( test1 ) then do j=1,niter ih1=accinf(j,1) ih2=accinf(j,12) ih3=accinf(j,13) ih4=accinf(j,14) write(60,2000) ih1,(accinf(j,k),k=2,11),ih2,ih3,ih4 enddo endif close(60) c short information in accinf: c it = stepnumber , f = function value, delx = norm of change in x c ngrad = norm of gradient, (nonmonotonic) , c n(a)= norm of quasi-newton-matrix, approxiamtes norm of Hessian c thet2 <<1 means that the function is strongly nonconvex c sig = stepsize x(new)=x(old)-sig*direction (d), c sig <<1 means that iterations stays far from solution c dscal < 1 if the computed d is considered too long (x changes to c drastically) c n(d)=norm of d, icf+ = number of f-evaluations for this step c In good cases 1,2. icf+ >>1 means a bad d c that is a bad (possibly restart would be better, may choose a greater eta) c or x far from solution. c restart = 1 means reinitialization a=einheitsmatrix*factor c update = -1 means no update computed this step c ( f stongly nonconvex at least in the current region) c termination : bfgsterm = 0 , 1 means desired accuracy attained c =2 directional derivative in the roundoff level of f c possibly epsx , epsg too small c -2 itermax steps needed without success c -1 no acceptable stepsize even after restart c -3 nonrecoverable error in the evaluation of f or gradf c (set by user) enddo 4000 continue if ( .not. single ) then write(40,fmt='('' number of failures = '',i10)') statfail write(40,fmt='('' sum of #call f= '',f10.2)') statf write(40,fmt='('' sum of #call g= '',f10.2)') statgf write(40,fmt='('' sum of cputime= '',f10.2)') statcpu write(40,fmt='('' sum of steps= '',i10)') statit write(40,fmt='('' sum of # updates = '',i10)') statup write(40,fmt='('' sum of # noupdates= '',i10)') statno write(40,fmt='('' sum of # restarts = '',i10)') statre close(30) close(40) endif if ( test4 ) close(50) stop 'ende' 500 format(/1x,'control parameters :', 1 /1x,'epsmac=',d12.5,/1x,'tolmac=',d12.5,/1x,' deldif=', 2 d12.5,/1x,' epsx=',d12.5,/1x,' epsg=',d12.5, 2 /1x,' epsf=',d12.5, 3 /1x,' eta=',d12.5,/1x,' sigsm=',d12.5,/1x,' sigla=', 4 d12.5,/1x,' beta=',d12.5,/1x,' alpha=',d12.5,/1x,' delta=', 5 d12.5,/1x,' theta=',d12.5) 2000 format(1x,'it=',i8,' f=',d10.4,' delx=',d10.4, 1 ' ngrad=',d10.4,/12x,' n(a)=',d10.4,' thet2=',d10.4, 2 ' dird=',d10.4,/12x,' cond=',d10.4, 3 ' sig=',d10.4,' dscal=',d10.4,/12x,' n(d)=', 4 d10.4,' icf+=',i10, 5 ' new=',i2,' upd=',i2) 3000 format(/1x,'function calls ',i11,' gradient calls ', 1 i10) end c**************************************************************** c determination of stepsize by powell-wolfe descnet test c albaali-fletcher like algorithm c**************************************************************** c copyright = spellucci subroutine step implicit none integer nx,nsmall,maxit parameter (nx=300,nsmall=34,maxit=4000) double precision x,d,gradx,fx,xnorm,dnorm,sig,x0,d0, 1 gradx0,fx0,x0norm,d0norm,sig0,difx,ftest,dirder common/xdat/x(nx),d(nx),gradx(nx),fx,xnorm,dnorm,sig, 1 x0(nx),d0(nx),gradx0(nx),fx0,x0norm,d0norm,sig0, 2 difx(nx),ftest,dirder double precision stepterm,sigsm,sigla,alphag,betag,theta,dscal, 1 delta,kappa,tl common/stepparam/stepterm,sigsm,sigla,alphag,betag,delta,theta, 1 dscal,kappa,tl double precision epsx,epsg,epsf logical test1,test2,test3,test4 common/param/epsx,epsg,epsf,test1,test2,test3,test4 double precision epsmac,tolmac,deldif common/mach/epsmac,tolmac,deldif integer n common/dim/n integer icf,icgf,iup,inoup,ires,imod common/count/icf,icgf,iup,inoup,ires,imod double precision a,diaga,a0,diaga0,accinf,z,y,tracea,traceb common/bfgsdat/a(nx,nx),diaga(nx),a0(nsmall,nsmall), 1 diaga0(nsmall),accinf(maxit,14),z(nx),y(nx) 2 ,tracea,traceb logical analyt common/gradinf/analyt integer niter,itermax double precision eta,thet2,bfgsterm common/bfgsparam/niter,itermax,eta,thet2,bfgsterm c******************************************************************** c n = dimension c x = current point c d = direction c gradx = gradient of f at x c x,d,gradx = input c x0,d0,gradx0 etc. information from previous step c (saved automatically) c xnorm,dnorm = length of x or d c stepterm = 1 success , =-1 otherwise c sig = computed stepsize c we assume that 1 is asymptotically exact c otherwise you must change initialization of sig c sigsm = minimal stepsize allowed c sigla = maximal stpsize allowed c alpha = maximal stepsize reduction allowed per step c delta = armijo parameter c beta = max. growth factor of xnorm for sig=1 c********************************************************************* logical err include 'bfgs_cons.inc' c**** computation of stepsize sig given x and direction of descent c**** d. gradf=gradient of f at x. values with index "0" are the corresponding c**** values at the previous point. there holds c**** x=x0+d0*sig0. we use fletcher's criterion c**** c**** f(x)-f(x+sig*d) >= delta*sig*gradf(x)*d c**** c**** and c**** c**** abs(gradf(x+sig*d)'d) <= -kappa*gradf(x)*d c**** c**** if no sig in [sigsm,sigla] is found , we return c**** with abb=-1 unsuccessfully. abb=3 means that including sig=sigla c**** the first, but not the second criterion has been satisfied. c**** the new point may be accepted, but no update takes place then. c**** this case is possible, if a direction of negative curvature has been used c**** and sigla was not large enough. c**** termination with abb=2 means, that the directional derivative is in c**** the order of magnitude of the ronding errors in evaluating f. c**** termination abb=1 means normal return with c**** a new acepted x. c**** if an error ocured during a call to f or gradf, c**** we replace d by .1*d and repeat. c**** if this is unsuccessful, we return with abb=-2. c**** abb = 0 means a too small d, which doesn't alter x much enough with sig=1. c**** abb = 4 means that we could not satisfy both conditions concerning sig c**** although an interval including a potential solution has been diminished c**** to length espx. possibly kappa is too small. remember kappa >> delta. c**** abb = -3 means directional derivative nonnegative. c**** facf = proportionalityfactor for requirement of descent c**** in the initial search-phase. =2 here c**** facf*dirder corresponds to ny in al baali & fletcher c**** smallx = describes zero-neighborhood of x c**** delta,kappa = parameter of the descent test. c**** tl = truncation-parameter during interval-reduction-phase c**** 0 < tl << .5 c**** dscal = scaling factor for d, 1. normally. c**** however if ||d|| is large compared with ||x||, we reduce ||d|| by dscal c**** dirder = gradf*d*dscal c**** dnorm = norm of d c**** xnorm = norm of x c**** xdifn = norm of x-x0, if abb=1 c**** abb = termination reason, s.a. c**** cont = t : continue interval reduction c**** desc = t : continue searching downhill c**** egb = t : directional derivative at c**** beta (b) has been evaluated. c**** [alpha,beta] is the current interval containing the potential solution sigma c**** bundef = t : beta is not yet defined (initial search phase) c**** err = t : erro in evaluating f or gradf occured. c**** here we assume this being the case because of ||d|| or sig too large. c**** part2 = t : evaluate the gradient at x+gamma*d c**** gundef = t : gradient at xloc not yet evaluated integer i,k logical cont,desc,egb,bundef,part2,gundef double precision dloc(nx),xloc(nx),gam0, f gam1,alpha,beta,gamma,dnmax double precision fgam0, fgam1,gradxl(nx),fgamma, f ggamma,falpha,galpha double precision fbeta,gbeta,gamnew,sigst,siglal double precision my,abb,gnorm,smallx double precision gam,fgam,ggam,xdifn common/trans/gam,fgam,ggam,abb parameter (smallx=one) double precision skal1,vecnor external skal1,vecnor save c**** c**** alpha=zero gamma=zero siglal=sigla my=one err=.false. abb=one sigst=one c******* to avoid erroneous termination with abb =4 beta=siglal c******* dnorm=vecnor(1,n,d) xnorm=vecnor(1,n,x) gnorm=vecnor(1,n,gradx) dscal=one dnmax=betag*(xnorm+smallx) if ( dnorm .le. sqrt(eta)*gnorm ) then siglal=siglal*(gnorm/dnorm)*sqrt(eta) endif dirder=skal1(1,n,gradx,d) if ( dirder .ge. zero ) then abb = -three goto 10000 endif if ( -dirder .le. tp3*epsmac*(abs(fx)+epsmac) ) then if ( err ) then abb = -two else abb = two endif goto 10000 endif if ( dnorm .gt. dnmax ) then c**** d too large sigst=min(sigst,dnmax/dnorm) dscal=dnmax/dnorm endif 100 continue c**** restartpoint if error in evaluating f or gradf occured if ( dnorm*dscal .le. epsmac*(xnorm+epsmac) ) then if ( err ) then abb = -two else abb = zero endif goto 10000 endif err=.false. cont=.true. egb=.false. falpha=fx galpha=dirder alpha=zero gamma=sigst bundef=.true. c********* main iteration for interval reduction do while (cont) gam0=alpha gam1=gamma fgam0=falpha k=1 desc=.true. gundef=.true. part2=.true. do while (desc) c**** search downhill with f decreasing fast c**** if ( k .gt. 1 ) then c**** choose a new gamma if ( bundef ) then c**** quadratic interpolation using falpha,galpha,fgam0 call quadl(alpha,falpha,galpha,gam1,fgam1,gamnew) gam1=max(gamnew,min(2.d0*gam1, f siglal)) else if ( fx-fbeta .ge. delta*beta*(-dirder) ) then call quadlr(alpha,falpha,gam0,fgam0,beta,fbeta,gamnew) else if ( egb ) then call quadr(gam0,fgam0,beta,fbeta,gbeta,gamnew) else call quadlr(alpha,falpha,gam0,fgam0,beta,fbeta, f gamnew) endif endif gam1=max(alpha+(beta-alpha)*tl, f min(gamnew,beta-(beta-alpha)*tl)) endif endif do i=1,n xloc(i)=x(i)+gam1*d(i) enddo c****** call f(xloc,fgam1,err) if ( err ) then sigst=min(sigst*tm1,gam1*tm1) siglal=min(siglal*tm1,gam1*p5) beta=siglal bundef=.true. gundef=.true. dscal=sigst goto 100 endif gam=gam1 fgam=fgam1 call info(6) c****** if ( fx-fgam1 .lt. delta*gam1*(-dirder) .or. fgam1 .ge. f falpha ) then beta=gam1 bundef=.false. egb=.false. fbeta=fgam1 desc=.false. part2=.false. call quadl(alpha,falpha,galpha,beta,fbeta,gamnew) gamma=max(alpha+(beta-alpha)*tl,min(gamnew, f beta-(beta-alpha)*tl)) elseif ( fgam1 .le. fgam0+dirder*2.d0 .and. gam1 .ne. sigst f ) then gam0=gam1 fgam0=fgam1 k=k+1 else desc=.false. endif if ( gam1 .eq. siglal ) then part2=.false. gundef=.true. gamma=gam1 fgamma=fgam1 abb=three desc=.false. cont=.false. endif enddo if ( part2 ) then gamma=gam1 fgamma=fgam1 c***** if ( analyt ) then call gradf(xloc,gradxl,err) else call gradfn(xloc,gradxl,err) endif if ( err ) then sigst=min(sigst*tm1,gamma*tm1) siglal=min(siglal*tm1,gamma*p5) beta=siglal bundef=.true. gundef=.true. dscal=sigst goto 100 endif gundef=.false. c***** ggamma=skal1(1,n,gradxl,d) c**** for a convex function the directional derivative is isotone and we c**** always have gamma>alpha my=ggamma/dirder gam=gamma ggam=my call info(7) if ( ggamma .ge. kappa*dirder ) then cont=.false. elseif ( ggamma .lt. kappa*dirder ) then if ( bundef ) then if ( gamma .ge. siglal ) then abb=three cont=.false. else call cubic(alpha,falpha,galpha,gamma,fgamma,ggamma, f gamnew) endif alpha=gamma falpha=fgamma galpha=ggamma gamma=min(siglal,max(two*gamma,gamnew)) else if ( egb ) then call cubic(gamma,fgamma,ggamma,beta,fbeta,gbeta,gamnew) else call quadl(gamma,fgamma,ggamma,beta,fbeta,gamnew) endif alpha=gamma falpha=fgamma galpha=ggamma gamma=max(gamma+(beta-gamma)*tl,min(gamnew, f beta-(beta-gamma)*tl)) endif else beta=gamma fbeta=fgamma gbeta=ggamma egb=.true. bundef=.false. call cubic(alpha,falpha,galpha,gamma,fgamma,ggamma,gamnew) gamma=max(alpha+(gamma-alpha)*tl,min(gamnew, f gamma-(gamma-alpha)*tl)) endif endif if ( abs(beta-alpha) .le. f epsx*(abs(alpha)+abs(beta)) .and. cont ) then cont=.false. abb=four endif if ( .not. bundef .and. beta .lt. sigsm ) then abb=-one goto 10000 endif enddo c******************************************* if ( abb .le. zero .or. abb .eq. four ) then stepterm=-one else stepterm=one endif 9000 continue call info(9) if ( gundef ) then if ( analyt ) then call gradf(xloc,gradxl,err) else call gradfn(xloc,gradxl,err) endif if ( err ) then stepterm=-1.d0 sig=zero my=one return endif endif sig0=sig sig=gamma fx0=fx fx=fgamma do i=1,n x0(i)=x(i) x(i)=xloc(i) difx(i)=x(i)-x0(i) d0(i)=d(i) gradx0(i)=gradx(i) gradx(i)=gradxl(i) dloc(i)=x(i)-x0(i) enddo xdifn=vecnor(1,n,dloc) return 10000 continue if (abb .le. zero .or. abb .eq. four ) then stepterm=-one sig=zero return endif stepterm=one goto 9000 end c******************************************************************* c scalar product of two vectors or parts of vectors c******************************************************************* c copyright = spellucci double precision function skal1(i,j,a,b) implicit none integer i,j,k double precision a(*),b(*) double precision s if ( i .gt. j ) then skal1=0.d0 return else s=0.d0 do k=i,j s=s+a(k)*b(k) enddo skal1=s return endif end c******************************************************************* c bfgs-method using modified update and cholesky-decomposition c inaccurate line search possible c******************************************************************* subroutine bfgs2(reenter) implicit none double precision reenter integer nx,nsmall,maxit,i,j parameter (nx=300,nsmall=34,maxit=4000) double precision x,d,gradx,fx,xnorm,dnorm,sig,x0,d0, 1 gradx0,fx0,x0norm,d0norm,sig0,difx,ftest,dirder common/xdat/x(nx),d(nx),gradx(nx),fx,xnorm,dnorm,sig, 1 x0(nx),d0(nx),gradx0(nx),fx0,x0norm,d0norm,sig0, 2 difx(nx),ftest,dirder double precision stepterm,sigsm,sigla,alpha,beta,theta,dscal, 1 delta,kappa,tl common/stepparam/stepterm,sigsm,sigla,alpha,beta,delta,theta, 1 dscal,kappa,tl double precision epsx,epsg,epsf logical test1,test2,test3,test4 common/param/epsx,epsg,epsf,test1,test2,test3,test4 double precision epsmac,tolmac,deldif common/mach/epsmac,tolmac,deldif integer n common/dim/n integer icf,icgf,iup,inoup,ires,imod common/count/icf,icgf,iup,inoup,ires,imod double precision a,diaga,a0,diaga0,accinf,z,y,tracea,traceb common/bfgsdat/a(nx,nx),diaga(nx),a0(nsmall,nsmall), 1 diaga0(nsmall),accinf(maxit,14),z(nx),y(nx) 2 ,tracea,traceb logical analyt common/gradinf/analyt integer niter,itermax double precision eta,thet2,bfgsterm common/bfgsparam/niter,itermax,eta,thet2,bfgsterm c for testing purposes only this kind of fdat integer icase double precision xst common/fdat/xst(nx),icase character*80 ident character*30 xifile,xofile common/ident/ident,xifile,xofile double precision v(nx),s1,s2,s3,s4,s5,difxin,wi,vecnor, 1 skal1,den1,h double precision dif1,dif2 logical restart,err integer icsf,icfold,icgold,ifail external vecnor,skal1 save c******* initialization ********************* icsf=0 if ( reenter .eq. 0.d0 ) then do i=1,n x(i)=xst(i) x0(i)=x(i) difx(i)=1.d0+dabs(x(i)) c************ to avoid termination in testx in step one enddo xnorm=vecnor(1,n,x) call setup1 c******** setup1 initializes the quasi-newton-matrix a endif c******** allows for modification of problem parameters sig0=1.d0 d0norm=1.d0 c******** no stepsize larger than 1 in the first step niter=0 thet2=1.d0 bfgsterm=0.d0 do i=1,n d0(i)=0.d0 enddo err=.false. restart=.false. call f(x,fx,err) if ( err ) then bfgsterm=-3.d0 return endif err=.false. call gradf(x,gradx,err) if ( err ) then bfgsterm=-3.d0 return endif c****************************************************** c main iteration loop starts here c****************************************************** 100 continue niter=niter+1 if ( niter .gt. itermax ) then bfgsterm=-2.d0 return endif do i=1,14 accinf(niter,i)=0.d0 enddo accinf(niter,1)=niter accinf(niter,2)=fx accinf(niter,3)=vecnor(1,n,difx) accinf(niter,4)=vecnor(1,n,gradx) call info(1) s1=diaga(1) s2=a(1,1) do i=2,n s1=s1+diaga(i) s2=dmin1(s2,a(i,i)) enddo c***** s1 = spur (a), s2=minimales element auf der diagonale des cholesky c***** faktors von a if ( niter .gt. 1 ) then if ( (tracea*epsmac .gt. 1.d0 .and. traceb*epsmac .gt. 1.d0) f .or. s1*epsmac .gt. s2**2 .or. f accinf(niter-1,8) .gt. 1.d0 /eta ) then call refresh(niter) if ( bfgsterm .eq. -3.d0 ) return restart=.true. else restart=.false. endif endif call testx c******** check accuracy if ( bfgsterm .eq. 1.d0 ) return accinf(niter,13)=-1.d0 200 continue call cholsol(a,n,gradx,d) c******** d = quasi-newton direction (uphill) dnorm=vecnor(1,n,d) if ( dnorm .eq. 0.d0 ) then bfgsterm=0.d0 return endif dirder=skal1(1,n,d,gradx) accinf(niter,7)=dirder if ( delta*dirder .le. epsmac*(dabs(fx)+1.d0)*100.d0 ) then bfgsterm=2.d0 return endif call info(2) c c icfold=icf icgold=icgf c**** this step-routine requires d downhill do i=1,n d(i)=-d(i) enddo call step do i=1,n d(i)=-d(i) enddo c c accinf(niter,9)=sig accinf(niter,10)=dscal accinf(niter,11)=dnorm accinf(niter,12)=icf-icfold if ( .not. analyt ) f accinf(niter,12)=accinf(niter,12)-(icgf-icgold)*6*n call info(3) if ( stepterm .ne. 1.d0 ) then accinf(niter,13)=1.d0 if ( restart ) then bfgsterm=-1.d0 return else call refresh(niter) if ( bfgsterm .eq. -3.d0 ) return restart=.true. goto 200 endif endif if ( fx0-fx .gt. epsf*(dabs(fx)+epsf)) then icsf=0 else icsf=icsf+1 endif if ( icsf .eq. 4 ) then bfgsterm=4.d0 return endif c********************************************* c compute update vectors z and y c********************************************* thet2=1.d0 s1=0.d0 s2=0.d0 s3=0.d0 den1=skal1(1,n,d,gradx0)/dscal h=1.d0/dsqrt(den1) do i=1,n y(i)=gradx0(i)*h dif1=difx(i) dif2=gradx(i)-gradx0(i) s1=dif1**2 + s1 s2=dif2**2 + s2 s3=dif1*dif2 + s3 z(i)=dif2 enddo if ( s3 .lt. eta*s2 ) then c********************************************** c bfgs-updating impossible from this information c f not uniformly convex of course c********************************************** s4=((sig*dscal)**2)*den1 c******** s(k)(transpose)*a(k)*s(k) in the usual notation if ( s3 .lt. .2d0*s4 ) thet2=.8d0*s4/(s4-s3) if ( thet2 .ne. 1.d0 ) imod=imod+1 s3=0.2d0*s4 h=(1.d0-thet2)*sig*dscal c******** powells modification c******** h*gradx0 = a(k)*s(k)*(1-thet2) do i=1,n z(i)=thet2*z(i)+h*gradx0(i) enddo endif accinf(niter,6)=thet2 s3=1.d0/dsqrt(s3) do i=1,n z(i)=z(i)*s3 enddo call info(4) accinf(niter,14)=1.d0 call cholsol(a,n,z,v) c**** v = b(k)*y(k)/sqrt(y's) in the usual notation s4=0.d0 s5=0.d0 do i=1,n difxin=difx(i)*s3 wi=difxin+v(i) c**** w=(sig(i)+(b*y)(i))/sqrt(y's) s5=s5+z(i)*wi s4=s4+v(i)*difxin enddo c******** update a now and compute its cholesky-decomposition if ( n .le. nsmall ) then c******** direct updating, save old information do i=1,n diaga0(i)=diaga(i) a0(i,i)=a(i,i) diaga(i)=diaga(i)-y(i)**2+z(i)**2 a(i,i)=diaga(i) c******** the strict lower part of a(a0) holds the quasi-newton c matrix, its diagonal being diaga (diaga0) c******** the upper triangular part holds the cholesky-factor do j=1,i-1 a0(i,j)=a(i,j) a0(j,i)=a(j,i) a(i,j)=a(i,j)-y(i)*y(j)+z(i)*z(j) enddo enddo call chol(a,n,a,ifail) if ( ifail .ne. 0 ) then c******** restore old matrices accinf(niter,14)=-1.d0 inoup=inoup+1 do i=1,n diaga(i)=diaga0(i) do j=1,n a(i,j)=a0(i,j) enddo enddo else iup=iup+1 endif else c********* large n: quasi-newton matrix given as its c cholesky-factor only. store/modify/restore c********* this information only do i=1,n diaga(i)=a(i,i) do j=1,i-1 a(i,j)=a(j,i) enddo enddo call aufdat(a,z,y,n,ifail) if ( ifail .ne. 0 ) then inoup=inoup+1 accinf(niter,14)=-1.d0 do i=1,n a(i,i)=diaga(i) do j=1,i-1 a(j,i)=a(i,j) enddo enddo else iup=iup+1 endif endif c**************************************************** c next step c**************************************************** if ( accinf (niter,14) .eq. 1.d0 ) then c**** recursive computation of trace(a) and trace(b) c**** it is stopped, if one trace becomes negative, which indicates c**** unreliability if ( tracea .gt. 0.d0 ) f tracea=tracea+vecnor(1,n,z)**2-vecnor(1,n,y)**2 if ( traceb .gt. 0.d0 ) f traceb=traceb-2.d0*s4+s1*s3**2*s5 call info(8) endif call info(5) goto 100 end c******************************************************************* c subprograms for computation of an updated cholesky- c decomposition, method of stewart c******************************************************************* c copyright = spellucci subroutine leftelim(a,b,y,yl,n) implicit none integer n,nx,i,j parameter (nx=300) double precision a(nx,nx),b(nx),y(nx),yl,h c leftelim assumes that the cholesky-factor of a c a=r(transpose)*r is stored in the upper half of a. c b is a right hand side. leftelim solves c y(transpose)*r = b(transpose) c yl=norm(y)**2 yl=0.d0 do i=1,n h=b(i) do j=1,i-1 h = h - a(j,i)*y(j) enddo h=h/a(i,i) y(i)=h yl = h**2 + yl enddo return end c********************************************************** double precision function dsq1(a,b) implicit none double precision a,b,a1,b1 a1=dabs(a) b1=dabs(b) if ( a1 .gt. b1 ) then dsq1=a1*dsqrt(1.d0+(b1/a1)**2) else if ( b1 .gt. a1 ) then dsq1=b1*dsqrt(1.d0+(a1/b1)**2) else dsq1=a1*dsqrt(2.d0) endif endif return end c******************************************************************* subroutine aufdat(r,z,y,n,ifail) implicit none integer n,ifail,nx,i,j,i1 c copyright = spellucci parameter (nx=300) double precision r(nx,nx),z(nx),y(nx) double precision sdiag(nx),rn1(nx),w(nx) double precision yl,zl,wl,wn1,ai,bi,h,dsq1 c aufdat computes the upper triangular cholesky-factor c r1 of c r(transpose)*r+z*z(transpose)-y*y(transpose) c and restores it in r. the strict lower triangle of r re- c mains unchanged. c ifail=1 if the decomposition does'nt exist, ifail=2 on dimension c error, ifail=0 on success. if ( n .gt. nx ) then ifail=2 return endif ifail=0 c save subdiagonal do i=1,n-1 sdiag(i)=r(i+1,i) r(i+1,i)=0.d0 enddo c step one: include z*z(transpose) zl=0.d0 do i=1,n zl = zl + z(i)**2 enddo if ( zl .ne. 0.d0 ) then c solve w(transpose)*r=z(transpose) call leftelim(r,z,w,wl,n) wl=dsqrt(wl+1.d0) c u(2)*u(3)*...*u(n)*w = ( norm(w),0,..,0)(transpose) c u(i) rotations do i=n,2,-1 if ( w(i) .ne. 0.d0 ) then i1=i-1 ai=w(i1) bi=w(i) w(i1)=dsq1(ai,bi) ai=ai/w(i1) bi=-bi/w(i1) r(i,i1)=bi*r(i1,i1) r(i1,i1)=ai*r(i1,i1) do j=i,n h = ai*r(i1,j) - bi*r(i,j) r(i,j) = bi*r(i1,j) + ai*r(i,j) r(i1,j) = h enddo endif enddo c r=d*r, d=diag(wl,1,...,1), r now hessenberg do i=1,n r(1,i)=r(1,i)*wl enddo c r=u(n-1)*...*u(1)*r now upper triangular, c u(i) givens-rotations do i=1,n-1 i1=i+1 ai=r(i,i) bi=-r(i1,i) h=dsq1(ai,bi) if ( h .ne. 0.d0 ) then ai=ai/h bi=bi/h r(i,i)=h r(i1,i)=0.d0 do j=i+1,n h = ai*r(i,j) - bi*r(i1,j) r(i1,j) = bi*r(i,j) + ai*r(i1,j) r(i,j) = h enddo endif enddo endif c step two : include -y*y(transpose) yl=0.d0 do i=1,n yl = yl + y(i)**2 enddo if ( yl .ne. 0.d0 ) then call leftelim(r,y,w,wl,n) if ( wl .ge. 1.d0 ) then ifail=1 else wl=dsqrt(1.d0-wl) wn1=wl c****************************************************** c ( r(new) ,0 ) ( r , w ) c (-----------) = u(1)*...u(n)*(-----------) c (y(transp),1) ((0,..,0),wl) c****************************************************** do i=n,1,-1 ai=wn1 bi=w(i) wn1=dsq1(ai,bi) ai=ai/wn1 bi=bi/wn1 rn1(i)=bi*r(i,i) r(i,i)=ai*r(i,i) do j=i+1,n h = ai*r(i,j) - bi*rn1(j) rn1(j) = bi*r(i,j) + ai*rn1(j) r(i,j) = h enddo enddo endif endif c restore subdiagonal do i=1,n-1 r(i+1,i)=sdiag(i) enddo return end c******************************************************************* c cholesky decomposition of a symmetric positive matrix a c the strict lower triangle remains unaffected c the upper triangle including the diagonal holds the cholesky- c factor. initially the strict upper triangle may be undefined c******************************************************************* subroutine chol(a,n,r,ifail) implicit none integer nx,ifail,n,i,j,k c copyright = spellucci parameter (nx=300) double precision a(nx,nx),r(nx,nx),h,s c computes the cholesky-factor r of a=r(transp)*r c and stores it in the upper triangle of r c a and r may be identical. the strict lower triangle c of a remains unchanged anyway. c ifail .ne. 0 if the decomposition does'nt exist, c =0 otherwise ifail=0 do i=1,n h=a(i,i) do j=1,i-1 h=h-r(j,i)**2 enddo if ( h .le. 0.d0 ) then ifail=1 return endif h=dsqrt(h) r(i,i)=h h=1.d0/h do k=i+1,n s=0.d0 do j=1,i-1 s=s+r(j,i)*r(j,k) enddo s=(a(k,i)-s)*h r(i,k)=s enddo enddo return end subroutine cholsol(a,n,b,x) implicit none integer nx,n,i,j c copyright = spellucci parameter (nx=300) double precision a(nx,nx),b(nx),c(nx),x(nx) double precision s c solves the linear system a*x=b, a symmetric positive definite c it is assumed that the cholesky-factor r of c a = r(transpose)*r c is stored in the upper triangle of a by chol. do i=1,n s=0.d0 do j=1,i-1 s=s+a(j,i)*c(j) enddo c(i)=(b(i)-s)/a(i,i) enddo do i=n,1,-1 s=0.d0 do j=i+1,n s=s+a(i,j)*x(j) enddo x(i)=(c(i)-s)/a(i,i) enddo return end c******************************************************************* c setup1 is used for initializing the quasi newton-update c******************************************************************* subroutine setup1 c copyright = spellucci implicit none integer nx,nsmall,maxit parameter (nx=300,nsmall=34,maxit=4000) double precision x,d,gradx,fx,xnorm,dnorm,sig,x0,d0, 1 gradx0,fx0,x0norm,d0norm,sig0,difx,ftest,dirder common/xdat/x(nx),d(nx),gradx(nx),fx,xnorm,dnorm,sig, 1 x0(nx),d0(nx),gradx0(nx),fx0,x0norm,d0norm,sig0, 2 difx(nx),ftest,dirder double precision stepterm,sigsm,sigla,alpha,beta,theta,dscal, 1 delta,kappa,tl common/stepparam/stepterm,sigsm,sigla,alpha,beta,delta,theta, 1 dscal,kappa,tl double precision epsx,epsg,epsf logical test1,test2,test3,test4 common/param/epsx,epsg,epsf,test1,test2,test3,test4 double precision epsmac,tolmac,deldif common/mach/epsmac,tolmac,deldif integer n common/dim/n integer icf,icgf,iup,inoup,ires,imod common/count/icf,icgf,iup,inoup,ires,imod double precision a,diaga,a0,diaga0,accinf,z,y,tracea,traceb common/bfgsdat/a(nx,nx),diaga(nx),a0(nsmall,nsmall), 1 diaga0(nsmall),accinf(maxit,14),z(nx),y(nx) 2 ,tracea,traceb logical analyt common/gradinf/analyt integer niter,itermax double precision eta,thet2,bfgsterm common/bfgsparam/niter,itermax,eta,thet2,bfgsterm c for testing purposes only this kind of fdat integer icase double precision xst common/fdat/xst(nx),icase character*80 ident character*30 xifile,xofile common/ident/ident,xifile,xofile double precision xd(nx),xst1(nx),gradx1(nx), f gradxs(nx),vecnor,gn logical err integer i,j c******* initialize with identity for testing purposes external vecnor save do i=1,n c*** xst1 should be in the range of definition of gradf xst1(i)=x(i)*(1.d0+epsx)+dsign(epsx**2,x(i)) xd(i)=x(i)-xst1(i) enddo err=.false. call gradf(x,gradxs,err) if ( err ) then bfgsterm=-3.d0 return endif gn=vecnor(1,n,gradxs) if ( accinf(1,1) .eq. 0.d0 ) then epsg=epsg*(gn+1.d0) write(60,fmt='('' epsg set to '', d11.4)') epsg endif call gradf(xst1,gradx1,err) if ( err ) then bfgsterm=-3.d0 return endif do i=1,n gradxs(i)=gradxs(i)-gradx1(i) enddo gn=vecnor(1,n,gradxs) c denominator n**2 is restart type I ccold gn=gn/vecnor(1,n,xd)/dble(n*n) gn=gn/vecnor(1,n,xd) if ( gn .le. dsqrt(epsmac) ) gn = 1.d0 c**** gn should be an estimate of the order of magnitude of the hessian c**** in case of a very small gn the computation of gn itself is expected c**** to be unreliable if ( n .le. nsmall ) then do i=1,n do j=1,n a(i,j)=0.d0 a0(i,j)=0.d0 enddo a(i,i)=dsqrt(gn) diaga(i)=gn diaga0(i)=gn enddo else do i=1,n do j=1,n a(j,i)=0.d0 enddo a(i,i)=dsqrt(gn) diaga(i)=gn enddo endif tracea=dble(n)*gn traceb=dble(n)/gn return end c******************************************************************* c refresh is used on restart in bfgs c restart is done only on occasion of a too small stepsize c sig .lt. sigsm c******************************************************************* subroutine refresh(niter) implicit none integer niter,nx,nsmall,maxit parameter (nx=300,nsmall=34,maxit=4000) double precision x,d,gradx,fx,xnorm,dnorm,sig,x0,d0, 1 gradx0,fx0,x0norm,d0norm,sig0,difx,ftest,dirder common/xdat/x(nx),d(nx),gradx(nx),fx,xnorm,dnorm,sig, 1 x0(nx),d0(nx),gradx0(nx),fx0,x0norm,d0norm,sig0, 2 difx(nx),ftest,dirder double precision stepterm,sigsm,sigla,alpha,beta,theta,dscal, 1 delta,kappa,tl common/stepparam/stepterm,sigsm,sigla,alpha,beta,delta,theta, 1 dscal,kappa,tl double precision epsx,epsg,epsf logical test1,test2,test3,test4 common/param/epsx,epsg,epsf,test1,test2,test3,test4 double precision epsmac,tolmac,deldif common/mach/epsmac,tolmac,deldif integer n common/dim/n integer icf,icgf,iup,inoup,ires,imod common/count/icf,icgf,iup,inoup,ires,imod double precision a,diaga,a0,diaga0,accinf,z,y,tracea,traceb common/bfgsdat/a(nx,nx),diaga(nx),a0(nsmall,nsmall), 1 diaga0(nsmall),accinf(maxit,14),z(nx),y(nx) 2 ,tracea,traceb logical analyt common/gradinf/analyt integer niterg,itermax double precision eta,thet2,bfgsterm common/bfgsparam/niterg,itermax,eta,thet2,bfgsterm c for testing purposes only this kind of fdat integer icase double precision xst common/fdat/xst(nx),icase character*80 ident character*30 xifile,xofile common/ident/ident,xifile,xofile save c copyright = spellucci ires=ires+1 accinf(niter,8)=1.d0 call setup1 c******** restart with scaled unit matrix return end c******************************************************************* c testx checks for sufficient accuracy c******************************************************************* subroutine testx implicit none integer nx,nsmall,maxit,i parameter (nx=300,nsmall=34,maxit=4000) double precision x,d,gradx,fx,xnorm,dnorm,sig,x0,d0, 1 gradx0,fx0,x0norm,d0norm,sig0,difx,ftest,dirder common/xdat/x(nx),d(nx),gradx(nx),fx,xnorm,dnorm,sig, 1 x0(nx),d0(nx),gradx0(nx),fx0,x0norm,d0norm,sig0, 2 difx(nx),ftest,dirder double precision stepterm,sigsm,sigla,alpha,beta,theta,dscal, 1 delta,kappa,tl common/stepparam/stepterm,sigsm,sigla,alpha,beta,delta,theta, 1 dscal,kappa,tl double precision epsx,epsg,epsf logical test1,test2,test3,test4 common/param/epsx,epsg,epsf,test1,test2,test3,test4 double precision epsmac,tolmac,deldif common/mach/epsmac,tolmac,deldif integer n common/dim/n integer icf,icgf,iup,inoup,ires,imod common/count/icf,icgf,iup,inoup,ires,imod double precision a,diaga,a0,diaga0,accinf,z,y,tracea,traceb common/bfgsdat/a(nx,nx),diaga(nx),a0(nsmall,nsmall), 1 diaga0(nsmall),accinf(maxit,14),z(nx),y(nx) 2 ,tracea,traceb logical analyt common/gradinf/analyt integer niter,itermax double precision eta,thet2,bfgsterm common/bfgsparam/niter,itermax,eta,thet2,bfgsterm c for testing purposes only this kind of fdat integer icase double precision xst common/fdat/xst(nx),icase character*80 ident character*30 xifile,xofile common/ident/ident,xifile,xofile logical xsatisfy,gsatisfy double precision an,sx,sg,damin,vecnor external vecnor save an=0.d0 sx=vecnor(1,n,difx) sg=vecnor(1,n,gradx) damin=a(1,1) do i=1,n an=an+diaga(i) damin=dmin1(damin,a(i,i)) enddo accinf(niter,5)=an c*** accinf(niter,5) norm of quasi newton update c*** accinf(niter,8) lower bound for condition number if ( tracea .gt. 0.d0 .and. traceb .gt. 0.d0 ) then accinf(niter,8)=tracea*traceb/dble(float(n*n)) else accinf(niter,8)=an/damin**2 endif xsatisfy = sx .le. epsx*(xnorm+epsx) gsatisfy = sg .le. dmax1(epsg,epsmac*accinf(niter,8)) if ( xsatisfy .and. gsatisfy ) bfgsterm=1.d0 return end c******************************************************************* c info provides for intermediate output for c testing purposes c info may have an empty body of course c******************************************************************* subroutine info(itest) c copyright = spellucci implicit double precision (a-h,o-z) parameter (nx=300,nsmall=34,maxit=4000) common/xdat/x(nx),d(nx),gradx(nx),fx,xnorm,dnorm,sig, 1 x0(nx),d0(nx),gradx0(nx),fx0,x0norm,d0norm,sig0, 2 difx(nx),ftest,dirder double precision delta,kappa,tl common/stepparam/stepterm,sigsm,sigla,alpha,beta,delta,theta, 1 dscal,kappa,tl common/param/epsx,epsg,epsf,test1,test2,test3,test4 common/mach/epsmac,tolmac,deldif common/dim/n common/count/icf,icgf,iup,inoup,ires,imod common/bfgsdat/a(nx,nx),diaga(nx),a0(nsmall,nsmall), 1 diaga0(nsmall),accinf(maxit,14),z(nx),y(nx) 2 ,tracea,traceb common/bfgsparam/niter,itermax,eta,thet2,bfgsterm logical test1,test2,test3,test4 double precision gam,fgam,ggam,abb common/trans/gam,fgam,ggam,abb integer ih1,ih2 character*40 heading save if ( itest .eq. 1 .and. test4 ) then write(50,fmt='(i4,2x,2(f11.7,1x),3(d10.3,1x))') f niter,x(1),x(2),accinf(niter,4),tracea,traceb endif if ( .not. test2 ) return goto ( 1000,2000,3000,4000,5000,6000,7000,8000,9000 ),itest c 1000 continue write(60,1100) niter,fx,accinf(niter,4) 1100 format(78('*'), 1 /1x,i3,'-th iteration step', 2 /1x,'f = ',d21.14,' norm(gradx)=',d21.14) write(60,fmt='('' x = '')') write(60,fmt='((1x,3(d21.14,2x)))') (x(j),j=1,n) if ( test3 .and. n .le. nsmall ) then heading='representation of qn-matrix' call matdr2(a,nx,nx,n,n,heading,60,.false.,1,8) write(60,fmt='('' diagonal of quasi-newton-matrix'')') write(60,1200) (diaga(j),j=1,n) 1200 format((1x,6(d11.4,2x))) endif return c 2000 continue write(60,2100) accinf(niter,5),accinf(niter,8) 2100 format(' norm(a)=',d11.4,' cond(a)=',d11.4) write(60,fmt='('' d ='')') write(60,fmt='((1x,6(d11.4,2x)))') (d(j),j=1,n) return c 3000 continue ih1=accinf(niter,12) write(60,3100) (accinf(niter,j),j=9,11),ih1 3100 format(' sigma=',d11.4,' dscal=',d11.4,' norm(d)=',d11.4, 1 ' calls (f) ',i3) return c 4000 continue write(60,4100) thet2 write(60,fmt='('' update vector z '')') write(60,4200) (z(i),i=1,n) write(60,fmt='('' update vector y '')') write(60,4200) (y(i),i=1,n) return 4100 format(' theta2=',d11.4) 4200 format((1x,6(d11.4,2x))) 5000 continue ih1=accinf(niter,13) ih2=accinf(niter,14) write(60,5100) ih1,ih2 5100 format(1x,'restart=',i2,' update=',i2) return 6000 continue write(60,fmt='(1x,''step*** gam ='',d11.4,'' fgam='',d11.4)') f gam, fgam return 7000 continue write(60,fmt='(1x,''step*** gam ='',d11.4,'' ggam='',d11.4)') f gam,ggam return 8000 continue write(60,fmt='(1x,'' tr(a)='',d11.4,'' tr(b)='',d11.4)') f tracea,traceb return 9000 continue write(60,fmt='(1x,'' stptrm='',d11.4)') abb end c***************************************************************** c subprogram for structured output of a matrix c uses a run time computed format string c***************************************************************** subroutine matdr2(a,me,ne,ma,na,heading,lognum,fix,nd1,nd2) c*** prints the matrix a on unit lognum using f-format if fix=.true. c and d-format otherwise with nd2 digits behind the point c if fix=.false. nd1 is not used and treated as nd1=0 c me and ne are the dimensions of the actual parameter a as defined c in the calling program unit and ma , na the dimensions actually used c ma<=me and na<=ne of course c output is for 80-column-devices. otherwise the parameter pwid below c may be changed implicit none integer me,ne,ma,na,nd1,nd2,lognum double precision a(me,ne) logical fix character*(*) heading c**** local integer anz,i,j,ju,jo,width,pwid,spacl,spacr parameter (pwid=80) character*2 widtxt,digtxt,anztxt,spaclt,spacrt character*30 form1 character*24 form2 if ( ma .gt. me .or. na .gt. ne .or. nd1 .lt. 0 .or. * nd2 .le. 0 ) then write(lognum,*) 'heading' write(lognum,*) 'erroneous call of matdr2' return endif if ( ma .gt. 999 .or. na .gt. 999 ) then write(lognum,*) heading write(lognum,*) 'called matdru with dim too large' return endif if ( fix ) then width=nd1+nd2+3 else width=nd2+7 endif write(lognum,fmt='(/1x,a)') heading anz=(pwid-11)/(width+1) write(anztxt,fmt='(i2)') anz write(digtxt,fmt='(i2)') nd2 write(widtxt,fmt='(i2)') width spacl=(width+1-3)/2 spacr=width+1-3-spacl write(spaclt,fmt='(i2)') spacl write(spacrt,fmt='(i2)') spacr form1='('' row/column'','//anztxt//'('//spaclt//'x,i3,' * //spacrt//'x))' if ( fix ) then form2='(4x,i3,4x,'//anztxt//'(f'//widtxt * //'.'//digtxt//',1x))' else form2='(4x,i3,4x,'//anztxt//'(d'//widtxt * //'.'//digtxt//',1x))' endif jo=0 do while ( jo .lt. na ) ju=jo+1 jo=min(ju+anz-1,na) write(lognum,fmt=form1) (j,j=ju,jo) do i=1,ma write(lognum,fmt=form2) i,(a(i,j),j=ju,jo) enddo enddo return end c******************************************************************* double precision function vecnor(i,j,x) implicit none c*** euclidean norm of x , avoid overflow integer i,j,k double precision x(*) double precision xm,h if ( i .gt. j ) then vecnor=0.d0 return endif xm=abs(x(i)) do k=i+1,j xm=max(xm,abs(x(k))) enddo if ( xm .eq. 0.d0 ) then vecnor=0.d0 return else h=0.d0 do k=i,j h=h+(x(k)/xm)**2 enddo vecnor=xm*sqrt(h) return endif end c*************************************************************** subroutine quadl(x,fx,f1x,y,fy,z) implicit none include 'bfgs_cons.inc' double precision x,fx,f1x,y,fy,z,u,v,t,va c****** quadratic interpolation of data (x,f(x),f'(x)), (y,f(y)) c****** and determination of zero z of p' c****** here z (should be) a minimum point. z not in (x,y) possible c****** in case of numerical uncertainty we let z=y c****** u=x-y v=fx-fy va=abs(v) t=va+abs(u) if ( t .eq. va ) then z=y return endif u=v/u-f1x v=-f1x*(y-x) va=abs(v) t=va+abs(u) if ( t .eq. va ) then z=y return endif z=x+v/u*p5 return end c********************************************************** subroutine quadr(x,fx,y,fy,f1y,z) implicit none include 'bfgs_cons.inc' double precision x,fx,f1y,y,fy,z,u,v,t,va c****** quadratic interpolation of data (x,f(x)), (y,f(y),f'(y)) c****** determination of the zero z of p' c****** z (should be) a minimum point here . z not in (x,y) is possible c****** in case of numerical uncertainty we let z=x c****** u=x-y v=fx-fy va=abs(v) t=va+abs(u) if ( t .eq. va ) then z=x return endif u=v/u-f1y v=-f1y*(x-y) va=abs(v) t=va+abs(u) if ( t .eq. va ) then z=x return endif z=y+v/u*p5 return end c********************************************************************** subroutine quadlr(x,fx,y,fy,z,fz,r) implicit none include 'bfgs_cons.inc' double precision x,fx,y,fy,z,fz,u,v,t,va,r,c1,c2 c****** quadratic interpolation of data (x,f(x)), (y,f(y)) c****** (z,f(z)) with x