921 REAL,
DIMENSION(2) :: ta
924 INTEGER(mpi) :: iopnmp
935 INTEGER(mpi) :: nsecnd
936 INTEGER(mpi) :: ntsec
938 CHARACTER (LEN=24) :: chdate
939 CHARACTER (LEN=24) :: chost
941 CHARACTER (LEN=6) :: c6
942 INTEGER major, minor, patch
966 CALL getenv(
'HOSTNAME',chost)
967 IF (chost(1:1) ==
' ')
CALL getenv(
'HOST',chost)
968 WRITE(*,*)
'($Id: c5a7342b3793f36f30ad6658bc1f72bc74250cf7 $)'
973 CALL ilaver( major,minor, patch )
974 WRITE(*,110) lapack64, major,minor, patch
975110
FORMAT(
' using LAPACK64 with ',(a),
', version ',i0,
'.',i0,
'.',i0)
977 WRITE(*,*)
'using Intel oneMKL PARDISO'
981 WRITE(*,111) __gnuc__ , __gnuc_minor__ , __gnuc_patchlevel__
982111
FORMAT(
' compiled with gcc ',i0,
'.',i0,
'.',i0)
985 WRITE(*,111) __pgic__ , __pgic_minor__ , __pgic_patchlevel__
986111
FORMAT(
' compiled with pgi ',i0,
'.',i0,
'.',i0)
989 WRITE(*,*)
' < Millepede II-P starting ... ',chdate
993 WRITE(8,*)
'($Id: c5a7342b3793f36f30ad6658bc1f72bc74250cf7 $)'
995 WRITE(8,*)
'Log-file Millepede II-P ', chdate
996 WRITE(8,*)
' ', chost
998 CALL peend(-1,
'Still running or crashed')
1007 WRITE(*,*)
'!!! Checking input only, no calculation of a solution !!!'
1008 WRITE(8,*)
'!!! Checking input only, no calculation of a solution !!!'
1027 CALL getenv(
'OMP_NUM_THREADS',c6)
1029 CALL getenv(lapack64//
'_NUM_THREADS',c6)
1031 IF (c6(1:1) ==
' ')
THEN
1033 WRITE(*,*)
'Number of threads for LAPACK: unkown (empty OMP_NUM_THREADS)'
1035 WRITE(*,*)
'Number of threads for LAPACK: unkown (empty ',lapack64//
'_NUM_THREADS)'
1038 WRITE(*,*)
'Number of threads for LAPACK: ', c6
1058 CALL mvopen(lun,
'millepede.his')
1064 CALL mvopen(1,
'mpdebug.txt')
1068 times(0)=rstext-rstp
1074 times(1)=rloop1-rstext
1078 WRITE(8,*)
'Chi square cut equiv 3 st.dev applied ...'
1079 WRITE(8,*)
' in first iteration with factor',
chicut
1080 WRITE(8,*)
' in second iteration with factor',
chirem
1081 WRITE(8,*)
' (reduced by sqrt in next iterations)'
1085 WRITE(8,*)
'Down-weighting of outliers in',
lhuber,
' iterations'
1086 WRITE(8,*)
'Cut on downweight fraction',
dwcut
1090 times(2)=rloop2-rloop1
1095 CALL peend(5,
'Ended without solution (empty constraints)')
1097 CALL peend(0,
'Ended normally')
1134 CALL gmpdef(6,1,
'log10(#records) vs file number')
1135 CALL gmpdef(7,1,
'final rejection fraction vs file number')
1137 'final <Chi^2/Ndf> from accepted local fits vs file number')
1138 CALL gmpdef(9,1,
'<Ndf> from accepted local fits vs file number')
1144 rej=real(nrc-
jfd(kfl),mps)/real(nrc,mps)
1145 CALL gmpxy(6,real(kfl,mps),log10(real(nrc,mps)))
1146 CALL gmpxy(7,real(kfl,mps),rej)
1148 IF (
jfd(kfl) > 0)
THEN
1149 c2ndf=
cfd(kfl)/real(
jfd(kfl),mps)
1150 CALL gmpxy(8,real(kfl,mps),c2ndf)
1151 andf=real(
dfd(kfl),mps)/real(
jfd(kfl),mps)
1152 CALL gmpxy(9,real(kfl,mps),andf)
1169 WRITE(*,*)
'Misalignment test wire chamber'
1172 CALL hmpdef( 9,-0.0015,+0.0015,
'True - fitted displacement')
1173 CALL hmpdef(10,-0.0015,+0.0015,
'True - fitted Vdrift')
1179 sums(1)=sums(1)+diff
1180 sums(2)=sums(2)+diff*diff
1182 sums(3)=sums(3)+diff
1183 sums(4)=sums(4)+diff*diff
1185 sums(1)=0.01_mpd*sums(1)
1186 sums(2)=sqrt(0.01_mpd*sums(2))
1187 sums(3)=0.01_mpd*sums(3)
1188 sums(4)=sqrt(0.01_mpd*sums(4))
1189 WRITE(*,143)
'Parameters 1 - 100: mean =',sums(1),
'rms =',sums(2)
1190 WRITE(*,143)
'Parameters 101 - 200: mean =',sums(3),
'rms =',sums(4)
1191143
FORMAT(6x,a28,f9.6,3x,a5,f9.6)
1194 WRITE(*,*)
' I label simulated fitted diff'
1195 WRITE(*,*)
' -------------------------------------------- '
1215 WRITE(*,*)
'Misalignment test Si tracker'
1218 CALL hmpdef( 9,-0.0025,+0.0025,
'True - fitted displacement X')
1219 CALL hmpdef(10,-0.025,+0.025,
'True - fitted displacement Y')
1230 sums(1)=sums(1)+1.0_mpd
1231 sums(2)=sums(2)+diff
1232 sums(3)=sums(3)+diff*diff
1237 err=sqrt(abs(gmati))
1239 sums(7)=sums(7)+1.0_mpd
1240 sums(8)=sums(8)+diff
1241 sums(9)=sums(9)+diff*diff
1244 IF (mod(i,3) == 1)
THEN
1248 sums(4)=sums(4)+1.0_mpd
1249 sums(5)=sums(5)+diff
1250 sums(6)=sums(6)+diff*diff
1255 err=sqrt(abs(gmati))
1257 sums(7)=sums(7)+1.0_mpd
1258 sums(8)=sums(8)+diff
1259 sums(9)=sums(9)+diff*diff
1264 sums(2)=sums(2)/sums(1)
1265 sums(3)=sqrt(sums(3)/sums(1))
1266 sums(5)=sums(5)/sums(4)
1267 sums(6)=sqrt(sums(6)/sums(4))
1268 WRITE(*,143)
'Parameters 1 - 500: mean =',sums(2),
'rms =',sums(3)
1269 WRITE(*,143)
'Parameters 501 - 700: mean =',sums(5),
'rms =',sums(6)
1270 IF (sums(7) > 0.5_mpd)
THEN
1271 sums(8)=sums(8)/sums(7)
1272 sums(9)=sqrt(sums(9)/sums(7))
1273 WRITE(*,143)
'Parameter pulls, all: mean =',sums(8),
'rms =',sums(9)
1277 WRITE(*,*)
' I label simulated fitted diff'
1278 WRITE(*,*)
' -------------------------------------------- '
1290 IF (mod(i,3) == 1)
THEN
1310 WRITE(8,*)
'Record',
nrec1,
' has largest residual:',
value1
1313 WRITE(8,*)
'Record',
nrec2,
' has largest Chi^2/Ndf:',
value2
1317 WRITE(8,*)
'Record',
nrec3,
' is first with error (rank deficit/NaN)'
1321 WRITE(8,*)
'In total 3 +',
nloopn,
' loops through the data files'
1323 WRITE(8,*)
'In total 2 +',
nloopn,
' loops through the data files'
1327 WRITE(8,*)
'In total ',
mnrsit,
' internal MINRES iterations'
1335 ntsec=nint(deltat,mpi)
1336 CALL sechms(deltat,nhour,minut,secnd)
1337 nsecnd=nint(secnd,mpi)
1338 WRITE(8,*)
'Total time =',ntsec,
' seconds =',nhour,
' h',minut, &
1339 ' m',nsecnd,
' seconds'
1341 WRITE(8,*)
'end ', chdate
1348 IF(sum(
nrejec) /= 0)
THEN
1350 WRITE(8,*)
'Data records rejected in last iteration: '
1357 WRITE(*,*)
' < Millepede II-P ending ... ', chdate
1364106
FORMAT(
' PARDISO dyn. memory allocation: ',f11.6,
' GB')
1369102
FORMAT(2x,i4,i10,2x,3f10.5)
1370103
FORMAT(
' Times [in sec] for text processing',f12.3/ &
1373 ' func. value ',f12.3,
' *',f4.0/ &
1374 ' func. value, global matrix, solution',f12.3,
' *',f4.0/ &
1375 ' new solution',f12.3,
' *',f4.0/)
1376105
FORMAT(
' Peak dynamic memory allocation: ',f11.6,
' GB')
1396 INTEGER(mpi) :: istop
1397 INTEGER(mpi) :: itgbi
1398 INTEGER(mpi) :: itgbl
1400 INTEGER(mpi) :: itnlim
1401 INTEGER(mpi) :: nout
1403 INTEGER(mpi),
INTENT(IN) :: ivgbi
1440 globalcorrections, itnlim, nout, rtol, istop, itn, anorm, acond, rnorm, arnorm, ynorm)
1444 globalcorrections, itnlim, nout, rtol, istop, itn, anorm, acond, rnorm, arnorm, ynorm)
1447 globalcorrections, itnlim, nout, rtol, istop, itn, anorm, acond, rnorm, arnorm, ynorm)
1453 err=sqrt(abs(real(gmati,mps)))
1454 IF(gmati < 0.0_mpd) err=-err
1455 diag=matij(ivgbi,ivgbi)
1456 gcor2=real(1.0_mpd-1.0_mpd/(gmati*diag),mps)
1458101
FORMAT(1x,
' label parameter presigma differ', &
1459 ' Error gcor^2 iit'/ 1x,
'---------',2x,5(
'-----------'),2x,
'----')
1460102
FORMAT(i10,2x,4g12.4,f7.4,i6,i4)
1480 INTEGER(mpi) :: istop
1481 INTEGER(mpi) :: itgbi
1482 INTEGER(mpi) :: itgbl
1484 INTEGER(mpi) :: itnlim
1485 INTEGER(mpi) :: nout
1487 INTEGER(mpi),
INTENT(IN) :: ivgbi
1516 mxxnrm = real(
nagb,mpd)/sqrt(epsilon(mxxnrm))
1518 trcond = 1.0_mpd/epsilon(trcond)
1519 ELSE IF(
mrmode == 2)
THEN
1527 itnlim=itnlim, rtol=rtol, maxxnorm=mxxnrm, trancond=trcond, &
1531 itnlim=itnlim, rtol=rtol, maxxnorm=mxxnrm, trancond=trcond, &
1535 itnlim=itnlim, rtol=rtol, maxxnorm=mxxnrm, trancond=trcond, &
1542 err=sqrt(abs(real(gmati,mps)))
1543 IF(gmati < 0.0_mpd) err=-err
1544 diag=matij(ivgbi,ivgbi)
1545 gcor2=real(1.0_mpd-1.0_mpd/(gmati*diag),mps)
1547101
FORMAT(1x,
' label parameter presigma differ', &
1548 ' Error gcor^2 iit'/ 1x,
'---------',2x,5(
'-----------'),2x,
'----')
1549102
FORMAT(i10,2x,4g12.4,f7.4,i6,i4)
1562 INTEGER(mpi) :: icgb
1563 INTEGER(mpi) :: irhs
1564 INTEGER(mpi) :: itgbi
1565 INTEGER(mpi) :: ivgb
1567 INTEGER(mpi) :: jcgb
1569 INTEGER(mpi) :: label
1571 INTEGER(mpi) :: inone
1574 REAL(mpd) :: drhs(4)
1575 INTEGER(mpi) :: idrh (4)
1600 IF(abs(rhs) > climit)
THEN
1606 WRITE(*,101) (idrh(l),drhs(l),l=1,irhs)
1615 WRITE(*,101) (idrh(l),drhs(l),l=1,irhs)
1618 WRITE(*,102)
' Constraints: only equation values >', climit,
' are printed'
1619101
FORMAT(
' ',4(i6,g11.3))
1620102
FORMAT(a,g11.2,a)
1633 INTEGER(mpi) :: icgb
1634 INTEGER(mpi) :: icgrp
1635 INTEGER(mpi) :: ioff
1636 INTEGER(mpi) :: itgbi
1638 INTEGER(mpi) :: jcgb
1639 INTEGER(mpi) :: label
1640 INTEGER(mpi) :: labelf
1641 INTEGER(mpi) :: labell
1642 INTEGER(mpi) :: last
1643 INTEGER(mpi) :: line1
1644 INTEGER(mpi) :: ncon
1645 INTEGER(mpi) :: ndiff
1646 INTEGER(mpi) :: npar
1647 INTEGER(mpi) :: inone
1648 INTEGER(mpi) :: itype
1649 INTEGER(mpi) :: ncgbd
1650 INTEGER(mpi) :: ncgbr
1651 INTEGER(mpi) :: ncgbw
1652 INTEGER(mpi) :: ncgrpd
1653 INTEGER(mpi) :: ncgrpr
1654 INTEGER(mpi) :: next
1656 INTEGER(mpl):: length
1657 INTEGER(mpl) :: rows
1659 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: vecParConsOffsets
1660 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: vecParConsList
1661 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: vecConsParOffsets
1662 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: vecConsParList
1663 INTEGER(mpi),
DIMENSION(:,:),
ALLOCATABLE :: matConsGroupIndex
1676 IF(last < 0.AND.label < 0)
THEN
1679 IF(itype == 2) ncgbw=ncgbw+1
1686 IF(label > 0.AND.itype == 2)
THEN
1693 IF (ncgbw == 0)
THEN
1694 WRITE(*,*)
'GRPCON:',
ncgb,
' constraints found in steering files'
1696 WRITE(*,*)
'GRPCON:',
ncgb,
' constraints found in steering files,',ncgbw,
'weighted'
1701 length=
ncgb+1; rows=3
1714 CALL mpalloc(matconsgroupindex,rows,length,
'group index for constraint (I)')
1717 CALL mpalloc(vecconsparoffsets,length,
'offsets for global par list for cons. (I)')
1719 CALL mpalloc(vecparconsoffsets,length,
'offsets for cons. list for global par. (I)')
1720 vecparconsoffsets(1)=0
1722 vecparconsoffsets(i+1)=vecparconsoffsets(i)+
globalparcons(i)
1726 length=vecparconsoffsets(
ntgb+1)
1727 CALL mpalloc(vecconsparlist,length,
'global par. list for constraint (I)')
1728 CALL mpalloc(vecparconslist,length,
'constraint list for global par. (I)')
1733 vecconsparoffsets(1)=ioff
1745 vecparconslist(vecparconsoffsets(itgbi)+
globalparcons(itgbi))=icgb
1747 vecconsparlist(ioff+npar)=itgbi
1753 CALL sort1k(vecconsparlist(ioff+1),npar)
1757 next=vecconsparlist(ioff+j)
1758 IF (next /= last)
THEN
1760 vecconsparlist(ioff+ndiff) = next
1769 vecconsparoffsets(icgb+1)=ioff
1782 print *,
' Constraint #parameters first par. last par. first line'
1790 npar=vecconsparoffsets(icgb+1)-vecconsparoffsets(icgb)
1794 print *, jcgb, npar, labelf, labell, line1
1797 icgrp=matconsgroupindex(1,icgb)
1798 IF (icgrp == 0)
THEN
1800 DO i=vecconsparoffsets(icgb)+1, vecconsparoffsets(icgb+1)
1801 itgbi=vecconsparlist(i)
1803 DO j=vecparconsoffsets(itgbi)+1,vecparconsoffsets(itgbi+1)
1804 icgrp=matconsgroupindex(1,vecparconslist(j))
1810 IF (icgrp == 0)
THEN
1817 matconsgroupindex(2,icgb)=jcgb
1818 matconsgroupindex(3,icgb)=icgb
1819 DO i=vecconsparoffsets(icgb)+1, vecconsparoffsets(icgb+1)
1820 itgbi=vecconsparlist(i)
1823 DO j=vecparconsoffsets(itgbi)+1,vecparconsoffsets(itgbi+1)
1824 matconsgroupindex(1,vecparconslist(j))=icgrp
1828 WRITE(*,*)
'GRPCON:',
ncgrp,
' disjoint constraints groups built'
1836 icgb=matconsgroupindex(3,jcgb)
1841 icgrp=matconsgroupindex(1,jcgb)
1861 print *,
' cons.group first con. first par. last par. #cons #par'
1872 print *, icgrp,
matconsgroups(1,icgrp), labelf, labell, ncon, npar
1875 IF (ncon == npar)
THEN
1882 print *, icgrp,
matconsgroups(1,icgrp), labelf, labell,
' : cons.group resolved'
1897 print *, icgrp,
matconsgroups(1,icgrp), labelf, labell,
' : cons.group redundant'
1902 IF (ncgrpr > 0)
THEN
1903 WRITE(*,*)
'GRPCON:',ncgbr,
' redundancy constraints in ', ncgrpr,
' groups resolved'
1907 IF (ncgrpd > 0)
THEN
1908 WRITE(*,*)
'GRPCON:',ncgbd,
' redundancy constraints in ', ncgrpd,
' groups detected'
1931 INTEGER(mpi) :: icgb
1932 INTEGER(mpi) :: icgrp
1933 INTEGER(mpi) :: ifrst
1934 INTEGER(mpi) :: ilast
1935 INTEGER(mpi) :: isblck
1936 INTEGER(mpi) :: itgbi
1937 INTEGER(mpi) :: ivgb
1939 INTEGER(mpi) :: jcgb
1940 INTEGER(mpi) :: jfrst
1941 INTEGER(mpi) :: label
1942 INTEGER(mpi) :: labelf
1943 INTEGER(mpi) :: labell
1944 INTEGER(mpi) :: ncon
1945 INTEGER(mpi) :: ngrp
1946 INTEGER(mpi) :: npar
1947 INTEGER(mpi) :: ncnmxb
1948 INTEGER(mpi) :: ncnmxg
1949 INTEGER(mpi) :: nprmxb
1950 INTEGER(mpi) :: nprmxg
1951 INTEGER(mpi) :: inone
1952 INTEGER(mpi) :: nvar
1954 INTEGER(mpl):: length
1955 INTEGER(mpl) :: rows
1957 INTEGER(mpi),
DIMENSION(:,:),
ALLOCATABLE :: matConsGroupIndex
1970 length=
ncgrp+1; rows=3
1975 CALL mpalloc(matconsgroupindex,rows,length,
'group index for constraint (I)')
2012 IF (nvar > 0 .OR.
iskpec == 0)
THEN
2015 matconsgroupindex(1,
ncgb)=ngrp+1
2016 matconsgroupindex(2,
ncgb)=icgb
2017 matconsgroupindex(3,
ncgb)=nvar
2020 IF (
ncgb > ncon) ngrp=ngrp+1
2026 WRITE(*,*)
'PRPCON:',
ncgbe,
' empty constraints skipped'
2028 WRITE(*,*)
'PRPCON:',
ncgbe,
' empty constraints detected, to be fixed !!!'
2029 WRITE(*,*)
' (use option "skipemptycons" to skip those)'
2034 WRITE(*,*)
'!!! Switch to "-C" (checking input only), no calculation of a solution !!!'
2035 WRITE(8,*)
'!!! Switch to "-C" (checking input only), no calculation of a solution !!!'
2040 WRITE(*,*)
'PRPCON:',
ncgb,
' constraints accepted'
2043 IF(
ncgb == 0)
RETURN
2050 icgb=matconsgroupindex(2,jcgb)
2055 icgrp=matconsgroupindex(1,jcgb)
2081 WRITE(*,*)
' Cons. sorted index #var.par. first line first label last label'
2082 WRITE(*,*)
' Cons. group index first cons. last cons. first label last label'
2083 WRITE(*,*)
' Cons. block index first group last group first label last label'
2089 nvar=matconsgroupindex(3,jcgb)
2093 WRITE(*,*)
' Cons. sorted', jcgb, nvar, &
2096 WRITE(*,*)
' Cons. sorted', jcgb,
' empty (0)', &
2105 WRITE(*,*)
' Cons. group ', icgrp,
matconsgroups(1,icgrp), &
2118 DO i=icgrp,isblck,-1
2132 WRITE(*,*)
' Cons. block ',
ncblck, isblck, icgrp, labelf, labell
2150 ncnmxb=max(ncnmxb,ncon)
2151 nprmxb=max(nprmxb,npar)
2176 ncnmxg=max(ncnmxg,ncon)
2177 nprmxg=max(nprmxg,npar)
2212 WRITE(*,*)
'PRPCON: constraints split into ',
ncgrp,
'(disjoint) groups,'
2213 WRITE(*,*)
' groups combined into ',
ncblck,
'(non overlapping) blocks'
2214 WRITE(*,*)
' max group size (cons., par.) ', ncnmxg, nprmxg
2215 WRITE(*,*)
' max block size (cons., par.) ', ncnmxb, nprmxb
2233 INTEGER(mpi) :: icgb
2234 INTEGER(mpi) :: icgrp
2236 INTEGER(mpi) :: ifirst
2237 INTEGER(mpi) :: ilast
2238 INTEGER(mpl) :: ioffc
2239 INTEGER(mpl) :: ioffp
2240 INTEGER(mpi) :: irank
2241 INTEGER(mpi) :: ipar0
2242 INTEGER(mpi) :: itgbi
2243 INTEGER(mpi) :: ivgb
2245 INTEGER(mpi) :: jcgb
2247 INTEGER(mpi) :: label
2248 INTEGER(mpi) :: ncon
2249 INTEGER(mpi) :: npar
2250 INTEGER(mpi) :: nrank
2251 INTEGER(mpi) :: inone
2256 INTEGER(mpl):: length
2257 REAL(mpd),
DIMENSION(:),
ALLOCATABLE :: matConstraintsT
2258 REAL(mpd),
DIMENSION(:),
ALLOCATABLE :: auxVectorD
2259 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: auxVectorI
2263 IF(
ncgb == 0)
RETURN
2272 CALL mpalloc(auxvectori,length,
'auxiliary array (I)')
2273 CALL mpalloc(auxvectord,length,
'auxiliary array (D)')
2276 CALL mpalloc(matconstraintst,length,
'transposed matrix of constraints (blocks)')
2277 matconstraintst=0.0_mpd
2290 WRITE(*,*)
' Constraint group, #con, rank', icgrp, ncon, 0,
' (empty)'
2293 DO jcgb=ifirst,ilast
2305 IF(ivgb > 0) matconstraintst(int(jcgb-ifirst,mpl)*int(npar,mpl)+ivgb-ipar0+ioffc)= &
2306 matconstraintst(int(jcgb-ifirst,mpl)*int(npar,mpl)+ivgb-ipar0+ioffc)+factr
2313 DO ll=ioffc+1,ioffc+npar
2319 matconstraintst(int(i-1,mpl)*int(npar,mpl)+ll)* &
2320 matconstraintst(int(j-1,mpl)*int(npar,mpl)+ll)
2326 IF (
icheck > 1 .OR. irank < ncon)
THEN
2327 WRITE(*,*)
' Constraint group, #con, rank', icgrp, ncon, irank
2328 IF (irank < ncon)
THEN
2329 WRITE(*,*)
' .. rank deficit !! '
2330 WRITE(*,*)
' E.g. fix all parameters and remove all constraints related to label ', &
2335 ioffc=ioffc+int(npar,mpl)*int(ncon,mpl)
2342 WRITE(*,*)
'Rank of product matrix of constraints is',nrank, &
2343 ' for',
ncgb,
' constraint equations'
2344 WRITE(8,*)
'Rank of product matrix of constraints is',nrank, &
2345 ' for',
ncgb,
' constraint equations'
2346 IF(nrank <
ncgb)
THEN
2347 WRITE(*,*)
'Warning: insufficient constraint equations!'
2348 WRITE(8,*)
'Warning: insufficient constraint equations!'
2351 WRITE(*,*)
' --> enforcing SUBITO mode'
2352 WRITE(8,*)
' --> enforcing SUBITO mode'
2359 print *,
'QL decomposition of constraints matrix'
2362 WRITE(
lunlog,*)
'QL decomposition of constraints matrix'
2374 CALL lpqldec(matconstraintst,evmin,evmax)
2378 print *,
' largest |eigenvalue| of L: ', evmax
2379 print *,
' smallest |eigenvalue| of L: ', evmin
2380 IF (evmin == 0.0_mpd.AND.
icheck == 0)
THEN
2381 CALL peend(27,
'Aborted, singular QL decomposition of constraints matrix')
2382 stop
'FEASMA: stopping due to singular QL decomposition of constraints matrix'
2408 INTEGER(mpi) :: icgb
2409 INTEGER(mpi) :: icgrp
2410 INTEGER(mpi) :: iter
2411 INTEGER(mpi) :: itgbi
2412 INTEGER(mpi) :: ivgb
2413 INTEGER(mpi) :: ieblck
2414 INTEGER(mpi) :: isblck
2415 INTEGER(mpi) :: ifirst
2416 INTEGER(mpi) :: ilast
2418 INTEGER(mpi) :: jcgb
2419 INTEGER(mpi) :: label
2420 INTEGER(mpi) :: inone
2421 INTEGER(mpi) :: ncon
2423 REAL(mps),
INTENT(IN) :: concut
2424 INTEGER(mpi),
INTENT(OUT) :: iact
2431 REAL(mpd),
DIMENSION(:),
ALLOCATABLE :: vecCorrections
2435 IF(
ncgb == 0)
RETURN
2465 sum1=sqrt(sum1/real(
ncgb,mpd))
2466 sum2=sum2/real(
ncgb,mpd)
2468 IF(iter == 1.AND.sum1 < concut)
RETURN
2470 IF(iter == 1.AND.
ncgb <= 12)
THEN
2472 WRITE(*,*)
'Constraint equation discrepancies:'
2474101
FORMAT(4x,4(i5,g12.4))
2476103
FORMAT(10x,
' Cut on rms value is',g8.1)
2481 WRITE(*,*)
'Improve constraints'
2485 WRITE(*,102) iter,sum1,sum2,sum3
2486102
FORMAT(i6,
' rms',g12.4,
' avrg_abs',g12.4,
' max_abs',g12.4)
2488 CALL mpalloc(veccorrections,int(
nvgb,mpl),
'constraint corrections')
2489 veccorrections=0.0_mpd
2497 ieblck=isblck+(ncon*(ncon+1))/2
2566 INTEGER(mpi) :: iact
2567 INTEGER(mpi) :: ierrc
2568 INTEGER(mpi) :: ierrf
2569 INTEGER(mpi) :: ioffp
2571 INTEGER(mpi) :: ithr
2572 INTEGER(mpi) :: jfile
2573 INTEGER(mpi) :: jrec
2575 INTEGER(mpi) :: kfile
2578 INTEGER(mpi) :: mpri
2580 INTEGER(mpi) :: nact
2581 INTEGER(mpi) :: nbuf
2582 INTEGER(mpi) :: ndata
2583 INTEGER(mpi) :: noff
2584 INTEGER(mpi) :: noffs
2585 INTEGER(mpi) :: npointer
2586 INTEGER(mpi) :: npri
2590 INTEGER(mpi) :: nrpr
2591 INTEGER(mpi) :: nthr
2592 INTEGER(mpi) :: ntot
2593 INTEGER(mpi) :: maxRecordSize
2594 INTEGER(mpi) :: maxRecordFile
2596 INTEGER(mpi),
INTENT(OUT) :: more
2606 CHARACTER (LEN=7) :: cfile
2611 SUBROUTINE readc(bufferD, bufferF, bufferI, bufferLength, lun, err)
BIND(c)
2613 REAL(c_double),
DIMENSION(*),
INTENT(OUT) :: bufferD
2614 REAL(c_float),
DIMENSION(*),
INTENT(OUT) :: bufferF
2615 INTEGER(c_int),
DIMENSION(*),
INTENT(OUT) :: bufferI
2616 INTEGER(c_int),
INTENT(INOUT) :: bufferLength
2617 INTEGER(c_int),
INTENT(IN),
VALUE :: lun
2618 INTEGER(c_int),
INTENT(OUT) :: err
2619 END SUBROUTINE readc
2625 DATA npri / 0 /, mpri / 1000 /
2674 files:
DO WHILE (jfile > 0)
2678 CALL binopn(kfile,ithr,ios)
2684 IF(kfile <=
nfilf)
THEN
2705 eof=(ierrc <= 0.AND.ierrc /= -4)
2706 IF(eof.AND.ierrc < 0)
THEN
2707 WRITE(*,*)
'Read error for binary Cfile', kfile,
'record', jrec+1,
':', ierrc
2708 WRITE(8,*)
'Read error for binary Cfile', kfile,
'record', jrec+1,
':', ierrc
2710 WRITE(cfile,
'(I7)') kfile
2711 CALL peend(18,
'Aborted, read error(s) for binary file ' // cfile)
2712 stop
'PEREAD: stopping due to read errors'
2714 IF (
kfd(1,jfile) == 1)
THEN
2720 IF(eof)
EXIT records
2725 xfd(jfile)=max(
xfd(jfile),n)
2742 IF ((noff+nr+2+
ndimbuf >= ndata*(iact+1)).OR.(nbuf >= npointer))
EXIT files
2757 IF (
kfd(1,jfile) == 1)
THEN
2758 print *,
'PEREAD: file ', kfile,
'read the first time, found',jrec,
' records'
2762 IF (-
kfd(1,jfile) /= jrec)
THEN
2763 WRITE(cfile,
'(I7)') kfile
2764 CALL peend(19,
'Aborted, binary file modified (length) ' // cfile)
2765 stop
'PEREAD: file modified (length)'
2816 IF (
kfd(1,jfile) == 1)
kfd(1,jfile)=-nrc
2835 DO WHILE (
nloopn == 0.AND.ntot >= nrpr)
2836 WRITE(*,*)
' Record ',nrpr
2837 IF (nrpr < 100000)
THEN
2846 IF (npri == 1)
WRITE(*,100)
2848100
FORMAT(/
' PeRead records active file' &
2849 /
' total block threads number')
2850101
FORMAT(
' PeRead',4i10)
2863 IF (
xfd(k) > maxrecordsize)
THEN
2864 maxrecordsize=
xfd(k)
2867 dw=real(-
kfd(1,k),mpd)
2870 ds1=ds1+dw*real(
wfd(k),mpd)
2871 ds2=ds2+dw*real(
wfd(k)**2,mpd)
2873 print *,
'PEREAD: file ', maxrecordfile,
'with max record size ', maxrecordsize
2874 IF (
nfilw > 0.AND.ds0 > 0.0_mpd)
THEN
2878 WRITE(lun,177)
nfilw,real(ds1,mps),real(ds2,mps)
2879177
FORMAT(/
' !!!!!',i4,
' weighted binary files', &
2880 /
' !!!!! mean, variance of weights =',2g12.4)
2898179
FORMAT(/
' Read cache usage (#blocks, #records, ', &
2899 'min,max records/block'/17x,i10,i12,2i10)
2917 INTEGER(mpi),
INTENT(IN) :: mode
2919 INTEGER(mpi) :: ibuf
2920 INTEGER(mpi) :: ichunk
2922 INTEGER(mpi) :: itgbi
2928 INTEGER(mpi),
PARAMETER :: maxbad = 100
2929 INTEGER(mpi) :: nbad
2930 INTEGER(mpi) :: nerr
2931 INTEGER(mpi) :: inone
2949 CALL isjajb(nst,ist,ja,jb,jsp)
2971 CALL pechk(ibuf,nerr)
2974 IF(nbad >= maxbad)
EXIT
2979 CALL isjajb(nst,ist,ja,jb,jsp)
2992 CALL peend(20,
'Aborted, bad binary records')
2993 stop
'PEREAD: stopping due to bad records'
3014 INTEGER(mpi) :: ioff
3021 INTEGER(mpi),
INTENT(IN) :: ibuf
3022 INTEGER(mpi),
INTENT(OUT) :: nerr
3031 outer:
DO WHILE(is < nst)
3036 IF(is > nst)
EXIT outer
3042 IF(is > nst)
EXIT outer
3059100
FORMAT(
' PEREAD: record ', i8,
' in file ',i6,
' is broken !!!')
3069101
FORMAT(
' PEREAD: record ', i8,
' in file ',i6,
' contains ', i6,
' NaNs !!!')
3085 INTEGER(mpi) :: ibuf
3086 INTEGER(mpi) :: ichunk
3087 INTEGER(mpi) :: iproc
3088 INTEGER(mpi) :: ioff
3089 INTEGER(mpi) :: ioffbi
3091 INTEGER(mpi) :: itgbi
3096 INTEGER(mpi) :: nalg
3097 INTEGER(mpi) :: neqna
3100 INTEGER(mpi) :: nzero
3101 INTEGER(mpi) :: inone
3102 INTEGER(mpl) :: length
3137 CALL isjajb(nst,ist,ja,jb,jsp)
3168 CALL isjajb(nst,ist,ja,jb,jsp)
3203 CALL pargrp(ist+1,ist+nnz)
3232 INTEGER(mpi) :: istart
3233 INTEGER(mpi) :: itgbi
3235 INTEGER(mpi) :: jstart
3236 INTEGER(mpi) :: jtgbi
3237 INTEGER(mpi) :: lstart
3238 INTEGER(mpi) :: ltgbi
3240 INTEGER(mpi),
INTENT(IN) :: inds
3241 INTEGER(mpi),
INTENT(IN) :: inde
3243 IF (inds > inde)
RETURN
3252 IF (istart == 0)
THEN
3253 IF (itgbi /= ltgbi+1)
THEN
3256 IF (lstart == 0)
THEN
3271 IF (istart /= jstart)
THEN
3282 IF (itgbi /= ltgbi+1)
THEN
3297 IF (istart /= jstart)
THEN
3349 INTEGER(mpi),
INTENT(IN) :: nst
3350 INTEGER(mpi),
INTENT(IN OUT) :: is
3351 INTEGER(mpi),
INTENT(OUT) :: ja
3352 INTEGER(mpi),
INTENT(OUT) :: jb
3353 INTEGER(mpi),
INTENT(OUT) :: jsp
3361 IF(is >= nst)
RETURN
3415 INTEGER(mpi) :: ioffb
3417 INTEGER(mpi) :: itgbi
3418 INTEGER(mpi) :: itgbij
3419 INTEGER(mpi) :: itgbik
3420 INTEGER(mpi) :: ivgb
3421 INTEGER(mpi) :: ivgbij
3422 INTEGER(mpi) :: ivgbik
3425 INTEGER(mpi) :: lastit
3427 INTEGER(mpi) :: ncrit
3428 INTEGER(mpi) :: ngras
3429 INTEGER(mpi) :: nparl
3431 INTEGER(mpl) :: nrej
3432 INTEGER(mpi) :: inone
3433 INTEGER(mpi) :: ilow
3434 INTEGER(mpi) :: nlow
3435 INTEGER(mpi) :: nzero
3455 CALL gmpdef(1,4,
'Function value in iterations')
3457 CALL gmpdef(2,3,
'Number of MINRES steps vs iteration nr')
3459 CALL hmpdef( 5,0.0,0.0,
'Number of degrees of freedom')
3460 CALL hmpdef(11,0.0,0.0,
'Number of local parameters')
3461 CALL hmpdef(16,0.0,24.0,
'LOG10(cond(band part decomp.)) local fit ')
3462 CALL hmpdef(23,0.0,0.0,
'SQRT of diagonal elements without presigma')
3463 CALL hmpdef(24,0.0,0.0,
'Log10 of off-diagonal elements')
3464 CALL hmpdef(25,0.0,0.0,
'Relative individual pre-sigma')
3465 CALL hmpdef(26,0.0,0.0,
'Relative global pre-sigma')
3470 'Normalized residuals of single (global) measurement')
3472 'Normalized residuals of single (local) measurement')
3474 'Pulls of single (global) measurement')
3476 'Pulls of single (local) measurement')
3477 CALL hmpdef(4,0.0,0.0,
'Chi^2/Ndf after local fit')
3478 CALL gmpdef(4,5,
'location, dispersion (res.) vs record nr')
3479 CALL gmpdef(5,5,
'location, dispersion (pull) vs record nr')
3493 IF(
nloopn == 2)
CALL hmpdef(6,0.0,0.0,
'Down-weight fraction')
3496 IF(
iterat /= lastit)
THEN
3504 ELSE IF(
iterat >= 1)
THEN
3552 WRITE(*,111) nparl,ncrit,usei,used,peaki,peakd
3553111
FORMAT(
' Write cache usage (#flush,#overrun,<levels>,', &
3554 'peak(levels))'/2i7,
',',4(f6.1,
'%'))
3562 ELSE IF (
mbandw > 0)
THEN
3594 print *,
" ... warning ..."
3595 print *,
" global parameters with too few (< MREQENA) accepted entries: ", nlow
3599 IF(
icalcm == 1 .AND. nzero > 0)
THEN
3601 WRITE(*,*)
'Warning: the rank defect of the symmetric',
nfgb, &
3602 '-by-',
nfgb,
' matrix is ',
ndefec,
' (should be zero).'
3603 WRITE(lun,*)
'Warning: the rank defect of the symmetric',
nfgb, &
3604 '-by-',
nfgb,
' matrix is ',
ndefec,
' (should be zero).'
3607 WRITE(*,*)
' --> enforcing SUBITO mode'
3608 WRITE(lun,*)
' --> enforcing SUBITO mode'
3618 elmt=real(matij(i,i),mps)
3619 IF(elmt > 0.0)
CALL hmpent(23,1.0/sqrt(elmt))
3653 CALL addsums(1, adder, 0, 1.0_mpl)
3667 IF(sgm > 0.0) weight=1.0/sgm**2
3683 IF(itgbij /= 0)
THEN
3689 IF (factrj == 0.0_mpd) cycle
3699 IF(
icalcm == 1.AND.ivgbij > 0)
THEN
3706 IF(ivgbij > 0.AND.ivgbik > 0)
THEN
3707 CALL mupdat(ivgbij,ivgbik,weight*factrj*factrk)
3713 adder=weight*dsum**2
3714 CALL addsums(1, adder, 1, 1.0_mpl)
3732 ratae=max(-99.9,ratae)
3741 WRITE(*,*)
'Data records rejected in initial loop:'
3786 WRITE(lun,*)
' === local fits have bordered band matrix structure ==='
3787 IF (
nbndr(1) > 0 )
WRITE(lun,101)
' NBNDR',
nbndr(1),
'number of records (upper/left border)'
3788 IF (
nbndr(2) > 0 )
WRITE(lun,101)
' NBNDR',
nbndr(2),
'number of records (lower/right border)'
3789 WRITE(lun,101)
' NBDRX',
nbdrx,
'max border size'
3790 WRITE(lun,101)
' NBNDX',
nbndx,
'max band width'
3799101
FORMAT(1x,a8,
' =',i14,
' = ',a)
3814 INTEGER(mpi),
INTENT(IN) :: lunp
3819101
FORMAT(
' it fc',
' fcn_value dfcn_exp slpr costh iit st', &
3820 ' ls step cutf',1x,
'rejects hhmmss FMS')
3821102
FORMAT(
' -- --',
' ----------- -------- ---- ----- --- --', &
3822 ' -- ----- ----',1x,
'------- ------ ---')
3838 INTEGER(mpl) :: nrej
3839 INTEGER(mpi) :: nsecnd
3843 REAL(mps) :: slopes(3)
3844 REAL(mps) :: steps(3)
3845 REAL,
DIMENSION(2) :: ta
3848 INTEGER(mpi),
INTENT(IN) :: lunp
3850 CHARACTER (LEN=4):: ccalcm(4)
3851 DATA ccalcm /
' end',
' S',
' F ',
' FMS' /
3855 IF(nrej > 9999999) nrej=9999999
3859 nsecnd=nint(secnd,mpi)
3868 CALL ptlopt(nfa,ma,slopes,steps)
3869 ratae=max(-99.9,min(99.9,slopes(2)/slopes(1)))
3875103
FORMAT(i3,i3,e12.5,38x,f5.1, 1x,i7, i3,i2.2,i2.2,a4)
3876104
FORMAT(i3,i3,e12.5,1x,e8.2,f6.3,f6.3,i5,2i3,f6.3,f5.1, &
3877 1x,i7, i3,i2.2,i2.2,a4)
3878105
FORMAT(i3,i3,e12.5,1x,e8.2,12x,i5,i3,9x,f5.1, &
3879 1x,i7, i3,i2.2,i2.2,a4)
3892 INTEGER(mpi) :: minut
3894 INTEGER(mpi) :: nhour
3895 INTEGER(mpl) :: nrej
3896 INTEGER(mpi) :: nsecnd
3900 REAL(mps) :: slopes(3)
3901 REAL(mps) :: steps(3)
3902 REAL,
DIMENSION(2) :: ta
3905 INTEGER(mpi),
INTENT(IN) :: lunp
3906 CHARACTER (LEN=4):: ccalcm(4)
3907 DATA ccalcm /
' end',
' S',
' F ',
' FMS' /
3911 IF(nrej > 9999999) nrej=9999999
3915 nsecnd=nint(secnd,mpi)
3919 CALL ptlopt(nfa,ma,slopes,steps)
3920 ratae=abs(slopes(2)/slopes(1))
3925104
FORMAT(3x,i3,e12.5,9x, 35x, i7, i3,i2.2,i2.2,a4)
3926105
FORMAT(3x,i3,e12.5,9x, f6.3,14x,i3,f6.3,6x, i7, i3,i2.2,i2.2,a4)
3940 INTEGER(mpi) :: nsecnd
3943 REAL,
DIMENSION(2) :: ta
3946 INTEGER(mpi),
INTENT(IN) :: lunp
3947 CHARACTER (LEN=4):: ccalcm(4)
3948 DATA ccalcm /
' end',
' S',
' F ',
' FMS' /
3953 nsecnd=nint(secnd,mpi)
3955 WRITE(lunp,106) nhour,minut,nsecnd,ccalcm(
lcalcm)
3956106
FORMAT(69x,i3,i2.2,i2.2,a4)
3966 INTEGER(mpi) :: lunit
3968 WRITE(lunit,102)
'Explanation of iteration table'
3969 WRITE(lunit,102)
'=============================='
3970 WRITE(lunit,101)
'it', &
3971 'iteration number. Global parameters are improved for it > 0.'
3972 WRITE(lunit,102)
'First function evaluation is called iteraton 0.'
3973 WRITE(lunit,101)
'fc',
'number of function evaluations.'
3974 WRITE(lunit,101)
'fcn_value',
'value of 2 x Likelihood function (LF).'
3975 WRITE(lunit,102)
'The final value is the chi^2 value of the fit and should'
3976 WRITE(lunit,102)
'be about equal to the NDF (see below).'
3977 WRITE(lunit,101)
'dfcn_exp', &
3978 'expected reduction of the value of the Likelihood function (LF)'
3979 WRITE(lunit,101)
'slpr',
'ratio of the actual slope to inital slope.'
3980 WRITE(lunit,101)
'costh', &
3981 'cosine of angle between search direction and -gradient'
3983 WRITE(lunit,101)
'iit', &
3984 'number of internal iterations in MINRES algorithm'
3985 WRITE(lunit,101)
'st',
'stop code of MINRES algorithm'
3986 WRITE(lunit,102)
'< 0: rhs is very special, with beta2 = 0'
3987 WRITE(lunit,102)
'= 0: rhs b = 0, i.e. the exact solution is x = 0'
3988 WRITE(lunit,102)
'= 1 requested accuracy achieved, as determined by rtol'
3989 WRITE(lunit,102)
'= 2 reasonable accuracy achieved, given eps'
3990 WRITE(lunit,102)
'= 3 x has converged to an eigenvector'
3991 WRITE(lunit,102)
'= 4 matrix ill-conditioned (Acond has exceeded 0.1/eps)'
3992 WRITE(lunit,102)
'= 5 the iteration limit was reached'
3993 WRITE(lunit,102)
'= 6 Matrix x vector does not define a symmetric matrix'
3994 WRITE(lunit,102)
'= 7 Preconditioner does not define a symmetric matrix'
3995 ELSEIF (
metsol == 5)
THEN
3996 WRITE(lunit,101)
'iit', &
3997 'number of internal iterations in MINRES-QLP algorithm'
3998 WRITE(lunit,101)
'st',
'stop code of MINRES-QLP algorithm'
3999 WRITE(lunit,102)
'= 1: beta_{k+1} < eps, iteration k is the final Lanczos step.'
4000 WRITE(lunit,102)
'= 2: beta2 = 0. If M = I, b and x are eigenvectors of A.'
4001 WRITE(lunit,102)
'= 3: beta1 = 0. The exact solution is x = 0.'
4002 WRITE(lunit,102)
'= 4: A solution to (poss. singular) Ax = b found, given rtol.'
4003 WRITE(lunit,102)
'= 5: A solution to (poss. singular) Ax = b found, given eps.'
4004 WRITE(lunit,102)
'= 6: Pseudoinverse solution for singular LS problem, given rtol.'
4005 WRITE(lunit,102)
'= 7: Pseudoinverse solution for singular LS problem, given eps.'
4006 WRITE(lunit,102)
'= 8: The iteration limit was reached.'
4007 WRITE(lunit,102)
'= 9: The operator defined by Aprod appears to be unsymmetric.'
4008 WRITE(lunit,102)
'=10: The operator defined by Msolve appears to be unsymmetric.'
4009 WRITE(lunit,102)
'=11: The operator defined by Msolve appears to be indefinite.'
4010 WRITE(lunit,102)
'=12: xnorm has exceeded maxxnorm or will exceed it next iteration.'
4011 WRITE(lunit,102)
'=13: Acond has exceeded Acondlim or 0.1/eps.'
4012 WRITE(lunit,102)
'=14: Least-squares problem but no converged solution yet.'
4013 WRITE(lunit,102)
'=15: A null vector obtained, given rtol.'
4015 WRITE(lunit,101)
'ls',
'line search info'
4016 WRITE(lunit,102)
'< 0 recalculate function'
4017 WRITE(lunit,102)
'= 0: N or STP lt 0 or step not descending'
4018 WRITE(lunit,102)
'= 1: Linesearch convergence conditions reached'
4019 WRITE(lunit,102)
'= 2: interval of uncertainty at lower limit'
4020 WRITE(lunit,102)
'= 3: max nr of line search calls reached'
4021 WRITE(lunit,102)
'= 4: step at the lower bound'
4022 WRITE(lunit,102)
'= 5: step at the upper bound'
4023 WRITE(lunit,102)
'= 6: rounding error limitation'
4024 WRITE(lunit,101)
'step', &
4025 'the factor for the Newton step during the line search. Usually'
4027 'a value of 1 gives a sufficient reduction of the LF. Oherwise'
4028 WRITE(lunit,102)
'other step values are tried.'
4029 WRITE(lunit,101)
'cutf', &
4030 'cut factor. Local fits are rejected, if their chi^2 value'
4032 'is larger than the 3-sigma chi^2 value times the cut factor.'
4033 WRITE(lunit,102)
'A cut factor of 1 is used finally, but initially a larger'
4034 WRITE(lunit,102)
'factor may be used. A value of 0.0 means no cut.'
4035 WRITE(lunit,101)
'rejects',
'total number of rejected local fits.'
4036 WRITE(lunit,101)
'hmmsec',
'the time in hours (h), minutes (mm) and seconds.'
4037 WRITE(lunit,101)
'FMS',
'calculation of Function value, Matrix, Solution.'
4040101
FORMAT(a9,
' = ',a)
4057 INTEGER(mpi),
INTENT(IN) :: i
4058 INTEGER(mpi),
INTENT(IN) :: j
4059 REAL(mpd),
INTENT(IN) :: add
4061 INTEGER(mpl):: ijadd
4062 INTEGER(mpl):: ijcsr3
4067 IF(i <= 0.OR.j <= 0.OR. add == 0.0_mpd)
RETURN
4088 ELSE IF(
matsto == 2)
THEN
4116 CALL peend(23,
'Aborted, bad matrix index')
4117 stop
'mupdat: bad index'
4142 INTEGER(mpi),
INTENT(IN) :: i
4143 INTEGER(mpi),
INTENT(IN) :: j1
4144 INTEGER(mpi),
INTENT(IN) :: j2
4145 INTEGER(mpi),
INTENT(IN) :: il
4146 INTEGER(mpi),
INTENT(IN) :: jl
4147 INTEGER(mpi),
INTENT(IN) :: n
4148 REAL(mpd),
INTENT(IN) :: sub((n*n+n)/2)
4155 INTEGER(mpi):: iblast
4156 INTEGER(mpi):: iblock
4162 INTEGER(mpi):: jblast
4163 INTEGER(mpi):: jblock
4173 IF(i <= 0.OR.j1 <= 0.OR.j2 > i)
RETURN
4187 print *,
' MGUPDT: ij=0', i,j1,j2,il,jl,ij,lr,iprc,
matsto
4194 ijl=(k*k-k)/2+jl+ir-ia1
4211 ijl=(k*k-k)/2+jl+ir-ia1
4215 IF (jblock /= jblast.OR.iblock /= iblast)
THEN
4239 CALL ijpgrp(i,j1,ij,lr,iprc)
4241 print *,
' MGUPDT: ij=0', i,j1,j2,il,jl,ij,lr,iprc,
matsto
4305SUBROUTINE loopbf(nrej,numfil,naccf,chi2f,ndff)
4332 INTEGER(mpi) :: ibuf
4333 INTEGER(mpi) :: ichunk
4334 INTEGER(mpl) :: icmn
4335 INTEGER(mpl) :: icost
4337 INTEGER(mpi) :: idiag
4339 INTEGER(mpi) :: iext
4347 INTEGER(mpi) :: imeas
4350 INTEGER(mpi) :: ioffb
4351 INTEGER(mpi) :: ioffc
4352 INTEGER(mpi) :: ioffd
4353 INTEGER(mpi) :: ioffe
4354 INTEGER(mpi) :: ioffi
4355 INTEGER(mpi) :: ioffq
4356 INTEGER(mpi) :: iprc
4357 INTEGER(mpi) :: iprcnx
4358 INTEGER(mpi) :: iprdbg
4359 INTEGER(mpi) :: iproc
4360 INTEGER(mpi) :: irbin
4361 INTEGER(mpi) :: isize
4363 INTEGER(mpi) :: iter
4364 INTEGER(mpi) :: itgbi
4365 INTEGER(mpi) :: ivgbj
4366 INTEGER(mpi) :: ivgbk
4367 INTEGER(mpi) :: ivpgrp
4377 INTEGER(mpi) :: joffd
4378 INTEGER(mpi) :: joffi
4379 INTEGER(mpi) :: jproc
4383 INTEGER(mpi) :: kbdr
4384 INTEGER(mpi) :: kbdrx
4385 INTEGER(mpi) :: kbnd
4388 INTEGER(mpi) :: lvpgrp
4389 INTEGER(mpi) :: mbdr
4390 INTEGER(mpi) :: mbnd
4391 INTEGER(mpi) :: mside
4392 INTEGER(mpi) :: nalc
4393 INTEGER(mpi) :: nalg
4397 INTEGER(mpi) :: ndown
4399 INTEGER(mpi) :: nfred
4400 INTEGER(mpi) :: nfrei
4402 INTEGER(mpi) :: nprdbg
4403 INTEGER(mpi) :: nrank
4406 INTEGER(mpi) :: nter
4407 INTEGER(mpi) :: nweig
4408 INTEGER(mpi) :: ngrp
4409 INTEGER(mpi) :: npar
4411 INTEGER(mpl),
INTENT(IN OUT) :: nrej(6)
4412 INTEGER(mpi),
INTENT(IN) :: numfil
4413 INTEGER(mpi),
INTENT(IN OUT) :: naccf(numfil)
4414 REAL(mps),
INTENT(IN OUT) :: chi2f(numfil)
4415 INTEGER(mpi),
INTENT(IN OUT) :: ndff(numfil)
4425 INTEGER(mpi) :: ijprec
4432 CHARACTER (LEN=3):: chast
4433 DATA chuber/1.345_mpd/
4434 DATA cauchy/2.3849_mpd/
4482 IF(
nloopn == 1.AND.mod(nrc,100000_mpl) == 0)
THEN
4483 WRITE(*,*)
'Record',nrc,
' ... still reading'
4484 IF(
monpg1>0)
WRITE(
lunlog,*)
'Record',nrc,
' ... still reading'
4494 IF(nrc ==
nrecpr) lprnt=.true.
4495 IF(nrc ==
nrecp2) lprnt=.true.
4496 IF(nrc ==
nrecer) lprnt=.true.
4501 IF (nprdbg == 1) iprdbg=iproc
4502 IF (iproc /= iprdbg) lprnt=.false.
4507 WRITE(1,*)
'------------------ Loop',
nloopn, &
4508 ': Printout for record',nrc,iproc
4519 CALL isjajb(nst,ist,ja,jb,jsp)
4521 IF(imeas == 0)
WRITE(1,1121)
45261121
FORMAT(/
'Measured value and local derivatives'/ &
4527 ' i measured std_dev index...derivative ...')
45281122
FORMAT(i3,2g12.4,3(i3,g12.4)/(27x,3(i3,g12.4)))
4534 CALL isjajb(nst,ist,ja,jb,jsp)
4536 IF(imeas == 0)
WRITE(1,1123)
45481123
FORMAT(/
'Global derivatives'/ &
4549 ' i label gindex vindex derivative ...')
45501124
FORMAT(i3,2(i9,i7,i7,g12.4)/(3x,2(i9,i7,i7,g12.4)))
45511125
FORMAT(i3,2(i9,i7,i7,g12.4))
4561 WRITE(1,*)
'Data corrections using values of global parameters'
4562 WRITE(1,*)
'=================================================='
4572 CALL isjajb(nst,ist,ja,jb,jsp)
4604101
FORMAT(
' index measvalue corrvalue sigma')
4605102
FORMAT(i6,2x,2g12.4,
' +-',g12.4)
4607 IF(nalc <= 0)
GO TO 90
4609 ngg=(nalg*nalg+nalg)/2
4622 IF (ivpgrp /= lvpgrp)
THEN
4630 IF (npar /= nalg)
THEN
4631 print *,
' mismatch of number of global parameters ', nrc, nalg, npar, ngrp
4641 CALL peend(35,
'Aborted, mismatch of number of global parameters')
4642 stop
' mismatch of number of global parameters '
4669 DO WHILE(iter < nter)
4674 WRITE(1,*)
'Outlier-suppression iteration',iter,
' of',nter
4675 WRITE(1,*)
'=========================================='
4685 DO i=1,(nalc*nalc+nalc)/2
4696 wght =1.0_mpd/rerr**2
4701 IF(abs(resid) > chuber*rerr)
THEN
4702 wght=wght*chuber*rerr/abs(resid)
4706 wght=wght/(1.0+(resid/rerr/cauchy)**2)
4710 IF(lprnt.AND.iter /= 1.AND.nter /= 1)
THEN
4712 IF(abs(resid) > chuber*rerr) chast=
'* '
4713 IF(abs(resid) > 3.0*rerr) chast=
'** '
4714 IF(abs(resid) > 6.0*rerr) chast=
'***'
4715 IF(imeas == 0)
WRITE(1,*)
'Second loop: accumulate'
4716 IF(imeas == 0)
WRITE(1,103)
4721 WRITE(1,104) imeas,rmeas,resid,rerr,r1,chast,r2
4723103
FORMAT(
' index corrvalue residuum sigma', &
4725104
FORMAT(i6,2x,2g12.4,
' +-',g12.4,f7.2,1x,a3,f8.2)
4740 im=min(nalc+1-ij,nalc+1-ik)
4747 IF (iter == 1.AND.nalc > 5.AND.
lfitbb > 0)
THEN
4750 icmn=int(nalc,mpl)**3
4756 icost=6*int(nalc-kbdr,mpl)*int(kbnd+kbdr+1,mpl)**2+2*int(kbdr,mpl)**3
4757 IF (icost < icmn)
THEN
4769 kbdr=max(kbdr,
ibandh(k+nalc))
4770 icost=6*int(nalc-kbdr,mpl)*int(kbnd+kbdr+1,mpl)**2+2*int(kbdr,mpl)**3
4771 IF (icost < icmn)
THEN
4795 IF (
icalcm == 1.OR.lprnt) inv=2
4796 IF (mside == 1)
THEN
4804 IF (evdmin > 0.0_mpl) cndl10=log10(real(evdmax/evdmin,mps))
4814 IF ((.NOT.(
blvec(k) <= 0.0_mpd)).AND. (.NOT.(
blvec(k) > 0.0_mpd))) nan=nan+1
4819 WRITE(1,*)
'Parameter determination:',nalc,
' parameters,',
' rank=',nrank
4820 WRITE(1,*)
'-----------------------'
4821 IF(ndown /= 0)
WRITE(1,*)
' ',ndown,
' data down-weighted'
4837 wght =1.0_mpd/rerr**2
4861 IF (0.999999_mpd/wght > dvar)
THEN
4862 pull=rmeas/sqrt(1.0_mpd/wght-dvar)
4865 CALL hmpent(13,real(pull,mps))
4866 CALL gmpms(5,rec,real(pull,mps))
4868 CALL hmpent(14,real(pull,mps))
4887 IF(iter == 1.AND.jb < ist.AND.lhist) &
4888 CALL gmpms(4,rec,real(rmeas/rerr,mps))
4890 dchi2=wght*rmeas*rmeas
4895 IF(abs(resid) > chuber*rerr)
THEN
4896 wght=wght*chuber*rerr/abs(resid)
4897 dchi2=2.0*chuber*(abs(resid)/rerr-0.5*chuber)
4900 wght=wght/(1.0_mpd+(resid/rerr/cauchy)**2)
4901 dchi2=log(1.0_mpd+(resid/rerr/cauchy)**2)*cauchy**2
4911 IF(abs(resid) > chuber*rerr) chast=
'* '
4912 IF(abs(resid) > 3.0*rerr) chast=
'** '
4913 IF(abs(resid) > 6.0*rerr) chast=
'***'
4914 IF(imeas == 0)
WRITE(1,*)
'Third loop: single residuals'
4915 IF(imeas == 0)
WRITE(1,105)
4919 IF(resid < 0.0) r1=-r1
4920 IF(resid < 0.0) r2=-r2
4923105
FORMAT(
' index corrvalue residuum sigma', &
4925106
FORMAT(i6,2x,2g12.4,
' +-',g12.4,f7.2,1x,a3,f8.2)
4927 IF(iter == nter)
THEN
4929 resmax=max(resmax,abs(rmeas)/rerr)
4932 IF(iter == 1.AND.lhist)
THEN
4934 CALL hmpent( 3,real(rmeas/rerr,mps))
4936 CALL hmpent(12,real(rmeas/rerr,mps))
4943 resing=(real(nweig,mps)-real(suwt,mps))/real(nweig,mps)
4945 IF(iter == 1)
CALL hmpent( 5,real(ndf,mps))
4946 IF(iter == 1)
CALL hmpent(11,real(nalc,mps))
4951 WRITE(1,*)
'Chi^2=',summ,
' at',ndf,
' degrees of freedom: ', &
4952 '3-sigma limit is',chindl(3,ndf)*real(ndf,mps)
4953 WRITE(1,*) suwt,
' is sum of factors, compared to',nweig, &
4954 ' Downweight fraction:',resing
4960 WRITE(1,*)
' NaNs ', nalc, nrank, nan
4961 WRITE(1,*)
' ---> rejected!'
4963 IF (
mdebug < 0.AND.
nloopn == 1) print *,
' bad local fit-1 ', kfl,jrc,nrc,mside,neq,nalc,nrank,nan,cndl10
4966 IF(nrank /= nalc)
THEN
4970 WRITE(1,*)
' rank deficit', nalc, nrank
4971 WRITE(1,*)
' ---> rejected!'
4973 IF (
mdebug < 0.AND.
nloopn == 1) print *,
' bad local fit-2 ', kfl,jrc,nrc,mside,neq,nalc,nrank,nan,cndl10
4980 WRITE(1,*)
' too large condition(band part) ', nalc, nrank, cndl10
4981 WRITE(1,*)
' ---> rejected!'
4983 IF (
mdebug < 0.AND.
nloopn == 1) print *,
' bad local fit-3 ', kfl,jrc,nrc,mside,neq,nalc,nrank,nan,cndl10
4989 WRITE(1,*)
' Ndf<=0', nalc, nrank, ndf
4990 WRITE(1,*)
' ---> rejected!'
4992 IF (
mdebug < 0.AND.
nloopn == 1) print *,
' bad local fit-4 ', kfl,jrc,nrc,mside,neq,nalc,nrank,nan,cndl10
4996 chndf=real(summ/real(ndf,mpd),mps)
4998 IF(iter == 1.AND.lhist)
CALL hmpent(4,chndf)
5011 chichi=chindl(3,ndf)*real(ndf,mps)
5015 IF(summ >
chhuge*chichi)
THEN
5018 WRITE(1,*)
' ---> rejected!'
5026 IF(summ > chlimt)
THEN
5028 WRITE(1,*)
' ---> rejected!'
5032 CALL addsums(iproc+1, dchi2, ndf, dw1)
5042 CALL addsums(iproc+1, dchi2, ndf, dw1)
5046 WRITE(1,*)
' ---> rejected!'
5059 naccf(kfl)=naccf(kfl)+1
5060 ndff(kfl) =ndff(kfl) +ndf
5061 chi2f(kfl)=chi2f(kfl)+chndf
5074 wght =1.0_mpd/rerr**2
5075 dchi2=wght*rmeas*rmeas
5079 IF(resid > chuber*rerr)
THEN
5080 wght=wght*chuber*rerr/resid
5081 dchi2=2.0*chuber*(resid/rerr-0.5*chuber)
5131 CALL addsums(iproc+1, summ, ndf, dw1)
5187 IF (nfred < 0.OR.nfrei < 0)
THEN
5214 iprcnx=ijprec(i,jnx)
5216 IF (.NOT.((jnx == j+1).AND.(iprc == iprcnx)))
THEN
5230 joffd=joffd+(il*il-il)/2
5242 WRITE(1,*)
'------------------ End of printout for record',nrc
5299 iprcnx=ijprec(i,jnx)
5301 IF (.NOT.((jnx == j+1).AND.(iprc == iprcnx)))
THEN
5316 joffd=joffd+(il*il-il)/2
5352 INTEGER(mpi),
INTENT(IN) :: lun
5354 IF (
nrejec(1)>0)
WRITE(lun,*)
nrejec(1),
' (local solution contains NaNs)'
5355 IF (
nrejec(2)>0)
WRITE(lun,*)
nrejec(2),
' (local matrix with rank deficit)'
5356 IF (
nrejec(3)>0)
WRITE(lun,*)
nrejec(3),
' (local matrix with ill condition)'
5357 IF (
nrejec(4)>0)
WRITE(lun,*)
nrejec(4),
' (local fit with Ndf=0)'
5358 IF (
nrejec(5)>0)
WRITE(lun,*)
nrejec(5),
' (local fit with huge Chi2(Ndf))'
5359 IF (
nrejec(6)>0)
WRITE(lun,*)
nrejec(6),
' (local fit with large Chi2(Ndf))'
5385 INTEGER(mpi) :: icom
5386 INTEGER(mpl) :: icount
5390 INTEGER(mpi) :: imin
5391 INTEGER(mpi) :: iprlim
5392 INTEGER(mpi) :: isub
5393 INTEGER(mpi) :: itgbi
5394 INTEGER(mpi) :: itgbl
5395 INTEGER(mpi) :: ivgbi
5397 INTEGER(mpi) :: label
5405 INTEGER(mpi) :: labele(3)
5406 REAL(mps):: compnt(3)
5411 CALL mvopen(lup,
'millepede.res')
5414 WRITE(*,*)
' Result of fit for global parameters'
5415 WRITE(*,*)
' ==================================='
5420 WRITE(lup,*)
'Parameter ! first 3 elements per line are', &
5421 ' significant (if used as input)'
5439 err=sqrt(abs(real(gmati,mps)))
5440 IF(gmati < 0.0_mpd) err=-err
5443 IF(gmati*diag > 0.0_mpd)
THEN
5444 gcor2=1.0_mpd-1.0_mpd/(gmati*diag)
5445 IF(gcor2 >= 0.0_mpd.AND.gcor2 <= 1.0_mpd) gcor=real(sqrt(gcor2),mps)
5450 IF(lowstat) icount=-(icount+1)
5452 IF(itgbi <= iprlim)
THEN
5466 ELSE IF(itgbi == iprlim+1)
THEN
5467 WRITE(* ,*)
'... (further printout suppressed, but see log file)'
5481 ELSE IF (
igcorr /= 0)
THEN
5499 CALL mvopen(lup,
'millepede.eve')
5509 DO isub=0,min(15,imin-1)
5531 WRITE(lup,103) (labele(ie),compnt(ie),ie=1,iev)
5535 IF(iev /= 0)
WRITE(lup,103) (labele(ie),compnt(ie),ie=1,iev)
5542101
FORMAT(1x,
' label parameter presigma differ', &
5543 ' error'/ 1x,
'-----------',4x,4(
'-------------'))
5544102
FORMAT(i10,2x,4g14.5,f8.3)
5545103
FORMAT(3(i11,f11.7,2x))
5546110
FORMAT(i10,2x,2g14.5,28x,i12)
5547111
FORMAT(i10,2x,3g14.5,14x,i12)
5548112
FORMAT(i10,2x,4g14.5,i12)
5570 INTEGER(mpi) :: icom
5571 INTEGER(mpl) :: icount
5572 INTEGER(mpi) :: ifrst
5573 INTEGER(mpi) :: ilast
5574 INTEGER(mpi) :: inext
5575 INTEGER(mpi) :: itgbi
5576 INTEGER(mpi) :: itgbl
5577 INTEGER(mpi) :: itpgrp
5578 INTEGER(mpi) :: ivgbi
5580 INTEGER(mpi) :: icgrp
5581 INTEGER(mpi) :: ipgrp
5583 INTEGER(mpi) :: jpgrp
5585 INTEGER(mpi) :: label1
5586 INTEGER(mpi) :: label2
5587 INTEGER(mpi) :: ncon
5588 INTEGER(mpi) :: npair
5589 INTEGER(mpi) :: nstep
5592 INTEGER(mpl):: length
5594 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: vecPairedParGroups
5597 SUBROUTINE ggbmap(ipgrp,npair,npgrp)
5599 INTEGER(mpi),
INTENT(IN) :: ipgrp
5600 INTEGER(mpi),
INTENT(OUT) :: npair
5601 INTEGER(mpi),
DIMENSION(:),
INTENT(OUT) :: npgrp
5609 CALL mvopen(lup,
'millepede.res')
5610 WRITE(lup,*)
'*** Results of checking input only, no solution performed ***'
5611 WRITE(lup,*)
'! === global parameters ==='
5612 WRITE(lup,*)
'! fixed-1: by pre-sigma, -2: by entries cut, -3: by iterated entries cut'
5614 WRITE(lup,*)
'! Label Value Pre-sigma SkippedEntries Cons. group Status '
5616 WRITE(lup,*)
'! Label Value Pre-sigma Entries Cons. group Status '
5632 IF (ivgbi <= 0)
THEN
5634 IF (ivgbi == -4)
THEN
5635 WRITE(lup,116) c1,itgbl,par,presig,icount,icgrp
5637 WRITE(lup,110) c1,itgbl,par,presig,icount,icgrp,ivgbi
5641 WRITE(lup,111) c1,itgbl,par,presig,icount,icgrp
5647 WRITE(lup,*)
'!.Appearance statistics '
5648 WRITE(lup,*)
'!. Label First file and record Last file and record #files #paired-par'
5651 IF (itpgrp > 0)
THEN
5659 WRITE(lup,*)
'* === constraint groups ==='
5661 WRITE(lup,*)
'* Group #Cons. Entries First label Last label'
5663 WRITE(lup,*)
'* Group #Cons. Entries First label Last label Paired label range'
5665 CALL mpalloc(vecpairedpargroups,length,
'paired global parameter groups (I)')
5677 IF (
icheck > 1 .AND. label1 > 0)
THEN
5681 vecpairedpargroups(npair+1)=0
5685 jpgrp=vecpairedpargroups(j)
5689 IF (ifrst /= 0.AND.inext /= (ilast+nstep))
THEN
5692 WRITE(lup,114) label1, label2
5697 IF (ifrst == 0) ifrst=inext
5704 IF (jpgrp == vecpairedpargroups(j+1)-1)
THEN
5709 IF (ifrst /= 0)
THEN
5712 WRITE(lup,114) label1, label2
5718 WRITE(lup,*)
'*.Appearance statistics '
5719 WRITE(lup,*)
'*. Group First file and record Last file and record #files'
5729110
FORMAT(
' !',a1,i10,2x,2g14.5,2i12,
' fixed',i2)
5730111
FORMAT(
' !',a1,i10,2x,2g14.5,2i12,
' variable')
5731112
FORMAT(
' !.',i10,6i11)
5732113
FORMAT(
' * ',i6,i8,3i12)
5733114
FORMAT(
' *:',48x,i12,
' ..',i12)
5734115
FORMAT(
' *.',i10,5i11)
5735116
FORMAT(
' !',a1,i10,2x,2g14.5,2i12,
' redundant')
5765 INTEGER(mpi) :: iproc
5775 INTEGER(mpi),
INTENT(IN) :: n
5776 INTEGER(mpl),
INTENT(IN) :: l
5777 REAL(mpd),
INTENT(IN) :: x(n)
5778 INTEGER(mpi),
INTENT(IN) :: is
5779 INTEGER(mpi),
INTENT(IN) :: ie
5780 REAL(mpd),
INTENT(OUT) :: b(n)
5785 INTEGER(mpl) :: indij
5786 INTEGER(mpl) :: indid
5788 INTEGER(mpi) :: ichunk
5824 CALL peend(24,
'Aborted, vector/matrix size mismatch')
5825 stop
'AVPRDS: mismatched vector and matrix'
5845 IF (ia2 <= ib2) b(ia2:ib2)=b(ia2:ib2)+
globalmatd(ia2:ib2)*x(ia2:ib2)
5856 IF (i <= ie.AND.i >= is)
THEN
5864 IF (j <= ie.AND.j >= is)
THEN
5880 IF (ja2 <= jb2)
THEN
5883 b(i)=b(i)+dot_product(
globalmatd(indij+lj+ja2-ja:indij+lj+jb2-ja),x(ja2:jb2))
5887 IF (
mextnd == 0.AND.ia2 <= ib2)
THEN
5890 b(j)=b(j)+dot_product(
globalmatd(indij+lj+jn*(ia2-ia):indij+lj+jn*(ib2-ia):jn),x(ia2:ib2))
5909 IF (i <= ie.AND.i >= is)
THEN
5917 IF (j <= ie.AND.j >= is)
THEN
5933 IF (ja2 <= jb2)
THEN
5936 b(i)=b(i)+dot_product(real(
globalmatf(indij+lj+ja2-ja:indij+lj+jb2-ja),mpd),x(ja2:jb2))
5940 IF (
mextnd == 0.AND.ia2 <= ib2)
THEN
5943 b(j)=b(j)+dot_product(real(
globalmatf(indij+lj+jn*(ia2-ia):indij+lj+jn*(ib2-ia):jn),mpd),x(ia2:ib2))
5977 INTEGER(mpi) :: iproc
5985 INTEGER(mpi),
INTENT(IN) :: n
5986 INTEGER(mpl),
INTENT(IN) :: l
5987 REAL(mpd),
INTENT(IN) :: x(n)
5988 REAL(mpd),
INTENT(OUT) :: b(n)
5993 INTEGER(mpl) :: indij
5994 INTEGER(mpl) :: indid
5996 INTEGER(mpi) :: ichunk
6025 CALL peend(24,
'Aborted, vector/matrix size mismatch')
6026 stop
'AVPRD0: mismatched vector and matrix'
6072 b(i)=b(i)+dot_product(
globalmatd(indij+lj:indij+lj+jn-1),x(ja:jb))
6078 b(j)=b(j)+dot_product(
globalmatd(indij+lj:indij+jn*in:jn),x(ia:ib))
6114 b(i)=b(i)+dot_product(real(
globalmatf(indij+lj:indij+lj+jn-1),mpd),x(ja:jb))
6120 b(j)=b(j)+dot_product(real(
globalmatf(indij+lj:indij+jn*in:jn),mpd),x(ia:ib))
6144 INTEGER(mpi) :: ispc
6155 INTEGER(mpl) :: indid
6156 INTEGER(mpl),
DIMENSION(12) :: icount
6163 icount(4)=huge(icount(4))
6164 icount(7)=huge(icount(7))
6165 icount(10)=huge(icount(10))
6175 ir=ipg+(ispc-1)*(
napgrp+1)
6182 icount(1)=icount(1)+in
6183 icount(2)=icount(2)+in*ku
6189 icount(3)=icount(3)+1
6190 icount(4)=min(icount(4),jn)
6191 icount(5)=icount(5)+jn
6192 icount(6)=max(icount(6),jn)
6193 icount(7)=min(icount(7),in)
6194 icount(8)=icount(8)+in
6195 icount(9)=max(icount(9),in)
6196 icount(10)=min(icount(10),in*jn)
6197 icount(11)=icount(11)+in*jn
6198 icount(12)=max(icount(12),in*jn)
6204 WRITE(*,*)
"analysis of sparsity structure"
6205 IF (icount(1) > 0)
THEN
6206 WRITE(*,101)
"rows without compression/blocks ", icount(1)
6207 WRITE(*,101)
" contained elements ", icount(2)
6209 WRITE(*,101)
"number of block matrices ", icount(3)
6210 avg=real(icount(5),mps)/real(icount(3),mps)
6211 WRITE(*,101)
"number of columns (min,mean,max) ", icount(4), avg, icount(6)
6212 avg=real(icount(8),mps)/real(icount(3),mps)
6213 WRITE(*,101)
"number of rows (min,mean,max) ", icount(7), avg, icount(9)
6214 avg=real(icount(11),mps)/real(icount(3),mps)
6215 WRITE(*,101)
"number of elements (min,mean,max) ", icount(10), avg, icount(12)
6216101
FORMAT(2x,a34,i10,f10.3,i10)
6235 INTEGER(mpi),
INTENT(IN) :: n
6236 REAL(mpd),
INTENT(IN) :: x(n)
6237 REAL(mpd),
INTENT(OUT) :: b(n)
6242 CALL peend(24,
'Aborted, vector/matrix size mismatch')
6243 stop
'AVPROD: mismatched vector and matrix'
6274 INTEGER(mpi) :: ispc
6275 INTEGER(mpi) :: item1
6276 INTEGER(mpi) :: item2
6277 INTEGER(mpi) :: itemc
6278 INTEGER(mpi) :: jtem
6279 INTEGER(mpi) :: jtemn
6282 INTEGER(mpi),
INTENT(IN) :: itema
6283 INTEGER(mpi),
INTENT(IN) :: itemb
6284 INTEGER(mpl),
INTENT(OUT) :: ij
6285 INTEGER(mpi),
INTENT(OUT) :: lr
6286 INTEGER(mpi),
INTENT(OUT) :: iprc
6297 item1=max(itema,itemb)
6298 item2=min(itema,itemb)
6299 IF(item2 <= 0.OR.item1 >
napgrp)
RETURN
6302 outer:
DO ispc=1,
nspc
6313 IF(ku < kl) cycle outer
6318 IF(item2 >= jtem.AND.item2 < jtemn)
THEN
6324 IF(item2 < jtem)
THEN
6326 ELSE IF(item2 >= jtemn)
THEN
6341 IF(ku < kl) cycle outer
6345 IF(itemc == jtem)
EXIT
6346 IF(itemc < jtem)
THEN
6348 ELSE IF(itemc > jtem)
THEN
6376 INTEGER(mpi),
INTENT(IN) :: itema
6377 INTEGER(mpi),
INTENT(IN) :: itemb
6402 INTEGER(mpi) :: item1
6403 INTEGER(mpi) :: item2
6404 INTEGER(mpi) :: ipg1
6405 INTEGER(mpi) :: ipg2
6407 INTEGER(mpi) :: iprc
6409 INTEGER(mpi),
INTENT(IN) :: itema
6410 INTEGER(mpi),
INTENT(IN) :: itemb
6412 INTEGER(mpl) ::
ijadd
6416 item1=max(itema,itemb)
6417 item2=min(itema,itemb)
6419 IF(item2 <= 0.OR.item1 >
nagb)
RETURN
6420 IF(item1 == item2)
THEN
6429 CALL ijpgrp(ipg1,ipg2,ij,lr,iprc)
6451 INTEGER(mpi) :: item1
6452 INTEGER(mpi) :: item2
6453 INTEGER(mpi) :: jtem
6455 INTEGER(mpi),
INTENT(IN) :: itema
6456 INTEGER(mpi),
INTENT(IN) :: itemb
6465 item1=max(itema,itemb)
6466 item2=min(itema,itemb)
6468 IF(item2 <= 0.OR.item1 >
nagb)
RETURN
6476 print *,
' IJCSR3 empty list ', item1, item2, ks, ke
6477 CALL peend(23,
'Aborted, bad matrix index')
6478 stop
'ijcsr3: empty list'
6483 IF(item1 == jtem)
EXIT
6484 IF(item1 < jtem)
THEN
6491 print *,
' IJCSR3 not found ', item1, item2, ks, ke
6492 CALL peend(23,
'Aborted, bad matrix index')
6493 stop
'ijcsr3: not found'
6509 INTEGER(mpi) :: item1
6510 INTEGER(mpi) :: item2
6514 INTEGER(mpl) ::
ijadd
6517 INTEGER(mpi),
INTENT(IN) :: itema
6518 INTEGER(mpi),
INTENT(IN) :: itemb
6523 item1=max(itema,itemb)
6524 item2=min(itema,itemb)
6525 IF(item2 <= 0.OR.item1 >
nagb)
RETURN
6534 ij=
ijadd(item1,item2)
6537 ELSE IF (ij < 0)
THEN
6569 INTEGER(mpi) :: ichunk
6573 INTEGER(mpi) :: ispc
6581 INTEGER(mpl) :: ijadd
6603 ir=ipg+(ispc-1)*(
napgrp+1)
6653 REAL(mps),
INTENT(IN) :: deltat
6654 INTEGER(mpi),
INTENT(OUT) :: minut
6655 INTEGER(mpi),
INTENT(OUT):: nhour
6656 REAL(mps),
INTENT(OUT):: secnd
6657 INTEGER(mpi) :: nsecd
6660 nsecd=nint(deltat,
mpi)
6662 minut=nsecd/60-60*nhour
6663 secnd=deltat-60*(minut+60*nhour)
6699 INTEGER(mpi),
INTENT(IN) :: item
6703 INTEGER(mpl) :: length
6704 INTEGER(mpl),
PARAMETER :: four = 4
6708 IF(item <= 0)
RETURN
6729 IF(j == 0)
EXIT inner
6751 IF (
lvllog > 1)
WRITE(
lunlog,*)
'INONE: array increased to', &
6772 INTEGER(mpi) :: iprime
6773 INTEGER(mpi) :: nused
6774 LOGICAL :: finalUpdate
6775 INTEGER(mpl) :: oldLength
6776 INTEGER(mpl) :: newLength
6777 INTEGER(mpl),
PARAMETER :: four = 4
6778 INTEGER(mpi),
DIMENSION(:,:),
ALLOCATABLE :: tempArr
6779 INTEGER(mpl),
DIMENSION(:),
ALLOCATABLE :: tempVec
6783 IF(finalupdate)
THEN
6792 CALL mpalloc(temparr,four,oldlength,
'INONE: temp array')
6794 CALL mpalloc(tempvec,oldlength,
'INONE: temp vector')
6818 IF(j == 0)
EXIT inner
6819 IF(j == i) cycle outer
6823 IF(.NOT.finalupdate)
RETURN
6828 WRITE(
lunlog,*)
'INONE: array reduced to',newlength,
' words'
6852 IF(j == 0)
EXIT inner
6853 IF(j == i) cycle outer
6870 INTEGER(mpi),
INTENT(IN) :: n
6871 INTEGER(mpi) :: nprime
6872 INTEGER(mpi) :: nsqrt
6877 IF(mod(nprime,2) == 0) nprime=nprime+1
6880 nsqrt=int(sqrt(real(nprime,
mps)),
mpi)
6882 IF(i*(nprime/i) == nprime) cycle outer
6904 INTEGER(mpi) :: idum
6906 INTEGER(mpi) :: indab
6907 INTEGER(mpi) :: itgbi
6908 INTEGER(mpi) :: itgbl
6909 INTEGER(mpi) :: ivgbi
6911 INTEGER(mpi) :: jgrp
6912 INTEGER(mpi) :: lgrp
6914 INTEGER(mpi) :: nc31
6916 INTEGER(mpi) :: nwrd
6917 INTEGER(mpi) :: inone
6922 INTEGER(mpl) :: length
6923 INTEGER(mpl) :: rows
6927 WRITE(
lunlog,*)
'LOOP1: starting'
6950 WRITE(
lunlog,*)
'LOOP1: reading data files'
6964 WRITE(*,*)
'Read all binary data files:'
6966 CALL hmpldf(1,
'Number of words/record in binary file')
6967 CALL hmpdef(8,0.0,60.0,
'not_stored data per record')
6997 WRITE(
lunlog,*)
'LOOP1: reading data files again'
7009 CALL peend(21,
'Aborted, no labels/parameters defined')
7010 stop
'LOOP1: no labels/parameters defined'
7015 ' is total number NTGB of labels/parameters'
7017 CALL hmpldf(2,
'Number of entries per label')
7061 WRITE(
lunlog,*)
'LOOP1:',
npresg,
' is number of pre-sigmas'
7062 WRITE(*,*)
'LOOP1:',
npresg,
' is number of pre-sigmas'
7063 IF(
npresg == 0)
WRITE(*,*)
'Warning: no pre-sigmas defined'
7085 WRITE(
lunlog,*)
'LOOP1:',
nvgb,
' is number NVGB of variable parameters'
7090 WRITE(
lunlog,*)
'LOOP1: counting records, NO iteration of entries cut !'
7096 CALL hmpdef(15,0.0,120.0,
'Number of parameters per group')
7131 ' is total number NTPGRP of label/parameter groups'
7147 WRITE(*,112)
' Default pre-sigma =',
regpre, &
7148 ' (if no individual pre-sigma defined)'
7149 WRITE(*,*)
'Pre-sigma factor is',
regula
7152 WRITE(*,*)
'No regularization will be done'
7154 WRITE(*,*)
'Regularization will be done, using factor',
regula
7158 CALL peend(22,
'Aborted, no variable global parameters')
7159 stop
'... no variable global parameters'
7166 IF(presg > 0.0)
THEN
7168 ELSE IF(presg == 0.0.AND.
regpre > 0.0)
THEN
7169 prewt=1.0/real(
regpre**2,mpd)
7185 WRITE(*,101)
'NTGB',
ntgb,
'total number of parameters'
7186 WRITE(*,101)
'NVGB',
nvgb,
'number of variable parameters'
7190 WRITE(*,101)
'Too many variable parameters for diagonalization, fallback is inversion'
7198 WRITE(*,101)
' NREC',
nrec,
'number of records'
7199 IF (
nrecd > 0)
WRITE(*,101)
' NRECD',
nrec,
'number of records containing doubles'
7200 WRITE(*,101)
' NEQN',
neqn,
'number of equations (measurements)'
7201 WRITE(*,101)
' NEGB',
negb,
'number of equations with global parameters'
7202 WRITE(*,101)
' NDGB',
ndgb,
'number of global derivatives'
7204 WRITE(*,101)
' NZGB',
nzgb,
'number of zero global der. (ignored in entry counts)'
7207 WRITE(*,101)
'MREQENF',
mreqenf,
'required number of entries (eqns in binary files)'
7209 WRITE(*,101)
'MREQENF',
mreqenf,
'required number of entries (recs in binary files)'
7212 WRITE(*,101)
'ITEREN',
iteren,
'iterate cut for parameters with less entries'
7213 WRITE(*,101)
'MREQENA',
mreqena,
'required number of entries (from accepted fits)'
7214 IF (
mreqpe > 1)
WRITE(*,101) &
7215 'MREQPE',
mreqpe,
'required number of pair entries'
7216 IF (
msngpe >= 1)
WRITE(*,101) &
7217 'MSNGPE',
msngpe,
'max pair entries single prec. storage'
7218 WRITE(*,101)
'NTGB',
ntgb,
'total number of parameters'
7219 WRITE(*,101)
'NVGB',
nvgb,
'number of variable parameters'
7222 WRITE(*,*)
'Global parameter labels:'
7229 mqi=((mqi-20)/20)*20+1
7237 WRITE(8,101)
' NREC',
nrec,
'number of records'
7238 IF (
nrecd > 0)
WRITE(8,101)
' NRECD',
nrec,
'number of records containing doubles'
7239 WRITE(8,101)
' NEQN',
neqn,
'number of equations (measurements)'
7240 WRITE(8,101)
' NEGB',
negb,
'number of equations with global parameters'
7241 WRITE(8,101)
' NDGB',
ndgb,
'number of global derivatives'
7243 WRITE(8,101)
'MREQENF',
mreqenf,
'required number of entries (eqns in binary files)'
7245 WRITE(8,101)
'MREQENF',
mreqenf,
'required number of entries (recs in binary files)'
7248 WRITE(8,101)
'ITEREN',
iteren,
'iterate cut for parameters with less entries'
7249 WRITE(8,101)
'MREQENA',
mreqena,
'required number of entries (from accepted fits)'
7251 WRITE(
lunlog,*)
'LOOP1: ending'
7255101
FORMAT(1x,a8,
' =',i14,
' = ',a)
7271 INTEGER(mpi) :: ibuf
7273 INTEGER(mpi) :: indab
7279 INTEGER(mpi) :: nc31
7281 INTEGER(mpi) :: nlow
7283 INTEGER(mpi) :: nwrd
7285 INTEGER(mpl) :: length
7286 INTEGER(mpl),
DIMENSION(:),
ALLOCATABLE :: newCounter
7291 WRITE(
lunlog,*)
'LOOP1: iterating'
7293 WRITE(*,*)
'LOOP1: iterating'
7296 CALL mpalloc(newcounter,length,
'new entries counter')
7320 CALL isjajb(nst,ist,ja,jb,jsp)
7321 IF(ja == 0.AND.jb == 0)
EXIT
7327 IF(ij == -2) nlow=nlow+1
7332 newcounter(ij)=newcounter(ij)+1
7361 WRITE(
lunlog,*)
'LOOP1:',
nvgb,
' is number NVGB of variable parameters'
7391 INTEGER(mpi) :: ibuf
7392 INTEGER(mpi) :: icblst
7393 INTEGER(mpi) :: icboff
7394 INTEGER(mpi) :: icgb
7395 INTEGER(mpi) :: icgrp
7396 INTEGER(mpi) :: icount
7397 INTEGER(mpi) :: iext
7398 INTEGER(mpi) :: ihis
7402 INTEGER(mpi) :: ioff
7403 INTEGER(mpi) :: ipoff
7404 INTEGER(mpi) :: iproc
7405 INTEGER(mpi) :: irecmm
7407 INTEGER(mpi) :: itgbi
7408 INTEGER(mpi) :: itgbij
7409 INTEGER(mpi) :: itgbik
7410 INTEGER(mpi) :: ivgbij
7411 INTEGER(mpi) :: ivgbik
7412 INTEGER(mpi) :: ivpgrp
7416 INTEGER(mpi) :: jcgrp
7417 INTEGER(mpi) :: jext
7418 INTEGER(mpi) :: jcgb
7419 INTEGER(mpi) :: jrec
7421 INTEGER(mpi) :: joff
7423 INTEGER(mpi) :: kcgrp
7424 INTEGER(mpi) :: kfile
7426 INTEGER(mpi) :: label
7427 INTEGER(mpi) :: labelf
7428 INTEGER(mpi) :: labell
7429 INTEGER(mpi) :: lvpgrp
7432 INTEGER(mpi) :: maeqnf
7433 INTEGER(mpi) :: nall
7434 INTEGER(mpi) :: naeqna
7435 INTEGER(mpi) :: naeqnf
7436 INTEGER(mpi) :: naeqng
7437 INTEGER(mpi) :: npdblk
7438 INTEGER(mpi) :: nc31
7439 INTEGER(mpi) :: ncachd
7440 INTEGER(mpi) :: ncachi
7441 INTEGER(mpi) :: ncachr
7442 INTEGER(mpi) :: ncon
7445 INTEGER(mpi) :: ndfmax
7446 INTEGER(mpi) :: nfixed
7447 INTEGER(mpi) :: nggd
7448 INTEGER(mpi) :: nggi
7449 INTEGER(mpi) :: nmatmo
7450 INTEGER(mpi) :: noff
7451 INTEGER(mpi) :: npair
7452 INTEGER(mpi) :: npar
7453 INTEGER(mpi) :: nparmx
7455 INTEGER(mpi) :: nrece
7456 INTEGER(mpi) :: nrecf
7457 INTEGER(mpi) :: nrecmm
7459 INTEGER(mpi) :: nwrd
7460 INTEGER(mpi) :: inone
7469 INTEGER(mpl):: nblock
7470 INTEGER(mpl):: nbwrds
7471 INTEGER(mpl):: noff8
7472 INTEGER(mpl):: ndimbi
7473 INTEGER(mpl):: ndimsa(4)
7475 INTEGER(mpl):: nnzero
7476 INTEGER(mpl):: matsiz(2)
7477 INTEGER(mpl):: matwords
7478 INTEGER(mpl):: mbwrds
7479 INTEGER(mpl):: length
7482 INTEGER(mpl),
PARAMETER :: two=2
7483 INTEGER(mpi) :: maxGlobalPar = 0
7484 INTEGER(mpi) :: maxLocalPar = 0
7485 INTEGER(mpi) :: maxEquations = 0
7487 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: vecConsGroupList
7488 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: vecConsGroupIndex
7489 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: vecPairedParGroups
7490 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: vecBlockCounts
7493 SUBROUTINE ndbits(npgrp,ndims,nsparr,ihst)
7495 INTEGER(mpi),
DIMENSION(:),
INTENT(IN) :: npgrp
7496 INTEGER(mpl),
DIMENSION(4),
INTENT(OUT) :: ndims
7497 INTEGER(mpl),
DIMENSION(:,:),
INTENT(OUT) :: nsparr
7498 INTEGER(mpi),
INTENT(IN) :: ihst
7500 SUBROUTINE ckbits(npgrp,ndims)
7502 INTEGER(mpi),
DIMENSION(:),
INTENT(IN) :: npgrp
7503 INTEGER(mpl),
DIMENSION(4),
INTENT(OUT) :: ndims
7505 SUBROUTINE spbits(npgrp,nsparr,nsparc)
7507 INTEGER(mpi),
DIMENSION(:),
INTENT(IN) :: npgrp
7508 INTEGER(mpl),
DIMENSION(:,:),
INTENT(IN) :: nsparr
7509 INTEGER(mpi),
DIMENSION(:),
INTENT(OUT) :: nsparc
7511 SUBROUTINE gpbmap(ngroup,npgrp,npair)
7513 INTEGER(mpi),
INTENT(IN) :: ngroup
7514 INTEGER(mpi),
DIMENSION(:,:),
INTENT(IN) :: npgrp
7515 INTEGER(mpi),
DIMENSION(:),
INTENT(OUT) :: npair
7517 SUBROUTINE ggbmap(ipgrp,npair,npgrp)
7519 INTEGER(mpi),
INTENT(IN) :: ipgrp
7520 INTEGER(mpi),
INTENT(OUT) :: npair
7521 INTEGER(mpi),
DIMENSION(:),
INTENT(OUT) :: npgrp
7523 SUBROUTINE pbsbits(npgrp,ibsize,nnzero,nblock,nbkrow)
7525 INTEGER(mpi),
DIMENSION(:),
INTENT(IN) :: npgrp
7526 INTEGER(mpi),
INTENT(IN) :: ibsize
7527 INTEGER(mpl),
INTENT(OUT) :: nnzero
7528 INTEGER(mpl),
INTENT(OUT) :: nblock
7529 INTEGER(mpi),
DIMENSION(:),
INTENT(OUT) :: nbkrow
7531 SUBROUTINE pblbits(npgrp,ibsize,nsparr,nsparc)
7533 INTEGER(mpi),
DIMENSION(:),
INTENT(IN) :: npgrp
7534 INTEGER(mpi),
INTENT(IN) :: ibsize
7535 INTEGER(mpl),
DIMENSION(:),
INTENT(IN) :: nsparr
7536 INTEGER(mpl),
DIMENSION(:),
INTENT(OUT) :: nsparc
7538 SUBROUTINE prbits(npgrp,nsparr)
7540 INTEGER(mpi),
DIMENSION(:),
INTENT(IN) :: npgrp
7541 INTEGER(mpl),
DIMENSION(:),
INTENT(OUT) :: nsparr
7543 SUBROUTINE pcbits(npgrp,nsparr,nsparc)
7545 INTEGER(mpi),
DIMENSION(:),
INTENT(IN) :: npgrp
7546 INTEGER(mpl),
DIMENSION(:),
INTENT(IN) :: nsparr
7547 INTEGER(mpl),
DIMENSION(:),
INTENT(OUT) :: nsparc
7557 WRITE(
lunlog,*)
'LOOP2: starting'
7580 WRITE(
lunlog,*)
' Elimination for constraints enforced for solution by decomposition!'
7585 WRITE(
lunlog,*)
' Lagrange multipliers enforced for solution by sparsePARDISO!'
7590 WRITE(
lunlog,*)
' Elimination for constraints with mpqldec enforced (LAPACK only for unpacked storage)!'
7607 noff8=int(
nagb,mpl)*int(
nagb-1,mpl)/2
7643 print *,
' Variable parameter groups ',
nvpgrp
7691 CALL mpalloc(vecconsgrouplist,length,
'constraint group list')
7692 CALL mpalloc(vecconsgroupindex,length,
'constraint group index')
7702 WRITE(
lunlog,*)
'LOOP2: start event reading'
7742 WRITE(*,*)
'Record number ',
nrec,
' from file ',kfile
7743 IF (wgh /= 1.0)
WRITE(*,*)
' weight ',wrec
7747 CALL isjajb(nst,ist,ja,jb,jsp)
7751 IF(nda ==
mdebg2+1)
WRITE(*,*)
'... and more data'
7756 WRITE(*,*)
'Local derivatives:'
7758107
FORMAT(6(i3,g12.4))
7760 WRITE(*,*)
'Global derivatives:'
7763108
FORMAT(3i11,g12.4)
7766 WRITE(*,*)
'total_par_label __label__ var_par_index derivative'
7781 CALL isjajb(nst,ist,ja,jb,jsp)
7782 IF(ja == 0.AND.jb == 0)
EXIT
7822 IF (kcgrp == jcgrp) cycle
7833 IF (vecconsgroupindex(k) == 0)
THEN
7836 vecconsgrouplist(icgrp)=k
7851 IF (vecconsgroupindex(k) < icount)
THEN
7853 vecconsgroupindex(k)=icount
7871 IF (nfixed > 0) naeqnf=naeqnf+1
7874 IF(ja /= 0.AND.jb /= 0)
THEN
7883 IF (naeqnf > maeqnf) nrecf=nrecf+1
7887 maxglobalpar=max(
nagbn,maxglobalpar)
7888 maxlocalpar=max(
nalcn,maxlocalpar)
7889 maxequations=max(
naeqn,maxequations)
7892 dstat(1)=dstat(1)+real((nwrd+2)*2,mpd)
7893 dstat(2)=dstat(2)+real(
nagbn+2,mpd)
7898 vecconsgroupindex(vecconsgrouplist(k))=0
7903 IF (
nagbn == 0)
THEN
7922 IF (ivpgrp /= lvpgrp)
THEN
7977 (irecmm >= nrecmm.OR.irecmm ==
mxrec))
THEN
7978 IF (nmatmo == 0)
THEN
7980 WRITE(*,*)
'Monitoring of sparse matrix construction'
7981 WRITE(*,*)
' records ........ off-diagonal elements ', &
7982 '....... compression memory'
7983 WRITE(*,*)
' non-zero used(double) used', &
7988 gbc=1.0e-9*real((mpi*ndimsa(2)+mpd*ndimsa(3)+mps*ndimsa(4))/mpi*(bit_size(1_mpi)/8),mps)
7989 gbu=1.0e-9*real(((mpi+mpd)*(ndimsa(3)+ndimsa(4)))/mpi*(bit_size(1_mpi)/8),mps)
7991 WRITE(*,1177) irecmm,ndimsa(1),ndimsa(3),ndimsa(4),cpr,gbc
79921177
FORMAT(i9,3i13,f10.2,f11.6)
7993 DO WHILE(irecmm >= nrecmm)
8012 WRITE(
lunlog,*)
'LOOP2: event reading ended - end of data'
8014 dstat(k)=dstat(k)/real(
nrec,mpd)
8028 CALL mpalloc(vecpairedpargroups,length,
'paired global parameter groups (I)')
8030 print *,
' Total parameter groups pairs',
ntpgrp
8033 CALL ggbmap(i,npair,vecpairedpargroups)
8101 IF(ivgbij > 0.AND.ivgbik > 0)
THEN
8103 IF (
mprint > 1)
WRITE(*,*)
'add index pair ',ivgbij,ivgbik
8180 nparmx=max(nparmx,int(rows,mpi))
8202 WRITE(*,*)
' Parameter block', i, ib-ia, jb-ja, labelf, labell
8206 WRITE(
lunlog,*)
'Detected',
npblck,
'(disjoint) parameter blocks, max size ', nparmx
8208 WRITE(*,*)
'Detected',
npblck,
'(disjoint) parameter blocks, max size ', nparmx
8210 WRITE(*,*)
'Using block diagonal storage mode'
8225 IF (
nagb >= 65536)
THEN
8226 noff=int(noff8/1000,mpi)
8236 CALL hmpdef(ihis,0.0,real(
mhispe,mps),
'NDBITS: #off-diagonal elements')
8241 ndgn=ndimsa(3)+ndimsa(4)
8242 matwords=ndimsa(2)+length*4
8257 length=int(npdblk,mpl)
8258 CALL mpalloc(vecblockcounts,length,
'sparse matrix row offsets (CSR3)')
8260 nbwrds=2*int(nblock,mpl)*int(
ipdbsz(i)*
ipdbsz(i)+1,mpl)
8261 IF ((i == 1).OR.(nbwrds < mbwrds))
THEN
8285 IF (
fcache(2) == 0.0)
THEN
8287 fcache(2)=real(dstat(2),mps)
8288 fcache(3)=real(dstat(3),mps)
8309 ncachd=
ncache-ncachr-ncachi
8340 WRITE(lu,101)
'NTGB',
ntgb,
'total number of parameters'
8341 WRITE(lu,102)
'(all parameters, appearing in binary files)'
8342 WRITE(lu,101)
'NVGB',
nvgb,
'number of variable parameters'
8343 WRITE(lu,102)
'(appearing in fit matrix/vectors)'
8344 WRITE(lu,101)
'NAGB',
nagb,
'number of all parameters'
8345 WRITE(lu,102)
'(including Lagrange multiplier or reduced)'
8346 WRITE(lu,101)
'NTPGRP',
ntpgrp,
'total number of parameter groups'
8347 WRITE(lu,101)
'NVPGRP',
nvpgrp,
'number of variable parameter groups'
8348 WRITE(lu,101)
'NFGB',
nfgb,
'number of fit parameters'
8350 WRITE(lu,101)
'MBANDW',
mbandw,
'band width of preconditioner matrix'
8351 WRITE(lu,102)
'(if <0, no preconditioner matrix)'
8353 IF (
nagb >= 65536)
THEN
8354 WRITE(lu,101)
'NOFF/K',noff,
'max number of off-diagonal elements'
8356 WRITE(lu,101)
'NOFF',noff,
'max number of off-diagonal elements'
8359 IF (
nagb >= 65536)
THEN
8360 WRITE(lu,101)
'NDGN/K',ndgn/1000,
'actual number of off-diagonal elements'
8362 WRITE(lu,101)
'NDGN',ndgn,
'actual number of off-diagonal elements'
8365 WRITE(lu,101)
'NCGB',
ncgb,
'number of constraints'
8366 WRITE(lu,101)
'NAGBN',
nagbn,
'max number of global parameters in an event'
8367 WRITE(lu,101)
'NALCN',
nalcn,
'max number of local parameters in an event'
8368 WRITE(lu,101)
'NAEQN',
naeqn,
'max number of equations in an event'
8370 WRITE(lu,101)
'NAEQNA',naeqna,
'number of equations'
8371 WRITE(lu,101)
'NAEQNG',naeqng, &
8372 'number of equations with global parameters'
8373 WRITE(lu,101)
'NAEQNF',naeqnf, &
8374 'number of equations with fixed global parameters'
8375 WRITE(lu,101)
'NRECF',nrecf, &
8376 'number of records with fixed global parameters'
8379 WRITE(lu,101)
'NRECE',nrece, &
8380 'number of records without variable parameters'
8383 WRITE(lu,101)
'NCACHE',
ncache,
'number of words for caching'
8384 WRITE(lu,111) (
fcache(k)*100.0,k=1,3)
8385111
FORMAT(22x,
'cache splitting ',3(f6.1,
' %'))
8390 WRITE(lu,*)
'Solution method and matrix-storage mode:'
8392 WRITE(lu,*)
' METSOL = 1: matrix inversion'
8393 ELSE IF(
metsol == 2)
THEN
8394 WRITE(lu,*)
' METSOL = 2: diagonalization'
8395 ELSE IF(
metsol == 3)
THEN
8396 WRITE(lu,*)
' METSOL = 3: decomposition'
8397 ELSE IF(
metsol == 4)
THEN
8398 WRITE(lu,*)
' METSOL = 4: MINRES (rtol',
mrestl,
')'
8399 ELSE IF(
metsol == 5)
THEN
8400 WRITE(lu,*)
' METSOL = 5: MINRES-QLP (rtol',
mrestl,
')'
8401 ELSE IF(
metsol == 6)
THEN
8402 WRITE(lu,*)
' METSOL = 6: GMRES'
8404 ELSE IF(
metsol == 7)
THEN
8405 WRITE(lu,*)
' METSOL = 7: LAPACK factorization'
8406 ELSE IF(
metsol == 8)
THEN
8407 WRITE(lu,*)
' METSOL = 8: LAPACK factorization'
8409 ELSE IF(
metsol == 9)
THEN
8410 WRITE(lu,*)
' METSOL = 9: Intel oneMKL PARDISO'
8414 WRITE(lu,*)
' with',
mitera,
' iterations'
8416 WRITE(lu,*)
' MATSTO = 0: unpacked symmetric matrix, ',
'n*n elements'
8417 ELSE IF(
matsto == 1)
THEN
8418 WRITE(lu,*)
' MATSTO = 1: full symmetric matrix, ',
'(n*n+n)/2 elements'
8419 ELSE IF(
matsto == 2)
THEN
8420 WRITE(lu,*)
' MATSTO = 2: sparse matrix (custom)'
8421 ELSE IF(
matsto == 3)
THEN
8423 WRITE(lu,*)
' MATSTO = 3: sparse matrix (upper triangle, CSR3)'
8425 WRITE(lu,*)
' MATSTO = 3: sparse matrix (upper triangle, BSR3)'
8426 WRITE(lu,*)
' block size',
matbsz
8430 WRITE(lu,*)
' block diagonal with',
npblck,
' blocks'
8432 IF(
mextnd>0)
WRITE(lu,*)
' with extended storage'
8433 IF(
dflim /= 0.0)
THEN
8434 WRITE(lu,103)
'Convergence assumed, if expected dF <',
dflim
8439 WRITE(lu,*)
'Constraints handled by elimination with LAPACK'
8441 WRITE(lu,*)
'Constraints handled by elimination'
8444 WRITE(lu,*)
'Constraints handled by Lagrange multipliers'
8451 CALL peend(28,
'Aborted, no local parameters')
8452 stop
'LOOP2: stopping due to missing local parameters'
8473105
FORMAT(
' Constants C1, C2 for Wolfe conditions:',g12.4,
', ',g12.4)
8481 matwords=(length+
nagb+1)*2
8488 ELSE IF (
matsto == 2)
THEN
8489 matsiz(1)=ndimsa(3)+
nagb
8504 IF (icblst > icboff)
THEN
8510 nall = npar;
IF (
icelim <= 0) nall=npar+ncon
8514 matsiz(1)=matsiz(1)+k
8516 matsiz(1)=matsiz(1)+nall
8521 matwords=matwords+matsiz(1)*2+matsiz(2)
8527 IF (matwords < 250000)
THEN
8528 WRITE(*,*)
'Size of global matrix: < 1 MB'
8530 WRITE(*,*)
'Size of global matrix:',int(real(matwords,mps)*4.0e-6,mpi),
' MB'
8536 WRITE(
lunlog,*)
' Cut values of Chi^2/Ndf and Chi2,'
8537 WRITE(
lunlog,*)
' corresponding to 2 and 3 standard deviations'
8538 WRITE(
lunlog,*)
' Ndf Chi^2/Ndf(2) Chi^2(2) ', &
8539 ' Chi^2/Ndf(3) Chi^2(3)'
8542 IF(ndf >
naeqn)
EXIT
8545 ELSE IF(ndf < 20)
THEN
8547 ELSE IF(ndf < 100)
THEN
8549 ELSE IF(ndf < 200)
THEN
8556 WRITE(
lunlog,106) ndf,chin2,chin2*real(ndf,mps),chin3, chin3*real(ndf,mps)
8559 WRITE(
lunlog,*)
'LOOP2: ending'
8563 IF (
ncgbe /= 0)
THEN
8566 WRITE(*,199)
'WarningWarningWarningWarningWarningWarningWarningWarningWar'
8567 WRITE(*,199)
'arningWarningWarningWarningWarningWarningWarningWarningWarn'
8568 WRITE(*,199)
'rningWarningWarningWarningWarningWarningWarningWarningWarni'
8569 WRITE(*,199)
'ningWarningWarningWarningWarningWarningWarningWarningWarnin'
8570 WRITE(*,199)
'ingWarningWarningWarningWarningWarningWarningWarningWarning'
8571 WRITE(*,199)
'ngWarningWarningWarningWarningWarningWarningWarningWarningW'
8572 WRITE(*,199)
'gWarningWarningWarningWarningWarningWarningWarningWarningWa'
8574 WRITE(*,*)
' Number of empty constraints =',abs(
ncgbe),
', should be 0'
8575 WRITE(*,*)
' => please check constraint definition, mille data'
8577 WRITE(*,199)
'WarningWarningWarningWarningWarningWarningWarningWarningWar'
8578 WRITE(*,199)
'arningWarningWarningWarningWarningWarningWarningWarningWarn'
8579 WRITE(*,199)
'rningWarningWarningWarningWarningWarningWarningWarningWarni'
8580 WRITE(*,199)
'ningWarningWarningWarningWarningWarningWarningWarningWarnin'
8581 WRITE(*,199)
'ingWarningWarningWarningWarningWarningWarningWarningWarning'
8582 WRITE(*,199)
'ngWarningWarningWarningWarningWarningWarningWarningWarningW'
8583 WRITE(*,199)
'gWarningWarningWarningWarningWarningWarningWarningWarningWa'
8588101
FORMAT(1x,a8,
' =',i14,
' = ',a)
8590103
FORMAT(1x,a,g12.4)
8591106
FORMAT(i6,2(3x,f9.3,f12.1,3x))
8606 INTEGER(mpi) :: imed
8609 INTEGER(mpi) :: nent
8610 INTEGER(mpi),
DIMENSION(measBins) :: isuml
8611 INTEGER(mpi),
DIMENSION(measBins) :: isums
8615 INTEGER(mpl) :: ioff
8618 DATA lfirst /.true./
8631 WRITE(
lunmon,
'(A)')
'*** Normalized residuals grouped by first global label (per local fit cycle) ***'
8633 WRITE(
lunmon,
'(A)')
'*** Pulls grouped by first global label (per local fit cycle) ***'
8635 WRITE(
lunmon,
'(A)')
'! LFC Label Entries Median RMS(MAD) <error>'
8654 IF (2*isuml(j) > nent)
EXIT
8658 IF (isuml(j) > isuml(j-1)) amed=amed+real(nent-2*isuml(j-1),mps)/real(2*isuml(j)-2*isuml(j-1),mps)
8671 isums(j)=isums(j)+isums(j-1)
8675 IF (2*isums(j) > nent)
EXIT
8678 IF (isums(j) > isums(j-1)) amad=amad+real(nent-2*isums(j-1),mps)/real(2*isums(j)-2*isums(j-1),mps)
8692110
FORMAT(i5,2i10,3g14.5)
8707 INTEGER(mpi) :: ioff
8708 INTEGER(mpi) :: ipar0
8709 INTEGER(mpi) :: ncon
8710 INTEGER(mpi) :: npar
8711 INTEGER(mpi) :: nextra
8713 INTEGER :: nbopt, nboptx, ILAENV
8716 INTEGER(mpl),
INTENT(IN) :: msize(2)
8718 INTEGER(mpl) :: length
8719 INTEGER(mpl) :: nwrdpc
8720 INTEGER(mpl),
PARAMETER :: three = 3
8733 CALL mpalloc(
aux,length,
' local fit scratch array: aux')
8734 CALL mpalloc(
vbnd,length,
' local fit scratch array: vbnd')
8735 CALL mpalloc(
vbdr,length,
' local fit scratch array: vbdr')
8738 CALL mpalloc(
vbk,length,
' local fit scratch array: vbk')
8741 CALL mpalloc(
vzru,length,
' local fit scratch array: vzru')
8742 CALL mpalloc(
scdiag,length,
' local fit scratch array: scdiag')
8743 CALL mpalloc(
scflag,length,
' local fit scratch array: scflag')
8744 CALL mpalloc(
ibandh,2*length,
' local fit band width hist.: ibandh')
8758 length=4*
ncblck;
IF(ncon == 0) length=0
8767 nwrdpc=nwrdpc+length
8779 length=(ncon*ncon+ncon)/2
8804 length=length+npar+nextra
8810 IF(
mbandw == 0) length=length+1
8818 nwrdpc=nwrdpc+2*length
8819 IF (nwrdpc > 250000)
THEN
8821 WRITE(*,*)
'Size of preconditioner matrix:',int(real(nwrdpc,mps)*4.0e-6,mpi),
' MB'
8843 CALL peend(23,
'Aborted, bad matrix index (will exceed 32bit)')
8844 stop
'vmprep: bad index (matrix to large for diagonalization)'
8866 nbopt = ilaenv( 1_mpl,
'DSYTRF',
'U', int(
nagb,mpl), int(
nagb,mpl), -1_mpl, -1_mpl )
8868 print *,
'LAPACK optimal block size for DSYTRF:', nbopt
8869 lplwrk=length*int(nbopt,mpl)
8877 nbopt = ilaenv( 1_mpl,
'DORMQL',
'RN', int(npar,mpl), int(npar,mpl), int(ncon,mpl), int(npar,mpl) )
8878 IF (int(npar,mpl)*int(nbopt,mpl) >
lplwrk)
THEN
8879 lplwrk=int(npar,mpl)*int(nbopt,mpl)
8884 print *,
'LAPACK optimal block size for DORMQL:', nboptx
8903 INTEGER(mpi) :: icoff
8904 INTEGER(mpi) :: ipoff
8907 INTEGER(mpi) :: ncon
8908 INTEGER(mpi) :: nfit
8909 INTEGER(mpi) :: npar
8910 INTEGER(mpi) :: nrank
8911 INTEGER(mpl) :: imoff
8912 INTEGER(mpl) :: ioff1
8930 WRITE(
lunlog,*)
'Shrinkage of global matrix (A->Q^t*A*Q)'
8945 nfit=npar+ncon;
IF (
icelim > 0) nfit=npar-ncon
8947 IF(nfit < npar)
THEN
8965 WRITE(
lunlog,*)
'Inversion of global matrix (A->A^-1)'
8972 IF(nfit /= nrank)
THEN
8973 WRITE(*,*)
'Warning: the rank defect of the symmetric',nfit, &
8974 '-by-',nfit,
' matrix is ',nfit-nrank,
' (should be zero).'
8975 WRITE(lun,*)
'Warning: the rank defect of the symmetric',nfit, &
8976 '-by-',nfit,
' matrix is ',nfit-nrank,
' (should be zero).'
8979 WRITE(*,*)
' --> enforcing SUBITO mode'
8980 WRITE(lun,*)
' --> enforcing SUBITO mode'
8982 ELSE IF(
ndefec == 0)
THEN
8984 WRITE(lun,*)
'No rank defect of the symmetric matrix'
8986 WRITE(lun,*)
'No rank defect of the symmetric block', ib,
' of size', npar
8997 IF(nfit < npar)
THEN
9016 INTEGER(mpi) :: icoff
9017 INTEGER(mpi) :: ipoff
9020 INTEGER(mpi) :: ncon
9021 INTEGER(mpi) :: nfit
9022 INTEGER(mpi) :: npar
9023 INTEGER(mpi) :: nrank
9024 INTEGER(mpl) :: imoff
9025 INTEGER(mpl) :: ioff1
9040 WRITE(
lunlog,*)
'Shrinkage of global matrix (A->Q^t*A*Q)'
9054 nfit=npar+ncon;
IF (
icelim > 0) nfit=npar-ncon
9056 IF(nfit < npar)
THEN
9074 WRITE(
lunlog,*)
'Decomposition of global matrix (A->L*D*L^t)'
9080 IF(nfit /= nrank)
THEN
9081 WRITE(*,*)
'Warning: the rank defect of the symmetric',nfit, &
9082 '-by-',nfit,
' matrix is ',nfit-nrank,
' (should be zero).'
9083 WRITE(lun,*)
'Warning: the rank defect of the symmetric',nfit, &
9084 '-by-',nfit,
' matrix is ',nfit-nrank,
' (should be zero).'
9087 WRITE(*,*)
' --> enforcing SUBITO mode'
9088 WRITE(lun,*)
' --> enforcing SUBITO mode'
9090 ELSE IF(
ndefec == 0)
THEN
9092 WRITE(lun,*)
'No rank defect of the symmetric matrix'
9094 WRITE(lun,*)
'No rank defect of the symmetric block', ib,
' of size', npar
9096 WRITE(lun,*)
' largest diagonal element (LDLt)', evmax
9097 WRITE(lun,*)
' smallest diagonal element (LDLt)', evmin
9106 IF(nfit < npar)
THEN
9128 INTEGER(mpi) :: icoff
9129 INTEGER(mpi) :: ipoff
9132 INTEGER(mpi) :: ncon
9133 INTEGER(mpi) :: nfit
9134 INTEGER(mpi) :: npar
9135 INTEGER(mpl) :: imoff
9136 INTEGER(mpl) :: ioff1
9137 INTEGER(mpi) :: infolp
9157 WRITE(
lunlog,*)
'Shrinkage of global matrix (A->Q^t*A*Q)'
9172 nfit=npar+ncon;
IF (
icelim > 0) nfit=npar-ncon
9174 IF(nfit < npar)
THEN
9191 IF (nfit > npar)
THEN
9194 WRITE(
lunlog,*)
'Factorization of global matrix (A->L*D*L^t)'
9204 WRITE(
lunlog,*)
'Factorization of global matrix (A->L*L^t)'
9208 CALL dpptrf(
'U',int(nfit,mpl),
globalmatd(imoff+1:),infolp)
9215 WRITE(lun,*)
'No rank defect of the symmetric matrix'
9217 WRITE(lun,*)
'No rank defect of the symmetric block', ib,
' of size', npar
9221 WRITE(*,*)
'Warning: factorization of the symmetric',nfit, &
9222 '-by-',nfit,
' failed at index ', infolp
9223 WRITE(lun,*)
'Warning: factorization of the symmetric',nfit, &
9224 '-by-',nfit,
' failed at index ', infolp
9225 CALL peend(29,
'Aborted, factorization of global matrix failed')
9226 stop
'mdptrf: bad matrix'
9231 IF (nfit > npar)
THEN
9234 IF(infolp /= 0) print *,
' DSPTRS failed: ', infolp
9236 CALL dpptrs(
'U',int(nfit,mpl),1_mpl,
globalmatd(imoff+1:),&
9238 IF(infolp /= 0) print *,
' DPPTRS failed: ', infolp
9242 IF(nfit < npar)
THEN
9263 INTEGER(mpi) :: icoff
9264 INTEGER(mpi) :: ipoff
9267 INTEGER(mpi) :: ncon
9268 INTEGER(mpi) :: nfit
9269 INTEGER(mpi) :: npar
9270 INTEGER(mpl) :: imoff
9271 INTEGER(mpl) :: ioff1
9272 INTEGER(mpl) :: iloff
9273 INTEGER(mpi) :: infolp
9294 WRITE(
lunlog,*)
'Shrinkage of global matrix (A->Q^t*A*Q)'
9314 nfit=npar+ncon;
IF (
icelim > 0) nfit=npar-ncon
9316 IF(nfit < npar)
THEN
9320 CALL dtrtrs(
'L',
'T',
'N',int(ncon,mpl),1_mpl,
lapackql(iloff+npar-ncon+1:),int(npar,mpl),&
9322 IF(infolp /= 0) print *,
' DTRTRS failed: ', infolp
9324 CALL dormql(
'L',
'T',int(npar,mpl),1_mpl,int(ncon,mpl),
lapackql(iloff+1:),int(npar,mpl),&
9326 IF(infolp /= 0) print *,
' DORMQL failed: ', infolp
9345 IF (nfit > npar)
THEN
9348 WRITE(
lunlog,*)
'Factorization of global matrix (A->L*D*L^t)'
9352 CALL dsytrf(
'U',int(nfit,mpl),
globalmatd(imoff+1:),int(nfit,mpl),&
9359 WRITE(
lunlog,*)
'Factorization of global matrix (A->L*L^t)'
9363 CALL dpotrf(
'U',int(nfit,mpl),
globalmatd(imoff+1:),int(npar,mpl),infolp)
9370 WRITE(lun,*)
'No rank defect of the symmetric matrix'
9372 WRITE(lun,*)
'No rank defect of the symmetric block', ib,
' of size', npar
9376 WRITE(*,*)
'Warning: factorization of the symmetric',nfit, &
9377 '-by-',nfit,
' failed at index ', infolp
9378 WRITE(lun,*)
'Warning: factorization of the symmetric',nfit, &
9379 '-by-',nfit,
' failed at index ', infolp
9380 CALL peend(29,
'Aborted, factorization of global matrix failed')
9381 stop
'mdutrf: bad matrix'
9386 IF (nfit > npar)
THEN
9387 CALL dsytrs(
'U',int(nfit,mpl),1_mpl,
globalmatd(imoff+1:),int(nfit,mpl),&
9389 IF(infolp /= 0) print *,
' DSYTRS failed: ', infolp
9391 CALL dpotrs(
'U',int(nfit,mpl),1_mpl,
globalmatd(imoff+1:),int(npar,mpl),&
9393 IF(infolp /= 0) print *,
' DPOTRS failed: ', infolp
9397 IF(nfit < npar)
THEN
9402 CALL dormql(
'L',
'N',int(npar,mpl),1_mpl,int(ncon,mpl),
lapackql(iloff+1:),int(npar,mpl),&
9404 IF(infolp /= 0) print *,
' DORMQL failed: ', infolp
9411 iloff=iloff+int(npar,mpl)*int(ncon,mpl)
9433 INTEGER(mpi) :: icboff
9434 INTEGER(mpi) :: icblst
9435 INTEGER(mpi) :: icoff
9436 INTEGER(mpi) :: icfrst
9437 INTEGER(mpi) :: iclast
9438 INTEGER(mpi) :: ipfrst
9439 INTEGER(mpi) :: iplast
9440 INTEGER(mpi) :: ipoff
9443 INTEGER(mpi) :: ncon
9444 INTEGER(mpi) :: npar
9446 INTEGER(mpl) :: imoff
9447 INTEGER(mpl) :: iloff
9448 INTEGER(mpi) :: infolp
9449 INTEGER :: nbopt, ILAENV
9451 REAL(mpd),
INTENT(IN) :: a(mszcon)
9452 REAL(mpd),
INTENT(OUT) :: emin
9453 REAL(mpd),
INTENT(OUT) :: emax
9464 iloff=iloff+int(npar,mpl)*int(ncon,mpl)
9483 DO icb=icboff+1,icboff+icblst
9490 lapackql(iloff+ipfrst:iloff+iplast)=a(imoff+1:imoff+npb)
9507 nbopt = ilaenv( 1_mpl,
'DGEQLF',
'', int(npar,mpl), int(ncon,mpl), int(npar,mpl), -1_mpl )
9508 print *,
'LAPACK optimal block size for DGEQLF:', nbopt
9509 lplwrk=int(ncon,mpl)*int(nbopt,mpl)
9512 CALL dgeqlf(int(npar,mpl),int(ncon,mpl),
lapackql(iloff+1:),int(npar,mpl),&
9514 IF(infolp /= 0) print *,
' DGEQLF failed: ', infolp
9517 iloff=iloff+int(npar,mpl)*int(ncon,mpl)
9520 IF(emax < emin)
THEN
9548 INTEGER(mpi) :: icoff
9549 INTEGER(mpi) :: ipoff
9551 INTEGER(mpi) :: ncon
9552 INTEGER(mpi) :: npar
9553 INTEGER(mpl) :: imoff
9554 INTEGER(mpl) :: iloff
9555 INTEGER(mpi) :: infolp
9556 CHARACTER (LEN=1) :: transr, transl
9558 LOGICAL,
INTENT(IN) :: t
9577 IF(ncon <= 0 ) cycle
9586 DO i=ipoff+1,ipoff+npar
9592 CALL dormql(
'R',transr,int(npar,mpl),int(npar,mpl),int(ncon,mpl),
lapackql(iloff+1:),&
9595 IF(infolp /= 0) print *,
' DORMQL failed: ', infolp
9597 CALL dormql(
'L',transl,int(npar,mpl),int(npar,mpl),int(ncon,mpl),
lapackql(iloff+1:),&
9600 IF(infolp /= 0) print *,
' DORMQL failed: ', infolp
9603 iloff=iloff+int(npar,mpl)*int(ncon,mpl)
9609include
'mkl_pardiso.f90'
9640 TYPE(mkl_pardiso_handle) :: pt(64)
9642 INTEGER(mpl),
PARAMETER :: maxfct =1
9643 INTEGER(mpl),
PARAMETER :: mnum = 1
9644 INTEGER(mpl),
PARAMETER :: nrhs = 1
9646 INTEGER(mpl) :: mtype
9647 INTEGER(mpl) :: phase
9648 INTEGER(mpl) :: error
9649 INTEGER(mpl) :: msglvl
9653 INTEGER(mpl) :: idum(1)
9655 INTEGER(mpl) :: length
9656 INTEGER(mpi) :: nfill
9657 INTEGER(mpi) :: npdblk
9658 REAL(mpd) :: adum(1)
9659 REAL(mpd) :: ddum(1)
9661 INTEGER(mpl) :: iparm(64)
9662 REAL(mpd),
ALLOCATABLE :: b( : )
9663 REAL(mpd),
ALLOCATABLE :: x( : )
9682 WRITE(*,*)
'MSPARDISO: number of rows to fill up = ', nfill
9695 CALL pardiso_64(pt, maxfct, mnum, mtype, phase, int(npdblk,mpl), adum, idum, idum, &
9696 idum, nrhs, iparm, msglvl, ddum, ddum, error)
9697 IF (error /= 0)
THEN
9698 WRITE(lun,*)
'The following ERROR was detected: ', error
9699 WRITE(*,
'(A,2I10)')
' PARDISO release failed (phase, error): ', phase, error
9700 IF (
ipddbg == 0)
WRITE(*,*)
' rerun with "debugPARDISO" for more info'
9701 CALL peend(40,
'Aborted, other error: PARDISO release')
9702 stop
'MSPARDISO: stopping due to error in PARDISO release'
9719 IF (iparm(1) == 0)
WRITE(lun,*)
'PARDISO using defaults '
9720 IF (iparm(43) /= 0)
THEN
9721 WRITE(lun,*)
'PARDISO: computation of the diagonal of inverse matrix not implemented !'
9735 WRITE(
lunlog,*)
'Decomposition of global matrix (A->L*D*L^t)'
9749 WRITE(lun,*)
' iparm(',i,
') =', iparm(i)
9752 CALL pardiso_64(pt, maxfct, mnum, mtype, phase, int(npdblk,mpl),
globalmatd,
csr3rowoffsets,
csr3columnlist, &
9753 idum, nrhs, iparm, msglvl, ddum, ddum, error)
9755 WRITE(lun,*)
'PARDISO reordering completed ... '
9756 WRITE(lun,*)
'PARDISO peak memory required (KB)', iparm(15)
9759 WRITE(lun,*)
' iparm(',i,
') =', iparm(i)
9762 IF (error /= 0)
THEN
9763 WRITE(lun,*)
'The following ERROR was detected: ', error
9764 WRITE(*,
'(A,2I10)')
' PARDISO decomposition failed (phase, error): ', phase, error
9765 IF (
ipddbg == 0)
WRITE(*,*)
' rerun with "debugPARDISO" for more info'
9766 CALL peend(40,
'Aborted, other error: PARDISO reordering')
9767 stop
'MSPARDISO: stopping due to error in PARDISO reordering'
9769 IF (iparm(60) == 0)
THEN
9774 WRITE(lun,*)
'Size (KB) of allocated memory = ',
ipdmem
9775 WRITE(lun,*)
'Number of nonzeros in factors = ',iparm(18)
9776 WRITE(lun,*)
'Number of factorization MFLOPS = ',iparm(19)
9781 CALL pardiso_64(pt, maxfct, mnum, mtype, phase, int(npdblk,mpl),
globalmatd,
csr3rowoffsets,
csr3columnlist, &
9782 idum, nrhs, iparm, msglvl, ddum, ddum, error)
9784 WRITE(lun,*)
'PARDISO factorization completed ... '
9787 WRITE(lun,*)
' iparm(',i,
') =', iparm(i)
9790 IF (error /= 0)
THEN
9791 WRITE(lun,*)
'The following ERROR was detected: ', error
9792 WRITE(*,
'(A,2I10)')
' PARDISO decomposition failed (phase, error): ', phase, error
9793 IF (
ipddbg == 0)
WRITE(*,*)
' rerun with "debugPARDISO" for more info'
9794 CALL peend(40,
'Aborted, other error: PARDISO factorization')
9795 stop
'MSPARDISO: stopping due to error in PARDISO factorization'
9798 IF (iparm(14) > 0) &
9799 WRITE(lun,*)
'Number of perturbed pivots = ',iparm(14)
9800 WRITE(lun,*)
'Number of positive eigenvalues = ',iparm(22)-nfill
9801 WRITE(lun,*)
'Number of negative eigenvalues = ',iparm(23)
9802 ELSE IF (iparm(30) > 0)
THEN
9803 WRITE(lun,*)
'Equation with bad pivot (<=0.) = ',iparm(30)
9812 CALL mpalloc(b,length,
' PARDISO r.h.s')
9813 CALL mpalloc(x,length,
' PARDISO solution')
9818 CALL pardiso_64(pt, maxfct, mnum, mtype, phase, int(npdblk,mpl),
globalmatd,
csr3rowoffsets,
csr3columnlist, &
9819 idum, nrhs, iparm, msglvl, b, x, error)
9824 WRITE(lun,*)
'PARDISO solve completed ... '
9825 IF (error /= 0)
THEN
9826 WRITE(lun,*)
'The following ERROR was detected: ', error
9827 WRITE(*,
'(A,2I10)')
' PARDISO decomposition failed (phase, error): ', phase, error
9828 IF (
ipddbg == 0)
WRITE(*,*)
' rerun with "debugPARDISO" for more info'
9829 CALL peend(40,
'Aborted, other error: PARDISO solve')
9830 stop
'MSPARDISO: stopping due to error in PARDISO solve'
9844 INTEGER(mpi) :: iast
9845 INTEGER(mpi) :: idia
9846 INTEGER(mpi) :: imin
9847 INTEGER(mpl) :: ioff1
9849 INTEGER(mpi) :: last
9851 INTEGER(mpi) :: nmax
9852 INTEGER(mpi) :: nmin
9853 INTEGER(mpi) :: ntop
9875 WRITE(
lunlog,*)
'Shrinkage of global matrix (A->Q^t*A*Q)'
9912 DO WHILE(ntop < nmax)
9916 CALL hmpdef(7,real(nmin,mps),real(ntop,mps),
'log10 of positive eigenvalues')
9927 CALL gmpdef(3,2,
'low-value end of eigenvalues')
9930 CALL gmpxy(3,real(i,mps),evalue)
9944 WRITE(lun,*)
'The first (largest) eigenvalues ...'
9947 WRITE(lun,*)
'The last eigenvalues ... up to',last
9951 WRITE(lun,*)
'The eigenvalues from',
nvgb+1,
' to',
nagb
9955 WRITE(lun,*)
'Log10 + 3 of ',
nfgb,
' eigenvalues in decreasing',
' order'
9956 WRITE(lun,*)
'(for Eigenvalue < 0.001 the value 0.0 is shown)'
9959 'printed for negative eigenvalues'
9962 WRITE(lun,*) last,
' significances: insignificant if ', &
9963 'compatible with N(0,1)'
9992 INTEGER(mpl) :: ioff1
9993 INTEGER(mpl) :: ioff2
10029 INTEGER(mpi) :: istop
10030 INTEGER(mpi) :: itn
10031 INTEGER(mpi) :: itnlim
10032 INTEGER(mpi) :: lun
10033 INTEGER(mpi) :: nout
10034 INTEGER(mpi) :: nrkd
10035 INTEGER(mpi) :: nrkd2
10041 REAL(mpd) :: arnorm
10080 WRITE(lun,*)
'MMINRS: PRECONS ended ', nrkd
10084 globalcorrections, itnlim, nout, rtol, istop, itn, anorm, acond, rnorm, arnorm, ynorm)
10085 ELSE IF(
mbandw > 0)
THEN
10091 WRITE(lun,*)
'MMINRS: EQUDECS ended ', nrkd, nrkd2
10095 globalcorrections, itnlim, nout, rtol, istop, itn, anorm, acond, rnorm, arnorm, ynorm)
10098 globalcorrections, itnlim, nout, rtol, istop, itn, anorm, acond, rnorm, arnorm, ynorm)
10112 IF (
istopa == 0) print *,
'MINRES: istop=0, exact solution x=0.'
10127 INTEGER(mpi) :: istop
10128 INTEGER(mpi) :: itn
10129 INTEGER(mpi) :: itnlim
10130 INTEGER(mpi) :: lun
10131 INTEGER(mpi) :: nout
10132 INTEGER(mpi) :: nrkd
10133 INTEGER(mpi) :: nrkd2
10136 REAL(mpd) :: mxxnrm
10137 REAL(mpd) :: trcond
10147 mxxnrm = real(
nagb,mpd)/sqrt(epsilon(mxxnrm))
10149 trcond = 1.0_mpd/epsilon(trcond)
10150 ELSE IF(
mrmode == 2)
THEN
10180 WRITE(lun,*)
'MMINRS: PRECONS ended ', nrkd
10184 itnlim=itnlim, rtol=rtol, maxxnorm=mxxnrm, trancond=trcond, &
10186 ELSE IF(
mbandw > 0)
THEN
10192 WRITE(lun,*)
'MMINRS: EQUDECS ended ', nrkd, nrkd2
10197 itnlim=itnlim, rtol=rtol, maxxnorm=mxxnrm, trancond=trcond, &
10201 itnlim=itnlim, rtol=rtol, maxxnorm=mxxnrm, trancond=trcond, &
10216 IF (
istopa == 3) print *,
'MINRES: istop=0, exact solution x=0.'
10232 INTEGER(mpi),
INTENT(IN) :: n
10233 REAL(mpd),
INTENT(IN) :: x(n)
10234 REAL(mpd),
INTENT(OUT) :: y(n)
10254 INTEGER(mpi),
INTENT(IN) :: n
10255 REAL(mpd),
INTENT(IN) :: x(n)
10256 REAL(mpd),
INTENT(OUT) :: y(n)
10287 REAL(mps) :: concu2
10288 REAL(mps) :: concut
10289 REAL,
DIMENSION(2) :: ta
10292 INTEGER(mpi) :: iact
10293 INTEGER(mpi) :: iagain
10294 INTEGER(mpi) :: idx
10295 INTEGER(mpi) :: info
10297 INTEGER(mpi) :: ipoff
10298 INTEGER(mpi) :: icoff
10299 INTEGER(mpl) :: ioff
10300 INTEGER(mpi) :: itgbi
10301 INTEGER(mpi) :: ivgbi
10302 INTEGER(mpi) :: jcalcm
10304 INTEGER(mpi) :: labelg
10305 INTEGER(mpi) :: litera
10306 INTEGER(mpl) :: lrej
10307 INTEGER(mpi) :: lun
10308 INTEGER(mpi) :: lunp
10309 INTEGER(mpi) :: minf
10310 INTEGER(mpi) :: mrati
10311 INTEGER(mpi) :: nan
10312 INTEGER(mpi) :: ncon
10313 INTEGER(mpi) :: nfaci
10314 INTEGER(mpi) :: nloopsol
10315 INTEGER(mpi) :: npar
10316 INTEGER(mpi) :: nrati
10317 INTEGER(mpl) :: nrej
10318 INTEGER(mpi) :: nsol
10319 INTEGER(mpi) :: inone
10321 INTEGER(mpi) :: infolp
10322 INTEGER(mpi) :: nfit
10323 INTEGER(mpl) :: imoff
10327 REAL(mpd) :: dratio
10328 REAL(mpd) :: dwmean
10337 LOGICAL :: warnerss
10338 LOGICAL :: warners3
10340 CHARACTER (LEN=7) :: cratio
10341 CHARACTER (LEN=7) :: cfacin
10342 CHARACTER (LEN=7) :: crjrat
10353 WRITE(lunp,*)
'Solution algorithm: '
10354 WRITE(lunp,121)
'=================================================== '
10357 WRITE(lunp,121)
'solution method:',
'matrix inversion'
10358 ELSE IF(
metsol == 2)
THEN
10359 WRITE(lunp,121)
'solution method:',
'diagonalization'
10360 ELSE IF(
metsol == 3)
THEN
10361 WRITE(lunp,121)
'solution method:',
'decomposition'
10362 ELSE IF(
metsol == 4)
THEN
10363 WRITE(lunp,121)
'solution method:',
'minres (Paige/Saunders)'
10364 ELSE IF(
metsol == 5)
THEN
10365 WRITE(lunp,121)
'solution method:',
'minres-qlp (Choi/Paige/Saunders)'
10367 WRITE(lunp,121)
' ',
' using QR factorization'
10368 ELSE IF(
mrmode == 2)
THEN
10369 WRITE(lunp,121)
' ',
' using QLP factorization'
10371 WRITE(lunp,121)
' ',
' using QR and QLP factorization'
10372 WRITE(lunp,123)
'transition condition',
mrtcnd
10374 ELSE IF(
metsol == 6)
THEN
10375 WRITE(lunp,121)
'solution method:', &
10376 'gmres (generalized minimzation of residuals)'
10378 ELSE IF(
metsol == 7)
THEN
10380 WRITE(lunp,121)
'solution method:',
'LAPACK factorization (DSPTRF)'
10382 WRITE(lunp,121)
'solution method:',
'LAPACK factorization (DPPTRF)'
10384 IF(
ilperr == 1)
WRITE(lunp,121)
' ',
'with error calculation (D??TRI)'
10385 ELSE IF(
metsol == 8)
THEN
10387 WRITE(lunp,121)
'solution method:',
'LAPACK factorization (DSYTRF)'
10389 WRITE(lunp,121)
'solution method:',
'LAPACK factorization (DPOTRF)'
10391 IF(
ilperr == 1)
WRITE(lunp,121)
' ',
'with error calculation (D??TRI)'
10393 ELSE IF(
metsol == 9)
THEN
10395 WRITE(lunp,121)
'solution method:',
'Intel oneMKL PARDISO (sparse matrix (CSR3))'
10397 WRITE(lunp,121)
'solution method:',
'Intel oneMKL PARDISO (sparse matrix (BSR3))'
10402 WRITE(lunp,123)
'convergence limit at Delta F=',
dflim
10403 WRITE(lunp,122)
'maximum number of iterations=',
mitera
10406 WRITE(lunp,122)
'matrix recalculation up to ',
matrit,
'. iteration'
10410 WRITE(lunp,121)
'matrix storage:',
'full'
10411 ELSE IF(
matsto == 2)
THEN
10412 WRITE(lunp,121)
'matrix storage:',
'sparse'
10414 WRITE(lunp,122)
'pre-con band-width parameter=',
mbandw
10416 WRITE(lunp,121)
'pre-conditioning:',
'default'
10417 ELSE IF(
mbandw < 0)
THEN
10418 WRITE(lunp,121)
'pre-conditioning:',
'none!'
10419 ELSE IF(
mbandw > 0)
THEN
10421 WRITE(lunp,121)
'pre-conditioning=',
'skyline-matrix (rank preserving)'
10423 WRITE(lunp,121)
'pre-conditioning=',
'band-matrix'
10428 WRITE(lunp,121)
'using pre-sigmas:',
'no'
10431 WRITE(lunp,124)
'pre-sigmas defined for', &
10432 REAL(100*npresg,mps)/
REAL(nvgb,mps),
' % of variable parameters'
10433 WRITE(lunp,123)
'default pre-sigma=',
regpre
10436 WRITE(lunp,121)
'regularization:',
'no'
10438 WRITE(lunp,121)
'regularization:',
'yes'
10439 WRITE(lunp,123)
'regularization factor=',
regula
10443 WRITE(lunp,121)
'Chi square cut equiv 3 st.dev applied'
10444 WRITE(lunp,123)
'... in first iteration with factor',
chicut
10445 WRITE(lunp,123)
'... in second iteration with factor',
chirem
10446 WRITE(lunp,121)
' (reduced by sqrt in next iterations)'
10449 WRITE(lunp,121)
'Scaling of measurement errors applied'
10450 WRITE(lunp,123)
'... factor for "global" measuements',
dscerr(1)
10451 WRITE(lunp,123)
'... factor for "local" measuements',
dscerr(2)
10454 WRITE(lunp,122)
'Down-weighting of outliers in',
lhuber,
' iterations'
10455 WRITE(lunp,123)
'Cut on downweight fraction',
dwcut
10459121
FORMAT(1x,a40,3x,a)
10460122
FORMAT(1x,a40,3x,i0,a)
10461123
FORMAT(1x,a40,2x,e9.2)
10462124
FORMAT(1x,a40,3x,f5.1,a)
10472 stepl =real(stp,mps)
10484 ELSE IF(
metsol == 2)
THEN
10487 ELSE IF(
metsol == 3)
THEN
10490 ELSE IF(
metsol == 4)
THEN
10493 ELSE IF(
metsol == 5)
THEN
10496 ELSE IF(
metsol == 6)
THEN
10508 WRITE(
lunlog,*)
'Checking feasibility of parameters:'
10509 WRITE(*,*)
'Checking feasibility of parameters:'
10510 CALL feasib(concut,iact)
10512 WRITE(*,102) concut
10513 WRITE(*,*)
' parameters are made feasible'
10514 WRITE(
lunlog,102) concut
10515 WRITE(
lunlog,*)
' parameters are made feasible'
10517 WRITE(*,*)
' parameters are feasible (i.e. satisfy constraints)'
10518 WRITE(
lunlog,*)
' parameters are feasible (i.e. satisfy constraints)'
10526 WRITE(*,*)
'Reading files and accumulating vectors/matrices ...'
10530 WRITE(
lunlog,*)
'Reading files and accumulating vectors/matrices ...'
10547 IF(jcalcm+1 /= 0)
THEN
10556 IF(
iterat /= litera)
THEN
10577 IF (iabs(jcalcm) <= 1)
THEN
10590 IF(3*nrej >
nrecal)
THEN
10592 WRITE(*,*)
'Data records rejected in previous loop: '
10594 WRITE(*,*)
'Too many rejects (>33.3%) - stop'
10595 CALL peend(26,
'Aborted, too many rejects')
10602 IF(abs(
icalcm) == 1)
THEN
10614 ELSE IF(
metsol == 2)
THEN
10616 ELSE IF(
metsol == 3)
THEN
10618 ELSE IF(
metsol == 4)
THEN
10620 ELSE IF(
metsol == 5)
THEN
10622 ELSE IF(
metsol == 6)
THEN
10623 WRITE(*,*)
'... reserved for GMRES (not yet!)'
10626 ELSE IF(
metsol == 7)
THEN
10628 ELSE IF(
metsol == 8)
THEN
10631 ELSE IF(
metsol == 9)
THEN
10645 CALL feasib(concut,iact)
10657 angras=real(db/sqrt(db1*db2),mps)
10658 dbsig=16.0_mpd*sqrt(max(db1,db2))*epsilon(db)
10664 lsflag=lsflag .AND. (db > dbsig)
10671 ELSE IF(
metsol == 2)
THEN
10674 ELSE IF(
metsol == 3)
THEN
10677 ELSE IF(
metsol == 4)
THEN
10680 ELSE IF(
metsol == 5)
THEN
10683 ELSE IF(
metsol == 6)
THEN
10693 IF(db <= -dbsig)
THEN
10694 WRITE(*,*)
'Function not decreasing:',db
10695 IF(db > -1.0e-3_mpd)
THEN
10697 IF (iagain <= 1)
THEN
10698 WRITE(*,*)
'... again matrix calculation'
10702 WRITE(*,*)
'... aborting iterations'
10706 WRITE(*,*)
'... stopping iterations'
10729 IF (
nloopn == nloopsol)
THEN
10735 stepl=real(stp,mps)
10745 WRITE(*,*)
'Result vector containes ', nan,
' NaNs - stop'
10746 CALL peend(25,
'Aborted, result vector contains NaNs')
10753 WRITE(*,*)
'Subito! Exit after first step.'
10758 WRITE(*,*)
'INFO=0 should not happen (line search input err)'
10759 IF (iagain <= 0)
THEN
10764 IF(info < 0 .OR.
nloopn == nloopsol) cycle
10768 CALL feasib(concut,iact)
10769 IF(iact /= 0.OR.
chicut > 1.0)
THEN
10784 IF(sum(
nrejec) /= 0)
THEN
10786 WRITE(*,*)
'Data records rejected in last loop: '
10808 nfit=npar+ncon;
IF (
icelim > 0) nfit=npar-ncon
10809 IF (nfit > npar)
THEN
10812 WRITE(
lunlog,*)
'Inverse of global matrix from LDLt factorization'
10818 IF(infolp /= 0) print *,
' DSPTRI failed: ', infolp
10823 CALL dsytri(
'U',int(nfit,mpl),
globalmatd(imoff+1:),int(nfit,mpl),&
10825 IF(infolp /= 0) print *,
' DSYTRI failed: ', infolp
10831 WRITE(
lunlog,*)
'Inverse of global matrix from LLt factorization'
10836 CALL dpptri(
'U',int(nfit,mpl),
globalmatd(imoff+1:),infolp)
10837 IF(infolp /= 0) print *,
' DPPTRI failed: ', infolp
10841 CALL dpotri(
'U',int(nfit,mpl),
globalmatd(imoff+1:),int(npar,mpl),infolp)
10842 IF(infolp /= 0) print *,
' DPOTRI failed: ', infolp
10859 DO i=npar-ncon+1,npar
10866 WRITE(
lunlog,*)
'Expansion of global matrix (A->Q*A*Q^t)'
10882 catio=real(dratio,mps)
10886 mrati=nint(100.0*catio,mpi)
10890 IF (
nfilw <= 0)
THEN
10891 WRITE(lunp,*)
'Sum(Chi^2)/Sum(Ndf) =',
fvalue
10893 WRITE(lunp,*)
' =',dratio
10895 WRITE(lunp,*)
'Sum(W*Chi^2)/Sum(Ndf)/<W> =',
fvalue
10897 WRITE(lunp,*)
' /',dwmean
10898 WRITE(lunp,*)
' =',dratio
10902 ' with correction for down-weighting ',catio
10923 nrati=nint(10000.0*real(nrej,mps)/real(
nrecal,mps),mpi)
10924 WRITE(crjrat,197) 0.01_mpd*real(nrati,mpd)
10925 nfaci=nint(100.0*sqrt(catio),mpi)
10927 WRITE(cratio,197) 0.01_mpd*real(mrati,mpd)
10928 WRITE(cfacin,197) 0.01_mpd*real(nfaci,mpd)
10931 IF(mrati < 90.OR.mrati > 110) warner=.true.
10932 IF(nrati > 100) warner=.true.
10933 IF(
ncgbe /= 0) warner=.true.
10935 IF(
nalow /= 0) warners=.true.
10937 IF(
nmiss1 /= 0) warnerss=.true.
10938 IF(iagain /= 0) warnerss=.true.
10939 IF(
ndefec /= 0) warnerss=.true.
10940 IF(
ndefpg /= 0) warnerss=.true.
10942 IF(
nrderr /= 0) warners3=.true.
10944 IF(warner.OR.warners.OR.warnerss.Or.warners3)
THEN
10947 WRITE(*,199)
'WarningWarningWarningWarningWarningWarningWarningWarningWar'
10948 WRITE(*,199)
'arningWarningWarningWarningWarningWarningWarningWarningWarn'
10949 WRITE(*,199)
'rningWarningWarningWarningWarningWarningWarningWarningWarni'
10950 WRITE(*,199)
'ningWarningWarningWarningWarningWarningWarningWarningWarnin'
10951 WRITE(*,199)
'ingWarningWarningWarningWarningWarningWarningWarningWarning'
10952 WRITE(*,199)
'ngWarningWarningWarningWarningWarningWarningWarningWarningW'
10953 WRITE(*,199)
'gWarningWarningWarningWarningWarningWarningWarningWarningWa'
10955 IF(mrati < 90.OR.mrati > 110)
THEN
10957 WRITE(*,*)
' Chi^2/Ndf = ',cratio,
' (should be close to 1)'
10958 WRITE(*,*)
' => multiply all input standard ', &
10959 'deviations by factor',cfacin
10962 IF(nrati > 100)
THEN
10964 WRITE(*,*)
' Fraction of rejects =',crjrat,
' %', &
10965 ' (should be far below 1 %)'
10966 WRITE(*,*)
' => please provide correct mille data'
10970 IF(iagain /= 0)
THEN
10972 WRITE(*,*)
' Matrix not positiv definite '// &
10973 '(function not decreasing)'
10974 WRITE(*,*)
' => please provide correct mille data'
10979 WRITE(*,*)
' Rank defect =',
ndefec, &
10980 ' for global matrix, should be 0'
10981 WRITE(*,*)
' => please provide correct mille data'
10986 WRITE(*,*)
' Rank defect for',
ndefpg, &
10987 ' parameter groups, should be 0'
10988 WRITE(*,*)
' => please provide correct mille data'
10993 WRITE(*,*)
' Rank defect =',
nmiss1, &
10994 ' for constraint equations, should be 0'
10995 WRITE(*,*)
' => please correct constraint definition'
10998 IF(
ncgbe /= 0)
THEN
11000 WRITE(*,*)
' Number of empty constraints =',
ncgbe,
', should be 0'
11001 WRITE(*,*)
' => please check constraint definition, mille data'
11004 IF(
nxlow /= 0)
THEN
11006 WRITE(*,*)
' Possible rank defects =',
nxlow,
' for global matrix'
11007 WRITE(*,*)
' (too few accepted entries)'
11008 WRITE(*,*)
' => please check mille data and ENTRIES cut'
11011 IF(
nalow /= 0)
THEN
11013 WRITE(*,*)
' Possible bad elements =',
nalow,
' in global vector'
11014 WRITE(*,*)
' (toos few accepted entries)'
11015 IF(
ipcntr > 0)
WRITE(*,*)
' (indicated in millepede.res by counts<0)'
11016 WRITE(*,*)
' => please check mille data and ENTRIES cut'
11021 WRITE(*,*)
' Binary file(s) with read errors =',
nrderr,
' (treated as EOF)'
11022 WRITE(*,*)
' => please check mille data'
11026 WRITE(*,199)
'WarningWarningWarningWarningWarningWarningWarningWarningWar'
11027 WRITE(*,199)
'arningWarningWarningWarningWarningWarningWarningWarningWarn'
11028 WRITE(*,199)
'rningWarningWarningWarningWarningWarningWarningWarningWarni'
11029 WRITE(*,199)
'ningWarningWarningWarningWarningWarningWarningWarningWarnin'
11030 WRITE(*,199)
'ingWarningWarningWarningWarningWarningWarningWarningWarning'
11031 WRITE(*,199)
'ngWarningWarningWarningWarningWarningWarningWarningWarningW'
11032 WRITE(*,199)
'gWarningWarningWarningWarningWarningWarningWarningWarningWa'
11043 ELSE IF(
metsol == 2)
THEN
11045 ELSE IF(
metsol == 3)
THEN
11051 IF(labelg == 0) cycle
11052 itgbi=inone(labelg)
11055 IF(ivgbi < 0) ivgbi=0
11056 IF(ivgbi == 0) cycle
11065 ELSE IF(
metsol == 6)
THEN
11068 ELSE IF(
metsol == 7)
THEN
11076 CALL peend(4,
'Ended with severe warnings (bad binary file(s))')
11077 ELSE IF (warnerss)
THEN
11078 CALL peend(3,
'Ended with severe warnings (bad global matrix)')
11079 ELSE IF (warners)
THEN
11080 CALL peend(2,
'Ended with severe warnings (insufficient measurements)')
11081 ELSE IF (warner)
THEN
11082 CALL peend(1,
'Ended with warnings (bad measurements)')
11084 CALL peend(0,
'Ended normally')
11087102
FORMAT(
' Call FEASIB with cut=',g10.3)
11105 INTEGER(mpi) :: kfl
11106 INTEGER(mpi) :: kmin
11107 INTEGER(mpi) :: kmax
11108 INTEGER(mpi) :: nrc
11109 INTEGER(mpl) :: nrej
11115 REAL(mpd) :: sumallw
11116 REAL(mpd) :: sumrejw
11118 sumallw=0.; sumrejw=0.;
11127 sumallw=sumallw+real(nrc,mpd)*
wfd(kfl)
11128 sumrejw=sumrejw+real(nrej,mpd)*
wfd(kfl)
11129 frac=real(nrej,mps)/real(nrc,mps)
11130 IF (frac > fmax)
THEN
11134 IF (frac < fmin)
THEN
11141 WRITE(*,
"(' Weighted fraction =',F8.2,' %')") 100.*sumrejw/sumallw
11142 IF (
nfilb > 1)
THEN
11143 WRITE(*,
"(' File with max. fraction ',I6,' :',F8.2,' %')") kmax, 100.*fmax
11144 WRITE(*,
"(' File with min. fraction ',I6,' :',F8.2,' %')") kmin, 100.*fmin
11170 INTEGER(mpi) :: iargc
11173 INTEGER(mpi) :: ierrf
11174 INTEGER(mpi) :: ieq
11175 INTEGER(mpi) :: ifilb
11176 INTEGER(mpi) :: ioff
11177 INTEGER(mpi) :: iopt
11178 INTEGER(mpi) :: ios
11179 INTEGER(mpi) :: iosum
11182 INTEGER(mpi) :: mat
11183 INTEGER(mpi) :: nab
11184 INTEGER(mpi) :: nline
11185 INTEGER(mpi) :: npat
11186 INTEGER(mpi) :: ntext
11188 INTEGER(mpi) :: nuf
11189 INTEGER(mpi) :: nums
11190 INTEGER(mpi) :: nufile
11191 INTEGER(mpi) :: lenfileInfo
11192 INTEGER(mpi) :: lenFileNames
11193 INTEGER(mpi) :: matint
11194 INTEGER(mpi),
DIMENSION(:,:),
ALLOCATABLE :: vecfileInfo
11195 INTEGER(mpi),
DIMENSION(:,:),
ALLOCATABLE :: tempArray
11196 INTEGER(mpl) :: rows
11197 INTEGER(mpl) :: cols
11198 INTEGER(mpl) :: newcols
11199 INTEGER(mpl) :: length
11201 CHARACTER (LEN=1024) :: text
11202 CHARACTER (LEN=1024) :: fname
11203 CHARACTER (LEN=14) :: bite(3)
11204 CHARACTER (LEN=32) :: keystx
11205 INTEGER(mpi),
PARAMETER :: mnum=100
11206 REAL(mpd) :: dnum(mnum)
11210 SUBROUTINE initc(nfiles)
BIND(c)
11212 INTEGER(c_int),
INTENT(IN),
VALUE :: nfiles
11213 END SUBROUTINE initc
11218 DATA bite/
'C_binary',
'text ',
'Fortran_binary'/
11233 WRITE(*,*)
'Command line options: '
11234 WRITE(*,*)
'--------------------- '
11236 CALL getarg(i,text)
11237 CALL rltext(text,ia,ib,nab)
11238 WRITE(*,101) i,text(1:nab)
11239 IF(text(ia:ia) /=
'-')
THEN
11240 nu=nufile(text(ia:ib))
11243 WRITE(*,*)
'Second text file in command line - stop'
11244 CALL peend(12,
'Aborted, second text file in command line')
11250 WRITE(*,*)
'Open error for file:',text(ia:ib),
' - stop'
11251 CALL peend(16,
'Aborted, open error for file')
11252 IF(text(ia:ia) /=
'/')
THEN
11253 CALL getenv(
'PWD',text)
11254 CALL rltext(text,ia,ib,nab)
11255 WRITE(*,*)
'PWD:',text(ia:ib)
11260 IF(index(text(ia:ib),
'b') /= 0)
THEN
11262 WRITE(*,*)
'Debugging requested'
11264 it=index(text(ia:ib),
't')
11267 ieq=index(text(ia+it:ib),
'=')+it
11268 IF (it /= ieq)
THEN
11269 IF (index(text(ia+ieq:ib),
'SL0' ) /= 0)
ictest=2
11270 IF (index(text(ia+ieq:ib),
'SLE' ) /= 0)
ictest=3
11271 IF (index(text(ia+ieq:ib),
'BP' ) /= 0)
ictest=4
11272 IF (index(text(ia+ieq:ib),
'BRLF') /= 0)
ictest=5
11273 IF (index(text(ia+ieq:ib),
'BRLC') /= 0)
ictest=6
11276 IF(index(text(ia:ib),
's') /= 0)
isubit=1
11277 IF(index(text(ia:ib),
'f') /= 0)
iforce=1
11278 IF(index(text(ia:ib),
'c') /= 0)
icheck=1
11279 IF(index(text(ia:ib),
'C') /= 0)
icheck=2
11281 IF(i == iargc())
WRITE(*,*)
'--------------------- '
11302 CALL rltext(text,ia,ib,nab)
11303 nu=nufile(text(ia:ib))
11307 CALL peend(10,
'Aborted, no steering file')
11308 stop
'in FILETC: no steering file. .'
11321 WRITE(*,*)
'Listing of steering file: ',
filnam(1:
nfnam)
11322 WRITE(*,*)
'-------------------------'
11325 WRITE(*,*)
'Open error for steering file - stop'
11326 CALL peend(11,
'Aborted, open error for steering file')
11327 IF(
filnam(1:1) /=
'/')
THEN
11328 CALL getenv(
'PWD',text)
11329 CALL rltext(text,ia,ib,nab)
11330 WRITE(*,*)
'PWD:',text(ia:ib)
11339 rows=6; cols=lenfileinfo
11340 CALL mpalloc(vecfileinfo,rows,cols,
'file info from steering')
11343 READ(10,102,iostat=ierrf) text
11344 IF (ierrf < 0)
EXIT
11345 CALL rltext(text,ia,ib,nab)
11347 IF(nline <= 50)
THEN
11348 WRITE(*,101) nline,text(1:nab)
11349 IF(nline == 50)
WRITE(*,*)
' ...'
11353 CALL rltext(text,ia,ib,nab)
11354 IF(ib == ia+2)
THEN
11355 mat=matint(text(ia:ib),
'end',npat,ntext)
11356 IF(mat == max(npat,ntext))
THEN
11359 WRITE(*,*)
' end-statement after',nline,
' text lines'
11364 keystx=
'fortranfiles'
11365 mat=matint(text(ia:ib),keystx,npat,ntext)
11366 IF(mat == max(npat,ntext))
THEN
11373 mat=matint(text(ia:ib),keystx,npat,ntext)
11374 IF(mat == max(npat,ntext))
THEN
11380 keystx=
'closeandreopen'
11381 mat=matint(text(ia:ib),keystx,npat,ntext)
11382 IF(mat == max(npat,ntext))
THEN
11390 iopt=index(text(ia:ib),
' -- ')
11391 IF (iopt > 0) ie=iopt-1
11394 nu=nufile(text(ia:ie))
11396 IF (
nfiles == lenfileinfo)
THEN
11397 CALL mpalloc(temparray,rows,cols,
'temp file info from steering')
11398 temparray=vecfileinfo
11400 lenfileinfo=lenfileinfo*2
11401 newcols=lenfileinfo
11402 CALL mpalloc(vecfileinfo,rows,newcols,
'file info from steering')
11403 vecfileinfo(:,1:cols)=temparray(:,1:cols)
11409 lenfilenames=lenfilenames+ie-ia+1
11410 vecfileinfo(1,
nfiles)=nline
11411 vecfileinfo(2,
nfiles)=nu
11412 vecfileinfo(3,
nfiles)=ia
11413 vecfileinfo(4,
nfiles)=ie
11414 vecfileinfo(5,
nfiles)=iopt
11415 vecfileinfo(6,
nfiles)=ib
11425 CALL mpalloc(
nfd,length,
'file line (in steering)')
11428 length=lenfilenames
11434 READ(10,102,iostat=ierrf) text
11435 IF (ierrf < 0)
EXIT
11437 IF (nline == vecfileinfo(1,i))
THEN
11438 nfd(i)=vecfileinfo(1,i)
11439 mfd(i)=vecfileinfo(2,i)
11440 ia=vecfileinfo(3,i)-1
11441 lfd(i)=vecfileinfo(4,i)-ia
11443 tfd(ioff+k)=text(ia+k:ia+k)
11448 IF (vecfileinfo(5,i) > 0)
THEN
11449 CALL ratext(text(vecfileinfo(5,i)+4:vecfileinfo(6,i)),nums,dnum,mnum)
11450 IF (nums > 0)
ofd(i)=real(dnum(1),mps)
11460 CALL mpalloc(
ifd,length,
'integrated record numbers (=offset)')
11461 CALL mpalloc(
jfd,length,
'number of accepted records')
11462 CALL mpalloc(
kfd,rows,length,
'number of records in file, file order')
11467 CALL mpalloc(
sfd,rows,length,
'start, end of file name in TFD')
11468 CALL mpalloc(
yfd,length,
'modification date')
11471 WRITE(*,*)
'-------------------------'
11477 WRITE(*,*)
'Table of files:'
11478 WRITE(*,*)
'---------------'
11481 WRITE(8,*)
'Text and data files:'
11485 fname(k:k)=
tfd(ioff+k)
11488 IF (
mprint > 1)
WRITE(*,103) i,bite(
mfd(i)),fname(1:
lfd(i))
11489 WRITE(8,103) i,bite(
mfd(i)),fname(1:
lfd(i))
11493 WRITE(*,*)
'---------------'
11507 IF(
mfd(i) == 3)
THEN
11531 IF(
mfd(i) == 1)
THEN
11552 WRITE(*,*)
'Opening of C-files not supported.'
11568 IF(iosum /= 0)
THEN
11569 CALL peend(15,
'Aborted, open error(s) for binary files')
11570 stop
'FILETC: open error '
11572 IF(
nfilb == 0)
THEN
11573 CALL peend(14,
'Aborted, no binary files')
11574 stop
'FILETC: no binary files '
11577 WRITE(*,*)
nfilb,
' binary files opened'
11579 WRITE(*,*)
nfilb,
' binary files opened and closed'
11583103
FORMAT(i3,2x,a14,3x,a)
11646 INTEGER(mpi) :: ierrf
11647 INTEGER(mpi) :: ioff
11648 INTEGER(mpi) :: ios
11649 INTEGER(mpi) :: iosum
11651 INTEGER(mpi) :: mat
11652 INTEGER(mpi) :: nab
11653 INTEGER(mpi) :: nfiln
11654 INTEGER(mpi) :: nline
11655 INTEGER(mpi) :: nlinmx
11656 INTEGER(mpi) :: npat
11657 INTEGER(mpi) :: ntext
11658 INTEGER(mpi) :: matint
11662 CHARACTER (LEN=1024) :: text
11663 CHARACTER (LEN=1024) :: fname
11666 WRITE(*,*)
'Processing text files ...'
11679 IF(
mfd(i) /= 2) cycle
11681 fname(k:k)=
tfd(ia+k)
11683 WRITE(*,*)
'File ',fname(1:
lfd(i))
11684 IF (
mprint > 1)
WRITE(*,*)
' '
11685 OPEN(10,file=fname(1:
lfd(i)),iostat=ios,form=
'FORMATTED')
11687 WRITE(*,*)
'Open error for file ',fname(1:
lfd(i))
11697 READ(10,102,iostat=ierrf) text
11698 IF (ierrf < 0)
THEN
11701 WRITE(*,*)
' end-of-file after',nline,
' text lines'
11705 IF(nline <= nlinmx.AND.
mprint > 1)
THEN
11706 CALL rltext(text,ia,ib,nab)
11708 WRITE(*,101) nline,text(1:nab)
11709 IF(nline == nlinmx)
WRITE(*,*)
' ...'
11712 CALL rltext(text,ia,ib,nab)
11713 IF(ib == ia+2)
THEN
11714 mat=matint(text(ia:ib),
'end',npat,ntext)
11715 IF(mat == max(npat,ntext))
THEN
11718 WRITE(*,*)
' end-statement after',nline,
' text lines'
11724 IF(nfiln <=
nfiles)
THEN
11725 IF(nline ==
nfd(nfiln))
THEN
11740 IF(iosum /= 0)
THEN
11741 CALL peend(16,
'Aborted, open error(s) for text files')
11742 stop
'FILETX: open error(s) in text files '
11745 WRITE(*,*)
'... end of text file processing.'
11750 WRITE(*,*)
lunkno,
' unknown keywords in steering files, ', &
11751 'or file non-existing,'
11752 WRITE(*,*)
' see above!'
11753 WRITE(*,*)
'------------> stop'
11755 CALL peend(13,
'Aborted, unknown keywords in steering file')
11764 ELSE IF(
matsto == 1)
THEN
11766 ELSE IF(
matsto == 2)
THEN
11769 ELSE IF(
metsol == 1)
THEN
11771 ELSE IF(
metsol == 2)
THEN
11773 ELSE IF(
metsol == 3)
THEN
11775 ELSE IF(
metsol == 4)
THEN
11777 ELSE IF(
metsol == 5)
THEN
11779 ELSE IF(
metsol == 6)
THEN
11782 ELSE IF(
metsol == 7)
THEN
11784 ELSE IF(
metsol == 8)
THEN
11787 ELSE IF(
metsol == 9)
THEN
11792 WRITE(*,*)
'MINRES forced with sparse matrix!'
11794 WRITE(*,*)
'MINRES forced with sparse matrix!'
11796 WRITE(*,*)
'MINRES forced with sparse matrix!'
11801 WRITE(*,*)
'MINRES forced with sparse matrix!'
11803 WRITE(*,*)
'MINRES forced with sparse matrix!'
11805 WRITE(*,*)
'MINRES forced with sparse matrix!'
11813 WRITE(*,*)
'Solution method and matrix-storage mode:'
11815 WRITE(*,*)
' METSOL = 1: matrix inversion'
11816 ELSE IF(
metsol == 2)
THEN
11817 WRITE(*,*)
' METSOL = 2: diagonalization'
11818 ELSE IF(
metsol == 3)
THEN
11819 WRITE(*,*)
' METSOL = 3: decomposition'
11820 ELSE IF(
metsol == 4)
THEN
11821 WRITE(*,*)
' METSOL = 4: MINRES'
11822 ELSE IF(
metsol == 5)
THEN
11823 WRITE(*,*)
' METSOL = 5: MINRES-QLP'
11824 ELSE IF(
metsol == 6)
THEN
11825 WRITE(*,*)
' METSOL = 6: GMRES (-> MINRES)'
11827 ELSE IF(
metsol == 7)
THEN
11828 WRITE(*,*)
' METSOL = 7: LAPACK factorization'
11829 ELSE IF(
metsol == 8)
THEN
11830 WRITE(*,*)
' METSOL = 8: LAPACK factorization'
11832 ELSE IF(
metsol == 9)
THEN
11833 WRITE(*,*)
' METSOL = 9: Intel oneMKL PARDISO'
11838 WRITE(*,*)
' with',
mitera,
' iterations'
11841 WRITE(*,*)
' MATSTO = 0: unpacked symmetric matrix, ',
'n*n elements'
11842 ELSEIF(
matsto == 1)
THEN
11843 WRITE(*,*)
' MATSTO = 1: full symmetric matrix, ',
'(n*n+n)/2 elements'
11844 ELSE IF(
matsto == 2)
THEN
11845 WRITE(*,*)
' MATSTO = 2: sparse matrix (custom)'
11846 ELSE IF(
matsto == 3)
THEN
11848 WRITE(*,*)
' MATSTO = 3: sparse matrix (upper triangle, CSR3)'
11850 WRITE(*,*)
' MATSTO = 3: sparse matrix (upper triangle, BSR3)'
11854 WRITE(*,*)
' and band matrix, width',
mbandw
11858 WRITE(*,*)
'Chi square cut equiv 3 st.dev applied ...'
11859 WRITE(*,*)
' in first iteration with factor',
chicut
11860 WRITE(*,*)
' in second iteration with factor',
chirem
11861 WRITE(*,*)
' (reduced by sqrt in next iterations)'
11865 WRITE(*,*)
' Down-weighting of outliers in',
lhuber,
' iterations'
11866 WRITE(*,*)
' Cut on downweight fraction',
dwcut
11869 WRITE(*,*)
'Iterations (solutions) with line search:'
11878 WRITE(*,*)
' All with Chi square cut scaling factor <= 1.'
11909 INTEGER(mpi) :: ios
11913 INTEGER(mpi) :: npat
11914 INTEGER(mpi) :: ntext
11915 INTEGER(mpi) :: nuprae
11918 CHARACTER (LEN=*),
INTENT(INOUT) :: fname
11924 IF(len(fname) > 5)
THEN
11925 IF(fname(1:5) ==
'rfio:') nuprae=1
11926 IF(fname(1:5) ==
'dcap:') nuprae=2
11927 IF(fname(1:5) ==
'root:') nuprae=3
11929 IF(nuprae == 0)
THEN
11930 INQUIRE(file=fname,iostat=ios,exist=ex)
11931 IF(ios /= 0)
nufile=-abs(ios)
11932 IF(ios /= 0)
RETURN
11933 ELSE IF(nuprae == 1)
THEN
11946 nm=
matint(
'xt',fname(l1:ll),npat,ntext)
11949 nm=
matint(
'tx',fname(l1:ll),npat,ntext)
11970 INTEGER(mpi) :: ier
11971 INTEGER(mpi) :: iomp
11974 INTEGER(mpi) :: kkey
11975 INTEGER(mpi) :: label
11976 INTEGER(mpi) :: lkey
11977 INTEGER(mpi) :: mat
11978 INTEGER(mpi) :: miter
11979 INTEGER(mpi) :: nab
11980 INTEGER(mpi) :: nkey
11981 INTEGER(mpi) :: nkeys
11983 INTEGER(mpi) :: nmeth
11984 INTEGER(mpi) :: npat
11985 INTEGER(mpi) :: ntext
11986 INTEGER(mpi) :: nums
11987 INTEGER(mpi) :: matint
11989 CHARACTER (LEN=*),
INTENT(IN) :: text
11990 INTEGER(mpi),
INTENT(IN) :: nline
11994 parameter(nkeys=7,nmeth=10)
11996 parameter(nkeys=6,nmeth=9)
11999 parameter(nkeys=6,nmeth=7)
12001 CHARACTER (LEN=16) :: methxt(nmeth)
12002 CHARACTER (LEN=16) :: keylst(nkeys)
12003 CHARACTER (LEN=32) :: keywrd
12004 CHARACTER (LEN=32) :: keystx
12005 CHARACTER (LEN=itemCLen) :: ctext
12006 INTEGER(mpi),
PARAMETER :: mnum=100
12007 REAL(mpd) :: dnum(mnum)
12010 INTEGER(mpi) :: ipvs
12013 INTEGER(mpi) :: lpvs
12017 SUBROUTINE additem(length,list,label,value)
12019 INTEGER(mpi),
INTENT(IN OUT) :: length
12020 TYPE(listitem),
DIMENSION(:),
INTENT(IN OUT),
ALLOCATABLE :: list
12021 INTEGER(mpi),
INTENT(IN) :: label
12022 REAL(mpd),
INTENT(IN) :: value
12024 SUBROUTINE additemc(length,list,label,text)
12026 INTEGER(mpi),
INTENT(IN OUT) :: length
12027 TYPE(listitemc),
DIMENSION(:),
INTENT(IN OUT),
ALLOCATABLE :: list
12028 INTEGER(mpi),
INTENT(IN) :: label
12029 CHARACTER(LEN = itemCLen),
INTENT(IN) :: text
12031 SUBROUTINE additemi(length,list,label,ivalue)
12033 INTEGER(mpi),
INTENT(IN OUT) :: length
12034 TYPE(listitemi),
DIMENSION(:),
INTENT(IN OUT),
ALLOCATABLE :: list
12035 INTEGER(mpi),
INTENT(IN) :: label
12036 INTEGER(mpi),
INTENT(IN) :: ivalue
12043 DATA keylst/
'unknown',
'parameter',
'constraint',
'measurement',
'method',
'comment',
'pardiso'/
12044 DATA methxt/
'diagonalization',
'inversion',
'fullMINRES',
'sparseMINRES', &
12045 'fullMINRES-QLP',
'sparseMINRES-QLP',
'decomposition',
'fullLAPACK',
'unpackedLAPACK', &
12048 DATA keylst/
'unknown',
'parameter',
'constraint',
'measurement',
'method',
'comment'/
12049 DATA methxt/
'diagonalization',
'inversion',
'fullMINRES',
'sparseMINRES', &
12050 'fullMINRES-QLP',
'sparseMINRES-QLP',
'decomposition',
'fullLAPACK',
'unpackedLAPACK'/
12053 DATA keylst/
'unknown',
'parameter',
'constraint',
'measurement',
'method',
'comment'/
12054 DATA methxt/
'diagonalization',
'inversion',
'fullMINRES',
'sparseMINRES', &
12055 'fullMINRES-QLP',
'sparseMINRES-QLP',
'decomposition'/
12061 CALL rltext(text,ia,ib,nab)
12062 IF(nab == 0)
GOTO 10
12063 CALL ratext(text(1:nab),nums,dnum,mnum)
12065 IF(nums /= 0) nkey=0
12073 keystx=keylst(nkey)
12074 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12075 IF(100*mat >= 80*max(npat,ntext))
GO TO 10
12081 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12082 IF(100*mat >= 80*max(npat,ntext))
THEN
12084 IF(nums > 0) mprint=nint(dnum(1),mpi)
12089 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12090 IF(100*mat >= 80*max(npat,ntext))
THEN
12093 IF(nums > 0) mdebug=nint(dnum(1),mpi)
12094 IF(nums > 1) mdebg2=nint(dnum(2),mpi)
12099 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12100 IF(100*mat >= 80*max(npat,ntext))
THEN
12101 IF(nums > 0 .AND. dnum(1) > 0.5) mreqenf=nint(dnum(1),mpi)
12102 IF(nums > 1 .AND. dnum(2) > 0.5) mreqena=nint(dnum(2),mpi)
12103 IF(nums > 2 .AND. dnum(3) > 0.5) iteren=nint(dnum(1)*dnum(3),mpi)
12107 keystx=
'printrecord'
12108 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12109 IF(100*mat >= 80*max(npat,ntext))
THEN
12110 IF(nums > 0) nrecpr=nint(dnum(1),mpi)
12111 IF(nums > 1) nrecp2=nint(dnum(2),mpi)
12116 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12117 IF(100*mat >= 80*max(npat,ntext))
THEN
12118 IF (nums > 0.AND.dnum(1) > 0.) mxrec=nint(dnum(1),mpi)
12123 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12124 IF(100*mat >= 80*max(npat,ntext))
THEN
12125 IF (nums > 0.AND.dnum(1) >= 0.) ncache=nint(dnum(1),mpi)
12126 IF (nums == 2.AND.dnum(2) > 0..AND.dnum(2) <= 1.0) &
12127 fcache(1)=real(dnum(2),mps)
12128 IF (nums >= 4)
THEN
12130 fcache(k)=real(dnum(k+1),mps)
12137 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12138 IF(100*mat >= 80*max(npat,ntext))
THEN
12143 chicut=real(dnum(1),mps)
12144 IF(chicut < 1.0) chicut=-1.0
12148 chirem=real(dnum(2),mps)
12149 IF(chirem < 1.0) chirem=1.0
12150 IF(chicut >= 1.0) chirem=min(chirem,chicut)
12158 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12159 IF(100*mat >= 80*max(npat,ntext))
THEN
12160 IF(nums > 0) chhuge=real(dnum(1),mps)
12161 IF(chhuge < 1.0) chhuge=1.0
12166 keystx=
'linesearch'
12167 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12168 IF(100*mat >= 80*max(npat,ntext))
THEN
12169 IF(nums > 0) lsearch=nint(dnum(1),mpi)
12174 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12175 IF(100*mat >= 80*max(npat,ntext))
THEN
12176 IF(nums > 0) lfitnp=nint(dnum(1),mpi)
12177 IF(nums > 1) lfitbb=nint(dnum(2),mpi)
12181 keystx=
'regularization'
12182 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12183 IF(100*mat >= 80*max(npat,ntext))
THEN
12185 regula=real(dnum(1),mps)
12186 IF(nums >= 2) regpre=real(dnum(2),mps)
12190 keystx=
'regularisation'
12191 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12192 IF(100*mat >= 80*max(npat,ntext))
THEN
12194 regula=real(dnum(1),mps)
12195 IF(nums >= 2) regpre=real(dnum(2),mps)
12200 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12201 IF(100*mat >= 80*max(npat,ntext))
THEN
12202 regpre=real(dnum(1),mps)
12207 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12208 IF(100*mat >= 80*max(npat,ntext))
THEN
12209 matrit=nint(dnum(1),mpi)
12214 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12215 IF(100*mat >= 80*max(npat,ntext))
THEN
12217 IF (nums > 0.AND.dnum(1) > 0.) matmon=nint(dnum(1),mpi)
12222 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12223 IF(100*mat >= 80*max(npat,ntext))
THEN
12224 IF(nums > 0) mbandw=nint(dnum(1),mpi)
12225 IF(mbandw < 0) mbandw=-1
12226 IF(nums > 1) lprecm=nint(dnum(2),mpi)
12250 keystx=
'outlierdownweighting'
12251 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12252 IF(100*mat >= 80*max(npat,ntext))
THEN
12253 lhuber=nint(dnum(1),mpi)
12254 IF(lhuber > 0.AND.lhuber <= 2) lhuber=2
12258 keystx=
'dwfractioncut'
12259 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12260 IF(100*mat >= 80*max(npat,ntext))
THEN
12261 dwcut=real(dnum(1),mps)
12262 IF(dwcut > 0.5) dwcut=0.5
12266 keystx=
'maxlocalcond'
12267 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12268 IF(100*mat >= 80*max(npat,ntext))
THEN
12269 IF (nums > 0.AND.dnum(1) > 0.0) cndlmx=real(dnum(1),mps)
12274 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12275 IF(100*mat >= 80*max(npat,ntext))
THEN
12276 prange=abs(real(dnum(1),mps))
12281 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12282 IF(100*mat >= 80*max(npat,ntext))
THEN
12288 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12289 IF(100*mat >= 80*max(npat,ntext))
THEN
12294 keystx=
'memorydebug'
12295 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12296 IF(100*mat >= 80*max(npat,ntext))
THEN
12298 IF (nums > 0.AND.dnum(1) > 0.0) memdbg=nint(dnum(1),mpi)
12302 keystx=
'globalcorr'
12303 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12304 IF(100*mat >= 80*max(npat,ntext))
THEN
12309 keystx=
'printcounts'
12310 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12311 IF(100*mat >= 80*max(npat,ntext))
THEN
12313 IF (nums > 0) ipcntr=nint(dnum(1),mpi)
12317 keystx=
'weightedcons'
12318 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12319 IF(100*mat >= 80*max(npat,ntext))
THEN
12321 IF (nums > 0) iwcons=nint(dnum(1),mpi)
12325 keystx=
'skipemptycons'
12326 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12327 IF(100*mat >= 80*max(npat,ntext))
THEN
12332 keystx=
'resolveredundancycons'
12333 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12334 IF(100*mat >= 80*max(npat,ntext))
THEN
12339 keystx=
'withelimination'
12340 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12341 IF(100*mat >= 80*max(npat,ntext))
THEN
12347 keystx=
'withLAPACKelimination'
12348 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12349 IF(100*mat >= 80*max(npat,ntext))
THEN
12355 keystx=
'withmultipliers'
12356 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12357 IF(100*mat >= 80*max(npat,ntext))
THEN
12362 keystx=
'checkinput'
12363 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12364 IF(100*mat >= 80*max(npat,ntext))
THEN
12366 IF (nums > 0) icheck=nint(dnum(1),mpi)
12370 keystx=
'checkparametergroups'
12371 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12372 IF(100*mat >= 80*max(npat,ntext))
THEN
12377 keystx=
'monitorresiduals'
12378 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12379 IF(100*mat >= 80*max(npat,ntext))
THEN
12381 IF (nums > 0) imonit=nint(dnum(1),mpi)
12382 IF (nums > 1) measbins=max(measbins,nint(dnum(2),mpi))
12386 keystx=
'monitorpulls'
12387 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12388 IF(100*mat >= 80*max(npat,ntext))
THEN
12391 IF (nums > 0) imonit=nint(dnum(1),mpi)
12392 IF (nums > 1) measbins=max(measbins,nint(dnum(2),mpi))
12396 keystx=
'monitorprogress'
12397 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12398 IF(100*mat >= 80*max(npat,ntext))
THEN
12401 IF (nums > 0) monpg1=max(1,nint(dnum(1),mpi))
12402 IF (nums > 1) monpg2=max(1,nint(dnum(2),mpi))
12406 keystx=
'scaleerrors'
12407 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12408 IF(100*mat >= 80*max(npat,ntext))
THEN
12410 IF (nums > 0) dscerr(1:2)=dnum(1)
12411 IF (nums > 1) dscerr(2)=dnum(2)
12415 keystx=
'iterateentries'
12416 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12417 IF(100*mat >= 80*max(npat,ntext))
THEN
12418 iteren=huge(iteren)
12419 IF (nums > 0) iteren=nint(dnum(1),mpi)
12424 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12425 IF(100*mat >= 80*max(npat,ntext))
THEN
12433 WRITE(*,*)
'WARNING: multithreading not available'
12439 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12440 IF(100*mat >= 80*max(npat,ntext))
THEN
12441 WRITE(*,*)
'WARNING: keyword COMPRESS is obsolete (compression is default)'
12453 keystx=
'countrecords'
12454 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12455 IF(100*mat >= 80*max(npat,ntext))
THEN
12461 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12462 IF(100*mat >= 80*max(npat,ntext).AND.mnrsel < 100)
THEN
12463 nl=min(nums,100-mnrsel)
12465 lbmnrs(mnrsel+k)=nint(dnum(k),mpi)
12471 keystx=
'pairentries'
12472 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12473 IF(100*mat >= 80*max(npat,ntext))
THEN
12476 IF (nums > 0.AND.dnum(1) > 0.0)
THEN
12477 mreqpe=nint(dnum(1),mpi)
12478 IF (nums >= 2.AND.dnum(2) >= dnum(1)) mhispe=nint(dnum(2),mpi)
12479 IF (nums >= 3.AND.dnum(3) >= dnum(1)) msngpe=nint(dnum(3),mpi)
12485 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12486 IF(100*mat >= 80*max(npat,ntext))
THEN
12487 wolfc1=real(dnum(1),mps)
12488 wolfc2=real(dnum(2),mps)
12495 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12496 IF(100*mat >= 80*max(npat,ntext))
THEN
12498 IF (dnum(1) < 1.0e-10_mpd.OR.dnum(1) > 1.0e-04_mpd)
THEN
12499 WRITE(*,*)
'ERROR: need 1.0D-10 <= MRESTL ', &
12500 '<= 1.0D-04, but get ', dnum(1)
12509 keystx=
'mrestranscond'
12510 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12511 IF(100*mat >= 80*max(npat,ntext))
THEN
12519 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12520 IF(100*mat >= 80*max(npat,ntext))
THEN
12522 mrmode = int(dnum(1),mpi)
12527 keystx=
'nofeasiblestart'
12528 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12529 IF(100*mat >= 80*max(npat,ntext))
THEN
12535 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12536 IF(100*mat >= 80*max(npat,ntext))
THEN
12541 keystx=
'readerroraseof'
12542 mat=matint(text(ia:ib),keystx,npat,ntext)
12543 IF(100*mat >= 80*max(npat,ntext))
THEN
12549 keystx=
'LAPACKwitherrors'
12550 mat=matint(text(ia:ib),keystx,npat,ntext)
12551 IF(100*mat >= 80*max(npat,ntext))
THEN
12556 keystx=
'debugPARDISO'
12557 mat=matint(text(ia:ib),keystx,npat,ntext)
12558 IF(100*mat >= 80*max(npat,ntext))
THEN
12563 keystx=
'blocksizePARDISO'
12564 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12565 IF(100*mat >= 80*max(npat,ntext).AND.mnrsel < 100)
THEN
12566 nl=min(nums,10-mpdbsz)
12568 IF (nint(dnum(k),mpi) > 0)
THEN
12569 IF (mpdbsz == 0)
THEN
12571 ipdbsz(mpdbsz)=nint(dnum(k),mpi)
12572 ELSE IF (nint(dnum(k),mpi) > ipdbsz(mpdbsz))
THEN
12574 ipdbsz(mpdbsz)=nint(dnum(k),mpi)
12582 keystx=
'fortranfiles'
12583 mat=matint(text(ia:ib),keystx,npat,ntext)
12584 IF(mat == max(npat,ntext))
RETURN
12587 mat=matint(text(ia:ib),keystx,npat,ntext)
12588 IF(mat == max(npat,ntext))
RETURN
12590 keystx=
'closeandreopen'
12591 mat=matint(text(ia:ib),keystx,npat,ntext)
12592 IF(mat == max(npat,ntext))
RETURN
12596 IF(nums /= 0) nkey=0
12599 WRITE(*,*)
'**************************************************'
12601 WRITE(*,*)
'Unknown keyword(s): ',text(1:min(nab,50))
12603 WRITE(*,*)
'**************************************************'
1262510
IF(nkey > 0)
THEN
12629 lpvs=nint(dnum(1),mpi)
12631 CALL additem(lenparameters,listparameters,lpvs,dnum(2))
12632 CALL additem(lenpresigmas,listpresigmas,lpvs,dnum(3))
12634 WRITE(*,*)
'Line',nline,
' error, label=',lpvs
12636 ELSE IF(nums /= 0)
THEN
12638 WRITE(*,*)
'Wrong text in line',nline
12639 WRITE(*,*)
'Status: new parameter'
12640 WRITE(*,*)
'> ',text(1:nab)
12642 ELSE IF(lkey == 3)
THEN
12644 IF(nums >= 1.AND.nums <= 2)
THEN
12646 CALL additem(lenconstraints,listconstraints,lpvs,dnum(1))
12648 IF(iwcons > 0) lpvs=-2
12650 IF(nums == 2) plvs=dnum(2)
12651 CALL additem(lenconstraints,listconstraints,lpvs,plvs)
12654 WRITE(*,*)
'Wrong text in line',nline
12655 WRITE(*,*)
'Status: new keyword constraint'
12656 WRITE(*,*)
'> ',text(1:nab)
12658 ELSE IF(lkey == 4)
THEN
12660 nummeasurements=nummeasurements+1
12662 CALL additem(lenmeasurements,listmeasurements,lpvs,dnum(1))
12664 CALL additem(lenmeasurements,listmeasurements,lpvs,dnum(2))
12667 WRITE(*,*)
'Wrong text in line',nline
12668 WRITE(*,*)
'Status: new keyword measurement'
12669 WRITE(*,*)
'> ',text(1:nab)
12671 ELSE IF(lkey == 5.AND.
keyb <
keyc)
THEN
12673 IF(nums >= 1) miter=nint(dnum(1),mpi)
12674 IF(miter >= 1) mitera=miter
12675 dflim=real(dnum(2),mps)
12679 mat=matint(text(
keyb+1:
keyc),keystx,npat,ntext)
12680 IF(100*mat >= 80*max(npat,ntext))
THEN
12684 ELSE IF(i == 2)
THEN
12687 ELSE IF(i == 3)
THEN
12690 ELSE IF(i == 4)
THEN
12693 ELSE IF(i == 5)
THEN
12696 ELSE IF(i == 6)
THEN
12699 ELSE IF(i == 7)
THEN
12703 ELSE IF(i == 8)
THEN
12706 ELSE IF(i == 9)
THEN
12710 ELSE IF(i == 10)
THEN
12719 ELSE IF(nkey == 0)
THEN
12722 lpvs=nint(dnum(1),mpi)
12724 CALL additem(lenparameters,listparameters,lpvs,dnum(2))
12725 CALL additem(lenpresigmas,listpresigmas,lpvs,dnum(3))
12727 WRITE(*,*)
'Line',nline,
' error, label=',lpvs
12729 ELSE IF(nums > 1.AND.nums < 3)
THEN
12731 WRITE(*,*)
'Wrong text in line',nline
12732 WRITE(*,*)
'Status continuation parameter'
12733 WRITE(*,*)
'> ',text(1:nab)
12736 ELSE IF(lkey == 3)
THEN
12739 label=nint(dnum(i),mpi)
12740 IF(label <= 0) ier=1
12742 IF(mod(nums,2) /= 0) ier=1
12745 lpvs=nint(dnum(i),mpi)
12747 CALL additem(lenconstraints,listconstraints,lpvs,plvs)
12751 WRITE(*,*)
'Wrong text in line',nline
12752 WRITE(*,*)
'Status continuation constraint'
12753 WRITE(*,*)
'> ',text(1:nab)
12756 ELSE IF(lkey == 4)
THEN
12760 label=nint(dnum(i),mpi)
12761 IF(label <= 0) ier=1
12763 IF(mod(nums,2) /= 0) ier=1
12767 lpvs=nint(dnum(i),mpi)
12769 CALL additem(lenmeasurements,listmeasurements,lpvs,plvs)
12773 WRITE(*,*)
'Wrong text in line',nline
12774 WRITE(*,*)
'Status continuation measurement'
12775 WRITE(*,*)
'> ',text(1:nab)
12777 ELSE IF(lkey == 6)
THEN
12779 lpvs=nint(dnum(1),mpi)
12783 IF (text(j:j) ==
' ')
EXIT
12786 CALL additemc(lencomments,listcomments,lpvs,ctext)
12788 WRITE(*,*)
'Line',nline,
' error, label=',lpvs
12790 ELSE IF(nums /= 0)
THEN
12792 WRITE(*,*)
'Wrong text in line',nline
12793 WRITE(*,*)
'Status: continuation comment'
12794 WRITE(*,*)
'> ',text(1:nab)
12798 ELSE IF(lkey == 7)
THEN
12801 label=nint(dnum(i),mpi)
12802 IF(label <= 0.OR.label > 64) ier=1
12804 IF(mod(nums,2) /= 0) ier=1
12808 lpvs=nint(dnum(i),mpi)
12809 ipvs=nint(dnum(i+1),mpi)
12810 CALL additemi(lenpardiso,listpardiso,lpvs,ipvs)
12814 WRITE(*,*)
'Wrong text in line',nline
12815 WRITE(*,*)
'Status continuation measurement'
12816 WRITE(*,*)
'> ',text(1:nab)
12835 INTEGER(mpi),
INTENT(IN OUT) :: length
12836 TYPE(
listitem),
DIMENSION(:),
INTENT(IN OUT),
ALLOCATABLE :: list
12837 INTEGER(mpi),
INTENT(IN) :: label
12838 REAL(mpd),
INTENT(IN) :: value
12840 INTEGER(mpl) :: newSize
12841 INTEGER(mpl) :: oldSize
12842 TYPE(
listitem),
DIMENSION(:),
ALLOCATABLE :: tempList
12844 IF (label > 0.AND.
value == 0.)
RETURN
12845 IF (length == 0 )
THEN
12847 CALL mpalloc(list,newsize,
' list ')
12849 oldsize=
size(list,kind=
mpl)
12850 IF (length >= oldsize)
THEN
12851 newsize = oldsize + oldsize/5 + 100
12852 CALL mpalloc(templist,oldsize,
' temp. list ')
12855 CALL mpalloc(list,newsize,
' list ')
12856 list(1:oldsize)=templist(1:oldsize)
12861 list(length)%label=label
12862 list(length)%value=
value
12877 INTEGER(mpi),
INTENT(IN OUT) :: length
12878 TYPE(
listitemc),
DIMENSION(:),
INTENT(IN OUT),
ALLOCATABLE :: list
12879 INTEGER(mpi),
INTENT(IN) :: label
12880 CHARACTER(len = itemCLen),
INTENT(IN) :: text
12882 INTEGER(mpl) :: newSize
12883 INTEGER(mpl) :: oldSize
12884 TYPE(
listitemc),
DIMENSION(:),
ALLOCATABLE :: tempList
12886 IF (label > 0.AND.text ==
'')
RETURN
12887 IF (length == 0 )
THEN
12889 CALL mpalloc(list,newsize,
' list ')
12891 oldsize=
size(list,kind=
mpl)
12892 IF (length >= oldsize)
THEN
12893 newsize = oldsize + oldsize/5 + 100
12894 CALL mpalloc(templist,oldsize,
' temp. list ')
12897 CALL mpalloc(list,newsize,
' list ')
12898 list(1:oldsize)=templist(1:oldsize)
12903 list(length)%label=label
12904 list(length)%text=text
12919 INTEGER(mpi),
INTENT(IN OUT) :: length
12920 TYPE(
listitemi),
DIMENSION(:),
INTENT(IN OUT),
ALLOCATABLE :: list
12921 INTEGER(mpi),
INTENT(IN) :: label
12922 INTEGER(mpi),
INTENT(IN) :: ivalue
12924 INTEGER(mpl) :: newSize
12925 INTEGER(mpl) :: oldSize
12926 TYPE(
listitemi),
DIMENSION(:),
ALLOCATABLE :: tempList
12928 IF (length == 0 )
THEN
12930 CALL mpalloc(list,newsize,
' list ')
12932 oldsize=
size(list,kind=
mpl)
12933 IF (length >= oldsize)
THEN
12934 newsize = oldsize + oldsize/5 + 100
12935 CALL mpalloc(templist,oldsize,
' temp. list ')
12938 CALL mpalloc(list,newsize,
' list ')
12939 list(1:oldsize)=templist(1:oldsize)
12944 list(length)%label=label
12945 list(length)%ivalue=ivalue
12959 CHARACTER (LEN=*),
INTENT(IN) :: text
12960 CHARACTER (LEN=16) :: textc
12969 textl(ka:kb)=text(1:l)
12973 textc=text(1:l)//
'-end'
12981 textl(ka:kb)=textc(1:l)
13008 INTEGER(mpi),
INTENT(IN) :: lun
13009 CHARACTER (LEN=*),
INTENT(IN) :: fname
13010 CHARACTER (LEN=33) :: nafile
13011 CHARACTER (LEN=33) :: nbfile
13017 CALL peend(17,
'Aborted, file name too long')
13018 stop
'File name too long '
13021 nafile(l+1:l+1)=
'~'
13023 INQUIRE(file=nafile(1:l),exist=ex)
13025 INQUIRE(file=nafile(1:l+1),exist=ex)
13027 CALL system(
'rm '//nafile)
13030 nafile(l+1:l+1)=
' '
13031 CALL system(
'mv '//nafile//nbfile)
13033 OPEN(unit=lun,file=fname)
13044 REAL,
DIMENSION(2) :: ta
13064 IF(ncount > 1)
THEN
13066 nsecd1=int(delta,
mpi)
13068 minut1=nsecd1/60-60*nhour1
13069 secnd1=delta-60*(minut1+60*nhour1)
13071 nsecd2=int(delta,
mpi)
13073 minut2=nsecd2/60-60*nhour2
13074 secnd2=delta-60*(minut2+60*nhour2)
13075 WRITE(*,101) nhour1,minut1,secnd1, nhour2,minut2,secnd2
13080101
FORMAT(i4,
' h',i3,
' min',f5.1,
' sec total',18x,
'elapsed', &
13081 i4,
' h',i3,
' min',f5.1,
' sec')
13095 INTEGER(mpi),
INTENT(IN) :: icode
13096 CHARACTER (LEN=*),
INTENT(IN) :: cmessage
13098 CALL mvopen(9,
'millepede.end')
13099 WRITE(9,101) icode, cmessage
13100101
FORMAT(1x,i4,3x,a)
13104END SUBROUTINE peend
13116 INTEGER(mpi),
INTENT(IN) :: kfile
13117 INTEGER(mpi),
INTENT(IN) :: ithr
13118 INTEGER(mpi),
INTENT(OUT) :: ierr
13120 INTEGER(mpi),
DIMENSION(13) :: ibuff
13121 INTEGER(mpi) :: ioff
13122 INTEGER(mpi) :: ios
13124 INTEGER(mpi) :: lfn
13125 INTEGER(mpi) :: lun
13126 INTEGER(mpi) :: moddate
13127 CHARACTER (LEN=1024) :: fname
13128 CHARACTER (LEN=7) :: cfile
13133 SUBROUTINE openc(filename, lfn, lun, ios)
BIND(c)
13135 CHARACTER(kind=c_char),
DIMENSION(*),
INTENT(IN) :: filename
13136 INTEGER(c_int),
INTENT(IN),
VALUE :: lfn
13137 INTEGER(c_int),
INTENT(IN),
VALUE :: lun
13138 INTEGER(c_int),
INTENT(INOUT) :: ios
13139 END SUBROUTINE openc
13151 fname(k:k)=
tfd(ioff+k)
13156 IF(kfile <=
nfilf)
THEN
13159 OPEN(lun,file=fname(1:lfn),iostat=ios, form=
'UNFORMATTED')
13160 print *,
' lun ', lun, ios
13164 CALL openc(fname(1:lfn),lfn,lun,ios)
13166 WRITE(*,*)
'Opening of C-files not supported.'
13173 WRITE(*,*)
'Open error for file ',fname(1:lfn), ios
13174 IF (moddate /= 0)
THEN
13175 WRITE(cfile,
'(I7)') kfile
13176 CALL peend(15,
'Aborted, open error(s) for binary file ' // cfile)
13177 stop
'PEREAD: open error'
13182 ios=stat(fname(1:lfn),ibuff)
13186 WRITE(*,*)
'STAT error for file ',fname(1:lfn), ios
13190 IF (moddate /= 0)
THEN
13191 IF (ibuff(10) /= moddate)
THEN
13192 WRITE(cfile,
'(I7)') kfile
13193 CALL peend(19,
'Aborted, binary file modified (date) ' // cfile)
13194 stop
'PEREAD: file modified'
13197 yfd(kfile)=ibuff(10)
13212 INTEGER(mpi),
INTENT(IN) :: kfile
13213 INTEGER(mpi),
INTENT(IN) :: ithr
13215 INTEGER(mpi) :: lun
13219 SUBROUTINE closec(lun)
BIND(c)
13221 INTEGER(c_int),
INTENT(IN),
VALUE :: lun
13228 IF(kfile <=
nfilf)
THEN
13247 INTEGER(mpi),
INTENT(IN) :: kfile
13249 INTEGER(mpi) :: lun
13253 SUBROUTINE resetc(lun)
BIND(c)
13255 INTEGER(c_int),
INTENT(IN),
VALUE :: lun
13261 IF (kfile <=
nfilf)
THEN
13280 INTEGER(mpi) :: ipgrp
13281 INTEGER(mpi) :: irank
13282 INTEGER(mpi) :: isize
13283 INTEGER(mpi) :: ivoff
13284 INTEGER(mpi) :: itgbi
13286 INTEGER(mpi) :: msize
13287 INTEGER(mpi),
PARAMETER :: mxsize = 1000
13289 INTEGER(mpl):: length
13291 REAL(mpd),
DIMENSION(:),
ALLOCATABLE :: auxVectorD
13292 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: auxVectorI
13293 REAL(mpd),
DIMENSION(:),
ALLOCATABLE :: resParGroup
13294 REAL(mpd),
DIMENSION(:),
ALLOCATABLE :: blockParGroup
13302 IF (isize <= mxsize)
THEN
13303 msize=max(msize,isize)
13305 print *,
' CKPGRP: par. group', ipgrp,
' not checked -- too large: ', isize
13308 IF (msize == 0)
RETURN
13311 length=int(msize,mpl)*(int(msize,mpl)+1)/2
13312 CALL mpalloc(blockpargroup,length,
'(matrix) block for parameter groups (D)')
13314 CALL mpalloc(respargroup,length,
'residuals for parameter groups (D)')
13315 CALL mpalloc(auxvectori,length,
'auxiliary array (I)')
13316 CALL mpalloc(auxvectord,length,
'auxiliary array (D)')
13320 print *,
' CKPGRP par. group first label size rank'
13323 IF (isize > mxsize) cycle
13330 blockpargroup(ij)=matij(ivoff+i,ivoff+j)
13334 CALL sqminv(blockpargroup,respargroup,isize,irank, auxvectord, auxvectori)
13337 IF (isize == irank)
THEN
13341 print *,
' CKPGRP ', ipgrp,
globalparlabelindex(1,itgbi), isize, irank,
' rank deficit !!!'
13359 INTEGER(mpl) :: nan
13360 INTEGER(mpl) :: neg
13362 print *,
' Checking global matrix(D) for NANs ',
size(
globalmatd,kind=mpl)
13367 print *,
' i, nan ', i, nan
13373 print *,
' Checking diagonal elements ',
nagb
13378 print *,
' i, neg ', i, neg
13382 print *,
' CHKMAT summary ', nan, neg
13404 REAL(mpd),
INTENT(IN) :: chi2
13405 INTEGER(mpi),
INTENT(IN) :: ithrd
13406 INTEGER(mpi),
INTENT(IN) :: ndf
13407 REAL(mpd),
INTENT(IN) :: dw
13409 INTEGER(mpl) ::nadd
13437 REAL(mpd),
INTENT(OUT) ::chi2
13438 INTEGER(mpl),
INTENT(OUT) ::ndf
13439 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 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 sqmibb(v, b, n, nbdr, nbnd, inv, nrank, vbnd, vbdr, aux, vbk, vzru, scdiag, scflag, evdmin, evdmax)
Bordered band matrix.
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 devinv(n, diag, u, v)
Inversion by diagonalization.
subroutine sqmibb2(v, b, n, nbdr, nbnd, inv, nrank, vbnd, vbdr, aux, vbk, vzru, scdiag, scflag, evdmin, evdmax)
Band bordered matrix.
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) nzgb
number of zero global derivatives read
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(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(6) nrejec
rejected records
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(mpl), dimension(:), allocatable globalparlabelzeros
global parameters label with zero derivative counters
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(mps) cndlmx
cut on log10(condition of band part) of local (bordered-band matrix) fit
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 prtrej(lun)
Print rejection statistics.
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