source: Sophya/trunk/SophyaExt/CodeMinuit/code/mnline.F@ 2929

Last change on this file since 2929 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: 7.4 KB
Line 
1*
2* $Id: mnline.F,v 1.1.1.1 2003-06-11 14:18:28 cmv Exp $
3*
4* $Log: not supported by cvs2svn $
5* Revision 1.1.1.1 1996/03/07 14:31:30 mclareni
6* Minuit
7*
8*
9#include "minuit/pilot.h"
10 SUBROUTINE MNLINE(FCN,START,FSTART,STEP,SLOPE,TOLER,FUTIL)
11#include "minuit/d506dp.inc"
12CC Perform a line search from position START
13CC along direction STEP, where the length of vector STEP
14CC gives the expected position of minimum.
15CC FSTART is value of function at START
16CC SLOPE (if non-zero) is df/dx along STEP at START
17CC TOLER is initial tolerance of minimum in direction STEP
18#include "minuit/d506cm.inc"
19 EXTERNAL FCN,FUTIL
20 DIMENSION START(*), STEP(*)
21 PARAMETER (MAXPT=12)
22 DIMENSION XPQ(MAXPT),YPQ(MAXPT)
23 CHARACTER*1 CHPQ(MAXPT)
24 DIMENSION XVALS(3),FVALS(3),COEFF(3)
25 CHARACTER*26 CHARAL
26 CHARACTER*60 CMESS
27 PARAMETER (SLAMBG=5.,ALPHA=2.)
28C SLAMBG and ALPHA control the maximum individual steps allowed.
29C The first step is always =1. The max length of second step is SLAMBG.
30C The max size of subsequent steps is the maximum previous successful
31C step multiplied by ALPHA + the size of most recent successful step,
32C but cannot be smaller than SLAMBG.
33 LOGICAL LDEBUG
34 DATA CHARAL / 'ABCDEFGHIJKLMNOPQRSTUVWXYZ' /
35 LDEBUG = (IDBG(1).GE.1)
36C starting values for overall limits on total step SLAM
37 OVERAL = 1000.
38 UNDRAL = -100.
39C debug check if start is ok
40 IF (LDEBUG) THEN
41 CALL MNINEX(START)
42 CALL FCN(NPARX,GIN,F1,U,4,FUTIL)
43 NFCN=NFCN+1
44 IF (F1 .NE. FSTART) THEN
45 WRITE (ISYSWR,'(A/2E14.5/2X,10F10.5)')
46 + ' MNLINE start point not consistent, F values, parameters=',
47 + (X(KK),KK=1,NPAR)
48 ENDIF
49 ENDIF
50C . set up linear search along STEP
51
52 FVMIN = FSTART
53 XVMIN = ZERO
54 NXYPT = 1
55 CHPQ(1) = CHARAL(1:1)
56 XPQ(1) = 0.
57 YPQ(1) = FSTART
58C SLAMIN = smallest possible value of ABS(SLAM)
59 SLAMIN = ZERO
60 DO 20 I= 1, NPAR
61 IF (STEP(I) .EQ. ZERO) GO TO 20
62 RATIO = ABS(START(I)/STEP(I))
63 IF (SLAMIN .EQ. ZERO) SLAMIN = RATIO
64 IF (RATIO .LT. SLAMIN) SLAMIN = RATIO
65 20 X(I) = START(I) + STEP(I)
66 IF (SLAMIN .EQ. ZERO) SLAMIN = EPSMAC
67 SLAMIN = SLAMIN*EPSMA2
68 NPARX = NPAR
69C
70 CALL MNINEX(X)
71 CALL FCN(NPARX,GIN,F1,U,4,FUTIL)
72 NFCN=NFCN+1
73 NXYPT = NXYPT + 1
74 CHPQ(NXYPT) = CHARAL(NXYPT:NXYPT)
75 XPQ(NXYPT) = 1.
76 YPQ(NXYPT) = F1
77 IF (F1 .LT. FSTART) THEN
78 FVMIN = F1
79 XVMIN = 1.0
80 ENDIF
81C . quadr interp using slope GDEL and two points
82 SLAM = 1.
83 TOLER8 = TOLER
84 SLAMAX = SLAMBG
85 FLAST = F1
86C can iterate on two-points (cut) if no imprvmnt
87 25 CONTINUE
88 DENOM = 2.0*(FLAST-FSTART-SLOPE*SLAM)/SLAM**2
89C IF (DENOM .EQ. ZERO) DENOM = -0.1*SLOPE
90 SLAM = 1.
91 IF (DENOM .NE. ZERO) SLAM = -SLOPE/DENOM
92 IF (SLAM .LT. ZERO) SLAM = SLAMAX
93 IF (SLAM .GT. SLAMAX) SLAM = SLAMAX
94 IF (SLAM .LT. TOLER8) SLAM = TOLER8
95 IF (SLAM .LT. SLAMIN) GO TO 80
96 IF (ABS(SLAM-1.0).LT.TOLER8 .AND. F1.LT.FSTART) GO TO 70
97 IF (ABS(SLAM-1.0).LT.TOLER8) SLAM = 1.0+TOLER8
98 IF (NXYPT .GE. MAXPT) GO TO 65
99 DO 30 I= 1, NPAR
100 30 X(I) = START(I) + SLAM*STEP(I)
101 CALL MNINEX(X)
102 CALL FCN(NPAR,GIN,F2,U,4,FUTIL)
103 NFCN = NFCN + 1
104 NXYPT = NXYPT + 1
105 CHPQ(NXYPT) = CHARAL(NXYPT:NXYPT)
106 XPQ(NXYPT) = SLAM
107 YPQ(NXYPT) = F2
108 IF (F2 .LT. FVMIN) THEN
109 FVMIN = F2
110 XVMIN = SLAM
111 ENDIF
112 IF (FSTART .EQ. FVMIN) THEN
113 FLAST = F2
114 TOLER8 = TOLER*SLAM
115 OVERAL = SLAM-TOLER8
116 SLAMAX = OVERAL
117 GO TO 25
118 ENDIF
119C . quadr interp using 3 points
120 XVALS(1) = XPQ(1)
121 FVALS(1) = YPQ(1)
122 XVALS(2) = XPQ(NXYPT-1)
123 FVALS(2) = YPQ(NXYPT-1)
124 XVALS(3) = XPQ(NXYPT)
125 FVALS(3) = YPQ(NXYPT)
126C begin iteration, calculate desired step
127 50 CONTINUE
128 SLAMAX = MAX(SLAMAX,ALPHA*ABS(XVMIN))
129 CALL MNPFIT(XVALS,FVALS,3,COEFF,SDEV)
130 IF (COEFF(3) .LE. ZERO) THEN
131 SLOPEM = 2.0*COEFF(3)*XVMIN + COEFF(2)
132 IF (SLOPEM .LE. ZERO) THEN
133 SLAM = XVMIN + SLAMAX
134 ELSE
135 SLAM = XVMIN - SLAMAX
136 ENDIF
137 ELSE
138 SLAM = -COEFF(2)/(2.0*COEFF(3))
139 IF (SLAM .GT. XVMIN+SLAMAX) SLAM = XVMIN+SLAMAX
140 IF (SLAM .LT. XVMIN-SLAMAX) SLAM = XVMIN-SLAMAX
141 ENDIF
142 IF (SLAM .GT. ZERO) THEN
143 IF (SLAM .GT. OVERAL) SLAM = OVERAL
144 ELSE
145 IF (SLAM .LT. UNDRAL) SLAM = UNDRAL
146 ENDIF
147C come here if step was cut below
148 52 CONTINUE
149 TOLER9 = MAX(TOLER8,ABS(TOLER8*SLAM))
150 DO 55 IPT= 1, 3
151 IF (ABS(SLAM-XVALS(IPT)) .LT. TOLER9) GO TO 70
152 55 CONTINUE
153C take the step
154 IF (NXYPT .GE. MAXPT) GO TO 65
155 DO 60 I= 1, NPAR
156 60 X(I) = START(I)+SLAM*STEP(I)
157 CALL MNINEX(X)
158 CALL FCN(NPARX,GIN,F3,U,4,FUTIL)
159 NFCN = NFCN + 1
160 NXYPT = NXYPT + 1
161 CHPQ(NXYPT) = CHARAL(NXYPT:NXYPT)
162 XPQ(NXYPT) = SLAM
163 YPQ(NXYPT) = F3
164C find worst previous point out of three
165 FVMAX = FVALS(1)
166 NVMAX = 1
167 IF (FVALS(2) .GT. FVMAX) THEN
168 FVMAX = FVALS(2)
169 NVMAX = 2
170 ENDIF
171 IF (FVALS(3) .GT. FVMAX) THEN
172 FVMAX = FVALS(3)
173 NVMAX = 3
174 ENDIF
175C if latest point worse than all three previous, cut step
176 IF (F3 .GE. FVMAX) THEN
177 IF (NXYPT .GE. MAXPT) GO TO 65
178 IF (SLAM .GT. XVMIN) OVERAL = MIN(OVERAL,SLAM-TOLER8)
179 IF (SLAM .LT. XVMIN) UNDRAL = MAX(UNDRAL,SLAM+TOLER8)
180 SLAM = 0.5*(SLAM+XVMIN)
181 GO TO 52
182 ENDIF
183C prepare another iteration, replace worst previous point
184 XVALS(NVMAX) = SLAM
185 FVALS(NVMAX) = F3
186 IF (F3 .LT. FVMIN) THEN
187 FVMIN = F3
188 XVMIN = SLAM
189 ELSE
190 IF (SLAM .GT. XVMIN) OVERAL = MIN(OVERAL,SLAM-TOLER8)
191 IF (SLAM .LT. XVMIN) UNDRAL = MAX(UNDRAL,SLAM+TOLER8)
192 ENDIF
193 IF (NXYPT .LT. MAXPT) GO TO 50
194C . . end of iteration . . .
195C stop because too many iterations
196 65 CMESS = ' LINE SEARCH HAS EXHAUSTED THE LIMIT OF FUNCTION CALLS '
197 IF (LDEBUG) THEN
198 WRITE (ISYSWR,'(A/(2X,6G12.4))') ' MNLINE DEBUG: steps=',
199 + (STEP(KK),KK=1,NPAR)
200 ENDIF
201 GO TO 100
202C stop because within tolerance
203 70 CONTINUE
204 CMESS = ' LINE SEARCH HAS ATTAINED TOLERANCE '
205 GO TO 100
206 80 CONTINUE
207 CMESS = ' STEP SIZE AT ARITHMETICALLY ALLOWED MINIMUM'
208 100 CONTINUE
209 AMIN = FVMIN
210 DO 120 I= 1, NPAR
211 DIRIN(I) = STEP(I)*XVMIN
212 120 X(I) = START(I) + DIRIN(I)
213 CALL MNINEX(X)
214 IF (XVMIN .LT. 0.) CALL MNWARN('D','MNLINE',
215 + ' LINE MINIMUM IN BACKWARDS DIRECTION')
216 IF (FVMIN .EQ. FSTART) CALL MNWARN('D','MNLINE',
217 + ' LINE SEARCH FINDS NO IMPROVEMENT ')
218 IF (LDEBUG) THEN
219 WRITE (ISYSWR,'('' AFTER'',I3,'' POINTS,'',A)') NXYPT,CMESS
220 CALL MNPLOT(XPQ,YPQ,CHPQ,NXYPT,ISYSWR,NPAGWD,NPAGLN)
221 ENDIF
222 RETURN
223 END
Note: See TracBrowser for help on using the repository browser.