Millepede-II  V04-03-09
isuper.f
Go to the documentation of this file.
1 *CMZ : 11/12/2006 16.01.29 by Claus Kleinwort
2 *-- Author : Claus Kleinwort 14/04/99
3 
4 * A Least Squares Method for Detector Alignment - Fortran code
5 *
6 * DALIGN short test program for detector aligment
7 *
8 * INITGL initialization
9 * PARGLO optional: initialize parameters with nonzero values
10 * PARSIG optional: define sigma for single parameter
11 * INITUN optional: unit for iterations
12 * INITCS optional: constraints
13 * EQULOC equations for local fit
14 * FITLOC local parameter fit
15 * FITGLO final global parameter fit
16 *
17 *
18 * Special matrix subprograms (all in Double Precision):
19 * SPMINV matrix inversion + solution
20 * SPAVAT special double matrix product
21 * SPAX special matrix product
22 *
23 * CHINDL limit of chi**2/nd
24 
25  SUBROUTINE dalign ! test program for detector alignment
26  REAL DERGB(8),DERLC(2),PAR(8)
27  parameter(nplan=10)
28  REAL Z(nplan),SIGMA(nplan),SL(nplan),BIAS(nplan),FONE(8),FROT(8)
29  DATA sigma/0.0020,0.0020,8*0.0300/
30  DATA sl/5.0,9.0,20.0,25.0,30.0,35.0,40.0,45.0,50.0,55.0/
31  DATA bias/0.0,0.0,+0.1,+0.05,-0.10,+0.05,-0.1,+0.1,-0.05,-0.05/
32  ncases=10000
33  CALL initgl(8,2,3,0) ! define dimension parameters
34 C CALL PARSIG(3,0.1)
35 C CALL PARSIG(4,0.0001)
36 C CALL PARGLO(BIAS(3))
37 * CALL INITUN(11,4.0)
38 C CALL PARSIG(8,0.0)
39  DO i=1,8
40  fone(i)=1.0
41  frot(i)=0.0
42  frot(i)=1.0/(sl(i+2)-37.5)
43  END DO
44 C CALL CONSTF(FONE,0.0)
45 C CALL CONSTF(FROT,0.0)
46  CALL zerloc(dergb,derlc)
47  DO nc=1,ncases
48 * generate straight line parameters
49  slope= 2.0*urand()-1.0 ! SLOPE = -1 ... +1
50  znull=20.0*urand()-10.0 ! ZNULL = -10 ... +10
51  DO i=1,nplan ! planes
52  z(i)=znull+slope*sl(i) + bias(i)+sigma(i)*unorm()
53  END DO
54  DO i=1,nplan
55 * calibrate Z(I) = ZNULL + SLOPE * SL(I)
56  derlc(1)= 1.0
57  derlc(2)=sl(i)
58  IF(i.GT.2) dergb(i-2)=1.0
59  CALL equloc(dergb,derlc,z(i),sigma(i)) ! single meas.
60  END DO
61  CALL fitloc
62  END DO
63  CALL fitglo(par)
64  CALL prtglo(66) ! print on unit 66 (-> file fort.66)
65  END
66  FUNCTION urand() ! return random number U(0,1)
67 * (simple generator, showing principle)
68  parameter(ia=205,ic=29573,im=139968)
69  DATA last/4711/
70  last=mod(ia*last+ic,im)
71  IF(last.EQ.0) last=mod(ia*last+ic,im)
72  urand=float(last)/float(im)
73  END
74  FUNCTION unorm() ! return random number N(0,1)
75 * (simple generator, showing principle)
76  parameter(ia=205,ic=29573,im=139968)
77  DATA last/4711/
78  unorm=-6.0
79  DO i=1,12
80  last=mod(ia*last+ic,im)
81  unorm=unorm+float(last)/float(im)
82  END DO
83  END
84 
85 
86  SUBROUTINE initgl(NAGBAR,NALCAR,NSTD,IPRLIM)
87 * initialization of package
88 * NAGB = number of global parameters
89 * DERGB(1) ... DERGB(NAGB) = derivatives w.r.t. global parameters
90 * NALC = number of local parameters (maximum)
91 * DERLC(1) ... DERLC(NALC) = derivatives w.r.t. local parameters
92 * ------------------------------------------------------------------
93 * Basic dimension parameters
94  parameter(mglobl=1400,mlocal=10,nstore=10000,mcs=30) ! dimensions
95  parameter(mgl=mglobl+mcs) ! derived parameter
96 * derived parameters
97  parameter(msymgb =(mglobl*mglobl+mglobl)/2,
98  + msym =(mgl*mgl+mgl)/2,
99  + msymlc=(mlocal*mlocal+mlocal)/2,
100  + mrecta= mglobl*mlocal,
101  + mglocs= mglobl*mcs,
102  + msymcs= (mcs*mcs+mcs)/2 )
103  DOUBLE PRECISION CGMAT,CLMAT,CLCMAT,BGVEC,BLVEC,
104  + corrm,corrv,summ,diag,scdiag,pparm,dparm,adercs,
105  + arhs
106  LOGICAL SCFLAG
107  common/lsqred/cgmat(msym),clmat(msymlc),clcmat(mrecta),
108  + diag(mgl),bgvec(mgl),blvec(mlocal),
109  + corrm(msymgb),corrv(mglobl),psigm(mglobl),
110  + pparm(mglobl),adercs(mglocs),arhs(mcs),
111  + dparm(mglobl),
112  + scdiag(mglobl),summ,scflag(mglobl),
113  + indgb(mglobl),indlc(mlocal),loctot,locrej,
114  + nagb,nalc,nsum,nhist,mhist(51),khist(51),lhist(51),
115  + nst,nfl,indst(nstore),arest(nstore),itert,lunit,ncs,
116  + nlnpa(mglobl),nstdev,cfactr,icnpr,icnlim
117 * ------------------------------------------------------------------
118  INTEGER NDR(7)
119  DATA ndr/1,2,5,10,20,50,100/
120 * ...
121  icnlim=iprlim
122  IF(icnlim.GE.0) WRITE(*,199)
123  199 FORMAT(
124  +' '/
125  +' * o o o '/
126  +' o o o '/
127  +' o ooooo o o o oo ooo oo ooo oo '/
128  +' o o o o o o o o o o o o o o o o '/
129  +' o o o o o o oooo o o oooo o o oooo '/
130  +' o o o o o o o ooo o o o o '/
131  +' o o o o oo oo oo o oo ooo oo starting'/
132  +' o ')
133  itert=0 ! reset iteration counter/flag
134  lunit=0 ! unit for scratch storage
135  cfact=1.0! factor for cut during iterations
136  ncs =0 ! number of constraints
137  icnpr=0 ! number of printouts
138  loctot=0 ! total number of local fits
139  locrej=0 ! number of rejected fits
140  nstdev=max(0,min(nstd,3)) ! 0 ... 4
141  cfactr=1.0
142  IF (nstd.GT.nstdev) cfactr=(float(nstd)/float(nstdev))**2
143  nagb=nagbar
144  nalc=nalcar
145  IF(icnlim.GE.0) THEN
146  WRITE(*,*) ' '
147  WRITE(*,*) 'Number of global parameters ',nagb
148  WRITE(*,*) 'Number of local parameters ',nalc
149  WRITE(*,*) ' Number of standard deviations ',nstdev
150  WRITE(*,*) ' Initial cut factor ',cfactr
151  WRITE(*,*) ' Number of test printouts ',icnlim
152  END IF
153  IF(nagb.GT.mglobl.OR.nalc.GT.mlocal) THEN
154  WRITE(*,*) 'Too many parameter - STOP'
155  stop
156  END IF
157  IF(nstdev.NE.0.AND.icnlim.GE.0) THEN
158  WRITE(*,*) 'Final cut corresponds to',nstdev,
159  + ' standard deviations.'
160  WRITE(*,*) 'The actual cuts are made in Chisquare/Ndf at:'
161  DO i=1,7
162  WRITE(*,101) ndr(i),chindl(nstdev,ndr(i))
163  101 FORMAT(20x,'Ndf =',i4,' Limit =',f8.2)
164  END DO
165  END IF
166 * reset matrices for global variables
167  DO i=1,nagb
168  bgvec(i)=0.0
169  pparm(i)=0.0 ! previous values of parameters set to zero
170  dparm(i)=0.0 ! corrections of parameters zero
171  psigm(i)=-1.0 ! no sigma defined for parameter I
172  nlnpa(i)=0 ! linear parameter by default
173  END DO
174  DO i=1,(nagb*nagb+nagb)/2
175  cgmat(i)=0.0
176  END DO
177 * reset histogram
178  nhist=0
179  DO i=1,51
180  mhist(i)=0
181  khist(i)=0
182  lhist(i)=0
183  END DO
184 * reset matrices for local variables
185  summ=0.0d0
186  nsum=0
187  DO i=1,nalc
188  blvec(i)=0.0
189  END DO
190  DO i=1,(nalc*nalc+nalc)/2
191  clmat(i)=0.0
192  END DO
193 * reset matrices for mixed variables
194 C DO I=1,NAGB*NALC
195 C CLCMAT(I)=0.0
196 C END DO
197  nst=0 ! reset counter for derivatives
198  nfl=0 ! reset flag for derivative storage
199  END
200 
201  SUBROUTINE parglo(PAR)
202 * optional: initialize parameters with nonzero values
203 * ------------------------------------------------------------------
204 * Basic dimension parameters
205  parameter(mglobl=1400,mlocal=10,nstore=10000,mcs=30) ! dimensions
206  parameter(mgl=mglobl+mcs) ! derived parameter
207 * derived parameters
208  parameter(msymgb =(mglobl*mglobl+mglobl)/2,
209  + msym =(mgl*mgl+mgl)/2,
210  + msymlc=(mlocal*mlocal+mlocal)/2,
211  + mrecta= mglobl*mlocal,
212  + mglocs= mglobl*mcs,
213  + msymcs= (mcs*mcs+mcs)/2 )
214  DOUBLE PRECISION CGMAT,CLMAT,CLCMAT,BGVEC,BLVEC,
215  + corrm,corrv,summ,diag,scdiag,pparm,dparm,adercs,
216  + arhs
217  LOGICAL SCFLAG
218  common/lsqred/cgmat(msym),clmat(msymlc),clcmat(mrecta),
219  + diag(mgl),bgvec(mgl),blvec(mlocal),
220  + corrm(msymgb),corrv(mglobl),psigm(mglobl),
221  + pparm(mglobl),adercs(mglocs),arhs(mcs),
222  + dparm(mglobl),
223  + scdiag(mglobl),summ,scflag(mglobl),
224  + indgb(mglobl),indlc(mlocal),loctot,locrej,
225  + nagb,nalc,nsum,nhist,mhist(51),khist(51),lhist(51),
226  + nst,nfl,indst(nstore),arest(nstore),itert,lunit,ncs,
227  + nlnpa(mglobl),nstdev,cfactr,icnpr,icnlim
228 * ------------------------------------------------------------------
229  REAL PAR(*)
230  DO i=1,nagb
231  pparm(i)=par(i)
232  END DO
233  END
234 
235  SUBROUTINE parsig(INDEX,SIGMA)
236 * optional: define sigma for single parameter
237 * ------------------------------------------------------------------
238 * Basic dimension parameters
239  parameter(mglobl=1400,mlocal=10,nstore=10000,mcs=30) ! dimensions
240  parameter(mgl=mglobl+mcs) ! derived parameter
241 * derived parameters
242  parameter(msymgb =(mglobl*mglobl+mglobl)/2,
243  + msym =(mgl*mgl+mgl)/2,
244  + msymlc=(mlocal*mlocal+mlocal)/2,
245  + mrecta= mglobl*mlocal,
246  + mglocs= mglobl*mcs,
247  + msymcs= (mcs*mcs+mcs)/2 )
248  DOUBLE PRECISION CGMAT,CLMAT,CLCMAT,BGVEC,BLVEC,
249  + corrm,corrv,summ,diag,scdiag,pparm,dparm,adercs,
250  + arhs
251  LOGICAL SCFLAG
252  common/lsqred/cgmat(msym),clmat(msymlc),clcmat(mrecta),
253  + diag(mgl),bgvec(mgl),blvec(mlocal),
254  + corrm(msymgb),corrv(mglobl),psigm(mglobl),
255  + pparm(mglobl),adercs(mglocs),arhs(mcs),
256  + dparm(mglobl),
257  + scdiag(mglobl),summ,scflag(mglobl),
258  + indgb(mglobl),indlc(mlocal),loctot,locrej,
259  + nagb,nalc,nsum,nhist,mhist(51),khist(51),lhist(51),
260  + nst,nfl,indst(nstore),arest(nstore),itert,lunit,ncs,
261  + nlnpa(mglobl),nstdev,cfactr,icnpr,icnlim
262 * ------------------------------------------------------------------
263  IF(index.LT.1.OR.index.GT.nagb) RETURN
264  IF(sigma.LT.0.0) RETURN
265  psigm(index)=sigma
266  END
267 
268 CHK
269  SUBROUTINE parget(INDEX,PAR,SIG)
270 * optional: get value, sigma for single parameter
271 * ------------------------------------------------------------------
272 * Basic dimension parameters
273  parameter(mglobl=1400,mlocal=10,nstore=10000,mcs=30) ! dimensions
274  parameter(mgl=mglobl+mcs) ! derived parameter
275 * derived parameters
276  parameter(msymgb =(mglobl*mglobl+mglobl)/2,
277  + msym =(mgl*mgl+mgl)/2,
278  + msymlc=(mlocal*mlocal+mlocal)/2,
279  + mrecta= mglobl*mlocal,
280  + mglocs= mglobl*mcs,
281  + msymcs= (mcs*mcs+mcs)/2 )
282  DOUBLE PRECISION CGMAT,CLMAT,CLCMAT,BGVEC,BLVEC,
283  + corrm,corrv,summ,diag,scdiag,pparm,dparm,adercs,
284  + arhs
285  LOGICAL SCFLAG
286  common/lsqred/cgmat(msym),clmat(msymlc),clcmat(mrecta),
287  + diag(mgl),bgvec(mgl),blvec(mlocal),
288  + corrm(msymgb),corrv(mglobl),psigm(mglobl),
289  + pparm(mglobl),adercs(mglocs),arhs(mcs),
290  + dparm(mglobl),
291  + scdiag(mglobl),summ,scflag(mglobl),
292  + indgb(mglobl),indlc(mlocal),loctot,locrej,
293  + nagb,nalc,nsum,nhist,mhist(51),khist(51),lhist(51),
294  + nst,nfl,indst(nstore),arest(nstore),itert,lunit,ncs,
295  + nlnpa(mglobl),nstdev,cfactr,icnpr,icnlim
296 * ------------------------------------------------------------------
297  par=0.0
298  sig=0.0
299  IF(index.LT.1.OR.index.GT.nagb) RETURN
300  par=pparm(index)
301  sig=psigm(index)
302  END
303 
304  SUBROUTINE nonlin(INDEX)
305 * optional: set nonlinear flag for single parameter
306 * ------------------------------------------------------------------
307 * Basic dimension parameters
308  parameter(mglobl=1400,mlocal=10,nstore=10000,mcs=30) ! dimensions
309  parameter(mgl=mglobl+mcs) ! derived parameter
310 * derived parameters
311  parameter(msymgb =(mglobl*mglobl+mglobl)/2,
312  + msym =(mgl*mgl+mgl)/2,
313  + msymlc=(mlocal*mlocal+mlocal)/2,
314  + mrecta= mglobl*mlocal,
315  + mglocs= mglobl*mcs,
316  + msymcs= (mcs*mcs+mcs)/2 )
317  DOUBLE PRECISION CGMAT,CLMAT,CLCMAT,BGVEC,BLVEC,
318  + corrm,corrv,summ,diag,scdiag,pparm,dparm,adercs,
319  + arhs
320  LOGICAL SCFLAG
321  common/lsqred/cgmat(msym),clmat(msymlc),clcmat(mrecta),
322  + diag(mgl),bgvec(mgl),blvec(mlocal),
323  + corrm(msymgb),corrv(mglobl),psigm(mglobl),
324  + pparm(mglobl),adercs(mglocs),arhs(mcs),
325  + dparm(mglobl),
326  + scdiag(mglobl),summ,scflag(mglobl),
327  + indgb(mglobl),indlc(mlocal),loctot,locrej,
328  + nagb,nalc,nsum,nhist,mhist(51),khist(51),lhist(51),
329  + nst,nfl,indst(nstore),arest(nstore),itert,lunit,ncs,
330  + nlnpa(mglobl),nstdev,cfactr,icnpr,icnlim
331 * ------------------------------------------------------------------
332  IF(index.LT.1.OR.index.GT.nagb) RETURN
333  nlnpa(index)=1
334  END
335 
336  SUBROUTINE initun(LUN,CUTFAC)
337 * optional: unit for iterations
338 * ------------------------------------------------------------------
339 * Basic dimension parameters
340  parameter(mglobl=1400,mlocal=10,nstore=10000,mcs=30) ! dimensions
341  parameter(mgl=mglobl+mcs) ! derived parameter
342 * derived parameters
343  parameter(msymgb =(mglobl*mglobl+mglobl)/2,
344  + msym =(mgl*mgl+mgl)/2,
345  + msymlc=(mlocal*mlocal+mlocal)/2,
346  + mrecta= mglobl*mlocal,
347  + mglocs= mglobl*mcs,
348  + msymcs= (mcs*mcs+mcs)/2 )
349  DOUBLE PRECISION CGMAT,CLMAT,CLCMAT,BGVEC,BLVEC,
350  + corrm,corrv,summ,diag,scdiag,pparm,dparm,adercs,
351  + arhs
352  LOGICAL SCFLAG
353  common/lsqred/cgmat(msym),clmat(msymlc),clcmat(mrecta),
354  + diag(mgl),bgvec(mgl),blvec(mlocal),
355  + corrm(msymgb),corrv(mglobl),psigm(mglobl),
356  + pparm(mglobl),adercs(mglocs),arhs(mcs),
357  + dparm(mglobl),
358  + scdiag(mglobl),summ,scflag(mglobl),
359  + indgb(mglobl),indlc(mlocal),loctot,locrej,
360  + nagb,nalc,nsum,nhist,mhist(51),khist(51),lhist(51),
361  + nst,nfl,indst(nstore),arest(nstore),itert,lunit,ncs,
362  + nlnpa(mglobl),nstdev,cfactr,icnpr,icnlim
363 * ------------------------------------------------------------------
364  LOGICAL EXS,OPN
365 * test file for existence
366  INQUIRE(unit=lun,opened=opn,iostat=ios)
367  IF(ios.NE.0) THEN
368  stop '<INITUN: Inquire error>'
369  END IF
370  iostat=0
371  IF(.NOT.opn) THEN
372  OPEN(unit=lun,form='UNFORMATTED',status='SCRATCH',iostat=ios)
373  IF(ios.NE.0) THEN
374  stop '<INITUN: Open error>'
375  END IF
376  IF(icnlim.GE.0) WRITE(*,*) 'Scratch file opened'
377  END IF
378  lunit=lun
379  cfactr=max(1.0,cutfac) ! > 1.0
380  IF(icnlim.GE.0) WRITE(*,*) 'Initial cut factor is',cfactr
381  itert=1 ! iteration 1 is first iteration
382  END
383 
384  SUBROUTINE constf(DERCS,RHS)
385 * optional: constraints
386 * ------------------------------------------------------------------
387 * Basic dimension parameters
388  parameter(mglobl=1400,mlocal=10,nstore=10000,mcs=30) ! dimensions
389  parameter(mgl=mglobl+mcs) ! derived parameter
390 * derived parameters
391  parameter(msymgb =(mglobl*mglobl+mglobl)/2,
392  + msym =(mgl*mgl+mgl)/2,
393  + msymlc=(mlocal*mlocal+mlocal)/2,
394  + mrecta= mglobl*mlocal,
395  + mglocs= mglobl*mcs,
396  + msymcs= (mcs*mcs+mcs)/2 )
397  DOUBLE PRECISION CGMAT,CLMAT,CLCMAT,BGVEC,BLVEC,
398  + corrm,corrv,summ,diag,scdiag,pparm,dparm,adercs,
399  + arhs
400  LOGICAL SCFLAG
401  common/lsqred/cgmat(msym),clmat(msymlc),clcmat(mrecta),
402  + diag(mgl),bgvec(mgl),blvec(mlocal),
403  + corrm(msymgb),corrv(mglobl),psigm(mglobl),
404  + pparm(mglobl),adercs(mglocs),arhs(mcs),
405  + dparm(mglobl),
406  + scdiag(mglobl),summ,scflag(mglobl),
407  + indgb(mglobl),indlc(mlocal),loctot,locrej,
408  + nagb,nalc,nsum,nhist,mhist(51),khist(51),lhist(51),
409  + nst,nfl,indst(nstore),arest(nstore),itert,lunit,ncs,
410  + nlnpa(mglobl),nstdev,cfactr,icnpr,icnlim
411 * ------------------------------------------------------------------
412  REAL DERCS(*)
413  IF(ncs.GE.mcs) stop '<INITCS> too many constraints'
414  DO i=1,nagb
415  adercs(nagb*ncs+i)=dercs(i)
416  END DO
417  ncs=ncs+1
418  arhs(ncs)=rhs
419  END
420 
421  SUBROUTINE equloc(DERGB,DERLC,RRMEAS,SIGMA)
422 * a single equation with its derivatives
423 * DERGB(1) ... DERGB(NAGB) = derivatives w.r.t. global parameters
424 * DERLC(1) ... DERLC(NALC) = derivatives w.r.t. local parameters
425 * RMEAS = measured value
426 * SIGMA = standard deviation
427 * (WGHT = weight = 1/SIGMA**2)
428 * ------------------------------------------------------------------
429 * Basic dimension parameters
430  parameter(mglobl=1400,mlocal=10,nstore=10000,mcs=30) ! dimensions
431  parameter(mgl=mglobl+mcs) ! derived parameter
432 * derived parameters
433  parameter(msymgb =(mglobl*mglobl+mglobl)/2,
434  + msym =(mgl*mgl+mgl)/2,
435  + msymlc=(mlocal*mlocal+mlocal)/2,
436  + mrecta= mglobl*mlocal,
437  + mglocs= mglobl*mcs,
438  + msymcs= (mcs*mcs+mcs)/2 )
439  DOUBLE PRECISION CGMAT,CLMAT,CLCMAT,BGVEC,BLVEC,
440  + corrm,corrv,summ,diag,scdiag,pparm,dparm,adercs,
441  + arhs
442  LOGICAL SCFLAG
443  common/lsqred/cgmat(msym),clmat(msymlc),clcmat(mrecta),
444  + diag(mgl),bgvec(mgl),blvec(mlocal),
445  + corrm(msymgb),corrv(mglobl),psigm(mglobl),
446  + pparm(mglobl),adercs(mglocs),arhs(mcs),
447  + dparm(mglobl),
448  + scdiag(mglobl),summ,scflag(mglobl),
449  + indgb(mglobl),indlc(mlocal),loctot,locrej,
450  + nagb,nalc,nsum,nhist,mhist(51),khist(51),lhist(51),
451  + nst,nfl,indst(nstore),arest(nstore),itert,lunit,ncs,
452  + nlnpa(mglobl),nstdev,cfactr,icnpr,icnlim
453 * ------------------------------------------------------------------
454  REAL DERGB(*),DERLC(*)
455 * ...
456  rmeas=rrmeas
457  IF(sigma.LE.0.0) THEN
458  DO i=1,nalc ! local parameters
459  derlc(i)=0.0 ! reset
460  END DO
461  DO i=1,nagb ! global parameters
462  dergb(i)=0.0 ! reset
463  END DO
464  RETURN
465  END IF
466  wght=1.0/sigma**2
467  nonzer=0
468  ialc=0
469  iblc=-1
470  DO i=1,nalc ! count number of local parameters
471  IF(derlc(i).NE.0.0) THEN
472  nonzer=nonzer+1
473  IF(ialc.EQ.0) ialc=i ! first and last index
474  iblc=i
475  END IF
476  END DO
477  iagb=0
478  ibgb=-1
479  DO i=1,nagb ! ... plus global parameters
480  IF(dergb(i).NE.0.0) THEN
481  nonzer=nonzer+1
482  IF(iagb.EQ.0) iagb=i ! first and last index
483  ibgb=i
484  END IF
485  END DO
486  IF(nst+nonzer+2.GE.nstore) THEN
487  nfl=1 ! set overflow flag
488  RETURN ! ignore data
489  END IF
490  nst=nst+1
491  indst(nst)=0
492  arest(nst)=rmeas
493  DO i=ialc,iblc ! local parameters
494  IF(derlc(i).NE.0.0) THEN
495  nst=nst+1
496  indst(nst)=i ! store index ...
497  arest(nst)=derlc(i) ! ... and value of nonzero derivative
498  derlc(i)=0.0 ! reset
499  END IF
500  END DO
501  nst=nst+1
502  indst(nst)=0
503  arest(nst)=wght
504  DO i=iagb,ibgb ! global parameters
505  IF(dergb(i).NE.0.0) THEN
506  nst=nst+1
507  indst(nst)=i ! store index ...
508  arest(nst)=dergb(i) ! ... and value of nonzero derivative
509  dergb(i)=0.0 ! reset
510  END IF
511  END DO
512  END
513 
514  SUBROUTINE zerloc(DERGB,DERLC)
515 * reset derivatives
516 * DERGB(1) ... DERGB(NAGB) = derivatives w.r.t. global parameters
517 * DERLC(1) ... DERLC(NALC) = derivatives w.r.t. local parameters
518 * ------------------------------------------------------------------
519 * Basic dimension parameters
520  parameter(mglobl=1400,mlocal=10,nstore=10000,mcs=30) ! dimensions
521  parameter(mgl=mglobl+mcs) ! derived parameter
522 * derived parameters
523  parameter(msymgb =(mglobl*mglobl+mglobl)/2,
524  + msym =(mgl*mgl+mgl)/2,
525  + msymlc=(mlocal*mlocal+mlocal)/2,
526  + mrecta= mglobl*mlocal,
527  + mglocs= mglobl*mcs,
528  + msymcs= (mcs*mcs+mcs)/2 )
529  DOUBLE PRECISION CGMAT,CLMAT,CLCMAT,BGVEC,BLVEC,
530  + corrm,corrv,summ,diag,scdiag,pparm,dparm,adercs,
531  + arhs
532  LOGICAL SCFLAG
533  common/lsqred/cgmat(msym),clmat(msymlc),clcmat(mrecta),
534  + diag(mgl),bgvec(mgl),blvec(mlocal),
535  + corrm(msymgb),corrv(mglobl),psigm(mglobl),
536  + pparm(mglobl),adercs(mglocs),arhs(mcs),
537  + dparm(mglobl),
538  + scdiag(mglobl),summ,scflag(mglobl),
539  + indgb(mglobl),indlc(mlocal),loctot,locrej,
540  + nagb,nalc,nsum,nhist,mhist(51),khist(51),lhist(51),
541  + nst,nfl,indst(nstore),arest(nstore),itert,lunit,ncs,
542  + nlnpa(mglobl),nstdev,cfactr,icnpr,icnlim
543 * ------------------------------------------------------------------
544  REAL DERGB(*),DERLC(*)
545  DO i=1,nalc ! local parameters
546  derlc(i)=0.0 ! reset
547  END DO
548  DO i=1,nagb ! global parameters
549  dergb(i)=0.0 ! reset
550  END DO
551  END
552 
553 
554  SUBROUTINE fitloc
555 * fit after end of local block - faster(?) version
556 * ------------------------------------------------------------------
557 * Basic dimension parameters
558  parameter(mglobl=1400,mlocal=10,nstore=10000,mcs=30) ! dimensions
559  parameter(mgl=mglobl+mcs) ! derived parameter
560 * derived parameters
561  parameter(msymgb =(mglobl*mglobl+mglobl)/2,
562  + msym =(mgl*mgl+mgl)/2,
563  + msymlc=(mlocal*mlocal+mlocal)/2,
564  + mrecta= mglobl*mlocal,
565  + mglocs= mglobl*mcs,
566  + msymcs= (mcs*mcs+mcs)/2 )
567  DOUBLE PRECISION CGMAT,CLMAT,CLCMAT,BGVEC,BLVEC,
568  + corrm,corrv,summ,diag,scdiag,pparm,dparm,adercs,
569  + arhs
570  LOGICAL SCFLAG
571  common/lsqred/cgmat(msym),clmat(msymlc),clcmat(mrecta),
572  + diag(mgl),bgvec(mgl),blvec(mlocal),
573  + corrm(msymgb),corrv(mglobl),psigm(mglobl),
574  + pparm(mglobl),adercs(mglocs),arhs(mcs),
575  + dparm(mglobl),
576  + scdiag(mglobl),summ,scflag(mglobl),
577  + indgb(mglobl),indlc(mlocal),loctot,locrej,
578  + nagb,nalc,nsum,nhist,mhist(51),khist(51),lhist(51),
579  + nst,nfl,indst(nstore),arest(nstore),itert,lunit,ncs,
580  + nlnpa(mglobl),nstdev,cfactr,icnpr,icnlim,
581  + indnz(mglobl),indbk(mglobl)
582 *
583 *KEEP,CSFSPY.
584  dimension parloc(7)
585  REAL*8 DCVLOC(7,7)
586  COMMON / csfspy / nrank, ndf, rms, parloc, dcvloc
587 * ------------------------------------------------------------
588 *KEND.
589 * ------------------------------------------------------------------
590  icnpr=icnpr+1
591  incrp=0
592  IF(itert.EQ.1) THEN
593  IF(nst.GT.0) THEN ! write to scratch file
594  WRITE(lunit) nst,(indst(i),i=1,nst),(arest(i),i=1,nst)
595  END IF
596  END IF
597 * reset matrices for local variables ooooooooooooooooooooooooooooooo
598  summ=0.0d0
599  nsum=0
600  DO i=1,nalc
601  blvec(i)=0.0
602  END DO
603  DO i=1,(nalc*nalc+nalc)/2
604  clmat(i)=0.0
605  END DO
606 * reset matrices for mixed variables
607 C DO I=1,NAGB*NALC
608 C CLCMAT(I)=0.0
609 C END DO
610 * reset pointer matrix for mixed variables
611  DO i=1,nagb
612  indnz(i)=0
613  END DO
614  nagbn=0 ! actual number of global parameters
615 C above code new
616 * normal equations - symmetric matrix for local parameters
617  IF(nst.LE.1) GOTO 100
618  ist=0
619  10 ja =0 ! new equation
620  jb =0
621  20 ist=ist+1
622  IF(ist.GT.nst.OR.indst(ist).EQ.0) THEN
623  IF(ja.EQ.0) THEN
624  ja=ist ! first zero: measured value
625  ELSE IF(jb.EQ.0) THEN
626  jb=ist ! second zero: weight
627  ELSE
628  ist=ist-1 ! end of equation
629  rmeas=arest(ja) ! use the data
630 * subtract global ... from measured value
631  DO j=1,ist-jb
632  ij=indst(jb+j) ! subtract from measured value ...
633  IF(nlnpa(ij).EQ.0) THEN ! ... for linear parameter
634  rmeas=rmeas-arest(jb+j)*(pparm(ij)+dparm(ij))
635  ELSE ! ... for nonlinear parameter
636  rmeas=rmeas-arest(jb+j)*dparm(ij)
637  END IF
638  END DO
639 * end-subtract
640  wght =arest(jb) ! ... and the weight
641  DO j=1,jb-ja-1 ! number of derivatives
642  ij=indst(ja+j)
643  blvec(ij)=blvec(ij)+wght*rmeas*arest(ja+j)
644  DO k=1,j
645  ik=indst(ja+k)
646  jk=(ij*ij-ij)/2+ik
647  clmat(jk)=clmat(jk)+wght*arest(ja+j)*arest(ja+k)
648  END DO
649  END DO
650 
651  IF(ist.EQ.nst) GOTO 21
652  GOTO 10 ! next equation
653  END IF
654  END IF
655  IF(ist.LE.nst) GOTO 20
656 * end of data - determine local fit parameters
657  21 CALL spminv(clmat,blvec,nalc,nrank,scdiag,scflag)
658  ! inversion and solution
659 CHK
660  DO kk=1,7
661  parloc(kk)=sngl(blvec(kk))
662  ENDDO
663  kk=0
664  DO k1=1,7
665  DO k2=1,k1
666  kk=kk+1
667  dcvloc(k1,k2)=clmat(kk)
668  dcvloc(k2,k1)=clmat(kk)
669  ENDDO
670  ENDDO
671 CHK
672  IF(icnpr.LE.icnlim) THEN
673  WRITE(*,*) ' '
674  WRITE(*,*) '__________________________________________________'
675  WRITE(*,*) 'Printout of local fit (FITLOC) with rank=',nrank
676  WRITE(*,*) ' Result of local fit: (Index/Parameter/error)'
677  WRITE(*,103) (j,blvec(j),sqrt(clmat((j*j+j/2))),j=1,jb-ja-1)
678  103 FORMAT(2(i6,2g12.4))
679  END IF
680 * calculate residuals
681  summ=0.0d0
682  nsum=0
683 * second loop for residual calculation -----------------------------
684  ist=0
685  30 ja =0 ! new equation
686  jb =0
687  40 ist=ist+1
688  IF(ist.GT.nst.OR.indst(ist).EQ.0) THEN
689  IF(ja.EQ.0) THEN
690  ja=ist ! first zero: measured value
691  ELSE IF(jb.EQ.0) THEN
692  jb=ist ! second zero: weight
693  ELSE
694  ist=ist-1 ! end of equation
695  IF(icnpr.LE.icnlim.AND.incrp.LE.11) THEN
696  nderlc=jb-ja-1
697  ndergl=ist-jb
698  incrp=incrp+1
699  WRITE(*,*) ' '
700  WRITE(*,*) incrp,'. equation: ',
701  + 'measured value ',arest(ja),' +-',1.0/sqrt(arest(jb))
702  IF(incrp.LE.11) THEN
703  WRITE(*,*)
704  + ' number of derivates (global, local):',ndergl,',',nderlc
705  WRITE(*,*)
706  + ' Global derivatives are: (index/derivative/parvalue)'
707  WRITE(*,101) (indst(jb+j),arest(jb+j),
708  + pparm(indst(jb+j)),j=1,ist-jb)
709  101 FORMAT(2(i6,2g12.4))
710  WRITE(*,*)
711  + ' Local derivatives are: (index/derivative)'
712  WRITE(*,102) (indst(ja+j),arest(ja+j),j=1,jb-ja-1)
713  102 FORMAT(3(i6,g14.6))
714  END IF
715  IF(incrp.EQ.11) THEN
716  WRITE(*,*) '... (+ further equations)'
717  END IF
718  END IF
719  rmeas=arest(ja) ! use the data
720 * subtract global ... from measured value
721  DO j=1,ist-jb
722  ij=indst(jb+j) ! subtract from measured value ...
723  IF(nlnpa(ij).EQ.0) THEN ! ... for linear parameter
724  rmeas=rmeas-arest(jb+j)*(pparm(ij)+dparm(ij))
725  ELSE ! ... for nonlinear parameter
726  rmeas=rmeas-arest(jb+j)*dparm(ij)
727  END IF
728  END DO
729 * end-subtract
730  wght =arest(jb) ! ... and the weight
731  DO j=1,jb-ja-1 ! number of derivatives
732  ij=indst(ja+j)
733  rmeas=rmeas-arest(ja+j)*blvec(ij)
734  END DO
735  IF(icnpr.LE.icnlim) THEN
736  WRITE(*,104) wght*rmeas**2,rmeas
737  104 FORMAT(' Chi square contribution=',f8.2,
738  + 5x,g12.4,'= residuum')
739  END IF
740 * residual histogram
741  ihbin=1.0+5.0*(rmeas*sqrt(wght)+5.0)
742  IF(ihbin.LT. 1) ihbin=51
743  IF(ihbin.GT.50) ihbin=51
744  lhist(ihbin)=lhist(ihbin)+1 ! residial histogram
745  summ=summ+wght*rmeas**2
746  nsum=nsum+1 ! ... and count equation
747  IF(ist.EQ.nst) GOTO 41
748  GOTO 30 ! next equation
749  END IF
750  END IF
751  IF(ist.LE.nst) GOTO 40
752  41 ndf=nsum-nrank
753  IF(icnpr.LE.icnlim) THEN
754  WRITE(*,*) 'Entry number',icnpr
755  WRITE(*,*) 'Final chi square, degrees of freedom',summ,ndf
756  END IF
757  rms=0.0
758  IF(ndf.GT.0) THEN ! histogram entry
759  rms=summ/float(ndf)
760  nhist=nhist+1
761  ihbin=1+5.0*rms ! bins width 0.2
762  IF(ihbin.LT. 1) ihbin= 1
763  IF(ihbin.GT.50) ihbin=51
764  mhist(ihbin)=mhist(ihbin)+1
765  END IF
766  loctot=loctot+1
767 * make eventual cut
768  IF(nstdev.NE.0.AND.ndf.GT.0) THEN
769  cutval=chindl(nstdev,ndf)*cfactr
770  IF(icnpr.LE.icnlim) THEN
771  WRITE(*,*) 'Reject if Chisq/Ndf=',rms,' > ',cutval
772  END IF
773  IF(rms.GT.cutval) THEN
774  locrej=locrej+1
775  rms=-rms
776  GOTO 100
777  END IF
778  END IF
779 * third loop for global parameters ---------------------------------
780  ist=0
781  50 ja =0 ! new equation
782  jb =0
783  60 ist=ist+1
784  IF(ist.GT.nst.OR.indst(ist).EQ.0) THEN
785  IF(ja.EQ.0) THEN
786  ja=ist ! first zero: measured value
787  ELSE IF(jb.EQ.0) THEN
788  jb=ist ! second zero: weight
789  ELSE
790  ist=ist-1 ! end of equation
791  rmeas=arest(ja) ! use the data
792 * subtract global ... from measured value
793  DO j=1,ist-jb
794  ij=indst(jb+j) ! subtract from measured value ...
795  IF(nlnpa(ij).EQ.0) THEN ! ... for linear parameter
796  rmeas=rmeas-arest(jb+j)*(pparm(ij)+dparm(ij))
797  ELSE ! ... for nonlinear parameter
798  rmeas=rmeas-arest(jb+j)*dparm(ij)
799  END IF
800  END DO
801 * end-subtract
802  wght =arest(jb) ! ... and the weight
803 * normal equations - symmetric matrix for global parameters
804  DO j=1,ist-jb
805  ij=indst(jb+j)
806  bgvec(ij)=bgvec(ij)+wght*rmeas*arest(jb+j)
807  DO k=1,j
808  ik=indst(jb+k)
809  jk=(ij*ij-ij)/2+ik
810  cgmat(jk)=cgmat(jk)+wght*arest(jb+j)*arest(jb+k)
811  END DO
812  END DO
813 * normal equations - rectangular matrix for global/local pars
814  DO j=1,ist-jb
815  ij=indst(jb+j) ! index of global variable
816 * new code start
817  ijn=indnz(ij) ! get index of index
818  IF(ijn.EQ.0) THEN
819 * new global variable - initialize matrix row
820  DO k=1,nalc
821  clcmat(nagbn*nalc+k)=0.0
822  END DO
823  nagbn=nagbn+1
824  indnz(ij)=nagbn ! insert pointer
825  indbk(nagbn)=ij ! ... and pointer back
826  ijn=nagbn
827  END IF
828 * new code end
829  DO k=1,jb-ja-1
830  ik=indst(ja+k)
831 C JK=IK+(IJ-1)*NALC ! old code
832  jk=ik+(ijn-1)*nalc ! new code
833  clcmat(jk)=clcmat(jk)+wght*arest(jb+j)*arest(ja+k)
834  END DO
835  END DO
836  IF(ist.EQ.nst) GOTO 70
837  GOTO 50 ! next equation
838  END IF
839  END IF
840  IF(ist.LE.nst) GOTO 60
841 * -------------------------------------------------------------
842 * update global matrices
843 C70 CALL SPAVAT(CLMAT,CLCMAT,CORRM,NALC,NAGB) ! correction matrix
844 C CALL SPAX(CLCMAT,BLVEC,CORRV,NAGB,NALC) ! correction vector
845 C IJ=0
846 C DO I=1,NAGB
847 C BGVEC(I)=BGVEC(I)-CORRV(I)
848 C DO J=1,I
849 C IJ=IJ+1
850 C CGMAT(IJ)=CGMAT(IJ)-CORRM(IJ)
851 C END DO
852 C END DO
853 * new code
854  70 CALL spavat(clmat,clcmat,corrm,nalc,nagbn) ! correction matrix
855  CALL spax(clcmat,blvec,corrv,nagbn,nalc) ! correction vector
856  ijn=0
857  DO in=1,nagbn
858  i=indbk(in) ! get pointer back to global index
859  bgvec(i)=bgvec(i)-corrv(in)
860  DO jn=1,in
861  j=indbk(jn)
862  ijn=ijn+1
863  IF(i.GE.j) THEN
864  ij=j+(i*i-i)/2
865  ELSE
866  ij=i+(j*j-j)/2
867  END IF
868  cgmat(ij)=cgmat(ij)-corrm(ijn)
869  END DO
870  END DO
871 * end new code
872  entry killoc
873  100 IF(itert.LE.1) THEN
874 * histogram of used store space
875  IF(nfl.EQ.0) THEN
876  ibin=1.0+50.0*float(nst)/float(nstore)
877  ibin=min(ibin,50)
878  khist(ibin)=khist(ibin)+1
879  ELSE
880  khist(51)=khist(51)+1
881  END IF
882  END IF
883  nst=0 ! reset counter
884  nfl=0 ! reset overflow flag
885  END
886 
887 
888  SUBROUTINE fitglo(PAR)
889 * final global fit
890  REAL PAR(*)
891  CHARACTER*2 PATEXT
892 * ------------------------------------------------------------------
893 * Basic dimension parameters
894  parameter(mglobl=1400,mlocal=10,nstore=10000,mcs=30) ! dimensions
895  parameter(mgl=mglobl+mcs) ! derived parameter
896 * derived parameters
897  parameter(msymgb =(mglobl*mglobl+mglobl)/2,
898  + msym =(mgl*mgl+mgl)/2,
899  + msymlc=(mlocal*mlocal+mlocal)/2,
900  + mrecta= mglobl*mlocal,
901  + mglocs= mglobl*mcs,
902  + msymcs= (mcs*mcs+mcs)/2 )
903  DOUBLE PRECISION CGMAT,CLMAT,CLCMAT,BGVEC,BLVEC,
904  + corrm,corrv,summ,diag,scdiag,pparm,dparm,adercs,
905  + arhs
906  LOGICAL SCFLAG
907  common/lsqred/cgmat(msym),clmat(msymlc),clcmat(mrecta),
908  + diag(mgl),bgvec(mgl),blvec(mlocal),
909  + corrm(msymgb),corrv(mglobl),psigm(mglobl),
910  + pparm(mglobl),adercs(mglocs),arhs(mcs),
911  + dparm(mglobl),
912  + scdiag(mglobl),summ,scflag(mglobl),
913  + indgb(mglobl),indlc(mlocal),loctot,locrej,
914  + nagb,nalc,nsum,nhist,mhist(51),khist(51),lhist(51),
915  + nst,nfl,indst(nstore),arest(nstore),itert,lunit,ncs,
916  + nlnpa(mglobl),nstdev,cfactr,icnpr,icnlim
917 * ------------------------------------------------------------------
918  IF(itert.LE.1) itelim=10 ! maximum number of iterations
919  IF(itert.NE.0) THEN
920  lthist=0
921  DO i=1,51
922  lthist=lthist+lhist(i)
923  END DO
924  IF(icnlim.GE.0) THEN
925  WRITE(*,*) ' Initial residual histogram:',
926  + ' Total',lthist,' entries with',lhist(51),
927  + ' overflows'
928  CALL pxhist(lhist,50,-5.0,+5.0) ! histogram printout
929  END IF
930  END IF
931  khtot=0
932  DO i=1,51
933  khtot=khtot+khist(i)
934  END DO
935  IF(icnlim.GE.0) THEN
936  WRITE(*,*) ' '
937  WRITE(*,*) 'Histogram of used local fit storage:',
938  + ' total',khtot,' local fits with',khist(51),
939  + ' overflows'
940  END IF
941  IF(khist(51).NE.0.AND.icnlim.GE.0) THEN
942  WRITE(*,*) 'Parameter NSTORE to small! Change and rerun!'
943  END IF
944  IF(icnlim.GE.0) THEN
945  CALL pxhist(khist,50,0.0,float(nstore)) ! histogram printout
946  END IF
947  10 IF(icnlim.GE.0) THEN
948  WRITE(*,*) ' '
949  WRITE(*,*) '... making global fit ...'
950  END IF
951  nn=0
952  ii=0 ! modify matrix acccording to PSIGM value
953  DO i=1,nagb
954  ll=ii
955  ii=ii+i
956  IF(psigm(i).EQ.0.0) THEN
957  nn=nn+1 ! count number of fixed parameters
958  DO j=1,nagb
959  ll=ll+1
960  IF(j.GT.i) ll=ll+j-2
961  cgmat(ll)=0.0 ! reset row and column for parameter
962  END DO
963  ELSE IF(psigm(i).GT.0.0) THEN
964  cgmat(ii)=cgmat(ii)+1.0/psigm(i)**2 ! add to diagonal
965  END IF
966  END DO
967  ii=0 ! start saving diagonal elements
968  DO i=1,nagb
969  ii=ii+i
970  diag(i)=cgmat(ii) ! save original diagonal elements
971  END DO
972  nvar=nagb
973  IF(ncs.NE.0) THEN ! add constraints
974  ii=(nagb*nagb+nagb)/2
975  DO i=1,ncs ! loop on constraints
976  sum=arhs(i)
977  DO j=1,nagb
978  cgmat(ii+j)=adercs(nagb*(i-1)+j)
979  sum=sum-adercs(nagb*(i-1)+j)*(pparm(j)+dparm(j))
980  END DO
981  DO j=1,i
982  cgmat(ii+nagb+j)=0.0
983  END DO
984  nvar=nvar+1
985  ii =ii+nvar
986  bgvec(nvar)=sum
987  END DO
988  END IF
989 * =================================================
990  CALL spminv(cgmat,bgvec,nvar,nrank,scdiag,scflag)!matrix inversion
991  ndefec=nvar-nrank-nn ! rank defect
992 * =====================================================
993  DO i=1,nagb
994  dparm(i)=dparm(i)+bgvec(i) ! accumulate corrections
995  END DO
996  CALL prtglo(66)
997  IF(icnlim.GE.0) THEN
998  WRITE(*,*) 'The rank defect of the symmetric',nvar,'-by-',nvar,
999  + ' matrix is ',ndefec,' (should be zero).'
1000  END IF
1001  IF(itert.EQ.0.OR.nstdev.EQ.0.OR.itert.GE.itelim) GOTO 90
1002 * iterations
1003  IF(icnlim.GE.0) THEN
1004  WRITE(*,*) ' '
1005  WRITE(*,*) ' Total',loctot,' local fits,',locrej,' rejected.'
1006  WRITE(*,*) ' Histogram of Chisq/Ndf:',
1007  + ' Total',nhist,' entries with',mhist(51),
1008  + ' overflows'
1009  CALL pxhist(mhist,50,0.0,10.0) ! histogram printout
1010  END IF
1011 * reset histogram
1012  nhist=0
1013  DO i=1,51
1014  mhist(i)=0
1015  lhist(i)=0
1016  END DO
1017  itert=itert+1
1018  loctot=0
1019  locrej=0
1020  IF(cfactr.NE.1.0) THEN
1021  cfactr=sqrt(cfactr)
1022  IF(cfactr.LT.1.2) THEN
1023  cfactr=1.0
1024  itelim=itert+1
1025  END IF
1026  END IF
1027  IF(icnlim.GE.0) WRITE(*,107) itert,cfactr
1028  107 FORMAT(' Iteration',i3,' with cut factor=',f6.2)
1029 * reset matrices for global variables
1030  DO i=1,nagb
1031  bgvec(i)=0.0
1032  END DO
1033  DO i=1,(nagb*nagb+nagb)/2
1034  cgmat(i)=0.0
1035  END DO
1036  rewind lunit
1037  20 READ(lunit,end=10) nst,(indst(i),i=1,nst),(arest(i),i=1,nst)
1038  CALL fitloc
1039  GOTO 20
1040 * ==================================================================
1041  90 IF(icnlim.GE.0) THEN
1042  WRITE(*,*) ' '
1043  WRITE(*,*) ' Result of fit for global parameters'
1044  WRITE(*,*) ' ==================================='
1045  WRITE(*,101)
1046  END IF
1047  ii=0
1048  DO i=1,nagb
1049  ii=ii+i
1050  err=sqrt(abs(cgmat(ii)))
1051  IF(cgmat((i*i+i)/2).LT.0.0) err=-err
1052  gcor=0.0
1053  IF(cgmat(ii)*diag(i).GT.0.0) THEN
1054 * global correlation
1055  gcor=sqrt(abs(1.0-1.0/(cgmat(ii)*diag(i))))
1056  END IF
1057  IF(i.LE.25.OR.nagb-i.LE.25) THEN
1058  patext=' '
1059  IF(nlnpa(i).NE.0) patext='nl'
1060  IF(icnlim.GE.0) WRITE(*,102) i,patext,
1061  + pparm(i),pparm(i)+dparm(i),dparm(i),bgvec(i),err,gcor
1062  END IF
1063  par(i)=pparm(i)+dparm(i) ! copy of result to array in argument
1064  END DO
1065  DO i=1,ncs ! constraints
1066  IF(i.EQ.1) WRITE(*,*) ' '
1067  sum=0.0
1068  DO j=1,nagb
1069  ppp=pparm(j)+dparm(j)
1070  sum=sum+ppp*adercs(nagb*(i-1)+j)
1071 CHK IF(ICNLIM.GE.0) WRITE(*,*) J,PPP,ADERCS(NAGB*(I-1)+J),SUM
1072  END DO
1073  WRITE(*,106) i,sum,arhs(i),sum-arhs(i)
1074  END DO
1075  IF(icnlim.GE.0) THEN
1076  WRITE(*,*) ' '
1077  WRITE(*,*) ' Total',loctot,' local fits,',locrej,' rejected.'
1078  WRITE(*,*) ' Histogram of Chisq/Ndf:',
1079  + ' Total',nhist,' entries with',mhist(51),
1080  + ' overflows'
1081  CALL pxhist(mhist,50,0.0,10.0) ! histogram printout
1082  END IF
1083  lthist=0
1084  DO i=1,51
1085  lthist=lthist+lhist(i)
1086  END DO
1087  IF(icnlim.GE.0) THEN
1088  WRITE(*,*) ' Residual histogram:',
1089  + ' Total',lthist,' entries with',lhist(51),
1090  + ' overflows'
1091  CALL pxhist(lhist,50,-5.0,+5.0) ! histogram printout
1092  WRITE(*,199)
1093  END IF
1094  199 FORMAT(
1095  +' '/
1096  +' * o o o '/
1097  +' o o o '/
1098  +' o ooooo o o o oo ooo oo ooo oo '/
1099  +' o o o o o o o o o o o o o o o o '/
1100  +' o o o o o o oooo o o oooo o o oooo '/
1101  +' o o o o o o o ooo o o o o '/
1102  +' o o o o oo oo oo o oo ooo oo ending.'/
1103  +' o ')
1104  101 FORMAT(1x,' I initial final differ',
1105  + ' lastcor Error glcor'/
1106  + 1x,' --- ----------- ----------- -----------',
1107  + ' ----------- -------- -----')
1108  102 FORMAT(1x,i4,1x,a2,1x,4f12.5,f9.5,f6.3)
1109  106 FORMAT(' Constraint',i2,' Sum - RHS =',g12.5,' -',g12.5,
1110  + ' = ',g12.5)
1111  END
1112 
1113  FUNCTION errpar(I)
1114 * return error for parameter I
1115 * ------------------------------------------------------------------
1116 * Basic dimension parameters
1117  parameter(mglobl=1400,mlocal=10,nstore=10000,mcs=30) ! dimensions
1118  parameter(mgl=mglobl+mcs) ! derived parameter
1119 * derived parameters
1120  parameter(msymgb =(mglobl*mglobl+mglobl)/2,
1121  + msym =(mgl*mgl+mgl)/2,
1122  + msymlc=(mlocal*mlocal+mlocal)/2,
1123  + mrecta= mglobl*mlocal,
1124  + mglocs= mglobl*mcs,
1125  + msymcs= (mcs*mcs+mcs)/2 )
1126  DOUBLE PRECISION CGMAT,CLMAT,CLCMAT,BGVEC,BLVEC,
1127  + corrm,corrv,summ,diag,scdiag,pparm,dparm,adercs,
1128  + arhs
1129  LOGICAL SCFLAG
1130  common/lsqred/cgmat(msym),clmat(msymlc),clcmat(mrecta),
1131  + diag(mgl),bgvec(mgl),blvec(mlocal),
1132  + corrm(msymgb),corrv(mglobl),psigm(mglobl),
1133  + pparm(mglobl),adercs(mglocs),arhs(mcs),
1134  + dparm(mglobl),
1135  + scdiag(mglobl),summ,scflag(mglobl),
1136  + indgb(mglobl),indlc(mlocal),loctot,locrej,
1137  + nagb,nalc,nsum,nhist,mhist(51),khist(51),lhist(51),
1138  + nst,nfl,indst(nstore),arest(nstore),itert,lunit,ncs,
1139  + nlnpa(mglobl),nstdev,cfactr,icnpr,icnlim
1140 * ------------------------------------------------------------------
1141  errpar=0.0
1142  IF(i.LE.0.OR.i.GT.nagb) RETURN
1143  ii=(i*i+i)/2
1144  errpar=sqrt(abs(cgmat(ii)))
1145  IF(cgmat(ii).LT.0.0) errpar=errpar
1146  END
1147 
1148  FUNCTION corpar(I,J)
1149 * return correlation between parameters I and J
1150 * ------------------------------------------------------------------
1151 * Basic dimension parameters
1152  parameter(mglobl=1400,mlocal=10,nstore=10000,mcs=30) ! dimensions
1153  parameter(mgl=mglobl+mcs) ! derived parameter
1154 * derived parameters
1155  parameter(msymgb =(mglobl*mglobl+mglobl)/2,
1156  + msym =(mgl*mgl+mgl)/2,
1157  + msymlc=(mlocal*mlocal+mlocal)/2,
1158  + mrecta= mglobl*mlocal,
1159  + mglocs= mglobl*mcs,
1160  + msymcs= (mcs*mcs+mcs)/2 )
1161  DOUBLE PRECISION CGMAT,CLMAT,CLCMAT,BGVEC,BLVEC,
1162  + corrm,corrv,summ,diag,scdiag,pparm,dparm,adercs,
1163  + arhs
1164  LOGICAL SCFLAG
1165  common/lsqred/cgmat(msym),clmat(msymlc),clcmat(mrecta),
1166  + diag(mgl),bgvec(mgl),blvec(mlocal),
1167  + corrm(msymgb),corrv(mglobl),psigm(mglobl),
1168  + pparm(mglobl),adercs(mglocs),arhs(mcs),
1169  + dparm(mglobl),
1170  + scdiag(mglobl),summ,scflag(mglobl),
1171  + indgb(mglobl),indlc(mlocal),loctot,locrej,
1172  + nagb,nalc,nsum,nhist,mhist(51),khist(51),lhist(51),
1173  + nst,nfl,indst(nstore),arest(nstore),itert,lunit,ncs,
1174  + nlnpa(mglobl),nstdev,cfactr,icnpr,icnlim
1175 * ------------------------------------------------------------------
1176  corpar=0.0
1177  IF(i.LE.0.OR.i.GT.nagb) RETURN
1178  IF(j.LE.0.OR.j.GT.nagb) RETURN
1179  IF(i.EQ.j) RETURN
1180  ii=(i*i+i)/2
1181  jj=(j*j+j)/2
1182  k=max(i,j)
1183  ij=(k*k-k)/2+min(i,j)
1184  err=sqrt(abs(cgmat(ii)*cgmat(jj)))
1185  IF(err.NE.0.0) corpar=cgmat(ij)/err
1186  END
1187 
1188  SUBROUTINE prtglo(LUN)
1189 * ------------------------------------------------------------------
1190 * Basic dimension parameters
1191  parameter(mglobl=1400,mlocal=10,nstore=10000,mcs=30) ! dimensions
1192  parameter(mgl=mglobl+mcs) ! derived parameter
1193 * derived parameters
1194  parameter(msymgb =(mglobl*mglobl+mglobl)/2,
1195  + msym =(mgl*mgl+mgl)/2,
1196  + msymlc=(mlocal*mlocal+mlocal)/2,
1197  + mrecta= mglobl*mlocal,
1198  + mglocs= mglobl*mcs,
1199  + msymcs= (mcs*mcs+mcs)/2 )
1200  DOUBLE PRECISION CGMAT,CLMAT,CLCMAT,BGVEC,BLVEC,
1201  + corrm,corrv,summ,diag,scdiag,pparm,dparm,adercs,
1202  + arhs
1203  LOGICAL SCFLAG
1204  common/lsqred/cgmat(msym),clmat(msymlc),clcmat(mrecta),
1205  + diag(mgl),bgvec(mgl),blvec(mlocal),
1206  + corrm(msymgb),corrv(mglobl),psigm(mglobl),
1207  + pparm(mglobl),adercs(mglocs),arhs(mcs),
1208  + dparm(mglobl),
1209  + scdiag(mglobl),summ,scflag(mglobl),
1210  + indgb(mglobl),indlc(mlocal),loctot,locrej,
1211  + nagb,nalc,nsum,nhist,mhist(51),khist(51),lhist(51),
1212  + nst,nfl,indst(nstore),arest(nstore),itert,lunit,ncs,
1213  + nlnpa(mglobl),nstdev,cfactr,icnpr,icnlim
1214 * ------------------------------------------------------------------
1215  CHARACTER*2 PATEXT
1216 * ...
1217  lup=lun
1218  IF(lup.EQ.0) lup=6
1219  WRITE(lup,*) ' Result of fit for global parameters'
1220  WRITE(lup,*) ' ==================================='
1221  WRITE(lup,101)
1222  ii=0
1223  DO i=1,nagb
1224  ii=ii+i
1225  err=sqrt(abs(cgmat(ii)))
1226  IF(cgmat(ii).LT.0.0) err=-err
1227  gcor=0.0
1228  IF(cgmat(ii)*diag(i).GT.0.0) THEN
1229 * global correlation
1230  gcor=sqrt(abs(1.0-1.0/(cgmat(ii)*diag(i))))
1231  END IF
1232  patext=' '
1233  IF(nlnpa(i).NE.0) patext='nl'
1234  WRITE(lup,102) i,patext,
1235  + pparm(i),pparm(i)+dparm(i),dparm(i),bgvec(i),err,gcor
1236  END DO
1237  101 FORMAT(1x,' I initial final differ',
1238  + ' lastcor Error glcor'/
1239  + 1x,' --- ----------- ----------- -----------',
1240  + ' ----------- -------- -----')
1241  102 FORMAT(1x,i4,1x,a2,1x,4f12.5,f9.5,f6.3)
1242  END
errpar
function errpar(I)
Definition: isuper.f:1114
dalign
subroutine dalign
Definition: isuper.f:26
prtglo
subroutine prtglo(LUN)
Definition: isuper.f:1189
fitloc
subroutine fitloc
Definition: isuper.f:555
initgl
subroutine initgl(NAGBAR, NALCAR, NSTD, IPRLIM)
Definition: isuper.f:87
parget
subroutine parget(INDEX, PAR, SIG)
Definition: isuper.f:270
urand
function urand()
Definition: isuper.f:67
initun
subroutine initun(LUN, CUTFAC)
Definition: isuper.f:337
parsig
subroutine parsig(INDEX, SIGMA)
Definition: isuper.f:236
chindl
real(mps) function chindl(n, nd)
Chi2/ndf cuts.
Definition: mpnum.f90:1700
zerloc
subroutine zerloc(DERGB, DERLC)
Definition: isuper.f:515
constf
subroutine constf(DERCS, RHS)
Definition: isuper.f:385
nonlin
subroutine nonlin(INDEX)
Definition: isuper.f:305
fitglo
subroutine fitglo(PAR)
Definition: isuper.f:889
corpar
function corpar(I, J)
Definition: isuper.f:1149
parglo
subroutine parglo(PAR)
Definition: isuper.f:202
unorm
function unorm()
Definition: isuper.f:75
equloc
subroutine equloc(DERGB, DERLC, RRMEAS, SIGMA)
Definition: isuper.f:422
Generated on Fri May 25 2018 13:13:40 for Millepede-II by   doxygen 1.8.13