CCC--- FIXME: Instead of subroutines forslv and baksol CCC--- ===== Use ./dtrsl.f !! [see dbksl() in ./S_compat.c !] CCC ~~~~~~~~~ CCC subroutines mvmlt[lsu] and sclmul should be REPLACED by BLAS ones! CCC CCC--- choldc(nr,n,a,diagmx,tol,addmax) is ``choleski + tolerance'' CCC ------ CCC it should make use of BLAS routines as [linkpack's dpofa!] c subroutine fdhess c c this subroutine calculates a numerical approximation to the upper c triangular portion of the second derivative matrix (the hessian). c algorithm a5.6.2 from dennis and schnable (1983), numerical methods c for unconstrained optimization and nonlinear equations, c prentice-hall, 321-322. c c programmed by richard h. jones, january 11, 1989 c c input to subroutine c c n.....the number of parameters c x.....vector of parameter values c fval..double precision value of function at x c fun...a function provided by the user which must be declared as c external in the calling program. its call must be of the c call fun(n,x,fval) where fval is the computed value of the c function c nfd...first dimension of h in the calling program c c output from subroutine c c h.....an n by n matrix of the approximate hessian c c work space c c step....a real array of length n c f.......a double precision array of length n c subroutine fdhess(n,x,fval,fun,h,nfd,step,f,ndigit,typx) implicit none external fun integer n,nfd,ndigit double precision x(n),fval, h(nfd,nfd),step(n),f(n),typx(n) c double precision eta,tempi,tempj,fii,fij integer i,j eta=(10.0d0**(-ndigit))**(1.0d0/3.0d0) do 10 i=1,n step(i)=eta*dmax1(x(i),typx(i)) if(typx(i).lt.0.0d0)step(i)=-step(i) tempi=x(i) x(i)=x(i)+step(i) step(i)=x(i)-tempi call fun(n,x,f(i)) x(i)=tempi 10 continue do 30 i=1,n tempi=x(i) x(i)=x(i)+2.0d0*step(i) call fun(n,x,fii) h(i,i)=((fval-f(i))+(fii-f(i)))/(step(i)*step(i)) x(i)=tempi+step(i) if(i.lt.n)then do 20 j=i+1,n tempj=x(j) x(j)=x(j)+step(j) call fun(n,x,fij) h(i,j)=((fval-f(i))+(fij-f(j)))/(step(i)*step(j)) x(j)=tempj 20 continue endif x(i)=tempi 30 continue return end c subroutine baksol c c purpose c c solve ax=b where a is upper triangular matrix. c note that a is input as a lower triangular matrix and c that this routine takes its transpose implicitly. c c parameters c c nr --> row dimension of matrix c n --> dimension of problem c a(n,n) --> lower triangular matrix (preserved) c x(n) <-- solution vector c b(n) --> right-hand side vector c c note c c if b is no longer required by calling routine, c then vectors b and x may share the same storage. c cstatic subroutine baksol(nr,n,a,x,b) implicit double precision (a-h,o-z) dimension a(nr,*),x(n),b(n) c c solve (l-transpose)x=b. (back solve) c i=n x(i)=b(i)/a(i,i) if(n.eq.1) return 30 ip1=i i=i-1 sum=0.0d0 do 40 j=ip1,n sum=sum+a(j,i)*x(j) 40 continue x(i)=(b(i)-sum)/a(i,i) if(i.gt.1) go to 30 return end c subroutine chlhsn c c purpose c c find the l(l-transpose) [written ll+] decomposition of the perturbed c model hessian matrix a+mu*i(where mu\0 and i is the identity matrix) c which is safely positive definite. if a is safely positive definite c upon entry, then mu=0. c c parameters c c nr --> row dimension of matrix c n --> dimension of problem c a(n,n) <--> on entry; "a" is model hessian (only lower c triangular part and diagonal stored) c on exit: a contains l of ll+ decomposition of c perturbed model hessian in lower triangular c part and diagonal and contains hessian in upper c triangular part and udiag c epsm --> machine epsilon c sx(n) --> diagonal scaling matrix for x c udiag(n) <-- on exit: contains diagonal of hessian c c internal variables c c tol tolerance c diagmn minimum element on diagonal of a c diagmx maximum element on diagonal of a c offmax maximum off-diagonal element of a c offrow sum of off-diagonal elements in a row of a c evmin minimum eigenvalue of a c evmax maximum eigenvalue of a c c description c c 1. if "a" has any negative diagonal elements, then choose mu>0 c such that the diagonal of a:=a+mu*i is all positive c with the ratio of its smallest to largest element on the c order of sqrt(epsm). c c 2. "a" undergoes a perturbed cholesky decomposition which c results in an ll+ decomposition of a+d, where d is a c non-negative diagonal matrix which is implicitly added to c "a" during the decomposition if "a" is not positive definite. c "a" is retained and not changed during this process by c copying l into the upper triangular part of "a" and the c diagonal into udiag. then the cholesky decomposition routine c is called. on return, addmax contains maximum element of d. c c 3. if addmax=0, "a" was positive definite going into step 2 c and return is made to calling program. otherwise, c the minimum number sdd which must be added to the c diagonal of a to make it safely strictly diagonally dominant c is calculated. since a+addmax*i and a+sdd*i are safely c positive definite, choose mu=min(addmax,sdd) and decompose c a+mu*i to obtain l. c cstatic subroutine chlhsn(nr,n,a,epsm,sx,udiag) implicit double precision (a-h,o-z) dimension a(nr,*),sx(n),udiag(n) c c scale hessian c pre- and post- multiply "a" by inv(sx) c do 20 j=1,n do 10 i=j,n a(i,j)=a(i,j)/(sx(i)*sx(j)) 10 continue 20 continue c c step1 c ----- c note: if a different tolerance is desired throughout this c algorithm, change tolerance here: c tol=sqrt(epsm) c diagmx=a(1,1) diagmn=a(1,1) if(n.eq.1) go to 35 do 30 i=2,n if(a(i,i).lt.diagmn) diagmn=a(i,i) if(a(i,i).gt.diagmx) diagmx=a(i,i) 30 continue 35 posmax=max(diagmx,0.0d0) c c diagmn .le. 0 c if(diagmn.gt.posmax*tol) go to 100 c if(diagmn.le.posmax*tol) c then amu=tol*(posmax-diagmn)-diagmn if(amu.ne.0.0d0) go to 60 c if(amu.eq.0.0d0) c then c c find largest off-diagonal element of a offmax=0.0d0 if(n.eq.1) go to 50 do 45 i=2,n im1=i-1 do 40 j=1,im1 if(abs(a(i,j)).gt.offmax) offmax=abs(a(i,j)) 40 continue 45 continue 50 amu=offmax if(amu.ne.0.0d0) go to 55 c if(amu.eq.0.0d0) c then amu=1.0d0 go to 60 c else 55 amu=amu*(1.0d0+tol) c endif c endif c c a=a + mu*i c 60 do 65 i=1,n a(i,i)=a(i,i)+amu 65 continue diagmx=diagmx+amu c endif c c step2 c c copy lower triangular part of "a" to upper triangular part c and diagonal of "a" to udiag c 100 continue do 110 j=1,n udiag(j)=a(j,j) if(j.eq.n) go to 110 jp1=j+1 do 105 i=jp1,n a(j,i)=a(i,j) 105 continue 110 continue c call choldc(nr,n,a,diagmx,tol,addmax) c c c step3 c c if addmax=0, "a" was positive definite going into step 2, c the ll+ decomposition has been done, and we return. c otherwise, addmax>0. perturb "a" so that it is safely c diagonally dominant and find ll+ decomposition c if(addmax.le.0.0d0) go to 170 c if(addmax.gt.0.0d0) c then c c restore original "a" (lower triangular part and diagonal) c do 120 j=1,n a(j,j)=udiag(j) if(j.eq.n) go to 120 jp1=j+1 do 115 i=jp1,n a(i,j)=a(j,i) 115 continue 120 continue c c find sdd such that a+sdd*i is safely positive definite c note: evmin<0 since a is not positive definite; c evmin=0.0d0 evmax=a(1,1) do 150 i=1,n offrow=0.0d0 if(i.eq.1) go to 135 im1=i-1 do 130 j=1,im1 offrow=offrow+abs(a(i,j)) 130 continue 135 if(i.eq.n) go to 145 ip1=i+1 do 140 j=ip1,n offrow=offrow+abs(a(j,i)) 140 continue 145 evmin=min(evmin,a(i,i)-offrow) evmax=max(evmax,a(i,i)+offrow) 150 continue sdd=tol*(evmax-evmin)-evmin c c perturb "a" and decompose again c amu=min(sdd,addmax) do 160 i=1,n a(i,i)=a(i,i)+amu udiag(i)=a(i,i) 160 continue c c "a" now guaranteed safely positive definite c call choldc(nr,n,a,0.0d0,tol,addmax) c endif c c unscale hessian and cholesky decomposition matrix c 170 do 190 j=1,n do 175 i=j,n a(i,j)=sx(i)*a(i,j) 175 continue if(j.eq.1) go to 185 jm1=j-1 do 180 i=1,jm1 a(i,j)=sx(i)*sx(j)*a(i,j) 180 continue 185 udiag(j)=udiag(j)*sx(j)*sx(j) 190 continue return end c subroutine choldc c c purpose c c find the perturbed l(l-transpose) [written ll+] decomposition c of a+d, where d is a non-negative diagonal matrix added to a if c necessary to allow the cholesky decomposition to continue. c c parameters c c nr --> row dimension of matrix c n --> dimension of problem c a(n,n) <--> on entry: matrix for which to find perturbed c cholesky decomposition c on exit: contains l of ll+ decomposition c in lower triangular part and diagonal of "a" c diagmx --> maximum diagonal element of "a" c tol --> tolerance c addmax <-- maximum amount implicitly added to diagonal of "a" c in forming the cholesky decomposition of a+d c internal variables c c aminl smallest element allowed on diagonal of l c amnlsq =aminl**2 c offmax maximum off-diagonal element in column of a c c c description c c the normal cholesky decomposition is performed. however, if at any c point the algorithm would attempt to set l(i,i)=sqrt(temp) c with temp < tol*diagmx, then l(i,i) is set to sqrt(tol*diagmx) c instead. this is equivalent to adding tol*diagmx-temp to a(i,i) c c cstatic subroutine choldc(nr,n,a,diagmx,tol,addmax) implicit double precision (a-h,o-z) dimension a(nr,*) c addmax=0.0d0 aminl=sqrt(diagmx*tol) amnlsq=aminl*aminl c c form column j of l c do 100 j=1,n c c find diagonal elements of l c sum=0.0d0 if(j.eq.1) go to 20 jm1=j-1 do 10 k=1,jm1 sum=sum + a(j,k)*a(j,k) 10 continue 20 temp=a(j,j)-sum if(temp.lt.amnlsq) go to 30 c if(temp.ge.aminl**2) c then a(j,j)=sqrt(temp) go to 40 c else c c find maximum off-diagonal element in column c 30 offmax=0.0d0 if(j.eq.n) go to 37 jp1=j+1 do 35 i=jp1,n if(abs(a(i,j)).gt.offmax) offmax=abs(a(i,j)) 35 continue 37 if(offmax.le.amnlsq) offmax=amnlsq c c add to diagonal element to allow cholesky decomposition to continue c a(j,j)=sqrt(offmax) addmax=max(addmax,offmax-temp) c endif c c find i,j element of lower triangular matrix 40 if(j.eq.n) go to 100 jp1=j+1 do 70 i=jp1,n sum=0.0d0 if(j.eq.1) go to 60 jm1=j-1 do 50 k=1,jm1 sum=sum+a(i,k)*a(j,k) 50 continue 60 a(i,j)=(a(i,j)-sum)/a(j,j) 70 continue 100 continue return end c subroutine d1fcn c c purpose c c dummy routine to prevent unsatisfied external diagnostic c when specific analytic gradient function not supplied. c cstatic subroutine d1fcn(n,x,g) implicit double precision (a-h,o-z) dimension x(n),g(n) g(n)=g(n) x(n)=x(n) end c subroutine d2fcn c c purpose c c dummy routine to prevent unsatisfied external diagnostic c when specific analytic hessian function not supplied. c cstatic subroutine d2fcn(nr,n,x,h) implicit double precision (a-h,o-z) dimension x(n),h(nr,*) h(nr,1)=h(nr,1) x(n)=x(n) end c subroutine dfault c c purpose c c set default values for each input variable to c minimization algorithm. c c parameters c c n --> dimension of problem c x(n) --> initial guess to solution (to compute max step size) c typsiz(n) <-- typical size for each component of x c fscale <-- estimate of scale of minimization function c method <-- algorithm to use to solve minimization problem c iexp <-- =0 if minimization function not expensive to evaluate c msg <-- message to inhibit certain automatic checks + output c ndigit <-- number of good digits in minimization function c itnlim <-- maximum number of allowable iterations c iagflg <-- =0 if analytic gradient not supplied c iahflg <-- =0 if analytic hessian not supplied c ipr <-- device to which to send output c dlt <-- trust region radius c gradtl <-- tolerance at which gradient considered close enough c to zero to terminate algorithm c stepmx <-- value of zero to trip default maximum in optchk c steptl <-- tolerance at which successive iterates considered c close enough to terminate algorithm c cstatic subroutine dfault(n,x,typsiz,fscale,method,iexp,msg,ndigit, + itnlim,iagflg,iahflg,ipr,dlt,gradtl,stepmx,steptl) implicit double precision (a-h,o-z) dimension typsiz(n),x(n) x(n)=x(n) c c set typical size of x and minimization function do 10 i=1,n typsiz(i)=1.0d0 10 continue fscale=1.0d0 c c set tolerances c dlt=-1.0d0 epsm=d1mach(4) gradtl=epsm**(1.0d0/3.0d0) stepmx=0.0d0 steptl=sqrt(epsm) c c set flags c method=1 iexp=1 msg=0 ndigit=-1 itnlim=150 iagflg=0 iahflg=0 ipr=i1mach(2) c return end c subroutine dogdrv c c purpose c c find a next newton iterate (xpls) by the double dogleg method c c parameters c c nr --> row dimension of matrix c n --> dimension of problem c x(n) --> old iterate x[k-1] c f --> function value at old iterate, f(x) c g(n) --> gradient at old iterate, g(x), or approximate c a(n,n) --> cholesky decomposition of hessian c in lower triangular part and diagonal c p(n) --> newton step c xpls(n) <-- new iterate x[k] c fpls <-- function value at new iterate, f(xpls) c fcn --> name of subroutine to evaluate function c sx(n) --> diagonal scaling matrix for x c stepmx --> maximum allowable step size c steptl --> relative step size at which successive iterates c considered close enough to terminate algorithm c dlt <--> trust region radius c [retain value between successive calls] c iretcd <-- return code c =0 satisfactory xpls found c =1 failed to find satisfactory xpls sufficiently c distinct from x c mxtake <-- boolean flag indicating step of maximum length used c sc(n) --> workspace [current step] c wrk1(n) --> workspace (and place holding argument to tregup) c wrk2(n) --> workspace c wrk3(n) --> workspace c ipr --> device to which to send output c cstatic subroutine dogdrv(nr,n,x,f,g,a,p,xpls,fpls,fcn,sx,stepmx, + steptl,dlt,iretcd,mxtake,sc,wrk1,wrk2,wrk3,ipr,itncnt) implicit double precision (a-h,o-z) dimension x(n),xpls(n),g(n),p(n) dimension sx(n) dimension sc(n),wrk1(n),wrk2(n),wrk3(n) dimension a(nr,*) logical fstdog,nwtake,mxtake integer itncnt external fcn c iretcd=4 fstdog=.true. tmp=0.0d0 do 5 i=1,n tmp=tmp+sx(i)*sx(i)*p(i)*p(i) 5 continue rnwtln=sqrt(tmp) c$ write(ipr,954) rnwtln c 100 continue c c find new step by double dogleg algorithm c call dogstp(nr,n,g,a,p,sx,rnwtln,dlt,nwtake,fstdog, + wrk1,wrk2,cln,eta,sc,ipr,stepmx) c c check new point and update trust region c call tregup(nr,n,x,f,g,a,fcn,sc,sx,nwtake,stepmx,steptl,dlt, + iretcd,wrk3,fplsp,xpls,fpls,mxtake,ipr,2,wrk1) if(iretcd.le.1) return go to 100 c%950 format(42h dogdrv initial trust region not given., c% + 22h compute cauchy step.) c%951 format(18h dogdrv alpha =,e20.13/ c% + 18h dogdrv beta =,e20.13/ c% + 18h dogdrv dlt =,e20.13/ c% + 18h dogdrv nwtake=,l1 ) c%952 format(28h dogdrv current step (sc)) c%954 format(18h0dogdrv rnwtln=,e20.13) c%955 format(14h dogdrv ,5(e20.13,3x)) end c subroutine dogstp c c purpose c c find new step by double dogleg algorithm c c parameters c c nr --> row dimension of matrix c n --> dimension of problem c g(n) --> gradient at current iterate, g(x) c a(n,n) --> cholesky decomposition of hessian in c lower part and diagonal c p(n) --> newton step c sx(n) --> diagonal scaling matrix for x c rnwtln --> newton step length c dlt <--> trust region radius c nwtake <--> boolean, =.true. if newton step taken c fstdog <--> boolean, =.true. if on first leg of dogleg c ssd(n) <--> workspace [cauchy step to the minimum of the c quadratic model in the scaled steepest descent c direction] [retain value between successive calls] c v(n) <--> workspace [retain value between successive calls] c cln <--> cauchy length c [retain value between successive calls] c eta [retain value between successive calls] c sc(n) <-- current step c ipr --> device to which to send output c stepmx --> maximum allowable step size c c internal variables c c cln length of cauchy step c cstatic subroutine dogstp(nr,n,g,a,p,sx,rnwtln,dlt,nwtake,fstdog, + ssd,v,cln,eta,sc,ipr,stepmx) implicit double precision (a-h,o-z) dimension g(n),p(n) dimension sx(n) dimension sc(n),ssd(n),v(n) dimension a(nr,*) logical nwtake,fstdog ipr=ipr c c can we take newton step c if(rnwtln.gt.dlt) go to 100 c if(rnwtln.le.dlt) c then nwtake=.true. do 10 i=1,n sc(i)=p(i) 10 continue dlt=rnwtln c$ write(ipr,951) go to 700 c else c c newton step too long c cauchy step is on double dogleg curve c 100 nwtake=.false. if(.not.fstdog) go to 200 c if(fstdog) c then c c calculate double dogleg curve (ssd) fstdog=.false. alpha=0.0d0 do 110 i=1,n alpha=alpha + (g(i)*g(i))/(sx(i)*sx(i)) 110 continue beta=0.0d0 do 130 i=1,n tmp=0.0d0 do 120 j=i,n tmp=tmp + (a(j,i)*g(j))/(sx(j)*sx(j)) 120 continue beta=beta+tmp*tmp 130 continue do 140 i=1,n ssd(i)=-(alpha/beta)*g(i)/sx(i) 140 continue cln=alpha*sqrt(alpha)/beta eta=.2 + (.8*alpha*alpha)/(-beta*ddot(n,g,1,p,1)) do 150 i=1,n v(i)=eta*sx(i)*p(i) - ssd(i) 150 continue if (dlt .eq. (-1.0d0)) dlt = min(cln, stepmx) c$ write(ipr,954) alpha,beta,cln,eta c$ write(ipr,955) c$ write(ipr,960) (ssd(i),i=1,n) c$ write(ipr,956) c$ write(ipr,960) (v(i),i=1,n) c endif 200 if(eta*rnwtln.gt.dlt) go to 220 c if(eta*rnwtln .le. dlt) c then c c take partial step in newton direction c do 210 i=1,n sc(i)=(dlt/rnwtln)*p(i) 210 continue c$ write(ipr,957) go to 700 c else 220 if(cln.lt.dlt) go to 240 c if(cln.ge.dlt) c then c take step in steepest descent direction c do 230 i=1,n sc(i)=(dlt/cln)*ssd(i)/sx(i) 230 continue c$ write(ipr,958) go to 700 c else c c calculate convex combination of ssd and eta*p c which has scaled length dlt c 240 dot1=ddot(n,v,1,ssd,1) dot2=ddot(n,v,1,v,1) alam=(-dot1+sqrt((dot1*dot1)-dot2*(cln*cln-dlt*dlt)))/dot2 do 250 i=1,n sc(i)=(ssd(i) + alam*v(i))/sx(i) 250 continue c$ write(ipr,959) c endif c endif c endif 700 continue c$ write(ipr,952) fstdog,nwtake,rnwtln,dlt c$ write(ipr,953) c$ write(ipr,960) (sc(i),i=1,n) return c c%951 format(27h0dogstp take newton step) c%952 format(18h dogstp fstdog=,l1/ c% + 18h dogstp nwtake=,l1/ c% + 18h dogstp rnwtln=,e20.13/ c% + 18h dogstp dlt =,e20.13) c%953 format(28h dogstp current step (sc)) c%954 format(18h dogstp alpha =,e20.13/ c% + 18h dogstp beta =,e20.13/ c% + 18h dogstp cln =,e20.13/ c% + 18h dogstp eta =,e20.13) c%955 format(28h dogstp cauchy step (ssd)) c%956 format(12h dogstp v) c%957 format(48h0dogstp take partial step in newton direction) c%958 format(50h0dogstp take step in steepest descent direction) c%959 format(39h0dogstp take convex combination step) c%960 format(14h dogstp ,5(e20.13,3x)) end c subroutine forslv c c purpose c c solve ax=b where a is lower triangular matrix c c parameters c c nr --> row dimension of matrix c n --> dimension of problem c a(n,n) --> lower triangular matrix (preserved) c x(n) <-- solution vector c b(n) --> right-hand side vector c c note c c if b is no longer required by calling routine, c then vectors b and x may share the same storage. c cstatic subroutine forslv(nr,n,a,x,b) implicit double precision (a-h,o-z) dimension a(nr,*),x(n),b(n) c c solve lx=b. (foreward solve) c x(1)=b(1)/a(1,1) if(n.eq.1) return do 20 i=2,n sum=0.0d0 im1=i-1 do 10 j=1,im1 sum=sum+a(i,j)*x(j) 10 continue x(i)=(b(i)-sum)/a(i,i) 20 continue return end c subroutine fstocd c c purpose c c find central difference approximation g to the first derivative c (gradient) of the function defined by fcn at the point x. c c parameters c c n --> dimension of problem c x --> point at which gradient is to be approximated. c fcn --> name of subroutine to evaluate function. c sx --> diagonal scaling matrix for x. c rnoise --> relative noise in fcn [f(x)]. c g <-- central difference approximation to gradient. c c cstatic subroutine fstocd (n, x, fcn, sx, rnoise, g) implicit none integer n external fcn double precision x(n), sx(n), g(n), rnoise c integer i double precision third, stepi,xtempi, fplus, fminus c c find i th stepsize, evaluate two neighbors in direction of i th c unit vector, and evaluate i th component of gradient. c third = 1.0d0/3.0d0 do 10 i = 1, n stepi = rnoise**third * max(abs(x(i)), 1.0d0/sx(i)) xtempi = x(i) x(i) = xtempi + stepi call fcn (n, x, fplus) x(i) = xtempi - stepi call fcn (n, x, fminus) x(i) = xtempi g(i) = (fplus - fminus)/(2.0d0*stepi) 10 continue return end c subroutine fstofd c c purpose c c find first order forward finite difference approximation "a" to the c first derivative of the function defined by the subprogram "fname" c evaluated at the new iterate "xpls". c c for optimization use this routine to estimate: c 1) the first derivative (gradient) of the optimization function "fcn c analytic user routine has been supplied; c 2) the second derivative (hessian) of the optimization function c if no analytic user routine has been supplied for the hessian but c one has been supplied for the gradient ("fcn") and if the c optimization function is inexpensive to evaluate c c note c c _m=1 (optimization) algorithm estimates the gradient of the function c (fcn). fcn(x) # f: r(n)-->r(1) c _m=n (systems) algorithm estimates the jacobian of the function c fcn(x) # f: r(n)-->r(n). c _m=n (optimization) algorithm estimates the hessian of the optimizatio c function, where the hessian is the first derivative of "fcn" c c parameters c c nr --> row dimension of matrix c m --> number of rows in a c n --> number of columns in a; dimension of problem c xpls(n) --> new iterate: x[k] c fcn --> name of subroutine to evaluate function c fpls(m) --> _m=1 (optimization) function value at new iterate: c fcn(xpls) c _m=n (optimization) value of first derivative c (gradient) given by user function fcn c _m=n (systems) function value of associated c minimization function c a(nr,n) <-- finite difference approximation (see note). only c lower triangular matrix and diagonal are returned c sx(n) --> diagonal scaling matrix for x c rnoise --> relative noise in fcn [f(x)] c fhat(m) --> workspace c icase --> =1 optimization (gradient) c =2 systems c =3 optimization (hessian) c c internal variables c c stepsz - stepsize in the j-th variable direction c cstatic subroutine fstofd(nr,m,n,xpls,fcn,fpls,a,sx,rnoise,fhat,icase) implicit double precision (a-h,o-z) dimension xpls(n),fpls(m) dimension fhat(m) dimension sx(n) dimension a(nr,*) c c find j-th column of a c each column is derivative of f(fcn) with respect to xpls(j) c do 30 j=1,n stepsz=sqrt(rnoise)*max(abs(xpls(j)),1.0d0/sx(j)) xtmpj=xpls(j) xpls(j)=xtmpj+stepsz call fcn(n,xpls,fhat) xpls(j)=xtmpj do 20 i=1,m a(i,j)=(fhat(i)-fpls(i))/stepsz 20 continue 30 continue if(icase.ne.3) return c c if computing hessian, a must be symmetric c if(n.eq.1) return nm1=n-1 do 50 j=1,nm1 jp1=j+1 do 40 i=jp1,m a(i,j)=(a(i,j)+a(j,i))/2.0d0 40 continue 50 continue return end c subroutine grdchk c c purpose c c check analytic gradient against estimated gradient c c parameters c c n --> dimension of problem c x(n) --> estimate to a root of fcn c fcn --> name of subroutine to evaluate optimization function c must be declared external in calling routine c fcn: r(n) --> r(1) c f --> function value: fcn(x) c g(n) --> gradient: g(x) c typsiz(n) --> typical size for each component of x c sx(n) --> diagonal scaling matrix: sx(i)=1./typsiz(i) c fscale --> estimate of scale of objective function fcn c rnf --> relative noise in optimization function fcn c analtl --> tolerance for comparison of estimated and c analytical gradients c wrk1(n) --> workspace c msg <-- message or error code c on output: =-21, probable coding error of gradient c ipr --> device to which to send output c cstatic subroutine grdchk(n,x,fcn,f,g,typsiz,sx,fscale,rnf, + analtl,wrk1,msg,ipr) implicit double precision (a-h,o-z) dimension x(n),g(n) dimension sx(n),typsiz(n) dimension wrk1(n) external fcn c c compute first order finite difference gradient and compare to c analytic gradient. c call fstofd(1,1,n,x,fcn,f,wrk1,sx,rnf,wrk,1) ker=0 do 5 i=1,n gs=max(abs(f),fscale)/max(abs(x(i)),typsiz(i)) if(abs(g(i)-wrk1(i)).gt.max(abs(g(i)),gs)*analtl) ker=1 5 continue if(ker.eq.0) go to 20 c% write(ipr,901) c% write(ipr,902) (i,g(i),wrk1(i),i=1,n) msg=-21 20 continue return c%901 format(47h0grdchk probable error in coding of analytic, c% + 19h gradient function./ c% + 16h grdchk comp,12x,8hanalytic,12x,8hestimate) c%902 format(11h grdchk ,i5,3x,e20.13,3x,e20.13) end c subroutine heschk c c purpose c c check analytic hessian against estimated hessian c (this may be done only if the user supplied analytic hessian c d2fcn fills only the lower triangular part and diagonal of a) c c parameters c c nr --> row dimension of matrix c n --> dimension of problem c x(n) --> estimate to a root of fcn c fcn --> name of subroutine to evaluate optimization function c must be declared external in calling routine c fcn: r(n) --> r(1) c d1fcn --> name of subroutine to evaluate gradient of fcn. c must be declared external in calling routine c d2fcn --> name of subroutine to evaluate hessian of fcn. c must be declared external in calling routine c f --> function value: fcn(x) c g(n) <-- gradient: g(x) c a(n,n) <-- on exit: hessian in lower triangular part and diag c typsiz(n) --> typical size for each component of x c sx(n) --> diagonal scaling matrix: sx(i)=1./typsiz(i) c rnf --> relative noise in optimization function fcn c analtl --> tolerance for comparison of estimated and c analytical gradients c iagflg --> =1 if analytic gradient supplied c udiag(n) --> workspace c wrk1(n) --> workspace c wrk2(n) --> workspace c msg <--> message or error code c on input : if =1xx do not compare anal + est hess c on output: =-22, probable coding error of hessian c ipr --> device to which to send output c cstatic subroutine heschk(nr,n,x,fcn,d1fcn,d2fcn,f,g,a,typsiz,sx,rnf, + analtl,iagflg,udiag,wrk1,wrk2,msg,ipr) implicit double precision (a-h,o-z) dimension x(n),g(n),a(nr,*) dimension typsiz(n),sx(n) dimension udiag(n),wrk1(n),wrk2(n) external fcn,d1fcn c c compute finite difference approximation a to the hessian. c if(iagflg.eq.1) call fstofd(nr,n,n,x,d1fcn,g,a,sx,rnf,wrk1,3) if(iagflg.ne.1) call sndofd(nr,n,x,fcn,f,a,sx,rnf,wrk1,wrk2) c ker=0 c c copy lower triangular part of "a" to upper triangular part c and diagonal of "a" to udiag c do 30 j=1,n udiag(j)=a(j,j) if(j.eq.n) go to 30 jp1=j+1 do 25 i=jp1,n a(j,i)=a(i,j) 25 continue 30 continue c c compute analytic hessian and compare to finite difference c approximation. c call d2fcn(nr,n,x,a) do 40 j=1,n hs=max(abs(g(j)),1.0d0)/max(abs(x(j)),typsiz(j)) if(abs(a(j,j)-udiag(j)).gt.max(abs(udiag(j)),hs)*analtl) + ker=1 if(j.eq.n) go to 40 jp1=j+1 do 35 i=jp1,n if(abs(a(i,j)-a(j,i)).gt.max(abs(a(i,j)),hs)*analtl) ker=1 35 continue 40 continue c if(ker.eq.0) go to 90 c% write(ipr,901) c% do 50 i=1,n c% if(i.eq.1) go to 45 c% im1=i-1 c% do 43 j=1,im1 c% write(ipr,902) i,j,a(i,j),a(j,i) c% 43 continue c% 45 write(ipr,902) i,i,a(i,i),udiag(i) c% 50 continue msg=-22 c endif 90 continue return c%901 format(47h heschk probable error in coding of analytic, c% + 18h hessian function./ c% + 21h heschk row col,14x,8hanalytic,14x,10h(estimate)) c%902 format(11h heschk ,2i5,2x,e20.13,2x,1h(,e20.13,1h)) end c subroutine hookdr c c purpose c c find a next newton iterate (xpls) by the more-hebdon method c c parameters c c nr --> row dimension of matrix c n --> dimension of problem c x(n) --> old iterate x[k-1] c f --> function value at old iterate, f(x) c g(n) --> gradient at old iterate, g(x), or approximate c a(n,n) --> cholesky decomposition of hessian in lower c triangular part and diagonal. c hessian in upper triangular part and udiag. c udiag(n) --> diagonal of hessian in a(.,.) c p(n) --> newton step c xpls(n) <-- new iterate x[k] c fpls <-- function value at new iterate, f(xpls) c fcn --> name of subroutine to evaluate function c sx(n) --> diagonal scaling matrix for x c stepmx --> maximum allowable step size c steptl --> relative step size at which successive iterates c considered close enough to terminate algorithm c dlt <--> trust region radius c iretcd <-- return code c =0 satisfactory xpls found c =1 failed to find satisfactory xpls sufficiently c distinct from x c mxtake <-- boolean flag indicating step of maximum length used c amu <--> [retain value between successive calls] c dltp <--> [retain value between successive calls] c phi <--> [retain value between successive calls] c phip0 <--> [retain value between successive calls] c sc(n) --> workspace c xplsp(n) --> workspace c wrk0(n) --> workspace c epsm --> machine epsilon c itncnt --> iteration count c ipr --> device to which to send output c cstatic subroutine hookdr(nr,n,x,f,g,a,udiag,p,xpls,fpls,fcn,sx,stepmx, + steptl,dlt,iretcd,mxtake,amu,dltp,phi,phip0, + sc,xplsp,wrk0,epsm,itncnt,ipr) implicit double precision (a-h,o-z) dimension x(n),g(n),p(n),xpls(n),sx(n) dimension a(nr,*),udiag(n) dimension sc(n),xplsp(n),wrk0(n) logical mxtake,nwtake logical fstime external fcn c iretcd=4 fstime=.true. tmp=0.0d0 do 5 i=1,n tmp=tmp+sx(i)*sx(i)*p(i)*p(i) 5 continue rnwtln=sqrt(tmp) c$ write(ipr,954) rnwtln c if(itncnt.gt.1) go to 100 c if(itncnt.eq.1) c then amu=0.0d0 c c if first iteration and trust region not provided by user, c compute initial trust region. c if(dlt.ne. (-1.0d0)) go to 100 c if(dlt.eq. (-1.0d0)) c then alpha=0.0d0 do 10 i=1,n alpha=alpha+(g(i)*g(i))/(sx(i)*sx(i)) 10 continue beta=0.0d0 do 30 i=1,n tmp=0.0d0 do 20 j=i,n tmp=tmp + (a(j,i)*g(j))/(sx(j)*sx(j)) 20 continue beta=beta+tmp*tmp 30 continue dlt=alpha*sqrt(alpha)/beta dlt = min(dlt, stepmx) c$ write(ipr,950) c$ write(ipr,951) alpha,beta,dlt c endif c endif c 100 continue c c find new step by more-hebdon algorithm c call hookst(nr,n,g,a,udiag,p,sx,rnwtln,dlt,amu, + dltp,phi,phip0,fstime,sc,nwtake,wrk0,epsm,ipr) dltp=dlt c c check new point and update trust region c call tregup(nr,n,x,f,g,a,fcn,sc,sx,nwtake,stepmx,steptl, + dlt,iretcd,xplsp,fplsp,xpls,fpls,mxtake,ipr,3,udiag) if(iretcd.le.1) return go to 100 c c%950 format(43h hookdr initial trust region not given. , c% + 21h compute cauchy step.) c%951 format(18h hookdr alpha =,e20.13/ c% + 18h hookdr beta =,e20.13/ c% + 18h hookdr dlt =,e20.13) c%952 format(28h hookdr current step (sc)) c%954 format(18h0hookdr rnwtln=,e20.13) c%955 format(14h hookdr ,5(e20.13,3x)) end c subroutine hookst c c purpose c c find new step by more-hebdon algorithm c c parameters c c nr --> row dimension of matrix c n --> dimension of problem c g(n) --> gradient at current iterate, g(x) c a(n,n) --> cholesky decomposition of hessian in c lower triangular part and diagonal. c hessian or approx in upper triangular part c udiag(n) --> diagonal of hessian in a(.,.) c p(n) --> newton step c sx(n) --> diagonal scaling matrix for n c rnwtln --> newton step length c dlt <--> trust region radius c amu <--> [retain value between successive calls] c dltp --> trust region radius at last exit from this routine c phi <--> [retain value between successive calls] c phip0 <--> [retain value between successive calls] c fstime <--> boolean. =.true. if first entry to this routine c during k-th iteration c sc(n) <-- current step c nwtake <-- boolean, =.true. if newton step taken c wrk0 --> workspace c epsm --> machine epsilon c ipr --> device to which to send output c cstatic subroutine hookst(nr,n,g,a,udiag,p,sx,rnwtln,dlt,amu, + dltp,phi,phip0,fstime,sc,nwtake,wrk0,epsm,ipr) implicit double precision (a-h,o-z) dimension g(n),p(n),sx(n),sc(n),wrk0(n) dimension a(nr,*),udiag(n) logical nwtake,done logical fstime c c hi and alo are constants used in this routine. c change here if other values are to be substituted. ipr=ipr hi=1.5d0 alo=0.75d0 c if(rnwtln.gt.hi*dlt) go to 15 c if(rnwtln.le.hi*dlt) c then c c take newton step c nwtake=.true. do 10 i=1,n sc(i)=p(i) 10 continue dlt=min(dlt,rnwtln) amu=0.0d0 c$ write(ipr,951) return c else c c newton step not taken c 15 continue c$ write(ipr,952) nwtake=.false. if(amu.le.0.0d0) go to 20 c if(amu.gt.0.0d0) c then amu=amu- (phi+dltp) *((dltp-dlt)+phi)/(dlt*phip) c$ write(ipr,956) amu c endif 20 continue phi=rnwtln-dlt if(.not.fstime) go to 28 c if(fstime) c then do 25 i=1,n wrk0(i)=sx(i)*sx(i)*p(i) 25 continue c c solve l*y = (sx**2)*p c call forslv(nr,n,a,wrk0,wrk0) phip0=-dnrm2(n,wrk0,1)**2/rnwtln fstime=.false. c endif 28 phip=phip0 amulo=-phi/phip amuup=0.0d0 do 30 i=1,n amuup=amuup+(g(i)*g(i))/(sx(i)*sx(i)) 30 continue amuup=sqrt(amuup)/dlt done=.false. c$ write(ipr,956) amu c$ write(ipr,959) phi c$ write(ipr,960) phip c$ write(ipr,957) amulo c$ write(ipr,958) amuup c c test value of amu; generate next amu if necessary c 100 continue if(done) return c$ write(ipr,962) if(amu.ge.amulo .and. amu.le.amuup) go to 110 c if(amu.lt.amulo .or. amu.gt.amuup) c then amu=max(sqrt(amulo*amuup),amuup*1.0d-3) c$ write(ipr,956) amu c endif 110 continue c c copy (h,udiag) to l c where h <-- h+amu*(sx**2) [do not actually change (h,udiag)] do 130 j=1,n a(j,j)=udiag(j) + amu*sx(j)*sx(j) if(j.eq.n) go to 130 jp1=j+1 do 120 i=jp1,n a(i,j)=a(j,i) 120 continue 130 continue c c factor h=l(l+) c call choldc(nr,n,a,0.0d0,sqrt(epsm),addmax) c c solve h*p = l(l+)*sc = -g c do 140 i=1,n wrk0(i)=-g(i) 140 continue call lltslv(nr,n,a,sc,wrk0) c$ write(ipr,955) c$ write(ipr,963) (sc(i),i=1,n) c c reset h. note since udiag has not been destroyed we need do c nothing here. h is in the upper part and in udiag, still intact c stepln=0.0d0 do 150 i=1,n stepln=stepln + sx(i)*sx(i)*sc(i)*sc(i) 150 continue stepln=sqrt(stepln) phi=stepln-dlt do 160 i=1,n wrk0(i)=sx(i)*sx(i)*sc(i) 160 continue call forslv(nr,n,a,wrk0,wrk0) phip=-dnrm2(n,wrk0,1)**2/stepln c$ write(ipr,961) dlt,stepln c$ write(ipr,959) phi c$ write(ipr,960) phip if((alo*dlt.gt.stepln .or. stepln.gt.hi*dlt) .and. + (amuup-amulo.gt.0.0d0)) go to 170 c if((alo*dlt.le.stepln .and. stepln.le.hi*dlt) .or. c (amuup-amulo.le.0.0d0)) c then c c sc is acceptable hookstep c c$ write(ipr,954) done=.true. go to 100 c else c c sc not acceptable hookstep. select new amu c 170 continue c$ write(ipr,953) amulo=max(amulo,amu-(phi/phip)) if(phi.lt.0.0d0) amuup=min(amuup,amu) amu=amu-(stepln*phi)/(dlt*phip) c$ write(ipr,956) amu c$ write(ipr,957) amulo c$ write(ipr,958) amuup go to 100 c endif c endif c c%951 format(27h0hookst take newton step) c%952 format(32h0hookst newton step not taken) c%953 format(31h hookst sc is not acceptable) c%954 format(27h hookst sc is acceptable) c%955 format(28h hookst current step (sc)) c%956 format(18h hookst amu =,e20.13) c%957 format(18h hookst amulo =,e20.13) c%958 format(18h hookst amuup =,e20.13) c%959 format(18h hookst phi =,e20.13) c%960 format(18h hookst phip =,e20.13) c%961 format(18h hookst dlt =,e20.13/ c% + 18h hookst stepln=,e20.13) c%962 format(23h0hookst find new amu) c%963 format(14h hookst ,5(e20.13,3x)) end c subroutine hsnint c c purpose c c provide initial hessian when using secant updates c c parameters c c nr --> row dimension of matrix c n --> dimension of problem c a(n,n) <-- initial hessian (lower triangular matrix) c sx(n) --> diagonal scaling matrix for x c method --> algorithm to use to solve minimization problem c =1,2 factored secant method used c =3 unfactored secant method used c cstatic subroutine hsnint(nr,n,a,sx,method) implicit double precision (a-h,o-z) dimension a(nr,*),sx(n) c do 100 j=1,n if(method.eq.3) a(j,j)=sx(j)*sx(j) if(method.ne.3) a(j,j)=sx(j) if(j.eq.n) go to 100 jp1=j+1 do 90 i=jp1,n a(i,j)=0.0d0 90 continue 100 continue return end c subroutine lltslv c c purpose c c solve ax=b where a has the form l(l-transpose) c but only the lower triangular part, l, is stored. c c parameters c c nr --> row dimension of matrix c n --> dimension of problem c a(n,n) --> matrix of form l(l-transpose). c on return a is unchanged. c x(n) <-- solution vector c b(n) --> right-hand side vector c c note c c if b is not required by calling program, then c b and x may share the same storage. c cstatic subroutine lltslv(nr,n,a,x,b) implicit double precision (a-h,o-z) dimension a(nr,*),x(n),b(n) c c forward solve, result in x c call forslv(nr,n,a,x,b) c c back solve, result in x c call baksol(nr,n,a,x,x) return end c subroutine lnsrch c c purpose c c find a next newton iterate by line search. c c parameters c c n --> dimension of problem c x(n) --> old iterate: x[k-1] c f --> function value at old iterate, f(x) c g(n) --> gradient at old iterate, g(x), or approximate c p(n) --> non-zero newton step c xpls(n) <-- new iterate x[k] c fpls <-- function value at new iterate, f(xpls) c fcn --> name of subroutine to evaluate function c iretcd <-- return code c mxtake <-- boolean flag indicating step of maximum length used c stepmx --> maximum allowable step size c steptl --> relative step size at which successive iterates c considered close enough to terminate algorithm c sx(n) --> diagonal scaling matrix for x c ipr --> device to which to send output c c internal variables c c sln newton length c rln relative length of newton step c cstatic subroutine lnsrch(n,x,f,g,p,xpls,fpls,fcn,mxtake, + iretcd,stepmx,steptl,sx,ipr) implicit double precision (a-h,o-z) integer n,iretcd dimension sx(n) dimension x(n),g(n),p(n) dimension xpls(n) logical mxtake c ipr=ipr mxtake=.false. iretcd=2 c$ write(ipr,954) c$ write(ipr,955) (p(i),i=1,n) tmp=0.0d0 do 5 i=1,n tmp=tmp+sx(i)*sx(i)*p(i)*p(i) 5 continue sln=sqrt(tmp) if(sln.le.stepmx) go to 10 c c newton step longer than maximum allowed scl=stepmx/sln call sclmul(n,scl,p,p) sln=stepmx c$ write(ipr,954) c$ write(ipr,955) (p(i),i=1,n) 10 continue slp=ddot(n,g,1,p,1) rln=0.0d0 do 15 i=1,n rln=max(rln,abs(p(i))/max(abs(x(i)),1.0d0/sx(i))) 15 continue rmnlmb=steptl/rln almbda=1.0d0 c$ write(ipr,952) sln,slp,rmnlmb,stepmx,steptl c c loop c check if new iterate satisfactory. generate new lambda if necessary. c 100 continue if(iretcd.lt.2) return do 105 i=1,n xpls(i)=x(i) + almbda*p(i) 105 continue call fcn(n,xpls,fpls) c$ write(ipr,950) almbda c$ write(ipr,951) c$ write(ipr,955) (xpls(i),i=1,n) c$ write(ipr,953) fpls if(fpls.gt. f+slp*1.d-4*almbda) go to 130 c if(fpls.le. f+slp*1.d-4*almbda) c then c c solution found c iretcd=0 if(almbda.eq.1.0d0 .and. sln.gt. 0.99d0*stepmx) mxtake=.true. go to 100 c c solution not (yet) found c c else 130 if(almbda .ge. rmnlmb) go to 140 c if(almbda .lt. rmnlmb) c then c c no satisfactory xpls found sufficiently distinct from x c iretcd=1 go to 100 c else c c calculate new lambda c 140 if(almbda.ne.1.0d0) go to 150 c if(almbda.eq.1.0d0) c then c c first backtrack: quadratic fit c tlmbda=-slp/(2.0d0*(fpls-f-slp)) go to 170 c else c c all subsequent backtracks: cubic fit c 150 t1=fpls-f-almbda*slp t2=pfpls-f-plmbda*slp t3=1.0d0/(almbda-plmbda) a=t3*(t1/(almbda*almbda) - t2/(plmbda*plmbda)) b=t3*(t2*almbda/(plmbda*plmbda) + - t1*plmbda/(almbda*almbda) ) disc=b*b-3.0d0*a*slp if(disc.le. b*b) go to 160 c if(disc.gt. b*b) c then c c only one positive critical point, must be minimum c tlmbda=(-b+sign(1.0d0,a)*sqrt(disc))/(3.0d0*a) go to 165 c else c c both critical points positive, first is minimum c 160 tlmbda=(-b-sign(1.0d0,a)*sqrt(disc))/(3.0d0*a) c endif 165 if(tlmbda.gt. .5*almbda) tlmbda=.5*almbda c endif 170 plmbda=almbda pfpls=fpls if(tlmbda.ge. almbda*.1) go to 180 c if(tlmbda.lt.almbda/10.0d0) c then almbda=almbda*.1 go to 190 c else 180 almbda=tlmbda c endif c endif c endif 190 go to 100 c%950 format(18h lnsrch almbda=,e20.13) c%951 format(29h lnsrch new iterate (xpls)) c%952 format(18h lnsrch sln =,e20.13/ c% + 18h lnsrch slp =,e20.13/ c% + 18h lnsrch rmnlmb=,e20.13/ c% + 18h lnsrch stepmx=,e20.13/ c% + 18h lnsrch steptl=,e20.13) c%953 format(19h lnsrch f(xpls)=,e20.13) c%954 format(26h0lnsrch newton step (p)) c%955 format(14h lnsrch ,5(e20.13,3x)) end c subroutine mvmltl c c purpose c c compute y=lx c where l is a lower triangular matrix stored in a c c parameters c c nr --> row dimension of matrix c n --> dimension of problem c a(n,n) --> lower triangular (n*n) matrix c x(n) --> operand vector c y(n) <-- result vector c c note c c x and y cannot share storage c cstatic subroutine mvmltl(nr,n,a,x,y) implicit double precision (a-h,o-z) dimension a(nr,*),x(n),y(n) do 30 i=1,n sum=0.0d0 do 10 j=1,i sum=sum+a(i,j)*x(j) 10 continue y(i)=sum 30 continue return end c subroutine mvmlts c c purpose c c compute y=ax c where "a" is a symmetric (n*n) matrix stored in its lower c triangular part and x,y are n-vectors c c parameters c c nr --> row dimension of matrix c n --> dimension of problem c a(n,n) --> symmetric (n*n) matrix stored in c lower triangular part and diagonal c x(n) --> operand vector c y(n) <-- result vector c c note c c x and y cannot share storage. c cstatic subroutine mvmlts(nr,n,a,x,y) implicit double precision (a-h,o-z) dimension a(nr,*),x(n),y(n) do 30 i=1,n sum=0.0d0 do 10 j=1,i sum=sum+a(i,j)*x(j) 10 continue if(i.eq.n) go to 25 ip1=i+1 do 20 j=ip1,n sum=sum+a(j,i)*x(j) 20 continue 25 y(i)=sum 30 continue return end c subroutine mvmltu c c purpose c c compute y=(l+)x c where l is a lower triangular matrix stored in a c (l-transpose (l+) is taken implicitly) c c parameters c c nr --> row dimension of matrix c n --> dimension of problem c a(nr,1) --> lower triangular (n*n) matrix c x(n) --> operand vector c y(n) <-- result vector c c note c c x and y cannot share storage c cstatic subroutine mvmltu(nr,n,a,x,y) implicit double precision (a-h,o-z) dimension a(nr,*),x(n),y(n) do 30 i=1,n sum=0.0d0 do 10 j=i,n sum=sum+a(j,i)*x(j) 10 continue y(i)=sum 30 continue return end c subroutine optchk c c purpose c c check input for reasonableness c c parameters c c n --> dimension of problem c x(n) --> on entry, estimate to root of fcn c typsiz(n) <--> typical size of each component of x c sx(n) <-- diagonal scaling matrix for x c fscale <--> estimate of scale of objective function fcn c gradtl --> tolerance at which gradient considered close c enough to zero to terminate algorithm c itnlim <--> maximum number of allowable iterations c ndigit <--> number of good digits in optimization function fcn c epsm --> machine epsilon c dlt <--> trust region radius c method <--> algorithm indicator c iexp <--> expense flag c iagflg <--> =1 if analytic gradient supplied c iahflg <--> =1 if analytic hessian supplied c stepmx <--> maximum step size c msg <--> message and error code c ipr --> device to which to send output c cstatic subroutine optchk(n,x,typsiz,sx,fscale,gradtl,itnlim,ndigit,epsm, + dlt,method,iexp,iagflg,iahflg,stepmx,msg,ipr) implicit double precision (a-h,o-z) dimension x(n),typsiz(n),sx(n) c c check that parameters only take on acceptable values. c if not, set them to default values. if(method.lt.1 .or. method.gt.3) method=1 if(iagflg.ne.1) iagflg=0 if(iahflg.ne.1) iahflg=0 if(iexp.ne.0) iexp=1 if(mod(msg/2,2).eq.1 .and. iagflg.eq.0) go to 830 if(mod(msg/4,2).eq.1 .and. iahflg.eq.0) go to 835 c c check dimension of problem c if(n.le.0) go to 805 if(n.eq.1 .and. mod(msg,2).eq.0) go to 810 c c compute scale matrix c do 10 i=1,n if(typsiz(i).eq.0.0d0) typsiz(i)=1.0 if(typsiz(i).lt.0.0d0) typsiz(i)=-typsiz(i) sx(i)=1.0d0/typsiz(i) 10 continue c c check maximum step size c if (stepmx .gt. 0.0d0) go to 20 stpsiz = 0.0d0 do 15 i = 1, n stpsiz = stpsiz + x(i)*x(i)*sx(i)*sx(i) 15 continue stpsiz = sqrt(stpsiz) stepmx = max(1.0d3*stpsiz, 1.0d3) 20 continue c check function scale if(fscale.eq.0.0d0) fscale=1.0d0 if(fscale.lt.0.0d0) fscale=-fscale c c check gradient tolerance if(gradtl.lt.0.0d0) go to 815 c c check iteration limit if(itnlim.le.0) go to 820 c c check number of digits of accuracy in function fcn if(ndigit.eq.0) go to 825 if(ndigit.lt.0) ndigit=-log10(epsm) c c check trust region radius if(dlt.le.0.0d0) dlt=-1.0d0 if (dlt .gt. stepmx) dlt = stepmx return c c error exits c c%805 write(ipr,901) n c% msg=-1 805 msg=-1 go to 895 c%810 write(ipr,902) c% msg=-2 810 msg=-2 go to 895 c%815 write(ipr,903) gradtl c% msg=-3 815 msg=-3 go to 895 c%820 write(ipr,904) itnlim c% msg=-4 820 msg=-4 go to 895 c%825 write(ipr,905) ndigit c% msg=-5 825 msg=-5 go to 895 c%830 write(ipr,906) msg,iagflg c% msg=-6 830 msg=-6 go to 895 c%835 write(ipr,907) msg,iahflg c% msg=-7 835 msg=-7 895 return c%901 format(32h0optchk illegal dimension, n=,i5) c%902 format(55h0optchk +++ warning +++ this package is inefficient, c% + 26h for problems of size n=1./ c% + 48h optchk check installation libraries for more, c% + 22h appropriate routines./ c% + 41h optchk if none, set msg and resubmit.) c%903 format(38h0optchk illegal tolerance. gradtl=,e20.13) c%904 format(44h0optchk illegal iteration limit. itnlim=,i5) c%905 format(52h0optchk minimization function has no good digits., c% + 9h ndigit=,i5) c%906 format(50h0optchk user requests that analytic gradient be, c% + 33h accepted as properly coded (msg=,i5, 2h),/ c% + 45h optchk but analytic gradient not supplied, c% + 9h (iagflg=,i5, 2h).) c%907 format(49h0optchk user requests that analytic hessian be, c% + 33h accepted as properly coded (msg=,i5, 2h),/ c% + 44h optchk but analytic hessian not supplied, c% + 9h (iahflg=,i5, 2h).) end c subroutine optdrv c c purpose c c driver for non-linear optimization problem c c parameters c c nr --> row dimension of matrix c n --> dimension of problem c x(n) --> on entry: estimate to a root of fcn c fcn --> name of subroutine to evaluate optimization function c must be declared external in calling routine c fcn: r(n) --> r(1) c d1fcn --> (optional) name of subroutine to evaluate gradient c of fcn. must be declared external in calling routine c d2fcn --> (optional) name of subroutine to evaluate hessian of c of fcn. must be declared external in calling routine c typsiz(n) --> typical size for each component of x c fscale --> estimate of scale of objective function c method --> algorithm to use to solve minimization problem c =1 line search c =2 double dogleg c =3 more-hebdon c iexp --> =1 if optimization function fcn is expensive to c evaluate, =0 otherwise. if set then hessian will c be evaluated by secant update instead of c analytically or by finite differences c msg <--> on input: (.gt.0) message to inhibit certain c automatic checks c on output: (.lt.0) error code; =0 no error c ndigit --> number of good digits in optimization function fcn c itnlim --> maximum number of allowable iterations c iagflg --> =1 if analytic gradient supplied c iahflg --> =1 if analytic hessian supplied c ipr --> device to which to send output c dlt --> trust region radius c gradtl --> tolerance at which gradient considered close c enough to zero to terminate algorithm c stepmx --> maximum allowable step size c steptl --> relative step size at which successive iterates c considered close enough to terminate algorithm c xpls(n) <--> on exit: xpls is local minimum c fpls <--> on exit: function value at solution, xpls c gpls(n) <--> on exit: gradient at solution xpls c itrmcd <-- termination code c a(n,n) --> workspace for hessian (or estimate) c and its cholesky decomposition c udiag(n) --> workspace [for diagonal of hessian] c g(n) --> workspace (for gradient at current iterate) c p(n) --> workspace for step c sx(n) --> workspace (for diagonal scaling matrix) c wrk0(n) --> workspace c wrk1(n) --> workspace c wrk2(n) --> workspace c wrk3(n) --> workspace c itncnt current iteration, k {{was `internal'}} c c c internal variables c c analtl tolerance for comparison of estimated and c analytical gradients and hessians c epsm machine epsilon c f function value: fcn(x) c rnf relative noise in optimization function fcn. c noise=10.**(-ndigit) c cstatic subroutine optdrv(nr,n,x,fcn,d1fcn,d2fcn,typsiz,fscale, + method,iexp,msg,ndigit,itnlim,iagflg,iahflg,ipr, + dlt,gradtl,stepmx,steptl, + xpls,fpls,gpls,itrmcd, + a,udiag,g,p,sx,wrk0,wrk1,wrk2,wrk3,itncnt) implicit double precision (a-h,o-z) dimension x(n),xpls(n),g(n),gpls(n),p(n) dimension typsiz(n),sx(n) dimension a(nr,*),udiag(n) dimension wrk0(n),wrk1(n),wrk2(n),wrk3(n) logical mxtake,noupdt integer itncnt external fcn,d1fcn,d2fcn external result c for iteration tracing [defined in ../main/optimize.c] c c initialization c do 10 i=1,n p(i)=0.0d0 10 continue itncnt=0 iretcd=-1 epsm=d1mach(4) call optchk(n,x,typsiz,sx,fscale,gradtl,itnlim,ndigit,epsm, + dlt,method,iexp,iagflg,iahflg,stepmx,msg,ipr) if(msg.lt.0) return rnf=max(10.0d0**(-ndigit),epsm) analtl=max(1.0d-2,sqrt(rnf)) c c% if(mod(msg/8,2).eq.1) go to 15 c% write(ipr,901) c% write(ipr,900) (typsiz(i),i=1,n) c% write(ipr,902) c% write(ipr,900) (sx(i),i=1,n) c% write(ipr,903) fscale c% write(ipr,904) ndigit,iagflg,iahflg,iexp,method,itnlim,epsm c% write(ipr,905) stepmx,steptl,gradtl,dlt,rnf,analtl c% 15 continue c c evaluate fcn(x) c call fcn(n,x,f) c c evaluate analytic or finite difference gradient and check analytic c gradient, if requested. c if (iagflg .eq. 1) go to 20 c if (iagflg .eq. 0) c then call fstofd (1, 1, n, x, fcn, f, g, sx, rnf, wrk, 1) go to 25 c 20 call d1fcn (n, x, g) if (mod(msg/2,2) .eq. 1) go to 25 c if (mod(msg/2,2).eq.0) c then call grdchk (n, x, fcn, f, g, typsiz, sx, fscale, 1 rnf, analtl, wrk1, msg, ipr) if (msg .lt. 0) return 25 continue c call optstp(n,x,f,g,wrk1,itncnt,icscmx, + itrmcd,gradtl,steptl,sx,fscale,itnlim,iretcd,mxtake, + ipr,msg) if(itrmcd.ne.0) go to 700 c if(iexp.ne.1) go to 80 c c if optimization function expensive to evaluate (iexp=1), then c hessian will be obtained by secant updates. get initial hessian. c call hsnint(nr,n,a,sx,method) go to 90 80 continue c c evaluate analytic or finite difference hessian and check analytic c hessian if requested (only if user-supplied analytic hessian c routine d2fcn fills only lower triangular part and diagonal of a). c if (iahflg .eq. 1) go to 82 c if (iahflg .eq. 0) c then if (iagflg .eq. 1) call fstofd (nr, n, n, x, d1fcn, g, a, sx, 1 rnf, wrk1, 3) if (iagflg .ne. 1) call sndofd (nr, n, x, fcn, f, a, sx, rnf, 1 wrk1, wrk2) go to 88 c c else 82 if (mod(msg/4,2).eq.0) go to 85 c if (mod(msg/4, 2) .eq. 1) c then call d2fcn (nr, n, x, a) go to 88 c c else 85 call heschk (nr, n, x, fcn, d1fcn, d2fcn, f, g, a, typsiz, 1 sx, rnf, analtl, iagflg, udiag, wrk1, wrk2, msg, ipr) c c heschk evaluates d2fcn and checks it against the finite c difference hessian which it calculates by calling fstofd c (if iagflg .eq. 1) or sndofd (otherwise). c if (msg .lt. 0) return 88 continue c 90 if(mod(msg/8,2).eq.0) + call result(nr,n,x,f,g,a,p,itncnt,1,ipr) c c c iteration c 100 itncnt=itncnt+1 c c find perturbed local model hessian and its ll+ decomposition c (skip this step if line search or dogstep techniques being used with c secant updates. cholesky decomposition l already obtained from c secfac.) c if(iexp.eq.1 .and. method.ne.3) go to 105 103 call chlhsn(nr,n,a,epsm,sx,udiag) 105 continue c c solve for newton step: ap=-g c do 110 i=1,n wrk1(i)=-g(i) 110 continue call lltslv(nr,n,a,p,wrk1) c c decide whether to accept newton step xpls=x + p c or to choose xpls by a global strategy. c if (iagflg .ne. 0 .or. method .eq. 1) go to 111 dltsav = dlt if (method .eq. 2) go to 111 amusav = amu dlpsav = dltp phisav = phi phpsav = phip0 111 if(method.eq.1) + call lnsrch(n,x,f,g,p,xpls,fpls,fcn,mxtake,iretcd, + stepmx,steptl,sx,ipr) if(method.eq.2) + call dogdrv(nr,n,x,f,g,a,p,xpls,fpls,fcn,sx,stepmx, + steptl,dlt,iretcd,mxtake,wrk0,wrk1,wrk2,wrk3,ipr,itncnt) if(method.eq.3) + call hookdr(nr,n,x,f,g,a,udiag,p,xpls,fpls,fcn,sx,stepmx, + steptl,dlt,iretcd,mxtake,amu,dltp,phi,phip0,wrk0, + wrk1,wrk2,epsm,itncnt,ipr) c c if could not find satisfactory step and forward difference c gradient was used, retry using central difference gradient. c if (iretcd .ne. 1 .or. iagflg .ne. 0) go to 112 c if (iretcd .eq. 1 .and. iagflg .eq. 0) c then c c set iagflg for central differences c iagflg = -1 c% write(ipr,906) itncnt c call fstocd (n, x, fcn, sx, rnf, g) if (method .eq. 1) go to 105 dlt = dltsav if (method .eq. 2) go to 105 amu = amusav dltp = dlpsav phi = phisav phip0 = phpsav go to 103 c endif c c calculate step for output c 112 continue do 114 i = 1, n p(i) = xpls(i) - x(i) 114 continue c c calculate gradient at xpls c if (iagflg .eq. (-1)) go to 116 if (iagflg .eq. 0) go to 118 c c analytic gradient call d1fcn (n, xpls, gpls) go to 120 c c central difference gradient 116 call fstocd (n, xpls, fcn, sx, rnf, gpls) go to 120 c c forward difference gradient 118 call fstofd (1, 1, n, xpls, fcn, fpls, gpls, sx, rnf, wrk, 1) 120 continue c c check whether stopping criteria satisfied c call optstp(n,xpls,fpls,gpls,x,itncnt,icscmx, + itrmcd,gradtl,steptl,sx,fscale,itnlim,iretcd,mxtake, + ipr,msg) if(itrmcd.ne.0) go to 690 c c evaluate hessian at xpls c if(iexp.eq.0) go to 130 if(method.eq.3) + call secunf(nr,n,x,g,a,udiag,xpls,gpls,epsm,itncnt,rnf, + iagflg,noupdt,wrk1,wrk2,wrk3) if(method.ne.3) + call secfac(nr,n,x,g,a,xpls,gpls,epsm,itncnt,rnf,iagflg, + noupdt,wrk0,wrk1,wrk2,wrk3) go to 150 130 if(iahflg.eq.1) go to 140 if(iagflg.eq.1) + call fstofd(nr,n,n,xpls,d1fcn,gpls,a,sx,rnf,wrk1,3) if(iagflg.ne.1) call sndofd(nr,n,xpls,fcn,fpls,a,sx,rnf,wrk1,wrk2) go to 150 140 call d2fcn(nr,n,xpls,a) 150 continue if(mod(msg/16,2).eq.1) + call result(nr,n,xpls,fpls,gpls,a,p,itncnt,1,ipr) c c x <-- xpls and g <-- gpls and f <-- fpls c f=fpls do 160 i=1,n x(i)=xpls(i) g(i)=gpls(i) 160 continue go to 100 c c termination c c reset xpls,fpls,gpls, if previous iterate solution c 690 if(itrmcd.ne.3) go to 710 700 continue fpls=f do 705 i=1,n xpls(i)=x(i) gpls(i)=g(i) 705 continue c c print results c 710 continue if(mod(msg/8,2).eq.0) + call result(nr,n,xpls,fpls,gpls,a,p,itncnt,0,ipr) msg=0 return c c%900 format(14h optdrv ,5(e20.13,3x)) c%901 format(20h0optdrv typical x) c%902 format(40h optdrv diagonal scaling matrix for x) c%903 format(22h optdrv typical f =,e20.13) c%904 format(40h0optdrv number of good digits in fcn=,i5/ c% + 27h optdrv gradient flag =,i5,18h (=1 if analytic, c% + 19h gradient supplied)/ c% + 27h optdrv hessian flag =,i5,18h (=1 if analytic, c% + 18h hessian supplied)/ c% + 27h optdrv expense flag =,i5, 9h (=1 if, c% + 45h minimization function expensive to evaluate)/ c% + 27h optdrv method to use =,i5,19h (=1,2,3 for line, c% + 49h search, double dogleg, more-hebdon respectively)/ c% + 27h optdrv iteration limit=,i5/ c% + 27h optdrv machine epsilon=,e20.13) c%905 format(30h0optdrv maximum step size =,e20.13/ c% + 30h optdrv step tolerance =,e20.13/ c% + 30h optdrv gradient tolerance=,e20.13/ c% + 30h optdrv trust reg radius =,e20.13/ c% + 30h optdrv rel noise in fcn =,e20.13/ c% + 30h optdrv anal-fd tolerance =,e20.13) c%906 format(52h optdrv shift from forward to central differences, c% 1 14h in iteration , i5) end c subroutine optif0 c c purpose c c provide simplest interface to minimization package. c user has no control over options. c c parameters c c nr --> row dimension of matrix c n --> dimension of problem c x(n) --> initial estimate of minimum c fcn --> name of routine to evaluate minimization function. c must be declared external in calling routine. c xpls(n) <-- local minimum c fpls <-- function value at local minimum xpls c gpls(n) <-- gradient at local minimum xpls c itrmcd <-- termination code c a(n,n) --> workspace c wrk(n,9) --> workspace c subroutine optif0(nr,n,x,fcn,xpls,fpls,gpls,itrmcd,a,wrk) implicit double precision (a-h,o-z) dimension x(n),xpls(n),gpls(n) dimension a(nr,*),wrk(nr,9) external fcn,d1fcn,d2fcn integer itncnt c c equivalence wrk(n,1) = udiag(n) c wrk(n,2) = g(n) c wrk(n,3) = p(n) c wrk(n,4) = typsiz(n) c wrk(n,5) = sx(n) c wrk(n,6) = wrk0(n) c wrk(n,7) = wrk1(n) c wrk(n,8) = wrk2(n) c wrk(n,9) = wrk3(n) c call dfault(n,x,wrk(1,4),fscale,method,iexp,msg,ndigit, + itnlim,iagflg,iahflg,ipr,dlt,gradtl,stepmx,steptl) call optdrv(nr,n,x,fcn,d1fcn,d2fcn,wrk(1,4),fscale, + method,iexp,msg,ndigit,itnlim,iagflg,iahflg,ipr, + dlt,gradtl,stepmx,steptl, + xpls,fpls,gpls,itrmcd, + a,wrk(1,1),wrk(1,2),wrk(1,3),wrk(1,5),wrk(1,6), + wrk(1,7),wrk(1,8),wrk(1,9),itncnt) return end C---- this one is called from ../main/optimize.c : --------------- c subroutine optif9 c c purpose c c provide complete interface to minimization package. c user has full control over options. c c parameters c c nr --> row dimension of matrix c n --> dimension of problem c x(n) --> on entry: estimate to a root of fcn c fcn --> name of subroutine to evaluate optimization function c must be declared external in calling routine c fcn: r(n) --> r(1) c d1fcn --> (optional) name of subroutine to evaluate gradient c of fcn. must be declared external in calling routine c d2fcn --> (optional) name of subroutine to evaluate hessian of c of fcn. must be declared external in calling routine c typsiz(n) --> typical size for each component of x c fscale --> estimate of scale of objective function c method --> algorithm to use to solve minimization problem c =1 line search c =2 double dogleg c =3 more-hebdon c iexp --> =1 if optimization function fcn is expensive to c evaluate, =0 otherwise. if set then hessian will c be evaluated by secant update instead of c analytically or by finite differences c msg <--> on input: (.gt.0) message to inhibit certain c automatic checks c on output: (.lt.0) error code; =0 no error c ndigit --> number of good digits in optimization function fcn c itnlim --> maximum number of allowable iterations c iagflg --> =1 if analytic gradient supplied c iahflg --> =1 if analytic hessian supplied c ipr --> device to which to send output c dlt --> trust region radius c gradtl --> tolerance at which gradient considered close c enough to zero to terminate algorithm c stepmx --> maximum allowable step size c steptl --> relative step size at which successive iterates c considered close enough to terminate algorithm c xpls(n) <--> on exit: xpls is local minimum c fpls <--> on exit: function value at solution, xpls c gpls(n) <--> on exit: gradient at solution xpls c itrmcd <-- termination code c a(n,n) --> workspace for hessian (or estimate) c and its cholesky decomposition c wrk(n,8) --> workspace c itncnt --> iteration count c subroutine optif9(nr,n,x,fcn,d1fcn,d2fcn,typsiz,fscale, + method,iexp,msg,ndigit,itnlim,iagflg,iahflg,ipr, + dlt,gradtl,stepmx,steptl, + xpls,fpls,gpls,itrmcd,a,wrk, itncnt) implicit double precision (a-h,o-z) dimension x(n),xpls(n),gpls(n),typsiz(n) dimension a(nr,*),wrk(nr,8) external fcn,d1fcn,d2fcn integer itncnt c c equivalence wrk(n,1) = udiag(n) c wrk(n,2) = g(n) c wrk(n,3) = p(n) c wrk(n,4) = sx(n) c wrk(n,5) = wrk0(n) c wrk(n,6) = wrk1(n) c wrk(n,7) = wrk2(n) c wrk(n,8) = wrk3(n) c call optdrv(nr,n,x,fcn,d1fcn,d2fcn,typsiz,fscale, + method,iexp,msg,ndigit,itnlim,iagflg,iahflg,ipr, + dlt,gradtl,stepmx,steptl, + xpls,fpls,gpls,itrmcd, + a,wrk(1,1),wrk(1,2),wrk(1,3),wrk(1,4),wrk(1,5), + wrk(1,6),wrk(1,7),wrk(1,8),itncnt) return end c subroutine optstp c c unconstrained minimization stopping criteria c c find whether the algorithm should terminate, due to any c of the following: c 1) problem solved within user tolerance c 2) convergence within user tolerance c 3) iteration limit reached c 4) divergence or too restrictive maximum step (stepmx) suspected c c parameters c c n --> dimension of problem c xpls(n) --> new iterate x[k] c fpls --> function value at new iterate f(xpls) c gpls(n) --> gradient at new iterate, g(xpls), or approximate c x(n) --> old iterate x[k-1] c itncnt --> current iteration k c icscmx <--> number consecutive steps .ge. stepmx c [retain value between successive calls] c itrmcd <-- termination code c gradtl --> tolerance at which relative gradient considered close c enough to zero to terminate algorithm c steptl --> relative step size at which successive iterates c considered close enough to terminate algorithm c sx(n) --> diagonal scaling matrix for x c fscale --> estimate of scale of objective function c itnlim --> maximum number of allowable iterations c iretcd --> return code c mxtake --> boolean flag indicating step of maximum length used c ipr --> device to which to send output c msg --> if msg includes a term 8, suppress output c c cstatic subroutine optstp(n,xpls,fpls,gpls,x,itncnt,icscmx, + itrmcd,gradtl,steptl,sx,fscale,itnlim,iretcd,mxtake,ipr,msg) implicit double precision (a-h,o-z) integer n,itncnt,icscmx,itrmcd,itnlim dimension sx(n) dimension xpls(n),gpls(n),x(n) logical mxtake c itrmcd=0 c c last global step failed to locate a point lower than x c if(iretcd.ne.1) go to 50 c if(iretcd.eq.1) c then jtrmcd=3 go to 600 c endif 50 continue c c find direction in which relative gradient maximum. c check whether within tolerance c d=max(abs(fpls),fscale) rgx=0.0d0 do 100 i=1,n relgrd=abs(gpls(i))*max(abs(xpls(i)),1.0d0/sx(i))/d rgx=max(rgx,relgrd) 100 continue jtrmcd=1 if(rgx.le.gradtl) go to 600 c if(itncnt.eq.0) return c c find direction in which relative stepsize maximum c check whether within tolerance. c rsx=0.0d0 do 120 i=1,n relstp=abs(xpls(i)-x(i))/max(abs(xpls(i)),1.0d0/sx(i)) rsx=max(rsx,relstp) 120 continue jtrmcd=2 if(rsx.le.steptl) go to 600 c c check iteration limit c jtrmcd=4 if(itncnt.ge.itnlim) go to 600 c c check number of consecutive steps \ stepmx c if(mxtake) go to 140 c if(.not.mxtake) c then icscmx=0 return c else 140 continue c% if (mod(msg/8,2) .eq. 0) write(ipr,900) icscmx=icscmx+1 if(icscmx.lt.5) return jtrmcd=5 c endif c c c print termination code c 600 itrmcd=jtrmcd c% if (mod(msg/8,2) .eq. 0) go to(601,602,603,604,605), itrmcd c% go to 700 c%601 write(ipr,901) c% go to 700 c%602 write(ipr,902) c% go to 700 c%603 write(ipr,903) c% go to 700 c%604 write(ipr,904) c% go to 700 c%605 write(ipr,905) c 700 return c c%900 format(48h0optstp step of maximum length (stepmx) taken) c%901 format(43h0optstp relative gradient close to zero./ c% + 48h optstp current iterate is probably solution.) c%902 format(48h0optstp successive iterates within tolerance./ c% + 48h optstp current iterate is probably solution.) c%903 format(52h0optstp last global step failed to locate a point, c% + 14h lower than x./ c% + 51h optstp either x is an approximate local minimum, c% + 17h of the function,/ c% + 50h optstp the function is too non-linear for this, c% + 11h algorithm,/ c% + 34h optstp or steptl is too large.) c%904 format(36h0optstp iteration limit exceeded./ c% + 28h optstp algorithm failed.) c%905 format(39h0optstp maximum step size exceeded 5, c% + 19h consecutive times./ c% + 50h optstp either the function is unbounded below,/ c% + 47h optstp becomes asymptotic to a finite value, c% + 30h from above in some direction,/ c% + 33h optstp or stepmx is too small) end c subroutine qraux1 c c purpose c c interchange rows i,i+1 of the upper hessenberg matrix r, c columns i to n c c parameters c c nr --> row dimension of matrix c n --> dimension of matrix c r(n,n) <--> upper hessenberg matrix c i --> index of row to interchange (i.lt.n) c cstatic subroutine qraux1(nr,n,r,i) implicit double precision (a-h,o-z) dimension r(nr,*) do 10 j=i,n tmp=r(i,j) r(i,j)=r(i+1,j) r(i+1,j)=tmp 10 continue return end c subroutine qraux2 c c purpose c c pre-multiply r by the jacobi rotation j(i,i+1,a,b) c c parameters c c nr --> row dimension of matrix c n --> dimension of matrix c r(n,n) <--> upper hessenberg matrix c i --> index of row c a --> scalar c b --> scalar c cstatic subroutine qraux2(nr,n,r,i,a,b) implicit double precision (a-h,o-z) dimension r(nr,*) den=sqrt(a*a + b*b) c=a/den s=b/den do 10 j=i,n y=r(i,j) z=r(i+1,j) r(i,j)=c*y - s*z r(i+1,j)=s*y + c*z 10 continue return end c subroutine qrupdt c c purpose c c find an orthogonal (n*n) matrix (q*) and an upper triangular (n*n) c matrix (r*) such that (q*)(r*)=r+u(v+) c c parameters c c nr --> row dimension of matrix c n --> dimension of problem c a(n,n) <--> on input: contains r c on output: contains (r*) c u(n) --> vector c v(n) --> vector c cstatic subroutine qrupdt(nr,n,a,u,v) implicit double precision (a-h,o-z) dimension a(nr,*) dimension u(n),v(n) c c determine last non-zero in u(.) c k=n 10 if(u(k).ne.0.0d0 .or. k.eq.1) go to 20 c if(u(k).eq.0.0d0 .and. k.gt.1) c then k=k-1 go to 10 c endif c c (k-1) jacobi rotations transform c r + u(v+) --> (r*) + (u(1)*e1)(v+) c which is upper hessenberg c 20 if(k.le.1) go to 40 km1=k-1 do 30 ii=1,km1 i=km1-ii+1 if(u(i).ne.0.0d0) go to 25 c if(u(i).eq.0.0d0) c then call qraux1(nr,n,a,i) u(i)=u(i+1) go to 30 c else 25 call qraux2(nr,n,a,i,u(i),-u(i+1)) u(i)=sqrt(u(i)*u(i) + u(i+1)*u(i+1)) c endif 30 continue c endif c c r <-- r + (u(1)*e1)(v+) c 40 do 50 j=1,n a(1,j)=a(1,j) +u(1)*v(j) 50 continue c c (k-1) jacobi rotations transform upper hessenberg r c to upper triangular (r*) c if(k.le.1) go to 100 km1=k-1 do 80 i=1,km1 if(a(i,i).ne.0.0d0) go to 70 c if(a(i,i).eq.0.0d0) c then call qraux1(nr,n,a,i) go to 80 c else 70 t1=a(i,i) t2=-a(i+1,i) call qraux2(nr,n,a,i,t1,t2) c endif 80 continue c endif 100 return end c subroutine sclmul c c purpose c c multiply vector by scalar c result vector may be operand vector c c parameters c c n --> dimension of vectors c s --> scalar c v(n) --> operand vector c z(n) <-- result vector cstatic subroutine sclmul(n,s,v,z) implicit double precision (a-h,o-z) dimension v(n),z(n) do 100 i=1,n z(i)=s*v(i) 100 continue return end c subroutine secfac c c purpose c c update hessian by the bfgs factored method c c parameters c c nr --> row dimension of matrix c n --> dimension of problem c x(n) --> old iterate, x[k-1] c g(n) --> gradient or approximate at old iterate c a(n,n) <--> on entry: cholesky decomposition of hessian in c lower part and diagonal. c on exit: updated cholesky decomposition of hessian c in lower triangular part and diagonal c xpls(n) --> new iterate, x[k] c gpls(n) --> gradient or approximate at new iterate c epsm --> machine epsilon c itncnt --> iteration count c rnf --> relative noise in optimization function fcn c iagflg --> =1 if analytic gradient supplied, =0 itherwise c noupdt <--> boolean: no update yet c [retain value between successive calls] c s(n) --> workspace c y(n) --> workspace c u(n) --> workspace c w(n) --> workspace c cstatic subroutine secfac(nr,n,x,g,a,xpls,gpls,epsm,itncnt,rnf, + iagflg,noupdt,s,y,u,w) implicit double precision (a-h,o-z) dimension x(n),xpls(n),g(n),gpls(n) dimension a(nr,*) dimension s(n),y(n),u(n),w(n) logical noupdt,skpupd c if(itncnt.eq.1) noupdt=.true. do 10 i=1,n s(i)=xpls(i)-x(i) y(i)=gpls(i)-g(i) 10 continue den1=ddot(n,s,1,y,1) snorm2=dnrm2(n,s,1) ynrm2=dnrm2(n,y,1) if(den1.lt.sqrt(epsm)*snorm2*ynrm2) go to 110 c if(den1.ge.sqrt(epsm)*snorm2*ynrm2) c then call mvmltu(nr,n,a,s,u) den2=ddot(n,u,1,u,1) c c l <-- sqrt(den1/den2)*l c alp=sqrt(den1/den2) if(.not.noupdt) go to 50 c if(noupdt) c then do 30 j=1,n u(j)=alp*u(j) do 20 i=j,n a(i,j)=alp*a(i,j) 20 continue 30 continue noupdt=.false. den2=den1 alp=1.0d0 c endif 50 skpupd=.true. c c w = l(l+)s = hs c call mvmltl(nr,n,a,u,w) i=1 if(iagflg.ne.0) go to 55 c if(iagflg.eq.0) c then reltol=sqrt(rnf) go to 60 c else 55 reltol=rnf c endif 60 if(i.gt.n .or. .not.skpupd) go to 70 c if(i.le.n .and. skpupd) c then if(abs(y(i)-w(i)) .lt. reltol*max(abs(g(i)),abs(gpls(i)))) + go to 65 c if(abs(y(i)-w(i)) .ge. reltol*dmax1(abs(g(i)),abs(gpls(i)))) c then skpupd=.false. go to 60 c else 65 i=i+1 go to 60 c endif c endif 70 if(skpupd) go to 110 c if(.not.skpupd) c then c c w=y-alp*l(l+)s c do 75 i=1,n w(i)=y(i)-alp*w(i) 75 continue c c alp=1/sqrt(den1*den2) c alp=alp/den1 c c u=(l+)/sqrt(den1*den2) = (l+)s/sqrt((y+)s * (s+)l(l+)s) c do 80 i=1,n u(i)=alp*u(i) 80 continue c c copy l into upper triangular part. zero l. c if(n.eq.1) go to 93 do 90 i=2,n im1=i-1 do 85 j=1,im1 a(j,i)=a(i,j) a(i,j)=0.0d0 85 continue 90 continue c c find q, (l+) such that q(l+) = (l+) + u(w+) c 93 call qrupdt(nr,n,a,u,w) c c upper triangular part and diagonal of a now contain updated c cholesky decomposition of hessian. copy back to lower c triangular part. c if(n.eq.1) go to 110 do 100 i=2,n im1=i-1 do 95 j=1,im1 a(i,j)=a(j,i) 95 continue 100 continue c endif c endif 110 return end c subroutine secunf c c purpose c c update hessian by the bfgs unfactored method c c parameters c c nr --> row dimension of matrix c n --> dimension of problem c x(n) --> old iterate, x[k-1] c g(n) --> gradient or approximate at old iterate c a(n,n) <--> on entry: approximate hessian at old iterate c in upper triangular part (and udiag) c on exit: updated approx hessian at new iterate c in lower triangular part and diagonal c [lower triangular part of symmetric matrix] c udiag --> on entry: diagonal of hessian c xpls(n) --> new iterate, x[k] c gpls(n) --> gradient or approximate at new iterate c epsm --> machine epsilon c itncnt --> iteration count c rnf --> relative noise in optimization function fcn c iagflg --> =1 if analytic gradient supplied, =0 otherwise c noupdt <--> boolean: no update yet c [retain value between successive calls] c s(n) --> workspace c y(n) --> workspace c t(n) --> workspace c cstatic subroutine secunf(nr,n,x,g,a,udiag,xpls,gpls,epsm,itncnt, + rnf,iagflg,noupdt,s,y,t) implicit double precision (a-h,o-z) dimension x(n),g(n),xpls(n),gpls(n) dimension a(nr,*) dimension udiag(n) dimension s(n),y(n),t(n) logical noupdt,skpupd c c copy hessian in upper triangular part and udiag to c lower triangular part and diagonal c do 5 j=1,n a(j,j)=udiag(j) if(j.eq.n) go to 5 jp1=j+1 do 4 i=jp1,n a(i,j)=a(j,i) 4 continue 5 continue c if(itncnt.eq.1) noupdt=.true. do 10 i=1,n s(i)=xpls(i)-x(i) y(i)=gpls(i)-g(i) 10 continue den1=ddot(n,s,1,y,1) snorm2=dnrm2(n,s,1) ynrm2=dnrm2(n,y,1) if(den1.lt.sqrt(epsm)*snorm2*ynrm2) go to 100 c if(den1.ge.sqrt(epsm)*snorm2*ynrm2) c then call mvmlts(nr,n,a,s,t) den2=ddot(n,s,1,t,1) if(.not. noupdt) go to 50 c if(noupdt) c then c c h <-- [(s+)y/(s+)hs]h c gam=den1/den2 den2=gam*den2 do 30 j=1,n t(j)=gam*t(j) do 20 i=j,n a(i,j)=gam*a(i,j) 20 continue 30 continue noupdt=.false. c endif 50 skpupd=.true. c c check update condition on row i c do 60 i=1,n tol=rnf*max(abs(g(i)),abs(gpls(i))) if(iagflg.eq.0) tol=tol/sqrt(rnf) if(abs(y(i)-t(i)).lt.tol) go to 60 c if(abs(y(i)-t(i)).ge.tol) c then skpupd=.false. go to 70 c endif 60 continue 70 if(skpupd) go to 100 c if(.not.skpupd) c then c c bfgs update c do 90 j=1,n do 80 i=j,n a(i,j)=a(i,j)+y(i)*y(j)/den1-t(i)*t(j)/den2 80 continue 90 continue c endif c endif 100 return end c subroutine sndofd c c purpose c c find second order forward finite difference approximation "a" c to the second derivative (hessian) of the function defined by the subp c "fcn" evaluated at the new iterate "xpls" c c for optimization use this routine to estimate c 1) the second derivative (hessian) of the optimization function c if no analytical user function has been supplied for either c the gradient or the hessian and if the optimization function c "fcn" is inexpensive to evaluate. c c parameters c c nr --> row dimension of matrix c n --> dimension of problem c xpls(n) --> new iterate: x[k] c fcn --> name of subroutine to evaluate function c fpls --> function value at new iterate, f(xpls) c a(n,n) <-- finite difference approximation to hessian c only lower triangular matrix and diagonal c are returned c sx(n) --> diagonal scaling matrix for x c rnoise --> relative noise in fname [f(x)] c stepsz(n) --> workspace (stepsize in i-th component direction) c anbr(n) --> workspace (neighbor in i-th direction) c c cstatic subroutine sndofd(nr,n,xpls,fcn,fpls,a,sx,rnoise,stepsz,anbr) implicit double precision (a-h,o-z) dimension xpls(n) dimension sx(n) dimension stepsz(n),anbr(n) dimension a(nr,*) c c find i-th stepsize and evaluate neighbor in direction c of i-th unit vector. c ov3 = 1.0d0/3.0d0 do 10 i=1,n stepsz(i)=rnoise**ov3 * max(abs(xpls(i)),1.0d0/sx(i)) xtmpi=xpls(i) xpls(i)=xtmpi+stepsz(i) call fcn(n,xpls,anbr(i)) xpls(i)=xtmpi 10 continue c c calculate column i of a c do 30 i=1,n xtmpi=xpls(i) xpls(i)=xtmpi+2.0d0*stepsz(i) call fcn(n,xpls,fhat) a(i,i)=((fpls-anbr(i))+(fhat-anbr(i)))/(stepsz(i)*stepsz(i)) c c calculate sub-diagonal elements of column if(i.eq.n) go to 25 xpls(i)=xtmpi+stepsz(i) ip1=i+1 do 20 j=ip1,n xtmpj=xpls(j) xpls(j)=xtmpj+stepsz(j) call fcn(n,xpls,fhat) a(j,i)=((fpls-anbr(i))+(fhat-anbr(j)))/(stepsz(i)*stepsz(j)) xpls(j)=xtmpj 20 continue 25 xpls(i)=xtmpi 30 continue return end c subroutine tregup c c purpose c c decide whether to accept xpls=x+sc as the next iterate and update the c trust region dlt. c c parameters c c nr --> row dimension of matrix c n --> dimension of problem c x(n) --> old iterate x[k-1] c f --> function value at old iterate, f(x) c g(n) --> gradient at old iterate, g(x), or approximate c a(n,n) --> cholesky decomposition of hessian in c lower triangular part and diagonal. c hessian or approx in upper triangular part c fcn --> name of subroutine to evaluate function c sc(n) --> current step c sx(n) --> diagonal scaling matrix for x c nwtake --> boolean, =.true. if newton step taken c stepmx --> maximum allowable step size c steptl --> relative step size at which successive iterates c considered close enough to terminate algorithm c dlt <--> trust region radius c iretcd <--> return code c =0 xpls accepted as next iterate; c dlt trust region for next iteration. c =1 xpls unsatisfactory but accepted as next iterate c because xpls-x .lt. smallest allowable c step length. c =2 f(xpls) too large. continue current iteration c with new reduced dlt. c =3 f(xpls) sufficiently small, but quadratic model c predicts f(xpls) sufficiently well to continue c current iteration with new doubled dlt. c xplsp(n) <--> workspace [value needs to be retained between c succesive calls of k-th global step] c fplsp <--> [retain value between successive calls] c xpls(n) <-- new iterate x[k] c fpls <-- function value at new iterate, f(xpls) c mxtake <-- boolean flag indicating step of maximum length used c ipr --> device to which to send output c method --> algorithm to use to solve minimization problem c =1 line search c =2 double dogleg c =3 more-hebdon c udiag(n) --> diagonal of hessian in a(.,.) c cstatic subroutine tregup(nr,n,x,f,g,a,fcn,sc,sx,nwtake,stepmx,steptl, + dlt,iretcd,xplsp,fplsp,xpls,fpls,mxtake,ipr,method,udiag) implicit double precision (a-h,o-z) dimension x(n),xpls(n),g(n) dimension sx(n),sc(n),xplsp(n) dimension a(nr,*) logical nwtake,mxtake dimension udiag(n) c ipr=ipr mxtake=.false. do 100 i=1,n xpls(i)=x(i)+sc(i) 100 continue call fcn(n,xpls,fpls) dltf=fpls-f slp=ddot(n,g,1,sc,1) c c next statement added for case of compilers which do not optimize c evaluation of next "if" statement (in which case fplsp could be c undefined). c if(iretcd.eq.4) fplsp=0.0d0 c$ write(ipr,961) iretcd,fpls,fplsp,dltf,slp if(iretcd.ne.3 .or. (fpls.lt.fplsp .and. dltf.le. 1.d-4*slp)) + go to 130 c if(iretcd.eq.3 .and. (fpls.ge.fplsp .or. dltf.gt. 1.d-4*slp)) c then c c reset xpls to xplsp and terminate global step c iretcd=0 do 110 i=1,n xpls(i)=xplsp(i) 110 continue fpls=fplsp dlt=.5*dlt c$ write(ipr,951) go to 230 c else c c fpls too large c 130 if(dltf.le. 1.d-4*slp) go to 170 c if(dltf.gt. 1.d-4*slp) c then c$ write(ipr,952) rln=0.0d0 do 140 i=1,n rln=max(rln,abs(sc(i))/max(abs(xpls(i)),1.0d0/sx(i))) 140 continue c$ write(ipr,962) rln if(rln.ge.steptl) go to 150 c if(rln.lt.steptl) c then c c cannot find satisfactory xpls sufficiently distinct from x c iretcd=1 c$ write(ipr,954) go to 230 c else c c reduce trust region and continue global step c 150 iretcd=2 dltmp=-slp*dlt/(2.0d0*(dltf-slp)) c$ write(ipr,963) dltmp if(dltmp.ge. .1*dlt) go to 155 c if(dltmp.lt. .1*dlt) c then dlt=.1*dlt go to 160 c else 155 dlt=dltmp c endif 160 continue c$ write(ipr,955) go to 230 c endif c else c c fpls sufficiently small c 170 continue c$ write(ipr,958) dltfp=0.0d0 if (method .eq. 3) go to 180 c c if (method .eq. 2) c then c do 177 i = 1, n temp = 0.0d0 do 173 j = i, n temp = temp + (a(j, i)*sc(j)) 173 continue dltfp = dltfp + temp*temp 177 continue go to 190 c c else c 180 do 187 i = 1, n dltfp = dltfp + udiag(i)*sc(i)*sc(i) if (i .eq. n) go to 187 temp = 0 ip1 = i + 1 do 183 j = ip1, n temp = temp + a(i, j)*sc(i)*sc(j) 183 continue dltfp = dltfp + 2.0d0*temp 187 continue c c end if c 190 dltfp = slp + dltfp/2.0d0 c$ write(ipr,964) dltfp,nwtake if(iretcd.eq.2 .or. (abs(dltfp-dltf).gt. .1*abs(dltf)) + .or. nwtake .or. (dlt.gt. .99*stepmx)) go to 210 c if(iretcd.ne.2 .and. (abs(dltfp-dltf) .le. .1*abs(dltf)) c + .and. (.not.nwtake) .and. (dlt.le. .99*stepmx)) c then c c double trust region and continue global step c iretcd=3 do 200 i=1,n xplsp(i)=xpls(i) 200 continue fplsp=fpls dlt=min(2.0d0*dlt,stepmx) c$ write(ipr,959) go to 230 c else c c accept xpls as next iterate. choose new trust region. c 210 continue c$ write(ipr,960) iretcd=0 if(dlt.gt. .99*stepmx) mxtake=.true. if(dltf.lt. .1*dltfp) go to 220 c if(dltf.ge. .1*dltfp) c then c c decrease trust region for next iteration c dlt=.5*dlt go to 230 c else c c check whether to increase trust region for next iteration c 220 if(dltf.le. .75*dltfp) dlt=min(2.0d0*dlt,stepmx) c endif c endif c endif c endif 230 continue c$ write(ipr,953) c$ write(ipr,956) iretcd,mxtake,dlt,fpls c$ write(ipr,957) c$ write(ipr,965) (xpls(i),i=1,n) return c c%951 format(55h tregup reset xpls to xplsp. termination global step) c%952 format(26h tregup fpls too large.) c%953 format(38h0tregup values after call to tregup) c%954 format(54h tregup cannot find satisfactory xpls distinct from, c% + 27h x. terminate global step.) c%955 format(53h tregup reduce trust region. continue global step.) c%956 format(21h tregup iretcd=,i3/ c% + 21h tregup mxtake=,l1/ c% + 21h tregup dlt =,e20.13/ c% + 21h tregup fpls =,e20.13) c%957 format(32h tregup new iterate (xpls)) c%958 format(35h tregup fpls sufficiently small.) c%959 format(54h tregup double trust region. continue global step.) c%960 format(50h tregup accept xpls as new iterate. choose new, c% + 38h trust region. terminate global step.) c%961 format(18h tregup iretcd=,i5/ c% + 18h tregup fpls =,e20.13/ c% + 18h tregup fplsp =,e20.13/ c% + 18h tregup dltf =,e20.13/ c% + 18h tregup slp =,e20.13) c%962 format(18h tregup rln =,e20.13) c%963 format(18h tregup dltmp =,e20.13) c%964 format(18h tregup dltfp =,e20.13/ c% + 18h tregup nwtake=,l1) c%965 format(14h tregup ,5(e20.13,3x)) end