source: Sophya/trunk/SophyaExt/CodeMinuit/code/mnsimp.F@ 2863

Last change on this file since 2863 was 2403, checked in by cmv, 22 years ago

Creation du module de code source de MINUIT (CERNLIB) extrait par CMV

cmv 11/06/2003
File size: 6.2 KB
RevLine 
[2403]1*
2* $Id: mnsimp.F,v 1.1.1.1 2003-06-11 14:18:29 cmv Exp $
3*
4* $Log: not supported by cvs2svn $
5* Revision 1.2 1996/03/15 18:02:54 james
6* Modified Files:
7* mnderi.F eliminate possible division by zero
8* mnexcm.F suppress print on STOP when print flag=-1
9* set FVAL3 to flag if FCN already called with IFLAG=3
10* mninit.F set version 96.03
11* mnlims.F remove arguments, not needed
12* mnmigr.F VLEN -> LENV in debug print statement
13* mnparm.F move call to MNRSET to after NPAR redefined, to zero all
14* mnpsdf.F eliminate possible division by zero
15* mnscan.F suppress printout when print flag =-1
16* mnset.F remove arguments in call to MNLIMS
17* mnsimp.F fix CSTATU so status is PROGRESS only if new minimum
18* mnvert.F eliminate possible division by zero
19*
20* Revision 1.1.1.1 1996/03/07 14:31:31 mclareni
21* Minuit
22*
23*
24#include "minuit/pilot.h"
25 SUBROUTINE MNSIMP(FCN,FUTIL)
26#include "minuit/d506dp.inc"
27CC Performs a minimization using the simplex method of Nelder
28CC and Mead (ref. -- Comp. J. 7,308 (1965)).
29CC
30#include "minuit/d506cm.inc"
31 EXTERNAL FCN,FUTIL
32 DIMENSION Y(MNI+1)
33 DATA ALPHA,BETA,GAMMA,RHOMIN,RHOMAX / 1.0, 0.5, 2.0, 4.0, 8.0/
34 IF (NPAR .LE. 0) RETURN
35 IF (AMIN .EQ. UNDEFI) CALL MNAMIN(FCN,FUTIL)
36 CFROM = 'SIMPLEX '
37 NFCNFR = NFCN
38 CSTATU= 'UNCHANGED '
39 NPFN=NFCN
40 NPARP1=NPAR+1
41 NPARX = NPAR
42 RHO1 = 1.0 + ALPHA
43 RHO2 = RHO1 + ALPHA*GAMMA
44 WG = 1.0/FLOAT(NPAR)
45 IF (ISW(5) .GE. 0) WRITE(ISYSWR,100) EPSI
46 100 FORMAT(' START SIMPLEX MINIMIZATION. CONVERGENCE WHEN EDM .LT.'
47 +,E10.2 )
48 DO 2 I= 1, NPAR
49 DIRIN(I) = WERR(I)
50 CALL MNDXDI(X(I),I,DXDI)
51 IF (DXDI .NE. ZERO) DIRIN(I)=WERR(I)/DXDI
52 DMIN = EPSMA2*ABS(X(I))
53 IF (DIRIN(I) .LT. DMIN) DIRIN(I)=DMIN
54 2 CONTINUE
55C** choose the initial simplex using single-parameter searches
56 1 CONTINUE
57 YNPP1 = AMIN
58 JL = NPARP1
59 Y(NPARP1) = AMIN
60 ABSMIN = AMIN
61 DO 10 I= 1, NPAR
62 AMING = AMIN
63 PBAR(I) = X(I)
64 BESTX = X(I)
65 KG = 0
66 NS = 0
67 NF = 0
68 4 X(I) = BESTX + DIRIN(I)
69 CALL MNINEX(X)
70 CALL FCN(NPARX,GIN, F, U, 4, FUTIL)
71 NFCN = NFCN + 1
72 IF (F .LT. AMING) GO TO 6
73C failure
74 IF (KG .EQ. 1) GO TO 8
75 KG = -1
76 NF = NF + 1
77 DIRIN(I) = DIRIN(I) * (-0.4)
78 IF (NF .LT. 3) GO TO 4
79C stop after three failures
80 BESTX = X(I)
81 DIRIN(I) = DIRIN(I) * 3.0
82 AMING = F
83 GO TO 8
84C
85C success
86 6 BESTX = X(I)
87 DIRIN(I) = DIRIN(I) * 3.0
88 AMING = F
89 CSTATU= 'PROGRESS '
90 KG = 1
91 NS = NS + 1
92 IF (NS .LT. 6) GO TO 4
93C
94C 3 failures or 6 successes or
95C local minimum found in ith direction
96 8 Y(I) = AMING
97 IF (AMING .LT. ABSMIN) JL = I
98 IF (AMING .LT. ABSMIN) ABSMIN = AMING
99 X(I) = BESTX
100 DO 9 K= 1, NPAR
101 9 P(K,I) = X(K)
102 10 CONTINUE
103 JH = NPARP1
104 AMIN=Y(JL)
105 CALL MNRAZZ(YNPP1,PBAR,Y,JH,JL)
106 DO 20 I= 1, NPAR
107 20 X(I) = P(I,JL)
108 CALL MNINEX(X)
109 IF (ISW(5) .GE. 1) CALL MNPRIN(5,AMIN)
110 EDM = BIGEDM
111 SIG2 = EDM
112 NCYCL=0
113C . . . . . start main loop
114 50 CONTINUE
115 IF (SIG2 .LT. EPSI .AND. EDM.LT.EPSI) GO TO 76
116 SIG2 = EDM
117 IF ((NFCN-NPFN) .GT. NFCNMX) GO TO 78
118C calculate new point * by reflection
119 DO 60 I= 1, NPAR
120 PB = 0.
121 DO 59 J= 1, NPARP1
122 59 PB = PB + WG * P(I,J)
123 PBAR(I) = PB - WG * P(I,JH)
124 60 PSTAR(I)=(1.+ALPHA)*PBAR(I)-ALPHA*P(I,JH)
125 CALL MNINEX(PSTAR)
126 CALL FCN(NPARX,GIN,YSTAR,U,4,FUTIL)
127 NFCN=NFCN+1
128 IF(YSTAR.GE.AMIN) GO TO 70
129C point * better than jl, calculate new point **
130 CSTATU = 'PROGRESS '
131 DO 61 I=1,NPAR
132 61 PSTST(I)=GAMMA*PSTAR(I)+(1.-GAMMA)*PBAR(I)
133 CALL MNINEX(PSTST)
134 CALL FCN(NPARX,GIN,YSTST,U,4,FUTIL)
135 NFCN=NFCN+1
136C try a parabola through ph, pstar, pstst. min = prho
137 Y1 = (YSTAR-Y(JH)) * RHO2
138 Y2 = (YSTST-Y(JH)) * RHO1
139 RHO = 0.5 * (RHO2*Y1 -RHO1*Y2) / (Y1 -Y2)
140 IF (RHO .LT. RHOMIN) GO TO 66
141 IF (RHO .GT. RHOMAX) RHO = RHOMAX
142 DO 64 I= 1, NPAR
143 64 PRHO(I) = RHO*PBAR(I) + (1.0-RHO)*P(I,JH)
144 CALL MNINEX(PRHO)
145 CALL FCN(NPARX,GIN,YRHO, U,4,FUTIL)
146 NFCN = NFCN + 1
147 IF (YRHO .LT. AMIN) CSTATU = 'PROGRESS '
148 IF (YRHO .LT. Y(JL) .AND. YRHO .LT. YSTST) GO TO 65
149 IF (YSTST .LT. Y(JL)) GO TO 67
150 IF (YRHO .GT. Y(JL)) GO TO 66
151C accept minimum point of parabola, PRHO
152 65 CALL MNRAZZ (YRHO,PRHO,Y,JH,JL)
153 GO TO 68
154 66 IF (YSTST .LT. Y(JL)) GO TO 67
155 CALL MNRAZZ(YSTAR,PSTAR,Y,JH,JL)
156 GO TO 68
157 67 CALL MNRAZZ(YSTST,PSTST,Y,JH,JL)
158 68 NCYCL=NCYCL+1
159 IF (ISW(5) .LT. 2) GO TO 50
160 IF (ISW(5) .GE. 3 .OR. MOD(NCYCL, 10) .EQ. 0) CALL MNPRIN(5,AMIN)
161 GO TO 50
162C point * is not as good as jl
163 70 IF (YSTAR .GE. Y(JH)) GO TO 73
164 JHOLD = JH
165 CALL MNRAZZ(YSTAR,PSTAR,Y,JH,JL)
166 IF (JHOLD .NE. JH) GO TO 50
167C calculate new point **
168 73 DO 74 I=1,NPAR
169 74 PSTST(I)=BETA*P(I,JH)+(1.-BETA)*PBAR(I)
170 CALL MNINEX (PSTST)
171 CALL FCN(NPARX,GIN,YSTST,U,4,FUTIL)
172 NFCN=NFCN+1
173 IF(YSTST.GT.Y(JH)) GO TO 1
174C point ** is better than jh
175 IF (YSTST .LT. AMIN) CSTATU = 'PROGRESS '
176 IF (YSTST .LT. AMIN) GO TO 67
177 CALL MNRAZZ(YSTST,PSTST,Y,JH,JL)
178 GO TO 50
179C . . . . . . end main loop
180 76 IF (ISW(5) .GE. 0) WRITE(ISYSWR,'(A)')
181 + ' SIMPLEX MINIMIZATION HAS CONVERGED.'
182 ISW(4) = 1
183 GO TO 80
184 78 IF (ISW(5) .GE. 0) WRITE(ISYSWR,'(A)')
185 + ' SIMPLEX TERMINATES WITHOUT CONVERGENCE.'
186 CSTATU= 'CALL LIMIT'
187 ISW(4) = -1
188 ISW(1) = 1
189 80 DO 82 I=1,NPAR
190 PB = 0.
191 DO 81 J=1,NPARP1
192 81 PB = PB + WG * P(I,J)
193 82 PBAR(I) = PB - WG * P(I,JH)
194 CALL MNINEX(PBAR)
195 CALL FCN(NPARX,GIN,YPBAR,U,4,FUTIL)
196 NFCN=NFCN+1
197 IF (YPBAR .LT. AMIN) CALL MNRAZZ(YPBAR,PBAR,Y,JH,JL)
198 CALL MNINEX(X)
199 IF (NFCNMX+NPFN-NFCN .LT. 3*NPAR) GO TO 90
200 IF (EDM .GT. 2.0*EPSI) GO TO 1
201 90 IF (ISW(5) .GE. 0) CALL MNPRIN(5, AMIN)
202 RETURN
203 END
Note: See TracBrowser for help on using the repository browser.