********************************************************** * Date created: August 8, 2007 * * * * UNCONSTRAINED OPTIMIZATION BY DIRECT SEARCHING * ================================================ * * * Neculai Andrei, * Research Institute for Informatics * Advanced Modeling and Optimization Laboratory * 8-10, Bdl. Maresal Alexandru Averescu * 71316, Bucharest 1, Romania * and * Academy of Romanian Scientists, * 54 Splaiul Independentei, Bucharest 5, Romania *-------------------------------------------------------------------- * * * This Fortran 77 package implements a number of 4 direct search * methods for unconstrained optimization. * * The methods considered are as follows: * 1) Nelder-Mead * 2) Rosenbrock * 3) Hook-Jeeves * 4) Powell * * Hook-Jeeves and Powell methods work in close relation with 5 * one dimensional search: * 1) Golden search * 2) Fibonacci search * 3) Quadratic interpolation by Powell * 4) Dichotomous search * 5) lambda=1 *-------------------------------------------------------------------- * * * * * Main program * DOUBLE PRECISION X(40) DOUBLE PRECISION EPS, vfi integer*4 timpexp real*8 time1,time2, timp integer method, itech character*40 numef, fnumef(200) COMMON /PREC/EPS c Input file open(unit=7,file='funcname.txt',status='old') c Output files open(unit=8,file='d1.out',status='unknown') open(unit=9,file='d1.rez',status='unknown') *------------------------------------------------------------------ eps = 0.000000001d0 nexpst = 78 nexptot = 78 c c method = 1 (Nelder - Mead) c 2 (Rosenbrock) c 3 (Hook - Jeeves) (itech) c 4 (Powell) method = 2 c c itech = 1 (golden search) c 2 (Fibonacci) c 3 (Powell quadratic interpolation) c 4 (dichotomous search) c 5 (1) itech = 1 do i=1,78 read(7,21) numef fnumef(i)=numef 21 format(a40) end do *---------------------------------------------- Start experiments do nexp = nexpst, nexptot itert = 0 nfunct= 0 timp = 0 numef = fnumef(nexp) write(8,41) nexp, numef write(*,41) nexp, numef 41 format(/,2x,i3,2x,'UNO Algorithm. Function: ',a40) if(method .eq. 1) write(8,42) 42 format(9x,'Nelder - Mead method',/) if(method .eq. 2) write(8,43) 43 format(9x,'Rosenbrock method',/) if(method .eq. 3) then if(itech .eq. 1) write(8,44) 44 format(9x,'Hook - Jeeves method with Golden search',/) if(itech .eq. 2) write(8,441) 441 format(9x,'Hook - Jeeves method with Fibonacci search',/) if(itech .eq. 3) write(8,442) 442 format(9x,'Hook - Jeeves method with Powell QI search',/) if(itech .eq. 4) write(8,443) 443 format(9x,'Hook - Jeeves method with dichotomous',/) if(itech .eq. 5) write(8,444) 444 format(9x,'Hook - Jeeves method with alpha=1',/) end if if(method .eq. 4) then if(itech .eq. 1) write(8,54) 54 format(9x,'Powell method with Golden search',/) if(itech .eq. 2) write(8,541) 541 format(9x,'Powell method with Fibonacci search',/) if(itech .eq. 3) write(8,542) 542 format(9x,'Powell method with Powell QI search',/) if(itech .eq. 4) write(8,543) 543 format(9x,'Powell method with dichotomous',/) if(itech .eq. 5) write(8,544) 544 format(9x,'Powell method with alpha=1',/) end if write(8,19) 19 format(9x,'n',5x,'iter',4x,'nfunc',4x,2x,'time(c)', * 15x,'f') write(8,20) 20 format(1x,94('-')) do n = 6,6,6 c Call INIPOINT (initial guess) call inipoint (n,x, nexp) c Call UNO for optimization call cpu_time(time1) call unos(n,x, method,itech, iter,nfunc,vfi, nexp) call cpu_time(time2) *-- timpexp = int(100*(time2-time1)) itert = itert + iter nfunct = nfunct + nfunc timp = timp + float(timpexp) write(8,51) n, iter, nfunc, timpexp, vfi write(*,51) n, iter, nfunc, timpexp, vfi 51 format(4x,i6,3x,i6,3x,i6,3x,i6,5x,e20.13) if(n .eq. 5) then write(9,611)nexp, n,iter, nfunc, timpexp,vfi 611 format(i2,i6,2x,i6,2x,i6,2x,i6,2x,e20.13) else write(9,61) n,iter, nfunc, timpexp,vfi 61 format(2x,i6,2x,i6,2x,i6,2x,i6,2x,e20.13) end if end do c--------------------------------------- End do n write(8,20) write(8,71) itert, nfunct, timp/100.d0 71 format(1x,'TOTAL',7x,i6,3x,i6,2x,f7.2,'(seconds)') end do c---------------------------------------------End do nexp stop end c c Last card MAINPROGRAM Neculai Andrei *============================================================================= subroutine unos(n,x,method,itech,iter,nfunc,vfi, nexp) real*8 X(40),S(40),H(41),P(1600) real*8 X1(40),X2(40),X3(40),Y(1600),Z(1600) real*8 EPS, vfi C COMMON /PREC/EPS c---------------------------------------------------------method=1 (Nelder-Mead method) if(method .eq. 1) then N1=N+1 NN=N*(N+1) C CALL NELMED(N,N1,NN,X,X1,X2,X3,Y,H, iter,nfunc,vfi,nexp) C write(*,884)(i,x(i),i=1,n) 884 format(4x,i4,4x,e20.13) write(*,886) vfi 886 format(4x,'Final value:',e20.13) write(8,884)(i,x(i),i=1,n) write(8,886) vfi end if c---------------------------------------------------------method=2 (Rosenbrock) if(method .eq. 2) then NN=N*N C CALL ROSE(N,NN,X,X1,X2,S,H,Y,Z,P,X3, iter,nfunc,vfi,nexp) C write(*,884)(i,x(i),i=1,n) write(*,886) vfi write(8,884)(i,x(i),i=1,n) write(8,886) vfi end if c---------------------------------------------------------method=3 (Hook - Jeeves) if(method .eq. 3) then C CALL HOOKJ(N,X,X2,Y,Z,S,ITECH, iter,nfunc,vfi, nexp) C write(*,884)(i,x(i),i=1,n) write(*,886) vfi write(8,884)(i,x(i),i=1,n) write(8,886) vfi end if c--------------------------------------------------------method=4 (Powell) if(method .eq. 4) then C NN=N*N CALL POWEL(N,NN,X,X1,X2,X3,S,Y,ITECH, * iter,nfunc,vfi, nexp) C write(*,884)(i,x(i),i=1,n) write(*,886) vfi write(8,884)(i,x(i),i=1,n) write(8,886) vfi end if return end C********************************************************* C C Nelder - Mead Optimization method. Nelder-Mead C Simplex Method. C C********************************************************* C Neculai Andrei, C Research Institute for Informatics C Advanced Modeling and Optimization Laboratory C********************************************************* C SUBROUTINE NELMED(N,N1,NN,X,X0,XR,XE,XD,VF, * iter,nfunc,vfi,nexp) C DOUBLE PRECISION X(N),X0(N),XR(N),XE(N) DOUBLE PRECISION XD(NN),VF(N1),VF1(100) C DOUBLE PRECISION ALFA,BETA,GAMA,EPS DOUBLE PRECISION FX,F0,FR,FE,VFI,VFM,S,T C COMMON /PREC/EPS C NFUNC=0 ITER=0 ALFA=1.0 BETA=0.5 GAMA=2.0 C CALL FUNC(N,X,FX,nexp) NFUNC=NFUNC+1 write(8,881) fx 881 format(4x,'Initial value:',e20.13) C C C------------- Build up the initial simplex C using the initial point X. C N2=N*N DO 1 J=1,N L=J+N2 1 XD(L)=X(J) C DO 3 I=1,N DO 2 J=1,N L=J+(I-1)*N 2 XD(L)=X(J) L=I+(I-1)*N 3 XD(L)=XD(L)+0.5 C C C-----------------Compute the values of the objective C function in vertices of the simplex C DO 5 I=1,N1 DO 4 J=1,N L=J+(I-1)*N 4 X(J)=XD(L) CALL FUNC(N,X,FX,nexp) NFUNC=NFUNC+1 VF(I)=FX 5 CONTINUE C C 6 ITER=ITER+1 C IF(ITER.GT.1000)THEN RETURN ENDIF C * do i=1,n1 do j=1,n l=j+(i-1)*n x(j)=xd(l) end do call func(n,x,fx,nexp) vf1(i)=fx end do C C VFI=VF1(1) VFM=VFI IMI=1 IMA=1 DO 7 I=2,N1 IF(VF1(I).LT.VFI) THEN VFI=VF1(I) IMI=I ENDIF IF(VF1(I).GT.VFM) THEN VFM=VF1(I) IMA=I ENDIF 7 CONTINUE C C C-----------------Compute the centroid of the simplex. C DO 9 J=1,N X0(J)=0.0 DO 8 I=1,N1 IF(I.EQ.IMA) GO TO 8 L=J+(I-1)*N X0(J)=X0(J)+XD(L) 8 CONTINUE 9 CONTINUE C DO 10 J=1,N 10 X0(J)=X0(J)/FLOAT(N) CALL FUNC(N,X0,F0,nexp) NFUNC=NFUNC+1 C C C----------------Compute the reflection point. C DO 11 J=1,N L=J+(IMA-1)*N 11 XR(J)=(1.+ALFA)*X0(J)-ALFA*XD(L) CALL FUNC(N,XR,FR,nexp) NFUNC=NFUNC+1 C C C----------------Decision on the way. C IF(FR.LT.VFI) THEN C DO 12 J=1,N 12 XE(J)=GAMA*XR(J)+(1.-GAMA)*X0(J) CALL FUNC(N,XE,FE,nexp) NFUNC=NFUNC+1 C C IF(FE.LT.VFI) THEN DO 13 J=1,N L=J+(IMA-1)*N 13 XD(L)=XE(J) VF1(IMA)=FE GO TO 23 ELSE DO 16 J=1,N L=J+(IMA-1)*N 16 XD(L)=XR(J) VF1(IMA)=FR GO TO 23 ENDIF C ELSE C IS=1 DO 14 I=1,N1 IF(I.EQ.IMA) GO TO 14 IF(FR.GT.VF1(I)) GO TO 14 IS=0 GO TO 15 14 CONTINUE 15 IF(IS.EQ.1) GO TO 18 DO 17 J=1,N L=J+(IMA-1)*N 17 XD(L)=XR(J) VF1(IMA)=FR GO TO 23 C ENDIF C C----------------Compute the contraction point. C 18 CONTINUE IF(FR.GT.VFM) THEN DO 19 J=1,N L=J+(IMA-1)*N 19 X(J)=BETA*XD(L)+(1.-BETA)*X0(J) ELSE DO 20 J=1,N 20 X(J)=BETA*XR(J)+(1.-BETA)*X0(J) ENDIF C CALL FUNC(N,X,FX,nexp) NFUNC=NFUNC+1 C C C----------------Reduce the simplex. C IF(FX.GT.VFM) THEN DO 21 I=1,N1 DO 21 J=1,N L=J+(I-1)*N K=J+(IMI-1)*N 21 XD(L)=(XD(L)+XD(K))/2. C GO TO 23 C ELSE C DO 22 J=1,N L=J+(IMA-1)*N 22 XD(L)=X(J) GO TO 23 C ENDIF C C--------------Test of continuation (accuracy). C 23 CONTINUE S=0.0 DO 24 I=1,N1 24 S=S+(VF1(I)-F0)*(VF1(I)-F0) S=S/FLOAT(N+1) T=SQRT(S) C C IF(T.LE.EPS) GO TO 30 GO TO 6 C 30 CONTINUE DO 31 J=1,N L=J+(IMI-1)*N 31 X(J)=XD(L) C RETURN END C C*****************************Last card of NELMED.FOR C********************************************************** C C Rosenbrock's method for unconstrained minimization. Rosenbrook C C********************************************************** C Neculai Andrei, C Research Institute for Informatics C Advanced Modeling and Optimization Laboratory C********************************************************** C SUBROUTINE ROSE(N,NN,X,X1,X2,Z,H,HM,XD,P,Y, * iter,nfunc,vfi,nexp) C DOUBLE PRECISION X(N),X1(N),X2(N),Z(N),H(N),Y(N) DOUBLE PRECISION HM(NN),XD(NN),P(NN) C DOUBLE PRECISION EPS,FX,FY,S,T,ALFA,BETA, vfi C COMMON /PREC/EPS C NFUNC=0 IR=0 C ALFA=3.0D0 BETA=5.0D-1 C DO 1 I=1,N 1 H(I)=8.0D-1 C ITER=1 C CALL FUNC(N,X,FX,nexp) NFUNC=NFUNC+1 write(8,881) fx 881 format(4x,'Initial value:',e20.13) C DO 3 I=1,N DO 2 J=1,N L=J+(I-1)*N 2 XD(L)=0.0 L=I+(I-1)*N 3 XD(L)=1.0 C DO 4 I=1,N X1(I)=0.0 X2(I)=0.0 4 Z(I)=0.0 C 5 CONTINUE C I=1 6 DO 7 J=1,N L=J+(I-1)*N 7 Y(J)=X(J)+H(I)*XD(L) C CALL FUNC(N,Y,FY,nexp) NFUNC=NFUNC+1 C IF(DABS(FX-FY).LE.EPS) IR=IR+1 IF(IR.GT.20) GO TO 30 C IF(FY.LT.FX) THEN C C--------------------------------------------------------Succes. C Z(I)=Z(I)+H(I) H(I)=ALFA*H(I) X1(I)=1.0 FX=FY DO 8 J=1,N 8 X(J)=Y(J) ITER=ITER+1 IF(ITER.GT.1000) THEN RETURN ENDIF GO TO 9 C ELSE C C-------------------------------------------------------Failure. C H(I)=-BETA*H(I) X2(I)=1.0 GO TO 9 C ENDIF C 9 I=I+1 IF(I.LE.N) GO TO 6 C IS=1 DO 10 J=1,N IF(X1(J)*X2(J).NE.0.0) GO TO 10 IS=0 GO TO 11 10 CONTINUE 11 IF(IS.EQ.0) GO TO 5 C IS=1 DO 12 J=1,N IF(DABS(Z(J)).LE.EPS) GO TO 12 IS=0 GO TO 13 12 CONTINUE 13 IF(IS.EQ.1) GO TO 30 C DO 15 I=1,N DO 14 J=1,N L=J+(I-1)*N 14 HM(L)=Z(J) Z(I)=0.0 15 CONTINUE C DO 16 I=1,N DO 16 J=1,N L=I+(J-1)*N P(L)=0.0 DO 16 K=1,N L1=I+(K-1)*N L2=K+(J-1)*N 16 P(L)=P(L)+XD(L1)*HM(L2) C S=0.0 DO 17 J=1,N 17 S=S+P(J)*P(J) T=SQRT(S) DO 18 J=1,N 18 XD(J)=P(J)/T C DO 25 I=2,N DO 19 K=1,N L=K+(I-1)*N 19 Y(K)=P(L) C DO 22 J=1,I-1 S=0.0 DO 20 K=1,N L1=K+(I-1)*N L2=K+(J-1)*N 20 S=S+P(L1)*XD(L2) C DO 21 K=1,N L=K+(J-1)*N 21 Y(K)=Y(K)-S*XD(L) C 22 CONTINUE C S=0.0 DO 23 K=1,N 23 S=S+Y(K)*Y(K) C T=SQRT(S) DO 24 K=1,N L=K+(I-1)*N 24 XD(L)=Y(K)/T 25 CONTINUE C DO 26 I=1,N X1(I)=0.0D0 X2(I)=0.0D0 Z(I)=0.0D0 26 H(I)=8.0D-1 C GO TO 5 C 30 CONTINUE C vfi=FX C RETURN END C C***************************Last card of ROSE.FOR C******************************************************* C C Searching method by Hook-Jeeves. Hook - Jeeves C C******************************************************* C Neculai Andrei, C Research Institute for Informatics C Advanced Modeling and Optimization Laboratory C******************************************************* C SUBROUTINE HOOKJ(N,X1,X2,Y,Z,S,ITECH, * iter,nfunc,vfi, nexp) C DOUBLE PRECISION X1(N),X2(N),Y(N),Z(N),S(N) DOUBLE PRECISION F,FP,FM,BL,FY,F2,H,EPS,A,vfi C COMMON /PREC/EPS COMMON /L3A/ICOD COMMON /NOF/NBL C C ICOD=0 H=8.0D-1 NFUNC=0 NBL=0 C K=1 CALL FUNC(N,X1,F,nexp) NFUNC=NFUNC+1 c write(8,881) f c881 format(4x,'Initial value:',e20.13) C 20 CONTINUE DO 1 J=1,N 1 Y(J)=X1(J) C 10 DO 5 I=1,N DO 2 J=1,N 2 Z(J)=Y(J) C CALL FUNC(N,Z,F,nexp) NFUNC=NFUNC+1 Z(I)=Z(I)+H C CALL FUNC(N,Z,FP,nexp) NFUNC=NFUNC+1 C IF(FP.LT.F) GO TO 3 Z(I)=Z(I)-2.*H C CALL FUNC(N,Z,FM,nexp) NFUNC=NFUNC+1 C IF(FM.LT.F.AND.F.LT.FP) GO TO 3 A=FP IF(FM.LT.A) A=FM IF(F.LT.A) GO TO 5 3 DO 4 J=1,N 4 Y(J)=Z(J) 5 CONTINUE C C IS=0 DO 6 I=1,N IF(X1(I).NE.Y(I)) IS=1 6 CONTINUE IF(IS.EQ.1) THEN DO 7 I=1,N 7 X2(I)=Y(I) GO TO 8 ELSE IF(H.LT.EPS) GO TO 100 H=H/2. GO TO 10 ENDIF C C 8 DO 9 I=1,N 9 S(I)=X2(I)-X1(I) C BL=0.d0 C---------------------------------------------------------------- GO TO (171,172,173,174,175) ITECH 171 CALL L1(N,X2,S,BL,nexp) GO TO 180 172 CALL L2(N,X2,S,BL,nexp) IF(EPS.EQ.0.0) RETURN GO TO 180 173 CALL L3(N,X2,S,BL,nexp) IF(ICOD.EQ.1) RETURN GO TO 180 174 CALL L4(N,X2,S,BL,nexp) GO TO 180 175 BL=1.0D0 IF(EPS.EQ.0.0) EPS=0.0000001 C 180 CONTINUE C--------------------------------------------------------------- C C NFUNC=NFUNC+NBL C DO 11 I=1,N 11 Y(I)=X2(I)+BL*S(I) K=K+1 IF(K.GT.1000) THEN RETURN ENDIF C CALL FUNC(N,Y,F,nexp) NFUNC=NFUNC+1 C DO 15 I=1,N DO 12 J=1,N 12 Z(J)=Y(J) C CALL FUNC(N,Z,F,nexp) NFUNC=NFUNC+1 C Z(I)=Z(I)+H C CALL FUNC(N,Z,FP,nexp) NFUNC=NFUNC+1 C IF(FP.LT.F) GO TO 13 Z(I)=Z(I)-2.*H C CALL FUNC(N,Z,FM,nexp) NFUNC=NFUNC+1 C IF(FM.LT.F.AND.F.LT.FP) GO TO 13 A=FP IF(FM.LT.A) A=FM IF(F.LT.A) GO TO 15 13 DO 14 J=1,N 14 Y(J)=Z(J) 15 CONTINUE C C CALL FUNC(N,X2,F2,nexp) CALL FUNC(N,Y,FY,nexp) NFUNC=NFUNC+2 C C IF(FY.LT.F2) THEN DO 16 I=1,N 16 X1(I)=X2(I) DO 17 I=1,N 17 X2(I)=Y(I) C IF(DABS(FY-F2).LE.EPS) THEN IF(H.LE.EPS) GO TO 100 H=H/2. ENDIF C GO TO 8 ELSE IF(H.LE.EPS) GO TO 100 H=H/2. DO 18 I=1,N 18 X1(I)=X2(I) GO TO 20 ENDIF C C 100 CONTINUE C iter=k vfi=f C c write(8,884)(i,x1(i),i=1,n) c884 format(4x,i4,4x,e20.13) c write(8,886) vfi c886 format(4x,'Final value:',e20.13) RETURN END C C******************************* Last card of HOOKJ.FOR C********************************************************** C C Powell's Method. Powell C C********************************************************** C Neculai Andrei, C Research Institute for Informatics C Advanced Modeling and Optimization Laboratory C********************************************************** C SUBROUTINE POWEL(N,NN,X,X1,X2,X3,S,XD,ITECH, * iter,nfunc,vfi, nexp) C DOUBLE PRECISION X(N),X1(N),X2(N),X3(N),S(N),XD(NN) DOUBLE PRECISION BL,EP,EPS,F1,F2,F3,FS,T,V,DELTA,DEL DOUBLE PRECISION FA,FB,vfi C COMMON /PREC/EPS COMMON /L3A/ICOD COMMON /NOF/NBL C C NFUNC=0 NBL=0 EP=1.0D-5 C ITER=0 CALL FUNC(N,X,F1,nexp) NFUNC=NFUNC+1 c write(*,881) f1 c881 format(4x,'Initial value:',e20.13) C DO 3 I=1,N DO 2 J=1,N L=J+(I-1)*N 2 XD(L)=0.0 L=I+(I-1)*N 3 XD(L)=1.0 C 4 CONTINUE ITER=ITER+1 C IF(ITER.GT.1000) THEN RETURN ENDIF C I=1 DO 5 J=1,N 5 S(J)=XD(J) C--------------------------- DO 25 J=1,N 25 X1(J)=X(J)+EP*S(J) C CALL FUNC(N,X1,FA, nexp) NFUNC=NFUNC+1 C--------------------------- DO 26 J=1,N 26 X1(J)=X(J)-EP*S(J) C CALL FUNC(N,X1,FB,nexp) NFUNC=NFUNC+1 C IF(FA.GT.FB) THEN DO 27 J=1,N 27 S(J)=-S(J) ENDIF C BL=0.d0 C--------------------------------------------------- GO TO(171,172,173,174) ITECH 171 CALL L1(N,X,S,BL,nexp) GO TO 180 172 CALL L2(N,X,S,BL,nexp) GO TO 180 173 CALL L3(N,X,S,BL,nexp) IF(ICOD.EQ.1) RETURN GO TO 180 174 CALL L4(N,X,S,BL,nexp) C 180 CONTINUE C----------------------------------------------------- C NFUNC=NFUNC+NBL C DO 6 J=1,N 6 X2(J)=X(J)+BL*S(J) C CALL FUNC(N,X2,F2,nexp) NFUNC=NFUNC+1 C DELTA=(F1-F2) M=1 C 7 I=I+1 IF(I.GT.N) GO TO 11 C DO 8 J=1,N L=J+(I-1)*N 8 S(J)=XD(L) C DO 30 J=1,N 30 X1(J)=X(J)+EP*S(J) C CALL FUNC(N,X1,FA,nexp) NFUNC=NFUNC+1 C DO 31 J=1,N 31 X1(J)=X(J)-EP*S(J) C CALL FUNC(N,X1,FB,nexp) NFUNC=NFUNC+1 C IF(FA.GT.FB) THEN DO 32 J=1,N 32 S(J)=-S(J) ENDIF C BL=0.d0 C----------------------------------------------------- GO TO(181,182,183,184) ITECH 181 CALL L1(N,X2,S,BL,nexp) GO TO 190 182 CALL L2(N,X2,S,BL,nexp) GO TO 190 183 CALL L3(N,X2,S,BL,nexp) IF(ICOD.EQ.1) RETURN GO TO 190 184 CALL L4(N,X2,S,BL,nexp) 190 CONTINUE C---------------------------------------------------- C NFUNC=NFUNC+NBL C DO 9 J=1,N 9 X3(J)=X2(J)+BL*S(J) C CALL FUNC(N,X3,F3,nexp) NFUNC=NFUNC+1 C DEL=(F2-F3) IF(DEL.GT.DELTA) THEN DELTA=DEL M=I ENDIF C DO 10 J=1,N 10 X2(J)=X3(J) F2=F3 GO TO 7 C 11 CONTINUE C IS=0 DO 12 I=1,N IF(DABS(X3(I)-X(I)).LE.EPS) GO TO 12 IS=1 GO TO 121 12 CONTINUE C 121 CONTINUE C IF(IS.EQ.0) GO TO 100 C DO 13 I=1,N 13 S(I)=2.*X3(I)-X(I) C CALL FUNC(N,S,FS,nexp) NFUNC=NFUNC+1 IF(FS.GE.F1) THEN DO 14 I=1,N 14 X(I)=X3(I) F1=F3 GO TO 4 ENDIF C T=(F1-2.*F3+FS)*(F1-F3-DELTA)**2 V=0.5*DELTA*(F1-FS)**2 C IF(T.GE.V) THEN c Here T > V: Consider the old directions C----------------------------------------------------------- DO 15 I=1,N 15 X(I)=X3(I) F1=F3 GO TO 4 ENDIF C DO 16 I=1,N 16 S(I)=X3(I)-X(I) C c Here T < V: Compute the direction X3-X C--------------------------------------------------------- BL=0.d0 GO TO(191,192,193,194) ITECH 191 CALL L1(N,X3,S,BL,nexp) GO TO 200 192 CALL L2(N,X3,S,BL,nexp) GO TO 200 193 CALL L3(N,X3,S,BL,nexp) IF(ICOD.EQ.1) RETURN GO TO 200 194 CALL L4(N,X3,S,BL,nexp) C 200 CONTINUE C--------------------------------------------- C NFUNC=NFUNC+NBL DO 17 I=1,N 17 X(I)=X3(I)+BL*S(I) C CALL FUNC(N,X,F1,nexp) NFUNC=NFUNC+1 IS=0 DO 18 I=1,N IF(DABS(BL*S(I)).LE.EPS) GO TO 18 IS=1 GO TO 118 18 CONTINUE C 118 CONTINUE C IF(IS.EQ.0) GO TO 100 C IF(M.EQ.N) THEN DO 19 J=1,N L=J+(N-1)*N 19 XD(L)=S(J) GO TO 4 ENDIF C DO 20 I=M+1,N DO 20 J=1,N J1=J+(I-2)*N J2=J+(I-1)*N 20 XD(J1)=XD(J2) C DO 21 J=1,N L=J+(N-1)*N 21 XD(L)=S(J) GO TO 4 C 100 CONTINUE C vfi=f1 C RETURN END C C*********************************Last card of POWEL.FOR C****************************************************** C C Subroutine for one dimensional searching. C Golden section with interval searching. Golden section C ============== C****************************************************** C Neculai Andrei, C Research Institute for Informatics C Advanced Modeling and Optimization Laboratory C****************************************************** C SUBROUTINE L1(N,X,S,BL,nexp) DOUBLE PRECISION X(N),S(N),BL DOUBLE PRECISION F0,FY,FY1,SL,SL1,G1,G2,G2P,D1,D2 DOUBLE PRECISION R,R1,EPS,EPS1,Y(40) C COMMON /PREC/EPS COMMON /NOF/NBL C EPS1=0.001 C R=0.618034 R1=1.0d0-R SL=1. NBL=0 C C CALL FUNC(N,X,F0,nexp) NBL=NBL+1 C DO 1 I=1,N 1 Y(I)=X(I)+SL*S(I) CALL FUNC(N,Y,FY,nexp) NBL=NBL+1 C SL1=SL-EPS1 DO 2 I=1,N 2 Y(I)=X(I)+SL1*S(I) CALL FUNC(N,Y,FY1,nexp) NBL=NBL+1 C IF(FY.LT.FY1) GO TO 6 C D1=0.0 D2=SL GO TO 7 C 6 K=1 G1=SL 5 G2=G1+SL C 15 CONTINUE C DO 3 I=1,N 3 Y(I)=X(I)+G2*S(I) CALL FUNC(N,Y,FY,nexp) NBL=NBL+1 C G2P=G2-EPS1 DO 4 I=1,N 4 Y(I)=X(I)+G2P*S(I) CALL FUNC(N,Y,FY1,nexp) NBL=NBL+1 C IF(FY.LE.FY1) THEN IF(K.GT.10) THEN K=1 EPS1=-EPS1 G1=0. G2=SL GO TO 15 ENDIF K=K+1 G1=G2 GO TO 5 ELSE D1=G1 D2=G2 GO TO 7 ENDIF C 7 SL=D2-D1 C IF(SL.LE.EPS) GO TO 10 C G1=D1+R1*SL G2=D1+R*SL C DO 8 I=1,N 8 Y(I)=X(I)+G1*S(I) CALL FUNC(N,Y,FY,nexp) NBL=NBL+1 C DO 9 I=1,N 9 Y(I)=X(I)+G2*S(I) CALL FUNC(N,Y,FY1,nexp) NBL=NBL+1 C IF(FY.LT.FY1) THEN D2=G2 ELSE D1=G1 ENDIF GO TO 7 C 10 BL=(D1+D2)/2 C RETURN END C C*********************************** Last card of L1.FOR C********************************************************** C C Subroutine for one dimensional searching. C Fibonacci searching technique. C ========= C********************************************************** C Neculai Andrei, C Research Institute for Informatics C Advanced Modeling and Optimization Laboratory C********************************************************** C SUBROUTINE L2(N,X,S,BL,nexp) DOUBLE PRECISION X(N),S(N),BL DOUBLE PRECISION X1(50),FIB(50) DOUBLE PRECISION D1,D2,R1,R2,F1,F2 DOUBLE PRECISION EPS,EPSR C COMMON /PREC/EPS COMMON /NOF/NBL C NBL=0 FIB(1)=1.0 FIB(2)=1.0 C DO 1 I=3,50 1 FIB(I)=FIB(I-1)+FIB(I-2) C EPSR=1.0/EPS C I=2 2 IF(FIB(I).GE.EPSR) GO TO 3 I=I+1 GO TO 2 C 3 D1=0.0 D2=FIB(I) C 4 R1=D1+FIB(I-2) R2=D1+FIB(I-1) C R11=R1*EPS R22=R2*EPS C DO 5 J=1,N 5 X1(J)=X(J)+R11*S(J) CALL FUNC(N,X1,F1,nexp) NBL=NBL+1 C DO 6 J=1,N 6 X1(J)=X(J)+R22*S(J) CALL FUNC(N,X1,F2,nexp) NBL=NBL+1 C IF(F2.GT.F1) THEN D2=R2 ELSE D1=R1 ENDIF C IF(R1.EQ.R2) GO TO 7 I=I-1 GO TO 4 C 7 BL=R11 C RETURN END C C***************************************Last card of L2.FOR C********************************************************* C C Subroutine for one dimensional searching. C Quadratic interpolation of Powell technique. C ======================= C********************************************************* C Neculai Andrei, C Research Institute for Informatics C Advanced Modeling and Optimization Laboratory C********************************************************* C SUBROUTINE L3(N,X,S,BL,nexp) DOUBLE PRECISION X(N),S(N),BL DOUBLE PRECISION V(3),P(3),XT(50) DOUBLE PRECISION H,HL,F0,E,F,D,R,DM,DMM,BLM,EPS C COMMON /PREC/EPS COMMON /L3A/ICOD COMMON /NOF/NBL C NBL=0 ICOD=0 ITER=0 H=2.0 HL=0.5 C CALL FUNC(N,X,F0,nexp) NBL=NBL+1 P(1)=0.0 V(1)=F0 C DO 1 I=1,N 1 XT(I)=X(I)+HL*S(I) CALL FUNC(N,XT,F0,nexp) NBL=NBL+1 P(2)=HL V(2)=F0 C IF(V(1).LT.V(2)) THEN DO 2 I=1,N 2 XT(I)=X(I)-HL*S(I) CALL FUNC(N,XT,F0,nexp) NBL=NBL+1 P(3)=-HL V(3)=F0 C ELSE C DO 3 I=1,N 3 XT(I)=X(I)+2.*HL*S(I) CALL FUNC(N,XT,F0,nexp) NBL=NBL+1 P(3)=2.*HL V(3)=F0 ENDIF C 4 CONTINUE C ITER=ITER+1 IF(ITER.GT.1000) THEN ICOD=1 BL=0.d0 RETURN ENDIF C E=(P(2)*P(2)-P(3)*P(3))*V(1) + + (P(3)*P(3)-P(1)*P(1))*V(2) + + (P(1)*P(1)-P(2)*P(2))*V(3) C F=(P(2)-P(3))*V(1)+(P(3)-P(1))*V(2)+(P(1)-P(2))*V(3) IF(F.EQ.0.0) THEN BL=1.0 GO TO 10 ENDIF C BLM=(0.5*E)/F C D=(P(1)-P(2))*(P(2)-P(3))*(P(3)-P(1)) C R=F/D C DM=DABS(BLM-P(1)) DMM=DM IMI=1 IMA=1 DO 5 I=2,3 D=DABS(BLM-P(I)) IF(DM.GT.D) THEN DM=D IMI=I ENDIF IF(DMM.LT.D) THEN DMM=D IMA=I ENDIF 5 CONTINUE C IF(DM.GT.H) THEN IF(R.LT.0.0) THEN IF(BLM.GT.P(IMI)) THEN P(IMA)=P(IMI)+H ELSE P(IMA)=P(IMI)-H ENDIF C ELSE C IF(BLM.LT.P(IMA)) THEN P(IMA)=P(IMA)-H ELSE P(IMA)=P(IMA)+H ENDIF C ENDIF C DO 6 I=1,N 6 XT(I)=X(I)+P(IMA)*S(I) CALL FUNC(N,XT,F0,nexp) NBL=NBL+1 V(IMA)=F0 GO TO 4 C ENDIF C DO 7 I=1,N 7 XT(I)=X(I)+BLM*S(I) CALL FUNC(N,XT,F0,nexp) NBL=NBL+1 C E=V(1) F=E IMI=1 IMA=1 DO 8 I=2,3 IF(E.GT.V(I)) THEN E=V(I) IMI=I ENDIF IF(F.LT.V(I)) THEN F=V(I) IMA=I ENDIF 8 CONTINUE C IF(R.LT.0.0.AND.DM.LT.EPS) THEN BL=BLM IF(E.LT.F0) BL=P(IMI) GO TO 10 ENDIF C IF(R.LT.0.0.AND.(DM.GT.EPS.OR.DM.LT.H)) THEN P(IMA)=BLM V(IMA)=F0 GO TO 4 ENDIF C 10 CONTINUE C RETURN END C C************************************ Last card of L3.FOR C****************************************************** C C Dichotomous search technique for C one-dimensional minimization. C C****************************************************** C Neculai Andrei, C Research Institute for Informatics C Advanced Modeling and Optimization Laboratory C****************************************************** C SUBROUTINE L4(N,X,S,BL,nexp) DOUBLE PRECISION X(N),S(N),BL DOUBLE PRECISION EPS,A,B,C,C1,C2 DOUBLE PRECISION XT(40),F1,F2 C COMMON /PREC/EPS COMMON /NOF/NBL C NBL=0 EPS1=EPS/10 A=0.0 B=1.0 C 1 C=A+(B-A)/2. C C1=C+EPS1/2. DO 2 I=1,N 2 XT(I)=X(I)+C1*S(I) CALL FUNC(N,XT,F1,nexp) NBL=NBL+1 C C2=C-EPS1/2. DO 3 I=1,N 3 XT(I)=X(I)+C2*S(I) CALL FUNC(N,XT,F2,nexp) NBL=NBL+1 C IF(F2.LT.F1) THEN B=C1 ELSE A=C2 ENDIF C IF(DABS(B-A).LE.EPS) GO TO 4 GO TO 1 4 BL=C C RETURN END C C*******************************Last card of L4.FOR * * * .-----------------------------------------------------------------. * | Subroutine inipoint | * | =================== | * | Final | * |Subroutine for initial point specification. | * |This is an user subroutine: | * | | * | The calling sequence is: | * | | * | call inipoint(n,x,nexp) | * |where: | * | n (integer) the number of variables, | * | x (double) array with the initial point. | * | nexp (integer) parameter specifying the number of the | * | problem considered in a train of experiments. | * | | * | Neculai Andrei | * .-----------------------------------------------------------------. ****************************************************************** subroutine inipoint(n,x, nexp) C This subroutine computes the initial point real*8 x(40) go to ( 1, 2, 3, 4, 5, 6, 7, 8, 9,10, * 11,12,13,14,15,16,17,18,19,20, * 21,22,23,24,25,26,27,28,29,30, * 31,32,33,34,35,36,37,38,39,40, * 41,42,43,44,45,46,47,48,49,50, * 51,52,53,54,55,56,57,58,59,60, * 61,62,63,64,65,66,67,68,69,70, * 71,72,73,74,75,76,77,78) nexp 1 continue c Freudenstein & Roth - FREUROTH (CUTE) i=1 991 x(i) = 0.5d0 x(i+1)= -2.d0 i=i+2 if(i.le.n) go to 991 return 2 continue c Trigonometric do i=1,n x(i) = 0.2d0 end do return 3 continue c Extended Rosenbrock SROSENBR (CUTE) i=1 993 x(i) = -1.2d0 x(i+1)= 1.d0 i=i+2 if(i.le.n) go to 993 return 4 continue c Extended White & Holst i=1 994 x(i) = -1.2d0 x(i+1)= 1.d0 i=i+2 if(i.le.n) go to 994 return 5 continue c Extended Beale i=1 995 x(i) = 1.d0 x(i+1)= 0.8d0 i=i+2 if(i.le.n) go to 995 return 6 continue c Penalty do i=1,n x(i) = float(i) end do return 7 continue c Perturbed Quadratic do i=1,n x(i) = 0.5d0 end do return 8 continue c Raydan 1 do i=1,n x(i) = 1.d0 end do return 9 continue c Raydan 2 do i=1,n x(i) = 1.d0 end do return 10 continue c Diagonal 1 do i=1,n x(i) = 10.d0 end do return 11 continue c Diagonal 2 do i=1,n x(i) = 1.d0/float(i) end do return 12 continue c Diagonal 3 do i=1,n x(i) = 1.d0 end do return 13 continue c Hager do i=1,n x(i) = 1.d0 end do return 14 continue c Generalized Tridiagonal 1 do i=1,n x(i) = 2.d0 end do return 15 continue c Extended Tridiagonal 1 do i=1,n x(i) = 2.d0 end do return 16 continue c Extended Three Expo Terms do i=1,n x(i) = 0.1d0 end do return 17 continue c Generalized Tridiagonal 2 do i=1,n x(i) = -1.d0 end do return 18 continue c Diagonal 4 do i=1,n x(i) = 1.d0 end do return 19 continue c Diagonal 5 do i=1,n x(i) = 1.1d0 end do return 20 continue c Extended Himmelblau do i=1,n x(i) = 1.d0 end do return 21 continue c Generalized PSC1 i=1 9921 x(i) = 3.d0 x(i+1)= 0.1d0 i=i+2 if(i.le.n) go to 9921 return 22 continue c Extended PSC1 i=1 9922 x(i) = 3.d0 x(i+1)= 0.1d0 i=i+2 if(i.le.n) go to 9922 return 23 continue c Extended Powell i=1 9923 x(i) = 3.d0 x(i+1)= -1.d0 x(i+2)= 0.d0 x(i+3)= 1.d0 i=i+4 if(i.le.n) go to 9923 return 24 continue c Extended BD1 do i=1,n x(i) = 0.1d0 end do return 25 continue c Extended Maratos i=1 9925 x(i) = 1.1d0 x(i+1)= 0.1d0 i=i+2 if(i.le.n) go to 9925 return 26 continue c Extended Cliff i=1 9926 x(i) = 0.d0 x(i+1)= -0.1d0 i=i+2 if(i.le.n) go to 9926 return 27 continue c Quadratic Diagonal Perturbed do i=1,n x(i) = 0.5d0 end do return 28 continue c Extended Wood WOODS (CUTE) i=1 9928 x(i) = -3.d0 x(i+1)= -1.d0 x(i+2)= -3.d0 x(i+3)= -1.d0 i=i+4 if(i.le.n) go to 9928 return 29 continue c Extended Hiebert do i=1,n x(i) = 0.0001d0 end do return 30 continue c Quadratic QF1 do i=1,n x(i) = 1.d0 end do return 31 continue c Extended Quadratic Penalty QP1 do i=1,n x(i) = 1.d0 end do return 32 continue c Extended Quadratic Penalty QP2 do i=1,n x(i) = 1.d0 end do return 33 continue c Quadratic QF2 do i=1,n x(i) = 0.5d0 end do return 34 continue c Extended EP1 do i=1,n x(i) = 1.5d0 end do return 35 continue c Extended Tridiagonal 2 do i=1,n x(i) = 1.d0 end do return 36 continue c BDQRTIC (CUTE) do i=1,n x(i) = 1.d0 end do return 37 continue c TRIDIA (CUTE) do i=1,n x(i) = 1.d0 end do return 38 continue c ARWHEAD (CUTE) do i=1,n x(i) = 1.d0 end do return 39 continue c NONDIA (CUTE) do i=1,n x(i) = -1.d0 end do return 40 continue c NONDQUAR (CUTE) i=1 9940 x(i) = 1.d0 x(i+1)= -1.d0 i=i+2 if(i.le.n) go to 9940 return 41 continue c DQDRTIC (CUTE) do i=1,n x(i) = 3.d0 end do return 42 continue c EG2 (CUTE) do i=1,n x(i) = 1.d0 end do return 43 continue c DIXMAANA (CUTE) do i=1,n x(i) = 2.d0 end do return 44 continue c DIXMAANB (CUTE) do i=1,n x(i) = 2.d0 end do return 45 continue c DIXMAANC (CUTE) do i=1,n x(i) = 2.d0 end do return 46 continue c DIXMAANE (CUTE) do i=1,n x(i) = 2.d0 end do return 47 continue c Partial Perturbed Quadratic do i=1,n x(i) = 0.5d0 end do return 48 continue c Broyden Tridiagonal do i=1,n x(i) = -1.d0 end do return 49 continue c Almost Perturbed Quadratic do i=1,n x(i) = 0.5d0 end do return 50 continue c Tridiagonal Perturbed Quadratic do i=1,n x(i) = 0.5d0 end do return 51 continue c EDENSCH (CUTE) do i=1,n x(i) = 0.0d0 end do return 52 continue c VARDIM (CUTE) do i=1,n x(i) = 1.d0- float(i)/float(n) end do return 53 continue c STAIRCASE S1 do i=1,n x(i) = 1.0d0 end do return 54 continue c LIARWHD (CUTE) do i=1,n x(i) = 4.d0 end do return 55 continue c DIAGONAL 6 do i=1,n x(i) = 1.d0 end do return 56 continue c DIXON3DQ (CUTE) do i=1,n x(i) = -1.d0 end do return 57 continue c ENGVAL1 (CUTE) do i=1,n x(i) = 2.d0 end do return 58 continue c DENSCHNA (CUTE) do i=1,n x(i) = 2.d0 end do return 59 continue c DENSCHNC (CUTE) do i=1,n x(i) = 2.d0 end do return 60 continue c DENSCHNB (CUTE) do i=1,n x(i) = 1.d0 end do return 61 continue c DENSCHNF (CUTE) i=1 9961 x(i) = 2.d0 x(i+1)= 0.d0 i=i+2 if(i.le.n) go to 9961 return 62 continue c SINQUAD (CUTE) do i=1,n x(i) = 0.1d0 end do return 63 continue c BIGGSB1 (CUTE) do i=1,n x(i) = 0.d0 end do return 64 continue c Extended Block-Diagonal BD2 i=1 9964 x(i) = 1.5d0 x(i+1)= 2.d0 i=i+2 if(i.le.n) go to 9964 return 65 continue c Generalized quartic GQ1 do i=1,n x(i) = 1.d0 end do return 66 continue c DIAGONAL 7 do i=1,n x(i) = 0.d0 end do return 67 continue c DIAGONAL 8 do i=1,n x(i) = 1.d0 end do return 68 continue c Full Hessian do i=1,n x(i) = 1.d0 end do return 69 continue c SINCOS i=1 9969 x(i) = 3.d0 x(i+1)= 0.1d0 i=i+2 if(i.le.n) go to 9969 return 70 continue c Generalized quartic GQ2 i=1 9970 x(i) = -1.2d0 x(i+1)= 1.d0 i=i+2 if(i.le.n) go to 9970 return 71 continue c EXTROSNB (CUTE) i=1 9971 x(i) = -1.2d0 x(i+1)= 1.d0 i=i+2 if(i.le.n) go to 9971 return 72 continue c ARGLINB (CUTE) i=1 9972 x(i) = 0.01d0 x(i+1)= 0.001d0 i=i+2 if(i.le.n) go to 9972 return 73 continue c FLETCHCR (CUTE) do i=1,n x(i) = 0.5d0 end do return 74 continue c HIMMELBG (CUTE) do i=1,n x(i) = 1.5d0 end do return 75 continue c HIMMELBH (CUTE) i=1 9975 x(i) = 0.0d0 x(i+1)= 2.0d0 i=i+2 if(i.le.n) go to 9975 return 76 continue c HIMMEL x(1)= 1.d0 x(2)= -0.5d0 77 continue c PROPAN x(1)=0.01d0 x(2)=30.d0 x(3)=0.05d0 x(4)=0.5d0 x(5)=0.05d0 return 78 continue c REACTOR x(1)=1.09d0 x(2)=1.05d0 x(3)=3.05d0 x(4)=0.99d0 x(5)=6.05d0 x(6)=1.09d0 return end c------------------------------------------------ End INIPOINT ************************************************************** * Date created: October 28, 2004 * * Final * * TEST FUNCTIONS FOR UNCONSTRAINED OPTIMIZATION * =============================================== * * 57 problems: October 28, 2004 * 66 problems: March 29, 2005 * 70 problems: April 28, 2005 * 75 problems: April 19, 2006 * Neculai Andrei ************************************************************** *cg_value(f,x,n,nexp) subroutine func(n,x,f, nexp) real*8 x(n), f real*8 t1,t2,t3,t4, c, d real*8 s, temp(500000), temp1, sum real*8 u(500000), v(500000), t(500000) real*8 u1, v1, c1, c2 real*8 alpha, beta, gamma, delta integer k1, k2, k3, k4 real*8 r5, r6, r7, r8, r9, r10, r, fx(20) real*8 k1d,k2d,k3d,r1, r2 * go to ( 1, 2, 3, 4, 5, 6, 7, 8, 9,10, * 11,12,13,14,15,16,17,18,19,20, * 21,22,23,24,25,26,27,28,29,30, * 31,32,33,34,35,36,37,38,39,40, * 41,42,43,44,45,46,47,48,49,50, * 51,52,53,54,55,56,57,58,59,60, * 61,62,63,64,65,66,67,68,69,70, * 71,72,73,74,75,76,77,78) nexp cF1 FREUROTH (CUTE) * Extended Freudenstein & Roth * * Initial Point: [0.5, -2, ...,0.5, -2]. * 1 continue f = 0.d0 do i=1,n/2 t1=-13.d0+x(2*i-1)+5.d0*x(2*i)*x(2*i)-x(2*i)**3-2.d0*x(2*i) t2=-29.d0+x(2*i-1)+x(2*i)**3+x(2*i)**2-14.d0*x(2*i) f = f + t1*t1 + t2*t2 end do return c******************************************************************** cF2 Extended Trigonometric Function * * Initial Point: [0.2, 0.2, ....,0.2]. 2 continue s= float(n) do i=1,n s = s - dcos(x(i)) end do do i=1,n temp(i) = s + float(i)*(1.d0-dcos(x(i))) - dsin(x(i)) end do f = 0.d0 do i=1,n f = f + temp(i)**2 end do return cF3 SROSENBR (CUTE) * Extended Rosenbrock function * * Initial point: [-1.2, 1, -1.2, 1, ..........., -1.2, 1] 3 continue c=100.d0 f=0.d0 do i=1,n/2 f = f + c*(x(2*i)-x(2*i-1)**2)**2 + (1.d0-x(2*i-1))**2 end do return c****************************************************************** cF4 Extended White & Holst function * * Initial point: [-1.2, 1, -1.2, 1, ..........., -1.2, 1] 4 continue c=100.d0 f=0.d0 do i=1,n/2 f = f + c*(x(2*i)-x(2*i-1)**3)**2 + (1.d0-x(2*i-1))**2 end do return c****************************************************************** * cF5 Extended Beale Function BEALE (CUTE) * * Initial Point: [1, 0.8, ...., 1, 0.8] *-------------------------------------------------------------- * 5 continue * f=0.d0 do i=1,n/2 t1=1.5d0 -x(2*i-1)+x(2*i-1)*x(2*i) t2=2.25d0 -x(2*i-1)+x(2*i-1)*x(2*i)*x(2*i) t3=2.625d0-x(2*i-1)+x(2*i-1)*x(2*i)*x(2*i)*x(2*i) f = f + t1*t1 + t2*t2 + t3*t3 end do * return cF6 Extended Penalty Function U52 (MatrixRom) * * Intial Point: [1,2,3,.....,n]. * 6 continue temp1=0.d0 do i=1,n temp1 = temp1 + x(i)**2 end do * f = (temp1 - 0.25d0)**2 do i=1,n-1 f = f + (x(i)-1.d0)**2 end do return cF7 Perturbed Quadratic function * * Initial Point: [0.5, 0.5, ......, 0.5]. * 7 continue temp1 = 0.d0 do i=1,n temp1 = temp1 + x(i) end do f = temp1*temp1/100.d0 do i=1,n f = f + float(i)*x(i)**2 end do return cF8 Raydan 1 Function * * Initial point: [1, 1, ..., 1] 8 continue f=0.d0 do i=1,n f = f + float(i) * (dexp(x(i))-x(i)) / 10.d0 end do return cF9 Raydan 2 Function * * Initial Point: [1, 1, .....,1] * 9 continue f=0.d0 do i=1,n f = f + dexp(x(i)) - x(i) end do return cF10 Diagonal1 Function * * Initial Point: [1/n, 1/n, ....., 1/n] * 10 continue f=0.d0 do i=1,n f = f + dexp(x(i)) - x(i)*float(i) end do return cF11 Diagonal2 Function * * Initial Point: [1/1, 1/2, ....., 1/n] * 11 continue f=0.d0 do i=1,n f = f + dexp(x(i)) - x(i)/float(i) end do return cF12 Diagonal3 Function * * Initial Point: [1,1,...,1] 12 continue f=0.d0 do i=1,n f = f + dexp(x(i)) - float(i)*dsin(x(i)) end do return cF13 Hager Function * * Initial Point: [1,1,...,1] 13 continue f=0.d0 do i=1,n f = f + dexp(x(i)) - x(i)*sqrt(float(i)) end do return cF14 Generalized Tridiagonal-1 Function * * Initial Point: [2,2,...,2] 14 continue f=0.d0 do i=1,n-1 u(i) = x(i) + x(i+1) - 3.d0 v(i) = x(i) - x(i+1) + 1.d0 end do do i=1,n-1 f = f + u(i)**2 + v(i)**4 end do return cF15 Extended Tridiagonal-1 Function * * Initial Point: [2,2,...,2] 15 continue f=0.d0 do i=1,n/2 u(i) = x(2*i-1) + x(2*i) - 3.d0 v(i) = x(2*i-1) - x(2*i) + 1.d0 end do do i=1,n/2 f = f + u(i)**2 + v(i)**4 end do return cF16 Extended Three Exponential Terms * * Intial Point: [0.1,0.1,.....,0.1]. * 16 continue f=0.d0 j=1 do i=1,n/2 t1=x(j) + 3.d0*x(j+1) - 0.1d0 t2=x(j) - 3.d0*x(j+1) - 0.1d0 t3=-x(j) - 0.1d0 f = f + dexp(t1) + dexp(t2) + dexp(t3) end do return cF17 Generalized Tridiagonal-2 * * Initial point: [-1, -1, .........., -1., -1] 17 continue f = 0.d0 u(1) = 5.d0*x(1)-3.d0*x(1)**2-x(1)**3-3.d0*x(2)+1.d0 do i=2,n-1 u(i)=5.d0*x(i)-3.d0*x(i)**2-x(i)**3-x(i-1)-3.d0*x(i+1)+1.d0 end do u(n)=5.d0*x(n)-3.d0*x(n)**2-x(n)**3-x(n-1)+1.d0 * do i=1,n f = f + u(i)**2 end do return cF18 Diagonal4 Function * * Initial point: [1, 1, .........., 1., 1] 18 continue c=100.d0 f = 0.d0 do i=1,n/2 f = f + (x(2*i-1)**2 + c*x(2*i)**2)/2.d0 end do return cF19 Diagonal5 Function (MatrixRom) * * Initial point: [1.1, 1.1, .........., 1.1] 19 continue f = 0.d0 do i=1,n f = f + dlog( dexp(x(i)) + dexp(-x(i)) ) end do return cF20 HIMMELBC (CUTE) * Extended Himmelblau Function * * Initial Point: [1, 1, ....., 1] 20 continue f=0.d0 do i=1,n/2 u1 = x(2*i-1)**2 + x(2*i) - 11.d0 v1 = x(2*i-1) + x(2*i)**2 - 7.d0 f = f + u1*u1 + v1*v1 end do return cF21 Generalized PSC1 Function * * Initial point: [3, 0.1, ..., 3, 0.1] 21 continue f = 0.d0 do i=1,n-1 f = f + (x(i)**2 +x(i+1)**2 +x(i)*x(i+1))**2 * + (dsin(x(i)))**2 + (dcos(x(i+1)))**2 end do return cF22 Extended PSC1 Function * * Initial point: [3, 0.1, ..., 3, 0.1] 22 continue f = 0.d0 do i=1,n/2 f = f + (x(2*i-1)**2 +x(2*i)**2 +x(2*i-1)*x(2*i))**2 * + (dsin(x(2*i-1)))**2 + (dcos(x(2*i)))**2 end do return cF23 Extended Powell * * Initial Point: [3, -1, 0, 1, ......]. * 23 continue f=0.d0 do i=1,n/4 t1= x(4*i-3) + 10.d0*x(4*i-2) t2= x(4*i-1) - x(4*i) t3= x(4*i-2) - 2.d0*x(4*i-1) t4= x(4*i-3) - x(4*i) f = f + t1*t1 + 5.d0*t2*t2 + t3**4 + 10.d0*t4**4 end do return cF24 Extended Block Diagonal BD1 Function * * Initial Point: [0.1, 0.1, ..., 0.1]. *------------------------------------------------------------------- * 24 continue f = 0.d0 do i=1,n/2 t1 = x(2*i-1)**2 + x(2*i)**2 - 2.d0 t2 = dexp(x(2*i-1)) - x(2*i) f = f + t1*t1 + t2*t2 end do return cF25 Extended Maratos Function * * Initial Point: [1.1, 0.1, ...,1.1, 0.1]. *------------------------------------------------------------------- * 25 continue c = 100.d0 f = 0.d0 do i=1,n/2 t1 = x(2*i-1)**2 + x(2*i)**2 - 1.d0 f = f + (x(2*i-1) + c*t1*t1) end do return cF26 Extended Cliff CLIFF (CUTE) * * Initial Point: [0, -1, ......, 0, -1]. * 26 continue f=0.d0 do i=1,n/2 temp1 = (x(2*i-1)-3.d0)/100.d0 f = f+temp1*temp1-(x(2*i-1)-x(2*i))+dexp(20.d0*(x(2*i-1)-x(2*i))) end do return cF27 Quadratic Diagonal Perturbed Function * * Initial Point: [0.5, 0.5, ......, 0.5]. * 27 continue temp1 = 0.d0 do i=1,n temp1 = temp1 + x(i) end do f = temp1*temp1 do i=1,n f = f + (float(i)/100.d0) * x(i)**2 end do return cF28 Extended Wood Function * WOODS (CUTE) * * Initial Point: [-3,-1,-3,-1,......] * 28 continue f=0.d0 do i=1,n/4 f = f + 100.d0*(x(4*i-3)**2-x(4*i-2))**2 * + (x(4*i-3)-1.d0)**2 * + 90.d0*(x(4*i-1)**2-x(4*i))**2 * + (1.d0-x(4*i-1))**2 * + 10.1d0*(x(4*i-2)-1.d0)**2 * + 10.1d0*(x(4*i) -1.d0)**2 * + 19.8d0*(x(4*i-2)-1.d0)*(x(4*i)-1.d0) end do return cF29 Extended Hiebert Function * * Initial Point: [0,0,...0]. 29 continue c1=10.d0 c2=500.d0 f = 0.d0 do i=1,n/2 f = f + (x(2*i-1)-c1)**2 + (x(2*i-1)*x(2*i)-c2)**2 end do return cF30 Quadratic Function QF1 * * Initial Point: [1,1,...1]. 30 continue f = 0.d0 do i=1,n f = f + float(i)*x(i)*x(i) end do f = f/2.d0 f = f - x(n) return cF31 Extended Quadratic Penalty QP1 Function * * Initial Point: [1, 1, ......,1]. * 31 continue t1=0.d0 do i=1,n t1 = t1 + x(i)*x(i) end do t1 = t1 - 0.5d0 f = 0.d0 do i=1,n-1 f = f + (x(i)*x(i) - 2.d0)**2 end do f = f + t1*t1 return cF32 Extended Quadratic Penalty QP2 Function * * Initial Point: [1, 1, ......,1]. * 32 continue t1=0.d0 do i=1,n t1 = t1 + x(i)*x(i) end do t1 = t1 - 100.d0 f = 0.d0 do i=1,n-1 f = f + (x(i)*x(i) - dsin(x(i)))**2 end do f = f + t1*t1 return cF33 A Quadratic Function QF2 * * Initial Point: [0.5, 0.5, ......,0.5]. * 33 continue f=0.d0 do i=1,n f = f + float(i)*(x(i)**2 - 1.d0)**2 end do f = f/2.d0 f = f - x(n) return cF34 Extended EP1 Function * * Initial Point: [1.5.,1.5.,...,1.5]. *------------------------------------------------------------------- * 34 continue * f=0.d0 do i=1,n/2 t1=dexp(x(2*i-1)-x(2*i)) - 5.d0 t2=x(2*i-1)-x(2*i) t3=x(2*i-1)-x(2*i)-11.d0 f = f + t1*t1 + (t2*t2)*(t3*t3) end do return cF35 Extended Tridiagonal-2 Function * * Initial Point: [1.,1.,...,1.]. *------------------------------------------------------------------- * 35 continue * c=0.1d0 f=0.d0 do i=1,n-1 f = f + (x(i)*x(i+1)-1.d0)**2 + c*(x(i)+1.d0)*(x(i+1)+1.d0) end do return cF36 BDQRTIC (CUTE) * * Initial point x0=[1.,1.,...,1.]. * * 36 continue * n4=n-4 f=0.d0 do i=1,n4 temp(i) = x(i)**2 + 2.d0*x(i+1)**2 + 3.d0*x(i+2)**2 * + 4.d0*x(i+3)**2 + 5.d0*x(n)**2 end do do i=1,n4 f = f + (-4.d0*x(i)+3.d0)**2 + temp(i)**2 end do return cF37 TRIDIA (CUTE) * * Initial point x0=[1,1,...,1]. * * 37 continue * alpha=2.d0 beta=1.d0 gamma=1.d0 delta=1.d0 f=gamma*(delta*x(1)-1.d0)**2 do i=2,n f = f + float(i)*(alpha*x(i)-beta*x(i-1))**2 end do return cF38 ARWHEAD (CUTE) * * Initial point x0=[1,1,...,1]. * * 38 continue f=0.d0 do i=1,n-1 f = f + (-4.d0*x(i)+3.d0) + (x(i)**2+x(n)**2)**2 end do return cF39 * NONDIA (Shanno-78) (CUTE) * * Initial point x0=[-1,-1,...,-1]. * * 39 continue c=100.d0 f=(x(1)-1.d0)**2 + c*(x(1)-x(1)**2)**2 do i=2,n f = f + c*(x(1)-x(i)**2)**2 end do return cF40 NONDQUAR (CUTE) * * Initial point x0=[1,-1,1,-1,...,]. * * 40 continue f = (x(1)-x(2))**2 + (x(n-1)+x(n))**2 do i=1,n-2 f = f + (x(i)+x(i+1)+x(n))**4 end do return cF41 DQDRTIC * * Initial point x0=[3,3,3...,3]. * * 41 continue c=100.d0 d=100.d0 f=0.d0 do i=1,n-2 f = f + (x(i)**2 + c*x(i+1)**2 + d*x(i+2)**2) end do return cF42 EG2 (CUTE) * * Initial point x0=[1,1,1...,1]. * * 42 continue f=0.5d0*dsin(x(n)*x(n)) do i=1,n-1 f = f + dsin(x(1)+x(i)*x(i)-1.d0) end do return cF43 DIXMAANA (CUTE) * * Initial point x0=[2.,2.,2...,2.]. * 43 continue * alpha = 1.d0 beta = 0.d0 gamma = 0.125d0 delta = 0.125d0 k1 = 0 k2 = 0 k3 = 0 k4 = 0 m = n/3 f = 1.d0 do i=1,n f = f + alpha * x(i)* x(i) * (float(i)/float(n))**k1 end do do i=1,n-1 f = f + beta*x(i)*x(i)*((x(i+1)+x(i+1)*x(i+1))**2) * * (float(i)/float(n))**k2 end do do i=1,2*m f = f + gamma * x(i)*x(i) * (x(i+m)**4) * * (float(i)/float(n))**k3 end do do i=1,m f = f + delta * x(i) * x(i+2*m) * * (float(i)/float(n))**k4 end do return cF44 DIXMAANB (CUTE) * * Initial point x0=[2.,2.,2...,2.]. * 44 continue * alpha = 1.d0 beta = 0.0625d0 gamma = 0.0625d0 delta = 0.0625d0 k1 = 0 k2 = 0 k3 = 0 k4 = 0 m = n/3 f = 1.d0 do i=1,n f = f + alpha * x(i)* x(i) * (float(i)/float(n))**k1 end do do i=1,n-1 f = f + beta*x(i)*x(i)*((x(i+1)+x(i+1)*x(i+1))**2) * * (float(i)/float(n))**k2 end do do i=1,2*m f = f + gamma * x(i)*x(i) * (x(i+m)**4) * * (float(i)/float(n))**k3 end do do i=1,m f = f + delta * x(i) * x(i+2*m) * * (float(i)/float(n))**k4 end do return cF45 DIXMAANC (CUTE) * * Initial point x0=[2.,2.,2...,2.]. * 45 continue * alpha = 1.d0 beta = 0.125d0 gamma = 0.125d0 delta = 0.125d0 k1 = 0 k2 = 0 k3 = 0 k4 = 0 m = n/3 f = 1.d0 do i=1,n f = f + alpha * x(i)* x(i) * (float(i)/float(n))**k1 end do do i=1,n-1 f = f + beta*x(i)*x(i)*((x(i+1)+x(i+1)*x(i+1))**2) * * (float(i)/float(n))**k2 end do do i=1,2*m f = f + gamma * x(i)*x(i) * (x(i+m)**4) * * (float(i)/float(n))**k3 end do do i=1,m f = f + delta * x(i) * x(i+2*m) * * (float(i)/float(n))**k4 end do return cF46 DIXMAANE (CUTE) * * Initial point x0=[2.,2.,2...,2.]. * 46 continue * alpha = 1.d0 beta = 0.d0 gamma = 0.125d0 delta = 0.125d0 k1 = 1 k2 = 0 k3 = 0 k4 = 1 m = n/3 f = 1.d0 do i=1,n f = f + alpha * x(i)* x(i) * (float(i)/float(n))**k1 end do do i=1,n-1 f = f + beta*x(i)*x(i)*((x(i+1)+x(i+1)*x(i+1))**2) * * (float(i)/float(n))**k2 end do do i=1,2*m f = f + gamma * x(i)*x(i) * (x(i+m)**4) * * (float(i)/float(n))**k3 end do do i=1,m f = f + delta * x(i) * x(i+2*m) * * (float(i)/float(n))**k4 end do return cF47 Partial Perturbed Quadratic * * Initial point x0=[0.5, 0.5, ..., 0.5]. * 47 continue * temp(1) = x(1) + x(2) do i=2,n-1 temp(i) = temp(i-1) + x(i+1) end do f=x(1)*x(1) do i=2,n f = f + float(i)*x(i)*x(i) + temp(i-1)*temp(i-1)/100.d0 end do return cF48 Broyden Tridiagonal * * Initial point x0=[-1., -1., ..., -1.]. * 48 continue * temp(1) = 3.d0*x(1) - 2.d0*x(1)*x(1) do i=2,n-1 temp(i) = 3.d0*x(i)-2.d0*x(i)*x(i)-x(i-1)-2.d0*x(i+1)+1.d0 end do temp(n) = 3.d0*x(n)-2.d0*x(n)*x(n)-x(n-1)+1.d0 f = 0.d0 do i=1,n f = f + temp(i)*temp(i) end do return cF49 Almost Perturbed Quadratic * * Initial point x0=[0.5, 0.5, ...,0.5]. * 49 continue * f = 0.01d0*(x(1)+x(n))**2 do i=1,n f = f + float(i)*x(i)*x(i) end do return cF50 Tridiagonal Perturbed Quadratic * * Initial point x0=[0.5, 0.5, ...,0.5]. * 50 continue * do i=1,n-2 temp(i) = x(i) + x(i+1) + x(i+2) end do f = x(1)*x(1) do i=2,n-1 f = f + float(i)*x(i)*x(i) + temp(i-1)**2 end do return cF51 EDENSCH Function (CUTE) * * Initial Point: [0., 0., ..., 0.]. 51 continue f = 16.d0 do i=1,n-1 f = f + (x(i)-2.d0)**4 + * (x(i)*x(i+1)-2.d0*x(i+1))**2 + * (x(i+1)+1.d0)**2 end do return cF52 VARDIM Function (CUTE) * * Initial Point: [1-1/n, 1-2/n, ..., 1-n/n]. 52 continue sum = 0.d0 do i=1,n sum = sum + float(i)*x(i) end do sum = sum - float(n)*float(n+1)/2.d0 f = sum**2 + sum**4 do i=1,n f = f + (x(i)-1.d0)**2 end do return cF53 STAIRCASE S1 * * Initial point x0=[1,1,...,1]. * 53 continue f=0.d0 do i=1,n-1 f = f + (x(i)+x(i+1)-float(i))**2 end do return cF54 LIARWHD (CUTE) * * Initial point x0=[4., 4., ....4.]. * 54 continue f=0.d0 do i=1,n f = f + 4.d0*(x(i)*x(i) - x(1))**2 + (x(i)-1.d0)**2 end do return cF55 DIAGONAL 6 * * Initial point x0=[1.,1., ..., 1.]. * 55 continue f = 0.d0 do i=1,n f = f + dexp(x(i)) - (1.d0+x(i)) end do return cF56 DIXON3DQ (CUTE) * * Initial Point x0=[-1, -1,..., -1] March 7, 2005 * 56 continue f=(x(1)-2.d0)**2 do i=1,n-1 f = f + (x(i)-x(i+1))**2 end do f = f + (x(n)-1.d0)**2 return cF57 ENGVAL1 (CUTE) * * Initial point x0=[2.,2.,2...,2.]. * 57 continue do i=1,n-1 t(i) = x(i)*x(i) + x(i+1)*x(i+1) end do f = 0.d0 do i=1,n-1 f = f + t(i)*t(i) + (-4.d0*x(i) + 3.d0) end do return cF58 DENSCHNA (CUTE) * * Initial point: [2,2,...,2] * 58 continue f=0.d0 do i=1,n/2 f = f + x(2*i-1)**4 + * (x(2*i-1)+x(2*i))**2 + * (-1.d0+dexp(x(2*i)))**2 end do return cF59 DENSCHNC (CUTE) * * Initial point: [2,2,...,2] 59 continue f=0.d0 do i=1,n/2 f = f + (-2.d0+x(2*i-1)**2+x(2*i)**2)**2 + * (-2.d0+dexp(x(2*i-1)-1.d0)+x(2*i)**3)**2 end do return cF60 DENSCHNB (CUTE) * * Initial point: [1,1,...,1] 60 continue f=0.d0 do i=1,n/2 f = f + (x(2*i-1)-2.d0)**2 + * ((x(2*i-1)-2.d0)**2)*(x(2*i)**2) + * (x(2*i)+1.d0)**2 end do return cF61 DENSCHNF (CUTE) * * Initial point: [2,0,2,0,...,2,0] 61 continue f=0.d0 do i=1,n/2 f=f+(2.d0*(x(2*i-1)+x(2*i))**2+(x(2*i-1)-x(2*i))**2-8.d0)**2+ * (5.d0*x(2*i-1)**2+(x(2*i)-3.d0)**2-9.d0)**2 end do return cF62 SINQUAD (CUTE) * * Initial Point: [0.1, 0.1, ..., 0.1] 62 continue f=(x(1)-1.d0)**4 + (x(n)**2-x(1)**2)**2 do i=1,n-2 t(i) = sin(x(i+1)-x(n)) - x(1)**2 + x(i+1)**2 f = f + t(i)*t(i) end do return cF63 BIGGSB1 (CUTE) * * Initial Point: [0., 0., ....,0.] 63 continue f=(x(1)-1.d0)**2 + (1.d0-x(n))**2 do i=2,n f = f + (x(i)-x(i-1))**2 end do return cF64 Extended Block-Diagonal BD2 * * Initial Point: [1.5, 2, ..., 1.5, 2] 64 continue f = 0.d0 do i=1,n/2 f = f + (x(2*i-1)**2 + x(2*i)**2 - 2.d0)**2 + * (dexp(x(2*i-1)-1.d0) + x(2*i)**3 - 2.d0)**2 end do return cF65 Generalized quartic GQ1 function * * Initial Point: [1,1,...,1] 65 continue f = 0.d0 do i=1,n-1 f = f + x(i)*x(i) + (x(i+1)+x(i)*x(i))**2 end do return cF66 Diagonal 7 * * Initial point: [1,1,...,1] 66 continue f=0.d0 do i=1,n f = f + dexp(x(i)) - 2.d0*x(i) - x(i)*x(i) end do return cF67 Diagonal 8 * * Initial Point: [1,1,...,1] 67 continue f=0.d0 do i=1,n f = f + x(i)*dexp(x(i)) - 2.d0*x(i) - x(i)*x(i) end do return cF68 Full Hessian * * Initial Point: [1,1,...,1] 68 continue sum=0.d0 do i=1,n sum = sum + x(i) end do f = sum*sum do i=1,n f = f + x(i)*dexp(x(i)) - 2.d0*x(i) - x(i)*x(i) end do return cF69 SINCOS * * Initial Point: [3, 0.1, 3, 0.1, ..., 3, 0.1] 69 continue f= 0.d0 do i=1,n/2 f = f + (x(2*i-1)**2 + x(2*i)**2 + x(2*i-1)*x(2*i))**2 + * (dsin(x(2*i-1)))**2 + (dcos(x(2*i)))**2 end do return cF70 Generalized quartic GQ2 function * * Initial Point: [1,1,...,1] 70 continue f= (x(1)*x(1)-1.d0)**2 do i=2,n f = f + (x(i)*x(i) - x(i-1) - 2.d0)**2 end do return cF71 EXTROSNB (CUTE) * * Initial Point: [1,1,...,1] 71 continue f = (x(1)-1.d0)**2 do i=2,n f = f + 100.d0*(x(i)-x(i-1)**2)**2 end do return cF72 ARGLINB (CUTE) * * Initial Point: [0.01,0.001,...,0.01,0.001] 72 continue m=2 f=0.d0 do i=1,m u(i)=0.d0 do j=1,n u(i) = u(i) + float(i)*float(j)*x(j) end do f = f + (u(i)-1.d0)**2 end do return cF73 FLETCHCR (CUTE) * * Initial Point: [0.5,0.5,...0.5] 73 continue f=0.d0 do i=1,n-1 f = f + 100.d0*(x(i+1)-x(i)+1.d0-x(i)*x(i))**2 end do return cF74 HIMMELBG (CUTE) * * Initial Point: [1.5,1.5,...,1.5] 74 continue f=0.d0 do i=1,n/2 f = f + (2.d0*x(2*i-1)**2+3.d0*x(2*i)**2)* * (dexp(-x(2*i-1)-x(2*i))) end do return cF75 HIMMELBH (CUTE) * * Initial Point: [1.5,1.5,...,1.5] 75 continue f= 0.d0 do i=1,n/2 f=f+(-3.d0*x(2*i-1)-2.d0*x(2*i)+2.d0+x(2*i-1)**3 + x(2*i)**2) end do return cF76 HIMMEL * * Initial Point: [1, 1, ..., 1] 76 continue t1 = (x(1)+x(2)+1.d0)**2 t2 = 19.d0-14.d0*x(1)+3.d0*x(1)*x(1)-14.d0*x(2)+ 1 6.d0*x(1)*x(2)+3.d0*x(2)*x(2) t3 = (2.d0*x(1)-3.d0*x(2))**2 t4 = 18.d0-32.d0*x(1)+12.d0*x(1)*x(1)+48.d0*x(2)- 1 36.d0*x(1)*x(2)+27.d0*x(2)*x(2) f = (1.d0+t1*t2)*(30.d0+t3*t4) return cF77 PROPAN * Identify the concentrations of the products of a hydrocarbon combustion * process at equilibrium. * * Meintjes, K., Morgan, A.P. (1990) Chemical-equilibrium systems as * numerical test problem. ACM Trans. Math. Soft., 16, 143-151. * 77 continue r5=0.193d0 r6=0.4106217541d-3 r7=0.5451766686d-3 r8=0.44975d-6 r9=0.3407354178d-4 r10=0.9615d-6 r=10.d0 fx(1)=x(1)*x(2) + x(1) - 3.d0*x(5) fx(2)=2.d0*x(1)*x(2) + x(1) + 2.d0*r10*x(2)**2 + x(2)*x(3)**2 + 1 r7*x(2)*x(3)+r9*x(2)*x(4)+r8*x(2) - r*x(5) fx(3)=2.d0*x(2)*x(3)**2 + r7*x(2)*x(3) + 2.d0*r5*x(3)**2 + 1 r6*x(3) - 8.d0*x(5) fx(4)=r9*x(2)*x(4) + 2.d0*x(4)**2 - 4.d0*r*x(5) fx(5)=x(1)*x(2) + x(1) + r10*x(2)**2 + x(2)*x(3)**2 + 1 r7*x(2)*x(3) + r9*x(2)*x(4) + r8*x(2) + r5*x(3)**2 + 1 r6*x(3) + x(4)**2 - 1.d0 f=0.d0 do i=1,5 f = f + fx(i)**2 end do return cF78 REACTOR * * Stationar solution of a chemical reactor., * * Shacham, M., (1986) Numerical solution of constrained nonlinear * algebraic equations. * International Journal for numerical methods in engineering, * vol. 123, pp. 1455-1481. * * Please see also: * Neculai Andrei, (2003) Models, Test Problems and Applications for * Mathematical Programming, Technical Press, Bucharest, 2003, * pp. 63 (problem U82) * 78 continue k1d=31.24d0 k2d=0.272d0 k3d=303.03d0 r1=2.062d0 r2=0.02d0 fx(1)=1.d0 - x(1) - k1d*x(1)*x(6) + r1*x(4) fx(2)=1.d0 - x(2) - k2d*x(2)*x(6) + r2*x(5) fx(3)=-x(3) + 2.d0*k3d*x(4)*x(5) fx(4)=k1d*x(1)*x(6) - r1*x(4) - k3d*x(4)*x(5) fx(5)=1.5d0*k2d*x(2)*x(6)-1.5d0*r2*x(5)-k3d*x(4)*x(5) fx(6)=1.d0 - x(4) - x(5) - x(6) f=0.d0 do i=1,6 f = f + fx(i)**2 end do return end c------------------------------------------------ END FUNC.FOR c c Last card FUNC Neculai Andrei