906#ifdef SCOREP_USER_ENABLE
907#include "scorep/SCOREP_User.inc"
932 REAL,
DIMENSION(2) :: ta
935 INTEGER(mpi) :: iopnmp
946 INTEGER(mpi) :: nsecnd
947 INTEGER(mpi) :: ntsec
949 CHARACTER (LEN=24) :: chdate
950 CHARACTER (LEN=24) :: chost
952 CHARACTER (LEN=6) :: c6
953 INTEGER major, minor, patch
977 CALL getenv(
'HOSTNAME',chost)
978 IF (chost(1:1) ==
' ')
CALL getenv(
'HOST',chost)
979 WRITE(*,*)
'($Id: dd0c569a1aafb6f6eb2f26b9b9537a685639ca25 $)'
984 CALL ilaver( major,minor, patch )
985 WRITE(*,110) lapack64, major,minor, patch
986110
FORMAT(
' using LAPACK64 with ',(a),
', version ',i0,
'.',i0,
'.',i0)
988 WRITE(*,*)
'using Intel oneMKL PARDISO'
992 WRITE(*,111) __gnuc__ , __gnuc_minor__ , __gnuc_patchlevel__
993111
FORMAT(
' compiled with gcc ',i0,
'.',i0,
'.',i0)
996 WRITE(*,111) __pgic__ , __pgic_minor__ , __pgic_patchlevel__
997111
FORMAT(
' compiled with pgi ',i0,
'.',i0,
'.',i0)
999#ifdef SCOREP_USER_ENABLE
1000 WRITE(*,*)
'instrumenting Score-P user regions'
1003 WRITE(*,*)
' < Millepede II-P starting ... ',chdate
1004 WRITE(*,*)
' ',chost
1007 WRITE(8,*)
'($Id: dd0c569a1aafb6f6eb2f26b9b9537a685639ca25 $)'
1009 WRITE(8,*)
'Log-file Millepede II-P ', chdate
1010 WRITE(8,*)
' ', chost
1012 CALL peend(-1,
'Still running or crashed')
1021 WRITE(*,*)
'!!! Checking input only, no calculation of a solution !!!'
1022 WRITE(8,*)
'!!! Checking input only, no calculation of a solution !!!'
1041 CALL getenv(
'OMP_NUM_THREADS',c6)
1043 CALL getenv(lapack64//
'_NUM_THREADS',c6)
1045 IF (c6(1:1) ==
' ')
THEN
1047 WRITE(*,*)
'Number of threads for LAPACK: unkown (empty OMP_NUM_THREADS)'
1049 WRITE(*,*)
'Number of threads for LAPACK: unkown (empty ',lapack64//
'_NUM_THREADS)'
1052 WRITE(*,*)
'Number of threads for LAPACK: ', c6
1072 CALL mvopen(lun,
'millepede.his')
1078 CALL mvopen(1,
'mpdebug.txt')
1082 times(0)=rstext-rstp
1088 times(1)=rloop1-rstext
1092 WRITE(8,*)
'Chi square cut equiv 3 st.dev applied ...'
1093 WRITE(8,*)
' in first iteration with factor',
chicut
1094 WRITE(8,*)
' in second iteration with factor',
chirem
1095 WRITE(8,*)
' (reduced by sqrt in next iterations)'
1099 WRITE(8,*)
'Down-weighting of outliers in',
lhuber,
' iterations'
1100 WRITE(8,*)
'Cut on downweight fraction',
dwcut
1104 times(2)=rloop2-rloop1
1109 CALL peend(5,
'Ended without solution (empty constraints)')
1111 CALL peend(0,
'Ended normally')
1148 CALL gmpdef(6,1,
'log10(#records) vs file number')
1149 CALL gmpdef(7,1,
'final rejection fraction vs file number')
1151 'final <Chi^2/Ndf> from accepted local fits vs file number')
1152 CALL gmpdef(9,1,
'<Ndf> from accepted local fits vs file number')
1158 rej=real(nrc-
jfd(kfl),mps)/real(nrc,mps)
1159 CALL gmpxy(6,real(kfl,mps),log10(real(nrc,mps)))
1160 CALL gmpxy(7,real(kfl,mps),rej)
1162 IF (
jfd(kfl) > 0)
THEN
1163 c2ndf=
cfd(kfl)/real(
jfd(kfl),mps)
1164 CALL gmpxy(8,real(kfl,mps),c2ndf)
1165 andf=real(
dfd(kfl),mps)/real(
jfd(kfl),mps)
1166 CALL gmpxy(9,real(kfl,mps),andf)
1183 WRITE(*,*)
'Misalignment test wire chamber'
1186 CALL hmpdef( 9,-0.0015,+0.0015,
'True - fitted displacement')
1187 CALL hmpdef(10,-0.0015,+0.0015,
'True - fitted Vdrift')
1193 sums(1)=sums(1)+diff
1194 sums(2)=sums(2)+diff*diff
1196 sums(3)=sums(3)+diff
1197 sums(4)=sums(4)+diff*diff
1199 sums(1)=0.01_mpd*sums(1)
1200 sums(2)=sqrt(0.01_mpd*sums(2))
1201 sums(3)=0.01_mpd*sums(3)
1202 sums(4)=sqrt(0.01_mpd*sums(4))
1203 WRITE(*,143)
'Parameters 1 - 100: mean =',sums(1),
'rms =',sums(2)
1204 WRITE(*,143)
'Parameters 101 - 200: mean =',sums(3),
'rms =',sums(4)
1205143
FORMAT(6x,a28,f9.6,3x,a5,f9.6)
1208 WRITE(*,*)
' I label simulated fitted diff'
1209 WRITE(*,*)
' -------------------------------------------- '
1229 WRITE(*,*)
'Misalignment test Si tracker'
1232 CALL hmpdef( 9,-0.0025,+0.0025,
'True - fitted displacement X')
1233 CALL hmpdef(10,-0.025,+0.025,
'True - fitted displacement Y')
1244 sums(1)=sums(1)+1.0_mpd
1245 sums(2)=sums(2)+diff
1246 sums(3)=sums(3)+diff*diff
1251 err=sqrt(abs(gmati))
1253 sums(7)=sums(7)+1.0_mpd
1254 sums(8)=sums(8)+diff
1255 sums(9)=sums(9)+diff*diff
1258 IF (mod(i,3) == 1)
THEN
1262 sums(4)=sums(4)+1.0_mpd
1263 sums(5)=sums(5)+diff
1264 sums(6)=sums(6)+diff*diff
1269 err=sqrt(abs(gmati))
1271 sums(7)=sums(7)+1.0_mpd
1272 sums(8)=sums(8)+diff
1273 sums(9)=sums(9)+diff*diff
1278 sums(2)=sums(2)/sums(1)
1279 sums(3)=sqrt(sums(3)/sums(1))
1280 sums(5)=sums(5)/sums(4)
1281 sums(6)=sqrt(sums(6)/sums(4))
1282 WRITE(*,143)
'Parameters 1 - 500: mean =',sums(2),
'rms =',sums(3)
1283 WRITE(*,143)
'Parameters 501 - 700: mean =',sums(5),
'rms =',sums(6)
1284 IF (sums(7) > 0.5_mpd)
THEN
1285 sums(8)=sums(8)/sums(7)
1286 sums(9)=sqrt(sums(9)/sums(7))
1287 WRITE(*,143)
'Parameter pulls, all: mean =',sums(8),
'rms =',sums(9)
1291 WRITE(*,*)
' I label simulated fitted diff'
1292 WRITE(*,*)
' -------------------------------------------- '
1304 IF (mod(i,3) == 1)
THEN
1324 WRITE(8,*)
'Record',
nrec1,
' has largest residual:',
value1
1327 WRITE(8,*)
'Record',
nrec2,
' has largest Chi^2/Ndf:',
value2
1331 WRITE(8,*)
'Record',
nrec3,
' is first with error (rank deficit/NaN)'
1335 WRITE(8,*)
'In total 3 +',
nloopn,
' loops through the data files'
1337 WRITE(8,*)
'In total 2 +',
nloopn,
' loops through the data files'
1341 WRITE(8,*)
'In total ',
mnrsit,
' internal MINRES iterations'
1349 ntsec=nint(deltat,mpi)
1350 CALL sechms(deltat,nhour,minut,secnd)
1351 nsecnd=nint(secnd,mpi)
1352 WRITE(8,*)
'Total time =',ntsec,
' seconds =',nhour,
' h',minut, &
1353 ' m',nsecnd,
' seconds'
1355 WRITE(8,*)
'end ', chdate
1362 IF(sum(
nrejec) /= 0)
THEN
1364 WRITE(8,*)
'Data records rejected in last iteration: '
1371 WRITE(*,*)
' < Millepede II-P ending ... ', chdate
1378106
FORMAT(
' PARDISO dyn. memory allocation: ',f11.6,
' GB')
1388 WRITE(*,*)
'Postprocessing:'
1398102
FORMAT(2x,i4,i10,2x,3f10.5)
1399103
FORMAT(
' Times [in sec] for text processing',f12.3/ &
1402 ' func. value ',f12.3,
' *',f4.0/ &
1403 ' func. value, global matrix, solution',f12.3,
' *',f4.0/ &
1404 ' new solution',f12.3,
' *',f4.0/)
1405105
FORMAT(
' Peak dynamic memory allocation: ',f11.6,
' GB')
1425 INTEGER(mpi) :: istop
1426 INTEGER(mpi) :: itgbi
1427 INTEGER(mpi) :: itgbl
1429 INTEGER(mpi) :: itnlim
1430 INTEGER(mpi) :: nout
1432 INTEGER(mpi),
INTENT(IN) :: ivgbi
1469 globalcorrections, itnlim, nout, rtol, istop, itn, anorm, acond, rnorm, arnorm, ynorm)
1473 globalcorrections, itnlim, nout, rtol, istop, itn, anorm, acond, rnorm, arnorm, ynorm)
1476 globalcorrections, itnlim, nout, rtol, istop, itn, anorm, acond, rnorm, arnorm, ynorm)
1482 err=sqrt(abs(real(gmati,mps)))
1483 IF(gmati < 0.0_mpd) err=-err
1484 diag=matij(ivgbi,ivgbi)
1485 gcor2=real(1.0_mpd-1.0_mpd/(gmati*diag),mps)
1487101
FORMAT(1x,
' label parameter presigma differ', &
1488 ' Error gcor^2 iit'/ 1x,
'---------',2x,5(
'-----------'),2x,
'----')
1489102
FORMAT(i10,2x,4g12.4,f7.4,i6,i4)
1509 INTEGER(mpi) :: istop
1510 INTEGER(mpi) :: itgbi
1511 INTEGER(mpi) :: itgbl
1513 INTEGER(mpi) :: itnlim
1514 INTEGER(mpi) :: nout
1516 INTEGER(mpi),
INTENT(IN) :: ivgbi
1545 mxxnrm = real(
nagb,mpd)/sqrt(epsilon(mxxnrm))
1547 trcond = 1.0_mpd/epsilon(trcond)
1548 ELSE IF(
mrmode == 2)
THEN
1556 itnlim=itnlim, rtol=rtol, maxxnorm=mxxnrm, trancond=trcond, &
1560 itnlim=itnlim, rtol=rtol, maxxnorm=mxxnrm, trancond=trcond, &
1564 itnlim=itnlim, rtol=rtol, maxxnorm=mxxnrm, trancond=trcond, &
1571 err=sqrt(abs(real(gmati,mps)))
1572 IF(gmati < 0.0_mpd) err=-err
1573 diag=matij(ivgbi,ivgbi)
1574 gcor2=real(1.0_mpd-1.0_mpd/(gmati*diag),mps)
1576101
FORMAT(1x,
' label parameter presigma differ', &
1577 ' Error gcor^2 iit'/ 1x,
'---------',2x,5(
'-----------'),2x,
'----')
1578102
FORMAT(i10,2x,4g12.4,f7.4,i6,i4)
1591 INTEGER(mpi) :: icgb
1592 INTEGER(mpi) :: irhs
1593 INTEGER(mpi) :: itgbi
1594 INTEGER(mpi) :: ivgb
1596 INTEGER(mpi) :: jcgb
1598 INTEGER(mpi) :: label
1600 INTEGER(mpi) :: inone
1603 REAL(mpd) :: drhs(4)
1604 INTEGER(mpi) :: idrh (4)
1629 IF(abs(rhs) > climit)
THEN
1635 WRITE(*,101) (idrh(l),drhs(l),l=1,irhs)
1644 WRITE(*,101) (idrh(l),drhs(l),l=1,irhs)
1647 WRITE(*,102)
' Constraints: only equation values >', climit,
' are printed'
1648101
FORMAT(
' ',4(i6,g11.3))
1649102
FORMAT(a,g11.2,a)
1662 INTEGER(mpi) :: icgb
1663 INTEGER(mpi) :: icgrp
1664 INTEGER(mpi) :: ioff
1665 INTEGER(mpi) :: itgbi
1667 INTEGER(mpi) :: jcgb
1668 INTEGER(mpi) :: label
1669 INTEGER(mpi) :: labelf
1670 INTEGER(mpi) :: labell
1671 INTEGER(mpi) :: last
1672 INTEGER(mpi) :: line1
1673 INTEGER(mpi) :: ncon
1674 INTEGER(mpi) :: ndiff
1675 INTEGER(mpi) :: npar
1676 INTEGER(mpi) :: inone
1677 INTEGER(mpi) :: itype
1678 INTEGER(mpi) :: ncgbd
1679 INTEGER(mpi) :: ncgbr
1680 INTEGER(mpi) :: ncgbw
1681 INTEGER(mpi) :: ncgrpd
1682 INTEGER(mpi) :: ncgrpr
1683 INTEGER(mpi) :: next
1685 INTEGER(mpl):: length
1686 INTEGER(mpl) :: rows
1688 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: vecParConsOffsets
1689 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: vecParConsList
1690 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: vecConsParOffsets
1691 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: vecConsParList
1692 INTEGER(mpi),
DIMENSION(:,:),
ALLOCATABLE :: matConsGroupIndex
1705 IF(last < 0.AND.label < 0)
THEN
1708 IF(itype == 2) ncgbw=ncgbw+1
1715 IF(label > 0.AND.itype == 2)
THEN
1722 IF (ncgbw == 0)
THEN
1723 WRITE(*,*)
'GRPCON:',
ncgb,
' constraints found in steering files'
1725 WRITE(*,*)
'GRPCON:',
ncgb,
' constraints found in steering files,',ncgbw,
'weighted'
1730 length=
ncgb+1; rows=3
1743 CALL mpalloc(matconsgroupindex,rows,length,
'group index for constraint (I)')
1746 CALL mpalloc(vecconsparoffsets,length,
'offsets for global par list for cons. (I)')
1748 CALL mpalloc(vecparconsoffsets,length,
'offsets for cons. list for global par. (I)')
1749 vecparconsoffsets(1)=0
1751 vecparconsoffsets(i+1)=vecparconsoffsets(i)+
globalparcons(i)
1755 length=vecparconsoffsets(
ntgb+1)
1756 CALL mpalloc(vecconsparlist,length,
'global par. list for constraint (I)')
1757 CALL mpalloc(vecparconslist,length,
'constraint list for global par. (I)')
1762 vecconsparoffsets(1)=ioff
1774 vecparconslist(vecparconsoffsets(itgbi)+
globalparcons(itgbi))=icgb
1776 vecconsparlist(ioff+npar)=itgbi
1782 CALL sort1k(vecconsparlist(ioff+1),npar)
1786 next=vecconsparlist(ioff+j)
1787 IF (next /= last)
THEN
1789 vecconsparlist(ioff+ndiff) = next
1798 vecconsparoffsets(icgb+1)=ioff
1811 print *,
' Constraint #parameters first par. last par. first line'
1819 npar=vecconsparoffsets(icgb+1)-vecconsparoffsets(icgb)
1823 print *, jcgb, npar, labelf, labell, line1
1826 icgrp=matconsgroupindex(1,icgb)
1827 IF (icgrp == 0)
THEN
1829 DO i=vecconsparoffsets(icgb)+1, vecconsparoffsets(icgb+1)
1830 itgbi=vecconsparlist(i)
1832 DO j=vecparconsoffsets(itgbi)+1,vecparconsoffsets(itgbi+1)
1833 icgrp=matconsgroupindex(1,vecparconslist(j))
1839 IF (icgrp == 0)
THEN
1846 matconsgroupindex(2,icgb)=jcgb
1847 matconsgroupindex(3,icgb)=icgb
1848 DO i=vecconsparoffsets(icgb)+1, vecconsparoffsets(icgb+1)
1849 itgbi=vecconsparlist(i)
1852 DO j=vecparconsoffsets(itgbi)+1,vecparconsoffsets(itgbi+1)
1853 matconsgroupindex(1,vecparconslist(j))=icgrp
1857 WRITE(*,*)
'GRPCON:',
ncgrp,
' disjoint constraints groups built'
1865 icgb=matconsgroupindex(3,jcgb)
1870 icgrp=matconsgroupindex(1,jcgb)
1890 print *,
' cons.group first con. first par. last par. #cons #par'
1901 print *, icgrp,
matconsgroups(1,icgrp), labelf, labell, ncon, npar
1904 IF (ncon == npar)
THEN
1911 print *, icgrp,
matconsgroups(1,icgrp), labelf, labell,
' : cons.group resolved'
1926 print *, icgrp,
matconsgroups(1,icgrp), labelf, labell,
' : cons.group redundant'
1931 IF (ncgrpr > 0)
THEN
1932 WRITE(*,*)
'GRPCON:',ncgbr,
' redundancy constraints in ', ncgrpr,
' groups resolved'
1936 IF (ncgrpd > 0)
THEN
1937 WRITE(*,*)
'GRPCON:',ncgbd,
' redundancy constraints in ', ncgrpd,
' groups detected'
1960 INTEGER(mpi) :: icgb
1961 INTEGER(mpi) :: icgrp
1962 INTEGER(mpi) :: ifrst
1963 INTEGER(mpi) :: ilast
1964 INTEGER(mpi) :: isblck
1965 INTEGER(mpi) :: itgbi
1966 INTEGER(mpi) :: ivgb
1968 INTEGER(mpi) :: jcgb
1969 INTEGER(mpi) :: jfrst
1970 INTEGER(mpi) :: label
1971 INTEGER(mpi) :: labelf
1972 INTEGER(mpi) :: labell
1973 INTEGER(mpi) :: ncon
1974 INTEGER(mpi) :: ngrp
1975 INTEGER(mpi) :: npar
1976 INTEGER(mpi) :: ncnmxb
1977 INTEGER(mpi) :: ncnmxg
1978 INTEGER(mpi) :: nprmxb
1979 INTEGER(mpi) :: nprmxg
1980 INTEGER(mpi) :: inone
1981 INTEGER(mpi) :: nvar
1983 INTEGER(mpl):: length
1984 INTEGER(mpl) :: rows
1986 INTEGER(mpi),
DIMENSION(:,:),
ALLOCATABLE :: matConsGroupIndex
1999 length=
ncgrp+1; rows=3
2004 CALL mpalloc(matconsgroupindex,rows,length,
'group index for constraint (I)')
2041 IF (nvar > 0 .OR.
iskpec == 0)
THEN
2044 matconsgroupindex(1,
ncgb)=ngrp+1
2045 matconsgroupindex(2,
ncgb)=icgb
2046 matconsgroupindex(3,
ncgb)=nvar
2049 IF (
ncgb > ncon) ngrp=ngrp+1
2055 WRITE(*,*)
'PRPCON:',
ncgbe,
' empty constraints skipped'
2057 WRITE(*,*)
'PRPCON:',
ncgbe,
' empty constraints detected, to be fixed !!!'
2058 WRITE(*,*)
' (use option "skipemptycons" to skip those)'
2063 WRITE(*,*)
'!!! Switch to "-C" (checking input only), no calculation of a solution !!!'
2064 WRITE(8,*)
'!!! Switch to "-C" (checking input only), no calculation of a solution !!!'
2069 WRITE(*,*)
'PRPCON:',
ncgb,
' constraints accepted'
2072 IF(
ncgb == 0)
RETURN
2079 icgb=matconsgroupindex(2,jcgb)
2084 icgrp=matconsgroupindex(1,jcgb)
2110 WRITE(*,*)
' Cons. sorted index #var.par. first line first label last label'
2111 WRITE(*,*)
' Cons. group index first cons. last cons. first label last label'
2112 WRITE(*,*)
' Cons. block index first group last group first label last label'
2118 nvar=matconsgroupindex(3,jcgb)
2122 WRITE(*,*)
' Cons. sorted', jcgb, nvar, &
2125 WRITE(*,*)
' Cons. sorted', jcgb,
' empty (0)', &
2134 WRITE(*,*)
' Cons. group ', icgrp,
matconsgroups(1,icgrp), &
2147 DO i=icgrp,isblck,-1
2161 WRITE(*,*)
' Cons. block ',
ncblck, isblck, icgrp, labelf, labell
2179 ncnmxb=max(ncnmxb,ncon)
2180 nprmxb=max(nprmxb,npar)
2205 ncnmxg=max(ncnmxg,ncon)
2206 nprmxg=max(nprmxg,npar)
2241 WRITE(*,*)
'PRPCON: constraints split into ',
ncgrp,
'(disjoint) groups,'
2242 WRITE(*,*)
' groups combined into ',
ncblck,
'(non overlapping) blocks'
2243 WRITE(*,*)
' max group size (cons., par.) ', ncnmxg, nprmxg
2244 WRITE(*,*)
' max block size (cons., par.) ', ncnmxb, nprmxb
2262 INTEGER(mpi) :: icgb
2263 INTEGER(mpi) :: icgrp
2265 INTEGER(mpi) :: ifirst
2266 INTEGER(mpi) :: ilast
2267 INTEGER(mpl) :: ioffc
2268 INTEGER(mpl) :: ioffp
2269 INTEGER(mpi) :: irank
2270 INTEGER(mpi) :: ipar0
2271 INTEGER(mpi) :: itgbi
2272 INTEGER(mpi) :: ivgb
2274 INTEGER(mpi) :: jcgb
2276 INTEGER(mpi) :: label
2277 INTEGER(mpi) :: ncon
2278 INTEGER(mpi) :: npar
2279 INTEGER(mpi) :: nrank
2280 INTEGER(mpi) :: inone
2285 INTEGER(mpl):: length
2286 REAL(mpd),
DIMENSION(:),
ALLOCATABLE :: matConstraintsT
2287 REAL(mpd),
DIMENSION(:),
ALLOCATABLE :: auxVectorD
2288 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: auxVectorI
2292 IF(
ncgb == 0)
RETURN
2301 CALL mpalloc(auxvectori,length,
'auxiliary array (I)')
2302 CALL mpalloc(auxvectord,length,
'auxiliary array (D)')
2305 CALL mpalloc(matconstraintst,length,
'transposed matrix of constraints (blocks)')
2306 matconstraintst=0.0_mpd
2319 WRITE(*,*)
' Constraint group, #con, rank', icgrp, ncon, 0,
' (empty)'
2322 DO jcgb=ifirst,ilast
2334 IF(ivgb > 0) matconstraintst(int(jcgb-ifirst,mpl)*int(npar,mpl)+ivgb-ipar0+ioffc)= &
2335 matconstraintst(int(jcgb-ifirst,mpl)*int(npar,mpl)+ivgb-ipar0+ioffc)+factr
2342 DO ll=ioffc+1,ioffc+npar
2348 matconstraintst(int(i-1,mpl)*int(npar,mpl)+ll)* &
2349 matconstraintst(int(j-1,mpl)*int(npar,mpl)+ll)
2355 IF (
icheck > 1 .OR. irank < ncon)
THEN
2356 WRITE(*,*)
' Constraint group, #con, rank', icgrp, ncon, irank
2357 IF (irank < ncon)
THEN
2358 WRITE(*,*)
' .. rank deficit !! '
2359 WRITE(*,*)
' E.g. fix all parameters and remove all constraints related to label ', &
2364 ioffc=ioffc+int(npar,mpl)*int(ncon,mpl)
2371 WRITE(*,*)
'Rank of product matrix of constraints is',nrank, &
2372 ' for',
ncgb,
' constraint equations'
2373 WRITE(8,*)
'Rank of product matrix of constraints is',nrank, &
2374 ' for',
ncgb,
' constraint equations'
2375 IF(nrank <
ncgb)
THEN
2376 WRITE(*,*)
'Warning: insufficient constraint equations!'
2377 WRITE(8,*)
'Warning: insufficient constraint equations!'
2380 WRITE(*,*)
' --> enforcing SUBITO mode'
2381 WRITE(8,*)
' --> enforcing SUBITO mode'
2388 print *,
'QL decomposition of constraints matrix'
2391 WRITE(
lunlog,*)
'QL decomposition of constraints matrix'
2403 CALL lpqldec(matconstraintst,evmin,evmax)
2407 print *,
' largest |eigenvalue| of L: ', evmax
2408 print *,
' smallest |eigenvalue| of L: ', evmin
2409 IF (evmin == 0.0_mpd.AND.
icheck == 0)
THEN
2410 CALL peend(27,
'Aborted, singular QL decomposition of constraints matrix')
2411 stop
'FEASMA: stopping due to singular QL decomposition of constraints matrix'
2437 INTEGER(mpi) :: icgb
2438 INTEGER(mpi) :: icgrp
2439 INTEGER(mpi) :: iter
2440 INTEGER(mpi) :: itgbi
2441 INTEGER(mpi) :: ivgb
2442 INTEGER(mpi) :: ieblck
2443 INTEGER(mpi) :: isblck
2444 INTEGER(mpi) :: ifirst
2445 INTEGER(mpi) :: ilast
2447 INTEGER(mpi) :: jcgb
2448 INTEGER(mpi) :: label
2449 INTEGER(mpi) :: inone
2450 INTEGER(mpi) :: ncon
2452 REAL(mps),
INTENT(IN) :: concut
2453 INTEGER(mpi),
INTENT(OUT) :: iact
2460 REAL(mpd),
DIMENSION(:),
ALLOCATABLE :: vecCorrections
2464 IF(
ncgb == 0)
RETURN
2494 sum1=sqrt(sum1/real(
ncgb,mpd))
2495 sum2=sum2/real(
ncgb,mpd)
2497 IF(iter == 1.AND.sum1 < concut)
RETURN
2499 IF(iter == 1.AND.
ncgb <= 12)
THEN
2501 WRITE(*,*)
'Constraint equation discrepancies:'
2503101
FORMAT(4x,4(i5,g12.4))
2505103
FORMAT(10x,
' Cut on rms value is',g8.1)
2510 WRITE(*,*)
'Improve constraints'
2514 WRITE(*,102) iter,sum1,sum2,sum3
2515102
FORMAT(i6,
' rms',g12.4,
' avrg_abs',g12.4,
' max_abs',g12.4)
2517 CALL mpalloc(veccorrections,int(
nvgb,mpl),
'constraint corrections')
2518 veccorrections=0.0_mpd
2526 ieblck=isblck+(ncon*(ncon+1))/2
2595 INTEGER(mpi) :: iact
2596 INTEGER(mpi) :: ierrc
2597 INTEGER(mpi) :: ierrf
2598 INTEGER(mpi) :: ioffp
2600 INTEGER(mpi) :: ithr
2601 INTEGER(mpi) :: jfile
2602 INTEGER(mpi) :: jrec
2604 INTEGER(mpi) :: kfile
2607 INTEGER(mpi) :: mpri
2609 INTEGER(mpi) :: nact
2610 INTEGER(mpi) :: nbuf
2611 INTEGER(mpi) :: ndata
2612 INTEGER(mpi) :: noff
2613 INTEGER(mpi) :: noffs
2614 INTEGER(mpi) :: npointer
2615 INTEGER(mpi) :: npri
2619 INTEGER(mpi) :: nrpr
2620 INTEGER(mpi) :: nthr
2621 INTEGER(mpi) :: ntot
2622 INTEGER(mpi) :: maxRecordSize
2623 INTEGER(mpi) :: maxRecordFile
2625 INTEGER(mpi),
INTENT(OUT) :: more
2635 CHARACTER (LEN=7) :: cfile
2640 SUBROUTINE readc(bufferD, bufferF, bufferI, bufferLength, lun, err)
BIND(c)
2642 REAL(c_double),
DIMENSION(*),
INTENT(OUT) :: bufferD
2643 REAL(c_float),
DIMENSION(*),
INTENT(OUT) :: bufferF
2644 INTEGER(c_int),
DIMENSION(*),
INTENT(OUT) :: bufferI
2645 INTEGER(c_int),
INTENT(INOUT) :: bufferLength
2646 INTEGER(c_int),
INTENT(IN),
VALUE :: lun
2647 INTEGER(c_int),
INTENT(OUT) :: err
2648 END SUBROUTINE readc
2654 DATA npri / 0 /, mpri / 1000 /
2702 files:
DO WHILE (jfile > 0)
2706 CALL binopn(kfile,ithr,ios)
2712 IF(kfile <=
nfilf)
THEN
2735 eof=(ierrc <= 0.AND.ierrc /= -4)
2736 IF(eof.AND.ierrc < 0)
THEN
2737 WRITE(*,*)
'Read error for binary Cfile', kfile,
'record', jrec+1,
':', ierrc
2738 WRITE(8,*)
'Read error for binary Cfile', kfile,
'record', jrec+1,
':', ierrc
2740 WRITE(cfile,
'(I7)') kfile
2741 CALL peend(18,
'Aborted, read error(s) for binary file ' // cfile)
2742 stop
'PEREAD: stopping due to read errors (bad record, wrong file type?)'
2744 IF (
kfd(1,jfile) == 1)
THEN
2750 IF(eof)
EXIT records
2755 xfd(jfile)=max(
xfd(jfile),n)
2772 IF ((noff+nr+2+
ndimbuf >= ndata*(iact+1)).OR.(nbuf >= npointer))
EXIT files
2787 IF (
kfd(1,jfile) == 1)
THEN
2788 print *,
'PEREAD: file ', kfile,
'read the first time, found',jrec,
' records'
2792 IF (-
kfd(1,jfile) /= jrec)
THEN
2793 WRITE(cfile,
'(I7)') kfile
2794 CALL peend(19,
'Aborted, binary file modified (length) ' // cfile)
2795 stop
'PEREAD: file modified (length)'
2846 IF (
kfd(1,jfile) == 1)
kfd(1,jfile)=-nrc
2865 DO WHILE (
nloopn == 0.AND.ntot >= nrpr)
2866 WRITE(*,*)
' Record ',nrpr
2867 IF (nrpr < 100000)
THEN
2876 IF (npri == 1)
WRITE(*,100)
2878100
FORMAT(/
' PeRead records active file' &
2879 /
' total block threads number')
2880101
FORMAT(
' PeRead',4i10)
2893 IF (
xfd(k) > maxrecordsize)
THEN
2894 maxrecordsize=
xfd(k)
2897 dw=real(-
kfd(1,k),mpd)
2900 ds1=ds1+dw*real(
wfd(k),mpd)
2901 ds2=ds2+dw*real(
wfd(k)**2,mpd)
2903 print *,
'PEREAD: file ', maxrecordfile,
'with max record size ', maxrecordsize
2904 IF (
nfilw > 0.AND.ds0 > 0.0_mpd)
THEN
2908 WRITE(lun,177)
nfilw,real(ds1,mps),real(ds2,mps)
2909177
FORMAT(/
' !!!!!',i4,
' weighted binary files', &
2910 /
' !!!!! mean, variance of weights =',2g12.4)
2928179
FORMAT(/
' Read cache usage (#blocks, #records, ', &
2929 'min,max records/block'/17x,i10,i12,2i10)
2947 INTEGER(mpi),
INTENT(IN) :: mode
2949 INTEGER(mpi) :: ibuf
2950 INTEGER(mpi) :: ichunk
2952 INTEGER(mpi) :: itgbi
2958 INTEGER(mpi),
PARAMETER :: maxbad = 100
2959 INTEGER(mpi) :: nbad
2960 INTEGER(mpi) :: nerr
2961 INTEGER(mpi) :: inone
2979 CALL isjajb(nst,ist,ja,jb,jsp)
2998#ifdef SCOREP_USER_ENABLE
2999 scorep_user_region_by_name_begin(
"UR_peprep", scorep_user_region_type_common)
3004 CALL pechk(ibuf,nerr)
3007 IF(nbad >= maxbad)
EXIT
3012 CALL isjajb(nst,ist,ja,jb,jsp)
3025 CALL peend(20,
'Aborted, bad binary records')
3026 stop
'PEREAD: stopping due to bad records'
3029#ifdef SCOREP_USER_ENABLE
3030 scorep_user_region_by_name_end(
"UR_peprep")
3050 INTEGER(mpi) :: ioff
3057 INTEGER(mpi),
INTENT(IN) :: ibuf
3058 INTEGER(mpi),
INTENT(OUT) :: nerr
3067 outer:
DO WHILE(is < nst)
3072 IF(is > nst)
EXIT outer
3078 IF(is > nst)
EXIT outer
3095100
FORMAT(
' PEREAD: record ', i8,
' in file ',i6,
' is broken !!!')
3105101
FORMAT(
' PEREAD: record ', i8,
' in file ',i6,
' contains ', i6,
' NaNs !!!')
3121 INTEGER(mpi) :: ibuf
3122 INTEGER(mpi) :: ichunk
3123 INTEGER(mpi) :: iproc
3124 INTEGER(mpi) :: ioff
3125 INTEGER(mpi) :: ioffbi
3127 INTEGER(mpi) :: itgbi
3132 INTEGER(mpi) :: nalg
3133 INTEGER(mpi) :: neqna
3136 INTEGER(mpi) :: nzero
3137 INTEGER(mpi) :: inone
3138 INTEGER(mpl) :: length
3173 CALL isjajb(nst,ist,ja,jb,jsp)
3204 CALL isjajb(nst,ist,ja,jb,jsp)
3232#ifdef SCOREP_USER_ENABLE
3233 scorep_user_region_by_name_begin(
"UR_pepgrp", scorep_user_region_type_common)
3242 CALL pargrp(ist+1,ist+nnz)
3254#ifdef SCOREP_USER_ENABLE
3255 scorep_user_region_by_name_end(
"UR_pepgrp")
3274 INTEGER(mpi) :: istart
3275 INTEGER(mpi) :: itgbi
3277 INTEGER(mpi) :: jstart
3278 INTEGER(mpi) :: jtgbi
3279 INTEGER(mpi) :: lstart
3280 INTEGER(mpi) :: ltgbi
3282 INTEGER(mpi),
INTENT(IN) :: inds
3283 INTEGER(mpi),
INTENT(IN) :: inde
3285 IF (inds > inde)
RETURN
3294 IF (istart == 0)
THEN
3295 IF (itgbi /= ltgbi+1)
THEN
3298 IF (lstart == 0)
THEN
3313 IF (istart /= jstart)
THEN
3324 IF (itgbi /= ltgbi+1)
THEN
3339 IF (istart /= jstart)
THEN
3391 INTEGER(mpi),
INTENT(IN) :: nst
3392 INTEGER(mpi),
INTENT(IN OUT) :: is
3393 INTEGER(mpi),
INTENT(OUT) :: ja
3394 INTEGER(mpi),
INTENT(OUT) :: jb
3395 INTEGER(mpi),
INTENT(OUT) :: jsp
3403 IF(is >= nst)
RETURN
3457 INTEGER(mpi) :: ioffb
3459 INTEGER(mpi) :: itgbi
3460 INTEGER(mpi) :: itgbij
3461 INTEGER(mpi) :: itgbik
3462 INTEGER(mpi) :: ivgb
3463 INTEGER(mpi) :: ivgbij
3464 INTEGER(mpi) :: ivgbik
3467 INTEGER(mpi) :: lastit
3469 INTEGER(mpi) :: ncrit
3470 INTEGER(mpi) :: ngras
3471 INTEGER(mpi) :: nparl
3473 INTEGER(mpl) :: nrej
3474 INTEGER(mpi) :: inone
3475 INTEGER(mpi) :: ilow
3476 INTEGER(mpi) :: nlow
3477 INTEGER(mpi) :: nzero
3497 CALL gmpdef(1,4,
'Function value in iterations')
3499 CALL gmpdef(2,3,
'Number of MINRES steps vs iteration nr')
3501 CALL hmpdef( 5,0.0,0.0,
'Number of degrees of freedom')
3502 CALL hmpdef(11,0.0,0.0,
'Number of local parameters')
3503 CALL hmpdef(16,0.0,24.0,
'LOG10(cond(band part decomp.)) local fit ')
3504 CALL hmpdef(23,0.0,0.0,
'SQRT of diagonal elements without presigma')
3505 CALL hmpdef(24,0.0,0.0,
'Log10 of off-diagonal elements')
3506 CALL hmpdef(25,0.0,0.0,
'Relative individual pre-sigma')
3507 CALL hmpdef(26,0.0,0.0,
'Relative global pre-sigma')
3512 'Normalized residuals of single (global) measurement')
3514 'Normalized residuals of single (local) measurement')
3516 'Pulls of single (global) measurement')
3518 'Pulls of single (local) measurement')
3519 CALL hmpdef(4,0.0,0.0,
'Chi^2/Ndf after local fit')
3520 CALL gmpdef(4,5,
'location, dispersion (res.) vs record nr')
3521 CALL gmpdef(5,5,
'location, dispersion (pull) vs record nr')
3535 IF(
nloopn == 2)
CALL hmpdef(6,0.0,0.0,
'Down-weight fraction')
3538 IF(
iterat /= lastit)
THEN
3546 ELSE IF(
iterat >= 1)
THEN
3594 WRITE(*,111) nparl,ncrit,usei,used,peaki,peakd
3595111
FORMAT(
' Write cache usage (#flush,#overrun,<levels>,', &
3596 'peak(levels))'/2i7,
',',4(f6.1,
'%'))
3604 ELSE IF (
mbandw > 0)
THEN
3636 print *,
" ... warning ..."
3637 print *,
" global parameters with too few (< MREQENA) accepted entries: ", nlow
3641 IF(
icalcm == 1 .AND. nzero > 0)
THEN
3643 WRITE(*,*)
'Warning: the rank defect of the symmetric',
nfgb, &
3644 '-by-',
nfgb,
' matrix is ',
ndefec,
' (should be zero).'
3645 WRITE(lun,*)
'Warning: the rank defect of the symmetric',
nfgb, &
3646 '-by-',
nfgb,
' matrix is ',
ndefec,
' (should be zero).'
3649 WRITE(*,*)
' --> enforcing SUBITO mode'
3650 WRITE(lun,*)
' --> enforcing SUBITO mode'
3660 elmt=real(matij(i,i),mps)
3661 IF(elmt > 0.0)
CALL hmpent(23,1.0/sqrt(elmt))
3695 CALL addsums(1, adder, 0, 1.0_mpl)
3709 IF(sgm > 0.0) weight=1.0/sgm**2
3725 IF(itgbij /= 0)
THEN
3731 IF (factrj == 0.0_mpd) cycle
3741 IF(
icalcm == 1.AND.ivgbij > 0)
THEN
3748 IF(ivgbij > 0.AND.ivgbik > 0)
THEN
3749 CALL mupdat(ivgbij,ivgbik,weight*factrj*factrk)
3755 adder=weight*dsum**2
3756 CALL addsums(1, adder, 1, 1.0_mpl)
3774 ratae=max(-99.9,ratae)
3783 WRITE(*,*)
'Data records rejected in initial loop:'
3828 WRITE(lun,*)
' === local fits have bordered band matrix structure ==='
3829 IF (
nbndr(1) > 0 )
WRITE(lun,101)
' NBNDR',
nbndr(1),
'number of records (upper/left border)'
3830 IF (
nbndr(2) > 0 )
WRITE(lun,101)
' NBNDR',
nbndr(2),
'number of records (lower/right border)'
3831 WRITE(lun,101)
' NBDRX',
nbdrx,
'max border size'
3832 WRITE(lun,101)
' NBNDX',
nbndx,
'max band width'
3841101
FORMAT(1x,a8,
' =',i14,
' = ',a)
3856 INTEGER(mpi),
INTENT(IN) :: lunp
3861101
FORMAT(
' it fc',
' fcn_value dfcn_exp slpr costh iit st', &
3862 ' ls step cutf',1x,
'rejects hhmmss FMS')
3863102
FORMAT(
' -- --',
' ----------- -------- ---- ----- --- --', &
3864 ' -- ----- ----',1x,
'------- ------ ---')
3880 INTEGER(mpl) :: nrej
3881 INTEGER(mpi) :: nsecnd
3885 REAL(mps) :: slopes(3)
3886 REAL(mps) :: steps(3)
3887 REAL,
DIMENSION(2) :: ta
3890 INTEGER(mpi),
INTENT(IN) :: lunp
3892 CHARACTER (LEN=4):: ccalcm(4)
3893 DATA ccalcm /
' end',
' S',
' F ',
' FMS' /
3897 IF(nrej > 9999999) nrej=9999999
3901 nsecnd=nint(secnd,mpi)
3910 CALL ptlopt(nfa,ma,slopes,steps)
3911 ratae=max(-99.9,min(99.9,slopes(2)/slopes(1)))
3917103
FORMAT(i3,i3,e12.5,38x,f5.1, 1x,i7, i3,i2.2,i2.2,a4)
3918104
FORMAT(i3,i3,e12.5,1x,e8.2,f6.3,f6.3,i5,2i3,f6.3,f5.1, &
3919 1x,i7, i3,i2.2,i2.2,a4)
3920105
FORMAT(i3,i3,e12.5,1x,e8.2,12x,i5,i3,9x,f5.1, &
3921 1x,i7, i3,i2.2,i2.2,a4)
3934 INTEGER(mpi) :: minut
3936 INTEGER(mpi) :: nhour
3937 INTEGER(mpl) :: nrej
3938 INTEGER(mpi) :: nsecnd
3942 REAL(mps) :: slopes(3)
3943 REAL(mps) :: steps(3)
3944 REAL,
DIMENSION(2) :: ta
3947 INTEGER(mpi),
INTENT(IN) :: lunp
3948 CHARACTER (LEN=4):: ccalcm(4)
3949 DATA ccalcm /
' end',
' S',
' F ',
' FMS' /
3953 IF(nrej > 9999999) nrej=9999999
3957 nsecnd=nint(secnd,mpi)
3961 CALL ptlopt(nfa,ma,slopes,steps)
3962 ratae=abs(slopes(2)/slopes(1))
3967104
FORMAT(3x,i3,e12.5,9x, 35x, i7, i3,i2.2,i2.2,a4)
3968105
FORMAT(3x,i3,e12.5,9x, f6.3,14x,i3,f6.3,6x, i7, i3,i2.2,i2.2,a4)
3982 INTEGER(mpi) :: nsecnd
3985 REAL,
DIMENSION(2) :: ta
3988 INTEGER(mpi),
INTENT(IN) :: lunp
3989 CHARACTER (LEN=4):: ccalcm(4)
3990 DATA ccalcm /
' end',
' S',
' F ',
' FMS' /
3995 nsecnd=nint(secnd,mpi)
3997 WRITE(lunp,106) nhour,minut,nsecnd,ccalcm(
lcalcm)
3998106
FORMAT(69x,i3,i2.2,i2.2,a4)
4008 INTEGER(mpi) :: lunit
4010 WRITE(lunit,102)
'Explanation of iteration table'
4011 WRITE(lunit,102)
'=============================='
4012 WRITE(lunit,101)
'it', &
4013 'iteration number. Global parameters are improved for it > 0.'
4014 WRITE(lunit,102)
'First function evaluation is called iteraton 0.'
4015 WRITE(lunit,101)
'fc',
'number of function evaluations.'
4016 WRITE(lunit,101)
'fcn_value',
'value of 2 x Likelihood function (LF).'
4017 WRITE(lunit,102)
'The final value is the chi^2 value of the fit and should'
4018 WRITE(lunit,102)
'be about equal to the NDF (see below).'
4019 WRITE(lunit,101)
'dfcn_exp', &
4020 'expected reduction of the value of the Likelihood function (LF)'
4021 WRITE(lunit,101)
'slpr',
'ratio of the actual slope to inital slope.'
4022 WRITE(lunit,101)
'costh', &
4023 'cosine of angle between search direction and -gradient'
4025 WRITE(lunit,101)
'iit', &
4026 'number of internal iterations in MINRES algorithm'
4027 WRITE(lunit,101)
'st',
'stop code of MINRES algorithm'
4028 WRITE(lunit,102)
'< 0: rhs is very special, with beta2 = 0'
4029 WRITE(lunit,102)
'= 0: rhs b = 0, i.e. the exact solution is x = 0'
4030 WRITE(lunit,102)
'= 1 requested accuracy achieved, as determined by rtol'
4031 WRITE(lunit,102)
'= 2 reasonable accuracy achieved, given eps'
4032 WRITE(lunit,102)
'= 3 x has converged to an eigenvector'
4033 WRITE(lunit,102)
'= 4 matrix ill-conditioned (Acond has exceeded 0.1/eps)'
4034 WRITE(lunit,102)
'= 5 the iteration limit was reached'
4035 WRITE(lunit,102)
'= 6 Matrix x vector does not define a symmetric matrix'
4036 WRITE(lunit,102)
'= 7 Preconditioner does not define a symmetric matrix'
4037 ELSEIF (
metsol == 5)
THEN
4038 WRITE(lunit,101)
'iit', &
4039 'number of internal iterations in MINRES-QLP algorithm'
4040 WRITE(lunit,101)
'st',
'stop code of MINRES-QLP algorithm'
4041 WRITE(lunit,102)
'= 1: beta_{k+1} < eps, iteration k is the final Lanczos step.'
4042 WRITE(lunit,102)
'= 2: beta2 = 0. If M = I, b and x are eigenvectors of A.'
4043 WRITE(lunit,102)
'= 3: beta1 = 0. The exact solution is x = 0.'
4044 WRITE(lunit,102)
'= 4: A solution to (poss. singular) Ax = b found, given rtol.'
4045 WRITE(lunit,102)
'= 5: A solution to (poss. singular) Ax = b found, given eps.'
4046 WRITE(lunit,102)
'= 6: Pseudoinverse solution for singular LS problem, given rtol.'
4047 WRITE(lunit,102)
'= 7: Pseudoinverse solution for singular LS problem, given eps.'
4048 WRITE(lunit,102)
'= 8: The iteration limit was reached.'
4049 WRITE(lunit,102)
'= 9: The operator defined by Aprod appears to be unsymmetric.'
4050 WRITE(lunit,102)
'=10: The operator defined by Msolve appears to be unsymmetric.'
4051 WRITE(lunit,102)
'=11: The operator defined by Msolve appears to be indefinite.'
4052 WRITE(lunit,102)
'=12: xnorm has exceeded maxxnorm or will exceed it next iteration.'
4053 WRITE(lunit,102)
'=13: Acond has exceeded Acondlim or 0.1/eps.'
4054 WRITE(lunit,102)
'=14: Least-squares problem but no converged solution yet.'
4055 WRITE(lunit,102)
'=15: A null vector obtained, given rtol.'
4057 WRITE(lunit,101)
'ls',
'line search info'
4058 WRITE(lunit,102)
'< 0 recalculate function'
4059 WRITE(lunit,102)
'= 0: N or STP lt 0 or step not descending'
4060 WRITE(lunit,102)
'= 1: Linesearch convergence conditions reached'
4061 WRITE(lunit,102)
'= 2: interval of uncertainty at lower limit'
4062 WRITE(lunit,102)
'= 3: max nr of line search calls reached'
4063 WRITE(lunit,102)
'= 4: step at the lower bound'
4064 WRITE(lunit,102)
'= 5: step at the upper bound'
4065 WRITE(lunit,102)
'= 6: rounding error limitation'
4066 WRITE(lunit,101)
'step', &
4067 'the factor for the Newton step during the line search. Usually'
4069 'a value of 1 gives a sufficient reduction of the LF. Oherwise'
4070 WRITE(lunit,102)
'other step values are tried.'
4071 WRITE(lunit,101)
'cutf', &
4072 'cut factor. Local fits are rejected, if their chi^2 value'
4074 'is larger than the 3-sigma chi^2 value times the cut factor.'
4075 WRITE(lunit,102)
'A cut factor of 1 is used finally, but initially a larger'
4076 WRITE(lunit,102)
'factor may be used. A value of 0.0 means no cut.'
4077 WRITE(lunit,101)
'rejects',
'total number of rejected local fits.'
4078 WRITE(lunit,101)
'hmmsec',
'the time in hours (h), minutes (mm) and seconds.'
4079 WRITE(lunit,101)
'FMS',
'calculation of Function value, Matrix, Solution.'
4082101
FORMAT(a9,
' = ',a)
4099 INTEGER(mpi),
INTENT(IN) :: i
4100 INTEGER(mpi),
INTENT(IN) :: j
4101 REAL(mpd),
INTENT(IN) :: add
4103 INTEGER(mpl):: ijadd
4104 INTEGER(mpl):: ijcsr3
4109 IF(i <= 0.OR.j <= 0.OR. add == 0.0_mpd)
RETURN
4130 ELSE IF(
matsto == 2)
THEN
4158 CALL peend(23,
'Aborted, bad matrix index')
4159 stop
'mupdat: bad index'
4184 INTEGER(mpi),
INTENT(IN) :: i
4185 INTEGER(mpi),
INTENT(IN) :: j1
4186 INTEGER(mpi),
INTENT(IN) :: j2
4187 INTEGER(mpi),
INTENT(IN) :: il
4188 INTEGER(mpi),
INTENT(IN) :: jl
4189 INTEGER(mpi),
INTENT(IN) :: n
4190 REAL(mpd),
INTENT(IN) :: sub((n*n+n)/2)
4197 INTEGER(mpi):: iblast
4198 INTEGER(mpi):: iblock
4204 INTEGER(mpi):: jblast
4205 INTEGER(mpi):: jblock
4215 IF(i <= 0.OR.j1 <= 0.OR.j2 > i)
RETURN
4229 print *,
' MGUPDT: ij=0', i,j1,j2,il,jl,ij,lr,iprc,
matsto
4236 ijl=(k*k-k)/2+jl+ir-ia1
4253 ijl=(k*k-k)/2+jl+ir-ia1
4257 IF (jblock /= jblast.OR.iblock /= iblast)
THEN
4281 CALL ijpgrp(i,j1,ij,lr,iprc)
4283 print *,
' MGUPDT: ij=0', i,j1,j2,il,jl,ij,lr,iprc,
matsto
4347SUBROUTINE loopbf(nrej,numfil,naccf,chi2f,ndff)
4374 INTEGER(mpi) :: ibuf
4375 INTEGER(mpi) :: ichunk
4376 INTEGER(mpl) :: icmn
4377 INTEGER(mpl) :: icost
4379 INTEGER(mpi) :: idiag
4381 INTEGER(mpi) :: iext
4389 INTEGER(mpi) :: imeas
4392 INTEGER(mpi) :: ioffb
4393 INTEGER(mpi) :: ioffc
4394 INTEGER(mpi) :: ioffd
4395 INTEGER(mpi) :: ioffe
4396 INTEGER(mpi) :: ioffi
4397 INTEGER(mpi) :: ioffq
4398 INTEGER(mpi) :: iprc
4399 INTEGER(mpi) :: iprcnx
4400 INTEGER(mpi) :: iprdbg
4401 INTEGER(mpi) :: iproc
4402 INTEGER(mpi) :: irbin
4403 INTEGER(mpi) :: isize
4405 INTEGER(mpi) :: iter
4406 INTEGER(mpi) :: itgbi
4407 INTEGER(mpi) :: ivgbj
4408 INTEGER(mpi) :: ivgbk
4409 INTEGER(mpi) :: ivpgrp
4419 INTEGER(mpi) :: joffd
4420 INTEGER(mpi) :: joffi
4421 INTEGER(mpi) :: jproc
4425 INTEGER(mpi) :: kbdr
4426 INTEGER(mpi) :: kbdrx
4427 INTEGER(mpi) :: kbnd
4430 INTEGER(mpi) :: lvpgrp
4431 INTEGER(mpi) :: mbdr
4432 INTEGER(mpi) :: mbnd
4433 INTEGER(mpi) :: mside
4434 INTEGER(mpi) :: nalc
4435 INTEGER(mpi) :: nalg
4439 INTEGER(mpi) :: ndown
4441 INTEGER(mpi) :: nfred
4442 INTEGER(mpi) :: nfrei
4444 INTEGER(mpi) :: nprdbg
4445 INTEGER(mpi) :: nrank
4448 INTEGER(mpi) :: nter
4449 INTEGER(mpi) :: nweig
4450 INTEGER(mpi) :: ngrp
4451 INTEGER(mpi) :: npar
4453 INTEGER(mpl),
INTENT(IN OUT) :: nrej(6)
4454 INTEGER(mpi),
INTENT(IN) :: numfil
4455 INTEGER(mpi),
INTENT(IN OUT) :: naccf(numfil)
4456 REAL(mps),
INTENT(IN OUT) :: chi2f(numfil)
4457 INTEGER(mpi),
INTENT(IN OUT) :: ndff(numfil)
4467 INTEGER(mpi) :: ijprec
4474 CHARACTER (LEN=3):: chast
4475 DATA chuber/1.345_mpd/
4476 DATA cauchy/2.3849_mpd/
4524 IF(
nloopn == 1.AND.mod(nrc,100000_mpl) == 0)
THEN
4525 WRITE(*,*)
'Record',nrc,
' ... still reading'
4526 IF(
monpg1>0)
WRITE(
lunlog,*)
'Record',nrc,
' ... still reading'
4536 IF(nrc ==
nrecpr) lprnt=.true.
4537 IF(nrc ==
nrecp2) lprnt=.true.
4538 IF(nrc ==
nrecer) lprnt=.true.
4543 IF (nprdbg == 1) iprdbg=iproc
4544 IF (iproc /= iprdbg) lprnt=.false.
4549 WRITE(1,*)
'------------------ Loop',
nloopn, &
4550 ': Printout for record',nrc,iproc
4561 CALL isjajb(nst,ist,ja,jb,jsp)
4563 IF(imeas == 0)
WRITE(1,1121)
45681121
FORMAT(/
'Measured value and local derivatives'/ &
4569 ' i measured std_dev index...derivative ...')
45701122
FORMAT(i3,2g12.4,3(i3,g12.4)/(27x,3(i3,g12.4)))
4576 CALL isjajb(nst,ist,ja,jb,jsp)
4578 IF(imeas == 0)
WRITE(1,1123)
45901123
FORMAT(/
'Global derivatives'/ &
4591 ' i label gindex vindex derivative ...')
45921124
FORMAT(i3,2(i9,i7,i7,g12.4)/(3x,2(i9,i7,i7,g12.4)))
45931125
FORMAT(i3,2(i9,i7,i7,g12.4))
4603 WRITE(1,*)
'Data corrections using values of global parameters'
4604 WRITE(1,*)
'=================================================='
4614 CALL isjajb(nst,ist,ja,jb,jsp)
4646101
FORMAT(
' index measvalue corrvalue sigma')
4647102
FORMAT(i6,2x,2g12.4,
' +-',g12.4)
4649 IF(nalc <= 0)
GO TO 90
4651 ngg=(nalg*nalg+nalg)/2
4664 IF (ivpgrp /= lvpgrp)
THEN
4672 IF (npar /= nalg)
THEN
4673 print *,
' mismatch of number of global parameters ', nrc, nalg, npar, ngrp
4683 CALL peend(35,
'Aborted, mismatch of number of global parameters')
4684 stop
' mismatch of number of global parameters '
4711 DO WHILE(iter < nter)
4716 WRITE(1,*)
'Outlier-suppression iteration',iter,
' of',nter
4717 WRITE(1,*)
'=========================================='
4727 DO i=1,(nalc*nalc+nalc)/2
4738 wght =1.0_mpd/rerr**2
4743 IF(abs(resid) > chuber*rerr)
THEN
4744 wght=wght*chuber*rerr/abs(resid)
4748 wght=wght/(1.0+(resid/rerr/cauchy)**2)
4752 IF(lprnt.AND.iter /= 1.AND.nter /= 1)
THEN
4754 IF(abs(resid) > chuber*rerr) chast=
'* '
4755 IF(abs(resid) > 3.0*rerr) chast=
'** '
4756 IF(abs(resid) > 6.0*rerr) chast=
'***'
4757 IF(imeas == 0)
WRITE(1,*)
'Second loop: accumulate'
4758 IF(imeas == 0)
WRITE(1,103)
4763 WRITE(1,104) imeas,rmeas,resid,rerr,r1,chast,r2
4765103
FORMAT(
' index corrvalue residuum sigma', &
4767104
FORMAT(i6,2x,2g12.4,
' +-',g12.4,f7.2,1x,a3,f8.2)
4782 im=min(nalc+1-ij,nalc+1-ik)
4789 IF (iter == 1.AND.nalc > 5.AND.
lfitbb > 0)
THEN
4792 icmn=int(nalc,mpl)**3
4798 icost=6*int(nalc-kbdr,mpl)*int(kbnd+kbdr+1,mpl)**2+2*int(kbdr,mpl)**3
4799 IF (icost < icmn)
THEN
4811 kbdr=max(kbdr,
ibandh(k+nalc))
4812 icost=6*int(nalc-kbdr,mpl)*int(kbnd+kbdr+1,mpl)**2+2*int(kbdr,mpl)**3
4813 IF (icost < icmn)
THEN
4837 IF (
icalcm == 1.OR.lprnt) inv=2
4838 IF (mside == 1)
THEN
4846 IF (evdmin > 0.0_mpl) cndl10=log10(real(evdmax/evdmin,mps))
4856 IF ((.NOT.(
blvec(k) <= 0.0_mpd)).AND. (.NOT.(
blvec(k) > 0.0_mpd))) nan=nan+1
4861 WRITE(1,*)
'Parameter determination:',nalc,
' parameters,',
' rank=',nrank
4862 WRITE(1,*)
'-----------------------'
4863 IF(ndown /= 0)
WRITE(1,*)
' ',ndown,
' data down-weighted'
4879 wght =1.0_mpd/rerr**2
4903 IF (0.999999_mpd/wght > dvar)
THEN
4904 pull=rmeas/sqrt(1.0_mpd/wght-dvar)
4907 CALL hmpent(13,real(pull,mps))
4908 CALL gmpms(5,rec,real(pull,mps))
4910 CALL hmpent(14,real(pull,mps))
4929 IF(iter == 1.AND.jb < ist.AND.lhist) &
4930 CALL gmpms(4,rec,real(rmeas/rerr,mps))
4932 dchi2=wght*rmeas*rmeas
4937 IF(abs(resid) > chuber*rerr)
THEN
4938 wght=wght*chuber*rerr/abs(resid)
4939 dchi2=2.0*chuber*(abs(resid)/rerr-0.5*chuber)
4942 wght=wght/(1.0_mpd+(resid/rerr/cauchy)**2)
4943 dchi2=log(1.0_mpd+(resid/rerr/cauchy)**2)*cauchy**2
4953 IF(abs(resid) > chuber*rerr) chast=
'* '
4954 IF(abs(resid) > 3.0*rerr) chast=
'** '
4955 IF(abs(resid) > 6.0*rerr) chast=
'***'
4956 IF(imeas == 0)
WRITE(1,*)
'Third loop: single residuals'
4957 IF(imeas == 0)
WRITE(1,105)
4961 IF(resid < 0.0) r1=-r1
4962 IF(resid < 0.0) r2=-r2
4965105
FORMAT(
' index corrvalue residuum sigma', &
4967106
FORMAT(i6,2x,2g12.4,
' +-',g12.4,f7.2,1x,a3,f8.2)
4969 IF(iter == nter)
THEN
4971 resmax=max(resmax,abs(rmeas)/rerr)
4974 IF(iter == 1.AND.lhist)
THEN
4976 CALL hmpent( 3,real(rmeas/rerr,mps))
4978 CALL hmpent(12,real(rmeas/rerr,mps))
4985 resing=(real(nweig,mps)-real(suwt,mps))/real(nweig,mps)
4987 IF(iter == 1)
CALL hmpent( 5,real(ndf,mps))
4988 IF(iter == 1)
CALL hmpent(11,real(nalc,mps))
4993 WRITE(1,*)
'Chi^2=',summ,
' at',ndf,
' degrees of freedom: ', &
4994 '3-sigma limit is',chindl(3,ndf)*real(ndf,mps)
4995 WRITE(1,*) suwt,
' is sum of factors, compared to',nweig, &
4996 ' Downweight fraction:',resing
5002 WRITE(1,*)
' NaNs ', nalc, nrank, nan
5003 WRITE(1,*)
' ---> rejected!'
5005 IF (
mdebug < 0.AND.
nloopn == 1) print *,
' bad local fit-1 ', kfl,jrc,nrc,mside,neq,nalc,nrank,nan,cndl10
5008 IF(nrank /= nalc)
THEN
5012 WRITE(1,*)
' rank deficit', nalc, nrank
5013 WRITE(1,*)
' ---> rejected!'
5015 IF (
mdebug < 0.AND.
nloopn == 1) print *,
' bad local fit-2 ', kfl,jrc,nrc,mside,neq,nalc,nrank,nan,cndl10
5022 WRITE(1,*)
' too large condition(band part) ', nalc, nrank, cndl10
5023 WRITE(1,*)
' ---> rejected!'
5025 IF (
mdebug < 0.AND.
nloopn == 1) print *,
' bad local fit-3 ', kfl,jrc,nrc,mside,neq,nalc,nrank,nan,cndl10
5031 WRITE(1,*)
' Ndf<=0', nalc, nrank, ndf
5032 WRITE(1,*)
' ---> rejected!'
5034 IF (
mdebug < 0.AND.
nloopn == 1) print *,
' bad local fit-4 ', kfl,jrc,nrc,mside,neq,nalc,nrank,nan,cndl10
5038 chndf=real(summ/real(ndf,mpd),mps)
5040 IF(iter == 1.AND.lhist)
CALL hmpent(4,chndf)
5053 chichi=chindl(3,ndf)*real(ndf,mps)
5057 IF(summ >
chhuge*chichi)
THEN
5060 WRITE(1,*)
' ---> rejected!'
5068 IF(summ > chlimt)
THEN
5070 WRITE(1,*)
' ---> rejected!'
5074 CALL addsums(iproc+1, dchi2, ndf, dw1)
5084 CALL addsums(iproc+1, dchi2, ndf, dw1)
5088 WRITE(1,*)
' ---> rejected!'
5101 naccf(kfl)=naccf(kfl)+1
5102 ndff(kfl) =ndff(kfl) +ndf
5103 chi2f(kfl)=chi2f(kfl)+chndf
5116 wght =1.0_mpd/rerr**2
5117 dchi2=wght*rmeas*rmeas
5121 IF(resid > chuber*rerr)
THEN
5122 wght=wght*chuber*rerr/resid
5123 dchi2=2.0*chuber*(resid/rerr-0.5*chuber)
5173 CALL addsums(iproc+1, summ, ndf, dw1)
5229 IF (nfred < 0.OR.nfrei < 0)
THEN
5256 iprcnx=ijprec(i,jnx)
5258 IF (.NOT.((jnx == j+1).AND.(iprc == iprcnx)))
THEN
5272 joffd=joffd+(il*il-il)/2
5284 WRITE(1,*)
'------------------ End of printout for record',nrc
5341 iprcnx=ijprec(i,jnx)
5343 IF (.NOT.((jnx == j+1).AND.(iprc == iprcnx)))
THEN
5358 joffd=joffd+(il*il-il)/2
5394 INTEGER(mpi),
INTENT(IN) :: lun
5396 IF (
nrejec(1)>0)
WRITE(lun,*)
nrejec(1),
' (local solution contains NaNs)'
5397 IF (
nrejec(2)>0)
WRITE(lun,*)
nrejec(2),
' (local matrix with rank deficit)'
5398 IF (
nrejec(3)>0)
WRITE(lun,*)
nrejec(3),
' (local matrix with ill condition)'
5399 IF (
nrejec(4)>0)
WRITE(lun,*)
nrejec(4),
' (local fit with Ndf=0)'
5400 IF (
nrejec(5)>0)
WRITE(lun,*)
nrejec(5),
' (local fit with huge Chi2(Ndf))'
5401 IF (
nrejec(6)>0)
WRITE(lun,*)
nrejec(6),
' (local fit with large Chi2(Ndf))'
5427 INTEGER(mpi) :: icom
5428 INTEGER(mpl) :: icount
5432 INTEGER(mpi) :: imin
5433 INTEGER(mpi) :: iprlim
5434 INTEGER(mpi) :: isub
5435 INTEGER(mpi) :: itgbi
5436 INTEGER(mpi) :: itgbl
5437 INTEGER(mpi) :: ivgbi
5439 INTEGER(mpi) :: label
5447 INTEGER(mpi) :: labele(3)
5448 REAL(mps):: compnt(3)
5453 CALL mvopen(lup,
'millepede.res')
5456 WRITE(*,*)
' Result of fit for global parameters'
5457 WRITE(*,*)
' ==================================='
5462 WRITE(lup,*)
'Parameter ! first 3 elements per line are', &
5463 ' significant (if used as input)'
5481 err=sqrt(abs(real(gmati,mps)))
5482 IF(gmati < 0.0_mpd) err=-err
5485 IF(gmati*diag > 0.0_mpd)
THEN
5486 gcor2=1.0_mpd-1.0_mpd/(gmati*diag)
5487 IF(gcor2 >= 0.0_mpd.AND.gcor2 <= 1.0_mpd) gcor=real(sqrt(gcor2),mps)
5492 IF(lowstat) icount=-(icount+1)
5494 IF(itgbi <= iprlim)
THEN
5508 ELSE IF(itgbi == iprlim+1)
THEN
5509 WRITE(* ,*)
'... (further printout suppressed, but see log file)'
5523 ELSE IF (
igcorr /= 0)
THEN
5541 CALL mvopen(lup,
'millepede.eve')
5551 DO isub=0,min(15,imin-1)
5573 WRITE(lup,103) (labele(ie),compnt(ie),ie=1,iev)
5577 IF(iev /= 0)
WRITE(lup,103) (labele(ie),compnt(ie),ie=1,iev)
5585101
FORMAT(1x,
' label parameter presigma differ', &
5586 ' error'/ 1x,
'-----------',4x,4(
'-------------'))
5587102
FORMAT(i10,2x,4g14.5,f8.3)
5588103
FORMAT(3(i11,f11.7,2x))
5589110
FORMAT(i10,2x,2g14.5,28x,i12)
5590111
FORMAT(i10,2x,3g14.5,14x,i12)
5591112
FORMAT(i10,2x,4g14.5,i12)
5613 INTEGER(mpi) :: icom
5614 INTEGER(mpl) :: icount
5615 INTEGER(mpi) :: ifrst
5616 INTEGER(mpi) :: ilast
5617 INTEGER(mpi) :: inext
5618 INTEGER(mpi) :: itgbi
5619 INTEGER(mpi) :: itgbl
5620 INTEGER(mpi) :: itpgrp
5621 INTEGER(mpi) :: ivgbi
5623 INTEGER(mpi) :: icgrp
5624 INTEGER(mpi) :: ipgrp
5626 INTEGER(mpi) :: jpgrp
5628 INTEGER(mpi) :: label1
5629 INTEGER(mpi) :: label2
5630 INTEGER(mpi) :: ncon
5631 INTEGER(mpi) :: npair
5632 INTEGER(mpi) :: nstep
5635 INTEGER(mpl):: length
5637 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: vecPairedParGroups
5640 SUBROUTINE ggbmap(ipgrp,npair,npgrp)
5642 INTEGER(mpi),
INTENT(IN) :: ipgrp
5643 INTEGER(mpi),
INTENT(OUT) :: npair
5644 INTEGER(mpi),
DIMENSION(:),
INTENT(OUT) :: npgrp
5652 CALL mvopen(lup,
'millepede.res')
5653 WRITE(lup,*)
'*** Results of checking input only, no solution performed ***'
5654 WRITE(lup,*)
'! === global parameters ==='
5655 WRITE(lup,*)
'! fixed-1: by pre-sigma, -2: by entries cut, -3: by iterated entries cut'
5657 WRITE(lup,*)
'! Label Value Pre-sigma SkippedEntries Cons. group Status '
5659 WRITE(lup,*)
'! Label Value Pre-sigma Entries Cons. group Status '
5675 IF (ivgbi <= 0)
THEN
5677 IF (ivgbi == -4)
THEN
5678 WRITE(lup,116) c1,itgbl,par,presig,icount,icgrp
5680 WRITE(lup,110) c1,itgbl,par,presig,icount,icgrp,ivgbi
5684 WRITE(lup,111) c1,itgbl,par,presig,icount,icgrp
5690 WRITE(lup,*)
'!.Appearance statistics '
5691 WRITE(lup,*)
'!. Label First file and record Last file and record #files #paired-par'
5694 IF (itpgrp > 0)
THEN
5702 WRITE(lup,*)
'* === constraint groups ==='
5704 WRITE(lup,*)
'* Group #Cons. Entries First label Last label'
5706 WRITE(lup,*)
'* Group #Cons. Entries First label Last label Paired label range'
5708 CALL mpalloc(vecpairedpargroups,length,
'paired global parameter groups (I)')
5720 IF (
icheck > 1 .AND. label1 > 0)
THEN
5724 vecpairedpargroups(npair+1)=0
5728 jpgrp=vecpairedpargroups(j)
5732 IF (ifrst /= 0.AND.inext /= (ilast+nstep))
THEN
5735 WRITE(lup,114) label1, label2
5740 IF (ifrst == 0) ifrst=inext
5747 IF (jpgrp == vecpairedpargroups(j+1)-1)
THEN
5752 IF (ifrst /= 0)
THEN
5755 WRITE(lup,114) label1, label2
5761 WRITE(lup,*)
'*.Appearance statistics '
5762 WRITE(lup,*)
'*. Group First file and record Last file and record #files'
5772110
FORMAT(
' !',a1,i10,2x,2g14.5,2i12,
' fixed',i2)
5773111
FORMAT(
' !',a1,i10,2x,2g14.5,2i12,
' variable')
5774112
FORMAT(
' !.',i10,6i11)
5775113
FORMAT(
' * ',i6,i8,3i12)
5776114
FORMAT(
' *:',48x,i12,
' ..',i12)
5777115
FORMAT(
' *.',i10,5i11)
5778116
FORMAT(
' !',a1,i10,2x,2g14.5,2i12,
' redundant')
5808 INTEGER(mpi) :: iproc
5818 INTEGER(mpi),
INTENT(IN) :: n
5819 INTEGER(mpl),
INTENT(IN) :: l
5820 REAL(mpd),
INTENT(IN) :: x(n)
5821 INTEGER(mpi),
INTENT(IN) :: is
5822 INTEGER(mpi),
INTENT(IN) :: ie
5823 REAL(mpd),
INTENT(OUT) :: b(n)
5828 INTEGER(mpl) :: indij
5829 INTEGER(mpl) :: indid
5831 INTEGER(mpi) :: ichunk
5867 CALL peend(24,
'Aborted, vector/matrix size mismatch')
5868 stop
'AVPRDS: mismatched vector and matrix'
5888 IF (ia2 <= ib2) b(ia2:ib2)=b(ia2:ib2)+
globalmatd(ia2:ib2)*x(ia2:ib2)
5899 IF (i <= ie.AND.i >= is)
THEN
5907 IF (j <= ie.AND.j >= is)
THEN
5923 IF (ja2 <= jb2)
THEN
5926 b(i)=b(i)+dot_product(
globalmatd(indij+lj+ja2-ja:indij+lj+jb2-ja),x(ja2:jb2))
5930 IF (
mextnd == 0.AND.ia2 <= ib2)
THEN
5933 b(j)=b(j)+dot_product(
globalmatd(indij+lj+jn*(ia2-ia):indij+lj+jn*(ib2-ia):jn),x(ia2:ib2))
5952 IF (i <= ie.AND.i >= is)
THEN
5960 IF (j <= ie.AND.j >= is)
THEN
5976 IF (ja2 <= jb2)
THEN
5979 b(i)=b(i)+dot_product(real(
globalmatf(indij+lj+ja2-ja:indij+lj+jb2-ja),mpd),x(ja2:jb2))
5983 IF (
mextnd == 0.AND.ia2 <= ib2)
THEN
5986 b(j)=b(j)+dot_product(real(
globalmatf(indij+lj+jn*(ia2-ia):indij+lj+jn*(ib2-ia):jn),mpd),x(ia2:ib2))
6020 INTEGER(mpi) :: iproc
6028 INTEGER(mpi),
INTENT(IN) :: n
6029 INTEGER(mpl),
INTENT(IN) :: l
6030 REAL(mpd),
INTENT(IN) :: x(n)
6031 REAL(mpd),
INTENT(OUT) :: b(n)
6036 INTEGER(mpl) :: indij
6037 INTEGER(mpl) :: indid
6039 INTEGER(mpi) :: ichunk
6068 CALL peend(24,
'Aborted, vector/matrix size mismatch')
6069 stop
'AVPRD0: mismatched vector and matrix'
6115 b(i)=b(i)+dot_product(
globalmatd(indij+lj:indij+lj+jn-1),x(ja:jb))
6121 b(j)=b(j)+dot_product(
globalmatd(indij+lj:indij+jn*in:jn),x(ia:ib))
6157 b(i)=b(i)+dot_product(real(
globalmatf(indij+lj:indij+lj+jn-1),mpd),x(ja:jb))
6163 b(j)=b(j)+dot_product(real(
globalmatf(indij+lj:indij+jn*in:jn),mpd),x(ia:ib))
6187 INTEGER(mpi) :: ispc
6198 INTEGER(mpl) :: indid
6199 INTEGER(mpl),
DIMENSION(12) :: icount
6206 icount(4)=huge(icount(4))
6207 icount(7)=huge(icount(7))
6208 icount(10)=huge(icount(10))
6218 ir=ipg+(ispc-1)*(
napgrp+1)
6225 icount(1)=icount(1)+in
6226 icount(2)=icount(2)+in*ku
6232 icount(3)=icount(3)+1
6233 icount(4)=min(icount(4),jn)
6234 icount(5)=icount(5)+jn
6235 icount(6)=max(icount(6),jn)
6236 icount(7)=min(icount(7),in)
6237 icount(8)=icount(8)+in
6238 icount(9)=max(icount(9),in)
6239 icount(10)=min(icount(10),in*jn)
6240 icount(11)=icount(11)+in*jn
6241 icount(12)=max(icount(12),in*jn)
6247 WRITE(*,*)
"analysis of sparsity structure"
6248 IF (icount(1) > 0)
THEN
6249 WRITE(*,101)
"rows without compression/blocks ", icount(1)
6250 WRITE(*,101)
" contained elements ", icount(2)
6252 WRITE(*,101)
"number of block matrices ", icount(3)
6253 avg=real(icount(5),mps)/real(icount(3),mps)
6254 WRITE(*,101)
"number of columns (min,mean,max) ", icount(4), avg, icount(6)
6255 avg=real(icount(8),mps)/real(icount(3),mps)
6256 WRITE(*,101)
"number of rows (min,mean,max) ", icount(7), avg, icount(9)
6257 avg=real(icount(11),mps)/real(icount(3),mps)
6258 WRITE(*,101)
"number of elements (min,mean,max) ", icount(10), avg, icount(12)
6259101
FORMAT(2x,a34,i10,f10.3,i10)
6278 INTEGER(mpi),
INTENT(IN) :: n
6279 REAL(mpd),
INTENT(IN) :: x(n)
6280 REAL(mpd),
INTENT(OUT) :: b(n)
6285 CALL peend(24,
'Aborted, vector/matrix size mismatch')
6286 stop
'AVPROD: mismatched vector and matrix'
6317 INTEGER(mpi) :: ispc
6318 INTEGER(mpi) :: item1
6319 INTEGER(mpi) :: item2
6320 INTEGER(mpi) :: itemc
6321 INTEGER(mpi) :: jtem
6322 INTEGER(mpi) :: jtemn
6325 INTEGER(mpi),
INTENT(IN) :: itema
6326 INTEGER(mpi),
INTENT(IN) :: itemb
6327 INTEGER(mpl),
INTENT(OUT) :: ij
6328 INTEGER(mpi),
INTENT(OUT) :: lr
6329 INTEGER(mpi),
INTENT(OUT) :: iprc
6340 item1=max(itema,itemb)
6341 item2=min(itema,itemb)
6342 IF(item2 <= 0.OR.item1 >
napgrp)
RETURN
6345 outer:
DO ispc=1,
nspc
6356 IF(ku < kl) cycle outer
6361 IF(item2 >= jtem.AND.item2 < jtemn)
THEN
6367 IF(item2 < jtem)
THEN
6369 ELSE IF(item2 >= jtemn)
THEN
6384 IF(ku < kl) cycle outer
6388 IF(itemc == jtem)
EXIT
6389 IF(itemc < jtem)
THEN
6391 ELSE IF(itemc > jtem)
THEN
6419 INTEGER(mpi),
INTENT(IN) :: itema
6420 INTEGER(mpi),
INTENT(IN) :: itemb
6445 INTEGER(mpi) :: item1
6446 INTEGER(mpi) :: item2
6447 INTEGER(mpi) :: ipg1
6448 INTEGER(mpi) :: ipg2
6450 INTEGER(mpi) :: iprc
6452 INTEGER(mpi),
INTENT(IN) :: itema
6453 INTEGER(mpi),
INTENT(IN) :: itemb
6455 INTEGER(mpl) ::
ijadd
6459 item1=max(itema,itemb)
6460 item2=min(itema,itemb)
6462 IF(item2 <= 0.OR.item1 >
nagb)
RETURN
6463 IF(item1 == item2)
THEN
6472 CALL ijpgrp(ipg1,ipg2,ij,lr,iprc)
6494 INTEGER(mpi) :: item1
6495 INTEGER(mpi) :: item2
6496 INTEGER(mpi) :: jtem
6498 INTEGER(mpi),
INTENT(IN) :: itema
6499 INTEGER(mpi),
INTENT(IN) :: itemb
6508 item1=max(itema,itemb)
6509 item2=min(itema,itemb)
6511 IF(item2 <= 0.OR.item1 >
nagb)
RETURN
6519 print *,
' IJCSR3 empty list ', item1, item2, ks, ke
6520 CALL peend(23,
'Aborted, bad matrix index')
6521 stop
'ijcsr3: empty list'
6526 IF(item1 == jtem)
EXIT
6527 IF(item1 < jtem)
THEN
6534 print *,
' IJCSR3 not found ', item1, item2, ks, ke
6535 CALL peend(23,
'Aborted, bad matrix index')
6536 stop
'ijcsr3: not found'
6552 INTEGER(mpi) :: item1
6553 INTEGER(mpi) :: item2
6557 INTEGER(mpl) ::
ijadd
6560 INTEGER(mpi),
INTENT(IN) :: itema
6561 INTEGER(mpi),
INTENT(IN) :: itemb
6566 item1=max(itema,itemb)
6567 item2=min(itema,itemb)
6568 IF(item2 <= 0.OR.item1 >
nagb)
RETURN
6577 ij=
ijadd(item1,item2)
6580 ELSE IF (ij < 0)
THEN
6612 INTEGER(mpi) :: ichunk
6616 INTEGER(mpi) :: ispc
6624 INTEGER(mpl) :: ijadd
6646 ir=ipg+(ispc-1)*(
napgrp+1)
6696 REAL(mps),
INTENT(IN) :: deltat
6697 INTEGER(mpi),
INTENT(OUT) :: minut
6698 INTEGER(mpi),
INTENT(OUT):: nhour
6699 REAL(mps),
INTENT(OUT):: secnd
6700 INTEGER(mpi) :: nsecd
6703 nsecd=nint(deltat,
mpi)
6705 minut=nsecd/60-60*nhour
6706 secnd=deltat-60*(minut+60*nhour)
6742 INTEGER(mpi),
INTENT(IN) :: item
6746 INTEGER(mpl) :: length
6747 INTEGER(mpl),
PARAMETER :: four = 4
6751 IF(item <= 0)
RETURN
6772 IF(j == 0)
EXIT inner
6794 IF (
lvllog > 1)
WRITE(
lunlog,*)
'INONE: array increased to', &
6815 INTEGER(mpi) :: iprime
6816 INTEGER(mpi) :: nused
6817 LOGICAL :: finalUpdate
6818 INTEGER(mpl) :: oldLength
6819 INTEGER(mpl) :: newLength
6820 INTEGER(mpl),
PARAMETER :: four = 4
6821 INTEGER(mpi),
DIMENSION(:,:),
ALLOCATABLE :: tempArr
6822 INTEGER(mpl),
DIMENSION(:),
ALLOCATABLE :: tempVec
6826 IF(finalupdate)
THEN
6835 CALL mpalloc(temparr,four,oldlength,
'INONE: temp array')
6837 CALL mpalloc(tempvec,oldlength,
'INONE: temp vector')
6861 IF(j == 0)
EXIT inner
6862 IF(j == i) cycle outer
6866 IF(.NOT.finalupdate)
RETURN
6871 WRITE(
lunlog,*)
'INONE: array reduced to',newlength,
' words'
6895 IF(j == 0)
EXIT inner
6896 IF(j == i) cycle outer
6913 INTEGER(mpi),
INTENT(IN) :: n
6914 INTEGER(mpi) :: nprime
6915 INTEGER(mpi) :: nsqrt
6920 IF(mod(nprime,2) == 0) nprime=nprime+1
6923 nsqrt=int(sqrt(real(nprime,
mps)),
mpi)
6925 IF(i*(nprime/i) == nprime) cycle outer
6947 INTEGER(mpi) :: idum
6949 INTEGER(mpi) :: indab
6950 INTEGER(mpi) :: itgbi
6951 INTEGER(mpi) :: itgbl
6952 INTEGER(mpi) :: ivgbi
6954 INTEGER(mpi) :: jgrp
6955 INTEGER(mpi) :: lgrp
6957 INTEGER(mpi) :: nc31
6959 INTEGER(mpi) :: nwrd
6960 INTEGER(mpi) :: inone
6965 INTEGER(mpl) :: length
6966 INTEGER(mpl) :: rows
6970 WRITE(
lunlog,*)
'LOOP1: starting'
6993 WRITE(
lunlog,*)
'LOOP1: reading data files'
7004 CALL hmpldf(1,
'Number of words/record in binary file')
7005 CALL hmpdef(8,0.0,60.0,
'not_stored data per record')
7010 CALL peend(20,
'Aborted, bad binary records')
7011 stop
'LOOP1: length of binary record exceeds cache size, wrong file type?'
7025 WRITE(*,*)
'Read all binary data files:'
7044 WRITE(
lunlog,*)
'LOOP1: reading data files again'
7056 CALL peend(21,
'Aborted, no labels/parameters defined')
7057 stop
'LOOP1: no labels/parameters defined'
7062 ' is total number NTGB of labels/parameters'
7064 CALL hmpldf(2,
'Number of entries per label')
7108 WRITE(
lunlog,*)
'LOOP1:',
npresg,
' is number of pre-sigmas'
7109 WRITE(*,*)
'LOOP1:',
npresg,
' is number of pre-sigmas'
7110 IF(
npresg == 0)
WRITE(*,*)
'Warning: no pre-sigmas defined'
7132 WRITE(
lunlog,*)
'LOOP1:',
nvgb,
' is number NVGB of variable parameters'
7137 WRITE(
lunlog,*)
'LOOP1: counting records, NO iteration of entries cut !'
7143 CALL hmpdef(15,0.0,120.0,
'Number of parameters per group')
7178 ' is total number NTPGRP of label/parameter groups'
7194 WRITE(*,112)
' Default pre-sigma =',
regpre, &
7195 ' (if no individual pre-sigma defined)'
7196 WRITE(*,*)
'Pre-sigma factor is',
regula
7199 WRITE(*,*)
'No regularization will be done'
7201 WRITE(*,*)
'Regularization will be done, using factor',
regula
7205 CALL peend(22,
'Aborted, no variable global parameters')
7206 stop
'... no variable global parameters'
7213 IF(presg > 0.0)
THEN
7215 ELSE IF(presg == 0.0.AND.
regpre > 0.0)
THEN
7216 prewt=1.0/real(
regpre**2,mpd)
7232 WRITE(*,101)
'NTGB',
ntgb,
'total number of parameters'
7233 WRITE(*,101)
'NVGB',
nvgb,
'number of variable parameters'
7237 WRITE(*,101)
'Too many variable parameters for diagonalization, fallback is inversion'
7245 WRITE(*,101)
' NREC',
nrec,
'number of records'
7246 IF (
nrecd > 0)
WRITE(*,101)
' NRECD',
nrec,
'number of records containing doubles'
7247 WRITE(*,101)
' NEQN',
neqn,
'number of equations (measurements)'
7248 WRITE(*,101)
' NEGB',
negb,
'number of equations with global parameters'
7249 WRITE(*,101)
' NDGB',
ndgb,
'number of global derivatives'
7251 WRITE(*,101)
' NZGB',
nzgb,
'number of zero global der. (ignored in entry counts)'
7254 WRITE(*,101)
'MREQENF',
mreqenf,
'required number of entries (eqns in binary files)'
7256 WRITE(*,101)
'MREQENF',
mreqenf,
'required number of entries (recs in binary files)'
7259 WRITE(*,101)
'ITEREN',
iteren,
'iterate cut for parameters with less entries'
7260 WRITE(*,101)
'MREQENA',
mreqena,
'required number of entries (from accepted fits)'
7261 IF (
mreqpe > 1)
WRITE(*,101) &
7262 'MREQPE',
mreqpe,
'required number of pair entries'
7263 IF (
msngpe >= 1)
WRITE(*,101) &
7264 'MSNGPE',
msngpe,
'max pair entries single prec. storage'
7265 WRITE(*,101)
'NTGB',
ntgb,
'total number of parameters'
7266 WRITE(*,101)
'NVGB',
nvgb,
'number of variable parameters'
7269 WRITE(*,*)
'Global parameter labels:'
7276 mqi=((mqi-20)/20)*20+1
7284 WRITE(8,101)
' NREC',
nrec,
'number of records'
7285 IF (
nrecd > 0)
WRITE(8,101)
' NRECD',
nrec,
'number of records containing doubles'
7286 WRITE(8,101)
' NEQN',
neqn,
'number of equations (measurements)'
7287 WRITE(8,101)
' NEGB',
negb,
'number of equations with global parameters'
7288 WRITE(8,101)
' NDGB',
ndgb,
'number of global derivatives'
7290 WRITE(8,101)
'MREQENF',
mreqenf,
'required number of entries (eqns in binary files)'
7292 WRITE(8,101)
'MREQENF',
mreqenf,
'required number of entries (recs in binary files)'
7295 WRITE(8,101)
'ITEREN',
iteren,
'iterate cut for parameters with less entries'
7296 WRITE(8,101)
'MREQENA',
mreqena,
'required number of entries (from accepted fits)'
7298 WRITE(
lunlog,*)
'LOOP1: ending'
7302101
FORMAT(1x,a8,
' =',i14,
' = ',a)
7318 INTEGER(mpi) :: ibuf
7320 INTEGER(mpi) :: indab
7326 INTEGER(mpi) :: nc31
7328 INTEGER(mpi) :: nlow
7330 INTEGER(mpi) :: nwrd
7332 INTEGER(mpl) :: length
7333 INTEGER(mpl),
DIMENSION(:),
ALLOCATABLE :: newCounter
7338 WRITE(
lunlog,*)
'LOOP1: iterating'
7340 WRITE(*,*)
'LOOP1: iterating'
7343 CALL mpalloc(newcounter,length,
'new entries counter')
7367 CALL isjajb(nst,ist,ja,jb,jsp)
7368 IF(ja == 0.AND.jb == 0)
EXIT
7374 IF(ij == -2) nlow=nlow+1
7379 newcounter(ij)=newcounter(ij)+1
7408 WRITE(
lunlog,*)
'LOOP1:',
nvgb,
' is number NVGB of variable parameters'
7438 INTEGER(mpi) :: ibuf
7439 INTEGER(mpi) :: icblst
7440 INTEGER(mpi) :: icboff
7441 INTEGER(mpi) :: icgb
7442 INTEGER(mpi) :: icgrp
7443 INTEGER(mpi) :: icount
7444 INTEGER(mpi) :: iext
7445 INTEGER(mpi) :: ihis
7449 INTEGER(mpi) :: ioff
7450 INTEGER(mpi) :: ipoff
7451 INTEGER(mpi) :: iproc
7452 INTEGER(mpi) :: irecmm
7454 INTEGER(mpi) :: itgbi
7455 INTEGER(mpi) :: itgbij
7456 INTEGER(mpi) :: itgbik
7457 INTEGER(mpi) :: ivgbij
7458 INTEGER(mpi) :: ivgbik
7459 INTEGER(mpi) :: ivpgrp
7463 INTEGER(mpi) :: jcgrp
7464 INTEGER(mpi) :: jext
7465 INTEGER(mpi) :: jcgb
7466 INTEGER(mpi) :: jrec
7468 INTEGER(mpi) :: joff
7470 INTEGER(mpi) :: kcgrp
7471 INTEGER(mpi) :: kfile
7473 INTEGER(mpi) :: label
7474 INTEGER(mpi) :: labelf
7475 INTEGER(mpi) :: labell
7476 INTEGER(mpi) :: lvpgrp
7479 INTEGER(mpi) :: maeqnf
7480 INTEGER(mpi) :: nall
7481 INTEGER(mpi) :: naeqna
7482 INTEGER(mpi) :: naeqnf
7483 INTEGER(mpi) :: naeqng
7484 INTEGER(mpi) :: npdblk
7485 INTEGER(mpi) :: nc31
7486 INTEGER(mpi) :: ncachd
7487 INTEGER(mpi) :: ncachi
7488 INTEGER(mpi) :: ncachr
7489 INTEGER(mpi) :: ncon
7492 INTEGER(mpi) :: ndfmax
7493 INTEGER(mpi) :: nfixed
7494 INTEGER(mpi) :: nggd
7495 INTEGER(mpi) :: nggi
7496 INTEGER(mpi) :: nmatmo
7497 INTEGER(mpi) :: noff
7498 INTEGER(mpi) :: npair
7499 INTEGER(mpi) :: npar
7500 INTEGER(mpi) :: nparmx
7502 INTEGER(mpi) :: nrece
7503 INTEGER(mpi) :: nrecf
7504 INTEGER(mpi) :: nrecmm
7506 INTEGER(mpi) :: nwrd
7507 INTEGER(mpi) :: inone
7516 INTEGER(mpl):: nblock
7517 INTEGER(mpl):: nbwrds
7518 INTEGER(mpl):: noff8
7519 INTEGER(mpl):: ndimbi
7520 INTEGER(mpl):: ndimsa(4)
7522 INTEGER(mpl):: nnzero
7523 INTEGER(mpl):: matsiz(2)
7524 INTEGER(mpl):: matwords
7525 INTEGER(mpl):: mbwrds
7526 INTEGER(mpl):: length
7529 INTEGER(mpl),
PARAMETER :: two=2
7530 INTEGER(mpi) :: maxGlobalPar = 0
7531 INTEGER(mpi) :: maxLocalPar = 0
7532 INTEGER(mpi) :: maxEquations = 0
7534 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: vecConsGroupList
7535 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: vecConsGroupIndex
7536 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: vecPairedParGroups
7537 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: vecBlockCounts
7540 SUBROUTINE ndbits(npgrp,ndims,nsparr,ihst)
7542 INTEGER(mpi),
DIMENSION(:),
INTENT(IN) :: npgrp
7543 INTEGER(mpl),
DIMENSION(4),
INTENT(OUT) :: ndims
7544 INTEGER(mpl),
DIMENSION(:,:),
INTENT(OUT) :: nsparr
7545 INTEGER(mpi),
INTENT(IN) :: ihst
7547 SUBROUTINE ckbits(npgrp,ndims)
7549 INTEGER(mpi),
DIMENSION(:),
INTENT(IN) :: npgrp
7550 INTEGER(mpl),
DIMENSION(4),
INTENT(OUT) :: ndims
7552 SUBROUTINE spbits(npgrp,nsparr,nsparc)
7554 INTEGER(mpi),
DIMENSION(:),
INTENT(IN) :: npgrp
7555 INTEGER(mpl),
DIMENSION(:,:),
INTENT(IN) :: nsparr
7556 INTEGER(mpi),
DIMENSION(:),
INTENT(OUT) :: nsparc
7558 SUBROUTINE gpbmap(ngroup,npgrp,npair)
7560 INTEGER(mpi),
INTENT(IN) :: ngroup
7561 INTEGER(mpi),
DIMENSION(:,:),
INTENT(IN) :: npgrp
7562 INTEGER(mpi),
DIMENSION(:),
INTENT(OUT) :: npair
7564 SUBROUTINE ggbmap(ipgrp,npair,npgrp)
7566 INTEGER(mpi),
INTENT(IN) :: ipgrp
7567 INTEGER(mpi),
INTENT(OUT) :: npair
7568 INTEGER(mpi),
DIMENSION(:),
INTENT(OUT) :: npgrp
7570 SUBROUTINE pbsbits(npgrp,ibsize,nnzero,nblock,nbkrow)
7572 INTEGER(mpi),
DIMENSION(:),
INTENT(IN) :: npgrp
7573 INTEGER(mpi),
INTENT(IN) :: ibsize
7574 INTEGER(mpl),
INTENT(OUT) :: nnzero
7575 INTEGER(mpl),
INTENT(OUT) :: nblock
7576 INTEGER(mpi),
DIMENSION(:),
INTENT(OUT) :: nbkrow
7578 SUBROUTINE pblbits(npgrp,ibsize,nsparr,nsparc)
7580 INTEGER(mpi),
DIMENSION(:),
INTENT(IN) :: npgrp
7581 INTEGER(mpi),
INTENT(IN) :: ibsize
7582 INTEGER(mpl),
DIMENSION(:),
INTENT(IN) :: nsparr
7583 INTEGER(mpl),
DIMENSION(:),
INTENT(OUT) :: nsparc
7585 SUBROUTINE prbits(npgrp,nsparr)
7587 INTEGER(mpi),
DIMENSION(:),
INTENT(IN) :: npgrp
7588 INTEGER(mpl),
DIMENSION(:),
INTENT(OUT) :: nsparr
7590 SUBROUTINE pcbits(npgrp,nsparr,nsparc)
7592 INTEGER(mpi),
DIMENSION(:),
INTENT(IN) :: npgrp
7593 INTEGER(mpl),
DIMENSION(:),
INTENT(IN) :: nsparr
7594 INTEGER(mpl),
DIMENSION(:),
INTENT(OUT) :: nsparc
7604 WRITE(
lunlog,*)
'LOOP2: starting'
7627 WRITE(
lunlog,*)
' Elimination for constraints enforced for solution by decomposition!'
7632 WRITE(
lunlog,*)
' Lagrange multipliers enforced for solution by sparsePARDISO!'
7637 WRITE(
lunlog,*)
' Elimination for constraints with mpqldec enforced (LAPACK only for unpacked storage)!'
7654 noff8=int(
nagb,mpl)*int(
nagb-1,mpl)/2
7690 print *,
' Variable parameter groups ',
nvpgrp
7738 CALL mpalloc(vecconsgrouplist,length,
'constraint group list')
7739 CALL mpalloc(vecconsgroupindex,length,
'constraint group index')
7749 WRITE(
lunlog,*)
'LOOP2: start event reading'
7789 WRITE(*,*)
'Record number ',
nrec,
' from file ',kfile
7790 IF (wgh /= 1.0)
WRITE(*,*)
' weight ',wrec
7794 CALL isjajb(nst,ist,ja,jb,jsp)
7798 IF(nda ==
mdebg2+1)
WRITE(*,*)
'... and more data'
7803 WRITE(*,*)
'Local derivatives:'
7805107
FORMAT(6(i3,g12.4))
7807 WRITE(*,*)
'Global derivatives:'
7810108
FORMAT(3i11,g12.4)
7813 WRITE(*,*)
'total_par_label __label__ var_par_index derivative'
7828 CALL isjajb(nst,ist,ja,jb,jsp)
7829 IF(ja == 0.AND.jb == 0)
EXIT
7869 IF (kcgrp == jcgrp) cycle
7880 IF (vecconsgroupindex(k) == 0)
THEN
7883 vecconsgrouplist(icgrp)=k
7898 IF (vecconsgroupindex(k) < icount)
THEN
7900 vecconsgroupindex(k)=icount
7918 IF (nfixed > 0) naeqnf=naeqnf+1
7921 IF(ja /= 0.AND.jb /= 0)
THEN
7930 IF (naeqnf > maeqnf) nrecf=nrecf+1
7934 maxglobalpar=max(
nagbn,maxglobalpar)
7935 maxlocalpar=max(
nalcn,maxlocalpar)
7936 maxequations=max(
naeqn,maxequations)
7939 dstat(1)=dstat(1)+real((nwrd+2)*2,mpd)
7940 dstat(2)=dstat(2)+real(
nagbn+2,mpd)
7945 vecconsgroupindex(vecconsgrouplist(k))=0
7950 IF (
nagbn == 0)
THEN
7969 IF (ivpgrp /= lvpgrp)
THEN
8024 (irecmm >= nrecmm.OR.irecmm ==
mxrec))
THEN
8025 IF (nmatmo == 0)
THEN
8027 WRITE(*,*)
'Monitoring of sparse matrix construction'
8028 WRITE(*,*)
' records ........ off-diagonal elements ', &
8029 '....... compression memory'
8030 WRITE(*,*)
' non-zero used(double) used', &
8035 gbc=1.0e-9*real((mpi*ndimsa(2)+mpd*ndimsa(3)+mps*ndimsa(4))/mpi*(bit_size(1_mpi)/8),mps)
8036 gbu=1.0e-9*real(((mpi+mpd)*(ndimsa(3)+ndimsa(4)))/mpi*(bit_size(1_mpi)/8),mps)
8038 WRITE(*,1177) irecmm,ndimsa(1),ndimsa(3),ndimsa(4),cpr,gbc
80391177
FORMAT(i9,3i13,f10.2,f11.6)
8040 DO WHILE(irecmm >= nrecmm)
8059 WRITE(
lunlog,*)
'LOOP2: event reading ended - end of data'
8061 dstat(k)=dstat(k)/real(
nrec,mpd)
8075 CALL mpalloc(vecpairedpargroups,length,
'paired global parameter groups (I)')
8077 print *,
' Total parameter groups pairs',
ntpgrp
8080 CALL ggbmap(i,npair,vecpairedpargroups)
8148 IF(ivgbij > 0.AND.ivgbik > 0)
THEN
8150 IF (
mprint > 1)
WRITE(*,*)
'add index pair ',ivgbij,ivgbik
8227 nparmx=max(nparmx,int(rows,mpi))
8249 WRITE(*,*)
' Parameter block', i, ib-ia, jb-ja, labelf, labell
8253 WRITE(
lunlog,*)
'Detected',
npblck,
'(disjoint) parameter blocks, max size ', nparmx
8255 WRITE(*,*)
'Detected',
npblck,
'(disjoint) parameter blocks, max size ', nparmx
8257 WRITE(*,*)
'Using block diagonal storage mode'
8272 IF (
nagb >= 65536)
THEN
8273 noff=int(noff8/1000,mpi)
8283 CALL hmpdef(ihis,0.0,real(
mhispe,mps),
'NDBITS: #off-diagonal elements')
8288 ndgn=ndimsa(3)+ndimsa(4)
8289 matwords=ndimsa(2)+length*4
8304 length=int(npdblk,mpl)
8305 CALL mpalloc(vecblockcounts,length,
'sparse matrix row offsets (CSR3)')
8307 nbwrds=2*int(nblock,mpl)*int(
ipdbsz(i)*
ipdbsz(i)+1,mpl)
8308 IF ((i == 1).OR.(nbwrds < mbwrds))
THEN
8332 IF (
fcache(2) == 0.0)
THEN
8334 fcache(2)=real(dstat(2),mps)
8335 fcache(3)=real(dstat(3),mps)
8356 ncachd=
ncache-ncachr-ncachi
8387 WRITE(lu,101)
'NTGB',
ntgb,
'total number of parameters'
8388 WRITE(lu,102)
'(all parameters, appearing in binary files)'
8389 WRITE(lu,101)
'NVGB',
nvgb,
'number of variable parameters'
8390 WRITE(lu,102)
'(appearing in fit matrix/vectors)'
8391 WRITE(lu,101)
'NAGB',
nagb,
'number of all parameters'
8392 WRITE(lu,102)
'(including Lagrange multiplier or reduced)'
8393 WRITE(lu,101)
'NTPGRP',
ntpgrp,
'total number of parameter groups'
8394 WRITE(lu,101)
'NVPGRP',
nvpgrp,
'number of variable parameter groups'
8395 WRITE(lu,101)
'NFGB',
nfgb,
'number of fit parameters'
8397 WRITE(lu,101)
'MBANDW',
mbandw,
'band width of preconditioner matrix'
8398 WRITE(lu,102)
'(if <0, no preconditioner matrix)'
8400 IF (
nagb >= 65536)
THEN
8401 WRITE(lu,101)
'NOFF/K',noff,
'max number of off-diagonal elements'
8403 WRITE(lu,101)
'NOFF',noff,
'max number of off-diagonal elements'
8406 IF (
nagb >= 65536)
THEN
8407 WRITE(lu,101)
'NDGN/K',ndgn/1000,
'actual number of off-diagonal elements'
8409 WRITE(lu,101)
'NDGN',ndgn,
'actual number of off-diagonal elements'
8412 WRITE(lu,101)
'NCGB',
ncgb,
'number of constraints'
8413 WRITE(lu,101)
'NAGBN',
nagbn,
'max number of global parameters in an event'
8414 WRITE(lu,101)
'NALCN',
nalcn,
'max number of local parameters in an event'
8415 WRITE(lu,101)
'NAEQN',
naeqn,
'max number of equations in an event'
8417 WRITE(lu,101)
'NAEQNA',naeqna,
'number of equations'
8418 WRITE(lu,101)
'NAEQNG',naeqng, &
8419 'number of equations with global parameters'
8420 WRITE(lu,101)
'NAEQNF',naeqnf, &
8421 'number of equations with fixed global parameters'
8422 WRITE(lu,101)
'NRECF',nrecf, &
8423 'number of records with fixed global parameters'
8426 WRITE(lu,101)
'NRECE',nrece, &
8427 'number of records without variable parameters'
8430 WRITE(lu,101)
'NCACHE',
ncache,
'number of words for caching'
8431 WRITE(lu,111) (
fcache(k)*100.0,k=1,3)
8432111
FORMAT(22x,
'cache splitting ',3(f6.1,
' %'))
8437 WRITE(lu,*)
'Solution method and matrix-storage mode:'
8439 WRITE(lu,*)
' METSOL = 1: matrix inversion'
8440 ELSE IF(
metsol == 2)
THEN
8441 WRITE(lu,*)
' METSOL = 2: diagonalization'
8442 ELSE IF(
metsol == 3)
THEN
8443 WRITE(lu,*)
' METSOL = 3: decomposition'
8444 ELSE IF(
metsol == 4)
THEN
8445 WRITE(lu,*)
' METSOL = 4: MINRES (rtol',
mrestl,
')'
8446 ELSE IF(
metsol == 5)
THEN
8447 WRITE(lu,*)
' METSOL = 5: MINRES-QLP (rtol',
mrestl,
')'
8448 ELSE IF(
metsol == 6)
THEN
8449 WRITE(lu,*)
' METSOL = 6: GMRES'
8451 ELSE IF(
metsol == 7)
THEN
8452 WRITE(lu,*)
' METSOL = 7: LAPACK factorization'
8453 ELSE IF(
metsol == 8)
THEN
8454 WRITE(lu,*)
' METSOL = 8: LAPACK factorization'
8456 ELSE IF(
metsol == 9)
THEN
8457 WRITE(lu,*)
' METSOL = 9: Intel oneMKL PARDISO'
8461 WRITE(lu,*)
' with',
mitera,
' iterations'
8463 WRITE(lu,*)
' MATSTO = 0: unpacked symmetric matrix, ',
'n*n elements'
8464 ELSE IF(
matsto == 1)
THEN
8465 WRITE(lu,*)
' MATSTO = 1: full symmetric matrix, ',
'(n*n+n)/2 elements'
8466 ELSE IF(
matsto == 2)
THEN
8467 WRITE(lu,*)
' MATSTO = 2: sparse matrix (custom)'
8468 ELSE IF(
matsto == 3)
THEN
8470 WRITE(lu,*)
' MATSTO = 3: sparse matrix (upper triangle, CSR3)'
8472 WRITE(lu,*)
' MATSTO = 3: sparse matrix (upper triangle, BSR3)'
8473 WRITE(lu,*)
' block size',
matbsz
8477 WRITE(lu,*)
' block diagonal with',
npblck,
' blocks'
8479 IF(
mextnd>0)
WRITE(lu,*)
' with extended storage'
8480 IF(
dflim /= 0.0)
THEN
8481 WRITE(lu,103)
'Convergence assumed, if expected dF <',
dflim
8486 WRITE(lu,*)
'Constraints handled by elimination with LAPACK'
8488 WRITE(lu,*)
'Constraints handled by elimination'
8491 WRITE(lu,*)
'Constraints handled by Lagrange multipliers'
8498 CALL peend(28,
'Aborted, no local parameters')
8499 stop
'LOOP2: stopping due to missing local parameters'
8520105
FORMAT(
' Constants C1, C2 for Wolfe conditions:',g12.4,
', ',g12.4)
8528 matwords=(length+
nagb+1)*2
8535 ELSE IF (
matsto == 2)
THEN
8536 matsiz(1)=ndimsa(3)+
nagb
8551 IF (icblst > icboff)
THEN
8557 nall = npar;
IF (
icelim <= 0) nall=npar+ncon
8561 matsiz(1)=matsiz(1)+k
8563 matsiz(1)=matsiz(1)+nall
8568 matwords=matwords+matsiz(1)*2+matsiz(2)
8574 IF (matwords < 250000)
THEN
8575 WRITE(*,*)
'Size of global matrix: < 1 MB'
8577 WRITE(*,*)
'Size of global matrix:',int(real(matwords,mps)*4.0e-6,mpi),
' MB'
8583 WRITE(
lunlog,*)
' Cut values of Chi^2/Ndf and Chi2,'
8584 WRITE(
lunlog,*)
' corresponding to 2 and 3 standard deviations'
8585 WRITE(
lunlog,*)
' Ndf Chi^2/Ndf(2) Chi^2(2) ', &
8586 ' Chi^2/Ndf(3) Chi^2(3)'
8589 IF(ndf >
naeqn)
EXIT
8592 ELSE IF(ndf < 20)
THEN
8594 ELSE IF(ndf < 100)
THEN
8596 ELSE IF(ndf < 200)
THEN
8603 WRITE(
lunlog,106) ndf,chin2,chin2*real(ndf,mps),chin3, chin3*real(ndf,mps)
8606 WRITE(
lunlog,*)
'LOOP2: ending'
8610 IF (
ncgbe /= 0)
THEN
8613 WRITE(*,199)
'WarningWarningWarningWarningWarningWarningWarningWarningWar'
8614 WRITE(*,199)
'arningWarningWarningWarningWarningWarningWarningWarningWarn'
8615 WRITE(*,199)
'rningWarningWarningWarningWarningWarningWarningWarningWarni'
8616 WRITE(*,199)
'ningWarningWarningWarningWarningWarningWarningWarningWarnin'
8617 WRITE(*,199)
'ingWarningWarningWarningWarningWarningWarningWarningWarning'
8618 WRITE(*,199)
'ngWarningWarningWarningWarningWarningWarningWarningWarningW'
8619 WRITE(*,199)
'gWarningWarningWarningWarningWarningWarningWarningWarningWa'
8621 WRITE(*,*)
' Number of empty constraints =',abs(
ncgbe),
', should be 0'
8622 WRITE(*,*)
' => please check constraint definition, mille data'
8624 WRITE(*,199)
'WarningWarningWarningWarningWarningWarningWarningWarningWar'
8625 WRITE(*,199)
'arningWarningWarningWarningWarningWarningWarningWarningWarn'
8626 WRITE(*,199)
'rningWarningWarningWarningWarningWarningWarningWarningWarni'
8627 WRITE(*,199)
'ningWarningWarningWarningWarningWarningWarningWarningWarnin'
8628 WRITE(*,199)
'ingWarningWarningWarningWarningWarningWarningWarningWarning'
8629 WRITE(*,199)
'ngWarningWarningWarningWarningWarningWarningWarningWarningW'
8630 WRITE(*,199)
'gWarningWarningWarningWarningWarningWarningWarningWarningWa'
8635101
FORMAT(1x,a8,
' =',i14,
' = ',a)
8637103
FORMAT(1x,a,g12.4)
8638106
FORMAT(i6,2(3x,f9.3,f12.1,3x))
8653 INTEGER(mpi) :: imed
8656 INTEGER(mpi) :: nent
8657 INTEGER(mpi),
DIMENSION(measBins) :: isuml
8658 INTEGER(mpi),
DIMENSION(measBins) :: isums
8662 INTEGER(mpl) :: ioff
8665 DATA lfirst /.true./
8678 WRITE(
lunmon,
'(A)')
'*** Normalized residuals grouped by first global label (per local fit cycle) ***'
8680 WRITE(
lunmon,
'(A)')
'*** Pulls grouped by first global label (per local fit cycle) ***'
8682 WRITE(
lunmon,
'(A)')
'! LFC Label Entries Median RMS(MAD) <error>'
8687#ifdef SCOREP_USER_ENABLE
8688 scorep_user_region_by_name_begin(
"UR_monres", scorep_user_region_type_common)
8704 IF (2*isuml(j) > nent)
EXIT
8708 IF (isuml(j) > isuml(j-1)) amed=amed+real(nent-2*isuml(j-1),mps)/real(2*isuml(j)-2*isuml(j-1),mps)
8721 isums(j)=isums(j)+isums(j-1)
8725 IF (2*isums(j) > nent)
EXIT
8728 IF (isums(j) > isums(j-1)) amad=amad+real(nent-2*isums(j-1),mps)/real(2*isums(j)-2*isums(j-1),mps)
8740#ifdef SCOREP_USER_ENABLE
8741 scorep_user_region_by_name_end(
"UR_monres")
8745110
FORMAT(i5,2i10,3g14.5)
8760 INTEGER(mpi) :: ioff
8761 INTEGER(mpi) :: ipar0
8762 INTEGER(mpi) :: ncon
8763 INTEGER(mpi) :: npar
8764 INTEGER(mpi) :: nextra
8766 INTEGER :: nbopt, nboptx, ILAENV
8769 INTEGER(mpl),
INTENT(IN) :: msize(2)
8771 INTEGER(mpl) :: length
8772 INTEGER(mpl) :: nwrdpc
8773 INTEGER(mpl),
PARAMETER :: three = 3
8786 CALL mpalloc(
aux,length,
' local fit scratch array: aux')
8787 CALL mpalloc(
vbnd,length,
' local fit scratch array: vbnd')
8788 CALL mpalloc(
vbdr,length,
' local fit scratch array: vbdr')
8791 CALL mpalloc(
vbk,length,
' local fit scratch array: vbk')
8794 CALL mpalloc(
vzru,length,
' local fit scratch array: vzru')
8795 CALL mpalloc(
scdiag,length,
' local fit scratch array: scdiag')
8796 CALL mpalloc(
scflag,length,
' local fit scratch array: scflag')
8797 CALL mpalloc(
ibandh,2*length,
' local fit band width hist.: ibandh')
8811 length=4*
ncblck;
IF(ncon == 0) length=0
8820 nwrdpc=nwrdpc+length
8832 length=(ncon*ncon+ncon)/2
8857 length=length+npar+nextra
8863 IF(
mbandw == 0) length=length+1
8871 nwrdpc=nwrdpc+2*length
8872 IF (nwrdpc > 250000)
THEN
8874 WRITE(*,*)
'Size of preconditioner matrix:',int(real(nwrdpc,mps)*4.0e-6,mpi),
' MB'
8896 CALL peend(23,
'Aborted, bad matrix index (will exceed 32bit)')
8897 stop
'vmprep: bad index (matrix to large for diagonalization)'
8919 nbopt = ilaenv( 1_mpl,
'DSYTRF',
'U', int(
nagb,mpl), int(
nagb,mpl), -1_mpl, -1_mpl )
8921 print *,
'LAPACK optimal block size for DSYTRF:', nbopt
8922 lplwrk=length*int(nbopt,mpl)
8930 nbopt = ilaenv( 1_mpl,
'DORMQL',
'RN', int(npar,mpl), int(npar,mpl), int(ncon,mpl), int(npar,mpl) )
8931 IF (int(npar,mpl)*int(nbopt,mpl) >
lplwrk)
THEN
8932 lplwrk=int(npar,mpl)*int(nbopt,mpl)
8937 print *,
'LAPACK optimal block size for DORMQL:', nboptx
8956 INTEGER(mpi) :: icoff
8957 INTEGER(mpi) :: ipoff
8960 INTEGER(mpi) :: ncon
8961 INTEGER(mpi) :: nfit
8962 INTEGER(mpi) :: npar
8963 INTEGER(mpi) :: nrank
8964 INTEGER(mpl) :: imoff
8965 INTEGER(mpl) :: ioff1
8983 WRITE(
lunlog,*)
'Shrinkage of global matrix (A->Q^t*A*Q)'
8998 nfit=npar+ncon;
IF (
icelim > 0) nfit=npar-ncon
9000 IF(nfit < npar)
THEN
9018 WRITE(
lunlog,*)
'Inversion of global matrix (A->A^-1)'
9025 IF(nfit /= nrank)
THEN
9026 WRITE(*,*)
'Warning: the rank defect of the symmetric',nfit, &
9027 '-by-',nfit,
' matrix is ',nfit-nrank,
' (should be zero).'
9028 WRITE(lun,*)
'Warning: the rank defect of the symmetric',nfit, &
9029 '-by-',nfit,
' matrix is ',nfit-nrank,
' (should be zero).'
9032 WRITE(*,*)
' --> enforcing SUBITO mode'
9033 WRITE(lun,*)
' --> enforcing SUBITO mode'
9035 ELSE IF(
ndefec == 0)
THEN
9037 WRITE(lun,*)
'No rank defect of the symmetric matrix'
9039 WRITE(lun,*)
'No rank defect of the symmetric block', ib,
' of size', npar
9050 IF(nfit < npar)
THEN
9069 INTEGER(mpi) :: icoff
9070 INTEGER(mpi) :: ipoff
9073 INTEGER(mpi) :: ncon
9074 INTEGER(mpi) :: nfit
9075 INTEGER(mpi) :: npar
9076 INTEGER(mpi) :: nrank
9077 INTEGER(mpl) :: imoff
9078 INTEGER(mpl) :: ioff1
9093 WRITE(
lunlog,*)
'Shrinkage of global matrix (A->Q^t*A*Q)'
9107 nfit=npar+ncon;
IF (
icelim > 0) nfit=npar-ncon
9109 IF(nfit < npar)
THEN
9127 WRITE(
lunlog,*)
'Decomposition of global matrix (A->L*D*L^t)'
9133 IF(nfit /= nrank)
THEN
9134 WRITE(*,*)
'Warning: the rank defect of the symmetric',nfit, &
9135 '-by-',nfit,
' matrix is ',nfit-nrank,
' (should be zero).'
9136 WRITE(lun,*)
'Warning: the rank defect of the symmetric',nfit, &
9137 '-by-',nfit,
' matrix is ',nfit-nrank,
' (should be zero).'
9140 WRITE(*,*)
' --> enforcing SUBITO mode'
9141 WRITE(lun,*)
' --> enforcing SUBITO mode'
9143 ELSE IF(
ndefec == 0)
THEN
9145 WRITE(lun,*)
'No rank defect of the symmetric matrix'
9147 WRITE(lun,*)
'No rank defect of the symmetric block', ib,
' of size', npar
9149 WRITE(lun,*)
' largest diagonal element (LDLt)', evmax
9150 WRITE(lun,*)
' smallest diagonal element (LDLt)', evmin
9159 IF(nfit < npar)
THEN
9181 INTEGER(mpi) :: icoff
9182 INTEGER(mpi) :: ipoff
9185 INTEGER(mpi) :: ncon
9186 INTEGER(mpi) :: nfit
9187 INTEGER(mpi) :: npar
9188 INTEGER(mpl) :: imoff
9189 INTEGER(mpl) :: ioff1
9190 INTEGER(mpi) :: infolp
9210 WRITE(
lunlog,*)
'Shrinkage of global matrix (A->Q^t*A*Q)'
9225 nfit=npar+ncon;
IF (
icelim > 0) nfit=npar-ncon
9227 IF(nfit < npar)
THEN
9244 IF (nfit > npar)
THEN
9247 WRITE(
lunlog,*)
'Factorization of global matrix (A->L*D*L^t)'
9251#ifdef SCOREP_USER_ENABLE
9252 scorep_user_region_by_name_begin(
"UR_dsptrf", scorep_user_region_type_common)
9255#ifdef SCOREP_USER_ENABLE
9256 scorep_user_region_by_name_end(
"UR_dsptrf")
9263 WRITE(
lunlog,*)
'Factorization of global matrix (A->L*L^t)'
9267#ifdef SCOREP_USER_ENABLE
9268 scorep_user_region_by_name_begin(
"UR_dpptrf", scorep_user_region_type_common)
9270 CALL dpptrf(
'U',int(nfit,mpl),
globalmatd(imoff+1:),infolp)
9271#ifdef SCOREP_USER_ENABLE
9272 scorep_user_region_by_name_end(
"UR_dpptrf")
9280 WRITE(lun,*)
'No rank defect of the symmetric matrix'
9282 WRITE(lun,*)
'No rank defect of the symmetric block', ib,
' of size', npar
9286 WRITE(*,*)
'Warning: factorization of the symmetric',nfit, &
9287 '-by-',nfit,
' failed at index ', infolp
9288 WRITE(lun,*)
'Warning: factorization of the symmetric',nfit, &
9289 '-by-',nfit,
' failed at index ', infolp
9290 CALL peend(29,
'Aborted, factorization of global matrix failed')
9291 stop
'mdptrf: bad matrix'
9296 IF (nfit > npar)
THEN
9299 IF(infolp /= 0) print *,
' DSPTRS failed: ', infolp
9301 CALL dpptrs(
'U',int(nfit,mpl),1_mpl,
globalmatd(imoff+1:),&
9303 IF(infolp /= 0) print *,
' DPPTRS failed: ', infolp
9307 IF(nfit < npar)
THEN
9328 INTEGER(mpi) :: icoff
9329 INTEGER(mpi) :: ipoff
9332 INTEGER(mpi) :: ncon
9333 INTEGER(mpi) :: nfit
9334 INTEGER(mpi) :: npar
9335 INTEGER(mpl) :: imoff
9336 INTEGER(mpl) :: ioff1
9337 INTEGER(mpl) :: iloff
9338 INTEGER(mpi) :: infolp
9359 WRITE(
lunlog,*)
'Shrinkage of global matrix (A->Q^t*A*Q)'
9379 nfit=npar+ncon;
IF (
icelim > 0) nfit=npar-ncon
9381 IF(nfit < npar)
THEN
9385 CALL dtrtrs(
'L',
'T',
'N',int(ncon,mpl),1_mpl,
lapackql(iloff+npar-ncon+1:),int(npar,mpl),&
9387 IF(infolp /= 0) print *,
' DTRTRS failed: ', infolp
9389 CALL dormql(
'L',
'T',int(npar,mpl),1_mpl,int(ncon,mpl),
lapackql(iloff+1:),int(npar,mpl),&
9391 IF(infolp /= 0) print *,
' DORMQL failed: ', infolp
9410 IF (nfit > npar)
THEN
9413 WRITE(
lunlog,*)
'Factorization of global matrix (A->L*D*L^t)'
9417#ifdef SCOREP_USER_ENABLE
9418 scorep_user_region_by_name_begin(
"UR_dsytrf", scorep_user_region_type_common)
9420 CALL dsytrf(
'U',int(nfit,mpl),
globalmatd(imoff+1:),int(nfit,mpl),&
9422#ifdef SCOREP_USER_ENABLE
9423 scorep_user_region_by_name_end(
"UR_dsytrf")
9430 WRITE(
lunlog,*)
'Factorization of global matrix (A->L*L^t)'
9434#ifdef SCOREP_USER_ENABLE
9435 scorep_user_region_by_name_begin(
"UR_dpotrf", scorep_user_region_type_common)
9437 CALL dpotrf(
'U',int(nfit,mpl),
globalmatd(imoff+1:),int(npar,mpl),infolp)
9438#ifdef SCOREP_USER_ENABLE
9439 scorep_user_region_by_name_end(
"UR_dpotrf")
9447 WRITE(lun,*)
'No rank defect of the symmetric matrix'
9449 WRITE(lun,*)
'No rank defect of the symmetric block', ib,
' of size', npar
9453 WRITE(*,*)
'Warning: factorization of the symmetric',nfit, &
9454 '-by-',nfit,
' failed at index ', infolp
9455 WRITE(lun,*)
'Warning: factorization of the symmetric',nfit, &
9456 '-by-',nfit,
' failed at index ', infolp
9457 CALL peend(29,
'Aborted, factorization of global matrix failed')
9458 stop
'mdutrf: bad matrix'
9463 IF (nfit > npar)
THEN
9464 CALL dsytrs(
'U',int(nfit,mpl),1_mpl,
globalmatd(imoff+1:),int(nfit,mpl),&
9466 IF(infolp /= 0) print *,
' DSYTRS failed: ', infolp
9468 CALL dpotrs(
'U',int(nfit,mpl),1_mpl,
globalmatd(imoff+1:),int(npar,mpl),&
9470 IF(infolp /= 0) print *,
' DPOTRS failed: ', infolp
9474 IF(nfit < npar)
THEN
9479 CALL dormql(
'L',
'N',int(npar,mpl),1_mpl,int(ncon,mpl),
lapackql(iloff+1:),int(npar,mpl),&
9481 IF(infolp /= 0) print *,
' DORMQL failed: ', infolp
9488 iloff=iloff+int(npar,mpl)*int(ncon,mpl)
9510 INTEGER(mpi) :: icboff
9511 INTEGER(mpi) :: icblst
9512 INTEGER(mpi) :: icoff
9513 INTEGER(mpi) :: icfrst
9514 INTEGER(mpi) :: iclast
9515 INTEGER(mpi) :: ipfrst
9516 INTEGER(mpi) :: iplast
9517 INTEGER(mpi) :: ipoff
9520 INTEGER(mpi) :: ncon
9521 INTEGER(mpi) :: npar
9523 INTEGER(mpl) :: imoff
9524 INTEGER(mpl) :: iloff
9525 INTEGER(mpi) :: infolp
9526 INTEGER :: nbopt, ILAENV
9528 REAL(mpd),
INTENT(IN) :: a(mszcon)
9529 REAL(mpd),
INTENT(OUT) :: emin
9530 REAL(mpd),
INTENT(OUT) :: emax
9541 iloff=iloff+int(npar,mpl)*int(ncon,mpl)
9560 DO icb=icboff+1,icboff+icblst
9567 lapackql(iloff+ipfrst:iloff+iplast)=a(imoff+1:imoff+npb)
9584 nbopt = ilaenv( 1_mpl,
'DGEQLF',
'', int(npar,mpl), int(ncon,mpl), int(npar,mpl), -1_mpl )
9585 print *,
'LAPACK optimal block size for DGEQLF:', nbopt
9586 lplwrk=int(ncon,mpl)*int(nbopt,mpl)
9589#ifdef SCOREP_USER_ENABLE
9590 scorep_user_region_by_name_begin(
"UR_dgeqlf", scorep_user_region_type_common)
9592 CALL dgeqlf(int(npar,mpl),int(ncon,mpl),
lapackql(iloff+1:),int(npar,mpl),&
9594 IF(infolp /= 0) print *,
' DGEQLF failed: ', infolp
9595#ifdef SCOREP_USER_ENABLE
9596 scorep_user_region_by_name_end(
"UR_dgeqlf")
9600 iloff=iloff+int(npar,mpl)*int(ncon,mpl)
9603 IF(emax < emin)
THEN
9631 INTEGER(mpi) :: icoff
9632 INTEGER(mpi) :: ipoff
9634 INTEGER(mpi) :: ncon
9635 INTEGER(mpi) :: npar
9636 INTEGER(mpl) :: imoff
9637 INTEGER(mpl) :: iloff
9638 INTEGER(mpi) :: infolp
9639 CHARACTER (LEN=1) :: transr, transl
9641 LOGICAL,
INTENT(IN) :: t
9660 IF(ncon <= 0 ) cycle
9663#ifdef SCOREP_USER_ENABLE
9664 scorep_user_region_by_name_begin(
"UR_dormql", scorep_user_region_type_common)
9672 DO i=ipoff+1,ipoff+npar
9678 CALL dormql(
'R',transr,int(npar,mpl),int(npar,mpl),int(ncon,mpl),
lapackql(iloff+1:),&
9681 IF(infolp /= 0) print *,
' DORMQL failed: ', infolp
9683 CALL dormql(
'L',transl,int(npar,mpl),int(npar,mpl),int(ncon,mpl),
lapackql(iloff+1:),&
9686 IF(infolp /= 0) print *,
' DORMQL failed: ', infolp
9687#ifdef SCOREP_USER_ENABLE
9688 scorep_user_region_by_name_end(
"UR_dormql")
9692 iloff=iloff+int(npar,mpl)*int(ncon,mpl)
9698include
'mkl_pardiso.f90'
9729 TYPE(mkl_pardiso_handle) :: pt(64)
9731 INTEGER(mpl),
PARAMETER :: maxfct =1
9732 INTEGER(mpl),
PARAMETER :: mnum = 1
9733 INTEGER(mpl),
PARAMETER :: nrhs = 1
9735 INTEGER(mpl) :: mtype
9736 INTEGER(mpl) :: phase
9737 INTEGER(mpl) :: error
9738 INTEGER(mpl) :: msglvl
9742 INTEGER(mpl) :: idum(1)
9744 INTEGER(mpl) :: length
9745 INTEGER(mpi) :: nfill
9746 INTEGER(mpi) :: npdblk
9747 REAL(mpd) :: adum(1)
9748 REAL(mpd) :: ddum(1)
9750 INTEGER(mpl) :: iparm(64)
9751 REAL(mpd),
ALLOCATABLE :: b( : )
9752 REAL(mpd),
ALLOCATABLE :: x( : )
9766#ifdef SCOREP_USER_ENABLE
9767 scorep_user_region_by_name_begin(
"UR_mspd00", scorep_user_region_type_common)
9774 WRITE(*,*)
'MSPARDISO: number of rows to fill up = ', nfill
9787 CALL pardiso_64(pt, maxfct, mnum, mtype, phase, int(npdblk,mpl), adum, idum, idum, &
9788 idum, nrhs, iparm, msglvl, ddum, ddum, error)
9789 IF (error /= 0)
THEN
9790 WRITE(lun,*)
'The following ERROR was detected: ', error
9791 WRITE(*,
'(A,2I10)')
' PARDISO release failed (phase, error): ', phase, error
9792 IF (
ipddbg == 0)
WRITE(*,*)
' rerun with "debugPARDISO" for more info'
9793 CALL peend(40,
'Aborted, other error: PARDISO release')
9794 stop
'MSPARDISO: stopping due to error in PARDISO release'
9811 IF (iparm(1) == 0)
WRITE(lun,*)
'PARDISO using defaults '
9812 IF (iparm(43) /= 0)
THEN
9813 WRITE(lun,*)
'PARDISO: computation of the diagonal of inverse matrix not implemented !'
9821#ifdef SCOREP_USER_ENABLE
9822 scorep_user_region_by_name_end(
"UR_mspd00")
9830 WRITE(
lunlog,*)
'Decomposition of global matrix (A->L*D*L^t)'
9837#ifdef SCOREP_USER_ENABLE
9838 scorep_user_region_by_name_begin(
"UR_mspd11", scorep_user_region_type_common)
9847 WRITE(lun,*)
' iparm(',i,
') =', iparm(i)
9850 CALL pardiso_64(pt, maxfct, mnum, mtype, phase, int(npdblk,mpl),
globalmatd,
csr3rowoffsets,
csr3columnlist, &
9851 idum, nrhs, iparm, msglvl, ddum, ddum, error)
9852#ifdef SCOREP_USER_ENABLE
9853 scorep_user_region_by_name_end(
"UR_mspd11")
9856 WRITE(lun,*)
'PARDISO reordering completed ... '
9857 WRITE(lun,*)
'PARDISO peak memory required (KB)', iparm(15)
9860 WRITE(lun,*)
' iparm(',i,
') =', iparm(i)
9863 IF (error /= 0)
THEN
9864 WRITE(lun,*)
'The following ERROR was detected: ', error
9865 WRITE(*,
'(A,2I10)')
' PARDISO decomposition failed (phase, error): ', phase, error
9866 IF (
ipddbg == 0)
WRITE(*,*)
' rerun with "debugPARDISO" for more info'
9867 CALL peend(40,
'Aborted, other error: PARDISO reordering')
9868 stop
'MSPARDISO: stopping due to error in PARDISO reordering'
9870 IF (iparm(60) == 0)
THEN
9875 WRITE(lun,*)
'Size (KB) of allocated memory = ',
ipdmem
9876 WRITE(lun,*)
'Number of nonzeros in factors = ',iparm(18)
9877 WRITE(lun,*)
'Number of factorization MFLOPS = ',iparm(19)
9881#ifdef SCOREP_USER_ENABLE
9882 scorep_user_region_by_name_begin(
"UR_mspd22", scorep_user_region_type_common)
9885 CALL pardiso_64(pt, maxfct, mnum, mtype, phase, int(npdblk,mpl),
globalmatd,
csr3rowoffsets,
csr3columnlist, &
9886 idum, nrhs, iparm, msglvl, ddum, ddum, error)
9887#ifdef SCOREP_USER_ENABLE
9888 scorep_user_region_by_name_end(
"UR_mspd22")
9891 WRITE(lun,*)
'PARDISO factorization completed ... '
9894 WRITE(lun,*)
' iparm(',i,
') =', iparm(i)
9897 IF (error /= 0)
THEN
9898 WRITE(lun,*)
'The following ERROR was detected: ', error
9899 WRITE(*,
'(A,2I10)')
' PARDISO decomposition failed (phase, error): ', phase, error
9900 IF (
ipddbg == 0)
WRITE(*,*)
' rerun with "debugPARDISO" for more info'
9901 CALL peend(40,
'Aborted, other error: PARDISO factorization')
9902 stop
'MSPARDISO: stopping due to error in PARDISO factorization'
9905 IF (iparm(14) > 0) &
9906 WRITE(lun,*)
'Number of perturbed pivots = ',iparm(14)
9907 WRITE(lun,*)
'Number of positive eigenvalues = ',iparm(22)-nfill
9908 WRITE(lun,*)
'Number of negative eigenvalues = ',iparm(23)
9909 ELSE IF (iparm(30) > 0)
THEN
9910 WRITE(lun,*)
'Equation with bad pivot (<=0.) = ',iparm(30)
9919 CALL mpalloc(b,length,
' PARDISO r.h.s')
9920 CALL mpalloc(x,length,
' PARDISO solution')
9923#ifdef SCOREP_USER_ENABLE
9924 scorep_user_region_by_name_begin(
"UR_mspd33", scorep_user_region_type_common)
9928 CALL pardiso_64(pt, maxfct, mnum, mtype, phase, int(npdblk,mpl),
globalmatd,
csr3rowoffsets,
csr3columnlist, &
9929 idum, nrhs, iparm, msglvl, b, x, error)
9930#ifdef SCOREP_USER_ENABLE
9931 scorep_user_region_by_name_end(
"UR_mspd33")
9937 WRITE(lun,*)
'PARDISO solve completed ... '
9938 IF (error /= 0)
THEN
9939 WRITE(lun,*)
'The following ERROR was detected: ', error
9940 WRITE(*,
'(A,2I10)')
' PARDISO decomposition failed (phase, error): ', phase, error
9941 IF (
ipddbg == 0)
WRITE(*,*)
' rerun with "debugPARDISO" for more info'
9942 CALL peend(40,
'Aborted, other error: PARDISO solve')
9943 stop
'MSPARDISO: stopping due to error in PARDISO solve'
9957 INTEGER(mpi) :: iast
9958 INTEGER(mpi) :: idia
9959 INTEGER(mpi) :: imin
9960 INTEGER(mpl) :: ioff1
9962 INTEGER(mpi) :: last
9964 INTEGER(mpi) :: nmax
9965 INTEGER(mpi) :: nmin
9966 INTEGER(mpi) :: ntop
9988 WRITE(
lunlog,*)
'Shrinkage of global matrix (A->Q^t*A*Q)'
10025 DO WHILE(ntop < nmax)
10029 CALL hmpdef(7,real(nmin,mps),real(ntop,mps),
'log10 of positive eigenvalues')
10039 iast=max(1,imin-60)
10040 CALL gmpdef(3,2,
'low-value end of eigenvalues')
10043 CALL gmpxy(3,real(i,mps),evalue)
10057 WRITE(lun,*)
'The first (largest) eigenvalues ...'
10060 WRITE(lun,*)
'The last eigenvalues ... up to',last
10064 WRITE(lun,*)
'The eigenvalues from',
nvgb+1,
' to',
nagb
10068 WRITE(lun,*)
'Log10 + 3 of ',
nfgb,
' eigenvalues in decreasing',
' order'
10069 WRITE(lun,*)
'(for Eigenvalue < 0.001 the value 0.0 is shown)'
10072 'printed for negative eigenvalues'
10075 WRITE(lun,*) last,
' significances: insignificant if ', &
10076 'compatible with N(0,1)'
10105 INTEGER(mpl) :: ioff1
10106 INTEGER(mpl) :: ioff2
10142 INTEGER(mpi) :: istop
10143 INTEGER(mpi) :: itn
10144 INTEGER(mpi) :: itnlim
10145 INTEGER(mpi) :: lun
10146 INTEGER(mpi) :: nout
10147 INTEGER(mpi) :: nrkd
10148 INTEGER(mpi) :: nrkd2
10154 REAL(mpd) :: arnorm
10193 WRITE(lun,*)
'MMINRS: PRECONS ended ', nrkd
10197 globalcorrections, itnlim, nout, rtol, istop, itn, anorm, acond, rnorm, arnorm, ynorm)
10198 ELSE IF(
mbandw > 0)
THEN
10204 WRITE(lun,*)
'MMINRS: EQUDECS ended ', nrkd, nrkd2
10208 globalcorrections, itnlim, nout, rtol, istop, itn, anorm, acond, rnorm, arnorm, ynorm)
10211 globalcorrections, itnlim, nout, rtol, istop, itn, anorm, acond, rnorm, arnorm, ynorm)
10225 IF (
istopa == 0) print *,
'MINRES: istop=0, exact solution x=0.'
10240 INTEGER(mpi) :: istop
10241 INTEGER(mpi) :: itn
10242 INTEGER(mpi) :: itnlim
10243 INTEGER(mpi) :: lun
10244 INTEGER(mpi) :: nout
10245 INTEGER(mpi) :: nrkd
10246 INTEGER(mpi) :: nrkd2
10249 REAL(mpd) :: mxxnrm
10250 REAL(mpd) :: trcond
10260 mxxnrm = real(
nagb,mpd)/sqrt(epsilon(mxxnrm))
10262 trcond = 1.0_mpd/epsilon(trcond)
10263 ELSE IF(
mrmode == 2)
THEN
10293 WRITE(lun,*)
'MMINRS: PRECONS ended ', nrkd
10297 itnlim=itnlim, rtol=rtol, maxxnorm=mxxnrm, trancond=trcond, &
10299 ELSE IF(
mbandw > 0)
THEN
10305 WRITE(lun,*)
'MMINRS: EQUDECS ended ', nrkd, nrkd2
10310 itnlim=itnlim, rtol=rtol, maxxnorm=mxxnrm, trancond=trcond, &
10314 itnlim=itnlim, rtol=rtol, maxxnorm=mxxnrm, trancond=trcond, &
10329 IF (
istopa == 3) print *,
'MINRES: istop=0, exact solution x=0.'
10345 INTEGER(mpi),
INTENT(IN) :: n
10346 REAL(mpd),
INTENT(IN) :: x(n)
10347 REAL(mpd),
INTENT(OUT) :: y(n)
10367 INTEGER(mpi),
INTENT(IN) :: n
10368 REAL(mpd),
INTENT(IN) :: x(n)
10369 REAL(mpd),
INTENT(OUT) :: y(n)
10400 REAL(mps) :: concu2
10401 REAL(mps) :: concut
10402 REAL,
DIMENSION(2) :: ta
10405 INTEGER(mpi) :: iact
10406 INTEGER(mpi) :: iagain
10407 INTEGER(mpi) :: idx
10408 INTEGER(mpi) :: info
10410 INTEGER(mpi) :: ipoff
10411 INTEGER(mpi) :: icoff
10412 INTEGER(mpl) :: ioff
10413 INTEGER(mpi) :: itgbi
10414 INTEGER(mpi) :: ivgbi
10415 INTEGER(mpi) :: jcalcm
10417 INTEGER(mpi) :: labelg
10418 INTEGER(mpi) :: litera
10419 INTEGER(mpl) :: lrej
10420 INTEGER(mpi) :: lun
10421 INTEGER(mpi) :: lunp
10422 INTEGER(mpi) :: minf
10423 INTEGER(mpi) :: mrati
10424 INTEGER(mpi) :: nan
10425 INTEGER(mpi) :: ncon
10426 INTEGER(mpi) :: nfaci
10427 INTEGER(mpi) :: nloopsol
10428 INTEGER(mpi) :: npar
10429 INTEGER(mpi) :: nrati
10430 INTEGER(mpl) :: nrej
10431 INTEGER(mpi) :: nsol
10432 INTEGER(mpi) :: inone
10434 INTEGER(mpi) :: infolp
10435 INTEGER(mpi) :: nfit
10436 INTEGER(mpl) :: imoff
10440 REAL(mpd) :: dratio
10441 REAL(mpd) :: dwmean
10450 LOGICAL :: warnerss
10451 LOGICAL :: warners3
10453 CHARACTER (LEN=7) :: cratio
10454 CHARACTER (LEN=7) :: cfacin
10455 CHARACTER (LEN=7) :: crjrat
10466 WRITE(lunp,*)
'Solution algorithm: '
10467 WRITE(lunp,121)
'=================================================== '
10470 WRITE(lunp,121)
'solution method:',
'matrix inversion'
10471 ELSE IF(
metsol == 2)
THEN
10472 WRITE(lunp,121)
'solution method:',
'diagonalization'
10473 ELSE IF(
metsol == 3)
THEN
10474 WRITE(lunp,121)
'solution method:',
'decomposition'
10475 ELSE IF(
metsol == 4)
THEN
10476 WRITE(lunp,121)
'solution method:',
'minres (Paige/Saunders)'
10477 ELSE IF(
metsol == 5)
THEN
10478 WRITE(lunp,121)
'solution method:',
'minres-qlp (Choi/Paige/Saunders)'
10480 WRITE(lunp,121)
' ',
' using QR factorization'
10481 ELSE IF(
mrmode == 2)
THEN
10482 WRITE(lunp,121)
' ',
' using QLP factorization'
10484 WRITE(lunp,121)
' ',
' using QR and QLP factorization'
10485 WRITE(lunp,123)
'transition condition',
mrtcnd
10487 ELSE IF(
metsol == 6)
THEN
10488 WRITE(lunp,121)
'solution method:', &
10489 'gmres (generalized minimzation of residuals)'
10491 ELSE IF(
metsol == 7)
THEN
10493 WRITE(lunp,121)
'solution method:',
'LAPACK factorization (DSPTRF)'
10495 WRITE(lunp,121)
'solution method:',
'LAPACK factorization (DPPTRF)'
10497 IF(
ilperr == 1)
WRITE(lunp,121)
' ',
'with error calculation (D??TRI)'
10498 ELSE IF(
metsol == 8)
THEN
10500 WRITE(lunp,121)
'solution method:',
'LAPACK factorization (DSYTRF)'
10502 WRITE(lunp,121)
'solution method:',
'LAPACK factorization (DPOTRF)'
10504 IF(
ilperr == 1)
WRITE(lunp,121)
' ',
'with error calculation (D??TRI)'
10506 ELSE IF(
metsol == 9)
THEN
10508 WRITE(lunp,121)
'solution method:',
'Intel oneMKL PARDISO (sparse matrix (CSR3))'
10510 WRITE(lunp,121)
'solution method:',
'Intel oneMKL PARDISO (sparse matrix (BSR3))'
10515 WRITE(lunp,123)
'convergence limit at Delta F=',
dflim
10516 WRITE(lunp,122)
'maximum number of iterations=',
mitera
10519 WRITE(lunp,122)
'matrix recalculation up to ',
matrit,
'. iteration'
10523 WRITE(lunp,121)
'matrix storage:',
'full'
10524 ELSE IF(
matsto == 2)
THEN
10525 WRITE(lunp,121)
'matrix storage:',
'sparse'
10527 WRITE(lunp,122)
'pre-con band-width parameter=',
mbandw
10529 WRITE(lunp,121)
'pre-conditioning:',
'default'
10530 ELSE IF(
mbandw < 0)
THEN
10531 WRITE(lunp,121)
'pre-conditioning:',
'none!'
10532 ELSE IF(
mbandw > 0)
THEN
10534 WRITE(lunp,121)
'pre-conditioning=',
'skyline-matrix (rank preserving)'
10536 WRITE(lunp,121)
'pre-conditioning=',
'band-matrix'
10541 WRITE(lunp,121)
'using pre-sigmas:',
'no'
10544 WRITE(lunp,124)
'pre-sigmas defined for', &
10545 REAL(100*npresg,mps)/
REAL(nvgb,mps),
' % of variable parameters'
10546 WRITE(lunp,123)
'default pre-sigma=',
regpre
10549 WRITE(lunp,121)
'regularization:',
'no'
10551 WRITE(lunp,121)
'regularization:',
'yes'
10552 WRITE(lunp,123)
'regularization factor=',
regula
10556 WRITE(lunp,121)
'Chi square cut equiv 3 st.dev applied'
10557 WRITE(lunp,123)
'... in first iteration with factor',
chicut
10558 WRITE(lunp,123)
'... in second iteration with factor',
chirem
10559 WRITE(lunp,121)
' (reduced by sqrt in next iterations)'
10562 WRITE(lunp,121)
'Scaling of measurement errors applied'
10563 WRITE(lunp,123)
'... factor for "global" measuements',
dscerr(1)
10564 WRITE(lunp,123)
'... factor for "local" measuements',
dscerr(2)
10567 WRITE(lunp,122)
'Down-weighting of outliers in',
lhuber,
' iterations'
10568 WRITE(lunp,123)
'Cut on downweight fraction',
dwcut
10572121
FORMAT(1x,a40,3x,a)
10573122
FORMAT(1x,a40,3x,i0,a)
10574123
FORMAT(1x,a40,2x,e9.2)
10575124
FORMAT(1x,a40,3x,f5.1,a)
10585 stepl =real(stp,mps)
10597 ELSE IF(
metsol == 2)
THEN
10600 ELSE IF(
metsol == 3)
THEN
10603 ELSE IF(
metsol == 4)
THEN
10606 ELSE IF(
metsol == 5)
THEN
10609 ELSE IF(
metsol == 6)
THEN
10621 WRITE(
lunlog,*)
'Checking feasibility of parameters:'
10622 WRITE(*,*)
'Checking feasibility of parameters:'
10623 CALL feasib(concut,iact)
10625 WRITE(*,102) concut
10626 WRITE(*,*)
' parameters are made feasible'
10627 WRITE(
lunlog,102) concut
10628 WRITE(
lunlog,*)
' parameters are made feasible'
10630 WRITE(*,*)
' parameters are feasible (i.e. satisfy constraints)'
10631 WRITE(
lunlog,*)
' parameters are feasible (i.e. satisfy constraints)'
10639 WRITE(*,*)
'Reading files and accumulating vectors/matrices ...'
10643 WRITE(
lunlog,*)
'Reading files and accumulating vectors/matrices ...'
10660 IF(jcalcm+1 /= 0)
THEN
10669 IF(
iterat /= litera)
THEN
10690 IF (iabs(jcalcm) <= 1)
THEN
10703 IF(3*nrej >
nrecal)
THEN
10705 WRITE(*,*)
'Data records rejected in previous loop: '
10707 WRITE(*,*)
'Too many rejects (>33.3%) - stop'
10708 CALL peend(26,
'Aborted, too many rejects')
10715 IF(abs(
icalcm) == 1)
THEN
10727 ELSE IF(
metsol == 2)
THEN
10729 ELSE IF(
metsol == 3)
THEN
10731 ELSE IF(
metsol == 4)
THEN
10733 ELSE IF(
metsol == 5)
THEN
10735 ELSE IF(
metsol == 6)
THEN
10736 WRITE(*,*)
'... reserved for GMRES (not yet!)'
10739 ELSE IF(
metsol == 7)
THEN
10741 ELSE IF(
metsol == 8)
THEN
10744 ELSE IF(
metsol == 9)
THEN
10758 CALL feasib(concut,iact)
10770 angras=real(db/sqrt(db1*db2),mps)
10771 dbsig=16.0_mpd*sqrt(max(db1,db2))*epsilon(db)
10777 lsflag=lsflag .AND. (db > dbsig)
10784 ELSE IF(
metsol == 2)
THEN
10787 ELSE IF(
metsol == 3)
THEN
10790 ELSE IF(
metsol == 4)
THEN
10793 ELSE IF(
metsol == 5)
THEN
10796 ELSE IF(
metsol == 6)
THEN
10806 IF(db <= -dbsig)
THEN
10807 WRITE(*,*)
'Function not decreasing:',db
10808 IF(db > -1.0e-3_mpd)
THEN
10810 IF (iagain <= 1)
THEN
10811 WRITE(*,*)
'... again matrix calculation'
10815 WRITE(*,*)
'... aborting iterations'
10819 WRITE(*,*)
'... stopping iterations'
10842 IF (
nloopn == nloopsol)
THEN
10848 stepl=real(stp,mps)
10858 WRITE(*,*)
'Result vector containes ', nan,
' NaNs - stop'
10859 CALL peend(25,
'Aborted, result vector contains NaNs')
10866 WRITE(*,*)
'Subito! Exit after first step.'
10871 WRITE(*,*)
'INFO=0 should not happen (line search input err)'
10872 IF (iagain <= 0)
THEN
10877 IF(info < 0 .OR.
nloopn == nloopsol) cycle
10881 CALL feasib(concut,iact)
10882 IF(iact /= 0.OR.
chicut > 1.0)
THEN
10897 IF(sum(
nrejec) /= 0)
THEN
10899 WRITE(*,*)
'Data records rejected in last loop: '
10921 nfit=npar+ncon;
IF (
icelim > 0) nfit=npar-ncon
10922 IF (nfit > npar)
THEN
10925 WRITE(
lunlog,*)
'Inverse of global matrix from LDLt factorization'
10930#ifdef SCOREP_USER_ENABLE
10931 scorep_user_region_by_name_begin(
"UR_dsptri", scorep_user_region_type_common)
10934 IF(infolp /= 0) print *,
' DSPTRI failed: ', infolp
10935#ifdef SCOREP_USER_ENABLE
10936 scorep_user_region_by_name_end(
"UR_dsptri")
10942#ifdef SCOREP_USER_ENABLE
10943 scorep_user_region_by_name_begin(
"UR_dsytri", scorep_user_region_type_common)
10945 CALL dsytri(
'U',int(nfit,mpl),
globalmatd(imoff+1:),int(nfit,mpl),&
10947 IF(infolp /= 0) print *,
' DSYTRI failed: ', infolp
10948#ifdef SCOREP_USER_ENABLE
10949 scorep_user_region_by_name_end(
"UR_dsytri")
10956 WRITE(
lunlog,*)
'Inverse of global matrix from LLt factorization'
10961#ifdef SCOREP_USER_ENABLE
10962 scorep_user_region_by_name_begin(
"UR_dpptri", scorep_user_region_type_common)
10964 CALL dpptri(
'U',int(nfit,mpl),
globalmatd(imoff+1:),infolp)
10965 IF(infolp /= 0) print *,
' DPPTRI failed: ', infolp
10966#ifdef SCOREP_USER_ENABLE
10967 scorep_user_region_by_name_end(
"UR_dpptri")
10972#ifdef SCOREP_USER_ENABLE
10973 scorep_user_region_by_name_begin(
"UR_dpotri", scorep_user_region_type_common)
10975 CALL dpotri(
'U',int(nfit,mpl),
globalmatd(imoff+1:),int(npar,mpl),infolp)
10976 IF(infolp /= 0) print *,
' DPOTRI failed: ', infolp
10977#ifdef SCOREP_USER_ENABLE
10978 scorep_user_region_by_name_end(
"UR_dpotri")
10996 DO i=npar-ncon+1,npar
11003 WRITE(
lunlog,*)
'Expansion of global matrix (A->Q*A*Q^t)'
11019 catio=real(dratio,mps)
11023 mrati=nint(100.0*catio,mpi)
11027 IF (
nfilw <= 0)
THEN
11028 WRITE(lunp,*)
'Sum(Chi^2)/Sum(Ndf) =',
fvalue
11030 WRITE(lunp,*)
' =',dratio
11032 WRITE(lunp,*)
'Sum(W*Chi^2)/Sum(Ndf)/<W> =',
fvalue
11034 WRITE(lunp,*)
' /',dwmean
11035 WRITE(lunp,*)
' =',dratio
11039 ' with correction for down-weighting ',catio
11060 nrati=nint(10000.0*real(nrej,mps)/real(
nrecal,mps),mpi)
11061 WRITE(crjrat,197) 0.01_mpd*real(nrati,mpd)
11062 nfaci=nint(100.0*sqrt(catio),mpi)
11064 WRITE(cratio,197) 0.01_mpd*real(mrati,mpd)
11065 WRITE(cfacin,197) 0.01_mpd*real(nfaci,mpd)
11068 IF(mrati < 90.OR.mrati > 110) warner=.true.
11069 IF(nrati > 100) warner=.true.
11070 IF(
ncgbe /= 0) warner=.true.
11072 IF(
nalow /= 0) warners=.true.
11074 IF(
nmiss1 /= 0) warnerss=.true.
11075 IF(iagain /= 0) warnerss=.true.
11076 IF(
ndefec /= 0) warnerss=.true.
11077 IF(
ndefpg /= 0) warnerss=.true.
11079 IF(
nrderr /= 0) warners3=.true.
11081 IF(warner.OR.warners.OR.warnerss.Or.warners3)
THEN
11084 WRITE(*,199)
'WarningWarningWarningWarningWarningWarningWarningWarningWar'
11085 WRITE(*,199)
'arningWarningWarningWarningWarningWarningWarningWarningWarn'
11086 WRITE(*,199)
'rningWarningWarningWarningWarningWarningWarningWarningWarni'
11087 WRITE(*,199)
'ningWarningWarningWarningWarningWarningWarningWarningWarnin'
11088 WRITE(*,199)
'ingWarningWarningWarningWarningWarningWarningWarningWarning'
11089 WRITE(*,199)
'ngWarningWarningWarningWarningWarningWarningWarningWarningW'
11090 WRITE(*,199)
'gWarningWarningWarningWarningWarningWarningWarningWarningWa'
11092 IF(mrati < 90.OR.mrati > 110)
THEN
11094 WRITE(*,*)
' Chi^2/Ndf = ',cratio,
' (should be close to 1)'
11095 WRITE(*,*)
' => multiply all input standard ', &
11096 'deviations by factor',cfacin
11099 IF(nrati > 100)
THEN
11101 WRITE(*,*)
' Fraction of rejects =',crjrat,
' %', &
11102 ' (should be far below 1 %)'
11103 WRITE(*,*)
' => please provide correct mille data'
11107 IF(iagain /= 0)
THEN
11109 WRITE(*,*)
' Matrix not positiv definite '// &
11110 '(function not decreasing)'
11111 WRITE(*,*)
' => please provide correct mille data'
11116 WRITE(*,*)
' Rank defect =',
ndefec, &
11117 ' for global matrix, should be 0'
11118 WRITE(*,*)
' => please provide correct mille data'
11123 WRITE(*,*)
' Rank defect for',
ndefpg, &
11124 ' parameter groups, should be 0'
11125 WRITE(*,*)
' => please provide correct mille data'
11130 WRITE(*,*)
' Rank defect =',
nmiss1, &
11131 ' for constraint equations, should be 0'
11132 WRITE(*,*)
' => please correct constraint definition'
11135 IF(
ncgbe /= 0)
THEN
11137 WRITE(*,*)
' Number of empty constraints =',
ncgbe,
', should be 0'
11138 WRITE(*,*)
' => please check constraint definition, mille data'
11141 IF(
nxlow /= 0)
THEN
11143 WRITE(*,*)
' Possible rank defects =',
nxlow,
' for global matrix'
11144 WRITE(*,*)
' (too few accepted entries)'
11145 WRITE(*,*)
' => please check mille data and ENTRIES cut'
11148 IF(
nalow /= 0)
THEN
11150 WRITE(*,*)
' Possible bad elements =',
nalow,
' in global vector'
11151 WRITE(*,*)
' (toos few accepted entries)'
11152 IF(
ipcntr > 0)
WRITE(*,*)
' (indicated in millepede.res by counts<0)'
11153 WRITE(*,*)
' => please check mille data and ENTRIES cut'
11158 WRITE(*,*)
' Binary file(s) with read errors =',
nrderr,
' (treated as EOF)'
11159 WRITE(*,*)
' => please check mille data'
11163 WRITE(*,199)
'WarningWarningWarningWarningWarningWarningWarningWarningWar'
11164 WRITE(*,199)
'arningWarningWarningWarningWarningWarningWarningWarningWarn'
11165 WRITE(*,199)
'rningWarningWarningWarningWarningWarningWarningWarningWarni'
11166 WRITE(*,199)
'ningWarningWarningWarningWarningWarningWarningWarningWarnin'
11167 WRITE(*,199)
'ingWarningWarningWarningWarningWarningWarningWarningWarning'
11168 WRITE(*,199)
'ngWarningWarningWarningWarningWarningWarningWarningWarningW'
11169 WRITE(*,199)
'gWarningWarningWarningWarningWarningWarningWarningWarningWa'
11180 ELSE IF(
metsol == 2)
THEN
11182 ELSE IF(
metsol == 3)
THEN
11188 IF(labelg == 0) cycle
11189 itgbi=inone(labelg)
11192 IF(ivgbi < 0) ivgbi=0
11193 IF(ivgbi == 0) cycle
11202 ELSE IF(
metsol == 6)
THEN
11205 ELSE IF(
metsol == 7)
THEN
11213 CALL peend(4,
'Ended with severe warnings (bad binary file(s))')
11214 ELSE IF (warnerss)
THEN
11215 CALL peend(3,
'Ended with severe warnings (bad global matrix)')
11216 ELSE IF (warners)
THEN
11217 CALL peend(2,
'Ended with severe warnings (insufficient measurements)')
11218 ELSE IF (warner)
THEN
11219 CALL peend(1,
'Ended with warnings (bad measurements)')
11221 CALL peend(0,
'Ended normally')
11224102
FORMAT(
' Call FEASIB with cut=',g10.3)
11242 INTEGER(mpi) :: kfl
11243 INTEGER(mpi) :: kmin
11244 INTEGER(mpi) :: kmax
11245 INTEGER(mpi) :: nrc
11246 INTEGER(mpl) :: nrej
11252 REAL(mpd) :: sumallw
11253 REAL(mpd) :: sumrejw
11255 sumallw=0.; sumrejw=0.;
11264 sumallw=sumallw+real(nrc,mpd)*
wfd(kfl)
11265 sumrejw=sumrejw+real(nrej,mpd)*
wfd(kfl)
11266 frac=real(nrej,mps)/real(nrc,mps)
11267 IF (frac > fmax)
THEN
11271 IF (frac < fmin)
THEN
11278 WRITE(*,
"(' Weighted fraction =',F8.2,' %')") 100.*sumrejw/sumallw
11279 IF (
nfilb > 1)
THEN
11280 WRITE(*,
"(' File with max. fraction ',I6,' :',F8.2,' %')") kmax, 100.*fmax
11281 WRITE(*,
"(' File with min. fraction ',I6,' :',F8.2,' %')") kmin, 100.*fmin
11307 INTEGER(mpi) :: iargc
11310 INTEGER(mpi) :: ierrf
11311 INTEGER(mpi) :: ieq
11312 INTEGER(mpi) :: ifilb
11313 INTEGER(mpi) :: ioff
11314 INTEGER(mpi) :: iopt
11315 INTEGER(mpi) :: ios
11316 INTEGER(mpi) :: iosum
11319 INTEGER(mpi) :: mat
11320 INTEGER(mpi) :: nab
11321 INTEGER(mpi) :: nline
11322 INTEGER(mpi) :: npat
11323 INTEGER(mpi) :: ntext
11325 INTEGER(mpi) :: nuf
11326 INTEGER(mpi) :: nums
11327 INTEGER(mpi) :: nufile
11328 INTEGER(mpi) :: lenfileInfo
11329 INTEGER(mpi) :: lenFileNames
11330 INTEGER(mpi) :: matint
11331 INTEGER(mpi),
DIMENSION(:,:),
ALLOCATABLE :: vecfileInfo
11332 INTEGER(mpi),
DIMENSION(:,:),
ALLOCATABLE :: tempArray
11333 INTEGER(mpl) :: rows
11334 INTEGER(mpl) :: cols
11335 INTEGER(mpl) :: newcols
11336 INTEGER(mpl) :: length
11338 CHARACTER (LEN=1024) :: text
11339 CHARACTER (LEN=1024) :: fname
11340 CHARACTER (LEN=14) :: bite(3)
11341 CHARACTER (LEN=32) :: keystx
11342 INTEGER(mpi),
PARAMETER :: mnum=100
11343 REAL(mpd) :: dnum(mnum)
11347 SUBROUTINE initc(nfiles)
BIND(c)
11349 INTEGER(c_int),
INTENT(IN),
VALUE :: nfiles
11350 END SUBROUTINE initc
11355 DATA bite/
'C_binary',
'text ',
'Fortran_binary'/
11370 WRITE(*,*)
'Command line options: '
11371 WRITE(*,*)
'--------------------- '
11373 CALL getarg(i,text)
11374 CALL rltext(text,ia,ib,nab)
11375 WRITE(*,101) i,text(1:nab)
11376 IF(text(ia:ia) /=
'-')
THEN
11377 nu=nufile(text(ia:ib))
11380 WRITE(*,*)
'Second text file in command line - stop'
11381 CALL peend(12,
'Aborted, second text file in command line')
11387 WRITE(*,*)
'Open error for file:',text(ia:ib),
' - stop'
11388 CALL peend(16,
'Aborted, open error for file')
11389 IF(text(ia:ia) /=
'/')
THEN
11390 CALL getenv(
'PWD',text)
11391 CALL rltext(text,ia,ib,nab)
11392 WRITE(*,*)
'PWD:',text(ia:ib)
11397 IF(index(text(ia:ib),
'b') /= 0)
THEN
11399 WRITE(*,*)
'Debugging requested'
11401 it=index(text(ia:ib),
't')
11404 ieq=index(text(ia+it:ib),
'=')+it
11405 IF (it /= ieq)
THEN
11406 IF (index(text(ia+ieq:ib),
'SL0' ) /= 0)
ictest=2
11407 IF (index(text(ia+ieq:ib),
'SLE' ) /= 0)
ictest=3
11408 IF (index(text(ia+ieq:ib),
'BP' ) /= 0)
ictest=4
11409 IF (index(text(ia+ieq:ib),
'BRLF') /= 0)
ictest=5
11410 IF (index(text(ia+ieq:ib),
'BRLC') /= 0)
ictest=6
11413 IF(index(text(ia:ib),
's') /= 0)
isubit=1
11414 IF(index(text(ia:ib),
'f') /= 0)
iforce=1
11415 IF(index(text(ia:ib),
'c') /= 0)
icheck=1
11416 IF(index(text(ia:ib),
'C') /= 0)
icheck=2
11418 IF(i == iargc())
WRITE(*,*)
'--------------------- '
11439 CALL rltext(text,ia,ib,nab)
11440 nu=nufile(text(ia:ib))
11444 CALL peend(10,
'Aborted, no steering file')
11445 stop
'in FILETC: no steering file. .'
11458 WRITE(*,*)
'Listing of steering file: ',
filnam(1:
nfnam)
11459 WRITE(*,*)
'-------------------------'
11462 WRITE(*,*)
'Open error for steering file - stop'
11463 CALL peend(11,
'Aborted, open error for steering file')
11464 IF(
filnam(1:1) /=
'/')
THEN
11465 CALL getenv(
'PWD',text)
11466 CALL rltext(text,ia,ib,nab)
11467 WRITE(*,*)
'PWD:',text(ia:ib)
11476 rows=6; cols=lenfileinfo
11477 CALL mpalloc(vecfileinfo,rows,cols,
'file info from steering')
11480 READ(10,102,iostat=ierrf) text
11481 IF (ierrf < 0)
EXIT
11482 CALL rltext(text,ia,ib,nab)
11484 IF(nline <= 50)
THEN
11485 WRITE(*,101) nline,text(1:nab)
11486 IF(nline == 50)
WRITE(*,*)
' ...'
11490 CALL rltext(text,ia,ib,nab)
11491 IF(ib == ia+2)
THEN
11492 mat=matint(text(ia:ib),
'end',npat,ntext)
11493 IF(mat == max(npat,ntext))
THEN
11496 WRITE(*,*)
' end-statement after',nline,
' text lines'
11501 keystx=
'fortranfiles'
11502 mat=matint(text(ia:ib),keystx,npat,ntext)
11503 IF(mat == max(npat,ntext))
THEN
11510 mat=matint(text(ia:ib),keystx,npat,ntext)
11511 IF(mat == max(npat,ntext))
THEN
11517 keystx=
'closeandreopen'
11518 mat=matint(text(ia:ib),keystx,npat,ntext)
11519 IF(mat == max(npat,ntext))
THEN
11527 iopt=index(text(ia:ib),
' -- ')
11528 IF (iopt > 0) ie=iopt-1
11531 nu=nufile(text(ia:ie))
11533 IF (
nfiles == lenfileinfo)
THEN
11534 CALL mpalloc(temparray,rows,cols,
'temp file info from steering')
11535 temparray=vecfileinfo
11537 lenfileinfo=lenfileinfo*2
11538 newcols=lenfileinfo
11539 CALL mpalloc(vecfileinfo,rows,newcols,
'file info from steering')
11540 vecfileinfo(:,1:cols)=temparray(:,1:cols)
11546 lenfilenames=lenfilenames+ie-ia+1
11547 vecfileinfo(1,
nfiles)=nline
11548 vecfileinfo(2,
nfiles)=nu
11549 vecfileinfo(3,
nfiles)=ia
11550 vecfileinfo(4,
nfiles)=ie
11551 vecfileinfo(5,
nfiles)=iopt
11552 vecfileinfo(6,
nfiles)=ib
11562 CALL mpalloc(
nfd,length,
'file line (in steering)')
11565 length=lenfilenames
11571 READ(10,102,iostat=ierrf) text
11572 IF (ierrf < 0)
EXIT
11574 IF (nline == vecfileinfo(1,i))
THEN
11575 nfd(i)=vecfileinfo(1,i)
11576 mfd(i)=vecfileinfo(2,i)
11577 ia=vecfileinfo(3,i)-1
11578 lfd(i)=vecfileinfo(4,i)-ia
11580 tfd(ioff+k)=text(ia+k:ia+k)
11585 IF (vecfileinfo(5,i) > 0)
THEN
11586 CALL ratext(text(vecfileinfo(5,i)+4:vecfileinfo(6,i)),nums,dnum,mnum)
11587 IF (nums > 0)
ofd(i)=real(dnum(1),mps)
11597 CALL mpalloc(
ifd,length,
'integrated record numbers (=offset)')
11598 CALL mpalloc(
jfd,length,
'number of accepted records')
11599 CALL mpalloc(
kfd,rows,length,
'number of records in file, file order')
11604 CALL mpalloc(
sfd,rows,length,
'start, end of file name in TFD')
11605 CALL mpalloc(
yfd,length,
'modification date')
11608 WRITE(*,*)
'-------------------------'
11614 WRITE(*,*)
'Table of files:'
11615 WRITE(*,*)
'---------------'
11618 WRITE(8,*)
'Text and data files:'
11622 fname(k:k)=
tfd(ioff+k)
11625 IF (
mprint > 1)
WRITE(*,103) i,bite(
mfd(i)),fname(1:
lfd(i))
11626 WRITE(8,103) i,bite(
mfd(i)),fname(1:
lfd(i))
11630 WRITE(*,*)
'---------------'
11644 IF(
mfd(i) == 3)
THEN
11668 IF(
mfd(i) == 1)
THEN
11689 WRITE(*,*)
'Opening of C-files not supported.'
11705 IF(iosum /= 0)
THEN
11706 CALL peend(15,
'Aborted, open error(s) for binary files')
11707 stop
'FILETC: open error '
11709 IF(
nfilb == 0)
THEN
11710 CALL peend(14,
'Aborted, no binary files')
11711 stop
'FILETC: no binary files '
11714 WRITE(*,*)
nfilb,
' binary files opened'
11716 WRITE(*,*)
nfilb,
' binary files opened and closed'
11720103
FORMAT(i3,2x,a14,3x,a)
11783 INTEGER(mpi) :: ierrf
11784 INTEGER(mpi) :: ioff
11785 INTEGER(mpi) :: ios
11786 INTEGER(mpi) :: iosum
11788 INTEGER(mpi) :: mat
11789 INTEGER(mpi) :: nab
11790 INTEGER(mpi) :: nfiln
11791 INTEGER(mpi) :: nline
11792 INTEGER(mpi) :: nlinmx
11793 INTEGER(mpi) :: npat
11794 INTEGER(mpi) :: ntext
11795 INTEGER(mpi) :: matint
11799 CHARACTER (LEN=1024) :: text
11800 CHARACTER (LEN=1024) :: fname
11803 WRITE(*,*)
'Processing text files ...'
11816 IF(
mfd(i) /= 2) cycle
11818 fname(k:k)=
tfd(ia+k)
11820 WRITE(*,*)
'File ',fname(1:
lfd(i))
11821 IF (
mprint > 1)
WRITE(*,*)
' '
11822 OPEN(10,file=fname(1:
lfd(i)),iostat=ios,form=
'FORMATTED')
11824 WRITE(*,*)
'Open error for file ',fname(1:
lfd(i))
11834 READ(10,102,iostat=ierrf) text
11835 IF (ierrf < 0)
THEN
11838 WRITE(*,*)
' end-of-file after',nline,
' text lines'
11842 IF(nline <= nlinmx.AND.
mprint > 1)
THEN
11843 CALL rltext(text,ia,ib,nab)
11845 WRITE(*,101) nline,text(1:nab)
11846 IF(nline == nlinmx)
WRITE(*,*)
' ...'
11849 CALL rltext(text,ia,ib,nab)
11850 IF(ib == ia+2)
THEN
11851 mat=matint(text(ia:ib),
'end',npat,ntext)
11852 IF(mat == max(npat,ntext))
THEN
11855 WRITE(*,*)
' end-statement after',nline,
' text lines'
11861 IF(nfiln <=
nfiles)
THEN
11862 IF(nline ==
nfd(nfiln))
THEN
11877 IF(iosum /= 0)
THEN
11878 CALL peend(16,
'Aborted, open error(s) for text files')
11879 stop
'FILETX: open error(s) in text files '
11882 WRITE(*,*)
'... end of text file processing.'
11887 WRITE(*,*)
lunkno,
' unknown keywords in steering files, ', &
11888 'or file non-existing,'
11889 WRITE(*,*)
' see above!'
11890 WRITE(*,*)
'------------> stop'
11892 CALL peend(13,
'Aborted, unknown keywords in steering file')
11901 ELSE IF(
matsto == 1)
THEN
11903 ELSE IF(
matsto == 2)
THEN
11906 ELSE IF(
metsol == 1)
THEN
11908 ELSE IF(
metsol == 2)
THEN
11910 ELSE IF(
metsol == 3)
THEN
11912 ELSE IF(
metsol == 4)
THEN
11914 ELSE IF(
metsol == 5)
THEN
11916 ELSE IF(
metsol == 6)
THEN
11919 ELSE IF(
metsol == 7)
THEN
11921 ELSE IF(
metsol == 8)
THEN
11924 ELSE IF(
metsol == 9)
THEN
11929 WRITE(*,*)
'MINRES forced with sparse matrix!'
11931 WRITE(*,*)
'MINRES forced with sparse matrix!'
11933 WRITE(*,*)
'MINRES forced with sparse matrix!'
11938 WRITE(*,*)
'MINRES forced with sparse matrix!'
11940 WRITE(*,*)
'MINRES forced with sparse matrix!'
11942 WRITE(*,*)
'MINRES forced with sparse matrix!'
11950 WRITE(*,*)
'Solution method and matrix-storage mode:'
11952 WRITE(*,*)
' METSOL = 1: matrix inversion'
11953 ELSE IF(
metsol == 2)
THEN
11954 WRITE(*,*)
' METSOL = 2: diagonalization'
11955 ELSE IF(
metsol == 3)
THEN
11956 WRITE(*,*)
' METSOL = 3: decomposition'
11957 ELSE IF(
metsol == 4)
THEN
11958 WRITE(*,*)
' METSOL = 4: MINRES'
11959 ELSE IF(
metsol == 5)
THEN
11960 WRITE(*,*)
' METSOL = 5: MINRES-QLP'
11961 ELSE IF(
metsol == 6)
THEN
11962 WRITE(*,*)
' METSOL = 6: GMRES (-> MINRES)'
11964 ELSE IF(
metsol == 7)
THEN
11965 WRITE(*,*)
' METSOL = 7: LAPACK factorization'
11966 ELSE IF(
metsol == 8)
THEN
11967 WRITE(*,*)
' METSOL = 8: LAPACK factorization'
11969 ELSE IF(
metsol == 9)
THEN
11970 WRITE(*,*)
' METSOL = 9: Intel oneMKL PARDISO'
11975 WRITE(*,*)
' with',
mitera,
' iterations'
11978 WRITE(*,*)
' MATSTO = 0: unpacked symmetric matrix, ',
'n*n elements'
11979 ELSEIF(
matsto == 1)
THEN
11980 WRITE(*,*)
' MATSTO = 1: full symmetric matrix, ',
'(n*n+n)/2 elements'
11981 ELSE IF(
matsto == 2)
THEN
11982 WRITE(*,*)
' MATSTO = 2: sparse matrix (custom)'
11983 ELSE IF(
matsto == 3)
THEN
11985 WRITE(*,*)
' MATSTO = 3: sparse matrix (upper triangle, CSR3)'
11987 WRITE(*,*)
' MATSTO = 3: sparse matrix (upper triangle, BSR3)'
11991 WRITE(*,*)
' and band matrix, width',
mbandw
11995 WRITE(*,*)
'Chi square cut equiv 3 st.dev applied ...'
11996 WRITE(*,*)
' in first iteration with factor',
chicut
11997 WRITE(*,*)
' in second iteration with factor',
chirem
11998 WRITE(*,*)
' (reduced by sqrt in next iterations)'
12002 WRITE(*,*)
' Down-weighting of outliers in',
lhuber,
' iterations'
12003 WRITE(*,*)
' Cut on downweight fraction',
dwcut
12006 WRITE(*,*)
'Iterations (solutions) with line search:'
12015 WRITE(*,*)
' All with Chi square cut scaling factor <= 1.'
12046 INTEGER(mpi) :: ios
12050 INTEGER(mpi) :: npat
12051 INTEGER(mpi) :: ntext
12052 INTEGER(mpi) :: nuprae
12055 CHARACTER (LEN=*),
INTENT(INOUT) :: fname
12061 IF(len(fname) > 5)
THEN
12062 IF(fname(1:5) ==
'rfio:') nuprae=1
12063 IF(fname(1:5) ==
'dcap:') nuprae=2
12064 IF(fname(1:5) ==
'root:') nuprae=3
12066 IF(nuprae == 0)
THEN
12067 INQUIRE(file=fname,iostat=ios,exist=ex)
12068 IF(ios /= 0)
nufile=-abs(ios)
12069 IF(ios /= 0)
RETURN
12070 ELSE IF(nuprae == 1)
THEN
12083 nm=
matint(
'xt',fname(l1:ll),npat,ntext)
12086 nm=
matint(
'tx',fname(l1:ll),npat,ntext)
12107 INTEGER(mpi) :: ier
12108 INTEGER(mpi) :: iomp
12111 INTEGER(mpi) :: kkey
12112 INTEGER(mpi) :: label
12113 INTEGER(mpi) :: lkey
12114 INTEGER(mpi) :: mat
12115 INTEGER(mpi) :: miter
12116 INTEGER(mpi) :: nab
12117 INTEGER(mpi) :: nkey
12118 INTEGER(mpi) :: nkeys
12120 INTEGER(mpi) :: nmeth
12121 INTEGER(mpi) :: npat
12122 INTEGER(mpi) :: ntext
12123 INTEGER(mpi) :: nums
12124 INTEGER(mpi) :: matint
12126 CHARACTER (LEN=*),
INTENT(IN) :: text
12127 INTEGER(mpi),
INTENT(IN) :: nline
12131 parameter(nkeys=7,nmeth=10)
12133 parameter(nkeys=6,nmeth=9)
12136 parameter(nkeys=6,nmeth=7)
12138 CHARACTER (LEN=16) :: methxt(nmeth)
12139 CHARACTER (LEN=16) :: keylst(nkeys)
12140 CHARACTER (LEN=32) :: keywrd
12141 CHARACTER (LEN=32) :: keystx
12142 CHARACTER (LEN=itemCLen) :: ctext
12143 INTEGER(mpi),
PARAMETER :: mnum=100
12144 REAL(mpd) :: dnum(mnum)
12147 INTEGER(mpi) :: ipvs
12150 INTEGER(mpi) :: lpvs
12154 SUBROUTINE additem(length,list,label,value)
12156 INTEGER(mpi),
INTENT(IN OUT) :: length
12157 TYPE(listitem),
DIMENSION(:),
INTENT(IN OUT),
ALLOCATABLE :: list
12158 INTEGER(mpi),
INTENT(IN) :: label
12159 REAL(mpd),
INTENT(IN) :: value
12161 SUBROUTINE additemc(length,list,label,text)
12163 INTEGER(mpi),
INTENT(IN OUT) :: length
12164 TYPE(listitemc),
DIMENSION(:),
INTENT(IN OUT),
ALLOCATABLE :: list
12165 INTEGER(mpi),
INTENT(IN) :: label
12166 CHARACTER(LEN = itemCLen),
INTENT(IN) :: text
12168 SUBROUTINE additemi(length,list,label,ivalue)
12170 INTEGER(mpi),
INTENT(IN OUT) :: length
12171 TYPE(listitemi),
DIMENSION(:),
INTENT(IN OUT),
ALLOCATABLE :: list
12172 INTEGER(mpi),
INTENT(IN) :: label
12173 INTEGER(mpi),
INTENT(IN) :: ivalue
12180 DATA keylst/
'unknown',
'parameter',
'constraint',
'measurement',
'method',
'comment',
'pardiso'/
12181 DATA methxt/
'diagonalization',
'inversion',
'fullMINRES',
'sparseMINRES', &
12182 'fullMINRES-QLP',
'sparseMINRES-QLP',
'decomposition',
'fullLAPACK',
'unpackedLAPACK', &
12185 DATA keylst/
'unknown',
'parameter',
'constraint',
'measurement',
'method',
'comment'/
12186 DATA methxt/
'diagonalization',
'inversion',
'fullMINRES',
'sparseMINRES', &
12187 'fullMINRES-QLP',
'sparseMINRES-QLP',
'decomposition',
'fullLAPACK',
'unpackedLAPACK'/
12190 DATA keylst/
'unknown',
'parameter',
'constraint',
'measurement',
'method',
'comment'/
12191 DATA methxt/
'diagonalization',
'inversion',
'fullMINRES',
'sparseMINRES', &
12192 'fullMINRES-QLP',
'sparseMINRES-QLP',
'decomposition'/
12198 CALL rltext(text,ia,ib,nab)
12199 IF(nab == 0)
GOTO 10
12200 CALL ratext(text(1:nab),nums,dnum,mnum)
12202 IF(nums /= 0) nkey=0
12210 keystx=keylst(nkey)
12211 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12212 IF(100*mat >= 80*max(npat,ntext))
GO TO 10
12218 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12219 IF(100*mat >= 80*max(npat,ntext))
THEN
12221 IF(nums > 0) mprint=nint(dnum(1),mpi)
12226 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12227 IF(100*mat >= 80*max(npat,ntext))
THEN
12230 IF(nums > 0) mdebug=nint(dnum(1),mpi)
12231 IF(nums > 1) mdebg2=nint(dnum(2),mpi)
12236 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12237 IF(100*mat >= 80*max(npat,ntext))
THEN
12238 IF(nums > 0 .AND. dnum(1) > 0.5) mreqenf=nint(dnum(1),mpi)
12239 IF(nums > 1 .AND. dnum(2) > 0.5) mreqena=nint(dnum(2),mpi)
12240 IF(nums > 2 .AND. dnum(3) > 0.5) iteren=nint(dnum(1)*dnum(3),mpi)
12244 keystx=
'printrecord'
12245 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12246 IF(100*mat >= 80*max(npat,ntext))
THEN
12247 IF(nums > 0) nrecpr=nint(dnum(1),mpi)
12248 IF(nums > 1) nrecp2=nint(dnum(2),mpi)
12253 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12254 IF(100*mat >= 80*max(npat,ntext))
THEN
12255 IF (nums > 0.AND.dnum(1) > 0.) mxrec=nint(dnum(1),mpi)
12260 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12261 IF(100*mat >= 80*max(npat,ntext))
THEN
12262 IF (nums > 0.AND.dnum(1) >= 0.) ncache=nint(dnum(1),mpi)
12263 IF (nums == 2.AND.dnum(2) > 0..AND.dnum(2) <= 1.0) &
12264 fcache(1)=real(dnum(2),mps)
12265 IF (nums >= 4)
THEN
12267 fcache(k)=real(dnum(k+1),mps)
12274 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12275 IF(100*mat >= 80*max(npat,ntext))
THEN
12280 chicut=real(dnum(1),mps)
12281 IF(chicut < 1.0) chicut=-1.0
12285 chirem=real(dnum(2),mps)
12286 IF(chirem < 1.0) chirem=1.0
12287 IF(chicut >= 1.0) chirem=min(chirem,chicut)
12295 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12296 IF(100*mat >= 80*max(npat,ntext))
THEN
12297 IF(nums > 0) chhuge=real(dnum(1),mps)
12298 IF(chhuge < 1.0) chhuge=1.0
12303 keystx=
'linesearch'
12304 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12305 IF(100*mat >= 80*max(npat,ntext))
THEN
12306 IF(nums > 0) lsearch=nint(dnum(1),mpi)
12311 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12312 IF(100*mat >= 80*max(npat,ntext))
THEN
12313 IF(nums > 0) lfitnp=nint(dnum(1),mpi)
12314 IF(nums > 1) lfitbb=nint(dnum(2),mpi)
12318 keystx=
'regularization'
12319 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12320 IF(100*mat >= 80*max(npat,ntext))
THEN
12322 regula=real(dnum(1),mps)
12323 IF(nums >= 2) regpre=real(dnum(2),mps)
12327 keystx=
'regularisation'
12328 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12329 IF(100*mat >= 80*max(npat,ntext))
THEN
12331 regula=real(dnum(1),mps)
12332 IF(nums >= 2) regpre=real(dnum(2),mps)
12337 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12338 IF(100*mat >= 80*max(npat,ntext))
THEN
12339 regpre=real(dnum(1),mps)
12344 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12345 IF(100*mat >= 80*max(npat,ntext))
THEN
12346 matrit=nint(dnum(1),mpi)
12351 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12352 IF(100*mat >= 80*max(npat,ntext))
THEN
12354 IF (nums > 0.AND.dnum(1) > 0.) matmon=nint(dnum(1),mpi)
12359 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12360 IF(100*mat >= 80*max(npat,ntext))
THEN
12361 IF(nums > 0) mbandw=nint(dnum(1),mpi)
12362 IF(mbandw < 0) mbandw=-1
12363 IF(nums > 1) lprecm=nint(dnum(2),mpi)
12387 keystx=
'outlierdownweighting'
12388 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12389 IF(100*mat >= 80*max(npat,ntext))
THEN
12390 lhuber=nint(dnum(1),mpi)
12391 IF(lhuber > 0.AND.lhuber <= 2) lhuber=2
12395 keystx=
'dwfractioncut'
12396 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12397 IF(100*mat >= 80*max(npat,ntext))
THEN
12398 dwcut=real(dnum(1),mps)
12399 IF(dwcut > 0.5) dwcut=0.5
12403 keystx=
'maxlocalcond'
12404 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12405 IF(100*mat >= 80*max(npat,ntext))
THEN
12406 IF (nums > 0.AND.dnum(1) > 0.0) cndlmx=real(dnum(1),mps)
12411 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12412 IF(100*mat >= 80*max(npat,ntext))
THEN
12413 prange=abs(real(dnum(1),mps))
12418 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12419 IF(100*mat >= 80*max(npat,ntext))
THEN
12425 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12426 IF(100*mat >= 80*max(npat,ntext))
THEN
12431 keystx=
'memorydebug'
12432 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12433 IF(100*mat >= 80*max(npat,ntext))
THEN
12435 IF (nums > 0.AND.dnum(1) > 0.0) memdbg=nint(dnum(1),mpi)
12439 keystx=
'globalcorr'
12440 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12441 IF(100*mat >= 80*max(npat,ntext))
THEN
12446 keystx=
'printcounts'
12447 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12448 IF(100*mat >= 80*max(npat,ntext))
THEN
12450 IF (nums > 0) ipcntr=nint(dnum(1),mpi)
12454 keystx=
'weightedcons'
12455 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12456 IF(100*mat >= 80*max(npat,ntext))
THEN
12458 IF (nums > 0) iwcons=nint(dnum(1),mpi)
12462 keystx=
'skipemptycons'
12463 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12464 IF(100*mat >= 80*max(npat,ntext))
THEN
12469 keystx=
'resolveredundancycons'
12470 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12471 IF(100*mat >= 80*max(npat,ntext))
THEN
12476 keystx=
'withelimination'
12477 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12478 IF(100*mat >= 80*max(npat,ntext))
THEN
12483 keystx=
'postprocessing'
12484 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12485 IF(100*mat >= 80*max(npat,ntext))
THEN
12486 lenpostproc=ib-
keyb-1
12487 cpostproc(1:lenpostproc)=text(
keyb+2:ib)
12492 keystx=
'withLAPACKelimination'
12493 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12494 IF(100*mat >= 80*max(npat,ntext))
THEN
12500 keystx=
'withmultipliers'
12501 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12502 IF(100*mat >= 80*max(npat,ntext))
THEN
12507 keystx=
'checkinput'
12508 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12509 IF(100*mat >= 80*max(npat,ntext))
THEN
12511 IF (nums > 0) icheck=nint(dnum(1),mpi)
12515 keystx=
'checkparametergroups'
12516 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12517 IF(100*mat >= 80*max(npat,ntext))
THEN
12522 keystx=
'monitorresiduals'
12523 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12524 IF(100*mat >= 80*max(npat,ntext))
THEN
12526 IF (nums > 0) imonit=nint(dnum(1),mpi)
12527 IF (nums > 1) measbins=max(measbins,nint(dnum(2),mpi))
12531 keystx=
'monitorpulls'
12532 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12533 IF(100*mat >= 80*max(npat,ntext))
THEN
12536 IF (nums > 0) imonit=nint(dnum(1),mpi)
12537 IF (nums > 1) measbins=max(measbins,nint(dnum(2),mpi))
12541 keystx=
'monitorprogress'
12542 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12543 IF(100*mat >= 80*max(npat,ntext))
THEN
12546 IF (nums > 0) monpg1=max(1,nint(dnum(1),mpi))
12547 IF (nums > 1) monpg2=max(1,nint(dnum(2),mpi))
12551 keystx=
'scaleerrors'
12552 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12553 IF(100*mat >= 80*max(npat,ntext))
THEN
12555 IF (nums > 0) dscerr(1:2)=dnum(1)
12556 IF (nums > 1) dscerr(2)=dnum(2)
12560 keystx=
'iterateentries'
12561 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12562 IF(100*mat >= 80*max(npat,ntext))
THEN
12563 iteren=huge(iteren)
12564 IF (nums > 0) iteren=nint(dnum(1),mpi)
12569 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12570 IF(100*mat >= 80*max(npat,ntext))
THEN
12578 WRITE(*,*)
'WARNING: multithreading not available'
12584 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12585 IF(100*mat >= 80*max(npat,ntext))
THEN
12586 WRITE(*,*)
'WARNING: keyword COMPRESS is obsolete (compression is default)'
12598 keystx=
'countrecords'
12599 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12600 IF(100*mat >= 80*max(npat,ntext))
THEN
12606 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12607 IF(100*mat >= 80*max(npat,ntext).AND.mnrsel < 100)
THEN
12608 nl=min(nums,100-mnrsel)
12610 lbmnrs(mnrsel+k)=nint(dnum(k),mpi)
12616 keystx=
'pairentries'
12617 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12618 IF(100*mat >= 80*max(npat,ntext))
THEN
12621 IF (nums > 0.AND.dnum(1) > 0.0)
THEN
12622 mreqpe=nint(dnum(1),mpi)
12623 IF (nums >= 2.AND.dnum(2) >= dnum(1)) mhispe=nint(dnum(2),mpi)
12624 IF (nums >= 3.AND.dnum(3) >= dnum(1)) msngpe=nint(dnum(3),mpi)
12630 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12631 IF(100*mat >= 80*max(npat,ntext))
THEN
12632 wolfc1=real(dnum(1),mps)
12633 wolfc2=real(dnum(2),mps)
12640 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12641 IF(100*mat >= 80*max(npat,ntext))
THEN
12643 IF (dnum(1) < 1.0e-10_mpd.OR.dnum(1) > 1.0e-04_mpd)
THEN
12644 WRITE(*,*)
'ERROR: need 1.0D-10 <= MRESTL ', &
12645 '<= 1.0D-04, but get ', dnum(1)
12654 keystx=
'mrestranscond'
12655 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12656 IF(100*mat >= 80*max(npat,ntext))
THEN
12664 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12665 IF(100*mat >= 80*max(npat,ntext))
THEN
12667 mrmode = int(dnum(1),mpi)
12672 keystx=
'nofeasiblestart'
12673 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12674 IF(100*mat >= 80*max(npat,ntext))
THEN
12680 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12681 IF(100*mat >= 80*max(npat,ntext))
THEN
12686 keystx=
'readerroraseof'
12687 mat=matint(text(ia:ib),keystx,npat,ntext)
12688 IF(100*mat >= 80*max(npat,ntext))
THEN
12694 keystx=
'LAPACKwitherrors'
12695 mat=matint(text(ia:ib),keystx,npat,ntext)
12696 IF(100*mat >= 80*max(npat,ntext))
THEN
12701 keystx=
'debugPARDISO'
12702 mat=matint(text(ia:ib),keystx,npat,ntext)
12703 IF(100*mat >= 80*max(npat,ntext))
THEN
12708 keystx=
'blocksizePARDISO'
12709 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12710 IF(100*mat >= 80*max(npat,ntext).AND.mnrsel < 100)
THEN
12711 nl=min(nums,10-mpdbsz)
12713 IF (nint(dnum(k),mpi) > 0)
THEN
12714 IF (mpdbsz == 0)
THEN
12716 ipdbsz(mpdbsz)=nint(dnum(k),mpi)
12717 ELSE IF (nint(dnum(k),mpi) > ipdbsz(mpdbsz))
THEN
12719 ipdbsz(mpdbsz)=nint(dnum(k),mpi)
12727 keystx=
'fortranfiles'
12728 mat=matint(text(ia:ib),keystx,npat,ntext)
12729 IF(mat == max(npat,ntext))
RETURN
12732 mat=matint(text(ia:ib),keystx,npat,ntext)
12733 IF(mat == max(npat,ntext))
RETURN
12735 keystx=
'closeandreopen'
12736 mat=matint(text(ia:ib),keystx,npat,ntext)
12737 IF(mat == max(npat,ntext))
RETURN
12741 IF(nums /= 0) nkey=0
12744 WRITE(*,*)
'**************************************************'
12746 WRITE(*,*)
'Unknown keyword(s): ',text(1:min(nab,50))
12748 WRITE(*,*)
'**************************************************'
1277010
IF(nkey > 0)
THEN
12774 lpvs=nint(dnum(1),mpi)
12776 CALL additem(lenparameters,listparameters,lpvs,dnum(2))
12777 CALL additem(lenpresigmas,listpresigmas,lpvs,dnum(3))
12779 WRITE(*,*)
'Line',nline,
' error, label=',lpvs
12781 ELSE IF(nums /= 0)
THEN
12783 WRITE(*,*)
'Wrong text in line',nline
12784 WRITE(*,*)
'Status: new parameter'
12785 WRITE(*,*)
'> ',text(1:nab)
12787 ELSE IF(lkey == 3)
THEN
12789 IF(nums >= 1.AND.nums <= 2)
THEN
12791 CALL additem(lenconstraints,listconstraints,lpvs,dnum(1))
12793 IF(iwcons > 0) lpvs=-2
12795 IF(nums == 2) plvs=dnum(2)
12796 CALL additem(lenconstraints,listconstraints,lpvs,plvs)
12799 WRITE(*,*)
'Wrong text in line',nline
12800 WRITE(*,*)
'Status: new keyword constraint'
12801 WRITE(*,*)
'> ',text(1:nab)
12803 ELSE IF(lkey == 4)
THEN
12805 nummeasurements=nummeasurements+1
12807 CALL additem(lenmeasurements,listmeasurements,lpvs,dnum(1))
12809 CALL additem(lenmeasurements,listmeasurements,lpvs,dnum(2))
12812 WRITE(*,*)
'Wrong text in line',nline
12813 WRITE(*,*)
'Status: new keyword measurement'
12814 WRITE(*,*)
'> ',text(1:nab)
12816 ELSE IF(lkey == 5.AND.
keyb <
keyc)
THEN
12818 IF(nums >= 1) miter=nint(dnum(1),mpi)
12819 IF(miter >= 1) mitera=miter
12820 dflim=real(dnum(2),mps)
12824 mat=matint(text(
keyb+1:
keyc),keystx,npat,ntext)
12825 IF(100*mat >= 80*max(npat,ntext))
THEN
12829 ELSE IF(i == 2)
THEN
12832 ELSE IF(i == 3)
THEN
12835 ELSE IF(i == 4)
THEN
12838 ELSE IF(i == 5)
THEN
12841 ELSE IF(i == 6)
THEN
12844 ELSE IF(i == 7)
THEN
12848 ELSE IF(i == 8)
THEN
12851 ELSE IF(i == 9)
THEN
12855 ELSE IF(i == 10)
THEN
12864 ELSE IF(nkey == 0)
THEN
12867 lpvs=nint(dnum(1),mpi)
12869 CALL additem(lenparameters,listparameters,lpvs,dnum(2))
12870 CALL additem(lenpresigmas,listpresigmas,lpvs,dnum(3))
12872 WRITE(*,*)
'Line',nline,
' error, label=',lpvs
12874 ELSE IF(nums > 1.AND.nums < 3)
THEN
12876 WRITE(*,*)
'Wrong text in line',nline
12877 WRITE(*,*)
'Status continuation parameter'
12878 WRITE(*,*)
'> ',text(1:nab)
12881 ELSE IF(lkey == 3)
THEN
12884 label=nint(dnum(i),mpi)
12885 IF(label <= 0) ier=1
12887 IF(mod(nums,2) /= 0) ier=1
12890 lpvs=nint(dnum(i),mpi)
12892 CALL additem(lenconstraints,listconstraints,lpvs,plvs)
12896 WRITE(*,*)
'Wrong text in line',nline
12897 WRITE(*,*)
'Status continuation constraint'
12898 WRITE(*,*)
'> ',text(1:nab)
12901 ELSE IF(lkey == 4)
THEN
12905 label=nint(dnum(i),mpi)
12906 IF(label <= 0) ier=1
12908 IF(mod(nums,2) /= 0) ier=1
12912 lpvs=nint(dnum(i),mpi)
12914 CALL additem(lenmeasurements,listmeasurements,lpvs,plvs)
12918 WRITE(*,*)
'Wrong text in line',nline
12919 WRITE(*,*)
'Status continuation measurement'
12920 WRITE(*,*)
'> ',text(1:nab)
12922 ELSE IF(lkey == 6)
THEN
12924 lpvs=nint(dnum(1),mpi)
12928 IF (text(j:j) ==
' ')
EXIT
12931 CALL additemc(lencomments,listcomments,lpvs,ctext)
12933 WRITE(*,*)
'Line',nline,
' error, label=',lpvs
12935 ELSE IF(nums /= 0)
THEN
12937 WRITE(*,*)
'Wrong text in line',nline
12938 WRITE(*,*)
'Status: continuation comment'
12939 WRITE(*,*)
'> ',text(1:nab)
12943 ELSE IF(lkey == 7)
THEN
12946 label=nint(dnum(i),mpi)
12947 IF(label <= 0.OR.label > 64) ier=1
12949 IF(mod(nums,2) /= 0) ier=1
12953 lpvs=nint(dnum(i),mpi)
12954 ipvs=nint(dnum(i+1),mpi)
12955 CALL additemi(lenpardiso,listpardiso,lpvs,ipvs)
12959 WRITE(*,*)
'Wrong text in line',nline
12960 WRITE(*,*)
'Status continuation measurement'
12961 WRITE(*,*)
'> ',text(1:nab)
12980 INTEGER(mpi),
INTENT(IN OUT) :: length
12981 TYPE(
listitem),
DIMENSION(:),
INTENT(IN OUT),
ALLOCATABLE :: list
12982 INTEGER(mpi),
INTENT(IN) :: label
12983 REAL(mpd),
INTENT(IN) :: value
12985 INTEGER(mpl) :: newSize
12986 INTEGER(mpl) :: oldSize
12987 TYPE(
listitem),
DIMENSION(:),
ALLOCATABLE :: tempList
12989 IF (label > 0.AND.
value == 0.)
RETURN
12990 IF (length == 0 )
THEN
12992 CALL mpalloc(list,newsize,
' list ')
12994 oldsize=
size(list,kind=
mpl)
12995 IF (length >= oldsize)
THEN
12996 newsize = oldsize + oldsize/5 + 100
12997 CALL mpalloc(templist,oldsize,
' temp. list ')
13000 CALL mpalloc(list,newsize,
' list ')
13001 list(1:oldsize)=templist(1:oldsize)
13006 list(length)%label=label
13007 list(length)%value=
value
13022 INTEGER(mpi),
INTENT(IN OUT) :: length
13023 TYPE(
listitemc),
DIMENSION(:),
INTENT(IN OUT),
ALLOCATABLE :: list
13024 INTEGER(mpi),
INTENT(IN) :: label
13025 CHARACTER(len = itemCLen),
INTENT(IN) :: text
13027 INTEGER(mpl) :: newSize
13028 INTEGER(mpl) :: oldSize
13029 TYPE(
listitemc),
DIMENSION(:),
ALLOCATABLE :: tempList
13031 IF (label > 0.AND.text ==
'')
RETURN
13032 IF (length == 0 )
THEN
13034 CALL mpalloc(list,newsize,
' list ')
13036 oldsize=
size(list,kind=
mpl)
13037 IF (length >= oldsize)
THEN
13038 newsize = oldsize + oldsize/5 + 100
13039 CALL mpalloc(templist,oldsize,
' temp. list ')
13042 CALL mpalloc(list,newsize,
' list ')
13043 list(1:oldsize)=templist(1:oldsize)
13048 list(length)%label=label
13049 list(length)%text=text
13064 INTEGER(mpi),
INTENT(IN OUT) :: length
13065 TYPE(
listitemi),
DIMENSION(:),
INTENT(IN OUT),
ALLOCATABLE :: list
13066 INTEGER(mpi),
INTENT(IN) :: label
13067 INTEGER(mpi),
INTENT(IN) :: ivalue
13069 INTEGER(mpl) :: newSize
13070 INTEGER(mpl) :: oldSize
13071 TYPE(
listitemi),
DIMENSION(:),
ALLOCATABLE :: tempList
13073 IF (length == 0 )
THEN
13075 CALL mpalloc(list,newsize,
' list ')
13077 oldsize=
size(list,kind=
mpl)
13078 IF (length >= oldsize)
THEN
13079 newsize = oldsize + oldsize/5 + 100
13080 CALL mpalloc(templist,oldsize,
' temp. list ')
13083 CALL mpalloc(list,newsize,
' list ')
13084 list(1:oldsize)=templist(1:oldsize)
13089 list(length)%label=label
13090 list(length)%ivalue=ivalue
13104 CHARACTER (LEN=*),
INTENT(IN) :: text
13105 CHARACTER (LEN=16) :: textc
13114 textl(ka:kb)=text(1:l)
13118 textc=text(1:l)//
'-end'
13126 textl(ka:kb)=textc(1:l)
13153 INTEGER(mpi),
INTENT(IN) :: lun
13154 CHARACTER (LEN=*),
INTENT(IN) :: fname
13155 CHARACTER (LEN=33) :: nafile
13156 CHARACTER (LEN=33) :: nbfile
13162 CALL peend(17,
'Aborted, file name too long')
13163 stop
'File name too long '
13166 nafile(l+1:l+1)=
'~'
13168 INQUIRE(file=nafile(1:l),exist=ex)
13170 INQUIRE(file=nafile(1:l+1),exist=ex)
13172 CALL system(
'rm '//nafile)
13175 nafile(l+1:l+1)=
' '
13176 CALL system(
'mv '//nafile//nbfile)
13178 OPEN(unit=lun,file=fname)
13189 REAL,
DIMENSION(2) :: ta
13209 IF(ncount > 1)
THEN
13211 nsecd1=int(delta,
mpi)
13213 minut1=nsecd1/60-60*nhour1
13214 secnd1=delta-60*(minut1+60*nhour1)
13216 nsecd2=int(delta,
mpi)
13218 minut2=nsecd2/60-60*nhour2
13219 secnd2=delta-60*(minut2+60*nhour2)
13220 WRITE(*,101) nhour1,minut1,secnd1, nhour2,minut2,secnd2
13225101
FORMAT(i4,
' h',i3,
' min',f5.1,
' sec total',18x,
'elapsed', &
13226 i4,
' h',i3,
' min',f5.1,
' sec')
13240 INTEGER(mpi),
INTENT(IN) :: icode
13241 CHARACTER (LEN=*),
INTENT(IN) :: cmessage
13243 CALL mvopen(9,
'millepede.end')
13244 WRITE(9,101) icode, cmessage
13245101
FORMAT(1x,i4,3x,a)
13249END SUBROUTINE peend
13261 INTEGER(mpi),
INTENT(IN) :: kfile
13262 INTEGER(mpi),
INTENT(IN) :: ithr
13263 INTEGER(mpi),
INTENT(OUT) :: ierr
13265 INTEGER(mpi),
DIMENSION(13) :: ibuff
13266 INTEGER(mpi) :: ioff
13267 INTEGER(mpi) :: ios
13269 INTEGER(mpi) :: lfn
13270 INTEGER(mpi) :: lun
13271 INTEGER(mpi) :: moddate
13272 CHARACTER (LEN=1024) :: fname
13273 CHARACTER (LEN=7) :: cfile
13278 SUBROUTINE openc(filename, lfn, lun, ios)
BIND(c)
13280 CHARACTER(kind=c_char),
DIMENSION(*),
INTENT(IN) :: filename
13281 INTEGER(c_int),
INTENT(IN),
VALUE :: lfn
13282 INTEGER(c_int),
INTENT(IN),
VALUE :: lun
13283 INTEGER(c_int),
INTENT(INOUT) :: ios
13284 END SUBROUTINE openc
13296 fname(k:k)=
tfd(ioff+k)
13301 IF(kfile <=
nfilf)
THEN
13304 OPEN(lun,file=fname(1:lfn),iostat=ios, form=
'UNFORMATTED')
13305 print *,
' lun ', lun, ios
13309 CALL openc(fname(1:lfn),lfn,lun,ios)
13311 WRITE(*,*)
'Opening of C-files not supported.'
13318 WRITE(*,*)
'Open error for file ',fname(1:lfn), ios
13319 IF (moddate /= 0)
THEN
13320 WRITE(cfile,
'(I7)') kfile
13321 CALL peend(15,
'Aborted, open error(s) for binary file ' // cfile)
13322 stop
'PEREAD: open error'
13327 ios=stat(fname(1:lfn),ibuff)
13331 WRITE(*,*)
'STAT error for file ',fname(1:lfn), ios
13335 IF (moddate /= 0)
THEN
13336 IF (ibuff(10) /= moddate)
THEN
13337 WRITE(cfile,
'(I7)') kfile
13338 CALL peend(19,
'Aborted, binary file modified (date) ' // cfile)
13339 stop
'PEREAD: file modified'
13342 yfd(kfile)=ibuff(10)
13357 INTEGER(mpi),
INTENT(IN) :: kfile
13358 INTEGER(mpi),
INTENT(IN) :: ithr
13360 INTEGER(mpi) :: lun
13364 SUBROUTINE closec(lun)
BIND(c)
13366 INTEGER(c_int),
INTENT(IN),
VALUE :: lun
13373 IF(kfile <=
nfilf)
THEN
13392 INTEGER(mpi),
INTENT(IN) :: kfile
13394 INTEGER(mpi) :: lun
13398 SUBROUTINE resetc(lun)
BIND(c)
13400 INTEGER(c_int),
INTENT(IN),
VALUE :: lun
13406 IF (kfile <=
nfilf)
THEN
13425 INTEGER(mpi) :: ipgrp
13426 INTEGER(mpi) :: irank
13427 INTEGER(mpi) :: isize
13428 INTEGER(mpi) :: ivoff
13429 INTEGER(mpi) :: itgbi
13431 INTEGER(mpi) :: msize
13432 INTEGER(mpi),
PARAMETER :: mxsize = 1000
13434 INTEGER(mpl):: length
13436 REAL(mpd),
DIMENSION(:),
ALLOCATABLE :: auxVectorD
13437 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: auxVectorI
13438 REAL(mpd),
DIMENSION(:),
ALLOCATABLE :: resParGroup
13439 REAL(mpd),
DIMENSION(:),
ALLOCATABLE :: blockParGroup
13447 IF (isize <= mxsize)
THEN
13448 msize=max(msize,isize)
13450 print *,
' CKPGRP: par. group', ipgrp,
' not checked -- too large: ', isize
13453 IF (msize == 0)
RETURN
13456 length=int(msize,mpl)*(int(msize,mpl)+1)/2
13457 CALL mpalloc(blockpargroup,length,
'(matrix) block for parameter groups (D)')
13459 CALL mpalloc(respargroup,length,
'residuals for parameter groups (D)')
13460 CALL mpalloc(auxvectori,length,
'auxiliary array (I)')
13461 CALL mpalloc(auxvectord,length,
'auxiliary array (D)')
13465 print *,
' CKPGRP par. group first label size rank'
13468 IF (isize > mxsize) cycle
13475 blockpargroup(ij)=matij(ivoff+i,ivoff+j)
13479 CALL sqminv(blockpargroup,respargroup,isize,irank, auxvectord, auxvectori)
13482 IF (isize == irank)
THEN
13486 print *,
' CKPGRP ', ipgrp,
globalparlabelindex(1,itgbi), isize, irank,
' rank deficit !!!'
13504 INTEGER(mpl) :: nan
13505 INTEGER(mpl) :: neg
13507 print *,
' Checking global matrix(D) for NANs ',
size(
globalmatd,kind=mpl)
13512 print *,
' i, nan ', i, nan
13518 print *,
' Checking diagonal elements ',
nagb
13523 print *,
' i, neg ', i, neg
13527 print *,
' CHKMAT summary ', nan, neg
13549 REAL(mpd),
INTENT(IN) :: chi2
13550 INTEGER(mpi),
INTENT(IN) :: ithrd
13551 INTEGER(mpi),
INTENT(IN) :: ndf
13552 REAL(mpd),
INTENT(IN) :: dw
13554 INTEGER(mpl) ::nadd
13582 REAL(mpd),
INTENT(OUT) ::chi2
13583 INTEGER(mpl),
INTENT(OUT) ::ndf
13584 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
character(len=1024) cpostproc
post processing string
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(mpi) lenpostproc
length of post processing string
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