      subroutine setup2
      implicit none
      integer nmax,i,n,case,iterma
      parameter(nmax=300)
      real*8 pi,epsx,epsg,epsf,epsmach,tolmach,deldif,xst
      common/const/pi
      logical test1,test2,test3,test4,analyt
      common/gradinf/analyt
      common/param/epsx,epsg,epsf,test1,test2,test3,test4
      common/mach/epsmach,tolmach,deldif
      common/dim/n
      common/fdat/xst(nmax),case
      character*30 xifile,xofile
      character*80 ident 
      common/ident/ident,xifile,xofile
      real*8 dt
      real*8 kl,ya,ye,rho
      common/kette/rho,kl,ya,ye
      real*8 idt
      integer niter
      real*8 eta,thet2,bfgsterm
      common/bfgsparam/niter,iterma,eta,thet2,bfgsterm
      iterma=100
      pi=datan(1.d17)*2.d0
      call setda1
      call setda2
      test1=.false.
      test2=.false.
      test3=.false.
      test4=.false.
      goto (100,200,300,400,500,600,700,800,900,1000,1100,1200,
     f      1300,1400,1500,1600,1700,1800,1900,2000,2100,2200,
     f      2300,2400,2500,2600,2700,2800,2900,3000,3100,3200,
     f      3300,3400,3500,3600,3700,3800,3900,4000,4100,
     f      4200,4300) case
c*******************************
  100 continue
      ident='rosenbrockbanana valley'
      test1=.true.
      test2=.true.
      test3=.true.
      n=2
      xst(1)=-1.2d0
      xst(2)=1.d0
      return
c*******************************
  200 continue
      ident='woodsquartic function'
      n=4
      xst(1)=-3.d0
      xst(2)=-1.d0
      xst(3)=-3.d0
      xst(4)=-1.d0
      return
c*******************************
  300 continue
      ident='bealesfunction'
      n=2
      xst(1)=5.d0
      xst(2)=5.d0
      return
c*******************************
  400 continue
      ident='quadrfunktion n=9 x*=0'
      n=9
      do i=1,n
        xst(i)=dble(float(i))
      enddo
      return
c*******************************
  500 continue
      ident='powellsquartic'
      epsx=1.d-5
      epsg=1.d-5
      n=4
      xst(1)=3.d0
      xst(2)=-1.d0
      xst(3)=0.d0
      xst(4)=1.d0            
      return
c*******************************
  600 continue
      ident='boxexponential, x*=(1,10,1)'
      xst(1)=0.d0
      xst(2)=20.d0
      xst(3)=1.d0
      n=3
      return
c*******************************
  700 continue
      ident='fletcherhelical function'
      n=3
      xst(1)=-1.d0
      xst(2)=0.d0
      xst(3)=0.d0
      return
c*******************************
  800 continue
      ident='powells second function, continuum of minima'
      epsx=1.d-5
      epsg=1.d-5
      n=3
      xst(1)=0.d0
      xst(2)=1.d0
      xst(3)=2.d0
      return
c*******************************
  900 continue
      ident='himmelblau 4 minima 4 sattelpunkte'
      n=2
      xst(1)=0.d0
      xst(2)=-2.d0
      return
c*******************************
 1000 continue
      ident='bsm8.5 betts modified rosenbrock'
      n=2
      xst(1)=-2.d0
      xst(2)=-2.d0
      return
c*******************************
 1100 continue
      ident='bsm8.6 betts modified rosenbrock '
      n=2
      xst(1)=1.d0
      xst(2)=5.d0
      return
c*******************************
 1200 continue
      ident='bsm8.7 betts modified rosenrock '
      n=2
      xst(1)=-1.2d0
      xst(2)=1.d0
      return
c*******************************
 1300 continue
      ident='pavianilikelihoodfunction'
      deldif=1.d-7
      n=10
      do i=1,n
        xst(i)=9.d0
      enddo
      return
c*******************************
 1400 continue
      ident='polynomleastsquares '
      n=6
      do i=1,n
        xst(i)=0.d0
      enddo
      epsg=1.d-18
      iterma=300
      return
c*******************************
 1500 continue
      ident='nonlinearlsq  least squares'
      n=2
      xst(1)=3.01d0
      xst(2)=-.51d0
      return
c*******************************
 1600 continue
      ident='shanno tridia'
      iterma=400
      n=50
      do i=1,n
       xst(i)=-1.d0
      enddo      
      return
c*******************************
 1700 continue
      ident='genrose1 function'
      iterma=500
      n=50
      do i=1,n
       xst(i)=-1.d0
      enddo
      return
c*******************************
 1800 continue
      ident='hs70himmelblau'
      deldif=1.d-7
      n=4
      xst(1)=2.d0
      xst(2)=4.d0
      xst(3)=4.d-2
      xst(4)=2.d0
      return
c*******************************
 1900 continue
      ident='gen2ros'
      dt=1.d0/dble(float(n+1))
      iterma=2000
      n=100
      do i=1,n
        xst(i)=dble(float(i))*dt
      enddo
      epsx=1.d-7
      epsf=1.d-9
      epsg=1.d-7
      return
c*******************************
 2000 continue
      ident='shubertmultimin '
      n=2
      xst(1)=1.d0
      xst(2)=1.d0
      return
c*******************************
 2100 continue
      ident='hs25holzmann'
      test2=.true.
      n=3
      xst(1)=55.d0
      xst(2)=12.5d0
      xst(3)=2.6d0
      return
c*******************************
 2200 continue
      ident='mccormickhs5 '
      n=2
      xst(1)=0.d0
      xst(2)=0.d0
      return
c*******************************
 2300 continue
      ident='collanihs68'
      deldif=1.d-7 
      n=2
      xst(1)=1.d0
      xst(2)=1.d0
      return
c*******************************
 2400 continue
      ident='collanihs69'
      deldif=1.d-7
      n=2
      xst(1)=1.d0
      xst(2)=1.d0
      return
c*******************************
 2500 continue
      ident='tridiaquad'
      iterma=1000
      n=50
      do i=1,n
        xst(i)=0.d0
      enddo
      return
c*******************************      
 2600 continue
      ident='kowalikosborne nlsq'
      n=4
      xst(1)=.25d0
      xst(2)=.39d0
      xst(3)=.415d0
      xst(4)=.39d0
      return
c*******************************
 2700 continue
      ident='bardnlsq70'
      n=3
      xst(1)=1.d0
      xst(2)=1.d0
      xst(3)=1.d0
      return
c*******************************
 2800 continue
      ident='browndennisnlsq'
      n=4
      xst(1)=1.d0
      xst(2)=(dexp(-4.d0)-1.d0)/4.d0
      xst(3)=0.d0
      xst(4)=0.d0
      return
c*******************************
 2900 continue
      ident='cal1fct nash bsp 1'
      iterma=1000
      n=50
      epsx=1.d-5
      epsf=1.d-9
      epsg=1.d-5
      do i=1,n
        xst(i)=0.d0
      enddo
      return
c*******************************
 3000 continue
      ident='cal2fct nash bsp 2'
      n=50
      iterma=2000
      epsx=1.d-5
      epsf=1.d-9
      epsg=1.d-5
      do i=1,n
        xst(i)=0.d0
      enddo
      return
c*******************************
 3100 continue
      ident='cal3fct nash bsp 3'
      n=50
      iterma=2000
      epsx=1.d-5
      epsf=1.d-9
      epsg=1.d-5
      do i=1,n
        xst(i)=0.d0
      enddo
      return
c*******************************
 3200 continue
      ident='kettenlinie'
      iterma=2000
      n=100
      kl=4.d0
      rho=100.d0
      ya=1.d0
      ye=3.d0
      dt=1.d0/dble(float(n+1))
c**** parabel-interpolation zu (0,ya),(.5,ye-kl),(1.,ye)
      do i=1,n
        idt=i*dt
        xst(i)=ya+(ye-ya-kl)*2.d0*idt
     f           +2.d0*(2.d0*kl-(ye-ya))*idt*(idt-.5d0)
      enddo
      epsx=1.d-9
      epsf=1.d-12
      epsg=1.d-9
      deldif=1.d-4
      return
c**********************************************
 3300 continue
      n=20
      iterma=300
      ident='chebyquad'
      do i=1,n
        xst(i)=dble(float(i))/dble(float(n+1))
      enddo
      return
 3400 continue
      ident='engva1fct(engvall)'
      n=2
      xst(1)=.5d0
      xst(2)=2.d0
      return
c************************************************
 3500 continue
      ident='engva2fct(engvall)'
      n=3
      xst(1)=1.d0
      xst(2)=2.d0
      xst(3)=0.d0
      return
c***********************************************
 3600 continue
      ident='cragglevy'
      n=4
      xst(1)=1.d0
      do i=2,4
        xst(i)=2.d0
      enddo
      return
c***********************************************
 3700 continue
      ident='charosfct'
      n=25
      do i=1,n
         xst(i)=-1.d0
      enddo
      return
c***********************************************
 3800 continue
      ident='pauly1fct'
      n=8
      do i=1,4
        xst(2*i-1)=.1d0
        xst(2*i)=.9d0
      enddo
      return
c***********************************************
 3900 continue
      ident='pauly2fct'
      n=30
      do i=1,15
        xst(2*i-1)=2.d0
        xst(2*i)=1.d0
      enddo
      return
c**********************************************
 4000 continue
      ident='mittelmann'
      iterma=2000
      n=100
      do i=1,n
        xst(i)=1/dble(float(i))
      enddo
      return
c***********************************************
 4100 continue
      ident='scalbeale mit x1-> x1/100'
      xst(1)=500.d0
      xst(2)=5.d0
      n=2
      test2=.true.
      return
c*******************************
 4200 continue
      ident='bsp3_1_6'
      n=2
      xst(1)=0.d0
      xst(2)=0.d0
      return
 4300 continue
      ident='modrosenbrock 10**4 statt 10**2'
      n=2
      xst(1)=-1.2d0
      xst(2)=1.d0
      iterma=2000
      return
      end
      subroutine f(x,fx,err)
      implicit none
      integer nmax,i,cf,cg,cu,csu,cre,cm,j,k,nn
      parameter (nmax=300)
      real*8 x(nmax),fx
      logical err
      integer n
      common/dim/n
      integer case
      real*8 xst
      common/fdat/xst(nmax),case
      character*80 ident 
      character*30 xifile,xofile
      common/ident/ident,xifile,xofile
      common/count/cf,cg,cu,csu,cre,cm
      real*8 x1,x2,x3,x4,f1,f2,f3,f4,f5,f6,x1h,x2h,t,thet
      real*8 vix,ui,arg,yint,yint1,yint2
      real*8 pi
      common/const/pi
      real*8 cb(3),ytrue,yest,sum,prod,u(13),v(13),w(13),y(13)
      real*8  xt(0:nmax-1),dxt(0:nmax-1),h,wa(0:nmax-1),a,b
      common /bspline/ xt,dxt,h,a,b
      real*8 d,c,cl,y18
      common /dat/ d(19),c(19),cl(19),y18(19)
      real*8 data,datb
      common /dat2/ data(99),datb(99)
      real*8 y26(11),u26(11),y27(15)
      real*8 kl,ya,ye,rho
      common/kette/rho,kl,ya,ye
      real*8 skal1,nue,ghn,ifx,ihx,ifgrad(nmax),ihgrad(nmax),e
      real*8 t0,t1,t2,sump(nmax),gam1,gam2,ksi1,ksi2,delcha(25)
      real*8 s,bet,betx4,w1,w2,dlbet,dlx4,fac1,fac2,yic,e1,e2,cl1,cl2
      real*8 norint
      data u/4*0.d0,4*1.d0,3*2.d0,2.5d0,2.9d0/
      data v/0.d0,1.d0,2.d0,3.d0,0.d0,1.d0,2.d0,2.d0,
     f       0.d0,1.d0,2.d0,2.d0,1.8d0/
      data w/4*1.d2,4*5.d1,3*2.5d1,2*1.d1/
      data y/2.93d0,1.95d0,.81d0,.58d0,5.9d0,4.74d0,4.18d0,4.05d0,
     f       9.03d0,7.85d0,7.22d0,8.5d0,9.81d0/
      data y26/.1957d0,.1947d0,.1735d0,.16d0,.0844d0,.0627d0,
     f         .0456d0,.0342d0,.0323d0,.0235d0,.0246d0/
      data u26/4.d0,2.d0,1.d0,.5d0,.25d0,.167d0,.125d0,.1d0,.0823d0,
     f         .0714d0,.0625d0/
      data y27/.14d0,.18d0,.22d0,.25d0,.29d0,.32d0,.35d0,.39d0,
     f         .37d0,.58d0,.73d0,.96d0,1.34d0,2.1d0,4.39d0/
      data delcha/1.25d0,1.4d0,2.4d0,1.4d0,1.75d0,1.2d0,2.25d0,
     f            1.2d0,1.0d0,1.1d0,1.5d0,1.6d0,1.25d0,
     f            1.25d0,1.2d0,1.2d0,1.4d0,.5d0,.5d0,1.25d0,1.8d0,
     f            0.75d0,1.25d0,1.4d0,1.6d0/
c******
      err=.false.
      cf=cf+1
      goto (100,200,300,400,500,600,700,800,900,1000,1100,1200,
     f      1300,1400,1500,1600,1700,1800,1900,2000,2100,2200,
     f      2300,2400,2500,2600,2700,2800,2900,3000,3100,3200,
     f      3300,3400,3500,3600,3700,3800,3900,4000,4100,4200
     f       ,4300) case
c*******************************
c******** rosenbrock      
  100 continue
      x1=x(1)
      x2=x(2)
      f1=x2-x1**2
      f2=x1-1.d0
      fx=100.d0*f1**2+f2**2
      return
c******** wood
  200 continue
      x1=x(1)
      x2=x(2)
      x3=x(3)
      x4=x(4)
      f1=x2-x1**2
      f2=x1-1.d0
      f3=x4-x3**2
      f4=x3-1.d0
      f5=x2-1.d0
      f6=x4-1.d0
      fx=1.d2*f1**2+f2**2+9.d1*f3**2+f4**2+10.1d0*(f5**2+f6**2)+
     f   19.8d0*f5*f6
      return
c****************** beale's function
  300 continue
      x1=x(1)
      x2=x(2)
      cb(1)=1.5d0
      cb(2)=2.25d0
      cb(3)=2.625d0
      fx=0.d0
      do i=1,3
        fx=(cb(i)-x1*(1.d0-x2**i))**2+fx
      enddo
      return
c**************** nazareth's quadratic
  400 continue
      fx=0.d0
      do i=1,n
        fx=fx+dble(float(i))*x(i)**2
      enddo
      return
c*************** powell's quartic
  500 continue
      x1=x(1)
      x2=x(2)
      x3=x(3)
      x4=x(4)
      f1=x1+10.d0*x2
      f2=x3-x4
      f3=x2-2.d0*x3
      f4=x1-x4
      fx=f1**2+5.d0*f2**2+f3**4+10.d0*f4**4
      return
c******************* box exponential
  600 continue
      x1=x(1)
      x2=x(2)
      x3=x(3)
      fx=0.d0
      do i=1,10
        f1=dble(float(i))
        f2=f1*.1d0
        fx=fx
     f   +((dexp(-x1*f2)-dexp(-x2*f2))-x3*(dexp(-f2)-dexp(-f1)))**2    
      enddo
      return
c*************************** fletcher's helical
  700 continue
      x1=x(1)
      x2=x(2)
      x3=x(3)
      x1h=dabs(x1)
      x2h=dabs(x2)
      t=x1h+x2h
      if ( t .eq. x2h ) then
       if ( x1*x2 .lt. 0.d0 ) then
         thet=-pi*.5d0
       else
         thet=pi*.5d0
       endif
      else
         thet=datan(x2/x1)
      endif
      thet=thet/(2.d0*pi)
      if ( x1 .eq. 0.d0 ) then
        thet=.25d0
      else
        if ( x1 .lt. 0.d0 )  thet=thet+.5d0
      endif
      fx=1.d2*((x3-1.d1*thet)**2+(dsqrt(x1**2+x2**2)-1.d0)**2)+x3**2
      return
c********************* powells second function
  800 continue
      x1=x(1)
      x2=x(2)
      x3=x(3)
      f1=1.d0+(x1-x2)**2
      f2=pi*x2*x3/2.d0
      f3=0.d0
      if ( dabs(x1) .le. 9.d0*dabs(x2) ) f3=dexp(-(x1/x2-1.d0)**2)    
      fx=3.d0-1.d0/f1-dsin(f2)-f3
      return
c******************** himmelblau 4 minima 4 sattelpunkte
  900 continue
      fx=(x(1)**2+x(2)-11.d0)**2+(x(1)+x(2)**2-7.d0)**2
      return
c******************** betts modified rosenbrock 8.5
 1000 continue
      x1=x(1)
      x2=x(2)
      f1=x2-x1**2 
      f2=1.d0-x1
      fx=f1**2+f2**2
      return
c******************** betts modified rosenbrock 8.6
 1100 continue
      x1=x(1)
      x2=x(2)
      f1=x2-x1**2
      f2=1.d0-x1
      fx=f1**2+1.d2*f2**2
      return
c******************** betts modified rosenbrock 8.7
 1200 continue
      x1=x(1)
      x2=x(2)
      f1=x2-x1**3
      f2=1.d0-x1
      fx=1.d2*f1**2+f2**2
      return
c******************** paviani maximum likelihood
 1300 continue
      sum=0.d0
      prod=1.d0
      do i=1,n
        if ( x(i).le.2.d0 .or. x(i) .ge. 10.d0 ) then
          fx=sum+prod
          err=.true.
          return
        endif
        prod=prod*x(i)
        sum=sum+dlog(x(i)-2.d0)**2+dlog(10.d0-x(i))**2
      enddo
      fx=sum-dexp(.2d0*dlog(prod))
      return
c******************** polynomial least squares
 1400 continue
      sum=0.d0
      do i=1,21
        t=dble(float(i))-1.d0
        ytrue=((((t+1.d0)*t+1.d0)*t+1.d0)*t+1.d0)*t+1.d0
        yest=((((x(6)*t+x(5))*t+x(4))*t+x(3))*t+x(2))*t+x(1)
        sum=(ytrue-yest)**2+sum
      enddo
      fx=sum
      return
c************************ nonlinear least squares
 1500 continue
      x1=x(1)
      x2=x(2)
      sum=0.d0
      do i=1,13
        sum=sum+w(i)*(y(i)-(x1*u(i)+dexp(x2*v(i))))**2
      enddo
      fx=sum
      return
c***************************** shanno tridia
 1600 continue
      sum=(1.d0-x(1))**2
      do i=2,n
        sum=sum+1.d2*(x(i)-x(i-1))**2
      enddo
      fx=sum
      return
c***************************** genrose
 1700 continue
      sum=(1.d0-x(1))**2
      do i=2,n
        sum=sum+(1.d0-x(i))**2+1.d2*(x(1)-x(i)**2)**2
      enddo
      fx=sum
      return
c***************************** himmelblau hs70
 1800 continue
      s=0.d0
      x1=x(1)
      x2=x(2)
      x3=x(3)
      x4=x(4)
      if ( x1 .le. 1.d-4 .or. x2 .le. 1.d-4 .or. x3 .le. 1.d-4
     f     .or. x4 .le. 1.d-4 .or. x3+(1.d0-x3)*x4 .le. 1.d-4 )
     f    then
          fx=0.d0
          err=.true.
          return
      endif
      bet=x3+(1.d0-x3)*x4
      w1 =dsqrt(dabs(x1))
      w2 =dsqrt(dabs(x2))
      betx4=0.d0
      dlx4=0.d0
      dlbet=0.d0
      dlbet=dlog(bet)
      betx4=bet/x4
      dlx4=dlog(x4)
      fac1=w2*x2*x3/(1.d0+12.d0*x2)
      fac2=w1*x1*(1.d0-x3)/(1.d0+12.d0*x1)
      do   i=1,19
       yic=0.d0
       e1=0.d0
       cl1=x2*(1.d0-c(i)*bet  +dlbet     +cl(i))
       if(cl1 .ge. -80.d0)   e1=dexp(cl1)
       e2=0.d0
       cl2=x1*(1.d0-c(i)*betx4+dlbet-dlx4+cl(i))
       if(cl2 .ge. -80.d0)   e2=dexp(cl2)
       yic=d(i)*(fac1*e1+fac2*e2)
       s=s+(yic-y18(i))**2
      enddo
      fx=s-1.d-3*(dlog(x1-1.d-4)+dlog(x2-1.d-4)+dlog(x3-1.d-4)
     f          +dlog(x4-1.d-4)+dlog(bet-1.d-4))
      return
c******************************* genrose2
 1900 continue
      fx=1.d0
      do  i=2,n
	 fx=fx+100.d0*((x(i)-x(i-1)**2)**2)+(1.d0-x(i))**2
      enddo
      return
c******************************* shubert multimin
 2000 continue
	fx = -(dcos(2.*x(1)+1.)+2.*dcos(3.*x(1)+2.)+3.*dcos(4.*x(1)+
     a	3.)+4.*dcos(5.*x(1)+4.)+5.*dcos(6.*x(1)+5.))*(dcos(2.*
     a	x(2)+1.)+2.*dcos(3.*x(2)+2.)+3.*dcos(4.*x(2)+3.)+4.*dcos(5.*
     a	x(2)+4.)+5.*dcos(6.*x(2)+5.))
      return
c******************************* holzmann
 2100 continue
c      if ( x(3) .lt. 1.d-2 .or. x(1) .lt. 1.d-2 
c     f     .or. x(2) .ge. 25.6d0 ) then
c        fx=0.d0
c        err=.true.
c        return
c      endif
      x(3)=max(1.d-2,x(3))
      x(1)=max(1.d-2,x(1))
      x(2)=min(25.6,x(2))
      fx=0.d0
      do   i=1,99
        ui=data(i)
        vix=0.d0
        arg=-(ui-x(2))**x(3)/x(1)
        if(arg .gt. -80.d0)    vix=dexp(arg)
        f1=datb(i)+vix
        fx=fx+f1**2
      enddo
      return
c******************************** hs5 mccormick
 2200 continue
      fx=dsin(x(1)+x(2))+(x(1)-x(2))**2-1.5d0*x(1)+2.5d0*x(2)+1.d0
      return
c******************************** collani 1 hs68
 2300 continue
      if ( x(1) .le. 1.d-4 ) then
       x(1)=1.d-4
      endif
      e=dexp(x(1))-1.d0
      f1=dsqrt(24.d0)
      yint=norint(-x(2))
c imsld      call mdnord(-x(2),yint)
      x3=2.d0*yint
      yint1=norint(-x(2)+f1)
      yint2=norint(-x(2)-f1)
c imsld      call mdnord(-x(2)+f1,yint1)
c imsld      call mdnord(-x(2)-f1,yint2)
      x4=yint1+yint2
      fx=(2.4d-3-(e-x3)*x4/(e+x4))/x(1)
     f     -1.d-1*dlog(x(1)-1.d-5)
      return
c******************************* collani 2 hs 69
 2400 continue
      if ( x(1) .le. 1.d-4 ) then
        x(1)=1.d-4
      endif
      e=dexp(x(1))-1.d0
      f1=2.d0
      yint=norint(-x(2))
c imsld      call mdnord(-x(2),yint)
      x3=2.d0*yint
      yint1=norint(-x(2)+f1)
      yint2=norint(-x(2)-f1)
c imsld      call mdnord(-x(2)+f1,yint1)
c imsld      call mdnord(-x(2)-f1,yint2)
      x4=yint1+yint2
      fx=(0.4d0-(1.d3*e-x3)*x4/(e+x4))/x(1)
     f     -1.d1*dlog(x(1)-1.d-5)
      return
c********************************* tridia quadratisch
 2500 continue
      sum=x(1)**2
      do i=2,n
        sum=sum+i*x(i)**2+x(i)*x(i-1)*2.d0
      enddo
      do i=1,n
        sum=-x(i)*(2*i+2)+sum
      enddo
      sum=sum+x(1)+x(n)
      fx=sum
      return
c********************************* kowalik & osborne
 2600 continue
      sum=0.d0
      x1=x(1)
      x2=x(2)
      x3=x(3)
      x4=x(4)
      do i=1,11
        sum=sum+
     f  (y26(i)-(x1*u26(i)*(x2+u26(i))
     f  /((u26(i)+x3)*u26(i)+x4)))**2
      enddo
      fx=sum
      return
c********************************* bard
 2700 continue
      sum=0.d0
      x1=x(1)
      x2=x(2)
      x3=x(3)
      do i=1,15
        f1=dble(float(i))
        f2=16.d0-f1
        f3=8.d0-dabs(8.d0-f1)
        sum=sum+(y27(i)-(x1+f1/(x2*f2+x3*f3)))**2
      enddo
      fx=sum
      return
c****************************** dennis & brown
 2800 continue
      sum=0.d0
      x1=x(1)
      x2=x(2)
      x3=x(3)
      x4=x(4)
      do i=1,20
        f1=dble(float(i))*.2d0
        sum=sum+
     f    ((x1+x2*f1-dexp(-f1))**2+(x3+x4*dsin(f1)-dcos(f1))**2)**2
      enddo
      fx=sum
      return
c*************************************** nash cal1
 2900 continue
      a=1.d0
      b=2.d0
      call ansatz(x)
      do  i=0,n-1
	 wa(i)=xt(i)**2+dxt(i)*datan(dxt(i))-.5d0*dlog(1.d0+dxt(i)**2)
      enddo
c**** berechnung von f(x) =integral(...) mit der summierten trapezregel
      fx=0.d0
      do  i=0,n-2
	 fx=fx+wa(i+1)+wa(i)
      enddo
      fx=fx*h*.5d0
      return 
c*************************************** nash cal2
 3000 continue
      a=0.d0
      b=0.d0
      call ansatz(x)
      do  i=0,n-1
	 wa(i)=100.d0*(xt(i)-dxt(i)**2)**2+(1.d0-dxt(i))**2
      enddo
c**** berechnung von f(x) =integral(...) mit der summierten trapezregel
      fx=0.d0
      do  i=0,n-2
	 fx=fx+wa(i+1)+wa(i)
      enddo
      fx=fx*h*.5d0
      return
c*************************************** nash cal3
 3100 continue
      a=1.d0
      b=0.d0
      call ansatz(x)
      do  i=0,n-1
	 f1=dexp(-2.d0*xt(i)**2)
         wa(i)=f1*(dxt(i)**2-1.d0)
      enddo
c**** berechnung von f(x) =integral(...) mit der summierten trapezregel
      fx=0.d0
      do  i=0,n-2
	 fx=fx+wa(i+1)+wa(i)
      enddo
      fx=fx*h*.5d0
      return
c*************************************** kettenlinie
 3200 continue
      call fandh(x,ifx,ihx,ifgrad,ihgrad)
      nue=skal1(1,n,ifgrad,ihgrad)
      ghn=skal1(1,n,ihgrad,ihgrad)
      fx=ifx-nue*ihx/ghn+rho*ihx**2/ghn
      return
c**********************************************
 3300 continue
      sum=0.d0
      do i=1,n
        sump(i)=0.d0
      enddo
      do j=1,n
        arg=x(j)
        t0=1.d0
        t1=2.d0*arg-1.d0
        sump(1)=sump(1)+t1
        do k=2,n
          t2=(4.d0*arg-2.d0)*t1-t0
          t0=t1
          t1=t2
          sump(k)=sump(k)+t1
        enddo
      enddo
      do i=1,n
        if ( (i/2)*2 .eq. i ) then
          sump(i)=1.d0/(1.d0-dble(float(i*i)))-sump(i)/dble(float(n))
        else
          sump(i)=sump(i)/dble(float(n))
        endif
        sum=sum+sump(i)**2
      enddo
      fx=sum
      return
 3400 continue
      f1=x(1)**4+x(2)**4
      f2=x(1)**2+x(2)**2
      fx=f1+2.d0*f2-4.d0*x(1)+3.d0
      return
c***********************************************
 3500 continue
      f1=x(1)**2+x(2)**2+x(3)**2-1.d0
      f2=x(1)**2+x(2)**2+(x(3)-2.d0)**2-1.d0
      f3=x(1)+x(2)+x(3)-1.d0
      f4=x(1)+x(2)-x(3)+1.d0
      f5=x(1)**3+3.d0*x(2)**2+(5.d0*x(3)-x(1)+1.d0)**2-36.d0
      fx=f1**2+f2**2+f3**2+f4**2+f5**2
      return
c***********************************************
 3600 continue
      f1=dexp(x(1))-x(2)
      f2=x(2)-x(3)
      f3=dtan(x(3)-x(4))
      f4=x(4)-1.d0
      fx=f1**4+1.d2*f2**6+f3**4+x(1)**8+f4**2
      return
c***********************************************
 3700 continue
      sum=0.d0
      do i=2,n
        f1=x(i-1)-x(i)**2
        f2=1.d0-x(i)
        sum=sum+4.d0*delcha(i)*f1**2+f2**2
      enddo
      fx=sum
      return
c***********************************************
 3800 continue
      prod=1.d0
      do i=1,4
        f1=x(2*i-1)
        f2=x(2*i)
        f3=f1**2+f2**2
        f4=f3**2.5d0*(dsqrt(f3)/6.d0-.4d0)+(.5d0-.2d0*dsqrt(f3))*
     f          f3*f2-47.d0/15.d0
        prod=prod*f4
      enddo
      fx=prod
      return
c***********************************************
 3900 continue
      sum=0.d0
      do i=1,n-1
        sum=sum+x(i)**2*(x(i+1)**2+1.d0)
      enddo
      sum=sum+x(n)**2
      fx=sum
      return
c**********************************************
 4000 continue
      nn=(n+2)/2
      sum=0.d0
      gam1=(1.d0-1.d0/dsqrt(5.d0))/2.d0
      gam2=1.d0-gam1
      do i=0,nn-1
c***    f1=x(2*i-1)  f2=x(2*i)  f3=x(2*i+1) f4=x(2*i+2)
        if ( i .eq. 0 ) then
          f1=0.d0
          f2=0.d0
          f3=x(1)
          f4=x(2)
        elseif ( i .eq. nn-1 ) then
          f1=x(n-1)
          f2=x(n)
          f3=0.d0
          f4=0.d0
        else
          j=2*i-1
          f1=x(j)
          f2=x(j+1)
          f3=x(j+2)
          f4=x(j+3)
        endif
        ksi1=3.d0*nn*(f3-f1)-2.d0*f2-f3
        ksi2=2.d0*nn*(f1-f3)+f2+f4
        f5=f2+2.d0*gam1*ksi1+3.d0*gam1**2*ksi2
        f6=f2+2.d0*gam2*ksi1+3.d0*gam2**2*ksi2
        sum=sum+(2.d0*dsqrt(1.d0+f2**2)+
     f          5.d0*(dsqrt(1.d0+f5**2)+dsqrt(1.d0+f6**2))-12.d0)
        enddo
        sum=-sum/(12.d0*nn)
        s=0.d0
        do i=1,n
          s=s+(x(i)-1.d0/dble(float(i)))**8
        enddo
        fx=sum+s*.01d0
        return
c************************************************************
 4100   continue
	x1=x(1)/100.d0
	x2=x(2)
	cb(1)=1.5d0
	cb(2)=2.25d0
	cb(3)=2.625d0
	fx=0.d0
	do i=1,3
	 fx=(cb(i)-x1*(1.d0-x2**i))**2+fx
	enddo
	return
c**********************************************************
 4200   continue
	fx=1.1*x(1)**2+1.2*x(2)**2-2.d0*x(1)*x(2)-7.d0*x(1)
     *     -3.d0*x(2)+dsqrt(1.d0+x(1)**2+x(2)**2)
	return
 4300   continue
        fx=1.d4*(x(2)-x(1)**2)**2+(1.d0-x(1))**2
        return
      end    
c***************************************************************
      block data
      implicit real*8 (a-h,o-z)
      parameter (nmax=300)
      logical test1,test2,test3,analyt,test4
      common/gradinf/analyt
      common/param/epsx,epsg,epsf,test1,test2,test3,test4
      common/mach/epsmach,tolmach,deldif
      common/dim/n
      character*30 xifile,xofile
      character*80 ident
      common/ident/ident,xifile,xofile
      common/fdat/xst(nmax),icase
      data ident/'xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx'/
      data n/0/,xst/300*0.d0/,test1,test2,test3,test4/4*.false./
      data analyt/.false./
      end
c******************************************************************
      subroutine setda1
      real*8 d,c,cl,y,y0(19),h
      common /dat/ d(19),c(19),cl(19),y(19)
      data y0/1.89d-3,.1038d0,.268d0,.506d0,.577d0,.604d0,.725d0,.898d0,
     1        .947d0,.845d0,.702d0,.528d0,.385d0,.257d0,.159d0,.869d-1,
     2        .453d-1,.1509d-1,.189d-2/
      c(1)=.1d0/7.658d0
      cl(1)=dlog(c(1))
      h=12.d0/dsqrt(6.2832d0)
      d(1)=h/c(1)
      y(1)=y0(1)
      do      100      i=2,19
      c(i)=dble(float(i-1))/7.658d0
      cl(i)=dlog(c(i))
      d(i)=h/c(i)
       y(i)=y0(i)
  100  continue
      return
      end
c************************************************************
      subroutine setda2
      implicit none
      real*8 data,datb
      integer i
      common /dat2/ data(99),datb(99)
      do i=1,99
      data(i)=25.d0+
     f   (-50.d0*dlog(.01d0*dble(float(i))))**.666666666666666d0
      datb(i)=-dble(float(i))*1.d-2
      enddo
      return
      end
c**********************************************************
c     fuer cal1 , cal2 , cal3 
      subroutine ansatz(x)
c***  berechnung der ansatzfunktion kubische b_splines
      implicit real*8 (a-h,o-z)

      parameter (nmax=300)
      integer j,m,n
      real*8  h,a,b,c,c1,c2,x(nmax),xt(0:nmax-1),dxt(0:nmax-1)
      common /bspline/ xt,dxt,h,a,b
      common/dim/n

      m=n-1
      h=1.d0/m
      c=b-a
c**** berechnung von xt(j)=x(tj)=x(1)*phi(tj)+...x(n+1)*phi(tj)+a*(1-tj)
c**** + b*tj und dxt(j)=x'(tj)=x(1)*phi'(tj)+...x(n+1)*phi'(tj)-a+b
      xt(0)=a
      xt(m)=b
      dxt(0)=m*(x(1)*2.d0+x(2))+c
      dxt(m)=-m*(x(m+1)*2.d0+x(m))+c
      c1=c*h
      c2=0.d0
      do  j=1,m-1
        c2=c2+c1
	xt(j)=(x(j)+x(j+2))/6.d0+x(j+1)*2.d0/3.d0+a+c2
	dxt(j)=m*(x(j+2)-x(j))*.5d0+c
      enddo
      return
      end
c****************************************************************************
c fuer var2fct = beispiel 32
      subroutine fandh(x,fx,hx,gradfx,gradhx)

      implicit real*8 (a-h,o-z)
      parameter (nmax=300,nmax1=301)

      real*8 kl,x(nmax),gradfx(nmax),gradhx(nmax)

      logical test1,test2,test3,test4
      common/param/epsx,epsg,epsf,test1,test2,test3,test4
      common/mach/epsmach,tolmach,deldif
      common/dim/n
      common/kette/rho,kl,ya,ye
      common/funcdat/fxgl,hxgl
      real*8 y(0:nmax1),ys(nmax1),ym(nmax1),dt,dth,w(nmax1),
     f       wr(nmax1)
      n1=n+1
      dt=1.d0/dble(float(n1))
      dth=.5d0*dt
      y(0)=ya
      y(n1)=ye
      do i=1,n
        y(i)=x(i)
        gradhx(i)=0.d0
        gradfx(i)=0.d0
      enddo
      do i=1,n1
        ym(i)=(y(i)+y(i-1))*.5d0
        ys(i)=(y(i)-y(i-1))/dt
        w(i)=dsqrt(1.d0+ys(i)**2)
        wr(i)=1.d0/w(i)
      enddo
      hx=w(1)+w(n1)
      fx=ym(1)*w(1)+ym(n1)*w(n1)
      gradhx(1)=wr(1)*ys(1)
      gradfx(1)=wr(1)*ys(1)*ym(1)+w(1)*dth
      gradhx(n)=-wr(n1)*ys(n1)
      gradfx(n)=-wr(n1)*ys(n1)*ym(n1)+w(n1)*dth
      do i=2,n
        fx=fx+ym(i)*w(i)
        hx=hx+w(i)
        term1=wr(i)*ys(i)
        term2=term1*ym(i)
        term3=w(i)*dth
        gradhx(i)=gradhx(i)+term1
        gradhx(i-1)=gradhx(i-1)-term1
        gradfx(i)=gradfx(i)+term2+term3
        gradfx(i-1)=gradfx(i-1)-term2+term3
      enddo
      fx=fx*dt
      hx=hx*dt-kl
      fxgl=fx
      hxgl=hx
c************** multiplikation der gradienten mit dt entfaellt,
c************** da innere ableitung 1/dt nicht beruecksichtigt
      return
      end      
c****************************************************************
c     computes gradient by high precision numerical differentiation
c*******************************************************************
c copyright = spellucci 
      subroutine gradf(x,gradx,err)
      parameter (nx=300)
      implicit real*8 (a-h,o-z)
      logical err
      common/count/icf,icgf,iup,inoup,ires,imod
      real*8 x(nx),gradx(nx),fx(6)
      logical test1,test2,test3,test4
      common/param/epsx,epsg,epsf,test1,test2,test3,test4
      common/mach/epsmach,tolmach,deldif
      common/dim/n
c****************************************************
c     high precision numerical differentiation
c     by sixth order extrapolation
c****************************************************
      err=.false.
      do   i=1,n
        z=x(i)
        xincr=deldif*(dabs(x(i))+1.d0)
        x(i)=z-xincr
        call f(x,fx(1),err)
        if ( err ) return
        x(i)=z+xincr
        call f(x,fx(2),err)
        if ( err ) return
      	xincr=xincr+xincr
        d1=xincr
        x(i)=z-xincr
        call f(x,fx(3),err)
        if ( err ) return
        x(i)=z+xincr
        call f(x,fx(4),err)
        if ( err ) return
        xincr=xincr+xincr
        d2=xincr
        x(i)=z-xincr
        call f(x,fx(5),err)
        if ( err ) return
        x(i)=z+xincr
        call f(x,fx(6),err)
        if ( err ) return
        x(i)=z
        d3=xincr+xincr
        sd1=(fx(2)-fx(1))/d1
        sd2=(fx(4)-fx(3))/d2
        sd3=(fx(6)-fx(5))/d3
        sd3=sd2-sd3
        sd2=sd1-sd2
        sd3=sd2-sd3
        gradx(i)=sd1+.4d0*sd2+sd3/45.d0
      enddo
      icgf=icgf+1
      return
      end
c*******************************************************************
c*******************************************************************
      double precision function norint(x)
      implicit none
      real*8 x
c********* computes the gaussian normal distribution integral
c********* precision about 16 digits
      real*8 p1(0:8),q1(0:9),p2(0:5),q2(0:6)
      real*8 sqrt2,rsqrtpi,arg,arg2,xabs,erf,erfc
c********** coefficients taken from hart computer approximations
c********** no erfc 5708 and erfc 5725
      data p1/.37235079815548067d4, .71136632469540499d4,
     f        .67582169641104859d4, .40322670108300497d4,
     f        .16317602687537147d4, .45626145870609263d3,
     f        .86082762211948595d2, .10064858974909542d2,
     f        .56418958676181361/
      data q1/.37235079815548065d4, .11315192081854405d5,
     f        .15802535999402043d5, .13349346561284457d5,
     f        .75424795101934758d4, .29680049014823087d4,
     f        .81762238630454408d3, .15307771075036222d3,
     f        .17839498439139557d2, .1d1/
      data p2/.29788656263939929d1, .74097406059647418d1,
     f        .61602098531096305d1, .50190497267842675d1,
     f        .12753666447299660d1, .56418958354775507/
      data q2/.33690752069827528d1, .96089653271927879d1,
     f        .17081440747466004d2, .12048951927855129d2,
     f        .93960340162350542d1, .22605285207673270d1, .1d1/
      data sqrt2/1.41421356237390505d0/
      data rsqrtpi/.56418958354775629d0/
      xabs=dabs(x)
      if ( xabs .gt. .5d0 ) then
	if ( xabs .gt. 8.d0 ) then
c******* take approximation 5725
	  if ( xabs .gt. 100.d0 ) then
	    erfc=0.d0
	  else
	    arg=xabs/sqrt2
	    erfc=(((((p2(5)*arg+p2(4))*arg+p2(3))*arg+p2(2))*arg+p2(1)
     f             )*arg+p2(0))/
     f           ((((((arg+q2(5))*arg+q2(4))*arg+q2(3))*arg+q2(2))*arg
     f                 +q2(1))*arg+q2(0))*dexp(-arg**2)
	  endif
	else
c******** take approximation 5708
	  arg=xabs/sqrt2
	  erfc=((((((((p1(8)*arg+p1(7))*arg+p1(6))*arg+p1(5))*arg
     f                 +p1(4))*arg+p1(3))*arg+p1(2))*arg+p1(1))*arg
     f                 +p1(0))/
     f         (((((((((arg+q1(8))*arg+q1(7))*arg+q1(6))*arg+q1(5))
     f                 *arg+q1(4))*arg+q1(3))*arg+q1(2))*arg+q1(1))
     f                 *arg+q1(0))*dexp(-arg**2)
	endif
	if ( x .lt. 0.d0 ) then
	  norint=erfc*.5d0
	  return
	else
	  norint=(2.d0-erfc)*.5d0
	  return
	endif
      else
c************ take power series approximation
	arg=xabs/sqrt2
	arg2=arg**2
	erf=arg*2.d0*rsqrtpi*((((((((((arg2/210.d0-1.d0/19.d0)*arg2
     f     /9.d0+1.d0/17.d0)*arg2/8.d0-1.d0/15.d0)*arg2/7.d0+1.d0/13.d0)
     f     *arg2/6.d0-1.d0/11.d0)*arg2/5.d0+1.d0/9.d0)*arg2/4.d0
     f     -1.d0/7.d0)*arg2/3.d0+1.d0/5.d0)*arg2/2.d0-1.d0/3.d0)*arg2
     f     +1.d0)
	if ( x .ge. 0.d0 ) then
	  norint=(1.d0+erf)*.5d0
	  return
	else
	  norint=(1.d0-erf)*.5d0
	  return
	endif
      endif
      end
