您所在位置：主题：请高手看看我的程序。请多多指教
shandong2005
[专家分：0] 发布于 2005-01-10 04:17:00
这是一个混合惩罚函数法的优化计算程序，不知是何原因 ，运算时总是出现除数为零的错误。请高手指教 。用fortran77编的。
─────── YOUHUA.FOR ─────────────────────────────────┐
DIMENSION X(25),GX(50), HX20) ↑│ COMMON /ONE/ ITE,KTE,ILI,NPE,NFX,NGR │ COMMON /TWO/ KCHECK,KDEDI,KWR
COMMON /THR/ EPS1,EPS2,R,T0 │ READ(*,*) N, KG, KH
WRITE(*,1000)
WRITE(*,1005) N, KG, KH
1000 FORMAT(18X,' PRIMARY DATA '/)
1005 FORMAT(5X,'N=',I4,5X,'KG=',I4,5X,'KH=',I4)
│ CALL MAISUB(N,KG,KH,X,GX,HX)
│ STOP
│
│ SUBROUTINE MAISUB(N,KG,KH,X,GX,HX)
│ DIMENSION X(N),X0(25),GX(KG),HX(KH),BL(25),BU(25)
│ COMMON /ONE/ ITE,KTE,ILI,NPE,NFX,NGR
│ COMMON /TWO/ KCHECK,KDEDI,KWR
│ COMMON /THR/ EPS1,EPS2,R,T0
│ READ(*,*) (X(I),I=1,N)
│ READ(*,*) R,C,T0,EPS1,EPS2 ↓ READ(*,*) KCHECK,KDEDI,KWR ↑ READ(*,*) (BL(I),I=1,N),(BU(I),I=1,N)
│ WRITE(*,1100) (X(I),I=1,N)
│ CALL FFX(N,X,FX)
│ CALL GGX(N,KG,X,GX)
│ WRITE(*,1101) FX
│ WRITE(*,1105) (GX(I),I=1,KG)
│ DO 1000 K=1,KG
│ IF(GX(K).GE.0.0) GOTO 1005
1000 CONTINUE
│ GOTO 1010
1005 CALL XRANDO(N,KG,X,GX)
1010 IF(R.LE.0.0) CALL RVALUL(N,KG,X,GX)
│ CALL FUNCTI(N,KG,KH,X,FX,GX,HX,PEN)
1015 WRITE(*,1100) (X(I),I=1,N)
│ WRITE(*,1101) FX
│ WRITE(*,1105) (GX(I),I=1,KG)
│ WRITE(*,1106) (HX(I),I=1,KH)
│ WRITE(*,1075) PEN
│ WRITE(*,1080) R,C,T0,EPS1,EPS2
IRC=0
│ ITE=0
ILI=0
│ NPE=0
│ NFX=0
│ NGR=0
│ DO 1020 I=1,N
1020 X0(I)=X(I)
│ PEN0=PEN
│ WRITE(*,1025)
1025 FORMAT(/18X,' ITERATION COMPUTE '/)
1030 IRC=IRC+1
│ WRITE(*,1035) IRC,R,PEN
1035 FORMAT(/1X,'***** IRC=',I2,4X,' R=',E15.7,4X,'PEN=',E15.7)
│ IF(KCHECK.GT.0) GOTO 1045
│ CALL POWELL(N,KG,KH,X,FX,GX,HX,PEN)
│ ITE=ITE+KTE
│ GOTO 1050
1045 CALL DFPBFG(N,KG,KH,X,FX,GX,HX,PEN)
│ ITE=ITE+KTE
1050 IF(IRC.EQ.1) GOTO 1060
│ IF(ABS(PEN-PEN0).GT.EPS1*ABS(PEN)+EPS2) GOTO 1060
│ DO 1055 I=1,N
│ IF(ABS(X(I)-X0(I)).GT.EPS1*ABS(X(I))+EPS2) GOTO 1060
1055 CONTINUE
│ GOTO 1070
1060 R=R*C
│ DO 1065 I=1,N
1065 X0(I)=X(I)
PEN0=PEN
│ GOTO 1030
1070 WRITE(*,1085)
│ WRITE(*,1090) IRC,ITE,ILI,NPE,NFX,NGR
│ WRITE(*,1095) R,PEN
│ WRITE(*,1100) (X(I),I=1,N)
│ WRITE(*,1101) FX
│ WRITE(*,1105) (GX(I),I=1,KG)
│ WRITE(*,1106) (HX(I),I=1,KH)
1075 FORMAT(5X,'PEN=',E15.7)
1080 FORMAT(5X,'R=',E15.7,5X,'C=',E15.7,' T0=',E15.7/
* 5X,'EPS1=',E15.7,5X,'EPS2=',E15.7/)
1085 FORMAT(/18X,' OPTIMUM SOLUTION '/)
1090 FORMAT(5X,'IRC=',I5,' ITE=',I5,' ILI=',I5,
│ * 'NPE=',I5,'NFX=',I5,'NGR=',I5)
1095 FORMAT(5X,'R=',E15.7,' PEN=',E15.7)
1100 FORMAT(5X,'X :'/(5X,5E15.7))
1101 FORMAT(5X,'FX:'/5X,E15.7)
1105 FORMAT(5X,'GX:'/(5X,5E15.7))
1106 FORMAT(5X,'HX:'/(5X,5E15.7))
│ RETURN
│ END
│
│ SUBROUTINE XRANDO(N,KG,X,GX)
│ DIMENSION X(N),BL(25),BU(25),GX(KG)
│ READ(*,*) (BL(I),I=1,N),(BU(I),I=1,N)
│ RM=2657863.0
100 CALL GGX(N,KG,X,GX)
│ DO 110 K=1,KG
│ IF(GX(K).LT.0.0) GOTO 110
DO 105 I=1,N
│ CALL RANDOM(RM,Q)
105 X(I)=BL(I)+Q*(BU(I)-BL(I))
│ GOTO 100 110 CONTINUE
│ RETURN
│ END
│
│ SUBROUTINE RVALUL(N,KG,X,GX)
│ DIMENSION X(N),GX(KG)
│ COMMON /THR/ EPS1,EPS2,R,T0
│ CALL FFX(N,X,FX)
│ CALL GGX(N,KG,X,GX)
│ PIN=0.0
│ DO 120 I=1,KG
120 PIN=PIN+1.0/ABS(GX(I))
│ R=ABS(FX/PIN)
│ RETURN
│ END
│
│ SUBROUTINE RANDOM(RM,Q)
│ RM35=2.0**35
│ RM36=2.0*RM35
│ RM37=2.0*RM36
│ RM =5.0*RM
│ IF(RM.GE.RM37) RM=RM-RM37
│ IF(RM.GE.RM36) RM=RM-RM36
│ IF(RM.GE.RM35) RM=RM-RM35
│ Q=RM/RM35
│ RETURN
│ END
│ SUBROUTINE POWELL(N,KG,KH,X,FX,GX,HX,F)
│ DIMENSION X(N),X1(25),X2(25),X3(25),S(25)
│ DIMENSION GX(KG),HX(KH),DIRECT(25,26)
│ COMMON /ONE/ ITE,KTE,ILI,NPE,NFX,NGR
│ COMMON /TWO/ KCHECK,KDEDI,KWR
│ COMMON /THR/ EPS1,EPS2,R,T0
│ KTE=0
│ DO 155 I=1,N │ ↓
DO 150 K=1,N
150 DIRECT(K,I)=0.0
155 DIRECT(I,I)=1.0
160 KTE=KTE+1
F1=F
│ F2=F
│ DO 165 I=1,N
│ X1(I)=X(I)
165 X2(I)=X(I)
│ IDF=1
│ DFM=0.0
│ DO 175 I=1,N
│ DO 170 K=1,N
170 S(K)=DIRECT(K,I)
│ CALL LINMIN(N,KG,KH,X2,F2,S,GX,HX)
│ DF=F-F2
│ F=F2
│ IF(DF.LE.DFM) GOTO 175
│ DFM=DF
│ IDF=I
175 CONTINUE
│ DO 180 I=1,N
180 X3(I)=2.0*X2(I)-X1(I)
│ CALL GGX(N,KG,X3,GX)
│ DO 185 K=1,KG
│ IF(GX(K).GT.0.0) GOTO 190 185 CONTINUE
│ CALL FUNCTI(N,KG,KH,X3,FX,GX,HX,F3)
│ AF=(F1-2.0*F2+F3)*(F1-F2-DFM)**2
│ BF=0.5*DFM*(F1-F3)**2
│ IF(F3.LT.FI.AND.AF.LT.BF) GOTO 210
│ IF(F3.LT.F2) GOTO 200 190 DO 195 I=1,N
195 X(I)=X2(I)
│ F=F2
│ GOTO 245
200 DO 205 I=1,N
205 X(I)=X3(I)
│ F=F3
│ GOTO 245
210 N1=N-1
│ IF(IDF.EQ.N) GOTO 220
│ DO 215 I=IDF,N1
│ I1=I+1
│ DO 215 K=1,N
215 DIRECT(K,I)=DIRECT(K,I1) 220 SUM=0.0
│ DO 225 K=1,N
│ DIRECT(K,N)=X3(K)-X1(K)
225 SUM=SUM+DIRECT(K,N)**2
│ SUM=1.0/SQRT(SUM)
│ DO 230 K=1,N
230 DIRECT(K,N)=DIRECT(K,N)*SUM
│ DO 235 K=1,N
235 S(K)=DIRECT(K,N)
│ CALL LINMIN(N,KG,KH,X2,F2,S,GX,HX)
│ DO 240 I=1,N
240 X(I)=X2(I)
│ F=F2
245 CALL FFX(N,X,FX)
CALL GGX(N,KG,X,GX)
│ CALL HHX(N,KH,X,HX)
│ IF(KWR.EQ.0) GOTO 265
│ WRITE(*,250) KTE,F
│ WRITE(*,255) (X(I),I=1,N)
│ WRITE(*,251) FX
│ WRITE(*,260) (GX(I),I=1,KG)
│ WRITE(*,270) (HX(I),I=1,KH)
250 FORMAT(5X,'KTE=',I3,4X,'PEN=',E15.7)
251 FORMAT(5X,'FX:'/5X,E15.7)
255 FORMAT(5X,'X :'/(5X,5E15.7))
260 FORMAT(5X,'GX:'/(5X,5E15.7))
270 FORMAT(5X,'HX:'/(5X,5E15.7))
265 IF(ABS(F-F1).GT.EPS1*ABS(F)+EPS2) GOTO 160
│ RETURN
│ END
│
│ SUBROUTINE DFPBFG(N,KG,KH,X,FX,GX,HX,F)
│ DIMENSION X(N),GX(KG),HX(KH),DX(25),DG(25),HDG(25)
│ DIMENSION DPDX(25),DFDX(25),DGDX(1250),DHDX(500)
DIMENSION S(25),H(25,25)
│ COMMON /ONE/ ITE,KTE,ILI,NPE,NFX,NGR
│ COMMON /TWO/ KCHECK,KDEDI,KWR
│ COMMON /THR/ EPS1,EPS2,R,T0
│ KTE=0
270 DO 280 I=1,N
│ DO 275 K=1,N
275 H(I,K)=0.0
280 H(I,I)=1.0
285 CALL GRAPEN(N,KG,KH,X,FX,GX,HX,DPDX,DFDX,DGDX,DHDX)
│ F0=F
│ GOTO 330 295 DO 300 I=1,N
│ DX(I)=X(I)-DX(I)
300 DG(I)=DPDX(I)-DG(I)
│ IF(KCHECK.EQ.2) GOTO 315
│ CALL ABC1(N,DX,DG,DXTDG)
│ CALL ABC2(N,H,DG,HDG)
│ CALL ABC1(N,DG,HDG,DGTHDG)
│ DO 310 I=1,N
DO 305 K=1,N
│ H(I,K)=H(I,K)+DX(I)*DX(K)/DXTDG-HDG(I)*HDG(K)/DGTHDG
│ IF(I.EQ.K) GOTO 310
305 H(K,I)=H(I,K)
310 CONTINUE
│ GOTO 330 315 CALL ABC2(N,H,DG,HDG)
│ CALL ABC1(N,DG,HDG,DGTHDG)
│ CALL ABC1(N,DX,DG,DXTDG)
│ PSI=1.0+DGTHDG/DXTDG
│ DO 325 I=1,N
│ DO 320 K=1,N
│ H(I,K)=H(I,K)+((DX(I)*PSI-HDG(I))*DX(K)-HDG(K)*DX(I))/DXTDG
│ IF(I.EQ.K) GOTO 325
320 H(K,I)=H(I,K)
325 CONTINUE
330 KTE=KTE+1
F1=F
│ DO 335 I=1,N
│ DX(I)=X(I)
335 DG(I)=DPDX(I)
340 CALL ABC2(N,H,DPDX,S)
│ SUM=0.0
│ DO 345 I=1,N
345 SUM=SUM+S(I)*S(I)
│ SUM=1.0/SQRT(SUM)
│ CALL ABC1(N,DPDX,S,A)
│ IF(A.GT.0.0) GOTO 351
│ A=1.0
│ GOTO 352 351 A=-1.0
352 DO 350 I=1,N
350 S(I)=A*S(I)*SUM
│ CALL LINMIN(N,KG,KH,X,F,S,GX,HX)
│ CALL GRAPEN(N,KG,KH,X,FX,GX,HX,DPDX,DFDX,DGDX,DHDX)
│ IF(KWR.EQ.0) GOTO 372
│ WRITE(*,355) KTE,F
│ WRITE(*,360) (X(I),I=1,N)
│ WRITE(*,365) (DPDX(I),I=1,N)
│ WRITE(*,366) FX
WRITE(*,370) (GX(I),I=1,KG)
│ WRITE(*,371) (HX(I),I=1,KH)
355 FORMAT(5X,'KET='I3,4X,'PEN=',E15.7)
360 FORMAT(5X,'X :'/(5X,5E15.7))
365 FORMAT(5X,'DPDX:'/(5X,E15.7))
366 FORMAT(5X,'FX:'/5X,E15.7)
370 FORMAT(5X,'GX:'/(5X,5E15.7))
371 FORMAT(5X,'HX:'/(5X,5E15.7))
372 IF(F.GT.F0) GOTO 270
│ IF(ABS(F-F1).LE.EPS1*ABS(F)+EPS2) GOTO 375
│ GOTO 295
375 RETURN
END
│
│ SUBROUTINE ABC1(N,A,B,ATB)
│ DIMENSION A(N),B(N)
│ ATB=0.0
│ DO 10 I=1,N
10 ATB=ATB+A(I)*B(I)
│ RETURN
END
│
│ SUBROUTINE ABC2(N,H,A,HA)
│ DIMENSION H(25,25),A(N),HA(N)
│ DO 20 I=1,N
│ HA(I)=0.0
│ DO 20 J=1,N
│20 HA(I)=HA(I)+H(I,J)*A(J)
│ RETURN
│ END │
│ SUBROUTINE LINMIN(N,KG,KH,X,F,S,GX,HX)
│ DIMENSION X(N),XX(25),S(N),GX(KG),HX(KH)
│ COMMON /ONE/ ITE,KTE,ILI,NPE,NFX,NGR
│ COMMON /THR/ EPS1,EPS2,R,T0
│ ILI=ILI+1
│ T1=0.0
│ TH=T0
T2=T0
│ F1=F
400 DO 405 I=1,N
405 XX(I)=X(I)+T2*S(I)
│ CALL GGX(N,KG,XX,GX)
│ DO 410 K=1,KG
│ IF(GX(K).LT.0.0) GOTO 410
│ T2=0.5*T2
│ GOTO 400
410 CONTINUE
│ CALL FUNCTI(N,KG,KH,XX,FX,GX,HX,F2)
│ IF(F2.LE.F1) GOTO 420
│ TH=-TH
│ T3=T1
│ F3=F1
415 T1=T2
│ F1=F2
│ T2=T3
│ F2=F3
420 T3=T2+TH
425 DO 430 I=1,N
430 XX(I)=X(I)+T3*S(I)
│ CALL GGX(N,KG,XX,GX)
│ DO 435 K=1,KG
│ IF(GX(K).LT.0.0) GOTO 435
│ T3=T3-0.5*(T3-T2)
│ IF(ABS(T3-T2).LE.1E-7) GOTO 475
│ GOTO 425
435 CONTINUE
│ CALL FUNCTI(N,KG,KH,XX,FX,GX,HX,F3)
│ IF(F2.LE.F3) GOTO 440
│ TH=TH+TH
│ GOTO 415
440 C1=(F3-F1)/(T3-T1)
│ IF(ABS(T2-T3).LE.1E-7) GOTO 475
│ C2=((F2-F1)/(T2-T1)-C1)/(T2-T3)
│ IF(ABS(C2).LE.1E-7) GOTO 475
│ T4=0.5*(T1+T3-C1/C2)
│ IF((T4-T1)*(T3-T4).LE.0.0) GOTO 475
│ DO 445 I=1,N
445 XX(I)=X(I)+T4*S(I)
│ CALL FUNCTI(N,KG,KH,XX,FX,GX,HX,F4)
│ IF(ABS(F2-F4).LT.EPS1*ABS(F2)+EPS2) GOTO 470
│ IF(ABS(T2-T1).LE.1E-7) GOTO 470
│ IF(ABS(T2-T3).LT.1E-7) GOTO 470
│ IF((T4-T2)*TH.GT.0.0.AND.F2.GT.F4) K=1
│ IF((T4-T2)*TH.GT.0.0.AND.F2.LE.F4) K=2
│ IF((T4-T2)*TH.LE.0.0.AND.F2.GT.F4) K=3
│ IF((T4-T2)*TH.LE.0.0.AND.F2.LE.F4) K=4
│ GOTO (450,455,460,465),K
450 T1=T2
│ F1=F2
│ F2=T4
│ F2=F4
│ GOTO 440
455 T3=T4
│ F3=F4
│ GOTO 440
460 T3=T2
│ F3=F2
T2=T4
│ F2=F4
│ GOTO 440
465 T1=T4
│ F1=F4
│ GOTO 440 470 IF(F4.LT.F2) GOTO 485
475 DO 480 I=1,N
480 X(I)=X(I)+T2*S(I)
│ F=F2
│ RETURN
485 DO 490 I=1,N
490 X(I)=X(I)+T4*S(I)
│ F=F4
│ RETURN
│ END
│
│ SUBROUTINE FUNCTI(N,KG,KH,X,FX,GX,HX,PEN)
│ DIMENSION X(N),GX(KG),HX(KH)
│ COMMON /ONE/ ITE,KTE,ILI,NPE,NFX,NGR
COMMON /THR/ EPS1,EPS2,R,T0
│ NPE=NPE+1
│ PIN=0.0
│ PEX=0.0
│ CALL FFX(N,X,FX)
│ CALL GGX(N,KG,X,GX)
│ CALL HHX(N,KH,X,HX)
│ DO 515 K=1,KG
515 PIN=PIN+1.0/GX(K)
│ DO 510 K=1,KH
510 PEX=PEX+HX(K)*HX(K)
│ PEN=FX-R*PIN+PEX/SQRT(R)
│ RETURN
│ END
│
│ SUBROUTINE GRAPEN(N,KG,KH,X,FX,GX,HX,DPDX,DFDX,DGDX,DHDX)
│ DIMENSION X(N),GX(KG),HX(KH),DPDX(N),DFDX(N)
│ DIMENSION DGDX(N,KG),DHDX(N,KH)
│ COMMON /ONE/ ITE,KTE,ILI,NPE,NFX,NGR
│ COMMON /TWO/ KCHECK,KDEDI,KWR
COMMON /THR/ EPS1,EPS2,R,T0
│ NGR=NGR+1
│ IF(KDEDI) 525,525,530
525 CALL DERIVA(N,KG,KH,X,DFDX,DGDX,DHDX)
│ GOTO 535
530 CALL DIFFER(N,KG,KH,X,FX,GX,HX,DFDX,DGDX,DHDX)
535 DO 560 I=1,N
│ DIN=0.0
│ DEX=0.0
│ DFU=DFDX(I)
│ DO 545 K=1,KG 545 DIN=DIN-DGDX(I,K)/GX(K)**2
│ DO 550 K=1,KH
550 DEX=DEX+2.0*DHDX(I,K)*HX(K)
│ DPDX(I)=DFU+R*DIN+DEX/SQRT(R)
560 CONTINUE
│ RETURN
│ END
│
│ SUBROUTINE DIFFER(N,KG,KH,X,FX,GX,HX,DFDX,DGDX,DHDX)
DIMENSION X(N),GX(KG),HX(KH),DFDX(N),DGDX(N,KG),DHDX(N,KH)
│ COMMON /ONE/ ITE,KTE,ILI,NPE,NFX,NGR
│ COMMON /TWO/ KCHECK,KDEDI,KWR
│ CALL FFX(N,X,GX)
│ CALL GGX(N,KG,X,GX)
│ CALL HHX(N,KH,X,HX)
│ DO 585 I=1,N
│ DJ=1E-6*ABS(X(I))
│ IF(DJ.GT.ABS(FX)) DJ=ABS(FX)
│ DO 570 K=1,KG
│ IF(DJ.GT.ABS(GX(K))) DJ=ABS(GX(K))
570 CONTINUE
│ DO 571 K=1,KH
│ IF(DJ.GT.ABS(HX(K))) DJ=ABS(HX(K))
571 CONTINUE
│ IF(DJ.LE.1E-5) DJ=1E-5
│ XI=X(I)
│ X(I)=XI+DJ
│ CALL FFX(N,X,FX)
│ CALL GGX(N,KG,X,GX)
CALL HHX(N,KH,X,HX)
│ DFDX(I)=FX
│ DO 586 K=1,KG
586 DGDX(I,K)=GX(K)
│ DO 587 K=1,KH
587 DHDX(I,K)=HX(K)
│ X(I)=XI-DJ
│ CALL FFX(N,X,FX)
│ CALL GGX(N,KG,X,GX)
│ CALL HHX(N,KH,X,HX)
│ DFDX(I)=(DFDX(I)-FX)/(2.0*DJ)
│ DO 588 K=1,KG
588 DGDX(I,K)=(DGDX(I,K)-GX(K))/(2.0*DJ)
│ DO 589 K=1,KH
589 DHDX(I,K)=(DHDX(I,K)-HX(K))/(2.0*DJ)
│ X(I)=XI 585 CONTINUE
│ RETURN
│ END
│ SUBROUTINE DERIVA(N,KG,KH,X,DFDX,DGDX,DHDX)
│ DIMENSION X(N),DFDX(N),DGDX(N,KG),DHDX(N,KH)
│ X(1)=X(1)
│ RETURN
│ END
│
│ SUBROUTINE FFX(N,X,FX)
│ DIMENSION X(N)
│ COMMON /ONE/ ITE,KTE,ILI,NPE,NFX,NGR
│ NFX=NFX+1
│ FX=4.0*X(1)-X(2)*X(2)-12.0
│ RETURN
│ END
│ SUBROUTINE GGX(N,KG,X,GX)
│ DIMENSION X(N),GX(KG)
│ GX(1)=34.0-10.0*X(1)-10.0*X(2)+X(1)*X(1)+X(2)*X(2)
│ GX(2)=-X(1)
│ GX(3)=-X(2)
│ RETURN
│ END │ ↓
SUBROUTINE HHX(N,KH,X,HX)
│ DIMENSION X(N),HX(KH)
│ HX(1)=X(1)*X(1)+X(2)*X(2)-25.0
│ RETURN
│ END
─────── YOUHUA.FOR ─────────────────────────────────┐
DIMENSION X(25),GX(50), HX20) ↑│ COMMON /ONE/ ITE,KTE,ILI,NPE,NFX,NGR │ COMMON /TWO/ KCHECK,KDEDI,KWR
COMMON /THR/ EPS1,EPS2,R,T0 │ READ(*,*) N, KG, KH
WRITE(*,1000)
WRITE(*,1005) N, KG, KH
1000 FORMAT(18X,' PRIMARY DATA '/)
1005 FORMAT(5X,'N=',I4,5X,'KG=',I4,5X,'KH=',I4)
│ CALL MAISUB(N,KG,KH,X,GX,HX)
│ STOP
│
│ SUBROUTINE MAISUB(N,KG,KH,X,GX,HX)
│ DIMENSION X(N),X0(25),GX(KG),HX(KH),BL(25),BU(25)
│ COMMON /ONE/ ITE,KTE,ILI,NPE,NFX,NGR
│ COMMON /TWO/ KCHECK,KDEDI,KWR
│ COMMON /THR/ EPS1,EPS2,R,T0
│ READ(*,*) (X(I),I=1,N)
│ READ(*,*) R,C,T0,EPS1,EPS2 ↓ READ(*,*) KCHECK,KDEDI,KWR ↑ READ(*,*) (BL(I),I=1,N),(BU(I),I=1,N)
│ WRITE(*,1100) (X(I),I=1,N)
│ CALL FFX(N,X,FX)
│ CALL GGX(N,KG,X,GX)
│ WRITE(*,1101) FX
│ WRITE(*,1105) (GX(I),I=1,KG)
│ DO 1000 K=1,KG
│ IF(GX(K).GE.0.0) GOTO 1005
1000 CONTINUE
│ GOTO 1010
1005 CALL XRANDO(N,KG,X,GX)
1010 IF(R.LE.0.0) CALL RVALUL(N,KG,X,GX)
│ CALL FUNCTI(N,KG,KH,X,FX,GX,HX,PEN)
1015 WRITE(*,1100) (X(I),I=1,N)
│ WRITE(*,1101) FX
│ WRITE(*,1105) (GX(I),I=1,KG)
│ WRITE(*,1106) (HX(I),I=1,KH)
│ WRITE(*,1075) PEN
│ WRITE(*,1080) R,C,T0,EPS1,EPS2
IRC=0
│ ITE=0
ILI=0
│ NPE=0
│ NFX=0
│ NGR=0
│ DO 1020 I=1,N
1020 X0(I)=X(I)
│ PEN0=PEN
│ WRITE(*,1025)
1025 FORMAT(/18X,' ITERATION COMPUTE '/)
1030 IRC=IRC+1
│ WRITE(*,1035) IRC,R,PEN
1035 FORMAT(/1X,'***** IRC=',I2,4X,' R=',E15.7,4X,'PEN=',E15.7)
│ IF(KCHECK.GT.0) GOTO 1045
│ CALL POWELL(N,KG,KH,X,FX,GX,HX,PEN)
│ ITE=ITE+KTE
│ GOTO 1050
1045 CALL DFPBFG(N,KG,KH,X,FX,GX,HX,PEN)
│ ITE=ITE+KTE
1050 IF(IRC.EQ.1) GOTO 1060
│ IF(ABS(PEN-PEN0).GT.EPS1*ABS(PEN)+EPS2) GOTO 1060
│ DO 1055 I=1,N
│ IF(ABS(X(I)-X0(I)).GT.EPS1*ABS(X(I))+EPS2) GOTO 1060
1055 CONTINUE
│ GOTO 1070
1060 R=R*C
│ DO 1065 I=1,N
1065 X0(I)=X(I)
PEN0=PEN
│ GOTO 1030
1070 WRITE(*,1085)
│ WRITE(*,1090) IRC,ITE,ILI,NPE,NFX,NGR
│ WRITE(*,1095) R,PEN
│ WRITE(*,1100) (X(I),I=1,N)
│ WRITE(*,1101) FX
│ WRITE(*,1105) (GX(I),I=1,KG)
│ WRITE(*,1106) (HX(I),I=1,KH)
1075 FORMAT(5X,'PEN=',E15.7)
1080 FORMAT(5X,'R=',E15.7,5X,'C=',E15.7,' T0=',E15.7/
* 5X,'EPS1=',E15.7,5X,'EPS2=',E15.7/)
1085 FORMAT(/18X,' OPTIMUM SOLUTION '/)
1090 FORMAT(5X,'IRC=',I5,' ITE=',I5,' ILI=',I5,
│ * 'NPE=',I5,'NFX=',I5,'NGR=',I5)
1095 FORMAT(5X,'R=',E15.7,' PEN=',E15.7)
1100 FORMAT(5X,'X :'/(5X,5E15.7))
1101 FORMAT(5X,'FX:'/5X,E15.7)
1105 FORMAT(5X,'GX:'/(5X,5E15.7))
1106 FORMAT(5X,'HX:'/(5X,5E15.7))
│ RETURN
│ END
│
│ SUBROUTINE XRANDO(N,KG,X,GX)
│ DIMENSION X(N),BL(25),BU(25),GX(KG)
│ READ(*,*) (BL(I),I=1,N),(BU(I),I=1,N)
│ RM=2657863.0
100 CALL GGX(N,KG,X,GX)
│ DO 110 K=1,KG
│ IF(GX(K).LT.0.0) GOTO 110
DO 105 I=1,N
│ CALL RANDOM(RM,Q)
105 X(I)=BL(I)+Q*(BU(I)-BL(I))
│ GOTO 100 110 CONTINUE
│ RETURN
│ END
│
│ SUBROUTINE RVALUL(N,KG,X,GX)
│ DIMENSION X(N),GX(KG)
│ COMMON /THR/ EPS1,EPS2,R,T0
│ CALL FFX(N,X,FX)
│ CALL GGX(N,KG,X,GX)
│ PIN=0.0
│ DO 120 I=1,KG
120 PIN=PIN+1.0/ABS(GX(I))
│ R=ABS(FX/PIN)
│ RETURN
│ END
│
│ SUBROUTINE RANDOM(RM,Q)
│ RM35=2.0**35
│ RM36=2.0*RM35
│ RM37=2.0*RM36
│ RM =5.0*RM
│ IF(RM.GE.RM37) RM=RM-RM37
│ IF(RM.GE.RM36) RM=RM-RM36
│ IF(RM.GE.RM35) RM=RM-RM35
│ Q=RM/RM35
│ RETURN
│ END
│ SUBROUTINE POWELL(N,KG,KH,X,FX,GX,HX,F)
│ DIMENSION X(N),X1(25),X2(25),X3(25),S(25)
│ DIMENSION GX(KG),HX(KH),DIRECT(25,26)
│ COMMON /ONE/ ITE,KTE,ILI,NPE,NFX,NGR
│ COMMON /TWO/ KCHECK,KDEDI,KWR
│ COMMON /THR/ EPS1,EPS2,R,T0
│ KTE=0
│ DO 155 I=1,N │ ↓
DO 150 K=1,N
150 DIRECT(K,I)=0.0
155 DIRECT(I,I)=1.0
160 KTE=KTE+1
F1=F
│ F2=F
│ DO 165 I=1,N
│ X1(I)=X(I)
165 X2(I)=X(I)
│ IDF=1
│ DFM=0.0
│ DO 175 I=1,N
│ DO 170 K=1,N
170 S(K)=DIRECT(K,I)
│ CALL LINMIN(N,KG,KH,X2,F2,S,GX,HX)
│ DF=F-F2
│ F=F2
│ IF(DF.LE.DFM) GOTO 175
│ DFM=DF
│ IDF=I
175 CONTINUE
│ DO 180 I=1,N
180 X3(I)=2.0*X2(I)-X1(I)
│ CALL GGX(N,KG,X3,GX)
│ DO 185 K=1,KG
│ IF(GX(K).GT.0.0) GOTO 190 185 CONTINUE
│ CALL FUNCTI(N,KG,KH,X3,FX,GX,HX,F3)
│ AF=(F1-2.0*F2+F3)*(F1-F2-DFM)**2
│ BF=0.5*DFM*(F1-F3)**2
│ IF(F3.LT.FI.AND.AF.LT.BF) GOTO 210
│ IF(F3.LT.F2) GOTO 200 190 DO 195 I=1,N
195 X(I)=X2(I)
│ F=F2
│ GOTO 245
200 DO 205 I=1,N
205 X(I)=X3(I)
│ F=F3
│ GOTO 245
210 N1=N-1
│ IF(IDF.EQ.N) GOTO 220
│ DO 215 I=IDF,N1
│ I1=I+1
│ DO 215 K=1,N
215 DIRECT(K,I)=DIRECT(K,I1) 220 SUM=0.0
│ DO 225 K=1,N
│ DIRECT(K,N)=X3(K)-X1(K)
225 SUM=SUM+DIRECT(K,N)**2
│ SUM=1.0/SQRT(SUM)
│ DO 230 K=1,N
230 DIRECT(K,N)=DIRECT(K,N)*SUM
│ DO 235 K=1,N
235 S(K)=DIRECT(K,N)
│ CALL LINMIN(N,KG,KH,X2,F2,S,GX,HX)
│ DO 240 I=1,N
240 X(I)=X2(I)
│ F=F2
245 CALL FFX(N,X,FX)
CALL GGX(N,KG,X,GX)
│ CALL HHX(N,KH,X,HX)
│ IF(KWR.EQ.0) GOTO 265
│ WRITE(*,250) KTE,F
│ WRITE(*,255) (X(I),I=1,N)
│ WRITE(*,251) FX
│ WRITE(*,260) (GX(I),I=1,KG)
│ WRITE(*,270) (HX(I),I=1,KH)
250 FORMAT(5X,'KTE=',I3,4X,'PEN=',E15.7)
251 FORMAT(5X,'FX:'/5X,E15.7)
255 FORMAT(5X,'X :'/(5X,5E15.7))
260 FORMAT(5X,'GX:'/(5X,5E15.7))
270 FORMAT(5X,'HX:'/(5X,5E15.7))
265 IF(ABS(F-F1).GT.EPS1*ABS(F)+EPS2) GOTO 160
│ RETURN
│ END
│
│ SUBROUTINE DFPBFG(N,KG,KH,X,FX,GX,HX,F)
│ DIMENSION X(N),GX(KG),HX(KH),DX(25),DG(25),HDG(25)
│ DIMENSION DPDX(25),DFDX(25),DGDX(1250),DHDX(500)
DIMENSION S(25),H(25,25)
│ COMMON /ONE/ ITE,KTE,ILI,NPE,NFX,NGR
│ COMMON /TWO/ KCHECK,KDEDI,KWR
│ COMMON /THR/ EPS1,EPS2,R,T0
│ KTE=0
270 DO 280 I=1,N
│ DO 275 K=1,N
275 H(I,K)=0.0
280 H(I,I)=1.0
285 CALL GRAPEN(N,KG,KH,X,FX,GX,HX,DPDX,DFDX,DGDX,DHDX)
│ F0=F
│ GOTO 330 295 DO 300 I=1,N
│ DX(I)=X(I)-DX(I)
300 DG(I)=DPDX(I)-DG(I)
│ IF(KCHECK.EQ.2) GOTO 315
│ CALL ABC1(N,DX,DG,DXTDG)
│ CALL ABC2(N,H,DG,HDG)
│ CALL ABC1(N,DG,HDG,DGTHDG)
│ DO 310 I=1,N
DO 305 K=1,N
│ H(I,K)=H(I,K)+DX(I)*DX(K)/DXTDG-HDG(I)*HDG(K)/DGTHDG
│ IF(I.EQ.K) GOTO 310
305 H(K,I)=H(I,K)
310 CONTINUE
│ GOTO 330 315 CALL ABC2(N,H,DG,HDG)
│ CALL ABC1(N,DG,HDG,DGTHDG)
│ CALL ABC1(N,DX,DG,DXTDG)
│ PSI=1.0+DGTHDG/DXTDG
│ DO 325 I=1,N
│ DO 320 K=1,N
│ H(I,K)=H(I,K)+((DX(I)*PSI-HDG(I))*DX(K)-HDG(K)*DX(I))/DXTDG
│ IF(I.EQ.K) GOTO 325
320 H(K,I)=H(I,K)
325 CONTINUE
330 KTE=KTE+1
F1=F
│ DO 335 I=1,N
│ DX(I)=X(I)
335 DG(I)=DPDX(I)
340 CALL ABC2(N,H,DPDX,S)
│ SUM=0.0
│ DO 345 I=1,N
345 SUM=SUM+S(I)*S(I)
│ SUM=1.0/SQRT(SUM)
│ CALL ABC1(N,DPDX,S,A)
│ IF(A.GT.0.0) GOTO 351
│ A=1.0
│ GOTO 352 351 A=-1.0
352 DO 350 I=1,N
350 S(I)=A*S(I)*SUM
│ CALL LINMIN(N,KG,KH,X,F,S,GX,HX)
│ CALL GRAPEN(N,KG,KH,X,FX,GX,HX,DPDX,DFDX,DGDX,DHDX)
│ IF(KWR.EQ.0) GOTO 372
│ WRITE(*,355) KTE,F
│ WRITE(*,360) (X(I),I=1,N)
│ WRITE(*,365) (DPDX(I),I=1,N)
│ WRITE(*,366) FX
WRITE(*,370) (GX(I),I=1,KG)
│ WRITE(*,371) (HX(I),I=1,KH)
355 FORMAT(5X,'KET='I3,4X,'PEN=',E15.7)
360 FORMAT(5X,'X :'/(5X,5E15.7))
365 FORMAT(5X,'DPDX:'/(5X,E15.7))
366 FORMAT(5X,'FX:'/5X,E15.7)
370 FORMAT(5X,'GX:'/(5X,5E15.7))
371 FORMAT(5X,'HX:'/(5X,5E15.7))
372 IF(F.GT.F0) GOTO 270
│ IF(ABS(F-F1).LE.EPS1*ABS(F)+EPS2) GOTO 375
│ GOTO 295
375 RETURN
END
│
│ SUBROUTINE ABC1(N,A,B,ATB)
│ DIMENSION A(N),B(N)
│ ATB=0.0
│ DO 10 I=1,N
10 ATB=ATB+A(I)*B(I)
│ RETURN
END
│
│ SUBROUTINE ABC2(N,H,A,HA)
│ DIMENSION H(25,25),A(N),HA(N)
│ DO 20 I=1,N
│ HA(I)=0.0
│ DO 20 J=1,N
│20 HA(I)=HA(I)+H(I,J)*A(J)
│ RETURN
│ END │
│ SUBROUTINE LINMIN(N,KG,KH,X,F,S,GX,HX)
│ DIMENSION X(N),XX(25),S(N),GX(KG),HX(KH)
│ COMMON /ONE/ ITE,KTE,ILI,NPE,NFX,NGR
│ COMMON /THR/ EPS1,EPS2,R,T0
│ ILI=ILI+1
│ T1=0.0
│ TH=T0
T2=T0
│ F1=F
400 DO 405 I=1,N
405 XX(I)=X(I)+T2*S(I)
│ CALL GGX(N,KG,XX,GX)
│ DO 410 K=1,KG
│ IF(GX(K).LT.0.0) GOTO 410
│ T2=0.5*T2
│ GOTO 400
410 CONTINUE
│ CALL FUNCTI(N,KG,KH,XX,FX,GX,HX,F2)
│ IF(F2.LE.F1) GOTO 420
│ TH=-TH
│ T3=T1
│ F3=F1
415 T1=T2
│ F1=F2
│ T2=T3
│ F2=F3
420 T3=T2+TH
425 DO 430 I=1,N
430 XX(I)=X(I)+T3*S(I)
│ CALL GGX(N,KG,XX,GX)
│ DO 435 K=1,KG
│ IF(GX(K).LT.0.0) GOTO 435
│ T3=T3-0.5*(T3-T2)
│ IF(ABS(T3-T2).LE.1E-7) GOTO 475
│ GOTO 425
435 CONTINUE
│ CALL FUNCTI(N,KG,KH,XX,FX,GX,HX,F3)
│ IF(F2.LE.F3) GOTO 440
│ TH=TH+TH
│ GOTO 415
440 C1=(F3-F1)/(T3-T1)
│ IF(ABS(T2-T3).LE.1E-7) GOTO 475
│ C2=((F2-F1)/(T2-T1)-C1)/(T2-T3)
│ IF(ABS(C2).LE.1E-7) GOTO 475
│ T4=0.5*(T1+T3-C1/C2)
│ IF((T4-T1)*(T3-T4).LE.0.0) GOTO 475
│ DO 445 I=1,N
445 XX(I)=X(I)+T4*S(I)
│ CALL FUNCTI(N,KG,KH,XX,FX,GX,HX,F4)
│ IF(ABS(F2-F4).LT.EPS1*ABS(F2)+EPS2) GOTO 470
│ IF(ABS(T2-T1).LE.1E-7) GOTO 470
│ IF(ABS(T2-T3).LT.1E-7) GOTO 470
│ IF((T4-T2)*TH.GT.0.0.AND.F2.GT.F4) K=1
│ IF((T4-T2)*TH.GT.0.0.AND.F2.LE.F4) K=2
│ IF((T4-T2)*TH.LE.0.0.AND.F2.GT.F4) K=3
│ IF((T4-T2)*TH.LE.0.0.AND.F2.LE.F4) K=4
│ GOTO (450,455,460,465),K
450 T1=T2
│ F1=F2
│ F2=T4
│ F2=F4
│ GOTO 440
455 T3=T4
│ F3=F4
│ GOTO 440
460 T3=T2
│ F3=F2
T2=T4
│ F2=F4
│ GOTO 440
465 T1=T4
│ F1=F4
│ GOTO 440 470 IF(F4.LT.F2) GOTO 485
475 DO 480 I=1,N
480 X(I)=X(I)+T2*S(I)
│ F=F2
│ RETURN
485 DO 490 I=1,N
490 X(I)=X(I)+T4*S(I)
│ F=F4
│ RETURN
│ END
│
│ SUBROUTINE FUNCTI(N,KG,KH,X,FX,GX,HX,PEN)
│ DIMENSION X(N),GX(KG),HX(KH)
│ COMMON /ONE/ ITE,KTE,ILI,NPE,NFX,NGR
COMMON /THR/ EPS1,EPS2,R,T0
│ NPE=NPE+1
│ PIN=0.0
│ PEX=0.0
│ CALL FFX(N,X,FX)
│ CALL GGX(N,KG,X,GX)
│ CALL HHX(N,KH,X,HX)
│ DO 515 K=1,KG
515 PIN=PIN+1.0/GX(K)
│ DO 510 K=1,KH
510 PEX=PEX+HX(K)*HX(K)
│ PEN=FX-R*PIN+PEX/SQRT(R)
│ RETURN
│ END
│
│ SUBROUTINE GRAPEN(N,KG,KH,X,FX,GX,HX,DPDX,DFDX,DGDX,DHDX)
│ DIMENSION X(N),GX(KG),HX(KH),DPDX(N),DFDX(N)
│ DIMENSION DGDX(N,KG),DHDX(N,KH)
│ COMMON /ONE/ ITE,KTE,ILI,NPE,NFX,NGR
│ COMMON /TWO/ KCHECK,KDEDI,KWR
COMMON /THR/ EPS1,EPS2,R,T0
│ NGR=NGR+1
│ IF(KDEDI) 525,525,530
525 CALL DERIVA(N,KG,KH,X,DFDX,DGDX,DHDX)
│ GOTO 535
530 CALL DIFFER(N,KG,KH,X,FX,GX,HX,DFDX,DGDX,DHDX)
535 DO 560 I=1,N
│ DIN=0.0
│ DEX=0.0
│ DFU=DFDX(I)
│ DO 545 K=1,KG 545 DIN=DIN-DGDX(I,K)/GX(K)**2
│ DO 550 K=1,KH
550 DEX=DEX+2.0*DHDX(I,K)*HX(K)
│ DPDX(I)=DFU+R*DIN+DEX/SQRT(R)
560 CONTINUE
│ RETURN
│ END
│
│ SUBROUTINE DIFFER(N,KG,KH,X,FX,GX,HX,DFDX,DGDX,DHDX)
DIMENSION X(N),GX(KG),HX(KH),DFDX(N),DGDX(N,KG),DHDX(N,KH)
│ COMMON /ONE/ ITE,KTE,ILI,NPE,NFX,NGR
│ COMMON /TWO/ KCHECK,KDEDI,KWR
│ CALL FFX(N,X,GX)
│ CALL GGX(N,KG,X,GX)
│ CALL HHX(N,KH,X,HX)
│ DO 585 I=1,N
│ DJ=1E-6*ABS(X(I))
│ IF(DJ.GT.ABS(FX)) DJ=ABS(FX)
│ DO 570 K=1,KG
│ IF(DJ.GT.ABS(GX(K))) DJ=ABS(GX(K))
570 CONTINUE
│ DO 571 K=1,KH
│ IF(DJ.GT.ABS(HX(K))) DJ=ABS(HX(K))
571 CONTINUE
│ IF(DJ.LE.1E-5) DJ=1E-5
│ XI=X(I)
│ X(I)=XI+DJ
│ CALL FFX(N,X,FX)
│ CALL GGX(N,KG,X,GX)
CALL HHX(N,KH,X,HX)
│ DFDX(I)=FX
│ DO 586 K=1,KG
586 DGDX(I,K)=GX(K)
│ DO 587 K=1,KH
587 DHDX(I,K)=HX(K)
│ X(I)=XI-DJ
│ CALL FFX(N,X,FX)
│ CALL GGX(N,KG,X,GX)
│ CALL HHX(N,KH,X,HX)
│ DFDX(I)=(DFDX(I)-FX)/(2.0*DJ)
│ DO 588 K=1,KG
588 DGDX(I,K)=(DGDX(I,K)-GX(K))/(2.0*DJ)
│ DO 589 K=1,KH
589 DHDX(I,K)=(DHDX(I,K)-HX(K))/(2.0*DJ)
│ X(I)=XI 585 CONTINUE
│ RETURN
│ END
│ SUBROUTINE DERIVA(N,KG,KH,X,DFDX,DGDX,DHDX)
│ DIMENSION X(N),DFDX(N),DGDX(N,KG),DHDX(N,KH)
│ X(1)=X(1)
│ RETURN
│ END
│
│ SUBROUTINE FFX(N,X,FX)
│ DIMENSION X(N)
│ COMMON /ONE/ ITE,KTE,ILI,NPE,NFX,NGR
│ NFX=NFX+1
│ FX=4.0*X(1)-X(2)*X(2)-12.0
│ RETURN
│ END
│ SUBROUTINE GGX(N,KG,X,GX)
│ DIMENSION X(N),GX(KG)
│ GX(1)=34.0-10.0*X(1)-10.0*X(2)+X(1)*X(1)+X(2)*X(2)
│ GX(2)=-X(1)
│ GX(3)=-X(2)
│ RETURN
│ END │ ↓
SUBROUTINE HHX(N,KH,X,HX)
│ DIMENSION X(N),HX(KH)
│ HX(1)=X(1)*X(1)+X(2)*X(2)-25.0
│ RETURN
│ END