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

Go to the source code of this file.

Functions/Subroutines

subroutine mnsimp (FCN, FUTIL)
 

Function/Subroutine Documentation

◆ mnsimp()

subroutine mnsimp ( external  FCN,
external  FUTIL 
)

Definition at line 25 of file mnsimp.f.

References a, b, c, d, e, e10, f, g, h, softenpik::half, i, j, m, maxp, mnamin(), mndxdi(), mninex(), mnprin(), mnrazz(), n, o, p, parameter(), up(), x, xs, y, and z.

Referenced by mnexcm(), and mnimpr().

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 minimization using the simplex method of Nelder
40 CC and Mead (ref. -- Comp. J. 7,308 (1965)).
41 CC
42 *
43 * $Id: d506cm.inc,v 1.1.1.1 1996/03/07 14:31:32 mclareni Exp $
44 *
45 * $Log: d506cm.inc,v $
46 * Revision 1.1.1.1 1996/03/07 14:31:32 mclareni
47 * Minuit
48 *
49 *
50 *
51 *
52 * d506cm.inc
53 *
54  parameter(mne=100 , mni=50)
55  parameter(mnihl=mni*(mni+1)/2)
56  CHARACTER*10 cpnam
57  COMMON
58  1/mn7nam/ cpnam(mne)
59  2/mn7ext/ u(mne) ,alim(mne) ,blim(mne)
60  3/mn7err/ erp(mni) ,ern(mni) ,werr(mni) ,globcc(mni)
61  4/mn7inx/ nvarl(mne) ,niofex(mne),nexofi(mni)
62  5/mn7int/ x(mni) ,xt(mni) ,dirin(mni)
63  6/mn7fx2/ xs(mni) ,xts(mni) ,dirins(mni)
64  7/mn7der/ grd(mni) ,g2(mni) ,gstep(mni) ,gin(mne) ,dgrd(mni)
65  8/mn7fx3/ grds(mni) ,g2s(mni) ,gsteps(mni)
66  9/mn7fx1/ ipfix(mni) ,npfix
67  a/mn7var/ vhmat(mnihl)
68  b/mn7vat/ vthmat(mnihl)
69  c/mn7sim/ p(mni,mni+1),pstar(mni),pstst(mni) ,pbar(mni),prho(mni)
70 C
71  parameter(maxdbg=10, maxstk=10, maxcwd=20, maxp=30, maxcpt=101)
72  parameter(zero=0.0, one=1.0, half=0.5)
73  COMMON
74  d/mn7npr/ maxint ,npar ,maxext ,nu
75  e/mn7iou/ isysrd ,isyswr ,isyssa ,npagwd ,npagln ,newpag
76  e/mn7io2/ istkrd(maxstk) ,nstkrd ,istkwr(maxstk) ,nstkwr
77  f/mn7tit/ cfrom ,cstatu ,ctitl ,cword ,cundef ,cvrsn ,covmes
78  g/mn7flg/ isw(7) ,idbg(0:maxdbg) ,nblock ,icomnd
79  h/mn7min/ amin ,up ,edm ,fval3 ,epsi ,apsi ,dcovar
80  i/mn7cnv/ nfcn ,nfcnmx ,nfcnlc ,nfcnfr ,itaur,istrat,nwrmes(2)
81  j/mn7arg/ word7(maxp)
82  k/mn7log/ lwarn ,lrepor ,limset ,lnolim ,lnewmn ,lphead
83  l/mn7cns/ epsmac ,epsma2 ,vlimlo ,vlimhi ,undefi ,bigedm,updflt
84  m/mn7rpt/ xpt(maxcpt) ,ypt(maxcpt)
85  n/mn7cpt/ chpt(maxcpt)
86  o/mn7xcr/ xmidcr ,ymidcr ,xdircr ,ydircr ,ke1cr ,ke2cr
87  CHARACTER ctitl*50, cword*(maxcwd), cundef*10, cfrom*8,
88  + cvrsn*6, covmes(0:3)*22, cstatu*10, chpt*1
89  LOGICAL lwarn, lrepor, limset, lnolim, lnewmn, lphead
90  EXTERNAL fcn,futil
91  dimension y(mni+1)
92  DATA alpha,beta,gamma,rhomin,rhomax / 1.0, 0.5, 2.0, 4.0, 8.0/
93  IF (npar .LE. 0) RETURN
94  IF (amin .EQ. undefi) CALL mnamin(fcn,futil)
95  cfrom = 'SIMPLEX '
96  nfcnfr = nfcn
97  cstatu= 'UNCHANGED '
98  npfn=nfcn
99  nparp1=npar+1
100  nparx = npar
101  rho1 = 1.0 + alpha
102  rho2 = rho1 + alpha*gamma
103  wg = 1.0/float(npar)
104  IF (isw(5) .GE. 0) WRITE(isyswr,100) epsi
105  100 FORMAT(.LT.' START SIMPLEX MINIMIZATION. CONVERGENCE WHEN EDM '
106  +,e10.2 )
107  DO 2 i= 1, npar
108  dirin(i) = werr(i)
109  CALL mndxdi(x(i),i,dxdi)
110  IF (dxdi .NE. zero) dirin(i)=werr(i)/dxdi
111  dmin = epsma2*abs(x(i))
112  IF (dirin(i) .LT. dmin) dirin(i)=dmin
113  2 CONTINUE
114 C** choose the initial simplex using single-parameter searches
115  1 CONTINUE
116  ynpp1 = amin
117  jl = nparp1
118  y(nparp1) = amin
119  absmin = amin
120  DO 10 i= 1, npar
121  aming = amin
122  pbar(i) = x(i)
123  bestx = x(i)
124  kg = 0
125  ns = 0
126  nf = 0
127  4 x(i) = bestx + dirin(i)
128  CALL mninex(x)
129  CALL fcn(nparx,gin, f, u, 4, futil)
130  nfcn = nfcn + 1
131  IF (f .LT. aming) GO TO 6
132 C failure
133  IF (kg .EQ. 1) GO TO 8
134  kg = -1
135  nf = nf + 1
136  dirin(i) = dirin(i) * (-0.4)
137  IF (nf .LT. 3) GO TO 4
138 C stop after three failures
139  bestx = x(i)
140  dirin(i) = dirin(i) * 3.0
141  aming = f
142  GO TO 8
143 C
144 C success
145  6 bestx = x(i)
146  dirin(i) = dirin(i) * 3.0
147  aming = f
148  cstatu= 'PROGRESS '
149  kg = 1
150  ns = ns + 1
151  IF (ns .LT. 6) GO TO 4
152 C
153 C 3 failures or 6 successes or
154 C local minimum found in ith direction
155  8 y(i) = aming
156  IF (aming .LT. absmin) jl = i
157  IF (aming .LT. absmin) absmin = aming
158  x(i) = bestx
159  DO 9 k= 1, npar
160  9 p(k,i) = x(k)
161  10 CONTINUE
162  jh = nparp1
163  amin=y(jl)
164  CALL mnrazz(ynpp1,pbar,y,jh,jl)
165  DO 20 i= 1, npar
166  20 x(i) = p(i,jl)
167  CALL mninex(x)
168  IF (isw(5) .GE. 1) CALL mnprin(5,amin)
169  edm = bigedm
170  sig2 = edm
171  ncycl=0
172 C . . . . . start main loop
173  50 CONTINUE
174  IF (sig2 .LT. epsi .AND. edm.LT.epsi) GO TO 76
175  sig2 = edm
176  IF ((nfcn-npfn) .GT. nfcnmx) GO TO 78
177 C calculate new point * by reflection
178  DO 60 i= 1, npar
179  pb = 0.
180  DO 59 j= 1, nparp1
181  59 pb = pb + wg * p(i,j)
182  pbar(i) = pb - wg * p(i,jh)
183  60 pstar(i)=(1.+alpha)*pbar(i)-alpha*p(i,jh)
184  CALL mninex(pstar)
185  CALL fcn(nparx,gin,ystar,u,4,futil)
186  nfcn=nfcn+1
187  IF(ystar.GE.amin) GO TO 70
188 C point * better than jl, calculate new point **
189  cstatu = 'PROGRESS '
190  DO 61 i=1,npar
191  61 pstst(i)=gamma*pstar(i)+(1.-gamma)*pbar(i)
192  CALL mninex(pstst)
193  CALL fcn(nparx,gin,ystst,u,4,futil)
194  nfcn=nfcn+1
195 C try a parabola through ph, pstar, pstst. min = prho
196  y1 = (ystar-y(jh)) * rho2
197  y2 = (ystst-y(jh)) * rho1
198  rho = 0.5 * (rho2*y1 -rho1*y2) / (y1 -y2)
199  IF (rho .LT. rhomin) GO TO 66
200  IF (rho .GT. rhomax) rho = rhomax
201  DO 64 i= 1, npar
202  64 prho(i) = rho*pbar(i) + (1.0-rho)*p(i,jh)
203  CALL mninex(prho)
204  CALL fcn(nparx,gin,yrho, u,4,futil)
205  nfcn = nfcn + 1
206  IF (yrho .LT. amin) cstatu = 'PROGRESS '
207  IF (yrho .LT. y(jl) .AND. yrho .LT. ystst) GO TO 65
208  IF (ystst .LT. y(jl)) GO TO 67
209  IF (yrho .GT. y(jl)) GO TO 66
210 C accept minimum point of parabola, PRHO
211  65 CALL mnrazz (yrho,prho,y,jh,jl)
212  GO TO 68
213  66 IF (ystst .LT. y(jl)) GO TO 67
214  CALL mnrazz(ystar,pstar,y,jh,jl)
215  GO TO 68
216  67 CALL mnrazz(ystst,pstst,y,jh,jl)
217  68 ncycl=ncycl+1
218  IF (isw(5) .LT. 2) GO TO 50
219  IF (isw(5) .GE. 3 .OR. mod(ncycl, 10) .EQ. 0) CALL mnprin(5,amin)
220  GO TO 50
221 C point * is not as good as jl
222  70 IF (ystar .GE. y(jh)) GO TO 73
223  jhold = jh
224  CALL mnrazz(ystar,pstar,y,jh,jl)
225  IF (jhold .NE. jh) GO TO 50
226 C calculate new point **
227  73 DO 74 i=1,npar
228  74 pstst(i)=beta*p(i,jh)+(1.-beta)*pbar(i)
229  CALL mninex (pstst)
230  CALL fcn(nparx,gin,ystst,u,4,futil)
231  nfcn=nfcn+1
232  IF(ystst.GT.y(jh)) GO TO 1
233 C point ** is better than jh
234  IF (ystst .LT. amin) cstatu = 'PROGRESS '
235  IF (ystst .LT. amin) GO TO 67
236  CALL mnrazz(ystst,pstst,y,jh,jl)
237  GO TO 50
238 C . . . . . . end main loop
239  76 IF (isw(5) .GE. 0) WRITE(isyswr,'(A)')
240  + ' SIMPLEX MINIMIZATION HAS CONVERGED.'
241  isw(4) = 1
242  GO TO 80
243  78 IF (isw(5) .GE. 0) WRITE(isyswr,'(A)')
244  + ' SIMPLEX TERMINATES WITHOUT CONVERGENCE.'
245  cstatu= 'CALL LIMIT'
246  isw(4) = -1
247  isw(1) = 1
248  80 DO 82 i=1,npar
249  pb = 0.
250  DO 81 j=1,nparp1
251  81 pb = pb + wg * p(i,j)
252  82 pbar(i) = pb - wg * p(i,jh)
253  CALL mninex(pbar)
254  CALL fcn(nparx,gin,ypbar,u,4,futil)
255  nfcn=nfcn+1
256  IF (ypbar .LT. amin) CALL mnrazz(ypbar,pbar,y,jh,jl)
257  CALL mninex(x)
258  IF (nfcnmx+npfn-nfcn .LT. 3*npar) GO TO 90
259  IF (edm .GT. 2.0*epsi) GO TO 1
260  90 IF (isw(5) .GE. 0) CALL mnprin(5, amin)
261  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
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
o
struct ob o[NpMax]
Definition: Zprivate.h:34
mndxdi
subroutine mndxdi(PINT, IPAR, DXDI)
Definition: mndxdi.f:10
mninex
subroutine mninex(PINT)
Definition: mninex.f:10
y
latitude latitude this system is used *****************************************************************! type coord sequence union map real y
Definition: Zcoord.h:25
d
dE dx *! Nuc Int sampling table d
Definition: cblkMuInt.h:130
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
nrl_atmos::a
real(4), save a
Definition: cNRLAtmos.f:20
g
dE dx *! Nuc Int sampling table g
Definition: cblkMuInt.h:130
mnrazz
subroutine mnrazz(YNEW, PNEW, Y, JH, JL)
Definition: mnrazz.f:10
e10
*************************block data cblkTracking *************************implicit none data *ExactThick *Freec *RatioToE0 *MagChgDist *TimeStructure *Truncn *Truncx data *IncMuonPolari *KEminObs *ThinSampling *EthinRatio *Generate *LpmEffect *MagPairEmin e10
Definition: cblkTracking.h:9
nrl_atmos::b
real(4), save b
Definition: cNRLAtmos.f:21
maxp
!onst int maxp
Definition: Zprivate.h:3
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: