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