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