Three dimensional max-min optimization

c three dimensional max-min  optimization
c See Numerical Analysis Burden & Faires Chapter 9
c exercise set 9.4
c the code presented is used to minimize
c i.e.  sm= sm=amin1(f1,f2,f3,f4,f5,f6)
c to maximize use sm=amax1(f1,f2,f3,f4,f5,f6)
      equivalence(dx,dy,dz,delta),(x,x1) ,(y,x2),(z,x3)
      data nstep ,delta, epsi,tol/ 1000, .01 ,1.e-4, 1.e-2/
      data x0,y0,z0/0.,0.,0./
      eq1(x,y,z)=x+cos(x*y*z)-1.
      eq2(x,y,z)=  (1.-x)**.25+y+.05*z**2-.15*z-1.
      eq3(x,y,z) = -x**2  -.1*y**2+.01*y+z-1.
      f(x,y,z)= eq1(x,y,z)**2  +eq2(x,y,z)**2 + eq3(x,y,z)**2
c      f(x,y,z)=(x-1.)**2 +(y-2.)**2 +(z-3.)**2
      dfdx1(x1,x2,x3)=(f(x1+epsi,x2,x3)-f(x1,x2,x3))/epsi
      dfdx2(x1,x2,x3)=(f(x1,x2+epsi,x3)-f(x1,x2,x3))/epsi
      dfdx3(x1,x2,x3)=(f(x1,x2,x3+epsi)-f(x1,x2,x3))/epsi
      pi=2.*asin(1.)
      kp=int(float(nstep)/20.)
      kount=kp
      s= f(x0,y0,z0)
      do 10 i=1,nstep
      f1=f(x0+dx,y0,z0)
      f2=f(x0-dx,y0,z0)
      f3=f(x0,y0+dy,z0)
      f4=f(x0,y0-dy,z0)
      f5=f(x0,y0,z0+dz)
      f6=f(x0,y0,z0-dz)
      sm=amin1(f1,f2,f3,f4,f5,f6)
      if(sm.eq.f1)then
      x0=x0+dx
      endif
      if(sm.eq.f2)then
      x0=x0-dx
      endif
      if(sm.eq.f3)then
      y0=y0+dy
      endif
      if(sm.eq.f4)then
      y0=y0-dy
      endif
      if(sm.eq.f5)then
      z0=z0+dz
      endif
      if(sm.eq.f6)then
      z0=z0-dz
      endif
      if(i.eq.kount)then
      print 100,x0,y0,z0,f(x0,y0,z0)
      print*,'  '
      kount=kount+kp
      endif
      if(abs(dfdx1(x0,y0,z0)).le.tol.and.abs(dfdx2(x0,y0,z0)).le.tol
     $ .and.abs(dfdx3(x0,y0,z0)).le.tol)goto 300
      if(i.gt.int(.5*float(nstep)))delta=.002
      if(i.gt.int(.8*float(nstep)))delta=.0005
      s=f(x0,y0,z0)
10    continue
100   format(1x,'x0,y0,z0,f=',3(3x,e13.6),2x,e11.4)
300   print*,'number of iterations=',i
      print 100,x0,y0,z0,f(x0,y0,z0)
      print*,'dfdx1,dfdx2,dfdx3=',dfdx1(x0,y0,z0), dfdx2(x0,y0,z0),
     $ dfdx3(x0,y0,z0)
      stop
      end

      RUN
      The initial value is (0., 0., 0.5)
      the root value is  (0., .1, 1.)
       x0,y0,z0,f=    0.000000E+00    0.000000E+00    0.500000E+00   0.2539E+00

 x0,y0,z0,f=    0.000000E+00    0.500000E-01    0.949999E+00   0.4720E-02

 number of iterations= 110
 x0,y0,z0,f=    0.000000E+00    0.100000E+00    0.999999E+00   0.4300E-12
 dfdx1,dfdx2,dfdx3=  0.000106230233  0.000100054858  9.88517713E-005