COSMOS v7.655  COSMOSv7655
(AirShowerMC)
mnline.f
Go to the documentation of this file.
1 *
2 * $Id: mnline.F,v 1.1.1.1 1996/03/07 14:31:30 mclareni Exp $
3 *
4 * $Log: mnline.F,v $
5 * Revision 1.1.1.1 1996/03/07 14:31:30 mclareni
6 * Minuit
7 *
8 *
9  SUBROUTINE mnline(FCN,START,FSTART,STEP,SLOPE,TOLER,FUTIL)
10 *
11 * $Id: d506dp.inc,v 1.1.1.1 1996/03/07 14:31:32 mclareni Exp $
12 *
13 * $Log: d506dp.inc,v $
14 * Revision 1.1.1.1 1996/03/07 14:31:32 mclareni
15 * Minuit
16 *
17 *
18 *
19 *
20 * d506dp.inc
21 *
22 C ************ DOUBLE PRECISION VERSION *************
23  IMPLICIT DOUBLE PRECISION (a-h,o-z)
24 CC Perform a line search from position START
25 CC along direction STEP, where the length of vector STEP
26 CC gives the expected position of minimum.
27 CC FSTART is value of function at START
28 CC SLOPE (if non-zero) is df/dx along STEP at START
29 CC TOLER is initial tolerance of minimum in direction STEP
30 *
31 * $Id: d506cm.inc,v 1.1.1.1 1996/03/07 14:31:32 mclareni Exp $
32 *
33 * $Log: d506cm.inc,v $
34 * Revision 1.1.1.1 1996/03/07 14:31:32 mclareni
35 * Minuit
36 *
37 *
38 *
39 *
40 * d506cm.inc
41 *
42  parameter(mne=100 , mni=50)
43  parameter(mnihl=mni*(mni+1)/2)
44  CHARACTER*10 CPNAM
45  COMMON
46  1/mn7nam/ cpnam(mne)
47  2/mn7ext/ u(mne) ,alim(mne) ,blim(mne)
48  3/mn7err/ erp(mni) ,ern(mni) ,werr(mni) ,globcc(mni)
49  4/mn7inx/ nvarl(mne) ,niofex(mne),nexofi(mni)
50  5/mn7int/ x(mni) ,xt(mni) ,dirin(mni)
51  6/mn7fx2/ xs(mni) ,xts(mni) ,dirins(mni)
52  7/mn7der/ grd(mni) ,g2(mni) ,gstep(mni) ,gin(mne) ,dgrd(mni)
53  8/mn7fx3/ grds(mni) ,g2s(mni) ,gsteps(mni)
54  9/mn7fx1/ ipfix(mni) ,npfix
55  a/mn7var/ vhmat(mnihl)
56  b/mn7vat/ vthmat(mnihl)
57  c/mn7sim/ p(mni,mni+1),pstar(mni),pstst(mni) ,pbar(mni),prho(mni)
58 C
59  parameter(maxdbg=10, maxstk=10, maxcwd=20, maxp=30, maxcpt=101)
60  parameter(zero=0.0, one=1.0, half=0.5)
61  COMMON
62  d/mn7npr/ maxint ,npar ,maxext ,nu
63  e/mn7iou/ isysrd ,isyswr ,isyssa ,npagwd ,npagln ,newpag
64  e/mn7io2/ istkrd(maxstk) ,nstkrd ,istkwr(maxstk) ,nstkwr
65  f/mn7tit/ cfrom ,cstatu ,ctitl ,cword ,cundef ,cvrsn ,covmes
66  g/mn7flg/ isw(7) ,idbg(0:maxdbg) ,nblock ,icomnd
67  h/mn7min/ amin ,up ,edm ,fval3 ,epsi ,apsi ,dcovar
68  i/mn7cnv/ nfcn ,nfcnmx ,nfcnlc ,nfcnfr ,itaur,istrat,nwrmes(2)
69  j/mn7arg/ word7(maxp)
70  k/mn7log/ lwarn ,lrepor ,limset ,lnolim ,lnewmn ,lphead
71  l/mn7cns/ epsmac ,epsma2 ,vlimlo ,vlimhi ,undefi ,bigedm,updflt
72  m/mn7rpt/ xpt(maxcpt) ,ypt(maxcpt)
73  n/mn7cpt/ chpt(maxcpt)
74  o/mn7xcr/ xmidcr ,ymidcr ,xdircr ,ydircr ,ke1cr ,ke2cr
75  CHARACTER CTITL*50, CWORD*(maxcwd), CUNDEF*10, CFROM*8,
76  + cvrsn*6, covmes(0:3)*22, cstatu*10, chpt*1
77  LOGICAL LWARN, LREPOR, LIMSET, LNOLIM, LNEWMN, LPHEAD
78  EXTERNAL fcn,futil
79  dimension start(*), step(*)
80  parameter(maxpt=12)
81  dimension xpq(maxpt),ypq(maxpt)
82  CHARACTER*1 CHPQ(maxpt)
83  dimension xvals(3),fvals(3),coeff(3)
84  CHARACTER*26 CHARAL
85  CHARACTER*60 CMESS
86  parameter(slambg=5.,alpha=2.)
87 C SLAMBG and ALPHA control the maximum individual steps allowed.
88 C The first step is always =1. The max length of second step is SLAMBG.
89 C The max size of subsequent steps is the maximum previous successful
90 C step multiplied by ALPHA + the size of most recent successful step,
91 C but cannot be smaller than SLAMBG.
92  LOGICAL LDEBUG
93  DATA charal / 'ABCDEFGHIJKLMNOPQRSTUVWXYZ' /
94  ldebug = (idbg(1).GE.1)
95 C starting values for overall limits on total step SLAM
96  overal = 1000.
97  undral = -100.
98 C debug check if start is ok
99  IF (ldebug) THEN
100  CALL mninex(start)
101  CALL fcn(nparx,gin,f1,u,4,futil)
102  nfcn=nfcn+1
103  IF (f1 .NE. fstart) THEN
104  WRITE (isyswr,'(A/2E14.5/2X,10F10.5)')
105  + ' MNLINE start point not consistent, F values, parameters=',
106  + (x(kk),kk=1,npar)
107  ENDIF
108  ENDIF
109 C . set up linear search along STEP
110  fvmin = fstart
111  xvmin = zero
112  nxypt = 1
113  chpq(1) = charal(1:1)
114  xpq(1) = 0.
115  ypq(1) = fstart
116 C SLAMIN = smallest possible value of ABS(SLAM)
117  slamin = zero
118  DO 20 i= 1, npar
119  IF (step(i) .EQ. zero) GO TO 20
120  ratio = abs(start(i)/step(i))
121  IF (slamin .EQ. zero) slamin = ratio
122  IF (ratio .LT. slamin) slamin = ratio
123  20 x(i) = start(i) + step(i)
124  IF (slamin .EQ. zero) slamin = epsmac
125  slamin = slamin*epsma2
126  nparx = npar
127 C
128  CALL mninex(x)
129  CALL fcn(nparx,gin,f1,u,4,futil)
130  nfcn=nfcn+1
131  nxypt = nxypt + 1
132  chpq(nxypt) = charal(nxypt:nxypt)
133  xpq(nxypt) = 1.
134  ypq(nxypt) = f1
135  IF (f1 .LT. fstart) THEN
136  fvmin = f1
137  xvmin = 1.0
138  ENDIF
139 C . quadr interp using slope GDEL and two points
140  slam = 1.
141  toler8 = toler
142  slamax = slambg
143  flast = f1
144 C can iterate on two-points (cut) if no imprvmnt
145  25 CONTINUE
146  denom = 2.0*(flast-fstart-slope*slam)/slam**2
147 C IF (DENOM .EQ. ZERO) DENOM = -0.1*SLOPE
148  slam = 1.
149  IF (denom .NE. zero) slam = -slope/denom
150  IF (slam .LT. zero) slam = slamax
151  IF (slam .GT. slamax) slam = slamax
152  IF (slam .LT. toler8) slam = toler8
153  IF (slam .LT. slamin) GO TO 80
154  IF (abs(slam-1.0).LT.toler8 .AND. f1.LT.fstart) GO TO 70
155  IF (abs(slam-1.0).LT.toler8) slam = 1.0+toler8
156  IF (nxypt .GE. maxpt) GO TO 65
157  DO 30 i= 1, npar
158  30 x(i) = start(i) + slam*step(i)
159  CALL mninex(x)
160  CALL fcn(npar,gin,f2,u,4,futil)
161  nfcn = nfcn + 1
162  nxypt = nxypt + 1
163  chpq(nxypt) = charal(nxypt:nxypt)
164  xpq(nxypt) = slam
165  ypq(nxypt) = f2
166  IF (f2 .LT. fvmin) THEN
167  fvmin = f2
168  xvmin = slam
169  ENDIF
170  IF (fstart .EQ. fvmin) THEN
171  flast = f2
172  toler8 = toler*slam
173  overal = slam-toler8
174  slamax = overal
175  GO TO 25
176  ENDIF
177 C . quadr interp using 3 points
178  xvals(1) = xpq(1)
179  fvals(1) = ypq(1)
180  xvals(2) = xpq(nxypt-1)
181  fvals(2) = ypq(nxypt-1)
182  xvals(3) = xpq(nxypt)
183  fvals(3) = ypq(nxypt)
184 C begin iteration, calculate desired step
185  50 CONTINUE
186  slamax = max(slamax,alpha*abs(xvmin))
187  CALL mnpfit(xvals,fvals,3,coeff,sdev)
188  IF (coeff(3) .LE. zero) THEN
189  slopem = 2.0*coeff(3)*xvmin + coeff(2)
190  IF (slopem .LE. zero) THEN
191  slam = xvmin + slamax
192  ELSE
193  slam = xvmin - slamax
194  ENDIF
195  ELSE
196  slam = -coeff(2)/(2.0*coeff(3))
197  IF (slam .GT. xvmin+slamax) slam = xvmin+slamax
198  IF (slam .LT. xvmin-slamax) slam = xvmin-slamax
199  ENDIF
200  IF (slam .GT. zero) THEN
201  IF (slam .GT. overal) slam = overal
202  ELSE
203  IF (slam .LT. undral) slam = undral
204  ENDIF
205 C come here if step was cut below
206  52 CONTINUE
207  toler9 = max(toler8,abs(toler8*slam))
208  DO 55 ipt= 1, 3
209  IF (abs(slam-xvals(ipt)) .LT. toler9) GO TO 70
210  55 CONTINUE
211 C take the step
212  IF (nxypt .GE. maxpt) GO TO 65
213  DO 60 i= 1, npar
214  60 x(i) = start(i)+slam*step(i)
215  CALL mninex(x)
216  CALL fcn(nparx,gin,f3,u,4,futil)
217  nfcn = nfcn + 1
218  nxypt = nxypt + 1
219  chpq(nxypt) = charal(nxypt:nxypt)
220  xpq(nxypt) = slam
221  ypq(nxypt) = f3
222 C find worst previous point out of three
223  fvmax = fvals(1)
224  nvmax = 1
225  IF (fvals(2) .GT. fvmax) THEN
226  fvmax = fvals(2)
227  nvmax = 2
228  ENDIF
229  IF (fvals(3) .GT. fvmax) THEN
230  fvmax = fvals(3)
231  nvmax = 3
232  ENDIF
233 C if latest point worse than all three previous, cut step
234  IF (f3 .GE. fvmax) THEN
235  IF (nxypt .GE. maxpt) GO TO 65
236  IF (slam .GT. xvmin) overal = min(overal,slam-toler8)
237  IF (slam .LT. xvmin) undral = max(undral,slam+toler8)
238  slam = 0.5*(slam+xvmin)
239  GO TO 52
240  ENDIF
241 C prepare another iteration, replace worst previous point
242  xvals(nvmax) = slam
243  fvals(nvmax) = f3
244  IF (f3 .LT. fvmin) THEN
245  fvmin = f3
246  xvmin = slam
247  ELSE
248  IF (slam .GT. xvmin) overal = min(overal,slam-toler8)
249  IF (slam .LT. xvmin) undral = max(undral,slam+toler8)
250  ENDIF
251  IF (nxypt .LT. maxpt) GO TO 50
252 C . . end of iteration . . .
253 C stop because too many iterations
254  65 cmess = ' LINE SEARCH HAS EXHAUSTED THE LIMIT OF FUNCTION CALLS '
255  IF (ldebug) THEN
256  WRITE (isyswr,'(A/(2X,6G12.4))') ' MNLINE DEBUG: steps=',
257  + (step(kk),kk=1,npar)
258  ENDIF
259  GO TO 100
260 C stop because within tolerance
261  70 CONTINUE
262  cmess = ' LINE SEARCH HAS ATTAINED TOLERANCE '
263  GO TO 100
264  80 CONTINUE
265  cmess = ' STEP SIZE AT ARITHMETICALLY ALLOWED MINIMUM'
266  100 CONTINUE
267  amin = fvmin
268  DO 120 i= 1, npar
269  dirin(i) = step(i)*xvmin
270  120 x(i) = start(i) + dirin(i)
271  CALL mninex(x)
272  IF (xvmin .LT. 0.) CALL mnwarn('D','MNLINE',
273  + ' LINE MINIMUM IN BACKWARDS DIRECTION')
274  IF (fvmin .EQ. fstart) CALL mnwarn('D','MNLINE',
275  + ' LINE SEARCH FINDS NO IMPROVEMENT ')
276  IF (ldebug) THEN
277  WRITE (isyswr,'('' AFTER'',I3,'' POINTS,'',A)') nxypt,cmess
278  CALL mnplot(xpq,ypq,chpq,nxypt,isyswr,npagwd,npagln)
279  ENDIF
280  RETURN
281  END
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
mnline
subroutine mnline(FCN, START, FSTART, STEP, SLOPE, TOLER, FUTIL)
Definition: mnline.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
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
mninex
subroutine mninex(PINT)
Definition: mninex.f:10
d
dE dx *! Nuc Int sampling table d
Definition: cblkMuInt.h:130
a
nodes a
Definition: cstdatmos0.fNew2.h:93
mnpfit
subroutine mnpfit(PARX2P, PARY2P, NPAR2P, COEF2P, SDEV2P)
Definition: mnpfit.f:10
b
dE dx *! Nuc Int sampling table b
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
h
dE dx *! Nuc Int sampling table h
Definition: cblkMuInt.h:130
g
dE dx *! Nuc Int sampling table g
Definition: cblkMuInt.h:130
maxp
!onst int maxp
Definition: Zprivate.h:3
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
mnplot
subroutine mnplot(XPT, YPT, CHPT, NXYPT, NUNIT, NPAGWD, NPAGLN)
Definition: mnplot.f:10
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