COSMOS v7.655  COSMOSv7655
(AirShowerMC)
Functions/Subroutines
mnmigr.f File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine mnmigr (FCN, FUTIL)
 

Function/Subroutine Documentation

◆ mnmigr()

subroutine mnmigr ( external  FCN,
external  FUTIL 
)

Definition at line 25 of file mnmigr.f.

References a, b, c, d, e, f, g, h, softenpik::half, i, j, m, maxp, mnamin(), mnderi(), mnhess(), mninex(), mnline(), mnmatu(), mnprin(), mnpsdf(), mnwarn(), mnwerr(), n, o, p, parameter(), up(), vg, x, xs, and z.

Referenced by mncuve(), mneval(), and mnexcm().

25 *
26 * $Id: d506dp.inc,v 1.1.1.1 1996/03/07 14:31:32 mclareni Exp $
27 *
28 * $Log: d506dp.inc,v $
29 * Revision 1.1.1.1 1996/03/07 14:31:32 mclareni
30 * Minuit
31 *
32 *
33 *
34 *
35 * d506dp.inc
36 *
37 C ************ DOUBLE PRECISION VERSION *************
38  IMPLICIT DOUBLE PRECISION (a-h,o-z)
39 CC Performs a local function minimization using basically the
40 CC method of Davidon-Fletcher-Powell as modified by Fletcher
41 CC ref. -- Fletcher, Comp.J. 13,317 (1970) "switching method"
42 CC
43 *
44 * $Id: d506cm.inc,v 1.1.1.1 1996/03/07 14:31:32 mclareni Exp $
45 *
46 * $Log: d506cm.inc,v $
47 * Revision 1.1.1.1 1996/03/07 14:31:32 mclareni
48 * Minuit
49 *
50 *
51 *
52 *
53 * d506cm.inc
54 *
55  parameter(mne=100 , mni=50)
56  parameter(mnihl=mni*(mni+1)/2)
57  CHARACTER*10 cpnam
58  COMMON
59  1/mn7nam/ cpnam(mne)
60  2/mn7ext/ u(mne) ,alim(mne) ,blim(mne)
61  3/mn7err/ erp(mni) ,ern(mni) ,werr(mni) ,globcc(mni)
62  4/mn7inx/ nvarl(mne) ,niofex(mne),nexofi(mni)
63  5/mn7int/ x(mni) ,xt(mni) ,dirin(mni)
64  6/mn7fx2/ xs(mni) ,xts(mni) ,dirins(mni)
65  7/mn7der/ grd(mni) ,g2(mni) ,gstep(mni) ,gin(mne) ,dgrd(mni)
66  8/mn7fx3/ grds(mni) ,g2s(mni) ,gsteps(mni)
67  9/mn7fx1/ ipfix(mni) ,npfix
68  a/mn7var/ vhmat(mnihl)
69  b/mn7vat/ vthmat(mnihl)
70  c/mn7sim/ p(mni,mni+1),pstar(mni),pstst(mni) ,pbar(mni),prho(mni)
71 C
72  parameter(maxdbg=10, maxstk=10, maxcwd=20, maxp=30, maxcpt=101)
73  parameter(zero=0.0, one=1.0, half=0.5)
74  COMMON
75  d/mn7npr/ maxint ,npar ,maxext ,nu
76  e/mn7iou/ isysrd ,isyswr ,isyssa ,npagwd ,npagln ,newpag
77  e/mn7io2/ istkrd(maxstk) ,nstkrd ,istkwr(maxstk) ,nstkwr
78  f/mn7tit/ cfrom ,cstatu ,ctitl ,cword ,cundef ,cvrsn ,covmes
79  g/mn7flg/ isw(7) ,idbg(0:maxdbg) ,nblock ,icomnd
80  h/mn7min/ amin ,up ,edm ,fval3 ,epsi ,apsi ,dcovar
81  i/mn7cnv/ nfcn ,nfcnmx ,nfcnlc ,nfcnfr ,itaur,istrat,nwrmes(2)
82  j/mn7arg/ word7(maxp)
83  k/mn7log/ lwarn ,lrepor ,limset ,lnolim ,lnewmn ,lphead
84  l/mn7cns/ epsmac ,epsma2 ,vlimlo ,vlimhi ,undefi ,bigedm,updflt
85  m/mn7rpt/ xpt(maxcpt) ,ypt(maxcpt)
86  n/mn7cpt/ chpt(maxcpt)
87  o/mn7xcr/ xmidcr ,ymidcr ,xdircr ,ydircr ,ke1cr ,ke2cr
88  CHARACTER ctitl*50, cword*(maxcwd), cundef*10, cfrom*8,
89  + cvrsn*6, covmes(0:3)*22, cstatu*10, chpt*1
90  LOGICAL lwarn, lrepor, limset, lnolim, lnewmn, lphead
91  EXTERNAL fcn,futil
92  dimension gs(mni), step(mni), xxs(mni), flnu(mni), vg(mni)
93  LOGICAL ldebug
94  parameter(toler=0.05)
95  IF (npar .LE. 0) RETURN
96  IF (amin .EQ. undefi) CALL mnamin(fcn,futil)
97  ldebug = (idbg(4) .GE. 1)
98  cfrom = 'MIGRAD '
99  nfcnfr = nfcn
100  nfcnmg = nfcn
101  cstatu= 'INITIATE '
102  iswtr = isw(5) - 2*itaur
103  npfn = nfcn
104  nparx = npar
105  lenv = npar*(npar+1)/2
106  nrstrt = 0
107  npsdf = 0
108  lined2 = 0
109  isw(4) = -1
110  rhotol = 1.0e-3*apsi
111  IF (iswtr .GE. 1) WRITE (isyswr,470) istrat,rhotol
112  470 FORMAT (' START MIGRAD MINIMIZATION. STRATEGY',i2,
113  +.LT.'. CONVERGENCE WHEN EDM ',e9.2)
114 C initialization strategy
115  IF (istrat.LT.2 .OR. isw(2).GE.3) GO TO 2
116 C come (back) here to restart completely
117  1 CONTINUE
118  IF (nrstrt .GT. istrat) THEN
119  cstatu= 'FAILED '
120  isw(4) = -1
121  GO TO 230
122  ENDIF
123 C . get full covariance and gradient
124  CALL mnhess(fcn,futil)
125  CALL mnwerr
126  npsdf = 0
127  IF (isw(2) .GE. 1) GO TO 10
128 C . get gradient at start point
129  2 CONTINUE
130  CALL mninex(x)
131  IF (isw(3) .EQ. 1) THEN
132  CALL fcn(nparx,gin,fzero,u,2,futil)
133  nfcn = nfcn + 1
134  ENDIF
135  CALL mnderi(fcn,futil)
136  IF (isw(2) .GE. 1) GO TO 10
137 C sometimes start with diagonal matrix
138  DO 3 i= 1, npar
139  xxs(i) = x(i)
140  step(i) = zero
141  3 CONTINUE
142 C do line search if second derivative negative
143  lined2 = lined2 + 1
144  IF (lined2 .LT. (istrat+1)*npar) THEN
145  DO 5 i= 1, npar
146  IF (g2(i) .GT. zero) GO TO 5
147  step(i) = -sign(gstep(i),grd(i))
148  gdel = step(i)*grd(i)
149  fs = amin
150  CALL mnline(fcn,xxs,fs,step,gdel,toler,futil)
151  CALL mnwarn('D','MNMIGR','Negative G2 line search')
152  iext = nexofi(i)
153  IF (ldebug) WRITE (isyswr,'(A,I3,2G13.3)')
154  + ' Negative G2 line search, param ',iext,fs,amin
155  GO TO 2
156  5 CONTINUE
157  ENDIF
158 C make diagonal error matrix
159  DO 8 i=1,npar
160  ndex = i*(i-1)/2
161  DO 7 j=1,i-1
162  ndex = ndex + 1
163  7 vhmat(ndex) = 0.
164  ndex = ndex + 1
165  IF (g2(i) .LE. zero) g2(i) = 1.
166  vhmat(ndex) = 2./g2(i)
167  8 CONTINUE
168  dcovar = 1.
169  IF (ldebug) WRITE (isyswr,'(A,A/(1X,10G10.2))') ' DEBUG MNMIGR,',
170  + ' STARTING MATRIX DIAGONAL, VHMAT=', (vhmat(kk),kk=1,lenv)
171 C ready to start first iteration
172  10 CONTINUE
173  nrstrt = nrstrt + 1
174  IF (nrstrt .GT. istrat+1) THEN
175  cstatu= 'FAILED '
176  GO TO 230
177  ENDIF
178  fs = amin
179 C . . . get EDM and set up loop
180  edm = 0.
181  DO 18 i= 1, npar
182  gs(i) = grd(i)
183  xxs(i) = x(i)
184  ndex = i*(i-1)/2
185  DO 17 j= 1, i-1
186  ndex = ndex + 1
187  17 edm = edm + gs(i)*vhmat(ndex)*gs(j)
188  ndex = ndex + 1
189  18 edm = edm + 0.5 * gs(i)**2 *vhmat(ndex)
190  edm = edm * 0.5 * (1.0+3.0*dcovar)
191  IF (edm .LT. zero) THEN
192  CALL mnwarn('W','MIGRAD','STARTING MATRIX NOT POS-DEFINITE.')
193  isw(2) = 0
194  dcovar = 1.
195  GO TO 2
196  ENDIF
197  IF (isw(2) .EQ. 0) edm=bigedm
198  iter = 0
199  CALL mninex(x)
200  CALL mnwerr
201  IF (iswtr .GE. 1) CALL mnprin(3,amin)
202  IF (iswtr .GE. 2) CALL mnmatu(0)
203 C . . . . . start main loop
204  24 CONTINUE
205  IF (nfcn-npfn .GE. nfcnmx) GO TO 190
206  gdel = 0.
207  gssq = 0.
208  DO 30 i=1,npar
209  ri = 0.
210  gssq = gssq + gs(i)**2
211  DO 25 j=1,npar
212  m = max(i,j)
213  n = min(i,j)
214  ndex = m*(m-1)/2 + n
215  25 ri = ri + vhmat(ndex) *gs(j)
216  step(i) = -0.5*ri
217  30 gdel = gdel + step(i)*gs(i)
218  IF (gssq .EQ. zero) THEN
219  CALL mnwarn('D','MIGRAD',
220  + ' FIRST DERIVATIVES OF FCN ARE ALL ZERO')
221  GO TO 300
222  ENDIF
223 C if gdel positive, V not posdef
224  IF (gdel .GE. zero) THEN
225  CALL mnwarn('D','MIGRAD',' NEWTON STEP NOT DESCENT.')
226  IF (npsdf .EQ. 1) GO TO 1
227  CALL mnpsdf
228  npsdf = 1
229  GO TO 24
230  ENDIF
231 C . . . . do line search
232  CALL mnline(fcn,xxs,fs,step,gdel,toler,futil)
233  IF (amin .EQ. fs) GO TO 200
234  cfrom = 'MIGRAD '
235  nfcnfr = nfcnmg
236  cstatu= 'PROGRESS '
237 C . get gradient at new point
238  CALL mninex(x)
239  IF (isw(3) .EQ. 1) THEN
240  CALL fcn(nparx,gin,fzero,u,2,futil)
241  nfcn = nfcn + 1
242  ENDIF
243  CALL mnderi(fcn,futil)
244 C . calculate new EDM
245  npsdf = 0
246  81 edm = 0.
247  gvg = 0.
248  delgam = 0.
249  gdgssq = 0.
250  DO 100 i= 1, npar
251  ri = 0.
252  vgi = 0.
253  DO 90 j= 1, npar
254  m = max(i,j)
255  n = min(i,j)
256  ndex = m*(m-1)/2 + n
257  vgi = vgi + vhmat(ndex)*(grd(j)-gs(j))
258  90 ri = ri + vhmat(ndex)* grd(j)
259  vg(i) = vgi*0.5
260  gami = grd(i) - gs(i)
261  gdgssq = gdgssq + gami**2
262  gvg = gvg + gami*vg(i)
263  delgam = delgam + dirin(i)*gami
264  100 edm = edm + grd(i)*ri*0.5
265  edm = edm * 0.5 * (1.0 + 3.0*dcovar)
266 C . if EDM negative, not positive-definite
267  IF (edm .LT. zero .OR. gvg .LE. zero) THEN
268  CALL mnwarn('D','MIGRAD','NOT POS-DEF. EDM OR GVG NEGATIVE.')
269  cstatu = 'NOT POSDEF'
270  IF (npsdf .EQ. 1) GO TO 230
271  CALL mnpsdf
272  npsdf = 1
273  GO TO 81
274  ENDIF
275 C print information about this iteration
276  iter = iter + 1
277  IF (iswtr.GE.3 .OR. (iswtr.EQ.2.AND.mod(iter,10).EQ.1)) THEN
278  CALL mnwerr
279  CALL mnprin(3,amin)
280  ENDIF
281  IF (gdgssq .EQ. zero) CALL mnwarn('D','MIGRAD',
282  + 'NO CHANGE IN FIRST DERIVATIVES OVER LAST STEP')
283  IF (delgam .LT. zero) CALL mnwarn('D','MIGRAD',
284  + 'FIRST DERIVATIVES INCREASING ALONG SEARCH LINE')
285 C . update covariance matrix
286  cstatu = 'IMPROVEMNT'
287  IF (ldebug) WRITE (isyswr,'(A,(1X,10G10.3))') ' VHMAT 1 =',
288  + (vhmat(kk),kk=1,10)
289  dsum = 0.
290  vsum = 0.
291  DO 120 i=1, npar
292  DO 120 j=1, i
293  d = dirin(i)*dirin(j)/delgam - vg(i)*vg(j)/gvg
294  dsum = dsum + abs(d)
295  ndex = i*(i-1)/2 + j
296  vhmat(ndex) = vhmat(ndex) + 2.0*d
297  vsum = vsum + abs(vhmat(ndex))
298  120 CONTINUE
299 C smooth local fluctuations by averaging DCOVAR
300  dcovar = 0.5*(dcovar + dsum/vsum)
301  IF (iswtr.GE.3 .OR. ldebug) WRITE (isyswr,'(A,F5.1,A)')
302  + ' RELATIVE CHANGE IN COV. MATRIX=',dcovar*100.,'%'
303  IF (ldebug) WRITE (isyswr,'(A,(1X,10G10.3))') ' VHMAT 2 =',
304  + (vhmat(kk),kk=1,10)
305  IF (delgam .LE. gvg) GO TO 135
306  DO 125 i= 1, npar
307  125 flnu(i) = dirin(i)/delgam - vg(i)/gvg
308  DO 130 i= 1, npar
309  DO 130 j= 1, i
310  ndex = i*(i-1)/2 + j
311  130 vhmat(ndex) = vhmat(ndex) + 2.0*gvg*flnu(i)*flnu(j)
312  135 CONTINUE
313 C and see if converged
314  IF (edm .LT. 0.1*rhotol) GO TO 300
315 C if not, prepare next iteration
316  DO 140 i= 1, npar
317  xxs(i) = x(i)
318  gs(i) = grd(i)
319  140 CONTINUE
320  fs = amin
321  IF (isw(2) .EQ. 0 .AND. dcovar.LT. 0.5 ) isw(2) = 1
322  IF (isw(2) .EQ. 3 .AND. dcovar.GT. 0.1 ) isw(2) = 1
323  IF (isw(2) .EQ. 1 .AND. dcovar.LT. 0.05) isw(2) = 3
324  GO TO 24
325 C . . . . . end main loop
326 C . . call limit in MNMIGR
327  190 isw(1) = 1
328  IF (isw(5) .GE. 0)
329  + WRITE (isyswr,'(A)') ' CALL LIMIT EXCEEDED IN MIGRAD.'
330  cstatu = 'CALL LIMIT'
331  GO TO 230
332 C . . fails to improve . .
333  200 IF (iswtr .GE. 1) WRITE (isyswr,'(A)')
334  + ' MIGRAD FAILS TO FIND IMPROVEMENT'
335  DO 210 i= 1, npar
336  210 x(i) = xxs(i)
337  IF (edm .LT. rhotol) GO TO 300
338  IF (edm .LT. abs(epsma2*amin)) THEN
339  IF (iswtr .GE. 0) WRITE (isyswr, '(A)')
340  + ' MACHINE ACCURACY LIMITS FURTHER IMPROVEMENT.'
341  GO TO 300
342  ENDIF
343  IF (istrat .LT. 1) THEN
344  IF (isw(5) .GE. 0) WRITE (isyswr, '(A)')
345  + ' MIGRAD FAILS WITH STRATEGY=0. WILL TRY WITH STRATEGY=1.'
346  istrat = 1
347  ENDIF
348  GO TO 1
349 C . . fails to converge
350  230 IF (iswtr .GE. 0) WRITE (isyswr,'(A)')
351  + ' MIGRAD TERMINATED WITHOUT CONVERGENCE.'
352  IF (isw(2) .EQ. 3) isw(2) = 1
353  isw(4) = -1
354  GO TO 400
355 C . . apparent convergence
356  300 IF (iswtr .GE. 0) WRITE(isyswr,'(/A)')
357  + ' MIGRAD MINIMIZATION HAS CONVERGED.'
358  IF (itaur .EQ. 0) THEN
359  IF (istrat .GE. 2 .OR. (istrat.EQ.1.AND.isw(2).LT.3)) THEN
360  IF (isw(5) .GE. 0) WRITE (isyswr, '(/A)')
361  + ' MIGRAD WILL VERIFY CONVERGENCE AND ERROR MATRIX.'
362  CALL mnhess(fcn,futil)
363  CALL mnwerr
364  npsdf = 0
365  IF (edm .GT. rhotol) GO TO 10
366  ENDIF
367  ENDIF
368  cstatu='CONVERGED '
369  isw(4) = 1
370 C come here in any case
371  400 CONTINUE
372  cfrom = 'MIGRAD '
373  nfcnfr = nfcnmg
374  CALL mninex(x)
375  CALL mnwerr
376  IF (iswtr .GE. 0) CALL mnprin (3,amin)
377  IF (iswtr .GE. 1) CALL mnmatu(1)
378  RETURN
parameter
integer npitbl real *nx parameter(n=101, npitbl=46, nx=n-1) real *8 uconst
z
nodes z
Definition: cstdatmos0.fNew2.h:93
up
block data include Zlatfit h c fitting region data data data data data d0 data data d0 data data h g *is for param c g data up(2, 1)/7.0d0/
e
dE dx *! Nuc Int sampling table e
Definition: cblkMuInt.h:130
i
nodes i
Definition: cstdatmos0.fNew2.h:93
mnwerr
subroutine mnwerr
Definition: mnwerr.f:10
mnline
subroutine mnline(FCN, START, FSTART, STEP, SLOPE, TOLER, FUTIL)
Definition: mnline.f:10
nrl_atmos::h
real(4), dimension(:), allocatable, save h
Definition: cNRLAtmos.f:28
mnprin
subroutine mnprin(INKODE, FVAL)
Definition: mnprin.f:10
xs
! constants thru Cosmos real ! if multiplied to deg radian Torad ! light velocity m sec ! infinty ! kg m2 *Togpcm2 g cm2 ! g cm2 *Tokgpm2 kg m2 ! cm *Tom m ! m *Tocm cm ! g cm3 *Tokgpm3 kg m3 ! kg m3 *Togpcm3 g cm3 ! sec *Tonsec nsec ! Tesla m ! Avogadro *A2deninv ! mfp *n * xs
Definition: Zglobalc.h:18
mnamin
subroutine mnamin(FCN, FUTIL)
Definition: mnamin.f:10
p
********************block data cblkHeavy ********************integer j data *HeavyG2symbol p
Definition: cblkHeavy.h:7
j
********************block data cblkHeavy ********************integer j data *HeavyG2symbol *data *HeavyG2code kiron data j
Definition: cblkHeavy.h:36
mnhess
subroutine mnhess(FCN, FUTIL)
Definition: mnhess.f:10
o
struct ob o[NpMax]
Definition: Zprivate.h:34
mninex
subroutine mninex(PINT)
Definition: mninex.f:10
mnpsdf
subroutine mnpsdf
Definition: mnpsdf.f:25
d
dE dx *! Nuc Int sampling table d
Definition: cblkMuInt.h:130
mnmatu
subroutine mnmatu(KODE)
Definition: mnmatu.f:10
m
block data include Zlatfit h c fitting region data data data data data d0 data data d0 data data m
Definition: ZlatfitBD.h:35
vg
! Zkcele h ! unit here is aunit parameter pi real vg
Definition: Zkcele.h:6
nrl_atmos::a
real(4), save a
Definition: cNRLAtmos.f:20
g
dE dx *! Nuc Int sampling table g
Definition: cblkMuInt.h:130
nrl_atmos::b
real(4), save b
Definition: cNRLAtmos.f:21
maxp
!onst int maxp
Definition: Zprivate.h:3
mnderi
subroutine mnderi(FCN, FUTIL)
Definition: mnderi.f:25
mnwarn
subroutine mnwarn(COPT, CORG, CMES)
Definition: mnwarn.f:10
n
integer n
Definition: Zcinippxc.h:1
softenpik::half
integer, parameter half
Definition: csoftenPiK.f:108
x
! structure defining a particle at production ! Basic idea of what is to be contained in ! the particle structue is that dynamical ones should be included those derivable from the particle code ! is not included ******************************************************type fmom momentum sequence union map real e endmap map real * x
Definition: Zptcl.h:21
f
dE dx *! Nuc Int sampling table f
Definition: cblkMuInt.h:130
c
dE dx *! Nuc Int sampling table c
Definition: cblkMuInt.h:130
Here is the call graph for this function:
Here is the caller graph for this function: