50SUBROUTINE initgl(nagbar,nalcar,nstd,iprlim)
54 INTEGER,
INTENT(IN) :: nagbar
55 INTEGER,
INTENT(IN) :: nalcar
56 INTEGER,
INTENT(IN OUT) :: nstd
57 INTEGER,
INTENT(IN) :: iprlim
58 INTEGER,
PARAMETER :: mglobl=1400
59 INTEGER,
PARAMETER :: mlocal=10
60 INTEGER,
PARAMETER :: nstore=10000
61 INTEGER,
PARAMETER :: mcs=10
62 INTEGER,
PARAMETER :: mgl=mglobl+mcs
64 INTEGER,
PARAMETER :: msymgb =(mglobl*mglobl+mglobl)/2
65 INTEGER,
PARAMETER :: msym =(mgl*mgl+mgl)/2
66 INTEGER,
PARAMETER :: msymlc=(mlocal*mlocal+mlocal)/2
67 INTEGER,
PARAMETER :: mrecta= mglobl*mlocal
68 INTEGER,
PARAMETER :: mglocs= mglobl*mcs
69 INTEGER,
PARAMETER :: msymcs= (mcs*mcs+mcs)/2
70 DOUBLE PRECISION :: cgmat,clmat,clcmat,bgvec,blvec, &
71 corrm,corrv,summ,diag,scdiag,pparm,dparm,adercs, arhs
73 common/lsqred/cgmat(msym),clmat(msymlc),clcmat(mrecta), &
74 diag(mgl),bgvec(mgl),blvec(mlocal), &
75 corrm(msymgb),corrv(mglobl),psigm(mglobl), &
76 pparm(mglobl),adercs(mglocs),arhs(mcs), dparm(mglobl), &
77 scdiag(mglobl),summ,scflag(mglobl), &
78 indgb(mglobl),indlc(mlocal),loctot,locrej, &
79 nagb,nalc,nsum,nhist,mhist(51),khist(51),lhist(51), &
80 nst,nfl,indst(nstore),arest(nstore),itert,lunit,ncs, &
81 nlnpa(mglobl),nstdev,cfactr,icnpr,icnlim
84 DATA ndr/1,2,5,10,20,50,100/
87 IF(icnlim >= 0)
WRITE(*,199)
88199
FORMAT(
' MP-I to MP-II interface: create steering, binary files'/ &
92 ' o ooooo o o o oo ooo oo ooo oo '/ &
93 ' o o o o o o o o o o o o o o o o '/ &
94 ' o o o o o o oooo o o oooo o o oooo '/ &
95 ' o o o o o o o ooo o o o o '/ &
96 ' o o o o oo oo oo o oo ooo oo starting'/ &
104 WRITE(*,*)
'Number of global parameters ',nagb
105 WRITE(*,*)
'Number of local parameters ',nalc
109 WRITE(*,*)
'MP-II uses nstd=3 istead of ',nstd
113 IF(nagb > mglobl.OR.nalc > mlocal)
THEN
114 WRITE(*,*)
'Too many parameter - STOP'
123 OPEN(unit=lunit,access=
'SEQUENTIAL',form=
'UNFORMATTED', file=
'mp2tst.bin')
137 REAL,
INTENT(IN) :: par(*)
138 INTEGER,
PARAMETER :: mglobl=1400
139 INTEGER,
PARAMETER :: mlocal=10
140 INTEGER,
PARAMETER :: nstore=10000
141 INTEGER,
PARAMETER :: mcs=10
142 INTEGER,
PARAMETER :: mgl=mglobl+mcs
144 INTEGER,
PARAMETER :: msymgb =(mglobl*mglobl+mglobl)/2
145 INTEGER,
PARAMETER :: msym =(mgl*mgl+mgl)/2
146 INTEGER,
PARAMETER :: msymlc=(mlocal*mlocal+mlocal)/2
147 INTEGER,
PARAMETER :: mrecta= mglobl*mlocal
148 INTEGER,
PARAMETER :: mglocs= mglobl*mcs
149 INTEGER,
PARAMETER :: msymcs= (mcs*mcs+mcs)/2
150 DOUBLE PRECISION :: cgmat,clmat,clcmat,bgvec,blvec, &
151 corrm,corrv,summ,diag,scdiag,pparm,dparm,adercs, arhs
153 common/lsqred/cgmat(msym),clmat(msymlc),clcmat(mrecta), &
154 diag(mgl),bgvec(mgl),blvec(mlocal), &
155 corrm(msymgb),corrv(mglobl),psigm(mglobl), &
156 pparm(mglobl),adercs(mglocs),arhs(mcs), dparm(mglobl), &
157 scdiag(mglobl),summ,scflag(mglobl), &
158 indgb(mglobl),indlc(mlocal),loctot,locrej, &
159 nagb,nalc,nsum,nhist,mhist(51),khist(51),lhist(51), &
160 nst,nfl,indst(nstore),arest(nstore),itert,lunit,ncs, &
161 nlnpa(mglobl),nstdev,cfactr,icnpr,icnlim
175 INTEGER,
INTENT(IN) :: INDEX
176 REAL,
INTENT(IN) :: sigma
177 INTEGER,
PARAMETER :: mglobl=1400
178 INTEGER,
PARAMETER :: mlocal=10
179 INTEGER,
PARAMETER :: nstore=10000
180 INTEGER,
PARAMETER :: mcs=10
181 INTEGER,
PARAMETER :: mgl=mglobl+mcs
183 INTEGER,
PARAMETER :: msymgb =(mglobl*mglobl+mglobl)/2
184 INTEGER,
PARAMETER :: msym =(mgl*mgl+mgl)/2
185 INTEGER,
PARAMETER :: msymlc=(mlocal*mlocal+mlocal)/2
186 INTEGER,
PARAMETER :: mrecta= mglobl*mlocal
187 INTEGER,
PARAMETER :: mglocs= mglobl*mcs
188 INTEGER,
PARAMETER :: msymcs= (mcs*mcs+mcs)/2
189 DOUBLE PRECISION :: cgmat,clmat,clcmat,bgvec,blvec, &
190 corrm,corrv,summ,diag,scdiag,pparm,dparm,adercs, arhs
192 common/lsqred/cgmat(msym),clmat(msymlc),clcmat(mrecta), &
193 diag(mgl),bgvec(mgl),blvec(mlocal), &
194 corrm(msymgb),corrv(mglobl),psigm(mglobl), &
195 pparm(mglobl),adercs(mglocs),arhs(mcs), dparm(mglobl), &
196 scdiag(mglobl),summ,scflag(mglobl), &
197 indgb(mglobl),indlc(mlocal),loctot,locrej, &
198 nagb,nalc,nsum,nhist,mhist(51),khist(51),lhist(51), &
199 nst,nfl,indst(nstore),arest(nstore),itert,lunit,ncs, &
200 nlnpa(mglobl),nstdev,cfactr,icnpr,icnlim
202 IF(index < 1.OR.index > nagb)
RETURN
203 IF(sigma < 0.0)
RETURN
212 INTEGER,
INTENT(IN) :: INDEX
213 INTEGER,
PARAMETER :: mglobl=1400
214 INTEGER,
PARAMETER :: mlocal=10
215 INTEGER,
PARAMETER :: nstore=10000
216 INTEGER,
PARAMETER :: mcs=10
217 INTEGER,
PARAMETER :: mgl=mglobl+mcs
219 INTEGER,
PARAMETER :: msymgb =(mglobl*mglobl+mglobl)/2
220 INTEGER,
PARAMETER :: msym =(mgl*mgl+mgl)/2
221 INTEGER,
PARAMETER :: msymlc=(mlocal*mlocal+mlocal)/2
222 INTEGER,
PARAMETER :: mrecta= mglobl*mlocal
223 INTEGER,
PARAMETER :: mglocs= mglobl*mcs
224 INTEGER,
PARAMETER :: msymcs= (mcs*mcs+mcs)/2
225 DOUBLE PRECISION :: cgmat,clmat,clcmat,bgvec,blvec, &
226 corrm,corrv,summ,diag,scdiag,pparm,dparm,adercs, arhs
228 common/lsqred/cgmat(msym),clmat(msymlc),clcmat(mrecta), &
229 diag(mgl),bgvec(mgl),blvec(mlocal), &
230 corrm(msymgb),corrv(mglobl),psigm(mglobl), &
231 pparm(mglobl),adercs(mglocs),arhs(mcs), dparm(mglobl), &
232 scdiag(mglobl),summ,scflag(mglobl), &
233 indgb(mglobl),indlc(mlocal),loctot,locrej, &
234 nagb,nalc,nsum,nhist,mhist(51),khist(51),lhist(51), &
235 nst,nfl,indst(nstore),arest(nstore),itert,lunit,ncs, &
236 nlnpa(mglobl),nstdev,cfactr,icnpr,icnlim
238 IF(index < 1.OR.index > nagb)
RETURN
244 WRITE(*,*)
' INITUN is dummy !', lun, cutfac
252 REAL,
INTENT(IN) :: dercs(*)
253 REAL,
INTENT(IN) :: rhs
254 INTEGER,
PARAMETER :: mglobl=1400
255 INTEGER,
PARAMETER :: mlocal=10
256 INTEGER,
PARAMETER :: nstore=10000
257 INTEGER,
PARAMETER :: mcs=10
258 INTEGER,
PARAMETER :: mgl=mglobl+mcs
260 INTEGER,
PARAMETER :: msymgb =(mglobl*mglobl+mglobl)/2
261 INTEGER,
PARAMETER :: msym =(mgl*mgl+mgl)/2
262 INTEGER,
PARAMETER :: msymlc=(mlocal*mlocal+mlocal)/2
263 INTEGER,
PARAMETER :: mrecta= mglobl*mlocal
264 INTEGER,
PARAMETER :: mglocs= mglobl*mcs
265 INTEGER,
PARAMETER :: msymcs= (mcs*mcs+mcs)/2
266 DOUBLE PRECISION :: cgmat,clmat,clcmat,bgvec,blvec, &
267 corrm,corrv,summ,diag,scdiag,pparm,dparm,adercs, arhs
269 common/lsqred/cgmat(msym),clmat(msymlc),clcmat(mrecta), &
270 diag(mgl),bgvec(mgl),blvec(mlocal), &
271 corrm(msymgb),corrv(mglobl),psigm(mglobl), &
272 pparm(mglobl),adercs(mglocs),arhs(mcs), dparm(mglobl), &
273 scdiag(mglobl),summ,scflag(mglobl), &
274 indgb(mglobl),indlc(mlocal),loctot,locrej, &
275 nagb,nalc,nsum,nhist,mhist(51),khist(51),lhist(51), &
276 nst,nfl,indst(nstore),arest(nstore),itert,lunit,ncs, &
277 nlnpa(mglobl),nstdev,cfactr,icnpr,icnlim
281 IF(ncs >= mcs) stop
'<INITCS> too many constraints'
283 adercs(nagb*ncs+i)=dercs(i)
296SUBROUTINE equloc(dergb,derlc,rrmeas,sigma)
298 REAL,
INTENT(OUT) :: dergb(*)
299 REAL,
INTENT(OUT) :: derlc(*)
300 REAL,
INTENT(IN) :: rrmeas
301 REAL,
INTENT(IN) :: sigma
304 INTEGER,
PARAMETER :: mglobl=1400
305 INTEGER,
PARAMETER :: mlocal=10
306 INTEGER,
PARAMETER :: nstore=10000
307 INTEGER,
PARAMETER :: mcs=10
308 INTEGER,
PARAMETER :: mgl=mglobl+mcs
310 INTEGER,
PARAMETER :: msymgb =(mglobl*mglobl+mglobl)/2
311 INTEGER,
PARAMETER :: msym =(mgl*mgl+mgl)/2
312 INTEGER,
PARAMETER :: msymlc=(mlocal*mlocal+mlocal)/2
313 INTEGER,
PARAMETER :: mrecta= mglobl*mlocal
314 INTEGER,
PARAMETER :: mglocs= mglobl*mcs
315 INTEGER,
PARAMETER :: msymcs= (mcs*mcs+mcs)/2
316 DOUBLE PRECISION :: cgmat,clmat,clcmat,bgvec,blvec, &
317 corrm,corrv,summ,diag,scdiag,pparm,dparm,adercs, arhs
319 common/lsqred/cgmat(msym),clmat(msymlc),clcmat(mrecta), &
320 diag(mgl),bgvec(mgl),blvec(mlocal), &
321 corrm(msymgb),corrv(mglobl),psigm(mglobl), &
322 pparm(mglobl),adercs(mglocs),arhs(mcs), dparm(mglobl), &
323 scdiag(mglobl),summ,scflag(mglobl), &
324 indgb(mglobl),indlc(mlocal),loctot,locrej, &
325 nagb,nalc,nsum,nhist,mhist(51),khist(51),lhist(51), &
326 nst,nfl,indst(nstore),arest(nstore),itert,lunit,ncs, &
327 nlnpa(mglobl),nstdev,cfactr,icnpr,icnlim
332 IF(sigma <= 0.0)
THEN
345 IF(derlc(i) /= 0.0)
THEN
354 IF(dergb(i) /= 0.0)
THEN
360 IF(nst+nonzer+2 >= nstore)
THEN
368 IF(derlc(i) /= 0.0)
THEN
379 IF(dergb(i) /= 0.0)
THEN
394 REAL,
INTENT(OUT) :: dergb(*)
395 REAL,
INTENT(OUT) :: derlc(*)
398 INTEGER,
PARAMETER :: mglobl=1400
399 INTEGER,
PARAMETER :: mlocal=10
400 INTEGER,
PARAMETER :: nstore=10000
401 INTEGER,
PARAMETER :: mcs=10
402 INTEGER,
PARAMETER :: mgl=mglobl+mcs
404 INTEGER,
PARAMETER :: msymgb =(mglobl*mglobl+mglobl)/2
405 INTEGER,
PARAMETER :: msym =(mgl*mgl+mgl)/2
406 INTEGER,
PARAMETER :: msymlc=(mlocal*mlocal+mlocal)/2
407 INTEGER,
PARAMETER :: mrecta= mglobl*mlocal
408 INTEGER,
PARAMETER :: mglocs= mglobl*mcs
409 INTEGER,
PARAMETER :: msymcs= (mcs*mcs+mcs)/2
410 DOUBLE PRECISION :: cgmat,clmat,clcmat,bgvec,blvec, &
411 corrm,corrv,summ,diag,scdiag,pparm,dparm,adercs, arhs
413 common/lsqred/cgmat(msym),clmat(msymlc),clcmat(mrecta), &
414 diag(mgl),bgvec(mgl),blvec(mlocal), &
415 corrm(msymgb),corrv(mglobl),psigm(mglobl), &
416 pparm(mglobl),adercs(mglocs),arhs(mcs), dparm(mglobl), &
417 scdiag(mglobl),summ,scflag(mglobl), &
418 indgb(mglobl),indlc(mlocal),loctot,locrej, &
419 nagb,nalc,nsum,nhist,mhist(51),khist(51),lhist(51), &
420 nst,nfl,indst(nstore),arest(nstore),itert,lunit,ncs, &
421 nlnpa(mglobl),nstdev,cfactr,icnpr,icnlim
438 INTEGER,
PARAMETER :: mglobl=1400
439 INTEGER,
PARAMETER :: mlocal=10
440 INTEGER,
PARAMETER :: nstore=10000
441 INTEGER,
PARAMETER :: mcs=10
442 INTEGER,
PARAMETER :: mgl=mglobl+mcs
444 INTEGER,
PARAMETER :: msymgb =(mglobl*mglobl+mglobl)/2
445 INTEGER,
PARAMETER :: msym =(mgl*mgl+mgl)/2
446 INTEGER,
PARAMETER :: msymlc=(mlocal*mlocal+mlocal)/2
447 INTEGER,
PARAMETER :: mrecta= mglobl*mlocal
448 INTEGER,
PARAMETER :: mglocs= mglobl*mcs
449 INTEGER,
PARAMETER :: msymcs= (mcs*mcs+mcs)/2
450 DOUBLE PRECISION :: cgmat,clmat,clcmat,bgvec,blvec, &
451 corrm,corrv,summ,diag,scdiag,pparm,dparm,adercs, arhs
453 common/lsqred/cgmat(msym),clmat(msymlc),clcmat(mrecta), &
454 diag(mgl),bgvec(mgl),blvec(mlocal), &
455 corrm(msymgb),corrv(mglobl),psigm(mglobl), &
456 pparm(mglobl),adercs(mglocs),arhs(mcs), dparm(mglobl), &
457 scdiag(mglobl),summ,scflag(mglobl), &
458 indgb(mglobl),indlc(mlocal),loctot,locrej, &
459 nagb,nalc,nsum,nhist,mhist(51),khist(51),lhist(51), &
460 nst,nfl,indst(nstore),arest(nstore),itert,lunit,ncs, &
461 nlnpa(mglobl),nstdev,cfactr,icnpr,icnlim
466 WRITE(lunit) nst*2,(arest(i),i=1,nst),(indst(i),i=1,nst)
473 ibin=int(1.0+50.0*float(nst)/float(nstore))
475 khist(ibin)=khist(ibin)+1
477 khist(51)=khist(51)+1
487 REAL,
INTENT(IN OUT) :: par(*)
491 INTEGER,
PARAMETER :: mglobl=1400
492 INTEGER,
PARAMETER :: mlocal=10
493 INTEGER,
PARAMETER :: nstore=10000
494 INTEGER,
PARAMETER :: mcs=10
495 INTEGER,
PARAMETER :: mgl=mglobl+mcs
497 INTEGER,
PARAMETER :: msymgb =(mglobl*mglobl+mglobl)/2
498 INTEGER,
PARAMETER :: msym =(mgl*mgl+mgl)/2
499 INTEGER,
PARAMETER :: msymlc=(mlocal*mlocal+mlocal)/2
500 INTEGER,
PARAMETER :: mrecta= mglobl*mlocal
501 INTEGER,
PARAMETER :: mglocs= mglobl*mcs
502 INTEGER,
PARAMETER :: msymcs= (mcs*mcs+mcs)/2
503 DOUBLE PRECISION :: cgmat,clmat,clcmat,bgvec,blvec, &
504 corrm,corrv,summ,diag,scdiag,pparm,dparm,adercs, arhs
506 common/lsqred/cgmat(msym),clmat(msymlc),clcmat(mrecta), &
507 diag(mgl),bgvec(mgl),blvec(mlocal), &
508 corrm(msymgb),corrv(mglobl),psigm(mglobl), &
509 pparm(mglobl),adercs(mglocs),arhs(mcs), dparm(mglobl), &
510 scdiag(mglobl),summ,scflag(mglobl), &
511 indgb(mglobl),indlc(mlocal),loctot,locrej, &
512 nagb,nalc,nsum,nhist,mhist(51),khist(51),lhist(51), &
513 nst,nfl,indst(nstore),arest(nstore),itert,lunit,ncs, &
514 nlnpa(mglobl),nstdev,cfactr,icnpr,icnlim
521 IF (pparm(i) /= 0.0.OR.psigm(i) >= 0.0) npar=npar+1
526 OPEN(unit=luns,access=
'SEQUENTIAL',form=
'FORMATTED', file=
'mp2str.txt')
527 WRITE(luns,101)
'* Default test steering file'
528 WRITE(luns,101)
'fortranfiles ! following bin files are fortran'
529 WRITE(luns,101)
'mp2tst.bin ! binary data file'
532 WRITE(luns,101)
'mp2con.txt ! constraints text file '
536 WRITE(luns,101)
'Parameter ! with start values or pre-sigma'
538 IF (pparm(i) /= 0.0.OR.psigm(i) >= 0.0)
THEN
541 IF (presig == 0.0) presig=-1.0
542 WRITE(luns,103) i, pparm(i), presig
547 WRITE(luns,101)
'chisqcut 30.0 6.0 ! Chi2/ndf cut factors'
548 WRITE(luns,101)
'method inversion 3 0.001 ! Gauss matrix inversion'
550 WRITE(luns,101)
'end ! optional for end-of-data'
555 OPEN(unit=lunc,access=
'SEQUENTIAL',form=
'FORMATTED', file=
'mp2con.txt')
557 WRITE(lunc,101)
'Constraint 0.0'
559 IF (adercs(nagb*i+j) /= 0.0)
WRITE(lunc,102) j, adercs(nagb*i+j)
568 ' o ooooo o o o oo ooo oo ooo oo '/ &
569 ' o o o o o o o o o o o o o o o o '/ &
570 ' o o o o o o oooo o o oooo o o oooo '/ &
571 ' o o o o o o o ooo o o o o '/ &
572 ' o o o o oo oo oo o oo ooo oo ending.'/ &
575 ' MP-I to MP-II interface: create steering, binary files')
585 WRITE(*,*)
' ERRPAR is dummy !', i
591 WRITE(*,*)
' CORPAR is dummy !', i, j
596 WRITE(*,*)
' PRTGLO is dummy !', lun
subroutine initgl(nagbar, nalcar, nstd, iprlim)
Initialization of package.
function errpar(i)
Return error for parameter I.
subroutine initun(lun, cutfac)
Define unit for iterations (optional).
subroutine equloc(dergb, derlc, rrmeas, sigma)
Add single equation with its derivatives.
subroutine fitglo(par)
Final global fit.
subroutine constf(dercs, rhs)
Add constraint (optional).
function corpar(i, j)
Return correlation between parameters I and J.
subroutine fitloc
Fit after end of local block.
subroutine zerloc(dergb, derlc)
Reset derivatives.
subroutine nonlin(INDEX)
Set nonlinear flag for single parameter (optional).
subroutine parglo(par)
Initialize global parameters.
subroutine parsig(INDEX, sigma)
Define sigma for single parameter (optional).
subroutine prtglo
Print final log file.