EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
lmsimp.F
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file lmsimp.F
1 
2 C **********************************************************************
3 
4  SUBROUTINE lmsimp
5 
6  IMPLICIT NONE
7 
8 C...This is the MINUIT routine SIMPLEX.
9 CC PERFORMS A MINIMIZATION USING THE SIMPLEX METHOD OF NELDER
10 CC AND MEAD (REF. -- COMP. J. 7,308 (1965)).
11  COMMON /lminui/ xkin(4),ukin(4),wkin(4),ain(4),bin(4),
12  &maxfin,relup,relerr,reler2,fcnmax
13  REAL xkin,ukin,wkin,ain,bin,relerr,relup,reler2,fcnmax
14  INTEGER maxfin
15  SAVE /lminui/
16 
17  COMMON /lpflag/ lst3
18  INTEGER lst3
19  SAVE /lpflag/
20 
21  COMMON
22  1/lmmine/ erp(30) ,ern(30)
23  2/lmpari/ x(15) ,xt(15) ,dirin(15) ,maxint ,npar
24  3/lmpare/ u(30) ,werr(30) ,maxext ,nu
25  4/lmlimi/ alim(30) ,blim(30) ,lcode(30) ,lcorsp(30) ,limset
26  5/lmvari/ v(15,15)
27  7/lmfix / ipfix(15),xs(15) ,xts(15) ,dirins(15) ,npfix
28  7/lmfix2/ grds(15) ,g2s(15) ,gsteps(15),aberfs(15)
29  c/lmcasc/ y(16) ,jh ,jl
30  f/lmderi/ gin(30) ,grd(15) ,g2(15) ,gstep(15) ,aberf(15)
31  g/lmsimv/ p(15,16) ,pstar(15),pstst(15) ,pbar(15) ,prho(15)
32  j/lmvart/ vt(15,15)
33  COMMON
34  6/lmunit/ isysrd ,isyswr ,isyspu
35  8/lmtitl/ title(13),date(2) ,isw(7) ,nblock
36  9/lmconv/ epsi ,apsi ,vtest ,nstepq ,nfcn ,nfcnmx
37  a/lmcard/ cword ,cword2 ,cword3 ,word7(7)
38  b/lmmini/ amin ,up ,newmin ,itaur ,sigma,epsmac
39 
40  REAL erp,ern
41  INTEGER maxint,npar
42  REAL x,xt,dirin
43  INTEGER maxext ,nu
44  REAL u,werr
45  INTEGER lcode,lcorsp,limset
46  REAL alim,blim
47  REAL v
48  INTEGER ipfix,npfix
49  REAL xs,xts,dirins
50  REAL grds,g2s,gsteps,aberfs
51  REAL y
52  INTEGER jh,jl
53  REAL gin,grd,g2,gstep,aberf
54  REAL p,pstar,pstst,pbar,prho
55  REAL vt
56 ***** COMMON
57  INTEGER isysrd ,isyswr ,isyspu
58  REAL title,date
59  INTEGER isw,nblock
60  REAL epsi ,apsi ,vtest
61  INTEGER nstepq ,nfcn ,nfcnmx
62  REAL cword ,cword2 ,cword3 ,word7
63  REAL amin ,up ,sigma,epsmac
64  INTEGER newmin ,itaur
65  INTEGER i,iflag,nparp1,npfn,kg,ns,nf,k,ncycl,j,jhold
66  REAL alpha,beta,gamma,rhomin,rhomax,rho1,rho2,wg,ynpp1,absmin,
67  + aming,bestx,f,sig2,pb,ystar,ystst,y1,y2,rho,yrho,ypbar
68 
69  DATA alpha,beta,gamma,rhomin,rhomax / 1.0, 0.5, 2.0, 4.0, 8.0/
70  IF (npar .LE. 0) RETURN
71  npfn=nfcn
72  nparp1=npar+1
73  rho1 = 1.0 + alpha
74  rho2 = rho1 + alpha*gamma
75  wg = 1.0/float(npar)
76  iflag=4
77  IF(lst3.GE.5) WRITE(isyswr,100) epsi
78  DO 2 i= 1, npar
79  IF (isw(2) .GE. 1) dirin(i) = sqrt(v(i,i)*up)
80  IF (abs(dirin(i)) .LT. 1.0e-10*abs(x(i))) dirin(i)=1.0e-8*x(i)
81  IF(itaur.LT. 1) v(i,i) = dirin(i)**2/up
82  2 CONTINUE
83  IF (itaur .LT. 1) isw(2) = 1
84 C** CHOOSE THE INITIAL SIMPLEX USING SINGLE-PARAMETER SEARCHES
85  1 CONTINUE
86  ynpp1 = amin
87  jl = nparp1
88  y(nparp1) = amin
89  absmin = amin
90  DO 10 i= 1, npar
91  aming = amin
92  pbar(i) = x(i)
93  bestx = x(i)
94  kg = 0
95  ns = 0
96  nf = 0
97  4 x(i) = bestx + dirin(i)
98  CALL lminto(x)
99  CALL lsigmx(npar,gin, f, u, 4)
100  nfcn = nfcn + 1
101  IF (f .LE. aming) go to 6
102 C FAILURE
103  IF (kg .EQ. 1) go to 8
104  kg = -1
105  nf = nf + 1
106  dirin(i) = dirin(i) * (-0.4)
107  IF (nf .LT. 3) go to 4
108  ns = 6
109 C SUCCESS
110  6 bestx = x(i)
111  dirin(i) = dirin(i) * 3.0
112  aming = f
113  kg = 1
114  ns = ns + 1
115  IF (ns .LT. 6) go to 4
116 C LOCAL MINIMUM FOUND IN ITH DIRECTION
117  8 y(i) = aming
118  IF (aming .LT. absmin) jl = i
119  IF (aming .LT. absmin) absmin = aming
120  x(i) = bestx
121  DO 9 k= 1, npar
122  9 p(k,i) = x(k)
123  10 CONTINUE
124  jh = nparp1
125  amin=y(jl)
126  CALL lmrazz(ynpp1,pbar)
127  DO 20 i= 1, npar
128  20 x(i) = p(i,jl)
129  CALL lminto(x)
130  DO 30 i=1,npar
131  30 IF(abs(dirin(i)).LE.abs(epsmac*x(i))) dirin(i)=4.*epsmac*x(i)
132  IF (isw(5) .GE. 1) CALL lmprin(0,amin)
133  sigma = sigma * 10.
134  sig2 = sigma
135  ncycl=0
136 C . . . . . START MAIN LOOP
137  50 CONTINUE
138 C...Change in SIMPLX; error redefined for second call to LMSIMP.
139  up=relup*abs(amin)
140  epsi=relerr*up
141  IF (sig2 .LT. epsi .AND. sigma.LT.epsi) go to 76
142  sig2 = sigma
143  IF ((nfcn-npfn) .GT. nfcnmx) go to 78
144 C CALCULATE NEW POINT * BY REFLECTION
145  DO 60 i= 1, npar
146  pb = 0.
147  DO 59 j= 1, nparp1
148  59 pb = pb + wg * p(i,j)
149  pbar(i) = pb - wg * p(i,jh)
150  60 pstar(i)=(1.+alpha)*pbar(i)-alpha*p(i,jh)
151  CALL lminto(pstar)
152  CALL lsigmx(npar,gin,ystar,u,4)
153  nfcn=nfcn+1
154  IF(ystar.GE.amin) go to 70
155 C POINT * BETTER THAN JL, CALCULATE NEW POINT **
156  DO 61 i=1,npar
157  61 pstst(i)=gamma*pstar(i)+(1.-gamma)*pbar(i)
158  CALL lminto(pstst)
159  CALL lsigmx(npar,gin,ystst,u,4)
160  nfcn=nfcn+1
161 C TRY A PARABOLA THROUGH PH, PSTAR, PSTST. MIN = PRHO
162  y1 = (ystar-y(jh)) * rho2
163  y2 = (ystst-y(jh)) * rho1
164  rho = 0.5 * (rho2*y1 -rho1*y2) / (y1 -y2)
165  IF (rho .LT. rhomin) go to 66
166  IF (rho .GT. rhomax) rho = rhomax
167  DO 64 i= 1, npar
168  64 prho(i) = rho*pbar(i) + (1.0-rho)*p(i,jh)
169  CALL lminto(prho)
170  CALL lsigmx(npar,gin,yrho, u,4)
171  nfcn = nfcn + 1
172  IF (yrho .LT. y(jl) .AND. yrho .LT. ystst) go to 65
173  IF (ystst .LT. y(jl)) go to 67
174  IF (yrho .GT. y(jl)) go to 66
175 C ACCEPT MINIMUM POINT OF PARABOLA, PRHO
176  65 CALL lmrazz(yrho,prho)
177  go to 68
178  66 IF (ystst .LT. y(jl)) go to 67
179  CALL lmrazz(ystar,pstar)
180  go to 68
181  67 CALL lmrazz(ystst,pstst)
182  68 ncycl=ncycl+1
183  IF (isw(5) .LT. 2) go to 50
184  IF (isw(5) .GE. 3 .OR. mod(ncycl, 10) .EQ. 0) CALL lmprin(0,amin)
185  go to 50
186 C POINT * IS NOT AS GOOD AS JL
187  70 IF (ystar .GE. y(jh)) go to 73
188  jhold = jh
189  CALL lmrazz(ystar,pstar)
190  IF (jhold .NE. jh) go to 50
191 C CALCULATE NEW POINT **
192  73 DO 74 i=1,npar
193  74 pstst(i)=beta*p(i,jh)+(1.-beta)*pbar(i)
194  CALL lminto(pstst)
195  CALL lsigmx(npar,gin,ystst,u,4)
196  nfcn=nfcn+1
197  IF(ystst.GT.y(jh)) go to 1
198 C POINT ** IS BETTER THAN JH
199  IF (ystst .LT. amin) go to 67
200  CALL lmrazz(ystst,pstst)
201  go to 50
202 C . . . . . . END MAIN LOOP
203  76 IF(lst3.GE.5) WRITE(isyswr,120)
204  go to 80
205  78 IF(lst3.GE.5) WRITE(isyswr,130)
206  isw(1) = 1
207  80 DO 82 i=1,npar
208  pb = 0.
209  DO 81 j=1,nparp1
210  81 pb = pb + wg * p(i,j)
211  82 pbar(i) = pb - wg * p(i,jh)
212  CALL lminto(pbar)
213  CALL lsigmx(npar,gin,ypbar,u,iflag)
214  nfcn=nfcn+1
215  IF (ypbar .LT. amin) CALL lmrazz(ypbar,pbar)
216  CALL lminto(x)
217  IF (nfcnmx+npfn-nfcn .LT. 3*npar) go to 90
218  IF (sigma .GT. 2.0*epsi) go to 1
219  90 CALL lmprin(1-itaur, amin)
220  RETURN
221  100 FORMAT(' START SIMPLEX MINIMIZATION ',8x ,'CON',
222  +.LT.'VERGENCE CRITERION -- ESTIMATED DISTANCE TO MINIMUM (EDM) ',
223  +e10.2 )
224 120 FORMAT(1h ,'SIMPLEX MINIMIZATION HAS CONVERGED')
225 130 FORMAT(1h ,'SIMPLEX TERMINATES WITHOUT CONVERGENCE')
226  END