c This Fortran version of PRAXIS c (supplied by Linda Kaufman, lck@research.att.com) c calls Port utilities to get scratch storage. c c R. Brent (rpb@phys4.anu.oz.au), 31 July 1990. c subroutine fminn(t0,h0,n,x,f,fmin) real t0,h0,x(n),f,fmin external f c c fminn returns the minimum of the function f(x,n) of n variables c using the principal axis method. The gradient of the function is c not required. c c For a description of the algorithm, see chapter 7 of c "Algorithms for finding zeros and extrema of functions without c calculating derivatives" by Richard P Brent (Prentice-Hall, 1973). c c The original Algol W version of praxis was translated into Fortran c by Linda Kaufman. c c The input parameters are c c t0 is a tolerance. fminn attempts to return fmin=f(x) c such that if x0 is the true local minimum near x, then c norm(x-x0) ) t0 + squareroot(machep)*norm(x). c where machep is the machine precision c h0 is the maximum step size. h0 should be set to about the c maximum distance from the initial guess to the minimum. c (if h0 is set too large or too small, the initial rate of c convergence may be slow.) c n (at least two) is the number of variables upon which c the function depends. c x is an array containing on entry a guess of the point of c minimum, on return the estimated point of minimum. c f(x,n) is the function to be minimized. f should be a real c function declared external in the calling program. c c output parameters are c c fmin the value of the function at the minimum c x the solution of the problem c the approximating quadratic form is c q(x*) = f(x,n) + (1/2) * (x*-x)-transpose * a * (x*-x) c where x is the best estimate of the minimum and a is c inverse(v-transpose) * d * inverse(v) c (v(*,*) is the matrix of search directions\ d(*) is the array c of second differences). if f has continuous second derivatives c near x0, a will tend to the hessian of f at x0 as x approaches x0. c c It is assumed that on floating-point underflow the result is set c to zero. c The user should observe the comment on heuristic numbers c c Error condition* c n.lt.1 fatal c c Extra storage required: 7n+n*n real storage locations c logical illc integer ktm real scbd common /cstak/ d double precision d(500) real r(1000) equivalence (d(1),r(1)) call enter(1) if (n.lt.1) call seterr(12hfminn-n.lt.1,12,1,2) id=istkgt(n,3) iy=istkgt(n,3) iz= istkgt(n,3) iq1=istkgt(n,3) iq0=istkgt(n,3) ie=istkgt(n,3) iv=istkgt(n*n,3) c Heuristic numbers c If the axes are badly scaled,set scbd to 10, otherwise c set it to 1. If the problem is known to be ill conditioned, c set illc to true, otherwise false. c ktm is the number of iterations without improvement c before the algorithm terminates. ktm is 4 is c very cautious/usually ktm=1 is sufficient. scbd=1.0 illc=.false. ktm=1 eps= r1mach(4) call praxis(t0,eps,h0,n,x,f,fmin,scbd,ktm,illc, 1 r(id),r(iy),r(iz),r(iq0),r(iq1),r(iv),r(ie)) call leave return end c subroutine praxis(t0,machep,h0,n,x,f,fmin,scbd,ktm, 1 illc,d,y,z,q0,q1,v,e) real t0,machep,h0,x(n),f,fmin external f logical illc integer nl,nf,kl,kt,ktm real s,sl,dn,dmin,fx,f1,lds,ldt,t,h,sf,df,qf1,qd0,qd1,qa,qb,qc real m2,m4,small,vsmall,large,vlarge,scbd,ldfac,t2,dni,value c real d(n),y(n),z(n),q0(n),q1(n),v(n,n) real e(n) logical kenplt common /kenpt/ kenplt common /global/ fx,ldt,dmin,nf,nl . /q/ qa,qb,qc,qd0,qd1,qf1 idim=n kenplt=.false. c c.....initialization..... c machine dependent numbers[ c small=machep*machep vsmall=small*small large=1.0/small vlarge=1.0/vsmall m2=sqrt(machep) m4=sqrt(m2) c c ldfac=0.01 if (illc) ldfac=0.1 kt=0 nl=0 nf=1 fx=f(x,n) qf1=fx t=small+abs(t0) t2=t dmin=small h=h0 if (h.lt.100.0*t) h=100.0*t ldt=h c.....the first set of search directions v is the identity matrix..... do 20 i=1,n do 10 j=1,n 10 v(i,j)=0.0 20 v(i,i)=1.0 d(1)=0.0 qd0=0.0 do 30 i=1,n q0(i)=x(i) 30 q1(i)=x(i) c c.....the main loop starts here..... 40 sf=d(1) d(1)=0.0 s=0.0 c c.....minimize along the first direction v(*,1). c fx must be passed to min by value. value=fx call min(n,1,2,d(1),s,value,.false.,f,x,t,machep,h, 1v,q0,q1) if (s.gt.0.0) go to 50 do 45 i=1,n 45 v(i,1)=-v(i,1) 50 if (sf.gt.0.9*d(1).and.0.9*sf.lt.d(1)) go to 70 do 60 i=2,n 60 d(i)=0.0 c c.....the inner loop starts here..... 70 do 170 k=2,n do 75 i=1,n 75 y(i)=x(i) sf=fx if (kt.gt.0) illc=.true. 80 kl=k df=0.0 c c.....a random step follows (to avoid resolution valleys). c praxis assumes that random returns a random number uniformly c distributed in (0,1). c if(.not.illc) go to 95 do 90 i=1,n s=(0.1*ldt+t2*(10.0**kt))*(uni(0)-0.5) z(i)=s do 85 j=1,n 85 x(j)=x(j)+s*v(j,i) 90 continue fx=f(x,n) nf=nf+1 c c.....minimize along the ?non-conjugate? directions v(*,k),...,v(*,n) c 95 do 105 k2=k,n sl=fx s = 0.0 value=fx call min(n,k2,2,d(k2),s,value,.false.,f,x,t,machep,h, 1 v,q0,q1) write(6,900)(x(i),i=1,n) write(6,901)fx,nf 900 format(" current x",5e15.7) 901 format(" functions value is",e15.7," after evaluation",i5) if (illc) go to 97 s=sl-fx go to 99 97 s=d(k2)*((s+z(k2))**2) 99 if (df.gt.s) go to 105 df=s kl=k2 105 continue if (illc.or.(df.ge.abs((100.*machep)*fx))) go to 110 c c.....if there was not much improvement on the first try, set c illc=true and start the inner loop again..... c illc=.true. go to 80 110 continue c c.....minimize along the ?conjugate? directions v(*,1),...,v(*,k-1) c km1=k-1 do 120 k2=1,km1 s=0 value=fx call min(n,k2,2,d(k2),s,value,.false.,f,x,t,machep,h, 1 v,q0,q1) write(6,900)(x(i),i=1,n) write(6,901)fx,nf 120 continue f1=fx fx=sf lds=0 do 130 i=1,n sl=x(i) x(i)=y(i) sl=sl-y(i) y(i)=sl 130 lds=lds+sl*sl lds=sqrt(lds) if (lds.le.small) go to 160 c c.....discard direction v(*,kl). c if no random step was taken, v(*,kl) is the ?non-conjugate? c direction along which the greatest improvement was made..... c klmk=kl-k if (klmk.lt.1) go to 141 do 140 ii=1,klmk i=kl-ii do 135 j=1,n 135 v(j,i+1)=v(j,i) 140 d(i+1)=d(i) 141 d(k)=0 do 145 i=1,n 145 v(i,k)=y(i)/lds c c.....minimize along the new ?conjugate? direction v(*,k), which is c the normalized vector[ (new x) - (0ld x)..... c value=f1 call min(n,k,4,d(k),lds,value,.true.,f,x,t,machep,h, 1 v,q0,q1) write(6,902)(x(i),i=1,n) write(6,901)fx,nf 902 format(" x at end of inner loop",4e15.7) kenplt=.true. fff=f(x,n) kenplt=.false. if (lds.gt.0.0) go to 160 lds=-lds do 150 i=1,n 150 v(i,k)=-v(i,k) 160 ldt=ldfac*ldt if (ldt.lt.lds) ldt=lds t2 = 0.0 do 165 i=1,n 165 t2=t2+x(i)**2 t2=m2*sqrt(t2)+t c c.....see whether the length of the step taken since starting the c inner loop exceeds half the tolerance..... c if (ldt.gt.(0.5*t2)) kt=-1 kt=kt+1 if (kt.gt.ktm) go to 400 170 continue c.....the inner loop ends here. c c try quadratic extrapolation in case we are in a curved valley. c 171 call quad(n,f,x,t,machep,h,v,q0,q1) dn=0.0 do 175 i=1,n d(i)=1.0/sqrt(d(i)) if (dn.lt.d(i)) dn=d(i) 175 continue do 180 j=1,n s=d(j)/dn do 180 i=1,n 180 v(i,j)=s*v(i,j) c c.....scale the axes to try to reduce the condition number..... c if (scbd.le.1.0) go to 200 s=vlarge do 185 i=1,n sl=0.0 do 182 j=1,n 182 sl=sl+v(i,j)*v(i,j) z(i)=sqrt(sl) if (z(i).lt.m4) z(i)=m4 if (s.gt.z(i)) s=z(i) 185 continue do 195 i=1,n sl=s/z(i) z(i) = 1.0/sl if (z(i).le.scbd) go to 189 sl = 1.0/scbd z(i)=scbd 189 do 190 j=1,n 190 v(i,j)=sl*v(i,j) 195 continue c c.....calculate a new set of orthogonal directions before repeating c the main loop. c first transpose v for minfit[ c 200 do 220 i=2,n im1=i-1 do 210 j=1,im1 s=v(i,j) v(i,j)=v(j,i) 210 v(j,i)=s 220 continue c c.....call minfit to find the singular value decomposition of v. c this gives the principal values and principal directions of the c approximating quadratic form without squaring the condition c number..... c call minfit(idim,n,machep,vsmall,v,d,e) c c.....unscale the axes..... c if (scbd.le.1.0) go to 250 do 230 i=1,n s=z(i) do 225 j=1,n 225 v(i,j)=s*v(i,j) 230 continue do 245 i=1,n s=0.0 do 235 j=1,n 235 s=s+v(j,i)**2 s=sqrt(s) d(i)=s*d(i) s=1.0/s do 240 j=1,n 240 v(j,i)=s*v(j,i) 245 continue c c 250 do 270 i=1,n dni=dn*d(i) if (dni.gt.large) go to 265 if (dni.lt.small) go to 260 d(i)=1.0/(dni*dni) go to 270 260 d(i)=vlarge go to 270 265 d(i)=vsmall 270 continue c c.....sort the eigenvalues and eigenvectors..... c call sort(idim,n,d,v) dmin=d(n) if (dmin.lt.small) dmin=small illc=.false. if (m2*d(1).gt.dmin) illc=.true. c.....the main loop ends here..... c go to 40 c c.....return..... c 400 continue fmin=fx return end subroutine minfit(m,n,machep,tol,ab,q,e) real machep dimension ab(m,n),q(n) c...an improved version of minfit (see golub and reinsch, 1969) c restricted to m=n,p=0. c the singular values of the array ab are returned in q and ab is c overwritten with the orthogonal matrix v such that u.diag(q) = ab.v, c where u is another orthogonal matrix. real e(n) c...householder*s reduction to bidiagonal form... if (n.eq.1) go to 200 eps = machep g = 0.0 x = 0.0 do 11 i=1,n e(i) = g s = 0.0 l = i + 1 do 1 j=i,n 1 s = s + ab(j,i)**2 g = 0.0 if (s.lt.tol) go to 4 f = ab(i,i) g = sqrt(s) if (f.ge.0.0) g = -g h = f*g - s ab(i,i)=f-g if (l.gt.n) go to 4 do 3 j=l,n f = 0.0 do 2 k=i,n 2 f = f + ab(k,i)*ab(k,j) f = f/h do 3 k=i,n 3 ab(k,j) = ab(k,j) + f*ab(k,i) 4 q(i) = g s = 0.0 if (i.eq.n) go to 6 do 5 j=l,n 5 s = s + ab(i,j)*ab(i,j) 6 g = 0.0 if (s.lt.tol) go to 10 if (i.eq.n) go to 16 f = ab(i,i+1) 16 g = sqrt(s) if (f.ge.0.0) g = -g h = f*g - s if (i.eq.n) go to 10 ab(i,i+1) = f - g do 7 j=l,n 7 e(j) = ab(i,j)/h do 9 j=l,n s = 0.0 do 8 k=l,n 8 s = s + ab(j,k)*ab(i,k) do 9 k=l,n 9 ab(j,k) = ab(j,k) + s*e(k) 10 y = abs(q(i)) + abs(e(i)) 11 if (y.gt.x) x = y c...accumulation of right-hand transformations... ab(n,n) = 1.0 g = e(n) l = n do 25 ii=2,n i = n - ii + 1 if (g.eq.0.0) go to 23 h = ab(i,i+1)*g do 20 j=l,n 20 ab(j,i) = ab(i,j)/h do 22 j=l,n s = 0.0 do 21 k=l,n 21 s = s + ab(i,k)*ab(k,j) do 22 k=l,n 22 ab(k,j) = ab(k,j) + s*ab(k,i) 23 do 24 j=l,n ab(i,j) = 0.0 24 ab(j,i) = 0.0 ab(i,i) = 1.0 g = e(i) 25 l = i c...diagonalization of the bidiagonal form... 100 eps = eps*x do 150 kk=1,n k = n - kk + 1 kt = 0 101 kt = kt + 1 if (kt.le.30) go to 102 e(k) = 0.0 call seterr(22hfminn-problem with svd,22,2,2) return 102 do 103 ll2=1,k l2 = k - ll2 + 1 l = l2 if (abs(e(l)).le.eps) go to 120 if (l.eq.1) go to 103 if (abs(q(l-1)).le.eps) go to 110 103 continue c...cancellation of e(l) if l'1... 110 c = 0.0 s = 1.0 do 116 i=l,k f = s*e(i) e(i) = c*e(i) if (abs(f).le.eps) go to 120 g = q(i) c...q(i) = h = dsqrt(g*g + f*f)... if (abs(f).lt.abs(g)) go to 113 if (f) 112,111,112 111 h = 0.0 go to 114 112 h = abs(f)*sqrt(1.0 + (g/f)**2) go to 114 113 h = abs(g)*sqrt(1.0+ (f/g)**2) 114 q(i) = h if (h.ne.0.0) go to 115 g = 1.0 h = 1.0 115 c = g/h 116 s = -f/h c...test for convergence... 120 z = q(k) if (l.eq.k) go to 140 c...shift from bottom 2*2 minor... x = q(l) y = q(k-1) g = e(k-1) h = e(k) f = ((y - z)*(y + z) + (g - h)*(g + h))/(2.0*h*y) g = sqrt(f*f + 1.0) temp = f - g if (f.ge.0.0) temp = f + g f = ((x - z)*(x + z) + h*(y/temp - h))/x c...next qr transformation... c = 1.0 s = 1.0 lp1 = l + 1 if (lp1.gt.k) go to 133 do 132 i=lp1,k g = e(i) y = q(i) h = s*g g = g*c if (abs(f).lt.abs(h)) go to 123 if (f) 122,121,122 121 z = 0.0 go to 124 122 z = abs(f)*sqrt(1.0+ (h/f)**2) go to 124 123 z = abs(h)*sqrt(1.0+ (f/h)**2) 124 e(i-1) = z if (z.ne.0.0) go to 125 f = 1.0 z = 1.0 125 c = f/z s = h/z f = x*c + g*s g = -x*s + g*c h = y*s y = y*c do 126 j=1,n x = ab(j,i-1) z = ab(j,i) ab(j,i-1) = x*c + z*s 126 ab(j,i) = -x*s + z*c if (abs(f).lt.abs(h)) go to 129 if (f) 128,127,128 127 z = 0.0 go to 130 128 z = abs(f)*sqrt(1.0+ (h/f)**2) go to 130 129 z = abs(h)*sqrt(1.0+ (f/h)**2) 130 q(i-1) = z if (z.ne.0.0) go to 131 f = 1.0 z = 1.0 131 c = f/z s = h/z f = c*g + s*y 132 x = -s*g + c*y 133 e(l) = 0.0 e(k) = f q(k) = x go to 101 c...convergence[ q(k) is made non-negative... 140 if (z.ge.0.0) go to 150 q(k) = -z do 141 j=1,n 141 ab(j,k) = -ab(j,k) 150 continue return 200 q(1) = ab(1,1) ab(1,1) = 1.0 return end subroutine min(n,j,nits,d2,x1,f1,fk,f,x,t,machep,h, 1v,q0,q1) external f logical fk real machep,x(n),ldt dimension v(n,n),q0(n),q1(n) common /global/ fx,ldt,dmin,nf,nl . /q/ qa,qb,qc,qd0,qd1,qf1 c...the subroutine min minimizes f from x in the direction v(*,j) unless c j is less than 1, when a quadratic search is made in the plane c defined by q0,q1,x. c d2 is either zero or an approximation to half f?. c on entry, x1 is an estimate of the distance from x to the minimum c along v(*,j) (or, if j=0, a curve). on return, x1 is the distance c found. c if fk=.true., then f1 is flin(x1). otherwise x1 and f1 are ignored c on entry unless final fx is greater than f1. c nits controls the number of times an attempt will be made to halve c the interval. logical dz real m2,m4 small = machep**2 m2=sqrt(machep) m4 = sqrt(m2) sf1 = f1 sx1 = x1 k = 0 xm = 0.0 fm = fx f0 = fx dz = d2.lt.machep c...find the step size... s = 0.0 do 1 i=1,n 1 s = s + x(i)**2 s = sqrt(s) temp = d2 if (dz) temp = dmin t2 = m4*sqrt(abs(fx)/temp + s*ldt) + m2*ldt s = m4*s + t if (dz.and.t2.gt.s) t2 = s t2=amax1(t2,small) t2=amin1(t2,.01*h) if (.not.fk.or.f1.gt.fm) go to 2 xm = x1 fm = f1 2 if (fk.and.abs(x1).ge.t2) go to 3 temp = 1.0 if (x1.lt.0.0) temp = -1.0 x1=temp*t2 f1 = flin(n,j,x1,f,x,nf,v,q0,q1) 3 if (f1.gt.fm) go to 4 xm = x1 fm = f1 4 if (.not.dz) go to 6 c...evaluate flin at another point and estimate the second derivative... x2 = -x1 if (f0.ge.f1) x2 = 2.0*x1 f2 = flin(n,j,x2,f,x,nf,v,q0,q1) if (f2.gt.fm) go to 5 xm = x2 fm = f2 5 d2 = (x2*(f1 - f0)-x1*(f2 - f0))/((x1*x2)*(x1 - x2)) c...estimate the first derivative at 0... 6 d1 = (f1 - f0)/x1 - x1*d2 dz = .true. c...predict the minimum... if (d2.gt.small) go to 7 x2 = h if (d1.ge.0.0) x2 = -x2 go to 8 7 x2 = (-0.5*d1)/d2 8 if (abs(x2).le.h) go to 11 if (x2) 9,9,10 9 x2 = -h go to 11 10 x2 = h c...evaluate f at the predicted minimum... 11 f2 = flin(n,j,x2,f,x,nf,v,q0,q1) if (k.ge.nits.or.f2.le.f0) go to 12 c...no success, so try again... k = k + 1 if (f0.lt.f1.and.(x1*x2).gt.0.0) go to 4 x2 = 0.5*x2 go to 11 c...increment the one-dimensional search counter... 12 nl = nl + 1 if (f2.le.fm) go to 13 x2 = xm go to 14 13 fm = f2 c...get a new estimate of the second derivative... 14 if (abs(x2*(x2-x1)).le.small) go to 15 d2 = (x2*(f1-f0) - x1*(fm-f0))/((x1*x2)*(x1 - x2)) go to 16 15 if (k.gt.0) d2 = 0.0 16 if (d2.le.small) d2 = small x1 = x2 fx = fm if (sf1.ge.fx) go to 17 fx = sf1 x1 = sx1 c...update x for linear but not parabolic search... 17 if (j.eq.0) return do 18 i=1,n 18 x(i) = x(i) + x1*v(i,j) return end real function flin(n,j,l,f,x,nf,v,q0,q1) real l,x(n) dimension v(n,n),q0(n),q1(n) c...flin is the function of one real variable l that is minimized c by the subroutine min... common /q/ qa,qb,qc,qd0,qd1,qf1 common /cstack/ d(500) double precision d(500) real t(1000) equivalence(d(1),t(1)) call enter(1) it =istkgt(n,3) if (j .eq. 0) go to 2 c...the search is linear... it1=it-1 do 1 i=1,n itt=it1+i 1 t(itt) = x(i) + l*v(i,j) go to 4 c...the search is along a parabolic space curve... 2 qa = (l*(l - qd1))/(qd0*(qd0 + qd1)) qb = ((l + qd0)*(qd1 - l))/(qd0*qd1) qc = (l*(l + qd0))/(qd1*(qd0 + qd1)) it1=it-1 do 3 i=1,n itt=it1+i 3 t(itt) = (qa*q0(i) + qb*x(i)) + qc*q1(i) c...the function evaluation counter nf is incremented... 4 nf = nf + 1 flin = f(t(it),n) call leave return end subroutine sort(m,n,d,v) real d(n),v(m,n) c...sorts the elements of d(n) into descending order and moves the c corresponding columns of v(n,n). c m is the row dimension of v as declared in the calling program. real s if (n.eq.1) return nm1 = n - 1 do 3 i = 1,nm1 k=i s = d(i) ip1 = i + 1 do 1 j = ip1,n if (d(j) .le. s) go to 1 k = j s = d(j) 1 continue if (k .le. i) go to 3 d(k) = d(i) d(i) = s do 2 j = 1,n s = v(j,i) v(j,i) = v(j,k) 2 v(j,k) = s 3 continue return end subroutine quad(n,f,x,t,machep,h,v,q0,q1) external f c...quad looks for the minimum of f along a curve defined by q0,q1,x... real x(n),machep,ldt,l dimension v(n,n),q0(n),q1(n) common /global/ fx,ldt,dmin,nf,nl . /q/ qa,qb,qc,qd0,qd1,qf1 s = fx fx = qf1 qf1 = s qd1 = 0.0 do 1 i=1,n s = x(i) l = q1(i) x(i) = l q1(i) = s 1 qd1 = qd1+ (s-l)**2 qd1 = sqrt(qd1) l = qd1 s = 0.0 if (qd0.le.0.0.or.qd1.le.0.0.or.nl.lt.3*n*n) go to 2 value=qf1 call min(n,0,2,s,l,value,.true.,f,x,t,machep,h,v,q0,q1) qa = (l*(l-qd1))/(qd0*(qd0+qd1)) qb = ((l+qd0)*(qd1-l))/(qd0*qd1) qc = (l*(l+qd0))/(qd1*(qd0+qd1)) go to 3 2 fx = qf1 qa = 0.0 qb = qa qc = 1.0 3 qd0 = qd1 do 4 i=1,n s = q0(i) q0(i) = x(i) 4 x(i) = (qa*s + qb*x(i)) + qc*q1(i) return end