906 REAL,
DIMENSION(2) :: ta
909 INTEGER(mpi) :: iopnmp
920 INTEGER(mpi) :: nsecnd
921 INTEGER(mpi) :: ntsec
923 CHARACTER (LEN=24) :: chdate
924 CHARACTER (LEN=24) :: chost
926 CHARACTER (LEN=6) :: c6
927 INTEGER major, minor, patch
951 CALL getenv(
'HOSTNAME',chost)
952 IF (chost(1:1) ==
' ')
CALL getenv(
'HOST',chost)
953 WRITE(*,*)
'($Id: c6be33cb87ef25e4226d5088ec53ebf6d39f4712 $)'
958 CALL ilaver( major,minor, patch )
959 WRITE(*,110) lapack64, major,minor, patch
960110
FORMAT(
' using LAPACK64 with ',(a),
', version ',i0,
'.',i0,
'.',i0)
962 WRITE(*,*)
'using Intel oneMKL PARDISO'
966 WRITE(*,111) __gnuc__ , __gnuc_minor__ , __gnuc_patchlevel__
967111
FORMAT(
' compiled with gcc ',i0,
'.',i0,
'.',i0)
970 WRITE(*,111) __pgic__ , __pgic_minor__ , __pgic_patchlevel__
971111
FORMAT(
' compiled with pgi ',i0,
'.',i0,
'.',i0)
974 WRITE(*,*)
' < Millepede II-P starting ... ',chdate
978 WRITE(8,*)
'($Id: c6be33cb87ef25e4226d5088ec53ebf6d39f4712 $)'
980 WRITE(8,*)
'Log-file Millepede II-P ', chdate
981 WRITE(8,*)
' ', chost
983 CALL peend(-1,
'Still running or crashed')
992 WRITE(*,*)
'!!! Checking input only, no calculation of a solution !!!'
993 WRITE(8,*)
'!!! Checking input only, no calculation of a solution !!!'
1012 CALL getenv(
'OMP_NUM_THREADS',c6)
1014 CALL getenv(lapack64//
'_NUM_THREADS',c6)
1016 IF (c6(1:1) ==
' ')
THEN
1018 WRITE(*,*)
'Number of threads for LAPACK: unkown (empty OMP_NUM_THREADS)'
1020 WRITE(*,*)
'Number of threads for LAPACK: unkown (empty ',lapack64//
'_NUM_THREADS)'
1023 WRITE(*,*)
'Number of threads for LAPACK: ', c6
1043 CALL mvopen(lun,
'millepede.his')
1049 CALL mvopen(1,
'mpdebug.txt')
1053 times(0)=rstext-rstp
1059 times(1)=rloop1-rstext
1063 WRITE(8,*)
'Chi square cut equiv 3 st.dev applied ...'
1064 WRITE(8,*)
' in first iteration with factor',
chicut
1065 WRITE(8,*)
' in second iteration with factor',
chirem
1066 WRITE(8,*)
' (reduced by sqrt in next iterations)'
1070 WRITE(8,*)
'Down-weighting of outliers in',
lhuber,
' iterations'
1071 WRITE(8,*)
'Cut on downweight fraction',
dwcut
1075 times(2)=rloop2-rloop1
1080 CALL peend(5,
'Ended without solution (empty constraints)')
1082 CALL peend(0,
'Ended normally')
1119 CALL gmpdef(6,1,
'log10(#records) vs file number')
1120 CALL gmpdef(7,1,
'final rejection fraction vs file number')
1122 'final <Chi^2/Ndf> from accepted local fits vs file number')
1123 CALL gmpdef(9,1,
'<Ndf> from accepted local fits vs file number')
1129 rej=real(nrc-
jfd(kfl),mps)/real(nrc,mps)
1130 CALL gmpxy(6,real(kfl,mps),log10(real(nrc,mps)))
1131 CALL gmpxy(7,real(kfl,mps),rej)
1133 IF (
jfd(kfl) > 0)
THEN
1134 c2ndf=
cfd(kfl)/real(
jfd(kfl),mps)
1135 CALL gmpxy(8,real(kfl,mps),c2ndf)
1136 andf=real(
dfd(kfl),mps)/real(
jfd(kfl),mps)
1137 CALL gmpxy(9,real(kfl,mps),andf)
1154 WRITE(*,*)
'Misalignment test wire chamber'
1157 CALL hmpdef( 9,-0.0015,+0.0015,
'True - fitted displacement')
1158 CALL hmpdef(10,-0.0015,+0.0015,
'True - fitted Vdrift')
1164 sums(1)=sums(1)+diff
1165 sums(2)=sums(2)+diff*diff
1167 sums(3)=sums(3)+diff
1168 sums(4)=sums(4)+diff*diff
1170 sums(1)=0.01_mpd*sums(1)
1171 sums(2)=sqrt(0.01_mpd*sums(2))
1172 sums(3)=0.01_mpd*sums(3)
1173 sums(4)=sqrt(0.01_mpd*sums(4))
1174 WRITE(*,143)
'Parameters 1 - 100: mean =',sums(1),
'rms =',sums(2)
1175 WRITE(*,143)
'Parameters 101 - 200: mean =',sums(3),
'rms =',sums(4)
1176143
FORMAT(6x,a28,f9.6,3x,a5,f9.6)
1179 WRITE(*,*)
' I label simulated fitted diff'
1180 WRITE(*,*)
' -------------------------------------------- '
1200 WRITE(*,*)
'Misalignment test Si tracker'
1203 CALL hmpdef( 9,-0.0025,+0.0025,
'True - fitted displacement X')
1204 CALL hmpdef(10,-0.025,+0.025,
'True - fitted displacement Y')
1215 sums(1)=sums(1)+1.0_mpd
1216 sums(2)=sums(2)+diff
1217 sums(3)=sums(3)+diff*diff
1222 err=sqrt(abs(gmati))
1224 sums(7)=sums(7)+1.0_mpd
1225 sums(8)=sums(8)+diff
1226 sums(9)=sums(9)+diff*diff
1229 IF (mod(i,3) == 1)
THEN
1233 sums(4)=sums(4)+1.0_mpd
1234 sums(5)=sums(5)+diff
1235 sums(6)=sums(6)+diff*diff
1240 err=sqrt(abs(gmati))
1242 sums(7)=sums(7)+1.0_mpd
1243 sums(8)=sums(8)+diff
1244 sums(9)=sums(9)+diff*diff
1249 sums(2)=sums(2)/sums(1)
1250 sums(3)=sqrt(sums(3)/sums(1))
1251 sums(5)=sums(5)/sums(4)
1252 sums(6)=sqrt(sums(6)/sums(4))
1253 WRITE(*,143)
'Parameters 1 - 500: mean =',sums(2),
'rms =',sums(3)
1254 WRITE(*,143)
'Parameters 501 - 700: mean =',sums(5),
'rms =',sums(6)
1255 IF (sums(7) > 0.5_mpd)
THEN
1256 sums(8)=sums(8)/sums(7)
1257 sums(9)=sqrt(sums(9)/sums(7))
1258 WRITE(*,143)
'Parameter pulls, all: mean =',sums(8),
'rms =',sums(9)
1262 WRITE(*,*)
' I label simulated fitted diff'
1263 WRITE(*,*)
' -------------------------------------------- '
1275 IF (mod(i,3) == 1)
THEN
1295 WRITE(8,*)
'Record',
nrec1,
' has largest residual:',
value1
1298 WRITE(8,*)
'Record',
nrec2,
' has largest Chi^2/Ndf:',
value2
1302 WRITE(8,*)
'Record',
nrec3,
' is first with error (rank deficit/NaN)'
1306 WRITE(8,*)
'In total 3 +',
nloopn,
' loops through the data files'
1308 WRITE(8,*)
'In total 2 +',
nloopn,
' loops through the data files'
1312 WRITE(8,*)
'In total ',
mnrsit,
' internal MINRES iterations'
1320 ntsec=nint(deltat,mpi)
1321 CALL sechms(deltat,nhour,minut,secnd)
1322 nsecnd=nint(secnd,mpi)
1323 WRITE(8,*)
'Total time =',ntsec,
' seconds =',nhour,
' h',minut, &
1324 ' m',nsecnd,
' seconds'
1326 WRITE(8,*)
'end ', chdate
1335 WRITE(8,*)
'Data rejected in last iteration: '
1337 nrejec(0),
' (rank deficit/NaN) ',
nrejec(1),
' (Ndf=0) ', &
1344 WRITE(*,*)
' < Millepede II-P ending ... ', chdate
1351106
FORMAT(
' PARDISO dyn. memory allocation: ',f11.6,
' GB')
1356102
FORMAT(2x,i4,i10,2x,3f10.5)
1357103
FORMAT(
' Times [in sec] for text processing',f12.3/ &
1360 ' func. value ',f12.3,
' *',f4.0/ &
1361 ' func. value, global matrix, solution',f12.3,
' *',f4.0/ &
1362 ' new solution',f12.3,
' *',f4.0/)
1363105
FORMAT(
' Peak dynamic memory allocation: ',f11.6,
' GB')
1383 INTEGER(mpi) :: istop
1384 INTEGER(mpi) :: itgbi
1385 INTEGER(mpi) :: itgbl
1387 INTEGER(mpi) :: itnlim
1388 INTEGER(mpi) :: nout
1390 INTEGER(mpi),
INTENT(IN) :: ivgbi
1427 globalcorrections, itnlim, nout, rtol, istop, itn, anorm, acond, rnorm, arnorm, ynorm)
1431 globalcorrections, itnlim, nout, rtol, istop, itn, anorm, acond, rnorm, arnorm, ynorm)
1434 globalcorrections, itnlim, nout, rtol, istop, itn, anorm, acond, rnorm, arnorm, ynorm)
1440 err=sqrt(abs(real(gmati,mps)))
1441 IF(gmati < 0.0_mpd) err=-err
1442 diag=matij(ivgbi,ivgbi)
1443 gcor2=real(1.0_mpd-1.0_mpd/(gmati*diag),mps)
1445101
FORMAT(1x,
' label parameter presigma differ', &
1446 ' Error gcor^2 iit'/ 1x,
'---------',2x,5(
'-----------'),2x,
'----')
1447102
FORMAT(i10,2x,4g12.4,f7.4,i6,i4)
1467 INTEGER(mpi) :: istop
1468 INTEGER(mpi) :: itgbi
1469 INTEGER(mpi) :: itgbl
1471 INTEGER(mpi) :: itnlim
1472 INTEGER(mpi) :: nout
1474 INTEGER(mpi),
INTENT(IN) :: ivgbi
1503 mxxnrm = real(
nagb,mpd)/sqrt(epsilon(mxxnrm))
1505 trcond = 1.0_mpd/epsilon(trcond)
1506 ELSE IF(
mrmode == 2)
THEN
1514 itnlim=itnlim, rtol=rtol, maxxnorm=mxxnrm, trancond=trcond, &
1518 itnlim=itnlim, rtol=rtol, maxxnorm=mxxnrm, trancond=trcond, &
1522 itnlim=itnlim, rtol=rtol, maxxnorm=mxxnrm, trancond=trcond, &
1529 err=sqrt(abs(real(gmati,mps)))
1530 IF(gmati < 0.0_mpd) err=-err
1531 diag=matij(ivgbi,ivgbi)
1532 gcor2=real(1.0_mpd-1.0_mpd/(gmati*diag),mps)
1534101
FORMAT(1x,
' label parameter presigma differ', &
1535 ' Error gcor^2 iit'/ 1x,
'---------',2x,5(
'-----------'),2x,
'----')
1536102
FORMAT(i10,2x,4g12.4,f7.4,i6,i4)
1549 INTEGER(mpi) :: icgb
1550 INTEGER(mpi) :: irhs
1551 INTEGER(mpi) :: itgbi
1552 INTEGER(mpi) :: ivgb
1554 INTEGER(mpi) :: jcgb
1556 INTEGER(mpi) :: label
1558 INTEGER(mpi) :: inone
1561 REAL(mpd) :: drhs(4)
1562 INTEGER(mpi) :: idrh (4)
1587 IF(abs(rhs) > climit)
THEN
1593 WRITE(*,101) (idrh(l),drhs(l),l=1,irhs)
1602 WRITE(*,101) (idrh(l),drhs(l),l=1,irhs)
1605 WRITE(*,102)
' Constraints: only equation values >', climit,
' are printed'
1606101
FORMAT(
' ',4(i6,g11.3))
1607102
FORMAT(a,g11.2,a)
1620 INTEGER(mpi) :: icgb
1621 INTEGER(mpi) :: icgrp
1622 INTEGER(mpi) :: ioff
1623 INTEGER(mpi) :: itgbi
1625 INTEGER(mpi) :: jcgb
1626 INTEGER(mpi) :: label
1627 INTEGER(mpi) :: labelf
1628 INTEGER(mpi) :: labell
1629 INTEGER(mpi) :: last
1630 INTEGER(mpi) :: line1
1631 INTEGER(mpi) :: ncon
1632 INTEGER(mpi) :: ndiff
1633 INTEGER(mpi) :: npar
1634 INTEGER(mpi) :: inone
1635 INTEGER(mpi) :: itype
1636 INTEGER(mpi) :: ncgbd
1637 INTEGER(mpi) :: ncgbr
1638 INTEGER(mpi) :: ncgbw
1639 INTEGER(mpi) :: ncgrpd
1640 INTEGER(mpi) :: ncgrpr
1641 INTEGER(mpi) :: next
1643 INTEGER(mpl):: length
1644 INTEGER(mpl) :: rows
1646 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: vecParConsOffsets
1647 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: vecParConsList
1648 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: vecConsParOffsets
1649 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: vecConsParList
1650 INTEGER(mpi),
DIMENSION(:,:),
ALLOCATABLE :: matConsGroupIndex
1663 IF(last < 0.AND.label < 0)
THEN
1666 IF(itype == 2) ncgbw=ncgbw+1
1673 IF(label > 0.AND.itype == 2)
THEN
1680 IF (ncgbw == 0)
THEN
1681 WRITE(*,*)
'GRPCON:',
ncgb,
' constraints found in steering files'
1683 WRITE(*,*)
'GRPCON:',
ncgb,
' constraints found in steering files,',ncgbw,
'weighted'
1688 length=
ncgb+1; rows=3
1701 CALL mpalloc(matconsgroupindex,rows,length,
'group index for constraint (I)')
1704 CALL mpalloc(vecconsparoffsets,length,
'offsets for global par list for cons. (I)')
1706 CALL mpalloc(vecparconsoffsets,length,
'offsets for cons. list for global par. (I)')
1707 vecparconsoffsets(1)=0
1709 vecparconsoffsets(i+1)=vecparconsoffsets(i)+
globalparcons(i)
1713 length=vecparconsoffsets(
ntgb+1)
1714 CALL mpalloc(vecconsparlist,length,
'global par. list for constraint (I)')
1715 CALL mpalloc(vecparconslist,length,
'constraint list for global par. (I)')
1720 vecconsparoffsets(1)=ioff
1732 vecparconslist(vecparconsoffsets(itgbi)+
globalparcons(itgbi))=icgb
1734 vecconsparlist(ioff+npar)=itgbi
1740 CALL sort1k(vecconsparlist(ioff+1),npar)
1744 next=vecconsparlist(ioff+j)
1745 IF (next /= last)
THEN
1747 vecconsparlist(ioff+ndiff) = next
1756 vecconsparoffsets(icgb+1)=ioff
1769 print *,
' Constraint #parameters first par. last par. first line'
1777 npar=vecconsparoffsets(icgb+1)-vecconsparoffsets(icgb)
1781 print *, jcgb, npar, labelf, labell, line1
1784 icgrp=matconsgroupindex(1,icgb)
1785 IF (icgrp == 0)
THEN
1787 DO i=vecconsparoffsets(icgb)+1, vecconsparoffsets(icgb+1)
1788 itgbi=vecconsparlist(i)
1790 DO j=vecparconsoffsets(itgbi)+1,vecparconsoffsets(itgbi+1)
1791 icgrp=matconsgroupindex(1,vecparconslist(j))
1797 IF (icgrp == 0)
THEN
1804 matconsgroupindex(2,icgb)=jcgb
1805 matconsgroupindex(3,icgb)=icgb
1806 DO i=vecconsparoffsets(icgb)+1, vecconsparoffsets(icgb+1)
1807 itgbi=vecconsparlist(i)
1810 DO j=vecparconsoffsets(itgbi)+1,vecparconsoffsets(itgbi+1)
1811 matconsgroupindex(1,vecparconslist(j))=icgrp
1815 WRITE(*,*)
'GRPCON:',
ncgrp,
' disjoint constraints groups built'
1823 icgb=matconsgroupindex(3,jcgb)
1828 icgrp=matconsgroupindex(1,jcgb)
1848 print *,
' cons.group first con. first par. last par. #cons #par'
1859 print *, icgrp,
matconsgroups(1,icgrp), labelf, labell, ncon, npar
1862 IF (ncon == npar)
THEN
1869 print *, icgrp,
matconsgroups(1,icgrp), labelf, labell,
' : cons.group resolved'
1884 print *, icgrp,
matconsgroups(1,icgrp), labelf, labell,
' : cons.group redundant'
1889 IF (ncgrpr > 0)
THEN
1890 WRITE(*,*)
'GRPCON:',ncgbr,
' redundancy constraints in ', ncgrpr,
' groups resolved'
1894 IF (ncgrpd > 0)
THEN
1895 WRITE(*,*)
'GRPCON:',ncgbd,
' redundancy constraints in ', ncgrpd,
' groups detected'
1918 INTEGER(mpi) :: icgb
1919 INTEGER(mpi) :: icgrp
1920 INTEGER(mpi) :: ifrst
1921 INTEGER(mpi) :: ilast
1922 INTEGER(mpi) :: isblck
1923 INTEGER(mpi) :: itgbi
1924 INTEGER(mpi) :: ivgb
1926 INTEGER(mpi) :: jcgb
1927 INTEGER(mpi) :: jfrst
1928 INTEGER(mpi) :: label
1929 INTEGER(mpi) :: labelf
1930 INTEGER(mpi) :: labell
1931 INTEGER(mpi) :: ncon
1932 INTEGER(mpi) :: ngrp
1933 INTEGER(mpi) :: npar
1934 INTEGER(mpi) :: ncnmxb
1935 INTEGER(mpi) :: ncnmxg
1936 INTEGER(mpi) :: nprmxb
1937 INTEGER(mpi) :: nprmxg
1938 INTEGER(mpi) :: inone
1939 INTEGER(mpi) :: nvar
1941 INTEGER(mpl):: length
1942 INTEGER(mpl) :: rows
1944 INTEGER(mpi),
DIMENSION(:,:),
ALLOCATABLE :: matConsGroupIndex
1957 length=
ncgrp+1; rows=3
1962 CALL mpalloc(matconsgroupindex,rows,length,
'group index for constraint (I)')
1999 IF (nvar > 0 .OR.
iskpec == 0)
THEN
2002 matconsgroupindex(1,
ncgb)=ngrp+1
2003 matconsgroupindex(2,
ncgb)=icgb
2004 matconsgroupindex(3,
ncgb)=nvar
2007 IF (
ncgb > ncon) ngrp=ngrp+1
2013 WRITE(*,*)
'PRPCON:',
ncgbe,
' empty constraints skipped'
2015 WRITE(*,*)
'PRPCON:',
ncgbe,
' empty constraints detected, to be fixed !!!'
2016 WRITE(*,*)
' (use option "skipemptycons" to skip those)'
2021 WRITE(*,*)
'!!! Switch to "-C" (checking input only), no calculation of a solution !!!'
2022 WRITE(8,*)
'!!! Switch to "-C" (checking input only), no calculation of a solution !!!'
2027 WRITE(*,*)
'PRPCON:',
ncgb,
' constraints accepted'
2030 IF(
ncgb == 0)
RETURN
2037 icgb=matconsgroupindex(2,jcgb)
2042 icgrp=matconsgroupindex(1,jcgb)
2068 WRITE(*,*)
' Cons. sorted index #var.par. first line first label last label'
2069 WRITE(*,*)
' Cons. group index first cons. last cons. first label last label'
2070 WRITE(*,*)
' Cons. block index first group last group first label last label'
2076 nvar=matconsgroupindex(3,jcgb)
2080 WRITE(*,*)
' Cons. sorted', jcgb, nvar, &
2083 WRITE(*,*)
' Cons. sorted', jcgb,
' empty (0)', &
2092 WRITE(*,*)
' Cons. group ', icgrp,
matconsgroups(1,icgrp), &
2105 DO i=icgrp,isblck,-1
2119 WRITE(*,*)
' Cons. block ',
ncblck, isblck, icgrp, labelf, labell
2137 ncnmxb=max(ncnmxb,ncon)
2138 nprmxb=max(nprmxb,npar)
2163 ncnmxg=max(ncnmxg,ncon)
2164 nprmxg=max(nprmxg,npar)
2199 WRITE(*,*)
'PRPCON: constraints split into ',
ncgrp,
'(disjoint) groups,'
2200 WRITE(*,*)
' groups combined into ',
ncblck,
'(non overlapping) blocks'
2201 WRITE(*,*)
' max group size (cons., par.) ', ncnmxg, nprmxg
2202 WRITE(*,*)
' max block size (cons., par.) ', ncnmxb, nprmxb
2220 INTEGER(mpi) :: icgb
2221 INTEGER(mpi) :: icgrp
2223 INTEGER(mpi) :: ifirst
2224 INTEGER(mpi) :: ilast
2225 INTEGER(mpl) :: ioffc
2226 INTEGER(mpl) :: ioffp
2227 INTEGER(mpi) :: irank
2228 INTEGER(mpi) :: ipar0
2229 INTEGER(mpi) :: itgbi
2230 INTEGER(mpi) :: ivgb
2232 INTEGER(mpi) :: jcgb
2234 INTEGER(mpi) :: label
2235 INTEGER(mpi) :: ncon
2236 INTEGER(mpi) :: npar
2237 INTEGER(mpi) :: nrank
2238 INTEGER(mpi) :: inone
2243 INTEGER(mpl):: length
2244 REAL(mpd),
DIMENSION(:),
ALLOCATABLE :: matConstraintsT
2245 REAL(mpd),
DIMENSION(:),
ALLOCATABLE :: auxVectorD
2246 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: auxVectorI
2250 IF(
ncgb == 0)
RETURN
2259 CALL mpalloc(auxvectori,length,
'auxiliary array (I)')
2260 CALL mpalloc(auxvectord,length,
'auxiliary array (D)')
2263 CALL mpalloc(matconstraintst,length,
'transposed matrix of constraints (blocks)')
2264 matconstraintst=0.0_mpd
2277 WRITE(*,*)
' Constraint group, #con, rank', icgrp, ncon, 0,
' (empty)'
2280 DO jcgb=ifirst,ilast
2292 IF(ivgb > 0) matconstraintst(int(jcgb-ifirst,mpl)*int(npar,mpl)+ivgb-ipar0+ioffc)= &
2293 matconstraintst(int(jcgb-ifirst,mpl)*int(npar,mpl)+ivgb-ipar0+ioffc)+factr
2300 DO ll=ioffc+1,ioffc+npar
2306 matconstraintst(int(i-1,mpl)*int(npar,mpl)+ll)* &
2307 matconstraintst(int(j-1,mpl)*int(npar,mpl)+ll)
2313 IF (
icheck > 1 .OR. irank < ncon)
THEN
2314 WRITE(*,*)
' Constraint group, #con, rank', icgrp, ncon, irank
2315 IF (irank < ncon)
THEN
2316 WRITE(*,*)
' .. rank deficit !! '
2317 WRITE(*,*)
' E.g. fix all parameters and remove all constraints related to label ', &
2322 ioffc=ioffc+int(npar,mpl)*int(ncon,mpl)
2329 WRITE(*,*)
'Rank of product matrix of constraints is',nrank, &
2330 ' for',
ncgb,
' constraint equations'
2331 WRITE(8,*)
'Rank of product matrix of constraints is',nrank, &
2332 ' for',
ncgb,
' constraint equations'
2333 IF(nrank <
ncgb)
THEN
2334 WRITE(*,*)
'Warning: insufficient constraint equations!'
2335 WRITE(8,*)
'Warning: insufficient constraint equations!'
2338 WRITE(*,*)
' --> enforcing SUBITO mode'
2339 WRITE(8,*)
' --> enforcing SUBITO mode'
2346 print *,
'QL decomposition of constraints matrix'
2349 WRITE(
lunlog,*)
'QL decomposition of constraints matrix'
2361 CALL lpqldec(matconstraintst,evmin,evmax)
2365 print *,
' largest |eigenvalue| of L: ', evmax
2366 print *,
' smallest |eigenvalue| of L: ', evmin
2367 IF (evmin == 0.0_mpd.AND.
icheck == 0)
THEN
2368 CALL peend(27,
'Aborted, singular QL decomposition of constraints matrix')
2369 stop
'FEASMA: stopping due to singular QL decomposition of constraints matrix'
2395 INTEGER(mpi) :: icgb
2396 INTEGER(mpi) :: icgrp
2397 INTEGER(mpi) :: iter
2398 INTEGER(mpi) :: itgbi
2399 INTEGER(mpi) :: ivgb
2400 INTEGER(mpi) :: ieblck
2401 INTEGER(mpi) :: isblck
2402 INTEGER(mpi) :: ifirst
2403 INTEGER(mpi) :: ilast
2405 INTEGER(mpi) :: jcgb
2406 INTEGER(mpi) :: label
2407 INTEGER(mpi) :: inone
2408 INTEGER(mpi) :: ncon
2410 REAL(mps),
INTENT(IN) :: concut
2411 INTEGER(mpi),
INTENT(OUT) :: iact
2418 REAL(mpd),
DIMENSION(:),
ALLOCATABLE :: vecCorrections
2422 IF(
ncgb == 0)
RETURN
2452 sum1=sqrt(sum1/real(
ncgb,mpd))
2453 sum2=sum2/real(
ncgb,mpd)
2455 IF(iter == 1.AND.sum1 < concut)
RETURN
2457 IF(iter == 1.AND.
ncgb <= 12)
THEN
2459 WRITE(*,*)
'Constraint equation discrepancies:'
2461101
FORMAT(4x,4(i5,g12.4))
2463103
FORMAT(10x,
' Cut on rms value is',g8.1)
2468 WRITE(*,*)
'Improve constraints'
2472 WRITE(*,102) iter,sum1,sum2,sum3
2473102
FORMAT(i6,
' rms',g12.4,
' avrg_abs',g12.4,
' max_abs',g12.4)
2475 CALL mpalloc(veccorrections,int(
nvgb,mpl),
'constraint corrections')
2476 veccorrections=0.0_mpd
2484 ieblck=isblck+(ncon*(ncon+1))/2
2553 INTEGER(mpi) :: iact
2554 INTEGER(mpi) :: ierrc
2555 INTEGER(mpi) :: ierrf
2556 INTEGER(mpi) :: ioffp
2558 INTEGER(mpi) :: ithr
2559 INTEGER(mpi) :: jfile
2560 INTEGER(mpi) :: jrec
2562 INTEGER(mpi) :: kfile
2565 INTEGER(mpi) :: mpri
2567 INTEGER(mpi) :: nact
2568 INTEGER(mpi) :: nbuf
2569 INTEGER(mpi) :: ndata
2570 INTEGER(mpi) :: noff
2571 INTEGER(mpi) :: noffs
2572 INTEGER(mpi) :: npointer
2573 INTEGER(mpi) :: npri
2577 INTEGER(mpi) :: nrpr
2578 INTEGER(mpi) :: nthr
2579 INTEGER(mpi) :: ntot
2580 INTEGER(mpi) :: maxRecordSize
2581 INTEGER(mpi) :: maxRecordFile
2583 INTEGER(mpi),
INTENT(OUT) :: more
2593 CHARACTER (LEN=7) :: cfile
2598 SUBROUTINE readc(bufferD, bufferF, bufferI, bufferLength, lun, err)
BIND(c)
2600 REAL(c_double),
DIMENSION(*),
INTENT(OUT) :: bufferD
2601 REAL(c_float),
DIMENSION(*),
INTENT(OUT) :: bufferF
2602 INTEGER(c_int),
DIMENSION(*),
INTENT(OUT) :: bufferI
2603 INTEGER(c_int),
INTENT(INOUT) :: bufferLength
2604 INTEGER(c_int),
INTENT(IN),
VALUE :: lun
2605 INTEGER(c_int),
INTENT(OUT) :: err
2606 END SUBROUTINE readc
2612 DATA npri / 0 /, mpri / 1000 /
2661 files:
DO WHILE (jfile > 0)
2665 CALL binopn(kfile,ithr,ios)
2671 IF(kfile <=
nfilf)
THEN
2692 eof=(ierrc <= 0.AND.ierrc /= -4)
2693 IF(eof.AND.ierrc < 0)
THEN
2694 WRITE(*,*)
'Read error for binary Cfile', kfile,
'record', jrec+1,
':', ierrc
2695 WRITE(8,*)
'Read error for binary Cfile', kfile,
'record', jrec+1,
':', ierrc
2697 WRITE(cfile,
'(I7)') kfile
2698 CALL peend(18,
'Aborted, read error(s) for binary file ' // cfile)
2699 stop
'PEREAD: stopping due to read errors'
2701 IF (
kfd(1,jfile) == 1)
THEN
2707 IF(eof)
EXIT records
2712 xfd(jfile)=max(
xfd(jfile),n)
2729 IF ((noff+nr+2+
ndimbuf >= ndata*(iact+1)).OR.(nbuf >= npointer))
EXIT files
2744 IF (
kfd(1,jfile) == 1)
THEN
2745 print *,
'PEREAD: file ', kfile,
'read the first time, found',jrec,
' records'
2749 IF (-
kfd(1,jfile) /= jrec)
THEN
2750 WRITE(cfile,
'(I7)') kfile
2751 CALL peend(19,
'Aborted, binary file modified (length) ' // cfile)
2752 stop
'PEREAD: file modified (length)'
2803 IF (
kfd(1,jfile) == 1)
kfd(1,jfile)=-nrc
2822 DO WHILE (
nloopn == 0.AND.ntot >= nrpr)
2823 WRITE(*,*)
' Record ',nrpr
2824 IF (nrpr < 100000)
THEN
2833 IF (npri == 1)
WRITE(*,100)
2835100
FORMAT(/
' PeRead records active file' &
2836 /
' total block threads number')
2837101
FORMAT(
' PeRead',4i10)
2850 IF (
xfd(k) > maxrecordsize)
THEN
2851 maxrecordsize=
xfd(k)
2854 dw=real(-
kfd(1,k),mpd)
2857 ds1=ds1+dw*real(
wfd(k),mpd)
2858 ds2=ds2+dw*real(
wfd(k)**2,mpd)
2860 print *,
'PEREAD: file ', maxrecordfile,
'with max record size ', maxrecordsize
2861 IF (
nfilw > 0.AND.ds0 > 0.0_mpd)
THEN
2865 WRITE(lun,177)
nfilw,real(ds1,mps),real(ds2,mps)
2866177
FORMAT(/
' !!!!!',i4,
' weighted binary files', &
2867 /
' !!!!! mean, variance of weights =',2g12.4)
2885179
FORMAT(/
' Read cache usage (#blocks, #records, ', &
2886 'min,max records/block'/17x,i10,i12,2i10)
2904 INTEGER(mpi),
INTENT(IN) :: mode
2906 INTEGER(mpi) :: ibuf
2907 INTEGER(mpi) :: ichunk
2909 INTEGER(mpi) :: itgbi
2915 INTEGER(mpi),
PARAMETER :: maxbad = 100
2916 INTEGER(mpi) :: nbad
2917 INTEGER(mpi) :: nerr
2918 INTEGER(mpi) :: inone
2936 CALL isjajb(nst,ist,ja,jb,jsp)
2958 CALL pechk(ibuf,nerr)
2961 IF(nbad >= maxbad)
EXIT
2966 CALL isjajb(nst,ist,ja,jb,jsp)
2979 CALL peend(20,
'Aborted, bad binary records')
2980 stop
'PEREAD: stopping due to bad records'
3001 INTEGER(mpi) :: ioff
3008 INTEGER(mpi),
INTENT(IN) :: ibuf
3009 INTEGER(mpi),
INTENT(OUT) :: nerr
3018 outer:
DO WHILE(is < nst)
3023 IF(is > nst)
EXIT outer
3029 IF(is > nst)
EXIT outer
3046100
FORMAT(
' PEREAD: record ', i8,
' in file ',i6,
' is broken !!!')
3056101
FORMAT(
' PEREAD: record ', i8,
' in file ',i6,
' contains ', i6,
' NaNs !!!')
3072 INTEGER(mpi) :: ibuf
3073 INTEGER(mpi) :: ichunk
3074 INTEGER(mpi) :: iproc
3075 INTEGER(mpi) :: ioff
3076 INTEGER(mpi) :: ioffbi
3078 INTEGER(mpi) :: itgbi
3083 INTEGER(mpi) :: nalg
3085 INTEGER(mpi) :: inone
3086 INTEGER(mpl) :: length
3120 CALL isjajb(nst,ist,ja,jb,jsp)
3144 CALL isjajb(nst,ist,ja,jb,jsp)
3165 CALL isjajb(nst,ist,ja,jb,jsp)
3195 INTEGER(mpi) :: istart
3196 INTEGER(mpi) :: itgbi
3198 INTEGER(mpi) :: jstart
3199 INTEGER(mpi) :: jtgbi
3200 INTEGER(mpi) :: lstart
3201 INTEGER(mpi) :: ltgbi
3203 INTEGER(mpi),
INTENT(IN) :: inds
3204 INTEGER(mpi),
INTENT(IN) :: inde
3206 IF (inds > inde)
RETURN
3216 IF (istart == 0)
THEN
3217 IF (itgbi /= ltgbi+1)
THEN
3220 IF (lstart == 0)
THEN
3235 IF (istart /= jstart)
THEN
3246 IF (itgbi /= ltgbi+1)
THEN
3261 IF (istart /= jstart)
THEN
3313 INTEGER(mpi),
INTENT(IN) :: nst
3314 INTEGER(mpi),
INTENT(IN OUT) :: is
3315 INTEGER(mpi),
INTENT(OUT) :: ja
3316 INTEGER(mpi),
INTENT(OUT) :: jb
3317 INTEGER(mpi),
INTENT(OUT) :: jsp
3325 IF(is >= nst)
RETURN
3379 INTEGER(mpi) :: ioffb
3381 INTEGER(mpi) :: itgbi
3382 INTEGER(mpi) :: itgbij
3383 INTEGER(mpi) :: itgbik
3384 INTEGER(mpi) :: ivgb
3385 INTEGER(mpi) :: ivgbij
3386 INTEGER(mpi) :: ivgbik
3389 INTEGER(mpi) :: lastit
3391 INTEGER(mpi) :: ncrit
3392 INTEGER(mpi) :: ngras
3393 INTEGER(mpi) :: nparl
3395 INTEGER(mpi) :: nrej
3396 INTEGER(mpi) :: inone
3397 INTEGER(mpi) :: ilow
3398 INTEGER(mpi) :: nlow
3399 INTEGER(mpi) :: nzero
3419 CALL gmpdef(1,4,
'Function value in iterations')
3421 CALL gmpdef(2,3,
'Number of MINRES steps vs iteration nr')
3423 CALL hmpdef( 5,0.0,0.0,
'Number of degrees of freedom')
3424 CALL hmpdef(11,0.0,0.0,
'Number of local parameters')
3425 CALL hmpdef(23,0.0,0.0,
'SQRT of diagonal elements without presigma')
3426 CALL hmpdef(24,0.0,0.0,
'Log10 of off-diagonal elements')
3427 CALL hmpdef(25,0.0,0.0,
'Relative individual pre-sigma')
3428 CALL hmpdef(26,0.0,0.0,
'Relative global pre-sigma')
3433 'Normalized residuals of single (global) measurement')
3435 'Normalized residuals of single (local) measurement')
3437 'Pulls of single (global) measurement')
3439 'Pulls of single (local) measurement')
3440 CALL hmpdef(4,0.0,0.0,
'Chi^2/Ndf after local fit')
3441 CALL gmpdef(4,5,
'location, dispersion (res.) vs record nr')
3442 CALL gmpdef(5,5,
'location, dispersion (pull) vs record nr')
3456 IF(
nloopn == 2)
CALL hmpdef(6,0.0,0.0,
'Down-weight fraction')
3459 IF(
iterat /= lastit)
THEN
3467 ELSE IF(
iterat >= 1)
THEN
3517 WRITE(*,111) nparl,ncrit,usei,used,peaki,peakd
3518111
FORMAT(
' Write cache usage (#flush,#overrun,<levels>,', &
3519 'peak(levels))'/2i7,
',',4(f6.1,
'%'))
3527 ELSE IF (
mbandw > 0)
THEN
3559 print *,
" ... warning ..."
3560 print *,
" global parameters with too few (< MREQENA) accepted entries: ", nlow
3564 IF(
icalcm == 1 .AND. nzero > 0)
THEN
3566 WRITE(*,*)
'Warning: the rank defect of the symmetric',
nfgb, &
3567 '-by-',
nfgb,
' matrix is ',
ndefec,
' (should be zero).'
3568 WRITE(lun,*)
'Warning: the rank defect of the symmetric',
nfgb, &
3569 '-by-',
nfgb,
' matrix is ',
ndefec,
' (should be zero).'
3572 WRITE(*,*)
' --> enforcing SUBITO mode'
3573 WRITE(lun,*)
' --> enforcing SUBITO mode'
3583 elmt=real(matij(i,i),mps)
3584 IF(elmt > 0.0)
CALL hmpent(23,1.0/sqrt(elmt))
3618 CALL addsums(1, adder, 0, 1.0_mpl)
3632 IF(sgm > 0.0) weight=1.0/sgm**2
3648 IF(itgbij /= 0)
THEN
3663 IF(
icalcm == 1.AND.ivgbij > 0)
THEN
3670 IF(ivgbij > 0.AND.ivgbik > 0)
THEN
3671 CALL mupdat(ivgbij,ivgbik,weight*factrj*factrk)
3677 adder=weight*dsum**2
3678 CALL addsums(1, adder, 1, 1.0_mpl)
3696 ratae=max(-99.9,ratae)
3705 WRITE(*,*)
'Data rejected in initial loop:'
3707 nrejec(0),
' (rank deficit/NaN) ',
nrejec(1),
' (Ndf=0) ', &
3756 WRITE(lun,*)
' === local fits have bordered band matrix structure ==='
3757 IF (
nbndr(1) > 0 )
WRITE(lun,101)
' NBNDR',
nbndr(1),
'number of records (upper/left border)'
3758 IF (
nbndr(2) > 0 )
WRITE(lun,101)
' NBNDR',
nbndr(2),
'number of records (lower/right border)'
3759 WRITE(lun,101)
' NBDRX',
nbdrx,
'max border size'
3760 WRITE(lun,101)
' NBNDX',
nbndx,
'max band width'
3769101
FORMAT(1x,a8,
' =',i14,
' = ',a)
3784 INTEGER(mpi),
INTENT(IN) :: lunp
3789101
FORMAT(
' it fc',
' fcn_value dfcn_exp slpr costh iit st', &
3790 ' ls step cutf',1x,
'rejects hhmmss FMS')
3791102
FORMAT(
' -- --',
' ----------- -------- ---- ----- --- --', &
3792 ' -- ----- ----',1x,
'------- ------ ---')
3808 INTEGER(mpi) :: nrej
3809 INTEGER(mpi) :: nsecnd
3813 REAL(mps) :: slopes(3)
3814 REAL(mps) :: steps(3)
3815 REAL,
DIMENSION(2) :: ta
3818 INTEGER(mpi),
INTENT(IN) :: lunp
3820 CHARACTER (LEN=4):: ccalcm(4)
3821 DATA ccalcm /
' end',
' S',
' F ',
' FMS' /
3825 IF(nrej > 9999999) nrej=9999999
3829 nsecnd=nint(secnd,mpi)
3838 CALL ptlopt(nfa,ma,slopes,steps)
3839 ratae=max(-99.9,min(99.9,slopes(2)/slopes(1)))
3845103
FORMAT(i3,i3,e12.5,38x,f5.1, 1x,i7, i3,i2.2,i2.2,a4)
3846104
FORMAT(i3,i3,e12.5,1x,e8.2,f6.3,f6.3,i5,2i3,f6.3,f5.1, &
3847 1x,i7, i3,i2.2,i2.2,a4)
3848105
FORMAT(i3,i3,e12.5,1x,e8.2,12x,i5,i3,9x,f5.1, &
3849 1x,i7, i3,i2.2,i2.2,a4)
3862 INTEGER(mpi) :: minut
3864 INTEGER(mpi) :: nhour
3865 INTEGER(mpi) :: nrej
3866 INTEGER(mpi) :: nsecnd
3870 REAL(mps) :: slopes(3)
3871 REAL(mps) :: steps(3)
3872 REAL,
DIMENSION(2) :: ta
3875 INTEGER(mpi),
INTENT(IN) :: lunp
3876 CHARACTER (LEN=4):: ccalcm(4)
3877 DATA ccalcm /
' end',
' S',
' F ',
' FMS' /
3881 IF(nrej > 9999999) nrej=9999999
3885 nsecnd=nint(secnd,mpi)
3889 CALL ptlopt(nfa,ma,slopes,steps)
3890 ratae=abs(slopes(2)/slopes(1))
3895104
FORMAT(3x,i3,e12.5,9x, 35x, i7, i3,i2.2,i2.2,a4)
3896105
FORMAT(3x,i3,e12.5,9x, f6.3,14x,i3,f6.3,6x, i7, i3,i2.2,i2.2,a4)
3910 INTEGER(mpi) :: nsecnd
3913 REAL,
DIMENSION(2) :: ta
3916 INTEGER(mpi),
INTENT(IN) :: lunp
3917 CHARACTER (LEN=4):: ccalcm(4)
3918 DATA ccalcm /
' end',
' S',
' F ',
' FMS' /
3923 nsecnd=nint(secnd,mpi)
3925 WRITE(lunp,106) nhour,minut,nsecnd,ccalcm(
lcalcm)
3926106
FORMAT(69x,i3,i2.2,i2.2,a4)
3936 INTEGER(mpi) :: lunit
3938 WRITE(lunit,102)
'Explanation of iteration table'
3939 WRITE(lunit,102)
'=============================='
3940 WRITE(lunit,101)
'it', &
3941 'iteration number. Global parameters are improved for it > 0.'
3942 WRITE(lunit,102)
'First function evaluation is called iteraton 0.'
3943 WRITE(lunit,101)
'fc',
'number of function evaluations.'
3944 WRITE(lunit,101)
'fcn_value',
'value of 2 x Likelihood function (LF).'
3945 WRITE(lunit,102)
'The final value is the chi^2 value of the fit and should'
3946 WRITE(lunit,102)
'be about equal to the NDF (see below).'
3947 WRITE(lunit,101)
'dfcn_exp', &
3948 'expected reduction of the value of the Likelihood function (LF)'
3949 WRITE(lunit,101)
'slpr',
'ratio of the actual slope to inital slope.'
3950 WRITE(lunit,101)
'costh', &
3951 'cosine of angle between search direction and -gradient'
3953 WRITE(lunit,101)
'iit', &
3954 'number of internal iterations in MINRES algorithm'
3955 WRITE(lunit,101)
'st',
'stop code of MINRES algorithm'
3956 WRITE(lunit,102)
'< 0: rhs is very special, with beta2 = 0'
3957 WRITE(lunit,102)
'= 0: rhs b = 0, i.e. the exact solution is x = 0'
3958 WRITE(lunit,102)
'= 1 requested accuracy achieved, as determined by rtol'
3959 WRITE(lunit,102)
'= 2 reasonable accuracy achieved, given eps'
3960 WRITE(lunit,102)
'= 3 x has converged to an eigenvector'
3961 WRITE(lunit,102)
'= 4 matrix ill-conditioned (Acond has exceeded 0.1/eps)'
3962 WRITE(lunit,102)
'= 5 the iteration limit was reached'
3963 WRITE(lunit,102)
'= 6 Matrix x vector does not define a symmetric matrix'
3964 WRITE(lunit,102)
'= 7 Preconditioner does not define a symmetric matrix'
3965 ELSEIF (
metsol == 5)
THEN
3966 WRITE(lunit,101)
'iit', &
3967 'number of internal iterations in MINRES-QLP algorithm'
3968 WRITE(lunit,101)
'st',
'stop code of MINRES-QLP algorithm'
3969 WRITE(lunit,102)
'= 1: beta_{k+1} < eps, iteration k is the final Lanczos step.'
3970 WRITE(lunit,102)
'= 2: beta2 = 0. If M = I, b and x are eigenvectors of A.'
3971 WRITE(lunit,102)
'= 3: beta1 = 0. The exact solution is x = 0.'
3972 WRITE(lunit,102)
'= 4: A solution to (poss. singular) Ax = b found, given rtol.'
3973 WRITE(lunit,102)
'= 5: A solution to (poss. singular) Ax = b found, given eps.'
3974 WRITE(lunit,102)
'= 6: Pseudoinverse solution for singular LS problem, given rtol.'
3975 WRITE(lunit,102)
'= 7: Pseudoinverse solution for singular LS problem, given eps.'
3976 WRITE(lunit,102)
'= 8: The iteration limit was reached.'
3977 WRITE(lunit,102)
'= 9: The operator defined by Aprod appears to be unsymmetric.'
3978 WRITE(lunit,102)
'=10: The operator defined by Msolve appears to be unsymmetric.'
3979 WRITE(lunit,102)
'=11: The operator defined by Msolve appears to be indefinite.'
3980 WRITE(lunit,102)
'=12: xnorm has exceeded maxxnorm or will exceed it next iteration.'
3981 WRITE(lunit,102)
'=13: Acond has exceeded Acondlim or 0.1/eps.'
3982 WRITE(lunit,102)
'=14: Least-squares problem but no converged solution yet.'
3983 WRITE(lunit,102)
'=15: A null vector obtained, given rtol.'
3985 WRITE(lunit,101)
'ls',
'line search info'
3986 WRITE(lunit,102)
'< 0 recalculate function'
3987 WRITE(lunit,102)
'= 0: N or STP lt 0 or step not descending'
3988 WRITE(lunit,102)
'= 1: Linesearch convergence conditions reached'
3989 WRITE(lunit,102)
'= 2: interval of uncertainty at lower limit'
3990 WRITE(lunit,102)
'= 3: max nr of line search calls reached'
3991 WRITE(lunit,102)
'= 4: step at the lower bound'
3992 WRITE(lunit,102)
'= 5: step at the upper bound'
3993 WRITE(lunit,102)
'= 6: rounding error limitation'
3994 WRITE(lunit,101)
'step', &
3995 'the factor for the Newton step during the line search. Usually'
3997 'a value of 1 gives a sufficient reduction of the LF. Oherwise'
3998 WRITE(lunit,102)
'other step values are tried.'
3999 WRITE(lunit,101)
'cutf', &
4000 'cut factor. Local fits are rejected, if their chi^2 value'
4002 'is larger than the 3-sigma chi^2 value times the cut factor.'
4003 WRITE(lunit,102)
'A cut factor of 1 is used finally, but initially a larger'
4004 WRITE(lunit,102)
'factor may be used. A value of 0.0 means no cut.'
4005 WRITE(lunit,101)
'rejects',
'total number of rejected local fits.'
4006 WRITE(lunit,101)
'hmmsec',
'the time in hours (h), minutes (mm) and seconds.'
4007 WRITE(lunit,101)
'FMS',
'calculation of Function value, Matrix, Solution.'
4010101
FORMAT(a9,
' = ',a)
4027 INTEGER(mpi),
INTENT(IN) :: i
4028 INTEGER(mpi),
INTENT(IN) :: j
4029 REAL(mpd),
INTENT(IN) :: add
4031 INTEGER(mpl):: ijadd
4032 INTEGER(mpl):: ijcsr3
4037 IF(i <= 0.OR.j <= 0.OR. add == 0.0_mpd)
RETURN
4058 ELSE IF(
matsto == 2)
THEN
4086 CALL peend(23,
'Aborted, bad matrix index')
4087 stop
'mupdat: bad index'
4112 INTEGER(mpi),
INTENT(IN) :: i
4113 INTEGER(mpi),
INTENT(IN) :: j1
4114 INTEGER(mpi),
INTENT(IN) :: j2
4115 INTEGER(mpi),
INTENT(IN) :: il
4116 INTEGER(mpi),
INTENT(IN) :: jl
4117 INTEGER(mpi),
INTENT(IN) :: n
4118 REAL(mpd),
INTENT(IN) :: sub((n*n+n)/2)
4125 INTEGER(mpi):: iblast
4126 INTEGER(mpi):: iblock
4132 INTEGER(mpi):: jblast
4133 INTEGER(mpi):: jblock
4143 IF(i <= 0.OR.j1 <= 0.OR.j2 > i)
RETURN
4157 print *,
' MGUPDT: ij=0', i,j1,j2,il,jl,ij,lr,iprc,
matsto
4164 ijl=(k*k-k)/2+jl+ir-ia1
4181 ijl=(k*k-k)/2+jl+ir-ia1
4185 IF (jblock /= jblast.OR.iblock /= iblast)
THEN
4209 CALL ijpgrp(i,j1,ij,lr,iprc)
4211 print *,
' MGUPDT: ij=0', i,j1,j2,il,jl,ij,lr,iprc,
matsto
4275SUBROUTINE loopbf(nrej,numfil,naccf,chi2f,ndff)
4302 INTEGER(mpi) :: ibuf
4303 INTEGER(mpi) :: ichunk
4304 INTEGER(mpl) :: icmn
4305 INTEGER(mpl) :: icost
4307 INTEGER(mpi) :: idiag
4309 INTEGER(mpi) :: iext
4317 INTEGER(mpi) :: imeas
4320 INTEGER(mpi) :: ioffb
4321 INTEGER(mpi) :: ioffc
4322 INTEGER(mpi) :: ioffd
4323 INTEGER(mpi) :: ioffe
4324 INTEGER(mpi) :: ioffi
4325 INTEGER(mpi) :: ioffq
4326 INTEGER(mpi) :: iprc
4327 INTEGER(mpi) :: iprcnx
4328 INTEGER(mpi) :: iprdbg
4329 INTEGER(mpi) :: iproc
4330 INTEGER(mpi) :: irbin
4331 INTEGER(mpi) :: isize
4333 INTEGER(mpi) :: iter
4334 INTEGER(mpi) :: itgbi
4335 INTEGER(mpi) :: ivgbj
4336 INTEGER(mpi) :: ivgbk
4337 INTEGER(mpi) :: ivpgrp
4347 INTEGER(mpi) :: joffd
4348 INTEGER(mpi) :: joffi
4349 INTEGER(mpi) :: jproc
4353 INTEGER(mpi) :: kbdr
4354 INTEGER(mpi) :: kbdrx
4355 INTEGER(mpi) :: kbnd
4358 INTEGER(mpi) :: lvpgrp
4359 INTEGER(mpi) :: mbdr
4360 INTEGER(mpi) :: mbnd
4361 INTEGER(mpi) :: mside
4362 INTEGER(mpi) :: nalc
4363 INTEGER(mpi) :: nalg
4367 INTEGER(mpi) :: ndown
4369 INTEGER(mpi) :: nfred
4370 INTEGER(mpi) :: nfrei
4372 INTEGER(mpi) :: nprdbg
4373 INTEGER(mpi) :: nrank
4376 INTEGER(mpi) :: nter
4377 INTEGER(mpi) :: nweig
4378 INTEGER(mpi) :: ngrp
4379 INTEGER(mpi) :: npar
4381 INTEGER(mpi),
INTENT(IN OUT) :: nrej(0:3)
4382 INTEGER(mpi),
INTENT(IN) :: numfil
4383 INTEGER(mpi),
INTENT(IN OUT) :: naccf(numfil)
4384 REAL(mps),
INTENT(IN OUT) :: chi2f(numfil)
4385 INTEGER(mpi),
INTENT(IN OUT) :: ndff(numfil)
4392 INTEGER(mpi) :: ijprec
4399 CHARACTER (LEN=3):: chast
4400 DATA chuber/1.345_mpd/
4401 DATA cauchy/2.3849_mpd/
4449 IF(
nloopn == 1.AND.mod(nrc,100000_mpl) == 0)
THEN
4450 WRITE(*,*)
'Record',nrc,
' ... still reading'
4451 IF(
monpg1>0)
WRITE(
lunlog,*)
'Record',nrc,
' ... still reading'
4461 IF(nrc ==
nrecpr) lprnt=.true.
4462 IF(nrc ==
nrecp2) lprnt=.true.
4463 IF(nrc ==
nrecer) lprnt=.true.
4468 IF (nprdbg == 1) iprdbg=iproc
4469 IF (iproc /= iprdbg) lprnt=.false.
4474 WRITE(1,*)
'------------------ Loop',
nloopn, &
4475 ': Printout for record',nrc,iproc
4486 CALL isjajb(nst,ist,ja,jb,jsp)
4488 IF(imeas == 0)
WRITE(1,1121)
44931121
FORMAT(/
'Measured value and local derivatives'/ &
4494 ' i measured std_dev index...derivative ...')
44951122
FORMAT(i3,2g12.4,3(i3,g12.4)/(27x,3(i3,g12.4)))
4501 CALL isjajb(nst,ist,ja,jb,jsp)
4503 IF(imeas == 0)
WRITE(1,1123)
45151123
FORMAT(/
'Global derivatives'/ &
4516 ' i label gindex vindex derivative ...')
45171124
FORMAT(i3,2(i9,i7,i7,g12.4)/(3x,2(i9,i7,i7,g12.4)))
45181125
FORMAT(i3,2(i9,i7,i7,g12.4))
4528 WRITE(1,*)
'Data corrections using values of global parameters'
4529 WRITE(1,*)
'=================================================='
4539 CALL isjajb(nst,ist,ja,jb,jsp)
4571101
FORMAT(
' index measvalue corrvalue sigma')
4572102
FORMAT(i6,2x,2g12.4,
' +-',g12.4)
4574 IF(nalc <= 0)
GO TO 90
4576 ngg=(nalg*nalg+nalg)/2
4589 IF (ivpgrp /= lvpgrp)
THEN
4597 IF (npar /= nalg)
THEN
4598 print *,
' mismatch of number of global parameters ', nrc, nalg, npar, ngrp
4608 CALL peend(35,
'Aborted, mismatch of number of global parameters')
4609 stop
' mismatch of number of global parameters '
4636 DO WHILE(iter < nter)
4641 WRITE(1,*)
'Outlier-suppression iteration',iter,
' of',nter
4642 WRITE(1,*)
'=========================================='
4652 DO i=1,(nalc*nalc+nalc)/2
4662 wght =1.0_mpd/rerr**2
4667 IF(abs(resid) > chuber*rerr)
THEN
4668 wght=wght*chuber*rerr/abs(resid)
4672 wght=wght/(1.0+(resid/rerr/cauchy)**2)
4676 IF(lprnt.AND.iter /= 1.AND.nter /= 1)
THEN
4678 IF(abs(resid) > chuber*rerr) chast=
'* '
4679 IF(abs(resid) > 3.0*rerr) chast=
'** '
4680 IF(abs(resid) > 6.0*rerr) chast=
'***'
4681 IF(imeas == 0)
WRITE(1,*)
'Second loop: accumulate'
4682 IF(imeas == 0)
WRITE(1,103)
4687 WRITE(1,104) imeas,rmeas,resid,rerr,r1,chast,r2
4689103
FORMAT(
' index corrvalue residuum sigma', &
4691104
FORMAT(i6,2x,2g12.4,
' +-',g12.4,f7.2,1x,a3,f8.2)
4706 im=min(nalc+1-ij,nalc+1-ik)
4713 IF (iter == 1.AND.nalc > 5.AND.
lfitbb > 0)
THEN
4716 icmn=int(nalc,mpl)**3
4722 icost=6*int(nalc-kbdr,mpl)*int(kbnd+kbdr+1,mpl)**2+2*int(kbdr,mpl)**3
4723 IF (icost < icmn)
THEN
4735 kbdr=max(kbdr,
ibandh(k+nalc))
4736 icost=6*int(nalc-kbdr,mpl)*int(kbnd+kbdr+1,mpl)**2+2*int(kbdr,mpl)**3
4737 IF (icost < icmn)
THEN
4761 IF (
icalcm == 1.OR.lprnt) inv=2
4762 IF (mside == 1)
THEN
4777 IF ((.NOT.(
blvec(k) <= 0.0_mpd)).AND. (.NOT.(
blvec(k) > 0.0_mpd))) nan=nan+1
4782 WRITE(1,*)
'Parameter determination:',nalc,
' parameters,',
' rank=',nrank
4783 WRITE(1,*)
'-----------------------'
4784 IF(ndown /= 0)
WRITE(1,*)
' ',ndown,
' data down-weighted'
4800 wght =1.0_mpd/rerr**2
4824 IF (0.999999_mpd/wght > dvar)
THEN
4825 pull=rmeas/sqrt(1.0_mpd/wght-dvar)
4828 CALL hmpent(13,real(pull,mps))
4829 CALL gmpms(5,rec,real(pull,mps))
4831 CALL hmpent(14,real(pull,mps))
4850 IF(iter == 1.AND.jb < ist.AND.lhist) &
4851 CALL gmpms(4,rec,real(rmeas/rerr,mps))
4853 dchi2=wght*rmeas*rmeas
4858 IF(abs(resid) > chuber*rerr)
THEN
4859 wght=wght*chuber*rerr/abs(resid)
4860 dchi2=2.0*chuber*(abs(resid)/rerr-0.5*chuber)
4863 wght=wght/(1.0_mpd+(resid/rerr/cauchy)**2)
4864 dchi2=log(1.0_mpd+(resid/rerr/cauchy)**2)*cauchy**2
4874 IF(abs(resid) > chuber*rerr) chast=
'* '
4875 IF(abs(resid) > 3.0*rerr) chast=
'** '
4876 IF(abs(resid) > 6.0*rerr) chast=
'***'
4877 IF(imeas == 0)
WRITE(1,*)
'Third loop: single residuals'
4878 IF(imeas == 0)
WRITE(1,105)
4882 IF(resid < 0.0) r1=-r1
4883 IF(resid < 0.0) r2=-r2
4886105
FORMAT(
' index corrvalue residuum sigma', &
4888106
FORMAT(i6,2x,2g12.4,
' +-',g12.4,f7.2,1x,a3,f8.2)
4890 IF(iter == nter)
THEN
4892 resmax=max(resmax,abs(rmeas)/rerr)
4895 IF(iter == 1.AND.lhist)
THEN
4897 CALL hmpent( 3,real(rmeas/rerr,mps))
4899 CALL hmpent(12,real(rmeas/rerr,mps))
4906 resing=(real(nweig,mps)-real(suwt,mps))/real(nweig,mps)
4908 IF(iter == 1)
CALL hmpent( 5,real(ndf,mps))
4909 IF(iter == 1)
CALL hmpent(11,real(nalc,mps))
4914 WRITE(1,*)
'Chi^2=',summ,
' at',ndf,
' degrees of freedom: ', &
4915 '3-sigma limit is',chindl(3,ndf)*real(ndf,mps)
4916 WRITE(1,*) suwt,
' is sum of factors, compared to',nweig, &
4917 ' Downweight fraction:',resing
4919 IF(nrank /= nalc.OR.nan > 0)
THEN
4923 WRITE(1,*)
' rank deficit/NaN ', nalc, nrank, nan
4924 WRITE(1,*)
' ---> rejected!'
4931 WRITE(1,*)
' ---> rejected!'
4936 chndf=real(summ/real(ndf,mpd),mps)
4938 IF(iter == 1.AND.lhist)
CALL hmpent(4,chndf)
4951 chichi=chindl(3,ndf)*real(ndf,mps)
4955 IF(summ >
chhuge*chichi)
THEN
4958 WRITE(1,*)
' ---> rejected!'
4966 IF(summ > chlimt)
THEN
4968 WRITE(1,*)
' ---> rejected!'
4972 CALL addsums(iproc+1, dchi2, ndf, dw1)
4982 CALL addsums(iproc+1, dchi2, ndf, dw1)
4986 WRITE(1,*)
' ---> rejected!'
4999 naccf(kfl)=naccf(kfl)+1
5000 ndff(kfl) =ndff(kfl) +ndf
5001 chi2f(kfl)=chi2f(kfl)+chndf
5014 wght =1.0_mpd/rerr**2
5015 dchi2=wght*rmeas*rmeas
5019 IF(resid > chuber*rerr)
THEN
5020 wght=wght*chuber*rerr/resid
5021 dchi2=2.0*chuber*(resid/rerr-0.5*chuber)
5070 CALL addsums(iproc+1, summ, ndf, dw1)
5126 IF (nfred < 0.OR.nfrei < 0)
THEN
5153 iprcnx=ijprec(i,jnx)
5155 IF (.NOT.((jnx == j+1).AND.(iprc == iprcnx)))
THEN
5169 joffd=joffd+(il*il-il)/2
5181 WRITE(1,*)
'------------------ End of printout for record',nrc
5238 iprcnx=ijprec(i,jnx)
5240 IF (.NOT.((jnx == j+1).AND.(iprc == iprcnx)))
THEN
5255 joffd=joffd+(il*il-il)/2
5308 INTEGER(mpi) :: icom
5309 INTEGER(mpl) :: icount
5313 INTEGER(mpi) :: imin
5314 INTEGER(mpi) :: iprlim
5315 INTEGER(mpi) :: isub
5316 INTEGER(mpi) :: itgbi
5317 INTEGER(mpi) :: itgbl
5318 INTEGER(mpi) :: ivgbi
5320 INTEGER(mpi) :: label
5328 INTEGER(mpi) :: labele(3)
5329 REAL(mps):: compnt(3)
5334 CALL mvopen(lup,
'millepede.res')
5337 WRITE(*,*)
' Result of fit for global parameters'
5338 WRITE(*,*)
' ==================================='
5343 WRITE(lup,*)
'Parameter ! first 3 elements per line are', &
5344 ' significant (if used as input)'
5362 err=sqrt(abs(real(gmati,mps)))
5363 IF(gmati < 0.0_mpd) err=-err
5366 IF(gmati*diag > 0.0_mpd)
THEN
5367 gcor2=1.0_mpd-1.0_mpd/(gmati*diag)
5368 IF(gcor2 >= 0.0_mpd.AND.gcor2 <= 1.0_mpd) gcor=real(sqrt(gcor2),mps)
5373 IF(lowstat) icount=-(icount+1)
5374 IF(itgbi <= iprlim)
THEN
5388 ELSE IF(itgbi == iprlim+1)
THEN
5389 WRITE(* ,*)
'... (further printout suppressed, but see log file)'
5403 ELSE IF (
igcorr /= 0)
THEN
5421 CALL mvopen(lup,
'millepede.eve')
5431 DO isub=0,min(15,imin-1)
5453 WRITE(lup,103) (labele(ie),compnt(ie),ie=1,iev)
5457 IF(iev /= 0)
WRITE(lup,103) (labele(ie),compnt(ie),ie=1,iev)
5464101
FORMAT(1x,
' label parameter presigma differ', &
5465 ' error'/ 1x,
'-----------',4x,4(
'-------------'))
5466102
FORMAT(i10,2x,4g14.5,f8.3)
5467103
FORMAT(3(i11,f11.7,2x))
5468110
FORMAT(i10,2x,2g14.5,28x,i12)
5469111
FORMAT(i10,2x,3g14.5,14x,i12)
5470112
FORMAT(i10,2x,4g14.5,i12)
5492 INTEGER(mpi) :: icom
5493 INTEGER(mpl) :: icount
5494 INTEGER(mpi) :: ifrst
5495 INTEGER(mpi) :: ilast
5496 INTEGER(mpi) :: inext
5497 INTEGER(mpi) :: itgbi
5498 INTEGER(mpi) :: itgbl
5499 INTEGER(mpi) :: itpgrp
5500 INTEGER(mpi) :: ivgbi
5502 INTEGER(mpi) :: icgrp
5503 INTEGER(mpi) :: ipgrp
5505 INTEGER(mpi) :: jpgrp
5507 INTEGER(mpi) :: label1
5508 INTEGER(mpi) :: label2
5509 INTEGER(mpi) :: ncon
5510 INTEGER(mpi) :: npair
5511 INTEGER(mpi) :: nstep
5514 INTEGER(mpl):: length
5516 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: vecPairedParGroups
5519 SUBROUTINE ggbmap(ipgrp,npair,npgrp)
5521 INTEGER(mpi),
INTENT(IN) :: ipgrp
5522 INTEGER(mpi),
INTENT(OUT) :: npair
5523 INTEGER(mpi),
DIMENSION(:),
INTENT(OUT) :: npgrp
5531 CALL mvopen(lup,
'millepede.res')
5532 WRITE(lup,*)
'*** Results of checking input only, no solution performed ***'
5533 WRITE(lup,*)
'! === global parameters ==='
5534 WRITE(lup,*)
'! fixed-1: by pre-sigma, -2: by entries cut, -3: by iterated entries cut'
5535 WRITE(lup,*)
'! Label Value Pre-sigma Entries Cons. group Status '
5549 IF (ivgbi <= 0)
THEN
5551 IF (ivgbi == -4)
THEN
5552 WRITE(lup,116) c1,itgbl,par,presig,icount,icgrp
5554 WRITE(lup,110) c1,itgbl,par,presig,icount,icgrp,ivgbi
5558 WRITE(lup,111) c1,itgbl,par,presig,icount,icgrp
5564 WRITE(lup,*)
'!.Appearance statistics '
5565 WRITE(lup,*)
'!. Label First file and record Last file and record #files #paired-par'
5568 IF (itpgrp > 0)
THEN
5576 WRITE(lup,*)
'* === constraint groups ==='
5578 WRITE(lup,*)
'* Group #Cons. Entries First label Last label'
5580 WRITE(lup,*)
'* Group #Cons. Entries First label Last label Paired label range'
5582 CALL mpalloc(vecpairedpargroups,length,
'paired global parameter groups (I)')
5594 IF (
icheck > 1 .AND. label1 > 0)
THEN
5598 vecpairedpargroups(npair+1)=0
5602 jpgrp=vecpairedpargroups(j)
5606 IF (ifrst /= 0.AND.inext /= (ilast+nstep))
THEN
5609 WRITE(lup,114) label1, label2
5614 IF (ifrst == 0) ifrst=inext
5621 IF (jpgrp == vecpairedpargroups(j+1)-1)
THEN
5626 IF (ifrst /= 0)
THEN
5629 WRITE(lup,114) label1, label2
5635 WRITE(lup,*)
'*.Appearance statistics '
5636 WRITE(lup,*)
'*. Group First file and record Last file and record #files'
5646110
FORMAT(
' !',a1,i10,2x,2g14.5,2i12,
' fixed',i2)
5647111
FORMAT(
' !',a1,i10,2x,2g14.5,2i12,
' variable')
5648112
FORMAT(
' !.',i10,6i11)
5649113
FORMAT(
' * ',i6,i8,3i12)
5650114
FORMAT(
' *:',48x,i12,
' ..',i12)
5651115
FORMAT(
' *.',i10,5i11)
5652116
FORMAT(
' !',a1,i10,2x,2g14.5,2i12,
' redundant')
5682 INTEGER(mpi) :: iproc
5692 INTEGER(mpi),
INTENT(IN) :: n
5693 INTEGER(mpl),
INTENT(IN) :: l
5694 REAL(mpd),
INTENT(IN) :: x(n)
5695 INTEGER(mpi),
INTENT(IN) :: is
5696 INTEGER(mpi),
INTENT(IN) :: ie
5697 REAL(mpd),
INTENT(OUT) :: b(n)
5702 INTEGER(mpl) :: indij
5703 INTEGER(mpl) :: indid
5705 INTEGER(mpi) :: ichunk
5741 CALL peend(24,
'Aborted, vector/matrix size mismatch')
5742 stop
'AVPRDS: mismatched vector and matrix'
5762 IF (ia2 <= ib2) b(ia2:ib2)=b(ia2:ib2)+
globalmatd(ia2:ib2)*x(ia2:ib2)
5773 IF (i <= ie.AND.i >= is)
THEN
5781 IF (j <= ie.AND.j >= is)
THEN
5797 IF (ja2 <= jb2)
THEN
5800 b(i)=b(i)+dot_product(
globalmatd(indij+lj+ja2-ja:indij+lj+jb2-ja),x(ja2:jb2))
5804 IF (
mextnd == 0.AND.ia2 <= ib2)
THEN
5807 b(j)=b(j)+dot_product(
globalmatd(indij+lj+jn*(ia2-ia):indij+lj+jn*(ib2-ia):jn),x(ia2:ib2))
5826 IF (i <= ie.AND.i >= is)
THEN
5834 IF (j <= ie.AND.j >= is)
THEN
5850 IF (ja2 <= jb2)
THEN
5853 b(i)=b(i)+dot_product(real(
globalmatf(indij+lj+ja2-ja:indij+lj+jb2-ja),mpd),x(ja2:jb2))
5857 IF (
mextnd == 0.AND.ia2 <= ib2)
THEN
5860 b(j)=b(j)+dot_product(real(
globalmatf(indij+lj+jn*(ia2-ia):indij+lj+jn*(ib2-ia):jn),mpd),x(ia2:ib2))
5894 INTEGER(mpi) :: iproc
5902 INTEGER(mpi),
INTENT(IN) :: n
5903 INTEGER(mpl),
INTENT(IN) :: l
5904 REAL(mpd),
INTENT(IN) :: x(n)
5905 REAL(mpd),
INTENT(OUT) :: b(n)
5910 INTEGER(mpl) :: indij
5911 INTEGER(mpl) :: indid
5913 INTEGER(mpi) :: ichunk
5942 CALL peend(24,
'Aborted, vector/matrix size mismatch')
5943 stop
'AVPRD0: mismatched vector and matrix'
5989 b(i)=b(i)+dot_product(
globalmatd(indij+lj:indij+lj+jn-1),x(ja:jb))
5995 b(j)=b(j)+dot_product(
globalmatd(indij+lj:indij+jn*in:jn),x(ia:ib))
6031 b(i)=b(i)+dot_product(real(
globalmatf(indij+lj:indij+lj+jn-1),mpd),x(ja:jb))
6037 b(j)=b(j)+dot_product(real(
globalmatf(indij+lj:indij+jn*in:jn),mpd),x(ia:ib))
6061 INTEGER(mpi) :: ispc
6072 INTEGER(mpl) :: indid
6073 INTEGER(mpl),
DIMENSION(12) :: icount
6080 icount(4)=huge(icount(4))
6081 icount(7)=huge(icount(7))
6082 icount(10)=huge(icount(10))
6092 ir=ipg+(ispc-1)*(
napgrp+1)
6099 icount(1)=icount(1)+in
6100 icount(2)=icount(2)+in*ku
6106 icount(3)=icount(3)+1
6107 icount(4)=min(icount(4),jn)
6108 icount(5)=icount(5)+jn
6109 icount(6)=max(icount(6),jn)
6110 icount(7)=min(icount(7),in)
6111 icount(8)=icount(8)+in
6112 icount(9)=max(icount(9),in)
6113 icount(10)=min(icount(10),in*jn)
6114 icount(11)=icount(11)+in*jn
6115 icount(12)=max(icount(12),in*jn)
6121 WRITE(*,*)
"analysis of sparsity structure"
6122 IF (icount(1) > 0)
THEN
6123 WRITE(*,101)
"rows without compression/blocks ", icount(1)
6124 WRITE(*,101)
" contained elements ", icount(2)
6126 WRITE(*,101)
"number of block matrices ", icount(3)
6127 avg=real(icount(5),mps)/real(icount(3),mps)
6128 WRITE(*,101)
"number of columns (min,mean,max) ", icount(4), avg, icount(6)
6129 avg=real(icount(8),mps)/real(icount(3),mps)
6130 WRITE(*,101)
"number of rows (min,mean,max) ", icount(7), avg, icount(9)
6131 avg=real(icount(11),mps)/real(icount(3),mps)
6132 WRITE(*,101)
"number of elements (min,mean,max) ", icount(10), avg, icount(12)
6133101
FORMAT(2x,a34,i10,f10.3,i10)
6152 INTEGER(mpi),
INTENT(IN) :: n
6153 REAL(mpd),
INTENT(IN) :: x(n)
6154 REAL(mpd),
INTENT(OUT) :: b(n)
6159 CALL peend(24,
'Aborted, vector/matrix size mismatch')
6160 stop
'AVPROD: mismatched vector and matrix'
6191 INTEGER(mpi) :: ispc
6192 INTEGER(mpi) :: item1
6193 INTEGER(mpi) :: item2
6194 INTEGER(mpi) :: itemc
6195 INTEGER(mpi) :: jtem
6196 INTEGER(mpi) :: jtemn
6199 INTEGER(mpi),
INTENT(IN) :: itema
6200 INTEGER(mpi),
INTENT(IN) :: itemb
6201 INTEGER(mpl),
INTENT(OUT) :: ij
6202 INTEGER(mpi),
INTENT(OUT) :: lr
6203 INTEGER(mpi),
INTENT(OUT) :: iprc
6214 item1=max(itema,itemb)
6215 item2=min(itema,itemb)
6216 IF(item2 <= 0.OR.item1 >
napgrp)
RETURN
6219 outer:
DO ispc=1,
nspc
6230 IF(ku < kl) cycle outer
6235 IF(item2 >= jtem.AND.item2 < jtemn)
THEN
6241 IF(item2 < jtem)
THEN
6243 ELSE IF(item2 >= jtemn)
THEN
6258 IF(ku < kl) cycle outer
6262 IF(itemc == jtem)
EXIT
6263 IF(itemc < jtem)
THEN
6265 ELSE IF(itemc > jtem)
THEN
6293 INTEGER(mpi),
INTENT(IN) :: itema
6294 INTEGER(mpi),
INTENT(IN) :: itemb
6319 INTEGER(mpi) :: item1
6320 INTEGER(mpi) :: item2
6321 INTEGER(mpi) :: ipg1
6322 INTEGER(mpi) :: ipg2
6324 INTEGER(mpi) :: iprc
6326 INTEGER(mpi),
INTENT(IN) :: itema
6327 INTEGER(mpi),
INTENT(IN) :: itemb
6329 INTEGER(mpl) ::
ijadd
6333 item1=max(itema,itemb)
6334 item2=min(itema,itemb)
6336 IF(item2 <= 0.OR.item1 >
nagb)
RETURN
6337 IF(item1 == item2)
THEN
6346 CALL ijpgrp(ipg1,ipg2,ij,lr,iprc)
6368 INTEGER(mpi) :: item1
6369 INTEGER(mpi) :: item2
6370 INTEGER(mpi) :: jtem
6372 INTEGER(mpi),
INTENT(IN) :: itema
6373 INTEGER(mpi),
INTENT(IN) :: itemb
6382 item1=max(itema,itemb)
6383 item2=min(itema,itemb)
6385 IF(item2 <= 0.OR.item1 >
nagb)
RETURN
6393 print *,
' IJCSR3 empty list ', item1, item2, ks, ke
6394 CALL peend(23,
'Aborted, bad matrix index')
6395 stop
'ijcsr3: empty list'
6400 IF(item1 == jtem)
EXIT
6401 IF(item1 < jtem)
THEN
6408 print *,
' IJCSR3 not found ', item1, item2, ks, ke
6409 CALL peend(23,
'Aborted, bad matrix index')
6410 stop
'ijcsr3: not found'
6426 INTEGER(mpi) :: item1
6427 INTEGER(mpi) :: item2
6431 INTEGER(mpl) ::
ijadd
6434 INTEGER(mpi),
INTENT(IN) :: itema
6435 INTEGER(mpi),
INTENT(IN) :: itemb
6440 item1=max(itema,itemb)
6441 item2=min(itema,itemb)
6442 IF(item2 <= 0.OR.item1 >
nagb)
RETURN
6451 ij=
ijadd(item1,item2)
6454 ELSE IF (ij < 0)
THEN
6486 INTEGER(mpi) :: ichunk
6490 INTEGER(mpi) :: ispc
6498 INTEGER(mpl) :: ijadd
6520 ir=ipg+(ispc-1)*(
napgrp+1)
6570 REAL(mps),
INTENT(IN) :: deltat
6571 INTEGER(mpi),
INTENT(OUT) :: minut
6572 INTEGER(mpi),
INTENT(OUT):: nhour
6573 REAL(mps),
INTENT(OUT):: secnd
6574 INTEGER(mpi) :: nsecd
6577 nsecd=nint(deltat,
mpi)
6579 minut=nsecd/60-60*nhour
6580 secnd=deltat-60*(minut+60*nhour)
6616 INTEGER(mpi),
INTENT(IN) :: item
6620 INTEGER(mpl) :: length
6621 INTEGER(mpl),
PARAMETER :: four = 4
6625 IF(item <= 0)
RETURN
6646 IF(j == 0)
EXIT inner
6668 IF (
lvllog > 1)
WRITE(
lunlog,*)
'INONE: array increased to', &
6689 INTEGER(mpi) :: iprime
6690 INTEGER(mpi) :: nused
6691 LOGICAL :: finalUpdate
6692 INTEGER(mpl) :: oldLength
6693 INTEGER(mpl) :: newLength
6694 INTEGER(mpl),
PARAMETER :: four = 4
6695 INTEGER(mpi),
DIMENSION(:,:),
ALLOCATABLE :: tempArr
6696 INTEGER(mpl),
DIMENSION(:),
ALLOCATABLE :: tempVec
6700 IF(finalupdate)
THEN
6709 CALL mpalloc(temparr,four,oldlength,
'INONE: temp array')
6711 CALL mpalloc(tempvec,oldlength,
'INONE: temp vector')
6735 IF(j == 0)
EXIT inner
6736 IF(j == i) cycle outer
6740 IF(.NOT.finalupdate)
RETURN
6745 WRITE(
lunlog,*)
'INONE: array reduced to',newlength,
' words'
6769 IF(j == 0)
EXIT inner
6770 IF(j == i) cycle outer
6787 INTEGER(mpi),
INTENT(IN) :: n
6788 INTEGER(mpi) :: nprime
6789 INTEGER(mpi) :: nsqrt
6794 IF(mod(nprime,2) == 0) nprime=nprime+1
6797 nsqrt=int(sqrt(real(nprime,
mps)),
mpi)
6799 IF(i*(nprime/i) == nprime) cycle outer
6821 INTEGER(mpi) :: idum
6823 INTEGER(mpi) :: indab
6824 INTEGER(mpi) :: itgbi
6825 INTEGER(mpi) :: itgbl
6826 INTEGER(mpi) :: ivgbi
6828 INTEGER(mpi) :: jgrp
6829 INTEGER(mpi) :: lgrp
6831 INTEGER(mpi) :: nc31
6833 INTEGER(mpi) :: nwrd
6834 INTEGER(mpi) :: inone
6839 INTEGER(mpl) :: length
6840 INTEGER(mpl) :: rows
6844 WRITE(
lunlog,*)
'LOOP1: starting'
6867 WRITE(
lunlog,*)
'LOOP1: reading data files'
6880 WRITE(*,*)
'Read all binary data files:'
6882 CALL hmpldf(1,
'Number of words/record in binary file')
6883 CALL hmpdef(8,0.0,60.0,
'not_stored data per record')
6913 WRITE(
lunlog,*)
'LOOP1: reading data files again'
6925 CALL peend(21,
'Aborted, no labels/parameters defined')
6926 stop
'LOOP1: no labels/parameters defined'
6931 ' is total number NTGB of labels/parameters'
6933 CALL hmpldf(2,
'Number of entries per label')
6977 WRITE(
lunlog,*)
'LOOP1:',
npresg,
' is number of pre-sigmas'
6978 WRITE(*,*)
'LOOP1:',
npresg,
' is number of pre-sigmas'
6979 IF(
npresg == 0)
WRITE(*,*)
'Warning: no pre-sigmas defined'
7001 WRITE(
lunlog,*)
'LOOP1:',
nvgb,
' is number NVGB of variable parameters'
7006 WRITE(
lunlog,*)
'LOOP1: counting records, NO iteration of entries cut !'
7012 CALL hmpdef(15,0.0,120.0,
'Number of parameters per group')
7047 ' is total number NTPGRP of label/parameter groups'
7063 WRITE(*,112)
' Default pre-sigma =',
regpre, &
7064 ' (if no individual pre-sigma defined)'
7065 WRITE(*,*)
'Pre-sigma factor is',
regula
7068 WRITE(*,*)
'No regularization will be done'
7070 WRITE(*,*)
'Regularization will be done, using factor',
regula
7074 CALL peend(22,
'Aborted, no variable global parameters')
7075 stop
'... no variable global parameters'
7082 IF(presg > 0.0)
THEN
7084 ELSE IF(presg == 0.0.AND.
regpre > 0.0)
THEN
7085 prewt=1.0/real(
regpre**2,mpd)
7101 WRITE(*,101)
'NTGB',
ntgb,
'total number of parameters'
7102 WRITE(*,101)
'NVGB',
nvgb,
'number of variable parameters'
7106 WRITE(*,101)
'Too many variable parameters for diagonalization, fallback is inversion'
7114 WRITE(*,101)
' NREC',
nrec,
'number of records'
7115 IF (
nrecd > 0)
WRITE(*,101)
' NRECD',
nrec,
'number of records containing doubles'
7116 WRITE(*,101)
' NEQN',
neqn,
'number of equations (measurements)'
7117 WRITE(*,101)
' NEGB',
negb,
'number of equations with global parameters'
7118 WRITE(*,101)
' NDGB',
ndgb,
'number of global derivatives'
7120 WRITE(*,101)
'MREQENF',
mreqenf,
'required number of entries (eqns in binary files)'
7122 WRITE(*,101)
'MREQENF',
mreqenf,
'required number of entries (recs in binary files)'
7125 WRITE(*,101)
'ITEREN',
iteren,
'iterate cut for parameters with less entries'
7126 WRITE(*,101)
'MREQENA',
mreqena,
'required number of entries (from accepted fits)'
7127 IF (
mreqpe > 1)
WRITE(*,101) &
7128 'MREQPE',
mreqpe,
'required number of pair entries'
7129 IF (
msngpe >= 1)
WRITE(*,101) &
7130 'MSNGPE',
msngpe,
'max pair entries single prec. storage'
7131 WRITE(*,101)
'NTGB',
ntgb,
'total number of parameters'
7132 WRITE(*,101)
'NVGB',
nvgb,
'number of variable parameters'
7135 WRITE(*,*)
'Global parameter labels:'
7142 mqi=((mqi-20)/20)*20+1
7150 WRITE(8,101)
' NREC',
nrec,
'number of records'
7151 IF (
nrecd > 0)
WRITE(8,101)
' NRECD',
nrec,
'number of records containing doubles'
7152 WRITE(8,101)
' NEQN',
neqn,
'number of equations (measurements)'
7153 WRITE(8,101)
' NEGB',
negb,
'number of equations with global parameters'
7154 WRITE(8,101)
' NDGB',
ndgb,
'number of global derivatives'
7156 WRITE(8,101)
'MREQENF',
mreqenf,
'required number of entries (eqns in binary files)'
7158 WRITE(8,101)
'MREQENF',
mreqenf,
'required number of entries (recs in binary files)'
7161 WRITE(8,101)
'ITEREN',
iteren,
'iterate cut for parameters with less entries'
7162 WRITE(8,101)
'MREQENA',
mreqena,
'required number of entries (from accepted fits)'
7164 WRITE(
lunlog,*)
'LOOP1: ending'
7168101
FORMAT(1x,a8,
' =',i14,
' = ',a)
7184 INTEGER(mpi) :: ibuf
7186 INTEGER(mpi) :: indab
7192 INTEGER(mpi) :: nc31
7194 INTEGER(mpi) :: nlow
7196 INTEGER(mpi) :: nwrd
7198 INTEGER(mpl) :: length
7199 INTEGER(mpl),
DIMENSION(:),
ALLOCATABLE :: newCounter
7204 WRITE(
lunlog,*)
'LOOP1: iterating'
7206 WRITE(*,*)
'LOOP1: iterating'
7209 CALL mpalloc(newcounter,length,
'new entries counter')
7233 CALL isjajb(nst,ist,ja,jb,jsp)
7234 IF(ja == 0.AND.jb == 0)
EXIT
7240 IF(ij == -2) nlow=nlow+1
7245 newcounter(ij)=newcounter(ij)+1
7274 WRITE(
lunlog,*)
'LOOP1:',
nvgb,
' is number NVGB of variable parameters'
7304 INTEGER(mpi) :: ibuf
7305 INTEGER(mpi) :: icblst
7306 INTEGER(mpi) :: icboff
7307 INTEGER(mpi) :: icgb
7308 INTEGER(mpi) :: icgrp
7309 INTEGER(mpi) :: icount
7310 INTEGER(mpi) :: iext
7311 INTEGER(mpi) :: ihis
7315 INTEGER(mpi) :: ioff
7316 INTEGER(mpi) :: ipoff
7317 INTEGER(mpi) :: iproc
7318 INTEGER(mpi) :: irecmm
7320 INTEGER(mpi) :: itgbi
7321 INTEGER(mpi) :: itgbij
7322 INTEGER(mpi) :: itgbik
7323 INTEGER(mpi) :: ivgbij
7324 INTEGER(mpi) :: ivgbik
7325 INTEGER(mpi) :: ivpgrp
7329 INTEGER(mpi) :: jcgrp
7330 INTEGER(mpi) :: jext
7331 INTEGER(mpi) :: jcgb
7332 INTEGER(mpi) :: jrec
7334 INTEGER(mpi) :: joff
7336 INTEGER(mpi) :: kcgrp
7337 INTEGER(mpi) :: kfile
7339 INTEGER(mpi) :: label
7340 INTEGER(mpi) :: labelf
7341 INTEGER(mpi) :: labell
7342 INTEGER(mpi) :: lvpgrp
7345 INTEGER(mpi) :: maeqnf
7346 INTEGER(mpi) :: nall
7347 INTEGER(mpi) :: naeqna
7348 INTEGER(mpi) :: naeqnf
7349 INTEGER(mpi) :: naeqng
7350 INTEGER(mpi) :: npdblk
7351 INTEGER(mpi) :: nc31
7352 INTEGER(mpi) :: ncachd
7353 INTEGER(mpi) :: ncachi
7354 INTEGER(mpi) :: ncachr
7355 INTEGER(mpi) :: ncon
7358 INTEGER(mpi) :: ndfmax
7359 INTEGER(mpi) :: nfixed
7360 INTEGER(mpi) :: nggd
7361 INTEGER(mpi) :: nggi
7362 INTEGER(mpi) :: nmatmo
7363 INTEGER(mpi) :: noff
7364 INTEGER(mpi) :: npair
7365 INTEGER(mpi) :: npar
7366 INTEGER(mpi) :: nparmx
7368 INTEGER(mpi) :: nrece
7369 INTEGER(mpi) :: nrecf
7370 INTEGER(mpi) :: nrecmm
7372 INTEGER(mpi) :: nwrd
7373 INTEGER(mpi) :: inone
7382 INTEGER(mpl):: nblock
7383 INTEGER(mpl):: nbwrds
7384 INTEGER(mpl):: noff8
7385 INTEGER(mpl):: ndimbi
7386 INTEGER(mpl):: ndimsa(4)
7388 INTEGER(mpl):: nnzero
7389 INTEGER(mpl):: matsiz(2)
7390 INTEGER(mpl):: matwords
7391 INTEGER(mpl):: mbwrds
7392 INTEGER(mpl):: length
7395 INTEGER(mpl),
PARAMETER :: two=2
7396 INTEGER(mpi) :: maxGlobalPar = 0
7397 INTEGER(mpi) :: maxLocalPar = 0
7398 INTEGER(mpi) :: maxEquations = 0
7400 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: vecConsGroupList
7401 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: vecConsGroupIndex
7402 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: vecPairedParGroups
7403 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: vecBlockCounts
7406 SUBROUTINE ndbits(npgrp,ndims,nsparr,ihst)
7408 INTEGER(mpi),
DIMENSION(:),
INTENT(IN) :: npgrp
7409 INTEGER(mpl),
DIMENSION(4),
INTENT(OUT) :: ndims
7410 INTEGER(mpl),
DIMENSION(:,:),
INTENT(OUT) :: nsparr
7411 INTEGER(mpi),
INTENT(IN) :: ihst
7413 SUBROUTINE ckbits(npgrp,ndims)
7415 INTEGER(mpi),
DIMENSION(:),
INTENT(IN) :: npgrp
7416 INTEGER(mpl),
DIMENSION(4),
INTENT(OUT) :: ndims
7418 SUBROUTINE spbits(npgrp,nsparr,nsparc)
7420 INTEGER(mpi),
DIMENSION(:),
INTENT(IN) :: npgrp
7421 INTEGER(mpl),
DIMENSION(:,:),
INTENT(IN) :: nsparr
7422 INTEGER(mpi),
DIMENSION(:),
INTENT(OUT) :: nsparc
7424 SUBROUTINE gpbmap(ngroup,npgrp,npair)
7426 INTEGER(mpi),
INTENT(IN) :: ngroup
7427 INTEGER(mpi),
DIMENSION(:,:),
INTENT(IN) :: npgrp
7428 INTEGER(mpi),
DIMENSION(:),
INTENT(OUT) :: npair
7430 SUBROUTINE ggbmap(ipgrp,npair,npgrp)
7432 INTEGER(mpi),
INTENT(IN) :: ipgrp
7433 INTEGER(mpi),
INTENT(OUT) :: npair
7434 INTEGER(mpi),
DIMENSION(:),
INTENT(OUT) :: npgrp
7436 SUBROUTINE pbsbits(npgrp,ibsize,nnzero,nblock,nbkrow)
7438 INTEGER(mpi),
DIMENSION(:),
INTENT(IN) :: npgrp
7439 INTEGER(mpi),
INTENT(IN) :: ibsize
7440 INTEGER(mpl),
INTENT(OUT) :: nnzero
7441 INTEGER(mpl),
INTENT(OUT) :: nblock
7442 INTEGER(mpi),
DIMENSION(:),
INTENT(OUT) :: nbkrow
7444 SUBROUTINE pblbits(npgrp,ibsize,nsparr,nsparc)
7446 INTEGER(mpi),
DIMENSION(:),
INTENT(IN) :: npgrp
7447 INTEGER(mpi),
INTENT(IN) :: ibsize
7448 INTEGER(mpl),
DIMENSION(:),
INTENT(IN) :: nsparr
7449 INTEGER(mpl),
DIMENSION(:),
INTENT(OUT) :: nsparc
7451 SUBROUTINE prbits(npgrp,nsparr)
7453 INTEGER(mpi),
DIMENSION(:),
INTENT(IN) :: npgrp
7454 INTEGER(mpl),
DIMENSION(:),
INTENT(OUT) :: nsparr
7456 SUBROUTINE pcbits(npgrp,nsparr,nsparc)
7458 INTEGER(mpi),
DIMENSION(:),
INTENT(IN) :: npgrp
7459 INTEGER(mpl),
DIMENSION(:),
INTENT(IN) :: nsparr
7460 INTEGER(mpl),
DIMENSION(:),
INTENT(OUT) :: nsparc
7470 WRITE(
lunlog,*)
'LOOP2: starting'
7489 WRITE(
lunlog,*)
' Elimination for constraints enforced for solution by decomposition!'
7494 WRITE(
lunlog,*)
' Lagrange multipliers enforced for solution by sparsePARDISO!'
7499 WRITE(
lunlog,*)
' Elimination for constraints with mpqldec enforced (LAPACK only for unpacked storage)!'
7516 noff8=int(
nagb,mpl)*int(
nagb-1,mpl)/2
7552 print *,
' Variable parameter groups ',
nvpgrp
7600 CALL mpalloc(vecconsgrouplist,length,
'constraint group list')
7601 CALL mpalloc(vecconsgroupindex,length,
'constraint group index')
7611 WRITE(
lunlog,*)
'LOOP2: start event reading'
7651 WRITE(*,*)
'Record number ',
nrec,
' from file ',kfile
7652 IF (wgh /= 1.0)
WRITE(*,*)
' weight ',wrec
7656 CALL isjajb(nst,ist,ja,jb,jsp)
7660 IF(nda ==
mdebg2+1)
WRITE(*,*)
'... and more data'
7665 WRITE(*,*)
'Local derivatives:'
7667107
FORMAT(6(i3,g12.4))
7669 WRITE(*,*)
'Global derivatives:'
7672108
FORMAT(3i11,g12.4)
7675 WRITE(*,*)
'total_par_label __label__ var_par_index derivative'
7690 CALL isjajb(nst,ist,ja,jb,jsp)
7691 IF(ja == 0.AND.jb == 0)
EXIT
7727 IF (kcgrp == jcgrp) cycle
7738 IF (vecconsgroupindex(k) == 0)
THEN
7741 vecconsgrouplist(icgrp)=k
7756 IF (vecconsgroupindex(k) < icount)
THEN
7758 vecconsgroupindex(k)=icount
7776 IF (nfixed > 0) naeqnf=naeqnf+1
7779 IF(ja /= 0.AND.jb /= 0)
THEN
7788 IF (naeqnf > maeqnf) nrecf=nrecf+1
7792 maxglobalpar=max(
nagbn,maxglobalpar)
7793 maxlocalpar=max(
nalcn,maxlocalpar)
7794 maxequations=max(
naeqn,maxequations)
7797 dstat(1)=dstat(1)+real((nwrd+2)*2,mpd)
7798 dstat(2)=dstat(2)+real(
nagbn+2,mpd)
7803 vecconsgroupindex(vecconsgrouplist(k))=0
7808 IF (
nagbn == 0)
THEN
7827 IF (ivpgrp /= lvpgrp)
THEN
7882 (irecmm >= nrecmm.OR.irecmm ==
mxrec))
THEN
7883 IF (nmatmo == 0)
THEN
7885 WRITE(*,*)
'Monitoring of sparse matrix construction'
7886 WRITE(*,*)
' records ........ off-diagonal elements ', &
7887 '....... compression memory'
7888 WRITE(*,*)
' non-zero used(double) used', &
7893 gbc=1.0e-9*real((mpi*ndimsa(2)+mpd*ndimsa(3)+mps*ndimsa(4))/mpi*(bit_size(1_mpi)/8),mps)
7894 gbu=1.0e-9*real(((mpi+mpd)*(ndimsa(3)+ndimsa(4)))/mpi*(bit_size(1_mpi)/8),mps)
7896 WRITE(*,1177) irecmm,ndimsa(1),ndimsa(3),ndimsa(4),cpr,gbc
78971177
FORMAT(i9,3i13,f10.2,f11.6)
7898 DO WHILE(irecmm >= nrecmm)
7917 WRITE(
lunlog,*)
'LOOP2: event reading ended - end of data'
7919 dstat(k)=dstat(k)/real(
nrec,mpd)
7933 CALL mpalloc(vecpairedpargroups,length,
'paired global parameter groups (I)')
7935 print *,
' Total parameter groups pairs',
ntpgrp
7938 CALL ggbmap(i,npair,vecpairedpargroups)
8006 IF(ivgbij > 0.AND.ivgbik > 0)
THEN
8008 IF (
mprint > 1)
WRITE(*,*)
'add index pair ',ivgbij,ivgbik
8085 nparmx=max(nparmx,int(rows,mpi))
8107 WRITE(*,*)
' Parameter block', i, ib-ia, jb-ja, labelf, labell
8111 WRITE(
lunlog,*)
'Detected',
npblck,
'(disjoint) parameter blocks, max size ', nparmx
8113 WRITE(*,*)
'Detected',
npblck,
'(disjoint) parameter blocks, max size ', nparmx
8115 WRITE(*,*)
'Using block diagonal storage mode'
8130 IF (
nagb >= 65536)
THEN
8131 noff=int(noff8/1000,mpi)
8141 CALL hmpdef(ihis,0.0,real(
mhispe,mps),
'NDBITS: #off-diagonal elements')
8146 ndgn=ndimsa(3)+ndimsa(4)
8147 matwords=ndimsa(2)+length*4
8162 length=int(npdblk,mpl)
8163 CALL mpalloc(vecblockcounts,length,
'sparse matrix row offsets (CSR3)')
8165 nbwrds=2*int(nblock,mpl)*int(
ipdbsz(i)*
ipdbsz(i)+1,mpl)
8166 IF ((i == 1).OR.(nbwrds < mbwrds))
THEN
8190 IF (
fcache(2) == 0.0)
THEN
8192 fcache(2)=real(dstat(2),mps)
8193 fcache(3)=real(dstat(3),mps)
8214 ncachd=
ncache-ncachr-ncachi
8245 WRITE(lu,101)
'NTGB',
ntgb,
'total number of parameters'
8246 WRITE(lu,102)
'(all parameters, appearing in binary files)'
8247 WRITE(lu,101)
'NVGB',
nvgb,
'number of variable parameters'
8248 WRITE(lu,102)
'(appearing in fit matrix/vectors)'
8249 WRITE(lu,101)
'NAGB',
nagb,
'number of all parameters'
8250 WRITE(lu,102)
'(including Lagrange multiplier or reduced)'
8251 WRITE(lu,101)
'NTPGRP',
ntpgrp,
'total number of parameter groups'
8252 WRITE(lu,101)
'NVPGRP',
nvpgrp,
'number of variable parameter groups'
8253 WRITE(lu,101)
'NFGB',
nfgb,
'number of fit parameters'
8255 WRITE(lu,101)
'MBANDW',
mbandw,
'band width of preconditioner matrix'
8256 WRITE(lu,102)
'(if <0, no preconditioner matrix)'
8258 IF (
nagb >= 65536)
THEN
8259 WRITE(lu,101)
'NOFF/K',noff,
'max number of off-diagonal elements'
8261 WRITE(lu,101)
'NOFF',noff,
'max number of off-diagonal elements'
8264 IF (
nagb >= 65536)
THEN
8265 WRITE(lu,101)
'NDGN/K',ndgn/1000,
'actual number of off-diagonal elements'
8267 WRITE(lu,101)
'NDGN',ndgn,
'actual number of off-diagonal elements'
8270 WRITE(lu,101)
'NCGB',
ncgb,
'number of constraints'
8271 WRITE(lu,101)
'NAGBN',
nagbn,
'max number of global parameters in an event'
8272 WRITE(lu,101)
'NALCN',
nalcn,
'max number of local parameters in an event'
8273 WRITE(lu,101)
'NAEQN',
naeqn,
'max number of equations in an event'
8275 WRITE(lu,101)
'NAEQNA',naeqna,
'number of equations'
8276 WRITE(lu,101)
'NAEQNG',naeqng, &
8277 'number of equations with global parameters'
8278 WRITE(lu,101)
'NAEQNF',naeqnf, &
8279 'number of equations with fixed global parameters'
8280 WRITE(lu,101)
'NRECF',nrecf, &
8281 'number of records with fixed global parameters'
8284 WRITE(lu,101)
'NRECE',nrece, &
8285 'number of records without variable parameters'
8288 WRITE(lu,101)
'NCACHE',
ncache,
'number of words for caching'
8289 WRITE(lu,111) (
fcache(k)*100.0,k=1,3)
8290111
FORMAT(22x,
'cache splitting ',3(f6.1,
' %'))
8295 WRITE(lu,*)
'Solution method and matrix-storage mode:'
8297 WRITE(lu,*)
' METSOL = 1: matrix inversion'
8298 ELSE IF(
metsol == 2)
THEN
8299 WRITE(lu,*)
' METSOL = 2: diagonalization'
8300 ELSE IF(
metsol == 3)
THEN
8301 WRITE(lu,*)
' METSOL = 3: decomposition'
8302 ELSE IF(
metsol == 4)
THEN
8303 WRITE(lu,*)
' METSOL = 4: MINRES (rtol',
mrestl,
')'
8304 ELSE IF(
metsol == 5)
THEN
8305 WRITE(lu,*)
' METSOL = 5: MINRES-QLP (rtol',
mrestl,
')'
8306 ELSE IF(
metsol == 6)
THEN
8307 WRITE(lu,*)
' METSOL = 6: GMRES'
8309 ELSE IF(
metsol == 7)
THEN
8310 WRITE(lu,*)
' METSOL = 7: LAPACK factorization'
8311 ELSE IF(
metsol == 8)
THEN
8312 WRITE(lu,*)
' METSOL = 8: LAPACK factorization'
8314 ELSE IF(
metsol == 9)
THEN
8315 WRITE(lu,*)
' METSOL = 9: Intel oneMKL PARDISO'
8319 WRITE(lu,*)
' with',
mitera,
' iterations'
8321 WRITE(lu,*)
' MATSTO = 0: unpacked symmetric matrix, ',
'n*n elements'
8322 ELSE IF(
matsto == 1)
THEN
8323 WRITE(lu,*)
' MATSTO = 1: full symmetric matrix, ',
'(n*n+n)/2 elements'
8324 ELSE IF(
matsto == 2)
THEN
8325 WRITE(lu,*)
' MATSTO = 2: sparse matrix (custom)'
8326 ELSE IF(
matsto == 3)
THEN
8328 WRITE(lu,*)
' MATSTO = 3: sparse matrix (upper triangle, CSR3)'
8330 WRITE(lu,*)
' MATSTO = 3: sparse matrix (upper triangle, BSR3)'
8331 WRITE(lu,*)
' block size',
matbsz
8335 WRITE(lu,*)
' block diagonal with',
npblck,
' blocks'
8337 IF(
mextnd>0)
WRITE(lu,*)
' with extended storage'
8338 IF(
dflim /= 0.0)
THEN
8339 WRITE(lu,103)
'Convergence assumed, if expected dF <',
dflim
8344 WRITE(lu,*)
'Constraints handled by elimination with LAPACK'
8346 WRITE(lu,*)
'Constraints handled by elimination'
8349 WRITE(lu,*)
'Constraints handled by Lagrange multipliers'
8356 CALL peend(28,
'Aborted, no local parameters')
8357 stop
'LOOP2: stopping due to missing local parameters'
8378105
FORMAT(
' Constants C1, C2 for Wolfe conditions:',g12.4,
', ',g12.4)
8386 matwords=(length+
nagb+1)*2
8393 ELSE IF (
matsto == 2)
THEN
8394 matsiz(1)=ndimsa(3)+
nagb
8409 IF (icblst > icboff)
THEN
8415 nall = npar;
IF (
icelim <= 0) nall=npar+ncon
8419 matsiz(1)=matsiz(1)+k
8421 matsiz(1)=matsiz(1)+nall
8426 matwords=matwords+matsiz(1)*2+matsiz(2)
8432 IF (matwords < 250000)
THEN
8433 WRITE(*,*)
'Size of global matrix: < 1 MB'
8435 WRITE(*,*)
'Size of global matrix:',int(real(matwords,mps)*4.0e-6,mpi),
' MB'
8441 WRITE(
lunlog,*)
' Cut values of Chi^2/Ndf and Chi2,'
8442 WRITE(
lunlog,*)
' corresponding to 2 and 3 standard deviations'
8443 WRITE(
lunlog,*)
' Ndf Chi^2/Ndf(2) Chi^2(2) ', &
8444 ' Chi^2/Ndf(3) Chi^2(3)'
8447 IF(ndf >
naeqn)
EXIT
8450 ELSE IF(ndf < 20)
THEN
8452 ELSE IF(ndf < 100)
THEN
8454 ELSE IF(ndf < 200)
THEN
8461 WRITE(
lunlog,106) ndf,chin2,chin2*real(ndf,mps),chin3, chin3*real(ndf,mps)
8464 WRITE(
lunlog,*)
'LOOP2: ending'
8468 IF (
ncgbe /= 0)
THEN
8471 WRITE(*,199)
'WarningWarningWarningWarningWarningWarningWarningWarningWar'
8472 WRITE(*,199)
'arningWarningWarningWarningWarningWarningWarningWarningWarn'
8473 WRITE(*,199)
'rningWarningWarningWarningWarningWarningWarningWarningWarni'
8474 WRITE(*,199)
'ningWarningWarningWarningWarningWarningWarningWarningWarnin'
8475 WRITE(*,199)
'ingWarningWarningWarningWarningWarningWarningWarningWarning'
8476 WRITE(*,199)
'ngWarningWarningWarningWarningWarningWarningWarningWarningW'
8477 WRITE(*,199)
'gWarningWarningWarningWarningWarningWarningWarningWarningWa'
8479 WRITE(*,*)
' Number of empty constraints =',abs(
ncgbe),
', should be 0'
8480 WRITE(*,*)
' => please check constraint definition, mille data'
8482 WRITE(*,199)
'WarningWarningWarningWarningWarningWarningWarningWarningWar'
8483 WRITE(*,199)
'arningWarningWarningWarningWarningWarningWarningWarningWarn'
8484 WRITE(*,199)
'rningWarningWarningWarningWarningWarningWarningWarningWarni'
8485 WRITE(*,199)
'ningWarningWarningWarningWarningWarningWarningWarningWarnin'
8486 WRITE(*,199)
'ingWarningWarningWarningWarningWarningWarningWarningWarning'
8487 WRITE(*,199)
'ngWarningWarningWarningWarningWarningWarningWarningWarningW'
8488 WRITE(*,199)
'gWarningWarningWarningWarningWarningWarningWarningWarningWa'
8493101
FORMAT(1x,a8,
' =',i14,
' = ',a)
8495103
FORMAT(1x,a,g12.4)
8496106
FORMAT(i6,2(3x,f9.3,f12.1,3x))
8511 INTEGER(mpi) :: imed
8514 INTEGER(mpi) :: nent
8515 INTEGER(mpi),
DIMENSION(measBins) :: isuml
8516 INTEGER(mpi),
DIMENSION(measBins) :: isums
8520 INTEGER(mpl) :: ioff
8523 DATA lfirst /.true./
8536 WRITE(
lunmon,
'(A)')
'*** Normalized residuals grouped by first global label (per local fit cycle) ***'
8538 WRITE(
lunmon,
'(A)')
'*** Pulls grouped by first global label (per local fit cycle) ***'
8540 WRITE(
lunmon,
'(A)')
'! LFC Label Entries Median RMS(MAD) <error>'
8559 IF (2*isuml(j) > nent)
EXIT
8563 IF (isuml(j) > isuml(j-1)) amed=amed+real(nent-2*isuml(j-1),mps)/real(2*isuml(j)-2*isuml(j-1),mps)
8576 isums(j)=isums(j)+isums(j-1)
8580 IF (2*isums(j) > nent)
EXIT
8583 IF (isums(j) > isums(j-1)) amad=amad+real(nent-2*isums(j-1),mps)/real(2*isums(j)-2*isums(j-1),mps)
8597110
FORMAT(i5,2i10,3g14.5)
8612 INTEGER(mpi) :: ioff
8613 INTEGER(mpi) :: ipar0
8614 INTEGER(mpi) :: ncon
8615 INTEGER(mpi) :: npar
8616 INTEGER(mpi) :: nextra
8618 INTEGER :: nbopt, nboptx, ILAENV
8621 INTEGER(mpl),
INTENT(IN) :: msize(2)
8623 INTEGER(mpl) :: length
8624 INTEGER(mpl) :: nwrdpc
8625 INTEGER(mpl),
PARAMETER :: three = 3
8638 CALL mpalloc(
aux,length,
' local fit scratch array: aux')
8639 CALL mpalloc(
vbnd,length,
' local fit scratch array: vbnd')
8640 CALL mpalloc(
vbdr,length,
' local fit scratch array: vbdr')
8643 CALL mpalloc(
vbk,length,
' local fit scratch array: vbk')
8646 CALL mpalloc(
vzru,length,
' local fit scratch array: vzru')
8647 CALL mpalloc(
scdiag,length,
' local fit scratch array: scdiag')
8648 CALL mpalloc(
scflag,length,
' local fit scratch array: scflag')
8649 CALL mpalloc(
ibandh,2*length,
' local fit band width hist.: ibandh')
8663 length=4*
ncblck;
IF(ncon == 0) length=0
8672 nwrdpc=nwrdpc+length
8684 length=(ncon*ncon+ncon)/2
8709 length=length+npar+nextra
8715 IF(
mbandw == 0) length=length+1
8723 nwrdpc=nwrdpc+2*length
8724 IF (nwrdpc > 250000)
THEN
8726 WRITE(*,*)
'Size of preconditioner matrix:',int(real(nwrdpc,mps)*4.0e-6,mpi),
' MB'
8748 CALL peend(23,
'Aborted, bad matrix index (will exceed 32bit)')
8749 stop
'vmprep: bad index (matrix to large for diagonalization)'
8771 nbopt = ilaenv( 1_mpl,
'DSYTRF',
'U', int(
nagb,mpl), int(
nagb,mpl), -1_mpl, -1_mpl )
8773 print *,
'LAPACK optimal block size for DSYTRF:', nbopt
8774 lplwrk=length*int(nbopt,mpl)
8782 nbopt = ilaenv( 1_mpl,
'DORMQL',
'RN', int(npar,mpl), int(npar,mpl), int(ncon,mpl), int(npar,mpl) )
8783 IF (int(npar,mpl)*int(nbopt,mpl) >
lplwrk)
THEN
8784 lplwrk=int(npar,mpl)*int(nbopt,mpl)
8789 print *,
'LAPACK optimal block size for DORMQL:', nboptx
8808 INTEGER(mpi) :: icoff
8809 INTEGER(mpi) :: ipoff
8812 INTEGER(mpi) :: ncon
8813 INTEGER(mpi) :: nfit
8814 INTEGER(mpi) :: npar
8815 INTEGER(mpi) :: nrank
8816 INTEGER(mpl) :: imoff
8817 INTEGER(mpl) :: ioff1
8835 WRITE(
lunlog,*)
'Shrinkage of global matrix (A->Q^t*A*Q)'
8850 nfit=npar+ncon;
IF (
icelim > 0) nfit=npar-ncon
8852 IF(nfit < npar)
THEN
8870 WRITE(
lunlog,*)
'Inversion of global matrix (A->A^-1)'
8877 IF(nfit /= nrank)
THEN
8878 WRITE(*,*)
'Warning: the rank defect of the symmetric',nfit, &
8879 '-by-',nfit,
' matrix is ',nfit-nrank,
' (should be zero).'
8880 WRITE(lun,*)
'Warning: the rank defect of the symmetric',nfit, &
8881 '-by-',nfit,
' matrix is ',nfit-nrank,
' (should be zero).'
8884 WRITE(*,*)
' --> enforcing SUBITO mode'
8885 WRITE(lun,*)
' --> enforcing SUBITO mode'
8887 ELSE IF(
ndefec == 0)
THEN
8889 WRITE(lun,*)
'No rank defect of the symmetric matrix'
8891 WRITE(lun,*)
'No rank defect of the symmetric block', ib,
' of size', npar
8902 IF(nfit < npar)
THEN
8921 INTEGER(mpi) :: icoff
8922 INTEGER(mpi) :: ipoff
8925 INTEGER(mpi) :: ncon
8926 INTEGER(mpi) :: nfit
8927 INTEGER(mpi) :: npar
8928 INTEGER(mpi) :: nrank
8929 INTEGER(mpl) :: imoff
8930 INTEGER(mpl) :: ioff1
8945 WRITE(
lunlog,*)
'Shrinkage of global matrix (A->Q^t*A*Q)'
8959 nfit=npar+ncon;
IF (
icelim > 0) nfit=npar-ncon
8961 IF(nfit < npar)
THEN
8979 WRITE(
lunlog,*)
'Decomposition of global matrix (A->L*D*L^t)'
8985 IF(nfit /= nrank)
THEN
8986 WRITE(*,*)
'Warning: the rank defect of the symmetric',nfit, &
8987 '-by-',nfit,
' matrix is ',nfit-nrank,
' (should be zero).'
8988 WRITE(lun,*)
'Warning: the rank defect of the symmetric',nfit, &
8989 '-by-',nfit,
' matrix is ',nfit-nrank,
' (should be zero).'
8992 WRITE(*,*)
' --> enforcing SUBITO mode'
8993 WRITE(lun,*)
' --> enforcing SUBITO mode'
8995 ELSE IF(
ndefec == 0)
THEN
8997 WRITE(lun,*)
'No rank defect of the symmetric matrix'
8999 WRITE(lun,*)
'No rank defect of the symmetric block', ib,
' of size', npar
9001 WRITE(lun,*)
' largest diagonal element (LDLt)', evmax
9002 WRITE(lun,*)
' smallest diagonal element (LDLt)', evmin
9011 IF(nfit < npar)
THEN
9033 INTEGER(mpi) :: icoff
9034 INTEGER(mpi) :: ipoff
9037 INTEGER(mpi) :: ncon
9038 INTEGER(mpi) :: nfit
9039 INTEGER(mpi) :: npar
9040 INTEGER(mpl) :: imoff
9041 INTEGER(mpl) :: ioff1
9042 INTEGER(mpi) :: infolp
9062 WRITE(
lunlog,*)
'Shrinkage of global matrix (A->Q^t*A*Q)'
9077 nfit=npar+ncon;
IF (
icelim > 0) nfit=npar-ncon
9079 IF(nfit < npar)
THEN
9096 IF (nfit > npar)
THEN
9099 WRITE(
lunlog,*)
'Factorization of global matrix (A->L*D*L^t)'
9109 WRITE(
lunlog,*)
'Factorization of global matrix (A->L*L^t)'
9113 CALL dpptrf(
'U',int(nfit,mpl),
globalmatd(imoff+1:),infolp)
9120 WRITE(lun,*)
'No rank defect of the symmetric matrix'
9122 WRITE(lun,*)
'No rank defect of the symmetric block', ib,
' of size', npar
9126 WRITE(*,*)
'Warning: factorization of the symmetric',nfit, &
9127 '-by-',nfit,
' failed at index ', infolp
9128 WRITE(lun,*)
'Warning: factorization of the symmetric',nfit, &
9129 '-by-',nfit,
' failed at index ', infolp
9130 CALL peend(29,
'Aborted, factorization of global matrix failed')
9131 stop
'mdptrf: bad matrix'
9136 IF (nfit > npar)
THEN
9139 IF(infolp /= 0) print *,
' DSPTRS failed: ', infolp
9141 CALL dpptrs(
'U',int(nfit,mpl),1_mpl,
globalmatd(imoff+1:),&
9143 IF(infolp /= 0) print *,
' DPPTRS failed: ', infolp
9147 IF(nfit < npar)
THEN
9168 INTEGER(mpi) :: icoff
9169 INTEGER(mpi) :: ipoff
9172 INTEGER(mpi) :: ncon
9173 INTEGER(mpi) :: nfit
9174 INTEGER(mpi) :: npar
9175 INTEGER(mpl) :: imoff
9176 INTEGER(mpl) :: ioff1
9177 INTEGER(mpl) :: iloff
9178 INTEGER(mpi) :: infolp
9199 WRITE(
lunlog,*)
'Shrinkage of global matrix (A->Q^t*A*Q)'
9219 nfit=npar+ncon;
IF (
icelim > 0) nfit=npar-ncon
9221 IF(nfit < npar)
THEN
9225 CALL dtrtrs(
'L',
'T',
'N',int(ncon,mpl),1_mpl,
lapackql(iloff+npar-ncon+1:),int(npar,mpl),&
9227 IF(infolp /= 0) print *,
' DTRTRS failed: ', infolp
9229 CALL dormql(
'L',
'T',int(npar,mpl),1_mpl,int(ncon,mpl),
lapackql(iloff+1:),int(npar,mpl),&
9231 IF(infolp /= 0) print *,
' DORMQL failed: ', infolp
9250 IF (nfit > npar)
THEN
9253 WRITE(
lunlog,*)
'Factorization of global matrix (A->L*D*L^t)'
9257 CALL dsytrf(
'U',int(nfit,mpl),
globalmatd(imoff+1:),int(nfit,mpl),&
9264 WRITE(
lunlog,*)
'Factorization of global matrix (A->L*L^t)'
9268 CALL dpotrf(
'U',int(nfit,mpl),
globalmatd(imoff+1:),int(npar,mpl),infolp)
9275 WRITE(lun,*)
'No rank defect of the symmetric matrix'
9277 WRITE(lun,*)
'No rank defect of the symmetric block', ib,
' of size', npar
9281 WRITE(*,*)
'Warning: factorization of the symmetric',nfit, &
9282 '-by-',nfit,
' failed at index ', infolp
9283 WRITE(lun,*)
'Warning: factorization of the symmetric',nfit, &
9284 '-by-',nfit,
' failed at index ', infolp
9285 CALL peend(29,
'Aborted, factorization of global matrix failed')
9286 stop
'mdutrf: bad matrix'
9291 IF (nfit > npar)
THEN
9292 CALL dsytrs(
'U',int(nfit,mpl),1_mpl,
globalmatd(imoff+1:),int(nfit,mpl),&
9294 IF(infolp /= 0) print *,
' DSYTRS failed: ', infolp
9296 CALL dpotrs(
'U',int(nfit,mpl),1_mpl,
globalmatd(imoff+1:),int(npar,mpl),&
9298 IF(infolp /= 0) print *,
' DPOTRS failed: ', infolp
9302 IF(nfit < npar)
THEN
9307 CALL dormql(
'L',
'N',int(npar,mpl),1_mpl,int(ncon,mpl),
lapackql(iloff+1:),int(npar,mpl),&
9309 IF(infolp /= 0) print *,
' DORMQL failed: ', infolp
9316 iloff=iloff+int(npar,mpl)*int(ncon,mpl)
9338 INTEGER(mpi) :: icboff
9339 INTEGER(mpi) :: icblst
9340 INTEGER(mpi) :: icoff
9341 INTEGER(mpi) :: icfrst
9342 INTEGER(mpi) :: iclast
9343 INTEGER(mpi) :: ipfrst
9344 INTEGER(mpi) :: iplast
9345 INTEGER(mpi) :: ipoff
9348 INTEGER(mpi) :: ncon
9349 INTEGER(mpi) :: npar
9351 INTEGER(mpl) :: imoff
9352 INTEGER(mpl) :: iloff
9353 INTEGER(mpi) :: infolp
9354 INTEGER :: nbopt, ILAENV
9356 REAL(mpd),
INTENT(IN) :: a(mszcon)
9357 REAL(mpd),
INTENT(OUT) :: emin
9358 REAL(mpd),
INTENT(OUT) :: emax
9369 iloff=iloff+int(npar,mpl)*int(ncon,mpl)
9388 DO icb=icboff+1,icboff+icblst
9395 lapackql(iloff+ipfrst:iloff+iplast)=a(imoff+1:imoff+npb)
9412 nbopt = ilaenv( 1_mpl,
'DGEQLF',
'', int(npar,mpl), int(ncon,mpl), int(npar,mpl), -1_mpl )
9413 print *,
'LAPACK optimal block size for DGEQLF:', nbopt
9414 lplwrk=int(ncon,mpl)*int(nbopt,mpl)
9417 CALL dgeqlf(int(npar,mpl),int(ncon,mpl),
lapackql(iloff+1:),int(npar,mpl),&
9419 IF(infolp /= 0) print *,
' DGEQLF failed: ', infolp
9422 iloff=iloff+int(npar,mpl)*int(ncon,mpl)
9425 IF(emax < emin)
THEN
9453 INTEGER(mpi) :: icoff
9454 INTEGER(mpi) :: ipoff
9456 INTEGER(mpi) :: ncon
9457 INTEGER(mpi) :: npar
9458 INTEGER(mpl) :: imoff
9459 INTEGER(mpl) :: iloff
9460 INTEGER(mpi) :: infolp
9461 CHARACTER (LEN=1) :: transr, transl
9463 LOGICAL,
INTENT(IN) :: t
9482 IF(ncon <= 0 ) cycle
9491 DO i=ipoff+1,ipoff+npar
9497 CALL dormql(
'R',transr,int(npar,mpl),int(npar,mpl),int(ncon,mpl),
lapackql(iloff+1:),&
9500 IF(infolp /= 0) print *,
' DORMQL failed: ', infolp
9502 CALL dormql(
'L',transl,int(npar,mpl),int(npar,mpl),int(ncon,mpl),
lapackql(iloff+1:),&
9505 IF(infolp /= 0) print *,
' DORMQL failed: ', infolp
9508 iloff=iloff+int(npar,mpl)*int(ncon,mpl)
9514include
'mkl_pardiso.f90'
9545 TYPE(mkl_pardiso_handle) :: pt(64)
9547 INTEGER(mpl),
PARAMETER :: maxfct =1
9548 INTEGER(mpl),
PARAMETER :: mnum = 1
9549 INTEGER(mpl),
PARAMETER :: nrhs = 1
9551 INTEGER(mpl) :: mtype
9552 INTEGER(mpl) :: phase
9553 INTEGER(mpl) :: error
9554 INTEGER(mpl) :: msglvl
9558 INTEGER(mpl) :: idum(1)
9560 INTEGER(mpl) :: length
9561 INTEGER(mpi) :: nfill
9562 INTEGER(mpi) :: npdblk
9563 REAL(mpd) :: adum(1)
9564 REAL(mpd) :: ddum(1)
9566 INTEGER(mpl) :: iparm(64)
9567 REAL(mpd),
ALLOCATABLE :: b( : )
9568 REAL(mpd),
ALLOCATABLE :: x( : )
9587 WRITE(*,*)
'MSPARDISO: number of rows to fill up = ', nfill
9600 CALL pardiso_64(pt, maxfct, mnum, mtype, phase, int(npdblk,mpl), adum, idum, idum, &
9601 idum, nrhs, iparm, msglvl, ddum, ddum, error)
9602 IF (error /= 0)
THEN
9603 WRITE(lun,*)
'The following ERROR was detected: ', error
9604 WRITE(*,
'(A,2I10)')
' PARDISO release failed (phase, error): ', phase, error
9605 IF (
ipddbg == 0)
WRITE(*,*)
' rerun with "debugPARDISO" for more info'
9606 CALL peend(40,
'Aborted, other error: PARDISO release')
9607 stop
'MSPARDISO: stopping due to error in PARDISO release'
9624 IF (iparm(1) == 0)
WRITE(lun,*)
'PARDISO using defaults '
9625 IF (iparm(43) /= 0)
THEN
9626 WRITE(lun,*)
'PARDISO: computation of the diagonal of inverse matrix not implemented !'
9640 WRITE(
lunlog,*)
'Decomposition of global matrix (A->L*D*L^t)'
9654 WRITE(lun,*)
' iparm(',i,
') =', iparm(i)
9657 CALL pardiso_64(pt, maxfct, mnum, mtype, phase, int(npdblk,mpl),
globalmatd,
csr3rowoffsets,
csr3columnlist, &
9658 idum, nrhs, iparm, msglvl, ddum, ddum, error)
9660 WRITE(lun,*)
'PARDISO reordering completed ... '
9661 WRITE(lun,*)
'PARDISO peak memory required (KB)', iparm(15)
9664 WRITE(lun,*)
' iparm(',i,
') =', iparm(i)
9667 IF (error /= 0)
THEN
9668 WRITE(lun,*)
'The following ERROR was detected: ', error
9669 WRITE(*,
'(A,2I10)')
' PARDISO decomposition failed (phase, error): ', phase, error
9670 IF (
ipddbg == 0)
WRITE(*,*)
' rerun with "debugPARDISO" for more info'
9671 CALL peend(40,
'Aborted, other error: PARDISO reordering')
9672 stop
'MSPARDISO: stopping due to error in PARDISO reordering'
9674 IF (iparm(60) == 0)
THEN
9679 WRITE(lun,*)
'Size (KB) of allocated memory = ',
ipdmem
9680 WRITE(lun,*)
'Number of nonzeros in factors = ',iparm(18)
9681 WRITE(lun,*)
'Number of factorization MFLOPS = ',iparm(19)
9686 CALL pardiso_64(pt, maxfct, mnum, mtype, phase, int(npdblk,mpl),
globalmatd,
csr3rowoffsets,
csr3columnlist, &
9687 idum, nrhs, iparm, msglvl, ddum, ddum, error)
9689 WRITE(lun,*)
'PARDISO factorization completed ... '
9692 WRITE(lun,*)
' iparm(',i,
') =', iparm(i)
9695 IF (error /= 0)
THEN
9696 WRITE(lun,*)
'The following ERROR was detected: ', error
9697 WRITE(*,
'(A,2I10)')
' PARDISO decomposition failed (phase, error): ', phase, error
9698 IF (
ipddbg == 0)
WRITE(*,*)
' rerun with "debugPARDISO" for more info'
9699 CALL peend(40,
'Aborted, other error: PARDISO factorization')
9700 stop
'MSPARDISO: stopping due to error in PARDISO factorization'
9703 IF (iparm(14) > 0) &
9704 WRITE(lun,*)
'Number of perturbed pivots = ',iparm(14)
9705 WRITE(lun,*)
'Number of positive eigenvalues = ',iparm(22)-nfill
9706 WRITE(lun,*)
'Number of negative eigenvalues = ',iparm(23)
9707 ELSE IF (iparm(30) > 0)
THEN
9708 WRITE(lun,*)
'Equation with bad pivot (<=0.) = ',iparm(30)
9717 CALL mpalloc(b,length,
' PARDISO r.h.s')
9718 CALL mpalloc(x,length,
' PARDISO solution')
9723 CALL pardiso_64(pt, maxfct, mnum, mtype, phase, int(npdblk,mpl),
globalmatd,
csr3rowoffsets,
csr3columnlist, &
9724 idum, nrhs, iparm, msglvl, b, x, error)
9729 WRITE(lun,*)
'PARDISO solve completed ... '
9730 IF (error /= 0)
THEN
9731 WRITE(lun,*)
'The following ERROR was detected: ', error
9732 WRITE(*,
'(A,2I10)')
' PARDISO decomposition failed (phase, error): ', phase, error
9733 IF (
ipddbg == 0)
WRITE(*,*)
' rerun with "debugPARDISO" for more info'
9734 CALL peend(40,
'Aborted, other error: PARDISO solve')
9735 stop
'MSPARDISO: stopping due to error in PARDISO solve'
9749 INTEGER(mpi) :: iast
9750 INTEGER(mpi) :: idia
9751 INTEGER(mpi) :: imin
9752 INTEGER(mpl) :: ioff1
9754 INTEGER(mpi) :: last
9756 INTEGER(mpi) :: nmax
9757 INTEGER(mpi) :: nmin
9758 INTEGER(mpi) :: ntop
9780 WRITE(
lunlog,*)
'Shrinkage of global matrix (A->Q^t*A*Q)'
9817 DO WHILE(ntop < nmax)
9821 CALL hmpdef(7,real(nmin,mps),real(ntop,mps),
'log10 of positive eigenvalues')
9832 CALL gmpdef(3,2,
'low-value end of eigenvalues')
9835 CALL gmpxy(3,real(i,mps),evalue)
9849 WRITE(lun,*)
'The first (largest) eigenvalues ...'
9852 WRITE(lun,*)
'The last eigenvalues ... up to',last
9856 WRITE(lun,*)
'The eigenvalues from',
nvgb+1,
' to',
nagb
9860 WRITE(lun,*)
'Log10 + 3 of ',
nfgb,
' eigenvalues in decreasing',
' order'
9861 WRITE(lun,*)
'(for Eigenvalue < 0.001 the value 0.0 is shown)'
9864 'printed for negative eigenvalues'
9867 WRITE(lun,*) last,
' significances: insignificant if ', &
9868 'compatible with N(0,1)'
9897 INTEGER(mpl) :: ioff1
9898 INTEGER(mpl) :: ioff2
9934 INTEGER(mpi) :: istop
9936 INTEGER(mpi) :: itnlim
9938 INTEGER(mpi) :: nout
9939 INTEGER(mpi) :: nrkd
9940 INTEGER(mpi) :: nrkd2
9985 WRITE(lun,*)
'MMINRS: PRECONS ended ', nrkd
9989 globalcorrections, itnlim, nout, rtol, istop, itn, anorm, acond, rnorm, arnorm, ynorm)
9996 WRITE(lun,*)
'MMINRS: EQUDECS ended ', nrkd, nrkd2
10000 globalcorrections, itnlim, nout, rtol, istop, itn, anorm, acond, rnorm, arnorm, ynorm)
10003 globalcorrections, itnlim, nout, rtol, istop, itn, anorm, acond, rnorm, arnorm, ynorm)
10017 IF (
istopa == 0) print *,
'MINRES: istop=0, exact solution x=0.'
10032 INTEGER(mpi) :: istop
10033 INTEGER(mpi) :: itn
10034 INTEGER(mpi) :: itnlim
10035 INTEGER(mpi) :: lun
10036 INTEGER(mpi) :: nout
10037 INTEGER(mpi) :: nrkd
10038 INTEGER(mpi) :: nrkd2
10041 REAL(mpd) :: mxxnrm
10042 REAL(mpd) :: trcond
10052 mxxnrm = real(
nagb,mpd)/sqrt(epsilon(mxxnrm))
10054 trcond = 1.0_mpd/epsilon(trcond)
10055 ELSE IF(
mrmode == 2)
THEN
10085 WRITE(lun,*)
'MMINRS: PRECONS ended ', nrkd
10089 itnlim=itnlim, rtol=rtol, maxxnorm=mxxnrm, trancond=trcond, &
10091 ELSE IF(
mbandw > 0)
THEN
10097 WRITE(lun,*)
'MMINRS: EQUDECS ended ', nrkd, nrkd2
10102 itnlim=itnlim, rtol=rtol, maxxnorm=mxxnrm, trancond=trcond, &
10106 itnlim=itnlim, rtol=rtol, maxxnorm=mxxnrm, trancond=trcond, &
10121 IF (
istopa == 3) print *,
'MINRES: istop=0, exact solution x=0.'
10137 INTEGER(mpi),
INTENT(IN) :: n
10138 REAL(mpd),
INTENT(IN) :: x(n)
10139 REAL(mpd),
INTENT(OUT) :: y(n)
10159 INTEGER(mpi),
INTENT(IN) :: n
10160 REAL(mpd),
INTENT(IN) :: x(n)
10161 REAL(mpd),
INTENT(OUT) :: y(n)
10192 REAL(mps) :: concu2
10193 REAL(mps) :: concut
10194 REAL,
DIMENSION(2) :: ta
10197 INTEGER(mpi) :: iact
10198 INTEGER(mpi) :: iagain
10199 INTEGER(mpi) :: idx
10200 INTEGER(mpi) :: info
10202 INTEGER(mpi) :: ipoff
10203 INTEGER(mpi) :: icoff
10204 INTEGER(mpl) :: ioff
10205 INTEGER(mpi) :: itgbi
10206 INTEGER(mpi) :: ivgbi
10207 INTEGER(mpi) :: jcalcm
10209 INTEGER(mpi) :: labelg
10210 INTEGER(mpi) :: litera
10211 INTEGER(mpi) :: lrej
10212 INTEGER(mpi) :: lun
10213 INTEGER(mpi) :: lunp
10214 INTEGER(mpi) :: minf
10215 INTEGER(mpi) :: mrati
10216 INTEGER(mpi) :: nan
10217 INTEGER(mpi) :: ncon
10218 INTEGER(mpi) :: nfaci
10219 INTEGER(mpi) :: nloopsol
10220 INTEGER(mpi) :: npar
10221 INTEGER(mpi) :: nrati
10222 INTEGER(mpi) :: nrej
10223 INTEGER(mpi) :: nsol
10224 INTEGER(mpi) :: inone
10226 INTEGER(mpi) :: infolp
10227 INTEGER(mpi) :: nfit
10228 INTEGER(mpl) :: imoff
10232 REAL(mpd) :: dratio
10233 REAL(mpd) :: dwmean
10242 LOGICAL :: warnerss
10243 LOGICAL :: warners3
10245 CHARACTER (LEN=7) :: cratio
10246 CHARACTER (LEN=7) :: cfacin
10247 CHARACTER (LEN=7) :: crjrat
10258 WRITE(lunp,*)
'Solution algorithm: '
10259 WRITE(lunp,121)
'=================================================== '
10262 WRITE(lunp,121)
'solution method:',
'matrix inversion'
10263 ELSE IF(
metsol == 2)
THEN
10264 WRITE(lunp,121)
'solution method:',
'diagonalization'
10265 ELSE IF(
metsol == 3)
THEN
10266 WRITE(lunp,121)
'solution method:',
'decomposition'
10267 ELSE IF(
metsol == 4)
THEN
10268 WRITE(lunp,121)
'solution method:',
'minres (Paige/Saunders)'
10269 ELSE IF(
metsol == 5)
THEN
10270 WRITE(lunp,121)
'solution method:',
'minres-qlp (Choi/Paige/Saunders)'
10272 WRITE(lunp,121)
' ',
' using QR factorization'
10273 ELSE IF(
mrmode == 2)
THEN
10274 WRITE(lunp,121)
' ',
' using QLP factorization'
10276 WRITE(lunp,121)
' ',
' using QR and QLP factorization'
10277 WRITE(lunp,123)
'transition condition',
mrtcnd
10279 ELSE IF(
metsol == 6)
THEN
10280 WRITE(lunp,121)
'solution method:', &
10281 'gmres (generalized minimzation of residuals)'
10283 ELSE IF(
metsol == 7)
THEN
10285 WRITE(lunp,121)
'solution method:',
'LAPACK factorization (DSPTRF)'
10287 WRITE(lunp,121)
'solution method:',
'LAPACK factorization (DPPTRF)'
10289 IF(
ilperr == 1)
WRITE(lunp,121)
' ',
'with error calculation (D??TRI)'
10290 ELSE IF(
metsol == 8)
THEN
10292 WRITE(lunp,121)
'solution method:',
'LAPACK factorization (DSYTRF)'
10294 WRITE(lunp,121)
'solution method:',
'LAPACK factorization (DPOTRF)'
10296 IF(
ilperr == 1)
WRITE(lunp,121)
' ',
'with error calculation (D??TRI)'
10298 ELSE IF(
metsol == 9)
THEN
10300 WRITE(lunp,121)
'solution method:',
'Intel oneMKL PARDISO (sparse matrix (CSR3))'
10302 WRITE(lunp,121)
'solution method:',
'Intel oneMKL PARDISO (sparse matrix (BSR3))'
10307 WRITE(lunp,123)
'convergence limit at Delta F=',
dflim
10308 WRITE(lunp,122)
'maximum number of iterations=',
mitera
10311 WRITE(lunp,122)
'matrix recalculation up to ',
matrit,
'. iteration'
10315 WRITE(lunp,121)
'matrix storage:',
'full'
10316 ELSE IF(
matsto == 2)
THEN
10317 WRITE(lunp,121)
'matrix storage:',
'sparse'
10319 WRITE(lunp,122)
'pre-con band-width parameter=',
mbandw
10321 WRITE(lunp,121)
'pre-conditioning:',
'default'
10322 ELSE IF(
mbandw < 0)
THEN
10323 WRITE(lunp,121)
'pre-conditioning:',
'none!'
10324 ELSE IF(
mbandw > 0)
THEN
10326 WRITE(lunp,121)
'pre-conditioning=',
'skyline-matrix (rank preserving)'
10328 WRITE(lunp,121)
'pre-conditioning=',
'band-matrix'
10333 WRITE(lunp,121)
'using pre-sigmas:',
'no'
10336 WRITE(lunp,124)
'pre-sigmas defined for', &
10337 REAL(100*npresg,mps)/
REAL(nvgb,mps),
' % of variable parameters'
10338 WRITE(lunp,123)
'default pre-sigma=',
regpre
10341 WRITE(lunp,121)
'regularization:',
'no'
10343 WRITE(lunp,121)
'regularization:',
'yes'
10344 WRITE(lunp,123)
'regularization factor=',
regula
10348 WRITE(lunp,121)
'Chi square cut equiv 3 st.dev applied'
10349 WRITE(lunp,123)
'... in first iteration with factor',
chicut
10350 WRITE(lunp,123)
'... in second iteration with factor',
chirem
10351 WRITE(lunp,121)
' (reduced by sqrt in next iterations)'
10354 WRITE(lunp,121)
'Scaling of measurement errors applied'
10355 WRITE(lunp,123)
'... factor for "global" measuements',
dscerr(1)
10356 WRITE(lunp,123)
'... factor for "local" measuements',
dscerr(2)
10359 WRITE(lunp,122)
'Down-weighting of outliers in',
lhuber,
' iterations'
10360 WRITE(lunp,123)
'Cut on downweight fraction',
dwcut
10364121
FORMAT(1x,a40,3x,a)
10365122
FORMAT(1x,a40,3x,i0,a)
10366123
FORMAT(1x,a40,2x,e9.2)
10367124
FORMAT(1x,a40,3x,f5.1,a)
10377 stepl =real(stp,mps)
10389 ELSE IF(
metsol == 2)
THEN
10392 ELSE IF(
metsol == 3)
THEN
10395 ELSE IF(
metsol == 4)
THEN
10398 ELSE IF(
metsol == 5)
THEN
10401 ELSE IF(
metsol == 6)
THEN
10413 WRITE(
lunlog,*)
'Checking feasibility of parameters:'
10414 WRITE(*,*)
'Checking feasibility of parameters:'
10415 CALL feasib(concut,iact)
10417 WRITE(*,102) concut
10418 WRITE(*,*)
' parameters are made feasible'
10419 WRITE(
lunlog,102) concut
10420 WRITE(
lunlog,*)
' parameters are made feasible'
10422 WRITE(*,*)
' parameters are feasible (i.e. satisfy constraints)'
10423 WRITE(
lunlog,*)
' parameters are feasible (i.e. satisfy constraints)'
10431 WRITE(*,*)
'Reading files and accumulating vectors/matrices ...'
10435 WRITE(
lunlog,*)
'Reading files and accumulating vectors/matrices ...'
10452 IF(jcalcm+1 /= 0)
THEN
10461 IF(
iterat /= litera)
THEN
10482 IF (iabs(jcalcm) <= 1)
THEN
10495 IF(3*nrej >
nrecal)
THEN
10497 WRITE(*,*)
'Data rejected in previous loop: '
10499 nrejec(0),
' (rank deficit/NaN) ',
nrejec(1),
' (Ndf=0) ', &
10501 WRITE(*,*)
'Too many rejects (>33.3%) - stop'
10502 CALL peend(26,
'Aborted, too many rejects')
10509 IF(abs(
icalcm) == 1)
THEN
10521 ELSE IF(
metsol == 2)
THEN
10523 ELSE IF(
metsol == 3)
THEN
10525 ELSE IF(
metsol == 4)
THEN
10527 ELSE IF(
metsol == 5)
THEN
10529 ELSE IF(
metsol == 6)
THEN
10530 WRITE(*,*)
'... reserved for GMRES (not yet!)'
10533 ELSE IF(
metsol == 7)
THEN
10535 ELSE IF(
metsol == 8)
THEN
10538 ELSE IF(
metsol == 9)
THEN
10552 CALL feasib(concut,iact)
10564 angras=real(db/sqrt(db1*db2),mps)
10565 dbsig=16.0_mpd*sqrt(max(db1,db2))*epsilon(db)
10571 lsflag=lsflag .AND. (db > dbsig)
10578 ELSE IF(
metsol == 2)
THEN
10581 ELSE IF(
metsol == 3)
THEN
10584 ELSE IF(
metsol == 4)
THEN
10587 ELSE IF(
metsol == 5)
THEN
10590 ELSE IF(
metsol == 6)
THEN
10600 IF(db <= -dbsig)
THEN
10601 WRITE(*,*)
'Function not decreasing:',db
10602 IF(db > -1.0e-3_mpd)
THEN
10604 IF (iagain <= 1)
THEN
10605 WRITE(*,*)
'... again matrix calculation'
10609 WRITE(*,*)
'... aborting iterations'
10613 WRITE(*,*)
'... stopping iterations'
10636 IF (
nloopn == nloopsol)
THEN
10642 stepl=real(stp,mps)
10652 WRITE(*,*)
'Result vector containes ', nan,
' NaNs - stop'
10653 CALL peend(25,
'Aborted, result vector contains NaNs')
10660 WRITE(*,*)
'Subito! Exit after first step.'
10665 WRITE(*,*)
'INFO=0 should not happen (line search input err)'
10666 IF (iagain <= 0)
THEN
10671 IF(info < 0 .OR.
nloopn == nloopsol) cycle
10675 CALL feasib(concut,iact)
10676 IF(iact /= 0.OR.
chicut > 1.0)
THEN
10693 WRITE(*,*)
'Data rejected in last loop: '
10695 nrejec(0),
' (rank deficit/NaN) ',
nrejec(1),
' (Ndf=0) ', &
10717 nfit=npar+ncon;
IF (
icelim > 0) nfit=npar-ncon
10718 IF (nfit > npar)
THEN
10721 WRITE(
lunlog,*)
'Inverse of global matrix from LDLt factorization'
10727 IF(infolp /= 0) print *,
' DSPTRI failed: ', infolp
10732 CALL dsytri(
'U',int(nfit,mpl),
globalmatd(imoff+1:),int(nfit,mpl),&
10734 IF(infolp /= 0) print *,
' DSYTRI failed: ', infolp
10740 WRITE(
lunlog,*)
'Inverse of global matrix from LLt factorization'
10745 CALL dpptri(
'U',int(nfit,mpl),
globalmatd(imoff+1:),infolp)
10746 IF(infolp /= 0) print *,
' DPPTRI failed: ', infolp
10750 CALL dpotri(
'U',int(nfit,mpl),
globalmatd(imoff+1:),int(npar,mpl),infolp)
10751 IF(infolp /= 0) print *,
' DPOTRI failed: ', infolp
10768 DO i=npar-ncon+1,npar
10775 WRITE(
lunlog,*)
'Expansion of global matrix (A->Q*A*Q^t)'
10791 catio=real(dratio,mps)
10795 mrati=nint(100.0*catio,mpi)
10799 IF (
nfilw <= 0)
THEN
10800 WRITE(lunp,*)
'Sum(Chi^2)/Sum(Ndf) =',
fvalue
10802 WRITE(lunp,*)
' =',dratio
10804 WRITE(lunp,*)
'Sum(W*Chi^2)/Sum(Ndf)/<W> =',
fvalue
10806 WRITE(lunp,*)
' /',dwmean
10807 WRITE(lunp,*)
' =',dratio
10811 ' with correction for down-weighting ',catio
10832 nrati=nint(10000.0*real(nrej,mps)/real(
nrecal,mps),mpi)
10833 WRITE(crjrat,197) 0.01_mpd*real(nrati,mpd)
10834 nfaci=nint(100.0*sqrt(catio),mpi)
10836 WRITE(cratio,197) 0.01_mpd*real(mrati,mpd)
10837 WRITE(cfacin,197) 0.01_mpd*real(nfaci,mpd)
10840 IF(mrati < 90.OR.mrati > 110) warner=.true.
10841 IF(nrati > 100) warner=.true.
10842 IF(
ncgbe /= 0) warner=.true.
10844 IF(
nalow /= 0) warners=.true.
10846 IF(
nmiss1 /= 0) warnerss=.true.
10847 IF(iagain /= 0) warnerss=.true.
10848 IF(
ndefec /= 0) warnerss=.true.
10849 IF(
ndefpg /= 0) warnerss=.true.
10851 IF(
nrderr /= 0) warners3=.true.
10853 IF(warner.OR.warners.OR.warnerss.Or.warners3)
THEN
10856 WRITE(*,199)
'WarningWarningWarningWarningWarningWarningWarningWarningWar'
10857 WRITE(*,199)
'arningWarningWarningWarningWarningWarningWarningWarningWarn'
10858 WRITE(*,199)
'rningWarningWarningWarningWarningWarningWarningWarningWarni'
10859 WRITE(*,199)
'ningWarningWarningWarningWarningWarningWarningWarningWarnin'
10860 WRITE(*,199)
'ingWarningWarningWarningWarningWarningWarningWarningWarning'
10861 WRITE(*,199)
'ngWarningWarningWarningWarningWarningWarningWarningWarningW'
10862 WRITE(*,199)
'gWarningWarningWarningWarningWarningWarningWarningWarningWa'
10864 IF(mrati < 90.OR.mrati > 110)
THEN
10866 WRITE(*,*)
' Chi^2/Ndf = ',cratio,
' (should be close to 1)'
10867 WRITE(*,*)
' => multiply all input standard ', &
10868 'deviations by factor',cfacin
10871 IF(nrati > 100)
THEN
10873 WRITE(*,*)
' Fraction of rejects =',crjrat,
' %', &
10874 ' (should be far below 1 %)'
10875 WRITE(*,*)
' => please provide correct mille data'
10879 IF(iagain /= 0)
THEN
10881 WRITE(*,*)
' Matrix not positiv definite '// &
10882 '(function not decreasing)'
10883 WRITE(*,*)
' => please provide correct mille data'
10888 WRITE(*,*)
' Rank defect =',
ndefec, &
10889 ' for global matrix, should be 0'
10890 WRITE(*,*)
' => please provide correct mille data'
10895 WRITE(*,*)
' Rank defect for',
ndefpg, &
10896 ' parameter groups, should be 0'
10897 WRITE(*,*)
' => please provide correct mille data'
10902 WRITE(*,*)
' Rank defect =',
nmiss1, &
10903 ' for constraint equations, should be 0'
10904 WRITE(*,*)
' => please correct constraint definition'
10907 IF(
ncgbe /= 0)
THEN
10909 WRITE(*,*)
' Number of empty constraints =',
ncgbe,
', should be 0'
10910 WRITE(*,*)
' => please check constraint definition, mille data'
10913 IF(
nxlow /= 0)
THEN
10915 WRITE(*,*)
' Possible rank defects =',
nxlow,
' for global matrix'
10916 WRITE(*,*)
' (too few accepted entries)'
10917 WRITE(*,*)
' => please check mille data and ENTRIES cut'
10920 IF(
nalow /= 0)
THEN
10922 WRITE(*,*)
' Possible bad elements =',
nalow,
' in global vector'
10923 WRITE(*,*)
' (toos few accepted entries)'
10924 IF(
ipcntr > 0)
WRITE(*,*)
' (indicated in millepede.res by counts<0)'
10925 WRITE(*,*)
' => please check mille data and ENTRIES cut'
10930 WRITE(*,*)
' Binary file(s) with read errors =',
nrderr,
' (treated as EOF)'
10931 WRITE(*,*)
' => please check mille data'
10935 WRITE(*,199)
'WarningWarningWarningWarningWarningWarningWarningWarningWar'
10936 WRITE(*,199)
'arningWarningWarningWarningWarningWarningWarningWarningWarn'
10937 WRITE(*,199)
'rningWarningWarningWarningWarningWarningWarningWarningWarni'
10938 WRITE(*,199)
'ningWarningWarningWarningWarningWarningWarningWarningWarnin'
10939 WRITE(*,199)
'ingWarningWarningWarningWarningWarningWarningWarningWarning'
10940 WRITE(*,199)
'ngWarningWarningWarningWarningWarningWarningWarningWarningW'
10941 WRITE(*,199)
'gWarningWarningWarningWarningWarningWarningWarningWarningWa'
10952 ELSE IF(
metsol == 2)
THEN
10954 ELSE IF(
metsol == 3)
THEN
10960 IF(labelg == 0) cycle
10961 itgbi=inone(labelg)
10964 IF(ivgbi < 0) ivgbi=0
10965 IF(ivgbi == 0) cycle
10974 ELSE IF(
metsol == 6)
THEN
10977 ELSE IF(
metsol == 7)
THEN
10985 CALL peend(4,
'Ended with severe warnings (bad binary file(s))')
10986 ELSE IF (warnerss)
THEN
10987 CALL peend(3,
'Ended with severe warnings (bad global matrix)')
10988 ELSE IF (warners)
THEN
10989 CALL peend(2,
'Ended with severe warnings (insufficient measurements)')
10990 ELSE IF (warner)
THEN
10991 CALL peend(1,
'Ended with warnings (bad measurements)')
10993 CALL peend(0,
'Ended normally')
10996102
FORMAT(
' Call FEASIB with cut=',g10.3)
11014 INTEGER(mpi) :: kfl
11015 INTEGER(mpi) :: kmin
11016 INTEGER(mpi) :: kmax
11017 INTEGER(mpi) :: nrc
11018 INTEGER(mpi) :: nrej
11024 REAL(mpd) :: sumallw
11025 REAL(mpd) :: sumrejw
11027 sumallw=0.; sumrejw=0.;
11036 sumallw=sumallw+real(nrc,mpd)*
wfd(kfl)
11037 sumrejw=sumrejw+real(nrej,mpd)*
wfd(kfl)
11038 frac=real(nrej,mps)/real(nrc,mps)
11039 IF (frac > fmax)
THEN
11043 IF (frac < fmin)
THEN
11050 WRITE(*,
"(' Weighted fraction =',F8.2,' %')") 100.*sumrejw/sumallw
11051 IF (
nfilb > 1)
THEN
11052 WRITE(*,
"(' File with max. fraction ',I6,' :',F8.2,' %')") kmax, 100.*fmax
11053 WRITE(*,
"(' File with min. fraction ',I6,' :',F8.2,' %')") kmin, 100.*fmin
11079 INTEGER(mpi) :: iargc
11082 INTEGER(mpi) :: ierrf
11083 INTEGER(mpi) :: ieq
11084 INTEGER(mpi) :: ifilb
11085 INTEGER(mpi) :: ioff
11086 INTEGER(mpi) :: iopt
11087 INTEGER(mpi) :: ios
11088 INTEGER(mpi) :: iosum
11091 INTEGER(mpi) :: mat
11092 INTEGER(mpi) :: nab
11093 INTEGER(mpi) :: nline
11094 INTEGER(mpi) :: npat
11095 INTEGER(mpi) :: ntext
11097 INTEGER(mpi) :: nuf
11098 INTEGER(mpi) :: nums
11099 INTEGER(mpi) :: nufile
11100 INTEGER(mpi) :: lenfileInfo
11101 INTEGER(mpi) :: lenFileNames
11102 INTEGER(mpi) :: matint
11103 INTEGER(mpi),
DIMENSION(:,:),
ALLOCATABLE :: vecfileInfo
11104 INTEGER(mpi),
DIMENSION(:,:),
ALLOCATABLE :: tempArray
11105 INTEGER(mpl) :: rows
11106 INTEGER(mpl) :: cols
11107 INTEGER(mpl) :: newcols
11108 INTEGER(mpl) :: length
11110 CHARACTER (LEN=1024) :: text
11111 CHARACTER (LEN=1024) :: fname
11112 CHARACTER (LEN=14) :: bite(3)
11113 CHARACTER (LEN=32) :: keystx
11114 INTEGER(mpi),
PARAMETER :: mnum=100
11115 REAL(mpd) :: dnum(mnum)
11119 SUBROUTINE initc(nfiles)
BIND(c)
11121 INTEGER(c_int),
INTENT(IN),
VALUE :: nfiles
11122 END SUBROUTINE initc
11127 DATA bite/
'C_binary',
'text ',
'Fortran_binary'/
11142 WRITE(*,*)
'Command line options: '
11143 WRITE(*,*)
'--------------------- '
11145 CALL getarg(i,text)
11146 CALL rltext(text,ia,ib,nab)
11147 WRITE(*,101) i,text(1:nab)
11148 IF(text(ia:ia) /=
'-')
THEN
11149 nu=nufile(text(ia:ib))
11152 WRITE(*,*)
'Second text file in command line - stop'
11153 CALL peend(12,
'Aborted, second text file in command line')
11159 WRITE(*,*)
'Open error for file:',text(ia:ib),
' - stop'
11160 CALL peend(16,
'Aborted, open error for file')
11161 IF(text(ia:ia) /=
'/')
THEN
11162 CALL getenv(
'PWD',text)
11163 CALL rltext(text,ia,ib,nab)
11164 WRITE(*,*)
'PWD:',text(ia:ib)
11169 IF(index(text(ia:ib),
'b') /= 0)
THEN
11171 WRITE(*,*)
'Debugging requested'
11173 it=index(text(ia:ib),
't')
11176 ieq=index(text(ia+it:ib),
'=')+it
11177 IF (it /= ieq)
THEN
11178 IF (index(text(ia+ieq:ib),
'SL0' ) /= 0)
ictest=2
11179 IF (index(text(ia+ieq:ib),
'SLE' ) /= 0)
ictest=3
11180 IF (index(text(ia+ieq:ib),
'BP' ) /= 0)
ictest=4
11181 IF (index(text(ia+ieq:ib),
'BRLF') /= 0)
ictest=5
11182 IF (index(text(ia+ieq:ib),
'BRLC') /= 0)
ictest=6
11185 IF(index(text(ia:ib),
's') /= 0)
isubit=1
11186 IF(index(text(ia:ib),
'f') /= 0)
iforce=1
11187 IF(index(text(ia:ib),
'c') /= 0)
icheck=1
11188 IF(index(text(ia:ib),
'C') /= 0)
icheck=2
11190 IF(i == iargc())
WRITE(*,*)
'--------------------- '
11211 CALL rltext(text,ia,ib,nab)
11212 nu=nufile(text(ia:ib))
11216 CALL peend(10,
'Aborted, no steering file')
11217 stop
'in FILETC: no steering file. .'
11230 WRITE(*,*)
'Listing of steering file: ',
filnam(1:
nfnam)
11231 WRITE(*,*)
'-------------------------'
11234 WRITE(*,*)
'Open error for steering file - stop'
11235 CALL peend(11,
'Aborted, open error for steering file')
11236 IF(
filnam(1:1) /=
'/')
THEN
11237 CALL getenv(
'PWD',text)
11238 CALL rltext(text,ia,ib,nab)
11239 WRITE(*,*)
'PWD:',text(ia:ib)
11248 rows=6; cols=lenfileinfo
11249 CALL mpalloc(vecfileinfo,rows,cols,
'file info from steering')
11252 READ(10,102,iostat=ierrf) text
11253 IF (ierrf < 0)
EXIT
11254 CALL rltext(text,ia,ib,nab)
11256 IF(nline <= 50)
THEN
11257 WRITE(*,101) nline,text(1:nab)
11258 IF(nline == 50)
WRITE(*,*)
' ...'
11262 CALL rltext(text,ia,ib,nab)
11263 IF(ib == ia+2)
THEN
11264 mat=matint(text(ia:ib),
'end',npat,ntext)
11265 IF(mat == max(npat,ntext))
THEN
11268 WRITE(*,*)
' end-statement after',nline,
' text lines'
11273 keystx=
'fortranfiles'
11274 mat=matint(text(ia:ib),keystx,npat,ntext)
11275 IF(mat == max(npat,ntext))
THEN
11282 mat=matint(text(ia:ib),keystx,npat,ntext)
11283 IF(mat == max(npat,ntext))
THEN
11289 keystx=
'closeandreopen'
11290 mat=matint(text(ia:ib),keystx,npat,ntext)
11291 IF(mat == max(npat,ntext))
THEN
11299 iopt=index(text(ia:ib),
' -- ')
11300 IF (iopt > 0) ie=iopt-1
11303 nu=nufile(text(ia:ie))
11305 IF (
nfiles == lenfileinfo)
THEN
11306 CALL mpalloc(temparray,rows,cols,
'temp file info from steering')
11307 temparray=vecfileinfo
11309 lenfileinfo=lenfileinfo*2
11310 newcols=lenfileinfo
11311 CALL mpalloc(vecfileinfo,rows,newcols,
'file info from steering')
11312 vecfileinfo(:,1:cols)=temparray(:,1:cols)
11318 lenfilenames=lenfilenames+ie-ia+1
11319 vecfileinfo(1,
nfiles)=nline
11320 vecfileinfo(2,
nfiles)=nu
11321 vecfileinfo(3,
nfiles)=ia
11322 vecfileinfo(4,
nfiles)=ie
11323 vecfileinfo(5,
nfiles)=iopt
11324 vecfileinfo(6,
nfiles)=ib
11334 CALL mpalloc(
nfd,length,
'file line (in steering)')
11337 length=lenfilenames
11343 READ(10,102,iostat=ierrf) text
11344 IF (ierrf < 0)
EXIT
11346 IF (nline == vecfileinfo(1,i))
THEN
11347 nfd(i)=vecfileinfo(1,i)
11348 mfd(i)=vecfileinfo(2,i)
11349 ia=vecfileinfo(3,i)-1
11350 lfd(i)=vecfileinfo(4,i)-ia
11352 tfd(ioff+k)=text(ia+k:ia+k)
11357 IF (vecfileinfo(5,i) > 0)
THEN
11358 CALL ratext(text(vecfileinfo(5,i)+4:vecfileinfo(6,i)),nums,dnum,mnum)
11359 IF (nums > 0)
ofd(i)=real(dnum(1),mps)
11369 CALL mpalloc(
ifd,length,
'integrated record numbers (=offset)')
11370 CALL mpalloc(
jfd,length,
'number of accepted records')
11371 CALL mpalloc(
kfd,rows,length,
'number of records in file, file order')
11376 CALL mpalloc(
sfd,rows,length,
'start, end of file name in TFD')
11377 CALL mpalloc(
yfd,length,
'modification date')
11380 WRITE(*,*)
'-------------------------'
11386 WRITE(*,*)
'Table of files:'
11387 WRITE(*,*)
'---------------'
11390 WRITE(8,*)
'Text and data files:'
11394 fname(k:k)=
tfd(ioff+k)
11397 IF (
mprint > 1)
WRITE(*,103) i,bite(
mfd(i)),fname(1:
lfd(i))
11398 WRITE(8,103) i,bite(
mfd(i)),fname(1:
lfd(i))
11402 WRITE(*,*)
'---------------'
11416 IF(
mfd(i) == 3)
THEN
11440 IF(
mfd(i) == 1)
THEN
11461 WRITE(*,*)
'Opening of C-files not supported.'
11477 IF(iosum /= 0)
THEN
11478 CALL peend(15,
'Aborted, open error(s) for binary files')
11479 stop
'FILETC: open error '
11481 IF(
nfilb == 0)
THEN
11482 CALL peend(14,
'Aborted, no binary files')
11483 stop
'FILETC: no binary files '
11486 WRITE(*,*)
nfilb,
' binary files opened'
11488 WRITE(*,*)
nfilb,
' binary files opened and closed'
11492103
FORMAT(i3,2x,a14,3x,a)
11555 INTEGER(mpi) :: ierrf
11556 INTEGER(mpi) :: ioff
11557 INTEGER(mpi) :: ios
11558 INTEGER(mpi) :: iosum
11560 INTEGER(mpi) :: mat
11561 INTEGER(mpi) :: nab
11562 INTEGER(mpi) :: nfiln
11563 INTEGER(mpi) :: nline
11564 INTEGER(mpi) :: nlinmx
11565 INTEGER(mpi) :: npat
11566 INTEGER(mpi) :: ntext
11567 INTEGER(mpi) :: matint
11571 CHARACTER (LEN=1024) :: text
11572 CHARACTER (LEN=1024) :: fname
11575 WRITE(*,*)
'Processing text files ...'
11588 IF(
mfd(i) /= 2) cycle
11590 fname(k:k)=
tfd(ia+k)
11592 WRITE(*,*)
'File ',fname(1:
lfd(i))
11593 IF (
mprint > 1)
WRITE(*,*)
' '
11594 OPEN(10,file=fname(1:
lfd(i)),iostat=ios,form=
'FORMATTED')
11596 WRITE(*,*)
'Open error for file ',fname(1:
lfd(i))
11606 READ(10,102,iostat=ierrf) text
11607 IF (ierrf < 0)
THEN
11610 WRITE(*,*)
' end-of-file after',nline,
' text lines'
11614 IF(nline <= nlinmx.AND.
mprint > 1)
THEN
11615 CALL rltext(text,ia,ib,nab)
11617 WRITE(*,101) nline,text(1:nab)
11618 IF(nline == nlinmx)
WRITE(*,*)
' ...'
11621 CALL rltext(text,ia,ib,nab)
11622 IF(ib == ia+2)
THEN
11623 mat=matint(text(ia:ib),
'end',npat,ntext)
11624 IF(mat == max(npat,ntext))
THEN
11627 WRITE(*,*)
' end-statement after',nline,
' text lines'
11633 IF(nfiln <=
nfiles)
THEN
11634 IF(nline ==
nfd(nfiln))
THEN
11649 IF(iosum /= 0)
THEN
11650 CALL peend(16,
'Aborted, open error(s) for text files')
11651 stop
'FILETX: open error(s) in text files '
11654 WRITE(*,*)
'... end of text file processing.'
11659 WRITE(*,*)
lunkno,
' unknown keywords in steering files, ', &
11660 'or file non-existing,'
11661 WRITE(*,*)
' see above!'
11662 WRITE(*,*)
'------------> stop'
11664 CALL peend(13,
'Aborted, unknown keywords in steering file')
11673 ELSE IF(
matsto == 1)
THEN
11675 ELSE IF(
matsto == 2)
THEN
11678 ELSE IF(
metsol == 1)
THEN
11680 ELSE IF(
metsol == 2)
THEN
11682 ELSE IF(
metsol == 3)
THEN
11684 ELSE IF(
metsol == 4)
THEN
11686 ELSE IF(
metsol == 5)
THEN
11688 ELSE IF(
metsol == 6)
THEN
11691 ELSE IF(
metsol == 7)
THEN
11693 ELSE IF(
metsol == 8)
THEN
11696 ELSE IF(
metsol == 9)
THEN
11701 WRITE(*,*)
'MINRES forced with sparse matrix!'
11703 WRITE(*,*)
'MINRES forced with sparse matrix!'
11705 WRITE(*,*)
'MINRES forced with sparse matrix!'
11710 WRITE(*,*)
'MINRES forced with sparse matrix!'
11712 WRITE(*,*)
'MINRES forced with sparse matrix!'
11714 WRITE(*,*)
'MINRES forced with sparse matrix!'
11722 WRITE(*,*)
'Solution method and matrix-storage mode:'
11724 WRITE(*,*)
' METSOL = 1: matrix inversion'
11725 ELSE IF(
metsol == 2)
THEN
11726 WRITE(*,*)
' METSOL = 2: diagonalization'
11727 ELSE IF(
metsol == 3)
THEN
11728 WRITE(*,*)
' METSOL = 3: decomposition'
11729 ELSE IF(
metsol == 4)
THEN
11730 WRITE(*,*)
' METSOL = 4: MINRES'
11731 ELSE IF(
metsol == 5)
THEN
11732 WRITE(*,*)
' METSOL = 5: MINRES-QLP'
11733 ELSE IF(
metsol == 6)
THEN
11734 WRITE(*,*)
' METSOL = 6: GMRES (-> MINRES)'
11736 ELSE IF(
metsol == 7)
THEN
11737 WRITE(*,*)
' METSOL = 7: LAPACK factorization'
11738 ELSE IF(
metsol == 8)
THEN
11739 WRITE(*,*)
' METSOL = 8: LAPACK factorization'
11741 ELSE IF(
metsol == 9)
THEN
11742 WRITE(*,*)
' METSOL = 9: Intel oneMKL PARDISO'
11747 WRITE(*,*)
' with',
mitera,
' iterations'
11750 WRITE(*,*)
' MATSTO = 0: unpacked symmetric matrix, ',
'n*n elements'
11751 ELSEIF(
matsto == 1)
THEN
11752 WRITE(*,*)
' MATSTO = 1: full symmetric matrix, ',
'(n*n+n)/2 elements'
11753 ELSE IF(
matsto == 2)
THEN
11754 WRITE(*,*)
' MATSTO = 2: sparse matrix (custom)'
11755 ELSE IF(
matsto == 3)
THEN
11757 WRITE(*,*)
' MATSTO = 3: sparse matrix (upper triangle, CSR3)'
11759 WRITE(*,*)
' MATSTO = 3: sparse matrix (upper triangle, BSR3)'
11763 WRITE(*,*)
' and band matrix, width',
mbandw
11767 WRITE(*,*)
'Chi square cut equiv 3 st.dev applied ...'
11768 WRITE(*,*)
' in first iteration with factor',
chicut
11769 WRITE(*,*)
' in second iteration with factor',
chirem
11770 WRITE(*,*)
' (reduced by sqrt in next iterations)'
11774 WRITE(*,*)
' Down-weighting of outliers in',
lhuber,
' iterations'
11775 WRITE(*,*)
' Cut on downweight fraction',
dwcut
11778 WRITE(*,*)
'Iterations (solutions) with line search:'
11787 WRITE(*,*)
' All with Chi square cut scaling factor <= 1.'
11818 INTEGER(mpi) :: ios
11822 INTEGER(mpi) :: npat
11823 INTEGER(mpi) :: ntext
11824 INTEGER(mpi) :: nuprae
11827 CHARACTER (LEN=*),
INTENT(INOUT) :: fname
11833 IF(len(fname) > 5)
THEN
11834 IF(fname(1:5) ==
'rfio:') nuprae=1
11835 IF(fname(1:5) ==
'dcap:') nuprae=2
11836 IF(fname(1:5) ==
'root:') nuprae=3
11838 IF(nuprae == 0)
THEN
11839 INQUIRE(file=fname,iostat=ios,exist=ex)
11840 IF(ios /= 0)
nufile=-abs(ios)
11841 IF(ios /= 0)
RETURN
11842 ELSE IF(nuprae == 1)
THEN
11855 nm=
matint(
'xt',fname(l1:ll),npat,ntext)
11858 nm=
matint(
'tx',fname(l1:ll),npat,ntext)
11879 INTEGER(mpi) :: ier
11880 INTEGER(mpi) :: iomp
11883 INTEGER(mpi) :: kkey
11884 INTEGER(mpi) :: label
11885 INTEGER(mpi) :: lkey
11886 INTEGER(mpi) :: mat
11887 INTEGER(mpi) :: miter
11888 INTEGER(mpi) :: nab
11889 INTEGER(mpi) :: nkey
11890 INTEGER(mpi) :: nkeys
11892 INTEGER(mpi) :: nmeth
11893 INTEGER(mpi) :: npat
11894 INTEGER(mpi) :: ntext
11895 INTEGER(mpi) :: nums
11896 INTEGER(mpi) :: matint
11898 CHARACTER (LEN=*),
INTENT(IN) :: text
11899 INTEGER(mpi),
INTENT(IN) :: nline
11903 parameter(nkeys=7,nmeth=10)
11905 parameter(nkeys=6,nmeth=9)
11908 parameter(nkeys=6,nmeth=7)
11910 CHARACTER (LEN=16) :: methxt(nmeth)
11911 CHARACTER (LEN=16) :: keylst(nkeys)
11912 CHARACTER (LEN=32) :: keywrd
11913 CHARACTER (LEN=32) :: keystx
11914 CHARACTER (LEN=itemCLen) :: ctext
11915 INTEGER(mpi),
PARAMETER :: mnum=100
11916 REAL(mpd) :: dnum(mnum)
11919 INTEGER(mpi) :: ipvs
11922 INTEGER(mpi) :: lpvs
11926 SUBROUTINE additem(length,list,label,value)
11928 INTEGER(mpi),
INTENT(IN OUT) :: length
11929 TYPE(listitem),
DIMENSION(:),
INTENT(IN OUT),
ALLOCATABLE :: list
11930 INTEGER(mpi),
INTENT(IN) :: label
11931 REAL(mpd),
INTENT(IN) :: value
11933 SUBROUTINE additemc(length,list,label,text)
11935 INTEGER(mpi),
INTENT(IN OUT) :: length
11936 TYPE(listitemc),
DIMENSION(:),
INTENT(IN OUT),
ALLOCATABLE :: list
11937 INTEGER(mpi),
INTENT(IN) :: label
11938 CHARACTER(LEN = itemCLen),
INTENT(IN) :: text
11940 SUBROUTINE additemi(length,list,label,ivalue)
11942 INTEGER(mpi),
INTENT(IN OUT) :: length
11943 TYPE(listitemi),
DIMENSION(:),
INTENT(IN OUT),
ALLOCATABLE :: list
11944 INTEGER(mpi),
INTENT(IN) :: label
11945 INTEGER(mpi),
INTENT(IN) :: ivalue
11952 DATA keylst/
'unknown',
'parameter',
'constraint',
'measurement',
'method',
'comment',
'pardiso'/
11953 DATA methxt/
'diagonalization',
'inversion',
'fullMINRES',
'sparseMINRES', &
11954 'fullMINRES-QLP',
'sparseMINRES-QLP',
'decomposition',
'fullLAPACK',
'unpackedLAPACK', &
11957 DATA keylst/
'unknown',
'parameter',
'constraint',
'measurement',
'method',
'comment'/
11958 DATA methxt/
'diagonalization',
'inversion',
'fullMINRES',
'sparseMINRES', &
11959 'fullMINRES-QLP',
'sparseMINRES-QLP',
'decomposition',
'fullLAPACK',
'unpackedLAPACK'/
11962 DATA keylst/
'unknown',
'parameter',
'constraint',
'measurement',
'method',
'comment'/
11963 DATA methxt/
'diagonalization',
'inversion',
'fullMINRES',
'sparseMINRES', &
11964 'fullMINRES-QLP',
'sparseMINRES-QLP',
'decomposition'/
11970 CALL rltext(text,ia,ib,nab)
11971 IF(nab == 0)
GOTO 10
11972 CALL ratext(text(1:nab),nums,dnum,mnum)
11974 IF(nums /= 0) nkey=0
11982 keystx=keylst(nkey)
11983 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
11984 IF(100*mat >= 80*max(npat,ntext))
GO TO 10
11990 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
11991 IF(100*mat >= 80*max(npat,ntext))
THEN
11993 IF(nums > 0) mprint=nint(dnum(1),mpi)
11998 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
11999 IF(100*mat >= 80*max(npat,ntext))
THEN
12002 IF(nums > 0) mdebug=nint(dnum(1),mpi)
12003 IF(nums > 1) mdebg2=nint(dnum(2),mpi)
12008 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12009 IF(100*mat >= 80*max(npat,ntext))
THEN
12010 IF(nums > 0 .AND. dnum(1) > 0.5) mreqenf=nint(dnum(1),mpi)
12011 IF(nums > 1 .AND. dnum(2) > 0.5) mreqena=nint(dnum(2),mpi)
12012 IF(nums > 2 .AND. dnum(3) > 0.5) iteren=nint(dnum(1)*dnum(3),mpi)
12016 keystx=
'printrecord'
12017 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12018 IF(100*mat >= 80*max(npat,ntext))
THEN
12019 IF(nums > 0) nrecpr=nint(dnum(1),mpi)
12020 IF(nums > 1) nrecp2=nint(dnum(2),mpi)
12025 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12026 IF(100*mat >= 80*max(npat,ntext))
THEN
12027 IF (nums > 0.AND.dnum(1) > 0.) mxrec=nint(dnum(1),mpi)
12032 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12033 IF(100*mat >= 80*max(npat,ntext))
THEN
12034 IF (nums > 0.AND.dnum(1) >= 0.) ncache=nint(dnum(1),mpi)
12035 IF (nums == 2.AND.dnum(2) > 0..AND.dnum(2) <= 1.0) &
12036 fcache(1)=real(dnum(2),mps)
12037 IF (nums >= 4)
THEN
12039 fcache(k)=real(dnum(k+1),mps)
12046 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12047 IF(100*mat >= 80*max(npat,ntext))
THEN
12052 chicut=real(dnum(1),mps)
12053 IF(chicut < 1.0) chicut=-1.0
12057 chirem=real(dnum(2),mps)
12058 IF(chirem < 1.0) chirem=1.0
12059 IF(chicut >= 1.0) chirem=min(chirem,chicut)
12067 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12068 IF(100*mat >= 80*max(npat,ntext))
THEN
12069 IF(nums > 0) chhuge=real(dnum(1),mps)
12070 IF(chhuge < 1.0) chhuge=1.0
12075 keystx=
'linesearch'
12076 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12077 IF(100*mat >= 80*max(npat,ntext))
THEN
12078 IF(nums > 0) lsearch=nint(dnum(1),mpi)
12083 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12084 IF(100*mat >= 80*max(npat,ntext))
THEN
12085 IF(nums > 0) lfitnp=nint(dnum(1),mpi)
12086 IF(nums > 1) lfitbb=nint(dnum(2),mpi)
12090 keystx=
'regularization'
12091 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12092 IF(100*mat >= 80*max(npat,ntext))
THEN
12094 regula=real(dnum(1),mps)
12095 IF(nums >= 2) regpre=real(dnum(2),mps)
12099 keystx=
'regularisation'
12100 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12101 IF(100*mat >= 80*max(npat,ntext))
THEN
12103 regula=real(dnum(1),mps)
12104 IF(nums >= 2) regpre=real(dnum(2),mps)
12109 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12110 IF(100*mat >= 80*max(npat,ntext))
THEN
12111 regpre=real(dnum(1),mps)
12116 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12117 IF(100*mat >= 80*max(npat,ntext))
THEN
12118 matrit=nint(dnum(1),mpi)
12123 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12124 IF(100*mat >= 80*max(npat,ntext))
THEN
12126 IF (nums > 0.AND.dnum(1) > 0.) matmon=nint(dnum(1),mpi)
12131 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12132 IF(100*mat >= 80*max(npat,ntext))
THEN
12133 IF(nums > 0) mbandw=nint(dnum(1),mpi)
12134 IF(mbandw < 0) mbandw=-1
12135 IF(nums > 1) lprecm=nint(dnum(2),mpi)
12159 keystx=
'outlierdownweighting'
12160 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12161 IF(100*mat >= 80*max(npat,ntext))
THEN
12162 lhuber=nint(dnum(1),mpi)
12163 IF(lhuber > 0.AND.lhuber <= 2) lhuber=2
12167 keystx=
'dwfractioncut'
12168 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12169 IF(100*mat >= 80*max(npat,ntext))
THEN
12170 dwcut=real(dnum(1),mps)
12171 IF(dwcut > 0.5) dwcut=0.5
12176 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12177 IF(100*mat >= 80*max(npat,ntext))
THEN
12178 prange=abs(real(dnum(1),mps))
12183 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12184 IF(100*mat >= 80*max(npat,ntext))
THEN
12190 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12191 IF(100*mat >= 80*max(npat,ntext))
THEN
12196 keystx=
'memorydebug'
12197 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12198 IF(100*mat >= 80*max(npat,ntext))
THEN
12200 IF (nums > 0.AND.dnum(1) > 0.0) memdbg=nint(dnum(1),mpi)
12204 keystx=
'globalcorr'
12205 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12206 IF(100*mat >= 80*max(npat,ntext))
THEN
12211 keystx=
'printcounts'
12212 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12213 IF(100*mat >= 80*max(npat,ntext))
THEN
12215 IF (nums > 0.AND.dnum(1) > 0.0) ipcntr=nint(dnum(1),mpi)
12219 keystx=
'weightedcons'
12220 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12221 IF(100*mat >= 80*max(npat,ntext))
THEN
12223 IF (nums > 0) iwcons=nint(dnum(1),mpi)
12227 keystx=
'skipemptycons'
12228 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12229 IF(100*mat >= 80*max(npat,ntext))
THEN
12234 keystx=
'resolveredundancycons'
12235 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12236 IF(100*mat >= 80*max(npat,ntext))
THEN
12241 keystx=
'withelimination'
12242 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12243 IF(100*mat >= 80*max(npat,ntext))
THEN
12249 keystx=
'withLAPACKelimination'
12250 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12251 IF(100*mat >= 80*max(npat,ntext))
THEN
12257 keystx=
'withmultipliers'
12258 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12259 IF(100*mat >= 80*max(npat,ntext))
THEN
12264 keystx=
'checkinput'
12265 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12266 IF(100*mat >= 80*max(npat,ntext))
THEN
12268 IF (nums > 0) icheck=nint(dnum(1),mpi)
12272 keystx=
'checkparametergroups'
12273 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12274 IF(100*mat >= 80*max(npat,ntext))
THEN
12279 keystx=
'monitorresiduals'
12280 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12281 IF(100*mat >= 80*max(npat,ntext))
THEN
12283 IF (nums > 0) imonit=nint(dnum(1),mpi)
12284 IF (nums > 1) measbins=max(measbins,nint(dnum(2),mpi))
12288 keystx=
'monitorpulls'
12289 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12290 IF(100*mat >= 80*max(npat,ntext))
THEN
12293 IF (nums > 0) imonit=nint(dnum(1),mpi)
12294 IF (nums > 1) measbins=max(measbins,nint(dnum(2),mpi))
12298 keystx=
'monitorprogress'
12299 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12300 IF(100*mat >= 80*max(npat,ntext))
THEN
12303 IF (nums > 0) monpg1=max(1,nint(dnum(1),mpi))
12304 IF (nums > 1) monpg2=max(1,nint(dnum(2),mpi))
12308 keystx=
'scaleerrors'
12309 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12310 IF(100*mat >= 80*max(npat,ntext))
THEN
12312 IF (nums > 0) dscerr(1:2)=dnum(1)
12313 IF (nums > 1) dscerr(2)=dnum(2)
12317 keystx=
'iterateentries'
12318 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12319 IF(100*mat >= 80*max(npat,ntext))
THEN
12320 iteren=huge(iteren)
12321 IF (nums > 0) iteren=nint(dnum(1),mpi)
12326 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12327 IF(100*mat >= 80*max(npat,ntext))
THEN
12335 WRITE(*,*)
'WARNING: multithreading not available'
12341 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12342 IF(100*mat >= 80*max(npat,ntext))
THEN
12343 WRITE(*,*)
'WARNING: keyword COMPRESS is obsolete (compression is default)'
12355 keystx=
'countrecords'
12356 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12357 IF(100*mat >= 80*max(npat,ntext))
THEN
12363 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12364 IF(100*mat >= 80*max(npat,ntext).AND.mnrsel < 100)
THEN
12365 nl=min(nums,100-mnrsel)
12367 lbmnrs(mnrsel+k)=nint(dnum(k),mpi)
12373 keystx=
'pairentries'
12374 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12375 IF(100*mat >= 80*max(npat,ntext))
THEN
12378 IF (nums > 0.AND.dnum(1) > 0.0)
THEN
12379 mreqpe=nint(dnum(1),mpi)
12380 IF (nums >= 2.AND.dnum(2) >= dnum(1)) mhispe=nint(dnum(2),mpi)
12381 IF (nums >= 3.AND.dnum(3) >= dnum(1)) msngpe=nint(dnum(3),mpi)
12387 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12388 IF(100*mat >= 80*max(npat,ntext))
THEN
12389 wolfc1=real(dnum(1),mps)
12390 wolfc2=real(dnum(2),mps)
12397 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12398 IF(100*mat >= 80*max(npat,ntext))
THEN
12400 IF (dnum(1) < 1.0e-10_mpd.OR.dnum(1) > 1.0e-04_mpd)
THEN
12401 WRITE(*,*)
'ERROR: need 1.0D-10 <= MRESTL ', &
12402 '<= 1.0D-04, but get ', dnum(1)
12411 keystx=
'mrestranscond'
12412 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12413 IF(100*mat >= 80*max(npat,ntext))
THEN
12421 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12422 IF(100*mat >= 80*max(npat,ntext))
THEN
12424 mrmode = int(dnum(1),mpi)
12429 keystx=
'nofeasiblestart'
12430 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12431 IF(100*mat >= 80*max(npat,ntext))
THEN
12437 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12438 IF(100*mat >= 80*max(npat,ntext))
THEN
12443 keystx=
'readerroraseof'
12444 mat=matint(text(ia:ib),keystx,npat,ntext)
12445 IF(100*mat >= 80*max(npat,ntext))
THEN
12451 keystx=
'LAPACKwitherrors'
12452 mat=matint(text(ia:ib),keystx,npat,ntext)
12453 IF(100*mat >= 80*max(npat,ntext))
THEN
12458 keystx=
'debugPARDISO'
12459 mat=matint(text(ia:ib),keystx,npat,ntext)
12460 IF(100*mat >= 80*max(npat,ntext))
THEN
12465 keystx=
'blocksizePARDISO'
12466 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12467 IF(100*mat >= 80*max(npat,ntext).AND.mnrsel < 100)
THEN
12468 nl=min(nums,10-mpdbsz)
12470 IF (nint(dnum(k),mpi) > 0)
THEN
12471 IF (mpdbsz == 0)
THEN
12473 ipdbsz(mpdbsz)=nint(dnum(k),mpi)
12474 ELSE IF (nint(dnum(k),mpi) > ipdbsz(mpdbsz))
THEN
12476 ipdbsz(mpdbsz)=nint(dnum(k),mpi)
12484 keystx=
'fortranfiles'
12485 mat=matint(text(ia:ib),keystx,npat,ntext)
12486 IF(mat == max(npat,ntext))
RETURN
12489 mat=matint(text(ia:ib),keystx,npat,ntext)
12490 IF(mat == max(npat,ntext))
RETURN
12492 keystx=
'closeandreopen'
12493 mat=matint(text(ia:ib),keystx,npat,ntext)
12494 IF(mat == max(npat,ntext))
RETURN
12498 IF(nums /= 0) nkey=0
12501 WRITE(*,*)
'**************************************************'
12503 WRITE(*,*)
'Unknown keyword(s): ',text(1:min(nab,50))
12505 WRITE(*,*)
'**************************************************'
1252710
IF(nkey > 0)
THEN
12531 lpvs=nint(dnum(1),mpi)
12533 CALL additem(lenparameters,listparameters,lpvs,dnum(2))
12534 CALL additem(lenpresigmas,listpresigmas,lpvs,dnum(3))
12536 WRITE(*,*)
'Line',nline,
' error, label=',lpvs
12538 ELSE IF(nums /= 0)
THEN
12540 WRITE(*,*)
'Wrong text in line',nline
12541 WRITE(*,*)
'Status: new parameter'
12542 WRITE(*,*)
'> ',text(1:nab)
12544 ELSE IF(lkey == 3)
THEN
12546 IF(nums >= 1.AND.nums <= 2)
THEN
12548 CALL additem(lenconstraints,listconstraints,lpvs,dnum(1))
12550 IF(iwcons > 0) lpvs=-2
12552 IF(nums == 2) plvs=dnum(2)
12553 CALL additem(lenconstraints,listconstraints,lpvs,plvs)
12556 WRITE(*,*)
'Wrong text in line',nline
12557 WRITE(*,*)
'Status: new keyword constraint'
12558 WRITE(*,*)
'> ',text(1:nab)
12560 ELSE IF(lkey == 4)
THEN
12562 nummeasurements=nummeasurements+1
12564 CALL additem(lenmeasurements,listmeasurements,lpvs,dnum(1))
12566 CALL additem(lenmeasurements,listmeasurements,lpvs,dnum(2))
12569 WRITE(*,*)
'Wrong text in line',nline
12570 WRITE(*,*)
'Status: new keyword measurement'
12571 WRITE(*,*)
'> ',text(1:nab)
12573 ELSE IF(lkey == 5.AND.
keyb <
keyc)
THEN
12575 IF(nums >= 1) miter=nint(dnum(1),mpi)
12576 IF(miter >= 1) mitera=miter
12577 dflim=real(dnum(2),mps)
12581 mat=matint(text(
keyb+1:
keyc),keystx,npat,ntext)
12582 IF(100*mat >= 80*max(npat,ntext))
THEN
12586 ELSE IF(i == 2)
THEN
12589 ELSE IF(i == 3)
THEN
12592 ELSE IF(i == 4)
THEN
12595 ELSE IF(i == 5)
THEN
12598 ELSE IF(i == 6)
THEN
12601 ELSE IF(i == 7)
THEN
12605 ELSE IF(i == 8)
THEN
12608 ELSE IF(i == 9)
THEN
12612 ELSE IF(i == 10)
THEN
12621 ELSE IF(nkey == 0)
THEN
12624 lpvs=nint(dnum(1),mpi)
12626 CALL additem(lenparameters,listparameters,lpvs,dnum(2))
12627 CALL additem(lenpresigmas,listpresigmas,lpvs,dnum(3))
12629 WRITE(*,*)
'Line',nline,
' error, label=',lpvs
12631 ELSE IF(nums > 1.AND.nums < 3)
THEN
12633 WRITE(*,*)
'Wrong text in line',nline
12634 WRITE(*,*)
'Status continuation parameter'
12635 WRITE(*,*)
'> ',text(1:nab)
12638 ELSE IF(lkey == 3)
THEN
12641 label=nint(dnum(i),mpi)
12642 IF(label <= 0) ier=1
12644 IF(mod(nums,2) /= 0) ier=1
12647 lpvs=nint(dnum(i),mpi)
12649 CALL additem(lenconstraints,listconstraints,lpvs,plvs)
12653 WRITE(*,*)
'Wrong text in line',nline
12654 WRITE(*,*)
'Status continuation constraint'
12655 WRITE(*,*)
'> ',text(1:nab)
12658 ELSE IF(lkey == 4)
THEN
12662 label=nint(dnum(i),mpi)
12663 IF(label <= 0) ier=1
12665 IF(mod(nums,2) /= 0) ier=1
12669 lpvs=nint(dnum(i),mpi)
12671 CALL additem(lenmeasurements,listmeasurements,lpvs,plvs)
12675 WRITE(*,*)
'Wrong text in line',nline
12676 WRITE(*,*)
'Status continuation measurement'
12677 WRITE(*,*)
'> ',text(1:nab)
12679 ELSE IF(lkey == 6)
THEN
12681 lpvs=nint(dnum(1),mpi)
12685 IF (text(j:j) ==
' ')
EXIT
12688 CALL additemc(lencomments,listcomments,lpvs,ctext)
12690 WRITE(*,*)
'Line',nline,
' error, label=',lpvs
12692 ELSE IF(nums /= 0)
THEN
12694 WRITE(*,*)
'Wrong text in line',nline
12695 WRITE(*,*)
'Status: continuation comment'
12696 WRITE(*,*)
'> ',text(1:nab)
12700 ELSE IF(lkey == 7)
THEN
12703 label=nint(dnum(i),mpi)
12704 IF(label <= 0.OR.label > 64) ier=1
12706 IF(mod(nums,2) /= 0) ier=1
12710 lpvs=nint(dnum(i),mpi)
12711 ipvs=nint(dnum(i+1),mpi)
12712 CALL additemi(lenpardiso,listpardiso,lpvs,ipvs)
12716 WRITE(*,*)
'Wrong text in line',nline
12717 WRITE(*,*)
'Status continuation measurement'
12718 WRITE(*,*)
'> ',text(1:nab)
12737 INTEGER(mpi),
INTENT(IN OUT) :: length
12738 TYPE(
listitem),
DIMENSION(:),
INTENT(IN OUT),
ALLOCATABLE :: list
12739 INTEGER(mpi),
INTENT(IN) :: label
12740 REAL(mpd),
INTENT(IN) :: value
12742 INTEGER(mpl) :: newSize
12743 INTEGER(mpl) :: oldSize
12744 TYPE(
listitem),
DIMENSION(:),
ALLOCATABLE :: tempList
12746 IF (label > 0.AND.
value == 0.)
RETURN
12747 IF (length == 0 )
THEN
12749 CALL mpalloc(list,newsize,
' list ')
12751 oldsize=
size(list,kind=
mpl)
12752 IF (length >= oldsize)
THEN
12753 newsize = oldsize + oldsize/5 + 100
12754 CALL mpalloc(templist,oldsize,
' temp. list ')
12757 CALL mpalloc(list,newsize,
' list ')
12758 list(1:oldsize)=templist(1:oldsize)
12763 list(length)%label=label
12764 list(length)%value=
value
12779 INTEGER(mpi),
INTENT(IN OUT) :: length
12780 TYPE(
listitemc),
DIMENSION(:),
INTENT(IN OUT),
ALLOCATABLE :: list
12781 INTEGER(mpi),
INTENT(IN) :: label
12782 CHARACTER(len = itemCLen),
INTENT(IN) :: text
12784 INTEGER(mpl) :: newSize
12785 INTEGER(mpl) :: oldSize
12786 TYPE(
listitemc),
DIMENSION(:),
ALLOCATABLE :: tempList
12788 IF (label > 0.AND.text ==
'')
RETURN
12789 IF (length == 0 )
THEN
12791 CALL mpalloc(list,newsize,
' list ')
12793 oldsize=
size(list,kind=
mpl)
12794 IF (length >= oldsize)
THEN
12795 newsize = oldsize + oldsize/5 + 100
12796 CALL mpalloc(templist,oldsize,
' temp. list ')
12799 CALL mpalloc(list,newsize,
' list ')
12800 list(1:oldsize)=templist(1:oldsize)
12805 list(length)%label=label
12806 list(length)%text=text
12821 INTEGER(mpi),
INTENT(IN OUT) :: length
12822 TYPE(
listitemi),
DIMENSION(:),
INTENT(IN OUT),
ALLOCATABLE :: list
12823 INTEGER(mpi),
INTENT(IN) :: label
12824 INTEGER(mpi),
INTENT(IN) :: ivalue
12826 INTEGER(mpl) :: newSize
12827 INTEGER(mpl) :: oldSize
12828 TYPE(
listitemi),
DIMENSION(:),
ALLOCATABLE :: tempList
12830 IF (length == 0 )
THEN
12832 CALL mpalloc(list,newsize,
' list ')
12834 oldsize=
size(list,kind=
mpl)
12835 IF (length >= oldsize)
THEN
12836 newsize = oldsize + oldsize/5 + 100
12837 CALL mpalloc(templist,oldsize,
' temp. list ')
12840 CALL mpalloc(list,newsize,
' list ')
12841 list(1:oldsize)=templist(1:oldsize)
12846 list(length)%label=label
12847 list(length)%ivalue=ivalue
12861 CHARACTER (LEN=*),
INTENT(IN) :: text
12862 CHARACTER (LEN=16) :: textc
12871 textl(ka:kb)=text(1:l)
12875 textc=text(1:l)//
'-end'
12883 textl(ka:kb)=textc(1:l)
12910 INTEGER(mpi),
INTENT(IN) :: lun
12911 CHARACTER (LEN=*),
INTENT(IN) :: fname
12912 CHARACTER (LEN=33) :: nafile
12913 CHARACTER (LEN=33) :: nbfile
12919 CALL peend(17,
'Aborted, file name too long')
12920 stop
'File name too long '
12923 nafile(l+1:l+1)=
'~'
12925 INQUIRE(file=nafile(1:l),exist=ex)
12927 INQUIRE(file=nafile(1:l+1),exist=ex)
12929 CALL system(
'rm '//nafile)
12932 nafile(l+1:l+1)=
' '
12933 CALL system(
'mv '//nafile//nbfile)
12935 OPEN(unit=lun,file=fname)
12946 REAL,
DIMENSION(2) :: ta
12966 IF(ncount > 1)
THEN
12968 nsecd1=int(delta,
mpi)
12970 minut1=nsecd1/60-60*nhour1
12971 secnd1=delta-60*(minut1+60*nhour1)
12973 nsecd2=int(delta,
mpi)
12975 minut2=nsecd2/60-60*nhour2
12976 secnd2=delta-60*(minut2+60*nhour2)
12977 WRITE(*,101) nhour1,minut1,secnd1, nhour2,minut2,secnd2
12982101
FORMAT(i4,
' h',i3,
' min',f5.1,
' sec total',18x,
'elapsed', &
12983 i4,
' h',i3,
' min',f5.1,
' sec')
12997 INTEGER(mpi),
INTENT(IN) :: icode
12998 CHARACTER (LEN=*),
INTENT(IN) :: cmessage
13000 CALL mvopen(9,
'millepede.end')
13001 WRITE(9,101) icode, cmessage
13002101
FORMAT(1x,i4,3x,a)
13006END SUBROUTINE peend
13018 INTEGER(mpi),
INTENT(IN) :: kfile
13019 INTEGER(mpi),
INTENT(IN) :: ithr
13020 INTEGER(mpi),
INTENT(OUT) :: ierr
13022 INTEGER(mpi),
DIMENSION(13) :: ibuff
13023 INTEGER(mpi) :: ioff
13024 INTEGER(mpi) :: ios
13026 INTEGER(mpi) :: lfn
13027 INTEGER(mpi) :: lun
13028 INTEGER(mpi) :: moddate
13029 CHARACTER (LEN=1024) :: fname
13030 CHARACTER (LEN=7) :: cfile
13035 SUBROUTINE openc(filename, lfn, lun, ios)
BIND(c)
13037 CHARACTER(kind=c_char),
DIMENSION(*),
INTENT(IN) :: filename
13038 INTEGER(c_int),
INTENT(IN),
VALUE :: lfn
13039 INTEGER(c_int),
INTENT(IN),
VALUE :: lun
13040 INTEGER(c_int),
INTENT(INOUT) :: ios
13041 END SUBROUTINE openc
13053 fname(k:k)=
tfd(ioff+k)
13058 IF(kfile <=
nfilf)
THEN
13061 OPEN(lun,file=fname(1:lfn),iostat=ios, form=
'UNFORMATTED')
13062 print *,
' lun ', lun, ios
13066 CALL openc(fname(1:lfn),lfn,lun,ios)
13068 WRITE(*,*)
'Opening of C-files not supported.'
13075 WRITE(*,*)
'Open error for file ',fname(1:lfn), ios
13076 IF (moddate /= 0)
THEN
13077 WRITE(cfile,
'(I7)') kfile
13078 CALL peend(15,
'Aborted, open error(s) for binary file ' // cfile)
13079 stop
'PEREAD: open error'
13084 ios=stat(fname(1:lfn),ibuff)
13088 WRITE(*,*)
'STAT error for file ',fname(1:lfn), ios
13092 IF (moddate /= 0)
THEN
13093 IF (ibuff(10) /= moddate)
THEN
13094 WRITE(cfile,
'(I7)') kfile
13095 CALL peend(19,
'Aborted, binary file modified (date) ' // cfile)
13096 stop
'PEREAD: file modified'
13099 yfd(kfile)=ibuff(10)
13114 INTEGER(mpi),
INTENT(IN) :: kfile
13115 INTEGER(mpi),
INTENT(IN) :: ithr
13117 INTEGER(mpi) :: lun
13121 SUBROUTINE closec(lun)
BIND(c)
13123 INTEGER(c_int),
INTENT(IN),
VALUE :: lun
13130 IF(kfile <=
nfilf)
THEN
13149 INTEGER(mpi),
INTENT(IN) :: kfile
13151 INTEGER(mpi) :: lun
13155 SUBROUTINE resetc(lun)
BIND(c)
13157 INTEGER(c_int),
INTENT(IN),
VALUE :: lun
13163 IF (kfile <=
nfilf)
THEN
13182 INTEGER(mpi) :: ipgrp
13183 INTEGER(mpi) :: irank
13184 INTEGER(mpi) :: isize
13185 INTEGER(mpi) :: ivoff
13186 INTEGER(mpi) :: itgbi
13188 INTEGER(mpi) :: msize
13189 INTEGER(mpi),
PARAMETER :: mxsize = 1000
13191 INTEGER(mpl):: length
13193 REAL(mpd),
DIMENSION(:),
ALLOCATABLE :: auxVectorD
13194 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: auxVectorI
13195 REAL(mpd),
DIMENSION(:),
ALLOCATABLE :: resParGroup
13196 REAL(mpd),
DIMENSION(:),
ALLOCATABLE :: blockParGroup
13204 IF (isize <= mxsize)
THEN
13205 msize=max(msize,isize)
13207 print *,
' CKPGRP: par. group', ipgrp,
' not checked -- too large: ', isize
13210 IF (msize == 0)
RETURN
13213 length=int(msize,mpl)*(int(msize,mpl)+1)/2
13214 CALL mpalloc(blockpargroup,length,
'(matrix) block for parameter groups (D)')
13216 CALL mpalloc(respargroup,length,
'residuals for parameter groups (D)')
13217 CALL mpalloc(auxvectori,length,
'auxiliary array (I)')
13218 CALL mpalloc(auxvectord,length,
'auxiliary array (D)')
13222 print *,
' CKPGRP par. group first label size rank'
13225 IF (isize > mxsize) cycle
13232 blockpargroup(ij)=matij(ivoff+i,ivoff+j)
13236 CALL sqminv(blockpargroup,respargroup,isize,irank, auxvectord, auxvectori)
13239 IF (isize == irank)
THEN
13243 print *,
' CKPGRP ', ipgrp,
globalparlabelindex(1,itgbi), isize, irank,
' rank deficit !!!'
13261 INTEGER(mpl) :: nan
13262 INTEGER(mpl) :: neg
13264 print *,
' Checking global matrix(D) for NANs ',
size(
globalmatd,kind=mpl)
13269 print *,
' i, nan ', i, nan
13275 print *,
' Checking diagonal elements ',
nagb
13280 print *,
' i, neg ', i, neg
13284 print *,
' CHKMAT summary ', nan, neg
13306 REAL(mpd),
INTENT(IN) :: chi2
13307 INTEGER(mpi),
INTENT(IN) :: ithrd
13308 INTEGER(mpi),
INTENT(IN) :: ndf
13309 REAL(mpd),
INTENT(IN) :: dw
13311 INTEGER(mpl) ::nadd
13339 REAL(mpd),
INTENT(OUT) ::chi2
13340 INTEGER(mpl),
INTENT(OUT) ::ndf
13341 REAL(mpd),
INTENT(OUT) ::wndf
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 pcbits(npgrp, nsparr, nsparc)
Analyze bit fields.
subroutine ndbits(npgrp, ndims, nsparr, ihst)
Analyze bit fields.
subroutine clbits(in, jreqpe, jhispe, jsngpe, jextnd, idimb, ispc)
Calculate bit (field) array size, encoding.
subroutine plbits(in, inar, inac, idimb)
Calculate bit field array size (PARDISO).
subroutine spbits(npgrp, nsparr, nsparc)
Create sparsity information.
subroutine irbits(i, j)
Fill bit fields (counters, rectangular part).
subroutine clbmap(in)
Clear (additional) bit map.
subroutine inbmap(im, jm)
Fill bit map.
subroutine ckbits(npgrp, ndims)
Check sparsity of matrix.
subroutine ggbmap(ipgrp, npair, npgrp)
Get paired (parameter) groups from map.
subroutine prbits(npgrp, nsparr)
Analyze bit fields.
subroutine gpbmap(ngroup, npgrp, npair)
Get pairs (statistic) from map.
subroutine pblbits(npgrp, ibsize, nsparr, nsparc)
Analyze bit fields.
subroutine pbsbits(npgrp, ibsize, nnzero, nblock, nbkrow)
Analyze bit fields.
subroutine inbits(im, jm, inc)
Fill bit fields (counters, triangular part).
subroutine hmplun(lunw)
unit for output
subroutine gmpdef(ig, ityp, text)
book, reset XY storage
subroutine gmpxy(ig, x, y)
add (X,Y) pair
subroutine hmpdef(ih, xa, xb, text)
book, reset histogram
subroutine gmplun(lunw)
unit for output
subroutine gmpxyd(ig, x, y, dx, dy)
add (X,Y,DX,DY)
subroutine hmpwrt(ih)
write histogram text file
subroutine gmpwrt(ig)
write XY text file
subroutine hmpldf(ih, text)
book, reset log histogram
subroutine gmprnt(ig)
print XY data
subroutine hmpent(ih, x)
entry flt.pt.
subroutine hmplnt(ih, ix)
entry integer
subroutine gmpms(ig, x, y)
mean sigma(X) from Y
subroutine hmprnt(ih)
print, content vert
subroutine monend()
End monitoring.
subroutine monini(l, n1, n2)
Initialize monitoring.
subroutine sqmibb2(v, b, n, nbdr, nbnd, inv, nrank, vbnd, vbdr, aux, vbk, vzru, scdiag, scflag)
Band bordered matrix.
subroutine dbavat(v, a, w, n, m, iopt)
A V AT product (similarity).
subroutine sqminl(v, b, n, nrank, diag, next, vk, mon)
Matrix inversion for LARGE matrices.
subroutine devsol(n, diag, u, b, x, work)
Solution by diagonalization.
subroutine dbsvxl(v, a, b, n)
Product LARGE symmetric matrix, vector.
subroutine devrot(n, diag, u, v, work, iwork)
Diagonalization.
subroutine sort22l(a, b, n)
Quick sort 2 with index.
subroutine dbavats(v, a, is, w, n, m, iopt, sc)
A V AT product (similarity, sparse).
subroutine chslv2(g, x, n)
Solve A*x=b using Cholesky decomposition.
subroutine sort1k(a, n)
Quick sort 1.
subroutine sqminv(v, b, n, nrank, diag, next)
Matrix inversion and solution.
subroutine presols(p, n, b, nm, cu, a, l, s, x, y)
Constrained (sparse) preconditioner, solution.
subroutine sqmibb(v, b, n, nbdr, nbnd, inv, nrank, vbnd, vbdr, aux, vbk, vzru, scdiag, scflag)
Bordered band matrix.
subroutine devinv(n, diag, u, v)
Inversion by diagonalization.
subroutine equdecs(n, m, b, nm, ls, c, india, l, nrkd, nrkd2)
Decomposition of (sparse) equilibrium systems.
subroutine chdec2(g, n, nrank, evmax, evmin, mon)
Cholesky decomposition (LARGE pos.
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 equslvs(n, m, b, nm, c, india, l, x)
Solution of (sparse) equilibrium systems (after decomposition).
subroutine precons(p, n, b, nm, c, cu, a, l, s, nrkd)
Constrained (sparse) preconditioner, decomposition.
subroutine sort2i(a, n)
Quick sort 2 with index.
subroutine qlpssq(aprod, B, m, t)
Partial similarity transformation by Q(t).
subroutine qldecb(a, bpar, bcon, rcon)
QL decomposition (for disjoint block matrix).
subroutine qlmlq(x, m, t)
Multiply left by Q(t) (per block).
subroutine qlsetb(ib)
Set block.
subroutine qlbsub(d, y)
Backward substitution (per block).
subroutine qlini(n, m, l, s, k)
Initialize QL decomposition.
subroutine qlgete(emin, emax)
Get eigenvalues.
subroutine qlssq(aprod, A, s, roff, t)
Similarity transformation by Q(t).
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, mnum)
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
integer, parameter mpl
long integer
integer, parameter mps
single precision
integer, parameter mpi
integer
Parameters, variables, dynamic arrays.
integer(mpl), dimension(:), allocatable csr3columnlist
list of columns for sparse matrix
integer(mpl) mszpcc
(integrated block) matrix size for constraint matrix for preconditioner
real(mpd), dimension(:), allocatable workspaceeigenvectors
workspace eigen vectors
real(mpd), dimension(:), allocatable globalparameter
global parameters (start values + sum(x_i))
integer(mpl) nrecal
number of records
integer(mpi), dimension(:), allocatable localglobalmap
matrix correlating local and global par, map (counts)
type(listitem), dimension(:), allocatable listparameters
list of parameters from steering file
integer(mpi), dimension(:), allocatable vecparblockconoffsets
global par block (constraint) offsets
real(mpd), dimension(:), allocatable lapacktau
LAPACK TAU (QL decomp.)
integer(mpl) mszprd
(integrated block) matrix size for (constraint) product matrix
integer(mpi) lunmon
unit for monitoring output file
real(mpd), dimension(:), allocatable vecconsresiduals
residuals of constraints
integer(mpl) nrec1
record number with largest residual
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)
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), dimension(10) ipdbsz
PARDISO, list of block sizes to be tried (by PBSBITS)
integer(mpi) metsol
solution method (1: inversion, 2: diagonalization, 3: decomposition, 4: MINRES, 5: 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(mpl) sumrecords
sum of records
integer(mpi) mreqpe
min number of pair entries
integer(mpi) memdbg
debug flag for memory management
integer(mpi), dimension(100) lbmnrs
MINRES error labels.
integer(mpi) ncgrp
number of (disjoint) constraint groups
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.
integer(mpi) matsto
(global) matrix storage mode (0: unpacked, 1: full = packed, 2: sparse(custom), 3: sparse(CSR3,...
integer(mpi), dimension(:,:), allocatable matconssort
keys and index for sorting
real(mpd), dimension(:), allocatable lapackwork
LAPACK work array.
integer(mpi) monpg1
progress monitoring, repetition rate start value
integer(mpi), dimension(:,:), allocatable readbufferinfo
buffer management (per thread)
integer(mpi) nhistp
flag for histogram printout
integer(mpl), dimension(:), allocatable csr3rowoffsets
row offsets for column list
real(mpd), dimension(:), allocatable globalparcopy
copy of global parameters
real(mpd), dimension(:), allocatable lapackql
LAPACK QL (QL decomp.)
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), dimension(:,:), allocatable matconsranges
parameter ranges for constraints
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) mdebug
debug flag (number of records to print)
integer(mpi) npblck
number of (disjoint) parameter blocks (>1: block diagonal storage)
real(mpd), dimension(:), allocatable matconsproduct
product matrix of constraints
integer(mpi), dimension(:), allocatable yfd
binary file: modification date
integer(mpi) nxlow
(max of) global parameters with too few accepted entries for icalcm=1
integer(mpl) ndgb
number of global derivatives read
real(mps) value1
largest residual
integer(mpi) ipddbg
flag for debugging Intel oneMKL PARDISO
real(mpd), dimension(:), allocatable localcorrections
local fit corrections (to residuals)
integer(mpl) nrec3
(1.) record number with error
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, content
integer(mpi) mhispe
upper bound for pair entry histogrammimg
integer(mpi) nfgb
number of fit parameters
integer(mpi), dimension(:,:), allocatable kfd
(1,.)= number of records in file, (2,..)= file order
real(mpd), dimension(:), allocatable globalchi2sumd
fractional part of Chi2 sum
integer(mpi) icheck
flag for checking input only (no solution determined)
integer(mpi), dimension(:), allocatable jfd
file: number of accepted records
integer(mpl) nrecer
record with error (rank deficit or Not-a-Number) for printout
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) nfilw
number of weighted binary files
integer(mpi) ndefpg
number of parameter groups with rank deficit (from inversion)
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(mpl) nrec2
record number with largest chi^2/Ndf
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 (var.
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
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(mpl), dimension(:), allocatable lapackipiv
LAPACK IPIV (pivot)
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), dimension(:), allocatable globalindexranges
global par ranges
integer(mpi) mthrd
number of (OpenMP) threads
integer(mpi) mbandw
band width of preconditioner matrix
integer(mpl) lplwrk
length of LAPACK WORK array
real(mps) dwcut
down-weight fraction cut
integer(mpl), dimension(:), allocatable globalcounter
global counter (entries in 'x')
real(mps), dimension(:), allocatable globalmatf
global matrix 'A' (float part for compressed sparse)
integer(mpi), dimension(:,:), allocatable matconsgroups
start of constraint groups, parameter range
real(mps), dimension(0:8) times
cpu time counters
integer(mpi) minrecordsinblock
min.
integer(mpi), dimension(:), allocatable localglobalstructure
matrix correlating local and global par, (sparsity) structure
real(mpd), dimension(:), allocatable globalndfsumw
weighted NDF sum
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
integer(mpi), dimension(:,:), allocatable globaltotindexgroups
integer(mpi), dimension(:), allocatable vecconsgroupcounts
counter for constraint groups
real(mps) deltim
cpu time difference
integer(mpi) igcorr
flag for output of global correlations for inversion, =0: none
integer(mpi), dimension(-8:0) globalparheader
global parameters (mapping) header
integer(mpi) lencomments
length of list of (global parameter) comments from steering file
integer(mpl), dimension(:), allocatable offprecond
preconditioner (block matrix) offsets
real(mpd), dimension(:), allocatable vecconssolution
solution for constraint elimination
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(mpl) negb
number of equations read with global parameters
integer(mpi) keepopen
flag for keeping binary files open
real(mpd), dimension(:), allocatable workspacediagonalization
workspace diag.
real(mps), dimension(:), allocatable wfd
binary file: weight
real(mpd), dimension(:), allocatable matprecond
preconditioner matrix (band and other parts)
integer(mpi) ntgb
total number of global parameters
real(mps) angras
angle between gradient and search direction
type(listitemc), dimension(:), allocatable listcomments
list of comments from steering file
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 (non overlapping) 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), dimension(:), allocatable workspacerow
(pivot) row of global matrix (for global corr.)
integer(mpi) maxrecordsinblock
max.
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) lenpardiso
length of list of Intel oneMKL PARDISO parameters (indices 1..64)
integer(mpi) nbdrx
max border size for local fit
integer(mpi), dimension(:,:), allocatable globalparlabelindex
global parameters label, total -> var.
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.
integer(mpi), dimension(:,:), allocatable sfd
offset (1,..), length (2,..) of binary file name in tfd
integer(mpi) lenconstraints
length of list of constraints from steering file
integer(mpi), dimension(:), allocatable blockprecond
preconditioner (constraint) blocks
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(mpl) nrecp2
record number with printout
integer(mpl) nrec
number of records read
integer(mpi), dimension(:,:), allocatable matparblockoffsets
global par block offsets (parameter, constraint blocks)
integer(mpl) nrecpr
record number with printout
integer(mpl), dimension(:), allocatable ifd
file: integrated record numbers (=offset)
integer(mpi) nofeas
flag for skipping making parameters feasible
integer(mpi) matbsz
(global) matrix (fixed) block size, only used for BSR3 storage mode (Intel oneMKL PARDISO)
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
integer(mpi) ntpgrp
number of parameter groups
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), dimension(3) nprecond
number of constraints (blocks), matrix size for preconditioner
integer(mpi) nalcn
max number of local paramters per record
integer(mpi), dimension(:), allocatable globalparcomments
global parameters comments
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
integer(mpi) mcount
flag for grouping and counting global parameters on equlation (0) or record (1) level
real(mps), dimension(:), allocatable ofd
file: option
integer(mpi) ireeof
flag for treating (binary file) read errors as end-of-file
integer(mpi) ifile
current file (index)
real(mps) delfun
expected function change
integer(mpi) iitera
MINRES iterations.
integer(mpl) skippedrecords
number of skipped records (buffer too small)
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) napgrp
number of all parameter groups (variable + Lagrange mult.)
integer(mpl) nrecd
number of records read containing doubles
integer(mpi), dimension(:,:), allocatable localequations
indices (ISJAJB) for local equations (measurements)
integer(mpi), dimension(:), allocatable globalallpartogroup
all parameters variable -> group index
integer(mpi), dimension(:), allocatable backindexusage
list of global par in record
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) mpdbsz
PARDISO, number of block sizes to be tried (by PBSBITS)
integer(mpi) lfitnp
local fit: number of iteration to calculate pulls
integer(mpl), dimension(:), allocatable globalparlabelcounter
global parameters label counters
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
integer(mpi) irslvrc
flag for resolving redundancy constraints (two equivalent parameter groups)
real(mpd), dimension(:), allocatable workspaced
(general) workspace (D)
integer(mpl) neqn
number of equations (measurements) read
integer(mpi) measbins
number of bins per measurement for monitoring
integer(mpl) mszcon
(integrated block) matrix size for constraint matrix
integer(mpi), dimension(:), allocatable nfd
index (line) in (steering) file
integer(mpi) ilperr
flag to calculate parameter errors with LAPACK
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), dimension(:), allocatable measindex
mapping of 1.
integer(mpi) istopa
MINRES istop (convergence)
integer(mpi), dimension(:), allocatable mfd
file mode: cbinary =1, text =2, fbinary=3
real(mpd), dimension(:), allocatable blvec
local fit vector 'b' (in A*x=b), replaced by 'x'
logical newite
flag for new iteration
integer(mpi) nrderr
number of binary files with read errors
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) monpg2
progress monitoring, repetition rate max increase
integer(mpl), dimension(:), allocatable globalrowoffsets
row offsets for full or unpacked matrix
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) nvpgrp
number of variable parameter groups
integer(mpi), dimension(:), allocatable xfd
file: max.
integer(mpi) mreqena
required number of entries (for variable global parameter from Accepted local fits)
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(mpl) ndfsum
sum(ndf)
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(mpl), dimension(:), allocatable globalchi2sumi
integer part of Chi2 sum
integer(mpl) ipdmem
memory (kB) used by Intel oneMKL PARDISO
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(mpl), dimension(:), allocatable globalndfsum
NDF sum.
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
type(listitemi), dimension(:), allocatable listpardiso
list of Intel oneMKL PARDISO parameters
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
integer(mpi) ichkpg
flag for checking (rank of) parameter groups
type(listitem), dimension(:), allocatable listpresigmas
list of pre-sgmas from steering file
integer(mpi), dimension(:), allocatable globalallindexgroups
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.
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 first keyword
integer(mpi) keya
start (position) of first keyword
integer(mpi) keyc
end (position) of last keyword
subroutine ploopb(lunp)
Print iteration line.
subroutine mchdec
Solution by Cholesky decomposition.
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 mgupdt(i, j1, j2, il, jl, n, sub)
Update global matrix for parameter group.
subroutine lpavat(t)
Similarity transformation by Q(t).
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 pargrp(inds, inde)
Parameter group info update for block of parameters.
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 (custom).
subroutine mdiags
Solution by diagonalization.
program mptwo
Millepede II main program Pede.
subroutine prtstat
Print input statistic.
real(mpd) function matij(itema, itemb)
Get matrix element at (i,j).
subroutine grpcon
Group constraints.
subroutine loopbf(nrej, numfil, naccf, chi2f, ndff)
Loop over records in read buffer (block), fits and sums.
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.
integer(mpl) function ijcsr3(itema, itemb)
Index for sparse storage (CSR3).
subroutine useone
Make usable (sort items and redefine hash indices).
subroutine mvopen(lun, fname)
Open file.
subroutine chkrej
Check rejection details.
subroutine avprd0(n, l, x, b)
Product symmetric (sub block) matrix times vector.
subroutine addsums(ithrd, chi2, ndf, dw)
Accurate summation.
subroutine solgloqlp(ivgbi)
Error for single global parameter from MINRES-QLP.
subroutine lpqldec(a, emin, emax)
QL decomposition.
subroutine addcst
Add constraint information to matrix and vector.
subroutine petime
Print times.
subroutine mstart(text)
Start of 'module' printout.
subroutine mend
End of 'module' printout.
subroutine anasps
Analyse sparsity structure.
subroutine minver
Solution by matrix inversion.
subroutine peprep(mode)
Prepare records.
integer(mpi) function ijprec(itema, itemb)
Precision for storage of parameter groups.
subroutine explfc(lunit)
Print explanation of iteration table.
subroutine getsums(chi2, ndf, wndf)
Get accurate sums.
subroutine chkmat
Check global matrix.
subroutine binopn(kfile, ithr, ierr)
Open binary file.
subroutine pepgrp
Parameter group info update.
subroutine sechms(deltat, nhour, minut, secnd)
Time conversion.
integer(mpi) function inone(item)
Translate labels to indices (for global parameters).
subroutine avprds(n, l, x, is, ie, b)
Product symmetric (sub block) matrix times sparse vector.
subroutine avprod(n, x, b)
Product symmetric matrix times vector.
subroutine ijpgrp(itema, itemb, ij, lr, iprc)
Index (region length and precision) for sparse storage of parameter groups.
subroutine loop1i
Iteration of first data loop.
subroutine mhalf2
Fill 2nd half of matrix for extended storage.
subroutine ckpgrp
Check (rank of) parameter groups.
subroutine additemi(length, list, label, ivalue)
add item to list
subroutine mminrsqlp
Solution with MINRES-QLP.
subroutine filetc
Interprete command line option, steering file.
subroutine feasib(concut, iact)
Make parameters feasible.
subroutine mspardiso
Solution with Intel(R) oneAPI Math Kernel Library (oneMKL) PARDISO.
subroutine mdutrf
Solution by factorization.
subroutine mdptrf
Solution by factorization.
subroutine mvsolv(n, x, y)
Solution for finite band width preconditioner.
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.
subroutine additemc(length, list, label, text)
add character item to list
void resetc(int nFileIn)
Rewind file.
void initc(int nFiles)
Initialises the 'global' variables used for file handling.
void closec(int nFileIn)
Close file.
void readc(double *bufferDouble, float *bufferFloat, int *bufferInt, int *lengthBuffers, int nFileIn, int *errorFlag)
Read record from file.
void openc(const char *fileName, int lfn, int nFileIn, int *errorFlag)
Open file.
list items from steering file
character list items from steering file
integer list items from steering file