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