574 REAL,
DIMENSION(2) :: ta
587 INTEGER(mpi) :: nsecnd
589 INTEGER(mpi) :: ntsec
591 CHARACTER (LEN=24) :: chdate
592 CHARACTER (LEN=24) :: chost
613 CALL getenv(
'HOSTNAME',chost)
614 IF (chost(1:1) ==
' ')
CALL getenv(
'HOST',chost)
615 WRITE(*,*)
'($Rev: 169 $)'
618 WRITE(*,111) __gnuc__ , __gnuc_minor__ , __gnuc_patchlevel__
619111
FORMAT(
' compiled with gcc ',i0,
'.',i0,
'.',i0)
622 WRITE(*,*)
' < Millepede II-P starting ... ',chdate
627 WRITE(8,*)
'Log-file Millepede II-P ', chdate
628 WRITE(8,*)
' ', chost
629 CALL peend(-1,
'Still running or crashed')
635 WRITE(*,*)
'!!! Checking input only, no calculation of a solution !!!'
636 WRITE(8,*)
'!!! Checking input only, no calculation of a solution !!!'
660 CALL mvopen(lun,
'millepede.his')
666 CALL mvopen(1,
'mpdebug.txt')
669 CALL etime(ta,rstext)
675 CALL etime(ta,rloop1)
676 times(1)=rloop1-rstext
680 WRITE(8,*)
'Chi square cut equiv 3 st.dev applied ...'
681 WRITE(8,*)
' in first iteration with factor',
chicut
682 WRITE(8,*)
' in second iteration with factor',
chirem
683 WRITE(8,*)
' (reduced by sqrt in next iterations)'
687 WRITE(8,*)
'Down-weighting of outliers in',
lhuber,
' iterations'
688 WRITE(8,*)
'Cut on downweight fraction',
dwcut
691 CALL etime(ta,rloop2)
692 times(2)=rloop2-rloop1
696 CALL peend(0,
'Ended normally')
732 CALL gmpdef(6,1,
'log10(#records) vs file number')
733 CALL gmpdef(7,1,
'final rejection fraction vs file number')
735 'final <Chi^2/Ndf> from accepted local fits vs file number')
736 CALL gmpdef(9,1,
'<Ndf> from accepted local fits vs file number')
742 rej=real(nrc-
jfd(kfl),mps)/real(nrc,mps)
743 CALL gmpxy(6,real(kfl,mps),log10(real(nrc,mps)))
744 CALL gmpxy(7,real(kfl,mps),rej)
746 IF (
jfd(kfl) > 0)
THEN
747 c2ndf=
cfd(kfl)/real(
jfd(kfl),mps)
748 CALL gmpxy(8,real(kfl,mps),c2ndf)
749 andf=real(
dfd(kfl),mps)/real(
jfd(kfl),mps)
750 CALL gmpxy(9,real(kfl,mps),andf)
767 WRITE(*,*)
'Misalignment test wire chamber'
770 CALL hmpdef( 9,-0.0015,+0.0015,
'True - fitted displacement')
771 CALL hmpdef(10,-0.0015,+0.0015,
'True - fitted Vdrift')
778 sums(2)=sums(2)+diff*diff
781 sums(4)=sums(4)+diff*diff
783 sums(1)=0.01_mpd*sums(1)
784 sums(2)=sqrt(0.01_mpd*sums(2))
785 sums(3)=0.01_mpd*sums(3)
786 sums(4)=sqrt(0.01_mpd*sums(4))
787 WRITE(*,143)
'Parameters 1 - 100: mean =',sums(1),
'rms =',sums(2)
788 WRITE(*,143)
'Parameters 101 - 200: mean =',sums(3),
'rms =',sums(4)
789143
FORMAT(6x,a28,f9.6,3x,a5,f9.6)
811 WRITE(*,*)
'Misalignment test Si tracker'
814 CALL hmpdef( 9,-0.0025,+0.0025,
'True - fitted displacement X')
815 CALL hmpdef(10,-0.025,+0.025,
'True - fitted displacement Y')
826 sums(1)=sums(1)+1.0_mpd
828 sums(3)=sums(3)+diff*diff
835 sums(7)=sums(7)+1.0_mpd
837 sums(9)=sums(9)+diff*diff
840 IF (mod(i,3) == 1)
THEN
844 sums(4)=sums(4)+1.0_mpd
846 sums(6)=sums(6)+diff*diff
853 sums(7)=sums(7)+1.0_mpd
855 sums(9)=sums(9)+diff*diff
860 sums(2)=sums(2)/sums(1)
861 sums(3)=sqrt(sums(3)/sums(1))
862 sums(5)=sums(5)/sums(4)
863 sums(6)=sqrt(sums(6)/sums(4))
864 WRITE(*,143)
'Parameters 1 - 500: mean =',sums(2),
'rms =',sums(3)
865 WRITE(*,143)
'Parameters 501 - 700: mean =',sums(5),
'rms =',sums(6)
866 IF (sums(7) > 0.5_mpd)
THEN
867 sums(8)=sums(8)/sums(7)
868 sums(9)=sqrt(sums(9)/sums(7))
869 WRITE(*,143)
'Parameter pulls, all: mean =',sums(8),
'rms =',sums(9)
886 IF (mod(i,3) == 1)
THEN
906 WRITE(8,*)
'Record',
nrec1,
' has largest residual:',
value1
909 WRITE(8,*)
'Record',
nrec2,
' has largest Chi^2/Ndf:',
value2
913 WRITE(8,*)
'Record',
nrec3,
' is first with error (rank deficit/NaN)'
917 WRITE(8,*)
'In total 3 +',
nloopn,
' loops through the data files'
919 WRITE(8,*)
'In total 2 +',
nloopn,
' loops through the data files'
923 WRITE(8,*)
'In total ',
mnrsit,
' internal MINRES iterations'
931 ntsec=nint(deltat,mpi)
932 CALL sechms(deltat,nhour,minut,secnd)
933 nsecnd=nint(secnd,mpi)
934 WRITE(8,*)
'Total time =',ntsec,
' seconds =',nhour,
' h',minut, &
935 ' m',nsecnd,
' seconds'
937 WRITE(8,*)
'end ', chdate
946 WRITE(8,*)
'Data rejected in last iteration: '
948 nrejec(0),
' (rank deficit/NaN) ',
nrejec(1),
' (Ndf=0) ', &
955 WRITE(*,*)
' < Millepede II-P ending ... ', chdate
961102
FORMAT(2x,i4,2x,3f10.5,2x,3f10.5)
962103
FORMAT(
' Times [in sec] for text processing',f12.3/ &
965 ' func. value ',f12.3,
' *',f4.0/ &
966 ' func. value, global matrix, solution',f12.3,
' *',f4.0/ &
967 ' new solution',f12.3,
' *',f4.0/)
968105
FORMAT(
' Peak dynamic memory allocation: ',f11.6,
' GB')
988 INTEGER(mpi) :: istop
989 INTEGER(mpi) :: itgbi
990 INTEGER(mpi) :: itgbl
992 INTEGER(mpi) :: itnlim
995 INTEGER(mpi),
INTENT(IN) :: ivgbi
1006 INTEGER(mpl) :: ijadd
1034 globalcorrections, itnlim, nout, rtol, istop, itn, anorm, acond, rnorm, arnorm, ynorm)
1038 globalcorrections, itnlim, nout, rtol, istop, itn, anorm, acond, rnorm, arnorm, ynorm)
1041 globalcorrections, itnlim, nout, rtol, istop, itn, anorm, acond, rnorm, arnorm, ynorm)
1047 err=sqrt(abs(real(gmati,mps)))
1048 IF(gmati < 0.0_mpd) err=-err
1052 ELSE IF(
matsto == 2)
THEN
1053 jk=ijadd(ivgbi,ivgbi)
1060 gcor2=real(1.0_mpd-1.0_mpd/(gmati*diag),mps)
1062101
FORMAT(1x,
' label parameter presigma differ', &
1063 ' Error gcor^2 iit'/ 1x,
'---------',2x,5(
'-----------'),2x,
'----')
1064102
FORMAT(i10,2x,4g12.4,f7.4,i6,i4)
1084 INTEGER(mpi) :: istop
1085 INTEGER(mpi) :: itgbi
1086 INTEGER(mpi) :: itgbl
1088 INTEGER(mpi) :: itnlim
1089 INTEGER(mpi) :: nout
1091 INTEGER(mpi),
INTENT(IN) :: ivgbi
1099 INTEGER(mpl) :: ijadd
1122 mxxnrm = real(
nagb,mpd)/sqrt(epsilon(mxxnrm))
1124 trcond = 1.0_mpd/epsilon(trcond)
1125 ELSE IF(
mrmode == 2)
THEN
1133 itnlim=itnlim, rtol=rtol, maxxnorm=mxxnrm, trancond=trcond, &
1137 itnlim=itnlim, rtol=rtol, maxxnorm=mxxnrm, trancond=trcond, &
1141 itnlim=itnlim, rtol=rtol, maxxnorm=mxxnrm, trancond=trcond, &
1148 err=sqrt(abs(real(gmati,mps)))
1149 IF(gmati < 0.0_mpd) err=-err
1153 ELSE IF(
matsto == 2)
THEN
1154 jk=ijadd(ivgbi,ivgbi)
1161 gcor2=real(1.0_mpd-1.0_mpd/(gmati*diag),mps)
1163101
FORMAT(1x,
' label parameter presigma differ', &
1164 ' Error gcor^2 iit'/ 1x,
'---------',2x,5(
'-----------'),2x,
'----')
1165102
FORMAT(i10,2x,4g12.4,f7.4,i6,i4)
1178 INTEGER(mpi) :: icgb
1179 INTEGER(mpi) :: irhs
1180 INTEGER(mpi) :: itgbi
1181 INTEGER(mpi) :: ivgb
1183 INTEGER(mpi) :: jcgb
1185 INTEGER(mpi) :: label
1187 INTEGER(mpi) :: inone
1190 REAL(mpd) :: drhs(4)
1191 INTEGER(mpi) :: idrh (4)
1216 IF(abs(rhs) > climit)
THEN
1222 WRITE(*,101) (idrh(l),drhs(l),l=1,irhs)
1231 WRITE(*,101) (idrh(l),drhs(l),l=1,irhs)
1234 WRITE(*,102)
' Constraints: only equation values >', climit,
' are printed'
1235101
FORMAT(
' ',4(i4,g11.3))
1236102
FORMAT(a,g11.2,a)
1249 INTEGER(mpi) :: icgb
1250 INTEGER(mpi) :: isblck
1251 INTEGER(mpi) :: ilast
1252 INTEGER(mpi) :: itgbi
1253 INTEGER(mpi) :: ivgb
1254 INTEGER(mpi) :: jcgb
1255 INTEGER(mpi) :: label
1256 INTEGER(mpi) :: labelf
1257 INTEGER(mpi) :: labell
1258 INTEGER(mpi) :: ncon
1259 INTEGER(mpi) :: npar
1260 INTEGER(mpi) :: nconmx
1261 INTEGER(mpi) :: nparmx
1262 INTEGER(mpi) :: inone
1263 INTEGER(mpi) :: itype
1264 INTEGER(mpi) :: ncgbw
1265 INTEGER(mpi) :: newlen
1266 INTEGER(mpi) :: nvar
1267 INTEGER(mpi) :: last
1268 INTEGER(mpi) :: lastlen
1270 INTEGER(mpl):: length
1271 INTEGER(mpl) :: rows
1288 IF(last == 0.AND.label < 0)
THEN
1289 IF (
ncgb > 0 .AND.
icheck>0)
WRITE(*,113)
ncgb, newlen-lastlen-3, nvar
1291 IF (nvar == 0 .AND.
iskpec > 0)
THEN
1302 IF(itype == 2) ncgbw=ncgbw+1
1309 IF (ivgb > 0) nvar=nvar+1
1311 IF(label > 0.AND.itype == 2)
THEN
1319 IF (
ncgb > 0 .AND.
icheck>0)
WRITE(*,113)
ncgb, newlen-lastlen-2, nvar
1321 IF (nvar == 0 .AND.
iskpec > 0) newlen=lastlen
1325 WRITE(*,*)
'PRPCON:',
ncgbe,
' empty constraints skipped'
1328 IF (ncgbw == 0)
THEN
1329 WRITE(*,*)
'PRPCON:',
ncgb,
' constraints accepted'
1331 WRITE(*,*)
'PRPCON:',
ncgb,
' constraints accepted,',ncgbw,
'weighted'
1338 length=
ncgb+1; rows=3
1389 WRITE(*,*)
' Cons. sorted', jcgb, icgb,
vecconsstart(icgb), labelf, labell
1398 nconmx=max(nconmx,ncon)
1399 nparmx=max(nparmx,npar)
1405 WRITE(*,*)
' Cons. block ',
ncblck, isblck, jcgb, labelf, labell
1415 WRITE(*,*)
'PRPCON: constraints split into ',
ncblck,
'(disjoint) blocks'
1416 WRITE(*,*)
' max block size (cons., par.) ', nconmx, nparmx
1419113
FORMAT(
' constraint',i6,
' : ',i9,
' parameters,',i9,
' variable')
1435 INTEGER(mpi) :: iblck
1436 INTEGER(mpi) :: icgb
1438 INTEGER(mpi) :: ifirst
1439 INTEGER(mpi) :: ilast
1440 INTEGER(mpi) :: ioffc
1441 INTEGER(mpi) :: ioffp
1442 INTEGER(mpi) :: irank
1443 INTEGER(mpi) :: ipar0
1444 INTEGER(mpi) :: itgbi
1445 INTEGER(mpi) :: ivgb
1447 INTEGER(mpi) :: jcgb
1449 INTEGER(mpi) :: label
1450 INTEGER(mpi) :: ncon
1451 INTEGER(mpi) :: npar
1452 INTEGER(mpi) :: nrank
1453 INTEGER(mpi) :: inone
1458 INTEGER(mpl):: length
1459 REAL(mpd),
DIMENSION(:),
ALLOCATABLE :: matConstraintsT
1460 REAL(mpd),
DIMENSION(:),
ALLOCATABLE :: auxVectorD
1461 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: auxVectorI
1465 IF(
ncgb == 0)
RETURN
1476 CALL mpalloc(auxvectori,length,
'auxiliary array (I)')
1477 CALL mpalloc(auxvectord,length,
'auxiliary array (D)')
1480 CALL mpalloc(matconstraintst,length,
'transposed matrix of constraints (blocks)')
1481 matconstraintst=0.0_mpd
1493 DO jcgb=ifirst,ilast
1505 IF(ivgb > 0) matconstraintst(ivgb-ipar0+ioffc+(jcgb-ifirst)*npar)=factr
1513 DO l=ioffc+1,ioffc+npar
1524 IF (
icheck > 1)
WRITE(*,*)
' Constraint block rank', iblck, ncon, irank
1526 ioffc=ioffc+npar*ncon
1533 WRITE(*,*)
'Rank of product matrix of constraints is',nrank, &
1534 ' for',
ncgb,
' constraint equations'
1535 WRITE(8,*)
'Rank of product matrix of constraints is',nrank, &
1536 ' for',
ncgb,
' constraint equations'
1537 IF(nrank <
ncgb)
THEN
1538 WRITE(*,*)
'Warning: insufficient constraint equations!'
1539 WRITE(8,*)
'Warning: insufficient constraint equations!'
1542 WRITE(*,*)
' --> enforcing SUBITO mode'
1543 WRITE(8,*)
' --> enforcing SUBITO mode'
1550 print *,
'QL decomposition of constraints matrix'
1554 print *,
' largest |eigenvalue| of L: ', evmax
1555 print *,
' smallest |eigenvalue| of L: ', evmin
1556 IF (evmin == 0.0_mpd)
THEN
1557 CALL peend(27,
'Aborted, singular QL decomposition of constraints matrix')
1558 stop
'FEASMA: stopping due to singular QL decomposition of constraints matrix'
1584 INTEGER(mpi) :: icgb
1585 INTEGER(mpi) :: iter
1586 INTEGER(mpi) :: itgbi
1587 INTEGER(mpi) :: ivgb
1588 INTEGER(mpi) :: iblck
1589 INTEGER(mpi) :: ieblck
1590 INTEGER(mpi) :: isblck
1591 INTEGER(mpi) :: ifirst
1592 INTEGER(mpi) :: ilast
1594 INTEGER(mpi) :: jcgb
1595 INTEGER(mpi) :: label
1596 INTEGER(mpi) :: inone
1597 INTEGER(mpi) :: ncon
1599 REAL(mps),
INTENT(IN) :: concut
1600 INTEGER(mpi),
INTENT(OUT) :: iact
1607 REAL(mpd),
DIMENSION(:),
ALLOCATABLE :: vecCorrections
1641 sum1=sqrt(sum1/real(
ncgb,mpd))
1642 sum2=sum2/real(
ncgb,mpd)
1644 IF(iter == 1.AND.sum1 < concut)
RETURN
1646 IF(iter == 1.AND.
ncgb <= 12)
THEN
1648 WRITE(*,*)
'Constraint equation discrepancies:'
1650101
FORMAT(4x,4(i5,g12.4))
1652103
FORMAT(10x,
' Cut on rms value is',g8.1)
1657 WRITE(*,*)
'Improve constraints'
1661 WRITE(*,102) iter,sum1,sum2,sum3
1662102
FORMAT(i6,
' rms',g12.4,
' avrg_abs',g12.4,
' max_abs',g12.4)
1664 CALL mpalloc(veccorrections,int(
nvgb,mpl),
'constraint corrections')
1665 veccorrections=0.0_mpd
1673 ieblck=isblck+(ncon*(ncon+1))/2
1742 INTEGER(mpi) :: iact
1743 INTEGER(mpi) :: ierrc
1744 INTEGER(mpi) :: ierrf
1745 INTEGER(mpi) :: inder
1746 INTEGER(mpi) :: ioffp
1748 INTEGER(mpi) :: ithr
1749 INTEGER(mpi) :: jfile
1750 INTEGER(mpi) :: jrec
1752 INTEGER(mpi) :: kfile
1755 INTEGER(mpi) :: mpri
1757 INTEGER(mpi) :: nact
1758 INTEGER(mpi) :: nbuf
1759 INTEGER(mpi) :: ndata
1760 INTEGER(mpi) :: noff
1761 INTEGER(mpi) :: noffs
1762 INTEGER(mpi) :: npointer
1763 INTEGER(mpi) :: npri
1767 INTEGER(mpi) :: nrpr
1768 INTEGER(mpi) :: nthr
1769 INTEGER(mpi) :: ntot
1770 INTEGER(mpi) :: maxRecordSize
1771 INTEGER(mpi) :: maxRecordFile
1773 INTEGER(mpi),
INTENT(OUT) :: more
1783 CHARACTER (LEN=7) :: cfile
1790 DATA npri / 0 /, mpri / 1000 /
1839 files:
DO WHILE (jfile > 0)
1843 CALL binopn(kfile,ithr,ios)
1849 IF(kfile <=
nfilf)
THEN
1870 eof=(ierrc <= 0.AND.ierrc /= -4)
1871 IF(eof.AND.ierrc < 0)
THEN
1872 WRITE(*,*)
'Read error for binary Cfile', kfile,
'record', jrec+1,
':', ierrc
1873 WRITE(8,*)
'Read error for binary Cfile', kfile,
'record', jrec+1,
':', ierrc
1875 WRITE(cfile,
'(I7)') kfile
1876 CALL peend(18,
'Aborted, read error(s) for binary file ' // cfile)
1877 stop
'PEREAD: stopping due to read errors'
1881 IF(eof)
EXIT records
1886 xfd(jfile)=max(
xfd(jfile),n)
1889 IF(inder(noff+1) /= 0)
CALL hmpent(8,real(inder(noff+1),mps))
1903 IF ((noff+nr+2+
ndimbuf >= ndata*(iact+1)).OR.(nbuf >= npointer))
EXIT files
1918 IF (
kfd(1,jfile) == 1)
THEN
1919 print *,
'PEREAD: file ', kfile,
'read the first time, found',jrec,
' records'
1970 IF (
kfd(1,jfile) == 1)
kfd(1,jfile)=-nrc
1989 DO WHILE (
nloopn == 0.AND.ntot >= nrpr)
1990 WRITE(*,*)
' Record ',nrpr
1991 IF (nrpr < 100000)
THEN
2000 IF (npri == 1)
WRITE(*,100)
2002100
FORMAT(/
' PeRead records active file' &
2003 /
' total block threads number')
2004101
FORMAT(
' PeRead',4i10)
2017 IF (
xfd(k) > maxrecordsize)
THEN
2018 maxrecordsize=
xfd(k)
2021 dw=real(-
kfd(1,k),mpd)
2024 ds1=ds1+dw*real(
wfd(k),mpd)
2025 ds2=ds2+dw*real(
wfd(k)**2,mpd)
2027 print *,
'PEREAD: file ', maxrecordfile,
'with max record size ', maxrecordsize
2028 IF (
nfilw > 0.AND.ds0 > 0.0_mpd)
THEN
2032 WRITE(lun,177)
nfilw,real(ds1,mps),real(ds2,mps)
2033177
FORMAT(/
' !!!!!',i4,
' weighted binary files', &
2034 /
' !!!!! mean, variance of weights =',2g12.4)
2052179
FORMAT(/
' Read cache usage (#blocks, #records, ', &
2053 'min,max records/block'/17x,4i10)
2071 INTEGER(mpi),
INTENT(IN) :: mode
2073 INTEGER(mpi) :: ibuf
2074 INTEGER(mpi) :: ichunk
2075 INTEGER(mpi) :: iproc
2076 INTEGER(mpi) :: isfrst
2077 INTEGER(mpi) :: islast
2084 INTEGER(mpi),
PARAMETER :: maxbad = 100
2085 INTEGER(mpi) :: nbad
2086 INTEGER(mpi) :: nerr
2087 INTEGER(mpi) :: inone
2107 CALL isjajb(nst,ist,ja,jb,jsp)
2129 CALL pechk(ibuf,nerr)
2132 IF(nbad >= maxbad)
EXIT
2137 CALL isjajb(nst,ist,ja,jb,jsp)
2146 CALL peend(20,
'Aborted, bad binary records')
2147 stop
'PEREAD: stopping due to bad records'
2169 INTEGER(mpi) :: inder
2170 INTEGER(mpi) :: ioff
2171 INTEGER(mpi) :: isfrst
2172 INTEGER(mpi) :: islast
2179 INTEGER(mpi),
INTENT(IN) :: ibuf
2180 INTEGER(mpi),
INTENT(OUT) :: nerr
2193 outer:
DO WHILE(is < nst)
2198 IF(is > nst)
EXIT outer
2199 IF(inder(is) == 0)
EXIT inner1
2204 IF(is > nst)
EXIT outer
2205 IF(inder(is) == 0)
EXIT inner2
2208 IF(ja+1 == jb.AND.glder(jb) < 0.0_mpr8)
THEN
2211 is=is+nint(-glder(jb),mpi)
2214 DO WHILE(inder(is+1) /= 0.AND.is < nst)
2221100
FORMAT(
' PEREAD: record ', i8,
' in file ',i6,
' is broken !!!')
2231101
FORMAT(
' PEREAD: record ', i8,
' in file ',i6,
' contains ', i6,
' NaNs !!!')
2265 INTEGER(mpi) :: inder
2267 INTEGER(mpi),
INTENT(IN) :: nst
2268 INTEGER(mpi),
INTENT(IN OUT) :: is
2269 INTEGER(mpi),
INTENT(OUT) :: ja
2270 INTEGER(mpi),
INTENT(OUT) :: jb
2271 INTEGER(mpi),
INTENT(OUT) :: jsp
2281 IF(is >= nst)
RETURN
2284 IF(inder(is) == 0)
EXIT
2289 IF(inder(is) == 0)
EXIT
2292 IF(ja+1 == jb.AND.glder(jb) < 0.0_mpr8)
THEN
2295 is=is+nint(-glder(jb),mpi)
2298 DO WHILE(inder(is+1) /= 0.AND.is < nst)
2336 INTEGER(mpi) :: ibuf
2337 INTEGER(mpi) :: inder
2338 INTEGER(mpi) :: ioffb
2340 INTEGER(mpi) :: isfrst
2341 INTEGER(mpi) :: islast
2342 INTEGER(mpi) :: itgbi
2343 INTEGER(mpi) :: itgbij
2344 INTEGER(mpi) :: itgbik
2345 INTEGER(mpi) :: ivgb
2346 INTEGER(mpi) :: ivgbij
2347 INTEGER(mpi) :: ivgbik
2350 INTEGER(mpi) :: lastit
2352 INTEGER(mpi) :: ncrit
2353 INTEGER(mpi) :: ndfs
2354 INTEGER(mpi) :: ngras
2355 INTEGER(mpi) :: nparl
2357 INTEGER(mpi) :: nrej
2358 INTEGER(mpi) :: inone
2359 INTEGER(mpi) :: ilow
2360 INTEGER(mpi) :: nlow
2361 INTEGER(mpi) :: nzero
2387 CALL gmpdef(1,4,
'Function value in iterations')
2389 CALL gmpdef(2,3,
'Number of MINRES steps vs iteration nr')
2391 CALL hmpdef( 5,0.0,0.0,
'Number of degrees of freedom')
2392 CALL hmpdef(11,0.0,0.0,
'Number of local parameters')
2393 CALL hmpdef(23,0.0,0.0,
'SQRT of diagonal elements without presigma')
2394 CALL hmpdef(24,0.0,0.0,
'Log10 of off-diagonal elements')
2395 CALL hmpdef(25,0.0,0.0,
'Relative individual pre-sigma')
2396 CALL hmpdef(26,0.0,0.0,
'Relative global pre-sigma')
2401 'Normalized residuals of single (global) measurement')
2403 'Normalized residuals of single (local) measurement')
2405 'Pulls of single (global) measurement')
2407 'Pulls of single (local) measurement')
2408 CALL hmpdef(4,0.0,0.0,
'Chi^2/Ndf after local fit')
2409 CALL gmpdef(4,5,
'location, dispersion (res.) vs record nr')
2410 CALL gmpdef(5,5,
'location, dispersion (pull) vs record nr')
2424 IF(
nloopn == 2)
CALL hmpdef(6,0.0,0.0,
'Down-weight fraction')
2427 IF(
iterat /= lastit)
THEN
2435 ELSE IF(
iterat >= 1)
THEN
2491 WRITE(*,111) nparl,ncrit,usei,used,peaki,peakd
2492111
FORMAT(
' Write cache usage (#flush,#overrun,<levels>,', &
2493 'peak(levels))'/2i7,
',',4(f6.1,
'%'))
2513 print *,
" ... warning ..."
2514 print *,
" global parameters with too few (< MREQENA) accepted entries: ", nlow
2518 IF(
icalcm == 1 .AND. nzero > 0)
THEN
2520 WRITE(*,*)
'Warning: the rank defect of the symmetric',
nfgb, &
2521 '-by-',
nfgb,
' matrix is ',
ndefec,
' (should be zero).'
2522 WRITE(lun,*)
'Warning: the rank defect of the symmetric',
nfgb, &
2523 '-by-',
nfgb,
' matrix is ',
ndefec,
' (should be zero).'
2526 WRITE(*,*)
' --> enforcing SUBITO mode'
2527 WRITE(lun,*)
' --> enforcing SUBITO mode'
2546 IF(elmt > 0.0)
CALL hmpent(23,1.0/sqrt(elmt))
2595 IF(sgm > 0.0) weight=1.0/sgm**2
2611 IF(itgbij /= 0)
THEN
2622 IF(
icalcm == 1.AND.ivgbij > 0)
THEN
2629 IF(ivgbij > 0.AND.ivgbik > 0)
THEN
2630 CALL mupdat(ivgbij,ivgbik,weight*factrj*factrk)
2636 adder=weight*dsum**2
2653 ratae=max(-99.9,ratae)
2662 WRITE(*,*)
'Data rejected in initial loop:'
2664 nrejec(0),
' (rank deficit/NaN) ',
nrejec(1),
' (Ndf=0) ', &
2701 IF(
nloopn == 2)
CALL hmpwrt(6)
2713 WRITE(lun,*)
' === local fits have bordered band matrix structure ==='
2714 IF (
nbndr(1) > 0 )
WRITE(lun,101)
' NBNDR',
nbndr(1),
'number of records (upper/left border)'
2715 IF (
nbndr(2) > 0 )
WRITE(lun,101)
' NBNDR',
nbndr(2),
'number of records (lower/right border)'
2716 WRITE(lun,101)
' NBDRX',
nbdrx,
'max border size'
2717 WRITE(lun,101)
' NBNDX',
nbndx,
'max band width'
2726101
FORMAT(1x,a8,
' =',i10,
' = ',a)
2741 INTEGER(mpi),
INTENT(IN) :: lunp
2746101
FORMAT(
' it fc',
' fcn_value dfcn_exp slpr costh iit st', &
2747 ' ls step cutf',1x,
'rejects hhmmss FMS')
2748102
FORMAT(
' -- --',
' ----------- -------- ---- ----- --- --', &
2749 ' -- ----- ----',1x,
'------- ------ ---')
2765 INTEGER(mpi) :: nrej
2766 INTEGER(mpi) :: nsecnd
2770 REAL(mps) :: slopes(3)
2771 REAL(mps) :: steps(3)
2772 REAL,
DIMENSION(2) :: ta
2774 INTEGER(mpi),
INTENT(IN) :: lunp
2776 CHARACTER (LEN=4):: ccalcm(4)
2777 DATA ccalcm /
' end',
' S',
' F ',
' FMS' /
2781 IF(nrej > 9999999) nrej=9999999
2785 nsecnd=nint(secnd,mpi)
2794 CALL ptlopt(nfa,ma,slopes,steps)
2795 ratae=max(-99.9,min(99.9,slopes(2)/slopes(1)))
2801103
FORMAT(i3,i3,e12.5,38x,f5.1, 1x,i7, i3,i2.2,i2.2,a4)
2802104
FORMAT(i3,i3,e12.5,1x,e8.2,f6.3,f6.3,i5,2i3,f6.3,f5.1, &
2803 1x,i7, i3,i2.2,i2.2,a4)
2804105
FORMAT(i3,i3,e12.5,1x,e8.2,12x,i5,i3,9x,f5.1, &
2805 1x,i7, i3,i2.2,i2.2,a4)
2818 INTEGER(mpi) :: minut
2820 INTEGER(mpi) :: nhour
2821 INTEGER(mpi) :: nrej
2822 INTEGER(mpi) :: nsecnd
2826 REAL(mps) :: slopes(3)
2827 REAL(mps) :: steps(3)
2828 REAL,
DIMENSION(2) :: ta
2830 INTEGER(mpi),
INTENT(IN) :: lunp
2831 CHARACTER (LEN=4):: ccalcm(4)
2832 DATA ccalcm /
' end',
' S',
' F ',
' FMS' /
2836 IF(nrej > 9999999) nrej=9999999
2840 nsecnd=nint(secnd,mpi)
2844 CALL ptlopt(nfa,ma,slopes,steps)
2845 ratae=abs(slopes(2)/slopes(1))
2850104
FORMAT(3x,i3,e12.5,9x, 35x, i7, i3,i2.2,i2.2,a4)
2851105
FORMAT(3x,i3,e12.5,9x, f6.3,14x,i3,f6.3,6x, i7, i3,i2.2,i2.2,a4)
2865 INTEGER(mpi) :: nsecnd
2868 REAL,
DIMENSION(2) :: ta
2871 INTEGER(mpi),
INTENT(IN) :: lunp
2872 CHARACTER (LEN=4):: ccalcm(4)
2873 DATA ccalcm /
' end',
' S',
' F ',
' FMS' /
2878 nsecnd=nint(secnd,mpi)
2880 WRITE(lunp,106) nhour,minut,nsecnd,ccalcm(
lcalcm)
2881106
FORMAT(69x,i3,i2.2,i2.2,a4)
2891 INTEGER(mpi) :: lunit
2893 WRITE(lunit,102)
'Explanation of iteration table'
2894 WRITE(lunit,102)
'=============================='
2895 WRITE(lunit,101)
'it', &
2896 'iteration number. Global parameters are improved for it > 0.'
2897 WRITE(lunit,102)
'First function evaluation is called iteraton 0.'
2898 WRITE(lunit,101)
'fc',
'number of function evaluations.'
2899 WRITE(lunit,101)
'fcn_value',
'value of 2 x Likelihood function (LF).'
2900 WRITE(lunit,102)
'The final value is the chi^2 value of the fit and should'
2901 WRITE(lunit,102)
'be about equal to the NDF (see below).'
2902 WRITE(lunit,101)
'dfcn_exp', &
2903 'expected reduction of the value of the Likelihood function (LF)'
2904 WRITE(lunit,101)
'slpr',
'ratio of the actual slope to inital slope.'
2905 WRITE(lunit,101)
'costh', &
2906 'cosine of angle between search direction and -gradient'
2908 WRITE(lunit,101)
'iit', &
2909 'number of internal iterations in MINRES algorithm'
2910 WRITE(lunit,101)
'st',
'stop code of MINRES algorithm'
2911 WRITE(lunit,102)
'< 0: rhs is very special, with beta2 = 0'
2912 WRITE(lunit,102)
'= 0: rhs b = 0, i.e. the exact solution is x = 0'
2913 WRITE(lunit,102)
'= 1 requested accuracy achieved, as determined by rtol'
2914 WRITE(lunit,102)
'= 2 reasonable accuracy achieved, given eps'
2915 WRITE(lunit,102)
'= 3 x has converged to an eigenvector'
2916 WRITE(lunit,102)
'= 4 matrix ill-conditioned (Acond has exceeded 0.1/eps)'
2917 WRITE(lunit,102)
'= 5 the iteration limit was reached'
2918 WRITE(lunit,102)
'= 6 Matrix x vector does not define a symmetric matrix'
2919 WRITE(lunit,102)
'= 7 Preconditioner does not define a symmetric matrix'
2920 ELSEIF (
metsol == 4)
THEN
2921 WRITE(lunit,101)
'iit', &
2922 'number of internal iterations in MINRES-QLP algorithm'
2923 WRITE(lunit,101)
'st',
'stop code of MINRES-QLP algorithm'
2924 WRITE(lunit,102)
'= 1: beta_{k+1} < eps, iteration k is the final Lanczos step.'
2925 WRITE(lunit,102)
'= 2: beta2 = 0. If M = I, b and x are eigenvectors of A.'
2926 WRITE(lunit,102)
'= 3: beta1 = 0. The exact solution is x = 0.'
2927 WRITE(lunit,102)
'= 4: A solution to (poss. singular) Ax = b found, given rtol.'
2928 WRITE(lunit,102)
'= 5: A solution to (poss. singular) Ax = b found, given eps.'
2929 WRITE(lunit,102)
'= 6: Pseudoinverse solution for singular LS problem, given rtol.'
2930 WRITE(lunit,102)
'= 7: Pseudoinverse solution for singular LS problem, given eps.'
2931 WRITE(lunit,102)
'= 8: The iteration limit was reached.'
2932 WRITE(lunit,102)
'= 9: The operator defined by Aprod appears to be unsymmetric.'
2933 WRITE(lunit,102)
'=10: The operator defined by Msolve appears to be unsymmetric.'
2934 WRITE(lunit,102)
'=11: The operator defined by Msolve appears to be indefinite.'
2935 WRITE(lunit,102)
'=12: xnorm has exceeded maxxnorm or will exceed it next iteration.'
2936 WRITE(lunit,102)
'=13: Acond has exceeded Acondlim or 0.1/eps.'
2937 WRITE(lunit,102)
'=14: Least-squares problem but no converged solution yet.'
2938 WRITE(lunit,102)
'=15: A null vector obtained, given rtol.'
2940 WRITE(lunit,101)
'ls',
'line search info'
2941 WRITE(lunit,102)
'< 0 recalculate function'
2942 WRITE(lunit,102)
'= 0: N or STP lt 0 or step not descending'
2943 WRITE(lunit,102)
'= 1: Linesearch convergence conditions reached'
2944 WRITE(lunit,102)
'= 2: interval of uncertainty at lower limit'
2945 WRITE(lunit,102)
'= 3: max nr of line search calls reached'
2946 WRITE(lunit,102)
'= 4: step at the lower bound'
2947 WRITE(lunit,102)
'= 5: step at the upper bound'
2948 WRITE(lunit,102)
'= 6: rounding error limitation'
2949 WRITE(lunit,101)
'step', &
2950 'the factor for the Newton step during the line search. Usually'
2952 'a value of 1 gives a sufficient reduction of the LF. Oherwise'
2953 WRITE(lunit,102)
'other step values are tried.'
2954 WRITE(lunit,101)
'cutf', &
2955 'cut factor. Local fits are rejected, if their chi^2 value'
2957 'is larger than the 3-sigma chi^2 value times the cut factor.'
2958 WRITE(lunit,102)
'A cut factor of 1 is used finally, but initially a larger'
2959 WRITE(lunit,102)
'factor may be used. A value of 0.0 means no cut.'
2960 WRITE(lunit,101)
'rejects',
'total number of rejected local fits.'
2961 WRITE(lunit,101)
'hmmsec',
'the time in hours (h), minutes (mm) and seconds.'
2962 WRITE(lunit,101)
'FMS',
'calculation of Function value, Matrix, Solution.'
2965101
FORMAT(a9,
' = ',a)
2982 INTEGER(mpi),
INTENT(IN) :: i
2983 INTEGER(mpi),
INTENT(IN) :: j
2984 REAL(mpd),
INTENT(IN) :: add
2986 INTEGER(mpl):: ijadd
2991 IF(i <= 0.OR.j <= 0)
RETURN
2998 ELSE IF(
matsto == 2)
THEN
3014 CALL peend(23,
'Aborted, bad matrix index')
3015 stop
'mupdat: bad index'
3021 ELSE IF(
mbandw == 0)
THEN
3063SUBROUTINE loopbf(nrej,ndfs,sndf,dchi2s, numfil,naccf,chi2f,ndff)
3091 INTEGER(mpi) :: ibuf
3092 INTEGER(mpi) :: ichunk
3093 INTEGER(mpl) :: icmn
3094 INTEGER(mpl) :: icost
3096 INTEGER(mpi) :: idiag
3097 INTEGER(mpi) :: iext
3101 INTEGER(mpi) :: ijsym
3105 INTEGER(mpi) :: imeas
3107 INTEGER(mpi) :: inder
3109 INTEGER(mpi) :: ioffb
3110 INTEGER(mpi) :: ioffc
3111 INTEGER(mpi) :: ioffd
3112 INTEGER(mpi) :: ioffe
3113 INTEGER(mpi) :: ioffi
3114 INTEGER(mpi) :: iprdbg
3115 INTEGER(mpi) :: iproc
3116 INTEGER(mpi) :: irbin
3117 INTEGER(mpi) :: irow
3118 INTEGER(mpi) :: isfrst
3119 INTEGER(mpi) :: islast
3121 INTEGER(mpi) :: iter
3122 INTEGER(mpi) :: itgbi
3123 INTEGER(mpi) :: ivgbj
3124 INTEGER(mpi) :: ivgbk
3130 INTEGER(mpi) :: joffd
3131 INTEGER(mpi) :: joffi
3132 INTEGER(mpi) :: jproc
3135 INTEGER(mpi) :: kbdr
3136 INTEGER(mpi) :: kbdrx
3137 INTEGER(mpi) :: kbnd
3140 INTEGER(mpi) :: mbdr
3141 INTEGER(mpi) :: mbnd
3142 INTEGER(mpi) :: mside
3143 INTEGER(mpi) :: nalc
3144 INTEGER(mpi) :: nalg
3148 INTEGER(mpi) :: ndfsd
3149 INTEGER(mpi) :: ndown
3151 INTEGER(mpi) :: nfred
3152 INTEGER(mpi) :: nfrei
3154 INTEGER(mpi) :: nprdbg
3155 INTEGER(mpi) :: nrank
3158 INTEGER(mpi) :: nter
3159 INTEGER(mpi) :: nweig
3161 INTEGER(mpi),
INTENT(IN OUT) :: nrej(0:3)
3162 INTEGER(mpi),
INTENT(IN OUT) :: ndfs
3163 REAL(mpd),
INTENT(IN OUT) :: sndf
3164 REAL(mpd),
INTENT(IN OUT) :: dchi2s
3165 INTEGER(mpi),
INTENT(IN) :: numfil
3166 INTEGER(mpi),
INTENT(IN OUT) :: naccf(numfil)
3167 REAL(mps),
INTENT(IN OUT) :: chi2f(numfil)
3168 INTEGER(mpi),
INTENT(IN OUT) :: ndff(numfil)
3180 CHARACTER (LEN=3):: chast
3181 DATA chuber/1.345_mpd/
3182 DATA cauchy/2.3849_mpd/
3185 ijsym(i,j)=min(i,j)+(max(i,j)*max(i,j)-max(i,j))/2
3232 IF(
nloopn == 1.AND.mod(nrc,100000) == 0)
THEN
3233 WRITE(*,*)
'Record',nrc,
' ... still reading'
3243 IF(nrc ==
nrecpr) lprnt=.true.
3244 IF(nrc ==
nrecp2) lprnt=.true.
3245 IF(nrc ==
nrecer) lprnt=.true.
3250 IF (nprdbg == 1) iprdbg=iproc
3251 IF (iproc /= iprdbg) lprnt=.false.
3256 WRITE(1,*)
'------------------ Loop',
nloopn, &
3257 ': Printout for record',nrc,iproc
3268 CALL isjajb(nst,ist,ja,jb,jsp)
3270 IF(imeas == 0)
WRITE(1,1121)
3272 WRITE(1,1122) imeas,glder(ja),glder(jb), &
3273 (inder(ja+j),glder(ja+j),j=1,jb-ja-1)
32751121
FORMAT(/
'Measured value and local derivatives'/ &
3276 ' i measured std_dev index...derivative ...')
32771122
FORMAT(i3,2g12.4,3(i3,g12.4)/(27x,3(i3,g12.4)))
3283 CALL isjajb(nst,ist,ja,jb,jsp)
3285 IF(imeas == 0)
WRITE(1,1123)
32971123
FORMAT(/
'Global derivatives'/ &
3298 ' i label gindex vindex derivative ...')
32991124
FORMAT(i3,2(i9,i7,i7,g12.4)/(3x,2(i9,i7,i7,g12.4)))
33001125
FORMAT(i3,2(i9,i7,i7,g12.4))
3310 WRITE(1,*)
'Data corrections using values of global parameters'
3311 WRITE(1,*)
'=================================================='
3320 CALL isjajb(nst,ist,ja,jb,jsp)
3322 rmeas=real(glder(ja),mpd)
3341 IF (jb < ist)
WRITE(1,102) neq,glder(ja),rmeas,glder(jb)
3349101
FORMAT(
' index measvalue corrvalue sigma')
3350102
FORMAT(i6,2x,2g12.4,
' +-',g12.4)
3352 IF(nalc <= 0)
GO TO 90
3354 ngg=(nalg*nalg+nalg)/2
3390 DO WHILE(iter < nter)
3395 WRITE(1,*)
'Outlier-suppression iteration',iter,
' of',nter
3396 WRITE(1,*)
'=========================================='
3406 DO i=1,(nalc*nalc+nalc)/2
3415 CALL isjajb(nst,ist,ja,jb,jsp)
3417 rmeas=real(glder(ja),mpd)
3418 rerr =real(glder(jb),mpd)
3419 wght =1.0_mpd/rerr**2
3425 IF(abs(resid) > chuber*rerr)
THEN
3426 wght=wght*chuber*rerr/abs(resid)
3430 wght=wght/(1.0+(resid/rerr/cauchy)**2)
3434 IF(lprnt.AND.iter /= 1.AND.nter /= 1)
THEN
3436 IF(abs(resid) > chuber*rerr) chast=
'* '
3437 IF(abs(resid) > 3.0*rerr) chast=
'** '
3438 IF(abs(resid) > 6.0*rerr) chast=
'***'
3439 IF(imeas == 0)
WRITE(1,*)
'Second loop: accumulate'
3440 IF(imeas == 0)
WRITE(1,103)
3445 WRITE(1,104) imeas,rmeas,resid,rerr,r1,chast,r2
3447103
FORMAT(
' index corrvalue residuum sigma', &
3449104
FORMAT(i6,2x,2g12.4,
' +-',g12.4,f7.2,1x,a3,f8.2)
3453 blvec(ij)=
blvec(ij)+wght*rmeas*real(glder(ja+j),mpd)
3458 +wght*real(glder(ja+j),mpd)*real(glder(ja+k),mpd)
3464 im=min(nalc+1-ij,nalc+1-ik)
3471 IF (iter == 1.AND.nalc > 5.AND.
lfitbb > 0)
THEN
3474 icmn=int(nalc,mpl)**3
3480 icost=6*int(nalc-kbdr,mpl)*int(kbnd+kbdr+1,mpl)**2+2*int(kbdr,mpl)**3
3481 IF (icost < icmn)
THEN
3493 kbdr=max(kbdr,
ibandh(k+nalc))
3494 icost=6*int(nalc-kbdr,mpl)*int(kbnd+kbdr+1,mpl)**2+2*int(kbdr,mpl)**3
3495 IF (icost < icmn)
THEN
3519 IF (
icalcm == 1.OR.lprnt) inv=2
3520 IF (mside == 1)
THEN
3535 IF ((.NOT.(
blvec(k) <= 0.0_mpd)).AND. (.NOT.(
blvec(k) > 0.0_mpd))) nan=nan+1
3540 WRITE(1,*)
'Parameter determination:',nalc,
' parameters,',
' rank=',nrank
3541 WRITE(1,*)
'-----------------------'
3542 IF(ndown /= 0)
WRITE(1,*)
' ',ndown,
' data down-weighted'
3556 CALL isjajb(nst,ist,ja,jb,jsp)
3558 rmeas=real(glder(ja),mpd)
3559 rerr =real(glder(jb),mpd)
3560 wght =1.0_mpd/rerr**2
3565 rmloc=rmloc+real(glder(ja+j),mpd)*
blvec(ij)
3578 dvar=dvar+
clmat(jk)*real(glder(ja+j),mpd)*real(glder(ja+k),mpd)
3582 IF (0.999999_mpd/wght > dvar)
THEN
3583 pull=rmeas/sqrt(1.0_mpd/wght-dvar)
3586 CALL hmpent(13,real(pull,mps))
3587 CALL gmpms(5,rec,real(pull,mps))
3589 CALL hmpent(14,real(pull,mps))
3608 IF(iter == 1.AND.jb < ist.AND.lhist) &
3609 CALL gmpms(4,rec,real(rmeas/rerr,mps))
3611 dchi2=wght*rmeas*rmeas
3616 IF(abs(resid) > chuber*rerr)
THEN
3617 wght=wght*chuber*rerr/abs(resid)
3618 dchi2=2.0*chuber*(abs(resid)/rerr-0.5*chuber)
3621 wght=wght/(1.0_mpd+(resid/rerr/cauchy)**2)
3622 dchi2=log(1.0_mpd+(resid/rerr/cauchy)**2)*cauchy**2
3632 IF(abs(resid) > chuber*rerr) chast=
'* '
3633 IF(abs(resid) > 3.0*rerr) chast=
'** '
3634 IF(abs(resid) > 6.0*rerr) chast=
'***'
3635 IF(imeas == 0)
WRITE(1,*)
'Third loop: single residuals'
3636 IF(imeas == 0)
WRITE(1,105)
3640 IF(resid < 0.0) r1=-r1
3641 IF(resid < 0.0) r2=-r2
3642 WRITE(1,106) imeas,glder(ja),rmeas,rerr,r1,chast,r2
3644105
FORMAT(
' index corrvalue residuum sigma', &
3646106
FORMAT(i6,2x,2g12.4,
' +-',g12.4,f7.2,1x,a3,f8.2)
3648 IF(iter == nter)
THEN
3650 resmax=max(resmax,abs(rmeas)/rerr)
3653 IF(iter == 1.AND.lhist)
THEN
3655 CALL hmpent( 3,real(rmeas/rerr,mps))
3657 CALL hmpent(12,real(rmeas/rerr,mps))
3663 resing=(real(nweig,mps)-real(suwt,mps))/real(nweig,mps)
3665 IF(iter == 1)
CALL hmpent( 5,real(ndf,mps))
3666 IF(iter == 1)
CALL hmpent(11,real(nalc,mps))
3667 IF(
nloopn == 2.AND.iter == nter)
CALL hmpent(6,resing)
3671 WRITE(1,*)
'Chi^2=',summ,
' at',ndf,
' degrees of freedom: ', &
3672 '3-sigma limit is',chindl(3,ndf)*real(ndf,mps)
3673 WRITE(1,*) suwt,
' is sum of factors, compared to',nweig, &
3674 ' Downweight fraction:',resing
3676 IF(nrank /= nalc.OR.nan > 0)
THEN
3680 WRITE(1,*)
' rank deficit/NaN ', nalc, nrank, nan
3681 WRITE(1,*)
' ---> rejected!'
3688 WRITE(1,*)
' ---> rejected!'
3693 chndf=real(summ/real(ndf,mpd),mps)
3695 IF(iter == 1.AND.lhist)
CALL hmpent(4,chndf)
3699 sndf=sndf+real(ndf,mpd)*dw1
3710 chichi=chindl(3,ndf)*real(ndf,mps)
3714 IF(summ >
chhuge*chichi)
THEN
3717 WRITE(1,*)
' ---> rejected!'
3725 IF(summ > chlimt)
THEN
3727 WRITE(1,*)
' ---> rejected!'
3731 dchi2s=dchi2s+dchi2*dw1
3741 dchi2s=dchi2s+dchi2*dw1
3745 WRITE(1,*)
' ---> rejected!'
3757 naccf(kfl)=naccf(kfl)+1
3758 ndff(kfl) =ndff(kfl) +ndf
3759 chi2f(kfl)=chi2f(kfl)+chndf
3768 CALL isjajb(nst,ist,ja,jb,jsp)
3771 rmeas=real(glder(ja),mpd)
3772 rerr =real(glder(jb),mpd)
3773 wght =1.0_mpd/rerr**2
3774 dchi2=wght*rmeas*rmeas
3778 IF(resid > chuber*rerr)
THEN
3779 wght=wght*chuber*rerr/resid
3780 dchi2=2.0*chuber*(resid/rerr-0.5*chuber)
3783 dchi2s=dchi2s+dchi2*dw1
3791 +dw1*wght*rmeas*real(glder(jb+j),mpd)
3803 -dw1*wght*real(glder(jb+j),mpd)*real(glder(jb+k),mpd)
3848 IF (nfred < 0.OR.nfrei < 0)
THEN
3883 WRITE(1,*)
'------------------ End of printout for record',nrc
3990 INTEGER(mpi) :: icount
3994 INTEGER(mpi) :: imin
3995 INTEGER(mpi) :: iprlim
3996 INTEGER(mpi) :: isub
3997 INTEGER(mpi) :: itgbi
3998 INTEGER(mpi) :: itgbl
3999 INTEGER(mpi) :: ivgbi
4001 INTEGER(mpi) :: label
4008 INTEGER(mpi) :: labele(3)
4010 REAL(mps):: compnt(3)
4015 CALL mvopen(lup,
'millepede.res')
4018 WRITE(*,*)
' Result of fit for global parameters'
4019 WRITE(*,*)
' ==================================='
4024 WRITE(lup,*)
'Parameter ! first 3 elements per line are', &
4025 ' significant (if used as input)'
4041 err=sqrt(abs(real(gmati,mps)))
4042 IF(gmati < 0.0_mpd) err=-err
4045 IF(gmati*diag > 0.0_mpd)
THEN
4046 gcor2=1.0_mpd-1.0_mpd/(gmati*diag)
4047 IF(gcor2 >= 0.0_mpd.AND.gcor2 <= 1.0_mpd) gcor=real(sqrt(gcor2),mps)
4052 IF(itgbi <= iprlim)
THEN
4066 ELSE IF(itgbi == iprlim+1)
THEN
4067 WRITE(* ,*)
'... (further printout suppressed, but see log file)'
4081 ELSE IF (
igcorr /= 0)
THEN
4099 CALL mvopen(lup,
'millepede.eve')
4109 DO isub=0,min(15,imin-1)
4131 WRITE(lup,103) (labele(ie),compnt(ie),ie=1,iev)
4135 IF(iev /= 0)
WRITE(lup,103) (labele(ie),compnt(ie),ie=1,iev)
4142101
FORMAT(1x,
' label parameter presigma differ', &
4143 ' error'/ 1x,
'-----------',4x,4(
'-------------'))
4144102
FORMAT(i10,2x,4g14.5,f8.3)
4145103
FORMAT(3(i11,f11.7,2x))
4146110
FORMAT(i10,2x,2g14.5,28x,i12)
4147111
FORMAT(i10,2x,3g14.5,14x,i12)
4148112
FORMAT(i10,2x,4g14.5,i12)
4168 INTEGER(mpi) :: icount
4169 INTEGER(mpi) :: itgbi
4170 INTEGER(mpi) :: itgbl
4171 INTEGER(mpi) :: ivgbi
4173 INTEGER(mpi) :: ncon
4180 CALL mvopen(lup,
'millepede.res')
4181 WRITE(lup,*)
'*** Results of checking input only, no solution performed ***'
4182 WRITE(lup,*)
'! fixed-1: by pre-sigma, -2: by entries cut, -3: by iterated entries cut'
4183 WRITE(lup,*)
'! Label Value Pre-sigma Entries Constraints Status '
4193 IF (ivgbi <= 0)
THEN
4194 WRITE(lup,110) itgbl,par,presig,icount,ncon,ivgbi
4196 WRITE(lup,111) itgbl,par,presig,icount,ncon
4202 WRITE(lup,*)
'! Appearance statistics '
4203 WRITE(lup,*)
'! Label First file and record Last file and record #files #paired-par'
4211110
FORMAT(
' ! ',i10,2x,2g14.5,2i12,
' fixed',i2)
4212111
FORMAT(
' ! ',i10,2x,2g14.5,2i12,
' variable')
4213112
FORMAT(
' !.',i10,6i11)
4231 INTEGER(mpi) :: iencdb
4232 INTEGER(mpi) :: iencdm
4233 INTEGER(mpi) :: iproc
4240 INTEGER(mpi),
INTENT(IN) :: n
4241 REAL(mpd),
INTENT(IN) :: x(n)
4242 REAL(mpd),
INTENT(OUT) :: b(n)
4249 INTEGER(mpl) :: indij
4250 INTEGER(mpl) :: indid
4252 INTEGER(mpi) :: ichunk
4282 CALL peend(24,
'Aborted, vector/matrix size mismatch')
4283 stop
'AVPRD0: mismatched vector and matrix'
4286 iencdm=ishft(1,iencdb)-1
4321 b(i)=b(i)+dot_product(
globalmatd(indij+lj:indij+lj+jn-1),x(j:j+jn-1))
4364 b(i)=b(i)+real(
globalmatf(indij+lj),mpd)*x(j)
4370 b(j)=b(j)+real(
globalmatf(indij+lj),mpd)*x(i)
4371 b(i)=b(i)+real(
globalmatf(indij+lj),mpd)*x(j)
4399 INTEGER(mpi),
INTENT(IN) :: n
4400 REAL(mpd),
INTENT(IN) :: x(n)
4401 REAL(mpd),
INTENT(OUT) :: b(n)
4406 CALL peend(24,
'Aborted, vector/matrix size mismatch')
4407 stop
'AVPROD: mismatched vector and matrix'
4436 INTEGER(mpi) :: iencdb
4437 INTEGER(mpi) :: iencdm
4438 INTEGER(mpi) :: isgn
4439 INTEGER(mpi) :: ispc
4440 INTEGER(mpi) :: item2
4441 INTEGER(mpi) :: jtem
4442 INTEGER(mpi) :: jtemc
4443 INTEGER(mpi) :: jtemn
4445 INTEGER(mpi),
INTENT(IN) :: itema
4446 INTEGER(mpi),
INTENT(IN) :: itemb
4448 INTEGER(mpl) ::
ijadd
4453 INTEGER(mpl) :: indid
4457 INTEGER(mpl) :: item1
4461 item1=max(itema,itemb)
4462 item2=min(itema,itemb)
4463 IF(item2 <= 0.OR.item1 > nd)
RETURN
4464 IF(item1 == item2)
THEN
4470 iencdm=ishft(1,iencdb)-1
4472 outer:
DO ispc=1,
nspc
4485 IF(ku < kl) cycle outer
4489 jtem =ishft(jtemc,-iencdb)
4490 jtemn=jtem+iand(jtemc,iencdm)
4491 IF(item2 >= jtem.AND.item2 < jtemn)
EXIT
4492 IF(item2 < jtem)
THEN
4494 ELSE IF(item2 >= jtemn)
THEN
4509 IF(ku < kl) cycle outer
4514 IF(item2 == jtem)
EXIT
4515 IF(item2 < jtem)
THEN
4517 ELSE IF(item2 > jtem)
THEN
4540 INTEGER(mpi) :: ichunk
4541 INTEGER(mpi) :: iencdb
4542 INTEGER(mpi) :: iencdm
4544 INTEGER(mpi) :: ispc
4546 INTEGER(mpi) :: jtem
4547 INTEGER(mpi) :: jtemc
4548 INTEGER(mpi) :: jtemn
4552 INTEGER(mpl) :: ijadd
4557 INTEGER(mpl) :: indid
4566 iencdm=ishft(1,iencdb)-1
4574 ir=i+(ispc-1)*(nd+1)
4590 jtem =ishft(jtemc,-iencdb)
4591 jtemn=jtem+iand(jtemc,iencdm)
4621 REAL(mps),
INTENT(IN) :: deltat
4622 INTEGER(mpi),
INTENT(OUT) :: minut
4623 INTEGER(mpi),
INTENT(OUT):: nhour
4624 REAL(mps),
INTENT(OUT):: secnd
4625 INTEGER(mpi) :: nsecd
4628 nsecd=nint(deltat,
mpi)
4630 minut=nsecd/60-60*nhour
4631 secnd=deltat-60*(minut+60*nhour)
4667 INTEGER(mpi),
INTENT(IN) :: item
4671 INTEGER(mpl) :: length
4672 INTEGER(mpl),
PARAMETER :: two = 2
4675 IF(item <= 0)
RETURN
4694 IF(j == 0)
EXIT inner
4713 IF (
lvllog > 1)
WRITE(
lunlog,*)
'INONE: array increased to', &
4733 INTEGER(mpi) :: iprime
4734 INTEGER(mpi) :: nused
4735 LOGICAL :: finalUpdate
4736 INTEGER(mpl) :: oldLength
4737 INTEGER(mpl) :: newLength
4738 INTEGER(mpl),
PARAMETER :: two = 2
4739 INTEGER(mpi),
DIMENSION(:,:),
ALLOCATABLE :: tempArr
4743 IF(finalupdate)
THEN
4749 CALL mpalloc(temparr,two,oldlength,
'INONE: temp array')
4769 IF(j == 0)
EXIT inner
4770 IF(j == i) cycle outer
4774 IF(.NOT.finalupdate)
RETURN
4779 WRITE(
lunlog,*)
'INONE: array reduced to',newlength,
' words'
4793 INTEGER(mpi),
INTENT(IN) :: n
4794 INTEGER(mpi) :: nprime
4795 INTEGER(mpi) :: nsqrt
4800 IF(mod(nprime,2) == 0) nprime=nprime+1
4803 nsqrt=int(sqrt(real(nprime,
mps)),
mpi)
4805 IF(i*(nprime/i) == nprime) cycle outer
4827 INTEGER(mpi) :: idum
4829 INTEGER(mpi) :: indab
4830 INTEGER(mpi) :: itgbi
4831 INTEGER(mpi) :: itgbl
4832 INTEGER(mpi) :: ivgbi
4835 INTEGER(mpi) :: nc31
4837 INTEGER(mpi) :: nwrd
4838 INTEGER(mpi) :: inone
4843 INTEGER(mpl) :: length
4847 WRITE(
lunlog,*)
'LOOP1: starting'
4866 WRITE(
lunlog,*)
'LOOP1: reading data files'
4876 WRITE(*,*)
'Read all binary data files:'
4878 CALL hmpldf(1,
'Number of words/record in binary file')
4879 CALL hmpdef(8,0.0,60.0,
'not_stored data per record')
4906 WRITE(
lunlog,*)
'LOOP1: reading data files again'
4918 CALL peend(21,
'Aborted, no labels/parameters defined')
4919 stop
'LOOP1: no labels/parameters defined'
4923 ' is total number NTGB of labels/parameters'
4925 CALL hmpldf(2,
'Number of entries per label')
4929 IF(
nhistp /= 0)
CALL hmprnt(2)
4963 WRITE(
lunlog,*)
'LOOP1:',
npresg,
' is number of pre-sigmas'
4964 WRITE(*,*)
'LOOP1:',
npresg,
' is number of pre-sigmas'
4965 IF(
npresg == 0)
WRITE(*,*)
'Warning: no pre-sigmas defined'
4983 WRITE(
lunlog,*)
'LOOP1:',
nvgb,
' is number NVGB of variable parameters'
4999 WRITE(*,112)
' Default pre-sigma =',
regpre, &
5000 ' (if no individual pre-sigma defined)'
5001 WRITE(*,*)
'Pre-sigma factor is',
regula
5004 WRITE(*,*)
'No regularization will be done'
5006 WRITE(*,*)
'Regularization will be done, using factor',
regula
5010 CALL peend(22,
'Aborted, no variable global parameters')
5011 stop
'... no variable global parameters'
5018 IF(presg > 0.0)
THEN
5020 ELSE IF(presg == 0.0.AND.
regpre > 0.0)
THEN
5021 prewt=1.0/real(
regpre**2,mpd)
5037 WRITE(*,101)
'NTGB',
ntgb,
'total number of parameters'
5038 WRITE(*,101)
'NVGB',
nvgb,
'number of variable parameters'
5045 WRITE(*,101)
' NREC',
nrec,
'number of records'
5046 IF (
nrecd > 0)
WRITE(*,101)
' NRECD',
nrec,
'number of records containing doubles'
5047 WRITE(*,101)
'MREQENF',
mreqenf,
'required number of entries (from binary files)'
5049 WRITE(*,101)
'ITEREN',
iteren,
'iterate cut for parameters with less entries'
5050 WRITE(*,101)
'MREQENA',
mreqena,
'required number of entries (from accepted fits)'
5051 IF (
mreqpe > 1)
WRITE(*,101) &
5052 'MREQPE',
mreqpe,
'required number of pair entries'
5053 IF (
msngpe >= 1)
WRITE(*,101) &
5054 'MSNGPE',
msngpe,
'max pair entries single prec. storage'
5055 WRITE(*,101)
'NTGB',
ntgb,
'total number of parameters'
5056 WRITE(*,101)
'NVGB',
nvgb,
'number of variable parameters'
5059 WRITE(*,*)
'Global parameter labels:'
5066 mqi=((mqi-20)/20)*20+1
5074 WRITE(8,101)
' NREC',
nrec,
'number of records'
5075 IF (
nrecd > 0)
WRITE(8,101)
' NRECD',
nrec,
'number of records containing doubles'
5076 WRITE(8,101)
'MREQENF',
mreqenf,
'required number of entries (from binary files)'
5078 WRITE(8,101)
'ITEREN',
iteren,
'iterate cut for parameters with less entries'
5079 WRITE(8,101)
'MREQENA',
mreqena,
'required number of entries (from accepted fits)'
5081 WRITE(
lunlog,*)
'LOOP1: ending'
5085101
FORMAT(1x,a8,
' =',i10,
' = ',a)
5101 INTEGER(mpi) :: ibuf
5103 INTEGER(mpi) :: indab
5104 INTEGER(mpi) :: inder
5105 INTEGER(mpi) :: isfrst
5106 INTEGER(mpi) :: islast
5112 INTEGER(mpi) :: nc31
5114 INTEGER(mpi) :: nlow
5116 INTEGER(mpi) :: nwrd
5119 INTEGER(mpl) :: length
5120 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: newCounter
5129 WRITE(
lunlog,*)
'LOOP1: iterating'
5131 WRITE(*,*)
'LOOP1: iterating'
5134 CALL mpalloc(newcounter,length,
'new entries counter')
5158 CALL isjajb(nst,ist,ja,jb,jsp)
5159 IF(ja == 0.AND.jb == 0)
EXIT
5165 IF(ij == -2) nlow=nlow+1
5170 newcounter(ij)=newcounter(ij)+1
5199 WRITE(
lunlog,*)
'LOOP1:',
nvgb,
' is number NVGB of variable parameters'
5230 INTEGER(mpi) :: ibuf
5231 INTEGER(mpi) :: icgb
5232 INTEGER(mpi) :: iext
5233 INTEGER(mpi) :: ihis
5236 INTEGER(mpi) :: inder
5237 INTEGER(mpi) :: ioff
5238 INTEGER(mpi) :: iproc
5239 INTEGER(mpi) :: irecmm
5240 INTEGER(mpi) :: isfrst
5241 INTEGER(mpi) :: islast
5243 INTEGER(mpi) :: itgbi
5244 INTEGER(mpi) :: itgbij
5245 INTEGER(mpi) :: itgbik
5246 INTEGER(mpi) :: ivgbij
5247 INTEGER(mpi) :: ivgbik
5251 INTEGER(mpi) :: jext
5252 INTEGER(mpi) :: jcgb
5254 INTEGER(mpi) :: joff
5256 INTEGER(mpi) :: kfile
5258 INTEGER(mpi) :: label
5261 INTEGER(mpi) :: maeqnf
5262 INTEGER(mpi) :: naeqna
5263 INTEGER(mpi) :: naeqnf
5264 INTEGER(mpi) :: naeqng
5265 INTEGER(mpi) :: nc31
5266 INTEGER(mpi) :: ncachd
5267 INTEGER(mpi) :: ncachi
5268 INTEGER(mpi) :: ncachr
5271 INTEGER(mpi) :: ndfmax
5272 INTEGER(mpi) :: nfixed
5273 INTEGER(mpi) :: nggd
5274 INTEGER(mpi) :: nggi
5275 INTEGER(mpi) :: nmatmo
5276 INTEGER(mpi) :: noff
5278 INTEGER(mpi) :: nrecf
5279 INTEGER(mpi) :: nrecmm
5281 INTEGER(mpi) :: nwrd
5282 INTEGER(mpi) :: inone
5291 INTEGER(mpl):: noff8
5292 INTEGER(mpl):: ndimbi
5293 INTEGER(mpl):: ndimsa(4)
5295 INTEGER(mpl):: matsiz(2)
5296 INTEGER(mpl):: matwords
5297 INTEGER(mpl):: length
5300 INTEGER(mpl),
PARAMETER :: two=2
5301 INTEGER(mpi) :: maxGlobalPar = 0
5302 INTEGER(mpi) :: maxLocalPar = 0
5303 INTEGER(mpi) :: maxEquations = 0
5306 SUBROUTINE ndbits(ndims,ncmprs,nsparr,ihst)
5308 INTEGER(mpl),
DIMENSION(4),
INTENT(OUT) :: ndims
5309 INTEGER(mpi),
DIMENSION(:),
INTENT(OUT) :: ncmprs
5310 INTEGER(mpl),
DIMENSION(:,:),
INTENT(OUT) :: nsparr
5311 INTEGER(mpi),
INTENT(IN) :: ihst
5313 SUBROUTINE spbits(nsparr,nsparc,ncmprs)
5315 INTEGER(mpl),
DIMENSION(:,:),
INTENT(IN) :: nsparr
5316 INTEGER(mpi),
DIMENSION(:),
INTENT(OUT) :: nsparc
5317 INTEGER(mpi),
DIMENSION(:),
INTENT(IN) :: ncmprs
5321 INTEGER(mpi),
DIMENSION(:),
INTENT(OUT) :: npair
5335 WRITE(
lunlog,*)
'LOOP2: starting'
5384 WRITE(
lunlog,*)
'LOOP2: start event reading'
5433 WRITE(*,*)
'Record number ',
nrec,
' from file ',kfile
5434 IF (wgh /= 1.0)
WRITE(*,*)
' weight ',wrec
5438 CALL isjajb(nst,ist,ja,jb,jsp)
5442 IF(nda ==
mdebg2+1)
WRITE(*,*)
'... and more data'
5446 WRITE(*,*) nda,
' Measured value =',glder(ja),
' +- ',glder(jb)
5447 WRITE(*,*)
'Local derivatives:'
5448 WRITE(*,107) (inder(ja+j),glder(ja+j),j=1,jb-ja-1)
5449107
FORMAT(6(i3,g12.4))
5451 WRITE(*,*)
'Global derivatives:'
5454108
FORMAT(3i11,g12.4)
5457 WRITE(*,*)
'total_par_label __label__ var_par_index derivative'
5471 CALL isjajb(nst,ist,ja,jb,jsp)
5472 IF(ja == 0.AND.jb == 0)
EXIT
5480 rerr =real(glder(jb),mpd)
5502 CALL inbmap(ij,inder(jb+k))
5518 IF (nfixed > 0) naeqnf=naeqnf+1
5521 IF(ja /= 0.AND.jb /= 0)
THEN
5530 IF (naeqnf > maeqnf) nrecf=nrecf+1
5534 maxglobalpar=max(
nagbn,maxglobalpar)
5535 maxlocalpar=max(
nalcn,maxlocalpar)
5536 maxequations=max(
naeqn,maxequations)
5539 dstat(1)=dstat(1)+real((nwrd+2)*2,mpd)
5540 dstat(2)=dstat(2)+real(
nagbn+2,mpd)
5580 (irecmm >= nrecmm.OR.irecmm ==
mxrec))
THEN
5581 IF (nmatmo == 0)
THEN
5583 WRITE(*,*)
'Monitoring of sparse matrix construction'
5584 WRITE(*,*)
' records ........ off-diagonal elements ', &
5585 '....... compression memory'
5586 WRITE(*,*)
' non-zero used(double) used', &
5591 gbc=1.0e-9*real((mpi*ndimsa(2)+mpd*ndimsa(3)+mps*ndimsa(4))/mpi*(bit_size(1_mpi)/8),mps)
5592 gbu=1.0e-9*real(((mpi+mpd)*(ndimsa(3)+ndimsa(4)))/mpi*(bit_size(1_mpi)/8),mps)
5594 WRITE(*,1177) irecmm,ndimsa(1),ndimsa(3),ndimsa(4),cpr,gbc
55951177
FORMAT(i9,3i13,f10.2,f11.6)
5596 DO WHILE(irecmm >= nrecmm)
5615 WRITE(
lunlog,*)
'LOOP2: event reading ended - end of data'
5617 dstat(k)=dstat(k)/real(
nrec,mpd)
5674 IF(ivgbij > 0.AND.ivgbik > 0)
THEN
5676 IF (
mprint > 1)
WRITE(*,*)
'add index pair ',ivgbij,ivgbik
5698 IF (
nagb >= 65536)
THEN
5699 noff=int(noff8/1000,mpi)
5709 CALL hmpdef(ihis,0.0,real(
mhispe,mps),
'NDBITS: #off-diagonal elements')
5717 ndgn=ndimsa(3)+ndimsa(4)
5718 matwords=ndimsa(2)+length
5721 IF (
nhistp /= 0)
CALL hmprnt(ihis)
5733 IF (
fcache(2) == 0.0)
THEN
5735 fcache(2)=real(dstat(2),mps)
5736 fcache(3)=real(dstat(3),mps)
5757 ncachd=
ncache-ncachr-ncachi
5785 WRITE(lu,101)
'NTGB',
ntgb,
'total number of parameters'
5786 WRITE(lu,102)
'(all parameters, appearing in binary files)'
5787 WRITE(lu,101)
'NVGB',
nvgb,
'number of variable parameters'
5788 WRITE(lu,102)
'(appearing in fit matrix/vectors)'
5789 WRITE(lu,101)
'NAGB',
nagb,
'number of all parameters'
5790 WRITE(lu,102)
'(including Lagrange multiplier or reduced)'
5791 WRITE(lu,101)
'NFGB',
nfgb,
'number of fit parameters'
5792 WRITE(lu,101)
'MBANDW',
mbandw,
'band width of band matrix'
5793 WRITE(lu,102)
'(if =0, no band matrix)'
5794 IF (
nagb >= 65536)
THEN
5795 WRITE(lu,101)
'NOFF/K',noff,
'max number of off-diagonal elements'
5797 WRITE(lu,101)
'NOFF',noff,
'max number of off-diagonal elements'
5800 IF (
nagb >= 65536)
THEN
5801 WRITE(lu,101)
'NDGN/K',ndgn/1000,
'actual number of off-diagonal elements'
5803 WRITE(lu,101)
'NDGN',ndgn,
'actual number of off-diagonal elements'
5806 WRITE(lu,101)
'NCGB',
ncgb,
'number of constraints'
5807 WRITE(lu,101)
'NAGBN',
nagbn,
'max number of global parameters in an event'
5808 WRITE(lu,101)
'NALCN',
nalcn,
'max number of local parameters in an event'
5809 WRITE(lu,101)
'NAEQN',
naeqn,
'max number of equations in an event'
5811 WRITE(lu,101)
'NAEQNA',naeqna,
'number of equations'
5812 WRITE(lu,101)
'NAEQNG',naeqng, &
5813 'number of equations with global derivatives'
5814 WRITE(lu,101)
'NAEQNF',naeqnf, &
5815 'number of equations with fixed global derivatives'
5816 WRITE(lu,101)
'NRECF',nrecf, &
5817 'number of records with fixed global derivatives'
5820 WRITE(lu,101)
'NCACHE',
ncache,
'number of words for caching'
5821 WRITE(lu,111) (
fcache(k)*100.0,k=1,3)
5822111
FORMAT(22x,
'cache splitting ',3(f6.1,
' %'))
5827 WRITE(lu,*)
'Solution method and matrix-storage mode:'
5829 WRITE(lu,*)
' METSOL = 1: matrix inversion'
5830 ELSE IF(
metsol == 2)
THEN
5831 WRITE(lu,*)
' METSOL = 2: diagonalization'
5832 ELSE IF(
metsol == 3)
THEN
5833 WRITE(lu,*)
' METSOL = 3: MINRES (rtol',
mrestl,
')'
5834 ELSE IF(
metsol == 4)
THEN
5835 WRITE(lu,*)
' METSOL = 4: MINRES-QLP (rtol',
mrestl,
')'
5836 ELSE IF(
metsol == 5)
THEN
5837 WRITE(lu,*)
' METSOL = 5: GMRES'
5839 WRITE(lu,*)
' with',
mitera,
' iterations'
5841 WRITE(lu,*)
' MATSTO = 1: symmetric matrix, ',
'(n*n+n)/2 elements'
5842 ELSE IF(
matsto == 2)
THEN
5843 WRITE(lu,*)
' MATSTO = 2: sparse matrix'
5845 IF(
mextnd>0)
WRITE(lu,*)
' with extended storage'
5846 IF(
dflim /= 0.0)
THEN
5847 WRITE(lu,103)
'Convergence assumed, if expected dF <',
dflim
5851 WRITE(lu,*)
'Constraints handled by elimination'
5853 WRITE(lu,*)
'Constraints handled by Lagrange multipliers'
5877105
FORMAT(
' Constants C1, C2 for Wolfe conditions:',g12.4,
', ',g12.4)
588132 matsiz(1)=int8(
nagb)*int8(
nagb+1)/2
5884 matsiz(1)=ndimsa(3)+
nagb
5889 matwords=matwords+matsiz(1)*2+matsiz(2)
5895 IF (matwords < 250000)
THEN
5896 WRITE(*,*)
'Size of global matrix: < 1 MB'
5898 WRITE(*,*)
'Size of global matrix:',int(real(matwords,mps)*4.0e-6,mpi),
' MB'
5904 WRITE(
lunlog,*)
' Cut values of Chi^2/Ndf and Chi2,'
5905 WRITE(
lunlog,*)
' corresponding to 2 and 3 standard deviations'
5906 WRITE(
lunlog,*)
' Ndf Chi^2/Ndf(2) Chi^2(2) ', &
5907 ' Chi^2/Ndf(3) Chi^2(3)'
5910 IF(ndf >
naeqn)
EXIT
5913 ELSE IF(ndf < 20)
THEN
5915 ELSE IF(ndf < 100)
THEN
5917 ELSE IF(ndf < 200)
THEN
5924 WRITE(
lunlog,106) ndf,chin2,chin2*real(ndf,mps),chin3, chin3*real(ndf,mps)
5927 WRITE(
lunlog,*)
'LOOP2: ending'
5930101
FORMAT(1x,a8,
' =',i10,
' = ',a)
5932103
FORMAT(1x,a,g12.4)
5933106
FORMAT(i6,2(3x,f9.3,f12.1,3x))
5947 INTEGER(mpi) :: imed
5950 INTEGER(mpi) :: nent
5951 INTEGER(mpi),
DIMENSION(measBins) :: isuml
5952 INTEGER(mpi),
DIMENSION(measBins) :: isums
5956 INTEGER(mpl) :: ioff
5959 DATA lfirst /.true./
5972 WRITE(
lunmon,
'(A)')
'*** Normalized residuals grouped by first global label (per local fit cycle) ***'
5974 WRITE(
lunmon,
'(A)')
'*** Pulls grouped by first global label (per local fit cycle) ***'
5976 WRITE(
lunmon,
'(A)')
'! LFC Label Entries Median RMS(MAD) <error>'
5993 IF (2*isuml(j) > nent)
EXIT
5997 IF (isuml(j) > isuml(j-1)) amed=amed+real(nent-2*isuml(j-1),mps)/real(2*isuml(j)-2*isuml(j-1),mps)
6010 isums(j)=isums(j)+isums(j-1)
6014 IF (2*isums(j) > nent)
EXIT
6017 IF (isums(j) > isums(j-1)) amad=amad+real(nent-2*isums(j-1),mps)/real(2*isums(j)-2*isums(j-1),mps)
6026110
FORMAT(i5,2i10,3g14.5)
6040 INTEGER(mpi) :: ncon
6042 INTEGER(mpl),
INTENT(IN) :: msize(2)
6044 INTEGER(mpl) :: length
6055 CALL mpalloc(
aux,length,
' local fit scratch array: aux')
6056 CALL mpalloc(
vbnd,length,
' local fit scratch array: vbnd')
6057 CALL mpalloc(
vbdr,length,
' local fit scratch array: vbdr')
6060 CALL mpalloc(
vbk,length,
' local fit scratch array: vbk')
6063 CALL mpalloc(
vzru,length,
' local fit scratch array: vzru')
6064 CALL mpalloc(
scdiag,length,
' local fit scratch array: scdiag')
6065 CALL mpalloc(
scflag,length,
' local fit scratch array: scflag')
6066 CALL mpalloc(
ibandh,2*length,
' local fit band width hist.: ibandh')
6092 length=
nvgb+ncon*
nvgb+(ncon*ncon+ncon)/2
6134 INTEGER(mpl) :: ioff1
6137 INTEGER(mpi) :: nrank
6176 IF(
nfgb /= nrank)
THEN
6177 WRITE(*,*)
'Warning: the rank defect of the symmetric',
nfgb, &
6178 '-by-',
nfgb,
' matrix is ',
nfgb-nrank,
' (should be zero).'
6179 WRITE(lun,*)
'Warning: the rank defect of the symmetric',
nfgb, &
6180 '-by-',
nfgb,
' matrix is ',
nfgb-nrank,
' (should be zero).'
6183 WRITE(*,*)
' --> enforcing SUBITO mode'
6184 WRITE(lun,*)
' --> enforcing SUBITO mode'
6186 ELSE IF(
ndefec == 0)
THEN
6187 WRITE(lun,*)
'No rank defect of the symmetric matrix'
6212 INTEGER(mpi) :: iast
6213 INTEGER(mpi) :: idia
6214 INTEGER(mpi) :: imin
6215 INTEGER(mpl) :: ioff1
6218 INTEGER(mpi) :: nmax
6219 INTEGER(mpi) :: nmin
6220 INTEGER(mpi) :: ntop
6274 DO WHILE(ntop < nmax)
6278 CALL hmpdef(7,real(nmin,mps),real(ntop,mps),
'log10 of positive eigenvalues')
6282 CALL hmpent(7,evalue)
6285 IF(
nhistp /= 0)
CALL hmprnt(7)
6289 CALL gmpdef(3,2,
'low-value end of eigenvalues')
6292 CALL gmpxy(3,real(i,mps),evalue)
6294 IF(
nhistp /= 0)
CALL gmprnt(3)
6305 WRITE(lun,*)
'The first (largest) eigenvalues ...'
6308 WRITE(lun,*)
'The last eigenvalues ... up to',
nvgb
6312 WRITE(lun,*)
'The eigenvalues from',
nvgb+1,
' to',
nagb
6316 WRITE(lun,*)
'Log10 + 3 of ',
nagb,
' eigenvalues in decreasing',
' order'
6317 WRITE(lun,*)
'(for Eigenvalue < 0.001 the value 0.0 is shown)'
6320 'printed for negative eigenvalues'
6323 WRITE(lun,*)
nvgb,
' significances: insignificant if ', &
6324 'compatible with N(0,1)'
6353 INTEGER(mpl) :: ioff1
6354 INTEGER(mpl) :: ioff2
6390 INTEGER(mpi) :: istop
6392 INTEGER(mpi) :: itnlim
6394 INTEGER(mpi) :: nout
6395 INTEGER(mpi) :: nrkd
6396 INTEGER(mpi) :: nrkd2
6442 globalcorrections, itnlim, nout, rtol, istop, itn, anorm, acond, rnorm, arnorm, ynorm)
6448 WRITE(lun,*)
'MMINRS: EQUDEC ended ', nrkd, nrkd2
6451 globalcorrections, itnlim, nout, rtol, istop, itn, anorm, acond, rnorm, arnorm, ynorm)
6454 globalcorrections, itnlim, nout, rtol, istop, itn, anorm, acond, rnorm, arnorm, ynorm)
6468 IF (
istopa == 0) print *,
'MINRES: istop=0, exact solution x=0.'
6483 INTEGER(mpi) :: istop
6485 INTEGER(mpi) :: itnlim
6487 INTEGER(mpi) :: nout
6488 INTEGER(mpi) :: nrkd
6489 INTEGER(mpi) :: nrkd2
6504 mxxnrm = real(
nagb,mpd)/sqrt(epsilon(mxxnrm))
6506 trcond = 1.0_mpd/epsilon(trcond)
6507 ELSE IF(
mrmode == 2)
THEN
6537 itnlim=itnlim, rtol=rtol, maxxnorm=mxxnrm, trancond=trcond, &
6544 WRITE(lun,*)
'MMINRS: EQUDEC ended ', nrkd, nrkd2
6548 itnlim=itnlim, rtol=rtol, maxxnorm=mxxnrm, trancond=trcond, &
6552 itnlim=itnlim, rtol=rtol, maxxnorm=mxxnrm, trancond=trcond, &
6567 IF (
istopa == 3) print *,
'MINRES: istop=0, exact solution x=0.'
6583 INTEGER(mpi),
INTENT(IN) :: n
6584 REAL(mpd),
INTENT(IN) :: x(n)
6585 REAL(mpd),
INTENT(OUT) :: y(n)
6604 INTEGER(mpi),
INTENT(IN) :: n
6605 REAL(mpd),
INTENT(IN) :: x(n)
6606 REAL(mpd),
INTENT(OUT) :: y(n)
6639 REAL,
DIMENSION(2) :: ta
6641 INTEGER(mpi) :: iact
6642 INTEGER(mpi) :: iagain
6644 INTEGER(mpi) :: info
6645 INTEGER(mpl) :: ioff
6646 INTEGER(mpi) :: itgbi
6647 INTEGER(mpi) :: ivgbi
6648 INTEGER(mpi) :: jcalcm
6650 INTEGER(mpi) :: labelg
6651 INTEGER(mpi) :: litera
6652 INTEGER(mpi) :: lrej
6654 INTEGER(mpi) :: lunp
6655 INTEGER(mpi) :: minf
6656 INTEGER(mpi) :: mrati
6658 INTEGER(mpi) :: nfaci
6659 INTEGER(mpi) :: nloopsol
6660 INTEGER(mpi) :: nrati
6661 INTEGER(mpi) :: nrej
6662 INTEGER(mpi) :: nsol
6663 INTEGER(mpi) :: inone
6677 CHARACTER (LEN=7) :: cratio
6678 CHARACTER (LEN=7) :: cfacin
6679 CHARACTER (LEN=7) :: crjrat
6691 WRITE(lunp,*)
'Solution algorithm: '
6692 WRITE(lunp,121)
'=================================================== '
6695 WRITE(lunp,121)
'solution method:',
'matrix inversion'
6696 ELSE IF(
metsol == 2)
THEN
6697 WRITE(lunp,121)
'solution method:',
'diagonalization'
6698 ELSE IF(
metsol == 3)
THEN
6699 WRITE(lunp,121)
'solution method:',
'minres (Paige/Saunders)'
6700 ELSE IF(
metsol == 4)
THEN
6701 WRITE(lunp,121)
'solution method:',
'minres-qlp (Choi/Paige/Saunders)'
6703 WRITE(lunp,121)
' ',
' using QR factorization'
6704 ELSE IF(
mrmode == 2)
THEN
6705 WRITE(lunp,121)
' ',
' using QLP factorization'
6707 WRITE(lunp,121)
' ',
' using QR and QLP factorization'
6708 WRITE(lunp,123)
'transition condition',
mrtcnd
6710 ELSE IF(
metsol == 5)
THEN
6711 WRITE(lunp,121)
'solution method:', &
6712 'gmres (generalized minimzation of residuals)'
6714 WRITE(lunp,123)
'convergence limit at Delta F=',
dflim
6715 WRITE(lunp,122)
'maximum number of iterations=',
mitera
6718 WRITE(lunp,122)
'matrix recalculation up to ',
matrit,
'. iteration'
6722 WRITE(lunp,121)
'matrix storage:',
'full'
6723 ELSE IF(
matsto == 2)
THEN
6724 WRITE(lunp,121)
'matrix storage:',
'sparse'
6726 WRITE(lunp,122)
'pre-con band-width parameter=',
mbandw
6728 WRITE(lunp,121)
'pre-conditioning:',
'default'
6730 WRITE(lunp,121)
'pre-conditioning:',
'none!'
6733 WRITE(lunp,121)
'pre-conditioning=',
'skyline-matrix (rank preserving)'
6735 WRITE(lunp,121)
'pre-conditioning=',
'band-matrix'
6740 WRITE(lunp,121)
'using pre-sigmas:',
'no'
6743 WRITE(lunp,124)
'pre-sigmas defined for', &
6744 REAL(100*npresg,mps)/
REAL(nvgb,mps),
' % of variable parameters'
6745 WRITE(lunp,123)
'default pre-sigma=',
regpre
6748 WRITE(lunp,121)
'regularization:',
'no'
6750 WRITE(lunp,121)
'regularization:',
'yes'
6751 WRITE(lunp,123)
'regularization factor=',
regula
6755 WRITE(lunp,121)
'Chi square cut equiv 3 st.dev applied'
6756 WRITE(lunp,123)
'... in first iteration with factor',
chicut
6757 WRITE(lunp,123)
'... in second iteration with factor',
chirem
6758 WRITE(lunp,121)
' (reduced by sqrt in next iterations)'
6761 WRITE(lunp,121)
'Scaling of measurement errors applied'
6762 WRITE(lunp,123)
'... factor for "global" measuements',
dscerr(1)
6763 WRITE(lunp,123)
'... factor for "local" measuements',
dscerr(2)
6766 WRITE(lunp,122)
'Down-weighting of outliers in',
lhuber,
' iterations'
6767 WRITE(lunp,123)
'Cut on downweight fraction',
dwcut
6771121
FORMAT(1x,a40,3x,a)
6772122
FORMAT(1x,a40,3x,i0,a)
6773123
FORMAT(1x,a40,2x,e9.2)
6774124
FORMAT(1x,a40,3x,f5.1,a)
6784 stepl =real(stp,mps)
6796 ELSE IF(
metsol == 2)
THEN
6799 ELSE IF(
metsol == 3)
THEN
6802 ELSE IF(
metsol == 4)
THEN
6805 ELSE IF(
metsol == 5)
THEN
6814 WRITE(
lunlog,*)
'Checking feasibility of parameters:'
6815 WRITE(*,*)
'Checking feasibility of parameters:'
6819 WRITE(*,*)
' parameters are made feasible'
6821 WRITE(
lunlog,*)
' parameters are made feasible'
6823 WRITE(*,*)
' parameters are feasible (i.e. satisfy constraints)'
6824 WRITE(
lunlog,*)
' parameters are feasible (i.e. satisfy constraints)'
6832 WRITE(*,*)
'Reading files and accumulating vectors/matrices ...'
6848 IF(jcalcm+1 /= 0)
THEN
6855 CALL gmpxyd(1,real(
nloopn,mps),real(
fvalue,mps),0.5,0.)
6857 IF(
iterat /= litera)
THEN
6869 CALL gmpxyd(1,real(
nloopn,mps),real(
fvalue,mps),0.5,0.)
6878 IF (iabs(jcalcm) <= 1)
THEN
6893 WRITE(*,*)
'Data rejected in previous loop: '
6895 nrejec(0),
' (rank deficit/NaN) ',
nrejec(1),
' (Ndf=0) ', &
6897 WRITE(*,*)
'Too many rejects (>33.3%) - stop'
6898 CALL peend(26,
'Aborted, too many rejects')
6904 IF(abs(
icalcm) == 1)
THEN
6915 ELSE IF(
metsol == 2)
THEN
6917 ELSE IF(
metsol == 3)
THEN
6919 ELSE IF(
metsol == 4)
THEN
6921 ELSE IF(
metsol == 5)
THEN
6922 WRITE(*,*)
'... reserved for GMRES (not yet!)'
6946 angras=real(db/sqrt(db1*db2),mps)
6958 ELSE IF(
metsol == 2)
THEN
6961 ELSE IF(
metsol == 3)
THEN
6964 ELSE IF(
metsol == 4)
THEN
6967 ELSE IF(
metsol == 5)
THEN
6974 IF(db <= -16.0_mpd*sqrt(max(db1,db2))*epsilon(db))
THEN
6975 WRITE(*,*)
'Function not decreasing:',db
6976 IF(db > -1.0e-3_mpd)
THEN
6978 IF (iagain <= 1)
THEN
6979 WRITE(*,*)
'... again matrix calculation'
6983 WRITE(*,*)
'... aborting iterations'
6987 WRITE(*,*)
'... stopping iterations'
7010 IF (
nloopn == nloopsol)
THEN
7026 WRITE(*,*)
'Result vector containes ', nan,
' NaNs - stop'
7027 CALL peend(25,
'Aborted, result vector contains NaNs')
7034 WRITE(*,*)
'Subito! Exit after first step.'
7039 WRITE(*,*)
'INFO=0 should not happen (line search input err)'
7040 IF (iagain <= 0)
THEN
7045 IF(info < 0 .OR.
nloopn == nloopsol) cycle
7050 IF(iact /= 0.OR.
chicut > 1.0)
THEN
7067 WRITE(*,*)
'Data rejected in last loop: '
7069 nrejec(0),
' (rank deficit/NaN) ',
nrejec(1),
' (Ndf=0) ', &
7079 catio=real(dratio,mps)
7083 mrati=nint(100.0*catio,mpi)
7087 IF (
nfilw <= 0)
THEN
7088 WRITE(lunp,*)
'Sum(Chi^2)/Sum(Ndf) =',
fvalue
7090 WRITE(lunp,*)
' =',dratio
7092 WRITE(lunp,*)
'Sum(W*Chi^2)/Sum(Ndf)/<W> =',
fvalue
7094 WRITE(lunp,*)
' /',dwmean
7095 WRITE(lunp,*)
' =',dratio
7099 ' with correction for down-weighting ',catio
7120 nrati=nint(10000.0*real(nrej,mps)/real(
nrecal,mps),mpi)
7121 WRITE(crjrat,197) 0.01_mpd*real(nrati,mpd)
7122 nfaci=nint(100.0*sqrt(catio),mpi)
7124 WRITE(cratio,197) 0.01_mpd*real(mrati,mpd)
7125 WRITE(cfacin,197) 0.01_mpd*real(nfaci,mpd)
7128 IF(mrati < 90.OR.mrati > 110) warner=.true.
7129 IF(nrati > 100) warner=.true.
7130 IF(
ncgbe /= 0) warner=.true.
7132 IF(
nalow /= 0) warners=.true.
7134 IF(
nmiss1 /= 0) warnerss=.true.
7135 IF(iagain /= 0) warnerss=.true.
7136 IF(
ndefec /= 0) warnerss=.true.
7138 IF(warner.OR.warners.OR.warnerss)
THEN
7141 WRITE(*,199)
'WarningWarningWarningWarningWarningWarningWarningWarningWar'
7142 WRITE(*,199)
'arningWarningWarningWarningWarningWarningWarningWarningWarn'
7143 WRITE(*,199)
'rningWarningWarningWarningWarningWarningWarningWarningWarni'
7144 WRITE(*,199)
'ningWarningWarningWarningWarningWarningWarningWarningWarnin'
7145 WRITE(*,199)
'ingWarningWarningWarningWarningWarningWarningWarningWarning'
7146 WRITE(*,199)
'ngWarningWarningWarningWarningWarningWarningWarningWarningW'
7147 WRITE(*,199)
'gWarningWarningWarningWarningWarningWarningWarningWarningWa'
7149 IF(mrati < 90.OR.mrati > 110)
THEN
7151 WRITE(*,*)
' Chi^2/Ndf = ',cratio,
' (should be close to 1)'
7152 WRITE(*,*)
' => multiply all input standard ', &
7153 'deviations by factor',cfacin
7156 IF(nrati > 100)
THEN
7158 WRITE(*,*)
' Fraction of rejects =',crjrat,
' %', &
7159 ' (should be far below 1 %)'
7160 WRITE(*,*)
' => please provide correct mille data'
7163 IF(iagain /= 0)
THEN
7165 WRITE(*,*)
' Matrix not positiv definite '// &
7166 '(function not decreasing)'
7167 WRITE(*,*)
' => please provide correct mille data'
7172 WRITE(*,*)
' Rank defect =',
ndefec, &
7173 ' for global matrix, should be 0'
7174 WRITE(*,*)
' => please provide correct mille data'
7179 WRITE(*,*)
' Rank defect =',
nmiss1, &
7180 ' for constraint equations, should be 0'
7181 WRITE(*,*)
' => please correct constraint definition'
7186 WRITE(*,*)
' Number of empty constraints =',
ncgbe,
', should be 0'
7187 WRITE(*,*)
' => please check constraint definition, mille data'
7192 WRITE(*,*)
' Possible rank defects =',
nalow, &
7193 ' for global vector (too few entries)'
7194 WRITE(*,*)
' => please check mille data and ENTRIES cut'
7198 WRITE(*,199)
'WarningWarningWarningWarningWarningWarningWarningWarningWar'
7199 WRITE(*,199)
'arningWarningWarningWarningWarningWarningWarningWarningWarn'
7200 WRITE(*,199)
'rningWarningWarningWarningWarningWarningWarningWarningWarni'
7201 WRITE(*,199)
'ningWarningWarningWarningWarningWarningWarningWarningWarnin'
7202 WRITE(*,199)
'ingWarningWarningWarningWarningWarningWarningWarningWarning'
7203 WRITE(*,199)
'ngWarningWarningWarningWarningWarningWarningWarningWarningW'
7204 WRITE(*,199)
'gWarningWarningWarningWarningWarningWarningWarningWarningWa'
7215 ELSE IF(
metsol == 2)
THEN
7221 IF(labelg == 0) cycle
7225 IF(ivgbi < 0) ivgbi=0
7226 IF(ivgbi == 0) cycle
7235 ELSE IF(
metsol == 5)
THEN
7254 CALL peend(3,
'Ended with severe warnings (bad global matrix)')
7255 ELSE IF (warners)
THEN
7256 CALL peend(2,
'Ended with severe warnings (insufficient measurements)')
7257 ELSE IF (warner)
THEN
7258 CALL peend(1,
'Ended with warnings (bad measurements)')
7260 CALL peend(0,
'Ended normally')
7263102
FORMAT(
' Call FEASIB with cut=',g10.3)
7290 INTEGER(mpi) :: iargc
7293 INTEGER(mpi) :: ierrf
7295 INTEGER(mpi) :: ifilb
7296 INTEGER(mpi) :: ioff
7297 INTEGER(mpi) :: iopt
7299 INTEGER(mpi) :: iosum
7304 INTEGER(mpi) :: nline
7305 INTEGER(mpi) :: npat
7306 INTEGER(mpi) :: ntext
7309 INTEGER(mpi) :: nums
7310 INTEGER(mpi) :: nufile
7311 INTEGER(mpi) :: lenfileInfo
7312 INTEGER(mpi) :: lenFileNames
7313 INTEGER(mpi) :: matint
7314 INTEGER(mpi),
DIMENSION(:,:),
ALLOCATABLE :: vecfileInfo
7315 INTEGER(mpi),
DIMENSION(:,:),
ALLOCATABLE :: tempArray
7316 INTEGER(mpl) :: rows
7317 INTEGER(mpl) :: cols
7318 INTEGER(mpl) :: newcols
7319 INTEGER(mpl) :: length
7321 CHARACTER (LEN=1024) :: text
7322 CHARACTER (LEN=1024) :: fname
7323 CHARACTER (LEN=14) :: bite(3)
7324 CHARACTER (LEN=32) :: keystx
7325 REAL(mpd) :: dnum(100)
7327 DATA bite/
'C_binary',
'text ',
'Fortran_binary'/
7342 WRITE(*,*)
'Command line options: '
7343 WRITE(*,*)
'--------------------- '
7346 CALL rltext(text,ia,ib,nab)
7347 WRITE(*,101) i,text(1:nab)
7348 IF(text(ia:ia) /=
'-')
THEN
7349 nu=nufile(text(ia:ib))
7352 WRITE(*,*)
'Second text file in command line - stop'
7353 CALL peend(12,
'Aborted, second text file in command line')
7359 WRITE(*,*)
'Open error for file:',text(ia:ib),
' - stop'
7360 CALL peend(16,
'Aborted, open error for file')
7364 IF(index(text(ia:ib),
'b') /= 0)
THEN
7366 WRITE(*,*)
'Debugging requested'
7368 it=index(text(ia:ib),
't')
7371 ieq=index(text(ia+it:ib),
'=')+it
7373 IF (index(text(ia+ieq:ib),
'SL0' ) /= 0)
ictest=2
7374 IF (index(text(ia+ieq:ib),
'SLE' ) /= 0)
ictest=3
7375 IF (index(text(ia+ieq:ib),
'BP' ) /= 0)
ictest=4
7376 IF (index(text(ia+ieq:ib),
'BRLF') /= 0)
ictest=5
7377 IF (index(text(ia+ieq:ib),
'BRLC') /= 0)
ictest=6
7380 IF(index(text(ia:ib),
's') /= 0)
isubit=1
7381 IF(index(text(ia:ib),
'f') /= 0)
iforce=1
7382 IF(index(text(ia:ib),
'c') /= 0)
icheck=1
7383 IF(index(text(ia:ib),
'C') /= 0)
icheck=2
7385 IF(i == iargc())
WRITE(*,*)
'--------------------- '
7406 CALL rltext(text,ia,ib,nab)
7407 nu=nufile(text(ia:ib))
7411 CALL peend(10,
'Aborted, no steering file')
7412 stop
'in FILETC: no steering file. .'
7425 WRITE(*,*)
'Listing of steering file: ',
filnam(1:
nfnam)
7426 WRITE(*,*)
'-------------------------'
7429 WRITE(*,*)
'Open error for steering file - stop'
7430 CALL peend(11,
'Aborted, open error for steering file')
7438 rows=6; cols=lenfileinfo
7439 CALL mpalloc(vecfileinfo,rows,cols,
'file info from steering')
7442 READ(10,102,iostat=ierrf) text
7444 CALL rltext(text,ia,ib,nab)
7446 IF(nline <= 50)
THEN
7447 WRITE(*,101) nline,text(1:nab)
7448 IF(nline == 50)
WRITE(*,*)
' ...'
7451 CALL rltext(text,ia,ib,nab)
7453 mat=matint(text(ia:ib),
'end',npat,ntext)
7457 WRITE(*,*)
' end-statement after',nline,
' text lines'
7462 keystx=
'fortranfiles'
7463 mat=matint(text(ia:ib),keystx,npat,ntext)
7464 IF(mat == ntext)
THEN
7471 mat=matint(text(ia:ib),keystx,npat,ntext)
7472 IF(mat == ntext)
THEN
7478 keystx=
'closeandreopen'
7479 mat=matint(text(ia:ib),keystx,npat,ntext)
7480 IF(mat == ntext)
THEN
7488 iopt=index(text(ia:ib),
' -- ')
7489 IF (iopt > 0) ie=iopt-1
7492 nu=nufile(text(ia:ie))
7494 IF (
nfiles == lenfileinfo)
THEN
7495 CALL mpalloc(temparray,rows,cols,
'temp file info from steering')
7496 temparray=vecfileinfo
7498 lenfileinfo=lenfileinfo*2
7500 CALL mpalloc(vecfileinfo,rows,newcols,
'file info from steering')
7501 vecfileinfo(:,1:cols)=temparray(:,1:cols)
7507 lenfilenames=lenfilenames+ie-ia+1
7508 vecfileinfo(1,
nfiles)=nline
7512 vecfileinfo(5,
nfiles)=iopt
7523 CALL mpalloc(
nfd,length,
'file line (in steering)')
7532 READ(10,102,iostat=ierrf) text
7535 IF (nline == vecfileinfo(1,i))
THEN
7536 nfd(i)=vecfileinfo(1,i)
7537 mfd(i)=vecfileinfo(2,i)
7538 ia=vecfileinfo(3,i)-1
7539 lfd(i)=vecfileinfo(4,i)-ia
7540 FORALL (k=1:
lfd(i))
tfd(ioff+k)=text(ia+k:ia+k)
7544 IF (vecfileinfo(5,i) > 0)
THEN
7545 CALL ratext(text(vecfileinfo(5,i)+4:vecfileinfo(6,i)),nums,dnum)
7546 IF (nums > 0)
ofd(i)=real(dnum(1),mps)
7556 CALL mpalloc(
ifd,length,
'integrated record numbers (=offset)')
7557 CALL mpalloc(
jfd,length,
'number of accepted records')
7558 CALL mpalloc(
kfd,rows,length,
'number of records in file, file order')
7563 CALL mpalloc(
sfd,rows,length,
'start, end of file name in TFD')
7567 WRITE(*,*)
'-------------------------'
7573 WRITE(*,*)
'Table of files:'
7574 WRITE(*,*)
'---------------'
7577 WRITE(8,*)
'Text and data files:'
7580 FORALL (k=1:
lfd(i)) fname(k:k)=
tfd(ioff+k)
7582 IF (
mprint > 1)
WRITE(*,103) i,bite(
mfd(i)),fname(1:
lfd(i))
7583 WRITE(8,103) i,bite(
mfd(i)),fname(1:
lfd(i))
7587 WRITE(*,*)
'---------------'
7601 IF(
mfd(i) == 3)
THEN
7625 IF(
mfd(i) == 1)
THEN
7646 WRITE(*,*)
'Opening of C-files not supported.'
7663 CALL peend(15,
'Aborted, open error(s) for binary files')
7664 stop
'FILETC: open error '
7667 CALL peend(14,
'Aborted, no binary files')
7668 stop
'FILETC: no binary files '
7671 WRITE(*,*)
nfilb,
' binary files opened'
7673 WRITE(*,*)
nfilb,
' binary files opened and closed'
7677103
FORMAT(i3,2x,a14,3x,a)
7740 INTEGER(mpi) :: ierrf
7741 INTEGER(mpi) :: ioff
7743 INTEGER(mpi) :: iosum
7747 INTEGER(mpi) :: nfiln
7748 INTEGER(mpi) :: nline
7749 INTEGER(mpi) :: nlinmx
7750 INTEGER(mpi) :: npat
7751 INTEGER(mpi) :: ntext
7752 INTEGER(mpi) :: matint
7756 CHARACTER (LEN=1024) :: text
7757 CHARACTER (LEN=1024) :: fname
7760 WRITE(*,*)
'Processing text files ...'
7773 IF(
mfd(i) /= 2) cycle
7774 FORALL (k=1:
lfd(i)) fname(k:k)=
tfd(ia+k)
7775 WRITE(*,*)
'File ',fname(1:
lfd(i))
7776 IF (
mprint > 1)
WRITE(*,*)
' '
7777 OPEN(10,file=fname(1:
lfd(i)),iostat=ios,form=
'FORMATTED')
7779 WRITE(*,*)
'Open error for file ',fname(1:
lfd(i))
7789 READ(10,102,iostat=ierrf) text
7793 WRITE(*,*)
' end-of-file after',nline,
' text lines'
7797 IF(nline <= nlinmx.AND.
mprint > 1)
THEN
7798 CALL rltext(text,ia,ib,nab)
7800 WRITE(*,101) nline,text(1:nab)
7801 IF(nline == nlinmx)
WRITE(*,*)
' ...'
7804 CALL rltext(text,ia,ib,nab)
7806 mat=matint(text(ia:ib),
'end',npat,ntext)
7810 WRITE(*,*)
' end-statement after',nline,
' text lines'
7816 IF(nfiln <=
nfiles.AND.nline ==
nfd(nfiln))
THEN
7831 CALL peend(16,
'Aborted, open error(s) for text files')
7832 stop
'FILETX: open error(s) in text files '
7835 WRITE(*,*)
'... end of text file processing.'
7840 WRITE(*,*)
lunkno,
' unknown keywords in steering files, ', &
7841 'or file non-existing,'
7842 WRITE(*,*)
' see above!'
7843 WRITE(*,*)
'------------> stop'
7845 CALL peend(13,
'Aborted, unknown keywords in steering file')
7855 ELSE IF(
matsto == 1)
THEN
7857 ELSE IF(
matsto == 2)
THEN
7860 ELSE IF(
metsol == 1)
THEN
7862 ELSE IF(
metsol == 2)
THEN
7864 ELSE IF(
metsol == 3)
THEN
7866 ELSE IF(
metsol == 4)
THEN
7868 ELSE IF(
metsol == 5)
THEN
7871 WRITE(*,*)
'MINRES forced with sparse matrix!'
7873 WRITE(*,*)
'MINRES forced with sparse matrix!'
7875 WRITE(*,*)
'MINRES forced with sparse matrix!'
7880 WRITE(*,*)
'MINRES forced with sparse matrix!'
7882 WRITE(*,*)
'MINRES forced with sparse matrix!'
7884 WRITE(*,*)
'MINRES forced with sparse matrix!'
7892 WRITE(*,*)
'Solution method and matrix-storage mode:'
7894 WRITE(*,*)
' METSOL = 1: matrix inversion'
7895 ELSE IF(
metsol == 2)
THEN
7896 WRITE(*,*)
' METSOL = 2: diagonalization'
7897 ELSE IF(
metsol == 3)
THEN
7898 WRITE(*,*)
' METSOL = 3: MINRES'
7899 ELSE IF(
metsol == 4)
THEN
7900 WRITE(*,*)
' METSOL = 4: MINRES-QLP'
7901 ELSE IF(
metsol == 5)
THEN
7902 WRITE(*,*)
' METSOL = 5: GMRES (-> MINRES)'
7906 WRITE(*,*)
' with',
mitera,
' iterations'
7909 WRITE(*,*)
' MATSTO = 1: symmetric matrix, ',
'(n*n+n)/2 elements'
7910 ELSE IF(
matsto == 2)
THEN
7911 WRITE(*,*)
' MATSTO = 2: sparse matrix'
7914 WRITE(*,*)
' and band matrix, width',
mbandw
7918 WRITE(*,*)
'Chi square cut equiv 3 st.dev applied ...'
7919 WRITE(*,*)
' in first iteration with factor',
chicut
7920 WRITE(*,*)
' in second iteration with factor',
chirem
7921 WRITE(*,*)
' (reduced by sqrt in next iterations)'
7925 WRITE(*,*)
' Down-weighting of outliers in',
lhuber,
' iterations'
7926 WRITE(*,*)
' Cut on downweight fraction',
dwcut
7929 WRITE(*,*)
'Iterations (solutions) with line search:'
7938 WRITE(*,*)
' All with Chi square cut scaling factor <= 1.'
7973 INTEGER(mpi) :: npat
7974 INTEGER(mpi) :: ntext
7975 INTEGER(mpi) :: nuprae
7978 CHARACTER (LEN=*),
INTENT(INOUT) :: fname
7983 IF(fname(1:5) ==
'rfio:') nuprae=1
7984 IF(fname(1:5) ==
'dcap:') nuprae=2
7985 IF(fname(1:5) ==
'root:') nuprae=3
7986 IF(nuprae == 0)
THEN
7987 INQUIRE(file=fname,iostat=ios,exist=ex)
7988 IF(ios /= 0)
nufile=-abs(ios)
7990 ELSE IF(nuprae == 1)
THEN
8003 nm=
matint(
'xt',fname(l1:ll),npat,ntext)
8006 nm=
matint(
'tx',fname(l1:ll),npat,ntext)
8028 INTEGER(mpi) :: iomp
8030 INTEGER(mpi) :: kkey
8031 INTEGER(mpi) :: label
8032 INTEGER(mpi) :: lkey
8034 INTEGER(mpi) :: miter
8036 INTEGER(mpi) :: nkey
8037 INTEGER(mpi) :: nkeys
8039 INTEGER(mpi) :: nmeth
8040 INTEGER(mpi) :: npat
8041 INTEGER(mpi) :: ntext
8042 INTEGER(mpi) :: nums
8043 INTEGER(mpi) :: matint
8045 CHARACTER (LEN=*),
INTENT(IN) :: text
8046 INTEGER(mpi),
INTENT(IN) :: nline
8048 parameter(nkeys=5,nmeth=6)
8049 CHARACTER (LEN=16) :: methxt(nmeth)
8050 CHARACTER (LEN=16) :: keylst(nkeys)
8051 CHARACTER (LEN=32) :: keywrd
8052 CHARACTER (LEN=32) :: keystx
8053 REAL(mpd) :: dnum(100)
8054 INTEGER(mpi) :: lpvs
8058 SUBROUTINE additem(length,list,label,value)
8060 INTEGER(mpi),
INTENT(IN OUT) :: length
8061 TYPE(listitem),
DIMENSION(:),
INTENT(IN OUT),
ALLOCATABLE :: list
8062 INTEGER(mpi),
INTENT(IN) :: label
8063 REAL(mpd),
INTENT(IN) :: value
8067 DATA keylst/
'unknown',
'parameter',
'constraint',
'measurement',
'method'/
8070 DATA methxt/
'diagonalization',
'inversion',
'fullMINRES',
'sparseMINRES', &
8071 'fullMINRES-QLP',
'sparseMINRES-QLP'/
8076 CALL rltext(text,ia,ib,nab)
8077 IF(nab == 0)
GOTO 10
8078 CALL ratext(text(1:nab),nums,dnum)
8080 IF(nums /= 0) nkey=0
8089 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
8090 IF(mat >= ntext-ntext/5)
GO TO 10
8096 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
8099 IF(mat >= (npat-npat/5))
THEN
8102 IF(nums > 0) mprint=nint(dnum(1),mpi)
8107 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
8108 IF(mat >= (npat-npat/5))
THEN
8112 IF(nums > 0) mdebug=nint(dnum(1),mpi)
8113 IF(nums > 1) mdebg2=nint(dnum(2),mpi)
8118 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
8119 IF(mat >= (npat-npat/5))
THEN
8121 IF(nums > 0 .AND. dnum(1) > 0.5) mreqenf=nint(dnum(1),mpi)
8122 IF(nums > 1 .AND. dnum(2) > 0.5) mreqena=nint(dnum(2),mpi)
8123 IF(nums > 2 .AND. dnum(3) > 0.5) iteren=nint(dnum(1)*dnum(3),mpi)
8127 keystx=
'printrecord'
8128 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
8129 IF(mat >= (npat-npat/5))
THEN
8130 IF(nums > 0) nrecpr=nint(dnum(1),mpi)
8131 IF(nums > 1) nrecp2=nint(dnum(2),mpi)
8136 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
8137 IF(mat >= (npat-npat/5))
THEN
8138 IF (nums > 0.AND.dnum(1) > 0.) mxrec=nint(dnum(1),mpi)
8143 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
8144 IF(mat >= (npat-npat/5))
THEN
8145 IF (nums > 0.AND.dnum(1) >= 0.) ncache=nint(dnum(1),mpi)
8146 IF (nums == 2.AND.dnum(2) > 0..AND.dnum(2) <= 1.0) &
8147 fcache(1)=real(dnum(2),mps)
8150 fcache(k)=real(dnum(k+1),mps)
8157 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
8158 IF(mat >= (npat-npat/5))
THEN
8163 chicut=real(dnum(1),mps)
8164 IF(chicut < 1.0) chicut=-1.0
8168 chirem=real(dnum(2),mps)
8169 IF(chirem < 1.0) chirem=1.0
8170 IF(chicut >= 1.0) chirem=min(chirem,chicut)
8178 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
8179 IF(mat >= (npat-npat/5))
THEN
8180 IF(nums > 0) chhuge=real(dnum(1),mps)
8181 IF(chhuge < 1.0) chhuge=1.0
8187 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
8188 IF(mat >= (npat-npat/5))
THEN
8189 IF(nums > 0) lsearch=nint(dnum(1),mpi)
8194 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
8195 IF(mat >= (npat-npat/5))
THEN
8196 IF(nums > 0) lfitnp=nint(dnum(1),mpi)
8197 IF(nums > 1) lfitbb=nint(dnum(2),mpi)
8201 keystx=
'regularization'
8202 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
8203 IF(mat >= (npat-npat/5))
THEN
8205 regula=real(dnum(1),mps)
8206 IF(nums >= 2) regpre=real(dnum(2),mps)
8210 keystx=
'regularisation'
8211 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
8212 IF(mat >= (npat-npat/5))
THEN
8214 regula=real(dnum(1),mps)
8215 IF(nums >= 2) regpre=real(dnum(2),mps)
8220 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
8221 IF(mat >= (npat-npat/5))
THEN
8222 regpre=real(dnum(1),mps)
8227 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
8228 IF(mat >= (npat-npat/5))
THEN
8229 matrit=nint(dnum(1),mpi)
8234 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
8235 IF(mat >= (npat-npat/5))
THEN
8237 IF (nums > 0.AND.dnum(1) > 0.) matmon=nint(dnum(1),mpi)
8242 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
8244 IF(mat >= (npat-npat/5))
THEN
8245 IF(nums > 0) mbandw=nint(dnum(1),mpi)
8246 IF(mbandw < 0) mbandw=-1
8247 IF(nums > 1) lprecm=nint(dnum(2),mpi)
8271 keystx=
'outlierdownweighting'
8272 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
8275 IF(mat >= (npat-npat/5))
THEN
8276 lhuber=nint(dnum(1),mpi)
8277 IF(lhuber > 0.AND.lhuber <= 2) lhuber=2
8281 keystx=
'dwfractioncut'
8282 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
8283 IF(mat >= (npat-npat/5))
THEN
8284 dwcut=real(dnum(1),mps)
8285 IF(dwcut > 0.5) dwcut=0.5
8290 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
8291 IF(mat >= (npat-npat/5))
THEN
8292 prange=abs(real(dnum(1),mps))
8297 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
8298 IF(mat >= (npat-npat/5))
THEN
8304 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
8305 IF(mat >= (npat-npat/5))
THEN
8310 keystx=
'memorydebug'
8311 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
8312 IF(mat >= (npat-npat/5))
THEN
8314 IF (nums > 0.AND.dnum(1) > 0.0) memdbg=nint(dnum(1),mpi)
8319 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
8320 IF(mat >= (npat-npat/5))
THEN
8325 keystx=
'printcounts'
8326 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
8327 IF(mat >= (npat-npat/5))
THEN
8329 IF (nums > 0.AND.dnum(1) > 0.0) ipcntr=nint(dnum(1),mpi)
8333 keystx=
'weightedcons'
8334 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
8335 IF(mat >= (npat-npat/5))
THEN
8337 IF (nums > 0) iwcons=nint(dnum(1),mpi)
8341 keystx=
'skipemptycons'
8342 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
8343 IF(mat >= (npat-npat/5))
THEN
8348 keystx=
'withelimination'
8349 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
8350 IF(mat >= (npat-npat/5))
THEN
8355 keystx=
'withmultipliers'
8356 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
8357 IF(mat >= (npat-npat/5))
THEN
8363 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
8364 IF(mat >= (npat-npat/5))
THEN
8366 IF (nums > 0) icheck=nint(dnum(1),mpi)
8370 keystx=
'monitorresiduals'
8371 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
8372 IF(mat >= (npat-npat/5))
THEN
8374 IF (nums > 0) imonit=nint(dnum(1),mpi)
8375 IF (nums > 1) measbins=max(measbins,nint(dnum(2),mpi))
8379 keystx=
'monitorpulls'
8380 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
8381 IF(mat >= (npat-npat/5))
THEN
8384 IF (nums > 0) imonit=nint(dnum(1),mpi)
8385 IF (nums > 1) measbins=max(measbins,nint(dnum(2),mpi))
8389 keystx=
'scaleerrors'
8390 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
8391 IF(mat >= (npat-npat/5))
THEN
8393 IF (nums > 0) dscerr(1:2)=dnum(1)
8394 IF (nums > 1) dscerr(2)=dnum(2)
8398 keystx=
'iterateentries'
8399 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
8400 IF(mat >= (npat-npat/5))
THEN
8402 IF (nums > 0) iteren=nint(dnum(1),mpi)
8407 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
8408 IF(mat >= (npat-npat/5))
THEN
8418 WRITE(*,*)
'WARNING: multithreading not available'
8424 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
8425 IF(mat >= (npat-npat/5))
THEN
8440 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
8441 IF(mat >= (npat-npat/5).AND.mnrsel < 100)
THEN
8442 nl=min(nums,100-mnrsel)
8444 lbmnrs(mnrsel+k)=nint(dnum(k),mpi)
8450 keystx=
'pairentries'
8451 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
8452 IF(mat >= (npat-npat/5))
THEN
8455 IF (nums > 0.AND.dnum(1) > 0.0)
THEN
8456 mreqpe=nint(dnum(1),mpi)
8457 IF (nums >= 2.AND.dnum(2) >= dnum(1)) mhispe=nint(dnum(2),mpi)
8458 IF (nums >= 3.AND.dnum(3) >= dnum(1)) msngpe=nint(dnum(3),mpi)
8464 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
8465 IF(mat >= (npat-npat/5))
THEN
8466 wolfc1=real(dnum(1),mps)
8467 wolfc2=real(dnum(2),mps)
8474 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
8475 IF(mat >= (npat-npat/5))
THEN
8477 IF (dnum(1) < 1.0e-10_mpd.OR.dnum(1) > 1.0e-04_mpd)
THEN
8478 WRITE(*,*)
'ERROR: need 1.0D-10 <= MRESTL ', &
8479 '<= 1.0D-04, but get ', dnum(1)
8488 keystx=
'mrestranscond'
8489 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
8490 IF(mat >= (npat-npat/5))
THEN
8498 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
8499 IF(mat >= (npat-npat/5))
THEN
8501 mrmode = int(dnum(1),mpi)
8506 keystx=
'nofeasiblestart'
8507 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
8508 IF(mat >= (npat-npat/5))
THEN
8514 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
8515 IF(mat >= ntext-ntext/10)
THEN
8521 keystx=
'fortranfiles'
8522 mat=matint(text(ia:ib),keystx,npat,ntext)
8523 IF(mat >= ntext-ntext/10)
RETURN
8526 mat=matint(text(ia:ib),keystx,npat,ntext)
8527 IF(mat >= ntext-ntext/10)
RETURN
8529 keystx=
'closeandreopen'
8530 mat=matint(text(ia:ib),keystx,npat,ntext)
8531 IF(mat >= ntext-ntext/10)
RETURN
8535 IF(nums /= 0) nkey=0
8538 WRITE(*,*)
'**************************************************'
8540 WRITE(*,*)
'Unknown keyword(s): ',text(1:min(nab,50))
8542 WRITE(*,*)
'**************************************************'
8568 lpvs=nint(dnum(1),mpi)
8570 CALL additem(lenparameters,listparameters,lpvs,dnum(2))
8571 CALL additem(lenpresigmas,listpresigmas,lpvs,dnum(3))
8573 WRITE(*,*)
'Line',nline,
' error, label=',lpvs
8575 ELSE IF(nums /= 0)
THEN
8577 WRITE(*,*)
'Wrong text in line',nline
8578 WRITE(*,*)
'Status: new parameter'
8579 WRITE(*,*)
'> ',text(1:nab)
8581 ELSE IF(lkey == 3)
THEN
8583 IF(nums >= 1.AND.nums <= 2)
THEN
8585 CALL additem(lenconstraints,listconstraints,lpvs,dnum(1))
8587 IF(iwcons > 0) lpvs=-2
8589 IF(nums == 2) plvs=dnum(2)
8590 CALL additem(lenconstraints,listconstraints,lpvs,plvs)
8593 WRITE(*,*)
'Wrong text in line',nline
8594 WRITE(*,*)
'Status: new keyword constraint'
8595 WRITE(*,*)
'> ',text(1:nab)
8597 ELSE IF(lkey == 4)
THEN
8599 nummeasurements=nummeasurements+1
8601 CALL additem(lenmeasurements,listmeasurements,lpvs,dnum(1))
8603 CALL additem(lenmeasurements,listmeasurements,lpvs,dnum(2))
8606 WRITE(*,*)
'Wrong text in line',nline
8607 WRITE(*,*)
'Status: new keyword measurement'
8608 WRITE(*,*)
'> ',text(1:nab)
8611 ELSE IF(lkey == 5)
THEN
8613 IF(nums >= 1) miter=nint(dnum(1),mpi)
8614 IF(miter >= 1) mitera=miter
8615 dflim=real(dnum(2),mps)
8619 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
8620 IF(mat >= ntext-ntext/5)
THEN
8624 ELSE IF(i == 2)
THEN
8627 ELSE IF(i == 3)
THEN
8630 ELSE IF(i == 4)
THEN
8633 ELSE IF(i == 5)
THEN
8636 ELSE IF(i == 6)
THEN
8643 ELSE IF(nkey == 0)
THEN
8646 lpvs=nint(dnum(1),mpi)
8648 CALL additem(lenparameters,listparameters,lpvs,dnum(2))
8649 CALL additem(lenpresigmas,listpresigmas,lpvs,dnum(3))
8651 WRITE(*,*)
'Line',nline,
' error, label=',lpvs
8653 ELSE IF(nums > 1.AND.nums < 3)
THEN
8655 WRITE(*,*)
'Wrong text in line',nline
8656 WRITE(*,*)
'Status continuation parameter'
8657 WRITE(*,*)
'> ',text(1:nab)
8660 ELSE IF(lkey == 3)
THEN
8663 label=nint(dnum(i),mpi)
8664 IF(label <= 0) ier=1
8666 IF(mod(nums,2) /= 0) ier=1
8669 lpvs=nint(dnum(i),mpi)
8671 CALL additem(lenconstraints,listconstraints,lpvs,plvs)
8675 WRITE(*,*)
'Wrong text in line',nline
8676 WRITE(*,*)
'Status continuation constraint'
8677 WRITE(*,*)
'> ',text(1:nab)
8680 ELSE IF(lkey == 4)
THEN
8684 label=nint(dnum(i),mpi)
8685 IF(label <= 0) ier=1
8687 IF(mod(nums,2) /= 0) ier=1
8691 lpvs=nint(dnum(i),mpi)
8693 CALL additem(lenmeasurements,listmeasurements,lpvs,plvs)
8697 WRITE(*,*)
'Wrong text in line',nline
8698 WRITE(*,*)
'Status continuation measurement'
8699 WRITE(*,*)
'> ',text(1:nab)
8717 INTEGER(mpi),
INTENT(IN OUT) :: length
8718 TYPE(
listitem),
DIMENSION(:),
INTENT(IN OUT),
ALLOCATABLE :: list
8719 INTEGER(mpi),
INTENT(IN) :: label
8720 REAL(mpd),
INTENT(IN) :: value
8722 INTEGER(mpl) :: newSize
8723 INTEGER(mpl) :: oldSize
8724 TYPE(
listitem),
DIMENSION(:),
ALLOCATABLE :: tempList
8726 IF (label > 0.AND.
value == 0.)
RETURN
8727 IF (length == 0 )
THEN
8729 CALL mpalloc(list,newsize,
' list ')
8731 oldsize=
size(list,kind=
mpl)
8732 IF (length >= oldsize)
THEN
8733 newsize = oldsize + oldsize/5 + 100
8734 CALL mpalloc(templist,oldsize,
' temp. list ')
8737 CALL mpalloc(list,newsize,
' list ')
8738 list(1:oldsize)=templist(1:oldsize)
8743 list(length)%label=label
8744 list(length)%value=
value
8758 CHARACTER (LEN=*),
INTENT(IN) :: text
8759 CHARACTER (LEN=16) :: textc
8768 textl(ka:kb)=text(1:l)
8772 textc=text(1:l)//
'-end'
8780 textl(ka:kb)=textc(1:l)
8807 INTEGER(mpi),
INTENT(IN) :: lun
8808 CHARACTER (LEN=*),
INTENT(IN) :: fname
8809 CHARACTER (LEN=33) :: nafile
8810 CHARACTER (LEN=33) :: nbfile
8816 CALL peend(17,
'Aborted, file name too long')
8817 stop
'File name too long '
8822 INQUIRE(file=nafile(1:l),exist=ex)
8824 INQUIRE(file=nafile(1:l+1),exist=ex)
8826 CALL system(
'rm '//nafile)
8830 CALL system(
'mv '//nafile//nbfile)
8832 OPEN(unit=lun,file=fname)
8843 REAL,
DIMENSION(2) :: ta
8864 nsecd1=int(delta,
mpi)
8866 minut1=nsecd1/60-60*nhour1
8867 secnd1=delta-60*(minut1+60*nhour1)
8869 nsecd2=int(delta,
mpi)
8871 minut2=nsecd2/60-60*nhour2
8872 secnd2=delta-60*(minut2+60*nhour2)
8873 WRITE(*,101) nhour1,minut1,secnd1, nhour2,minut2,secnd2
8878101
FORMAT(i4,
' h',i3,
' min',f5.1,
' sec total',18x,
'elapsed', &
8879 i4,
' h',i3,
' min',f5.1,
' sec')
8893 INTEGER(mpi),
INTENT(IN) :: icode
8894 CHARACTER (LEN=*),
INTENT(IN) :: cmessage
8896 CALL mvopen(9,
'millepede.end')
8897 WRITE(9,101) icode, cmessage
8898101
FORMAT(1x,i4,3x,a)
8913 INTEGER(mpi),
INTENT(IN) :: kfile
8914 INTEGER(mpi),
INTENT(IN) :: ithr
8915 INTEGER(mpi),
INTENT(OUT) :: ierr
8917 INTEGER(mpi),
DIMENSION(13) :: ibuff
8918 INTEGER(mpi) :: ioff
8923 INTEGER(mpi) :: moddate
8924 CHARACTER (LEN=1024) :: fname
8933 FORALL (k=1:lfn) fname(k:k)=
tfd(ioff+k)
8937 IF(kfile <=
nfilf)
THEN
8940 OPEN(lun,file=fname(1:lfn),iostat=ios, form=
'UNFORMATTED')
8941 print *,
' lun ', lun, ios
8945 CALL openc(fname(1:lfn),lun,ios)
8947 WRITE(*,*)
'Opening of C-files not supported.'
8954 WRITE(*,*)
'Open error for file ',fname(1:lfn), ios
8955 IF (moddate /= 0)
THEN
8956 CALL peend(15,
'Aborted, open error(s) for binary files')
8957 stop
'PEREAD: open error '
8962 CALL stat(fname(1:lfn),ibuff,ios)
8966 WRITE(*,*)
'STAT error for file ',fname(1:lfn), ios
8970 IF (moddate /= 0)
THEN
8971 IF (ibuff(10) /= moddate)
THEN
8972 CALL peend(19,
'Aborted, binary file(s) modified')
8973 stop
'PEREAD: file modified '
8976 yfd(kfile)=ibuff(10)
8991 INTEGER(mpi),
INTENT(IN) :: kfile
8992 INTEGER(mpi),
INTENT(IN) :: ithr
8998 IF(kfile <=
nfilf)
THEN
9017 INTEGER(mpi),
INTENT(IN) :: kfile
9022 IF (kfile <=
nfilf)
THEN
9070 REAL(mpd),
INTENT(OUT) ::asum
subroutine ptlopt(nf, m, slopes, steps)
Get details.
subroutine ptline(n, x, f, g, s, step, info)
Perform linesearch.
subroutine ptldef(gtole, stmax, minfe, maxfe)
Initialize line search.
subroutine ptlprt(lunp)
Print line search data.
subroutine clbits(in, jreqpe, jhispe, jsngpe, jcmprs, jextnd, idimb, iencdb, ispc)
Calculate bit (field) array size, encoding.
subroutine gpbmap(npair)
Get pairs (statistic) from map.
subroutine ckbits(ndims)
Check sparsity of matrix.
subroutine spbits(nsparr, nsparc, ncmprs)
Create sparsity information.
subroutine clbmap(in)
Clear (additional) bit map.
subroutine inbmap(im, jm)
Fill bit map.
subroutine inbits(im, jm, inc)
Fill bit fields (counters).
subroutine ndbits(ndims, ncmprs, nsparr, ihst)
Analyze bit fields.
subroutine gmpdef(ig, ityp, text)
subroutine hmpdef(ih, xa, xb, text)
subroutine sqmibb2(v, b, n, nbdr, nbnd, inv, nrank, vbnd, vbdr, aux, vbk, vzru, scdiag, scflag)
Band bordered matrix.
subroutine precon(p, n, c, cu, a, s, nrkd)
Constrained preconditioner, decomposition.
subroutine devsol(n, diag, u, b, x, work)
Solution by diagonalization.
subroutine sqminl(v, b, n, nrank, diag, next)
Matrix inversion for LARGE matrices.
subroutine dbsvxl(v, a, b, n)
Product LARGE symmetric matrix, vector.
subroutine devrot(n, diag, u, v, work, iwork)
Diagonalization.
subroutine equslv(n, m, c, india, x)
Solution of equilibrium systems (after decomposition).
subroutine equdec(n, m, ls, c, india, nrkd, nrkd2)
Decomposition of equilibrium systems.
subroutine sort1k(a, n)
Quick sort 1.
subroutine sqminv(v, b, n, nrank, diag, next)
Matrix inversion and solution.
subroutine sqmibb(v, b, n, nbdr, nbnd, inv, nrank, vbnd, vbdr, aux, vbk, vzru, scdiag, scflag)
Bordered band matrix.
subroutine presol(p, n, cu, a, s, x, y)
Constrained preconditioner, solution.
subroutine devinv(n, diag, u, v)
Inversion by diagonalization. Get inverse matrix V from DIAG and U.
subroutine sort2k(a, n)
Quick sort 2.
subroutine devsig(n, diag, u, b, coef)
Calculate significances.
subroutine dbsvx(v, a, b, n)
Product symmetric matrix, vector.
subroutine dbavat(v, a, w, n, ms)
A V AT product (similarity).
subroutine sort2i(a, n)
Quick sort 2 with index.
subroutine qldecb(a, nb, b)
QL decomposition (for disjoint block matrix).
subroutine qlmlq(x, m, t)
Multiply left by Q(t).
subroutine qlbsub(d, y)
Backward substitution.
subroutine qlgete(emin, emax)
Get eigenvalues.
subroutine qlssq(aprod, a, t)
Similarity transformation by Q(t).
subroutine qlpssq(aprod, b, m, t)
Partial similarity transformation by Q(t).
subroutine qlini(n, m)
Initialize QL decomposition.
subroutine mptest
Generate test files.
subroutine mptst2(imodel)
Generate test files.
integer(mpi) function matint(pat, text, npat, ntext)
Approximate string matching.
subroutine ratext(text, nums, dnum)
Translate text.
subroutine rltext(text, ia, ib, nab)
Analyse text range.
MINRES solves symmetric systems Ax = b or min ||Ax - b||_2, where the matrix A may be indefinite and/...
subroutine, public minres(n, aprod, msolve, b, shift, checka, precon, x, itnlim, nout, rtol, istop, itn, anorm, acond, rnorm, arnorm, ynorm)
Solution of linear equation system.
MINRESQLP solves symmetric systems Ax = b or min ||Ax - b||_2, where the matrix A may be indefinite a...
subroutine, public minresqlp(n, aprod, b, shift, msolve, disable, nout, itnlim, rtol, maxxnorm, trancond, acondlim, x, istop, itn, rnorm, arnorm, xnorm, anorm, acond)
Solution of linear equation system or least squares problem.
(De)Allocate vectors and arrays.
integer(mpl) maxwordsalloc
peak dynamic memory allocation (words)
integer(mpi) printflagalloc
print flag for dynamic allocations
Parameters, variables, dynamic arrays.
real(mpd), dimension(:), allocatable workspaceeigenvectors
workspace eigen vectors
real(mpd), dimension(:), allocatable globalparameter
global parameters (start values + sum(x_i))
type(listitem), dimension(:), allocatable listparameters
list of parameters from steering file
integer(mpi) lunmon
unit for monitoring output file
real(mpd), dimension(:), allocatable vecconsresiduals
residuals of constraints
integer(mpi) iskpec
flag for skipping empty constraints (no variable parameters)
integer(mpi) mnrsel
number of MINRES error labels in LBMNRS (calc err, corr with SOLGLO)
integer(mpi) nrecer
record with error (rank deficit or Not-a-Number) for printout
real(mps) actfun
actual function change
integer(mpi), dimension(:), allocatable globalindexusage
indices of global par in record
real(mps) regpre
default presigma
integer(mpi) mnrsit
total number of MINRES internal iterations
integer(mpi) metsol
solution method (1: inversion, 2: diagonalization, 3: MINRES-QLP)
integer(mpi) nagbn
max number of global paramters per record
character(len=74) textl
name of current MP 'module' (step)
integer(mpi) nloopn
number of data reading, fitting loops
integer(mpi) mreqpe
min number of pair entries
integer(mpi) memdbg
debug flag for memory management
integer(mpi), dimension(100) lbmnrs
MINRES error labels.
real(mpd) mrtcnd
transition (QR -> QLP) (matrix) condition for MINRES-QLP
real(mpd), dimension(:), allocatable vbk
local fit 'matrix for border solution'
real(mps) prange
range (-PRANGE..PRANGE) for histograms of pulls, norm. residuals
integer(mpi) matsto
(global) matrix storage mode (1: full, 2: sparse)
integer(mpi), dimension(:,:), allocatable matconssort
keys and index for sorting
integer(mpi), dimension(:,:), allocatable readbufferinfo
buffer management (per thread)
integer(mpi) nhistp
flag for histogram printout
real(mpd), dimension(:), allocatable globalparcopy
copy of global parameters
real(mpd), dimension(2) dscerr
scaling factors for errors of 'global' and 'local' measurement
real(mps) chhuge
cut in terms of 3-sigma for unreasonable data, all iterations
integer(mpi), dimension(:), allocatable sparsematrixcolumns
(compressed) list of columns for sparse matrix
integer(mpl), dimension(:,:), allocatable sparsematrixoffsets
row offsets for column list, sparse matrix elements
integer(mpi) iteren
entries cut is iterated for parameters with less entries (if > mreqenf)
integer(mpi) lunkno
flag for unkown keywords
integer(mpi), dimension(:), allocatable scflag
local fit workspace (I)
real(mpd), parameter measbinsize
bins size for monitoring
integer(mpi) nrecp2
record number with printout
integer(mpi) mdebug
debug flag (number of records to print)
real(mpd), dimension(:), allocatable matconsproduct
product matrix of constraints
integer(mpi), dimension(:), allocatable yfd
binary file: modification date
integer(mpi) nrecpr
record number with printout
real(mps) value1
largest residual
real(mpd), dimension(:), allocatable localcorrections
local fit corrections (to residuals)
real(mps) chirem
cut in terms of 3-sigma cut, other iterations, approaching 1.
real(mpd), dimension(:), allocatable localglobalmatrix
matrix correlating local and global par
integer(mpi) mhispe
upper bound for pair entry histogrammimg
integer(mpi) nfgb
number of fit parameters
integer(mpi) nrec3
(1.) record number with error
integer(mpi), dimension(:,:), allocatable kfd
(1,.)= number of records in file, (2,..)= file order
integer(mpi) icheck
flag for checking input only (no solution determined)
integer(mpi), parameter nexp20
integer(mpi), dimension(:), allocatable jfd
file: number of accepted records
integer(mpi) matmon
record interval for monitoring of (sparse) matrix construction
integer(mpi) nbndx
max band width for local fit
type(listitem), dimension(:), allocatable listconstraints
list of constraints from steering file
real(mpd), dimension(:), allocatable globalmatd
global matrix 'A' (double, full or sparse)
real(mpr8), dimension(:), allocatable readbufferdatad
double data
type(listitem), dimension(:), allocatable listmeasurements
list of (external) measurements from steering file
integer(mpi) lsinfo
line search: returned information
integer(mpi) nregul
regularization flag
integer(mpi), dimension(-7:0) globalparheader
global parameters (mapping) header
integer(mpi) nfilw
number of weighted binary files
integer(mpi), dimension(:), allocatable paircounter
number of paired parameters (in equations)
integer(mpi) nummeasurements
number of (external) measurements from steering file
integer(mpi), dimension(0:3) nrejec
rejected events
integer(mpi) ndimbuf
default read buffer size (I/F words, half record length)
real(mpd) fvalue
function value (chi2 sum) solution
real(mpd), dimension(:), allocatable globalcorrections
correction x_i (from A*x_i=b_i in iteration i)
real(mps), dimension(:), allocatable cfd
file: chi2 sum
real(mps) regula
regularization parameter, add regula * norm(global par.) to objective function
integer(mpi) nspc
number of precision for sparse global matrix (1=D, 2=D+F)
integer(mpi) nfilc
number of C binary files
integer(mpi) nagb
number of all parameters (global par. + Lagrange mult.)
integer(mpi) nmiss1
rank deficit for constraints
integer(mpi), dimension(:), allocatable globalparhashtable
global parameters hash table
integer(mpi) nalow
(sum of) global parameters with too few accepted entries
integer(mpi) iscerr
flag for scaling of errors
integer(mpi), dimension(:), allocatable globalcounter
global counter (entries in 'x')
real(mpd) sumndf
weighted sum(ndf)
integer(mpi), dimension(2) nbndr
number of records with bordered band matrix for local fit (upper/left, lower/right)
integer(mpi) iterat
iterations in solution
real(mpd) flines
function value line search
integer(mpi), dimension(:), allocatable meashists
measurement histograms (100 bins per thread)
integer(mpi) mthrd
number of (OpenMP) threads
integer(mpi) mbandw
band width of preconditioner matrix
real(mps) dwcut
down-weight fraction cut
real(mps), dimension(:), allocatable globalmatf
global matrix 'A' (float part for compressed sparse)
real(mps), dimension(0:8) times
cpu time counters
integer(mpi) minrecordsinblock
min. records in block
integer(mpi) naeqn
max number of equations (measurements) per record
integer(mpi) nfilb
number of binary files
real(mpd), dimension(:), allocatable vzru
local fit 'border solution'
real(mpd), dimension(:), allocatable globalparpreweight
weight from pre-sigma
integer(mpi) ictest
test mode '-t'
real(mpd), dimension(:), allocatable vbdr
local fit border part of 'A'
integer(mpi) mdebg2
number of measurements for record debug printout
real(mps) deltim
cpu time difference
integer(mpi) igcorr
flag for output of global correlations for inversion, =0: none
real(mpd), dimension(:), allocatable vecconssolution
solution for constraint elimination
integer(mpi), dimension(2) nprecond
number of constraints, matrix size for preconditioner
integer(mpi) nfiles
number of files
integer(mpi) ipcntr
flag for output of global parameter counts (entries), =0: none, =1: local fits, >1: binary files
integer(mpi) keepopen
flag for keeping binary files open
real(mpd), dimension(:), allocatable workspacediagonalization
workspace diag.
integer(mpi) nrec
number of records read
real(mps), dimension(:), allocatable wfd
binary file: weight
integer(mpi) accuratenexp
sum / 2**20
real(mpd), dimension(:), allocatable matprecond
preconditioner (band) matrix
integer(mpi) nrec1
record number with largest residual
integer(mpi) ntgb
total number of global parameters
integer(mpi) mszprd
(integrated block) matrix size for (constraint) product matrix
real(mps) angras
angle between gradient and search direction
integer(mpi) mthrdr
number of threads for reading binary files
integer(mpi) numreadbuffer
number of buffers (records) in (read) block
integer(mpi) imonmd
monitoring mode: 0:residuals (normalized to average error), 1:pulls
character(len=1024) filnam
name of steering file
integer(mpi) lunlog
unit for logfile
integer(mpi) ncblck
number of (disjoint) constraint blocks
real(mps), dimension(3) fcache
read cache, average fill level; write cache; dynamic size
real(mps) wolfc2
C_2 of strong Wolfe condition.
real(mpd) accuratedsum
fractional part of sum
integer(mpi) maxrecordsinblock
max. records in block
real(mpd) mrestl
tolerance criterion for MINRES-QLP
real(mpd), dimension(:), allocatable globalparpresigma
pre-sigma for global parameters
integer(mpi) icelim
flag for using elimination (instead of multipliers) for constraints
integer(mpi) mitera
number of iterations
integer(mpi) nbdrx
max border size for local fit
integer(mpi), dimension(:,:), allocatable globalparlabelindex
global parameters label, total -> var. index
real(mpd), dimension(:), allocatable scdiag
local fit workspace (D)
integer(mpi), dimension(:), allocatable readbufferdatai
integer data
integer(mpi) mextnd
flag for extended storage (both 'halves' of sym. mat. for improved access patterns)
integer(mpi), dimension(:,:), allocatable sfd
offset (1,..), length (2,..) of binary file name in tfd
integer(mpi) ndfsum
sum(ndf)
integer(mpi) mcmprs
compression flag for sparsity (column indices)
integer(mpi) lenconstraints
length of list of constraints from steering file
integer(mpi) lenparameters
list items from steering file
integer(mpi) lprecm
additional flag for preconditioner (band) matrix (>0: preserve rank by skyline matrix)
integer(mpi) ndefec
rank deficit for global matrix (from inversion)
integer(mpi), dimension(:), allocatable ifd
file: integrated record numbers (=offset)
integer(mpi) nofeas
flag for skipping making parameters feasible
integer(mpi) nfnam
length of sterring file name
real rstart
cpu start time for solution iterations
integer(mpi), dimension(:), allocatable writebufferindices
write buffer for indices
integer(mpi) iforce
switch to SUBITO for (global) rank defects if zero
real(mpd), dimension(:), allocatable workspacelinesearch
workspace line search
integer(mpi), dimension(:), allocatable globalparvartototal
global parameters variable -> total index
real(mpd), dimension(:), allocatable clmat
local fit matrix 'A' (in A*x=b)
integer(mpi), dimension(:), allocatable lfd
length of file name
character, dimension(:), allocatable tfd
file names (concatenation)
integer(mpi) ncgbe
number of empty constraints (no variable parameters)
integer(mpi) mprint
print flag (0: minimal, 1: normal, >1: more)
integer(mpi), dimension(:), allocatable vecconsstart
start of constraint in listConstraints (unsorted input)
integer(mpi) nummeas
number of measurement groups for monitoring
integer(mpi) lvllog
log level
integer(mpi) nalcn
max number of local paramters per record
integer(mpi) mreqenf
required number of entries (for variable global parameter from binary Files)
real(mps) value2
largest chi^2/Ndf
integer(mpi) icalcm
calculation mode (for XLOOPN) , >0: calculate matrix
real(mps), dimension(:), allocatable ofd
file: option
integer(mpi) ifile
current file (index)
real(mps) delfun
expected function change
integer(mpi) iitera
MINRES iterations.
integer(mpi) lenmeasurements
length of list of (external) measurements from steering file
real(mps) wolfc1
C_1 of strong Wolfe condition.
real(mpd), dimension(:), allocatable aux
local fit 'solutions for border rows'
integer(mpi) nrecal
number of records
integer(mpi) skippedrecords
number of skipped records (buffer too small)
integer(mpi), dimension(:), allocatable backindexusage
list of global par in record
integer(mpi), dimension(:), allocatable globalparcounts
global parameters counts (from binary files)
integer(mpi), dimension(:), allocatable ibandh
local fit 'band width histogram' (band size autodetection)
integer(mpi) isubit
subito flag '-s'
integer(mpi), dimension(:), allocatable indprecond
preconditioner pointer array
real(mps) dflim
convergence limit
integer(mpi) ncache
buffer size for caching (default 100MB per thread)
integer(mpi) mxrec
max number of records
integer(mpi) nrecd
number of records read containing doubles
integer(mpi) lfitnp
local fit: number of iteration to calculate pulls
integer(mpi) lcalcm
last calclation mode
real(mpd), dimension(:), allocatable globalvector
global vector 'x' (in A*x=b)
real(mpd), dimension(:), allocatable writebufferupdates
write buffer for update matrices
real(mpd), dimension(:), allocatable workspaced
(general) workspace (D)
integer(mpi) measbins
number of bins per measurement for monitoring
integer(mpi) accuratensum
sum mod 2**20
integer(mpi), dimension(:), allocatable nfd
index (line) in (steering) file
integer(mpi) numblocks
number of (read) blocks
integer(mpi) ncgb
number of constraints
integer(mpi), dimension(:,:), allocatable matconsblocks
start of constraint blocks, parameter range
real(mpd), dimension(:), allocatable workspaceeigenvalues
workspace eigen values
integer(mpi) lhuber
Huber down-weighting flag.
integer(mpi) nvgb
number of variable global parameters
integer(mpi) nfilf
number of Fortran binary files
integer(mpi) mszcon
(integrated block) matrix size for constraint matrix
integer(mpi), dimension(:), allocatable measindex
mapping of 1. global label to measurement index
integer(mpi) istopa
MINRES istop (convergence)
integer(mpi), dimension(:), allocatable mfd
file mode: cbinary =1, text =2, fbinary=3
integer(mpi) nrec2
record number with largest chi^2/Ndf
integer(mpi) nencdb
encoding info (number bits for column counter)
real(mpd), dimension(:), allocatable blvec
local fit vector 'b' (in A*x=b), replaced by 'x'
logical newite
flag for new iteration
real(mpd), dimension(:), allocatable measres
average measurement error
real(mpd), dimension(:), allocatable vecxav
vector x for AVPROD (A*x=b)
real(mpd), dimension(:), allocatable globalparstart
start value for global parameters
integer(mpi), dimension(-6:6) writebufferheader
write buffer header (-6..-1: updates, 1..6: indices)
integer(mpi) lenpresigmas
length of list of pre-sigmas from steering file
integer(mpi) npresg
number of pre-sigmas
integer(mpi), dimension(:), allocatable appearancecounter
appearance statistics for global par (first/last file,record)
integer(mpi), dimension(:), allocatable xfd
file: max. record size
integer(mpi) mreqena
required number of entries (for variable global parameter from Accepted local fits)
integer(mpi) sumrecords
sum of records
real(mps), dimension(:,:), allocatable writebufferdata
write buffer data (largest residual, Chi2/ndf, per thread)
real(mpd), dimension(:), allocatable workspacediag
diagonal of global matrix (for global corr.)
integer(mpi) lenglobalvec
length of global vector 'b' (A*x=b)
real(mps) stepl
step length (line search)
integer(mpi) msngpe
upper bound for pair entry single precision storage
real(mpd), dimension(:), allocatable vecbav
vector b for AVPROD (A*x=b)
integer(mpi), dimension(:), allocatable readbufferpointer
pointer to used buffers
integer(mpi), dimension(:), allocatable workspacei
(general) workspace (I)
integer(mpi), dimension(:), allocatable globalparcons
global parameters (number of) constraints
integer(mpi), dimension(:,:), allocatable writebufferinfo
write buffer management (per thread)
integer(mpi) matrit
matrix calculation up to iteration MATRIT
real(mpd), dimension(:), allocatable vbnd
local fit band part of 'A'
real(mpr4), dimension(:), allocatable readbufferdataf
float data
integer(mpi) lfitbb
local fit: check for bordered band matrix (if >0)
integer(mpi) lsearch
iterations (solutions) with line search: >2: all, =2: all with (next) Chi2 cut scaling factor =1....
integer(mpi), dimension(:), allocatable dfd
file: ndf sum
type(listitem), dimension(:), allocatable listpresigmas
list of pre-sgmas from steering file
integer(mpi), dimension(:), allocatable sparsematrixcompression
compression info (per row)
integer(mpi) mrmode
MINRES-QLP mode (0: QR+QLP, 1: only QR, 2: only QLP factorization)
real(mps) chicut
cut in terms of 3-sigma cut, first iteration
integer(mpi) imonit
flag for monitoring residuals per local fit cycle (=0: none, <0: all, bit 0: first,...
real(mps), dimension(nplan) dvd
rel. drift velocity deviation (calibration parameter)
real(mps), dimension(nplan) del
shift (position deviation) (alignment parameter)
integer(mpi), parameter nplan
integer(mpi), parameter nmx
number of modules in x direction
real(mps), dimension(ntot) sdevx
shift in x (alignment parameter)
real(mps), dimension(ntot) sdevy
shift in y (alignment parameter)
integer(mpi), parameter nmy
number of modules in y direction
integer(mpi), parameter nlyr
number of detector layers
integer(mpi), parameter ntot
total number of modules
integer(mpi) keyb
end (position) of keyword
integer(mpi) keya
start (position) of keyword
subroutine ploopb(lunp)
Print iteration line.
subroutine bincls(kfile, ithr)
Close binary file.
subroutine prpcon
Prepare constraints.
subroutine mminrs
Solution with MINRES.
subroutine mcsolv(n, x, y)
Solution for zero band width preconditioner.
subroutine mupdat(i, j, add)
Update element of global matrix.
subroutine peend(icode, cmessage)
Print exit code.
subroutine loopn
Loop with fits and sums.
subroutine loop1
First data loop (get global labels).
subroutine feasma
Matrix for feasible solution.
subroutine xloopn
Standard solution algorithm.
subroutine ploopa(lunp)
Print title for iteration.
subroutine isjajb(nst, is, ja, jb, jsp)
Decode Millepede record.
subroutine additem(length, list, label, value)
add item to list
subroutine binrwd(kfile)
Rewind binary file.
subroutine zdiags
Covariance matrix for diagonalization (,correction of eigenvectors).
subroutine solglo(ivgbi)
Error for single global parameter from MINRES.
subroutine upone
Update, redefine hash indices.
subroutine prtglo
Print final log file.
subroutine monres
Monitor input residuals.
subroutine intext(text, nline)
Interprete text.
integer(mpl) function ijadd(itema, itemb)
Index for sparse storage.
subroutine mdiags
Solution by diagonalization.
program mptwo
Millepede II main program Pede.
subroutine prtstat
Print input statistic.
subroutine getsum(asum)
Get accurate sum.
subroutine peread(more)
Read (block of) records from binary files.
subroutine filetx
Interprete text files.
integer(mpi) function iprime(n)
largest prime number < N.
subroutine ploopc(lunp)
Print sub-iteration line.
subroutine mvopen(lun, fname)
Open file.
subroutine solgloqlp(ivgbi)
Error for single global parameter from MINRES-QLP.
subroutine addcst
Add constraint information to matrix and vector.
subroutine avprd0(n, x, b)
Product symmetric matrix times vector.
subroutine petime
Print times.
subroutine mstart(text)
Start of 'module' printout.
subroutine mend
End of 'module' printout.
subroutine minver
Solution by matrix inversion.
subroutine peprep(mode)
Prepare records.
subroutine explfc(lunit)
Print explanation of iteration table.
subroutine binopn(kfile, ithr, ierr)
Open binary file.
subroutine addsum(add)
Accurate summation.
subroutine sechms(deltat, nhour, minut, secnd)
Time conversion.
integer(mpi) function inone(item)
Translate labels to indices (for global parameters).
subroutine avprod(n, x, b)
Product symmetric matrix times vector.
subroutine loop1i
Iteration of first data loop.
subroutine mhalf2
Fill 2nd half of matrix for extended storage.
subroutine mminrsqlp
Solution with MINRES-QLP.
subroutine filetc
Interprete command line option, steering file.
subroutine feasib(concut, iact)
Make parameters feasible.
subroutine mvsolv(n, x, y)
Solution for finite band width preconditioner.
subroutine loopbf(nrej, ndfs, sndf, dchi2s, numfil, naccf, chi2f, ndff)
Loop over records in read buffer (block), fits and sums.
subroutine vmprep(msize)
Prepare storage for vectors and matrices.
subroutine ploopd(lunp)
Print solution line.
subroutine pechk(ibuf, nerr)
Check Millepede record.
subroutine loop2
Second data loop (number of derivatives, global label pairs).
integer(mpi) function nufile(fname)
Inquire on file.
list items from steering file