source: Sophya/trunk/SophyaExt/CodeMinuit/code/mncont.F@ 3632

Last change on this file since 3632 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: 8.7 KB
Line 
1*
2* $Id: mncont.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 MNCONT(FCN,KE1,KE2,NPTU,XPTU,YPTU,IERRF,FUTIL)
11#include "minuit/d506dp.inc"
12CC Find NPTU points along a contour where the function
13CC FMIN (X(KE1),X(KE2)) = AMIN+UP
14CC where FMIN is the minimum of FCN with respect to all
15CC the other NPAR-2 variable parameters (if any).
16CC IERRF on return will be equal to the number of points found:
17CC NPTU if normal termination with NPTU points found
18CC -1 if errors in the calling sequence (KE1, KE2 not variable)
19CC 0 if less than four points can be found (using MNMNOT)
20CC n>3 if only n points can be found (n < NPTU)
21CC
22#include "minuit/d506cm.inc"
23 DIMENSION XPTU(NPTU), YPTU(NPTU), W(MNI),GCC(MNI)
24 CHARACTER CHERE*10
25 PARAMETER (CHERE='MNContour ')
26 LOGICAL LDEBUG
27 EXTERNAL FCN,FUTIL
28C input arguments: parx, pary, devs, ngrid
29 LDEBUG = (IDBG(6) .GE. 1)
30 IF (KE1.LE.0 .OR. KE2.LE.0) GO TO 1350
31 IF (KE1.GT.NU .OR. KE2.GT.NU) GO TO 1350
32 KI1 = NIOFEX(KE1)
33 KI2 = NIOFEX(KE2)
34 IF (KI1.LE.0 .OR. KI2.LE.0) GO TO 1350
35 IF (KI1 .EQ. KI2) GO TO 1350
36 IF (NPTU .LT. 4) GO TO 1400
37C
38 NFCNCO = NFCN
39 NFCNMX = 100*(NPTU+5)*(NPAR+1)
40C The minimum
41 CALL MNCUVE(FCN,FUTIL)
42 U1MIN = U(KE1)
43 U2MIN = U(KE2)
44 IERRF = 0
45 CFROM = CHERE
46 NFCNFR = NFCNCO
47 IF (ISW(5) .GE. 0) THEN
48 WRITE (ISYSWR,'(1X,A,I4,A)')
49 + 'START MNCONTOUR CALCULATION OF',NPTU,' POINTS ON CONTOUR.'
50 IF (NPAR .GT. 2) THEN
51 IF (NPAR .EQ. 3) THEN
52 KI3 = 6 - KI1 - KI2
53 KE3 = NEXOFI(KI3)
54 WRITE (ISYSWR,'(1X,A,I3,2X,A)')
55 + 'EACH POINT IS A MINIMUM WITH RESPECT TO PARAMETER ',
56 + KE3, CPNAM(KE3)
57 ELSE
58 WRITE (ISYSWR,'(1X,A,I3,A)')
59 + 'EACH POINT IS A MINIMUM WITH RESPECT TO THE OTHER',
60 + NPAR-2, ' VARIABLE PARAMETERS.'
61 ENDIF
62 ENDIF
63 ENDIF
64C
65C Find the first four points using MNMNOT
66C ........................ first two points
67 CALL MNMNOT(FCN,KE1,KE2,VAL2PL,VAL2MI,FUTIL)
68 IF (ERN(KI1) .EQ. UNDEFI) THEN
69 XPTU(1) = ALIM(KE1)
70 CALL MNWARN('W',CHERE,'Contour squeezed by parameter limits.')
71 ELSE
72 IF (ERN(KI1) .GE. ZERO) GO TO 1500
73 XPTU(1) = U1MIN+ERN(KI1)
74 ENDIF
75 YPTU(1) = VAL2MI
76C
77 IF (ERP(KI1) .EQ. UNDEFI) THEN
78 XPTU(3) = BLIM(KE1)
79 CALL MNWARN('W',CHERE,'Contour squeezed by parameter limits.')
80 ELSE
81 IF (ERP(KI1) .LE. ZERO) GO TO 1500
82 XPTU(3) = U1MIN+ERP(KI1)
83 ENDIF
84 YPTU(3) = VAL2PL
85 SCALX = 1.0/(XPTU(3) - XPTU(1))
86C ........................... next two points
87 CALL MNMNOT(FCN,KE2,KE1,VAL2PL,VAL2MI,FUTIL)
88 IF (ERN(KI2) .EQ. UNDEFI) THEN
89 YPTU(2) = ALIM(KE2)
90 CALL MNWARN('W',CHERE,'Contour squeezed by parameter limits.')
91 ELSE
92 IF (ERN(KI2) .GE. ZERO) GO TO 1500
93 YPTU(2) = U2MIN+ERN(KI2)
94 ENDIF
95 XPTU(2) = VAL2MI
96 IF (ERP(KI2) .EQ. UNDEFI) THEN
97 YPTU(4) = BLIM(KE2)
98 CALL MNWARN('W',CHERE,'Contour squeezed by parameter limits.')
99 ELSE
100 IF (ERP(KI2) .LE. ZERO) GO TO 1500
101 YPTU(4) = U2MIN+ERP(KI2)
102 ENDIF
103 XPTU(4) = VAL2PL
104 SCALY = 1.0/(YPTU(4) - YPTU(2))
105 NOWPTS = 4
106 NEXT = 5
107 IF (LDEBUG) THEN
108 WRITE (ISYSWR,'(A)') ' Plot of four points found by MINOS'
109 XPT(1) = U1MIN
110 YPT(1) = U2MIN
111 CHPT(1) = ' '
112 NALL = MIN(NOWPTS+1,MAXCPT)
113 DO 85 I= 2, NALL
114 XPT(I) = XPTU(I-1)
115 YPT(I) = YPTU(I-1)
116 85 CONTINUE
117 CHPT(2)= 'A'
118 CHPT(3)= 'B'
119 CHPT(4)= 'C'
120 CHPT(5)= 'D'
121 CALL MNPLOT(XPT,YPT,CHPT,NALL,ISYSWR,NPAGWD,NPAGLN)
122 ENDIF
123C
124C ..................... save some values before fixing
125 ISW2 = ISW(2)
126 ISW4 = ISW(4)
127 SIGSAV = EDM
128 ISTRAV = ISTRAT
129 DC = DCOVAR
130 APSI = EPSI*0.5
131 ABEST=AMIN
132 MPAR=NPAR
133 NFMXIN = NFCNMX
134 DO 125 I= 1, MPAR
135 125 XT(I) = X(I)
136 DO 130 J= 1, MPAR*(MPAR+1)/2
137 130 VTHMAT(J) = VHMAT(J)
138 DO 135 I= 1, MPAR
139 GCC(I) = GLOBCC(I)
140 135 W(I) = WERR(I)
141C fix the two parameters in question
142 KINTS = NIOFEX(KE1)
143 CALL MNFIXP (KINTS,IERR)
144 KINTS = NIOFEX(KE2)
145 CALL MNFIXP (KINTS,IERR)
146C ......................Fill in the rest of the points
147 DO 900 INEW= NEXT, NPTU
148C find the two neighbouring points with largest separation
149 BIGDIS = 0.
150 DO 200 IOLD = 1, INEW-1
151 I2 = IOLD + 1
152 IF (I2 .EQ. INEW) I2 = 1
153 DIST = (SCALX*(XPTU(IOLD)-XPTU(I2)))**2 +
154 + (SCALY*(YPTU(IOLD)-YPTU(I2)))**2
155 IF (DIST .GT. BIGDIS) THEN
156 BIGDIS = DIST
157 IDIST = IOLD
158 ENDIF
159 200 CONTINUE
160 I1 = IDIST
161 I2 = I1 + 1
162 IF (I2 .EQ. INEW) I2 = 1
163C next point goes between I1 and I2
164 A1 = HALF
165 A2 = HALF
166 300 XMIDCR = A1*XPTU(I1) + A2*XPTU(I2)
167 YMIDCR = A1*YPTU(I1) + A2*YPTU(I2)
168 XDIR = YPTU(I2) - YPTU(I1)
169 YDIR = XPTU(I1) - XPTU(I2)
170 SCLFAC = MAX(ABS(XDIR*SCALX), ABS(YDIR*SCALY))
171 XDIRCR = XDIR/SCLFAC
172 YDIRCR = YDIR/SCLFAC
173 KE1CR = KE1
174 KE2CR = KE2
175C Find the contour crossing point along DIR
176 AMIN = ABEST
177 CALL MNCROS(FCN,AOPT,IERCR,FUTIL)
178 IF (IERCR .GT. 1) THEN
179C If cannot find mid-point, try closer to point 1
180 IF (A1 .GT. HALF) THEN
181 IF (ISW(5) .GE. 0)
182 + WRITE (ISYSWR,'(A,A,I3,A)') ' MNCONT CANNOT FIND NEXT',
183 + ' POINT ON CONTOUR. ONLY ',NOWPTS,' POINTS FOUND.'
184 GO TO 950
185 ENDIF
186 CALL MNWARN('W',CHERE,'Cannot find midpoint, try closer.')
187 A1 = 0.75
188 A2 = 0.25
189 GO TO 300
190 ENDIF
191C Contour has been located, insert new point in list
192 DO 830 MOVE= NOWPTS,I1+1,-1
193 XPTU(MOVE+1) = XPTU(MOVE)
194 YPTU(MOVE+1) = YPTU(MOVE)
195 830 CONTINUE
196 NOWPTS = NOWPTS + 1
197 XPTU(I1+1) = XMIDCR + XDIRCR*AOPT
198 YPTU(I1+1) = YMIDCR + YDIRCR*AOPT
199 900 CONTINUE
200 950 CONTINUE
201C
202 IERRF = NOWPTS
203 CSTATU = 'SUCCESSFUL'
204 IF (NOWPTS .LT. NPTU) CSTATU = 'INCOMPLETE'
205C make a lineprinter plot of the contour
206 IF (ISW(5) .GE. 0) THEN
207 XPT(1) = U1MIN
208 YPT(1) = U2MIN
209 CHPT(1) = ' '
210 NALL = MIN(NOWPTS+1,MAXCPT)
211 DO 1000 I= 2, NALL
212 XPT(I) = XPTU(I-1)
213 YPT(I) = YPTU(I-1)
214 CHPT(I)= 'X'
215 1000 CONTINUE
216 WRITE (ISYSWR,'(A,I3,2X,A)') ' Y-AXIS: PARAMETER ',KE2,
217 + CPNAM(KE2)
218 CALL MNPLOT(XPT,YPT,CHPT,NALL,ISYSWR,NPAGWD,NPAGLN)
219 WRITE (ISYSWR,'(25X,A,I3,2X,A)') 'X-AXIS: PARAMETER ',
220 + KE1,CPNAM(KE1)
221 ENDIF
222C print out the coordinates around the contour
223 IF (ISW(5) .GE. 1) THEN
224 NPCOL = (NOWPTS+1)/2
225 NFCOL = NOWPTS/2
226 WRITE (ISYSWR,'(/I5,A,G13.5,A,G11.3)') NOWPTS,
227 + ' POINTS ON CONTOUR. FMIN=',ABEST,' ERRDEF=',UP
228 WRITE (ISYSWR,'(9X,A,3X,A,18X,A,3X,A)')
229 + CPNAM(KE1),CPNAM(KE2),CPNAM(KE1),CPNAM(KE2)
230 DO 1050 LINE = 1, NFCOL
231 LR = LINE + NPCOL
232 WRITE (ISYSWR,'(1X,I5,2G13.5,10X,I5,2G13.5)')
233 + LINE,XPTU(LINE),YPTU(LINE),LR,XPTU(LR),YPTU(LR)
234 1050 CONTINUE
235 IF (NFCOL .LT. NPCOL) WRITE (ISYSWR,'(1X,I5,2G13.5)')
236 + NPCOL,XPTU(NPCOL),YPTU(NPCOL)
237 ENDIF
238C . . contour finished. reset v
239 ITAUR = 1
240 CALL MNFREE(1)
241 CALL MNFREE(1)
242 DO 1100 J= 1, MPAR*(MPAR+1)/2
243 1100 VHMAT(J) = VTHMAT(J)
244 DO 1120 I= 1, MPAR
245 GLOBCC(I) = GCC(I)
246 WERR(I) = W(I)
247 1120 X(I) = XT(I)
248 CALL MNINEX (X)
249 EDM = SIGSAV
250 AMIN = ABEST
251 ISW(2) = ISW2
252 ISW(4) = ISW4
253 DCOVAR = DC
254 ITAUR = 0
255 NFCNMX = NFMXIN
256 ISTRAT = ISTRAV
257 U(KE1) = U1MIN
258 U(KE2) = U2MIN
259 GO TO 2000
260C Error returns
261 1350 WRITE (ISYSWR,'(A)') ' INVALID PARAMETER NUMBERS.'
262 GO TO 1450
263 1400 WRITE (ISYSWR,'(A)') ' LESS THAN FOUR POINTS REQUESTED.'
264 1450 IERRF = -1
265 CSTATU = 'USER ERROR'
266 GO TO 2000
267 1500 WRITE (ISYSWR,'(A)') ' MNCONT UNABLE TO FIND FOUR POINTS.'
268 U(KE1) = U1MIN
269 U(KE2) = U2MIN
270 IERRF = 0
271 CSTATU = 'FAILED'
272 2000 CONTINUE
273 CFROM = CHERE
274 NFCNFR = NFCNCO
275 RETURN
276 END
Note: See TracBrowser for help on using the repository browser.