source: Sophya/trunk/SophyaExt/CodeMinuit/code/mnimpr.F@ 3005

Last change on this file since 3005 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: 5.8 KB
Line 
1*
2* $Id: mnimpr.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 MNIMPR(FCN,FUTIL)
11#include "minuit/d506dp.inc"
12CC Attempts to improve on a good local minimum by finding a
13CC better one. The quadratic part of FCN is removed by MNCALF
14CC and this transformed function is minimized using the simplex
15CC method from several random starting points.
16CC ref. -- Goldstein and Price, Math.Comp. 25, 569 (1971)
17CC
18#include "minuit/d506cm.inc"
19 EXTERNAL FCN,FUTIL
20 DIMENSION DSAV(MNI), Y(MNI+1)
21 PARAMETER (ALPHA=1.,BETA=0.5,GAMMA=2.0)
22 DATA RNUM/0./
23 IF (NPAR .LE. 0) RETURN
24 IF (AMIN .EQ. UNDEFI) CALL MNAMIN(FCN,FUTIL)
25 CSTATU = 'UNCHANGED '
26 ITAUR = 1
27 EPSI = 0.1*UP
28 NPFN=NFCN
29 NLOOP = WORD7(2)
30 IF (NLOOP .LE. 0) NLOOP = NPAR + 4
31 NPARX = NPAR
32 NPARP1=NPAR+1
33 WG = 1.0/FLOAT(NPAR)
34 SIGSAV = EDM
35 APSI = AMIN
36 DO 2 I= 1, NPAR
37 XT(I) = X(I)
38 DSAV(I) = WERR(I)
39 DO 2 J = 1, I
40 NDEX = I*(I-1)/2 + J
41 P(I,J) = VHMAT(NDEX)
42 2 P(J,I) = P(I,J)
43 CALL MNVERT(P,MAXINT,MAXINT,NPAR,IFAIL)
44 IF (IFAIL .GE. 1) GO TO 280
45C Save inverted matrix in VT
46 DO 12 I= 1, NPAR
47 NDEX = I*(I-1)/2
48 DO 12 J= 1, I
49 NDEX = NDEX + 1
50 12 VTHMAT(NDEX) = P(I,J)
51 LOOP = 0
52C
53 20 CONTINUE
54 DO 25 I= 1, NPAR
55 DIRIN(I) = 2.0*DSAV(I)
56 CALL MNRN15(RNUM,ISEED)
57 25 X(I) = XT(I) + 2.0*DIRIN(I)*(RNUM-0.5)
58 LOOP = LOOP + 1
59 REG = 2.0
60 IF (ISW(5) .GE. 0) WRITE (ISYSWR, 1040) LOOP
61 30 CALL MNCALF(FCN,X,YCALF,FUTIL)
62 AMIN = YCALF
63C . . . . set up random simplex
64 JL = NPARP1
65 JH = NPARP1
66 Y(NPARP1) = AMIN
67 AMAX = AMIN
68 DO 45 I= 1, NPAR
69 XI = X(I)
70 CALL MNRN15(RNUM,ISEED)
71 X(I) = XI - DIRIN(I) *(RNUM-0.5)
72 CALL MNCALF(FCN,X,YCALF,FUTIL)
73 Y(I) = YCALF
74 IF (Y(I) .LT. AMIN) THEN
75 AMIN = Y(I)
76 JL = I
77 ELSE IF (Y(I) .GT. AMAX) THEN
78 AMAX = Y(I)
79 JH = I
80 ENDIF
81 DO 40 J= 1, NPAR
82 40 P(J,I) = X(J)
83 P(I,NPARP1) = XI
84 X(I) = XI
85 45 CONTINUE
86C
87 EDM = AMIN
88 SIG2 = EDM
89C . . . . . . . start main loop
90 50 CONTINUE
91 IF (AMIN .LT. ZERO) GO TO 95
92 IF (ISW(2) .LE. 2) GO TO 280
93 EP = 0.1*AMIN
94 IF (SIG2 .LT. EP .AND. EDM.LT.EP ) GO TO 100
95 SIG2 = EDM
96 IF ((NFCN-NPFN) .GT. NFCNMX) GO TO 300
97C calculate new point * by reflection
98 DO 60 I= 1, NPAR
99 PB = 0.
100 DO 59 J= 1, NPARP1
101 59 PB = PB + WG * P(I,J)
102 PBAR(I) = PB - WG * P(I,JH)
103 60 PSTAR(I)=(1.+ALPHA)*PBAR(I)-ALPHA*P(I,JH)
104 CALL MNCALF(FCN,PSTAR,YCALF,FUTIL)
105 YSTAR = YCALF
106 IF(YSTAR.GE.AMIN) GO TO 70
107C point * better than jl, calculate new point **
108 DO 61 I=1,NPAR
109 61 PSTST(I)=GAMMA*PSTAR(I)+(1.-GAMMA)*PBAR(I)
110 CALL MNCALF(FCN,PSTST,YCALF,FUTIL)
111 YSTST = YCALF
112 66 IF (YSTST .LT. Y(JL)) GO TO 67
113 CALL MNRAZZ(YSTAR,PSTAR,Y,JH,JL)
114 GO TO 50
115 67 CALL MNRAZZ(YSTST,PSTST,Y,JH,JL)
116 GO TO 50
117C point * is not as good as jl
118 70 IF (YSTAR .GE. Y(JH)) GO TO 73
119 JHOLD = JH
120 CALL MNRAZZ(YSTAR,PSTAR,Y,JH,JL)
121 IF (JHOLD .NE. JH) GO TO 50
122C calculate new point **
123 73 DO 74 I=1,NPAR
124 74 PSTST(I)=BETA*P(I,JH)+(1.-BETA)*PBAR(I)
125 CALL MNCALF(FCN,PSTST,YCALF,FUTIL)
126 YSTST = YCALF
127 IF(YSTST.GT.Y(JH)) GO TO 30
128C point ** is better than jh
129 IF (YSTST .LT. AMIN) GO TO 67
130 CALL MNRAZZ(YSTST,PSTST,Y,JH,JL)
131 GO TO 50
132C . . . . . . end main loop
133 95 IF (ISW(5) .GE. 0) WRITE (ISYSWR,1000)
134 REG = 0.1
135C . . . . . ask if point is new
136 100 CALL MNINEX(X)
137 CALL FCN(NPARX,GIN,AMIN,U,4,FUTIL)
138 NFCN = NFCN + 1
139 DO 120 I= 1, NPAR
140 DIRIN(I) = REG*DSAV(I)
141 IF (ABS(X(I)-XT(I)) .GT. DIRIN(I)) GO TO 150
142 120 CONTINUE
143 GO TO 230
144 150 NFCNMX = NFCNMX + NPFN - NFCN
145 NPFN = NFCN
146 CALL MNSIMP(FCN,FUTIL)
147 IF (AMIN .GE. APSI) GO TO 325
148 DO 220 I= 1, NPAR
149 DIRIN(I) = 0.1 *DSAV(I)
150 IF (ABS(X(I)-XT(I)) .GT. DIRIN(I)) GO TO 250
151 220 CONTINUE
152 230 IF (AMIN .LT. APSI) GO TO 350
153 GO TO 325
154C . . . . . . truly new minimum
155 250 LNEWMN = .TRUE.
156 IF (ISW(2) .GE. 1) THEN
157 ISW(2) = 1
158 DCOVAR = MAX(DCOVAR,HALF)
159 ELSE
160 DCOVAR = 1.
161 ENDIF
162 ITAUR = 0
163 NFCNMX = NFCNMX + NPFN - NFCN
164 CSTATU = 'NEW MINIMU'
165 IF (ISW(5) .GE. 0) WRITE (ISYSWR,1030)
166 RETURN
167C . . . return to previous region
168 280 IF (ISW(5) .GT. 0) WRITE (ISYSWR,1020)
169 GO TO 325
170 300 ISW(1) = 1
171 325 DO 330 I= 1, NPAR
172 DIRIN(I) = 0.01*DSAV(I)
173 330 X(I) = XT(I)
174 AMIN = APSI
175 EDM = SIGSAV
176 350 CALL MNINEX(X)
177 IF (ISW(5) .GT. 0) WRITE (ISYSWR,1010)
178 CSTATU= 'UNCHANGED '
179 CALL MNRSET(0)
180 IF (ISW(2) .LT. 2) GO TO 380
181 IF (LOOP .LT. NLOOP .AND. ISW(1) .LT. 1) GO TO 20
182 380 CALL MNPRIN (5,AMIN)
183 ITAUR = 0
184 RETURN
185 1000 FORMAT (54H AN IMPROVEMENT ON THE PREVIOUS MINIMUM HAS BEEN FOUND)
186 1010 FORMAT (51H IMPROVE HAS RETURNED TO REGION OF ORIGINAL MINIMUM)
187 1020 FORMAT (/44H COVARIANCE MATRIX WAS NOT POSITIVE-DEFINITE)
188 1030 FORMAT (/38H IMPROVE HAS FOUND A TRULY NEW MINIMUM/1H ,37(1H*)/)
189 1040 FORMAT (/18H START ATTEMPT NO.,I2, 20H TO FIND NEW MINIMUM)
190 END
Note: See TracBrowser for help on using the repository browser.