896 REAL,
DIMENSION(2) :: ta
899 INTEGER(mpi) :: iopnmp
910 INTEGER(mpi) :: nsecnd
911 INTEGER(mpi) :: ntsec
913 CHARACTER (LEN=24) :: chdate
914 CHARACTER (LEN=24) :: chost
916 CHARACTER (LEN=6) :: c6
917 INTEGER major, minor, patch
941 CALL getenv(
'HOSTNAME',chost)
942 IF (chost(1:1) ==
' ')
CALL getenv(
'HOST',chost)
943 WRITE(*,*)
'($Id: c6be33cb87ef25e4226d5088ec53ebf6d39f4712 $)'
948 CALL ilaver( major,minor, patch )
949 WRITE(*,110) lapack64, major,minor, patch
950110
FORMAT(
' using LAPACK64 with ',(a),
', version ',i0,
'.',i0,
'.',i0)
952 WRITE(*,*)
'using Intel oneMKL PARDISO'
956 WRITE(*,111) __gnuc__ , __gnuc_minor__ , __gnuc_patchlevel__
957111
FORMAT(
' compiled with gcc ',i0,
'.',i0,
'.',i0)
960 WRITE(*,111) __pgic__ , __pgic_minor__ , __pgic_patchlevel__
961111
FORMAT(
' compiled with pgi ',i0,
'.',i0,
'.',i0)
964 WRITE(*,*)
' < Millepede II-P starting ... ',chdate
968 WRITE(8,*)
'($Id: c6be33cb87ef25e4226d5088ec53ebf6d39f4712 $)'
970 WRITE(8,*)
'Log-file Millepede II-P ', chdate
971 WRITE(8,*)
' ', chost
973 CALL peend(-1,
'Still running or crashed')
982 WRITE(*,*)
'!!! Checking input only, no calculation of a solution !!!'
983 WRITE(8,*)
'!!! Checking input only, no calculation of a solution !!!'
1002 CALL getenv(
'OMP_NUM_THREADS',c6)
1004 CALL getenv(lapack64//
'_NUM_THREADS',c6)
1006 IF (c6(1:1) ==
' ')
THEN
1008 WRITE(*,*)
'Number of threads for LAPACK: unkown (empty OMP_NUM_THREADS)'
1010 WRITE(*,*)
'Number of threads for LAPACK: unkown (empty ',lapack64//
'_NUM_THREADS)'
1013 WRITE(*,*)
'Number of threads for LAPACK: ', c6
1033 CALL mvopen(lun,
'millepede.his')
1039 CALL mvopen(1,
'mpdebug.txt')
1043 times(0)=rstext-rstp
1049 times(1)=rloop1-rstext
1053 WRITE(8,*)
'Chi square cut equiv 3 st.dev applied ...'
1054 WRITE(8,*)
' in first iteration with factor',
chicut
1055 WRITE(8,*)
' in second iteration with factor',
chirem
1056 WRITE(8,*)
' (reduced by sqrt in next iterations)'
1060 WRITE(8,*)
'Down-weighting of outliers in',
lhuber,
' iterations'
1061 WRITE(8,*)
'Cut on downweight fraction',
dwcut
1065 times(2)=rloop2-rloop1
1070 CALL peend(5,
'Ended without solution (empty constraints)')
1072 CALL peend(0,
'Ended normally')
1109 CALL gmpdef(6,1,
'log10(#records) vs file number')
1110 CALL gmpdef(7,1,
'final rejection fraction vs file number')
1112 'final <Chi^2/Ndf> from accepted local fits vs file number')
1113 CALL gmpdef(9,1,
'<Ndf> from accepted local fits vs file number')
1119 rej=real(nrc-
jfd(kfl),mps)/real(nrc,mps)
1120 CALL gmpxy(6,real(kfl,mps),log10(real(nrc,mps)))
1121 CALL gmpxy(7,real(kfl,mps),rej)
1123 IF (
jfd(kfl) > 0)
THEN
1124 c2ndf=
cfd(kfl)/real(
jfd(kfl),mps)
1125 CALL gmpxy(8,real(kfl,mps),c2ndf)
1126 andf=real(
dfd(kfl),mps)/real(
jfd(kfl),mps)
1127 CALL gmpxy(9,real(kfl,mps),andf)
1144 WRITE(*,*)
'Misalignment test wire chamber'
1147 CALL hmpdef( 9,-0.0015,+0.0015,
'True - fitted displacement')
1148 CALL hmpdef(10,-0.0015,+0.0015,
'True - fitted Vdrift')
1154 sums(1)=sums(1)+diff
1155 sums(2)=sums(2)+diff*diff
1157 sums(3)=sums(3)+diff
1158 sums(4)=sums(4)+diff*diff
1160 sums(1)=0.01_mpd*sums(1)
1161 sums(2)=sqrt(0.01_mpd*sums(2))
1162 sums(3)=0.01_mpd*sums(3)
1163 sums(4)=sqrt(0.01_mpd*sums(4))
1164 WRITE(*,143)
'Parameters 1 - 100: mean =',sums(1),
'rms =',sums(2)
1165 WRITE(*,143)
'Parameters 101 - 200: mean =',sums(3),
'rms =',sums(4)
1166143
FORMAT(6x,a28,f9.6,3x,a5,f9.6)
1169 WRITE(*,*)
' I label simulated fitted diff'
1170 WRITE(*,*)
' -------------------------------------------- '
1190 WRITE(*,*)
'Misalignment test Si tracker'
1193 CALL hmpdef( 9,-0.0025,+0.0025,
'True - fitted displacement X')
1194 CALL hmpdef(10,-0.025,+0.025,
'True - fitted displacement Y')
1205 sums(1)=sums(1)+1.0_mpd
1206 sums(2)=sums(2)+diff
1207 sums(3)=sums(3)+diff*diff
1212 err=sqrt(abs(gmati))
1214 sums(7)=sums(7)+1.0_mpd
1215 sums(8)=sums(8)+diff
1216 sums(9)=sums(9)+diff*diff
1219 IF (mod(i,3) == 1)
THEN
1223 sums(4)=sums(4)+1.0_mpd
1224 sums(5)=sums(5)+diff
1225 sums(6)=sums(6)+diff*diff
1230 err=sqrt(abs(gmati))
1232 sums(7)=sums(7)+1.0_mpd
1233 sums(8)=sums(8)+diff
1234 sums(9)=sums(9)+diff*diff
1239 sums(2)=sums(2)/sums(1)
1240 sums(3)=sqrt(sums(3)/sums(1))
1241 sums(5)=sums(5)/sums(4)
1242 sums(6)=sqrt(sums(6)/sums(4))
1243 WRITE(*,143)
'Parameters 1 - 500: mean =',sums(2),
'rms =',sums(3)
1244 WRITE(*,143)
'Parameters 501 - 700: mean =',sums(5),
'rms =',sums(6)
1245 IF (sums(7) > 0.5_mpd)
THEN
1246 sums(8)=sums(8)/sums(7)
1247 sums(9)=sqrt(sums(9)/sums(7))
1248 WRITE(*,143)
'Parameter pulls, all: mean =',sums(8),
'rms =',sums(9)
1252 WRITE(*,*)
' I label simulated fitted diff'
1253 WRITE(*,*)
' -------------------------------------------- '
1265 IF (mod(i,3) == 1)
THEN
1285 WRITE(8,*)
'Record',
nrec1,
' has largest residual:',
value1
1288 WRITE(8,*)
'Record',
nrec2,
' has largest Chi^2/Ndf:',
value2
1292 WRITE(8,*)
'Record',
nrec3,
' is first with error (rank deficit/NaN)'
1296 WRITE(8,*)
'In total 3 +',
nloopn,
' loops through the data files'
1298 WRITE(8,*)
'In total 2 +',
nloopn,
' loops through the data files'
1302 WRITE(8,*)
'In total ',
mnrsit,
' internal MINRES iterations'
1310 ntsec=nint(deltat,mpi)
1311 CALL sechms(deltat,nhour,minut,secnd)
1312 nsecnd=nint(secnd,mpi)
1313 WRITE(8,*)
'Total time =',ntsec,
' seconds =',nhour,
' h',minut, &
1314 ' m',nsecnd,
' seconds'
1316 WRITE(8,*)
'end ', chdate
1325 WRITE(8,*)
'Data rejected in last iteration: '
1327 nrejec(0),
' (rank deficit/NaN) ',
nrejec(1),
' (Ndf=0) ', &
1334 WRITE(*,*)
' < Millepede II-P ending ... ', chdate
1341106
FORMAT(
' PARDISO dyn. memory allocation: ',f11.6,
' GB')
1346102
FORMAT(2x,i4,i10,2x,3f10.5)
1347103
FORMAT(
' Times [in sec] for text processing',f12.3/ &
1350 ' func. value ',f12.3,
' *',f4.0/ &
1351 ' func. value, global matrix, solution',f12.3,
' *',f4.0/ &
1352 ' new solution',f12.3,
' *',f4.0/)
1353105
FORMAT(
' Peak dynamic memory allocation: ',f11.6,
' GB')
1373 INTEGER(mpi) :: istop
1374 INTEGER(mpi) :: itgbi
1375 INTEGER(mpi) :: itgbl
1377 INTEGER(mpi) :: itnlim
1378 INTEGER(mpi) :: nout
1380 INTEGER(mpi),
INTENT(IN) :: ivgbi
1417 globalcorrections, itnlim, nout, rtol, istop, itn, anorm, acond, rnorm, arnorm, ynorm)
1421 globalcorrections, itnlim, nout, rtol, istop, itn, anorm, acond, rnorm, arnorm, ynorm)
1424 globalcorrections, itnlim, nout, rtol, istop, itn, anorm, acond, rnorm, arnorm, ynorm)
1430 err=sqrt(abs(real(gmati,mps)))
1431 IF(gmati < 0.0_mpd) err=-err
1432 diag=matij(ivgbi,ivgbi)
1433 gcor2=real(1.0_mpd-1.0_mpd/(gmati*diag),mps)
1435101
FORMAT(1x,
' label parameter presigma differ', &
1436 ' Error gcor^2 iit'/ 1x,
'---------',2x,5(
'-----------'),2x,
'----')
1437102
FORMAT(i10,2x,4g12.4,f7.4,i6,i4)
1457 INTEGER(mpi) :: istop
1458 INTEGER(mpi) :: itgbi
1459 INTEGER(mpi) :: itgbl
1461 INTEGER(mpi) :: itnlim
1462 INTEGER(mpi) :: nout
1464 INTEGER(mpi),
INTENT(IN) :: ivgbi
1493 mxxnrm = real(
nagb,mpd)/sqrt(epsilon(mxxnrm))
1495 trcond = 1.0_mpd/epsilon(trcond)
1496 ELSE IF(
mrmode == 2)
THEN
1504 itnlim=itnlim, rtol=rtol, maxxnorm=mxxnrm, trancond=trcond, &
1508 itnlim=itnlim, rtol=rtol, maxxnorm=mxxnrm, trancond=trcond, &
1512 itnlim=itnlim, rtol=rtol, maxxnorm=mxxnrm, trancond=trcond, &
1519 err=sqrt(abs(real(gmati,mps)))
1520 IF(gmati < 0.0_mpd) err=-err
1521 diag=matij(ivgbi,ivgbi)
1522 gcor2=real(1.0_mpd-1.0_mpd/(gmati*diag),mps)
1524101
FORMAT(1x,
' label parameter presigma differ', &
1525 ' Error gcor^2 iit'/ 1x,
'---------',2x,5(
'-----------'),2x,
'----')
1526102
FORMAT(i10,2x,4g12.4,f7.4,i6,i4)
1539 INTEGER(mpi) :: icgb
1540 INTEGER(mpi) :: irhs
1541 INTEGER(mpi) :: itgbi
1542 INTEGER(mpi) :: ivgb
1544 INTEGER(mpi) :: jcgb
1546 INTEGER(mpi) :: label
1548 INTEGER(mpi) :: inone
1551 REAL(mpd) :: drhs(4)
1552 INTEGER(mpi) :: idrh (4)
1577 IF(abs(rhs) > climit)
THEN
1583 WRITE(*,101) (idrh(l),drhs(l),l=1,irhs)
1592 WRITE(*,101) (idrh(l),drhs(l),l=1,irhs)
1595 WRITE(*,102)
' Constraints: only equation values >', climit,
' are printed'
1596101
FORMAT(
' ',4(i6,g11.3))
1597102
FORMAT(a,g11.2,a)
1610 INTEGER(mpi) :: icgb
1611 INTEGER(mpi) :: icgrp
1612 INTEGER(mpi) :: ioff
1613 INTEGER(mpi) :: itgbi
1615 INTEGER(mpi) :: jcgb
1616 INTEGER(mpi) :: label
1617 INTEGER(mpi) :: labelf
1618 INTEGER(mpi) :: labell
1619 INTEGER(mpi) :: last
1620 INTEGER(mpi) :: line1
1621 INTEGER(mpi) :: ncon
1622 INTEGER(mpi) :: ndiff
1623 INTEGER(mpi) :: npar
1624 INTEGER(mpi) :: inone
1625 INTEGER(mpi) :: itype
1626 INTEGER(mpi) :: ncgbd
1627 INTEGER(mpi) :: ncgbr
1628 INTEGER(mpi) :: ncgbw
1629 INTEGER(mpi) :: ncgrpd
1630 INTEGER(mpi) :: ncgrpr
1631 INTEGER(mpi) :: next
1633 INTEGER(mpl):: length
1634 INTEGER(mpl) :: rows
1636 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: vecParConsOffsets
1637 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: vecParConsList
1638 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: vecConsParOffsets
1639 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: vecConsParList
1640 INTEGER(mpi),
DIMENSION(:,:),
ALLOCATABLE :: matConsGroupIndex
1653 IF(last < 0.AND.label < 0)
THEN
1656 IF(itype == 2) ncgbw=ncgbw+1
1663 IF(label > 0.AND.itype == 2)
THEN
1670 IF (ncgbw == 0)
THEN
1671 WRITE(*,*)
'GRPCON:',
ncgb,
' constraints found in steering files'
1673 WRITE(*,*)
'GRPCON:',
ncgb,
' constraints found in steering files,',ncgbw,
'weighted'
1678 length=
ncgb+1; rows=3
1691 CALL mpalloc(matconsgroupindex,rows,length,
'group index for constraint (I)')
1694 CALL mpalloc(vecconsparoffsets,length,
'offsets for global par list for cons. (I)')
1696 CALL mpalloc(vecparconsoffsets,length,
'offsets for cons. list for global par. (I)')
1697 vecparconsoffsets(1)=0
1699 vecparconsoffsets(i+1)=vecparconsoffsets(i)+
globalparcons(i)
1703 length=vecparconsoffsets(
ntgb+1)
1704 CALL mpalloc(vecconsparlist,length,
'global par. list for constraint (I)')
1705 CALL mpalloc(vecparconslist,length,
'constraint list for global par. (I)')
1710 vecconsparoffsets(1)=ioff
1722 vecparconslist(vecparconsoffsets(itgbi)+
globalparcons(itgbi))=icgb
1724 vecconsparlist(ioff+npar)=itgbi
1730 CALL sort1k(vecconsparlist(ioff+1),npar)
1734 next=vecconsparlist(ioff+j)
1735 IF (next /= last)
THEN
1737 vecconsparlist(ioff+ndiff) = next
1746 vecconsparoffsets(icgb+1)=ioff
1759 print *,
' Constraint #parameters first par. last par. first line'
1767 npar=vecconsparoffsets(icgb+1)-vecconsparoffsets(icgb)
1771 print *, jcgb, npar, labelf, labell, line1
1774 icgrp=matconsgroupindex(1,icgb)
1775 IF (icgrp == 0)
THEN
1777 DO i=vecconsparoffsets(icgb)+1, vecconsparoffsets(icgb+1)
1778 itgbi=vecconsparlist(i)
1780 DO j=vecparconsoffsets(itgbi)+1,vecparconsoffsets(itgbi+1)
1781 icgrp=matconsgroupindex(1,vecparconslist(j))
1787 IF (icgrp == 0)
THEN
1794 matconsgroupindex(2,icgb)=jcgb
1795 matconsgroupindex(3,icgb)=icgb
1796 DO i=vecconsparoffsets(icgb)+1, vecconsparoffsets(icgb+1)
1797 itgbi=vecconsparlist(i)
1800 DO j=vecparconsoffsets(itgbi)+1,vecparconsoffsets(itgbi+1)
1801 matconsgroupindex(1,vecparconslist(j))=icgrp
1805 WRITE(*,*)
'GRPCON:',
ncgrp,
' disjoint constraints groups built'
1813 icgb=matconsgroupindex(3,jcgb)
1818 icgrp=matconsgroupindex(1,jcgb)
1838 print *,
' cons.group first con. first par. last par. #cons #par'
1849 print *, icgrp,
matconsgroups(1,icgrp), labelf, labell, ncon, npar
1852 IF (ncon == npar)
THEN
1859 print *, icgrp,
matconsgroups(1,icgrp), labelf, labell,
' : cons.group resolved'
1874 print *, icgrp,
matconsgroups(1,icgrp), labelf, labell,
' : cons.group redundant'
1879 IF (ncgrpr > 0)
THEN
1880 WRITE(*,*)
'GRPCON:',ncgbr,
' redundancy constraints in ', ncgrpr,
' groups resolved'
1884 IF (ncgrpd > 0)
THEN
1885 WRITE(*,*)
'GRPCON:',ncgbd,
' redundancy constraints in ', ncgrpd,
' groups detected'
1908 INTEGER(mpi) :: icgb
1909 INTEGER(mpi) :: icgrp
1910 INTEGER(mpi) :: ifrst
1911 INTEGER(mpi) :: ilast
1912 INTEGER(mpi) :: isblck
1913 INTEGER(mpi) :: itgbi
1914 INTEGER(mpi) :: ivgb
1916 INTEGER(mpi) :: jcgb
1917 INTEGER(mpi) :: jfrst
1918 INTEGER(mpi) :: label
1919 INTEGER(mpi) :: labelf
1920 INTEGER(mpi) :: labell
1921 INTEGER(mpi) :: ncon
1922 INTEGER(mpi) :: ngrp
1923 INTEGER(mpi) :: npar
1924 INTEGER(mpi) :: ncnmxb
1925 INTEGER(mpi) :: ncnmxg
1926 INTEGER(mpi) :: nprmxb
1927 INTEGER(mpi) :: nprmxg
1928 INTEGER(mpi) :: inone
1929 INTEGER(mpi) :: nvar
1931 INTEGER(mpl):: length
1932 INTEGER(mpl) :: rows
1934 INTEGER(mpi),
DIMENSION(:,:),
ALLOCATABLE :: matConsGroupIndex
1947 length=
ncgrp+1; rows=3
1952 CALL mpalloc(matconsgroupindex,rows,length,
'group index for constraint (I)')
1989 IF (nvar > 0 .OR.
iskpec == 0)
THEN
1992 matconsgroupindex(1,
ncgb)=ngrp+1
1993 matconsgroupindex(2,
ncgb)=icgb
1994 matconsgroupindex(3,
ncgb)=nvar
1997 IF (
ncgb > ncon) ngrp=ngrp+1
2003 WRITE(*,*)
'PRPCON:',
ncgbe,
' empty constraints skipped'
2005 WRITE(*,*)
'PRPCON:',
ncgbe,
' empty constraints detected, to be fixed !!!'
2006 WRITE(*,*)
' (use option "skipemptycons" to skip those)'
2011 WRITE(*,*)
'!!! Switch to "-C" (checking input only), no calculation of a solution !!!'
2012 WRITE(8,*)
'!!! Switch to "-C" (checking input only), no calculation of a solution !!!'
2017 WRITE(*,*)
'PRPCON:',
ncgb,
' constraints accepted'
2020 IF(
ncgb == 0)
RETURN
2027 icgb=matconsgroupindex(2,jcgb)
2032 icgrp=matconsgroupindex(1,jcgb)
2058 WRITE(*,*)
' Cons. sorted index #var.par. first line first label last label'
2059 WRITE(*,*)
' Cons. group index first cons. last cons. first label last label'
2060 WRITE(*,*)
' Cons. block index first group last group first label last label'
2066 nvar=matconsgroupindex(3,jcgb)
2070 WRITE(*,*)
' Cons. sorted', jcgb, nvar, &
2073 WRITE(*,*)
' Cons. sorted', jcgb,
' empty (0)', &
2082 WRITE(*,*)
' Cons. group ', icgrp,
matconsgroups(1,icgrp), &
2095 DO i=icgrp,isblck,-1
2109 WRITE(*,*)
' Cons. block ',
ncblck, isblck, icgrp, labelf, labell
2127 ncnmxb=max(ncnmxb,ncon)
2128 nprmxb=max(nprmxb,npar)
2153 ncnmxg=max(ncnmxg,ncon)
2154 nprmxg=max(nprmxg,npar)
2189 WRITE(*,*)
'PRPCON: constraints split into ',
ncgrp,
'(disjoint) groups,'
2190 WRITE(*,*)
' groups combined into ',
ncblck,
'(non overlapping) blocks'
2191 WRITE(*,*)
' max group size (cons., par.) ', ncnmxg, nprmxg
2192 WRITE(*,*)
' max block size (cons., par.) ', ncnmxb, nprmxb
2210 INTEGER(mpi) :: icgb
2211 INTEGER(mpi) :: icgrp
2213 INTEGER(mpi) :: ifirst
2214 INTEGER(mpi) :: ilast
2215 INTEGER(mpl) :: ioffc
2216 INTEGER(mpl) :: ioffp
2217 INTEGER(mpi) :: irank
2218 INTEGER(mpi) :: ipar0
2219 INTEGER(mpi) :: itgbi
2220 INTEGER(mpi) :: ivgb
2222 INTEGER(mpi) :: jcgb
2224 INTEGER(mpi) :: label
2225 INTEGER(mpi) :: ncon
2226 INTEGER(mpi) :: npar
2227 INTEGER(mpi) :: nrank
2228 INTEGER(mpi) :: inone
2233 INTEGER(mpl):: length
2234 REAL(mpd),
DIMENSION(:),
ALLOCATABLE :: matConstraintsT
2235 REAL(mpd),
DIMENSION(:),
ALLOCATABLE :: auxVectorD
2236 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: auxVectorI
2240 IF(
ncgb == 0)
RETURN
2249 CALL mpalloc(auxvectori,length,
'auxiliary array (I)')
2250 CALL mpalloc(auxvectord,length,
'auxiliary array (D)')
2253 CALL mpalloc(matconstraintst,length,
'transposed matrix of constraints (blocks)')
2254 matconstraintst=0.0_mpd
2267 WRITE(*,*)
' Constraint group, #con, rank', icgrp, ncon, 0,
' (empty)'
2270 DO jcgb=ifirst,ilast
2282 IF(ivgb > 0) matconstraintst(int(jcgb-ifirst,mpl)*int(npar,mpl)+ivgb-ipar0+ioffc)= &
2283 matconstraintst(int(jcgb-ifirst,mpl)*int(npar,mpl)+ivgb-ipar0+ioffc)+factr
2290 DO ll=ioffc+1,ioffc+npar
2296 matconstraintst(int(i-1,mpl)*int(npar,mpl)+ll)* &
2297 matconstraintst(int(j-1,mpl)*int(npar,mpl)+ll)
2303 IF (
icheck > 1 .OR. irank < ncon)
THEN
2304 WRITE(*,*)
' Constraint group, #con, rank', icgrp, ncon, irank
2305 IF (irank < ncon)
THEN
2306 WRITE(*,*)
' .. rank deficit !! '
2307 WRITE(*,*)
' E.g. fix all parameters and remove all constraints related to label ', &
2312 ioffc=ioffc+int(npar,mpl)*int(ncon,mpl)
2319 WRITE(*,*)
'Rank of product matrix of constraints is',nrank, &
2320 ' for',
ncgb,
' constraint equations'
2321 WRITE(8,*)
'Rank of product matrix of constraints is',nrank, &
2322 ' for',
ncgb,
' constraint equations'
2323 IF(nrank <
ncgb)
THEN
2324 WRITE(*,*)
'Warning: insufficient constraint equations!'
2325 WRITE(8,*)
'Warning: insufficient constraint equations!'
2328 WRITE(*,*)
' --> enforcing SUBITO mode'
2329 WRITE(8,*)
' --> enforcing SUBITO mode'
2336 print *,
'QL decomposition of constraints matrix'
2339 WRITE(
lunlog,*)
'QL decomposition of constraints matrix'
2351 CALL lpqldec(matconstraintst,evmin,evmax)
2355 print *,
' largest |eigenvalue| of L: ', evmax
2356 print *,
' smallest |eigenvalue| of L: ', evmin
2357 IF (evmin == 0.0_mpd.AND.
icheck == 0)
THEN
2358 CALL peend(27,
'Aborted, singular QL decomposition of constraints matrix')
2359 stop
'FEASMA: stopping due to singular QL decomposition of constraints matrix'
2385 INTEGER(mpi) :: icgb
2386 INTEGER(mpi) :: icgrp
2387 INTEGER(mpi) :: iter
2388 INTEGER(mpi) :: itgbi
2389 INTEGER(mpi) :: ivgb
2390 INTEGER(mpi) :: ieblck
2391 INTEGER(mpi) :: isblck
2392 INTEGER(mpi) :: ifirst
2393 INTEGER(mpi) :: ilast
2395 INTEGER(mpi) :: jcgb
2396 INTEGER(mpi) :: label
2397 INTEGER(mpi) :: inone
2398 INTEGER(mpi) :: ncon
2400 REAL(mps),
INTENT(IN) :: concut
2401 INTEGER(mpi),
INTENT(OUT) :: iact
2408 REAL(mpd),
DIMENSION(:),
ALLOCATABLE :: vecCorrections
2412 IF(
ncgb == 0)
RETURN
2442 sum1=sqrt(sum1/real(
ncgb,mpd))
2443 sum2=sum2/real(
ncgb,mpd)
2445 IF(iter == 1.AND.sum1 < concut)
RETURN
2447 IF(iter == 1.AND.
ncgb <= 12)
THEN
2449 WRITE(*,*)
'Constraint equation discrepancies:'
2451101
FORMAT(4x,4(i5,g12.4))
2453103
FORMAT(10x,
' Cut on rms value is',g8.1)
2458 WRITE(*,*)
'Improve constraints'
2462 WRITE(*,102) iter,sum1,sum2,sum3
2463102
FORMAT(i6,
' rms',g12.4,
' avrg_abs',g12.4,
' max_abs',g12.4)
2465 CALL mpalloc(veccorrections,int(
nvgb,mpl),
'constraint corrections')
2466 veccorrections=0.0_mpd
2474 ieblck=isblck+(ncon*(ncon+1))/2
2543 INTEGER(mpi) :: iact
2544 INTEGER(mpi) :: ierrc
2545 INTEGER(mpi) :: ierrf
2546 INTEGER(mpi) :: ioffp
2548 INTEGER(mpi) :: ithr
2549 INTEGER(mpi) :: jfile
2550 INTEGER(mpi) :: jrec
2552 INTEGER(mpi) :: kfile
2555 INTEGER(mpi) :: mpri
2557 INTEGER(mpi) :: nact
2558 INTEGER(mpi) :: nbuf
2559 INTEGER(mpi) :: ndata
2560 INTEGER(mpi) :: noff
2561 INTEGER(mpi) :: noffs
2562 INTEGER(mpi) :: npointer
2563 INTEGER(mpi) :: npri
2567 INTEGER(mpi) :: nrpr
2568 INTEGER(mpi) :: nthr
2569 INTEGER(mpi) :: ntot
2570 INTEGER(mpi) :: maxRecordSize
2571 INTEGER(mpi) :: maxRecordFile
2573 INTEGER(mpi),
INTENT(OUT) :: more
2583 CHARACTER (LEN=7) :: cfile
2588 SUBROUTINE readc(bufferD, bufferF, bufferI, bufferLength, lun, err)
BIND(c)
2590 REAL(c_double),
DIMENSION(*),
INTENT(OUT) :: bufferD
2591 REAL(c_float),
DIMENSION(*),
INTENT(OUT) :: bufferF
2592 INTEGER(c_int),
DIMENSION(*),
INTENT(OUT) :: bufferI
2593 INTEGER(c_int),
INTENT(INOUT) :: bufferLength
2594 INTEGER(c_int),
INTENT(IN),
VALUE :: lun
2595 INTEGER(c_int),
INTENT(OUT) :: err
2596 END SUBROUTINE readc
2602 DATA npri / 0 /, mpri / 1000 /
2651 files:
DO WHILE (jfile > 0)
2655 CALL binopn(kfile,ithr,ios)
2661 IF(kfile <=
nfilf)
THEN
2682 eof=(ierrc <= 0.AND.ierrc /= -4)
2683 IF(eof.AND.ierrc < 0)
THEN
2684 WRITE(*,*)
'Read error for binary Cfile', kfile,
'record', jrec+1,
':', ierrc
2685 WRITE(8,*)
'Read error for binary Cfile', kfile,
'record', jrec+1,
':', ierrc
2687 WRITE(cfile,
'(I7)') kfile
2688 CALL peend(18,
'Aborted, read error(s) for binary file ' // cfile)
2689 stop
'PEREAD: stopping due to read errors'
2691 IF (
kfd(1,jfile) == 1)
THEN
2697 IF(eof)
EXIT records
2702 xfd(jfile)=max(
xfd(jfile),n)
2719 IF ((noff+nr+2+
ndimbuf >= ndata*(iact+1)).OR.(nbuf >= npointer))
EXIT files
2734 IF (
kfd(1,jfile) == 1)
THEN
2735 print *,
'PEREAD: file ', kfile,
'read the first time, found',jrec,
' records'
2739 IF (-
kfd(1,jfile) /= jrec)
THEN
2740 WRITE(cfile,
'(I7)') kfile
2741 CALL peend(19,
'Aborted, binary file modified (length) ' // cfile)
2742 stop
'PEREAD: file modified (length)'
2793 IF (
kfd(1,jfile) == 1)
kfd(1,jfile)=-nrc
2812 DO WHILE (
nloopn == 0.AND.ntot >= nrpr)
2813 WRITE(*,*)
' Record ',nrpr
2814 IF (nrpr < 100000)
THEN
2823 IF (npri == 1)
WRITE(*,100)
2825100
FORMAT(/
' PeRead records active file' &
2826 /
' total block threads number')
2827101
FORMAT(
' PeRead',4i10)
2840 IF (
xfd(k) > maxrecordsize)
THEN
2841 maxrecordsize=
xfd(k)
2844 dw=real(-
kfd(1,k),mpd)
2847 ds1=ds1+dw*real(
wfd(k),mpd)
2848 ds2=ds2+dw*real(
wfd(k)**2,mpd)
2850 print *,
'PEREAD: file ', maxrecordfile,
'with max record size ', maxrecordsize
2851 IF (
nfilw > 0.AND.ds0 > 0.0_mpd)
THEN
2855 WRITE(lun,177)
nfilw,real(ds1,mps),real(ds2,mps)
2856177
FORMAT(/
' !!!!!',i4,
' weighted binary files', &
2857 /
' !!!!! mean, variance of weights =',2g12.4)
2875179
FORMAT(/
' Read cache usage (#blocks, #records, ', &
2876 'min,max records/block'/17x,i10,i12,2i10)
2894 INTEGER(mpi),
INTENT(IN) :: mode
2896 INTEGER(mpi) :: ibuf
2897 INTEGER(mpi) :: ichunk
2899 INTEGER(mpi) :: itgbi
2905 INTEGER(mpi),
PARAMETER :: maxbad = 100
2906 INTEGER(mpi) :: nbad
2907 INTEGER(mpi) :: nerr
2908 INTEGER(mpi) :: inone
2926 CALL isjajb(nst,ist,ja,jb,jsp)
2948 CALL pechk(ibuf,nerr)
2951 IF(nbad >= maxbad)
EXIT
2956 CALL isjajb(nst,ist,ja,jb,jsp)
2969 CALL peend(20,
'Aborted, bad binary records')
2970 stop
'PEREAD: stopping due to bad records'
2991 INTEGER(mpi) :: ioff
2998 INTEGER(mpi),
INTENT(IN) :: ibuf
2999 INTEGER(mpi),
INTENT(OUT) :: nerr
3008 outer:
DO WHILE(is < nst)
3013 IF(is > nst)
EXIT outer
3019 IF(is > nst)
EXIT outer
3036100
FORMAT(
' PEREAD: record ', i8,
' in file ',i6,
' is broken !!!')
3046101
FORMAT(
' PEREAD: record ', i8,
' in file ',i6,
' contains ', i6,
' NaNs !!!')
3062 INTEGER(mpi) :: ibuf
3063 INTEGER(mpi) :: ichunk
3064 INTEGER(mpi) :: iproc
3065 INTEGER(mpi) :: ioff
3066 INTEGER(mpi) :: ioffbi
3068 INTEGER(mpi) :: itgbi
3073 INTEGER(mpi) :: nalg
3075 INTEGER(mpi) :: inone
3076 INTEGER(mpl) :: length
3110 CALL isjajb(nst,ist,ja,jb,jsp)
3134 CALL isjajb(nst,ist,ja,jb,jsp)
3155 CALL isjajb(nst,ist,ja,jb,jsp)
3185 INTEGER(mpi) :: istart
3186 INTEGER(mpi) :: itgbi
3188 INTEGER(mpi) :: jstart
3189 INTEGER(mpi) :: jtgbi
3190 INTEGER(mpi) :: lstart
3191 INTEGER(mpi) :: ltgbi
3193 INTEGER(mpi),
INTENT(IN) :: inds
3194 INTEGER(mpi),
INTENT(IN) :: inde
3196 IF (inds > inde)
RETURN
3206 IF (istart == 0)
THEN
3207 IF (itgbi /= ltgbi+1)
THEN
3210 IF (lstart == 0)
THEN
3225 IF (istart /= jstart)
THEN
3236 IF (itgbi /= ltgbi+1)
THEN
3251 IF (istart /= jstart)
THEN
3303 INTEGER(mpi),
INTENT(IN) :: nst
3304 INTEGER(mpi),
INTENT(IN OUT) :: is
3305 INTEGER(mpi),
INTENT(OUT) :: ja
3306 INTEGER(mpi),
INTENT(OUT) :: jb
3307 INTEGER(mpi),
INTENT(OUT) :: jsp
3315 IF(is >= nst)
RETURN
3369 INTEGER(mpi) :: ioffb
3371 INTEGER(mpi) :: itgbi
3372 INTEGER(mpi) :: itgbij
3373 INTEGER(mpi) :: itgbik
3374 INTEGER(mpi) :: ivgb
3375 INTEGER(mpi) :: ivgbij
3376 INTEGER(mpi) :: ivgbik
3379 INTEGER(mpi) :: lastit
3381 INTEGER(mpi) :: ncrit
3382 INTEGER(mpi) :: ngras
3383 INTEGER(mpi) :: nparl
3385 INTEGER(mpi) :: nrej
3386 INTEGER(mpi) :: inone
3387 INTEGER(mpi) :: ilow
3388 INTEGER(mpi) :: nlow
3389 INTEGER(mpi) :: nzero
3409 CALL gmpdef(1,4,
'Function value in iterations')
3411 CALL gmpdef(2,3,
'Number of MINRES steps vs iteration nr')
3413 CALL hmpdef( 5,0.0,0.0,
'Number of degrees of freedom')
3414 CALL hmpdef(11,0.0,0.0,
'Number of local parameters')
3415 CALL hmpdef(23,0.0,0.0,
'SQRT of diagonal elements without presigma')
3416 CALL hmpdef(24,0.0,0.0,
'Log10 of off-diagonal elements')
3417 CALL hmpdef(25,0.0,0.0,
'Relative individual pre-sigma')
3418 CALL hmpdef(26,0.0,0.0,
'Relative global pre-sigma')
3423 'Normalized residuals of single (global) measurement')
3425 'Normalized residuals of single (local) measurement')
3427 'Pulls of single (global) measurement')
3429 'Pulls of single (local) measurement')
3430 CALL hmpdef(4,0.0,0.0,
'Chi^2/Ndf after local fit')
3431 CALL gmpdef(4,5,
'location, dispersion (res.) vs record nr')
3432 CALL gmpdef(5,5,
'location, dispersion (pull) vs record nr')
3446 IF(
nloopn == 2)
CALL hmpdef(6,0.0,0.0,
'Down-weight fraction')
3449 IF(
iterat /= lastit)
THEN
3457 ELSE IF(
iterat >= 1)
THEN
3507 WRITE(*,111) nparl,ncrit,usei,used,peaki,peakd
3508111
FORMAT(
' Write cache usage (#flush,#overrun,<levels>,', &
3509 'peak(levels))'/2i7,
',',4(f6.1,
'%'))
3517 ELSE IF (
mbandw > 0)
THEN
3549 print *,
" ... warning ..."
3550 print *,
" global parameters with too few (< MREQENA) accepted entries: ", nlow
3554 IF(
icalcm == 1 .AND. nzero > 0)
THEN
3556 WRITE(*,*)
'Warning: the rank defect of the symmetric',
nfgb, &
3557 '-by-',
nfgb,
' matrix is ',
ndefec,
' (should be zero).'
3558 WRITE(lun,*)
'Warning: the rank defect of the symmetric',
nfgb, &
3559 '-by-',
nfgb,
' matrix is ',
ndefec,
' (should be zero).'
3562 WRITE(*,*)
' --> enforcing SUBITO mode'
3563 WRITE(lun,*)
' --> enforcing SUBITO mode'
3573 elmt=real(matij(i,i),mps)
3574 IF(elmt > 0.0)
CALL hmpent(23,1.0/sqrt(elmt))
3608 CALL addsums(1, adder, 0, 1.0_mpl)
3622 IF(sgm > 0.0) weight=1.0/sgm**2
3638 IF(itgbij /= 0)
THEN
3653 IF(
icalcm == 1.AND.ivgbij > 0)
THEN
3660 IF(ivgbij > 0.AND.ivgbik > 0)
THEN
3661 CALL mupdat(ivgbij,ivgbik,weight*factrj*factrk)
3667 adder=weight*dsum**2
3668 CALL addsums(1, adder, 1, 1.0_mpl)
3686 ratae=max(-99.9,ratae)
3695 WRITE(*,*)
'Data rejected in initial loop:'
3697 nrejec(0),
' (rank deficit/NaN) ',
nrejec(1),
' (Ndf=0) ', &
3746 WRITE(lun,*)
' === local fits have bordered band matrix structure ==='
3747 IF (
nbndr(1) > 0 )
WRITE(lun,101)
' NBNDR',
nbndr(1),
'number of records (upper/left border)'
3748 IF (
nbndr(2) > 0 )
WRITE(lun,101)
' NBNDR',
nbndr(2),
'number of records (lower/right border)'
3749 WRITE(lun,101)
' NBDRX',
nbdrx,
'max border size'
3750 WRITE(lun,101)
' NBNDX',
nbndx,
'max band width'
3759101
FORMAT(1x,a8,
' =',i14,
' = ',a)
3774 INTEGER(mpi),
INTENT(IN) :: lunp
3779101
FORMAT(
' it fc',
' fcn_value dfcn_exp slpr costh iit st', &
3780 ' ls step cutf',1x,
'rejects hhmmss FMS')
3781102
FORMAT(
' -- --',
' ----------- -------- ---- ----- --- --', &
3782 ' -- ----- ----',1x,
'------- ------ ---')
3798 INTEGER(mpi) :: nrej
3799 INTEGER(mpi) :: nsecnd
3803 REAL(mps) :: slopes(3)
3804 REAL(mps) :: steps(3)
3805 REAL,
DIMENSION(2) :: ta
3808 INTEGER(mpi),
INTENT(IN) :: lunp
3810 CHARACTER (LEN=4):: ccalcm(4)
3811 DATA ccalcm /
' end',
' S',
' F ',
' FMS' /
3815 IF(nrej > 9999999) nrej=9999999
3819 nsecnd=nint(secnd,mpi)
3828 CALL ptlopt(nfa,ma,slopes,steps)
3829 ratae=max(-99.9,min(99.9,slopes(2)/slopes(1)))
3835103
FORMAT(i3,i3,e12.5,38x,f5.1, 1x,i7, i3,i2.2,i2.2,a4)
3836104
FORMAT(i3,i3,e12.5,1x,e8.2,f6.3,f6.3,i5,2i3,f6.3,f5.1, &
3837 1x,i7, i3,i2.2,i2.2,a4)
3838105
FORMAT(i3,i3,e12.5,1x,e8.2,12x,i5,i3,9x,f5.1, &
3839 1x,i7, i3,i2.2,i2.2,a4)
3852 INTEGER(mpi) :: minut
3854 INTEGER(mpi) :: nhour
3855 INTEGER(mpi) :: nrej
3856 INTEGER(mpi) :: nsecnd
3860 REAL(mps) :: slopes(3)
3861 REAL(mps) :: steps(3)
3862 REAL,
DIMENSION(2) :: ta
3865 INTEGER(mpi),
INTENT(IN) :: lunp
3866 CHARACTER (LEN=4):: ccalcm(4)
3867 DATA ccalcm /
' end',
' S',
' F ',
' FMS' /
3871 IF(nrej > 9999999) nrej=9999999
3875 nsecnd=nint(secnd,mpi)
3879 CALL ptlopt(nfa,ma,slopes,steps)
3880 ratae=abs(slopes(2)/slopes(1))
3885104
FORMAT(3x,i3,e12.5,9x, 35x, i7, i3,i2.2,i2.2,a4)
3886105
FORMAT(3x,i3,e12.5,9x, f6.3,14x,i3,f6.3,6x, i7, i3,i2.2,i2.2,a4)
3900 INTEGER(mpi) :: nsecnd
3903 REAL,
DIMENSION(2) :: ta
3906 INTEGER(mpi),
INTENT(IN) :: lunp
3907 CHARACTER (LEN=4):: ccalcm(4)
3908 DATA ccalcm /
' end',
' S',
' F ',
' FMS' /
3913 nsecnd=nint(secnd,mpi)
3915 WRITE(lunp,106) nhour,minut,nsecnd,ccalcm(
lcalcm)
3916106
FORMAT(69x,i3,i2.2,i2.2,a4)
3926 INTEGER(mpi) :: lunit
3928 WRITE(lunit,102)
'Explanation of iteration table'
3929 WRITE(lunit,102)
'=============================='
3930 WRITE(lunit,101)
'it', &
3931 'iteration number. Global parameters are improved for it > 0.'
3932 WRITE(lunit,102)
'First function evaluation is called iteraton 0.'
3933 WRITE(lunit,101)
'fc',
'number of function evaluations.'
3934 WRITE(lunit,101)
'fcn_value',
'value of 2 x Likelihood function (LF).'
3935 WRITE(lunit,102)
'The final value is the chi^2 value of the fit and should'
3936 WRITE(lunit,102)
'be about equal to the NDF (see below).'
3937 WRITE(lunit,101)
'dfcn_exp', &
3938 'expected reduction of the value of the Likelihood function (LF)'
3939 WRITE(lunit,101)
'slpr',
'ratio of the actual slope to inital slope.'
3940 WRITE(lunit,101)
'costh', &
3941 'cosine of angle between search direction and -gradient'
3943 WRITE(lunit,101)
'iit', &
3944 'number of internal iterations in MINRES algorithm'
3945 WRITE(lunit,101)
'st',
'stop code of MINRES algorithm'
3946 WRITE(lunit,102)
'< 0: rhs is very special, with beta2 = 0'
3947 WRITE(lunit,102)
'= 0: rhs b = 0, i.e. the exact solution is x = 0'
3948 WRITE(lunit,102)
'= 1 requested accuracy achieved, as determined by rtol'
3949 WRITE(lunit,102)
'= 2 reasonable accuracy achieved, given eps'
3950 WRITE(lunit,102)
'= 3 x has converged to an eigenvector'
3951 WRITE(lunit,102)
'= 4 matrix ill-conditioned (Acond has exceeded 0.1/eps)'
3952 WRITE(lunit,102)
'= 5 the iteration limit was reached'
3953 WRITE(lunit,102)
'= 6 Matrix x vector does not define a symmetric matrix'
3954 WRITE(lunit,102)
'= 7 Preconditioner does not define a symmetric matrix'
3955 ELSEIF (
metsol == 5)
THEN
3956 WRITE(lunit,101)
'iit', &
3957 'number of internal iterations in MINRES-QLP algorithm'
3958 WRITE(lunit,101)
'st',
'stop code of MINRES-QLP algorithm'
3959 WRITE(lunit,102)
'= 1: beta_{k+1} < eps, iteration k is the final Lanczos step.'
3960 WRITE(lunit,102)
'= 2: beta2 = 0. If M = I, b and x are eigenvectors of A.'
3961 WRITE(lunit,102)
'= 3: beta1 = 0. The exact solution is x = 0.'
3962 WRITE(lunit,102)
'= 4: A solution to (poss. singular) Ax = b found, given rtol.'
3963 WRITE(lunit,102)
'= 5: A solution to (poss. singular) Ax = b found, given eps.'
3964 WRITE(lunit,102)
'= 6: Pseudoinverse solution for singular LS problem, given rtol.'
3965 WRITE(lunit,102)
'= 7: Pseudoinverse solution for singular LS problem, given eps.'
3966 WRITE(lunit,102)
'= 8: The iteration limit was reached.'
3967 WRITE(lunit,102)
'= 9: The operator defined by Aprod appears to be unsymmetric.'
3968 WRITE(lunit,102)
'=10: The operator defined by Msolve appears to be unsymmetric.'
3969 WRITE(lunit,102)
'=11: The operator defined by Msolve appears to be indefinite.'
3970 WRITE(lunit,102)
'=12: xnorm has exceeded maxxnorm or will exceed it next iteration.'
3971 WRITE(lunit,102)
'=13: Acond has exceeded Acondlim or 0.1/eps.'
3972 WRITE(lunit,102)
'=14: Least-squares problem but no converged solution yet.'
3973 WRITE(lunit,102)
'=15: A null vector obtained, given rtol.'
3975 WRITE(lunit,101)
'ls',
'line search info'
3976 WRITE(lunit,102)
'< 0 recalculate function'
3977 WRITE(lunit,102)
'= 0: N or STP lt 0 or step not descending'
3978 WRITE(lunit,102)
'= 1: Linesearch convergence conditions reached'
3979 WRITE(lunit,102)
'= 2: interval of uncertainty at lower limit'
3980 WRITE(lunit,102)
'= 3: max nr of line search calls reached'
3981 WRITE(lunit,102)
'= 4: step at the lower bound'
3982 WRITE(lunit,102)
'= 5: step at the upper bound'
3983 WRITE(lunit,102)
'= 6: rounding error limitation'
3984 WRITE(lunit,101)
'step', &
3985 'the factor for the Newton step during the line search. Usually'
3987 'a value of 1 gives a sufficient reduction of the LF. Oherwise'
3988 WRITE(lunit,102)
'other step values are tried.'
3989 WRITE(lunit,101)
'cutf', &
3990 'cut factor. Local fits are rejected, if their chi^2 value'
3992 'is larger than the 3-sigma chi^2 value times the cut factor.'
3993 WRITE(lunit,102)
'A cut factor of 1 is used finally, but initially a larger'
3994 WRITE(lunit,102)
'factor may be used. A value of 0.0 means no cut.'
3995 WRITE(lunit,101)
'rejects',
'total number of rejected local fits.'
3996 WRITE(lunit,101)
'hmmsec',
'the time in hours (h), minutes (mm) and seconds.'
3997 WRITE(lunit,101)
'FMS',
'calculation of Function value, Matrix, Solution.'
4000101
FORMAT(a9,
' = ',a)
4017 INTEGER(mpi),
INTENT(IN) :: i
4018 INTEGER(mpi),
INTENT(IN) :: j
4019 REAL(mpd),
INTENT(IN) :: add
4021 INTEGER(mpl):: ijadd
4022 INTEGER(mpl):: ijcsr3
4027 IF(i <= 0.OR.j <= 0.OR. add == 0.0_mpd)
RETURN
4048 ELSE IF(
matsto == 2)
THEN
4076 CALL peend(23,
'Aborted, bad matrix index')
4077 stop
'mupdat: bad index'
4102 INTEGER(mpi),
INTENT(IN) :: i
4103 INTEGER(mpi),
INTENT(IN) :: j1
4104 INTEGER(mpi),
INTENT(IN) :: j2
4105 INTEGER(mpi),
INTENT(IN) :: il
4106 INTEGER(mpi),
INTENT(IN) :: jl
4107 INTEGER(mpi),
INTENT(IN) :: n
4108 REAL(mpd),
INTENT(IN) :: sub((n*n+n)/2)
4115 INTEGER(mpi):: iblast
4116 INTEGER(mpi):: iblock
4122 INTEGER(mpi):: jblast
4123 INTEGER(mpi):: jblock
4133 IF(i <= 0.OR.j1 <= 0.OR.j2 > i)
RETURN
4147 print *,
' MGUPDT: ij=0', i,j1,j2,il,jl,ij,lr,iprc,
matsto
4154 ijl=(k*k-k)/2+jl+ir-ia1
4171 ijl=(k*k-k)/2+jl+ir-ia1
4175 IF (jblock /= jblast.OR.iblock /= iblast)
THEN
4199 CALL ijpgrp(i,j1,ij,lr,iprc)
4201 print *,
' MGUPDT: ij=0', i,j1,j2,il,jl,ij,lr,iprc,
matsto
4265SUBROUTINE loopbf(nrej,numfil,naccf,chi2f,ndff)
4292 INTEGER(mpi) :: ibuf
4293 INTEGER(mpi) :: ichunk
4294 INTEGER(mpl) :: icmn
4295 INTEGER(mpl) :: icost
4297 INTEGER(mpi) :: idiag
4299 INTEGER(mpi) :: iext
4307 INTEGER(mpi) :: imeas
4310 INTEGER(mpi) :: ioffb
4311 INTEGER(mpi) :: ioffc
4312 INTEGER(mpi) :: ioffd
4313 INTEGER(mpi) :: ioffe
4314 INTEGER(mpi) :: ioffi
4315 INTEGER(mpi) :: ioffq
4316 INTEGER(mpi) :: iprc
4317 INTEGER(mpi) :: iprcnx
4318 INTEGER(mpi) :: iprdbg
4319 INTEGER(mpi) :: iproc
4320 INTEGER(mpi) :: irbin
4321 INTEGER(mpi) :: isize
4323 INTEGER(mpi) :: iter
4324 INTEGER(mpi) :: itgbi
4325 INTEGER(mpi) :: ivgbj
4326 INTEGER(mpi) :: ivgbk
4327 INTEGER(mpi) :: ivpgrp
4337 INTEGER(mpi) :: joffd
4338 INTEGER(mpi) :: joffi
4339 INTEGER(mpi) :: jproc
4343 INTEGER(mpi) :: kbdr
4344 INTEGER(mpi) :: kbdrx
4345 INTEGER(mpi) :: kbnd
4348 INTEGER(mpi) :: lvpgrp
4349 INTEGER(mpi) :: mbdr
4350 INTEGER(mpi) :: mbnd
4351 INTEGER(mpi) :: mside
4352 INTEGER(mpi) :: nalc
4353 INTEGER(mpi) :: nalg
4357 INTEGER(mpi) :: ndown
4359 INTEGER(mpi) :: nfred
4360 INTEGER(mpi) :: nfrei
4362 INTEGER(mpi) :: nprdbg
4363 INTEGER(mpi) :: nrank
4366 INTEGER(mpi) :: nter
4367 INTEGER(mpi) :: nweig
4368 INTEGER(mpi) :: ngrp
4369 INTEGER(mpi) :: npar
4371 INTEGER(mpi),
INTENT(IN OUT) :: nrej(0:3)
4372 INTEGER(mpi),
INTENT(IN) :: numfil
4373 INTEGER(mpi),
INTENT(IN OUT) :: naccf(numfil)
4374 REAL(mps),
INTENT(IN OUT) :: chi2f(numfil)
4375 INTEGER(mpi),
INTENT(IN OUT) :: ndff(numfil)
4382 INTEGER(mpi) :: ijprec
4389 CHARACTER (LEN=3):: chast
4390 DATA chuber/1.345_mpd/
4391 DATA cauchy/2.3849_mpd/
4439 IF(
nloopn == 1.AND.mod(nrc,100000_mpl) == 0)
THEN
4440 WRITE(*,*)
'Record',nrc,
' ... still reading'
4441 IF(
monpg1>0)
WRITE(
lunlog,*)
'Record',nrc,
' ... still reading'
4451 IF(nrc ==
nrecpr) lprnt=.true.
4452 IF(nrc ==
nrecp2) lprnt=.true.
4453 IF(nrc ==
nrecer) lprnt=.true.
4458 IF (nprdbg == 1) iprdbg=iproc
4459 IF (iproc /= iprdbg) lprnt=.false.
4464 WRITE(1,*)
'------------------ Loop',
nloopn, &
4465 ': Printout for record',nrc,iproc
4476 CALL isjajb(nst,ist,ja,jb,jsp)
4478 IF(imeas == 0)
WRITE(1,1121)
44831121
FORMAT(/
'Measured value and local derivatives'/ &
4484 ' i measured std_dev index...derivative ...')
44851122
FORMAT(i3,2g12.4,3(i3,g12.4)/(27x,3(i3,g12.4)))
4491 CALL isjajb(nst,ist,ja,jb,jsp)
4493 IF(imeas == 0)
WRITE(1,1123)
45051123
FORMAT(/
'Global derivatives'/ &
4506 ' i label gindex vindex derivative ...')
45071124
FORMAT(i3,2(i9,i7,i7,g12.4)/(3x,2(i9,i7,i7,g12.4)))
45081125
FORMAT(i3,2(i9,i7,i7,g12.4))
4518 WRITE(1,*)
'Data corrections using values of global parameters'
4519 WRITE(1,*)
'=================================================='
4529 CALL isjajb(nst,ist,ja,jb,jsp)
4561101
FORMAT(
' index measvalue corrvalue sigma')
4562102
FORMAT(i6,2x,2g12.4,
' +-',g12.4)
4564 IF(nalc <= 0)
GO TO 90
4566 ngg=(nalg*nalg+nalg)/2
4579 IF (ivpgrp /= lvpgrp)
THEN
4587 IF (npar /= nalg)
THEN
4588 print *,
' mismatch of number of global parameters ', nrc, nalg, npar, ngrp
4598 CALL peend(35,
'Aborted, mismatch of number of global parameters')
4599 stop
' mismatch of number of global parameters '
4626 DO WHILE(iter < nter)
4631 WRITE(1,*)
'Outlier-suppression iteration',iter,
' of',nter
4632 WRITE(1,*)
'=========================================='
4642 DO i=1,(nalc*nalc+nalc)/2
4652 wght =1.0_mpd/rerr**2
4657 IF(abs(resid) > chuber*rerr)
THEN
4658 wght=wght*chuber*rerr/abs(resid)
4662 wght=wght/(1.0+(resid/rerr/cauchy)**2)
4666 IF(lprnt.AND.iter /= 1.AND.nter /= 1)
THEN
4668 IF(abs(resid) > chuber*rerr) chast=
'* '
4669 IF(abs(resid) > 3.0*rerr) chast=
'** '
4670 IF(abs(resid) > 6.0*rerr) chast=
'***'
4671 IF(imeas == 0)
WRITE(1,*)
'Second loop: accumulate'
4672 IF(imeas == 0)
WRITE(1,103)
4677 WRITE(1,104) imeas,rmeas,resid,rerr,r1,chast,r2
4679103
FORMAT(
' index corrvalue residuum sigma', &
4681104
FORMAT(i6,2x,2g12.4,
' +-',g12.4,f7.2,1x,a3,f8.2)
4696 im=min(nalc+1-ij,nalc+1-ik)
4703 IF (iter == 1.AND.nalc > 5.AND.
lfitbb > 0)
THEN
4706 icmn=int(nalc,mpl)**3
4712 icost=6*int(nalc-kbdr,mpl)*int(kbnd+kbdr+1,mpl)**2+2*int(kbdr,mpl)**3
4713 IF (icost < icmn)
THEN
4725 kbdr=max(kbdr,
ibandh(k+nalc))
4726 icost=6*int(nalc-kbdr,mpl)*int(kbnd+kbdr+1,mpl)**2+2*int(kbdr,mpl)**3
4727 IF (icost < icmn)
THEN
4751 IF (
icalcm == 1.OR.lprnt) inv=2
4752 IF (mside == 1)
THEN
4767 IF ((.NOT.(
blvec(k) <= 0.0_mpd)).AND. (.NOT.(
blvec(k) > 0.0_mpd))) nan=nan+1
4772 WRITE(1,*)
'Parameter determination:',nalc,
' parameters,',
' rank=',nrank
4773 WRITE(1,*)
'-----------------------'
4774 IF(ndown /= 0)
WRITE(1,*)
' ',ndown,
' data down-weighted'
4790 wght =1.0_mpd/rerr**2
4814 IF (0.999999_mpd/wght > dvar)
THEN
4815 pull=rmeas/sqrt(1.0_mpd/wght-dvar)
4818 CALL hmpent(13,real(pull,mps))
4819 CALL gmpms(5,rec,real(pull,mps))
4821 CALL hmpent(14,real(pull,mps))
4840 IF(iter == 1.AND.jb < ist.AND.lhist) &
4841 CALL gmpms(4,rec,real(rmeas/rerr,mps))
4843 dchi2=wght*rmeas*rmeas
4848 IF(abs(resid) > chuber*rerr)
THEN
4849 wght=wght*chuber*rerr/abs(resid)
4850 dchi2=2.0*chuber*(abs(resid)/rerr-0.5*chuber)
4853 wght=wght/(1.0_mpd+(resid/rerr/cauchy)**2)
4854 dchi2=log(1.0_mpd+(resid/rerr/cauchy)**2)*cauchy**2
4864 IF(abs(resid) > chuber*rerr) chast=
'* '
4865 IF(abs(resid) > 3.0*rerr) chast=
'** '
4866 IF(abs(resid) > 6.0*rerr) chast=
'***'
4867 IF(imeas == 0)
WRITE(1,*)
'Third loop: single residuals'
4868 IF(imeas == 0)
WRITE(1,105)
4872 IF(resid < 0.0) r1=-r1
4873 IF(resid < 0.0) r2=-r2
4876105
FORMAT(
' index corrvalue residuum sigma', &
4878106
FORMAT(i6,2x,2g12.4,
' +-',g12.4,f7.2,1x,a3,f8.2)
4880 IF(iter == nter)
THEN
4882 resmax=max(resmax,abs(rmeas)/rerr)
4885 IF(iter == 1.AND.lhist)
THEN
4887 CALL hmpent( 3,real(rmeas/rerr,mps))
4889 CALL hmpent(12,real(rmeas/rerr,mps))
4896 resing=(real(nweig,mps)-real(suwt,mps))/real(nweig,mps)
4898 IF(iter == 1)
CALL hmpent( 5,real(ndf,mps))
4899 IF(iter == 1)
CALL hmpent(11,real(nalc,mps))
4904 WRITE(1,*)
'Chi^2=',summ,
' at',ndf,
' degrees of freedom: ', &
4905 '3-sigma limit is',chindl(3,ndf)*real(ndf,mps)
4906 WRITE(1,*) suwt,
' is sum of factors, compared to',nweig, &
4907 ' Downweight fraction:',resing
4909 IF(nrank /= nalc.OR.nan > 0)
THEN
4913 WRITE(1,*)
' rank deficit/NaN ', nalc, nrank, nan
4914 WRITE(1,*)
' ---> rejected!'
4921 WRITE(1,*)
' ---> rejected!'
4926 chndf=real(summ/real(ndf,mpd),mps)
4928 IF(iter == 1.AND.lhist)
CALL hmpent(4,chndf)
4941 chichi=chindl(3,ndf)*real(ndf,mps)
4945 IF(summ >
chhuge*chichi)
THEN
4948 WRITE(1,*)
' ---> rejected!'
4956 IF(summ > chlimt)
THEN
4958 WRITE(1,*)
' ---> rejected!'
4962 CALL addsums(iproc+1, dchi2, ndf, dw1)
4972 CALL addsums(iproc+1, dchi2, ndf, dw1)
4976 WRITE(1,*)
' ---> rejected!'
4989 naccf(kfl)=naccf(kfl)+1
4990 ndff(kfl) =ndff(kfl) +ndf
4991 chi2f(kfl)=chi2f(kfl)+chndf
5004 wght =1.0_mpd/rerr**2
5005 dchi2=wght*rmeas*rmeas
5009 IF(resid > chuber*rerr)
THEN
5010 wght=wght*chuber*rerr/resid
5011 dchi2=2.0*chuber*(resid/rerr-0.5*chuber)
5060 CALL addsums(iproc+1, summ, ndf, dw1)
5116 IF (nfred < 0.OR.nfrei < 0)
THEN
5143 iprcnx=ijprec(i,jnx)
5145 IF (.NOT.((jnx == j+1).AND.(iprc == iprcnx)))
THEN
5159 joffd=joffd+(il*il-il)/2
5171 WRITE(1,*)
'------------------ End of printout for record',nrc
5228 iprcnx=ijprec(i,jnx)
5230 IF (.NOT.((jnx == j+1).AND.(iprc == iprcnx)))
THEN
5245 joffd=joffd+(il*il-il)/2
5298 INTEGER(mpi) :: icom
5299 INTEGER(mpl) :: icount
5303 INTEGER(mpi) :: imin
5304 INTEGER(mpi) :: iprlim
5305 INTEGER(mpi) :: isub
5306 INTEGER(mpi) :: itgbi
5307 INTEGER(mpi) :: itgbl
5308 INTEGER(mpi) :: ivgbi
5310 INTEGER(mpi) :: label
5318 INTEGER(mpi) :: labele(3)
5319 REAL(mps):: compnt(3)
5324 CALL mvopen(lup,
'millepede.res')
5327 WRITE(*,*)
' Result of fit for global parameters'
5328 WRITE(*,*)
' ==================================='
5333 WRITE(lup,*)
'Parameter ! first 3 elements per line are', &
5334 ' significant (if used as input)'
5352 err=sqrt(abs(real(gmati,mps)))
5353 IF(gmati < 0.0_mpd) err=-err
5356 IF(gmati*diag > 0.0_mpd)
THEN
5357 gcor2=1.0_mpd-1.0_mpd/(gmati*diag)
5358 IF(gcor2 >= 0.0_mpd.AND.gcor2 <= 1.0_mpd) gcor=real(sqrt(gcor2),mps)
5363 IF(lowstat) icount=-(icount+1)
5364 IF(itgbi <= iprlim)
THEN
5378 ELSE IF(itgbi == iprlim+1)
THEN
5379 WRITE(* ,*)
'... (further printout suppressed, but see log file)'
5393 ELSE IF (
igcorr /= 0)
THEN
5411 CALL mvopen(lup,
'millepede.eve')
5421 DO isub=0,min(15,imin-1)
5443 WRITE(lup,103) (labele(ie),compnt(ie),ie=1,iev)
5447 IF(iev /= 0)
WRITE(lup,103) (labele(ie),compnt(ie),ie=1,iev)
5454101
FORMAT(1x,
' label parameter presigma differ', &
5455 ' error'/ 1x,
'-----------',4x,4(
'-------------'))
5456102
FORMAT(i10,2x,4g14.5,f8.3)
5457103
FORMAT(3(i11,f11.7,2x))
5458110
FORMAT(i10,2x,2g14.5,28x,i12)
5459111
FORMAT(i10,2x,3g14.5,14x,i12)
5460112
FORMAT(i10,2x,4g14.5,i12)
5482 INTEGER(mpi) :: icom
5483 INTEGER(mpl) :: icount
5484 INTEGER(mpi) :: ifrst
5485 INTEGER(mpi) :: ilast
5486 INTEGER(mpi) :: inext
5487 INTEGER(mpi) :: itgbi
5488 INTEGER(mpi) :: itgbl
5489 INTEGER(mpi) :: itpgrp
5490 INTEGER(mpi) :: ivgbi
5492 INTEGER(mpi) :: icgrp
5493 INTEGER(mpi) :: ipgrp
5495 INTEGER(mpi) :: jpgrp
5497 INTEGER(mpi) :: label1
5498 INTEGER(mpi) :: label2
5499 INTEGER(mpi) :: ncon
5500 INTEGER(mpi) :: npair
5501 INTEGER(mpi) :: nstep
5504 INTEGER(mpl):: length
5506 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: vecPairedParGroups
5509 SUBROUTINE ggbmap(ipgrp,npair,npgrp)
5511 INTEGER(mpi),
INTENT(IN) :: ipgrp
5512 INTEGER(mpi),
INTENT(OUT) :: npair
5513 INTEGER(mpi),
DIMENSION(:),
INTENT(OUT) :: npgrp
5521 CALL mvopen(lup,
'millepede.res')
5522 WRITE(lup,*)
'*** Results of checking input only, no solution performed ***'
5523 WRITE(lup,*)
'! === global parameters ==='
5524 WRITE(lup,*)
'! fixed-1: by pre-sigma, -2: by entries cut, -3: by iterated entries cut'
5525 WRITE(lup,*)
'! Label Value Pre-sigma Entries Cons. group Status '
5539 IF (ivgbi <= 0)
THEN
5541 IF (ivgbi == -4)
THEN
5542 WRITE(lup,116) c1,itgbl,par,presig,icount,icgrp
5544 WRITE(lup,110) c1,itgbl,par,presig,icount,icgrp,ivgbi
5548 WRITE(lup,111) c1,itgbl,par,presig,icount,icgrp
5554 WRITE(lup,*)
'!.Appearance statistics '
5555 WRITE(lup,*)
'!. Label First file and record Last file and record #files #paired-par'
5558 IF (itpgrp > 0)
THEN
5566 WRITE(lup,*)
'* === constraint groups ==='
5568 WRITE(lup,*)
'* Group #Cons. Entries First label Last label'
5570 WRITE(lup,*)
'* Group #Cons. Entries First label Last label Paired label range'
5572 CALL mpalloc(vecpairedpargroups,length,
'paired global parameter groups (I)')
5584 IF (
icheck > 1 .AND. label1 > 0)
THEN
5588 vecpairedpargroups(npair+1)=0
5592 jpgrp=vecpairedpargroups(j)
5596 IF (ifrst /= 0.AND.inext /= (ilast+nstep))
THEN
5599 WRITE(lup,114) label1, label2
5604 IF (ifrst == 0) ifrst=inext
5611 IF (jpgrp == vecpairedpargroups(j+1)-1)
THEN
5616 IF (ifrst /= 0)
THEN
5619 WRITE(lup,114) label1, label2
5625 WRITE(lup,*)
'*.Appearance statistics '
5626 WRITE(lup,*)
'*. Group First file and record Last file and record #files'
5636110
FORMAT(
' !',a1,i10,2x,2g14.5,2i12,
' fixed',i2)
5637111
FORMAT(
' !',a1,i10,2x,2g14.5,2i12,
' variable')
5638112
FORMAT(
' !.',i10,6i11)
5639113
FORMAT(
' * ',i6,i8,3i12)
5640114
FORMAT(
' *:',48x,i12,
' ..',i12)
5641115
FORMAT(
' *.',i10,5i11)
5642116
FORMAT(
' !',a1,i10,2x,2g14.5,2i12,
' redundant')
5672 INTEGER(mpi) :: iproc
5682 INTEGER(mpi),
INTENT(IN) :: n
5683 INTEGER(mpl),
INTENT(IN) :: l
5684 REAL(mpd),
INTENT(IN) :: x(n)
5685 INTEGER(mpi),
INTENT(IN) :: is
5686 INTEGER(mpi),
INTENT(IN) :: ie
5687 REAL(mpd),
INTENT(OUT) :: b(n)
5692 INTEGER(mpl) :: indij
5693 INTEGER(mpl) :: indid
5695 INTEGER(mpi) :: ichunk
5731 CALL peend(24,
'Aborted, vector/matrix size mismatch')
5732 stop
'AVPRDS: mismatched vector and matrix'
5752 IF (ia2 <= ib2) b(ia2:ib2)=b(ia2:ib2)+
globalmatd(ia2:ib2)*x(ia2:ib2)
5763 IF (i <= ie.AND.i >= is)
THEN
5771 IF (j <= ie.AND.j >= is)
THEN
5787 IF (ja2 <= jb2)
THEN
5790 b(i)=b(i)+dot_product(
globalmatd(indij+lj+ja2-ja:indij+lj+jb2-ja),x(ja2:jb2))
5794 IF (
mextnd == 0.AND.ia2 <= ib2)
THEN
5797 b(j)=b(j)+dot_product(
globalmatd(indij+lj+jn*(ia2-ia):indij+lj+jn*(ib2-ia):jn),x(ia2:ib2))
5816 IF (i <= ie.AND.i >= is)
THEN
5824 IF (j <= ie.AND.j >= is)
THEN
5840 IF (ja2 <= jb2)
THEN
5843 b(i)=b(i)+dot_product(real(
globalmatf(indij+lj+ja2-ja:indij+lj+jb2-ja),mpd),x(ja2:jb2))
5847 IF (
mextnd == 0.AND.ia2 <= ib2)
THEN
5850 b(j)=b(j)+dot_product(real(
globalmatf(indij+lj+jn*(ia2-ia):indij+lj+jn*(ib2-ia):jn),mpd),x(ia2:ib2))
5884 INTEGER(mpi) :: iproc
5892 INTEGER(mpi),
INTENT(IN) :: n
5893 INTEGER(mpl),
INTENT(IN) :: l
5894 REAL(mpd),
INTENT(IN) :: x(n)
5895 REAL(mpd),
INTENT(OUT) :: b(n)
5900 INTEGER(mpl) :: indij
5901 INTEGER(mpl) :: indid
5903 INTEGER(mpi) :: ichunk
5932 CALL peend(24,
'Aborted, vector/matrix size mismatch')
5933 stop
'AVPRD0: mismatched vector and matrix'
5979 b(i)=b(i)+dot_product(
globalmatd(indij+lj:indij+lj+jn-1),x(ja:jb))
5985 b(j)=b(j)+dot_product(
globalmatd(indij+lj:indij+jn*in:jn),x(ia:ib))
6021 b(i)=b(i)+dot_product(real(
globalmatf(indij+lj:indij+lj+jn-1),mpd),x(ja:jb))
6027 b(j)=b(j)+dot_product(real(
globalmatf(indij+lj:indij+jn*in:jn),mpd),x(ia:ib))
6051 INTEGER(mpi) :: ispc
6062 INTEGER(mpl) :: indid
6063 INTEGER(mpl),
DIMENSION(12) :: icount
6070 icount(4)=huge(icount(4))
6071 icount(7)=huge(icount(7))
6072 icount(10)=huge(icount(10))
6082 ir=ipg+(ispc-1)*(
napgrp+1)
6089 icount(1)=icount(1)+in
6090 icount(2)=icount(2)+in*ku
6096 icount(3)=icount(3)+1
6097 icount(4)=min(icount(4),jn)
6098 icount(5)=icount(5)+jn
6099 icount(6)=max(icount(6),jn)
6100 icount(7)=min(icount(7),in)
6101 icount(8)=icount(8)+in
6102 icount(9)=max(icount(9),in)
6103 icount(10)=min(icount(10),in*jn)
6104 icount(11)=icount(11)+in*jn
6105 icount(12)=max(icount(12),in*jn)
6111 WRITE(*,*)
"analysis of sparsity structure"
6112 IF (icount(1) > 0)
THEN
6113 WRITE(*,101)
"rows without compression/blocks ", icount(1)
6114 WRITE(*,101)
" contained elements ", icount(2)
6116 WRITE(*,101)
"number of block matrices ", icount(3)
6117 avg=real(icount(5),mps)/real(icount(3),mps)
6118 WRITE(*,101)
"number of columns (min,mean,max) ", icount(4), avg, icount(6)
6119 avg=real(icount(8),mps)/real(icount(3),mps)
6120 WRITE(*,101)
"number of rows (min,mean,max) ", icount(7), avg, icount(9)
6121 avg=real(icount(11),mps)/real(icount(3),mps)
6122 WRITE(*,101)
"number of elements (min,mean,max) ", icount(10), avg, icount(12)
6123101
FORMAT(2x,a34,i10,f10.3,i10)
6142 INTEGER(mpi),
INTENT(IN) :: n
6143 REAL(mpd),
INTENT(IN) :: x(n)
6144 REAL(mpd),
INTENT(OUT) :: b(n)
6149 CALL peend(24,
'Aborted, vector/matrix size mismatch')
6150 stop
'AVPROD: mismatched vector and matrix'
6181 INTEGER(mpi) :: ispc
6182 INTEGER(mpi) :: item1
6183 INTEGER(mpi) :: item2
6184 INTEGER(mpi) :: itemc
6185 INTEGER(mpi) :: jtem
6186 INTEGER(mpi) :: jtemn
6189 INTEGER(mpi),
INTENT(IN) :: itema
6190 INTEGER(mpi),
INTENT(IN) :: itemb
6191 INTEGER(mpl),
INTENT(OUT) :: ij
6192 INTEGER(mpi),
INTENT(OUT) :: lr
6193 INTEGER(mpi),
INTENT(OUT) :: iprc
6204 item1=max(itema,itemb)
6205 item2=min(itema,itemb)
6206 IF(item2 <= 0.OR.item1 >
napgrp)
RETURN
6209 outer:
DO ispc=1,
nspc
6220 IF(ku < kl) cycle outer
6225 IF(item2 >= jtem.AND.item2 < jtemn)
THEN
6231 IF(item2 < jtem)
THEN
6233 ELSE IF(item2 >= jtemn)
THEN
6248 IF(ku < kl) cycle outer
6252 IF(itemc == jtem)
EXIT
6253 IF(itemc < jtem)
THEN
6255 ELSE IF(itemc > jtem)
THEN
6283 INTEGER(mpi),
INTENT(IN) :: itema
6284 INTEGER(mpi),
INTENT(IN) :: itemb
6309 INTEGER(mpi) :: item1
6310 INTEGER(mpi) :: item2
6311 INTEGER(mpi) :: ipg1
6312 INTEGER(mpi) :: ipg2
6314 INTEGER(mpi) :: iprc
6316 INTEGER(mpi),
INTENT(IN) :: itema
6317 INTEGER(mpi),
INTENT(IN) :: itemb
6319 INTEGER(mpl) ::
ijadd
6323 item1=max(itema,itemb)
6324 item2=min(itema,itemb)
6326 IF(item2 <= 0.OR.item1 >
nagb)
RETURN
6327 IF(item1 == item2)
THEN
6336 CALL ijpgrp(ipg1,ipg2,ij,lr,iprc)
6358 INTEGER(mpi) :: item1
6359 INTEGER(mpi) :: item2
6360 INTEGER(mpi) :: jtem
6362 INTEGER(mpi),
INTENT(IN) :: itema
6363 INTEGER(mpi),
INTENT(IN) :: itemb
6372 item1=max(itema,itemb)
6373 item2=min(itema,itemb)
6375 IF(item2 <= 0.OR.item1 >
nagb)
RETURN
6383 print *,
' IJCSR3 empty list ', item1, item2, ks, ke
6384 CALL peend(23,
'Aborted, bad matrix index')
6385 stop
'ijcsr3: empty list'
6390 IF(item1 == jtem)
EXIT
6391 IF(item1 < jtem)
THEN
6398 print *,
' IJCSR3 not found ', item1, item2, ks, ke
6399 CALL peend(23,
'Aborted, bad matrix index')
6400 stop
'ijcsr3: not found'
6416 INTEGER(mpi) :: item1
6417 INTEGER(mpi) :: item2
6421 INTEGER(mpl) ::
ijadd
6424 INTEGER(mpi),
INTENT(IN) :: itema
6425 INTEGER(mpi),
INTENT(IN) :: itemb
6430 item1=max(itema,itemb)
6431 item2=min(itema,itemb)
6432 IF(item2 <= 0.OR.item1 >
nagb)
RETURN
6441 ij=
ijadd(item1,item2)
6444 ELSE IF (ij < 0)
THEN
6476 INTEGER(mpi) :: ichunk
6480 INTEGER(mpi) :: ispc
6488 INTEGER(mpl) :: ijadd
6510 ir=ipg+(ispc-1)*(
napgrp+1)
6560 REAL(mps),
INTENT(IN) :: deltat
6561 INTEGER(mpi),
INTENT(OUT) :: minut
6562 INTEGER(mpi),
INTENT(OUT):: nhour
6563 REAL(mps),
INTENT(OUT):: secnd
6564 INTEGER(mpi) :: nsecd
6567 nsecd=nint(deltat,
mpi)
6569 minut=nsecd/60-60*nhour
6570 secnd=deltat-60*(minut+60*nhour)
6606 INTEGER(mpi),
INTENT(IN) :: item
6610 INTEGER(mpl) :: length
6611 INTEGER(mpl),
PARAMETER :: four = 4
6615 IF(item <= 0)
RETURN
6636 IF(j == 0)
EXIT inner
6658 IF (
lvllog > 1)
WRITE(
lunlog,*)
'INONE: array increased to', &
6679 INTEGER(mpi) :: iprime
6680 INTEGER(mpi) :: nused
6681 LOGICAL :: finalUpdate
6682 INTEGER(mpl) :: oldLength
6683 INTEGER(mpl) :: newLength
6684 INTEGER(mpl),
PARAMETER :: four = 4
6685 INTEGER(mpi),
DIMENSION(:,:),
ALLOCATABLE :: tempArr
6686 INTEGER(mpl),
DIMENSION(:),
ALLOCATABLE :: tempVec
6690 IF(finalupdate)
THEN
6699 CALL mpalloc(temparr,four,oldlength,
'INONE: temp array')
6701 CALL mpalloc(tempvec,oldlength,
'INONE: temp vector')
6725 IF(j == 0)
EXIT inner
6726 IF(j == i) cycle outer
6730 IF(.NOT.finalupdate)
RETURN
6735 WRITE(
lunlog,*)
'INONE: array reduced to',newlength,
' words'
6759 IF(j == 0)
EXIT inner
6760 IF(j == i) cycle outer
6777 INTEGER(mpi),
INTENT(IN) :: n
6778 INTEGER(mpi) :: nprime
6779 INTEGER(mpi) :: nsqrt
6784 IF(mod(nprime,2) == 0) nprime=nprime+1
6787 nsqrt=int(sqrt(real(nprime,
mps)),
mpi)
6789 IF(i*(nprime/i) == nprime) cycle outer
6811 INTEGER(mpi) :: idum
6813 INTEGER(mpi) :: indab
6814 INTEGER(mpi) :: itgbi
6815 INTEGER(mpi) :: itgbl
6816 INTEGER(mpi) :: ivgbi
6818 INTEGER(mpi) :: jgrp
6819 INTEGER(mpi) :: lgrp
6821 INTEGER(mpi) :: nc31
6823 INTEGER(mpi) :: nwrd
6824 INTEGER(mpi) :: inone
6829 INTEGER(mpl) :: length
6830 INTEGER(mpl) :: rows
6834 WRITE(
lunlog,*)
'LOOP1: starting'
6857 WRITE(
lunlog,*)
'LOOP1: reading data files'
6870 WRITE(*,*)
'Read all binary data files:'
6872 CALL hmpldf(1,
'Number of words/record in binary file')
6873 CALL hmpdef(8,0.0,60.0,
'not_stored data per record')
6903 WRITE(
lunlog,*)
'LOOP1: reading data files again'
6915 CALL peend(21,
'Aborted, no labels/parameters defined')
6916 stop
'LOOP1: no labels/parameters defined'
6921 ' is total number NTGB of labels/parameters'
6923 CALL hmpldf(2,
'Number of entries per label')
6967 WRITE(
lunlog,*)
'LOOP1:',
npresg,
' is number of pre-sigmas'
6968 WRITE(*,*)
'LOOP1:',
npresg,
' is number of pre-sigmas'
6969 IF(
npresg == 0)
WRITE(*,*)
'Warning: no pre-sigmas defined'
6991 WRITE(
lunlog,*)
'LOOP1:',
nvgb,
' is number NVGB of variable parameters'
6996 WRITE(
lunlog,*)
'LOOP1: counting records, NO iteration of entries cut !'
7002 CALL hmpdef(15,0.0,120.0,
'Number of parameters per group')
7037 ' is total number NTPGRP of label/parameter groups'
7053 WRITE(*,112)
' Default pre-sigma =',
regpre, &
7054 ' (if no individual pre-sigma defined)'
7055 WRITE(*,*)
'Pre-sigma factor is',
regula
7058 WRITE(*,*)
'No regularization will be done'
7060 WRITE(*,*)
'Regularization will be done, using factor',
regula
7064 CALL peend(22,
'Aborted, no variable global parameters')
7065 stop
'... no variable global parameters'
7072 IF(presg > 0.0)
THEN
7074 ELSE IF(presg == 0.0.AND.
regpre > 0.0)
THEN
7075 prewt=1.0/real(
regpre**2,mpd)
7091 WRITE(*,101)
'NTGB',
ntgb,
'total number of parameters'
7092 WRITE(*,101)
'NVGB',
nvgb,
'number of variable parameters'
7096 WRITE(*,101)
'Too many variable parameters for diagonalization, fallback is inversion'
7104 WRITE(*,101)
' NREC',
nrec,
'number of records'
7105 IF (
nrecd > 0)
WRITE(*,101)
' NRECD',
nrec,
'number of records containing doubles'
7106 WRITE(*,101)
' NEQN',
neqn,
'number of equations (measurements)'
7107 WRITE(*,101)
' NEGB',
negb,
'number of equations with global parameters'
7108 WRITE(*,101)
' NDGB',
ndgb,
'number of global derivatives'
7110 WRITE(*,101)
'MREQENF',
mreqenf,
'required number of entries (eqns in binary files)'
7112 WRITE(*,101)
'MREQENF',
mreqenf,
'required number of entries (recs in binary files)'
7115 WRITE(*,101)
'ITEREN',
iteren,
'iterate cut for parameters with less entries'
7116 WRITE(*,101)
'MREQENA',
mreqena,
'required number of entries (from accepted fits)'
7117 IF (
mreqpe > 1)
WRITE(*,101) &
7118 'MREQPE',
mreqpe,
'required number of pair entries'
7119 IF (
msngpe >= 1)
WRITE(*,101) &
7120 'MSNGPE',
msngpe,
'max pair entries single prec. storage'
7121 WRITE(*,101)
'NTGB',
ntgb,
'total number of parameters'
7122 WRITE(*,101)
'NVGB',
nvgb,
'number of variable parameters'
7125 WRITE(*,*)
'Global parameter labels:'
7132 mqi=((mqi-20)/20)*20+1
7140 WRITE(8,101)
' NREC',
nrec,
'number of records'
7141 IF (
nrecd > 0)
WRITE(8,101)
' NRECD',
nrec,
'number of records containing doubles'
7142 WRITE(8,101)
' NEQN',
neqn,
'number of equations (measurements)'
7143 WRITE(8,101)
' NEGB',
negb,
'number of equations with global parameters'
7144 WRITE(8,101)
' NDGB',
ndgb,
'number of global derivatives'
7146 WRITE(8,101)
'MREQENF',
mreqenf,
'required number of entries (eqns in binary files)'
7148 WRITE(8,101)
'MREQENF',
mreqenf,
'required number of entries (recs in binary files)'
7151 WRITE(8,101)
'ITEREN',
iteren,
'iterate cut for parameters with less entries'
7152 WRITE(8,101)
'MREQENA',
mreqena,
'required number of entries (from accepted fits)'
7154 WRITE(
lunlog,*)
'LOOP1: ending'
7158101
FORMAT(1x,a8,
' =',i14,
' = ',a)
7174 INTEGER(mpi) :: ibuf
7176 INTEGER(mpi) :: indab
7182 INTEGER(mpi) :: nc31
7184 INTEGER(mpi) :: nlow
7186 INTEGER(mpi) :: nwrd
7188 INTEGER(mpl) :: length
7189 INTEGER(mpl),
DIMENSION(:),
ALLOCATABLE :: newCounter
7194 WRITE(
lunlog,*)
'LOOP1: iterating'
7196 WRITE(*,*)
'LOOP1: iterating'
7199 CALL mpalloc(newcounter,length,
'new entries counter')
7223 CALL isjajb(nst,ist,ja,jb,jsp)
7224 IF(ja == 0.AND.jb == 0)
EXIT
7230 IF(ij == -2) nlow=nlow+1
7235 newcounter(ij)=newcounter(ij)+1
7264 WRITE(
lunlog,*)
'LOOP1:',
nvgb,
' is number NVGB of variable parameters'
7294 INTEGER(mpi) :: ibuf
7295 INTEGER(mpi) :: icblst
7296 INTEGER(mpi) :: icboff
7297 INTEGER(mpi) :: icgb
7298 INTEGER(mpi) :: icgrp
7299 INTEGER(mpi) :: icount
7300 INTEGER(mpi) :: iext
7301 INTEGER(mpi) :: ihis
7305 INTEGER(mpi) :: ioff
7306 INTEGER(mpi) :: ipoff
7307 INTEGER(mpi) :: iproc
7308 INTEGER(mpi) :: irecmm
7310 INTEGER(mpi) :: itgbi
7311 INTEGER(mpi) :: itgbij
7312 INTEGER(mpi) :: itgbik
7313 INTEGER(mpi) :: ivgbij
7314 INTEGER(mpi) :: ivgbik
7315 INTEGER(mpi) :: ivpgrp
7319 INTEGER(mpi) :: jcgrp
7320 INTEGER(mpi) :: jext
7321 INTEGER(mpi) :: jcgb
7322 INTEGER(mpi) :: jrec
7324 INTEGER(mpi) :: joff
7326 INTEGER(mpi) :: kcgrp
7327 INTEGER(mpi) :: kfile
7329 INTEGER(mpi) :: label
7330 INTEGER(mpi) :: labelf
7331 INTEGER(mpi) :: labell
7332 INTEGER(mpi) :: lvpgrp
7335 INTEGER(mpi) :: maeqnf
7336 INTEGER(mpi) :: nall
7337 INTEGER(mpi) :: naeqna
7338 INTEGER(mpi) :: naeqnf
7339 INTEGER(mpi) :: naeqng
7340 INTEGER(mpi) :: npdblk
7341 INTEGER(mpi) :: nc31
7342 INTEGER(mpi) :: ncachd
7343 INTEGER(mpi) :: ncachi
7344 INTEGER(mpi) :: ncachr
7345 INTEGER(mpi) :: ncon
7348 INTEGER(mpi) :: ndfmax
7349 INTEGER(mpi) :: nfixed
7350 INTEGER(mpi) :: nggd
7351 INTEGER(mpi) :: nggi
7352 INTEGER(mpi) :: nmatmo
7353 INTEGER(mpi) :: noff
7354 INTEGER(mpi) :: npair
7355 INTEGER(mpi) :: npar
7356 INTEGER(mpi) :: nparmx
7358 INTEGER(mpi) :: nrece
7359 INTEGER(mpi) :: nrecf
7360 INTEGER(mpi) :: nrecmm
7362 INTEGER(mpi) :: nwrd
7363 INTEGER(mpi) :: inone
7372 INTEGER(mpl):: nblock
7373 INTEGER(mpl):: nbwrds
7374 INTEGER(mpl):: noff8
7375 INTEGER(mpl):: ndimbi
7376 INTEGER(mpl):: ndimsa(4)
7378 INTEGER(mpl):: nnzero
7379 INTEGER(mpl):: matsiz(2)
7380 INTEGER(mpl):: matwords
7381 INTEGER(mpl):: mbwrds
7382 INTEGER(mpl):: length
7385 INTEGER(mpl),
PARAMETER :: two=2
7386 INTEGER(mpi) :: maxGlobalPar = 0
7387 INTEGER(mpi) :: maxLocalPar = 0
7388 INTEGER(mpi) :: maxEquations = 0
7390 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: vecConsGroupList
7391 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: vecConsGroupIndex
7392 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: vecPairedParGroups
7393 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: vecBlockCounts
7396 SUBROUTINE ndbits(npgrp,ndims,nsparr,ihst)
7398 INTEGER(mpi),
DIMENSION(:),
INTENT(IN) :: npgrp
7399 INTEGER(mpl),
DIMENSION(4),
INTENT(OUT) :: ndims
7400 INTEGER(mpl),
DIMENSION(:,:),
INTENT(OUT) :: nsparr
7401 INTEGER(mpi),
INTENT(IN) :: ihst
7403 SUBROUTINE ckbits(npgrp,ndims)
7405 INTEGER(mpi),
DIMENSION(:),
INTENT(IN) :: npgrp
7406 INTEGER(mpl),
DIMENSION(4),
INTENT(OUT) :: ndims
7408 SUBROUTINE spbits(npgrp,nsparr,nsparc)
7410 INTEGER(mpi),
DIMENSION(:),
INTENT(IN) :: npgrp
7411 INTEGER(mpl),
DIMENSION(:,:),
INTENT(IN) :: nsparr
7412 INTEGER(mpi),
DIMENSION(:),
INTENT(OUT) :: nsparc
7414 SUBROUTINE gpbmap(ngroup,npgrp,npair)
7416 INTEGER(mpi),
INTENT(IN) :: ngroup
7417 INTEGER(mpi),
DIMENSION(:,:),
INTENT(IN) :: npgrp
7418 INTEGER(mpi),
DIMENSION(:),
INTENT(OUT) :: npair
7420 SUBROUTINE ggbmap(ipgrp,npair,npgrp)
7422 INTEGER(mpi),
INTENT(IN) :: ipgrp
7423 INTEGER(mpi),
INTENT(OUT) :: npair
7424 INTEGER(mpi),
DIMENSION(:),
INTENT(OUT) :: npgrp
7426 SUBROUTINE pbsbits(npgrp,ibsize,nnzero,nblock,nbkrow)
7428 INTEGER(mpi),
DIMENSION(:),
INTENT(IN) :: npgrp
7429 INTEGER(mpi),
INTENT(IN) :: ibsize
7430 INTEGER(mpl),
INTENT(OUT) :: nnzero
7431 INTEGER(mpl),
INTENT(OUT) :: nblock
7432 INTEGER(mpi),
DIMENSION(:),
INTENT(OUT) :: nbkrow
7434 SUBROUTINE pblbits(npgrp,ibsize,nsparr,nsparc)
7436 INTEGER(mpi),
DIMENSION(:),
INTENT(IN) :: npgrp
7437 INTEGER(mpi),
INTENT(IN) :: ibsize
7438 INTEGER(mpl),
DIMENSION(:),
INTENT(IN) :: nsparr
7439 INTEGER(mpl),
DIMENSION(:),
INTENT(OUT) :: nsparc
7441 SUBROUTINE prbits(npgrp,nsparr)
7443 INTEGER(mpi),
DIMENSION(:),
INTENT(IN) :: npgrp
7444 INTEGER(mpl),
DIMENSION(:),
INTENT(OUT) :: nsparr
7446 SUBROUTINE pcbits(npgrp,nsparr,nsparc)
7448 INTEGER(mpi),
DIMENSION(:),
INTENT(IN) :: npgrp
7449 INTEGER(mpl),
DIMENSION(:),
INTENT(IN) :: nsparr
7450 INTEGER(mpl),
DIMENSION(:),
INTENT(OUT) :: nsparc
7460 WRITE(
lunlog,*)
'LOOP2: starting'
7479 WRITE(
lunlog,*)
' Elimination for constraints enforced for solution by decomposition!'
7484 WRITE(
lunlog,*)
' Lagrange multipliers enforced for solution by sparsePARDISO!'
7489 WRITE(
lunlog,*)
' Elimination for constraints with mpqldec enforced (LAPACK only for unpacked storage)!'
7506 noff8=int(
nagb,mpl)*int(
nagb-1,mpl)/2
7542 print *,
' Variable parameter groups ',
nvpgrp
7590 CALL mpalloc(vecconsgrouplist,length,
'constraint group list')
7591 CALL mpalloc(vecconsgroupindex,length,
'constraint group index')
7601 WRITE(
lunlog,*)
'LOOP2: start event reading'
7641 WRITE(*,*)
'Record number ',
nrec,
' from file ',kfile
7642 IF (wgh /= 1.0)
WRITE(*,*)
' weight ',wrec
7646 CALL isjajb(nst,ist,ja,jb,jsp)
7650 IF(nda ==
mdebg2+1)
WRITE(*,*)
'... and more data'
7655 WRITE(*,*)
'Local derivatives:'
7657107
FORMAT(6(i3,g12.4))
7659 WRITE(*,*)
'Global derivatives:'
7662108
FORMAT(3i11,g12.4)
7665 WRITE(*,*)
'total_par_label __label__ var_par_index derivative'
7680 CALL isjajb(nst,ist,ja,jb,jsp)
7681 IF(ja == 0.AND.jb == 0)
EXIT
7717 IF (kcgrp == jcgrp) cycle
7728 IF (vecconsgroupindex(k) == 0)
THEN
7731 vecconsgrouplist(icgrp)=k
7746 IF (vecconsgroupindex(k) < icount)
THEN
7748 vecconsgroupindex(k)=icount
7766 IF (nfixed > 0) naeqnf=naeqnf+1
7769 IF(ja /= 0.AND.jb /= 0)
THEN
7778 IF (naeqnf > maeqnf) nrecf=nrecf+1
7782 maxglobalpar=max(
nagbn,maxglobalpar)
7783 maxlocalpar=max(
nalcn,maxlocalpar)
7784 maxequations=max(
naeqn,maxequations)
7787 dstat(1)=dstat(1)+real((nwrd+2)*2,mpd)
7788 dstat(2)=dstat(2)+real(
nagbn+2,mpd)
7793 vecconsgroupindex(vecconsgrouplist(k))=0
7798 IF (
nagbn == 0)
THEN
7817 IF (ivpgrp /= lvpgrp)
THEN
7872 (irecmm >= nrecmm.OR.irecmm ==
mxrec))
THEN
7873 IF (nmatmo == 0)
THEN
7875 WRITE(*,*)
'Monitoring of sparse matrix construction'
7876 WRITE(*,*)
' records ........ off-diagonal elements ', &
7877 '....... compression memory'
7878 WRITE(*,*)
' non-zero used(double) used', &
7883 gbc=1.0e-9*real((mpi*ndimsa(2)+mpd*ndimsa(3)+mps*ndimsa(4))/mpi*(bit_size(1_mpi)/8),mps)
7884 gbu=1.0e-9*real(((mpi+mpd)*(ndimsa(3)+ndimsa(4)))/mpi*(bit_size(1_mpi)/8),mps)
7886 WRITE(*,1177) irecmm,ndimsa(1),ndimsa(3),ndimsa(4),cpr,gbc
78871177
FORMAT(i9,3i13,f10.2,f11.6)
7888 DO WHILE(irecmm >= nrecmm)
7907 WRITE(
lunlog,*)
'LOOP2: event reading ended - end of data'
7909 dstat(k)=dstat(k)/real(
nrec,mpd)
7923 CALL mpalloc(vecpairedpargroups,length,
'paired global parameter groups (I)')
7925 print *,
' Total parameter groups pairs',
ntpgrp
7928 CALL ggbmap(i,npair,vecpairedpargroups)
7996 IF(ivgbij > 0.AND.ivgbik > 0)
THEN
7998 IF (
mprint > 1)
WRITE(*,*)
'add index pair ',ivgbij,ivgbik
8075 nparmx=max(nparmx,int(rows,mpi))
8097 WRITE(*,*)
' Parameter block', i, ib-ia, jb-ja, labelf, labell
8101 WRITE(
lunlog,*)
'Detected',
npblck,
'(disjoint) parameter blocks, max size ', nparmx
8103 WRITE(*,*)
'Detected',
npblck,
'(disjoint) parameter blocks, max size ', nparmx
8105 WRITE(*,*)
'Using block diagonal storage mode'
8120 IF (
nagb >= 65536)
THEN
8121 noff=int(noff8/1000,mpi)
8131 CALL hmpdef(ihis,0.0,real(
mhispe,mps),
'NDBITS: #off-diagonal elements')
8136 ndgn=ndimsa(3)+ndimsa(4)
8137 matwords=ndimsa(2)+length*4
8152 length=int(npdblk,mpl)
8153 CALL mpalloc(vecblockcounts,length,
'sparse matrix row offsets (CSR3)')
8155 nbwrds=2*int(nblock,mpl)*int(
ipdbsz(i)*
ipdbsz(i)+1,mpl)
8156 IF ((i == 1).OR.(nbwrds < mbwrds))
THEN
8180 IF (
fcache(2) == 0.0)
THEN
8182 fcache(2)=real(dstat(2),mps)
8183 fcache(3)=real(dstat(3),mps)
8204 ncachd=
ncache-ncachr-ncachi
8235 WRITE(lu,101)
'NTGB',
ntgb,
'total number of parameters'
8236 WRITE(lu,102)
'(all parameters, appearing in binary files)'
8237 WRITE(lu,101)
'NVGB',
nvgb,
'number of variable parameters'
8238 WRITE(lu,102)
'(appearing in fit matrix/vectors)'
8239 WRITE(lu,101)
'NAGB',
nagb,
'number of all parameters'
8240 WRITE(lu,102)
'(including Lagrange multiplier or reduced)'
8241 WRITE(lu,101)
'NTPGRP',
ntpgrp,
'total number of parameter groups'
8242 WRITE(lu,101)
'NVPGRP',
nvpgrp,
'number of variable parameter groups'
8243 WRITE(lu,101)
'NFGB',
nfgb,
'number of fit parameters'
8245 WRITE(lu,101)
'MBANDW',
mbandw,
'band width of preconditioner matrix'
8246 WRITE(lu,102)
'(if <0, no preconditioner matrix)'
8248 IF (
nagb >= 65536)
THEN
8249 WRITE(lu,101)
'NOFF/K',noff,
'max number of off-diagonal elements'
8251 WRITE(lu,101)
'NOFF',noff,
'max number of off-diagonal elements'
8254 IF (
nagb >= 65536)
THEN
8255 WRITE(lu,101)
'NDGN/K',ndgn/1000,
'actual number of off-diagonal elements'
8257 WRITE(lu,101)
'NDGN',ndgn,
'actual number of off-diagonal elements'
8260 WRITE(lu,101)
'NCGB',
ncgb,
'number of constraints'
8261 WRITE(lu,101)
'NAGBN',
nagbn,
'max number of global parameters in an event'
8262 WRITE(lu,101)
'NALCN',
nalcn,
'max number of local parameters in an event'
8263 WRITE(lu,101)
'NAEQN',
naeqn,
'max number of equations in an event'
8265 WRITE(lu,101)
'NAEQNA',naeqna,
'number of equations'
8266 WRITE(lu,101)
'NAEQNG',naeqng, &
8267 'number of equations with global parameters'
8268 WRITE(lu,101)
'NAEQNF',naeqnf, &
8269 'number of equations with fixed global parameters'
8270 WRITE(lu,101)
'NRECF',nrecf, &
8271 'number of records with fixed global parameters'
8274 WRITE(lu,101)
'NRECE',nrece, &
8275 'number of records without variable parameters'
8278 WRITE(lu,101)
'NCACHE',
ncache,
'number of words for caching'
8279 WRITE(lu,111) (
fcache(k)*100.0,k=1,3)
8280111
FORMAT(22x,
'cache splitting ',3(f6.1,
' %'))
8285 WRITE(lu,*)
'Solution method and matrix-storage mode:'
8287 WRITE(lu,*)
' METSOL = 1: matrix inversion'
8288 ELSE IF(
metsol == 2)
THEN
8289 WRITE(lu,*)
' METSOL = 2: diagonalization'
8290 ELSE IF(
metsol == 3)
THEN
8291 WRITE(lu,*)
' METSOL = 3: decomposition'
8292 ELSE IF(
metsol == 4)
THEN
8293 WRITE(lu,*)
' METSOL = 4: MINRES (rtol',
mrestl,
')'
8294 ELSE IF(
metsol == 5)
THEN
8295 WRITE(lu,*)
' METSOL = 5: MINRES-QLP (rtol',
mrestl,
')'
8296 ELSE IF(
metsol == 6)
THEN
8297 WRITE(lu,*)
' METSOL = 6: GMRES'
8299 ELSE IF(
metsol == 7)
THEN
8300 WRITE(lu,*)
' METSOL = 7: LAPACK factorization'
8301 ELSE IF(
metsol == 8)
THEN
8302 WRITE(lu,*)
' METSOL = 8: LAPACK factorization'
8304 ELSE IF(
metsol == 9)
THEN
8305 WRITE(lu,*)
' METSOL = 9: Intel oneMKL PARDISO'
8309 WRITE(lu,*)
' with',
mitera,
' iterations'
8311 WRITE(lu,*)
' MATSTO = 0: unpacked symmetric matrix, ',
'n*n elements'
8312 ELSE IF(
matsto == 1)
THEN
8313 WRITE(lu,*)
' MATSTO = 1: full symmetric matrix, ',
'(n*n+n)/2 elements'
8314 ELSE IF(
matsto == 2)
THEN
8315 WRITE(lu,*)
' MATSTO = 2: sparse matrix (custom)'
8316 ELSE IF(
matsto == 3)
THEN
8318 WRITE(lu,*)
' MATSTO = 3: sparse matrix (upper triangle, CSR3)'
8320 WRITE(lu,*)
' MATSTO = 3: sparse matrix (upper triangle, BSR3)'
8321 WRITE(lu,*)
' block size',
matbsz
8325 WRITE(lu,*)
' block diagonal with',
npblck,
' blocks'
8327 IF(
mextnd>0)
WRITE(lu,*)
' with extended storage'
8328 IF(
dflim /= 0.0)
THEN
8329 WRITE(lu,103)
'Convergence assumed, if expected dF <',
dflim
8334 WRITE(lu,*)
'Constraints handled by elimination with LAPACK'
8336 WRITE(lu,*)
'Constraints handled by elimination'
8339 WRITE(lu,*)
'Constraints handled by Lagrange multipliers'
8346 CALL peend(28,
'Aborted, no local parameters')
8347 stop
'LOOP2: stopping due to missing local parameters'
8368105
FORMAT(
' Constants C1, C2 for Wolfe conditions:',g12.4,
', ',g12.4)
8376 matwords=(length+
nagb+1)*2
8383 ELSE IF (
matsto == 2)
THEN
8384 matsiz(1)=ndimsa(3)+
nagb
8399 IF (icblst > icboff)
THEN
8405 nall = npar;
IF (
icelim <= 0) nall=npar+ncon
8409 matsiz(1)=matsiz(1)+k
8411 matsiz(1)=matsiz(1)+nall
8416 matwords=matwords+matsiz(1)*2+matsiz(2)
8422 IF (matwords < 250000)
THEN
8423 WRITE(*,*)
'Size of global matrix: < 1 MB'
8425 WRITE(*,*)
'Size of global matrix:',int(real(matwords,mps)*4.0e-6,mpi),
' MB'
8431 WRITE(
lunlog,*)
' Cut values of Chi^2/Ndf and Chi2,'
8432 WRITE(
lunlog,*)
' corresponding to 2 and 3 standard deviations'
8433 WRITE(
lunlog,*)
' Ndf Chi^2/Ndf(2) Chi^2(2) ', &
8434 ' Chi^2/Ndf(3) Chi^2(3)'
8437 IF(ndf >
naeqn)
EXIT
8440 ELSE IF(ndf < 20)
THEN
8442 ELSE IF(ndf < 100)
THEN
8444 ELSE IF(ndf < 200)
THEN
8451 WRITE(
lunlog,106) ndf,chin2,chin2*real(ndf,mps),chin3, chin3*real(ndf,mps)
8454 WRITE(
lunlog,*)
'LOOP2: ending'
8458 IF (
ncgbe /= 0)
THEN
8461 WRITE(*,199)
'WarningWarningWarningWarningWarningWarningWarningWarningWar'
8462 WRITE(*,199)
'arningWarningWarningWarningWarningWarningWarningWarningWarn'
8463 WRITE(*,199)
'rningWarningWarningWarningWarningWarningWarningWarningWarni'
8464 WRITE(*,199)
'ningWarningWarningWarningWarningWarningWarningWarningWarnin'
8465 WRITE(*,199)
'ingWarningWarningWarningWarningWarningWarningWarningWarning'
8466 WRITE(*,199)
'ngWarningWarningWarningWarningWarningWarningWarningWarningW'
8467 WRITE(*,199)
'gWarningWarningWarningWarningWarningWarningWarningWarningWa'
8469 WRITE(*,*)
' Number of empty constraints =',abs(
ncgbe),
', should be 0'
8470 WRITE(*,*)
' => please check constraint definition, mille data'
8472 WRITE(*,199)
'WarningWarningWarningWarningWarningWarningWarningWarningWar'
8473 WRITE(*,199)
'arningWarningWarningWarningWarningWarningWarningWarningWarn'
8474 WRITE(*,199)
'rningWarningWarningWarningWarningWarningWarningWarningWarni'
8475 WRITE(*,199)
'ningWarningWarningWarningWarningWarningWarningWarningWarnin'
8476 WRITE(*,199)
'ingWarningWarningWarningWarningWarningWarningWarningWarning'
8477 WRITE(*,199)
'ngWarningWarningWarningWarningWarningWarningWarningWarningW'
8478 WRITE(*,199)
'gWarningWarningWarningWarningWarningWarningWarningWarningWa'
8483101
FORMAT(1x,a8,
' =',i14,
' = ',a)
8485103
FORMAT(1x,a,g12.4)
8486106
FORMAT(i6,2(3x,f9.3,f12.1,3x))
8501 INTEGER(mpi) :: imed
8504 INTEGER(mpi) :: nent
8505 INTEGER(mpi),
DIMENSION(measBins) :: isuml
8506 INTEGER(mpi),
DIMENSION(measBins) :: isums
8510 INTEGER(mpl) :: ioff
8513 DATA lfirst /.true./
8526 WRITE(
lunmon,
'(A)')
'*** Normalized residuals grouped by first global label (per local fit cycle) ***'
8528 WRITE(
lunmon,
'(A)')
'*** Pulls grouped by first global label (per local fit cycle) ***'
8530 WRITE(
lunmon,
'(A)')
'! LFC Label Entries Median RMS(MAD) <error>'
8549 IF (2*isuml(j) > nent)
EXIT
8553 IF (isuml(j) > isuml(j-1)) amed=amed+real(nent-2*isuml(j-1),mps)/real(2*isuml(j)-2*isuml(j-1),mps)
8566 isums(j)=isums(j)+isums(j-1)
8570 IF (2*isums(j) > nent)
EXIT
8573 IF (isums(j) > isums(j-1)) amad=amad+real(nent-2*isums(j-1),mps)/real(2*isums(j)-2*isums(j-1),mps)
8587110
FORMAT(i5,2i10,3g14.5)
8602 INTEGER(mpi) :: ioff
8603 INTEGER(mpi) :: ipar0
8604 INTEGER(mpi) :: ncon
8605 INTEGER(mpi) :: npar
8606 INTEGER(mpi) :: nextra
8608 INTEGER :: nbopt, nboptx, ILAENV
8611 INTEGER(mpl),
INTENT(IN) :: msize(2)
8613 INTEGER(mpl) :: length
8614 INTEGER(mpl) :: nwrdpc
8615 INTEGER(mpl),
PARAMETER :: three = 3
8628 CALL mpalloc(
aux,length,
' local fit scratch array: aux')
8629 CALL mpalloc(
vbnd,length,
' local fit scratch array: vbnd')
8630 CALL mpalloc(
vbdr,length,
' local fit scratch array: vbdr')
8633 CALL mpalloc(
vbk,length,
' local fit scratch array: vbk')
8636 CALL mpalloc(
vzru,length,
' local fit scratch array: vzru')
8637 CALL mpalloc(
scdiag,length,
' local fit scratch array: scdiag')
8638 CALL mpalloc(
scflag,length,
' local fit scratch array: scflag')
8639 CALL mpalloc(
ibandh,2*length,
' local fit band width hist.: ibandh')
8653 length=4*
ncblck;
IF(ncon == 0) length=0
8662 nwrdpc=nwrdpc+length
8674 length=(ncon*ncon+ncon)/2
8699 length=length+npar+nextra
8705 IF(
mbandw == 0) length=length+1
8713 nwrdpc=nwrdpc+2*length
8714 IF (nwrdpc > 250000)
THEN
8716 WRITE(*,*)
'Size of preconditioner matrix:',int(real(nwrdpc,mps)*4.0e-6,mpi),
' MB'
8738 CALL peend(23,
'Aborted, bad matrix index (will exceed 32bit)')
8739 stop
'vmprep: bad index (matrix to large for diagonalization)'
8761 nbopt = ilaenv( 1_mpl,
'DSYTRF',
'U', int(
nagb,mpl), int(
nagb,mpl), -1_mpl, -1_mpl )
8763 print *,
'LAPACK optimal block size for DSYTRF:', nbopt
8764 lplwrk=length*int(nbopt,mpl)
8772 nbopt = ilaenv( 1_mpl,
'DORMQL',
'RN', int(npar,mpl), int(npar,mpl), int(ncon,mpl), int(npar,mpl) )
8773 IF (int(npar,mpl)*int(nbopt,mpl) >
lplwrk)
THEN
8774 lplwrk=int(npar,mpl)*int(nbopt,mpl)
8779 print *,
'LAPACK optimal block size for DORMQL:', nboptx
8798 INTEGER(mpi) :: icoff
8799 INTEGER(mpi) :: ipoff
8802 INTEGER(mpi) :: ncon
8803 INTEGER(mpi) :: nfit
8804 INTEGER(mpi) :: npar
8805 INTEGER(mpi) :: nrank
8806 INTEGER(mpl) :: imoff
8807 INTEGER(mpl) :: ioff1
8825 WRITE(
lunlog,*)
'Shrinkage of global matrix (A->Q^t*A*Q)'
8840 nfit=npar+ncon;
IF (
icelim > 0) nfit=npar-ncon
8842 IF(nfit < npar)
THEN
8860 WRITE(
lunlog,*)
'Inversion of global matrix (A->A^-1)'
8867 IF(nfit /= nrank)
THEN
8868 WRITE(*,*)
'Warning: the rank defect of the symmetric',nfit, &
8869 '-by-',nfit,
' matrix is ',nfit-nrank,
' (should be zero).'
8870 WRITE(lun,*)
'Warning: the rank defect of the symmetric',nfit, &
8871 '-by-',nfit,
' matrix is ',nfit-nrank,
' (should be zero).'
8874 WRITE(*,*)
' --> enforcing SUBITO mode'
8875 WRITE(lun,*)
' --> enforcing SUBITO mode'
8877 ELSE IF(
ndefec == 0)
THEN
8879 WRITE(lun,*)
'No rank defect of the symmetric matrix'
8881 WRITE(lun,*)
'No rank defect of the symmetric block', ib,
' of size', npar
8892 IF(nfit < npar)
THEN
8911 INTEGER(mpi) :: icoff
8912 INTEGER(mpi) :: ipoff
8915 INTEGER(mpi) :: ncon
8916 INTEGER(mpi) :: nfit
8917 INTEGER(mpi) :: npar
8918 INTEGER(mpi) :: nrank
8919 INTEGER(mpl) :: imoff
8920 INTEGER(mpl) :: ioff1
8935 WRITE(
lunlog,*)
'Shrinkage of global matrix (A->Q^t*A*Q)'
8949 nfit=npar+ncon;
IF (
icelim > 0) nfit=npar-ncon
8951 IF(nfit < npar)
THEN
8969 WRITE(
lunlog,*)
'Decomposition of global matrix (A->L*D*L^t)'
8975 IF(nfit /= nrank)
THEN
8976 WRITE(*,*)
'Warning: the rank defect of the symmetric',nfit, &
8977 '-by-',nfit,
' matrix is ',nfit-nrank,
' (should be zero).'
8978 WRITE(lun,*)
'Warning: the rank defect of the symmetric',nfit, &
8979 '-by-',nfit,
' matrix is ',nfit-nrank,
' (should be zero).'
8982 WRITE(*,*)
' --> enforcing SUBITO mode'
8983 WRITE(lun,*)
' --> enforcing SUBITO mode'
8985 ELSE IF(
ndefec == 0)
THEN
8987 WRITE(lun,*)
'No rank defect of the symmetric matrix'
8989 WRITE(lun,*)
'No rank defect of the symmetric block', ib,
' of size', npar
8991 WRITE(lun,*)
' largest diagonal element (LDLt)', evmax
8992 WRITE(lun,*)
' smallest diagonal element (LDLt)', evmin
9001 IF(nfit < npar)
THEN
9023 INTEGER(mpi) :: icoff
9024 INTEGER(mpi) :: ipoff
9027 INTEGER(mpi) :: ncon
9028 INTEGER(mpi) :: nfit
9029 INTEGER(mpi) :: npar
9030 INTEGER(mpl) :: imoff
9031 INTEGER(mpl) :: ioff1
9032 INTEGER(mpi) :: infolp
9052 WRITE(
lunlog,*)
'Shrinkage of global matrix (A->Q^t*A*Q)'
9067 nfit=npar+ncon;
IF (
icelim > 0) nfit=npar-ncon
9069 IF(nfit < npar)
THEN
9086 IF (nfit > npar)
THEN
9089 WRITE(
lunlog,*)
'Factorization of global matrix (A->L*D*L^t)'
9099 WRITE(
lunlog,*)
'Factorization of global matrix (A->L*L^t)'
9103 CALL dpptrf(
'U',int(nfit,mpl),
globalmatd(imoff+1:),infolp)
9110 WRITE(lun,*)
'No rank defect of the symmetric matrix'
9112 WRITE(lun,*)
'No rank defect of the symmetric block', ib,
' of size', npar
9116 WRITE(*,*)
'Warning: factorization of the symmetric',nfit, &
9117 '-by-',nfit,
' failed at index ', infolp
9118 WRITE(lun,*)
'Warning: factorization of the symmetric',nfit, &
9119 '-by-',nfit,
' failed at index ', infolp
9120 CALL peend(29,
'Aborted, factorization of global matrix failed')
9121 stop
'mdptrf: bad matrix'
9126 IF (nfit > npar)
THEN
9129 IF(infolp /= 0) print *,
' DSPTRS failed: ', infolp
9131 CALL dpptrs(
'U',int(nfit,mpl),1_mpl,
globalmatd(imoff+1:),&
9133 IF(infolp /= 0) print *,
' DPPTRS failed: ', infolp
9137 IF(nfit < npar)
THEN
9158 INTEGER(mpi) :: icoff
9159 INTEGER(mpi) :: ipoff
9162 INTEGER(mpi) :: ncon
9163 INTEGER(mpi) :: nfit
9164 INTEGER(mpi) :: npar
9165 INTEGER(mpl) :: imoff
9166 INTEGER(mpl) :: ioff1
9167 INTEGER(mpl) :: iloff
9168 INTEGER(mpi) :: infolp
9189 WRITE(
lunlog,*)
'Shrinkage of global matrix (A->Q^t*A*Q)'
9209 nfit=npar+ncon;
IF (
icelim > 0) nfit=npar-ncon
9211 IF(nfit < npar)
THEN
9215 CALL dtrtrs(
'L',
'T',
'N',int(ncon,mpl),1_mpl,
lapackql(iloff+npar-ncon+1:),int(npar,mpl),&
9217 IF(infolp /= 0) print *,
' DTRTRS failed: ', infolp
9219 CALL dormql(
'L',
'T',int(npar,mpl),1_mpl,int(ncon,mpl),
lapackql(iloff+1:),int(npar,mpl),&
9221 IF(infolp /= 0) print *,
' DORMQL failed: ', infolp
9240 IF (nfit > npar)
THEN
9243 WRITE(
lunlog,*)
'Factorization of global matrix (A->L*D*L^t)'
9247 CALL dsytrf(
'U',int(nfit,mpl),
globalmatd(imoff+1:),int(nfit,mpl),&
9254 WRITE(
lunlog,*)
'Factorization of global matrix (A->L*L^t)'
9258 CALL dpotrf(
'U',int(nfit,mpl),
globalmatd(imoff+1:),int(npar,mpl),infolp)
9265 WRITE(lun,*)
'No rank defect of the symmetric matrix'
9267 WRITE(lun,*)
'No rank defect of the symmetric block', ib,
' of size', npar
9271 WRITE(*,*)
'Warning: factorization of the symmetric',nfit, &
9272 '-by-',nfit,
' failed at index ', infolp
9273 WRITE(lun,*)
'Warning: factorization of the symmetric',nfit, &
9274 '-by-',nfit,
' failed at index ', infolp
9275 CALL peend(29,
'Aborted, factorization of global matrix failed')
9276 stop
'mdutrf: bad matrix'
9281 IF (nfit > npar)
THEN
9282 CALL dsytrs(
'U',int(nfit,mpl),1_mpl,
globalmatd(imoff+1:),int(nfit,mpl),&
9284 IF(infolp /= 0) print *,
' DSYTRS failed: ', infolp
9286 CALL dpotrs(
'U',int(nfit,mpl),1_mpl,
globalmatd(imoff+1:),int(npar,mpl),&
9288 IF(infolp /= 0) print *,
' DPOTRS failed: ', infolp
9292 IF(nfit < npar)
THEN
9297 CALL dormql(
'L',
'N',int(npar,mpl),1_mpl,int(ncon,mpl),
lapackql(iloff+1:),int(npar,mpl),&
9299 IF(infolp /= 0) print *,
' DORMQL failed: ', infolp
9306 iloff=iloff+int(npar,mpl)*int(ncon,mpl)
9328 INTEGER(mpi) :: icboff
9329 INTEGER(mpi) :: icblst
9330 INTEGER(mpi) :: icoff
9331 INTEGER(mpi) :: icfrst
9332 INTEGER(mpi) :: iclast
9333 INTEGER(mpi) :: ipfrst
9334 INTEGER(mpi) :: iplast
9335 INTEGER(mpi) :: ipoff
9338 INTEGER(mpi) :: ncon
9339 INTEGER(mpi) :: npar
9341 INTEGER(mpl) :: imoff
9342 INTEGER(mpl) :: iloff
9343 INTEGER(mpi) :: infolp
9344 INTEGER :: nbopt, ILAENV
9346 REAL(mpd),
INTENT(IN) :: a(mszcon)
9347 REAL(mpd),
INTENT(OUT) :: emin
9348 REAL(mpd),
INTENT(OUT) :: emax
9359 iloff=iloff+int(npar,mpl)*int(ncon,mpl)
9378 DO icb=icboff+1,icboff+icblst
9385 lapackql(iloff+ipfrst:iloff+iplast)=a(imoff+1:imoff+npb)
9402 nbopt = ilaenv( 1_mpl,
'DGEQLF',
'', int(npar,mpl), int(ncon,mpl), int(npar,mpl), -1_mpl )
9403 print *,
'LAPACK optimal block size for DGEQLF:', nbopt
9404 lplwrk=int(ncon,mpl)*int(nbopt,mpl)
9407 CALL dgeqlf(int(npar,mpl),int(ncon,mpl),
lapackql(iloff+1:),int(npar,mpl),&
9409 IF(infolp /= 0) print *,
' DGEQLF failed: ', infolp
9412 iloff=iloff+int(npar,mpl)*int(ncon,mpl)
9415 IF(emax < emin)
THEN
9443 INTEGER(mpi) :: icoff
9444 INTEGER(mpi) :: ipoff
9446 INTEGER(mpi) :: ncon
9447 INTEGER(mpi) :: npar
9448 INTEGER(mpl) :: imoff
9449 INTEGER(mpl) :: iloff
9450 INTEGER(mpi) :: infolp
9451 CHARACTER (LEN=1) :: transr, transl
9453 LOGICAL,
INTENT(IN) :: t
9472 IF(ncon <= 0 ) cycle
9481 DO i=ipoff+1,ipoff+npar
9487 CALL dormql(
'R',transr,int(npar,mpl),int(npar,mpl),int(ncon,mpl),
lapackql(iloff+1:),&
9490 IF(infolp /= 0) print *,
' DORMQL failed: ', infolp
9492 CALL dormql(
'L',transl,int(npar,mpl),int(npar,mpl),int(ncon,mpl),
lapackql(iloff+1:),&
9495 IF(infolp /= 0) print *,
' DORMQL failed: ', infolp
9498 iloff=iloff+int(npar,mpl)*int(ncon,mpl)
9504include
'mkl_pardiso.f90'
9535 TYPE(mkl_pardiso_handle) :: pt(64)
9537 INTEGER(mpl),
PARAMETER :: maxfct =1
9538 INTEGER(mpl),
PARAMETER :: mnum = 1
9539 INTEGER(mpl),
PARAMETER :: nrhs = 1
9541 INTEGER(mpl) :: mtype
9542 INTEGER(mpl) :: phase
9543 INTEGER(mpl) :: error
9544 INTEGER(mpl) :: msglvl
9548 INTEGER(mpl) :: idum(1)
9550 INTEGER(mpl) :: length
9551 INTEGER(mpi) :: nfill
9552 INTEGER(mpi) :: npdblk
9553 REAL(mpd) :: adum(1)
9554 REAL(mpd) :: ddum(1)
9556 INTEGER(mpl) :: iparm(64)
9557 REAL(mpd),
ALLOCATABLE :: b( : )
9558 REAL(mpd),
ALLOCATABLE :: x( : )
9577 WRITE(*,*)
'MSPARDISO: number of rows to fill up = ', nfill
9590 CALL pardiso_64(pt, maxfct, mnum, mtype, phase, int(npdblk,mpl), adum, idum, idum, &
9591 idum, nrhs, iparm, msglvl, ddum, ddum, error)
9592 IF (error /= 0)
THEN
9593 WRITE(lun,*)
'The following ERROR was detected: ', error
9594 WRITE(*,
'(A,2I10)')
' PARDISO release failed (phase, error): ', phase, error
9595 IF (
ipddbg == 0)
WRITE(*,*)
' rerun with "debugPARDISO" for more info'
9596 CALL peend(40,
'Aborted, other error: PARDISO release')
9597 stop
'MSPARDISO: stopping due to error in PARDISO release'
9614 IF (iparm(1) == 0)
WRITE(lun,*)
'PARDISO using defaults '
9615 IF (iparm(43) /= 0)
THEN
9616 WRITE(lun,*)
'PARDISO: computation of the diagonal of inverse matrix not implemented !'
9630 WRITE(
lunlog,*)
'Decomposition of global matrix (A->L*D*L^t)'
9644 WRITE(lun,*)
' iparm(',i,
') =', iparm(i)
9647 CALL pardiso_64(pt, maxfct, mnum, mtype, phase, int(npdblk,mpl),
globalmatd,
csr3rowoffsets,
csr3columnlist, &
9648 idum, nrhs, iparm, msglvl, ddum, ddum, error)
9650 WRITE(lun,*)
'PARDISO reordering completed ... '
9651 WRITE(lun,*)
'PARDISO peak memory required (KB)', iparm(15)
9654 WRITE(lun,*)
' iparm(',i,
') =', iparm(i)
9657 IF (error /= 0)
THEN
9658 WRITE(lun,*)
'The following ERROR was detected: ', error
9659 WRITE(*,
'(A,2I10)')
' PARDISO decomposition failed (phase, error): ', phase, error
9660 IF (
ipddbg == 0)
WRITE(*,*)
' rerun with "debugPARDISO" for more info'
9661 CALL peend(40,
'Aborted, other error: PARDISO reordering')
9662 stop
'MSPARDISO: stopping due to error in PARDISO reordering'
9664 IF (iparm(60) == 0)
THEN
9669 WRITE(lun,*)
'Size (KB) of allocated memory = ',
ipdmem
9670 WRITE(lun,*)
'Number of nonzeros in factors = ',iparm(18)
9671 WRITE(lun,*)
'Number of factorization MFLOPS = ',iparm(19)
9676 CALL pardiso_64(pt, maxfct, mnum, mtype, phase, int(npdblk,mpl),
globalmatd,
csr3rowoffsets,
csr3columnlist, &
9677 idum, nrhs, iparm, msglvl, ddum, ddum, error)
9679 WRITE(lun,*)
'PARDISO factorization completed ... '
9682 WRITE(lun,*)
' iparm(',i,
') =', iparm(i)
9685 IF (error /= 0)
THEN
9686 WRITE(lun,*)
'The following ERROR was detected: ', error
9687 WRITE(*,
'(A,2I10)')
' PARDISO decomposition failed (phase, error): ', phase, error
9688 IF (
ipddbg == 0)
WRITE(*,*)
' rerun with "debugPARDISO" for more info'
9689 CALL peend(40,
'Aborted, other error: PARDISO factorization')
9690 stop
'MSPARDISO: stopping due to error in PARDISO factorization'
9693 IF (iparm(14) > 0) &
9694 WRITE(lun,*)
'Number of perturbed pivots = ',iparm(14)
9695 WRITE(lun,*)
'Number of positive eigenvalues = ',iparm(22)-nfill
9696 WRITE(lun,*)
'Number of negative eigenvalues = ',iparm(23)
9697 ELSE IF (iparm(30) > 0)
THEN
9698 WRITE(lun,*)
'Equation with bad pivot (<=0.) = ',iparm(30)
9707 CALL mpalloc(b,length,
' PARDISO r.h.s')
9708 CALL mpalloc(x,length,
' PARDISO solution')
9713 CALL pardiso_64(pt, maxfct, mnum, mtype, phase, int(npdblk,mpl),
globalmatd,
csr3rowoffsets,
csr3columnlist, &
9714 idum, nrhs, iparm, msglvl, b, x, error)
9719 WRITE(lun,*)
'PARDISO solve completed ... '
9720 IF (error /= 0)
THEN
9721 WRITE(lun,*)
'The following ERROR was detected: ', error
9722 WRITE(*,
'(A,2I10)')
' PARDISO decomposition failed (phase, error): ', phase, error
9723 IF (
ipddbg == 0)
WRITE(*,*)
' rerun with "debugPARDISO" for more info'
9724 CALL peend(40,
'Aborted, other error: PARDISO solve')
9725 stop
'MSPARDISO: stopping due to error in PARDISO solve'
9739 INTEGER(mpi) :: iast
9740 INTEGER(mpi) :: idia
9741 INTEGER(mpi) :: imin
9742 INTEGER(mpl) :: ioff1
9744 INTEGER(mpi) :: last
9746 INTEGER(mpi) :: nmax
9747 INTEGER(mpi) :: nmin
9748 INTEGER(mpi) :: ntop
9770 WRITE(
lunlog,*)
'Shrinkage of global matrix (A->Q^t*A*Q)'
9807 DO WHILE(ntop < nmax)
9811 CALL hmpdef(7,real(nmin,mps),real(ntop,mps),
'log10 of positive eigenvalues')
9822 CALL gmpdef(3,2,
'low-value end of eigenvalues')
9825 CALL gmpxy(3,real(i,mps),evalue)
9839 WRITE(lun,*)
'The first (largest) eigenvalues ...'
9842 WRITE(lun,*)
'The last eigenvalues ... up to',last
9846 WRITE(lun,*)
'The eigenvalues from',
nvgb+1,
' to',
nagb
9850 WRITE(lun,*)
'Log10 + 3 of ',
nfgb,
' eigenvalues in decreasing',
' order'
9851 WRITE(lun,*)
'(for Eigenvalue < 0.001 the value 0.0 is shown)'
9854 'printed for negative eigenvalues'
9857 WRITE(lun,*) last,
' significances: insignificant if ', &
9858 'compatible with N(0,1)'
9887 INTEGER(mpl) :: ioff1
9888 INTEGER(mpl) :: ioff2
9924 INTEGER(mpi) :: istop
9926 INTEGER(mpi) :: itnlim
9928 INTEGER(mpi) :: nout
9929 INTEGER(mpi) :: nrkd
9930 INTEGER(mpi) :: nrkd2
9975 WRITE(lun,*)
'MMINRS: PRECONS ended ', nrkd
9979 globalcorrections, itnlim, nout, rtol, istop, itn, anorm, acond, rnorm, arnorm, ynorm)
9986 WRITE(lun,*)
'MMINRS: EQUDECS ended ', nrkd, nrkd2
9990 globalcorrections, itnlim, nout, rtol, istop, itn, anorm, acond, rnorm, arnorm, ynorm)
9993 globalcorrections, itnlim, nout, rtol, istop, itn, anorm, acond, rnorm, arnorm, ynorm)
10007 IF (
istopa == 0) print *,
'MINRES: istop=0, exact solution x=0.'
10022 INTEGER(mpi) :: istop
10023 INTEGER(mpi) :: itn
10024 INTEGER(mpi) :: itnlim
10025 INTEGER(mpi) :: lun
10026 INTEGER(mpi) :: nout
10027 INTEGER(mpi) :: nrkd
10028 INTEGER(mpi) :: nrkd2
10031 REAL(mpd) :: mxxnrm
10032 REAL(mpd) :: trcond
10042 mxxnrm = real(
nagb,mpd)/sqrt(epsilon(mxxnrm))
10044 trcond = 1.0_mpd/epsilon(trcond)
10045 ELSE IF(
mrmode == 2)
THEN
10075 WRITE(lun,*)
'MMINRS: PRECONS ended ', nrkd
10079 itnlim=itnlim, rtol=rtol, maxxnorm=mxxnrm, trancond=trcond, &
10081 ELSE IF(
mbandw > 0)
THEN
10087 WRITE(lun,*)
'MMINRS: EQUDECS ended ', nrkd, nrkd2
10092 itnlim=itnlim, rtol=rtol, maxxnorm=mxxnrm, trancond=trcond, &
10096 itnlim=itnlim, rtol=rtol, maxxnorm=mxxnrm, trancond=trcond, &
10111 IF (
istopa == 3) print *,
'MINRES: istop=0, exact solution x=0.'
10127 INTEGER(mpi),
INTENT(IN) :: n
10128 REAL(mpd),
INTENT(IN) :: x(n)
10129 REAL(mpd),
INTENT(OUT) :: y(n)
10149 INTEGER(mpi),
INTENT(IN) :: n
10150 REAL(mpd),
INTENT(IN) :: x(n)
10151 REAL(mpd),
INTENT(OUT) :: y(n)
10182 REAL(mps) :: concu2
10183 REAL(mps) :: concut
10184 REAL,
DIMENSION(2) :: ta
10187 INTEGER(mpi) :: iact
10188 INTEGER(mpi) :: iagain
10189 INTEGER(mpi) :: idx
10190 INTEGER(mpi) :: info
10192 INTEGER(mpi) :: ipoff
10193 INTEGER(mpi) :: icoff
10194 INTEGER(mpl) :: ioff
10195 INTEGER(mpi) :: itgbi
10196 INTEGER(mpi) :: ivgbi
10197 INTEGER(mpi) :: jcalcm
10199 INTEGER(mpi) :: labelg
10200 INTEGER(mpi) :: litera
10201 INTEGER(mpi) :: lrej
10202 INTEGER(mpi) :: lun
10203 INTEGER(mpi) :: lunp
10204 INTEGER(mpi) :: minf
10205 INTEGER(mpi) :: mrati
10206 INTEGER(mpi) :: nan
10207 INTEGER(mpi) :: ncon
10208 INTEGER(mpi) :: nfaci
10209 INTEGER(mpi) :: nloopsol
10210 INTEGER(mpi) :: npar
10211 INTEGER(mpi) :: nrati
10212 INTEGER(mpi) :: nrej
10213 INTEGER(mpi) :: nsol
10214 INTEGER(mpi) :: inone
10216 INTEGER(mpi) :: infolp
10217 INTEGER(mpi) :: nfit
10218 INTEGER(mpl) :: imoff
10222 REAL(mpd) :: dratio
10223 REAL(mpd) :: dwmean
10232 LOGICAL :: warnerss
10233 LOGICAL :: warners3
10235 CHARACTER (LEN=7) :: cratio
10236 CHARACTER (LEN=7) :: cfacin
10237 CHARACTER (LEN=7) :: crjrat
10248 WRITE(lunp,*)
'Solution algorithm: '
10249 WRITE(lunp,121)
'=================================================== '
10252 WRITE(lunp,121)
'solution method:',
'matrix inversion'
10253 ELSE IF(
metsol == 2)
THEN
10254 WRITE(lunp,121)
'solution method:',
'diagonalization'
10255 ELSE IF(
metsol == 3)
THEN
10256 WRITE(lunp,121)
'solution method:',
'decomposition'
10257 ELSE IF(
metsol == 4)
THEN
10258 WRITE(lunp,121)
'solution method:',
'minres (Paige/Saunders)'
10259 ELSE IF(
metsol == 5)
THEN
10260 WRITE(lunp,121)
'solution method:',
'minres-qlp (Choi/Paige/Saunders)'
10262 WRITE(lunp,121)
' ',
' using QR factorization'
10263 ELSE IF(
mrmode == 2)
THEN
10264 WRITE(lunp,121)
' ',
' using QLP factorization'
10266 WRITE(lunp,121)
' ',
' using QR and QLP factorization'
10267 WRITE(lunp,123)
'transition condition',
mrtcnd
10269 ELSE IF(
metsol == 6)
THEN
10270 WRITE(lunp,121)
'solution method:', &
10271 'gmres (generalized minimzation of residuals)'
10273 ELSE IF(
metsol == 7)
THEN
10275 WRITE(lunp,121)
'solution method:',
'LAPACK factorization (DSPTRF)'
10277 WRITE(lunp,121)
'solution method:',
'LAPACK factorization (DPPTRF)'
10279 IF(
ilperr == 1)
WRITE(lunp,121)
' ',
'with error calculation (D??TRI)'
10280 ELSE IF(
metsol == 8)
THEN
10282 WRITE(lunp,121)
'solution method:',
'LAPACK factorization (DSYTRF)'
10284 WRITE(lunp,121)
'solution method:',
'LAPACK factorization (DPOTRF)'
10286 IF(
ilperr == 1)
WRITE(lunp,121)
' ',
'with error calculation (D??TRI)'
10288 ELSE IF(
metsol == 9)
THEN
10290 WRITE(lunp,121)
'solution method:',
'Intel oneMKL PARDISO (sparse matrix (CSR3))'
10292 WRITE(lunp,121)
'solution method:',
'Intel oneMKL PARDISO (sparse matrix (BSR3))'
10297 WRITE(lunp,123)
'convergence limit at Delta F=',
dflim
10298 WRITE(lunp,122)
'maximum number of iterations=',
mitera
10301 WRITE(lunp,122)
'matrix recalculation up to ',
matrit,
'. iteration'
10305 WRITE(lunp,121)
'matrix storage:',
'full'
10306 ELSE IF(
matsto == 2)
THEN
10307 WRITE(lunp,121)
'matrix storage:',
'sparse'
10309 WRITE(lunp,122)
'pre-con band-width parameter=',
mbandw
10311 WRITE(lunp,121)
'pre-conditioning:',
'default'
10312 ELSE IF(
mbandw < 0)
THEN
10313 WRITE(lunp,121)
'pre-conditioning:',
'none!'
10314 ELSE IF(
mbandw > 0)
THEN
10316 WRITE(lunp,121)
'pre-conditioning=',
'skyline-matrix (rank preserving)'
10318 WRITE(lunp,121)
'pre-conditioning=',
'band-matrix'
10323 WRITE(lunp,121)
'using pre-sigmas:',
'no'
10326 WRITE(lunp,124)
'pre-sigmas defined for', &
10327 REAL(100*npresg,mps)/
REAL(nvgb,mps),
' % of variable parameters'
10328 WRITE(lunp,123)
'default pre-sigma=',
regpre
10331 WRITE(lunp,121)
'regularization:',
'no'
10333 WRITE(lunp,121)
'regularization:',
'yes'
10334 WRITE(lunp,123)
'regularization factor=',
regula
10338 WRITE(lunp,121)
'Chi square cut equiv 3 st.dev applied'
10339 WRITE(lunp,123)
'... in first iteration with factor',
chicut
10340 WRITE(lunp,123)
'... in second iteration with factor',
chirem
10341 WRITE(lunp,121)
' (reduced by sqrt in next iterations)'
10344 WRITE(lunp,121)
'Scaling of measurement errors applied'
10345 WRITE(lunp,123)
'... factor for "global" measuements',
dscerr(1)
10346 WRITE(lunp,123)
'... factor for "local" measuements',
dscerr(2)
10349 WRITE(lunp,122)
'Down-weighting of outliers in',
lhuber,
' iterations'
10350 WRITE(lunp,123)
'Cut on downweight fraction',
dwcut
10354121
FORMAT(1x,a40,3x,a)
10355122
FORMAT(1x,a40,3x,i0,a)
10356123
FORMAT(1x,a40,2x,e9.2)
10357124
FORMAT(1x,a40,3x,f5.1,a)
10367 stepl =real(stp,mps)
10379 ELSE IF(
metsol == 2)
THEN
10382 ELSE IF(
metsol == 3)
THEN
10385 ELSE IF(
metsol == 4)
THEN
10388 ELSE IF(
metsol == 5)
THEN
10391 ELSE IF(
metsol == 6)
THEN
10403 WRITE(
lunlog,*)
'Checking feasibility of parameters:'
10404 WRITE(*,*)
'Checking feasibility of parameters:'
10405 CALL feasib(concut,iact)
10407 WRITE(*,102) concut
10408 WRITE(*,*)
' parameters are made feasible'
10409 WRITE(
lunlog,102) concut
10410 WRITE(
lunlog,*)
' parameters are made feasible'
10412 WRITE(*,*)
' parameters are feasible (i.e. satisfy constraints)'
10413 WRITE(
lunlog,*)
' parameters are feasible (i.e. satisfy constraints)'
10421 WRITE(*,*)
'Reading files and accumulating vectors/matrices ...'
10425 WRITE(
lunlog,*)
'Reading files and accumulating vectors/matrices ...'
10442 IF(jcalcm+1 /= 0)
THEN
10451 IF(
iterat /= litera)
THEN
10472 IF (iabs(jcalcm) <= 1)
THEN
10485 IF(3*nrej >
nrecal)
THEN
10487 WRITE(*,*)
'Data rejected in previous loop: '
10489 nrejec(0),
' (rank deficit/NaN) ',
nrejec(1),
' (Ndf=0) ', &
10491 WRITE(*,*)
'Too many rejects (>33.3%) - stop'
10492 CALL peend(26,
'Aborted, too many rejects')
10499 IF(abs(
icalcm) == 1)
THEN
10511 ELSE IF(
metsol == 2)
THEN
10513 ELSE IF(
metsol == 3)
THEN
10515 ELSE IF(
metsol == 4)
THEN
10517 ELSE IF(
metsol == 5)
THEN
10519 ELSE IF(
metsol == 6)
THEN
10520 WRITE(*,*)
'... reserved for GMRES (not yet!)'
10523 ELSE IF(
metsol == 7)
THEN
10525 ELSE IF(
metsol == 8)
THEN
10528 ELSE IF(
metsol == 9)
THEN
10542 CALL feasib(concut,iact)
10554 angras=real(db/sqrt(db1*db2),mps)
10555 dbsig=16.0_mpd*sqrt(max(db1,db2))*epsilon(db)
10561 lsflag=lsflag .AND. (db > dbsig)
10568 ELSE IF(
metsol == 2)
THEN
10571 ELSE IF(
metsol == 3)
THEN
10574 ELSE IF(
metsol == 4)
THEN
10577 ELSE IF(
metsol == 5)
THEN
10580 ELSE IF(
metsol == 6)
THEN
10590 IF(db <= -dbsig)
THEN
10591 WRITE(*,*)
'Function not decreasing:',db
10592 IF(db > -1.0e-3_mpd)
THEN
10594 IF (iagain <= 1)
THEN
10595 WRITE(*,*)
'... again matrix calculation'
10599 WRITE(*,*)
'... aborting iterations'
10603 WRITE(*,*)
'... stopping iterations'
10626 IF (
nloopn == nloopsol)
THEN
10632 stepl=real(stp,mps)
10642 WRITE(*,*)
'Result vector containes ', nan,
' NaNs - stop'
10643 CALL peend(25,
'Aborted, result vector contains NaNs')
10650 WRITE(*,*)
'Subito! Exit after first step.'
10655 WRITE(*,*)
'INFO=0 should not happen (line search input err)'
10656 IF (iagain <= 0)
THEN
10661 IF(info < 0 .OR.
nloopn == nloopsol) cycle
10665 CALL feasib(concut,iact)
10666 IF(iact /= 0.OR.
chicut > 1.0)
THEN
10683 WRITE(*,*)
'Data rejected in last loop: '
10685 nrejec(0),
' (rank deficit/NaN) ',
nrejec(1),
' (Ndf=0) ', &
10707 nfit=npar+ncon;
IF (
icelim > 0) nfit=npar-ncon
10708 IF (nfit > npar)
THEN
10711 WRITE(
lunlog,*)
'Inverse of global matrix from LDLt factorization'
10717 IF(infolp /= 0) print *,
' DSPTRI failed: ', infolp
10722 CALL dsytri(
'U',int(nfit,mpl),
globalmatd(imoff+1:),int(nfit,mpl),&
10724 IF(infolp /= 0) print *,
' DSYTRI failed: ', infolp
10730 WRITE(
lunlog,*)
'Inverse of global matrix from LLt factorization'
10735 CALL dpptri(
'U',int(nfit,mpl),
globalmatd(imoff+1:),infolp)
10736 IF(infolp /= 0) print *,
' DPPTRI failed: ', infolp
10740 CALL dpotri(
'U',int(nfit,mpl),
globalmatd(imoff+1:),int(npar,mpl),infolp)
10741 IF(infolp /= 0) print *,
' DPOTRI failed: ', infolp
10758 DO i=npar-ncon+1,npar
10765 WRITE(
lunlog,*)
'Expansion of global matrix (A->Q*A*Q^t)'
10781 catio=real(dratio,mps)
10785 mrati=nint(100.0*catio,mpi)
10789 IF (
nfilw <= 0)
THEN
10790 WRITE(lunp,*)
'Sum(Chi^2)/Sum(Ndf) =',
fvalue
10792 WRITE(lunp,*)
' =',dratio
10794 WRITE(lunp,*)
'Sum(W*Chi^2)/Sum(Ndf)/<W> =',
fvalue
10796 WRITE(lunp,*)
' /',dwmean
10797 WRITE(lunp,*)
' =',dratio
10801 ' with correction for down-weighting ',catio
10822 nrati=nint(10000.0*real(nrej,mps)/real(
nrecal,mps),mpi)
10823 WRITE(crjrat,197) 0.01_mpd*real(nrati,mpd)
10824 nfaci=nint(100.0*sqrt(catio),mpi)
10826 WRITE(cratio,197) 0.01_mpd*real(mrati,mpd)
10827 WRITE(cfacin,197) 0.01_mpd*real(nfaci,mpd)
10830 IF(mrati < 90.OR.mrati > 110) warner=.true.
10831 IF(nrati > 100) warner=.true.
10832 IF(
ncgbe /= 0) warner=.true.
10834 IF(
nalow /= 0) warners=.true.
10836 IF(
nmiss1 /= 0) warnerss=.true.
10837 IF(iagain /= 0) warnerss=.true.
10838 IF(
ndefec /= 0) warnerss=.true.
10839 IF(
ndefpg /= 0) warnerss=.true.
10841 IF(
nrderr /= 0) warners3=.true.
10843 IF(warner.OR.warners.OR.warnerss.Or.warners3)
THEN
10846 WRITE(*,199)
'WarningWarningWarningWarningWarningWarningWarningWarningWar'
10847 WRITE(*,199)
'arningWarningWarningWarningWarningWarningWarningWarningWarn'
10848 WRITE(*,199)
'rningWarningWarningWarningWarningWarningWarningWarningWarni'
10849 WRITE(*,199)
'ningWarningWarningWarningWarningWarningWarningWarningWarnin'
10850 WRITE(*,199)
'ingWarningWarningWarningWarningWarningWarningWarningWarning'
10851 WRITE(*,199)
'ngWarningWarningWarningWarningWarningWarningWarningWarningW'
10852 WRITE(*,199)
'gWarningWarningWarningWarningWarningWarningWarningWarningWa'
10854 IF(mrati < 90.OR.mrati > 110)
THEN
10856 WRITE(*,*)
' Chi^2/Ndf = ',cratio,
' (should be close to 1)'
10857 WRITE(*,*)
' => multiply all input standard ', &
10858 'deviations by factor',cfacin
10861 IF(nrati > 100)
THEN
10863 WRITE(*,*)
' Fraction of rejects =',crjrat,
' %', &
10864 ' (should be far below 1 %)'
10865 WRITE(*,*)
' => please provide correct mille data'
10869 IF(iagain /= 0)
THEN
10871 WRITE(*,*)
' Matrix not positiv definite '// &
10872 '(function not decreasing)'
10873 WRITE(*,*)
' => please provide correct mille data'
10878 WRITE(*,*)
' Rank defect =',
ndefec, &
10879 ' for global matrix, should be 0'
10880 WRITE(*,*)
' => please provide correct mille data'
10885 WRITE(*,*)
' Rank defect for',
ndefpg, &
10886 ' parameter groups, should be 0'
10887 WRITE(*,*)
' => please provide correct mille data'
10892 WRITE(*,*)
' Rank defect =',
nmiss1, &
10893 ' for constraint equations, should be 0'
10894 WRITE(*,*)
' => please correct constraint definition'
10897 IF(
ncgbe /= 0)
THEN
10899 WRITE(*,*)
' Number of empty constraints =',
ncgbe,
', should be 0'
10900 WRITE(*,*)
' => please check constraint definition, mille data'
10903 IF(
nxlow /= 0)
THEN
10905 WRITE(*,*)
' Possible rank defects =',
nxlow,
' for global matrix'
10906 WRITE(*,*)
' (too few accepted entries)'
10907 WRITE(*,*)
' => please check mille data and ENTRIES cut'
10910 IF(
nalow /= 0)
THEN
10912 WRITE(*,*)
' Possible bad elements =',
nalow,
' in global vector'
10913 WRITE(*,*)
' (toos few accepted entries)'
10914 IF(
ipcntr > 0)
WRITE(*,*)
' (indicated in millepede.res by counts<0)'
10915 WRITE(*,*)
' => please check mille data and ENTRIES cut'
10920 WRITE(*,*)
' Binary file(s) with read errors =',
nrderr,
' (treated as EOF)'
10921 WRITE(*,*)
' => please check mille data'
10925 WRITE(*,199)
'WarningWarningWarningWarningWarningWarningWarningWarningWar'
10926 WRITE(*,199)
'arningWarningWarningWarningWarningWarningWarningWarningWarn'
10927 WRITE(*,199)
'rningWarningWarningWarningWarningWarningWarningWarningWarni'
10928 WRITE(*,199)
'ningWarningWarningWarningWarningWarningWarningWarningWarnin'
10929 WRITE(*,199)
'ingWarningWarningWarningWarningWarningWarningWarningWarning'
10930 WRITE(*,199)
'ngWarningWarningWarningWarningWarningWarningWarningWarningW'
10931 WRITE(*,199)
'gWarningWarningWarningWarningWarningWarningWarningWarningWa'
10942 ELSE IF(
metsol == 2)
THEN
10944 ELSE IF(
metsol == 3)
THEN
10950 IF(labelg == 0) cycle
10951 itgbi=inone(labelg)
10954 IF(ivgbi < 0) ivgbi=0
10955 IF(ivgbi == 0) cycle
10964 ELSE IF(
metsol == 6)
THEN
10967 ELSE IF(
metsol == 7)
THEN
10975 CALL peend(4,
'Ended with severe warnings (bad binary file(s))')
10976 ELSE IF (warnerss)
THEN
10977 CALL peend(3,
'Ended with severe warnings (bad global matrix)')
10978 ELSE IF (warners)
THEN
10979 CALL peend(2,
'Ended with severe warnings (insufficient measurements)')
10980 ELSE IF (warner)
THEN
10981 CALL peend(1,
'Ended with warnings (bad measurements)')
10983 CALL peend(0,
'Ended normally')
10986102
FORMAT(
' Call FEASIB with cut=',g10.3)
11004 INTEGER(mpi) :: kfl
11005 INTEGER(mpi) :: kmin
11006 INTEGER(mpi) :: kmax
11007 INTEGER(mpi) :: nrc
11008 INTEGER(mpi) :: nrej
11014 REAL(mpd) :: sumallw
11015 REAL(mpd) :: sumrejw
11017 sumallw=0.; sumrejw=0.;
11026 sumallw=sumallw+real(nrc,mpd)*
wfd(kfl)
11027 sumrejw=sumrejw+real(nrej,mpd)*
wfd(kfl)
11028 frac=real(nrej,mps)/real(nrc,mps)
11029 IF (frac > fmax)
THEN
11033 IF (frac < fmin)
THEN
11040 WRITE(*,
"(' Weighted fraction =',F8.2,' %')") 100.*sumrejw/sumallw
11041 IF (
nfilb > 1)
THEN
11042 WRITE(*,
"(' File with max. fraction ',I6,' :',F8.2,' %')") kmax, 100.*fmax
11043 WRITE(*,
"(' File with min. fraction ',I6,' :',F8.2,' %')") kmin, 100.*fmin
11069 INTEGER(mpi) :: iargc
11072 INTEGER(mpi) :: ierrf
11073 INTEGER(mpi) :: ieq
11074 INTEGER(mpi) :: ifilb
11075 INTEGER(mpi) :: ioff
11076 INTEGER(mpi) :: iopt
11077 INTEGER(mpi) :: ios
11078 INTEGER(mpi) :: iosum
11081 INTEGER(mpi) :: mat
11082 INTEGER(mpi) :: nab
11083 INTEGER(mpi) :: nline
11084 INTEGER(mpi) :: npat
11085 INTEGER(mpi) :: ntext
11087 INTEGER(mpi) :: nuf
11088 INTEGER(mpi) :: nums
11089 INTEGER(mpi) :: nufile
11090 INTEGER(mpi) :: lenfileInfo
11091 INTEGER(mpi) :: lenFileNames
11092 INTEGER(mpi) :: matint
11093 INTEGER(mpi),
DIMENSION(:,:),
ALLOCATABLE :: vecfileInfo
11094 INTEGER(mpi),
DIMENSION(:,:),
ALLOCATABLE :: tempArray
11095 INTEGER(mpl) :: rows
11096 INTEGER(mpl) :: cols
11097 INTEGER(mpl) :: newcols
11098 INTEGER(mpl) :: length
11100 CHARACTER (LEN=1024) :: text
11101 CHARACTER (LEN=1024) :: fname
11102 CHARACTER (LEN=14) :: bite(3)
11103 CHARACTER (LEN=32) :: keystx
11104 INTEGER(mpi),
PARAMETER :: mnum=100
11105 REAL(mpd) :: dnum(mnum)
11109 SUBROUTINE initc(nfiles)
BIND(c)
11111 INTEGER(c_int),
INTENT(IN),
VALUE :: nfiles
11112 END SUBROUTINE initc
11117 DATA bite/
'C_binary',
'text ',
'Fortran_binary'/
11132 WRITE(*,*)
'Command line options: '
11133 WRITE(*,*)
'--------------------- '
11135 CALL getarg(i,text)
11136 CALL rltext(text,ia,ib,nab)
11137 WRITE(*,101) i,text(1:nab)
11138 IF(text(ia:ia) /=
'-')
THEN
11139 nu=nufile(text(ia:ib))
11142 WRITE(*,*)
'Second text file in command line - stop'
11143 CALL peend(12,
'Aborted, second text file in command line')
11149 WRITE(*,*)
'Open error for file:',text(ia:ib),
' - stop'
11150 CALL peend(16,
'Aborted, open error for file')
11151 IF(text(ia:ia) /=
'/')
THEN
11152 CALL getenv(
'PWD',text)
11153 CALL rltext(text,ia,ib,nab)
11154 WRITE(*,*)
'PWD:',text(ia:ib)
11159 IF(index(text(ia:ib),
'b') /= 0)
THEN
11161 WRITE(*,*)
'Debugging requested'
11163 it=index(text(ia:ib),
't')
11166 ieq=index(text(ia+it:ib),
'=')+it
11167 IF (it /= ieq)
THEN
11168 IF (index(text(ia+ieq:ib),
'SL0' ) /= 0)
ictest=2
11169 IF (index(text(ia+ieq:ib),
'SLE' ) /= 0)
ictest=3
11170 IF (index(text(ia+ieq:ib),
'BP' ) /= 0)
ictest=4
11171 IF (index(text(ia+ieq:ib),
'BRLF') /= 0)
ictest=5
11172 IF (index(text(ia+ieq:ib),
'BRLC') /= 0)
ictest=6
11175 IF(index(text(ia:ib),
's') /= 0)
isubit=1
11176 IF(index(text(ia:ib),
'f') /= 0)
iforce=1
11177 IF(index(text(ia:ib),
'c') /= 0)
icheck=1
11178 IF(index(text(ia:ib),
'C') /= 0)
icheck=2
11180 IF(i == iargc())
WRITE(*,*)
'--------------------- '
11201 CALL rltext(text,ia,ib,nab)
11202 nu=nufile(text(ia:ib))
11206 CALL peend(10,
'Aborted, no steering file')
11207 stop
'in FILETC: no steering file. .'
11220 WRITE(*,*)
'Listing of steering file: ',
filnam(1:
nfnam)
11221 WRITE(*,*)
'-------------------------'
11224 WRITE(*,*)
'Open error for steering file - stop'
11225 CALL peend(11,
'Aborted, open error for steering file')
11226 IF(
filnam(1:1) /=
'/')
THEN
11227 CALL getenv(
'PWD',text)
11228 CALL rltext(text,ia,ib,nab)
11229 WRITE(*,*)
'PWD:',text(ia:ib)
11238 rows=6; cols=lenfileinfo
11239 CALL mpalloc(vecfileinfo,rows,cols,
'file info from steering')
11242 READ(10,102,iostat=ierrf) text
11243 IF (ierrf < 0)
EXIT
11244 CALL rltext(text,ia,ib,nab)
11246 IF(nline <= 50)
THEN
11247 WRITE(*,101) nline,text(1:nab)
11248 IF(nline == 50)
WRITE(*,*)
' ...'
11252 CALL rltext(text,ia,ib,nab)
11253 IF(ib == ia+2)
THEN
11254 mat=matint(text(ia:ib),
'end',npat,ntext)
11255 IF(mat == max(npat,ntext))
THEN
11258 WRITE(*,*)
' end-statement after',nline,
' text lines'
11263 keystx=
'fortranfiles'
11264 mat=matint(text(ia:ib),keystx,npat,ntext)
11265 IF(mat == max(npat,ntext))
THEN
11272 mat=matint(text(ia:ib),keystx,npat,ntext)
11273 IF(mat == max(npat,ntext))
THEN
11279 keystx=
'closeandreopen'
11280 mat=matint(text(ia:ib),keystx,npat,ntext)
11281 IF(mat == max(npat,ntext))
THEN
11289 iopt=index(text(ia:ib),
' -- ')
11290 IF (iopt > 0) ie=iopt-1
11293 nu=nufile(text(ia:ie))
11295 IF (
nfiles == lenfileinfo)
THEN
11296 CALL mpalloc(temparray,rows,cols,
'temp file info from steering')
11297 temparray=vecfileinfo
11299 lenfileinfo=lenfileinfo*2
11300 newcols=lenfileinfo
11301 CALL mpalloc(vecfileinfo,rows,newcols,
'file info from steering')
11302 vecfileinfo(:,1:cols)=temparray(:,1:cols)
11308 lenfilenames=lenfilenames+ie-ia+1
11309 vecfileinfo(1,
nfiles)=nline
11310 vecfileinfo(2,
nfiles)=nu
11311 vecfileinfo(3,
nfiles)=ia
11312 vecfileinfo(4,
nfiles)=ie
11313 vecfileinfo(5,
nfiles)=iopt
11314 vecfileinfo(6,
nfiles)=ib
11324 CALL mpalloc(
nfd,length,
'file line (in steering)')
11327 length=lenfilenames
11333 READ(10,102,iostat=ierrf) text
11334 IF (ierrf < 0)
EXIT
11336 IF (nline == vecfileinfo(1,i))
THEN
11337 nfd(i)=vecfileinfo(1,i)
11338 mfd(i)=vecfileinfo(2,i)
11339 ia=vecfileinfo(3,i)-1
11340 lfd(i)=vecfileinfo(4,i)-ia
11342 tfd(ioff+k)=text(ia+k:ia+k)
11347 IF (vecfileinfo(5,i) > 0)
THEN
11348 CALL ratext(text(vecfileinfo(5,i)+4:vecfileinfo(6,i)),nums,dnum,mnum)
11349 IF (nums > 0)
ofd(i)=real(dnum(1),mps)
11359 CALL mpalloc(
ifd,length,
'integrated record numbers (=offset)')
11360 CALL mpalloc(
jfd,length,
'number of accepted records')
11361 CALL mpalloc(
kfd,rows,length,
'number of records in file, file order')
11366 CALL mpalloc(
sfd,rows,length,
'start, end of file name in TFD')
11367 CALL mpalloc(
yfd,length,
'modification date')
11370 WRITE(*,*)
'-------------------------'
11376 WRITE(*,*)
'Table of files:'
11377 WRITE(*,*)
'---------------'
11380 WRITE(8,*)
'Text and data files:'
11384 fname(k:k)=
tfd(ioff+k)
11387 IF (
mprint > 1)
WRITE(*,103) i,bite(
mfd(i)),fname(1:
lfd(i))
11388 WRITE(8,103) i,bite(
mfd(i)),fname(1:
lfd(i))
11392 WRITE(*,*)
'---------------'
11406 IF(
mfd(i) == 3)
THEN
11430 IF(
mfd(i) == 1)
THEN
11451 WRITE(*,*)
'Opening of C-files not supported.'
11467 IF(iosum /= 0)
THEN
11468 CALL peend(15,
'Aborted, open error(s) for binary files')
11469 stop
'FILETC: open error '
11471 IF(
nfilb == 0)
THEN
11472 CALL peend(14,
'Aborted, no binary files')
11473 stop
'FILETC: no binary files '
11476 WRITE(*,*)
nfilb,
' binary files opened'
11478 WRITE(*,*)
nfilb,
' binary files opened and closed'
11482103
FORMAT(i3,2x,a14,3x,a)
11545 INTEGER(mpi) :: ierrf
11546 INTEGER(mpi) :: ioff
11547 INTEGER(mpi) :: ios
11548 INTEGER(mpi) :: iosum
11550 INTEGER(mpi) :: mat
11551 INTEGER(mpi) :: nab
11552 INTEGER(mpi) :: nfiln
11553 INTEGER(mpi) :: nline
11554 INTEGER(mpi) :: nlinmx
11555 INTEGER(mpi) :: npat
11556 INTEGER(mpi) :: ntext
11557 INTEGER(mpi) :: matint
11561 CHARACTER (LEN=1024) :: text
11562 CHARACTER (LEN=1024) :: fname
11565 WRITE(*,*)
'Processing text files ...'
11578 IF(
mfd(i) /= 2) cycle
11580 fname(k:k)=
tfd(ia+k)
11582 WRITE(*,*)
'File ',fname(1:
lfd(i))
11583 IF (
mprint > 1)
WRITE(*,*)
' '
11584 OPEN(10,file=fname(1:
lfd(i)),iostat=ios,form=
'FORMATTED')
11586 WRITE(*,*)
'Open error for file ',fname(1:
lfd(i))
11596 READ(10,102,iostat=ierrf) text
11597 IF (ierrf < 0)
THEN
11600 WRITE(*,*)
' end-of-file after',nline,
' text lines'
11604 IF(nline <= nlinmx.AND.
mprint > 1)
THEN
11605 CALL rltext(text,ia,ib,nab)
11607 WRITE(*,101) nline,text(1:nab)
11608 IF(nline == nlinmx)
WRITE(*,*)
' ...'
11611 CALL rltext(text,ia,ib,nab)
11612 IF(ib == ia+2)
THEN
11613 mat=matint(text(ia:ib),
'end',npat,ntext)
11614 IF(mat == max(npat,ntext))
THEN
11617 WRITE(*,*)
' end-statement after',nline,
' text lines'
11623 IF(nfiln <=
nfiles)
THEN
11624 IF(nline ==
nfd(nfiln))
THEN
11639 IF(iosum /= 0)
THEN
11640 CALL peend(16,
'Aborted, open error(s) for text files')
11641 stop
'FILETX: open error(s) in text files '
11644 WRITE(*,*)
'... end of text file processing.'
11649 WRITE(*,*)
lunkno,
' unknown keywords in steering files, ', &
11650 'or file non-existing,'
11651 WRITE(*,*)
' see above!'
11652 WRITE(*,*)
'------------> stop'
11654 CALL peend(13,
'Aborted, unknown keywords in steering file')
11663 ELSE IF(
matsto == 1)
THEN
11665 ELSE IF(
matsto == 2)
THEN
11668 ELSE IF(
metsol == 1)
THEN
11670 ELSE IF(
metsol == 2)
THEN
11672 ELSE IF(
metsol == 3)
THEN
11674 ELSE IF(
metsol == 4)
THEN
11676 ELSE IF(
metsol == 5)
THEN
11678 ELSE IF(
metsol == 6)
THEN
11681 ELSE IF(
metsol == 7)
THEN
11683 ELSE IF(
metsol == 8)
THEN
11686 ELSE IF(
metsol == 9)
THEN
11691 WRITE(*,*)
'MINRES forced with sparse matrix!'
11693 WRITE(*,*)
'MINRES forced with sparse matrix!'
11695 WRITE(*,*)
'MINRES forced with sparse matrix!'
11700 WRITE(*,*)
'MINRES forced with sparse matrix!'
11702 WRITE(*,*)
'MINRES forced with sparse matrix!'
11704 WRITE(*,*)
'MINRES forced with sparse matrix!'
11712 WRITE(*,*)
'Solution method and matrix-storage mode:'
11714 WRITE(*,*)
' METSOL = 1: matrix inversion'
11715 ELSE IF(
metsol == 2)
THEN
11716 WRITE(*,*)
' METSOL = 2: diagonalization'
11717 ELSE IF(
metsol == 3)
THEN
11718 WRITE(*,*)
' METSOL = 3: decomposition'
11719 ELSE IF(
metsol == 4)
THEN
11720 WRITE(*,*)
' METSOL = 4: MINRES'
11721 ELSE IF(
metsol == 5)
THEN
11722 WRITE(*,*)
' METSOL = 5: MINRES-QLP'
11723 ELSE IF(
metsol == 6)
THEN
11724 WRITE(*,*)
' METSOL = 6: GMRES (-> MINRES)'
11726 ELSE IF(
metsol == 7)
THEN
11727 WRITE(*,*)
' METSOL = 7: LAPACK factorization'
11728 ELSE IF(
metsol == 8)
THEN
11729 WRITE(*,*)
' METSOL = 8: LAPACK factorization'
11731 ELSE IF(
metsol == 9)
THEN
11732 WRITE(*,*)
' METSOL = 9: Intel oneMKL PARDISO'
11737 WRITE(*,*)
' with',
mitera,
' iterations'
11740 WRITE(*,*)
' MATSTO = 0: unpacked symmetric matrix, ',
'n*n elements'
11741 ELSEIF(
matsto == 1)
THEN
11742 WRITE(*,*)
' MATSTO = 1: full symmetric matrix, ',
'(n*n+n)/2 elements'
11743 ELSE IF(
matsto == 2)
THEN
11744 WRITE(*,*)
' MATSTO = 2: sparse matrix (custom)'
11745 ELSE IF(
matsto == 3)
THEN
11747 WRITE(*,*)
' MATSTO = 3: sparse matrix (upper triangle, CSR3)'
11749 WRITE(*,*)
' MATSTO = 3: sparse matrix (upper triangle, BSR3)'
11753 WRITE(*,*)
' and band matrix, width',
mbandw
11757 WRITE(*,*)
'Chi square cut equiv 3 st.dev applied ...'
11758 WRITE(*,*)
' in first iteration with factor',
chicut
11759 WRITE(*,*)
' in second iteration with factor',
chirem
11760 WRITE(*,*)
' (reduced by sqrt in next iterations)'
11764 WRITE(*,*)
' Down-weighting of outliers in',
lhuber,
' iterations'
11765 WRITE(*,*)
' Cut on downweight fraction',
dwcut
11768 WRITE(*,*)
'Iterations (solutions) with line search:'
11777 WRITE(*,*)
' All with Chi square cut scaling factor <= 1.'
11808 INTEGER(mpi) :: ios
11812 INTEGER(mpi) :: npat
11813 INTEGER(mpi) :: ntext
11814 INTEGER(mpi) :: nuprae
11817 CHARACTER (LEN=*),
INTENT(INOUT) :: fname
11823 IF(len(fname) > 5)
THEN
11824 IF(fname(1:5) ==
'rfio:') nuprae=1
11825 IF(fname(1:5) ==
'dcap:') nuprae=2
11826 IF(fname(1:5) ==
'root:') nuprae=3
11828 IF(nuprae == 0)
THEN
11829 INQUIRE(file=fname,iostat=ios,exist=ex)
11830 IF(ios /= 0)
nufile=-abs(ios)
11831 IF(ios /= 0)
RETURN
11832 ELSE IF(nuprae == 1)
THEN
11845 nm=
matint(
'xt',fname(l1:ll),npat,ntext)
11848 nm=
matint(
'tx',fname(l1:ll),npat,ntext)
11869 INTEGER(mpi) :: ier
11870 INTEGER(mpi) :: iomp
11873 INTEGER(mpi) :: kkey
11874 INTEGER(mpi) :: label
11875 INTEGER(mpi) :: lkey
11876 INTEGER(mpi) :: mat
11877 INTEGER(mpi) :: miter
11878 INTEGER(mpi) :: nab
11879 INTEGER(mpi) :: nkey
11880 INTEGER(mpi) :: nkeys
11882 INTEGER(mpi) :: nmeth
11883 INTEGER(mpi) :: npat
11884 INTEGER(mpi) :: ntext
11885 INTEGER(mpi) :: nums
11886 INTEGER(mpi) :: matint
11888 CHARACTER (LEN=*),
INTENT(IN) :: text
11889 INTEGER(mpi),
INTENT(IN) :: nline
11893 parameter(nkeys=7,nmeth=10)
11895 parameter(nkeys=6,nmeth=9)
11898 parameter(nkeys=6,nmeth=7)
11900 CHARACTER (LEN=16) :: methxt(nmeth)
11901 CHARACTER (LEN=16) :: keylst(nkeys)
11902 CHARACTER (LEN=32) :: keywrd
11903 CHARACTER (LEN=32) :: keystx
11904 CHARACTER (LEN=itemCLen) :: ctext
11905 INTEGER(mpi),
PARAMETER :: mnum=100
11906 REAL(mpd) :: dnum(mnum)
11909 INTEGER(mpi) :: ipvs
11912 INTEGER(mpi) :: lpvs
11916 SUBROUTINE additem(length,list,label,value)
11918 INTEGER(mpi),
INTENT(IN OUT) :: length
11919 TYPE(listitem),
DIMENSION(:),
INTENT(IN OUT),
ALLOCATABLE :: list
11920 INTEGER(mpi),
INTENT(IN) :: label
11921 REAL(mpd),
INTENT(IN) :: value
11923 SUBROUTINE additemc(length,list,label,text)
11925 INTEGER(mpi),
INTENT(IN OUT) :: length
11926 TYPE(listitemc),
DIMENSION(:),
INTENT(IN OUT),
ALLOCATABLE :: list
11927 INTEGER(mpi),
INTENT(IN) :: label
11928 CHARACTER(LEN = itemCLen),
INTENT(IN) :: text
11930 SUBROUTINE additemi(length,list,label,ivalue)
11932 INTEGER(mpi),
INTENT(IN OUT) :: length
11933 TYPE(listitemi),
DIMENSION(:),
INTENT(IN OUT),
ALLOCATABLE :: list
11934 INTEGER(mpi),
INTENT(IN) :: label
11935 INTEGER(mpi),
INTENT(IN) :: ivalue
11942 DATA keylst/
'unknown',
'parameter',
'constraint',
'measurement',
'method',
'comment',
'pardiso'/
11943 DATA methxt/
'diagonalization',
'inversion',
'fullMINRES',
'sparseMINRES', &
11944 'fullMINRES-QLP',
'sparseMINRES-QLP',
'decomposition',
'fullLAPACK',
'unpackedLAPACK', &
11947 DATA keylst/
'unknown',
'parameter',
'constraint',
'measurement',
'method',
'comment'/
11948 DATA methxt/
'diagonalization',
'inversion',
'fullMINRES',
'sparseMINRES', &
11949 'fullMINRES-QLP',
'sparseMINRES-QLP',
'decomposition',
'fullLAPACK',
'unpackedLAPACK'/
11952 DATA keylst/
'unknown',
'parameter',
'constraint',
'measurement',
'method',
'comment'/
11953 DATA methxt/
'diagonalization',
'inversion',
'fullMINRES',
'sparseMINRES', &
11954 'fullMINRES-QLP',
'sparseMINRES-QLP',
'decomposition'/
11960 CALL rltext(text,ia,ib,nab)
11961 IF(nab == 0)
GOTO 10
11962 CALL ratext(text(1:nab),nums,dnum,mnum)
11964 IF(nums /= 0) nkey=0
11972 keystx=keylst(nkey)
11973 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
11974 IF(100*mat >= 80*max(npat,ntext))
GO TO 10
11980 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
11981 IF(100*mat >= 80*max(npat,ntext))
THEN
11983 IF(nums > 0) mprint=nint(dnum(1),mpi)
11988 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
11989 IF(100*mat >= 80*max(npat,ntext))
THEN
11992 IF(nums > 0) mdebug=nint(dnum(1),mpi)
11993 IF(nums > 1) mdebg2=nint(dnum(2),mpi)
11998 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
11999 IF(100*mat >= 80*max(npat,ntext))
THEN
12000 IF(nums > 0 .AND. dnum(1) > 0.5) mreqenf=nint(dnum(1),mpi)
12001 IF(nums > 1 .AND. dnum(2) > 0.5) mreqena=nint(dnum(2),mpi)
12002 IF(nums > 2 .AND. dnum(3) > 0.5) iteren=nint(dnum(1)*dnum(3),mpi)
12006 keystx=
'printrecord'
12007 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12008 IF(100*mat >= 80*max(npat,ntext))
THEN
12009 IF(nums > 0) nrecpr=nint(dnum(1),mpi)
12010 IF(nums > 1) nrecp2=nint(dnum(2),mpi)
12015 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12016 IF(100*mat >= 80*max(npat,ntext))
THEN
12017 IF (nums > 0.AND.dnum(1) > 0.) mxrec=nint(dnum(1),mpi)
12022 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12023 IF(100*mat >= 80*max(npat,ntext))
THEN
12024 IF (nums > 0.AND.dnum(1) >= 0.) ncache=nint(dnum(1),mpi)
12025 IF (nums == 2.AND.dnum(2) > 0..AND.dnum(2) <= 1.0) &
12026 fcache(1)=real(dnum(2),mps)
12027 IF (nums >= 4)
THEN
12029 fcache(k)=real(dnum(k+1),mps)
12036 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12037 IF(100*mat >= 80*max(npat,ntext))
THEN
12042 chicut=real(dnum(1),mps)
12043 IF(chicut < 1.0) chicut=-1.0
12047 chirem=real(dnum(2),mps)
12048 IF(chirem < 1.0) chirem=1.0
12049 IF(chicut >= 1.0) chirem=min(chirem,chicut)
12057 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12058 IF(100*mat >= 80*max(npat,ntext))
THEN
12059 IF(nums > 0) chhuge=real(dnum(1),mps)
12060 IF(chhuge < 1.0) chhuge=1.0
12065 keystx=
'linesearch'
12066 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12067 IF(100*mat >= 80*max(npat,ntext))
THEN
12068 IF(nums > 0) lsearch=nint(dnum(1),mpi)
12073 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12074 IF(100*mat >= 80*max(npat,ntext))
THEN
12075 IF(nums > 0) lfitnp=nint(dnum(1),mpi)
12076 IF(nums > 1) lfitbb=nint(dnum(2),mpi)
12080 keystx=
'regularization'
12081 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12082 IF(100*mat >= 80*max(npat,ntext))
THEN
12084 regula=real(dnum(1),mps)
12085 IF(nums >= 2) regpre=real(dnum(2),mps)
12089 keystx=
'regularisation'
12090 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12091 IF(100*mat >= 80*max(npat,ntext))
THEN
12093 regula=real(dnum(1),mps)
12094 IF(nums >= 2) regpre=real(dnum(2),mps)
12099 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12100 IF(100*mat >= 80*max(npat,ntext))
THEN
12101 regpre=real(dnum(1),mps)
12106 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12107 IF(100*mat >= 80*max(npat,ntext))
THEN
12108 matrit=nint(dnum(1),mpi)
12113 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12114 IF(100*mat >= 80*max(npat,ntext))
THEN
12116 IF (nums > 0.AND.dnum(1) > 0.) matmon=nint(dnum(1),mpi)
12121 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12122 IF(100*mat >= 80*max(npat,ntext))
THEN
12123 IF(nums > 0) mbandw=nint(dnum(1),mpi)
12124 IF(mbandw < 0) mbandw=-1
12125 IF(nums > 1) lprecm=nint(dnum(2),mpi)
12149 keystx=
'outlierdownweighting'
12150 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12151 IF(100*mat >= 80*max(npat,ntext))
THEN
12152 lhuber=nint(dnum(1),mpi)
12153 IF(lhuber > 0.AND.lhuber <= 2) lhuber=2
12157 keystx=
'dwfractioncut'
12158 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12159 IF(100*mat >= 80*max(npat,ntext))
THEN
12160 dwcut=real(dnum(1),mps)
12161 IF(dwcut > 0.5) dwcut=0.5
12166 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12167 IF(100*mat >= 80*max(npat,ntext))
THEN
12168 prange=abs(real(dnum(1),mps))
12173 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12174 IF(100*mat >= 80*max(npat,ntext))
THEN
12180 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12181 IF(100*mat >= 80*max(npat,ntext))
THEN
12186 keystx=
'memorydebug'
12187 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12188 IF(100*mat >= 80*max(npat,ntext))
THEN
12190 IF (nums > 0.AND.dnum(1) > 0.0) memdbg=nint(dnum(1),mpi)
12194 keystx=
'globalcorr'
12195 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12196 IF(100*mat >= 80*max(npat,ntext))
THEN
12201 keystx=
'printcounts'
12202 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12203 IF(100*mat >= 80*max(npat,ntext))
THEN
12205 IF (nums > 0.AND.dnum(1) > 0.0) ipcntr=nint(dnum(1),mpi)
12209 keystx=
'weightedcons'
12210 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12211 IF(100*mat >= 80*max(npat,ntext))
THEN
12213 IF (nums > 0) iwcons=nint(dnum(1),mpi)
12217 keystx=
'skipemptycons'
12218 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12219 IF(100*mat >= 80*max(npat,ntext))
THEN
12224 keystx=
'resolveredundancycons'
12225 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12226 IF(100*mat >= 80*max(npat,ntext))
THEN
12231 keystx=
'withelimination'
12232 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12233 IF(100*mat >= 80*max(npat,ntext))
THEN
12239 keystx=
'withLAPACKelimination'
12240 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12241 IF(100*mat >= 80*max(npat,ntext))
THEN
12247 keystx=
'withmultipliers'
12248 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12249 IF(100*mat >= 80*max(npat,ntext))
THEN
12254 keystx=
'checkinput'
12255 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12256 IF(100*mat >= 80*max(npat,ntext))
THEN
12258 IF (nums > 0) icheck=nint(dnum(1),mpi)
12262 keystx=
'checkparametergroups'
12263 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12264 IF(100*mat >= 80*max(npat,ntext))
THEN
12269 keystx=
'monitorresiduals'
12270 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12271 IF(100*mat >= 80*max(npat,ntext))
THEN
12273 IF (nums > 0) imonit=nint(dnum(1),mpi)
12274 IF (nums > 1) measbins=max(measbins,nint(dnum(2),mpi))
12278 keystx=
'monitorpulls'
12279 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12280 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=
'monitorprogress'
12289 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12290 IF(100*mat >= 80*max(npat,ntext))
THEN
12293 IF (nums > 0) monpg1=max(1,nint(dnum(1),mpi))
12294 IF (nums > 1) monpg2=max(1,nint(dnum(2),mpi))
12298 keystx=
'scaleerrors'
12299 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12300 IF(100*mat >= 80*max(npat,ntext))
THEN
12302 IF (nums > 0) dscerr(1:2)=dnum(1)
12303 IF (nums > 1) dscerr(2)=dnum(2)
12307 keystx=
'iterateentries'
12308 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12309 IF(100*mat >= 80*max(npat,ntext))
THEN
12310 iteren=huge(iteren)
12311 IF (nums > 0) iteren=nint(dnum(1),mpi)
12316 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12317 IF(100*mat >= 80*max(npat,ntext))
THEN
12325 WRITE(*,*)
'WARNING: multithreading not available'
12331 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12332 IF(100*mat >= 80*max(npat,ntext))
THEN
12333 WRITE(*,*)
'WARNING: keyword COMPRESS is obsolete (compression is default)'
12345 keystx=
'countrecords'
12346 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12347 IF(100*mat >= 80*max(npat,ntext))
THEN
12353 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12354 IF(100*mat >= 80*max(npat,ntext).AND.mnrsel < 100)
THEN
12355 nl=min(nums,100-mnrsel)
12357 lbmnrs(mnrsel+k)=nint(dnum(k),mpi)
12363 keystx=
'pairentries'
12364 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12365 IF(100*mat >= 80*max(npat,ntext))
THEN
12368 IF (nums > 0.AND.dnum(1) > 0.0)
THEN
12369 mreqpe=nint(dnum(1),mpi)
12370 IF (nums >= 2.AND.dnum(2) >= dnum(1)) mhispe=nint(dnum(2),mpi)
12371 IF (nums >= 3.AND.dnum(3) >= dnum(1)) msngpe=nint(dnum(3),mpi)
12377 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12378 IF(100*mat >= 80*max(npat,ntext))
THEN
12379 wolfc1=real(dnum(1),mps)
12380 wolfc2=real(dnum(2),mps)
12387 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12388 IF(100*mat >= 80*max(npat,ntext))
THEN
12390 IF (dnum(1) < 1.0e-10_mpd.OR.dnum(1) > 1.0e-04_mpd)
THEN
12391 WRITE(*,*)
'ERROR: need 1.0D-10 <= MRESTL ', &
12392 '<= 1.0D-04, but get ', dnum(1)
12401 keystx=
'mrestranscond'
12402 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12403 IF(100*mat >= 80*max(npat,ntext))
THEN
12411 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12412 IF(100*mat >= 80*max(npat,ntext))
THEN
12414 mrmode = int(dnum(1),mpi)
12419 keystx=
'nofeasiblestart'
12420 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12421 IF(100*mat >= 80*max(npat,ntext))
THEN
12427 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12428 IF(100*mat >= 80*max(npat,ntext))
THEN
12433 keystx=
'readerroraseof'
12434 mat=matint(text(ia:ib),keystx,npat,ntext)
12435 IF(100*mat >= 80*max(npat,ntext))
THEN
12441 keystx=
'LAPACKwitherrors'
12442 mat=matint(text(ia:ib),keystx,npat,ntext)
12443 IF(100*mat >= 80*max(npat,ntext))
THEN
12448 keystx=
'debugPARDISO'
12449 mat=matint(text(ia:ib),keystx,npat,ntext)
12450 IF(100*mat >= 80*max(npat,ntext))
THEN
12455 keystx=
'blocksizePARDISO'
12456 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12457 IF(100*mat >= 80*max(npat,ntext).AND.mnrsel < 100)
THEN
12458 nl=min(nums,10-mpdbsz)
12460 IF (nint(dnum(k),mpi) > 0)
THEN
12461 IF (mpdbsz == 0)
THEN
12463 ipdbsz(mpdbsz)=nint(dnum(k),mpi)
12464 ELSE IF (nint(dnum(k),mpi) > ipdbsz(mpdbsz))
THEN
12466 ipdbsz(mpdbsz)=nint(dnum(k),mpi)
12474 keystx=
'fortranfiles'
12475 mat=matint(text(ia:ib),keystx,npat,ntext)
12476 IF(mat == max(npat,ntext))
RETURN
12479 mat=matint(text(ia:ib),keystx,npat,ntext)
12480 IF(mat == max(npat,ntext))
RETURN
12482 keystx=
'closeandreopen'
12483 mat=matint(text(ia:ib),keystx,npat,ntext)
12484 IF(mat == max(npat,ntext))
RETURN
12488 IF(nums /= 0) nkey=0
12491 WRITE(*,*)
'**************************************************'
12493 WRITE(*,*)
'Unknown keyword(s): ',text(1:min(nab,50))
12495 WRITE(*,*)
'**************************************************'
1251710
IF(nkey > 0)
THEN
12521 lpvs=nint(dnum(1),mpi)
12523 CALL additem(lenparameters,listparameters,lpvs,dnum(2))
12524 CALL additem(lenpresigmas,listpresigmas,lpvs,dnum(3))
12526 WRITE(*,*)
'Line',nline,
' error, label=',lpvs
12528 ELSE IF(nums /= 0)
THEN
12530 WRITE(*,*)
'Wrong text in line',nline
12531 WRITE(*,*)
'Status: new parameter'
12532 WRITE(*,*)
'> ',text(1:nab)
12534 ELSE IF(lkey == 3)
THEN
12536 IF(nums >= 1.AND.nums <= 2)
THEN
12538 CALL additem(lenconstraints,listconstraints,lpvs,dnum(1))
12540 IF(iwcons > 0) lpvs=-2
12542 IF(nums == 2) plvs=dnum(2)
12543 CALL additem(lenconstraints,listconstraints,lpvs,plvs)
12546 WRITE(*,*)
'Wrong text in line',nline
12547 WRITE(*,*)
'Status: new keyword constraint'
12548 WRITE(*,*)
'> ',text(1:nab)
12550 ELSE IF(lkey == 4)
THEN
12552 nummeasurements=nummeasurements+1
12554 CALL additem(lenmeasurements,listmeasurements,lpvs,dnum(1))
12556 CALL additem(lenmeasurements,listmeasurements,lpvs,dnum(2))
12559 WRITE(*,*)
'Wrong text in line',nline
12560 WRITE(*,*)
'Status: new keyword measurement'
12561 WRITE(*,*)
'> ',text(1:nab)
12563 ELSE IF(lkey == 5.AND.
keyb <
keyc)
THEN
12565 IF(nums >= 1) miter=nint(dnum(1),mpi)
12566 IF(miter >= 1) mitera=miter
12567 dflim=real(dnum(2),mps)
12571 mat=matint(text(
keyb+1:
keyc),keystx,npat,ntext)
12572 IF(100*mat >= 80*max(npat,ntext))
THEN
12576 ELSE IF(i == 2)
THEN
12579 ELSE IF(i == 3)
THEN
12582 ELSE IF(i == 4)
THEN
12585 ELSE IF(i == 5)
THEN
12588 ELSE IF(i == 6)
THEN
12591 ELSE IF(i == 7)
THEN
12595 ELSE IF(i == 8)
THEN
12598 ELSE IF(i == 9)
THEN
12602 ELSE IF(i == 10)
THEN
12611 ELSE IF(nkey == 0)
THEN
12614 lpvs=nint(dnum(1),mpi)
12616 CALL additem(lenparameters,listparameters,lpvs,dnum(2))
12617 CALL additem(lenpresigmas,listpresigmas,lpvs,dnum(3))
12619 WRITE(*,*)
'Line',nline,
' error, label=',lpvs
12621 ELSE IF(nums > 1.AND.nums < 3)
THEN
12623 WRITE(*,*)
'Wrong text in line',nline
12624 WRITE(*,*)
'Status continuation parameter'
12625 WRITE(*,*)
'> ',text(1:nab)
12628 ELSE IF(lkey == 3)
THEN
12631 label=nint(dnum(i),mpi)
12632 IF(label <= 0) ier=1
12634 IF(mod(nums,2) /= 0) ier=1
12637 lpvs=nint(dnum(i),mpi)
12639 CALL additem(lenconstraints,listconstraints,lpvs,plvs)
12643 WRITE(*,*)
'Wrong text in line',nline
12644 WRITE(*,*)
'Status continuation constraint'
12645 WRITE(*,*)
'> ',text(1:nab)
12648 ELSE IF(lkey == 4)
THEN
12652 label=nint(dnum(i),mpi)
12653 IF(label <= 0) ier=1
12655 IF(mod(nums,2) /= 0) ier=1
12659 lpvs=nint(dnum(i),mpi)
12661 CALL additem(lenmeasurements,listmeasurements,lpvs,plvs)
12665 WRITE(*,*)
'Wrong text in line',nline
12666 WRITE(*,*)
'Status continuation measurement'
12667 WRITE(*,*)
'> ',text(1:nab)
12669 ELSE IF(lkey == 6)
THEN
12671 lpvs=nint(dnum(1),mpi)
12675 IF (text(j:j) ==
' ')
EXIT
12678 CALL additemc(lencomments,listcomments,lpvs,ctext)
12680 WRITE(*,*)
'Line',nline,
' error, label=',lpvs
12682 ELSE IF(nums /= 0)
THEN
12684 WRITE(*,*)
'Wrong text in line',nline
12685 WRITE(*,*)
'Status: continuation comment'
12686 WRITE(*,*)
'> ',text(1:nab)
12690 ELSE IF(lkey == 7)
THEN
12693 label=nint(dnum(i),mpi)
12694 IF(label <= 0.OR.label > 64) ier=1
12696 IF(mod(nums,2) /= 0) ier=1
12700 lpvs=nint(dnum(i),mpi)
12701 ipvs=nint(dnum(i+1),mpi)
12702 CALL additemi(lenpardiso,listpardiso,lpvs,ipvs)
12706 WRITE(*,*)
'Wrong text in line',nline
12707 WRITE(*,*)
'Status continuation measurement'
12708 WRITE(*,*)
'> ',text(1:nab)
12727 INTEGER(mpi),
INTENT(IN OUT) :: length
12728 TYPE(
listitem),
DIMENSION(:),
INTENT(IN OUT),
ALLOCATABLE :: list
12729 INTEGER(mpi),
INTENT(IN) :: label
12730 REAL(mpd),
INTENT(IN) :: value
12732 INTEGER(mpl) :: newSize
12733 INTEGER(mpl) :: oldSize
12734 TYPE(
listitem),
DIMENSION(:),
ALLOCATABLE :: tempList
12736 IF (label > 0.AND.
value == 0.)
RETURN
12737 IF (length == 0 )
THEN
12739 CALL mpalloc(list,newsize,
' list ')
12741 oldsize=
size(list,kind=
mpl)
12742 IF (length >= oldsize)
THEN
12743 newsize = oldsize + oldsize/5 + 100
12744 CALL mpalloc(templist,oldsize,
' temp. list ')
12747 CALL mpalloc(list,newsize,
' list ')
12748 list(1:oldsize)=templist(1:oldsize)
12753 list(length)%label=label
12754 list(length)%value=
value
12769 INTEGER(mpi),
INTENT(IN OUT) :: length
12770 TYPE(
listitemc),
DIMENSION(:),
INTENT(IN OUT),
ALLOCATABLE :: list
12771 INTEGER(mpi),
INTENT(IN) :: label
12772 CHARACTER(len = itemCLen),
INTENT(IN) :: text
12774 INTEGER(mpl) :: newSize
12775 INTEGER(mpl) :: oldSize
12776 TYPE(
listitemc),
DIMENSION(:),
ALLOCATABLE :: tempList
12778 IF (label > 0.AND.text ==
'')
RETURN
12779 IF (length == 0 )
THEN
12781 CALL mpalloc(list,newsize,
' list ')
12783 oldsize=
size(list,kind=
mpl)
12784 IF (length >= oldsize)
THEN
12785 newsize = oldsize + oldsize/5 + 100
12786 CALL mpalloc(templist,oldsize,
' temp. list ')
12789 CALL mpalloc(list,newsize,
' list ')
12790 list(1:oldsize)=templist(1:oldsize)
12795 list(length)%label=label
12796 list(length)%text=text
12811 INTEGER(mpi),
INTENT(IN OUT) :: length
12812 TYPE(
listitemi),
DIMENSION(:),
INTENT(IN OUT),
ALLOCATABLE :: list
12813 INTEGER(mpi),
INTENT(IN) :: label
12814 INTEGER(mpi),
INTENT(IN) :: ivalue
12816 INTEGER(mpl) :: newSize
12817 INTEGER(mpl) :: oldSize
12818 TYPE(
listitemi),
DIMENSION(:),
ALLOCATABLE :: tempList
12820 IF (length == 0 )
THEN
12822 CALL mpalloc(list,newsize,
' list ')
12824 oldsize=
size(list,kind=
mpl)
12825 IF (length >= oldsize)
THEN
12826 newsize = oldsize + oldsize/5 + 100
12827 CALL mpalloc(templist,oldsize,
' temp. list ')
12830 CALL mpalloc(list,newsize,
' list ')
12831 list(1:oldsize)=templist(1:oldsize)
12836 list(length)%label=label
12837 list(length)%ivalue=ivalue
12851 CHARACTER (LEN=*),
INTENT(IN) :: text
12852 CHARACTER (LEN=16) :: textc
12861 textl(ka:kb)=text(1:l)
12865 textc=text(1:l)//
'-end'
12873 textl(ka:kb)=textc(1:l)
12900 INTEGER(mpi),
INTENT(IN) :: lun
12901 CHARACTER (LEN=*),
INTENT(IN) :: fname
12902 CHARACTER (LEN=33) :: nafile
12903 CHARACTER (LEN=33) :: nbfile
12909 CALL peend(17,
'Aborted, file name too long')
12910 stop
'File name too long '
12913 nafile(l+1:l+1)=
'~'
12915 INQUIRE(file=nafile(1:l),exist=ex)
12917 INQUIRE(file=nafile(1:l+1),exist=ex)
12919 CALL system(
'rm '//nafile)
12922 nafile(l+1:l+1)=
' '
12923 CALL system(
'mv '//nafile//nbfile)
12925 OPEN(unit=lun,file=fname)
12936 REAL,
DIMENSION(2) :: ta
12956 IF(ncount > 1)
THEN
12958 nsecd1=int(delta,
mpi)
12960 minut1=nsecd1/60-60*nhour1
12961 secnd1=delta-60*(minut1+60*nhour1)
12963 nsecd2=int(delta,
mpi)
12965 minut2=nsecd2/60-60*nhour2
12966 secnd2=delta-60*(minut2+60*nhour2)
12967 WRITE(*,101) nhour1,minut1,secnd1, nhour2,minut2,secnd2
12972101
FORMAT(i4,
' h',i3,
' min',f5.1,
' sec total',18x,
'elapsed', &
12973 i4,
' h',i3,
' min',f5.1,
' sec')
12987 INTEGER(mpi),
INTENT(IN) :: icode
12988 CHARACTER (LEN=*),
INTENT(IN) :: cmessage
12990 CALL mvopen(9,
'millepede.end')
12991 WRITE(9,101) icode, cmessage
12992101
FORMAT(1x,i4,3x,a)
12996END SUBROUTINE peend
13008 INTEGER(mpi),
INTENT(IN) :: kfile
13009 INTEGER(mpi),
INTENT(IN) :: ithr
13010 INTEGER(mpi),
INTENT(OUT) :: ierr
13012 INTEGER(mpi),
DIMENSION(13) :: ibuff
13013 INTEGER(mpi) :: ioff
13014 INTEGER(mpi) :: ios
13016 INTEGER(mpi) :: lfn
13017 INTEGER(mpi) :: lun
13018 INTEGER(mpi) :: moddate
13019 CHARACTER (LEN=1024) :: fname
13020 CHARACTER (LEN=7) :: cfile
13025 SUBROUTINE openc(filename, lfn, lun, ios)
BIND(c)
13027 CHARACTER(kind=c_char),
DIMENSION(*),
INTENT(IN) :: filename
13028 INTEGER(c_int),
INTENT(IN),
VALUE :: lfn
13029 INTEGER(c_int),
INTENT(IN),
VALUE :: lun
13030 INTEGER(c_int),
INTENT(INOUT) :: ios
13031 END SUBROUTINE openc
13043 fname(k:k)=
tfd(ioff+k)
13048 IF(kfile <=
nfilf)
THEN
13051 OPEN(lun,file=fname(1:lfn),iostat=ios, form=
'UNFORMATTED')
13052 print *,
' lun ', lun, ios
13056 CALL openc(fname(1:lfn),lfn,lun,ios)
13058 WRITE(*,*)
'Opening of C-files not supported.'
13065 WRITE(*,*)
'Open error for file ',fname(1:lfn), ios
13066 IF (moddate /= 0)
THEN
13067 WRITE(cfile,
'(I7)') kfile
13068 CALL peend(15,
'Aborted, open error(s) for binary file ' // cfile)
13069 stop
'PEREAD: open error'
13074 ios=stat(fname(1:lfn),ibuff)
13078 WRITE(*,*)
'STAT error for file ',fname(1:lfn), ios
13082 IF (moddate /= 0)
THEN
13083 IF (ibuff(10) /= moddate)
THEN
13084 WRITE(cfile,
'(I7)') kfile
13085 CALL peend(19,
'Aborted, binary file modified (date) ' // cfile)
13086 stop
'PEREAD: file modified'
13089 yfd(kfile)=ibuff(10)
13104 INTEGER(mpi),
INTENT(IN) :: kfile
13105 INTEGER(mpi),
INTENT(IN) :: ithr
13107 INTEGER(mpi) :: lun
13111 SUBROUTINE closec(lun)
BIND(c)
13113 INTEGER(c_int),
INTENT(IN),
VALUE :: lun
13120 IF(kfile <=
nfilf)
THEN
13139 INTEGER(mpi),
INTENT(IN) :: kfile
13141 INTEGER(mpi) :: lun
13145 SUBROUTINE resetc(lun)
BIND(c)
13147 INTEGER(c_int),
INTENT(IN),
VALUE :: lun
13153 IF (kfile <=
nfilf)
THEN
13172 INTEGER(mpi) :: ipgrp
13173 INTEGER(mpi) :: irank
13174 INTEGER(mpi) :: isize
13175 INTEGER(mpi) :: ivoff
13176 INTEGER(mpi) :: itgbi
13178 INTEGER(mpi) :: msize
13179 INTEGER(mpi),
PARAMETER :: mxsize = 1000
13181 INTEGER(mpl):: length
13183 REAL(mpd),
DIMENSION(:),
ALLOCATABLE :: auxVectorD
13184 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: auxVectorI
13185 REAL(mpd),
DIMENSION(:),
ALLOCATABLE :: resParGroup
13186 REAL(mpd),
DIMENSION(:),
ALLOCATABLE :: blockParGroup
13194 IF (isize <= mxsize)
THEN
13195 msize=max(msize,isize)
13197 print *,
' CKPGRP: par. group', ipgrp,
' not checked -- too large: ', isize
13200 IF (msize == 0)
RETURN
13203 length=int(msize,mpl)*(int(msize,mpl)+1)/2
13204 CALL mpalloc(blockpargroup,length,
'(matrix) block for parameter groups (D)')
13206 CALL mpalloc(respargroup,length,
'residuals for parameter groups (D)')
13207 CALL mpalloc(auxvectori,length,
'auxiliary array (I)')
13208 CALL mpalloc(auxvectord,length,
'auxiliary array (D)')
13212 print *,
' CKPGRP par. group first label size rank'
13215 IF (isize > mxsize) cycle
13222 blockpargroup(ij)=matij(ivoff+i,ivoff+j)
13226 CALL sqminv(blockpargroup,respargroup,isize,irank, auxvectord, auxvectori)
13229 IF (isize == irank)
THEN
13233 print *,
' CKPGRP ', ipgrp,
globalparlabelindex(1,itgbi), isize, irank,
' rank deficit !!!'
13251 INTEGER(mpl) :: nan
13252 INTEGER(mpl) :: neg
13254 print *,
' Checking global matrix(D) for NANs ',
size(
globalmatd,kind=mpl)
13259 print *,
' i, nan ', i, nan
13265 print *,
' Checking diagonal elements ',
nagb
13270 print *,
' i, neg ', i, neg
13274 print *,
' CHKMAT summary ', nan, neg
13296 REAL(mpd),
INTENT(IN) :: chi2
13297 INTEGER(mpi),
INTENT(IN) :: ithrd
13298 INTEGER(mpi),
INTENT(IN) :: ndf
13299 REAL(mpd),
INTENT(IN) :: dw
13301 INTEGER(mpl) ::nadd
13329 REAL(mpd),
INTENT(OUT) ::chi2
13330 INTEGER(mpl),
INTENT(OUT) ::ndf
13331 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