912#ifdef SCOREP_USER_ENABLE
913#include "scorep/SCOREP_User.inc"
938 REAL,
DIMENSION(2) :: ta
941 INTEGER(mpi) :: iopnmp
952 INTEGER(mpi) :: nsecnd
953 INTEGER(mpi) :: ntsec
955 CHARACTER (LEN=24) :: chdate
956 CHARACTER (LEN=24) :: chost
958 CHARACTER (LEN=6) :: c6
959 INTEGER major, minor, patch
983 CALL getenv(
'HOSTNAME',chost)
984 IF (chost(1:1) ==
' ')
CALL getenv(
'HOST',chost)
985 WRITE(*,*)
'($Id: dd0c569a1aafb6f6eb2f26b9b9537a685639ca25 $)'
990 CALL ilaver( major,minor, patch )
991 WRITE(*,110) lapack64, major,minor, patch
992110
FORMAT(
' using LAPACK64 with ',(a),
', version ',i0,
'.',i0,
'.',i0)
994 WRITE(*,*)
'using Intel oneMKL PARDISO'
998 WRITE(*,111) __gnuc__ , __gnuc_minor__ , __gnuc_patchlevel__
999111
FORMAT(
' compiled with gcc ',i0,
'.',i0,
'.',i0)
1002 WRITE(*,111) __pgic__ , __pgic_minor__ , __pgic_patchlevel__
1003111
FORMAT(
' compiled with pgi ',i0,
'.',i0,
'.',i0)
1005#ifdef SCOREP_USER_ENABLE
1006 WRITE(*,*)
'instrumenting Score-P user regions'
1009 WRITE(*,*)
' < Millepede II-P starting ... ',chdate
1010 WRITE(*,*)
' ',chost
1013 WRITE(8,*)
'($Id: dd0c569a1aafb6f6eb2f26b9b9537a685639ca25 $)'
1015 WRITE(8,*)
'Log-file Millepede II-P ', chdate
1016 WRITE(8,*)
' ', chost
1018 CALL peend(-1,
'Still running or crashed')
1027 WRITE(*,*)
'!!! Checking input only, no calculation of a solution !!!'
1028 WRITE(8,*)
'!!! Checking input only, no calculation of a solution !!!'
1047 CALL getenv(
'OMP_NUM_THREADS',c6)
1049 CALL getenv(lapack64//
'_NUM_THREADS',c6)
1051 IF (c6(1:1) ==
' ')
THEN
1053 WRITE(*,*)
'Number of threads for LAPACK: unkown (empty OMP_NUM_THREADS)'
1055 WRITE(*,*)
'Number of threads for LAPACK: unkown (empty ',lapack64//
'_NUM_THREADS)'
1058 WRITE(*,*)
'Number of threads for LAPACK: ', c6
1078 CALL mvopen(lun,
'millepede.his')
1084 CALL mvopen(1,
'mpdebug.txt')
1088 times(0)=rstext-rstp
1094 times(1)=rloop1-rstext
1098 WRITE(8,*)
'Chi square cut equiv 3 st.dev applied ...'
1099 WRITE(8,*)
' in first iteration with factor',
chicut
1100 WRITE(8,*)
' in second iteration with factor',
chirem
1101 WRITE(8,*)
' (reduced by sqrt in next iterations)'
1105 WRITE(8,*)
'Down-weighting of outliers in',
lhuber,
' iterations'
1106 WRITE(8,*)
'Cut on downweight fraction',
dwcut
1110 times(2)=rloop2-rloop1
1115 CALL peend(5,
'Ended without solution (empty constraints)')
1117 CALL peend(0,
'Ended normally')
1154 CALL gmpdef(6,1,
'log10(#records) vs file number')
1155 CALL gmpdef(7,1,
'final rejection fraction vs file number')
1157 'final <Chi^2/Ndf> from accepted local fits vs file number')
1158 CALL gmpdef(9,1,
'<Ndf> from accepted local fits vs file number')
1164 rej=real(nrc-
jfd(kfl),mps)/real(nrc,mps)
1165 CALL gmpxy(6,real(kfl,mps),log10(real(nrc,mps)))
1166 CALL gmpxy(7,real(kfl,mps),rej)
1168 IF (
jfd(kfl) > 0)
THEN
1169 c2ndf=
cfd(kfl)/real(
jfd(kfl),mps)
1170 CALL gmpxy(8,real(kfl,mps),c2ndf)
1171 andf=real(
dfd(kfl),mps)/real(
jfd(kfl),mps)
1172 CALL gmpxy(9,real(kfl,mps),andf)
1189 WRITE(*,*)
'Misalignment test wire chamber'
1192 CALL hmpdef( 9,-0.0015,+0.0015,
'True - fitted displacement')
1193 CALL hmpdef(10,-0.0015,+0.0015,
'True - fitted Vdrift')
1199 sums(1)=sums(1)+diff
1200 sums(2)=sums(2)+diff*diff
1202 sums(3)=sums(3)+diff
1203 sums(4)=sums(4)+diff*diff
1205 sums(1)=0.01_mpd*sums(1)
1206 sums(2)=sqrt(0.01_mpd*sums(2))
1207 sums(3)=0.01_mpd*sums(3)
1208 sums(4)=sqrt(0.01_mpd*sums(4))
1209 WRITE(*,143)
'Parameters 1 - 100: mean =',sums(1),
'rms =',sums(2)
1210 WRITE(*,143)
'Parameters 101 - 200: mean =',sums(3),
'rms =',sums(4)
1211143
FORMAT(6x,a28,f9.6,3x,a5,f9.6)
1214 WRITE(*,*)
' I label simulated fitted diff'
1215 WRITE(*,*)
' -------------------------------------------- '
1235 WRITE(*,*)
'Misalignment test Si tracker'
1238 CALL hmpdef( 9,-0.0025,+0.0025,
'True - fitted displacement X')
1239 CALL hmpdef(10,-0.025,+0.025,
'True - fitted displacement Y')
1250 sums(1)=sums(1)+1.0_mpd
1251 sums(2)=sums(2)+diff
1252 sums(3)=sums(3)+diff*diff
1257 err=sqrt(abs(gmati))
1259 sums(7)=sums(7)+1.0_mpd
1260 sums(8)=sums(8)+diff
1261 sums(9)=sums(9)+diff*diff
1264 IF (mod(i,3) == 1)
THEN
1268 sums(4)=sums(4)+1.0_mpd
1269 sums(5)=sums(5)+diff
1270 sums(6)=sums(6)+diff*diff
1275 err=sqrt(abs(gmati))
1277 sums(7)=sums(7)+1.0_mpd
1278 sums(8)=sums(8)+diff
1279 sums(9)=sums(9)+diff*diff
1284 sums(2)=sums(2)/sums(1)
1285 sums(3)=sqrt(sums(3)/sums(1))
1286 sums(5)=sums(5)/sums(4)
1287 sums(6)=sqrt(sums(6)/sums(4))
1288 WRITE(*,143)
'Parameters 1 - 500: mean =',sums(2),
'rms =',sums(3)
1289 WRITE(*,143)
'Parameters 501 - 700: mean =',sums(5),
'rms =',sums(6)
1290 IF (sums(7) > 0.5_mpd)
THEN
1291 sums(8)=sums(8)/sums(7)
1292 sums(9)=sqrt(sums(9)/sums(7))
1293 WRITE(*,143)
'Parameter pulls, all: mean =',sums(8),
'rms =',sums(9)
1297 WRITE(*,*)
' I label simulated fitted diff'
1298 WRITE(*,*)
' -------------------------------------------- '
1310 IF (mod(i,3) == 1)
THEN
1330 WRITE(8,*)
'Record',
nrec1,
' has largest residual:',
value1
1333 WRITE(8,*)
'Record',
nrec2,
' has largest Chi^2/Ndf:',
value2
1337 WRITE(8,*)
'Record',
nrec3,
' is first with error (rank deficit/NaN)'
1341 WRITE(8,*)
'In total 3 +',
nloopn,
' loops through the data files'
1343 WRITE(8,*)
'In total 2 +',
nloopn,
' loops through the data files'
1347 WRITE(8,*)
'In total ',
mnrsit,
' internal MINRES iterations'
1355 ntsec=nint(deltat,mpi)
1356 CALL sechms(deltat,nhour,minut,secnd)
1357 nsecnd=nint(secnd,mpi)
1358 WRITE(8,*)
'Total time =',ntsec,
' seconds =',nhour,
' h',minut, &
1359 ' m',nsecnd,
' seconds'
1361 WRITE(8,*)
'end ', chdate
1368 IF(sum(
nrejec) /= 0)
THEN
1370 WRITE(8,*)
'Data records rejected in last iteration: '
1377 WRITE(*,*)
' < Millepede II-P ending ... ', chdate
1384106
FORMAT(
' PARDISO dyn. memory allocation: ',f11.6,
' GB')
1394 WRITE(*,*)
'Postprocessing:'
1404102
FORMAT(2x,i4,i10,2x,3f10.5)
1405103
FORMAT(
' Times [in sec] for text processing',f12.3/ &
1408 ' func. value ',f12.3,
' *',f4.0/ &
1409 ' func. value, global matrix, solution',f12.3,
' *',f4.0/ &
1410 ' new solution',f12.3,
' *',f4.0/)
1411105
FORMAT(
' Peak dynamic memory allocation: ',f11.6,
' GB')
1431 INTEGER(mpi) :: istop
1432 INTEGER(mpi) :: itgbi
1433 INTEGER(mpi) :: itgbl
1435 INTEGER(mpi) :: itnlim
1436 INTEGER(mpi) :: nout
1438 INTEGER(mpi),
INTENT(IN) :: ivgbi
1475 globalcorrections, itnlim, nout, rtol, istop, itn, anorm, acond, rnorm, arnorm, ynorm)
1479 globalcorrections, itnlim, nout, rtol, istop, itn, anorm, acond, rnorm, arnorm, ynorm)
1482 globalcorrections, itnlim, nout, rtol, istop, itn, anorm, acond, rnorm, arnorm, ynorm)
1488 err=sqrt(abs(real(gmati,mps)))
1489 IF(gmati < 0.0_mpd) err=-err
1490 diag=matij(ivgbi,ivgbi)
1491 gcor2=real(1.0_mpd-1.0_mpd/(gmati*diag),mps)
1493101
FORMAT(1x,
' label parameter presigma differ', &
1494 ' Error gcor^2 iit'/ 1x,
'---------',2x,5(
'-----------'),2x,
'----')
1495102
FORMAT(i10,2x,4g12.4,f7.4,i6,i4)
1515 INTEGER(mpi) :: istop
1516 INTEGER(mpi) :: itgbi
1517 INTEGER(mpi) :: itgbl
1519 INTEGER(mpi) :: itnlim
1520 INTEGER(mpi) :: nout
1522 INTEGER(mpi),
INTENT(IN) :: ivgbi
1551 mxxnrm = real(
nagb,mpd)/sqrt(epsilon(mxxnrm))
1553 trcond = 1.0_mpd/epsilon(trcond)
1554 ELSE IF(
mrmode == 2)
THEN
1562 itnlim=itnlim, rtol=rtol, maxxnorm=mxxnrm, trancond=trcond, &
1566 itnlim=itnlim, rtol=rtol, maxxnorm=mxxnrm, trancond=trcond, &
1570 itnlim=itnlim, rtol=rtol, maxxnorm=mxxnrm, trancond=trcond, &
1577 err=sqrt(abs(real(gmati,mps)))
1578 IF(gmati < 0.0_mpd) err=-err
1579 diag=matij(ivgbi,ivgbi)
1580 gcor2=real(1.0_mpd-1.0_mpd/(gmati*diag),mps)
1582101
FORMAT(1x,
' label parameter presigma differ', &
1583 ' Error gcor^2 iit'/ 1x,
'---------',2x,5(
'-----------'),2x,
'----')
1584102
FORMAT(i10,2x,4g12.4,f7.4,i6,i4)
1597 INTEGER(mpi) :: icgb
1598 INTEGER(mpi) :: irhs
1599 INTEGER(mpi) :: itgbi
1600 INTEGER(mpi) :: ivgb
1602 INTEGER(mpi) :: jcgb
1604 INTEGER(mpi) :: label
1606 INTEGER(mpi) :: inone
1609 REAL(mpd) :: drhs(4)
1610 INTEGER(mpi) :: idrh (4)
1635 IF(abs(rhs) > climit)
THEN
1641 WRITE(*,101) (idrh(l),drhs(l),l=1,irhs)
1650 WRITE(*,101) (idrh(l),drhs(l),l=1,irhs)
1653 WRITE(*,102)
' Constraints: only equation values >', climit,
' are printed'
1654101
FORMAT(
' ',4(i6,g11.3))
1655102
FORMAT(a,g11.2,a)
1668 INTEGER(mpi) :: icgb
1669 INTEGER(mpi) :: icgrp
1670 INTEGER(mpi) :: ioff
1671 INTEGER(mpi) :: itgbi
1673 INTEGER(mpi) :: jcgb
1674 INTEGER(mpi) :: label
1675 INTEGER(mpi) :: labelf
1676 INTEGER(mpi) :: labell
1677 INTEGER(mpi) :: last
1678 INTEGER(mpi) :: line1
1679 INTEGER(mpi) :: ncon
1680 INTEGER(mpi) :: ndiff
1681 INTEGER(mpi) :: npar
1682 INTEGER(mpi) :: inone
1683 INTEGER(mpi) :: itype
1684 INTEGER(mpi) :: ncgbd
1685 INTEGER(mpi) :: ncgbr
1686 INTEGER(mpi) :: ncgbw
1687 INTEGER(mpi) :: ncgrpd
1688 INTEGER(mpi) :: ncgrpr
1689 INTEGER(mpi) :: next
1691 INTEGER(mpl):: length
1692 INTEGER(mpl) :: rows
1694 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: vecParConsOffsets
1695 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: vecParConsList
1696 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: vecConsParOffsets
1697 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: vecConsParList
1698 INTEGER(mpi),
DIMENSION(:,:),
ALLOCATABLE :: matConsGroupIndex
1711 IF(last < 0.AND.label < 0)
THEN
1714 IF(itype == 2) ncgbw=ncgbw+1
1721 IF(label > 0.AND.itype == 2)
THEN
1728 IF (ncgbw == 0)
THEN
1729 WRITE(*,*)
'GRPCON:',
ncgb,
' constraints found in steering files'
1731 WRITE(*,*)
'GRPCON:',
ncgb,
' constraints found in steering files,',ncgbw,
'weighted'
1736 length=
ncgb+1; rows=3
1749 CALL mpalloc(matconsgroupindex,rows,length,
'group index for constraint (I)')
1752 CALL mpalloc(vecconsparoffsets,length,
'offsets for global par list for cons. (I)')
1754 CALL mpalloc(vecparconsoffsets,length,
'offsets for cons. list for global par. (I)')
1755 vecparconsoffsets(1)=0
1757 vecparconsoffsets(i+1)=vecparconsoffsets(i)+
globalparcons(i)
1761 length=vecparconsoffsets(
ntgb+1)
1762 CALL mpalloc(vecconsparlist,length,
'global par. list for constraint (I)')
1763 CALL mpalloc(vecparconslist,length,
'constraint list for global par. (I)')
1768 vecconsparoffsets(1)=ioff
1780 vecparconslist(vecparconsoffsets(itgbi)+
globalparcons(itgbi))=icgb
1782 vecconsparlist(ioff+npar)=itgbi
1788 CALL sort1k(vecconsparlist(ioff+1),npar)
1792 next=vecconsparlist(ioff+j)
1793 IF (next /= last)
THEN
1795 vecconsparlist(ioff+ndiff) = next
1804 vecconsparoffsets(icgb+1)=ioff
1817 print *,
' Constraint #parameters first par. last par. first line'
1825 npar=vecconsparoffsets(icgb+1)-vecconsparoffsets(icgb)
1829 print *, jcgb, npar, labelf, labell, line1
1832 icgrp=matconsgroupindex(1,icgb)
1833 IF (icgrp == 0)
THEN
1835 DO i=vecconsparoffsets(icgb)+1, vecconsparoffsets(icgb+1)
1836 itgbi=vecconsparlist(i)
1838 DO j=vecparconsoffsets(itgbi)+1,vecparconsoffsets(itgbi+1)
1839 icgrp=matconsgroupindex(1,vecparconslist(j))
1845 IF (icgrp == 0)
THEN
1852 matconsgroupindex(2,icgb)=jcgb
1853 matconsgroupindex(3,icgb)=icgb
1854 DO i=vecconsparoffsets(icgb)+1, vecconsparoffsets(icgb+1)
1855 itgbi=vecconsparlist(i)
1858 DO j=vecparconsoffsets(itgbi)+1,vecparconsoffsets(itgbi+1)
1859 matconsgroupindex(1,vecparconslist(j))=icgrp
1863 WRITE(*,*)
'GRPCON:',
ncgrp,
' disjoint constraints groups built'
1871 icgb=matconsgroupindex(3,jcgb)
1876 icgrp=matconsgroupindex(1,jcgb)
1896 print *,
' cons.group first con. first par. last par. #cons #par'
1907 print *, icgrp,
matconsgroups(1,icgrp), labelf, labell, ncon, npar
1910 IF (ncon == npar)
THEN
1917 print *, icgrp,
matconsgroups(1,icgrp), labelf, labell,
' : cons.group resolved'
1932 print *, icgrp,
matconsgroups(1,icgrp), labelf, labell,
' : cons.group redundant'
1937 IF (ncgrpr > 0)
THEN
1938 WRITE(*,*)
'GRPCON:',ncgbr,
' redundancy constraints in ', ncgrpr,
' groups resolved'
1942 IF (ncgrpd > 0)
THEN
1943 WRITE(*,*)
'GRPCON:',ncgbd,
' redundancy constraints in ', ncgrpd,
' groups detected'
1966 INTEGER(mpi) :: icgb
1967 INTEGER(mpi) :: icgrp
1968 INTEGER(mpi) :: ifrst
1969 INTEGER(mpi) :: ilast
1970 INTEGER(mpi) :: isblck
1971 INTEGER(mpi) :: itgbi
1972 INTEGER(mpi) :: ivgb
1974 INTEGER(mpi) :: jcgb
1975 INTEGER(mpi) :: jfrst
1976 INTEGER(mpi) :: label
1977 INTEGER(mpi) :: labelf
1978 INTEGER(mpi) :: labell
1979 INTEGER(mpi) :: ncon
1980 INTEGER(mpi) :: ngrp
1981 INTEGER(mpi) :: npar
1982 INTEGER(mpi) :: ncnmxb
1983 INTEGER(mpi) :: ncnmxg
1984 INTEGER(mpi) :: nprmxb
1985 INTEGER(mpi) :: nprmxg
1986 INTEGER(mpi) :: inone
1987 INTEGER(mpi) :: nvar
1989 INTEGER(mpl):: length
1990 INTEGER(mpl) :: rows
1992 INTEGER(mpi),
DIMENSION(:,:),
ALLOCATABLE :: matConsGroupIndex
2005 length=
ncgrp+1; rows=3
2010 CALL mpalloc(matconsgroupindex,rows,length,
'group index for constraint (I)')
2047 IF (nvar > 0 .OR.
iskpec == 0)
THEN
2050 matconsgroupindex(1,
ncgb)=ngrp+1
2051 matconsgroupindex(2,
ncgb)=icgb
2052 matconsgroupindex(3,
ncgb)=nvar
2055 IF (
ncgb > ncon) ngrp=ngrp+1
2061 WRITE(*,*)
'PRPCON:',
ncgbe,
' empty constraints skipped'
2063 WRITE(*,*)
'PRPCON:',
ncgbe,
' empty constraints detected, to be fixed !!!'
2064 WRITE(*,*)
' (use option "skipemptycons" to skip those)'
2069 WRITE(*,*)
'!!! Switch to "-C" (checking input only), no calculation of a solution !!!'
2070 WRITE(8,*)
'!!! Switch to "-C" (checking input only), no calculation of a solution !!!'
2075 WRITE(*,*)
'PRPCON:',
ncgb,
' constraints accepted'
2078 IF(
ncgb == 0)
RETURN
2085 icgb=matconsgroupindex(2,jcgb)
2090 icgrp=matconsgroupindex(1,jcgb)
2116 WRITE(*,*)
' Cons. sorted index #var.par. first line first label last label'
2117 WRITE(*,*)
' Cons. group index first cons. last cons. first label last label'
2118 WRITE(*,*)
' Cons. block index first group last group first label last label'
2124 nvar=matconsgroupindex(3,jcgb)
2128 WRITE(*,*)
' Cons. sorted', jcgb, nvar, &
2131 WRITE(*,*)
' Cons. sorted', jcgb,
' empty (0)', &
2140 WRITE(*,*)
' Cons. group ', icgrp,
matconsgroups(1,icgrp), &
2153 DO i=icgrp,isblck,-1
2167 WRITE(*,*)
' Cons. block ',
ncblck, isblck, icgrp, labelf, labell
2185 ncnmxb=max(ncnmxb,ncon)
2186 nprmxb=max(nprmxb,npar)
2211 ncnmxg=max(ncnmxg,ncon)
2212 nprmxg=max(nprmxg,npar)
2247 WRITE(*,*)
'PRPCON: constraints split into ',
ncgrp,
'(disjoint) groups,'
2248 WRITE(*,*)
' groups combined into ',
ncblck,
'(non overlapping) blocks'
2249 WRITE(*,*)
' max group size (cons., par.) ', ncnmxg, nprmxg
2250 WRITE(*,*)
' max block size (cons., par.) ', ncnmxb, nprmxb
2268 INTEGER(mpi) :: icgb
2269 INTEGER(mpi) :: icgrp
2271 INTEGER(mpi) :: ifirst
2272 INTEGER(mpi) :: ilast
2273 INTEGER(mpl) :: ioffc
2274 INTEGER(mpl) :: ioffp
2275 INTEGER(mpi) :: irank
2276 INTEGER(mpi) :: ipar0
2277 INTEGER(mpi) :: itgbi
2278 INTEGER(mpi) :: ivgb
2280 INTEGER(mpi) :: jcgb
2282 INTEGER(mpi) :: label
2283 INTEGER(mpi) :: ncon
2284 INTEGER(mpi) :: npar
2285 INTEGER(mpi) :: nrank
2286 INTEGER(mpi) :: inone
2291 INTEGER(mpl):: length
2292 REAL(mpd),
DIMENSION(:),
ALLOCATABLE :: matConstraintsT
2293 REAL(mpd),
DIMENSION(:),
ALLOCATABLE :: auxVectorD
2294 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: auxVectorI
2298 IF(
ncgb == 0)
RETURN
2307 CALL mpalloc(auxvectori,length,
'auxiliary array (I)')
2308 CALL mpalloc(auxvectord,length,
'auxiliary array (D)')
2311 CALL mpalloc(matconstraintst,length,
'transposed matrix of constraints (blocks)')
2312 matconstraintst=0.0_mpd
2325 WRITE(*,*)
' Constraint group, #con, rank', icgrp, ncon, 0,
' (empty)'
2328 DO jcgb=ifirst,ilast
2340 IF(ivgb > 0) matconstraintst(int(jcgb-ifirst,mpl)*int(npar,mpl)+ivgb-ipar0+ioffc)= &
2341 matconstraintst(int(jcgb-ifirst,mpl)*int(npar,mpl)+ivgb-ipar0+ioffc)+factr
2348 DO ll=ioffc+1,ioffc+npar
2354 matconstraintst(int(i-1,mpl)*int(npar,mpl)+ll)* &
2355 matconstraintst(int(j-1,mpl)*int(npar,mpl)+ll)
2361 IF (
icheck > 1 .OR. irank < ncon)
THEN
2362 WRITE(*,*)
' Constraint group, #con, rank', icgrp, ncon, irank
2363 IF (irank < ncon)
THEN
2364 WRITE(*,*)
' .. rank deficit !! '
2365 WRITE(*,*)
' E.g. fix all parameters and remove all constraints related to label ', &
2370 ioffc=ioffc+int(npar,mpl)*int(ncon,mpl)
2377 WRITE(*,*)
'Rank of product matrix of constraints is',nrank, &
2378 ' for',
ncgb,
' constraint equations'
2379 WRITE(8,*)
'Rank of product matrix of constraints is',nrank, &
2380 ' for',
ncgb,
' constraint equations'
2381 IF(nrank <
ncgb)
THEN
2382 WRITE(*,*)
'Warning: insufficient constraint equations!'
2383 WRITE(8,*)
'Warning: insufficient constraint equations!'
2386 WRITE(*,*)
' --> enforcing SUBITO mode'
2387 WRITE(8,*)
' --> enforcing SUBITO mode'
2394 print *,
'QL decomposition of constraints matrix'
2397 WRITE(
lunlog,*)
'QL decomposition of constraints matrix'
2409 CALL lpqldec(matconstraintst,evmin,evmax)
2413 print *,
' largest |eigenvalue| of L: ', evmax
2414 print *,
' smallest |eigenvalue| of L: ', evmin
2415 IF (evmin == 0.0_mpd.AND.
icheck == 0)
THEN
2416 CALL peend(27,
'Aborted, singular QL decomposition of constraints matrix')
2417 stop
'FEASMA: stopping due to singular QL decomposition of constraints matrix'
2443 INTEGER(mpi) :: icgb
2444 INTEGER(mpi) :: icgrp
2445 INTEGER(mpi) :: iter
2446 INTEGER(mpi) :: itgbi
2447 INTEGER(mpi) :: ivgb
2448 INTEGER(mpi) :: ieblck
2449 INTEGER(mpi) :: isblck
2450 INTEGER(mpi) :: ifirst
2451 INTEGER(mpi) :: ilast
2453 INTEGER(mpi) :: jcgb
2454 INTEGER(mpi) :: label
2455 INTEGER(mpi) :: inone
2456 INTEGER(mpi) :: ncon
2458 REAL(mps),
INTENT(IN) :: concut
2459 INTEGER(mpi),
INTENT(OUT) :: iact
2466 REAL(mpd),
DIMENSION(:),
ALLOCATABLE :: vecCorrections
2470 IF(
ncgb == 0)
RETURN
2500 sum1=sqrt(sum1/real(
ncgb,mpd))
2501 sum2=sum2/real(
ncgb,mpd)
2503 IF(iter == 1.AND.sum1 < concut)
RETURN
2505 IF(iter == 1.AND.
ncgb <= 12)
THEN
2507 WRITE(*,*)
'Constraint equation discrepancies:'
2509101
FORMAT(4x,4(i5,g12.4))
2511103
FORMAT(10x,
' Cut on rms value is',g8.1)
2516 WRITE(*,*)
'Improve constraints'
2520 WRITE(*,102) iter,sum1,sum2,sum3
2521102
FORMAT(i6,
' rms',g12.4,
' avrg_abs',g12.4,
' max_abs',g12.4)
2523 CALL mpalloc(veccorrections,int(
nvgb,mpl),
'constraint corrections')
2524 veccorrections=0.0_mpd
2532 ieblck=isblck+(ncon*(ncon+1))/2
2601 INTEGER(mpi) :: iact
2602 INTEGER(mpi) :: ierrc
2603 INTEGER(mpi) :: ierrf
2604 INTEGER(mpi) :: ioffp
2606 INTEGER(mpi) :: ithr
2607 INTEGER(mpi) :: jfile
2608 INTEGER(mpi) :: jrec
2610 INTEGER(mpi) :: kfile
2613 INTEGER(mpi) :: mpri
2615 INTEGER(mpi) :: nact
2616 INTEGER(mpi) :: nbuf
2617 INTEGER(mpi) :: ndata
2618 INTEGER(mpi) :: noff
2619 INTEGER(mpi) :: noffs
2620 INTEGER(mpi) :: npointer
2621 INTEGER(mpi) :: npri
2625 INTEGER(mpi) :: nrpr
2626 INTEGER(mpi) :: nthr
2627 INTEGER(mpi) :: ntot
2628 INTEGER(mpi) :: maxRecordSize
2629 INTEGER(mpi) :: maxRecordFile
2631 INTEGER(mpi),
INTENT(OUT) :: more
2641 CHARACTER (LEN=7) :: cfile
2646 SUBROUTINE readc(bufferD, bufferF, bufferI, bufferLength, lun, err)
BIND(c)
2648 REAL(c_double),
DIMENSION(*),
INTENT(OUT) :: bufferD
2649 REAL(c_float),
DIMENSION(*),
INTENT(OUT) :: bufferF
2650 INTEGER(c_int),
DIMENSION(*),
INTENT(OUT) :: bufferI
2651 INTEGER(c_int),
INTENT(INOUT) :: bufferLength
2652 INTEGER(c_int),
INTENT(IN),
VALUE :: lun
2653 INTEGER(c_int),
INTENT(OUT) :: err
2654 END SUBROUTINE readc
2660 DATA npri / 0 /, mpri / 1000 /
2708 files:
DO WHILE (jfile > 0)
2712 CALL binopn(kfile,ithr,ios)
2718 IF(kfile <=
nfilf)
THEN
2741 eof=(ierrc <= 0.AND.ierrc /= -4)
2742 IF(eof.AND.ierrc < 0)
THEN
2743 WRITE(*,*)
'Read error for binary Cfile', kfile,
'record', jrec+1,
':', ierrc
2744 WRITE(8,*)
'Read error for binary Cfile', kfile,
'record', jrec+1,
':', ierrc
2746 WRITE(cfile,
'(I7)') kfile
2747 CALL peend(18,
'Aborted, read error(s) for binary file ' // cfile)
2748 stop
'PEREAD: stopping due to read errors (bad record, wrong file type?)'
2750 IF (
kfd(1,jfile) == 1)
THEN
2756 IF(eof)
EXIT records
2761 xfd(jfile)=max(
xfd(jfile),n)
2778 IF ((noff+nr+2+
ndimbuf >= ndata*(iact+1)).OR.(nbuf >= npointer))
EXIT files
2793 IF (
kfd(1,jfile) == 1)
THEN
2794 print *,
'PEREAD: file ', kfile,
'read the first time, found',jrec,
' records'
2798 IF (-
kfd(1,jfile) /= jrec)
THEN
2799 WRITE(cfile,
'(I7)') kfile
2800 CALL peend(19,
'Aborted, binary file modified (length) ' // cfile)
2801 stop
'PEREAD: file modified (length)'
2852 IF (
kfd(1,jfile) == 1)
kfd(1,jfile)=-nrc
2871 DO WHILE (
nloopn == 0.AND.ntot >= nrpr)
2872 WRITE(*,*)
' Record ',nrpr
2873 IF (nrpr < 100000)
THEN
2882 IF (npri == 1)
WRITE(*,100)
2884100
FORMAT(/
' PeRead records active file' &
2885 /
' total block threads number')
2886101
FORMAT(
' PeRead',4i10)
2899 IF (
xfd(k) > maxrecordsize)
THEN
2900 maxrecordsize=
xfd(k)
2903 dw=real(-
kfd(1,k),mpd)
2906 ds1=ds1+dw*real(
wfd(k),mpd)
2907 ds2=ds2+dw*real(
wfd(k)**2,mpd)
2909 print *,
'PEREAD: file ', maxrecordfile,
'with max record size ', maxrecordsize
2910 IF (
nfilw > 0.AND.ds0 > 0.0_mpd)
THEN
2914 WRITE(lun,177)
nfilw,real(ds1,mps),real(ds2,mps)
2915177
FORMAT(/
' !!!!!',i4,
' weighted binary files', &
2916 /
' !!!!! mean, variance of weights =',2g12.4)
2934179
FORMAT(/
' Read cache usage (#blocks, #records, ', &
2935 'min,max records/block'/17x,i10,i12,2i10)
2953 INTEGER(mpi),
INTENT(IN) :: mode
2955 INTEGER(mpi) :: ibuf
2956 INTEGER(mpi) :: ichunk
2958 INTEGER(mpi) :: itgbi
2964 INTEGER(mpi),
PARAMETER :: maxbad = 100
2965 INTEGER(mpi) :: nbad
2966 INTEGER(mpi) :: nerr
2967 INTEGER(mpi) :: inone
2985 CALL isjajb(nst,ist,ja,jb,jsp)
3004#ifdef SCOREP_USER_ENABLE
3005 scorep_user_region_by_name_begin(
"UR_peprep", scorep_user_region_type_common)
3010 CALL pechk(ibuf,nerr)
3013 IF(nbad >= maxbad)
EXIT
3018 CALL isjajb(nst,ist,ja,jb,jsp)
3031 CALL peend(20,
'Aborted, bad binary records')
3032 stop
'PEREAD: stopping due to bad records'
3035#ifdef SCOREP_USER_ENABLE
3036 scorep_user_region_by_name_end(
"UR_peprep")
3056 INTEGER(mpi) :: ioff
3063 INTEGER(mpi),
INTENT(IN) :: ibuf
3064 INTEGER(mpi),
INTENT(OUT) :: nerr
3073 outer:
DO WHILE(is < nst)
3078 IF(is > nst)
EXIT outer
3084 IF(is > nst)
EXIT outer
3101100
FORMAT(
' PEREAD: record ', i8,
' in file ',i6,
' is broken !!!')
3111101
FORMAT(
' PEREAD: record ', i8,
' in file ',i6,
' contains ', i6,
' NaNs !!!')
3127 INTEGER(mpi) :: ibuf
3128 INTEGER(mpi) :: ichunk
3129 INTEGER(mpi) :: iproc
3130 INTEGER(mpi) :: ioff
3131 INTEGER(mpi) :: ioffbi
3133 INTEGER(mpi) :: itgbi
3138 INTEGER(mpi) :: nalg
3139 INTEGER(mpi) :: neqna
3142 INTEGER(mpi) :: nzero
3143 INTEGER(mpi) :: inone
3144 INTEGER(mpl) :: length
3179 CALL isjajb(nst,ist,ja,jb,jsp)
3210 CALL isjajb(nst,ist,ja,jb,jsp)
3238#ifdef SCOREP_USER_ENABLE
3239 scorep_user_region_by_name_begin(
"UR_pepgrp", scorep_user_region_type_common)
3248 CALL pargrp(ist+1,ist+nnz)
3260#ifdef SCOREP_USER_ENABLE
3261 scorep_user_region_by_name_end(
"UR_pepgrp")
3280 INTEGER(mpi) :: istart
3281 INTEGER(mpi) :: itgbi
3283 INTEGER(mpi) :: jstart
3284 INTEGER(mpi) :: jtgbi
3285 INTEGER(mpi) :: lstart
3286 INTEGER(mpi) :: ltgbi
3288 INTEGER(mpi),
INTENT(IN) :: inds
3289 INTEGER(mpi),
INTENT(IN) :: inde
3291 IF (inds > inde)
RETURN
3300 IF (istart == 0)
THEN
3301 IF (itgbi /= ltgbi+1)
THEN
3304 IF (lstart == 0)
THEN
3319 IF (istart /= jstart)
THEN
3330 IF (itgbi /= ltgbi+1)
THEN
3345 IF (istart /= jstart)
THEN
3397 INTEGER(mpi),
INTENT(IN) :: nst
3398 INTEGER(mpi),
INTENT(IN OUT) :: is
3399 INTEGER(mpi),
INTENT(OUT) :: ja
3400 INTEGER(mpi),
INTENT(OUT) :: jb
3401 INTEGER(mpi),
INTENT(OUT) :: jsp
3409 IF(is >= nst)
RETURN
3463 INTEGER(mpi) :: ioffb
3465 INTEGER(mpi) :: itgbi
3466 INTEGER(mpi) :: itgbij
3467 INTEGER(mpi) :: itgbik
3468 INTEGER(mpi) :: ivgb
3469 INTEGER(mpi) :: ivgbij
3470 INTEGER(mpi) :: ivgbik
3473 INTEGER(mpi) :: lastit
3475 INTEGER(mpi) :: ncrit
3476 INTEGER(mpi) :: ngras
3477 INTEGER(mpi) :: nparl
3479 INTEGER(mpl) :: nrej
3480 INTEGER(mpi) :: inone
3481 INTEGER(mpi) :: ilow
3482 INTEGER(mpi) :: nlow
3483 INTEGER(mpi) :: nzero
3503 CALL gmpdef(1,4,
'Function value in iterations')
3505 CALL gmpdef(2,3,
'Number of MINRES steps vs iteration nr')
3507 CALL hmpdef( 5,0.0,0.0,
'Number of degrees of freedom')
3508 CALL hmpdef(11,0.0,0.0,
'Number of local parameters')
3509 CALL hmpdef(16,0.0,24.0,
'LOG10(cond(band part decomp.)) local fit ')
3510 CALL hmpdef(23,0.0,0.0,
'SQRT of diagonal elements without presigma')
3511 CALL hmpdef(24,0.0,0.0,
'Log10 of off-diagonal elements')
3512 CALL hmpdef(25,0.0,0.0,
'Relative individual pre-sigma')
3513 CALL hmpdef(26,0.0,0.0,
'Relative global pre-sigma')
3518 'Normalized residuals of single (global) measurement')
3520 'Normalized residuals of single (local) measurement')
3522 'Pulls of single (global) measurement')
3524 'Pulls of single (local) measurement')
3525 CALL hmpdef(4,0.0,0.0,
'Chi^2/Ndf after local fit')
3526 CALL gmpdef(4,5,
'location, dispersion (res.) vs record nr')
3527 CALL gmpdef(5,5,
'location, dispersion (pull) vs record nr')
3541 IF(
nloopn == 2)
CALL hmpdef(6,0.0,0.0,
'Down-weight fraction')
3544 IF(
iterat /= lastit)
THEN
3552 ELSE IF(
iterat >= 1)
THEN
3600 WRITE(*,111) nparl,ncrit,usei,used,peaki,peakd
3601111
FORMAT(
' Write cache usage (#flush,#overrun,<levels>,', &
3602 'peak(levels))'/2i7,
',',4(f6.1,
'%'))
3610 ELSE IF (
mbandw > 0)
THEN
3642 print *,
" ... warning ..."
3643 print *,
" global parameters with too few (< MREQENA) accepted entries: ", nlow
3647 IF(
icalcm == 1 .AND. nzero > 0)
THEN
3649 WRITE(*,*)
'Warning: the rank defect of the symmetric',
nfgb, &
3650 '-by-',
nfgb,
' matrix is ',
ndefec,
' (should be zero).'
3651 WRITE(lun,*)
'Warning: the rank defect of the symmetric',
nfgb, &
3652 '-by-',
nfgb,
' matrix is ',
ndefec,
' (should be zero).'
3655 WRITE(*,*)
' --> enforcing SUBITO mode'
3656 WRITE(lun,*)
' --> enforcing SUBITO mode'
3666 elmt=real(matij(i,i),mps)
3667 IF(elmt > 0.0)
CALL hmpent(23,1.0/sqrt(elmt))
3701 CALL addsums(1, adder, 0, 1.0_mpl)
3715 IF(sgm > 0.0) weight=1.0/sgm**2
3731 IF(itgbij /= 0)
THEN
3737 IF (factrj == 0.0_mpd) cycle
3747 IF(
icalcm == 1.AND.ivgbij > 0)
THEN
3754 IF(ivgbij > 0.AND.ivgbik > 0)
THEN
3755 CALL mupdat(ivgbij,ivgbik,weight*factrj*factrk)
3761 adder=weight*dsum**2
3762 CALL addsums(1, adder, 1, 1.0_mpl)
3780 ratae=max(-99.9,ratae)
3789 WRITE(*,*)
'Data records rejected in initial loop:'
3834 WRITE(lun,*)
' === local fits have bordered band matrix structure ==='
3835 IF (
nbndr(1) > 0 )
WRITE(lun,101)
' NBNDR',
nbndr(1),
'number of records (upper/left border)'
3836 IF (
nbndr(2) > 0 )
WRITE(lun,101)
' NBNDR',
nbndr(2),
'number of records (lower/right border)'
3837 WRITE(lun,101)
' NBDRX',
nbdrx,
'max border size'
3838 WRITE(lun,101)
' NBNDX',
nbndx,
'max band width'
3847101
FORMAT(1x,a8,
' =',i14,
' = ',a)
3862 INTEGER(mpi),
INTENT(IN) :: lunp
3867101
FORMAT(
' it fc',
' fcn_value dfcn_exp slpr costh iit st', &
3868 ' ls step cutf',1x,
'rejects hhmmss FMS')
3869102
FORMAT(
' -- --',
' ----------- -------- ---- ----- --- --', &
3870 ' -- ----- ----',1x,
'------- ------ ---')
3886 INTEGER(mpl) :: nrej
3887 INTEGER(mpi) :: nsecnd
3891 REAL(mps) :: slopes(3)
3892 REAL(mps) :: steps(3)
3893 REAL,
DIMENSION(2) :: ta
3896 INTEGER(mpi),
INTENT(IN) :: lunp
3898 CHARACTER (LEN=4):: ccalcm(4)
3899 DATA ccalcm /
' end',
' S',
' F ',
' FMS' /
3903 IF(nrej > 9999999) nrej=9999999
3907 nsecnd=nint(secnd,mpi)
3916 CALL ptlopt(nfa,ma,slopes,steps)
3917 ratae=max(-99.9,min(99.9,slopes(2)/slopes(1)))
3923103
FORMAT(i3,i3,e12.5,38x,f5.1, 1x,i7, i3,i2.2,i2.2,a4)
3924104
FORMAT(i3,i3,e12.5,1x,e8.2,f6.3,f6.3,i5,2i3,f6.3,f5.1, &
3925 1x,i7, i3,i2.2,i2.2,a4)
3926105
FORMAT(i3,i3,e12.5,1x,e8.2,12x,i5,i3,9x,f5.1, &
3927 1x,i7, i3,i2.2,i2.2,a4)
3940 INTEGER(mpi) :: minut
3942 INTEGER(mpi) :: nhour
3943 INTEGER(mpl) :: nrej
3944 INTEGER(mpi) :: nsecnd
3948 REAL(mps) :: slopes(3)
3949 REAL(mps) :: steps(3)
3950 REAL,
DIMENSION(2) :: ta
3953 INTEGER(mpi),
INTENT(IN) :: lunp
3954 CHARACTER (LEN=4):: ccalcm(4)
3955 DATA ccalcm /
' end',
' S',
' F ',
' FMS' /
3959 IF(nrej > 9999999) nrej=9999999
3963 nsecnd=nint(secnd,mpi)
3967 CALL ptlopt(nfa,ma,slopes,steps)
3968 ratae=abs(slopes(2)/slopes(1))
3973104
FORMAT(3x,i3,e12.5,9x, 35x, i7, i3,i2.2,i2.2,a4)
3974105
FORMAT(3x,i3,e12.5,9x, f6.3,14x,i3,f6.3,6x, i7, i3,i2.2,i2.2,a4)
3988 INTEGER(mpi) :: nsecnd
3991 REAL,
DIMENSION(2) :: ta
3994 INTEGER(mpi),
INTENT(IN) :: lunp
3995 CHARACTER (LEN=4):: ccalcm(4)
3996 DATA ccalcm /
' end',
' S',
' F ',
' FMS' /
4001 nsecnd=nint(secnd,mpi)
4003 WRITE(lunp,106) nhour,minut,nsecnd,ccalcm(
lcalcm)
4004106
FORMAT(69x,i3,i2.2,i2.2,a4)
4014 INTEGER(mpi) :: lunit
4016 WRITE(lunit,102)
'Explanation of iteration table'
4017 WRITE(lunit,102)
'=============================='
4018 WRITE(lunit,101)
'it', &
4019 'iteration number. Global parameters are improved for it > 0.'
4020 WRITE(lunit,102)
'First function evaluation is called iteraton 0.'
4021 WRITE(lunit,101)
'fc',
'number of function evaluations.'
4022 WRITE(lunit,101)
'fcn_value',
'value of 2 x Likelihood function (LF).'
4023 WRITE(lunit,102)
'The final value is the chi^2 value of the fit and should'
4024 WRITE(lunit,102)
'be about equal to the NDF (see below).'
4025 WRITE(lunit,101)
'dfcn_exp', &
4026 'expected reduction of the value of the Likelihood function (LF)'
4027 WRITE(lunit,101)
'slpr',
'ratio of the actual slope to inital slope.'
4028 WRITE(lunit,101)
'costh', &
4029 'cosine of angle between search direction and -gradient'
4031 WRITE(lunit,101)
'iit', &
4032 'number of internal iterations in MINRES algorithm'
4033 WRITE(lunit,101)
'st',
'stop code of MINRES algorithm'
4034 WRITE(lunit,102)
'< 0: rhs is very special, with beta2 = 0'
4035 WRITE(lunit,102)
'= 0: rhs b = 0, i.e. the exact solution is x = 0'
4036 WRITE(lunit,102)
'= 1 requested accuracy achieved, as determined by rtol'
4037 WRITE(lunit,102)
'= 2 reasonable accuracy achieved, given eps'
4038 WRITE(lunit,102)
'= 3 x has converged to an eigenvector'
4039 WRITE(lunit,102)
'= 4 matrix ill-conditioned (Acond has exceeded 0.1/eps)'
4040 WRITE(lunit,102)
'= 5 the iteration limit was reached'
4041 WRITE(lunit,102)
'= 6 Matrix x vector does not define a symmetric matrix'
4042 WRITE(lunit,102)
'= 7 Preconditioner does not define a symmetric matrix'
4043 ELSEIF (
metsol == 5)
THEN
4044 WRITE(lunit,101)
'iit', &
4045 'number of internal iterations in MINRES-QLP algorithm'
4046 WRITE(lunit,101)
'st',
'stop code of MINRES-QLP algorithm'
4047 WRITE(lunit,102)
'= 1: beta_{k+1} < eps, iteration k is the final Lanczos step.'
4048 WRITE(lunit,102)
'= 2: beta2 = 0. If M = I, b and x are eigenvectors of A.'
4049 WRITE(lunit,102)
'= 3: beta1 = 0. The exact solution is x = 0.'
4050 WRITE(lunit,102)
'= 4: A solution to (poss. singular) Ax = b found, given rtol.'
4051 WRITE(lunit,102)
'= 5: A solution to (poss. singular) Ax = b found, given eps.'
4052 WRITE(lunit,102)
'= 6: Pseudoinverse solution for singular LS problem, given rtol.'
4053 WRITE(lunit,102)
'= 7: Pseudoinverse solution for singular LS problem, given eps.'
4054 WRITE(lunit,102)
'= 8: The iteration limit was reached.'
4055 WRITE(lunit,102)
'= 9: The operator defined by Aprod appears to be unsymmetric.'
4056 WRITE(lunit,102)
'=10: The operator defined by Msolve appears to be unsymmetric.'
4057 WRITE(lunit,102)
'=11: The operator defined by Msolve appears to be indefinite.'
4058 WRITE(lunit,102)
'=12: xnorm has exceeded maxxnorm or will exceed it next iteration.'
4059 WRITE(lunit,102)
'=13: Acond has exceeded Acondlim or 0.1/eps.'
4060 WRITE(lunit,102)
'=14: Least-squares problem but no converged solution yet.'
4061 WRITE(lunit,102)
'=15: A null vector obtained, given rtol.'
4063 WRITE(lunit,101)
'ls',
'line search info'
4064 WRITE(lunit,102)
'< 0 recalculate function'
4065 WRITE(lunit,102)
'= 0: N or STP lt 0 or step not descending'
4066 WRITE(lunit,102)
'= 1: Linesearch convergence conditions reached'
4067 WRITE(lunit,102)
'= 2: interval of uncertainty at lower limit'
4068 WRITE(lunit,102)
'= 3: max nr of line search calls reached'
4069 WRITE(lunit,102)
'= 4: step at the lower bound'
4070 WRITE(lunit,102)
'= 5: step at the upper bound'
4071 WRITE(lunit,102)
'= 6: rounding error limitation'
4072 WRITE(lunit,101)
'step', &
4073 'the factor for the Newton step during the line search. Usually'
4075 'a value of 1 gives a sufficient reduction of the LF. Oherwise'
4076 WRITE(lunit,102)
'other step values are tried.'
4077 WRITE(lunit,101)
'cutf', &
4078 'cut factor. Local fits are rejected, if their chi^2 value'
4080 'is larger than the 3-sigma chi^2 value times the cut factor.'
4081 WRITE(lunit,102)
'A cut factor of 1 is used finally, but initially a larger'
4082 WRITE(lunit,102)
'factor may be used. A value of 0.0 means no cut.'
4083 WRITE(lunit,101)
'rejects',
'total number of rejected local fits.'
4084 WRITE(lunit,101)
'hmmsec',
'the time in hours (h), minutes (mm) and seconds.'
4085 WRITE(lunit,101)
'FMS',
'calculation of Function value, Matrix, Solution.'
4088101
FORMAT(a9,
' = ',a)
4105 INTEGER(mpi),
INTENT(IN) :: i
4106 INTEGER(mpi),
INTENT(IN) :: j
4107 REAL(mpd),
INTENT(IN) :: add
4109 INTEGER(mpl):: ijadd
4110 INTEGER(mpl):: ijcsr3
4115 IF(i <= 0.OR.j <= 0.OR. add == 0.0_mpd)
RETURN
4136 ELSE IF(
matsto == 2)
THEN
4164 CALL peend(23,
'Aborted, bad matrix index')
4165 stop
'mupdat: bad index'
4190 INTEGER(mpi),
INTENT(IN) :: i
4191 INTEGER(mpi),
INTENT(IN) :: j1
4192 INTEGER(mpi),
INTENT(IN) :: j2
4193 INTEGER(mpi),
INTENT(IN) :: il
4194 INTEGER(mpi),
INTENT(IN) :: jl
4195 INTEGER(mpi),
INTENT(IN) :: n
4196 REAL(mpd),
INTENT(IN) :: sub((n*n+n)/2)
4203 INTEGER(mpi):: iblast
4204 INTEGER(mpi):: iblock
4210 INTEGER(mpi):: jblast
4211 INTEGER(mpi):: jblock
4221 IF(i <= 0.OR.j1 <= 0.OR.j2 > i)
RETURN
4235 print *,
' MGUPDT: ij=0', i,j1,j2,il,jl,ij,lr,iprc,
matsto
4242 ijl=(k*k-k)/2+jl+ir-ia1
4259 ijl=(k*k-k)/2+jl+ir-ia1
4263 IF (jblock /= jblast.OR.iblock /= iblast)
THEN
4287 CALL ijpgrp(i,j1,ij,lr,iprc)
4289 print *,
' MGUPDT: ij=0', i,j1,j2,il,jl,ij,lr,iprc,
matsto
4353SUBROUTINE loopbf(nrej,numfil,naccf,chi2f,ndff)
4380 INTEGER(mpi) :: ibuf
4381 INTEGER(mpi) :: ichunk
4382 INTEGER(mpl) :: icmn
4383 INTEGER(mpl) :: icost
4385 INTEGER(mpi) :: idiag
4387 INTEGER(mpi) :: iext
4395 INTEGER(mpi) :: imeas
4398 INTEGER(mpi) :: ioffb
4399 INTEGER(mpi) :: ioffc
4400 INTEGER(mpi) :: ioffd
4401 INTEGER(mpi) :: ioffe
4402 INTEGER(mpi) :: ioffi
4403 INTEGER(mpi) :: ioffq
4404 INTEGER(mpi) :: iprc
4405 INTEGER(mpi) :: iprcnx
4406 INTEGER(mpi) :: iprdbg
4407 INTEGER(mpi) :: iproc
4408 INTEGER(mpi) :: irbin
4409 INTEGER(mpi) :: isize
4411 INTEGER(mpi) :: iter
4412 INTEGER(mpi) :: itgbi
4413 INTEGER(mpi) :: ivgbj
4414 INTEGER(mpi) :: ivgbk
4415 INTEGER(mpi) :: ivpgrp
4425 INTEGER(mpi) :: joffd
4426 INTEGER(mpi) :: joffi
4427 INTEGER(mpi) :: jproc
4431 INTEGER(mpi) :: kbdr
4432 INTEGER(mpi) :: kbdrx
4433 INTEGER(mpi) :: kbnd
4436 INTEGER(mpi) :: lvpgrp
4437 INTEGER(mpi) :: mbdr
4438 INTEGER(mpi) :: mbnd
4439 INTEGER(mpi) :: mside
4440 INTEGER(mpi) :: nalc
4441 INTEGER(mpi) :: nalg
4445 INTEGER(mpi) :: ndown
4447 INTEGER(mpi) :: nfred
4448 INTEGER(mpi) :: nfrei
4450 INTEGER(mpi) :: nprdbg
4451 INTEGER(mpi) :: nrank
4454 INTEGER(mpi) :: nter
4455 INTEGER(mpi) :: nweig
4456 INTEGER(mpi) :: ngrp
4457 INTEGER(mpi) :: npar
4459 INTEGER(mpl),
INTENT(IN OUT) :: nrej(6)
4460 INTEGER(mpi),
INTENT(IN) :: numfil
4461 INTEGER(mpi),
INTENT(IN OUT) :: naccf(numfil)
4462 REAL(mps),
INTENT(IN OUT) :: chi2f(numfil)
4463 INTEGER(mpi),
INTENT(IN OUT) :: ndff(numfil)
4473 INTEGER(mpi) :: ijprec
4480 CHARACTER (LEN=3):: chast
4481 DATA chuber/1.345_mpd/
4482 DATA cauchy/2.3849_mpd/
4530 IF(
nloopn == 1.AND.mod(nrc,100000_mpl) == 0)
THEN
4531 WRITE(*,*)
'Record',nrc,
' ... still reading'
4532 IF(
monpg1>0)
WRITE(
lunlog,*)
'Record',nrc,
' ... still reading'
4542 IF(nrc ==
nrecpr) lprnt=.true.
4543 IF(nrc ==
nrecp2) lprnt=.true.
4544 IF(nrc ==
nrecer) lprnt=.true.
4549 IF (nprdbg == 1) iprdbg=iproc
4550 IF (iproc /= iprdbg) lprnt=.false.
4555 WRITE(1,*)
'------------------ Loop',
nloopn, &
4556 ': Printout for record',nrc,iproc
4567 CALL isjajb(nst,ist,ja,jb,jsp)
4569 IF(imeas == 0)
WRITE(1,1121)
45741121
FORMAT(/
'Measured value and local derivatives'/ &
4575 ' i measured std_dev index...derivative ...')
45761122
FORMAT(i3,2g12.4,3(i3,g12.4)/(27x,3(i3,g12.4)))
4582 CALL isjajb(nst,ist,ja,jb,jsp)
4584 IF(imeas == 0)
WRITE(1,1123)
45961123
FORMAT(/
'Global derivatives'/ &
4597 ' i label gindex vindex derivative ...')
45981124
FORMAT(i3,2(i9,i7,i7,g12.4)/(3x,2(i9,i7,i7,g12.4)))
45991125
FORMAT(i3,2(i9,i7,i7,g12.4))
4609 WRITE(1,*)
'Data corrections using values of global parameters'
4610 WRITE(1,*)
'=================================================='
4620 CALL isjajb(nst,ist,ja,jb,jsp)
4652101
FORMAT(
' index measvalue corrvalue sigma')
4653102
FORMAT(i6,2x,2g12.4,
' +-',g12.4)
4655 IF(nalc <= 0)
GO TO 90
4657 ngg=(nalg*nalg+nalg)/2
4670 IF (ivpgrp /= lvpgrp)
THEN
4678 IF (npar /= nalg)
THEN
4679 print *,
' mismatch of number of global parameters ', nrc, nalg, npar, ngrp
4689 CALL peend(35,
'Aborted, mismatch of number of global parameters')
4690 stop
' mismatch of number of global parameters '
4717 DO WHILE(iter < nter)
4722 WRITE(1,*)
'Outlier-suppression iteration',iter,
' of',nter
4723 WRITE(1,*)
'=========================================='
4733 DO i=1,(nalc*nalc+nalc)/2
4744 wght =1.0_mpd/rerr**2
4749 IF(abs(resid) > chuber*rerr)
THEN
4750 wght=wght*chuber*rerr/abs(resid)
4754 wght=wght/(1.0+(resid/rerr/cauchy)**2)
4758 IF(lprnt.AND.iter /= 1.AND.nter /= 1)
THEN
4760 IF(abs(resid) > chuber*rerr) chast=
'* '
4761 IF(abs(resid) > 3.0*rerr) chast=
'** '
4762 IF(abs(resid) > 6.0*rerr) chast=
'***'
4763 IF(imeas == 0)
WRITE(1,*)
'Second loop: accumulate'
4764 IF(imeas == 0)
WRITE(1,103)
4769 WRITE(1,104) imeas,rmeas,resid,rerr,r1,chast,r2
4771103
FORMAT(
' index corrvalue residuum sigma', &
4773104
FORMAT(i6,2x,2g12.4,
' +-',g12.4,f7.2,1x,a3,f8.2)
4788 im=min(nalc+1-ij,nalc+1-ik)
4795 IF (iter == 1.AND.nalc > 5.AND.
lfitbb > 0)
THEN
4798 icmn=int(nalc,mpl)**3
4804 icost=6*int(nalc-kbdr,mpl)*int(kbnd+kbdr+1,mpl)**2+2*int(kbdr,mpl)**3
4805 IF (icost < icmn)
THEN
4817 kbdr=max(kbdr,
ibandh(k+nalc))
4818 icost=6*int(nalc-kbdr,mpl)*int(kbnd+kbdr+1,mpl)**2+2*int(kbdr,mpl)**3
4819 IF (icost < icmn)
THEN
4843 IF (
icalcm == 1.OR.lprnt) inv=2
4844 IF (mside == 1)
THEN
4852 IF (evdmin > 0.0_mpl) cndl10=log10(real(evdmax/evdmin,mps))
4862 IF ((.NOT.(
blvec(k) <= 0.0_mpd)).AND. (.NOT.(
blvec(k) > 0.0_mpd))) nan=nan+1
4867 WRITE(1,*)
'Parameter determination:',nalc,
' parameters,',
' rank=',nrank
4868 WRITE(1,*)
'-----------------------'
4869 IF(ndown /= 0)
WRITE(1,*)
' ',ndown,
' data down-weighted'
4885 wght =1.0_mpd/rerr**2
4909 IF (0.999999_mpd/wght > dvar)
THEN
4910 pull=rmeas/sqrt(1.0_mpd/wght-dvar)
4913 CALL hmpent(13,real(pull,mps))
4914 CALL gmpms(5,rec,real(pull,mps))
4916 CALL hmpent(14,real(pull,mps))
4935 IF(iter == 1.AND.jb < ist.AND.lhist) &
4936 CALL gmpms(4,rec,real(rmeas/rerr,mps))
4938 dchi2=wght*rmeas*rmeas
4943 IF(abs(resid) > chuber*rerr)
THEN
4944 wght=wght*chuber*rerr/abs(resid)
4945 dchi2=2.0*chuber*(abs(resid)/rerr-0.5*chuber)
4948 wght=wght/(1.0_mpd+(resid/rerr/cauchy)**2)
4949 dchi2=log(1.0_mpd+(resid/rerr/cauchy)**2)*cauchy**2
4959 IF(abs(resid) > chuber*rerr) chast=
'* '
4960 IF(abs(resid) > 3.0*rerr) chast=
'** '
4961 IF(abs(resid) > 6.0*rerr) chast=
'***'
4962 IF(imeas == 0)
WRITE(1,*)
'Third loop: single residuals'
4963 IF(imeas == 0)
WRITE(1,105)
4967 IF(resid < 0.0) r1=-r1
4968 IF(resid < 0.0) r2=-r2
4971105
FORMAT(
' index corrvalue residuum sigma', &
4973106
FORMAT(i6,2x,2g12.4,
' +-',g12.4,f7.2,1x,a3,f8.2)
4975 IF(iter == nter)
THEN
4977 resmax=max(resmax,abs(rmeas)/rerr)
4980 IF(iter == 1.AND.lhist)
THEN
4982 CALL hmpent( 3,real(rmeas/rerr,mps))
4984 CALL hmpent(12,real(rmeas/rerr,mps))
4991 resing=(real(nweig,mps)-real(suwt,mps))/real(nweig,mps)
4993 IF(iter == 1)
CALL hmpent( 5,real(ndf,mps))
4994 IF(iter == 1)
CALL hmpent(11,real(nalc,mps))
4999 WRITE(1,*)
'Chi^2=',summ,
' at',ndf,
' degrees of freedom: ', &
5000 '3-sigma limit is',chindl(3,ndf)*real(ndf,mps)
5001 WRITE(1,*) suwt,
' is sum of factors, compared to',nweig, &
5002 ' Downweight fraction:',resing
5008 WRITE(1,*)
' NaNs ', nalc, nrank, nan
5009 WRITE(1,*)
' ---> rejected!'
5011 IF (
mdebug < 0.AND.
nloopn == 1) print *,
' bad local fit-1 ', kfl,jrc,nrc,mside,neq,nalc,nrank,nan,cndl10
5014 IF(nrank /= nalc)
THEN
5018 WRITE(1,*)
' rank deficit', nalc, nrank
5019 WRITE(1,*)
' ---> rejected!'
5021 IF (
mdebug < 0.AND.
nloopn == 1) print *,
' bad local fit-2 ', kfl,jrc,nrc,mside,neq,nalc,nrank,nan,cndl10
5028 WRITE(1,*)
' too large condition(band part) ', nalc, nrank, cndl10
5029 WRITE(1,*)
' ---> rejected!'
5031 IF (
mdebug < 0.AND.
nloopn == 1) print *,
' bad local fit-3 ', kfl,jrc,nrc,mside,neq,nalc,nrank,nan,cndl10
5037 WRITE(1,*)
' Ndf<=0', nalc, nrank, ndf
5038 WRITE(1,*)
' ---> rejected!'
5040 IF (
mdebug < 0.AND.
nloopn == 1) print *,
' bad local fit-4 ', kfl,jrc,nrc,mside,neq,nalc,nrank,nan,cndl10
5044 chndf=real(summ/real(ndf,mpd),mps)
5046 IF(iter == 1.AND.lhist)
CALL hmpent(4,chndf)
5059 chichi=chindl(3,ndf)*real(ndf,mps)
5063 IF(summ >
chhuge*chichi)
THEN
5066 WRITE(1,*)
' ---> rejected!'
5074 IF(summ > chlimt)
THEN
5076 WRITE(1,*)
' ---> rejected!'
5080 CALL addsums(iproc+1, dchi2, ndf, dw1)
5090 CALL addsums(iproc+1, dchi2, ndf, dw1)
5094 WRITE(1,*)
' ---> rejected!'
5107 naccf(kfl)=naccf(kfl)+1
5108 ndff(kfl) =ndff(kfl) +ndf
5109 chi2f(kfl)=chi2f(kfl)+chndf
5122 wght =1.0_mpd/rerr**2
5123 dchi2=wght*rmeas*rmeas
5127 IF(resid > chuber*rerr)
THEN
5128 wght=wght*chuber*rerr/resid
5129 dchi2=2.0*chuber*(resid/rerr-0.5*chuber)
5179 CALL addsums(iproc+1, summ, ndf, dw1)
5235 IF (nfred < 0.OR.nfrei < 0)
THEN
5262 iprcnx=ijprec(i,jnx)
5264 IF (.NOT.((jnx == j+1).AND.(iprc == iprcnx)))
THEN
5278 joffd=joffd+(il*il-il)/2
5290 WRITE(1,*)
'------------------ End of printout for record',nrc
5347 iprcnx=ijprec(i,jnx)
5349 IF (.NOT.((jnx == j+1).AND.(iprc == iprcnx)))
THEN
5364 joffd=joffd+(il*il-il)/2
5400 INTEGER(mpi),
INTENT(IN) :: lun
5402 IF (
nrejec(1)>0)
WRITE(lun,*)
nrejec(1),
' (local solution contains NaNs)'
5403 IF (
nrejec(2)>0)
WRITE(lun,*)
nrejec(2),
' (local matrix with rank deficit)'
5404 IF (
nrejec(3)>0)
WRITE(lun,*)
nrejec(3),
' (local matrix with ill condition)'
5405 IF (
nrejec(4)>0)
WRITE(lun,*)
nrejec(4),
' (local fit with Ndf=0)'
5406 IF (
nrejec(5)>0)
WRITE(lun,*)
nrejec(5),
' (local fit with huge Chi2(Ndf))'
5407 IF (
nrejec(6)>0)
WRITE(lun,*)
nrejec(6),
' (local fit with large Chi2(Ndf))'
5433 INTEGER(mpi) :: icom
5434 INTEGER(mpl) :: icount
5438 INTEGER(mpi) :: imin
5439 INTEGER(mpi) :: iprlim
5440 INTEGER(mpi) :: isub
5441 INTEGER(mpi) :: itgbi
5442 INTEGER(mpi) :: itgbl
5443 INTEGER(mpi) :: ivgbi
5445 INTEGER(mpi) :: label
5453 INTEGER(mpi) :: labele(3)
5454 REAL(mps):: compnt(3)
5459 CALL mvopen(lup,
'millepede.res')
5462 WRITE(*,*)
' Result of fit for global parameters'
5463 WRITE(*,*)
' ==================================='
5468 WRITE(lup,*)
'Parameter ! first 3 elements per line are', &
5469 ' significant (if used as input)'
5487 err=sqrt(abs(real(gmati,mps)))
5488 IF(gmati < 0.0_mpd) err=-err
5491 IF(gmati*diag > 0.0_mpd)
THEN
5492 gcor2=1.0_mpd-1.0_mpd/(gmati*diag)
5493 IF(gcor2 >= 0.0_mpd.AND.gcor2 <= 1.0_mpd) gcor=real(sqrt(gcor2),mps)
5498 IF(lowstat) icount=-(icount+1)
5500 IF(itgbi <= iprlim)
THEN
5514 ELSE IF(itgbi == iprlim+1)
THEN
5515 WRITE(* ,*)
'... (further printout suppressed, but see log file)'
5529 ELSE IF (
igcorr /= 0)
THEN
5547 CALL mvopen(lup,
'millepede.eve')
5557 DO isub=0,min(15,imin-1)
5579 WRITE(lup,103) (labele(ie),compnt(ie),ie=1,iev)
5583 IF(iev /= 0)
WRITE(lup,103) (labele(ie),compnt(ie),ie=1,iev)
5591101
FORMAT(1x,
' label parameter presigma differ', &
5592 ' error'/ 1x,
'-----------',4x,4(
'-------------'))
5593102
FORMAT(i10,2x,4g14.5,f8.3)
5594103
FORMAT(3(i11,f11.7,2x))
5595110
FORMAT(i10,2x,2g14.5,28x,i12)
5596111
FORMAT(i10,2x,3g14.5,14x,i12)
5597112
FORMAT(i10,2x,4g14.5,i12)
5619 INTEGER(mpi) :: icom
5620 INTEGER(mpl) :: icount
5621 INTEGER(mpi) :: ifrst
5622 INTEGER(mpi) :: ilast
5623 INTEGER(mpi) :: inext
5624 INTEGER(mpi) :: itgbi
5625 INTEGER(mpi) :: itgbl
5626 INTEGER(mpi) :: itpgrp
5627 INTEGER(mpi) :: ivgbi
5629 INTEGER(mpi) :: icgrp
5630 INTEGER(mpi) :: ipgrp
5632 INTEGER(mpi) :: jpgrp
5634 INTEGER(mpi) :: label1
5635 INTEGER(mpi) :: label2
5636 INTEGER(mpi) :: ncon
5637 INTEGER(mpi) :: npair
5638 INTEGER(mpi) :: nstep
5641 INTEGER(mpl):: length
5643 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: vecPairedParGroups
5646 SUBROUTINE ggbmap(ipgrp,npair,npgrp)
5648 INTEGER(mpi),
INTENT(IN) :: ipgrp
5649 INTEGER(mpi),
INTENT(OUT) :: npair
5650 INTEGER(mpi),
DIMENSION(:),
INTENT(OUT) :: npgrp
5658 CALL mvopen(lup,
'millepede.res')
5659 WRITE(lup,*)
'*** Results of checking input only, no solution performed ***'
5660 WRITE(lup,*)
'! === global parameters ==='
5661 WRITE(lup,*)
'! fixed-1: by pre-sigma, -2: by entries cut, -3: by iterated entries cut'
5663 WRITE(lup,*)
'! Label Value Pre-sigma SkippedEntries Cons. group Status '
5665 WRITE(lup,*)
'! Label Value Pre-sigma Entries Cons. group Status '
5681 IF (ivgbi <= 0)
THEN
5683 IF (ivgbi == -4)
THEN
5684 WRITE(lup,116) c1,itgbl,par,presig,icount,icgrp
5686 WRITE(lup,110) c1,itgbl,par,presig,icount,icgrp,ivgbi
5690 WRITE(lup,111) c1,itgbl,par,presig,icount,icgrp
5696 WRITE(lup,*)
'!.Appearance statistics '
5697 WRITE(lup,*)
'!. Label First file and record Last file and record #files #paired-par'
5700 IF (itpgrp > 0)
THEN
5708 WRITE(lup,*)
'* === constraint groups ==='
5710 WRITE(lup,*)
'* Group #Cons. Entries First label Last label'
5712 WRITE(lup,*)
'* Group #Cons. Entries First label Last label Paired label range'
5714 CALL mpalloc(vecpairedpargroups,length,
'paired global parameter groups (I)')
5726 IF (
icheck > 1 .AND. label1 > 0)
THEN
5730 vecpairedpargroups(npair+1)=0
5734 jpgrp=vecpairedpargroups(j)
5738 IF (ifrst /= 0.AND.inext /= (ilast+nstep))
THEN
5741 WRITE(lup,114) label1, label2
5746 IF (ifrst == 0) ifrst=inext
5753 IF (jpgrp == vecpairedpargroups(j+1)-1)
THEN
5758 IF (ifrst /= 0)
THEN
5761 WRITE(lup,114) label1, label2
5767 WRITE(lup,*)
'*.Appearance statistics '
5768 WRITE(lup,*)
'*. Group First file and record Last file and record #files'
5778110
FORMAT(
' !',a1,i10,2x,2g14.5,2i12,
' fixed',i2)
5779111
FORMAT(
' !',a1,i10,2x,2g14.5,2i12,
' variable')
5780112
FORMAT(
' !.',i10,6i11)
5781113
FORMAT(
' * ',i6,i8,3i12)
5782114
FORMAT(
' *:',48x,i12,
' ..',i12)
5783115
FORMAT(
' *.',i10,5i11)
5784116
FORMAT(
' !',a1,i10,2x,2g14.5,2i12,
' redundant')
5814 INTEGER(mpi) :: iproc
5824 INTEGER(mpi),
INTENT(IN) :: n
5825 INTEGER(mpl),
INTENT(IN) :: l
5826 REAL(mpd),
INTENT(IN) :: x(n)
5827 INTEGER(mpi),
INTENT(IN) :: is
5828 INTEGER(mpi),
INTENT(IN) :: ie
5829 REAL(mpd),
INTENT(OUT) :: b(n)
5834 INTEGER(mpl) :: indij
5835 INTEGER(mpl) :: indid
5837 INTEGER(mpi) :: ichunk
5873 CALL peend(24,
'Aborted, vector/matrix size mismatch')
5874 stop
'AVPRDS: mismatched vector and matrix'
5894 IF (ia2 <= ib2) b(ia2:ib2)=b(ia2:ib2)+
globalmatd(ia2:ib2)*x(ia2:ib2)
5905 IF (i <= ie.AND.i >= is)
THEN
5913 IF (j <= ie.AND.j >= is)
THEN
5929 IF (ja2 <= jb2)
THEN
5932 b(i)=b(i)+dot_product(
globalmatd(indij+lj+ja2-ja:indij+lj+jb2-ja),x(ja2:jb2))
5936 IF (
mextnd == 0.AND.ia2 <= ib2)
THEN
5939 b(j)=b(j)+dot_product(
globalmatd(indij+lj+jn*(ia2-ia):indij+lj+jn*(ib2-ia):jn),x(ia2:ib2))
5958 IF (i <= ie.AND.i >= is)
THEN
5966 IF (j <= ie.AND.j >= is)
THEN
5982 IF (ja2 <= jb2)
THEN
5985 b(i)=b(i)+dot_product(real(
globalmatf(indij+lj+ja2-ja:indij+lj+jb2-ja),mpd),x(ja2:jb2))
5989 IF (
mextnd == 0.AND.ia2 <= ib2)
THEN
5992 b(j)=b(j)+dot_product(real(
globalmatf(indij+lj+jn*(ia2-ia):indij+lj+jn*(ib2-ia):jn),mpd),x(ia2:ib2))
6026 INTEGER(mpi) :: iproc
6034 INTEGER(mpi),
INTENT(IN) :: n
6035 INTEGER(mpl),
INTENT(IN) :: l
6036 REAL(mpd),
INTENT(IN) :: x(n)
6037 REAL(mpd),
INTENT(OUT) :: b(n)
6042 INTEGER(mpl) :: indij
6043 INTEGER(mpl) :: indid
6045 INTEGER(mpi) :: ichunk
6074 CALL peend(24,
'Aborted, vector/matrix size mismatch')
6075 stop
'AVPRD0: mismatched vector and matrix'
6121 b(i)=b(i)+dot_product(
globalmatd(indij+lj:indij+lj+jn-1),x(ja:jb))
6127 b(j)=b(j)+dot_product(
globalmatd(indij+lj:indij+jn*in:jn),x(ia:ib))
6163 b(i)=b(i)+dot_product(real(
globalmatf(indij+lj:indij+lj+jn-1),mpd),x(ja:jb))
6169 b(j)=b(j)+dot_product(real(
globalmatf(indij+lj:indij+jn*in:jn),mpd),x(ia:ib))
6193 INTEGER(mpi) :: ispc
6204 INTEGER(mpl) :: indid
6205 INTEGER(mpl),
DIMENSION(12) :: icount
6212 icount(4)=huge(icount(4))
6213 icount(7)=huge(icount(7))
6214 icount(10)=huge(icount(10))
6224 ir=ipg+(ispc-1)*(
napgrp+1)
6231 icount(1)=icount(1)+in
6232 icount(2)=icount(2)+in*ku
6238 icount(3)=icount(3)+1
6239 icount(4)=min(icount(4),jn)
6240 icount(5)=icount(5)+jn
6241 icount(6)=max(icount(6),jn)
6242 icount(7)=min(icount(7),in)
6243 icount(8)=icount(8)+in
6244 icount(9)=max(icount(9),in)
6245 icount(10)=min(icount(10),in*jn)
6246 icount(11)=icount(11)+in*jn
6247 icount(12)=max(icount(12),in*jn)
6253 WRITE(*,*)
"analysis of sparsity structure"
6254 IF (icount(1) > 0)
THEN
6255 WRITE(*,101)
"rows without compression/blocks ", icount(1)
6256 WRITE(*,101)
" contained elements ", icount(2)
6258 WRITE(*,101)
"number of block matrices ", icount(3)
6259 avg=real(icount(5),mps)/real(icount(3),mps)
6260 WRITE(*,101)
"number of columns (min,mean,max) ", icount(4), avg, icount(6)
6261 avg=real(icount(8),mps)/real(icount(3),mps)
6262 WRITE(*,101)
"number of rows (min,mean,max) ", icount(7), avg, icount(9)
6263 avg=real(icount(11),mps)/real(icount(3),mps)
6264 WRITE(*,101)
"number of elements (min,mean,max) ", icount(10), avg, icount(12)
6265101
FORMAT(2x,a34,i10,f10.3,i10)
6284 INTEGER(mpi),
INTENT(IN) :: n
6285 REAL(mpd),
INTENT(IN) :: x(n)
6286 REAL(mpd),
INTENT(OUT) :: b(n)
6291 CALL peend(24,
'Aborted, vector/matrix size mismatch')
6292 stop
'AVPROD: mismatched vector and matrix'
6323 INTEGER(mpi) :: ispc
6324 INTEGER(mpi) :: item1
6325 INTEGER(mpi) :: item2
6326 INTEGER(mpi) :: itemc
6327 INTEGER(mpi) :: jtem
6328 INTEGER(mpi) :: jtemn
6331 INTEGER(mpi),
INTENT(IN) :: itema
6332 INTEGER(mpi),
INTENT(IN) :: itemb
6333 INTEGER(mpl),
INTENT(OUT) :: ij
6334 INTEGER(mpi),
INTENT(OUT) :: lr
6335 INTEGER(mpi),
INTENT(OUT) :: iprc
6346 item1=max(itema,itemb)
6347 item2=min(itema,itemb)
6348 IF(item2 <= 0.OR.item1 >
napgrp)
RETURN
6351 outer:
DO ispc=1,
nspc
6362 IF(ku < kl) cycle outer
6367 IF(item2 >= jtem.AND.item2 < jtemn)
THEN
6373 IF(item2 < jtem)
THEN
6375 ELSE IF(item2 >= jtemn)
THEN
6390 IF(ku < kl) cycle outer
6394 IF(itemc == jtem)
EXIT
6395 IF(itemc < jtem)
THEN
6397 ELSE IF(itemc > jtem)
THEN
6425 INTEGER(mpi),
INTENT(IN) :: itema
6426 INTEGER(mpi),
INTENT(IN) :: itemb
6451 INTEGER(mpi) :: item1
6452 INTEGER(mpi) :: item2
6453 INTEGER(mpi) :: ipg1
6454 INTEGER(mpi) :: ipg2
6456 INTEGER(mpi) :: iprc
6458 INTEGER(mpi),
INTENT(IN) :: itema
6459 INTEGER(mpi),
INTENT(IN) :: itemb
6461 INTEGER(mpl) ::
ijadd
6465 item1=max(itema,itemb)
6466 item2=min(itema,itemb)
6468 IF(item2 <= 0.OR.item1 >
nagb)
RETURN
6469 IF(item1 == item2)
THEN
6478 CALL ijpgrp(ipg1,ipg2,ij,lr,iprc)
6500 INTEGER(mpi) :: item1
6501 INTEGER(mpi) :: item2
6502 INTEGER(mpi) :: jtem
6504 INTEGER(mpi),
INTENT(IN) :: itema
6505 INTEGER(mpi),
INTENT(IN) :: itemb
6514 item1=max(itema,itemb)
6515 item2=min(itema,itemb)
6517 IF(item2 <= 0.OR.item1 >
nagb)
RETURN
6525 print *,
' IJCSR3 empty list ', item1, item2, ks, ke
6526 CALL peend(23,
'Aborted, bad matrix index')
6527 stop
'ijcsr3: empty list'
6532 IF(item1 == jtem)
EXIT
6533 IF(item1 < jtem)
THEN
6540 print *,
' IJCSR3 not found ', item1, item2, ks, ke
6541 CALL peend(23,
'Aborted, bad matrix index')
6542 stop
'ijcsr3: not found'
6558 INTEGER(mpi) :: item1
6559 INTEGER(mpi) :: item2
6563 INTEGER(mpl) ::
ijadd
6566 INTEGER(mpi),
INTENT(IN) :: itema
6567 INTEGER(mpi),
INTENT(IN) :: itemb
6572 item1=max(itema,itemb)
6573 item2=min(itema,itemb)
6574 IF(item2 <= 0.OR.item1 >
nagb)
RETURN
6583 ij=
ijadd(item1,item2)
6586 ELSE IF (ij < 0)
THEN
6618 INTEGER(mpi) :: ichunk
6622 INTEGER(mpi) :: ispc
6630 INTEGER(mpl) :: ijadd
6652 ir=ipg+(ispc-1)*(
napgrp+1)
6702 REAL(mps),
INTENT(IN) :: deltat
6703 INTEGER(mpi),
INTENT(OUT) :: minut
6704 INTEGER(mpi),
INTENT(OUT):: nhour
6705 REAL(mps),
INTENT(OUT):: secnd
6706 INTEGER(mpi) :: nsecd
6709 nsecd=nint(deltat,
mpi)
6711 minut=nsecd/60-60*nhour
6712 secnd=deltat-60*(minut+60*nhour)
6748 INTEGER(mpi),
INTENT(IN) :: item
6752 INTEGER(mpl) :: length
6753 INTEGER(mpl),
PARAMETER :: four = 4
6757 IF(item <= 0)
RETURN
6778 IF(j == 0)
EXIT inner
6800 IF (
lvllog > 1)
WRITE(
lunlog,*)
'INONE: array increased to', &
6821 INTEGER(mpi) :: iprime
6822 INTEGER(mpi) :: nused
6823 LOGICAL :: finalUpdate
6824 INTEGER(mpl) :: oldLength
6825 INTEGER(mpl) :: newLength
6826 INTEGER(mpl),
PARAMETER :: four = 4
6827 INTEGER(mpi),
DIMENSION(:,:),
ALLOCATABLE :: tempArr
6828 INTEGER(mpl),
DIMENSION(:),
ALLOCATABLE :: tempVec
6832 IF(finalupdate)
THEN
6841 CALL mpalloc(temparr,four,oldlength,
'INONE: temp array')
6843 CALL mpalloc(tempvec,oldlength,
'INONE: temp vector')
6867 IF(j == 0)
EXIT inner
6868 IF(j == i) cycle outer
6872 IF(.NOT.finalupdate)
RETURN
6877 WRITE(
lunlog,*)
'INONE: array reduced to',newlength,
' words'
6901 IF(j == 0)
EXIT inner
6902 IF(j == i) cycle outer
6919 INTEGER(mpi),
INTENT(IN) :: n
6920 INTEGER(mpi) :: nprime
6921 INTEGER(mpi) :: nsqrt
6926 IF(mod(nprime,2) == 0) nprime=nprime+1
6929 nsqrt=int(sqrt(real(nprime,
mps)),
mpi)
6931 IF(i*(nprime/i) == nprime) cycle outer
6953 INTEGER(mpi) :: idum
6955 INTEGER(mpi) :: indab
6956 INTEGER(mpi) :: itgbi
6957 INTEGER(mpi) :: itgbl
6958 INTEGER(mpi) :: ivgbi
6960 INTEGER(mpi) :: jgrp
6961 INTEGER(mpi) :: lgrp
6963 INTEGER(mpi) :: nc31
6965 INTEGER(mpi) :: nwrd
6966 INTEGER(mpi) :: inone
6971 INTEGER(mpl) :: length
6972 INTEGER(mpl) :: rows
6976 WRITE(
lunlog,*)
'LOOP1: starting'
6999 WRITE(
lunlog,*)
'LOOP1: reading data files'
7010 CALL hmpldf(1,
'Number of words/record in binary file')
7011 CALL hmpdef(8,0.0,60.0,
'not_stored data per record')
7016 CALL peend(20,
'Aborted, bad binary records')
7017 stop
'LOOP1: length of binary record exceeds cache size, wrong file type?'
7031 WRITE(*,*)
'Read all binary data files:'
7050 WRITE(
lunlog,*)
'LOOP1: reading data files again'
7062 CALL peend(21,
'Aborted, no labels/parameters defined')
7063 stop
'LOOP1: no labels/parameters defined'
7068 ' is total number NTGB of labels/parameters'
7070 CALL hmpldf(2,
'Number of entries per label')
7114 WRITE(
lunlog,*)
'LOOP1:',
npresg,
' is number of pre-sigmas'
7115 WRITE(*,*)
'LOOP1:',
npresg,
' is number of pre-sigmas'
7116 IF(
npresg == 0)
WRITE(*,*)
'Warning: no pre-sigmas defined'
7138 WRITE(
lunlog,*)
'LOOP1:',
nvgb,
' is number NVGB of variable parameters'
7143 WRITE(
lunlog,*)
'LOOP1: counting records, NO iteration of entries cut !'
7149 CALL hmpdef(15,0.0,120.0,
'Number of parameters per group')
7184 ' is total number NTPGRP of label/parameter groups'
7200 WRITE(*,112)
' Default pre-sigma =',
regpre, &
7201 ' (if no individual pre-sigma defined)'
7202 WRITE(*,*)
'Pre-sigma factor is',
regula
7205 WRITE(*,*)
'No regularization will be done'
7207 WRITE(*,*)
'Regularization will be done, using factor',
regula
7211 CALL peend(22,
'Aborted, no variable global parameters')
7212 stop
'... no variable global parameters'
7219 IF(presg > 0.0)
THEN
7221 ELSE IF(presg == 0.0.AND.
regpre > 0.0)
THEN
7222 prewt=1.0/real(
regpre**2,mpd)
7238 WRITE(*,101)
'NTGB',
ntgb,
'total number of parameters'
7239 WRITE(*,101)
'NVGB',
nvgb,
'number of variable parameters'
7243 WRITE(*,101)
'Too many variable parameters for diagonalization, fallback is inversion'
7251 WRITE(*,101)
' NREC',
nrec,
'number of records'
7252 IF (
nrecd > 0)
WRITE(*,101)
' NRECD',
nrec,
'number of records containing doubles'
7253 WRITE(*,101)
' NEQN',
neqn,
'number of equations (measurements)'
7254 WRITE(*,101)
' NEGB',
negb,
'number of equations with global parameters'
7255 WRITE(*,101)
' NDGB',
ndgb,
'number of global derivatives'
7257 WRITE(*,101)
' NZGB',
nzgb,
'number of zero global der. (ignored in entry counts)'
7260 WRITE(*,101)
'MREQENF',
mreqenf,
'required number of entries (eqns in binary files)'
7262 WRITE(*,101)
'MREQENF',
mreqenf,
'required number of entries (recs in binary files)'
7265 WRITE(*,101)
'ITEREN',
iteren,
'iterate cut for parameters with less entries'
7266 WRITE(*,101)
'MREQENA',
mreqena,
'required number of entries (from accepted fits)'
7267 IF (
mreqpe > 1)
WRITE(*,101) &
7268 'MREQPE',
mreqpe,
'required number of pair entries'
7269 IF (
msngpe >= 1)
WRITE(*,101) &
7270 'MSNGPE',
msngpe,
'max pair entries single prec. storage'
7271 WRITE(*,101)
'NTGB',
ntgb,
'total number of parameters'
7272 WRITE(*,101)
'NVGB',
nvgb,
'number of variable parameters'
7275 WRITE(*,*)
'Global parameter labels:'
7282 mqi=((mqi-20)/20)*20+1
7290 WRITE(8,101)
' NREC',
nrec,
'number of records'
7291 IF (
nrecd > 0)
WRITE(8,101)
' NRECD',
nrec,
'number of records containing doubles'
7292 WRITE(8,101)
' NEQN',
neqn,
'number of equations (measurements)'
7293 WRITE(8,101)
' NEGB',
negb,
'number of equations with global parameters'
7294 WRITE(8,101)
' NDGB',
ndgb,
'number of global derivatives'
7296 WRITE(8,101)
'MREQENF',
mreqenf,
'required number of entries (eqns in binary files)'
7298 WRITE(8,101)
'MREQENF',
mreqenf,
'required number of entries (recs in binary files)'
7301 WRITE(8,101)
'ITEREN',
iteren,
'iterate cut for parameters with less entries'
7302 WRITE(8,101)
'MREQENA',
mreqena,
'required number of entries (from accepted fits)'
7304 WRITE(
lunlog,*)
'LOOP1: ending'
7308101
FORMAT(1x,a8,
' =',i14,
' = ',a)
7324 INTEGER(mpi) :: ibuf
7326 INTEGER(mpi) :: indab
7332 INTEGER(mpi) :: nc31
7334 INTEGER(mpi) :: nlow
7336 INTEGER(mpi) :: nwrd
7338 INTEGER(mpl) :: length
7339 INTEGER(mpl),
DIMENSION(:),
ALLOCATABLE :: newCounter
7344 WRITE(
lunlog,*)
'LOOP1: iterating'
7346 WRITE(*,*)
'LOOP1: iterating'
7349 CALL mpalloc(newcounter,length,
'new entries counter')
7373 CALL isjajb(nst,ist,ja,jb,jsp)
7374 IF(ja == 0.AND.jb == 0)
EXIT
7380 IF(ij == -2) nlow=nlow+1
7385 newcounter(ij)=newcounter(ij)+1
7414 WRITE(
lunlog,*)
'LOOP1:',
nvgb,
' is number NVGB of variable parameters'
7444 INTEGER(mpi) :: ibuf
7445 INTEGER(mpi) :: icblst
7446 INTEGER(mpi) :: icboff
7447 INTEGER(mpi) :: icgb
7448 INTEGER(mpi) :: icgrp
7449 INTEGER(mpi) :: icount
7450 INTEGER(mpi) :: iext
7451 INTEGER(mpi) :: ihis
7455 INTEGER(mpi) :: ioff
7456 INTEGER(mpi) :: ipoff
7457 INTEGER(mpi) :: iproc
7458 INTEGER(mpi) :: irecmm
7460 INTEGER(mpi) :: itgbi
7461 INTEGER(mpi) :: itgbij
7462 INTEGER(mpi) :: itgbik
7463 INTEGER(mpi) :: ivgbij
7464 INTEGER(mpi) :: ivgbik
7465 INTEGER(mpi) :: ivpgrp
7469 INTEGER(mpi) :: jcgrp
7470 INTEGER(mpi) :: jext
7471 INTEGER(mpi) :: jcgb
7472 INTEGER(mpi) :: jrec
7474 INTEGER(mpi) :: joff
7476 INTEGER(mpi) :: kcgrp
7477 INTEGER(mpi) :: kfile
7479 INTEGER(mpi) :: label
7480 INTEGER(mpi) :: labelf
7481 INTEGER(mpi) :: labell
7482 INTEGER(mpi) :: lvpgrp
7485 INTEGER(mpi) :: maeqnf
7486 INTEGER(mpi) :: nall
7487 INTEGER(mpi) :: naeqna
7488 INTEGER(mpi) :: naeqnf
7489 INTEGER(mpi) :: naeqng
7490 INTEGER(mpi) :: npdblk
7491 INTEGER(mpi) :: nc31
7492 INTEGER(mpi) :: ncachd
7493 INTEGER(mpi) :: ncachi
7494 INTEGER(mpi) :: ncachr
7495 INTEGER(mpi) :: ncon
7498 INTEGER(mpi) :: ndfmax
7499 INTEGER(mpi) :: nfixed
7500 INTEGER(mpi) :: nggd
7501 INTEGER(mpi) :: nggi
7502 INTEGER(mpi) :: nmatmo
7503 INTEGER(mpi) :: noff
7504 INTEGER(mpi) :: npair
7505 INTEGER(mpi) :: npar
7506 INTEGER(mpi) :: nparmx
7508 INTEGER(mpi) :: nrece
7509 INTEGER(mpi) :: nrecf
7510 INTEGER(mpi) :: nrecmm
7512 INTEGER(mpi) :: nwrd
7513 INTEGER(mpi) :: inone
7522 INTEGER(mpl):: nblock
7523 INTEGER(mpl):: nbwrds
7524 INTEGER(mpl):: noff8
7525 INTEGER(mpl):: ndimbi
7526 INTEGER(mpl):: ndimsa(4)
7528 INTEGER(mpl):: nnzero
7529 INTEGER(mpl):: matsiz(2)
7530 INTEGER(mpl):: matwords
7531 INTEGER(mpl):: mbwrds
7532 INTEGER(mpl):: length
7535 INTEGER(mpl),
PARAMETER :: two=2
7536 INTEGER(mpi) :: maxGlobalPar = 0
7537 INTEGER(mpi) :: maxLocalPar = 0
7538 INTEGER(mpi) :: maxEquations = 0
7540 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: vecConsGroupList
7541 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: vecConsGroupIndex
7542 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: vecPairedParGroups
7543 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: vecBlockCounts
7546 SUBROUTINE ndbits(npgrp,ndims,nsparr,ihst)
7548 INTEGER(mpi),
DIMENSION(:),
INTENT(IN) :: npgrp
7549 INTEGER(mpl),
DIMENSION(4),
INTENT(OUT) :: ndims
7550 INTEGER(mpl),
DIMENSION(:,:),
INTENT(OUT) :: nsparr
7551 INTEGER(mpi),
INTENT(IN) :: ihst
7553 SUBROUTINE ckbits(npgrp,ndims)
7555 INTEGER(mpi),
DIMENSION(:),
INTENT(IN) :: npgrp
7556 INTEGER(mpl),
DIMENSION(4),
INTENT(OUT) :: ndims
7558 SUBROUTINE spbits(npgrp,nsparr,nsparc)
7560 INTEGER(mpi),
DIMENSION(:),
INTENT(IN) :: npgrp
7561 INTEGER(mpl),
DIMENSION(:,:),
INTENT(IN) :: nsparr
7562 INTEGER(mpi),
DIMENSION(:),
INTENT(OUT) :: nsparc
7564 SUBROUTINE gpbmap(ngroup,npgrp,npair)
7566 INTEGER(mpi),
INTENT(IN) :: ngroup
7567 INTEGER(mpi),
DIMENSION(:,:),
INTENT(IN) :: npgrp
7568 INTEGER(mpi),
DIMENSION(:),
INTENT(OUT) :: npair
7570 SUBROUTINE ggbmap(ipgrp,npair,npgrp)
7572 INTEGER(mpi),
INTENT(IN) :: ipgrp
7573 INTEGER(mpi),
INTENT(OUT) :: npair
7574 INTEGER(mpi),
DIMENSION(:),
INTENT(OUT) :: npgrp
7576 SUBROUTINE pbsbits(npgrp,ibsize,nnzero,nblock,nbkrow)
7578 INTEGER(mpi),
DIMENSION(:),
INTENT(IN) :: npgrp
7579 INTEGER(mpi),
INTENT(IN) :: ibsize
7580 INTEGER(mpl),
INTENT(OUT) :: nnzero
7581 INTEGER(mpl),
INTENT(OUT) :: nblock
7582 INTEGER(mpi),
DIMENSION(:),
INTENT(OUT) :: nbkrow
7584 SUBROUTINE pblbits(npgrp,ibsize,nsparr,nsparc)
7586 INTEGER(mpi),
DIMENSION(:),
INTENT(IN) :: npgrp
7587 INTEGER(mpi),
INTENT(IN) :: ibsize
7588 INTEGER(mpl),
DIMENSION(:),
INTENT(IN) :: nsparr
7589 INTEGER(mpl),
DIMENSION(:),
INTENT(OUT) :: nsparc
7591 SUBROUTINE prbits(npgrp,nsparr)
7593 INTEGER(mpi),
DIMENSION(:),
INTENT(IN) :: npgrp
7594 INTEGER(mpl),
DIMENSION(:),
INTENT(OUT) :: nsparr
7596 SUBROUTINE pcbits(npgrp,nsparr,nsparc)
7598 INTEGER(mpi),
DIMENSION(:),
INTENT(IN) :: npgrp
7599 INTEGER(mpl),
DIMENSION(:),
INTENT(IN) :: nsparr
7600 INTEGER(mpl),
DIMENSION(:),
INTENT(OUT) :: nsparc
7610 WRITE(
lunlog,*)
'LOOP2: starting'
7633 WRITE(
lunlog,*)
' Elimination for constraints enforced for solution by decomposition!'
7638 WRITE(
lunlog,*)
' Lagrange multipliers enforced for solution by sparsePARDISO!'
7643 WRITE(
lunlog,*)
' Elimination for constraints with mpqldec enforced (LAPACK only for unpacked storage)!'
7660 noff8=int(
nagb,mpl)*int(
nagb-1,mpl)/2
7696 print *,
' Variable parameter groups ',
nvpgrp
7744 CALL mpalloc(vecconsgrouplist,length,
'constraint group list')
7745 CALL mpalloc(vecconsgroupindex,length,
'constraint group index')
7755 WRITE(
lunlog,*)
'LOOP2: start event reading'
7795 WRITE(*,*)
'Record number ',
nrec,
' from file ',kfile
7796 IF (wgh /= 1.0)
WRITE(*,*)
' weight ',wrec
7800 CALL isjajb(nst,ist,ja,jb,jsp)
7804 IF(nda ==
mdebg2+1)
WRITE(*,*)
'... and more data'
7809 WRITE(*,*)
'Local derivatives:'
7811107
FORMAT(6(i3,g12.4))
7813 WRITE(*,*)
'Global derivatives:'
7816108
FORMAT(3i11,g12.4)
7819 WRITE(*,*)
'total_par_label __label__ var_par_index derivative'
7834 CALL isjajb(nst,ist,ja,jb,jsp)
7835 IF(ja == 0.AND.jb == 0)
EXIT
7875 IF (kcgrp == jcgrp) cycle
7886 IF (vecconsgroupindex(k) == 0)
THEN
7889 vecconsgrouplist(icgrp)=k
7904 IF (vecconsgroupindex(k) < icount)
THEN
7906 vecconsgroupindex(k)=icount
7924 IF (nfixed > 0) naeqnf=naeqnf+1
7927 IF(ja /= 0.AND.jb /= 0)
THEN
7936 IF (naeqnf > maeqnf) nrecf=nrecf+1
7940 maxglobalpar=max(
nagbn,maxglobalpar)
7941 maxlocalpar=max(
nalcn,maxlocalpar)
7942 maxequations=max(
naeqn,maxequations)
7945 dstat(1)=dstat(1)+real((nwrd+2)*2,mpd)
7946 dstat(2)=dstat(2)+real(
nagbn+2,mpd)
7951 vecconsgroupindex(vecconsgrouplist(k))=0
7956 IF (
nagbn == 0)
THEN
7975 IF (ivpgrp /= lvpgrp)
THEN
8030 (irecmm >= nrecmm.OR.irecmm ==
mxrec))
THEN
8031 IF (nmatmo == 0)
THEN
8033 WRITE(*,*)
'Monitoring of sparse matrix construction'
8034 WRITE(*,*)
' records ........ off-diagonal elements ', &
8035 '....... compression memory'
8036 WRITE(*,*)
' non-zero used(double) used', &
8041 gbc=1.0e-9*real((mpi*ndimsa(2)+mpd*ndimsa(3)+mps*ndimsa(4))/mpi*(bit_size(1_mpi)/8),mps)
8042 gbu=1.0e-9*real(((mpi+mpd)*(ndimsa(3)+ndimsa(4)))/mpi*(bit_size(1_mpi)/8),mps)
8044 WRITE(*,1177) irecmm,ndimsa(1),ndimsa(3),ndimsa(4),cpr,gbc
80451177
FORMAT(i9,3i13,f10.2,f11.6)
8046 DO WHILE(irecmm >= nrecmm)
8065 WRITE(
lunlog,*)
'LOOP2: event reading ended - end of data'
8067 dstat(k)=dstat(k)/real(
nrec,mpd)
8081 CALL mpalloc(vecpairedpargroups,length,
'paired global parameter groups (I)')
8083 print *,
' Total parameter groups pairs',
ntpgrp
8086 CALL ggbmap(i,npair,vecpairedpargroups)
8154 IF(ivgbij > 0.AND.ivgbik > 0)
THEN
8156 IF (
mprint > 1)
WRITE(*,*)
'add index pair ',ivgbij,ivgbik
8233 nparmx=max(nparmx,int(rows,mpi))
8255 WRITE(*,*)
' Parameter block', i, ib-ia, jb-ja, labelf, labell
8259 WRITE(
lunlog,*)
'Detected',
npblck,
'(disjoint) parameter blocks, max size ', nparmx
8261 WRITE(*,*)
'Detected',
npblck,
'(disjoint) parameter blocks, max size ', nparmx
8263 WRITE(*,*)
'Using block diagonal storage mode'
8278 IF (
nagb >= 65536)
THEN
8279 noff=int(noff8/1000,mpi)
8289 CALL hmpdef(ihis,0.0,real(
mhispe,mps),
'NDBITS: #off-diagonal elements')
8294 ndgn=ndimsa(3)+ndimsa(4)
8295 matwords=ndimsa(2)+length*4
8310 length=int(npdblk,mpl)
8311 CALL mpalloc(vecblockcounts,length,
'sparse matrix row offsets (CSR3)')
8313 nbwrds=2*int(nblock,mpl)*int(
ipdbsz(i)*
ipdbsz(i)+1,mpl)
8314 IF ((i == 1).OR.(nbwrds < mbwrds))
THEN
8338 IF (
fcache(2) == 0.0)
THEN
8340 fcache(2)=real(dstat(2),mps)
8341 fcache(3)=real(dstat(3),mps)
8362 ncachd=
ncache-ncachr-ncachi
8393 WRITE(lu,101)
'NTGB',
ntgb,
'total number of parameters'
8394 WRITE(lu,102)
'(all parameters, appearing in binary files)'
8395 WRITE(lu,101)
'NVGB',
nvgb,
'number of variable parameters'
8396 WRITE(lu,102)
'(appearing in fit matrix/vectors)'
8397 WRITE(lu,101)
'NAGB',
nagb,
'number of all parameters'
8398 WRITE(lu,102)
'(including Lagrange multiplier or reduced)'
8399 WRITE(lu,101)
'NTPGRP',
ntpgrp,
'total number of parameter groups'
8400 WRITE(lu,101)
'NVPGRP',
nvpgrp,
'number of variable parameter groups'
8401 WRITE(lu,101)
'NFGB',
nfgb,
'number of fit parameters'
8403 WRITE(lu,101)
'MBANDW',
mbandw,
'band width of preconditioner matrix'
8404 WRITE(lu,102)
'(if <0, no preconditioner matrix)'
8406 IF (
nagb >= 65536)
THEN
8407 WRITE(lu,101)
'NOFF/K',noff,
'max number of off-diagonal elements'
8409 WRITE(lu,101)
'NOFF',noff,
'max number of off-diagonal elements'
8412 IF (
nagb >= 65536)
THEN
8413 WRITE(lu,101)
'NDGN/K',ndgn/1000,
'actual number of off-diagonal elements'
8415 WRITE(lu,101)
'NDGN',ndgn,
'actual number of off-diagonal elements'
8418 WRITE(lu,101)
'NCGB',
ncgb,
'number of constraints'
8419 WRITE(lu,101)
'NAGBN',
nagbn,
'max number of global parameters in an event'
8420 WRITE(lu,101)
'NALCN',
nalcn,
'max number of local parameters in an event'
8421 WRITE(lu,101)
'NAEQN',
naeqn,
'max number of equations in an event'
8423 WRITE(lu,101)
'NAEQNA',naeqna,
'number of equations'
8424 WRITE(lu,101)
'NAEQNG',naeqng, &
8425 'number of equations with global parameters'
8426 WRITE(lu,101)
'NAEQNF',naeqnf, &
8427 'number of equations with fixed global parameters'
8428 WRITE(lu,101)
'NRECF',nrecf, &
8429 'number of records with fixed global parameters'
8432 WRITE(lu,101)
'NRECE',nrece, &
8433 'number of records without variable parameters'
8436 WRITE(lu,101)
'NCACHE',
ncache,
'number of words for caching'
8437 WRITE(lu,111) (
fcache(k)*100.0,k=1,3)
8438111
FORMAT(22x,
'cache splitting ',3(f6.1,
' %'))
8443 WRITE(lu,*)
'Solution method and matrix-storage mode:'
8445 WRITE(lu,*)
' METSOL = 1: matrix inversion'
8446 ELSE IF(
metsol == 2)
THEN
8447 WRITE(lu,*)
' METSOL = 2: diagonalization'
8448 ELSE IF(
metsol == 3)
THEN
8449 WRITE(lu,*)
' METSOL = 3: decomposition'
8450 ELSE IF(
metsol == 4)
THEN
8451 WRITE(lu,*)
' METSOL = 4: MINRES (rtol',
mrestl,
')'
8452 ELSE IF(
metsol == 5)
THEN
8453 WRITE(lu,*)
' METSOL = 5: MINRES-QLP (rtol',
mrestl,
')'
8454 ELSE IF(
metsol == 6)
THEN
8455 WRITE(lu,*)
' METSOL = 6: GMRES'
8457 ELSE IF(
metsol == 7)
THEN
8458 WRITE(lu,*)
' METSOL = 7: LAPACK factorization'
8459 ELSE IF(
metsol == 8)
THEN
8460 WRITE(lu,*)
' METSOL = 8: LAPACK factorization'
8462 ELSE IF(
metsol == 9)
THEN
8463 WRITE(lu,*)
' METSOL = 9: Intel oneMKL PARDISO'
8467 WRITE(lu,*)
' with',
mitera,
' iterations'
8469 WRITE(lu,*)
' MATSTO = 0: unpacked symmetric matrix, ',
'n*n elements'
8470 ELSE IF(
matsto == 1)
THEN
8471 WRITE(lu,*)
' MATSTO = 1: full symmetric matrix, ',
'(n*n+n)/2 elements'
8472 ELSE IF(
matsto == 2)
THEN
8473 WRITE(lu,*)
' MATSTO = 2: sparse matrix (custom)'
8474 ELSE IF(
matsto == 3)
THEN
8476 WRITE(lu,*)
' MATSTO = 3: sparse matrix (upper triangle, CSR3)'
8478 WRITE(lu,*)
' MATSTO = 3: sparse matrix (upper triangle, BSR3)'
8479 WRITE(lu,*)
' block size',
matbsz
8483 WRITE(lu,*)
' block diagonal with',
npblck,
' blocks'
8485 IF(
mextnd>0)
WRITE(lu,*)
' with extended storage'
8486 IF(
dflim /= 0.0)
THEN
8487 WRITE(lu,103)
'Convergence assumed, if expected dF <',
dflim
8492 WRITE(lu,*)
'Constraints handled by elimination with LAPACK'
8494 WRITE(lu,*)
'Constraints handled by elimination'
8497 WRITE(lu,*)
'Constraints handled by Lagrange multipliers'
8504 CALL peend(28,
'Aborted, no local parameters')
8505 stop
'LOOP2: stopping due to missing local parameters'
8526105
FORMAT(
' Constants C1, C2 for Wolfe conditions:',g12.4,
', ',g12.4)
8534 matwords=(length+
nagb+1)*2
8541 ELSE IF (
matsto == 2)
THEN
8542 matsiz(1)=ndimsa(3)+
nagb
8557 IF (icblst > icboff)
THEN
8563 nall = npar;
IF (
icelim <= 0) nall=npar+ncon
8567 matsiz(1)=matsiz(1)+k
8569 matsiz(1)=matsiz(1)+nall
8574 matwords=matwords+matsiz(1)*2+matsiz(2)
8580 IF (matwords < 250000)
THEN
8581 WRITE(*,*)
'Size of global matrix: < 1 MB'
8583 WRITE(*,*)
'Size of global matrix:',int(real(matwords,mps)*4.0e-6,mpi),
' MB'
8589 WRITE(
lunlog,*)
' Cut values of Chi^2/Ndf and Chi2,'
8590 WRITE(
lunlog,*)
' corresponding to 2 and 3 standard deviations'
8591 WRITE(
lunlog,*)
' Ndf Chi^2/Ndf(2) Chi^2(2) ', &
8592 ' Chi^2/Ndf(3) Chi^2(3)'
8595 IF(ndf >
naeqn)
EXIT
8598 ELSE IF(ndf < 20)
THEN
8600 ELSE IF(ndf < 100)
THEN
8602 ELSE IF(ndf < 200)
THEN
8609 WRITE(
lunlog,106) ndf,chin2,chin2*real(ndf,mps),chin3, chin3*real(ndf,mps)
8612 WRITE(
lunlog,*)
'LOOP2: ending'
8616 IF (
ncgbe /= 0)
THEN
8619 WRITE(*,199)
'WarningWarningWarningWarningWarningWarningWarningWarningWar'
8620 WRITE(*,199)
'arningWarningWarningWarningWarningWarningWarningWarningWarn'
8621 WRITE(*,199)
'rningWarningWarningWarningWarningWarningWarningWarningWarni'
8622 WRITE(*,199)
'ningWarningWarningWarningWarningWarningWarningWarningWarnin'
8623 WRITE(*,199)
'ingWarningWarningWarningWarningWarningWarningWarningWarning'
8624 WRITE(*,199)
'ngWarningWarningWarningWarningWarningWarningWarningWarningW'
8625 WRITE(*,199)
'gWarningWarningWarningWarningWarningWarningWarningWarningWa'
8627 WRITE(*,*)
' Number of empty constraints =',abs(
ncgbe),
', should be 0'
8628 WRITE(*,*)
' => please check constraint definition, mille data'
8630 WRITE(*,199)
'WarningWarningWarningWarningWarningWarningWarningWarningWar'
8631 WRITE(*,199)
'arningWarningWarningWarningWarningWarningWarningWarningWarn'
8632 WRITE(*,199)
'rningWarningWarningWarningWarningWarningWarningWarningWarni'
8633 WRITE(*,199)
'ningWarningWarningWarningWarningWarningWarningWarningWarnin'
8634 WRITE(*,199)
'ingWarningWarningWarningWarningWarningWarningWarningWarning'
8635 WRITE(*,199)
'ngWarningWarningWarningWarningWarningWarningWarningWarningW'
8636 WRITE(*,199)
'gWarningWarningWarningWarningWarningWarningWarningWarningWa'
8641101
FORMAT(1x,a8,
' =',i14,
' = ',a)
8643103
FORMAT(1x,a,g12.4)
8644106
FORMAT(i6,2(3x,f9.3,f12.1,3x))
8659 INTEGER(mpi) :: imed
8662 INTEGER(mpi) :: nent
8663 INTEGER(mpi),
DIMENSION(measBins) :: isuml
8664 INTEGER(mpi),
DIMENSION(measBins) :: isums
8668 INTEGER(mpl) :: ioff
8671 DATA lfirst /.true./
8684 WRITE(
lunmon,
'(A)')
'*** Normalized residuals grouped by first global label (per local fit cycle) ***'
8686 WRITE(
lunmon,
'(A)')
'*** Pulls grouped by first global label (per local fit cycle) ***'
8688 WRITE(
lunmon,
'(A)')
'! LFC Label Entries Median RMS(MAD) <error>'
8693#ifdef SCOREP_USER_ENABLE
8694 scorep_user_region_by_name_begin(
"UR_monres", scorep_user_region_type_common)
8710 IF (2*isuml(j) > nent)
EXIT
8714 IF (isuml(j) > isuml(j-1)) amed=amed+real(nent-2*isuml(j-1),mps)/real(2*isuml(j)-2*isuml(j-1),mps)
8727 isums(j)=isums(j)+isums(j-1)
8731 IF (2*isums(j) > nent)
EXIT
8734 IF (isums(j) > isums(j-1)) amad=amad+real(nent-2*isums(j-1),mps)/real(2*isums(j)-2*isums(j-1),mps)
8746#ifdef SCOREP_USER_ENABLE
8747 scorep_user_region_by_name_end(
"UR_monres")
8751110
FORMAT(i5,2i10,3g14.5)
8766 INTEGER(mpi) :: ioff
8767 INTEGER(mpi) :: ipar0
8768 INTEGER(mpi) :: ncon
8769 INTEGER(mpi) :: npar
8770 INTEGER(mpi) :: nextra
8772 INTEGER :: nbopt, nboptx, ILAENV
8775 INTEGER(mpl),
INTENT(IN) :: msize(2)
8777 INTEGER(mpl) :: length
8778 INTEGER(mpl) :: nwrdpc
8779 INTEGER(mpl),
PARAMETER :: three = 3
8792 CALL mpalloc(
aux,length,
' local fit scratch array: aux')
8793 CALL mpalloc(
vbnd,length,
' local fit scratch array: vbnd')
8794 CALL mpalloc(
vbdr,length,
' local fit scratch array: vbdr')
8797 CALL mpalloc(
vbk,length,
' local fit scratch array: vbk')
8800 CALL mpalloc(
vzru,length,
' local fit scratch array: vzru')
8801 CALL mpalloc(
scdiag,length,
' local fit scratch array: scdiag')
8802 CALL mpalloc(
scflag,length,
' local fit scratch array: scflag')
8803 CALL mpalloc(
ibandh,2*length,
' local fit band width hist.: ibandh')
8817 length=4*
ncblck;
IF(ncon == 0) length=0
8826 nwrdpc=nwrdpc+length
8838 length=(ncon*ncon+ncon)/2
8863 length=length+npar+nextra
8869 IF(
mbandw == 0) length=length+1
8877 nwrdpc=nwrdpc+2*length
8878 IF (nwrdpc > 250000)
THEN
8880 WRITE(*,*)
'Size of preconditioner matrix:',int(real(nwrdpc,mps)*4.0e-6,mpi),
' MB'
8902 CALL peend(23,
'Aborted, bad matrix index (will exceed 32bit)')
8903 stop
'vmprep: bad index (matrix to large for diagonalization)'
8925 nbopt = ilaenv( 1_mpl,
'DSYTRF',
'U', int(
nagb,mpl), int(
nagb,mpl), -1_mpl, -1_mpl )
8927 print *,
'LAPACK optimal block size for DSYTRF:', nbopt
8928 lplwrk=length*int(nbopt,mpl)
8936 nbopt = ilaenv( 1_mpl,
'DORMQL',
'RN', int(npar,mpl), int(npar,mpl), int(ncon,mpl), int(npar,mpl) )
8937 IF (int(npar,mpl)*int(nbopt,mpl) >
lplwrk)
THEN
8938 lplwrk=int(npar,mpl)*int(nbopt,mpl)
8943 print *,
'LAPACK optimal block size for DORMQL:', nboptx
8962 INTEGER(mpi) :: icoff
8963 INTEGER(mpi) :: ipoff
8966 INTEGER(mpi) :: ncon
8967 INTEGER(mpi) :: nfit
8968 INTEGER(mpi) :: npar
8969 INTEGER(mpi) :: nrank
8970 INTEGER(mpl) :: imoff
8971 INTEGER(mpl) :: ioff1
8989 WRITE(
lunlog,*)
'Shrinkage of global matrix (A->Q^t*A*Q)'
9004 nfit=npar+ncon;
IF (
icelim > 0) nfit=npar-ncon
9006 IF(nfit < npar)
THEN
9024 WRITE(
lunlog,*)
'Inversion of global matrix (A->A^-1)'
9031 IF(nfit /= nrank)
THEN
9032 WRITE(*,*)
'Warning: the rank defect of the symmetric',nfit, &
9033 '-by-',nfit,
' matrix is ',nfit-nrank,
' (should be zero).'
9034 WRITE(lun,*)
'Warning: the rank defect of the symmetric',nfit, &
9035 '-by-',nfit,
' matrix is ',nfit-nrank,
' (should be zero).'
9038 WRITE(*,*)
' --> enforcing SUBITO mode'
9039 WRITE(lun,*)
' --> enforcing SUBITO mode'
9041 ELSE IF(
ndefec == 0)
THEN
9043 WRITE(lun,*)
'No rank defect of the symmetric matrix'
9045 WRITE(lun,*)
'No rank defect of the symmetric block', ib,
' of size', npar
9056 IF(nfit < npar)
THEN
9075 INTEGER(mpi) :: icoff
9076 INTEGER(mpi) :: ipoff
9079 INTEGER(mpi) :: ncon
9080 INTEGER(mpi) :: nfit
9081 INTEGER(mpi) :: npar
9082 INTEGER(mpi) :: nrank
9083 INTEGER(mpl) :: imoff
9084 INTEGER(mpl) :: ioff1
9099 WRITE(
lunlog,*)
'Shrinkage of global matrix (A->Q^t*A*Q)'
9113 nfit=npar+ncon;
IF (
icelim > 0) nfit=npar-ncon
9115 IF(nfit < npar)
THEN
9133 WRITE(
lunlog,*)
'Decomposition of global matrix (A->L*D*L^t)'
9139 IF(nfit /= nrank)
THEN
9140 WRITE(*,*)
'Warning: the rank defect of the symmetric',nfit, &
9141 '-by-',nfit,
' matrix is ',nfit-nrank,
' (should be zero).'
9142 WRITE(lun,*)
'Warning: the rank defect of the symmetric',nfit, &
9143 '-by-',nfit,
' matrix is ',nfit-nrank,
' (should be zero).'
9146 WRITE(*,*)
' --> enforcing SUBITO mode'
9147 WRITE(lun,*)
' --> enforcing SUBITO mode'
9149 ELSE IF(
ndefec == 0)
THEN
9151 WRITE(lun,*)
'No rank defect of the symmetric matrix'
9153 WRITE(lun,*)
'No rank defect of the symmetric block', ib,
' of size', npar
9155 WRITE(lun,*)
' largest diagonal element (LDLt)', evmax
9156 WRITE(lun,*)
' smallest diagonal element (LDLt)', evmin
9165 IF(nfit < npar)
THEN
9187 INTEGER(mpi) :: icoff
9188 INTEGER(mpi) :: ipoff
9191 INTEGER(mpi) :: ncon
9192 INTEGER(mpi) :: nfit
9193 INTEGER(mpi) :: npar
9194 INTEGER(mpl) :: imoff
9195 INTEGER(mpl) :: ioff1
9196 INTEGER(mpi) :: infolp
9216 WRITE(
lunlog,*)
'Shrinkage of global matrix (A->Q^t*A*Q)'
9231 nfit=npar+ncon;
IF (
icelim > 0) nfit=npar-ncon
9233 IF(nfit < npar)
THEN
9250 IF (nfit > npar)
THEN
9253 WRITE(
lunlog,*)
'Factorization of global matrix (A->L*D*L^t)'
9257#ifdef SCOREP_USER_ENABLE
9258 scorep_user_region_by_name_begin(
"UR_dsptrf", scorep_user_region_type_common)
9261#ifdef SCOREP_USER_ENABLE
9262 scorep_user_region_by_name_end(
"UR_dsptrf")
9269 WRITE(
lunlog,*)
'Factorization of global matrix (A->L*L^t)'
9273#ifdef SCOREP_USER_ENABLE
9274 scorep_user_region_by_name_begin(
"UR_dpptrf", scorep_user_region_type_common)
9276 CALL dpptrf(
'U',int(nfit,mpl),
globalmatd(imoff+1:),infolp)
9277#ifdef SCOREP_USER_ENABLE
9278 scorep_user_region_by_name_end(
"UR_dpptrf")
9286 WRITE(lun,*)
'No rank defect of the symmetric matrix'
9288 WRITE(lun,*)
'No rank defect of the symmetric block', ib,
' of size', npar
9292 WRITE(*,*)
'Warning: factorization of the symmetric',nfit, &
9293 '-by-',nfit,
' failed at index ', infolp
9294 WRITE(lun,*)
'Warning: factorization of the symmetric',nfit, &
9295 '-by-',nfit,
' failed at index ', infolp
9296 CALL peend(29,
'Aborted, factorization of global matrix failed')
9297 stop
'mdptrf: bad matrix'
9302 IF (nfit > npar)
THEN
9305 IF(infolp /= 0) print *,
' DSPTRS failed: ', infolp
9307 CALL dpptrs(
'U',int(nfit,mpl),1_mpl,
globalmatd(imoff+1:),&
9309 IF(infolp /= 0) print *,
' DPPTRS failed: ', infolp
9313 IF(nfit < npar)
THEN
9334 INTEGER(mpi) :: icoff
9335 INTEGER(mpi) :: ipoff
9338 INTEGER(mpi) :: ncon
9339 INTEGER(mpi) :: nfit
9340 INTEGER(mpi) :: npar
9341 INTEGER(mpl) :: imoff
9342 INTEGER(mpl) :: ioff1
9343 INTEGER(mpl) :: iloff
9344 INTEGER(mpi) :: infolp
9365 WRITE(
lunlog,*)
'Shrinkage of global matrix (A->Q^t*A*Q)'
9385 nfit=npar+ncon;
IF (
icelim > 0) nfit=npar-ncon
9387 IF(nfit < npar)
THEN
9391 CALL dtrtrs(
'L',
'T',
'N',int(ncon,mpl),1_mpl,
lapackql(iloff+npar-ncon+1:),int(npar,mpl),&
9393 IF(infolp /= 0) print *,
' DTRTRS failed: ', infolp
9395 CALL dormql(
'L',
'T',int(npar,mpl),1_mpl,int(ncon,mpl),
lapackql(iloff+1:),int(npar,mpl),&
9397 IF(infolp /= 0) print *,
' DORMQL failed: ', infolp
9416 IF (nfit > npar)
THEN
9419 WRITE(
lunlog,*)
'Factorization of global matrix (A->L*D*L^t)'
9423#ifdef SCOREP_USER_ENABLE
9424 scorep_user_region_by_name_begin(
"UR_dsytrf", scorep_user_region_type_common)
9426 CALL dsytrf(
'U',int(nfit,mpl),
globalmatd(imoff+1:),int(nfit,mpl),&
9428#ifdef SCOREP_USER_ENABLE
9429 scorep_user_region_by_name_end(
"UR_dsytrf")
9436 WRITE(
lunlog,*)
'Factorization of global matrix (A->L*L^t)'
9440#ifdef SCOREP_USER_ENABLE
9441 scorep_user_region_by_name_begin(
"UR_dpotrf", scorep_user_region_type_common)
9443 CALL dpotrf(
'U',int(nfit,mpl),
globalmatd(imoff+1:),int(npar,mpl),infolp)
9444#ifdef SCOREP_USER_ENABLE
9445 scorep_user_region_by_name_end(
"UR_dpotrf")
9453 WRITE(lun,*)
'No rank defect of the symmetric matrix'
9455 WRITE(lun,*)
'No rank defect of the symmetric block', ib,
' of size', npar
9459 WRITE(*,*)
'Warning: factorization of the symmetric',nfit, &
9460 '-by-',nfit,
' failed at index ', infolp
9461 WRITE(lun,*)
'Warning: factorization of the symmetric',nfit, &
9462 '-by-',nfit,
' failed at index ', infolp
9463 CALL peend(29,
'Aborted, factorization of global matrix failed')
9464 stop
'mdutrf: bad matrix'
9469 IF (nfit > npar)
THEN
9470 CALL dsytrs(
'U',int(nfit,mpl),1_mpl,
globalmatd(imoff+1:),int(nfit,mpl),&
9472 IF(infolp /= 0) print *,
' DSYTRS failed: ', infolp
9474 CALL dpotrs(
'U',int(nfit,mpl),1_mpl,
globalmatd(imoff+1:),int(npar,mpl),&
9476 IF(infolp /= 0) print *,
' DPOTRS failed: ', infolp
9480 IF(nfit < npar)
THEN
9485 CALL dormql(
'L',
'N',int(npar,mpl),1_mpl,int(ncon,mpl),
lapackql(iloff+1:),int(npar,mpl),&
9487 IF(infolp /= 0) print *,
' DORMQL failed: ', infolp
9494 iloff=iloff+int(npar,mpl)*int(ncon,mpl)
9516 INTEGER(mpi) :: icboff
9517 INTEGER(mpi) :: icblst
9518 INTEGER(mpi) :: icoff
9519 INTEGER(mpi) :: icfrst
9520 INTEGER(mpi) :: iclast
9521 INTEGER(mpi) :: ipfrst
9522 INTEGER(mpi) :: iplast
9523 INTEGER(mpi) :: ipoff
9526 INTEGER(mpi) :: ncon
9527 INTEGER(mpi) :: npar
9529 INTEGER(mpl) :: imoff
9530 INTEGER(mpl) :: iloff
9531 INTEGER(mpi) :: infolp
9532 INTEGER :: nbopt, ILAENV
9534 REAL(mpd),
INTENT(IN) :: a(mszcon)
9535 REAL(mpd),
INTENT(OUT) :: emin
9536 REAL(mpd),
INTENT(OUT) :: emax
9547 iloff=iloff+int(npar,mpl)*int(ncon,mpl)
9566 DO icb=icboff+1,icboff+icblst
9573 lapackql(iloff+ipfrst:iloff+iplast)=a(imoff+1:imoff+npb)
9590 nbopt = ilaenv( 1_mpl,
'DGEQLF',
'', int(npar,mpl), int(ncon,mpl), int(npar,mpl), -1_mpl )
9591 print *,
'LAPACK optimal block size for DGEQLF:', nbopt
9592 lplwrk=int(ncon,mpl)*int(nbopt,mpl)
9595#ifdef SCOREP_USER_ENABLE
9596 scorep_user_region_by_name_begin(
"UR_dgeqlf", scorep_user_region_type_common)
9598 CALL dgeqlf(int(npar,mpl),int(ncon,mpl),
lapackql(iloff+1:),int(npar,mpl),&
9600 IF(infolp /= 0) print *,
' DGEQLF failed: ', infolp
9601#ifdef SCOREP_USER_ENABLE
9602 scorep_user_region_by_name_end(
"UR_dgeqlf")
9606 iloff=iloff+int(npar,mpl)*int(ncon,mpl)
9609 IF(emax < emin)
THEN
9637 INTEGER(mpi) :: icoff
9638 INTEGER(mpi) :: ipoff
9640 INTEGER(mpi) :: ncon
9641 INTEGER(mpi) :: npar
9642 INTEGER(mpl) :: imoff
9643 INTEGER(mpl) :: iloff
9644 INTEGER(mpi) :: infolp
9645 CHARACTER (LEN=1) :: transr, transl
9647 LOGICAL,
INTENT(IN) :: t
9666 IF(ncon <= 0 ) cycle
9669#ifdef SCOREP_USER_ENABLE
9670 scorep_user_region_by_name_begin(
"UR_dormql", scorep_user_region_type_common)
9678 DO i=ipoff+1,ipoff+npar
9684 CALL dormql(
'R',transr,int(npar,mpl),int(npar,mpl),int(ncon,mpl),
lapackql(iloff+1:),&
9687 IF(infolp /= 0) print *,
' DORMQL failed: ', infolp
9689 CALL dormql(
'L',transl,int(npar,mpl),int(npar,mpl),int(ncon,mpl),
lapackql(iloff+1:),&
9692 IF(infolp /= 0) print *,
' DORMQL failed: ', infolp
9693#ifdef SCOREP_USER_ENABLE
9694 scorep_user_region_by_name_end(
"UR_dormql")
9698 iloff=iloff+int(npar,mpl)*int(ncon,mpl)
9704include
'mkl_pardiso.f90'
9735 TYPE(mkl_pardiso_handle) :: pt(64)
9737 INTEGER(mpl),
PARAMETER :: maxfct =1
9738 INTEGER(mpl),
PARAMETER :: mnum = 1
9739 INTEGER(mpl),
PARAMETER :: nrhs = 1
9741 INTEGER(mpl) :: mtype
9742 INTEGER(mpl) :: phase
9743 INTEGER(mpl) :: error
9744 INTEGER(mpl) :: msglvl
9748 INTEGER(mpl) :: idum(1)
9750 INTEGER(mpl) :: length
9751 INTEGER(mpi) :: nfill
9752 INTEGER(mpi) :: npdblk
9753 REAL(mpd) :: adum(1)
9754 REAL(mpd) :: ddum(1)
9756 INTEGER(mpl) :: iparm(64)
9757 REAL(mpd),
ALLOCATABLE :: b( : )
9758 REAL(mpd),
ALLOCATABLE :: x( : )
9772#ifdef SCOREP_USER_ENABLE
9773 scorep_user_region_by_name_begin(
"UR_mspd00", scorep_user_region_type_common)
9780 WRITE(*,*)
'MSPARDISO: number of rows to fill up = ', nfill
9793 CALL pardiso_64(pt, maxfct, mnum, mtype, phase, int(npdblk,mpl), adum, idum, idum, &
9794 idum, nrhs, iparm, msglvl, ddum, ddum, error)
9795 IF (error /= 0)
THEN
9796 WRITE(lun,*)
'The following ERROR was detected: ', error
9797 WRITE(*,
'(A,2I10)')
' PARDISO release failed (phase, error): ', phase, error
9798 IF (
ipddbg == 0)
WRITE(*,*)
' rerun with "debugPARDISO" for more info'
9799 CALL peend(40,
'Aborted, other error: PARDISO release')
9800 stop
'MSPARDISO: stopping due to error in PARDISO release'
9817 IF (iparm(1) == 0)
WRITE(lun,*)
'PARDISO using defaults '
9818 IF (iparm(43) /= 0)
THEN
9819 WRITE(lun,*)
'PARDISO: computation of the diagonal of inverse matrix not implemented !'
9827#ifdef SCOREP_USER_ENABLE
9828 scorep_user_region_by_name_end(
"UR_mspd00")
9836 WRITE(
lunlog,*)
'Decomposition of global matrix (A->L*D*L^t)'
9843#ifdef SCOREP_USER_ENABLE
9844 scorep_user_region_by_name_begin(
"UR_mspd11", scorep_user_region_type_common)
9853 WRITE(lun,*)
' iparm(',i,
') =', iparm(i)
9856 CALL pardiso_64(pt, maxfct, mnum, mtype, phase, int(npdblk,mpl),
globalmatd,
csr3rowoffsets,
csr3columnlist, &
9857 idum, nrhs, iparm, msglvl, ddum, ddum, error)
9858#ifdef SCOREP_USER_ENABLE
9859 scorep_user_region_by_name_end(
"UR_mspd11")
9862 WRITE(lun,*)
'PARDISO reordering completed ... '
9863 WRITE(lun,*)
'PARDISO peak memory required (KB)', iparm(15)
9866 WRITE(lun,*)
' iparm(',i,
') =', iparm(i)
9869 IF (error /= 0)
THEN
9870 WRITE(lun,*)
'The following ERROR was detected: ', error
9871 WRITE(*,
'(A,2I10)')
' PARDISO decomposition failed (phase, error): ', phase, error
9872 IF (
ipddbg == 0)
WRITE(*,*)
' rerun with "debugPARDISO" for more info'
9873 CALL peend(40,
'Aborted, other error: PARDISO reordering')
9874 stop
'MSPARDISO: stopping due to error in PARDISO reordering'
9876 IF (iparm(60) == 0)
THEN
9881 WRITE(lun,*)
'Size (KB) of allocated memory = ',
ipdmem
9882 WRITE(lun,*)
'Number of nonzeros in factors = ',iparm(18)
9883 WRITE(lun,*)
'Number of factorization MFLOPS = ',iparm(19)
9887#ifdef SCOREP_USER_ENABLE
9888 scorep_user_region_by_name_begin(
"UR_mspd22", scorep_user_region_type_common)
9891 CALL pardiso_64(pt, maxfct, mnum, mtype, phase, int(npdblk,mpl),
globalmatd,
csr3rowoffsets,
csr3columnlist, &
9892 idum, nrhs, iparm, msglvl, ddum, ddum, error)
9893#ifdef SCOREP_USER_ENABLE
9894 scorep_user_region_by_name_end(
"UR_mspd22")
9897 WRITE(lun,*)
'PARDISO factorization completed ... '
9900 WRITE(lun,*)
' iparm(',i,
') =', iparm(i)
9903 IF (error /= 0)
THEN
9904 WRITE(lun,*)
'The following ERROR was detected: ', error
9905 WRITE(*,
'(A,2I10)')
' PARDISO decomposition failed (phase, error): ', phase, error
9906 IF (
ipddbg == 0)
WRITE(*,*)
' rerun with "debugPARDISO" for more info'
9907 CALL peend(40,
'Aborted, other error: PARDISO factorization')
9908 stop
'MSPARDISO: stopping due to error in PARDISO factorization'
9911 IF (iparm(14) > 0) &
9912 WRITE(lun,*)
'Number of perturbed pivots = ',iparm(14)
9913 WRITE(lun,*)
'Number of positive eigenvalues = ',iparm(22)-nfill
9914 WRITE(lun,*)
'Number of negative eigenvalues = ',iparm(23)
9915 ELSE IF (iparm(30) > 0)
THEN
9916 WRITE(lun,*)
'Equation with bad pivot (<=0.) = ',iparm(30)
9925 CALL mpalloc(b,length,
' PARDISO r.h.s')
9926 CALL mpalloc(x,length,
' PARDISO solution')
9929#ifdef SCOREP_USER_ENABLE
9930 scorep_user_region_by_name_begin(
"UR_mspd33", scorep_user_region_type_common)
9934 CALL pardiso_64(pt, maxfct, mnum, mtype, phase, int(npdblk,mpl),
globalmatd,
csr3rowoffsets,
csr3columnlist, &
9935 idum, nrhs, iparm, msglvl, b, x, error)
9936#ifdef SCOREP_USER_ENABLE
9937 scorep_user_region_by_name_end(
"UR_mspd33")
9943 WRITE(lun,*)
'PARDISO solve completed ... '
9944 IF (error /= 0)
THEN
9945 WRITE(lun,*)
'The following ERROR was detected: ', error
9946 WRITE(*,
'(A,2I10)')
' PARDISO decomposition failed (phase, error): ', phase, error
9947 IF (
ipddbg == 0)
WRITE(*,*)
' rerun with "debugPARDISO" for more info'
9948 CALL peend(40,
'Aborted, other error: PARDISO solve')
9949 stop
'MSPARDISO: stopping due to error in PARDISO solve'
9963 INTEGER(mpi) :: iast
9964 INTEGER(mpi) :: idia
9965 INTEGER(mpi) :: imin
9966 INTEGER(mpl) :: ioff1
9968 INTEGER(mpi) :: last
9970 INTEGER(mpi) :: nmax
9971 INTEGER(mpi) :: nmin
9972 INTEGER(mpi) :: ntop
9994 WRITE(
lunlog,*)
'Shrinkage of global matrix (A->Q^t*A*Q)'
10031 DO WHILE(ntop < nmax)
10035 CALL hmpdef(7,real(nmin,mps),real(ntop,mps),
'log10 of positive eigenvalues')
10045 iast=max(1,imin-60)
10046 CALL gmpdef(3,2,
'low-value end of eigenvalues')
10049 CALL gmpxy(3,real(i,mps),evalue)
10063 WRITE(lun,*)
'The first (largest) eigenvalues ...'
10066 WRITE(lun,*)
'The last eigenvalues ... up to',last
10070 WRITE(lun,*)
'The eigenvalues from',
nvgb+1,
' to',
nagb
10074 WRITE(lun,*)
'Log10 + 3 of ',
nfgb,
' eigenvalues in decreasing',
' order'
10075 WRITE(lun,*)
'(for Eigenvalue < 0.001 the value 0.0 is shown)'
10078 'printed for negative eigenvalues'
10081 WRITE(lun,*) last,
' significances: insignificant if ', &
10082 'compatible with N(0,1)'
10111 INTEGER(mpl) :: ioff1
10112 INTEGER(mpl) :: ioff2
10148 INTEGER(mpi) :: istop
10149 INTEGER(mpi) :: itn
10150 INTEGER(mpi) :: itnlim
10151 INTEGER(mpi) :: lun
10152 INTEGER(mpi) :: nout
10153 INTEGER(mpi) :: nrkd
10154 INTEGER(mpi) :: nrkd2
10160 REAL(mpd) :: arnorm
10199 WRITE(lun,*)
'MMINRS: PRECONS ended ', nrkd
10203 globalcorrections, itnlim, nout, rtol, istop, itn, anorm, acond, rnorm, arnorm, ynorm)
10204 ELSE IF(
mbandw > 0)
THEN
10210 WRITE(lun,*)
'MMINRS: EQUDECS ended ', nrkd, nrkd2
10214 globalcorrections, itnlim, nout, rtol, istop, itn, anorm, acond, rnorm, arnorm, ynorm)
10217 globalcorrections, itnlim, nout, rtol, istop, itn, anorm, acond, rnorm, arnorm, ynorm)
10231 IF (
istopa == 0) print *,
'MINRES: istop=0, exact solution x=0.'
10246 INTEGER(mpi) :: istop
10247 INTEGER(mpi) :: itn
10248 INTEGER(mpi) :: itnlim
10249 INTEGER(mpi) :: lun
10250 INTEGER(mpi) :: nout
10251 INTEGER(mpi) :: nrkd
10252 INTEGER(mpi) :: nrkd2
10255 REAL(mpd) :: mxxnrm
10256 REAL(mpd) :: trcond
10266 mxxnrm = real(
nagb,mpd)/sqrt(epsilon(mxxnrm))
10268 trcond = 1.0_mpd/epsilon(trcond)
10269 ELSE IF(
mrmode == 2)
THEN
10299 WRITE(lun,*)
'MMINRS: PRECONS ended ', nrkd
10303 itnlim=itnlim, rtol=rtol, maxxnorm=mxxnrm, trancond=trcond, &
10305 ELSE IF(
mbandw > 0)
THEN
10311 WRITE(lun,*)
'MMINRS: EQUDECS ended ', nrkd, nrkd2
10316 itnlim=itnlim, rtol=rtol, maxxnorm=mxxnrm, trancond=trcond, &
10320 itnlim=itnlim, rtol=rtol, maxxnorm=mxxnrm, trancond=trcond, &
10335 IF (
istopa == 3) print *,
'MINRES: istop=0, exact solution x=0.'
10351 INTEGER(mpi),
INTENT(IN) :: n
10352 REAL(mpd),
INTENT(IN) :: x(n)
10353 REAL(mpd),
INTENT(OUT) :: y(n)
10373 INTEGER(mpi),
INTENT(IN) :: n
10374 REAL(mpd),
INTENT(IN) :: x(n)
10375 REAL(mpd),
INTENT(OUT) :: y(n)
10406 REAL(mps) :: concu2
10407 REAL(mps) :: concut
10408 REAL,
DIMENSION(2) :: ta
10411 INTEGER(mpi) :: iact
10412 INTEGER(mpi) :: iagain
10413 INTEGER(mpi) :: idx
10414 INTEGER(mpi) :: info
10416 INTEGER(mpi) :: ipoff
10417 INTEGER(mpi) :: icoff
10418 INTEGER(mpl) :: ioff
10419 INTEGER(mpi) :: itgbi
10420 INTEGER(mpi) :: ivgbi
10421 INTEGER(mpi) :: jcalcm
10423 INTEGER(mpi) :: labelg
10424 INTEGER(mpi) :: litera
10425 INTEGER(mpl) :: lrej
10426 INTEGER(mpi) :: lun
10427 INTEGER(mpi) :: lunp
10428 INTEGER(mpi) :: minf
10429 INTEGER(mpi) :: mrati
10430 INTEGER(mpi) :: nan
10431 INTEGER(mpi) :: ncon
10432 INTEGER(mpi) :: nfaci
10433 INTEGER(mpi) :: nloopsol
10434 INTEGER(mpi) :: npar
10435 INTEGER(mpi) :: nrati
10436 INTEGER(mpl) :: nrej
10437 INTEGER(mpi) :: nsol
10438 INTEGER(mpi) :: inone
10440 INTEGER(mpi) :: infolp
10441 INTEGER(mpi) :: nfit
10442 INTEGER(mpl) :: imoff
10446 REAL(mpd) :: dratio
10447 REAL(mpd) :: dwmean
10456 LOGICAL :: warnerss
10457 LOGICAL :: warners3
10459 CHARACTER (LEN=7) :: cratio
10460 CHARACTER (LEN=7) :: cfacin
10461 CHARACTER (LEN=7) :: crjrat
10472 WRITE(lunp,*)
'Solution algorithm: '
10473 WRITE(lunp,121)
'=================================================== '
10476 WRITE(lunp,121)
'solution method:',
'matrix inversion'
10477 ELSE IF(
metsol == 2)
THEN
10478 WRITE(lunp,121)
'solution method:',
'diagonalization'
10479 ELSE IF(
metsol == 3)
THEN
10480 WRITE(lunp,121)
'solution method:',
'decomposition'
10481 ELSE IF(
metsol == 4)
THEN
10482 WRITE(lunp,121)
'solution method:',
'minres (Paige/Saunders)'
10483 ELSE IF(
metsol == 5)
THEN
10484 WRITE(lunp,121)
'solution method:',
'minres-qlp (Choi/Paige/Saunders)'
10486 WRITE(lunp,121)
' ',
' using QR factorization'
10487 ELSE IF(
mrmode == 2)
THEN
10488 WRITE(lunp,121)
' ',
' using QLP factorization'
10490 WRITE(lunp,121)
' ',
' using QR and QLP factorization'
10491 WRITE(lunp,123)
'transition condition',
mrtcnd
10493 ELSE IF(
metsol == 6)
THEN
10494 WRITE(lunp,121)
'solution method:', &
10495 'gmres (generalized minimzation of residuals)'
10497 ELSE IF(
metsol == 7)
THEN
10499 WRITE(lunp,121)
'solution method:',
'LAPACK factorization (DSPTRF)'
10501 WRITE(lunp,121)
'solution method:',
'LAPACK factorization (DPPTRF)'
10503 IF(
ilperr == 1)
WRITE(lunp,121)
' ',
'with error calculation (D??TRI)'
10504 ELSE IF(
metsol == 8)
THEN
10506 WRITE(lunp,121)
'solution method:',
'LAPACK factorization (DSYTRF)'
10508 WRITE(lunp,121)
'solution method:',
'LAPACK factorization (DPOTRF)'
10510 IF(
ilperr == 1)
WRITE(lunp,121)
' ',
'with error calculation (D??TRI)'
10512 ELSE IF(
metsol == 9)
THEN
10514 WRITE(lunp,121)
'solution method:',
'Intel oneMKL PARDISO (sparse matrix (CSR3))'
10516 WRITE(lunp,121)
'solution method:',
'Intel oneMKL PARDISO (sparse matrix (BSR3))'
10521 WRITE(lunp,123)
'convergence limit at Delta F=',
dflim
10522 WRITE(lunp,122)
'maximum number of iterations=',
mitera
10525 WRITE(lunp,122)
'matrix recalculation up to ',
matrit,
'. iteration'
10529 WRITE(lunp,121)
'matrix storage:',
'full'
10530 ELSE IF(
matsto == 2)
THEN
10531 WRITE(lunp,121)
'matrix storage:',
'sparse'
10533 WRITE(lunp,122)
'pre-con band-width parameter=',
mbandw
10535 WRITE(lunp,121)
'pre-conditioning:',
'default'
10536 ELSE IF(
mbandw < 0)
THEN
10537 WRITE(lunp,121)
'pre-conditioning:',
'none!'
10538 ELSE IF(
mbandw > 0)
THEN
10540 WRITE(lunp,121)
'pre-conditioning=',
'skyline-matrix (rank preserving)'
10542 WRITE(lunp,121)
'pre-conditioning=',
'band-matrix'
10547 WRITE(lunp,121)
'using pre-sigmas:',
'no'
10550 WRITE(lunp,124)
'pre-sigmas defined for', &
10551 REAL(100*npresg,mps)/
REAL(nvgb,mps),
' % of variable parameters'
10552 WRITE(lunp,123)
'default pre-sigma=',
regpre
10555 WRITE(lunp,121)
'regularization:',
'no'
10557 WRITE(lunp,121)
'regularization:',
'yes'
10558 WRITE(lunp,123)
'regularization factor=',
regula
10562 WRITE(lunp,121)
'Chi square cut equiv 3 st.dev applied'
10563 WRITE(lunp,123)
'... in first iteration with factor',
chicut
10564 WRITE(lunp,123)
'... in second iteration with factor',
chirem
10565 WRITE(lunp,121)
' (reduced by sqrt in next iterations)'
10568 WRITE(lunp,121)
'Scaling of measurement errors applied'
10569 WRITE(lunp,123)
'... factor for "global" measuements',
dscerr(1)
10570 WRITE(lunp,123)
'... factor for "local" measuements',
dscerr(2)
10573 WRITE(lunp,122)
'Down-weighting of outliers in',
lhuber,
' iterations'
10574 WRITE(lunp,123)
'Cut on downweight fraction',
dwcut
10578121
FORMAT(1x,a40,3x,a)
10579122
FORMAT(1x,a40,3x,i0,a)
10580123
FORMAT(1x,a40,2x,e9.2)
10581124
FORMAT(1x,a40,3x,f5.1,a)
10591 stepl =real(stp,mps)
10603 ELSE IF(
metsol == 2)
THEN
10606 ELSE IF(
metsol == 3)
THEN
10609 ELSE IF(
metsol == 4)
THEN
10612 ELSE IF(
metsol == 5)
THEN
10615 ELSE IF(
metsol == 6)
THEN
10627 WRITE(
lunlog,*)
'Checking feasibility of parameters:'
10628 WRITE(*,*)
'Checking feasibility of parameters:'
10629 CALL feasib(concut,iact)
10631 WRITE(*,102) concut
10632 WRITE(*,*)
' parameters are made feasible'
10633 WRITE(
lunlog,102) concut
10634 WRITE(
lunlog,*)
' parameters are made feasible'
10636 WRITE(*,*)
' parameters are feasible (i.e. satisfy constraints)'
10637 WRITE(
lunlog,*)
' parameters are feasible (i.e. satisfy constraints)'
10645 WRITE(*,*)
'Reading files and accumulating vectors/matrices ...'
10649 WRITE(
lunlog,*)
'Reading files and accumulating vectors/matrices ...'
10666 IF(jcalcm+1 /= 0)
THEN
10675 IF(
iterat /= litera)
THEN
10696 IF (iabs(jcalcm) <= 1)
THEN
10709 IF(3*nrej >
nrecal)
THEN
10711 WRITE(*,*)
'Data records rejected in previous loop: '
10713 WRITE(*,*)
'Too many rejects (>33.3%) - stop'
10714 CALL peend(26,
'Aborted, too many rejects')
10721 IF(abs(
icalcm) == 1)
THEN
10733 ELSE IF(
metsol == 2)
THEN
10735 ELSE IF(
metsol == 3)
THEN
10737 ELSE IF(
metsol == 4)
THEN
10739 ELSE IF(
metsol == 5)
THEN
10741 ELSE IF(
metsol == 6)
THEN
10742 WRITE(*,*)
'... reserved for GMRES (not yet!)'
10745 ELSE IF(
metsol == 7)
THEN
10747 ELSE IF(
metsol == 8)
THEN
10750 ELSE IF(
metsol == 9)
THEN
10764 CALL feasib(concut,iact)
10776 angras=real(db/sqrt(db1*db2),mps)
10777 dbsig=16.0_mpd*sqrt(max(db1,db2))*epsilon(db)
10783 lsflag=lsflag .AND. (db > dbsig)
10790 ELSE IF(
metsol == 2)
THEN
10793 ELSE IF(
metsol == 3)
THEN
10796 ELSE IF(
metsol == 4)
THEN
10799 ELSE IF(
metsol == 5)
THEN
10802 ELSE IF(
metsol == 6)
THEN
10812 IF(db <= -dbsig)
THEN
10813 WRITE(*,*)
'Function not decreasing:',db
10814 IF(db > -1.0e-3_mpd)
THEN
10816 IF (iagain <= 1)
THEN
10817 WRITE(*,*)
'... again matrix calculation'
10821 WRITE(*,*)
'... aborting iterations'
10825 WRITE(*,*)
'... stopping iterations'
10848 IF (
nloopn == nloopsol)
THEN
10854 stepl=real(stp,mps)
10864 WRITE(*,*)
'Result vector containes ', nan,
' NaNs - stop'
10865 CALL peend(25,
'Aborted, result vector contains NaNs')
10872 WRITE(*,*)
'Subito! Exit after first step.'
10877 WRITE(*,*)
'INFO=0 should not happen (line search input err)'
10878 IF (iagain <= 0)
THEN
10883 IF(info < 0 .OR.
nloopn == nloopsol) cycle
10887 CALL feasib(concut,iact)
10888 IF(iact /= 0.OR.
chicut > 1.0)
THEN
10903 IF(sum(
nrejec) /= 0)
THEN
10905 WRITE(*,*)
'Data records rejected in last loop: '
10927 nfit=npar+ncon;
IF (
icelim > 0) nfit=npar-ncon
10928 IF (nfit > npar)
THEN
10931 WRITE(
lunlog,*)
'Inverse of global matrix from LDLt factorization'
10936#ifdef SCOREP_USER_ENABLE
10937 scorep_user_region_by_name_begin(
"UR_dsptri", scorep_user_region_type_common)
10940 IF(infolp /= 0) print *,
' DSPTRI failed: ', infolp
10941#ifdef SCOREP_USER_ENABLE
10942 scorep_user_region_by_name_end(
"UR_dsptri")
10948#ifdef SCOREP_USER_ENABLE
10949 scorep_user_region_by_name_begin(
"UR_dsytri", scorep_user_region_type_common)
10951 CALL dsytri(
'U',int(nfit,mpl),
globalmatd(imoff+1:),int(nfit,mpl),&
10953 IF(infolp /= 0) print *,
' DSYTRI failed: ', infolp
10954#ifdef SCOREP_USER_ENABLE
10955 scorep_user_region_by_name_end(
"UR_dsytri")
10962 WRITE(
lunlog,*)
'Inverse of global matrix from LLt factorization'
10967#ifdef SCOREP_USER_ENABLE
10968 scorep_user_region_by_name_begin(
"UR_dpptri", scorep_user_region_type_common)
10970 CALL dpptri(
'U',int(nfit,mpl),
globalmatd(imoff+1:),infolp)
10971 IF(infolp /= 0) print *,
' DPPTRI failed: ', infolp
10972#ifdef SCOREP_USER_ENABLE
10973 scorep_user_region_by_name_end(
"UR_dpptri")
10978#ifdef SCOREP_USER_ENABLE
10979 scorep_user_region_by_name_begin(
"UR_dpotri", scorep_user_region_type_common)
10981 CALL dpotri(
'U',int(nfit,mpl),
globalmatd(imoff+1:),int(npar,mpl),infolp)
10982 IF(infolp /= 0) print *,
' DPOTRI failed: ', infolp
10983#ifdef SCOREP_USER_ENABLE
10984 scorep_user_region_by_name_end(
"UR_dpotri")
11002 DO i=npar-ncon+1,npar
11009 WRITE(
lunlog,*)
'Expansion of global matrix (A->Q*A*Q^t)'
11025 catio=real(dratio,mps)
11029 mrati=nint(100.0*catio,mpi)
11033 IF (
nfilw <= 0)
THEN
11034 WRITE(lunp,*)
'Sum(Chi^2)/Sum(Ndf) =',
fvalue
11036 WRITE(lunp,*)
' =',dratio
11038 WRITE(lunp,*)
'Sum(W*Chi^2)/Sum(Ndf)/<W> =',
fvalue
11040 WRITE(lunp,*)
' /',dwmean
11041 WRITE(lunp,*)
' =',dratio
11045 ' with correction for down-weighting ',catio
11066 nrati=nint(10000.0*real(nrej,mps)/real(
nrecal,mps),mpi)
11067 WRITE(crjrat,197) 0.01_mpd*real(nrati,mpd)
11068 nfaci=nint(100.0*sqrt(catio),mpi)
11070 WRITE(cratio,197) 0.01_mpd*real(mrati,mpd)
11071 WRITE(cfacin,197) 0.01_mpd*real(nfaci,mpd)
11074 IF(mrati < 90.OR.mrati > 110) warner=.true.
11075 IF(nrati > 100) warner=.true.
11076 IF(
ncgbe /= 0) warner=.true.
11078 IF(
nalow /= 0) warners=.true.
11080 IF(
nmiss1 /= 0) warnerss=.true.
11081 IF(iagain /= 0) warnerss=.true.
11082 IF(
ndefec /= 0) warnerss=.true.
11083 IF(
ndefpg /= 0) warnerss=.true.
11085 IF(
nrderr /= 0) warners3=.true.
11087 IF(warner.OR.warners.OR.warnerss.Or.warners3)
THEN
11090 WRITE(*,199)
'WarningWarningWarningWarningWarningWarningWarningWarningWar'
11091 WRITE(*,199)
'arningWarningWarningWarningWarningWarningWarningWarningWarn'
11092 WRITE(*,199)
'rningWarningWarningWarningWarningWarningWarningWarningWarni'
11093 WRITE(*,199)
'ningWarningWarningWarningWarningWarningWarningWarningWarnin'
11094 WRITE(*,199)
'ingWarningWarningWarningWarningWarningWarningWarningWarning'
11095 WRITE(*,199)
'ngWarningWarningWarningWarningWarningWarningWarningWarningW'
11096 WRITE(*,199)
'gWarningWarningWarningWarningWarningWarningWarningWarningWa'
11098 IF(mrati < 90.OR.mrati > 110)
THEN
11100 WRITE(*,*)
' Chi^2/Ndf = ',cratio,
' (should be close to 1)'
11101 WRITE(*,*)
' => multiply all input standard ', &
11102 'deviations by factor',cfacin
11105 IF(nrati > 100)
THEN
11107 WRITE(*,*)
' Fraction of rejects =',crjrat,
' %', &
11108 ' (should be far below 1 %)'
11109 WRITE(*,*)
' => please provide correct mille data'
11113 IF(iagain /= 0)
THEN
11115 WRITE(*,*)
' Matrix not positiv definite '// &
11116 '(function not decreasing)'
11117 WRITE(*,*)
' => please provide correct mille data'
11122 WRITE(*,*)
' Rank defect =',
ndefec, &
11123 ' for global matrix, should be 0'
11124 WRITE(*,*)
' => please provide correct mille data'
11129 WRITE(*,*)
' Rank defect for',
ndefpg, &
11130 ' parameter groups, should be 0'
11131 WRITE(*,*)
' => please provide correct mille data'
11136 WRITE(*,*)
' Rank defect =',
nmiss1, &
11137 ' for constraint equations, should be 0'
11138 WRITE(*,*)
' => please correct constraint definition'
11141 IF(
ncgbe /= 0)
THEN
11143 WRITE(*,*)
' Number of empty constraints =',
ncgbe,
', should be 0'
11144 WRITE(*,*)
' => please check constraint definition, mille data'
11147 IF(
nxlow /= 0)
THEN
11149 WRITE(*,*)
' Possible rank defects =',
nxlow,
' for global matrix'
11150 WRITE(*,*)
' (too few accepted entries)'
11151 WRITE(*,*)
' => please check mille data and ENTRIES cut'
11154 IF(
nalow /= 0)
THEN
11156 WRITE(*,*)
' Possible bad elements =',
nalow,
' in global vector'
11157 WRITE(*,*)
' (toos few accepted entries)'
11158 IF(
ipcntr > 0)
WRITE(*,*)
' (indicated in millepede.res by counts<0)'
11159 WRITE(*,*)
' => please check mille data and ENTRIES cut'
11164 WRITE(*,*)
' Binary file(s) with read errors =',
nrderr,
' (treated as EOF)'
11165 WRITE(*,*)
' => please check mille data'
11169 WRITE(*,199)
'WarningWarningWarningWarningWarningWarningWarningWarningWar'
11170 WRITE(*,199)
'arningWarningWarningWarningWarningWarningWarningWarningWarn'
11171 WRITE(*,199)
'rningWarningWarningWarningWarningWarningWarningWarningWarni'
11172 WRITE(*,199)
'ningWarningWarningWarningWarningWarningWarningWarningWarnin'
11173 WRITE(*,199)
'ingWarningWarningWarningWarningWarningWarningWarningWarning'
11174 WRITE(*,199)
'ngWarningWarningWarningWarningWarningWarningWarningWarningW'
11175 WRITE(*,199)
'gWarningWarningWarningWarningWarningWarningWarningWarningWa'
11186 ELSE IF(
metsol == 2)
THEN
11188 ELSE IF(
metsol == 3)
THEN
11194 IF(labelg == 0) cycle
11195 itgbi=inone(labelg)
11198 IF(ivgbi < 0) ivgbi=0
11199 IF(ivgbi == 0) cycle
11208 ELSE IF(
metsol == 6)
THEN
11211 ELSE IF(
metsol == 7)
THEN
11219 CALL peend(4,
'Ended with severe warnings (bad binary file(s))')
11220 ELSE IF (warnerss)
THEN
11221 CALL peend(3,
'Ended with severe warnings (bad global matrix)')
11222 ELSE IF (warners)
THEN
11223 CALL peend(2,
'Ended with severe warnings (insufficient measurements)')
11224 ELSE IF (warner)
THEN
11225 CALL peend(1,
'Ended with warnings (bad measurements)')
11227 CALL peend(0,
'Ended normally')
11230102
FORMAT(
' Call FEASIB with cut=',g10.3)
11248 INTEGER(mpi) :: kfl
11249 INTEGER(mpi) :: kmin
11250 INTEGER(mpi) :: kmax
11251 INTEGER(mpi) :: nrc
11252 INTEGER(mpl) :: nrej
11258 REAL(mpd) :: sumallw
11259 REAL(mpd) :: sumrejw
11261 sumallw=0.; sumrejw=0.;
11270 sumallw=sumallw+real(nrc,mpd)*
wfd(kfl)
11271 sumrejw=sumrejw+real(nrej,mpd)*
wfd(kfl)
11272 frac=real(nrej,mps)/real(nrc,mps)
11273 IF (frac > fmax)
THEN
11277 IF (frac < fmin)
THEN
11284 WRITE(*,
"(' Weighted fraction =',F8.2,' %')") 100.*sumrejw/sumallw
11285 IF (
nfilb > 1)
THEN
11286 WRITE(*,
"(' File with max. fraction ',I6,' :',F8.2,' %')") kmax, 100.*fmax
11287 WRITE(*,
"(' File with min. fraction ',I6,' :',F8.2,' %')") kmin, 100.*fmin
11313 INTEGER(mpi) :: iargc
11316 INTEGER(mpi) :: ierrf
11317 INTEGER(mpi) :: ieq
11318 INTEGER(mpi) :: ifilb
11319 INTEGER(mpi) :: ioff
11320 INTEGER(mpi) :: iopt
11321 INTEGER(mpi) :: ios
11322 INTEGER(mpi) :: iosum
11325 INTEGER(mpi) :: mat
11326 INTEGER(mpi) :: nab
11327 INTEGER(mpi) :: nline
11328 INTEGER(mpi) :: npat
11329 INTEGER(mpi) :: ntext
11331 INTEGER(mpi) :: nuf
11332 INTEGER(mpi) :: nums
11333 INTEGER(mpi) :: nufile
11334 INTEGER(mpi) :: lenfileInfo
11335 INTEGER(mpi) :: lenFileNames
11336 INTEGER(mpi) :: matint
11337 INTEGER(mpi),
DIMENSION(:,:),
ALLOCATABLE :: vecfileInfo
11338 INTEGER(mpi),
DIMENSION(:,:),
ALLOCATABLE :: tempArray
11339 INTEGER(mpl) :: rows
11340 INTEGER(mpl) :: cols
11341 INTEGER(mpl) :: newcols
11342 INTEGER(mpl) :: length
11344 CHARACTER (LEN=1024) :: text
11345 CHARACTER (LEN=1024) :: fname
11346 CHARACTER (LEN=14) :: bite(3)
11347 CHARACTER (LEN=32) :: keystx
11348 INTEGER(mpi),
PARAMETER :: mnum=100
11349 REAL(mpd) :: dnum(mnum)
11353 SUBROUTINE initc(nfiles)
BIND(c)
11355 INTEGER(c_int),
INTENT(IN),
VALUE :: nfiles
11356 END SUBROUTINE initc
11361 DATA bite/
'C_binary',
'text ',
'Fortran_binary'/
11376 WRITE(*,*)
'Command line options: '
11377 WRITE(*,*)
'--------------------- '
11379 CALL getarg(i,text)
11380 CALL rltext(text,ia,ib,nab)
11381 WRITE(*,101) i,text(1:nab)
11382 IF(text(ia:ia) /=
'-')
THEN
11383 nu=nufile(text(ia:ib))
11386 WRITE(*,*)
'Second text file in command line - stop'
11387 CALL peend(12,
'Aborted, second text file in command line')
11393 WRITE(*,*)
'Open error for file:',text(ia:ib),
' - stop'
11394 CALL peend(16,
'Aborted, open error for file')
11395 IF(text(ia:ia) /=
'/')
THEN
11396 CALL getenv(
'PWD',text)
11397 CALL rltext(text,ia,ib,nab)
11398 WRITE(*,*)
'PWD:',text(ia:ib)
11403 IF(index(text(ia:ib),
'b') /= 0)
THEN
11405 WRITE(*,*)
'Debugging requested'
11407 it=index(text(ia:ib),
't')
11410 ieq=index(text(ia+it:ib),
'=')+it
11411 IF (it /= ieq)
THEN
11412 IF (index(text(ia+ieq:ib),
'SL0' ) /= 0)
ictest=2
11413 IF (index(text(ia+ieq:ib),
'SLE' ) /= 0)
ictest=3
11414 IF (index(text(ia+ieq:ib),
'BP' ) /= 0)
ictest=4
11415 IF (index(text(ia+ieq:ib),
'BRLF') /= 0)
ictest=5
11416 IF (index(text(ia+ieq:ib),
'BRLC') /= 0)
ictest=6
11419 IF(index(text(ia:ib),
's') /= 0)
isubit=1
11420 IF(index(text(ia:ib),
'f') /= 0)
iforce=1
11421 IF(index(text(ia:ib),
'c') /= 0)
icheck=1
11422 IF(index(text(ia:ib),
'C') /= 0)
icheck=2
11424 IF(i == iargc())
WRITE(*,*)
'--------------------- '
11445 CALL rltext(text,ia,ib,nab)
11446 nu=nufile(text(ia:ib))
11450 CALL peend(10,
'Aborted, no steering file')
11451 stop
'in FILETC: no steering file. .'
11464 WRITE(*,*)
'Listing of steering file: ',
filnam(1:
nfnam)
11465 WRITE(*,*)
'-------------------------'
11468 WRITE(*,*)
'Open error for steering file - stop'
11469 CALL peend(11,
'Aborted, open error for steering file')
11470 IF(
filnam(1:1) /=
'/')
THEN
11471 CALL getenv(
'PWD',text)
11472 CALL rltext(text,ia,ib,nab)
11473 WRITE(*,*)
'PWD:',text(ia:ib)
11482 rows=6; cols=lenfileinfo
11483 CALL mpalloc(vecfileinfo,rows,cols,
'file info from steering')
11486 READ(10,102,iostat=ierrf) text
11487 IF (ierrf < 0)
EXIT
11488 CALL rltext(text,ia,ib,nab)
11490 IF(nline <= 50)
THEN
11491 WRITE(*,101) nline,text(1:nab)
11492 IF(nline == 50)
WRITE(*,*)
' ...'
11496 CALL rltext(text,ia,ib,nab)
11497 IF(ib == ia+2)
THEN
11498 mat=matint(text(ia:ib),
'end',npat,ntext)
11499 IF(mat == max(npat,ntext))
THEN
11502 WRITE(*,*)
' end-statement after',nline,
' text lines'
11507 keystx=
'fortranfiles'
11508 mat=matint(text(ia:ib),keystx,npat,ntext)
11509 IF(mat == max(npat,ntext))
THEN
11516 mat=matint(text(ia:ib),keystx,npat,ntext)
11517 IF(mat == max(npat,ntext))
THEN
11523 keystx=
'closeandreopen'
11524 mat=matint(text(ia:ib),keystx,npat,ntext)
11525 IF(mat == max(npat,ntext))
THEN
11533 iopt=index(text(ia:ib),
' -- ')
11534 IF (iopt > 0) ie=iopt-1
11537 nu=nufile(text(ia:ie))
11539 IF (
nfiles == lenfileinfo)
THEN
11540 CALL mpalloc(temparray,rows,cols,
'temp file info from steering')
11541 temparray=vecfileinfo
11543 lenfileinfo=lenfileinfo*2
11544 newcols=lenfileinfo
11545 CALL mpalloc(vecfileinfo,rows,newcols,
'file info from steering')
11546 vecfileinfo(:,1:cols)=temparray(:,1:cols)
11552 lenfilenames=lenfilenames+ie-ia+1
11553 vecfileinfo(1,
nfiles)=nline
11554 vecfileinfo(2,
nfiles)=nu
11555 vecfileinfo(3,
nfiles)=ia
11556 vecfileinfo(4,
nfiles)=ie
11557 vecfileinfo(5,
nfiles)=iopt
11558 vecfileinfo(6,
nfiles)=ib
11568 CALL mpalloc(
nfd,length,
'file line (in steering)')
11571 length=lenfilenames
11577 READ(10,102,iostat=ierrf) text
11578 IF (ierrf < 0)
EXIT
11580 IF (nline == vecfileinfo(1,i))
THEN
11581 nfd(i)=vecfileinfo(1,i)
11582 mfd(i)=vecfileinfo(2,i)
11583 ia=vecfileinfo(3,i)-1
11584 lfd(i)=vecfileinfo(4,i)-ia
11586 tfd(ioff+k)=text(ia+k:ia+k)
11591 IF (vecfileinfo(5,i) > 0)
THEN
11592 CALL ratext(text(vecfileinfo(5,i)+4:vecfileinfo(6,i)),nums,dnum,mnum)
11593 IF (nums > 0)
ofd(i)=real(dnum(1),mps)
11603 CALL mpalloc(
ifd,length,
'integrated record numbers (=offset)')
11604 CALL mpalloc(
jfd,length,
'number of accepted records')
11605 CALL mpalloc(
kfd,rows,length,
'number of records in file, file order')
11610 CALL mpalloc(
sfd,rows,length,
'start, end of file name in TFD')
11611 CALL mpalloc(
yfd,length,
'modification date')
11614 WRITE(*,*)
'-------------------------'
11620 WRITE(*,*)
'Table of files:'
11621 WRITE(*,*)
'---------------'
11624 WRITE(8,*)
'Text and data files:'
11628 fname(k:k)=
tfd(ioff+k)
11631 IF (
mprint > 1)
WRITE(*,103) i,bite(
mfd(i)),fname(1:
lfd(i))
11632 WRITE(8,103) i,bite(
mfd(i)),fname(1:
lfd(i))
11636 WRITE(*,*)
'---------------'
11650 IF(
mfd(i) == 3)
THEN
11674 IF(
mfd(i) == 1)
THEN
11695 WRITE(*,*)
'Opening of C-files not supported.'
11711 IF(iosum /= 0)
THEN
11712 CALL peend(15,
'Aborted, open error(s) for binary files')
11713 stop
'FILETC: open error '
11715 IF(
nfilb == 0)
THEN
11716 CALL peend(14,
'Aborted, no binary files')
11717 stop
'FILETC: no binary files '
11720 WRITE(*,*)
nfilb,
' binary files opened'
11722 WRITE(*,*)
nfilb,
' binary files opened and closed'
11726103
FORMAT(i3,2x,a14,3x,a)
11789 INTEGER(mpi) :: ierrf
11790 INTEGER(mpi) :: ioff
11791 INTEGER(mpi) :: ios
11792 INTEGER(mpi) :: iosum
11794 INTEGER(mpi) :: mat
11795 INTEGER(mpi) :: nab
11796 INTEGER(mpi) :: nfiln
11797 INTEGER(mpi) :: nline
11798 INTEGER(mpi) :: nlinmx
11799 INTEGER(mpi) :: npat
11800 INTEGER(mpi) :: ntext
11801 INTEGER(mpi) :: matint
11805 CHARACTER (LEN=1024) :: text
11806 CHARACTER (LEN=1024) :: fname
11809 WRITE(*,*)
'Processing text files ...'
11822 IF(
mfd(i) /= 2) cycle
11824 fname(k:k)=
tfd(ia+k)
11826 WRITE(*,*)
'File ',fname(1:
lfd(i))
11827 IF (
mprint > 1)
WRITE(*,*)
' '
11828 OPEN(10,file=fname(1:
lfd(i)),iostat=ios,form=
'FORMATTED')
11830 WRITE(*,*)
'Open error for file ',fname(1:
lfd(i))
11840 READ(10,102,iostat=ierrf) text
11841 IF (ierrf < 0)
THEN
11844 WRITE(*,*)
' end-of-file after',nline,
' text lines'
11848 IF(nline <= nlinmx.AND.
mprint > 1)
THEN
11849 CALL rltext(text,ia,ib,nab)
11851 WRITE(*,101) nline,text(1:nab)
11852 IF(nline == nlinmx)
WRITE(*,*)
' ...'
11855 CALL rltext(text,ia,ib,nab)
11856 IF(ib == ia+2)
THEN
11857 mat=matint(text(ia:ib),
'end',npat,ntext)
11858 IF(mat == max(npat,ntext))
THEN
11861 WRITE(*,*)
' end-statement after',nline,
' text lines'
11867 IF(nfiln <=
nfiles)
THEN
11868 IF(nline ==
nfd(nfiln))
THEN
11883 IF(iosum /= 0)
THEN
11884 CALL peend(16,
'Aborted, open error(s) for text files')
11885 stop
'FILETX: open error(s) in text files '
11888 WRITE(*,*)
'... end of text file processing.'
11893 WRITE(*,*)
lunkno,
' unknown keywords in steering files, ', &
11894 'or file non-existing,'
11895 WRITE(*,*)
' see above!'
11896 WRITE(*,*)
'------------> stop'
11898 CALL peend(13,
'Aborted, unknown keywords in steering file')
11907 ELSE IF(
matsto == 1)
THEN
11909 ELSE IF(
matsto == 2)
THEN
11912 ELSE IF(
metsol == 1)
THEN
11914 ELSE IF(
metsol == 2)
THEN
11916 ELSE IF(
metsol == 3)
THEN
11918 ELSE IF(
metsol == 4)
THEN
11920 ELSE IF(
metsol == 5)
THEN
11922 ELSE IF(
metsol == 6)
THEN
11925 ELSE IF(
metsol == 7)
THEN
11927 ELSE IF(
metsol == 8)
THEN
11930 ELSE IF(
metsol == 9)
THEN
11935 WRITE(*,*)
'MINRES forced with sparse matrix!'
11937 WRITE(*,*)
'MINRES forced with sparse matrix!'
11939 WRITE(*,*)
'MINRES forced with sparse matrix!'
11944 WRITE(*,*)
'MINRES forced with sparse matrix!'
11946 WRITE(*,*)
'MINRES forced with sparse matrix!'
11948 WRITE(*,*)
'MINRES forced with sparse matrix!'
11956 WRITE(*,*)
'Solution method and matrix-storage mode:'
11958 WRITE(*,*)
' METSOL = 1: matrix inversion'
11959 ELSE IF(
metsol == 2)
THEN
11960 WRITE(*,*)
' METSOL = 2: diagonalization'
11961 ELSE IF(
metsol == 3)
THEN
11962 WRITE(*,*)
' METSOL = 3: decomposition'
11963 ELSE IF(
metsol == 4)
THEN
11964 WRITE(*,*)
' METSOL = 4: MINRES'
11965 ELSE IF(
metsol == 5)
THEN
11966 WRITE(*,*)
' METSOL = 5: MINRES-QLP'
11967 ELSE IF(
metsol == 6)
THEN
11968 WRITE(*,*)
' METSOL = 6: GMRES (-> MINRES)'
11970 ELSE IF(
metsol == 7)
THEN
11971 WRITE(*,*)
' METSOL = 7: LAPACK factorization'
11972 ELSE IF(
metsol == 8)
THEN
11973 WRITE(*,*)
' METSOL = 8: LAPACK factorization'
11975 ELSE IF(
metsol == 9)
THEN
11976 WRITE(*,*)
' METSOL = 9: Intel oneMKL PARDISO'
11981 WRITE(*,*)
' with',
mitera,
' iterations'
11984 WRITE(*,*)
' MATSTO = 0: unpacked symmetric matrix, ',
'n*n elements'
11985 ELSEIF(
matsto == 1)
THEN
11986 WRITE(*,*)
' MATSTO = 1: full symmetric matrix, ',
'(n*n+n)/2 elements'
11987 ELSE IF(
matsto == 2)
THEN
11988 WRITE(*,*)
' MATSTO = 2: sparse matrix (custom)'
11989 ELSE IF(
matsto == 3)
THEN
11991 WRITE(*,*)
' MATSTO = 3: sparse matrix (upper triangle, CSR3)'
11993 WRITE(*,*)
' MATSTO = 3: sparse matrix (upper triangle, BSR3)'
11997 WRITE(*,*)
' and band matrix, width',
mbandw
12001 WRITE(*,*)
'Chi square cut equiv 3 st.dev applied ...'
12002 WRITE(*,*)
' in first iteration with factor',
chicut
12003 WRITE(*,*)
' in second iteration with factor',
chirem
12004 WRITE(*,*)
' (reduced by sqrt in next iterations)'
12008 WRITE(*,*)
' Down-weighting of outliers in',
lhuber,
' iterations'
12009 WRITE(*,*)
' Cut on downweight fraction',
dwcut
12012 WRITE(*,*)
'Iterations (solutions) with line search:'
12021 WRITE(*,*)
' All with Chi square cut scaling factor <= 1.'
12052 INTEGER(mpi) :: ios
12056 INTEGER(mpi) :: npat
12057 INTEGER(mpi) :: ntext
12058 INTEGER(mpi) :: nuprae
12061 CHARACTER (LEN=*),
INTENT(INOUT) :: fname
12067 IF(len(fname) > 5)
THEN
12068 IF(fname(1:5) ==
'rfio:') nuprae=1
12069 IF(fname(1:5) ==
'dcap:') nuprae=2
12070 IF(fname(1:5) ==
'root:') nuprae=3
12072 IF(nuprae == 0)
THEN
12073 INQUIRE(file=fname,iostat=ios,exist=ex)
12074 IF(ios /= 0)
nufile=-abs(ios)
12075 IF(ios /= 0)
RETURN
12076 ELSE IF(nuprae == 1)
THEN
12089 nm=
matint(
'xt',fname(l1:ll),npat,ntext)
12092 nm=
matint(
'tx',fname(l1:ll),npat,ntext)
12113 INTEGER(mpi) :: ier
12114 INTEGER(mpi) :: iomp
12117 INTEGER(mpi) :: kkey
12118 INTEGER(mpi) :: label
12119 INTEGER(mpi) :: lkey
12120 INTEGER(mpi) :: mat
12121 INTEGER(mpi) :: miter
12122 INTEGER(mpi) :: nab
12123 INTEGER(mpi) :: nkey
12124 INTEGER(mpi) :: nkeys
12126 INTEGER(mpi) :: nmeth
12127 INTEGER(mpi) :: npat
12128 INTEGER(mpi) :: ntext
12129 INTEGER(mpi) :: nums
12130 INTEGER(mpi) :: matint
12132 CHARACTER (LEN=*),
INTENT(IN) :: text
12133 INTEGER(mpi),
INTENT(IN) :: nline
12137 parameter(nkeys=7,nmeth=10)
12139 parameter(nkeys=6,nmeth=9)
12142 parameter(nkeys=6,nmeth=7)
12144 CHARACTER (LEN=16) :: methxt(nmeth)
12145 CHARACTER (LEN=16) :: keylst(nkeys)
12146 CHARACTER (LEN=32) :: keywrd
12147 CHARACTER (LEN=32) :: keystx
12148 CHARACTER (LEN=itemCLen) :: ctext
12149 INTEGER(mpi),
PARAMETER :: mnum=100
12150 REAL(mpd) :: dnum(mnum)
12153 INTEGER(mpi) :: ipvs
12156 INTEGER(mpi) :: lpvs
12160 SUBROUTINE additem(length,list,label,value)
12162 INTEGER(mpi),
INTENT(IN OUT) :: length
12163 TYPE(listitem),
DIMENSION(:),
INTENT(IN OUT),
ALLOCATABLE :: list
12164 INTEGER(mpi),
INTENT(IN) :: label
12165 REAL(mpd),
INTENT(IN) :: value
12167 SUBROUTINE additemc(length,list,label,text)
12169 INTEGER(mpi),
INTENT(IN OUT) :: length
12170 TYPE(listitemc),
DIMENSION(:),
INTENT(IN OUT),
ALLOCATABLE :: list
12171 INTEGER(mpi),
INTENT(IN) :: label
12172 CHARACTER(LEN = itemCLen),
INTENT(IN) :: text
12174 SUBROUTINE additemi(length,list,label,ivalue)
12176 INTEGER(mpi),
INTENT(IN OUT) :: length
12177 TYPE(listitemi),
DIMENSION(:),
INTENT(IN OUT),
ALLOCATABLE :: list
12178 INTEGER(mpi),
INTENT(IN) :: label
12179 INTEGER(mpi),
INTENT(IN) :: ivalue
12186 DATA keylst/
'unknown',
'parameter',
'constraint',
'measurement',
'method',
'comment',
'pardiso'/
12187 DATA methxt/
'diagonalization',
'inversion',
'fullMINRES',
'sparseMINRES', &
12188 'fullMINRES-QLP',
'sparseMINRES-QLP',
'decomposition',
'fullLAPACK',
'unpackedLAPACK', &
12191 DATA keylst/
'unknown',
'parameter',
'constraint',
'measurement',
'method',
'comment'/
12192 DATA methxt/
'diagonalization',
'inversion',
'fullMINRES',
'sparseMINRES', &
12193 'fullMINRES-QLP',
'sparseMINRES-QLP',
'decomposition',
'fullLAPACK',
'unpackedLAPACK'/
12196 DATA keylst/
'unknown',
'parameter',
'constraint',
'measurement',
'method',
'comment'/
12197 DATA methxt/
'diagonalization',
'inversion',
'fullMINRES',
'sparseMINRES', &
12198 'fullMINRES-QLP',
'sparseMINRES-QLP',
'decomposition'/
12204 CALL rltext(text,ia,ib,nab)
12205 IF(nab == 0)
GOTO 10
12206 CALL ratext(text(1:nab),nums,dnum,mnum)
12208 IF(nums /= 0) nkey=0
12216 keystx=keylst(nkey)
12217 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12218 IF(100*mat >= 80*max(npat,ntext))
GO TO 10
12224 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12225 IF(100*mat >= 80*max(npat,ntext))
THEN
12227 IF(nums > 0) mprint=nint(dnum(1),mpi)
12232 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12233 IF(100*mat >= 80*max(npat,ntext))
THEN
12236 IF(nums > 0) mdebug=nint(dnum(1),mpi)
12237 IF(nums > 1) mdebg2=nint(dnum(2),mpi)
12242 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12243 IF(100*mat >= 80*max(npat,ntext))
THEN
12244 IF(nums > 0 .AND. dnum(1) > 0.5) mreqenf=nint(dnum(1),mpi)
12245 IF(nums > 1 .AND. dnum(2) > 0.5) mreqena=nint(dnum(2),mpi)
12246 IF(nums > 2 .AND. dnum(3) > 0.5) iteren=nint(dnum(1)*dnum(3),mpi)
12250 keystx=
'printrecord'
12251 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12252 IF(100*mat >= 80*max(npat,ntext))
THEN
12253 IF(nums > 0) nrecpr=nint(dnum(1),mpi)
12254 IF(nums > 1) nrecp2=nint(dnum(2),mpi)
12259 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12260 IF(100*mat >= 80*max(npat,ntext))
THEN
12261 IF (nums > 0.AND.dnum(1) > 0.) mxrec=nint(dnum(1),mpi)
12266 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12267 IF(100*mat >= 80*max(npat,ntext))
THEN
12268 IF (nums > 0.AND.dnum(1) >= 0.) ncache=nint(dnum(1),mpi)
12269 IF (nums == 2.AND.dnum(2) > 0..AND.dnum(2) <= 1.0) &
12270 fcache(1)=real(dnum(2),mps)
12271 IF (nums >= 4)
THEN
12273 fcache(k)=real(dnum(k+1),mps)
12280 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12281 IF(100*mat >= 80*max(npat,ntext))
THEN
12286 chicut=real(dnum(1),mps)
12287 IF(chicut < 1.0) chicut=-1.0
12291 chirem=real(dnum(2),mps)
12292 IF(chirem < 1.0) chirem=1.0
12293 IF(chicut >= 1.0) chirem=min(chirem,chicut)
12301 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12302 IF(100*mat >= 80*max(npat,ntext))
THEN
12303 IF(nums > 0) chhuge=real(dnum(1),mps)
12304 IF(chhuge < 1.0) chhuge=1.0
12309 keystx=
'linesearch'
12310 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12311 IF(100*mat >= 80*max(npat,ntext))
THEN
12312 IF(nums > 0) lsearch=nint(dnum(1),mpi)
12317 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12318 IF(100*mat >= 80*max(npat,ntext))
THEN
12319 IF(nums > 0) lfitnp=nint(dnum(1),mpi)
12320 IF(nums > 1) lfitbb=nint(dnum(2),mpi)
12324 keystx=
'regularization'
12325 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12326 IF(100*mat >= 80*max(npat,ntext))
THEN
12328 regula=real(dnum(1),mps)
12329 IF(nums >= 2) regpre=real(dnum(2),mps)
12333 keystx=
'regularisation'
12334 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12335 IF(100*mat >= 80*max(npat,ntext))
THEN
12337 regula=real(dnum(1),mps)
12338 IF(nums >= 2) regpre=real(dnum(2),mps)
12343 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12344 IF(100*mat >= 80*max(npat,ntext))
THEN
12345 regpre=real(dnum(1),mps)
12350 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12351 IF(100*mat >= 80*max(npat,ntext))
THEN
12352 matrit=nint(dnum(1),mpi)
12357 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12358 IF(100*mat >= 80*max(npat,ntext))
THEN
12360 IF (nums > 0.AND.dnum(1) > 0.) matmon=nint(dnum(1),mpi)
12365 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12366 IF(100*mat >= 80*max(npat,ntext))
THEN
12367 IF(nums > 0) mbandw=nint(dnum(1),mpi)
12368 IF(mbandw < 0) mbandw=-1
12369 IF(nums > 1) lprecm=nint(dnum(2),mpi)
12393 keystx=
'outlierdownweighting'
12394 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12395 IF(100*mat >= 80*max(npat,ntext))
THEN
12396 lhuber=nint(dnum(1),mpi)
12397 IF(lhuber > 0.AND.lhuber <= 2) lhuber=2
12401 keystx=
'dwfractioncut'
12402 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12403 IF(100*mat >= 80*max(npat,ntext))
THEN
12404 dwcut=real(dnum(1),mps)
12405 IF(dwcut > 0.5) dwcut=0.5
12409 keystx=
'maxlocalcond'
12410 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12411 IF(100*mat >= 80*max(npat,ntext))
THEN
12412 IF (nums > 0.AND.dnum(1) > 0.0) cndlmx=real(dnum(1),mps)
12417 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12418 IF(100*mat >= 80*max(npat,ntext))
THEN
12419 prange=abs(real(dnum(1),mps))
12424 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12425 IF(100*mat >= 80*max(npat,ntext))
THEN
12431 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12432 IF(100*mat >= 80*max(npat,ntext))
THEN
12437 keystx=
'memorydebug'
12438 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12439 IF(100*mat >= 80*max(npat,ntext))
THEN
12441 IF (nums > 0.AND.dnum(1) > 0.0) memdbg=nint(dnum(1),mpi)
12445 keystx=
'globalcorr'
12446 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12447 IF(100*mat >= 80*max(npat,ntext))
THEN
12452 keystx=
'printcounts'
12453 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12454 IF(100*mat >= 80*max(npat,ntext))
THEN
12456 IF (nums > 0) ipcntr=nint(dnum(1),mpi)
12460 keystx=
'weightedcons'
12461 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12462 IF(100*mat >= 80*max(npat,ntext))
THEN
12464 IF (nums > 0) iwcons=nint(dnum(1),mpi)
12468 keystx=
'skipemptycons'
12469 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12470 IF(100*mat >= 80*max(npat,ntext))
THEN
12475 keystx=
'resolveredundancycons'
12476 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12477 IF(100*mat >= 80*max(npat,ntext))
THEN
12482 keystx=
'withelimination'
12483 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12484 IF(100*mat >= 80*max(npat,ntext))
THEN
12489 keystx=
'postprocessing'
12490 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12491 IF(100*mat >= 80*max(npat,ntext))
THEN
12492 lenpostproc=ib-
keyb-1
12493 cpostproc(1:lenpostproc)=text(
keyb+2:ib)
12498 keystx=
'withLAPACKelimination'
12499 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12500 IF(100*mat >= 80*max(npat,ntext))
THEN
12506 keystx=
'withmultipliers'
12507 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12508 IF(100*mat >= 80*max(npat,ntext))
THEN
12513 keystx=
'checkinput'
12514 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12515 IF(100*mat >= 80*max(npat,ntext))
THEN
12517 IF (nums > 0) icheck=nint(dnum(1),mpi)
12521 keystx=
'checkparametergroups'
12522 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12523 IF(100*mat >= 80*max(npat,ntext))
THEN
12528 keystx=
'monitorresiduals'
12529 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12530 IF(100*mat >= 80*max(npat,ntext))
THEN
12532 IF (nums > 0) imonit=nint(dnum(1),mpi)
12533 IF (nums > 1) measbins=max(measbins,nint(dnum(2),mpi))
12537 keystx=
'monitorpulls'
12538 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12539 IF(100*mat >= 80*max(npat,ntext))
THEN
12542 IF (nums > 0) imonit=nint(dnum(1),mpi)
12543 IF (nums > 1) measbins=max(measbins,nint(dnum(2),mpi))
12547 keystx=
'monitorprogress'
12548 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12549 IF(100*mat >= 80*max(npat,ntext))
THEN
12552 IF (nums > 0) monpg1=max(1,nint(dnum(1),mpi))
12553 IF (nums > 1) monpg2=max(1,nint(dnum(2),mpi))
12557 keystx=
'scaleerrors'
12558 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12559 IF(100*mat >= 80*max(npat,ntext))
THEN
12561 IF (nums > 0) dscerr(1:2)=dnum(1)
12562 IF (nums > 1) dscerr(2)=dnum(2)
12566 keystx=
'iterateentries'
12567 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12568 IF(100*mat >= 80*max(npat,ntext))
THEN
12569 iteren=huge(iteren)
12570 IF (nums > 0) iteren=nint(dnum(1),mpi)
12575 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12576 IF(100*mat >= 80*max(npat,ntext))
THEN
12584 WRITE(*,*)
'WARNING: multithreading not available'
12590 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12591 IF(100*mat >= 80*max(npat,ntext))
THEN
12592 WRITE(*,*)
'WARNING: keyword COMPRESS is obsolete (compression is default)'
12604 keystx=
'countrecords'
12605 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12606 IF(100*mat >= 80*max(npat,ntext))
THEN
12612 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12613 IF(100*mat >= 80*max(npat,ntext).AND.mnrsel < 100)
THEN
12614 nl=min(nums,100-mnrsel)
12616 lbmnrs(mnrsel+k)=nint(dnum(k),mpi)
12622 keystx=
'pairentries'
12623 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12624 IF(100*mat >= 80*max(npat,ntext))
THEN
12627 IF (nums > 0.AND.dnum(1) > 0.0)
THEN
12628 mreqpe=nint(dnum(1),mpi)
12629 IF (nums >= 2.AND.dnum(2) >= dnum(1)) mhispe=nint(dnum(2),mpi)
12630 IF (nums >= 3.AND.dnum(3) >= dnum(1)) msngpe=nint(dnum(3),mpi)
12636 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12637 IF(100*mat >= 80*max(npat,ntext))
THEN
12638 wolfc1=real(dnum(1),mps)
12639 wolfc2=real(dnum(2),mps)
12646 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12647 IF(100*mat >= 80*max(npat,ntext))
THEN
12649 IF (dnum(1) < 1.0e-10_mpd.OR.dnum(1) > 1.0e-04_mpd)
THEN
12650 WRITE(*,*)
'ERROR: need 1.0D-10 <= MRESTL ', &
12651 '<= 1.0D-04, but get ', dnum(1)
12660 keystx=
'mrestranscond'
12661 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12662 IF(100*mat >= 80*max(npat,ntext))
THEN
12670 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12671 IF(100*mat >= 80*max(npat,ntext))
THEN
12673 mrmode = int(dnum(1),mpi)
12678 keystx=
'nofeasiblestart'
12679 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12680 IF(100*mat >= 80*max(npat,ntext))
THEN
12686 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12687 IF(100*mat >= 80*max(npat,ntext))
THEN
12692 keystx=
'readerroraseof'
12693 mat=matint(text(ia:ib),keystx,npat,ntext)
12694 IF(100*mat >= 80*max(npat,ntext))
THEN
12700 keystx=
'LAPACKwitherrors'
12701 mat=matint(text(ia:ib),keystx,npat,ntext)
12702 IF(100*mat >= 80*max(npat,ntext))
THEN
12707 keystx=
'debugPARDISO'
12708 mat=matint(text(ia:ib),keystx,npat,ntext)
12709 IF(100*mat >= 80*max(npat,ntext))
THEN
12714 keystx=
'blocksizePARDISO'
12715 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12716 IF(100*mat >= 80*max(npat,ntext).AND.mnrsel < 100)
THEN
12717 nl=min(nums,10-mpdbsz)
12719 IF (nint(dnum(k),mpi) > 0)
THEN
12720 IF (mpdbsz == 0)
THEN
12722 ipdbsz(mpdbsz)=nint(dnum(k),mpi)
12723 ELSE IF (nint(dnum(k),mpi) > ipdbsz(mpdbsz))
THEN
12725 ipdbsz(mpdbsz)=nint(dnum(k),mpi)
12733 keystx=
'fortranfiles'
12734 mat=matint(text(ia:ib),keystx,npat,ntext)
12735 IF(mat == max(npat,ntext))
RETURN
12738 mat=matint(text(ia:ib),keystx,npat,ntext)
12739 IF(mat == max(npat,ntext))
RETURN
12741 keystx=
'closeandreopen'
12742 mat=matint(text(ia:ib),keystx,npat,ntext)
12743 IF(mat == max(npat,ntext))
RETURN
12747 IF(nums /= 0) nkey=0
12750 WRITE(*,*)
'**************************************************'
12752 WRITE(*,*)
'Unknown keyword(s): ',text(1:min(nab,50))
12754 WRITE(*,*)
'**************************************************'
1277610
IF(nkey > 0)
THEN
12780 lpvs=nint(dnum(1),mpi)
12782 CALL additem(lenparameters,listparameters,lpvs,dnum(2))
12783 CALL additem(lenpresigmas,listpresigmas,lpvs,dnum(3))
12785 WRITE(*,*)
'Line',nline,
' error, label=',lpvs
12787 ELSE IF(nums /= 0)
THEN
12789 WRITE(*,*)
'Wrong text in line',nline
12790 WRITE(*,*)
'Status: new parameter'
12791 WRITE(*,*)
'> ',text(1:nab)
12793 ELSE IF(lkey == 3)
THEN
12795 IF(nums >= 1.AND.nums <= 2)
THEN
12797 CALL additem(lenconstraints,listconstraints,lpvs,dnum(1))
12799 IF(iwcons > 0) lpvs=-2
12801 IF(nums == 2) plvs=dnum(2)
12802 CALL additem(lenconstraints,listconstraints,lpvs,plvs)
12805 WRITE(*,*)
'Wrong text in line',nline
12806 WRITE(*,*)
'Status: new keyword constraint'
12807 WRITE(*,*)
'> ',text(1:nab)
12809 ELSE IF(lkey == 4)
THEN
12811 nummeasurements=nummeasurements+1
12813 CALL additem(lenmeasurements,listmeasurements,lpvs,dnum(1))
12815 CALL additem(lenmeasurements,listmeasurements,lpvs,dnum(2))
12818 WRITE(*,*)
'Wrong text in line',nline
12819 WRITE(*,*)
'Status: new keyword measurement'
12820 WRITE(*,*)
'> ',text(1:nab)
12822 ELSE IF(lkey == 5.AND.
keyb <
keyc)
THEN
12824 IF(nums >= 1) miter=nint(dnum(1),mpi)
12825 IF(miter >= 1) mitera=miter
12826 dflim=real(dnum(2),mps)
12830 mat=matint(text(
keyb+1:
keyc),keystx,npat,ntext)
12831 IF(100*mat >= 80*max(npat,ntext))
THEN
12835 ELSE IF(i == 2)
THEN
12838 ELSE IF(i == 3)
THEN
12841 ELSE IF(i == 4)
THEN
12844 ELSE IF(i == 5)
THEN
12847 ELSE IF(i == 6)
THEN
12850 ELSE IF(i == 7)
THEN
12854 ELSE IF(i == 8)
THEN
12857 ELSE IF(i == 9)
THEN
12861 ELSE IF(i == 10)
THEN
12870 ELSE IF(nkey == 0)
THEN
12873 lpvs=nint(dnum(1),mpi)
12875 CALL additem(lenparameters,listparameters,lpvs,dnum(2))
12876 CALL additem(lenpresigmas,listpresigmas,lpvs,dnum(3))
12878 WRITE(*,*)
'Line',nline,
' error, label=',lpvs
12880 ELSE IF(nums > 1.AND.nums < 3)
THEN
12882 WRITE(*,*)
'Wrong text in line',nline
12883 WRITE(*,*)
'Status continuation parameter'
12884 WRITE(*,*)
'> ',text(1:nab)
12887 ELSE IF(lkey == 3)
THEN
12890 label=nint(dnum(i),mpi)
12891 IF(label <= 0) ier=1
12893 IF(mod(nums,2) /= 0) ier=1
12896 lpvs=nint(dnum(i),mpi)
12898 CALL additem(lenconstraints,listconstraints,lpvs,plvs)
12902 WRITE(*,*)
'Wrong text in line',nline
12903 WRITE(*,*)
'Status continuation constraint'
12904 WRITE(*,*)
'> ',text(1:nab)
12907 ELSE IF(lkey == 4)
THEN
12911 label=nint(dnum(i),mpi)
12912 IF(label <= 0) ier=1
12914 IF(mod(nums,2) /= 0) ier=1
12918 lpvs=nint(dnum(i),mpi)
12920 CALL additem(lenmeasurements,listmeasurements,lpvs,plvs)
12924 WRITE(*,*)
'Wrong text in line',nline
12925 WRITE(*,*)
'Status continuation measurement'
12926 WRITE(*,*)
'> ',text(1:nab)
12928 ELSE IF(lkey == 6)
THEN
12930 lpvs=nint(dnum(1),mpi)
12934 IF (text(j:j) ==
' ')
EXIT
12937 CALL additemc(lencomments,listcomments,lpvs,ctext)
12939 WRITE(*,*)
'Line',nline,
' error, label=',lpvs
12941 ELSE IF(nums /= 0)
THEN
12943 WRITE(*,*)
'Wrong text in line',nline
12944 WRITE(*,*)
'Status: continuation comment'
12945 WRITE(*,*)
'> ',text(1:nab)
12949 ELSE IF(lkey == 7)
THEN
12952 label=nint(dnum(i),mpi)
12953 IF(label <= 0.OR.label > 64) ier=1
12955 IF(mod(nums,2) /= 0) ier=1
12959 lpvs=nint(dnum(i),mpi)
12960 ipvs=nint(dnum(i+1),mpi)
12961 CALL additemi(lenpardiso,listpardiso,lpvs,ipvs)
12965 WRITE(*,*)
'Wrong text in line',nline
12966 WRITE(*,*)
'Status continuation measurement'
12967 WRITE(*,*)
'> ',text(1:nab)
12986 INTEGER(mpi),
INTENT(IN OUT) :: length
12987 TYPE(
listitem),
DIMENSION(:),
INTENT(IN OUT),
ALLOCATABLE :: list
12988 INTEGER(mpi),
INTENT(IN) :: label
12989 REAL(mpd),
INTENT(IN) :: value
12991 INTEGER(mpl) :: newSize
12992 INTEGER(mpl) :: oldSize
12993 TYPE(
listitem),
DIMENSION(:),
ALLOCATABLE :: tempList
12995 IF (label > 0.AND.
value == 0.)
RETURN
12996 IF (length == 0 )
THEN
12998 CALL mpalloc(list,newsize,
' list ')
13000 oldsize=
size(list,kind=
mpl)
13001 IF (length >= oldsize)
THEN
13002 newsize = oldsize + oldsize/5 + 100
13003 CALL mpalloc(templist,oldsize,
' temp. list ')
13006 CALL mpalloc(list,newsize,
' list ')
13007 list(1:oldsize)=templist(1:oldsize)
13012 list(length)%label=label
13013 list(length)%value=
value
13028 INTEGER(mpi),
INTENT(IN OUT) :: length
13029 TYPE(
listitemc),
DIMENSION(:),
INTENT(IN OUT),
ALLOCATABLE :: list
13030 INTEGER(mpi),
INTENT(IN) :: label
13031 CHARACTER(len = itemCLen),
INTENT(IN) :: text
13033 INTEGER(mpl) :: newSize
13034 INTEGER(mpl) :: oldSize
13035 TYPE(
listitemc),
DIMENSION(:),
ALLOCATABLE :: tempList
13037 IF (label > 0.AND.text ==
'')
RETURN
13038 IF (length == 0 )
THEN
13040 CALL mpalloc(list,newsize,
' list ')
13042 oldsize=
size(list,kind=
mpl)
13043 IF (length >= oldsize)
THEN
13044 newsize = oldsize + oldsize/5 + 100
13045 CALL mpalloc(templist,oldsize,
' temp. list ')
13048 CALL mpalloc(list,newsize,
' list ')
13049 list(1:oldsize)=templist(1:oldsize)
13054 list(length)%label=label
13055 list(length)%text=text
13070 INTEGER(mpi),
INTENT(IN OUT) :: length
13071 TYPE(
listitemi),
DIMENSION(:),
INTENT(IN OUT),
ALLOCATABLE :: list
13072 INTEGER(mpi),
INTENT(IN) :: label
13073 INTEGER(mpi),
INTENT(IN) :: ivalue
13075 INTEGER(mpl) :: newSize
13076 INTEGER(mpl) :: oldSize
13077 TYPE(
listitemi),
DIMENSION(:),
ALLOCATABLE :: tempList
13079 IF (length == 0 )
THEN
13081 CALL mpalloc(list,newsize,
' list ')
13083 oldsize=
size(list,kind=
mpl)
13084 IF (length >= oldsize)
THEN
13085 newsize = oldsize + oldsize/5 + 100
13086 CALL mpalloc(templist,oldsize,
' temp. list ')
13089 CALL mpalloc(list,newsize,
' list ')
13090 list(1:oldsize)=templist(1:oldsize)
13095 list(length)%label=label
13096 list(length)%ivalue=ivalue
13110 CHARACTER (LEN=*),
INTENT(IN) :: text
13111 CHARACTER (LEN=16) :: textc
13120 textl(ka:kb)=text(1:l)
13124 textc=text(1:l)//
'-end'
13132 textl(ka:kb)=textc(1:l)
13159 INTEGER(mpi),
INTENT(IN) :: lun
13160 CHARACTER (LEN=*),
INTENT(IN) :: fname
13161 CHARACTER (LEN=33) :: nafile
13162 CHARACTER (LEN=33) :: nbfile
13168 CALL peend(17,
'Aborted, file name too long')
13169 stop
'File name too long '
13172 nafile(l+1:l+1)=
'~'
13174 INQUIRE(file=nafile(1:l),exist=ex)
13176 INQUIRE(file=nafile(1:l+1),exist=ex)
13178 CALL system(
'rm '//nafile)
13181 nafile(l+1:l+1)=
' '
13182 CALL system(
'mv '//nafile//nbfile)
13184 OPEN(unit=lun,file=fname)
13195 REAL,
DIMENSION(2) :: ta
13215 IF(ncount > 1)
THEN
13217 nsecd1=int(delta,
mpi)
13219 minut1=nsecd1/60-60*nhour1
13220 secnd1=delta-60*(minut1+60*nhour1)
13222 nsecd2=int(delta,
mpi)
13224 minut2=nsecd2/60-60*nhour2
13225 secnd2=delta-60*(minut2+60*nhour2)
13226 WRITE(*,101) nhour1,minut1,secnd1, nhour2,minut2,secnd2
13231101
FORMAT(i4,
' h',i3,
' min',f5.1,
' sec total',18x,
'elapsed', &
13232 i4,
' h',i3,
' min',f5.1,
' sec')
13246 INTEGER(mpi),
INTENT(IN) :: icode
13247 CHARACTER (LEN=*),
INTENT(IN) :: cmessage
13249 CALL mvopen(9,
'millepede.end')
13250 WRITE(9,101) icode, cmessage
13251101
FORMAT(1x,i4,3x,a)
13255END SUBROUTINE peend
13267 INTEGER(mpi),
INTENT(IN) :: kfile
13268 INTEGER(mpi),
INTENT(IN) :: ithr
13269 INTEGER(mpi),
INTENT(OUT) :: ierr
13271 INTEGER(mpi),
DIMENSION(13) :: ibuff
13272 INTEGER(mpi) :: ioff
13273 INTEGER(mpi) :: ios
13275 INTEGER(mpi) :: lfn
13276 INTEGER(mpi) :: lun
13277 INTEGER(mpi) :: moddate
13278 CHARACTER (LEN=1024) :: fname
13279 CHARACTER (LEN=7) :: cfile
13284 SUBROUTINE openc(filename, lfn, lun, ios)
BIND(c)
13286 CHARACTER(kind=c_char),
DIMENSION(*),
INTENT(IN) :: filename
13287 INTEGER(c_int),
INTENT(IN),
VALUE :: lfn
13288 INTEGER(c_int),
INTENT(IN),
VALUE :: lun
13289 INTEGER(c_int),
INTENT(INOUT) :: ios
13290 END SUBROUTINE openc
13302 fname(k:k)=
tfd(ioff+k)
13307 IF(kfile <=
nfilf)
THEN
13310 OPEN(lun,file=fname(1:lfn),iostat=ios, form=
'UNFORMATTED')
13311 print *,
' lun ', lun, ios
13315 CALL openc(fname(1:lfn),lfn,lun,ios)
13317 WRITE(*,*)
'Opening of C-files not supported.'
13324 WRITE(*,*)
'Open error for file ',fname(1:lfn), ios
13325 IF (moddate /= 0)
THEN
13326 WRITE(cfile,
'(I7)') kfile
13327 CALL peend(15,
'Aborted, open error(s) for binary file ' // cfile)
13328 stop
'PEREAD: open error'
13333 ios=stat(fname(1:lfn),ibuff)
13337 WRITE(*,*)
'STAT error for file ',fname(1:lfn), ios
13341 IF (moddate /= 0)
THEN
13342 IF (ibuff(10) /= moddate)
THEN
13343 WRITE(cfile,
'(I7)') kfile
13344 CALL peend(19,
'Aborted, binary file modified (date) ' // cfile)
13345 stop
'PEREAD: file modified'
13348 yfd(kfile)=ibuff(10)
13363 INTEGER(mpi),
INTENT(IN) :: kfile
13364 INTEGER(mpi),
INTENT(IN) :: ithr
13366 INTEGER(mpi) :: lun
13370 SUBROUTINE closec(lun)
BIND(c)
13372 INTEGER(c_int),
INTENT(IN),
VALUE :: lun
13379 IF(kfile <=
nfilf)
THEN
13398 INTEGER(mpi),
INTENT(IN) :: kfile
13400 INTEGER(mpi) :: lun
13404 SUBROUTINE resetc(lun)
BIND(c)
13406 INTEGER(c_int),
INTENT(IN),
VALUE :: lun
13412 IF (kfile <=
nfilf)
THEN
13431 INTEGER(mpi) :: ipgrp
13432 INTEGER(mpi) :: irank
13433 INTEGER(mpi) :: isize
13434 INTEGER(mpi) :: ivoff
13435 INTEGER(mpi) :: itgbi
13437 INTEGER(mpi) :: msize
13438 INTEGER(mpi),
PARAMETER :: mxsize = 1000
13440 INTEGER(mpl):: length
13442 REAL(mpd),
DIMENSION(:),
ALLOCATABLE :: auxVectorD
13443 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: auxVectorI
13444 REAL(mpd),
DIMENSION(:),
ALLOCATABLE :: resParGroup
13445 REAL(mpd),
DIMENSION(:),
ALLOCATABLE :: blockParGroup
13453 IF (isize <= mxsize)
THEN
13454 msize=max(msize,isize)
13456 print *,
' CKPGRP: par. group', ipgrp,
' not checked -- too large: ', isize
13459 IF (msize == 0)
RETURN
13462 length=int(msize,mpl)*(int(msize,mpl)+1)/2
13463 CALL mpalloc(blockpargroup,length,
'(matrix) block for parameter groups (D)')
13465 CALL mpalloc(respargroup,length,
'residuals for parameter groups (D)')
13466 CALL mpalloc(auxvectori,length,
'auxiliary array (I)')
13467 CALL mpalloc(auxvectord,length,
'auxiliary array (D)')
13471 print *,
' CKPGRP par. group first label size rank'
13474 IF (isize > mxsize) cycle
13481 blockpargroup(ij)=matij(ivoff+i,ivoff+j)
13485 CALL sqminv(blockpargroup,respargroup,isize,irank, auxvectord, auxvectori)
13488 IF (isize == irank)
THEN
13492 print *,
' CKPGRP ', ipgrp,
globalparlabelindex(1,itgbi), isize, irank,
' rank deficit !!!'
13510 INTEGER(mpl) :: nan
13511 INTEGER(mpl) :: neg
13513 print *,
' Checking global matrix(D) for NANs ',
size(
globalmatd,kind=mpl)
13518 print *,
' i, nan ', i, nan
13524 print *,
' Checking diagonal elements ',
nagb
13529 print *,
' i, neg ', i, neg
13533 print *,
' CHKMAT summary ', nan, neg
13555 REAL(mpd),
INTENT(IN) :: chi2
13556 INTEGER(mpi),
INTENT(IN) :: ithrd
13557 INTEGER(mpi),
INTENT(IN) :: ndf
13558 REAL(mpd),
INTENT(IN) :: dw
13560 INTEGER(mpl) ::nadd
13588 REAL(mpd),
INTENT(OUT) ::chi2
13589 INTEGER(mpl),
INTENT(OUT) ::ndf
13590 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