905#ifdef SCOREP_USER_ENABLE
906#include "scorep/SCOREP_User.inc"
931 REAL,
DIMENSION(2) :: ta
934 INTEGER(mpi) :: iopnmp
945 INTEGER(mpi) :: nsecnd
946 INTEGER(mpi) :: ntsec
948 CHARACTER (LEN=24) :: chdate
949 CHARACTER (LEN=24) :: chost
951 CHARACTER (LEN=6) :: c6
952 INTEGER major, minor, patch
976 CALL getenv(
'HOSTNAME',chost)
977 IF (chost(1:1) ==
' ')
CALL getenv(
'HOST',chost)
978 WRITE(*,*)
'($Id: 2abb4a313debe0a8077a24e1c46a3c0875f9041b $)'
983 CALL ilaver( major,minor, patch )
984 WRITE(*,110) lapack64, major,minor, patch
985110
FORMAT(
' using LAPACK64 with ',(a),
', version ',i0,
'.',i0,
'.',i0)
987 WRITE(*,*)
'using Intel oneMKL PARDISO'
991 WRITE(*,111) __gnuc__ , __gnuc_minor__ , __gnuc_patchlevel__
992111
FORMAT(
' compiled with gcc ',i0,
'.',i0,
'.',i0)
995 WRITE(*,111) __pgic__ , __pgic_minor__ , __pgic_patchlevel__
996111
FORMAT(
' compiled with pgi ',i0,
'.',i0,
'.',i0)
998#ifdef SCOREP_USER_ENABLE
999 WRITE(*,*)
'instrumenting Score-P user regions'
1002 WRITE(*,*)
' < Millepede II-P starting ... ',chdate
1003 WRITE(*,*)
' ',chost
1006 WRITE(8,*)
'($Id: 2abb4a313debe0a8077a24e1c46a3c0875f9041b $)'
1008 WRITE(8,*)
'Log-file Millepede II-P ', chdate
1009 WRITE(8,*)
' ', chost
1011 CALL peend(-1,
'Still running or crashed')
1020 WRITE(*,*)
'!!! Checking input only, no calculation of a solution !!!'
1021 WRITE(8,*)
'!!! Checking input only, no calculation of a solution !!!'
1040 CALL getenv(
'OMP_NUM_THREADS',c6)
1042 CALL getenv(lapack64//
'_NUM_THREADS',c6)
1044 IF (c6(1:1) ==
' ')
THEN
1046 WRITE(*,*)
'Number of threads for LAPACK: unkown (empty OMP_NUM_THREADS)'
1048 WRITE(*,*)
'Number of threads for LAPACK: unkown (empty ',lapack64//
'_NUM_THREADS)'
1051 WRITE(*,*)
'Number of threads for LAPACK: ', c6
1071 CALL mvopen(lun,
'millepede.his')
1077 CALL mvopen(1,
'mpdebug.txt')
1081 times(0)=rstext-rstp
1087 times(1)=rloop1-rstext
1091 WRITE(8,*)
'Chi square cut equiv 3 st.dev applied ...'
1092 WRITE(8,*)
' in first iteration with factor',
chicut
1093 WRITE(8,*)
' in second iteration with factor',
chirem
1094 WRITE(8,*)
' (reduced by sqrt in next iterations)'
1098 WRITE(8,*)
'Down-weighting of outliers in',
lhuber,
' iterations'
1099 WRITE(8,*)
'Cut on downweight fraction',
dwcut
1103 times(2)=rloop2-rloop1
1108 CALL peend(5,
'Ended without solution (empty constraints)')
1110 CALL peend(0,
'Ended normally')
1147 CALL gmpdef(6,1,
'log10(#records) vs file number')
1148 CALL gmpdef(7,1,
'final rejection fraction vs file number')
1150 'final <Chi^2/Ndf> from accepted local fits vs file number')
1151 CALL gmpdef(9,1,
'<Ndf> from accepted local fits vs file number')
1157 rej=real(nrc-
jfd(kfl),mps)/real(nrc,mps)
1158 CALL gmpxy(6,real(kfl,mps),log10(real(nrc,mps)))
1159 CALL gmpxy(7,real(kfl,mps),rej)
1161 IF (
jfd(kfl) > 0)
THEN
1162 c2ndf=
cfd(kfl)/real(
jfd(kfl),mps)
1163 CALL gmpxy(8,real(kfl,mps),c2ndf)
1164 andf=real(
dfd(kfl),mps)/real(
jfd(kfl),mps)
1165 CALL gmpxy(9,real(kfl,mps),andf)
1182 WRITE(*,*)
'Misalignment test wire chamber'
1185 CALL hmpdef( 9,-0.0015,+0.0015,
'True - fitted displacement')
1186 CALL hmpdef(10,-0.0015,+0.0015,
'True - fitted Vdrift')
1192 sums(1)=sums(1)+diff
1193 sums(2)=sums(2)+diff*diff
1195 sums(3)=sums(3)+diff
1196 sums(4)=sums(4)+diff*diff
1198 sums(1)=0.01_mpd*sums(1)
1199 sums(2)=sqrt(0.01_mpd*sums(2))
1200 sums(3)=0.01_mpd*sums(3)
1201 sums(4)=sqrt(0.01_mpd*sums(4))
1202 WRITE(*,143)
'Parameters 1 - 100: mean =',sums(1),
'rms =',sums(2)
1203 WRITE(*,143)
'Parameters 101 - 200: mean =',sums(3),
'rms =',sums(4)
1204143
FORMAT(6x,a28,f9.6,3x,a5,f9.6)
1207 WRITE(*,*)
' I label simulated fitted diff'
1208 WRITE(*,*)
' -------------------------------------------- '
1228 WRITE(*,*)
'Misalignment test Si tracker'
1231 CALL hmpdef( 9,-0.0025,+0.0025,
'True - fitted displacement X')
1232 CALL hmpdef(10,-0.025,+0.025,
'True - fitted displacement Y')
1243 sums(1)=sums(1)+1.0_mpd
1244 sums(2)=sums(2)+diff
1245 sums(3)=sums(3)+diff*diff
1250 err=sqrt(abs(gmati))
1252 sums(7)=sums(7)+1.0_mpd
1253 sums(8)=sums(8)+diff
1254 sums(9)=sums(9)+diff*diff
1257 IF (mod(i,3) == 1)
THEN
1261 sums(4)=sums(4)+1.0_mpd
1262 sums(5)=sums(5)+diff
1263 sums(6)=sums(6)+diff*diff
1268 err=sqrt(abs(gmati))
1270 sums(7)=sums(7)+1.0_mpd
1271 sums(8)=sums(8)+diff
1272 sums(9)=sums(9)+diff*diff
1277 sums(2)=sums(2)/sums(1)
1278 sums(3)=sqrt(sums(3)/sums(1))
1279 sums(5)=sums(5)/sums(4)
1280 sums(6)=sqrt(sums(6)/sums(4))
1281 WRITE(*,143)
'Parameters 1 - 500: mean =',sums(2),
'rms =',sums(3)
1282 WRITE(*,143)
'Parameters 501 - 700: mean =',sums(5),
'rms =',sums(6)
1283 IF (sums(7) > 0.5_mpd)
THEN
1284 sums(8)=sums(8)/sums(7)
1285 sums(9)=sqrt(sums(9)/sums(7))
1286 WRITE(*,143)
'Parameter pulls, all: mean =',sums(8),
'rms =',sums(9)
1290 WRITE(*,*)
' I label simulated fitted diff'
1291 WRITE(*,*)
' -------------------------------------------- '
1303 IF (mod(i,3) == 1)
THEN
1323 WRITE(8,*)
'Record',
nrec1,
' has largest residual:',
value1
1326 WRITE(8,*)
'Record',
nrec2,
' has largest Chi^2/Ndf:',
value2
1330 WRITE(8,*)
'Record',
nrec3,
' is first with error (rank deficit/NaN)'
1334 WRITE(8,*)
'In total 3 +',
nloopn,
' loops through the data files'
1336 WRITE(8,*)
'In total 2 +',
nloopn,
' loops through the data files'
1340 WRITE(8,*)
'In total ',
mnrsit,
' internal MINRES iterations'
1348 ntsec=nint(deltat,mpi)
1349 CALL sechms(deltat,nhour,minut,secnd)
1350 nsecnd=nint(secnd,mpi)
1351 WRITE(8,*)
'Total time =',ntsec,
' seconds =',nhour,
' h',minut, &
1352 ' m',nsecnd,
' seconds'
1354 WRITE(8,*)
'end ', chdate
1361 IF(sum(
nrejec) /= 0)
THEN
1363 WRITE(8,*)
'Data records rejected in last iteration: '
1370 WRITE(*,*)
' < Millepede II-P ending ... ', chdate
1377106
FORMAT(
' PARDISO dyn. memory allocation: ',f11.6,
' GB')
1387 WRITE(*,*)
'Postprocessing:'
1397102
FORMAT(2x,i4,i10,2x,3f10.5)
1398103
FORMAT(
' Times [in sec] for text processing',f12.3/ &
1401 ' func. value ',f12.3,
' *',f4.0/ &
1402 ' func. value, global matrix, solution',f12.3,
' *',f4.0/ &
1403 ' new solution',f12.3,
' *',f4.0/)
1404105
FORMAT(
' Peak dynamic memory allocation: ',f11.6,
' GB')
1424 INTEGER(mpi) :: istop
1425 INTEGER(mpi) :: itgbi
1426 INTEGER(mpi) :: itgbl
1428 INTEGER(mpi) :: itnlim
1429 INTEGER(mpi) :: nout
1431 INTEGER(mpi),
INTENT(IN) :: ivgbi
1468 globalcorrections, itnlim, nout, rtol, istop, itn, anorm, acond, rnorm, arnorm, ynorm)
1472 globalcorrections, itnlim, nout, rtol, istop, itn, anorm, acond, rnorm, arnorm, ynorm)
1475 globalcorrections, itnlim, nout, rtol, istop, itn, anorm, acond, rnorm, arnorm, ynorm)
1481 err=sqrt(abs(real(gmati,mps)))
1482 IF(gmati < 0.0_mpd) err=-err
1483 diag=matij(ivgbi,ivgbi)
1484 gcor2=real(1.0_mpd-1.0_mpd/(gmati*diag),mps)
1486101
FORMAT(1x,
' label parameter presigma differ', &
1487 ' Error gcor^2 iit'/ 1x,
'---------',2x,5(
'-----------'),2x,
'----')
1488102
FORMAT(i10,2x,4g12.4,f7.4,i6,i4)
1508 INTEGER(mpi) :: istop
1509 INTEGER(mpi) :: itgbi
1510 INTEGER(mpi) :: itgbl
1512 INTEGER(mpi) :: itnlim
1513 INTEGER(mpi) :: nout
1515 INTEGER(mpi),
INTENT(IN) :: ivgbi
1544 mxxnrm = real(
nagb,mpd)/sqrt(epsilon(mxxnrm))
1546 trcond = 1.0_mpd/epsilon(trcond)
1547 ELSE IF(
mrmode == 2)
THEN
1555 itnlim=itnlim, rtol=rtol, maxxnorm=mxxnrm, trancond=trcond, &
1559 itnlim=itnlim, rtol=rtol, maxxnorm=mxxnrm, trancond=trcond, &
1563 itnlim=itnlim, rtol=rtol, maxxnorm=mxxnrm, trancond=trcond, &
1570 err=sqrt(abs(real(gmati,mps)))
1571 IF(gmati < 0.0_mpd) err=-err
1572 diag=matij(ivgbi,ivgbi)
1573 gcor2=real(1.0_mpd-1.0_mpd/(gmati*diag),mps)
1575101
FORMAT(1x,
' label parameter presigma differ', &
1576 ' Error gcor^2 iit'/ 1x,
'---------',2x,5(
'-----------'),2x,
'----')
1577102
FORMAT(i10,2x,4g12.4,f7.4,i6,i4)
1590 INTEGER(mpi) :: icgb
1591 INTEGER(mpi) :: irhs
1592 INTEGER(mpi) :: itgbi
1593 INTEGER(mpi) :: ivgb
1595 INTEGER(mpi) :: jcgb
1597 INTEGER(mpi) :: label
1599 INTEGER(mpi) :: inone
1602 REAL(mpd) :: drhs(4)
1603 INTEGER(mpi) :: idrh (4)
1628 IF(abs(rhs) > climit)
THEN
1634 WRITE(*,101) (idrh(l),drhs(l),l=1,irhs)
1643 WRITE(*,101) (idrh(l),drhs(l),l=1,irhs)
1646 WRITE(*,102)
' Constraints: only equation values >', climit,
' are printed'
1647101
FORMAT(
' ',4(i6,g11.3))
1648102
FORMAT(a,g11.2,a)
1661 INTEGER(mpi) :: icgb
1662 INTEGER(mpi) :: icgrp
1663 INTEGER(mpi) :: ioff
1664 INTEGER(mpi) :: itgbi
1666 INTEGER(mpi) :: jcgb
1667 INTEGER(mpi) :: label
1668 INTEGER(mpi) :: labelf
1669 INTEGER(mpi) :: labell
1670 INTEGER(mpi) :: last
1671 INTEGER(mpi) :: line1
1672 INTEGER(mpi) :: ncon
1673 INTEGER(mpi) :: ndiff
1674 INTEGER(mpi) :: npar
1675 INTEGER(mpi) :: inone
1676 INTEGER(mpi) :: itype
1677 INTEGER(mpi) :: ncgbd
1678 INTEGER(mpi) :: ncgbr
1679 INTEGER(mpi) :: ncgbw
1680 INTEGER(mpi) :: ncgrpd
1681 INTEGER(mpi) :: ncgrpr
1682 INTEGER(mpi) :: next
1684 INTEGER(mpl):: length
1685 INTEGER(mpl) :: rows
1687 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: vecParConsOffsets
1688 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: vecParConsList
1689 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: vecConsParOffsets
1690 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: vecConsParList
1691 INTEGER(mpi),
DIMENSION(:,:),
ALLOCATABLE :: matConsGroupIndex
1704 IF(last < 0.AND.label < 0)
THEN
1707 IF(itype == 2) ncgbw=ncgbw+1
1714 IF(label > 0.AND.itype == 2)
THEN
1721 IF (ncgbw == 0)
THEN
1722 WRITE(*,*)
'GRPCON:',
ncgb,
' constraints found in steering files'
1724 WRITE(*,*)
'GRPCON:',
ncgb,
' constraints found in steering files,',ncgbw,
'weighted'
1729 length=
ncgb+1; rows=3
1742 CALL mpalloc(matconsgroupindex,rows,length,
'group index for constraint (I)')
1745 CALL mpalloc(vecconsparoffsets,length,
'offsets for global par list for cons. (I)')
1747 CALL mpalloc(vecparconsoffsets,length,
'offsets for cons. list for global par. (I)')
1748 vecparconsoffsets(1)=0
1750 vecparconsoffsets(i+1)=vecparconsoffsets(i)+
globalparcons(i)
1754 length=vecparconsoffsets(
ntgb+1)
1755 CALL mpalloc(vecconsparlist,length,
'global par. list for constraint (I)')
1756 CALL mpalloc(vecparconslist,length,
'constraint list for global par. (I)')
1761 vecconsparoffsets(1)=ioff
1773 vecparconslist(vecparconsoffsets(itgbi)+
globalparcons(itgbi))=icgb
1775 vecconsparlist(ioff+npar)=itgbi
1781 CALL sort1k(vecconsparlist(ioff+1),npar)
1785 next=vecconsparlist(ioff+j)
1786 IF (next /= last)
THEN
1788 vecconsparlist(ioff+ndiff) = next
1797 vecconsparoffsets(icgb+1)=ioff
1810 print *,
' Constraint #parameters first par. last par. first line'
1818 npar=vecconsparoffsets(icgb+1)-vecconsparoffsets(icgb)
1822 print *, jcgb, npar, labelf, labell, line1
1825 icgrp=matconsgroupindex(1,icgb)
1826 IF (icgrp == 0)
THEN
1828 DO i=vecconsparoffsets(icgb)+1, vecconsparoffsets(icgb+1)
1829 itgbi=vecconsparlist(i)
1831 DO j=vecparconsoffsets(itgbi)+1,vecparconsoffsets(itgbi+1)
1832 icgrp=matconsgroupindex(1,vecparconslist(j))
1838 IF (icgrp == 0)
THEN
1845 matconsgroupindex(2,icgb)=jcgb
1846 matconsgroupindex(3,icgb)=icgb
1847 DO i=vecconsparoffsets(icgb)+1, vecconsparoffsets(icgb+1)
1848 itgbi=vecconsparlist(i)
1851 DO j=vecparconsoffsets(itgbi)+1,vecparconsoffsets(itgbi+1)
1852 matconsgroupindex(1,vecparconslist(j))=icgrp
1856 WRITE(*,*)
'GRPCON:',
ncgrp,
' disjoint constraints groups built'
1864 icgb=matconsgroupindex(3,jcgb)
1869 icgrp=matconsgroupindex(1,jcgb)
1889 print *,
' cons.group first con. first par. last par. #cons #par'
1900 print *, icgrp,
matconsgroups(1,icgrp), labelf, labell, ncon, npar
1903 IF (ncon == npar)
THEN
1910 print *, icgrp,
matconsgroups(1,icgrp), labelf, labell,
' : cons.group resolved'
1925 print *, icgrp,
matconsgroups(1,icgrp), labelf, labell,
' : cons.group redundant'
1930 IF (ncgrpr > 0)
THEN
1931 WRITE(*,*)
'GRPCON:',ncgbr,
' redundancy constraints in ', ncgrpr,
' groups resolved'
1935 IF (ncgrpd > 0)
THEN
1936 WRITE(*,*)
'GRPCON:',ncgbd,
' redundancy constraints in ', ncgrpd,
' groups detected'
1959 INTEGER(mpi) :: icgb
1960 INTEGER(mpi) :: icgrp
1961 INTEGER(mpi) :: ifrst
1962 INTEGER(mpi) :: ilast
1963 INTEGER(mpi) :: isblck
1964 INTEGER(mpi) :: itgbi
1965 INTEGER(mpi) :: ivgb
1967 INTEGER(mpi) :: jcgb
1968 INTEGER(mpi) :: jfrst
1969 INTEGER(mpi) :: label
1970 INTEGER(mpi) :: labelf
1971 INTEGER(mpi) :: labell
1972 INTEGER(mpi) :: ncon
1973 INTEGER(mpi) :: ngrp
1974 INTEGER(mpi) :: npar
1975 INTEGER(mpi) :: ncnmxb
1976 INTEGER(mpi) :: ncnmxg
1977 INTEGER(mpi) :: nprmxb
1978 INTEGER(mpi) :: nprmxg
1979 INTEGER(mpi) :: inone
1980 INTEGER(mpi) :: nvar
1982 INTEGER(mpl):: length
1983 INTEGER(mpl) :: rows
1985 INTEGER(mpi),
DIMENSION(:,:),
ALLOCATABLE :: matConsGroupIndex
1998 length=
ncgrp+1; rows=3
2003 CALL mpalloc(matconsgroupindex,rows,length,
'group index for constraint (I)')
2040 IF (nvar > 0 .OR.
iskpec == 0)
THEN
2043 matconsgroupindex(1,
ncgb)=ngrp+1
2044 matconsgroupindex(2,
ncgb)=icgb
2045 matconsgroupindex(3,
ncgb)=nvar
2048 IF (
ncgb > ncon) ngrp=ngrp+1
2054 WRITE(*,*)
'PRPCON:',
ncgbe,
' empty constraints skipped'
2056 WRITE(*,*)
'PRPCON:',
ncgbe,
' empty constraints detected, to be fixed !!!'
2057 WRITE(*,*)
' (use option "skipemptycons" to skip those)'
2062 WRITE(*,*)
'!!! Switch to "-C" (checking input only), no calculation of a solution !!!'
2063 WRITE(8,*)
'!!! Switch to "-C" (checking input only), no calculation of a solution !!!'
2068 WRITE(*,*)
'PRPCON:',
ncgb,
' constraints accepted'
2071 IF(
ncgb == 0)
RETURN
2078 icgb=matconsgroupindex(2,jcgb)
2083 icgrp=matconsgroupindex(1,jcgb)
2109 WRITE(*,*)
' Cons. sorted index #var.par. first line first label last label'
2110 WRITE(*,*)
' Cons. group index first cons. last cons. first label last label'
2111 WRITE(*,*)
' Cons. block index first group last group first label last label'
2117 nvar=matconsgroupindex(3,jcgb)
2121 WRITE(*,*)
' Cons. sorted', jcgb, nvar, &
2124 WRITE(*,*)
' Cons. sorted', jcgb,
' empty (0)', &
2133 WRITE(*,*)
' Cons. group ', icgrp,
matconsgroups(1,icgrp), &
2146 DO i=icgrp,isblck,-1
2160 WRITE(*,*)
' Cons. block ',
ncblck, isblck, icgrp, labelf, labell
2178 ncnmxb=max(ncnmxb,ncon)
2179 nprmxb=max(nprmxb,npar)
2204 ncnmxg=max(ncnmxg,ncon)
2205 nprmxg=max(nprmxg,npar)
2240 WRITE(*,*)
'PRPCON: constraints split into ',
ncgrp,
'(disjoint) groups,'
2241 WRITE(*,*)
' groups combined into ',
ncblck,
'(non overlapping) blocks'
2242 WRITE(*,*)
' max group size (cons., par.) ', ncnmxg, nprmxg
2243 WRITE(*,*)
' max block size (cons., par.) ', ncnmxb, nprmxb
2261 INTEGER(mpi) :: icgb
2262 INTEGER(mpi) :: icgrp
2264 INTEGER(mpi) :: ifirst
2265 INTEGER(mpi) :: ilast
2266 INTEGER(mpl) :: ioffc
2267 INTEGER(mpl) :: ioffp
2268 INTEGER(mpi) :: irank
2269 INTEGER(mpi) :: ipar0
2270 INTEGER(mpi) :: itgbi
2271 INTEGER(mpi) :: ivgb
2273 INTEGER(mpi) :: jcgb
2275 INTEGER(mpi) :: label
2276 INTEGER(mpi) :: ncon
2277 INTEGER(mpi) :: npar
2278 INTEGER(mpi) :: nrank
2279 INTEGER(mpi) :: inone
2284 INTEGER(mpl):: length
2285 REAL(mpd),
DIMENSION(:),
ALLOCATABLE :: matConstraintsT
2286 REAL(mpd),
DIMENSION(:),
ALLOCATABLE :: auxVectorD
2287 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: auxVectorI
2291 IF(
ncgb == 0)
RETURN
2300 CALL mpalloc(auxvectori,length,
'auxiliary array (I)')
2301 CALL mpalloc(auxvectord,length,
'auxiliary array (D)')
2304 CALL mpalloc(matconstraintst,length,
'transposed matrix of constraints (blocks)')
2305 matconstraintst=0.0_mpd
2318 WRITE(*,*)
' Constraint group, #con, rank', icgrp, ncon, 0,
' (empty)'
2321 DO jcgb=ifirst,ilast
2333 IF(ivgb > 0) matconstraintst(int(jcgb-ifirst,mpl)*int(npar,mpl)+ivgb-ipar0+ioffc)= &
2334 matconstraintst(int(jcgb-ifirst,mpl)*int(npar,mpl)+ivgb-ipar0+ioffc)+factr
2341 DO ll=ioffc+1,ioffc+npar
2347 matconstraintst(int(i-1,mpl)*int(npar,mpl)+ll)* &
2348 matconstraintst(int(j-1,mpl)*int(npar,mpl)+ll)
2354 IF (
icheck > 1 .OR. irank < ncon)
THEN
2355 WRITE(*,*)
' Constraint group, #con, rank', icgrp, ncon, irank
2356 IF (irank < ncon)
THEN
2357 WRITE(*,*)
' .. rank deficit !! '
2358 WRITE(*,*)
' E.g. fix all parameters and remove all constraints related to label ', &
2363 ioffc=ioffc+int(npar,mpl)*int(ncon,mpl)
2370 WRITE(*,*)
'Rank of product matrix of constraints is',nrank, &
2371 ' for',
ncgb,
' constraint equations'
2372 WRITE(8,*)
'Rank of product matrix of constraints is',nrank, &
2373 ' for',
ncgb,
' constraint equations'
2374 IF(nrank <
ncgb)
THEN
2375 WRITE(*,*)
'Warning: insufficient constraint equations!'
2376 WRITE(8,*)
'Warning: insufficient constraint equations!'
2379 WRITE(*,*)
' --> enforcing SUBITO mode'
2380 WRITE(8,*)
' --> enforcing SUBITO mode'
2387 print *,
'QL decomposition of constraints matrix'
2390 WRITE(
lunlog,*)
'QL decomposition of constraints matrix'
2402 CALL lpqldec(matconstraintst,evmin,evmax)
2406 print *,
' largest |eigenvalue| of L: ', evmax
2407 print *,
' smallest |eigenvalue| of L: ', evmin
2408 IF (evmin == 0.0_mpd.AND.
icheck == 0)
THEN
2409 CALL peend(27,
'Aborted, singular QL decomposition of constraints matrix')
2410 stop
'FEASMA: stopping due to singular QL decomposition of constraints matrix'
2436 INTEGER(mpi) :: icgb
2437 INTEGER(mpi) :: icgrp
2438 INTEGER(mpi) :: iter
2439 INTEGER(mpi) :: itgbi
2440 INTEGER(mpi) :: ivgb
2441 INTEGER(mpi) :: ieblck
2442 INTEGER(mpi) :: isblck
2443 INTEGER(mpi) :: ifirst
2444 INTEGER(mpi) :: ilast
2446 INTEGER(mpi) :: jcgb
2447 INTEGER(mpi) :: label
2448 INTEGER(mpi) :: inone
2449 INTEGER(mpi) :: ncon
2451 REAL(mps),
INTENT(IN) :: concut
2452 INTEGER(mpi),
INTENT(OUT) :: iact
2459 REAL(mpd),
DIMENSION(:),
ALLOCATABLE :: vecCorrections
2463 IF(
ncgb == 0)
RETURN
2493 sum1=sqrt(sum1/real(
ncgb,mpd))
2494 sum2=sum2/real(
ncgb,mpd)
2496 IF(iter == 1.AND.sum1 < concut)
RETURN
2498 IF(iter == 1.AND.
ncgb <= 12)
THEN
2500 WRITE(*,*)
'Constraint equation discrepancies:'
2502101
FORMAT(4x,4(i5,g12.4))
2504103
FORMAT(10x,
' Cut on rms value is',g8.1)
2509 WRITE(*,*)
'Improve constraints'
2513 WRITE(*,102) iter,sum1,sum2,sum3
2514102
FORMAT(i6,
' rms',g12.4,
' avrg_abs',g12.4,
' max_abs',g12.4)
2516 CALL mpalloc(veccorrections,int(
nvgb,mpl),
'constraint corrections')
2517 veccorrections=0.0_mpd
2525 ieblck=isblck+(ncon*(ncon+1))/2
2594 INTEGER(mpi) :: iact
2595 INTEGER(mpi) :: ierrc
2596 INTEGER(mpi) :: ierrf
2597 INTEGER(mpi) :: ioffp
2599 INTEGER(mpi) :: ithr
2600 INTEGER(mpi) :: jfile
2601 INTEGER(mpi) :: jrec
2603 INTEGER(mpi) :: kfile
2606 INTEGER(mpi) :: mpri
2608 INTEGER(mpi) :: nact
2609 INTEGER(mpi) :: nbuf
2610 INTEGER(mpi) :: ndata
2611 INTEGER(mpi) :: noff
2612 INTEGER(mpi) :: noffs
2613 INTEGER(mpi) :: npointer
2614 INTEGER(mpi) :: npri
2618 INTEGER(mpi) :: nrpr
2619 INTEGER(mpi) :: nthr
2620 INTEGER(mpi) :: ntot
2621 INTEGER(mpi) :: maxRecordSize
2622 INTEGER(mpi) :: maxRecordFile
2624 INTEGER(mpi),
INTENT(OUT) :: more
2634 CHARACTER (LEN=7) :: cfile
2639 SUBROUTINE readc(bufferD, bufferF, bufferI, bufferLength, lun, err)
BIND(c)
2641 REAL(c_double),
DIMENSION(*),
INTENT(OUT) :: bufferD
2642 REAL(c_float),
DIMENSION(*),
INTENT(OUT) :: bufferF
2643 INTEGER(c_int),
DIMENSION(*),
INTENT(OUT) :: bufferI
2644 INTEGER(c_int),
INTENT(INOUT) :: bufferLength
2645 INTEGER(c_int),
INTENT(IN),
VALUE :: lun
2646 INTEGER(c_int),
INTENT(OUT) :: err
2647 END SUBROUTINE readc
2653 DATA npri / 0 /, mpri / 1000 /
2701 files:
DO WHILE (jfile > 0)
2705 CALL binopn(kfile,ithr,ios)
2711 IF(kfile <=
nfilf)
THEN
2732 eof=(ierrc <= 0.AND.ierrc /= -4)
2733 IF(eof.AND.ierrc < 0)
THEN
2734 WRITE(*,*)
'Read error for binary Cfile', kfile,
'record', jrec+1,
':', ierrc
2735 WRITE(8,*)
'Read error for binary Cfile', kfile,
'record', jrec+1,
':', ierrc
2737 WRITE(cfile,
'(I7)') kfile
2738 CALL peend(18,
'Aborted, read error(s) for binary file ' // cfile)
2739 stop
'PEREAD: stopping due to read errors'
2741 IF (
kfd(1,jfile) == 1)
THEN
2747 IF(eof)
EXIT records
2752 xfd(jfile)=max(
xfd(jfile),n)
2769 IF ((noff+nr+2+
ndimbuf >= ndata*(iact+1)).OR.(nbuf >= npointer))
EXIT files
2784 IF (
kfd(1,jfile) == 1)
THEN
2785 print *,
'PEREAD: file ', kfile,
'read the first time, found',jrec,
' records'
2789 IF (-
kfd(1,jfile) /= jrec)
THEN
2790 WRITE(cfile,
'(I7)') kfile
2791 CALL peend(19,
'Aborted, binary file modified (length) ' // cfile)
2792 stop
'PEREAD: file modified (length)'
2843 IF (
kfd(1,jfile) == 1)
kfd(1,jfile)=-nrc
2862 DO WHILE (
nloopn == 0.AND.ntot >= nrpr)
2863 WRITE(*,*)
' Record ',nrpr
2864 IF (nrpr < 100000)
THEN
2873 IF (npri == 1)
WRITE(*,100)
2875100
FORMAT(/
' PeRead records active file' &
2876 /
' total block threads number')
2877101
FORMAT(
' PeRead',4i10)
2890 IF (
xfd(k) > maxrecordsize)
THEN
2891 maxrecordsize=
xfd(k)
2894 dw=real(-
kfd(1,k),mpd)
2897 ds1=ds1+dw*real(
wfd(k),mpd)
2898 ds2=ds2+dw*real(
wfd(k)**2,mpd)
2900 print *,
'PEREAD: file ', maxrecordfile,
'with max record size ', maxrecordsize
2901 IF (
nfilw > 0.AND.ds0 > 0.0_mpd)
THEN
2905 WRITE(lun,177)
nfilw,real(ds1,mps),real(ds2,mps)
2906177
FORMAT(/
' !!!!!',i4,
' weighted binary files', &
2907 /
' !!!!! mean, variance of weights =',2g12.4)
2925179
FORMAT(/
' Read cache usage (#blocks, #records, ', &
2926 'min,max records/block'/17x,i10,i12,2i10)
2944 INTEGER(mpi),
INTENT(IN) :: mode
2946 INTEGER(mpi) :: ibuf
2947 INTEGER(mpi) :: ichunk
2949 INTEGER(mpi) :: itgbi
2955 INTEGER(mpi),
PARAMETER :: maxbad = 100
2956 INTEGER(mpi) :: nbad
2957 INTEGER(mpi) :: nerr
2958 INTEGER(mpi) :: inone
2976 CALL isjajb(nst,ist,ja,jb,jsp)
2995#ifdef SCOREP_USER_ENABLE
2996 scorep_user_region_by_name_begin(
"UR_peprep", scorep_user_region_type_common)
3001 CALL pechk(ibuf,nerr)
3004 IF(nbad >= maxbad)
EXIT
3009 CALL isjajb(nst,ist,ja,jb,jsp)
3022 CALL peend(20,
'Aborted, bad binary records')
3023 stop
'PEREAD: stopping due to bad records'
3026#ifdef SCOREP_USER_ENABLE
3027 scorep_user_region_by_name_end(
"UR_peprep")
3047 INTEGER(mpi) :: ioff
3054 INTEGER(mpi),
INTENT(IN) :: ibuf
3055 INTEGER(mpi),
INTENT(OUT) :: nerr
3064 outer:
DO WHILE(is < nst)
3069 IF(is > nst)
EXIT outer
3075 IF(is > nst)
EXIT outer
3092100
FORMAT(
' PEREAD: record ', i8,
' in file ',i6,
' is broken !!!')
3102101
FORMAT(
' PEREAD: record ', i8,
' in file ',i6,
' contains ', i6,
' NaNs !!!')
3118 INTEGER(mpi) :: ibuf
3119 INTEGER(mpi) :: ichunk
3120 INTEGER(mpi) :: iproc
3121 INTEGER(mpi) :: ioff
3122 INTEGER(mpi) :: ioffbi
3124 INTEGER(mpi) :: itgbi
3129 INTEGER(mpi) :: nalg
3130 INTEGER(mpi) :: neqna
3133 INTEGER(mpi) :: nzero
3134 INTEGER(mpi) :: inone
3135 INTEGER(mpl) :: length
3170 CALL isjajb(nst,ist,ja,jb,jsp)
3201 CALL isjajb(nst,ist,ja,jb,jsp)
3229#ifdef SCOREP_USER_ENABLE
3230 scorep_user_region_by_name_begin(
"UR_pepgrp", scorep_user_region_type_common)
3239 CALL pargrp(ist+1,ist+nnz)
3251#ifdef SCOREP_USER_ENABLE
3252 scorep_user_region_by_name_end(
"UR_pepgrp")
3271 INTEGER(mpi) :: istart
3272 INTEGER(mpi) :: itgbi
3274 INTEGER(mpi) :: jstart
3275 INTEGER(mpi) :: jtgbi
3276 INTEGER(mpi) :: lstart
3277 INTEGER(mpi) :: ltgbi
3279 INTEGER(mpi),
INTENT(IN) :: inds
3280 INTEGER(mpi),
INTENT(IN) :: inde
3282 IF (inds > inde)
RETURN
3291 IF (istart == 0)
THEN
3292 IF (itgbi /= ltgbi+1)
THEN
3295 IF (lstart == 0)
THEN
3310 IF (istart /= jstart)
THEN
3321 IF (itgbi /= ltgbi+1)
THEN
3336 IF (istart /= jstart)
THEN
3388 INTEGER(mpi),
INTENT(IN) :: nst
3389 INTEGER(mpi),
INTENT(IN OUT) :: is
3390 INTEGER(mpi),
INTENT(OUT) :: ja
3391 INTEGER(mpi),
INTENT(OUT) :: jb
3392 INTEGER(mpi),
INTENT(OUT) :: jsp
3400 IF(is >= nst)
RETURN
3454 INTEGER(mpi) :: ioffb
3456 INTEGER(mpi) :: itgbi
3457 INTEGER(mpi) :: itgbij
3458 INTEGER(mpi) :: itgbik
3459 INTEGER(mpi) :: ivgb
3460 INTEGER(mpi) :: ivgbij
3461 INTEGER(mpi) :: ivgbik
3464 INTEGER(mpi) :: lastit
3466 INTEGER(mpi) :: ncrit
3467 INTEGER(mpi) :: ngras
3468 INTEGER(mpi) :: nparl
3470 INTEGER(mpl) :: nrej
3471 INTEGER(mpi) :: inone
3472 INTEGER(mpi) :: ilow
3473 INTEGER(mpi) :: nlow
3474 INTEGER(mpi) :: nzero
3494 CALL gmpdef(1,4,
'Function value in iterations')
3496 CALL gmpdef(2,3,
'Number of MINRES steps vs iteration nr')
3498 CALL hmpdef( 5,0.0,0.0,
'Number of degrees of freedom')
3499 CALL hmpdef(11,0.0,0.0,
'Number of local parameters')
3500 CALL hmpdef(16,0.0,24.0,
'LOG10(cond(band part decomp.)) local fit ')
3501 CALL hmpdef(23,0.0,0.0,
'SQRT of diagonal elements without presigma')
3502 CALL hmpdef(24,0.0,0.0,
'Log10 of off-diagonal elements')
3503 CALL hmpdef(25,0.0,0.0,
'Relative individual pre-sigma')
3504 CALL hmpdef(26,0.0,0.0,
'Relative global pre-sigma')
3509 'Normalized residuals of single (global) measurement')
3511 'Normalized residuals of single (local) measurement')
3513 'Pulls of single (global) measurement')
3515 'Pulls of single (local) measurement')
3516 CALL hmpdef(4,0.0,0.0,
'Chi^2/Ndf after local fit')
3517 CALL gmpdef(4,5,
'location, dispersion (res.) vs record nr')
3518 CALL gmpdef(5,5,
'location, dispersion (pull) vs record nr')
3532 IF(
nloopn == 2)
CALL hmpdef(6,0.0,0.0,
'Down-weight fraction')
3535 IF(
iterat /= lastit)
THEN
3543 ELSE IF(
iterat >= 1)
THEN
3591 WRITE(*,111) nparl,ncrit,usei,used,peaki,peakd
3592111
FORMAT(
' Write cache usage (#flush,#overrun,<levels>,', &
3593 'peak(levels))'/2i7,
',',4(f6.1,
'%'))
3601 ELSE IF (
mbandw > 0)
THEN
3633 print *,
" ... warning ..."
3634 print *,
" global parameters with too few (< MREQENA) accepted entries: ", nlow
3638 IF(
icalcm == 1 .AND. nzero > 0)
THEN
3640 WRITE(*,*)
'Warning: the rank defect of the symmetric',
nfgb, &
3641 '-by-',
nfgb,
' matrix is ',
ndefec,
' (should be zero).'
3642 WRITE(lun,*)
'Warning: the rank defect of the symmetric',
nfgb, &
3643 '-by-',
nfgb,
' matrix is ',
ndefec,
' (should be zero).'
3646 WRITE(*,*)
' --> enforcing SUBITO mode'
3647 WRITE(lun,*)
' --> enforcing SUBITO mode'
3657 elmt=real(matij(i,i),mps)
3658 IF(elmt > 0.0)
CALL hmpent(23,1.0/sqrt(elmt))
3692 CALL addsums(1, adder, 0, 1.0_mpl)
3706 IF(sgm > 0.0) weight=1.0/sgm**2
3722 IF(itgbij /= 0)
THEN
3728 IF (factrj == 0.0_mpd) cycle
3738 IF(
icalcm == 1.AND.ivgbij > 0)
THEN
3745 IF(ivgbij > 0.AND.ivgbik > 0)
THEN
3746 CALL mupdat(ivgbij,ivgbik,weight*factrj*factrk)
3752 adder=weight*dsum**2
3753 CALL addsums(1, adder, 1, 1.0_mpl)
3771 ratae=max(-99.9,ratae)
3780 WRITE(*,*)
'Data records rejected in initial loop:'
3825 WRITE(lun,*)
' === local fits have bordered band matrix structure ==='
3826 IF (
nbndr(1) > 0 )
WRITE(lun,101)
' NBNDR',
nbndr(1),
'number of records (upper/left border)'
3827 IF (
nbndr(2) > 0 )
WRITE(lun,101)
' NBNDR',
nbndr(2),
'number of records (lower/right border)'
3828 WRITE(lun,101)
' NBDRX',
nbdrx,
'max border size'
3829 WRITE(lun,101)
' NBNDX',
nbndx,
'max band width'
3838101
FORMAT(1x,a8,
' =',i14,
' = ',a)
3853 INTEGER(mpi),
INTENT(IN) :: lunp
3858101
FORMAT(
' it fc',
' fcn_value dfcn_exp slpr costh iit st', &
3859 ' ls step cutf',1x,
'rejects hhmmss FMS')
3860102
FORMAT(
' -- --',
' ----------- -------- ---- ----- --- --', &
3861 ' -- ----- ----',1x,
'------- ------ ---')
3877 INTEGER(mpl) :: nrej
3878 INTEGER(mpi) :: nsecnd
3882 REAL(mps) :: slopes(3)
3883 REAL(mps) :: steps(3)
3884 REAL,
DIMENSION(2) :: ta
3887 INTEGER(mpi),
INTENT(IN) :: lunp
3889 CHARACTER (LEN=4):: ccalcm(4)
3890 DATA ccalcm /
' end',
' S',
' F ',
' FMS' /
3894 IF(nrej > 9999999) nrej=9999999
3898 nsecnd=nint(secnd,mpi)
3907 CALL ptlopt(nfa,ma,slopes,steps)
3908 ratae=max(-99.9,min(99.9,slopes(2)/slopes(1)))
3914103
FORMAT(i3,i3,e12.5,38x,f5.1, 1x,i7, i3,i2.2,i2.2,a4)
3915104
FORMAT(i3,i3,e12.5,1x,e8.2,f6.3,f6.3,i5,2i3,f6.3,f5.1, &
3916 1x,i7, i3,i2.2,i2.2,a4)
3917105
FORMAT(i3,i3,e12.5,1x,e8.2,12x,i5,i3,9x,f5.1, &
3918 1x,i7, i3,i2.2,i2.2,a4)
3931 INTEGER(mpi) :: minut
3933 INTEGER(mpi) :: nhour
3934 INTEGER(mpl) :: nrej
3935 INTEGER(mpi) :: nsecnd
3939 REAL(mps) :: slopes(3)
3940 REAL(mps) :: steps(3)
3941 REAL,
DIMENSION(2) :: ta
3944 INTEGER(mpi),
INTENT(IN) :: lunp
3945 CHARACTER (LEN=4):: ccalcm(4)
3946 DATA ccalcm /
' end',
' S',
' F ',
' FMS' /
3950 IF(nrej > 9999999) nrej=9999999
3954 nsecnd=nint(secnd,mpi)
3958 CALL ptlopt(nfa,ma,slopes,steps)
3959 ratae=abs(slopes(2)/slopes(1))
3964104
FORMAT(3x,i3,e12.5,9x, 35x, i7, i3,i2.2,i2.2,a4)
3965105
FORMAT(3x,i3,e12.5,9x, f6.3,14x,i3,f6.3,6x, i7, i3,i2.2,i2.2,a4)
3979 INTEGER(mpi) :: nsecnd
3982 REAL,
DIMENSION(2) :: ta
3985 INTEGER(mpi),
INTENT(IN) :: lunp
3986 CHARACTER (LEN=4):: ccalcm(4)
3987 DATA ccalcm /
' end',
' S',
' F ',
' FMS' /
3992 nsecnd=nint(secnd,mpi)
3994 WRITE(lunp,106) nhour,minut,nsecnd,ccalcm(
lcalcm)
3995106
FORMAT(69x,i3,i2.2,i2.2,a4)
4005 INTEGER(mpi) :: lunit
4007 WRITE(lunit,102)
'Explanation of iteration table'
4008 WRITE(lunit,102)
'=============================='
4009 WRITE(lunit,101)
'it', &
4010 'iteration number. Global parameters are improved for it > 0.'
4011 WRITE(lunit,102)
'First function evaluation is called iteraton 0.'
4012 WRITE(lunit,101)
'fc',
'number of function evaluations.'
4013 WRITE(lunit,101)
'fcn_value',
'value of 2 x Likelihood function (LF).'
4014 WRITE(lunit,102)
'The final value is the chi^2 value of the fit and should'
4015 WRITE(lunit,102)
'be about equal to the NDF (see below).'
4016 WRITE(lunit,101)
'dfcn_exp', &
4017 'expected reduction of the value of the Likelihood function (LF)'
4018 WRITE(lunit,101)
'slpr',
'ratio of the actual slope to inital slope.'
4019 WRITE(lunit,101)
'costh', &
4020 'cosine of angle between search direction and -gradient'
4022 WRITE(lunit,101)
'iit', &
4023 'number of internal iterations in MINRES algorithm'
4024 WRITE(lunit,101)
'st',
'stop code of MINRES algorithm'
4025 WRITE(lunit,102)
'< 0: rhs is very special, with beta2 = 0'
4026 WRITE(lunit,102)
'= 0: rhs b = 0, i.e. the exact solution is x = 0'
4027 WRITE(lunit,102)
'= 1 requested accuracy achieved, as determined by rtol'
4028 WRITE(lunit,102)
'= 2 reasonable accuracy achieved, given eps'
4029 WRITE(lunit,102)
'= 3 x has converged to an eigenvector'
4030 WRITE(lunit,102)
'= 4 matrix ill-conditioned (Acond has exceeded 0.1/eps)'
4031 WRITE(lunit,102)
'= 5 the iteration limit was reached'
4032 WRITE(lunit,102)
'= 6 Matrix x vector does not define a symmetric matrix'
4033 WRITE(lunit,102)
'= 7 Preconditioner does not define a symmetric matrix'
4034 ELSEIF (
metsol == 5)
THEN
4035 WRITE(lunit,101)
'iit', &
4036 'number of internal iterations in MINRES-QLP algorithm'
4037 WRITE(lunit,101)
'st',
'stop code of MINRES-QLP algorithm'
4038 WRITE(lunit,102)
'= 1: beta_{k+1} < eps, iteration k is the final Lanczos step.'
4039 WRITE(lunit,102)
'= 2: beta2 = 0. If M = I, b and x are eigenvectors of A.'
4040 WRITE(lunit,102)
'= 3: beta1 = 0. The exact solution is x = 0.'
4041 WRITE(lunit,102)
'= 4: A solution to (poss. singular) Ax = b found, given rtol.'
4042 WRITE(lunit,102)
'= 5: A solution to (poss. singular) Ax = b found, given eps.'
4043 WRITE(lunit,102)
'= 6: Pseudoinverse solution for singular LS problem, given rtol.'
4044 WRITE(lunit,102)
'= 7: Pseudoinverse solution for singular LS problem, given eps.'
4045 WRITE(lunit,102)
'= 8: The iteration limit was reached.'
4046 WRITE(lunit,102)
'= 9: The operator defined by Aprod appears to be unsymmetric.'
4047 WRITE(lunit,102)
'=10: The operator defined by Msolve appears to be unsymmetric.'
4048 WRITE(lunit,102)
'=11: The operator defined by Msolve appears to be indefinite.'
4049 WRITE(lunit,102)
'=12: xnorm has exceeded maxxnorm or will exceed it next iteration.'
4050 WRITE(lunit,102)
'=13: Acond has exceeded Acondlim or 0.1/eps.'
4051 WRITE(lunit,102)
'=14: Least-squares problem but no converged solution yet.'
4052 WRITE(lunit,102)
'=15: A null vector obtained, given rtol.'
4054 WRITE(lunit,101)
'ls',
'line search info'
4055 WRITE(lunit,102)
'< 0 recalculate function'
4056 WRITE(lunit,102)
'= 0: N or STP lt 0 or step not descending'
4057 WRITE(lunit,102)
'= 1: Linesearch convergence conditions reached'
4058 WRITE(lunit,102)
'= 2: interval of uncertainty at lower limit'
4059 WRITE(lunit,102)
'= 3: max nr of line search calls reached'
4060 WRITE(lunit,102)
'= 4: step at the lower bound'
4061 WRITE(lunit,102)
'= 5: step at the upper bound'
4062 WRITE(lunit,102)
'= 6: rounding error limitation'
4063 WRITE(lunit,101)
'step', &
4064 'the factor for the Newton step during the line search. Usually'
4066 'a value of 1 gives a sufficient reduction of the LF. Oherwise'
4067 WRITE(lunit,102)
'other step values are tried.'
4068 WRITE(lunit,101)
'cutf', &
4069 'cut factor. Local fits are rejected, if their chi^2 value'
4071 'is larger than the 3-sigma chi^2 value times the cut factor.'
4072 WRITE(lunit,102)
'A cut factor of 1 is used finally, but initially a larger'
4073 WRITE(lunit,102)
'factor may be used. A value of 0.0 means no cut.'
4074 WRITE(lunit,101)
'rejects',
'total number of rejected local fits.'
4075 WRITE(lunit,101)
'hmmsec',
'the time in hours (h), minutes (mm) and seconds.'
4076 WRITE(lunit,101)
'FMS',
'calculation of Function value, Matrix, Solution.'
4079101
FORMAT(a9,
' = ',a)
4096 INTEGER(mpi),
INTENT(IN) :: i
4097 INTEGER(mpi),
INTENT(IN) :: j
4098 REAL(mpd),
INTENT(IN) :: add
4100 INTEGER(mpl):: ijadd
4101 INTEGER(mpl):: ijcsr3
4106 IF(i <= 0.OR.j <= 0.OR. add == 0.0_mpd)
RETURN
4127 ELSE IF(
matsto == 2)
THEN
4155 CALL peend(23,
'Aborted, bad matrix index')
4156 stop
'mupdat: bad index'
4181 INTEGER(mpi),
INTENT(IN) :: i
4182 INTEGER(mpi),
INTENT(IN) :: j1
4183 INTEGER(mpi),
INTENT(IN) :: j2
4184 INTEGER(mpi),
INTENT(IN) :: il
4185 INTEGER(mpi),
INTENT(IN) :: jl
4186 INTEGER(mpi),
INTENT(IN) :: n
4187 REAL(mpd),
INTENT(IN) :: sub((n*n+n)/2)
4194 INTEGER(mpi):: iblast
4195 INTEGER(mpi):: iblock
4201 INTEGER(mpi):: jblast
4202 INTEGER(mpi):: jblock
4212 IF(i <= 0.OR.j1 <= 0.OR.j2 > i)
RETURN
4226 print *,
' MGUPDT: ij=0', i,j1,j2,il,jl,ij,lr,iprc,
matsto
4233 ijl=(k*k-k)/2+jl+ir-ia1
4250 ijl=(k*k-k)/2+jl+ir-ia1
4254 IF (jblock /= jblast.OR.iblock /= iblast)
THEN
4278 CALL ijpgrp(i,j1,ij,lr,iprc)
4280 print *,
' MGUPDT: ij=0', i,j1,j2,il,jl,ij,lr,iprc,
matsto
4344SUBROUTINE loopbf(nrej,numfil,naccf,chi2f,ndff)
4371 INTEGER(mpi) :: ibuf
4372 INTEGER(mpi) :: ichunk
4373 INTEGER(mpl) :: icmn
4374 INTEGER(mpl) :: icost
4376 INTEGER(mpi) :: idiag
4378 INTEGER(mpi) :: iext
4386 INTEGER(mpi) :: imeas
4389 INTEGER(mpi) :: ioffb
4390 INTEGER(mpi) :: ioffc
4391 INTEGER(mpi) :: ioffd
4392 INTEGER(mpi) :: ioffe
4393 INTEGER(mpi) :: ioffi
4394 INTEGER(mpi) :: ioffq
4395 INTEGER(mpi) :: iprc
4396 INTEGER(mpi) :: iprcnx
4397 INTEGER(mpi) :: iprdbg
4398 INTEGER(mpi) :: iproc
4399 INTEGER(mpi) :: irbin
4400 INTEGER(mpi) :: isize
4402 INTEGER(mpi) :: iter
4403 INTEGER(mpi) :: itgbi
4404 INTEGER(mpi) :: ivgbj
4405 INTEGER(mpi) :: ivgbk
4406 INTEGER(mpi) :: ivpgrp
4416 INTEGER(mpi) :: joffd
4417 INTEGER(mpi) :: joffi
4418 INTEGER(mpi) :: jproc
4422 INTEGER(mpi) :: kbdr
4423 INTEGER(mpi) :: kbdrx
4424 INTEGER(mpi) :: kbnd
4427 INTEGER(mpi) :: lvpgrp
4428 INTEGER(mpi) :: mbdr
4429 INTEGER(mpi) :: mbnd
4430 INTEGER(mpi) :: mside
4431 INTEGER(mpi) :: nalc
4432 INTEGER(mpi) :: nalg
4436 INTEGER(mpi) :: ndown
4438 INTEGER(mpi) :: nfred
4439 INTEGER(mpi) :: nfrei
4441 INTEGER(mpi) :: nprdbg
4442 INTEGER(mpi) :: nrank
4445 INTEGER(mpi) :: nter
4446 INTEGER(mpi) :: nweig
4447 INTEGER(mpi) :: ngrp
4448 INTEGER(mpi) :: npar
4450 INTEGER(mpl),
INTENT(IN OUT) :: nrej(6)
4451 INTEGER(mpi),
INTENT(IN) :: numfil
4452 INTEGER(mpi),
INTENT(IN OUT) :: naccf(numfil)
4453 REAL(mps),
INTENT(IN OUT) :: chi2f(numfil)
4454 INTEGER(mpi),
INTENT(IN OUT) :: ndff(numfil)
4464 INTEGER(mpi) :: ijprec
4471 CHARACTER (LEN=3):: chast
4472 DATA chuber/1.345_mpd/
4473 DATA cauchy/2.3849_mpd/
4521 IF(
nloopn == 1.AND.mod(nrc,100000_mpl) == 0)
THEN
4522 WRITE(*,*)
'Record',nrc,
' ... still reading'
4523 IF(
monpg1>0)
WRITE(
lunlog,*)
'Record',nrc,
' ... still reading'
4533 IF(nrc ==
nrecpr) lprnt=.true.
4534 IF(nrc ==
nrecp2) lprnt=.true.
4535 IF(nrc ==
nrecer) lprnt=.true.
4540 IF (nprdbg == 1) iprdbg=iproc
4541 IF (iproc /= iprdbg) lprnt=.false.
4546 WRITE(1,*)
'------------------ Loop',
nloopn, &
4547 ': Printout for record',nrc,iproc
4558 CALL isjajb(nst,ist,ja,jb,jsp)
4560 IF(imeas == 0)
WRITE(1,1121)
45651121
FORMAT(/
'Measured value and local derivatives'/ &
4566 ' i measured std_dev index...derivative ...')
45671122
FORMAT(i3,2g12.4,3(i3,g12.4)/(27x,3(i3,g12.4)))
4573 CALL isjajb(nst,ist,ja,jb,jsp)
4575 IF(imeas == 0)
WRITE(1,1123)
45871123
FORMAT(/
'Global derivatives'/ &
4588 ' i label gindex vindex derivative ...')
45891124
FORMAT(i3,2(i9,i7,i7,g12.4)/(3x,2(i9,i7,i7,g12.4)))
45901125
FORMAT(i3,2(i9,i7,i7,g12.4))
4600 WRITE(1,*)
'Data corrections using values of global parameters'
4601 WRITE(1,*)
'=================================================='
4611 CALL isjajb(nst,ist,ja,jb,jsp)
4643101
FORMAT(
' index measvalue corrvalue sigma')
4644102
FORMAT(i6,2x,2g12.4,
' +-',g12.4)
4646 IF(nalc <= 0)
GO TO 90
4648 ngg=(nalg*nalg+nalg)/2
4661 IF (ivpgrp /= lvpgrp)
THEN
4669 IF (npar /= nalg)
THEN
4670 print *,
' mismatch of number of global parameters ', nrc, nalg, npar, ngrp
4680 CALL peend(35,
'Aborted, mismatch of number of global parameters')
4681 stop
' mismatch of number of global parameters '
4708 DO WHILE(iter < nter)
4713 WRITE(1,*)
'Outlier-suppression iteration',iter,
' of',nter
4714 WRITE(1,*)
'=========================================='
4724 DO i=1,(nalc*nalc+nalc)/2
4735 wght =1.0_mpd/rerr**2
4740 IF(abs(resid) > chuber*rerr)
THEN
4741 wght=wght*chuber*rerr/abs(resid)
4745 wght=wght/(1.0+(resid/rerr/cauchy)**2)
4749 IF(lprnt.AND.iter /= 1.AND.nter /= 1)
THEN
4751 IF(abs(resid) > chuber*rerr) chast=
'* '
4752 IF(abs(resid) > 3.0*rerr) chast=
'** '
4753 IF(abs(resid) > 6.0*rerr) chast=
'***'
4754 IF(imeas == 0)
WRITE(1,*)
'Second loop: accumulate'
4755 IF(imeas == 0)
WRITE(1,103)
4760 WRITE(1,104) imeas,rmeas,resid,rerr,r1,chast,r2
4762103
FORMAT(
' index corrvalue residuum sigma', &
4764104
FORMAT(i6,2x,2g12.4,
' +-',g12.4,f7.2,1x,a3,f8.2)
4779 im=min(nalc+1-ij,nalc+1-ik)
4786 IF (iter == 1.AND.nalc > 5.AND.
lfitbb > 0)
THEN
4789 icmn=int(nalc,mpl)**3
4795 icost=6*int(nalc-kbdr,mpl)*int(kbnd+kbdr+1,mpl)**2+2*int(kbdr,mpl)**3
4796 IF (icost < icmn)
THEN
4808 kbdr=max(kbdr,
ibandh(k+nalc))
4809 icost=6*int(nalc-kbdr,mpl)*int(kbnd+kbdr+1,mpl)**2+2*int(kbdr,mpl)**3
4810 IF (icost < icmn)
THEN
4834 IF (
icalcm == 1.OR.lprnt) inv=2
4835 IF (mside == 1)
THEN
4843 IF (evdmin > 0.0_mpl) cndl10=log10(real(evdmax/evdmin,mps))
4853 IF ((.NOT.(
blvec(k) <= 0.0_mpd)).AND. (.NOT.(
blvec(k) > 0.0_mpd))) nan=nan+1
4858 WRITE(1,*)
'Parameter determination:',nalc,
' parameters,',
' rank=',nrank
4859 WRITE(1,*)
'-----------------------'
4860 IF(ndown /= 0)
WRITE(1,*)
' ',ndown,
' data down-weighted'
4876 wght =1.0_mpd/rerr**2
4900 IF (0.999999_mpd/wght > dvar)
THEN
4901 pull=rmeas/sqrt(1.0_mpd/wght-dvar)
4904 CALL hmpent(13,real(pull,mps))
4905 CALL gmpms(5,rec,real(pull,mps))
4907 CALL hmpent(14,real(pull,mps))
4926 IF(iter == 1.AND.jb < ist.AND.lhist) &
4927 CALL gmpms(4,rec,real(rmeas/rerr,mps))
4929 dchi2=wght*rmeas*rmeas
4934 IF(abs(resid) > chuber*rerr)
THEN
4935 wght=wght*chuber*rerr/abs(resid)
4936 dchi2=2.0*chuber*(abs(resid)/rerr-0.5*chuber)
4939 wght=wght/(1.0_mpd+(resid/rerr/cauchy)**2)
4940 dchi2=log(1.0_mpd+(resid/rerr/cauchy)**2)*cauchy**2
4950 IF(abs(resid) > chuber*rerr) chast=
'* '
4951 IF(abs(resid) > 3.0*rerr) chast=
'** '
4952 IF(abs(resid) > 6.0*rerr) chast=
'***'
4953 IF(imeas == 0)
WRITE(1,*)
'Third loop: single residuals'
4954 IF(imeas == 0)
WRITE(1,105)
4958 IF(resid < 0.0) r1=-r1
4959 IF(resid < 0.0) r2=-r2
4962105
FORMAT(
' index corrvalue residuum sigma', &
4964106
FORMAT(i6,2x,2g12.4,
' +-',g12.4,f7.2,1x,a3,f8.2)
4966 IF(iter == nter)
THEN
4968 resmax=max(resmax,abs(rmeas)/rerr)
4971 IF(iter == 1.AND.lhist)
THEN
4973 CALL hmpent( 3,real(rmeas/rerr,mps))
4975 CALL hmpent(12,real(rmeas/rerr,mps))
4982 resing=(real(nweig,mps)-real(suwt,mps))/real(nweig,mps)
4984 IF(iter == 1)
CALL hmpent( 5,real(ndf,mps))
4985 IF(iter == 1)
CALL hmpent(11,real(nalc,mps))
4990 WRITE(1,*)
'Chi^2=',summ,
' at',ndf,
' degrees of freedom: ', &
4991 '3-sigma limit is',chindl(3,ndf)*real(ndf,mps)
4992 WRITE(1,*) suwt,
' is sum of factors, compared to',nweig, &
4993 ' Downweight fraction:',resing
4999 WRITE(1,*)
' NaNs ', nalc, nrank, nan
5000 WRITE(1,*)
' ---> rejected!'
5002 IF (
mdebug < 0.AND.
nloopn == 1) print *,
' bad local fit-1 ', kfl,jrc,nrc,mside,neq,nalc,nrank,nan,cndl10
5005 IF(nrank /= nalc)
THEN
5009 WRITE(1,*)
' rank deficit', nalc, nrank
5010 WRITE(1,*)
' ---> rejected!'
5012 IF (
mdebug < 0.AND.
nloopn == 1) print *,
' bad local fit-2 ', kfl,jrc,nrc,mside,neq,nalc,nrank,nan,cndl10
5019 WRITE(1,*)
' too large condition(band part) ', nalc, nrank, cndl10
5020 WRITE(1,*)
' ---> rejected!'
5022 IF (
mdebug < 0.AND.
nloopn == 1) print *,
' bad local fit-3 ', kfl,jrc,nrc,mside,neq,nalc,nrank,nan,cndl10
5028 WRITE(1,*)
' Ndf<=0', nalc, nrank, ndf
5029 WRITE(1,*)
' ---> rejected!'
5031 IF (
mdebug < 0.AND.
nloopn == 1) print *,
' bad local fit-4 ', kfl,jrc,nrc,mside,neq,nalc,nrank,nan,cndl10
5035 chndf=real(summ/real(ndf,mpd),mps)
5037 IF(iter == 1.AND.lhist)
CALL hmpent(4,chndf)
5050 chichi=chindl(3,ndf)*real(ndf,mps)
5054 IF(summ >
chhuge*chichi)
THEN
5057 WRITE(1,*)
' ---> rejected!'
5065 IF(summ > chlimt)
THEN
5067 WRITE(1,*)
' ---> rejected!'
5071 CALL addsums(iproc+1, dchi2, ndf, dw1)
5081 CALL addsums(iproc+1, dchi2, ndf, dw1)
5085 WRITE(1,*)
' ---> rejected!'
5098 naccf(kfl)=naccf(kfl)+1
5099 ndff(kfl) =ndff(kfl) +ndf
5100 chi2f(kfl)=chi2f(kfl)+chndf
5113 wght =1.0_mpd/rerr**2
5114 dchi2=wght*rmeas*rmeas
5118 IF(resid > chuber*rerr)
THEN
5119 wght=wght*chuber*rerr/resid
5120 dchi2=2.0*chuber*(resid/rerr-0.5*chuber)
5170 CALL addsums(iproc+1, summ, ndf, dw1)
5226 IF (nfred < 0.OR.nfrei < 0)
THEN
5253 iprcnx=ijprec(i,jnx)
5255 IF (.NOT.((jnx == j+1).AND.(iprc == iprcnx)))
THEN
5269 joffd=joffd+(il*il-il)/2
5281 WRITE(1,*)
'------------------ End of printout for record',nrc
5338 iprcnx=ijprec(i,jnx)
5340 IF (.NOT.((jnx == j+1).AND.(iprc == iprcnx)))
THEN
5355 joffd=joffd+(il*il-il)/2
5391 INTEGER(mpi),
INTENT(IN) :: lun
5393 IF (
nrejec(1)>0)
WRITE(lun,*)
nrejec(1),
' (local solution contains NaNs)'
5394 IF (
nrejec(2)>0)
WRITE(lun,*)
nrejec(2),
' (local matrix with rank deficit)'
5395 IF (
nrejec(3)>0)
WRITE(lun,*)
nrejec(3),
' (local matrix with ill condition)'
5396 IF (
nrejec(4)>0)
WRITE(lun,*)
nrejec(4),
' (local fit with Ndf=0)'
5397 IF (
nrejec(5)>0)
WRITE(lun,*)
nrejec(5),
' (local fit with huge Chi2(Ndf))'
5398 IF (
nrejec(6)>0)
WRITE(lun,*)
nrejec(6),
' (local fit with large Chi2(Ndf))'
5424 INTEGER(mpi) :: icom
5425 INTEGER(mpl) :: icount
5429 INTEGER(mpi) :: imin
5430 INTEGER(mpi) :: iprlim
5431 INTEGER(mpi) :: isub
5432 INTEGER(mpi) :: itgbi
5433 INTEGER(mpi) :: itgbl
5434 INTEGER(mpi) :: ivgbi
5436 INTEGER(mpi) :: label
5444 INTEGER(mpi) :: labele(3)
5445 REAL(mps):: compnt(3)
5450 CALL mvopen(lup,
'millepede.res')
5453 WRITE(*,*)
' Result of fit for global parameters'
5454 WRITE(*,*)
' ==================================='
5459 WRITE(lup,*)
'Parameter ! first 3 elements per line are', &
5460 ' significant (if used as input)'
5478 err=sqrt(abs(real(gmati,mps)))
5479 IF(gmati < 0.0_mpd) err=-err
5482 IF(gmati*diag > 0.0_mpd)
THEN
5483 gcor2=1.0_mpd-1.0_mpd/(gmati*diag)
5484 IF(gcor2 >= 0.0_mpd.AND.gcor2 <= 1.0_mpd) gcor=real(sqrt(gcor2),mps)
5489 IF(lowstat) icount=-(icount+1)
5491 IF(itgbi <= iprlim)
THEN
5505 ELSE IF(itgbi == iprlim+1)
THEN
5506 WRITE(* ,*)
'... (further printout suppressed, but see log file)'
5520 ELSE IF (
igcorr /= 0)
THEN
5538 CALL mvopen(lup,
'millepede.eve')
5548 DO isub=0,min(15,imin-1)
5570 WRITE(lup,103) (labele(ie),compnt(ie),ie=1,iev)
5574 IF(iev /= 0)
WRITE(lup,103) (labele(ie),compnt(ie),ie=1,iev)
5582101
FORMAT(1x,
' label parameter presigma differ', &
5583 ' error'/ 1x,
'-----------',4x,4(
'-------------'))
5584102
FORMAT(i10,2x,4g14.5,f8.3)
5585103
FORMAT(3(i11,f11.7,2x))
5586110
FORMAT(i10,2x,2g14.5,28x,i12)
5587111
FORMAT(i10,2x,3g14.5,14x,i12)
5588112
FORMAT(i10,2x,4g14.5,i12)
5610 INTEGER(mpi) :: icom
5611 INTEGER(mpl) :: icount
5612 INTEGER(mpi) :: ifrst
5613 INTEGER(mpi) :: ilast
5614 INTEGER(mpi) :: inext
5615 INTEGER(mpi) :: itgbi
5616 INTEGER(mpi) :: itgbl
5617 INTEGER(mpi) :: itpgrp
5618 INTEGER(mpi) :: ivgbi
5620 INTEGER(mpi) :: icgrp
5621 INTEGER(mpi) :: ipgrp
5623 INTEGER(mpi) :: jpgrp
5625 INTEGER(mpi) :: label1
5626 INTEGER(mpi) :: label2
5627 INTEGER(mpi) :: ncon
5628 INTEGER(mpi) :: npair
5629 INTEGER(mpi) :: nstep
5632 INTEGER(mpl):: length
5634 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: vecPairedParGroups
5637 SUBROUTINE ggbmap(ipgrp,npair,npgrp)
5639 INTEGER(mpi),
INTENT(IN) :: ipgrp
5640 INTEGER(mpi),
INTENT(OUT) :: npair
5641 INTEGER(mpi),
DIMENSION(:),
INTENT(OUT) :: npgrp
5649 CALL mvopen(lup,
'millepede.res')
5650 WRITE(lup,*)
'*** Results of checking input only, no solution performed ***'
5651 WRITE(lup,*)
'! === global parameters ==='
5652 WRITE(lup,*)
'! fixed-1: by pre-sigma, -2: by entries cut, -3: by iterated entries cut'
5654 WRITE(lup,*)
'! Label Value Pre-sigma SkippedEntries Cons. group Status '
5656 WRITE(lup,*)
'! Label Value Pre-sigma Entries Cons. group Status '
5672 IF (ivgbi <= 0)
THEN
5674 IF (ivgbi == -4)
THEN
5675 WRITE(lup,116) c1,itgbl,par,presig,icount,icgrp
5677 WRITE(lup,110) c1,itgbl,par,presig,icount,icgrp,ivgbi
5681 WRITE(lup,111) c1,itgbl,par,presig,icount,icgrp
5687 WRITE(lup,*)
'!.Appearance statistics '
5688 WRITE(lup,*)
'!. Label First file and record Last file and record #files #paired-par'
5691 IF (itpgrp > 0)
THEN
5699 WRITE(lup,*)
'* === constraint groups ==='
5701 WRITE(lup,*)
'* Group #Cons. Entries First label Last label'
5703 WRITE(lup,*)
'* Group #Cons. Entries First label Last label Paired label range'
5705 CALL mpalloc(vecpairedpargroups,length,
'paired global parameter groups (I)')
5717 IF (
icheck > 1 .AND. label1 > 0)
THEN
5721 vecpairedpargroups(npair+1)=0
5725 jpgrp=vecpairedpargroups(j)
5729 IF (ifrst /= 0.AND.inext /= (ilast+nstep))
THEN
5732 WRITE(lup,114) label1, label2
5737 IF (ifrst == 0) ifrst=inext
5744 IF (jpgrp == vecpairedpargroups(j+1)-1)
THEN
5749 IF (ifrst /= 0)
THEN
5752 WRITE(lup,114) label1, label2
5758 WRITE(lup,*)
'*.Appearance statistics '
5759 WRITE(lup,*)
'*. Group First file and record Last file and record #files'
5769110
FORMAT(
' !',a1,i10,2x,2g14.5,2i12,
' fixed',i2)
5770111
FORMAT(
' !',a1,i10,2x,2g14.5,2i12,
' variable')
5771112
FORMAT(
' !.',i10,6i11)
5772113
FORMAT(
' * ',i6,i8,3i12)
5773114
FORMAT(
' *:',48x,i12,
' ..',i12)
5774115
FORMAT(
' *.',i10,5i11)
5775116
FORMAT(
' !',a1,i10,2x,2g14.5,2i12,
' redundant')
5805 INTEGER(mpi) :: iproc
5815 INTEGER(mpi),
INTENT(IN) :: n
5816 INTEGER(mpl),
INTENT(IN) :: l
5817 REAL(mpd),
INTENT(IN) :: x(n)
5818 INTEGER(mpi),
INTENT(IN) :: is
5819 INTEGER(mpi),
INTENT(IN) :: ie
5820 REAL(mpd),
INTENT(OUT) :: b(n)
5825 INTEGER(mpl) :: indij
5826 INTEGER(mpl) :: indid
5828 INTEGER(mpi) :: ichunk
5864 CALL peend(24,
'Aborted, vector/matrix size mismatch')
5865 stop
'AVPRDS: mismatched vector and matrix'
5885 IF (ia2 <= ib2) b(ia2:ib2)=b(ia2:ib2)+
globalmatd(ia2:ib2)*x(ia2:ib2)
5896 IF (i <= ie.AND.i >= is)
THEN
5904 IF (j <= ie.AND.j >= is)
THEN
5920 IF (ja2 <= jb2)
THEN
5923 b(i)=b(i)+dot_product(
globalmatd(indij+lj+ja2-ja:indij+lj+jb2-ja),x(ja2:jb2))
5927 IF (
mextnd == 0.AND.ia2 <= ib2)
THEN
5930 b(j)=b(j)+dot_product(
globalmatd(indij+lj+jn*(ia2-ia):indij+lj+jn*(ib2-ia):jn),x(ia2:ib2))
5949 IF (i <= ie.AND.i >= is)
THEN
5957 IF (j <= ie.AND.j >= is)
THEN
5973 IF (ja2 <= jb2)
THEN
5976 b(i)=b(i)+dot_product(real(
globalmatf(indij+lj+ja2-ja:indij+lj+jb2-ja),mpd),x(ja2:jb2))
5980 IF (
mextnd == 0.AND.ia2 <= ib2)
THEN
5983 b(j)=b(j)+dot_product(real(
globalmatf(indij+lj+jn*(ia2-ia):indij+lj+jn*(ib2-ia):jn),mpd),x(ia2:ib2))
6017 INTEGER(mpi) :: iproc
6025 INTEGER(mpi),
INTENT(IN) :: n
6026 INTEGER(mpl),
INTENT(IN) :: l
6027 REAL(mpd),
INTENT(IN) :: x(n)
6028 REAL(mpd),
INTENT(OUT) :: b(n)
6033 INTEGER(mpl) :: indij
6034 INTEGER(mpl) :: indid
6036 INTEGER(mpi) :: ichunk
6065 CALL peend(24,
'Aborted, vector/matrix size mismatch')
6066 stop
'AVPRD0: mismatched vector and matrix'
6112 b(i)=b(i)+dot_product(
globalmatd(indij+lj:indij+lj+jn-1),x(ja:jb))
6118 b(j)=b(j)+dot_product(
globalmatd(indij+lj:indij+jn*in:jn),x(ia:ib))
6154 b(i)=b(i)+dot_product(real(
globalmatf(indij+lj:indij+lj+jn-1),mpd),x(ja:jb))
6160 b(j)=b(j)+dot_product(real(
globalmatf(indij+lj:indij+jn*in:jn),mpd),x(ia:ib))
6184 INTEGER(mpi) :: ispc
6195 INTEGER(mpl) :: indid
6196 INTEGER(mpl),
DIMENSION(12) :: icount
6203 icount(4)=huge(icount(4))
6204 icount(7)=huge(icount(7))
6205 icount(10)=huge(icount(10))
6215 ir=ipg+(ispc-1)*(
napgrp+1)
6222 icount(1)=icount(1)+in
6223 icount(2)=icount(2)+in*ku
6229 icount(3)=icount(3)+1
6230 icount(4)=min(icount(4),jn)
6231 icount(5)=icount(5)+jn
6232 icount(6)=max(icount(6),jn)
6233 icount(7)=min(icount(7),in)
6234 icount(8)=icount(8)+in
6235 icount(9)=max(icount(9),in)
6236 icount(10)=min(icount(10),in*jn)
6237 icount(11)=icount(11)+in*jn
6238 icount(12)=max(icount(12),in*jn)
6244 WRITE(*,*)
"analysis of sparsity structure"
6245 IF (icount(1) > 0)
THEN
6246 WRITE(*,101)
"rows without compression/blocks ", icount(1)
6247 WRITE(*,101)
" contained elements ", icount(2)
6249 WRITE(*,101)
"number of block matrices ", icount(3)
6250 avg=real(icount(5),mps)/real(icount(3),mps)
6251 WRITE(*,101)
"number of columns (min,mean,max) ", icount(4), avg, icount(6)
6252 avg=real(icount(8),mps)/real(icount(3),mps)
6253 WRITE(*,101)
"number of rows (min,mean,max) ", icount(7), avg, icount(9)
6254 avg=real(icount(11),mps)/real(icount(3),mps)
6255 WRITE(*,101)
"number of elements (min,mean,max) ", icount(10), avg, icount(12)
6256101
FORMAT(2x,a34,i10,f10.3,i10)
6275 INTEGER(mpi),
INTENT(IN) :: n
6276 REAL(mpd),
INTENT(IN) :: x(n)
6277 REAL(mpd),
INTENT(OUT) :: b(n)
6282 CALL peend(24,
'Aborted, vector/matrix size mismatch')
6283 stop
'AVPROD: mismatched vector and matrix'
6314 INTEGER(mpi) :: ispc
6315 INTEGER(mpi) :: item1
6316 INTEGER(mpi) :: item2
6317 INTEGER(mpi) :: itemc
6318 INTEGER(mpi) :: jtem
6319 INTEGER(mpi) :: jtemn
6322 INTEGER(mpi),
INTENT(IN) :: itema
6323 INTEGER(mpi),
INTENT(IN) :: itemb
6324 INTEGER(mpl),
INTENT(OUT) :: ij
6325 INTEGER(mpi),
INTENT(OUT) :: lr
6326 INTEGER(mpi),
INTENT(OUT) :: iprc
6337 item1=max(itema,itemb)
6338 item2=min(itema,itemb)
6339 IF(item2 <= 0.OR.item1 >
napgrp)
RETURN
6342 outer:
DO ispc=1,
nspc
6353 IF(ku < kl) cycle outer
6358 IF(item2 >= jtem.AND.item2 < jtemn)
THEN
6364 IF(item2 < jtem)
THEN
6366 ELSE IF(item2 >= jtemn)
THEN
6381 IF(ku < kl) cycle outer
6385 IF(itemc == jtem)
EXIT
6386 IF(itemc < jtem)
THEN
6388 ELSE IF(itemc > jtem)
THEN
6416 INTEGER(mpi),
INTENT(IN) :: itema
6417 INTEGER(mpi),
INTENT(IN) :: itemb
6442 INTEGER(mpi) :: item1
6443 INTEGER(mpi) :: item2
6444 INTEGER(mpi) :: ipg1
6445 INTEGER(mpi) :: ipg2
6447 INTEGER(mpi) :: iprc
6449 INTEGER(mpi),
INTENT(IN) :: itema
6450 INTEGER(mpi),
INTENT(IN) :: itemb
6452 INTEGER(mpl) ::
ijadd
6456 item1=max(itema,itemb)
6457 item2=min(itema,itemb)
6459 IF(item2 <= 0.OR.item1 >
nagb)
RETURN
6460 IF(item1 == item2)
THEN
6469 CALL ijpgrp(ipg1,ipg2,ij,lr,iprc)
6491 INTEGER(mpi) :: item1
6492 INTEGER(mpi) :: item2
6493 INTEGER(mpi) :: jtem
6495 INTEGER(mpi),
INTENT(IN) :: itema
6496 INTEGER(mpi),
INTENT(IN) :: itemb
6505 item1=max(itema,itemb)
6506 item2=min(itema,itemb)
6508 IF(item2 <= 0.OR.item1 >
nagb)
RETURN
6516 print *,
' IJCSR3 empty list ', item1, item2, ks, ke
6517 CALL peend(23,
'Aborted, bad matrix index')
6518 stop
'ijcsr3: empty list'
6523 IF(item1 == jtem)
EXIT
6524 IF(item1 < jtem)
THEN
6531 print *,
' IJCSR3 not found ', item1, item2, ks, ke
6532 CALL peend(23,
'Aborted, bad matrix index')
6533 stop
'ijcsr3: not found'
6549 INTEGER(mpi) :: item1
6550 INTEGER(mpi) :: item2
6554 INTEGER(mpl) ::
ijadd
6557 INTEGER(mpi),
INTENT(IN) :: itema
6558 INTEGER(mpi),
INTENT(IN) :: itemb
6563 item1=max(itema,itemb)
6564 item2=min(itema,itemb)
6565 IF(item2 <= 0.OR.item1 >
nagb)
RETURN
6574 ij=
ijadd(item1,item2)
6577 ELSE IF (ij < 0)
THEN
6609 INTEGER(mpi) :: ichunk
6613 INTEGER(mpi) :: ispc
6621 INTEGER(mpl) :: ijadd
6643 ir=ipg+(ispc-1)*(
napgrp+1)
6693 REAL(mps),
INTENT(IN) :: deltat
6694 INTEGER(mpi),
INTENT(OUT) :: minut
6695 INTEGER(mpi),
INTENT(OUT):: nhour
6696 REAL(mps),
INTENT(OUT):: secnd
6697 INTEGER(mpi) :: nsecd
6700 nsecd=nint(deltat,
mpi)
6702 minut=nsecd/60-60*nhour
6703 secnd=deltat-60*(minut+60*nhour)
6739 INTEGER(mpi),
INTENT(IN) :: item
6743 INTEGER(mpl) :: length
6744 INTEGER(mpl),
PARAMETER :: four = 4
6748 IF(item <= 0)
RETURN
6769 IF(j == 0)
EXIT inner
6791 IF (
lvllog > 1)
WRITE(
lunlog,*)
'INONE: array increased to', &
6812 INTEGER(mpi) :: iprime
6813 INTEGER(mpi) :: nused
6814 LOGICAL :: finalUpdate
6815 INTEGER(mpl) :: oldLength
6816 INTEGER(mpl) :: newLength
6817 INTEGER(mpl),
PARAMETER :: four = 4
6818 INTEGER(mpi),
DIMENSION(:,:),
ALLOCATABLE :: tempArr
6819 INTEGER(mpl),
DIMENSION(:),
ALLOCATABLE :: tempVec
6823 IF(finalupdate)
THEN
6832 CALL mpalloc(temparr,four,oldlength,
'INONE: temp array')
6834 CALL mpalloc(tempvec,oldlength,
'INONE: temp vector')
6858 IF(j == 0)
EXIT inner
6859 IF(j == i) cycle outer
6863 IF(.NOT.finalupdate)
RETURN
6868 WRITE(
lunlog,*)
'INONE: array reduced to',newlength,
' words'
6892 IF(j == 0)
EXIT inner
6893 IF(j == i) cycle outer
6910 INTEGER(mpi),
INTENT(IN) :: n
6911 INTEGER(mpi) :: nprime
6912 INTEGER(mpi) :: nsqrt
6917 IF(mod(nprime,2) == 0) nprime=nprime+1
6920 nsqrt=int(sqrt(real(nprime,
mps)),
mpi)
6922 IF(i*(nprime/i) == nprime) cycle outer
6944 INTEGER(mpi) :: idum
6946 INTEGER(mpi) :: indab
6947 INTEGER(mpi) :: itgbi
6948 INTEGER(mpi) :: itgbl
6949 INTEGER(mpi) :: ivgbi
6951 INTEGER(mpi) :: jgrp
6952 INTEGER(mpi) :: lgrp
6954 INTEGER(mpi) :: nc31
6956 INTEGER(mpi) :: nwrd
6957 INTEGER(mpi) :: inone
6962 INTEGER(mpl) :: length
6963 INTEGER(mpl) :: rows
6967 WRITE(
lunlog,*)
'LOOP1: starting'
6990 WRITE(
lunlog,*)
'LOOP1: reading data files'
7004 WRITE(*,*)
'Read all binary data files:'
7006 CALL hmpldf(1,
'Number of words/record in binary file')
7007 CALL hmpdef(8,0.0,60.0,
'not_stored data per record')
7037 WRITE(
lunlog,*)
'LOOP1: reading data files again'
7049 CALL peend(21,
'Aborted, no labels/parameters defined')
7050 stop
'LOOP1: no labels/parameters defined'
7055 ' is total number NTGB of labels/parameters'
7057 CALL hmpldf(2,
'Number of entries per label')
7101 WRITE(
lunlog,*)
'LOOP1:',
npresg,
' is number of pre-sigmas'
7102 WRITE(*,*)
'LOOP1:',
npresg,
' is number of pre-sigmas'
7103 IF(
npresg == 0)
WRITE(*,*)
'Warning: no pre-sigmas defined'
7125 WRITE(
lunlog,*)
'LOOP1:',
nvgb,
' is number NVGB of variable parameters'
7130 WRITE(
lunlog,*)
'LOOP1: counting records, NO iteration of entries cut !'
7136 CALL hmpdef(15,0.0,120.0,
'Number of parameters per group')
7171 ' is total number NTPGRP of label/parameter groups'
7187 WRITE(*,112)
' Default pre-sigma =',
regpre, &
7188 ' (if no individual pre-sigma defined)'
7189 WRITE(*,*)
'Pre-sigma factor is',
regula
7192 WRITE(*,*)
'No regularization will be done'
7194 WRITE(*,*)
'Regularization will be done, using factor',
regula
7198 CALL peend(22,
'Aborted, no variable global parameters')
7199 stop
'... no variable global parameters'
7206 IF(presg > 0.0)
THEN
7208 ELSE IF(presg == 0.0.AND.
regpre > 0.0)
THEN
7209 prewt=1.0/real(
regpre**2,mpd)
7225 WRITE(*,101)
'NTGB',
ntgb,
'total number of parameters'
7226 WRITE(*,101)
'NVGB',
nvgb,
'number of variable parameters'
7230 WRITE(*,101)
'Too many variable parameters for diagonalization, fallback is inversion'
7238 WRITE(*,101)
' NREC',
nrec,
'number of records'
7239 IF (
nrecd > 0)
WRITE(*,101)
' NRECD',
nrec,
'number of records containing doubles'
7240 WRITE(*,101)
' NEQN',
neqn,
'number of equations (measurements)'
7241 WRITE(*,101)
' NEGB',
negb,
'number of equations with global parameters'
7242 WRITE(*,101)
' NDGB',
ndgb,
'number of global derivatives'
7244 WRITE(*,101)
' NZGB',
nzgb,
'number of zero global der. (ignored in entry counts)'
7247 WRITE(*,101)
'MREQENF',
mreqenf,
'required number of entries (eqns in binary files)'
7249 WRITE(*,101)
'MREQENF',
mreqenf,
'required number of entries (recs in binary files)'
7252 WRITE(*,101)
'ITEREN',
iteren,
'iterate cut for parameters with less entries'
7253 WRITE(*,101)
'MREQENA',
mreqena,
'required number of entries (from accepted fits)'
7254 IF (
mreqpe > 1)
WRITE(*,101) &
7255 'MREQPE',
mreqpe,
'required number of pair entries'
7256 IF (
msngpe >= 1)
WRITE(*,101) &
7257 'MSNGPE',
msngpe,
'max pair entries single prec. storage'
7258 WRITE(*,101)
'NTGB',
ntgb,
'total number of parameters'
7259 WRITE(*,101)
'NVGB',
nvgb,
'number of variable parameters'
7262 WRITE(*,*)
'Global parameter labels:'
7269 mqi=((mqi-20)/20)*20+1
7277 WRITE(8,101)
' NREC',
nrec,
'number of records'
7278 IF (
nrecd > 0)
WRITE(8,101)
' NRECD',
nrec,
'number of records containing doubles'
7279 WRITE(8,101)
' NEQN',
neqn,
'number of equations (measurements)'
7280 WRITE(8,101)
' NEGB',
negb,
'number of equations with global parameters'
7281 WRITE(8,101)
' NDGB',
ndgb,
'number of global derivatives'
7283 WRITE(8,101)
'MREQENF',
mreqenf,
'required number of entries (eqns in binary files)'
7285 WRITE(8,101)
'MREQENF',
mreqenf,
'required number of entries (recs in binary files)'
7288 WRITE(8,101)
'ITEREN',
iteren,
'iterate cut for parameters with less entries'
7289 WRITE(8,101)
'MREQENA',
mreqena,
'required number of entries (from accepted fits)'
7291 WRITE(
lunlog,*)
'LOOP1: ending'
7295101
FORMAT(1x,a8,
' =',i14,
' = ',a)
7311 INTEGER(mpi) :: ibuf
7313 INTEGER(mpi) :: indab
7319 INTEGER(mpi) :: nc31
7321 INTEGER(mpi) :: nlow
7323 INTEGER(mpi) :: nwrd
7325 INTEGER(mpl) :: length
7326 INTEGER(mpl),
DIMENSION(:),
ALLOCATABLE :: newCounter
7331 WRITE(
lunlog,*)
'LOOP1: iterating'
7333 WRITE(*,*)
'LOOP1: iterating'
7336 CALL mpalloc(newcounter,length,
'new entries counter')
7360 CALL isjajb(nst,ist,ja,jb,jsp)
7361 IF(ja == 0.AND.jb == 0)
EXIT
7367 IF(ij == -2) nlow=nlow+1
7372 newcounter(ij)=newcounter(ij)+1
7401 WRITE(
lunlog,*)
'LOOP1:',
nvgb,
' is number NVGB of variable parameters'
7431 INTEGER(mpi) :: ibuf
7432 INTEGER(mpi) :: icblst
7433 INTEGER(mpi) :: icboff
7434 INTEGER(mpi) :: icgb
7435 INTEGER(mpi) :: icgrp
7436 INTEGER(mpi) :: icount
7437 INTEGER(mpi) :: iext
7438 INTEGER(mpi) :: ihis
7442 INTEGER(mpi) :: ioff
7443 INTEGER(mpi) :: ipoff
7444 INTEGER(mpi) :: iproc
7445 INTEGER(mpi) :: irecmm
7447 INTEGER(mpi) :: itgbi
7448 INTEGER(mpi) :: itgbij
7449 INTEGER(mpi) :: itgbik
7450 INTEGER(mpi) :: ivgbij
7451 INTEGER(mpi) :: ivgbik
7452 INTEGER(mpi) :: ivpgrp
7456 INTEGER(mpi) :: jcgrp
7457 INTEGER(mpi) :: jext
7458 INTEGER(mpi) :: jcgb
7459 INTEGER(mpi) :: jrec
7461 INTEGER(mpi) :: joff
7463 INTEGER(mpi) :: kcgrp
7464 INTEGER(mpi) :: kfile
7466 INTEGER(mpi) :: label
7467 INTEGER(mpi) :: labelf
7468 INTEGER(mpi) :: labell
7469 INTEGER(mpi) :: lvpgrp
7472 INTEGER(mpi) :: maeqnf
7473 INTEGER(mpi) :: nall
7474 INTEGER(mpi) :: naeqna
7475 INTEGER(mpi) :: naeqnf
7476 INTEGER(mpi) :: naeqng
7477 INTEGER(mpi) :: npdblk
7478 INTEGER(mpi) :: nc31
7479 INTEGER(mpi) :: ncachd
7480 INTEGER(mpi) :: ncachi
7481 INTEGER(mpi) :: ncachr
7482 INTEGER(mpi) :: ncon
7485 INTEGER(mpi) :: ndfmax
7486 INTEGER(mpi) :: nfixed
7487 INTEGER(mpi) :: nggd
7488 INTEGER(mpi) :: nggi
7489 INTEGER(mpi) :: nmatmo
7490 INTEGER(mpi) :: noff
7491 INTEGER(mpi) :: npair
7492 INTEGER(mpi) :: npar
7493 INTEGER(mpi) :: nparmx
7495 INTEGER(mpi) :: nrece
7496 INTEGER(mpi) :: nrecf
7497 INTEGER(mpi) :: nrecmm
7499 INTEGER(mpi) :: nwrd
7500 INTEGER(mpi) :: inone
7509 INTEGER(mpl):: nblock
7510 INTEGER(mpl):: nbwrds
7511 INTEGER(mpl):: noff8
7512 INTEGER(mpl):: ndimbi
7513 INTEGER(mpl):: ndimsa(4)
7515 INTEGER(mpl):: nnzero
7516 INTEGER(mpl):: matsiz(2)
7517 INTEGER(mpl):: matwords
7518 INTEGER(mpl):: mbwrds
7519 INTEGER(mpl):: length
7522 INTEGER(mpl),
PARAMETER :: two=2
7523 INTEGER(mpi) :: maxGlobalPar = 0
7524 INTEGER(mpi) :: maxLocalPar = 0
7525 INTEGER(mpi) :: maxEquations = 0
7527 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: vecConsGroupList
7528 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: vecConsGroupIndex
7529 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: vecPairedParGroups
7530 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: vecBlockCounts
7533 SUBROUTINE ndbits(npgrp,ndims,nsparr,ihst)
7535 INTEGER(mpi),
DIMENSION(:),
INTENT(IN) :: npgrp
7536 INTEGER(mpl),
DIMENSION(4),
INTENT(OUT) :: ndims
7537 INTEGER(mpl),
DIMENSION(:,:),
INTENT(OUT) :: nsparr
7538 INTEGER(mpi),
INTENT(IN) :: ihst
7540 SUBROUTINE ckbits(npgrp,ndims)
7542 INTEGER(mpi),
DIMENSION(:),
INTENT(IN) :: npgrp
7543 INTEGER(mpl),
DIMENSION(4),
INTENT(OUT) :: ndims
7545 SUBROUTINE spbits(npgrp,nsparr,nsparc)
7547 INTEGER(mpi),
DIMENSION(:),
INTENT(IN) :: npgrp
7548 INTEGER(mpl),
DIMENSION(:,:),
INTENT(IN) :: nsparr
7549 INTEGER(mpi),
DIMENSION(:),
INTENT(OUT) :: nsparc
7551 SUBROUTINE gpbmap(ngroup,npgrp,npair)
7553 INTEGER(mpi),
INTENT(IN) :: ngroup
7554 INTEGER(mpi),
DIMENSION(:,:),
INTENT(IN) :: npgrp
7555 INTEGER(mpi),
DIMENSION(:),
INTENT(OUT) :: npair
7557 SUBROUTINE ggbmap(ipgrp,npair,npgrp)
7559 INTEGER(mpi),
INTENT(IN) :: ipgrp
7560 INTEGER(mpi),
INTENT(OUT) :: npair
7561 INTEGER(mpi),
DIMENSION(:),
INTENT(OUT) :: npgrp
7563 SUBROUTINE pbsbits(npgrp,ibsize,nnzero,nblock,nbkrow)
7565 INTEGER(mpi),
DIMENSION(:),
INTENT(IN) :: npgrp
7566 INTEGER(mpi),
INTENT(IN) :: ibsize
7567 INTEGER(mpl),
INTENT(OUT) :: nnzero
7568 INTEGER(mpl),
INTENT(OUT) :: nblock
7569 INTEGER(mpi),
DIMENSION(:),
INTENT(OUT) :: nbkrow
7571 SUBROUTINE pblbits(npgrp,ibsize,nsparr,nsparc)
7573 INTEGER(mpi),
DIMENSION(:),
INTENT(IN) :: npgrp
7574 INTEGER(mpi),
INTENT(IN) :: ibsize
7575 INTEGER(mpl),
DIMENSION(:),
INTENT(IN) :: nsparr
7576 INTEGER(mpl),
DIMENSION(:),
INTENT(OUT) :: nsparc
7578 SUBROUTINE prbits(npgrp,nsparr)
7580 INTEGER(mpi),
DIMENSION(:),
INTENT(IN) :: npgrp
7581 INTEGER(mpl),
DIMENSION(:),
INTENT(OUT) :: nsparr
7583 SUBROUTINE pcbits(npgrp,nsparr,nsparc)
7585 INTEGER(mpi),
DIMENSION(:),
INTENT(IN) :: npgrp
7586 INTEGER(mpl),
DIMENSION(:),
INTENT(IN) :: nsparr
7587 INTEGER(mpl),
DIMENSION(:),
INTENT(OUT) :: nsparc
7597 WRITE(
lunlog,*)
'LOOP2: starting'
7620 WRITE(
lunlog,*)
' Elimination for constraints enforced for solution by decomposition!'
7625 WRITE(
lunlog,*)
' Lagrange multipliers enforced for solution by sparsePARDISO!'
7630 WRITE(
lunlog,*)
' Elimination for constraints with mpqldec enforced (LAPACK only for unpacked storage)!'
7647 noff8=int(
nagb,mpl)*int(
nagb-1,mpl)/2
7683 print *,
' Variable parameter groups ',
nvpgrp
7731 CALL mpalloc(vecconsgrouplist,length,
'constraint group list')
7732 CALL mpalloc(vecconsgroupindex,length,
'constraint group index')
7742 WRITE(
lunlog,*)
'LOOP2: start event reading'
7782 WRITE(*,*)
'Record number ',
nrec,
' from file ',kfile
7783 IF (wgh /= 1.0)
WRITE(*,*)
' weight ',wrec
7787 CALL isjajb(nst,ist,ja,jb,jsp)
7791 IF(nda ==
mdebg2+1)
WRITE(*,*)
'... and more data'
7796 WRITE(*,*)
'Local derivatives:'
7798107
FORMAT(6(i3,g12.4))
7800 WRITE(*,*)
'Global derivatives:'
7803108
FORMAT(3i11,g12.4)
7806 WRITE(*,*)
'total_par_label __label__ var_par_index derivative'
7821 CALL isjajb(nst,ist,ja,jb,jsp)
7822 IF(ja == 0.AND.jb == 0)
EXIT
7862 IF (kcgrp == jcgrp) cycle
7873 IF (vecconsgroupindex(k) == 0)
THEN
7876 vecconsgrouplist(icgrp)=k
7891 IF (vecconsgroupindex(k) < icount)
THEN
7893 vecconsgroupindex(k)=icount
7911 IF (nfixed > 0) naeqnf=naeqnf+1
7914 IF(ja /= 0.AND.jb /= 0)
THEN
7923 IF (naeqnf > maeqnf) nrecf=nrecf+1
7927 maxglobalpar=max(
nagbn,maxglobalpar)
7928 maxlocalpar=max(
nalcn,maxlocalpar)
7929 maxequations=max(
naeqn,maxequations)
7932 dstat(1)=dstat(1)+real((nwrd+2)*2,mpd)
7933 dstat(2)=dstat(2)+real(
nagbn+2,mpd)
7938 vecconsgroupindex(vecconsgrouplist(k))=0
7943 IF (
nagbn == 0)
THEN
7962 IF (ivpgrp /= lvpgrp)
THEN
8017 (irecmm >= nrecmm.OR.irecmm ==
mxrec))
THEN
8018 IF (nmatmo == 0)
THEN
8020 WRITE(*,*)
'Monitoring of sparse matrix construction'
8021 WRITE(*,*)
' records ........ off-diagonal elements ', &
8022 '....... compression memory'
8023 WRITE(*,*)
' non-zero used(double) used', &
8028 gbc=1.0e-9*real((mpi*ndimsa(2)+mpd*ndimsa(3)+mps*ndimsa(4))/mpi*(bit_size(1_mpi)/8),mps)
8029 gbu=1.0e-9*real(((mpi+mpd)*(ndimsa(3)+ndimsa(4)))/mpi*(bit_size(1_mpi)/8),mps)
8031 WRITE(*,1177) irecmm,ndimsa(1),ndimsa(3),ndimsa(4),cpr,gbc
80321177
FORMAT(i9,3i13,f10.2,f11.6)
8033 DO WHILE(irecmm >= nrecmm)
8052 WRITE(
lunlog,*)
'LOOP2: event reading ended - end of data'
8054 dstat(k)=dstat(k)/real(
nrec,mpd)
8068 CALL mpalloc(vecpairedpargroups,length,
'paired global parameter groups (I)')
8070 print *,
' Total parameter groups pairs',
ntpgrp
8073 CALL ggbmap(i,npair,vecpairedpargroups)
8141 IF(ivgbij > 0.AND.ivgbik > 0)
THEN
8143 IF (
mprint > 1)
WRITE(*,*)
'add index pair ',ivgbij,ivgbik
8220 nparmx=max(nparmx,int(rows,mpi))
8242 WRITE(*,*)
' Parameter block', i, ib-ia, jb-ja, labelf, labell
8246 WRITE(
lunlog,*)
'Detected',
npblck,
'(disjoint) parameter blocks, max size ', nparmx
8248 WRITE(*,*)
'Detected',
npblck,
'(disjoint) parameter blocks, max size ', nparmx
8250 WRITE(*,*)
'Using block diagonal storage mode'
8265 IF (
nagb >= 65536)
THEN
8266 noff=int(noff8/1000,mpi)
8276 CALL hmpdef(ihis,0.0,real(
mhispe,mps),
'NDBITS: #off-diagonal elements')
8281 ndgn=ndimsa(3)+ndimsa(4)
8282 matwords=ndimsa(2)+length*4
8297 length=int(npdblk,mpl)
8298 CALL mpalloc(vecblockcounts,length,
'sparse matrix row offsets (CSR3)')
8300 nbwrds=2*int(nblock,mpl)*int(
ipdbsz(i)*
ipdbsz(i)+1,mpl)
8301 IF ((i == 1).OR.(nbwrds < mbwrds))
THEN
8325 IF (
fcache(2) == 0.0)
THEN
8327 fcache(2)=real(dstat(2),mps)
8328 fcache(3)=real(dstat(3),mps)
8349 ncachd=
ncache-ncachr-ncachi
8380 WRITE(lu,101)
'NTGB',
ntgb,
'total number of parameters'
8381 WRITE(lu,102)
'(all parameters, appearing in binary files)'
8382 WRITE(lu,101)
'NVGB',
nvgb,
'number of variable parameters'
8383 WRITE(lu,102)
'(appearing in fit matrix/vectors)'
8384 WRITE(lu,101)
'NAGB',
nagb,
'number of all parameters'
8385 WRITE(lu,102)
'(including Lagrange multiplier or reduced)'
8386 WRITE(lu,101)
'NTPGRP',
ntpgrp,
'total number of parameter groups'
8387 WRITE(lu,101)
'NVPGRP',
nvpgrp,
'number of variable parameter groups'
8388 WRITE(lu,101)
'NFGB',
nfgb,
'number of fit parameters'
8390 WRITE(lu,101)
'MBANDW',
mbandw,
'band width of preconditioner matrix'
8391 WRITE(lu,102)
'(if <0, no preconditioner matrix)'
8393 IF (
nagb >= 65536)
THEN
8394 WRITE(lu,101)
'NOFF/K',noff,
'max number of off-diagonal elements'
8396 WRITE(lu,101)
'NOFF',noff,
'max number of off-diagonal elements'
8399 IF (
nagb >= 65536)
THEN
8400 WRITE(lu,101)
'NDGN/K',ndgn/1000,
'actual number of off-diagonal elements'
8402 WRITE(lu,101)
'NDGN',ndgn,
'actual number of off-diagonal elements'
8405 WRITE(lu,101)
'NCGB',
ncgb,
'number of constraints'
8406 WRITE(lu,101)
'NAGBN',
nagbn,
'max number of global parameters in an event'
8407 WRITE(lu,101)
'NALCN',
nalcn,
'max number of local parameters in an event'
8408 WRITE(lu,101)
'NAEQN',
naeqn,
'max number of equations in an event'
8410 WRITE(lu,101)
'NAEQNA',naeqna,
'number of equations'
8411 WRITE(lu,101)
'NAEQNG',naeqng, &
8412 'number of equations with global parameters'
8413 WRITE(lu,101)
'NAEQNF',naeqnf, &
8414 'number of equations with fixed global parameters'
8415 WRITE(lu,101)
'NRECF',nrecf, &
8416 'number of records with fixed global parameters'
8419 WRITE(lu,101)
'NRECE',nrece, &
8420 'number of records without variable parameters'
8423 WRITE(lu,101)
'NCACHE',
ncache,
'number of words for caching'
8424 WRITE(lu,111) (
fcache(k)*100.0,k=1,3)
8425111
FORMAT(22x,
'cache splitting ',3(f6.1,
' %'))
8430 WRITE(lu,*)
'Solution method and matrix-storage mode:'
8432 WRITE(lu,*)
' METSOL = 1: matrix inversion'
8433 ELSE IF(
metsol == 2)
THEN
8434 WRITE(lu,*)
' METSOL = 2: diagonalization'
8435 ELSE IF(
metsol == 3)
THEN
8436 WRITE(lu,*)
' METSOL = 3: decomposition'
8437 ELSE IF(
metsol == 4)
THEN
8438 WRITE(lu,*)
' METSOL = 4: MINRES (rtol',
mrestl,
')'
8439 ELSE IF(
metsol == 5)
THEN
8440 WRITE(lu,*)
' METSOL = 5: MINRES-QLP (rtol',
mrestl,
')'
8441 ELSE IF(
metsol == 6)
THEN
8442 WRITE(lu,*)
' METSOL = 6: GMRES'
8444 ELSE IF(
metsol == 7)
THEN
8445 WRITE(lu,*)
' METSOL = 7: LAPACK factorization'
8446 ELSE IF(
metsol == 8)
THEN
8447 WRITE(lu,*)
' METSOL = 8: LAPACK factorization'
8449 ELSE IF(
metsol == 9)
THEN
8450 WRITE(lu,*)
' METSOL = 9: Intel oneMKL PARDISO'
8454 WRITE(lu,*)
' with',
mitera,
' iterations'
8456 WRITE(lu,*)
' MATSTO = 0: unpacked symmetric matrix, ',
'n*n elements'
8457 ELSE IF(
matsto == 1)
THEN
8458 WRITE(lu,*)
' MATSTO = 1: full symmetric matrix, ',
'(n*n+n)/2 elements'
8459 ELSE IF(
matsto == 2)
THEN
8460 WRITE(lu,*)
' MATSTO = 2: sparse matrix (custom)'
8461 ELSE IF(
matsto == 3)
THEN
8463 WRITE(lu,*)
' MATSTO = 3: sparse matrix (upper triangle, CSR3)'
8465 WRITE(lu,*)
' MATSTO = 3: sparse matrix (upper triangle, BSR3)'
8466 WRITE(lu,*)
' block size',
matbsz
8470 WRITE(lu,*)
' block diagonal with',
npblck,
' blocks'
8472 IF(
mextnd>0)
WRITE(lu,*)
' with extended storage'
8473 IF(
dflim /= 0.0)
THEN
8474 WRITE(lu,103)
'Convergence assumed, if expected dF <',
dflim
8479 WRITE(lu,*)
'Constraints handled by elimination with LAPACK'
8481 WRITE(lu,*)
'Constraints handled by elimination'
8484 WRITE(lu,*)
'Constraints handled by Lagrange multipliers'
8491 CALL peend(28,
'Aborted, no local parameters')
8492 stop
'LOOP2: stopping due to missing local parameters'
8513105
FORMAT(
' Constants C1, C2 for Wolfe conditions:',g12.4,
', ',g12.4)
8521 matwords=(length+
nagb+1)*2
8528 ELSE IF (
matsto == 2)
THEN
8529 matsiz(1)=ndimsa(3)+
nagb
8544 IF (icblst > icboff)
THEN
8550 nall = npar;
IF (
icelim <= 0) nall=npar+ncon
8554 matsiz(1)=matsiz(1)+k
8556 matsiz(1)=matsiz(1)+nall
8561 matwords=matwords+matsiz(1)*2+matsiz(2)
8567 IF (matwords < 250000)
THEN
8568 WRITE(*,*)
'Size of global matrix: < 1 MB'
8570 WRITE(*,*)
'Size of global matrix:',int(real(matwords,mps)*4.0e-6,mpi),
' MB'
8576 WRITE(
lunlog,*)
' Cut values of Chi^2/Ndf and Chi2,'
8577 WRITE(
lunlog,*)
' corresponding to 2 and 3 standard deviations'
8578 WRITE(
lunlog,*)
' Ndf Chi^2/Ndf(2) Chi^2(2) ', &
8579 ' Chi^2/Ndf(3) Chi^2(3)'
8582 IF(ndf >
naeqn)
EXIT
8585 ELSE IF(ndf < 20)
THEN
8587 ELSE IF(ndf < 100)
THEN
8589 ELSE IF(ndf < 200)
THEN
8596 WRITE(
lunlog,106) ndf,chin2,chin2*real(ndf,mps),chin3, chin3*real(ndf,mps)
8599 WRITE(
lunlog,*)
'LOOP2: ending'
8603 IF (
ncgbe /= 0)
THEN
8606 WRITE(*,199)
'WarningWarningWarningWarningWarningWarningWarningWarningWar'
8607 WRITE(*,199)
'arningWarningWarningWarningWarningWarningWarningWarningWarn'
8608 WRITE(*,199)
'rningWarningWarningWarningWarningWarningWarningWarningWarni'
8609 WRITE(*,199)
'ningWarningWarningWarningWarningWarningWarningWarningWarnin'
8610 WRITE(*,199)
'ingWarningWarningWarningWarningWarningWarningWarningWarning'
8611 WRITE(*,199)
'ngWarningWarningWarningWarningWarningWarningWarningWarningW'
8612 WRITE(*,199)
'gWarningWarningWarningWarningWarningWarningWarningWarningWa'
8614 WRITE(*,*)
' Number of empty constraints =',abs(
ncgbe),
', should be 0'
8615 WRITE(*,*)
' => please check constraint definition, mille data'
8617 WRITE(*,199)
'WarningWarningWarningWarningWarningWarningWarningWarningWar'
8618 WRITE(*,199)
'arningWarningWarningWarningWarningWarningWarningWarningWarn'
8619 WRITE(*,199)
'rningWarningWarningWarningWarningWarningWarningWarningWarni'
8620 WRITE(*,199)
'ningWarningWarningWarningWarningWarningWarningWarningWarnin'
8621 WRITE(*,199)
'ingWarningWarningWarningWarningWarningWarningWarningWarning'
8622 WRITE(*,199)
'ngWarningWarningWarningWarningWarningWarningWarningWarningW'
8623 WRITE(*,199)
'gWarningWarningWarningWarningWarningWarningWarningWarningWa'
8628101
FORMAT(1x,a8,
' =',i14,
' = ',a)
8630103
FORMAT(1x,a,g12.4)
8631106
FORMAT(i6,2(3x,f9.3,f12.1,3x))
8646 INTEGER(mpi) :: imed
8649 INTEGER(mpi) :: nent
8650 INTEGER(mpi),
DIMENSION(measBins) :: isuml
8651 INTEGER(mpi),
DIMENSION(measBins) :: isums
8655 INTEGER(mpl) :: ioff
8658 DATA lfirst /.true./
8671 WRITE(
lunmon,
'(A)')
'*** Normalized residuals grouped by first global label (per local fit cycle) ***'
8673 WRITE(
lunmon,
'(A)')
'*** Pulls grouped by first global label (per local fit cycle) ***'
8675 WRITE(
lunmon,
'(A)')
'! LFC Label Entries Median RMS(MAD) <error>'
8680#ifdef SCOREP_USER_ENABLE
8681 scorep_user_region_by_name_begin(
"UR_monres", scorep_user_region_type_common)
8697 IF (2*isuml(j) > nent)
EXIT
8701 IF (isuml(j) > isuml(j-1)) amed=amed+real(nent-2*isuml(j-1),mps)/real(2*isuml(j)-2*isuml(j-1),mps)
8714 isums(j)=isums(j)+isums(j-1)
8718 IF (2*isums(j) > nent)
EXIT
8721 IF (isums(j) > isums(j-1)) amad=amad+real(nent-2*isums(j-1),mps)/real(2*isums(j)-2*isums(j-1),mps)
8733#ifdef SCOREP_USER_ENABLE
8734 scorep_user_region_by_name_end(
"UR_monres")
8738110
FORMAT(i5,2i10,3g14.5)
8753 INTEGER(mpi) :: ioff
8754 INTEGER(mpi) :: ipar0
8755 INTEGER(mpi) :: ncon
8756 INTEGER(mpi) :: npar
8757 INTEGER(mpi) :: nextra
8759 INTEGER :: nbopt, nboptx, ILAENV
8762 INTEGER(mpl),
INTENT(IN) :: msize(2)
8764 INTEGER(mpl) :: length
8765 INTEGER(mpl) :: nwrdpc
8766 INTEGER(mpl),
PARAMETER :: three = 3
8779 CALL mpalloc(
aux,length,
' local fit scratch array: aux')
8780 CALL mpalloc(
vbnd,length,
' local fit scratch array: vbnd')
8781 CALL mpalloc(
vbdr,length,
' local fit scratch array: vbdr')
8784 CALL mpalloc(
vbk,length,
' local fit scratch array: vbk')
8787 CALL mpalloc(
vzru,length,
' local fit scratch array: vzru')
8788 CALL mpalloc(
scdiag,length,
' local fit scratch array: scdiag')
8789 CALL mpalloc(
scflag,length,
' local fit scratch array: scflag')
8790 CALL mpalloc(
ibandh,2*length,
' local fit band width hist.: ibandh')
8804 length=4*
ncblck;
IF(ncon == 0) length=0
8813 nwrdpc=nwrdpc+length
8825 length=(ncon*ncon+ncon)/2
8850 length=length+npar+nextra
8856 IF(
mbandw == 0) length=length+1
8864 nwrdpc=nwrdpc+2*length
8865 IF (nwrdpc > 250000)
THEN
8867 WRITE(*,*)
'Size of preconditioner matrix:',int(real(nwrdpc,mps)*4.0e-6,mpi),
' MB'
8889 CALL peend(23,
'Aborted, bad matrix index (will exceed 32bit)')
8890 stop
'vmprep: bad index (matrix to large for diagonalization)'
8912 nbopt = ilaenv( 1_mpl,
'DSYTRF',
'U', int(
nagb,mpl), int(
nagb,mpl), -1_mpl, -1_mpl )
8914 print *,
'LAPACK optimal block size for DSYTRF:', nbopt
8915 lplwrk=length*int(nbopt,mpl)
8923 nbopt = ilaenv( 1_mpl,
'DORMQL',
'RN', int(npar,mpl), int(npar,mpl), int(ncon,mpl), int(npar,mpl) )
8924 IF (int(npar,mpl)*int(nbopt,mpl) >
lplwrk)
THEN
8925 lplwrk=int(npar,mpl)*int(nbopt,mpl)
8930 print *,
'LAPACK optimal block size for DORMQL:', nboptx
8949 INTEGER(mpi) :: icoff
8950 INTEGER(mpi) :: ipoff
8953 INTEGER(mpi) :: ncon
8954 INTEGER(mpi) :: nfit
8955 INTEGER(mpi) :: npar
8956 INTEGER(mpi) :: nrank
8957 INTEGER(mpl) :: imoff
8958 INTEGER(mpl) :: ioff1
8976 WRITE(
lunlog,*)
'Shrinkage of global matrix (A->Q^t*A*Q)'
8991 nfit=npar+ncon;
IF (
icelim > 0) nfit=npar-ncon
8993 IF(nfit < npar)
THEN
9011 WRITE(
lunlog,*)
'Inversion of global matrix (A->A^-1)'
9018 IF(nfit /= nrank)
THEN
9019 WRITE(*,*)
'Warning: the rank defect of the symmetric',nfit, &
9020 '-by-',nfit,
' matrix is ',nfit-nrank,
' (should be zero).'
9021 WRITE(lun,*)
'Warning: the rank defect of the symmetric',nfit, &
9022 '-by-',nfit,
' matrix is ',nfit-nrank,
' (should be zero).'
9025 WRITE(*,*)
' --> enforcing SUBITO mode'
9026 WRITE(lun,*)
' --> enforcing SUBITO mode'
9028 ELSE IF(
ndefec == 0)
THEN
9030 WRITE(lun,*)
'No rank defect of the symmetric matrix'
9032 WRITE(lun,*)
'No rank defect of the symmetric block', ib,
' of size', npar
9043 IF(nfit < npar)
THEN
9062 INTEGER(mpi) :: icoff
9063 INTEGER(mpi) :: ipoff
9066 INTEGER(mpi) :: ncon
9067 INTEGER(mpi) :: nfit
9068 INTEGER(mpi) :: npar
9069 INTEGER(mpi) :: nrank
9070 INTEGER(mpl) :: imoff
9071 INTEGER(mpl) :: ioff1
9086 WRITE(
lunlog,*)
'Shrinkage of global matrix (A->Q^t*A*Q)'
9100 nfit=npar+ncon;
IF (
icelim > 0) nfit=npar-ncon
9102 IF(nfit < npar)
THEN
9120 WRITE(
lunlog,*)
'Decomposition of global matrix (A->L*D*L^t)'
9126 IF(nfit /= nrank)
THEN
9127 WRITE(*,*)
'Warning: the rank defect of the symmetric',nfit, &
9128 '-by-',nfit,
' matrix is ',nfit-nrank,
' (should be zero).'
9129 WRITE(lun,*)
'Warning: the rank defect of the symmetric',nfit, &
9130 '-by-',nfit,
' matrix is ',nfit-nrank,
' (should be zero).'
9133 WRITE(*,*)
' --> enforcing SUBITO mode'
9134 WRITE(lun,*)
' --> enforcing SUBITO mode'
9136 ELSE IF(
ndefec == 0)
THEN
9138 WRITE(lun,*)
'No rank defect of the symmetric matrix'
9140 WRITE(lun,*)
'No rank defect of the symmetric block', ib,
' of size', npar
9142 WRITE(lun,*)
' largest diagonal element (LDLt)', evmax
9143 WRITE(lun,*)
' smallest diagonal element (LDLt)', evmin
9152 IF(nfit < npar)
THEN
9174 INTEGER(mpi) :: icoff
9175 INTEGER(mpi) :: ipoff
9178 INTEGER(mpi) :: ncon
9179 INTEGER(mpi) :: nfit
9180 INTEGER(mpi) :: npar
9181 INTEGER(mpl) :: imoff
9182 INTEGER(mpl) :: ioff1
9183 INTEGER(mpi) :: infolp
9203 WRITE(
lunlog,*)
'Shrinkage of global matrix (A->Q^t*A*Q)'
9218 nfit=npar+ncon;
IF (
icelim > 0) nfit=npar-ncon
9220 IF(nfit < npar)
THEN
9237 IF (nfit > npar)
THEN
9240 WRITE(
lunlog,*)
'Factorization of global matrix (A->L*D*L^t)'
9244#ifdef SCOREP_USER_ENABLE
9245 scorep_user_region_by_name_begin(
"UR_dsptrf", scorep_user_region_type_common)
9248#ifdef SCOREP_USER_ENABLE
9249 scorep_user_region_by_name_end(
"UR_dsptrf")
9256 WRITE(
lunlog,*)
'Factorization of global matrix (A->L*L^t)'
9260#ifdef SCOREP_USER_ENABLE
9261 scorep_user_region_by_name_begin(
"UR_dpptrf", scorep_user_region_type_common)
9263 CALL dpptrf(
'U',int(nfit,mpl),
globalmatd(imoff+1:),infolp)
9264#ifdef SCOREP_USER_ENABLE
9265 scorep_user_region_by_name_end(
"UR_dpptrf")
9273 WRITE(lun,*)
'No rank defect of the symmetric matrix'
9275 WRITE(lun,*)
'No rank defect of the symmetric block', ib,
' of size', npar
9279 WRITE(*,*)
'Warning: factorization of the symmetric',nfit, &
9280 '-by-',nfit,
' failed at index ', infolp
9281 WRITE(lun,*)
'Warning: factorization of the symmetric',nfit, &
9282 '-by-',nfit,
' failed at index ', infolp
9283 CALL peend(29,
'Aborted, factorization of global matrix failed')
9284 stop
'mdptrf: bad matrix'
9289 IF (nfit > npar)
THEN
9292 IF(infolp /= 0) print *,
' DSPTRS failed: ', infolp
9294 CALL dpptrs(
'U',int(nfit,mpl),1_mpl,
globalmatd(imoff+1:),&
9296 IF(infolp /= 0) print *,
' DPPTRS failed: ', infolp
9300 IF(nfit < npar)
THEN
9321 INTEGER(mpi) :: icoff
9322 INTEGER(mpi) :: ipoff
9325 INTEGER(mpi) :: ncon
9326 INTEGER(mpi) :: nfit
9327 INTEGER(mpi) :: npar
9328 INTEGER(mpl) :: imoff
9329 INTEGER(mpl) :: ioff1
9330 INTEGER(mpl) :: iloff
9331 INTEGER(mpi) :: infolp
9352 WRITE(
lunlog,*)
'Shrinkage of global matrix (A->Q^t*A*Q)'
9372 nfit=npar+ncon;
IF (
icelim > 0) nfit=npar-ncon
9374 IF(nfit < npar)
THEN
9378 CALL dtrtrs(
'L',
'T',
'N',int(ncon,mpl),1_mpl,
lapackql(iloff+npar-ncon+1:),int(npar,mpl),&
9380 IF(infolp /= 0) print *,
' DTRTRS failed: ', infolp
9382 CALL dormql(
'L',
'T',int(npar,mpl),1_mpl,int(ncon,mpl),
lapackql(iloff+1:),int(npar,mpl),&
9384 IF(infolp /= 0) print *,
' DORMQL failed: ', infolp
9403 IF (nfit > npar)
THEN
9406 WRITE(
lunlog,*)
'Factorization of global matrix (A->L*D*L^t)'
9410#ifdef SCOREP_USER_ENABLE
9411 scorep_user_region_by_name_begin(
"UR_dsytrf", scorep_user_region_type_common)
9413 CALL dsytrf(
'U',int(nfit,mpl),
globalmatd(imoff+1:),int(nfit,mpl),&
9415#ifdef SCOREP_USER_ENABLE
9416 scorep_user_region_by_name_end(
"UR_dsytrf")
9423 WRITE(
lunlog,*)
'Factorization of global matrix (A->L*L^t)'
9427#ifdef SCOREP_USER_ENABLE
9428 scorep_user_region_by_name_begin(
"UR_dpotrf", scorep_user_region_type_common)
9430 CALL dpotrf(
'U',int(nfit,mpl),
globalmatd(imoff+1:),int(npar,mpl),infolp)
9431#ifdef SCOREP_USER_ENABLE
9432 scorep_user_region_by_name_end(
"UR_dpotrf")
9440 WRITE(lun,*)
'No rank defect of the symmetric matrix'
9442 WRITE(lun,*)
'No rank defect of the symmetric block', ib,
' of size', npar
9446 WRITE(*,*)
'Warning: factorization of the symmetric',nfit, &
9447 '-by-',nfit,
' failed at index ', infolp
9448 WRITE(lun,*)
'Warning: factorization of the symmetric',nfit, &
9449 '-by-',nfit,
' failed at index ', infolp
9450 CALL peend(29,
'Aborted, factorization of global matrix failed')
9451 stop
'mdutrf: bad matrix'
9456 IF (nfit > npar)
THEN
9457 CALL dsytrs(
'U',int(nfit,mpl),1_mpl,
globalmatd(imoff+1:),int(nfit,mpl),&
9459 IF(infolp /= 0) print *,
' DSYTRS failed: ', infolp
9461 CALL dpotrs(
'U',int(nfit,mpl),1_mpl,
globalmatd(imoff+1:),int(npar,mpl),&
9463 IF(infolp /= 0) print *,
' DPOTRS failed: ', infolp
9467 IF(nfit < npar)
THEN
9472 CALL dormql(
'L',
'N',int(npar,mpl),1_mpl,int(ncon,mpl),
lapackql(iloff+1:),int(npar,mpl),&
9474 IF(infolp /= 0) print *,
' DORMQL failed: ', infolp
9481 iloff=iloff+int(npar,mpl)*int(ncon,mpl)
9503 INTEGER(mpi) :: icboff
9504 INTEGER(mpi) :: icblst
9505 INTEGER(mpi) :: icoff
9506 INTEGER(mpi) :: icfrst
9507 INTEGER(mpi) :: iclast
9508 INTEGER(mpi) :: ipfrst
9509 INTEGER(mpi) :: iplast
9510 INTEGER(mpi) :: ipoff
9513 INTEGER(mpi) :: ncon
9514 INTEGER(mpi) :: npar
9516 INTEGER(mpl) :: imoff
9517 INTEGER(mpl) :: iloff
9518 INTEGER(mpi) :: infolp
9519 INTEGER :: nbopt, ILAENV
9521 REAL(mpd),
INTENT(IN) :: a(mszcon)
9522 REAL(mpd),
INTENT(OUT) :: emin
9523 REAL(mpd),
INTENT(OUT) :: emax
9534 iloff=iloff+int(npar,mpl)*int(ncon,mpl)
9553 DO icb=icboff+1,icboff+icblst
9560 lapackql(iloff+ipfrst:iloff+iplast)=a(imoff+1:imoff+npb)
9577 nbopt = ilaenv( 1_mpl,
'DGEQLF',
'', int(npar,mpl), int(ncon,mpl), int(npar,mpl), -1_mpl )
9578 print *,
'LAPACK optimal block size for DGEQLF:', nbopt
9579 lplwrk=int(ncon,mpl)*int(nbopt,mpl)
9582#ifdef SCOREP_USER_ENABLE
9583 scorep_user_region_by_name_begin(
"UR_dgeqlf", scorep_user_region_type_common)
9585 CALL dgeqlf(int(npar,mpl),int(ncon,mpl),
lapackql(iloff+1:),int(npar,mpl),&
9587 IF(infolp /= 0) print *,
' DGEQLF failed: ', infolp
9588#ifdef SCOREP_USER_ENABLE
9589 scorep_user_region_by_name_end(
"UR_dgeqlf")
9593 iloff=iloff+int(npar,mpl)*int(ncon,mpl)
9596 IF(emax < emin)
THEN
9624 INTEGER(mpi) :: icoff
9625 INTEGER(mpi) :: ipoff
9627 INTEGER(mpi) :: ncon
9628 INTEGER(mpi) :: npar
9629 INTEGER(mpl) :: imoff
9630 INTEGER(mpl) :: iloff
9631 INTEGER(mpi) :: infolp
9632 CHARACTER (LEN=1) :: transr, transl
9634 LOGICAL,
INTENT(IN) :: t
9653 IF(ncon <= 0 ) cycle
9656#ifdef SCOREP_USER_ENABLE
9657 scorep_user_region_by_name_begin(
"UR_dormql", scorep_user_region_type_common)
9665 DO i=ipoff+1,ipoff+npar
9671 CALL dormql(
'R',transr,int(npar,mpl),int(npar,mpl),int(ncon,mpl),
lapackql(iloff+1:),&
9674 IF(infolp /= 0) print *,
' DORMQL failed: ', infolp
9676 CALL dormql(
'L',transl,int(npar,mpl),int(npar,mpl),int(ncon,mpl),
lapackql(iloff+1:),&
9679 IF(infolp /= 0) print *,
' DORMQL failed: ', infolp
9680#ifdef SCOREP_USER_ENABLE
9681 scorep_user_region_by_name_end(
"UR_dormql")
9685 iloff=iloff+int(npar,mpl)*int(ncon,mpl)
9691include
'mkl_pardiso.f90'
9722 TYPE(mkl_pardiso_handle) :: pt(64)
9724 INTEGER(mpl),
PARAMETER :: maxfct =1
9725 INTEGER(mpl),
PARAMETER :: mnum = 1
9726 INTEGER(mpl),
PARAMETER :: nrhs = 1
9728 INTEGER(mpl) :: mtype
9729 INTEGER(mpl) :: phase
9730 INTEGER(mpl) :: error
9731 INTEGER(mpl) :: msglvl
9735 INTEGER(mpl) :: idum(1)
9737 INTEGER(mpl) :: length
9738 INTEGER(mpi) :: nfill
9739 INTEGER(mpi) :: npdblk
9740 REAL(mpd) :: adum(1)
9741 REAL(mpd) :: ddum(1)
9743 INTEGER(mpl) :: iparm(64)
9744 REAL(mpd),
ALLOCATABLE :: b( : )
9745 REAL(mpd),
ALLOCATABLE :: x( : )
9759#ifdef SCOREP_USER_ENABLE
9760 scorep_user_region_by_name_begin(
"UR_mspd00", scorep_user_region_type_common)
9767 WRITE(*,*)
'MSPARDISO: number of rows to fill up = ', nfill
9780 CALL pardiso_64(pt, maxfct, mnum, mtype, phase, int(npdblk,mpl), adum, idum, idum, &
9781 idum, nrhs, iparm, msglvl, ddum, ddum, error)
9782 IF (error /= 0)
THEN
9783 WRITE(lun,*)
'The following ERROR was detected: ', error
9784 WRITE(*,
'(A,2I10)')
' PARDISO release failed (phase, error): ', phase, error
9785 IF (
ipddbg == 0)
WRITE(*,*)
' rerun with "debugPARDISO" for more info'
9786 CALL peend(40,
'Aborted, other error: PARDISO release')
9787 stop
'MSPARDISO: stopping due to error in PARDISO release'
9804 IF (iparm(1) == 0)
WRITE(lun,*)
'PARDISO using defaults '
9805 IF (iparm(43) /= 0)
THEN
9806 WRITE(lun,*)
'PARDISO: computation of the diagonal of inverse matrix not implemented !'
9814#ifdef SCOREP_USER_ENABLE
9815 scorep_user_region_by_name_end(
"UR_mspd00")
9823 WRITE(
lunlog,*)
'Decomposition of global matrix (A->L*D*L^t)'
9830#ifdef SCOREP_USER_ENABLE
9831 scorep_user_region_by_name_begin(
"UR_mspd11", scorep_user_region_type_common)
9840 WRITE(lun,*)
' iparm(',i,
') =', iparm(i)
9843 CALL pardiso_64(pt, maxfct, mnum, mtype, phase, int(npdblk,mpl),
globalmatd,
csr3rowoffsets,
csr3columnlist, &
9844 idum, nrhs, iparm, msglvl, ddum, ddum, error)
9845#ifdef SCOREP_USER_ENABLE
9846 scorep_user_region_by_name_end(
"UR_mspd11")
9849 WRITE(lun,*)
'PARDISO reordering completed ... '
9850 WRITE(lun,*)
'PARDISO peak memory required (KB)', iparm(15)
9853 WRITE(lun,*)
' iparm(',i,
') =', iparm(i)
9856 IF (error /= 0)
THEN
9857 WRITE(lun,*)
'The following ERROR was detected: ', error
9858 WRITE(*,
'(A,2I10)')
' PARDISO decomposition failed (phase, error): ', phase, error
9859 IF (
ipddbg == 0)
WRITE(*,*)
' rerun with "debugPARDISO" for more info'
9860 CALL peend(40,
'Aborted, other error: PARDISO reordering')
9861 stop
'MSPARDISO: stopping due to error in PARDISO reordering'
9863 IF (iparm(60) == 0)
THEN
9868 WRITE(lun,*)
'Size (KB) of allocated memory = ',
ipdmem
9869 WRITE(lun,*)
'Number of nonzeros in factors = ',iparm(18)
9870 WRITE(lun,*)
'Number of factorization MFLOPS = ',iparm(19)
9874#ifdef SCOREP_USER_ENABLE
9875 scorep_user_region_by_name_begin(
"UR_mspd22", scorep_user_region_type_common)
9878 CALL pardiso_64(pt, maxfct, mnum, mtype, phase, int(npdblk,mpl),
globalmatd,
csr3rowoffsets,
csr3columnlist, &
9879 idum, nrhs, iparm, msglvl, ddum, ddum, error)
9880#ifdef SCOREP_USER_ENABLE
9881 scorep_user_region_by_name_end(
"UR_mspd22")
9884 WRITE(lun,*)
'PARDISO factorization completed ... '
9887 WRITE(lun,*)
' iparm(',i,
') =', iparm(i)
9890 IF (error /= 0)
THEN
9891 WRITE(lun,*)
'The following ERROR was detected: ', error
9892 WRITE(*,
'(A,2I10)')
' PARDISO decomposition failed (phase, error): ', phase, error
9893 IF (
ipddbg == 0)
WRITE(*,*)
' rerun with "debugPARDISO" for more info'
9894 CALL peend(40,
'Aborted, other error: PARDISO factorization')
9895 stop
'MSPARDISO: stopping due to error in PARDISO factorization'
9898 IF (iparm(14) > 0) &
9899 WRITE(lun,*)
'Number of perturbed pivots = ',iparm(14)
9900 WRITE(lun,*)
'Number of positive eigenvalues = ',iparm(22)-nfill
9901 WRITE(lun,*)
'Number of negative eigenvalues = ',iparm(23)
9902 ELSE IF (iparm(30) > 0)
THEN
9903 WRITE(lun,*)
'Equation with bad pivot (<=0.) = ',iparm(30)
9912 CALL mpalloc(b,length,
' PARDISO r.h.s')
9913 CALL mpalloc(x,length,
' PARDISO solution')
9916#ifdef SCOREP_USER_ENABLE
9917 scorep_user_region_by_name_begin(
"UR_mspd33", scorep_user_region_type_common)
9921 CALL pardiso_64(pt, maxfct, mnum, mtype, phase, int(npdblk,mpl),
globalmatd,
csr3rowoffsets,
csr3columnlist, &
9922 idum, nrhs, iparm, msglvl, b, x, error)
9923#ifdef SCOREP_USER_ENABLE
9924 scorep_user_region_by_name_end(
"UR_mspd33")
9930 WRITE(lun,*)
'PARDISO solve completed ... '
9931 IF (error /= 0)
THEN
9932 WRITE(lun,*)
'The following ERROR was detected: ', error
9933 WRITE(*,
'(A,2I10)')
' PARDISO decomposition failed (phase, error): ', phase, error
9934 IF (
ipddbg == 0)
WRITE(*,*)
' rerun with "debugPARDISO" for more info'
9935 CALL peend(40,
'Aborted, other error: PARDISO solve')
9936 stop
'MSPARDISO: stopping due to error in PARDISO solve'
9950 INTEGER(mpi) :: iast
9951 INTEGER(mpi) :: idia
9952 INTEGER(mpi) :: imin
9953 INTEGER(mpl) :: ioff1
9955 INTEGER(mpi) :: last
9957 INTEGER(mpi) :: nmax
9958 INTEGER(mpi) :: nmin
9959 INTEGER(mpi) :: ntop
9981 WRITE(
lunlog,*)
'Shrinkage of global matrix (A->Q^t*A*Q)'
10018 DO WHILE(ntop < nmax)
10022 CALL hmpdef(7,real(nmin,mps),real(ntop,mps),
'log10 of positive eigenvalues')
10032 iast=max(1,imin-60)
10033 CALL gmpdef(3,2,
'low-value end of eigenvalues')
10036 CALL gmpxy(3,real(i,mps),evalue)
10050 WRITE(lun,*)
'The first (largest) eigenvalues ...'
10053 WRITE(lun,*)
'The last eigenvalues ... up to',last
10057 WRITE(lun,*)
'The eigenvalues from',
nvgb+1,
' to',
nagb
10061 WRITE(lun,*)
'Log10 + 3 of ',
nfgb,
' eigenvalues in decreasing',
' order'
10062 WRITE(lun,*)
'(for Eigenvalue < 0.001 the value 0.0 is shown)'
10065 'printed for negative eigenvalues'
10068 WRITE(lun,*) last,
' significances: insignificant if ', &
10069 'compatible with N(0,1)'
10098 INTEGER(mpl) :: ioff1
10099 INTEGER(mpl) :: ioff2
10135 INTEGER(mpi) :: istop
10136 INTEGER(mpi) :: itn
10137 INTEGER(mpi) :: itnlim
10138 INTEGER(mpi) :: lun
10139 INTEGER(mpi) :: nout
10140 INTEGER(mpi) :: nrkd
10141 INTEGER(mpi) :: nrkd2
10147 REAL(mpd) :: arnorm
10186 WRITE(lun,*)
'MMINRS: PRECONS ended ', nrkd
10190 globalcorrections, itnlim, nout, rtol, istop, itn, anorm, acond, rnorm, arnorm, ynorm)
10191 ELSE IF(
mbandw > 0)
THEN
10197 WRITE(lun,*)
'MMINRS: EQUDECS ended ', nrkd, nrkd2
10201 globalcorrections, itnlim, nout, rtol, istop, itn, anorm, acond, rnorm, arnorm, ynorm)
10204 globalcorrections, itnlim, nout, rtol, istop, itn, anorm, acond, rnorm, arnorm, ynorm)
10218 IF (
istopa == 0) print *,
'MINRES: istop=0, exact solution x=0.'
10233 INTEGER(mpi) :: istop
10234 INTEGER(mpi) :: itn
10235 INTEGER(mpi) :: itnlim
10236 INTEGER(mpi) :: lun
10237 INTEGER(mpi) :: nout
10238 INTEGER(mpi) :: nrkd
10239 INTEGER(mpi) :: nrkd2
10242 REAL(mpd) :: mxxnrm
10243 REAL(mpd) :: trcond
10253 mxxnrm = real(
nagb,mpd)/sqrt(epsilon(mxxnrm))
10255 trcond = 1.0_mpd/epsilon(trcond)
10256 ELSE IF(
mrmode == 2)
THEN
10286 WRITE(lun,*)
'MMINRS: PRECONS ended ', nrkd
10290 itnlim=itnlim, rtol=rtol, maxxnorm=mxxnrm, trancond=trcond, &
10292 ELSE IF(
mbandw > 0)
THEN
10298 WRITE(lun,*)
'MMINRS: EQUDECS ended ', nrkd, nrkd2
10303 itnlim=itnlim, rtol=rtol, maxxnorm=mxxnrm, trancond=trcond, &
10307 itnlim=itnlim, rtol=rtol, maxxnorm=mxxnrm, trancond=trcond, &
10322 IF (
istopa == 3) print *,
'MINRES: istop=0, exact solution x=0.'
10338 INTEGER(mpi),
INTENT(IN) :: n
10339 REAL(mpd),
INTENT(IN) :: x(n)
10340 REAL(mpd),
INTENT(OUT) :: y(n)
10360 INTEGER(mpi),
INTENT(IN) :: n
10361 REAL(mpd),
INTENT(IN) :: x(n)
10362 REAL(mpd),
INTENT(OUT) :: y(n)
10393 REAL(mps) :: concu2
10394 REAL(mps) :: concut
10395 REAL,
DIMENSION(2) :: ta
10398 INTEGER(mpi) :: iact
10399 INTEGER(mpi) :: iagain
10400 INTEGER(mpi) :: idx
10401 INTEGER(mpi) :: info
10403 INTEGER(mpi) :: ipoff
10404 INTEGER(mpi) :: icoff
10405 INTEGER(mpl) :: ioff
10406 INTEGER(mpi) :: itgbi
10407 INTEGER(mpi) :: ivgbi
10408 INTEGER(mpi) :: jcalcm
10410 INTEGER(mpi) :: labelg
10411 INTEGER(mpi) :: litera
10412 INTEGER(mpl) :: lrej
10413 INTEGER(mpi) :: lun
10414 INTEGER(mpi) :: lunp
10415 INTEGER(mpi) :: minf
10416 INTEGER(mpi) :: mrati
10417 INTEGER(mpi) :: nan
10418 INTEGER(mpi) :: ncon
10419 INTEGER(mpi) :: nfaci
10420 INTEGER(mpi) :: nloopsol
10421 INTEGER(mpi) :: npar
10422 INTEGER(mpi) :: nrati
10423 INTEGER(mpl) :: nrej
10424 INTEGER(mpi) :: nsol
10425 INTEGER(mpi) :: inone
10427 INTEGER(mpi) :: infolp
10428 INTEGER(mpi) :: nfit
10429 INTEGER(mpl) :: imoff
10433 REAL(mpd) :: dratio
10434 REAL(mpd) :: dwmean
10443 LOGICAL :: warnerss
10444 LOGICAL :: warners3
10446 CHARACTER (LEN=7) :: cratio
10447 CHARACTER (LEN=7) :: cfacin
10448 CHARACTER (LEN=7) :: crjrat
10459 WRITE(lunp,*)
'Solution algorithm: '
10460 WRITE(lunp,121)
'=================================================== '
10463 WRITE(lunp,121)
'solution method:',
'matrix inversion'
10464 ELSE IF(
metsol == 2)
THEN
10465 WRITE(lunp,121)
'solution method:',
'diagonalization'
10466 ELSE IF(
metsol == 3)
THEN
10467 WRITE(lunp,121)
'solution method:',
'decomposition'
10468 ELSE IF(
metsol == 4)
THEN
10469 WRITE(lunp,121)
'solution method:',
'minres (Paige/Saunders)'
10470 ELSE IF(
metsol == 5)
THEN
10471 WRITE(lunp,121)
'solution method:',
'minres-qlp (Choi/Paige/Saunders)'
10473 WRITE(lunp,121)
' ',
' using QR factorization'
10474 ELSE IF(
mrmode == 2)
THEN
10475 WRITE(lunp,121)
' ',
' using QLP factorization'
10477 WRITE(lunp,121)
' ',
' using QR and QLP factorization'
10478 WRITE(lunp,123)
'transition condition',
mrtcnd
10480 ELSE IF(
metsol == 6)
THEN
10481 WRITE(lunp,121)
'solution method:', &
10482 'gmres (generalized minimzation of residuals)'
10484 ELSE IF(
metsol == 7)
THEN
10486 WRITE(lunp,121)
'solution method:',
'LAPACK factorization (DSPTRF)'
10488 WRITE(lunp,121)
'solution method:',
'LAPACK factorization (DPPTRF)'
10490 IF(
ilperr == 1)
WRITE(lunp,121)
' ',
'with error calculation (D??TRI)'
10491 ELSE IF(
metsol == 8)
THEN
10493 WRITE(lunp,121)
'solution method:',
'LAPACK factorization (DSYTRF)'
10495 WRITE(lunp,121)
'solution method:',
'LAPACK factorization (DPOTRF)'
10497 IF(
ilperr == 1)
WRITE(lunp,121)
' ',
'with error calculation (D??TRI)'
10499 ELSE IF(
metsol == 9)
THEN
10501 WRITE(lunp,121)
'solution method:',
'Intel oneMKL PARDISO (sparse matrix (CSR3))'
10503 WRITE(lunp,121)
'solution method:',
'Intel oneMKL PARDISO (sparse matrix (BSR3))'
10508 WRITE(lunp,123)
'convergence limit at Delta F=',
dflim
10509 WRITE(lunp,122)
'maximum number of iterations=',
mitera
10512 WRITE(lunp,122)
'matrix recalculation up to ',
matrit,
'. iteration'
10516 WRITE(lunp,121)
'matrix storage:',
'full'
10517 ELSE IF(
matsto == 2)
THEN
10518 WRITE(lunp,121)
'matrix storage:',
'sparse'
10520 WRITE(lunp,122)
'pre-con band-width parameter=',
mbandw
10522 WRITE(lunp,121)
'pre-conditioning:',
'default'
10523 ELSE IF(
mbandw < 0)
THEN
10524 WRITE(lunp,121)
'pre-conditioning:',
'none!'
10525 ELSE IF(
mbandw > 0)
THEN
10527 WRITE(lunp,121)
'pre-conditioning=',
'skyline-matrix (rank preserving)'
10529 WRITE(lunp,121)
'pre-conditioning=',
'band-matrix'
10534 WRITE(lunp,121)
'using pre-sigmas:',
'no'
10537 WRITE(lunp,124)
'pre-sigmas defined for', &
10538 REAL(100*npresg,mps)/
REAL(nvgb,mps),
' % of variable parameters'
10539 WRITE(lunp,123)
'default pre-sigma=',
regpre
10542 WRITE(lunp,121)
'regularization:',
'no'
10544 WRITE(lunp,121)
'regularization:',
'yes'
10545 WRITE(lunp,123)
'regularization factor=',
regula
10549 WRITE(lunp,121)
'Chi square cut equiv 3 st.dev applied'
10550 WRITE(lunp,123)
'... in first iteration with factor',
chicut
10551 WRITE(lunp,123)
'... in second iteration with factor',
chirem
10552 WRITE(lunp,121)
' (reduced by sqrt in next iterations)'
10555 WRITE(lunp,121)
'Scaling of measurement errors applied'
10556 WRITE(lunp,123)
'... factor for "global" measuements',
dscerr(1)
10557 WRITE(lunp,123)
'... factor for "local" measuements',
dscerr(2)
10560 WRITE(lunp,122)
'Down-weighting of outliers in',
lhuber,
' iterations'
10561 WRITE(lunp,123)
'Cut on downweight fraction',
dwcut
10565121
FORMAT(1x,a40,3x,a)
10566122
FORMAT(1x,a40,3x,i0,a)
10567123
FORMAT(1x,a40,2x,e9.2)
10568124
FORMAT(1x,a40,3x,f5.1,a)
10578 stepl =real(stp,mps)
10590 ELSE IF(
metsol == 2)
THEN
10593 ELSE IF(
metsol == 3)
THEN
10596 ELSE IF(
metsol == 4)
THEN
10599 ELSE IF(
metsol == 5)
THEN
10602 ELSE IF(
metsol == 6)
THEN
10614 WRITE(
lunlog,*)
'Checking feasibility of parameters:'
10615 WRITE(*,*)
'Checking feasibility of parameters:'
10616 CALL feasib(concut,iact)
10618 WRITE(*,102) concut
10619 WRITE(*,*)
' parameters are made feasible'
10620 WRITE(
lunlog,102) concut
10621 WRITE(
lunlog,*)
' parameters are made feasible'
10623 WRITE(*,*)
' parameters are feasible (i.e. satisfy constraints)'
10624 WRITE(
lunlog,*)
' parameters are feasible (i.e. satisfy constraints)'
10632 WRITE(*,*)
'Reading files and accumulating vectors/matrices ...'
10636 WRITE(
lunlog,*)
'Reading files and accumulating vectors/matrices ...'
10653 IF(jcalcm+1 /= 0)
THEN
10662 IF(
iterat /= litera)
THEN
10683 IF (iabs(jcalcm) <= 1)
THEN
10696 IF(3*nrej >
nrecal)
THEN
10698 WRITE(*,*)
'Data records rejected in previous loop: '
10700 WRITE(*,*)
'Too many rejects (>33.3%) - stop'
10701 CALL peend(26,
'Aborted, too many rejects')
10708 IF(abs(
icalcm) == 1)
THEN
10720 ELSE IF(
metsol == 2)
THEN
10722 ELSE IF(
metsol == 3)
THEN
10724 ELSE IF(
metsol == 4)
THEN
10726 ELSE IF(
metsol == 5)
THEN
10728 ELSE IF(
metsol == 6)
THEN
10729 WRITE(*,*)
'... reserved for GMRES (not yet!)'
10732 ELSE IF(
metsol == 7)
THEN
10734 ELSE IF(
metsol == 8)
THEN
10737 ELSE IF(
metsol == 9)
THEN
10751 CALL feasib(concut,iact)
10763 angras=real(db/sqrt(db1*db2),mps)
10764 dbsig=16.0_mpd*sqrt(max(db1,db2))*epsilon(db)
10770 lsflag=lsflag .AND. (db > dbsig)
10777 ELSE IF(
metsol == 2)
THEN
10780 ELSE IF(
metsol == 3)
THEN
10783 ELSE IF(
metsol == 4)
THEN
10786 ELSE IF(
metsol == 5)
THEN
10789 ELSE IF(
metsol == 6)
THEN
10799 IF(db <= -dbsig)
THEN
10800 WRITE(*,*)
'Function not decreasing:',db
10801 IF(db > -1.0e-3_mpd)
THEN
10803 IF (iagain <= 1)
THEN
10804 WRITE(*,*)
'... again matrix calculation'
10808 WRITE(*,*)
'... aborting iterations'
10812 WRITE(*,*)
'... stopping iterations'
10835 IF (
nloopn == nloopsol)
THEN
10841 stepl=real(stp,mps)
10851 WRITE(*,*)
'Result vector containes ', nan,
' NaNs - stop'
10852 CALL peend(25,
'Aborted, result vector contains NaNs')
10859 WRITE(*,*)
'Subito! Exit after first step.'
10864 WRITE(*,*)
'INFO=0 should not happen (line search input err)'
10865 IF (iagain <= 0)
THEN
10870 IF(info < 0 .OR.
nloopn == nloopsol) cycle
10874 CALL feasib(concut,iact)
10875 IF(iact /= 0.OR.
chicut > 1.0)
THEN
10890 IF(sum(
nrejec) /= 0)
THEN
10892 WRITE(*,*)
'Data records rejected in last loop: '
10914 nfit=npar+ncon;
IF (
icelim > 0) nfit=npar-ncon
10915 IF (nfit > npar)
THEN
10918 WRITE(
lunlog,*)
'Inverse of global matrix from LDLt factorization'
10923#ifdef SCOREP_USER_ENABLE
10924 scorep_user_region_by_name_begin(
"UR_dsptri", scorep_user_region_type_common)
10927 IF(infolp /= 0) print *,
' DSPTRI failed: ', infolp
10928#ifdef SCOREP_USER_ENABLE
10929 scorep_user_region_by_name_end(
"UR_dsptri")
10935#ifdef SCOREP_USER_ENABLE
10936 scorep_user_region_by_name_begin(
"UR_dsytri", scorep_user_region_type_common)
10938 CALL dsytri(
'U',int(nfit,mpl),
globalmatd(imoff+1:),int(nfit,mpl),&
10940 IF(infolp /= 0) print *,
' DSYTRI failed: ', infolp
10941#ifdef SCOREP_USER_ENABLE
10942 scorep_user_region_by_name_end(
"UR_dsytri")
10949 WRITE(
lunlog,*)
'Inverse of global matrix from LLt factorization'
10954#ifdef SCOREP_USER_ENABLE
10955 scorep_user_region_by_name_begin(
"UR_dpptri", scorep_user_region_type_common)
10957 CALL dpptri(
'U',int(nfit,mpl),
globalmatd(imoff+1:),infolp)
10958 IF(infolp /= 0) print *,
' DPPTRI failed: ', infolp
10959#ifdef SCOREP_USER_ENABLE
10960 scorep_user_region_by_name_end(
"UR_dpptri")
10965#ifdef SCOREP_USER_ENABLE
10966 scorep_user_region_by_name_begin(
"UR_dpotri", scorep_user_region_type_common)
10968 CALL dpotri(
'U',int(nfit,mpl),
globalmatd(imoff+1:),int(npar,mpl),infolp)
10969 IF(infolp /= 0) print *,
' DPOTRI failed: ', infolp
10970#ifdef SCOREP_USER_ENABLE
10971 scorep_user_region_by_name_end(
"UR_dpotri")
10989 DO i=npar-ncon+1,npar
10996 WRITE(
lunlog,*)
'Expansion of global matrix (A->Q*A*Q^t)'
11012 catio=real(dratio,mps)
11016 mrati=nint(100.0*catio,mpi)
11020 IF (
nfilw <= 0)
THEN
11021 WRITE(lunp,*)
'Sum(Chi^2)/Sum(Ndf) =',
fvalue
11023 WRITE(lunp,*)
' =',dratio
11025 WRITE(lunp,*)
'Sum(W*Chi^2)/Sum(Ndf)/<W> =',
fvalue
11027 WRITE(lunp,*)
' /',dwmean
11028 WRITE(lunp,*)
' =',dratio
11032 ' with correction for down-weighting ',catio
11053 nrati=nint(10000.0*real(nrej,mps)/real(
nrecal,mps),mpi)
11054 WRITE(crjrat,197) 0.01_mpd*real(nrati,mpd)
11055 nfaci=nint(100.0*sqrt(catio),mpi)
11057 WRITE(cratio,197) 0.01_mpd*real(mrati,mpd)
11058 WRITE(cfacin,197) 0.01_mpd*real(nfaci,mpd)
11061 IF(mrati < 90.OR.mrati > 110) warner=.true.
11062 IF(nrati > 100) warner=.true.
11063 IF(
ncgbe /= 0) warner=.true.
11065 IF(
nalow /= 0) warners=.true.
11067 IF(
nmiss1 /= 0) warnerss=.true.
11068 IF(iagain /= 0) warnerss=.true.
11069 IF(
ndefec /= 0) warnerss=.true.
11070 IF(
ndefpg /= 0) warnerss=.true.
11072 IF(
nrderr /= 0) warners3=.true.
11074 IF(warner.OR.warners.OR.warnerss.Or.warners3)
THEN
11077 WRITE(*,199)
'WarningWarningWarningWarningWarningWarningWarningWarningWar'
11078 WRITE(*,199)
'arningWarningWarningWarningWarningWarningWarningWarningWarn'
11079 WRITE(*,199)
'rningWarningWarningWarningWarningWarningWarningWarningWarni'
11080 WRITE(*,199)
'ningWarningWarningWarningWarningWarningWarningWarningWarnin'
11081 WRITE(*,199)
'ingWarningWarningWarningWarningWarningWarningWarningWarning'
11082 WRITE(*,199)
'ngWarningWarningWarningWarningWarningWarningWarningWarningW'
11083 WRITE(*,199)
'gWarningWarningWarningWarningWarningWarningWarningWarningWa'
11085 IF(mrati < 90.OR.mrati > 110)
THEN
11087 WRITE(*,*)
' Chi^2/Ndf = ',cratio,
' (should be close to 1)'
11088 WRITE(*,*)
' => multiply all input standard ', &
11089 'deviations by factor',cfacin
11092 IF(nrati > 100)
THEN
11094 WRITE(*,*)
' Fraction of rejects =',crjrat,
' %', &
11095 ' (should be far below 1 %)'
11096 WRITE(*,*)
' => please provide correct mille data'
11100 IF(iagain /= 0)
THEN
11102 WRITE(*,*)
' Matrix not positiv definite '// &
11103 '(function not decreasing)'
11104 WRITE(*,*)
' => please provide correct mille data'
11109 WRITE(*,*)
' Rank defect =',
ndefec, &
11110 ' for global matrix, should be 0'
11111 WRITE(*,*)
' => please provide correct mille data'
11116 WRITE(*,*)
' Rank defect for',
ndefpg, &
11117 ' parameter groups, should be 0'
11118 WRITE(*,*)
' => please provide correct mille data'
11123 WRITE(*,*)
' Rank defect =',
nmiss1, &
11124 ' for constraint equations, should be 0'
11125 WRITE(*,*)
' => please correct constraint definition'
11128 IF(
ncgbe /= 0)
THEN
11130 WRITE(*,*)
' Number of empty constraints =',
ncgbe,
', should be 0'
11131 WRITE(*,*)
' => please check constraint definition, mille data'
11134 IF(
nxlow /= 0)
THEN
11136 WRITE(*,*)
' Possible rank defects =',
nxlow,
' for global matrix'
11137 WRITE(*,*)
' (too few accepted entries)'
11138 WRITE(*,*)
' => please check mille data and ENTRIES cut'
11141 IF(
nalow /= 0)
THEN
11143 WRITE(*,*)
' Possible bad elements =',
nalow,
' in global vector'
11144 WRITE(*,*)
' (toos few accepted entries)'
11145 IF(
ipcntr > 0)
WRITE(*,*)
' (indicated in millepede.res by counts<0)'
11146 WRITE(*,*)
' => please check mille data and ENTRIES cut'
11151 WRITE(*,*)
' Binary file(s) with read errors =',
nrderr,
' (treated as EOF)'
11152 WRITE(*,*)
' => please check mille data'
11156 WRITE(*,199)
'WarningWarningWarningWarningWarningWarningWarningWarningWar'
11157 WRITE(*,199)
'arningWarningWarningWarningWarningWarningWarningWarningWarn'
11158 WRITE(*,199)
'rningWarningWarningWarningWarningWarningWarningWarningWarni'
11159 WRITE(*,199)
'ningWarningWarningWarningWarningWarningWarningWarningWarnin'
11160 WRITE(*,199)
'ingWarningWarningWarningWarningWarningWarningWarningWarning'
11161 WRITE(*,199)
'ngWarningWarningWarningWarningWarningWarningWarningWarningW'
11162 WRITE(*,199)
'gWarningWarningWarningWarningWarningWarningWarningWarningWa'
11173 ELSE IF(
metsol == 2)
THEN
11175 ELSE IF(
metsol == 3)
THEN
11181 IF(labelg == 0) cycle
11182 itgbi=inone(labelg)
11185 IF(ivgbi < 0) ivgbi=0
11186 IF(ivgbi == 0) cycle
11195 ELSE IF(
metsol == 6)
THEN
11198 ELSE IF(
metsol == 7)
THEN
11206 CALL peend(4,
'Ended with severe warnings (bad binary file(s))')
11207 ELSE IF (warnerss)
THEN
11208 CALL peend(3,
'Ended with severe warnings (bad global matrix)')
11209 ELSE IF (warners)
THEN
11210 CALL peend(2,
'Ended with severe warnings (insufficient measurements)')
11211 ELSE IF (warner)
THEN
11212 CALL peend(1,
'Ended with warnings (bad measurements)')
11214 CALL peend(0,
'Ended normally')
11217102
FORMAT(
' Call FEASIB with cut=',g10.3)
11235 INTEGER(mpi) :: kfl
11236 INTEGER(mpi) :: kmin
11237 INTEGER(mpi) :: kmax
11238 INTEGER(mpi) :: nrc
11239 INTEGER(mpl) :: nrej
11245 REAL(mpd) :: sumallw
11246 REAL(mpd) :: sumrejw
11248 sumallw=0.; sumrejw=0.;
11257 sumallw=sumallw+real(nrc,mpd)*
wfd(kfl)
11258 sumrejw=sumrejw+real(nrej,mpd)*
wfd(kfl)
11259 frac=real(nrej,mps)/real(nrc,mps)
11260 IF (frac > fmax)
THEN
11264 IF (frac < fmin)
THEN
11271 WRITE(*,
"(' Weighted fraction =',F8.2,' %')") 100.*sumrejw/sumallw
11272 IF (
nfilb > 1)
THEN
11273 WRITE(*,
"(' File with max. fraction ',I6,' :',F8.2,' %')") kmax, 100.*fmax
11274 WRITE(*,
"(' File with min. fraction ',I6,' :',F8.2,' %')") kmin, 100.*fmin
11300 INTEGER(mpi) :: iargc
11303 INTEGER(mpi) :: ierrf
11304 INTEGER(mpi) :: ieq
11305 INTEGER(mpi) :: ifilb
11306 INTEGER(mpi) :: ioff
11307 INTEGER(mpi) :: iopt
11308 INTEGER(mpi) :: ios
11309 INTEGER(mpi) :: iosum
11312 INTEGER(mpi) :: mat
11313 INTEGER(mpi) :: nab
11314 INTEGER(mpi) :: nline
11315 INTEGER(mpi) :: npat
11316 INTEGER(mpi) :: ntext
11318 INTEGER(mpi) :: nuf
11319 INTEGER(mpi) :: nums
11320 INTEGER(mpi) :: nufile
11321 INTEGER(mpi) :: lenfileInfo
11322 INTEGER(mpi) :: lenFileNames
11323 INTEGER(mpi) :: matint
11324 INTEGER(mpi),
DIMENSION(:,:),
ALLOCATABLE :: vecfileInfo
11325 INTEGER(mpi),
DIMENSION(:,:),
ALLOCATABLE :: tempArray
11326 INTEGER(mpl) :: rows
11327 INTEGER(mpl) :: cols
11328 INTEGER(mpl) :: newcols
11329 INTEGER(mpl) :: length
11331 CHARACTER (LEN=1024) :: text
11332 CHARACTER (LEN=1024) :: fname
11333 CHARACTER (LEN=14) :: bite(3)
11334 CHARACTER (LEN=32) :: keystx
11335 INTEGER(mpi),
PARAMETER :: mnum=100
11336 REAL(mpd) :: dnum(mnum)
11340 SUBROUTINE initc(nfiles)
BIND(c)
11342 INTEGER(c_int),
INTENT(IN),
VALUE :: nfiles
11343 END SUBROUTINE initc
11348 DATA bite/
'C_binary',
'text ',
'Fortran_binary'/
11363 WRITE(*,*)
'Command line options: '
11364 WRITE(*,*)
'--------------------- '
11366 CALL getarg(i,text)
11367 CALL rltext(text,ia,ib,nab)
11368 WRITE(*,101) i,text(1:nab)
11369 IF(text(ia:ia) /=
'-')
THEN
11370 nu=nufile(text(ia:ib))
11373 WRITE(*,*)
'Second text file in command line - stop'
11374 CALL peend(12,
'Aborted, second text file in command line')
11380 WRITE(*,*)
'Open error for file:',text(ia:ib),
' - stop'
11381 CALL peend(16,
'Aborted, open error for file')
11382 IF(text(ia:ia) /=
'/')
THEN
11383 CALL getenv(
'PWD',text)
11384 CALL rltext(text,ia,ib,nab)
11385 WRITE(*,*)
'PWD:',text(ia:ib)
11390 IF(index(text(ia:ib),
'b') /= 0)
THEN
11392 WRITE(*,*)
'Debugging requested'
11394 it=index(text(ia:ib),
't')
11397 ieq=index(text(ia+it:ib),
'=')+it
11398 IF (it /= ieq)
THEN
11399 IF (index(text(ia+ieq:ib),
'SL0' ) /= 0)
ictest=2
11400 IF (index(text(ia+ieq:ib),
'SLE' ) /= 0)
ictest=3
11401 IF (index(text(ia+ieq:ib),
'BP' ) /= 0)
ictest=4
11402 IF (index(text(ia+ieq:ib),
'BRLF') /= 0)
ictest=5
11403 IF (index(text(ia+ieq:ib),
'BRLC') /= 0)
ictest=6
11406 IF(index(text(ia:ib),
's') /= 0)
isubit=1
11407 IF(index(text(ia:ib),
'f') /= 0)
iforce=1
11408 IF(index(text(ia:ib),
'c') /= 0)
icheck=1
11409 IF(index(text(ia:ib),
'C') /= 0)
icheck=2
11411 IF(i == iargc())
WRITE(*,*)
'--------------------- '
11432 CALL rltext(text,ia,ib,nab)
11433 nu=nufile(text(ia:ib))
11437 CALL peend(10,
'Aborted, no steering file')
11438 stop
'in FILETC: no steering file. .'
11451 WRITE(*,*)
'Listing of steering file: ',
filnam(1:
nfnam)
11452 WRITE(*,*)
'-------------------------'
11455 WRITE(*,*)
'Open error for steering file - stop'
11456 CALL peend(11,
'Aborted, open error for steering file')
11457 IF(
filnam(1:1) /=
'/')
THEN
11458 CALL getenv(
'PWD',text)
11459 CALL rltext(text,ia,ib,nab)
11460 WRITE(*,*)
'PWD:',text(ia:ib)
11469 rows=6; cols=lenfileinfo
11470 CALL mpalloc(vecfileinfo,rows,cols,
'file info from steering')
11473 READ(10,102,iostat=ierrf) text
11474 IF (ierrf < 0)
EXIT
11475 CALL rltext(text,ia,ib,nab)
11477 IF(nline <= 50)
THEN
11478 WRITE(*,101) nline,text(1:nab)
11479 IF(nline == 50)
WRITE(*,*)
' ...'
11483 CALL rltext(text,ia,ib,nab)
11484 IF(ib == ia+2)
THEN
11485 mat=matint(text(ia:ib),
'end',npat,ntext)
11486 IF(mat == max(npat,ntext))
THEN
11489 WRITE(*,*)
' end-statement after',nline,
' text lines'
11494 keystx=
'fortranfiles'
11495 mat=matint(text(ia:ib),keystx,npat,ntext)
11496 IF(mat == max(npat,ntext))
THEN
11503 mat=matint(text(ia:ib),keystx,npat,ntext)
11504 IF(mat == max(npat,ntext))
THEN
11510 keystx=
'closeandreopen'
11511 mat=matint(text(ia:ib),keystx,npat,ntext)
11512 IF(mat == max(npat,ntext))
THEN
11520 iopt=index(text(ia:ib),
' -- ')
11521 IF (iopt > 0) ie=iopt-1
11524 nu=nufile(text(ia:ie))
11526 IF (
nfiles == lenfileinfo)
THEN
11527 CALL mpalloc(temparray,rows,cols,
'temp file info from steering')
11528 temparray=vecfileinfo
11530 lenfileinfo=lenfileinfo*2
11531 newcols=lenfileinfo
11532 CALL mpalloc(vecfileinfo,rows,newcols,
'file info from steering')
11533 vecfileinfo(:,1:cols)=temparray(:,1:cols)
11539 lenfilenames=lenfilenames+ie-ia+1
11540 vecfileinfo(1,
nfiles)=nline
11541 vecfileinfo(2,
nfiles)=nu
11542 vecfileinfo(3,
nfiles)=ia
11543 vecfileinfo(4,
nfiles)=ie
11544 vecfileinfo(5,
nfiles)=iopt
11545 vecfileinfo(6,
nfiles)=ib
11555 CALL mpalloc(
nfd,length,
'file line (in steering)')
11558 length=lenfilenames
11564 READ(10,102,iostat=ierrf) text
11565 IF (ierrf < 0)
EXIT
11567 IF (nline == vecfileinfo(1,i))
THEN
11568 nfd(i)=vecfileinfo(1,i)
11569 mfd(i)=vecfileinfo(2,i)
11570 ia=vecfileinfo(3,i)-1
11571 lfd(i)=vecfileinfo(4,i)-ia
11573 tfd(ioff+k)=text(ia+k:ia+k)
11578 IF (vecfileinfo(5,i) > 0)
THEN
11579 CALL ratext(text(vecfileinfo(5,i)+4:vecfileinfo(6,i)),nums,dnum,mnum)
11580 IF (nums > 0)
ofd(i)=real(dnum(1),mps)
11590 CALL mpalloc(
ifd,length,
'integrated record numbers (=offset)')
11591 CALL mpalloc(
jfd,length,
'number of accepted records')
11592 CALL mpalloc(
kfd,rows,length,
'number of records in file, file order')
11597 CALL mpalloc(
sfd,rows,length,
'start, end of file name in TFD')
11598 CALL mpalloc(
yfd,length,
'modification date')
11601 WRITE(*,*)
'-------------------------'
11607 WRITE(*,*)
'Table of files:'
11608 WRITE(*,*)
'---------------'
11611 WRITE(8,*)
'Text and data files:'
11615 fname(k:k)=
tfd(ioff+k)
11618 IF (
mprint > 1)
WRITE(*,103) i,bite(
mfd(i)),fname(1:
lfd(i))
11619 WRITE(8,103) i,bite(
mfd(i)),fname(1:
lfd(i))
11623 WRITE(*,*)
'---------------'
11637 IF(
mfd(i) == 3)
THEN
11661 IF(
mfd(i) == 1)
THEN
11682 WRITE(*,*)
'Opening of C-files not supported.'
11698 IF(iosum /= 0)
THEN
11699 CALL peend(15,
'Aborted, open error(s) for binary files')
11700 stop
'FILETC: open error '
11702 IF(
nfilb == 0)
THEN
11703 CALL peend(14,
'Aborted, no binary files')
11704 stop
'FILETC: no binary files '
11707 WRITE(*,*)
nfilb,
' binary files opened'
11709 WRITE(*,*)
nfilb,
' binary files opened and closed'
11713103
FORMAT(i3,2x,a14,3x,a)
11776 INTEGER(mpi) :: ierrf
11777 INTEGER(mpi) :: ioff
11778 INTEGER(mpi) :: ios
11779 INTEGER(mpi) :: iosum
11781 INTEGER(mpi) :: mat
11782 INTEGER(mpi) :: nab
11783 INTEGER(mpi) :: nfiln
11784 INTEGER(mpi) :: nline
11785 INTEGER(mpi) :: nlinmx
11786 INTEGER(mpi) :: npat
11787 INTEGER(mpi) :: ntext
11788 INTEGER(mpi) :: matint
11792 CHARACTER (LEN=1024) :: text
11793 CHARACTER (LEN=1024) :: fname
11796 WRITE(*,*)
'Processing text files ...'
11809 IF(
mfd(i) /= 2) cycle
11811 fname(k:k)=
tfd(ia+k)
11813 WRITE(*,*)
'File ',fname(1:
lfd(i))
11814 IF (
mprint > 1)
WRITE(*,*)
' '
11815 OPEN(10,file=fname(1:
lfd(i)),iostat=ios,form=
'FORMATTED')
11817 WRITE(*,*)
'Open error for file ',fname(1:
lfd(i))
11827 READ(10,102,iostat=ierrf) text
11828 IF (ierrf < 0)
THEN
11831 WRITE(*,*)
' end-of-file after',nline,
' text lines'
11835 IF(nline <= nlinmx.AND.
mprint > 1)
THEN
11836 CALL rltext(text,ia,ib,nab)
11838 WRITE(*,101) nline,text(1:nab)
11839 IF(nline == nlinmx)
WRITE(*,*)
' ...'
11842 CALL rltext(text,ia,ib,nab)
11843 IF(ib == ia+2)
THEN
11844 mat=matint(text(ia:ib),
'end',npat,ntext)
11845 IF(mat == max(npat,ntext))
THEN
11848 WRITE(*,*)
' end-statement after',nline,
' text lines'
11854 IF(nfiln <=
nfiles)
THEN
11855 IF(nline ==
nfd(nfiln))
THEN
11870 IF(iosum /= 0)
THEN
11871 CALL peend(16,
'Aborted, open error(s) for text files')
11872 stop
'FILETX: open error(s) in text files '
11875 WRITE(*,*)
'... end of text file processing.'
11880 WRITE(*,*)
lunkno,
' unknown keywords in steering files, ', &
11881 'or file non-existing,'
11882 WRITE(*,*)
' see above!'
11883 WRITE(*,*)
'------------> stop'
11885 CALL peend(13,
'Aborted, unknown keywords in steering file')
11894 ELSE IF(
matsto == 1)
THEN
11896 ELSE IF(
matsto == 2)
THEN
11899 ELSE IF(
metsol == 1)
THEN
11901 ELSE IF(
metsol == 2)
THEN
11903 ELSE IF(
metsol == 3)
THEN
11905 ELSE IF(
metsol == 4)
THEN
11907 ELSE IF(
metsol == 5)
THEN
11909 ELSE IF(
metsol == 6)
THEN
11912 ELSE IF(
metsol == 7)
THEN
11914 ELSE IF(
metsol == 8)
THEN
11917 ELSE IF(
metsol == 9)
THEN
11922 WRITE(*,*)
'MINRES forced with sparse matrix!'
11924 WRITE(*,*)
'MINRES forced with sparse matrix!'
11926 WRITE(*,*)
'MINRES forced with sparse matrix!'
11931 WRITE(*,*)
'MINRES forced with sparse matrix!'
11933 WRITE(*,*)
'MINRES forced with sparse matrix!'
11935 WRITE(*,*)
'MINRES forced with sparse matrix!'
11943 WRITE(*,*)
'Solution method and matrix-storage mode:'
11945 WRITE(*,*)
' METSOL = 1: matrix inversion'
11946 ELSE IF(
metsol == 2)
THEN
11947 WRITE(*,*)
' METSOL = 2: diagonalization'
11948 ELSE IF(
metsol == 3)
THEN
11949 WRITE(*,*)
' METSOL = 3: decomposition'
11950 ELSE IF(
metsol == 4)
THEN
11951 WRITE(*,*)
' METSOL = 4: MINRES'
11952 ELSE IF(
metsol == 5)
THEN
11953 WRITE(*,*)
' METSOL = 5: MINRES-QLP'
11954 ELSE IF(
metsol == 6)
THEN
11955 WRITE(*,*)
' METSOL = 6: GMRES (-> MINRES)'
11957 ELSE IF(
metsol == 7)
THEN
11958 WRITE(*,*)
' METSOL = 7: LAPACK factorization'
11959 ELSE IF(
metsol == 8)
THEN
11960 WRITE(*,*)
' METSOL = 8: LAPACK factorization'
11962 ELSE IF(
metsol == 9)
THEN
11963 WRITE(*,*)
' METSOL = 9: Intel oneMKL PARDISO'
11968 WRITE(*,*)
' with',
mitera,
' iterations'
11971 WRITE(*,*)
' MATSTO = 0: unpacked symmetric matrix, ',
'n*n elements'
11972 ELSEIF(
matsto == 1)
THEN
11973 WRITE(*,*)
' MATSTO = 1: full symmetric matrix, ',
'(n*n+n)/2 elements'
11974 ELSE IF(
matsto == 2)
THEN
11975 WRITE(*,*)
' MATSTO = 2: sparse matrix (custom)'
11976 ELSE IF(
matsto == 3)
THEN
11978 WRITE(*,*)
' MATSTO = 3: sparse matrix (upper triangle, CSR3)'
11980 WRITE(*,*)
' MATSTO = 3: sparse matrix (upper triangle, BSR3)'
11984 WRITE(*,*)
' and band matrix, width',
mbandw
11988 WRITE(*,*)
'Chi square cut equiv 3 st.dev applied ...'
11989 WRITE(*,*)
' in first iteration with factor',
chicut
11990 WRITE(*,*)
' in second iteration with factor',
chirem
11991 WRITE(*,*)
' (reduced by sqrt in next iterations)'
11995 WRITE(*,*)
' Down-weighting of outliers in',
lhuber,
' iterations'
11996 WRITE(*,*)
' Cut on downweight fraction',
dwcut
11999 WRITE(*,*)
'Iterations (solutions) with line search:'
12008 WRITE(*,*)
' All with Chi square cut scaling factor <= 1.'
12039 INTEGER(mpi) :: ios
12043 INTEGER(mpi) :: npat
12044 INTEGER(mpi) :: ntext
12045 INTEGER(mpi) :: nuprae
12048 CHARACTER (LEN=*),
INTENT(INOUT) :: fname
12054 IF(len(fname) > 5)
THEN
12055 IF(fname(1:5) ==
'rfio:') nuprae=1
12056 IF(fname(1:5) ==
'dcap:') nuprae=2
12057 IF(fname(1:5) ==
'root:') nuprae=3
12059 IF(nuprae == 0)
THEN
12060 INQUIRE(file=fname,iostat=ios,exist=ex)
12061 IF(ios /= 0)
nufile=-abs(ios)
12062 IF(ios /= 0)
RETURN
12063 ELSE IF(nuprae == 1)
THEN
12076 nm=
matint(
'xt',fname(l1:ll),npat,ntext)
12079 nm=
matint(
'tx',fname(l1:ll),npat,ntext)
12100 INTEGER(mpi) :: ier
12101 INTEGER(mpi) :: iomp
12104 INTEGER(mpi) :: kkey
12105 INTEGER(mpi) :: label
12106 INTEGER(mpi) :: lkey
12107 INTEGER(mpi) :: mat
12108 INTEGER(mpi) :: miter
12109 INTEGER(mpi) :: nab
12110 INTEGER(mpi) :: nkey
12111 INTEGER(mpi) :: nkeys
12113 INTEGER(mpi) :: nmeth
12114 INTEGER(mpi) :: npat
12115 INTEGER(mpi) :: ntext
12116 INTEGER(mpi) :: nums
12117 INTEGER(mpi) :: matint
12119 CHARACTER (LEN=*),
INTENT(IN) :: text
12120 INTEGER(mpi),
INTENT(IN) :: nline
12124 parameter(nkeys=7,nmeth=10)
12126 parameter(nkeys=6,nmeth=9)
12129 parameter(nkeys=6,nmeth=7)
12131 CHARACTER (LEN=16) :: methxt(nmeth)
12132 CHARACTER (LEN=16) :: keylst(nkeys)
12133 CHARACTER (LEN=32) :: keywrd
12134 CHARACTER (LEN=32) :: keystx
12135 CHARACTER (LEN=itemCLen) :: ctext
12136 INTEGER(mpi),
PARAMETER :: mnum=100
12137 REAL(mpd) :: dnum(mnum)
12140 INTEGER(mpi) :: ipvs
12143 INTEGER(mpi) :: lpvs
12147 SUBROUTINE additem(length,list,label,value)
12149 INTEGER(mpi),
INTENT(IN OUT) :: length
12150 TYPE(listitem),
DIMENSION(:),
INTENT(IN OUT),
ALLOCATABLE :: list
12151 INTEGER(mpi),
INTENT(IN) :: label
12152 REAL(mpd),
INTENT(IN) :: value
12154 SUBROUTINE additemc(length,list,label,text)
12156 INTEGER(mpi),
INTENT(IN OUT) :: length
12157 TYPE(listitemc),
DIMENSION(:),
INTENT(IN OUT),
ALLOCATABLE :: list
12158 INTEGER(mpi),
INTENT(IN) :: label
12159 CHARACTER(LEN = itemCLen),
INTENT(IN) :: text
12161 SUBROUTINE additemi(length,list,label,ivalue)
12163 INTEGER(mpi),
INTENT(IN OUT) :: length
12164 TYPE(listitemi),
DIMENSION(:),
INTENT(IN OUT),
ALLOCATABLE :: list
12165 INTEGER(mpi),
INTENT(IN) :: label
12166 INTEGER(mpi),
INTENT(IN) :: ivalue
12173 DATA keylst/
'unknown',
'parameter',
'constraint',
'measurement',
'method',
'comment',
'pardiso'/
12174 DATA methxt/
'diagonalization',
'inversion',
'fullMINRES',
'sparseMINRES', &
12175 'fullMINRES-QLP',
'sparseMINRES-QLP',
'decomposition',
'fullLAPACK',
'unpackedLAPACK', &
12178 DATA keylst/
'unknown',
'parameter',
'constraint',
'measurement',
'method',
'comment'/
12179 DATA methxt/
'diagonalization',
'inversion',
'fullMINRES',
'sparseMINRES', &
12180 'fullMINRES-QLP',
'sparseMINRES-QLP',
'decomposition',
'fullLAPACK',
'unpackedLAPACK'/
12183 DATA keylst/
'unknown',
'parameter',
'constraint',
'measurement',
'method',
'comment'/
12184 DATA methxt/
'diagonalization',
'inversion',
'fullMINRES',
'sparseMINRES', &
12185 'fullMINRES-QLP',
'sparseMINRES-QLP',
'decomposition'/
12191 CALL rltext(text,ia,ib,nab)
12192 IF(nab == 0)
GOTO 10
12193 CALL ratext(text(1:nab),nums,dnum,mnum)
12195 IF(nums /= 0) nkey=0
12203 keystx=keylst(nkey)
12204 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12205 IF(100*mat >= 80*max(npat,ntext))
GO TO 10
12211 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12212 IF(100*mat >= 80*max(npat,ntext))
THEN
12214 IF(nums > 0) mprint=nint(dnum(1),mpi)
12219 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12220 IF(100*mat >= 80*max(npat,ntext))
THEN
12223 IF(nums > 0) mdebug=nint(dnum(1),mpi)
12224 IF(nums > 1) mdebg2=nint(dnum(2),mpi)
12229 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12230 IF(100*mat >= 80*max(npat,ntext))
THEN
12231 IF(nums > 0 .AND. dnum(1) > 0.5) mreqenf=nint(dnum(1),mpi)
12232 IF(nums > 1 .AND. dnum(2) > 0.5) mreqena=nint(dnum(2),mpi)
12233 IF(nums > 2 .AND. dnum(3) > 0.5) iteren=nint(dnum(1)*dnum(3),mpi)
12237 keystx=
'printrecord'
12238 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12239 IF(100*mat >= 80*max(npat,ntext))
THEN
12240 IF(nums > 0) nrecpr=nint(dnum(1),mpi)
12241 IF(nums > 1) nrecp2=nint(dnum(2),mpi)
12246 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12247 IF(100*mat >= 80*max(npat,ntext))
THEN
12248 IF (nums > 0.AND.dnum(1) > 0.) mxrec=nint(dnum(1),mpi)
12253 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12254 IF(100*mat >= 80*max(npat,ntext))
THEN
12255 IF (nums > 0.AND.dnum(1) >= 0.) ncache=nint(dnum(1),mpi)
12256 IF (nums == 2.AND.dnum(2) > 0..AND.dnum(2) <= 1.0) &
12257 fcache(1)=real(dnum(2),mps)
12258 IF (nums >= 4)
THEN
12260 fcache(k)=real(dnum(k+1),mps)
12267 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12268 IF(100*mat >= 80*max(npat,ntext))
THEN
12273 chicut=real(dnum(1),mps)
12274 IF(chicut < 1.0) chicut=-1.0
12278 chirem=real(dnum(2),mps)
12279 IF(chirem < 1.0) chirem=1.0
12280 IF(chicut >= 1.0) chirem=min(chirem,chicut)
12288 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12289 IF(100*mat >= 80*max(npat,ntext))
THEN
12290 IF(nums > 0) chhuge=real(dnum(1),mps)
12291 IF(chhuge < 1.0) chhuge=1.0
12296 keystx=
'linesearch'
12297 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12298 IF(100*mat >= 80*max(npat,ntext))
THEN
12299 IF(nums > 0) lsearch=nint(dnum(1),mpi)
12304 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12305 IF(100*mat >= 80*max(npat,ntext))
THEN
12306 IF(nums > 0) lfitnp=nint(dnum(1),mpi)
12307 IF(nums > 1) lfitbb=nint(dnum(2),mpi)
12311 keystx=
'regularization'
12312 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12313 IF(100*mat >= 80*max(npat,ntext))
THEN
12315 regula=real(dnum(1),mps)
12316 IF(nums >= 2) regpre=real(dnum(2),mps)
12320 keystx=
'regularisation'
12321 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12322 IF(100*mat >= 80*max(npat,ntext))
THEN
12324 regula=real(dnum(1),mps)
12325 IF(nums >= 2) regpre=real(dnum(2),mps)
12330 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12331 IF(100*mat >= 80*max(npat,ntext))
THEN
12332 regpre=real(dnum(1),mps)
12337 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12338 IF(100*mat >= 80*max(npat,ntext))
THEN
12339 matrit=nint(dnum(1),mpi)
12344 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12345 IF(100*mat >= 80*max(npat,ntext))
THEN
12347 IF (nums > 0.AND.dnum(1) > 0.) matmon=nint(dnum(1),mpi)
12352 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12353 IF(100*mat >= 80*max(npat,ntext))
THEN
12354 IF(nums > 0) mbandw=nint(dnum(1),mpi)
12355 IF(mbandw < 0) mbandw=-1
12356 IF(nums > 1) lprecm=nint(dnum(2),mpi)
12380 keystx=
'outlierdownweighting'
12381 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12382 IF(100*mat >= 80*max(npat,ntext))
THEN
12383 lhuber=nint(dnum(1),mpi)
12384 IF(lhuber > 0.AND.lhuber <= 2) lhuber=2
12388 keystx=
'dwfractioncut'
12389 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12390 IF(100*mat >= 80*max(npat,ntext))
THEN
12391 dwcut=real(dnum(1),mps)
12392 IF(dwcut > 0.5) dwcut=0.5
12396 keystx=
'maxlocalcond'
12397 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12398 IF(100*mat >= 80*max(npat,ntext))
THEN
12399 IF (nums > 0.AND.dnum(1) > 0.0) cndlmx=real(dnum(1),mps)
12404 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12405 IF(100*mat >= 80*max(npat,ntext))
THEN
12406 prange=abs(real(dnum(1),mps))
12411 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12412 IF(100*mat >= 80*max(npat,ntext))
THEN
12418 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12419 IF(100*mat >= 80*max(npat,ntext))
THEN
12424 keystx=
'memorydebug'
12425 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12426 IF(100*mat >= 80*max(npat,ntext))
THEN
12428 IF (nums > 0.AND.dnum(1) > 0.0) memdbg=nint(dnum(1),mpi)
12432 keystx=
'globalcorr'
12433 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12434 IF(100*mat >= 80*max(npat,ntext))
THEN
12439 keystx=
'printcounts'
12440 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12441 IF(100*mat >= 80*max(npat,ntext))
THEN
12443 IF (nums > 0) ipcntr=nint(dnum(1),mpi)
12447 keystx=
'weightedcons'
12448 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12449 IF(100*mat >= 80*max(npat,ntext))
THEN
12451 IF (nums > 0) iwcons=nint(dnum(1),mpi)
12455 keystx=
'skipemptycons'
12456 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12457 IF(100*mat >= 80*max(npat,ntext))
THEN
12462 keystx=
'resolveredundancycons'
12463 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12464 IF(100*mat >= 80*max(npat,ntext))
THEN
12469 keystx=
'withelimination'
12470 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12471 IF(100*mat >= 80*max(npat,ntext))
THEN
12476 keystx=
'postprocessing'
12477 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12478 IF(100*mat >= 80*max(npat,ntext))
THEN
12479 lenpostproc=ib-
keyb-1
12480 cpostproc(1:lenpostproc)=text(
keyb+2:ib)
12485 keystx=
'withLAPACKelimination'
12486 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12487 IF(100*mat >= 80*max(npat,ntext))
THEN
12493 keystx=
'withmultipliers'
12494 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12495 IF(100*mat >= 80*max(npat,ntext))
THEN
12500 keystx=
'checkinput'
12501 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12502 IF(100*mat >= 80*max(npat,ntext))
THEN
12504 IF (nums > 0) icheck=nint(dnum(1),mpi)
12508 keystx=
'checkparametergroups'
12509 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12510 IF(100*mat >= 80*max(npat,ntext))
THEN
12515 keystx=
'monitorresiduals'
12516 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12517 IF(100*mat >= 80*max(npat,ntext))
THEN
12519 IF (nums > 0) imonit=nint(dnum(1),mpi)
12520 IF (nums > 1) measbins=max(measbins,nint(dnum(2),mpi))
12524 keystx=
'monitorpulls'
12525 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12526 IF(100*mat >= 80*max(npat,ntext))
THEN
12529 IF (nums > 0) imonit=nint(dnum(1),mpi)
12530 IF (nums > 1) measbins=max(measbins,nint(dnum(2),mpi))
12534 keystx=
'monitorprogress'
12535 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12536 IF(100*mat >= 80*max(npat,ntext))
THEN
12539 IF (nums > 0) monpg1=max(1,nint(dnum(1),mpi))
12540 IF (nums > 1) monpg2=max(1,nint(dnum(2),mpi))
12544 keystx=
'scaleerrors'
12545 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12546 IF(100*mat >= 80*max(npat,ntext))
THEN
12548 IF (nums > 0) dscerr(1:2)=dnum(1)
12549 IF (nums > 1) dscerr(2)=dnum(2)
12553 keystx=
'iterateentries'
12554 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12555 IF(100*mat >= 80*max(npat,ntext))
THEN
12556 iteren=huge(iteren)
12557 IF (nums > 0) iteren=nint(dnum(1),mpi)
12562 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12563 IF(100*mat >= 80*max(npat,ntext))
THEN
12571 WRITE(*,*)
'WARNING: multithreading not available'
12577 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12578 IF(100*mat >= 80*max(npat,ntext))
THEN
12579 WRITE(*,*)
'WARNING: keyword COMPRESS is obsolete (compression is default)'
12591 keystx=
'countrecords'
12592 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12593 IF(100*mat >= 80*max(npat,ntext))
THEN
12599 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12600 IF(100*mat >= 80*max(npat,ntext).AND.mnrsel < 100)
THEN
12601 nl=min(nums,100-mnrsel)
12603 lbmnrs(mnrsel+k)=nint(dnum(k),mpi)
12609 keystx=
'pairentries'
12610 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12611 IF(100*mat >= 80*max(npat,ntext))
THEN
12614 IF (nums > 0.AND.dnum(1) > 0.0)
THEN
12615 mreqpe=nint(dnum(1),mpi)
12616 IF (nums >= 2.AND.dnum(2) >= dnum(1)) mhispe=nint(dnum(2),mpi)
12617 IF (nums >= 3.AND.dnum(3) >= dnum(1)) msngpe=nint(dnum(3),mpi)
12623 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12624 IF(100*mat >= 80*max(npat,ntext))
THEN
12625 wolfc1=real(dnum(1),mps)
12626 wolfc2=real(dnum(2),mps)
12633 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12634 IF(100*mat >= 80*max(npat,ntext))
THEN
12636 IF (dnum(1) < 1.0e-10_mpd.OR.dnum(1) > 1.0e-04_mpd)
THEN
12637 WRITE(*,*)
'ERROR: need 1.0D-10 <= MRESTL ', &
12638 '<= 1.0D-04, but get ', dnum(1)
12647 keystx=
'mrestranscond'
12648 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12649 IF(100*mat >= 80*max(npat,ntext))
THEN
12657 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12658 IF(100*mat >= 80*max(npat,ntext))
THEN
12660 mrmode = int(dnum(1),mpi)
12665 keystx=
'nofeasiblestart'
12666 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12667 IF(100*mat >= 80*max(npat,ntext))
THEN
12673 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12674 IF(100*mat >= 80*max(npat,ntext))
THEN
12679 keystx=
'readerroraseof'
12680 mat=matint(text(ia:ib),keystx,npat,ntext)
12681 IF(100*mat >= 80*max(npat,ntext))
THEN
12687 keystx=
'LAPACKwitherrors'
12688 mat=matint(text(ia:ib),keystx,npat,ntext)
12689 IF(100*mat >= 80*max(npat,ntext))
THEN
12694 keystx=
'debugPARDISO'
12695 mat=matint(text(ia:ib),keystx,npat,ntext)
12696 IF(100*mat >= 80*max(npat,ntext))
THEN
12701 keystx=
'blocksizePARDISO'
12702 mat=matint(text(
keya:
keyb),keystx,npat,ntext)
12703 IF(100*mat >= 80*max(npat,ntext).AND.mnrsel < 100)
THEN
12704 nl=min(nums,10-mpdbsz)
12706 IF (nint(dnum(k),mpi) > 0)
THEN
12707 IF (mpdbsz == 0)
THEN
12709 ipdbsz(mpdbsz)=nint(dnum(k),mpi)
12710 ELSE IF (nint(dnum(k),mpi) > ipdbsz(mpdbsz))
THEN
12712 ipdbsz(mpdbsz)=nint(dnum(k),mpi)
12720 keystx=
'fortranfiles'
12721 mat=matint(text(ia:ib),keystx,npat,ntext)
12722 IF(mat == max(npat,ntext))
RETURN
12725 mat=matint(text(ia:ib),keystx,npat,ntext)
12726 IF(mat == max(npat,ntext))
RETURN
12728 keystx=
'closeandreopen'
12729 mat=matint(text(ia:ib),keystx,npat,ntext)
12730 IF(mat == max(npat,ntext))
RETURN
12734 IF(nums /= 0) nkey=0
12737 WRITE(*,*)
'**************************************************'
12739 WRITE(*,*)
'Unknown keyword(s): ',text(1:min(nab,50))
12741 WRITE(*,*)
'**************************************************'
1276310
IF(nkey > 0)
THEN
12767 lpvs=nint(dnum(1),mpi)
12769 CALL additem(lenparameters,listparameters,lpvs,dnum(2))
12770 CALL additem(lenpresigmas,listpresigmas,lpvs,dnum(3))
12772 WRITE(*,*)
'Line',nline,
' error, label=',lpvs
12774 ELSE IF(nums /= 0)
THEN
12776 WRITE(*,*)
'Wrong text in line',nline
12777 WRITE(*,*)
'Status: new parameter'
12778 WRITE(*,*)
'> ',text(1:nab)
12780 ELSE IF(lkey == 3)
THEN
12782 IF(nums >= 1.AND.nums <= 2)
THEN
12784 CALL additem(lenconstraints,listconstraints,lpvs,dnum(1))
12786 IF(iwcons > 0) lpvs=-2
12788 IF(nums == 2) plvs=dnum(2)
12789 CALL additem(lenconstraints,listconstraints,lpvs,plvs)
12792 WRITE(*,*)
'Wrong text in line',nline
12793 WRITE(*,*)
'Status: new keyword constraint'
12794 WRITE(*,*)
'> ',text(1:nab)
12796 ELSE IF(lkey == 4)
THEN
12798 nummeasurements=nummeasurements+1
12800 CALL additem(lenmeasurements,listmeasurements,lpvs,dnum(1))
12802 CALL additem(lenmeasurements,listmeasurements,lpvs,dnum(2))
12805 WRITE(*,*)
'Wrong text in line',nline
12806 WRITE(*,*)
'Status: new keyword measurement'
12807 WRITE(*,*)
'> ',text(1:nab)
12809 ELSE IF(lkey == 5.AND.
keyb <
keyc)
THEN
12811 IF(nums >= 1) miter=nint(dnum(1),mpi)
12812 IF(miter >= 1) mitera=miter
12813 dflim=real(dnum(2),mps)
12817 mat=matint(text(
keyb+1:
keyc),keystx,npat,ntext)
12818 IF(100*mat >= 80*max(npat,ntext))
THEN
12822 ELSE IF(i == 2)
THEN
12825 ELSE IF(i == 3)
THEN
12828 ELSE IF(i == 4)
THEN
12831 ELSE IF(i == 5)
THEN
12834 ELSE IF(i == 6)
THEN
12837 ELSE IF(i == 7)
THEN
12841 ELSE IF(i == 8)
THEN
12844 ELSE IF(i == 9)
THEN
12848 ELSE IF(i == 10)
THEN
12857 ELSE IF(nkey == 0)
THEN
12860 lpvs=nint(dnum(1),mpi)
12862 CALL additem(lenparameters,listparameters,lpvs,dnum(2))
12863 CALL additem(lenpresigmas,listpresigmas,lpvs,dnum(3))
12865 WRITE(*,*)
'Line',nline,
' error, label=',lpvs
12867 ELSE IF(nums > 1.AND.nums < 3)
THEN
12869 WRITE(*,*)
'Wrong text in line',nline
12870 WRITE(*,*)
'Status continuation parameter'
12871 WRITE(*,*)
'> ',text(1:nab)
12874 ELSE IF(lkey == 3)
THEN
12877 label=nint(dnum(i),mpi)
12878 IF(label <= 0) ier=1
12880 IF(mod(nums,2) /= 0) ier=1
12883 lpvs=nint(dnum(i),mpi)
12885 CALL additem(lenconstraints,listconstraints,lpvs,plvs)
12889 WRITE(*,*)
'Wrong text in line',nline
12890 WRITE(*,*)
'Status continuation constraint'
12891 WRITE(*,*)
'> ',text(1:nab)
12894 ELSE IF(lkey == 4)
THEN
12898 label=nint(dnum(i),mpi)
12899 IF(label <= 0) ier=1
12901 IF(mod(nums,2) /= 0) ier=1
12905 lpvs=nint(dnum(i),mpi)
12907 CALL additem(lenmeasurements,listmeasurements,lpvs,plvs)
12911 WRITE(*,*)
'Wrong text in line',nline
12912 WRITE(*,*)
'Status continuation measurement'
12913 WRITE(*,*)
'> ',text(1:nab)
12915 ELSE IF(lkey == 6)
THEN
12917 lpvs=nint(dnum(1),mpi)
12921 IF (text(j:j) ==
' ')
EXIT
12924 CALL additemc(lencomments,listcomments,lpvs,ctext)
12926 WRITE(*,*)
'Line',nline,
' error, label=',lpvs
12928 ELSE IF(nums /= 0)
THEN
12930 WRITE(*,*)
'Wrong text in line',nline
12931 WRITE(*,*)
'Status: continuation comment'
12932 WRITE(*,*)
'> ',text(1:nab)
12936 ELSE IF(lkey == 7)
THEN
12939 label=nint(dnum(i),mpi)
12940 IF(label <= 0.OR.label > 64) ier=1
12942 IF(mod(nums,2) /= 0) ier=1
12946 lpvs=nint(dnum(i),mpi)
12947 ipvs=nint(dnum(i+1),mpi)
12948 CALL additemi(lenpardiso,listpardiso,lpvs,ipvs)
12952 WRITE(*,*)
'Wrong text in line',nline
12953 WRITE(*,*)
'Status continuation measurement'
12954 WRITE(*,*)
'> ',text(1:nab)
12973 INTEGER(mpi),
INTENT(IN OUT) :: length
12974 TYPE(
listitem),
DIMENSION(:),
INTENT(IN OUT),
ALLOCATABLE :: list
12975 INTEGER(mpi),
INTENT(IN) :: label
12976 REAL(mpd),
INTENT(IN) :: value
12978 INTEGER(mpl) :: newSize
12979 INTEGER(mpl) :: oldSize
12980 TYPE(
listitem),
DIMENSION(:),
ALLOCATABLE :: tempList
12982 IF (label > 0.AND.
value == 0.)
RETURN
12983 IF (length == 0 )
THEN
12985 CALL mpalloc(list,newsize,
' list ')
12987 oldsize=
size(list,kind=
mpl)
12988 IF (length >= oldsize)
THEN
12989 newsize = oldsize + oldsize/5 + 100
12990 CALL mpalloc(templist,oldsize,
' temp. list ')
12993 CALL mpalloc(list,newsize,
' list ')
12994 list(1:oldsize)=templist(1:oldsize)
12999 list(length)%label=label
13000 list(length)%value=
value
13015 INTEGER(mpi),
INTENT(IN OUT) :: length
13016 TYPE(
listitemc),
DIMENSION(:),
INTENT(IN OUT),
ALLOCATABLE :: list
13017 INTEGER(mpi),
INTENT(IN) :: label
13018 CHARACTER(len = itemCLen),
INTENT(IN) :: text
13020 INTEGER(mpl) :: newSize
13021 INTEGER(mpl) :: oldSize
13022 TYPE(
listitemc),
DIMENSION(:),
ALLOCATABLE :: tempList
13024 IF (label > 0.AND.text ==
'')
RETURN
13025 IF (length == 0 )
THEN
13027 CALL mpalloc(list,newsize,
' list ')
13029 oldsize=
size(list,kind=
mpl)
13030 IF (length >= oldsize)
THEN
13031 newsize = oldsize + oldsize/5 + 100
13032 CALL mpalloc(templist,oldsize,
' temp. list ')
13035 CALL mpalloc(list,newsize,
' list ')
13036 list(1:oldsize)=templist(1:oldsize)
13041 list(length)%label=label
13042 list(length)%text=text
13057 INTEGER(mpi),
INTENT(IN OUT) :: length
13058 TYPE(
listitemi),
DIMENSION(:),
INTENT(IN OUT),
ALLOCATABLE :: list
13059 INTEGER(mpi),
INTENT(IN) :: label
13060 INTEGER(mpi),
INTENT(IN) :: ivalue
13062 INTEGER(mpl) :: newSize
13063 INTEGER(mpl) :: oldSize
13064 TYPE(
listitemi),
DIMENSION(:),
ALLOCATABLE :: tempList
13066 IF (length == 0 )
THEN
13068 CALL mpalloc(list,newsize,
' list ')
13070 oldsize=
size(list,kind=
mpl)
13071 IF (length >= oldsize)
THEN
13072 newsize = oldsize + oldsize/5 + 100
13073 CALL mpalloc(templist,oldsize,
' temp. list ')
13076 CALL mpalloc(list,newsize,
' list ')
13077 list(1:oldsize)=templist(1:oldsize)
13082 list(length)%label=label
13083 list(length)%ivalue=ivalue
13097 CHARACTER (LEN=*),
INTENT(IN) :: text
13098 CHARACTER (LEN=16) :: textc
13107 textl(ka:kb)=text(1:l)
13111 textc=text(1:l)//
'-end'
13119 textl(ka:kb)=textc(1:l)
13146 INTEGER(mpi),
INTENT(IN) :: lun
13147 CHARACTER (LEN=*),
INTENT(IN) :: fname
13148 CHARACTER (LEN=33) :: nafile
13149 CHARACTER (LEN=33) :: nbfile
13155 CALL peend(17,
'Aborted, file name too long')
13156 stop
'File name too long '
13159 nafile(l+1:l+1)=
'~'
13161 INQUIRE(file=nafile(1:l),exist=ex)
13163 INQUIRE(file=nafile(1:l+1),exist=ex)
13165 CALL system(
'rm '//nafile)
13168 nafile(l+1:l+1)=
' '
13169 CALL system(
'mv '//nafile//nbfile)
13171 OPEN(unit=lun,file=fname)
13182 REAL,
DIMENSION(2) :: ta
13202 IF(ncount > 1)
THEN
13204 nsecd1=int(delta,
mpi)
13206 minut1=nsecd1/60-60*nhour1
13207 secnd1=delta-60*(minut1+60*nhour1)
13209 nsecd2=int(delta,
mpi)
13211 minut2=nsecd2/60-60*nhour2
13212 secnd2=delta-60*(minut2+60*nhour2)
13213 WRITE(*,101) nhour1,minut1,secnd1, nhour2,minut2,secnd2
13218101
FORMAT(i4,
' h',i3,
' min',f5.1,
' sec total',18x,
'elapsed', &
13219 i4,
' h',i3,
' min',f5.1,
' sec')
13233 INTEGER(mpi),
INTENT(IN) :: icode
13234 CHARACTER (LEN=*),
INTENT(IN) :: cmessage
13236 CALL mvopen(9,
'millepede.end')
13237 WRITE(9,101) icode, cmessage
13238101
FORMAT(1x,i4,3x,a)
13242END SUBROUTINE peend
13254 INTEGER(mpi),
INTENT(IN) :: kfile
13255 INTEGER(mpi),
INTENT(IN) :: ithr
13256 INTEGER(mpi),
INTENT(OUT) :: ierr
13258 INTEGER(mpi),
DIMENSION(13) :: ibuff
13259 INTEGER(mpi) :: ioff
13260 INTEGER(mpi) :: ios
13262 INTEGER(mpi) :: lfn
13263 INTEGER(mpi) :: lun
13264 INTEGER(mpi) :: moddate
13265 CHARACTER (LEN=1024) :: fname
13266 CHARACTER (LEN=7) :: cfile
13271 SUBROUTINE openc(filename, lfn, lun, ios)
BIND(c)
13273 CHARACTER(kind=c_char),
DIMENSION(*),
INTENT(IN) :: filename
13274 INTEGER(c_int),
INTENT(IN),
VALUE :: lfn
13275 INTEGER(c_int),
INTENT(IN),
VALUE :: lun
13276 INTEGER(c_int),
INTENT(INOUT) :: ios
13277 END SUBROUTINE openc
13289 fname(k:k)=
tfd(ioff+k)
13294 IF(kfile <=
nfilf)
THEN
13297 OPEN(lun,file=fname(1:lfn),iostat=ios, form=
'UNFORMATTED')
13298 print *,
' lun ', lun, ios
13302 CALL openc(fname(1:lfn),lfn,lun,ios)
13304 WRITE(*,*)
'Opening of C-files not supported.'
13311 WRITE(*,*)
'Open error for file ',fname(1:lfn), ios
13312 IF (moddate /= 0)
THEN
13313 WRITE(cfile,
'(I7)') kfile
13314 CALL peend(15,
'Aborted, open error(s) for binary file ' // cfile)
13315 stop
'PEREAD: open error'
13320 ios=stat(fname(1:lfn),ibuff)
13324 WRITE(*,*)
'STAT error for file ',fname(1:lfn), ios
13328 IF (moddate /= 0)
THEN
13329 IF (ibuff(10) /= moddate)
THEN
13330 WRITE(cfile,
'(I7)') kfile
13331 CALL peend(19,
'Aborted, binary file modified (date) ' // cfile)
13332 stop
'PEREAD: file modified'
13335 yfd(kfile)=ibuff(10)
13350 INTEGER(mpi),
INTENT(IN) :: kfile
13351 INTEGER(mpi),
INTENT(IN) :: ithr
13353 INTEGER(mpi) :: lun
13357 SUBROUTINE closec(lun)
BIND(c)
13359 INTEGER(c_int),
INTENT(IN),
VALUE :: lun
13366 IF(kfile <=
nfilf)
THEN
13385 INTEGER(mpi),
INTENT(IN) :: kfile
13387 INTEGER(mpi) :: lun
13391 SUBROUTINE resetc(lun)
BIND(c)
13393 INTEGER(c_int),
INTENT(IN),
VALUE :: lun
13399 IF (kfile <=
nfilf)
THEN
13418 INTEGER(mpi) :: ipgrp
13419 INTEGER(mpi) :: irank
13420 INTEGER(mpi) :: isize
13421 INTEGER(mpi) :: ivoff
13422 INTEGER(mpi) :: itgbi
13424 INTEGER(mpi) :: msize
13425 INTEGER(mpi),
PARAMETER :: mxsize = 1000
13427 INTEGER(mpl):: length
13429 REAL(mpd),
DIMENSION(:),
ALLOCATABLE :: auxVectorD
13430 INTEGER(mpi),
DIMENSION(:),
ALLOCATABLE :: auxVectorI
13431 REAL(mpd),
DIMENSION(:),
ALLOCATABLE :: resParGroup
13432 REAL(mpd),
DIMENSION(:),
ALLOCATABLE :: blockParGroup
13440 IF (isize <= mxsize)
THEN
13441 msize=max(msize,isize)
13443 print *,
' CKPGRP: par. group', ipgrp,
' not checked -- too large: ', isize
13446 IF (msize == 0)
RETURN
13449 length=int(msize,mpl)*(int(msize,mpl)+1)/2
13450 CALL mpalloc(blockpargroup,length,
'(matrix) block for parameter groups (D)')
13452 CALL mpalloc(respargroup,length,
'residuals for parameter groups (D)')
13453 CALL mpalloc(auxvectori,length,
'auxiliary array (I)')
13454 CALL mpalloc(auxvectord,length,
'auxiliary array (D)')
13458 print *,
' CKPGRP par. group first label size rank'
13461 IF (isize > mxsize) cycle
13468 blockpargroup(ij)=matij(ivoff+i,ivoff+j)
13472 CALL sqminv(blockpargroup,respargroup,isize,irank, auxvectord, auxvectori)
13475 IF (isize == irank)
THEN
13479 print *,
' CKPGRP ', ipgrp,
globalparlabelindex(1,itgbi), isize, irank,
' rank deficit !!!'
13497 INTEGER(mpl) :: nan
13498 INTEGER(mpl) :: neg
13500 print *,
' Checking global matrix(D) for NANs ',
size(
globalmatd,kind=mpl)
13505 print *,
' i, nan ', i, nan
13511 print *,
' Checking diagonal elements ',
nagb
13516 print *,
' i, neg ', i, neg
13520 print *,
' CHKMAT summary ', nan, neg
13542 REAL(mpd),
INTENT(IN) :: chi2
13543 INTEGER(mpi),
INTENT(IN) :: ithrd
13544 INTEGER(mpi),
INTENT(IN) :: ndf
13545 REAL(mpd),
INTENT(IN) :: dw
13547 INTEGER(mpl) ::nadd
13575 REAL(mpd),
INTENT(OUT) ::chi2
13576 INTEGER(mpl),
INTENT(OUT) ::ndf
13577 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