dscpow.for

共轭梯度法源代码

	SUBROUTINE DSCPOW(X,FX,Y,FY,S,DELX,N,INDIC,TOL,ITOL,ITER,IEX)
	DIMENSION X(N),Y(N),S(N),DELX(N)

	IF(ITOL.EQ.1)GO TO 110
	TOL=.001
110	FTOL2=TOL/100
	IEX=1
	K=-2
	M=0
	FA=FX
	DA=0.
	STEP=1.0
	D=STEP
	IF(INDIC.EQ.2)GO TO 1
	IF(ITER.EQ.0)GO TO 1
	DXNORM=0.
	SNORM=0.
	DO 102 I=1,N
	DXNORM=DXNORM+DELX(I)*DELX(I)
102	SNORM=SNORM+S(I)*S(I)
	IF(INDIC.NE.1)GO TO 103
	IF(DXNORM.GE.SNORM)GO TO 1
103	RATIO=DXNORM/SNORM
	STEP=SQRT(RATIO)
	D=STEP
1	DO 2 I=1,N
2	Y(I)=X(I)+D*S(I)
	CALL OBF(Y,F)
	K=K+1
	IF(F-FA)5,3,6
3	D=0.5*(DA+D)
	DO 4 I=1,N
4	Y(I)=X(I)+D*S(I)
	CALL OBF(Y,F)
	IF(F-FA)204,202,205
202	DO 203 I=1,N
203	Y(I)=X(I)+DA*S(I)
	FY=FA
	GO TO 326
204	FC=FA
	DC=DA+2.*(D-DA)
	FB=F
	DB=D
	GO TO 21
205	STEP=0.5*STEP
	IF(K)7,8,206
206	DC=D
	D=DA
	DA=DB
	DB=D
	FC=F
	F=FA
	FA=FB
	FB=F
	GO TO 10
5	FC=FB
	FB=FA
	FA=F
	DC=DB
	DB=DA
	DA=D
	D=2.0*D+STEP
	GO TO 1
6	IF(k)7,8,9
7	FB=F
	DB=D
	D=-D
	STEP=-STEP
	GO TO 1
8	FC=FB
	FB=FA
	FA=F
	DC=DB
	DB=DA
	DA=D
	GO TO 21
9	DC=DB
	DB=DA
	DA=D
	FC=FB
	FB=FA
	FA=F
10	D=0.5*(DA+DB)
	DO 11 I=1,N
11	Y(I)=X(I)+D*S(I)
	CALL OBF(Y,F)
12	IF((DC-D)*(D-DB))15,13,18
13	DO 14 I=1,N
14	Y(I)=X(I)+DB*S(I)
	FY=FB
	IF(IEX.EQ.0)GO TO 32
	GO TO 326
15	IF(F-FB)16,13,17
16	FC=FB
	FB=F
	DC=DB
	DB=D
	GO TO 21
17	FA=F
	DA=D
	GO TO 21
18	IF(F-FB)19,13,20
19	FA=FB
	FB=F
	DA=DB
	DB=D
	GO TO 21
20	FC=F
	DC=D
21	IF(ABS(DA-DB).GE.ABS(DC-DB))GO TO 50
	D=DA
	DA=DC
	DC=D
	F=FA
	FA=FC
	FC=F
50	IF(ABS((DC-DB)/(DA-DB)).GT.0.25)GO TO 26
	D=0.5*(DA+DB)
	DO 51 I=1,N
51	Y(I)=X(I)+D*S(I)
	CALL OBF(Y,F)
	IF(F-FB)52,13,53
52	FC=FB
	FB=F
	DC=DB
	DB=D
	GO TO 26
53	FA=F
	DA=D
26	A=FA*(DB-DC)+FB*(DC-DA)+FC*(DA-DB)
	IF(A)22,30,22
22	D=0.5*((DB*DB-DC*DC)*FA+(DC*DC-DA*DA)*FB+(DA*DA-DB*DB)*FC)/A
	IF((DA-D)*(D-DC))13,13,23
23	DO 24 I=1,N
24	Y(I)=X(I)+D*S(I)
	CALL OBF(Y,F)
	IF(ABS(FB-F)-(ABS(FB)+FTOL2)*TOL)28,28,12
28	IEX=0
	IF(F-FB)29,13,13
29	FY=F
	GO TO 32
30	IF(M)31,31,13
31	M=M+1
	GO TO 10
32	DO 99 I=1,N
	IF(Y(I).NE.X(I))GO TO 326
99	CONTINUE
33	IF(NTOL.EQ.5)GO TO 34
	IEX=1
	NTOL=NTOL+1
	TOL=TOL/10.
	GO TO 12
34	IEX=1
	GO TO 35
326	IF(FY.LT.FX)GO TO 36
	IEX=1
	GO TO 35
36	IF(IEX.EQ.0)GO TO 35
	IEX=2
35	RETURN
	END