sumt.for

来自「内点惩罚法,fortran77编制可以在visual fortran运行」· FOR 代码 · 共 300 行

	PROGRAM SUMT
	DIMENSION X(60),X0(60),XX(60),X3(60),BL(60),BU(60),
     1  GX(150),HX(30),SS(60),S(60,61)
	WRITE(*,123)
 123	FORMAT(5X,1HN,5X,2HKG,5X,2HKH)
	READ(*,*)N,KG,KH
	WRITE(*,1002) 
 	WRITE(*,1003) N,KG,KH
1002    FORMAT(25X,'-----PRIMARY DATA-----',/)
1003 	FORMAT(5X,2HN=,I4,5X,3HKG=,I4,5X,3HKH=,I4)
	N1=N+1
	CALL MAISUB(N,N1,KG,KH,X,X0,XX,X3,BL,BU,FX,GX,HX,SS,S)
	STOP
	END
C
	SUBROUTINE MAISUB(N,N1,KG,KH,X,X0,XX,X3,BL,BU,FX,GX,HX,SS,S)
	DIMENSION X(N),X0(N),XX(N),X3(N),BL(N),BU(N),GX(KG),
     1  HX(KH),SS(N),S(N,N1)
	WRITE(*,1113) 
 1113	FORMAT(5X,2HX:)
	READ(*,*) (X(I),I=1,N)
        WRITE(*,1123) (X(I),I=1,N)                                      
        WRITE(*,1116)
1116	FORMAT(5X,3HBL:)
	READ(*,*) (BL(I),I=1,N)
        WRITE(*,1123) (BL(I),I=1,N)                                     
        WRITE(*,1115)
 1115	FORMAT(5X,3HBU:)
	READ(*,*)(BU(I),I=1,N)
        WRITE(*,1123) (BU(I),I=1,N)                                     
1123    FORMAT(10X,4(F10.4,8X))                                          
        WRITE(*,1114)
 1114   FORMAT(5X,'EP:',11X,'EPS:',13X,'C:',8X,'H0:',8X,'R:')           
	READ(*,*) EP,EPS,C,H0,R
        WRITE(*,1124) EP,EPS,C,H0,R                                     
 1124   FORMAT(6X,2(F10.7,5X),3(F6.3,5X))                               
2000	CALL XRAN(N,X,BL,BU,GX,KG)
	IXF=0
	IF(R.LE.0.0) CALL RVALUE(R,N,KG,KH,FX,GX,HX,X,IXF)
	DO 1006 I=1,N
	DO 1005 K=1,N1
1005	S(I,K)=0.
1006	S(I,I)=1.
	IRC=0
	IQU=0
	DO 1007 K=1,N
1007	X0(K)=X(K)
	CALL FUNCT(N,KG,KH,X,FX,GX,HX,F0,R,IXF)
	FOM=1.E15
1013	IRC=IRC+1
	CALL POWELL(N,N1,X,X0,X3,S,SS,EPS,F0,IQU,IXF,
     1  FX,GX,HX,KG,KH,XX,R,H0)
	IF(ABS((FOM-F0)/F0).LT.EP) GOTO 1016
	FOM=F0
	R=R*C
	GOTO 1013
1016	WRITE(*,1017)
1017	FORMAT(2X,/25X,'-----OPTIMUM SOLUTION-----',/)
	WRITE(*,1018)IRC,IQU,IXF
	WRITE(*,1019)R,F0
	WRITE(*,1020) (X(I),I=1,N)
1018	FORMAT(5X,4HIRC=,I5,5X,4HIQU=,I5,5X,4HIXF=,I5)
1019	FORMAT(5X,3HR =,E15.6/5X,3HF0=,E15.6)
1020	FORMAT(5X,3HX :/(5X,5E15.6))
	CALL FFX(N,X,FX,IXF)
	CALL GGX(N,KG,X,GX)
	CALL HHX(N,KH,X,HX)
	WRITE(*,1021) FX
	WRITE(*,1022)(GX(I),I=1,KG)
 	WRITE(*,1023)(HX(I),I=1,KH)
1021	FORMAT(5X,3HFX=,E15.6)
1022 	FORMAT(5X,3HGX:/(5X,5E15.6))
1023	FORMAT(5X,3HHX:/(5X,5E15.6))
	RETURN
	END
C
	SUBROUTINE XRAN(N,X,BL,BU,GX,KG)
	DIMENSION X(N),BL(N),BU(N),GX(KG)
	RM=2657863.
401     CALL GGX(N,KG,X,GX)
	DO 403 I=1,KG
	IF(GX(I).GT.0.0) GOTO 403
	DO 402 K=1,N
	CALL RANDOM(RM,Q)
402     X(K)=BL(K)+Q*(BU(K)-BL(K))
	GOTO 401
403     CONTINUE
	RETURN
	END
C
	SUBROUTINE RVALUE(R,N,KG,KH,FX,GX,HX,X,IXF)
	DIMENSION X(N),GX(KG),HX(KH)
	CALL FFX(N,X,FX,IXF)
	CALL GGX(N,KG,X,GX)
	SG=0.
	DO 501 I=1,KG
501     SG=SG+1./GX(I)
	CALL HHX(N,KH,X,HX)
	SH=0.
	DO 502 I=1,KH
502     SH=SH+HX(I)*HX(I)
	IF(ABS(SH).LE.1.E-15) GOTO 503
	R=ABS(FX/SG/SH)
	RETURN
503     R=ABS(FX/SG)
	RETURN
	END
C
	SUBROUTINE RANDOM(RM,Q)
	RM35=2.**35
	RM36=2.*RM35
	RM37=2.*RM36
	RM =5.*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
C
	SUBROUTINE QUAINS(HO,YO,Y,N,X,XX,SS,EPS,IQU,
     1  IXF,FX,GX,HX,KG,KH,R)
	DIMENSION X(N),XX(N),SS(N),GX(KG),HX(KH)
	IQU=IQU+1
	HT=HO
	T2=HO
	T1=0.
	Y0=F0
	Y1=YO
100     DO 101 K=1,N
101     XX(K)=X(K)+T2*SS(K)
	CALL GGX(N,KG,XX,GX)
	DO 102 K=1,KG
	IF(GX(K).GE.1.E-15) GOTO 102
	T2=0.5*T2
	GOTO 100
102     CONTINUE
	CALL FUNCT(N,KG,KH,XX,FX,GX,HX,Y2,R,IXF)
	IF(Y2.LE.Y1) GOTO 104
	HT=-HT
	T3=T1
	Y3=Y1
103     T1=T2
	Y1=Y2
	T2=T3
	Y2=Y3
104     T3=T2+HT
105     DO 106 K=1,N
106     XX(K)=X(K)+T3*SS(K)
	CALL GGX(N,KG,XX,GX)
	DO 107 K=1,KG
	IF(GX(K).GE.1E-15) GOTO 107
	T3=T3-0.5*(T3-T2)
	GOTO 105
107     CONTINUE
	CALL FUNCT(N,KG,KH,XX,FX,GX,HX,Y3,R,IXF)
	IF(Y2.GT.Y3) GOTO 116
108	C1=(Y3-Y1)/(T3-T1)
	C2=((Y2-Y1)/(T2-T1)-C1)/(T2-T3)
	IF(ABS(C2).LE.1.E-15) GOTO 118
	T4=0.5*(T1+T3-C1/C2)
	IF((T4-T1)*(T3-T4).LE.0.) GOTO 118
	DO 109 K=1,N
109	XX(K)=X(K)+T4*SS(K)
	CALL FUNCT(N,KG,KH,XX,FX,GX,HX,Y4,R,IXF)
	IF(ABS(Y2).LT.1.) GOTO 110
	A=Y2
	GOTO 111
110	A=1.
111	IF(ABS((Y2-Y4)/A).LT.EPS) GOTO 117
	IF(ABS(T2-T1).LT.1E-5) GOTO 117
	IF(ABS(T2-T3).LT.1E-5) GOTO 117
	IF((T4-T2)*HT.GT.0.0.AND.Y2.GT.Y4) K=1
	IF((T4-T2)*HT.GT.0.0.AND.Y2.LE.Y4) K=2
	IF((T4-T2)*HT.LE.0.0.AND.Y2.GT.Y4) K=3
	IF((T4-T2)*HT.LE.0.0.AND.Y2.LE.Y4) K=4
	GOTO (112,113,114,115),K
112	T1=T2
	T2=T4
	Y1=Y2
	Y2=Y4
	GOTO 108
113 	T3=T4
	Y3=Y4
	GOTO 108
114	T3=T2
	Y3=Y2
	T2=T4
	Y2=Y4
	GOTO 108
115	T1=T4
	Y1=Y4
	GOTO 108
116	HT=HT+HT
	GOTO 103
117	IF(Y4.LT.Y2) GOTO 120
118 	DO 119 K=1,N
119	X(K)=X(K)+T2*SS(K)
	Y=Y2
	RETURN
120	DO 121 K=1,N
121	X(K)=X(K)+T4*SS(K)
	Y=Y4
	RETURN
	END
C
	SUBROUTINE POWELL(N,N1,X,X0,X3,S,SS,EPS,F0,IQU,
     1  IXF,FX,GX,HX,KG,KH,XX,R,H0)
	DIMENSION X(N),X0(N),X3(N),S(N,N1),SS(N),GX(KG),HX(KH),XX(N)
	ITE=0
201	ITE=ITE+1
	F1=F0
	IDF=1
	DFM=0.0
	DO 203 I=1,N
	DO 202 K=1,N
202	SS(K)=S(K,I)
	CALL QUAINS(H0,F1,F2,N,X,XX,SS,EPS,IQU,IXF,FX,GX,HX,
     1  KG,KH,R)
	DF=F1-F2
	F1=F2
	IF(DF.LE.DFM) GOTO 203
	DFM=DF
	IDF=I
203	CONTINUE
	SDX=0.0
	DO 204 K=1,N
204 	SDX=SDX+(X(K)-X0(K))**2
	IF(SDX.LE.EPS) GOTO 214
	DO 205 K=1,N
	X3(K)=2.0*X(K)-X0(K)
205	S(K,N1)=X(K)-X0(K)
	CALL GGX(N,KG,X3,GX)
	DO 2051 K=1,KG
	IF(GX(K).LE.1E-15) GOTO 2052
2051 	CONTINUE
	CALL FUNCT(N,KG,KH,X3,FX,GX,HX,F3,R,IXF)
	AF=(F0-2.*F2+F3)*(F0-F2-DFM)**2
	BF=0.5*DFM*(F0-F3)**2
	IF(F3.LT.F0.AND.AF.LT.BF) GOTO 209
	IF(F3.LT.F2) GOTO 207
2052	DO 206 K=1,N
206	X0(K)=X(K)
	F0=F2
 	GOTO 213
207	F0=F3
	DO 208 K=1,N
208     X0(K)=X3(K)
	X(K)=X3(K)
	GOTO 213
209 	DO 210 I=IDF,N
	DO 210 K=1,N
	M=I+1
210	S(K,I)=S(K,M)
	DO 211 K=1,N
211	SS(K)=S(K,N)
	CALL QUAINS(H0,F2,F0,N,X,XX,SS,EPS,IQU,IXF,FX,GX,HX,
     1  KG,KH,R)
	DO 212 K=1,N
212	X0(K)=X(K)
	CALL FUNCT(N,KG,KH,X,FX,GX,HX,F0,R,IXF)
213	GOTO 201
214	RETURN
	END
C
	SUBROUTINE  FUNCT(N,KG,KH,X,FX,GX,HX,Y,R,IXF)
	DIMENSION GX(KG),HX(KH),X(N)
	SG=0.0
	SH=0.
	CALL FFX(N,X,FX,IXF)
	CALL GGX(N,KG,X,GX)
	CALL HHX(N,KH,X,HX)
	DO 701 I=1,KG
701	SG=SG+1./GX(I)
	DO 702 I=1,KH
702 	SH=SH+HX(I)*HX(I)
	Y=FX+R*SG+SH/SQRT(R)
	RETURN
	END

	SUBROUTINE FFX(N,X,FX,IXF)
	DIMENSION X(N)
	IXF=IXF+1
	FX=120.*X(1)+X(2)
	RETURN
	END
C
	SUBROUTINE GGX(N,KG,X,GX)
	DIMENSION X(N),GX(KG)
	GX(1)=X(2)-4.
	GX(2)=X(1)*X(2)-6.429
	GX(3)=X(2)*X(1)**3-6.429
	GX(4)=X(2)**2*X(1)-321.14
	RETURN
	END
C
	SUBROUTINE HHX(N,KH,X,HX)
	RETURN
	END