FORTRAN版阻尼最小二乘法

thinco

最优化程序汇编  常用程序集
工人出版社

C**********************************************************************
 &nbspROGRAM TAY
C**********************************************************************
C  DIMENSION X(N),H(N),F(M),DF(M,N),AA(N,N)
  DIMENSION X(2),H(2),F(479),DF(479,2),AA(2,2),TT(479),U(479)
  COMMON NF
  COMMON TT,U
  OPEN(2,FILE='CL.DAT')
  OPEN(3,FILE='RESULT2.DAT')
  DO 5 I=1,479
5  READ(2,*,END=6) TT( I ),U( I )
6  N=2
  M=479
  DO 10 I=1,N
  H( I )=0.001
10  CONTINUE
  NF=0
  ITMAX=1000
  S=1.0E20
  X(1)=60
  X(2)=4
  CALL TAYLOR(N,M,X,H,F,ITMAX,0.5E-1,0.5E-1,S,KENN,DF,AA)
  IF(KENN.EQ.1) WRITE(*,*)'NO REACH CONVERGENCE PRECISION WHEN
     *                         ITMAX=1000'
  IF(KENN.EQ.-1) WRITE(*,*)'MINIZATION FAILED'  
  IF(KENN.EQ.0) WRITE(*,30)
30  FORMAT(/1X, 'REACH CONVERGENCE PRECISION '/)
  WRITE(*,*) 'X( I ),I=1,N'
  WRITE(*,40) (I,X( I ),I=1,N)
  WRITE(*,*) 'F( I ),I=1,N'
  WRITE(*,40) (I,F( I ),I=1,N)
  WRITE(3,*) 'X( I ),I=1,N'
  WRITE(3,40) (I,X( I ),I=1,N)
  WRITE(3,*) 'F( I ) ,I=1,N'
  WRITE(3,40) (I,F( I ),I=1,N)
40  FORMAT(1X,5(3X,I3,2X,E15.8))
  WRITE(*,50) ITMAX,NF,S
  WRITE(6,50) ITMAX,NF,S
50  FORMAT(1X,'ITMAX=',I3,3X,'NF=',I4,3X,'S=',E15.8)  
  STOP
  END

C**********************************************************************
  SUBROUTINE TAYLOR(N,M,X,H,F,ITMAX,EPS1,EPS2,S,KENN,DFDX,AA)
C**********************************************************************
  DIMENSION FP(500),FM(500),B(30),DX(30),DFDX(M,N),AA(N,N),
     * X(N),F(M),H(N)
  HS=S
  KENN=0
  IZ=0
23  IZ=IZ+1
  IF(IZ-ITMAX) 1,1,2
2  KENN=1
  GO TO 25
1  L=0
  HL=1.0
24  L=L+1
  IF(L-16) 3,3,4
4  KENN=-1
  GO TO 25
3  CALL FUNCT(X,F,N,M)
  HF=0
  DO 5 I=1,M
5  HF=HF+F( I )*F( I )
  IF(HF-HS*(1.0-0.2*HL)) 6,6,7
7  HL=0.5*HL
  DO 8 K=1,N
8  X(K)=X(K)+HL*DX(K)
  GO TO 24
6  HS=HF
  IF(HS.LT.EPS1) GO TO 25
  DO 12 I=1,N
  HF=H( I )
  HZ=2.0*HF
  HH=X( I )
  X( I )=X( I )+HF
  CALL FUNCT(X,FP,N,M)
  X( I )=X( I )-HZ
  CALL FUNCT(X,FM,N,M)
  X( I )=HH
  HZ=1.0/HZ
  DO 12 K=1,M
12  DFDX(K,I)=HZ*(FP(K)-FM(K))
  IF(M-N) 14,15,14
15  CALL GAUSS(N,DFDX,F,DX,1.0E-10,N)
  GO TO 20
14  DO 17 I=1,N
  HF=0.0
  DO 18 K=1,M
18  HF=HF+DFDX(K,I)*F(K)
  B( I )=HF
  DO 17 K=1,N
  HF=0.0
  DO 16 J=1,M
16  HF=HF+DFDX(J,I)*DFDX(J,K)
  AA(I,K)=HF
17  AA(K,I)=HF
  CALL GAUSS(N,AA,B,DX,1.E-10,N)
20  HZ=0.0
  HF=0.0
  DO 21 I=1,N
  X( I )=X( I )-DX( I )
  HZ=HZ+ABS(X( I ))
21  HF=HF+ABS(DX( I ))
  IF(HF.GE.EPS2*HZ) GO TO 23
25  CALL FUNCT(X,F,N,M)
  S=0.0
  ITMAX=IZ
  DO 22 I=1,M
22  S=S+F( I )*F( I )
  RETURN
  END
C**********************************************************************
  SUBROUTINE GAUSS(N,A,B,X,EP,M)
C**********************************************************************
  DIMENSION A(M,M),B(N),X(N)
  DO 10 K=1,N
  C=0.0
  DO 9 I=K,N
  IF(ABS(A(I,K)).LE.ABS(C)) GO TO 9
  C=A(I,K)
  I0=I
9  CONTINUE
  IF(ABS(C).LE.EP) PAUSE 1
  IF(I0.EQ.K) GO TO 13
  DO 12 J=K,N
  T=A(K,J)
  A(K,J)=A(I0,J)
12  A(I0,J)=T
  T=B(K)
  B(K)=B(I0)
  B(I0)=T
13  C=1.0/C
  IK=K+1
  DO 15 J=IK,N
  A(K,J)=A(K,J)*C
  DO 16 I=IK,N
16  A(I,J)=A(I,J)-A(I,K)*A(K,J)
15  CONTINUE
  B(K)=B(K)*C
  DO 8 I=IK,N
8  B( I )=B( I )-B(K)*A(I,K)
10  CONTINUE
  DO 18 I=N,1,-1
  DO 18 J=I+1,N
  B( I )=B( I )-A(I,J)*B(J)
18  CONTINUE
  DO 11 I=1,N
11  X( I )=B( I )
  RETURN
  END
C**********************************************************************
  SUBROUTINE FUNCT(X,F,N,M)
C**********************************************************************
  DIMENSION X(N),F(M),TT(479),U(479),FF(479)
  COMMON NF
  COMMON TT,U
  P0=316.77
  WRITE(*,*) P0,X(1),X(2)
  WRITE(3,5) P0,X(1),X(2)
  PAUSE 2
  DO 10 I=1,M
10  FF( I )=P0-X(1)*LOG(TT( I )/(X(2)+TT( I ))+1)
  DO 20 I=1,M
20  F( I )=U( I )-FF( I )
  NF=NF+1
5  FORMAT(F8.4,6X,F10.6,4X,F10.6)
  RETURN
  END

编辑过了，请大家禁用笑脸标记

iwtry

好多心呀！hoho^

sbn

这有什么可隐藏的啊，numeircal recipes里不就有现成的嘛

Faith165

这么差的东西还拿出来吓人。