source: Sophya/trunk/SophyaExt/CodeMinuit/code/mncros.F@ 3056

Last change on this file since 3056 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: 11.0 KB
Line 
1*
2* $Id: mncros.F,v 1.1.1.1 2003-06-11 14:18:26 cmv Exp $
3*
4* $Log: not supported by cvs2svn $
5* Revision 1.1.1.1 1996/03/07 14:31:29 mclareni
6* Minuit
7*
8*
9#include "minuit/pilot.h"
10 SUBROUTINE MNCROS(FCN,AOPT,IERCR,FUTIL)
11#include "minuit/d506dp.inc"
12CC Find point where MNEVAL=AMIN+UP, along the line through
13CC XMIDCR,YMIDCR with direction XDIRCR,YDIRCR, where X and Y
14CC are parameters KE1CR and KE2CR. If KE2CR=0 (from MINOS),
15CC only KE1CR is varied. From MNCONT, both are varied.
16CC Crossing point is at
17CC (U(KE1),U(KE2)) = (XMID,YMID) + AOPT*(XDIR,YDIR)
18CC
19#include "minuit/d506cm.inc"
20 CHARACTER CHERE*10, CHARAL*28, CHSIGN*4
21 PARAMETER (CHERE='MNCROS ', MLSB=3, MAXITR=15, TLR=0.01)
22 DIMENSION FLSB(MLSB),ALSB(MLSB), COEFF(3)
23 LOGICAL LDEBUG
24 EXTERNAL FCN,FUTIL
25 DATA CHARAL/' .ABCDEFGHIJKLMNOPQRSTUVWXYZ'/
26 LDEBUG = (IDBG(6) .GE. 1)
27 AMINSV = AMIN
28C convergence when F is within TLF of AIM and next prediction
29C of AOPT is within TLA of previous value of AOPT
30 AIM = AMIN + UP
31 TLF = TLR*UP
32 TLA = TLR
33 XPT(1) = 0.0
34 YPT(1) = AIM
35 CHPT(1) = ' '
36 IPT = 1
37 IF (KE2CR .EQ. 0) THEN
38 XPT(2) = -1.0
39 YPT(2) = AMIN
40 CHPT(2) = '.'
41 IPT = 2
42 ENDIF
43C find the largest allowed A
44 AULIM = 100.
45 DO 100 IK= 1, 2
46 IF (IK .EQ. 1) THEN
47 KEX = KE1CR
48 ZMID = XMIDCR
49 ZDIR = XDIRCR
50 ELSE
51 IF (KE2CR .EQ. 0) GO TO 100
52 KEX = KE2CR
53 ZMID = YMIDCR
54 ZDIR = YDIRCR
55 ENDIF
56 IF (NVARL(KEX) .LE. 1) GO TO 100
57 IF (ZDIR .EQ. ZERO) GO TO 100
58 ZLIM = ALIM(KEX)
59 IF (ZDIR .GT. ZERO) ZLIM = BLIM(KEX)
60 AULIM = MIN(AULIM,(ZLIM-ZMID)/ZDIR)
61 100 CONTINUE
62C LSB = Line Search Buffer
63C first point
64 ANEXT = 0.
65 AOPT = ANEXT
66 LIMSET = .FALSE.
67 IF (AULIM .LT. AOPT+TLA) LIMSET = .TRUE.
68 CALL MNEVAL(FCN,ANEXT,FNEXT,IEREV,FUTIL)
69C debug printout:
70 IF (LDEBUG) WRITE (ISYSWR,'(A,I8,A,F10.5,A,2F10.5)')
71 + ' MNCROS: calls=',NFCN,' AIM=',AIM,' F,A=',FNEXT,AOPT
72 IF (IEREV .GT. 0) GO TO 900
73 IF (LIMSET .AND. FNEXT .LE. AIM) GO TO 930
74 IPT = IPT + 1
75 XPT(IPT) = ANEXT
76 YPT(IPT) = FNEXT
77 CHPT(IPT)= CHARAL(IPT:IPT)
78 ALSB(1) = ANEXT
79 FLSB(1) = FNEXT
80 FNEXT = MAX(FNEXT,AMINSV+0.1*UP)
81 AOPT = SQRT((UP)/(FNEXT-AMINSV)) - 1.0
82 IF (ABS(FNEXT-AIM) .LT. TLF) GO TO 800
83C
84 IF (AOPT .LT. -HALF) AOPT = -HALF
85 IF (AOPT .GT. ONE) AOPT = ONE
86 LIMSET = .FALSE.
87 IF (AOPT .GT. AULIM) THEN
88 AOPT = AULIM
89 LIMSET = .TRUE.
90 ENDIF
91 CALL MNEVAL(FCN,AOPT,FNEXT,IEREV,FUTIL)
92C debug printout:
93 IF (LDEBUG) WRITE (ISYSWR,'(A,I8,A,F10.5,A,2F10.5)')
94 + ' MNCROS: calls=',NFCN,' AIM=',AIM,' F,A=',FNEXT,AOPT
95 IF (IEREV .GT. 0) GO TO 900
96 IF (LIMSET .AND. FNEXT .LE. AIM) GO TO 930
97 ALSB(2) = AOPT
98 IPT = IPT + 1
99 XPT(IPT) = ALSB(2)
100 YPT(IPT) = FNEXT
101 CHPT(IPT)= CHARAL(IPT:IPT)
102 FLSB(2) = FNEXT
103 DFDA = (FLSB(2)-FLSB(1))/ (ALSB(2)-ALSB(1))
104C DFDA must be positive on the contour
105 IF (DFDA .GT. ZERO) GO TO 460
106 300 CALL MNWARN('D',CHERE,'Looking for slope of the right sign')
107 MAXLK = MAXITR - IPT
108 DO 400 IT= 1, MAXLK
109 ALSB(1) = ALSB(2)
110 FLSB(1) = FLSB(2)
111 AOPT = ALSB(1) + 0.2*REAL(IT)
112 LIMSET = .FALSE.
113 IF (AOPT .GT. AULIM) THEN
114 AOPT = AULIM
115 LIMSET = .TRUE.
116 ENDIF
117 CALL MNEVAL(FCN,AOPT,FNEXT,IEREV,FUTIL)
118C debug printout:
119 IF (LDEBUG) WRITE (ISYSWR,'(A,I8,A,F10.5,A,2F10.5)')
120 + ' MNCROS: calls=',NFCN,' AIM=',AIM,' F,A=',FNEXT,AOPT
121 IF (IEREV .GT. 0) GO TO 900
122 IF (LIMSET .AND. FNEXT .LE. AIM) GO TO 930
123 ALSB(2) = AOPT
124 IPT = IPT + 1
125 XPT(IPT) = ALSB(2)
126 YPT(IPT) = FNEXT
127 CHPT(IPT)= CHARAL(IPT:IPT)
128 FLSB(2) = FNEXT
129 DFDA = (FLSB(2)-FLSB(1))/ (ALSB(2)-ALSB(1))
130 IF (DFDA .GT. ZERO) GO TO 450
131 400 CONTINUE
132 CALL MNWARN('W',CHERE,'Cannot find slope of the right sign')
133 GO TO 950
134 450 CONTINUE
135C we have two points with the right slope
136 460 AOPT = ALSB(2) + (AIM-FLSB(2))/DFDA
137 FDIST = MIN(ABS(AIM -FLSB(1)),ABS(AIM -FLSB(2)))
138 ADIST = MIN(ABS(AOPT-ALSB(1)),ABS(AOPT-ALSB(2)))
139 TLA = TLR
140 IF (ABS(AOPT) .GT. ONE) TLA = TLR*ABS(AOPT)
141 IF (ADIST .LT. TLA .AND. FDIST .LT. TLF) GO TO 800
142 IF (IPT .GE. MAXITR) GO TO 950
143 BMIN = MIN(ALSB(1),ALSB(2)) - 1.0
144 IF (AOPT .LT. BMIN) AOPT = BMIN
145 BMAX = MAX(ALSB(1),ALSB(2)) + 1.0
146 IF (AOPT .GT. BMAX) AOPT = BMAX
147C Try a third point
148 LIMSET = .FALSE.
149 IF (AOPT .GT. AULIM) THEN
150 AOPT = AULIM
151 LIMSET = .TRUE.
152 ENDIF
153 CALL MNEVAL(FCN,AOPT,FNEXT,IEREV,FUTIL)
154C debug printout:
155 IF (LDEBUG) WRITE (ISYSWR,'(A,I8,A,F10.5,A,2F10.5)')
156 + ' MNCROS: calls=',NFCN,' AIM=',AIM,' F,A=',FNEXT,AOPT
157 IF (IEREV .GT. 0) GO TO 900
158 IF (LIMSET .AND. FNEXT .LE. AIM) GO TO 930
159 ALSB(3) = AOPT
160 IPT = IPT + 1
161 XPT(IPT) = ALSB(3)
162 YPT(IPT) = FNEXT
163 CHPT(IPT)= CHARAL(IPT:IPT)
164 FLSB(3) = FNEXT
165 INEW = 3
166C now we have three points, ask how many <AIM
167 ECARMN = ABS(FNEXT-AIM)
168 IBEST = 3
169 ECARMX = 0.
170 NOLESS = 0
171 DO 480 I= 1, 3
172 ECART = ABS(FLSB(I) - AIM)
173 IF (ECART .GT. ECARMX) THEN
174 ECARMX = ECART
175 IWORST = I
176 ENDIF
177 IF (ECART .LT. ECARMN) THEN
178 ECARMN = ECART
179 IBEST = I
180 ENDIF
181 IF (FLSB(I) .LT. AIM) NOLESS = NOLESS + 1
182 480 CONTINUE
183 INEW = IBEST
184C if at least one on each side of AIM, fit a parabola
185 IF (NOLESS.EQ.1 .OR. NOLESS.EQ.2) GO TO 500
186C if all three are above AIM, third must be closest to AIM
187 IF (NOLESS .EQ. 0 .AND. IBEST .NE. 3) GO TO 950
188C if all three below, and third is not best, then slope
189C has again gone negative, look for positive slope.
190 IF (NOLESS .EQ. 3 .AND. IBEST .NE. 3) THEN
191 ALSB(2) = ALSB(3)
192 FLSB(2) = FLSB(3)
193 GO TO 300
194 ENDIF
195C in other cases, new straight line thru last two points
196 ALSB(IWORST) = ALSB(3)
197 FLSB(IWORST) = FLSB(3)
198 DFDA = (FLSB(2)-FLSB(1))/ (ALSB(2)-ALSB(1))
199 GO TO 460
200C parabola fit
201 500 CALL MNPFIT(ALSB,FLSB,3,COEFF,SDEV)
202 IF (COEFF(3) .LE. ZERO) CALL MNWARN ('D',CHERE,
203 + 'Curvature is negative near contour line.')
204 DETERM = COEFF(2)**2 - 4.*COEFF(3)*(COEFF(1)-AIM)
205 IF (DETERM .LE. ZERO) THEN
206 CALL MNWARN('D',CHERE,'Problem 2, impossible determinant')
207 GO TO 950
208 ENDIF
209C Find which root is the right one
210 RT = SQRT(DETERM)
211 X1 = (-COEFF(2) + RT)/(2.*COEFF(3))
212 X2 = (-COEFF(2) - RT)/(2.*COEFF(3))
213 S1 = COEFF(2) + 2.*X1*COEFF(3)
214 S2 = COEFF(2) + 2.*X2*COEFF(3)
215 IF (S1*S2 .GT. ZERO) WRITE (ISYSWR,'(A)') ' MNCONTour problem 1'
216 AOPT = X1
217 SLOPE = S1
218 IF (S2 .GT. ZERO) THEN
219 AOPT = X2
220 SLOPE = S2
221 ENDIF
222C ask if converged
223 TLA = TLR
224 IF (ABS(AOPT) .GT. ONE) TLA = TLR*ABS(AOPT)
225 IF (ABS(AOPT-ALSB(IBEST)) .LT. TLA .AND.
226 & ABS(FLSB(IBEST)-AIM) .LT. TLF) GO TO 800
227 IF (IPT .GE. MAXITR) GO TO 950
228C see if proposed point is in acceptable zone between L and R
229C first find ILEFT, IRIGHT, IOUT and IBEST
230 ILEFT = 0
231 IRIGHT = 0
232 IBEST = 1
233 ECARMX = 0.
234 ECARMN = ABS(AIM-FLSB(1))
235 DO 550 I= 1, 3
236 ECART = ABS(FLSB(I) - AIM)
237 IF (ECART .LT. ECARMN) THEN
238 ECARMN = ECART
239 IBEST = I
240 ENDIF
241 IF (ECART .GT. ECARMX) ECARMX = ECART
242 IF (FLSB(I) .GT. AIM) THEN
243 IF (IRIGHT .EQ. 0) THEN
244 IRIGHT = I
245 ELSE IF (FLSB(I) .GT. FLSB(IRIGHT)) THEN
246 IOUT = I
247 ELSE
248 IOUT = IRIGHT
249 IRIGHT = I
250 ENDIF
251 ELSE IF (ILEFT .EQ. 0) THEN
252 ILEFT = I
253 ELSE IF (FLSB(I) .LT. FLSB(ILEFT)) THEN
254 IOUT = I
255 ELSE
256 IOUT = ILEFT
257 ILEFT = I
258 ENDIF
259 550 CONTINUE
260C avoid keeping a very bad point next time around
261 IF (ECARMX .GT. 10.*ABS(FLSB(IOUT)-AIM))
262 & AOPT = HALF*AOPT + HALF*HALF*(ALSB(IRIGHT)+ALSB(ILEFT))
263C knowing ILEFT and IRIGHT, get acceptable window
264 SMALLA = 0.1*TLA
265 IF (SLOPE*SMALLA .GT. TLF) SMALLA = TLF/SLOPE
266 ALEFT = ALSB(ILEFT) + SMALLA
267 ARIGHT = ALSB(IRIGHT) - SMALLA
268C move proposed point AOPT into window if necessary
269 IF (AOPT .LT. ALEFT) AOPT = ALEFT
270 IF (AOPT .GT. ARIGHT) AOPT = ARIGHT
271 IF (ALEFT .GT. ARIGHT)AOPT = HALF*(ALEFT + ARIGHT)
272C see if proposed point outside limits (should be impossible!)
273 LIMSET = .FALSE.
274 IF (AOPT .GT. AULIM) THEN
275 AOPT = AULIM
276 LIMSET = .TRUE.
277 ENDIF
278C Evaluate function at new point AOPT
279 CALL MNEVAL(FCN,AOPT,FNEXT,IEREV,FUTIL)
280C debug printout:
281 IF (LDEBUG) WRITE (ISYSWR,'(A,I8,A,F10.5,A,2F10.5)')
282 + ' MNCROS: calls=',NFCN,' AIM=',AIM,' F,A=',FNEXT,AOPT
283 IF (IEREV .GT. 0) GO TO 900
284 IF (LIMSET .AND. FNEXT .LE. AIM) GO TO 930
285 IPT = IPT + 1
286 XPT(IPT) = AOPT
287 YPT(IPT) = FNEXT
288 CHPT(IPT)= CHARAL(IPT:IPT)
289C Replace odd point by new one
290 ALSB(IOUT) = AOPT
291 FLSB(IOUT) = FNEXT
292C the new point may not be the best, but it is the only one
293C which could be good enough to pass convergence criteria
294 IBEST = IOUT
295 GO TO 500
296C
297C Contour has been located, return point to MNCONT OR MINOS
298 800 CONTINUE
299 IERCR = 0
300 GO TO 1000
301C error in the minimization
302 900 IF (IEREV .EQ. 1) GO TO 940
303 GO TO 950
304C parameter up against limit
305 930 IERCR = 1
306 GO TO 1000
307C too many calls to FCN
308 940 IERCR = 2
309 GO TO 1000
310C cannot find next point
311 950 IERCR = 3
312C in any case
313 1000 CONTINUE
314 IF (LDEBUG) THEN
315 ITOOHI = 0
316 DO 1100 I= 1, IPT
317 IF (YPT(I) .GT. AIM+UP) THEN
318 YPT(I) = AIM+UP
319 CHPT(I) = '+'
320 ITOOHI = 1
321 ENDIF
322 1100 CONTINUE
323 CHSIGN = 'POSI'
324 IF (XDIRCR .LT. ZERO) CHSIGN = 'NEGA'
325 IF (KE2CR .EQ. 0) WRITE (ISYSWR, '(2X,A,A,I3)')
326 + CHSIGN,'TIVE MINOS ERROR, PARAMETER ',KE1CR
327 IF (ITOOHI .EQ. 1) WRITE (ISYSWR, '(10X,A)')
328 + 'POINTS LABELLED "+" WERE TOO HIGH TO PLOT.'
329 IF (IERCR .EQ. 1) WRITE (ISYSWR,'(10X,A)')
330 + 'RIGHTMOST POINT IS UP AGAINST LIMIT.'
331 CALL MNPLOT(XPT,YPT,CHPT,IPT,ISYSWR,NPAGWD,NPAGLN)
332 ENDIF
333 RETURN
334 END
Note: See TracBrowser for help on using the repository browser.