Millepede-II V04-17-03
millepede1.f90
Go to the documentation of this file.
1! Code converted using TO_F90 by Alan Miller
2! Date: 2024-04-29 Time: 13:19:49
3
63
68SUBROUTINE testmp(iarg) ! test program IARG = 0 or 1
69
70 INTEGER, INTENT(IN) :: iarg
71 REAL :: dergb(10),derlc(2),par(10)
72 INTEGER, PARAMETER :: nplan=10
73 REAL :: x(nplan),y(nplan),sigma(nplan),bias(nplan),heff(nplan)
74 REAL :: frot(nplan)
75 DATA x/5.0,9.0,20.0,25.0,30.0,35.0,40.0,45.0,50.0,55.0/
76 DATA sigma/2*0.0020,8*0.0300/
77 DATA bias /0.00,-0.04,0.15,0.03,-0.075,0.045,0.035,-0.08,0.09,-0.05/
78 DATA heff/4*0.9,0.50,5*0.9/
79 ! ...
80 ncases=1000
81 CALL initgl(10,2,3,1) ! define dimension parameters
82 CALL parsig(1,0.0) ! parameter 1 is fixed
83 IF(iarg == 0) THEN
84 frot(1)=0.0
85 DO i=2,10
86 frot(i)=1.0/(x(i)-x(1))
87 END DO
88 CALL constf(frot,0.0) ! constraint: total rotation zero
89 ELSE
90 DO i=3,10
91 CALL parsig(i,0.03) ! parameters 3...8
92 END DO
93 END IF
94 CALL initun(11,10.0) ! option iterations
95 CALL zerloc(dergb,derlc) ! initialization of arrays to zero
96 ! -------------- loop
97 DO nc=1,ncases
98 ! generate straight line parameters
99 CALL gener(y,a,b,x,sigma,bias,heff)
100 DO i=1,nplan
101 ! calibrate Y(I) = A + B * X(I)
102 IF(y(i) /= 0.0) THEN
103 derlc(1)= 1.0
104 derlc(2)= x(i)
105 dergb(i)= 1.0
106 CALL equloc(dergb,derlc,y(i),sigma(i)) ! single measurement
107 END IF
108 END DO
109 CALL fitloc
110 END DO ! end loop cases
111! -------------- loop end
112CONTINUE
113CALL fitglo(par) ! final solution
114END SUBROUTINE testmp
115
117SUBROUTINE gener(y,a,b,x,sigma,bias,heff)
118
119 INTEGER, PARAMETER :: nplan=10
120
121 REAL, INTENT(OUT) :: y(*)
122 REAL, INTENT(OUT) :: a
123 REAL, INTENT(OUT) :: b
124 REAL, INTENT(IN) :: x(nplan)
125 REAL, INTENT(IN) :: sigma(nplan)
126 REAL, INTENT(IN) :: bias(nplan)
127 REAL, INTENT(IN) :: heff(nplan)
128
129 ! ...
130 ! generate straight line parameters
131 a=20.0*zrand()-10.0 ! A = -10 ... +10 uniform
132 b= 0.3*znorm() ! B = -0.3 ... +0.3 gaussian
133 DO i=1,nplan ! planes
134 IF(zrand() < heff(i)) THEN
135 y(i)=a+b*x(i) + bias(i)+sigma(i)*znorm()
136 ELSE
137 y(i)=0.0 ! unmeasured
138 END IF
139 END DO
140END SUBROUTINE gener
141
143FUNCTION zrand()
144 ! (simple generator, showing principle)
145 parameter(ia=205,ic=29573,im=139968)
146 DATA last/4711/
147 last=mod(ia*last+ic,im)
148 IF(last == 0) last=mod(ia*last+ic,im)
149 zrand=float(last)/float(im)
150END
151
153FUNCTION znorm()
154 ! (simple generator, showing principle)
155 parameter(ia=205,ic=29573,im=139968)
156 DATA last/4711/
157 znorm=-6.0
158 DO i=1,12
159 last=mod(ia*last+ic,im)
160 znorm=znorm+float(last)/float(im)
161 END DO
162END
163
170SUBROUTINE initgl(nagbar,nalcar,nstd,iprlim)
171 ! ------------------------------------------------------------------
172 ! Basic dimension parameters
173 parameter(mglobl=1400,mlocal=10,nstore=10000,mcs=10) ! dimensions
174 parameter(mgl=mglobl+mcs) ! derived parameter
175 ! derived parameters
176 parameter(msymgb =(mglobl*mglobl+mglobl)/2, msym =(mgl*mgl+mgl)/2, &
177 msymlc=(mlocal*mlocal+mlocal)/2, mrecta= mglobl*mlocal, &
178 mglocs= mglobl*mcs, msymcs= (mcs*mcs+mcs)/2 )
179 DOUBLE PRECISION :: cgmat,clmat,clcmat,bgvec,blvec, &
180 corrm,corrv,summ,diag,scdiag,pparm,dparm,adercs, arhs
181 LOGICAL :: scflag
182 common/lsqred/cgmat(msym),clmat(msymlc),clcmat(mrecta), &
183 diag(mgl),bgvec(mgl),blvec(mlocal), &
184 corrm(msymgb),corrv(mglobl),psigm(mglobl), &
185 pparm(mglobl),adercs(mglocs),arhs(mcs), dparm(mglobl), &
186 scdiag(mglobl),summ,scflag(mglobl), &
187 indgb(mglobl),indlc(mlocal),loctot,locrej, &
188 nagb,nalc,nsum,nhist,mhist(51),khist(51),lhist(51), &
189 nst,nfl,indst(nstore),arest(nstore),itert,lunit,ncs, &
190 nlnpa(mglobl),nstdev,cfactr,icnpr,icnlim, indnz(mglobl),indbk(mglobl)
191 ! ------------------------------------------------------------------
192 INTEGER :: ndr(7)
193 DATA ndr/1,2,5,10,20,50,100/
194 ! ...
195 icnlim=iprlim
196 IF(icnlim >= 0) WRITE(*,199)
197199 FORMAT( ' '/ &
198 ' * o o o '/ &
199 ' o o o '/ &
200 ' o ooooo o o o oo ooo oo ooo oo '/ &
201 ' o o o o o o o o o o o o o o o o '/ &
202 ' o o o o o o oooo o o oooo o o oooo '/ &
203 ' o o o o o o o ooo o o o o '/ &
204 ' o o o o oo oo oo o oo ooo oo starting'/ &
205 ' o ')
206 itert=0 ! reset iteration counter/flag
207 lunit=0 ! unit for scratch storage
208 cfact=1.0! factor for cut during iterations
209 ncs =0 ! number of constraints
210 icnpr=0 ! number of printouts
211 loctot=0 ! total number of local fits
212 locrej=0 ! number of rejected fits
213 nstdev=max(0,min(nstd,3)) ! 0 ... 4
214 cfactr=1.0
215 nagb=nagbar
216 nalc=nalcar
217 IF(icnlim >= 0) THEN
218 WRITE(*,*) ' '
219 WRITE(*,*) 'Number of global parameters ',nagb
220 WRITE(*,*) 'Number of local parameters ',nalc
221 WRITE(*,*) ' Number of standard deviations ',nstdev
222 WRITE(*,*) ' Number of test printouts ',icnlim
223 END IF
224 IF(nagb > mglobl.OR.nalc > mlocal) THEN
225 WRITE(*,*) 'Too many parameter - STOP'
226 stop
227 END IF
228 IF(nstdev /= 0.AND.icnlim >= 0) THEN
229 WRITE(*,*) 'Final cut corresponds to',nstdev, ' standard deviations.'
230 WRITE(*,*) 'The actual cuts are made in Chisquare/Ndf at:'
231 DO i=1,7
232 WRITE(*,101) ndr(i),chindl(nstdev,ndr(i))
233101 FORMAT(20x,'Ndf =',i4,' Limit =',f8.2)
234 END DO
235 END IF
236 ! reset matrices for global variables
237 DO i=1,nagb
238 bgvec(i)=0.0
239 pparm(i)=0.0 ! previous values of parameters set to zero
240 dparm(i)=0.0 ! corrections of parameters zero
241 psigm(i)=-1.0 ! no sigma defined for parameter I
242 nlnpa(i)=0 ! linear parameter by default
243 END DO
244 DO i=1,(nagb*nagb+nagb)/2
245 cgmat(i)=0.0
246 END DO
247 ! reset histogram
248 nhist=0
249 DO i=1,51
250 mhist(i)=0
251 khist(i)=0
252 lhist(i)=0
253 END DO
254 ! reset matrices for local variables
255 summ=0.0d0
256 nsum=0
257 DO i=1,nalc
258 blvec(i)=0.0
259 END DO
260 DO i=1,(nalc*nalc+nalc)/2
261 clmat(i)=0.0
262 END DO
263 nst=0 ! reset counter for derivatives
264 nfl=0 ! reset flag for derivative storage
265END SUBROUTINE initgl
266
270SUBROUTINE parglo(par)
271 ! ------------------------------------------------------------------
272 ! Basic dimension parameters
273 parameter(mglobl=1400,mlocal=10,nstore=10000,mcs=10) ! dimensions
274 parameter(mgl=mglobl+mcs) ! derived parameter
275 ! derived parameters
276 parameter(msymgb =(mglobl*mglobl+mglobl)/2, msym =(mgl*mgl+mgl)/2, &
277 msymlc=(mlocal*mlocal+mlocal)/2, mrecta= mglobl*mlocal, &
278 mglocs= mglobl*mcs, msymcs= (mcs*mcs+mcs)/2 )
279 DOUBLE PRECISION :: cgmat,clmat,clcmat,bgvec,blvec, &
280 corrm,corrv,summ,diag,scdiag,pparm,dparm,adercs, arhs
281 LOGICAL :: scflag
282 common/lsqred/cgmat(msym),clmat(msymlc),clcmat(mrecta), &
283 diag(mgl),bgvec(mgl),blvec(mlocal), &
284 corrm(msymgb),corrv(mglobl),psigm(mglobl), &
285 pparm(mglobl),adercs(mglocs),arhs(mcs), dparm(mglobl), &
286 scdiag(mglobl),summ,scflag(mglobl), &
287 indgb(mglobl),indlc(mlocal),loctot,locrej, &
288 nagb,nalc,nsum,nhist,mhist(51),khist(51),lhist(51), &
289 nst,nfl,indst(nstore),arest(nstore),itert,lunit,ncs, &
290 nlnpa(mglobl),nstdev,cfactr,icnpr,icnlim, indnz(mglobl),indbk(mglobl)
291 ! ------------------------------------------------------------------
292 REAL :: par(*)
293 DO i=1,nagb
294 pparm(i)=par(i)
295 END DO
296END SUBROUTINE parglo
297
299SUBROUTINE parsig(INDEX,sigma)
300 ! ------------------------------------------------------------------
301 ! Basic dimension parameters
302 parameter(mglobl=1400,mlocal=10,nstore=10000,mcs=10) ! dimensions
303 parameter(mgl=mglobl+mcs) ! derived parameter
304 ! derived parameters
305 parameter(msymgb =(mglobl*mglobl+mglobl)/2, msym =(mgl*mgl+mgl)/2, &
306 msymlc=(mlocal*mlocal+mlocal)/2, mrecta= mglobl*mlocal, &
307 mglocs= mglobl*mcs, msymcs= (mcs*mcs+mcs)/2 )
308 DOUBLE PRECISION :: cgmat,clmat,clcmat,bgvec,blvec, &
309 corrm,corrv,summ,diag,scdiag,pparm,dparm,adercs, arhs
310 LOGICAL :: scflag
311 common/lsqred/cgmat(msym),clmat(msymlc),clcmat(mrecta), &
312 diag(mgl),bgvec(mgl),blvec(mlocal), &
313 corrm(msymgb),corrv(mglobl),psigm(mglobl), &
314 pparm(mglobl),adercs(mglocs),arhs(mcs), dparm(mglobl), &
315 scdiag(mglobl),summ,scflag(mglobl), &
316 indgb(mglobl),indlc(mlocal),loctot,locrej, &
317 nagb,nalc,nsum,nhist,mhist(51),khist(51),lhist(51), &
318 nst,nfl,indst(nstore),arest(nstore),itert,lunit,ncs, &
319 nlnpa(mglobl),nstdev,cfactr,icnpr,icnlim, indnz(mglobl),indbk(mglobl)
320 ! ------------------------------------------------------------------
321 IF(index < 1.OR.index > nagb) RETURN
322 IF(sigma < 0.0) RETURN
323 psigm(index)=sigma
324END SUBROUTINE parsig
325
327SUBROUTINE nonlin(INDEX)
328 ! ------------------------------------------------------------------
329 ! Basic dimension parameters
330 parameter(mglobl=1400,mlocal=10,nstore=10000,mcs=10) ! dimensions
331 parameter(mgl=mglobl+mcs) ! derived parameter
332 ! derived parameters
333 parameter(msymgb =(mglobl*mglobl+mglobl)/2, msym =(mgl*mgl+mgl)/2, &
334 msymlc=(mlocal*mlocal+mlocal)/2, mrecta= mglobl*mlocal, &
335 mglocs= mglobl*mcs, msymcs= (mcs*mcs+mcs)/2 )
336 DOUBLE PRECISION :: cgmat,clmat,clcmat,bgvec,blvec, &
337 corrm,corrv,summ,diag,scdiag,pparm,dparm,adercs, arhs
338 LOGICAL :: scflag
339 common/lsqred/cgmat(msym),clmat(msymlc),clcmat(mrecta), &
340 diag(mgl),bgvec(mgl),blvec(mlocal), &
341 corrm(msymgb),corrv(mglobl),psigm(mglobl), &
342 pparm(mglobl),adercs(mglocs),arhs(mcs), dparm(mglobl), &
343 scdiag(mglobl),summ,scflag(mglobl), &
344 indgb(mglobl),indlc(mlocal),loctot,locrej, &
345 nagb,nalc,nsum,nhist,mhist(51),khist(51),lhist(51), &
346 nst,nfl,indst(nstore),arest(nstore),itert,lunit,ncs, &
347 nlnpa(mglobl),nstdev,cfactr,icnpr,icnlim, indnz(mglobl),indbk(mglobl)
348 ! ------------------------------------------------------------------
349 IF(index < 1.OR.index > nagb) RETURN
350 nlnpa(index)=1
351END SUBROUTINE nonlin
352
354SUBROUTINE initun(lun,cutfac)
355 ! ------------------------------------------------------------------
356 ! Basic dimension parameters
357 parameter(mglobl=1400,mlocal=10,nstore=10000,mcs=10) ! dimensions
358 parameter(mgl=mglobl+mcs) ! derived parameter
359 ! derived parameters
360 parameter(msymgb =(mglobl*mglobl+mglobl)/2, msym =(mgl*mgl+mgl)/2, &
361 msymlc=(mlocal*mlocal+mlocal)/2, mrecta= mglobl*mlocal, &
362 mglocs= mglobl*mcs, msymcs= (mcs*mcs+mcs)/2 )
363 DOUBLE PRECISION :: cgmat,clmat,clcmat,bgvec,blvec, &
364 corrm,corrv,summ,diag,scdiag,pparm,dparm,adercs, arhs
365 LOGICAL :: scflag
366 common/lsqred/cgmat(msym),clmat(msymlc),clcmat(mrecta), &
367 diag(mgl),bgvec(mgl),blvec(mlocal), &
368 corrm(msymgb),corrv(mglobl),psigm(mglobl), &
369 pparm(mglobl),adercs(mglocs),arhs(mcs), dparm(mglobl), &
370 scdiag(mglobl),summ,scflag(mglobl), &
371 indgb(mglobl),indlc(mlocal),loctot,locrej, &
372 nagb,nalc,nsum,nhist,mhist(51),khist(51),lhist(51), &
373 nst,nfl,indst(nstore),arest(nstore),itert,lunit,ncs, &
374 nlnpa(mglobl),nstdev,cfactr,icnpr,icnlim, indnz(mglobl),indbk(mglobl)
375 ! ------------------------------------------------------------------
376 LOGICAL :: opn
377 ! test file for existence
378 INQUIRE(unit=lun,opened=opn,iostat=ios)
379 IF(ios /= 0) THEN
380 stop '<INITUN: Inquire error>'
381 END IF
382 iostat=0
383 IF(.NOT.opn) THEN
384 OPEN(unit=lun,form='UNFORMATTED',status='SCRATCH',iostat=ios)
385 IF(ios /= 0) THEN
386 stop '<INITUN: Open error>'
387 END IF
388 IF(icnlim >= 0) WRITE(*,*) 'Scratch file opened'
389 END IF
390 lunit=lun
391 cfactr=max(1.0,cutfac) ! > 1.0
392 IF(icnlim >= 0) WRITE(*,*) 'Initial cut factor is',cfactr
393 itert=1 ! iteration 1 is first iteration
394END SUBROUTINE initun
395
397SUBROUTINE constf(dercs,rhs)
398 ! ------------------------------------------------------------------
399 ! Basic dimension parameters
400 parameter(mglobl=1400,mlocal=10,nstore=10000,mcs=10) ! dimensions
401 parameter(mgl=mglobl+mcs) ! derived parameter
402 ! derived parameters
403 parameter(msymgb =(mglobl*mglobl+mglobl)/2, msym =(mgl*mgl+mgl)/2, &
404 msymlc=(mlocal*mlocal+mlocal)/2, mrecta= mglobl*mlocal, &
405 mglocs= mglobl*mcs, msymcs= (mcs*mcs+mcs)/2 )
406 DOUBLE PRECISION :: cgmat,clmat,clcmat,bgvec,blvec, &
407 corrm,corrv,summ,diag,scdiag,pparm,dparm,adercs, arhs
408 LOGICAL :: scflag
409 common/lsqred/cgmat(msym),clmat(msymlc),clcmat(mrecta), &
410 diag(mgl),bgvec(mgl),blvec(mlocal), &
411 corrm(msymgb),corrv(mglobl),psigm(mglobl), &
412 pparm(mglobl),adercs(mglocs),arhs(mcs), dparm(mglobl), &
413 scdiag(mglobl),summ,scflag(mglobl), &
414 indgb(mglobl),indlc(mlocal),loctot,locrej, &
415 nagb,nalc,nsum,nhist,mhist(51),khist(51),lhist(51), &
416 nst,nfl,indst(nstore),arest(nstore),itert,lunit,ncs, &
417 nlnpa(mglobl),nstdev,cfactr,icnpr,icnlim, indnz(mglobl),indbk(mglobl)
418 ! ------------------------------------------------------------------
419 REAL :: dercs(*)
420 IF(ncs >= mcs) stop '<INITCS> too many constraints'
421 DO i=1,nagb
422 adercs(nagb*ncs+i)=dercs(i)
423 END DO
424 ncs=ncs+1
425 arhs(ncs)=rhs
426 WRITE(*,*) 'Number of constraints increased to ',ncs
427END SUBROUTINE constf
428
436SUBROUTINE equloc(dergb,derlc,rrmeas,sigma)
437 ! ------------------------------------------------------------------
438 ! Basic dimension parameters
439 parameter(mglobl=1400,mlocal=10,nstore=10000,mcs=10) ! dimensions
440 parameter(mgl=mglobl+mcs) ! derived parameter
441 ! derived parameters
442 parameter(msymgb =(mglobl*mglobl+mglobl)/2, msym =(mgl*mgl+mgl)/2, &
443 msymlc=(mlocal*mlocal+mlocal)/2, mrecta= mglobl*mlocal, &
444 mglocs= mglobl*mcs, msymcs= (mcs*mcs+mcs)/2 )
445 DOUBLE PRECISION :: cgmat,clmat,clcmat,bgvec,blvec, &
446 corrm,corrv,summ,diag,scdiag,pparm,dparm,adercs, arhs
447 LOGICAL :: scflag
448 common/lsqred/cgmat(msym),clmat(msymlc),clcmat(mrecta), &
449 diag(mgl),bgvec(mgl),blvec(mlocal), &
450 corrm(msymgb),corrv(mglobl),psigm(mglobl), &
451 pparm(mglobl),adercs(mglocs),arhs(mcs), dparm(mglobl), &
452 scdiag(mglobl),summ,scflag(mglobl), &
453 indgb(mglobl),indlc(mlocal),loctot,locrej, &
454 nagb,nalc,nsum,nhist,mhist(51),khist(51),lhist(51), &
455 nst,nfl,indst(nstore),arest(nstore),itert,lunit,ncs, &
456 nlnpa(mglobl),nstdev,cfactr,icnpr,icnlim, indnz(mglobl),indbk(mglobl)
457 ! ------------------------------------------------------------------
458 REAL :: dergb(*),derlc(*)
459 ! ...
460 rmeas=rrmeas
461 IF(sigma <= 0.0) THEN
462 DO i=1,nalc ! local parameters
463 derlc(i)=0.0 ! reset
464 END DO
465 DO i=1,nagb ! global parameters
466 dergb(i)=0.0 ! reset
467 END DO
468 RETURN
469 END IF
470 wght=1.0/sigma**2
471 nonzer=0
472 ialc=0
473 iblc=-1
474 DO i=1,nalc ! count number of local parameters
475 IF(derlc(i) /= 0.0) THEN
476 nonzer=nonzer+1
477 IF(ialc == 0) ialc=i ! first and last index
478 iblc=i
479 END IF
480 END DO
481 iagb=0
482 ibgb=-1
483 DO i=1,nagb ! ... plus global parameters
484 IF(dergb(i) /= 0.0) THEN
485 nonzer=nonzer+1
486 IF(iagb == 0) iagb=i ! first and last index
487 ibgb=i
488 END IF
489 END DO
490 IF(nst+nonzer+2 >= nstore) THEN
491 nfl=1 ! set overflow flag
492 RETURN ! ignore data
493 END IF
494 nst=nst+1
495 indst(nst)=0
496 arest(nst)=rmeas
497 DO i=ialc,iblc ! local parameters
498 IF(derlc(i) /= 0.0) THEN
499 nst=nst+1
500 indst(nst)=i ! store index ...
501 arest(nst)=derlc(i) ! ... and value of nonzero derivative
502 derlc(i)=0.0 ! reset
503 END IF
504 END DO
505 nst=nst+1
506 indst(nst)=0
507 arest(nst)=wght
508 DO i=iagb,ibgb ! global parameters
509 IF(dergb(i) /= 0.0) THEN
510 nst=nst+1
511 indst(nst)=i ! store index ...
512 arest(nst)=dergb(i) ! ... and value of nonzero derivative
513 dergb(i)=0.0 ! reset
514 END IF
515 END DO
516END SUBROUTINE equloc
517
522SUBROUTINE zerloc(dergb,derlc)
523 ! ------------------------------------------------------------------
524 ! Basic dimension parameters
525 parameter(mglobl=1400,mlocal=10,nstore=10000,mcs=10) ! dimensions
526 parameter(mgl=mglobl+mcs) ! derived parameter
527 ! derived parameters
528 parameter(msymgb =(mglobl*mglobl+mglobl)/2, msym =(mgl*mgl+mgl)/2, &
529 msymlc=(mlocal*mlocal+mlocal)/2, mrecta= mglobl*mlocal, &
530 mglocs= mglobl*mcs, msymcs= (mcs*mcs+mcs)/2 )
531 DOUBLE PRECISION :: cgmat,clmat,clcmat,bgvec,blvec, &
532 corrm,corrv,summ,diag,scdiag,pparm,dparm,adercs, arhs
533 LOGICAL :: scflag
534 common/lsqred/cgmat(msym),clmat(msymlc),clcmat(mrecta), &
535 diag(mgl),bgvec(mgl),blvec(mlocal), &
536 corrm(msymgb),corrv(mglobl),psigm(mglobl), &
537 pparm(mglobl),adercs(mglocs),arhs(mcs), 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, indnz(mglobl),indbk(mglobl)
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
551END SUBROUTINE zerloc
552
554SUBROUTINE fitloc
555 ! faster(?) version
556 ! ------------------------------------------------------------------
557 ! Basic dimension parameters
558 parameter(mglobl=1400,mlocal=10,nstore=10000,mcs=10) ! dimensions
559 parameter(mgl=mglobl+mcs) ! derived parameter
560 ! derived parameters
561 parameter(msymgb =(mglobl*mglobl+mglobl)/2, msym =(mgl*mgl+mgl)/2, &
562 msymlc=(mlocal*mlocal+mlocal)/2, mrecta= mglobl*mlocal, &
563 mglocs= mglobl*mcs, msymcs= (mcs*mcs+mcs)/2 )
564 DOUBLE PRECISION :: cgmat,clmat,clcmat,bgvec,blvec, &
565 corrm,corrv,summ,diag,scdiag,pparm,dparm,adercs, arhs
566 LOGICAL :: scflag
567 common/lsqred/cgmat(msym),clmat(msymlc),clcmat(mrecta), &
568 diag(mgl),bgvec(mgl),blvec(mlocal), &
569 corrm(msymgb),corrv(mglobl),psigm(mglobl), &
570 pparm(mglobl),adercs(mglocs),arhs(mcs), dparm(mglobl), &
571 scdiag(mglobl),summ,scflag(mglobl), &
572 indgb(mglobl),indlc(mlocal),loctot,locrej, &
573 nagb,nalc,nsum,nhist,mhist(51),khist(51),lhist(51), &
574 nst,nfl,indst(nstore),arest(nstore),itert,lunit,ncs, &
575 nlnpa(mglobl),nstdev,cfactr,icnpr,icnlim, indnz(mglobl),indbk(mglobl)
576 ! ------------------------------------------------------------------
577 icnpr=icnpr+1
578 incrp=0
579 IF(itert == 1) THEN
580 IF(nst > 0) THEN ! write to scratch file
581 WRITE(lunit) nst,(indst(i),i=1,nst),(arest(i),i=1,nst)
582 END IF
583 END IF
584 ! reset matrices for local variables ooooooooooooooooooooooooooooooo
585 summ=0.0d0
586 nsum=0
587 DO i=1,nalc
588 blvec(i)=0.0
589 END DO
590 DO i=1,(nalc*nalc+nalc)/2
591 clmat(i)=0.0
592 END DO
593 ! reset pointer matrix for mixed variables
594 DO i=1,nagb
595 indnz(i)=0
596 END DO
597 nagbn=0 ! actual number of global parameters
598 ! normal equations - symmetric matrix for local parameters
599 IF(nst <= 1) GO TO 100
600 ist=0
60110 ja =0 ! new equation
602 jb =0
60320 ist=ist+1
604 IF(ist > nst.OR.indst(ist) == 0) THEN
605 IF(ja == 0) THEN
606 ja=ist ! first zero: measured value
607 ELSE IF(jb == 0) THEN
608 jb=ist ! second zero: weight
609 ELSE
610 ist=ist-1 ! end of equation
611 rmeas=arest(ja) ! use the data
612 ! subtract global ... from measured value
613 DO j=1,ist-jb
614 ij=indst(jb+j) ! subtract from measured value ...
615 IF(nlnpa(ij) == 0) THEN ! ... for linear parameter
616 rmeas=rmeas-arest(jb+j)*real(pparm(ij)+dparm(ij))
617 ELSE ! ... for nonlinear parameter
618 rmeas=rmeas-arest(jb+j)*real(dparm(ij))
619 END IF
620 END DO
621 ! end-subtract
622 wght =arest(jb) ! ... and the weight
623 DO j=1,jb-ja-1 ! number of derivatives
624 ij=indst(ja+j)
625 blvec(ij)=blvec(ij)+wght*rmeas*arest(ja+j)
626 DO k=1,j
627 ik=indst(ja+k)
628 jk=(ij*ij-ij)/2+ik
629 clmat(jk)=clmat(jk)+wght*arest(ja+j)*arest(ja+k)
630 END DO
631 END DO
632
633 IF(ist == nst) GO TO 21
634 GO TO 10 ! next equation
635 END IF
636 END IF
637 IF(ist <= nst) GO TO 20
638 ! end of data - determine local fit parameters
63921 CALL spminv(clmat,blvec,nalc,nrank,scdiag,scflag)
640 ! inversion and solution
641 IF(icnpr <= icnlim) THEN
642 WRITE(*,*) ' '
643 WRITE(*,*) '__________________________________________________'
644 WRITE(*,*) 'Printout of local fit (FITLOC) with rank=',nrank
645 WRITE(*,*) ' Result of local fit: (Index/Parameter/error)'
646 WRITE(*,103) (j,blvec(j),sqrt(clmat((j*j+j)/2)),j=1,jb-ja-1)
647103 FORMAT(2(i6,2g12.4))
648 END IF
649 ! calculate residuals
650 summ=0.0d0
651 nsum=0
652 ! second loop for residual calculation -----------------------------
653 ist=0
65430 ja =0 ! new equation
655 jb =0
65640 ist=ist+1
657 IF(ist > nst.OR.indst(ist) == 0) THEN
658 IF(ja == 0) THEN
659 ja=ist ! first zero: measured value
660 ELSE IF(jb == 0) THEN
661 jb=ist ! second zero: weight
662 ELSE
663 ist=ist-1 ! end of equation
664 IF(icnpr <= icnlim.AND.incrp <= 11) THEN
665 nderlc=jb-ja-1
666 ndergl=ist-jb
667 incrp=incrp+1
668 WRITE(*,*) ' '
669 WRITE(*,*) incrp,'. equation: ', &
670 'measured value ',arest(ja),' +-',1.0/sqrt(arest(jb))
671 IF(incrp <= 11) THEN
672 WRITE(*,*) ' number of derivates (global, local):',ndergl,',',nderlc
673 WRITE(*,*) ' Global derivatives are: (index/derivative/parvalue)'
674 WRITE(*,101) (indst(jb+j),arest(jb+j), pparm(indst(jb+j)),j=1,ist-jb)
675101 FORMAT(2(i6,2g12.4))
676 WRITE(*,*) ' Local derivatives are: (index/derivative)'
677 WRITE(*,102) (indst(ja+j),arest(ja+j),j=1,jb-ja-1)
678102 FORMAT(3(i6,g14.6))
679 END IF
680 IF(incrp == 11) THEN
681 WRITE(*,*) '... (+ further equations)'
682 END IF
683 END IF
684 rmeas=arest(ja) ! use the data
685 ! subtract global ... from measured value
686 DO j=1,ist-jb
687 ij=indst(jb+j) ! subtract from measured value ...
688 IF(nlnpa(ij) == 0) THEN ! ... for linear parameter
689 rmeas=rmeas-arest(jb+j)*real(pparm(ij)+dparm(ij))
690 ELSE ! ... for nonlinear parameter
691 rmeas=rmeas-arest(jb+j)*real(dparm(ij))
692 END IF
693 END DO
694 ! end-subtract
695 wght =arest(jb) ! ... and the weight
696 DO j=1,jb-ja-1 ! number of derivatives
697 ij=indst(ja+j)
698 rmeas=rmeas-arest(ja+j)*real(blvec(ij))
699 END DO
700 IF(icnpr <= icnlim) THEN
701 WRITE(*,104) wght*rmeas**2,rmeas
702104 FORMAT(' Chi square contribution=',f8.2, &
703 5x,g12.4,'= residuum')
704 END IF
705 ! residual histogram
706 ihbin=int(1.0+5.0*(rmeas*sqrt(wght)+5.0))
707 IF(ihbin < 1) ihbin=51
708 IF(ihbin > 50) ihbin=51
709 lhist(ihbin)=lhist(ihbin)+1 ! residial histogram
710 summ=summ+wght*rmeas**2
711 nsum=nsum+1 ! ... and count equation
712 IF(ist == nst) GO TO 41
713 GO TO 30 ! next equation
714 END IF
715 END IF
716 IF(ist <= nst) GO TO 40
71741 ndf=nsum-nrank
718 IF(icnpr <= icnlim) THEN
719 WRITE(*,*) 'Entry number',icnpr
720 WRITE(*,*) 'Final chi square, degrees of freedom',summ,ndf
721 END IF
722 rms=0.0
723 IF(ndf > 0) THEN ! histogram entry
724 rms=real(summ)/real(ndf)
725 nhist=nhist+1
726 ihbin=int(1+5.0*rms) ! bins width 0.2
727 IF(ihbin < 1) ihbin= 1
728 IF(ihbin > 50) ihbin=51
729 mhist(ihbin)=mhist(ihbin)+1
730 END IF
731 loctot=loctot+1
732 ! make eventual cut
733 IF(nstdev /= 0.AND.ndf > 0) THEN
734 cutval=chindl(nstdev,ndf)*cfactr
735 IF(icnpr <= icnlim) THEN
736 WRITE(*,*) 'Reject if Chisq/Ndf=',rms,' > ',cutval
737 END IF
738 IF(rms > cutval) THEN
739 locrej=locrej+1
740 GO TO 100
741 END IF
742 END IF
743 ! third loop for global parameters ---------------------------------
744 ist=0
74550 ja =0 ! new equation
746 jb =0
74760 ist=ist+1
748 IF(ist > nst.OR.indst(ist) == 0) THEN
749 IF(ja == 0) THEN
750 ja=ist ! first zero: measured value
751 ELSE IF(jb == 0) THEN
752 jb=ist ! second zero: weight
753 ELSE
754 ist=ist-1 ! end of equation
755 rmeas=arest(ja) ! use the data
756 ! subtract global ... from measured value
757 DO j=1,ist-jb
758 ij=indst(jb+j) ! subtract from measured value ...
759 IF(nlnpa(ij) == 0) THEN ! ... for linear parameter
760 rmeas=rmeas-arest(jb+j)*real(pparm(ij)+dparm(ij))
761 ELSE ! ... for nonlinear parameter
762 rmeas=rmeas-arest(jb+j)*real(dparm(ij))
763 END IF
764 END DO
765 ! end-subtract
766 wght =arest(jb) ! ... and the weight
767 ! normal equations - symmetric matrix for global parameters
768 DO j=1,ist-jb
769 ij=indst(jb+j)
770 bgvec(ij)=bgvec(ij)+wght*rmeas*arest(jb+j)
771 DO k=1,j
772 ik=indst(jb+k)
773 jk=(ij*ij-ij)/2+ik
774 cgmat(jk)=cgmat(jk)+wght*arest(jb+j)*arest(jb+k)
775 END DO
776 END DO
777 ! normal equations - rectangular matrix for global/local pars
778 DO j=1,ist-jb
779 ij=indst(jb+j) ! index of global variable
780 ! new code start
781 ijn=indnz(ij) ! get index of index
782 IF(ijn == 0) THEN
783 ! new global variable - initialize matrix row
784 DO k=1,nalc
785 clcmat(nagbn*nalc+k)=0.0
786 END DO
787 nagbn=nagbn+1
788 indnz(ij)=nagbn ! insert pointer
789 indbk(nagbn)=ij ! ... and pointer back
790 ijn=nagbn
791 END IF
792 ! new code end
793 DO k=1,jb-ja-1
794 ik=indst(ja+k)
795 ! JK=IK+(IJ-1)*NALC ! old code
796 jk=ik+(ijn-1)*nalc ! new code
797 clcmat(jk)=clcmat(jk)+wght*arest(jb+j)*arest(ja+k)
798 END DO
799 END DO
800 IF(ist == nst) GO TO 70
801 GO TO 50 ! next equation
802 END IF
803 END IF
804 IF(ist <= nst) GO TO 60
805 ! -------------------------------------------------------------
806 ! update global matrices
80770 CALL spavat(clmat,clcmat,corrm,nalc,nagbn) ! correction matrix
808 CALL spax(clcmat,blvec,corrv,nagbn,nalc) ! correction vector
809 ijn=0
810 DO in=1,nagbn
811 i=indbk(in) ! get pointer back to global index
812 bgvec(i)=bgvec(i)-corrv(in)
813 DO jn=1,in
814 j=indbk(jn)
815 ijn=ijn+1
816 IF(i >= j) THEN
817 ij=j+(i*i-i)/2
818 ELSE
819 ij=i+(j*j-j)/2
820 END IF
821 cgmat(ij)=cgmat(ij)-corrm(ijn)
822 END DO
823 END DO
824 entry killoc
825100 IF(itert <= 1) THEN
826 ! histogram of used store space
827 IF(nfl == 0) THEN
828 ibin=int(1.0+50.0*float(nst)/float(nstore))
829 ibin=min(ibin,50)
830 khist(ibin)=khist(ibin)+1
831 ELSE
832 khist(51)=khist(51)+1
833 END IF
834 END IF
835 nst=0 ! reset counter
836 nfl=0 ! reset overflow flag
837END SUBROUTINE fitloc
838
840SUBROUTINE fitglo(par)
841 REAL :: par(*)
842 CHARACTER (LEN=2) :: patext
843 ! ------------------------------------------------------------------
844 ! Basic dimension parameters
845 parameter(mglobl=1400,mlocal=10,nstore=10000,mcs=10) ! 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, indnz(mglobl),indbk(mglobl)
863 ! ------------------------------------------------------------------
864 DOUBLE PRECISION :: dsum
865 IF(itert <= 1) itelim=10 ! maximum number of iterations
866 IF(itert /= 0) THEN
867 lthist=0
868 DO i=1,51
869 lthist=lthist+lhist(i)
870 END DO
871 IF(icnlim >= 0) THEN
872 WRITE(*,*) ' Initial residual histogram:', &
873 ' Total',lthist,' entries with',lhist(51), ' overflows'
874 CALL pxhist(lhist,50,-5.0,+5.0) ! histogram printout
875 END IF
876 END IF
877 khtot=0
878 DO i=1,51
879 khtot=khtot+khist(i)
880 END DO
881 IF(icnlim >= 0) THEN
882 WRITE(*,*) ' '
883 WRITE(*,*) 'Histogram of used local fit storage:', &
884 ' total',khtot,' local fits with',khist(51), ' overflows'
885 END IF
886 IF(khist(51) /= 0.AND.icnlim >= 0) THEN
887 WRITE(*,*) 'Parameter NSTORE to small! Change and rerun!'
888 END IF
889 IF(icnlim >= 0) THEN
890 CALL pxhist(khist,50,0.0,float(nstore)) ! histogram printout
891 END IF
89210 IF(icnlim >= 0) THEN
893 WRITE(*,*) ' '
894 WRITE(*,*) '... making global fit ...'
895 END IF
896 nn=0
897 ii=0 ! modify matrix acccording to PSIGM value
898 DO i=1,nagb
899 ll=ii
900 ii=ii+i
901 IF(psigm(i) == 0.0) THEN
902 nn=nn+1 ! count number of fixed parameters
903 DO j=1,nagb
904 ll=ll+1
905 IF(j > i) ll=ll+j-2
906 cgmat(ll)=0.0 ! reset row and column for parameter
907 END DO
908 ELSE IF(psigm(i) > 0.0) THEN
909 cgmat(ii)=cgmat(ii)+1.0/psigm(i)**2 ! add to diagonal
910 END IF
911 END DO
912 ii=0 ! start saving diagonal elements
913 DO i=1,nagb
914 ii=ii+i
915 diag(i)=cgmat(ii) ! save original diagonal elements
916 END DO
917 nvar=nagb
918 WRITE(*,*) 'Number of constraints ',ncs
919 IF(ncs /= 0) THEN ! add constraints
920 ii=(nagb*nagb+nagb)/2
921 DO i=1,ncs ! loop on constraints
922 dsum=arhs(i)
923 DO j=1,nagb
924 cgmat(ii+j)=adercs(nagb*(i-1)+j)
925 dsum=dsum-adercs(nagb*(i-1)+j)*(pparm(j)+dparm(j))
926 END DO
927 DO j=1,i
928 cgmat(ii+nagb+j)=0.0
929 END DO
930 nvar=nvar+1
931 ii =ii+nvar
932 bgvec(nvar)=dsum
933 END DO
934 END IF
935 ! =================================================
936 CALL spminv(cgmat,bgvec,nvar,nrank,scdiag,scflag)!matrix inversion
937 ndefec=nvar-nrank-nn ! rank defect
938 ! =====================================================
939 DO i=1,nagb
940 dparm(i)=dparm(i)+bgvec(i) ! accumulate corrections
941 END DO
942 CALL prtglo(66)
943 IF(icnlim >= 0) THEN
944 WRITE(*,*) 'The rank defect of the symmetric',nvar,'-by-',nvar, &
945 ' matrix is ',ndefec,' (should be zero).'
946 END IF
947 IF(itert == 0.OR.nstdev == 0.OR.itert >= itelim) GO TO 90
948 ! iterations
949 IF(icnlim >= 0) THEN
950 WRITE(*,*) ' '
951 WRITE(*,*) ' Total',loctot,' local fits,',locrej,' rejected.'
952 WRITE(*,*) ' Histogram of Chisq/Ndf:', &
953 ' Total',nhist,' entries with',mhist(51), ' overflows'
954 CALL pxhist(mhist,50,0.0,10.0) ! histogram printout
955 END IF
956 ! reset histogram
957 nhist=0
958 DO i=1,51
959 mhist(i)=0
960 lhist(i)=0
961 END DO
962 itert=itert+1
963 loctot=0
964 locrej=0
965 IF(cfactr /= 1.0) THEN
966 cfactr=sqrt(cfactr)
967 IF(cfactr < 1.2) THEN
968 cfactr=1.0
969 itelim=itert+1
970 END IF
971 END IF
972 IF(icnlim >= 0) WRITE(*,107) itert,cfactr
973107 FORMAT(' Iteration',i3,' with cut factor=',f6.2)
974 ! reset matrices for global variables
975 DO i=1,nagb
976 bgvec(i)=0.0
977 END DO
978 DO i=1,(nagb*nagb+nagb)/2
979 cgmat(i)=0.0
980 END DO
981 rewind lunit
98220 READ(lunit,END=10) nst,(indst(i),i=1,nst),(arest(i),i=1,nst)
983 CALL fitloc
984 GO TO 20
985 ! ==================================================================
98690 IF(icnlim >= 0) THEN
987 WRITE(*,*) ' '
988 WRITE(*,*) ' Result of fit for global parameters'
989 WRITE(*,*) ' ==================================='
990 WRITE(*,101)
991 END IF
992 ii=0
993 DO i=1,nagb
994 ii=ii+i
995 err=real(sqrt(abs(cgmat(ii))))
996 IF(cgmat((i*i+i)/2) < 0.0) err=-err
997 gcor=0.0
998 IF(cgmat(ii)*diag(i) > 0.0) THEN
999 ! global correlation
1000 gcor=real(sqrt(abs(1.0-1.0/(cgmat(ii)*diag(i)))))
1001 END IF
1002 IF(i <= 25.OR.nagb-i <= 25) THEN
1003 patext=' '
1004 IF(nlnpa(i) /= 0) patext='nl'
1005 IF(icnlim >= 0) WRITE(*,102) i,patext, &
1006 pparm(i),pparm(i)+dparm(i),dparm(i),bgvec(i),err,gcor
1007 END IF
1008 par(i)=real(pparm(i)+dparm(i)) ! copy of result to array in argument
1009 END DO
1010 DO i=1,ncs ! constraints
1011 IF(i == 1.AND.icnlim >= 0) WRITE(*,*) ' '
1012 dsum=0.0
1013 DO j=1,nagb
1014 dsum=dsum+adercs(nagb*(i-1)+j)*(pparm(j)+dparm(j))
1015 END DO
1016 IF(icnlim >= 0) WRITE(*,106) i,dsum,arhs(i),dsum-arhs(i)
1017 END DO
1018 IF(icnlim >= 0) THEN
1019 WRITE(*,*) ' '
1020 WRITE(*,*) ' Total',loctot,' local fits,',locrej,' rejected.'
1021 WRITE(*,*) ' Histogram of RMS:', &
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
1035199 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 ')
1044101 FORMAT(1x,' I initial final differ', &
1045 ' lastcor Error glcor'/ &
1046 1x,' --- ----------- ----------- -----------', &
1047 ' ----------- -------- -----')
1048102 FORMAT(1x,i4,1x,a2,1x,4f12.5,f9.5,f6.3)
1049106 FORMAT(' Constraint',i2,' Sum - RHS =',g12.5,' -',g12.5, ' = ',g12.5)
1050END SUBROUTINE fitglo
1051
1053FUNCTION errpar(i)
1054 ! ------------------------------------------------------------------
1055 ! Basic dimension parameters
1056 parameter(mglobl=1400,mlocal=10,nstore=10000,mcs=10) ! 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, indnz(mglobl),indbk(mglobl)
1074 ! ------------------------------------------------------------------
1075 errpar=0.0
1076 IF(i <= 0.OR.i > nagb) RETURN
1077 ii=(i*i+i)/2
1078 errpar=real(sqrt(abs(cgmat(ii))))
1079 IF(cgmat(ii) < 0.0) errpar=errpar
1080END FUNCTION errpar
1081
1083FUNCTION corpar(i,j)
1084 ! ------------------------------------------------------------------
1085 ! Basic dimension parameters
1086 parameter(mglobl=1400,mlocal=10,nstore=10000,mcs=10) ! 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, indnz(mglobl),indbk(mglobl)
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=real(sqrt(abs(cgmat(ii)*cgmat(jj))))
1114 IF(err /= 0.0) corpar=real(cgmat(ij))/err
1115END FUNCTION corpar
1116
1118SUBROUTINE prtglo(lun)
1119 ! ------------------------------------------------------------------
1120 ! Basic dimension parameters
1121 parameter(mglobl=1400,mlocal=10,nstore=10000,mcs=10) ! dimensions
1122 parameter(mgl=mglobl+mcs) ! derived parameter
1123 ! derived parameters
1124 parameter(msymgb =(mglobl*mglobl+mglobl)/2, msym =(mgl*mgl+mgl)/2, &
1125 msymlc=(mlocal*mlocal+mlocal)/2, mrecta= mglobl*mlocal, &
1126 mglocs= mglobl*mcs, msymcs= (mcs*mcs+mcs)/2 )
1127 DOUBLE PRECISION :: cgmat,clmat,clcmat,bgvec,blvec, &
1128 corrm,corrv,summ,diag,scdiag,pparm,dparm,adercs, 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), dparm(mglobl), &
1134 scdiag(mglobl),summ,scflag(mglobl), &
1135 indgb(mglobl),indlc(mlocal),loctot,locrej, &
1136 nagb,nalc,nsum,nhist,mhist(51),khist(51),lhist(51), &
1137 nst,nfl,indst(nstore),arest(nstore),itert,lunit,ncs, &
1138 nlnpa(mglobl),nstdev,cfactr,icnpr,icnlim, indnz(mglobl),indbk(mglobl)
1139 ! ------------------------------------------------------------------
1140 CHARACTER (LEN=2) :: patext
1141 ! ...
1142 lup=lun
1143 IF(lup == 0) lup=6
1144 WRITE(lup,*) ' Result of fit for global parameters'
1145 WRITE(lup,*) ' ==================================='
1146 WRITE(lup,101)
1147 ii=0
1148 DO i=1,nagb
1149 ii=ii+i
1150 err=real(sqrt(abs(cgmat(ii))))
1151 IF(cgmat(ii) < 0.0) err=-err
1152 gcor=0.0
1153 IF(cgmat(ii)*diag(i) > 0.0) THEN
1154 ! global correlation
1155 gcor=real(sqrt(abs(1.0-1.0/(cgmat(ii)*diag(i)))))
1156 END IF
1157 patext=' '
1158 IF(nlnpa(i) /= 0) patext='nl'
1159 WRITE(lup,102) i,patext, &
1160 pparm(i),pparm(i)+dparm(i),dparm(i),bgvec(i),err,gcor
1161 END DO
1162101 FORMAT(1x,' I initial final differ', &
1163 ' lastcor Error glcor'/ &
1164 1x,' --- ----------- ----------- -----------', &
1165 ' ----------- -------- -----')
1166102 FORMAT(1x,i4,1x,a2,1x,4f12.5,f9.5,f6.3)
1167END SUBROUTINE prtglo
1168
1169
1188SUBROUTINE spminv(v,b,n,nrank,diag,flag)
1189
1190 DOUBLE PRECISION :: v(*),b(n),diag(n),vkk,vjk,eps
1191 LOGICAL :: flag(*)
1192 parameter(eps=1.0d-10)
1193 ! ...
1194 DO i=1,n
1195 flag(i)=.true. ! reset flags
1196 diag(i)=abs(v((i*i+i)/2)) ! save abs of diagonal elements
1197 END DO
1198 nrank=0
1199 DO i=1,n ! start of loop
1200 k =0
1201 jj =0
1202 kk =0
1203 vkk=0.0d0
1204 DO j=1,n ! search for pivot
1205 jj=jj+j
1206 IF(flag(j)) THEN ! not used so far
1207 IF(abs(v(jj)) > max(abs(vkk),eps*diag(j))) THEN
1208 vkk=v(jj) ! pivot (candidate)
1209 k =j ! index of pivot
1210 kk =jj ! index of diagonal element
1211 END IF
1212 END IF
1213 END DO
1214 IF(k /= 0) THEN ! pivot found
1215 nrank=nrank+1 ! increase rank and ...
1216 flag(k)=.false. ! ... reset flag
1217 vkk =1.0/vkk
1218 v(kk) =-vkk
1219 b(k) =b(k)*vkk
1220 jk =kk-k
1221 jl =0
1222 DO j=1,n ! elimination
1223 IF(j == k) THEN
1224 jk=kk
1225 jl=jl+j
1226 ELSE
1227 IF(j < k) THEN
1228 jk=jk+1
1229 ELSE
1230 jk=jk+j-1
1231 END IF
1232 vjk =v(jk)
1233 v(jk)=vkk*vjk
1234 b(j) =b(j)-b(k)*vjk
1235 lk =kk-k
1236 DO l=1,j
1237 jl=jl+1
1238 IF(l == k) THEN
1239 lk=kk
1240 ELSE
1241 IF(l < k) THEN
1242 lk=lk+1
1243 ELSE
1244 lk=lk+l-1
1245 END IF
1246 v(jl)=v(jl)-v(lk)*vjk
1247 END IF
1248 END DO
1249 END IF
1250 END DO
1251 ELSE
1252 DO k=1,n
1253 IF(flag(k)) THEN
1254 b(k)=0.0d0 ! clear vector element
1255 DO j=1,k
1256 IF(flag(j)) v((k*k-k)/2+j)=0.0d0 ! clear matrix row/col
1257 END DO
1258 END IF
1259 END DO
1260 GO TO 10
1261 END IF
1262 END DO ! end of loop
1263 10 DO ij=1,(n*n+n)/2
1264 v(ij)=-v(ij) ! finally reverse sign of all matrix elements
1265 END DO
1266END SUBROUTINE spminv
1267
1282SUBROUTINE spavat(v,a,w,n,m)
1283
1284 DOUBLE PRECISION :: v,a,w,cik
1285 dimension v(*),a(*),w(*)
1286 ! ...
1287 DO i=1,(m*m+m)/2
1288 w(i)=0.0 ! reset output matrix
1289 END DO
1290 il=-n
1291 ijs=0
1292 DO i=1,m ! do I
1293 ijs=ijs+i-1 !
1294 il=il+n !
1295 lkl=0 !
1296 DO k=1,n ! do K
1297 cik=0.0d0 !
1298 lkl=lkl+k-1 !
1299 lk=lkl !
1300 DO l=1,k ! do L
1301 lk=lk+1 ! .
1302 cik=cik+a(il+l)*v(lk) ! .
1303 END DO ! end do L
1304 DO l=k+1,n ! do L
1305 lk=lk+l-1 ! .
1306 cik=cik+a(il+l)*v(lk) ! .
1307 END DO ! end do L
1308 jk=k !
1309 ij=ijs !
1310 DO j=1,i ! do J
1311 ij=ij+1 ! .
1312 w(ij)=w(ij)+cik*a(jk) ! .
1313 jk=jk+n ! .
1314 END DO ! end do J
1315 END DO ! end do K
1316 END DO ! end do I
1317END SUBROUTINE spavat
1318
1328SUBROUTINE spax(a,x,y,m,n)
1329
1330 DOUBLE PRECISION :: a(*),x(*),y(*)
1331 ! ...
1332 ij=0
1333 DO i=1,m
1334 y(i)=0.0d0
1335 DO j=1,n
1336 ij=ij+1
1337 y(i)=y(i)+a(ij)*x(j)
1338 END DO
1339 END DO
1340END SUBROUTINE spax
1341
1343SUBROUTINE pxhist(inc,n,xa,xb)
1344 parameter(maxpl=70)
1345 INTEGER :: inc(n),num(maxpl)
1346 parameter(np=maxpl/10+1)
1347 REAL :: p(np)
1348 equivalence(nval,fval)
1349 CHARACTER (LEN=1) :: text(10)*130 ! 10 rows of X's
1350 CHARACTER (LEN=4) :: echar
1351 CHARACTER (LEN=8) :: xchar(maxpl)
1352 CHARACTER (LEN=1) :: chn(0:10)*1
1353 DATA chn/'0','1','2','3','4','5','6','7','8','9','0'/
1354 ! ...
1355 IF(n < 10) RETURN
1356 ntyp=0
1357 DO i=1,n
1358 IF(inc(i) /= 0) THEN
1359 nval=inc(i)
1360 IF(abs(fval) < 1.0e-30) THEN
1361 ntyp=ntyp+1 ! integer
1362 ELSE
1363 ntyp=ntyp-1 ! floating point
1364 END IF
1365 END IF
1366 END DO
1367 m=1 ! M bins combined to 1
136810 nred=n/m
1369 IF(nred > maxpl) THEN
1370 m=m+m ! M = power of 2
1371 GO TO 10
1372 END IF
1373 nred=n/m ! reduced number of bins
1374 nk=0
1375 IF(ntyp >= 0) THEN ! integer
1376 DO i=1,n,m ! add M bins together
1377 mum=0
1378 DO k=i,min(i+m-1,n)
1379 mum=mum+inc(k)
1380 END DO
1381 nk=nk+1
1382 num(nk)=mum ! copy to array NUM(.)
1383 END DO
1384 ELSE ! floating point
1385 DO i=1,n,m ! add M bins together
1386 sum=0
1387 DO k=i,min(i+m-1,n)
1388 nval=inc(k)
1389 sum=sum+fval
1390 END DO
1391 nk=nk+1
1392 num(nk)=int(sum+0.5) ! copy to array NUM(.)
1393 END DO
1394 END IF
1395 inmax=1 ! find maximum bin
1396 DO i=1,nk
1397 IF(num(i) > num(inmax)) inmax=i
1398 END DO
1399 inh=num(inmax) ! maximum bin content
1400 idiv=1+(inh-1)/10 ! X equivalent
1401 nr=inh/idiv ! number of X lines
1402 DO l=1,nr
1403 text(l)=' ' ! blank text line
1404 END DO
1405 DO k=1,nred
1406 lr=num(k)/idiv
1407 IF(lr /= 0) THEN
1408 DO l=1,lr
1409 text(l)(k:k)='X'
1410 END DO
1411 ELSE
1412 IF(num(k) /= 0) text(1)(i:i)='.'
1413 END IF
1414 END DO
1415 DO l=nr,1,-1 ! print X's
1416 WRITE(*,103) text(l)(1:nred)
1417 END DO
1418 n10=1+(nred-1)/10
1419 text(1)=' '
1420 DO i=1,n10
1421 iup=min(10*i,nred)
1422 text(1)(10*i-9:iup)='----+----+'
1423 IF(iup == 10*i) text(1)(iup:iup)=chn(i)
1424 END DO
1425 WRITE(*,103) text(1)(1:nred)
1426 nap=nred/10+1
1427 IF(nap*10-10 > nred) nap=nap-1
1428 DO l=1,nap
1429 p(l)=xa+float(10*l-10)*(xb-xa)/float(nred)
1430 END DO
1431 CALL chfmt(p,nap,xchar,echar)
1432 WRITE(*,104) (xchar(l),l=1,nap)
1433 IF(echar /= ' ') WRITE(*,105) (echar,l=1,nap)
1434 WRITE(*,103) ' '
1435 nt=0
1436 DO k=1,10
1437 text(k)=' '
1438 DO i=1,nred
1439 IF(num(i) /= 0) THEN
1440 nt=k
1441 l=mod(num(i),10)
1442 num(i)=num(i)/10
1443 text(k)(i:i)=chn(l)
1444 END IF
1445 END DO
1446 END DO
1447 DO l=nt,1,-1
1448 WRITE(*,103) text(l)(1:nred)
1449 END DO
1450 WRITE(*,103) ' '
1451103 FORMAT(5x,a)
1452104 FORMAT(12(a8,2x)/)
1453105 FORMAT(3x,12(a4,6x))
1454END SUBROUTINE pxhist
1455
1472SUBROUTINE chfmt(x,n,xchar,echar)
1473
1474 REAL :: x(*)
1475 CHARACTER (LEN=1) :: ch(10)
1476 CHARACTER (LEN=4) :: echar
1477 CHARACTER (LEN=8) :: xchar(n),fxf,sxf,null
1478 DATA null/'00000000'/,ch/'0','1','2','3','4','5','6','7','8','9'/
1479 ! ...
1480 ! determine factor fc, so that fc*xmax is 5 digit number
1481 ip=0
1482 xm=0.0
1483 DO i=1,n
1484 xm=amax1(xm,abs(x(i)))
1485 END DO
1486 IF(xm /= 0.0) THEN
1487 jp=104-ifix(alog10(abs(xm))+100.04)
1488 fc=10.0**jp
1489 ELSE
1490 jp=5
1491 fc=1.0
1492 END IF
1493 ! store digits as characters and find jm = first nonzero digit
1494 im=6
1495 DO i=1,n
1496 fxf=null
1497 ij=int(fc*abs(x(i))+0.5)
1498 jm=6
1499 IF(ij /= 0) THEN
1500 DO j=1,5
1501 jn=mod(ij,10)
1502 ij=ij/10
1503 IF(jn /= 0.AND.jm == 6) jm=j
1504 fxf(j:j)=ch(jn+1)
1505 END DO
1506 im=min(im,jm)
1507 END IF
1508 xchar(i)=fxf
1509 END DO
1510 jm=im
1511 ! determine exponent as a multiple of 3
151232 IF(jp < 1) THEN
1513 jp=jp+3
1514 ip=ip+3
1515 GO TO 32
1516 END IF
151734 IF(jp > jm+4.OR.jp >= 8) THEN
1518 jp=jp-3
1519 ip=ip-3
1520 GO TO 34
1521 END IF
1522 ! loop to convert to print format
1523 ja=min(jm,jp)
1524 jb=max(6,jp+1)
1525 DO i=1,n
1526 fxf=xchar(i)
1527 sxf=' '
1528 ib=7+(jb-ja)/2
1529 DO j=ja,jb
1530 IF(fxf(j:j) /= ch(1)) GO TO 70
1531 IF(j > jp+1) GO TO 50
1532 IF(fxf /= null.OR.j >= jp) GO TO 70
1533 ib=ib-1
1534 cycle
1535 50 DO k=j,jb
1536 IF(fxf(k:k) /= ch(1)) GO TO 70
1537 END DO
1538 cycle
1539 ! insert digit
154070 ib=ib-1
1541 sxf(ib:ib)=fxf(j:j)
1542 IF(j == jp) THEN
1543 ! insert decimal dot
1544 ib=ib-1
1545 sxf(ib:ib)='.'
1546 END IF
1547 END DO
1548 ! insert - sign
1549 IF(x(i) < 0.0) sxf(ib-1:ib-1)='-'
1550 xchar(i)=sxf
1551 END DO
1552 ! prepare print format for exponent
1553 echar=' '
1554 IF(ip /= 0) THEN
1555 echar='E 0 '
1556 IF(ip <= 0) THEN
1557 echar(2:2)='-'
1558 ip=iabs(ip)
1559 END IF
1560 j=mod(ip,10)
1561 echar(4:4)=ch(j+1)
1562 ip=(ip-j)/10
1563 IF(ip /= 0) THEN
1564 j=mod(ip,10)
1565 echar(3:3)=ch(j+1)
1566 END IF
1567 END IF
1568END SUBROUTINE chfmt
1569
1571FUNCTION chindl(n,nd)
1572 REAL :: pn(3),sn(3),table(30,3)
1573 DATA pn/0.31731,0.0455002785,2.69985e-3/ ! probabilities
1574 DATA sn/0.47523,1.690140,2.782170/
1575 DATA table/ 1.0000, 1.1479, 1.1753, 1.1798, 1.1775, 1.1730, 1.1680, 1.1630, &
1576 1.1581, 1.1536, 1.1493, 1.1454, 1.1417, 1.1383, 1.1351, 1.1321, &
1577 1.1293, 1.1266, 1.1242, 1.1218, 1.1196, 1.1175, 1.1155, 1.1136, &
1578 1.1119, 1.1101, 1.1085, 1.1070, 1.1055, 1.1040, &
1579 4.0000, 3.0900, 2.6750, 2.4290, 2.2628, 2.1415, 2.0481, 1.9736, &
1580 1.9124, 1.8610, 1.8171, 1.7791, 1.7457, 1.7161, 1.6897, 1.6658, &
1581 1.6442, 1.6246, 1.6065, 1.5899, 1.5745, 1.5603, 1.5470, 1.5346, &
1582 1.5230, 1.5120, 1.5017, 1.4920, 1.4829, 1.4742, &
1583 9.0000, 5.9146, 4.7184, 4.0628, 3.6410, 3.3436, 3.1209, 2.9468, &
1584 2.8063, 2.6902, 2.5922, 2.5082, 2.4352, 2.3711, 2.3143, 2.2635, &
1585 2.2178, 2.1764, 2.1386, 2.1040, 2.0722, 2.0428, 2.0155, 1.9901, &
1586 1.9665, 1.9443, 1.9235, 1.9040, 1.8855, 1.8681/
1587 ! ...
1588 IF(nd < 1) THEN
1589 chindl=0.0
1590 ELSE
1591 m=max(1,min(n,3)) ! 1, 2 or 3 sigmas
1592 IF(nd <= 30) THEN
1593 chindl=table(nd,m) ! from table
1594 ELSE ! approximation for ND > 30
1595 chindl=(sn(m)+sqrt(float(nd+nd-1)))**2/float(nd+nd)
1596 END IF
1597 END IF
1598END FUNCTION chindl
1599
1601SUBROUTINE fitmut(nvec,vec)
1602 ! ------------------------------------------------------------------
1603 ! Basic dimension parameters
1604 parameter(mglobl=1400,mlocal=10,nstore=10000,mcs=10) ! dimensions
1605 parameter(mgl=mglobl+mcs) ! derived parameter
1606 ! derived parameters
1607 parameter(msymgb =(mglobl*mglobl+mglobl)/2, msym =(mgl*mgl+mgl)/2, &
1608 msymlc=(mlocal*mlocal+mlocal)/2, mrecta= mglobl*mlocal, &
1609 mglocs= mglobl*mcs, msymcs= (mcs*mcs+mcs)/2 )
1610 DOUBLE PRECISION :: cgmat,clmat,clcmat,bgvec,blvec, &
1611 corrm,corrv,summ,diag,scdiag,pparm,dparm,adercs, arhs
1612 LOGICAL :: scflag
1613 common/lsqred/cgmat(msym),clmat(msymlc),clcmat(mrecta), &
1614 diag(mgl),bgvec(mgl),blvec(mlocal), &
1615 corrm(msymgb),corrv(mglobl),psigm(mglobl), &
1616 pparm(mglobl),adercs(mglocs),arhs(mcs), dparm(mglobl), &
1617 scdiag(mglobl),summ,scflag(mglobl), &
1618 indgb(mglobl),indlc(mlocal),loctot,locrej, &
1619 nagb,nalc,nsum,nhist,mhist(51),khist(51),lhist(51), &
1620 nst,nfl,indst(nstore),arest(nstore),itert,lunit,ncs, &
1621 nlnpa(mglobl),nstdev,cfactr,icnpr,icnlim, indnz(mglobl),indbk(mglobl)
1622 ! ------------------------------------------------------------------
1623 REAL :: vec(*)
1624 ! ...
1625 nvec=(nagb*nagb+nagb)/2+nagb
1626 DO i=1,nvec-nagb
1627 vec(i)=real(cgmat(i))
1628 END DO
1629 DO i=1,nagb
1630 vec(nvec-nagb+i)=real(bgvec(i))
1631 END DO
1632 WRITE(*,*) 'FITMUT ',nvec,(vec(i),i=1,nvec)
1633 RETURN
1634 entry fitmin(nvec,vec)
1635 ! insert information
1636 IF(nvec /= (nagb*nagb+nagb)/2+nagb) THEN
1637 WRITE(*,*) ' Wrong dimensions in FITMUT/FITMIN'
1638 WRITE(*,*) ' Argument NVEC =',nvec
1639 WRITE(*,*) ' Expected NVEC =',(nagb*nagb+nagb)/2+nagb
1640 stop 'FITMIN'
1641 END IF
1642 DO i=1,nvec-nagb
1643 cgmat(i)=vec(i)
1644 END DO
1645 DO i=1,nagb
1646 bgvec(i)=vec(nvec-nagb+i)
1647 END DO
1648END SUBROUTINE fitmut
1649
1650
1651
1652
1653
1654
1655
1656
1657
1658
1659
znorm
function znorm()
Return random number N(0,1).
Definition: millepede1.f90:154
initgl
subroutine initgl(nagbar, nalcar, nstd, iprlim)
Initialization of package.
Definition: millepede1.f90:171
errpar
function errpar(i)
Return error for parameter I.
Definition: millepede1.f90:1054
initun
subroutine initun(lun, cutfac)
Define unit for iterations (optional).
Definition: millepede1.f90:355
equloc
subroutine equloc(dergb, derlc, rrmeas, sigma)
Add single equation with its derivatives.
Definition: millepede1.f90:437
fitmut
subroutine fitmut(nvec, vec)
Get matrix information out.
Definition: millepede1.f90:1602
fitglo
subroutine fitglo(par)
Final global fit.
Definition: millepede1.f90:841
pxhist
subroutine pxhist(inc, n, xa, xb)
Print X histogram.
Definition: millepede1.f90:1344
constf
subroutine constf(dercs, rhs)
Add constraint (optional).
Definition: millepede1.f90:398
chfmt
subroutine chfmt(x, n, xchar, echar)
Prepare printout of array of real numbers as character strings.
Definition: millepede1.f90:1473
testmp
subroutine testmp(iarg)
Test of millepede.
Definition: millepede1.f90:69
gener
subroutine gener(y, a, b, x, sigma, bias, heff)
Generate Y values.
Definition: millepede1.f90:118
spavat
subroutine spavat(v, a, w, n, m)
Similarity operation A*V*A^t.
Definition: millepede1.f90:1283
corpar
function corpar(i, j)
Return correlation between parameters I and J.
Definition: millepede1.f90:1084
spax
subroutine spax(a, x, y, m, n)
Multiply general M-by-N matrix A and N-vector X.
Definition: millepede1.f90:1329
zrand
function zrand()
Return random number U(0,1).
Definition: millepede1.f90:144
fitloc
subroutine fitloc
Fit after end of local block.
Definition: millepede1.f90:555
zerloc
subroutine zerloc(dergb, derlc)
Reset derivatives.
Definition: millepede1.f90:523
nonlin
subroutine nonlin(INDEX)
Set nonlinear flag for single parameter (optional).
Definition: millepede1.f90:328
parglo
subroutine parglo(par)
Initialize global parameters.
Definition: millepede1.f90:271
parsig
subroutine parsig(INDEX, sigma)
Define sigma for single parameter (optional).
Definition: millepede1.f90:300
spminv
subroutine spminv(v, b, n, nrank, diag, flag)
Obtain solution of a system of linear equations with symmetric matrix and the inverse.
Definition: millepede1.f90:1189
chindl
real(mps) function chindl(n, nd)
Chi2/ndf cuts.
Definition: mpnum.f90:2079
prtglo
subroutine prtglo
Print final log file.
Definition: pede.f90:5417
Generated on Wed Sep 18 2024 14:36:52 for Millepede-II by doxygen 1.9.3